Numerical Study on Heat-Transfer Characteristics of Convection Melting in Metal Foam under Sinusoidal Temperature Boundary Conditions

Abstract

Convection melting in metal foam under sinusoidal temperature boundary conditions is numerically studied in the present study. A multiple-relaxation-time lattice Boltzmann method, in conjunction with the enthalpy approach, is constructed to model the melting process without iteration steps. The effects of the porosity, phase deviation, and periodicity parameter on the heat-transfer characteristics are investigated. For the cases considered in this work, it is found that the effects of the phase deviation and periodicity parameter on the melting rate are weak, but the melting front can be significantly affected by the sinusoidal temperature boundary conditions.

1. Introduction

Latent heat storage (LHS), which uses solid–liquid phase-change materials (PCMs) as thermal-energy storage media, has been widely employed in industrial waste heat utilization to build energy saving systems, solar thermal utilization systems, etc. LHS with solid–liquid PCMs has become an important research topic during the past 30 years, and numerous reviews about this topic have been published. Zalba et al. [1] carried out a comprehensive review of the materials, the heat transfer process, and applications of LHS using solid–liquid PCMs. Farid et al. [2] reviewed the efforts in developing new PCMs for LHS applications. In a recent review by Nazir et al. [3], the applications of various PCMs, based on their thermophysical properties, were summarized, and the strategies for improving the characteristics of thermal-energy storage through nanomaterial additives, as well as encapsulation, were discussed in detail.

LHS with the use of solid–liquid PCMs has gradually become the preferred thermal-energy storage pattern, as solid–liquid PCMs have some outstanding features, such as the energy storage density being very high and the temperature fluctuation being small. However, the thermal conductivities for most of the solid–liquid PCMs are low (0.1~0.6 W/(m·K) [4]). This serious shortcoming strongly slows down the charging and discharging rates of thermal energy. To improve the LHS system’s thermal performance, three main kinds of enhancement approaches have been employed: improving the uniformity of heat-transfer process, enhancing the thermal conductivity of PCMs, and extending the heat-transfer surface [5]. Among these enhancement approaches, enhancing the thermal conductivity performance of PCMs is an efficient way to improve the LHS system’s thermal performance. High-porosity metal foams attract great attention for LHS applications because of their attractive advantages, such as high thermal conductivity and large specific surface areas.

In recent decades, numerous numerical studies on the characteristics of solid–liquid phase change in metal foams (porous media) have been performed. Weaver and Viskanta [6] numerically and experimentally investigated the melting process of ice in a cylindrical capsule filled with glass or aluminum beads. Beckermann and Viskanta [7] studied the melting and solidification processes of gallium in a square cavity filled with glass beads. They found that the shape of the interface can be considerably influenced by the convection effect in the liquid region. Tong et al. [8] performed a numerical study on the melting and freezing of a water–aluminum matrix system in a cylindrical annulus. They found that the heat-transfer rates of enhanced cases were increased by one order of magnitude, compared with that of the base case without an aluminum matrix.

In the numerical studies [6,7,8], the local thermal equilibrium (LTE) assumption is adopted, as the thermal conductivity of the solid matrix is low. However, for high-thermal-conductivity metal foams, such as copper or aluminum foam, the local thermal non-equilibrium (LTNE) effect (temperature difference) between a PCM and a metal matrix during the melting process should be considered. Harris et al. [9] developed an approximate theoretical enthalpy model (LTNE model) in which a temperature difference between the PCM and the walls of the pores was maintained. Based on the approximate model, the conditions for the occurrence of LTE were analyzed. Mesalhy et al. [10] developed a two-temperature model to analyze the LTNE effect between the PCM and the metal matrix, and a parametric study was performed to investigate the effects of thermal conductivity and porosity. Krishnan et al. [11] also proposed an LTNE model for simulating convection melting in metal foams, and the merits of using metal foam for enhancing thermal storage systems’ effective thermal conductivity were discussed. An LTNE model regarding the volume change of the PCM was proposed by Yang and Garimella [12], and the effects of volume expansion/shrinkage were analyzed. Li et al. [13] investigated the melting of paraffin embedded in open-cell copper foam, and the effects of the morphology parameters of the metal foam on the temperature distributions were investigated. Zhao et al. [14] performed a numerical investigation on melting and solidification in copper foam, and the kinetic undercooling of solidification was analyzed. Wang et al. [15] studied the pore-scale melting in metal foams; different metal foams combined with paraffin and other PCMs were investigated to obtain the composite materials’ effective thermal conductivity.

