Contact melting of a phase change material inside a heated parallelepedic capsule

Energy Conversion & Management 42 (2001) 35±47

www.elsevier.com/locate/enconman

Contact melting of a phase change material inside a heated parallelepedic capsule Marcel Lacroix DeÂpartement de geÂnie meÂcanique, Faculte de GeÂnie, Universite de Sherbrooke, Sherbrooke, QueÂ., Canada J1K 2R1 Received 30 August 1999; accepted 14 March 2000

Abstract A mathematical model is presented for the contact melting of a subcooled phase change material (PCM) inside a heated parallelepipedic capsule. Results indicate that the melted fraction from close contact melting at the bottom of the capsule is larger, by an order of magnitude, than that from the conduction dominated melting at the top. The melting process is chie¯y governed by the magnitude of the Stefan number Ste. It is also strongly in¯uenced by the side dimensions of the capsule Lx and Lz. The total melting time increases linearly with the degree of subcooling R, but its e€ect is marginal for low and moderate subcooling (R R 1.5) and for Ste R 0.1. Finally, an analytical expression is determined for the duration of close contact melting of a parallelepipedic PCM block at the melting temperature resting on a heated plate. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Contact melting; Modelling; Parallelepipedic capsule

1. Introduction Close contact melting occurs when a heated surface and a solid phase change material (PCM) are pressed against each other while the solid is being melted. The resulting thin liquid layer that exists between the solid PCM and the heated plate is then continuously squeezed out by the force that presses them against each other. The problem of contact melting has been the subject of many investigations related to the E-mail address: [email protected] (M. Lacroix). 0196-8904/01/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 1 9 6 - 8 9 0 4 ( 0 0 ) 0 0 0 4 7 - 9

36

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

Nomenclature Ar Cf C1 F g hfs k Lx, Lz P Pr Q R S Ste Ste S1 S2 t T u,v,w U V x, y, z a d DT DT l r m n

Archimedes number, ˆ …rs ÿrf †L3x g=rs nf2 heat capacity of melt (J kgÿ1 Kÿ1) constant, Eq. (19) constant Eq. (20) acceleration of gravity (m sÿ2) latent heat of fusion (J kgÿ1) thermal conductivity (W mÿ1Kÿ1) side dimensions of capsule (m) pressure (Pa)  Prandtl number ˆ nf af ¯ow rates, Eqs. (13) and (14) (m2 sÿ1) degree of subcooling, ˆ ÿSte Ste thickness of liquid phase (m) Stefan number, ˆ Cf DT hfs  dimensionless number, ˆ Cf DT hfs thickness of liquid layer induced by contact melting (m) thickness of liquid phase induced by conduction melting at top (m) time (s) temperature (K) velocities (m sÿ1) characteristic velocity (m sÿ1) melting speed (m sÿ1) coordinates (m) thermal di€usivity of melt (m2 sÿ1) molten layer thickness (m) temperature di€erence, ˆ TH ÿ Tm temperature di€erence, ˆ Tc ÿ Tm constant, Eq. (17) density (kg mÿ3) dynamic viscosity (N s mÿ2) kinematic viscosity (m2 sÿ1)

Subscripts c cold H hot s solid m melting point f liquid Superscripts  indicates dimensionless quantity

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

37

fundamentals of heat transfer [1±4], lubrication [5] and latent heat thermal energy storage [6±8]. The interest for promoting direct contact melting in latent heat storage devices lies in the fact that the heat ¯ux across the thin melt layer separating the heated surface from the solid PCM is much higher than the heat transfer dominated by natural convection, which generally occurs in much thicker layers of molten material. As a result of the higher heat ¯uxes, the storage period of latent heat storage devices is considerably reduced. Recently, Hirata et al. [7] and Chen et al. [8] have examined the problem of close contact melting in rectangular capsules. Although these studies provide valuable guidance for engineering designs, they are nevertheless limited to two-dimensional capsules. The present paper overcomes this limitation by looking at the problem of melting inside parallelepipedic capsules. Moreover, it examines the e€ects of solid subcooling on the melting process. Both contact melting from the bottom of the capsule and conduction dominated melting from the top are considered simultaneously. A mathematical model is ®rst derived and next employed to predict the melting rates in terms of the key parameters.

