Latent Heat Thermal Storage of Nano ‐ Enhanced Phase Change Material Filled by Copper Foam with Linear Porosity Variation in Vertical Direction

: The melting flow and heat transfer of copper ‐ oxide coconut oil in thermal energy storage filled with a nonlinear copper metal foam are addressed. The porosity of the copper foam changes linearly from bottom to top. The phase change material (PCM) is filled into the metal foam pores, which form a composite PCM. The natural convection effect is also taken into account. The effect of average porosity; porosity distribution; pore size density; the inclination angle of enclosure; and nanoparticles’ concentration on the isotherms, melting maps, and the melting rate are investigated. The results show that the average porosity is the most important parameter on the melting behavior. The variation in porosity from 0.825 to 0.9 changes the melting time by about 116%. The natural convection flows are weak in the metal foam, and hence, the impact of each of the other parameters on the melting time is insignificant (less than 5%).


Introduction
Heat transfer enhancement is a major challenge in various engineering fields, including heat exchangers, solar collectors, building heating, ventilation, conditioning, bio-and chemical reactors, and others [1][2][3]. The solution to this aforementioned problem can be achieved using different approaches. One of them is an introduction to the considered system, metal foam saturated with nano-enhanced phase change material. Such an approach includes different passive heat transfer intensifiers such as porous metal foam, phase change material, and low concentration nanoparticles. Numerical and experimental analyses of the industrial systems with such intensifiers separately have been presented by different authors [4][5][6][7][8]. Thus, in an experimental study, Ryu et al. [4] considered a nanoliquid supply layer with a metal foam for heat transfer augmentation in micro heat spreaders. Delgado-Diaz et al. [9] proposed several key performance indicators such as the normalized heat transfer performance coefficient and compactness degree. They then examined the performance of several enhancement techniques, including fins, composite metal foams, and macro capsulation techniques. The results showed that the composites show promising potential for latent heat thermal energy storage applications.
Feng et al. [5] experimentally studied thermal convection from a heat sink based on the copper foam fins under constant pore density and porosity. It has been revealed that the heat transfer augmentation can be achieved with a rise in foam height but with a limit. Numerical analysis of the metal foam influence on heat transfer and fluid flow in a finned heat sink under the heat-generating unit was conducted by Astanina et al. [6]. They ascertained that using three fins is an optimal number reflecting the reduction in heater temperature for the analyzed case. Liang et al. [7] scrutinized experimentally and numerically the latent heat storage combined with a flat micro-heat pipe array-metal foam system. Such a combination allowed them to improve the heat transfer performance of conventional heat storage systems. Moreover, the authors showed that heat conduction has a crucial influence on the process during the initial level. Experimental investigation of the opportunities of metal foam in a finned heat sink was performed by Shen et al. [8]. It was demonstrated that the addition of metal foam essentially decreases the total heat resistance by 25.5% compared to conventional finned heat pipe radiators.
It is well-known that the addition of solid nano-sized particles to conventional heat transfer liquids (water, engine oil, ethylene glycol, and others) allows for an increase in thermal conductivity that can be the reason for energy transport augmentation [10,11]. Some research has been performed for the case of nanofluid circulation in metal foam inserts [12][13][14][15]. Xu and Xing [12] examined the nanoliquid behavior in a chamber filled with a metal foam using the lattice Boltzmann method. They found that energy transport augmentation can be achieved by employing the highly conductive metal foam and nanosuspension.
Qi et al. [13] investigated the ferroliquid circulation experimentally in an inclined chamber with a metal foam under Lorentz force influence. The authors showed that growth of the nanoparticle mass fraction characterizes an augmentation of the Nusselt number but only up to a critical value; after that, this number reduces with the nanoparticles mass fraction. For the considered case, the optimal value of the nanoparticles mass fraction was 0.3%. Analysis of the nanoliquid behavior in a metal foam using the local non-equilibrium approach with the Brownian diffusion and thermophoresis influence was performed by Xu et al. [14]. It was ascertained that the Nusselt number can be increased with a rise in the nanoparticle concentration or diminution of the metal foam porosity. Ghaneifar et al. [15] studied the hybrid nanofluid influence on heat sink performance containing the multilayered copper foam using constant porosity and various particle diameter model and vice versa. The authors revealed that using such a layer allowed them to enhance the energy transport performance and pumping power.
Nowadays, there are many published papers on phase change material applications in different engineering systems due to moderate melting temperature, high latent melting energy, and high heat capacity [16,17]. In addition, some papers have been published on the phase change material behavior in a metal foam structure. El Idi and Karkri [18] investigated the natural convection of phase change material melting/solidification in aluminum foam within a square enclosure under the influence of temperature-dependent wall heat flux. They showed that energy removal on the wall has an essential impact in the case of sinusoidal time law for the heat flux compared with the case of fixed heat flux. Such an influence becomes greater in the case of the solidification stage in comparison with the melting one. Yang et al. [19] analyzed the effect of phase change natural convection in metal foam-filled enclosures by altering the enclosure inclination angle and aspect ratio. The results showed that the impact of inclination angle on the melting rate was minimal. However, the aspect ratio could significantly affect the melting rate. Iasiello et al. [20] examined the impact of various design parameters such as number of pores per inch (PPIs) and porosity on the melting heat transfer in an enclosure filled with aluminum foam. The results showed that a decrease in the porosity of metal foam could drastically enhance the melting process, but the pore-density could alter the melting process marginally. Zhang et al. [21] explored the melting of paraffin in copper foam in an enclosure containing a heated solid block. They performed a nondimensional analysis and controlled the strength of natural convection effects by Rayleigh number. They found that an increase in Rayleigh number boosts the natural convection at upper parts of the cavity, but it has a minimal effect in the bottom parts. Thus, the Rayleigh number could not change the overall melting time notably.
This brief review shows that analysis of nano-enhanced phase change material melting/solidification within the chamber filled with metal foam plays a crucial role in many engineering systems. Therefore, the aim of the present research is to numerically simulate nano-enhanced phase change material in a porous metal foam enclosure under linear porosity variation in the vertical direction. The variable porosity could contribute to conduction heat transfer in solid regions while allowing some degree of convection heat transfer in liquid regions.

