Latent Heat Thermal Storage in Non-Uniform Metal Foam Filled with Nano-Enhanced Phase Change Material

Abstract

The melting heat transfer of CuO—coconut oil embedded in a non-uniform copper metal foam—was addressed. Copper foam is placed in a channel-shaped Thermal Energy Storage (TES) unit heated from one side. The foam is non-uniform with a linear porosity gradient in a direction perpendicular to the heated surface. The finite element method was applied to simulate natural convection flow and phase change heat transfer in the TES unit. The results showed that the porosity gradient could significantly boost the melting rate and stored energy rate in the TES unit. The best non-uniform porosity corresponds to a case in which the maximum porosity is next to a heated surface. The variation of the unit placement’s inclination angle is only important in the final stage of charging, where there is a dominant natural convection flow. The variation of porous pore size induces minimal impact on the phase change rate, except in the case of a large pore size of 30 pore density (PPI). The presence of nanoparticles could increase or decrease the charging time. However, using a 4% volume fraction of nanoparticles could mainly reduce the charging time.

1. Introduction

The rise in the participation of renewable energy in the prevailing cooling, heating, and power infrastructures demands new encounters to be overcome. Thermal energy storage systems are capable of storing/releasing a notable amount of heat in a compact volume. Thermal energy storage (TES) systems could also help shift heat loads or damp transient heat load fluctuations. Thus, such systems are crucial to heat load distribution time and reduction in the design size of thermal systems. For instance, the presence of a TES unit could act as a thermal buffer and allow a design of a heating system for average thermal load instead of maximum thermal load.

Latent Heat Storage (LHS) is an advantageous type of thermal energy storage due to its large melting enthalpy and the sharpness of the working temperature. From a thermodynamic point of view, a pure substance can phase change at a constant fusion temperature and release/absorb a significant amount of heat. During the last two decades, latent heat storage using Phase Change Materials (PCMs) has been employed in many engineering applications [1]. Besides the high thermal capacity features of LHS, practical application of PCMs encounters some heat transfer drawbacks, which are mainly due to the low thermal conductivity of PCMs [2,3]. Low thermal conductivity acts as a barrier to heat transfer and drastically increases the charging/discharging times of TES units. A fast charging/discharging TES is essential for many industrial applications so that energy storage units can quickly supply transient demanded thermal loads. Thus, improving the heat transfer performance of TES units is a hot topic. An increase in heat transfer performance could be achieved by increasing the thermal conductivity enhancement (TCE) of PCMs or the effective thermal conductivity of the overall TES.

Different techniques have been used to enhance the conductivity of PCMs. Among these techniques, we can quote the embedding of fins or extended surfaces [4,5,6], encapsulated PCMs [7,8], insertion of particles/nanoparticles of high thermal conductivity within the PCM [9,10,11], and inclusion of high thermal conductivity material structures [12,13,14,15]. One promising TCE technique consists of the use of porous open-cell metal foam [16,17,18]. Metal foams have high thermal conductivity, large porosity, and large surface area per unit of volume, which makes them potentially good candidates for promoting latent heat storage performance. Additionally, they are characterized by a strong mechanical structure and stable thermo-physics properties.

Huang et al. [19] conducted an experimental protocol to make a composite of PCM/ metal foam using vacuum melting infiltration. Myristyl alcohol (MA) was used as PCM and copper/nickel as foams. Measurements showed that thermal conductivity was improved by 1.8 with nickel foam (40 PPI) and 7.51 with copper foam (40 PPI). Martenelli et al. [20] experimented with the phase change heat transfer of a PCM in a copper-foam-filled shell-and-tube heat exchanger. These researchers found that the charging and discharging times were reduced remarkably.

Xiao et al. [21] investigated the effect of porosity on the effective thermal conductivity of a paraffin/nickel foam composite and paraffin/copper foam composite experimentally and found that thermal conductivity was drastically improved. Moreover, they found that pore size induces minimal impact on phase change performance.

Some researchers have tried to integrate two or more TCE techniques to further improve the heat transfer performance of TES units. For instance, a combination of fins and Nano-enhanced PCMs (NePCM) [22,23], a combination of Nano-Encapsulated PCMs (NEPCMs) and fins [24], and a combination of NePCM and metal foams have been employed.

