THE MATHEMATICAL MODEL AND NUMERICAL SIMULATION OF THE HEAT PUMP SYSTEM

Abstract

In this paper water-water heat pump is analyzed using a computer simulation. The simulation program is based upon steady state mathematical models of the refrigeration system components. The numerical simulation based off on mass, momentum and energy conservation is developed. A fully distributed parameter method was used for analyzing the condenser and evaporator. Lumped parameter method was applied to simulate the compressor and expansion device. The Runge-Kutta method is used to solve the differential and algebraic equations. [PUBLICATION ABSTRACT]

Headnote

ABSTRACT: In this paper water-water heat pump is analyzed using a computer simulation. The simulation program is based upon steady state mathematical models of the refrigeration system components. The numerical simulation based off on mass, momentum and energy conservation is developed. A fully distributed parameter method was used for analyzing the condenser and evaporator. Lumped parameter method was applied to simulate the compressor and expansion device. The Runge-Kutta method is used to solve the differential and algebraic equations.

KEYWORDS:

(ProQuest: ... denotes formulae omitted.)

INTRODUCTION

The European Commission accepts a proposal package [1] of future indicators at the beginning of 2008. The aim of this is to decrease the rapidly increasing greenhouse gas emissions based on the 2007 data. Furthermore to increase the renewable energy sources in the total energy consumption in the proportion of 20% by the year 2020.

The utilization of renewable energy sources is influenced by several factors. Besides the natural environment, the economic conditions are also major factors affecting the case of renewable energies. The fossil fuel prices and conditions of other energy costs are significantly determined by the demand for renewable as well as the amount of state aid and government fiscal policy application. According to the current standpoint of the "Új Széchenyi Terv" - An Economic Development Plan [2] - all conditions exist and could be provided so that the energy efficient method of heat pump technology of heating and cooling buildings could be spread in Hungary.

The energy efficiency improvement of the heat pump to improve the operation quality makes it unavoidable to strive for the heat pump operation. Nevertheless, the processes taking place in it should be as accurate as possible to describe the underlying physical and mathematical model development and refinement.

The experimental investigation of any refrigeration system is usually very complicated, mainly due to the financial costs and the large number of variables involved. The use of numerical models can reduce the costs and also facilitate understanding the phenomena related to the problem. Refrigeration systems models are divided in two broad classes: steady-state models and transient models. Dynamic phenomena have been widely observed and studied by several investigators including MacArthur [3], Chi and Didion [4], P. Welsby [5], Wang [6], Yasuda et al. [7] and Nyers and Gisbert [8].

Many researchers dealt with the description of the steady state behavior of a vapor compression refrigeration system such as Koury et al. [9] proposed a model for a refrigeration system with distributed parameter model for heat exchanger, Jong Won Choi et al. [10], Belman et al. [11], Yang Zhao et al. [12] and Garbai et al. [13].

The main objective of this work is to present the mathematical model and numerical simulation and to simulate the steady state behavior of a vapor compression refrigeration system. A distributed parameter method was used for analyzing the evaporators and condenser. Lumped parameter method was applied to simulate the compressor and expansion device. The simulations were performed with the refrigerants R134a.

DESCRIPTION OF PHYSICAL SYSTEM

The cooling circuit, i.e. the heat pump, is in the centre of the heat pump heating - cooling system. A heat pump consists of the four main components: evaporator, compressor, condenser and expansion valves. The refrigerant is the working fluid of heat pumps.

In the evaporator, the refrigerant takes over the heat from the low-temperature primary fluid by vaporization. In the evaporator, the refrigerant is superheated and vapor is sucked in by compressor. With the invested mechanical work, it is brought higher level of energy. On the discharge side of the compressor, the now hot and highly pressurized vapor is cooled in the condenser.

In the condenser, vapor provides the heat to the secondary fluid and then the refrigeration condenses. Condensation in horizontal tubes may involve partial or total condensation of the vapor. Depending on the application, the inlet vapor may be superheated, equal to 1.0 or below 1.0.

Superheated vapor enters the horizontal tube, which has a temperature below the saturation temperature of the vapor. The flow at this point in the tube is single-phase vapor flow. After the vapor cools and becomes saturated, condensation starts to occur on the inner wall of the tube.

The flow pattern at the beginning of condensation is annular, because the velocity of the vapor is much higher than that of the condensate. The dominant force in the annular flow regime is shear stress at the liquid-vapor interface, with gravity playing a less important role. As condensation continues, the velocity of the vapor phase decreases and the dominant force shifts from shear force at the interface to gravitational force. Liquid accumulates at the bottom of the tube, while condensation takes place mainly at the top portion of the tube where the liquid film is thin.

