Lattice Boltzmann modeling of melting of phase change materials in porous media with conducting fins

Accepted Manuscript Lattice Boltzmann modeling of melting of phase change materials in porous media with conducting fins Dongyan Gao, Zhenqian Chen, Dongliang Zhang, Linghai Chen PII: DOI: Reference:

S1359-4311(16)32362-6 http://dx.doi.org/10.1016/j.applthermaleng.2017.03.002 ATE 10003

To appear in:

Applied Thermal Engineering

Received Date: Revised Date: Accepted Date:

11 October 2016 12 January 2017 1 March 2017

Please cite this article as: D. Gao, Z. Chen, D. Zhang, L. Chen, Lattice Boltzmann modeling of melting of phase change materials in porous media with conducting fins, Applied Thermal Engineering (2017), doi: http://dx.doi.org/ 10.1016/j.applthermaleng.2017.03.002

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Lattice Boltzmann modeling of melting of phase change materials in porous media with conducting fins

Dongyan Gao a , Zhenqian Chen b, Dongliang Zhang a , Linghai Chen a a

School of Energy and Power Engineering, Nanjing Institute of Technology, Nanjing 211167, China b

School of Energy and Environment, Southeast University, Nanjing 210096, China

Abstract: Based on the averaged energy equation in terms of enthalpy, a modified lattice Boltzmann (LB) model is proposed for simulating melting of phase change materials in porous media with a conducting fin. Different from previous LB models, the present model incorporates the total enthalpy and a free parameter into the equilibrium distribution function, and thus makes it suitable for modeling conjugate heat transfer and has high computational efficiency by avoiding iteration procedure for phase change. To ensure the numerical stability, the multiple-relaxation-time collision scheme is adopted in the model. The model is validated by three benchmark problems. It is found that the numerical results are in good agreement with the analytical solutions and numerical results. In addition, the effect of the conducting fin on the melting processes in a porous cavity is investigated. The numerical results indicate that the melting heat transfer can be further enhanced by adding a conducting fin into the porous medium. The melting speed increases as the length of the fin increases and the heat capacity of the fin decreases. However, the vertical position of the fin has no remarkable impact on the melting speed.

Keywords: lattice Boltzmann method; porous media; melting; conjugate heat transfer; fins



Corresponding author.

Tel.: +86 25 8611 8355; E-mail address: [email protected] (D. GAO).

1. Introduction Phase change materials (PCMs) are used in numerous engineering applications such as latent thermal energy storage system and thermal management system due to their large latent heats and capabilities of maintaining nearly constant temperature [1-3]. However, the performance of most PCMs in thermal systems has been seriously limited by their low thermal conductivities and poor heat transfer performance. Low thermal conductivity of the PCM will lead to small heat transfer rate, resulting in slow heat charging and discharging in thermal energy storage system and inability to maintain constant temperature in the thermal management systems. Therefore, many methods have been proposed to improve the thermal conductivities and heat transfer performance of PCMs [1-2]. A common and promising method to enhance the thermal conductivity of PCMs is using porous media materials with high thermal conductivity as enhancer, such as porous metal foam [4-7] and porous carbon materials [8], which can form the effective heat transfer network in the fabricated composite PCM. Besides porous media, the use of solid fins as extended surface is the most straightforward approach [9-11]. It is evident that the addition of both fins and different porous media structures can enhance the internal heat transfer of PCMs. However, few studies involving the heat transfer enhancement of PCMs by combining these two technologies have been conducted in the literature. Recently, Feng et al. [12] have numerically investigated the heat transfer performance of PCMs in finned metal foams. Their results indicated that the heat transfer of PCMs in porous metal foam could be further improved by incorporating solid fins into the metal foam. Thus, there is a need for deeper understanding of heat transfer processes of PCMs in porous media with conducting fins. This study aims to deal with these issues using numerical method. To the present, it is still a challenging task to efficiently simulate the melting in porous media because of the inherent nonlinearities of melting process coupled with the fluid flow, as well as the physical complexities of porous media system. In the past several decades, various traditional numerical methods, such as the finite volume method, the finite element method and the finite difference method, have been adopted to investigate melting problems in porous media [13-17]. As an alternative and powerful numerical method for simulating fluid flows and complex physical processes [18-26], the lattice Boltzmann (LB) method has also been applied to the simulation of melting processes in porous media

[27-33]. Generally, the existing LB methods for melting in porous media can be classified into two major categories, i.e., the pore scale approach [27] and the representative elementary volume (REV) scale approach [28-33]. At the pore scale, Chen et al. applied thermal LB model with doubled populations to simulate the two-dimensional melting processes of phase change material in open-cell metal foams, the numerical results agree reasonably well with the experimental observation [27]. However, the pore-scale simulation requires much more computational resources. Thus more LB numerical studies for melting in porous media are based on the REV scale. Gao and Chen proposed a LB model to simulate the melting with natural convection in porous media based on the local thermal equilibrium (LTE) assumption [28]. Gao et al. developed a thermal LB model with multi-distribution function for natural convection with phase change under the local thermal non-equilibrium condition [29]. Gao et al. have investigated the melting processes of phase change materials in metal foam using LB method [30]. To improve the numerical stability, Liu and He developed a multiple-relaxation time (MRT) LB model for simulating solid-liquid phase change in porous media [31]. Tao et al. applied LB method to model the melting processes in the metal foam/paraffin composite PCMs [32]. To avoid iteration procedure for treatment phase change, Wu et al. proposed a LB model to simulate solid-liquid phase change problems in porous media [33]. Among the above mentioned studies, most of them are based on the assumption of local thermal equilibrium between solid matrix and PCM phase, because this assumption of LTE may be applicable for many practical applications [12]. The above literature survey shows that the LB method have gained increasing attention on simulation of the melting problems in porous media. However, all of the above existing LB studies [28-33] are focused on the melting of phase change materials saturated in porous media without solid fins, while simulation of the melting of PCMs in porous media with conducting fins using LB method has not been considered in the literature. In the presence of conducting fins, the melting phenomenon in porous media becomes more complicated. In fact, this is a conjugate heat transfer problem for the case with fins. When dealing with conjugate heat transfer, one needs to consider continuity conditions at the interface of two phases. The objective of the present study is to extend the lattice Boltzmann method (LBM) to deal with

melting in porous media with conducting fins at the REV scale, as well as investigating the effect of fins on the melting heat transfer of PCMs in porous media. In this study, the assumption of local thermal equilibrium between solid matrix and PCM phase is used, and thus the one-temperature model describing heat transfer in porous media is adopted. To more effectively simulate conjugate heat transfer as well as phase change, the total enthalpy of the mixture system, the volumetric heat capacity and a free parameter are introduced into the equilibrium distribution function and thus a modified LB formulation is proposed to handle the difference in both thermal conductivity and volumetric heat capacity at the porous/solid interface. In addition, by using a modified linear equilibrium distribution and an additional discrete source term in the evolution equation for the temperature field, a general form of energy equation can be correctly recovered through the Chapman-Enskog procedure without the deviation terms. To improve the numerical stability, the multiple-relaxation-time collision scheme is adopted. 2. Problem description and macroscopic governing equations As shown in Fig. 1, the melting problem in porous media with a conducting fin is considered. A two dimensional cavity of height H and width L is filled with a porous medium that is saturated by a phase change material (PCM). A conducting fin is attached along a wall of the cavity. The top and bottom walls are adiabatic. The left and right walls are kept at the uniform temperature, Th and Tc, respectively. At all walls of the cavity, the zero velocity boundary condition is imposed. Initially, the PCM and porous matrix and solid fin are at equilibrium at temperature Tinit. At the time t=0, the temperature of the left surface is raised to Th, the PCM begins to melt. In the melting process, the fluid flow in the liquid PCM induced by the buoyancy force is considered. Under the assumptions that the properties of porous medium and PCM are homogeneous and isotropic, porous matrix and the solid PCM are rigid, the liquid PCM is Boussinesq incompressible and the properties of the liquid and solid phases are constant, and the densities of the solid and liquid phases of the PCM are nearly the same, the continuity equation and the Brinkmann-Forchheimer equation for fluid flow in porous media at REV scale can be written as [13] u  0 ,

