Three-dimensional numerical simulation of melting characteristics of phase change materials embedded with various TPMS skeletons

Abstract

Phase change energy storage technology has great potential for enhancing the efficient conversion and storage of energy. While triply periodic minimal surface (TPMS) structures have shown promise in improving heat transfer, research on their application in phase change heat transfer remains limited. This paper presents numerical simulations of composite phase change materials (PCMs) featuring TPMS skeletons, specifically gyroid, diamond, primitive, and I-graph and wrapped package-graph (I-WP) utilizing the lattice Boltzmann method (LBM). A comparative analysis of the effects of four TPMS skeletons on enhancing the phase change process reveals that the PCM containing the gyroid skeleton melts the fastest, with a complete melting time of 24.1% shorter than that of the PCM containing the I-WP skeleton. The PCM containing the gyroid skeleton is further simulated to explore the effects of the Rayleigh (Ra) number, Prandtl (Pr) number, and Stefan (Ste) number on the melting characteristics. Notably, the complete melting time is reduced by 60.44% when Ra is increased to 10⁶ compared to the case with Ra at 10⁴. Increasing the Pr number accelerates the migration of the mushy zone, resulting in fast melting. Conversely, the convective heat transfer effect from the heating surface decreases as the Ste number increases. The temperature differences caused by the local thermal non-equilibrium (LTNE) effect over time are significant and complex, with peaks becoming more pronounced nearer the heating surface. This study intends to provide theoretical support for the further development of TPMS skeletons in enhancing the phase change process.

Graphical abstract

Keywords

solid–liquid phase change / lattice Boltzmann method (LBM) / triply periodic minimal surface (TPMS) / mushy zone / local thermal non-equilibrium effect (LTNE)

1 Introduction

Today, the world is facing a severe energy crisis. Climate-related reductions in renewable energy production have led to soaring energy prices, posing significant risks to the energy supply security in many countries. Therefore, it is essential to build a stable energy system. Thermal energy storage technology can address mismatches in energy supply and demand regarding time, space, or intensity [1,2]. It can be classified into three categories: sensible heat, latent heat, and thermochemical energy storage [3]. Latent heat storage, in particular, offers high energy storage density and maintains a near-constant temperature in the heat storage and exothermic processes, making it a key focus of research [4]. However, the low thermal conductivity and poor thermophysical stability of phase change materials (PCMs) have constrained the development and application of latent heat energy storage technology [5].

Therefore, advancing methods and techniques to enhance phase change heat transfer has become a major focus of current research. The application of composite PCMs in the thermal management of electronic devices provides significant advantages [6,7]. Li et al. [8], based on lauric acid, using styrene ethylene butylene styrene as the supporting material, and carboxylation multiwalled carbon nanotubes as the high thermal conductivity additive, generated a more cost-effective high thermal conductivity composite PCM. The performance of battery thermal management under high temperature conditions was found to be outstanding. Mo et al. [9] coupled the units-assembled composite PCM module with forced air convection found good insulation, cooling, and heat dissipation performance through experimental tests. Additionally, the energy density of the battery module was increased from 75.6 to 94.4 Wh/kg. Yao et al. [10] designed a novel leak-proof composite PCM by constructing a nano-level polymer framework with high adsorption strength. This design showed outstanding performance in leakage test, mechanical stability assessments, and temperature control evaluations of the battery module. Moreover, numerous studies have demonstrated that the addition of metal fins [11,12], foam metals [13,14], nanoparticles [15], graphite [16], and honeycomb structures [17] could effectively improve the thermal performance of phase change energy storage systems.

Porous metal structures have the advantages of high porosity, large specific surface area, low density, and lightweight [18], which are widely used to enhance the phase change heat transfer process. The melting behavior of the paraffin packaged in a heat sink was experimentally investigated by Feng et al. [19]. The extent and characteristics of the influence of copper foam on the phase change process were analyzed by the effective protection time of the PCM heat sink. A comparative analysis of the foam metal filling ratios for five different porosity conditions was conducted by Li et al. [20] to obtain a filling ratio that maximizes the melting process of PCM, and it was found that the optimum filling ratio decreases with increasing porosity. Moreover, the porosity and pore per inch (PPI) density are important parameters affecting the porous metal structures. Zhuang et al. [21] performed numerical simulations of the heat transfer properties of a metal foam structure with a three-layer gradient porosity in PCM to investigate the melting and heat storage performance under the influence of the gradient change direction. Marri and Balaji [22] combined experiments and numerical simulations to study the thermal performance of composite foam metal PCM heat sinks with gradient porosity and gradient PPI density.

For fluid motion in porous media structures, the size, orientation, and distribution of the pores are extremely complex. Traditional numerical computational methods often encounter challenges such as low parallel computational efficiency and difficulties in resolving complex fluid systems, particularly when describing changes in interface motion [23]. The lattice Boltzmann method (LBM) effectively overcomes this series of difficulties [24]. The LBM is a mesoscopic simulation method that can effectively simulate the flow and heat transfer of the medium within porous structures. He et al. [25] reviewed the research progress of applying the LBM to the solid–liquid phase change process at pore scale and representative elementary volume (REV) scale. Chen et al. [26] simulated the melting process of paraffin in metal foams of two-dimensional at pore scale by LBM to analyze the effects on the phase distribution, temperature, and flow fields. They also verified the results experimentally, with good agreement. Based on the enthalpy-porous medium model, Chen et al. [27] proposed a two-zone composite model of the multiphase flow-porous medium to describe the flow characteristics in the mushy zone. Then, they investigated the flow and heat transfer characteristics in the mushy zone of PCM in detail. Gaedtke et al. [28] proposed an LBM with two-relaxation time collision, which was applied to the melting and conjugate heat transfer processes of metal foam composite PCMs. By varying the combination of three Reynolds (Re) numbers, three porosities, and three PPI densities, Mabrouk et al. [29] simulated the phase change process in a two-dimensional rectangular space containing metal foam to compare and analyze the effects of different parameters on the enhancement of phase change.

The triply periodic minimal surface (TPMS) is a kind of special surface structure that repeats periodically in three independent directions and has zero mean curvature [30]. Additionally, the node coordinates could be described by level-set approximation equations [31]. TPMS structures have good mechanical properties, a high porosity, and a high specific surface area, which have important applications in bone tissue engineering, biological cell growth, heat and mass transfer, reverse osmosis, and ultrafiltration technologies [32–36]. Sadeghi et al. [37] selected eight TPMS structures and developed a model of TPMS metamaterials based on shape memory polymers to assess their thermal, force, and thermomechanical properties at four different volume fractions. Khaleghii et al. [38] numerically modeled the directionality of the elastic properties of 7 TPMS structures across solid phase volume fractions ranging from 10% to 90%, and designed hybrid structures to effectively weaken the anisotropy of elastic modulus compared to the single TPMS structure. Iyer et al. [39] used numerical simulations to evaluate the hydraulic characteristics and thermal performance of 6 TPMS-structured heat exchangers. The results indicated that the TPMS structure enhanced convective heat transfer, significantly reducing the size of heat exchanger. Alteneiji et al. [40] developed three-dimensional models of compact cross-flow heat exchangers featuring gyroid and primitive structures. Their simulations of heat transfer efficiency, convective heat transfer coefficient, and inlet/outlet fluid temperature revealed that the gyroid structure could increase the heat transfer efficiency by 35%. This shows that TPMS structures have significantly superior qualities for enhancing heat transfer.

TPMS structures applied in the phase change process are effective in enhancing the thermal conductivity of PCM [41], which will accelerate the phase change melting process of PCM. Qureshi et al. [42,43] numerically simulated the effects of gyroid, primitive, and I-WP structures compared to the Kelvin structure on the heat transfer characteristics of composite PCM using computational fluid dynamics (CFD). It was found that the ability of TPMS structures to enhance heat transfer was superior to the effect of Kelvin structure to enhance natural convective heat transfer, which could effectively shorten the melting time of PCM. Moreover, the TPMS structures resulted in a better temperature homogeneity of PCM during the melting process under the isothermal flow condition. Qureshi et al. [44] further designed TPMS structures with both positively and negatively graded porosity, maintaining the same average porosity, and compared them to TPMS structures with uniform porosity. The results show that the positively graded structure is more favorable for heat transfer during phase change process for the same structure. Tian et al. [45] combined experiments and simulations to determine the effective thermal conductivity of PCMs embedded with four TPMS structures, as well as those embedded with conventional lattice structures and pure PCMs. Their findings demonstrated that the TPMS structures could significantly accelerate the melting rate of the phase change process. Based on the primitive and I-WP structures, Fan et al. [46] constructed a P-IWP structure (a composite TPMS structure for primitive and I-WP) for improving the localized heat transfer performance of PCMs, which can effectively reduce the temperature rise of the battery when applied to the battery system. The temperature is reduced by 11.3% compared with the battery components with only PCMs. Zhang et al. [47] constructed TPMS structures with three gradient-varying porosities and combined experiments with simulations to compare their effects on the melting performance of PCM. It is found that the skeleton with linear porosity has a higher heat storage rate and that with Boltzmann gradient provides a more homogeneous internal temperature distribution than the uniform skeleton. The studies mentioned above show that the TPMS structures can provide a more stable temperature environment and enhance the heat transfer rate, thereby accelerating the melting process of PCMs. Moreover, the effects of different TPMS structures on the phase change process vary significantly. Therefore, it is important to explore the influence of different TPMS structures on the melting process of PCM.

