Dynamic Modelling and Validation of an Air-to-Water Reversible R744 Heat Pump for High Energy Demand Buildings

This paper presents a reversible heat pump based on CO2 as the refrigerant, able to provide heating, cooling, and domestic hot water to high energy demand buildings. The unit was developed and tested under the EU H2020 project MultiPACK, which has the main goal of assuring the market about the feasibility, reliability, and energy efficiency of CO2 integrated systems for heating and cooling and promoting a fast transition to low environmental impact solutions. Within the project, the confidence raising was performed by installation and monitoring of fully integrated state-of-the art CO2 systems in the Southern European Climate. With the aim of predicting the unit behaviour under variable load and boundary conditions, a dynamic model of the entire unit was developed with commercial software, considering actual components and the implemented control system and it was validated with experimental data, collected at the factory’s lab before commissioning. The validation against experimental data collected during operation as a heat pump demonstrated a maximum percentage difference between the experimental and predicted value of gas–cooler heat flow rate equal to +5.0%. A preliminary comparison with the experimental data in chiller configuration is reported, however further development was required to achieve a satisfactory validation. Lastly, the numerical model was utilized to simulate a typical operation in heat pump configuration with the system coupled with a hot water tank storage for the production of domestic hot water and space heating; the model predicts higher COP when operating in domestic hot water operation due to the lower water inlet temperature.


Introduction
The European Commission [1] has recently confirmed that a renovation wave is needed to green the EU building stock, buildings being responsible for about 40% of the EU's total energy consumption and for 36% of its greenhouse gas emissions from energy. One of the key principles for building renovation towards 2030 and 2050 is energy efficiency, which is classified by EU as the "horizontal guiding principle of European climate and energy governance and beyond to make sure we only produce the energy we really need". According to the same communication, one of the seven areas of intervention is the decarbonisation of heating and cooling. The heating sector offers significant opportunities to decarbonise the energy business [2] and heat pumps appear as one of the most promising way to lower the carbon footprint of the heating system.
The push towards heat pump technology in place of fossil fuels for heating purposes, as previously described, together with the need for a long-term environmentally sustainable refrigerant solution, inspired by regulations, as F-gas regulation in Europe, experimental validation of such models is essential and concluding with a more useful perspective for physical models. CO 2 cycles introduce a higher complexity in the layout design and the need for the further control of high pressure. Rasmussen, 2005 [17] presented the transient model of a carbon dioxide transcritical cycle, where the heat exchangers adopted the moving-boundary approach; the model was validated experimentally. Shi et al., 2010 [18], developed a dynamic model of a transcritical CO 2 cooling unit using the commercial software Dymola to assess the system performance in the absence of physical prototypes. In the same year, Gräber et al., 2010 [19] presented a new library for Dymola that was used to simulate and validate the dynamic model of a CO 2 heat pump. Zheng et al., 2015 [20] developed a transient model of a CO 2 transcritical cycle equipped with an ejector, assuming a 1-D two-phase ejector model, and they evaluated the response of the system to different expansion valve opening degree and area ratios of the ejector. Barta et al., 2019 [21] recently presented the dynamic simulation of a three-stage compression, two-level evaporation CO 2 cycle for military applications, under extreme environmental conditions. The model was developed with Dymola/Modelica language and used to investigate the control system.
Artuso et al., 2020 [22] developed a numerical model of a CO 2 system for road transport refrigeration, using the software AMESim; the energy performance was discussed for different values of internal space temperature and external ambient temperature.
In this paper, pre-commissioning tests of the MultiPACK unit for hotels, carried out at the factory's facilities, are used to validate a dynamic model of the unit, including all its components, i.e., two compressors in parallel, four CO 2 /water heat exchangers, and a CO 2 /air heat exchanger. Considering the complexity and the size of the MultiPACK units, which make experimental activity expensive and time consuming, the model can offer a useful tool to different purposes. Forecasting the performance of the system in a reliable and affordable way, once boundaries are defined, such as thermal load profiles and outdoor temperature, or testing different control logic for operation optimisation and energy consumption minimisation are crucial to promote the successful diffusion of CO 2 based, environmentally friendly HVAC units for high energy consumption buildings. The model was developed with the software AMESim.