The above literature review indicates that many numerical studies have been carried out on the heat-transfer performance of PCMs in metal foams based on an LTNE model. Moreover, our literature survey with respect to improving an LHS system’s thermal performance using metal foams found that nearly all of the numerical studies were conducted under constant wall heat flux or constant wall temperatures (uniform thermal boundary conditions). A fundamental understanding of the heat-transfer characteristics of melting in metal foams under non-uniform thermal boundary conditions is still lacking, and more studies are required. For natural convection in enclosures, previous studies indicated that non-uniform thermal boundary conditions (e.g., sinusoidal temperature boundary conditions) can significantly affect the flow structures and heat-transfer characteristics [16,17,18]. As expected, new heat-transfer characteristics can be created in the solid–liquid phase change of PCMs under sinusoidal temperature boundary conditions. Hence, this work aimed to study the heat-transfer characteristics of convection melting in metal foam under sinusoidal temperature boundary conditions. A multiple-relaxation-time (MRT) lattice Boltzmann (LB) method, in conjunction with the enthalpy approach, was constructed to model the melting process without iteration steps. This work will help in providing a valuable reference for improving the thermal performance of LHS systems.

2. Model Description

2.1. Physical Model

The problem considered in this work is shown in Figure 1. Initially, the temperatures of the PCM and metal foam are equal to $T_i$ ($T_i<T_{\mathrm{melt}}$). At $t=0$, a sinusoidally varying temperature $T=T_h+\Delta T\sin \left(2k\pi y/L+\phi \right)$ ($T_h>T_{\mathrm{melt}}$) is imposed on the left wall, and then the PCM begins to melt. Note that the average temperature of the left wall is $T_h$, $\Delta T=T_h-T_{\mathrm{melt}}$ is the characteristics temperature, $k$ is the periodicity parameter, and $\phi$ is the phase deviation (phase of the sinusoidal profile).

2.2. Governing Equations

For convection melting of solid–liquid PCMs embedded in metal foams, the following assumptions are made: (1) metal foam ($m$) is homogeneous and the pore diameter is uniform; (2) the flow (liquid region) is incompressible and laminar; (3) the volume change is neglected, i.e., ${\rho }_f = {\rho }_l= {\rho }_s$ (the subscript $f$ denotes PCM, $l$ denotes liquid PCM and $s$ denotes solid PCM). Based on the LTNE model, the governing equations are provided by [11,19,20,21]

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

(1)

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

(2)

$$\frac{\partial }{\partial t}\left(\varphi {\rho }_f{c}_{pf}{T}_f\right)+\nabla \cdot \left({\rho }_l{c}_{pl }{T}_f\mathbf{u}\right)=\nabla \cdot \left(k_{ef}\nabla T_f\right)+h_{mf}{a}_{mf}\left(T_m-T_f\right)-\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]=\nabla \cdot \left(k_{em}\nabla T_m\right)+h_{mf}{a}_{mf}\left(T_f-T_m\right)$$

(4)

where u and $p$ are the velocity and pressure in the liquid region, respectively; ${\rho }_f$ is the density; $T$ is the temperature; $\varphi$ is the metal foam’s porosity; $c_p$ is the specific heat; $v_e$ is the effective kinematic viscosity; $k_e$ is the effective thermal conductivity; $h_{mf}$ is the interfacial heat-transfer coefficient; $a_{mf}$ is the specific surface area of metal matrix; $f_l$ is the liquid fraction; and $L_a$ is the PCM’s latent heat.

The total body force $\mathbf{F}$ is determined by [22,23]

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

(5)

where $v$ is the liquid PCM’s kinematic viscosity and $K$ and $F_{\varphi }$ are the metal foam’s permeability and inertia coefficient, respectively. $\mathbf{G}$ is the buoyancy force approximated by

