Lattice Boltzmann simulation on thermal performance of composite phase change material based on Voronoi models∗

Meng-Yue Guo(郭孟月) Qun Han(韩群) Xiang-Dong Liu(刘向东) and Bo Zhou(周博)

1Key Laboratory of Energy Thermal Conversion and Control of Ministry of Education,School of Energy and Environment,Southeast University,Nanjing 210096,China

2College of Electrical,Energy,and Power Engineering,Yangzhou University,Yangzhou 225127,China

3Jiangsu Key Laboratory of Micro and Nano Heat Fluid Flow Technology and Energy Application,School of Environmental Science and Engineering,Suzhou University of Science and Technology,Suzhou 215009,China

Keywords: metal foam,Voronoi,melting,heat transfer enhancement

1. Introduction

Clean and renewable energy has received considerable attention since the fossil energy shortage and environmental pollution problems. However, the renewable energy including wind energy, wave energy and ocean thermal energy,the renewable energy has drawbacks of instability and nonuniformity in space and time distribution. To make up for these deficiencies,the latent heat storage(LHS)technique,in which the PCM is used to store energy in fields of the building energy saving,[1]thermal storage system,[2]PV-efficiency enhancement,[3]new materials,[4]and electronic devices cooling,is adopted,[5]in that it has lager energy density and temperature stability.[6]

However, most of PCMs have drawbacks of low heat storage rate. To improve thermal performance, some methods have been developed in the past years. Among those methods, adding additives to form a composite PCM is always a direct and effective way, including phase change microcapsule,[7,8]metallic finned tubes,[9,10]metal powder,[11,12]nanoparticles,[13]and metal foam.[14]The metal foam has the merit of high surface-to-volume ratio, which is widely adopted. For example, long fins at the lower part and shorter fins at the upper part show a better melting behaviour in the triplex-tube.[10]Laraibet al.[1]made the NePCMs,and thus the most of the energy and smart textiles improve their thermal stability due to the application of the NePCMs.Ah[13]obtained the Al2O3nanomaterial by physically shaking and tested the thermal properties of the Al2O3nanoparticles based nanofluid. As expected,the heat conductivity of the improved nanofluid was enhanced by 20%and its heat capacity was increased by 7%. Sardariet al.[14]made a comparison between the PCM with and without metal,where discharging time decreases and the temperature is more uniform for composite PCM.These attempts have enhanced the efficiency of the LHS unit to some extent. Nevertheless, it is significant to notice that the immature production process and particle flow instability restrict the broad applications in metallic finned tubes,nanoparticles and metal powder. Relatively, the metal foam has a well-rounded production technology, large specific surface and light weight,which makes it attractive in the fields of phase-change heat transfer enhancement.

From a comprehensive review, numerous researches on the key factors of the metal foam/PCM composite have been conducted. In the investigation, the experimental and numerical work conducted by Zhanget al.[15]showed that the treeshaped structure has better thermal performances: higher heat transfer rate and more uniform temperature distribution. A review of recent work reveals that in the current researches,significant attention was paid to the influences of foam metals’macro-physical properties on melting and solidification performance. It is important to note that most of the researches are conducted primarily through experiment because the internal structure of foam metal is complex and difficult to construct by numerical modeling. Moreover, the detailed microscopic heat transfer mechanism between pores of the metal foam is still not clear. For the reconstruction of the metal foam,some appropriate methods have been provided,including pore overlap criterion,[16]fractal geometry,[17,18]and the Voronoi random model. The pore overlap criterion[16]generates the metal foam geometry by using random spheres with different pore radii, porosities, and spacings, which is unable to present the actual continuity between the metal skeletons. Furthermore,the fractal geometry[17,19]reconstructs the porous metal foam with random and multiscale characterization, while it fails to describe the irregular properties of the actual metal foam. Besides, in Ref. [11] great efforts were made to reveal the microstructure of the composite PCM and analyze some effects of the geometric parameters on the thermal performance. But the discontinuous structure of the metal additives is different from the actual metal frame. More significantly,the Voronoi random model[20]has been used to investigate the influences of density,thickness and imperfection of pore edges on the metal foam for the composite PCM due to the merits of flexibility, controllability, diversity, and high quality. Therefore,by using the Voronoi random model,metal foam/PCM composite with microstructure is generated.

