Enthalpic lattice Boltzmann formulation for heat conduction during melting of PCMs with embedded solid blocks with different thermophysical properties

Enthalpic lattice Boltzmann formulation for heat conduction during melting of PCMs with embedded solid blocks with different thermophysical properties

Rihab Hamila*Abdelmajid Jemni | Sassi Ben Nasrallah | Patrick Perré

Laboratoire d’Études des Systèmes Thermiques et Énergétiques, École Nationale d’Ingénieurs de Monastir, Université de Monastir, Avenue Ibn El Jazzar, 5019 Monastir, Tunisia

Laboratoire de Génie des Procédés et Matériaux, Centrale Supélec, Université de Paris-Saclay Grande Voie des Vignes 92295 Châtenay-Malabry Cedex, France

Corresponding Author Email: 
hamilara7@gmail.com
Page: 
330-338
|
DOI: 
https://doi.org/10.18280/ijht.350214
Received: 
| |
Accepted: 
| | Citation

OPEN ACCESS

Abstract: 

In this paper, we extend the enthalpic formulation introduced in previous work of Hamila et al. [1], to solve transient conduction heat transfer during melting of phase change material. The volume fraction of liquid phase is described by a Heaviside step function smoothly approximated by a logistic curve having a common"S" shape, also called sigmoid curve. Several conduction test problems during melting of some specific phase change material are simulated and comparisons are made with the control volume solutions. The simulation of conduction heat transfer during melting of commercial PCM filling the cavity, in which one or more circular solid blocks (with enhanced thermophysical properties) are embedded, is performed. The proposed formulation can handle accurately conduction heat transfer problem with phase change including interface evolving over time. Special attention is given to distinction between thermophysical properties of solid and liquid phase which are not equal.

Keywords: 

phase change material, lattice boltzmann method, diffusion, melting

1. Introduction

Solid-liquid phase change is a basic physical phenomenon that has wide applications in engineering, such as the casting of metals, melting and solidification of alloys, artificial crystal growth from melts, green building, thermal energy storage, etc. When dealing with phase change problem, the main issue that has to be worked on is the propagation of the moving phase change interface, causing inherent nonlinearities. During the phase change a large amount of latent heat is released or absorbed on the phase interface, which gives rise to a significant and sharp variation into thermophysical properties near the phase change interface. That's why accurate simulation of solid-liquid phase change is a challenging task.