u u 1  (u  )( )   (p)  e2u  F , t  

(1) (2)

where ρ, u and p represent the volume-average density, velocity and pressure temperature, respectively. νe

is the effective viscosity. F is the total body force due to the presence of porous medium and other external force, and can be given as

F 

 f K

u

F K

u u  G ,

(3)

where νf is the kinematic viscosity of the fluid, K is the permeability, and Fε is the Forchheimer form coefficient. G is the external body force and further with Boussinesq approximation, which can be written as G   g (T  Tref )  a ,

(4)

where T and Tref refer to the average temperature and the reference temperature, respectively, β is thermal expansion coefficient,

g and a are the gravity acceleration and the acceleration induced by other

external force, respectively. Based on LTE assumption, the averaged energy equation in terms of the enthalpy for the porous matrix/solid/liquid mixture can be written as [13]

 [f l  f hf   (1  f l )  s hs  (1   )  m hm ]    (  f hf u)    (k e T ) , t

(5)

where, the subscripts f, s and m refer to properties of the fluid phase and solid phase of interstitial phase change substance, and the solid matrix, respectively. ke and T are the effective thermal conductivity and the temperature of the mixture, respectively. fl is the liquid volume fraction of the PCM. h is the total enthalpy of each component, which is approximately considered as the sum of sensible enthalpy and latent heat enthalpy for the phase change problems [13], and thus can be written as:

hf  c p,f T  La , hs  c p,sT , hm  c p,mT ,

(6)

where cp and La are the specific heat and the latent heat of constituent, respectively. By substituting equation (6) into equation (5), after a few steps, we get the following result [13]  [ ( f l  f c p,f T  (1  f l )  s c p,sT  f l  f La )  (1   )  m c p,mTm ] t    (  f c p,f Tu)    (k e T )

.

(7)

In order to simplify the form of equation (7), we introduce the following two parameters: hsys  hpcm  (1   ) mc p,mT , h pcm  [ f l  f c p,f  (1  f l )  s c p,s ]T  f l  f La .

(8)

And thus obtain

hs y s    ( f c p,f Tu)    (keT ) , t

(9)

Here, hsys and hpcm represent the total enthalpies of the porous matrix/solid/liquid mixture system and PCM, respectively. Note that there is an alternative form of the energy equation (7) frequently adopted in the previous studies and written as



T h f l    (Tu)    ( e T )   , t c p,f t

(10)

where αe is the apparent thermal diffusivity defined as  e  ke /(c p ) f . σ is the ratio of the mean thermal capacitance of the mixture to thermal capacitance of the melted fluid, which is defined as



[f l  f c p,f   (1  f l )  s c p,s  (1   )  m c p,m ]

 f c p ,f

.

(11)

Meanwhile, the transport process within the solid fin can be represented as

hfin    (k finT ) , t

(12)

where the enthalpy hfin=(ρcp)finT and the subscripts fin refers to properties of the solid fin. Note that equation (12) can be considered as a special case of equation (9) with the zero velocity. In order to close equations (1), (2) and (9), some parameters including ke and Fε need to be determined. Since the geometry of porous medium is considerably complex, workable approximations have been derived based on simpler models. The following different empirical equations for the effective conductivities ke are used in this work [13, 34],

ke  

k m  k pcm k pcm1/ 3

k e1/ 3  k m  0 ,

(13a)

ke  k p c m (1   )k m , k e  0.35[k pcm  (1   )k m ] 

(13b)

1  0.35 , 1  (  ) k pcm k m



(13c)

where k pcm  f l kf  (1  f l )ks . The correlations for the geometric function Fε are used as follows [13, 34]:

Fε 

Fε  0.00212(1   ) 0.132[1.18

1.75 150 3

,

1  1 ]1.63 . 3 1  e (1 ) / 0.04

(14a)

(14b)

The fluid flow and heat transfer governed by equations (1), (2) and (9), can be characterized by some dimensionless parameters: the Prandtl number Pr, the viscosity ratio J, the Rayleigh number Ra, the Grashof number Gr, the Darcy number Da and the Stefan number Ste, Pr 

vf

f

, J

c p ,f T gTL3 gTL3 e K , Ra  , Gr  , Da  2 , Ste  , 2 La  ff f f L

(15)

where ΔT represents the characteristic temperature difference and L denotes the characteristic length. In addition, conjugate thermal boundary condition should be employed at the interface between solid fins and porous mixture to ensure continuity of temperature and normal heat flux: Tfin  Tporous,

(16)

n  (kT ) fin  n  (kT ) porous ,

(17)

where n is normal to the interface. 3. MRT LB model for solid-liquid phase change in porous media In general, LB method for thermal flows can be grouped into three major categories: (1) the multispeed (MS) approach, (2) the double distribution function (DDF) approach and (3) the hybrid approach. Among three approaches, the DDF approach is commonly applied in the simulation of thermal flows due to its good numerical stability and simplicity. In this section, a total enthalpy-based MRT-LB model for melting in porous media is developed in the DDF approach framework, without iteration procedure for phase change treatment. 3.1 MRT LB equation for the flow field For the velocity field governed by

, the MRT-LB equation can be written as

S f i (r  eit , t  t )  f i (r , t )  (M 1SM)ij [ f j (r , t )  f jeq (r , t )]  t[M 1 (I  )M]ij F j , 2

(18)

where ei is the discrete velocity in direction i, fi (r, t) is the density distribution function with velocity ei at position r and time t, δt is the time increment, f i eq is the equilibrium distribution function and is defined as [35]

fi eq  i [1 

ei  u (ei  u) 2 u2   ], 2 4 cs 2cs 2cs2

(19)

where the weight coefficients ωi for the D2Q9 model can be given by ω0 = 4/9, ω1-4 = 1/9, and ω5-8 = 1/36.

cs is the sound speed. The discrete body force term Fi in Eq. (18) can be described as Fi  i [

ei  F (uF : ei ei ) u  F   ]. cs2 cs4 cs2

(20)

The macroscopic fluid density and velocity can be computed by [35]

   f i , u   ei f i /   i

i

t F. 2

(21)

The fluid velocity u is calculated using a temporal velocity V and is given by [35]

u

V d 0  d  d1 V 2 0

, V   ei f i /  

t

i

2

a ,

(22)

where the two parameters d0 and d1 are given by [35] d0 

t F 1 t  . (1   ) , d1   2 K 2 2 K

(23)

The transformation matrix M in eq. (13) for the D2Q9 model is given by 1 1 1 1 1 1  1   4 1 1 1 1 2 2  4 2 2 2 2 1 1   0 1 0 1 0 1 1  M 0 2 0 2 0 1 1  0 0 1 0 1 1 1  0 2 0 2 1 1  0  0 1 1 1 1 0 0  0 0 0 0 1 1  0

1 1   2 2 1 1   1 1   1 1  .  1  1   1  1 0 0  1  1

(24)

