release in a plate-fin unit

Applied Thermal Engineering 31 (2011) 3871e3884

Contents lists available at ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

Numerical simulation on phase-change thermal storage/release in a plate-fin unit Wei-Biao Ye a, Dong-Sheng Zhu a, b, *, Nan Wang a a

Key Lab. of Enhanced Heat Transfer and Energy Conservation, Ministry of Education, School of Chemistry and Chemical Engineering, South China University of Technology, Guangzhou 510640, PR China b School of Mechanical and Power Engineering, East China University of Science and Technology, Shanghai 200237, PR China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 11 April 2011 Accepted 19 July 2011 Available online 26 July 2011

The fluid flow and heat transfer in a plate-fin unit with a characteristic length of 2 mm used for rapid heat storage/release by paraffin (phase change material, PCM) are investigated numerically. Transient simulations are performed using a commercial computational fluid dynamics (CFD) package, Fluent, based on the finite volume method. The effect of temperature differences on the fluid flow and heat transfer in the energy storage/release system is analyzed. It is found that temperature differences play a key role in the performances of energy storage when temperature differences are less than 20
C. It is noted that part of not solidified PCM can be observed clearly during energy release, and a vortex in the air region is formed remarkably at the moment of complete thermal energy release. The correlations are developed as a function of the associated variables. The obtained correlations are useful for future component design and system optimization. Additionally, Experimental data taken from the literatures are conducted to validate the model. The numerical results show a good agreement with the experimental ones. Crown Copyright Ó 2011 Published by Elsevier Ltd. All rights reserved.

Keywords: Numerical prediction Latent thermal storage/release Plate-fin unit Temperature differences Correlations

1. Introduction Thermal energy storage/release systems can be used as either a heat source or heat sink for energy conversion and heat rejection systems. Further, it plays an important role in the rationale engineering use of energy and makes the system more cost-effective by reducing the wastage of energy. How to use thermal energy storage/release device effectively have expanded notably in recent decades. It allows the decoupling between supply and consumption of energy. In the applications with intermittent energy generation, such as waster heat recovery, ground heat thermal and solar thermal systems, an appropriate thermal storage/release system is essential. The thermal storage/release technology based on the use of phase change materials (PCMs), which possess a great capacity of accumulation energy for consideration as heat storage media, has raised an important practical interest [1]. This is mainly due to the high energy storage density during phase change process within a very narrow temperature range. These materials are used in

* Corresponding author. Key Lab. of Enhanced Heat Transfer and Energy Conservation, Ministry of Education, School of Chemistry and Chemical Engineering, South China University of Technology, Guangzhou 510640, PR China. Tel./ fax: þ86 20 87114568. E-mail address: [email protected] (D.-S. Zhu).

applications where it is necessary to store/release energy thanks to the temporary phase change between the offer and demand of thermal energy. The use of PCMs in energy storage/release has been widely studied numerically and experimentally. Ho and Viskanta [2] reported an experimental study on the evolution of melt front profiles and heat transfer. The effect of initial subcooling on the solid n-Octadecane in a two-dimensional rectangular cavity consisting of a heated bottom and two conducting vertical sidewalls was performed. The research showed that buoyancy force plays an important role in energy storage process. However, heat conduction was the dominant mode during energy release. Additionally, a correlation was developed between Nusselt number and Rayleigh number. Short extended surfaces were more effective than longer ones for solideliquid phase change. Gadgil and Gobin [3] simulated numerically two-dimensional thermal energy storage process of a solideliquid PCM in a rectangular enclosure heated from one side wall. They predicted the shape and motion of the solideliquid boundary by divided the process in a large number of quasi-static steps. Natural convection of steady state, in each quasi-static step, was calculated by directly solving the governing equations of motion with a finite difference method (FDM). The calculated results showed good agreement with the experiments. At the same time, Okada [4] analyzed natural convection of energy storage in transient melt region by means of FDM using n-Octadecane as a PCM. Analytical results found that the

1359-4311/$ e see front matter Crown Copyright Ó 2011 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2011.07.035

3872

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

Nomenclature A a C cp Fo g h H k L Nu P PCM q’ R2 ! S Ste t T u, v U ! V

porosity function thermal diffusivity coefficient, m2/s constant specific heat, J/kg K Fourier number gravity acceleration, m/s2 specific enthalpy, J/kg height, m thermal conductivity, W/m K latent heat capacity, J/kg Nusselt number pressure, Pa phase change material heat flux, W/m2 coefficient of determination momentum source term, N/m3 Stefan number jtime, s temperature, K velocity magnitude in x and y coordinates, respectively, m/s heat transfer coefficient, W/m2 K velocity vector, m/s

profile of the melting front and the temperature distribution in the melt region agreed well with the experimental results. Additionally, a simple relationship between the average Nusselt number along the vertical wall and dimensionless time was obtained. Inaba et al. [5] investigated experimentally in dealing with natural convective heat transfer characteristics of PCM composed of water and surfactant, in rectangular enclosures heated from below and cooled from above. The results showed that there was not much differences in natural convective heat transfer characteristics of the PCM slurry with low PCM concentration (<10% in mass). However, the difference became great with increasing the PCM concentration due to contribution of latent heat generated by melting and solidification. Pal and Joshi [6] performed, experimentally and computationally, melting of n-Triacontane in a side heated tall enclosure, by a uniformly dissipating heat source. An implicit enthalpy-porosity approach was used for computational modeling of the thermal energy storage process. The predictions were found to be in good agreement with experiment. The increasing importance of energy managements and use of energy attracted interest in PCMs, mainly concerning optimization of heat conduction and convection in reservoirs with planar, cylindrical or spherical configuration. However, almost all organic PCMs have low thermal conductivity which frequently makes an anticipated level of thermal storage untenable within an acceptable time period. The rate of heat transfer can be enhanced by incorporating high thermal conductivity materials, known as thermal conductivity enhancer (TCE) into the PCM. TCE is distributed in the PCM in the form of metal matrix [7,8], fins [9] and uniformly dispersed high thermal conductivity particles (e.g. graphitecompound matrices [10,11] and carbon fiber [12]). TCE also discussed in detail by Agyenim et al. [13], Velraj [14], Jegadheeswaran and Pohekar [15], and Zalba et al. [16]. Cabeza et al. [17] experimentally performed three different heat transfer enhancement methods, i.e., stainless steel pieces, copper pieces, and graphite matrix impregnated with water (PCM). The results found that the PCM-graphite composite material shows a very effective heat flux

W x, y, z

width, m Cartesian coordinates

Greek symbols volume fraction in multiphase media thermal expansion coefficient, 1/K Δ difference Q dimensionless mean temperature m dynamic viscosity, kg/m s s dimensionless time ɸ liquid fraction r density, kg/m3 ε constant

a b

Subscripts 0 initial value avg average c complete value f fin L liquid PCM-phase n nth phase pcm phase change material ref reference value w wall of plate

