Numerical investigation on latent thermal energy storage in shell and corrugated internal tube with PCM and metal foam

A numerical investigation on Latent Heat Thermal Energy Storage System (LHTESS) based on an aluminum foam totally filled with phase change material (PCM) is accomplished. The PCM used is a pure paraffin wax with melting over a range of temperature and a high latent heat of fusion. The LHTESS geometry under investigation is a vertical shell and tube. The corrugated internal surface of the hollow cylinder is assumed at a constant temperature above the PCM melting temperature. The other external surfaces are assumed adiabatic. The paraffin wax phase change process is modelled with the enthalpy-porosity theory, while the metal foam is considered as a porous media obeying to the DarcyForchheimer law. Local thermal non-equilibrium (LTNE) model is assumed to analyze the heat transfer in the metal foam. The governing equations are solved employing the Ansys-Fluent code. The numerical simulations results, reported as a function of time, and concerning the LHTESS charging phase, are compared in terms of melting time, average temperature and energy storage rate. The corrugated internal surface effect is analyzed with respect to the wavelength and wave amplitude of the corrugation.


Introduction
Greenhouse gases emissions have significantly increased in the last decades, especially owing to the widespread exploitation of fossil fuels. The main solution to slow down the mentioned phenomenon is definitely the use of renewable resources, whose availability is however discontinuous. Storage systems are certainly a remarkable tool to face the time mismatch between energy demand and supply. Specifically, Thermal Energy Storage (TES) units have been widely studied and numerous researches give evidence of their benefits in terms of heat handling and saving both for civil and industrial applications. Phase change materials (PCMs) are characterized by high melting latent heat, which allows them to embed large amounts of thermal energy in reduced volumes. Paraffin wax is finding a large use in Latent Heat Thermal Energy Storage System (LHTESS), since its melting temperature range makes it a useful heat transfer fluid (HTF) for several applications. The main weakness of this material is its low thermal conductivity, which negatively affects the system heat storage rate. In order to overcome this flaw, several solutions have been investigated and proposed. Joshi and Rathod [1] performed numerical simulations based on Local Thermal Non-Equilibrium (LTNE) method to evaluate the effects of a copper foam on the charging phase of a LHTESS. Their study highlights a remarkable thermal transient reduction can be obtained by increasing the fill height ratio of the metal foam and decreasing the volume porosity. Chen et al [2] studied a cylindrical TESS containing a paraffin wax of Rubitherm, RT 58, embedded in a copper foam in the annular section, as hot water flowing through the internal one is used as HTF. Their research not only shows that the melting time is shortened by 10 times by totally filling the annular section with the metal foam respect to the clean configuration (only PCM), but its application in the HTF side enhances the heat transfer, so reducing additionally the melting time and the thermal transient. Moreover, their simulations prove that initial hot fluid temperature strongly affects the heating phase. Sardari et al [3] investigated the thermal efficiency of a storage system made up of RT 35 embedded in a high-conductive copper foam. Their study gives evidence that porosity and the heat source position considerably affect the PCM melting and the charging time. Specifically, when thermal energy is provided to the fluid by the cross-sectional faces, the thermal transient significantly decreases owing to the heat transfer area augmentation. Furthermore, a porosity reduction leads to a system overall diffusivity increase, so determining a higher storage rate. Nie et al [4] employed the enthalpy-porosity model to perform a numerical study in LTNE assumption, in which both charging and discharging processes have been evaluated for three different shell and tube configurations of a TESS heated with water: cylindrical, with conical shell and with frustum tube (the last two with variable tilt angles). The simulations show that non-cylindrical geometries guarantee an enhanced natural convection (and a more efficient conductive transfer in the case of frustum tube), culminating in a higher melting rate. Finally, their work reveals that the advantages brought by the alternative configurations are less significant when air is employed as HTF instead of water. Zhang et al [5] performed a 3D numerical study on a plane geometry, focusing on the evaluation of the melting rate during the heat transfer at variable operating conditions. Their study shows that higher Rayleigh numbers produce an increased buoyancy effect especially in the upper side occupied by the PCM, so generating different temperature distributions. A RT 82 composed with copper foam storage system employed in a ventilation unit for dwelling applications has been studied by Sardari et al [6]. Both 2D and 3D simulations have been carried out to assess the PCM solidification and the simultaneous heating phase of an air stream of pre-determined inlet temperature and mass flow rate. The conductivity enhancement due to the presence of the copper foam results in a significantly higher PCM solidification rate and a more efficient heat transfer. One more research by Buonomo et al [7] concerns a totally filled LHTESS consisting in paraffin wax and high-conductive aluminum foam. The enthalpy-porosity model in LTNE assumptions was employed to numerically investigate the time evolution of the average liquid fraction of a vertical shell and tube cylindrical storage unit with the inner surface as heated wall. The totally filled configuration allows to reduce by two orders of magnitude the melting rate compared to the clean case. Moreover, the thermal diffusivity increase generates a uniform cross temperature distribution through the entire simulation. The present study proposes to numerically investigate the charging phase of a LHTESS made up of aluminum foam totally filled with RT 50. The main indicators describing the process have been evaluated under varying geometric and structural configurations of the system.