The lattice Boltzmann method(LBM)has attracted many scholars’attention due to the merits of high calculation accuracy, stability, simple treatment for the complex boundaries,and high parallel efficiency compared with traditional computational fluid dynamics(CFD)method.[21–24]Heet al.[21]has reviewed the applications of the LBM in the field of the solid–liquid phase change in the past 30 years. Huo and Rao[23]developed a quasi-enthalpy-based LBM,which obtains temperature without iterations through two key points of“prediction”and“consumption”. For simplification of the thermal boundary treatments, Wanget al.[25]proposed an improved LBM to avoid treating special interfaces for the contact boundary conditions. Zhaoet al.[26]creatively proposed a generalized two-relaxation-time LBM which can reduce iterative calculation and improve the accuracy for the phase change problems.Moreover, the theoretical development of LBM provides an effective mean to study the melting behavior. Luet al.[27]presented an optimal lattice Boltzmann (LB) model which is effective in reducing the unphysical diffusion for the melting problems. Besides,numerical and experimental studies on the melting phenomenon,performed by Chenet al.,[28]indicated that the LBM results match experimental results well. The above researches have provided a great mathematical theory for studying metal foam/PCM composite. However, few researches are conducted on a pore scale and microscopic heat transport mechanism between pores in metal foams has not yet been clear. Therefore, it is significant to explore the melting performance of the composite PCM and relevant research results are helpful in further understanding the microscopic heat transfer mechanism in metal foams.

In summary,there have been several attempts to improve the melting performances by using metal foam. However,it is still unclear how the pore structure affects the melting heat transfer mechanism on a microscopic scale. In addition, the quantitative characterization of geometry structure of metal foam is far from satisfactory due to the irregularity,randomness and complexity of porous structure. To reveal the microscopic melting mechanism of metal foam/PCM,the porous structure is quantitatively characterized by the Voronoi method,and a transient melting model with natural convection is developed and simulated based on the LBM in this work.The LB model is presented and followed by the model validation in terms of phase change problems. Moreover,the role of metal foam in melting phase change process of the PCM is examined and analyzed. Using the Voronoi method, the effects of geometric parameters(including porosity,pore density,and irregularity)on the melting heat transfer characteristics are investigated quantitatively,which provides a useful insight into the melting heat transfer enhancement.

2. Characterization of composite PCM

To show the microstructure of metal foams in detail and reveal the microscopic melting mechanism on a pore scare,the Voronoi random model is used here to reconstruct the composite PCM and detailed construction algorithm is given below.

(i)Takingnoriginal points at different positions in a twodimensional (2D) area (i.e., the composite PCM zone) withN×Nnodes,the original points are generated by the random numbers and are not distributed on a single line.Then,Voronoi polygons withnareas centered at the above original points are generated. In this context,the pore density and irregularity are determined by the number and positional distribution of original points.

(ii) Constructing an irregular triangle network based on discrete original points,i.e., the three nearest points of the above original points form a triangle. As shown in Fig. 1(a),constructing a super triangle which includes all the original points firstly.Then inserting a scattered point(the point P)and making circumcircles which include the new scattered point(the triangle is called effective triangle).Furthermore,deleting the common edges of effective triangles,moreover,generating new edges with the new point optimally,finally,repeating the above steps until all points are inserted,the generated triangle network is the Delaunay triangle network and the apex of the triangle is the original point.

(iii)Making perpendicular bisectors of the three sides of all Delaunay triangles, the center of the circumscribed circle is obtained. As depicted in Fig.1(b),the Voronoi structures are generated by connecting centers of the circumscribed circles between adjacent triangles. After eliminating center points,denoising,filtering and constructing binary Voronoi polygons,the composite PCM is reconstructed. In Fig. 1(c), the blue area denotes PCM and the silver area refers to the metal foam network.

Fig.1. Building process of composite PCM by Voronoi method: (a)delaunay triangle network,(b)Voronoi structure,(c)reconstructed composite PCM.

Fig. 2. Construction of metal foam/PCM with different geometric parameters: (a)PPI=15,ε =0.90,c=0.22;(b)PPI=30,ε =0.90,c=0.22;(c)PPI=44,ε=0.90,c=0.22;(d)PPI=44,ε=0.75,c=0.22;(e)PPI=30,ε =0.85,c=0.11;(f)PPI=30,ε =0.85,c=0.32.