The lattice Boltzmann method (LBM) is a widely used numerical technique which has been successfully applied to various hydrodynamics problems [2, 3], fluid flows and complex physical processes in fluids [4-7]. Its main asset is  its simple algebraic manipulation, its easy solution procedure and implementation of boundary conditions, together with its ability of dealing with complex fluids [8]. This explains why LBM is a very promising and competitive numerical tool in solving heat transfer process [9-11], turbulence flows [12, 13],micro-flows [14, 15], porous media [16-18], multiphase flow [19, 20] and simulation of solid-liquid phase change  [21-23].Among authors who have used the LBM to solve phase change heat conduction problems, Jiuang et al. [22] have extended the lattice Boltzmann method to solve a phase change problem governed by the heat conduction equation using an enthalpic formulation. The proposed formulation is based on bringing out an extra source term added to the lattice Boltzmann equation, in which calculation of the derivative of the liquid phase volume phase fraction is performed. To validate the feasibility and accuracy of the numerical scheme presented, some test problems have been simulated, namely melting in a half-space, solidification in a half-space, benchmark three region problem with solidification phenomena, and finally solidification from a corner in a quarter spaces. When solving solidification in a half space, they set both thermal diffusivity and relaxation time to be unity, next when they simulated three regions problem, they specified thermal conductivities and heat capacities ratio between solid phase, mushy zone and liquid phase. Their validations have been done with analytical solutions. Huang et al. [23] used an immersed boundary thermal lattice Boltzmann method to simulate solid-liquid phase change problems. Tracking the liquid-solid interface was made through a Lagrangian grid making this interface represented as a sharp interface rather than a diffusive one. Numerical examples, including conduction-induced melting in a semi-infinite space and melting in a square cavity, were carried out and good agreements were found. Beside the complexity of the proposed LB formulation, the authors assumed that the thermophysical properties of both liquid and solid phases are equal and constant, limiting their results to a specific case study. In their work, Li et al. [24] simulated two melting problems by a lattice Boltzmann method with an interfacial tracking method. Both conduction- and convection-controlled melting problems were solved in order to validate the proposed formulation. Two different lattice schemes were used to compute thermal and velocity fields. The authors assumed that the thermophysical properties of both phases are equal and constants. Kang et al. [25] opted for a double-population lattice Boltzmann method to solve a convection-diffusion problem involving solid-liquid phase transition process. To obtain the melt fraction, the authors linearly interpolate the enthalpy and subdivide the expression into three domains describing the solid phase, the liquid phase and the melting front. Validation of their model was performed by simulating first a square cavity initially filled with a solid pure substance and melting with convection transfer. Next they extended their work by demonstrating the ability of LBM to simulate complex systems by solving a Shell with different numbers of tubes. They also investigated the effect of the number and arrangement of inner tubes on multi-tube thermal energy storage systems. However during all their simulations they made assumptions stating that densities of both phases are equal, specific heats of both solid and liquid phases are constant and equal and both thermal conductivities are constant and equal. Huang et al. [26] developed a formulation of lattice Boltzmann method able to solve solid- liquid phase change problems, by modifying the equilibrium distribution function for the temperature. The phase interface was traced by updating the total enthalpy, and the moving interface was treated by the immersed boundary scheme. The authors assumed that thermophysical properties of both liquid and solid phases of pure substance are equal and constant. Huang et al. [27] adapted a total enthalpy-based lattice Boltzmann method with a mesh refinement process to efficiently simulate solid-liquid phase change problem. Using an MRT lattice Boltzmann formulation to compute the thermal and fluid flow, the authors simulate the thermal field by a modified relaxation matrix. When it comes to setting the thermophysical properties of solid and liquid phase, even if the authors claim that the density and heat capacity of the solid and the liquid are different, they assumed during their simulations, equal densities and heat capacities of both phases.Liu et al. [28] used a double multiple-relaxation-time lattice Boltzmann model to simulate at macroscopic scale transient solid-liquid phase change problems in porous media.The solid-liquid phase change interface is traced through the liquid fraction which is determined by the enthalpy method.Numerical simulations were performed, namely conduction with melting in a semi-infinite space, solidification in a semi-infinite corner, and convection melting in a square cavity filled with porous media. In their simulation, the authors imposed equal values of liquid and solid phase heat capacities.Chen et al. [29] experimentally investigated the melting behavior of a phase change material in metal foams. The authors used Paraffin wax as a phase change material in which aluminum foams were embedded to enhance the heat transfer. A thermal lattice Boltzmann model with doubled populations was implemented to simulate, at the pore scale,the two-dimensional melting of the phase change material in metal foams. The authors assumed that thermo-physical properties of both liquid and solid phases of the PCM are equal and constant. By doing so, the authors omitted afundamental requirement needed in performing heat transfer problems in heterogeneous media, which is the continuity of the normal flux at the interface between the considered phase change material and the foam material. Knowing that, there is obvious difference between the heat capacities of the two considered materials. The conventional lattice Boltzmann formulation is unable to take into account this difference and thereby cannot ensure continuity of the normal heat flux at interfaces.

In this work, we bring a new alternative method based on an enthalpic lattice Boltzmann formulation able to accurately solve phase change material melting by conduction heat transfer and able to implement different thermophysical properties of both involved liquid and solid phases. We start by solving a conduction heat transfer problem during melting of phase change material, and then we investigate the effect of including a discrete number of circular solid blocks inside an enclosure filled with PCM. In all simulated test problems, an enthalpy-based lattice Boltzmann formulation is adopted and comparisons with Control Volume solutions are done.

The present paper is organized as follows; section II presents the proposed enthalpic formulation of the lattice Boltzmann Method. Section III provides the validation of the proposed method through comparisons with Control Volume Method solutions of several test problems. Section IV concludes the paper.

2. Enthalpic Formulation

Considering a two dimensional enclosure filled with phase change material, the governing heat conduction is expressed as [22]:

where:

$\rho=\left(1-f_{l}\right) \rho_{s d}+f_{l} \rho_{l i q}$

$C_{p}=\left(1-f_{l}\right) C_{p, s d}+f_{l} C_{p, l i q}$