$$\mathbf{G}=-\mathbf{g}\beta \left(T_f-T_0\right)f_l$$

(6)

where $\mathbf{g}$ is the gravitational acceleration, $T_0$ is the reference temperature, and $\beta$ is the thermal-expansion coefficient.

For metal foams (e.g., aluminum or copper foam), the correlations of $F_{\varphi }$ and $K$ can be found in [24,25]. The effective thermal conductivities $k_{ef}$ and $k_{em}$ can be determined by analytical models [25,26]. In previous studies, the correlation for convection heat transfer through a bank of staggered cylinders proposed by Churchill and Chu [27] was widely employed to determine $h_{mf}$. The empirical formula for $a_{mf}$ can be found in [25]. To determine the temperature-dependent (thermodynamic or dynamic mechanical) properties of complex materials, such as the cross-linking of polymers, the methods proposed by Likozar and Krajnc [28,29,30] can be employed.

3. Numerical Method

As a mesoscopic approach evolved from the lattice-gas automata [31], the LB method [32,33,34] has become an efficient numerical methodology for modeling solid–liquid phase-change problems [35,36]. In this section, the MRT-LB method, in conjunction with the enthalpy approach, is introduced to model the melting process without iteration steps.

3.1. MRT-LB Equation for Flow Field

For the 2D problem considered in this work, the D2Q9 lattice is employed [34]