Geometrical and physical models
The LHTESS under investigation is a vertical shell and tube unit with a corrugated internal cross section. The space between the external shell and the corrugated internal tube contains the porous matrix, consisting in aluminum foam, which is totally filled with paraffin wax RT50. The 3D axisymmetric geometrical configuration is reported in Figure 1. Figure 2 outlines the 2D configuration of the system, whose main geometric parameters are listed below: H is the cylinder height, equal to 100 mm, ri,avg and re correspond, to the average internal radius and the external one, and measure, respectively, 6 mm and 50 mm. In addition, in this investigation, an aluminum foam featured by a pore density of 20 PPI has been considered. The inner surface, supposed to be the heat source, has been modeled as a cosine wave, having the following analytic expression: where r is the radial coordinate, A is the wave half-amplitude and n stands as the waves number.
In the basic geometric configuration, the parameters identifying the cosine wave have been given the values below: A=2 mm, n=6. Furthermore, the previously mentioned basic configuration includes an aluminum foam with porosity of 0.9546. In this study, different foam structures as well as different corrugated wall geometric parameters have been analyzed. To this proposal, different porosities, wavelengths and waves amplitudes have been considered. The thermophysical properties are considered temperature independent and are summarized in Table 1 [7,8]. 1) An unsteady and laminar two-dimensional axisymmetric flow is considered.
2) The internal cylindrical surface represents the heat source and is supposed to be at the uniform temperature of 350 K, whereas the external and cross-sectional ones are assumed adiabatic.
3) At the beginning of the thermal transient, the temperature is equal to 300 K. 4) The Boussinesq approximation is used to simulate the PCM natural convection due to the buoyancy phenomenon, which is considered to trigger when the temperature exceeds the RT50 solidus temperature. 5) The PCM thermal properties are assumed constant and evaluated in the liquid phase, at 300 K; only its density is assumed variable with temperature, just according to the Boussinesq approximation, previously mentioned. 6) The metal foam is isotropic and homogenous and simulated as a porous medium.
Furthermore, the LTNE hypothesis and the Forchheimer-extended Darcy model have been adopted to model the heat transfer between the paraffin and the aluminum foam, and to simulate the presence of the metal foam. As a consequence, the equations governing the physical model are the following. Continuity equation [7], with u and w being, respectively, the radial and axial components of the fluid velocity in the cylindrical coordinates r and z.
Momentum equations [7], Where ρp and μp stand for the paraffin density and dynamic viscosity, p represents the static pressure, ε is the foam porosity. The buoyancy effect is taken into account through the last term in equation (4), in which γ is the volumetric expansion coefficient and T0 is a reference temperature, corresponding to the RT 50 solidus temperature, equals to 318 K. The Darcy and Forchheimer terms, both contained on the right-hand side of the momentum equations simulate the pressure losses, respectively, due to viscous and inertial effects related to the coexistence of the PCM and the metal foam [9]. The variable K indicates the foam permeability, which has been calculated through the following equation [9]: Furthermore, the analytic expression of the variable CF, known as inertial drag factor, is given below [9]: The parameters governing the Darcy and Forchheimer terms depend on df and dp, respectively being the metal foam ligament diameter and the pore diameter. These parameters describing the foam structure are related through the following ratio [9]: In the simulations performed in this study, dp has been assigned the experimental values provided in references [9,10]. In equations (3) and (4), the term Si is associated both to the presence of the solid matrix in the mixed zone and to the enthalpy-porosity model, and can be expressed as follows [11]: where β represents the local liquid fraction, 0.001 is added to prevent the division to zero when β is next to zero, Amush is a constant which governs the amplitude of the velocity dumping owing to phase change processes. A fixed Amush value equal to 5x10 5 [kg m -3 s -1 ] has been set in the present investigation. The volume liquid fraction β, previously mentioned, is an index measuring the ratio between the liquid volume and the PCM total volume. Hence, it can take any value between 0 and 1, according to the local temperature T, as outlined below [ Under LTNE assumption, PCM and metal foam are featured by distinct temperature fields, as a consequence, two energy equations have to be solved.
The next-to-last term on the right-hand side of equation (10) is the unsteady term associated with the phase change; indeed, it's proportional to the PCM melting latent heat HL. In this investigation, the enthalpy-porosity model has been used to describe the physical correlation between melting process and heat storage. Generally, the overall specific enthalpy content of a material involved in a melting process can be expressed through the formula below [12]: H h H    (12) As summarized below, h is the sensible contribute per mass unit, intended at a specific temperature T, whereas ΔH is the latent enthalpy content, proportional to HL through the volume liquid fraction β [12]: Furthermore, keff,p and keff,s correspond, respectively, to the PCM and metal thermal conductivity. These properties have been assessed using the formulations below [14]: where kRT50 and kAl are the RT 50 and aluminum conductivity at 300 K (provided in Table  1). Furthermore, the contribution of the dispersion thermal conductivity kd should be considered in such an investigation. This coefficient linearly increases with velocity, as highlighted by the correlation given by Calmidi and Mahajan [9]. Since the velocity field is featured by a very low intensity in this case study, according to the previous formulation, kd is by two orders of magnitude lower than kRT50, hence it's been neglected. The right-hand sides of both energy equations contain a term associated to the convective heat transfer between MF and PCM. Several correlations are available in literature to determine the convective heat transfer coefficient hsf. In the present investigation, the one given by Churchill, widely used in natural convection problems, has been employed [15]. This correlation has been implemented in an User-Defined Function in C++ language. Furthermore, the interfacial area asf per volume unit depends on the metal foam structure and has been calculated as it follows [9]:

NUMERICAL MODEL
The governing equations (2) -(4) and (11) - (12) have been solved by means of Ansys-Fluent commercial code. The calculation domain has been discretized according to the finite volume method. The SIMPLE algorithm has been employed for the pressure-velocity coupling, while PRESTO algorithm has been used for the pressure calculation. The high order term relaxation is set equal to 0.75 for all variables, to avoid the numerical instability of the system. The unsteady simulations have been performed with a time step size set to 0.2 s. The convergence errors have the following values: 10 −3 , for continuity and momentum equation and 10 −6 for energy equation. The numerical grid is made up of quadrilateral cells having smaller size near the boundary and a larger size in the center of the grid, to improve the solution accuracy on the boundaries of the geometry. Cell size is uniform along the axial direction, and variable along the radius with a maximum size ratio equal to 12. To check the accuracy of the results, a grid independency analysis has been performed with three discretization levels: 738, 2916 and 11431 cells. The grids have been compared in terms of time evolution of liquid fraction and average temperature. These evaluations have been carried out on the basic geometrical configuration of the system, already introduced in Paragraph 2. The discretization analysis has highlighted the 2916 cells grid represents the best compromise between solution accuracy and computing cost; hence, it has been employed to perform the simulations. Indeed, as shown in Figure 3c  Figure 3d, a modest difference in terms of both average temperature and liquid fraction (less than 1%) can be observed by the comparison between the 738 and 2916 cells grid during the thermal transient. Furthermore, a comparison between the Local Thermal Equilibrium (LTE) and LTNE assumptions has been carried out in the present investigation. Figure 4 just highlights that a good agreement exists between the two approaches given the negligible maximum difference between MF and PCM temperature and melting rate (respectively reported in Figure 4c and Figure 4d).

RESULTS AND DISCUSSION
In the present investigation, the basic configuration, introduced in the previous paragraph, has been compared to systems characterized by both different geometrical features (wave amplitude and wavelenght) and structural parameters (porosity). In order to highlight how these variables affect the simulations results, the following indicators, which exhaustively describe the charging phase of a LHTESS, have been assessed: average temperature (both for the PCM and MF zone), average liquid fraction, overall enthalpy stored, storage rate.

Porosity variations
The foam structure has remarkable effects on the paraffin melting process, since the higher is the MF volume fraction, the more efficient is the heating process. In order to highlight this correlation, three porosity values have been analyzed: ε=0.9005, ε=0.9546 (basic configuration), ε=0.978. The time trends of the indicators under examination, showed in Figure 5, highlight that the presence of a high conductive metal foam improves the heat transfer inside the computational domain, in spite of a low increase of the overall thermal inertia. The quick melting phase in the lowest porosity foam (ε=0.9005) determines a greater latent energy storage rate, leading to a reduced thermal transient. It is worth specifying different values for cylinder height H have been assigned in each case, in order to take into account the PCM amount variation with different porosities. Specifically, as previously mentioned, for the basic configuration, corresponding to the intermediate porosity value (ε=0.9546), the cylinder height has been set equal to 100 mm. As a result, the PCM volume VPCM is calculated as VPCM = π (re 2 -ri,avg 2 ) H ε = 738.57 cm 3 . In order to prevent the PCM volume from being modified as the foam structure varies, the present value has been imposed for all the studied systems. As a consequence, according to the evaluation of the volume, a cylinder height equals to 106 mm has been assigned to the storage unit featured by the lowest porosity (ε=0.9005); on the contrary, for the most porous foam (ε=0.978), the investigation has been carried out with H=97.6 mm.

Wave amplitude variations
Several simulations have been performed to investigate the influence of the cosine wave amplitude on the thermal transient. The basic case, having A=2, has been compared to the following configurations: A=0 (not corrugated channel), A=4. The results are provided in Figure 6. A noticeable heat storage rate improvement occurs for greater amplitudes, owing to the increase of the inner surface area.

Waves number variations
The last parameter whose influence on the charging phase of the system under investigation has been analyzed is the wavelenght. To this purpose, in the present investigation the TESS performances have been assessed with the following waves numbers: n=0 (not corrugated channel), n=3, n=6 (basic configuration), n=9, n=12. This parameter increase leads to a heat transfer surface augmentation; as a consequence, as shown in Figure 7, greater waves numbers generate a higher PCM melting rate, which results in a more efficient charging time for the TESS.

CONCLUSIONS
The study aims to assess the influence of both structural and geometric parameters on the charging phase of a LHTESS consisting of a high-conductive aluminum foam totally filled with RT 50. The investigation has been carried out with a constant pore density equal to 20 PPI. The heat source is represented by the wavy inner surface of the cylinder, whose initial temperature is supposed to be higher than the PCM melting temperature range. As several experimental and numerical researches have given evidence of the foam porosity effect on the melting phase of paraffin wax in smooth storage units, few studies have been lead dealing with corrugated configurations. The present investigation proves, by increasing the waves amplitude and reducing the wavelength, the phase change process is enhanced and, consequently, the same does the storage rate.