$\lambda=\left(1-f_{l}\right) \lambda_{p, s d}+f_{l} \lambda_{p, l i q}$

$S=-\rho_{l i q} L_{M} \frac{\partial f_{l}}{\partial t}$       (1)

where, $T \quad$ and $\quad \lambda$ are the temperature and thermal conductivity, respectively. denotes the volume fraction of the liquid phase.$\rho_{s d}, \rho_{l i q}, C_{p, s d}, C_{p, l i q}$ and $L_{M}$  are density of the solid phase, density of the liquid phase, specific heat of solid phase, specific heat of liquid phase and latent heat of melting, respectively.

In most works investigating phase change heat transfer problem, the volume-phase fraction of the liquid phase, is a function defined on slice intervals marked out by liquid-and solid phase enthalpy [24-27]. In the present study the expression of is taken as a Heaviside step function of temperature, smoothly approximated by a logistic curve having a common "S" shape (also called sigmoid curve):

$f_{l}(T)=\frac{1}{1+e^{-2 \gamma\left(T-T_{M}\right)}}$       (2)

g is the steepness of the so-called sigmoid curve [28]. TM is the melting temperature.The single relaxation time lattice Boltzmann equation for heat conduction problem is [31]

$\begin{aligned} f_{k}\left(\vec{r}+\vec{e}_{k} \Delta t, t+\Delta t\right) &=f_{k}(\vec{r}, t)-\frac{\Delta t}{\tau}\left[f_{k}(\vec{r}, t)-f_{k}^{e q}(\vec{r}, t)\right] \\ &+w_{k} \Delta t \frac{S}{\rho C_{p}} \end{aligned}$       (3)

where, $f_{k}$ is the distribution function,$\tau$  is the relaxation time,$f_{k}^{e q}$ is the equilibrium distribution function and $\Delta t$ is the time step.

In two dimensional problems, the D2Q9 lattice scheme is usually used. and $f_{k}^{e q}$ are as below [31]:

$\tau=\frac{3 \alpha}{\left|\vec{e}_{k}\right|^{2}}+\frac{\Delta t}{2}$       (4)

$f_{k}^{e q}(\vec{r}, t)=w_{k} T(\vec{r}, t)$       (5)

where is the thermal diffusivity. The nine velocities and their corresponding weights $w_{k}$ in the D2Q9 lattice scheme are the following:

$e_{0}=(0,0), e_{1}=( \pm 1,0) \cdot U$

$e_{2,4}=(0, \pm 1) \cdot U, e_{5,6,7,8}=( \pm 1, \pm 1) \cdot U$       (6)

$w_{0}=\frac{4}{9}, w_{1,2,3,4}=\frac{1}{9}, w_{5,6,7,8}=\frac{1}{36}$       (7)

where is the grid size and $U=\frac{\Delta x}{\Delta t}$.Once the values of fk over all directions are known,temperature T is calculated as follows [31]:

$T(\vec{r}, t)=\sum_{k} f_{k}(\vec{r}, t)$       (8)

By using Chapmann-Enskog expansion, the lattice Boltzmann model (Equation 1) retrieves the following heat conduction equation:

$\frac{\partial T}{\partial t}=\vec{\nabla} \cdot(\alpha \vec{\nabla} \cdot T)+\frac{S}{\rho C_{p}}$       (9)

The LBM works with Equation (10), while Equation (1)can be solved by lattice Boltzmann method, only in the case of constant heat capacitance Cp. In our case, pCp is not constant and vary with time and in space.  Among authors who investigated heat transfer problems with phase change, a large number used the conventional lattice Boltzmann method: the diffusivity is evaluated by using pCp even when it is not constant. They often omitted differences in thermophysical properties of liquid and solid phase of the considered phase change material [23-27]. In order to solve Equation (1), using lattice Boltzmann method, we transform it to take the same form as Equation (10). So, we use the enthalpy h rnCp ,nT as a new variable (=sd or liq). By taking h rsd Cp ,sd T , Equation (1) becomes:

$\frac{\partial h}{\partial t}=\vec{\nabla} \cdot(\alpha \vec{\nabla} \cdot h)+S_{t o t}$       (10)

where