$${\mathbf{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 \\ \left(\cos \left[\left(2i-9\right)\pi /4\right],\sin \left[\left(2i-9\right)\pi /4\right]\right)\sqrt{2}c,\hfill & i=5~8\hfill \end{array}$$

(7)

where $c={\delta }_x/{\delta }_t$ is the lattice speed (${\delta }_x$ is the lattice step and ${\delta }_t$ is the time step). In this work, $c$ is set to $1$ (${\delta }_x={\delta }_t$).

The MRT-LB equation for the flow field can be written as [37,38,39]

$$f_i\left(\mathbf{x}+{\mathbf{e}}_i{\delta }_t, t+{\delta }_t\right)=f_i\left(\mathbf{x}, t\right)-{\overline{\mathsf{\Lambda }}}_{ij}\left(f_j-f_j^{eq}\right)|{}_{\left(x, t\right)}+{\delta }_t\left({\tilde{S}}_i-0.5{\overline{\mathsf{\Lambda }}}_{ij}{\tilde{S}}_j\right)$$

(8)

where $f_i\left(\mathbf{x}, t\right)$ is the density distribution function, $f_i^{eq}\left(\mathbf{x}, t\right)$ is the equilibrium value of $f_i\left(\mathbf{x}, t\right)$, $\overline{\mathsf{\Lambda }}=[{\overline{\mathsf{\Lambda }}}_{ij}]$ is the collision matrix, and ${\tilde{S}}_i$ is the forcing term.

The MRT-LB Equation (8) can be divided into two parts: a collision part and a streaming part. By multiplying a transformation matrix $\mathsf{M}$, the collision part can be carried out in moment space as

$${\mathbf{m}}^{*}\left(\mathbf{x}, t\right)=\mathbf{m}\left(\mathbf{x}, t\right)-\mathsf{\Lambda }\left(\mathbf{m}-{\mathbf{m}}^{eq}\right)|{}_{\left(x, t\right)}+{\delta }_t\left(\mathbf{I}-\frac{\mathsf{\Lambda }}{2}\right)\mathbf{S}$$

(9)

The streaming part is performed in velocity space as

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

(10)

where $\mathsf{\Lambda }$ is the relaxation matrix ($\mathsf{\Lambda }=\mathbf{M}\overline{\mathsf{\Lambda }}{\mathbf{M}}^{-1}=\mathrm{diag}\left(1, 1, 1,s_e, s_v, s_v, s_q, s_q,s_{\epsilon }\right)$), $\mathbf{m}=|m\rangle =\mathbf{M}\mathbf{f}$, ${\mathbf{m}}^{eq}=|m_{}^{eq}\rangle =\mathbf{M}{\mathbf{f}}^{eq}$, $\mathbf{S}=|S\rangle =\mathbf{M}\tilde{\mathbf{S}}$, in which $\mathbf{f}=|f\rangle$, ${\mathbf{f}}^{eq}=|f_{}^{eq}\rangle$, and $\tilde{S}=|\tilde{S}\rangle$. Here, Dirac notation $|\cdot \rangle$ denotes a nine-dimensional column vector, e.g., $|m\rangle ={\left(m_0,m_1,\dots , m_8\right)}^T$. $f_i^{*}$ is determined by ${\mathbf{f}}^{*}=|f_{}^{*}\rangle ={\mathbf{M}}^{-1}{\mathbf{m}}^{*}$.

$\mathbf{M}$ is a non-orthogonal transformation matrix [39]

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

(11)

The equilibrium moments $\left\{m_i^{eq}\right\}$ are determined by

$$\begin{gathered} m_0^{eq}={\rho }_f, m_1^{eq}={\rho }_f{u}_x, m_2^{eq}={\rho }_f{u}_y, m_3^{eq}=\frac{2}{3}{\rho }_f+\frac{{\rho }_f\left(u_x^2+u_y^2\right)}{\varphi }, m_4^{eq}=\frac{{\rho }_f\left(u_x^2-u_y^2\right)}{\varphi } \\ {m}_5^{eq}={\rho }_f{u}_x^{}{u}_y, m_6^{eq}=\frac{1}{3}{\rho }_f{u}_y, m_7^{eq}=\frac{1}{3}{\rho }_f{u}_x, m_8^{eq}=\frac{1}{9}{\rho }_f+\frac{1}{3}\frac{{\rho }_f\left(u_x^2+u_y^2\right)}{\varphi } \end{gathered}$$

(12)

The source terms $\left\{S_i^{}\right\}$ are determined by

$$\begin{gathered} S_0=0, S_1={\rho }_f{F}_x, S_2={\rho }_f{F}_y, S_3=\frac{2{\rho }_f\left(u_x{F}_x+u_y{F}_y\right)}{\varphi }, S_4=\frac{2{\rho }_f\left(u_x{F}_x-u_y{F}_y\right)}{\varphi } \\ {S}_5=\frac{{\rho }_f\left(u_x{F}_y+u_y{F}_x\right)}{\varphi }, S_6=\frac{1}{3}{\rho }_f{F}_y, S_7=\frac{1}{3}{\rho }_f{F}_x, S_8=\frac{2}{3}\frac{{\rho }_f\left(u_x{F}_x+u_y{F}_y\right)}{\varphi } \end{gathered}$$

(13)

To implement the non-slip velocity boundary condition on the phase interface accurately, the volumetric LB scheme [40] is employed; then, a new density distribution function is defined:

$$f_i^{+}=f_l{f}_i+\left(1-f_l\right)f_i^{eq}\left({\rho }_f,{\mathbf{u}}_s\right)$$

(14)

In Equation (14), the superscript “$+$” denotes that the solid-phase effect has been considered, and ${\mathbf{u}}_s =0$ (the solid phase is static). Accordingly, ${\rho }_f$ and $\mathbf{u}$ are defined as

$${\rho }_f={\displaystyle \sum _{i=0}^8{f}_i^{}}$$

(15)

$${\rho }_f\mathbf{u}={\displaystyle \sum _{i=0}^8{e}_i{f}_i^{+}}+\frac{{\delta }_t}{2}{\rho }_f\mathbf{F}$$

(16)

$p$ is defined as $p=\rho c_s^2/\varphi$ ($c_s=1/\sqrt{3}$ is the sound speed). Explicitly, $\mathbf{u}$ can be calculated via [41]

$$\mathbf{u}=\frac{\mathbf{v}}{l_0+\sqrt{l_0^2+l_1\left|\mathbf{v}\right|}}$$

(17)

where

$${\rho }_f\mathbf{v}={\displaystyle \sum _{i=0}^8{\mathbf{e}}_i{f}_i^{+}}+\frac{{\delta }_t}{2}{\rho }_f\varphi \mathbf{G}$$

(18)

$$l_0=\frac{1}{2}\left(1+\varphi \frac{{\delta }_t}{2}\frac{v}{K}\right), l_1=\varphi \frac{{\delta }_t}{2}\frac{F_{\varphi }}{\sqrt{K}}$$

(19)

The kinetic viscosity $v=c_s^2\left(s_v^{-1}-0.5\right){\delta }_t$ and the bulk viscosity $\xi =c_s^2\left(s_e^{-1}-0.5\right){\delta }_t$.

3.2. MRT-LB Equation for the Temperature Field of the PCM

Equation (3) can be rewritten as

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

(20)

where $H_f=\sigma c_{pl }{T}_f+L_a{f}_l$ is the effective 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. When $f_l=1$ (liquid region), $H_f=c_{pl }{T}_f+L_a$ and ${\sigma }_l=1$; when $f_l=0$ (solid region), $H_f={\sigma }_s{c}_{pl }{T}_f$ and ${\sigma }_s=\frac{{\rho }_s{c}_{ps}}{{\rho }_l{c}_{pl }}$.

For the temperature field of the PCM, governed by Equation (20), the D2Q5 lattice is adopted and $\left\{{\mathbf{e}}_i|i=0, \dots , 4\right\}$ are provided in Equation (7). The enthalpy-based MRT-LB equation is determined by

$$\mathbf{g}\left(\mathbf{x}+\mathbf{e}{\delta }_t, t+{\delta }_t\right)=\mathbf{g}\left(\mathbf{x}, t\right)-{\mathbf{N}}^{-1}\mathsf{\Theta }\left({\mathbf{n}}_g-{\mathbf{n}}_g^{eq}\right)|{}_{\left(\mathbf{x}, t\right)}+{\delta }_t{\mathbf{N}}^{-1}{\mathbf{S}}_{\mathrm{PCM}}$$

(21)

where $g_i$ is the enthalpy distribution function and $\mathsf{\Theta }=\mathrm{diag}\left(1, {\zeta }_{\alpha },{\zeta }_{\alpha },{\zeta }_e,{\zeta }_e\right)$ is the relaxation matrix.

Through the transformation matrix $\mathbf{N}$, the collision part of the MRT-LB Equation (21) is carried out in moment space as

$${\mathbf{n}}_g^{*}\left(\mathbf{x}, t\right)={\mathbf{n}}_g\left(\mathbf{x}, t\right)-\mathsf{\Theta }\left({\mathbf{n}}_g-{\mathbf{n}}_g^{eq}\right)|{}_{\left(\mathbf{x}, t\right)}+{\delta }_t{\mathbf{S}}_{\mathrm{PCM}}$$

(22)

The streaming part is performed in velocity space as

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

(23)

where ${\mathbf{n}}_g =\mathbf{N}\mathbf{g}$ is the moment, and ${\mathbf{n}}_g^{eq} =\mathbf{N}{\mathbf{g}}^{eq}$ is the corresponding equilibrium moment. Here, $g_i^{eq}$ is the equilibrium value of $g_i$, and ${\mathbf{g}}^{*} ={\mathbf{N}}^{-1}{\mathbf{n}}_g^{*}$.

$\mathbf{N}$ is a non-orthogonal transformation matrix [39]

$$\mathbf{N}=\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 1& \hfill 1& \hfill 1\\ \hfill 0& \hfill 1& \hfill -1& \hfill 1& \hfill -1\end{array}\right]$$

(24)

The equilibrium moment ${\mathbf{n}}_g^{eq}$ is

$${\mathbf{n}}_g^{eq}={\left(H_f,\frac{c_{pl }{T}_f{u}_x}{\varphi },\frac{c_{pl }{T}_f{u}_y}{\varphi },{\varpi }_1{c}_{pl }{T}_f ,0\right)}^T$$

(25)

where ${\varpi }_1\in \left(0, 1\right)$. Correspondingly, $g_i^{eq}$ is given by

$$g_i^{eq}=\{\begin{array}{ll}{H}_f-{\varpi }_1{c}_{pl }{T}_f,\hfill & i=0\hfill \\ \frac{1}{4}{\varpi }_1{c}_{pl }{T}_f\left(1+\frac{e_i\cdot \mathbf{u}}{c_{sf}^2\varphi }\right),\hfill & i=1~4\hfill \end{array}$$

(26)

where $c_{sf}^{}=\sqrt{{\varpi }_1/2}$ is the sound speed.

The source term ${\mathbf{S}}_{\mathrm{PCM}}$ is chosen as

$${\mathbf{S}}_{\mathrm{PCM}}=S_{\mathrm{PCM}}{\left(1,0, 0,0,0\right)}^T$$

(27)

where $S_{\mathrm{PCM}}=S{r}_f +\frac{1}{2}{\delta }_t{\partial }_{t}S{r}_f$ and $S{r}_f=h_{mf}{a}_{mf}\left(T_m-T_f\right)/\left(\varphi {\rho }_l\right)$.

$H_f$ is defined as

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

(28)

$T_f$ can be determined via the following equation:

$$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}$$

(29)

where $T_{fs}$ is solidus temperature and $T_{fl}$ is liquidus temperature ($T_{fs}\le T_{fl}$); $H_{fs}$ ($H_{fl}$) is the effective enthalpy corresponding to $T_{fs}$ ($T_{fl}$).

$f_l$ is determined by

$$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}$$

