On Numerical Modeling of Thermal Performance Enhancementof a Heat Thermal Energy Storage System Using a Phase Change Material and a Porous Foam

: In this investigation, a comprehensive numerical analysis of the ﬂow involved in an open-ended straight channel fully ﬁlled with a porous metal foam saturated and a phase change material (parafﬁn) has been performed using a single relaxation time lattice Boltzmann method (SRT-LBM) at the representative elementary volume (REV) scale. The enthalpy-based approach with three density functions has been employed to cope with the governing equations under the local thermal non-equilibrium (LTNE) condition. The in-house code has been validated through a comparison with a previous case in literature. The pore per inch density (10 ≤ PPI ≤ 60) and porosity (0.7 ≤ ε ≤ 0.9) effects of the metal structure were analyzed during melting/solidifying phenomena at two Reynolds numbers (Re = 200 and 400). The relevant ﬁndings are discussed for the LTNE intensity and the entropy generation rate (Ns). Through the simulations, the LTNE hypothesis turned out to be secure and valid. In addition, it is maximum for small PPI value (=10) whatever the parameters deemed. On the other hand, high porosity (=0.9) is advised to reduce the system’s irreversibility. However, at a moderate Re (=200), a small PPI (=10) would be appropriate to mitigate the system irreversibility during the charging case, while a large value (PPI = 60) might be advised for the discharging case. In this context, it can be stated that during the melting period, low porosity (=0.7) with low PPI (=10) improves thermal performance, reduces the system irreversibility and speeds up the melting rate, while for high porosity (=0.9), a moderate PPI (=30) should be used during the melting process to achieve an optimal system.


