FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 16, N o 2, 2018, pp. 249 - 259 https://doi.org/10.22190/FUME180424019N © 2018 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper DYNAMICAL BEHAVIOR OF A HEAT PUMP COAXIAL EVAPORATOR CONSIDERING THE PHASE BORDER’S IMPACT ON CONVERGENCE UDC 621.1 Arpad Nyers 1 , Zoltan Pek 2 , Jozsef Nyers 2,3 1 University of Pécs, Faculty of Engineering and Information Technology, Hungary 2 Obuda University, Doctoral School of Applied Informatics and Applied Mathematics, Hungary 3 Szent István University, Mechanical Engineering PhD School, Hungary Abstract. Using a dynamical mathematical model, we investigated transient behavior of a water-water heat pump’s evaporator. The model consists of time and space dependent partial differential equations of water, pipe wall and refrigerant. Mathematically the thermal expansion valve (TEV) and compressor are described with lumped parameters and represent the boundary conditions. During the numerical solution of this system of equations the problem emerged of divergence of solutions. It was determined that the cause of the divergence solution was in the location of phase change of the refrigerant. The aim of this paper is, firstly, to display and propose a new approach to the elimination of divergence. In addition, it examines dynamic behavior of the heat pumps’ coaxial evaporator. Key Words: Heat Pump, Coaxial Evaporator, Dynamical Behavior, Phase Border, Convergence, Numerical Simulation 1. INTRODUCTION The evaporator is probably the most important part of the heat pump since it represents the place where heat is recovered; this, to a great extent, influences energy efficiency of the system. The process efficiency is significantly influenced by the amount and quality of injected refrigerant into the evaporator. The injection process must be controlled and for the design of the controller the transient behavior of the evaporator must be known. We tried to take into account all those physical phenomena which significantly influence the process. Received April 24, 2018 / Accepted June 10, 2018 Corresponding author: Jozsef Nyers Obuda University, Doctoral School of Applied Informatics and Applied Mathematics, Budapest, Hungary E-mail: jnyers1@gmail.com 250 A. NYERS, Z. PEK, J. NYERS Several dynamic models of heat pumps have been developed but few of them dealt with the heat pump’s coaxial evaporator. MacArthur and Grald [1] presented a model of vapor- compression heat pumps. The heat exchangers were modeled with detailed distributed formulations while the expansion device was modeled as a simple fixed orifice. Pavkovic et al. [2] dealt with a fully distributed dynamic numerical model of a shell and a tube type of the refrigerant evaporator, suitable for control and design purposes. Avoidance of the fraction model is used to account for a slip effect. The solution is achieved by the finite volume method. Fu et al. [3] presented a dynamic model of air-to-water dual-mode heat pump with screw compressor having four-step capacities. The dynamic responses of adding additional compressor capacity in a step-wise manner were studied. Kima et al. [4] presented a dynamic model of the water heater system driven by the heat pump; the finite volume method was applied to describe the heat exchangers while the lumped parameter models were used to analyze the compressor and the hot water reservoir where dynamic simulations were carried out for various reservoir sizes. Rasmussen and Alleyne [5] developed an air to air heat pump system using a moving boundary-lumped parameter approach to the heat exchangers. Đorđević et al. [6] experimentally studied the heat transfer of the Archimedean spiral coil made of a transversely corrugated tube as a heat absorber of the parabolic dish solar concentrator. The results showed enhancement of the heat transfer rate in comparison with the spirally-coil smooth tubes, up to 240%. Nyers et al. presented in their earlier works [7-10] the structure and functioning of the heat pump’s evaporator, its dynamic and discrete mathematical model as well as the numerical method for solving the model as well. Szamarszkij et al. [11] described in detail the methods for solution of gas dynamics equations. Kajtar and Kassai [12-15] dealt with computerized simulation of the energy consumption. Some further works [16-21] dealt with the same topics as in this paper, including mathematical modeling, numerical procedures, heat pump and heat pump heating-cooling system. In this study we propose a new numerical method for solving the divergence of solution at the phase change border which can also be applied to the most complicated model. The obtained results presented by graphs and diagrams visibly prove the solution improvement. The improvement is achieved by means of the proposed method for calculating the time step. 2. PHYSICAL MODEL OF EVAPORATOR Constructively, the considered evaporator is a coaxial tube heat exchanger. From thermal aspect, in the evaporator the heat transfer is performed between the refrigerant and the well water. In the observed case, the refrigerant Freon R134a flows inside the evaporator pipes while the cooled mass, the well water, flows inside the evaporator’s shell. In the evaporator, the parallel pipes are connected with the baffles. The baffles are placed perpendicularly to the pipes. The pipes bundle at a distance of about 150 mm in the direction of the tube axis. Water flows between the baffles approximately as a sinusoid. The evaporator works most efficiently when the full length of parallel pipes is filled with vapor-liquid mix phase of the refrigerant. The quantity of the evaporated refrigerant depends on the compressor capacity, the temperature and the mass flow rate of well water, respectively. The thermo-expansion (TEX) valve doses an adequate quantity and quality of the refrigerant into the evaporator. This valve uses the sensors for monitoring Dynamical Behavior of a Heat Pump Coaxial Evaporator Considering the Phase Border’s Impact ... 251 temperature and pressure of the outlet refrigerant from the evaporator. Based on the measured temperature and pressure, the valve realizes optimal dosing of the refrigerant. In case of a refrigerant quantity deficit, the liquid flows for a shorter length of the evaporator’s pipes while in the remaining part of the evaporator flows the superheated vapor. The heat transfer of the refrigerant's vapor is multiple times lower than of the vapor-liquid mix phase. The refrigerant vapor in the dry superheated section superheats. The refrigerant over-dosing means that the liquid phase cannot evaporate and some quantity of the liquid leaves the evaporator. This phenomenon most reduces the refrigeration effect, i.e. the performance of the heat pumps' evaporator. The non-evaporated liquid evaporates in the compressor and the consequences of the hydro hummer can happen in the compressor. In the correct operation mode the superheated, dry vapor flows out of the evaporator. The level of the superheating is up to 4 °C. The evaporator is connected to the compressor. The compressor sucks out the superheated vapor from the evaporator. The compressor using mechanical work compresses the vapor which primarily increases the temperature on the adequate level. The compression is unfortunately accompanied with an intensive increasing of the vapor pressure as well. Fig. 1 The physical model of heat pump's evaporator with the system parameters Fig. 2 Photo of the heat pump with coaxial evaporator and shell-tube condenser 3. DYNAMIC MATHEMATIC MODEL OF EVAPORATOR In this paper, both the mathematical model and its numerical solution represent an improvement over those reported in our earlier work [7, 10]. The mathematical model is built up with following assumptions:  evaporator consists of horizontal bundled pipes, the refrigerant is distributed uniformly between them, the heat resistance of the pipe wall neglected,  refrigerating flow inside the pipes is homogeneous and one-dimensional,  heat resistance of the pipe wall is neglected, and,  axial heat conduction is neglected everywhere. 252 A. NYERS, Z. PEK, J. NYERS The mathematical model built on these assumptions contains the equations of the refrigerant, water, pipe wall, thermal expansion valve TEV and compressor. The mathematical description of the refrigerant’s dynamic behavior is the most complicated part of the model. It contains the conservation equations of mass, impulse and energy which are valid in the evaporation region, at the phase change border and in the superheated region as well. Conservation equations of mass ( ) 0 , p ρ w 1 , ρ t z v         (1) where p is the pressure, w the velocity, ρ the density, v the specific volume, z the space increment while t denotes the time. Conservation equations of impulse ( ) ( ) ( ) 0, 2 ρ w ρ w p f x t z           (2) where f(x) is the friction factor in flowing refrigerant as a function of vapor quality, x. The equation of energy conservation yields: ( / 2 ) ( ( / 2 )) ( ) 0 2 2 ρ w ρ i p w ρ w ρ i q x , t z                (3) where i is the refrigerant enthalpy and q(x) is the heat flux as a function of vapor quality. Equation (3) is written in a divergent form since this is a solid basis for the construction of a reasonable discretization of the equations: it assures the validity of discrete conservation laws [11]. The above equations are complemented by equations of mixture and vapor state, equations of heat transfer, friction coefficients and heat flux. The equations of state differ in the evaporator and the superheated region (the phase change point will be determined by the vapor quality): Evaporation region: xo <= x <= xmax Vapor-liquid equilibrium equation as a function of pressure, p, and temperature, T: 1( ) 0f p,T , (4) Specific enthalpy of refrigerant mixture: ( )t g ti i x i i ,    (5) where the subscripting t and g denote that the quantity is related to liquid and vapor, respectively. Specific volume of refrigerant mixture: ( )t g tv v x v v ,    (6) Dynamical Behavior of a Heat Pump Coaxial Evaporator Considering the Phase Border’s Impact ... 253 The heat flux between the pipe wall (index c) and the refrigerant (index f) as a function of vapor quality: ( ) ( ) ( )u c fq x α x 4/d T T    (7) where du is the inner diameter. The coefficients of heat transfer and friction as a function of vapor quality are denoted by (x) and f(x), respectively. Phase change border x = xmax =1 Vapor-liquid equilibrium equation: 1( ) 0f p,T , (8) Equation of the superheated vapor state: 2 ( ) 0f p,T, v , (9) Fictive equation, (due to the number of equations): 3 ( ) 0maxf x, x , (10) Again, the coefficients of heat transfer and friction as a function of vapor quality are denoted by (x) and f(x), respectively. Superheated region: x = xmax =1 Equation of superheated vapor state: 2 ( ) 0f p,T, v , (11) Equation of superheated vapor enthalpy: 4 ( ) 0f p,T, v,i , (12) Fictive equation (due to the number of equations): 5 ( 1) 0f x, , (13) Also, in this case, the coefficients of heat transfer and friction as a function of vapor quality are denoted by (x) and f(x), respectively. The water through the pipe wall transfers heat to the refrigerant. The dynamical heat transient behavior of the water and the wall are described by the energy conservation law: ( ) 0 v v v v v c T T w C T T t z           (14) where the subscripting v means that the quantity is related to water, while the sign (+) is used in the case of concurrent flow and the sign (-) in the case of countercurrent flow. ( ) ( ) 0 c cv v c cf c f T C T T C T T t          (15) where Ccv and Ccf are the coefficients of interaction between the tube and water and the tube and the refrigerant, respectively. 254 A. NYERS, Z. PEK, J. NYERS The refrigerant flows through the thermal expansion valve TEV and the compressor. These components of the system are described using algebraic equations with lumped parameters, as follows: Thermal expansion valve’s equation, TEV based on Darcy – Weissbach formula: 2 ( ) ( ) ( ) 0conp p C t w     (16) where the subscripting con means that the quantity is related to the compressor and C(t) is the throttle coefficient of the TEV valve as a function of time. The compressor`s steady-state equation for the piston mechanism reads: 1 /( ) 0 n con c h c c p C A w V n p C            (17) where n is the exponent of polytrophic, 1≤ n ≤ 1.26, Cc is the coefficient of clearance volume, A is the surface area of piston, Vh the work volume of compressor and nc the number of revolutions of compressor. In practice, the piston speed is high so that the vapor heat loss to the environment is low; therefore, the process of compression can be approximately regarded as an adiabatic, consequently n=1.26. 4. NUMERICAL SOLUTION METHOD The above mathematical model can be solved only approximately, using a numerical approach. To obtain the corresponding formula, the conservation laws are integrated over small sub regions (cells) of the z, t-domain in which the solution is to be determined. The area integrals can be transformed to line integrals using Green’s integral theorem; the line integrals are approximated by simple formula like the trapezium and rectangle rule. Weighting parameters are introduced for stability reasons. In this way a system of coupled nonlinear, algebraic equations arises for the solution of which the Newton method is applied. The linear algebraic system to be solved in every step of the Newton method has a distinct structure; it is of a block-tri diagonal form. This is a result of the discretization described above and of the fact that we write the mass conservation equation not for the density but for the specific volume and then discretize it over the cell {xi ≤ x ≤ xi+1, ti ≤ t ≤ ti+1}, whereas the remaining equations are discretized over the cell {xi-1 ≤ x ≤ xi, ti ≤ t ≤ ti+1}. The block-tri diagonal linear equations are solved by block-Gauss elimination. Special treatment needs the point of phase change where the vapor quality reaches the value x=1. At this point the approximation is written down not on a rectangular cell but on a triangular one. By appropriate choice of the time step we manage the phase change point to move from time level to time level by at most one spatial interval only (a shift of this point by more than one interval has an adverse influence). The solution of this problem is obtained as follows: in the first iteration on a given time level, using old time step ∆to, the number m of ∆z-intervals is found by which the phase change point moves. Dynamical Behavior of a Heat Pump Coaxial Evaporator Considering the Phase Border’s Impact ... 255 Fig. 3 Phase boundary displacement after ot time The solution of this problem is obtained as follows: in the first iteration on a given time level, using old time step ∆to, the number m of ∆z-intervals is found by which the phase change point moves. The number of space intervals m, zm z1 t t ,k-km o t        (18) m t t o    (19) where k is the number of increments and kt is the number of increment at phase boundary when the vapor quality x=1. ∆t is taken as a new time step for the next iteration. When the iteration on a fixed time level has converged, ∆to is taken as the first time step on the next time level since in our problem, after a period of considerable changes in the solution, convergence to a new stationary state follows. In that latter r time interval, the phase change point settles down and time step (much smaller than ∆t is not necessary for reaching sufficient accuracy) but would mean a waste of computing time. 5. RESULTS AND DISCUSSION Using real-life physical data we simulated the evaporation in a pipe, taking into account the different processes listed above. The evaporator consists of a 13x10/8 mm copper pipe of 10 m length. The compressor supplies 10.m f  kg/s R134a refrigerant whereas the mass flow of water is 50.mv  kg/s. The system’s transient behavior is affected by a jump in the water temperature for - 4°C (from 14°C to 6°C) or for +6°C (from 14°C to 20°C). The transient processes induced by these temperature changes are well depicted by our computations, see Figs. 4-7. Also visible is the advantage of the described time step choice over an arbitrarily specified time step, especially if looking at the computed refrigerant temperatures. In the case of an arbitrary time step the temperature at the phase change point attains too low or high values and oscillates in the superheated region, significantly differing from reasonable values; see the Figs. 4 (-4°C) and 6 (+6°C) displaying the Tf values. This non-physical effect practically disappears if using the proposed time step choice; see Figs. 5 (-4°C) and 7 (+6°C). 256 A. NYERS, Z. PEK, J. NYERS Fig. 4 Refrigerant temperature changes in time and pipe length affected by water temperature jump of -4 °C (from 14 °C to 10 °C). In the superheated section no natural oscillations are visible because of the wrong time step. Fig. 5 Refrigerant temperature changes in time and the pipe length is affected by a water temperature jump of -4 °C (from 14 °C to 10 °C). In the superheated section the unnatural oscillations disappear as a result of the appropriate time step. Dynamical Behavior of a Heat Pump Coaxial Evaporator Considering the Phase Border’s Impact ... 257 Fig. 6 Refrigerant temperature changes in time and the pipe length is affected by water temperature jump of +6 °C (from 14 °C to 20 °C). In the superheated section unnatural oscillations are visible because of the wrong time step. Fig. 7 Refrigerant temperature changes in time and the pipe length is affected by water temperature jump of +6 °C (from 14 °C to 20 °C). In the superheated section unnatural oscillations disappear as a result of the appropriate time step. 6. CONCLUSION In this section some remarks about the determination of the phase change point are given. Without an accurate calculation of the location of this point no numerical convergence could be obtained in the Newton method combined with the block-Gauss elimination for linear systems. In our model the phase change point was determined by vapor quality x=1. The choice of the mathematical model was also influenced by our experiments with solving this problem; finally, as reported above, we have found it necessary to separate the evaporation region and 258 A. NYERS, Z. PEK, J. NYERS the superheated region and to write at the phase change point the equilibrium equation as well as the state equation of the superheated gas together with the other leading equations. Using the recommended algorithm it is possible to determine the appropriate time step to ensure the convergence of solution. In the given case the convergence time step was 0.5s. As mentioned, the appropriate time step was obtained iteratively by using Eq. (19). The advantage of the recommended method is that it is very simple and efficient. The disadvantage is that the appropriate time step cannot be calculated in an exact mathematical manner. Rather, it can be defined only numerically in real time using the recommended iterative procedure. The diagrams in Figs. 4-5 show the evaporator dynamic behavior for -4 °C unit jump of water temperature. As a response: the refrigerant input temperature decreases from 10.3 °C to 7.2 °C, the output vapor temperature barely varies whilst the evaporation length significantly reduces from 7.3m to 5.7m. The time duration of the transient process is 10 s. In the case of water temperature increase of +6 °C (Fig 6-7), the refrigerant's input temperature increases from 10.5 °C to 12.6 °C, the vapor output temperature barely varies while the evaporation length reduces from 8.6 m to 7.5 m. REFERENCES 1. MacArthur, J.W., Grald, E.W., 1989, Unsteady compressible two-phase flow model for predicting cyclic heat pump performance and a comparison with experimental data, International Journal of Refrigeration, 12(1), pp. 29-41. 2. Pavkovic, B., Vilicic, I., Medica, V., 2000, Numerical Simulation of Transients in Shell and Tube Refrigerant vaporator, Advanced Computational Methods in Heat Transfer VI, ed. B. Sunden & C. A. Brebbia, WIT Press: Southampton, pp. 457-466. 3. Fu, L., Ding, G. and Zhang, C., 2003, Dynamic simulation of air-to-water dual-mode heat pump with screw compressor, Applied Thermal Engineering, 23, pp. 1629-1645. 4. Kima, M., Kim, M.S. and Chung, J.D., 2004, Transient thermal behavior of a water heater system driven by a heat pump, International Journal of Refrigeration, 27, pp. 415-421. 5. Rasmussen, B. P. and Alleyne, A. G., 2006, Dynamic Modeling and Advanced Control of Air Conditioning and Refrigeration Systems, Air Conditioning and Refrigeration Center University of Illinois TR-244. 6. Đorđević, M., Stefanović, V., Vukić, M., Mančić, M., 2017, Experimental investigation of the convective heat transfer in a spirally coiled corrugated tube with radiant heating , Facta Universitatis- Series Mechanical Engineering, 15( 3), pp. 495-506. 7. Nyers, J., Stoyan, G., 1994, A Dynamical Model Adequate for Controlling the Evaporator of Heat Pump. International Journal of Refrigeration, 17(2), pp.101-108. 8. Nyers, J., Pek, Z., 2014, Mathematical Model of Heat Pumps' Coaxial Evaporator with Distributed Parameters, Acta Polytechnica Hungarica, 11(10), pp. 41-54. 9. Nyers, J., 2016, COP and Economic Analysis of the Heat Recovery from Waste Water using Heat Pumps, International J. Acta Polytechnica Hungarica, 13(5), pp. 135-154. 10. Nyers, J., Nyers, A., 2013, Hydraulic analysis of heat pump's heating circuit using mathematical model, 9 th ICCC International Conference, Proceedings-USB, Tihany, Hungary, pp. 349-353,. 11. Szamarszkij, A., Popov A., Yu, P., 1980, Difference methods for solution of gas dynamics equations, 2nd ed., Nauka, Moscow. 12. Kajtár, L., Kassai, M., Bánhidi, L., 2011, Computerized simulation of the energy consumption of air handling units, Energy and Buildings, (45), pp. 54-59. 13. Kassai, M., 2016, Energy demand and consumption investigation of a single family house, International symposium, Subotica, Serbia, Proceedings EXPRES 2016, pp. 30-35. 14. Kajtár, L., Kassai, M., 2010, A new calculation procedure to analyze the energy consumption of air handling units, Periodica polytechnica-mechanical engineering, 54(1), pp. 21-26. Dynamical Behavior of a Heat Pump Coaxial Evaporator Considering the Phase Border’s Impact ... 259 15. Kassai, M., Ge, G., Simonson, J. C., 2016, Dehumidification performance investigation of run-around membrane energy exchanger system, Thermal Science, 20(6), pp. 1927-1938. 16. Wu, W., Skye, H. M., Domanski, P. A., 2018, Selecting HVAC systems to achieve comfortable and cost- effective residential net-zero energy buildings, Applied Energy, 212, pp. 577-591. 17. Kalmár, F., 2016, Interrelation between glazing and summer operative temperature in buildings, International Review of Applied Sciences and Engineering, 7(1), pp. 51–60. 18. Szabó, G. L., Kalmár, F., 2017, Investigation of subjective and objective thermal comfort in the case of ceiling and wall cooling systems, International Review of Applied Sciences and Engineering, 8(2), pp. 153-162. 19. Takács, J., 2015, Enhance of the efficiency of exploitation of geothermal energy, International symposium, Subotica, Serbia, Proceedings EXPRES 2015, pp. 46-49. 20. Poos, T., Szabo, V., 2017, Determination of entry length of a fluidized bed dryer using volumetric heat transfer coefficient, International Review of Applied Sciences and Engineering, 8(1), pp. 57-65. 21. Odry, A., Kecskés, I., Burkus, E., Odry, P., 2017, Protective Fuzzy Control of a Two-Wheeled Mobile Pendulum Robot: Design and Optimization, WSEAS transactions on systems and control, 12(32), pp. 297-306.