enhancement with an increase in heat flux bigger than that with any of the other techniques due to its generations of a consolidated block with many uninterrupted heat flux paths. Recently, Medrano et al. [18] experimentally investigated five different tube heat exchangers as heat storage systems. Results also showed that heat exchanger with PCM embedded in a graphite matrix has the higher heat transfer coefficient (700e800 W/m2 K). Furthermore, Akgun et al. [19,20] experimentally investigated the melting and solidification processes of paraffin as a PCM in a novel tube in shell heat exchanger system. Nevertheless, modeling of melting/solidification in plate-fin geometry is a considerable challenge. In addition to all phase-change problems, we encounter here such phenomena as convection in the melt, volumetric expansion, motion of the solid in the melt, and opposite effect will be found in solidification. In case of fins, TCE can be distributed in two ways, (I) plate type fins and (II) pin type fins. Stritih [21] tested the heat transfer characteristics of a latent heat storage unit with a finned surface in terms of the energy storage/release processes by comparing them with a plain surface heat storage unit. The heat storage was not a problem during thermal storage applications. The thermal release can be effectively enhanced with fins. The effect of natural convection, which was dominant during thermal energy storage and increased the heat transfer, was reduced due to the fins. If the fins existed, the heat transfer was greater during energy release and 40% reduction in release time. However, fluid flow and heat transfer in latent heat storage/ release system is a very complex problem. The process is nonstationary, there is also a phase-change problem, and the process is non-homogeneous due to the fins. Time-dependent melting and solidification of the process is the most important for thermal energy applications. Therefore, numerical methods provide us more convenient and detailed information to investigate phasechange problems due to the complexity of nonlinear governing equations. So far, there have been a few reports about the phase change thermal energy storage/release process (paraffin as PCM) in a plate-

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

fin unit by visual study numerically. Furthermore, a long time for energy storage/release was extensively evaluated by Refs. [22e26]. However, in some practical applications, an objective requires to store/release as much heat within a few minutes, i.e., fast charge/ discharge. Therefore, the latent thermal storage/release system with small characteristic size can be used for rapidly supplying a large amount of heat to reach the objective. In the present study, a thermal energy storage/release system employing PCM having a characteristic length of 2 mm for rapid heat storage/release is investigated numerically using computational fluid dynamics (CFD) approach. Recently, modeling of the fluid flow and phase-change heat transfer in thermal storage/release processes has attracted considerable attention of the researchers. The paper reports quantitatively the detailed results of numerical modeling of the phase-change thermal storage/release in a plate-fin unit, and paraffin is used as PCM. Especially, the performances of thermal release process are proposed for the first time in present research. The purposes of this research are two aspects. One is to supply detail data/figure of thermal energy storage/release process for

3873

a real practical engineering application, and the other one is to figure out the correlations between influential variables. This study would provide useful information not only for recognizing the effect of temperature differences on the thermal energy storage/ release process but also for the design of plate-fin latent thermal system. 2. Numerical analysis A physical model of the system is presented at first. Then, the conservation equations, boundary and initial conditions, CFD solution procedure, time-step and grid independence test, and experimental validation are discussed in detail. 2.1. Physical model The three- and two-dimensional map of the thermal energy storage/release system is presented in Fig. 1. The thermal energy storage/release system consists of square channel with inner plates, outer plates and fins space filled with PCM inside the cavities (see

Fig. 1. Geometry of the thermal energy storage/release system: (a) Actual geometry; (b) geometry used for three-dimensional unit; (c) two-dimensional computational domain of front view; (d) boundary conditions for two-dimensional simulation.

3874

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

Table 1 Parameters for ΔT used in the study. case

ΔT for energy storage,
C

ΔT for energy release,
C

1 2 3 4 5 6 7

5 10 15 20 25 30 35

20

2.2. Governing equations

Fig. 1a). The heat transfer fluid (HTF), water, flows through the inner plates and exchanges heat with PCM through heating/cooling plate wall. The schematic of the computational domain for three- and two-dimensional are illustrated in Fig. 1b and c, respectively, due to the symmetric arrangement. To take an account of volume expansion during solideliquid phase change due to a large difference in solid and liquid density of PCM that exists in reality. The height of the fins is larger than PCM, which fills 85% of the enclosed space. In other words, the other 15% is filled up with air above the PCM. Geometrical dimensions are explored in this work. The PCM layer of width Wpcm ¼ 4 mm, and the vertical fins of height Hf ¼ 12 mm and width Wf ¼ 1.6 mm are included, as shown in Fig. 1c. These values are scaled according to the practical application in order to obtain the fluid flow and heat transfer performances more accurately. The fins and plates are all made of aluminum to ensure high thermal conductivity and low cost. The density of air is obtained as a polynomial function of temperature while constant thermophysical properties are specified for aluminum. Following the same geometry parameters, numerical calculations are performed for various temperature differences (ΔT) between the temperature of heating/cooling plate wall (Tw) and the PCM mean melting temperature (Tpcm,avg) during the process of energy storage/release, i.e., ΔT ¼ jTweTpcm, avgj. However, the temperature difference is kept all the same for all cases on energy release, namely 20
C as reflected in Table 1. The properties of the PCM, paraffin C19-C20 [27], with a melting interval of 33e35
C as discussed below in detail. Material density of the PCM-phase depends on their temperature. For the liquid PCM-phase, considering computational continuity during phase change, the density is calculated with the following relations:

rPCM ¼

rL






(1)

b  T  Tpcm þ 1

where rL is the density of liquid PCM-phase at the phase change temperature Tpcm, and b is the thermal expansion coefficient. Here b ¼ 0.001 K1 based on the analysis of the detailed date followed by Humphries and Griggs [28]. The dynamic viscosity of liquid PCM-phase is expressed as:



m ¼ exp A þ

B T

Tpcm, latent heat L and specific heat cp, are given respectively. The detailed properties of materials, used for computation for aluminum and air, are listed in Table 2.



The PCM and air between the neighboring fins are assumed laminar owing to the low velocity and small fin pitch. In order to describe the PCM-air system with a moving internal interface but without inter-penetration of the two phase fluids, a so-called “volume-of-fluid” (VOF) model has been used successfully [30,31]. In the VOF model, a single set of conservation equations is shared by the phases, the each phases volume fraction is tracked in each computational cell throughout the domain by interface reconstruction [32]. Accordingly, conservation equations, employed for modeling the PCM-air system in a Cartesian coordinate system by enthalpyporosity approach [6,33], are as follow:

Continuity

Momentum

Energy

(3)

D! v ! ! v þrg þ S ¼ VP þ mV2 ! Dt

r

Dh ¼ kV2 T Dt

r

(4)

(5)

Where an is denotes the nth fluid’s volume fraction (i.e., the following three conditions are possible: if an ¼ 0 the cell does not contain the nth fluid; if an ¼ 1 the cell is full of the nth fluid; if 0 < an < 1 the cell contains the interface between the nth fluid and one or more other fluids.), k is the thermal conductivity, t is the time, ! ! v is the velocity vector. S is the momentum source term for PCM ! v , where A (f) is the “porosity funcphase, given by S ¼ AðfÞ! tion” defined by Brent et al. [34] to make the momentum equation mimic CarmaneKozeny equations for flow in porous media:

AðfÞ ¼

Cð1  fÞ2

(6)