The System Layout and the Experimental Setup
The scheme of the refrigerating system is provided in Figure 1, while the main component characteristics are listed in Table 1. The unit has two semi-hermetic compressors (compressor (1) and compressor (2) ), two brazed-plate gas-coolers (gas-cooler (1) and gas-cooler (2) ), a finned-coil heat exchanger (gas-cooler/evaporator (3) ), an internal heat exchanger (IHX) between the compressor suction line and gas-cooler outlet line, a backpressure valve (BPV), a multi-ejector rack, a low-pressure receiver (LPR), and an electronic expansion valve (EEV (1) ). The operation of the system can be switched between two main configurations, i.e., heat pump and chiller, by means of circulation valves located in the circuit. As the refrigerating system is utilized for the production of chilled and hot water in different operating conditions, a more complex configuration was presented compared to the traditional solution. In particular, the system is capable of operating according to a heat pump configuration, for the production of domestic hot water and space heating, and as a chiller configuration for the production of chilled water and eventually hot water. The heat pump configuration consists of a simple configuration air/water heat pump, where the heat is absorbed by the finned coil evaporator and rejected to the brazed plate gas-cooler in order to heat up the incoming water heat flow rate. On the other hand, the chiller configuration has a more complicated operating scheme: a parallel between the back-pressure valve and a multi-ejector system can be employed for the expansion of the refrigerant, to feed an ejector driven evaporator (Evaporator (1) ) and a natural circulation evaporator (Evaporator (2) ), which cool down a water stream matching the temperature profile. An extra brazed-plate gas-cooler is included for domestic hot water production.  Table 1. Geometric dimensions of the system main components as named in Figure 1.

Component Geometric Dimensions/Characteristics
Compressors (1) Compressor (2) Reciprocating, semi-hermetic, single stage compressor.  Table 2 reports the accuracy of the instruments used to monitor the system.

Heat Pump Configuration
During operation as a heat pump, CO 2 is compressed by the two semi-hermetic compressors. One of the compressors is inverter controlled, based on water outlet temperature, while the other compressor operates at a fixed speed. After the compression, the refrigerant in a supercritical state flows in a single-pass brazed-plate gas-cooler before entering the high-pressure side of the internal heat exchanger where it is further cooled down. The refrigerant is then expanded to the evaporation pressure by a back-pressure valve, whose cross-section is modulated by a PI controller to maintain the design high-pressure value p set−point HP , before flowing into the low pressure receiver (LPR) after the evaporation. The internal heat exchanger (IHX) provides saturated vapour superheating before compression. The simple heat pump configuration is reported in Figure 2.

Chiller Configuration
When operating according to the chiller configuration, the system features two evaporators, according to a patented design (EP3486580). One evaporator is fed from the liquid separator, thanks to natural circulation; the flow in the second evaporator is driven by the two-phase ejector block, thus bringing a lower saturation temperature. Purely for testing purposes, the unit can be forced to operate with a single compressor and the ejector driven evaporator, closing the valve BV (1) , thus resulting in the reference ejector layout [4]. The unit standard configuration includes two evaporators and compressors, whose capacity is controlled on cold water delivery temperature, by acting on compressors on/off and inverter frequency. Figure 3 reports the schematic of the refrigerating system operating according to the chiller configuration with a single compressor and a single evaporator. When operating in this configuration the refrigerant is compressed to the transcritical state with a single compressor operating at fixed speed. A brazed plate gas-cooler (gas-cooler (2) ) can be utilized for the hot water production (though it can be bypassed) and the refrigerant is further cooled down in the finned coil gas-cooler. The finned coil fans rotation speed is modulated by a PI controller in order to optimise the temperature at the gas-cooler's outlet, T set−point out−gc . After flowing into the high pressure side of the internal heat exchanger, the refrigerant streams through the motive nozzle, and is expanded in the multi-ejector rack, which has also the role of maintaining the high pressure set point p set−point HP , while the back-pressure valve is kept closed. The refrigerant in the state of saturated liquid flows through the expansion valve (EEV (1) ), whose opening degree is decided by a PI controller in order to maintain a certain set-point pressure at the low-pressure receiver p set−point LPR , and is expanded to the evaporation pressure, before flowing into the brazed-plate evaporator (Evaporator (1) ). The refrigerant exiting the evaporator is recirculated by the ejector rack back to the LPR.