$S_{t o t}=\left[1-\frac{\rho_{l i q} C_{p, l i q}}{\rho_{s d} C_{p, s d}}\right] \frac{\partial\left(f_{l} h\right)}{\partial t}+S$       (11)

$\alpha=\frac{\lambda}{\rho_{s d} C_{p, s d}}$       (12)

Consequently, Equation (11) can be solved using the lattice Boltzmann equation with an extra term source[31] Fk wk Stot , by

$\begin{aligned} f_{k}\left(\vec{r}+\vec{e}_{k} \Delta t, t+\Delta t\right) &=f_{k}(\vec{r}, t)-\frac{\Delta t}{\tau}\left[f_{k}(\vec{r}, t)-f_{k}^{e q}(\vec{r}, t)\right] \\ &+\Delta t F_{k} \end{aligned}$       (13)

Once the values of fk over all directions are known, temperature T is calculated as follows [31]:

$T(\vec{r}, t)=\frac{1}{\rho_{s d} C_{p, s d}} \sum_{k} f_{k}(\vec{r}, t)$       (14)

3. Numerical Results

To prove the validity of the proposed formulation of LBM,several conduction heat transfer problems in melting process of a specific phase change material have been investigated,and comparison of obtained isotherms with Control volume solutions, have been done. All the proposed problems were considered inside a square enclosure fulfilled with a selected PCM. In the first test problem, we consider the most common phase change material, water. Two different configurations are adopted, we start by an initial cold temperature equal to $T_{i}=-0.3^{\circ} \mathrm{C}$, and keep heating up the West and South wall of the enclosure by imposing hot temperature $T_{\text {west}}=T_{\text {south}}=+1^{\circ} \mathrm{C}$. The remaining walls (namely, North and East walls) are maintained adiabatic (Figure 1(a)). In the second configuration of melting ice, shown in (Figure 1(b)),only the east wall is transformed into a heated wall, initial and boundary condition are $T_{i}=-0.3^{\circ} \mathrm{C}$ $T_{\text {west}}=T_{\text {south}}=T_{\text {east}}=+1^{\circ} \mathrm{C}$, and the North wall is insulated.The second test problem consists of an enclosure containing a selected commercial PCM [32], in which we insert a circular solid block of natural Graphite [33]. Initially the hole enclosure is kept at $T_{i}=0^{\circ} \mathrm{C}$. North and East wall are insulated and $T_{\text {west}}=T_{\text {south}}=+26^{\circ} \mathrm{C}$. In the last test example performed in this study initial and boundary condition of the previous test is the same, we multiply the number of circular  natural Graphite blocks inside the enclosure filled with the same PCM used in the previous test example. By doing so,we tend to approach a porous medium configuration

3.1 Melting of a phase change material (water)

To demonstrate the feasibility and accuracy of the presented formulation, we begin by solving conduction heat transfer in melting of ice inside a square enclosure. Initially the medium is at cold temperature. Coordinates and boundary conditions used in this test problem are illustrated in Figure 1. Two different set of boundary condition are investigated as shown in Figure 1(a) and (b) respectively.

In Table 1, we provide the thermophysical properties used in this simulation.

Table 1. Thermophysical properties of water

Phase

r

[Kg.m-3 ]

C p

[KJ.Kg-1 K -1 ]

l

[W.m-1 K -1 ]

LM

[KJ.Kg -1 ]

Liquid

1000

4.186

1.6

 

Solid

917

2.05

2.22

333.6

Figure 1. Illustration of heat conduction in melting of ice (a) 1st configuration and (b) 2nd configuration

In Figure 2, comparisons of the obtained lattice Boltzmann results with the Control volume solutions are shown and excellent agreements are found. Isotherms were plotted at several time steps, highlighting the melting evolution of ice during time. In this configuration melting was progressing from a specific corner, while in the Figure 3, where isotherms at two different time steps are plotted, evolution of the melting process is progressing from three sides of the enclosure.Simulation of the same test problems were conducted with the "Conventional" lattice Boltzmann method, where solution of the conduction heat transfer during melting is performed by solving Equation (4) where the diffusivity is evaluated by using pCp even when it is not constant. When compared with the proposed lattice Boltzmann formulation the conventional LBM results show some discrepancy with the control volume solutions literally in all performed test problems.