f3 þ ε

where ε ¼ 0.001 is a small computational constant used to avoid division by zero, and C ¼ 1.0  105 is a mushy zone constant reflecting the morphology of the melting front. RT Also, h is the specific enthalpy defined as h ¼ href þ Tref Cp dT, and the enthalpy change due to the phase change fL, where href is the reference enthalpy at the reference temperature Tref, cp is the specific heat, L is the latent heat capacity of the PCM, and f is the liquid fraction during the phase-change which occurs over a range of temperatures Tsolid ＜ T ＜ Tliquid, defined by the follow relations:

f ¼ 0;

if T < Tsolid

(2)

where A ¼ 4.25 and B ¼ 1790 are constant coefficients as suggested by Reid et al. [29]. In addition, the thermophysical properties of the PCM-phase, including thermal conductivity k, melting point

Dan ¼ 0 Dt

if Tsolid < T < Tliquid if T>Tliquid

f ¼

T  Tsolid Tliquid  Tsolid

(7)

f¼ 1

Table 2 Properties of paraffin, air and aluminum used for simulation. Materials

r (kg/m3)

k (W/m K)

cp (J/kg K)

Tpcm (
C)

L (J/kg)

m (kg/m s)

Paraffin Liquid Solid Air Aluminum

750 0:001  ðT  308:15Þ þ 1 850 1.2  105T20.01134T þ 3.498 2719

0.15 0.22 0.0242 202.4

2630

33e35

176000

1006.43 871

e e

e e

0:001  exp  4:25 þ e 1.7894  105 e

1790 T

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

The numerical simulation model mentioned above can be found in Assis et al. [35] and Shatikian et al. [36]. They simulated the performances of PCM in spherical geometry and in internal fins heat storage unit, respectively. It can be observed that the obtained results were in adequate agreement with experimental approach presented by Assis et al. [35], Ghosh et al. [36], and Shmueli et al. [37]. Thus, this method is now employed to simulate the latent thermal energy storage/release system in plate-fin heat exchanger unit, which based on the assumption that liquid PCM and air are compressible inside the plate-fin thermal energy storage/release systems. 2.3. Boundary and initial conditions In the present study, characterization of phase-change thermal energy storage and release processes in a plate-fin unit for uniform heating/cooling from the bottom plate wall are presented. The diagram of boundary condition for two-dimensional plate-fin thermal energy storage/release unit is shown in Fig. 1d. The boundary conditions are described as follow: The coupled-wall boundary condition is specified at the interface between PCM and fin to allow for conjugate heat transfer (Tf,wall ¼ Tpcm/air,A constant wall,kf ðvTf =vxÞjwall ¼ kpcm =airððvTpcm =airÞ=vxÞjwall ). temperature, which represents the temperature of heat transfer fluid, is prescribed at the heating/cooling wall (Twall). Therefore, the heat is transferred through fins by conduction first and then through the PCM by convection. Insulated boundary condition is specified for outer plates wall ðkwall vT=vy ¼ 0Þ. Symmetric boundary condition are used in fins and PCM, for fins: u ¼ v ¼ 0, vT/ vx ¼ 0; for PCM: vT/vx ¼ 0, vv/vx ¼ 0, u ¼ 0. Because of the proposed latent thermal storage/release system in a plate-fin unit is transient, initial conditions must be set appropriately before calculation. In the simulations, the initial condition of the whole system is supposed to be an average surrounding ambient temperature of 30
C, i.e., the PCM is slightly sub-cooling during the process of thermal storage. However, the instant of complete thermal storage will be considered as the initial conditions of thermal release process.

3875

simulations, the time-step and grid size independence study have been carefully made in preliminary calculations. The numerical results can be almost regarded as time-step and grid size independence as far as the heat flux of heating/cooling plate wall; liquid fraction and solid fraction are converged. In particular, the time-step is as small as 0.004, 0.006 (chosen), and 0.008 s for energy storage, and 0.002, 0.003 (chosen), and 0.004 s for energy release, respectively. As can be seen from Fig. 2, the results obtained for heat flux, liquid fraction and solid fraction became practically independent of the time-step, through the whole process, i.e., from t ¼ 0 to the complete energy storage/ release of the plate-fin unit. Here, the “þ”sign in the heat flux represents energy storage and the “” sign represents energy release (see Fig. 2a). Four quadrilateral grid systems are studied, namely 13  67, 16  80, 20  100, and 24  120, with special attention paid to the thermal bounder layer of the heating/cooling plate wall. Therefore, a fine structured mesh near the heating/cooling plate wall to resolve the thermal boundary layer. As appears from the results of these different grid sizes of Fig. 3, it can be seen that the differences of heat flux are rather small for the whole simulation process. The computational result of liquid fraction and solid fraction of the latent thermal energy system varies within 1% and 1.2%, respectively, as the number of grids increasing from 20  100 to 24  120. To save computing time with no penalty of poor accuracy, the grid system of 20  100 can be regarded as grid independent and sufficed for the numerical investigation purposes.

2.4. Numerical solution Because of the nonlinear variation of rPCM, mPCM and rair with temperature, no analytic solutions of Eqs. (3e5) are available. Therefore, the conservation of mass, momentum and energy equations expressed by Eqs. (3e5), together with appropriate boundary conditions and initial conditions, governing heat transfer in the latent thermal energy storage/release system, are solved using commercial CFD code, Fluent [38], based on finite volume method described by Patankar [39]. Based on the local value of an, the appropriate properties and variables will be assigned to each control volume within the domain. User-defined functions (UDF) are written in C language to account for temperature-depended of the thermo-physical properties of PCM. The coupling between pressure and velocity is conducted by the well-known SemiImplicit Pressure-Linked Equation (SIMPLE) algorithm. The convergence is checked at each time step, with the convergence criterion is that the scaled residuals are less than 104 for the continuity equation, 105 for the momentum equation and 108 for the energy equation. 2.5. Time-step and grid independence study The software Gambit is a geometric modeling and mesh generation tool. The model drawing are created and meshed by using Gambit, which is packaged with Fluent. As a first step of all

Fig. 2. Time-step dependence of the numerical solution for ΔT equal to 20
C (energy storage) and 25
C (energy release) for grid size 20  100.

3876

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

a + 200

0

+ 180

-40 Energy release

+ 160

-80 -120

+ 120 + 100

13 16 20 24

+ 80 + 60

-160

67 80 100 120

-200

q' (kW/m2)

2

q' (kW/m )

+ 140 Energy storage

-240

+ 40

-280

+ 20 0

60

80

100

120 t (s)

140

160

180

200

-320 220

1.0

1.0 0.9

Liquid fraction

40

Energy storage

0.9

0.8

0.8

0.7

0.7 0.6

0.6 Energy release

0.5

0.5 13 16 20 24

0.4 0.3

67 80 100 120

0.4 0.3

0.2

0.2

0.1

0.1

0.0 0

20

40

60

80

100 120 t (s)

140

160

180

Solid fraction

b

20

200

0.0 220

Fig. 3. Grid dependence of the numerical solution for ΔT ¼ 20
C, Δt ¼ 0.006 s (energy storage), and ΔT ¼ 25
C, Δt ¼ 0.003 s (energy release).

