Simulating Solid-Liquid Phase-Change Heat Transfer in Metal Foams via a Cascaded Lattice Boltzmann Model

Abstract

A cascaded lattice Boltzmann (CLB) model is constructed for simulating heat transfer in metal-foam-based solid-liquid phase change materials (PCMs). The present model captures the phase interface implicitly via the enthalpy methodology, and to avoid iterations in simulations, the CLB equation of the PCM employs the enthalpy as the basic evolution variable through modifying the cascaded collision process. Numerical results demonstrate the effectiveness and practicability of the CLB model for investigating heat transfer in solid-liquid PCMs with metal foams. The effects of the inertial coefficient, permeability and porosity on the melting process are investigated. The results indicate that the empirical correlations of inertial coefficient and permeability based on packed beds overestimate the melting rate at high porosities. Moreover, the porosity has significant impact on phase-change processes. The melting rate increases as the porosity of the metal foam decreases since heat conduction through high thermal conductive metal foam dominates the total heat transfer.

1. Introduction

Latent heat storage (LHS) using solid-liquid PCMs attracts wide attention from the academic research community and industry field due to its importance for solving energy and environment issues [1,2,3,4]. Although solid-liquid PCMs have some outstanding advantages (i.e., high energy storage density), they usually suffer from low thermal conductivities (in the range of 0.1~0.6 W/(m·K) [5]), which strongly affects the LHS system’s thermal efficiency. To enhance PCMs’ thermal conductivities, embedding high thermal conductive metal foams within PCMs to form composite PCMs has long been practiced [6,7]. In the past several decades, various conventional numerical methods (e.g., FVM [8,9]) have been developed to investigate thermal behaviors of heat transfer in solid-liquid PCMs with metal foams. The lattice Boltzmann (LB) method [10], as a particle-based numerical tool evolved from lattice-gas automata [11], has also attracted great attention in studying heat transfer in solid-liquid PCMs [12].

The LB method is a mesoscopic modeling method sitting in the intermediate zone between macroscopic continuum-based methods and microscopic molecular dynamics method. This scale-bridging nature is a fundamental feature of the LB method, and consequently, it has several distinctive advantages, such as local nature of operations, ideal for parallel computing, avoiding solving Poisson equation for incompressible flows, and easy implementation of complex boundary conditions with elementary mechanical rules [13,14]. Due to its inherent advantage, the LB method is particularly effective in studying transient phase-change processes with complex interfacial dynamics and moving phase interface. Moreover, for transient phase-change processes, the efficiency of the LB method is higher than that of the conventional numerical methods if time steps are equal. In 1998, Fabritiis et al. [15] developed the first LB scheme for investigating solid-liquid phase change. Since then, many LB models have been constructed for transient solid-liquid phase-change problems, and these models can be divided into three major groups: phase-field approach [16,17], enthalpy-based approach [18,19,20], and interface-tracking approach [21,22]. Among the three types of phase-change LB models, the enthalpy-based approach has attracted great attention owing to its effectiveness and easy implementation in modeling solid-liquid phase-change processes. In the enthalpy-based approach, the diffusive phase interface as well as the liquid and solid phases are distinguished via liquid fraction, which makes this approach simple and effective. In the past decade, the enthalpy-based approach has also been successfully adopted to investigate solid-liquid phase change in porous media by many researchers [23,24,25,26,27,28]. These enthalpy-based LB models are mainly constructed based on Bhatnagar–Gross–Krook (BGK) [23,24,25] and multiple-relaxation-time (MRT) [26,27,28] models. As reported in Ref. [27], the effect of unphysical numerical diffusion on numerical simulations can be serious in the BGK model, while the MRT model can effectively eliminate such effect since it has sufficient tunable relaxation parameters.

In addition to BGK and MRT collision models, the cascaded collision model [29,30,31,32] also attracts significant attention in the LB community. The cascaded collision model was first developed by Geier et al. [29] in the year of 2006. In the framework of this model, the collision step is executed based on central moments, and different central moments are relaxed to their equilibria with different relaxation rates. In the CLB method, a high degree of Galilean invariance can be preserved since the influence between different orders of moments has been eliminated [29,30]. For heat transfer in solid-liquid PCMs with metal foams, the interfacial behaviors in phase-change region are very important, and interfacial heat transfer (thermal non-equilibrium effect) plays a significant role because metal foam’s thermal conductivity is usually three orders of magnitude higher than PCM’s thermal conductivity. An in-depth understanding of the complex interfacial dynamics and thermal behaviors requires efficient and powerful numerical tool. Since the cascaded collision model has tunable relaxation parameters as well as a high degree of Galilean invariance, it is expected that the complex interfacial dynamics and thermal behaviors of heat transfer in metal-foam-based solid-liquid PCMs can be well captured by using the cascaded collision model. Hence, this work aims to propose a CLB model for simulating heat transfer in solid-liquid PCMs with metal foams, in which the enthalpy methodology is adopted to capture the phase interface implicitly. Moreover, the effects of the inertial coefficient, permeability and porosity on melting process will be investigated.

2. Macroscopic Governing Equations

For heat transfer in metal-foam-based solid-liquid PCMs, the macroscopic governing equations are given by [9,12,23]:

$$\mathbf{\nabla }·\mathbf{u}=0$$

(1)