Figure 2. Comparisons between present LBM source formulation (black solid line), conventional LBM (red solid line) and CVM (symbol) prediction of isotherms at (a)t = 0.1 and (b) t = 0.5 for g = 7 in case of melting of ice from a corner

Figure 3. Comparisons between present LBM source formulation (black solid line) conventional LB (red solid line) and CVM (symbol) prediction of isotherms at (a) t = 0.1 and(b)t = 0.5 for g = 7 in case of melting of ice in a square enclosure

To give an idea of the volume fraction of the liquid phaseused in this paper, we plot the melting front of each conduction test problem performed. In Figure 4, plots of fl for melting ice for the two different configurations mentioned in Figure 1 (a) and (b) are presented.

Since the considered PCM is water, melting temperature is obviously set to $T_{M}=0^{\circ} C$, which is not the case of the next test problem, in which the PCM investigated has a melting temperature equal to $T_{M}=25^{\circ} \mathrm{C}$ [32].

Figure 4. melting front position for g = 7 at t = 0.5 for (a) melting of ice from a corner and (b) melting of ice from three sides at time, present LBM source formulation (black solid line) conventional LB (red solid line) and CVM (symbol)

3.2 Melting from a corner of a PCM containing another solid with circular conducting solid

The second test problem investigated is a conduction heat transfer problem in a square enclosure filled with a phase change material in which we insert a circular conducting solid with different thermophysical properties than those of the PCM. The circular solid block inserted for the current simulated test problem has radius equal to R=0.5.In Table 2, thermophysical properties of PCM and solids inserted are given

Table 2. Thermophysical properties of used materials

 

 

PCM

 

Phase

r

[Kg.m-3 ]

C p

[KJ.Kg-1 K -1 ]

l L

M

[W.m-1 K -1 ]

[KJ.Kg -1 ]

Liqui

d

1510

1

0.55

 

 

185

Solid

1800

2.3

1.5

 

Solid (Natural Graphite)

r

[Kg.m-3 ]

 

1100

C p

[KJ.Kg-1 K -1 ]

0.846

l

[W.m-1 K -1 ]

 

3

The previous enthalpic formulation detailed in section II, have to be modified to take into account the presence of one or more conducting solids of circular form, inside the PCM medium.

The same source term formulation is adopted to solve the current and next conduction problems. The enthalpy of the solid phase change material h rsd Cp,sd T is used as a new variable. In this test problem we consider two different media, the first is the PCM, and the second is the solid block inserted in the center of the PCM region (Figure 5).

By rewriting governing equation for heat conduction in both considered media, we get:

➢ 1 st medium: solid block

$\frac{\partial h}{\partial t}=\vec{\nabla} \cdot(\alpha \vec{\nabla} \cdot h)+S_{t o t}$

$S_{t o t}=\left[1-\frac{\rho_{s b} C_{p, s b}}{\rho_{s d} C_{p, s d}}\right] \frac{\partial(h)}{\partial t}$

$\alpha=\frac{\lambda_{s b}}{\rho_{s d} C_{p, s d}}$     (15)

➢ 2 nd medium: Phase Change Material:

$\frac{\partial h}{\partial t}=\vec{\nabla} \cdot(\alpha \vec{\nabla} \cdot h)+S_{t o t}$

$S_{t o t}=\left[\begin{array}{c}{\rho_{s d} C_{p, b i q}} \\ {\frac{\lambda}{\rho_{s d} C_{p, s d}}}\end{array}\right] \frac{\partial\left(f_{l} h\right)}{\partial t}+S$

$\alpha=\frac{\lambda}{\rho_{s d} C_{p, s d}}$     (16)

For multi-domain problems, continuity of both temperature and normal heat flux should be ensured. The present enthalpic lattice Boltzmann formulation, ensure temperature and normal flux continuity [1]

Figure 5. Schematic diagram of heat conduction in melting of a PCM with one encapsulated circular solid block

Figure 6. Comparisons between present LBM source formulation (black solid line), conventional LBM (red solid line) and CVM (symbol) prediction of isotherms at (a)t = 0.2 , (b) t = 5 for melting of PCM with one conducting circular solid inside a square enclosure