Currently, the enhancement of the phase change process by adding structures such as foam metals and fins to PCMs has been widely studied. However, research on embedding skeletal structures to enhance the thermal characteristics of PCMs is relatively limited. In particular, TPMS skeletons featuring twisted and complex pore structures, distinguish them from conventional cylindrical or rectangular skeletons. This enables the TPMS skeletons to have a larger heat exchange area at the same porosity. Unlike fins, TPMS structures can enhance thermal conductivity across the entire PCM, rather than being limited to localized areas. Additionally, compared to foam metals, TPMS physical models are easier to generate, and the skeletons can be conveniently produced for experimentation using 3D printing technology. Moreover, the TPMS structures are designed with a high degree of flexibility and parameter controllability. However, there are comparatively few studies on the application of TPMS structures in phase change heat transfer. It is important to further explore the effects on heat transfer and flow characteristics of TPMS structures during the phase change process. In addition, there is a lack of research on the impact of relevant dimensionless parameters on the phase change process and the local thermal non-equilibrium (LTNE) effect during the phase change heat transfer process. In this paper, the solid–liquid phase change processes of composite PCMs embedding various TPMS skeletons are numerically simulated using LBM. Four TPMS skeletal structures with a porosity of 0.9, i.e., gyroid, diamond, primitive, and I-WP are constructed to comparatively analyze the heat transfer and flow characteristics of composite PCMs during the melting process. The composite PCM with a gyroid skeletal structure is selected for further investigation into the influence of dimensionless parameters on the enhancement of the melting process, with a detailed exploration of the LTNE effect. The findings of this research aim to provide theoretical support for the development and application of TPMS structures in enhancing phase change heat transfer.

2 Model building

2.1 Physical model

Four common types of TPMS structures are gyroid, diamond, primitive, and I-WP [48], which can be modeled quickly by implicit trigonometric functions. The higher porosity of the skeleton indicates the smaller volume of the skeleton. At the same volume of the composite PCM, the embedded skeleton with high porosity has the minimal effect on the heat storage capacity of the PCM. Therefore, skeletal structures with a porosity of 0.9 are selected for the simulation investigation.

Fewer studies have focused on the effect of TPMS skeletons on phase change heat transfer and flow under high porosity conditions, highlighting the need for further exploration. The gyroid and diamond skeletons are high-porosity continuous structures that maintain continuity with a porosity of 0.9. The trigonometric functions described in Yoo [49] are used for modeling. In contrast, the primitive and I-WP skeletons already show discontinuities at a porosity above 0.8. To ensure that the skeletons can be generated continuously under the same porosity conditions, optimized trigonometric functions from Zhang [50] are employed for modeling the primitive and I-WP skeletons.

The trigonometric functions for the four TPMS skeletal structures are outlined below.

Gyroid:

(1)\begin{array}{rl}{\Phi }_{\mathrm{G}}=& \sin \left(ax\right)\cos \left(by\right)+\sin \left(by\right)\cos \left(cz\right)+\\ & \sin \left(cz\right)\cos \left(ax\right)-1.17,\end{array}

Diamond:

(2)\begin{array}{rl}{\Phi }_{\mathrm{D}}=& \cos \left(ax\right)\cos \left(by\right)\cos \left(cz\right)-\\ & \sin (ax)\sin (by)\sin (cz)-0.94,\end{array}

Primitive:

(3)\begin{array}{rl}{\Phi }_{\mathrm{P}}=& 10[\cos \left(ax\right)+\cos \left(by\right)+\cos \left(cz\right)]-\\ & 5[\cos (ax)\cos \left(by\right)+\cos (by)\cos \left(cz\right)+\\ & \cos \left(cz\right)\cos \left(ax\right)]-13.5,\end{array}

I-WP:

(4)\begin{array}{rl}{\Phi }_{\text{I-WP}}=& 10[\cos \left(ax\right)\cos \left(by\right)+\cos (by)\cos \left(cz\right)+\\ & \cos \left(cz\right)\cos \left(ax\right)]-5[\cos \left(2ax\right)+\cos \left(2by\right)+\\ & \cos \left(2cz\right)]-13.4,\end{array}

where x, y, and z denote the Cartesian coordinate system; a, b, and c can control the variation of unit pore size in the x, y, and z coordinate directions. Usually enlarging the values of a, b, and c increases the pore density of the skeleton.

Gyroid, diamond, primitive, and I-WP skeletal structures with a porosity of 0.9 are generated using the above formulas. Fig.1 shows the physical models of the cubic cavities of the PCM embedded with skeletons. The blue part is the solid skeletons, and the space outside the skeletons is filled with PCM.

2.2 Mathematical model

The solid–liquid phase change process has an extremely complex heat transfer behavior. Fig.2 illustrates the physical process of the melting pure PCM. Initially, the left surface (xOz plane) of the cubic cavity is suddenly heated and maintained at θ_h = 1. At this point, the heat starts to transfer through the PCM toward the right of the cubic cavity. As the heating continuous, the melted PCM influenced by natural convection, forms a clockwise flow within the cubic cavity, depicted schematically by the white line segment with arrows. As the natural convection is enhanced, the degree of the bending at the phase change interface is increased. During the phase change process, a liquid-phase zone, a solid-phase zone, and a mushy zone are formed in the cubic cavity. The solid and liquid PCM coexist in the mushy zone. The local enlargement of the mushy zone is shown in Fig.2(c), where the black part is the solid PCM while the white area is the liquid PCM. The solid PCM near the liquid-phase zone is small, becoming larger and more extensive as it moves further away from the heating surface. Once the melting process is complete, only the liquid PCM remains, and the internal temperature gradually homogenizes.

This paper is based on the LBM and enthalpy method which can avoid complex iterations in dealing with latent heat source terms. The two-zone composite model of the multiphase flow-porous medium proposed by Chen et al. [27] is used to develop the mathematical model for numerical analysis. To simplify the model, it is assumed that the thermal properties of the PCM are to remain constant during the phase change. In addition, the liquid PCM is regarded as a Newtonian incompressible fluid and the density change satisfy the Boussinesq assumption. Moreover, the liquid PCM flows as a non-constant laminar flow in the cubic cavity and the viscous dissipation during the flow is neglected. Furthermore, the volume change during the melting process of PCM is neglected.

The dual temperature distribution is used to establish the energy equations. Based on the above-mentioned assumptions, the control equation [51] in the phase change process are established.

Continuity equation:

(5)\mathrm{\nabla }\cdot u=0.

Momentum equation:

(6)\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+\left(u\cdot \mathrm{\nabla }\right)\frac{u}{\gamma }=-\frac{1}{\rho {}_{\mathrm{f}\mathrm{l}}}\mathrm{\nabla }\left(\gamma p\right)+{\nu }_{\mathrm{e}\mathrm{f}\mathrm{f}}{\mathrm{\nabla }}^{2}u +F.

Energy equation for PCM:

(7)\begin{array}{rl}& \left[\gamma {\left(\rho c{}_{\mathrm{p}}\right)}_{\mathrm{f}\mathrm{l}}+\left(1-\gamma \right){\left(\rho c{}_{\mathrm{p}}\right)}_{\mathrm{f}\mathrm{s}}\right]\frac{\mathrm{\partial }{T}_{\mathrm{f}}}{\mathrm{\partial }t}+{\left(\rho c{}_{\mathrm{p}}\right)}_{\mathrm{f}\mathrm{l}}\left(u\cdot \mathrm{\nabla }{T}_{\mathrm{f}}\right)=\\ & k_{\mathrm{e}\mathrm{f}\mathrm{f}}{\mathrm{\nabla }}^2{T}_{\mathrm{f}}-\gamma {\rho }_{\mathrm{f}\mathrm{l}}La\frac{\mathrm{\partial }\gamma }{\mathrm{\partial }t}.\end{array}

