Alissar Yehya
email: alissar_yehya@seas.harvard.edu

Hassane Naji
email: hassane.naji@univ-artois.fr

Laurent Zalewski
email: laurent.zalewski@univ-artois.fr

Experimental and numerical characterization of an impure phase change material using a thermal lattice Boltzmann method

Keywords: Phase change materials (PCMs), Lattice Boltzmann method (LBM), Convection, Supercooling, Impurities

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The global population growth has led to a high energy demand and increased CO2 emissions triggering a major environmental crisis. The energy use in buildings accounts for a large share of the total end use of energy, with an average consumption of about 50% of all raw materials owned in the world [START_REF] Hegger | Energy Manual, Sustainable architecture[END_REF]. This percentage would differ depending on the region, and is equal to 40% in Europe, 30% in USA, and reaches 70% in undeveloped countries [START_REF]Energy Efficiency Requirements in Building Codes, Energy Efficiency Policies for New Buildings[END_REF]. Energy consumption for thermal comfort in buildings rises with increasing in user demand for comfort conditions. This rise in, both, energy consumption and CO2 emissions thus promoted a new policy aimed at creating more sustainable buildings [START_REF] Castell | Experimental Study of PCM in Mediterranean Buildings, Effstock 2009 Conference[END_REF]. Scientists all over the world are in search of new and renewable energy sources. One of the options is to develop energy storage devices, which are as important as developing new sources of energy. Thermal energy storage does not only reduce the temporal gap between supply and demand, but it also improves the performance of energy systems and plays a significant role in conserving energy [START_REF] Sharma | Review on thermal energy storage with phase change materials and applications[END_REF]. Such storage can be achieved either by using sensible heat or thermo-chemical heat or latent heat. Among these different methods, the latent thermal energy storage (LTES) method, using phase change materials (PCMs), is the most preferred and the most widely used at present [START_REF] Soaresa | Review of passive PCM latent heat thermal energy storage systems towards buildings' energy efficiency[END_REF][START_REF] Baetens | Phase change materials for building applications: a state-of-the-art review[END_REF] mainly for their application for heating and cooling of buildings [START_REF] Sharma | Review on thermal energy storage with phase change materials and applications[END_REF][START_REF] Zalba | Review on thermal energy storage with phase change materials: Heat transfer analysis and applications[END_REF][START_REF] Farid | A review on phase change energy storage: materials and applications[END_REF][START_REF] Zhang | Application of latent heat thermal energy storage in buildings: State of the art and outlook[END_REF][START_REF] Dutil | A review on phase-change materials: Mathematical modeling and simulations[END_REF][START_REF] Voller | Fixed grid techniques for phase change problems: a review[END_REF]. PCMs can be used for temperature regulation and their high storage capacity can reduce energy consumption in buildings, and thereby contribute to the reduction of CO2 emissions.

Recently, interest in the use of Octadecane in latent heat systems and building materials has increased. Its choice is mainly due to its transparency properties and its melting temperature, which lies within the thermal comfort of human beings. The original idea is to fill Plexiglas enclosures with PCM having transparency properties, and thereby create a wall of a building façade or windows that benefits from solar gains. For PCMs in windows, paraffin-based organic materials are the most interesting since they are transparent in the liquid state and translucent in the solid state. For that, it is interesting to investigate more on Octadecane as a PCM. The main task throughout this research is to develop a numerical model that matches the experimental results and helps infer some thermo-physical properties that are difficult to predict via single experiments, example: the thermal contact resistance, conductivity in liquid phase, the degree of supercooling, heat transfer coefficient, etc. Following this, the model can be also used to enhance and optimize PCM performance in practical applications. Note that, at this stage, the heat transfer by radiation is not considered since our experimental apparatus induces heat through attached plates controlled by thermo-regulated baths (see Section 4). Actually, the study of such materials will involve the characterization part, as our present study, and an application part where the experimental setup will be different and will follow appropriate climatic and external conditions and adequate geometric scales. Both tests and studies are required for this material to reach a mature state for practical applications.

Thus, prior to the large-scale practical application of this technology, it is crucial to consider certain issues in the stage of research and development [START_REF] Sharma | Review on thermal energy storage with phase change materials and applications[END_REF]. The experimental investigation of some PCMs showed the presence of some phenomena that could affect the performance if not considered in the design, namely supercooling, convection and impurities. When supercooling occurs, the start of solidification is delayed, and the liquid solidifies at a temperature below its freezing temperature.

This can lead to a mismatch between the design and the real behavior of the PCM [START_REF] White | An enthalpy formulation of the Stefan problem[END_REF]. Most of the theoretical heat transfer modeling of PCM energy storage systems is based on the conjecture that the supercooling degree is negligible or even completely absent. However, recently published researches [START_REF] White | An enthalpy formulation of the Stefan problem[END_REF][START_REF] Bonacina | Numerical solution of phase-change problems[END_REF][START_REF] Zeng | Temperature-transforming model for binary solid-liquid phasechange problems, Part I: Mathematical modeling and numerical methodology[END_REF] indicate that even a relatively small degree of supercooling might result in a significant undesirable effect on the heat release rate as a function of time, reducing the thermal efficiency of the system. Therefore, it is important to provide reliable theoretical modeling of supercooling in PCM for a more adequate design.

Another important phenomenon, ignored in some numerical models, is the natural convection in the liquid phase. If PCMs with high Prandtl number are used, e.g. Octadecane, it is very likely to have convection even in small samples. This is because the relative thickness of the momentum boundary layer is much more than the thermal one. Hence, convection will be dominant with momentum diffusing more quickly compared to heat. In addition, soluble impurities, if they exist, may lead to depression in melting point thus altering the temperature-enthalpy curve, so that the traditional piecewise function can no longer be used. Also, the thermal contact between different materials in the system can have similar effects on the enthalpy curve. Therefore, using a simplified model may fail to successfully predict the phase change process, if the aforementioned phenomena are present. Thereby, the main aim targeted here is to provide an improved model that accounts for supercooling, convection and the modified enthalpy-temperature curve. To carry out this work, we consider the thermal lattice Boltzmann method (TLBM) [START_REF] Talati | Lattice Boltzmann application on the PCM solidification within a rectangular finned container[END_REF][START_REF] Taghilou | Numerical investigation on the natural convection effects in the melting process of PCM in a finned container using lattice Boltzmann[END_REF][START_REF] Darzi | Melting and solidification of PCM enhanced by radial conductive fins and nanoparticles in cylindrical annulus[END_REF][START_REF] Su | A new mesoscopic scale timestep adjustable non-dimensional lattice Boltzmann method for melting and solidification heat[END_REF][START_REF] Jourabian | Melting of nanoparticles-enhanced phase change material (NEPCM) in vertical semicircle enclosure: numerical study[END_REF][START_REF] Jourabian | On the expedited melting of phase change material (PCM) through dispersion of nanoparticles in the thermal storage unit[END_REF] based on a partial bounce-back (PBB) approach [START_REF] Yehya | Thermal lattice Boltzmann simulation of entropy generation within a square enclosure for sensible and latent heat transfers[END_REF][START_REF] Yehya | Assessment of a Lattice Boltzmann Model to simulate fluid flows with complex geometries[END_REF][23] supplemented with an enthalpy-based model. This enhanced model provides a more accurate simulation of phase-change materials, can better predict the solidification and melting processes, and thus match the experimental results.

The paper is organized as follows: Section 2 provides a literature review on challenges in modeling phase-change materials. Section 3 is devoted to the lattice Boltzmann method adopted herein to achieve the numerical solution while supplementing the appropriate enhancements. The experimental setup with equipment and the measurement method is outlined in Section 4. Section 5 presents the experimental determination of the thermo-physical properties of Octadecane: conductivities, specific capacities, latent heat of fusion, and the analysis of the heat flux curves.

Thorough discussions of the confrontation between numerical and experimental results are presented in Section 6. Finally, Section 7 concludes the paper by providing key findings based on the dual experimental and numerical approach.

Challenges of modeling of phase-change materials