Mahdi et al. [25] studied the discharge of LHS numerically during the solidification of multiple PCMs with nanoparticles in a shell-and-tube containment with a cascade of multiple foams. They showed that with the application of multiple cascade foams and multiple PCMs, the solidification time can be enhanced by 17 times compared with the reference case of one PCM with no foam. An experiment was performed by Al Jethelah et al. [26] to examine the charging of latent heat during the melting of bio-based NePCM, embedded with aluminum foam. The NePCM was made of CuO nanoparticles dispersed in coconut oil. The authors examined the melting rate of NePCM in aluminum foams with different porosities of 88%, 92%, and 96%. It was found that the insertion of metal foam boosts the melting rate of PCM significantly. Thus, the performance of the storage unit after employing nanoadditives and metal foam was much better than only using nanoadditives. These researchers also reported that a dense metal foam with a porosity of 88% suppresses natural convection flow and leads to uniform melting, but natural convection plays a clear role when porosity increases to 96%.

Li et al. [27] analyzed the effectiveness of a triple-tube heat exchanger by filling the middle tube with NePCM/copper foam composite. Li et al. tested different nanoparticle concentrations and metal foam porosities. The outcomes showed that the melting/solidification time was shortened by 25.9/28.2% by adding 5% copper nanoparticles. Additionally, the melting/solidification time was reduced by 83.7/88.2% due to the presence of metal foam with 95% porosity. This reduction in charging/discharging time confirms the advantage of adding metal foam compared with adding nanoparticles.

Utilizing metal foam with non-uniform porosity or employing a cascade of metal foams with different porosities is a novel idea that could adjust local thermal conductivity to the initial melting times and allow degrees of natural convection in the final stages. There is a scarce number of publications on the phase change of PCMs in non-uniform metal foams. One of the main reasons could be the difficulties in producing non-uniform metal foams. Recently, with the advancement of 3D metal printers and other production methods, synthesis of non-uniform metal foams is feasible. Regarding the phase change of non-uniform metal foams, Yang et al. [28] proposed a model of the vertical porosity gradient of the metal foam, in which the porosity increases linearly from bottom to top. Such non-uniform porosity could shorten the full melting time (charging time) by enhancing natural convection and improving heat transfer performance, compared with a case of uniform porosity. Mahdi et al. [29] investigated the melting and solidification of PCM-filled metal foam in a shell-and-tube thermal storage system theoretically. The space between the two concentric cylinders was filled with multiple-segments of different porosities. It was found that the performance of the storage system in terms of charging and discharging times was improved by cascading the porosity. Recently, Girish et al. [30] investigated the phase change heat transfer in a cylindrical heatsink filled with composite PCM-open metal foam. They used a vertical porosity gradient from the bottom to the top of the cylinder with uniform PPI density. They showed an increase in performance by using non-uniform porosity metal foam.

2. Mathematical Model

2.1. Physical Model

A schematic view of an inclined porous rectangular enclosure saturated by a nano-enhanced phase change material, studied in this work, is depicted in Figure 1. The active wall of the enclosed medium, which is along the y-axis and located at x = 0, is kept at a high temperature, i.e. T_h. However, the other walls are fully insulated. The PCM and the nanoadditives on the NePCM are coconut oil and CuO nanoparticles, with the thermophysical characteristics listed in

Table 1. Thermophysical properties of the coconut oil and nanoadditives [26,31].

Properties	Coconut Oil (Measured)		Metal Foam	CuO Nano Additives
	Liquid (32 C)	Solid (15 C)
Cp (JkgK−1)	2010	3750	386	540
${\rho }_{PCM}$(kgm−1)	914	920	8900	6500
k (Wm−1K−1)	0.166	0.228	380	18
${\mu }_{PCM}$(Nsm−2)	0.0326	-	-	-
hf (kJkg−1)	103	-	-	-
Tme (ΔTme)	-	24 C	-	-
Pr	394.73	-	-	-

.

The length of the domian, i.e. L, is 40 mm, and the width, i.e. W, is 20 mm. The porosity varies as:

$${\epsilon }_{avg}=ax+b$$

(1)

${\epsilon }_{avg}$ is the average porosiy. a as the porosity gradient is the input changeable parameter, and the b parameter is defined based on the values of ${\epsilon }_{avg}$ and a.

2.2. Convective Phase Change Heat Transfer in NePCM

NePCM melts over time, and the natural convection of molten NePCM occurs in the liquid phase. The transient phase change process with natural convection in the porous domain is described using the volume-averaged approach:

Mass conservation

$$\nabla \cdot \overrightarrow{\mathit{q}}=0$$

(2a)

Momentum equations

