- Email: [email protected]
Numerical investigation of the phase change process of low melting point metal
International Journal of Heat and Mass Transfer 100 (2016) 899–907
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Numerical investigation of the phase change process of low melting point metal Xiao-Hu Yang a, Si-Cong Tan a, Jing Liu a,b,⇑ a b
Beijing Key Lab of Cryo-Biomedical Engineering and Key Lab of Cryogenics, Technical Institute of Physics and Chemistry, Chinese Academy of Sciences, Beijing 100190, China Department of Biomedical Engineering, School of Medicine, Tsinghua University, Beijing 100084, China
a r t i c l e
i n f o
Article history: Received 1 April 2016 Received in revised form 28 April 2016 Accepted 30 April 2016 Available online 23 May 2016 Keywords: Phase change material Low melting point metal Parametric numerical investigation Dimensionless correlation
a b s t r a c t Low melting point metals (LMPMs) are a new kind of phase change materials (PCMs) which exhibit excellent heat extraction capability due to their high thermal conductivity. Such a characteristic renders it distinctive melting mode which hardly exists in conventional paraffin PCMs. In this paper, comprehensive numerical investigation is carried out to reveal the melting process of LMPMs. First, comparative simulations of typical paraffin eicosane and typical LMPM gallium is conducted. Then, dimensionless analysis is performed based on a series of parametric studies and corresponding dimensionless correlations are given under both constant wall temperature condition and constant heat flux condition. Further, the critical dimensionless values within which the correlations are applicable are suggested. Finally, a preliminary exploration is performed to illustrate the influence of natural convection under relatively high Ra number condition. This study can serve as a valuable reference for future practical thermal design of LMPMs based thermal management systems and thermal energy storage systems. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction Phase change materials (PCMs) are materials which undergo melting/solidification process at a constant or nearly constant temperature and absorb/release thermal energy from/to the surroundings. PCMs are widely used in the field of thermal management and thermal energy storage [1], such as electronics cooling [2–4], energy storage for buildings [5], solar energy systems [6], thermal comfortable textile design [7] and space systems [8]. Conventionally, PCMs are classified into organic PCMs and inorganic PCMs. Organic PCMs are generally used in relatively low temperature situations (from 35 °C to about 70 °C). Paraffin is a typical kind of organic PCM and has been extensively researched and used. The main drawback of paraffin PCM lies in its low thermal conductivity, which seriously hinders the heat conduction inside the PCM and thus decreases its efficiency. Many efforts have been made to improve such a situation, such as internal fin [3,9,10], metallic foam [11], and nano-particle inclusion [12,13]. However, such improvement is still limited. Ge and Liu [14] had ever proposed that using low melting point metals (LMPMs), typically gallium and its alloys and bismuth based ⇑ Corresponding author at: Technical Institute of Physics and Chemistry, Chinese Academy of Sciences, Beijing 100190, China. Tel.: +86 10 82543765; fax: +86 10 82543767. E-mail address: [email protected] (J. Liu). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.04.109 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
and indium based alloys, as a new kind of PCMs and successfully used it in the thermal management of smart phone [15,16] which exhibited good performance. The main advantages of LMPMs lie in their high thermal conductivity and high volumetric latent heat. Since the thermophysical properties of LMPMs are much different from that of conventional paraffin PCMs, some distinguishing features exist during their phase change process. However, up till now, there are rather limited literatures to reveal such a difference. This paper is dedicated to provide a comprehensive understanding about the phase change characteristics of LMPMs via parametric numerical investigation. The phase change process of PCMs includes melting and solidification. In this paper, attentions will be focused on the former one. First, comparative investigation of the melting process of typical paraffin eicosane and typical LMPM gallium is conducted, which reveals the superiority of using LMPMs as PCMs and their special characteristics. Then, dimensionless analysis of the melting process of LMPM is carried out based on a series of parametric study and corresponding dimensionless correlations are given under both constant wall temperature condition and constant heat flux condition. Further, the critical dimensionless values within which the correlations are applicable are suggested. Finally, preliminary exploration is performed to reveal the melting process of LMPM under relatively high Ra
900
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
Nomenclature A C cp Fo g H DH h k Nu p q00 Ra Ste T t u v W x, y
Greek letters thermal diffusivity (m2/s) b volume thermal expansivity (1/K) c local liquid fraction e small constant g dimensionless length along y direction h dimensionless temperature l dynamic viscosity (kg/s/m) n dimensionless length along x direction q mass density (kg/m3) u melting fraction s dimensionless time
source item coefficient in momentum equation constant coefficient of A specific heat capacity (J/kg/K) Fourier number gravitational acceleration (m/s2) height (m) fusion latent heat (J/kg) enthalpy (J/kg) thermal conductivity (W/m/K) Nusselt number pressure (N/m2) heat flux (W/m2) Rayleigh number Stefan number temperature (K) time (s) velocity along x direction (m/s) velocity along y direction (m/s) width (m) rectangular coordinates (m)
a
Subscripts c cold h hot l liquid phase m melting point ref reference s solid phase w wall
number condition where the influence of natural convection becomes significant.
avoiding zero denominator. In the current simulation, C = 105 and
e ¼ 0:001.
Energy equation: 2. Numerical method 2.1. Enthalpy-porosity method Enthalpy-porosity method, which was proposed by Voller et al. [17], is widely used for the numerical simulation of solid–liquid phase change problems [18–22]. The main features of this method are that: (1) The energy equation is based on the variable enthalpy but not temperature, which can conveniently take the phase change process into consideration; (2) the phase change region (mushy region) is treated as a porous domain, and the local liquid fraction of the PCM is regarded equal to the porosity. In enthalpyporosity method, the governing equations in two-dimensional Cartesian coordinate are given by: Continuity equation:
@ q @ðquÞ @ðqv Þ þ ¼0 þ @x @y @t
ð1Þ
Momentum equations with Boussinesq approximation invoked:
@ðquÞ @ðquuÞ @ðqv uÞ @p @ þ þ ¼ þ @t @x @y @x @x
l
@u @ @u þ l @x @y @y
þ Au @ðqv Þ @ðquv Þ @ðqvv Þ @p @ þ þ ¼ þ @t @x @y @x @x
l
@v @x
þ
þ Av þ qgbðT T m Þ
@ @y
l
@v @y
ð4Þ
8RT cp dT; T < T solidus > T > > < R ref T solidus cp dT þ cDH; T solidus 6 T 6 T liquidus h¼ T ref > > > : R T solidus c dT þ DH þ R T c dT; T > T liquidus p T liquidus p T ref
ð5Þ
where k and cp are the thermal conductivity and specific heat capacity of the PCM, respectively. h denotes the enthalpy of the PCM and is defined by Eq. (5), where, Tref is a reference temperature and the enthalpy at this temperature is set to be zero, DH is the fusion latent heat, Tsolidus and Tliquidus are the temperature limit at which phase change starts and ends up, respectively. Eq. (5) indicates that, in the mushy region, if Tsolidus = Tliquidus = Tm, then T = Tm; if T solidus – T liquidus , T can be calculated by the following definition
c¼
T T solidus T liquidus T solidus
ð6Þ
RT where c can be obtained from Eq. (5) as c ¼ h T refsolidus cp dT =DH. It
ð2Þ
@ðqhÞ @ðquhÞ @ðqv hÞ @ k @h @ k @h þ þ þ ¼ @t @x @y @x cp @x @y cp @y
ð3Þ
where g is the gravitational acceleration; q, l, b and Tm are the density, dynamic viscosity, volume thermal expansivity and melting point of the PCM; Au and Av are the momentum dissipation source items, which are used for suppressing the velocity in the mushy and solid region, hence it is a function of local liquid fraction c. Brent et al. [23] defined such a function in the form of A ¼ Cð1 cÞ2 =ðc3 þ eÞ, where C is a large constant number to enhance the suppression effect and e is a small constant set for
is worthy to note here that the term local liquid fraction c should be distinguished from the term melting fraction u, the latter of which means the ratio of the volume of the molten PCM to the volume of the initial solid PCM and this definition will be widely used in the following. Mathematically, u is the volume weighted average of c R over the whole cavity defined as u ¼ V cdV=V. 2.2. Validation of the numerical method Before parametric study, the numerical method is validated first by simulating the melting process of pure gallium in a rectangular cavity (47.6 mm 47.6 mm) which had ever been experimentally investigated by Beckermann et al. [24]. The cavity is heated at its left wall with overheating of 10.2 °C and cooled at the right wall with subcooling of 4.8 °C. Here, overheating and subcooling are relative
901
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907 Table 1 Thermophysical properties of gallium. Parameter
Value
Parameter
Value
Melting point Tm (K) Fusion latent heat DH (J/kg) Density (solid) qs (kg/m3) Density (liquid) ql (kg/m3) Specific heat capacity cps (solid) (J/kg/K)
302.93 80,160 5903.7 6094.7 372.3
Specific heat capacity (liquid) cpl (J/kg/K) Thermal conductivity (solid) ks (W/m/K) Thermal conductivity (liquid) kl (W/m/K) Expansion coefficient b (1/K) Viscosity l (kg/m/s)
397.6 33.49 33.6767 1.2 104 1.75 103
(a)
within which the natural convection is laminar [27]. As can be seen later that all of the cases investigated in this paper are laminar and hence laminar viscous flow mode is adopted in the numerical method for all of the simulations. SIMPLE algorithm is selected for velocity–pressure coupling and body force weighted scheme is used for the spatial discretization of pressure. Momentum and energy equations use second order upwind scheme to ensure accuracy. Convergence criteria of 104 is set for the momentum equation and 107 for the energy equation. Time step of 0.01 s is chosen for the transient computation with compromise of enough accuracy and time consumption. Grid independence study is also conducted with 6 sets of grid sizes (20 20, 30 30, 40 40, 50 50, 60 60, 80 80), the result of which indicate that 60 60 grid is accurate enough for this simulation since the difference of the calculated melting fraction and heat flux based on 60 60 grid and 80 80 grid are less than 0.5%. Typical computational results are presented in Fig. 1 in comparison with the experimental data in Ref. [24] in the form of dimensionless parameters. Fig. 1(a) shows the temperature distribution on three specific horizontal lines in the cavity and Fig. 1(b) exhibits the solid–liquid interface under 5 different boundary conditions. Herein, hl and hs are dimensionless temperature in the liquid and solid region of the PCM defined as hl ¼ ðT T m Þ=ðT h T m Þ and hs ¼ ðT T m Þ=ðT m T c Þ, where Th and Tc are the temperature of the left hot wall and right cold wall of the cavity, respectively; n and g are dimensionless length in horizontal and vertical direction defined as n ¼ x=W and g ¼ y=H, where W and H are the width and height of the cavity. It can be seen from Fig. 1(a) and (b) that good agreement between the current computational results and the experimental data is achieved, which implies the reliability of the current numerical method.
3. Results and discussions 3.1. Comparison with paraffin
(b) Fig. 1. Comparison of the current computational results and the experimental results in literature [24]. (a) Temperature profiles on 3 specified lines. (b) Liquid– solid interfaces under 5 different boundary conditions. Case 1: DT h ¼ T h T m ¼ 10:2 C, DT c ¼ T m T c ¼ 4:8 C; case 2: DT h ¼ 10:2 C, DT c ¼ 9:8 C; case 3: DT h ¼ 10:2 C, DT c ¼ 19:8 C; case 4: DT h ¼ 5:2 C, DT c ¼ 14:8 C; case 5: DT h ¼ 15:2 C, DT c ¼ 14:8 C.
to the melting point of gallium 29.78 °C (302.93 K). The top and bottom wall of the cavity are set to be adiabatic. Thermophysical properties of gallium refer to Ref. [25] and are listed in Table 1. Many literatures [21,26] have verified the accuracy of using two-dimensional simulation to compute the realistic threedimensional problem and in this paper, 2-D simulation is performed. Commercial codes Fluent 6.3 is used to carry out the simulation, which is running based on the aforementioned enthalpy-porosity method. Boussinesq approximation is applied to activate natural convection. The Ra number (Eq. (9)) of this simulation case is 3.24 105, which is much smaller than the critical value 109
Based on the aforementioned numerical method, a series of parametric study are performed. At first, a comparative investigation is conducted to illustrate the different melting process of LMPMs and conventional paraffin PCMs. Typical LMPM pure gallium and typical paraffin eicosane are chosen here for comparative simulation and the main thermophysical properties of which are listed in Table 2. The melting process happened in a rectangular cavity with dimension of 5 mm in width and 20 mm in height. The top, right and bottom wall of the cavity are set to be adiabatic. Constant temperature condition is set for the left wall of the cavity with 3 different over-temperature sets (i.e. 5 K, 15 K and 25 K). Here, over-temperature indicates the temperature relative to the melting point of the PCM. Transient simulation is carried out and the initial temperature of the whole PCM is set to be at its melting point. To intuitively illustrate the different melting modes of eicosane and gallium, liquid fraction contours and over-temperature contours at different melting fractions (i.e. u ¼ 0:3; 0:6 and 0:9) for the two PCMs are presented, in Fig. 2(a) for eicosane and
902
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
Table 2 Thermophysical properties of gallium and eicosane.
a
PCM
qa (kg/m3)
cp (J/kg/K)
l (103 kg/m/s)
k (W/m/K)
b (1/K)
T m (K)
DH (J/kg)
Gallium Eicosane [19]
6094.7 785
397.6 2460
1.75 exp(1790/T 4.25) [21]
33.6767 0.15
0.00012 0.001
302.93 309.15
80,160 247,000
The density of liquid phase PCM at its melting point.
Fig. 2. Liquid fraction and temperature distribution in eicosane and gallium under different melting fractions. Herein, over-temperature indicates temperature relative to the melting point of the PCM. (a) Eicosane; (b) gallium.
Fig. 2(b) for gallium. The streamline and flow direction of the natural convective flow of liquid PCM at melting fraction of 0.6 and 0.9 are also inserted in the liquid fraction contours. It can be seen that for eicosane, the liquid–solid interface is inclined from bottom-left to top-right, especially under large melting fraction conditions, which is caused by the clockwise natural convective flow and thus induced heat transfer enhancement in the upper cavity region. While in the lower region, heat transfer is mainly dependent on heat conduction which is relatively weak thus the PCM there melts more slowly. The temperature distribution in the liquid region of eicosane exhibits horizontal slicing which is a typical feature of natural convection. While for gallium the interface is nearly vertical all the time, as well as its temperature distribution, which is a typical feature of heat conduction. These phenomena qualitatively indicate that for eicosane, natural convection plays an important role in the total heat transfer, while for gallium heat conduction is dominant. Fig. 2 also presents the velocity distribution in the liquid PCM region. It can be found that in fact, the natural convective flow in liquid gallium is more intensive since the velocity of which is several times larger than that in eicosane. However, such an intensive flow does not make noticeable contribution to the total heat transfer. The reason causing the different melting modes between gallium and eicosane mainly lies in their large difference in thermal conductivity. As shown in Table 2, the thermal conductivity of gallium is more than 200 times larger than that of eicosane. Hence, in eicosane, the heat conduction is very weak, and natural convection will significantly enhance the heat transfer. While in gallium, heat
conduction is so strong that the natural convective heat transfer becomes relatively weak in this case though the convection is more intensive. Fig. 3 quantitatively presents the melting fraction curves and heat flux curves over time of eicosane (Fig. 3(a) and (b)) and gallium (Fig. 3(c) and (d)) at different over-temperatures. It can be seen from Fig. 3(a)–(d) that, for all cases, a very high heat flux exists at the beginning of the melting due to the large temperature difference between the cavity wall and adjacent PCM. Quickly, a liquid film is formed between the wall and the solid PCM. Such a liquid region is in fact a thermal resistance for the heat transfer from the wall to the solid PCM. With time going on, the liquid region becomes thicker and thicker and corresponding thermal resistance keeps getting larger, which leads to the gradual decrease of the heat flux as well as the growth rate of the melting fraction. At the same time, natural convection in the liquid region is developed gradually due to the temperature difference and thus induced density difference of the liquid PCM, which promotes the heat transfer in some degree. However, at the end of the melting process, there generally exists a sharp drop of the heat flux, as can be seen in Fig. 3(b) and (d). Such a drop is inferred to also be caused by the natural convection. The reason lies in that at the end, only small amount of solid PCM exist at the bottom-right corner of the cavity, while the liquid flows clockwise at the bottom, namely from right to left (Fig. 2), which is opposite to the heat flux direction and thus the heat transfer is adversely hindered at this point. The larger the over-temperature is, the more sharply the heat flux decreases.
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
903
Fig. 3. Melting fraction and heat flux versus time of different PCMs under different overheating temperatures. Dimension of the cavity is 5 mm 20 mm (width height). (a) Eicosane, melting fraction; (b) eicosane, heat flux; (c) gallium, melting fraction; (d) gallium, heat flux; (e) melting fraction versus FoSte; (f) Nu versus FoSte.
It can also be seen that the time needed to completely melt eicosane are 515 s, 910 s and 3080 s respectively under overtemperatures of 25 K, 15 K and 5 K. While for gallium, the time are only 8.5 s, 13.5 s and 38 s, respectively. The heat flux in the cases with gallium as PCM is nearly two orders of magnitude larger than that with eicosane as PCM at the same over-temperature. Hence it can be concluded that gallium has much superior heat extraction efficiency as a PCM than eicosane. Dimensionless analysis for the melting process is conducted for both eicosane and gallium, as shown in Fig. 3(e) and (f). Here, Fourier number is dimensionless time, defined as Fo ¼ at=W 2 , where a is the thermal diffusivity of the PCM and W denotes the width of the cavity which is regarded as the characteristic length here. Stefan number is generally used for phase change problem in the form of Ste ¼ cp DT=DH, where DT is the characteristic temperature difference and defined here as DT ¼ T w T m , T w is the left
wall temperature of the cavity. To take into account both the transient heat conduction and phase change, the product of Fo number and Ste number can serve as a dimensionless time [21,24], namely
s ¼ FoSte ¼
acp DTt DHW 2
ð7Þ
Nu number is used to describe the dimensionless heat transfer intensity, defined as
Nu ¼
W q00 k DT
ð8Þ
where q00 denotes the heat flux on the left wall. Fig. 3(e) shows the melting fraction curve versus FoSte and Fig. 3(f) is Nu versus FoSte. It can be seen that both in Fig. 3(e) and (f) the dimensionless curves of gallium under 3 different overtemperatures are almost coincident. While for eicosane, significant
904
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
Table 3 Sets of cavity sizes and boundary conditions and corresponding Ra numbers. Case
Case 1 Case 2 Cases 3/4/5 Case 6 Case 7 Case 8
Width (mm)
3 5
8 10
Height (mm)
10 10 20 40 30 20
Constant wall temperature
Constant heat flux
Over-temperature (K)
Ra
Heat flux (kW/m2)
Ra
15 15 5/15/25 15 15 15
119 552 184/552/921 552 2263 4420
50 50 30/50/100 50 50 50
35 273 164/273/547 273 1792 4376
Fig. 4. Melting process of gallium in different cavities under different wall temperatures. W and H mean the width (mm) and height (mm) of the cavity and the temperature means over-temperature of the wall relative to the melting point of the PCM, respectively. (a) Melting fraction versus time; (b) melting fraction versus FoSte; (c) heat flux versus time; (d) Nu versus FoSte.
difference exists. The higher the wall temperature is, the shorter the dimensionless time for complete fusion is and the higher Nu number is. Such a difference is essentially a quantitative reflection of the aforementioned influence of natural convection. Rayleigh number is generally used in natural convection problems which can be viewed as the ratio of buoyancy and viscosity forces multiplied by the ratio of momentum and thermal diffusivities, defined as
Ra ¼
gbDTcp q2 W 3 lk
ð9Þ
In fact, Ra number can also be understood as the ratio between the rate of convective heat transfer to the rate of heat conduction, up to a numerical factor. Under the current cavity geometry and an overtemperature of 15 K, the Ra number for eicosane and gallium are 45,659 and 552 respectively, which quantitatively indicates that significant natural convective heat transfer relative to heat conduction exists in eicosane while it is rather weak in gallium. Natural convection enhances the heat transfer and thus increases Nu
number and decreases the dimensionless time for complete fusion. It is worthy to note that, although gallium exhibits inferior heat transfer capacity than eicosane under the dimensionless scale (Fig. 3(e) and (f)), it is much superior under real spatiotemporal scale (Fig. 3(a) to (d)). Up to this point, the reason why the dimensionless curves of gallium under different conditions coincide is very clear: the natural convective heat transfer in gallium is relatively weak compared with the heat conduction and it can be neglected, thus the dimensionless analysis has no relationship with Ra number. In fact for gallium, there exist tremendous situations where the influence of natural convection can be neglected due to its high thermal conductivity which hardly exists in conventional paraffin PCMs. In the following, attentions will be focused on these situations. 3.2. Dimensionless correlations: constant wall temperature Parametric study is performed to investigate the melting process of gallium under relatively low Ra number conditions. In this
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
section, boundary condition with constant wall temperature is set for the left cavity wall. Sets of cavity sizes and over-temperatures and corresponding Ra numbers are listed in Table 3. In which, the density and viscosity of the PCM is taken at the mean temperature of its melting point and the temperature of the wall. The melting processes of all the cases are very similar to what have been shown in Fig. 2(b). Quantitative results are plotted in Fig. 4, which presents the melting fraction and heat flux versus time and their corresponding dimensionless forms. It should be clarified here that, in the legend of Fig. 4, W and H denote the width (mm) and height (mm) of the cavity and the temperature means over-temperature of the wall relative to the melting point of gallium. Comparing the computational results of cases 2, 4 and 6, which have the same width and over-temperature but different height, the melting fraction and heat flux curves of the 3 cases under real spatiotemporal scale almost coincide and thus conclusion can be drawn that the height of the cavity has no significant influence on the melting process of gallium. Expectedly, the higher the over-temperature, the faster the melting takes places and the larger the heat flux is, which is clearly shown by the curves of cases 3, 4 and 5. For situations which have the same over-temperature, (i.e. all the cases except cases 3 and 5), their heat flux variations almost fall on a same curve if the sharp drop at the end is neglected. Dimensionless analysis is performed for the 8 cases, as shown in Fig. 4(b) and (d). It is pleasing to find that good coincidence of all the cases is obtained. Hence there exist generalized dimensionless correlations to quantitatively describe such a melting process, which are fitted as
u ¼ 1:4ðFoSteÞ0:5
ð10Þ
Nu ¼ 0:727ðFoSteÞ0:5
905
ð11Þ
The fitting curves are also presented in Fig. 4(b) and (d). Noting that at the end of the melting process, significant deviation exists between the simulation results and the fitting curves, this is mainly caused by the aforementioned heat transfer inhibition by natural convection and thus induced sharp drop of heat flux at the end. 3.3. Dimensionless correlations: constant heat flux Boundary condition of constant heat flux for the left cavity wall is also studied. The geometry parameters kept the same as before and the corresponding Ra numbers are also listed in Table 3. In this situation, the characteristic temperature difference DT ¼ T w T m as appeared in the former definitions (Eqs. (7) and (9)) under constant wall temperature condition is now unknown, and it is replaced by q00 W=k in which q00 is known in current situation [18,22]. Hence, the dimensionless parameters are redefined as
s ¼ FoSte ¼ Ra ¼
q00 t DHqW
gbq2 cp W 4 q00 k
2
l
ð12Þ
ð13Þ
Here, Nu number is in the same form of Eq. (8). However, under current condition the heat flux q00 is already known while the characteristic temperature difference DT is obtained from the computational results. The main simulation results are presented in Fig. 5, in which Fig. 5(c) shows the mean temperature rise of the left cavity wall. Very similar conclusions can be drawn in this situation as for that
Fig. 5. Melting process of gallium in different cavities under different heat fluxes. (a) Melting fraction versus time; (b) melting fraction versus FoSte; (c) left wall temperature rise versus time; (d) Nu versus FoSte.
906
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
Fig. 6. Melting process of gallium in a rectangular cavity (20 mm 20 mm) under different over-temperatures. (a) Melting fraction versus time; (b) Nu versus FoSte; (c) melting fraction versus FoSteRa0.1; (d) Nu versus FoSteRa0.1.
of constant wall temperature. Some distinguishing features are that under constant heat flux condition, the melting fraction grows nearly linearly, as well as the wall temperature. As discussed in Section 3.1, at the end of the melting process, the heat transfer will be hindered. At this point, the growth of melting fraction slows down and much heat is stored in the liquid PCM in the form of sensible heat which rapidly raises its temperature as well as the wall, as shown in Fig. 5(c). Attentions should be paid to such a phenomenon in practical thermal design to avoid overheating of the component requiring cooling. Dimensionless analysis gives the corresponding correlations for constant heat flux conditions as
u ¼ 0:985FoSte
ð14Þ
Nu ¼ 1:08ðFoSteÞ0:958
ð15Þ
It is worthy to note that in fact, Eq. (14) can be approximately obtained from energy conservation equation. Take the PCM in the cavity as a whole, and neglect the sensible heat stored in the liquid phase, energy conservation law gives q00 t ¼ uqW DH, where the item on the left denotes the heat transferred from the left wall to the PCM at unit cavity height and length, while the right item presents the energy absorbed by the PCM in the form of latent heat. Rearrange this equation, one can get
u¼
q00 t at Wq00 ¼ 2 ¼ FoSte qW DH W aqDH
ð16Þ
The coefficient 0.985 in Eq. (14) is in fact a correction for Eq. (16) to take the sensible heat absorbed by the liquid PCM into consideration.
If carefully reviewing the dimensionless analysis in Figs. 4 (b), (d) and 5(b), (d), it can be found that significant derivation exists between case 8 and the remaining 7 cases. In fact, all the fitting curves above are based on the data of the first 7 cases and case 8 is excluded. Comparing the Ra number of all the 8 cases in Table 3, case 8 has the largest value, 4420 for constant wall temperature and 4376 for constant heat flux. Further increase the Ra number by increasing the width of the cavity or increasing the over-temperature or heat flux, more significant derivation will appear and thus the fitted dimensionless correlations are no longer applicable. The reason lies in that in the dimensionless analysis above, Ra number is out of account due to its negligible influence on the heat transfer process when it is at relatively low value. However, with the increase of Ra, the natural convection becomes more and more intensive, which will significantly affect the melting process and thus it cannot be ignored. The current numerical investigation suggests that the fitted correlations are applicable when the Ra number is below 2000–4000, both for constant wall temperature condition and constant heat flux condition. In the next section, a preliminary exploration will be performed to investigate the melting process of gallium with relatively high Ra number. 3.4. Discussion on high Ra number condition In this section, a larger cavity with dimension of 20 mm 20 mm is selected as an example for investigation. Constant wall temperature condition is set with 6 different overtemperatures (i.e. 5 K, 10 K, 15 K, 20 K, 25 K and 30 K), the Ra number of which are 11,788, 23,577, 35,365, 47,153, 58,941, 70,730, respectively. Dimensionless results of current configuration are
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
presented in Fig. 6(a) and (b). It can be found that the curves have poor coincidence when the Ra number is out of account. The larger the Ra number is, the shorter the melting time is and the larger Nu number is before the sharp drop at the end, which is much similar to that of eicosane. Similar to Refs. [18] and [22], in which traditional paraffin PCM is used for investigation of PCM based heat sink, the dimensionless time is redefined in the form of s ¼ FoSteRan to take the natural convection into consideration. It is found that n = 0.1 gives best coincidence of the dimensionless curves, as shown in Fig. 6(c) and (d).
4. Conclusion LMPMs are a new kind of PCMs which exhibit excellent heat extraction capability due to their high thermal conductivity, which are generally two orders of magnitude larger than that of conventional paraffin PCMs. Such a feature renders it distinctive melting mode which hardly exists in paraffin PCMs. Under relatively low Ra number (below 2000– 4000) conditions, the melting process of LMPMs shows great similarity and corresponding fitting dimensionless correlations can be given by u ¼ 1:4ðFoSteÞ0:5 , Nu ¼ 0:727ðFoSteÞ0:5 for constant wall temperature condition and u ¼ 0:985FoSte, Nu ¼ 1:08ðFoSteÞ0:958 for constant heat flux condition. As for higher Ra number condition, natural convection will significantly influence the melting process. A preliminary exploration is conducted in this paper to illustrate such an influence which takes the Ra number into consideration in the dimensionless analysis. More works are needed in the near future to further figure out the melting characteristics of LMPMs at high Ra number. The conclusions in this paper can serve as a valuable reference for future practical thermal design of LMPMs based thermal management systems and thermal energy storage systems. Acknowledgment This work is partially supported by the funding from Tsinghua University (No. 7131538) and Chinese Academy of Sciences. References [1] A.S. Fleischer, Thermal Energy Storage Using Phase Change Materials: Fundamentals and Applications, Springer, 2015. [2] F.L. Tan, C.P. Tso, Cooling of mobile electronic devices using phase change materials, Appl. Therm. Eng. 24 (2–3) (2004) 159–169. [3] R. Kandasamy, X.-Q. Wang, A.S. Mujumdar, Transient cooling of electronics using phase change material (PCM)-based heat sinks, Appl. Therm. Eng. 28 (8–9) (2008) 1047–1057.
907
[4] G. Setoh, F. Tan, S. Fok, Experimental studies on the use of a phase change material for cooling mobile phones, Int. Commun. Heat Mass Transfer 37 (9) (2010) 1403–1410. [5] H. Zhang, Q. Xu, Z. Zhao, J. Zhang, Y. Sun, L. Sun, F. Xu, Y. Sawada, Preparation and thermal performance of gypsum boards incorporated with microencapsulated phase change materials for thermal regulation, Sol. Energy Mater. Sol. Cells 102 (2012) 93–102. [6] R. Adinberg, D. Zvegilsky, M. Epstein, Heat transfer efficient thermal energy storage for steam generation, Energy Convers. Manage. 51 (1) (2010) 9–15. [7] N. Sarier, E. Onder, Organic phase change materials and their textile applications: an overview, Thermochim. Acta 540 (2012) 7–60. [8] T.D. Swanson, G.C. Birur, NASA thermal control technologies for robotic spacecraft, Appl. Therm. Eng. 23 (9) (2003) 1055–1065. [9] S. Krishnan, S.V. Garimella, S.S. Kang, A novel hybrid heat sink using phase change materials for transient thermal management of electronics, IEEE Trans. Compon. Packag. Technol. 28 (2) (2005) 281–289. [10] R. Akhilesh, A. Narasimhan, C. Balaji, Method to improve geometry for heat transfer enhancement in PCM composite heat sinks, Int. J. Heat Mass Transfer 48 (13) (2005) 2759–2770. [11] C.-Y. Zhao, W. Lu, Y. Tian, Heat transfer enhancement for thermal energy storage using metal foams embedded within phase change materials (PCMs), Sol. Energy 84 (8) (2010) 1402–1412. [12] Q. He, S. Wang, M. Tong, Y. Liu, Experimental study on thermophysical properties of nanofluids as phase-change material (PCM) in low temperature cool storage, Energy Convers. Manage. 64 (2012) 199–205. [13] M. Abdollahzadeh, M. Esmaeilpour, Enhancement of phase change material (PCM) based latent heat storage system with nano fluid and wavy surface, Int. J. Heat Mass Transfer 80 (2015) 376–385. [14] H. Ge, H. Li, S. Mei, J. Liu, Low melting point liquid metal as a new class of phase change material: an emerging frontier in energy area, Renew. Sustain. Energy Rev. 21 (2013) 331–346. [15] H. Ge, J. Liu, Keeping smartphones cool with gallium phase change material, J. Heat Transfer 135 (2013) 054503. [16] H. Ge, J. Liu, Cooling capacity of metal phase change material for thermal management, in: Proceedings of the ASME 2013 International Mechanical Engineering Congress and Exposition, 2013. [17] V. Voller, M. Cross, N. Markatos, An enthalpy method for convection/diffusion phase change, Int. J. Numer. Meth. Eng. 24 (1) (1987) 271–284. [18] S.K. Saha, P. Dutta, Heat transfer correlations for PCM-based heat sinks with plate fins, Appl. Therm. Eng. 30 (16) (2010) 2485–2491. [19] X.-Q. Wang, C. Yap, A.S. Mujumdar, A parametric study of phase change material (PCM)-based heat sinks, Int. J. Therm. Sci. 47 (8) (2008) 1055–1068. [20] E. Assis, L. Katsman, G. Ziskind, R. Letan, Numerical and experimental study of melting in a spherical shell, Int. J. Heat Mass Transfer 50 (9–10) (2007) 1790– 1804. [21] V. Shatikian, G. Ziskind, R. Letan, Numerical investigation of a PCM-based heat sink with internal fins, Int. J. Heat Mass Transfer 48 (17) (2005) 3689–3706. [22] V. Shatikian, G. Ziskind, R. Letan, Numerical investigation of a PCM-based heat sink with internal fins: constant heat flux, Int. J. Heat Mass Transfer 51 (5–6) (2008) 1488–1493. [23] A. Brent, V. Voller, K.t.J. Reid, Enthalpy-porosity technique for modeling convection-diffusion phase change: application to the melting of a pure metal, Numer. Heat Transfer, Part A: Appl. 13 (3) (1988) 297–318. [24] C. Beckermann, R. Viskanta, Effect of solid subcooling on natural convection melting of a pure metal, J. Heat Transfer 111 (2) (1989) 416–424. [25] C. Gau, R. Viskanta, Melting and solidification of a pure metal on a vertical wall, J. Heat Transfer 108 (1) (1986) 174–181. [26] M. Huang, P. Eames, B. Norton, Comparison of a small-scale 3D PCM thermal control model with a validated 2D PCM thermal control model, Sol. Energy Mater. Sol. Cells 90 (13) (2006) 1961–1972. [27] T.L. Bergman, F.P. Incropera, A.S. Lavine, Fundamentals of Heat and Mass Transfer, John Wiley & Sons, 2011.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Numerical investigation of the phase change process of low melting point metal Xiao-Hu Yang a, Si-Cong Tan a, Jing Liu a,b,⇑ a b
Beijing Key Lab of Cryo-Biomedical Engineering and Key Lab of Cryogenics, Technical Institute of Physics and Chemistry, Chinese Academy of Sciences, Beijing 100190, China Department of Biomedical Engineering, School of Medicine, Tsinghua University, Beijing 100084, China
a r t i c l e
i n f o
Article history: Received 1 April 2016 Received in revised form 28 April 2016 Accepted 30 April 2016 Available online 23 May 2016 Keywords: Phase change material Low melting point metal Parametric numerical investigation Dimensionless correlation
a b s t r a c t Low melting point metals (LMPMs) are a new kind of phase change materials (PCMs) which exhibit excellent heat extraction capability due to their high thermal conductivity. Such a characteristic renders it distinctive melting mode which hardly exists in conventional paraffin PCMs. In this paper, comprehensive numerical investigation is carried out to reveal the melting process of LMPMs. First, comparative simulations of typical paraffin eicosane and typical LMPM gallium is conducted. Then, dimensionless analysis is performed based on a series of parametric studies and corresponding dimensionless correlations are given under both constant wall temperature condition and constant heat flux condition. Further, the critical dimensionless values within which the correlations are applicable are suggested. Finally, a preliminary exploration is performed to illustrate the influence of natural convection under relatively high Ra number condition. This study can serve as a valuable reference for future practical thermal design of LMPMs based thermal management systems and thermal energy storage systems. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction Phase change materials (PCMs) are materials which undergo melting/solidification process at a constant or nearly constant temperature and absorb/release thermal energy from/to the surroundings. PCMs are widely used in the field of thermal management and thermal energy storage [1], such as electronics cooling [2–4], energy storage for buildings [5], solar energy systems [6], thermal comfortable textile design [7] and space systems [8]. Conventionally, PCMs are classified into organic PCMs and inorganic PCMs. Organic PCMs are generally used in relatively low temperature situations (from 35 °C to about 70 °C). Paraffin is a typical kind of organic PCM and has been extensively researched and used. The main drawback of paraffin PCM lies in its low thermal conductivity, which seriously hinders the heat conduction inside the PCM and thus decreases its efficiency. Many efforts have been made to improve such a situation, such as internal fin [3,9,10], metallic foam [11], and nano-particle inclusion [12,13]. However, such improvement is still limited. Ge and Liu [14] had ever proposed that using low melting point metals (LMPMs), typically gallium and its alloys and bismuth based ⇑ Corresponding author at: Technical Institute of Physics and Chemistry, Chinese Academy of Sciences, Beijing 100190, China. Tel.: +86 10 82543765; fax: +86 10 82543767. E-mail address: [email protected] (J. Liu). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.04.109 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
and indium based alloys, as a new kind of PCMs and successfully used it in the thermal management of smart phone [15,16] which exhibited good performance. The main advantages of LMPMs lie in their high thermal conductivity and high volumetric latent heat. Since the thermophysical properties of LMPMs are much different from that of conventional paraffin PCMs, some distinguishing features exist during their phase change process. However, up till now, there are rather limited literatures to reveal such a difference. This paper is dedicated to provide a comprehensive understanding about the phase change characteristics of LMPMs via parametric numerical investigation. The phase change process of PCMs includes melting and solidification. In this paper, attentions will be focused on the former one. First, comparative investigation of the melting process of typical paraffin eicosane and typical LMPM gallium is conducted, which reveals the superiority of using LMPMs as PCMs and their special characteristics. Then, dimensionless analysis of the melting process of LMPM is carried out based on a series of parametric study and corresponding dimensionless correlations are given under both constant wall temperature condition and constant heat flux condition. Further, the critical dimensionless values within which the correlations are applicable are suggested. Finally, preliminary exploration is performed to reveal the melting process of LMPM under relatively high Ra
900
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
Nomenclature A C cp Fo g H DH h k Nu p q00 Ra Ste T t u v W x, y
Greek letters thermal diffusivity (m2/s) b volume thermal expansivity (1/K) c local liquid fraction e small constant g dimensionless length along y direction h dimensionless temperature l dynamic viscosity (kg/s/m) n dimensionless length along x direction q mass density (kg/m3) u melting fraction s dimensionless time
source item coefficient in momentum equation constant coefficient of A specific heat capacity (J/kg/K) Fourier number gravitational acceleration (m/s2) height (m) fusion latent heat (J/kg) enthalpy (J/kg) thermal conductivity (W/m/K) Nusselt number pressure (N/m2) heat flux (W/m2) Rayleigh number Stefan number temperature (K) time (s) velocity along x direction (m/s) velocity along y direction (m/s) width (m) rectangular coordinates (m)
a
Subscripts c cold h hot l liquid phase m melting point ref reference s solid phase w wall
number condition where the influence of natural convection becomes significant.
avoiding zero denominator. In the current simulation, C = 105 and
e ¼ 0:001.
Energy equation: 2. Numerical method 2.1. Enthalpy-porosity method Enthalpy-porosity method, which was proposed by Voller et al. [17], is widely used for the numerical simulation of solid–liquid phase change problems [18–22]. The main features of this method are that: (1) The energy equation is based on the variable enthalpy but not temperature, which can conveniently take the phase change process into consideration; (2) the phase change region (mushy region) is treated as a porous domain, and the local liquid fraction of the PCM is regarded equal to the porosity. In enthalpyporosity method, the governing equations in two-dimensional Cartesian coordinate are given by: Continuity equation:
@ q @ðquÞ @ðqv Þ þ ¼0 þ @x @y @t
ð1Þ
Momentum equations with Boussinesq approximation invoked:
@ðquÞ @ðquuÞ @ðqv uÞ @p @ þ þ ¼ þ @t @x @y @x @x
l
@u @ @u þ l @x @y @y
þ Au @ðqv Þ @ðquv Þ @ðqvv Þ @p @ þ þ ¼ þ @t @x @y @x @x
l
@v @x
þ
þ Av þ qgbðT T m Þ
@ @y
l
@v @y
ð4Þ
8RT cp dT; T < T solidus > T > > < R ref T solidus cp dT þ cDH; T solidus 6 T 6 T liquidus h¼ T ref > > > : R T solidus c dT þ DH þ R T c dT; T > T liquidus p T liquidus p T ref
ð5Þ
where k and cp are the thermal conductivity and specific heat capacity of the PCM, respectively. h denotes the enthalpy of the PCM and is defined by Eq. (5), where, Tref is a reference temperature and the enthalpy at this temperature is set to be zero, DH is the fusion latent heat, Tsolidus and Tliquidus are the temperature limit at which phase change starts and ends up, respectively. Eq. (5) indicates that, in the mushy region, if Tsolidus = Tliquidus = Tm, then T = Tm; if T solidus – T liquidus , T can be calculated by the following definition
c¼
T T solidus T liquidus T solidus
ð6Þ
RT where c can be obtained from Eq. (5) as c ¼ h T refsolidus cp dT =DH. It
ð2Þ
@ðqhÞ @ðquhÞ @ðqv hÞ @ k @h @ k @h þ þ þ ¼ @t @x @y @x cp @x @y cp @y
ð3Þ
where g is the gravitational acceleration; q, l, b and Tm are the density, dynamic viscosity, volume thermal expansivity and melting point of the PCM; Au and Av are the momentum dissipation source items, which are used for suppressing the velocity in the mushy and solid region, hence it is a function of local liquid fraction c. Brent et al. [23] defined such a function in the form of A ¼ Cð1 cÞ2 =ðc3 þ eÞ, where C is a large constant number to enhance the suppression effect and e is a small constant set for
is worthy to note here that the term local liquid fraction c should be distinguished from the term melting fraction u, the latter of which means the ratio of the volume of the molten PCM to the volume of the initial solid PCM and this definition will be widely used in the following. Mathematically, u is the volume weighted average of c R over the whole cavity defined as u ¼ V cdV=V. 2.2. Validation of the numerical method Before parametric study, the numerical method is validated first by simulating the melting process of pure gallium in a rectangular cavity (47.6 mm 47.6 mm) which had ever been experimentally investigated by Beckermann et al. [24]. The cavity is heated at its left wall with overheating of 10.2 °C and cooled at the right wall with subcooling of 4.8 °C. Here, overheating and subcooling are relative
901
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907 Table 1 Thermophysical properties of gallium. Parameter
Value
Parameter
Value
Melting point Tm (K) Fusion latent heat DH (J/kg) Density (solid) qs (kg/m3) Density (liquid) ql (kg/m3) Specific heat capacity cps (solid) (J/kg/K)
302.93 80,160 5903.7 6094.7 372.3
Specific heat capacity (liquid) cpl (J/kg/K) Thermal conductivity (solid) ks (W/m/K) Thermal conductivity (liquid) kl (W/m/K) Expansion coefficient b (1/K) Viscosity l (kg/m/s)
397.6 33.49 33.6767 1.2 104 1.75 103
(a)
within which the natural convection is laminar [27]. As can be seen later that all of the cases investigated in this paper are laminar and hence laminar viscous flow mode is adopted in the numerical method for all of the simulations. SIMPLE algorithm is selected for velocity–pressure coupling and body force weighted scheme is used for the spatial discretization of pressure. Momentum and energy equations use second order upwind scheme to ensure accuracy. Convergence criteria of 104 is set for the momentum equation and 107 for the energy equation. Time step of 0.01 s is chosen for the transient computation with compromise of enough accuracy and time consumption. Grid independence study is also conducted with 6 sets of grid sizes (20 20, 30 30, 40 40, 50 50, 60 60, 80 80), the result of which indicate that 60 60 grid is accurate enough for this simulation since the difference of the calculated melting fraction and heat flux based on 60 60 grid and 80 80 grid are less than 0.5%. Typical computational results are presented in Fig. 1 in comparison with the experimental data in Ref. [24] in the form of dimensionless parameters. Fig. 1(a) shows the temperature distribution on three specific horizontal lines in the cavity and Fig. 1(b) exhibits the solid–liquid interface under 5 different boundary conditions. Herein, hl and hs are dimensionless temperature in the liquid and solid region of the PCM defined as hl ¼ ðT T m Þ=ðT h T m Þ and hs ¼ ðT T m Þ=ðT m T c Þ, where Th and Tc are the temperature of the left hot wall and right cold wall of the cavity, respectively; n and g are dimensionless length in horizontal and vertical direction defined as n ¼ x=W and g ¼ y=H, where W and H are the width and height of the cavity. It can be seen from Fig. 1(a) and (b) that good agreement between the current computational results and the experimental data is achieved, which implies the reliability of the current numerical method.
3. Results and discussions 3.1. Comparison with paraffin
(b) Fig. 1. Comparison of the current computational results and the experimental results in literature [24]. (a) Temperature profiles on 3 specified lines. (b) Liquid– solid interfaces under 5 different boundary conditions. Case 1: DT h ¼ T h T m ¼ 10:2 C, DT c ¼ T m T c ¼ 4:8 C; case 2: DT h ¼ 10:2 C, DT c ¼ 9:8 C; case 3: DT h ¼ 10:2 C, DT c ¼ 19:8 C; case 4: DT h ¼ 5:2 C, DT c ¼ 14:8 C; case 5: DT h ¼ 15:2 C, DT c ¼ 14:8 C.
to the melting point of gallium 29.78 °C (302.93 K). The top and bottom wall of the cavity are set to be adiabatic. Thermophysical properties of gallium refer to Ref. [25] and are listed in Table 1. Many literatures [21,26] have verified the accuracy of using two-dimensional simulation to compute the realistic threedimensional problem and in this paper, 2-D simulation is performed. Commercial codes Fluent 6.3 is used to carry out the simulation, which is running based on the aforementioned enthalpy-porosity method. Boussinesq approximation is applied to activate natural convection. The Ra number (Eq. (9)) of this simulation case is 3.24 105, which is much smaller than the critical value 109
Based on the aforementioned numerical method, a series of parametric study are performed. At first, a comparative investigation is conducted to illustrate the different melting process of LMPMs and conventional paraffin PCMs. Typical LMPM pure gallium and typical paraffin eicosane are chosen here for comparative simulation and the main thermophysical properties of which are listed in Table 2. The melting process happened in a rectangular cavity with dimension of 5 mm in width and 20 mm in height. The top, right and bottom wall of the cavity are set to be adiabatic. Constant temperature condition is set for the left wall of the cavity with 3 different over-temperature sets (i.e. 5 K, 15 K and 25 K). Here, over-temperature indicates the temperature relative to the melting point of the PCM. Transient simulation is carried out and the initial temperature of the whole PCM is set to be at its melting point. To intuitively illustrate the different melting modes of eicosane and gallium, liquid fraction contours and over-temperature contours at different melting fractions (i.e. u ¼ 0:3; 0:6 and 0:9) for the two PCMs are presented, in Fig. 2(a) for eicosane and
902
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
Table 2 Thermophysical properties of gallium and eicosane.
a
PCM
qa (kg/m3)
cp (J/kg/K)
l (103 kg/m/s)
k (W/m/K)
b (1/K)
T m (K)
DH (J/kg)
Gallium Eicosane [19]
6094.7 785
397.6 2460
1.75 exp(1790/T 4.25) [21]
33.6767 0.15
0.00012 0.001
302.93 309.15
80,160 247,000
The density of liquid phase PCM at its melting point.
Fig. 2. Liquid fraction and temperature distribution in eicosane and gallium under different melting fractions. Herein, over-temperature indicates temperature relative to the melting point of the PCM. (a) Eicosane; (b) gallium.
Fig. 2(b) for gallium. The streamline and flow direction of the natural convective flow of liquid PCM at melting fraction of 0.6 and 0.9 are also inserted in the liquid fraction contours. It can be seen that for eicosane, the liquid–solid interface is inclined from bottom-left to top-right, especially under large melting fraction conditions, which is caused by the clockwise natural convective flow and thus induced heat transfer enhancement in the upper cavity region. While in the lower region, heat transfer is mainly dependent on heat conduction which is relatively weak thus the PCM there melts more slowly. The temperature distribution in the liquid region of eicosane exhibits horizontal slicing which is a typical feature of natural convection. While for gallium the interface is nearly vertical all the time, as well as its temperature distribution, which is a typical feature of heat conduction. These phenomena qualitatively indicate that for eicosane, natural convection plays an important role in the total heat transfer, while for gallium heat conduction is dominant. Fig. 2 also presents the velocity distribution in the liquid PCM region. It can be found that in fact, the natural convective flow in liquid gallium is more intensive since the velocity of which is several times larger than that in eicosane. However, such an intensive flow does not make noticeable contribution to the total heat transfer. The reason causing the different melting modes between gallium and eicosane mainly lies in their large difference in thermal conductivity. As shown in Table 2, the thermal conductivity of gallium is more than 200 times larger than that of eicosane. Hence, in eicosane, the heat conduction is very weak, and natural convection will significantly enhance the heat transfer. While in gallium, heat
conduction is so strong that the natural convective heat transfer becomes relatively weak in this case though the convection is more intensive. Fig. 3 quantitatively presents the melting fraction curves and heat flux curves over time of eicosane (Fig. 3(a) and (b)) and gallium (Fig. 3(c) and (d)) at different over-temperatures. It can be seen from Fig. 3(a)–(d) that, for all cases, a very high heat flux exists at the beginning of the melting due to the large temperature difference between the cavity wall and adjacent PCM. Quickly, a liquid film is formed between the wall and the solid PCM. Such a liquid region is in fact a thermal resistance for the heat transfer from the wall to the solid PCM. With time going on, the liquid region becomes thicker and thicker and corresponding thermal resistance keeps getting larger, which leads to the gradual decrease of the heat flux as well as the growth rate of the melting fraction. At the same time, natural convection in the liquid region is developed gradually due to the temperature difference and thus induced density difference of the liquid PCM, which promotes the heat transfer in some degree. However, at the end of the melting process, there generally exists a sharp drop of the heat flux, as can be seen in Fig. 3(b) and (d). Such a drop is inferred to also be caused by the natural convection. The reason lies in that at the end, only small amount of solid PCM exist at the bottom-right corner of the cavity, while the liquid flows clockwise at the bottom, namely from right to left (Fig. 2), which is opposite to the heat flux direction and thus the heat transfer is adversely hindered at this point. The larger the over-temperature is, the more sharply the heat flux decreases.
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
903
Fig. 3. Melting fraction and heat flux versus time of different PCMs under different overheating temperatures. Dimension of the cavity is 5 mm 20 mm (width height). (a) Eicosane, melting fraction; (b) eicosane, heat flux; (c) gallium, melting fraction; (d) gallium, heat flux; (e) melting fraction versus FoSte; (f) Nu versus FoSte.
It can also be seen that the time needed to completely melt eicosane are 515 s, 910 s and 3080 s respectively under overtemperatures of 25 K, 15 K and 5 K. While for gallium, the time are only 8.5 s, 13.5 s and 38 s, respectively. The heat flux in the cases with gallium as PCM is nearly two orders of magnitude larger than that with eicosane as PCM at the same over-temperature. Hence it can be concluded that gallium has much superior heat extraction efficiency as a PCM than eicosane. Dimensionless analysis for the melting process is conducted for both eicosane and gallium, as shown in Fig. 3(e) and (f). Here, Fourier number is dimensionless time, defined as Fo ¼ at=W 2 , where a is the thermal diffusivity of the PCM and W denotes the width of the cavity which is regarded as the characteristic length here. Stefan number is generally used for phase change problem in the form of Ste ¼ cp DT=DH, where DT is the characteristic temperature difference and defined here as DT ¼ T w T m , T w is the left
wall temperature of the cavity. To take into account both the transient heat conduction and phase change, the product of Fo number and Ste number can serve as a dimensionless time [21,24], namely
s ¼ FoSte ¼
acp DTt DHW 2
ð7Þ
Nu number is used to describe the dimensionless heat transfer intensity, defined as
Nu ¼
W q00 k DT
ð8Þ
where q00 denotes the heat flux on the left wall. Fig. 3(e) shows the melting fraction curve versus FoSte and Fig. 3(f) is Nu versus FoSte. It can be seen that both in Fig. 3(e) and (f) the dimensionless curves of gallium under 3 different overtemperatures are almost coincident. While for eicosane, significant
904
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
Table 3 Sets of cavity sizes and boundary conditions and corresponding Ra numbers. Case
Case 1 Case 2 Cases 3/4/5 Case 6 Case 7 Case 8
Width (mm)
3 5
8 10
Height (mm)
10 10 20 40 30 20
Constant wall temperature
Constant heat flux
Over-temperature (K)
Ra
Heat flux (kW/m2)
Ra
15 15 5/15/25 15 15 15
119 552 184/552/921 552 2263 4420
50 50 30/50/100 50 50 50
35 273 164/273/547 273 1792 4376
Fig. 4. Melting process of gallium in different cavities under different wall temperatures. W and H mean the width (mm) and height (mm) of the cavity and the temperature means over-temperature of the wall relative to the melting point of the PCM, respectively. (a) Melting fraction versus time; (b) melting fraction versus FoSte; (c) heat flux versus time; (d) Nu versus FoSte.
difference exists. The higher the wall temperature is, the shorter the dimensionless time for complete fusion is and the higher Nu number is. Such a difference is essentially a quantitative reflection of the aforementioned influence of natural convection. Rayleigh number is generally used in natural convection problems which can be viewed as the ratio of buoyancy and viscosity forces multiplied by the ratio of momentum and thermal diffusivities, defined as
Ra ¼
gbDTcp q2 W 3 lk
ð9Þ
In fact, Ra number can also be understood as the ratio between the rate of convective heat transfer to the rate of heat conduction, up to a numerical factor. Under the current cavity geometry and an overtemperature of 15 K, the Ra number for eicosane and gallium are 45,659 and 552 respectively, which quantitatively indicates that significant natural convective heat transfer relative to heat conduction exists in eicosane while it is rather weak in gallium. Natural convection enhances the heat transfer and thus increases Nu
number and decreases the dimensionless time for complete fusion. It is worthy to note that, although gallium exhibits inferior heat transfer capacity than eicosane under the dimensionless scale (Fig. 3(e) and (f)), it is much superior under real spatiotemporal scale (Fig. 3(a) to (d)). Up to this point, the reason why the dimensionless curves of gallium under different conditions coincide is very clear: the natural convective heat transfer in gallium is relatively weak compared with the heat conduction and it can be neglected, thus the dimensionless analysis has no relationship with Ra number. In fact for gallium, there exist tremendous situations where the influence of natural convection can be neglected due to its high thermal conductivity which hardly exists in conventional paraffin PCMs. In the following, attentions will be focused on these situations. 3.2. Dimensionless correlations: constant wall temperature Parametric study is performed to investigate the melting process of gallium under relatively low Ra number conditions. In this
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
section, boundary condition with constant wall temperature is set for the left cavity wall. Sets of cavity sizes and over-temperatures and corresponding Ra numbers are listed in Table 3. In which, the density and viscosity of the PCM is taken at the mean temperature of its melting point and the temperature of the wall. The melting processes of all the cases are very similar to what have been shown in Fig. 2(b). Quantitative results are plotted in Fig. 4, which presents the melting fraction and heat flux versus time and their corresponding dimensionless forms. It should be clarified here that, in the legend of Fig. 4, W and H denote the width (mm) and height (mm) of the cavity and the temperature means over-temperature of the wall relative to the melting point of gallium. Comparing the computational results of cases 2, 4 and 6, which have the same width and over-temperature but different height, the melting fraction and heat flux curves of the 3 cases under real spatiotemporal scale almost coincide and thus conclusion can be drawn that the height of the cavity has no significant influence on the melting process of gallium. Expectedly, the higher the over-temperature, the faster the melting takes places and the larger the heat flux is, which is clearly shown by the curves of cases 3, 4 and 5. For situations which have the same over-temperature, (i.e. all the cases except cases 3 and 5), their heat flux variations almost fall on a same curve if the sharp drop at the end is neglected. Dimensionless analysis is performed for the 8 cases, as shown in Fig. 4(b) and (d). It is pleasing to find that good coincidence of all the cases is obtained. Hence there exist generalized dimensionless correlations to quantitatively describe such a melting process, which are fitted as
u ¼ 1:4ðFoSteÞ0:5
ð10Þ
Nu ¼ 0:727ðFoSteÞ0:5
905
ð11Þ
The fitting curves are also presented in Fig. 4(b) and (d). Noting that at the end of the melting process, significant deviation exists between the simulation results and the fitting curves, this is mainly caused by the aforementioned heat transfer inhibition by natural convection and thus induced sharp drop of heat flux at the end. 3.3. Dimensionless correlations: constant heat flux Boundary condition of constant heat flux for the left cavity wall is also studied. The geometry parameters kept the same as before and the corresponding Ra numbers are also listed in Table 3. In this situation, the characteristic temperature difference DT ¼ T w T m as appeared in the former definitions (Eqs. (7) and (9)) under constant wall temperature condition is now unknown, and it is replaced by q00 W=k in which q00 is known in current situation [18,22]. Hence, the dimensionless parameters are redefined as
s ¼ FoSte ¼ Ra ¼
q00 t DHqW
gbq2 cp W 4 q00 k
2
l
ð12Þ
ð13Þ
Here, Nu number is in the same form of Eq. (8). However, under current condition the heat flux q00 is already known while the characteristic temperature difference DT is obtained from the computational results. The main simulation results are presented in Fig. 5, in which Fig. 5(c) shows the mean temperature rise of the left cavity wall. Very similar conclusions can be drawn in this situation as for that
Fig. 5. Melting process of gallium in different cavities under different heat fluxes. (a) Melting fraction versus time; (b) melting fraction versus FoSte; (c) left wall temperature rise versus time; (d) Nu versus FoSte.
906
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
Fig. 6. Melting process of gallium in a rectangular cavity (20 mm 20 mm) under different over-temperatures. (a) Melting fraction versus time; (b) Nu versus FoSte; (c) melting fraction versus FoSteRa0.1; (d) Nu versus FoSteRa0.1.
of constant wall temperature. Some distinguishing features are that under constant heat flux condition, the melting fraction grows nearly linearly, as well as the wall temperature. As discussed in Section 3.1, at the end of the melting process, the heat transfer will be hindered. At this point, the growth of melting fraction slows down and much heat is stored in the liquid PCM in the form of sensible heat which rapidly raises its temperature as well as the wall, as shown in Fig. 5(c). Attentions should be paid to such a phenomenon in practical thermal design to avoid overheating of the component requiring cooling. Dimensionless analysis gives the corresponding correlations for constant heat flux conditions as
u ¼ 0:985FoSte
ð14Þ
Nu ¼ 1:08ðFoSteÞ0:958
ð15Þ
It is worthy to note that in fact, Eq. (14) can be approximately obtained from energy conservation equation. Take the PCM in the cavity as a whole, and neglect the sensible heat stored in the liquid phase, energy conservation law gives q00 t ¼ uqW DH, where the item on the left denotes the heat transferred from the left wall to the PCM at unit cavity height and length, while the right item presents the energy absorbed by the PCM in the form of latent heat. Rearrange this equation, one can get
u¼
q00 t at Wq00 ¼ 2 ¼ FoSte qW DH W aqDH
ð16Þ
The coefficient 0.985 in Eq. (14) is in fact a correction for Eq. (16) to take the sensible heat absorbed by the liquid PCM into consideration.
If carefully reviewing the dimensionless analysis in Figs. 4 (b), (d) and 5(b), (d), it can be found that significant derivation exists between case 8 and the remaining 7 cases. In fact, all the fitting curves above are based on the data of the first 7 cases and case 8 is excluded. Comparing the Ra number of all the 8 cases in Table 3, case 8 has the largest value, 4420 for constant wall temperature and 4376 for constant heat flux. Further increase the Ra number by increasing the width of the cavity or increasing the over-temperature or heat flux, more significant derivation will appear and thus the fitted dimensionless correlations are no longer applicable. The reason lies in that in the dimensionless analysis above, Ra number is out of account due to its negligible influence on the heat transfer process when it is at relatively low value. However, with the increase of Ra, the natural convection becomes more and more intensive, which will significantly affect the melting process and thus it cannot be ignored. The current numerical investigation suggests that the fitted correlations are applicable when the Ra number is below 2000–4000, both for constant wall temperature condition and constant heat flux condition. In the next section, a preliminary exploration will be performed to investigate the melting process of gallium with relatively high Ra number. 3.4. Discussion on high Ra number condition In this section, a larger cavity with dimension of 20 mm 20 mm is selected as an example for investigation. Constant wall temperature condition is set with 6 different overtemperatures (i.e. 5 K, 10 K, 15 K, 20 K, 25 K and 30 K), the Ra number of which are 11,788, 23,577, 35,365, 47,153, 58,941, 70,730, respectively. Dimensionless results of current configuration are
X.-H. Yang et al. / International Journal of Heat and Mass Transfer 100 (2016) 899–907
presented in Fig. 6(a) and (b). It can be found that the curves have poor coincidence when the Ra number is out of account. The larger the Ra number is, the shorter the melting time is and the larger Nu number is before the sharp drop at the end, which is much similar to that of eicosane. Similar to Refs. [18] and [22], in which traditional paraffin PCM is used for investigation of PCM based heat sink, the dimensionless time is redefined in the form of s ¼ FoSteRan to take the natural convection into consideration. It is found that n = 0.1 gives best coincidence of the dimensionless curves, as shown in Fig. 6(c) and (d).
4. Conclusion LMPMs are a new kind of PCMs which exhibit excellent heat extraction capability due to their high thermal conductivity, which are generally two orders of magnitude larger than that of conventional paraffin PCMs. Such a feature renders it distinctive melting mode which hardly exists in paraffin PCMs. Under relatively low Ra number (below 2000– 4000) conditions, the melting process of LMPMs shows great similarity and corresponding fitting dimensionless correlations can be given by u ¼ 1:4ðFoSteÞ0:5 , Nu ¼ 0:727ðFoSteÞ0:5 for constant wall temperature condition and u ¼ 0:985FoSte, Nu ¼ 1:08ðFoSteÞ0:958 for constant heat flux condition. As for higher Ra number condition, natural convection will significantly influence the melting process. A preliminary exploration is conducted in this paper to illustrate such an influence which takes the Ra number into consideration in the dimensionless analysis. More works are needed in the near future to further figure out the melting characteristics of LMPMs at high Ra number. The conclusions in this paper can serve as a valuable reference for future practical thermal design of LMPMs based thermal management systems and thermal energy storage systems. Acknowledgment This work is partially supported by the funding from Tsinghua University (No. 7131538) and Chinese Academy of Sciences. References [1] A.S. Fleischer, Thermal Energy Storage Using Phase Change Materials: Fundamentals and Applications, Springer, 2015. [2] F.L. Tan, C.P. Tso, Cooling of mobile electronic devices using phase change materials, Appl. Therm. Eng. 24 (2–3) (2004) 159–169. [3] R. Kandasamy, X.-Q. Wang, A.S. Mujumdar, Transient cooling of electronics using phase change material (PCM)-based heat sinks, Appl. Therm. Eng. 28 (8–9) (2008) 1047–1057.
907
[4] G. Setoh, F. Tan, S. Fok, Experimental studies on the use of a phase change material for cooling mobile phones, Int. Commun. Heat Mass Transfer 37 (9) (2010) 1403–1410. [5] H. Zhang, Q. Xu, Z. Zhao, J. Zhang, Y. Sun, L. Sun, F. Xu, Y. Sawada, Preparation and thermal performance of gypsum boards incorporated with microencapsulated phase change materials for thermal regulation, Sol. Energy Mater. Sol. Cells 102 (2012) 93–102. [6] R. Adinberg, D. Zvegilsky, M. Epstein, Heat transfer efficient thermal energy storage for steam generation, Energy Convers. Manage. 51 (1) (2010) 9–15. [7] N. Sarier, E. Onder, Organic phase change materials and their textile applications: an overview, Thermochim. Acta 540 (2012) 7–60. [8] T.D. Swanson, G.C. Birur, NASA thermal control technologies for robotic spacecraft, Appl. Therm. Eng. 23 (9) (2003) 1055–1065. [9] S. Krishnan, S.V. Garimella, S.S. Kang, A novel hybrid heat sink using phase change materials for transient thermal management of electronics, IEEE Trans. Compon. Packag. Technol. 28 (2) (2005) 281–289. [10] R. Akhilesh, A. Narasimhan, C. Balaji, Method to improve geometry for heat transfer enhancement in PCM composite heat sinks, Int. J. Heat Mass Transfer 48 (13) (2005) 2759–2770. [11] C.-Y. Zhao, W. Lu, Y. Tian, Heat transfer enhancement for thermal energy storage using metal foams embedded within phase change materials (PCMs), Sol. Energy 84 (8) (2010) 1402–1412. [12] Q. He, S. Wang, M. Tong, Y. Liu, Experimental study on thermophysical properties of nanofluids as phase-change material (PCM) in low temperature cool storage, Energy Convers. Manage. 64 (2012) 199–205. [13] M. Abdollahzadeh, M. Esmaeilpour, Enhancement of phase change material (PCM) based latent heat storage system with nano fluid and wavy surface, Int. J. Heat Mass Transfer 80 (2015) 376–385. [14] H. Ge, H. Li, S. Mei, J. Liu, Low melting point liquid metal as a new class of phase change material: an emerging frontier in energy area, Renew. Sustain. Energy Rev. 21 (2013) 331–346. [15] H. Ge, J. Liu, Keeping smartphones cool with gallium phase change material, J. Heat Transfer 135 (2013) 054503. [16] H. Ge, J. Liu, Cooling capacity of metal phase change material for thermal management, in: Proceedings of the ASME 2013 International Mechanical Engineering Congress and Exposition, 2013. [17] V. Voller, M. Cross, N. Markatos, An enthalpy method for convection/diffusion phase change, Int. J. Numer. Meth. Eng. 24 (1) (1987) 271–284. [18] S.K. Saha, P. Dutta, Heat transfer correlations for PCM-based heat sinks with plate fins, Appl. Therm. Eng. 30 (16) (2010) 2485–2491. [19] X.-Q. Wang, C. Yap, A.S. Mujumdar, A parametric study of phase change material (PCM)-based heat sinks, Int. J. Therm. Sci. 47 (8) (2008) 1055–1068. [20] E. Assis, L. Katsman, G. Ziskind, R. Letan, Numerical and experimental study of melting in a spherical shell, Int. J. Heat Mass Transfer 50 (9–10) (2007) 1790– 1804. [21] V. Shatikian, G. Ziskind, R. Letan, Numerical investigation of a PCM-based heat sink with internal fins, Int. J. Heat Mass Transfer 48 (17) (2005) 3689–3706. [22] V. Shatikian, G. Ziskind, R. Letan, Numerical investigation of a PCM-based heat sink with internal fins: constant heat flux, Int. J. Heat Mass Transfer 51 (5–6) (2008) 1488–1493. [23] A. Brent, V. Voller, K.t.J. Reid, Enthalpy-porosity technique for modeling convection-diffusion phase change: application to the melting of a pure metal, Numer. Heat Transfer, Part A: Appl. 13 (3) (1988) 297–318. [24] C. Beckermann, R. Viskanta, Effect of solid subcooling on natural convection melting of a pure metal, J. Heat Transfer 111 (2) (1989) 416–424. [25] C. Gau, R. Viskanta, Melting and solidification of a pure metal on a vertical wall, J. Heat Transfer 108 (1) (1986) 174–181. [26] M. Huang, P. Eames, B. Norton, Comparison of a small-scale 3D PCM thermal control model with a validated 2D PCM thermal control model, Sol. Energy Mater. Sol. Cells 90 (13) (2006) 1961–1972. [27] T.L. Bergman, F.P. Incropera, A.S. Lavine, Fundamentals of Heat and Mass Transfer, John Wiley & Sons, 2011.