PCMs have various possible applications ranging from temperature stabilization to energy storage. Though many pure chemical elements can be used as PCMs, the majority often consists of a mixture, or contains at least soluble impurities. The main reason for creating a PCM as a mixture of various substances is to achieve a desirable melting temperature for a particular application.

However, these mixed or impure PCMs require accurate and reliable methods for determining their thermo-physical properties, which represents a significant issue that considerably affects the accuracy and credibility of the corresponding studies. Numerical modeling can help in determining some of these properties. However, other serious complexities interfere in the modeling of phase change materials, which makes it a real challenge.

Many macroscopic mathematical modeling schemes for solidification/melting problems can be found in the literature [23][START_REF] Voller | Fast implicit finite-difference method for the analysis of phase change problems[END_REF][START_REF] Swaminathan | Towards a general numerical scheme for solidification systems[END_REF]. Early efforts initiated with the moving/deforming grid approach [START_REF] Voller | Fast implicit finite-difference method for the analysis of phase change problems[END_REF][START_REF] Swaminathan | Towards a general numerical scheme for solidification systems[END_REF] in which independent conservation equations for each phase are initially formulated and are subsequently coupled with appropriate boundary conditions at the interfaces. However, such multiple region solutions require the existence of discrete interfaces between the respective phases.

In fact, a major difficulty in their implementation [START_REF] Chatterjee | Lattice Boltzmann Modeling for Melting/Solidification Processes[END_REF] is associated with tracking the phase interfaces (which are generally unknown functions of space and time). Additionally, a serious limitation exists for modeling phase change behavior of multi-component systems, since they do not exhibit a sharp interface between solid and liquid phases. Moreover, solidification occurs over extended temperature ranges and solid formation often occurs as a permeable crystalline-like matrix, which coexists with the liquid phase. Otherwise, there are the fixed-grid models, using finite element (FE) or finite volume (FV) methods [START_REF] Voller | Fast implicit finite-difference method for the analysis of phase change problems[END_REF][START_REF] Al-Saadi | Modeling phase change materials embedded in building enclosure: A review[END_REF]. A separate equation for the liquid fraction evolution is solved, which implicitly specifies and updates the interfacial locations with respect to space and time. The fixed grid method is relatively simple, versatile, practical, adaptable and easily programmable [START_REF] Voller | Fixed grid techniques for phase change problems: a review[END_REF]. It can easily handle melting or solidifying materials over a range of temperatures. The latent heat evolution is accounted for in the governing equation by using, either enthalpy method [START_REF] White | An enthalpy formulation of the Stefan problem[END_REF]23], heat capacity method [START_REF] Bonacina | Numerical solution of phase-change problems[END_REF][START_REF] Morgan | An improved algorithm for heat conduction problems with phase change[END_REF], temperature transforming model [START_REF] Zeng | Temperature-transforming model for binary solid-liquid phasechange problems, Part I: Mathematical modeling and numerical methodology[END_REF],

heat source method [START_REF] Voller | Fast implicit finite-difference method for the analysis of phase change problems[END_REF][START_REF] Swaminathan | Towards a general numerical scheme for solidification systems[END_REF], or other methods.

On the other hand, the multi-scale mesoscopic lattice Boltzmann method (LBM) has emerged to offer huge potential for solving complex thermo-fluidic problems involving morphological development of complicated phase boundaries [START_REF] Yehya | Thermal lattice Boltzmann simulation of entropy generation within a square enclosure for sensible and latent heat transfers[END_REF][START_REF] Yehya | Assessment of a Lattice Boltzmann Model to simulate fluid flows with complex geometries[END_REF][23][START_REF] Huber | Lattice Boltzmann model for melting with natural convection[END_REF][START_REF] Jourabian | Convection-dominated melting of phase change material in partially heated cavity: lattice Boltzmann study[END_REF][START_REF] Mezrhab | Double MRT thermal lattice Boltzmann method for simulating convective flows[END_REF][START_REF] Chen | Lattice Boltzmann Method for Fluid Flows[END_REF][START_REF] Xu | Lattice Boltzmann modeling of transport phenomena in fuel cells and flow batteries[END_REF] and recovers the Navier-Stokes and energy equations. In addition, it has been witnessed that the LBM stands out in major fields of classical fluid dynamics: multiple-scale flows and heat transfers with or without phase change process. Its advantage, compared to a classical based-continuum formulation, is that it operates at a mesoscopic level, which incorporates micro and meso-scale physics of phase transitions, and bypasses the explicit calculation of the pressure equation, leading to time-efficient computational simulations. Further, LBMs are inherently transient and parallelizable, which renders their suitability to address phase change processes over large-scale computational domains [START_REF] Xu | Accelerated lattice Boltzmann simulation using GPU and OpenACC with data management[END_REF].

Moreover, by adopting partial bounce-back approach with enthalpy formulation [START_REF] Talati | Lattice Boltzmann application on the PCM solidification within a rectangular finned container[END_REF], LBM treats phase-change as flow in porous medium. Extensive work has been done to simulate fluid and heat flow in porous media using LBM [START_REF] Gao | Lattice Boltzmann simulation of natural convection dominated melting in a rectangular cavity filled with porous media[END_REF][START_REF] Gao | A thermal lattice Boltzmann model for natural convection in porous media under local thermal non-equilibrium conditions[END_REF]. The liquid fraction calculated through an enthalpy formulation designates the porosity in the mushy zone and helps tracking fusion front and simulating convection. The local nature of the method allows adding local complexities where necessary.

Numerical Model

Enthalpy-temperature relationship

The most common enthalpy/temperature relationship is often translated as a piecewise function in the case of pure materials. However, even a small amount of impurities can lead to a significant change in the enthalpy curve shape, with a depression of the melting temperature followed up by a melting range. Hence, melting point analysis can provide information on the sample purity [START_REF] Schoff | Purity determinations by thermal methods[END_REF]. In other words, a substance including soluble impurities will be prone to "melting point depression". This is related to the intermolecular forces within the material. As the solid impurity increases, its structure will be more disrupted, and hence will result in a greater variation in intermolecular forces throughout different areas of the solid. The effect will be that the melting temperature is lower than that of a pure material, and the solid melts over a wider range of temperatures. Note that this is not the case for insoluble impurities. The latter will have no effect on the compound's melting point.

When soluble impurities exist, the model can be dealt as a binary mixture, with one component having very small portion compared to the other. Thereby, we can calculate the depression in the melting temperature (ΔTam) as follows [START_REF] Schoff | Purity determinations by thermal methods[END_REF]:

∆ = -= 2 / am a m a i f T T T RT X L ( 1
)
where Ta is the melting point of a 100% pure material, Tm is the lowered melting point of the impure material, Lf (J/mol) is the molar heat of fusion, R (= 8.314 J/mol/K) is the ideal gases universal constant, and Xi is impurities' fraction in the considered material. In this work, this can be used as a pre-estimate of the real melting temperature of the studied sample with respect to the expected melting temperature in literature.

It seems obvious that impurities' percentage in a material alters the enthalpy curve of a substance and consequently affects the resulting heat flux. The peak of the heat flux curves becomes wider and shorter as the impurity increases while maintaining the same underlying area (under it), which is related to the melting enthalpy.

• Enthalpy-temperature relationship for pure PCMs

The variation of enthalpy vs. temperature of a pure material is:

<   =  >   , , s pcm a pcm l pcm a c for T T dh c for T T dT (2)
It should be noted that / pcm dh dT tends to infinity when T approaches Ta.

• Enthalpy-temperature relationship for non-pure PCMs, i.e. with soluble impurities For an impure material, the variation of the enthalpy vs. temperature can be written as [START_REF] Franquet | Experimental and theoretical analysis of a cement mortar containing microencapsulated PCM[END_REF][START_REF] Gibout | Comparison of different modelings of pure substances during melting in a DSC experiment[END_REF][START_REF] Dumas | Model for the DSC thermograms of the melting of ideal binary solutions[END_REF]:

- -  - - - ≤  - - =   >  , , , 2 , ( ) ( )
m a m a s pcm s pcm l pcm f m pcm a a l pcm m T T T T c c c L for T T dh T T T T dT c for T T (3)
The resulting expression of the enthalpy becomes:

     - - - + - - + + ≤      - - =       - + + >  , , , 0 , 0 ( ) ( )( )ln (
)