Introduction
Latent Heat Thermal Energy Storage (LHTES) mode continues to gain great attention due to its large amount of energy stored in a small volume at almost constant temperature.
The LHTES systems using phase change materials (PCMs) have been applied in several engineering applications such as concentrating solar plants, thermal isolation, groundwater infiltration, solar energy, building sectors, nuclear reactors, and so on thanks to the phase transition process [1,2]. However, the only issue of these systems is the low thermal conductivity of PCM. Thereby, the integration of porous metallic structures with PCMs is one of the most suitable solutions and widely used owing to their large thermal conductivity, high specific surface area and their low cost [3]. In this area, representative studies can be found such as that of Mabrouk et al. [2] who performed a numerical study on the heat transfer improvement in a rectangular channel including a metallic foam/paraffin. They examined the porosity and pore density (PPI) effects on the thermal performance of the system under consideration. They demonstrated that high porosity with large PPI speeds up the melting phenomenon. In addition, an increase in PPI restricts the forced convection and improves the heat conduction. In a further numerical investigation of Mabrouk et al. [4], different parameters have been considered such as the porosity, Reynolds and Eckert numbers to deal with the melting/solidifying processes. They pointed out that the quantity and quality of stored energy is maximum for a low porosity value, while at moderate Re (=400) the energetic and exergetic efficiencies are optimum for a high porosity. Sardari et al. [5] have numerically conducted some morphological parameters of metal foams such as the porosity and pore size on the PCM melting progression and demonstrated that lowering porosity enhances the system performance and accelerates the melting rate compared to a pure PCM. They stated that the metal pore size has sparsely effect.Tao et al. [3] investigated the influence of the same parameters, viz., PPI and porosity on the melting process in a square cavity filled with metal foams and paraffins subjected to the natural convection. They indicated that, for a greater porosity (=0.98) with a decrease in PPI, natural convection strongly dominates the heat transfer and that the melt front evolvement advances rapidly. Zadeh et al. [6] surveyed the phase change transition in a circular thermal energy storage unit partially filled porous copper foam and nano-additives. They found that an increase in the porosity improves latent heat capacity due to the PCM volume rises, slows the melting rate and intensifies heat transfer by convection. Zhang et al. [7] numerically analyzed the heat transfer enhancement of PCM melting/solidification through a combination of heat pipe-fins-copper foam (HP-Fin-CF). They reported that when the porosity of the HP-CF combination decreases, the heat transfer rate of phase change processes of PCM rises obviously at fixed PCM volume. Also, they found the same as the number of fins increases. Nevertheless, this combination strongly restricts natural convection and strengthens conduction. Esapour et al. [8] presented a model of multi-tube heat exchanger based on a metal foam/PCM composite. They found that the melting time is reduced by 14% and 55% for porosities 0.9 and 0.7, respectively indicating thereby that a porous metal foam insertion is more efficient during PCM solidification period than that of melting. Yang et al. [9] developed an experimental and numerical study on the PCM (water) solidification features impregnated into metal foam samples with different pore parameters such as porosity (0.93 and 0.97) and pore density (8 and 30 PPI) for cold storage systems. They showed that the time interval of solidification process was reduced by 87.5% and 76.7% with porosities of 0.93 and 0.97, respectively compared with pure PCM. Furthermore, they demonstrated that porosity rather than pore density represents the key parameter that influenced and dominated heat transfer enhancement of the phenomenon.
As a relevant numerical approach for coping with complex physics and heat transfer problems, the lattice Boltzmann method (LBM) showed a sound ability over the past 30 years compared to traditional computational fluid dynamics (CFD) methods (volumes and finite elements, etc.). It is worth noting that, when dealing with transport problems in porous medium, the LBM approach is mostly split into two approaches, viz., the LBM at the pore scale and the LBM at the representative elementary volume (REV) scale [2,4]. Various studies [2][3][4], to name a few, have demonstrated that LBMs have become prominent approaches in handling fluid-fluid and fluid-solid interactions within complex models such as porous structures.
Through the works briefly mentioned herein up, it can be stated that the novelty of this paper is the study of the effects of PPI and low porosities under forced convection with phase change transition for latent heat thermal energy storage system. To the best of our knowledge, these effects have been scarcely dealt. Besides, numerical and experimental investigations available in the literature that deal with the system irreversibility are sorely lacking.
This numerical study leans on an enthalpy-based TLBM approach to investigate the influence of the metal foam pore density and porosity on heat transfer enhancement during the phase transition phenomenon of PCMs. The main aim targeted here is to analyse the entropy generation rate (Ns) and the LTNE criteria under the effect of considered parameters and two Reynolds numbers. This paper is organized as follow. After Section 1, the physical and mathematical models are exhibited in Section 2. Afterwards, the numerical method (TLBM) is described and the in-house code is validated with available cases in literature in Section 3. In Section 4, numerical findings are illustrated and analysed. Finally, the main outcomes are highlighted in Section 5.

Physical Model
The physical model with two-dimensional system deemed herein is illustrated in Figure 1.
It is an open-ended straight channel with height H and length L fully filled with a porous structure (metal foam) and saturated with a phase change material: PCM (paraffin). The top and bottom walls are insulated (adiabatic) and non-slip. During the charging process, the fluid flow enters through the porous channel at a high temperature T h and a velocity U 0 , so that the paraffin begins to melt. Then the fluid leaves the channel through the east with lower constant temperature T c . Inversely, during the discharge process, the cold fluid enters the channel from the east at −U 0 . At this point, the liquid paraffin releases more heat energy and begins to solidify. So, the fluid heats up and set (force) a flow which leaves the channel at hot temperature T h . This numerical study leans on an enthalpy-based TLBM approach to investigate the influence of the metal foam pore density and porosity on heat transfer enhancement during the phase transition phenomenon of PCMs. The main aim targeted here is to analyse the entropy generation rate (Ns) and the LTNE criteria under the effect of considered parameters and two Reynolds numbers. This paper is organized as follow. After Section 1, the physical and mathematical models are exhibited in Section 2. Afterwards, the numerical method (TLBM) is described and the in-house code is validated with available cases in literature in Section 3. In Section 4, numerical findings are illustrated and analysed. Finally, the main outcomes are highlighted in Section 5.

Physical Model
The physical model with two-dimensional system deemed herein is illustrated in