Additionally, three-dimensional model is also explored. No considerable difference has been found between threedimensional and two-dimensional model. A similar treatment method was also presented by Agyenim et al. [40], they realized a comparative study on latent heat storage using bare and finned tubes. They found to vindicate most models (two-dimensional symmetry) that ignore the thermal conductivity in the axial direction. In addition, Lacroix [41] and Ghoneim [42] developed the solutions of phase change problems with respect to ignoring the thermal conductivity or the temperature gradient of the PCM in the direction of the heat transfer fluid flow based on assumption that the phase change problems were two-dimensional symmetric. 2.6. Experimental validation The present numerical model validation is conducted by comparison of the numerical prediction with the experimental data proposed by Gau et al. [43] and the implicit finite difference results of Lacroix [44] melting of a pure metal (gallium) inside a twodimensional rectangular cavity with aspect ratio Hy/Wx ¼ 0.5. The top and bottom walls are insulated. At the beginning, the enclosure contained solid gallium slightly below its melting point, i.e., T0 ¼ 301.45 K. The temperature of the left vertical wall is suddenly raised to 311.15 K, while the right wall is maintained at 301.45 K. The values of the governing dimensionless numbers Rayleigh number Ra, Prandtl number Pr, and Stefan number Ste are equal to 2.2  105, 0.021, and 0.042, respectively.

The predicted phase interface position with both the experimental results from Ref. [43] and the finite difference prediction from Ref. [44] are compared, respectively, as shown in Fig.4a. It is noted that the results of the present model are in good agreement with the experimental data and Ref. [44], respectively. The discrepancy between the predicted phase interface position of the present model and the experimental results may be explained as follows. The solid PCM shows an initial sub-cooling of 274.63 K in the experiment. This degree of sub-cooling is significant in the light of the fact that the heating wall is only 281.15 K higher than the melting temperature of gallium. Furthermore, it is difficult to heat the vertical wall to a desired temperature in reality due to twodimensional simplicity. However, the effect of the z-direction is neglected in our numerical simulations, and thus leads to a deviation. The discrepancy of predicted phase interface position between the present model and Lacroix’s model may be attributed the difference of the numerical methods used. Lacroix used a front tracking method, while this model uses an enthalpy-porosity approach to model the phase change effects. Another test case investigated is transient natural convection from a finned vertical heating wall for thermal storage in threedimensional rectangular cavity. The experimental data proposed by Lacroix et al. [45]. The rectangular test cell, containing the PCM n-Octadecane with Tpcm ¼ 300.7 K, has inside dimensions of 0.60 m in z-direction, aspect ratio Hy/Wx ¼ 5.0, and horizontal equidistant length of fins (one fin and five fins are considered) equal to 0.02 m. A small air gap is maintained between the top of the PCM due to the volume expansion associated with the phase change process from the solid to the liquid. The right vertical wall is heated by a constant temperature (333 K) water bath. The other walls of the cavity are well insulated by a thick layer of Styrofoam insulation. The comparison of the calculated liquid volume fraction with experimental data by Ref. [45], reaching a very good agreement except for a slight deviation in the final stage, as shown in Fig.4b. The possible reasons can be attributed to numerical errors caused by long calculation periods. Furthermore, they do not take into account the volume expansion due to the phase change and temperature-dependent thermophysical in their model. However, the predicted results can be still acceptable as an engineering computation. Additionally, Zivkovic et al. [46] performed a comparison between the measured and predicted values of the temperature at the center of the three-dimensional rectangular container (0.1 m in length, 0.1 m in width, and 0.02 m in height) for the melting of CaCl2$6H2O. The mathematical model proposed by Ref. [46] was solved using finite volume method [39]. The agreement between experiment and numerical simulations is reasonable, as displayed in Fig.4b. However, the simulated curve is higher slope than that from experimental result. The reasons may be included in the follows: Firstly, the natural convection in the liquid PCM is ignored. Secondly, the convective heat transfer coefficient between the air and the container wall in this numerical solution (16 W/m2 K) is higher than the really. Thirdly, the conduction resistance of the container wall is neglected in mathematical model. It implies that the present numerical simulations are reasonable and the validity of the present model is verified. Therefore, the numerical model can be used successfully to predict the fluid flow and heat transfer performances of thermal energy storage/release system. 3. Results and discussion of the numerical simulation This section presents the detailed analyses fluid flow and heat transfer of the thermal energy storage and thermal energy release processes. Then, dimensional analyses in the process of thermal storage/release are performed.

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

a

3877

1.0 0.9 0.8 0.7

y/H

0.6 0.5 Experiment (Ref. [43])

0.4

5 min 8 min 12 min

0.3 0.2

Prediction (Ref. [44])

0.1

Simulation 0.0 0.0

0.2

0.3

0.4

0.5 x/W

0.6

0.7

0.8

0.9

330 1.0 0.9

325

0.8

320

0.7 Liquid fraction

1.0

315

0.6 310

0.5 0.4

T (K)

b

0.1

305

0.3

1 fin Experimental (Ref. [45]) 5 fins Experimental (Ref. [45]) Experimental (Ref. [46]) Simulations

0.2 0.1 0.0 0

10

20

30

40

50

60 70 80 t (min)

300 295

290 90 100 110 120 130 Fig. 5. Liquid fraction (a) and transient heat fluxes q’ (b) from the heating plate wall during energy storage.

Fig. 4. Comparison of the numerical predictions and experimental results.

3.1. The process of thermal energy storage A multiphase computational analysis is performed by numerical simulation in detail, including the thermal and flow characterization, solideliquid interface moving, interaction between the PCM and air, and phase distribution in PCM. Fig. 5 illustrates the calculated liquid fractions (also melting amount) and heat fluxes q’ on the thermal parameters of the system as a function of time. The curves are displayed for ΔT ¼ 5
C, 10
C, 15
C, 20
C, 25
C, 30
C, and 35
C. The x-coordinate is time (s). The y-coordinate in Fig. 5a is the liquid fraction. Non-melting to complete melting is expressed from 0 to 1.0. As expected, the larger is the temperature difference, ΔT, the more rapid is the growth of the liquid fraction, and shorter is the complete time. Consequently, the larger is the ΔT, the steeper in liquid fraction of increase. That’s to say, ΔT significantly dominates thermal storage speed. Since it has been found that the behavior of the system is similar for various temperature differences as shown in Fig. 5, only representative correlations are given. Therefore, the relation between liquid fraction and the thermal storage time (seconds), t, for case 6 can be written as:

Liquid fraction ¼ 0:0672 þ 0:0222t    0:0001t 2 R2 ¼ 0:9985

(8)

From Fig. 5b, one can be observed that transient heat flux q’, i.e., q’ ¼ kpcm (vT/vy), is considerably decreases at the early stages of