a m a s pcm m l pcm s pcm a m f m m a a pcm l pcm m f m

T T T T c T T c c T T L H for T T T T T T h c T T L H for T T (4)

where Ta and Tm are the melting temperatures of pure and impure PCMs, respectively. However, Lf is the latent heat of fusion, H0 is the enthalpy at TE, which denotes the onset of melting.

Thermal lattice Boltzmann model (TLBM)

The Lattice Boltzmann method (LBM) is a mesoscale discrete model that has become an increasingly popular tool for simulating fluid flows with and without heat transfers [START_REF] Talati | Lattice Boltzmann application on the PCM solidification within a rectangular finned container[END_REF][START_REF] Taghilou | Numerical investigation on the natural convection effects in the melting process of PCM in a finned container using lattice Boltzmann[END_REF][START_REF] Darzi | Melting and solidification of PCM enhanced by radial conductive fins and nanoparticles in cylindrical annulus[END_REF][START_REF] Su | A new mesoscopic scale timestep adjustable non-dimensional lattice Boltzmann method for melting and solidification heat[END_REF][START_REF] Jourabian | Melting of nanoparticles-enhanced phase change material (NEPCM) in vertical semicircle enclosure: numerical study[END_REF][START_REF] Jourabian | On the expedited melting of phase change material (PCM) through dispersion of nanoparticles in the thermal storage unit[END_REF]. It consists of simulating the statistical behavior of a set of particles on a lattice with finite velocities.

This evolution is carried out in a cycle of "streaming" (advection) and "collision" steps. The essential interpretation of such an approach is that it is a special finite difference form of the continuous Boltzmann equation. Furthermore, it allows providing macroscopic fluid properties, such as density, velocity, pressure, etc., through weighted averages, or moments, of the particle distribution for all discrete lattice velocities. The single relaxation time (SRT) lattice Boltzmann model (also called the lattice Bhatnagar-Gross-Krook (LBGK) model) [START_REF] Chapman | The mathematical theory of non-uniform Gases[END_REF][START_REF] Guo | A coupled lattice BGK model for the Boussinesq equations[END_REF] for incompressible thermal flows leans on two distribution functions (DFs) with their corresponding evolution equations to solve the evolution of the two mesoscopic particle distribution functions, i f and i g , via the discretized lattice Boltzmann equation as following:

( ) ( ) ( ) ( ) ( ) ( )
1 eq i i i f i i f 0 f x e t 1 f x t f x t f x t g T T T τ τ β ρ - + + = - - - - ∆ , , , , / (5)
( ) ( ) ( ) ( ) ( )
1 eq g i i i h i i h g x e t 1 g x t g x t g x t S τ - + + = - - + , , , , (6)
where i e is the microscopic particle velocity in the i -direction, f τ and h τ are the dimensionless relaxation times, and eq i f and eq i g are local equilibrium distributions functions that can be computed from:

=   = =   =  i / for i / for i ,
 =  ± =  =  ± =   ± ± =  (9)
are, respectively, the weight coefficient and the velocity vector of the D 2 Q 9 model; is the macroscopic velocity, with, u and v representing velocities in the x -and y -directions, respectively. Note that the relaxation times f Fτ ρ as proposed by Shan and Chen [START_REF] Chen | Simulation of nonideal gases and gas-liquid phase transitions by the lattice Boltzmann equation[END_REF],

where b F ( ) ( ) = - - ∆ 0 g T T / T β
. By this treatment, there is no need to add a force term to the collision operator. On the other hand, in g-distribution function, the source term is treated as per the method proposed by Luo [START_REF] Luo | Lattice-Gas automata and lattice Boltzmann equations for two-Dimensional hydrodynamics[END_REF]. Hence, the resulting force in the LBM frame will be:

g h i h S S ω = - with h S ( ) - = ∂ ∂ 1 Ste / t
ε being the source (or sink) term that handles the phase-change. Note that, we performed tests on the temperature and liquid fraction to insure that they converge within an acceptable tolerance after one iteration for the parameters set considered in this study.

It is worth recalling that the description and validation of the model adopted herein can be found in [START_REF] Yehya | Thermal lattice Boltzmann simulation of entropy generation within a square enclosure for sensible and latent heat transfers[END_REF][START_REF] Yehya | Simulating flows in multi-layered and spatially-variable permeability media via a new gray lattice Boltzmann model[END_REF][START_REF] Yehya | Contribution to the experimental and numerical characterization of phase change materials: Consideration of convection, supercooling, and soluble impurities[END_REF]. The numerical model will calculate the temperature and velocity fields for both conduction and convection heat transfer modes. Afterwards, the liquid fraction will be calculated from the enthalpy while taking into consideration the impurities' presence. Hence, the enthalpy h is calculated from equation ( 4).

The liquid fraction ( ε ) is computed as:

( )

1 2 1 1 if 1 h if 0 ε ε ε ε =   =  =   (10)
where 1

ε is an indicator parameter that is equal to 1 as long as nucleation did not occur for the solidification phase. For nucleation to occur the temperature should fall below the nucleation temperature Tnuc, where the degree of supercooling is the estimated as (Tm -Tnuc). Note that this happens only in solidification, so that the solidification and melting curves are defined according to different indicators. The above numeric indicator controls what enthalpy curve is used for the calculation of liquid fraction according to Fig. 3, and then, ( )

< = =   -  = ≤ ≤ = +  -   >  s s pcm E 0 s 2 s l s f l s l 0 for h h c T H h h h for h h h h L h h 1 for h h ε , ( 11
)
It is useful to recall that, in the current model, the solid and liquid phases are defined according to the liquid fraction value. Therefore, a mushy zone state is assigned when the value of ε is between zero and one. In this case, the velocity field is partially bounced back and the macroscopic velocity is modified [START_REF] Talati | Lattice Boltzmann application on the PCM solidification within a rectangular finned container[END_REF][START_REF] Yehya | Simulating flows in multi-layered and spatially-variable permeability media via a new gray lattice Boltzmann model[END_REF][START_REF] Yehya | Contribution to the experimental and numerical characterization of phase change materials: Consideration of convection, supercooling, and soluble impurities[END_REF]. The procedure for implementing the partial bounce-back approach is described in [START_REF] Yehya | Simulating flows in multi-layered and spatially-variable permeability media via a new gray lattice Boltzmann model[END_REF][START_REF] Yehya | Contribution to the experimental and numerical characterization of phase change materials: Consideration of convection, supercooling, and soluble impurities[END_REF].

For the velocity field, the non-slip BCs are used for all the cavity walls. These are performed by the on-grid bounce-back (BB) boundary conditions:

+ = ( , 1) ( , ) i w j w f x t f x t (12) 13 w
x being the fluid node adjacent to the wall, and i and j represent two opposite lattice directions on the boundary site. Note that the BB conditions apply to the DF in non-parallel directions at a solid wall.

To specify a constant temperature at the left and right walls, we use the method proposed by Inamuro et al. [START_REF] Inamuro | Lattice Boltzmann method for a binary miscible fluid mixture and its application to a heat-transfer problem[END_REF]. Its principle is to substitute unknown DFs for a boundary point with local equilibrium values using an adjusted temperature to set the defined temperature at that point.

Specifically, the adjusted temperature on the left side can be expressed as: In the nearly incompressible formulation, the basic thermo-hydrodynamic properties, such as density ρ , momentum density, U ρ , and temperature, T , are defined as moments of the DFs, i f and i g , as follows,

