Computational Modeling of Latent Heat Thermal Energy Storage in a Shell-Tube Unit: Using Neural Networks and Anisotropic Metal Foam

Abstract

Latent heat storage in a shell-tube is a promising method to store excessive solar heat for later use. The shell-tube unit is filled with a phase change material PCM combined with a high porosity anisotropic copper metal foam (FM) of high thermal conductivity. The PCM-MF composite was modeled as an anisotropic porous medium. Then, a two-heat equation mathematical model, a local thermal non-equilibrium approach LTNE, was adopted to consider the effects of the difference between the thermal conductivities of the PCM and the copper foam. The Darcy–Brinkman–Forchheimer formulation was employed to model the natural convection circulations in the molten PCM region. The thermal conductivity and the permeability of the porous medium were a function of an anisotropic angle. The finite element method was employed to integrate the governing equations. A neural network model was successfully applied to learn the transient physical behavior of the storage unit. The neural network was trained using 4998 sample data. Then, the trained neural network was utilized to map the relationship between control parameters and melting behavior to optimize the storage design. The impact of the anisotropic angle and the inlet pressure of heat transfer fluid (HTF) was addressed on the thermal energy storage of the storage unit. Moreover, an artificial neural network was successfully utilized to learn the transient behavior of the thermal storage unit for various combinations of control parameters and map the storage behavior. The results showed that the anisotropy angle significantly affects the energy storage time. The melting volume fraction MVF was maximum for a zero anisotropic angle where the local thermal conductivity was maximum perpendicular to the heated tube. An optimum storage rate could be obtained for an anisotropic angle smaller than 45°. Compared to a uniform MF, utilizing an optimum anisotropic angle could reduce the melting time by about 7% without impacting the unit’s thermal energy storage capacity or adding weight.

1. Introduction

Thermal energy storage (TES) has found essential applications in the automotive, power generation, buildings, and food industries [1]. TES technology has shown great potential to regulate the difference between the supply and demand chain. TES technology promotes the flexibility, efficiency, and reliability of thermal energy systems by regulating real-world heat sources’ fluctuation and intermittency nature [2]. Latent thermal energy storage (LHTES) utilizes the latent heat of fusion and store/release thermal energy at a constant temperature in a compact space. Thus, LHTES systems are promising for large- and small-scale applications. LHTES units have been utilized in buildings [3,4], solar systems [5], solar air heating [6], solar dryers [7], and waste heat recovery [8].

Although LHTES systems can store/release a considerable amount of latent heat at a constant fusion temperature, they suffer from poor thermal conductivity of conventional phase change material (PCM) [5,9]. Thus various methods, such as using fins [10,11,12], nanoparticles [13,14], metal foams [15,16], or a combination of these methods [17]. Embedding PCMs in open-cell metal foams is an excellent approach for improving the heat transfer performance of LHTES systems, which has attracted the attention of recent researchers [18]. However, metal foams do not undergo a phase change and cannot contribute to latent heat TES. Moreover, metal foams are made of thermally conductive metals such as copper, aluminum, and nickel, which increase the cost and weight of the LHTES unit. Thus, some researchers tried to utilize layers of metal foam to improve the system’s thermal performance by using a minimum amount of metal foam.

For example, Zhang et al. [19] experimentally investigated the influence of using layers of metal foams in an enclosure. They also investigated the impact of utilizing different layers of metal foam on the heat transfer behavior of an LHTES unit. The results showed that filling the enclosure with metal foam up to 5/6 of its capacity improves heat transfer. Using different types of foam layers could also improve heat transfer. A metal foam layer with a high porosity should be placed next to the heated wall or at the top regions of the enclosure unit. Buonomo et al. [20] studied the effect of a metal foam layer on the thermal performance of a shell-tube shape LHTES unit where the foam layer was placed on the tube. The authors examined the effect of various porosities in the range of 0.9–0.97 and various foam thicknesses on the melting time. An increase in the thickness of the foam layer could change the heat transfer regimes and reduce the melting time. Ge et al. [21] utilized several optimized shapes of fins in a horizontal shell-tube shape thermal energy storage. They compared the advantage of fins with metal foams. They reported that using fins or metal foams (fully or half-filled enclosures) could improve the heat transfer rate by 3.3–5.8 times. They demonstrated that the topology-optimized fins could provide better thermal enhancement compared to metal foams. However, the production cost of optimized fins was higher. Hashem Zadeh et al. [22] utilized different configurations of fin shape layers of metal foams to improve heat transfer in a horizontal shell-tube shape LHTES unit. A Y shape design could produce the best charging performance.

Some researchers tried to change the structure of metal foams and produce anisotropic metal foams with improved thermal conductivity in a direction. In this approach, the amount of metal foam is fixed, and the changes are applied to the internal structure of the foam. Bamdezh et al. [23] examined the metal foam anisotropy on the thermal behavior of a thermal management system. The authors reported that an increase in the tangential conductance raises the phase change rate. Moreover, promoting axial thermal conductivity enhances heat transfer. Yu et al. [24] explored the anisotropic properties of porous media with two different pore structures and found anisotropic characteristics important. Ren et al. [25] reported the anisotropic thermal conductivity effects in a woven metal foam PCM composite. They reported a better thermal performance for an anisotropic than an isotropic for thermal energy storage applications. They found that a metal foam with engineered anisotropic properties could notably promote heat transfer performance and required material.

The few available literature works on anisotropic metal foams reveal the potential advantage of these materials for LHTES. However, there are only a few experimental or numerical investigations on the engineering design of LHTES units with anisotropic metal foams. The literature review showed that shell-tube type LHTES has many applications in solar energy storage systems. Thus, the present study aims to address the impact of foam anisotropic characteristics on the phase change behavior of shell-tube LHTES units for the first time.

2. Mathematical Model

A shell-tube type latent heat thermal energy storage stores/releases thermal energy in thermodynamic cycles. For example, Figure 1 illustrates a view of a solar water heating system. In this system, a cold-water source supplies the building’s water demand. The cold water enters a hot water storage tank, providing hot water to the house. The water in the tank circulates in a cycle and goes through solar collectors in a solar heating cycle as water absorbs solar heat and its temperature rises. The excess solar energy can be stored in a compact space, an LHTES unit, in the form of sensible and latent heat. Thus, an LHTES cycle circulates the tanks’ hot water in the LHTES units and charges the PCMs inside the unit. Later, at any point, some of the tank’s water can be consumed by the building, and cold fresh water will enter the tank, causing a significant temperature drop in the hot water tank. In this case, the LHTES cycle circulates the middle-temperature water of the tank in the LHTES unit. Thus, the PCM in the unit releases its energy to the water and undergoes a solidification process.