Energy equation for solid skeleton:

(8){\left(\rho c_{\mathrm{p}}\right)}_{\mathrm{s}}\frac{\mathrm{\partial }{T}_{\mathrm{s}}}{\mathrm{\partial }t}=k_{\mathrm{s}}{\mathrm{\nabla }}^2{T}_{\mathrm{s}},

where γ is the liquid fraction, indicating the ratio at which the volume of liquid PCM accounts for the total volume of PCM, k_s is the ratio of thermal conductivity of solid skeleton, and F is the external source term, expressed as

(9)F=-\frac{\gamma {\nu }_{\mathrm{f}\mathrm{l}}}{K}u-\frac{\gamma F_{\epsilon }}{\sqrt{K}}u|u|+\gamma g\beta (T- T_{\mathrm{r}\mathrm{e}\mathrm{f}}),

where K is the permeability of the porous skeleton, F_ε is the shape factor of the porous skeleton, and T_{\text{ref}} is the reference temperature.

The porous skeleton permeability is expressed as

(10)K=\frac{{\gamma }^3}{C(1-\gamma {)}^2} ,

where C is the constant of the mushy zone. This paper sets C = 1600 according to Voller and Prakash [52].

At γ = 0, no liquid PCM exists, indicating that melting has not yet started. When γ = 1, the PCM has completely melted and exists entirely in the liquid state. When γ is between 0 and 1, the PCM absorbed the heat gained as latent heat. At this stage, solid and liquid coexist in a specific region known as the mushy zone. The two-zone composite model divides the mushy zone into high and low liquid-content zones based on the liquid fraction γ = 0.7. The Brinkmann-Forchheimer-Darcy model is applied in the low liquid-content zone (0 < γ < 0.7), while the multiphase flow model is utilized in the high liquid-content zone (0.7 < γ < 1).

In the low liquid-content zone, the effective thermal conductivity is expressed as

(11)k_{\mathrm{e}\mathrm{f}\mathrm{f}}=\left(1-\gamma \right)k_{\mathrm{f}\mathrm{s}}+\gamma k_{\mathrm{f}\mathrm{l}}.

In the high liquid-content zone, the effective thermal conductivity is expressed [53] as

(12)k_{\mathrm{e}\mathrm{f}\mathrm{f}}=k_{\mathrm{f}\mathrm{l}}\cdot {\displaystyle \frac{2+{\displaystyle \frac{k_{\mathrm{f}\mathrm{s}}}{k_{\mathrm{f}\mathrm{l}}}}+2(1-\gamma )\left({\displaystyle \frac{k_{\mathrm{f}\mathrm{s}}}{k_{\mathrm{f}\mathrm{l}}}}-1\right)}{2+{\displaystyle \frac{k_{\mathrm{f}\mathrm{s}}}{k_{\mathrm{f}\mathrm{l}}}}-(1-\gamma )\left({\displaystyle \frac{k_{\mathrm{f}\mathrm{s}}}{k_{\mathrm{f}\mathrm{l}}}}-1\right)}}.

The expression for the characterized kinematic viscosity [54] is given by

(13)\frac{{\nu }_{\mathrm{e}\mathrm{f}\mathrm{f}}}{{\nu }_{\mathrm{f}}}={\left(\gamma -1.16(1-\gamma {)}^2\right)}^{-2.5}.

The liquid fraction of PCM is denoted by

(14)En=( c_{p}T{)}_{\mathrm{f}}+\gamma La,

where En is the enthalpy of the PCM and the relationship with the liquid PCM temperature T_fl is as follows [55]:

(15)

where E{n}_{\mathrm{f}\mathrm{s}} is the enthalpy corresponding to the phase change onset temperature T_{\mathrm{f}\mathrm{s}}, E{n}_{\mathrm{f}\mathrm{l}} is the enthalpy corresponding to the phase change termination temperature T_{\mathrm{f}\mathrm{l}}, and T_{\mathrm{R}} is the phase change radius, expressed as

(16)T_{\mathrm{R}}=\frac{T_{\mathrm{f}\mathrm{l}}-T_{\mathrm{f}\mathrm{s}}}{2}.

To reduce the number of variables in the simulation and to derive more general laws, dimensionless parameters are introduced to transform the control equations into dimensionless form.

\begin{array}{rl}& X=\frac{x}{L},Y=\frac{y}{L},Z=\frac{z}{L},U=\frac{uL}{{\alpha }_{\mathrm{f}\mathrm{l}}},\theta =\frac{T-T_{\mathrm{c}}}{T_{\mathrm{h}}-T_{\mathrm{c}}},\\ & P=\frac{p{L}^2}{{\rho }_{\mathrm{f}\mathrm{l}}{\nu }_{\mathrm{f}\mathrm{l}}{\alpha }_{\mathrm{f}\mathrm{l}}},Fo=\frac{t{\alpha }_{\mathrm{f}\mathrm{l}}}{L^2},Pr=\frac{{\nu }_{\mathrm{f}\mathrm{l}}}{{\alpha }_{\mathrm{f}\mathrm{l}}},\end{array}

\begin{array}{rl}& Ra=\frac{g\beta (T_{\mathrm{h}}-T_{\mathrm{c}})L^3}{{\nu }_{\mathrm{f}\mathrm{l}}{\alpha }_{\mathrm{f}\mathrm{l}}},Ste=\frac{c_{p_{\mathrm{f}\mathrm{l}}}(T_{\mathrm{h}}-T_{\mathrm{c}})}{La},\\ & \sigma =\frac{\rho c_p}{{(\rho c_p)}_{\mathrm{f}\mathrm{l}}},R=\frac{k}{k_{\mathrm{f}\mathrm{l}}},Da=\frac{K}{L^2},\end{array}

where L is the feature length; and Fo, Pr, Ra, Ste, and Da are the Fourier number, Prandtl number, Rayleigh number, Stefan number, and Darcy number, correspondingly. Fo number can represent the dimensionless time of the phase change. The physical relationship between the PCM and the skeleton material is determined by the thermal conductivity ratio R_{\mathrm{s}}= k_{\mathrm{s}}/k_{\mathrm{f}\mathrm{l}}.

The dimensionless control equations are obtained as follows.

Continuity equation:

(17)\mathrm{\nabla }\cdot U=0.

Momentum equation:

(18)\begin{array}{rl}\frac{1}{Pr}& \left[\frac{1}{\gamma }\frac{\mathrm{\partial }U}{\mathrm{\partial }Fo}+\frac{1}{{\gamma }^2}\left(U\cdot \mathrm{\nabla }\right)U\right]=-\mathrm{\nabla }P+\frac{1}{\gamma }{\mathrm{\nabla }}^{2}U-\\ & \left(\frac{1}{Da}+\frac{F_{\epsilon }}{Pr\sqrt{Da}}\left|U\right|\right)U+Ra{\theta }_{\mathrm{f}\mathrm{l}}\frac{g}{\left|g\right|}.\end{array}

Energy equation for PCM:

(19)\sigma \frac{\mathrm{\partial }\theta }{\mathrm{\partial }Fo}+U\cdot \mathrm{\nabla }\theta =R_{\mathrm{f}}{\mathrm{\nabla }}^2{\theta }_{\mathrm{f}}-\frac{1}{Ste}\frac{\mathrm{\partial }\gamma }{\mathrm{\partial }Fo}.

Energy equation for solid skeleton:

(20){\sigma }_{\mathrm{s}}\frac{\mathrm{\partial }{\theta }_{\mathrm{s}}}{\mathrm{\partial }Fo}=R_{\mathrm{s}}{\mathrm{\nabla }}^2{\theta }_{\mathrm{s}}.

2.3 LBM

In this paper, the LBM is applied, utilizing the D3Q19 model proposed by Qian and Zhou [56]. Double-distribution functions are employed to describe the velocity and temperature fields during the phase change process. The boundary conditions are handled using the non-equilibrium extrapolation method [57]. The following equations [58] are formulated to study the heat transfer and flow in the solid–liquid phase change process of the PCM containing porous skeletons.

The evolution equation of the velocity distribution function is given by

(21)f_i\left(r+e_i\delta t,t+\delta t\right)-f_i\left(r,t\right)=-\frac{1}{{\tau }_{\mathrm{f}}} \left[f_i\left(r,t\right)-f_i^{\mathrm{e}\mathrm{q}}\left(r,t\right)\right]+\delta t{F}_i,

