Dynamic Model of a Transcritical CO 2 Heat Pump for Residential Water Heating

: This paper presents a distributed mathematical model for a carbon dioxide direct expansion solar-assisted heat pump used to heat bath water. The main components are a gas cooler, a needle valve, an evaporator/collector, and a compressor. To develop the heat exchange models, mass, energy, and momentum balances were used. The model was validated for transient as well as steady state conditions using experimental data. A reasonably good agreement was observed between the predicted temperatures and experimental data. The simulations showed that the time step required to demonstrate the behavior of the heat pump in the transient regime is greater than the time step required for the steady state. The results obtained with the mathematical model revealed that a reduction in the water mass ﬂow rate results in an increase in the water outlet temperature. In addition, when the carbon dioxide mass ﬂow rate is reduced, the compressor inlet and outlet temperatures increase as well as the water outlet temperature. software: R.N.O.; supervision: J.J.G.P., A.A.T.M., and R.N.O.; validation: R.N.O.; visualization: R.N.O.; writ-ing—original draft: H.A.G.D., T.F.P., and R.N.O.; writing—review and editing: J.J.G.P. and A.A.T.M. All authors have read and agreed to the published version of the manuscript.


Introduction
The increase in the electrical energy consumption reflects the improvement in the quality of life as well as economic growth. Between 1974 and 2017, world gross production of electricity increased annually at an average rate of 3.3%, and it has grown continuously since 1974. In addition, the total electricity produced from fossil fuels accounted for 64.5% in 2017 [1]. The emission of greenhouse gases into the atmosphere caused an increase in the Earth's average temperature with a consequent effect on climate change [2]. According to the World Health Organization, more than seven million people die every year due to poor air quality, with three million of these deaths being premature [3]. The concern with global warming and the reduction of CO 2 emissions observed in recent years have become research topics worldwide [4]. Currently, different researchers have worked on solutions to increase the energy efficiency of the equipment used in industrial, commercial, and domestic applications. In many cases, the aim of these studies was to contribute to the sustainable development, and in this context, the heat pump appears as an alternative for heating water for residential use. In developed countries, there is an increasing demand for heat pumps [5]. However, in Brazil, the electric shower [6] is still extensively used to heat water. The use of a heat pump to produce hot water could result in a significant electrical energy saving, since it absorbs the environment thermal energy and consumes only a fraction of the electricity consumed by the electric shower [7][8][9].
The vapor compression cycle was developed in 1834; since then, different refrigerant fluids have been used, in both refrigerating and heat pump systems. In the beginning, better COPs than a non-solar assisted heat pump [18]. A study about refrigerants selection for a DX-SAHP for domestic hot water was presented in [19] by the authors evaluated alternatives to replace R134a, and the results showed that propane followed by CO 2 are a good option.
The use of the CO 2 in a DX-SAHP was recently reported in the literature by Oliveira et al. [20], who developed a distributed mathematical model to study the transient condition in the gas cooler of a CO 2 DX-SAHP. That mathematical model was validated using experimental data. The difference noted between the experimental and theoretical data of the inlet and outlet CO 2 temperatures at the gas cooler was 1 • C. Faria et al. [21] also developed a distributed mathematical model to study the transient condition in the evaporator and expansion device of a CO 2 DX-SAHP. Their mathematical model was validated using experimental data with good agreement. Paulino et al. [22] studied the effects of an abrupt change in the solar radiation on a CO 2 DX-SAHP. The distributed model was, once again, validated using experimental data. The results showed that after a variation in the solar radiation, the expansion device needs to act almost immediately to maintain the superheat in the correct value range. In addition, Rabelo et al. [23,24] presented two experimental studies with a CO 2 DX-SAHP. In the first one, the authors discussed the influence of water mass flow rate on the performance of the equipment. The results showed that the increase in the mass flow rate causes a reduction in the outlet water temperature, in the CO 2 outlet temperature, and in the evaporating temperature. In the second study, the authors demonstrated that there is an expansion valve opening at which the heat pump operates at an optimum high pressure and, consequently, at a maximum COP, regardless of the environmental conditions.
The development of mathematical models to simulate the transient and steady state regimes in CO 2 DX-SAHPs is a topic that still demands continued scientific effort. It should be highlighted that due to the fact that the CO 2 heat pump operates in a transcritical cycle, its thermodynamic and transport properties are different from other refrigerants [21]. Additionally, a DX-SAHP is more affected by a solar radiation variation than a non-solar assisted heat pump.
Based on a literature review, it was observed that some theoretical and experimental studies have already been made about CO 2 DX-SAHP. However, none of them have presented a distributed model to study, in a same simulation configuration, for the behavior of all components of a CO 2 DX-SAHP. Thus, the objective of this paper is to study the Sustainability 2021, 13, 3464 3 of 17 dynamic behavior of the CO 2 DX-SAHP with a more inclusive approach. For this purpose, a distributed model is used to find, as a typical application, the steady state of the system as well as the transient response after a variation in the water mass flow rate and in the internal area of the expansion device, followed by an experimental validation. This paper presented here also differs from the others in that it uses the water outlet temperature as the model input data. This type of treatment is most suitable for heating bath water and allows for a better analysis of COP variation under different environmental conditions. Furthermore, the computational effort of the model is assessed.

