Numerical solutions for mixed controlled solidification of phase change materials

International Journal of Heat and Mass Transfer 53 (2010) 5335–5342

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Numerical solutions for mixed controlled solidification of phase change materials N. Vitorino a, J.C.C. Abrantes a,b,*, J.R. Frade a a b

Ceramics Dep., (CICECO), University of Aveiro, 3810 Aveiro, Portugal UIDM, ESTG, Polytechnic Institute of Viana do Castelo, 4900 Viana do Castelo, Portugal

a r t i c l e

i n f o

Article history: Available online 9 August 2010 Keywords: Phase change materials Heat storage Mixed control Numerical method Moving boundaries

a b s t r a c t Kinetics of solidification of phase change materials (PCM) was analyzed for combined heat conduction in the PCM and container wall, and convection in the cold fluid. Stable and convergent numerical methods were derived after transformation to normalized size scale, and corresponding immobilization of the moving boundary. The accuracy was confirmed by comparing numerical solutions and corresponding analytical solutions for control by heat conduction in solidifying layers. The proposed method was used to assess when solidification is controlled solely by conduction in solidified layer, and to analyze relative roles of conduction through the container and convection in the cold fluid. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Solidification is very important in classical materials processing, with emphasis on metallurgy of metals and alloys, purification by zone melting, etc. More recently, one finds extensive literature on heat storage based on phase change materials [1–6]. This heat storage is based on highly endothermic melting, combined with heat discharge upon exothermic solidification. The high enthalpy changes upon melting/solidification inevitably requires proper understanding and description of the underlying heat transfer processes, to optimize materials solidification processes, and to design efficient heat storage appliances. Most prospective applications of phase change materials are for melting temperatures in the range from 0 to 100 °C, as required for domestic heating/air conditioning or hot water appliances, and for temperature regulation to prevent overheating or overcooling of food, electronic equipments, etc. In most of these cases, the selection of suitable phase change materials and appliance designs must take into account multiple requirements which include affordable price and availability, high heat storage density (per unit mass or per unit volume), reversibility of melting/solidification with reduced overheating/overcooling temperature differences, suitable heat transfer parameters, stability and inertness towards prospective container materials, low volume changes upon melting/cooling, etc. It is widely recognized that most phase change materials for near room temperature applications possess limited thermal con-

* Corresponding author at: Ceramics Dep., (CICECO), University of Aveiro, 3810 Aveiro, Portugal. Tel.: +351 258 819 700; fax: +351 258 827 636. E-mail addresses: [email protected] (N. Vitorino), [email protected] (J.C.C. Abrantes), [email protected] (J.R. Frade). 0017-9310/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2010.07.023

ductivity. The dynamics of heat conduction upon phase changes is relatively complex due to the presence of a moving boundary, and this class of problems is usually denoted Stefan-type, to recognize his pioneering contributions [7]. A variety of Stefan-type processes have been dealt with in the literature [8–19], with emphasis on heat or diffusion controlled processes [8–10], including solidification of phase change materials or alloys, solid state reactions [11], etc. By resorting to dimensionless analysis one recognizes that different types of Stefan problems are parametrically identical, which allows one to share methods developed independently for different technologies with heat conduction or diffusion control. Analytical solutions are available only for particular Stefan problems, with emphasis on planar geometry [8,9]. Alternative numerical or approximate solutions are needed for most other cases, thus requiring stable and converging computing codes, and criteria to assess the convergence of these numerical or alternative approximate solutions. Indeed, exact solutions available for particular cases may be very useful to establish those convergence criteria. A variety of numerical methods has been proposed for moving boundary problems [12], including the enthalpy method [13], boundary immobilization [14], and other approaches based on variable grids to avoid instabilities and enhance the convergence of numerical methods. Some methods adjust the time step in order to move the solid/liquid interface through predefined grid positions [15,16]. Other methods rely on moving meshes [17,18], coordinate transformations [19], etc. [20]. The increasing importance of energy managements and use of solar energy attracted renewed interest in phase change materials, mainly concerning optimization of heat conduction and storage in single reservoirs with planar, cylindrical or spherical configuration. However, prospective appliances rely on arrays of small PCM

5336

N. Vitorino et al. / International Journal of Heat and Mass Transfer 53 (2010) 5335–5342

Nomenclature cp cp,w h h* k kw  kw L L* T Tc Tm T W x X