whose equilibrium distribution function is expressed as

(22)f_i^{\mathrm{e}\mathrm{q}}=\rho {\omega }_i\left[1+\frac{e_i\cdot u}{c{}_{\mathrm{s}}^2}+\frac{{\left(e_i\cdot u\right)}^2}{2c{}_{\mathrm{s}}^4}-\frac{u^2}{2c{}_{\mathrm{s}}^2}\right],

where {\omega }_i is the weighting factor, e_i is the the lattice speed, and c_s is the lattice sound speed.

The external force source term is expressed as

(23)F_i=\left(1-\frac{1}{2{\tau }_{\mathrm{f}}}\right)\rho {\omega }_i\left[\frac{e_i\cdot F}{c_{\mathrm{s}}^2}+\frac{uF:(e_i{e}_i-c_{\mathrm{s}}^{2}I)}{\gamma c_{\mathrm{s}}^4}\right],

The evolution equation of the temperature distribution function of the PCM is given by

(24)\begin{array}{rl}& g_{i,\mathrm{f}}\left(r+e_i\delta t,t+\delta t\right)-g_{i,\mathrm{f}}\left(r,t\right)=\\ -& \frac{1}{{\tau }_{T,\mathrm{f}}}\left[g_{i,\mathrm{f}}\left(r,t\right)-g_{i,\mathrm{f}}^{\mathrm{e}\mathrm{q}}\left(r,t\right)\right]+\delta tS{r}_{i,\mathrm{f}},\end{array}

whose equilibrium distribution function is expressed as

(25)g_{i,\mathrm{f}}^{\mathrm{e}\mathrm{q}}= T_{\mathrm{f}}{\omega }_i \left[1+\frac{e_i\cdot u}{c_{\mathrm{s}}^2}+\frac{{\left(e_i\cdot u\right)}^2}{2c_{\mathrm{s}}^4}-\frac{u^2}{2c_{\mathrm{s}}^2}\right],

The discrete source term is expressed as

(26)S{r}_i= -{\omega }_i \frac{La}{c_{p_{\mathrm{f}\mathrm{l}}}}\frac{{\gamma }_{\mathrm{f}\mathrm{l}}\left(t+\delta t\right)-{\gamma }_{\mathrm{f}\mathrm{l}}\left(t\right)}{\delta t}\left[1+\left(1-\frac{1}{2{\tau }_T}\right)\frac{e_i\cdot u}{c_{\mathrm{s}}^2}\right],

The evolution equation of the temperature distribution function of the solid skeleton is given by

(27)g_{i,\mathrm{s}}\left(r+e_i\delta t,t+\delta t\right)- g_{i,\mathrm{s}}\left(r,t\right)=-\frac{1}{{\tau }_{T,\mathrm{s}}}\left[g_{i,\mathrm{s}}\left(r,t\right)-g_{i,\mathrm{s}}^{\mathrm{e}\mathrm{q}}\left(r,t\right)\right],

whose equilibrium distribution function is expressed as

(28)g_{i,\mathrm{s}}^{\text{eq}}=T_{\mathrm{s}}{\omega }_i,

where {\omega }_i is the weight coefficient as a function of particle velocity.

The dimensionless relaxation time expressions are given by

(29){\tau }_{\mathrm{f}}={\nu }_{\mathrm{f}\mathrm{l}}\frac{1}{c_{\mathrm{s}}^2\delta t}+0.5,

(30){\tau }_{T,\mathrm{f}}={\alpha }_f\frac{1}{c_{\mathrm{s}}^2\delta t}+0.5,

(31){\tau }_{T,\mathrm{s}}={\alpha }_{\mathrm{s}}\frac{1}{c_{\mathrm{s}}^2\delta t}+0.5,

where {\tau }_{\text{f}} is the velocity dimensionless relaxation time and {\tau }_T is the temperature dimensionless relaxation time.

2.4 Conditions for numerical simulation

In this paper, a dimensionless representation is utilized and the simulation space consists of cubic cavities with a dimensionless length of 1. By applying dimensionless treatment of the parameters, the resulting findings can be generalized for broader applications. Fig.3 shows a three-dimensional physical model of PCM embedded with the gyroid skeleton. The yOz plane of the cubic cavity is the heating surface, and the dimensionless temperature is set to θ_h = 1. The rest of the walls are adiabatic. The red plane indicates the heating surface while the blue planes represent the adiabatic surfaces. The phase change cubic cavity is filled with solid skeleton and PCM. The initial temperature of the phase change process is θ₀ = 0. The thermal conductivity ratio of the solid PCM to the liquid PCM is set to be k_fs = 2, while the thermal conductivity ratio of the skeleton and the liquid PCM is set to be k_s = 100. The PCM has a phase change center temperature of θ_c = 0.2 and a phase change radius of θ_R = 0.1. This means that the PCM starts to melt when the phase change temperature is θ_fs = 0.1. When the temperature reaches θ_fl = 0.3, the PCM completely melts and exists in liquid form. When the temperature is between 0.1 and 0.3, the PCM absorbs heat and stores it as latent heat of phase change. In this range, the PCM is in the phase change mushy zone, where both liquid and solid states coexist.

3 Model validation

3.1 Mesh independence verification

The quantity of mesh division for the model directly affects the accuracy of the calculation results, making it essential to select a suitable mesh criterion. A tetrahedral meshing pattern is employed for the models. The PCM embedded in the gyroid skeleton with a porosity of 0.9 is divided into three mesh levels of 60 × 60 × 60, 80 × 80 × 80 and 100 × 100 × 100, respectively. As shown in Fig.4, the Fo number for the complete melting of PCM obtained from the simulation at the three mesh levels are 0.165, 0.163, and 0.162. Therefore, the model used in this study has a good mesh independence. In comparison to the 80 × 80 × 80 mesh, the difference in complete melting time for the PCM is 1.2% with the 60 × 60 × 60 mesh and 0.6% with the 100 × 100 × 100 mesh. Additionally, the calculation time for the 80 × 80 × 80 mesh is shorter. To ensure a high computational accuracy while maintaining a better numerical simulation efficiency, the 80 × 80 × 80 mesh is selected for the working condition models in this paper.

3.2 Validation of phase change process for pure PCM

The phase change processes are often accompanied by complex and nonlinear natural convection phenomena. These processes are often influenced by various aspects such as initial conditions, boundary treatments, modeling methods, and model shapes during simulation. For these reasons, the model used in this paper must be validated with relevant experiments.

To properly validate the numerical model, the experiment of phase change of lauric acid in a rectangular cavity in Kamkari et al. [59] is referred to. The physical parameters of lauric acid are normalized and the dimensionless parameters obtained are input into the model to define the material properties of the PCM. A rectangular structure with lattice counts of 50 × 120 × 120 is constructed. The corresponding boundary and initial conditions are established through dimensionless quantization of experimental conditions from Kamkari et al. [59]. The right wall serves as a thermostatically heating surface. The simulation results are converted into quantitative data and compared with the experimental findings in Kamkari et al. [59]. Fig.5 presents a comparison of liquid fraction from numerical simulation and experimental results, with a maximum error of 5.33%. Therefore, the results of the current model align well with the experimental findings. This indicates that the dimensionless lattice Boltzmann model employed in this paper effectively captures the solid–liquid phase change process.

3.3 Verification of fluid-solid coupling

Compared to the phase change process of pure PCM, the PCM embedded with porous skeletons exhibits more complex fluid–solid coupling phenomena during the melting process. A square cavity with side length L/2 nested in the middle of a square cavity with side length L is used as a validation model. The parameters of the simulation are set to Pr = 0.71 and Ra = 10³, 10⁴, 10⁵, and 10⁶. The results of the simulation are compared with the findings in Das and Reddy [60], as shown in Fig.6. It can be found that the numerical simulation results are well consistent with the findings in Das and Reddy [60], with a maximum error of 2.03%. Therefore, the theoretical model developed in this paper accurately represents the fluid-solid coupling between the solid skeleton and the PCM.

4 Results and discussion

In this paper, the above-mentioned dimensionless lattice Boltzmann model is employed to study the melting process of the PCMs containing TPMS skeletons. The distributions and characteristics of the phase change interface, velocity, and temperature are compared in composite PCMs embedded with different TPMS skeletons. The phase change cubic cavity containing the gyroid skeleton is selected to study the influence of dimensionless parameters Ra, Ste, and Pr numbers on the melting process. The LTNE in the horizontal direction is analyzed to investigate the differences in heat transfer between PCM and TPMS skeletons.