Prototype Description
The prototype used in this study was previously presented by Faria et al. [21] and by Oliveira et al. [20]. In this study, the equipment needs to produce hot water for a family of four people. For this purpose, it is necessary to heat 200 L of water up to 40 • C [25,26]. A schematic of the prototype heat pump is given in Figure 1, where "T" and "P" are the thermocouples and pressure sensors, respectively. The temperature is measured with Ttype thermocouples, with an accuracy of ±0.5 • C. Three FuelTech pressure sensors with an accuracy of ±0.5% (from 0 to 102 bar) were used to measure the pressure at the compressor suction, the compressor discharge, and the gas cooler outlet. A mechanical pressure gauge ZURICH, with an accuracy of ±0.25% (from 0.1 to 200 bar), was installed at the compressor inlet. The water volumetric flow rate is measured with a beaker and a chronometer with an accuracy of ±2.5%. Based on a literature review, it was observed that some theoretical and experimental studies have already been made about CO2 DX-SAHP. However, none of them have presented a distributed model to study, in a same simulation configuration, for the behavior of all components of a CO2 DX-SAHP. Thus, the objective of this paper is to study the dynamic behavior of the CO2 DX-SAHP with a more inclusive approach. For this purpose, a distributed model is used to find, as a typical application, the steady state of the system as well as the transient response after a variation in the water mass flow rate and in the internal area of the expansion device, followed by an experimental validation. This paper presented here also differs from the others in that it uses the water outlet temperature as the model input data. This type of treatment is most suitable for heating bath water and allows for a better analysis of COP variation under different environmental conditions. Furthermore, the computational effort of the model is assessed.

