Insight into Foam Pore Effect on Phase Change Process in a Plane Channel under Forced Convection Using the Thermal Lattice Boltzmann Method

Abstract

In this work, the two-dimensional laminar flow and the heat transfer in an open-ended rectangular porous channel (metal foam) including a phase change material (PCM; paraffin) under forced convection were numerically investigated. To gain further insight into the foam pore effect on charging/discharging processes, the Darcy–Brinkmann–Forchheimer (DBF) unsteady flow model and that with two temperature equations based on the local thermal non-equilibrium (LTNE) were solved at the representative elementary volume (REV) scale. The enthalpy-based thermal lattice Boltzmann method (TLBM) with triple distribution function (TDF) was employed at the REV scale to perform simulations for different porosities ($0.7\le \mathsf{\epsilon }\le 0.9$) and pore per inch (PPI) density ($10\le \mathrm{P}\mathrm{P}\mathrm{I}\le 60$) at Reynolds numbers (Re) of $200$ and $400$. It turned out that increasing Re with high porosity and PPI (0.9 and 60) speeds up the melting process, while, at low PPI and porosity (10 and 0.7), the complete melting time increases. In addition, during the charging process, increasing the PPI with a small porosity (0.7) weakens the forced convection in the first two-thirds of the channel. However, the increase in PPI with large porosity and high Re number limits the forced convection while improving the heat transfer. To sum up, the study findings clearly evidence the foam pore effect on the phase change process under unsteady forced convection in a PCM-saturated porous channel under local thermal non-equilibrium (LTNE).

1. Introduction

More attention is being paid to latent heat thermal energy storage systems (LHTESS) due to their large storage capacity in a small volume and low temperature fluctuation [1]. This allows removing temperature peaks due to the intermittent high power demand and storing a large amount of solar energy without significantly increasing the system temperature. Furthermore, LHTESS using phase change materials (PCMs) are widely used in various engineering applications such as building temperature regulation, waste heat recovery, compact heat exchangers, solar energy storage, concentrating solar plants, etc. [2,3]. However, their low thermal conductivity is one of their weaknesses because it restricts thermal transport by giving rise to slow heat dissipation during charging/discharging periods. Therefore, infiltrated metallic foams with PCMs are now increasingly used as porous media to overcome this defect due to their high thermal conductivity, their large specific surface area, and their low cost [4]. Indeed, this infiltration approach recently attracted more attention due to the high thermal conductivity of metallic struts which further spreads heat faster through the PCM matrix.

In terms of research dealing with the effectiveness of these foams, Zhao et al. [5] experimentally demonstrated that the metallic foam addition improves heat transfer via phase change by up to 10 times more than pure paraffin. Chen et al. [6] found that incorporating metallic foam in the heat transfer fluid (water) including the PCM raises the energy storage performance during charging/discharging processes while shortening the time of the 84.9% cycle compared to pure PCM. Cui [7] incorporated copper foam with PCM (paraffin) and found that the heat transfer rate increased (by 36%) more than that without the foam while the charging time decreased. Yang et al. [8] established physical and numerical models for a novel LHTESS using two gradient kinds (positive and negative) to simulate the melting front evolvement and thermal and dynamic fields during the melting process. They indicated that an optimized positive porosity gradient (0.89–0.95–0.98) is better suited to further decrease the melting time than a uniform gradient. Sardari et al. [9] investigated the porosity and pore size effects of metal foams on the PCM melting process. They stated that inserting porous foam with low porosity provides high performance and decreases the melting time by 85% compared to PCM alone. Furthermore, Lafdi et al. [10] experimentally studied the effects of the same parameters and found that a larger pore size and a greater porosity lead to the steady state being reached more rapidly and the heat transfer being intensified. Zhao et al. [11] showed that lower porosity speeds up PCM melting and solidification, and that a low pore density reduces the melting rate with a higher Rayleigh number (Ra). Zhu et al. [12] pointed out the porosity and pore density effects on thermal performance during PCM melting in a metallic fin foam. They stated that increasing the pore per inch (PPI) restricts convective heat transfer while dropping the energy storage rate. Yang et al. [13] deemed the effects of gradient in porosity, pore density, and matrix conductivity of water-saturated cell foams on the solidification process. They found that the foam property related to the gradient significantly affects the solidification rate and the complete solidification time. More explicitly, a positive porosity gradient accelerates solidification compared to cases of negative gradient or uniform porosity. Tao et al. [4] studied the influence of PPI and porosity on the melting process of a composite phase change material (metallic foam/paraffin) in natural convection. They showed that, for greater porosity (=0.98) and when the PPI decreases, natural convection dominates strongly, thereby increasing the PCM melting rate. A year earlier, Tao et al. [14] studied the effects of natural convection on latent heat storage performance. They noted that natural convection improves the PCM heat transfer and gives rise to a non-uniform evolvement of the temperature field and phase change interface during the melting process. Longeon et al. [15] explored the performance of an annular latent storage unit under natural convection. They concluded that, during the melting period, natural convection has a more meaningful effect on the heat transfer rate than that during solidification mode. Wang et al. [16] examined the influence of natural convection on PCM melting in a sleeve tube with internal fins. They argued that natural convection greatly affected melting mode upon considering a small fin ratio and angle between neighboring fins.

From the literature review conducted above, it can be noted that the number of studies (experimental and/or numerical) related to the storage of thermal energy at latent heat is still limited and deserves to be further investigated. Moreover, most of these studies were based on classical numerical approaches while being subject to various factors and time-consuming. However, heat transfer problems caused by forced convection in porous media filled with PCM were hardly addressed. Indeed, modeling and simulating such a problem remain arduous tasks because of the porous metallic medium complexity and the phase change process involved. Initially, research on these problems invoked the local thermal equilibrium (LTE) assumption. However, it turned out that it is rather the local thermal non-equilibrium (LTNE) hypothesis which should be considered, thereby requiring the use of the so-called LTNE models (two-equation energy models), one for each phase (solid and liquid), which can spot LTE zones, if any, in forced convection. Such models are advised in the literature because of their ability to deal with these issues, in particular, at the representative elementary volume (REV) scale.

Generally, lattice Boltzmann (LB) simulations in porous media involve the microscopic (pore scale) and mesoscopic (REV) scales. The first approach requires detailed knowledge of the porous medium structure, while the second settles for a simplified description avoiding the need to deal with large computation domains at the pore scale. Increasingly, the REV–LB model, as a computationally efficient numerical method [17,18], is being successfully used to predict flows with or without phase change in porous media. However, to the best of our knowledge, work on the application of this method (REV–LB) to problems under forced convection in porous media saturated with a phase change material remains scarce. Therefore, the present work leans on such an approach to scrutinize the effects of metallic foam pores on the phase change process in a rectangular channel under forced convection at the REV scale. Here, this method allows handling the porous medium’s existence and the transport phenomena involved using an additional term which translates the effects of drag forces, buoyancy, etc. [19].

In this context, Wu et al. [20] employed a new multiple relaxation time (MRT)–LB method to examine the transient fusion process in porous media at the REV scale. Their findings demonstrated the capability and reliability of the method for handling solid–liquid phase change issues in porous media. Gao and Chen [21] simulated freezing and solidification behavior in porous media under natural convection using an enthalpy-based LBM at REV scale. From the results obtained over an ample range of dimensionless parameters, they stated that the model can handle freezing and solidification processes in porous media. Liu and He [22] succeeded in simulating conduction fusion in a semi-infinite space, solidification in a semi-infinite corner, and convection fusion in a square cavity filled with porous media using a double MRT–LBM. Most recently, Mabrouk et al. [23] analyzed the heat transfer under forced convection through a PCM-saturated porous channel under LTNE during the melting and solidification cycles using an LBM with a single relaxation time (SRT) at the REV scale and relevant parameters such as the Reynolds number, Eckert number, and porosity. They showed that the storage energy is optimal for a critical Reynolds number in the range examined. At the REV scale, Jourabian et al. [24] also numerically studied the fusion cycle in a porous rectangular cavity including two hot cylinders using an enthalpy-based LBM with double distribution functions (DDF). They showed that porosity diminution leads to a complete decrease in melting time and system thermal storage capacity.