Conjectures
To facilitate the mathematical model resolution of the problem, the following hypotheses are taken into consideration: The fluid flow is laminar, Newtonian and incompressible where forced convection mode prevails in the channel. The fluid and solid phases are in local thermal non equilibrium (LTNE) condition. The thermo-physical characteristics are assumed to be constant, homogeneous and isotropic. For more details, the reader may consult Nield and Bejan [10].

Mathematical Formulation
Based on the above conjectures, the two-dimensional governing equations for mass, momentum, and energy are as follows [2,4,11]:

Conjectures
To facilitate the mathematical model resolution of the problem, the following hypotheses are taken into consideration: The fluid flow is laminar, Newtonian and incompressible where forced convection mode prevails in the channel. The fluid and solid phases are in local thermal non equilibrium (LTNE) condition. The thermo-physical characteristics are assumed to be constant, homogeneous and isotropic. For more details, the reader may consult Nield and Bejan [10].

Mathematical Formulation
Based on the above conjectures, the two-dimensional governing equations for mass, momentum, and energy are as follows [2,4,11]: Computation 2022, 10, 3 Subscripts f and s indicate the fluid and solid phases, respectively. The variables → u , P, T f , T s , ε, ρ, ν f , C p , λ eff , Γ, L a , a sf , h sf , K and F ε represent the velocity vector, pressure, fluid and solid temperatures, metal foam porosity, density, paraffin kinematic viscosity, specific heat capacity, effective thermal conductivity, liquid fraction, latent heat, specific surface area, interfacial heat transfer coefficient, permeability and the Forchheimer coefficient, respectively.
The liquid fraction Γ and the empirical quantities a sf , h sf and λ eff can be calculated as follows [11][12][13]: where d f , d p , and ω denote the ligament diameter, pore size, and the pore density, respectively. Under the LTNE assumption, the effective thermal conductivity can be correlated using the following relationships [11,14,15]: with e = 0.16 and σ = As for the medium permeability and the Forchheimer coefficient (Equation (2)), they have been computed as [11][12][13]:

Boundary and Initial Conditions (B & ICs)
The boundary and initial conditions of the physical model deemed are the following:  Note that the fluid flow evolves from right to left, as depicted in Figure 1, thereby initiating the discharging process.

Entropy Generation Rate
As the entropy generation reflected system irreversibility causing by heat transfer irreversibility (HTI) and fluid friction irreversibility (FFI), the local and average entropy generation rate under the LTNE assumption can be formulated as follows [16,17]: D and τ being the strain and stress tensors, respectively. The average entropy generation Ns ave per unit of volume can be calculated as [17]:

LTNE Intensity
The LTNE condition is computed via the following criterion [17]: N and Θ f;s are the total number of nodes and the dimensionless temperature, respectively.
It should be noted that if LTNE > 5%, it is the LTNE condition that prevails, and if LTNE < 5%, the local thermal equilibrium (LTE) condition takes place.

LB Equations
The thermal lattice Boltzmann method (TLBM) controls the propagation of particles on a lattice x at time t with a discrete distribution velocities. In this study, three distribution functions are used to handle the dynamic, f i (x, t), and thermal, g i,f,s (x, t), fields using the single relaxation time SRT approach at the representative elementary volume (REV) scale. Note that the use of the SRT model could cause numerical diffusion when the collision time is not the unity [11]. To deal with this problem, the multiple relaxation time (MRT) approach should be considered [18]. Despite this weakness, the SRT approach remains the most popular and has shown its ability to handle phase change in porous media due to its simplicity, numerical precision and computational efficiency while mitigating the numerical diffusion [11,19] across phase interface. Thereby, the LBEs are expressed as [11,19]: where τ v (= 3ν + 0.5) and τ T,f;s are the dimensionless collision time for velocity and for temperatures [17,19]: τ T,f = 3α e,f /(δtc 2 ) + 0.5 with α e,f = λ e,f / ε(ρC p ) f (19) τ T,s = 3α e,s /(δtc 2 ) + 0.5 with α e,s = λ e,s / (1 − ε)(ρC p ) s (20) δt, c (= δx/δt = 1; δx = δt) and α e,f,s being the lattice time step, the streaming speed and the effective diffusivity, respectively. Under the D2Q9 model (Figure 2), the local equilibrium distribution functions f eq i (x, t) and g eq i,f;s can be expressed as follows [4,11,17]: g eq f,i = w i T f 1 + e i u/(εc 2 s ) and g eq s,i = w i T s (22) → e i and w i being the particle streaming velocity and the equilibrium weighting coefficients, which are written, respectively, as follows:   The source terms → F e i , Sr i,f;s and q i are computed as follows [4,11,17]: Finally, the macroscopic quantities ρ, → u , T f , T s are computed as: The solution procedure is shown through the flow chart (see Figure 3).
Finally, the macroscopic quantities f s , u, T , T   are computed as: The solution procedure is shown through the flow chart (see Figure 3).