heat capacity of PCM (J/(kg °C)) heat capacity of wall material (J/(kg °C)) heat transfer coefficient (W/(m2 K)) relative contribution in heat transfer of convection conduction thermal conductivity of PCM (W/(m K)) thermal conductivity of wall material (W/(m K)) relation between thermal conductivity of wall material and PCM wall thickness (m) relation between wall reservoir thickness and maximum thickness of PCM layer temperature (°C) temperature of cold fluid (°C) melting temperature of PCM (°C) time (s) dimensionless thickness distance from the inner wall of the reservoir (m) thickness of the solid phase layer (m)

reservoirs, often with spherical shape, for fast heat charge/discharge. Under these conditions, the complexity of such systems may require simpler solutions for the required geometries, often based on quasi steady state approximations [21]. The common assumption that solidification is controlled by heat conduction through the solid phase change layer is mostly correct for organic or inorganic materials with high melting enthalpy and relatively low thermal conductivity, separated from the cooling fluid by a good mixed conducting wall. However, this may not be the case if one considers the prospective use of phase change materials with significantly high thermal conductivity (e.g. metallic materials with low melting temperatures), and/or the need to use corrosion resistant container materials with lower thermal conductivity. The actual work is a development of a numerical method for mixed control, which accounts for the contributions of heat conduction through the growing layer of solid phase change material, heat conduction through the container wall with fixed thickness, and heat transfer from the container wall to the cooling fluid. Similar mixed control methods for solid state reactions may be found in the literature [22].

Xmax z

maximum thickness of PCM solid layer (m) dimensionless thickness of wall

Greek symbols a thermal diffusivity of PCM (m2/s) relation between thermal diffusivity of wall material a* and PCM aw thermal diffusivity of wall material (m2/s) b rate constant n fraction of solid layer k latent heat of PCM (J/kg) q density of PCM (kg/m3) qw density of wall material (kg/m3) / relative contribution between heat capacity and latent heat h dimensionless temperature s dimensionless time

– The heat exchanging reservoir wall separates the solid PCM from a cold fluid and heat transfer is controlled by conduction through the solid PCM layer, combined with conduction through the heat exchanging wall and heat transfer to the cold fluid. – The temperature at the reservoir wall/PCM interface is established by continuity of heat conduction fluxes in the solid PCM layer and reservoir wall. – The temperature at the wall/cold fluid interface is established by continuity of heat conduction in the wall and heat transfer in the cold fluid. Under these circumstances, the heat balance in the solid PCM layer can be expressed by:

@2T a @x2

!

 ¼

@T @t

 ð1Þ

Solid

2. Formulation of physical models

Solid

2.1. Control by heat conduction through the PCM solid layer Fig. 1 shows a schematic representation of solidification controlled by heat transfer to a cold fluid. The following assumptions are used to define a model behaviour for heat transfer in a solid PCM layer, formed upon highly exothermic solidification, and for a planar reservoir: – – – –

The PCM phases are contained inside a planar reservoir. Heat is being discharged from both surfaces. Edge effects are neglected. Symmetry is assumed throughout the entire process, corresponding to series association of solid–liquid–solid PCM layers inside the planar reservoir. – Solid–liquid interfaces are nearly flat. – Under-cooling is negligible, implying that the solid–liquid interfaces remain close to the melting temperature of the PCM. – The inner liquid PCM layer remains close to melting temperature, with negligible overheating.

Liquid

Cold fluid

Cold fluid

L

X X max

Fig. 1. Schematic representation of a PCM layer undergoing solidification by cooling under symmetrical ideal conditions.

5337

N. Vitorino et al. / International Journal of Heat and Mass Transfer 53 (2010) 5335–5342

where T is temperature, x is distance from the inner wall of the reservoir, t is time, and a is thermal diffusivity

a¼

k

qc p

ð2Þ

which combines the values of thermal conductivity k, density q and specific heat cp of PCM. Solutions for heat balance described by Eq. (1) must take into account the inner moving boundary, due to solidification. Under the assumptions formulated above, the moving boundary conditions are set by the rate of solidification, controlled by heat conduction from the solid/liquid interface. In other words, the heat flux tends to the following condition at the solid/liquid boundary (x = X):

qk

    @X @T ¼k @t @x X

ð3Þ

