Numerical Study of Melting Process of a High-Temperature Phase Change Material Including Natural Convection and Turbulence

Thermal energy storage (TES) systems including sensible heat (SHTES) and latent heat (LHTES) can be integrated in concentrated solar power (CSP) plants, delivering continuous and cheap electricity. SHTES systems have been studied and integrated into CSP plants for sometimes as reported by Kuravi et al. [1]. However, there have been less research efforts in the field of LHTES system, which is considered to deliver higher energy storage density compared to SHTES, resulting in a smaller scale storage system as discussed by Wang et al. [2] which is also cost effective. With SHTES, heat is stored or released due to the heat capacity ( $\dot{\mathrm{mC}} \mathrm{C}_{\mathrm{p}}$ ) of the storage medium (fluid) whereas a LHTES unit stores heat via melting a PCM and releases heat through freezing the PCM during a nearly isothermal process. In both systems, natural convection heat transfer is a major heat transfer mechanism. In a LHTES system, the thermal conductivity of potential PCMs are low as reported in a review study by Kenisarin [3], and the natural convection role is important to enhance the heat transfer rate. In most cases natural convection is laminar, particularly in low temperature PCMs (e.g. $\mathrm{T}<50^{\circ} \mathrm{C}$ ) resulting in a process with $\mathrm{Ra}<10^9$. However, in high-temperature PCMs (e.g. $300^{\circ} \mathrm{C}<\mathrm{T}<700^{\circ} \mathrm{C}$ ) with higher temperature differences between a heat source (e.g., a heat transfer fluid (HTF)) and the melting temperature of PCM ( $\Delta \mathrm{T}=\mathrm{T}_{\mathrm{h}}-\mathrm{T}_{\mathrm{m}}$ ), the flow becomes turbulent for Ra~10 ${ }^{12}$ as reported by Bejan [4]. Turbulence makes the natural convection more effective in mixing the storage medium.

High-temperature PCMs are promising heat storage materials for LHTES systems to be used in CSP plants. High-temperature differences (e.g., $\Delta \mathrm{T}>50^{\circ} \mathrm{C}$ ) are expected in using a LHTES system, particularly for protection purposes when a fast response is required, e.g. protecting the receiver from high temperatures as proposed by Verdier et al. [5].

This study aims to investigate the impact of natural convection and turbulence in a numerical study of the melting process of a high-temperature PCM. To provide an insight to the subject, a brief review of the relevant research work is reported here.

In a two dimensional numerical study, Markatos and Periclbous [6] developed a model to capture the changes in temperature and velocity fields while flow changes from laminar to turbulent in an enclosed cavity, i.e., from $\mathrm{Ra}=10^3$ to $10^{16}$.

The results of the study revealed higher Nu and heat transfer rate for higher Ra for which correlation between Nu and Ra were proposed.

Basal and Ünal [7] developed a computational code to study a triple concentric-tube LHTES, ignoring natural convection. Results showed that melting duration does not change with the mass flow rate of HTF, however, it decreases exponentially as the temperature of the HTF rises. Including natural convection calculation might change the result of this study, which was not addressed by the authors.

In a numerical investigation, Belusko et al. [8] studied the impact of different shell and tube configurations, e.g. parallel and counter flow on the charging and discharging process of a LHTES system while ignoring natural convection. They found that the high heat transfer area for a longer period of the discharging process occurs in the counter flow configuration due to a one dimensional phase change front. This study could be conclusive for the discharge process as natural convection was not the dominant heat transfer mode, however, this was not the case for the charge processes due to the major role of natural convection.

In numerical and analytical work, Lakeh et al. [9] studied the impact of turbulent natural convection in a supercritical TES in a vertical tube with fixed outer temperature. Comparing different aspect ratios of the system, they found that the charging time is significantly shorter due to the heat transfer enhancement effect of turbulence.

In other numerical work, Verdier et al. [5] studied the charging and discharging process in a rectangular LHTES system, with a high-temperature difference between the heat source and sink, e.g. $\Delta \mathrm{T}=200^{\circ} \mathrm{C}$. In that study, the results of the numerical calculation were not compared with the experimental data, hence the accuracy of the predicted data was not determined.

There are studies that natural convection was included or ignored, however limited research has investigated the natural convection impact on the mass transfer, heat transfer area (melt front in a melting process), the temperature gradient in space and time for the melting process of a high-temperature PCM. Moreover, there is lack of investigations on the impact of turbulence on the melting process in a LHTES system and the importance of this phenomenon in high-temperature PCM systems. Due to the importance of natural convection and turbulence in the charging process of heat storage systems, this study aims to numerically investigate firstly the effect of including natural convection and secondly the impact of turbulent natural convection on the temperature and velocity fields, heat transfer area and rate, and duration of the phase change process.