Near the outlet of the horizontal tube, the vapor quality reduces to zero and the flow in the tube becomes single-phase liquid flow.

The condensed refrigerant then passes through a pressure-lowering expansion device.

The expansion device is to reduce the pressure and to regulate the refrigerant mass flow rate. The widely utilized expansion device is the thermostatic expansion valves. The thermostatic expansion is a valve for controlling the refrigerant flow by a sensor bulb placed in the evaporator discharge line and hence controls the mass rate by the degree of superheat. The low pressure, refrigerant leaving the expansion device enters the evaporator, in which the refrigerant absorbs heat and boils. The refrigerant then returns to the compressor and the cycle is repeated.

The observed heat exchangers are counter-cross flow, shell and tube type. Tubes are made of copper and have a staggered layout. In the current case, the refrigerant R134a flows through in finned tube bundle of heat exchangers, while primary and secondary fluid flow in the shell across the tube bundle. In the shell side of heat exchangers baffles were applied primarily used for supporting the tubes and for inducing cross flow over the tubes, resulting in improved heat transfer performance. The baffles are perpendicular to the tube bundle. The examined heat exchangers have segmental baffles. The segment baffles force the (main stream) primary and secondary fluid to change direction. The rate of change in direction is affected by the baffles insertion length. In this case, the compression occurs on the principle of displacement reciprocating compressor. While the throttle is isenthalpic and occurs with variable cross-section of expansion valve.

MATHEMATICAL MODEL OF REFRIGERATION SYSTEM

The simulation program is based upon steady state mathematical models of the components of the refrigeration circuit including the compressor, heat exchangers, thermostatic expansion valve. What follows the model corresponding to each component of the refrigeration circuit will be detailed. It must be noted that mathematical models of the component such as the compressor is based upon manufacturer data.

Heat exchangers (Shell and tube)

The evaporator and the condenser were assumed to have the same size and type. The refrigerant flows through the inside of the inner tube and the primary/secondary fluid flows through the annulus. The evaporator is divided into two zones (evaporation and superheating zones), and the condenser is also divided into two zones, such as superheating, condensation zones. Each zone the equations of continuity, energy and momentum are described. The general forms of the governing equations are shown in Eqs. (1), (2), and (3).

Continuity equation:

... (1)

Momentum equation:

... (2)

Energy equation:

... (3)

Pipe wall energy conservation:

... (4)

Water (secondary fluid) energy conservation:

... (5)

Heat transfer

Since phase change of the refrigerant occurs as it flows in the system, the correlations for heat transfer coefficient for the refrigerant should be selected considering its phase.

Kandlikar (1998) [14] correlation was applied to obtain the two-phase heat transfer coefficient for the evaporator while for the condensation heat transfer coefficient for the condenser, Shah correlation was used.

For the single phase regions, vapor phase in the evaporator and condenser, Dittus -Boelter correlation (Eqs. 10) was selected. For calculating water (primary/secondary) side heat transfer coefficient correlation proposed by J. Taborek (1983) [15] was used.

Two-phase flow boiling heat transfer

Kandlikar correlation was used for the prediction of the heat transfer coefficient in flow boiling. The final correlation consists of two sets of constants. One for the convective evaporation dominated regime and the other for the nucletae boiling dominated regime. The correlation is:

... (6)

and the constants are given in the table below.

Heat transfer coefficients of single-phase refrigerant in eqs. (7) was calculated by the Petukhov and Popov (1963), [16].

Interval of applicability: 0.5 < Prlo < 2000 and 104 < Relo < 5 ?106

... (7)

The correlation is calculated twice using each set of constants and the greater of the two values is used as the heat transfer coefficient.

... (8)

Two-phase condensation heat transfer

For the two phase regions, in the condenser, Shah correlation [17] was selected to proposed heat transfer coefficient. The Shah correlation is a modified version of Dittus-Boelter single-phase heat transfer correlation (10). The two-phase model the reduced pressure refrigerant takes into account.

... (9)

In Eqs. (9) the single phase heat transfer coefficients are determined by the Dittus-Boelter correlation, where the reduced pressure is:

... (10)

Single-phase shell side heat transfer

J. Taborek [18] correlation proposes the heat transfer coefficient of the secondary fluid. The single-phase correlation is taking the evaporator and condenser tube bundle and baffle plates of the spatial distribution into account.

... (11)

* Baffle cut correction factor: ...

* Baffle leakage correction factor: ...

* Bundle bypass correction factor: ...