where X is the thickness of the solid phase layer and k is the latent heat of the PCM. One of the boundary conditions was based on the assumption that the solid/liquid interface remains at melting temperature, T(X, t) = Tm. The other condition depends on conditions at the wall/PCM interface and/or at the wall/cold fluid interface. The simplest condition is for cases when the thermal diffusivity of the reservoir wall material is much higher than for the PCM material, and for sufficiently high heat transfer coefficient (h) (i.e. h  k/X), in which case T(0, t)  Tc. Under these boundary conditions one can find relatively simple analytical solutions for this type of Stefan problems [9–11], and the thickness of the solidifying layer varies as: 1=2 X=X max ¼ b  X 1 max ðatÞ

ð4Þ

where the maximum thickness Xmax is attained when the PCM reaches complete solidification, and the rate constant b is given by [9]:

u¼

c p DT ¼ p1=2 b expðb2 Þerf ðbÞ k

ð5Þ

with DT = Tm  Tc. This dependence was used to obtain the time dependence of X/Xmax shown in Fig. 2 (continuous line), and to com-

pare these analytical solutions with corresponding numerical results obtained by a finite difference method (described below). The dimensionless parameter / = cpDT/k is a measure of the relative importance of sensible heat, associated to temperature changes in the solid PCM layer, relative to latent heat. In other words, high / values indicate that sensible heat (cpDT) is much higher than latent heat (k). In the opposite case, low / values means that latent heat prevails. For planar geometry one also finds analytical solutions for the relevant temperature or concentration profiles in Stefan problems, as follows [9]:

T  Tc erf ðx=ð4atÞ1=2 Þ ¼ T m  T c erf ðX=ð4atÞ1=2 Þ

ð6Þ

This is shown in Fig. 3, for different values of b = X/(at)1/2. For sufficiently low values of b, these solutions converge to the quasi steady state linear dependence given by (T  Tc)/(Tm  Tc)  (x/X), as predicted on solving Eq. (1) for oT/ox  0. On combining the corresponding quasi steady state gradient oT/ox  (Tm/Tc)/X with Eq. (3) one also obtains b  (2/)1/2, which is the limiting trend shown in Fig. 2. This shows that quasi steady state solutions are nearly correct for cpDT/k  1, i.e. when latent heat is much higher than sensible heat. This assumption is reasonably correct for many potential phase change materials, as shown in Table 1. Phase change materials may have potential applications requiring thermal inertia. In this case, one can establish a performance criterion based on the limiting trend shown in Fig. 3. Note that the time scale is given by t  0.5X2/(a/) = 0.5X2qk/(kDT), thus showing that qk/k is a performance criterion for prospective applications based on thermal inertia, i.e. high qk/k values minimize the impact of temperature fluctuations by maximizing the latent heat per unit volume, and minimizing thermal conductivity. Actually, the opposite limitation is expected in prospective applications requiring fast response time, i.e. quick discharge of latent heat, as for hot water appliances, requiring sufficiently high conductivity or higher temperature gradients. Thus, thermal conductivity limitations are often discussed in the literature [23–25]. 2.2. Heat transfer to a cooling fluid through the heat exchanging wall The container wall may also exert some effects on the kinetics of solidification. This can be analyzed by series association of com-

0.5

log[X/(α α t) 1/2]

β ≈(2φ )0,5 0

-0.5

-1 -2 Fig. 2. Time dependence of the thickness of solidified layer computed for conditions when the process is controlled by heat conduction through the solidified PCM layer (see inserted labels). Symbols represent finite difference predictions and solid lines are the corresponding analytic solutions (Eqs. (4) and (5)).

-1

log ( φ )

0

1

Fig. 3. Predictions of rate constant b = X/(at) obtained from numerical (symbols) and analytical solutions (solid line), computed for the conditions shown in Fig. 2.

5338

N. Vitorino et al. / International Journal of Heat and Mass Transfer 53 (2010) 5335–5342

Table 1 Relevant physical properties of potential PCM materials. The values of / = cPDT/k were obtained a temperature difference DT = Tm  15. Tf (°C) Parafin 1 Parafin 2 Parafin 3 Naftalene Stearic acid Palmític acid Lauric acid Na2SO410H2O Na3PO410H2O TH29 TH 21 Aluminum Bismuth Cadmium Sodium Lead Zinc

52 67 59 80 69 64 42 33 70 29 21 660 271 321 98 327 419

k (kJ/kg) 244 189 189 148 202 180 212 254 184 190 150 395 50 55 115 25 102

qe (MJ/m3) 197 176 174 169 172 153 212 377 331 325 222 672 492 470 111 280 730