In Figure 6, isotherms are plotted for different time steps, and comparisons are made between the present enthalpic lattice Boltzmann formulation and the control volume solution. An overview of Figure 6 shows good agreements. The circular solid block inserted for the current simulated test problem has radius equal to R=0.5. Likewise the previous heat conduction problems, results given by "conventional"lattice Boltzmann method are also plotted in Figure 6 and clearly they are mismatching with control volume solutions.

3.3 Melting from a corner of a PCM in a quarter space containing discrete circular conducting solids

In this case, precisely, we tend to approach a porous media configuration by including inside the phase change material fulfilling the enclosure, a discrete number of circular conducting solids owing more efficient thermophysical properties than those of the considered PCM. In this context we chose to simulate melting of a commercial PCM inside a square enclosure in which we insert 16 solid blocks (Figure 7), with radius equal to 0.2, whose thermophysical properties are same as in Graphite (natural Graphite [33]) with melting temperature  of  the  PCM  equal  to = 25°C . Initially at (t = 0) the medium is at T = 0°C , West and South wall of  the enclosure are kept at Twest Tsouth = +26°C .

Figure 7. Schematic diagram of heat conduction in melting of a CPM with encapsulated circular solid blocks

Figure 8. melting front position for melting of PCM with 16 conducting circular solids inside a square enclosure

In Figure 8, plot of the volume fraction of the liquid phase is established, and as time evolves, the process of PCM melting is on progress starting from the South-West corner. In Figure 9, we compare results from CVM solutions with those obtained by the proposed LB formulation. Very good agreements are observed

Figure 9. Comparisons between present LBM source formulation (black solid line) Conventional LBM (red solid line) and CVM (symbol) prediction of isotherms at (a)t = 0.2 and (b) t = 5 for melting of PCM with 16 conducting circular solids inside a square enclosure

Figure 9 provides evolution of melting process of the selected PCM with 16 solid blocks inserted, performed by three formulations, namely "conventional" LB, proposed enthalpic LB formulation and control volume solutions. At all plotted time steps, the proposed enthalpic LB formulation agrees well with the CVM solutions, while the conventional LB shows noticeable discrepancies with the CVM solutions.

Considering the fact that the conventional lattice Boltzmann formulation causes notable discrepancies with the control volume solution, solving heat transfer problems with phase change in porous media with the “conventional” LB formulation as in [29] seems to be inappropriate.

4. Conclusions

In this paper, an enthalpy-based lattice Boltzmann formulation was adopted to solve melting problems governed by the heat conduction equation. Test examples including melting of ice from a corner and from three sides of a square enclosure, melting of a commercial PCM with an inserted circular conducting solid at the centre of the enclosure, melting of the same PCM with 16 circular conducting solids, were illustrated. All numerical results obtained by the present method agree very well with the control volume method solution. Obtained results also show the inability of the usually used "conventional" lattice Boltzmann method to accurately simulate conduction heat transfer in melting process of phase change material. The proposed enthalpic LB formulation possesses thereby a great potential to solve heat conduction problem with irregular and evolving interface over time. It can naturally be extended to multiple relaxation time lattice Boltzmann model.

Nomenclature

Cp

specific heat capacity, (J / KgK )

ek

lattice velocity vector in direction k

fk

particle distribution function in the kthdirection

f eq k

equilibrium distribution function in the kthdirection

fl

Volume fraction of liquid phase

h

enthalpy

LM

Latent heat of melting

r

position vector (x, y), m

S

Source term

t

time, s

T

Temperature

TM

melting temperature

wk

weight factor in the kth direction

x, y

axial coordinates

Greek symbols

a

thermal diffusivity, (m2 s)

g

steepness of the sigmoid curve

l

thermal conductivity, W / mK

r

density, Kg / m3

t

relaxation time, s

Dx

space step

Dt

time step

Subscripts

liq

liquid phase

sb

Solid block

sd

solid phase

  References

[1] Hamila R., Nouri M., Ben Nasrallah S., Perré P.(2016). Enthalpic lattice Boltzmann formulation for unsteady heat conduction in heterogeneous media,International Journal of Heat and Mass Transfer, Vol.100, pp. 728–736. DOI:10.1016/j.ijheatmasstransfer.2016.05.001

[2] Bespalko D., Pollard A., Uddin M. (2012). Analysis of the pressure fluctuations from an LBM simulation of turbulent channel flow, Computational Fluids, Vol. 54,pp. 143-146. DOI: 10.1016/j.compfluid.2011.10.008