Validation
Figure 4 depicts our in-house code through the study by Krishnan et al. [20]. As it can observed, the evolution dimensionless fluid (Θ f ) and solid (Θ s ) temperatures for various dimensionless times corroborate the results of the Ref. [20].

Validation
Figure 4 depicts our in-house code through the study by Krishnan et al. [20]. As it can observed, the evolution dimensionless fluid ( f  ) and solid ( s  ) temperatures for various dimensionless times corroborate the results of the Ref. [20].

Grid Test
To ensure the stability of the in-house code, grid independence test for four mesh

Grid Test
To ensure the stability of the in-house code, grid independence test for four mesh size, viz., 25 × 50, 50 × 100, 100 × 200 and 150 × 300 has been performed for Reynolds number of 400. The dimensionless fluid temperature (Θ f ) is portrayed in Figure 5. As can be seen, the profiles are close to each other. Thereby, the difference between the two last meshes is approximately 0.2%. However, it is about 0.8% between the last grid and the two first. Therefore, the 100 × 200 grid was selected for all subsequent simulations.

Grid Test
To ensure the stability of the in-house code, grid independence test for four mesh size, viz., 25 50  , 50 100  ,100 200  and150 300  has been performed for Reynolds number of 400. The dimensionless fluid temperature ( f  ) is portrayed in Figure 5. As can be seen, the profiles are close to each other. Thereby, the difference between the two last meshes is approximately 0.2%. However, it is about 0.8% between the last grid and the two first. Therefore, the 100 200  grid was selected for all subsequent simulations. Pr 50, Kr 10 , Rc 1, Da 10 ,Ec 0, 0.7, 10 PPI, Re 400 This section is devoted to analyzing the effects of pore density  ( PPI 10, 30  and 60) and Re number (200 and 400) while keeping the other parameters set at Pr 50, Ste 1, Ec 0, Rc 1     . This section is devoted to analyzing the effects of pore density ω (PPI = 10, 30 and 60) and Re number (200 and 400) while keeping the other parameters set at Pr = 50, Ste = 1, Ec = 0, Rc = 1. Figure 6 depicts the effects of the pore density (=10, 30 and 60) and porosity on the LTNE parameter for Re = 200 during the charging and discharging processes. As shown, it can be stated the LTNE condition is secured (LTNE > 0.05) for Re=200. For the charging case, it decreases for ε = 0.7. However, it increases for ε = 0.8 and 0.9 and PPI ≤ 30 then, it decreases ( Figure 6a). For the discharging case (Figure 6b), it decreases whatever ε. It is observable that for Re = 200, the LTNE is maximum for small porosity (=0.7) and PPI (=10) values during both charging/discharging periods. Figure 7 illustrates the effects of the pore density (=10, 30 and 60) and porosity on the LTNE parameter for Re =400 during the charging and discharging processes. During the charging period (Figure 7a), the LTNE increases proportionally with the PPI for ε = 0.8 and 0.9 and inversely, it decreases for low porosity value (=0.7) with increasing PPI. However, during the discharging period, it decreases with the increasing PPI regardless of the porosity. For the case of Re = 400, The LTNE is maximum for lower PPI (=10) during two processes. In addition, it turns out that the increase in Re secures the LTNE validity.