Two Compressors and Two Evaporators
The standard operation model of the chiller includes both semi-hermetic compressors and two evaporators, including the flooded evaporator (Evaporator (2) ). The operation of this configuration is similar to the previous one, however, the additional semi-hermetic compressor is inverter controlled (compressor (2) ) on the water outlet temperature. Furthermore, the expansion of the refrigerant is achieved with a parallel between a back-pressure valve (BPV) and the multi-ejector system. The multi-ejector and the BPV are controlled to maintain the set-point high pressure p set−point HP1 . The two compressors, two evaporators system configuration is reported in Figure 4.

Data Reduction
Experimental data were collected during the operation in steady-state conditions of the refrigerating system in heat pump and chiller configurations. When operating according to the heat pump configuration, the heat flow rate at the gas-cooler was derived on refrigerant side (Equation (1)).
(1) Figure 4. Schematic of the refrigerating system operating according to the chiller configuration, two compressors, and two evaporators.
When the system was operating with both semi-hermetic compressors, CO 2 mass flow rate was not measured by the Coriolis meter, due to its limited FS. The performance of each compressor was then validated against the manufacturer's polynomials independently from the other, with dedicated tests considering different compressor ratio and frequencies.
In Figure 5a, the mass flow rate is reported for both the fixed speed and the inverter controlled compressor. Percentage errors were on average 3%, which is within the tolerance admitted by EN12900:2013 for mass flow declared data and the measurement accuracy. These tests were also used to verify the polynomial curves describing the overall compressor efficiency, thus, the compressor power consumption. When the fixed speed compressor was considered, the experimental reading could be directly compared to the manufacturer's data. On the other side, in the case of inverter controlled compressor, we have to account for the inverter efficiency. This was assumed to be 95%, according to the inverter data sheet. Figure 5b presents the same operating points of Figure 5a in terms of nominal against experimental electrical power consumption. For both the compressors, the power was systematically underestimated by the nominal data. For each frequency lower than 50 Hz, the average error was 5.8%, which is consistent with the tolerance allowed for data declaration (5% according to EN12900:2013) and the measurement accuracy. The 60 Hz actual power was 10.5% higher than the expected. This can be ascribed to a drop in the compressor efficiency at higher speeds and/or a decrease in the inverter performances. This aspect will need further investigation.
The heat flow rate at the internal heat exchanger (Equation (4)), was evaluated as the average value of the heat flow rate in the high (Equation (2)) and low pressure (Equation (3)) side of the heat exchanger. . .
The COP of the heat pump was calculated according to Equation (5).
When the system was operating in chiller configuration, the heat flow rate at the gas-cooler (2) was evaluated according to (Equation (6)): The heat flow rate at the finned coil gas-cooler was evaluated according to Equation (7): Lastly, the heat flow rate at the brazed-plate evaporator (1) was evaluated according to (Equation (8)), when the refrigerant at the evaporator (1) outlet was in a superheated state.
When the refrigerant was in a two-phase state, the heat flow rate at the evaporator was evaluated with (Equation (9)), as the refrigerant enthalpy at evaporator outlet was not defined using pressure and temperature at point 9: As for the evaporator (2) , the heat flow rate at the evaporator was given only on the water side (Equation (10)), as the refrigerant mass flow was not available.
The coefficient of performance of the system operating according to the chiller configuration was then given by Equation (11). Table 3 summarises the maximum uncertainties evaluated during the experiments presented in Section 4; on the high pressure side, whenever the heat flow was calculated on CO 2 side, the uncertainty related to the CO 2 mass flow calculation with compressor curves was always applied, in favour of safety. In the case of the calculation performed on the water side, which occurred whenever the evaporator exit was not superheated, the uncertainty was heavily affected by the very limited water temperature difference.