Through the brief review performed above, it appears that porosity and PPI are relevant parameters. Thereby, framed in such a general background, the primary intent of this study is to thoroughly assess the effects of PPI and porosity on heat transfer during melting/solidifying processes in a PCM-saturated open-ended rectangular porous channel under forced convection. Furthermore, as far as the numerical method is concerned, the enthalpy–SRT–LB method with triple distribution functions (TDF) is implemented herein.

This paper is laid out as follows: after the literature review (Section 1), the physical problem and the mathematical model, along with boundary conditions, are presented in Section 2. Afterward, a brief description of the numerical method is provided in Section 3. In Section 4, mesh control and code validation are performed. Numerical findings are presented and commented on in Section 5. Finally, major conclusions from this investigation are drawn in Section 6.

2. Model Description and Mathematical Formulation

In this section, firstly, the model geometry, the boundary conditions, and the simplifying hypotheses are outlined. Then, the governing equations with the adopted numerical resolution method are presented.

2.1. Considered Configuration

The physical model with a Cartesian coordinate system considered in this study is depicted in Figure 1. It is an open-ended rectangular channel of dimensions $\mathrm{L}\times \mathrm{H}$ completely filled with a porous medium (metallic foam) and a PCM (paraffin) embedded in its pores. Its lower and upper walls are adiabatic, impermeable, and non-slip. The uniform velocity and the hot temperature of the fluid at the channel entrance (west) were set to ${\mathrm{U}}_{\mathrm{i}\mathrm{n}}$ and ${\mathrm{T}}_{\mathrm{f},\mathrm{h}}$, respectively, to start the charging (melting) process. At the channel exit (east), the velocity and the temperature were set to $-{\mathrm{U}}_{\mathrm{i}\mathrm{n}}$ and ${\mathrm{T}}_{\mathrm{f},\mathrm{c}}$ ($<{\mathrm{T}}_{\mathrm{f},\mathrm{h}}$) to initiate the discharging process. The thermal input condition was wisely set such that the PCM could start to melt (${\mathrm{T}}_{\mathrm{f},\mathrm{h}}>{\mathrm{T}}_{\mathrm{m}}$) and the output conditions were assumed to be fully developed since the geometry is two-dimensional.

2.2. Computational Assumptions

Some assumptions were set out to simplify and establish the mathematical model to be dealt numerically. Forced convection prevailed in the channel compared to other heat transfer modes regardless of the porous medium interface or surface. Buoyancy and radiation effects were not considered. The flow was deemed to be laminar and unsteady. The working fluid was Newtonian, incompressible, and viscous, while it formed a second phase, i.e., the LTNE, with the porous matrix (first phase). It is worth indicating that the medium’s thermo-physical properties were constant, homogeneous, and isotropic within the deemed temperature and Re values. In addition, at the REV scale, homogenization of the governing equations reflects the porous medium effect. Hence, the Darcy–Brinkmann–Forchheimer (DBF) model was adopted here to model and investigate the convective fluid flow in the porous medium at Reynolds numbers of 200 and 400. Furthermore, the solid matrix was assumed to be in an LTNE state with the fluid. For more details, interested readers may refer to, for example, Nield and Bejan [25].

2.3. Mathematical Model

In light of the above-mentioned assumptions, the two-dimensional macroscopic conservation equations for mass, momentum, and energy at the REV scale are as follows [26,27,28]:

$$\nabla \cdot \overrightarrow{\mathrm{u}}=0,$$

(1)

$$\frac{\partial \overrightarrow{\mathrm{u}}}{\partial \mathrm{t}}+\left(\overrightarrow{\mathrm{u}}·\nabla \right)\left({\mathsf{\epsilon }}^{-1}\overrightarrow{\mathrm{u}}\right)=-\nabla \left(\mathsf{\epsilon }\mathrm{P}\right)+{\mathsf{\nu }}_{\mathrm{f}}{\nabla }^2\overrightarrow{\mathrm{u}}-\mathsf{\epsilon }\left(\frac{{\mathsf{\nu }}_{\mathrm{f}}}{\mathrm{K}}+\frac{{\mathrm{F}}_{\mathsf{\epsilon }}}{\sqrt{\mathrm{K}}}\Vert \overrightarrow{\mathrm{u}}\Vert \right)\overrightarrow{\mathrm{u}},$$

(2)

$$\mathsf{\epsilon }{\mathsf{\rho }}_{\mathrm{f}}{\mathrm{C}}_{\mathrm{p}\mathrm{f}}\left(\frac{\partial {\mathrm{T}}_{\mathrm{f}}}{\partial \mathrm{t}}+\overrightarrow{\mathrm{u}}·\nabla {\mathrm{T}}_{\mathrm{f}}\right)=\nabla ·\left({\mathsf{\lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f},\mathrm{f}}\nabla {\mathrm{T}}_{\mathrm{f}}\right)+{\mathrm{h}}_{\mathrm{s}\mathrm{f}}{\mathrm{a}}_{\mathrm{s}\mathrm{f}}\left({\mathrm{T}}_{\mathrm{s}}-{\mathrm{T}}_{\mathrm{f}}\right)-\mathsf{\epsilon }{\mathsf{\rho }}_{\mathrm{f}}{\mathrm{L}}_{\mathrm{a}}\frac{\partial \mathsf{\Gamma }}{\partial \mathrm{t}}+\mathsf{\Phi },$$

(3)

$$\left(1-\mathsf{\epsilon }\right){\mathsf{\rho }}_{\mathrm{s}}{\mathrm{C}}_{\mathrm{p}\mathrm{s}}\frac{\partial {\mathrm{T}}_{\mathrm{s}}}{\partial \mathrm{t}}=\nabla ·\left({\mathsf{\lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f},\mathrm{s}}\nabla {\mathrm{T}}_{\mathrm{s}}\right)+{\mathrm{h}}_{\mathrm{s}\mathrm{f}}{\mathrm{a}}_{\mathrm{s}\mathrm{f}}\left({\mathrm{T}}_{\mathrm{f}}-{\mathrm{T}}_{\mathrm{s}}\right)·$$

(4)