q (g/cm3) 0.81 0.93 0.92 1.14 0.85 0.85 1.00 1.48 1.80 1.71 1.48 1.70 9.80 8.65 9.70 11.34 7.14

cp (kJ kg1 K1) 2.4 2.6 2.6 2.6 1.7 0.7 1.8 3.5 3.1 1.4 0.7 0.9 0.1 0.2 1.2 0.1 0.4

bined heat conduction in the heat exchanging wall, with constant thickness, and in the growing solidified layer. Heat conduction through the heat exchanging wall can be described as:

aw

@2T @x2

! ¼

  @T @t

ð7Þ

where

aw ¼

kw

qw cp;w

ð8Þ

is the thermal diffusivity of the wall material, which is related to the corresponding thermal conductivity kw density qw, and specific heat cp,w of the wall material. The thermal balance through the wall applies for L < x < 0, where L denotes the wall thickness, where the origin (x = 0) is located at the wall/PCM interface. The relevant relations between the dynamics of heat conduction in the exchanging wall and solidified layer can be set by assuming continuity of temperature and heat flux at the interface, i.e. T(L, t)  T(0, t) and:

    @T @T kw ¼k @x L @x 0

ð9Þ

One must also take into account heat transfer to the cooling fluid, outside the heat exchanging wall. Upon assuming continuity of the heat flux:

kw

  @T ¼ h½TðL; tÞ  T c  @x L

k (W m1 K1)

ð10Þ

3. Numerical solutions

0.15 0.21 0.21 0.31 0.17 0.16 1.60 0.54 0.60 1.09 0.43 222 7.04 95 141 34 108

8

8 10 8.7 108 8.8 108 1.1 107 1.2 107 3 107 9 107 107 1.1 107 4.6 107 4.3 107 1.44 104 5.9 106 4.8 105 1.19 104 2.3 105 3.9 105

/ 0.364 0.715 0.605 1.142 0.454 0.191 0.229 0.248 0.927 0.103 0.028 1.470 0.512 1.113 0.866 1.25 1.58

qk/k (sK/m2) 7

3.43 10 1.61 107 1.88 107 0.80 107 1.83 107 1.75 107 4.8 106 4.0 107 9.8 106 2.11 107 8.3 107 4.7 103 3.3 105 1.89 104 9.7 103 3.5 104 1.62 104

Ref. [27] [1] [1] [1] [1] [1] [28] [1] [1] [29] [29] [26] [26] [26] [26] [26] [26]

Boundary immobilization was achieved by the following change:

w¼

x XðtÞ

ð11Þ

thereby normalizing the size scale 0 6 w 6 1. This implies transformation of the differential equations by replacing the dependence on x and t, by the dependence on w and t. In addition, the relevant equations were rewritten in dimensionless terms, to minimize the number of truly independent parameters; this transforms Eq. (1) as follows:

!       @2h @h @h 2 @h ¼ n þ w/ @w2 @w 1 @w @s     @n / @h ¼ @s n @w 1

ð12Þ ð13Þ

where w = x/X(t) and

h ¼ ðT  T c Þ=DT

ð14Þ

/ ¼ cp DT=k n ¼ X=X max ta s¼ 2 X max

ð15Þ ð16Þ ð17Þ

with DT = (Tm  Tc), and maximum thickness Xmax for the growing solid slab, as illustrated in Fig. 1; this determines the actual time scale t max ¼ smax X 2max =a. If one assumes that the thermal diffusivity of the wall material is much higher than for the PCM layer, and for relatively high heat transfer coefficient to the cooling fluid (i.e. h  k/X), the boundary conditions reduce to h(0, s) = 0 and h(1, s) = 1. Otherwise, heat conduction through the heat exchanging wall with fixed thickness L can be written:

3.1. Immobilization of moving boundaries and dimensionless analysis

a Moving boundaries usually cause instabilities and poor convergence of numerical methods. Our approach to circumvent these difficulties is based on boundary immobilization by changing the actual size scale to a corresponding normalized scale, and performing the corresponding transformation of the relevant differential equations. The stability and convergence of this method can be demonstrated by comparing the numerical results to known analytical solutions available for the simplest case (Eqs. (4)–(6)). The convergence of numerical methods is also tested by varying the effects of mesh sizes, and analyzing their effects on computed results.

a  108 (m2/s)

!   @2h @h ¼ @z2 @s

ð18Þ

for L* 6 z 6 0, where

a ¼

a  w

a

ð19Þ