Numerical Model Development
After the experimental characterization of the refrigerating system operating according to the different configurations, the mathematical model of the system was developed with the commercial software Simcenter AMESim v.17. The libraries of the software are composed of basic elements designed to model the transient behaviour of the internal and external flows, i.e., refrigerant, air and water.
On the refrigerant side, components were modelled based on a homogeneous fluid model, i.e., the two-phase flow was considered as a single phase fluid, displaying the mean fluid properties. Refrigerant and water pressure losses and heat transfer coefficients were evaluated with empirical correlations from the literature. In addition, the model neglected the effect of gravity; the natural convection heat exchanger was thus neglected in this work and preliminary results of the chiller configuration will be presented including only the ejector driven evaporator. The two-phase flow library provided the thermophysical properties of carbon dioxide, water, and air; the thermodynamic properties of CO 2 were evaluated with the Modified Benedict Webb Rubin (MBWR) equation [23], considering the pressure as a function of the fluid's density and temperature. Water was assumed to be pure and its thermodynamic properties were calculated with Bode's formulation, and air was considered as a mixture, composed of a combination of pure fluids weighted by their concentration.
Basic elements were coupled to develop the submodel of the main components of the system, which were connected to form the complete numerical model of the refrigerating unit. The system included mass flow devices (compressor, expansion devices, and twophase ejector) and energy flow devices (evaporators, gas-coolers/condensers, and internal heat exchangers); as the dynamics of the mass flow devices are generally an order of magnitude faster than those of the heat exchangers, according to Rasmussen et al., 2005 [17], the compressor, expansion devices, and two-phase ejector were considered static components and were modelled using steady-state empirical equations. On the other hand, the modelling of the heat exchangers, which are characterized by a more complex nature, was conceived with a discretization of the devices into smaller volumes, each treated with a lumped parameters approach. Finally, the pressure drops of the refrigerant flowing through the pipes connecting different components was neglected, as the assumption of short pipes was made.

Low Pressure Receiver Numerical Model
The low-pressure receiver is a device used to separate liquid from the vapour, storing excess of refrigerant mass to ensure system capacity over a wide range of operating conditions while preventing liquid flow back and consequent damage to the compressor. The submodel of the liquid-vapour separator was conceived assuming an adiabatic constant cross-sectional area; pressure was assumed to be homogeneous in the entire volume and the densities of both the liquid and vapour phase were considered homogeneous as well in their respective volumes. The average refrigerant pressure p LPR and specific enthalpy h LPR were used to compute the state of the refrigerant leaving the receiver (in the state of saturated liquid and vapour, respectively) and the percentage of liquid volume r LV inside of the receiver, which was used to determine the height of the liquid-vapor interface H LV .