where $\overrightarrow{\mathrm{u}}$, $\mathrm{P}$, ${\mathrm{T}}_{\mathrm{f}}$, ${\mathrm{T}}_{\mathrm{s}}$, $\mathsf{\epsilon }$, $\mathsf{\rho }$, ${\mathsf{\nu }}_{\mathrm{f}}$, ${\mathrm{C}}_{\mathrm{p}}$, and ${\mathsf{\lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ are the velocity vector field, pressure, fluid and porous medium temperatures, copper porosity, density, kinematic viscosity of PCM, specific heat capacity, and effective thermal conductivity, respectively. Subscripts $\mathrm{f}$ and $\mathrm{s}$ denote the fluid and solid phases, respectively. Moreover, $\mathsf{\Gamma }$, ${\mathrm{L}}_{\mathrm{a}}$, ${\mathrm{a}}_{\mathrm{s}\mathrm{f}}$, and ${\mathrm{h}}_{\mathrm{s}\mathrm{f}}$ are the liquid fraction in pore space, the PCM latent heat, specific surface area of the porous matrix, and the local interfacial heat transfer coefficient between copper and paraffin, respectively. Note that the mutual establishment of heat transfer between the porous medium and the PCM is reflected via the two energy equations (Equations (3) and (4)). Simply stated, the melting heat transfer was modeled differently in the PCM and the metallic foam.

The liquid fraction $\mathsf{\Gamma }$ appearing in the third term on the right-hand side of Equation (3) could be determined according to the enthalpy method as follows [29]:

$$\mathsf{\Gamma }=\{\begin{array}{lll}0& & \mathrm{T}<{\mathrm{T}}_{\mathrm{m}}-\Delta \mathrm{T}\\ \left(\mathrm{T}-{\mathrm{T}}_{\mathrm{m}}+\Delta \mathrm{T}\right)/2\Delta \mathrm{T}& \mathrm{if}& {\mathrm{T}}_{\mathrm{m}}-\Delta \mathrm{T}\le \mathrm{T}\le {\mathrm{T}}_{\mathrm{m}}+\Delta \mathrm{T}\\ 1& & \mathrm{T}>{\mathrm{T}}_{\mathrm{m}}+\Delta \mathrm{T}\end{array}.$$

(5)

Since the metal and the PCM were in thermal equilibrium at the channel entrance (west), the matrix temperature was next raised to initiate the PCM melting. $\Delta \mathrm{T}$ is the PCM’s melting temperature range whose value was chosen to be small to ensure numerical stability. For the computations carried out here, it was set to 0.04.

The last term of Equation (3) points out the PCM viscous dissipation, which can be determined as follows [30]:

$$\mathsf{\Phi }=\mathsf{\epsilon }\left(\frac{\mathsf{\mu }}{\mathrm{K}}+\frac{{\mathrm{F}}_{\mathsf{\epsilon }}}{\sqrt{\mathrm{K}}}\Vert \overrightarrow{\mathrm{u}}\Vert \right){\Vert \overrightarrow{\mathrm{u}}\Vert }^2+\mathsf{\mu }\left(2\left[{\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{v}}{\partial \mathrm{y}}\right)}^2\right]+{\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}+\frac{\partial \mathrm{v}}{\partial \mathrm{x}}\right)}^2\right)·$$

(6)

2.4. Metal Foam Parameters