Prototype Description
The prototype used in this study was previously presented by Faria et al. [21] and by Oliveira et al. [20]. In this study, the equipment needs to produce hot water for a family of four people. For this purpose, it is necessary to heat 200 L of water up to 40 °C [25,26]. A schematic of the prototype heat pump is given in Figure 1, where "T" and "P" are the thermocouples and pressure sensors, respectively. The temperature is measured with Ttype thermocouples, with an accuracy of ±0.5 °C. Three FuelTech pressure sensors with an accuracy of ±0.5% (from 0 to 102 bar) were used to measure the pressure at the compressor suction, the compressor discharge, and the gas cooler outlet. A mechanical pressure gauge ZURICH, with an accuracy of ±0.25% (from 0.1 to 200 bar), was installed at the compressor inlet. The water volumetric flow rate is measured with a beaker and a chronometer with an accuracy of ±2.5%.  The CO 2 compressor model chosen is the EK 6210 CD manufactured by Embraco. The compressor is the hermetic and reciprocating type, has constant speed, operates with a voltage of 110-127 V, and has a displacement of 1.75 cm 3 .rev −1 and a maximum operating pressure of 12 MPa. This equipment was chosen because of its capability of working at high pressures and its high cooling capacity, which consequently generates a high heating capacity.
The gas cooler was designed according to the methodology proposed by Cavallini et al. [27], which suggests temperature differences between the fluids at the inlet and outlet of 20 and 5 • C, respectively. The use of energy conservation equations resulted in a gas cooler tube length of 21 m. The tube was arranged in a spiral shape around the water storage tank. This configuration made the heat pump assembly more compact.
The expansion device is a needle valve model SS-31RS4, manufactured by Swagelok with an orifice area of 1.6 mm 2 .
The evaporator/collector is a flat solar plate panel made of copper with an area of 1.57 m 2 . It receives energy from solar radiation, from natural convection, as well as from the condensation of water vapor present in the atmospheric air. The total length of the tube (16.3 m) was obtained according to the methodology presented by Cheng et al. [28]. The data acquisition system receives the signals generated by the thermocouples and by the pressure sensors. A LabVIEW code was used to save and to show in real time the temperatures and pressures. Figure 2 shows a photograph of the test facility. The CO2 compressor model chosen is the EK 6210 CD manufactured by Embraco. The compressor is the hermetic and reciprocating type, has constant speed, operates with a voltage of 110-127 V, and has a displacement of 1.75 cm³.rev −1 and a maximum operating pressure of 12 MPa. This equipment was chosen because of its capability of working at high pressures and its high cooling capacity, which consequently generates a high heating capacity.
The gas cooler was designed according to the methodology proposed by Cavallini et al. [27], which suggests temperature differences between the fluids at the inlet and outlet of 20 and 5 °C, respectively. The use of energy conservation equations resulted in a gas cooler tube length of 21 m. The tube was arranged in a spiral shape around the water storage tank. This configuration made the heat pump assembly more compact.
The expansion device is a needle valve model SS-31RS4, manufactured by Swagelok with an orifice area of 1.6 mm².
The evaporator/collector is a flat solar plate panel made of copper with an area of 1.57 m². It receives energy from solar radiation, from natural convection, as well as from the condensation of water vapor present in the atmospheric air. The total length of the tube (16.3 m) was obtained according to the methodology presented by Cheng et al. [28]. The data acquisition system receives the signals generated by the thermocouples and by the pressure sensors. A LabVIEW code was used to save and to show in real time the temperatures and pressures. Figure 2 shows a photograph of the test facility.

Model Descriptions
To simulate the transient and the steady-state behavior of the CO2 heat pump, a distributed dynamic model was written in Fortran. This section presents the mathematical model of the evaporator, compressor, gas cooler, and expansion device. Furthermore, it explains how these models are coupled to obtain the overall CO2 heat pump dynamic model.

Compressor Model
The compressor model has the function of providing the outlet mass flow rate for the evaporator model, and the mass flow rate and enthalpy at the gas cooler inlet. In order to formulate the mathematical model for the compressor, the following premises were considered: (i) an irreversible adiabatic compression, (ii) a negligible pressure loss at the suction and discharge, (iii) a constant mass flow rate, and (iv) a refrigerant without any contamination of lubricating oil.

Model Descriptions
To simulate the transient and the steady-state behavior of the CO 2 heat pump, a distributed dynamic model was written in Fortran. This section presents the mathematical model of the evaporator, compressor, gas cooler, and expansion device. Furthermore, it explains how these models are coupled to obtain the overall CO 2 heat pump dynamic model.