2. Analysis A schematic representation of the capsule is depicted in Fig. 1. The solid PCM of initial height H, and width 2Lx and 2Lz is isothermal at the temperature Tc below the melting point

Fig. 1. Schematic of parallelepipedic capsule.

38

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

Tm. At the time t = 0, a constant temperature TH ˆ Tm ‡ DT is imposed at both the top and bottom surfaces of the cavity, while the lateral vertical walls remain adiabatic. Melting is eventually triggered at the top and bottom surfaces of the PCM, and as it proceeds, the solid is assumed to descend vertically while squeezing the melt out of the thin gap separating it from the bottom heated wall. In the present paper, ``liquid layer'' refers to the thin ®lm that ®lls the contact gap and ``liquid phase'' denotes the liquid above the solid PCM.

3. Contact melting at the bottom The focus is ®rst on the thin liquid ®lm. The vertical dimension of the thin gap d is assumed to be very small and is, in general, a function of both x and z. Neglecting the e€ect of thermal di€usion in the x and z directions, we recognize the energy conservation equation: u

@T @T @T @ 2T m ‡v ‡w ˆ af 2 ‡ f F r f Cf @x @y @z @y

…1†

in which the last term represents the e€ect of internal heat generation by ¯uid friction. On the left side of Eq. (1), the convection terms, whose representative scale is UDT=L, are collected. The scale of the vertical di€usion e€ect (®rst term on rightside of the equation) is af DT=d 2 , while the scale of the volumetric heating e€ect is …m rC†f …U d† 2 : A comparison of these scales reveals  that the e€ects of  convection and frictional heating can be neglected if, simultaneously, 2 2 Ud af L  1 and mf U kf DT  1: If these inequalities are true with respect to their order of magnitude, then the only term remaining in Eq. (1) is the thermal di€usion term, and the temperature distribution in the liquid layer becomes   y DT …2† Tf …y † ˆ Tm ‡ 1 ÿ d The unknown gap thickness d follows from the melting condition that controls the position of the upper surface of the liquid gap     @Tf @Ts ÿkf ˆ ÿks …3† ‡ hfs rs V1 @y yˆd @y Eq. (3) stipulates that some of the heat transferred by conduction through the liquid layer di€uses into the solid subcooled core (®rst term on right side of Eq. (3)), and only part of it is used to actually melt the PCM (second term on right side of Eq. (3)). Combining Eqs. (2) and (3) yields   DT @Ts ˆ ÿks ‡hfs rs V1 …4† kf d @y B  The gap thickness is still unknown, because the temperature gradient …@Ts @y†B at the bottom of the  solid core and the melting speed V1 are not speci®ed. The temperature gradient …@Ts @y†B will be evaluated later from the temperature distribution inside the solid core. On

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

39

the other hand, the speed V1 is related to the normal force applied by the weight of the sinking block on the ®lm layer: … Lx … Lz P…x, z † dz dx ˆ …rs ÿ rf †Lx Lz …H ÿ S †g …5† 0

0

The pressure ®eld P(x, z ) is determined from the longitudinal and sideways momentum equations, i.e. ! @u @u @u ÿ1 @P mf @ 2 u @ 2 u @ 2 u ˆ ‡ …6† ‡ ‡ u ‡v ‡w @x @y @z rf @x rf @x 2 @y 2 @z 2 @w @w @w ÿ1 @P mf @ 2 w @ 2 w @ 2 w u ‡v ‡w ˆ ‡ ‡ ‡ 2 @x @y @z rf @z @y 2 @z rf @x 2

! …7†

The scale of the convection terms (left side of Eqs. (6) and (7)) is U 2/L,while that of the di€usion terms (second term on right side of Eqs. (6) and (7)) is …mf rf †…U d 2 †: The e€ect of  convection can be neglected if …d L† 2 …ULrf mf †  1: If this inequality holds true, Eqs. (6) and (7) become @ 2 u…x, y, z † @P ˆ mf @y 2 @x

…8†

@ 2 w…x, y, z † @P ˆ mf @y 2 @z

…9†

and

Integrated twice in the vertical direction and subjected to the velocity boundary conditions u=w = 0 at y = 0 and d, Eqs. (8) and (9) read u…x, y, z † ˆ

1 @P y…y ÿ d† 2m @x

…10†

w…x, y, z † ˆ

1 @P y…y ÿ d† 2m @z