The empirical quantities ${\mathrm{a}}_{\mathrm{s}\mathrm{f}}$, ${\mathrm{h}}_{\mathrm{s}\mathrm{f}}$, and ${\mathsf{\lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ (Equations (1)–(4)) characterizing the thermal properties and the metallic foam geometry can be correlated as follows [31,32,33]:

$${\mathrm{a}}_{\mathrm{s}\mathrm{f}}=3\mathsf{\pi }{\mathrm{d}}_{\mathrm{f}}\mathrm{G}/0.59/{\mathrm{d}}_{\mathrm{p}}^2,\mathrm{where}\mathrm{G}=1-{\mathrm{e}}^{-(1-\mathsf{\epsilon })/0.004}\left(\mathrm{shape}\mathrm{function}\right),$$

(7)

with

$${\mathrm{d}}_{\mathrm{f}}=1.18{\left(\left(1-\mathsf{\epsilon }\right)/3\mathsf{\pi }\right)}^{1/2}{\mathrm{d}}_{\mathrm{p}}\mathrm{and}{\mathrm{d}}_{\mathrm{p}}=22.4\times {10}^{-3}/\mathsf{\omega },$$

(8)

$${\mathrm{h}}_{\mathrm{s}\mathrm{f}}=\{\begin{array}{lll}0.76.{\mathrm{Re}}_{\mathrm{d}}^{0.4}{\mathrm{Pr}}^{0.37}\frac{{\mathsf{\lambda }}_{\mathrm{f}}}{{\mathrm{d}}_{\mathrm{f}}}& & 1\le {\mathrm{Re}}_{\mathrm{d}}\le 40\\ 0.52.{\mathrm{Re}}_{\mathrm{d}}^{0.5}{\mathrm{Pr}}^{0.37}\frac{{\mathsf{\lambda }}_{\mathrm{f}}}{{\mathrm{d}}_{\mathrm{f}}}& \mathrm{for}& 140\le {\mathrm{Re}}_{\mathrm{d}}\le {10}^3\\ 0.26.{\mathrm{Re}}_{\mathrm{d}}^{0.6}{\mathrm{Pr}}^{0.37}\frac{{\mathsf{\lambda }}_{\mathrm{f}}}{{\mathrm{d}}_{\mathrm{f}}}& & {10}^3\le {\mathrm{Re}}_{\mathrm{d}}\le {2.10}^5\end{array},$$

(9)

where ${\mathrm{Re}}_{\mathrm{d}}\left(={\mathrm{d}}_{\mathrm{f}}{\mathrm{U}}_{\mathrm{i}\mathrm{n}}/\mathsf{\epsilon }{\mathsf{\nu }}_{\mathrm{f}}\right)$ is the pore Reynolds number, and ${\mathrm{d}}_{\mathrm{f}}$, ${\mathrm{d}}_{\mathrm{p}}$, and $\mathsf{\omega }$ are the ligament diameter, the pore size, and the pore density, respectively.

As the LTNE assumption was invoked in present the study, the effective thermal conductivity could be calculated as follows [9,34,35]:

$${\mathsf{\lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f}}=1/\sqrt{2}\left({\mathrm{R}}_{\mathrm{A}}+{\mathrm{R}}_{\mathrm{B}}+{\mathrm{R}}_{\mathrm{C}}+{\mathrm{R}}_{\mathrm{D}}\right),$$

(10)

where

$${\mathrm{R}}_{\mathrm{A}}=4\mathsf{\sigma }/\left(\left(2{\mathrm{e}}^2+\mathsf{\pi }\mathsf{\sigma }\left(1-\mathrm{e}\right){\mathsf{\lambda }}_{\mathrm{s}}+(4-2{\mathrm{e}}^2-\mathsf{\pi }\mathsf{\sigma }\left(1-\mathrm{e}\right)\right){\mathsf{\lambda }}_{\mathrm{f}}\right),$$

(11)

$${\mathrm{R}}_{\mathrm{B}}={\left(\mathrm{e}-2\mathsf{\sigma }\right)}^2/\left(\left(\mathrm{e}-2\mathsf{\sigma }\right){\mathrm{e}}^2{\mathsf{\lambda }}_{\mathrm{s}}+\left(2\mathrm{e}-4\mathsf{\sigma }-\left(\mathrm{e}-2\mathsf{\sigma }\right){\mathrm{e}}^2\right){\mathsf{\lambda }}_{\mathrm{f}}\right),$$

(12)

$${\mathrm{R}}_{\mathrm{C}}=\left(\sqrt{2}-2\mathrm{e}\right)/\left(\sqrt{2}\mathsf{\pi }{\mathsf{\sigma }}^2{\mathsf{\lambda }}_{\mathrm{s}}+\left(2-\sqrt{2}\mathsf{\pi }{\mathsf{\sigma }}^2\right){\mathsf{\lambda }}_{\mathrm{f}}\right),$$

(13)

$${\mathrm{R}}_{\mathrm{D}}=2\mathrm{e}/\left({\mathrm{e}}^2{\mathsf{\lambda }}_{\mathrm{s}}+\left(4-{\mathrm{e}}^2\right){\mathsf{\lambda }}_{\mathrm{f}}\right),$$

(14)

with

$$\mathrm{e}=0.16\mathrm{and}\mathsf{\sigma }={\left(\frac{\sqrt{2}(2-\frac{3\sqrt{2}}{4}{\mathrm{e}}^3-2\mathsf{\epsilon })}{\mathsf{\pi }(3-2\sqrt{2}\mathrm{e}-\mathrm{e})}\right)}^{1/2}·$$

(15)

Therefore,

$${\lambda }_{eff,f}={\lambda }_{eff}{}_{|{\lambda }_s=0},$$

(16)

$${\lambda }_{eff,s}={\lambda }_{eff}{}_{|{\lambda }_f=0}·$$

(17)

The permeability ($\mathrm{K}$) and the Forchheimer coefficient (${\mathrm{F}}_{\mathsf{\epsilon }}$) (see Equation (2)) were computed with the following correlations [31,32,33]:

$$\mathrm{K}=7.3\times {10}^{-4}{\mathrm{d}}_{\mathrm{p}}^2{(1-\mathsf{\epsilon })}^{-0.224}{({\mathrm{d}}_{\mathrm{f}}/{\mathrm{d}}_{\mathrm{p}})}^{-1.11},$$

(18)

$${\mathrm{F}}_{\mathsf{\epsilon }}=2.12\times {10}^{-3}{(1-\mathsf{\epsilon })}^{-0.132}{({\mathrm{d}}_{\mathrm{f}}/{\mathrm{d}}_{\mathrm{p}})}^{-1.63}·$$

(19)

2.5. Boundary and Initial Conditions

The necessary boundary conditions (BCs) and initial condition (IC) for the governing equations were as follows:

• $\mathrm{u}=0$; $\mathrm{v}=0$; ${\mathrm{T}}_{\mathrm{f}}={\mathrm{T}}_{\mathrm{c}}={\mathrm{T}}_0$ $\forall$ $\mathrm{x}$ and $\mathrm{y}$ (IC);
• $\mathrm{u}={\mathrm{u}}_{\mathrm{i}\mathrm{n}}$; $\mathrm{v}=0$; ${\mathrm{T}}_{\mathrm{f}}={\mathrm{T}}_{\mathrm{h}}$ at $\mathrm{x}=0$ and $0\le \mathrm{y}\le \mathrm{H}$ (west boundary);
• ${\nabla }_{\mathrm{x}}\mathrm{u}=0$; $\mathrm{v}=0$; ${\mathrm{T}}_{\mathrm{f}}={\mathrm{T}}_{\mathrm{c}}$ at $\mathrm{x}=\mathrm{L}$ and $0\le \mathrm{y}\le \mathrm{H}$ (east boundary);
• $\mathrm{u}=0$; $\mathrm{v}=0$ (no-slip condition); ${\nabla }_{\mathrm{y}}{\mathrm{T}}_{\mathrm{f}}={\nabla }_{\mathrm{y}}{\mathrm{T}}_{\mathrm{s}}=0$ (perfect-insulated) at $0\le \mathrm{x}\le \mathrm{L}$ and $\mathrm{y}=\mathrm{H}$ (upper boundary);
• $\mathrm{u}=0$; $\mathrm{v}=0$ (no-slip condition); ${\nabla }_{\mathrm{y}}{\mathrm{T}}_{\mathrm{f}}={\nabla }_{\mathrm{y}}{\mathrm{T}}_{\mathrm{s}}=0$ (perfect-insulated) at $0\le \mathrm{x}\le \mathrm{L}$ and $\mathrm{y}=0$ (lower boundary).

It should be pointed out that the fluid flow was inverted from right to left (see Figure 1) for the discharging cycle.

2.6. Dimensionless Mathematic Model and Key Parameters

Using the dimensionless variables below,

$$\begin{gathered} \left(\mathrm{X},\mathrm{Y}\right)=\left(\mathrm{x},\mathrm{y}\right)/\mathrm{H},\mathrm{U}=\mathrm{u}/{\mathrm{U}}_0,\mathrm{P}=\mathrm{p}/\mathsf{\rho }{\mathrm{U}}_0^2,\tilde{\mathrm{t}}=\mathrm{t}{\mathrm{U}}_0/\mathrm{H}, \\ \mathsf{\Theta }=\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{c}}\right)/\left({\mathrm{T}}_{\mathrm{h}}-{\mathrm{T}}_{\mathrm{c}}\right),\Delta {\mathrm{T}}_{\mathrm{r}\mathrm{e}\mathrm{f}}={\mathrm{T}}_{\mathrm{h}}-{\mathrm{T}}_{\mathrm{c}}, \end{gathered}$$

(20)

$$\begin{gathered} \mathrm{D}\mathrm{a}=\mathrm{K}/{\mathrm{H}}^2,\mathrm{Pr}={\mathsf{\nu }}_{\mathrm{f}}/{\mathsf{\alpha }}_{\mathrm{f}},\mathrm{Re}={\mathrm{U}}_0\mathrm{H}/{\mathsf{\nu }}_{\mathrm{f}},\mathrm{R}\mathrm{c}={(\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}})}_{\mathrm{s}}/{(\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}})}_{\mathrm{f}},\mathrm{K}\mathrm{r}={\mathsf{\lambda }}_{\mathrm{s}}/{\mathsf{\lambda }}_{\mathrm{f}}, \\ \mathrm{B}\mathrm{i}={\mathrm{h}}_{\mathrm{s}\mathrm{f}}{\mathrm{a}}_{\mathrm{s}\mathrm{f}}{\mathrm{H}}^2/{\mathsf{\lambda }}_{\mathrm{s}},Ste=C_{pf}(T_h-T_m)/L_a,\mathrm{E}\mathrm{c}={\mathrm{U}}_0^2/\left({\mathrm{C}}_{\mathrm{f}}\Delta {\mathrm{T}}_{\mathrm{r}\mathrm{e}\mathrm{f}}\right), \end{gathered}$$

(21)

the governing Equations (1)–(4) were converted to the following dimensionless forms [25,27,36]:

$$\nabla \cdot \overrightarrow{\mathrm{U}}=0,$$

(22)

$$\frac{\partial \overrightarrow{\mathrm{U}}}{\partial \tilde{\mathrm{t}}}+\left(\overrightarrow{\mathrm{U}}·\nabla \right)\left({\mathsf{\epsilon }}^{-1}\overrightarrow{\mathrm{U}}\right)=-\nabla \left(\mathsf{\epsilon }\mathrm{P}\right)+\frac{1}{\mathrm{Re}}{\nabla }^2\overrightarrow{\mathrm{U}}-\mathsf{\epsilon }\left(\frac{1}{\mathrm{Re}\mathrm{D}\mathrm{a}}+\frac{{\mathrm{F}}_{\mathsf{\epsilon }}}{\sqrt{\mathrm{D}\mathrm{a}}}\Vert \overrightarrow{\mathrm{U}}\Vert \right)\overrightarrow{\mathrm{U}},$$

(23)

$$\frac{\partial {\mathsf{\Theta }}_{\mathrm{f}}}{\partial \tilde{\mathrm{t}}}+\overrightarrow{\mathrm{U}}·\nabla {\mathsf{\Theta }}_{\mathrm{f}}=\frac{1}{\mathrm{Re}·\mathrm{Pr}}\nabla ·\left(\frac{{\mathsf{\lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f},\mathrm{f}}}{{\mathsf{\lambda }}_{\mathrm{f}}}\nabla \frac{{\mathsf{\Theta }}_{\mathrm{f}}}{\mathsf{\epsilon }}\right)+\mathrm{K}\mathrm{r}·\frac{\mathrm{B}\mathrm{i}}{\mathrm{Re}·\mathrm{Pr}}\left(\frac{{\mathsf{\Theta }}_{\mathrm{s}}-{\mathsf{\Theta }}_{\mathrm{f}}}{\mathsf{\epsilon }}\right)-\frac{1}{\mathrm{S}\mathrm{t}\mathrm{e}}\frac{\partial \mathsf{\Gamma }}{\partial \tilde{\mathrm{t}}}+\tilde{\mathsf{\Phi }},$$

