Comparison of enthalpy-porosity and lattice Boltzmann-phase field techniques for the simulation of the heat transfer and melting processes in LHTES devices

Thermal energy torage (TES) is a key enabling technology for the efficient exploitation of distributed generation systems based on renewable energy sources. Among the available options, research on latent heat TES (LHTES) solutions has been particularly active in the last decade, due to their ability to store and release high amounts of thermal energy in a very narrow temperature range. LHTES devices are based on phase change materials (PCMs), which act as thermal sinks or sources during their solid-to-liquid transition and vice-versa. As such, the development of reliable numerical tools for the prediction of the heat transfer and phase change characteristics of PCMs is of foremost importance, to help designing innovative and efficiently integrated LHTES implementations. In the present paper, the consolidated enthalpy-porosity (EP) method is compared to a novel lattice Boltzmann-phase field (LB-PF) algorithm in the simulation of a standard numerical benchmark for paraffin-like PCM melting problems. Performances and limitations of the two approaches are discussed, including the influence of model-related and purely numerical parameters. Outcomes from this study are used to confirm general guidelines for the application of well established methodologies, as well as to suggest new pathways for outof-standard modeling techniques.


Introduction
In 2019, the share of renewable energy sources (RES) in gross final energy consumption was at 19.7% in the EU-27, which is very close to the 2020 EU target of 20 % [1]. Compared with the 9.6% share in 2004, this might be regarded as a positive development, but the average value hinders very different performances at a national level. Sweden, for instance, had a 56.4% renewables share (with a 49% 2020 national target), compared to the 8.8% share of the Netherlands (with a 14% 2020 national target). Furthermore, a much higher effort is needed to increase the penetration of renewables and set the path towards the ambitious goal behind the European Green Deal -becoming the world's first climate-neutral continent by 2050 [2]. To this aim, the integration of distributed RES technologies with thermal energy storage (TES) systems has proven to be an effective solution to both reduce the energy consumption and mitigate the inherent RES local intermittency. For instance, studies from 2011 had already estimated an energy efficiency growth of 7.5% and CO 2 emission reduction of 5.5% if TES were to be deployed at a large scale in EU [3]. Recent examples of TES-based optimization of heating/cooling appliances for residential applications can be found in [4,5].
For low and moderately high temperature applications (below 300°C), latent heat TES (LHTES) technologies have attracted the attention of a continuously growing worldwide community of researchers in the last two decades [6]. All LHTES systems feature a Phase Changing Material (PCM) to accumulate/release large amounts of thermal energy in a narrow temperature range. Compared to sensible heat TES (SHTES), LHTES generally possess higher storage densities and allow for a better control of thermal fluctuations within the energy system. Typical limitations of the PCMs used in LHTES are [7]: low thermal conductivity (low charging/discharging power), subcooling, chemical stability and cost.
Concerning modeling of PCMs for TES applications, a typical issue is represented by the different physical behavior between the melting (charging) and solidification (discharging) stages [8]. More specifically, while solidification is conduction-dominated, the melting stage generates complex convective motions within the liquid phase which strongly affects the melting front position in time and thus the overall duration of the phase change process. As such, analytical solutions for the melting/solidification problem are available only for a few simplified configurations, thus increasing the need for developing and consolidating numerical tools as a key factor for the efficient LHTES design and optimization.
A common way to approach the numerical modeling of PCMs is a single-domain method called enthalpy-porosity (EP) technique [9][10][11], in which the melting front is not tracked explicitly. Instead, a scalar liquid fraction is associated with each computational cell or node in the domain, while a small semi-solid porous region called mushy zone is created at the solid/liquid interface as a function of the liquid fraction value. Although it has proven to be reliable for many PCM applications, the EP method performance is known to be significantly affected by the choice of two modeling parameters: the mushy zone constant and the phase change temperature interval (∆T ls = T liquidus − T solidus ). While the first parameter influence has been studied quite extensively for different PCM types [12], the interlink with the phase transition temperature range is currently less obvious.
As such, the first goal of the present work is to assess the combined effects of the mushy zone constant and the phase change temperature range for a reference 2D melting case [13,14], including also the impact of purely numerical parameters such as grid density. In addition to that, an original lattice-Boltzmann phase field (LB-PF) formulation [15,16] is proposed for the modeling and simulation of the same class of PCM melting problems. Lattice boltzmann algorithms for fluid dynamic simulations are known to be computationally efficient and easily modifiable/extendable to incorporate additional complex physical and chemical phenomena [17,18]. Therefore, reasons for considering an LB-based methods for melting/solidification modeling include: i) faster simulation times for a given set of thermopysical properties and initial and boundary conditions: ii) lesser dependence on model-specific parameters, to encompass a wider class of applications. To assess the proposed LB-PF methodology, the same 2D reference case has been simulated, highlighting differences, stengths and weaknesses compared to EP.
The remainder of the paper is organized as follows: first, the EP technique (as implemented in the ANSYS ® Fluent finite-volume software package) and the LB-PF method are briefly described; after that, the core section of the paper is dedicated to the reference 2D melting case, including a discussion on the results obtained with the two simulation approaches; finally, the main outcomes from the present study are summarizes in the Conclusions section.

