Numerical Study on Phase Change Material with Metal Foam in Shell and Convergent/Divergent Tube Thermal Energy Storage Systems with External Heat Losses

The utilization of renewable energy has grown recently. By the end of 2021, renewable energy sources made about 38% of the capacity installed globally [1]. An increase in renewable energy also means an increase in the gap between energy supply and demand. The Thermal Energy Storage systems (TESs) must be improved to cover this gap. The TESs can essentially be separated into sensible, thermochemical, and latent heat. Sensible heat storage (SHS) is a method for storing thermal energy without undergoing any phase changes in the substance that is being employed. The energy is kept through material phase transitions with latent heat thermal energy storage (LHTES), whereas thermochemical storage (TCS) necessitates a reversible physio-chemical condition upon heating to store the thermal energy. The most promising of the three thermal energy storage techniques (TES) is LHTES because of its high energy storage density and isothermal phase transitions [2, 3]. Phase change materials (PCMs) are used in LHTESs because they can store thermal energy during phase changes from solid to liquid caused by temperature increases, and to release thermal energy during phase changes from liquid to solid caused by temperature decreases. Due to their distinct thermal characteristics, PCMs can nevertheless be effective in energy storage applications despite having limited thermal conductivity. Due to phase transition, PCMs have a high energy storage density, and the process of storing heat takes place in an almost isothermal operating environment. For these reasons, PCMs in LHTESs may store more energy per volume and are more practical than SHS materials. Additionally, the most crucial factor in choosing the optimum material for any application is the melting temperature of PCMs [4]. Organic PCMs can be used in residential applications including drying, heating, and hot water because of their low melting temperature less than 100℃, and they also offer safe working temperatures for research and development [5, 6]. A state-of-the-art review of SHS and LHTES on porous media was made by Mabrouk et al. [7]. They demonstrate some common methods for improving the PCMs' poor heat conductivity. As they demonstrated to deal with this issue, the suggested solutions include creating a composite by embedding the PCM in a porous matrix (metal foam or expanded graphite), adding metal spheres, fins, and wool to an existing PCM to create a new one with improved thermal conductivity, or incorporating some nanomaterials inside the PCM [7].

Truncated cone shell-and-tube TESs were proposed by Mao et al. [8]. Comparing the cylindrical and truncated cone TES tanks, they found that the truncated cone tank performed better than the conventional cylindrical tank under the same operating conditions and that the total melting time was reduced by about 30.69%. Alhusseny et al. [9] investigated a TES made of staggered coiled tubes filled with paraffin wax compounded to open-cell copper foam to address the issue of the PCMs' inadequate conductivity. A moderately hot or cold-water stream that is running across the tube-bundle units is used to charge or discharge the PCM unit. An enhanced LHTES was presented by Karimi et al. [10]. The energy storage storage container was separated into two parts, and these portions were in contact with a helical coil through which the heat transfer fluid (HTF) flowed. The coil's outside portion included pure PCM, while its interior region contained PCM embedded within open-cell metal foams. Some different configurations were studied in studies [11-13].

In their experimental study of a novel group of combined sensible-latent heat configurations for medium temperature, Zahid et al. [14] published their findings. The objective of the research was to suggest a more cost-effective and highly heat efficient TES technology.

The topology optimized fin, metal foam, longitudinal fin, and composite PCM were the four alternative improved designs that Li et al. [15] assessed and analyzed.

Shafiei Ghazani and Gholamzadeh [16] work examined the numerical impact of the conicality of either the shell or tube on the storage rate of a latent heat thermal energy storage tank. The study concluded that the melting time is not an appropriate criterion for evaluating the performance of a storage tank. Instead, a new parameter was introduced, which is defined as the time taken by the system to reach 99% of the maximum storable thermal energy. The study also found that, when the height of the tank is kept constant, the shape factor, defined as the cross-sectional area ratio of the PCM at the top to the bottom of the tank, is the most important determining parameter. The optimal shape factor in this study was around 10, and it can be achieved by increasing the tilting angle of either the shell or tube wall or a combination of both.

Nie et al. [17] have studied the effect the effect of geometry on the melting and solidification behavior of PCM on a vertical shell and tube configuration considering a conical shell or a frustum tube.