The heat transfer capabilities of the LHTES unit it is an important factor in thermodynamic cycles since it dictates the required time to regulate the tank’s temperature. Here, it assumed that the thermal charging of LHTES is set at 70 °C, and the LHTES loops activate at this temperature. Thus, by neglecting the line energy losses, the inlet temperature to LHTES can be considered 70 °C.

The storage capacity of the cycle can be adjusted by using several parallel LHTES units. Since the units can be mounted in parallel, their thermal behavior can be similar, and thus, the present study focuses on the modeling and analysis of just one LHTES unit. A schematic view of the LHTES unit is depicted in Figure 2. A shell-tube storage unit is filled by the anisotropic copper foam of porosity ε. The thermal conductivity (k) and permeability (κ) of metal foam (MF) can be changed in perpendicular directions. The MF is filled by paraffin with a melting temperature of T_f. The tube is made of copper with a nominal diameter of D= (½)” or D = (2.54/2) cm [26]. A relative pressure P_in drives the HTF water into the HTF tube at an entrance temperature of T_h. The LHTES unit is initially at a supper cold temperature T_c. The tube’s outside diameter, inside diameter, and wall thickness are 1.59 cm, 1.34 cm, and 0.13 cm, respectively. The geometrical details of the shell are reported in Figure 1 as a function of tube nominal radius (R = D/2).

2.1. Governing Equations

There are three regions of HTF, copper tube wall, and composite MF-PCM, which are physically coupled. The governing equation for the copper wall is the conduction heat transfer which can be written as in a cylindrical coordinate system:

$${(\rho C_p)}_{\mathrm{Wall}}\frac{\partial T}{\partial t}=k_{\mathrm{Wall}}\frac{{\partial }^{2}T}{\partial z^2}+\frac{k_{\mathrm{Wall}}}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)$$

(1)

where ρ, C_p, and k are the tube wall’s density, specific heat capacity, and thermal conductivity, respectively. Moreover, T and t are the temperature field and time, respectively. The cylindrical coordinate system is introduced by the r- and z-axis, as indicated in Figure 1. The subscript wall indicates the copper wall.

The governing equations for the HTF region are the continuity, momentum, and convection heat transfer as follows:

$$\frac{\partial u_z}{\partial z}+\frac{1}{r}\frac{\partial \left(r{u}_r\right)}{\partial r}=0$$

(2)

$$\begin{gathered} {\rho }_{\mathrm{HTF}}\left(\frac{\partial u_r}{\partial t}+u_z\frac{\partial u_r}{\partial z}+u_r\frac{\partial u_r}{\partial r}\right)=-\left(\frac{\partial p}{\partial r}\right)+{\mu }_{\mathrm{HTF}}\frac{{\partial }^2{u}_r}{\partial z^2}+\frac{{\mu }_{\mathrm{HTF}}}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_r}{\partial r}\right)-\frac{{\mu }_{\mathrm{HTF}}{u}_r}{r^2} \\ {\rho }_{\mathrm{HTF}}\left(\frac{\partial u_z}{\partial t}+u_z\frac{\partial u_z}{\partial z}+u_r\frac{\partial u_z}{\partial r}\right)=-\left(\frac{\partial p}{\partial z}\right)+{\mu }_{\mathrm{HTF}}\frac{{\partial }^2{u}_z}{\partial z^2}+\frac{{\mu }_{\mathrm{HTF}}}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_z}{\partial r}\right) \end{gathered}$$

(3)

$${(\rho C_p)}_{\mathrm{HTF}}\left(\frac{\partial T}{\partial t}+u_z\frac{\partial T}{\partial z}+u_r\frac{\partial T}{\partial r}\right)=k_{\mathrm{HTF}}\frac{{\partial }^{2}T}{\partial z^2}+\frac{k_{HTF}}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)$$

(4)

where the velocity components in the r and z directions are denoted by u_r and u_z, respectively. The fluid’s dynamic viscosity is μ, and the subscript HTF represents the heat transfer liquid. Here, it was assumed that the thermophysical properties are constant. Moreover, the fluid flow is laminar (Re_HTF < 2000) where Re_HTF = ρ_HTF × u × D/μ_HTF, and u =$\sqrt{u_r^2+u_z^2}$. It should be noted that there is no gravity effect in the HTF governing equations since it was assumed that the HTF works in a closed loop, and the static pressure does not contribute to the fluid movement.

Regarding the composite MF-PCM, the metal foam is made of open cells, and thus the liquid PCM can move inside the pores. In the molten PCM, there are temperature gradients; hence, the buoyancy forces can induce natural convection flow circulation. Thus, the continuity and momentum equations should be solved along with the phase change heat equation to simulate the natural convection heat transfer with phase change. Here, the enthalpy porosity formulation is used to simulate phase change. Source terms as a function of melt volume fraction (φ) are added to the momentum equation to force velocities to zero in solid regions. Here, a difference between the temperature of MF porous matrix (T_MF) and the temperature of the PCM (T_PCM) inside the pores was considered. Thus, a local-thermal non-equilibrium (LTNE) model was used to explain the temperature differences. Considering the non-Darcy effects and controlling source terms, the governing equations for continuity, momentum, and phase change heat transfer can be expressed as follows [27,28,29]:

$$\frac{\partial u_z}{\partial z}+\frac{1}{r}\frac{\partial \left(r{u}_r\right)}{\partial r}=0$$

(5)

$$\begin{array}{l}\frac{{\rho }_{\mathrm{PCM}}}{\epsilon }\frac{\partial u_r}{\partial t}+\frac{{\rho }_{\mathrm{PCM}}}{{\epsilon }^2}\left(u_z\frac{\partial u_r}{\partial z}+u_r\frac{\partial u_r}{\partial r}\right)=-\left(\frac{\partial p}{\partial r}\right)+\frac{{\mu }_{\mathrm{PCM}}}{\epsilon }\left[\frac{{\partial }^2{u}_r}{\partial z^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_r}{\partial r}\right)-\frac{u_r}{r^2}\right]\\ -\frac{{\mu }_{\mathrm{PCM}}}{\kappa }{u}_r-{\rho }_{\mathrm{PCM}}\frac{C_F}{\sqrt{\kappa }}\left|u\right|u_r+A_{mush}\frac{{\left(1-\phi \left(T\right)\right)}^2}{{\lambda }_{mush}+{\phi }^3\left(T\right)}{u}_r\end{array}$$

(6a)

