Heat Transfer Enhancement Technique of PCMs and Its Lattice Boltzmann Modeling

Abstract

Phase change materials (PCMs) have several advantages for thermal energy storage due to their high energy storage density and nearly constant working temperature. Unfortunately, the low thermal conductivity of PCM impedes its efficiency of charging and discharging processes. To solve this issue, different techniques are developed to enhance the heat transfer capability of PCMs. In this chapter, the common approaches, which include the use of extended internal fins, porous matrices or metal foams, high thermal conductivity nanoparticles, and heat pipes for enhancing the heat transfer rate of PCMs, are presented in details. In addition, mathematical modeling plays a significant role in clarifying the PCM melting and solidification mechanisms and directs the experiments. As a powerful mesoscopic numerical approach, the enthalpy-based lattice Boltzmann method (LBM), which is robust to investigate the solid-liquid phase change phenomenon without iteration of source terms, is also introduced in this chapter, and its applications in latent heat thermal energy storage (LHTES) unit using different heat transfer enhancement techniques are discussed.

1. Introduction to heat transfer enhancement techniques of PCMs

The development of renewable energy such as solar energy and wind energy has attracted lots of attention during the past decades due to the gap between the increased global energy demand and the decreased amount of fossil fuel in the world. However, one of the major drawbacks for renewable energy is its territorial, time-dependent, and intermittent characteristics. Under this circumstance, the energy storage techniques play an indispensable role for achieving a continuous and reliable supply of renewable energy [1, 2]. Thermal energy storage (TES), which stores the heat in the materials and generates the electricity with heat engine cycles later, is a promising energy storage technique. In general, TES could be classified into three different categories, namely latent, sensible, and thermochemical. With the advantages of high energy storage density and nearly constant charging/discharging temperature, latent heat thermal energy storage (LHTES) using phase change materials (PCMs) is widely used in several renewable energy applications. However, a major issue of LHTES system is the low thermal conductivity of most PCMs, which seriously impedes the energy storage efficiency. To handle this challenge, several heat transfer enhancement techniques are developed and discussed by researches during the past years. The existing effective approaches to ameliorate the thermal performance of PCMs include using extended internal fins, filling porous matrix or metal foams, adding high thermal conductivity nanoparticles, and using heat pipes [3].

With the characteristics of simple fabrication, low cost construction, and large heat transfer surfaces, fins are used in a majority of PCM-based LHTES systems. There are several different configurations of fins such as longitudinal, annular, circular, plates, pins, tree shape, and other novel geometries as shown in Figure 1 [4]. By applying extended internal fins, the average thermal conductivity and heat transfer depth of LHTES system are improved, so that the melting and solidification rate of PCMs are accelerated. However, there exists a tradeoff between the increased PCM charging/discharging rate in LHTES unit and its corresponding reduced energy storage capacity with the existence of internal fins. An optimum design of fin configuration and arrangement becomes significant for achieving high energy storage efficiency of LHTES unit with PCMs. Under this circumstance, lots of researches are numerically and experimentally carried out to investigate the conjugate heat transfer between fins and PCMs, and the enhancement of PCM thermal performance with different type of fins is deeply understood and optimized during the recent years [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21].

Figure 1.

Different configuration of fins used in LHTES system with PCMs [4].