(30)

${\alpha }_{ef}$ is defined as

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

(31)

3.3. MRT-LB Equation for the Temperature Field of Metal Foam

Equation (4) can be rewritten as

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

(32)

For the temperature field of metal foam, governed by Equation (32), the MRT-LB equation based on D2Q5 lattice is as follows:

$$\mathbf{h}\left(\mathbf{x}+\mathbf{e}{\delta }_t, t+{\delta }_t\right)=\mathbf{h}\left(\mathbf{x}, t\right)-{\mathbf{N}}^{-1}\mathbf{Q}\left({\mathbf{n}}_h-{\mathbf{n}}_h^{eq}\right)|{}_{\left(x, t\right)}+{\delta }_t{\mathbf{N}}^{-1}{\mathbf{S}}_{\mathrm{metal}}$$

(33)

where $h_i\left(\mathbf{x}, t\right)$ is the temperature distribution function, $\mathbf{Q}=\mathrm{diag}\left(1,{\eta }_{\alpha }, {\eta }_{\alpha }, {\eta }_e, {\eta }_e\right)$ is the relaxation matrix, and $\mathbf{N}$ is given by Equation (24).

The collision part of the MRT-LB Equation (33) is performed in moment space as

$${\mathbf{n}}_h^{*}\left(\mathbf{x}, t\right)={\mathbf{n}}_h\left(\mathbf{x}, t\right)-\mathbf{Q}\left({\mathbf{n}}_h-{\mathbf{n}}_h^{eq}\right)|{}_{\left(x, t\right)}+{\delta }_t{\mathbf{S}}_{\mathrm{metal}}$$

