Pore-scale investigation on natural convection melting in a square cavity with gradient porous media

In this paper, a numerical study on the melting behavior of phase change material (PCM) with gradient porous media has been carried out at the pore scales. In order to solve the governing equations, a pore-scale lattice Boltzmann method with the double distribution functions is used, in which a volumetric LB scheme is employed to handle the boundary. The Monte Carlo random sampling is adopted to generate a microstructure of two-dimensional gradient foam metal which are then used to simulate the solid-liquid phase transition in the cavity. The effect of several factors, such as gradient porosity structure, gradient direction, Rayleigh number and particle diameters on the liquid fraction of PCM are systematically investigated. It is observed that the presence of gradient media affect significantly the melting rate and shortens full melting time compared to that for constant porosity by enhancing natural convection. The melting time of positive and negative gradients will change with Rayleigh number, and there is a critical value for Rayleigh number. Specifically, when Rayleigh number is below the critical value, the positive gradient is more advantageous, and when Rayleigh number exceeds the critical value, the negative gradient is more advantageous. Moreover, smaller particle diameters would lead to lower permeability and larger internal surfaces for heat transfer.

In order to improve the performance of LHTES, scholars have proposed many heat transfer enhancement techniques, including adding highly conductive particles [6,7], adding highly conductive metallic fins [8,9,10], using multiple PCM methods [11,12] , embedding PCMs in highly conductive porous media [13,14,15,16,17,18,19] and so on. Among these heat transfer enhancement methods, highly conductive porous media with high thermal conductivity, high heat penetration, high porosity and high specific surface area has been showing progressive potentials.
Metal foam is a typical porous media with high thermal conductivity. Metal foam perfectly maintains the excellent properties of base metals, such as high stability, light weight, high thermal conductivity, ductility, and improved heat transfer characteristics. Therefore, metal foam is a widely used thermal conductivity enhancers. Some studies have explored the problems of solid-liquid phase change heat transfer in PCM-filled metal foam [13,14,15,16,17,18,19].
Chen et al. [13] used an infrared camera and microscope to experimentally study the heat transfer phenomenon of PCM in foamed metal and used the double-distributed lattice Boltzmann method to carry out a numerical simulation.
It is found that the numerical value is consistent with the experimental results and the foam metal has a significant enhancement effect on the melting of PCM. Yang et al. [14] experimentally studied the effects of metal foam to melting and found that completely melting the PCM in metal foam takes over 1/3 less time than that of pure paraffin under the same conditions. Zhao et al. [15] numerically studied the melting behavior of paraffin in foam metal, it is noted that the Rayleigh number, porosity and pore density have a significant influence on on the melting and solidification process. Tao et al. [16] used the lattice Boltzmann method to study the heat storage performance of paraffin foam metal composites. The author studied the influence of porosity and pore density on the melting rate, and proposed the best metal foam structure with a porosity of 0.94 and a PPI of 45. Zhu et al. [17] used the finite volume method to explore the influence of the three strengthening methods of foam metal porosity, cold wall shape and discrete heat source on the thermal properties of foam metal/phase change materials. Yao et al. [18] conducted a visualized experiment to study the melting of paraffin in high porosity copper foam at pore scale. They concluded that the copper foam with a high porosity of 0.974 effectively extends the phase change interface and improves the heat storage of paraffin, while the reduction in the amount of latent heat is only 2.6%. Yang et al. [19] numerically and experimentally explored the influence of the inclination angle of the inclined cavity containing foamed metal and the cavity aspect ratio on the melting of PCM. It is found that the tilt angle has little effect when the aspect ratio is given, and the smaller aspect ratio is better than the larger aspect ratio when the aspect ratio is not fixed.
All the studies mentioned above only consider metal foams with fixed pore parameters. Recently, many investigations have found that gradient metal foams can further enhance melting heat transfer under the same conditions. Yang et al. [20] numerically investigated the melting process of sodium nitrate inside porous copper foam with linearly changed porosity. The numerical results show that porosity linearly increased from bottom to top could improve the heat transfer performance and shorten the completely melted time compared to that for constant porosity by enhancing natural convection. Yang et al. [21] numerically studied the solidification behavior of saturated distilled water in open-cell foam metal and compared it with ungraded foam metal. The authors found that the positive gradient porosity and negative gradient pore density structures have a faster solidification rate compared to the non-graded foams.
Zhu et al. [22] proposed an improved metal foam structure, which is composed of metal foam and finned metal foam with gradient pores. The finite volume method is also used to analyze the influence of structural parameters on energy storage performance. It is found that this structure can shorten the melting time by changing the melting sequence of the phase change material. Zhang et al. [23] numerically studied the melting behavior in gradient foam metal which was consist of three different homogeneous porosity slices, and results show that the gradient porosity structure can overcome the corner phenomenon at the bottom to increase the heat storage rate. Yang et al. [24] experimentally and numerically studied the effect of gradient porosity and gradient density in tube latent heat thermal energy storage. Results indicated that the positive gradient design of porosity can significantly reduce the melting time of PCMs filled in the pore space and obtain simultaneously a better temperature uniformity. Ghahremannezhad et al. [25] used a finite volume approach to numerically simulate the melting behavior of PCMs in gradient foam metal under different heating modes. They found that the direction of gradient porosity and PPI can affect the heat transfer rate. Hu et al. [26] use a three-dimensional model to numerically simulate the melting behavior in a gradient metal foam saturated with PCM and quantitatively explored the effect of gradient size and gradient difference on the heat storage characteristics of PCMs. It found that gradient metal foam effectively improves and accelerates the heat storage efficiency and there is an optimal gradient difference under a fixed average porosity. Marri et al. [27] experimentally numerically studied the thermal function of cylindrical foam metal/PCM composite heat sinks with gradients of porosity and density. The results show that the three-level gradient and the two-level gradient have comparable thermal performance, which are 4.4 times and 4 times stronger than the thermal performance of the uniform structure, respectively.
From the above experiments and numerical studies, it can be seen that the gradient pore structure has a significant impact on the heat transfer characteristics of PCM, and it is an effective means to enhance phase change heat transfer.
However, most of the research on the influence mechanism of the gradient metal foam on the phase change process of the phase change material (PCM) is based on the representative elementary volume (REV) models, and seldom is studied by the pore-scale models. Compared with the REV model, the pore-scale model can provide detailed local information of the fluid flow and heat transfer inside the gradient metal foam. As far as the author knows, the effect of some key parameters such as Reynolds number on the solid-liquid phase transition in gradient metal foam has not been systematically studied. Therefore, we studied the melting behavior of PCM under a gradient porosity structure at the pore scale.
In this work, a pore-scale LB model was used to study the melting of composite metal foam PCM with porosity gradient. Lattice Boltzmann method (LBM) has been developed as an effective numerical method to solve the problem of solid-liquid phase transition [28,29,30], and has solved the problem of solid-liquid phase transition of PCM in porous media [13,16,31,32]. The remainder of this paper is structured as follows. The governing equations for the melting problem and physical problem are presented in Section 2, followed by the LBM for the enthalpy equation and the fluid flow in Section 3.1 and 3.2. Subsequently, the melting model is verified and boundary treatment in Section 3.3 and Section 3.4, respectively. Numerical results are presented and discussed in Section 4, and finally, some conclusions from the present investigation are summarized in Section 5.