Due to the high surface area to volume ratio generated by the tortuosity of metal foams as shown in Figure 2, the PCM charging/discharging rates could be highly improved by inserting metal foams into the LHTES unit [22]. As air inevitably exists in the porous structure of metal foams, the infiltration of PCMs into metal foams is hindered, which correspondingly reduces the impregnation ratio of PCMs. As a consequence, the energy storage capacity of LHTES system using combination of PCMs and metal foams is affected. To handle this difficulty, a vacuum impregnation method is generally used to prepare the PCM/metal foam composite materials. Figure 3 shows the apparatus and procedures of impregnation treatment for PCMs with vacuum assistance, and the detailed steps could be found in Ref. [23]. As the PCM melting and solidification are actually accelerated due to the interconnected heat transfer channel inside the metal foams, the porosity and pore size of metal foams are the most significant factors, which affect the energy storage efficiency of LHTES. The copper foam with different porosities and pore sizes is displayed in Figure 4 as an example [24]. The conduction heat transfer in the LHTES system could be consolidated with the decrement of porosity and pore sizes because of the increased density of high speed heat transfer channels inside the metal foams. However, natural convective heat transfer of liquid PCMs through the metal foams is hampered due to the reduced void space caused by the decreased metal foam porosity and pore size. In addition, when the porosity of metal foams decreases, the amount of pure PCMs in the LHTES unit is reduced, which decreases the energy storage capacity. Due to the above reasons, metal foams with appropriate porosity, pore size, and filling ratio, which balance the conduction and natural convection, are essential for achieving the optimum heat transfer rate of PCMs and the most efficient energy storage of LHTES unit. Hence, the mechanisms of PCM melting and solidification processes inserted with various metal foams are studied by several researches at both macroscopic and pore scales [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. Besides, the heat transfer rate of PCMs could also be ameliorated by applying other additives such as graphite [45, 46], carbon nanotubes [47], and graphene [48].

Figure 2.

A piece of FeCrAlY foam sample manufactured with the sintering route [22].

Figure 3.

Schematic of the apparatus and procedure for the preparation of composite PCMs using a vacuum impregnation method [23].

Figure 4.

Copper foam samples with different pore sizes and porosities [24].

As the nanotechnology has achieved rapid development during the past decades, adding high thermal conductivity nanoparticles becomes a new technique to improve the low heat transfer rate of PCMs [49]. Khodadadi and Hosseinizadeh first investigated the enhancement of PCM heat transfer capability using nanoparticles [50], and their results demonstrated that nanoparticle-enhanced PCMs (NEPCMs) have a great potential in TES applications. The SEM micrographs of nanoparticle-enhanced PCM (NaNO₃-KNO₃) with different nanoparticles and mass fractions are shown in Figure 5 [51]. With the existence of nanoparticles, the thermophysical properties of PCMs such as thermal conductivity and latent heat capacity are varied. The mechanism of the effects of surface, chemical, and physical properties of nanoparticles on the thermal properties of PCMs is reviewed in Ref. [52]. In the recent years, many researches related to NEPCMs are carried out, which mainly focus on enhancing the charging/discharging speed of PCMs with nanoparticles [53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73]. However, although the effective thermal conductivity of PCMs is ameliorated by adding nanoparticles, the energy storage capacity of LHTES unit is decreased. Furthermore, the use of nanoparticles increases the viscosity of PCMs, which impedes the development of natural convective heat transfer. Under this circumstance, the total heat transfer rate of PCMs may decrease especially for the cases under high temperature with dominant convective heat transfer. Compared with melting/solidification rate of PCMs, the energy storage rate of LHTES system is the essential goal of storing heat using PCMs. Hence, more investigations, which concentrate on the energy storage rate of NEPCM, should be completed in the future research, and the technique of adding nanoparticles for enhancing the thermal performance of LHTES needs to be further compared with other heat transfer improvement technologies in order to realize the optimum energy storage efficiency.

Figure 5.

SEM photos of nanoparticle-enhanced NaNO₃-KNO₃ with 0.5 wt.% (a–d), 1.0 wt% (e–h), and 1.5 wt% (i–l) of nanoparticles. Silica (a, e, i), alumina (b, f, j), titania (c, g, k), and alumina-silica nanoparticles (d, h, l) [51].

As the most commonly used heat exchanger devices, heat pipes are widely used to amplify the charging/discharging processes of PCMs and to transfer heat from a source to the storage or from the storage to a sink with heat transfer fluid (HTF) [74]. Although increasing the heat transfer area on the PCM side using extended fins or metal foams is the most efficient and simple method to ameliorate the energy storage rate of LHTES system as previously discussed, when there exists high temperature HTF passing through the LHTES tank such as the waste heat recovery, heat pipes are indispensable for achieving the high efficiency energy storage. The transient charging/discharging processes of PCMs in a LHTES unit with heat pipes are shown in Figure 6 [75]. The configuration and arrangement of heat pipes play a significant role in the energy storage rate of PCMs. To optimize the thermal performance of heat pipe–assisted LHTES systems, lots of experimental and numerical works are carried out during the past few years [76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96].

Figure 6.

The transient charging/discharging processes of PCM in a LHTES tank with HTF passing through heat pipes [75]. (a) Charging process and (b) discharging process.

The significant research progress of PCM heat transfer enhancement using a single technique is discussed and summarized in the above paragraphs. Recently, to further improve the heat transfer capability of PCM and compare the effectiveness of different approaches (use of fins, metal foams, nanoparticles, or heat pipes), the charging/discharging processes of PCMs with hybrid heat transfer enhancement techniques are investigated. Although adding nanoparticles could ameliorate the effective thermal conductivity of PCMs, the heat transfer area on the PCM side is not improved. Based on this, the extended fins are considered to be used for enhancing the heat transfer depth of NEPCM in a LHTES unit [97, 98, 99, 100]. Darzi et al. studied the melting and solidification of PCM enhanced with radial fins and nanoparticles in cylindrical annulus, and they found that adding fins on the hot or cold tubes is the best approach to expedite the heat transfer rate [98]. Lohrasbi et al. optimized the copper nanoparticle-PCM solidification process in a fin-assisted LHTES system [99]. The results indicate that immersing fin in LHTES unit increases the solidification rate more significantly than dispersing nanoparticles. Parsazadeh and Duan investigated the effects of fins and nanoparticle in a shell and tube LHTES unit [100]. They found that adding Al₂O₃ nanoparticles even decreases the overall heat transfer rate because the thermal conductivity enhancement with nanoparticles could not compensate for the natural convection reduction. Similarly, porous matrices are also inserted into the LHTES unit to improve the thermal performance of NEPCM [101, 102]. Hossain et al. studied the melting process of NEPCM inside the porous medium [101], and it is observed that the movement of PCM melting front is more significant under the influence of porous medium than that of nanoparticles. Tasnim et al. investigated the convection effect on the melting process of NEPCM filled in porous enclosure [102]. The results showed that both the conduction and convection heat transfer are degraded by the presence of nanoparticles. From these researches, it could be found that using extended fins or porous matrices is more effective than adding nanoparticles for enhancing the charging/discharging rate of LHTES system. Besides, other hybrid heat transfer enhancement techniques for enhancing the energy storage rate of LHTES with PCMs such as combination of fins and metal foams [103] or using combined three techniques [104, 105, 106, 107] are also recently studied and analyzed.

In this chapter, the mathematical models for PCM charging and discharging processes with different heat transfer enhancement techniques are shown. In addition, the lattice Boltzmann method (LBM) for solid-liquid phase change phenomenon is reviewed with some classical analytical and numerical validation cases, and the implementation of graphic processor unit (GPU) computing is also presented. Furthermore, the applications of LBM modeling for LHTES system with various heat transfer improvement approaches are discussed.

2. Mathematical models

2.1 Governing equations of fluid flow and solid-liquid phase change

To simulate the charging and discharging processes of PCMs, the following assumptions are usually made to simplify the mathematical models: (1) the fluid flow of liquid PCMs is Newtonian, laminar, and incompressible with negligible viscous dissipation and (2) the thermophysical properties of PCMs, nanoparticles, fins, and metal foams are constant, except the PCM density $\rho$, which is a linear function of temperature $T$ using the Boussinesq approximation. Based on the above assumptions, the flow of liquid PCMs is governed by the continuity equation and the momentum equation expressed in the Cartesian coordinate as:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$$E1

$$\rho \left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}\right)=-\frac{\partial p}{\partial x}+\mu \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)$$E2

$$\rho \left(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}\right)=-\frac{\partial p}{\partial y}+\mu \left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}+\frac{{\partial }^{2}v}{\partial z^2}\right)+\left(\mathit{\rho \beta }\right)g\left(T-T_r\right)$$E3

$$\rho \left(\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}\right)=-\frac{\partial p}{\partial z}+\mu \left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}+\frac{{\partial }^{2}w}{\partial z^2}\right)$$E4

where $\rho$ is the density; $p$ is the pressure; $\mu$ is the dynamic viscosity; $t$ is the time; $\beta$ is the thermal expansion coefficient of PCMs; $g$ is the magnitude of gravitational acceleration; $T$ is the temperature; $T_r$ is the reference temperature; $x$ is the horizontal coordinate; $y$ is the vertical coordinate; $z$ is the coordinate, which is orthogonal with $x$ and $y$ coordinates; and $u$, $v$, and $w$ are the fluid velocities in the $x$, $y$, and $z$ directions, respectively.

The solid-liquid phase change process of PCMs is governed by the enthalpy equation as:

$$\frac{\partial \left(\mathit{\rho H}\right)}{\partial t}+\rho c_p\left(u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}\right)=k\left(\frac{{\partial }^{2}T}{{\partial x}^2}+\frac{{\partial }^{2}T}{{\partial y}^2}+\frac{{\partial }^{2}T}{{\partial z}^2}\right)$$E5

where $c_p$ and $k$ are the specific heat and thermal conductivity of PCMs. The enthalpy $H$ of PCMs is defined as:

$$H=c_p\left(T-T_r\right)+f_{l}L$$E6

where $f_l$ is the PCM liquid fraction, and $L$ is the latent heat of PCMs. By calculating the enthalpy $H$ of PCMs, the liquid fraction $f_l$ and temperature $T$ could be updated by the following equations:

$$f_l=\left\{\begin{array}{cc}0,& H\le H_s\\ \frac{H-H_s}{H_l-H_s},& H_s<H<H_l\\ 1& H\ge H_l\end{array}\right.$$E7

$$T=\left\{\begin{array}{c}{T}_m-\frac{H_s-H}{c_p},H\le H_s\\ T_m,H_s<H<H_l\\ T_m+\frac{H-H_l}{c_p}H\ge H_l\end{array}\right.$$E8

where $H_s$ is the enthalpy of solid state PCMs, $H_l$ is the enthalpy of liquid state PCMs, and $T_m$ is the melting/solidification temperature of PCMs.

2.2 Nanofluid models

For simulating the fluid flow and heat transfer of nanoparticle-enhanced PCMs (NEPCMs), the nanoparticle is assumed to be spherical in shape, so that the Brinkman model and the Maxwell model for nanofluid are valid. In addition, the NEPCMs are considered as continuous media with the thermal dispersion being neglected. Based on this, the effective viscosity of NEPCMs is computed using the Brinkman model as [108]:

$${\mu }_{\mathit{nf}}=\frac{{\mu }_{\mathit{PCM}}}{{\left(1-\Phi \right)}^{2.5}}$$E9

where $\Phi$ is the volume fraction of nanoparticles, ${\mu }_{\mathit{PCM}}$ is the dynamic viscosity of pure PCM, and ${\mu }_{\mathit{nf}}$ is the dynamic viscosity of NEPCMs. The thermal conductivity of NEPCMs is calculated according to the Maxwell model as [109]:

$$\frac{k_{\mathit{nf}}}{k_{\mathit{PCM}}}=\frac{k_p+2k_{\mathit{PCM}}-2\left(k_{\mathit{PCM}}-k_p\right)\Phi }{k_p+2k_{\mathit{PCM}}+\left(k_{\mathit{PCM}}-k_p\right)\Phi }$$E10

where $k_{\mathit{PCM}}$, $k_p$, and $k_{\mathit{nf}}$ are the thermal conductivities of pure PCMs, nanoparticles, and NEPCMs, respectively. Furthermore, the density of nanofluid ${\rho }_{\mathit{nf}}$ is computed using interpolation as:

$${\rho }_{\mathit{nf}}=\left(1-\Phi \right){\rho }_{\mathit{PCM}}+\Phi {\rho }_p$$E11

where ${\rho }_{\mathit{PCM}}$ and ${\rho }_p$ are the densities of pure PCM and nanoparticles. The heat capacitance of NEPCMs ${\left(\rho c_p\right)}_{\mathit{nf}}$ is defined as:

$${\left(\rho c_p\right)}_{\mathit{nf}}=\left(1-\Phi \right){\left(\rho c_p\right)}_{\mathit{PCM}}+\Phi {\left(\rho c_p\right)}_p$$E12

where ${\left(\rho c_p\right)}_{\mathit{PCM}}$ is the heat capacitance of pure PCM, and ${\left(\rho c_p\right)}_p$ is the heat capacitance of nanoparticles. The thermal expansion volume of NEPCMs ${\left(\mathit{\rho \beta }\right)}_{\mathit{nf}}$ is given as:

$${\left(\mathit{\rho \beta }\right)}_{\mathit{nf}}=\left(1-\Phi \right){\left(\mathit{\rho \beta }\right)}_{\mathit{PCM}}+\Phi {\left(\mathit{\rho \beta }\right)}_p$$E13

where ${\left(\mathit{\rho \beta }\right)}_{\mathit{PCM}}$ and ${\left(\mathit{\rho \beta }\right)}_p$ are the thermal expansion volume of pure PCM and nanoparticles, respectively. The latent heat of NEPCMs ${\left(\mathit{\rho L}\right)}_{\mathit{nf}}$ is computed as:

$${\left(\mathit{\rho L}\right)}_{\mathit{nf}}=\left(1-\Phi \right){\left(\mathit{\rho L}\right)}_{\mathit{PCM}}$$E14

where ${\left(\mathit{\rho L}\right)}_{\mathit{PCM}}$ is the latent heat of pure PCM. Then, the corresponding enthalpy of NEPCMs $H_{\mathit{nf}}$ is given as:

$$H_{\mathit{nf}}={c_p}_{\mathit{nf}}\left(T-T_r\right)+f_l{L}_{\mathit{nf}}$$E15

2.3 Conjugate heat transfer between PCMs and fins or metal foams

When the extended fins or metal foams are used as the heat transfer enhancement techniques for LHTES unit, the conjugate heat transfer occurs between the PCMs and the fins or metal foams. The heat transfer inside the fins or metal foams is governed by the conduction equation as:

$${\rho }_s{c_p}_s\frac{\partial T}{\partial t}=k_s\left(\frac{{\partial }^{2}T}{{\partial x}^2}+\frac{{\partial }^{2}T}{{\partial y}^2}+\frac{{\partial }^{2}T}{{\partial z}^2}\right)$$E16

${\rho }_s$, ${c_p}_s$, and $k_s$ are the density, specific heat, and thermal conductivity of fins or metal foams, respectively. On the interface of PCMs and fins or metal foams, the coupled Dirichlet-Neumann boundary conditions for conjugate heat transfer should be satisfied:

$$-k_{\mathit{PCM}}\frac{\partial T_{\mathit{PCM}}}{\partial n}=-k_s\frac{\partial T_s}{\partial n}$$E17

$$T_{\mathit{PCM}}=T_s$$E18

$T_{\mathit{PCM}}$ and $T_s$ are the temperature of PCMs and fins or metal foams, $k_{\mathit{PCM}}$ is the thermal conductivity of PCMs, and $n$ is the normal unit of the interface. The above boundary conditions could be satisfied in the lattice Boltzmann method automatically by using the harmonic mean value of thermophysical properties of PCMs and fins or metal foams as discussed in the following section. To investigate the melting and solidification processes of PCMs filled with metal foams at pore scale, the morphology of the real metal foam structures could be reconstructed using the quartet structure generation set (QSGS) as shown in Figure 7 [110, 111], where ${\epsilon }_{\mathit{ave}}$ and $d_p$ are the porosity and the pore size of metal foams. Besides, the metal foams could also be reconstructed with scanning electron microscopy (SEM).

Figure 7.

Reconstructed PCM filled with metal foams using QSGS [111].

3. Lattice Boltzmann method for solid-liquid phase change

3.1 The development history of enthalpy-based LBM for solid-liquid phase change

As a mesoscopic numerical approach developed during the past more than two decades, the lattice Boltzmann method (LBM) has become a powerful tool for simulating complex heat transfer and fluid flow problems such as single-phase flows, multiphase flows, turbulence flows, flows in porous media, heat transfer, phase change, and transport in microfluidics [112, 113, 114, 115]. The mesoscopic nature of LBM makes it appropriate for tackling the evolvement of solid-liquid interface during the phase change process. The existing LBM for solid-liquid phase change problems could be generally categorized as follows: (1) the phase-field based method, (2) the enthalpy-based method, and (3) the immersed boundary method. For the phase-field method, an auxiliary parameter, which varies smoothly across the diffusive phase interface, is used to track the solid-liquid interface implicitly [116, 117]. Unfortunately, the extremely finer grids are indispensable in the interfacial region, which increases the computational effort. In addition, Huang and Wu developed the immersed boundary-thermal lattice Boltzmann method for modeling the solid-liquid phase change phenomenon [118]. The melting interface is explicitly tracked by the Lagrangian grids, and the temperature and fluid flow are solved on the Eulerian nodes. However, an interpolation between the Lagrangian nodes and the Eulerian nodes should be carried out using the Dirac delta function, which reduces the computation accuracy. Compared with other methods, the enthalpy-based method becomes attractive for solid-liquid phase change due to its high efficiency and robustness. Jiang et al. first used LBM with enthalpy formulation to investigate the phase change problems [119]. However, the convective heat transfer is not considered in this work. Huber et al. developed a LBM model using double distribution functions to simulate the solid-liquid phase change process with convective heat transfer [120], and the numerical results are analyzed and compared with scaling laws and previous numerical work. Eshraghi and Felicelli developed an implicit LBM scheme to investigate the conduction with phase change [121]. Different from the previous enthalpy-based LBM, the iteration of nonlinear latent heat source term in the energy equation is avoided by solving a linear system of equations. Feng et al. further extended this implicit model with a consideration of natural convection to investigate the melting process of nanoparticle-enhanced PCMs [122]. To improve the computational efficiency, Huang et al. modified the equilibrium distribution function of temperature for enthalpy, so that the iteration of heat source term or solving a linear system of equations is not indispensable [123]. Unfortunately, the stability range of relaxation time for this model is narrow, which limits its applications for real scientific and engineering problems. To handle this drawback, Huang and Wu improved their work by using a multiple-relaxation-time (MRT) collision scheme instead of the single-relaxation-time (SRT) one, so that the numerical stability is highly ameliorated [124]. In addition, the thermal conductivity and the specific heat capacity are decoupled from the relaxation time and equilibrium distribution function, which make this model appropriate for satisfying the Dirichlet-Neumann boundary conditions for conjugate heat transfer. As the model developed by Huang et al. is limited in two-dimensional cases, Li et al. recently developed the SRT and MRT models for axisymmetric and three-dimensional solid-liquid phase change problems [125, 126]. Besides, to speed up the computation, Su and Davidson worked out a mesoscopic scale timestep adjustable non-dimensional LBM [127], and the time steps can be adjusted independent of mesh size by changing the transient mesoscopic Mach number. Furthermore, although the key advantage of LBM is to carry out the pore-scale numerical modeling of heat transfer in porous media, a few lattice Boltzmann models for solid-liquid phase change in porous media at the representative elementary volume (REV) scale are also developed [128, 129, 130]. The classical enthalpy-based MRT lattice Boltzmann model developed by Huang and Wu is reviewed in this section because of its simplicity and wide application by researchers for latent heat thermal energy storage problems.

3.2 Multiple-relaxation-time (MRT) method

3.2.1 MRT LBM for fluid flow

The general two dimensional nine-velocity (D2Q9) MRT LBM model is presented in this part for simulating the fluid flow. In the D2Q9 model, lattice velocities ${\mathit{e}}_i$ are given by:

$${\mathit{e}}_i=\left\{\begin{array}{c}{\mathit{e}}_0=\left(0,0\right)\hfill \\ {\mathit{e}}_i=c\left(cos\left[\left(i-1\right)\pi /2\right],sin\left[\left(i-1\right)\pi /2\right]\right),i=1,2,3,4\hfill \\ {\mathit{e}}_i=\sqrt{2}c\left(cos\left[\left(2i-9\right)\pi /4\right],sin\left[\left(2i-9\right)\pi /4\right]\right),i=5,6,7,8\hfill \end{array}\right.$$E19

where $c$ is the lattice speed, and the collision step carried out in the momentum space is given as:

$${\mathit{m}}_{\mathit{f}}\left(\mathit{x},t+{\delta }_t\right)={\mathit{m}}_{\mathit{f}}\left(\mathit{x},t\right)-\mathit{S}\left[{\mathit{m}}_{\mathit{f}}\left(\mathit{x},t\right)-{\mathit{m}}_{\mathit{f}}^{\mathit{eq}}\left(\mathit{x},t\right)\right]+{\delta }_t\left(\mathit{I}-\frac{\mathit{S}}{2}\right){\mathit{F}}_{\mathit{m}}\left(\mathit{x},t\right)$$E20

where $\mathit{x}$ is the location vector, $t$ is the time, ${\delta }_t$ is the time step, $\mathit{I}$ is the unit matrix, and ${\mathit{m}}_{\mathit{f}}$ is the distribution function in momentum space defined as:

$${\mathit{m}}_{\mathit{f}}\left(\mathit{x},t\right)={\left[m_{f0}\left(\mathit{x},t\right),m_{f1}\left(\mathit{x},t\right),\dots ,m_{f8}\left(\mathit{x},t\right)\right]}^T$$E21

The equilibrium distribution function in momentum space ${\mathit{m}}_{\mathit{f}}^{\mathit{eq}}$ is [131, 132, 133]:

$${\mathit{m}}_{\mathit{f}}^{\mathit{eq}}\mathbf{=}{\left(\rho ,-2\rho +3\rho \frac{{\left|\mathit{u}\right|}^2}{c^2},\rho -3\rho \frac{{\left|\mathit{u}\right|}^2}{c^2},\rho \frac{u_x}{c},-\rho \frac{u_x}{c},\rho \frac{u_y}{c},-\rho \frac{u_y}{c},\rho \frac{u_x^2-u_y^2}{c^2},\rho \frac{u_x{u}_y}{c^2}\right)}^T$$E22

where $\rho$ is the density, $\mathit{u}$ is the fluid velocity vector, $u_x$ is the horizontal velocity, and $u_y$ is the vertical velocity. The diagonal relaxation matrix $\mathit{S}$ is defined as:

$$\mathit{S}\mathbf{=}\mathit{diag}\left(s_0,s_e,s_{\epsilon },s_j,s_q,s_j,s_q,s_p,s_p\right)$$E23

where $s_0$, $s_e$, $s_{\epsilon }$, $s_j$, $s_q$, and $s_p$ are the parameters related to relaxation time, which could be chosen as described in Ref. [134]. The discrete force term in the momentum space ${\mathit{F}}_{\mathit{m}}\left(\mathit{x},t\right)$ is given by [135, 136]:

$${\mathit{F}}_{\mathit{m}}\left(\mathit{x},t\right)={\left(0,6\frac{\mathit{F}\bullet \mathit{u}}{c^2},-6\frac{\mathit{F}\bullet \mathit{u}}{c^2},\frac{F_x}{c},-\frac{F_x}{c},\frac{F_y}{c},-\frac{F_y}{c},2\frac{F_x{u}_x-F_y{u}_y}{c^2},\frac{F_x{u}_y+F_y{u}_x}{c^2}\right)}^T$$E24

where $\mathit{F}$ is the body force of fluid flow, $F_x$ is the body force in the horizontal direction, and $F_y$ is the body force in the vertical direction. After the collision process, the post collision distribution function in the velocity space $f_i\left(\mathit{x},t+{\delta }_t\right)$ is calculated through inverse transformation:

$$f_i\left(\mathit{x},t+{\delta }_t\right)={\mathit{M}}^{\mathbf{-}\mathbf{1}}{\mathit{m}}_{\mathit{f}}\left(\mathit{x},t+{\delta }_t\right)$$E25

where the dimensionless orthogonal transformation matrix $\mathit{M}$ is chosen as [137]:

$$\mathit{M}=\left\{\begin{array}{ccccccccc}1& 1& 1& 1& 1& 1& 1& 1& 1\\ -4& -1& -1& -1& -1& 2& 2& 2& 2\\ 4& -2& -2& -2& -2& 1& 1& 1& 1\\ 0& 1& 0& -1& 0& 1& -1& -1& 1\\ 0& -2& 0& 2& 0& 1& -1& -1& 1\\ 0& 0& 1& 0& -1& 1& 1& -1& -1\\ 0& 0& -2& 0& 2& 1& 1& -1& -1\\ 0& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& -1& 1& -1\end{array}\right\}$$E26

Then, the streaming step is carried out as:

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

The nonslip velocity condition on the diffusive interface and in the solid phase is tackled by recalculating the density distribution function $f_i$ through a linear interpolation as [138]:

$$f_i=f_l{f}_i+\left(1-f_l\right)f_i^{\mathit{eq}}\left(\rho ,{\mathit{u}}_{\mathit{s}}\right)$$E28

$f_l$ is the liquid fraction of PCM, $f_i^{\mathit{eq}}$ is the equilibrium distribution function, which could be calculated by inverse transformation as $f_i^{\mathit{eq}}={\mathit{M}}^{\mathit{-}\mathbf{1}}{\mathit{m}}_{\mathit{f}}^{\mathit{eq}}$. For the solid phase, there is ${\mathit{u}}_{\mathit{s}}\mathit{=}\mathbf{0}$. Hence, the macroscopic variables, density $\rho$ and velocity $\mathit{u}$, are defined as:

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

As mentioned in Ref. [138], the density $\rho$ in the term $f_i^{\mathit{eq}}\left(\rho ,{\mathit{u}}_{\mathit{s}}\right)$ in Eq. (28) should be first calculated by Eq. (29) in order to ensure the mass conservation. Then, for the liquid phase of $f_l=1$, the above lattice Boltzmann model recovers to the standard scheme for incompressible flow. On the other hand, for the solid phase of $f_l=0$, the above model could satisfy that $f_i=f_i^{\mathit{eq}}\left(\rho ,{\mathit{u}}_{\mathit{s}}\right)$ indicating that the nonslip velocity $\mathit{u}\mathit{=}{\mathit{u}}_{\mathit{s}}$ is ensured.

3.2.2 MRT LBM for solid-liquid phase change

The MRT lattice Boltzmann equation (LBE) for the total enthalpy $H$ distribution function ${\mathit{m}}_{\mathit{g}}\left(\mathit{x},t\right)$ is expressed as [124]:

$${\mathit{m}}_{\mathit{g}}\left(\mathit{x},t+{\delta }_t\right)={\mathit{m}}_{\mathit{g}}\left(\mathit{x},t\right)-\mathit{S}\left[{\mathit{m}}_{\mathit{g}}\left(\mathit{x},t\right)-{\mathit{m}}_{\mathit{g}}^{\mathit{eq}}\left(\mathit{x},t\right)\right]$$E30

where ${\mathit{m}}_{\mathit{g}}$ is the distribution in momentum space given as:

$${\mathit{m}}_{\mathit{g}}\left(\mathit{x},t\right)={\left[m_{g0}\left(\mathit{x},t\right),m_{g1}\left(\mathit{x},t\right),\dots ,m_{g8}\left(\mathit{x},t\right)\right]}^T$$E31

The equilibrium moment ${\mathit{m}}_{\mathit{g}}^{\mathit{eq}}$ is given by:

$$\begin{array}{c}{m}^{eq}=(H\text{,}-4H+2c_{p,ref}T+3c_{p}T\frac{u^2}{c^2}\text{,}4H-3c_{p,ref}T-3c_{p}T\frac{u^2}{c^2}\text{,}{c}_{p}T\frac{u_x}{c}\text{,}\\ {-c_{p}T\frac{u_x}{c}\text{,}{c}_{p}T\frac{u_y}{c}\text{,}-c_{p}T\frac{u_y}{c}\text{,}{c}_{p}T\frac{u_x^2-u_y^2}{c^2}\text{,}{c}_{p}T\frac{u_x{u}_y}{c^2})}^T\end{array}$$E32

where $T$ is the temperature, and $c_p$ is the specific heat calculated by interpolation as:

$$c_p=\left(1-f_l\right)c_{p,s}+f_l{c}_{p,l}$$E33

$c_{p,s}$ is the specific heat of PCM at solid state, and $c_{p,l}$ is the specific heat of PCM at liquid state. To achieve good numerical stability, the reference specific heat $c_{p,\mathit{ref}}$ is defined by the harmonic mean of $c_{p,s}$ and $c_{p,l}$ as:

$$c_{p,\mathit{ref}}=\frac{2c_{p,s}{c}_{p,l}}{c_{p,s}+c_{p,l}}$$E34

The parameters in the diagonal relaxation matrix $\mathit{S}$ satisfy $s_0=1$, $s_e=s_p$, $s_j=\frac{1}{\tau }$, and $0<s_{e,\epsilon ,q}<2$, where the relaxation time $\tau$ is given as:

$$\frac{k}{\rho c_{p,\mathit{ref}}}=c_s^2\left(\tau -0.5\right){\delta }_t$$E35

where $k$ is the thermal conductivity, and $c_s=\frac{c}{\sqrt{3}}$ is the sound speed. To reduce the numerical diffusion, a “magic” relationship is found by Huang and Wu [124] as:

$$\left(\frac{1}{s_e}-\frac{1}{2}\right)\left(\frac{1}{s_j}-\frac{1}{2}\right)=\frac{1}{4}$$E36

Similar to the computation of fluid flow, the post-collision distribution function in the velocity space $g_i$ could be calculated by inverse transformation as:

$$g_i\left(\mathit{x},t+{\delta }_t\right)={\mathit{M}}^{\mathbf{-}\mathbf{1}}{\mathit{m}}_{\mathit{g}}\left(\mathit{x},t+{\delta }_t\right)$$E37

The streaming process is completed as:

$$g_i\left(\mathit{x}+{\mathit{e}}_i{\delta }_t,t+{\delta }_t\right)=g_i\left(\mathit{x},t+{\delta }_t\right)$$E38

Then, the enthalpy $H$ is calculated as:

$$H=\sum _{i=0}^8{g}_i$$E39

Furthermore, it should be pointed out that the nonequilibrium extrapolation method developed by Guo et al. could be applied for the boundary conditions of fluid flow and enthalpy on the surfaces of LHTES unit [139, 140].

3.3 Classical examples for code validation

3.3.1 One-dimensional transient conjugate heat transfer

To validate the capability of MRT LBM for tackling the differences in thermophysical properties, the one-dimensional transient conjugate heat transfer in two regions without phase change is used to compare the numerical results with analytical solutions. Initially, T is equal to 1 in the region A at $x>0$, and T is equal to 0 in the region B at $x<0$. The analytical solution for this problem is given as [141]:

$$T^A\left(x,t\right)=\frac{1}{1+\sqrt{{\left(\rho C_p\right)}^B{k}^B/{\left(\rho C_p\right)}^A{k}^A}}[1+\sqrt{{\left(\rho C_p\right)}^B{k}^B/{\left(\rho C_p\right)}^A{k}^A}\mathit{erf}\left(\frac{x}{2\sqrt{k^{A}t/{\left(\rho C_p\right)}^A}}\right)]$$E40

$$T^B\left(x,t\right)=\frac{1}{1+\sqrt{{\left(\rho C_p\right)}^B{k}^B/{\left(\rho C_p\right)}^A{k}^A}}\mathit{erfc}\left(-\frac{x}{2\sqrt{k^{B}t/{\left(\rho C_p\right)}^B}}\right)$$E41

When the aluminum is used for region A while the liquid water is chosen as the material in region B, the comparison between LBM results and analytical solutions is presented in Figure 8. A good agreement is observed, and the maximum L2 error with 200 × 200 lattice grids is only $6.7983\times {10}^{-5}$ [142].

Figure 8.

One-dimensional transient conjugate heat transfer [142].

3.3.2 One-dimensional melting by conduction

In order to verify the MRT LBM for solid-liquid phase change phenomenon, the one-dimensional melting by conduction at a constant phase change temperature $T_m$ is simulated. Initially, the substance is uniformly solid at a temperature $T_0$ ($T_0<T_m$). The melting process begins at time $t=0$ when the temperature of the left wall is at a high temperature of $T_h$ ($T_h$ > $T_m$). Then, the analytical solution for the temperature in this problem is given as [143]:

$$T\left(x,t\right)=T_h-\frac{T_h-T_m}{erf\left(\lambda \right)}\mathit{erf}\left(\frac{x}{2\sqrt{k_{l}t/\left(\rho c_{p,l}\right)}}\right),0\le x\le X_i\left(t\right),\text{liquid}$$E42

$$T\left(x,t\right)=T_0+\frac{T_m-T_0}{\mathit{erfc}\left(\lambda /\sqrt{R_{\mathit{ac}}}\right)}\mathit{erfc}\left(\frac{x}{2\sqrt{k_{s}t/\left(\rho c_{p,s}\right)}}\right),x>X_i\left(t\right),\text{solid}$$E43

where $k_l$ and $k_s$ are the thermal conductivities of liquid PCM and solid PCM, respectively. The parameter $R_{\mathit{ac}}=\frac{\left[k_{s}t/\left(\rho c_{p,s}\right)\right]}{\left[k_{l}t/\left(\rho c_{p,l}\right)\right]}$ is the ratio of thermal diffusivity between the solid and liquid phases of PCM, and the location of phase interface $X_i$ is calculated as:

$$X_i=2\lambda \sqrt{k_{l}t/\left(\rho c_{p,l}\right)}$$E44

The parameter $\lambda$ is the root of the transcendental equation:

$$\frac{{\mathit{Ste}}_l}{\mathit{exp}\left({\lambda }^2\right)\mathit{erf}\left(\lambda \right)}-\frac{{\mathit{Ste}}_s\sqrt{R_{\mathit{ac}}}}{\mathit{exp}\left({\lambda }^2/R_{\mathit{ac}}\right)\mathrm{erfc}\left(\lambda /\sqrt{R_{\mathit{ac}}}\right)}=\lambda \sqrt{\pi }$$E45

where the Stefan numbers, ${\mathit{Ste}}_l$ and ${\mathit{Ste}}_s$, are defined as:

$${\mathit{Ste}}_l=\frac{c_{p,l}\left(T_h-T_m\right)}{L},{\mathit{Ste}}_s=\frac{c_{p,s}\left(T_m-T_0\right)}{L}$$E46

The comparison of temperature $T$ between the analytical solutions and the SRT or MRT LBM with different relaxation times is presented in Figure 9. It could be observed that the numerical diffusion exists for the small $\left({\tau }_l\text{close to}0.5\right)$ and large $\left({\tau }_l>2\right)$ relaxation times when the SRT model is applied. However, the numerical diffusion is highly reduced once the MRT model with a “magic” parameter relation shown in Eq. (36) is used.

Figure 9.

(a) SRT LBM (b) MRT LBM. Comparison of temperature $T$ between analytical solutions and SRT or MRT LBM results with different relaxation times [124].

3.3.3 Two-dimensional melting by convection

The natural convection with melting in a square cavity heated from the side wall is usually used to validate the code for solid-liquid phase change. First, the solid-liquid phase change with convection in a cavity at the Rayleigh number $\mathit{Ra}=25000$, the Prandtl number $\mathit{Pr}=0.02$, and the Stefan number $\mathit{Ste}=0.01$ is compared with the results by Mencinger [144]. The average Nusselt number ${\mathit{Nu}}_{\mathit{ave}}$ at the hot wall in terms of the Fourier numbers $\mathit{Fo}$ is plotted in Figure 10(a), and the melting interface positions at $\mathit{Fo}=10$ and $\mathit{Fo}=20$ are shown in Figure 10(b). In addition, the average melting fraction $f_l$ is also presented in Figure 10(c). It is obvious that the current LBM results are consistent with the work of Mencinger [142]. Besides, for the cases of $\mathit{Pr}>1$, the MRT LBM for natural convection with melting in cavity could be calibrated by the scaling laws and correlations derived by Jany and Bejan [145]. The average Nusselt number ${\mathit{Nu}}_{\mathit{ave}}$ at hot wall is given as:

$${\mathit{Nu}}_{\mathit{ave}}={\left(2\mathit{\text{FoSte}}\right)}^{-1/2}+\left[0.33{\mathit{Ra}}^{1/4}-{\left(2\mathit{\text{FoSte}}\right)}^{-1/2}\right]{\left\{1+{\left[0.0175{\mathit{Ra}}^{3/4}{\left(\mathit{\text{FoSte}}\right)}^{3/2}\right]}^{-2}\right\}}^{-1/2}$$E47

Figure 10.

Natural convection with melting in a cavity [142].

As shown in Figure 10(d), the average Nusselt number ${\mathit{Nu}}_{\mathit{ave}}$ at $\mathit{Pr}=6.1989$ and $\mathit{Ste}=0.1$ agrees well with the results of scaling law correlations at different Rayleigh number $\mathit{Ra}$.

3.4 GPU acceleration

The characteristic of highly parallel nature is a significant advantage of LBM over other traditional macroscopic numerical methods. In the CUDA programming platform, the CPU and GPU work as the host and the devices, respectively, as shown in Figure 11 [146]. It means that the parallel tasks are executed on GPU, while the CPU is responsible for the initial conditions and all the sequential commands. First, the initial conditions are set up in the host memory, and then the data are moved to the memory of GPU. The threads are grouped into blocks, which are the component of grids as displayed in Figure 11. For instance, in the real simulation, if the lattice grid number is (M_x, M_y, M_z) and the block size is (N_x, N_y, N_z), the corresponding grid size is $\left(M_x/N_x,M_y/N_y,M_z/N_z\right)$. Furthermore, a kernel is a function, which is executed on the concurrent threads on GPU. The collision step, streaming step, tackling of boundary conditions, and computation of macroscopic variables are completed in different kernels. It is important to note that all the kernels should be synchronized between CPU and GPU. Finally, the data on the GPU should be copied to the host memory of CPU for printing at any specific time when the results are needed. During the recent years, LBM is demonstrated to be appropriate for GPU computing [147, 148, 149, 150, 151], and it has been applied to solve several different physical problems [152, 153, 154, 155].

Figure 11.

Schematic diagram of CUDA platform [146].

4. Applications of LBM modeling in latent heat thermal energy storage

4.1 Enhancement of PCM performance with fins

As presented in the previous sections, the extended fins could be used to increase the heat transfer depth in LHTES system and accelerate the energy storage efficiency. In the recent years, LBM is applied by several researches to study the solid-liquid phase change of PCMs with extended fins [142, 156, 157]. Jourabian et al. studied the melting process of PCM in a cavity with a horizontal fin heated from the sidewall by enthalpy-based D2Q5 LBM [156]. The results indicated that adding a fin enhances the melting rate for all positions and different lengths compared with the LHTES cavity without fin. They also found that although varying the position of the fin from the bottom to the middle has a negligible effect on melting rate, the melting time is increased once the fin is mounted on the top of LHTES cavity. Talati and Taghilou used an implicit LBM to study the PCM solidification in a rectangular finned container [157]. It was found that the optimum aspect ratio of container for solidification equals to $0.5$, and changing the fin material from aluminum to copper has no significant influence on the solidification rate. Ren and Chan applied enthalpy-based MRT LBM to investigate the PCM melting performance in an enclosure with internal fins and finned thick walls with GPU acceleration [142]. The transient PCM melting process with different number of fins at the Fourier number $\mathit{Fo}=0.15$ is shown in Figure 12 in terms of temperature contours, liquid fraction, and streamlines. It could be found that the PCM melting rate is obviously enhanced by adding more internal fins in the cavity. However, the energy storage capacity of LHTES system is reduced when the number of internal fins increases, so that the appropriate fin configuration and number should be designed for engineering applications. They found that using a less number of longer fins is more effective than applying shorter fins for enhancing the thermal performance of PCMs. Besides, compared with the LHTES cavity with horizontal fins heated from side walls, the LHTES enclosure using vertical fins heated from the bottom surface has a better charging rate. From the above researches, it should be concluded that the enthalpy-based LBM is successful for simulating the conjugate heat transfer with solid-liquid phase change for melting and solidification of PCMs accelerated with fins.

Figure 12.

Transient PCM melting process with different number of fins [142].

4.2 Nanoparticle-enhanced PCM

The nanoparticles are commonly used to ameliorate the low thermal conductivity of most PCMs. With some appropriate assumptions as presented in Section 2, the thermal performance of PCMs with nanoparticles could be modeled using enthalpy-based LBM. Jourabian et al. studied the convective melting process of Cu-nanoparticle enhanced water-ice in a cylindrical-horizontal annulus [158, 159]. It was demonstrated that the melting rate of NEPCM is accelerated due to the improvement of thermal conductivity and the reduced latent heat. When the heated cylinder located in the bottom section of the annulus, the increment effect on NEPCM melting rate by adding more nanoparticles decreases because of the augmentation of NEPCM viscosity, which weakens the convective heat transfer. Feng et al. investigated the melting of water (ice)-copper nanoparticle NEPCM in a bottom-heated rectangular cavity by treating the latent heat source term with an implicit scheme [122]. They also found that the heat transfer rate of NEPCM increases with respect to the increment of nanoparticle volume fractions. However, the energy storage rate is the most significant parameter for a LHTES unit. Although adding nanoparticles into PCMs could increase their thermal conductivity, it increases the viscosity of PCM, which weakens the convective heat transfer. Due to this reason, the energy storage rate of LHTES system may even decrease with the increasing nanoparticle volume fractions, especially for the case with large temperature gradient. On the other hand, the increment of nanoparticle volume fraction decreases the energy storage capacity of LHTES unit, so that the energy storage rate could also be affected. Under this circumstance, more future research attentions should be paid on the influences of nanoparticles on the energy storage rate of NEPCM. In addition, the solid-liquid phase change model of NEPCM using LBM could be extended to study the charging and discharging of NEPCM under hybrid heat transfer enhancement techniques such as nanoparticle-fin or nanoparticle-metal foam combinations.

4.3 PCM filled with metal foams

Due to its intersected and connected heat transfer channels, metal foam is inserted into the LHTES system for enhancing the melting and solidification rates of PCMs. The numerical modeling of PCM solid-liquid phase change phenomenon filled in metal foams using enthalpy-based LBM could be categorized as the representative elementary volume (REV) scale modeling [128, 129, 130, 160] and pore-scale modeling [111, 161]. Tao et al. investigated the performance of metal foams/paraffin composite PCM in a LHTES cavity using LBM at REV scale [160]. They found that increasing the metal foam PPI (number of pores per inch) could enhance the conduction heat transfer, while the convective heat transfer is weakened. In addition, although decreasing the metal foam porosity could accelerate the PCM melting rate, the energy storage density of LHTES unit is dramatically reduced. Due to the above two tradeoffs, the optimum metal foam structure with the porosity of 0.94 and PPI of 45 is highly recommended. However, the REV scale modeling requires the use of some empirical relations, and the influences of metal foam morphology on energy storage rate are difficult to be analyzed. As a comparison, with the advantage of LBM for tackling complex boundary conditions, pore-scale modeling without using any empirical equations achieves the detailed fluid flow and heat transfer information inside the metal foams, so that the energy storage efficiency of LHTES unit could be further optimized. Ren et al. investigated the effects of metal foam characteristics on the PCM melting rate through pore-scale LBM modeling [111]. The transient PCM melting process inside the metal foams of different pore sizes at the Fourier number $\mathit{Fo}=0.04$ is plotted with respect to temperature field, melting interface, and fluid velocity vectors in Figure 13. It could be observed that the charging rate of PCM is accelerated when the pore size decreases from $d_p=1\mathrm{mm}$ to $d_p=0.75\mathrm{mm}$ at a porosity of ${\epsilon }_{\mathit{ave}}=0.95$. Furthermore, the temperature field in the melted region of LHTES with smaller pore size $d_p=0.75\mathrm{mm}$ is found to be more uniform than that of LHTES with pore size $d_p=1.0\mathrm{mm}$, which means that the thermal conductivity and its corresponding conduction heat transfer in LHTES unit are improved with a decreasing pore size (increasing PPI). Besides, they also concluded that an appropriate metal foam porosity should be chosen in the real engineering applications in order to balance the PCM melting speed and the energy storage density of LHTES unit. By reconstructing the microstructure of metal foam with X-ray computed tomography, Li et al. investigated the solid-liquid phase change phenomenon of PCM inserted with metal foams under different gravitational acceleration conditions [161]. The results indicated that the transition of the dominant heat transfer mechanism from convection to conduction occurs when the gravity gradually decreases. Due to the attenuated convection effect with decreasing gravity, the PCM melting development is dramatically hindered. Besides, they also concluded that the decreasing metal foam porosity could enhance the effective thermal conductivity of LHTES unit because of the extended heat transfer area. However, the above pore-scale modeling of PCM charging and discharging processes enhanced by metal foams is limited in two-dimensional cases. The three-dimensional pore-scale modeling, which is indispensable for optimizing the complicated metal foam structures, should be carried out in the future work.

Figure 13.

Transient PCM melting process filled in metal foams with different pore sizes [111]. (a) ${\epsilon }_{\mathit{ave}}=0.95$ and $d_p=1\mathrm{mm}$ and (b) ${\epsilon }_{\mathit{ave}}=0.95$ and $d_p=0.75\mathrm{mm}$.

4.4 PCM with heat pipes

The heat pipes are used to transfer heat between PCMs and heat transfer fluid (HTF), so that the charging and discharging of LHTES system could be accelerated. The configuration, arrangement, and number of heat pipes in the LHTES unit are essential for the PCM melting and solidification speed. Luo et al. applied LBM to study the convection melting in complex LHTES system with heat tubes [162]. The effects of inner heat pipe arrangement on PCM melting process are illustrated in Figure 14 with respect to temperature, flow, and phase fields at the Fourier number $\mathit{Fo}=3$, the Prandtl number $\mathit{Pr}=0.2$, the Stefan number $\mathit{Ste}=0.02$, and the Rayleigh number $\mathit{Ra}=5\times {10}^4$. For the case using centrosymmetric inner heat pipes, the conduction heat transfer is dominant because the inner heat pipe is surrounded by solid PCM at the melting temperature, so that its melting rate is faster than the LHTES system with inline or staggered heat tubes as displayed in Figure 15. In addition, there is no obvious difference between the melting rates of LHTES systems using inline or staggered heat pipes. Although using heat pipe individually could enhance the thermal performance of LHTES to some extent, other heat transfer enhancement techniques are usually coupled with heat pipes such as fins or nanoparticles to further improve the PCM charging and discharging rate.

Figure 14.

Temperature, flow, and phase fields at $\mathit{Fo}=3$ with different arrangements of inner tubes: (a) centrosymmetric, (b) inline, and (c) staggered for the Prandtl number $\mathit{Pr}=0.2$, the Stefan number $\mathit{Ste}=0.02$, and the Rayleigh number $\mathit{Ra}=5\times {10}^4$ [162].

Figure 15.

The total liquid fraction $f_l$ versus the Fourier number $\mathit{Fo}$ for different arrangements of inner tubes at the Prandtl number $\mathit{Pr}=0.2$, the Stefan number $\mathit{Ste}=0.02$, and the Rayleigh number $\mathit{Ra}=5\times {10}^4$ [162].

4.5 PCM with hybrid heat transfer enhancement techniques

In this section, the melting and solidification of PCMs enhanced using combined heat transfer enhancement techniques modeled by LBM are discussed, and the effectiveness of different approaches for improving the heat transfer capability of LHTES unit is compared. Huo and Rao carried out an investigation of NEPCM solid-liquid phase change process in a cavity with a heat pipe and separate plates using LBM [163]. It was found that the NEPCM of the case with separate located in the middle of cavity melts fastest because of the weakened heat accumulation. The results also showed that when the location of separate is less than 0.3, the melting rate of NEPCM is even slowed down due to the heat accumulation around the separate plate. Gao et al. developed an enthalpy-based MRT LBM model with a free parameter in the equilibrium distribution function for solid-liquid phase change in porous media and conjugate heat transfer with high-computational efficiency and stability [164]. Then, they investigated the PCM melting process in porous media with a conducting fin, and the results indicated that the heat transfer rate could be further improved by adding a fin into the porous medium. It was also observed that the vertical position of the fin has no remarkable impact on the PCM melting speed when there exists porous media. Jourabian et al. investigated the constrained ice melting around a cylinder using three heat transfer enhancement techniques by LBM [165]. They pointed out that adding nanoparticles may increase the dynamic viscosity of the base PCM, which has a negative effect on natural convective heat transfer. Besides, the tradeoff between consolidating the conduction heat transfer and weakening the convective heat transfer through decreasing the metal foam porosity needs further attention and investigation. Ren et al. completed a pore-scale comparative study of nanoparticle-enhanced PCM melting in a heat pipe–assisted LHTES unit with metal foams [166]. The melting process of NEPCM in a LHTES cavity filled with metal foams using a heat pipe at the Fourier number $\mathit{Fo}=0.06$, pore size $d_p=1\mathrm{mm}$, and heat pipe radius $R=2\mathrm{mm}$ is shown in Figure 16. Although different combinations of nanoparticle volume fraction Φ and metal foam porosity ${\epsilon }_{\mathit{ave}}$ are used for the LHTES unit, it should be noted that the volume fraction of pure PCM is kept at 90%, so that the energy storage capacity is unchanged for a fair comparison. It could be observed that the NEPCM with less nanoparticle volume fraction and lower metal foam porosity melts faster than that with more nanoparticles in the higher porosity metal foam. This finding indicated that using metal foams is more effective than adding nanoparticles for enhancing the charging rate of NEPCMs. They also found that there exists an optimum heat pipe radius for achieving the best energy storage rate in LHTES unit. As discussed in the above sections, LBM is demonstrated to be appropriate for simulating the charging and discharging processes of PCMs in LHTES system with different kinds of heat transfer enhancement technologies. Under this circumstance, the enthalpy-based LBM will definitely play a more significant role in the future research of thermal energy storage using PCMs.

Figure 16.

Transient nanoparticle-enhanced PCM melting process filled in metal foams of pore size $d_p=1\mathrm{mm}$ and heat pipe radius $R=2\mathrm{mm}$ [166]. (a) $\mathrm{\Phi }=0.09$, ${\epsilon }_{\mathit{ave}}=0.99$ and (b) $\mathrm{\Phi }=0.01$, ${\epsilon }_{\mathit{ave}}=0.91$.

5. Summary

In this chapter, different heat transfer enhancement techniques of PCMs for LHTES unit are discussed and compared. As the numerical modeling plays a significant role in clarifying the mechanism of complicated physical processes, the mathematical models for fluid flow of liquid PCMs and the PCM solid-liquid phase change phenomenon are presented. In order to investigate the PCM charging and discharging processes enhanced by nanoparticles, fins, or metal foams, the empirical relations for nanofluids and the mathematical formula for conjugate heat transfer are shown. The development history of the lattice Boltzmann method for solid-liquid phase change problems is carefully reviewed, and the enthalpy-based multiple-relaxation-time LBM is discussed in detail due to its simplicity and robustness. Besides, the implementation of GPU computing is briefly discussed to accelerate the computational efficiency of LBM modeling. Then, the applications of LBM modeling in LHTES system with different heat transfer enhancement approaches are presented, which demonstrate that the mesoscopic and highly parallel LBM is powerful for understanding the melting and solidification processes of PCMs.

Acknowledgments

This work was financially sponsored by the National Key Research and Development Program of China (No. 2017YFB0102703), the National Natural Science Foundation of China (No. 51536003), 111 Project (grant number B16038), and the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (No. 51721004).