[3] Yu D., Mei R., Luo L.S., Shyy W. (2003). Viscous flow computations with the method of lattice Boltzmann equation, Progress in Aerospace sciences,Vol. 39, pp. 329-367. DOI: 10.1016/S0376-0421(03)00003-4

[4] Benzi R., Succi S., Vergassola M. (1992). The lattice Boltzmann equation: theory and applications, Physics Reports, Vol. 222, pp. 145-197. DOI: 10.1016/0370-1573(92)90090-M

[5] Qian Y.H., Succi S., Orszag S.A. (1995). Recent advances in lattice Boltzmann computing, in: Annual Reviews of Computational Physics III, pp. 195-242.DOI: 10.1142/9789812830647_0006

[6] Chen S.Y., Doolen G.D. (1998). Lattice Boltzmann method for fluids flows, Annual Review of Fluid Mechanics, Vol. 30, pp. 329-364. DOI:10.1146/annurev.fluid.30.1.329

[7] Aidun C.K., Clausen J.R. (2010). Lattice Boltzmann method for complex flows, Annual Review of Fluid Mechanics, Vol. 42, pp. 439-472. DOI:10.1146/annurev-fluid-121108-145519

[8] Huang X.P., Chen Z.Q., Shi J. (2016). Simulation of solid-liquid phase transition process in aluminum foams using the Lattice Boltzmann method,International Journal of Heat and Technology, Vol.34, No. 4, pp. 694-700. DOI: 10.18280/ijht.340420

[9] Yao S.G., Lei L., Deng J.W., Lu S., Zhang W. (2015).Heat transfer mechanism in porous copper foam wick heat pipes using nanofluids, International Journal of Heat and Technology, Vol. 33, No. 3, pp. 133-138.DOI: 10.18280/ijht.330320

[10] Yuang H.Z., Niu X.D., Shu S., Li M.J., Yamaguchi H.(2014). A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating a flexible filament in an incompressible flow, Computers & Mathematics with Applications, Vol. 67, pp. 1039-1056. DOI: 10.1016/j.camwa.2014.01.006

[11] Mishra S.C., Lankadasu A. (2005). Application of the lattice Boltzmann method for solving the energy equation of a 2-D transient conduction-radiation problem, Numerical Heat Transfer A, Vol. 47, pp.935-954. DOI:10.1016/j.ijheatmasstransfer.2004.10.041

[12] Massailoi F., Benzi R., Succi S. (1993). Exponential tails in two-dimensional Rayleigh-Bénard convection,Europhysics Letters, Vol. 212, pp. 303-310. DOI:iopscience.iop.org/0295-5075/21/3/009)

[13] Mishra S.C., Roy H.K. (2007), Solving transient conduction and radiation heat transfer problems using the lattice Boltzmann method and the finite volume method, Journal of Computational Physics, Vol. 223,pp. 89-107. DOI: 10.1016/j.jcp.2006.08.021

[14] Teixeira C.M. (1998). Incorporating Turbulence Models into the Lattice-Boltzmann Method,International Journal of Modern Physics C, Vol. 09,No. 08, p. 1159. DOI: 10.1142/S0129183198001060

[15] Yu H., Girimaji S.S., Luo L.S. (2005). DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method,Journal of Computational Physics, Vol. 209, pp. 599-616. DOI: 10.1016/j.jcp.2005.03.022

[16] Nie X., Doolen G.D, Chen S. (2002). Lattice Boltzmann simulation of fluid flows in MEMS,Journal of Statistical Physics, Vol. 107, pp. 279-289.DOI: 10.1023/A:1014523007427

[17] Lim C.Y., Shu C., Niu X.D., Chew Y.T. (2002).Application of lattice Boltzmann method to simulate microchannel flows, Phys. Fluids, Vol. 14, pp. 2299-2308. DOI: 10.1063/1.1483841

[18] Guo Z., Zhao T.S. (2002). Lattice Boltzmann model for incompressible flows through porous media,Physical Review E, Vol. 66, p. 036304. DOI:10.1103/PhysRevE.66.036304