$$\begin{array}{l}\frac{1}{\epsilon }{\rho }_{npcm,l}\frac{\partial \overrightarrow{\mathit{q}}}{\partial t}+\frac{1}{\epsilon }{\rho }_{npcm,l}\left(\overrightarrow{\mathit{q}}\cdot \nabla \right)\overrightarrow{\mathit{q}}\frac{1}{\epsilon }=\nabla \cdot \left[-p\mathbf{I}+{\mu }_{npcm,l}\frac{1}{\epsilon }\left(\nabla \overrightarrow{\mathit{q}}+{\left(\nabla \overrightarrow{\mathit{q}}\right)}^T\right)\right]-\frac{{\mu }_{npcm,l}}{\kappa }\overrightarrow{\mathit{q}}\\ + {\left(\rho \beta \right)}_{npcm,l}\overrightarrow{\mathit{g}}\left(\theta -{\theta }_{mel}\right) + F\left(\theta \right)\overrightarrow{\mathit{q}}\end{array}$$

(2b)

in which,

$$\overrightarrow{\mathit{g}}=g\left(\cos \gamma \overrightarrow{i} + \sin \gamma \overrightarrow{j}\right)$$

(3a)

$$F\left(\theta \right) = A_{mush}\frac{1-2l\left(\theta \right) + l^2\left(\theta \right)}{\Gamma +l^3\left(\theta \right)};\hspace{1em}{A}_{mush} = 5\times {10}^5\hspace{1em}and\hspace{1em}\mathsf{\Gamma }= {10}^{-3}$$

(3b)