Semi-Hermetic Compressor Numerical Model
The semi-hermetic compressor was modelled as a fixed-displacement compressor. The parameters required as input were the value of volumetric efficiency η vol and overall compression efficiency η comp , which were computed from data supplied by the compressor manufacture, the compressor's displacement V d , and the nominal rotatory speed n.
The adherence of the nominal maps has been already discussed in Section 2.2.3, showing a good agreement in terms of mass flow rate and a systematic underestimation of the power consumption of 5%, which increased to 10% when the inverter was set to 60 Hz. As a complete experimental characterization of the compressors in all the working conditions is beyond the scope of this research, the nominal data were used to model the compressor, still being aware of the findings illustrated in Section 2.2.3. for the critical evaluation of the model accuracy. Given the state of the refrigerant at the inlet, the numerical model of the compressor provided the mass-flowrate developed (Equations (12) and (13)), electric power input (Equation (14)), and state of the refrigerant at the compressor outlet. In the case of the fixed speed compressor, which operated at the constant speed of 1450 rpm, the developed mass flow rate was evaluated with Equation (12). In the case of the variable speed compressor, the developed mass flow rate was evaluated with Equation (13) instead, where the value of frequency modulation between 30 and 60 Hz provided by the inverter, Inv ∈ [0, 1] was considered. .