The diagonal relaxation matrix S is S  diag(s0 , s1 , s2 , s3 , s4 , s5 , s6 , s7 , s8 ) ,

(25)

where 0
In the past years, some approaches for treating conjugate heat transfer by using LB method have been proposed [37-38]. One of the common ways is to apply the correction of the incoming distribution function at interface nodes [37]. However, this kind of approach needs the identification of interface points and the angle between the lattice velocity components and the normal to the interface. An alternative way for treating conjugate conditions by introducing an additional source term in the evolution equation is developed to avoid any correction of distribution functions neighboring the interface [38]. However, in this approach, a relative complex source term requires to be calculated. In the LB method for solid-liquid phase change, the enthalpy-based LB method has been frequently employed due to its simplicity and computational efficiency [29-32, 39]. However, some of them are limited to treating phase change problem using iteration enthalpy method, which makes computation time-consuming. The efforts of simulating solid-liquid phase change heat transfer using LB method without iteration operations have been also conducted in some open literature [40-41]. Huang et al. [40] proposed the LB model based on the general energy equation involving the enthalpy term by combining the latent heat source term into the transient term of the conventional energy equation, avoiding iteration operations. In addition, it is worth pointing out that, for most of the existing LB models for solid-liquid phase change in porous media, there are the deviation terms in the recovered macroscopic equation from the LB equation, which are ignored with some assumptions. However, these deviation terms may still exist for most practical situations and affect the accuracy of the LB model [42-43]. In this paper, in order to eliminate this imperfectness, as well as inspired by the idea of Chopard et al. [43], an additional source term should be added to the LB evolution equation. Hence, the lattice Boltzmann equation of the modified SRT model is expressed as gi (r  eit , t  t )  gi (r , t )  

1

T

[ gi (r , t )  gieq (r , t )]  tSui ,

(26)

where g ieq is the equilibrium temperature distribution functions and is defined by g ieq

 hsys  T  i T  e u   T (  c p i 2 )  i cs 

i0 i0

,

(27)

where γ is a free parameter, which keeps unvaried over the entire space. Adding a free parameter, γ, into the equilibrium distribution function is to further enhance the numerical ability of the model [44-45].

g i and g ieq satisfy

g

i

e e g

  gieq  hsys ,  ei g ieq  c pTu ,

i i

eq i

 Tcs2  ,

(28)

where I is the unit tensor. The additional discrete source term Sui is written as

Sui  i (1 

1 2 T

)

e i ( c pTu)  , t cs2

(29)

and it satisfies

i Sui

 0,  ei Sui  (1  i

1 ( c pTu) ) . 2 T t

(30)

To recover the macroscopic equation (9) using equation (26), the following Chapman-Enskog analysis is performed. We first expand the distribution, the time and space derivatives and the discrete source as:

gi  gieq  gi(1)   2 gi( 2) , t  t1   2t 2 ,   1, Sui  Sui(1) ,

(31)

where ξ is a small expansion parameter. The Taylor series expansion for Eq. (26) to O(δt3) leads to Di gi 

t 2

Di gi  

1

2

 T t

( gi  g eq )  Sui , i

(32)

where Di   t  ei   . Based on this equation, we can derive the equations at different orders of parameter ξ as

 0 : gi(0)  gieq ,  1 : D1i gieq  

(33)

1

 T t

 2 :  t 2 gieq  D1i gi(1) 

gi(1)  Sui(1) ,

t 2

D12i gieq  

(34) 1

 T t

gi( 2) .

(35)

Substituting Eq. (34) into Eq. (35), one can rewrite Eq. (35) in another form

 2 :  t 2 gieq  (1 

1 t 1 ( 2) ) D1i gi(1)  D1i Sui(1)   gi . 2 T 2  T t

(36)

Summing Eqs. (34) and (36) over i and with the help of Eqs. (28) and (30), we can have

 t1 (hsys )  1  ( c pTu)  0 ,

 t 2 hsys  (1 

  1 t 1 )1  [ ei g i(1) ]  1  [1  ] t1 ( c pTu)  0 . 2 T 2 2  i T  

Using equations (34), (28) and (30), we can get

(37)

(38)

e g

(1) i i

  T t

 e (D i

i

eq 1i g i

 Sui(1) )

i

  T t[ t1 ( c pTu)  1  Tcs2 I  (1    T tcs21T 

t 2

1 2 T

) t1 ( c pTu)] ,

(39)

 t1 ( c pTu)

Substituting Eq. (39) into Eq. (38), we obtain

1    t 2 hsys  1  [cs2 ( T  )t ]1T  . 2  

(40)

Combining equations (37) and (40), and taking 1 k e  c s2 ( T  )t , 2

(41)

We can derive the recovered the macroscopic equation, i.e.,

 t hsys    ( c pTu)    [keT ] ,

(42)

which is the equation (9). Note that the above model is valid for the energy equation of solid fin. So far, the macroscopic variable, the total enthalpy hsys, can be obtained by hsys 

g

i

, but the

other three variables, namely, the total enthalpy hpcm, the temperature and the liquid fraction, are still unknown. As we all known, in the enthalpy methodology, if the total enthalpy hpcm is known, the temperature T and the liquid fraction γ can be determined by the thermodynamic relations as  Ts  (hs  hpcm ) /(  s c p ,s ) h  h h pcm  hs  l pcm T  Ts  Tl hl  hs  hl  hs   Tl  (hpcm  hl ) /(  f c p ,f ) 0   hpcm  hs fl    hl  hs  1

hpcm  hs hl  hpcm  hs ,

(43a)

hpcm  hl

hpcm  hs hl  hpcm  hs ,

(43b)

hpcm  hl

where hs  s c p ,sTs is the total enthalpy at the solidus temperature Ts, hl   f c p,f Tl  La is the total enthalpy at liquidus temperature Tl. Thus, the next key step is how to get the value of hpcm . Combining equations (8), (43a) and (43b), we have

h pcm

  ( c p ) s hsys    s c p , s  (1   )  m c p , m  hsys  (hl  hs )  (1   )  m c p , m  (hlTs  hsTl )   (hl  hs )    (1   )  m c p , m  (Tl  Ts )   f c p ,f (hsys   f La )    f c p , f  (1   )  m c p , m

hsys  hsys, s hsys, s  hsys  hsys, l , hsys  hsys, l

(44)

where hsys, s   s c p,sTs  (1   )  m c p, mTs is the corresponding total enthalpy of the mixture system at the PCM solidus temperature Ts. hsys, l   f c p,f Tl  f La  (1   )  m c p,mTl , is the total enthalpy of the mixture system at the PCM liquidus temperature Tl. Hence, the numerical steps for the temperature resolution are given below: (1) Perform streaming process, (2) Firstly calculate the total enthalpy hsys of the mixture system by hsys   g i , and then calculate the total enthalpy hpcm of PCM using equation (44), and at last compute the temperature T and the liquid fraction fl by equations (43a) and (43b), (3) Apply thermal boundary conditions, (4) Perform collision process. 3.3 Total enthalpy-based lattice Boltzmann model: MRT collision scheme The MRT LB equation for the total enthalpy hsys distribution function can be written as gi (r  eit , t  t )  gi (r , t )  (M 1ST )ij [n j (r , t )  neq j (r , t )]  tSui ,

(45)

where ni (r,t)=Mgi (r,t) is the rescaled moment with n ieq being the corresponding equilibrium moment that is written as n ieq  (hsys ,4hsys  2T ,4hsys  3T , c pT  c pT

ux , c

uy uy ux , c pT ,  c p T ,0,0) T c c c

(46)

