**1. Introduction**

The design of thermal energy storage (TES) systems has improved significantly in recent years. Phase change materials (PCMs) have been implemented in various practical applications such as buildings materials [1], air heating and cooling in buildings [2], cooling of electronic components [3], recovering low-temperature industrial waste heat [4], and automotive applications [5].

Most of the advances in TES system design have been around improving the response time of storage units during the charging and discharging process and the synthesis of new composite PCM materials. For instance, using expanded graphite additives [6], tree-like

**Citation:** Ghalambaz, M.; Mehryan, S.A.M.; Shirivand, H.; Shalbafi, F.; Younis, O.; Inthavong, K.; Ahmadi, G.; Talebizadehsardari, P. Simulation of a Fast-Charging Porous Thermal Energy Storage System Saturated with a Nano-Enhanced Phase Change Material. *Energies* **2021**, *14*, 1575. https://doi.org/10.3390/en14061575

Academic Editor: Alessandro Mauro

Received: 27 January 2021 Accepted: 6 March 2021 Published: 12 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

fins [7], stepped fins [8], and heat pipes [9] are some of the recent techniques. Moreover, thermal conductive structures such as three-dimensional (3D) porous diamond foams [10] and continuous diamond–carbon nanotube foams [11] have also shown promising performance.

Many researchers attempted using a frame of thermal conductive metals and porous metal foams, saturating them with phase change material. The thermal conductive structure of the metal foam carries out the heat to and from PCM inside the pores. Thus, the composite PCM–metal foam channels the heat from the heat sources and improves the thermal conductivity. Open metal foams allow some degree of natural convection heat transfer in molten areas of the TES unit. The natural convection heat transfer is an important phenomenon that transfers thermal energy by advection in a molten region. Considering open metal foams and natural convection effects, Talebizadehsardari [12] utilized a composite of metal foam and PCMs and designed an air heater for domestic application. The results indicated that the geometrical shape of air passages induces dominant effects on the PCM unit's thermal behavior. Sardari et al. [13] used an aluminum foam–PCM composite TES unit to absorb a wall-mounted radiator's heat during off-peak loads. Later, the TES unit releases the absorbed heat when the primary heating system turns off. The authors showed that using the metal foam shortens the charging time by 95% compared to a simple PCM. Interestingly, the increase in metal foam porosity (97%) produces a positive influence and reduces the TES unit's charging/discharging time.

Zhao et al. [14] compared the advantages of using fins with those of using metal foams in reducing the melting time of a PCM in a shell-and-tube shape TES unit. They found an optimum design for the fins. They then compared the melting time of using the optimum fin with a case made using the same amount of metal foam instead of the fins. They reported that properly designed fins could be as advantageous as metal foam.

Using nanoparticles is another approach to synthesize nano-enhanced phase change materials (NePCMs). Zhang et al. [15] dispersed copper oxide nanoparticles in RT28 PCM to improve the charging heat transfer rate in a wavy-channel TES unit. The presence of nanoparticles improved the heat transfer rate and reduced the charging time. Bondareva et al. [16] examined the melting of an alumina–paraffin NePCM inside a copper radiator. This investigation showed that the nanoparticle presence increases the viscosity of molten PCM, limiting the mobility of liquid PCM.

On the other hand, nano-additives enhance the thermal conductivity of the PCM. Thus, using NePCM was most advantageous at an early melting heat transfer stage, where the conduction regime was dominant. It can be concluded that the advantage of using nanoparticles mainly depends on the internal structure of the TES enclosure and dominant mechanisms of heat transfer.

All of the above studies used a uniform metal foam to improve the thermal conductivity of composite PCM. Mahdi et al. [17] employed a cascade (multiple segments) of PCM–metal foam composites to enhance a TES unit's discharging time. They also employed 5% nanoparticles to further enhance the thermal conductivity of the system. Using the multiple segments and nanoparticles, they reduced the discharging time by 94%.