z ¼ x=X max

ð20Þ

L ¼ L=X max

ð21Þ

The temperature at the wall/PCM interface h(0, s) is given by continuity of heat flux, kw(oT/ox)d = k(oT/ox)+d and thus:

5339

N. Vitorino et al. / International Journal of Heat and Mass Transfer 53 (2010) 5335–5342 

nkw



   @h @h ¼ @z L @w o



kw ¼ kw =k

ð22Þ ð23Þ

Temperature at the cold fluid/wall interface may also be established by assuming conservation of the heat flux kw(oT/ ox)L = h[T(L, t)/Tc] and thus:

  @h   ¼ ðh =k ÞhðL ; sÞ @z L 

h ¼ X max h=k

ð24Þ ð25Þ

One thus used these relations to develop numerical algorithms for series association of heat conduction through the growing PCM layer, heat conduction through the reservoir wall and heat transfer into the cold fluid, as detailed in Appendix. 3.2. Stability and convergence

Fig. 4. Temperature profiles computed numerically for the inserted conditions (symbols), and corresponding analytical solutions (Eq. (6)).

tance contribution of the metallic container wall can be neglected even for cases when the thermal conductivity and thermal diffusivity values of the PCM and container material are similar (i.e. kw/k = 1 and aw/a = 1), provided that the thickness of the wall remains much thinner than the PCM layer (i.e. L*  1). Thus, one may assume that solidification is solely controlled by conduction through the PCM  layer for conditions when kw =L P 102 , and for sufficiently high heat transfer coefficient in the cold fluid. Some solidification processes may still be affected by poorly conducting container materials. For example, this may be the case if one considers the combination of low melting metallic PCM materials within ceramic containers. In these cases, one may expect that kw/k  1 and aw/a  1, and kinetics solidification will be dependent on heat conduction through the ceramic wall.

1

0.8

0.6

X/Xmax

The proposed algorithms yield very stable computing for wide ranges of parameters without undue computing requirements. The stability is not lost on varying the mesh sizes or the relative scale of time increments. One thus performed typical convergence tests by varying the mesh sizes (dw and dz), and also the relative time step (ds)/s. Convergence is easily recognized by verifying that predictions of time scale required to reach complete solidification (i.e. X/Xmax = 1) become nearly insensitive to changes in mesh sizes for relatively coarse mesh sizes. Typical examples were computed for spatial mesh sizes ranging from 0.25% (dw = 0.0025 and dz/ L* = 0.0025) to 1% (dw = 0.04 and dz/L* = 0.01) while retaining the mesh time (ds)/s = 0.005 (i.e. 0.5%) yielding differences in total time scale for complete solidification in the order of 0.1–0.2%. Similar or even better convergence was found on varying the mesh time. An even more conclusive demonstration of convergence is obtained by comparing numerical predictions and exact analytical solutions for conditions when the overall resistance to heat transfer is solely controlled by conduction through the solid PCM layer, i.e. for very large differences in thermal conductivity between the wall and PCM materials, and high heat transfer coefficient in the cold fluid. Under these conditions, one can use analytical solutions (Eqs. (4), (5), and (7)) to demonstrate the accuracy of our numerical method, as shown in Figs. 2–4 for aw/a = 103; kw/k = 103, (L/Xmax) = 0.02, and h/(k/Xmax) = 103 (Figs. 2–4). Fig. 2 clearly shows the expected parabolic time dependence for the thickness of solidified layer and close agreement between the finite difference results (symbols) and analytical solutions given by Eqs. (4) and (5) (solid lines). The slope of results shown in Fig. 2 was also used to obtain the values of rate constant b = X/(at)1/2, as shown in Fig. 3. Note that these solutions converge to a rather simple solution b = (2/)1/2, as expected for the limiting trend of Eq. (5) when b  1. In addition, Fig. 4 shows that the numerical predictions of temperature profiles are also nearly indistinguishable from the analytical solutions given by Eq. (6).

0.02 0.4

4. Mixed controlled kinetics

0.2 One may be tempted to ignore the contribution of the container wall material to the overall heat transfer resistance which controls solidification. This is indeed the case if one considers solidification of typical organic or inorganic non metallic PCM materials, and compares their thermal properties to those of copper or other metallic heat exchange materials. Note that the thermal conductivity and thermal diffusivity values of metallic heat exchange materials are about 3–4 orders of magnitude higher than for those PCMs (Table 1). Actually numerical solutions (Fig. 5) show that the thermal resis-