* Unequal baffle spacing correction factor: ... * Unequal baffle spacing correction factor: ...

* Ideal tube bank heat transfer coefficient: ...

Pressure drop

* The total pressure drop of the refrigerant

The total pressure drop of a fluid is due to the variation of kinetic and potential energy of the refrigerant and that is due to the friction on the flow channel wall. The total pressure drop is the sum of the momentum pressure drop and the frictional pressure drop. The two-phase friction factors are estimated by using the Friedel correlation [19] while the Blasius correlation was adopted for the single phase flow.

* Momentum pressure drop

In the present work the change in momentum caused by the pressure drop is taken into account.

... (12)

* Frictional pressure drop

The calculation of the frictional pressure drop is based on a method from Friedel Esq. (13)-(17) The two-phase pressure gradient is given by:

... (13)

* The resulting two-phase multiplier is now calculated as follows:

... (14)

* The dimensionless parameters E, F, and H are evaluated as follows:

...

* The two-phase mixture density is calculated as follows:

... (15)

* The corresponding single-phase pressure gradients are given by:

... (16)

* For smooth tubes the friction factor, f, can be determined as following Eq. (17)

... (17)

* Single-phase shell side pressure drop

For calculating water side pressure drop, correlation proposed by J. Tagore (1983) was used. The pressure drop for shell-side flow is equal to the sum of inlet nozzle pressure drop, the bundle pressure drop and the outlet of the nozzle pressure drop.

The bundle pressure drop is equal to the sum of the cross-flow pressure drops ?pC, the window pressure drops ?pw , and the two end zone pressure drops (first and last baffle compartments) ?pE.

... (18)

* The pressure drop in all the central baffle compartments is:

... (19)

* The pressure drop and mass velocity in all Nb window zones for turbulent flow is:

... (20)

* The pressure drop in the two end zone of the tube bundle is:

... (21)

Void fraction

The void fraction was computed by the Rouhani-Axelsson correlation. Rouhani-Axelsson [20] (1970) presented a drift flux model based on that of Zuber and Findlay [21] (1965) in the form of

... (22)

The weighted mean drift flux velocity and the distribution parameter based on a horizontal configuration are in agreement with experimental data proposed by VDI-Warmeatlas [22] (1993) using the additional expressions:

... (23)

After substitution into equation, this results in the expression:

... (24)

This form of the Rouhani-Axelsson void fraction model was proposed for ge neral use by Steiner [22] (1993) irrespective of flow pattern for horizontal flows. This method takes the effects of concentration distribution and relative velocity into account and includes the important effects of mass velocity, viscosity and surface tension.

Compressor model

A compressor model is developed from the data provided by the manufacturer. Performance characteristic are given using three equations fitting the compressor energy consumption (25), cooling capacity (26) and the mass flow rate of refrigeration (27) as a function of the evaporation temperature Te and the condensation temperature Tc for a given suction temperature:

... (25)

... (26)

... (27)

where Pi, Qi and Mi are constant coefficients.

* Thermostatic expansion valve model

Thermostatic expansion valve is the valve that controls the refrigerants mass flow rate by sensing the degree of suction vapor superheat temperature. The enthalpy is assumed to be constant. The refrigerant mass flow is calculated by the following equation

... (28)

where C is the characteristic constant of the valve.

* Thermodynamics equations of state refrigerant

Refrigerant R134a is used in this study as a working fluid. A computer program was developed to compute various properties of R134a. The P-R state equation is used here [23]

... (29)

The predicted transport properties were compared with the values given in Solkane R134a database. The transport properties of secondary fluid include viscosity, thermal conductivity; surface tension and diffusion coefficient are generally calculated from the existing correlations.

NUMERICAL CONSIDERATIONS

The governing equations for the refrigerant, the water and the pipe wall and the auxiliary equations are solved numerically by casting them in finite difference form using the control volume approach, as shown in figure 4. Because the steady state terms are considered, the mass flux of each cell should be the same.

Because of the fact that the refrigerant and the water flow are in opposite directions, the simultaneous solution of the five equations (1) - (5) is very complicated. Hence, Eqs. (1)-(3) were solved separately by the fourth order Runge-Kutta method, along the z axis, while the temperature profile was determined from Eqs. (4) and (5) by the finite difference implicit method. This sequence was repeated until achieving convergence. The second part of the model consists of explicit and implicit non-linear algebraic equations. For solving the implicit non-linear algebraic equations the iterative procedure with Newton linearization is used.