The literature studies showed that the natural convection effects and the thermal conductivity of composite PCMs could significantly influence a TES unit's phase change behavior. A dense porous structure (low porosity) increases the composite thermal conductivity but suppresses convection flows. This is while a high porosity metal foam allows convection flows. Thus, a nonuniform porous metal foam could be advantageous from both points of view. The present study investigates the melting process of a PCM embedded in a nonuniform metal foam. The influence of using various concentrations of nanoparticles, porous densities, pores sizes, and porous gradients on the full charging time of the TES unit are addressed. The Taguchi optimization approach was employed to systematically find an optimum design of nonuniform porosity. In particular, the study addressed the following research questions: (1) Does a nonuniform porosity reduce the melting time for a fixed amount of metal foam? (2) What is the impact of nanoparticles and pore size on the melting time of TES?

#### **2. Governing Equations and Boundary Conditions**

A schematic view of the physical and computational domains with dimensions of the latent heat thermal energy storage unit is depicted in Figure 1a,b. The TES unit is a cylinder with a longitudinal inner tube containing heat transfer fluid (HTF). The cylinder is filled with an inhomogeneous porous medium in which porosity and permeability vary along the *z*-direction. The porosity in the direction of the *z*-axis could be increasing, constant, or decreasing. Biobased coconut oil containing the CuO nanoparticles fills the pores. Table 1 lists the properties of coconut oil, CuO, copper metal foam, and water.

**Figure 1.** The physcial model of the storage unit: (**a**) cylindrical thermal energy storage (TES) system filled with phase change material (PCM) and metal foam, along with a heat transfer fluid tube in the center; (**b**) the computational domain.



The modeling was done under the following assumptions: (1) the melted NePCM behaves as an incompressible Newtonian fluid, and the flow in the pores is laminar; (2) nanoparticles do not settle, and they are homogeneously dispersed in the PCM; (3) expansion of the NePCM is neglected; (4) the properties of the HTF in tube are considered to be constant and Newtonian. Here, a linear porosity profile is assumed as

$$
\kappa(z) = az + b,\tag{1}
$$

where Equation (1) shows the linear variation of porosity which is proposed in the present study. As seen, the porosity varies linearly in the vertical direction. The values of porosity at the lower and upper walls are

$$
\varepsilon(z=0) = \varepsilon\_0 = b,\tag{2}
$$

$$
\varepsilon(z=L) = \varepsilon\_L.\tag{3}
$$

According to the values of porosity on the lower and upper walls, and the average porosity of the domain, i.e., *εavg*, the *b* parameter is calculated as follows:

$$b = \varepsilon\_{a\overline{v}\xi} - \frac{aL}{2}.\tag{4}$$

The equations for the transient melting process include mass conservation, momentum balance, and energy conservation equations [20–22].

Mass conservation:

$$\nabla \cdot \overrightarrow{\boldsymbol{U}} = 0;\tag{5}$$

Momentum equations:

$$\begin{aligned} \frac{1}{\varepsilon} \frac{\partial \left(\rho\_{\text{lnpcm}} \overrightarrow{\boldsymbol{U}}\right)}{\partial t} + \frac{1}{\varepsilon} \rho\_{\text{lnpcm}} \left(\overrightarrow{\boldsymbol{U}} \cdot \nabla\right) \overrightarrow{\boldsymbol{U}} \frac{1}{\varepsilon} + \nabla \cdot \boldsymbol{p} \boldsymbol{I} &= \nabla \cdot \left(\mu\_{\text{lnpcm}} \frac{1}{\varepsilon} \nabla \overrightarrow{\boldsymbol{U}}\right) \\ - \rho\_{\text{lnpcm}} \rho\_{\text{lnpcm}} \overrightarrow{\boldsymbol{g}} \left(\theta\_{\text{melting}} - \theta\right) + \left(M(\theta) - \frac{\mu\_{\text{lnpcm}}}{\sigma}\right) \overrightarrow{\boldsymbol{U}} \end{aligned} \tag{6}$$

where <sup>→</sup> *U*, *p*, and *θ* are the velocity vector, pressure, and temperature, respectively. Here, *ε*, *β*, *μ*, and *g* are the porosity, volume expansion coefficient, viscosity, and gravity acceleration, respectively. The subscript *lnpcm* denotes the molten PCM, and *θmelting* is the phase change temperature. The linear Boussinesq model is employed to apply the buoyancy force, and the Forchheimer term was neglected since the natural convection velocities are very small. The source term *M*(*θ*) is a function of temperature (*θ*) and is zero in liquid PCM, but it rises to large values in the solid PCM. The large values of *M*(*θ*) induce significant resistance forces to the fluid motion and force the velocities to zero in a solid PCM. The source term *M*(*θ*) is defined as

