Flexible dynamic model of PHEX for transient simulations in Matlab/Simulink using finite control volume method

1 Flexible dynamic model of PHEX for transient simulations in Matlab/Simulink using finite control volume method Erik Salazar-Herran, Koldobika Martin-Escudero, Luis A. del Portillo-Valdes, Ivan Flores-Abascal, Naiara Romero-Anton ENEDI Research group, Department of Thermal Engineering, University of the Basque Country (UPV/EHU), Plaza Torres Quevedo 1, 48013 Bilbao, Spain Email: [email protected] Telephone number: + (34) 94 601 7322 ABSTRACT In order to improve the energy efficiency and control of heat pump systems, it is necessary to develop dynamic models that accurately simulate their real performance. In addition, these models will help to carry out future works of research, such as new low carbon refrigerant testing. Physics-based models follow a set of physics laws that characterize the model as the most accurate, versatile and robust to simulate different heat pump systems. Taking into account the fact that the dynamics of the elements that regulate mass flow (compressors and valves) are much faster than the dynamics of the components that regulate heat transfer (heat exchangers), the model complexity usually resides in the latter. This paper provides a detailed explanation of the physics-based dynamic model in Matlab/Simulink using the finite-control volume approach applied to a refrigerant-to-liquid plate heat exchanger. Dynamic experimental tests were developed to validate the model under four possible situations: condenser and evaporator heat exchangers working in both counterand parallel-flow. In addition, an approximation of the number of finite control volumes required to reach a good accuracy, while maintaining a reasonable simulation time is presented. Simulation results show great accuracy when compared to experimental tests. It was proved by calculating the Normalized Residual Error, which is between 1.1 E-04 and 1.0 E-03 in all cases. It was also concluded that using twenty finite control volumes, there is good agreement between the accuracy of the results and the computational time. Keywords: Plate heat exchangers; Dynamic modeling; Experimental validation; Finite control volume approach; Reversible heat pump systems; Matlab/Simulink. This is the accepted manuscript of the article that appeared in final form in International Journal of Refrigeration 110 : 83-94 (2020) , which has been published in final form at https://doi.org/10.1016/j.ijrefrig.2019.11.003. © 2019 Elsevier Ltd and IIR. under CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) 2 Nomenclature Symbols A area [m2] c specific heat [kJ kg-1 K-1]] h enthalpy [kJ kg-1] ℎ󰇗 enthalpy time derivative [kJ kg-1 s-1] L width of the plates [m] m mass [kg] 𝑚󰇗 mass flow rate [kg s-1] N number of finite volumes [-] P pressure [kPa] 𝑃󰇗 pressure time derivative [kPa s-1] 𝑞󰇗 heat transfer rate [kW s-1] t time [s] T temperature [K] V volume [m3] x length of the plates in the flow direction [m] Greek letters 𝛼 convection heat transfer coefficient [kW m-2 K-1] 𝜕 partial derivative 𝜌 density [kg m-3] 𝜌󰇗 density time derivative [kg m-3 s-1] Subscripts ave average property c cross-sectional area h at constant enthalpy i control volume index in inlet out outlet p at constant pressure R refrigerant s heat transfer area S secondary fluid W wall/plate Acronyms COP coefficient of performance [-] FCV Finite Control Volume GSHP Ground-Source Heat Pump systems HEX Heat Exchangers HTC Heat Transfer Coefficient [kW m-2 K-1] MB Moving Boundary NRE Normalized Residual Error PHEX Plate Heat Exchanger RTF Real Time Factor 3 1. Introduction Energy sustainability is a global challenge that needs to be achieved in the coming years. Different steps have been taken towards this goal over the last few years, such as the Paris Agreement [1] and the Horizon 2020 program [2]. Within this challenge, Heat Pump systems are destined to be one of the systems that will help to achieve these internationally proposed goals. In fact, since 2009 in the European Union, a part of the energy produced by these systems can be considered renewable [3]. Ground source heat pump systems (GSHP) are a type of vapor compression system that, driven by electric energy, can heat or cool water or air throughout the year by transferring heat from or to a source, such as the subsoil or a water reservoir [4]. These systems have been proved to have better energy performance than air source heat pump systems [5]. Improving the performance and control of these systems is a key issue in making GSHP an increasingly competitive technology as compared to conventional heating and cooling systems, as well as for reaching the above-mentioned energetic sustainability. In this issue, heat pump manufacturers play a very important role in developing innovative GSHP with upgraded predictive control systems and better performance. Nevertheless, heat pump prototypes and experimental tests with new equipment require a high investment and this is time-consuming. Dynamic simulation models and software, which could reduce the costs of new products and improve their performance [6], would be beneficial all over the world, but such models need to be accurate, versatile and robust. In this line, physics-based models are an appropriate solution. These kinds of model use physical laws and partial derivatives to simulate the dynamic behavior of the system, providing great accuracy and the opportunity to simulate different components and configurations. Underwood [7] reviewed GSHP simulations, comparing models based on transfer functions and physics-based models, emphasizing the accuracy of the latter. Similarly, Rasmussen [8] compared physics-based and data-based models, highlighting the versatility of the former. The challenge when modeling a dynamic system usually lies in modeling the components with fast dynamics. In the case of heat pump systems, the heat exchangers have much slower dynamics than the components that regulate the mass flow rate, such as compressors or valves. For modeling heat exchangers (HEX) with physics-based models, different approaches and resolution methods can be found. Among others, two of the most commonly used approaches are the moving-boundary (MB) [9] and finite-control volume (FCV) [10]. While FCV divides the total length of the HEX into a finite number of volumes of equal size; the MB approach divides it into a reduced number of unequal and variable size finite volumes. The sizes of those finite volumes varies along the simulation and depends on a chosen property, such as the mean void fraction [11]. Some comparisons between both approaches have been presented over the last few years. Bendapudi et al. [12] compared them in a shell-and-tube heat exchanger of a centrifugal chiller system. The results were compared with experimental data and it was concluded that, during steady-states, MB was faster reaching equal accuracy. Nevertheless, during transient-states, the FCV formulation was found to be more robust. Herschel et al. [13] made a comparison between modeling an air-to-liquid HEX with both FVC and MB approaches, concluding that the computational time of MB was lower than 4 that of FCV, maintaining a good order of accuracy in both cases. Desideri et al. [14] compared MB and FCV approaches to simulate a plate heat exchanger (PHEX) concluding the similar good grade of accuracy of both model and highlighting the lower computational time of the MB approach. Similar conclusions were reached by Bonilla et al. [15] simulating an evaporator of a direct solar steam facility. Rodriguez et al. [16] compared four different dynamic modeling paradigms for an air-to-refrigerant evaporator. It was compared FCV approach with three MB methods: enthalpy, void fraction and density based MB. They highlighted the advantages of the FCV method due to the no need of switching triggers, thresholds or tolerances to reach an accurate results, although it is done at the expenses of a higher computational time. The MB approach has been widely used in air-to-refrigerant HEX dynamic modeling, as described in [8,17,18]. For instance, Ibrahim et al. [19] used it to predict the energy performance of an air source heat pump water heater. Similarly, MB approach has been used to simulate liquid-to-refrigerant HEX, such as [14,20]. Bell et al. [20] developed an MB model for counter flow HEX taking into account the fluid phase change, but was only able to predict stationary states. Chu et al. [21] were able to model a counter flow PHEX by using a MB-FCV coupling algorithm. Regarding FCV approach, Ozana et al. [22] implemented a dynamic model of a steam superheater HEX of an industrial boiler using the FCV approach. It was validated for both parallel and counter flow connections. Nevertheless, as for a superheated HEX, it was not necessary to model the phase change of the fluids. Bendapudi et al. [23] developed a centrifugal chiller system model using the FCV approach for the shell-and-tube heat exchanger, focusing on such aspects as mesh dependence, integrator order and stepsize. It was concluded that an increase in the number of finite control volumes increases the simulation accuracy up to a limit. Srihari et al. [24] analyzed the effect of a bad distribution of the flow into the plate heat exchanger (PHEX) channels for single-phase fluids using an FCV model. The transient performance of the airside heat exchanger of an air-water heat pump system under conditions of frosting was analyzed using a dynamic FCV model by Gao et al. [25]. Two of the most common modeling environments of physic-based models for HEX are Modelica [26] and Matlab/Simulink [27]. In Modelica many researches have been carried out, both with MB and FCV approach [14,15,28,29]. Moreover, some drawbacks of FCV such as chattering [30] have been studied concluding that that phenomenon would be reduced by using MB approach [31]. Regarding Matlab/Simulink, among others [32,33], one of the most known thermal systems simulation models are in the toolbox Thermosys [34]. It allows the user to simulate different refrigeration systems with different components such as air-torefrigerant condensers and evaporators or liquid-to-refrigerant condensers and evaporators. However, this last ones only in parallel flow configuration. These components are developed by using MB approach, which make difficult to develop a liquid-to-refrigerant HEX with counter-flow configuration. Taking all this into consideration and to the best of authors knowledge, a refrigerant-toliquid HEX dynamic model that includes phase-change phenomena (both condensation and evaporation) and that can be used indistinctly for both counter-flow and parallel-flow configurations in Matlab/Simulink environment cannot be found in the literature. This kind of model allows calculations to be carried out in GSHP, having the possibility to simulate switching mode operation. In addition, the extensive use of Matlab/Simulink software 5 must be taken into account in the world of thermal facilities manufacturers, especially for the development of control loops. The main goal of this document is to present a high accuracy dynamic model to simulate transitory changes in PHEX that can be connected in both parallel and counter-flow. Additionally, the model allows the same PHEX to work as a condenser or evaporator. The governing equations of the model and its matrix form are developed in detail. Moreover, an example of the formulation for a three finite control volumes PHEX for both counter and parallel-flow is presented. Finally, the validation of the model under different experimental situations is carried out and a parametric analysis to assess the influence of the number of FVC in the simulation speed is done. 2. Model development The model presented is a dynamic model based on physical equations using a distributed-parameter finite volume modeling method developed in a Matlab/Simulink environment. As we have already said, the developed model is able to simulate the behavior of a refrigerant-to-liquid PHEX under different situations. All these situations are grouped under the cases shown in Table 1. The governing equations presented and dealt with here are equal for all of the said cases. Nevertheless, when the equations are discretized into different finite control volumes and brought together in a matrix form, depending on whether the PHEX connections are in counter-flow or parallel-flow, the resulting matrices are different. As the direction of the refrigerant is inverted, what was the inlet port of the refrigerant becomes the outlet port and vice versa, changing the sign of some variables of the matrix. Additionally, the correlations for calculating the heat transfer coefficient (HTC) to obtain the heat transfer rates are subordinate to the working mode of the PHEX, and these are different for the condenser or evaporator working modes. Table 1. Cases analyzed for PHEX working modes and connections Case number PHEX connections and working mode Case 1 Counter-flow condenser Case 2 Parallel-flow condenser Case 3 Counter-flow evaporator Case 4 Parallel-flow evaporator 2.1. Model equations The equations used to describe the dynamic, physical behavior of the fluids inside the PHEX are now presented. The following assumptions have been taken into account for modeling the PHEX. • Pressure drops through the PHEX are negligible. • The maldistribution of the fluids is neglected. • Axial heat conduction is negligible, so the properties inside each finite volume are considered to be uniform. • The PHEX is supposed to be perfectly insulated from the surrounding environment. • The fluid flow is modeled as one-dimensional in the longitudinal direction. 6 The FCV approach is used to carry out the simulations. This model can capture the behavior and dynamics inside heat exchangers in great detail. The thermo-physical gradients and distributed parameters for each small discretization of the heat exchanger can be obtained. 2.1.1. Refrigerant equations The refrigerant fluid undergoes a phase change when it circulates through the PHEX. It can be condensed, going from superheated gas to sub-cooled liquid, or evaporated, going from two-phase fluid to superheated gas. Additionally, it can leave the PHEX as a two-phase fluid if the heat transferred between fluids is insufficient to conclude the condensation or evaporation processes. During transient states, the properties of the refrigerant fluid change over time. In order to model the changes in the fluid properties, energy and mass conservation equations are applied to the refrigerant fluid. The momentum equation is not applied because, as said before, the pressure drop through the PHEX is negligible. Eq. (1) and (2) present the mass and energy conservation equation applied to the refrigerant fluid respectively. The following steps are followed to transform both equations: • Apply 𝑐𝑅𝑇𝑅=ℎ𝑅 equality to calculate the enthalpy variation instead of temperature variation. • Integrate over the length of the control volume, which in this case is the total length of the heat exchanger. • Apply the Leibniz integral rule [36]. • Assume average values of fluid and plate properties inside the PHEX. • Solve time derivatives. • Take the enthalpy and the pressure as the independent variables to calculate the density time derivative. Now, the mass conservation equation and energy conservation equation applied to the refrigerant fluid has been modified in order to be implemented in a simulation code. where |𝜕𝜌 𝜕𝑃|ℎrepresents the refrigerant density variation with respect to the pressure variation at constant enthalpy and |𝜕𝜌 𝜕ℎ|𝑃 represents the refrigerant density variation with respect to the enthalpy variation at constant pressure [37]. 𝜕(𝜌𝑅𝐴𝑐) 𝜕𝑡 +𝜕(𝑚󰇗𝑅) 𝜕𝑥 =0 (1) 𝜕(𝜌𝑅𝑐𝑅𝐴𝑐𝑇𝑅) 𝜕𝑡 −𝜕(𝐴𝑐𝑃𝑅) 𝜕𝑡 +𝜕(𝑚󰇗𝑅𝑐𝑅𝑇𝑅) 𝜕𝑥 =𝛼𝑅𝐿(𝑇𝑅−𝑇𝑊) (2) (|𝜕𝜌 𝜕𝑃|ℎ𝑉)𝑃󰇗𝑅+(|𝜕𝜌 𝜕ℎ|𝑃𝑉)ℎ󰇗𝑅+𝑚󰇗𝑅,𝑜𝑢𝑡−𝑚󰇗𝑅,𝑖𝑛=0 (3) (|𝜕𝜌 𝜕𝑃|ℎℎ𝑅,𝑎𝑣𝑒−1)𝑉𝑃󰇗𝑅+(|𝜕𝜌 𝜕ℎ|𝑃ℎ𝑅,𝑎𝑣𝑒+𝜌𝑅,𝑎𝑣𝑒)𝑉ℎ󰇗𝑅+𝑚󰇗𝑅,𝑜𝑢𝑡ℎ𝑅,𝑜𝑢𝑡−𝑚󰇗𝑅,𝑖𝑛ℎ𝑅,𝑖𝑛 = 𝛼𝑅𝐴𝑠(𝑇𝑅−𝑇𝑊) (4) 7 Implementing Eq. (3) and Eq. (4) together in a simulation code, the variation of the refrigerant properties during transient states can be predicted. 2.1.2. Secondary fluid equations The velocity and the density of the secondary fluid, generally water, are assumed to be constant all along the heat exchanger. This means the equation of the mass conservation does not need to be applied for the secondary fluid. Eq. (5) presents the energy conservation equation for the secondary fluid. Unlike the refrigerant fluid, the secondary fluid will remain in a liquid phase throughout the PHEX. Therefore, the temperature variation can be directly calculated instead of enthalpy variation. Since the secondary fluid is assumed to be an incompressible fluid, the variation in the pressure is neglected. Therefore, the partial derivative of the pressure in terms of the time does not appear in Eq. (5). 𝜕(𝜌𝑆𝐴𝑐𝑐𝑆𝑇𝑆) 𝜕𝑡 +𝜕(𝑚󰇗𝑆𝑐𝑆𝑇𝑆) 𝜕𝑥 =𝛼𝑆𝐿(𝑇𝑊−𝑇𝑆) (5) Integrating each term of the equation over the length of the PHEX, assuming average properties inside the PHEX and solving the time derivatives, Eq (6) is obtained. 𝑉𝜌𝑆,𝑎𝑣𝑒𝑐𝑆,𝑎𝑣𝑒𝑇󰇗𝑆+𝑚󰇗𝑆𝑐𝑆,𝑎𝑣𝑒(𝑇𝑆,𝑜𝑢𝑡−𝑇𝑆,𝑖𝑛)= 𝛼𝑆𝐴𝑠(𝑇𝑊−𝑇𝑆) (6) By implementing Eq. (6) in a simulation code, the secondary fluid temperature variation in a HEX can be calculated with respect to the time. 2.1.3. Heat exchanger plate equations The temperature variation of the heat exchanger plate is described by applying the energy conservation equation. It is presented in Eq. (7) and no transformations are needed. 𝑐𝑊𝑚𝑊𝑇󰇗𝑊=𝛼𝑅𝐴𝑠(𝑇𝑅−𝑇𝑊)−𝛼𝑆𝐴𝑠(𝑇𝑆−𝑇𝑊) (7) The temperature of the plate depends on the mass, the thermal properties of the plate material and the heat transferred to or from both refrigerant and secondary fluids. 2.2. Model implementation Once the necessary equations have been transformed into temporal derivative equations, they are then discretized for use with the FCV approach. For an 𝑁 finite number of control volumes, Eqs. (3), (4), (6) and (7) are transformed into Eqs. (8) - (11). 𝑉𝑐𝑣[(|𝜕𝜌 𝜕𝑃|ℎ)𝑖ℎ𝑅,𝑖−1]𝑃󰇗𝑅+𝑉𝑐𝑣[(|𝜕𝜌 𝜕ℎ|𝑃)𝑖ℎ𝑅,𝑖+𝜌𝑅,𝑖]ℎ󰇗𝑅,𝑖+𝑚󰇗𝑅,𝑜𝑢𝑡,𝑖ℎ𝑅,𝑜𝑢𝑡,𝑖 −𝑚󰇗𝑅,𝑖𝑛,𝑖ℎ𝑅,𝑖𝑛,𝑖= 𝛼𝑅,𝑖𝐴𝑠,𝑐𝑣(𝑇𝑅,𝑖−𝑇𝑊,𝑖) (8) 𝑉𝑐𝑣(|𝜕𝜌 𝜕𝑃|ℎ)𝑖𝑃󰇗𝑅+𝑉𝑐𝑣(|𝜕𝜌 𝜕ℎ|𝑃)𝑖ℎ󰇗𝑅,𝑖+𝑚󰇗𝑅,𝑜𝑢𝑡,𝑖−𝑚󰇗𝑅,𝑖𝑛,𝑖=0 (9) 𝑉𝑐𝑣𝜌𝑆,𝑖𝑐𝑆,𝑖𝑇󰇗𝑆,𝑖+𝑚󰇗𝑆,𝑖𝑐𝑆,𝑖(𝑇𝑆,𝑜𝑢𝑡,𝑖−𝑇𝑆,𝑖𝑛,𝑖)= 𝛼𝑆,𝑖𝐴𝑠,𝑐𝑣(𝑇𝑊,𝑖−𝑇𝑆,𝑖) (10) 8 𝑐𝑊,𝑐𝑣𝑚𝑊,𝑐𝑣𝑇󰇗𝑊,𝑖=𝛼𝑅,𝑖𝐴𝑠,𝑐𝑣(𝑇𝑅,𝑖−𝑇𝑊,𝑖)−𝛼𝑆,𝑖𝐴𝑠,𝑐𝑣(𝑇𝑆,𝑖−𝑇𝑊,𝑖) (11) Then, in order to implement them in a Matlab code, they are reorganized into a matrix equation. [ 𝑉𝑐𝑣[(|𝜕𝜌 𝜕𝑃|ℎ)𝑖ℎ𝑅,𝑖−1]𝑉𝑐𝑣[(|𝜕𝜌 𝜕ℎ|𝑃)𝑖ℎ𝑅,𝑖+𝜌𝑅,𝑖]0 0 𝑉𝑐𝑣(|𝜕𝜌 𝜕𝑃|ℎ)𝑖𝑉𝑐𝑣(|𝜕𝜌 𝜕ℎ|𝑃)𝑖0 0 0 0 𝑉𝑐𝑣𝜌𝑆,𝑖𝑐𝑝,𝑆,𝑖 0 00 0 𝑐𝑊𝑚𝑊,𝑐𝑣 ] × [ 𝑃󰇗𝑅 ℎ󰇗𝑅,𝑖 𝑇󰇗𝑆,𝑖 𝑇󰇗𝑊,𝑖 ] = [ 𝑚󰇗𝑅,𝑖𝑛,𝑖ℎ𝑅,𝑖𝑛,𝑖−𝑚󰇗𝑅,𝑜𝑢𝑡,𝑖ℎ𝑅,𝑜𝑢𝑡,𝑖+2𝑞󰇗𝑅,𝑖 𝑚󰇗𝑅,𝑖𝑛,𝑖−𝑚󰇗𝑅,𝑜𝑢𝑡,𝑖 −𝑚󰇗𝑆,𝑖𝑐𝑆,𝑖(𝑇𝑆,𝑜𝑢𝑡,𝑖−𝑇𝑆,𝑖𝑛,𝑖)+2𝑞󰇗𝑆,𝑖 𝑞󰇗𝑅,𝑖−𝑞󰇗𝑆,𝑖 ] (12) Both heat transfer rates (𝑞󰇗𝑅 and 𝑞󰇗𝑆) are multiplied by two because, in the actual PHEX, every refrigerant or secondary fluid are in contact with two intermediate plates, except for the boundary secondary fluid channels. Empirical correlations are used in the calculation of the heat transfer rates. For the one-phase convection HTC, the DittusBoelter correlation is used [38]. In addition, the correlations presented in Han et al. [39,40] are used for the two-phase convection HTC. In order to solve Eq. (12), the inlet and outlet mass flow rates of the refrigerant fluid must be known. When a steady-state is calculated, both are equal. Nevertheless, during a transient-state, expansion and compression devices regulate the refrigerant mass flow and the one can be different from the other. On the other hand, the inlet enthalpy of the refrigerant fluid and the inlet temperature and mass flow rate of the secondary fluid must also be known. As will be shown later, when Eq. (12) is extended to an arbitrary number of control volumes and the refrigerant inlet and outlet mass flow rates are known, the intermediate mass flow rates must be calculated at each time step. These intermediate refrigerant mass flow rates will be moved to the second term of the matrix equation, forming another state vector. In the aggregate, for each time step, 4𝑁 derivative states and integrations must be solved. Nevertheless, only the refrigerant energy and mass equations have to be solved together. The energy equation of plates and secondary fluid can be removed from the matrix and solved separately, thus reducing the computational time. On the other hand, at the time of applying the presented formulation to a fixed number of FCV, counter-flow and parallel-flow connections must be differentiated. Figure 1 and Figure 2 show a graph description of a PHEX discretized in N number of FCV for counterflow and parallel-flow connections, respectively. 9 Figure 1. Finite control volume counter-flow PHEX model. Figure 2. Finite control volume parallel-flow PHEX model. As can be seen, the outlets and inlets of the refrigerant fluid are inverted when the connection of the PHEX is in counter-flow or parallel-flow. The inversion of the refrigerant fluid allows the model to be used in reversible heat pump systems which can heat and cool-down water, depending on the user demand. 2.3. Three FCV example cases Eq. (12) will now be extended to a three control volume formulation for both a counterflow and a parallel-flow connection. The differences in the formulation between both connections are shown below. Eqs. (13) and (14) show the matrix equations that solve the refrigerant energy and mass equations for counter-flow and parallel-flow PHEX, respectively. As the energy equation of secondary fluid and plates can be solved separately and are equal regardless of the connection of the PHEX, they have been removed from the main matrix. Eq. (15) shows the secondary fluid energy matrix equation and Eq. (16) shows the intermediate plate energy matrix equation for the three control volume example. Secondary fluid Refrigerant fluid Plate Control volumen 1 Control volumen 2 Control volumen n-1 Control volumen n . . . 𝑚󰇗𝑆𝑚󰇗𝑆 𝑚󰇗𝑅,𝑜𝑢𝑡 𝑚󰇗𝑅,𝑖𝑛 𝑞󰇗𝑅, 𝑞󰇗𝑆, 𝑇𝑆,𝑜𝑢𝑡 ℎ𝑅,𝑜𝑢𝑡 𝑇𝑆,𝑖𝑛 ℎ𝑅,𝑖𝑛 𝑇𝑆, 𝑇𝑆,𝑛 𝑛 𝑇𝑆, 𝑇𝑆, 𝑇𝑆,𝑛 𝑇𝑆,𝑛 ℎ𝑅, ℎ𝑅,𝑛 𝑛 ℎ𝑅, ℎ𝑅, ℎ𝑅,𝑛 ℎ󰇗𝑅,𝑛 𝑞󰇗𝑅, 𝑞󰇗𝑆, 𝑞󰇗𝑅,𝑛 𝑞󰇗𝑆,𝑛 𝑞󰇗𝑅,𝑛 𝑞󰇗𝑆,𝑛 𝑞󰇗𝑆, 𝑞󰇗𝑆, 𝑞󰇗𝑆,𝑛 𝑞󰇗𝑆,𝑛 𝑞󰇗𝑅, 𝑞󰇗𝑅, 𝑞󰇗𝑅,𝑛 𝑞󰇗𝑅,𝑛 Plate Plate 𝑚󰇗𝑅, 𝑚󰇗𝑅,𝑛 𝑛 Secondary fluid Refrigerant fluid Plate Control volumen 1 Control volumen 2 Control volumen n-1 Control volumen n . . . 𝑚󰇗𝑆𝑚󰇗𝑆 𝑚󰇗𝑅,𝑜𝑢𝑡 𝑚󰇗𝑅,𝑖𝑛 𝑞󰇗𝑅, 𝑞󰇗𝑆, 𝑇𝑆,𝑜𝑢𝑡 ℎ𝑅,𝑜𝑢𝑡 𝑇𝑆,𝑖𝑛 ℎ𝑅,𝑖𝑛 𝑇𝑆, 𝑇𝑆,𝑛 𝑛 𝑇𝑆, 𝑇𝑆, 𝑇𝑆,𝑛 𝑇𝑆,𝑛 ℎ𝑅, ℎ𝑅,𝑛 𝑛 ℎ𝑅, ℎ𝑅, ℎ𝑅,𝑛 ℎ𝑅,𝑛 𝑞󰇗𝑅, 𝑞󰇗𝑆, 𝑞󰇗𝑅,𝑛 𝑞󰇗𝑆,𝑛 𝑞󰇗𝑅,𝑛 𝑞󰇗𝑆,𝑛 𝑞󰇗𝑆, 𝑞󰇗𝑆, 𝑞󰇗𝑆,𝑛 𝑞󰇗𝑆,𝑛 𝑞󰇗𝑅, 𝑞󰇗𝑅, 𝑞󰇗𝑅,𝑛 𝑞󰇗𝑅,𝑛 Plate Plate 𝑚󰇗𝑅, 𝑚󰇗𝑅,𝑛 𝑛 16 (a) (b) (c) Figure 8. Parallel-flow condenser (Case 2) results of refrigerant pressure (a), water temperatures (b) and temperature distribution at time 800 seconds (c). Figure 8(c) depicts the temperature distribution in the parallel-flow condenser. As can be seen, the refrigerant fluid just reaches the sub-cooled saturation line in the final part of the PHEX. On the other hand, Figure 9(c) and Figure 10(c) represent the temperature distribution when the PHEX works as an evaporator for the counter-flow and the parallel-flow, respectively. As can be seen, the refrigerant fluid goes into the PHEX as a two-phase fluid. Inside the PHEX, it is evaporated and superheated. (a) (b) 2400 2500 2600 2700 2800 2900 0 200 400 600 800 1000 Condensation pressure [kPa] Time [s] Test data Model results 25 30 35 40 45 0 200 400 600 800 1000 Water temperature [ºC] Time [s] Outlet (model) Outlet (test) Inlet 30 40 50 60 70 80 010 20 30 40 50 Temperature [ºC] PHEX length [cm] Refrigerant Intermediate plate Water 700 800 900 1000 1100 0 200 400 600 800 1000 Evaporation pressure [kPa] Time [s] Test data Model results 0 5 10 15 20 0 200 400 600 800 1000 Water temperature [ºC] Time [s] Outlet (model) Outlet (test) Inlet 17 (c) Figure 9. Counter-flow evaporator (Case 3) results of refrigerant pressure (a), water temperatures (b) and temperature distribution at time 800 seconds (c). (a) (b) (c) Figure 10. Parallel-flow evaporator (Case 4) results of refrigerant pressure (a), water temperatures (b) and temperature distribution at time 800 seconds (c). Regarding the PHEX connections, in both Case 2 and Case 3, the water inside the PHEX is cooled down by around 5 ºC. Nevertheless, in Case 2, the refrigerant fluid goes out from the PHEX with a superheating of around 5ºC, while in Case 3, the superheating is of approximately 2 ºC. Such a difference is related to the PHEX connection, the heat transfer being higher in a counter-flow connection than in a parallel-flow connection. In Table 6, the maximum and mean numerical differences between the model results and the data from the tests for refrigerant pressure and water outlet temperature are presented. Table 6. Maximum and mean numerical differences between model results and tests data for the four studied cases. Refrigerant pressure [kPa] Water outlet temperature [ºC] 5 7 9 11 13 15 17 010 20 30 40 50 Temperature [ºC] PHEX length [cm] Refrigerant Intermediate plate Water 900 950 1000 1050 1100 1150 0 200 400 600 800 1000 Evaporation pressure [kPa] Time [s] Test data Model results 0 5 10 15 20 0 200 400 600 800 1000 Water temperature [ºC] Time [s] Outlet (model) Outlet (test) Inlet 4 6 8 10 12 14 010 20 30 40 50 Temperature [ºC] PHEX length [cm] Refrigerant Intermediate plate Water 18 Maximum Mean Maximum Mean Case 1 51.5 20.3 1.1 0.3 Case 2 52.6 30.4 1.3 0.6 Case 3 46.9 12.4 0.5 0.1 Case 4 51.9 21.8 0.8 0.2 The results show that the maximum numerical difference of the refrigerant pressure is around 50 kPa for any studied case. Regarding the comparison for the water outlet temperature, the greatest difference is found in Case 2. Here, the maximum difference increases up to 1.3 ºC. Additionally, the NRE for the refrigerant pressures and water outlet temperatures for the four cases are shown in Table 7. As can be seen, in all the cases studied, the NRE is smaller than the accepted value of 0.05. Table 7. Normalized Residual Error for the four studied cases. Refrigerant pressure Water outlet temperature Case 1 1.1 E-04 1.2 E-04 Case 2 1.6 E-04 2.5 E-04 Case 3 2.2 E-04 2.5 E-04 Case 4 5.8 E-04 1.0 E-03 4.3. FVC number vs simulation speed Finally, the simulation speed and accuracy have been studied by increasing the number of FCV used to simulate the described model from five to one hundred. The counter-flow condenser (Case 1) has been simulated for 20 minutes and the Real Time Factor (RTF) of the simulation has been calculated, dividing the computational time by the simulated time. Figure 11 (a) shows the differences in the condensing pressure calculation using a different number of finite volumes. Although the simulation time was 1200 seconds, in the image only the period where the transient occurs has been shown, so the differences can be better perceived. Figure 11 (b) shows the RTF for different numbers of FCV. As can be appreciated in Figure 11 (a), the shape of the system dynamic is well followed, even with a low number of FCV. Nevertheless, the numerical results do not accurately fit with a low number of FCV. With ten FCV, the results can be improved. The numerical results are accurate enough once the number of FCV is equal to or higher than twenty. Regarding the computational time, the RTF remains low up to 20 FCV and increases greatly with 100 FCV. These results agree with the conclusions given in [42], where it is stated that at least 15 FCV are required to achieve good accuracy. 19 (a) (b) Figure 11. Condensing pressure with different numbers of finite control volumes (a) and Real Time Factor for different numbers of finite control volumes (b). 5. Conclusions In order to improve the control systems of reversible heat pump systems, models that can predict the dynamic behavior and transient states of these systems are needed. In this task, physics-based models are the most accurate ones. In this paper, a physics-based model using the finite volume control approach to simulate a PHEX is presented. The model has been developed in Matlab/Simulink environment. This kind of heat exchanger is commonly used in residential reversible ground source heat pump systems. As it is based on physical equations, this model can be used to model different types of PHEX, achieving a good level of accuracy. The model has been validated with experimental dynamic tests for a PHEX working as condenser and as evaporator. Moreover, for each working mode, it has been simulated under counter-flow and parallel-flow arrangements. The accuracy of the model has been proved by calculating the NRE of the simulated results with respect to the test data of the refrigerant pressure and the water outlet temperature. The minimum value of the NRE was 1.1 E-04, while the maximum was 1.0 E-03, which is lower than the accepted value for this kind of model. It has also been demonstrated that discretizing the model in twenty finite volumes allows a high degree of accuracy to be reached, while maintaining the computational time low. This paper provides a useful tool for heat pump systems and PHEX design as it allows the performance of a PHEX to be compared, varying the working mode of the PHEX, the fluid connections, the PHEX size, the number of PHEX plates, the working fluids, etc. Additionally, as the model can simulate a PHEX working as a condenser and as an evaporator, it can be used to simulate the behavior of reversible heat pump systems. REFERENCES [1] United Nations. Paris Agreement 2017;21:1–7. [2] European Union. Horizon 2020 EU Research and Innovation programme n.d. https://ec.europa.eu/programmes/horizon2020/en/what-horizon-2020. [3] European Union. Directive 2009/28/Ce. Off Diary Eur Union 2009;140:16–62. [4] Lucia U, Simonetti M, Chiesa G, Grisolia G. Ground-source pump system for heating and cooling: Review and thermodynamic approach. Renew Sustain Energy Rev 2016;70:0–1. doi:10.1016/j.rser.2016.11.268. 1850 1950 2050 2150 2250 2350 350 450 550 650 750 850 Condensing pressure [kPa] Time [s] 5 FCV 10 FCV 20 FCV 50 FCV 100 FCV 0.0044 0.0053 0.0079 0.0177 0.0738 0.00 0.02 0.04 0.06 0.08 5 FCV 10 FCV 20 FCV 50 FCV 100 FCV Real Time Factor [-] 20 [5] Salazar-Herran E, Martin-Escudero K, Del Portillo-Valdes LA, Flores-Abascal I, Picallo-Perez A. Residential Heat Pumps as Renewable Energy. 8th Eur. Conf. Energy Effic. Sustain. Archit. Plan., 2017, p. 163–70. [6] Romero-Anton N, Martin-Escudero K, Portillo-Valdés LA, Gómez-Elvira I, SalazarHerran E. Improvement of auxiliary BI-DRUM boiler operation by dynamic simulation. Energy 2018;148. doi:10.1016/j.energy.2018.01.160. [7] Underwood CP. Heat pump modelling. Elsevier Ltd; 2016. doi:10.1016/B978-008-100311-4.00014-5. [8] Rasmussen BP. Dynamic modeling for vapor compression systems — Part I : Literature review Review Article Dynamic modeling for vapor compression systems — Part I : Literature review 2012;9669. doi:10.1080/10789669.2011.582916. [9] Grald EW, MacArthur JW. A moving-boundary formulation for modeling timedependent two-phase flows. Int J Heat Fluid Flow 1992;13:266–72. doi:10.1016/0142-727X(92)90040-G. [10] MacArthur J. Transient heat pump behaviour: a theoretical investigation. Int J Refrig 1984;7:123–32. doi:10.1016/0140-7007(84)90025-2. [11] Wedekind GL, Bhatt BL, Beck BT. A system mean void fraction model for predicting various transient phenomena associated with two-phase evaporating and condensing flows. Int J Multiph Flow 1978;4:97–114. doi:10.1016/03019322(78)90029-0. [12] Bendapudi S, Braun JE, Groll EA. A comparison of moving-boundary and finitevolume formulations for transients in centrifugal chillers. Int J Refrig 2008;31:1437–52. doi:10.1016/J.IJREFRIG.2008.03.006. [13] Pangborn H, Alleyne AG, Wu N. A comparison between finite volume and switched moving boundary approaches for dynamic vapor compression system modeling. Int J Refrig 2015;53:101–14. doi:10.1016/j.ijrefrig.2015.01.009. [14] Desideri A, Dechesne B, Wronski J, van den Broek M, Gusev S, Lemort V, et al. Comparison of Moving Boundary and Finite-Volume Heat Exchanger Models in the Modelica Language. Energies 2016;9:339. doi:10.3390/en9050339. [15] Bonilla J, Dormido S, Cellier FE. Switching moving boundary models for twophase flow evaporators and condensers. Commun Nonlinear Sci Numer Simul 2015;20:743–68. doi:10.1016/J.CNSNS.2014.06.035. [16] Rodriguez E, Rasmussen B. A comparison of modeling paradigms for dynamic evaporator simulations with variable fluid phases. Appl Therm Eng 2017;112:1326–42. doi:10.1016/J.APPLTHERMALENG.2016.10.131. [17] Rasmussen BP, Shenoy B, Rasmussen BP, Shenoy B. Dynamic modeling for vapor compression systems — Part II : Simulation tutorial Review Article Dynamic modeling for vapor compression systems — Part II : Simulation tutorial 2012;9669. doi:10.1080/10789669.2011.582917. [18] Qiao H, Laughman CR, Aute V, Radermacher R. An advanced switching moving boundary heat exchanger model with pressure drop. Int J Refrig 2016;65:154–71. doi:10.1016/J.IJREFRIG.2016.01.026. [19] Ibrahim O, Fardoun F, Louahlia-gualous H. Air source heat pump water heater : Dynamic modeling , optimal energy management and mini-tubes condensers 2014;64:1102–16. doi:10.1016/j.energy.2013.11.017. 21 [20] Bell IH, Quoilin S, Georges E, Braun JE, Groll EA, Horton WT, et al. A generalized moving-boundary algorithm to predict the heat transfer rate of counterflow heat exchangers for any phase configuration. Appl Therm Eng 2015;79:192–201. doi:10.1016/J.APPLTHERMALENG.2014.12.028. [21] Chu Z, Zhang W. Moving-boundary and finite volume coupling algorithm for heat exchanger with fluid phase change. Int J Heat Mass Transf 2019;131:313–28. doi:10.1016/J.IJHEATMASSTRANSFER.2018.11.066. [22] Ozana S, Pies M. Using Simulink S-Functions with Finite Difference Method Applied for Heat Exchangers n.d.:210–5. [23] Bendapudi S, Braun JE, Groll EA. Dynamic model of a centrifugal chiller system - Model development, numerical study, and validation. ASHRAE Trans., vol. 111 PART 1, 2005, p. 132–48. [24] Srihari N, Prabhakara Rao B, Sunden B, Das SK. Transient response of plate heat exchangers considering effect of flow maldistribution. Int J Heat Mass Transf 2005;48:3231–43. doi:10.1016/j.ijheatmasstransfer.2005.02.032. [25] Gao T, Gong J. Modeling the airside dynamic behavior of a heat exchanger under frosting conditions. J Mech Sci Technol 2011;25:2719–28. doi:10.1007/s12206011-0615-5. [26] Modelica n.d. [27] MathWorks, Moler C. Matlab/Simulink 1984. [28] Bonilla J, Yebra LJ, Dormido S, Cellier FE. Object-Oriented Modeling of Switching Moving Boundary Models for Two-phase Flow Evaporators. IFAC Proc Vol 2012;45:1069–74. doi:10.3182/20120215-3-AT-3016.00189. [29] Sangi R, Jahangiri P, Müller D. A combined moving boundary and discretized approach for dynamic modeling and simulation of geothermal heat pump systems. Therm Sci Eng Prog 2019. doi:10.1016/j.tsep.2018.11.015. [30] Quoilin S, Bell I, Desideri A, Dewallef P, Lemort V. Methods to Increase the Robustness of Finite-Volume Flow Models in Thermodynamic Systems. Energies 2014;7:1621–40. doi:10.3390/en7031621. [31] Bonilla J, Yebra LJ, Dormido S. Chattering in dynamic mathematical two-phase flow models. Appl Math Model 2012;36:2067–81. doi:10.1016/J.APM.2011.08.013. [32] Underwood CP. Fuzzy multivariable control of domestic heat pumps. Appl Therm Eng 2015;90:957–69. doi:10.1016/j.applthermaleng.2015.07.068. [33] Al-Dawery SK, Alrahawi AM, Al-Zobai KM. Dynamic modeling and control of plate heat exchanger. Int J Heat Mass Transf 2012;55:6873–80. doi:10.1016/j.ijheatmasstransfer.2012.06.094. [34] Alleyne A, Rasmussen B, Jeddeloh J. Thermosys Toolbox n.d. [35] Rasmussen, Bryan P; Alleyne A. Dynamic modeling and advanced control of air conditioning and refrigeration systems. University of Illinois at UrbanaChampaign, 2005. [36] Protter MH, Morrey CB. Differentiation under the Integral Sign. Improper Integrals. The Gamma Function. Intermed. Calc., New York, NY: Springer New York; 1985, p. 421–53. doi:10.1007/978-1-4612-1086-3_8. 22 [37] Thorade M, Saadat A. Partial derivatives of thermodynamic state properties for dynamic simulation. Environ Earth Sci 2013;70:3497–503. doi:10.1007/s12665013-2394-z. [38] Dittus FW, Boelter LMK. Heat transfer in automobile radiators of the tubular type. Int Commun Heat Mass Transf 1985;12:3–22. doi:10.1016/0735-1933(85)90003X. [39] Han DH, Lee KJ, Kim YH. The characteristics of brazed plate heat exchangers with different chevron angles. J Korean Phys Soc 2003;43:66–73. doi:10.1016/S1359-4311(03)00061-9. [40] Han DH, Lee KJ, Kim YH. Experiments on the characteristics of evaporation of R410A in brazed plate heat exchangers with different geometric configuration. Appl Therm Eng 2003;23:1209–25. doi:10.1016/S1359-4311(03)00061-9. [41] Salazar-Herran E, Martin-Escudero K, Alleyne AG, del Portillo-Valdes LA, Romero-Anton N. Numerical model for liquid-to-liquid heat pumps implementing switching mode. Appl Therm Eng 2019;160:114054. doi:10.1016/j.applthermaleng.2019.114054. [42] Bendapudi S, Braun JE, Groll EA. Dynamic model of a centrifugal chiller system - Model development, numerical study, and validation. ASHRAE Trans., vol. 111 PART 1, 2005, p. 132–48.