(24)

$$\frac{\partial {\mathsf{\Theta }}_{\mathrm{s}}}{\partial \tilde{\mathrm{t}}}=\frac{\mathrm{K}\mathrm{r}}{\mathrm{R}\mathrm{c}}\frac{1}{\mathrm{Re}·\mathrm{Pr}}\nabla ·\left(\frac{{\mathsf{\lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f},\mathrm{s}}}{{\mathsf{\lambda }}_{\mathrm{s}}}\nabla \frac{{\mathsf{\Theta }}_{\mathrm{s}}}{1-\mathsf{\epsilon }}\right)-\frac{\mathrm{K}\mathrm{r}}{\mathrm{R}\mathrm{c}}·\frac{\mathrm{B}\mathrm{i}}{\mathrm{Re}·\mathrm{Pr}}\left(\frac{{\mathsf{\Theta }}_{\mathrm{s}}-{\mathsf{\Theta }}_{\mathrm{f}}}{1-\mathsf{\epsilon }}\right),$$

(25)

where

$$\tilde{\mathsf{\Phi }}=\mathsf{\epsilon }·\mathrm{E}\mathrm{c}\left\{\frac{1}{\mathrm{D}\mathrm{a}·\mathrm{Re}}+\frac{{\mathrm{F}}_{\mathsf{\epsilon }}}{\sqrt{\mathrm{D}\mathrm{a}}}\Vert \overrightarrow{\mathrm{U}}\Vert \right\}{\Vert \overrightarrow{\mathrm{U}}\Vert }^2+\frac{\mathrm{E}\mathrm{c}}{\mathrm{Re}}\left\{2\left[{\left(\frac{\partial \overrightarrow{\mathrm{U}}}{\partial \mathrm{X}}\right)}^2+{\left(\frac{\partial \overrightarrow{\mathrm{V}}}{\partial \mathrm{Y}}\right)}^2\right]+{\left(\frac{\partial \overrightarrow{\mathrm{U}}}{\partial \mathrm{Y}}+\frac{\partial \overrightarrow{\mathrm{V}}}{\partial \mathrm{X}}\right)}^2\right\}·$$

(26)

The corresponding dimensionless BCs and IC were as follows:

• $\mathrm{U}=0$; $\mathrm{V}=0$; ${\mathsf{\Theta }}_{\mathrm{f}}=0$ $\forall$ $\mathrm{X}$ and $\mathrm{Y}$ (IC);
• $\mathrm{U}=1$; $\mathrm{V}=0$; ${\mathsf{\Theta }}_{\mathrm{f},\mathrm{h}}=1$ at $\mathrm{X}=0$ and $0\le \mathrm{Y}\le 1$ (west boundary);
• ${\nabla }_{\mathrm{X}}\mathrm{U}=0$; $\mathrm{V}=0$; ${\mathsf{\Theta }}_{\mathrm{f},\mathrm{c}}=0$ at $\mathrm{X}=\mathrm{L}/\mathrm{H}$ and $0\le \mathrm{Y}\le 1$ (east boundary);
• $\mathrm{U}=0$; $\mathrm{V}=0$; ${\nabla }_{\mathrm{Y}}{\mathsf{\Theta }}_{\mathrm{f},\mathrm{h}}={\nabla }_{\mathrm{Y}}{\mathsf{\Theta }}_{\mathrm{f},\mathrm{s}}=0$ at $0\le \mathrm{X}\le \mathrm{L}/\mathrm{H}$ and $\mathrm{Y}=1$ (upper boundary);
• $\mathrm{U}=0$; $\mathrm{V}=0$; ${\nabla }_{\mathrm{Y}}{\mathsf{\Theta }}_{\mathrm{f},\mathrm{h}}={\nabla }_{\mathrm{Y}}{\mathsf{\Theta }}_{\mathrm{f},\mathrm{s}}=0$ at $0\le \mathrm{X}\le \mathrm{L}/\mathrm{H}$ and $\mathrm{Y}=0$ (lower boundary).

3. Lattice Boltzmann Methodology

For a long time now, the LBM aroused much interest in CFD-related fields. Unlike traditional CFD approaches (volumes, finite elements, etc.), such a method models a fluid using the kinetics of discrete particles which propagate (streaming step) and collide (relaxation step) on a discrete lattice. Generally, its mechanism consists of handling the streaming of fluid particle clusters via a monitor called the density distribution function (DDF) which controls the particle evolution position $\mathrm{x}$ at each time $\mathrm{t}$ with a certain discrete distribution velocity ${\overrightarrow{\mathrm{e}}}_{\mathrm{i}}$ to move regularly in the lattice toward bordering sights. Explicitly, at each time step, the particle circulation is achieved by continuous diffusion where the particles maintain their velocity direction, and collision occurs when particles change velocity. In addition, it has advantages in terms of geometry flexibility, a simple computation procedure, high parallelism, an ability to deal complex boundary conditions, and an intrinsic transient feature.

Recently, research corroborated the LBM’s success at the pore or REV scales for handling flows with and without phase change in porous media (e.g., References [4,26]). The first approach is often avoided owing to the geometry and the computational domain size, while the second approach REV is frequently adopted because of the averaging method. However, so far, and to the best of our knowledge, none of these approaches were used to deal with phase change, let alone forced convection under LTNE. In these circumstances, we studied the phase change process with forced convection under LTN in the REV-scale SRT–LBM framework [26].

3.1. Lattice Boltzmann Equation (LBE) for Dynamic Field

In this study, the LB model adopted to simulate the flow of the fluid was the SRT–LBM (also called the Bhatnagar–Gross–Krook (BGK) model) which leans on the evolution equation of the distribution function of the particle velocity density ${\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)$ [37].

$$\underset{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{a}\min \mathrm{g}}{\underbrace{{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}+{\mathrm{e}}_{\mathrm{i}}\mathsf{\delta }\mathrm{t},\mathrm{t}+{\mathrm{e}}_{\mathrm{i}}\mathsf{\delta }\mathrm{t}\right)-{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)}}=-\mathsf{\delta }\mathrm{t}{\mathsf{\omega }}_{\mathrm{i}}\underset{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}}{\underbrace{\left[{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)-{\mathrm{f}}_{\mathrm{i}}^{\mathrm{e}\mathrm{q}}\left(\mathrm{x},\mathrm{t}\right)\right]}}+\mathsf{\delta }\mathrm{t}·\underset{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}{\underbrace{{\mathrm{F}}_{{\mathrm{e}}_{\mathrm{i}}}}}·$$

(27)

The left-hand term represents advection (streaming), the first right-hand term is the discrete collision term, and the second right-hand term is the external force.

In Equation (27), ${\mathsf{\omega }}_{\mathrm{i}}(=1/{\mathsf{\tau }}_{\mathrm{v}};{\mathsf{\tau }}_{\mathrm{v}}=3\mathsf{\nu }+0.5)$, $\mathsf{\delta }\mathrm{t}$, and ${\mathrm{f}}_{\mathrm{i}}^{\mathrm{e}\mathrm{q}}(\mathrm{x},\mathrm{t})$ are the single relaxation collision frequency, the lattice time step, and the equilibrium fluid distribution function (EFDF), respectively. Under the well-known D2Q9 lattice model (Figure 2), the latter term can be expressed as follows (e.g., Reference [17]):