Enthalpy-porosity method
The ANSYS Fluent CFD package features the standard EP technique for the modeling of melting/solidification processes [19]. The liquid fraction is computed at each iteration, based on an enthalpy balance, while the mushy zone is defined as a region where the liquid fraction lies between 0 and 1 and, at the same time, as a pseudo porous medium in which the porosity decreases from 1 to 0 as the material solidifies. When the material has fully solidified in a computational cell, the porosity becomes zero and hence the velocities also drop to zero.
The enthalpy of the material is computed as the sum of the sensible enthalpy h and the latent heat content ∆H, where: with c p = specific heat at constant pressure. The liquid fraction γ can be defined as: The latent heat content can now be written in terms of the latent heat of the material Λ: Finally, the energy equation is written as: with H = h + ∆H, ρ = density, k = thermal conductivity, ⃗ V = fluid velocity and S T = source term. The solution for temperature is essentially an iteration between the energy equation and the liquid fraction equation. Directly using Equation (2) to update γ usually results in poor convergence of the energy equation: to avoid this, in ANSYS Fluent the method suggested by Voller and Swaminathan is used to update the liquid fraction when T solidus T liquidus [11]. For pure substances (T solidus = T liquidus ), the original implementation based on specific heat is used instead [9].
The momentum sink due to the reduced porosity in the mushy zone takes the following form: where ϵ is a small positive number (default is set to 10 −3 ) to prevent division by zero, A mush is the mushy zone constant, ⃗ V p and is the solid velocity due to the (enetual) pulling of solidified material out of the domain.
The method described above is general and is able to account for thermal and solutal buoyancy, turbulent flows (for the liquid phase) and multicomponent mixtures.

Lattice Boltzmann -phase field method
To track the fluid dynamics and the thermal evolution within a melting/solidification problem, a standard 2-population LB approach can be implemented [17,18,20]: where τ F is the (standard) fluid relaxation time, while τ T is the thermal relaxation time.
Within this approach, the temperature is a passive scalar, whose effects are transferred towards the fluid according to specifically tuned forces. The three forcing terms at the rhs of Equations (6) and (7) represent: F B the buoyancy effect; F Drag an empirical mesoscopic drag, which enforces zero-velocity inside the solid region [15]; F L the latent heat release/absorption during melting/solidification phenomena, respectively. All these forces depend on the time evolution of the phase field ϕ that governs the melting/solidification front. The discretized rate equation for ϕ is: in which, R is the reaction term expressed as in [15]: where f + and f − are frequency scales (the inverse of the time scale that rules melting and solidification, respectively) and K + and K − are melting/solidification switch functions. Further details on the implementation of K + and K − , as well as on the forces dependence on ϕ can be found in [15] and [16].
3 Reference 2D cavity test case

Problem statement and case setup
The physics and scaling laws of 2D melting in a rectangular cavity with fixed temperature side walls have been discussed in dertails in the 1988 paper from Jany and Bejan [13]. The general problem is depicted in Figure 1 More specifically, a general Nu correlation can be written as follows: where: for any Pr. Note that θ = FoS te is a characterisctic nondimensional time coordinate, with Fo = αt/H 2 being the Fourier number (α = thermal diffusivity, t = physical time coordinate and H = cavity height). For a sufficiently high Rayleigh number (say above 10 5 -10 6 ), the heat transfer at the hot (left) wall is characterized by three subsequent regimes, as shown in Figure 1: i) an initial pure-conduction regime, with a steady decrease of Nu and a vertical wall-parallel melting front; ii) a mixed regime, with a knee appearing in the upper part of the melting front and the Nu value reaching a minimum and then increasing towards the pure-convection limit; iii) the pure-convection limit, with well developed thermal boundary layers on the left wall and melting front interface and with Nu approaching Nu ∞ . Three configurations have been considered for the purposes of the present work, as shown in Table 1. Case A is based on a high-Ra/high-Pr configuration, with the octadecane assumed as a reference PCM. The octadecane is a reasonable surrogate of paraffin-like PCMs and thus Case A can be viewed as close to a realistic paraffin-based LHTES application. Cases B and C have been taken into account due to the availability of previous LB-based computational studies [21], which are functional to the LB-PF methodology assessment. Most of the EP-related part of the present study has been focused on Case A. The EP implementation within the ANSYS Fluent software has already been described in Section 2.1. Here some details are added on the computational setup for the Case A simulations. The unsteady SIMPLE algorithm has been selected for pressure-velocity coupling, with implicit first-order time integration and up to 50 time step sub-iterations. Second-order discretization schemes have been applied for spatial gradients and convective terms. In terms of nondimensional time units, a fixed time step length of ∆θ = 4 · 10 −7 has been set with a total simulation time of θ = 0.01. Since octadecane melts between 28°C and 30°C, to analyze the effects of the phase transition temperature range three values of ∆T ls have been considered: 1°C, 1.5°C and 2°C. Regarding A mush , previous reviews on PCM modeling with the EP methos have highlighted a practical range of 10 5 -10 8 for paraffin-like materials [12]: as such, here we start from an intermediate value, 5 · 10 6 , which is then raised up to 10 7 for each ∆T sl . The left and right walls are treated as fixed-temperature walls, while the top and bottom walls are assumed adiabatic. Note that the fixed temperature condition at both vertical walls is coherent with the numerical studies reported in [14]. The cavity is squared in shape (H = L = 0.1 m) and three computational grids have been considered to assess the effects of numerical accuracy: Grid #1, with 93 × 50 elements, Grid #2, with 143 × 100 elements and Grid #3, with 193 × 150 elements.
For the setup of the LB-PF simulations of Cases B and C, the reader is redirected towards reference [16].