= - - ∑ , 6 ' ( ) 1 3 h k w x T T g v ( 13
8 0 = = = ∑ i i i f ρ , 8 0 = = = ∑ i i i i U e f ρ , 8
0 i i i T g = = = ∑ (14)
To simulate the presence of a mushy zone (co-existence of liquid and solid), a step is added to mimic the bounce-back at a lattice node [23,[START_REF] Yehya | Simulating flows in multi-layered and spatially-variable permeability media via a new gray lattice Boltzmann model[END_REF][START_REF] Yehya | Contribution to the experimental and numerical characterization of phase change materials: Consideration of convection, supercooling, and soluble impurities[END_REF]. This condition redirects the incoming fluid packets prior to the collision step, as follows:

1 out col in i i i f ( x,t ) ( ) f f ( x,t ) ε ε = - +
where fi in and fi out denote the incoming and outgoing fluid packet densities at a node.

This ensures that for a completely liquid phase ( 1 ε = ), a normal collision is conducted.

However, when it is completely solid ( 0 ε = ), the flow is bounce-backed and blocked, whereas for values between 0 and 1, the flux is partially bounce-backed according to the solid fraction estimated at the node. Such a scheme has been verified in our previous work compared to analytical solutions and other numerical methods [21-23, 40, 41].

The velocity field is then modified in the mushy zone and expressed in terms of the in-coming fluid packet densities,

= ε * . U U ( 15
)
This modified velocity is in accordance with the requirements of the partial bounce-back approach [23,[START_REF] Yehya | Simulating flows in multi-layered and spatially-variable permeability media via a new gray lattice Boltzmann model[END_REF][START_REF] Yehya | Contribution to the experimental and numerical characterization of phase change materials: Consideration of convection, supercooling, and soluble impurities[END_REF]. As a result, the flow in the mushy zone will be dealt as flow in a porous medium and will therefore be governed by Darcy's law as proved by derived analytical solutions [23,[START_REF] Yehya | Simulating flows in multi-layered and spatially-variable permeability media via a new gray lattice Boltzmann model[END_REF][START_REF] Yehya | Contribution to the experimental and numerical characterization of phase change materials: Consideration of convection, supercooling, and soluble impurities[END_REF].

Experimental Setup

Recall that the main task targeted here is to develop a mutually compatible experimental and numerical method to characterize the studied PCM. The major enhancement in the current numerical model is to add natural convection and supercooling, adopt variable thermo-physical properties, and use the real enthalpy curves corresponding to the considered PCM (viz.

Octadecane). Let's point out that this PCM is an alkane hydrocarbon of chemical formula is CH3(CH2)16CH3 with Prandtl number is about 50. In general, for high Prandtl values, the momentum diffusivity dominates, and convection seems very efficient in transferring heat (energy).

To solve this physical problem, the TLBM, for phase change in multi-layered materials, is used herein. However, before confronting numerical simulations findings, we perform experiments on the Octadecane to characterize its main thermo-physical properties: conductivity (solid and liquid), heat capacity (solid and liquid), and latent heat of fusion. The numerical model, hereby, helps us to determine both the real enthalpy curve and the degree of supercooling in the material. Figure 2 shows the schema of the experimentally studied sample of Plexiglas containing Octadecane as PCM. Tables 1 and2 gather the properties of Octadecane and Plexiglas, used in the numerical model, as per the literature [START_REF] Messerly | Low-temperature thermal data for pentane, n-heptadecane, and n-octadecane. Revised thermodynamic functions for the nalkanes, C5-C18[END_REF][START_REF] Mehling | Heat and cold storage with PCM: An up to date introduction into basics and applications[END_REF] and verified experimentally with an acceptable tolerance.

The experimental setup and its associated equipment, used in this work, are schematically depicted in Figs. 2 and3(a, b). The studied phase-change material, viz. Octadecane (99% pure), is placed inside a parallelepiped Plexiglas enclosure and prone to a temperature gradient through exchange plates controlled by two thermo-regulated baths (Julabo Model 34 HE -1kW with a precision of about ±0.01 °C)) that define thermal conditions (see specifications in [START_REF] Joulin | Experimental investigation of thermal characteristics of a mortar with or without a micro-encapsulated phase change material[END_REF]). The exchange aluminum plates (500 mm × 500 mm × 19 mm) allow imposing temperature variations, with respect to selected time scales, on the two dominant faces (left and right faces). The upper and bottom faces, of the studied sample, are insulated by a 14.5 cm-thick polystyrene (a thermally insulating material). Such insulation also serves to minimize the heat transfer from the sample lateral faces into the surrounding medium. This is the guarded hot plate approach principle, often used when characterizing the thermal performance of PCM samples. The flux-meters used herein are "tangential gradients flux meters". The inserted sensors (Captecv, France) on the sample's both sides were pre-calibrated (using the comparative method) with a precision of about 2%, using the device described by [START_REF] Joulin | Experimental and numerical investigation of a phase change material: thermal-energy storage and release[END_REF]. Their surface, thickness and sensitivity are 210×140 mm 2 , 0.2 mm and 124 μV.W -1 .m -2 , respectively. The various sensors are connected to a multichannel multimeter (Keithley 2700) adapted to low level signal measurements. Experimental data are scanned at regular and adjustable time steps of 10 s. The adopted experimental setup permits the simultaneous measurement of heat fluxes and temperatures through the different faces. Note that data acquisition is facilitated by a PC running LabVIEW for further analysis.

According to the thermal program applied to the sample, it is possible to characterize the apparent thermal conductivity, the specific heat capacity, the phase change temperature and latent heat of the sample. Indeed, the experimental set-up permits to impose temperature loading on each face of the studied sample. Simultaneous measurements of temperatures and heat fluxes exchanged during heating and cooling processes allow the determination of the apparent thermo-physical properties. The characterization also yields to the determination of conductivities and heat capacities when the PCM is in solid or liquid state. During the phase change, temperature and the latent heat can be determined via proven methods [50 -52], to name few.

Obtaining Octadecane thermo-physical properties

Thermal conductivity and thermal contact resistance in solid phase

For the solid phase test, we impose, using the previous experimental setup, temperatures of TG = 10 °C and TD = 15 °C (below melting temperature) on the left and right faces, respectively, until a thermal steady state is reached. The heat fluxes are initially constant. At a particular time, an increase of 5°C in temperature is imposed in both thermo-regulated baths linked to the plate heat exchangers. This will induce an increase of temperature in the material, before reaching again a thermal equilibrium. The sample would thus store energy (sensible heat) between these two permanent steady states. A similar experiment is performed to determine the thermal conductivity in the liquid phase. The sample is also subjected to a temperature increase of 5°C starting from TG = 30 °C and TD = 35 °C, insuring temperatures higher than the melting temperature of the Octadecane. The subjected thermal loading ramp and the resulting heat fluxes are shown in Fig.

4(a, b).

The thermal conductivity λ can be determined depending on the sum of heat fluxes and the temperature difference, using the following relation [START_REF] Zalewski | A novel technique for experimental thermophysical characterization of phase-change materials[END_REF],

2 L / T λ ϕ = ∆ ∑ ( 16
)
where φ is the heat flux density in W/m 2 and L is the thickness of the sample.

It is useful to mention that the used method involved simultaneous measurements of the heat fluxes and temperatures on both sides of the differentially heated sample. Compared to conventional methods of determining the thermal conductivity in steady state, this method is four times faster [START_REF] Zalewski | A novel technique for experimental thermophysical characterization of phase-change materials[END_REF]. The exact term for the quantity measured is "thermal transmission", which, depending on the material being measured, can have components of convective, radiative and conductive heat transfer; it is commonly referred to as the effective or apparent thermal conductivity.

The heat sensor is made of a thin material with a high thermal conductivity (copper). Thermal contact resistance (or contact resistivity) may cause huge errors of thermal conductivity measurements if it is not taken into account. Moreover, when there are two different material layers, a thermal resistance exists on their interface. It is more obvious in the case of two solids. In case of thermal insulation materials (small conductivity) the sample's thermal resistance is large and thermal contact resistance can be neglected. But, in case of higher conductivity materials (> 0.1 W/mK), the thermal contact resistance becomes significant compared to the sample's thermal resistance and cannot be neglected [START_REF] Tleoubaev | Fast measurements of absolute thermal conductivity excluding thermal contact resistance errors[END_REF]. Hence, it is important to estimate this effect for accurate results.

The resistance of the sample can be defined by:

eq c pl oct R 2R 2R R = + + ( 17
)
where Rc, Rpl, and Roct, are respectively the thermal resistivity contact, resistance of Plexiglas and Octadecane. We first neglect the thermal contact resistivity. Then,

pl oct eq eq pl oct 2 L L L R λ λ λ = = + (18)
Here, we apply Eq. 18 to calculate the equivalent conductivity of the sample (Plexiglas + Octadecane). The literature properties of Plexiglas are gathered in Table 1, and those for Octadecane are in Table 2. The resultant conductivity in the solid case is: λeq = 0.167 W/(m.K).

Hence, by applying Eq. 18, we get λoct = 0.152 W/(m.K). This value is very far from the conductivity of Octadecane reported in the literature [START_REF] Messerly | Low-temperature thermal data for pentane, n-heptadecane, and n-octadecane. Revised thermodynamic functions for the nalkanes, C5-C18[END_REF][START_REF] Mehling | Heat and cold storage with PCM: An up to date introduction into basics and applications[END_REF]. The error is probably due to ignoring the thermal contact resistivity. As stated before, when testing samples of moderate thermal conductivity (λ~ 0.1-10 W/(m.K) the thermal contact (or interface) resistance must be addressed, otherwise significant errors will result. Since we are not sure about Rc, we will suppose range of Plexiglas and solid PCM. This is applicable for the 4 faces of the enclosure. However, since we take a 2D model, we ignore the effect of resistance on the upper and lower faces.

Thermal conductivity, thermal contact resistance and convection in liquid phase

As in the previous procedure, we calculate the conductivity in liquid Octadecane. The Flux density and surface temperatures are plotted in Fig 4b . The estimated value by Eq. 18 is λapp (liquid) = 0.22 W/(m.K), which can be seen as an apparent conductivity. It is important to recall here that the quantity measured may have a convection component. Therefore, we seek to, theoretically, prove the existence of convection in the liquid state. This would check the high value of the apparent conductivity in liquid phase, and then to confirm this experimentally and numerically.

Theoretical evaluation of the convection

We need to calculate RaH = βgΔTH 3 /αυ for the above problem (A=H/L ≈ 20), corresponding to the considered sample. This gives us an idea about the mechanism of heat transfer, within the sample, and on the effect of the fluid flow, following the scaling laws [START_REF] Faghri | Advanced Heat and Mass Transfer[END_REF]. In such a tall enclosure the convection is usually ignored. Here, our aim is to check the importance of considering convection theoretically, numerically, and experimentally.

According to the scaling laws [START_REF] Faghri | Advanced Heat and Mass Transfer[END_REF], the parameter RaH needs to be calculated and to H/L:

(RaH) 1/4 = (βgΔTH 3 /να) 1/4 = (1.415×10 10 ) 1/4 ≈ 345. (RaH) -1/4 = (βgΔTH 3 /να) -1/4 = (1.415×10 10 ) -1/4 ≈ 0.003 H/L = 20.2, so (RaH) -1/4 < H/L < (RaH) 1/4 .
So, the flow pattern is expected to show a boundary layer on all four walls, with one core remaining stagnant. The dominant heat transfer mechanism is "boundary layer convection" with a significant effect of the fluid flow. To calculate the Nusselt number, we can use an adequate empirical formula. Since RaL = βgΔTH 3 /να = 3.43×10 5 , with 10< H/L < 40, 1 <Pr < 2×10 4 , 10 4 <RaL <10 7 , then:

( )

L 0.3 1 / 4 0.012 L Nu 0.42 xRa x Pr x H / L - = (19)
This relationship allows to deduce with RaL = βgΔTL 3 /να. This reveals that the convective heat transfer is four times stronger than the conductive one. Once again, this shows that the convection is the dominant mechanism in this problem. The apparent conductivity measured by Eq. 18 has a convection contribution, and we can state that it is not the equivalent conductivity of the sample.

Experimental investigation of the convection

To achieve this, we perform a new experiment. We inject in the Plexiglas enclosure of Octadecane, described in Fig. 2, silver coated hollow glass spheres, and allow them to settle down.

Using Particle Image Velocimetry (PIV) system, we trace the flow of the particles inside the enclosure. In the case of liquid Octadecane, due to the presence of convection, and as expected, the L Nu 4.33 = 20 powder circulates inside liquid Octadecane. In contrary, if conduction was the dominant mode, the powder would have stayed at the bottom, and no motion would have been noticed.

After proving the presence of convection, we will estimate the value of the conductivity of liquid Octadecane, by the aid of the numerical simulations.

Numerical assessment of the convection

The implemented LBM numerical model is used to evaluate the convection in liquid Octadecane within the Plexiglas enclosure. Thus, the average Nusselt number ( L Nu ) value on the Octadecane interface is 4.81. Such a value seems close to that estimated empirically by the order of 10%.

The average equivalent conductivity for natural convection can be estimated via the following relationship [START_REF] Faghri | Advanced Heat and Mass Transfer[END_REF]:

app LBM ,conv eq LBM ,cond Q Q λ λ ≈ (20)
where λapp and λeq are the apparent and equivalent conductivities of liquid Octadecane. Note that, by apparent here, we mean that it contains the convection contribution. As for the equivalent, we point to the conductivity calculated for the sample (Plexiglas + Octadecane). To compute the ratio, we implement two models. The first is for convection (with a Boussinesq force evaluated by the corresponding Ra value), and the second is a pure conduction model (without the Boussinesq force, insuring a zero-velocity field). The simulation results do not depend on the selected value of conductivity, since in LBM we assign a fixed Prandtl number of 50. The main aim here is to evaluate the effect of including the Boussinesq force that designates the convection. This method is only an approximation. The calculated ratio is λapp/λeq The conductivity of liquid Octadecane is λoct = 0.148 W/(m.K) in literature. So in both liquid and solid cases the closest calculated value of conductivity is obtained when considering Rc = 1/60.

Fig. 5b shows the estimated values of λoct for the liquid Octadecane versus the thermal contact resistance. As deduced, the effect of thermal resistivity contact is much more in the case of the solid state, where there is a solid-solid contact. Thereby an interface of the air gap may rise in between. From our approach, the conductivity of solid and liquid Octadecane are λoct (solid) = 0.31 W/(m.K) and λoct (liquid) = 0.151 W/(m.K), respectively with Rc = 1/60. However, we should still consider border effects for more accuracy.

Specific heat capacity

The temperatures and the heat fluxes evolution on both sides of the sample, when Octadecane is in the solid and the liquid phases are represented in Fig. 4. A symmetrical behavior of heat fluxes and temperatures measured on both faces of the sample can be observed, which correspond to the results classically obtained with a solid material without phase change.

Initially, the sample is isothermal (at a thermal equilibrium state), and then a temperature change is imposed on its both sides. This induces a thermal evolution of the system in an asymptotic way towards a second thermal steady state. Also, it should be noted that the heat fluxes evolve quickly when the temperatures of the sample are increased and then converge to a second thermal equilibrium state obtained at the end of the test. The specific heat capacities of the sample are obtained starting from the determination of the sensible heat accumulated by the material between the imposed temperatures. The sensible heat is calculated by integrating the difference in heat fluxes between the initial and the final state using the following relationship:

( )

1 end init T p,i end init T i Q C T T L ϕ ρ = ∆ = - ∑ ∫ (21)
where Δφ represents the cumulative heat rate entering the sample and Cp is the apparent specific heat capacity (kJ/(kg.°C)).

The calculated specific heat capacity of Octadecane in solid and liquid phases are: coct,solid = 1925 kJ/(kg.°C) and coct,liquid = 2365 kJ/(kg.°C). It is noted that these values are consistent with those available in the literature.

Analysis of the heat flux with phase-change

To analyze the heat flux in the presence of a phase change, we use a four-hour ramp. The sample is initially at 15.8 °C, then the both plates' temperature Tp is increased from 15.8 °C to 40 °C in 4 hours (4h). Then a constant loading of 40 °C is applied followed by a decreasing ramp until reaching 15.8 °C again, as illustrated in Fig. 6. As noted, the sample is initially in thermal equilibrium, where the heat flux is zero. This is followed by sensible heat flow when the PCM is still in its solid state. The curve then starts increasing gradually at a time corresponding to the onset of melting. The latent heat absorbed also increases gradually and the largest quantity is absorbed in the vicinity of the melting point. After complete melting, and reaching the constant heating phase, the PCM goes back to a thermal equilibrium, with no more flux variation, until the cooling ramp starts.

When cooling starts, the liquid PCM releases sensible heat until reaching the solidification temperature. The difference between the melting and solidification is that the majority of freezing happens directly, and the largest amount of latent heat is released at first, then the rest will solidify gradually. This behavior can be also regarded to the expected shape of the enthalpy curve. The portion around the curve at the solidification temperature is discontinuous. There is a direct energy jump from a temperature point in the liquid state to the solidification point. We note also the existence of supercooling, which is illustrated by the formation of a vertical portion or discontinuity in the curve (Fig. 1). After complete freezing, the PCM returns back to its equilibrium state with almost zero flux.

Results and discussion

As previously stated, the major aim is to first estimate the main thermo-physical properties and to highlight their effect on the numerical results if the supercooling, convection, and soluble impurities are ignored or not. To achieve this, we compare the numerical results, using an improved model, with the experimental ones. The properties in Tables 1 and2 and the verified thermophysical properties (conductivities, specific heat capacities, latent heat of fusion, thermal contact resistance, etc.) are fixed in the model. It should be noted that the Boussinesq force has been varied to account for convection, the melting temperature to calibrate the enthalpy curve, and the degree of supercooling to estimate the nucleation temperature. In the following, some significant results showing effects of physical phenomena are comprehensively presented and discussed.

Effect of enthalpy curve shape

As explained before, the presence of soluble impurities in the PCM results in a melting temperature depression. Recall that, the used Octadecane is thought to be 99% pure. So, the 1% of soluble impurities can have an impact on the heat flux curve. This will be investigated theoretically and numerically. First, the depression in the melting temperature is calculated from Eq. ( 1) to obtain Tm = 28.03 °C. Thus, the PCM melts with a depression of about 0.1 °C compared to the case of pure Octadecane. Then, this melting temperature will be used in the numerical model, to check if this the numerical and experimental fluxes match. If not, Tm will be changed until they well match. The enthalpy-temperature curves are plotted in Fig. 7 for different values of (Ta -Tm). As demonstrated, the curve becomes less steep as Tm deviates from Ta. This induces a melting temperature range that causes a gradual increase in latent heat instead of absorbing it at once. Moreover, to show the important impact of enthalpy curve on the flux, we compare the heat flux obtained numerically for the case pure Octadecane to the experimental one (cf. Fig. 8a). In pure PCM, there is no melting range. Thus, the immediate absorption of latent heat at the melting temperature results in a sharp and high rate of heat flow at this moment, as shown in the comparison curves. However, when there exists a temperature depression, the rate of heat absorption decreases more, since there is no instant absorption and release. The accurate prediction of the onset of melting and freeing is highly dependent on taking into account the real enthalpy-temperature relation.

By comparing the experimental and numerical curves, the best fit was obtained for Tm = 27.6 °C, which corresponds to a variation of about 0.5°C from the supposed pure melting temperature found in literature (Fig. 7). This shows that using Eq. ( 1) may not be adequate in our case but can be used as pre-estimation. Also, the supposed melting temperature [START_REF] Messerly | Low-temperature thermal data for pentane, n-heptadecane, and n-octadecane. Revised thermodynamic functions for the nalkanes, C5-C18[END_REF][START_REF] Mehling | Heat and cold storage with PCM: An up to date introduction into basics and applications[END_REF] of pure Octadecane is also questionable. If Octadecane is considered pure, the heat flux curve will show a high discrepancy when compared to the experimental one, as shown in Fig. 8a. The curve obtained numerically comes from considering the enthalpy curve as the traditional piece-wise function. This shows that such a curve cannot be used for our case.

The final numerical curve for melting is plotted in Fig. 11 despite the fact that this is the result from the enhanced numerical model, i.e. including also convection. In addition, we highlight that the heating rate does not influence the prediction of the melting point Tm, since the enthalpy curve is unique for a material.

Effect of presence of supercooling and convection

We show, in what follows, the discrepancy that would have resulted if we ignored the presence of convection and supercooling. The major effect of ignoring convection (cf. Fig. 9a) is on estimating the time of melting and mainly that of solidification. At first, the material is solid. So convection will start affecting the curve as soon as the liquid fraction starts increasing. It is important here to address the effectiveness of our partial bounce-back LBM model. By this method, we are able to treat the mushy zone as a porous medium, where the velocity of the fluid is restrained but not completely blocked. The Rayleigh number (Ra), corresponding to the real convection, is searched by a fitting method, and is found to be Ra = 1.2×10 4 . Higher values would lead to an over estimation (Fig. 9b).

On the other hand, if supercooling is ignored (see Fig. 10a), the solidification starts before the real instant. So, the experimental and numerical results will show a lag in time. To overcome this, we have sought the degree of supercooling through a tuning (fitting) method. Hereafter, the best estimated degree of supercooling corresponds to 1.4 °C. If higher supercooling degree is used, the supercooling phenomenon will be overestimated (see Fig. 10b).

By this, we can combine all the enhancements to the present model: using the real enthalpy curve, including convection and supercooling, using variable thermo-physical properties for solid and liquid phases, and using a partial bounce-back LBM for phase change to better simulate the low effect of convection. The enhanced model calculated for Ra = 1.2×10 4 , degree of super cooling of 1.4 °C, and a melting point depression of 0.43 °C is displayed in Fig. 11, and shows a great improvement, and matches the experimental results.

The simplifications taken in some numerical models can result in an overall discrepancy in the real behavior of PCM. These discrepancies may lead to wrong estimation of the fusion times and the energy amount stored. Consequently, the PCM will not give the desired performance. The major enhancement to the proposed model was to take into account natural convection, supercooling, to adopt variable thermo-physical properties for solid and liquid phases, and to use the real enthalpy curves corresponding to the considered PCM. These improvements have led to better agreement with the experimental results, and showed that if the above phenomena were present, they cannot be ignored.

Conclusion

The main aim of this work is to characterize the thermo-physical properties and performance of Octadecane. This PCM may be promising for the usage in glass facades or other building materials and latent heat systems. To achieve such an objective, we adopted an approach based on experimental and numerical techniques. The enhanced numerical model, used here, can help estimate some properties, and can be also used as an optimization tool for any practical application.

The thermo-physical properties obtained experimentally and, using the numerical model, match the expected values from literature. The adopted numerical model is a thermal LBM with a partial

  being the lattice sound speed. It should be noted that, the lattice viscosity and thermal diffusivity are selected so as to conform to the intended Prandtl number Pr the additional force term related to Boussinesq force Fb is incorporated in the model by shifting the velocity field by a term of b f

  ) , w x v being the computed near-wall velocity, and p g represents a known distribution function. Hence, the unknown DFs p g are computed by . As for the adiabatic BCs, the Neumann BCs are achieved using the BB boundary conditions for the distribution i g , as prescribed for i f .

values of 1 /

 1 Rc = 50 -100. Thus, as we change Rc , we obtain different thermal conductivities of solid Octadecane: (Rc = 1/100; λoct = 0.219 W/(m.K); (Rc = 1/70; λoct = 0.27 W/(m.K); (Rc = 1/60; λoct = 0.31 W/(m.K); (Rc = 1/50; λoct = 0.39 W/(m.K). From here, we can realize the major impact of considering the thermal contact in calculating the conductivity of Octadecane, which has a relatively moderate conductivity. The estimated values of λoct are plotted versus Rc in Fig 5a. The value in the literature for the conductivity of Octadecane is λoct (solid) = 0.356 W/(m.K) which corresponds mostly to around Rc = 1/60. It is worth mentioning that the contact resistances here occur between the flux sensors and Plexiglas, and also between

  = 1.76, hence the value of the equivalent conductivity is λeq = 0.125 W/mK. Here, the thermal contact resistance occurs only between the flux sensors and the Plexiglas. Now we use Eq. 18 with different trials on the thermal contact resistance to get: Rc = 1/100; λoct = 0.123 W/(m.K); Rc = 1/70, λoct = 0.137 W/(m.K);Rc = 1/60; λoct = 0.151 W/(m.K); Rc = 1/50; λoct = 0.162 W/(m.K).

Figure 1 .

 1 Figure 1. Enthalpy-temperature relationship for pure PCMs with soluble impurities in the presence of supercooling.

Figure 2 .

 2 Figure 2. Two-dimensional schematic of the studied sample.

Figure 3 .

 3 Figure 3. (a) Schematic illustration of the experimental setup, (b) Schematic illustration of the studied sample.

Figure 4 .

 4 Figure 4. Heat flux curves with no phase-change for (a) solid state, (b) liquid state.

Figure 5 .

 5 Figure 5. Estimated thermal conductivity vs thermal contact resistivity: (a) solid state, (b) liquid state.

Figure 6 .

 6 Figure 6. Heat flux curves indicating the presence of supercooling.

Figure 7 .

 7 Figure 7. Enthalpy-temperature relationship for different (Ta -Tm) values.

Figure 8 .

 8 Figure 8. (a) Effect of the presence of impurities on the heat flux, (b) Numerical and experimental shapes with Tm = 27.6 °C.

Figure 9 .Figure 10 .Figure 11 . 35 Figure 1 . 2 . 36 Figure 3a .Figure 4 .Figure 5 . 37 Figure 6 .

 910113512363a45376 Figure 9. Effect of convection on the shape of the heat flux curve: (a) ignoring convection, (b) overestimating convection Figure 10. Effect of supercooling on the shape of the heat flux curve: (a) ignoring supercooling, (b) overestimating supercooling Figure 11. Experimental vs. numerical results when using the enhanced model.

Figure 7 .Figure 8 .Figure 9 .Figure 10 .

 78910 Figure 7. Enthalpy-temperature relationship for different (Ta -Tm) values.

bounce-back approach and an enthalpy-based model. We considered the natural convection, supercooling, variable thermo-physical properties (when necessary), and the real enthalpy curve of Octadecane. Simulations performed using the improved model corroborate the experimental results and show that ignoring these phenomena may lead to wrong estimation of the fusion times and amount of energy stored and the PCM will not give the desired performance.

Disclosure statement

The authors declare no potential conflicts of interest regarding authorship and/or publication of this paper. Table 1. Properties of Plexiglas [START_REF] Messerly | Low-temperature thermal data for pentane, n-heptadecane, and n-octadecane. Revised thermodynamic functions for the nalkanes, C5-C18[END_REF][START_REF] Mehling | Heat and cold storage with PCM: An up to date introduction into basics and applications[END_REF] Property Value Heat capacity 1470 J.Kg -

List of Tables


Moudhaffar Nouri
email: moudhaffar.nouri@centralesupelec.fr

Rihab Hamila

Patrick Perre

A double distribution lattice Boltzmann scheme for unsteady Conjugate Heat Transfer : the DD-CHT LB method

In this paper, we propose an innovative formulation of the "well-known" thermal lattice Boltzmann method for solving transient conduction heat transfer problems in heterogeneous media : the DD-CHT LB method for Double Distribution -Conjugate Heat Transfer. The main idea of this new formulation is to bring out a second distribution function that represent the effect of thermal inertia and its variation in a heterogeneous media. By adopting this original formulation we achieved several goals, mainly conservation of local properties of the studied media and ensuring the stability of the proposed formulation by avoiding any alteration of the "well-known" TLBM collision step. Validation of the formulation is fulfilled through three different benchmark conduction heat transfer problems in heterogeneous media. All results are compared with those obtained by the Control Volume Method (CVM) and excellent agreements were found.

Introduction

Conjugate Heat Transfer (CHT) takes place when heat is transferred between two materials with different transport properties. The major difficulty when solving CHT is the necessity of satisfying two fundamental conditions: the temperature continuity and the normal heat flux continuity.

Generally when adopting conventional Computational Fluid Dynamic (CFD) solvers in solving CHT problems, extrapolations and iterative schemes have to be applied. As a result extra computational cost and complex codes are to be expected.

In earlier researches trying to solve CHT with the LB method, most of the proposed solutions are based on solving the convection-diffusion equation in heterogeneous media based on temperature as a dependent variable [START_REF] Chen | A note on the solution of conjugate heat transfer problems using simple-like algorithms[END_REF][START_REF] Meng | Lattice boltzmann simulations of conjugate heat transfer in high-frequency oscillating flows[END_REF][START_REF] Hu | Simulation of steady fluid-solid conjugate heat transfer problems via immersed boundary-lattice boltzmann method[END_REF]. In these works precautions must be taken regarding tracking and correcting for the presence of interfaces. By doing so, the overall solution procedure become more complex especially in case of curved interfaces, moving or evolving interfaces over time. [START_REF] Wang | Counter-extrapolation method for conjugate heat and mass transfer with interfacial discontinuity[END_REF] extended a previous work of [START_REF] Le | Counter-extrapolation method for conjugate interfaces in computational heat and mass transfer[END_REF] in order to include the interface discontinuity. By presenting their model they proposed a solution for treatment of discontinuity at interfaces. Their solution is opting for the Counterextrapolation method which provides extrapolations of the computed transport scalar along the local normal direction of the conjugate interface and combining that with bilinear interpolations. [START_REF] Patel | A spectral-element discontinuous galerkin thermal lattice boltzmann method for conjugate heat transfer applications[END_REF] Another attempt to solve CHT problems with the LB method was performed by Hamila et al. [START_REF] Hamila | An enthalpy-based lattice boltzmann formulation for unsteady convection-diffusion heat transfer problems in heterogeneous media[END_REF]. In their work, the authors proposed a general formulation of the convection-diffusion equation in multicomponent media at its enthalpic form. The main idea of the their work is to introduce an extra source term to the enthalpic form of the con-vection diffusion equation. The same idea was also tested in the three dimensional case [START_REF] Nouri | A three-dimensional enthalpic lattice boltzmann formulation for convection-diffusion heat transfer problems in heterogeneous media[END_REF].

In this paper we propose a completely different approach from what we mentioned above.

The basic idea of the proposed formulation is a double distribution approach of thermal transfer, in other words the concept is an internal reformulation of the lattice Boltzmann equation. Without any further calculations of term source or specific interface tracking method, the proposed model aim to solve CHT problems in heterogeneous media with local properties conservation and stability of the model proposed.

Mathematical formulation

Thermal Lattice Boltzmann formulation

The Thermal Lattice Boltzmann formulation for heat transfer problems based on BGK collision operator takes the following form:

f k (r + e k δt, t + δt) = f k (r, t) - δt τ [f k (r, t) -f (eq) k (r, t)] (1)
where f k is the distribution function, denoting the number of particles at lattice node r and time t, τ is the relaxation time, f (eq) k

is the equilibrium distribution function, and δt is the time step. In two dimensional problems, when using the D2Q9 lattice scheme τ and f (eq) i are as below :

τ = 2α |e k | 2 + δt 2 (2) f (eq) k = w k T (r, t) 1 + e k • U c 2 s ( 3
)
where α = k ρCp is the thermal diffusivity, k and ρCp are thermal conductivity and heat capacitance, respectively. The nine velocities e k and their corresponding weights w k in the D2Q9 lattice scheme are the following:

e 0 = (0, 0) • U e 1,3 = (±1, 0) • U e 2,4 = (0, ±1) • U ; e 5,6,7,8 = (±1, ±1) • U (4) w 0 = 4 9 , w 1,2,3,4 = 1 9 w 5,6,7,8 = 1 36 (5)
3 once the values of population f k are known, temperature T is calculated as follows:

T ( r, t) = k f k ( r, t) (6)
Equation 1 based on Chapmann-Ensgok expansion and BGK collision model recover the following macroscopic equation

∂T ∂t = ∇ • (α∇T ) (7)
When analyzing the macroscopic energy equation recovred by thermal lattice Boltzmann formulation, we notice that the method can only provide correct results in homogenous media or in heterogenous media with identical heat capacity in either sides of the existing interfaces [START_REF] Hamila | An enthalpy-based lattice boltzmann formulation for unsteady convection-diffusion heat transfer problems in heterogeneous media[END_REF][START_REF] Karani | Lattice boltzmann formulation for conjugate heat transfer in heterogeneous media[END_REF]. Indeed, as the distribution function is used to compute both the flux between points of the lattice and the temperature via equation ( 6), temperature and enthalpy are equal in this approach, which implicitly imposes ρCp = 1.

The double distribution lattice Boltzmann equation for unsteady heat conduction in heterogeneous media

In this paper we propose a more generalized formulation of the thermal LBM. The formulation allows to distinguish temperature and enthalpy when solving heat transfer problems.

To do so, a second distribution function is added. This function is very simple to compute as it contains only one component. Through these changes we will be able to correctly simulate thermal transfer in heterogeneous media without further considerations of heat flux and temperature continuities at the interface.

Writing the energy equation in the enthalpy form is as following:

∂h ∂t = ∇ • (k∇T ) (8)
with h = ρCpT , and k represent the conductivity.

By operating a slight change on Eq. 8 , we aim to express this equation with temperature as a dependent variable. So by adding in both sides of equation ∂T ∂t the equation takes the following form:

∂(ρCpT -T ) ∂t II + ∂(T ) ∂t = ∇ • (k∇T ) I (9)
Eq. 9 can be seen as the sum of two parts. Part I resolves the energy equation for a reference medium in which the volumetric heat capacity is equal to unity ρCp = 1 and Part II account for the actual thermal inertia.

In order to distinguish temperature and enthalpy, we introduce a second distribution g which contains a unique component g 0 with a zero velocity. g is defined to obtain the enthalpy :

T = f i h = ρCpT = g 0 + f i (10)
These expressions tell us that a simple relation relies g 0 and f i :

g 0 = (ρCp -1) f i (11)
At each time-iteration, from t to t + 1, solving Part I using Eq. 1 accounts for the right heat flux, driven by the temperature field, but estimates a wrong temperature change as this equation assumes the heat capacity to be equal to unit ρCp = 1. Let's name the new temperature fields obtained after the streaming-collision steps T * (t+1) = f * (t+1) i

. At this stage, g didn't change as the single component of this distribution is at rest and g is not involved in the collision step. For this reason, the set g

(t)
0 and f * (t+1) i

do not respect equation 13. A correction is therefore needed at constant enthalpy :

h (t+1) = g (t) 0 + i f * (t+1) i = g (t+1) 0 + i f (t+1) i (12)
Together with the definition of the enthalpy :

g (t+1) 0 = (ρCp -1) i f (t+1) i (13) 5
This leads to the following expression :

T (t+1) = β × T * (t+1) (14)
Where :

β = 1 + g (t) 0 /T * (t+1) ρC p ( 15
)
This correction is applied to each component of distribution f (Fig. 1):

f (t+1) i = β × f * (t+1) i (16)
The value of g 0 at t + 1, g (t+1) 0 is then updated using equation 12. The algorithm to be implemented is detailed hereafter :

Algorithm 1 The DD-CHT LBM algorithm Initialize Lattice Boltzmann distribution f (t=0) i Initialize thermal inertia g (t=0) 0 ← (ρCp -1) × i f (t=0) i loop Streaming step 1) Collision step 1) Correct temperature distribution f (t+1) i ← β × f * (t+1) i Update inertia distribution g (t+1) 0 ← g (t) 0 + T * × (1 -β) end loop
1) Boundary nodes need the same treatment as for the classical LBM.

Numerical simulations

Transient heat conduction in a two-layer medium with inclined interface

To demonstrate the ability of the proposed formulation to solve conjugate conduction heat transfer problems in heterogeneous media, a transient heat transfer problem in a two-layer medium with an inclined interface is proposed first. Initially the medium is uniformly at the initial temperature T int . At (t>0), a Dirichlet condition at the hot temperature T 0 is applied to all sides (Fig. 2a). Simulation is carried out on a (250 × 250) mesh. The same problem was solved using the control volume method for validation purposes. Different sets of thermophysical properties were tested for the two sub-domains 1 and 2 to obtain contrasted values of the thermal diffusivity either due to the heat capacity or to the thermal conductivity. The values used are : ρCp 1 = 20, ρCp 2 = 1, k 1 = 1, k 2 = 6 for the first simulation (Fig. 2b) and ρCp 1 = ρCp 2 = 2, k 1 = 1, k 2 = 4 for the second simulation (Fig. 2c). The results obtained with the new LBM formulation and the CVM solutions are depicted as isotherms after 30 seconds. We also plotted in (Fig. 2d) and (Fig. 2e)the dimensionless temperature profile along the diagonal for a set of selected times for case 1 and case 2, respectively. Excellent agreements were found for the two case studies.

Unsteady conduction heat transfer problem inside a square medium with a circular inclusion

The problem investigated in this part of the paper is a 2D unsteady conduction heat transfer problem inside a square medium with circular inclusion. The two mediums have different values of thermophysical properties. Initially the temperature of the entire medium is set to T ini . At (t>0), the left-hand side of the square medium is heated to T 0 while the other sides are kept insulated. In this test case we choose to perform two different simulations, in other word in the first simulation we imposed similar thermal conductivity values of the sub-domains and imposed non unitary ratio between the heat capacity's values, while in the second simulation we chose to impose non unitary ration between the heat capacity's values and between the thermal conductivity's values (Fig. 3a) . Both simulations were carried out on a (150 × 150) mesh.

(Fig. 3b -3c) present the results as computed with the new LBM method for the first simulation case. In this figure, the isotherm lines are compared with the same solution computed with a the control volumes method(CVM). We can see that excellent agreement are found both in transient period (750 s) and at steady-state.

In the second configuration, the conduction heat transfer problem is solved in a square medium with a circular inclusion, the values of thermophysical properties are k 1 /k 2 = 20 and ρCp 2 /ρCp 1 = 20. We intentionally selected the same ratio of properties to obtain the same thermal diffusivity even though the medium is heterogeneous. The first validation for this set of thermophysical properties is plotted as isotherms at t=100s (Fig. 3d).

  presented a spectral element discontinious Galerkin thermal LB model for fluid-solid conjugate heat transfer problems. The authors employed a double distribution thermal LB model for solving flow field and CHT problems on unstructured non-uniform grids. The results presented are mostly at the steady state. Overall comparisons are made with finite difference and continuous Spectral elements solutions. Recently Pareschi et al. (Pareschi et al., 2016) proposed a CHT model based on the two populations entropic LB method. In their study, the authors chose the entropic multi-relaxation LB method and validated their formulation by solving a pure conduction transfer problem, a CHT problem with a backward facing step flow and finally tested their model with a three dimensional turbulent flow around a heated mounted cube. Some other researchers opted for an enthalpy based LB formulation to recover the convectiondiffusion equation under its enthalpic form. Lately Huang et al. (Huang and Wu, 2016) proposed a total enthalpy based MRT LB model associated with an adaptive mesh refinement near the solid-liquid phase change interface where thermodynamic properties vary significantly. Several melting problems in a square cavity with different values of Prandtl number and Rayleigh number were simulated. More recently Liu et al. (Liu et al., 2017) simulated solid-liquid phase change heat transfer problem in metal foam. Their formulation comes under the form of an enthalpic MRT LB model used with the local thermal non-equilibrium (LTNE) condition. The authors implicitly tracked the moving liquid-solid interface by mean of the liquid fraction expression.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Modified distribution functions with the DD-CHT LBM algorithm.

Acknowledgments

The authors would like to thank Department of Marne, Great East Region and Grand Reims for their financial support

Conclusion

In this paper, a rigorous reformulation of the "well-known" thermal lattice Boltzmann method is detailed. The proposed approach aim to correctly simulate thermal transfer in heterogeneous media which obviously mean dealing with both different conductivities and heat capacities on both sides of the interface. Mainly the basic idea of the paper involve a second distribution function, with a single component, that allows distinguishing temperature and enthalpy in the TLBM algorithm. The new formulation offer an innovative method for treatment of unsteady conduction heat transfer problems in heterogeneous media. The method is presented as a novel LBM algorithm which is easy to implement and works without need to track the interface of the heterogeneities. After showing the ability of the presented formulation to simulate conduction heat transfer in homogeneous media, the method has been tested and validated with three unsteady conduction heat transfer problems in heterogeneous media.