It appears that there is a lack of information about shell and tube TESs, specifically for PCMs embedded in metal foam having a variable internal section. Due to this, the shell and converging or diverging internal tube topologies with a PCM embedded in a metal foam are examined in this article. PCM and metal foam are used to fill the TES completely. Moreover, different heat transfer coefficients are assigned to the external surface of the TES to simulate a heat loss towards this surface. This work is novel in that it investigates the possibility of higher buoyancy effects, which could decrease the time the PCM takes to melt inside the tank, if the inner section is not constant. According to the writers' knowledge, this represents the first paper on this arrangement that also consider the PCM embedded in a metal foam and the heat losses toward the external surface of the shell and tube configuration.

The aim of this study is to investigate how the inclination of the interior tube affects the performance of thermal energy storage (TES) systems, taking into account the heat loss that occurs within the environment where the TES is installed. The study compares the effectiveness of two different setups: one with a vertical shell and tube configuration that has a convergent portion at its top, and the other with a top part that diverges. The inclination angle of the inside wall is varied between these two setups, with angles of 0°, 2.5°, 5°, and 10° selected to determine whether the melting time is shortened compared to the case without inclination. It is assumed that the heat transfer fluid flows from the bottom to the top, and external heat transfer coefficients of 0, 0.1, 1, 0.5, and 5 W/m²K are considered for each analyzed configuration. The study also assumes that the surrounding temperature is approximately 293 K.

In this section the mathematical model of the above physical problem is described. Considering all the assumptions, the governing equations have the following form starting from the Continuity equation [19]:

$\nabla u=0$     (1)

u denotes the velocity vector of the PCM in the liquid phase.

The following momentum equation is referred to the study [19]:

$\begin{align}  & \frac{{{\rho }_{p}}}{\varepsilon }\frac{\partial u}{\partial t}+\frac{{{\rho }_{p}}}{{{\varepsilon }^{2}}}u\cdot \nabla u=-\nabla p+\frac{{{\mu }_{p}}}{\varepsilon }{{\nabla }^{2}}u+ \\ & -\frac{{{(1-\beta )}^{2}}}{({{\beta }^{3}}+\omega )}{{A}_{mush}}u+ \\ & -\frac{{{\mu }_{p}}}{K}u-\frac{{{C}_{F}}}{\sqrt{K}}{{\rho }_{p}}u\left| u \right|+{{\rho }_{p}}g\gamma \Delta {{T}_{0}} \\\end{align}$     (2)

In this context, ρ_p represents the density of PCM and is dependent on temperature, assuming the Boussinesq approximation. μ_p is the viscosity of the PCM in its liquid phase, and p represents pressure. A_mush is a constant termed as the "mushy constant," which has a value of 1×10⁶ kg m⁻³ s⁻¹. This value was derived from Buonomo et al. [20] paper, which carried out an analysis on this constant. β represents the liquid fraction, and a constant ω is added to prevent division by zero when β is close to zero. Additionally, γ represents the thermal expansion coefficient of the PCM.

The PCM energy equation and the metal foam energy equation must both be solved under the LTNE assumption.

Under the LTNE assumption two energy equations must be solved, the first for the PCM and the second for the metal foam. These equations are [19]:

$\begin{align}  & \frac{{{\rho }_{p}}{{c}_{p}}\partial (\varepsilon {{T}_{p}})}{\partial t}+{{\rho }_{p}}{{c}_{p}}u\cdot \nabla {{T}_{p}}={{k}_{eff,p}}{{\nabla }^{2}}{{T}_{p}}+ \\ & -\varepsilon {{\rho }_{p}}{{H}_{L}}\frac{\partial \beta }{\partial t}+{{h}_{sf}}{{a}_{s}}_{f}({{T}_{s}}-{{T}_{f}}) \\\end{align}$     (3)

$\frac{{{\rho }_{s}}{{c}_{s}}\partial \,[(1-\varepsilon ){{T}_{s}}]}{\partial t}={{k}_{eff,s}}{{\nabla }^{2}}{{T}_{s}}-{{h}_{sf}}{{a}_{s}}_{f}({{T}_{s}}-{{T}_{f}})$     (4)

The energy equations incorporate several variables, including ε, which refers to the porosity of the metal foam, and ρ_p and c_p, which denote the density and specific heat of the PCM. The latent heat of the PCM is represented by H_L, while time is denoted by t. The effective thermal conductivity k_eff is determined by the conductivity of the metal foam (k_eff,s) and the conductivity of the PCM (k_eff,p). It can be defined using the following equations [21]:

${{k}_{eff,p}}=\varepsilon {{k}_{RT50}}+{{k}_{d}}$      (5)