Results
Figure 2 summarizes Nusselt number EP results. In Figure 2a, effects of the numerical grid resolution are highlighted: the average Nu variation across the simulation time span is of about 10% passing from Grid #1 to Grid #2 and about 5% passing from Grid #2 to Grid #3. Further refinements did not produce significant variations, so Grid #3 has been selected as the reference for all the subsequent simulations. Figures 2b-d show the cross-effects of A mush and ∆T ls . Apparently, graphs for ∆T ls = 1°C return the best fit with the theoretical Nu scaling, while for ∆T ls = 2°C the influence of A mush becomes small but the Nu trend is less consistent, especially during the mixed regime (Nu recovery is too slow).
To extend the analysis, contours of the liquid fraction γ have been analyzed at θ = 0.01 for different A mush /∆T ls configurations. For A mush = 5 · 10 6 and ∆T ls = 1°C there is a very good agreement with the theoretical Nu trend, but at θ = 0.01 the melting front has already been reached the right wall (Figure 3a), which is not consistent with previous numerical studies of the same Case [14]. Raising A mush improves the melting front position, but introduces also numerical wiggles (see Figure 3b) due to instabilities in the energy equation.
Based on this, the setup with A mush = 5 · 10 6 and ∆T ls = 1.5°C has been selected as the optimal one for the correct representation of the Nu scaling and melting front evolution. A snapshot of the liquid fraction and temperature distributions at θ = 0.01 is shown in Figure  4 for the optimal setup: the presence of the left wall and melting interface thermal boundary layers is apparent from the temperature field. Figure 5 collects Nu trends predicted by LB-PF for the Cases B and C, including also the reference correlation and simulations from [21]. Interestingly, both numerical data sets return a slower Nu dynamics compared to what theoretically expected. However, the Nu recovery computed by LB-PF is faster and, additionally, the predicted asymptotic convective Nu value is much more accurate.
Starting from these considerations, a further comparison has been made between LB-PF and EP, based on Case C. The EP simulation has been performed implementing the optimal setup previously found for Case A and adapting the thermophysical parameters to match the Pr, S te and Ra parameters. Surprisingly, Figure 6 reveals that the Nu dynamics predicted by LB-PF is still the faster one, although both approaches are in excellent agreement with the theoretical asymptotic Nu value. A plausible explanation for the slower initial heat transfer dynamics in both numerical approaches might be found in the Pr number: one of the basic assumptions behind Equation (10) is Pr ≥ 1, meaning that Pr = 1 is at the limit of its validity range. It should be also noted that Pr ∼ 1 is not a value of interest for actual LHTES applications, since it is typical of gaseous susbtances.

Conclusions
Outcomes from the study presented in the previous Sections allow to draw the following concluding remarks: • the combination of standard EP methodology and state-of-the art commercial CFD packages essentially confirms its reliablity in the simulation of high Ra/high Pr PCM melting problems; 8 E3S Web of Conferences 312, 01002 (2021) https://doi.org/10.1051/e3sconf/202131201002 76° Italian National Congress ATI • using relatively large ∆T ls values (1.5°C or more) lessens the impact of A mush , but might slighlty impair the results' consistency if the reference PCM is a pure substance (i. e. with very narrow phase transition range); • using large A mush seems to improve some aspects of the numerical predictions, but is also very likely to introduce melting front oscillations due to instabilities in the energy equation: as such, more conservative values (5 · 10 6 or even less) are preferable; • in moderate Rayleigh number cases, the proposed LB-PF algorithm has disclosed a clear potential for added accuracy and flexibility compared to standard EP.
On the LB-PF side, since the explicit LB fluid dynamic algorithm requires very small time steps (of the order of ∼ 1 · 10 −10 in nondimensional θ units), future developments will be mainly focused on: i) the efficient parallelization of the simulation code; ii) the expansion of the stability constraints on the fluid and thermal relaxation times. Both development lines will allow to extend the LB-PF applicability towards high Ra/high Pr regimes of practical interest. Figure 6: LBPF vs. EP comparison at Ra = 6.8 · 10 6 .