α* =kw* = 103; 1; 0.2; 0.1; 0.04; 0.02 L*=0,02; h * =103; φ=1

0 0

0.4

0.8

1.2

(tα α φ)1/2/Xmax Fig. 5. Predictions of the role of container wall for cases when the thermal properties of the container material may be inferior to those of the phase change material. The simulated conditions are shown in the insert.

5340

N. Vitorino et al. / International Journal of Heat and Mass Transfer 53 (2010) 5335–5342

Fig. 6. Predictions for conditions when the container wall acts as isolation, to delay heat losses and solidification. The simulated conditions are shown in the insert.

Fig. 8. Temperature profiles illustrating the mixed controlled behaviour for conditions inserted in the figure.

Actual values of heat transfer coefficient depend on the type of cooling fluid (liquid or gas) and its flow dynamics. Typical cooling fluids in heat storage applications are air or water, circulating under forced or even natural convection. If one considers latent heat storage and discharge in heat hot water appliances, one may expect typical values of heat transfer coefficient in the order of 103 W m2, and on combining this with a typical value for thermal conductivity of organic PCM materials (k  0.3 W m1 K1), and typical values of reservoir thickness in the range of 1 mm to 1 cm, one expects 0.3 < k/(Xmaxh) < 0.03. Thus, one cannot assume that the rate of heat transfer is always controlled by the poorly conducting PCM. This is even more likely to occur when latent heat is being transferred to air (e.g. heat storage for central heating), because typical values of heat transfer coefficient in air are often below 102 W m2, and thus quite lower than in water. During heat discharging one must also take into account that the relative roles of heat conduction through the PCM layer and heat convection in the adjacent fluid change with discharging time, i.e. the overall kinetics becomes more dependent on heat conduction through the solid phase change material as its thickness increases. This is illustrated by the simulations in Fig. 8. Fig. 7. Predictions of the role of heat transfer in the cold fluid for the conditions shown in the insert

5. Conclusions

Thermal isolation corresponds to the most extreme case of control by the thermal properties of the wall materials. This is, indeed, needed for long term heat storage in PCM materials (e.g. temperature regulation), without undue losses. In this case, one should  seek the lowest conductivity and kw =L  1, to ensure that heat losses are slow and controlled by the isolating container, as simulated in Fig. 6. The time scale then tends to t  XLqk/(kwDT). Heat transfer in the cold fluid may also play a significant role in the kinetics of solidification (i.e. heat discharge). The conditions when this occurs are simulated in Fig. 7. These calculations show that one may neglect limitations of heat transfer in the cold fluid when the heat transfer coefficient sufficiently high heat transfer coefficient, i.e. h* = Xmaxh/k P 100. Otherwise, one should consider mixed control, mainly for intermediate values of heat transfer coefficient, 0.1 < h* < 10.

A finite difference method has been developed to obtain solutions for the kinetics of solidification, combining the effects of heat conduction in the growing solidified layer, conduction through the container wall and heat transfer into the cold fluid. The risks of instability and divergence of numerical codes for this type of moving boundary problems has been eliminated by a suitable transformation of variables, which is based on normalization of the scale w = x/X(t). The accuracy of the numerical method was verified by comparing numerical and analytical solutions for the case when the process is controlled solely by conduction through the solidified layer. Stability and convergence of the numerical solutions was also demonstrated for wider ranges of conditions by varying the mesh sizes dw and relative increments of ds per time step. These solutions were then used to obtain typical criteria for mixed control

N. Vitorino et al. / International Journal of Heat and Mass Transfer 53 (2010) 5335–5342

and conditions when the dynamics of the process is controlled mainly by isolating walls or by the thermal inertia of the phase change material, due to combination of highly exothermic solidification and low thermal conductivity. Appendix A An implicit finite difference method was used to solve the partial differential equations for the heat balances, combined with the motion of the solid/liquid interface. The generic algorithms are based on implicit finite difference methods, as previously derived for diffusion controlled processes [22]. These algorithms are based on the following formulae for the time derivative, first and second spatial derivatives:



 @h hi;jþ1  hi;j ¼ @s i ds   @h hiþ1;jþ1  hi1;jþ1 ¼ @w i;jþ1 2dw ! 2 @ h hiþ1;jþ1  2hi;jþ1 þ hi1;jþ1 ¼ @w2 ðdwÞ2

ða1Þ ða2Þ