…11†

De®ning the height integrated ¯ow rates …d …d Qx …x, z † ˆ u…x, y, z † dy and Qz …x, z † ˆ w…x, y, z † dy 0

0

…12†

we get from Eqs. (10) and (11) Qx ˆ

ÿ1 @P 3 d 12mf @x

…13†

40

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

and Qz ˆ

ÿ1 @P 3 d 12mf @z

…14†

Now, integration of the mass conservation equation in the liquid layer yields …d @u @v @w dy ‡ dy ‡ dy ˆ 0 @y @z 0 @x or d dx

…d

d u dy ‡ njyˆd ÿ vjyˆ0 ‡ dz 0

…d 0

w dy ˆ 0

which becomes @Qx …x, z † @Qz …x, z † ÿ V1 ÿ 0 ‡ ˆ0 @x @z

…15†

Di€erentiating Eqs. (13) and (14) and combining them with Eq. (15) for the purpose of eliminating Qx and Qz results in a single equation for P(x, z ): @ 2P @ 2P ÿ12mf V1 ‡ 2 ˆ @x 2 @z d3

…16†

  subjected to the boundary conditions P…Lx , z† ˆ 0, @P…0; z† @x ˆ 0, P…x, Lz † ˆ 0, @P…x; 0† @z ˆ0: Eq. (16) is solved using the method of separation of variables. Its solution is 1
… ÿ 1 †n‡1 24mf V1 X 6mf V1 ÿ 2 2 L …17† cos…l x† cosh…l z† ‡ ÿ x P…x, y † ˆ n n x Lx d3 nˆ0 l3n cosh…ln Lz † d3  where ln ˆ …2n ‡ 1†p 2Lx : By substituting Eq. (17) into Eq. (5) and integrating over the contact area, the pressure distribution yields the following dimensionless expression for V1:

V1 ˆ

r Ard 4C1

3

…H  ÿ S  †

…18†

where V1 ˆ V1 Lx =nf , r ˆ rs =rf and Ar ˆ …rs ÿ rf †L3x g=rs nf2 : d, H and S have been nondimensionalized with respect to Lx : C1 is a dimensionless coecient that depends only on Lz : h  i  1 … † 2n ‡ 1 p 2 Lz X 192 C1 ˆ 5  tanh …19† p Lz nˆ0 …2n ‡ 1 †5 This function is plotted in Fig. 2. Of particular interest are C1 11 when Lz  1 and C1 1L2 z when Lz  1:

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

41