Problem statement and governing equations
In this work, we consider the numerical simulation of solid-liquid phase change in a square cavity equipped filled with a gradient porous media. The schematic diagram of the two-dimensional physical model considered in this study is depicted in Fig.1. A square cavity with side length δ encloses the gradient foam metal, which is filled with solid PCMs at the melting temperature T m . In the solid-liquid phase change, the left wall is raised to a constant temperature T h , higher than T m , the temperature of the right wall is kept at the  melting temperature T m , while the other two walls are assumed to be adiabatic boundary. Additionally, the detailed description of the microstructure of porous media is a necessary prerequisite for the pore scale model. In the present work, the stochastic simulation method of porous media proposed by Chen et al. [33] is employed to generate porous structure. Since the gradient porosity structure created with tri-layer media have almost the same performance as a bilayer media [27], so we adopt a two-layer structure for simplicity. Three structures of porous media have been defined with fixed average porosity as follows: positive gradient structure (porosity increases from 0.8 to 0.95), negative gradient structure (porosity decreases from 0.8 to 0.95), and uniform gradient structure (porosity keep at 0.875). The positive porosity gradient is obtained by increasing the porosity of porous media along the positive x-axis or y-axis direction, while the negative one meant the decreased porosity along the positive x-axis or y-axis direction. Case A and Case B study horizontal gradient porosity and vertical gradient porosity, respectively. Furthermore, It is worth noting that the particle diameter with 7.5 l.u. (lattice unit) and the volume of PCM that can be filled in the gradient structure is the same, which lays the foundation for comparisons of melting behaviors during phase change process.
and the energy equation for the solid phase is: where the subscript l and s denote the liquid and solid phases, respectively. g, u, T, p are the gravity acceleration vector, velocity, temperature, and pressure, respectively. β f is the thermal volumetric expansion coefficient of the fluid. ρ, C p , λ represent the density, heat capacity, and thermal conductivity coefficient, respectively.
Moreover, thermal conductivity coefficient of PCM is: In two phase regions, the enthalpy H can be solved by the enthalpy-based method and is given by here f l denotes the liquid fraction and values within 0 and 1. L represent the latent heat of PCM.
This problem can be characterized by the following three main dimensionless parameters, i.e. Rayleigh number Ra, Prandtl number Pr, Stefan number Ste, which are defined by where ν and α are the kinetic viscosity and thermal diffusivity, respectively