Mathematical Model
A representative schematic of the tilted inhomogeneous porous enclosed medium is displayed in Figure 1. As depicted in Figure 1, the porosity is a linear function in the ydirection with positive or negative gradients. The inhomogeneous pores are occupied by a mixture of bio-based coconut oil and homogeneous dispersed CuO nanoparticles. The solid matrix of the porous medium is copper metal foam. Table 1 tabulates the properties of the PCM, CuO, and solid matrix.
The change of density from solid (920 kg/m 3 ) to liquid (914 kg/m 3 ) is just 0.65%. Thus, coconut oil has a very small density variation. This variation could be important for mechanical issues and may require a small void gap in the enclosure design to allow for liquid expansion during melting. However, the effect of density variation is neglected in the fluid motion and heat transfer model. The size of the 2D enclosure is 40 mm. The left side surface of the enclosure is set at a high constant-temperature of θh, and the other surfaces are insulated. The porosity function is as follows: where 0  and L  are the porosities at y = 0 and L, respectively. The average porosity, avg  , can be defined by integrating the porosity all over the domain. The melted nanoenhanced phase change material (NePCM) is considered laminar and Newtonian. In addition, the volume change in the NePCM during the phase change process is assumed to be zero. The enthalpy-porosity technique is employed to simulate melting flow. The balance equations for mass, momentum, and energy are as follows: where Xi is the ith coordinate such that X1 = x and X1 = y. Vi is the velocity in the ith coordinate.
in which Since the porosity is a spatial function, the corresponding permeability, i.e., η, is also a spatial variable [24]: where Ol is the diameter of fiber and Op is computed using the pore density. Considering a local thermal equation for PCM and metal foam, the conservation of energy can be written as follows: in which where lnp and snp refer to the liquid and solid NePCMs. cmf introduces the properties of the copper metal foam, and eff denotes the effective properties of the porous medium saturated with NePCM.
The thermal conductivity of the liquid/solid NePCM is defined as the following [25,26]: where z = (1 − ε)/3. In order to calculate the properties of the mixture of coconut oil and dispersed nanoparticles, the below-expressed relations are employed: Density: where np and pp introduce the properties of the NePCM and pure PCM, respectively. na denotes the properties of the nano-additives, and VFna is the nano-additives concentration. Viscosity: Thermal expansion coefficient: Thermal conductivity: Effective heat capacity: The following boundary and initial conditions are applied to the problem: On the hot surface: On the insulated vertical surfaces: On the insulated horizontal surfaces: Initial conditions: The total energy stored capacity in the NePCM domain during melting flow is expressed as follows: The first term of the above equation denotes the sensible energy, and the second term shows the latent energy stored in the PCM. The liquid fraction of the NePCM is as follows:

Numerical Approach and Grid Check
The finite element method (FEM) is known as a potent technique for solving partial differential and integral equations. Equations (2), (3) and (5) along with the pertinent initial and boundary conditions are solved by employing user-defined FEM codes. The Galerkin-weighted residual approach is employed to create the weak form of the controlling equations. This method applies the Lagrange finite elements with various orders to obtain the object variables in the non-overlapping zones of the computational region. The weighted mean of the residuals, Rm, while calculating the object variables is zero for the whole domain.
where Wf is the weight function. Nonlinear residual formulas are integrated by utilizing the second-order Gaussian-quadrature technique. The backward differentiation formula (BDF) [27], as an automatic matter, is employed to control the time step. To solve the residual formulas, a Newtonian damping with the value of 0.8 in PARallel DIrect SOlver (PARDISO) is selected. Moreover, the relative tolerance of 10 −3 is set for the convergence criteria [28][29][30].
The whole region of NePCM is divided into rectangular subdomains, namely finite elements. To assess an appropriate mesh that can simultaneously ensure the accuracy of the numerical results and lower computational cost, the grid sizes are increased from 125 × 125, 150 × 150, 175 × 175, 200 × 200, and 250 × 250. Figure 2a depicts the structured grid 150 × 150 adopted for the computations. As seen in Figure 2b To use the present model for further calculation, the model is verified and validated with the numerical and experimental work in [22,31]. In the numerical verification, isotherms and streamlines in an inhomogeneous porous enclosed medium of the work conducted by Xiong et al. [32] and the current model are compared, as illustrated in Figure  3. In this comparative analysis, free convection heat transfer in an enclosure with the inhomogeneous porosity distribution is analyzed. High and low temperatures of θh and θc are imposed on the left and right surfaces of the domain. However, the other walls are adiabatic. In the comparison presented in Figure 3, the porosity distribution is a linear function in y-direction. Herein, the nondimensional parameters of Rayleigh are 10 4 , Darcy 10 −1 , and Prandtl 6.2. This verification clearly depicts that the present model is in excellent agreement with the computations of the study carried out by Xiong et al. [32]. In the second comparative analysis, the results of the experiments carried out by Al-Jethelah et al. [22] were re-stimulated with our developed model for code validation. In [22], while the other walls of the enclosure were well insulated, the left wall was heated continuously with constant heat flux. The copper foam with high porosity (ε = 0.92) was employed. The permeability of the copper foam calculated using Equation (4) was assessed to be 3.3142 × 10 −7 by the authors. As depicted in Figure 4, the obtained melt fraction of the present model and those presented by [22] are in good agreement, and similar trends for conduction and convection modes could be seen.