$${\mathrm{f}}_{\mathrm{i}}^{\mathrm{e}\mathrm{q}}=\mathsf{\rho }{\mathrm{w}}_{\mathrm{i}}\left(1+\frac{{\overrightarrow{\mathrm{e}}}_{\mathrm{i}}·\overrightarrow{\mathrm{u}}}{{\mathrm{c}}_{\mathrm{s}}^2}+\frac{\overrightarrow{\mathrm{u}}\otimes \overrightarrow{\mathrm{u}}:\left({\overrightarrow{\mathrm{e}}}_{\mathrm{i}}\otimes {\overrightarrow{\mathrm{e}}}_{\mathrm{i}}-{\mathrm{c}}_{\mathrm{s}}^2\mathrm{I}\right)}{2{\mathrm{c}}_{\mathrm{s}}^4\mathsf{\epsilon }}\right)·$$

(28)

Moreover, the model discrete velocity ${\overrightarrow{\mathrm{e}}}_{\mathrm{k}}$ can be given by

$${\overrightarrow{\mathrm{e}}}_{\mathrm{k}}=\{\begin{array}{ll}0\overrightarrow{\mathrm{i}}+0\overrightarrow{\mathrm{j}},& \mathrm{k}=0\\ \mathrm{c}\left(\cos \left((\mathrm{k}-1)\mathsf{\pi }/2\right)\overrightarrow{\mathrm{i}}+\sin \left((\mathrm{k}-1)\mathsf{\pi }/2\right)\overrightarrow{\mathrm{j}}\right),& \mathrm{k}=1,2,3,4,\\ \sqrt{2}\mathrm{c}\left[\cos \left((2\mathrm{k}-9)\mathsf{\pi }/4\right)\overrightarrow{\mathrm{i}}+\sin \left((2\mathrm{k}-9)\mathsf{\pi }/4\right)\overrightarrow{\mathrm{j}}\right],& \mathrm{k}=5,6,7,8\end{array}$$

(29)

where $\mathrm{c}(=\mathsf{\delta }\mathrm{x}/\mathsf{\delta }\mathrm{t}=1;\mathsf{\delta }\mathrm{x}=\mathsf{\delta }\mathrm{t})$, ${\mathrm{c}}_{\mathrm{s}}(=\mathrm{c}/\sqrt{3})$, and $\mathrm{I}$ are the streaming speed, the sound speed, and the unit tensor, respectively. The equilibrium weighting coefficient ${\mathrm{w}}_{\mathrm{i}}$ in Equation (28) was set as ${\mathrm{w}}_0=4/9$, ${\mathrm{w}}_{1,4}=1/9$, and ${\mathrm{w}}_{5,8}=1/36$.

The force term pointed out in Equation (27) can be expressed as follows [4]:

$${\mathrm{F}}_{{\mathrm{e}}_{\mathrm{i}}}={\mathrm{w}}_{\mathrm{i}}\mathsf{\rho }(1-\frac{1}{2{\mathsf{\tau }}_{\mathrm{v}}})\left[\frac{{\overrightarrow{\mathrm{e}}}_{\mathrm{i}}·\overrightarrow{\mathrm{F}}}{{\mathrm{c}}_{\mathrm{s}}^2}+\frac{\overrightarrow{\mathrm{u}}·\overrightarrow{\mathrm{F}}:\left({\overrightarrow{\mathrm{e}}}_{\mathrm{i}}{\overrightarrow{\mathrm{e}}}_{\mathrm{i}}-{\mathrm{c}}_{\mathrm{s}}^2\mathrm{I}\right)}{\mathsf{\epsilon }{\mathrm{c}}_{\mathrm{s}}^4}\right]\mathrm{with}\overrightarrow{\mathrm{F}}=-\left(\frac{{\mathsf{\nu }}_{\mathrm{f}}}{\mathrm{K}}+\frac{{\mathrm{F}}_{\mathsf{\epsilon }}}{\sqrt{\mathrm{K}}}\Vert \overrightarrow{\mathrm{U}}\Vert \right)\overrightarrow{\mathrm{U}}·$$

(30)

In the current model, the macroscopic density and velocity could be obtained by

$$\mathsf{\rho }={\displaystyle \sum _{\mathrm{i}}{\mathrm{f}}_{\mathrm{i}}}\mathrm{and}\mathsf{\rho }\overrightarrow{\mathrm{u}}={\displaystyle \sum _{\mathrm{i}}{\mathrm{f}}_{\mathrm{i}}}{\overrightarrow{\mathrm{e}}}_{\mathrm{i}}+\mathsf{\delta }\mathrm{t}\overrightarrow{\mathrm{F}}/2·$$

(31)

Note that the computation is typically split into a collision and streaming step as follows:

$${\mathrm{f}}_{\mathrm{i}}^{*}\left(\mathrm{x},\mathrm{t}\right)={\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)-\mathsf{\delta }\mathrm{t}{\mathsf{\omega }}_{\mathrm{i}}\left({\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)-{\mathrm{f}}_{\mathrm{i}}^{\mathrm{e}\mathrm{q}}\left(\mathrm{x},\mathrm{t}\right)\right)+\mathsf{\delta }\mathrm{t}{\mathrm{F}}_{{\mathrm{e}}_{\mathrm{i}}},$$

(32)

$${\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}+{\overrightarrow{\mathrm{e}}}_{\mathrm{i}}\mathsf{\delta }\mathrm{t},\mathrm{t}+{\mathrm{e}}_{\mathrm{i}}\mathsf{\delta }\mathrm{t}\right)={\mathrm{f}}_{\mathrm{i}}^{*}\left(\mathrm{x},\mathrm{t}\right)$$

(33)

where ${\mathrm{f}}_{\mathrm{i}}^{*}\left(\mathrm{x},\mathrm{t}\right)$ is the post-collision density distribution function.

It is noteworthy that the execution order of these two steps is arbitrary and may vary from one code to another for implementation reasons.

3.2. Lattice Boltzmann Equation (LBE) for Thermal Field

Here, we follow an approach similar to that which was just adopted for the dynamic field. Simply stated, we still use the LBGK with two thermal density distribution functions, ${\mathrm{g}}_{\mathrm{i},\mathrm{f};\mathrm{s}}$, to predict the thermal field for the PCM and porous matrix using the following scheme [4,28]:

$$\begin{array}{l}{\mathrm{g}}_{\mathrm{f},\mathrm{i}}\left(\mathrm{x}+{\overrightarrow{\mathrm{e}}}_{\mathrm{i}}\mathsf{\delta }\mathrm{t},\mathrm{t}+\mathsf{\delta }\mathrm{t}\right)-{\mathrm{g}}_{\mathrm{f},\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)=-{\mathsf{\omega }}_{\mathrm{T},\mathrm{f}}\left({\mathrm{g}}_{\mathrm{f},\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)-{\mathrm{g}}_{\mathrm{f},\mathrm{i}}^{\mathrm{e}\mathrm{q}}\left(\mathrm{x},\mathrm{t}\right)\right)+\left(1+\mathsf{\delta }\mathrm{t}{\partial }_{\mathrm{t}}/2\right)\mathsf{\delta }\mathrm{t}\mathrm{S}{\mathrm{r}}_{\mathrm{i},\mathrm{f}}\\ +\mathsf{\delta }\mathrm{t}{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right){\mathrm{q}}_{\mathrm{i}}\end{array}$$

(34)

