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

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.


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 solidliquid 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 enthalpybased 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 solidliquid 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.

Macroscopic Governing Equations
For heat transfer in metal-foam-based solid-liquid PCMs, the macroscopic governing equations are given by [9,12,23]: 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 F is determined by: According to the Boussinesq approximation, the body force G is given by G = −gβ(T − T 0 ). For flow through packed bed consisting of spherical particles, F φ and K are given by: d p in Equation (6) denotes the mean particle diameter.
For metal foams such as aluminum foam, F φ and K can be well predicted by the following empirical correlations [33]: where 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]: The effective thermal conductivities k e f 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 m f can be predicted by the following relation [35]: h m f depends on the foam structure, Reynolds number, and Prandtl number. For forced convective flow in metal foams, several empirical correlations for h m f have been proposed [8,34]. For natural convection flow through metal foams, the following empirical correlation is widely used [35]: where

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 {e i |i = 0, . . . , 4 } are given by: We first introduce the thermal CLB equation for solving Equation (3) without phasechange term based on the simplified CLB method [30]. The raw moments k T mn and central moments k T mn of the discrete temperature distribution function g f i of the fluid phase are defined as follows [30]: The simplified raw-moment n f i and central-moment n f i are given by [30] Through the transformation matrix M T and shift matrix N T , n f i and n f i can be determined by [30] The thermal CLB equation for solving Equation (12) can be divided into two parts: collision and streaming processes. Through shifting procedure ( n f i = N T n f i ), collision process is implemented in central-moment space as where n eq f i are equilibrium central moments, n * f i denotes the post-collision central moments, Θ= diag(ζ T , ζ α , ζ α , ζ e , ζ e ) is the relaxation matrix, S f i is the source term.
The streaming process is still implemented in velocity space as where . M T and N T are given by (c = 1): The equilibrium temperature distribution function g eq f i can be determined by g T n eq f i . Explicitly, the source term S f i is given by: where The temperature T f is defined as: The effective thermal diffusivity α e f is:

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 ∂ t (φρ l L a f l ) with the transient term ∂ t φρ f c p f T f in Equation (3), we can obtain: where H f = σc pl T f + L a f l is the enthalpy, and σ = is the heat capacity ratio. In the liquid phase ( f l = 1), H f = c pl T f + L a and σ l = 1; in the solid phase . 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 n eq f i and raw moments n eq f i should be chosen as follows: In the system, n f 0 is the only conserved quantity ( n f 0 = n f 0 = H f ). As the enthalpy H f is the basic evolution variable (g f i is now defined as the enthalpy distribution function), the shifting procedure n f i = N T n f i in the original thermal CLB equation should be modified. Explicitly, the modified shifting procedure is given by: 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 Explicitly, g eq f i in the velocity space is: Simultaneously, the temperature T f can be calculated by: where T f s and T f l (T f s ≤ T f l ) are solidus and liquidus temperatures, respectively; H f s and H f l are enthalpies corresponding to T f s and T f l , respectively. The liquid fraction f l is:

Thermal MRT-LB Equation
Equation (4) can be rewritten as: The MRT-LB equation for Equation (34) is: The collision process of the evolution Equation (35) is implemented in raw-moment space as: where n mi = M T g mi , n eq mi = M T g eq mi , and S mi = M T S mi . The streaming process is executed in velocity space as where g * mi = M −1 T n * mi . The equilibrium moment n eq mi is: The equilibrium temperature distribution function g eq mi (x, t) in the velocity space is given by: where 2 ∈ (0, 1) is a tunable parameter of the D2Q5 lattice model. Explicitly, the source term S mi is given by: The temperature T m is defined by: The effective thermal diffusivity

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 ≤ T 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: where is the characteristic temperature and α m = k m / ρ m c pm .
Entropy 2022, 24, x FOR PEER REVIEW 9 of 18 [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 CPU@2.80GHz; RAM: 8.0 GB). Ra 10 = , 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 (  In simulations, some required parameters are: δ x = δ y = δ t = 1 (c = 1), c pl = c ps = 1, σ = 1, Pr = 50, H v = 5.9, d p /L = 0.0135, λ = Γ = 10 3 , k f = 0.0005, J = 1, and 1 = 2 = 1/2. The effective thermal conductivities k e f and k em are simply determined by k e f = φk f (k f = k f l = k f s ) and k em = (1 − φ)k m , respectively. To eliminate the unphysical numerical diffusion, the relaxation parameters ζ 3 and ζ 4 (ζ 3 = ζ 4 = ζ e ) related to the second-order moments are determined by ζ e = 2 − ζ α [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 CPU@2.80GHz; RAM: 8.0 GB).
In Figure 2, the melting front at different values of Fo for Ra = 10 6 and 10 8 with Da = 10 −4 , φ = 0.8 and St = 1 are presented. For comparison purpose, the other parameters are chosen as F φ = 0.068, θ h = 1, θ i = θ c = θ 0 = 0, and θ melt = 0.3 (θ = (T − T c )/∆T) [9]. For Ra = 10 6 , N x × N y = 100 × 100 is employed, and for Ra = 10 8 , N x × N y = 300 × 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 , 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 ( 8 10 Ra = ), the effect of natural convection becomes stronger, then the melting front moves faster near the top wall (see Figure 2b). ), 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 0.001 Fo = , 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]. 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 , φ = 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 figure, at the very beginning ( 0.00005 Fo = ), 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 0.001 Fo = , 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].  . As Fo increases, na ural convection effect becomes stronger, resulting in more hot liquid to move upward and consequently, the PCM melts faster in the upper region of the cavity. Although th Rayleigh number is high, the natural convection effect is weak since it has been severel suppressed by the effect caused by the metal foam's ligament, and, during the meltin process, the heat conduction through the high thermal conductive metal foam dominate the heat transfer.  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. ural 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 0.002  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 , φ = 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 N T is an identity matrix and ζ i = 1/τ f l , the enthalpy-based CLB equation degrades into the enthalpy-based BGK-LB equation with α e f = c 2 s f τ f l − 0.5 δ 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 ζ e = 2 − ζ α , the effect of the numerical diffusion on numerical simulations can be effectively eliminated by the present CLB model.
). 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 ζ 2 ζ  F  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  and K given by Equations (6) and (7) are evaluated. In numerical simulations, the related parameters are set as  (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  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 F φ 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 φ and K given by Equations (6) and (7) are evaluated. In numerical simulations, the related parameters are set as θ h = 1, θ i = θ c = θ melt = θ 0 = 0. Numerical simulations are conducted based on N x × N y = 100 × 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 (φ = 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 φ 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 6 10 Ra = and 8 10 for different porosities with F φ and K given by Equation (7) are shown. As can be seen in the figure, when φ 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 6 10 Ra = and 8 10 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. In Figure 7, the total liquid fractions of Ra = 10 6 and 10 8 for different porosities with F φ and K given by Equation (7) are shown. As can be seen in the figure, when φ 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.

Conclusions
An enthalpy-based CLB model is constructed for investigating heat transfer in solidliquid 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 φ 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.