Compressor Model
The compressor model has the function of providing the outlet mass flow rate for the evaporator model, and the mass flow rate and enthalpy at the gas cooler inlet. In order to formulate the mathematical model for the compressor, the following premises were considered: (i) an irreversible adiabatic compression, (ii) a negligible pressure loss at the suction and discharge, (iii) a constant mass flow rate, and (iv) a refrigerant without any contamination of lubricating oil.
The isentropic efficiency (η is ), of the compressor presented in Equation (1) was obtained using a multiple regression as described by Sarkar et al. [29]. Then, using Equation (2), the compressor outlet enthalpy (h comp 2 ) was obtained, which is also the value of the gas cooler inlet enthalpy (h gc 1 ).
where a 1 to a 5 are the regression coefficients, P refers to pressure, and the subscripts comp, 1, 2, and s, refer to compressor, inlet, outlet, and isentropic process, respectively. Similarly, using a multiple regression, the compressor volumetric efficiency (η vol ) was obtained as presented Equation (3). Then, it was possible to find the compressor mass flow rate ( . m comp ) as shown in Equation (4).
where b 1 to b 3 are the regression coefficients, N is the rotation speed, V is the volumetric displacement, and v is the specific volume. It is important to highlight that the mass flow rate presented in Equation (4) is the same in evaporator outlet and in the gas cooler inlet, as previously discussed.

Expansion Device Model
The expansion valve model provides the mass flow rate at the gas cooler outlet and at the evaporator inlet. To obtain a simplified numerical simulation, the following assumptions were made: (i) the process is adiabatic, (ii) the flow is one-dimensional, (iii) the relation between orifice of the expansion device (d eq ) and number of turns of the valve (nv) is linear, and (iv) the refrigerant is without any contamination of lubricating oil.
The Equation (5) provides the function to calculate the mass flow rate through the expansion device. .
where T is temperature, ρ is the density, and the subscripts ed and c refer to expansion device and critical point, respectively. Using the Buckham-Pi theorem, it was possible to obtain Equation (6), for the expansion device mass flow rate [30]. The parameters k 1 to k 3 were obtained experimentally, which makes the mathematical model semi-empirical. Martin et al. [31] obtained these parameters for a CO 2 system operating under similar conditions for pressure and mass flow rate.

Heat Exchangers Models
The heat exchangers models have the function of providing the spatial profile of pressure, temperature, enthalpy, specific volume, and mass flow rate for each time step. The following simplifying assumptions were considered: (i) the flows of refrigerant and water were one dimensional, (ii) in the two-phase region, the liquid and vapor were in thermodynamic equilibrium, (iii) the refrigerant was without any contamination by lubricating oil, (iv) there was no axial heat conduction, (v) the radial temperature gradient along the tube wall was negligible, and (vi) the external heat transfer coefficient was uniform. In case of the gas cooler, there was no heat loss to the ambient, and the internal tube is exactly in center of the external tube. For the evaporator, the fin was considered well fixed to the tube. Therefore, there was no contact resistance.
The heat exchangers models can be separated in external and internal balance. For both of them, the solution is similar. The heat exchanger was divided into "n" control volumes and the balances of mass, energy, and momentum applied for the refrigerant in each coil section generate the Equations (7)- (9).
where G and P are the mass flux and the refrigerant modified pressure, respectively. The terms A i and p i are the cross-section area and the internal perimeter of the tube. H i is the heat transfer coefficient between the wall and the refrigerant. This coefficient was calculated for the gas cooler using the methodology proposed by Gnielinski [32], and for the evaporator, the methodology proposed by Cheng et al. [28,33]. (dP/dz) f is the pressure loss due to the friction. In the boiling region, the pressure loss was calculated by the correlation proposed by Friedel and recommended by Cheng et al. [28]. In the superheat region, the pressure loss was given by the Fanning equation, and the friction factor was determined by Blasius equation as described by Ozisik [34]. The subscripts f and p refer to the refrigerant and wall. The superscript "0" refers to the previous instant of time. Finally, ∆t is the time step. The input data for the first control volume are the mass flow rate and the enthalpy provided by the compressor model (for the gas cooler) or by the expansion device model (for the evaporator). Furthermore, the assumptions of the inlet pressure and the temperature of the tube wall profile are important as initial conditions. After the first control volume, the equations are solved by the finite difference method in which the output conditions for the current control volume are the input conditions for the next one. In the last control volume, the outlet mass flow rate is compared to that provided by the compressor (for the gas cooler) or by the expansion device (for the evaporator). If the maximum error is greater than 0.1%, the process is repeated until the convergence occurs. After the second iteration, the Newton-Raphson algorithm is used [6] to accelerate the convergence. Then, having achieved the first convergence, the model carries out the external energy balance, from which a tube wall temperature profile is obtained and compared with that was assumed initially. If this temperature profile presents a maximum error greater than 0.01 or 0.1 • C [6] for the phase change regions and for the region close to the critical point, the internal energy balance is restarted. The current temperature of the tube wall profile is used as the input data until convergence occurs. The thermodynamic and transport properties of carbon dioxide are obtained following the methodology proposed by Span and Wagner [35]. Figure 3 shows the flowchart for the heat exchanger models.