where wi+1 and wi are consecutive mesh points, dw = wi+1  wi, and si, and si+1 are consecutive time values, with ds = si+1  si, and hi,j denotes the values of the dependent function at spatial mesh point wi and time si. On combining these derivatives with Eq. (12) one obtains generic finite difference implicit formulae for heat conduction through the growing PCM layer:

þ

o b2;i hi;jþ1

þ b3;i hiþ1;jþ1 ¼

o b4;i

ða5Þ

b2;i ¼ 2  ðdw  nÞ2 =ds

o

ða6Þ

b3;i ¼ 1 þ 0:5uwj ðhn;j  hn1;j Þ

ða7Þ

o

ða8Þ

Similarly one can use Eq. (a4) as a finite difference algorithm for heat conduction through the wall, as hi1,j+1  [2 + (dz)2/ (a*ds)]hi,j+1 + hi+1,j+1 = hi,j(dz)2/(a*ds), which corresponds to: o

b1;i ¼ 1 o b2;i

ða9Þ 2



¼ 2  ðdzÞ =ða dsÞ

ða10Þ

b3;i ¼ 1

ða11Þ

o

b4;i ¼ hi;j ðdzÞ2 =ða dsÞ

ða12Þ

for i = 1, . . . , m  1. The continuity at the wall/PCM interface yields the additional finite difference algorithm (Eq. (a4)) for i = m, with: 



hm1;jþ1 k ndw=dz þ ð1 þ k ndw=dzÞhm;jþ1  hmþ1;jþ1 ¼ 0

ða13Þ

Finally, conservation at the cold fluid/wall interface yields: 



½1 þ dzðh =k ÞÞh0;jþ1  h1;jþ1 ¼ 0

or

1 hL hX þ1 þ kw k  1 kL k ðT m  T w;h Þ=ðT m  T c Þ  þ1þ kw X hX

ðT w;c  T c Þ=ðT m  T c Þ 



ða16Þ ða17Þ

This yields the corresponding relations in dimensionless terms as: 

 

h0;0 ¼ ðh þ nh kw =L þ 1Þ1 

ða18Þ ða19Þ

The remaining points for the initial temperature gradient are then estimated by assuming linear dependence on distance. One also needs a starting condition for the thickness of solidified layer (X0) and corresponding time (t0). This was based on assuming a quasi steady state approximation, i.e.:

½ðX=kÞ þ ðL=kw Þ þ ð1=hÞ@T=@t  DT=ðqkÞ

ða20Þ

or: 



ð0:5n2 þ nðL =kw Þ þ 1=h Þ=/  s

ða21Þ

References

o

b4;i ¼ hi;j ðdw  nÞ2 =ds

ða15Þ

ða4Þ

for i = m + 1, . . . , n  1, with

b1;i ¼ 1  0:5uwj ðhn;j  hn1;j Þ

      @T T w;h  T w;c @T  kw hðT w;c  T c Þ ¼ kw  kw @x @x L   T m  T w;h k X



i;jþ1

o b1;i hi1;jþ1

A quasi steady state approximation can be used to obtain starting conditions, i.e. the initial temperature profile; this is based on assuming series association of resistive contributions as follows:

hm;0 ¼ 1  ð1=ðnh Þ þ kw L =n þ 1Þ1 ða3Þ

5341

ða14Þ

Since one assumed that the PCM solid/liquid interface remains at melting temperature (hn,j+1 = 1), this combines with Eq. (a4) to obo o tain b1;n1 hn2;jþ1 þ b2;n1 hn1;jþ1 ¼ b4;i with b4;i ¼ b4;i  b3;n1 hn;jþ1 . This starts a series of backwards substitutions yielding a series of linear equations containing only two independent unknowns. The very o last in this series b1;1 h0;jþ1 þ b2;1 h1;jþ1 ¼ b4;i is then combined with Eq. (a14) to obtain the values of h0,j+1 and h1,j+1. This starts the series of forward substitution to compute the remaining unknowns.