$$l\left(\theta \right) = \left\{\begin{array}{cc}0\hfill & \theta <{\theta }_{mel}-\frac{\mathsf{\Delta }{\theta }_{mel}}{2}\hfill \\ \frac{\theta - {\theta }_{mel}}{\mathsf{\Delta }{\theta }_{mel}}\hfill & + \frac{1}{2}{\theta }_{mel}-\frac{\mathsf{\Delta }{\theta }_{mel}}{2}<\theta <{\theta }_{mel}+\frac{\mathsf{\Delta }{\theta }_{mel}}{2}\hfill \\ 1\hfill & \theta >{\theta }_{mel}+\frac{\mathsf{\Delta }{\theta }_{mel}}{2}\hfill \end{array}\right.$$

(3c)

Additionally, the permeability of the porous domain can be expressed as [32]:

$$\begin{array}{l}\kappa =\frac{73 \times {10}^{-5}{D}_p^2}{{\left(1 - \epsilon \right)}^{0.224}{\left(D_l{D}_p^{-1}\right)}^{1.11}}\\ D_l{D}_p^{-1}={\left(\frac{1 - \epsilon }{3\pi }\right)}^{0.5}\frac{1.18}{1-e^{-\left(1-\epsilon \right)/0.04}}\\ D_p=254\times {10}^{-4}PP{I}^{-1} \end{array}$$

(4)

Energy conservation

$${\left(\rho c_p\right)}_{eff}\frac{\partial \theta }{\partial t} + {\left(\rho c_p\right)}_{npcm,l}\overrightarrow{\mathit{q}} · \nabla \theta = \nabla · \left({\lambda }_{eff}\nabla \theta \right) - \left(1-V{F}_{na}\right){\rho }_{pcm,l}{h}_{f,pcm}\epsilon \frac{\partial l\left(\theta \right)}{\partial t}$$

(5)

in which,

$${\left(\rho c_p\right)}_{eff} = l\left(\theta \right)\left[{\left(\rho c_p\right)}_{eff,l} - {\left(\rho c_p\right)}_{eff,s}\right] + {\left(\rho c_p\right)}_{eff,s}$$

(6a)

$${\left(\rho c_p\right)}_{eff,l\left(s\right)} = \left(1-\epsilon \right){\left(\rho c_p\right)}_{mf}+\epsilon {\left(\rho c_p\right)}_{npcm,l\left(s\right)}$$

(6b)

The subscript mf denotes the metal foam. The thermal conductivity of PCM depends on the melt fraction of NePCM in an element. Thus, the melt volume fraction is used, and the effective thermal conductivity of metal foam and NePCM is evaluated as:

$${\lambda }_{eff}=l\left(\theta \right){\lambda }_{eff,l} + \left(1-l\left(\theta \right)\right){\lambda }_{eff,s}$$

(7)

Some researchers argued for possible relations for evaluation of the effective thermal conductivity of a composite metal foam and the material inside pores, as reviewed in Ranut [33]. The below relation has been generally adopted in many literature works. Thus, it has also been selected for the current investigation [13,34]:

$$\begin{array}{l}{\lambda }_{eff,l\left(s\right)}=\frac{\left[{\lambda }_{npcm,l\left(s\right)} + \pi \left(\sqrt{\iota }-\iota \right)\mathsf{\Delta }{\lambda }_{l\left(s\right)}\right]\left[{\lambda }_{npcm,l\left(s\right)} + \left(\iota \pi \right)\mathsf{\Delta }{\lambda }_{l\left(s\right)}\right]}{{\lambda }_{npcm,l\left(s\right)} + \left[\frac{4}{3}\sqrt{\iota }\left(1 - \epsilon \right) + \pi \sqrt{\iota } - \left(1 - \epsilon \right)\right]\mathsf{\Delta }{\lambda }_{l\left(s\right)}}\\ \iota =\frac{1 - \epsilon }{3\pi }, \mathsf{\Delta }{\lambda }_{l\left(s\right)}={\lambda }_{mf}-{\lambda }_{npcm,l\left(s\right)}\end{array}$$

(8)

2.3. NePCM Thermo-Physical Properties

The density of NePCM is the linear average of the PCM and nanoparticles [35]:

$${\rho }_{npcm}={\rho }_{pcm}+V{F}_{na}\left({\rho }_{na}-{\rho }_{pcm}\right)$$

(9a)

$${\rho }_{pcm}\left(\theta \right) = l\left(\theta \right)\left({\rho }_{pcm,l}-{\rho }_{pcm,s}\right) + {\rho }_{pcm,s}$$

(9b)

where VF_na represents the nanoparticles’ volume fraction. Here, the well-known Brinkman relation is applied to approximate the dynamic viscosity of molten NePCM:

$${\mu }_{npcm,l}={\mu }_{pcm,l}{\left(1-V{F}_{na}\right)}^{-2.5}$$

(10)

The coefficients of thermal expansion and effective heat capacities are also computed as the linear average of density and the expansion coefficient. Density is also involved in this relationship, since the production of these terms includes thermal terms and buoyancy force:

$${\rho }_{npcm,l}{\beta }_{npcm,l}={\rho }_{pcm,l}{\beta }_{pcm,l}+V{F}_{na}\left({\rho }_{na}{\beta }_{na}-{\rho }_{pcm,l}{\beta }_{pcm,l}\right)$$

(11)

$${\left(\rho c_p\right)}_{npcm}={\rho }_{pcm}{c}_{p,}{}_{pcm}+V{F}_{na}\left({\rho }_{na}{c}_p{}_{,na}-{\rho }_{pcm}{c}_{p,}{}_{pcm}\right)$$

(12a)

$${\rho }_{pcm}{c}_{p,}{}_{pcm}\left(\theta \right) = l\left(\theta \right)\left[{\rho }_{pcm,l}{c}_p{}_{,pcm,l}-{\rho }_{pcm,s}{c}_p{}_{,pcm,s}\right] + {\rho }_{pcm,s}{c}_p{}_{,pcm,s}$$

(12b)

The well-known Maxwell model is applied to estimate the thermal conductivity of NePCM:

$$\frac{{\lambda }_{npcm,l\left(s\right)}}{{\lambda }_{pcm,l\left(s\right)}}=\frac{\left({\lambda }_{na}+2{\lambda }_{pcm,l\left(s\right)}\right) - 2V{F}_{na}\left({\lambda }_{pcm,l\left(s\right)}-{\lambda }_{na}\right)}{\left({\lambda }_{na}+2{\lambda }_{pcm,l\left(s\right)}\right) + V{F}_{na}\left({\lambda }_{pcm,l\left(s\right)} - {\lambda }_{na}\right)}$$

(13)

2.4. Model Conditions

Boundary condition on the hot wall:

$$\theta ={\theta }_h, \overrightarrow{\mathit{q}}=0$$

(14a)

Boundary condition on insulated walls:

$$\frac{\partial \theta }{\partial n}=0, \overrightarrow{\mathit{q}}=0$$

(14b)

Initial condition for the computational domain:

$$\theta =295.15 \mathrm{K}, \overrightarrow{\mathit{q}}=0$$

(14c)

2.5. Characteristics Parameters

The cumulative heat transferred ($ES\left(t\right)$) from the hot wall to the NePCM is exactly equal to the stored energy in the unit. $ES\left(t\right)$ can be calculated as the following:

$$ES\left(t\right) = {\displaystyle \underset{\forall }{\int }{\left(\rho C_p\right)}_{eff}\left(\theta -{\theta }_{initial}\right)} d\forall +{\displaystyle \underset{\forall }{\int }\left(1-V{F}_{na}\right){\rho }_{pcm,l}\epsilon h_{pcm}d\forall }$$

(15)

The melting volume fraction is evaluated as:

$$MVF\left(t\right) = \frac{{\displaystyle \underset{\forall }{\int }l\left(\theta \right)d\forall }}{{\displaystyle \underset{\forall }{\int }d\forall }}$$

(16)

In the practical application of the LHS unit, the performance of the unit can be evaluated by charging power. This parameter shows the capacity of PCM to store energy and depends on the amount of the energy stored at 100% melting volume fraction:

$$P=\frac{ES \mathrm{when} \mathrm{the} \mathrm{melting} \mathrm{process} \mathrm{complete}}{\mathrm{complete} \mathrm{melting} \mathrm{time}}$$

(17)

To further discuss temperature distribution, the temperature uniformity index TU(t) was proposed, w as described below:

$$TU\left(t\right) = \sqrt{{\displaystyle \underset{\forall }{\int }{\left(\theta \left(x,t\right) - {\theta }_{ave}\right)}^{2}d\forall }}$$

(18a)

$${\theta }_{ave}=\frac{{\displaystyle \underset{\forall }{\int }\theta d\forall }}{{\displaystyle \underset{\forall }{\int }d\forall }}$$

(18b)

3. Numerical Approach and Grid Dependency

The Finite Element Method was employed to solve simultaneous dominant partial differential equations (PDEs). User-defined codes were developed to model phase change in a non-uniform porous medium. The PDEs are discretized by the use of rectangular elements in the computational region. As previously mentioned in detail, the thermophysical characteristics of NePCM are determined based on temperature, defining the liquid, mushy and solid zones, and the concentration of dispersed nanoadditives. Mesh quality is a key parameter influencing the precision of results for a numerical problem. To find an appropriate grid, five sets of grids are examined. The selected parameters in this test are L = 50 mm, VF_na = 0.06, ε_avg = 0.900, a = 4, γ = π/4, PPI = 10. As Figure 2a,b depict, no significant changes in the Melted Volume Fraction (MVF) and the temperature at the center of the enclosure could be seen when varying the mesh sizes. Based on these results, a uniform mesh of 150 × 150 intervals will be used for the present study. Additionally, the time-stepping with a maximum number of 8 iterations in each time step is automatically adopted.

The validation of the employed numerical approach can be confirmed by comparisons between the simulated results and the numerical/experimental data available in the literature [26,36]. In the first comparison, the melted liquid domains of the experimental study conducted by Al-Jethelah et al. [26] and the established numerical model are illustrated in Figure 3. Al-Jethelah et al. [26] carried out an experimental study to track the melting flow corresponding to a bio-based NePCM inserted in an aluminum foam enclosed medium with constant heat flux imposed on the left sidewall. The bio-based PCM and high thermal conductivity particles were coconut oil and CuO. The dimensions of the porous enclosed medium were 72 mm in height and 50 mm in width. In the comparison, illustrated in Figure 3, PPI = 5, ε = 92 and κ = 3.3142 × 10⁻⁷.

According to the comparative illustrations in Figure 4, the obtained melted liquid domains from the established model are compatible with the experimental observations. In the second comparative analysis, the reliability of the current adopted model is tested through the progressive melting fronts from this work and the experimental-numerical study discussed in [36]. The problem physics studied in [36] was a square enclosed medium with a size of 100 mm containing Octadecane with the fusion temperature of 30 °C. Several other researchers simulated the investigation conducted by Bertrand et al. [36] as a benchmark. As shown in Figure 4, the numerical approach agrees with numerical-experimental outcomes, demonstrating that melting flow is appropriately modeled.

4. Results and Discussion

In the present model, flow and heat transfer behavior in the enclosure are influenced by the following parameters: average porosity ($0.8\le {\epsilon }_{avg}\le 0.9$), the volume fraction of nanoadditives ($0\le V{F}_{na}\le 0.04$), the gradient of porosity, ($-4\le a\le 4$), the inclination angle of the enclosed medium ($0\le \gamma \le \pi$) and the pore per inch of the metal matrix ($10PPI\le \chi \le 30PPI$).

It is obvious that performing numerical experiments to test all possible combinations of the abovementioned parameters presents a substantial computational cost and is impracticable. In order to optimize such experiments, the L-25 orthogonal array statistical technique is employed as follows.

Each one of the five parameters, or control factors, is given five values, each representing a level (

Table 2. The range and levels of control parameters.

Factors	Control Parameters	Level 1	Level 2	Level 3	Level 4	Level 5
A	VFna	0.00	0.01	0.02	0.03	0.04
B	εavg	0.800	0.825	0.850	0.875	0.9
C	a	−4	−2	0	2	4
D	γ	0	π/4	π/2	3π/4	π
E	PPI	10	15	20	25	30

).

Testing all possible combinations would require 5⁵ experiments. However, The L-25 orthogonal array method is employed to explore the design space systematically. As the objective here is to select the design that would lead to the fastest melting,

Table 3. L-25 orthogonal table corresponding to range and levels of control parameters.

Experiment Number	Control Parameters					Time (s) When MVF = 0.9	ES (kJ/m)	P (kW/m)
	VFna	εavg	a	γ	PPI
1	0.00	0.800	−4	0	10	480	320.30	0.667
2	0.00	0.825	−2	π/4	15	488	328.86	0.674
3	0.00	0.850	0	2π/4	20	521	337.71	0.648
4	0.00	0.875	2	3π/4	25	627	347.45	0.554
5	0.00	0.900	4	π	30	773	359.82	0.466
6	0.01	0.800	−2	2π/4	25	397	320.90	0.808
7	0.01	0.825	0	3π/4	30	404	329.27	0.815
8	0.01	0.850	2	π	10	470	338.52	0.720
9	0.01	0.875	4	0	15	579	349.09	0.603
10	0.01	0.900	−4	π/4	20	4460	335.91	0.075
11	0.02	0.800	0	π	15	351	320.95	0.914
12	0.02	0.825	2	0	20	428	329.95	0.772
13	0.02	0.850	4	π/4	25	437	339.40	0.776
14	0.02	0.875	−4	2π/4	30	1288	332.20	0.258
15	0.02	0.900	−2	3π/4	10	1174	342.28	0.292
16	0.03	0.800	2	π/4	30	324	321.24	0.991
17	0.03	0.825	4	2π/4	10	358	330.39	0.922
18	0.03	0.850	−4	3π/4	15	774	325.11	0.420
19	0.03	0.875	−2	π	20	795	334.15	0.420
20	0.03	0.900	0	0	25	810	343.61	0.424
21 *	0.04	0.800	4	3π/4	20	318	321.79	1.012
22	0.04	0.825	−4	π	25	615	317.63	0.517
23	0.04	0.850	−2	0	30	606	326.23	0.538
24	0.04	0.875	0	π/4	10	830	335.39	0.404
25	0.04	0.900	2	2π/4	15	712	345.21	0.485

* The row 21 was selected as the best design.

presents the details of the L-25 orthogonal array and its corresponding parameters. It is shown that P is maximum for the following values of the control factors, which are considered optimal values:

$$V{F}_{na}= 0.04,{\epsilon }_{avg}= 0.8,a= 4,\gamma = 3\pi /4,PPI= 20$$

(19)

Further experiments are performed around the optimal values and summarized in

Table 4. A further 21 cases to explore the design space around the optimum design.

Experiment Number		Control Parameters					Time When MVF = 0.9	ES (kJ/m)	P (kW/m)
	Parameter Study	VFna	εavg	a	γ	PPI
1	VFna	0.00	0.800	4	3π/4	20	351.511	329.927	0.9386
2		0.01	0.800	4	3π/4	20	334.002	327.870	0.9816
3		0.02	0.800	4	3π/4	20	369.538	325.911	0.8819
4		0.03	0.800	4	3π/4	20	365.244	323.935	0.8869
5	εavg	0.04	0.825	4	3π/4	20	392.582	328.438	0.8366
6		0.04	0.850	4	3π/4	20	467.525	335.255	0.7171
7		0.04	0.875	4	3π/4	20	528.626	342.399	0.6477
8		0.04	0.900	4	3π/4	20	688.149	350.656	0.5096
9	a	0.04	0.800	−4	3π/4	20	469.522	312.043	0.6646
10		0.04	0.800	−2	3π/4	20	407.366	314.557	0.7722
11		0.04	0.800	0	3π/4	20	348.865	316.916	0.9084
12		0.04	0.800	2	3π/4	20	338.125	319.286	0.9443
13	γ	0.04	0.800	4	0	20	347.665	321.838	0.9257
14		0.04	0.800	4	π/4	20	315.240	321.787	1.0208
15		0.04	0.800	4	2π/4	20	353.908	334.902	0.9463
16		0.04	0.800	4	π	20	360.902	321.827	0.8917
17	PPI	0.04	0.800	4	3π/4	10	312.557	321.807	1.0296
18		0.04	0.800	4	3π/4	15	391.438	321.881	0.8223
19		0.04	0.800	4	3π/4	25	346.688	321.853	0.9284
20		0.04	0.800	4	3π/4	30	307.910	321.748	1.0449
21	Optimum case	0.04	0.800	4	3π/4	20	319.425	321.766	1.0073

in order to confirm the design optimization for these values. In each experiment, one control factor is varied among four values, while the other factors are kept constant. It can be seen that in most cases, the other experiments led to a greater melting duration compared with the one obtained with the optimal values, which confirms the optimization of the design. However, some minor exceptions exist in which the melting duration is slightly shorter than the optimal one, but the difference remains relatively low. In these exceptions, the values of the parameters VF_na, ε_avg, and a are equal to their optimal values, but those of γ and PPI are changed. This indicates that the influence of the latter parameters on optimization is relatively weak compared with the other parameters, as the power is still maximized when they are changed.

In any case, the flow and thermal patterns inside the enclosure should be studied to better understand how various control parameters affect NePCM behavior and clarify how optimal values lead to better results in melting duration and charging power.

Figure 5a,b show the variation of porosity and metal foam permeability for the optimum case, respectively. The optimum value of a was 4. Thus, as seen in Figure 5a, the porosity next to the left wall is low (due to the high amount of metal foam) and linearly increases toward the right wall. From Figure 5b, it can be seen that by moving from left to right the permeability increases. The increase in permeability is due to the increase in local porosity.

Figure 6 and Figure 7 depict the streamlines and isothermal contours in the cavity for the different nanoparticle volume fraction, VF_na. It should be noted that the streamlines indicate the presence of melted PCM, and the red line separates the liquid and solid phases. It can be seen that the behavior of NePCM is similar in all cases. The PCM starts melting initially in the left region near the heated wall. As time goes, the hot melt circulates upwards while the colder one goes downwards, and convective flow takes place, as illustrated by the single recirculation zone occupying the cavity. Nonetheless, the vertical shape of the isotherms indicates that heat transfer remains mainly dominated by conduction, expect in the final stages of melting. As the melted region grows, natural convection flows can be initiated, and isotherms would deflect due to the natural convection effects. The melting interface is more advanced toward the right when VF_na is raised, indicating that the volume of melted PCM is slightly higher when VF_na is increased. Furthermore, and following the development of the flow patterns, the isothermal contours reveal that the temperature is higher when a greater value of VF_na is employed in the same zones of the cavity. Indeed, dispersing highly conductive nanoparticles in pure PCM improves its thermal conductivity and intensifies heat transfer between the different regions, raising the overall temperature and enhancing PCM melting.

Figure 8 shows the time evolution of the melted volume fraction (MVF) and the stored energy (ES). It is clear that the MVF is at its maximum when VF_na = 0.04, which corresponds to the optimal value, and full melting (MVF = 1) is achieved faster for this value of VF_na. On the other hand, VF_na seems to have little effect on ES. However, considering the charging power, defined as the ratio of ES to full melting time, it becomes clear that the power is at its maximum when VF_na = 0.4 and confirms the optimization for that value.

The effect of the average porosity ε_avg on flow and temperature patterns in the enclosure is illustrated in Figure 9 and Figure 10. It is clear that the melting front is developed, and the size of the melted PCM zone is greater when ε_avg is decreased. Full melting of the PCM is reached faster in that case. In addition, the temperature in the cavity increases when a lower value of ε_avg is utilized. These observations are confirmed in the plots of MVF and ES in Figure 11. A substantial increase in both the MVF and ES is observed when ε_avg is reduced. Full melting time is at its minimum when ε_avg = 0.8, and ES is also at its maximum in this case, which indicates a maximization of the charging power and confirms the optimization when ε_avg = 0.8. Indeed, the porosity is defined as the ratio of the volume of voids to the total volume of the metal foam. When ε_avg is reduced, the volume of the solid matrix is increased relative to that of the fluid. Since this solid matrix is thermally conductive, heat transfer is enhanced in this case, which results in higher temperature and faster PCM melting.

The streamlines and isotherms are plotted in Figure 12 and Figure 13 for various values of the porosity gradient a. A positive value of an indicates that porosity is low near the left hot wall and increases in the right direction, while a negative value indicates a decreasing porosity starting from the left wall. Raising a clearly tends to intensify heat transfer and accelerate PCM melting. The fastest melting occurs when a = 4, which validates the optimal aspect of that value. The reason behind this behavior is that in the initial stages, heat transfer is dominated by conduction. Conductive effects are enhanced near the left wall for a higher value of a, i.e., low porosity in that area and consequent higher metal volume, which tends to melt PCM at a faster rate. Once the PCM has melted and the convective effects take place in the last stage in the right part of the cavity, the high porosity in that area facilitates liquid circulation, which further contributes to convective heat transfer. The opposite occurs when a is decreased, mainly to negative values, where the porosity near the wall is high; conductive heat transfer is inhibited, and porosity is low in the right region, which slows down the flow and diminishes the convective effects. For this reason, both MVF and ES increase when a is higher and are at their maximum when a = 4, as observed in Figure 14. The charging power is at its maximum when a = 4, and that value is therefore confirmed as optimal.

Figure 15 and Figure 16 illustrate the impact of the inclination angle of cavity γ on flow and thermal patterns. It can be seen that γ has a negligible effect on the streamlines and isotherms in the initial stages. As time goes by, a slight development in the melting front can be seen for γ = π/4 and γ = 3π/4 compared with the cases of γ = 0 and γ = π. In fact, tilting the cavity impacts the location of the hot wall and influences the effects of gravitational forces. This consequently affects the convective heat transfer, which is slightly enhanced when the cavity is tilted. This appears in the plots of MVF and ES in Figure 17, where it can be seen that both parameters have the same variation in the initial stages, but this trend changes over time, as both MVF and ES, and consequently the power P, become relatively higher for γ = π/4 and γ = 3π/4.

The flow patterns and isothermal contours are plotted in Figure 18 and Figure 19 for different values of the pores-per-inch of metal square PPI. It should be noted that even though the PPI is varied, the average porosity is constant, which means the total metal volume is the same in all the cases. When PPI is decreased, the number of pores is lower, so the thickness of the solid metal matrix is increased, which enhances thermal conductivity. However, in that case, the interface contact area between PCM and the metal structure is reduced, which may diminish heat transfer between the two components. For this reason, PPI seems to have a little effect on the thermal behavior of PCM due to the balance between positive and negative effects, as appears in the plots of MVF and ES in Figure 20. Nonetheless, it seems that the case PPI = 15 leads to lower MVF and ES compared with the other cases, which indicates that the improvement of thermal conductivity fails to balance the decrease in the interface area in this case.

5. Conclusions

The melting heat transfer of NePCMs in non-uniform metal foam was addressed systematically. Porosity changes linearly, in the x-direction. Non-linear porosity could improve local effective thermal conductivity in conduction-dominant regions by decreasing porosity, while low porosity could facilitate natural convection in convection-dominant regions. The FEM was applied to solve the governing equations. The impacts of the volume fraction of nanoparticles, average porosity of metal foam, porosity gradient, inclination angle, and pore size were investigated on the phase change process and melting rate. The results were reported in the form of isotherms, melting maps, streamlines, and graphs of melting rates and stored energy. The following main outcomes were found:

• Increasing the volume fraction of the nanoparticles improves the effective thermal conductivity of PCM, which enhances heat transfer and accelerates PCM melting. The value that leads to the fastest melting and maximum charging power P is VF_na = 0.04.
• Reducing the average porosity ε_avg of metal foam increased the volume of the metallic matrix and resulted in higher thermal conductivity. This led to faster PCM melting and higher stored energy and charging power. The lowest melting duration was obtained for ε_avg = 0.8.
• At a constant average porosity, using a higher positive porosity gradient a, defined from the left hot wall toward the right, enhances conductive heat transfer near the left region with low porosity in initial stages and the convective heat transfer in the right side of the cavity in the final stages. The opposite occurs when lower and negative values of a are used. Melting duration is minimized for a = 4.
• Raising the PPI presents two effects. First, it reduces the thickness of the solid matrix, which diminishes effective thermal conductivity. Second, it increases surface contact and improves the heat transfer between metal foam and PCM. The balance between these two effects leads to a slight overall impact of PPI on the charging power.
• Tilting the cavity has a limited effect on the flow and thermal behavior of NePCM, mainly in the initial stages. This effect becomes more apparent when the convective heat transfer starts taking place, and the charging power can be improved when the inclination angle is γ = π/4 or γ = 3π/4 compared with the cases of γ = 0 and γ = π.

The focus of the present study was thermal energy storage, which is concerned with the melting of PCM. Investigation of solidification (discharging) in a non-uniform metal foam could be subject to future studies.