- Email: [email protected]
The novel use of phase change materials in refrigeration plant. Part 2: Dynamic simulation model for the combined system
Applied Thermal Engineering 27 (2007) 2902–2910 www.elsevier.com/locate/apthermeng
The novel use of phase change materials in refrigeration plant. Part 2: Dynamic simulation model for the combined system Fuqiao Wang a
a,*
, Graeme Maidment a, John Missenden a, Robert Tozer
b
London South Bank University, Faculty of Engineering, Science and the Built Environment, 103 Borough Road, London SE1 0AA, UK b Waterman Gore-Mechanical and Electrical Consulting Engineers, Versailles Court, 3 Paris Garden, London SE1 8ND, UK Received 23 August 2004; accepted 17 June 2005 Available online 11 August 2005
Abstract A dynamic mathematical model for coupling the refrigeration system and PCMs has been developed in this paper. Overall the model consists of the following basic components: a compressor, a condenser, an expansion valve, an evaporator cooler and a PCM heat exchanger. The model developed here, is based on a lumped-parameter method. The condenser and evaporator were treated as storage tanks at different states, which have a superheat region, a two-phase region and a sub-cooled region. In the single-phase region the parameters are considered homogeneous whereas in the two-phase region, the intensive properties are considered as in thermal equilibrium. The compressor model is considered as an adiabatic process; an isentropic efficiency is employed in this process. The expansion process in the thermostatic expansion valve is considered as an isenthalpic process. The PCM is treated as a onedimensional heat transfer model. The mathematical simulation in this study predicts the refrigerant states and dynamic coefficient of performance in the system with respect to time. The dynamic validation shows good agreement with the test result. Ó 2005 Published by Elsevier Ltd. Keywords: Dynamic simulation; Refrigeration cycle; Phase change materials (PCMs); Thermal energy storage
1. Introduction The mathematical simulation of a refrigeration plant was first demonstrated in the 1970s. The first mathematical model of a refrigeration system used algebraic equations derived from the assumption of steady state flow [1]. During transient operation, the refrigeration system components experience unsteady state operation. The refrigerant mass flow rate is continuously changing, which causes spatial variations in the refrigerant distribution in the system components as well as variable refrigerant states at the inlet and outlet of each component. Of the four major components in the vapour compression system, the transients in the heat exchang*
Corresponding author. Tel.: +44 20 78157634; fax: +44 20 78157699. E-mail address: [email protected] (F. Wang). 1359-4311/$ - see front matter Ó 2005 Published by Elsevier Ltd. doi:10.1016/j.applthermaleng.2005.06.009
ers are usually the slowest to respond and have the largest impact on system performance [2]. It is necessary to consider the mass distribution within the heat exchangers as a function of time and space, and this requires transient mass balances to allow for local storage. Thermal capacitances of heat exchangers and the refrigerant have to be considered to account for local energy storage. James and James [3] and Cleland [4] have reported comprehensive surveys of mathematical models in the area of industrial refrigeration system simulation. Bendapudi and Braun [2] also have reviewed the dynamic simulation of the refrigeration system. Generally, heat exchanger models dealing with compressible twophase flow fall into two categories: the lumped-parameter approach or the spatially distributed approach [7]. The lumped-parameter approach is generally employed to integrated with thermal environment simulation
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
2903
Nomenclature
N P R r T U V v X
area (m2) mass flow factor enthalpy (J kg1) mass flow rate (kg s1) rotation speed of the compressor (rpm) heat flux (W m2) compressor power (W) fusion heat of PCM (J kg1) stroke of piston (m) heat transfer fluid temperature (K) specific internal energy (J kg1) volume rate (m3 s1) work (kJ kg1) coordinate (m) specific heat (kJ kg1 K1) diameter (m) latent heat (J kg1) the refrigerant flashing back into the vapour region (kg s1) wetted perimeter (m) pressure (Pa) thermal resistance (°C m2 W1) radius (m) temperature (K) overall heat transfer coefficient (W m2 K1) volume (m3) specific volume (m3 kg1) displacement of valve needle (m)
Greeks a gv q gi k s
heat transfer coefficient (W m2 K1) volumetric efficiency of the compressor density (kg m3) compressor isentropic efficiency thermal conductivity (W m1 K1) time (s)
A Cd h m n q Pe rh S Tf u V_ W x Cp D hfg Dm
because of less calculation and proper accuracy. The spatial distributed approach can reveal details of the flow and heat transfer in the heat exchanger in the refrigeration cycle, however, more calculation and lengthy time consumption are required. In this paper, the lumped parameters approach is used because the PCM heat exchanger model is integrated. For lumped parameters model, Marshal and James [5] were the first to simulate a complex system. They used 46 ordinary differential equations (ODEs) and approximately 100 algebraic equations to simulate a vegetable freezer and its dedicated refrigeration system. Cleland [4] is probably the first to simulate a thermal environment which involves a large meat industry refrigeration system at three different evaporating tempera-
Subscripts a air side cond1 gas region in the condenser cond3 sub-cool region in the condenser evap evaporator evap2 saturation region in the evaporator G gas L liquid pcm phase change material r33 the outlet to the saturation region in the condenser r31 r4 the outlet to the condenser r55 the inlet to the gas superheat region in the evaporator sub sub-cool V,sat saturation area in the evaporator 00 vapour parameters, 0 initial state 2 outlet, zone 2 cond condenser cond2 saturation region in the condenser comp compressor evap1 superheat region in the evaporator f heat transfer fluid i inside o outside r3 the outlet to the compressor the outlet to the gas region in the condenser r5 the inlet to the evaporator r6 the outlet to the evaporator TEV thermal expansion valve w pipe wall 0 liquid parameters 1 inlet, zone 1 20 compressor outlet at real condition
tures. For the refrigeration system simulation, much work can also be found in Grald and MacArthur [6,7], He et al. [8,9], Willatzen et al. [10]. From these studies, the main assumptions of the transient model can be assumed as follows: (1) Liquid and vapour refrigerant in the heat exchangers are in thermal equilibrium. (2) Effects of pressure wave dynamics are negligible. (3) Expansion is isenthalpic. (4) Compression is adiabatic. (5) Thermal resistances of metallic elements in the system are negligible in comparison with other thermal resistances, however, their capacitance is important.
2904
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
The process of building a system transient model consists primarily of constructing transient models of individual components that are integrated to model the complete system. The details of each model are described below. 2. Mathematical model of refrigeration plant 2.1. Model of semi-hermetic reciprocating compressor The compressor generates a mass flow from the lowpressure side of the cycle to the high-pressure side of the cycle. The heat generated by the electric motor is transferred to the refrigerant. Thus the compressive process is an adiabatic process and the isentropic efficiency is employed, which is defined as h2 h1 T2 g¼ ¼f h2 h1 T1 6 5 T2 T2 ¼ 590.32 þ 3895.8 T1 T1 4 3 T2 T2 10; 686 þ 15; 638 T1 T1 2 T2 T2 12; 925 þ 5743.3 1075 ð1Þ T1 T1 The refrigerant mass flow of the compressor is given by Eq. (2) m ¼ q V_ gv ð2Þ p V_ ¼ D2 S n=60 ð3Þ 4 where the volumetric efficiency of the compressor, gv, is determined by p gv ¼ 0.851 0.0241 2 ð4Þ p1 The coefficiency listed above related to individual compressor, which is derived according to the tested compressor. The power of the compressor, Pe, can be calculated as P e ¼ mðh2 h1 Þ
ð5Þ
2.2. Condenser model The condenser determines the pressure level on the high-pressure side of the cycle. In the condenser [2], there are mainly three different heat-dissipating zones: a cool down zone for the superheated gas (Vcond1), a condensation zone (Vcond2 = VcondG + VcondL) and a sub-cooling zone (Vcond3) for liquid refrigerant as shown in Fig. 1. According to initial calculation, the superheat region takes a small part of the condenser. This region will be neglected in our model, and as such only two regions are considered; two-phase region and sub-cool region. The mass and energy equations can be expressed as follows. In saturated region, mass and energy balance equation are dðqcond2 V cond2 Þ dt mr3 hr3 mr33 hr33 Qcond1 Qcond2 þ Dmr dðqcond2 V cond2 ucond2 Þ ¼ dt mr3 mr33 þ Dm ¼
ð6Þ
ð7Þ
In liquid region, mass and energy balance equation are dV cond3 dt mr33 hr33 mr4 hr4 Qcond3 Dmr dðV cond3 ucond3 Þ ¼ qcond3 dt V cond2 ¼ V condL þ V condG
mr33 mr4 Dm ¼ qcond3
ð8Þ
ð9Þ ð10Þ
where Dm is the refrigerant flashing back into the vapour region, which occurs when the compressor rotation speed is suddenly reduced. This term was discovered important in order to simulate accurately the influence of the compressor speed variation. The result without this modification is negative sub-cooling which is not possible. Qcond2 and Qcond3 are calculated using NTU-e method is a method using for heat exchanger performance calculation, which is base on a known heat exchanger surface and structure. The equations are
Fig. 1. The concept model of condenser.
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
2905
omitted here, the details can be referred in Ref. [19]. The pipe energy equation can be written as C w q w Aw
oT w ¼ ai pDi ðT T w Þ þ a0 Ao gf ðT a T Þ ot
ð11Þ
The refrigerant properties state equations are f ðP ; v; T Þ ¼ 0
ð12Þ
Here M–H equations for refrigerants properties are used [19], in this model 10 refrigerants include R11, R12, R13, R14, R134a, R500, R114, R22, R717, R502 are later developed. Here, R22 is simulated. hr4 hr33 ¼ C p T sub
Fig. 3. A typical TEV.
ð13Þ
where, Cp is the specific heat of the refrigerant. 2.3. Evaporator model Similar to the condenser model, the model for the evaporator is shown in Fig. 2. However, in the evaporator, only two regions exist: the two-phase region and the superheat regions. According to this assumption, the mass and energy equations can be expressed as the following: mr5 hr5 mr55 hr55 þ qevap1 ¼
dðmuÞV;sat ds
ð14Þ
dmV;sat ds ¼ ðmuÞV;evapL þ ðmuÞV;evapG
mr5 mr55 ¼
ð15Þ
ðmuÞV;sat
ð16Þ
ðmÞV;sat ¼ ðmÞV;evapL þ ðmÞV;evapG
ð17Þ
Heat and mass storage for the superheated zone are negligible. Therefore mass and energy balance for the superheated region should be mr55 mr6 ¼ 0
ð18Þ
mr55 hr55 mr6 hr6 þ Qevap2 ¼ 0
ð19Þ
2.4. Thermostatic expansion valve (TEV) model The TEV shown in Figs. 3 and 4 is a proportional controller with the control signal from the superheat
Fig. 4. Diagram showing the geometrical relationship between the value needle and seat.
of the vapour at the evaporator outlet. It responds to the difference between the pressure of the refrigerant at the position that the pressure-sensing connection is made and the pressure developed in the temperaturesensing remote phial located at the outlet of the evaporator. The phial is normally charged with the same refrigerant as the plant. A spring setting in the valve is provided to adjust the required superheat of the vapour leaving the evaporator. The mass flow rate can be calculated as pffiffiffiffiffiffiffiffiffi m ¼ C d A0 qDP ð20Þ A0 ¼ f ðD; X ; aÞ ð21Þ where, A0 is the refrigerant flow area of the TEV, a is the angle of the valve seat as shown in Fig. 3.
Fig. 2. The concept of evaporator.
2906
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
X ¼ C 0 ðP 1 P 3 DP 0 Þ
ð22Þ
0
where C is the parameter of the TEV, which equals to the area of the diaphragm divided by the spring stiffness. P1 is the pressure in the remote phial and P3 is the evaporator pressure at the position of the pressure-sensing connection. DP0 is the spring pre-tension and is known as the static pressure of the valve. 3. Phase change material heat exchanger model Many theoretical studies concerning PCM storage have been undertaken. Heat transfer to and from a phase change material involves both convective and conductive processes. Amongst the studies reported in the literature, Cao and Faghri [11] are two of the few authors who have modeled convective heat transfer as well as conductive in a thermal storage system. The method used is similar to a model proposed by Lacroix [12] who presented a theoretical model for predicting the transient behaviour of a shell-and-tube storage unit with the PCM on the shell side and the heat transfer fluid circulating inside the tubes. This phase change problem was tackled with an enthalpy-based method combined with convective heat transfer from the fluid. The numerical results were validated with experimental data. Ismail [13] presented a two-dimensional axis-symmetrical model to present the convection dominated fusion of PCM around a vertical cylinder submersed in the PCM. The numerical predictions were compared with available numerical and experimental results, which showed good agreement. Kang and Zhang [14] built a general model for analyzing the thermal characteristics of a latent heat storage system. The investigations reported above involved the use of CFD (computational fluid dynamics) and considerable computer power to converge their solution. Also, these investigations were only concerned with heat transfer and did not connect to a system model. Modeling the refrigeration system coupled with PCM stores is complex because of the moving liquid/solid interface boundary problem, which occurs within the PCM. The liquid and solid interface will move during the charging and discharging period of the PCM store. The moving interface will cause the thermal resistance to change over time and this causes the heat transfer coefficient to vary and that requires the use of non-steady state governing equations. When combined with a refrigeration system, these parameters will affect the refrigeration system performance, which means all the parameters both for the PCM storage and for the refrigeration cycle need to be solved with respect to time. There are only few studies that have been cited in the literature on the modeling of the thermal performance of low temperature storage combined with a refrigeration system. Jekel [15] developed a model of ice-storage tanks. Finer
[16] published a simple mathematical model for predicting the transient behavior of an ice-bank system. This model was improved by Lovatt et al. [17] who integrated it with a refrigeration system. However, these models assume convection within the PCM is not important and also that the refrigerant in contact with the PCM is in single-phase only. The new model addresses this by allowing two-phase refrigerant flow in the heat transfer process, and also a novel approach which improves the accuracy of this method has been validated for use in this application. This uses an equivalent conduction that includes the effects of convection. 3.1. The common structure of PCM storage system The PCM heat exchanger was modeled as a shelland-tube storage unit with the PCM on the shell side and the heat transfer fluid (HTF) circulating inside the tubes. When the phase transition temperature is selected properly the PCM will melt or solidify depending on the inner refrigerant temperature. The energy will be stored in or be released from the PCM during the phase transition process. The reason for the choice of this configuration is that it maximizes heat transfer and minimizes the quantity of refrigerant charge and thus is environmentally acceptable. Furthermore, using the refrigeration in the outside of the tubes would involve complex vessel design. A solid form of the PCM was assumed around the outside tube when the PCM store was charging. The interface of the liquid and solid was assumed to move evenly away from the tube surface during the charging process. 3.2. Mathematical description of PCM model A number of assumptions have been made in describing the PCM model mathematically. These are (1) Axial conduction in the PCM is negligible, axial temperature gradient is small when compared to radial temperature gradient. (2) The Stefan number Ste 1, which means that the sensible thermal storage capacity is very small compared with latent thermal storage capacity and therefore can be neglected. (3) The capacitance of the heat transfer fluid (refrigerant) is small and can be neglected. (4) Natural convective heat transfer in the liquid PCM is assumed to be part of the equivalent conductive heat transfer, which increases the thermal conductivity by 50% according to the test data. (5) Heat storage in the metal tubes is neglected. With the above assumption, the energy balance equation for the PCM is written as follows:
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
oAðs; xÞ ¼ UN jT f ðs; xÞ T pcm j U o N o ðT pcm T a Þ os ð23Þ where A(s, x) is the PCM transition interface cross-sectional area, N is wetted perimeter of heat transfer fluid. Uo is the outside air convective heat transfer coefficient around the PCM container, No is the perimeter of PCM container. The energy equation of the heat transfer fluid is described by oT f ðs; xÞ ¼ UN jT f ðs; xÞ T pcm j mc ð24Þ ox where, m is the mass flow rate of the refrigerant passing through the PCM heat exchanger, C is the specific heat of the refrigerant. U and N as stated above, Tf is the refrigerant temperature, Tpcm is the temperature of PCM material. rhq
ð25Þ
U ¼ af Rr Rr ¼
Rf Rf þ Rw þ Rpcm
The initial conditions are Aðx; s ¼ 0Þ ¼ A0 ðxÞ T f ðx; s ¼ 0Þ ¼ T f0
ð26Þ ð27Þ
T pcm ðx; s ¼ 0Þ ¼ T 0
ð28Þ
and the boundary conditions are T ðs; x ¼ 0Þ ¼ T f oT ðr; sÞ k ¼ a ½T ðr; sÞ T f jr¼rr or r¼ri oT ðr; sÞ ¼ a ½T ðr; sÞ T air jr¼rpcm k or r¼rpcm
ð29Þ ð30Þ ð31Þ
2907
4. Dynamic results validation An experimental facility has been designed and constructed to enable validation of results and system optimization. This has been described in part 1 of the paper. 4.1. Basic system without PCM heat exchanger The simulation results are compared with data from test results and do not include the PCM heat exchanger. Fig. 5 shows the pressure comparisons between the tested and calculated data from start to nearly an hour. Fig. 6 is the COP (coefficient of performance) comparison. The main discrepancies occur at the beginning of the test. After this the model and experiment are in agreement within the test error. 4.2. PCM heat exchanger positioned between condenser and TEV PCM heat exchanger between condenser and TEV is charged with the materials E21, made by EPSLtd. The physical properties of the latent heat and transition temperature are listed in Table 1. The heat exchangers are about 5 m long in total and measure 110 mm in diameter. Fig. 7 shows the pressure comparisons between the tested and calculated data from start to nearly 70 min. Fig. 8 is the COP comparison. These figures show that the model results are in good agreement with the measured results. However, both cases with or without PCM heat exchanger show that the main errors between the model and measured results occur at the beginning of the operation. This might be caused by the following reasons:
2000000
Pressure (Pa)
1700000
1400000 Pe–cal
Pc–cal
Pc–test
Pe–test
1100000
800000
500000 0
500
1000
1500
2000
2500
3000
Time (sec) Fig. 5. Pressure comparisons between test data and calculation data.
3500
4000
2908
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
3
2.5
COP
2
cop_test
cop_cal 1.5
1
0.5 0
500
1000
1500
2000
2500
3000
3500
4000
Time (sec) Fig. 6. COP comparisons between test data and calculation data.
(2) The heat transfer coefficient for the condenser and evaporator used in the model are the correlations under steady state conditions. (3) The time step for the model is 0.2 s. However, the measurement time step is 5 s roughly. This results in the input data (inlet air temperature of both condenser and evaporator) being not exactly the same as the tested data.
Table 1 PCM materials physical properties
E21
Latent heat (kJ kg1)
Transition temperature (°C)
Thermal conductivity (W m K1)
Density (kg m3)
150
21
0.43
1480
(1) The initial mass distribution of refrigerant at different components is assumed that initial twophase equilibrium existed both in condenser and evaporator at initial state, as initial dry fraction in evaporator was assumed as well.
Fig. 9 shows a good agreement between the temperatures differences either side of the PCM for both test and model cases. The small deviation at the start of the cycle is possibly due to the sensitivity of the moving boundary
1800000 1600000 1400000
Pressure (Pa)
1200000
Pe_cal
Pc_test
Pc_cal
Pe_test
1000000 800000 600000 400000 200000 0 0
500
1000
1500
2000
2500
3000
3500
Time (sec) Fig. 7. Pressure comparison between test data and calculation data.
4000
4500
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
2909
3
2.5
COP
2 cop_cal
cop_test
1.5
1
0.5 0
500
1000
1500
2000
2500
3000
3500
4000
4500
4000
4500
Time (sec)
Fig. 8. COP comparison between test data and calculation data.
12
temp_difference_ cal
temp_difference_test
10
Temperature (ºC)
8
6
4
2
0 0
500
1000
1500
2000
2500
3000
3500
Time (sec) Fig. 9. Comparison of the refrigerant temperature difference across the PCM heat exchanger.
as the solid PCM melts or the initial mass distribution in both cases being different.
5. Conclusion A dynamic model of the novel system has been developed which can be used to design and optimize the performance of the system. In the simulation, a lumped parameters approach with moving boundary is used to ensure greater accuracy. The dynamic model developed in this paper incorporates the following features:
(1) Calculation of the refrigeration properties extends to ten refrigerants, including R11, R12, R13, R14, R134a, R500, R114, R22, R717, R502. Here, R22. (2) Application of the model also extends to variable frequency compressors. (3) A one-dimensional PCM model has been developed to simulate unsteady state conditions. This novel PCM heat exchanger model is integrated with the system. Also, a new simulation concept of liquid refrigerant flash into the gas region in the condenser is introduced
2910
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
into the model. This improves the accuracy of the model when the compressor rotation speed is reduced. A transient void fraction is used instead of mean void fraction in the evaporator. Furthermore, a novel PCM heat exchanger model with two fluids in phase change has been developed and integrated into refrigeration system. The lumped parameters model has been validated with the test results. It can be found that the model has a good agreement with the tested data for parameters such as condenser pressure, evaporator pressure and COP. For these parameters the model predicted the test data within 8%. However, this model presented some weakness in predicting the superheat and subcooling. This is due to the lumped parameters assumption. In the model, all the parameters calculated are the average parameters rather than the spatial distributed parameters. This is the main drawback of the model. As a dynamic model, the system parameters are sensitive to its initial conditions. Model validation results are also keen to this condition. The total mass of the refrigerant can be monitored. However, the initial mass distribution is hardly to be traced. This makes the dynamic validation more difficult, and results in the initial errors. Moreover, the time step for recording the dada is up to the thermocouples responding to its thermal environment. Also, there is always a time lag existed in testing process. This time lag will also results in some errors. Acknowledgement The financial support of the Engineering and Physical Sciences Research Council (EPSRC) is gratefully acknowledged. References [1] K. Cornwell, Thermal optimisation of refrigeration evaporators, in: Proceedings of the Institution of Refrigeration, 1992, pp. 1–11. [2] S. Bendapudi, J.E. Braun, A review of literature on dynamic models of vapour compression equipment, Sponsored by ASHRAE, Report # 4036-5, 2002.
[3] K.A. James, R.W. James, A Critical Survey of Dynamic Mathematical Models of Refrigeration Systems and Heat Pumps and Their Components. Technical Memorandum Õ97, Institute of Environmental Engineering, South Bank Polytechnic. 1986. [4] A.C. Cleland, Food Refrigeration Processes—Analysis, Design and Simulation, Elsevier Science, London, 1990. [5] S.A. Marshal, R.W. James, Dynamic analysis of an industrial refrigeration system to investigate capacity control, Proceedings of the Institution of Mechanical Engineers 189 (44) (1975) 437–445. [6] J.W. MacArthur, E.W. Grald, Unsteady compressible two-phase flow model for predicting cyclic heat pump performance and a comparison with experimental data, Int. J. Refrigr. 12 (1989) 29–41. [7] E.W. Grald, J.W. MacArthur, A moving Boundary formulation for modelling time-dependent two-phase flows, Int. J. Heat Fluid Flow 13 (3) (1992) 266–272. [8] X.D. He, S. Liu, H. Asada, Modeling of vapour compression cycles for advanced control in HVAC systems. in: Proceedings of the American Control Conference, Seattle, Washington, 1986. [9] X.D. He, S. Liu, H. Asada, A moving interface model of twophase flow heat exchanger dynamics for control of vapour compression cycle heat pump and refrigeration system design, analysis and applications. ASE Vol. 32 ASME, 1994. [10] M. Willatzen, N. Petit, L. Ploug-Sorensen, A general dynamic simulation model for evaporators and condensers in refrigeration: Part 1-moving boundary formulation of two-phase flows with heat exchange, Int. J. Refrigr. 21 (5) (1998) 398–403. [11] Y. Cao, A. Faghri, PCM/forced convection conjugate transient analysis of energy storage systems with annular and countercurrent flows, ASME, vol. 129, Heat Transfer Division, HTD, (Publication), 1990, pp. 63–69. [12] M. Lacroix, Study of heat transfer behavior of a latent heat thermal energy unit with a finned tube, Int. J. Heat Mass Transfer 36 (1993) 2083–2092. [13] K. Ismail, Convection-based model for a PCM vertical storage unit, Int. J. Energy Res. 22 (14) (1998) 1249–1265. [14] Y.B. Kang, Y. Zhang, General model for analysing the thermal characteristics of a class of latent heat thermal energy storage systems, J. Sol. Energy Eng., Trans. ASME 121 (4) (1999) 185– 193. [15] T.B. Jekel, Modelling of ice storage tanks, ASHRAE Trans. 99 (1) (1993) 1016–1024. [16] S.I. Finer, Simple mathematical model for predicting the transient behaviour of an ice storage system, Int. J. Refrigr. 16 (5) (1992) 312–320. [17] S.J. Lovatt, M.P. Loeffen, A.C. Cleland, Improved dynamic simulation of multi-temperature industrial refrigeration systems for food chilling, freezing and cold storage, Int. J. Refrigr. 21 (3) (1998) 247–260. [19] F. Wang, 2003, the passive use of phase change materials (PCMs) in refrigeration system, PhD thesis, London South Bank University, 2003.
The novel use of phase change materials in refrigeration plant. Part 2: Dynamic simulation model for the combined system Fuqiao Wang a
a,*
, Graeme Maidment a, John Missenden a, Robert Tozer
b
London South Bank University, Faculty of Engineering, Science and the Built Environment, 103 Borough Road, London SE1 0AA, UK b Waterman Gore-Mechanical and Electrical Consulting Engineers, Versailles Court, 3 Paris Garden, London SE1 8ND, UK Received 23 August 2004; accepted 17 June 2005 Available online 11 August 2005
Abstract A dynamic mathematical model for coupling the refrigeration system and PCMs has been developed in this paper. Overall the model consists of the following basic components: a compressor, a condenser, an expansion valve, an evaporator cooler and a PCM heat exchanger. The model developed here, is based on a lumped-parameter method. The condenser and evaporator were treated as storage tanks at different states, which have a superheat region, a two-phase region and a sub-cooled region. In the single-phase region the parameters are considered homogeneous whereas in the two-phase region, the intensive properties are considered as in thermal equilibrium. The compressor model is considered as an adiabatic process; an isentropic efficiency is employed in this process. The expansion process in the thermostatic expansion valve is considered as an isenthalpic process. The PCM is treated as a onedimensional heat transfer model. The mathematical simulation in this study predicts the refrigerant states and dynamic coefficient of performance in the system with respect to time. The dynamic validation shows good agreement with the test result. Ó 2005 Published by Elsevier Ltd. Keywords: Dynamic simulation; Refrigeration cycle; Phase change materials (PCMs); Thermal energy storage
1. Introduction The mathematical simulation of a refrigeration plant was first demonstrated in the 1970s. The first mathematical model of a refrigeration system used algebraic equations derived from the assumption of steady state flow [1]. During transient operation, the refrigeration system components experience unsteady state operation. The refrigerant mass flow rate is continuously changing, which causes spatial variations in the refrigerant distribution in the system components as well as variable refrigerant states at the inlet and outlet of each component. Of the four major components in the vapour compression system, the transients in the heat exchang*
Corresponding author. Tel.: +44 20 78157634; fax: +44 20 78157699. E-mail address: [email protected] (F. Wang). 1359-4311/$ - see front matter Ó 2005 Published by Elsevier Ltd. doi:10.1016/j.applthermaleng.2005.06.009
ers are usually the slowest to respond and have the largest impact on system performance [2]. It is necessary to consider the mass distribution within the heat exchangers as a function of time and space, and this requires transient mass balances to allow for local storage. Thermal capacitances of heat exchangers and the refrigerant have to be considered to account for local energy storage. James and James [3] and Cleland [4] have reported comprehensive surveys of mathematical models in the area of industrial refrigeration system simulation. Bendapudi and Braun [2] also have reviewed the dynamic simulation of the refrigeration system. Generally, heat exchanger models dealing with compressible twophase flow fall into two categories: the lumped-parameter approach or the spatially distributed approach [7]. The lumped-parameter approach is generally employed to integrated with thermal environment simulation
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
2903
Nomenclature
N P R r T U V v X
area (m2) mass flow factor enthalpy (J kg1) mass flow rate (kg s1) rotation speed of the compressor (rpm) heat flux (W m2) compressor power (W) fusion heat of PCM (J kg1) stroke of piston (m) heat transfer fluid temperature (K) specific internal energy (J kg1) volume rate (m3 s1) work (kJ kg1) coordinate (m) specific heat (kJ kg1 K1) diameter (m) latent heat (J kg1) the refrigerant flashing back into the vapour region (kg s1) wetted perimeter (m) pressure (Pa) thermal resistance (°C m2 W1) radius (m) temperature (K) overall heat transfer coefficient (W m2 K1) volume (m3) specific volume (m3 kg1) displacement of valve needle (m)
Greeks a gv q gi k s
heat transfer coefficient (W m2 K1) volumetric efficiency of the compressor density (kg m3) compressor isentropic efficiency thermal conductivity (W m1 K1) time (s)
A Cd h m n q Pe rh S Tf u V_ W x Cp D hfg Dm
because of less calculation and proper accuracy. The spatial distributed approach can reveal details of the flow and heat transfer in the heat exchanger in the refrigeration cycle, however, more calculation and lengthy time consumption are required. In this paper, the lumped parameters approach is used because the PCM heat exchanger model is integrated. For lumped parameters model, Marshal and James [5] were the first to simulate a complex system. They used 46 ordinary differential equations (ODEs) and approximately 100 algebraic equations to simulate a vegetable freezer and its dedicated refrigeration system. Cleland [4] is probably the first to simulate a thermal environment which involves a large meat industry refrigeration system at three different evaporating tempera-
Subscripts a air side cond1 gas region in the condenser cond3 sub-cool region in the condenser evap evaporator evap2 saturation region in the evaporator G gas L liquid pcm phase change material r33 the outlet to the saturation region in the condenser r31 r4 the outlet to the condenser r55 the inlet to the gas superheat region in the evaporator sub sub-cool V,sat saturation area in the evaporator 00 vapour parameters, 0 initial state 2 outlet, zone 2 cond condenser cond2 saturation region in the condenser comp compressor evap1 superheat region in the evaporator f heat transfer fluid i inside o outside r3 the outlet to the compressor the outlet to the gas region in the condenser r5 the inlet to the evaporator r6 the outlet to the evaporator TEV thermal expansion valve w pipe wall 0 liquid parameters 1 inlet, zone 1 20 compressor outlet at real condition
tures. For the refrigeration system simulation, much work can also be found in Grald and MacArthur [6,7], He et al. [8,9], Willatzen et al. [10]. From these studies, the main assumptions of the transient model can be assumed as follows: (1) Liquid and vapour refrigerant in the heat exchangers are in thermal equilibrium. (2) Effects of pressure wave dynamics are negligible. (3) Expansion is isenthalpic. (4) Compression is adiabatic. (5) Thermal resistances of metallic elements in the system are negligible in comparison with other thermal resistances, however, their capacitance is important.
2904
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
The process of building a system transient model consists primarily of constructing transient models of individual components that are integrated to model the complete system. The details of each model are described below. 2. Mathematical model of refrigeration plant 2.1. Model of semi-hermetic reciprocating compressor The compressor generates a mass flow from the lowpressure side of the cycle to the high-pressure side of the cycle. The heat generated by the electric motor is transferred to the refrigerant. Thus the compressive process is an adiabatic process and the isentropic efficiency is employed, which is defined as h2 h1 T2 g¼ ¼f h2 h1 T1 6 5 T2 T2 ¼ 590.32 þ 3895.8 T1 T1 4 3 T2 T2 10; 686 þ 15; 638 T1 T1 2 T2 T2 12; 925 þ 5743.3 1075 ð1Þ T1 T1 The refrigerant mass flow of the compressor is given by Eq. (2) m ¼ q V_ gv ð2Þ p V_ ¼ D2 S n=60 ð3Þ 4 where the volumetric efficiency of the compressor, gv, is determined by p gv ¼ 0.851 0.0241 2 ð4Þ p1 The coefficiency listed above related to individual compressor, which is derived according to the tested compressor. The power of the compressor, Pe, can be calculated as P e ¼ mðh2 h1 Þ
ð5Þ
2.2. Condenser model The condenser determines the pressure level on the high-pressure side of the cycle. In the condenser [2], there are mainly three different heat-dissipating zones: a cool down zone for the superheated gas (Vcond1), a condensation zone (Vcond2 = VcondG + VcondL) and a sub-cooling zone (Vcond3) for liquid refrigerant as shown in Fig. 1. According to initial calculation, the superheat region takes a small part of the condenser. This region will be neglected in our model, and as such only two regions are considered; two-phase region and sub-cool region. The mass and energy equations can be expressed as follows. In saturated region, mass and energy balance equation are dðqcond2 V cond2 Þ dt mr3 hr3 mr33 hr33 Qcond1 Qcond2 þ Dmr dðqcond2 V cond2 ucond2 Þ ¼ dt mr3 mr33 þ Dm ¼
ð6Þ
ð7Þ
In liquid region, mass and energy balance equation are dV cond3 dt mr33 hr33 mr4 hr4 Qcond3 Dmr dðV cond3 ucond3 Þ ¼ qcond3 dt V cond2 ¼ V condL þ V condG
mr33 mr4 Dm ¼ qcond3
ð8Þ
ð9Þ ð10Þ
where Dm is the refrigerant flashing back into the vapour region, which occurs when the compressor rotation speed is suddenly reduced. This term was discovered important in order to simulate accurately the influence of the compressor speed variation. The result without this modification is negative sub-cooling which is not possible. Qcond2 and Qcond3 are calculated using NTU-e method is a method using for heat exchanger performance calculation, which is base on a known heat exchanger surface and structure. The equations are
Fig. 1. The concept model of condenser.
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
2905
omitted here, the details can be referred in Ref. [19]. The pipe energy equation can be written as C w q w Aw
oT w ¼ ai pDi ðT T w Þ þ a0 Ao gf ðT a T Þ ot
ð11Þ
The refrigerant properties state equations are f ðP ; v; T Þ ¼ 0
ð12Þ
Here M–H equations for refrigerants properties are used [19], in this model 10 refrigerants include R11, R12, R13, R14, R134a, R500, R114, R22, R717, R502 are later developed. Here, R22 is simulated. hr4 hr33 ¼ C p T sub
Fig. 3. A typical TEV.
ð13Þ
where, Cp is the specific heat of the refrigerant. 2.3. Evaporator model Similar to the condenser model, the model for the evaporator is shown in Fig. 2. However, in the evaporator, only two regions exist: the two-phase region and the superheat regions. According to this assumption, the mass and energy equations can be expressed as the following: mr5 hr5 mr55 hr55 þ qevap1 ¼
dðmuÞV;sat ds
ð14Þ
dmV;sat ds ¼ ðmuÞV;evapL þ ðmuÞV;evapG
mr5 mr55 ¼
ð15Þ
ðmuÞV;sat
ð16Þ
ðmÞV;sat ¼ ðmÞV;evapL þ ðmÞV;evapG
ð17Þ
Heat and mass storage for the superheated zone are negligible. Therefore mass and energy balance for the superheated region should be mr55 mr6 ¼ 0
ð18Þ
mr55 hr55 mr6 hr6 þ Qevap2 ¼ 0
ð19Þ
2.4. Thermostatic expansion valve (TEV) model The TEV shown in Figs. 3 and 4 is a proportional controller with the control signal from the superheat
Fig. 4. Diagram showing the geometrical relationship between the value needle and seat.
of the vapour at the evaporator outlet. It responds to the difference between the pressure of the refrigerant at the position that the pressure-sensing connection is made and the pressure developed in the temperaturesensing remote phial located at the outlet of the evaporator. The phial is normally charged with the same refrigerant as the plant. A spring setting in the valve is provided to adjust the required superheat of the vapour leaving the evaporator. The mass flow rate can be calculated as pffiffiffiffiffiffiffiffiffi m ¼ C d A0 qDP ð20Þ A0 ¼ f ðD; X ; aÞ ð21Þ where, A0 is the refrigerant flow area of the TEV, a is the angle of the valve seat as shown in Fig. 3.
Fig. 2. The concept of evaporator.
2906
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
X ¼ C 0 ðP 1 P 3 DP 0 Þ
ð22Þ
0
where C is the parameter of the TEV, which equals to the area of the diaphragm divided by the spring stiffness. P1 is the pressure in the remote phial and P3 is the evaporator pressure at the position of the pressure-sensing connection. DP0 is the spring pre-tension and is known as the static pressure of the valve. 3. Phase change material heat exchanger model Many theoretical studies concerning PCM storage have been undertaken. Heat transfer to and from a phase change material involves both convective and conductive processes. Amongst the studies reported in the literature, Cao and Faghri [11] are two of the few authors who have modeled convective heat transfer as well as conductive in a thermal storage system. The method used is similar to a model proposed by Lacroix [12] who presented a theoretical model for predicting the transient behaviour of a shell-and-tube storage unit with the PCM on the shell side and the heat transfer fluid circulating inside the tubes. This phase change problem was tackled with an enthalpy-based method combined with convective heat transfer from the fluid. The numerical results were validated with experimental data. Ismail [13] presented a two-dimensional axis-symmetrical model to present the convection dominated fusion of PCM around a vertical cylinder submersed in the PCM. The numerical predictions were compared with available numerical and experimental results, which showed good agreement. Kang and Zhang [14] built a general model for analyzing the thermal characteristics of a latent heat storage system. The investigations reported above involved the use of CFD (computational fluid dynamics) and considerable computer power to converge their solution. Also, these investigations were only concerned with heat transfer and did not connect to a system model. Modeling the refrigeration system coupled with PCM stores is complex because of the moving liquid/solid interface boundary problem, which occurs within the PCM. The liquid and solid interface will move during the charging and discharging period of the PCM store. The moving interface will cause the thermal resistance to change over time and this causes the heat transfer coefficient to vary and that requires the use of non-steady state governing equations. When combined with a refrigeration system, these parameters will affect the refrigeration system performance, which means all the parameters both for the PCM storage and for the refrigeration cycle need to be solved with respect to time. There are only few studies that have been cited in the literature on the modeling of the thermal performance of low temperature storage combined with a refrigeration system. Jekel [15] developed a model of ice-storage tanks. Finer
[16] published a simple mathematical model for predicting the transient behavior of an ice-bank system. This model was improved by Lovatt et al. [17] who integrated it with a refrigeration system. However, these models assume convection within the PCM is not important and also that the refrigerant in contact with the PCM is in single-phase only. The new model addresses this by allowing two-phase refrigerant flow in the heat transfer process, and also a novel approach which improves the accuracy of this method has been validated for use in this application. This uses an equivalent conduction that includes the effects of convection. 3.1. The common structure of PCM storage system The PCM heat exchanger was modeled as a shelland-tube storage unit with the PCM on the shell side and the heat transfer fluid (HTF) circulating inside the tubes. When the phase transition temperature is selected properly the PCM will melt or solidify depending on the inner refrigerant temperature. The energy will be stored in or be released from the PCM during the phase transition process. The reason for the choice of this configuration is that it maximizes heat transfer and minimizes the quantity of refrigerant charge and thus is environmentally acceptable. Furthermore, using the refrigeration in the outside of the tubes would involve complex vessel design. A solid form of the PCM was assumed around the outside tube when the PCM store was charging. The interface of the liquid and solid was assumed to move evenly away from the tube surface during the charging process. 3.2. Mathematical description of PCM model A number of assumptions have been made in describing the PCM model mathematically. These are (1) Axial conduction in the PCM is negligible, axial temperature gradient is small when compared to radial temperature gradient. (2) The Stefan number Ste 1, which means that the sensible thermal storage capacity is very small compared with latent thermal storage capacity and therefore can be neglected. (3) The capacitance of the heat transfer fluid (refrigerant) is small and can be neglected. (4) Natural convective heat transfer in the liquid PCM is assumed to be part of the equivalent conductive heat transfer, which increases the thermal conductivity by 50% according to the test data. (5) Heat storage in the metal tubes is neglected. With the above assumption, the energy balance equation for the PCM is written as follows:
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
oAðs; xÞ ¼ UN jT f ðs; xÞ T pcm j U o N o ðT pcm T a Þ os ð23Þ where A(s, x) is the PCM transition interface cross-sectional area, N is wetted perimeter of heat transfer fluid. Uo is the outside air convective heat transfer coefficient around the PCM container, No is the perimeter of PCM container. The energy equation of the heat transfer fluid is described by oT f ðs; xÞ ¼ UN jT f ðs; xÞ T pcm j mc ð24Þ ox where, m is the mass flow rate of the refrigerant passing through the PCM heat exchanger, C is the specific heat of the refrigerant. U and N as stated above, Tf is the refrigerant temperature, Tpcm is the temperature of PCM material. rhq
ð25Þ
U ¼ af Rr Rr ¼
Rf Rf þ Rw þ Rpcm
The initial conditions are Aðx; s ¼ 0Þ ¼ A0 ðxÞ T f ðx; s ¼ 0Þ ¼ T f0
ð26Þ ð27Þ
T pcm ðx; s ¼ 0Þ ¼ T 0
ð28Þ
and the boundary conditions are T ðs; x ¼ 0Þ ¼ T f oT ðr; sÞ k ¼ a ½T ðr; sÞ T f jr¼rr or r¼ri oT ðr; sÞ ¼ a ½T ðr; sÞ T air jr¼rpcm k or r¼rpcm
ð29Þ ð30Þ ð31Þ
2907
4. Dynamic results validation An experimental facility has been designed and constructed to enable validation of results and system optimization. This has been described in part 1 of the paper. 4.1. Basic system without PCM heat exchanger The simulation results are compared with data from test results and do not include the PCM heat exchanger. Fig. 5 shows the pressure comparisons between the tested and calculated data from start to nearly an hour. Fig. 6 is the COP (coefficient of performance) comparison. The main discrepancies occur at the beginning of the test. After this the model and experiment are in agreement within the test error. 4.2. PCM heat exchanger positioned between condenser and TEV PCM heat exchanger between condenser and TEV is charged with the materials E21, made by EPSLtd. The physical properties of the latent heat and transition temperature are listed in Table 1. The heat exchangers are about 5 m long in total and measure 110 mm in diameter. Fig. 7 shows the pressure comparisons between the tested and calculated data from start to nearly 70 min. Fig. 8 is the COP comparison. These figures show that the model results are in good agreement with the measured results. However, both cases with or without PCM heat exchanger show that the main errors between the model and measured results occur at the beginning of the operation. This might be caused by the following reasons:
2000000
Pressure (Pa)
1700000
1400000 Pe–cal
Pc–cal
Pc–test
Pe–test
1100000
800000
500000 0
500
1000
1500
2000
2500
3000
Time (sec) Fig. 5. Pressure comparisons between test data and calculation data.
3500
4000
2908
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
3
2.5
COP
2
cop_test
cop_cal 1.5
1
0.5 0
500
1000
1500
2000
2500
3000
3500
4000
Time (sec) Fig. 6. COP comparisons between test data and calculation data.
(2) The heat transfer coefficient for the condenser and evaporator used in the model are the correlations under steady state conditions. (3) The time step for the model is 0.2 s. However, the measurement time step is 5 s roughly. This results in the input data (inlet air temperature of both condenser and evaporator) being not exactly the same as the tested data.
Table 1 PCM materials physical properties
E21
Latent heat (kJ kg1)
Transition temperature (°C)
Thermal conductivity (W m K1)
Density (kg m3)
150
21
0.43
1480
(1) The initial mass distribution of refrigerant at different components is assumed that initial twophase equilibrium existed both in condenser and evaporator at initial state, as initial dry fraction in evaporator was assumed as well.
Fig. 9 shows a good agreement between the temperatures differences either side of the PCM for both test and model cases. The small deviation at the start of the cycle is possibly due to the sensitivity of the moving boundary
1800000 1600000 1400000
Pressure (Pa)
1200000
Pe_cal
Pc_test
Pc_cal
Pe_test
1000000 800000 600000 400000 200000 0 0
500
1000
1500
2000
2500
3000
3500
Time (sec) Fig. 7. Pressure comparison between test data and calculation data.
4000
4500
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
2909
3
2.5
COP
2 cop_cal
cop_test
1.5
1
0.5 0
500
1000
1500
2000
2500
3000
3500
4000
4500
4000
4500
Time (sec)
Fig. 8. COP comparison between test data and calculation data.
12
temp_difference_ cal
temp_difference_test
10
Temperature (ºC)
8
6
4
2
0 0
500
1000
1500
2000
2500
3000
3500
Time (sec) Fig. 9. Comparison of the refrigerant temperature difference across the PCM heat exchanger.
as the solid PCM melts or the initial mass distribution in both cases being different.
5. Conclusion A dynamic model of the novel system has been developed which can be used to design and optimize the performance of the system. In the simulation, a lumped parameters approach with moving boundary is used to ensure greater accuracy. The dynamic model developed in this paper incorporates the following features:
(1) Calculation of the refrigeration properties extends to ten refrigerants, including R11, R12, R13, R14, R134a, R500, R114, R22, R717, R502. Here, R22. (2) Application of the model also extends to variable frequency compressors. (3) A one-dimensional PCM model has been developed to simulate unsteady state conditions. This novel PCM heat exchanger model is integrated with the system. Also, a new simulation concept of liquid refrigerant flash into the gas region in the condenser is introduced
2910
F. Wang et al. / Applied Thermal Engineering 27 (2007) 2902–2910
into the model. This improves the accuracy of the model when the compressor rotation speed is reduced. A transient void fraction is used instead of mean void fraction in the evaporator. Furthermore, a novel PCM heat exchanger model with two fluids in phase change has been developed and integrated into refrigeration system. The lumped parameters model has been validated with the test results. It can be found that the model has a good agreement with the tested data for parameters such as condenser pressure, evaporator pressure and COP. For these parameters the model predicted the test data within 8%. However, this model presented some weakness in predicting the superheat and subcooling. This is due to the lumped parameters assumption. In the model, all the parameters calculated are the average parameters rather than the spatial distributed parameters. This is the main drawback of the model. As a dynamic model, the system parameters are sensitive to its initial conditions. Model validation results are also keen to this condition. The total mass of the refrigerant can be monitored. However, the initial mass distribution is hardly to be traced. This makes the dynamic validation more difficult, and results in the initial errors. Moreover, the time step for recording the dada is up to the thermocouples responding to its thermal environment. Also, there is always a time lag existed in testing process. This time lag will also results in some errors. Acknowledgement The financial support of the Engineering and Physical Sciences Research Council (EPSRC) is gratefully acknowledged. References [1] K. Cornwell, Thermal optimisation of refrigeration evaporators, in: Proceedings of the Institution of Refrigeration, 1992, pp. 1–11. [2] S. Bendapudi, J.E. Braun, A review of literature on dynamic models of vapour compression equipment, Sponsored by ASHRAE, Report # 4036-5, 2002.
[3] K.A. James, R.W. James, A Critical Survey of Dynamic Mathematical Models of Refrigeration Systems and Heat Pumps and Their Components. Technical Memorandum Õ97, Institute of Environmental Engineering, South Bank Polytechnic. 1986. [4] A.C. Cleland, Food Refrigeration Processes—Analysis, Design and Simulation, Elsevier Science, London, 1990. [5] S.A. Marshal, R.W. James, Dynamic analysis of an industrial refrigeration system to investigate capacity control, Proceedings of the Institution of Mechanical Engineers 189 (44) (1975) 437–445. [6] J.W. MacArthur, E.W. Grald, Unsteady compressible two-phase flow model for predicting cyclic heat pump performance and a comparison with experimental data, Int. J. Refrigr. 12 (1989) 29–41. [7] E.W. Grald, J.W. MacArthur, A moving Boundary formulation for modelling time-dependent two-phase flows, Int. J. Heat Fluid Flow 13 (3) (1992) 266–272. [8] X.D. He, S. Liu, H. Asada, Modeling of vapour compression cycles for advanced control in HVAC systems. in: Proceedings of the American Control Conference, Seattle, Washington, 1986. [9] X.D. He, S. Liu, H. Asada, A moving interface model of twophase flow heat exchanger dynamics for control of vapour compression cycle heat pump and refrigeration system design, analysis and applications. ASE Vol. 32 ASME, 1994. [10] M. Willatzen, N. Petit, L. Ploug-Sorensen, A general dynamic simulation model for evaporators and condensers in refrigeration: Part 1-moving boundary formulation of two-phase flows with heat exchange, Int. J. Refrigr. 21 (5) (1998) 398–403. [11] Y. Cao, A. Faghri, PCM/forced convection conjugate transient analysis of energy storage systems with annular and countercurrent flows, ASME, vol. 129, Heat Transfer Division, HTD, (Publication), 1990, pp. 63–69. [12] M. Lacroix, Study of heat transfer behavior of a latent heat thermal energy unit with a finned tube, Int. J. Heat Mass Transfer 36 (1993) 2083–2092. [13] K. Ismail, Convection-based model for a PCM vertical storage unit, Int. J. Energy Res. 22 (14) (1998) 1249–1265. [14] Y.B. Kang, Y. Zhang, General model for analysing the thermal characteristics of a class of latent heat thermal energy storage systems, J. Sol. Energy Eng., Trans. ASME 121 (4) (1999) 185– 193. [15] T.B. Jekel, Modelling of ice storage tanks, ASHRAE Trans. 99 (1) (1993) 1016–1024. [16] S.I. Finer, Simple mathematical model for predicting the transient behaviour of an ice storage system, Int. J. Refrigr. 16 (5) (1992) 312–320. [17] S.J. Lovatt, M.P. Loeffen, A.C. Cleland, Improved dynamic simulation of multi-temperature industrial refrigeration systems for food chilling, freezing and cold storage, Int. J. Refrigr. 21 (3) (1998) 247–260. [19] F. Wang, 2003, the passive use of phase change materials (PCMs) in refrigeration system, PhD thesis, London South Bank University, 2003.