$$\begin{aligned} M(\theta) &= B\_{\text{vc}} \frac{\xi^2(\theta) - 2\xi(\theta) + 1}{\xi^3(\theta) + 10^{-3}} \\ \xi(\theta) &= \begin{cases} 0 & \theta < \theta\_{\text{melting}} - \theta\_{\text{unlow}}/2 \\ 0.5 - \frac{\theta\_{\text{unlow}} - \theta}{\theta\_{\text{unlow}}} & \theta\_{\text{melting}} - 0.5\theta\_{\text{unlow}} < \theta < \theta\_{\text{unlow}} + 0.5\theta\_{\text{unlow}} \\ 1 & \theta > \theta\_{\text{multiting}} + 0.5\theta\_{\text{unlow}} \end{cases} \end{aligned} \tag{7}$$

where *Bvc* is a large value, which intensifies the magnitude of the temperature-dependent term. The porous permeability, *σ*, as a function of porosity, is calculated as follows [23,24]:

$$\begin{array}{l} \sigma = \frac{0.73 \times 10^{-3}}{(1 - \varepsilon)^{0.224} \frac{R\_m^{-2} \left(R\_z R\_m^{-1}\right)^{1.11}}{\pi}}\\ R\_z R\_m^{-1} = \sqrt{\frac{1.3924(1 - \varepsilon)}{3\pi}} \frac{1}{1 - \varepsilon^{25(\varepsilon - 1)}}\\ R\_m = 0.254 \times 10^{-7} \chi^{-1} \text{ (PPI)} \end{array} \tag{8}$$

where *Rz* and *Rm* are the pore characteristics, as introduced in [23,24], and *χ* is the pore density in pore size per inch (*PPI*).

Energy conservation [25]:

$$\begin{array}{c} \left[\xi(\theta)\left[\left(\rho\varepsilon\_{\mathcal{P}}\right)\_{\textit{eff\\_Inpcm}} - \left(\rho\varepsilon\_{\mathcal{P}}\right)\_{\textit{eff\\_supcun}}\right] + \left(\rho\varepsilon\_{\mathcal{P}}\right)\_{\textit{eff\\_supcun}}\right] \overset{\scriptstyle{\scriptstyle\Theta}}{\delta t} + \left(\rho\varepsilon\_{\mathcal{P}}\right)\_{\textit{Inpcm}} \overset{\scriptstyle{\scriptstyle\Theta}}{\mathcal{U}} \cdot \nabla\theta = \\ \qquad \qquad \qquad \qquad \nabla. \left(\kappa\_{\mathcal{eff},\textit{nepc}}\nabla\theta\right) + \left(V\mathcal{F}\_{\textit{na}} - 1\right)\rho\_{\textit{ppc}m}l\_{\textit{ppc}m}\varepsilon^{\underline{\bf{\mathcal{E}}}\underline{\mathcal{E}}(\theta)} \end{array} \tag{9}$$

where

$$\left(\rho c\_p\right)\_{\varepsilon f; \operatorname{lnpccm}(snpcm)} = \left(1-\varepsilon\right)\left(\rho c\_p\right)\_{\operatorname{sun}} + \varepsilon \left(\rho c\_p\right)\_{\operatorname{lnpccm}(snpcm)}.\tag{10}$$

The subscript *snepcm* denotes the characteristics of the NePCM in solid state, and *sm* refers to the solid matrix of metal foam. The term *ξ*(*θ*) denotes the liquid volume fraction and can be changed in each element. The variation of the liquid fraction represents the variation of the stored/released latent heat. The thermal conductivity of the porous medium saturated with the NePCM is considered as the weighted average of the two media and is computed as

$$
\kappa\_{eff, \text{perp}} = \mathfrak{z}(\theta) \kappa\_{eff, \text{perp}} + (1 - \mathfrak{z}(\theta)) \mathfrak{z}\_{eff, \text{perp}} \tag{11}
$$

in which,

$$
\kappa\_{eff,Impcm(supm)} = \frac{\Xi\_1}{\Xi\_1},
\tag{12}
$$