energy storage which indicative the transient heat conduction through the heating plate wall. The larger is the temperature difference, ΔT, the higher is the heat flux initially. Therefore, the complete time of thermal storage is shortening as mentioned in Fig. 5a. During this process, the surface heat flux decreases due to the increasing thermal resistance of the growing layer thickness of the melting medium. However, since buoyancy driven flow is sets in and develops, the decrease subsequently in heat transfer slows down and maintains for a long time. In the final stage, heat flux is followed by a decrease tendency. The increase of temperature differences leads to increase the temperature gradient between the heating plate wall and the PCM, and consequently increases heat flux. During thermal energy storage process, the relation between the transient heat flux q’ (kW/m2) and the melting period t (seconds) for case 5 can be written in the form as follows:

  q0 ¼ A0 þ A1et=a þ A2et=b R2 ¼ 0:9777

(9)

where A0 ¼ 44.9009  0.0538; A1 ¼ 200.9921  0.6777; a ¼ 2.3948  0.0104; A2 ¼ 638.6243  3.1028; b ¼ 0.0787  5.9353  104. A variety of numerical methods had been proposed for solideliquid interface problems [47], such as enthalpy method, moving meshes, coordinate transformations, the boundary immobilization, etc [48]. Nevertheless, the enthalpy method is used for present study as mentioned in Section 2.2.

3878

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

A typical evolution of the predicted PCM liquidesolid interface moving against time (see Fig. 6a) and the detailed liquid volume fraction/vector distribution (see Fig. 6b) for liquid fraction equal to 0.65 corresponding to every case are depicted in Fig. 6. This curve is obtained from the contour images of the interface during the simulation. Two parameters, x/W and y/H, are introduced for defining the interface position dimensionless. In Fig. 6a, one can see that the interface is nearly planar due to conduction during the early stage of melting (t ¼ 12 s). However, the effect of natural convection caused the interface shape to become oblique as the time goes. The instantaneous melting interface of cases 6, 7 are obviously more parallel to the fin wall than others, while are curved at the upper part of the interfaces for all cases, which indicated that all cases in the upper portion melted faster that in the lower portion. Additionally, the liquid phase starts to move inward due to the buoyancy force after a liquid layer is formed around the heating plate wall, i.e., the solideliquid interface moves unparallel away from the heat transfer surface. These phenomena differ distinctly from those derived from the thermal conduction theory but are in accordance with those dominated by the natural convection. Moreover, when the heating plate wall temperature is high, the temperature gradient between the working fluid and the PCM is big leading to high heat transfer rates and hence more PCM melted mass. Opposite effects are found when the heating plate wall

temperature is low. The curves also show longer melting intervals lead to more molten mass. Similar results are found for the other cases. It should be noted that in the figures the white line in Fig. 6b represents the interface between the air and PCM. The air above the liquid PCM (shows in red) is closed to stationary due to its lower density. Due to buoyancy, the liquid PCM near the fin and heated plate wall is driven to flow up. However, at the same time the gravity tries to drive the liquid PCM flowing downward from the top to the bottom. At last, the liquid PCM flows are driven by the stronger one. The weak circulation can be found in cases 1, 2 as the result of low temperature differences triggered the stronger buoyancy driven flows between the fins and phase interface. Ultimately, it is interesting to observe that the velocity vector distribution region decreases through cases 3e7. Furthermore, the profiles of the melt front for cases 1, 2 are different with cases 3e7, which display the similar solideliquid phase interface profile. As can be seen from cases 4e7, part of the heat is transferred through the melting from the bottom heating plate wall to the top cold phase front. This situation is unstable as layers of cold dense fluid adjacent to the solideliquid interface lie above layers of hot and lighter fluid near the bottom heating wall. As a result, circulating flow is eventually occurred near the bottom heating wall thereby increasing significantly the melting rate in these regions.

Fig. 6. The phase interface position (a) at various moment, and liquid fraction/vector distribution (b) for phase volume fraction equal to 0.65 for energy storage.

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

To measure how much heat flux is transferred through heating plate wall throughout the energy storage/release process, average heat transfer coefficient (Uavg) will provide further information about the influence of the various ΔT, is evaluated by:

Ztc Uavg ¼

q0avg

DT

¼

q0 dt

t0

tc  t0

=DT

(10)

where to and tc are the initial and complete time of energy storage, respectively. The above integral is taken over the entire surface of the heating wall. The dynamics of heat conduction upon phase changes is relatively complex due to the presence of a nonlinear moving boundary, and this class of problem is usually denoted Stefan-type [49]. Therefore, an important parameter has to be examined in latent thermal energy storage/release system about the energy density, Ste, which comprises both latent heat and sensible heat, is given by:

Ste ¼

cp DT L

(11)

i.e., it is based on specific heat cp, latent heat capacity of PCM L, and temperature difference ΔT. Fig. 7 plots the effect of Ste on the average heat transfer coefficient for all the cases. As can be seen from the curve, average heat transfer coefficient decreases rapidly with Ste for Ste  0.299, that is, the Ste number of case 4 is a critical value for heat transfer enhanced. However, the variation of the average heat transfer coefficient due to the Ste influence is small when Ste > 0.299. It remains near constant as Ste increases from 0.374 (case 5) to 0.523 (case 7), i.e., average heat transfer coefficient does not clearly vary when Ste up to 0.299. Therefore, it can be used to guide the designing and optimization of plate-fin unit for efficient thermal energy system. The variation of average heat transfer coefficient Uavg (W/m2K) with Ste can be predicted from the relation derived for all of cases as below:

Uavg ¼ 5:6537  30:4629Ste þ 82:6092Ste2    72:5821Ste3 R2 ¼ 0:9785

(12)

4.0 case 1

2

Uavg (kW/m K)

3.5

3.0 case 2 2.5

case 3 case 4

2.0 case 5 1.5 0.0

0.1

0.2

0.3 Ste

0.4

case 6

case 7

0.5

0.6

Fig. 7. The average heat transfer coefficient versus Stefan number for energy storage process.

3879