Gas Cooler Model
The models of the heat exchangers are similar to each other, the main difference occurs in the external energy balance. In the gas cooler, the CO2 and water flow take place in counter-current in concentric tubes, as shown in Figure 4. In the gas cooler, the external energy balance applied between water and the copper tube is represented by Equations (10) and (11). The first one represents the energy gain over time for the water and the second one for the copper tube.

Gas Cooler Model
The models of the heat exchangers are similar to each other, the main difference occurs in the external energy balance. In the gas cooler, the CO 2 and water flow take place in counter-current in concentric tubes, as shown in Figure 4.

Gas Cooler Model
The models of the heat exchangers are similar to each other, the main difference occurs in the external energy balance. In the gas cooler, the CO2 and water flow take place in counter-current in concentric tubes, as shown in Figure 4. In the gas cooler, the external energy balance applied between water and the copper tube is represented by Equations (10) and (11). The first one represents the energy gain over time for the water and the second one for the copper tube. In the gas cooler, the external energy balance applied between water and the copper tube is represented by Equations (10) and (11). The first one represents the energy gain over time for the water and the second one for the copper tube.
where c p w is the specific heat at constant pressure and the subscripts w and o are the water and the internal tube. The external convective coefficient H o is calculated using the Dittus-Boelter equations [36].

Evaporator Model
The evaporator model is different from gas cooler model with respect to the external energy balance, as previously discussed. The evaporator/collector is a flat solar plate panel. It receives energy from solar radiation, from natural convection, as well as from the condensation of the water vapor present in the atmospheric air. Equation (12) presents the external energy balance for the evaporator.
where S and U L are solar radiation and overall heat transfer coefficient, respectively. F and w are fin efficiency and distance between tubes. The subscript e indicates external. For more details on external heat transfer in the evaporator, see Faria et al. [21]. Moreover, the evaporator is shown in Figure 5.
where is the specific heat at constant pressure and the subscripts and are the water and the internal tube. The external convective coefficient is calculated using the Dittus-Boelter equations [36].

Evaporator Model
The evaporator model is different from gas cooler model with respect to the external energy balance, as previously discussed. The evaporator/collector is a flat solar plate panel. It receives energy from solar radiation, from natural convection, as well as from the condensation of the water vapor present in the atmospheric air. Equation (12) presents the external energy balance for the evaporator.
where and are solar radiation and overall heat transfer coefficient, respectively. and are fin efficiency and distance between tubes. The subscript indicates external. For more details on external heat transfer in the evaporator, see Faria et al. [21]. Moreover, the evaporator is shown in Figure 5.

Heat Pump Model
The heat pump model was written in Fortran language, and each component of the model was developed in different subroutines. The integration of these subroutines is carried out in the main program. First, the outlet pressure and the evaporator superheat are estimated. After that, the gas cooler model generates the spatial profile using the data provided by the expansion device and compressor models. Using the output data of the gas cooler, the evaporator model subroutine starts. Similarly, the evaporator model generates the spatial profile of temperature and pressure using the data provided by the expansion device and compressor models. The evaporator model provides the new values of pressure and superheat, and they are compared with the previously estimated values. If there is no convergence, the process starts again, using the new values for the outlet pressure and superheat. Once the convergence between them has been reached, the model