According to the above construction process,the detailed structure of the metal foam/PCM composite is distinguished by three primary parameters: pore density,porosity,and irregularity. Figures 2(a)and 2(b)show that the density of Voronoi polygon turns higher with pore density increasing. Besides,the porosity of metal foam is a significant parameter. Figure 2(c)and 2(d)show the metallic frameworks with different thickness result in different porsities.

Voronoi structures with different irregularities are obtained by controlling the positions of original points. Here,the irregularity of the metal foam can be expressed as[29]

whereAis the total area,cis the irregularity,dis equivalent diameter of the averaged area,anddiis the equivalent diameter of then-th cell. As depicted in Figs.2(e)and 2(f),with the increase in irregularity(i.e.,c),the structure of the metal foam becomes more random.

3. Melting model description

3.1. Physical problem

The practical metal foam possesses a complex structure which will consume plenty of time and budget in the process of the simulation. However,previous researches[25]show that the two-dimensional (2D) model can predict the phase change process appropriately with high computational efficiency. Figure 3 gives the 2D physical model for the investigation of the composite PCM. The size of the domain isL×L=15 mm×15 mm. The PCM(i.e.,the blue zone)is immersed in the metal foam(i.e.,the silver region). After comprehensively considering the economy cost and efficiency,the paraffin and aluminum are chosen.

Fig.3. Physical model of metal foam/PCM.

Table 1 gives the thermophysical parameters of the two materials. The initial temperature of the whole composite PCM is 20°C, while left side is maintained at 70°C. The other three surfaces are all regarded as adiabatic.

Table 1. Thermophysical parameters of materials.

3.2. Governing equations

The melting process is complicated due to the multiple heat transport mechanisms and the complex structure.To comprehensively study the melting process,some assumptions are made in this work below.

i)The liquid is incompressible and the Bossinesq hypothesis isapplicable;

ii)The viscous heat dissipation and the heat radiation are ignored;iii)The volume keeps unchanged under different phases;iv) The physical properties are constant for the different phases;

v)The melting temperature is almost constant.[28]

Based on the assumptions above,the governing equations can be expressed as

3.3. Lattice Boltzmann method

By using the LBM, the double distribution is employed.Through the force term and the equilibrium distribution, the velocity and temperature can be connected together. Among them,the density distribution function is expressed as

where the subscript f and s represent the PCM and the metal foam,respectively. The source termSriis described as

For the heat transfer in the melting process, there is a mutual coupling between the temperature and the liquid fraction and these parameters are calculated by the enthalpy-based method. At each time step, we calculate the enthalpy value from the macroscopic temperature. Then the new liquid fraction can be obtained and the temperature distribution function will be updated. Furthermore, the liquid fraction and temperature are calculated iteratively until the convergence has been satisfied. And the macroscopic parameters including the density,temperature and velocity are obtained fromρ=∑i fi,T=∑i gi,andu=∑i ei fi/ρ+Fδt/2.

The PCM with different phases can be qualitatively described by the liquid fraction,where the liquid fraction of the solid,liquid and mushy zone areγ=0,γ=1,and 0<γ<1,respectively.

Moreover, the treatment of the boundary is of significance for the simulation result. The bounce-back scheme is widely used.[33]However, the scheme is of first-order accuracy. Therefore, the non-equilibrium extrapolation scheme is adopted.[34]

For the velocity boundary nodex0, the velocity condition is given while the densityρ(x0,t)is unknown. And thexfis the nearest node (i.e. xf=xo+eiδt). For the nodexo,the distribution which includes the equilibrium and the nonequilibrium distribution functions is expressed as

In addition, to satisfy the conjugate heat transfer conditions,several available LB models have been proposed,which can recover the diffusion term in the energy equation without special treatments of the interface.[24,35]In this work, we use the interface treatment method with second-order accuracy.[36]And the interfacial temperature can be calculated from the microscopic distribution functions by this method.

3.4. Model validation

To ensure the accuracy,the grid independence test is performed first. Figure 4 shows the results of grid independence tests under grid systems with the numbers of nodes of 300×300, 600×600, and 900×900, respectively. The PCMs are under solid state first, the left side is kept at a high temperatureTh, and the right side stays at a low temperatureTc,whereTh=1.0 andTc=0. Besides, the horizontal sides are assumed to be adiabatic. The dimensionless Nusselt numbers of the left wall withPr=5.0,Ste=1.0,andRa=50000 under three grid systems are given,where the dimensionless Nusselt number is given by