$$\begin{array}{c} \boldsymbol{\Sigma}\_{1} = \begin{bmatrix} \boldsymbol{\kappa}\_{\mathrm{lnpcm}(\mathrm{snpcm})} + \left( \left( \frac{\boldsymbol{\pi} - \boldsymbol{\pi}\boldsymbol{\varepsilon}}{3} \right)^{0.5} - \frac{1 - \boldsymbol{\varepsilon}}{3} \right) \left( \boldsymbol{\kappa}\_{\mathrm{s}\boldsymbol{m}} - \boldsymbol{\kappa}\_{\mathrm{lnpcm}(\mathrm{snpcm})} \right) \\ \times \begin{bmatrix} \boldsymbol{\kappa}\_{\mathrm{lnpcm}(\mathrm{snpcm})} + \left( \frac{1 - \boldsymbol{\varepsilon}}{3} \right) \left( \boldsymbol{\kappa}\_{\mathrm{nf}} - \boldsymbol{\kappa}\_{\mathrm{lnpcm}(\mathrm{snpcm})} \right) \end{bmatrix} \end{array} \tag{13}$$

$$\Xi\_2 = \kappa\_{\text{lppcm}(\text{suppm})} + \left[ \frac{4}{3} \left( \frac{1-\varepsilon}{3\pi} \right)^{0.5} (1-\varepsilon) + \left( \frac{\pi-\pi\varepsilon}{3} \right)^{0.5} - (1-\varepsilon) \right] \left( \kappa\_{\text{mf}} - \kappa\_{\text{lppcm}(\text{suppm})} \right)\_{\text{'}} \tag{14}$$

The governing equations for flow and heat transfer of HTF passing the tube are laminar convection heat transfer.

Mass conservation:

$$\nabla \cdot \overrightarrow{\boldsymbol{U}} = 0;\tag{15}$$

Momentum Balance:

$$
\rho\_{HTF} \frac{\partial \overrightarrow{\boldsymbol{U}}}{\partial t} + \rho\_{HTF} \left(\overrightarrow{\boldsymbol{U}} \cdot \nabla\right) \overrightarrow{\boldsymbol{U}} + \nabla \cdot \boldsymbol{p} \boldsymbol{I} = \nabla \cdot \left(\mu\_{HTF} \overrightarrow{\boldsymbol{\nabla}} \overrightarrow{\boldsymbol{U}}\right);\tag{16}$$

Energy conservation:

$$\left(\rho c\_p\right)\_{HTF} \frac{\partial \theta}{\partial t} + \vec{l} \vec{l} . \nabla \left(\left(\rho c\_p\right)\_{HTF} \theta\right) = \nabla. \left(\kappa\_{HTF} \nabla \theta\right)\_t \tag{17}$$

where subscript HTF shows the heat transfer fluid.

The thermophysical properties of the nano-enhanced phase change material can be defined as weighted functions of the pure PCM and nanoparticle properties. These properties are tabulated in Table 2. The initial and boundary conditions are listed in Table 3.

**Table 2.** Properties of the nano-enhanced phase change material.


**Table 3.** Initial and boundary conditions.


Sensible energy, latent energy, and total energy stored in the unit are

$$ES(t) = \int \left( \underbrace{\underbrace{(1 - VF\_{\text{na}}) \rho\_{pp\text{r}m} \varepsilon(z) l\_{pp\text{r}m}}\_{\text{Latent energy}} + \underbrace{(\rho \mathbb{C}\_p)\_{eff\text{\\_ref}, \text{up}\text{cm}} (\theta - \theta\_{\text{initial}})}\_{\text{Catal\text{\\_energy}} \text{ stored}} \right) dV,\tag{29}$$

where *V* is the total volume of the NePCM domain. The melt volume fraction is

$$MVF(t) = \frac{S\_{\text{In }cpcm}}{S\_{\text{In }cpcm} + S\_{\text{sucpcum}}}.\tag{30}$$

The charging power of the TES unit is defined as

$$CP = \frac{ES|\_{MVF=1}}{t|\_{MVF=1}}.\tag{31}$$

Moreover, the uniformity of the temperature field can be measured using

$$T\mathcal{U}^2(t) = \left(\int\_V \left(\theta - \theta\_{\text{mean}}\right)^2 dV\right) \Big/ \left(\int\_V dV\right),\tag{32}$$

where *θmean* is the mean temperature of the NePCM domain.

#### **3. Numerical Approach and Grid Dependency**

The model equations described above were solved by employing the finite element approach. A mesh independence test was performed to determine a balance between the precision and performance of the numerical simulation. Structural meshes with rectangular elements were employed for both computational domains. Initially, a mesh with 75 × 75 elements was considered. Then, five levels of finer meshes were used to test the solution accuracy and grid optimization. The simulation results for different grids are depicted in Figure 2. Figure 2a,b show, respectively, the mesh size's influence on the melting volume fraction (*MVF*) and local temperature at a specified point. It is seen that the temperature field is more sensitive to the mesh size. On the basis of Figure 2b, a mesh size with 125 × 125 elements was selected to satisfy the balance between the solution accuracy and simulation efficiency.

**Figure 2.** The variation of (**a**) melting volume fraction (*MVF*), and (**b**) the temperature at selected points with coordinates (5*Ri*, 0.5*L*) for various cases of grid size for *εavg* = 0.84, *VFna* = 0.04, *a* = 0.6, and *PPI* = 30.

To guarantee the reliability of the implemented model and the employed numerical solution, we studied the liquid fraction field of biobased coconut oil contained in a uniform porous rectangular enclosure by conducting the two-dimensional (2D) numerical simulation. The predicted liquid fraction field was compared with the experimental observations of Al-Jethelah et al. [18]. In this validation test, a net heat flux was imposed on the cubic cavity's left side wall, and the insulated boundary conditions were set at the other surfaces. The porosity of copper metal foam used in the container was 0.92, and the corresponding permeability was 3.3142 × <sup>10</sup>−<sup>7</sup> <sup>m</sup>2. The liquid fraction field obtained from the present model is in good agreement with the experimental field observation of Al-Jethelah et al. [18], as shown in Figure 3.

**Figure 3.** A comparison of melting fields simulated in the present simulations and the experimental study conducted by Al-Jethelah et al. [18].

To verify the computational model for application to an inhomogeneous porous domain, a square cavity occupied with an inhomogeneous metal foam utilized in [26] was studied. Here, the insulation boundary conditions were assigned to the cavity's upper and lower surfaces. However, the left and right vertical side boundaries were kept at high and low temperatures. The isotherms and streamlines of the present study are in excellent agreement with those of Xiong et al. [26], as illustrated in Figure 4. In this analysis, porosity was a linear function of the *y*-coordinate such that *εavg*=0.7. In addition, the nondimensional parameters were as follows: *Ra* = 106 and *Da* = 10<sup>−</sup>1.

**Figure 4.** Comparisons of the predicted streamlines and isotherms with the simulations of Xiong et al. [26].

#### **4. Results and Discussion**

Here, the Taguchi optimization method was adopted to systematically investigate the impact of design parameters on the melting rate of PCM in the TES unit. The volume fraction of nanoparticles (*VFna*), average porosity of porous medium (*εavg*), intensity of porosity gradient (*a*), and pore size per inch (PPI) were selected as design parameters. These parameters are known as control factors in the context of the Taguchi method. The aim was to minimize the full melting time (charging time). Thus, the charging time was selected as the target parameter. Since the time should be minimized, the target function was selected as "the lower, the better" for the Taguchi method. Following the Taguchi method, the control factors were divided into several levels. Here, five levels were selected for each factor. The details of the levels and factors are summarized in Table 4. The porosity gradient was selected in the range of −0.6 to +0.6, which covered a negative and positive gradient distribution. The zero gradient case represents uniform porosity.

**Table 4.** The range and levels of control factors.


The total combination of four factors and five levels resulted in 5<sup>4</sup> combination variables. A full melting process should be computed for each combination, which is computationally impractical. The Taguchi method uses an orthogonal table to probe the possible solution space efficiently. Here, we adopted the L25 design. The L25 design selects only 25

unique combinations of possible designs out of 5<sup>4</sup> possibilities. These selected combinations are summarized in Table 5. The simulations were carried out for all 25 cases listed in Table 5, and the corresponding values for the time of full melting (*t*|*MVF*=<sup>1</sup> ), stored energy (*ES*), and TES unit power (CP) at full charge are reported in Table 5. The Taguchi method was then used to compute the signal to noise (S/N) ratio values. The S/N ratio indicates the robustness of a factor level to possible noises. Thus, a factor with the highest S/N value is a promising candidate to produce an optimum design (lowest charging time). Using the values of *t*|*MVF*=<sup>1</sup> from Table 5, the S/N of the Taguchi method was computed, and the results are summarized in Table 5.

**Table 5.** Taguchi L25 orthogonal table corresponding to range and levels of control parameters. PPI, pore size per inch; ES, stored energy; CP, unit power.


Following the standard Taguchi method, Table 6 shows the S/N and rank values of each level and control factor. Figure 5 shows a graphical representation of S/N values in Table 6. The highest value of S/N for each control parameter shows the promising level for the minimum charging time. Table 6 also shows the rank of each factor with respect to the variation of melting time. As seen, the average porosity (*εavg*) was the most effective parameter influencing the melting time, followed by the volume fraction of nanoparticles. The porosity gradient (*a*) and pore size (*PPI*) were ranked third and fourth. Here, δ indicates the maximum difference between the computed S/N values of each factor. A higher *δ* value denotes a more influential factor.


**Table 6.** The S/N and rank values of the control factors.

**Figure 5.** Mean values of the signal-to-noise (S/N) ratios for all the levels of the controlling parameters. Optimum case: *A* = 5, *B* = 1, *C* = 4, and *D* = 4.

Figure 5 shows that Levels 4, 1, 4, and 4 corresponded to the highest S/N ratios for nanoparticle volume fraction (*VFna*), average porosity (*εavg*), porosity gradient (*a*), and pore size (PPI), respectively. The approach for computing the S/N ratio was "smaller is better"; thus, a case with a lower melting time produced a larger S/N ratio value. Using Table 4, these levels can be read as *VFna* = 0.04, *εavg* = 0.8, *a* = 0.3, and PPI = 25. The Taguchi method estimated a melting time of 8784 s for this design parameter combination. A simulation was performed for this specific combination to confirm the reduction in the TES unit's charging time. The simulation outcomes showed a full melting (charging) time of 9097 s, which is slightly higher than the estimated time of 8784 s estimated using the Taguchi method. However, a comparison between the computed melting time and the 25 cases of Table 5 confirms that this melting time was the smallest. Thus, the proposed design according to the Taguchi method was selected as the optimum design. The details of the optimum design are shown in Table 7.

**Table 7.** The optimum values of the controlling parameters.


Table 5 also shows that case 5 resulted in the highest melting time (15,900 s). Considering the investigated range of control parameters, the optimum melting time of 9097s was 74% better than the highest melting time. It should be noted that case 5 contained no nanoparticles and had the maximum average porosity of 0.9. Moreover, case 21 with a full melting time of 9150 s provided a melting time that was close to the optimum case. The only difference between these two cases was the porosity gradient.

The results of Table 5 were also used to develop a linear relationship for design parameters and the charging time of the TES. That is,

$$\text{Time} + 8892 - 243.0 \text{ V} \\ \text{F}\_{\text{nat}} + 1047.0 \text{ } \varepsilon\_{\text{avg}} \\ - 9.0 \text{ a} + 45.0 \text{ PPI.} \tag{33}$$

Here, 16 more cases were adopted to further explore the impact of the variation of the design parameters on the optimum design. The details of selected cases are presented in Table 8. Figure 6a shows the time evaluations of MVF during the charging process for various values of nanoparticle volume fraction. It is seen that the increase in the concentration of nanoparticles slightly increased the MVF. The differences were only visible in the final stages of charging.

**Table 8.** Table of 16 cases for further analysis around the optimum design.


Figure 6b depicts the predicted temperature of point A in the enclosure. The temperature in the solid region increased sharply until it reached the fusion temperature. Then, the temperature slightly increased until the entire PCM changed to liquid at this point. The temperature then increased again. The increase in nanoparticle concentration slightly increased the temperature of point A in the liquid region due to the improved thermal conductivity of NePCM and better heat transfer between the hot liquid in the tube and the PCM inside the enclosure. The melting interfaces during the melting process were very close; thus, they were not plotted.

**Figure 6.** The time variation of (**a**) *MVF*, and (**b**) temperature distribution for various nano-additive volume fractions for *Ro* = 0.1 m, *Ri* = 0.01 m, *L* = 0.3 m, *εavg* = 0.800, *a* = 0.3, and *PPI* = 25.

Figure 7 shows the phase change interfaces (- = 0.5) for various average porosity (*εavg*) values and different melting times. It is seen that the melting interface advanced in the enclosure toward the right wall with time.

**Figure 7.** Melting interface at four steps of melting progress for various average porosities for *VFna* = 0.04, *a* = 0.3, and *PPI* = 25.

Interestingly, the melting started at the bottom since the HTF fluid at the entry was hot, and the convection heat transfer coefficients at the undeveloped region of the HTF tube was also larger. Thus, the heat transfer between HTF liquid and NePCM at the bottom was more intense. The advancement of the melting interface accelerated as the porosity decreased (mass of metal foam increased). This was due to the increase in thermal conductivity of metal foam, which enhanced the heat transfer rate. It should be noted that the inlet temperature of HTF fluid was constant; hence, a higher composite thermal conductivity at the PCM side increased the heat transfer rate. Thus, the amount of stored heat in an enclosure increased for a higher composite thermal conductivity. Therefore, as seen in Figure 7 at *t* = 7800 s, there was a clear difference between the melting interfaces. Figure 8a illustrates the time history of MVF for various average porosities. The average porosity induced notable changes in the values of MVF. In agreement with the melting interfaces of Figure 7, this figure shows that the increase in average porosity raised the MVF.

**Figure 8.** The variation of (**a**) *MVF*, and (**b**) temperature uniformity for the various average porosities as a function of time for *VFna* = 0.04, *a* = 0.3, and *PPI* = 25.

Figure 8b depicts the temperature uniformity inside the PCM domain. The enclosure temperature distribution experienced a slightly higher temperature gradient during the melting process when average porosity was maximum. The rise in temperature was due to the lower composite thermal conductivity and large temperature gradients. At the early stages of heat transfer, the NePCM was in a solid state, and there was no practical molten region or natural convection circulations. As the MVF increased and the molten region grew, natural convection flows occurred. During the phase change, the temperature uniformity stayed almost constant. When the molten region grew significantly, the natural convection circulation started, and the temperature gradients also increased.

Figure 9 is plotted to show the impact of pore sizes on MVF and temperature evaluations of point A. It is seen that the variation of pore size induced a negligible impact on the MVF. There was a slight temperature variation for point A around 6000 s. Half of the enclosure was in a molten state, and there was a dominant natural convection heat transfer flow. In a microscopic view, the overall heat, which diffused to the composite PCM (foam and PCM), first channeled through the pore walls, and then the PCM inside the pores could absorb it. A decrease in the pore size could provide a larger contact surface and a better melting rate. However, since the amount of foam was low (80% or higher porosity), the heat transfer rate was limited by the foam's capacity to conduct heat into the solid PCM region. Since the amount of foam was limited (constant porosity), the condition heat transfer through the foam was also limited. Thus, the change in pore size could impact the heat transfer, but it did not play a dominant role. Moreover, a smaller pore size led to a smaller permeability, which reduced the natural convection effects.

The variation of pore sizes could affect the permeability of the medium, according to Equation (4); thus, the impact of pore size on the flow and heat transfer can be boosted in convective dominant heat transfer regimes. Since the influence of the pore size (PPI) on the MVF was minimal, the melting interfaces were not plotted for the sake of brevity.

Figure 10 display the dependency of the porosity gradient on the evaluation of melting interfaces at various time steps. Interestingly, the increase in *a* shifted the melting interface toward the left wall and resulted in an increased melting rate. According to Equation (1), the increase in *a* meant more metal foam at the bottom and less metal foam at the top. Thus, the PCM's thermal conductivity at the bottom of the enclosure was higher than that at the top. At the bottom and in the initial stages of phase change, the dominant mechanism of heat transfer was conduction. Thus, the increase in thermal conductivity enhanced the heat transfer. As the phase change continued, the molten region grew, and convection heat transfer occurred in the enclosure's top regions. In natural convection flows, the increase in porosity of metal foams increased the medium's permeability. Consequently, the liquid PCM could circulate more easily in the top area, contributing to convection heat transfer. The relative difference in full charging of TES for *a* = 0.3 (optimum case) and *a* = 0.6 (case 12, Table 8) was only 0.5%, and Figure 10 shows the close competition of these two cases at 7800 s. Case *a* = 0.3, which was the adopted optimum design, led to a better melting at the final stage of thermal energy storage (above 90% MVF) compared to case 12 in Table 8.

**Figure 9.** The variation of (**a**) *MVF*, and (**b**) temperature distribution at the various cases of pore per inch of the metal matrix as a function of time for *VFna* = 0.04, *a* = 0.3, and *εavg* = 0.800.

**Figure 10.** Melting interface at four steps of melting for various porosity gradients for *VFna* = 0.04, *PPI* = 25, and *εavg* = 0.800.

Figure 11 illustrates the time history of MVF and temperature uniformity for various values of porosity gradients. Figure 11a shows that an increase in *a* improved the MVF during the middle stages of phase change. This finding is in agreement with the advancement of the melting interfaces, as observed in Figure 10. However, in the final stages of charging, the case *a* = 0.3 provided better performance than the case of *a* = 0.6. The temperature uniformity for the initial stages of the melting process was almost independent of the *a* parameter. However, in the middle stages, where the molten region grew, the increase in *a* parameter promoted the temperature nonuniformity. In the middle stages, the increase in *a* promoted temperature nonuniformities. A high value of *a* led to a stronger heat transfer flow and natural convection flows.

**Figure 11.** The variation of (**a**) *MVF*, and the (**b**) temperature distribution at various gradients of porosity as a function of time for *VFna* = 0.04, *PPI* = 25, and *εavg* = 0.800.

Nine points in the NePCM–metal foam region were selected to investigate the temperature variations during the melting process. These points were located at a1(2*Ri*, 0.1*L*/3), a2(2*Ri*, *L*/2), a3(2*Ri*, 2.9*L*/3), b1(5*Ri*, 0.1*L*/3), b2(5*Ri*, *L*/2), b3(5*Ri*, 2.9*L*/3), c1(9*Ri*, 0.1*L*/3), c2(9*Ri*, *L*/2), and c3(9*Ri*, 2.9*L*/3). Figure 12 shows the variation of temperatures at these nine points in the PCM domain during the melting process. All points started with an initial temperature of 293 K followed by a linear increase in temperature. The linear segment showed pure conduction in a super-cold solid PCM before it reached a melting temperature. Then, there was a flattened segment, in which the temperature slightly increased. This was the phase change stage, where the PCM absorbed the energy in the form of latent heat. Then, there was a semi-linear increase in the temperature, which corresponded to the heat transfer in the molten region with no phase change. The time variations of temperatures for all cases showed that the impact of the porosity gradient (*a*) on temperature profiles was minimal in the solid state and during melting phase change. However, after the phase change, the molten NePCM circulated in the enclosure, enhancing the effect of the porosity gradient.

Here, *a*1, *b*1, and *c*<sup>1</sup> were the points placed at the bottom of the enclosure. Point a1 was next to the tube wall, point b1 was in the middle, and c1 was next to the insulated shell wall. Since a1 was the closest point next to the tube, it experienced a sharp temperature rise at the beginning of the charging process. In contrast, c1 showed a smooth temperature variation since it was far from the tube. Point c3 was placed near the top-right corner of the enclosure, which was the last place to be melted. Thus, during the melting process, this point remained at the fusion temperature.

Figure 13 compares the full charging power of TES units for various values of design factors. Figure 13a shows that the variation of nanoparticle concentration slightly changed the unit power. This is because the presence of nanoparticles improved the heat transfer rate; they also reduced the heat capacity of the storage unit. Figure 13b depicts an increasing trend of power (*CP*) as the porosity gradient (*a*) increased. Both cases of *a* = 0.3 and a = 0.6 showed almost similar charging powers. The increase in average porosity (Figure 13c) reduced the charging power. As the average porosity increased, the heat capacity of the TES also increased. However, a reduction in the composite PCM's effective thermal conductivity weakened the heat transfer rate and consequently increased the charging time. The variation of pore sizes had a negligible impact on the TES power since PPI variation did not induce notable changes in either the unit's heat capacity or the melting time.

**Figure 12.** Variation of temperature at selected points during the charging process for optimum condition.

**Figure 13.** Variation of charging power during the charging process for optimum designs as a reference and the variation of (**a**) the volume fraction of the nano-additives, (**b**) the gradient of porosity for optimum conditions, (**c**) the average porosity, and (**d**) the number of pores per inch of the metal matrix.