Heat Pump Model
The heat pump model was written in Fortran language, and each component of the model was developed in different subroutines. The integration of these subroutines is carried out in the main program. First, the outlet pressure and the evaporator superheat are estimated. After that, the gas cooler model generates the spatial profile using the data provided by the expansion device and compressor models. Using the output data of the gas cooler, the evaporator model subroutine starts. Similarly, the evaporator model generates the spatial profile of temperature and pressure using the data provided by the expansion device and compressor models. The evaporator model provides the new values of pressure and superheat, and they are compared with the previously estimated values. If there is no convergence, the process starts again, using the new values for the outlet pressure and superheat. Once the convergence between them has been reached, the model moves on to the next time step.
The model runs until it reaches a steady state configuration. At this moment, the water outlet temperature in the gas cooler is compared to that suitable for bath (40 • C). If they do not converge, the model is restarted, and a correction in the water mass flow rate is made to obtain a new value of water outlet temperature. This methodology makes the water outlet temperature an input data instead of an output data. Furthermore, this convergence process also follows the Newton-Raphson methodology and represents a great differential of this model in relation to those that have been presented in the literature up to now. Figure 6 shows the flowchart for the heat pump distributed model.

2021, 13, x FOR PEER REVIEW 9 of 17
The model runs until it reaches a steady state configuration. At this moment, the water outlet temperature in the gas cooler is compared to that suitable for bath (40 °C). If they do not converge, the model is restarted, and a correction in the water mass flow rate is made to obtain a new value of water outlet temperature. This methodology makes the water outlet temperature an input data instead of an output data. Furthermore, this convergence process also follows the Newton-Raphson methodology and represents a great differential of this model in relation to those that have been presented in the literature up to now. Figure 6 shows the flowchart for the heat pump distributed model.

Volumetric Flow Rate Step-Model Validation
First of all, in order to validate the mathematical model, in transient and steady states, step change in the water volumetric flow rate was adopted. The model and experiment were running for 1500 s, while the water volumetric flow rate changed from 1.0 to 1.5 L/min (Step 1), as presented Figure 8. The results of the model and experiments are closer in the steady state as well as in the transient condition. Comparing the data, average and maximum deviations of 5% and 10% can be verified, respectively. Figure 9 shows the effects of the step in the following parameters: CO2 inlet ( 1 ) and outlet ( 2 ) temperatures at the gas cooler, CO2 inlet ( 1 ) and outlet ( 2 ) temperatures at the evaporator, and water outlet temperature ( 2 ).

Volumetric Flow Rate Step-Model Validation
First of all, in order to validate the mathematical model, in transient and steady states, step change in the water volumetric flow rate was adopted. The model and experiment were running for 1500 s, while the water volumetric flow rate changed from 1.0 to 1.5 L/min (Step 1), as presented Figure 8. The results of the model and experiments are closer in the steady state as well as in the transient condition. Comparing the data, average and maximum deviations of 5% and 10% can be verified, respectively. Figure 9 shows the effects of the step in the following parameters: CO 2 inlet (T gc 1 ) and outlet (T gc 2 ) temperatures at the gas cooler, CO 2 inlet (T ev 1 ) and outlet (T ev 2 ) temperatures at the evaporator, and water outlet temperature (T w 2 ).
In addition of the first step in the water volumetric flow rate at 1500 s, a second one at 2200 s (Step 2) occurred as presented Figure 9. In the first, the water volumetric flow rate increased, and in the second, it decreased. It should be noted that the variation of the water volumetric flow rate produced, as expected, an inverse effect on water outlet temperature. The deviations between model and experimental temperatures in Figure 9 are shown in Table 1.
were running for 1500 s, while the water volumetric flow rate changed from 1.0 to 1.5 L/min (Step 1), as presented Figure 8. The results of the model and experiments are closer in the steady state as well as in the transient condition. Comparing the data, average and maximum deviations of 5% and 10% can be verified, respectively. Figure 9 shows the effects of the step in the following parameters: CO2 inlet ( 1 ) and outlet ( 2 ) temperatures at the gas cooler, CO2 inlet ( 1 ) and outlet ( 2 ) temperatures at the evaporator, and water outlet temperature ( 2 ). In addition of the first step in the water volumetric flow rate at 1500 s, a second one at 2200 s (Step 2) occurred as presented Figure 9. In the first, the water volumetric flow rate increased, and in the second, it decreased. It should be noted that the variation of the water volumetric flow rate produced, as expected, an inverse effect on water outlet temperature. The deviations between model and experimental temperatures in Figure 9 are shown in Table 1.  Figure 10 presents the processing time for the model to perform the analysis presented in Figure 9. In the moments close to the steps (1 and 2) in the water volumetric flow rate and until the system reaches the steady state regime, the processing times were longer. In the steady state, the processing time is around 0.5 s. However, under the tran-Time (s) Figure 9. Variation of temperatures regarding the water volumetric flow rate steps. Table 1. Deviation between model and experiment temperatures presented in Figure 9.