It can be seen that the result of the case with nodes of 600×600 is in accordance with that of the case with nodes of 900×900. Considering the the computation cost, the grid system with nodes of 600×600 is adopted in the following simulation.

Fig.4. Grid independence test under various grid systems.

Furthermore,the results from the LBM and analytic solutions are compared by solving the conduction-dominated melting problem(i.e.the Stefan problem). Figure 5(a)shows that the left side is kept at a temperatureThwhile others are adiabatic boundaries. Moreover,initial temperature isT0,which is the melting temperature and lower thanTh. The analytical solution for conduction-dominated melting problem can be obtained by Rubinstein[37]and Huberet al.[38]

wherexmis the melting position. Figure 5(b) displays the melting front position evolution for the cases with different thermal diffusivities (α=0.033, 0.36, and 0.7) but the same Stefan number ofSte=1.0.Among them,the normalized timet∗=tα/L2,whereLrepresents the length of the system.

The numerical predictions by LBM and analytic solutions(AS)are in good agreement with each other while the relative error is only 4.9%.

Fig.5. Model validation of conduction melting problem: (a)Schematic diagram of the conduction melting,(b)comparison of the melting front.

For further validation,the results from the present model are compared with those from Ref. [39] in terms of the dimensionless Nusselt number. Figure 6 shows the physical model used for this comparison, the left boundary and right boundary are kept,respectively,at constant temperature ofThandTc, while other two walls are assumed to be adiabatic.Moreover, the computational domain is a grid system with 600×600 nodes and the numerical simulation is conducted withRa=105, andPr=1.0. Table 2 gives the results from the literature and the present model under various thermal conductivity ratios. The maximum relative error is 4.7%between the result from the present model and the Merrikhet al.’s result,[39]thereby indicating that present model is correct and reliable.

Fig.6. Schematic diagram of model validation in a square.

Table 2. Comparison between dimensionless Nusselt numbers from the present model and literature.

4. Results and discussion

4.1. Distribution of melting front and temperature

For analyzing the thermal performance,a comparison between cases with and without metal foam (e.g.porosityε=0.8)is carried out underRa=50000,Pr=5.0,andSte=1.0.Figure 7 displays the distributions of melting front that the solid-liquid interface of the pure PCM is a vertical line at the early stage and then presents an inclined shape with the melting time increasing (i.e.the upper part of PCM melts faster than the lower part). This phenomenon implies that the heat conduction plays a leading role at the early stage while the natural convection exerts a more significant influence at the later stage. However, the melting front interface of the composite PCM maintains a nearly straight line during the whole stage,which differs from that of the pure PCM.Moreover,the metal foam/PCM has a faster melting rate, and a melting rate enhanced by 69%compared with that of the PCM without metal foam.When the composite PCM melting process is almost finished,the pure PCM melting process is completed by less than half. It can be concluded that the existence of metal enhances melting capability and the heat conduction plays a major role in the composite PCM.

Fig.7. Distributions of melting front for PCM with and without metal foam.

Fig.8. Distributions of temperature for PCM with and without metal foam.

Figure 8 displays distributions of the temperature for PCM with and without metal foam, indicating that the metal foam accelerates the heat transport rate with respect to the heat transport rate of the pure PCM without any additives. Moreover,when the natural convection develops,the temperature in the upper part of the melting front is high,but low in the lower part for the pure PCM. Nevertheless, the metal matrix with high thermal conductivity enhances the heat transport and the heat first pass through the metal skeleton. The heat transports almost horizontally in the composite PCM since the heat conduction exerts a larger influence.

4.2. Analysis of metal foam configurations

4.2.1. Effects of pore density

The pore density is a significant structural parameter of metal foam which can affect the melting behaviors. Figure 9 exhibits the distributions of the melting front and temperature with different pore densities to discover the effect of the pore density,with parameters of porosityεand irregularityckeeping unchanged,i.e.ε=0.90 andc=0.22. It is observed that the melting front evolution and the thermal response become quicker as the pore density increases. It implies that the larger pore density results in higher contact surface area between metal frame and PCM. And larger contact area conduces to the enhancement of melting heat transfer performance.

The times taken to melt the metal foam/PCM with different pore densities are depicted in Fig. 10, indicating that the time decreases obviously as the pore density increases. Moreover,as the pore density increases fromPPI=15 toPPI=44,the melting duration time decreases by 17.8%. The results imply that larger pore density is beneficial to the heat transport of the composite PCM normally.