(34)

where ${\mathbf{n}}_h =\mathbf{N}\mathbf{h}$ is the moment, and ${\mathbf{n}}_h^{eq} =\mathbf{N}{\mathbf{h}}^{eq}$ is the corresponding equilibrium moment. Here, $h_i^{eq}$ is the equilibrium value of $h_i$. The streaming step is carried out in the velocity space as follows:

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

(35)

where ${\mathbf{h}}^{*} ={\mathbf{N}}^{-1}{\mathbf{n}}_h^{*}$. The equilibrium moment ${\mathbf{n}}_h^{eq}$ is defined as

$${\mathbf{n}}_h^{eq}={\left(c_{pm}{T}_m,0,0,{\varpi }_2{c}_{pm}{T}_m,0\right)}^T$$

(36)

where ${\varpi }_2\in (0,1)$. $h_i^{eq}$ is determined by

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

(37)

The source term ${\mathbf{S}}_{\mathrm{metal}}$ is chosen as

$${\mathbf{S}}_{\mathrm{metal}}=S_{\mathrm{metal}}{\left(1,0, 0,0,0\right)}^T$$

(38)

where $S_{\mathrm{metal}}=S{r}_m +\frac{1}{2}{\delta }_t{\partial }_{t}S{r}_m$ and $S{r}_m=h_{mf}{a}_{mf}\left(T_f-T_m\right)/\left[\left(1-\varphi \right){\rho }_m\right]$.