where c is equal to x / t and x is the lattice spacing. ST is the diagonal relaxation matrix and can be written as ST  diag(0 , 1, 2 , 3 , 4 , 5 , 6 , 7 , 8 )

(47)

for the D2Q9 model, where 0  1 , 3  5  1 /  T , 0  1, 2, 4,6,7,8  2 . In order to reduce the numerical diffusion across phase change interface, some relaxation parameters can be further adjusted as [41] (

1

1,7,8

1 1 1 1  )(  ) 2 3,5 2 4

(48)

4. Results and discussions In this section, we first validate the proposed LB model through some numerical tests. To confirm the

applicability and accuracy of the LB model proposed above, we applied it to three different problems: two-region melting with different thermal properties between solid and liquid phases, melting with natural convection in a porous cavity and conjugate natural convection in porous media with a conducting wall. Then we apply the present LB model to simulate the melting processes in porous media with a conducting fin, and effect of fin on the melting heat transfer is investigated. In addition, it is noted that the choice of the free parameter γ is relatively arbitrary. As pointed out in Ref. [44], if the thermal conductivity ke is very small, γ can be tuned to keep the relaxation time be not close to 0.5. In Ref. [41], γ is recommended to be the harmonic mean of two heat capacities. Our numerical experiments show that the value of γ may not be too larger than twice of the minimum value of volumetric heat capacity between different phases, otherwise, it frequently leads to numerical instability. In the following simulations, we set it as the minimum value of volumetric heat capacity between different phases, which is of better numerical stability. In this work, the non-equilibrium extrapolation scheme with second-order accuracy [46] is adopted to the present model for velocity and thermal boundaries. The basic idea of this approach is to decompose the distribution function at the boundary node into its equilibrium and non-equilibrium parts. The non-equilibrium part is approximated by extrapolating from the non-equilibrium distribution at the neighboring node. 4.1. Two-region melting in a semi-infinite space The first test problem is the two-region melting problem in a semi-infinite space, as shown in Fig. 2. Initially, the PCM is kept in the solidus state with the temperature Tinit below the melting point of the PCM, Tm. At time t=0, the temperature at the boundary x=0 is suddenly increased to a high temperature Th, which exceeds the melting point of the PCM, and melting begins to occur. The analytical solutions for temperature distribution in the liquid and solid phases are given as [47]

Tf ( x, t )  Th erf ( x /(2  f t ))  Tm  Th erf ( )

(49a)

Ts ( x, t )  Tinit erfc( x /(2  st ))  Tm  Tinit erfc(  f /  s )

(49b)

where, the superscripts f and s represent the liquid and solid phases of PCM. The melting coefficient η is determined from the following transcendental equation [47]:

k  La  T  Ti e  e  ( f /  s )  s ( f )1 / 2 m  erf ( ) k f  s Tm  Th erfc(  f /  s ) c p ,f (Tm  Th ) 2

2

(50)

In the computations, the following parameters are given: Th=1, Tinit=0, Tm=0.5, Ste=cp,f (Th- Tm)=0.01. A grid size of 500×5 is employed and 100 lattices were selected as the reference length scale L. The dimensionless time, Fo, is defined as Fo=tαs/L2. Comparison between LB simulation results with the analytical solution for two cases is presented in Fig.3. It can be seen from Fig. 3 (a) that all the comparisons show an excellent agreement between the LBM results and the analytical solution for the case of kf /ks=5 and αf /αs=5. On the other hand, for the case of k f /ks =5 and αf /αs =25 as shown in Fig. 3(b), the present LBM result shows very good agreements with the analytical solution, while the previous LBM [28] result is mismatching with the analytical solution, which indicates that the LBM with the present formulation for melting heat transfer is effective for all cases, However, the previous model is only effective when the ratio of thermal conductivities between solid and fluid phases matches the ratio of thermal diffusivities between solid and fluid phases. This indicates that the present model can be used for the melting problem with different thermo-physical properties between two phases. 4.2. Melting with natural convection in a porous cavity In this subsection, we apply the present LB model to simulate the natural convection melting in a porous cavity investigated by Beckermann and Viskanta [13]. The physical model and boundary conditions are shown in Fig. 4, where the top and bottom walls are adiabatic, while the left and right walls are isothermal but held at different temperatures, Th and Tc, respectively. Initially, the domain are kept at temperature Tinit. At the time t=0, the left wall temperature is raised to Th, the PCM is allowed to melt. The main parameters are set as follow: ε=0.385, Ra=8.409×105, Da=1.37×10-5, Ste=0.1241, Pr=0.0208, J=1, Th=1.0, Tc=0, Tinit=0, Tm=0.3912, ke/kf =0.2719, kf/ks =1.0, σ=0.8604 (in melt) and σ=0.8352 (in solid). The equation (13a) is adopted as the empirical equation for the effective conductivities ke. In the computations, a lattice size is 160×160. Fig. 5 shows the comparison of the melting front locations predicted by LBM with experimental and numerical results in the ref. [13]. The present numerical results agree reasonably well with the experimentally measured and predicted results in [13], particularly for the case of the time Fo=7.314. Nevertheless, there are some discrepancies between the numerical and experimental results, which can be explained by the uncertainties in the measurements and effective

parameters as pointed out by Beckermann and Viskanta [13]. In addition, one can see from Fig. 5 that at the early time Fo=0.366, the melting front is almost parallel to the hot wall indicating a conduction-dominated melting. From the time Fo=1.829 to Fo=7.314, the melting front gradually presents a typical shape of convection-dominated melting. These features can also be observed from isotherms of melting in a porous cavity at the two times as shown in Fig. 6. One can also found from Fig. 6 that there is a large single eddy in the melt region for two cases. All of these observations are qualitatively consistent with those of numerical results in [13]. To quantify the results, we further compare the temperature profiles at different heights of the domain predicted by the present LBM with the results of previous studies in Fig. 7. It is seen that the present results agree relatively well with previous results. 4.3. Steady conjugate natural convection in porous media with a conducting wall Although natural convection in porous media has been extensively investigated by the LB method [35, 48-52], little work has been conducted to simulate conjugate natural convection in porous media. To further demonstrate the capacity of the proposed method in simulating conjugate heat transfer in porous media, in this subsection, we chose to solve steady conjugate natural convection in porous media with a conducting wall as depicted in Fig. 8. The left wall of thickness D is maintained at a hot temperature Th and the right wall is kept at a cold temperature Tc. The upper and lower walls of the cavity are assumed to be adiabatic. In the simulations, some parameters are set as: Th=1.0, Tc=0, ε=0.9, L/H=1, D/H=0.2, km/kf=100, (ρcp)m/(ρcp)f =1, (ρcp)w/(ρcp)f =5, Pr=1, Da=10-3, Gr=105 and the system size is chosen to be 120×120 grid size. In addition, the following convergent criterion was used to check on the convergence to reach the steady-state solution,

 T ( x, t )  T ( x, t  t )  10  T ( x, t )

9

,

(51)

where the summation is over the entire domain. Steady-state streamlines and isotherms at different wall-to-fluid thermal conductivity ratio are presented in Fig. 9, where the stream function is normalized by gTH 3 and the average Nusselt numbers is given as [53] x=-D:

Nu 

L H



H

0

k w  T  Tc ( ) x   D dy kf x Th  Tc

(52)

where kw and kf are the thermal conductivity of the solid wall and the fluid, respectively. As shown, the