The water side and pipe wall energy equations (4), (5) are alternatively used to determine the pipe wall temperature and the leaving water temperature of each control volume where the former is assumed at the beginning of the iteration. The boundary conditions for the water side are the mass flow rate and the inlet water temperature.

A flow chart outlines the model structure and the order of computation in Fig. 5. Prior to cycle calculation, initial conditions must be specified together with a degree; this condition prevents performing the refrigerant inventory which usually leads to high uncertainties.

Sizing the compressor, geometrical parameters heat exchangers and transport properties of refrigerant and water is a preliminary step.

Initialization values are necessary to run the model: inlet water conditions (flow rate, supply temperature) are required to estimate initial evaporation and condensation pressures.

In the evaporator and condenser models, the enthalpies and mass flux values at inlet and outlet, were determined from the models of the expansion device and compressor. Starting from arbitrary values of evaporation and condensation pressures and utilizing the initial conditions of the problem, the spatial profiles of the refrigerant enthalpy, mass flux and density can be calculated at each point of the heat exchanger.

The transition between the single-phase zone to the two-phase-zone and vice versa are determined by computing the refrigerant vapor quality along the entire length of the heat exchanger.

The superheat is compared to the required value. Evaporation pressure is adjusted according to the difference. It seems that the condenser parameters mainly affect the discharge pressure and the compressor mass flow rate, while the evaporator parameters mainly affect the refrigerant mass flow rate and the heat transfer rates which are also noticed by Harms [25].

RESULT AND DISCUSSION

In this study, the steady-state behavior simulation of a refrigeration system of the water-water type was presented. During the work the changes of parameters of the heat pump system were analyzed. The following investigations occurred in dependence of changes of the mass flow rate of the cooled fluid in the evaporator shell side. The numerical results of the obtained simulation are presented in a visually graphical form, as shown in figure 5 - 19.

The simulations were carried out under the following circumstances:

* Refrigerant R134a

* Number and diameter of tubes in the heat exchagers:5 x Φ 8 x 1 [mm]

* Diameter of shell side of heat exchangers: 36 [mm]

* The heat exchangers were discredited into control volumes numbered from i = 1 to n=30

* At a constant temperature of superheated in evaporation ?TT = 4 [K]

* At a constant temperature of the input primary fluid Tpv,in =13 [oC]

* At a constant temperature of the input secondary fluid Tsv,in =27 [oC]

* At a constant mass flow rate of secondary fluid m? s =0.28 ???kgs ???

In the simulation the following was changed:

* mass flow rate of primary fluid mp = 0.28, 0.42, 0.55, 0.69 |^kg/]

CONCLUSIONS

Based on the simulations the following can be stated:

* The accuracy of the heat pump mathematical model mainly depends on the applied correlation for two phase boiling and condensation heat transfer, both primary and secondary fluid. The investigated five models of two phase boiling heat transfer correlation the Kandlikar's model has proven to be suitable.[24 At the beginning and at the end of the evaporation heat transfer coefficient corresponds to the expectation. The character, i.e. tendency to change the Kandlikar`s heat transfer coefficient along the evaporator is real. It resembles a parabola, as show in figure 9.

* The refrigerant flow in the tube of condenser at the phase change heat transfer coefficient values are changed by leaps and bounds, as show in figure 17.

* The friction pressure drop and void fraction correlation affects the accuracy of the mathematical model of the heat pump at the lower level of extension. The pressure drop refrigerant in the evaporator is of greater amount than in the condenser. In the condenser, the pressure drop of the refrigerant is practically constant. Shown in figure 14.

* The increase of the heating water temperature profiles is not even, the heat transfer of the two-phase refrigerant is higher, because of the latent heat than the single phase heat transfer of the refrigerant. Shown in figure 19.

Nomenclature / Greeks

D diameter (m)

x vapor quality(-)

G mass velocity (kg s ? m2 )

m& mass flow (kg s )

f friction factor

p pressure (Pa)

p* reduced pressure

Q heat flux (W/m2)

T temperature (K)

i entalpy (J/kg)

tp two phase

? heat transfer coefficient (W m2K)

λ thermal conductivity (W mK)

ρ density (kg m3 )

[varepsilon] void fraction

η dynamic viscosity (Ns m2 )

σ surface tension(-) Subscripts

f liquid

g vapor

References

REFERENCES

[1] Combating climate change The EU leads the way 2008 Edition, Catalogue number: NA-AB-08-128-EN-C.

[2] http://ujszechenyiterv.gov.hu/.

[3] MacArthur JW. Theoretical analysis of the dynamic interactions of vapour-compression heat pumps. Energy Conservation and Management 1984;24(1):49-66.