4.1 Influence of different TPMS skeletons on phase change characteristics

The phase change interfaces are visual representations to characterize the melting extent of the PCM. The parameters are set to Ra = 10⁵, Ste = 5, Pr = 1, and Rs = 100 to perform numerical simulations for pure PCM and composite PCMs embedding four TPMS skeletons. Fig.7 shows the migrations of phase change interface for pure PCM and PCM containing four different TPMS skeletons. The red area represents the fully melted zone while the blue area is the unmelted zone. The colored area between red and blue is the mushy zone. In the cubic cavity, the red surface defined by the dashed line is the phase change interface between the liquid-phase zone and the mushy zone, while the blue surface defined by the dashed line is the phase change interface between the mushy zone and the solid-phase zone. Looking laterally at Fig.7, the phase change interfaces are almost parallel to the heating surface at the dimensionless time of Fo = 0.01 in all cases. This indicates that heat is transferred uniformly from the heating surface to the interior of the cubic cavities, with thermal conduction playing a dominant role. As time progresses, the mushy zone tilts, which can be attributed to the natural convection of the liquid PCM driven by thermal buoyancy. The high-temperature liquid rises to the upper part of the cavities, and this heat convergence causes the upper PCM to melt more rapidly. As a result, the mushy zone in the upper part of the cubic cavities migrates faster. With increased natural convection, the tendency of the mushy zone to bend to the lower right corner of the cubic cavities becomes more pronounced. Consequently, the upper part of the mushy zone becomes narrower, while the lower part widens.

The migrations of the phase change interface with dimensionless time in different cases have significant features. A longitudinal comparison of Fig.7 shows that the melting process of the PCMs embedded with skeletal structures is faster than that of the pure PCM. The phase change interface of PCMs containing TPMS skeletons is not smooth, but forms prominences along the skeletons. At Fo = 0.05, the mushy zone of the PCM containing the gyroid skeleton is the widest. This indicates that more heat is absorbed by the PCM so that there is more PCM in a molten state within the mushy zone. The phase change interfaces of the PCM containing the gyroid skeleton migrate the fastest at every moment. With the same skeleton volume, it can be seen that the gyroid skeleton is the most favorable for the heat transfer development of the phase change. The second is the PCM containing the diamond skeleton. Compared to the primitive and I-WP skeletons, the gyroid and diamond skeletons exhibit more sparse distributions, resulting in less flow obstruction for the liquid PCM. The natural convection of liquid PCM enhances convective heat transfer, which accelerates the melting process of the PCM, causing the phase change interface migrates faster. In comparison to the primitive skeleton, the internal pore distribution of the I-WP skeleton is more complex and denser, which creates a greater resistance to flow and is not conducive to heat transfer. As a result, the PCM containing the I-WP skeleton melts the slowest. The weak junctions of the cells in the diamond skeleton and the large center volume of the cells increase the degree of non-uniformity in the distribution of the skeleton volume in the cubic cavity. The enhancement of thermal conductivity is weaker due to the lower volume of the skeleton at the weak joints. On the other hand, the gyroid skeleton has a more uniform density distribution in the cubic cavity, which is beneficial to enhance the overall thermal conductivity of the PCM and has a faster melting rate. The effectiveness of skeletal structures for phase change is ranked from highest to lowest as follows: gyroid, diamond, primitive, and I-WP.

The solid–liquid phase change produces an interaction between the velocity and temperature fields. Two X−Z planes perpendicular to the heating surface are selected in the cubic cavity at Y = 0.3 and 0.9 to explore this relationship. Fig.8 shows the velocity and temperature fields of the phase change for PCM containing different TPMS skeletons. The dark gray irregular closed figures are the skeleton cross-sections. The color gradient from red to blue indicates a temperature range from high to low. The lines between each color in the X−Z planes are isotherms with dimensionless temperatures labeled on them. The blue line segments with arrows denote the velocity vectors in the X−Z planes. The temperature cloud maps and isotherms illustrate the temperature field while the velocity vector lines depict the velocity field. The isotherms at the skeleton cross-sections distributed in the X−Z planes are more prominent. The reason for the bumpy shape of the isotherms is that the thermal conductivity of the skeleton is higher than that of the PCM, and the heat is transferred faster through the skeleton. For the same cubic cavity, the distributions of isotherms and velocity vectors are significantly different for the two X−Z planes with different shapes of skeleton cross-sections. For example, the velocity vectors of the PCM containing the gyroid skeleton at the X−Z plane with Y = 0.3 is significantly larger than those at the X−Z plane with Y = 0.9, and the temperature change occurs more rapidly. This difference is attributed to the greater distribution of the gyroid skeleton in the Y = 0.3 plane, which enhances local thermal conductivity and accelerates the heat conduction process, leading to a more extensive heat transfer. In the PCM containing the primitive skeleton, the excessive distribution of the skeleton in the Y = 0.9 plane, occupying almost one-half of the plane area, creates a great resistance to the flow of the liquid PCM. It significantly restricts the full development of natural convection, resulting in a weak velocity vector distribution. This indicates that the skeletal structure has a remarkable influence on the temperature and velocity fields during the phase change process.

At Fo = 0.01, the phase change process has just started, with thermal conductivity being the dominant factor; as a result, the isotherms are almost parallel to the heating surface. At this stage, the liquid PCM is minimal, and only a faint velocity distribution exists at the heating surface. As the phase change process progresses, the temperature and density differences between the melted PCM and the unmelted PCM lead to natural convection in the liquid PCM. The velocity vectors increase in magnitude and form clockwise circulations. As a result, the isotherms in the upper part of the cubic cavity migrate more rapidly, resulting in a steeper temperature gradient compared to the lower part of the cubic cavity, especially at Fo = 0.10 and 0.15. Moreover, the velocity is more intense at the heating surface and near the phase change interface between the liquid-phase zone and mushy zone.

Comparing different cases at the same Fo, the high-temperature zones in the X-Z planes of the PCMs containing the primitive and I-WP skeletons are much smaller than those of the other two cases. This suggests that the primitive and I-WP skeletons are less effective in enhancing heat transfer. The cross section of the primitive skeleton in the Y = 0.9 plane is significantly large, while the distribution of the skeleton in the Y = 0.3 plane is sparse. Similarly, there are also few distributions of I-WP skeleton in the selected planes. This results in a gentle flow with weak perturbations as the liquid PCM flows through the skeleton. The distribution of the gyroid and diamond skeletons is more pronounced in the Y = 0.3 and 0.9 planes. This enhanced distribution promotes natural convection driven not only by temperature differences but also by the stronger perturbation of the skeleton. As a result, convective heat transfer is further developed, accelerating the melting process. The PCM containing the gyroid skeleton has a wider range of natural convection, leading to a greater volume of liquid PCM and an increased melting rate. Additionally, its temperature field evolves more quickly, facilitating faster heat transfer.

The liquid fraction γ is an important indicator of the degree of phase change. The tangent slope at any point on the curve in Fig.9 reflects the melting rate at that moment. The PCM containing the gyroid skeleton has the fastest melting rate, reaching the final melting state first. There is no significant difference in the melting rate of the PCMs containing primitive and I-WP skeletons until Fo = 0.041. However, beyond this point, the PCM containing the primitive skeleton begins to melt faster. The analysis indicates that the contact areas of the primitive and I-WP skeletons with the heating surface are essentially equal, leading to similar effect on the heat transfer during the initial stage of the phase change. However, the different internal structures of the skeletons have an impact on the natural convection of the liquid PCM. While the skeletal structures enhance thermal conductivity, they also impede the development of natural convection. This result in varying melting rates among the four PCMs containing the skeletons. Notably, the TPMS skeletons significantly facilitate the melting process of the PCM compared to pure PCM.

As can be seen in Fig.10, the complete melting dimensionless time Fo of the pure PCM is 0.25. Under the same conditions, the composite PCMs containing gyroid, diamond, primitive, and I-WP skeletal structures show reductions in complete melting time of 35.6%, 29.6%, 22.4%, and 15.2%, respectively. The gyroid skeleton, with optimal enhancement, results in complete melting of the PCM 24.1% faster than the I-WP skeleton.

4.2 Effect of Ra, Pr, and Ste numbers on melting characteristics of PCM containing gyroid skeleton

The above results indicate that the gyroid skeleton has the best enhancement in phase change heat transfer. In this section, the PCM containing the gyroid skeleton will be selected to investigate the influence of different dimensionless parameters on the heat transfer and flow characteristics during the melting process.

4.2.1 Influence of Ra number on melting process of composite PCM