predicted flow patterns by the present model agree well with the numerical results obtained by Al-Amiri et al [53]. It can be observed that there is a central vortex within the porous medium for all cases and its intensity intensifies as the wall-to-fluid thermal conductivity ratio increases. It is also found that the average Nusselt number increases with the increasing of wall-to-fluid thermal conductivity ratio. All of these observations are consistent with those of previous numerical studies in the Ref. [53]. For quantitative comparison, Table 1 presents comparisons of the average Nuesslt number and the minimum dimensionless stream function at steady state predicted by the present model with results obtained by Al-Amiri et al [53]. It can be found that the predicted results by the present LB model are in excellent agreement with those reported data. 4.4. Melting in a porous cavity with a conducting fin In this subsection, we apply the present LB model to investigate the melting in a porous cavity with a conducting fin attached on the left wall as shown in Fig. 1. The effects of the fin-PCM heat capacity ratio, as well as the length and vertical position of the fin on the melting processes are investigated. In the following simulations, the required parameters are fixed as follows: Th=1.0, Tc=0, Tinit=0, Tm=0, ε=0.9, Pr=1.0, Ra=105, Da=10-3, ks/kf=1.0, km/kf=300, kfin/kf=500, (ρcp)s/(ρcp)f=1.0, (ρcp)m/(ρcp)f=1.2, L/H=1.0 and wf/H=0.05. The effective conductivities ke of porous media with high porosity is calculated by Eq. (13c) and the geometric function Fε of porous media is determined as Eq. (14b). A grid-independence study was firstly carried out to test the sensitivity of the chosen grids in the present simulations for the case without fin. Different grid sizes were selected and tested to ensure the independence of solution from the adopted grid size based on the comparison of melting volume fraction. Three grids-80×80 (grid #1), 120×120 (grid #2) and 160×160 (grid #3) were tested. For grid #1, a deviation of 1.81% was observed in the melt volume fraction with respect to grid #3 at time Fo=0.001, with this deviation between grid #2 and #3 reducing to 0.45%. Therefore grid #2 (120×120) was chosen for the following calculations considering both accuracy and computational cost. Fig. 10 shows temperature contours and streamlines for melting in a porous cavity with no fin and with fin at the dimensionless time Fo=tαf/L2=0.00255. It can be observed that

a low heat capacity heats more quickly than because in heat transfer, a lower

a high

heat capacity. This is

heat capacity means a higher thermal diffusivity, which

measures the rate of transfer of heat of a material from the hot side to the cold side. For further understanding of the effect of adding the fin on the melting processes in porous media, Fig. 11 compares the melt volume fraction of the molten PCM as a function of time for the cases without fin and with fin. Here, a volume fraction of the melted PCMs is defined as: fl  

Vf , Vtotal

(53)

where Vf and Vtotal denote the volume of the molten PCM and the total domain, respectively. As shown, the melt speed for the case with fin is higher than that for the case without fin, and a higher melt speed is also induced by adding a fin with a lower

heat capacity.

The length of the fin is generally considered as an important design parameter and thus the effect of the length of the fin on the melting processes in the porous cavity is investigated in the following part. Fig. 12 illustrates the melt volume fraction versus dimensionless time for hf/H=0.5 and different lengths of the fin. It can be seen that the melting speed of PCM is accelerated by increasing fin length from 0.2H to 0.8 H. For instance, at the dimensionless time of Fo=0.01, when the fin length is changed to 0.2H, 0.5H and 0.8H, the melt volume fraction increases 5.0%, 15.0% and 22.3% in comparison with the case without fin. The PCM region is fully melted at about Fo=0.02 with no fin, Fo=0.018 with lf =0.2H, Fo=0.015 with lf =0.5H and Fo=0.014 with lf =0.8H. The effect of vertical position of the fin on the melting processes is also investigated as shown in Fig. 13. It can be observed that varying the vertical position from the bottom to the top of the cavity has no remarkable impact on the melting speed. In particular, the melt history of the case with hf/H==0.25 is almost the same as that of the case with hf/H==0.75. This may be due to the fact that the melting process is dominated by conduction heat transfer for the cases considered. 5. Conclusions In this paper, a MRT-LB model has been developed for simulating transient melting heat transfer of phase change material in porous media coupled with conduction fins at the representative element volume scale. In this model, a total enthalpy-based LB equation is constructed to solve the temperature field in

addition to the MRT-LB equation with the density distribution function for the flow field. The general macroscopic energy equation can be correctly recovered from the proposed LB formulation by using an equilibrium temperature distribution function and an additional source term in the evolution of temperature field. Numerical results for two-region melting with different thermal properties between solid and liquid phases, melting in a porous cavity and conjugate natural convection in porous media with a solid wall, agree well with the analytical or numerical solutions in previous studies. Besides, the melting of PCMs in porous media with a conducting fin was numerically investigated. The numerical results indicate that the conducting fin added into the porous medium can further improve the melting heat transfer. Furthermore, the melting speed increases when the length of the fin becomes greater and the heat capacity of the fin becomes lower. However, it is found that varying the vertical position of the fin has no remarkable impact on the melting speed for the cases considered. The present model, as a generalized one, can be directly applied to simulate conjugate heat transfer in heterogeneous media and solidification processes in porous media. The present model is based on local thermal equilibrium assumption, which may not be valid in some practical application. Thus, extensions of the present model for modeling conjugate heat transfer with phase change in porous media under local thermal non-equilibrium condition will be conducted in our future studies.

Acknowledgement The authors thank Dr. Zhang Chengbin at the Southeast University of China for helpful discussions. This work was supported by the National Natural Science Foundation of China (No. 51206076, No. 51408302 and No. 51676037), and the Science Fund of Nanjing Institute of Technology of China (No. CKJA201604 and No. CKJC201503 ).

Nomenclature c

lattice speed, (m/s)

cp

specific heat at constant pressure (kJ/(kg K))

cs

sound speed (m/s)

Da

Darcy number, Da=K/L2

ei

discrete lattice velocity in direction i (m/s)

fi

density distribution function in direction i (kg/m3)

fi eq

equilibrium distribution function of density in direction i (kg/m3)

fl

liquid volume fraction in pore space

fl *

liquid volume fraction based on the entire domain

F

body force per unit mass (N/kg)

Fi

discrete body force in direction i (kg/( m3 s))

Fo

Fourier number, Fo=αt/L2

Fε

Forchheimer form coefficient

g

acceleration due to gravity (m/s2)

G

external body force per unit mass (N/kg)

gi

temperature distribution function in direction i (K)

gi eq

equilibrium temperature distribution function in direction i (K)

Gr

Grashof number

h

enthalpy (J/m3)

hf

vertical position of the fin (m)

H

height of the domain or characteristic length (m)

I

unit tensor

J

viscosity ratio, J=νe/νf

k

thermal conductivity (W/(m K))

K

permeability (m2)

lf

length of the fin (m)

L

Weight of the domain or characteristic length (m)

M

transformation matrix

n

moment of the distribution function

Nu

Nusselt number

p

pressure (Pa)

Pr

Prandtl number, Pr= νf/αf

r

space position (m)

Ra

Rayleigh number, Ra=|g|βΔTL3/(νfαf)

S, ST

diagonal relaxation matrix

Ste

Stefan number

Sui

discrete source term related to velocity (W/m3)

t

time (s)

T

temperature (K)

Tc

temperature of cold wall (K)

Th

temperature of hot wall (K)

Tinit

initial temperature (K)

Tm

melting temperature (K)

u

velocity (m/s)

V

temporal velocity (m/s)

wf

width of the fin

x, y

Cartesian coordinates (m)

Greek Symbols α

thermal diffusivity (m2/s)

β

coefficient of thermal expansion (1/K)

γ

a constant parameter

δx

lattice space (m)

δt

time step (s)

ε

porosity of porous media

νe

effective kinematic viscosity (m2/s)

ξ

small expansion parameter

ρ

density (kg/m3)

σ

ratio of volumetric heat capacities of mixture and fluid

τ, τT

dimensionless relaxation time

Ψ

stream function

ωi

weight number in direction i

Subscripts e

effective or equivalent

f

fluid

fin

fin

i

direction i in a lattice

l

liquidus

m

porous matrix

PCM

phase change material

s

solid matrix or solidus

sys

porous matrix/solid/liquid mixture system

w

wall

References [1] B. Zaiba, J.M. Marin, L.F. Cabeza, H. Mehling, Review on thermal energy storage with phase change: materials, heat transfer analysis and applications, Appl. Therm. Eng. 23 (2003) 251-283. [2] P. Zhang, X. Xiao, Z.W. Ma, A review of the composite phase change materials: Fabrication, characterization, mathematical modeling and application to performance enhancement, Appl. Energ. 165 (2016) 472-510. [3] Z. Rao, S. Wang, A review of power battery thermal energy management, Renew. Sust. Energ. Rev. 15 (2011) 4554-4571. [4] S. Hong, D. R. Herling, Open-cell aluminum foams filled with phase change materials as compact heat sinks, Scr. Mater. 55 (2006) 887-890. [5] C.Y. Zhao, W. Lu, Y. Tian, Heat transfer enhancement for thermal energy storage using metal foams embedded within phase change materials (PCMs), Sol. Energ. 84 (2010) 1402-1412. [6] X. Xiao, P. Zhang, M. Li, Effective thermal conductivity of open-cell metal foams impregnated with pure paraffin for latent heat storage, Int. J. Therm. Sci. 81 (2014) 94-105.

[7] S. Feng, Y. Zhang, M. Shi, T. Wen, T.J. Lu, Unidirectional freezing of phase change materials saturated in open-cell metal foams, Appl. Therm. Eng. 88 (2015) 315-321. [8] X. Xiao, P. Zhang, Morphologies and thermal characterization of paraffin/carbon foam composite phase change material, Sol. Energy Mater. Sol. Cells 117 (2013) 451-461. [9] Y.W. Zhang, A. Faghri, Heat transfer enhancement in latent heat thermal energy storage system by using the internally finned tube, Int. J. Heat Mass Transfer 39 (1996) 3165-3173. [10] A. Castell, C. Sole, M. Medrano, J. Roca, L. Cabeza, D. Garcia, Natural convection heat transfer coefficients in Phase Change Material (PCM) modules with external vertical fins, Appl. Therm. Eng. 28 (2008) 1676-1686. [11] P. Lamberg, K. Siren, Analytical model for melting in a semi-infinite PCM storage with an internal fin, Heat Mass Transfer 39 (2003) 167-176. [12] S. Feng, M. Shi, Y. Li, T. J. Lu, Pore-scale and volume-averaged numerical simulations of melting phase change heat transfer in finned metal foam, Int. J. Heat Mass Transfer 90 (2015) 838-847. [13] C. Beckermann, R. Viskanta, Natural convection solid/liquid phase change in porous media, Int. J. Heat Mass Transfer 31 (1988) 35-46. [14] O. Mesalhy, K. Lafdi, A. Elgafy, K. Bowman, Numerical study for enhancing the thermal conductivity of phase change material (PCM) storage using high thermal conductivity porous matrix, Energy Convers. Manage. 46 (2005) 847-867. [15] S. Krishnan, J.Y. Murthy, S.V. Garimella, A two-temperature model for solid–liquid phase change in metal foams, J. Heat Transfer 127 (2005) 995-1004. [16] Y. Tian, C.Y. Zhao, A numerical investigation of heat transfer in phase change materials (PCMs) embedded in porous metals, Energ. 36 (2011) 5539-5546. [17] W.Q. Li, Z.G. Qu, Y.L. He, W.Q. Tao, Experimental and numerical studies on melting phase change heat transfer in open-cell metallic foams filled with paraffin, Appl. Therm. Eng. 37 (2012) 1-9. [18] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222 (1992) 145-197. [19] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329-364.

[20] X. He, S. Chen, G.D. Doolen, A novel thermal model for the lattice Boltzmann method in incompressible limit, J. Comput. Phys. 146 (1) (1998) 282-300. [21] S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond, Clarendon Press, Oxford, 2001. [22] A.A. Mohamad, Lattice boltzmann method, fundamentals and engineering applications with computer codes, Springer, 2011. [23] Q. Li, Y.L. He, Y. Wang, W.Q. Tao, Coupled double-distribution-function lattice Boltzmann method for the compressible Navier–Stokes equations, Phys. Rev. E 76 (2007) 056705. [24] Z. Chai, C. Huang, B. Shi, Z. Guo, A comparative study on the lattice Boltzmann models for predicting effective diffusivity of porous media, Int. J. Heat Mass Transfer 98 (2016) 687-696. [25] D. Sun, M. Zhu, J. Wang, B. Sun, Lattice Boltzmann modeling of bubble formation and dendritic growth in solidification of binary alloys, Int. J. Heat Mass Transfer 94 (2016) 139-149. [26] A. Xu, T.S Zhao, L. Shi, X.H Yan, Three-dimensional lattice Boltzmann simulation of suspensions containing both micro- and nanoparticles, Int. J. Heat Fluid Flow 62 (2016) 560-567. [27] Z. Chen, D. Gao, J. Shi, Experimental and numerical study on melting of phase change materials in metal foams at pore scale, Int. J. Heat Mass Transfer 72 (2014) 646-655. [28] D. Gao, Z. Chen, Lattice Boltzmann simulation of natural convection dominated melting in a rectangular cavity filled with porous media, Int. J. Therm. Sci. 50 (2011) 493-501. [29] D. Gao, Z. Chen, L. Chen, A thermal lattice Boltzmann model for natural convection in porous media under local thermal non-equilibrium conditions, Int. J. Heat Mass Transfer 70 (2014) 979-989. [30] D. Gao, Z. Chen, M. Shi, Z. Wu, Study on the melting process of phase change materials in metal foams using lattice Boltzmann method, Sci. China Tech. Sci. 53 (2010) 3079-3087. [31] Q. Liu, Y.L. He, Double multiple-relaxation-time lattice Boltzmann model for solid-liquid phase change with natural convection in porous media, Phys. A 438 (2015) 94-106. [32] Y.B. Tao, Y. You, Y.L. He, Lattice Boltzmann simulation on phase change heat transfer in metal foams/paraffin composite phase change material, Appl. Therm. Eng. 93 (2016) 476-485. [33] W. Wu, S. Zhang, S. Wang, A novel lattice Boltzmann model for the solid-liquid phase change with the convection heat transfer in the porous media, Int. J. Heat Mass Transfer 104 (2017) 675-687.

[34] K. Bhattacharya, V.V. Calmidi, R.L. Mahajan, Thermo-physical properties of high porosity metal foams, Int. J. Heat Mass Transfer 45 (2002) 1016-1031. [35] Z. Guo, T.S. Zhao, A lattice Boltzmann model for convection heat transfer in porous media, Numer. Heat Transfer, Part B 47 (2005) 157-177. [36] C. Huber, A. Parmigiani, B. Chopard, M. Manga, O. Bachmann, Lattice Boltzmann model for melting with natural convection, Int. J. Heat Fluid Flow 29 (2008) 1469-1480. [37] L. Li, C. Chen, R. M, J.F. Klausner, Conjugate heat and mass transfer in the lattice Boltzmann equation method, Phys. Rev. E 89 (2014) 043308. [38] H. Karani, C. Huber, Lattice Boltzmann formulation for conjugate heat transfer in heterogeneous media, Phys. Rev. E 91 (2015) 023304. [39] W.S. Jiaung, J.R. Ho, C.P. Kuo, Lattice-Boltzmann method for the heat conduction problem with phase change, Numer. Heat Transfer: Part B 39 (2001)167-187. [40] R. Huang, H. Wu, P. Cheng, A new lattice Boltzmann model for solid-liquid phase change, Int. J. Heat Mass Transfer 59 (2013) 295-301. [41] R. Huang, H. Wu, Phase interface effects in the total enthalpy-based lattice Boltzmann model for solid–liquid phase change, J. Comput. Phys. 294 (2015) 346-362. [42] B. Chopard, J.L. Falcone, J. Latt, Lattice Boltzmann advection-diffusion model revisited, Euro. Phys. J. Special Topics 171 (2009) 245-249. [43] Z. Chai, T.S. Zhao, Lattice Boltzmann model for the convection-diffusion equation, Phys. Rev. E 87 (2013) 063309. [44] Z. Chai, B. Shi, Z. Guo, A multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropic convection-diffusion equations, J. Sci. Comput. 69 (2016) 355-390. [45] D. Gao, Z. Chen, L. Chen, D. Zhang, A modified lattice Boltzmann model for conjugate heat transfer in porous media, Int. J. Heat Mass Transfer 105 (2017) 673-683. [46] Z. Guo, C. Zheng, B. Shi, An extrapolation method for boundary conditions in lattice Boltzmann method, Phys. Fluids 14 (2002) 2006-2010. [47] M.N. Ozisik, Heat Conduction, 2nd edition, John Wiley & Sons Press, New York, 1993. [48] F. Rong, Z. Guo, Z. Chai, B. Shi, A lattice Boltzmann model for axisymmetric thermal flows through

porous media, Int. J. Heat Mass Transfer, 53 (2010) 519-5527. [49] Q. Liu, Y.L. He, Q. Li, W.Q. Tao, A multiple-relaxation-time lattice Boltzmann model for convection heat transfer in porous media, Int. J. Heat Mass Transfer, 73 (2014) 761-775. [50] L. Wang, Z. Zeng, L. Zhang, H. Xie, G. Liang, Y. Lu, A lattice Boltzmann model for thermal flows through porous media, Appl. Therm. Eng., 108 (2016) 66-75. [51] L. Wang, J. Mi, Z. Guo, A modified lattice Bhatnagar–Gross–Krook model for convection heat transfer in porous media, Int. J. Heat Mass Transfer, 94(2016) 269-291. [52] Y. Hu, D. Li, S. Shu, X. Niu, Immersed boundary-lattice Boltzmann simulation of natural convection in a square enclosure with a cylinder covered by porous layer, Int. J. Heat Mass Transfer 92 (2016) 1166-1170. [53] A. Al-Amiri, K. Khanafer b, I. Pop, Steady-state conjugate natural convection in a fluid-saturated porous cavity, Int. J. Heat Mass Transfer 51 (2008) 4260-4275.

Captions of figures Fig. 1. Schematic diagram for melting in porous medium with a conducting fin Fig. 2. Schematic of two-region melting problem Fig. 3. Comparisons of LBM solutions with the analytical solution of two-region melting for (a) kf /ks=5 and αf /αs=5, and (b) kf /ks=5 and αf /αs=25 Fig. 4. Schematic of melting with natural convection in a porous cavity heated from the side Fig. 5. Comparison of experimentally measured and predicted interface locations at various times Fig. 6. Isotherms and streamlines of melting in a porous cavity at two times (a) Fo=1.829 and (b) Fo=7.314 Fig. 7. Comparisons of experimentally measured and predicted temperature profiles at different height of the domain at two times (a) Fo=1.829 and (b) Fo=7.314 Fig. 8. Schematic diagram of conjugate heat transfer in porous medium with a conducting wall Fig. 9. Steady-state isotherms and streamlines at different wall-to-fluid thermal conductivity ratio Fig. 10. Temperature contours and streamlines for melting in a porous cavity without fin and with a fin at the time Fo=0.00255 Fig. 11. Comparison of melt volume fraction versus time for the cases without fin and with fin Fig. 12. Melt volume fraction versus time for different lengths of the fin Fig. 13. Melt volume fraction versus time for different vertical positions of the fin

Captions of tables Table 1 Comparisons of the average Nusselt number and the minimum stream function predicted by the present model with numerical results in Ref. [53] at steady state

y Adiabatic

ux=uy=0

H PCM + Porous media

Th

Tc

ux=uy=0

Fin

ux=uy=0

wf

lf

hf

g 0

Adiabatic

ux=uy=0

L

x

Fig. 1. Schematic diagram for melting in porous medium with a conducting fin

T Th

Liquid PCM Tl(x,t)

Solid PCM Ts(x,t)

Tm Tinit

Tinit

Interface Fig. 2 Schematic of two-region melting problem

x

1.00 Analytical Previous model [28] Present model

0.75

0.50

T

kf / ks=5, f /s=5

0.25 Ste=0.01, Fo=0.2 0.00 0.0

0.2

0.4

0.6

0.8

1.0

x/L (a) kf /ks=5 and αf /αs=5

1.00 Analytical Previous model [28] Present model

T

0.75

0.50

0.25

kf / ks=5, f /s=25 Ste=0.01, Fo=0.2

0.00 0.0

0.2

0.4

0.6

0.8

1.0

x/L (b) kf /ks=5 and αf /αs=25 Fig. 3 Comparisons of LBM solutions with the analytical solution of two-region melting for (a) kf /ks=5 and αf /αs=5, and (b) kf /ks=5 and αf /αs=25

ux=uy=0

Adiabatic L Interface

Th

Liquid PCM + Porous matrix

Solid PCM + Porous + matrix

Tc

ux=uy=0

ux=uy=0 Tinit

y

Tm 0

x Adiabatic

ux=uy=0

L

Fig. 4 Schematic of melting with natural convection in a porous cavity heated from the side

1.0 0.8

Present

0.6

y/H

Numer. [13]

0.4

Exp. [13] Fo 0.366 1.829 7.314

0.2 0.0 0.0

0.2

0.4

x/L

0.6

0.8

1.0

Fig. 5 Comparison of experimentally measured and predicted interface locations at various times

Frame 001  10 Sep 2016  field

Frame 001  09 Sep 2016  field

0.8

0.8

0.6

0.1 Y 0.05

0.2

Y

0.6

0.3

1

0.9 0.8 0.7 0.6 0.5 0.4

1

0.4

0.4

0.2

0.2 Frame 001  07 Sep 2016  field

Frame 001  07 Sep 2016  field

0

0

0

0.2

0.4

0.6

(a) 0 1.2 0.2 1 Fo=1.829

0.8

0.6

0.8

1

1.2

1

0.8

0.6 0. 5

0.7

0.8

0.8

0.9

1

0.6

Y 0.2

0.4

0.4

0.1

0.3

Y

0.4

0.6

0.2

0.2

0

0

-0.2

0.4

X

X

0

(b) Fo=7.314 -0.2 0.2 0.4 0.6 0.8 1 0 1.2 0.2 0.4 0.6 0.8 1 Fig. 6 Isotherms and streamlines of melting in a porous cavity at two times (a) Fo=1.829 and (b) X

X

Fo=7.314

1.2

1.0

y/H Exp. [13] Numer. [13] 0.1333 0.5000 0.8667 y/H Present 0.1333 0.5000 0.8667

0.8

T

0.6 0.4 0.2 0.0 0.0

Fo=1.829 0.2

0.4

0.6

0.8

1.0

x/L (a) Fo=1.829

1.0 Fo=7.314 0.8

y/H Present 0.1333 0.5000 0.8667

T

0.6 0.4 y/H Exp. [13] Numer. [13] 0.2 0.1333 0.5000 0.0 0.8667 0.0 0.2 0.4 0.6

0.8

1.0

x/L (b) Fo=7.314 Fig. 7 Comparisons of experimentally measured and predicted temperature profiles at different height of the domain at two times (a) Fo=1.829 and (b) Fo=7.314

y H Adiabatic ux=uy=0

ux=uy=0

ux=uy=0 Porous medium

Th

Tc

才c

solid wall

g

D

0

Adiabatic

ux=uy=0

L

x

Fig. 8 Schematic diagram of conjugate heat transfer in porous medium with a conducting wall

Frame 001  12 Sep 2016  field

Frame 001  12 Sep 2016  field

1.2

1.2

1

1 -6.12E-05 -1.64E-04 -2.66E-04 -3.68E-04

1.000 0.8 0.862

0.8 0.6

0.724 0.586 0.447 0.4 0.309 Frame 001  12 Sep 2016  field 0.171 0.2 0.033

0.6

Y

Y

0.033

0.4 0.2

Frame 001  12 Sep 2016  field

0

-8.80E-04

0 1.2 (a) kw/kf =0.1 (Nu=0.489) -0.2

1.2

-0.2 -0.4

-4.71E-04 -5.73E-04 -6.75E-04 -7.78E-04

1

-0.4

-0.2

0

0.8

0.2

0.4

0.6

X

0.8

1 -0.4 1 1.2 -0.41.4-0.2 1.000 0.8

Y

Y

0.6

Frame 001 0.2 12 Sep 2016  field

0.2

0.4

0.6

0.8

-4.06E-04 1.2 1.4 -1.12E-03 -1.83E-03 -2.54E-03

1

X

0.862

0.6 0.4

0

0.724

0.586 0.447 0.40.309 0.171

-3.25E-03 -3.96E-03 -4.68E-03 -5.39E-03 -6.10E-03

Frame 001  12 Sep 2016  field 0.2

0.033

0

0 1.2 (b) -0.2kw/kf =1.0 (Nu=4.488)

1.2 -0.2 1

0

0.5

X

1 -0.4

1

-0.4 1.000 0.8 0.862

0.8

-0.2

0

0.2

0.4

0.6

0.8

1

1.2-1.00E-03 1.4 -2.71E-03 -4.43E-03 -6.14E-03 -7.86E-03 -9.57E-03 -1.13E-02 -1.30E-02

X

0.724 0.6 0.586 0.447 0.4 0.309 0.171 Frame 0010.2  12 Sep 2016  field 0.033

Y

Y

0.6 0.4 Frame 0010.2  12 Sep 2016  field

0

0 1.2 -0.2

1.2 (c)-0.2 kw/kf =5 (Nu=7.745)

1

1

0

0.5

1

1.000 0.8 0.862

X

0.8

Y

Y

0.6 0.4 0.2 0

0

0.5

1

X

-1.00E-03 -3.00E-03 -5.00E-03

0.724 0.6 0.586 0.447 0.4 0.309 0.171 0.2 0.033

-7.00E-03 -9.00E-03 -1.10E-02 -1.30E-02 -1.50E-02

0

(d) -0.2 kw/kf =10 (Nu=9.172)

-0.2

Fig. 9 Steady-state isotherms and streamlines at different wall-to-fluid thermal conductivity ratio 0

0.5

X

1

0

0.5

X

1

Frame 001  09 Oct 2016  field Frame 001  09 Oct 2016  field

1.2 1

1

0.95 0.85 0.8 0.75 0.65 0.6 0.55 0.45

0.8

Y

0.6

0.4

0.4

0.35 0.25 Frame 001 0.2  09 Oct 2016  field 0.15 0.05

0.2 Frame 001  09 Oct 2016  field

0

0 1.2

(a) without fin 0

0.2

0.4

0.6

0.8

X

1

1 -0.2 1.2

1

0.8

0.6

0.6

0.4

0.4

Y

0.8

0.95 -0.2 0.85 0.75 0.65 0.55 0.45

0

0.2

0.4

0.6

0.8

1

1.2

0.4

0.6

0.8

1

1.2

0.35 0.25 Frame 0010.2  09 Oct 2016  field 0.15 0.05

0.2 Frame 001  09 Oct 2016  field

0 1.2

0

(b) (ρcp)fin/(ρcp)f =20, hf/H=0.5, lf/H=0.8 1

0

0.2

0.4

0.6

0.8

X

1-0.2

1

0.6

0.6

0.4

0.4

0.2

0.2

0

0

Y

0.8

0.8

1.2

0.95 0.85 0.75 0.65 0.55 0.45

-0.2

0

0.2

0.35 0.25 0.15 0.05

(c) (ρcp)fin/(ρcp)f =1.05, hf/H=0.5, lf/H=0.8 0

0.2

0.4

0.6

0.8

-0.2 1

1.2

X and streamlines for melting in a porous cavity without fin and with a fin at Fig. 10. Temperature contours 0 the time-0.2 Fo=0.00255

0.2

0.4

0.6

0.8

1

1.2

1.0 0.8 0.6 fl *

Without fin (cp)fin/(cp)f=20

0.4 0.2 0.0 0.000

(cp)fin/(cp)f=1.05 =0.9, Pr=1.0, Ra=10 , Da=10 , kfin/kf=500, lf/H=0.8, hf/H=0.5 5

0.005

0.010 Fo

-3

0.015

0.020

Fig. 11. Comparison of melt volume fraction versus time for the cases without fin and with fin

1.0 0.8 without fin lf/H=0.2 lf/H=0.5 lf/H=0.8

fl *

0.6 0.4 0.2 0.0 0.000

=0.9, Pr=1.0, Ra=10 , Da=10 , kfin/kf=500, (cp)fin/(cp)f=1.05, hf/H=0.5 5

0.005

0.010 Fo

-3

0.015

0.020

Fig. 12. Melt volume fraction versus time for different lengths of the fin

1.0 0.8 hf/H=0.25 hf/H=0.5 hf/H=0.75

fl *

0.6 0.4

=0.9, Pr=1.0, Ra=10 , Da=10 , kfin/kf=500, (cp)fin/(cp)f=1.05, lf/H=0.8 5

0.2 0.0 0.000

0.005

0.010

-3

0.015

Fo Fig. 13. Melt volume fraction versus time for different vertical positions of the fin

Table 1 Comparisons of the average Nusselt number and the minimum stream function predicted by the present model with numerical results in Ref. [53] at steady state Ψmin

Nu kw/kf Present

Ref. [53]

% Error

Present

Ref. [53]

% Error

0.1

0.489

0.478

2.30

-8.8755E-04

-8.8079E-04

0.76

1

3.488

3.433

1.60

-6.1661E-03

-6.1681E-03

0.03

5

7.745

7.75

0.06

-1.3084E-02

-1.3209E-02

0.95

10

9.172

9.168

0.04

-1.5338E-02

-1.5444E-02

0.69

Research highlights 

A modified model is proposed for melting in porous media with fins.



The model it suitable for modeling conjugate heat transfer in porous media.



The feasibility of the model are verified by several numerical tests.



Adding a conducting fin into porous medium can further enhance heat transfer.