Lattice Boltzmann method for velocity field
In the present study, the incompressible lattice Bhatnagar-Gross-Krook (LBGK) model proposed by Guo et al.
[34] is used to simulate the fluid flow, the evolution equation of the particle velocity field is where f i (x, t) is the probability density distribution functions with velocity eiat position x and time t, ∆t is the time In addition, f eq i is the equilibrium distribution function given by [34] where η i is the model parameter satisfying η 0 = (ω 0 − 1) /c 2 s + ρ 0 , η i = ω i /c 2 s (i 0) with the constant ρ 0 being the fluid average density. ω i is the weight coefficient, represented as The direction of the discrete velocity c i of the model in i direction are given by where i is velocity direction, and c is the lattice speed satisfying c = ∆x/∆t, ∆x is the lattice spacing. The body force where c s = c/ √ 3 is the sound speed of the model, F is the buoyancy force, and can be calculated according to Boussinesq assumption: After a serious of iterative computations, the macroscopic pressure and velocity can be obtained by

Lattice Boltzmann method for temperature field
For the heat transfer in the fluid and solid domains, the double relaxation time LB model proposed by Lu et al.
[36] is adopted to simulate the temperature. It can effectively suppress the non-physical diffusion phenomenon at the solid-liquid interface during the phase transition simulation process. The evolution equation of the temperature distribution function can be described as where g s i (x, t) and g a i (x, t) are the symmetric and anti-symmetric parts of the particle distribution function where the superscript s and a represent the symmetric and anti-symmetric parts of distribution function, respectively. τ g s and τ g a are the symmetric relaxation time and anti-symmetric relaxation time, respectively. The expressions of in which i represents the opposite direction of i. The equilibrium distribution function can be expressed as The dimensionless relaxation time τ s g and τ a g can be obtained by And the total enthalpy can be obtained as follows Finally, the liquid fraction f l and the temperature T can be calculated from total enthalpy, in which H s = C p T s and H l = C p T s + L are total enthalpy of the solid phase and the liquid phase, respectively. T s and T l denote the solidus and liquidus temperatures, respectively.

LBM boundary conditions
In this context, a volumetric LB scheme proposed by Huang et al. is adopted, which can be expressed as where the f * i = f i (x + c i ∆t, t + ∆t) is given by Eq. 8 and u s = 0 is the the velocity of solid phase. What's more, the nonequilibrium extrapolation scheme (NEES) proposed by Guo et al. [37] is applied to deal with all boundary of walls for its second order accuracy. In order to demonstrate the ability of the present code in simulation of fluid flow and heat transfer, we simulated the pure PCM melting at Ra = 25000.0, S te = 0.01, and Pr = 0.02 in a square cavity. Additionally, a grid resolution of 512 × 512 is employed in the simulations, which is fine enough to give the grid-independent solution. We compare the average Nusselt number of left wall, total liquid fraction and melting interface with Huang et.al [38]. Fig.2(a) shows

Validation of phase change model
Nu ave and f l at different Fourier numbers (Fo). It can be clearly seen that the heat transfer performance matches the Huang's results well at different Fo, a good agreement is obtained. Moreover, the comparison of the position of the solid-liquid interface at Fo = 4, Fo = 10 and Fo = 20 is shown in Fig.2(b), and shows good agreement. The present numerical results are in good agreement with the results of Huang's results, indicating that the present LBM model and simulation code are accurate and reliable.