$$\begin{array}{l}\frac{{\rho }_{\mathrm{PCM}}}{\epsilon }\frac{\partial u_z}{\partial t}+\frac{{\rho }_{\mathrm{PCM}}}{{\epsilon }^2}\left(u_z\frac{\partial u_z}{\partial z}+u_r\frac{\partial u_z}{\partial r}\right)=-\left(\frac{\partial p}{\partial z}\right)+\frac{{\mu }_{\mathrm{PCM}}}{\epsilon }\left[\frac{{\partial }^2{u}_z}{\partial z^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_z}{\partial z}\right)\right]\\ +g{\rho }_{\mathrm{PCM}}{\beta }_{\mathrm{PCM}}\left(T-T_0\right)-\frac{{\mu }_{\mathrm{PCM}}}{\kappa }{u}_z-{\rho }_{\mathrm{PCM}}\frac{C_F}{\sqrt{\kappa }}\left|u\right|u_z+A_{mush}\frac{{\left(1-\phi \left(T\right)\right)}^2}{{\lambda }_{mush}+{\phi }^3\left(T\right)}{u}_z\end{array}$$

(6b)

$$\begin{array}{l}\epsilon {\left(\rho C_p\right)}_{\mathrm{PCM}}\frac{\partial T_{\mathrm{PCM}}}{\partial t}+(\rho C_p{)}_{\mathrm{PCM}}\left(u_z\frac{\partial T_{\mathrm{PCM}}}{\partial z}+u_r\frac{\partial T_{\mathrm{PCM}}}{\partial r}\right)=\\ \frac{\partial }{\partial z}\left(k_{\mathrm{eff},\hspace{0.17em}\mathrm{PCM}}\frac{\partial T_{\mathrm{PCM}}}{\partial z}\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(k_{\mathrm{eff}, \mathrm{PCM}}r\frac{\partial T_{\mathrm{PCM}}}{\partial r}\right)+h_v\left(T_{\mathrm{MF}}-T_{\mathrm{PCM}}\right)-\epsilon {\rho }_{\mathrm{PCM}}{L}_{\mathrm{PCM}}\frac{\partial \phi \left(T\right)}{\partial t}\end{array}$$

(7a)

$$\left(1-\epsilon \right){\left(\rho C_p\right)}_{\mathrm{MF}}\frac{\partial T_{\mathrm{MF}}}{\partial t}=\frac{\partial }{\partial z}\left(k_{\mathrm{eff},\hspace{0.17em}\mathrm{MF}}\frac{\partial T_{\mathrm{MF}}}{\partial z}\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(k_{\mathrm{eff}, \mathrm{MF}}r\frac{\partial T_{\mathrm{MF}}}{\partial r}\right)-h_v\left(T_{\mathrm{MF}}-T_{\mathrm{PCM}}\right)$$

(7b)

where κ and C_F are the porous permeability and Frochheimer parameter, respectively. The volumetric interstitial heat transfer coefficient (h_v) is introduced later using the pore scale Nusselt number. The subscripts MF and eff denoted the metal foam and effective properties, respectively. The symbols β, g, and L denote the thermal volume expansion coefficient, gravity acceleration, and latent heat of fusion, respectively. The source terms containing A_mush and λ_mush force the fluid velocity to zero in the solid regions. Thus, a substantial value of A_mush = 10¹⁰ Pa.s/m² was adopted to adequately force the fluid velocities to zero by dominating the flow solution in solid regions.

Moreover, λ_mush = 0.001 was used to avoid a zero denominator in the source terms. The dynamic viscosity was considered as a function of melt fraction (φ) to better control the solver stability and force the velocity to zero in solid regions. Thus, the viscosity was introduced as μ_PCM = φ × μ_PCM,l + (1 − φ) × μ_a where μ_a is an artificially large viscosity value. It was adopted as 10⁴ Pa.s in the current study. Using this formula, the viscosity tends to the typical dynamic viscosity μ_PCM,l when φ = 1 (liquid region), and it artificially increases to a substantial value in a solid region (φ = 0). Since the aim of source terms is to force the velocities to zero in a solid region, the artificial increase in the dynamic viscosity helps force the velocities to zero in a solid region. Moreover, such a dynamic viscosity formulation helps with the solver stability and facilitates numerical computations without impacting the physical model.

Equation (7a,b) represent heat transfer in the PCM and MF phases. The term h_v (T_MF − T_PCM) considers the thermal interaction of two phases. The effective thermal conductivities for PCM and MF are used in these equations since the pore structures influence these thermal conductivities [30]. Melting volume fraction, φ, is a function of PCM temperature as [31]:

$$\phi \left(T\right)=\{\begin{array}{l}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}T<T_{\mathrm{f}}-0.5\mathsf{\Delta }{T}_{\mathrm{f}}\hspace{0.17em}(\mathrm{Solid}\hspace{0.17em}\mathrm{phase})\\ 0.5+\frac{\left(T-T_{\mathrm{f}}\right)}{\mathsf{\Delta }{T}_{\mathrm{f}}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{T}_{\mathrm{f}}-0.5\mathsf{\Delta }{T}_{\mathrm{f}}\le T\le T_{\mathrm{f}}+0.5\mathsf{\Delta }{T}_{\mathrm{f}}\hspace{0.17em}\hspace{0.17em}(\mathrm{Mushy}\hspace{0.17em}\mathrm{region})\\ 1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}T>T_{\mathrm{f}}+0.5\mathsf{\Delta }{T}_{\mathrm{f}}\hspace{0.17em}(\mathrm{Liquid}\hspace{0.17em}\mathrm{phase})\end{array}$$

(8)

where ΔT_f is the phase transition interval around the fusion temperature T_f. MF’s thermal conductivity and permeability are anisotropic, and they can change in perpendicular directions. Supporting the MF ligaments in one direction enhances the thermal conductivity with the cost of reducing MF permeability. Hence, they have been introduced using a second-order tensor follows [32]:

$$\kappa =\hspace{0.17em}\left[\begin{array}{l}{\kappa }_2{\left(\sin \theta \right)}^2+{\kappa }_1{\left(\cos \theta \right)}^2\hspace{1em}\hspace{1em}\left({\kappa }_1-{\kappa }_2\right)\left(\cos \theta \right)\left(\sin \theta \right)\\ \left({\kappa }_1-{\kappa }_2\right)\left(\cos \theta \right)\left(\sin \theta \right)\hspace{1em}\hspace{1em}{\kappa }_2{\left(\cos \theta \right)}^2+{\kappa }_1{\left(\sin \theta \right)}^2\end{array}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(9a)

$$k_{\mathrm{eff},\hspace{0.17em}\mathrm{MF}}=\hspace{0.17em}\left[\begin{array}{l}{k}_2{\left(\sin \theta \right)}^2+k_1{\left(\cos \theta \right)}^2\hspace{1em}\hspace{1em}\left(k_1-k_2\right)\left(\sin \theta \right)\left(\cos \theta \right)\\ \left(k_1-k_2\right)\left(\sin \theta \right)\left(\cos \theta \right)\hspace{1em}\hspace{1em}{k}_2{\left(\cos \theta \right)}^2+k_1{\left(\sin \theta \right)}^2\end{array}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(9b)

where κ₁ = (1 − ζ)×κ_a, κ₂ = (1 + ζ) × κ_a, k₁ = (1 + ζ)×k_a, and k₂ = (1 − ζ) × k_a. The anisotropic angle, θ, is defined in Figure 2b, and it represents the angle of anisotropic MF placement in the enclosure. Here, ζ is the anisotropic factor, where ζ = 0 represents an isotropic MF. The anisotropic factor (ζ) cannot be too high since the MF should stay connected and continuous. Thus, in the present study, the anisotropic factor is adopted as a moderate value of ζ = 0.2 [32]. The average permeability (κ_a) [33] and average MF effective thermal conductivity (k_a) of the MF are computed using [27,30]:

$$k_a=\frac{1-\epsilon }{3}{k}_{\mathrm{MF}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(10a)

$${\kappa }_a=\frac{{\epsilon }^2{\left(d_{fp}\sqrt{\frac{{\kappa }_{tor}}{3\epsilon }}\right)}^2}{36\left({\kappa }_{tor}-1\right){\kappa }_{tor}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(10b)

where k_MF is the bulk thermal conductivity of MF material and d_fp and κ_tor are introduced as [33]:

$$d_{fs}=1.18\sqrt{\frac{1-\epsilon }{3\pi }}\left\{\frac{1}{1-e^{\frac{\left(\epsilon -1\right)}{0.04}}}\right\}{d}_{fp}$$

(11a)

$$\frac{1}{{\kappa }_{tor}}=\frac{3}{4\epsilon }+\frac{\sqrt{9-8\epsilon }}{2\epsilon }\cos \left\{\frac{4\pi }{3}+\frac{1}{3}{\cos }^{-1}\left(\frac{8{\epsilon }^2-36\epsilon +27}{{\left(9-8\epsilon \right)}^{\frac{3}{2}}}\right)\right\}{d}_{fp}$$

(11b)

in which, d_fp is computed using MF pore per inch (PPI) property as [33]:

$$d_{fp}=\frac{0.0254}{PPI}$$

(11c)

The Frochheimer parameter is computed using [33]:

$$C_F=0.00212{\left(1-\epsilon \right)}^{-0.132}{\left(\frac{d_{fs}}{d_{fp}}\right)}^{-1.63}$$

(12)

The heat exchange between the PCM and MF was evaluated using the volumetric interstitial heat transfer coefficient as [34]:

$$h_v=N{u}_v\frac{k_{\mathrm{PCM}}}{d_{fs}^2}$$

(13)

where Nu_v is the interstitial heat transfer Nusselt number introduced in a volumetric form [35]. Then, Nu_v is computed for low pore scale Reynold numbers as [34]:

$$N{u}_v=\{\begin{array}{l}76.99-152.01\epsilon +75.04{\epsilon }^2\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0\le \mathrm{Re}\le 0.1\\ \left(1.72+1.71\epsilon -3.46{\epsilon }^2\right){\mathrm{Re}}^{0.26}{\mathrm{Pr}}^{0.28}\hspace{1em}\hspace{1em}\hspace{1em}0.1<\mathrm{Re}\le 1\end{array}$$

(14)

where Pr = ρ_PCM × μ_PCM/α_PCM and α_PCM = k_PCM/(ρC_p)_PCM. The pore scale Reynold number (Re) was also introduced as Re = ρ_PCM × u_PCM × d_fs/μ_PCM. Many past literature studies used different relations to estimate the Nusselt number. However, most of the available relationships are related to high values of the Reynolds number and typically Re > 1. However, for the natural convection of PCM in MF, the circulation velocities are small, and Re is typically smaller than unity. Thus, a relation for the specific range of Re < 1 should be used. Therefore, here, the relation reported by Yao et al. [34] was adopted. This relation was texted against the experimental observations and agreed well with the phase change in PCMs in MF. Equation (16) holds for 21 < Pr < 41 and 0.929 < ε < 0.974. Moreover, it should be noted that the heat transfer exchange (Nu_v) is not zero when the velocity is very small or zero since there is still an active thermal diffusion heat transfer which can play a significant role in the thermal exchange between the PCM and MF phases.

Finally, the effective thermal conductivity of PCM was also evaluated [27,30]:

$$k_{\mathrm{eff},\hspace{0.17em}\mathrm{PCM}}=\frac{2+\epsilon }{3}{k}_{\mathrm{PCM}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(15a)

A linear weight average was used to compute the thermophysical properties of the PCM in the mushy region:

$${\left(\mathsf{\rho }{C}_p\right)}_{\mathrm{PCM}}=\left(1-\phi \right){\left(\mathsf{\rho }{C}_p\right)}_{\mathrm{l}}+\phi {\left(\mathsf{\rho }{C}_p\right)}_{\mathrm{s}}$$

(15b)

$${\rho }_{\mathrm{PCM}}=\left(1-\phi \right){\rho }_{\mathrm{l}}+\phi {\rho }_{\mathrm{s}}$$

(15c)

where subscripts l and s indicate the PCM in the liquid and solid phases.

Table 1. Thermophysical properties of PCM, metal foam, and HTF fluid.

Materials	K (W/m.K)	ρ (kg/m3)	cp (J/kg.K)	Tm (°C)	L (kJ/kg)	μf (kg/m.s)	Β (1/K)
Paraffin (solid/liquid) [36,37,38]	0.21/0.12	916/790	2700/2900	49–54	176	0.0036	0.00091
Copper foam [39]	380	8900	386	-	-	-	-
Water [40]	0.613	997.1	4179	-	-	0.000957	0.00021

summarizes the thermophysical properties of the tube, MF, and paraffin.

2.2. Initial and Boundary Conditions

The continuity of temperature and heat flux was considered for all coupled interfaces. For the MF region with LTNE conditions, the continuity of temperature was applied. For the continuity of the heat flux (q), it was divided into two parts based on the volume fraction of each phase (porosity). Thus, the heat flux continuity at the interface of the tube wall and MF region was introduced as:

$$q_{\mathrm{Wall}}=\epsilon q_{\mathrm{MF}}+\left(1-\epsilon \right)q_{\mathrm{PCM}}$$

(16)

A uniform inlet temperature of T_h was applied as T = 70 °C. A constant relative pressure for p = P_in was applied at the HTF intel boundary. The outflow condition with a relative pressure zero and −n.q = 0 was applied at the outlet where n is the normal vector to a surface. The shell walls are assumed well insulated with −n.q = 0. The no-slip and impermeability boundary conditions were employed for all boundary surfaces except the inlet and outlet boundaries which were explained. The entire LHTES unit was assumed at an initial temperature of T_c = 14 °C. A reference pressure point with zero relative pressure was considered at the bottom left corner of the shell space. The thermophysical properties of HTF fluid (water) and copper foam can be varied slightly by temperature variations.

Since the temperature variation is limited in the present study, the thermophysical properties were assumed constant except for the phase change material which undergoes a phase change. Moreover, the impact of the temperature variations on the density change was included in the buoyancy effect using the Boussinesq model in the liquid PCM.

2.3. Characteristics Parameters

The melting volume fraction shows the average amount of molten PCM in the enclosure, and it was defined as:

$$MVF=\frac{{\displaystyle {\oint }_V\epsilon \phi dV}}{{\displaystyle {\oint }_V\epsilon dV}}$$

(17)

where dV is the volume element of the shell domain. The average outlet temperature difference is introduced as:

$$\mathsf{\Delta }T=T_h-\frac{{\displaystyle {\oint }_A{u}_{z}TdA}}{{\displaystyle {\oint }_A{u}_{z}dA}}$$

(18)

where A is the surface area of the outlet. The stored thermal energy was computed as the sum of sensible and latent heat:

$$\begin{array}{l}{Q}_{sensible}=\left[{\displaystyle {\oint }_V\left({\displaystyle {\int }_{T_0}^T\epsilon {\left(\rho C_P\right)}_{\mathrm{PCM}}(T)dT}\right)dV}\right]+\\ \left(T-T_0\right){\left(\rho C_P\right)}_{\mathrm{MF}}{\displaystyle {\oint }_V\left(1-\epsilon \right)dV}+\\ \left(T-T_0\right){\left(\rho C_P\right)}_{\mathrm{Wall}}{V}_{\mathrm{Wall}}+\left(T-T_0\right){\left(\rho C_P\right)}_{\mathrm{HTF}}{V}_{\mathrm{HTF}}\end{array}$$

(19)

$$Q_{latent}=\epsilon {\displaystyle {\oint }_V{\rho }_{\mathrm{PCM}}\phi L_{\mathrm{PCM}}dV}$$

(20)

Q_store = Q_sensible + Q_latent

(21)

The average velocity of HTF (u_a) leaving the unit is computed as $u_a={\displaystyle {\oint }_A{u}_{z}dA}/{\displaystyle {\oint }_{A}dA}$. It should be noted that the present study is an axis-symmetric model, and the 2D lines are axis-symmetric surfaces, and the 2D surfaces are indeed axis-symmetric volumes. Using the HTF average velocity, the HTF tube Reynolds number (Re_HTF) is introduced as Re_HTF = ρ_HTF × u_a(D_in)/μ_HTF where D_in is the inside diameter of the HTF tube.

Energy storage power

Power = Q_store/t

(22)

The energy storage power indicates the average power of thermal energy storage from the initial state (super cold PCM) to any given time during the charging process.

3. Solution Method

3.1. Numerical Method

The governing equations in each domain are integrated over a structured mesh using the finite element method [41] in their weak form to obtain a set of residual equations. In the present governing equations, there are strong multi-physics non-linear source/sink terms due to the phase change effects. Thus, a numerical method capable of dealing with such non-linearities should be utilized. Here, the finite element provides a smooth, accurate, and continuous approximate solution over mesh elements, which is suitable for the purpose of the present study.

In the finite element method, the Gauss quadrature integration over the linear shape function discretization was used. The field variables are introduced using shape functions at each time step as:

$$\begin{array}{l}{\mathrm{u}}_{\mathrm{r}}\approx {\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{u}}_{\mathrm{r},\mathrm{i}}{\mathsf{\gamma }}_{\mathrm{i}}\left(\mathrm{r},\mathrm{z}\right),}{\mathrm{u}}_{\mathrm{z}}\approx {\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{u}}_{\mathrm{z},}{}_{\mathrm{i}}{\mathsf{\gamma }}_{\mathrm{i}}\left(\mathrm{r},\mathrm{z}\right),}\mathrm{p}\approx {\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{p}}_{\mathrm{i}}{\mathsf{\gamma }}_{\mathrm{i}}\left(\mathrm{r},\mathrm{z}\right),}\\ \mathrm{T}\approx {\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{T}}_{\mathrm{i}}{\mathsf{\gamma }}_{\mathrm{i}}\left(\mathrm{r},\mathrm{z}\right),}\hspace{0.17em}{\mathrm{T}}_{\mathrm{M}\mathrm{F}}\approx {\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{T}}_{\mathrm{M}\mathrm{F},\mathrm{i}}{\mathsf{\gamma }}_{\mathrm{i}}\left(\mathrm{r},\mathrm{z}\right),}{\mathrm{T}}_{\mathrm{P}\mathrm{C}\mathrm{M}}\approx {\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{T}}_{\mathrm{P}\mathrm{C}\mathrm{M},\mathrm{i}}{\mathsf{\gamma }}_{\mathrm{i}}\left(\mathrm{r},\mathrm{z}\right),}\end{array}$$

(23)

where γ is the shape function and is the number of discretized regions. The finite element algebraic equations were computed by inserting the above-filed variables into the weak form of the governing equations. Then, the algebraic equations were solved using the Newton method to find the variable fields.

The PARDISO parallel solver [42,43] was invoked to solve the residual equations over multiple computational cores. A constant Newtonian damping factor of 0.9 and relative error tolerance of 10⁻⁴ were utilized. The solution time steps and convergence were automatically controlled to regulate solution errors within the predefined error tolerance (10⁻⁴) by the first-second order backward differential formula (BDF) [44]. The convergence of the numerical method was monitored through the evaluation of the residual equations. The BDF scheme monitors the convergence solution at each time step for a relative error to remain in the specific threshold and also to select the next time step. The computations were commenced by employing the initial conditions and continued with the computation of the model parameters. The continuity, momentum, and phase change heat equations are integrated and solved in a fully coupled method. A stop condition of MVF ≥ 0.999 was added to terminate the computations when a completed melting condition reaches.

The solution steps are summarized in a flow chart depicted in Figure 3.

Moreover, the momentum equations in the HTF tube were solved as a steady-state step since they are independent of the thermal operating conditions of the LHTES unit. Then, the computed velocities were used in the time-dependent simulations as constant initial values. This step reduces the computational costs without affecting the simulation results.

3.2. Mesh Study

A non-uniform unstructured mesh was created in the HTF domain. Moreover, the wall and PCM regions were meshed using a uniform, unstructured mesh. A view of the utilized mesh is displayed in Figure 4.

The number of mesh cells in each domain was 45N_m × 5N_m (HTF), 45N_m × N_m (Tube wall), and 45N_m × 30N_m (PCM). The mesh in the HTF domain was stretched from the tube wall toward the axis with an element ratio of 10. The simulations were repeated for four cases with N_m equal to 2, 3, 4, and 5 when P_in = 5 Pa, ζ = 0.2, and θ = 60°.

The values of the obtained results are summarized in

Table 2. The computed characteristics for various mesh sizes at t = 4 h when P_in = 5 Pa, ζ = 0.2, and θ = 60°. The computational time was reported for a full melting. A mesh with N_m ≥ 3 provides results with 0.5% error or smaller, which is acceptable for most engineering applications.

Nm	MVF	%err	Qstore (kJ)	%err	Computational Time (min)
2	0.6958	1.6	1422.7	0.9	34
3	0.6883	0.5	1414.2	0.3	74
4	0.6860	0.1	1411.5	0.1	135
5	0.6851	0.0	1410.6	0.0	217

. The error was computed with respect to the finest mesh. For example, the error for MVF was computed as err = 100 × (MVF-MVF_{@Nm = 5})/MVF_{@Nm = 5}. As seen, a mesh with N_m = 3 or higher can compute MVF and stored energy with 0.5% error or less. Figure 5 depicts the time history of MVF for various mesh sizes, which agrees with the results of Table 2 over the melting process. Thus, as a fair trade between computational cost and accuracy, the mesh with N_m = 3 was adopted for the computation of the results section.

3.3. Model Verification

Zheng et al. [39] empirically examined the melting of paraffin wax in a square enclosure filled with copper metal foam. The paraffin in the enclosure was initially at a super cold temperature of 14 °C and fully solid. The enclosure’s height and width were 100 mm. The left enclosure was heated by an electrical heater with a heat flux of 1150 W/m2. The PCM and MF in the enclosure absorb the heat, and the temperature in the enclosure increases gradually until the PCM reaches the melting temperature T_f = 55.3°. Figure 6 illustrates a comparison between the simulated results of the present study and the observations of Zheng et al. [39]. The porosity was ε = 0.95 with PPI = 5 and ζ = 0. The permeability and the effective thermal conductivity of metal foam (k_{MF, eff}) were computed at 8.3857 × 10⁻⁷ m² and 6.33 W/(m·K), respectively.

The temperatures were measured at the heated wall and over a vertical line with a 2.5 cm distance from the heated wall. Figure 7 illustrates a comparison between the time history of the measured ([39]) and computed (present study) temperatures. The simulations are computed for local thermal equilibrium (LTE) and LTNE models. LTE was obtained by considering a significant interaction between phases (h_v = 10¹² W/m³). The results show the LTNE is only significant near the heated wall while there is no significant difference between porous and PCM far from the wall. LTNE model provides a temperature prediction closer to experiments.

As another comparison, the free convection heat transfer in an enclosure with the temperature difference between side walls was simulated. The enclosure was filled with an anisotropic porous medium [45]. In [45], the thermal conductivity was a function of the material space, k(x, y), and it was introduced with a gradient in x- and y- directions depending on a distribution parameter, ζ. Figure 8 compares the non-dimensional temperatures of the present computations and those of [45] when ζ = 2.0. Figure 8 depicts the excellent proximity of the present computations and the literature results.

4. Results and Discussions

The impact of the MF anisotropic angle (θ) and the inlet pressure P_in were investigated on the charging time and the thermal behavior of the LHTES unit. As explained in the modeling section, the inlet temperature (Th) and anisotropic parameter (ζ) were kept constant. The anisotropic angle (θ) was addressed for −90°, −75°, −60°, −45°, −30°, −15°, and 0°, and positive values of these angles. The simulations were also performed for three inlet pressures (P_in) of 5 Pa, 7.5 Pa, and 10 Pa. Moreover, three reference cases of ζ = 0 and P_in = 5, 7.5, and 10 Pa were investigated. The investigated cases are summarized in

Table 3. Detail of simulated cases for various values of inlet pressure and anisotropic angles and reference cases when ζ = 0.2.

Case	ζ	Pin	θ	Case	ζ	Pin	θ
1	0.2	5	−90	22	0.2	7.5	30
2	0.2	5	−75	23	0.2	7.5	45
3	0.2	5	−60	24	0.2	7.5	60
4	0.2	5	−45	25	0.2	7.5	75
5	0.2	5	−30	26	0.2	7.5	90
6	0.2	5	−15	27	0.2	10	−90
7	0.2	5	0	28	0.2	10	−75
8	0.2	5	15	29	0.2	10	−60
9	0.2	5	30	30	0.2	10	−45
10	0.2	5	45	31	0.2	10	−30
11	0.2	5	60	32	0.2	10	−15
12	0.2	5	75	33	0.2	10	0
13	0.2	5	90	34	0.2	10	15
14	0.2	7.5	−90	35	0.2	10	30
15	0.2	7.5	−75	36	0.2	10	45
16	0.2	7.5	−60	37	0.2	10	60
17	0.2	7.5	−45	38	0.2	10	75
18	0.2	7.5	−30	39	0.2	10	90
19	0.2	7.5	−15	40	0	5	0
20	0.2	7.5	0	41	0	7.5	0
21	0.2	7.5	15	42	0	10	0

.

Figure 9 depicts the MVF and stored energy over time for various anisotropic angles. As seen, the decrease in the anisotropic angle from 90° to 0° accelerates the melting process and reduces the melting time. The results for negative angles are very close to similar positive angles. This is since the difference between the positive and negative angles is just the minor difference produced by sin(θ)cos(θ) terms in the thermal conductivity and permeability matrixes. An isotropic case, a reference case with the average porosity and permeability, was also simulated, and the results were added to this figure. The results show that the typical isotropic MF (reference case) can produce a melting process between the anisotropic MF when the anisotropic angle is |θ| > 45°. However, the decline of the anisotropic angle below 45° down to 0° provides an accelerated melting process.

In practice, reaching different anisotropic angles is easy since it can change by the mounting-direction of the anisotropic MF in an enclosure. The reason for the improvement of melting heat transfer by the increase in anisotropic angle is attributed to the conductivity matrix. Using Equation (2), the thermal conductivity and permeability matrix can be obtained as follows for θ = 0° and θ = 90°:

$$\begin{gathered} k_{\mathrm{eff},\hspace{0.17em}\mathrm{MF}}=\hspace{0.17em}\left[\begin{array}{l}{k}_1\hspace{1em}\hspace{1em}\hspace{1em}0\\ 0\hspace{1em}\hspace{1em}\hspace{1em}{k}_2\end{array}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{and} \kappa =\hspace{0.17em}\left[\begin{array}{l}{\kappa }_1\hspace{1em}\hspace{1em}\hspace{1em}0\\ 0\hspace{1em}\hspace{1em}\hspace{1em}{\kappa }_2\end{array}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{when}\theta =0° \\ {k}_{\mathrm{eff},\hspace{0.17em}\mathrm{MF}}=\hspace{0.17em}\left[\begin{array}{l}{k}_2\hspace{1em}\hspace{1em}\hspace{1em}0\\ 0\hspace{1em}\hspace{1em}\hspace{1em}{k}_1\end{array}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{and} \kappa =\hspace{0.17em}\left[\begin{array}{l}{\kappa }_2\hspace{1em}\hspace{1em}\hspace{1em}0\\ 0\hspace{1em}\hspace{1em}\hspace{1em}{\kappa }_1\end{array}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{when}\theta =90° \end{gathered}$$

(24)

where k₁ > k₂ and κ₁ < κ₂. Thus, for a case with θ = 0°, there is a better thermal conductivity in the r direction to transfer the heat from the HTF tube wall into the enclosure. Moreover, there would be a high permeability in the z-direction, allowing the molten liquid to move in an upward direction and disperse the heat by the advection mechanism. In contrast, a case with θ = 90° provides minimal thermal conductivity in the r-direction and minimal permeability in the z-direction, which declines the advection and conduction mechanisms in the molten PCM. The same conclusion is proper for negative angles.

The stored energy shows the same trend as MVF since the thermal energy storage is dominated by latent heat energy storage. The sensible energy is also most significant in the liquid region, closely following the MVF. Thus, the highest energy storage rate can be seen for an anisotropic angle θ = 0 °C.

Figure 10 shows the temperature distribution and the interface for phase change heat transfer in the unit when the anisotropic angle is 0°. The temperature distribution is reported for the PCM (T_PCM) and MF (T_MF) phases. The temperature labels have not been added to avoid congestion. The contour levels are 2°, and the minimum and maximum temperature of each image were also reported above the image. The time of each image is also indicated below the image. As seen, the temperature differences between MF and PCM phases are not much. The differences are mostly next to the insulated wall. There is no temperature difference at the HTF tube wall, indicating accurate implementation of the continuity of temperatures at the boundary. The difference between the PCM and MF at the insulated wall is due to MF’s superior thermal conductivity, which better channels the heat into the solid regions. There are also slight temperature differences at the bottom due to the natural convection effects. They are significant at the final times, with a notable amount of liquid PCM and convection circulation flows.

The melting interface starts from the bottom, where the HTF enters the tube at the highest average temperature. Moreover, the convection heat transfer between the HTF flow and tube wall is highest in the entrance region. The heat reaches the PCM, and melting commences. Thus, the melting front develops from the bottom and extends toward the top right. Due to the presence of MF in the unit, the natural convection circulation is weak, and the melting front mainly progresses from the bottom.

Figure 11 shows the isotherms and melting front for an isotropic case with no anisotropic properties. The temperature difference between T_PCM and T_MF is almost similar to the anisotropic case of θ = 0°, which was addressed in Figure 10. The difference between isotherm lines near the bottom is more pronounced. The MVF maps also show that the interface moves toward the top right of the enclosure more slowly compared to the anisotropic case of θ = 0°. A comparison between the melting interface of Figure 10 and Figure 11 at t = 4 h reveals a more solid area for the isotropic case (Figure 11).

Figure 12 shows the temperature and melting maps for an isotropic case with θ = 90°. The difference between T_PCM and T_MF lines is almost similar to previous images. However, the difference is more inclined to the vertical direction, where the thermal conductivity of MF is slightly higher in a vertical direction for this case.

Figure 13 plots the impact of anisotropic angle (θ) and inlet pressure (P_in) on the MVF, energy storage power, and stored energy when the melting fraction reaches 0.9. The value of MVF = 0.9 was selected for analysis since, at this number, most of the enclosure is in a liquid state, representing the unit’s thermal behavior next to a full charge level. When the enclosure reaches a fully charged level (MVF = 0.999), there is no solid PCM left, and the unit behavior is mainly dominated by a small variation of sensible heat. Thus, an appropriate state next to the fully charged state, i.e., MVF = 0.9, was adopted to investigate the impact of control parameters on the thermal behavior of the unit.

Figure 13a shows that the increase in inlet pressure notably reduces the charging time. Regardless of the inlet pressure, the charging time is minimal for an anisotropic angle θ = 0°. This finding agrees with the findings of Figure 9 for a fully melted unit when P_in = 5 Pa. For the optimum anisotropic angle θ = 0°, increasing the inlet pressure from 5 Pa to 7.5 Pa (a 50% increase) reduces the melting time from 273 min to 247 min. Thus a 50% increase in the inlet pressure could save the charging time by 9.5%. For an inlet pressure of 10 Pa, the optimum melting time was 230 (min), which is only 15.7% lower than that of P_in = 5 Pa. The Reynolds number for the inlet pressures of 5 Pa, 7.5 Pa, and 10 Pa was 1054, 1628, and 2233, respectively.

Moreover, a semi-symmetric behavior for negative and positive angles can be observed. The maximum charging time corresponds to anisotropic angles +90° and −90°. These trends agree with the results of Figure 9 for a fully charged enclosure.

The energy storage power (Figure 13b) shows a similar trend of behavior to the charging time but in reverse. Figure 13c shows the amount of stored energy. There are small energy fluctuations (<1%) for various inlet pressure and anisotropic angles when the enclosure reaches a charging state with MVF = 90%. These fluctuations are primarily due to the sensible heat energy, a function of the unit’s local temperature distribution. Most of the enclosure is in a liquid state, with slight temperature variations in the liquid phase. Thus, since there are slight variations in the amount of stored energy, the behavior of the energy storage power is almost proportional to the inverse of melting time (Figure 13b). The maximum energy storage power for the unit is 119.9 W, obtained for P_in = 10 Pa and θ = 0°. The reference lines for an enclosure filled by an isotropic metal foam are also plotted in all sub-figures of Figure 13. Figure 13a,b show that an anisotropic angle below 45° (|θ| < 45°) could reduce the melting time and increase the charging power compared to an isotropic metal foam design.

The charging time (MVF = 0.9) for the isotropic metal foam was computed as 292 min, 266 min, and 249 min for inlet pressures of 5 Pa, 7.5 Pa, and 10 Pa. An optimum anisotropic design (θ = 0°) could reduce the charging time to 273 min (6.5%), 247 min (7.1%), and 230 min (7.6%). Considering all inlet pressures, the improvements average is about 7%. The same improvement can be seen for energy storage power. Thus, compared to an isotropic metal foam design, an anisotropic property (ζ = 0.2) could improve the unit’s thermal performance by about 7% without impacting the unit’s thermal energy storage capacity or adding weight.

5. Artificial Intelligence and Optimization

A neural network (NN) is a set of algorithms that attempt to discover fundamental relationships in a group of data through a process that mimics how the human brain works. Thus, a NN can be explained as a system of organic or artificial neurons. A “neuron” in a NN is a function introduced mathematically, which renders information following its specific architecture. Thus, such NN is analogous to statistical methods such as regression and curve fitting [46]. Since the neurons in a mathematical neural network are introduced artificially, the resulting network is an Artificial Neural Network (ANN).

NNs are capable of adapting to variable inputs and thus providing accurate outputs. A NN contains layers of interconnected nodes (perceptron) and activation functions. A feed-forward neural network, which was utilized in the present study, sends information in one direction through its entry node and processes them in a unidirectional fashion until they reach the output node [46].

The computational time for each simulation was reported in Table 2. It was observed that the computational time for the numerical experiments was more than one hour. This means, each simulation for a set of the design parameters, i.e., P_in and θ, takes more than an hour. Therefore, finding an accurate numerical relationship between the design parameters and the melting behavior (MVF) was computationally expensive. Here, the ANN is used to find an accurate map between the design parameters and the design target at a low cost.

In order to train the ANN, the melting process data for cases 1–39, consisting of 4998 sample data, were fed to a neural network to learn the relationship between the input control variables of inlet pressure (P_in), anisotropic angle (θ), and melting time (t), and the target data (MVF). The melting instances were captured from simulation data at each 2.5 min instance until a full melting reaches. Then, the ANN was trained and used to learn the relationship between the input and target data. The size of the training database was 4998 samples which were divided into 70% training, 15% test, and 15% validation.

Table 4. A summary of the input and output variables for the ANN model and the data partitions.

Input Variables			Output Variable	Train	Test	Validation
Pin (Pa)	θ (deg)	t (s)	MVF	70%	15%	15%

shows the list of input and output variables.

The Levenberg–Marquardt optimization, which is a network training function that updates weight and bias values, was used. This method is one of the fastest supervised algorithms, which uses the backpropagation algorithms for training. Training record (epoch and performance) returned as a struct whose fields depend on the network training function. A trained neural network structure contains the following parts: (1) Training, data splitting, and performance functions and parameters; (2) Data split index for training, validation, and test sets; (3) Data split mask for training validation and test sets; (4) Number of epochs and best epoch; (5) A list of training state names (states); (6) A field for each state name that records the value during training; and (7) Best network performance such as train, validation, test and best.

Some studies, such as [47,48], utilized systematic a hyper-parameter method to select the ANN details. However, in the present study, some preliminary investigations on the number of neurons and the number of hidden layers were performed manually. It was noted that using two hidden layers with many neurons results in excellent training data, but results in an overfit and deteriorates the validation predictions. Thus, only ten neurons with a single hidden layer were used to avoid overfitting issues. Since the utilized neural network was fairly simple and could well predict the physical model no further investigations were seem necessary for fine-tuning the ANN parameters.

Figure 14 shows the structure of the utilized neural network, which consists of ten neurons in its first hidden layer and an output layer made of one neuron. Figure 15 depicts the neural network prediction accuracy for the test, training, and validation aspects. It should be noted that the validation data should not be just used for validation [49] and not be fed to the ANN during the training or other processes. Here, the validation data were independent of the training process, but it was just plotted in the same figure for sake of representation. As seen, the ANN learns to predict the results after a few epochs very well. The validation process consists of results for cases that the ANN has never seen during the training and test process. Thus, this index can show the reliability of the ANN in estimating the solution field. The best validation was achieved after 58 epochs; hence, this state was adopted as the best state of ANN for producing the maps of MVF as a function of control parameters. Figure 16 contains the subplots for the train, test, and validation data details. The y-axis shows the normalized estimated output while the x-axis shows the amount of normalized target (CFD) value for each sample. This figure shows that the ANN estimated most of the data with excellent accuracy.

The trained ANN was used to produce maps of MVF as a function of control parameters. Figure 17 depicts the contours of MVF for various anisotropic angles and inlet pressures. The MVF was simulated using the trained ANN. The inlet pressure was adopted in the range of 5 Pa–10 Pa with intervals of 1 Pa, while the anisotropic angle was in the range of −90° to +90° with an interval of 10°. The plots were created for hourly melting intervals of 2 to 5 h. Based on the previous data, an anisotropic angle of 45° provides similar data to a uniform metal foam with no anisotropic properties. A zero anisotropic angle can provide the melting volume fraction at a low inlet pressure comparable with a much higher inlet pressure but with anisotropic angles about 45° or higher. Thus, an optimum anisotropic angle can promote MVF without adding weight to an LHTES unit.

6. Conclusions

The impact of utilizing an anisotropic metal foam on the melting behavior of a PCM was addressed in a shell-tube latent heat storage unit. The permeability and thermal conductivity were a function of anisotropic intensity and angle. A two-equation heat transfer model, based on the local thermal non-equilibrium approach LTNE model, was employed to capture the metal foam and PCM temperatures. About 40 numerical experiments were executed, and the impact of the anisotropy angle (between −90° and 90° with an anisotropy factor ζ = 0.2) and the HTF inlet pressures on the melting behavior were examined. Results are compared to a reference case of isotropic permeability and conductivity. An artificial neural network was utilized to learn the transient behavior of the melting heat transfer in the LHTES unit. Then, the neural network model was used to discuss the critical impact of control parameters on the melting rate.

Results showed a decrease in the anisotropy angle from $\left|\theta \right|$= 90° to 0° accelerates the melting rate. Nevertheless, the typical isotropic MF (reference case) can produce a melting process between the anisotropic MF when the anisotropic angle is |θ| = 45°. However, increasing the anisotropic angle above 45° to 90° slows the melting process.

The increase in the inlet pressure of the HTF from 5 Pa to 10 Pa accelerates the melting process by about 16% for the optimum angle of anisotropy θ = 0° and increases the storage capacity. Hence the maximum storage power is achieved for θ = 0° and P_in = 10 Pa.

The feed-forward neural network successfully learned the transient behavior of MVF from 4998 sample data with an R-value greater than 0.99. Using the ANN, hourly maps of MVF as a function of pressure and anisotropic angles were produced. Using the MVF maps, the advantage of anisotropic angle versus pressure drop was discussed. During the whole melting process, using a zero-angle anisotropic MF could significantly accelerate the melting process compared to a uniform MF.