$$\frac{\partial \mathbf{u}}{\partial t}+\left(\mathbf{u}·\mathbf{\nabla }\right)\left(\frac{\mathbf{u}}{\varphi }\right)=-\frac{1}{{\rho }_f}\mathbf{\nabla }\left(\varphi p\right)+v_e{\mathbf{\nabla }}^2\mathbf{u}+\mathbf{F}$$

(2)

$$\frac{\partial }{\partial t}\left(\varphi {\rho }_f{c}_{pf}{T}_f\right)+\mathbf{\nabla }·\left({\rho }_l{c}_{pl}{T}_f\mathbf{u}\right)=\mathbf{\nabla }·\left[\left(k_{ef}+k_d\right)\mathbf{\nabla }{T}_f\right]+h_{mf}{a}_{mf}\left(T_m-T_f\right)-\underset{¯}{\frac{\partial }{\partial t}\left(\varphi {\rho }_l{L}_a{f}_l\right)}$$

(3)

$$\frac{\partial }{\partial t}\left[\left(1-\varphi \right){\rho }_m{c}_{pm}{T}_m\right]=\mathbf{\nabla }·\left(k_{em}\mathbf{\nabla }{T}_m\right)+h_{mf}{a}_{mf}\left(T_f-T_m\right)$$

(4)

Equations (1) and (2) are the so-called generalized non-Darcy equations. The last term of Equation (3) is the phase-change term ($f_l=0$ denotes the solid region, $0<f_l<1$ denotes the phase-change region, and $f_l=1$ denotes the liquid region).

The total body force $\mathbf{F}$ is determined by:

$$\mathbf{F}=-\frac{\varphi v_{fl}}{K}\mathbf{u}-\frac{\varphi F_{\varphi }}{\sqrt{K}}\left|\mathbf{u}\right|\mathbf{u}+\varphi \mathbf{G}$$

(5)

According to the Boussinesq approximation, the body force $\mathbf{G}$ is given by $\mathbf{G}=-\mathbf{g}\beta \left(T-T_0\right)$.

For flow through packed bed consisting of spherical particles, $F_{\varphi }$ and $K$ are given by:

$$F_{\varphi }=\frac{1.75}{\sqrt{150{\varphi }^3}}, K=\frac{{\varphi }^3{d}_p^2}{150{\left(1-\varphi \right)}^2}$$

(6)

$d_p^{}$ in Equation (6) denotes the mean particle diameter.

For metal foams such as aluminum foam, $F_{\varphi }$ and $K$ can be well predicted by the following empirical correlations [33]:

$$F_{\varphi }=A_1{\left(1-\varphi \right)}^{b_1}{\left(\frac{d_f}{d_p}\right)}^{c_1}, K=A_2{\left(1-\varphi \right)}^{b_2}{d}_p^2{\left(\frac{d_f}{d_p}\right)}^{c_2}$$

(7)

where $A_1=0.00212$, $b_1=-0.132$, $c_1=-1.63$, $A_2=0.00073$, $b_2=-0.224$, and $c_2=-1.11$. $d_p^{}$ in Equation (7) denotes the mean diameter of the pores.

An appropriate equation for the structure of the metal foam is given by [33]:

$$\frac{d_f}{d_p}=1.18\sqrt{\frac{1-\varphi }{3\pi }}\frac{1}{1-e^{-(1-\varphi )/0.04}}$$

(8)

The effective thermal conductivities $k_{ef}$ and $k_{em}$ can be determined by analytical models [33,34]. Since the metal foam’s thermal conductivity is rather high, the thermal dispersion can be neglected and $k_d$ is usually set to in Equation (3) [8].

$a_{mf}$ can be predicted by the following relation [35]:

$$a_{mf}=\frac{3\pi d_f}{{\left(0.59d_p\right)}^2}\left[1-e^{-(1-\varphi )/0.04}\right]$$

(9)

$h_{mf}$ depends on the foam structure, Reynolds number, and Prandtl number. For forced convective flow in metal foams, several empirical correlations for $h_{mf}$ have been proposed [8,34]. For natural convection flow through metal foams, the following empirical correlation is widely used [35]:

$$h_{mf}=\frac{k_f}{d_f}\left(0.36+\frac{0.518R{a}_d^{1/4}}{{\left[1+{\left(0.559/Pr\right)}^{9/16}\right]}^{4/9}}\right)$$

(10)

where $R{a}_d=g\beta \Delta T{d}_f^3/({\alpha }_{fl}{v}_{fl})$, in which ${\alpha }_{fl}=k_{fl}/\left({\rho }_l{c}_{pl}\right)$.

3. Enthalpy-Based CLB Model for Solid-Liquid Phase Change in Metal Foams

For the flow field described by Equations (1) and (2), the isothermal CLB equation in Ref. [32] is adopted. In what follows, the enthalpy-based CLB equation for PCM’s temperature field and the thermal MRT-LB equation for metal foam’s temperature field are presented in detail. The D2Q5 lattice is adopted, in which the discrete velocities $\left\{e_i|i=0, \dots , 4\right\}$ are given by:

$$e_i=\{\begin{array}{ll}\left(0,0\right),\hfill & i=0\hfill \\ \left(\cos \left[\left(i-1\right)\pi /2\right],\sin \left[\left(i-1\right)\pi /2\right]\right)c,\hfill & i=1~4\hfill \end{array}$$

(11)

3.1. Enthalpy-Based CLB Equation

3.1.1. Thermal CLB Equation for Fluid Phase without Phase-Change Term

Without the phase-change term, Equation (3) can be rewritten as:

$$\frac{\partial }{\partial t}\left(c_{pf}{T}_f\right)+\mathbf{\nabla }·\left(\frac{c_{pf}{T}_f\mathbf{u}}{\varphi }\right)=\mathbf{\nabla }·\left(\frac{k_{ef}}{\varphi {\rho }_f}\mathbf{\nabla }{T}_f\right)+\frac{h_{mf}{a}_{mf}\left(T_m-T_f\right)}{\varphi {\rho }_f}$$

(12)

We first introduce the thermal CLB equation for solving Equation (3) without phase-change term based on the simplified CLB method [30]. The raw moments $\left\{k_{mn}^T\right\}$ and central moments $\left\{{\tilde{k}}_{mn}^T\right\}$ of the discrete temperature distribution function $g_{fi}$ of the fluid phase are defined as follows [30]:

$$k_{mn}^T=\langle g_{fi}|e_{ix}^m{e}_{iy}^n\rangle , {\tilde{k}}_{mn}^T=\langle g_{fi}|{\left(e_{ix}^{}-u_x\right)}^m{\left(e_{iy}^{}-u_y\right)}^n\rangle$$

(13)

The simplified raw-moment $|n_{fi}^{}\rangle$ and central-moment $|{\tilde{n}}_{fi}^{}\rangle$ are given by [30]

$$|n_{fi}^{}\rangle ={\left[k_{00}^T,k_{10}^T,k_{01}^T,k_{20}^T,k_{02}^T\right]}^{{}^T}, |{\tilde{n}}_{fi}^{}\rangle ={\left[{\tilde{k}}_{00}^T,{\tilde{k}}_{10}^T,{\tilde{k}}_{01}^T,{\tilde{k}}_{20}^T,{\tilde{k}}_{02}^T\right]}^{{}^T}$$

(14)

Through the transformation matrix ${\mathbf{M}}_T$ and shift matrix ${\mathbf{N}}_T$, $\left\{n_{fi}^{}\right\}$ and $\left\{{\tilde{n}}_{fi}^{}\right\}$ can be determined by [30]

$$|{\tilde{n}}_{fi}^{}\rangle ={\mathbf{N}}_T|n_{fi}^{}\rangle , |n_{fi}^{}\rangle ={\mathbf{M}}_T|g_{fi}^{}\rangle$$

(15)

The thermal CLB equation for solving Equation (12) can be divided into two parts: collision and streaming processes. Through shifting procedure ($|{\tilde{n}}_{fi}^{}\rangle ={\mathbf{N}}_T|n_{fi}^{}\rangle$), collision process is implemented in central-moment space as

$$|{\tilde{n}}_{fi}^{*}\rangle =|{\tilde{n}}_{fi}^{}\rangle -\mathbf{\Theta }\left(|{\tilde{n}}_{fi}^{}\rangle -|{\tilde{n}}_{fi}^{eq}\rangle \right)+|{\tilde{S}}_{fi}^{}\rangle$$

(16)

where $\left\{{\tilde{n}}_{fi}^{eq}\right\}$ are equilibrium central moments, $\left\{{\tilde{n}}_{fi}^{*}\right\}$ denotes the post-collision central moments, $\mathbf{\Theta }=\mathrm{diag}\left({\zeta }_T,{\zeta }_{\alpha },{\zeta }_{\alpha },{\zeta }_e,{\zeta }_e\right)$ is the relaxation matrix, $|{\tilde{S}}_{fi}^{}\rangle$ is the source term.

The streaming process is still implemented in velocity space as

$$g_{fi}\left(\mathbf{x}+{\mathbf{e}}_i{\delta }_t, t+{\delta }_t\right)=g_{fi}^{*}\left(\mathbf{x}, t\right)$$

(17)

where $g_{fi}^{*}$ is determined by $|g_{fi}^{*}\rangle ={\mathbf{M}}_T^{-1}{\mathbf{N}}_T^{-1}|{\tilde{n}}_{fi}^{*}\rangle$ ($|g_{fi}^{*}\rangle ={\mathbf{M}}_T^{-1}|n_{fi}^{*}\rangle$ and $|n_{fi}^{*}\rangle ={\mathbf{N}}_T^{-1}|{\tilde{n}}_{fi}^{*}\rangle$).

${\mathbf{M}}_T$ and ${\mathbf{N}}_T$ are given by ($c=1$):

$${\mathbf{M}}_T=\left[\begin{array}{rrrrr}\hfill 1& \hfill 1& \hfill 1& \hfill 1& \hfill 1\\ \hfill 0& \hfill 1& \hfill 0& \hfill -1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1& \hfill 0& \hfill -1\\ \hfill 0& \hfill 1& \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1& \hfill 0& \hfill 1\end{array}\right]$$

(18)

$${\mathbf{N}}_T=\left[\begin{array}{rrrrr}\hfill 1& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill -u_x/\varphi & \hfill 1& \hfill 0& \hfill 0& \hfill 0\\ \hfill -u_y/\varphi & \hfill 0& \hfill 1& \hfill 0& \hfill 0\\ \hfill u_x^2& \hfill -2\varphi u_x& \hfill 0& \hfill 1& \hfill 0\\ \hfill u_y^2& \hfill 0& \hfill -2\varphi u_y& \hfill 0& \hfill 1\end{array}\right]$$

(19)

Explicitly, the equilibrium central moments $\left\{{\tilde{n}}_{fi}^{eq}\right\}$ and raw moments $\left\{n_{fi}^{eq}\right\}$ are given by:

$$|{\tilde{n}}_{fi}^{eq}\rangle ={\left[c_{pf}{T}_f,0,0,c_{sf}^2{c}_{pf}{T}_f,c_{sf}^2{c}_{pf}{T}_f\right]}^{\mathsf{T}}$$

(20)

$$|n_{fi}^{eq}\rangle ={\left[c_{pf}{T}_f,\frac{c_{pf}{T}_f{u}_x}{\varphi },\frac{c_{pf}{T}_f{u}_y}{\varphi },c_{sf}^2{c}_{pf}{T}_f+u_x^2{c}_{pf}{T}_f,c_{sf}^2{c}_{pf}{T}_f+u_y^2{c}_{pf}{T}_f\right]}^{\mathsf{T}}$$

(21)

The equilibrium temperature distribution function $g_{fi}^{eq}$ can be determined by $|g_{fi}^{eq}\rangle ={\mathbf{M}}_T^{-1}|n_{fi}^{eq}\rangle$.

Explicitly, the source term $|{\tilde{S}}_{fi}^{}\rangle$ is given by:

$$|{\tilde{S}}_{fi}^{}\rangle ={\left[S{r}_f+\frac{1}{2}{\delta }_t{\partial }_{t}S{r}_f,0,0,0,0\right]}^T$$

(22)

where $S{r}_f=h_{mf}{a}_{mf}\left(T_m-T_f\right)/\left(\varphi {\rho }_f\right)$.

The temperature $T_f$ is defined as:

$$c_{pf}{T}_f={\displaystyle \sum _{i=0}^4{g}_{fi}}$$

(23)

The effective thermal diffusivity ${\alpha }_{ef}$ is:

$${\alpha }_{ef}=\frac{k_{ef}}{\varphi {\rho }_f{c}_{pf}}=c_{sf}^2\left({\mathsf{\zeta }}_{\alpha }^{-1}-\frac{1}{2}\right){\delta }_t$$

(24)

3.1.2. Enthalpy-Based CLB Equation for PCM with Phase Change

The enthalpy-based CLB equation for solving Equation (3) is developed based on the thermal CLB equation in Section 3.1.1. By combining the phase-change term ${\partial }_t\left(\varphi {\rho }_l{L}_a{f}_l\right)$ with the transient term ${\partial }_t\left(\varphi {\rho }_f{c}_{pf}{T}_f\right)$ in Equation (3), we can obtain:

$$\frac{\partial H_f}{\partial t}+\mathbf{\nabla }·\left(\frac{c_{pl}{T}_f\mathbf{u}}{\varphi }\right)=\mathbf{\nabla }·\left(\frac{k_{ef}}{\varphi {\rho }_l}\mathbf{\nabla }{T}_f\right)+\frac{h_{mf}{a}_{mf}\left(T_m-T_f\right)}{\varphi {\rho }_l}$$

(25)

where $H_f=\sigma c_{pl}{T}_f+L_a{f}_l$ is the enthalpy, and $\sigma =\frac{{\rho }_f{c}_{pf}}{{\rho }_l{c}_{pl}}=\frac{f_l{\rho }_l{c}_{pl}+\left(1-f_l\right){\rho }_s{c}_{ps}}{{\rho }_l{c}_{pl}}$ is the heat capacity ratio. In the liquid phase ($f_l=1$), $H_f=c_{pl}{T}_f+L_a$ and ${\sigma }_l=1$; in the solid phase ($f_l=0$), $H_f={\sigma }_s{c}_{pl}{T}_f$ and ${\sigma }_s=\frac{{\rho }_s{c}_{ps}}{{\rho }_l{c}_{pl}}$. The source term $S{r}_f=h_{mf}{a}_{mf}\left(T_m-T_f\right)/\left(\varphi {\rho }_l\right)$.

Formally, the collision and streaming processes of the CLB equation for solving Equation (25) are still given by Equations (16) and (17), respectively. However, to match the enthalpy-based energy Equation (25), the equilibrium central moments $\left\{{\tilde{n}}_{fi}^{eq}\right\}$ and raw moments $\left\{n_{fi}^{eq}\right\}$ should be chosen as follows:

$$|{\tilde{n}}_{fi}^{eq}\rangle ={\left[H_f,0,0,c_{sf}^2{c}_{pl}{T}_f,c_{sf}^2{c}_{pl}{T}_f\right]}^T$$

(26)

$$|n_{fi}^{eq}\rangle ={\left[H_f,\frac{c_{pl}{T}_f{u}_x}{\varphi },\frac{c_{pl}{T}_f{u}_y}{\varphi },c_{sf}^2{c}_{pl}{T}_f+u_x^2{c}_{pl}{T}_f,c_{sf}^2{c}_{pl}{T}_f+u_y^2{c}_{pl}{T}_f\right]}^T$$

(27)

In the system, ${\tilde{n}}_{f0}^{}$ is the only conserved quantity (${\tilde{n}}_{f0}^{}=n_{f0}^{}=H_f$). As the enthalpy $H_f$ is the basic evolution variable ($g_{fi}$ is now defined as the enthalpy distribution function), the shifting procedure $|{\tilde{n}}_{fi}^{}\rangle ={\mathbf{N}}_T|n_{fi}^{}\rangle$ in the original thermal CLB equation should be modified. Explicitly, the modified shifting procedure is given by:

$$\begin{gathered} {\tilde{n}}_{f0}^{}=n_{f0}^{} \\ {\tilde{n}}_{f1}^{}=-\frac{u_x{c}_{pl}{T}_f}{\varphi }+n_{f1}^{} \\ {\tilde{n}}_{f2}^{}=-\frac{u_y{c}_{pl}{T}_f}{\varphi }+n_{f2}^{} \\ {\tilde{n}}_{f3}^{}=u_x^2{c}_{pl}{T}_f-2\varphi u_x{n}_{f1}^{}+n_{f3}^{} \\ {\tilde{n}}_{f4}^{}=u_y^2{c}_{pl}{T}_f-2\varphi u_y{n}_{f2}^{}+n_{f4}^{} \end{gathered}$$

(28)

Formally, the collision process is still given by Equation (16), but it should be noted that it has been reconstructed since the shifting procedure is modified as Equation (28). The post-collision enthalpy distribution function $g_{fi}^{*}$ is determined via $|g_{fi}^{*}\rangle ={\mathbf{M}}_T^{-1}|n_{fi}^{*}\rangle$, in which $|n_{fi}^{*}\rangle$ is given by:

$$\begin{gathered} n_{f0}^{*}={\tilde{n}}_{f0}^{*} \\ {n}_{f1}^{*}=\frac{u_x{c}_{pl}{T}_f}{\varphi }+{\tilde{n}}_{f1}^{*} \\ {n}_{f2}^{*}=\frac{u_y{c}_{pl}{T}_f}{\varphi }+{\tilde{n}}_{f2}^{*} \\ {n}_{f3}^{*}=u_x^2{c}_{pl}{T}_f+2\varphi u_x{\tilde{n}}_{f1}^{*}+{\tilde{n}}_{f3}^{*} \\ {n}_{f4}^{*}=u_y^2{c}_{pl}{T}_f+2\varphi u_y{\tilde{n}}_{f2}^{*}+{\tilde{n}}_{f4}^{*} \end{gathered}$$

(29)

Explicitly, $g_{fi}^{eq}$ in the velocity space is:

$$\begin{gathered} g_{f0}^{eq}=H_f-{\varpi }_1{c}_{pl}{T}_f-\left(u_x^2+u_y^2\right)c_{pl}{T}_f \\ {g}_{f1}^{eq}=\frac{{\varpi }_1}{4}{c}_{pl}{T}_f+\frac{1}{2}{u}_x{c}_{pl}{T}_f+\frac{1}{2}{u}_x^2{c}_{pl}{T}_f \\ {g}_{f2}^{eq}=\frac{{\varpi }_1}{4}{c}_{pl}{T}_f+\frac{1}{2}{u}_y{c}_{pl}{T}_f+\frac{1}{2}{u}_y^2{c}_{pl}{T}_f \\ {g}_{f3}^{eq}=\frac{{\varpi }_1}{4}{c}_{pl}{T}_f-\frac{1}{2}{u}_x{c}_{pl}{T}_f+\frac{1}{2}{u}_x^2{c}_{pl}{T}_f \\ {g}_{f4}^{eq}=\frac{{\varpi }_1}{4}{c}_{pl}{T}_f-\frac{1}{2}{u}_y{c}_{pl}{T}_f+\frac{1}{2}{u}_y^2{c}_{pl}{T}_f \end{gathered}$$

(30)

where ${\varpi }_1\in \left(0, 1\right)$ and $c_{sf}^{}=c\sqrt{{\varpi }_1/2}=\sqrt{{\varpi }_1/2}$.

The enthalpy $H_f$ is defined by:

$$H_f={\displaystyle \sum _{i=0}^4{g}_{fi}}$$

(31)

Simultaneously, the temperature $T_f$ can be calculated by:

$$T_f=\{\begin{array}{ll}{H}_f/\left({\sigma }_s{c}_{pl}\right),\hfill & H_f\le H_{fs}\hfill \\ T_{fs}+\frac{H_f-H_{fs}}{H_{fl}-H_{fs}}\left(T_{fl}-T_{fs}\right),\hfill & H_{fs}<H_f<H_{fl}\hfill \\ T_{fl}+\left(H_f-H_{fl}\right)/\left({\sigma }_l{c}_{pl}\right),\hfill & H_f\ge H_{fl}\hfill \end{array}$$

(32)

where $T_{fs}$ and $T_{fl}$ ($T_{fs}\le T_{fl}$) are solidus and liquidus temperatures, respectively; $H_{fs}$ and $H_{fl}$ are enthalpies corresponding to $T_{fs}$ and $T_{fl}$, respectively.

The liquid fraction $f_l$ is:

$$f_l=\{\begin{array}{ll}0,\hfill & H_f\le H_{fs}\hfill \\ \frac{H_f-H_{fs}}{H_{fl}-H_{fs}},\hfill & H_{fs}<H_f<H_{fl}\hfill \\ 1,\hfill & H_f\ge H_{fl}\hfill \end{array}$$

(33)

${\alpha }_{ef}$ is defined as ${\alpha }_{ef}=k_{ef}/\left(\varphi {\rho }_l{c}_{pl}\right)=c_{sf}^2\left({\mathsf{\zeta }}_{\alpha }^{-1}-0.5\right){\delta }_t$.

3.2. Thermal MRT-LB Equation

Equation (4) can be rewritten as:

$$\frac{\partial \left(c_{pm}{T}_m\right)}{\partial t}=\mathbf{\nabla }·\left(\frac{k_{em}}{\left(1-\varphi \right){\rho }_m}\mathbf{\nabla }{T}_m\right)+\frac{h_{mf}{a}_{mf}\left(T_f-T_m\right)}{\left(1-\varphi \right){\rho }_m}$$

(34)

The MRT-LB equation for Equation (34) is:

$$g_{mi}\left(\mathbf{x}+{\mathbf{e}}_i{\delta }_t, t+{\delta }_t\right)-g_{mi}\left(\mathbf{x}, t\right)=-{\left({\mathbf{M}}_T^{-1}\mathbf{Q}{\mathbf{M}}_T\right)}_{ij}\left(g_{mj}-g_{mj}^{eq}\right)|{}_{\left(x, t\right)}+{\delta }_t{S}_{mi}^{}$$

(35)

where $g_{mi}\left(\mathbf{x}, t\right)$ is the temperature distribution function of the metal foam, $g_{mi}^{eq}\left(\mathbf{x}, t\right)$ is the equilibrium of $g_{mi}\left(\mathbf{x}, t\right)$, and $\mathbf{Q}=\mathrm{diag}\left({\eta }_T,{\eta }_{\alpha },{\eta }_{\alpha },{\eta }_e,{\eta }_e\right)$ is the relaxation matrix.

The collision process of the evolution Equation (35) is implemented in raw-moment space as:

$$|n_{mi}^{*}\rangle =|n_{mi}^{}\rangle -\mathbf{Q}\left(|n_{mi}^{}\rangle -|n_{mi}^{eq}\rangle \right)+|{\tilde{S}}_{mi}^{}\rangle$$

(36)

where $|n_{mi}^{}\rangle ={\mathbf{M}}_T|g_{mi}^{}\rangle$, $|n_{mi}^{eq}\rangle ={\mathbf{M}}_T|g_{mi}^{eq}\rangle$, and $|{\tilde{S}}_{mi}^{}\rangle ={\mathbf{M}}_T|S_{mi}^{}\rangle$.

The streaming process is executed in velocity space as

$$g_{mi}\left(\mathbf{x}+{\mathbf{e}}_i{\delta }_t, t+{\delta }_t\right)=g_{mi}^{*}\left(\mathbf{x}, t\right)$$

(37)

where $|g_{mi}^{*}\rangle ={\mathbf{M}}_T^{-1}|n_{mi}^{*}\rangle$.

The equilibrium moment $|n_{mi}^{eq}\rangle$ is:

$$|n_{mi}^{eq}\rangle ={\left[c_{pm}{T}_m,0,0,\frac{\varpi c_{pm}{T}_m}{2},\frac{\varpi c_{pm}{T}_m}{2}\right]}^{\mathsf{T}}$$

(38)

The equilibrium temperature distribution function $g_{mi}^{eq}\left(\mathbf{x}, t\right)$ in the velocity space is given by:

$$g_{mi}^{eq}=\{\begin{array}{ll}\left(1-{\varpi }_2\right)c_{pm}{T}_m,\hfill & i=0\hfill \\ \frac{1}{4}{\varpi }_2{c}_{pm}{T}_m,\hfill & i=1~4\hfill \end{array}$$

(39)

where ${\varpi }_2\in \left(0, 1\right)$ is a tunable parameter of the D2Q5 lattice model.

Explicitly, the source term $|{\tilde{S}}_{mi}^{}\rangle$ is given by:

$$|{\tilde{S}}_{mi}^{}\rangle ={\left[S{r}_m+\frac{1}{2}{\delta }_t{\partial }_{t}S{r}_m,0,0,0,0\right]}^{\mathsf{T}}$$

(40)

where $S{r}_m=h_{mf}{a}_{mf}\left(T_f-T_m\right)/\left[\left(1-\varphi \right){\rho }_m\right]$.

The temperature $T_m$ is defined by:

$$c_{pm}{T}_m={\displaystyle \sum _{i=0}^4{g}_{mi}}$$

(41)

The effective thermal diffusivity ${\alpha }_{em}=k_{em}/\left[\left(1-\varphi \right){\rho }_m\right]=c_{sm}^2\left({\eta }_{\alpha }^{-1}-0.5\right){\delta }_t$, in which $c_{sm}=c\sqrt{{\varpi }_2/2}=\sqrt{{\varpi }_2/2}$ is the sound speed.

4. Numerical Simulations

In this section, the present enthalpy-based CLB model is employed to investigate solid-liquid phase change (melting) with natural convection in metal foams. The schematic of the problem under investigation is shown in Figure 1. Initially, the PCM (solid state) and metal foam are kept at equilibrium temperature $T_i$ ($T_i\le T_{\mathrm{melt}}$). At $t=0$, the temperature of the left wall is raised to $T_h$ (T_h > T_melt), and consequently, the PCM starts to melt. The characteristic parameters are defined as follows:

$$\begin{gathered} Da=\frac{K}{L^2}, Ra=\frac{g\beta \Delta T{L}^3}{v_{fl}{\alpha }_{fl}}, Pr=\frac{v_{fl}}{{\alpha }_{fl}}, J=\frac{v_e}{v_{fl}}, \widehat{\sigma }=\frac{{\rho }_m{c}_{pm}}{{\rho }_l{c}_{pl}}, \lambda =\frac{k_m}{k_{fl}} \\ \Gamma =\frac{{\alpha }_m}{{\alpha }_{fl}}, H_v=\frac{h_{mf}{a}_{mf}{d}_p^2}{k_f}, Fo=\frac{t{\alpha }_{fl}}{L^2}, St=\frac{c_{pl}\Delta T}{L_a} \end{gathered}$$

(42)

where $\Delta T=T_h-T_c$ ($T_h>T_c$) is the characteristic temperature and ${\alpha }_m=k_m/\left({\rho }_m{c}_{pm}\right)$.

In simulations, some required parameters are: ${\delta }_x={\delta }_y={\delta }_t=1$ ($c=1$), $c_{pl}=c_{ps}=1$, $\widehat{\sigma }=1$, $Pr=50$, $H_v=5.9$, $d_p/L=0.0135$, $\lambda =\Gamma ={10}^3$, $k_f=0.0005$, $J=1$, and ${\varpi }_1={\varpi }_2=1/2$. The effective thermal conductivities $k_{ef}$ and $k_{em}$ are simply determined by $k_{ef}=\varphi k_f$ ($k_f=k_{fl}=k_{fs}$) and $k_{em}=\left(1-\varphi \right)k_m$, respectively. To eliminate the unphysical numerical diffusion, the relaxation parameters ${\mathsf{\zeta }}_3$ and ${\mathsf{\zeta }}_4$ (${\mathsf{\zeta }}_3={\mathsf{\zeta }}_4={\mathsf{\zeta }}_e$) related to the second-order moments are determined by ${\mathsf{\zeta }}_e=2-{\mathsf{\zeta }}_{\alpha }$ [20,27]. The isothermal CLB equation [32] in conjunction with the volumetric LB scheme [36] is employed to solve the flow filed, and to realize the velocity and thermal boundary conditions, the non-equilibrium extrapolation scheme [37] is adopted. Each C++ serial program runs on a desktop computer (Inter(R) Core(TM) i5-8400 [email protected]; RAM: 8.0 GB).

In Figure 2, the melting front at different values of $Fo$ for $Ra={10}^6$ and 10⁸ with $Da={10}^{-4}$, $\varphi =0.8$ and $St=1$ are presented. For comparison purpose, the other parameters are chosen as $F_{\varphi }=0.068$, ${\theta }_h=1$, ${\theta }_i={\theta }_c={\theta }_0=0$, and ${\theta }_{\mathrm{melt}}=0.3$ ($\theta =\left(T-T_c\right)/\Delta T$) [9]. For $Ra={10}^6$, $N_x\times N_y=100\times 100$ is employed, and for $Ra={10}^8$, $N_x\times N_y=300\times 300$ is employed. From the figure it can be found that the CLB results match well with the results obtained by FVM [9]. As show in Figure 2a, at $Ra={10}^6$, the shape of the melting front is almost planar since the heat transfer inside the cavity is controlled by conduction due to the large value of metal foam-to-PCM thermal conductivity ratio. As $Ra$ increases ($Ra={10}^8$), the effect of natural convection becomes stronger, then the melting front moves faster near the top wall (see Figure 2b).

In Figure 3, the temperature profiles (at $Y=0.5$) at different values of $Fo$ for $Ra={10}^6$ and ${10}^8$ with $Da={10}^{-4}$, $\varphi =0.8$ and $St=1$ are shown. As can be seen from the figure, at the very beginning ($Fo=0.00005$), the thermal non-equilibrium effect is apparent (the temperature difference is very high). As $Fo$ increases, the temperature difference progressively decreases due to the interfacial heat transfer. At $Fo=0.001$, the temperature difference tends to $0$ except in the mushy zone, indicating that the thermal non-equilibrium effect at this stage is weak. As clearly illustrated in Figure 3, the temperature difference has a maximum value near the phase-change region (mushy zone). The FVM results [9] are also shown in Figure 3 for comparisons. It can be found that the CLB results are in good agreement with those results in Ref. [9].

The flow field with the phase field at different values of $Fo$ for $Ra={10}^8$ with $Da={10}^{-4}$, $\varphi =0.8$ and $St=1$ are presented in Figure 4. As shown in the figure, an elliptical shape vortex appears in the liquid region at $Fo=0.0004$. As $Fo$ increases, natural convection effect becomes stronger, resulting in more hot liquid to move upwards, and consequently, the PCM melts faster in the upper region of the cavity. Although the Rayleigh number is high, the natural convection effect is weak since it has been severely suppressed by the effect caused by the metal foam’s ligament, and, during the melting process, the heat conduction through the high thermal conductive metal foam dominates the heat transfer.

Since the cascaded collision model has tunable relaxation parameters as well as a high degree of Galilean invariance, it is expected that the unphysical numerical diffusion can be effectively eliminated in numerical simulations by using the CLB model. To confirm this statement, the liquid fraction distributions calculated by the CLB and BGK-LB models at $Fo=0.002$ for $Ra={10}^8$ with $Da={10}^{-4}$, $\varphi =0.8$ and $St=1$ are presented in Figure 5. The BGK-LB result means that the PCM’s temperature field is solved by the enthalpy-based BGK-LB equation (when ${\mathbf{N}}_T$ is an identity matrix and ${\mathsf{\zeta }}_i=1/{\tau }_{f_l}$, the enthalpy-based CLB equation degrades into the enthalpy-based BGK-LB equation with ${\alpha }_{ef}=c_{sf}^2\left({\tau }_{f_l}-0.5\right){\delta }_t$). As can be seen in Figure 5a, the liquid fraction distribution calculated by the BGK-LB model exhibits significant numerical oscillations because the BGK-LB model has no free relaxation parameters to eliminate the effect of numerical diffusion on numerical simulations. On the contrary, the effect of numerical diffusion is almost invisible in Figure 5b because the CLB model has free relaxation parameters to eliminate such an effect. By setting ${\mathsf{\zeta }}_e=2-{\mathsf{\zeta }}_{\alpha }$, the effect of the numerical diffusion on numerical simulations can be effectively eliminated by the present CLB model.

$F_{\varphi }$ and $K$ are functions of the structure of the porous media. Equation (6) is proposed for packed beds of spherical particles based on Ergun’s experimental study, while Equation (7) is proposed for high porosity metal foams such as aluminum foam [33,34]. In what follows, the empirical correlations for $F_{\varphi }$ and $K$ given by Equations (6) and (7) are evaluated. In numerical simulations, the related parameters are set as ${\theta }_h=1$, ${\theta }_i={\theta }_c={\theta }_{\mathrm{melt}}={\theta }_0=0$. Numerical simulations are conducted based on $N_x\times N_y=100\times 100$. The total liquid fractions $f_{l,total}$ for different porosities with $Ra={10}^8$ and $St=1$ are presented in Figure 6. For a low value of porosity ($\varphi =0.8$), the total liquid fraction $f_{l,total}$ predicted with Equation (6) is almost the same as that predicted with Equation (7). As the porosity increases, the two empirical correlations yield different results (see Figure 6b–d). When $F_{\varphi }$ and $K$ are given by Equation (6), the melting rate of the PCM is faster than that predicted with Equation (7). This can be expected because the permeability given by Equation (6) is much larger than that given by Equation (7) at high porosities, leading to a stronger convection effect which promotes the heat transfer during melting process. Hence, for high porosity metal foams, appropriate empirical correlations should be employed. Since Equation (6) is proposed for packed beds, it overestimates the melting rate of the PCM. The empirical correlations given by Equation (7) are expected to provide a good estimate of the practical situations.

In Figure 7, the total liquid fractions of $Ra={10}^6$ and ${10}^8$ for different porosities with $F_{\varphi }$ and $K$ given by Equation (7) are shown. As can be seen in the figure, when $\varphi$ increases, the PCM’s melting rate decreases. An increase in the metal foam’s porosity leads to stronger convection effect of the liquid; however, it also reduces the effect of conduction through the metal foam’s ligament. This observation denotes that the total heat transfer from hot (left) wall is dominated by heat conduction through the high thermal conductive metal foam. The comparison between the total liquid fractions of $Ra={10}^6$ and ${10}^8$ indicates that the effect of $Ra$ on melting is weak because the convection has been severely suppressed by the metal foam’s ligament. Although smaller porosity increases the PCM’s melting rate, it also denotes smaller volume of PCM for LHS. A further optimization investigation on porosity of the metal foam to balance the PCM’s melting rate and LHS capacity is quite necessary.

5. Conclusions

An enthalpy-based CLB model is constructed for investigating heat transfer in solid-liquid PCMs with metal foams. The effects of the inertial coefficient, permeability and porosity on phase-change processes are examined. The findings of this research are:

(1)

The melting front and temperature profiles at different Fourier numbers predicted by the CLB model match well with the available data in previous studies, demonstrating the effectiveness and practicability of the CLB model for investigating heat transfer in solid-liquid PCMs with metal foams.

(2)

The empirical correlations of $F_{\varphi }$ and $K$ given by Equation (6) based on packed beds overestimate the PCM’s melting rate when the metal foam’s porosity is high, while the empirical correlations for metal foams such as aluminum foam given by Equation (7) are expected to provide a good estimate of the practical situations.

(3)

The PCM’s melting rate increases as the metal foam’s porosity decreases since the total heat transfer from the hot wall is dominated by heat conduction through high thermal conductive metal foam. Moreover, the effect of the Rayleigh number on phase-change process is weak since it has been severely suppressed by the metal foam’s ligament.

(4)

Although smaller porosity increases the PCM’s melting rate, it also reduces the volume of PCM for LHS, leading to a lower LHS capability. A further optimization investigation on metal foam’s porosity to balance PCM’s melting rate and LHS capacity is quite necessary.