Results and Discussion
The average porosity ( 0 . Here, an orthogonal table for a combination of parameters is considered, as shown in Table 2. This table searches for the alterable combination of parameter efficiency. The numerical tests were done for each case of Table 2, and the required time for melting of 90% and corresponding stored energy and power are reported. As seen, case 6 of this table shows the minimum 90% charging time of 313 s. Thus, this case is considered the best design. Table 3 shows that the largest melting time is for case 5 with a melting time of 953 s. This case has a maximum porosity of 0.9. Thus, the relative difference between the best case (313 s) and the worst case (935) is 205% compared to the best case (100 × (935-313)/313).
Twenty more systematic cases were executed, according to Table 3, to explore the influence of volume fraction of nanoparticles, porosity, linear distribution of porosity, and pore size on the storage design's melting behavior. This table shows that the increase in concentration of nanoparticles reduces the melting time. Growth of the inclination angle or pores density (PPI) slightly raises the melting time, while a rise in porosity notably increases the melting time. The variation in porosity distribution (a) leads to a nonmonotonic variation of melting time.
The results of Table 3 show that the variation in volume fraction of nanoparticles from zero to 4% could change the melting time by 3.4%. The variation of porosity from 0.825 to 0.90 induces 116% changes in the melting time. Alteration of the nonuniform porosity parameter (a) from −4 to 4 could change the melting time by 2.3%. Finally, the variation in inclination angle (γ) and pore density (PPI) changed the melting time by 0.5% and 4.4%, respectively.  The streamlines can be seen in Figure 5. Interestingly, most of the PCM melted in the case εave = 0.8 and at 450 s, while for the case of εave = 0.9 and the same time-step, only about a half of the enclosure melted. The isotherms show almost a linear variation in the enclosure. Here, t = 50 s is the beginning of melting and most of the cavity is solid. Thus as mentioned, the heat transfer is mostly conduction dominant. For low-porosity cases, the isotherms are generally farther from each other than the cases with high porosity. This is since a low porosity means more metal foam and better composite effective thermal conductivity. The lower the composite thermal conductivity, the lower the temperature gradients. Figures 5 and 6 depict the melting maps and isotherms in the enclosure for various porosity values. The results are reported at three time-points: 50 s, 250 s, and 450 s. The first 50 s could be considered a conduction dominant heat transfer regime. As seen, the increase in porosity reduces the overall molten region. This is since a higher porosity means lower metal foam and smaller composite thermal conductivity. Moreover, for this case, the porosity parameter (a) is negative (a = −2), and hence, as height (y coordinate) increases, the porosity decreases locally. Thus, the amount of metal at the high cavity regions is more and the local composite thermal conductivity is also high. As seen, the top region melted faster than the bottom region. The results are for an inclination angle of π/2, where the gravity vector acts toward the hot wall. The natural convection flow tends to improve the melting heat transfer in the middle of the enclosure due to the direction of the buoyancy force. Figure 7 illustrates the MVF and stored energy (ES) in the enclosure for various porosity values. This figure is in agreement with the contours of Figures 5 and 6 showing, that the lower the porosity, the higher the melting rate and stored energy.  Figures 8 and 9 show the impact of the porosity parameter (a) in the melting map and isotherms. As seen, by variation of a, some degrees of isotherm deflections could be seen. For the cases of a = −4 and −2, the isotherms are close to each other at the bottom and have more distance at the top. This is because, for negative values of a, the porosity is high at the bottom (low composite thermal conductivity). A low composite thermal conductivity leads to high-temperature gradients. In contrast, for the cases with positive a, the behavior is reversed. It should be noted that, for these cases, gravity acts toward the hot wall, and hence, the buoyancy force is perpendicular to the hot wall and tends to melt the PCM mostly from the center. Hence, when a is small, the general impact of the negative and positive values of a on the MVF and ES could be almost the same. When a is large, it could lead to a significant geometrical variation in the melting interface and could change the melting behavior. Figure 10 shows the MVF and ES for various values of a. This figure shows that the variation in a could induce minimal impact on the MVF and stored energy.     The influence of inclination angle (γ) on the melting process is investigated in the contours of Figures 11 and 12. The inclination of the cavity changes the buoyancy forces and the direction of streamlines. For example, for the cases of γ = 0 and π/2, the streamlines are clockwise, while for the cases of γ = 3π/4 and π, the streamlines are counterclockwise. Figure 13 is plotted to depict the impact of the inclination angle on the melting fraction and stored energy. These figures show that the inclination angle induces a minimal impact on the meting rate and stored energy. The reason for this is that the convection heat transfer is limited in a porous medium and particularly in a case with low porosity. Thus, the variation in inclination angle, which changes the natural convection flows, could also induce minimal impacts. Maybe using a larger enclosure could pronounce such effects. However, for a large enclosure, the overall melting time would also increase, which could be impractical for fast-charging energy storage systems. The effect of variation in nanoparticle concentration and pores density (PPI) on the melting process was also investigated. However, due to a significant amount of metal foams and minimal importance of natural convection effects, the results showed that the nanoparticles could add minimal impact on the phase change rate. Such graphs were not reported for the sake of brevity.

The Key Findings
The flow and heat transfer in an enclosure filled with a nonuniform metal foam was investigated. The porosity was changed in the y-direction using a controlling parameter a, while the average porosity was fixed. The governing equations and phase change properties for phase change heat transfer in the metal foam were introduced in the form of partial differential equations and were solved using the finite element method. The storage unit's melting behavior was studied for various values of average porosity, nonuniformity parameters, and inclination angles. The presence of metal foams could significantly improve the melting rate. The variation of porosity from 0.825 to 0.9 could change the melting time by 116%

Practical Significance/Usefulness
The presence of metal foam can improve heat transfer by enhancing the conduction heat transfer mechanism. However, as the melting process advances, the number of liquid PCM rises and natural convection flows can be formed. This is while the presence of a metal foam resists fluid motion and tends to suppress the natural convection flows. The presence of nanoparticles can also improve conduction heat transfer, but they increase the dynamic viscosity of liquid PCM, which acts against the fluid motion. Here, a nonuniform porosity variation was proposed to support the conduction heat transfer mechanism in the early stages of heat trasfer, where the PCM is solid. Then, the porosity can be increased where the natural convection flows commences. It was found that the nonuniformity parameter could change the isotherm and melting behavior in the enclosure. However, due to the configuration of the enclosure, the variation in porosity in the y-direction did not induce a significant impact on the melting rate. The impact of volume fraction of nanoparticles, the pore density, and inclination angles on the melting heat transfer was also small. Indeed, these parameters could impact the natural convection flows directly. This is while the natural convection flows in the metal foam were limited. Hence, the control parameters, which could influence natural convection, could only induce minimal alteration of the final melting process. The impact of each of these design parameters on the melting time was less than 5%.
It can be concluded that the porosity is a major design parameter and should be taken into account for acceleration of melting heat transfer. Then, the heat transfer can be slightly improved by using other methods.
One can also conclude that the porosity variation and using nanoparticles could better impact the melting process for large-scale thermal energy storage systems. In a large-scale system, the molten region can be extended and the distance between the heated surface and solid PCM could be large. Therefore, heat should be conducted through long distances to reach the cold PCM. While there is an extended molten region, the hot liquid PCM can easily move from the hot surface and reach the cold PCM. Thus, in large-scale systems, the natural convection flows could be significant. The large-scale systems were not in the scope of the present research, and hence, they were not investigated here.
In the continuation of this study, some suggestions can be made. As discussed in detail, this study provides the melting process in a copper foam with linear porosity variation in the vertical direction. Moreover, it was assumed that the PCM and metal foam are in local thermal equilibrium. Hence, in an upcoming study, the local thermal nonequilibrium condition can be applied to the copper foam with linear porosity variation. In addition, the porosity variations can be considered nonlinear.