The Ra number characterizes the strength of natural convection. As the Ra number increases, the intensity of natural convection in the liquid PCM becomes greater. Under the same conditions of Ste = 5, Pr = 1, and Rs = 100. Ra = 10⁴, Ra = 10⁵, and Ra = 10⁶ are selected for comparative analysis. Fig.11 illustrates the variation of the liquid fraction during the phase change with different Ra numbers. The complete melting dimensionless times of the composite PCMs for Ra = 10⁴, Ra = 10⁵, and Ra = 10⁶ are 0.225, 0.163, and 0.089, respectively. Increasing the Ra number from 10⁴ to 10⁵ reduces the complete melting time by 27.56%, while the extending the Ra number from 10⁵ to 10⁶ decreases it by 45.40%. Moreover, when the Ra number is increased from 10⁴ to 10⁶, the complete melting of the PCM precedes by 60.44%. This enhancement is attributed to the increase in natural convection, which intensifies convective heat transfer in the cubic cavity, allowing the PCM to absorb more heat rapidly to facilitate the melting process.

A transversal line perpendicular to the center of the heating surface and parallel to the x-axis is selected at the Z/2 height. Uz is the velocity component in the z-direction, with upward movement specified as positive. The distribution of the Uz varying along the X position is investigated at Fo = 0.05. It is observed in Fig.12 that Uz first increases to a positive peak, then decreases to a negative peak, and finally converges to zero. At Fo = 0.05, the PCM near the heating surface has melted. The high-temperature liquid PCM flows upward in the cubic cavity due to the thermal buoyancy force, resulting in a positive Uz. The dimensionless positive peaks of velocity for Ra = 10⁴, 10⁵, and 10⁶ are 0.00487, 0.01169, and 0.01634, respectively. Compared to the Ra = 10⁴ and 10⁵, the positive velocity peak at the Z/2 height is increased by a factor of 1.40 and 0.40 when the Ra = 10⁶. In the low liquid-content zone, the temperature is lower than that of the liquid PCM on the left, leading to a density difference. As a result, the liquid PCM moves downward under the force of gravity, resulting in a negative velocity peak. The dimensionless negative peaks for velocity of Ra = 10⁴, 10⁵, and 10⁶ are −0.00234, −0.00768, and −0.01321, respectively, which are lower than the positive velocity peaks. The temperature and density differences within the liquid-phase zone are minimal, resulting in a relatively low value of Uz. The position where Uz finally stabilizes at zero is closer to the right side of the cubic cavity, indicating a larger liquid-phase zone. In the case of Ra = 10⁶, Uz at the height of Z/2 reaches zero at X = 0.785, which implies a more rapid evolution of the phase change interface and the fastest melting of the PCM.

4.2.2 Influence of Pr number on melting process of composite PCM

The Pr number indicates the relative importance of momentum diffusion to thermal diffusion, reflecting the relationship between the thickness of the boundary layers of flow and temperature. For the numerical simulations, the parameters are set to Ste = 5, Ra = 10⁵, and Rs = 100, with Pr values of 0.5, 1, and 10 selected for analysis. Fig.13 shows the variation of liquid fraction for phase change with different Pr numbers. A higher Pr number is more favorable to the melting process of the PCM. The complete melting dimensionless times for Pr = 0.5, 1, and 10 are 0.174, 0.162, and 0.126, respectively. It is worth noting that the Pr number increases by a factor of 20 from 0.5 to 10, resulting in a 27.59% reduction in the complete melting dimensionless time. The difference in the liquid fraction is relatively slight for small Pr numbers of 0.5 and 1. It is evident that as the Pr number increases, the enhancement in heat transfer becomes more pronounced.

A transversal line, perpendicular to the center of the heating surface and parallel to the x-axis, is selected at the height of Z/2. The variation in the width of the mushy zone along this line, as shown in Fig.14, indirectly indicates the melting rate of the PCM. A wider mushy zone means a faster heat transfer. The PCM absorbs more heat, which facilitates the melting process. At the onset of the phase change, the PCM first stores heat as latent heat, leading to the formation of the mushy zone. Once the phase change temperature is reached, liquid PCM begins to appear. With the phase change process proceeding, more heat is transferred to the interior of the cubic cavity through thermal conductivity and natural convection. The width of the mush zone gradually expands, ultimately reaching a similar peak. The dimensionless times to peak are 0.067, 0.06, and 0.045 for Pr = 0.5, 1, and 10, respectively. The peak disappears when the right boundary of the mush zone makes contact with the right wall at the height of Z/2, after which the width of the mushy zone begins to decrease. Subsequently, the left boundary of the mush zone continues to approach the right wall until the width reaches zero. At this point, all of the PCM at the Z/2 height melted. It is obvious that the width of the mushy zone peaks earliest for Pr = 10, indicating that the phase change interfaces move the fastest. Consequently, the solid-phase zone at the height of Z/2 is the first to disappear, allowing the melting process to be completed more quickly for Pr = 10.

4.2.3 Influence of Ste number on melting process of composite PCM

The Ste number indicates the ratio of sensible heat to the latent heat in the PCM, playing an important role in the evolution of the mushy zone and the heat storage capacity. For the numerical simulations, the parameters are set to Pr = 1, Ra = 10⁵, and Rs = 100, with Ste values of 1, 5, and 10, selected. The variation of the liquid fraction during phase change for different Ste numbers is shown in Fig.15. The melting rate of the PCM increases at a higher Ste number and the complete melting dimensionless times for Ste = 5 and 10 are 0.162 and 0.137, respectively. At this point, the liquid fraction for Ste = 1 has not yet reached 0.8, with a complete melting dimensionless time of 0.336. When the Ste number increases to 5 and 10 from 1, the melting process completes 51.79% and 59.23% earlier, respectively, compared to Ste = 1. This demonstrates that increasing the Ste number from a lower value can effectively accelerate the melting process. However, with a further increase in the Ste number, the melting rate grows at a reduced degree.

Fig.16 shows the variation of Nu_avg number on the heating surface during phase change with different Ste numbers. A larger Ste number means an increase in sensible heat, resulting in less heat required for the phase change process, which accelerates melting. As the PCM near the heating surface melts, heat transfer occurs at a faster rate, leading to a more rapid decrease in the temperature gradient at the heating surface. In contrast, the temperature gradient at the heating surface is relatively larger for a smaller Ste number, resulting in a higher Nu_avg value during leveling-off and a stronger convective heat transfer effect. The difference in Nu_avg number between Ste = 5 and 10 is not significant, indicating that larger Ste numbers have a limited influence on the convective heat transfer at the heating surface.

4.3 LTNE effect

The principle of adding skeletal structures to the PCM to enhance phase change heat transfer is to improve the overall thermal performance of the PCM by using high thermal conductivity skeletal materials. The heat can be transferred quickly through the solid skeletons while accelerating the melting of the PCM around the skeletons. However, there are differences in heat transfer between the skeleton and the PCM, especially for metal skeletons, where the thermal conductivity is much higher than that of the PCM. The skeleton transfers heat preferentially through itself, resulting in less heat being transferred to the surrounding PCM. Moreover, the PCM absorbs heat from the heating surface relatively slowly. As a result, a non-negligible temperature difference (Δθ = θ_s− θ_f) occurs between the skeleton surface and the immediately adjacent PCM which is used to represent the LTNE effect. Due to the notable temperature decrease that occurs in thin layer from the surface of the skeleton to the PCM, the heat transfer to PCM is impeded, affecting the rate of heat absorption and exothermicity of the PCM during the phase change process. This will prolong the thermal response time of the PCM and reduce the overall heat transfer rate. Moreover, the high degree of non-uniformity of heat transfer between the PCM and the high-thermal-conductivity skeleton can affect the overall thermal stability of the composite PCM.

The PCM containing the gyroid skeleton is selected to investigate the LTNE effect. The parameters are set to Ra = 10⁵, Pr = 1, Ste = 5, and Rs = 100. Investigation points for the LTNE effect are selected based on the special structural distribution of the TPMS skeleton. Four groups of investigation points are positioned perpendicular to the heating surface at a height of 0.525Z, labeled L1, L2, L3, and L4 in the x-direction, with L1 being close to the heating surface. As shown in Fig.17, each group consists of one point on the skeleton surface and another point in the PCM, with a distance of three grid units between the two points. The temperature difference at the investigation points in the same group is used to reflect the LTNE effect.