$${\mathrm{g}}_{\mathrm{s},\mathrm{i}}\left(\mathrm{x}+{\overrightarrow{\mathrm{e}}}_{\mathrm{i}}\mathsf{\delta }\mathrm{t},\mathrm{t}+\mathsf{\delta }\mathrm{t}\right)-{\mathrm{g}}_{\mathrm{s},\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)=-{\mathsf{\omega }}_{\mathrm{T},\mathrm{s}}\left({\mathrm{g}}_{\mathrm{s},\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)-{\mathrm{g}}_{\mathrm{s},\mathrm{i}}^{\mathrm{e}\mathrm{q}}\left(\mathrm{x},\mathrm{t}\right)\right)+\left(1+\mathsf{\delta }\mathrm{t}{\partial }_{\mathrm{t}}/2\right)\mathsf{\delta }\mathrm{t}\mathrm{S}{\mathrm{r}}_{\mathrm{i},\mathrm{s}},$$

(35)

where ${\mathsf{\omega }}_{\mathrm{T},\mathrm{f};\mathrm{s}}(=1/{\mathsf{\tau }}_{\mathrm{T},\mathrm{f};\mathrm{s}})$ is the single relaxation collision frequency for the solid and fluid temperature distribution functions which can be estimated through the dimensionless relaxation times ${\mathsf{\tau }}_{\mathrm{T},\mathrm{f};\mathrm{s}}$ as follows [28,38]:

$${\mathsf{\tau }}_{\mathrm{T},\mathrm{f}}=3{\mathsf{\alpha }}_{\mathrm{e},\mathrm{f}}/(\mathsf{\delta }\mathrm{t}{\mathrm{c}}^2)+0.5\mathrm{with}{\mathsf{\alpha }}_{\mathrm{e},\mathrm{f}}={\mathrm{k}}_{\mathrm{e},\mathrm{f}}/\left(\mathsf{\epsilon }{(\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}})}_{\mathrm{f}}\right),$$

(36)

$${\mathsf{\tau }}_{\mathrm{T},\mathrm{s}}=3{\mathsf{\alpha }}_{\mathrm{e},\mathrm{s}}/(\mathsf{\delta }\mathrm{t}{\mathrm{c}}^2)+0.5\mathrm{with}{\mathsf{\alpha }}_{\mathrm{e},\mathrm{s}}={\mathrm{k}}_{\mathrm{e},\mathrm{s}}/\left((1-\mathsf{\epsilon }){(\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}})}_{\mathrm{s}}\right),$$

(37)

where ${\mathsf{\alpha }}_{\mathrm{e},\mathrm{f},\mathrm{s}}$ is the effective diffusivity.

The equilibrium distribution function ${\mathrm{g}}_{\mathrm{i},\mathrm{f};\mathrm{s}}^{\mathrm{e}\mathrm{q}}$ (Equations (34) and (35)) can be expressed as

$${\mathrm{g}}_{\mathrm{f},\mathrm{i}}^{\mathrm{e}\mathrm{q}}={\mathrm{w}}_{\mathrm{i}}{\mathrm{T}}_{\mathrm{f}}\left(1+{\overrightarrow{\mathrm{e}}}_{\mathrm{i}}·\overrightarrow{\mathrm{u}}/(\mathsf{\epsilon }{\mathrm{c}}_{\mathrm{s}}^2)\right)\mathrm{and}{\mathrm{g}}_{\mathrm{s},\mathrm{i}}^{\mathrm{e}\mathrm{q}}={\mathrm{w}}_{\mathrm{i}}{\mathrm{T}}_{\mathrm{s}}·$$

(38)

The penultimate terms (source terms) $\mathrm{S}{\mathrm{r}}_{\mathrm{i},\mathrm{f};\mathrm{s}}$ can be computed as follows [28,38]:

$$\mathrm{S}{\mathrm{r}}_{\mathrm{i},\mathrm{f}}={\mathrm{w}}_{\mathrm{i}}\left(\frac{{\mathrm{L}}_{\mathrm{a}}}{{\mathrm{C}}_{\mathrm{p},\mathrm{f}}}\left[\frac{\mathsf{\gamma }(\mathrm{t}+\mathsf{\delta }\mathrm{t})-\mathsf{\gamma }(\mathrm{t})}{\mathsf{\delta }\mathrm{t}}\right]+\frac{\mathrm{h}({\mathrm{T}}_{\mathrm{s}}-{\mathrm{T}}_{\mathrm{f}})}{\mathsf{\epsilon }{(\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}})}_{\mathrm{f}}}\right),$$

(39)

$$\mathrm{S}{\mathrm{r}}_{\mathrm{i},\mathrm{s}}={\mathrm{w}}_{\mathrm{i}}\left(\frac{\mathrm{h}({\mathrm{T}}_{\mathrm{s}}-{\mathrm{T}}_{\mathrm{f}})}{(1-\mathsf{\epsilon }){(\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}})}_{\mathrm{s}}}\right)·$$

(40)

As for the last term (Equation (34)), it can be computed as follows [39]:

$${\mathrm{q}}_{\mathrm{i}}=-({\mathrm{f}}_{\mathrm{i}}-{\mathrm{f}}_{\mathrm{i}}^{\mathrm{e}\mathrm{q}})({\overrightarrow{\mathrm{e}}}_{\mathrm{i}}-\overrightarrow{\mathrm{u}})({\overrightarrow{\mathrm{e}}}_{\mathrm{i}}-\overrightarrow{\mathrm{u}}):\Delta \overrightarrow{\mathrm{u}}·$$

(41)

The macroscopic temperature can be computed from

$${\mathrm{T}}_{\mathrm{f}}={\displaystyle \sum {\mathrm{g}}_{\mathrm{f}\mathrm{i}}}\mathrm{and}{\mathrm{T}}_{\mathrm{s}}={\displaystyle \sum {\mathrm{g}}_{\mathrm{s}\mathrm{i}}}·$$

(42)

4. Mesh Control and Code Validation

For the sake of brevity, validation and grid dependency analyses were limited to two cases that were deemed relevant.

4.1. Mesh Control

To achieve the grid-mesh test, simulations were performed for Reynolds numbers of 200 and 400 using four grids, i.e., $100\times 50$, $200\times 100$, $300\times 150$, and $350\times 175$. The dimensionless $\mathrm{U}$-velocity is depicted in Figure 3 using these grids. As can be seen, the predicted values were very close to each other. As illustrated, the maximum difference between the first two meshes was approximately 0.8%. Then, it was reduced to reach a deviation of 0.2% between the last three grids. This suggests that the solution became independent from the $200\times 100$ grid. Therefore, such a grid was selected for the remaining computations.

4.2. Code Validation

It should be stressed that our in-house code was already validated [23]. However, before presenting and commenting on the results in the dedicated section, we ensured to further validate it using other previously published benchmarks.

For this ultimate validation, we compared the results obtained using the present thermal SRT–LBM with those of References [40,41]. Figure 4 shows the comparison of the dimensionless velocity and temperature results generated by forced convection in a rectangular porous channel differentially heated via its two lower and upper walls along the line X = 0.5 for three selected Da numbers of 0.01, 0.05, and 0.1, under the LTE conjecture. As shown in Figure 4a, the permeability increased with the Da number, thereby reducing the velocity profile flatness. As for the thermal field, the viscous effects decreased when the Da number increased (Figure 4b). As noted, our predictions corroborated these selected benchmarks. Furthermore, we supported our findings with maximum errors (%), which are gathered in

Table 1. Model validation (error %).

$\mathrm{Da}$	Error (%)/ $\mathrm{U}$-Velocity	Error (%)/ $\mathsf{\Theta }$
0.01	2.01	0.81
0.05	1.18	0.47
0.1	1.15	0.42

.

Figure 5 depicts the comparison of fluid and solid temperature evolution (${\mathsf{\Theta }}_f,{\mathsf{\Theta }}_s$) in a square domain including a porous metallic foam filled with PCM under natural convection with results of Reference [41] at Y = 0.5 and $\mathrm{D}\mathrm{a}={10}^{-2},\mathrm{Pr}=50,\mathrm{R}\mathrm{a}={10}^6,\mathrm{S}\mathrm{t}\mathrm{e}=1,\mathrm{N}{\mathrm{u}}_{\mathrm{i}}=0$. As can be seen, excellent matching was observed.