The LHTES system assumed using sodium nitrate as the PCM and air as the HTF in a shell and tube configuration with parallel flow. In such a system, the HTF flows through tubes in parallel and melts/solidifies PCM which is confined in the shell. The melting temperature of the PCM, $\mathrm{T}_{\mathrm{m}}=306.8^{\circ} \mathrm{C}$, and all other thermophysical properties (e.g. $\operatorname{Pr}=9.2$ ) were acquired from a study by Lan and Yang [10]. Using this system, three melting processes were studied to compare the impact of convection and turbulence on the thermal behaviour of PCM.

Numerical modeling of the melting/solidification of a high-temperature PCM with natural convection was performed, using ANSYS Fluent [11]. The time dependent Navier-Stokes equations were solved including the latent heat transfer and phase change in the mushy zone, and zero velocity in the solid region as introduced in a previous work by Riahi et al. [12]. The following simplifications were assumed: (i) density is constant apart from the linear density-temperature relation in the Boussinesq approximation, (ii) same density for the solid and liquid phases, (iii) the natural convection flow is laminar (for the third modeling, turbulent flow was considered), incompressible and inviscid flow, (iv) internal thermal radiation is negligible as the Boltzmann number $\gg 1$, i.e. $\sim 10^4$, as suggested by Modest [13].

The enthalpy-porosity method by Voller et al. [14] with a fixed grid was used in fluent [11] to track the phase front by the evolution of temperature-enthalpy profiles. The mushy region where the phase fraction is between 0 to 1, was treated as a porous medium in which the porosity in each cell follows the phase fraction during a phase change process. The total enthalpy changes were equal to the latent heat of the fraction of the PCM changing phase, e.g. from solid to liquid. The melt fraction is defined to be 1 for the temperatures above Tm and 0 for the lower temperatures. A sink term is included in the momentum equation to keep solid velocity at zero. Details can be found in [11].

A model was developed in the previous work by Riahi et al. [12], using the solidification and melting option in Fluent [11] and was validated with experimental data from a study by Jones et al. [15]. In this work, the same model was used for the investigation of convection and turbulence impact on the temperature and velocity field, heat flux and melting process duration. The pressure based coupled algorithm was selected to solve the momentum and continuity equations. The second order upwind scheme for the advection term, the central differencing for the diffusion term and the second order implicit as a discretization scheme for the transient term were used. As a pressure interpolation scheme for the transient calculations, PRESTO was implemented. Realisable k-ε model with the standard wall functions as the near-wall treatment method were implemented for the turbulence calculations. The convergence threshold for all calculations were at least 10⁻⁴ for the continuity equation and 10⁻⁶ for other equations.

A dual-tube of a parallel tube bundle immersed in a shell and tube LHTES system was selected to study the thermal behaviour of PCM in a melting process under Ra $\sim 10^{12}$ where turbulent natural convection is involved. The specification of this model was acquired from a lab scale LHTES prototype which was built for the test purposes. The height of the tubes and PCM was 55 cm and the width of the geometry (centre to centre of tubes) was 6 cm . Three domains of HTF, tube thickness and PCM were defined. Boundary conditions were set as symmetry on both sides of the geometry, adiabatic at the bottom and top of PCM, flow inlet at HTF top and flow outlet at HTF bottom, as shown in Fig. 1. The HTF inlet temperature was set at $456.8^{\circ} \mathrm{C}$ resulting in $\mathrm{T}_{\mathrm{h}}-\mathrm{T}_{\mathrm{m}}=150^{\circ} \mathrm{C}$. The initial temperature was set $5^{\circ} \mathrm{C}$ colder than the PCM melting temperature ( $\mathrm{T}_{\mathrm{m}}=306.8^{\circ} \mathrm{C}$ ) to capture the latent heat transfer after the first stage of sensible heat transfer from the HTF to the PCM.

A two dimensional symmetric grid with 33,000 cells was generated from half of the geometry (Fig. 1) with a maximum aspect ratio = 2 and orthogonal quality = 1. A fine grid with 78,000 cells was generated to investigate the grid dependency of the results. The results of melting rate for the melting without turbulence with fine grid (cw-fg) is included in section 4.3 (Fig. 3b) and shows minor differences in the results of coarse grid for the same case. The time step was set to 0.1 seconds to capture details in the temperature and velocity fields in a feasible computation time.