Variable Average Deviation ( • C) Maximum Deviation ( • C)
T gc 1 0.9 2.1 T gc 2 1.7 1.9 T ev 1 1.4 1.5 T ev 2 1.9 2.7 T w 2 −0.3 1.3 Figure 10 presents the processing time for the model to perform the analysis presented in Figure 9. In the moments close to the steps (1 and 2) in the water volumetric flow rate and until the system reaches the steady state regime, the processing times were longer. In the steady state, the processing time is around 0.5 s. However, under the transient condition, the process time varied between 0.5 and 7.5 s. This is due to the greater computational effort required to represent the heat pump behavior in the transient regime. Finally, the processing time returns to the previous values when the system returns to the steady state condition.

Expansion Valve Opening Area
Step Figure 11 shows the variation of temperatures regarding the expansion valve opening area step. The step occurs at 1500 s, when the internal area of the expansion device is reduced. The reduction of the internal area of the expansion device produces a reduction in the CO2 mass flow rate. This leads to an increase in the superheat as previously discussed by Paulino et al. [22]. In other words, the compressor inlet temperature (evaporator outlet temperature, 2 ) increases. Consequently, the compressor outlet temperature (gas cooler inlet temperature, 1 ) increases as well. The gas cooler is a counter-current type. Therefore, if the CO2 gas cooler inlet temperature rises, the water outlet temperature must, of course, increase ( 2 ). If, however, the water mass flow rate and water inlet temperature are kept constant, then the gas cooler outlet temperature ( 2 ) maintained almost the same value. The energy received in the evaporator is the same during the simulation procedure, and it is

Expansion Valve Opening Area
Step Figure 11 shows the variation of temperatures regarding the expansion valve opening area step. The step occurs at 1500 s, when the internal area of the expansion device is reduced. The reduction of the internal area of the expansion device produces a reduction in the CO 2 mass flow rate. This leads to an increase in the superheat as previously discussed by Paulino et al. [22]. In other words, the compressor inlet temperature (evaporator outlet temperature, T ev 2 ) increases. Consequently, the compressor outlet temperature (gas cooler inlet temperature, T gc 1 ) increases as well. The gas cooler is a counter-current type. Therefore, if the CO 2 gas cooler inlet temperature rises, the water outlet temperature must, of course, increase (T w 2 ).