3.2. The process of thermal energy release The process of thermal energy release or solidification of PCM is based on after highly endothermic heat storage and combined with heat release upon exothermic solidification. The high enthalpy changes upon solidification inevitably requires proper understanding and description of the underlying heat transfer process, and to design efficient heat release appliances. The example of detailed phase fields throughout the energy release process, i.e., the solid shape (shows in blue) evolution for case 6, as shown in Fig. 8a. One can see that the liquid PCM sinks in the solidified PCM, part of not solidified PCM can be seen clearly and continued for a long time (t ¼ 374 s, solid fraction of 0.886) because of the amorphous structure of PCM. From the position of white line, which denote the interface between air and PCM. The volume of the PCM contracts obviously associated with solidifying between 2 and 38 s due to the energy release (the first four graphs as shown in Fig. 8a). However, the decreases in volume are slight after 38 s. In other words, the energy release is mainly focus on the first 38 s. The temperature distribution of PCM-air system at different moment during energy release for case 6 is given in Fig. 8b. From the five graphs, the temperature decreases from the left to the right side of symmetry and from the middle to the exterior. It is because that the temperature of liquid PCM is higher than solideliquid interface, the heat transfer by natural convection from liquid PCM to interface first. Then heat transfer from interface to solid PCM layer and fin by conduction, the temperature of solid PCM is lower than interface. That’s to say the overall kinetics becomes more dependent on heat conduction through the solid PCM as its thickness increase [50]. The representative velocity field of PCM-air system at the moment of complete thermal release in comparison with the end of thermal storage (namely at the beginning of energy release) for case 2 is depicted in Fig. 9, i.e., the time of energy release equal to 102 s (solid fraction of 0.883). It can be seen from Fig. 9a; a vortex is formed in the air region due to the natural convection of the air carried thermal energy from PCM. As a result, it makes large temperature differences between the top and the bottom PCM. However, the magnitude of velocity is very small in the region of air at the end of thermal storage. Comparison of Fig. 9a with Fig. 9b, when PCM solidifies and volume contracts, a gap at the end surface of solid PCM forms due to the negative pressure caused by volume contraction. Furthermore, part of not solidified PCM can be seen clearly in comparison with the process of energy storage (see Fig. 6b). The phenomena of not solidified PCM also similarly display in Fig. 8a. Fig. 10 plots a typical variation of solid fraction, transient heat flux q’ against time t. As can be seen from Fig. 10a the solid fraction increases rapidly at the initial time due to the small thermal resistance between the cooling plate wall and the PCM. Low temperature differences of melting lead to more solidified mass and short complete solidification time. The influence of the heat absorbed, during thermal storage process, can be seen clearly. Higher the absorbed heat during the thermal storage process, more dominant is the latent heat affect the thermal release process. As the time passes, more PCM is solidified, which causing the increase of the thermal resistance to heat transfer due to the growing layer of solid PCM, and consequently the decrease of the heat transfer rate as indicated by the continuously decreasing inclination of the curve. When the energy release process is trended to be completed, the solid fraction is infinitely closed to 0.890 not 1.0 due to existing part of discharging PCM as simulated in Figs. 8 and 9a. Therefore, the instant of time of thermal energy release is considered to be finished as long as solid fraction is converged to 0.883, such as

3880

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

Fig. 8. a) Example of detailed phase fields throughout the energy release process, for case 6. b) The temperature distribution of PCM-air system at various times during energy release process, for case 6.

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

3881

t ¼ 164 s for case 6, t ¼ 102 s for case 2. Similar treatments can be used for other cases. Towards the end of the thermal release process the solid fraction is nearly unchanged and heat flux nearly close to zero (see Fig. 10b), that is, thermal release is completed as indicated by the end on the curves. The results display a rapid decrease in heat flux at the early stages of energy release process which is indicative the transient heat convection as can be found by Fig. 10b. Due to the growing layer thickness of the solidified medium, heat flux decreases slowly. The major parts of energy release are focus on the first 38 s (approximate) and then heat flux will nearly keep constant after 38 s. These finds can be also qualitatively obtained as mentioned from Fig. 8a. The variation of solid fraction with the solidification time (seconds), t, can be predicted from the relation derived for case 4 as be given by:

solid fraction ¼ 0:037 þ 0:0268t  0:0003t 2   þ 106 t 3 R2 ¼ 0:9936

(13)

The relation between the transient heat flux (kW/m2), q’, and the solidification period (seconds), t, for case 4 can be written as:

  q0 ¼ A0 þ A1et=a þ A2et=b R2 ¼ 0:9902 Fig. 9. Examples of the flow field at the end of thermal release (a) in comparison with the end of thermal storage (b), for case 2.

(14)

where A0 ¼ 9.0084  0.0301; A1 ¼ 129.3904  0.1500; a ¼ 18.3937  0.0323; A2 ¼ 320.3457  4.1382; b ¼ 0.5775  0.0050. Fig. 11 displays the average heat transfer coefficient variations of all the simulated cases during energy release process. As can be seen from Fig. 11, the average heat transfer coefficient decreases when the cases vary from case 1 to case 7, i.e., temperature differences are increasing gradually. Nevertheless, the average heat transfer coefficient of case 1 is much larger than other cases. With certain dimensions of the thermal storage/release system the natural convection effect can increase the heat transfer coefficient. However, the natural convection effect is limited when the temperature differences are larger than 20
C, i.e., case 4 is a critical case for heat transfer enhanced or not. 3.3. Dimensional analysis Some relevant non-dimensional parameters are chosen, as reported in ref. [36] for Nu  Ste  Fo, liquid fraction  Ste  Fo and

5.0 4.5 4.0

2

Uavg (kW/m K)

3.5 3.0 2.5 2.0 1.5 1.0 0.5 Fig. 10. Solid fraction and heat fluxes q’ profiles during energy release: (a) solid fraction; (b) heat fluxes from the cooling plate wall.

case1

case2

case3

case4

case5

case6

case7

Fig. 11. The average heat transfer coefficient during energy release process.

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

solid fraction  Ste  Fo for cases selected. Heat transfer is described with the dimensionless Nusselt number, Nu, which is a function combined with various parameters, is defined as:

Nu ¼

Uavg  Wpcm Wpcm q0 ¼  DT 2kpcm 2kpcm

(15)

where Uavg is the average heat transfer coefficient of the plate-fin unit, Wpcm is the thickness of the PCM, i.e., Wpcm =2 is chosen as the characteristic length scale, kpcm denotes the thermal conductivity of the PCM, q’ is the mean heat flux from the plate wall. Nu is a function of time in dimensionless form, s, where dimensionless time is written as a product of the Stefan number, Ste, and Fourier number, Fo, and reads:

s ¼ Ste  Fo ¼

cp DT at 
2 L Wpcm =2

(16)

  NuðstorageÞ ¼ A0 þ A1eas R2 ¼ 0:780




A1 ^ NuðreleaseÞ ¼ A0 þ 0:7979   e2 sa A2 2 A2   R2 ¼ 0:9667

(18)

where A0 ¼ 4.6207  0.0273; A1 ¼ 1.059  1.3826; A2 ¼ 1.0929  0.0667; a ¼ 1.8361  0.3011. Fig. 13 plots liquid fraction and solid fraction for various dimensionless times for typical cases. From Fig. 13, it can be seen that both liquid fraction and solid fraction increase following the dimensionless time increase. As seen from Figs. 13, 5a and 10a, it is indicated that similar tendency is found for these simulated cases. It is higher of liquid fraction and lower of solid fraction in case 1 than case 4, 7, respectively, for a given dimensionless time. However, opposite effects are found form Figs. 5 and 10a. For case 4 liquid fraction and solid fraction increase with dimensionless time s increasing can be predicted from following relation, respectively, as below:

Liquid fraction ¼ 0:0460 þ 3:2589s    2:2419s2 R2 ¼ 0:9994

(19)

Solid fraction ¼ 3:3790 þ 9:8443s  7:7293s2   þ 2:0469s3 R2 ¼ 0:9920

(20)

The dimensionless mean temperature is used to evaluate the PCM-air system during energy release, is given by:

Q¼

Tavg  Tw Tpcm:avg  Tw

(21)

which combines the values of the mean temperature of PCM-air system during energy release Tavg, the mean phase change temperature of the PCM Tpcm, avg, and the heating/cooling plate wall temperature Tw.

1.0

1.0

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

case 1 case 4 case 7

0.4 0.3

0.4 0.3

0.2

0.2