The variation of temperature difference with dimensionless time for each group of investigation points is shown in Fig.18. Each curve shows a distinct positive peak. The dimensionless times for the appearance of positive peaks in the L1, L2, L3, and L4 cases are 0.02055, 0.08064, 0.117, and 0.12017, respectively. The corresponding dimensionless peak values are 0.0375, 0.01082, 0.00619, and 0.00383. It can be noticed that the peak of the temperature difference for the groups further away from the heating surface occurs later and gradually decreases. The thermal conductivity of the skeleton is much higher than that of the PCM, causing the skeleton heat up faster than the PCM. Consequently, the temperature difference increases in the early stage of phase change. The investigation points in group L1, being closer to the heating surface, are predominantly affected by the thermal conductivity, leading to a rapid increase in temperature difference. For the other three groups of investigation points, the PCM absorbs heat from the heating surface, while the surrounding high-temperature skeleton transfers some of the heat to the PCM. This interaction results in a brief fluctuation before the temperature difference reaches its maximum peak.

As the melting process continues, the PCM stores heat as latent heat during the phase change, resulting in no further temperature increase of the PCM. However, the skeleton continues to warm up, leading to a peak in the temperature difference. The PCM takes more time to reach the phase change temperature as the distance between the group of investigation points and the heating surface increases, causing the temperature peak to occur much later. For the same Fo, the investigation points farther away from the heating surface receive less heat, leading to a smaller overall temperature difference. The temperature difference reaches the peak and then starts to decrease quickly. The reason for this is that the temperature of the liquid PCM is increasing, and natural convection facilitates the rapid heat transfer among the liquid PCMs. The investigation points of groups L2, L3, and L4 are located in the middle and rear sections of the cubic cavity, where the liquid PCM exhibits a strong natural convection. This enhances the heat flow in the liquid PCM, causing the temperature of the PCM point to increase continuously, eventually exceeding the temperature of the skeleton point. As a consequence, the temperature differences of these three groups become negative in the later stage of the phase change process. However, since L1 is closer to the heating surface, the thermal conductivity of the skeleton has a greater effect on the temperature rise. Therefore, only brief temperature fluctuations occur, and the temperature differences remain above zero.

5 Conclusions

In this paper, the dimensionless lattice Boltzmann model is developed to numerically simulate the phase change process of the PCMs containing gyroid, diamond, primitive, and I-WP skeletons. The Ra, Pr, and Ste numbers are varied to further study the flow and heat transfer characteristics of composite PCM containing the gyroid skeleton during the melting process. The LTNE effect is analyzed due to the embedding of the high thermal conductivity solid skeleton in the PCM. The numerical simulations aim to provide a research basis for the application and technical development of various TPMS structures in the field of phase change heat transfer. The conclusions drawn are as follows:

1) The addition of TPMS skeletons significantly accelerates the melting rate of the PCM, with the gyroid skeleton demonstrating the most pronounced enhancement effect on phase change heat transfer. The PCM containing the gyroid skeleton achieves complete melting 35.6% earlier than the pure PCM.

2) Different types of TPMS skeletons significantly affect the velocity distribution of the phase change process. The velocity field of the PCM containing the gyroid skeleton varies rapidly, enhancing convective heat transfer and thereby accelerating the melting process. In contrast, the I-WP skeleton has the weakest effect on enhancing the heat transfer of PCM among the four TPMS skeletons studied.

3) The variation of dimensionless parameters has a noticeable effect on the heat transfer and flow characteristics during the solid–liquid phase change process. A higher Ra number results in a stronger natural convection effect which accelerates heat transfer. For Pr = 10, the mushy zone migrates the fastest, indicating that the phase change process advances more rapidly. The composite PCMs have faster melting rates at a higher Ste number, with complete melting dimensionless time reduced by 51.79% and 59.23% when the Ste number is increased from 1 to 5 and 10, respectively.

4) The variation of temperature difference caused by the LTNE effect during the phase change process is both significant and complex. As the horizontal position of the groups of investigation points changes, the peaks of temperature difference decrease sequentially and occur later. The group of investigation points closest to the heating surface has the largest dimensionless peak of temperature difference of 0.0375. In contrast, the groups of investigation points distributed in the rear part of the cubic cavity are strongly affected by the natural convection of the melted liquid PCM, resulting in negative temperature differences.

References

[1] Jiang L L, Chen G Y, Zhao L, et al. Research on the low-temperature waste heat storage performance of a novel microencapsulated phase change material with a double-shell encapsulation of hexadecylamine. Journal of Energy Storage, 2024, 80: 110253 .

[2] Weng M M, Lin J H, Yang Y J. . MXene-based phase change materials for multi-source driven energy storage, conversion and applications. Solar Energy Materials and Solar Cells, 2024, 272: 112915

[3] Jouhara H, Żabnieńska-Góra A, Khordehgah N. . Latent thermal energy storage technologies and applications: A review. International Journal of Thermofluids, 2020, 5–6: 100039

[4] Yu N, Chen C, Zhang Y. . Discussion on optimization method of the wall in PC component solar-steam curing building based on phase change energy storage technology. Energy Procedia, 2019, 159: 16–21

[5] Khan Z, Khan Z, Ghafoor A. A review of performance enhancement of PCM based latent heat storage system within the context of materials, thermal stability and compatibility. Energy Conversion and Management, 2016, 115: 132–158

[6] Xu W H, Yang W S, Su J T. . Diatom-based biomass composites phase change materials with high thermal conductivity for battery thermal management. Journal of Energy Storage, 2024, 96: 112737

[7] Yan R H, Huang Z, Chen Y. . Phase change composite based on lignin carbon aerogel/nickel foam dual-network for multisource energy harvesting and superb EMI shielding. International Journal of Biological Macromolecules, 2024, 277: 134233

[8] Li Z C, Zhang Y, Wang X. . Thermally-induced flexible composite phase change material with enhanced thermal conductivity. Journal of Power Sources, 2024, 603: 234447

[9] Mo C M, Xie J K, Zhang G Q. . All-climate battery thermal management system integrating units-assembled phase change material module with forced air convection. Energy, 2024, 294: 130642

[10] Yao Z B, Xie J K, Fu T Q. . One-pot preparation of phase change material employing nano-scaled resorcinol-furfural frameworks. Chemical Engineering Journal, 2024, 484: 149553

[11] Bondareva N S, Sheremet M A. Influence of metal foam and nanoparticles on convective melting in a rectangular cavity under heating of bottom wall. International Journal of Thermofluids, 2023, 19: 100374

[12] Yang C, Xu Y, Cai X. . Melting behavior of the latent heat thermal energy storage unit with fins and graded metal foam. Applied Thermal Engineering, 2021, 198: 117462

[13] Guo J F, Liu Z, Du Z. . Effect of fin-metal foam structure on thermal energy storage: An experimental study. Renewable Energy, 2021, 172: 57–70

[14] Aramesh M, Shabani B. Metal foams application to enhance the thermal performance of phase change materials: A review of experimental studies to understand the mechanisms. Journal of Energy Storage, 2022, 50: 104650

[15] Ortaç B, Mutlu S, Baskan T. . Thermal conductivity and phase-change properties of boron nitride-lead oxide nanoparticle-doped polymer nanocomposites. Polymers, 2023, 15(10): 2326

[16] Jadal M, Soto J, Delaunay D. Thermal conductivity evolution of a compressed expanded natural graphite—Phase change material composite after thermal cycling. Thermal Science and Engineering Progress, 2022, 28: 101047

[17] Xie B, Cheng W L, Xu Z M. Studies on the effect of shape-stabilized PCM filled aluminum honeycomb composite material on thermal control. International Journal of Heat and Mass Transfer, 2015, 91: 135–143

[18] Nematpourkeshteli A, Iasiello M, Langella G. . Enhancing PCMs thermal conductivity: A comparison among porous metal foams, nanoparticles and finned surfaces in triplex tube heat exchangers. Applied Thermal Engineering, 2022, 212: 118623

[19] Feng W C, Ding B, Zhang Y. . How can copper foam better promote the melting process of phase change materials. International Journal of Thermal Sciences, 2023, 187: 108199

[20] Li Y J, Niu Z Y, Gao X Y. . Effect of filling height of metal foam on improving energy storage for a thermal storage tank. Applied Thermal Engineering, 2023, 229: 120584

[21] Zhuang Y J, Liu Z B, Xu W B. Effects of gradient porous metal foam on the melting performance and energy storage of composite phase change materials subjected to an internal heater: A numerical study and PIV experimental validation. International Journal of Heat and Mass Transfer, 2022, 183: 122081

[22] Marri G K, Balaji C. Experimental and numerical investigations on the effect of porosity and PPI gradients of metal foams on the thermal performance of a composite phase change material heat sink. International Journal of Heat and Mass Transfer, 2021, 164: 120454