$T_m$ is defined by

$$T_m=\frac{1}{c_{pm}}{\displaystyle \sum _{i=0}^4{h}_i}$$

(39)

${\alpha }_{em}$ is given by

$${\alpha }_{em}=\frac{k_{em}}{\left(1-\varphi \right){\rho }_m{c}_{pm}}=c_{sm}^2\left({\eta }_{\alpha }^{-1}-\frac{1}{2}\right){\delta }_t$$

(40)

where $c_{sm}=\sqrt{{\varpi }_2/2}$ is the sound speed.

4. Numerical Results

In this section, numerical simulations are performed to investigate the effects of the porosity, the phase deviation, and the periodicity parameter on the heat-transfer performance of the convection melting of solid–liquid PCM embedded in metal foam under sinusoidal temperature boundary conditions. The characteristic parameters include $Pr=v_{fl}/{\alpha }_{fl}$ (Prandtl number), $Ra=g\beta \Delta T{L}^3/\left(v_{fl}{\alpha }_{fl}\right)$ (Rayleigh number), $Da=K/L^2$ (Darcy number), $J=v_e/v_{fl}$ (viscosity ratio), $\lambda =k_m/k_{fl}$ (thermal conductivity ratio), $\mathsf{\Gamma }={\alpha }_m/{\alpha }_{fl}$ (thermal diffusivity ratio), $\widehat{\sigma }={\rho }_m{c}_{pm}/\left({\rho }_l{c}_{pl}\right)$ (metal foam to liquid PCM heat capacity ratio), $H_v=h_{mf}{a}_{mf}{d}_p^2/k_f$ (volumetric heat transfer coefficient), $Fo=t{\alpha }_{fl}/L^2$ (Fourier number), and $St=c_{pl}\Delta T/L_a$ (Stefan number), where ${\alpha }_f=k_f/{\left(\rho c_{p }\right)}_f$ is thermal diffusivity of PCM, ${\alpha }_m=k_m/{\left(\rho c_{p }\right)}_m$ is thermal diffusivity of metal foam, and $d_p^{}$ is mean pore diameter.

In simulations, the required parameters are chosen as follows: $Pr=50$, $F_{\varphi }=0.068$, $Da={10}^{-4}$, $St=1$, ${\delta }_x={\delta }_y={\delta }_t=1$ ($c=1$), $c_{pl}=c_{ps} =1$, $J=\widehat{\sigma }=1$, $H_v=5.9$, $\lambda =\mathsf{\Gamma }={10}^3$, $d_p/L=0.0135$, $k_f=0.0005$ and ${\varpi }_1={\varpi }_2=1/2$. The relaxation rate ${\zeta }_e$ is determined by ${\mathsf{\zeta }}_e =2-{\mathsf{\zeta }}_{\alpha }$ to reduce the unphysical numerical diffusion [21,42]. The non-equilibrium extrapolation scheme [43] is adopted to realize the velocity and thermal boundary conditions. Numerical simulations are performed based on a grid size of $N_x\times N_y=150\times 150$. First, comparisons between the results predicted by the finite-volume method (FVM) [11] and the present method are made to validate the reliability of the present method. The predicted results are shown in Figure 2 and Figure 3, where the melting front (solid–liquid interface) and temperature profiles ($\theta =\left(T-T_{melt}\right)/\Delta T$) at different times for $Ra={10}^6$ and ${10}^8$ with $\varphi =0.8$ are presented. In the figures, it can be seen that the present results match well with the results in [11]. In what follows, the effects of the porosity, the phase deviation, and the periodicity parameter on the heat-transfer performance are investigated.

4.1. Effects of the Porosity and Phase Deviation