0.1

0.1

0.0 0.0 Fig. 12. The transient Nusselt number, Nu, versus dimensionless time, for typical examples.

(17)

where A0 ¼ 23.4713  0.0166; A1 ¼ 87.1843  0.2195; a ¼ 36.9897  0.1106.

Liquid fraction

where a is thermal diffusivity coefficient, and the time t. The variation of the Nusselt number, Nu, with dimensionless time, s, for cases 1, 4 and 7 is presented in Fig. 12. Nu (storage) and Nu (release) denotes the Nu values of energy storage and release, respectively. For cases 4 and 7, Nu (storage) decreases as dimensionless time increases. Nu (storage) decreases for case 4 similar to case 7 due to large temperature differences of melting leads to melt rapidly, the PCM near the heating plate wall is in mushy zone and the temperature remains practically unchanged. However, a difference is observed between case 1 with a low temperature difference, and the other cases. In the latter, the Nu (storage) decreases nearly monotonically throughout the thermal storage process. In the former, an increase in the Nu (storage) is observed, but its duration is shorter. This increase may be attributed to motion initiated in the liquid phase. The motion is due to temperature-driven natural convection. Nu(release) of case 1 (low ΔT of melting) is higher than that of cases 4 and 7 for a given dimensionless time. After approximately dimensionless time of 1.0, a stable Nu (release) is reached for all cases due to the solidification of PCM are commonly controlled by heat conduction through the solid phase layer. Furthermore, Nu (release) trends in the end stage are similar with flat in the three curves, while significant difference in the initial stage with steep slopes that indicated the heat conduction rate is high with good heat release. In other words, parameter Ste remarkably affects the initial Nu (release).

The relation between Nu and dimensionless time s, for case 4 can be written as:

0.2

0.4

0.6

0.8 τ 1.0

1.2

Solid fraction

3882

1.4

1.6

0.0 1.8

Fig. 13. Liquid fraction, solid fraction versus the dimensionless time, for typical examples.

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

3883

3. The correlations are established as a function of associated variables. These relationships can be used for further system optimization to design systems to meet the required conditions for application plate-fin heat storage/release system in practical engineering applications, particularly those involving solideliquid phase change materials. Acknowledgements The authors are grateful to Dr. Simin Hung of South China University of Technology for his helpful discussion and valuable advice. This work was financially supported by Key Science and Technology Park Support Program of Guangzhou Municipal of China for this research project (No. 2008Z1-1061). The authors also take this opportunity to express sincere respect to the peer reviewers for their valuable comments and suggestions.

Fig. 14. Generalized results for all simulated case during energy release, the dimensionless average temperature of PCM-air system versus dimensionless time.

The effect of temperature differences on dimensionless average temperature of PCM-air system Q is given in Fig. 14. The figure shows that these curves have the similar shapes and tendencies i.e., the dimensionless mean temperature diminishes as dimensionless time is increasing. It can be seen from cases 1e7, the dimensionless mean temperature increases as ΔT increases at the same dimensionless time due to ΔT of energy storage increasing gradually. However, when the dimensionless time is larger than 0.4, the dimensionless average temperature decreases linearly. This is caused by the dimensionless time just covers over the critical point (dimensionless time of 0.4), as shown in Fig. 13. The relation between the dimensionless mean temperature Q and the dimensionless time s for case 4 can be predicted from the relation as:

Q ¼ 1:3046  2:7325s þ 3:9267s2    2:0528s3 R2 ¼ 0:9839

(22)

4. Conclusions A computational study has been conducted to investigate the processes of the latent thermal energy storage/release in a plate-fin unit with the uniform temperature on heating/cooling plate wall. The effect of temperature differences between the heating/cooling plate wall and the mean melting temperature of PCM on the performances of phase-change fluid flow and heat transfer are investigated in detail. Following results can be found: 1. During the process of energy storage, the weakened circulation can be found at the solideliquid interface as the result of lower temperature differences (10
C) triggered the stronger buoyancy driven flows. The performances of the energy storage are enhanced when temperature differences are less than 20
C. Further, circulating flow is occurred near the bottom heating wall for high temperatures difference. Hence, the melting rate increases remarkably in these regions. 2. During the process of thermal energy release, part of not solidified PCM can be observed clearly at the stage of energy release and continued for a long time. Therefore, results in the solid fraction less than 1.0. At last, a vortex appears in the air region at the moment of complete thermal release.

References [1] A. Barba, M. Spiga, Discharge mode for encapsulated PCMs in storage tanks, Solar Energy 74 (2) (2003) 141e148. [2] C.J. Ho, R. Viskanta, Inward solid-liquid phase-change heat transfer in a rectangular cavity with conducting vertical walls, International Journal of Heat Mass Transfer 27 (7) (1984) 1055e1065. [3] A. Gadgil, D. Gobin, Analysis of two-dimensional melting in rectangular enclosures in presence of convection, Journal of Heat Transfer 106 (1) (1984) 20e26. [4] M. Okada, Analysis of heat transfer during melting from a vertical wall, International Journal of Heat and Mass Transfer 27 (11) (1984) 2057e2066. [5] H. Inaba, C. Dai, A. Horibe, Natural convection heat transfer of microemulsion phase-change-material slurry in rectangular cavities heated from below and cooled from above, International Journal of Heat Mass Transfer 46 (23) (2003) 4427e4438. [6] D. Pal, Y.K. Joshi, Melting in a side heated tall enclosure by a uniformly dissipating heat source, International Journal of Heat Mass Transfer 44 (2) (2001) 375e387. [7] H. Ettouney, I. Alatiqi, M. Al-Sahali, K. Al-Hajirie, Heat transfer enhancement in energy storage in spherical capsules filled with paraffin wax and metal beads, Energy Conversion and Management 47 (2) (2006) 211e228. [8] D. Zhou, C.Y. Zhao, Experimental investigations on heat transfer in phase change materials (PCMs) embedded in porous materials, Applied Thermal Engineering 31 (5) (2011) 970e977. [9] K.A.R. Ismail, F.A.M. Lino, Fins and turbulence promoters for heat transfer enhancement in latent heat storage systems, Experimental Thermal and Fluid Science 35 (6) (2011) 1010e1018. [10] A. Sari, A. Karaipekli, Thermal conductivity and latent heat thermal energy storage characteristics of paraffin/expanded graphite composite as phase change material, Applied Thermal Engineering 27 (8e9) (2007) 1271e1277. [11] A. Al-Hallaj, M. Mills, J.R. Farid, S. Selman, Thermal conductivity enhancement of phase change materials using a graphite matrix, Applied Thermal Engineering 26 (14e15) (2006) 1652e1661. [12] F. Frusteri, V. Leonardi, S. Vasta, G. Restuccia, Thermal conductivity measurement of a PCM based storage system containing carbon fibers, Applied Thermal Engineering 25 (11e12) (2005) 1623e1633. [13] F. Agyenim, N. Hewitt, P. Eames, M. Smyth, A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (LHTESS), Renewable and Sustainable Energy Reviews 14 (2) (2010) 615e628. [14] R. Velraj, R.V. Seeniraj, B. Hafner, C. Faber, K. Schwarzer, Heat transfer enhancement in a latent heat storage system, Solar Energy 65 (3) (1999) 171e180. [15] S. Jegadheeswaran, S.D. Pohekar, Performance enhancement in latent heat thermal storage system: a review, Renewable and Sustainable Energy Reviews 13 (9) (2009) 2225e2244. [16] B. Zalba, J.M. Marin, L.F. Cabeza, H. Mehling, Review on thermal energy storage with phase change: materials, heat transfer analysis and applications, Applied Thermal Engineering 23 (3) (2003) 251e283. [17] L.F. Cabeza, H. Mehling, S. Hiebler, F. Ziegler, Heat transfer enhancement in water when used as PCM in thermal energy storage, Applied Thermal Engineering 22 (10) (2002) 1141e1151. [18] M. Medrano, M.O. Yilmaz, M. Nogues, I. Martorell, J. Roca, L.F. Cabeza, Experimental evaluation of commercial heat exchangers for use as PCM thermal storage systems, Applied Energy 86 (10) (2009) 2047e2055. [19] M. Akgun, O. Aydin, K. Kaygusuz, Thermal energy storage behavior of a paraffin during melting and solidification, Energy Sources Part A-Recovery Utilization and Environmental Effects 29 (14) (2007) 1315e1326. [20] M. Akgun, O. Aydin, K. Kaygusuz, Experimental study on melting/solidification characteristics of a paraffin as PCM, Energy Conversion and Management 48 (2) (2007) 669e678.