[1] B. Zalba, J.M. Martin, L.F. Cabeza, H. Mehling, Review on thermal energy storage with thermal storage by latent heat using phase change materials, Energy Convers. Manage. 45 (2004) 263–275. [2] A. Abhat, Low temperature latent heat thermal energy storage: heat storage materials, Sol. Energy 30 (1983) 313–332. [3] G.A. Lane, Solar Heat Storage: Latent Heat Material Technology, vol. II, CRC Press, Florida, 1986. [4] I. Dincer, M.A. Rosen, Thermal Energy Storage, Systems and Applications, John Wiley & Sons, Chichester, England, 2002. [5] A.F. Regin, S.C. Solanki, J.S. Saini, Heat transfer characteristics of thermal energy storage system using PCM capsules: a review, Renew. Sust. Energy 12 (2008) 2458–2483. [6] M. Hadjieva, S. Kanev, J. Argirov, Thermophysical properties of some paraffins applicable to thermal energy storage, Sol. Energy Mater. Sol. Cells 27 (1992) 181–187. [7] J. Stefan, Sber. Akad. Wiss. Wien. 79 (1879) 161; J. Stefan, Sber. Akad. Wiss. Wien. 98 (1890) 965. [8] J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1984. [9] J. Crank, The Mathematics of Diffusion, second ed., Clarendon Press, Oxford, 1975. pp. 305–306. [10] V.R. Voller, An implicit enthalpy solution for phase change problems: with application to a binary alloy solidification, Appl. Math. Model. 11 (1987) 110– 116. [11] J.R. Frade, M. Cable, Re-examination of the basic theoretical model for the kinetics of solid state reactions, J. Am. Ceram. Soc. 385 (1992) 1949–1958. [12] N. Shamsundar, E. Rooz, Numerical methods for moving boundary problems, Handbook of Numerical Heat Transfer, Wiley-Interscience, 1988. [13] V.R. Voller, M. Cross, Accurate solutions of moving boundary problems using the enthalpy method, Int. J. Heat Mass Transfer 24 (1981) 545–556. [14] J. Crank, Diffusion with immobilizing reaction, Trans. Faraday Soc. 53 (1957) 1147. [15] J. Douglas, T.M. Gallie, On the numerical integration of a parabolic differential equation subject to a moving boundary conditions, Duke Math. J. 22 (1955) 557. [16] R.S. Gupta, D. Kumar, A modified variable time step method for the onedimensional Stefan problem, Comput. Methods Appl. Mech. Eng. 23 (1980) 101–109. [17] J. Crank, R. Gupta, Isothermal migration in two dimensions, Int. J. Heat Mass Transfer 18 (1975) 1101–1117. [18] S. Savovic, J. Caldwell, Numerical solution of Stefan problem with timedependent boundary conditions by variable space grid method, Therm. Sci. 13 (2009) 165–174.

5342

N. Vitorino et al. / International Journal of Heat and Mass Transfer 53 (2010) 5335–5342

[19] J.U. Brackbill, J.S. Salzman, Adaptive zoning for singular problems in two dimensions, J. Comput. Phys. 46 (1982) 342–368. [20] J. Caldwell, Y. Kwan, A brief review of several numerical methods for onedimensional Stefan problems, Therm. Sci. 13 (2) (2009) 61–72. [21] J.R. Frade, M. Cable, Numerical solutions for mixed control of solid state reactions, J. Am. Ceram. Soc. 78 (1995) 90. [22] S. Lin, Z. Jiang, An improved quasi-steady analysis for solving freezing problems in a plate, a cylinder and a sphere, ASME J. Heat Transfer 125 (2003) 1123–1238. [23] J.F. Wang, H.Q. Xie, et al., Thermal properties of paraffin based composites containing multi-walled carbon nanotubes, Thermochim. Acta 488 (2009) 39–42. [24] Z.G. Zhang, X.M. Fang, Study on paraffin/expanded graphite composite phase change thermal energy storage material, Energy Convers. Manage. 47 (2006) 303–310.

[25] W.L. Wang, X.X. Yang, et al., Enhanced thermal conductivity and thermal performance of form-stable composite phase change materials by using betaaluminum nitride, Appl. Energy 86 (2009) 1196–1200. [26] F.P. Incropera, D.P. DeWitt, Fundamentos de Transferência de Calor e de Massa, Rio de Janeiro, LTC, 1998. [27] P.D. Silva, L.C. Gonçalves, et al., Transient behaviour of a latent-heat thermalenergy store: numerical and experimental studies, Appl. Energy 73 (2002) 83– 98. [28] A. Sari, K. Kaygusuz, Thermal and heat transfer characteristics in a latent heat storage system using lauric acid, Energy Convers. Manage. 43 (2002) 2493– 2507. [29] A. Najjar, A. Hasan, Modeling of greenhouse with PCM energy storage, Energy Convers. Manage. 11 (2008) 3338–3342.