PPI's Effect on the LTNE Condition
To sum up, a decrease in PPI for a low porosity value (=0.7) gives a maximum LTNE value during the melting case. On the other hand, during the discharging (solidifying) period, decreasing the PPI for a low porosity (ε ≤ 0.8 ) increases the LTNE criteria whatever Re. Figure 6 depicts the effects of the pore density (=10, 30 and 60) and porosity o LTNE parameter for Re =200 during the charging and discharging processes. As sh it can be stated the LTNE condition is secured (  (Figure 6a). For the discharging case (Figure 6b), it decr whatever  .It is observable that for Re=200, the LTNE is maximum for small por (=0.7) and PPI (=10) values during both charging/discharging periods.
(a) (b) Figure 6. Pore density effect on LTNE criteria for Re=200 during: (a) charging case; (b) discha case. Figure 7 illustrates the effects of the pore density (=10, 30 and 60) and porosity o LTNE parameter for Re =400 during the charging and discharging processes. Durin charging period (Figure 7a), the LTNE increases proportionally with the PPI for   and 0.9 and inversely, it decreases for low porosity value (=0.7) with increasing However, during the discharging period, it decreases with the increasing PPI regar of the porosity. For the case of Re=400, The LTNE is maximum for lower PPI (=10) du two processes. In addition, it turns out that the increase in Re secures the LTNE valid To sum up, a decrease in PPI for a low porosity value (=0.7) gives a maximum L value during the melting case. On the other hand, during the discharging (solidif period, decreasing the PPI for a low porosity ( 0.8   ) increases the LTNE cr whatever Re. (a) (b) Figure 6. Pore density effect on LTNE criteria for Re=200 during: (a) charging case; (b) disch case. Figure 7 illustrates the effects of the pore density (=10, 30 and 60) and porosity o LTNE parameter for Re =400 during the charging and discharging processes. Durin charging period (Figure 7a), the LTNE increases proportionally with the PPI for  and 0.9 and inversely, it decreases for low porosity value (=0.7) with increasing However, during the discharging period, it decreases with the increasing PPI rega of the porosity. For the case of Re=400, The LTNE is maximum for lower PPI (=10) d two processes. In addition, it turns out that the increase in Re secures the LTNE val To sum up, a decrease in PPI for a low porosity value (=0.7) gives a maximum value during the melting case. On the other hand, during the discharging (solidi period, decreasing the PPI for a low porosity ( 0.8   ) increases the LTNE c whatever Re.

PPI's Effect on the Dimensionless Entropy Generation Rate
The entropy generation rate is an indicator used to assess the system stability. Figures 8-10 outline the pore density effect on the dimensionless entropy generation rate during changing and discharging processes for Re = 200 and 400 at ε = 0.7, 0.8 and 0.9.
For both cases charging/discharging and the three porosities, it can be seen that the entropy generation increases proportionally with Re number. Hence, lower Re should be used to reduce the system irreversibility, which refers to the dominance of the heat transfer irreversibility.
As seen in some cases in Figure 10 for example for 0.9   , the Ns's amplitude d ing the discharging period exceeds that of the charging for Re=400, indicating the portant role of the metal structure characteristics on the heat transfer irreversibility.
To sum up, to reduce the irreversibility of the system, lower PPI with lower should be used during the charging period, while a larger PPI could be adopted dur the discharging case. In addition, high porosity (=0.9) should be used to mitigate system irreversibility whatever the parameters deemed.

Conclusions
The current numerical investigation reported the porosity, pore density and number effects on heat transfer subjected to forced convection in an open-ended recta gular channel fully filled with porous metal foam saturated with paraffin. The compu tions have been performed using an enthalpy based SRT-TLBM approach at the RE scale. Validation of the built code with previous casein literature has exhibited go

Conclusions
The current numerical investigation reported the porosity, pore density and number effects on heat transfer subjected to forced convection in an open-ended rec gular channel fully filled with porous metal foam saturated with paraffin. The comp tions have been performed using an enthalpy based SRT-TLBM approach at the R During charging process, for ε = 0.7 and 0.8 (Figures 8 and 9), the entropy generation rate is minimal for low PPI value (=10) regardless of Re. It is minimal for PPI = 30 for higher porosity (=0.9) as exhibited in Figure 8. Indeed, it is maximal for larger PPI (=60) whatever Re and ε. Furthermore, as Re increases (=400) and PPI is stronger (=60), the dimensionless entropy generation rate reaches it is maximum indicating that the forced convection dominates the heat transfer irreversibility. As for discharging process, at the flow entrance, the dimensionless entropy increases with the increase of PPI showing that the forced convection dominates the channel. Then, from the channel's 1/5, it decreases as PPI increases due to the weakness of forced convection and subsequently the dominance of the heat conduction on the global heat transfer.
As seen in some cases in Figure 10 for example for ε = 0.9, the Ns's amplitude during the discharging period exceeds that of the charging for Re = 400, indicating the important role of the metal structure characteristics on the heat transfer irreversibility.
To sum up, to reduce the irreversibility of the system, lower PPI with lower Re should be used during the charging period, while a larger PPI could be adopted during the discharging case. In addition, high porosity (=0.9) should be used to mitigate the system irreversibility whatever the parameters deemed.

Conclusions
The current numerical investigation reported the porosity, pore density and Re number effects on heat transfer subjected to forced convection in an open-ended rectangular channel fully filled with porous metal foam saturated with paraffin. The computations have been performed using an enthalpy based SRT-TLBM approach at the REV scale. Validation of the built code with previous casein literature has exhibited good agreement. To cope with the present study, three distribution functions for dynamic and thermal fields have been applied under the LTNE condition. The relevance of the pore density effect is highlighted and analyzed during the melting and solidifying processes of the PCM.Based on the numerical outcomes obtained, the main conclusions are as:

•
Enthalpy-based REV-TLBM method has a robust ability for handling the phase change phenomena in a porous channel subjected to steady forced fluid flow.

•
During charging and discharging process, the LTNE hypothesis is valid for the case of Re = 200 and 400 whatever the PPI and ε. • For the melting and solidifying periods, small PPI (=10) and porosity (ε ≤ 0.8) give a maximum LTNE regardless of Re.

•
The system irreversibility can be reduced via a small PPI (=10) and Re (=200) during the charging case, while a large value (PPI = 60) can be used during the discharging period.

•
High porosity (=0.9) is recommended to mitigate the system irreversibility whatever the parameters deemed. • From the channel's 1/5 and during the solidification process, thermal conduction dominates the overall heat transfer as the PPI increases. • During the charging process, a low porosity value (=0.7) with a low PPI (=10) improves heat transfer, reduces the system irreversibility and speeds up the melting rate, while for a high porosity (=0.9), a moderate PPI (=30) should be considered an optimal value during the fusion loop. This was supported in Ref. [2]. • During the solidification process, an increase in PPI gives rise to a thermal conduction improvement. For more details, one could consult Ref. [3].
The use of renewable energies such as phase change in porous media, however, faces serious challenges, one of the main ones being its small volume at almost constant temperature and fitting consideration of pore geometry. Recent advances (especially in PCMs for TES) should help better overcome these barriers. To keep future research directions consistent with major issues, the following future research can be addressed:

•
Tracks such as the use of porous media with inexpensive MCPs can be followed to better deepen our knowledge of such systems.
• PCMs and latent heat energy storage systems should be also economically analysed to assess their feasibility on commercial scale, as the initial cost for setting up such systems may turn out to be high.
To sum up, this study was initiated to better understand the nature and scope of a real or artificial system able of presenting storage performances to meet the needs. Funding: This study was funded by a scholarship that was granted by the National School of Engineers of Monastir/Monastir University as part of the doctoral thesis of the first author.
Data Availability Statement: Data sharing is not applicable to this article.

Acknowledgments:
The authors would like to thank to the National School of Engineers of Monastir/Monastir University for supporting this work. All authors have read and approved the present manuscript version. In addition, the authors' order was approved by mutual agreement between the authors.

Conflicts of Interest:
The authors declare no conflict of interest.