[19] Chen L., Fang W., Kang Q., De'Haven H.J.,Viswanathan H.S., Tao W.Q. (2015). Generalized lattice Boltzmann model for flow through tight porous media with Klinkenberg's effect, Physical Review. E,Vol. 91, p. 033004. DOI:10.1103/PhysRevE.91.033004

[20] Hu Y., Li D.C., Shu S., Niu X.D. (2016). Immersed boundary-lattice Boltzmann simulation of natural convection in a square enclosure with a cylinder covered by porous layer, International Journal of Heat and Mass transfer, Vol. 92, pp. 1166-1170. DOI:10.1016/j.ijheatmasstransfer.2015.09.034

[21] Chen X., Chen H. (1993). Lattice Boltzmann model for simulating flows with multiple phases and components, Physical Review E, Vol. 47, p. 1815.DOI: 10.1103/PhysRevE.47.1815

[22] He X.Y., Chen S.Y., Zhang R. (1999). A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability, Journal of Computational Physics, Vol.152, p. 642-663. DOI: 10.1006/jcph.1999.6257

[23] Miller W., Succi S., Mansutti D. (2001). Lattice Boltzmann model for anisotropic liquid-solid phase transition, Physical Review Letter, Vol. 86, pp. 3578-3581. DOI: 10.1103/PhysRevLett.86.3578

[24] Jiaung W.S., Ho J.R., Kuo C.P. (2001). Lattice Boltzmann method for the heat conduction problem with phase change, Numerical heat transfer, Part B, Fundamentals, Vol. 39, pp. 167-187. DOI:10.1080/10407790150503495

[25] Huang R.Z., Wu H.Y. (2014). An immersed boundarythermal lattice Boltzmann method for solid-liquid phase change, Journal of Computational Physics, Vol.277, pp. 305-319. DOI: 10.1016/j.jcp.2014.08.020

[26] Li Z., Yang M., Zhang Y. (2015). Numerical simulation of melting problems using the lattice Boltzmann method with the interfacial tracking method,Numerical Heat Transfer, Part A, Vol. 68, pp. 1175-1197. DOI: 10.1080/10407782.2015.1037126

[27] Luo K., Yao F.J., Yi H.L., Tan H.P. (2015). Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials, Applied Thermal Engineering, Vol. 86, pp.238-250. DOI: 10.1016/j.applthermaleng.2015.04.059

[28] Huang R., Wu H., Cheng P. (2013). A new lattice Boltzmann model for solid-liquid phase change,International Journal of Heat and Mass Transfer, Vol.59, pp. 295-301. DOI:10.1016/j.ijheatmasstransfer.2012.12.027

[29] Huang R., Wu H. (2016). Total enthalpy-based lattice Boltzmann method with adaptive mesh refinement for solid-liquid phase change, Journal of Computational  Physics, Vol. 315, pp. 65-83. DOI:10.1016/j.jcp.2016.03.043

[30] Liu Q., He Y.L. (2015). Double multiple-relaxationtime lattice Boltzmann model for solid-liquid phase change with natural convection in porous media,Physica A: Statistical Mechanics and its Applications,Vol. 438, pp. 94–106. DOI:10.1016/j.physa.2015.06.018

[31] Chen Z., Gao D., Shi J. (2014). Experimental and numerical study on melting of phase change materials in metal foams at pore scale, International Journal of Heat and Mass Transfer, Vol. 72, pp. 646-655. DOI:10.1016/j.ijheatmasstransfer.2014.01.003

[32] Guichard S. (2016). Contribution a l'Etude des Parois Complexes intégrant des Matériaux a Changement de Phase: Modélisation, Expérimentation et Evaluation de la performance énergétique globale, Thesis dissertation, Réunion University, France, tel.archivesouvertes.fr/tel-01379518

[33] Mohamad A.A. (2011). Lattice Boltzmann method:fundamentals and engineering applications with computer codes, Springer-Verlag London Limited.DOI: 10.1007/978-0-85729-455-5

[34] http://www.pluss.co.in/technical-datasheets/Doc303- TDS-HS-24.pdf, acceded on Avril 2016.

[35] Smalc M., Shives G., Chen G., Guggari S., Norley J.,Andy Reynolds III R. (2005). Thermal performance of natural graphite heat spreaders, Proceedings of IPACK 2005, ASME InterPACK ‘05 July17-22, San Francisco,California USA.