To sum up, such comparisons highlight good agreement with previous studies, further confirming the reliability of the built code.

5. Results and Comments

It is worth pointing out that that it is the numerical simulation of unsteady forced convection flow in an open channel filled with porous medium/PCM (see Section 2 and Section 3) via a thermal SRT–LBM which was targeted. In this study, it was the effect of pore density $\mathsf{\omega }$ ($\mathrm{P}\mathrm{P}\mathrm{I}=10,30$, and $60$) for three porosity values ($\mathsf{\epsilon }=0.7,0.8$, and $0.9$) and two Re numbers (200 and 400) which was mainly investigated, while the other parameters such as the numbers of Prandtl, Stephan, Eckert, and heat capacity ratio were held fixed to $\mathrm{Pr}=50,\mathrm{S}\mathrm{t}\mathrm{e}=1,\mathrm{E}\mathrm{c}=0,\mathrm{R}\mathrm{c}=1$.

5.1. PPI Effect on Velocity Field

Figure 6a–f illustrate the PPI effect on the streamwise velocity ($\mathrm{U}$) during charging and discharging processes. Already and as expected, PPI had an impact on such a velocity whatever the porosity. As shown, during the charging case, a low porosity (0.7) with a high PPI value (60) had a strong effect on such a velocity, which was slowed and undeveloped. It only started to exhibit an almost developed profile in the channel’s last third for the two Reynolds numbers (200 and 400) and at low PPI (=10), demonstrating that the force convection was weak. Likewise, it appeared that forced convection acted in the first third of the channel for the same porosity ($\mathsf{\epsilon }=0.7$) at PPI = 10 and 30 when Re increased. At PPI = 30 and at low porosity (0.7), forced convection mainly affected the channel middle (red zone). As for PPI = 60 and $\mathsf{\epsilon }=0.7$, the $\mathrm{U}$-velocity presented an increasingly twisted profile.

During the discharging process (Figure 6a–c), for the two Re values (Re = 200 and 400) and three porosity values, the increase in PPI slowed down the flow, revealing the convection’s role as being lesser.

Figure 6d–f presents the dynamic U-contours. They show that PPI had a manifest influence on such a velocity during the charging process. It stands to reason that this comment remains in the event of the discharging process even if it is not depicted here for the sake of brevity.

5.2. PPI Effect on the PCM and Metal Foam Temperature

To better determine the effect induced by the PPI and the porosity, the PCM temperature (${\mathsf{\Theta }}_f$) during charging and discharging is plotted in Figure 7a–c. As can be seen, the metal foam temperature (${\mathsf{\Theta }}_s$) was generally higher than that of the PCM temperature (${\mathsf{\Theta }}_f$) in the first two-thirds of the channel regardless of the porosity, the pore density, and Re. In this same region, the PCM temperature decreased with the increase in PPI. This demonstrates that, for a larger PPI, forced convection is limited due to the large specific surface of the metallic structure and its high conductivity compared to the PCM (paraffin). It can be observed that, as time elapsed, the forced convection strengthened, thereby raising the liquid PCM temperature, which went beyond the metal foam temperature in the channel’s last third. In addition, as shown in Figure 7a–c, as the porosity decreased and the PPI increased, the difference $\left|{\mathsf{\Theta }}_f-{\mathsf{\Theta }}_s\right|$ became larger, proving once more that LTNE is indeed involved. On the other hand, during the discharge process, ${\mathsf{\Theta }}_s$ was always higher than ${\mathsf{\Theta }}_f$. However, the difference $\left|{\mathsf{\Theta }}_f-{\mathsf{\Theta }}_s\right|$ decreased as the PPI increased, restricting forced convection and showing that LTNE decreased and approximated (tended toward) LTE. In this case, only one energy equation (instead of two) would suffice.

5.3. PPI Effect on Melting Front and PCM Temperature Field

The melting evolution with the PCM temperature distribution is exhibited in Figure 8a–i at a Reynolds number (Re) of $200$. Note that the melt front mushy zone is highlighted with colors other than red and blue, denoting the liquid and solid PCM regions, respectively. It was found that, as time elapsed, the melting front advanced more quickly when PPI = 30 than with other values, while the melting approached its end after 600 s. However, a higher PPI (60) slowed the heat transfer progress and, therefore, the melting front evolution. For $\mathsf{\epsilon }=0.9$ (Figure 8g–i), the increase in PPI improved the heat transfer between the PCM and metal foam and sped up the melting time due to the metal foam permeability influence.

The same trend can be seen in Figure 8d–f at $\mathsf{\epsilon }=0.8$. Simply put, the melting rate sped up for fairly low PPI values (30), while a larger PPI (60) mitigated and slowed the heat transfer rate and, therefore, the melting front evolution. However, the melting evolution increased with PPI and became stronger for higher porosity (0.9) as presented in Figure 8g–i owing to the increase in metal foam heat transfer. More specifically, a higher PPI and porosity (60 and 0.9) improved the heat transfer between the PCM and the metal foam for Re = 200 while reducing the melting time.

The distribution of streamlines, flow, and temperature fields during the charging process at Re = 400 are depicted in Figure 9a–i, for the same values of PPI and porosity as shown in the previous figures, i.e., $10\le \mathrm{P}\mathrm{P}\mathrm{I}\le 30$ and $0.7\le \mathsf{\epsilon }\le 0.9$.

For $\mathsf{\epsilon }=0.7$, the melting front advanced faster as time elapsed for small PPI (10 or 30) than larger PPI (60). As mentioned above, a larger PPI (60) associated with a small porosity (0.7) continued to slow the melting front evolution and limit forced convection even upon doubling Re. Meanwhile, considering $\mathsf{\epsilon }=0.8$ and upon raising the PPI, it can be observed that the PCM melting front sped up even more quickly over time suggesting that the metallic foam has the capacity to further improve the transfer heat. As for $\mathsf{\epsilon }=0.9$ with a PPI of 60 (the largest value deemed here) (Figure 9g–i), the melting front advanced rapidly and melted the whole PCM in a short time just above 600 s (not illustrated here).

6. Conclusions

This paper conducted a numerical study of the effect of pore density and porosity of a metallic foam on the phase change in an open channel under forced convection. Simulations were performed under LTNE conditions using the REV enthalpy–SRT–LB method. Effects of the studied parameters on the flow, thermal field, and melt front evolution were surveyed during the charging (melting) and discharging (solidifying) cycles, and demonstrated that the composite material insert has an undeniable effect on the phase change via its key parameters (PPI and porosity). Based on the investigation results, the following main conclusions are proposed:

• During the charging process, increasing PPI decelerates the velocity distribution whatever the porosity, while lower values of PPI speed up the velocity field whatever the porosity.
• During the charging process, it is the thermal conduction which seems to be involved in the first two-thirds of the channel, then it is the forced convection which prevails within the channel’s last third.
• During the discharging process, the increase in PPI improves thermal conduction and restricts forced convection.
• At a moderate Reynolds number (200), an increase in PPI (60) slows down the melting front evolution for $\mathsf{\epsilon }\ge 0.8$ while reducing the melting time at high porosity ($\mathsf{\epsilon }=0.9$).
• On the other hand, at Re = 400, an increase in PPI (60) slows down the melting front evolution at low porosity ($\mathsf{\epsilon }=0.7$) while dropping the melting time for $\mathsf{\epsilon }\ge 0.8$.