Modeling of three melting processes with the same setting was performed, firstly calculating conduction heat transfer in PCM, then adding convection calculation, and finally including turbulent natural convection in the third melting process.

The temperature distribution in space and three time slots are presented in Fig. 2a–k. Figure 2a–c presents the temperature field of the melting process during pure conduction heat transfer while isotherms are mainly vertical. Comparing to Fig. 2d–f shows that involving natural convection results in horizontal isotherms with higher peak temperatures and temperature gradients. Due to the convection heat transfer at the highest levels of PCM, melting starts earlier during the process. However, involving turbulence (Fig. 2g–k) with a higher rate of mixing due to the higher rate of mass transfer, peak temperatures and temperature gradients are lower in comparison to the melting process including natural convection (Fig. 2d–f). In all melting processes, isotherms evolution being vertical in pure conduction and horizontal with convection and turbulence is consistent with the melting fronts evolution, as shown in Fig. 4a and b.

For the pure conduction melting process no mass transfer was predicted in the liquid PCM. The mass transfer rate for melting with convection and for melting with convection and turbulence is calculated by stream functions in the numerical modeling, however the profiles are not shown here due to space limitations. In general, involving turbulence results in a higher rate of mass transfer and increases by melting fraction, as shown in Table 1. In other word, higher liquid PCM results in more turbulent convection cells with higher velocity to carry and mix the fluid. Without turbulence, the mass transfer rate decreases during the melting process as the main driving force (difference between average temperatures of HTF and PCM) diminishes and slows down convection cells.

The heat transfer area being the melt front (Fig. 4a and b) and heat flux (Fig. 3a) can be compared to realise the impact of convection and turbulence on the melting process. In the case of pure conduction melting, the melt front is one dimensional and evolves parallel to the tube up to about 65% melt fraction and thereafter starts to shrink downward. In contrast, melt fronts in the two other cases with convection and turbulence are two dimensional, moving perpendicular to the tube after about 20% melt fraction. The two dimensional melt front before the 20% melt fraction means a higher heat transfer area with convection and turbulence, however, for most of the process (after 20% melt fraction) the melt front shrinks downward faster compared to pure conduction melting where the melt front holds higher heat transfer area between solid and liquid up to the end of the process.

For pure conduction melting, the heat fluxes are lower and decrease sharply during the process whereas including convection and turbulence presents more uniform processes due to the mixing effect of convection cells which facilitates a higher rate of heat transfer. Comparing turbulent convection with the case of ignoring turbulence (Fig. 3a), clearly shows a higher rate of heat flux with a more uniform process; providing a higher rate of heat transfer for a longer period and starts to decrease sharply at the later stages of the process. At about 90% melt fraction, pure conduction melting shows a higher heat transfer rate compared to the melting with convection due to the lower heat transfer area and also the lower velocity of convection cells which is the result of the lower driving force mentioned in 4.2. However, including turbulence compensates for the losses as a higher rate of turbulence (i.e. more convection cells with higher velocity) occurs at higher melt fraction (reported in section 4.2) resulting in higher heat transfer in comparison to the pure conduction melting and melting with convection.

Despite the lower heat flux for the case of pure conduction, the melting time is shorter in comparison to the melting with convection, e.g. 97% of PCM melts after about 5600s and 6500s, respectively. This is due to the higher heat transfer area of the one dimensional melting front in pure conduction which is effective at the later stages of process, e.g. after 60% melt fraction (Fig. 3b). Comparing the calculations with and without turbulence, it has been found that involving turbulence results in a shorter melting time (5600s in comparison to 6500s to melt 97% of PCM) due to the higher rate of convection heat transfer which increases with a higher rate of melt fraction and compensates for the lower average temperature difference between the heat source (HTF) and sink (PCM).

A numerical study of turbulent melting process for a high-temperature PCM was conducted, aiming to investigate the impact of natural convection and turbulence on the temperature and velocity fields, heat transfer area (melting front) and rate, and the melting time.

The results show that a one dimensional melting front, which moves parallel to the tube in pure conduction melting, provides higher heat transfer area. Involving convection and turbulence results in a two dimensional melt front which shrinks downward after 20% melt fraction and demonstrates a lower heat transfer area compared to the pure conduction melting. However, higher heat transfer rates are observed by including convection and turbulence due to mass and heat transfer by convection cells. Including turbulence results in a more uniform process due to the mixing effect of convection cells, resulting in lower peak temperatures and temperature gradients in comparison with melting with convection and ignoring turbulence.

Including natural convection and turbulence in numerical calculations provides a more realistic view of the charging process for the purposes of design and optimisation of a high-temperature LHTES.