[23] Shirbani M, Siavashi M, Bidabadi M. Phase change materials energy storage enhancement schemes and implementing the lattice Boltzmann method for simulations: A review. Energies, 2023, 16(3): 1059

[24] Rui Z L, Li J, Ding T X. . Research on thermal performance and optimization design of phase change pin fin heat sink based on Lattice Boltzmann method. Journal of Energy Storage, 2022, 55: 105659

[25] He Y L, Liu Q, Li Q. . Lattice Boltzmann methods for single-phase and solid-liquid phase-change heat transfer in porous media: A review. International Journal of Heat and Mass Transfer, 2019, 129: 160–197

[26] Chen Z Q, Gao D Y, Shi J. Experimental and numerical study on melting of phase change materials in metal foams at pore scale. International Journal of Heat and Mass Transfer, 2014, 72: 646–655

[27] Chen B M, Gao K K, Jiang H. Lattice Boltzmann study of the solid-liquid phase change mushy region. Journal of Engineering Thermophysics, 2017, 38(11): 2453–2460 (in Chinese)

[28] Gaedtke M, Abishek S, Mead-Hunter R. . Total enthalpy-based lattice Boltzmann simulations of melting in paraffin/metal foam composite phase change materials. International Journal of Heat and Mass Transfer, 2020, 155: 119870

[29] Mabrouk R, Naji H, Dhahri H. . Insight into foam pore effect on phase change process in a plane channel under forced convection using the thermal lattice Boltzmann method. Energies, 2020, 13(15): 3979

[30] Zeng C Q, Wang W. Modeling method for variable and isotropic permeability design of porous material based on TPMS lattices. Tribology International, 2022, 176: 107913

[31] Al-Ketan O, Lee D W, Rowshan R. . Functionally graded and multi-morphology sheet TPMS lattices: Design, manufacturing, and mechanical properties. Journal of the Mechanical Behavior of Biomedical Materials, 2020, 102: 103520

[32] Montazerian H, Davoodi E, Asadi-Eydivand M. . Porous scaffold internal architecture design based on minimal surfaces: A compromise between permeability and elastic properties. Materials & Design, 2017, 126: 98–114

[33] Lehder E F, Ashcroft I A, Wildman R D. . A multiscale optimisation method for bone growth scaffolds based on triply periodic minimal surfaces. Biomechanics and Modeling in Mechanobiology, 2021, 20(6): 2085–2096

[34] Liu F, Mao Z F, Zhang P. . Functionally graded porous scaffolds in multiple patterns: New design method, physical and mechanical properties. Materials & Design, 2018, 160: 849–860

[35] Hassan Ali M I, Al-Ketan O, Baobaid N. . A study on the fluid flow and heat transfer for a porous architected heat sink using the idea of CFD modelling. In: Proceedings of theASME 2019 International Mechanical Engineering Congress and Exposition. Salt Lake City: American Society of Mechanical Engineers, 2019, V008T09A060

[36] Sreedhar N, Thomas N, Al-Ketan O. . 3D printed feed spacers based on triply periodic minimal surfaces for flux enhancement and biofouling mitigation in RO and UF. Desalination, 2018, 425: 12–21

[37] Sadeghi F, Baniassadi M, Shahidi A. . TPMS metamaterial structures based on shape memory polymers: Mechanical, thermal and thermomechanical assessment. Journal of Materials Research and Technology, 2023, 23: 3726–3743

[38] Khaleghi S, Dehnavi F N, Baghani M. . On the directional elastic modulus of the TPMS structures and a novel hybridization method to control anisotropy. Materials & Design, 2021, 210: 110074

[39] Iyer J, Moore T, Nguyen D. . Heat transfer and pressure drop characteristics of heat exchangers based on triply periodic minimal and periodic nodal surfaces. Applied Thermal Engineering, 2022, 209: 118192

[40] Alteneiji M, Ali M I H, Khan K A. . Heat transfer effectiveness characteristics maps for additively manufactured TPMS compact heat exchangers. Energy Storage and Saving, 2022, 1(3): 153–161

[41] Qureshi Z A, Al Omari S A B, Elnajjar E. . Thermal characterization of 3D-printed lattices based on triply periodic minimal surfaces embedded with organic phase change material. Case Studies in Thermal Engineering, 2021, 27: 101315

[42] Qureshi Z A, Al-Omari S A B, Elnajjar E. . Using triply periodic minimal surfaces (TPMS)-based metal foams structures as skeleton for metal-foam-PCM composites for thermal energy storage and energy management applications. International Communications in Heat and Mass Transfer, 2021, 124: 105265

[43] Qureshi Z A, Elnajjar E, Al-Ketan O. . Heat transfer performance of a finned metal foam-phase change material (FMF-PCM) system incorporating triply periodic minimal surfaces (TPMS). International Journal of Heat and Mass Transfer, 2021, 170: 121001

[44] Qureshi Z A, Addin Burhan Al-Omari S, Elnajjar E. . On the effect of porosity and functional grading of 3D printable triply periodic minimal surface (TPMS) based architected lattices embedded with a phase change material. International Journal of Heat and Mass Transfer, 2022, 183: 122111

[45] Tian Y, Liu X L, Luo Q Y. . Sea urchin skeleton-inspired triply periodic foams for fast latent heat storage. International Journal of Heat and Mass Transfer, 2023, 206: 123944

[46] Fan Z H, Fu Y J, Gao R J. . Investigation on heat transfer enhancement of phase change material for battery thermal energy storage system based on composite triply periodic minimal surface. Journal of Energy Storage, 2023, 57: 106222

[47] Zhang T, Zhang K F, Liu F. . Analysis of thermal storage behavior of composite phase change materials embedded with gradient-designed TPMS thermal conductivity enhancers: A numerical and experimental study. Applied Energy, 2024, 358: 122630

[48] Feng J W, Fu J Z, Shang C. . Porous scaffold design by solid T-splines and triply periodic minimal surfaces. Computer Methods in Applied Mechanics and Engineering, 2018, 336: 333–352

[49] Yoo D J. Computer-aided porous scaffold design for tissue engineering using triply periodic minimal surfaces. International Journal of Precision Engineering and Manufacturing, 2011, 12(1): 61–71

[50] Zhang M K. Optimal design and mechanical properties of implicit surface gradient porous structure manufactured by selective laser melting. Guangzhou: South China University of Technology, 2020,

[51] Hamidi E, Ganesan P B, Sharma R K. . Computational study of heat transfer enhancement using porous foams with phase change materials: A comparative review. Renewable and Sustainable Energy Reviews, 2023, 176: 113196

[52] Voller V R, Prakash C. A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems. International Journal of Heat and Mass Transfer, 1987, 30(8): 1709–1719

[53] Hu X X, Zhang Y P. Novel insight and numerical analysis of convective heat transfer enhancement with microencapsulated phase change material slurries: Laminar flow in a circular tube with constant heat flux. International Journal of Heat and Mass Transfer, 2002, 45(15): 3163–3172

[54] Goel M, Roy S K, Sengupta S. Laminar forced convection heat transfer in microcapsulated phase change material suspensions. International Journal of Heat and Mass Transfer, 1994, 37(4): 593–604

[55] Liu X, Tong Z X, He Y L. Enthalpy-based immersed boundary-lattice Boltzmann model for solid–liquid phase change in porous media under local thermal non-equilibrium condition. International Journal of Thermal Sciences, 2022, 182: 107786

[56] Qian Y H, Zhou Y. Complete Galilean-invariant lattice BGK models for the Navier-Stokes equation. Europhysics Letters, 1998, 42(4): 359–364

[57] Guo Z L, Zheng C G, Shi B C. Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method. Chinese Physics, 2002, 11(4): 366–374

[58] Zhou W N, Li S, Feng Y H. . Evaluation and optimization on heat transfer performance of a composite phase change material embedded in porous ceramic skeleton: A lattice Boltzmann study. International Journal of Thermal Sciences, 2023, 188: 108214

[59] Kamkari B, Shokouhmand H, Bruno F. Experimental investigation of the effect of inclination angle on convection-driven melting of phase change material in a rectangular enclosure. International Journal of Heat and Mass Transfer, 2014, 72: 186–200

[60] Das M K, Reddy K S K. Conjugate natural convection heat transfer in an inclined square cavity containing a conducting block. International Journal of Heat and Mass Transfer, 2006, 49(25–26): 4987–5000

RIGHTS & PERMISSIONS

Higher Education Press 2024