4. Melting from the top The focus is now shifted to the liquid phase above the solid PCM. Since the liquid at the top of the cavity is nearly stagnant, i.e. layers of hot and light ¯uid rest on layers of cold and dense ¯uid, heat conduction prevails and the heat balance equation reads   kf DT @Ts ˆ ks ‡hfs rs V2 …20† S @y T Eq. (20) indicates that some of the heat transferred through the liquid phase at the top di€uses into the solid core (®rst term on right side of Eq. (20) and the rest is used to melt the PCM (second term on right side of Eq. (20)).

5. Temperature distribution inside the solid core In order to evaluate V1 and V2 and, therefore, S1 and S2 from Eqs. (4), (18) and (20), the temperature gradients …@Ts @y† at the bottom and at the top of the solid subcooled core must be determined. Neglecting the e€ect of thermal di€usion in the x and z directions, we recognize the conduction equation @Ts @ 2 Ts ˆ as @t @y 2

…21†

subjected to the boundary conditions, Ts …0, t† ˆ Tm and Ts …H ÿ S, t† ˆ Tm where the origin of the y-axis has been set, for the present equation, at the bottom surface of the solid core. The

Fig. 2. C1 versus Lz Eq. (19).

42

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

initial condition is Ts …y, 0† ˆ Tc : Once again, the method of separation of variables is invoked and the solution of Eq. (21) reads   1 X … ÿ 1 †n ÿ1 ÿas g 2 t n sin g y …22† e Ts …y, t † ˆ Tm ÿ 2DT n np nˆ0  where DT ˆ Tc ÿ Tm and gn ˆ np …H ÿ S†: Di€erentiating Eq. (22) with respect to y, the temperature gradients at the boundaries of the solid core become   1   @Ts 2DT X 2 … ÿ 1 †n ÿ1 e ÿas gn t …23† ˆ2 …H ÿ S † nˆ0 @y The minus sign is for the gradient at the bottom of the solid core ( y = 0), and the plus sign is for the top ( y=H ÿ S ).

6. Total melting The simultaneous e€ect of both contact melting at the bottom of the capsule and conduction dominated melting at the top may now be obtained by substituting Eq. (23) into Eqs. (4) and (20). Combining the dimensionless form of these equations with Eq. (18) for the purpose of eliminating the gap thickness d results in the following set of coupled nonlinear di€erential equations for S1 and S2: "   #1=3   4C1 dS1 dS1 1=3   r Pr  ‡ G ˆ Ste…H  ÿ S  † …24† r Ar dt dt r Pr

dS2 Ste ˆ  ÿG  dt S

…25†

where

and

S  ˆ S1 ‡ S2

…26†

" # 1   2 2   X  2k Ste ÿn p a t n … ÿ 1 † ÿ1 exp G ˆ 2 …H  ÿ S  † nˆ0 Pr…H  ÿ S  †

…27†

k ˆ ks =kf and a ˆ as =af : Ste and Pr are the Stefan and Prandtl numbers, respectively. Ste is a dimensionless number used to de®ne the degree of subcooling R ˆ ÿSte=Ste: Eqs. (24) and (25) are subjected to the initial conditions S1=S2=0 at time t=0. Of particular interest is the case of no subcooling, i.e., Tc ˆ Tm : In this case, Ste ˆ 0 and Eqs. (24) and (25) may be combined to yield a single di€erential equation for S:

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

"   #1=4 dS  Arr Ste 3 Ste 1 …H  ÿ S  †1=4 ‡ ˆ   Prr S  dt 4C1 Prr

43

…28†

subjected to the initial condition S=0 at time t=0. Eqs. (24) and (25) are solved for S1 and S2 using a fourth-order Runge±Kutta method. At a given time step, both equations are solved iteratively and convergence is declared when     dS  dS1 1 Re1 ÿ  dt i‡1 dt i     dS  dS2 2 Re2 ÿ dt i‡1 dt i where i stands for the iteration number and e1 ˆ e2 ˆ 10 ÿ8 : For all cases studied in the present paper, convergence is reached, at a given time step, in less than 20 iterations. Simulations were performed on a PC 200 MHz, and the CPU time never exceeded 3 s with a time step of Dt ˆ 0:005: 7. Results and discussion Simulations were conducted for the melting of n-octadecane con®ned to an enclosure for which the height H and the width Lx were ®xed at 0.03 and 0.075 m, respectively. The width Lz ranged from 0.0225 to 0.225 m (0.3 R Lz R 3.0), and the Stefan number, Ste, varied from 9.26  10ÿ3 to 9.26  10ÿ2 (1 K R DT R 10 K). The Archimedes number, Ar, was set equal to  K). 7.73  106. The degree of subcooling R varied from 5 to 0 (ÿ50 K R DTR0  Fig. 3 shows the temporal variation of S for two di€erent Stefan numbers Ste 0.0926 and 0.00926, and three di€erent widths Lz . The degree of subcooling R is set equal to zero. As expected, the melting process is chie¯y governed by the magnitude of the Stefan number, i.e., the temperature di€erence DT: The high melting rates (the slopes dS/dt) prevailing at early

Fig. 3. Thickness of liquid phase versus time.

44

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

times decrease monotonically and reach a nearly constant value after a while. These rates are higher for larger values of Ste. S is also strongly in¯uenced by the width Lz . As Lz is augmented, so is the path in the z direction through which the liquid between the bottom heated surface and the solid PCM is squeezed out, and as a result, the contact melting process is slowed. Fig. 4 depicts the timewise variation of S1, S2 and S as given by Eqs. (24)±(26), respectively, for the case Ste = 9. 26  10ÿ2 and Lz =1. This example clearly indicates that the melting fraction of contact melting S1 is larger, by an order of magnitude, than that from the conduction dominated melting S2 at the top, and the contact melting is dominant during the entire melting process. The temporal variation of the gap thickness d corresponding to the test cases of Fig. 3 is elucidated in Fig. 5. The increase in d with time is initially slow and then rapid in the ®nal stages of the melting process, i.e., when H ÿ S tends to zero. The e€ect of the degree of subcooling on the melting process is illustrated in Fig. 6. For the three cases shown in this ®gure, Ste was set equal to 9.26  10ÿ2 and Lz was ®xed at 1.0. As expected, the melting rate decreases as the degree of subcooling increases. Fig. 7 shows that the ratio of the total melting time of a subcooled PCM block over that of a block at the melting point (RMT) increases linearly with the degree of subcooling (R ). The slope of the RMT versus R relation gets steeper with the Stefan number but it was determined that it remains independent of the width Lz . For Ste = 9.26  10ÿ3, the slope is equal to 0.00315, and for Ste = 9.26  10ÿ2, it is equal to 0.028. Moreover, Fig. 7 reveals that for low and moderate subcooling, i.e., R R 1.5 and for Ste R 0.1, the increase in the melting time remains marginal. Finally, an interesting feature of the foregoing model is that it can readily be applied to the problem of contact melting of a parallelepipedic block at the melting temperature resting on an isothermally heated plate by setting the second term on the right side of Eq. (28) equal to zero and by rede®ning the dimensionless number Ar as Ar=gL3x/n2f . In this case, the exact solution of Eq. (28) is given by  4=3 3    3=4 H ÿS ˆ H ÿ Ft …29† 4

Fig. 4. Temporal variation of thicknesses S1, S2 and S .

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

45

Fig. 5. Temporal variation of gap thickness.

  where F ˆ ‰…Arr 4C1 †…Ste Prr †3 Š1=4 : The time required for completion of the contact melting process of a pure substance tt may then be evaluated by setting H  ÿ S  ˆ 0 in Eq. (29). This yields tt ˆ

4 H 3=4 3 F

…30†

As an example, for a rectangular block of n-octadecane (Lx=Lz), the dimensional form of Eq. (30) becomes 3=4  Ste tt ˆ 8160L1=2 …31† x H where tt is expressed in seconds, and H and Lx are in meters. A plot of this relation is provided in Fig. 8.

Fig. 6. E€ect of the degree of subcooling.

46

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

Fig. 7. Ratio of the melting times versus the degree of subcooling.

Fig. 8. Time required to melt a parallelepipedic block resting on a heated plate.

8. Conclusion A mathematical model was presented for the contact melting of a subcooled PCM inside a heated parallelepipedic capsule. Results indicate that the melted fraction from close contact melting at the bottom of the capsule is larger, by an order of magnitude, than that from the conduction dominated melting at the top. The melting process is chie¯y governed by the magnitude of the Stefan number Ste. It is also strongly in¯uenced by the side dimensions of the capsule Lx and Lz. The total melting time increases linearly with the degree of subcooling R, but its e€ect is marginal for low and moderate subcooling (R R 1.5) and for Ste R 0.1. Finally, an analytical expression was determined for the duration of close contact melting of a parallelepipedic PCM block at the melting temperature resting on a heated plate. Acknowledgements The author is grateful to the Natural Sciences and Engineering Research Council of Canada for their ®nancial support.

M. Lacroix / Energy Conversion & Management 42 (2001) 35±47

47

References [1] Saito A, Hong H, Hirokane O. Heat transfer enhancement in the direct contact melting process. Int J Heat Mass Transfer 1992;35(2):295±305. [2] Hong H, Saito A. Numerical method for direct contact melting in transient process. Int J Heat Mass Transfer 1993;36(8):2093±103. [3] Moallemi M, Webb B, Viskanta R. An experimental and analytical study of close-contact melting. J Heat Transfer 1986;108:894±9. [4] Bareiss M, Beer H. An analytical solution of the heat transfer process during melting of an un®xed solid phase change material inside a horizontal tube. Int J Heat Mass Transfer 1984;27(5):739±46. [5] Bejan A. Contact melting heat transfer and lubrication. In: Advances in Heat Transfer, vol. 24. New York: Academic Press, 1994. [6] Roy K, Sengupta S. The melting process within spherical enclosures. J Heat Transfer 1987;109:460±2. [7] Hirata T, Makino Y. Analysis of close-contact melting for octadecane and ice inside isothermally heated rectangular capsule. Int J Heat Mass Transfer 1991;34(12):3097±106. [8] Chen W, Cheng S, Luo Z, Gu W. Analysis of contact melting of phase change materials inside a heated rectangular capsules. Int J Research 1995;19:337±45.