In this subsection, the effects of the porosity and the phase deviation are investigated. In Figure 4, the total liquid fractions for different $\varphi$ with $Ra={10}^6$, $\phi =\pi /4$ and $k=0$ are shown. As can be seen in Figure 4, the melting rate decreases as $\varphi$ increases. When $\varphi =0.8$, the completely melting time $Fo=0.00485$. As $\varphi$ increases to 0.9 and 0.95, the completely melting time $Fo$ augments to 0.0103 and 0.0209, respectively. When $\varphi$ increases from 0.8 to 0.9, the completely melting time increases by 112.37%; when $\varphi$ increases from 0.9 to 0.95, the completely melting time increases by 102.91%. The influence of the porosity on the melting rate is induced by two factors: one is that the mass of the metal foam decreases as the porosity increases, which reduces the effective thermal conductivity, and consequently, the performance of heat transfer is deteriorated; the other is that the mass of the PCM increases as the porosity increases, which results in melting-time augmentation.

In Figure 5, the total liquid fractions for $\varphi =0.8$ and $0.95$ under non-uniform ($T=T_h+\Delta T\sin \left(2\pi y/L+\pi /4\right)$ at $t=0$) and uniform ($T=T_h$ at $t=0$) thermal boundary conditions are presented. One can observe that the melting rate of the uniform case is only a little faster than that of the non-uniform case with the given parameters, as the characteristics temperatures are equal for the cases considered. In Figure 6, the total liquid fractions for different $\phi$ with $Ra={10}^6$ and $k=0$ under non-uniform thermal boundary conditions are shown. It can be observed that the melting rate increases as $\phi$ increases from 0 to $\pi /2$. As shown in Figure 5 and Figure 6, it seems that the effects of the phase deviation on the total liquid fraction are not very strong. This is because the average temperature of the left wall (with sinusoidally varying temperature) equals a constant.

In Figure 7 and Figure 8, the liquid-fraction fields for different values of $\phi$ with $\varphi =0.8$ and $0.95$ under non-uniform thermal boundary conditions are shown. As mentioned above, the effects of the phase deviation on the total liquid fraction are weak. However, in Figure 7 and Figure 8 it can be clearly observed that the melting process can be significantly affected by the phase deviation. For the cases under uniform heating (Figure 7a and Figure 8a), the melting front is almost parallel to the vertical walls, as the conduction effect dominates the heat-transfer process. For the cases under sinusoidal temperature boundary conditions, the phase interface is in a bending shape. This is because under the non-uniform thermal boundary condition, the convection effect in the related region is much stronger than that in the rest of the region. As shown in the figures, for $0<\phi <\pi /2$, the convective effect near the bottom wall is stronger and the melting front moves faster near the bottom wall. Obviously, this feature is rather valuable for practical LHS applications, as it offers a possible tool for controlling the melting front.

4.2. Effects of the Periodicity Parameter

In this subsection, the effects of the periodicity parameter $k$ on the performance of heat transfer are studied. In Figure 9, the total liquid fractions for different values of the periodicity parameter $k$ with $Ra={10}^6$, $\varphi =0.9$ and $\phi =0$ are shown. As presented in the figure, the effects of the periodicity parameter on the total liquid fraction are weak. As $k$ increases, the melting rate slightly increases, and approaches that of the uniform heating case. The liquid-fraction fields for different $k$ at $Fo=0.002$ under non-uniform thermal boundary conditions are presented in Figure 10, and one can observe that the melting front can also be affected by the periodicity parameter. As $k$ increases to 4, the melting front is almost parallel to the vertical walls, which is similar to the situation of the uniform case.

5. Conclusions

An MRT-LB method in conjunction with the enthalpy approach was constructed for simulating convection melting in metal foam under sinusoidal temperature boundary conditions. The effects of the porosity, the phase deviation, and the periodicity parameter on the heat-transfer characteristics were investigated. The main conclusions are listed as follows:

(1)

The melting rate decreases as $\varphi$ increases. The influence of the porosity on the melting rate is induced by two factors: one is that the mass of the metal foam decreases as the porosity increases, which reduces the effective thermal conductivity; the other is that the mass of the PCM increases as the porosity increases, which results in melting time augmentation.

(2)

The melting rate increases as the phase deviation increases from 0 to $\pi /2$. Although the effects of the phase deviation on the melting rate (total liquid fraction) are weak, the melting front can be significantly affected by the phase deviation.

(3)

The effects of the periodicity parameter on the total liquid fraction are weak. However, the melting process can also be affected by the periodicity parameter. The above results provide a valuable reference for practical applications of LHS systems.