3884

W.-B. Ye et al. / Applied Thermal Engineering 31 (2011) 3871e3884

[21] U. Stritih, An experimental study of enhanced heat transfer in rectangular PCM thermal storage, International Journal of Heat Mass Transfer 47 (12e13) (2004) 2841e2847. [22] K.A.R. Ismail, J.R. Henriquez, T.M.D. Silva, A parametric study on ice formation inside a spherical capsule, International Journal of Thermal Science 42 (9) (2003) 881e887. [23] K.A.R. Ismail, J.R. Henriquez, Numerical and experimental study of spherical capsules packed bed latent heat storage system, Applied Thermal Engineering 22 (15) (2002) 1705e1716. [24] D. Buddhi, S.D. Sharma, A. Sharma, Thermal performance evaluation of a latent heat storage unit for late evening cooking in a solar cooker having three reflectors, Energy Conversion and Management 44 (6) (2003) 809e817. [25] M. Akgun, O. Aydin, K. Kaygusuz, Thermal energy storage performance of paraffin in a novel tube-in-shell system, Applied Thermal Engineering 28 (5e6) (2008) 405e413. [26] J.-Y. Long, D.-S. Zhu, Numerical and experimental study on heat pump water heater with PCM for thermal storage, Energy and Buildings 40 (4) (2008) 666e672. [27] S.-Y. Wu, Enhanced heat transfer experimental and simulation research of nanocomposite phase change materials, Doctor of Philosophy Thesis, South China University of Technology, Guangzhou, China, 2010. [28] W.R. Humphries, E.I. Griggs, A design handbook for phase change thermal control and energy storage devices. Technical Report, 1074NASA Scientific and Technical Information Office, 1977. [29] R.C. Reid, J.M. Prausnitz, B.E. Poling, The properties of gases and liquids. McGraw-Hill, New York, 1987. [30] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics 39 (1) (1981) 201e225. [31] K. Ghosh, A. Mukhopadhyay, S. Sen, D. Sanyal, A sphericosymmetric VOF approach for investigating immiscible two-phase systems with one liquid phase, Numerical Heat Transfer, Part A: Applications 50 (10) (2006) 949e974. [32] V.V. Ranade, in: , Computational Flow Modeling for Chemical Reactor Engineering, vol. 5, Academic Press, St. Louis, MO, USA, 2002. [33] V.R. Voller, C. Prakash, . A fixed grid numerical modeling methodology for convection-diffusion mushy region phase-change problems, International Journal of Heat Mass Transfer 30 (8) (1987) 1709e1719. [34] A.D. Brent, V.R. Voller, K.J. Reid, Enthalpy-porosity technique for modeling convection-diffusion phase change: application to the melting of a pure metal, Numerical Heat Transfer 13 (3) (1988) 297e318.

[35] E. Assis, L. Katsman, G. Ziskind, R. Letan, Numerical and experimental study of melting in a spherical shell, International Journal of Heat and Mass Transfer 50 (9e10) (2007) 1790e1804. [36] V. Shatikian, G. Ziskind, R. Letan, Numerical investigation of a PCM-based heat sink with internal fins, International Journal of Heat and Mass Transfer 48 (17) (2005) 3689e3706. [37] H. Shmueli, G. Ziskind, R. Letan, Melting in a vertical cylindrical tube: numerical investigation and comparison with experiments, International Journal of Heat and Mass Transfer 53 (19e20) (2010) 4082e4091. [38] Fluent 6 User’s Guide. Fluent Inc, New Hampshire, USA, 2003. [39] S.V. Patankar, Numerical heat transfer and fluid flow. McGraw Hill, New York, 1980. [40] F. Agyenim, P. Eames, M. Smyth, A comparison of heat transfer enhancement in a medium temperature thermal energy storage heat exchanger using fins, Solar Energy 83 (9) (2009) 1509e1520. [41] M. Lacroix, Numerical simulation of a shell-and-tube latent heat thermal energy storage unit, Solar Energy 50 (4) (1993) 357e367. [42] A.A. Ghoneim, Comparison of theoretical models of phase-change and sensible heat storage for air and water-based solar heating systems, Solar Energy 42 (3) (1989) 209e220. [43] C. Gau, R. Viskanta, Melting and solidification of a pure metal on a vertical wall, Journal of Heat Transfer 108 (1) (1986) 174e181. [44] M. Lacroix, Predictions of natural convection-dominated phase change problems by the vorticity-velocity formulation of the Navier-Stokes equations, numerical heat transfer, Part B: Fundamentals 22 (1) (1992) 79e93. [45] M. Lacroix, M. Benmadda, Numerical simulation of natural convection dominated melting and solidification from a finned vertical wall, numerical heat transfer, Part A: Applications 31 (1) (1997) 71e86. [46] B. Zivkovic, I. Fujii, An analysis of isothermal phase change of phase change material within rectangular and cylindrical containers, Solar Energy 70 (1) (2001) 51e61. [47] N. Shamsundar, E. Rooz, Numerical methods for moving boundary problems, Handbook of Numerical Heat Transfer. Wiley-Interscience, 1988. [48] J. Caldwell, Y.Y. Kwan, A brief review of several numerical methods for onedimensional Stefan problems, Thermal Science 13 (2) (2009) 61e72. [49] J. Stefan. Sber, Akad. Wiss. Wien. 79 (1879) 161; J. Stefan, Sber. Akad. Wiss. Wien. 98 (1890) 965. [50] N. Vitorino, J.C.C. Abrantes, J.R. Frade, Numerical solutions for mixed controlled solidification of phase change materials, International Journal of Heat and Mass Transfer 53 (2010) 5335e5342.