[4] Chi J, Didion D. A simulation model of the transient performance of a heat pump. Int J Refrigeration 1982; 5(3):176-84.

[5] P. Welsby, M. Pezzani, S. Devotta, P. J. Diggory, J. J. Guy: Steady- and dynamic-state simulations of heat-pumps. Part II: Modelling of a motor driven water-to-water heat-pump Original Research Article Applied Energy, Volume 31, Issue 4, 1988, Pages 239-262

[6] Wang H, Touber S. Distributed and non-steady state modeling of an air cooler. International Journal of Refrigeration 1991, 14, 98-111.

[7] Yasuda H, Touber S, Machielsen CHM. Simulation model of a vapor compression refrigeration system. ASHRAE Transactions 1983, 89(2A), 408-24

[8] J. Nyers and S. Gisbert, A dynamical model adequate for controlling the evaporator of heat pump, Int. J. Of Refrigeration 17 (2), 101-108 (1994).

[9] R.N.N. Koury, L. Machado, K.A.R. Ismail: Numerical simulation of a variable speed refrigeration system, International Journal of Refrigeration 24 (2001) 192-200, PII: S0140-7007(00)00014-1.

[10] Jong Won Choi, Gilbong Lee, Min Soo Kim:Numerical study on the steady state and transient performance of a multi-type heat pump system, International Journal of Refrigeration 34 (2011) 1157-1172, doi:10.1016/j.ijrefrig.2010.09.021.

[11] Belman, J. M., Navarro-Esbrí, J., Ginestar, D. and Milian V., Steady-state model of a variable speed vapor compression system using R134a as working fluid. International Journal of Energy Research, 34: 933-945. (2010) doi: 10.1002/er.1606.

[12] Yang Zhao, Zhao Haibo, Fang Zheng: Modeling and dynamic control simulation of unitary gas engine heat pump, Energy Conversion and Management 48 (2007) 3146-3153, doi:10.1016/j.enconman.2007.01.033.

[13] L. Garbai, Sz. Mehes: System Theory Modell of Heat Pumps, Gépészet 2008. pp.12., Budapest, 2008, Hungary.

[14] Satish G. Kandlikar,Heat transfer and fluid low in minichannels and microchannels, Mechanical Engineering Department, Rochester Institute of Technology, Elsevier Science (2005), ISBN: 0-0804-4527-6 , USA.

[15] Taborek J., Shell-and-Tube Heat Exchangers, Single-Phase Flow, Heat Exchanger Design Handbook (1983), Chapter 3.3, Hemisphere, New York, USA.

[16] B.S. Pethukov, V.N. Popov, Theoretical calculation of heat exchange in turbulent flow in tubes of an incompressible fluid with variable physical properties, High Temp., 1 (1963), pp. 69-83.

[17] M.M. Shah. A general correlation for heat transfer during film condensation in tubes. International Journal of Heat and Mass Transfer, 22(4):547-556, 1979.

[18] John R. Thome, Engineering Data Book III, Single-Phase Shell-Side Flows and Heat Transfer, Chapter 3, Swiss Federal Institute of Technology Lausanne CH-1015 Lausanne, 2004, Switzerland.

[19] Friedel L. 1979. Improved friction pressure drop correlations for horizontal and vertical two phase pipe flow. Paper E2, European Two Phase Flow Group Meeting, Ispra, Italy.

[20] S.Z. Rouhani and E.Axelsson, Calculation of volume void fraction in subcooled and quality region, Int. J. Heat Mass Transfer 13 (1970), 383 - 393.

[21] Verein Detscher Ingenieure VDI - Warmeatlas (VDI Heat Atlas), Chapter HBB, VDI - Gesellschaft Verfahrrenstechnik und Chemieingenieurwesen (GVC), Düsseldorf, 1993.

[22] N.Zuber and J.A. Findlay, Average volumetric concentration in two-phase flow systems, J.Heat Transfer 87 (1965), 453 - 468.

[23] Su Changsun. Advanced engineering thermodynamics. Higher Education Press; 1987.

[24] Jong Won Choi, Gilbong Lee, Min Soo Kim:Numerical study on the steady state and transient performance of a multi-type heat pump system, International Journal of Refrigeration 34 (2011) 1157-1172, doi:10.1016/j.ijrefrig.2010.09.021

[25] T.M. Harms, J.E. Braun, E.A. Groll, The impact of modeling complexity and two-phase flow parameters on the accuracy of system modeling for unitary air conditioners, HVAC&R Research Journal 10 (1) (2004) 5-20.