Results and discussion
In this section, numerical simulations are carried out for the solid-liquid phase transition in gradient porous media via the LBM. And we mainly investigate the influence of the geometrical structure and arrangement of gradient porous media on the melting improvement of PCM, which is measured by the nondimensional Fourier numaber (Fo), liquid phase fraction ( f l ). Unless otherwise stated, the Stefan numbe and the Prandtl number are selected as 0.1 and 0.2, which have been studied [38]. In order to show the main physical difference between porous media and PCM, we set the thermal conductivity of porous media to 500 times that of PCM.
We first examine the effect of the gradient porosity on melting behavior at Ra = 10 6 . Fig. 3 illustrates the liquid fraction f l as a function of dimensionless time Fo for Case A and Case B with different gradient structures: uniform, positive and negative gradients. Based on the difference in the slope of the melting curve, the melting process can be divided into three stages. At the early stage, heat is mainly transferred through conduction so that the curves of the three gradient porosity structure completely overlap. With time elapsed, f l incresaes as Fo incresses, and the trend is firstly steep and then gradually weaker, which indicates thermal conductivity performance degradation for melting interface. Specifically, the liquid fraction of the negative gradient in Case A is higher than the other two structures, and this trend is maintained until completely melted, which may due to the fact that the left part with relatively low porosity of negative gradient leads to a growth in natural convection strength.
However, the situation of Case B can get a little complicated, a turning point is observed and the liquid fraction of negative gradient increases with a higher rate compared to the other cases while the melting rate of other cases decreases. Eventually, since the fact that the main heat conduction mechanism has changed from the natural convection to the heat conduction, the melting rate of the three gradients structure is reduced. The negative gradient ends up with a relatively short period of time in Case A and Case B, followed by the uniform and the negative cases.
In order to intuitively understand the flow transition process in different gradient structures. Fig. 4 shows streamline and total liquid fraction distributions with different gradient structure at Ra = 10 6 condition, where  It can be seen from Eqs. 7 that the buoyancy force, as the driving force of natural convection, gradually increases with the Ra increasing. Consequently, the intensity of natural convection became increasing with the increase of the Ra. Finally, the influence of the particle diameters Ra is explored, and the corresponding numerical results are presented in Fig. 11 and Fig. 12. It can be seen from Fig. 11(a) that the liquid fraction cannot be distinguished at the beginning of melting process, which indicates that the effect of particle diameters is not such significant on conduction. About Fo = 3.5 later, the influence of particle diameters appears gradually, which is caused by the decrease of the particle size leads to a growth in the internal surfaces for heat transfer in porous media. However, as Rayleigh number increases further, the situation begins to change and one can observed from Fig. 11(b) and Fig. 12(b) that owing to the increasing of the particle size, the toal melting time of PCM decreases, which is caused by the enhancement of convection. The reason why this phenomenon occurs is that smaller particle diameters indeed enhances the heat transfer by providing larger heat transfer surfaces, but as Rayleigh number increases, natural convection tends to play an increasingly stronger role and smaller particle diameters would lead to lower permeability which suppresses natural convection and then weakens the overall heat transfer. Further, one can deduce that the variation of total melting time depend little on the particle diameters since the  curves in Fig. 11 are almost irregular.

Conclusion
In this paper, the total enthalpy-based lattice Boltzmann model is adopted to simulate the solid-liquid phase change in a square cavity equipped filled with a gradient porous structure. In order to improve computational efficiency, the proposed algorithm is programmed in parallel by using CUDA. The influence of the gradient porosity, Stefan number, gradient direction and Rayleigh number on the heat transfer in the composite are investigated in detail.
According to the present numerical results, it turns out that porous media with gradient porosity have different effects on the solid-liquid phase change process. The positive gradient porosity shows a further reduction for melting time when the Rayleigh number is small. As the Rayleigh numbers continues to increase and when it exceeds a certain critical value, positive gradient showed further enhancement for melting heat transfer, while the negative gradient deteriorated the heat transfer performance. The melting time of the positive gradient increases as the number of Rayleigh increases, while total melting time of the negative gradient increases as the number of Rayleigh increases when the Rayleigh numbers is greater than 1.0 × 10 5 . Furthermore, whether positive gradient or negative gradient, the decreasing of the partical size leads to a growth in effective thermal conductivity, while can create an obstacle to natural convection.

Conflict of interest
We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.