Expansion Valve Opening Area
Step Figure 11 shows the variation of temperatures regarding the expansion valve opening area step. The step occurs at 1500 s, when the internal area of the expansion device is reduced. The reduction of the internal area of the expansion device produces a reduction in the CO2 mass flow rate. This leads to an increase in the superheat as previously discussed by Paulino et al. [22]. In other words, the compressor inlet temperature (evaporator outlet temperature, 2 ) increases. Consequently, the compressor outlet temperature (gas cooler inlet temperature, 1 ) increases as well. The gas cooler is a counter-current type. Therefore, if the CO2 gas cooler inlet temperature rises, the water outlet temperature must, of course, increase ( 2 ). If, however, the water mass flow rate and water inlet temperature are kept constant, then the gas cooler outlet temperature ( 2 ) maintained almost the same value. The energy received in the evaporator is the same during the simulation procedure, and it is If, however, the water mass flow rate and water inlet temperature are kept constant, then the gas cooler outlet temperature (T gc 2 ) maintained almost the same value. The energy received in the evaporator is the same during the simulation procedure, and it is possible to notice that the evaporator inlet temperature (T ev 1 ) was almost unchanged due to the reduction of the internal area of expansion device.
The COP of the CO 2 DX-SAHP varies because of the adjustment of the expansion valve opening area, as shown in Figure 12. When the internal area of the expansion device changes, the COP also changes. This behavior is discussed by Rabelo et al. [24]. However, to analyze the COP trend (if it will increase or decrease for a specific operational point), it is important to know about the optimum high pressure.
Sustainability 2021, 13, x FOR PEER REVIEW 13 of 17 possible to notice that the evaporator inlet temperature ( 1 ) was almost unchanged due to the reduction of the internal area of expansion device. The COP of the CO2 DX-SAHP varies because of the adjustment of the expansion valve opening area, as shown in Figure 12. When the internal area of the expansion device changes, the COP also changes. This behavior is discussed by Rabelo et al. [24]. However, to analyze the COP trend (if it will increase or decrease for a specific operational point), it is important to know about the optimum high pressure.  Figures  11 and 12. However, as the step was higher, the effects on the evaporator outlet temperature, gas cooler inlet temperature and water outlet temperature were more pronounced. A higher step in the expansion valve opening area was executed. The effects of the greater reduction of the internal area are shown in Figure 13 for temperatures and in Figure 14 for the COP. The trends of all variables are the same as those presented in Figures 11 and 12. However, as the step was higher, the effects on the evaporator outlet temperature, gas cooler inlet temperature and water outlet temperature were more pronounced. possible to notice that the evaporator inlet temperature ( 1 ) was almost unchanged due to the reduction of the internal area of expansion device. The COP of the CO2 DX-SAHP varies because of the adjustment of the expansion valve opening area, as shown in Figure 12. When the internal area of the expansion device changes, the COP also changes. This behavior is discussed by Rabelo et al. [24]. However, to analyze the COP trend (if it will increase or decrease for a specific operational point), it is important to know about the optimum high pressure.  Figures  11 and 12. However, as the step was higher, the effects on the evaporator outlet temperature, gas cooler inlet temperature and water outlet temperature were more pronounced.

Conclusions
In this paper, a distributed mathematical model for a CO2 DX-SAHP is presented and validated experimentally. Based on the analysis of the dynamic behavior of the system, the following information was observed:

Conclusions
In this paper, a distributed mathematical model for a CO 2 DX-SAHP is presented and validated experimentally. Based on the analysis of the dynamic behavior of the system, the following information was observed: the comparison between the experimental results and those obtained with the model show small deviations for the analyzed temperatures and volumetric flow rate; -the computational processing time is longer for the transient condition. To show the variables with good representativeness, time steps of 25 and 5 s are appropriate for steady state and transient conditions, respectively; -an increase in the volumetric water flow rate produces a decrease in the water outlet temperature. On the other hand, a reduction in the volumetric water flow rate produces a rise in the water outlet temperature, as one would expect; -the reduction of the internal area of the expansion device shows an increase of inlet and outlet compressor temperatures as well as in the water outlet temperature; -the reduction of the internal area of the expansion device shows a variation in the CO 2 DX-SAHP coefficient of performance, as one would expect; -in the model presented in this paper, differently from what is usual in the literature for the CO 2 DX-SAHP, the water outlet temperature is an input data. This makes it possible to use the model under real operating conditions for which this variable is controlled.
At least, the presented model is a useful tool to analyze the CO 2 DX-SAHP in the transient and steady state regimes, demonstrating its behavior under different operating conditions.

Acknowledgments:
The authors appreciate the support of CAPES, CNPq, and FAPEMIG. We would also like to express our heartfelt gratitude to the late Ricardo Koury, at Federal University of Minas Gerais and our dear friend, who provided his exceptional insight and expertise while assisting us in our research.

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