Heat Exchanger Numerical Model
The modelling of the heat exchangers was achieved with an approximation to a simple straight tube in perfect counter flow configuration in a similar manner as discussed in Artuso et al., 2020 [22]. The equivalent tube was characterized by the hydraulic diameter of the heat exchanger (D h = D i for the finned coil heat exchanger and internal heat exchanger and D h = 2b p φ p for the brazed-plate heat exchanger) and the equivalent length which was calculated according to Equation (15).
The equivalent pipe was then discretized into N disc elements (Figure 6), which was decided with a preliminary sensitivity study where N disc was increased until the solution demonstrated insensitivity to the variation of number of elements, and the length of the pipe and heat transfer area were equally distributed in each element. For each discretized volume, it was possible to identify three nodes: a node referring to the refrigerant flow (C r ), a node referring to the state of the heat exchanger's walls (C w ), and a node referring to the state of the air or water, respectively, (C a ). The refrigerant and water flow were considered one-dimensional and the fluid properties varied only on the direction of the flow. In each discretized volume, the refrigerant pressure and density in the associated volume were computed according to Equations (16) and (17), which are given by the mass and energy balance equations applied to the j-th element: From the above equations, the trend of the mean thermodynamic properties of the refrigerant flowing inside the j-th tube element (T r ) were obtained. The internal convective heat flux between the refrigerant and the wall of the heat exchanger in each element is given by Equation (18): where α i is the internal convective heat transfer coefficient, computed from empirical correlations available in the literature, which will be explained later in this section. The thermal-capacity component accounting for the thermal structure of the heat exchanger quantifies the mean wall temperature T (j) w through Equation (19), obtained by the energy conservation equation with the assumption of monodimensional heat flow, as conduction heat flow through wall elements was neglected: The external flow component assessed the mean thermodynamic properties of the air or water flowing outside of the j-th tube element (T a ). In the case of water, the flow was assumed to be one dimensional as with the refrigerant flow, while in the case of air, zero-dimensional flow and negligible conduction in air flow direction was assumed. The convective heat flow of the external flow was computed through Equation (20).
All the empirical correlations utilized in the refrigerant, air side, and water side in the heat exchangers of the system are listed in Tables 4-6.    Friedel, 1979 [27] The inputs required by the submodel were the geometric dimensions of the heat exchanger, and the state of the two fluids, i.e., refrigerant and air/water, at the inlet. The model evaluated, for each discretized element: the state of the refrigerant (T The total heat flow rate and pressure losses were then evaluated with Equations (21) and (22):

Expansion Device Numerical Model
The expansion devices of the system are the back-pressure valve, performing the expansion after the gas cooler, and the expansion valves located between the low pressure receiver and the evaporators. The throttling process was considered isenthalpic, while the effect of choked flow was neglected. The mass flow rate through the throttle was determined by the flow coefficient C q , the actual cross flow area A op , the inlet density, and the pressure difference across the valve ∆p VALVE , according to Equation (23): The model required as inputs the opening ratio of the expansion valve, expressed as the ratio between the actual cross flow area A op and the maximum opening area A op,MAX , previously provided to the model as a constant parameter.

Two-Phase Ejector Numerical Model
The ejector performance was estimated on the basis of the data supplied by the manufacturer and an extrapolation map, relying on experimental data.
The nominal working point, expressed in terms of the entrainment ratio φ nom and pressure lift ∆p li f t nom given by the manufacturer as a function of the motive nozzle and the suction nozzle, stated: where the pressure lift was defined as: An example of the nominal performance working points of the ejector is reported in Figure 7, for a motive nozzle pressure of 90 bar. In order to extend the numerical model to non-nominal conditions, the polynomial mode proposed by Banasiak et al., 2015 [30] was used as a reference. In this paper, a complete polynomial model of an ejector geometry named VEJ1 was developed based on an extensive experimental database. Rearranging the formulation proposed, the entrainment ratio was obtained as a function of the motive and suction thermodynamic condition and of the pressure lift: Equation (28) was used in the model in order to predict the correlation between the lift and the entrainment ratio for values other than the nominal ones: This extrapolation was accurate only when the equilibrium point of the system was close to the nominal point, and decreased as the actual working condition differed from the nominal. This will be discussed in the Results and Validation session.
Given the entrainment ratio from Equation (28) and the motive mass flow rate, the suction nozzle mass flow rate was given by Equation (29). The mass flow rate and specific enthalpy at the diffuser were given by Equations (30) and (31), which enforce the conservation of the mass and the energy: .

Boundary Conditions and Control Systems
Due to the model being a dynamic one, each boundary condition had to be stated as a function of time. For stationary tests, inputs were maintained as constants. The required inputs were the water inlet conditions on both high and low pressure sides, the air inlet conditions at the finned coil heat exchanger, the compressors' status (ON/OFF and, in case, the inverter status), and the mass flow circulating in the natural circulation flooded evaporator.
In the validation phase, these inputs were set according to the corresponding experimental values.

Solver
To evaluate the numerical solution, the software Simcenter AMESim employed a default optimized solver, based on an optimized variable integration method. When the numerical solution was evaluated, to examine its accuracy, a tolerance equal to tol = 10 −7 was chosen; after each computation step, the software integrator checked for the convergence and estimated the error: for each state, variable y i , the error must satisfy the conditions ε i < tol(1 + y i ). When the iterations converged, an error test was applied to determine if the results were accurate enough, however if either of these tests failed, the integration step was replaced with a smaller step size.

Results and Validation
The numerical model was intended to investigate the thermal performance of the refrigerating system operating according to heat pump configuration and chiller configuration and to extend the results coming from the experimental campaign. In order to obtain the model validation, the results obtained from the model were compared against experimental data.

Validation of the Heat Pump Configuration
The validation of the numerical model was provided in both steady state and transient operation. The inputs given to the numerical model were the state of the air and the water, which were collected experimentally during the test, which were m water,1 = 2 kgs −1 , T water,1 = 23.9 • C for the water.
During the steady state operation, the numerical model demonstrated good agreement with the experimental data (Figure 8), while the comparison between the main heat flow rates of the system is reported in Table 7.  The numerical validation in transient conditions (Figure 8) was obtained by suddenly shutting down the compressor (1a) (7.5 min) after the system had previously reached steady state operating condition. Figure 9 reports the trend of . Q gc (1) as a function of time during the steady state condition prior to the compressor shutdown, the transient evolution (7.5-24.5 min), and the following steady state operation (24.5-30 min).
The thermodynamic response of the model was found to be in accordance with the experimental data.
The difference between the experimental and numerical compressor power input was consistent with the data previously presented in Figure 5b, regarding the compressors' characterization. The compressor nominal power consumption has been proven to be in fact 5% to 10% lower than the experimental, leading in this case to an underestimation of power input and consequent and consistent overestimation of the COP.

Preliminary Results in Chiller Configuration
The validation of the numerical model of the refrigerating system operating according to the chiller configuration is discussed here presenting two different cases. In both cases the system runs with a single compressor and a single evaporator: in the first one, the hot water production at the brazed plate gas-cooler is taking place while in the second one it is not.
Differently from the heat pump configuration, where each component of the system was validated ensuring a satisfying accuracy of the numerical model, in the chiller configuration, the model accuracy will depend strongly on the accuracy of the numerical model of the multi-phase ejector.
When the brazed plate gas-cooler is bypassed, the input are given by T I N air = 34.8 • C, ϕ I N air = 84 % for the air side and . m water = 2.4 kgs −1 , T I N water = 12.3 • C for the water side at the evaporator.
Experimental data were used to critically discuss the thermodynamic cycle (Figure 9), the energy flows of the system (Electrical power and heat flows), as well as the operating point of the ejector (Table 8). During the test, the high pressure set point value was set equal to 88.5 bar while the low-pressure receiver set point temperature was set equal to 39.8 bar.
The results of the validation point were consistent with the ones presented for the heat pump operations in the previous section. The ejector equilibrium point was predicted with reasonable accuracy in terms of the recirculation factor but with a lower pressure lift. This caused the increase in the evaporation pressure visible in Figure 10. Except for this, the model obtained a good accuracy in terms of the thermodynamic description of the cycle. The compressor's power consumption was underestimated (in the expected range 5-10%) with respect to the measured one as already discussed previously, while the heat transfer was correctly modelled. The underestimation of the power consumption led to the COP overestimation.  Figure 10. Operation of the chiller configuration with a single evaporator without DHW production.
When the brazed plate heat exchanger was included ( Figure 11, Table 9), the inputs to the numerical model were given by T I N air = 36.3 • C, ϕ I N air = 84 % at the finned coil gas-cooler, . m water = 0.75 kgs −1 , T I N water = 36.8 • C at the brazed-plate gas-cooler, and . m water = 2.6 kgs −1 , T I N water = 13.5 • C at the brazed plate evaporator. In this case, the high pressure set point was set equal to 99.8 bar, while the low-pressure receiver set point operating pressure was equal to 42.3 bar.
As the water flow rate in the gas-cooler was not measured, the boundary value given to the numerical model was derived indirectly from the experimental data according to Equation (32): . m water = . Q EXP gc (2) c p,water T OUT water − T I N water (32) While the high pressure side was well predicted, in this case, the evaporator heat flow was over predicted by 20%. This could be related to the ejector model, which converged for this operating point to a much higher entrainment ratio, leading to a corresponding increase in the evaporator flow rate and then heat exchange. The reported data demonstrated the high impact of the ejector model on the overall system. In contrast to the heat pump configuration, the chiller configuration cannot be considered validated. The following session will then focus on the heat pump, while the improvement of the ejector model will be the object of further research. Figure 11. Operation of the chiller configuration with a single evaporator and DHW production.

Simulation of the System Operating According to a Heat Pump Configuration
The numerical model was utilized to discuss a typical application of the refrigerating unit operating in heat pump configuration with a 24 h simulation, where the system was designed to supply heat to two different water tanks for domestic hot water (DHW) and space heating (SH), as reported in Figure 12.
The inputs into the numerical model were the state of the air at the evaporator inlet (T air,I N , ϕ air,I N ), assumed to be equal to the mean ambient conditions, and the hourly profile of the hot water consumption n . Q load DHW , . Q load SH at the domestic hot water and space heating tanks, respectively.
The reference profiles for the mean ambient temperature and the thermal loads were taken from the work of Smitt et al., 2020 [31], whereas the heating demand values were scaled by a factor of 0.6, in order to adjust them to the nominal heating power of the heat pump presented in this study. In the 24 h simulation, the water flow rate inside the gas-cooler was modulated by a PI controller in order to ensure a certain water outlet temperature depending on SH or DHW production ( Table 10). The controller gains first guess values were obtained applying the Ziegler-Nichols method. Then slight changes were applied to obtain a good compromise between controller error and stability for the specific application. The compressors were turned on and off depending on DHW and SH demand, in particular, during DHW demand both compressors were turned on and the inverter signal was kept at 100%. However, during the SH demand the inverter signal was modulated as the DHW demand decreased, providing energy saving at the compressor. The return water temperature boundary condition changed for SH and DHW and is reported in Table 10. As for the high-pressure control of the heat pump, during the DHW production the high-pressure was kept to a constant value of 105 bar, while during the SH production its value was evaluated depending on the refrigerant temperature at the gas-cooler outlet, according to the study of Chen and Gu., 2005 [32].
The 24 h simulation of the system was solved on a low performance notebook (AMD A8-7410 processor) in approximatively 100 s. Figure 13 reports the trend of the heat flow rate provided at the gas-cooler, . Q gc , as a function of simulation time and DHW/SH demand, while Figure 14 provides the trend of both compressors rotating speed as a function of the same simulation time and DHW/SH demand. It can be observed from both figures that the inverter regulated compressor modulated the rotating speed during the SH production, consequently regulating the heat flow rate at the gas-cooler and the electric power input. On the other hand, Figure 15 reports the trend of high pressure and evaporation pressure during the simulation time. It can be observed that during the SH request time, the value of evaporation pressure stayed between 32 and 35 bar, as the high pressure was kept at approximatively 95 bar, as the refrigerant temperature at the gas-cooler outlet was at a mean value of 37.5 • C, while during DHW production the evaporation temperature was around 30 bar.   As presented in Figure 16, the numerical model demonstrated that, during DHW production, the heat pump presented a higher COP even if the mean pressure ratio between high pressure and evaporation pressure was higher than the one during SH production. The benefit is of course related to the low water temperature at the gas cooler inlet, which results in a higher useful heat flow rate at the gas-cooler. The mean values of high-pressure, evaporation pressure, water mass flow rate developed by the variable speed pump, gascooler heat flow rate, and COP during DHW and SH production are presented in Table 11.

Conclusions
This paper presented the numerical model of a CO 2 unit, which can operate according to a chiller or heat pump configuration. The numerical model was developed with the software Simcenter AMESim v.17, and its results were compared against experimental data in both heat pump and chiller configuration.
The heat pump configuration of the numerical model was validated in both steadystate and dynamic conditions. In the steady-state condition, the model was in accordance with the experimental data in operating points and gas-cooler heat flow rate. The model was also validated in dynamic conditions by suddenly shutting off one of the compressors; the model was able to correctly replicate the heat pump dynamic behaviour, with a maximum difference between the experimentally collected and predicted gas-cooler heat flow rate of 5.0%.
The numerical model assessment against experimental data was also performed in chiller configuration. The validation process demonstrated inconsistent accuracy between the two tested conditions. Furthermore, as the chiller configuration was much more complex than the heat pump one, due to the use of a flooded evaporator and two-phase ejector, further investigation is needed to provide a robust validation under different operating conditions.
Lastly, the numerical model was utilized to investigate the dynamic performance of the refrigerating system operating in heat pump configuration in a typical application where the gas-cooler was coupled with hot water tanks. The model demonstrated a higher COP when operating in DHW operations due to a higher value of gas-cooler heat flow rate.
Future work will consider: • The validation of the numerical model of the system operating in chiller configuration in dynamic conditions and the development of a predictive model of the natural circulation evaporator; • An experimental campaign in order to increase the accuracy of the model in the prediction of the performance of the system. • The coupling of the system with a validated numerical model of water tanks, in order to correctly quantify the dynamic behaviour of the coupled system and discuss the optimal control strategy; • The prediction of the system performance in a yearly simulation.

•
The use of the numerical model in a theoretical study to increase the efficiency and optimize the system during operation in chiller and heat pump configuration.

Acknowledgments:
The experimental activity presented in this manuscript was carried out within the project MultiPACK.

Conflicts of Interest:
The authors declare no conflict of interest.