${{k}_{eff,s}}=0.181\left[ {{(1-\varepsilon )}^{0}}^{.763} \right]\cdot {{k}_{Al}}$     (6)

The symbols k_RT50 and k_Al represent the conductivity of RT50 and aluminum at a temperature of 293 K. In this paper, the contribution of dispersion thermal conductivity (k_d), which increases linearly with velocity as per the correlation provided by Calmidi and Mahajan [22], is ignored due to the very low intensity of the velocity field. It is important to note that k_d has a thermal conductivity two orders of magnitude lower than k_RT50.

The term h_sf refers to the convective heat transfer coefficient, which represents the convective heat transfer between metal foam and PCM. This coefficient is obtained using the Churchill correlation [23], which is implemented with a User-Defined Function in C++, and then used in Ansys-Fluent. On the other hand, the term a_sf represents the surface area of the metal foam per unit volume and is calculated as per [22]:

${{a}_{sf}}=\frac{3\pi {{d}_{f}}(1-{{e}^{-((1-\varepsilon )/0.04)}})}{\begin{align}  & {{(0.59{{d}_{p}})}^{2}} \\\end{align}}$     (7)

In Eq. (2), the porous zone is modeled using the Darcy and Forchheimer laws. These terms account for the pressure losses caused by viscous and inertial effects resulting from the presence of PCM and metal foam, where the permeability of the porous media is denoted by K and the inertial drag factor is denoted by C_F. The value of K can be determined using the following equation [22]:

$\frac{K}{{{({{d}_{p}})}^{2}}}=~0.00073{{(1-\varepsilon )}^{-0.224}}{{\left( \frac{{{d}_{f}}}{{{d}_{p}}} \right)}^{-1.11}}$     (8)

The inertial drag factor, C_F, is equal to [22]:

${{C}_{F}}=0.00212{{(1-\varepsilon )}^{-0.132}}{{\left( \frac{{{d}_{f}}}{{{d}_{p}}} \right)}^{-}}^{1.63}$     (9)

Both parameters are influenced by the foam ligament diameter, d_f, and the pore diameter, d_p, which describe the structure of the foam. The two parameters are related to each other by the following ratio [22]:

$\frac{d_f}{d_p}=1.18 \sqrt{\frac{(1-\varepsilon)}{3 \pi}}\left(\frac{1}{1-(e)^{-((1-\varepsilon) / 0.04)}}\right)$     (10)

In this paper, the value of d_p is obtained from the studies of Calmidi and Mahajan [22], and Bhattacharya et al. [24]. The physical properties of PCM and metal foam used in numerical simulations are presented in Table 1 and Table 2, respectively. Table 1 summarizes the physical properties of aluminum and Paraffin wax RT 50 [25], assuming that the thermal conductivity and specific heat of Paraffin wax remain the same in both liquid and solid phases. Aluminum foam is selected as the porous medium due to its low cost, density, and weight compared to other materials. The parameters of aluminum foam, as reported in Table 2, are obtained following the work of Battacharya et al. [24], with a particular reference to sample 6. The momentum equation introduces a parameter known as the liquid fraction, β, which describes the mixed region of the enthalpy-porosity model and represents the ratio of the liquid volume to the total volume of a cell in the mixed region. The β value varies between 0 and 1 and depends on the local temperature in the cell:

$\left\{ \begin{align}  & \beta =0\quad \quad \quad \quad \quad \quad for\ \ \ T<{{T}_{solidus}} \\ & \beta =\frac{T-{{T}_{solidus}}}{{{T}_{liquidus}}-{{T}_{solidus}}}\quad for\ \ \ {{T}_{solidus}}<T<{{T}_{liquidus}} \\ & \beta =1\quad \quad \quad \quad \quad \quad for\ \ \ T>{{T}_{liquidus}} \\\end{align} \right.$      (11)

When the region is entirely solid, the value of β is 0, and vice versa, when the region is entirely liquid, β is 1. In Eq. (11), T represents the local temperature, T_solidus represents the temperature in the solid phase below which the domain is entirely solid, and T_liquidus represents the temperature in the liquid phase above which the domain is completely liquid. Melting occurs when the temperature falls between T_liquidus and T_solidus. The enthalpy-porosity method includes a source term in the momentum equations during the melting or solidification phases, accounting for the solid portion of the mixed zone. This solid portion is modeled as a pseudo-porous zone, with the liquid fraction representing the porosity.

Table 1. Thermal properties

Table 2. Foam parameters

Firstly, a comparison is made between the temperature of the fluid and solid in the porous zone. As shown in Figure 5, the two curves match well, indicating that the Local Thermal Equilibrium (LTE) hypothesis is a valid assumption in this case. Therefore, the corresponding term can be considered negligible in the energy equation.

5.png

Figure 5. Time evolution of solid and fluid average temperature

In Figure 6 is shown the average fluid temperature for different configurations and different inclinations angle using the same value of the convective heat transfer coefficient on the external surface. Figure 6 shows that the converging configuration has the shortest charging time. In fact, the temperature at thermal equilibrium is reached faster by converging configurations than by diverging configurations with the same tilt angle. Also, the higher the tilt angle, the faster the TES charges.

Figure 7 displays the time evolution of the average fluid temperature for the convergent configuration at different inclinations angle comparing the cases without and with heat losses. Figure 7 shows how in configurations where a convective heat transfer coefficient is assigned to the external surface, the temperature at thermal equilibrium decreases compared to the adiabatic case. In addition, since the volume varies as the tilt angle varies, it is seen that the temperature at thermal equilibrium is slightly different as the tilt angle varies for the cases with outward heat loss. For the adiabatic cases, the temperature at thermal equilibrium is the same.

6.png

Figure 6. Time evolution of the fluid temperature for different configurations at constant value of the convective heat transfer coefficient

7.png

Figure 7. Time evolution of the fluid temperature for the convergent configuration with different inclination angles and for different values of heat transfer coefficient on the external surface

Figure 8 shows the liquid fraction as a function of time for different geometric configurations as the tilt angle changes with a convective heat transfer coefficient value of 5 W/m²K. As previously mentioned, the higher the tilt angle the faster the PCM melts inside the TES. Likewise, the convergent configuration is always faster than the divergent configuration at the same tilt angle. This is due to the greater buoyancy effects present in the convergent configuration.

8.png

Figure 8. Time evolution of the average liquid fraction for different geometries

Figure 9 shows the liquid fraction as a function of time for the convergent configuration as the convective heat transfer coefficient changes. In this case, melting is faster in adiabatic configurations than in those with outward heat losses.

9.png

Figure 9. Time evolution of the average liquid fraction of the convergent configuration for different heat losses

Although the configuration with a tilt angle of 10° is the one with a larger volume, it is still faster than the others. In addition, in Figure 10, it can be seen that the configuration with a tilt angle of 10° is also the one with higher external heat losses. Obviously, the higher the convective heat transfer coefficient, the larger the value of heat loss to the outside.

10.png

Figure 10. Time evolution of the heat losses for different inclination angles and different convective heat transfer coefficient

Figure 11 shows the trend of outward heat loss as a function of time considering the same value of convective heat transfer coefficient. Converging and diverging configurations as the angle of inclination changes are shown in this figure.

11.png

Figure 11. Time evolution of the heat losses for different configurations and different inclination angles using the same heat transfer coefficient on the external surface

Converging configurations undergo a rapid increase in outward heat loss around 1200 seconds, this effect is surely due to the fact that the melting is totally completed in the upper part of the TES, this results in lower thermal resistance in that region. In addition, for such configurations the effects due to natural convection are more pronounced. Divergent configurations, on the other hand, have a more continuous pattern. Initially, it seems that the heat losses are greater for the diverging configurations, but instead, it is shown that the converging ones have the greatest losses. This is due to the fact that at the top of the TES the thermal resistances of the two configurations are different. For the divergent case, the top section is less than the converging case and therefore melts earlier initially decreasing its thermal resistance. While for the converging case that has larger top section than the diverging case, the such section melts slower, so the thermal resistance is initially higher. When such a section melts totally, the thermal resistance decreases very quickly.

Figures 12 and 13 show the specific enthalpy. It can be seen that by introducing heat losses to the outside, the energy stored by the thermal reservoir decreases. Also, as mentioned earlier, the higher the angle of inclination, the faster the maximum storable energy is reached.

12.png

Figure 12. Specific enthalpy for the convergent configuration with different inclination angles and heat transfer coefficient values

13.png

Figure 13. Time evolution of the specific enthalpy for different configurations with the same heat transfer coefficient values

The Authors would like to thank you very much Professor Enrico Lorenzini for the very useful discussion and suggestions on the paper.

This research was partially funded by MIUR (Ministero dell’Istruzione, dell’Università e della Ricerca), grant number PRIN-2017F7KZWS.