Fig.9. Melting behaviors of metal foam/PCM with different pore densities.

Fig. 10. Times taken to melt composite PCM with different pore densities.

4.2.2. Effects of porosity

The porosity seriously influences the pore structure of the metal foam. Figure 11 shows the distribution of phase change region and temperature with various porosites to study the effects of porosity on the melting behaviors. Here, the pore density isPPI=30 and the irregularity isc=0.22. The result shows that the melting rate decreases as the porosity increases. The reason is that the conduction exerts a major influence in the heat transfer mechanism, and a lower porosity contributes to a larger metal foam fraction, which enhances the melting heat conduction performance. However,it should be noted that the high porosity reduces the heat capacity of metal foam/PCM.For this reason,it is important to construct the metal foam with a suitable porosity for making a balance between melting rate and heat storage capacity.

Fig.11. Melting behaviors of metal foam/PCM with different porosities.

Fig.12. Times taken to melt composite PCM with different porosities.

Figure 12 shows the times taken to melt the metal foam/PCM with different porosities. It shows that the porosity plays a major role in the melting performance: as the porosity increases, the melting duration time increases evidently.Moreover,the melting time is increased by 6.2%as the porosity increases from 0.75 to 0.8, while there is an increase of 16.1% as the porosity changes from 0.8 to 0.85. It indicates that a higher porosity weakens the melting rate more obviously. Based on the requirement for the balance between the heat transfer rate and heat storage density, a suitable porosity is seriously important in improving the thermal performance of the composite PCM.

4.2.3. Effects of irregularity

The irregularity is a crucial characteristic parameter for the metal foam. Figure 13 exhibits the melting behaviors of metal foam/PCM with different irregularities. In this subsubsection, the pore density isPPI=30 and the porosity isε=0.85. As shown in this figure, the lower irregularity results in a more uniform metal foam in terms of sizes and shapes, where the metal foam with a larger irregularity has cells with different sizes. Moreover, it is clear that melting rate decreases as irregularity increases under the same pore density and porosity. The cause is that the hampered and circuitous circled heat transfer path reduces the thermal performance in the composite PCM.Besides,interestingly,the metal foam/PCM withc=0.32 melts quicker in the middle part than in the upper and lower part in Fig. 13(a), where the distribution of metal framework is non-uniform and the middle part has a dense skeleton which enhances the heat conduction in this part.

Fig.13. Melting behaviors of metal foam/PCM with different irregularities.

Figure 14 shows the times taken to melt metal foam/PCM with different irregularities. From this figure, it can be obtained that the low irregularity can improve the ability of metal foam/PCM to transfer heat. And interestingly, the irregularity decreases from 0.32 to 0.21, the melting time is reduced by 20%. Therefore, under the same conditions including the porosity and pore density,the metal foam/PCM with low irregularity is recommended to be used to enhance the heat storage rate.

Fig.14. Times taken to melt metal foam/PCM with different irregularities.

5. Conclusions

In this paper, the melting behavior of PCM with metal foam generated by the Voronoi method is studied. An LBM is implemented and validated on a pore scale. Moreover, the evolution of phase and temperature distribution of the PCM with and without metal foam are analyzed. Besides,the influences of metal foam structure parameters on the melting performance are studied,including the pore density,porosity and irregularity. The main conclusions are drawn below.

(I) The metal foam is beneficial to the strengthening of the heat transport. Compared with pure PCM, the metal foam/PCM possesses quick evolution of melting front morphology and heat transport rate.

(II)The major heat transport mechanism of pure PCM is different from that of metal foam/PCM.For the pure PCM,the major heat transport mechanism is heat conduction first, and convection later. Nevertheless,the conduction exerts a crucial influence on the whole melting process since the convection is weakened for the metal foam/PCM.

(III) The structural characterization of metal foam remarkably influences the melting behaviors. The decrease of the porosity can improve the heat transport rate of composite PCM,while heat storage density will decrease accordingly.With the pore density increasing fromPPI=15 toPPI=44,the melting time decreases by 17.8%;with the irregularity decreasing fromc=0.32 toc=0.21,the melting time decreases by 20%.Therefore,metal foam is recommended to have a high pore density and thus presenting more uniform pore structure.