Aina Rakotondrandisa

Je profite également de ce moment pour remercier le projet M2NUM 1 pour le financement des trois ans de travaux de thèse. Cela m'a offert l'opportunité d'assister à de nombreux congrès et conférences nationaux et internationaux. Merci également au centre de calcul régional CRIANN 2 d'avoir accordé des heures de calculs, indispensables notamment pour les calculs tri-dimensionnels.

Je témoigne de ma reconnaissance à Pierre Calka, le directeur du laboratoire LMRS 3 . Sans son accord et sa confiance, ce travail n'aurait pu être mené à son terme. Merci au laboratoire pour le financement de nombreux déplacements pour ne citer que "l'atelier de résolution de problèmes industriels" à Montréal ou encore "l'école de printemps inter-disciplinaire en simulation numérique et optimisations de code" à Lyon. Merci à tout le personnel administratif : Edwige Auvray, Hamed Smail, Sandrine Halé, Marc Jolly, et Saada Benamar, pour leur coopération et leur aide au quotidien. Merci à l'équipe du midi, pour votre bonne humeur et les moments conviviaux. Un grand merci spécialement à toute l'équipe calcul scientifique, ceux qui y sont encore, ceux qui étaient de passage, ou ceux qui sont partis mais que j'ai pu côtoyer durant cette thèse : Francky Luddens, Guillaume Vergez, Lei Bian, Corentin Lothodé, et Geoges Sadaka. Je tiens à mentionner spécialement Georges Sadaka pour le grand travail qu'il a effectué sur la partie 3D du code. J'ai conscience d'avoir eu la chance de travailler auprès de personnes exceptionnelles -pas que sur le plan scientifique.

Je remercie aussi les développeurs du logiciel éléments finis FreeFem++ et des libraires FFDDM et HPDDM, Frédéric Hecht, Pierre-Henry Tournier et Pierre Jolivet, pour toute l'aide dont ils ont fait preuve pour me guider dans l'utilisation du logiciel. J'adresse en particulier mes remerciements à Luminita Danaila, que j'ai eu l'opportunité d'avoir en cours de mécanique des fluides et de turbulence depuis les années de Master à l'Université de Rouen Normandie, pour sa grande pédagogie et sa gentillesse. Je suis reconnaissant pour toute la patience et la compréhension dont elle a fait preuve, notamment durant la rédaction du papier sur le cycle complet fusion-solidification.

Je tiens aussi à exprimer ma gratitude à l'église protestante malgache de Rouen (FPMA) pour leur soutien indéfectible depuis ces sept années que je suis à Rouen. Merci pour l'amour

since it generates high temperatures and pressures. Storing thermal energy is hence an efficient and fundamental way to store the energy. This can be realised by rising a substance's temperature (sensible heat energy storage) or by changing a substance's phase (latent heat energy storage).

To this end, research on passive heat storage systems has generated many studies lastly, namely on latent heat energy storage like Phase Change Material (PCM). PCMs are used to store heat during the melting of the material and release later the stored energy during the solidification process. At the present time, the latent heat storage technologies are proven as an effective solution to decrease the use of fossil fuel and in the same time to increase the efficiency of energy usage.

Aside from energy storage technologies, PCMs are also widely used in building applications, to decrease temperature fluctuations. The latest announcement of the French ministry of ecological transition, detailing the repartition of the energy consumption in different domains, indicates that more than 35% of the total energy is consumed by residential and commercial buildings. Details about energy consumptions in 2017, reported by the French ministry of ecological transition are shown in Tab. More than 60% of the total energy consumption in the residential sector is used to space heating (see Fig. 1.1b). Research towards an energy-efficient building to reduce heating and cooling demands is of principal interest nowadays. Taking advantage of the high value of the latent heat of solid-liquid transformations, PCMs are extensively encountered in building's thermal regulation to reduce overheating. In summer, PCMs are used to absorb the excessive solar radiation heat and maintain a bracing indoor ambience. PCM with a temperature of fusion close to the ambient temperature is generally used to ensure melting during the daytime and solidification during the night-time. During winter however, PCMs can be used to store heat generated by electrical heating system during the night and then release it in the daytime.

PCMs are also encountered in other applications ranging from metal casting and passive temperature control devices for modern portable electronics, to food freezing and cryosurgery. In many of these applications, the choice of an appropriate material for a specific end depends of many criteria, such as the operating temperature range, the thermal conductivity, the cost, etc. They are generally classified into three classes: organic, inorganic and eutectic. A target application area for some PCMs is drawn in Tab. 1.2. The main operating temperature range can be assorted in three groups of melting temperature range. First, 0 to 65 o C for thermal storage used in domestic heating/cooling or for thermal storage of solar energy. Paraffins and water are used for such applications. Second, 80 to 120 o C for the cooling of systems with temperature of fusion less than 120 o C. Finally, temperatures greater than 150 o C for the heat storage in solar power plants based on parabolic trough collectors and direct steam generation.

For more details about these applications, see the review by [START_REF] Agyenim | A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (LHTESS)[END_REF].

Melting and solidification are also fundamental processes in geophysical problem such as Earth's mantle formations, lava lakes [START_REF] Davaille | Thermal convection in lava lakes[END_REF], thermal convection in magma chambers [START_REF] Brandeis | The convective liquidus in a solidifying magma chamber: a fluid dynamic investigation[END_REF] or ice-melt lakes [START_REF] Polashenski | The mechanisms of sea ice melt pond formation and evolution[END_REF]. Ice-melt ponds that form during the summer season in the Arctic are known for example to display natural convection coupled 2 1.2 Theoretical and experimental studies of phase-change systems Storage for solar power plants based on parabolic trough collectors and direct steam generation.

Table 1.2: Target application area for some PCMs studied in the literature [START_REF] Agyenim | A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (LHTESS)[END_REF].

to a phase-change process on the bottom side [START_REF] Polashenski | The mechanisms of sea ice melt pond formation and evolution[END_REF][START_REF] Esfahani | Basal melting driven by turbulent thermal convection[END_REF]. In this case, Rayleigh-Bénard like convection cells are observed in the liquid phase.

The solid-liquid phase-change phenomenon is known to involve complex processes. It couples the natural convection in the liquid phase induced by the buoyancy force to the non-linear evolution of the phase-change interface and the heat transfer process, which could be different from one configuration to another. The coupling of all these physical phenomena induces a strong non-linear process in the solid-liquid problems, difficult to analyse, excepting for simple and ideal test cases. As an example, Fig. 1.2 shows the experimental investigation of the melting of n-octadecane inside a transparent brick, used as a smart material to control the indoor environment of a building [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice Boltzmann method and performed with CUDA enabled GPU[END_REF]. The investigated material is heated from the right and melts consequently from the right to the left. Figs. 1.2a and 1.2b display the evolution of the liquid fraction and the vector field in the fluid obtained by particle image velocimetry (PIV) method. Capturing accurately the non-linear evolution of the solidliquid interface and the strong convection in the melting PCM is clearly a challenging numerical task. Furthermore, the existence of boundary layers near the walls and the interface suggests that the mesh resolution should adapt to capture these structures.

Theoretical and experimental studies of phase-change systems

Solid-liquid phase-change problems have attracted attention since the study of crust formation of earth by [START_REF] Lamé | Mémoire sur la solidification par refroidissement d'un globe liquide[END_REF] or the ice formation problem by [START_REF] Stefan | Uber einige probleme der theorie der warmeletung[END_REF]. From a mathematical point of view, the phase transition is referred to a moving boundary problem. The Stefan problem considers a linear evolution of the interface by assuming only conductive heat transfer. Under this as-Figure 1.2: Experimental observation of the melting of n-octadecane within a transparent brick of plexiglas heated vertically from the right: (a) melted fraction (b) velocity field in the liquid region obtained by particle image velocimetry (PIV) method [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice Boltzmann method and performed with CUDA enabled GPU[END_REF].

sumption, the existence and the uniqueness of the solution for the one-dimensional problem was proven by [START_REF] Rubinstein | On the solution of Stefan's problem[END_REF] using constant properties. For a historical review detailing the evolution of the Stefan problem and the progressive consideration of convection in the liquid, one can refer to [START_REF] Yao | Melting and freezing[END_REF]. Over the past 30 years, solid-liquid phase-change involving natural convection has been widely studied. Rectangular and cylindrical geometries were the most investigated configurations. [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]; [START_REF] Gau | Flow visualization during solid-liquid phase change heat transfer II. Melting in a rectangular cavity[END_REF]; [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice Boltzmann method and performed with CUDA enabled GPU[END_REF] have investigated experimentally the melting of noctadecane and gallium PCMs inside rectangular containers. [START_REF] Ho | Inward solid-liquid phase-change heat transfer in a rectangular cavity with conducting vertical walls[END_REF]; [START_REF] Liu | Experimental study of the phase change heat transfer inside a horizontal cylindrical latent heat energy storage system[END_REF]; [START_REF] Omojaro | Study on solid liquid interface heat transfer of PCM under simultaneous charging and discharging (SCD) in horizontal cylinder annulus[END_REF] have studied experimentally the melting of different paraffins, within cylindrical enclosures. These experimental investigations have been extensively used for numerical validations (see [START_REF] Bertrand | Melting driven by natural convection A comparison exercise: first results[END_REF]; [START_REF] Gobin | Melting from an isothermal vertical wall. Synthesis of numerical comparison exercise[END_REF]; Wang et al. (2010b); [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]; [START_REF] Rakotondrandisa | Numerical modelling of a melting-solidification cycle of a phase-change material with complete or partial melting[END_REF]) and have permitted a better comprehension of the physical phenomena that occur during melting and solidification, mainly the heat transfer process. The experimental work of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] has permitted for example to formulate correlations for the variation of dimensionless thermal energy stored as latent heat. Later, [START_REF] Ho | Inward solid-liquid phase-change heat transfer in a rectangular cavity with conducting vertical walls[END_REF] have analyzed the solid-liquid interface position for different configurations and [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] have suggested Nu-Ra correlations through a scaling analysis. [START_REF] Bejan | Analysis of melting by natural convection in an enclosure[END_REF] gathered the previous observations and described analytically the solid-liquid melting process in the case of a vertical heating. Le Quéré and Gobin (1999) have performed instability analysis during the melting of low and large Prandtl numbers. Gallium and n-octadecane were the most investigated materials. Besides, their physical properties are equal in both liquid and solid phases (differences less than 3%). Gallium and n-octadecane melt near the room temperature, making them the preferred materials for experimentally investigating melting and solidification of low or high Prandtl-number PCMs.

The previously mentioned works were mostly focused on studying separately melting or solidification and only recently alternate melting and solidification cycles were addressed (Wang et al., 2010b;[START_REF] Rakotondrandisa | Numerical modelling of a melting-solidification cycle of a phase-change material with complete or partial melting[END_REF]. Prior to these studies on complete cycles, periodic melting and solidi-1.3 Physical models for phase-change materials fication problems were the most studied in the literature. [START_REF] Ho | Periodic melting within a square enclosure with an oscillatory surface temperature[END_REF] and [START_REF] Voller | Cyclic phase change with fluid flow[END_REF] studied numerically periodic melting in square enclosures. Recently, [START_REF] Hosseini | Experimental and computational evolution of a shell and tube heat exchanger as a PCM thermal storage system[END_REF] studied melting and solidification of a cylindrical PCM during charging and discharging processes and [START_REF] Chabot | Solid-liquid phase change around a tube with periodic heating and cooling: Scale analysis, numerical simulations and correlations[END_REF] studied analytically the effect of an alternate heating and cooling in a cylindrical PCM, with periodical boundary conditions. We have also contributed recently to the understanding of the governing mechanism during the melting and the solidification process in the paper [START_REF] Rakotondrandisa | Numerical modelling of a melting-solidification cycle of a phase-change material with complete or partial melting[END_REF]. We analysed in detail the difference between solidification occuring after a partial melting or a complete melting by providing temporal evolutions of solid-liquid interface, liquid fraction, Nusselt number and accumulated heat input.

While the literature about melting and solidification of PCM heated from the side is very abundant, research on melting of PCM heated from below is quite rare. [START_REF] Diaz | Visualization of the solid-liquid interface morphology formed by natural convection during melting of a solid from below[END_REF]; [START_REF] Hale | Solid-liquid phase-change heat transfer and interface motion in materials cooled or heated from above or below[END_REF] have studied experimentally the solid-liquid interface morphology of the PCM during basal heating. [START_REF] Gong | Flow and heat transfer in convection-dominated melting in a rectangular cavity heated from below[END_REF] studied numerically the flow and heat transfer during the melting of pure n-octadecane in a rectangular cavity heated from below. Recent numerical simulations have investigated different boundary conditions, such as periodic configurations along the horizontal axis [START_REF] Esfahani | Basal melting driven by turbulent thermal convection[END_REF][START_REF] Madruga | Dynamic of plumes and scaling during the melting of a Phase Change Material heated from below[END_REF][START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF] or wavy surface in a rectangular cavity heated from below [START_REF] Kousksou | Melting over a wavy surface in a rectangular cavity heated from below[END_REF]. Concerning theoretical works, [START_REF] Vasil | Dynamic bifurcations and pattern formation in melting-boundary convection[END_REF]; [START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF] have studied the hydrodynamic instabilities at the onset of convection and compared their observations with the classical Rayleigh-Bénard instability mechanism [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF]. [START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF] have focused on the effect of the non-planar topography of the interface on the convection flow. On the other hand, [START_REF] Gong | Flow and heat transfer in convection-dominated melting in a rectangular cavity heated from below[END_REF]; [START_REF] Esfahani | Basal melting driven by turbulent thermal convection[END_REF]; [START_REF] Madruga | Dynamic of plumes and scaling during the melting of a Phase Change Material heated from below[END_REF] mostly focused on global quantities, such as the heat flux and statistical properties of the interface.

Physical models for phase-change materials

Temperature gradients induce buoyancy forces in the liquid (melted) phase and generate a significant convective flow. The appropriate mathematical description of the liquid phase is thus the usual model for the natural convection flow: the incompressible Navier-Stokes system of equations with Boussinesq approximation for thermal (buoyancy) effects (e. g. [START_REF] Viskanta | Natural convection in melting and solidification[END_REF]). In this model, the energy conservation equation is written as a convection-diffusion equation for the temperature. In the solid phase, conduction is the main phenomenon and the appropriate model is the classical heat equation. The main modelling difficulty is to link these two models by taking into account the separation of the two phases by a sharp interface, across which thermodynamic properties are discontinuous.

We offer in this section a description of the two main approaches suggested in the literature to deal with solid-liquid phase change problem. For a comprehensive review of models for phase-change problems with convection, see [START_REF] Kowalewski | Phase change with convection: modelling and validation[END_REF]. Note that a different category of models was recently suggested in the literature, based on the Lattice Boltzmann Method [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF][START_REF] Gong | Numerical simulation of melting with natural convection based on lattice Boltzmann method and performed with CUDA enabled GPU[END_REF] or meshless methods [START_REF] Atluri | The meshless local Petrov-Galerkin (MLPG) method[END_REF]. Such methods based on non-deterministic models are not discussed in this introduction.

A first modelling approach, usually referred to as the multi-domain (or deforming-grid) method, is based on the classical Stefan two-phase model. Solid and liquid domains are separated and the corresponding conservation equations are solved in each domain. Boundary conditions at the interface between domains are obtained by imposing the Stefan condition (balance of heat fluxes at the interface). The position of the solid-liquid interface is tracked and moved explicitly using either front tracking or front fixing methods. The former method uses deforming grids to reconstruct the interface, while the latter is based on a time-depending coordinate transform, mapping the physical domain into a fixed computational domain. For a detailed description of these methods, see e. g. [START_REF] Sparrow | Analysis of melting in the presence of natural convection in the melt region[END_REF]; [START_REF] Unverdi | A front-tracking method for viscous incompressible, multifluid flows[END_REF]; [START_REF] Gupta | A moving grid numerical scheme for multi-dimensional solidification with transition temperature range[END_REF]; [START_REF] Tenchev | Finite element moving mesh analysis of phase change problems with natural convection[END_REF]. The main drawback of deforming-grid methods is their algorithmic complexity, which makes difficult to accurately capture solid-liquid interfaces of complicated shape or structure (e. g. with mushy regions between solid and liquid phases). Configurations with multiple interacting interfaces are also difficult to simulate with these methods (see also [START_REF] Stella | Modelling Methodologies for Convection-Diffusion Phase-Change Problems[END_REF]).

The second modelling approach avoids to impose explicitly the Stefan condition at the solid-liquid interface and therefore uses a single-domain (or fixed-grid) model. The same system of equations is solved in both liquid and solid phases. The energy balance at the interface is implicitly taken into account by the model. Consequently, the position of the interface is obtained a posteriori by postprocessing the computed temperature field. Phase-field methods [START_REF] Fabbri | The phase-field method in the sharp interface limit: a comparison between model potential[END_REF] and enthalpy methods [START_REF] Voller | A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[END_REF][START_REF] Cao | A numerical analysis of Stefan problems for generalized multidimensional phase-change structures using the enthalpy transforming model[END_REF] are the most commonly used single-domain models. In phase-field methods, a supplementary partial differential equation for the evolution of the order parameter (a continuous variable taking the value 0 in the solid and 1 in the liquid) has to be solved, coupled with the conservation laws [START_REF] Shyy | Computational fluid dynamics with moving boundaries[END_REF]. This new equation is model dependent and its numerical solution could lead to diffuse solid-liquid interfaces. For recent contributions in this area, see [START_REF] Boettinger | Phase-field simulation of solidification[END_REF]; [START_REF] Singer-Loginova | The phase field technique for modeling multiphase materials[END_REF]; [START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF]. We focus below on enthalpy methods, which are the most widely used single-domain models due to their algorithmic simplicity.

The main idea behind enthalpy models is to formulate the energy conservation law in terms of enthalpy and temperature, and include latent heat effects in the definition of the enthalpy. The obtained equation applies to both liquid and solid phases and implicitly takes into account the separation of the phases. Another advantage of enthalpy methods, when compared to previously described models, is to remove the limitation of the phase-change occurring at a fixed temperature. The presence of mushy regions can be easily modelled with these methods. Two types of formulations of enthalpy methods exist in the literature, depending on the main variable used to solve the energy equation: enthalpy or temperaturebased methods.

In enthalpy-based formulations the main variable is the enthalpy [START_REF] Eyres | The calculation of variable heat flow in solids[END_REF][START_REF] Rose | A method for calculating solutions of parabolic equations with a free boundary[END_REF][START_REF] Bhattacharya | An enthalpy method for modeling eutectic solidification[END_REF]. Temperature is computed from the temperature-enthalpy coupling model. An iterative loop is necessary to solve the energy equation, formulated with both enthalpy and temperature variables. For a review of different iterative techniques to solve the energy equation, see [START_REF] Voller | Cyclic phase change with fluid flow[END_REF]. A second variety of enthalpy-based formulations consists in rewriting the energy equation with enthalpy terms only [START_REF] Rady | Natural convection during melting and solidification of pure metals in a cavity[END_REF][START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF].

In temperature-based formulations, the energy equation is formulated in terms of temperature only. The latent heat is treated either by deriving an apparent heat capacity coefficient to define the total enthalpy [START_REF] Morgan | An improved algrorithm for heat conduction problems with phase change[END_REF][START_REF] Chiesa | Natural convective heat transfer rates during the solidification and melting of metals and alloy systems[END_REF][START_REF] Gau | Melting and solidification of a metal system in a rectangular cavity[END_REF] or by introducing a source term in the energy equation [START_REF] Voller | Cyclic phase change with fluid flow[END_REF][START_REF] Swaminathan | Towards a general numerical scheme for solidification systems[END_REF]. Advantages and drawbacks of each approach are discussed in detail in [START_REF] König-Haagen | A comprehensive benchmark of fixedgrid methods for the modeling of melting[END_REF].

Numerical approaches

Single-domain methods are very appealing for numerical implementations. The same system of equations is solved in the entire computational domain, making possible algorithmic or computer-architecture optimisations. If enthalpy models offer an elegant solution to deal with the same energy conservation 6 1.4 Numerical approaches equation in both phases, a last modelling problem has to be solved. It concerns the extension of the Navier-Stokes-Boussinesq equations from the liquid to the solid phase. Different techniques to bring the velocity to zero in the solid region were suggested.

The most straightforward is the switch-off technique, which decouples solid and liquid computational points and overwrites the value of the velocity by setting it to zero in the solid region. Different implementations of this technique with finite-volume methods are presented in [START_REF] Ma | Solid velocity correction schemes for a temperature transforming model for convection phase change[END_REF]; Wang et al. (2010a).

In variable viscosity techniques (Gartling, 1980;Voller et al., 1987;[START_REF] Cao | A numerical analysis of phase change problem including natural convection[END_REF], the fluid viscosity depends on the temperature and is artificially increased to huge values in the solid regions through a regularisation or mushy zone. To avoid blow-up or numerical inconsistencies, the large gradients of viscosity must be correctly resolved in the mushy region. This is naturally achieved in finite-element methods with dynamical mesh adaptivity [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF], while in finite-volume methods with fixed grids, the time step has to be adapted to the space resolution [START_REF] Ma | Solid velocity correction schemes for a temperature transforming model for convection phase change[END_REF]. Versions of the variable viscosity approach suggested in [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] were further studied by [START_REF] Aldbaissy | A full discretisation of the time-dependent Boussinesq (buoyancy) model with nonlinear viscosity[END_REF]; [START_REF] Woodfield | Stability and finite element approximation of phase change models for natural convection in porous media[END_REF] and implemented in a different finite-element framework (FEniCS) by [START_REF] Zimmerman | Monolithic Simulation of Convection-Coupled Phase-Change: Verification and Reproducibility[END_REF].

A third technique used to ensure a zero velocity field in the solid phase is the so-called enthalpyporosity model [START_REF] Brent | Enthalpy-porosity technique for modeling convectiondiffusion phase change: Application to the melting of a pure metal[END_REF]. A penalisation source term is introduced in the momentum equation to allow the switch from the full Navier-Stokes equations in the liquid phase to a Darcy equation for porous media. The mushy region is thus regarded as a very dense porous medium that sharply brings the velocity to zero in the solid region. The expression of the penalisation source term generally follows the Carman-Kozeny model for the permeability of a porous medium [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF][START_REF] Hannoun | A reference solution for phase change with convection[END_REF][START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF], but other mathematically equivalent expressions were suggested [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF]. Different formulations and implementations of the enthalpy-porosity model are presented in [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]; Giangi et al. (2000); [START_REF] Stella | Modelling Methodologies for Convection-Diffusion Phase-Change Problems[END_REF].

Concerning the space discretization of these models, finite difference (FD) or finite volume (FV) methods are generally used in the literature. When single-mesh models are used, the general strategy to capture the solid-liquid interface is to dramatically increase the mesh resolution in the whole domain. This results in a considerable increase of the computational time, even for two dimensional cases. [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF] have reported that the simulation of the melting of tin with a 200 × 200 fixed grid required 2, 400 CPU hours, 111 runs (restarts), and 3 months of calculation. Dynamical mesh adaptivity becomes in this context a valuable tool to concentrate the grid refinement effort only in regions displaying high gradients of the computed variables (melting-solidification fronts, thermal or viscous boundary layers, recirculation zones).

Finite element (FE) methods offer the possibility to dynamically refine the mesh only in specific regions of the domain. Different FE approaches were suggested, from enthalpy-type methods (e. g. [START_REF] Elliott | Error analysis of the enthalpy method for the Stefan problem[END_REF]) to front-tracking methods (e. g. [START_REF] Li | A finite-element front-tracking enthalpy method for Stefan Problems[END_REF]). Recently, adaptive FE methods were proposed for classical two-phase Stefan problem (Belhamadia et al., 2004a) using an anisotropic mesh adaptation algorithm based on solution-dependent metrics. The authors extended their algorithm for the threedimensional simulation of the same problem (Belhamadia et al., 2004b) and showed that the use of locally adapted meshes with strong anisotropy proved to be very effective in reducing the number of computational nodes for such phase-change systems without convection.

To simulate melting or solidification problems with convection, [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] recently suggested a dynamical mesh adaptation algorithm based on metrics control and implemented with the FreeFem++ software [START_REF] Hecht | FreeFem++ (manual)[END_REF][START_REF] Hecht | New developments in Freefem++[END_REF]. The advantage of this adaptive finite-element method, which will be also used in the present study, is to make possible, with reasonable computational cost, the re-meshing of the computational domain at each time step. A very refined discretization of the regularization zone between solid and liquid phases is thus obtained, while regions with low gradients are coarsened in order to balance the overall computational effort.

Present numerical approach

The present study is based on a single-domain enthalpy-porosity model for solid-liquid phase change problems with convection. For the energy conservation equation, a temperature-based formulation takes into account the latent heat by introducing a discontinuous source term. For the mass and momentum conservation equations, we solve in the entire domain the incompressible Navier-Stokes equations with Boussinesq approximation for buoyancy effects. To bring the velocity to zero in the solid phase, we introduce in the momentum equation a penalty term following the Carman-Kozeny model. The coupled system of momentum and energy equations is integrated in time using a second-order Gear scheme. All the terms are treated implicitly and the resulting discretized equations are solved using a Newton method [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF].

The advantage of this formulation is to permit a straightforward implementation of different types of non-linearities. For the space discretization we use Taylor-Hood triangular finite elements, i. e. P 2 for the velocity and P 1 for the pressure. Temperature is discretized using P 2 or P 1 finite elements. Discontinuous variables (latent heat, thermal diffusivity, etc) at the solid-liquid interface are regularized through an intermediate artificial mushy region.

Single domain methods require a refined mesh near the interface, where large enthalpy gradients have to be correctly resolved. An optimized dynamical mesh adaptivity algorithm allows us to adapt the mesh every time step and thus accurately capture the evolution of the interface. Mesh adaptivity, a feature of the current method, offers the possibility to deal with complicated phase-change cases, involving multiple solid-liquid interfaces.

There are three main novelties in the present numerical approach, when compared to [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]:

(i) we use the Carman-Kozeny model to bring the velocity to zero inside the solid phase, instead of a viscosity penalty method (imposing a large value of the viscosity in the solid),

(ii) we increase the time accuracy of the algorithm by replacing the first-order Euler scheme with the second-order Gear (BDF2) scheme (see also [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF]),

(iii) we improve the metric calculation procedure for the mesh adaptivity.

The programs were built and organized as a toolbox for FreeFem++ [START_REF] Hecht | New developments in Freefem++[END_REF][START_REF] Hecht | FreeFem++ (manual)[END_REF], which is a free software (under LGPL license). FreeFem++ 1 offers a large variety of triangular finite elements (linear and quadratic Lagrangian elements, discontinuous P 1 , Raviart-Thomas elements, etc.) to solve partial differential equations. It is an integrated product with its own high level programming language and a syntax close to mathematical formulations, making the implementation of numerical algorithms very easy. Among the features making FreeFem++ an easy-to-use and highly adaptive software we recall the advanced automatic mesh generator, mesh adaptation, problem description by its variational formulation, automatic interpolation of data, colour display on line, postscript printouts, etc. The FreeFem++ programming framework offers the advantage to hide all technical issues related to the implementation of the finite element method. It becomes then easy to use the present toolbox to 1.6 Purpose of the thesis code new numerical algorithms for similar problems with phase-change.

Purpose of the thesis

The purpose of the present work was to investigate numerically solid-liquid phase-change systems. The investigation tool used in this thesis is the open-source software FreeFem++ [START_REF] Hecht | FreeFem++ (manual)[END_REF][START_REF] Hecht | New developments in Freefem++[END_REF]. A high accuracy numerical model using a Newton method with adaptive finite elements is used to simulate phase-change problems with natural convection.

A first numerical investigation of convective phase-change problems using an adaptive finite element method has been carried out by [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]. Their method solved the Navier-Stokes-Boussinesq equations by the mean of a single domain approach using first order scheme in time. The technique of variable viscosity approach was applied to bring to zero the velocity in the solid phase. The study was focused on two-dimensional square cavity configurations.

A first objective of this thesis was to improve the existing code, developed by the Numerical methods and applications group of the LMRS Laboratory2 , and to organize the program as a toolbox for the software FreeFem++. To this end, the following objectives were addressed:

(i) increase the accuracy by using a second order scheme for the time discretization and P 2 finite elements for the temperature, (ii) implement a Carman-Kozeny model, in addition to the viscosity-based approach, to bring to zero the velocity in the solid region, (iii) investigate challenging cases by simulating complex geometries (highly distorted meshes, cylindrical PCM with inner heated tubes) and computationally demanding cases (high Rayleigh numbers),

(iv) simulate a complete melting-solidification cycle of a PCM.

A second objective was to extend the program to three-dimensional configurations. The two dimensional assumption is indeed not valid for high Rayleigh problems, especially for basal melting cases. Moreover, three dimensional adaptive finite element methods for convective melting problems are less present in the literature.

A third objective was to provide a thorough analysis of both melting and solidification processes and compare with numerical results.

Thesis plan

Chapter 2 sets the mathematical and physical basis of the numerical system used to simulate phasechange problems involving natural convection. We present in detail the incompressible Navier-Stokes-Boussinesq equation and the single-domain approach to solve the same system of equations in the whole domain. The enthalpy method, modeling the phase-change phenomenon, is presented first. The Navier-Stokes equation with Boussinesq approximation to simulate the natural convection in the liquid flow is then recalled. Finally, the non-dimensional system of equations is described in detail, with a discussion on the Carman-Kozeny penalty term in the momentum equation.

Chapter 3 is devoted to the numerical algorithm for solving the numerical system presented previously. The finite element algorithm that we have developed in this work is first presented: the time integration scheme, the finite element discretization and the Newton method. Then, we describe the mesh adaptivity by metric control, which is a standard function offered by FreeFem++. Some theoretical tests to assess for the accuracy of our numerical method are also presented. The space and time convergence orders are demonstrated using the Burggraf flow and a manufactured solution designed for the incompressible Navier-Stokes equation. The structure of the new finite-element toolbox for the simulation of PCMs is also described. The program architecture and the parameter details are presented. Finally, the domain decomposition method used for large simulations is described.

A first validation of our numerical method is presented in Chapter 4. The capability of the code to deal with linear and non-linear forms of the buoyancy force in the Boussinesq approximation is first tested. Test cases are presented by increasing progressively the difficulty. We simulate first the natural convection of air in square enclosures, differentially heated from the vertical walls. Then, a heated obstacle is added in the center of the domain. Finally, we add the non-linearity in the Boussinesq force to simulate the natural convection of water.

Once the capability of our algorithm to deal with natural convection problems was demonstrated, much more attentions to phase-change problems is paid in Chapter 5. Numerical simulations of melting and solidification of PCM are presented. The melting processes for octadecane and Gallium inside rectangular enclosures are investigated first since they were extensively used for numerical method validations. Their physical properties are approximately equal in both solid and liquid phases. Finally, several PCM container geometries are simulated to prove the robustness of our numerical method. In particular, cylindrical and irregular domains are computed.

A scale correlation analysis is presented in Chapter 6. We compare the behavior of the PCM for lateral and basal heating. The melting of n-octadecane heated from the side is carried out first. Then, numerical results for the melting from below are presented. For both cases, we provide a comprehensive description of the melting process through a scale analysis for a better understanding of the heat transfer mechanism. The time evolution of some physical parameters, such as the Nusselt number, the liquid fraction, the accumulated heat input and the time evolution of the melting are discussed and analysed.

Chapter 7 presents the numerical simulation of a full melting/solidification cycle of a PCM. Two solidification fronts have to be tracked, which makes the case very challenging. A differentially heated cavity case and a circular PCM with inner cooled tube are studied. The influence of the Rayleigh number during the melting and the solidification processes is emphasized, since both cycles are not driven by the same mechanism.

Chapter 8 is devoted to three-dimensional configurations simulated using parallel algorithms. 3D simulations of PCMs are indeed less investigated in the literature.

Finally, Chapter 9 draws the conclusion of this study and some ideas for future developments.

Navier-Stokes-Boussinesq and enthalpy-porosity models

We consider a solid-liquid system placed in a 2D or 3D domain Ω of characteristic length L ref .

In the following, subscripts s and l will refer to the solid and liquid phases, respectively.

The single domain approach, using the same system of equations in both phases, is described first in detail. The model is based on the Navier-Stokes equations with Boussinesq approximation, which is the natural description of the fluid flow with natural convection. A penalty term is added to the momentum equations to bring the velocity to zero inside the solid region. For the energy conservation equation, an enthalpy-porosity method is used to model the phase change process.

Enthalpy-porosity model

The phase change process is modeled using an enthalpy method with temperature-based formulation (Voller et al., 1987;[START_REF] Cao | A numerical analysis of Stefan problems for generalized multidimensional phase-change structures using the enthalpy transforming model[END_REF][START_REF] Cao | A numerical analysis of phase change problem including natural convection[END_REF]. We start from the classical energy equation:

∂(ρh) ∂t ϕ + ∇ • (ρhU ) -∇ • (k∇T ) = 0, (2.1)
where t ϕ is the physical time, h the enthalpy, ρ the density, U the velocity vector, T the temperature and k the thermal conductivity. The total enthalpy h is transformed as the sum of the sensible heat and the latent heat:

h = h sen + h lat = c(T + s(T )), (2.2)
with c the local specific heat. The function s is introduced to model the jump of the enthalpy during the solid-liquid transition. For pure materials, s is theoretically a Heaviside step function depending on the temperature: it takes the zero value in the solid region and a large value in the liquid, equal to h sl /c l , with h sl the latent heat of fusion.

If the phase-change is assumed to occur within a mushy zone defined by a small temperature interval

T ∈ [T f -T ε , T f + T ε ]
around the temperature of fusion T f , a model for s(T ) is necessary. Linear (Voller et al., 1987;Wang et al., 2010a) or more smooth functions [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] can be used to regularize s(T ) and thus model the jump of material properties from solid to liquid. In the current work we use a regularization of all step-functions (latent heat source, specific heat, thermal diffusivity or conductivity) by a continuous and differentiable hyperbolic-tangent function suggested by [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF].

We assume moreover that the undercooling phenomenon is negligible during the solidification stage (see also Wang et al. (2010b); [START_REF] Kowalewski | Phase change with convection: modelling and validation[END_REF]). The solid-liquid interface is identified by the isoline T = T f , while the Gibbs-Thomson effect due to the surface energy of the solid-liquid interface is assumed to be negligible.

Equation (2.1) can be further simplified by considering the following assumptions:

(i) The density difference between solid and liquid phases is negligible, i. e. ρ l = ρ s = ρ is constant. We note however that this is not strictly true for all substances, but it serves here as a convenient simplification.

(ii) The regularization zone is narrow and the velocity inside this zone is negligible.

Consequently, the final form of the energy equation is obtained by combining Eqs. (2.2) and (2.1) and neglecting the convection term ∇ • (csU )1 :

∂ (cT ) ∂t ϕ + ∇ • (cT U ) -∇ • k ρ ∇T + ∂ (cs) ∂t ϕ = 0. (2.3)

Navier-Stokes equations with Boussinesq approximation

The natural convection in the liquid part of the system is modeled using the incompressible Navier-Stokes equations, with Boussinesq approximation for buoyancy effects. To make this model valid for both liquid and solid phases, the momentum equation is modified as following:

∂U ∂t ϕ + (U • ∇)U + 1 ρ ∇P - µ l ρ ∇ 2 U + ρge y = A(T )U , (2.4)
where P denotes the pressure, µ l the dynamic viscosity of the liquid (assumed to be constant), and g the gravitational acceleration.

The penalty term A(T )U is artificially introduced in Eq. (2.4) to extend this equation in the solid phase, where the velocity, pressure, viscosity and Boussinesq force are meaningless. Consequently, A(T ) is modelled to vanish in the liquid, where the Navier-Stokes-Boussinesq momentum equation is recovered. A large value of A(T ) is imposed in the solid, reducing the momentum Eq. (2.4) to A(T )U = 0, equivalent to U = 0. The exact expression for A will be given in the next section.

The density is assumed to be constant everywhere except for the body force term ρge y in Eq. (2.4).

Under the assumption of a small variations of density and temperature, the Boussinesq approximation allows to linearize the density as follows:

ρ = ρ ref (1 -β(T -T ref )), (2.5) with β = -(1/ρ ref )(∂ρ/∂T ) ρ ref the thermal expansion coefficient and (ρ ref , T ref ) a reference state.
It is worth noting that this approximation is valid for β(T -T ref ) considerably smaller than the unity. Therefore, the momentum Eq. (2.4) can be written as Finally, the conservation of mass in the liquid phase is expressed by the continuity equation in the frame of incompressible fluids:

∂U ∂t ϕ + (U • ∇)U + ∇p -ν l ∇ 2 U -f B (T )e y = A(T )U , ( 2
∇ • U = 0.
(2.7)

The final system of equations for the single-domain approach is thus:

∇ • U = 0, (2.8) ∂U ∂t ϕ + (U • ∇)U + ∇p -ν l ∇ 2 U -f B (T )e y -A(T )U = 0, (2.9) ∂ (cT ) ∂t ϕ + ∇ • (cT U ) -∇ • k ρ ∇T + ∂ (cs) ∂t ϕ = 0.
(2.10)

Dimensionless system of equations for the single-domain approach

It is convenient to numerically solve a dimensionless form of the previous equations. After choosing a reference length L ref (usually the height of the cavity when a rectangular domain is considered) and a reference state (ρ, V ref , T ref ), we can define the following scaling for the space, velocity, temperature and time variables:

x = X L ref , u = U V ref , θ = T -T ref δT , t = V ref L ref t ϕ .
(2.11)

T ref is the reference temperature and in most cases T ref = T f (the temperature of fusion), unless otherwise specified. Consequently, the non-dimensional temperature of fusion is set to θ f = 0. The temperature difference δT , defines a temperature scale, that will be set differently for melting and solidification cases. δT is considered as the representative temperature scale for the natural convection onset in the liquid region. For the classical natural convection problem without phase-change, δT is generally defined as δT = T h -T c since the flow in the fluid is driven by the temperature difference between the "hot" and the "cold" temperature. However, for the melting PCM, the convection is driven by the temperature difference δT = T h -T f , with T f the temperature of fusion. As far as the solidification process is concerned, a distinct discussion will be provided in Chapter 7.

The dimensionless system of equations to be solved in both liquid and solid regions can be finally written as:

∇ • u = 0, (2.12) ∂u ∂t + (u • ∇)u + ∇p - 1 Re ∇ 2 u -f B (θ) e y -A(θ)u = 0, (2.13) ∂ (Cθ) ∂t + ∇ • (Cθu) -∇ • K RePr ∇θ + ∂ (CS) ∂t = 0. (2.14)
f B is the linearised Boussinesq buoyancy force, Re is the Reynolds number relating the inertia forces to the viscous forces, and Pr the Prandtl number expressing the ratio of the molecular momentum and thermal diffusivities:

f B (θ) = Ra PrRe 2 θ, Re = ρV ref L ref µ l = V ref L ref ν l , Pr = ν l α l , (2.15)
with α = k/(ρc) the thermal diffusivity. In the expression of f B , the Rayleigh number of the flow, which indicates the relative importance of the buoyancy term with respect to viscous effects, is defined as

Ra = gβL 3 ref δT ν l α l .
(2.16)

If previous non-dimensional numbers are pertinent only in the liquid phase, the non-dimensional conductivity and specific heat are defined in both phases

K(θ) = k k l = 1, θ ≥ θ f , k s /k l , θ < θ f . , C(θ) = c c l = 1, θ ≥ θ f , c s /c l , θ < θ f .
(2.17)

The non-dimensional function S = s/s l in the energy Eq. (2.14) takes a similar non-dimensional form:

S(θ) = s s l =    h sl /c l δT = 1 Ste , θ ≥ θ f , 0, θ < θ f , (2.18)
with Ste the Stefan number.

Discontinuous step-functions defined in Eqs. (2.17) and (2.18) are replaced by continuous and differentiable hyperbolic-tangent functions, generically defined for all θ by the formula [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF])

F (θ; a r , θ r , R r ) = f l + f s -f l 2 1 + tanh a r θ r -θ R r , (2.19)
where f l , f s are the imposed values in the liquid and solid phases, a r a smoothing parameter, θ r the central value (around which we regularize), and R r the smoothing radius. For example, we use for the non-dimensional source term in Eq. (2.14) the following regularisation over the artificial mushy region θ ∈ [-ε, ε]:

S(θ) = 1 Ste - 1 2 Ste 1 + tanh θ s -θ R s , (2.20)
with θ s = θ f = 0 and R s = ε for the melting case.

Finally, the penalty term in the momentum Eq. (2.13) is derived from the Darcy's law, by modeling the fluid flow within the mushy region as a flow through a porous medium. In fact, the Darcy's law states that the velocity of flow in porous medium is proportional to the pressure gradient:

u = - ζ * µ ∇p. (2.21)
where ζ * is the permeability, which is a function of the porosity. As the porosity decreases, the permeability (and the velocity) also decreases, down to the limiting value of zero when the mushy zone becomes completely solid. This behavior can be accounted in a numerical model by adding a source term Au in the momentum equation. The well-known equation derived from the Darcy law is the Carman-Kozeny model:

∇p = - C CK (1 -λ) 2 λ 3 u.
(2.22) 14

Boundary layer approximation and scale analysis

Consequently, A takes the form [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF]Kheirabadi and Groulx, 2015)

A(θ) = - C CK (1 -L f (θ)) 2 L f (θ) 3 + b , (2.23)
where L f (θ) is the local liquid fraction, which is 1 in the fluid region and 0 in the solid. L f is regularized inside the artificial mushy-region using a hyperbolic-tangent similar to Eq. (2.20). The Carman-Kozeny constant C CK is set to a large value (as discussed below) and the constant b = 10 -6 is introduced to avoid divisions by zero.

Boundary layer approximation and scale analysis

Regardless to the practical use of PCM (energy storage, building insulation, or for other purposes), one would necessarily assess the heat transfer during the phase-change process. It was extensively proven that the convective heat transfer plays a significant role during the melting stage. Therefore, before solving numerically Eqs. (2.12) -(2.14), we first rely on a scale analysis to predict theoretically the fluid flow and heat transfer patterns that can develop in the fluid part. The idea behind the scaling analysis is to identify the proper scales of the phenomenon, in order to understand the evolution of the heat transfer and the melting rate. We consider only the liquid phase without phase-change, inside a two-dimensional enclosure of height H filled with Newtonian fluid, differentially heated from the vertical walls and insulated at the horizontal walls. A homogeneous Dirichlet boundary condition is prescribed for the velocity. One may refer to the book by [START_REF] Bejan | Convection heat transfer[END_REF] for a more detailed presentation.

Since no external force is applied to our system, the fluid flow is mainly driven by natural convection, induced by the temperature difference from the vertical walls. The heat transfer through the fluid layer immediately adjacent to the wall is assumed to be driven by pure conduction. We therefore define the average Nusselt number to quantify the heat transfer rate at the hot wall (placed at x = 0):

Nu = - 1 0 ∂θ ∂x x=0
dy.

(2.24)

When a steady state could be reached, the fluid near each sidewall is characterized by two boundary layers: a thermal boundary layer of thickness δ θ and a viscous boundary layer of thickness δ ν . The boundary layer approximation assumes that the flow and the energy transfer are restricted predominantly to the boundary layer region. The main consequences of the boundary layer approximation are that:

(i) the normal part of the momentum has a negligible importance, (ii) the downstream diffusion terms in the momentum and energy equations are neglected in comparison with the normal diffusion terms, (iii) the pressure distribution is purely hydrostatic, (iv) the thermal and the viscous boundary layer thickness are related by following expression:

δ ν /δ θ = o Pr 1/2 .
The scaling laws suggested by [START_REF] Bejan | Convection heat transfer[END_REF] are presented in Fig. 2.1, and rely on boundary layer approximation:

δ T ∼ HRa -1/4 if Pr ≫ 1, HPr -1/4 Ra -1/4 if Pr ≪ 1.
(2.25)

3 A finite element method for the Navier-Stokes-Boussinesq equations using a Newton algorithm This chapter sets the numerical algorithm for solving the system of equations described previously in Chapter 2 for both two and three-dimensional configurations. The finite element method and the Newton algorithm are presented first in detail. Then, the domain decomposition method for the threedimensional configuration is discussed.

Motivation for the choice of the numerical method

In many in-house or commercial codes used to simulate phase-change problems (ANSYS CFX, Fluent, etc.), finite difference (FD) or finite volume (FV) methods are used most of the time with a fixed-mesh approach. Accurate simulations are therefore carried out by refining considerably the mesh resolution in the whole domain, and thus increasing dramatically the computational time. In the example of [START_REF] Voller | Cyclic phase change with fluid flow[END_REF] who used insufficient grid resolution to compute the melting of Gallium, Wang et al. (2010b) had to consider thinner meshes, resulting in a strong increase, of a factor of three, of the total mesh nodes to capture correctly the boundary layer structure.

Our choice for the finite element (FE) method is motivated by its capability to adapt dynamically the mesh and to deal with several geometrical domains. The adaptive capabilities of the FE discretization to use finer meshes where sharp phenomena take place (solid-liquid interfaces, boundary layers, recirculation zones) and coarser meshes elsewhere (in the solid, in the bulk of the fluid region where the gradient are lower than in the boundary layer region) is helpful to reduce the degrees of freedom involved in the numerical resolution and thus reduce the computational time. We use a finite-element method that was implemented using the open-source software FreeFem++ [START_REF] Hecht | FreeFem++ (manual)[END_REF][START_REF] Hecht | New developments in Freefem++[END_REF].

In this study, the use of mesh adaptivity proved mandatory to obtain accurate results within reasonable computational time. The numerical code was optimized to afford the mesh refinement every time step: the mesh density was increased around the phase change interfaces, offering an optimal resolution of the large gradients of all regularized functions (S, K, L f ), while the mesh was de-refined (larger triangles) in the solid part, where a coarser mesh could be used. A simulation using a globally refined mesh would require a prohibitive computational time for an equivalent accuracy of the melting front resolution. Similar algorithms based on FreeFem++ were successfully used for solving different systems of equations with locally sharp variation of the solution, such as Gross-Pitaevskii equation [START_REF] Danaila | A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose-Einstein condensates[END_REF][START_REF] Vergez | A finite-element toolbox for the stationary Gross-Pitaevskii equation with rotation[END_REF] or Laplace equations with nonlinear source terms [START_REF] Zhang | Existence and numerical modelling of vortex rings with elliptic boundaries[END_REF].

The space discretization is based on Taylor-Hood finite elements, approximating the velocity with P 2 Lagrange finite elements (piecewise quadratic), and the the pressure with the P 1 finite elements (piecewise linear). The temperature and the enthalpy are discretized using P 2 finite elements. This discretization is second order accurate in space. We also use a second order accurate discretization in time. A fully implicit backward second order scheme (BDF2 or GEAR) is employed in the present study. The time derivative of a variable φ is approximated by:

dφ dt ≃ 3φ n+1 -4φ n + φ n-1 2δt , (3.1)
computing the solution φ n+1 at time t n+1 = (n + 1)δt by using two previous states (φ n , φ n-1 ). We use this scheme to advance in time both velocity (φ = u) and temperature fields (φ = θ). The other terms in Eqs. (2.12) -(2.14) are treated implicitly (i. e. taken at time t n+1 ). The resulting non-linear equations are solved using a Newton algorithm.

Some authors used a Richardson extrapolation [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF]Wang et al., 2010b) in the momentum equation by extrapolating U from previous time steps. We have tested this approach, but the results exhibited less accurate solutions and the computations requested small δt. In addition, a projection algorithm with an explicit discretization of the Navier-Stokes-Boussinesq equations was also investigated, using the Adams-Bashforth and Crank-Nicolson second order schemes. The main drawback is the very small δt (∼ 10 -6 ) needed to ensure convergence at each time steps. A last alternative we tested for the treatment of non-linear terms in the momentum equation is the characteristics Galerkin method [START_REF] Pironneau | Characteristic-Galerkin and Galerkin/least-squares space-time formulations for the advection-diffusion equations with time-dependent domains[END_REF], but this method is not convenient for high Ra simulations and requires very small time steps.

Finally, a supplementary difficulty comes from the one-domain method which needs techniques to bring the velocity to zero in the solid region. The switch-off technique, the variable viscosity approach and the Carman-Kozeny penalty term are the most used in the literature. We note that the physical meaning of the variable viscosity formulation and the Carman-Kozeny penalty term is fundamentally different. The Carman-Kozeny approach considers the mushy zone as a porous media, i. e. the solid is stationary and the liquid flows through the porous structure, while the variable viscosity formulation treats the mushy zone as a mixture of solid crystals and liquid, permitting thus a movement of both the solid and the liquid. The viscosity-based method was investigated by [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] and the Carman-Kozeny penalty method is investigated in the present work. We note however that techniques in FV methods based on the modification of the numerical algorithm to switch-off the velocity in the solid by a relaxation scheme also exist.

Finite element algorithm

To solve the system of Eqs. (2.12) -(2.14) we use a finite-element method. Finite-element methods for solving Navier-Stokes type systems of equations are generally based on a separate discretization of the temporal derivative (using finite difference, splitting or characteristics methods) and the generalization of the Stokes problem for the resulting system [START_REF] Temam | Navier-Stokes equations and nonlinear functional analysis[END_REF][START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]. We use the second-order implicit finite-difference discretization (3.1) of the temporal derivative and obtain the time semi-discretization of the single-domain model (2.12) -(2.14): 18

Finite element algorithm

∇ • u n+1 + γp n+1 = 0, (3.2) 3 2 u n+1 δt + (u n+1 • ∇)u n+1 + ∇p n+1 - 1 Re ∇ 2 u n+1 (3.3) -A(θ n+1 )u n+1 -f B (θ n+1 ) e y = 2 u n δt - u n-1 2δt , 3 2 θ n+1 + S(θ n+1 ) δt + ∇ • u n+1 θ n+1 -∇ • K RePr ∇θ n+1 = (3.4) 2 θ n + S(θ n ) δt - θ n-1 + S(θ n-1 ) 2δt .
The penalty parameter γ is introduced to fill the diagonal of the pressure term at the algebraic level, when a Dirichlet boundary condition is applied for the velocity. The average value of the pressure field and the divergence-free velocity field, computed for different values of the penalty parameter γ, for the case of the natural convection of air, are given in Tab. 3.1. γ 0.1 10 -3 10 -7 10 -10 One can observe a decreasing values of Ω ∇ • u for very low values of γ. However, larger values of Ω p are observed for decreasing values of the penalty parameter. To ensure a good approximation of the pressure field with zero average, a value of γ = 10 -7 is set for all simulations. This system of non-linear equations is solved at time t n+1 = (n + 1)δt, using two previous states at t n and t n-1 .

Ω ∇ • u 1.
The space discretization of variables over the domain Ω uses a finite-element method based on a weak formulation of the system of Eqs. (3.2) -(3.4). We consider homogeneous Dirichlet boundary conditions for the velocity, i. e. u = 0 on ∂Ω, and set the classical Hilbert spaces for the velocity and pressure:

V = V × V, V = H 1 0 (Ω), Q = q ∈ L 2 (Ω) Ω q = 0 . (3.5)
Following the generalization of the Stokes problem [START_REF] Temam | Navier-Stokes equations and nonlinear functional analysis[END_REF][START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF], the variational formulation of the system of Eqs. (3.2) -(3.4) can be written as:

find (u n+1 , p n+1 , θ n+1 ) ∈ V × Q × V , such that b u n+1 , q -γ(p n+1 , q) = 0, ∀ q ∈ Q (3.6) 3 2δt u n+1 , v + c u n+1 ; u n+1 , v + 1 Re a u n+1 , v -(A(θ n+1 ) u n+1 , v) + b v, p n+1 -f B (θ n+1 ) e y , v (3.7) = 2 δt (u n , v) - 1 2δt u n-1 , v , ∀ v ∈ V 3 2δt θ n+1 + S(θ n+1 ), φ + u n+1 • ∇θ n+1 , φ + K(θ n+1 ) ReP r ∇θ n+1 , ∇φ (3.8) = 2 δt (θ n + S(θ n ), φ) - 1 2δt θ n-1 + S(θ n-1 ), φ , ∀ φ ∈ V,
where (u, v) = Ω u • v denotes the scalar product in L 2 (Ω) or L 2 (Ω) 2 ; the bilinear forms a, b and trilinear form c are defined as [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]:

a : V × V → R, a(u, v) = Ω ∇ t u : ∇v = 2 i,j=1 Ω ∂ j u j • ∂ j v i , b : V × Q → R, b(u, q) = - Ω ∇ • u q = - 2 i=1 Ω ∂ i u i • q, c : V × V × V → R, c(w; z, v) = Ω [(w • ∇) z] • v = 2 i,j=1 Ω w j (∂ j z i )v i .
The system of non-linear Eqs. (3.6) -(3.8), regarded as a general equation F(w) = 0, with w = (u n+1 , p n+1 , θ n+1 ), is solved using a Newton method. To advance the solution from time t n to t n+1 , we start from an initial guess w 0 = (u n , p n , θ n ) (which is the solution at t n ), and construct the Newton sequence w k = (u k , p k , θ k ) by solving for each inner iteration k:

D w F(w k )w k+1 = D w F(w k )w k -F(w k ).
(3.9) D w F is the linear operator representing the differential of F. Equation (3.9) can be rewritten as follows:

b (u k+1 , q)γ(p k+1 , q) = 0, (3.10)

3 2δt (u k+1 , v) + c (u k+1 ; u k , v) + c (u k ; u k+1 , v) (3.11) + 1 Re a (u k+1 , v) - dA dθ (θ k ) θ k+1 u k , v -(A(θ k ) u k+1 , v) + b (v, p k+1 ) - df B dθ (θ k ) θ k+1 e y , v = 1 δt 2u n - 1 2 u n-1 , v +c (u k ; u k , v) - dA dθ (θ k ) θ k u k , v - df B dθ (θ k ) θ k -f B (θ k ) e y , v ,
3 2δt θ k+1 + dS dθ (θ k ) θ k+1 , φ + (u k • ∇θ k+1 , φ) + (u k+1 • ∇θ k , φ) (3.12) + K(θ k ) RePr ∇θ k+1 , ∇φ + dK dθ (θ k ) θ k+1 ReP r ∇θ k , ∇φ = 2 δt (θ n + S(θ n ), φ) + (u k • ∇θ k , φ) + 3 2δt dS dθ (θ k ) θ k -S(θ k ), φ - 1 2δt θ n-1 + S(θ n-1 ), φ + dK dθ (θ k ) θ k ReP r ∇θ k , ∇φ .
Note that the last term of Eq. (3.11) cancels in the case of a linear Boussinesq force f B (see Eq. (2.15)); this is not the case when non-linear variations of the density of the liquid are considered (convection or solidification of water). Note also that the previous system of Eqs. (3.10) -(3.12) depends only on u n , u n-1 , θ n and θ n-1 and is independent of p n , the pressure being in this approach a Lagrange multiplier for the divergence free constraint. We underline the fact that the Newton loop (following k) has to be 20

Mesh adaptivity

iterated until convergence for each time step δt following the algorithm Navier-Stokes time loop following n set w 0 = (u n , p n , θ n ) Newton iterations following k solve Eqs. (3.10) -(3.12) to get w k+1 stop when w k+1w k < ξ N actualize (u n+1 , p n+1 , θ n+1 ) = w k+1 .

(3.13)

In some cases, it is more relevant to solve directly the steady equation corresponding to Eqs. (3.6) -(3.8) (i. e. without the temporal derivatives), when a steady state could be achieved. It is, for example, the case for the natural convection of air or water in differentially heated configurations and low Ra numbers. In the stationary case, the weak formulation of the problem becomes: b u n+1 , qγ(p n+1 , q) = 0, (3.14)

u n+1 , v + c u n+1 ; u n+1 , v + a u n+1 , v (3.15) -(A(θ n+1 ) u n+1 , v) + b v, p n+1 -α f B (θ n+1 ) e y , v = 0, u n+1 • ∇θ n+1 , φ + K RePr ∇θ n+1 , ∇φ = 0. (3.16)
The parameter α is introduced in front of the buoyancy force term f B as a continuation parameter of the value of the Ra number. The idea is to solve successively steady problems by increasing α from 0 to 1. At each stage, the previously computed solution is used as initial condition. These steps could be assimilated to a smooth increase of the Ra number: (3.17) with N the number of steps of the continuation algorithm. This approach allows to compute the steady natural convection of air in 6 iterations for Ra = 10 6 instead of 150 iterations when the unsteady system of equations is solved until a steady state is reached.

Loops for m = 1 : N steps set w 0 = (u n , p n , θ n ) Newton iterations following k solve (3.14) -(3.16) with α = m N 4 stop when w k+1 -w k < ξ N actualize (u n+1 , p n+1 , θ n+1 ) = w k+1 ,

Mesh adaptivity

We use the standard mesh adaptivity function (adaptmesh) offered by FreeFem++ [START_REF] Hecht | New developments in Freefem++[END_REF].

The key idea implemented in this function (see also [START_REF] Borouchaki | Anisotropic adaptive mesh generation in two dimensions for CFD[END_REF]; Castro- [START_REF] Castro-Diaz | Anisotropic Grid Adaptation for Inviscid and Viscous Flows Simulations[END_REF]; [START_REF] Hecht | Mesh Adaptation by Metric Control for Multi-scale Phenomena and Turbulence[END_REF]; [START_REF] George | Delaunay triangulation and meshing[END_REF]; [START_REF] Frey | Maillages. Hermès[END_REF]; [START_REF] Mohammadi | Applied Shape Design for Fluids[END_REF]) is to modify the scalar product used in the automatic mesh generator to evaluate distance and volume. Equilateral elements are thus constructed, accordingly to the new metric. The scalar product is based on the evaluation of the Hessian H of the variables of the problem. For example, for a P 1 discretization of a variable χ, the interpolation error is bounded by:

E = |χ -Π h χ| 0 ≤ c sup T ∈T h sup x,y,z∈T |H(x)|(y -z, y -z), (3.18)
where Π h χ is the P 1 interpolate of χ, |H(x)| is the Hessian of χ at point x after being made positive definite. Using the Delaunay algorithm (e. g. [START_REF] George | Delaunay triangulation and meshing[END_REF]) to generate a triangular mesh with edges close to the unit length in the metric M = |H| (cE) will result in a equally distributed interpolation error E over the edges a i of the mesh. More precisely, we get

1 cE a T i Ma i ≤ 1. (3.19)
The previous approach could be generalized for a vector variable χ = [χ 1 , χ 2 ]. After computing the metrics M 1 and M 2 for each variable, we define a metric intersection M = M 1 ∩ M 2 , such that the unit ball of M is included in the intersection of the two unit balls of metrics M 2 and M 1 . For details, see the procedure defined in [START_REF] Frey | Maillages. Hermès[END_REF].

The adaptmesh function offers the possibility to take into account several metrics computed from different variables monitoring the evolution of the phase-change systems. For natural convection systems, the mesh will be adapted using the values of the two velocity components and the temperature. For phase-change systems, to accurately track the solid-liquid interface we add the variation of the enthalphy source term in the adaptivity criterion. For water systems (convection or freezing), we also add an extra function tracking the anomalous change of density around 4 o C. To reduce the impact of the interpolation on the global accuracy for time-depending problems, we consider, for each variable used for adaptivity, the metrics computed at actual (t n+1 ) and previous (t n ) time instants (see also Belhamadia et al. (2004a)). The anisotropy of the mesh is a parameter of the algorithm and it was set to values close to 1. This is an inevitable limitation since we also impose the minimum edge-length of triangles to avoid too large meshes. The capabilities of the mesh adaptivity algorithm are illustrated in Chapter 5.

A finite-element toolbox for the simulation of phase-change systems with natural convection

The methods described previously were implemented in a 2D toolbox based on FreeFem++ . The syntax to implement the Newton algorithm is very close to the mathematical formulation given above.

After defining a vectorial finite-element space fespace Wh(Th,[P2,P2,P1,P1]);, associated to the mesh Th, we define the velocity, pressure and temperature variables in a compact manner by Wh [u1,u2,p,T];. Corresponding test functions are defined similarly. It is then very easy to define a problem formulation in FreeFem++ and include all the terms of the algorithm (3.10) -(3.12). This makes the reading of the programs very intuitive by comparing each term to its mathematical expression. New terms could be added to the variational formulation expressed in the problem structure, without affecting other parts of the program. Consequently, the implementation of new models or numerical methods for this problem is greatly facilitated by this modular structure of programs.

In this section we first describe the architecture of the programs and the organisation of files. Then we focus on the list of input parameters and the structure of output files. Test cases and common macros are separated into two folders.

Program architecture

This directory is organized as follows:

1. The directory Common-Macros contains five files:

• Macro operator.idp includes macros and functions defining mathematical operators,

• Macro problem.idp: macros defining the variational formulation of the problem,

• Macro restart.idp: macros used to start a new simulation from a saved field,

• Macro output.idp: macros used to save the solution with different formats,

• Macro system.idp: macros identifying the OS and defining specific OS-commands.

2. The directory Test-Cases contains four subdirectories, each of them defining one of the following applications:

• natural convection of air or water in a differentially heated square cavity,

• melting of a PCM stored in containers of different shapes,

• melting followed by solidification of a rectangular PCM,

• freezing of pure water in a square cavity.

Each subdirectory contains three files: NEWTON $case.edp is the main FreeFem++ script file, param phys.inc defines the physical parameters and param num.inc the numerical parameters. For example, to run the natural convection case of air in a square cavity, one can use the following command in a terminal window:

FreeFem ++ NEWTON stat natconv . edp

The folder structure of each test case is illustrated in Fig.

Input parameters

Physical parameters and parameters related to the run are separated into two files.

(1) The file param phys.inc contains the physical descriptions of the problem:

• typeT: is the finite-element type for the temperature, with possible values P2 or P1,

• Torder: is the accuracy order of the time integration scheme, with possible values 1 (Euler scheme) or 2 (Gear scheme),

• scalAdim: defines the characteristic scales of the problem, see Eq. (2.11). Possible values 1, 2 or 3 correspond to the following choice of the characteristic scales [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]:

(1) :

V (1) ref = ν l H =⇒ t (1) ref = H 2 ν l =⇒ Re = 1, (3.20) (2) : V (2) ref = α H =⇒ t (2) ref = t (1) ref Pr =⇒ Re = 1/Pr, (3.21) (3) : V (3) ref = ν l H Ra Pr =⇒ t (3) ref = t (1) ref Pr Ra =⇒ Re = Ra Pr , (3.22)
• x l , x r , y l , y r : are the values defining the dimensions of the cavity [x l , x r ] × [y l , y r ],

• Pr, Ra, Ste: are the Prandtl, Rayleigh and Stefan numbers, see Eq. (2.16) and (2.15),

• T hot , T cold : are dimensionless temperatures according to Eq. (2.11),

• bcu 1 , bcu 2 , bcT: are macros defining the velocity (u) and the temperature (T ) boundary conditions.

• epsi: is the half width ε of the mushy region. Default value = 0.01,

• dt: is the dimensionless time step,

• t max : is the dimensionless final time,

• Parameters for regularization functions:

The parameters of the hyperbolic-tangent function in Eq. ( 2.19), used to regularize discontinuous functions are set by default as follows:

f s f l a s θ s R s C CK b Enthalpy 0 1/Ste 1 0.01 0.01 - -
Carman -Kozeny 0 1 1 0.01 0.01 10 6 10 -7

Conductivity (water) 1 2.26/0.578 1 θ f 0.015 --

• rho(T) and Drho(T): (water cases only) define the density and its derivative as functions of the temperature, following the model [START_REF] Gebhart | A new density relation for pure and saline water[END_REF]:

ρ(T ) = ρ m (1 -ω|T -T m | q ), ρ m [kg/m 3 ] ω [ o C -q ] q T m [ o C]
999.972 9.2793 • 10 -6 1.894816 4.0293

• f B (T), df B (T): define the buoyancy force and its derivative.

(2) The file param num.inc contains the parameters controlling the run.

Restart parameters:

• Nsave: the solution is saved every Nsave time steps in the Data folder (see Fig. 3.2). The temperature and the velocity fields are saved in Tecplot and Medit folders, while the liquid fraction, the Nusselt number, and the accumulated heat input are saved in the Gnuplot folder.

• Nrestart: restart files (mesh and solution) are saved every Nrestart time steps. Solutions at current and previous iterations, the CPU time, the accumulated heat input Q 0 , and the time step dt are saved in the folder RST.

• Ncondt: allows the user to stop the run and save the solution properly. The file OUTPUT/zz.condt is read every Ncondt time steps: if the user replaces the value "0" in this file by "1" the run is stopped. This is a simple solution for a clean stop of the job by the user. Default value = 20.

• Nremesh: the mesh is adapted every Nremesh iterations. If this parameter is set to "1" the mesh is adapted every time step.

• IFrestart: is a boolean controlling the set up of the initial field.

IF restart = 0, the initial condition is built in the code for each test case. For the PCM melting cases, the PCM is initially motionless at isothermal temperature. To set-up a smooth initial field, a few time steps (with very small δt) are computed by increasing progressively the boundary temperature at the hot wall and the Rayleigh number (by continuation). Outputs are saved in OUTPUT/Data-RST-0.

IF restart > 0, (positive integer values) the solution field previously computed at iteration IF restart is loaded from the folder OUTPUT/Data-RST-filenameRST/RST, with filenameRST a variable selecting the restart folder.

IF restart < 0, (negative integer values), the same principle for loading a solution is used, but from the folder INIT (see Fig. 3.2). The solution fields stored in this folder could come from different previous calculations (e. g. a steady state solution or, for the water, the natural convection field before freezing).

Newton parameters:

• epsconv: is the value of the stopping criterion for steady cases,

• gamma: is the penalty parameter in Eq. (3.2). Default value = 10 -7 ,

• tolNewton: is the Newton tolerance ξ N (see Eq. (3.13)). Default value = 10 -6 ,

• newtonMax: limits the maximum number of iterations in the Newton algorithm (3.13).

Default value = 50, Mesh parameters:

• nbseg: is the number of segments for the discretisation along the x and y directions,

• errh: is the interpolation error level. Default value = 0.02,

• hmin, hmax: are the minimum and maximum edge size, respectively,

• adaptratio: is the ratio for a prescribed smoothing of the metric. For a value less than 1.1 no smoothing is done. Default value = 1.5,

• nbvx: is the maximum number of vertices allowed in the mesh generator. Default value = 50000.

Output parameters:

• dircase: is the name of the output folder,

• fcase: is the prefix-name for ouput files.

• Tecplot, Medit, Gnu: correspond to the name of the visualisation software to be used; the format of the outputs written in OUTPUT/Data (see Fig. 3.2) is accordingly set. The files from the Tecplot folder can be easily read also with Paraview.

Outputs

When a computation starts, the OUTPUT directory is created (see Fig. 3.2). It contains two folders storing the output data and the echo of the run parameters. The folder Data contains four subdirectories with different output format files of the computed solution. File names are created using the prefix defined by the parameter fcase, the current iteration and the current dimensionless time t. Solution files can be visualized using either Tecplot or any other CFD Visualization tools (Paraview, Visit, etc.). Moreover, .gmsh (mesh) and .rst (fields) files are generated in the folder RST to enable restarts of the computation. Note that the folder FFglut contains FreeFem++ scripts that re-read and visualize the RST-files to facilitate the selection of a restart field. An .echo file with a summary of the main parameters, informations on the run and the names of the output files is saved in the folder RUNPARAM. This directory additionally contains a copy of the .inc parameter files, allowing an easy identification of each case and preparing an eventual rerun of the same case.

Numerical tests of the accuracy of the numerical method

We start by presenting tests of the accuracy of our numerical method. We used the technique of manufactured solutions (e. g. Roache ( 1998)) which has the advantage of providing an exact solution to a modified problem, related to the initial one. The general idea is to modify the original system of equations by introducing an extra source term, such that the new system admits an exact solution given by a convenient analytic expression. Even though in most cases exact solutions constructed in this way are not physically realistic, this approach allows one to rigorously verify computations. We tested the space and time accuracy using manufactured solutions for the system of Eqs. (2.12) -(2.14) for a stationary case (Burggraf flow) and a time-dependent one [START_REF] Nourgaliev | Fullyimplicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change[END_REF]. For both cases, we computed the global error ε for different norms in space:

ε = Φ e -φ h , (3.23)
with Φ e the exact solution and φ h the numerical solution. Computations were performed for the convection of air (C = K = 1, A(θ) = S(θ) = 0), with a Rayleigh number Ra = 10 4 and a Prandtl number Pr = 0.71.

Space accuracy: Burggraf stationary flow with thermal effects

The Burggraf manufactured solution is a time-independent recirculating flow inside a square cavity

[0, 1] × [0, 1].
It is similar to the well-known entrained cavity flow, with the difference that the velocity singularity at the top corners of the cavity is avoided. We added to the classical Burggraf flow [START_REF] Shih | Effect of grid staggering on numerical schemes[END_REF][START_REF] Laizet | High-order compact schemes for incompressible flows: A simple and efficient method with quasi-spectral accuracy[END_REF] a manufactured solution for the temperature, with constant temperature imposed at the top and the bottom walls. Vertical walls are assumed to be adiabatic. The exact solution of the new flow with thermal effects is with σ > 0 a scaling parameter and functions

u 1 (x, y) = σg ′ (x)h ′ (y), (3.24) u 2 (x, y) = -σg ′′ (x)h(y), p(x, y) = σ Re h (3) (y)g(x) + g ′′ (x)h ′ (y) + σ 2 2 g ′ (x) 2 h(y)h ′′ (y) -h ′ (y) 2 , T (x, y) = T c + (T h -T c )y + a(x)b(y),
g(x) = x 5 5 - x 4 2 + x 3 3 , (3.25) h(y) = y 4 -y 2 , a(x) = cos(πx), b(y) = y(1 -y).
Note that the velocity at the top border of the cavity is

u 1 (0, 1) = 2σ(x 4 -2x 3 + x 2 ), u 2 (x, 1) = 0, (3.26)
which ensures the continuity of the velocity at corners (u(0, 1) = u(1, 1) = 0), since no-slip walls are imposed for the other borders:

u(x, 0) = u(0, y) = u(1, y) = 0.
The forcing terms that have to be added to the momentum and energy (temperature) equation are derived by injecting the exact solution (3.24) into the system of Eqs. (2.12) -(2.14):

f u 1 = 0, (3.27) f u 2 = σ 2 h(y)h ′ (y) g ′′ (x) 2 -g ′ (x)g (3) (x) + σ Re g (4) (x)h(y) + 2g ′′ (x)h ′′ (y) + g(x)h (4) (y) + σ 2 2 g ′ (x) 2 h(y)h (3) (y) -h ′ (y)h ′′ (y) - Ra PrRe 2 T (x, y), f T = u 1 (x, y)a ′ (x)b(y) + u 2 (x, y) T h -T c + a(x)b ′ (y) - K RePr a ′′ (x)b(y) + a(x)b ′′ (y) .
We used the Taylor-Hood finite element (P 2 for the velocity and P 1 for the pressure) and tested P 1 or P 2 finite elements for the temperature. Figures 3.3a and 3.3b illustrate the streamlines and the temperature field, respectively.

Figure 3.4 plots the discretization error ε as a function of the grid size h = δx = δy for the temperature and the velocity. Both L 2 and L ∞ norms are displayed. The expected second order accuracy in L 2norm is obtained with P 1 finite elements (Figs. 3.4a and 3.4b), while an order exceeding three is observed when P 2 finite elements are used (Fig. 3.4c and 3.4d).

The optimal convergence for the velocity, the pressure and the temperature is shown in Fig. 3.5, for their natural norms: H 1 norms for the velocity and the temperature, and L 2 norm for the pressure.

A seconder order convergence is observed for the velocity and the pressure, whatever P 1 or P 2 finite elements are considered for the temperature (Figs. 3.5a -3.5c), consistent with the Taylor-Hood finite elements employed. Moreover, as expected, a first order accuracy is obtained for the temperature using P 1 finite elements (Fig. 3.5b) and a second order using P 2 finite elements (Fig. 3.5d).

Time accuracy: manufactured unsteady solution

To test the time accuracy of the Gear (BDF2) scheme, we used the manufactured time-dependent solution suggested in [START_REF] Nourgaliev | Fullyimplicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change[END_REF]:

u 1 (x, y, t) = (δU 0 + α u sin(t)) cos(x + γ 1 t) sin(y + γ 2 t), (3.28) u 2 (x, y, t) = -(δU 0 + α u sin(t)) sin(x + γ 1 t) cos(y + γ 2 t), T (x, y, t) = T + (δT 0 + α t sin(t)) cos(x + γ 1 t) sin(y + γ 2 t), p(x, y, t) = P + (δP 0 + α p sin(t)) sin(x + γ 1 t) cos(y + γ 2 t),
The values of the constants are reported in Table 3.2.

γ 1 γ 2 P T δP 0 δT 0 δU 0 α p α u α t 0.1 0.1 0 1.0 0.1 1.0 1.0 0.05 0.4 0.1 Table 3.2: Parameter for the time-dependent manufactured solution (3.28).

The corresponding forcing source terms are where a = (x + γ 1 t), b = (y + γ 2 t) and U c = (δU 0 + α u sin(t)), T c = (δT 0 + α u sin(t)), P c = (δP 0 + α u sin(t)).

f u 1 = α u cos(t) cos(a) sin(b) -U c γ 1 sin(a) sin(b) + U c γ 2 cos(a) cos(b) (3.29) -U c u 1 (x, y, t) sin(a) sin(b) + U c u 2 (x, y, t) cos(a) cos(b) + P c cos(a) cos(b) + 2 Re u 1 (x, y, t), f u 2 = -α u cos(t) sin(a) cos(b) -U c γ 1 cos(a) cos(b) + U c γ 2 sin(a) sin(b) -U c u 1 (x, y, t) cos(a) cos(b) + U c u 2 (x, y, t) sin(a) sin(b) -P c sin(a) sin(b) + 2 Re u 2 (x, y, t) - Ra PrRe 2 T (x, y, t), f T = α t cos(t) cos(a) sin(b) -T c γ 1 sin(a) sin(b) + T c γ 2 cos(a) cos(b) -T c u 1 (x, y, t) sin(a) sin(b) + T c u 2 (x,
Guided by the results obtained in Sec. 3.5.1 for the space accuracy, we fixed the grid size to h = dx = 0.01 and we applied P 2 discretization for the temperature, to ensure small spatial discretization errors. For diminishing values of the time step δt, the solution was evolved in time up to the time instant t max = π at which the error (3.23) was computed. The time convergence is displayed in Fig. 3.6 for the temperature (panel a) and the velocity (panel b) variables, for which the expected second order convergence in time is obtained.

Domain decomposition method with FreeFem++: FFDDM

Solving the Navier-Stokes-Boussinesq systems of Eqs. (3.6) -(3.8) in three-dimensional configurations could generate large-size problems. The natural convection of air in a cube of dimensions [0, 1] 3 with 40 × 40 × 40 uniform grid involves 3 millions of unknowns (d.o.f.) in the linear system, when a P 1 finite element is considered for the temperature. For such a large-size problem, memory limitation can rapidly arise with sequential algorithms. It is thus essential to distribute data among several processors.

A natural approach is the domain decomposition method (DDM). DDM aims at dividing the computational domain in many subdomains on which we solve local problems with adequate interface conditions. Two families of DDM exists: non-overlapping methods and overlapping methods such as the Schwarz method.

We use in this work a domain decomposition Schwarz method, enhanced by coarse space corrections 32 3.6 Domain decomposition method with FreeFem++: FFDDM through the FreeFem++ library ffddm. This enrichment with a coarse space is mandatory in the present study to avoid the lack of robustness related to the one-level method. The main drawback of one-level methods is their convergence rate that depends on the number of subdomains, resulting in a poor scaling for large problems. This is mostly due to a lack of global communication between subdomains that exchange informations with their direct neighbors only. The additional coarse space has consequently the task to spread the information to all subdomains at each iteration.

The data distribution among processors is done via an overlapping domain decomposition and involves a related linear algebra. To enable parallel computing, the mesh is first split into subdomains using Scotch or Metis libraries. Figure 3.7 illustrates the domain decomposition of a cube of dimension [0, 1] 3 into 8 subdomains with Metis graph partitioner. Mesh adaptivity using metrics control makes possible the optimisation of the distribution of mesh elements. The final linear system of Eqs. (3.10) -(3.12) is solved in parallel using a GMRES Krylov method, with an Optimized Restricted Additive Schwarz (ORAS) preconditioner. To solve the linear equation Ax = b, the ORAS preconditionner reads

M -1 RAS = N j=1 R T j D j (R j AR T j ) -1 R j , (3.30)
R j denote the restriction operators and D j are square diagonal matrices. Local matrices are defined as:

A j = R j AR T j . (3.31)
The duplicated unknowns due to the overlap between subdomains are coupled via a partition of unity

I = N i=1 R T i D i R i . (3.32)
Thus, the global solution U is defined as:

U = N i=1 R T i D i R i U = N i=1 R T i D i U i . (3.33)
4 Numerical simulation of natural convection flow in a 2D square cavity

We first focus on the capability of our code to deal with natural convection flows in square enclosures.

The convection flow is a relevant case to validate rigorously the Navier-Stokes-Boussinesq solver. A large number of benchmarks on natural convection problems induced by temperature difference could be found in the literature, with connection to engineering (air circulation in buildings, double-wall insulations, solar collectors, etc.) and geophysical systems.

In this chapter, we are interested on the natural convection of fluids in a square cavity differentially heated from the vertical walls. The fluid temperature rises and its density decreases along the heated wall, convecting the fluid up to the point where it reaches the cold wall, along which the reverse process occurs. These two simultaneous opposed effects create a recirculation cell within a stationary zone in the center.

We solve the system of Eqs. (2.12) -(2.14), with A(θ) = 0 in the momentum equation and S(θ) = 0 in the energy equation. Linear and non-linear expressions of the buoyancy force f B (T ) are investigated, by simulating the natural convection of air and water, respectively. Natural convection of water exhibits a non-linear variation of the density with a maximum value around T = 4 o C. A linear variation with temperature is generally assumed for the natural convection of air in the Boussinesq approximation. We consider a square enclosure of height H. Physical properties of air and water used in our simulations are listed in Tab. 4.1. Isothermal boundary conditions are applied at the vertical walls and adiabatic boundary conditions to the upper and lower walls. Quantitative and qualitative validations are carried out as following. (i) The velocity profile along symmetry lines, the maximum value of u at mid-domain (x = 0.5 and y = 0.5), locations x max and y max of this maximum and the Nusselt number are compared with the spectral-accurate simulations by Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF] and [START_REF] De Vahl | Natural convection of air in a square cavity: a benchmark numerical solution[END_REF] in Sec. 4.1.1 for natural convection of air within a square cavity. (ii) A heated obstacle is included in the center of the domain and the transversal velocity profile along the horizontal symmetry lines is compared with numerical results of [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF] in Sec. 4.1.2. (iii) The temperature profile along the horizontal symmetry line is compared with numerical results of [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF] in Sec. 4.2 for the natural convection of water.

ρ µ c p k α β kg/m 3 kg/(m s) J/(kg K) W/(m K) m 2 /s 1/

Natural convection of air

We start by testing the Newton algorithm (3.10) with linear expression of f B (θ) (see Eq. (2.15)).

The classical problem of the thermally driven square cavity with adiabatic top and bottom walls is simulated. We consider a cavity of height H = 0.1 m, initially filled with motionless air and a linear distribution of the temperature. The dimensionless parameters describing the investigated configuration are based on the fluid properties presented in Tab. 4.1. Pr = 0.71 and Rayleigh number ranges from Ra = 10 4 to 10 8 . The characteristic scales of the problem are:

L ref = H, T ref = T h + T c 2 , (4.1)
and,

V ref = ν l H Ra Pr =⇒ t ref = H 2 ν l Pr Ra =⇒ Re = Ra Pr . (4.2)
The left wall is heated with dimensionless hot temperature θ h = 0.5 and the right wall is cooled with θ c = -0.5 (resulting from Eq. (4.1)). A homogenous Dirichlet boundary condition (u = 0) is applied for the velocity. It has been shown by Le Quéré and [START_REF] Quéré | From onset of unsteadiness to chaos in a differentially heated square cavity[END_REF] that the solution of the 2-D Boussinesq equation in this configuration becomes unsteady at critical Rayleigh number Ra = 1.82 • 10 8 due to the instability mode that breaks the usual central-symmetry of the solution. Therefore, a steady state can be achieved for the chosen Rayleigh numbers.

We solve the unsteady and steady Navier-Stokes-Boussinesq equations, related to algorithms (3.13) and (3.17), respectively. The unsteady case is computed until the steady state with a single convection cell is reached, with a numerical tolerance of 10 -9 . In parallel, the steady case is performed using a Rayleigh number continuation: a smaller value of the Rayleigh number is set initially and increased smoothly until reaching the correct value. At each stage, the computation starts from the solution obtained from the previous Rayleigh number simulation. Two cases are simulated: i) a differentially heated square cavity and ii) a differentially heated cavity with an inner heated obstacle. For each of them, the horizontal and the vertical velocity profiles u(y) and v(x) at mid-domain (y = 0.5 and x = 0.5, respectively) are plotted and compared with numerical results by Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF] and [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF]. Moreover, for i), we compare with the benchmark solutions by [START_REF] De Vahl | Natural convection of air in a square cavity: a benchmark numerical solution[END_REF] and Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF]: the maximum value u max and v max with their location y max and x max , the Nusselt number at the heated (left) wall and at the mid-domain, and the maximum and the minimum value of the local Nusselt number at x = 0 with their y-locations. [START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF], who used a fourth-order finite difference method for the spatial discretization and a third-order backward finite difference scheme for the time integration, for Ra ranging from 10 4 to 10 6 . The temperature distribution and the streamlines at the steady state display a good agreement with the benchmark solution. The higher the Ra number, the stronger is the natural convection flow in the cavity. A single convection cell is observed in the center of the cavity for Ra = 10 4 (Fig. 4.1d) during which the heat transfer is merely bulk heat transfer. Figures 4.1e and 4.1f however exhibit stronger convection with more rolls, when the boundary layer heat transfer dominates.

Differentially heated square cavity

A more accurate validation is performed with respect to the spectral-accurate results of Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF]. We plot in Fig. 4.2 the horizontal (panel a) and the vertical (panel b) profiles of the velocity and compare the current simulation with data extracted from Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF] for each of the three computed Ra numbers. Results from Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF] are represented by solid lines and the current simulation by symbols. A very good agreement can be noticed for each case. The thickness of the boundary layer is decreasing with increasing value of the temperature difference δT . Following the correlation of [START_REF] Bejan | Convection heat transfer[END_REF], a viscous boundary layer with a dimensionless thickness of order of δ ν ∼ 0.02 should be present close to the vertical walls for Ra = 10 6 . Accordingly, the mesh resolution should allow to capture these structures. We assess on the influence of the mesh resolution on the solution in Tab. 4.2. We report the L 2 and L ∞ norms of simulation is compared with the accurate solutions, V LQ and U LQ , of Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF]. Differences decrease when the mesh resolution is increased. We observe that errors of order of 10 -5 are obtained for the 80 × 80 grid resolution, with a reasonable computational time of 14 seconds. All computations in this section are performed with a fixed triangular mesh, generated by the Delaunay algorithm starting with nbseg = 80 points on each side of the square. Higher Ra number simulations were performed to test the robustness of our algorithm. We depict the temperature isolines at the steady state for Ra = 10 7 , 4 • 10 7 , and 10 8 in Fig. 4.3. One can observe the onset of detached regions along the upper and lower adiabatic walls for Ra = 4 • 10 7 , which is in agreement with observations of Le Quéré and De [START_REF] Quéré | Computation of natural convection in two-dimensional cavities with Chebyshev polynomials[END_REF]. According to Le Quéré and De [START_REF] Quéré | Computation of natural convection in two-dimensional cavities with Chebyshev polynomials[END_REF], this phenomenon is in part induced by the Neumann boundary conditions at the horizontal (top and bottom) walls. To assess for the steadiness of the solution, we plot in Fig. 4.4 the time evolution of the Nusselt number Nu 0 at the heated wall (x = 0) and Nu 1/2 at the mid-domain (x = 0.5). The Nusselt number in the center region of the cavity is defined as

ε V = |V -V LQ | (along x-direction) and nbseg ||ε U || 2 ||ε U || ∞ ||ε V || 2 ||ε V || ∞ 40 1.777477 • 10 -4 5.
Nu 1/2 = √ PrRa 1 0 uθ - ∂θ ∂x x=0.5 dy. (4.3)
The Nusselt at the hot wall becomes constant very soon when compared to Nu 1/2 , which oscillates during a longer period and tends to the constant value of Nu 0 = 16.5237. Table 4.3 offers a more quantitative assessment of the accuracy of the present Newton method. The values of (u max , y max ), (v max , x max ) are compared to reference values from [START_REF] De Vahl | Natural convection of air in a square cavity: a benchmark numerical solution[END_REF] and Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF] for Ra ranging from 10 4 to 10 8 . The Newton method gives results identical to reference values. Relative errors lower than 0.01% are obtained for Ra = 10 5 and 10 6 . A relative difference of 0.43% is observed for Ra = 10 6 with respect to the reference solution of [START_REF] De Vahl | Natural convection of air in a square cavity: a benchmark numerical solution[END_REF], and a relative difference of 0.003% with the spectral accurate solution of Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF]. The highest discrepancy is noted for Ra = 10 8 , but still offers reasonable agreement with benchmark values, within 0.52% error for u max 5.57% for v max . We compare moreover the average Nusselt number at the left vertical wall (the mean, the max, and the min) and at the mid-domain, since the Nusselt number will be largely used to compare the optimized configuration of a PCM, either for melting or solidification cycles in the next chapters. It is thus essential to ensure the accuracy of the computed value of this parameter. Excellent agreement with [START_REF] De Vahl | Natural convection of air in a square cavity: a benchmark numerical solution[END_REF] and Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF] is obtained for all simulated Ra numbers, with a relative error lower than 0.038%. [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF] for Ra = 10 4 to 10 8 .

Differentially heated cavity with an inner heated square

Thermally driven cavity including a heated square obstacle is computed in this section. We consider the same configuration presented in Sec. 4.1.1 and a square domain with isothermal boundary conditions is added in the initial set-up. This kind of basic configuration could be representative of telecommunication outdoor cabinet applications, in which the use of passive cooling solutions have begun to be more and more investigated. Indeed, inside an outdoor cabinet, electronic equipments generate heat when active. The study of the flow structures within the enclosure with obstacles was studied experimentally and numerically by [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF]. The simplified model of cavity with rectangular heated obstacles investigated by [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF] will be reproduced in this section to test the robustness of our numerical algorithm.

A linear distribution of the temperature is imposed initially in the motionless air inside the cavity. The obstacle is maintained at a dimensionless hot temperature θ h = 0.8, with a no-slip boundary condition for the velocity. The solutions for Ra = 10 4 , 10 5 , 10 6 and Pr = 0.71 are compared with the result obtained by [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF], who used an immersed boundary method with a FD code using high order schemes for time and spatial discretization.

The temperature distribution in the cavity when the steady-state is reached, for each of the three computed Ra numbers, are reported in panels (a) to (c) of Fig. 4.5. The temperature gradient gives rise to a clockwise circulation and when Ra is increased, vertical thermal boundary layers form distinctly along the differentially heated sidewalls and the obstacle. Consequently, higher is the Rayleigh number, the hot temperature in the center of the domain is more rapidly advected by the natural convection flow into the cold part of the cavity. It is worth noting that for Ra = 10 6 in panels (c) and (f), a stagnant fluid with a stratified temperature forms in a small portion of the fluid between the cold wall and the obstacle.

We present a more accurate validation in panels (d) to (f) of Fig. 4.5. The transversal velocity profiles along the x-axis are plotted and compared with the numerical data of [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF] for each cases. A good agreement can be observed, with a relatively small difference between the extremum of the velocity, while the shape of the velocity profile match well.

Natural convection of water inside a square cavity

We simulate in this section the natural convection of water in a differentially heated cavity. A further difficulty is introduced compared to the previous validations by taking into account the non-linear variation of the density in the buoyancy force. Pure water exhibits a non-linear density variation for T < 10.2 • C with a maximum at T m = 4.0293 • C. We use below the following density-temperature relationship suggested in [START_REF] Gebhart | A new density relation for pure and saline water[END_REF]:

ρ(T ) = ρ m (1 -w |T -T m | q ) , (4.4) with ρ m = 999.972 [kg/m 3 ], w = 9.2793 • 10 -6 [( • C) -q ]
, and q = 1.894816. The bouyancy term

f B = g(ρ ref -ρ)/ρ ref in Eq.
(2.6) becomes after scaling:

f B (θ) = Ra Pr Re 2 1 βδT ρ(θ f ) -ρ(θ) ρ(θ f ) , (4.5)
where β = (1/ρ m ) (dρ/dT ) is the thermal expansion coefficient with the value β = 6.91 [START_REF] Scanlon | A numerical analysis of buoyancy-driven melting and freezing[END_REF].

• 10 -5 [(K) -1 ] (Scanlon
We consider a cavity of height H = 0.38 m filled with liquid pure distilled water. This problem was investigated experimentally and numerically by Giangi et al. (2000); [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]; [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF]. The height H of the cavity is considered as a the length scale of the problem, L ref = H. We choose T ref = T h -T c = 10K in order to compare our simulation with the numerical results of [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF], and define the following scaling:

V ref = ν l H =⇒ t ref = H 2 ν l =⇒ Re = 1. (4.6)
The non-dimensional parameters describing the problem result from the physical properties of water given in Tab. 4.1: Ra = 2.518084 • 10 6 and Pr = 6.99.

The initial temperature is linearly distributed with a hot (non-dimensional) temperature θ h = 1 at the left wall and a cold temperature θ c = 0 at the right wall. The top and the bottom of the cavity are adiabatic and no-slip boundary condition u = 0 is applied for the velocity.

The temperature field of the steady state is presented in Fig. 4.6a. Unlike the natural convection of air, in which two distinct boundary layers along the vertical walls and a stagnant and thermally stratified fluid in the core of the fluid flow were observed, an anomalous variation of the temperature is pointed out around the iso-line θ = θ m = 0.40293 for the natural convection of water. This anomalous thermal variation of water density, is clearly discernible in the streamlines of the steady flow in Fig. 4.6b. Two recirculating zones are formed in the flow: a lower (abnormal) recirculation in the vicinity of the cold wall where θ < θ m and an upper (normal) one where the density decreases with temperature (θ > θ m ).

Following the above observations, a higher mesh resolution should be applied around θ m . We thus define a P 1 function Φ(θ) by the following hyperbolic-tangent function similar to Eq. (2.20): layer regions. To reduce the impact of the interpolation on the global accuracy, since our algorithm is optimized to afford the mesh refinement every time step, we use both Φ(θ n ) and Φ(θ n+1 ) in the adaptivity procedure.

Φ(θ) = 1 2 1 + tanh θ m -θ R Φ , ( 4
The final mesh is displayed in Fig. 4.6c. The mesh is clearly refined along the line θ = θ m , where the structure and the extent of the two recirculating zones should be captured, and along the vertical walls where the heat transfer is dominated by the boundary layer transfer. Furthermore, as expected in the relatively stagnant fluid region, a coarser mesh is applied.

A more accurate comparison is given in Fig. 4.6d. The temperature profile θ(x) along the horizontal symmetry line of the cavity is in good agreement with the numerical results of [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF] obtained with FV and FD codes (FLUENT and FRECONV3V), commonly used in the heat transfer community. Differences are visible in the vicinity of the maximum density line, region where our mesh is well refined to capture the separation line between the two recirculation zones. It should be noted that the FLUENT simulations in [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF] are performed with a fixed uniform grid with 380 × 380 nodes, while our adapted grid has only 3422 triangles.

Concluding remarks

We can conclude from sections. 4.1 and 4.2 that our Newton method is able to deal efficiently with the two-dimensional Navier-Stokes-Boussinesq problem, either with linear or non-linear buoyancy forces.

Our simulations of natural convection of air including or not obstacles are in good agreement with the numerical solutions of Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF] and [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF]. Excellent agreement with [START_REF] De Vahl | Natural convection of air in a square cavity: a benchmark numerical solution[END_REF] and Le [START_REF] Queré | Accuratte solutions to the square thermally driven cavity at high Rayleigh number[END_REF] has also been observed for the value of the average Nusselt number at the heated wall.

The challenging case of the natural convection of water was used to demonstrate the robustness of the 2D code. Good agreement with the numerical simulations of [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF] was noticed. The advantage of the mesh adaptivity was clearly shown by using a total number of grid points 40 times lower than in the simulation of [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF]. The refined mesh along the line θ m have permitted to solve more accurately the structures and the extent of the two recirculating zones.

The properties of the run presented in this chapter are summarized in Tab. 4.4 for the highest value of the Ra number, Ra = 10 6 . The interest of solving the steady equation, when steady state could be reached, is clearly emphasised when regarding the number of iterations and the CPU time. 1, 164 iterations and 3162.86 CPU seconds are necessary to achieve the steady state with a tolerance of 10 -6 for the natural convection of water by solving the unsteady equation, while only 20 iterations and 51.1662 CPU seconds are needed with the steady algorithm. We simulate in the next chapter the melting and the solidification of a PCM. Aside from the nonlinear convection terms in the momentum and the energy Eqs. (2.12) -(2.14), a further difficulty will arise from the non-linearity introduced by the source term ∂(CS)/∂t in Eq. (2.14).

Numerical simulations of phase change materials in 2D configurations

The purpose of this chapter is to validate our numerical method by simulating convective phase-change problems. A discussion on numerical parameters is given first, followed by a sequence of validations computing different physical and geometrical configurations.

The Newton algorithm (3.10) is now solved throughout the whole domain containing liquid and solid phases. When compared to the previous validation, two new non-linearities are now added: the Carman-Kozeny penalty term in the momentum equation and the enthalpy source term in the energy equation. The Carman-Kozeny penalty term is used to ensure zero velocity in the solid region and the enthalpy non-linear source term is included to model the phase-change in the energy equation. Following the same idea as in the natural convection validation cases in Chapter 4, the non-linearity in the body force is gradually added. We consider first the linear form of f B (θ) by investigating the melting of various PCMs within several shapes of the containers. Then, the challenging case of water freezing, characterized by a non-linear variation of the density (see Eq. (4.4) in Chapter 4), is simulated.

Details about the parameter settings are first given in Sec. 5.1: the characteristics of the mushyzone, the mesh adaptivity and the initial condition are discussed in detail. Second, the melting of n-octadecane and Gallium are presented in Secs. 5.2 and 5.3. The physical properties of n-octadecane and Gallium used in our simulations are reported in Tab. 5.1.

ρ µ h SL c p k T f β kg/m 3 kg/(m s) kJ/kg J/(kg K) W/(m K) K 1/K Octadecane 774
3.9 Then, in Secs. 5.4 and 5.5 we demonstrate the capability of our code to deal with different shapes of the domain. The melting of a cylindrical PCM with inner heated tubes is presented in Sec. 5.4 and the solid crust formation in a highly distorted mesh is simulated in Sec. 5.5. Finally, the water freezing case is performed in Sec. 5.6. Besides the non-linear definition of the body force, the striking feature of the water freezing simulation is the tracking of several interfaces, namely the solidification front and the anomalous thermal variation of density. Our result is qualitatively compared with the experimental results of [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF].

Parameter settings

We present in this section the setting of the main numerical parameters. The influence of the penalty term A(θ) in the momentum equation is discussed first, namely the value of the Carman-Kozeny constant C CK . Second, details about the mesh adaptivity parameters are given. Finally, the initialization process is described.

The Carman-Kozeny model

When a single domain method is applied, the Carman-Kozeny penalty term is used to bring to zero the velocity in the solid region, through an intermediate mushy zone. In the mushy zone, at the leading term, the momentum Eq. (2.13) becomes ∂u ∂t = Au.

(5.1)

with solution u = u 0 exp(At).

(5.2)

Since A < 0, see Eq. (2.23), we infer that lim A→∞ u = 0. It is thus expected that, when C CK is large enough, the value of the velocity in the solid region will tend to zero. The effect of this constant on the melting and the solidification processes has attracted some attention in the literature.

The influence of C CK in the location of the solid-liquid front was investigated by Kheirabadi and Groulx (2015). They have concluded that high C CK values induce a slower melting rates, and conversely small values result in unphysical predictions of the melting front evolution. They also noticed that C CK and the melting temperature range ∆T ε = T ε2 -T ε1 are dependent of one another: different values of ∆T ε require different values of C CK to obtain the same melting front location. [START_REF] Mathura | Influence of Mushy Zone Constant on Thermohydraulics of a PCM[END_REF] studied the effect of the mushy zone constant and its influence on the melt fraction, vortex strength and the amount of heat storage. They have shown that increasing C CK leads to a decrease of the convection strength and consequently to a decrease of the heat transfer rate. These studies show the necessity of choosing an appropriate value of this parameter. Because of the semi-solid state and porous nature of the mushy zone, the choice of the value of this constant is still an open problem.

To assess on the influence of the Carman-Kozeny constant, we simulate the melting of a n-octadecane PCM within a square enclosure with the following dimensionless parameters: Ra = 3.27 • 10 5 , Pr = 56.2 and Ste = 0.045. See Sec. 5.2 for a detailed description. Different values of C CK are used ranging from 10 6 to 10 20 . The local liquid fraction function L f (θ) in Eq. (2.23), is regularized as:

L f (θ) = 1 - 1 2 1 + tanh θ f -θ R ck .
(5.

3)

The value of θ ck is set in order to have the sharp variation of the derivative of dL f /dθ near the new phase appearing in the system, i. e. liquid for melting (θ ck = ε), and solid for freezing (θ ck = θ f ) (see also [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]). This is to ensure that the velocity in the solid region is correctly set to u = 0.

The computed velocity in the solid phase is reported in Tab. 5.2 for different values of C CK , ranging from 10 6 to 10 20 . Increasing the constant C CK leads to lower values of the velocity in the solid, as it was expected. However, even though it is generally assumed that a large value for C CK must be set, the exact value of this constant could influence the accuracy of the results. Table 5.2: The influence of the constant of Carman-Kozeny on the values of the velocity in the solid region.

C CK A max u solid v solid 10
We compare the location of the melting front with the experimental data of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]. Figure 5.1a displays the location of the phase-change interface for four values of C CK : 10 6 , 10 10 , 10 15 , and 10 20 . Very good agreement with the experimental result of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] is obtained for C CK = 10 6 . Increasing the constant C CK , results in artificially slowing the propagation of the melting front. This behavior could be explained by the profile of the velocity U = Σ i u 2 i , plotted as a function of temperature θ following the line y = 0.9 (Fig. 5.1b). For a value of C CK = 10 6 , the velocity is zero in the expected temperature interval [-0.01, 0.01], corresponding to the characteristic length of the regularization zone R ck = ε = 0.01. However, imposing too large values of C CK , induces a larger temperature interval in which the velocity is forced to be zero. For C CK = 10 20 , the velocity is canceled up to θ = 0.1, i. e. an artificial mushy zone 10 times larger than the initial case C CK = 10 6 . [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]. (b) Velocity profile following the line y = 0.9 with respect to the temperature.

To obtain the same location of the melting front using different values of C CK , one should thus decrease the smoothing radius parameter R ck . Several values of R ck ranging from ε to ε/4 were tested and the results are plotted in Fig. 5.2. Decreasing R ck by a factor of 2 (Fig. 5.2 a) and 4 (Fig. 5.2 b), results in a merging of the location of the solid-liquid interface for different values of R ck . It is however worth noting that decreasing the regularization interval R ck implies a greater mesh resolution to solve accurately the strong gradient induced by discontinuities at the solid-liquid interface. We thus set for all subsequent simulations C CK = 10 6 and R ck = ε, which corresponds to the best fit with the experimental result of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and the less time-consuming configuration.

Mesh adaptativity

The use of mesh adaptation proved mandatory in the present simulation of PCM melting to obtain accurate results within reasonable computational time. For the melting case, we used five metrics intersection to adapt the mesh, based on S n+1 , S n , T n+1 , T n , and u n+1 . To reduce the impact of the interpolation on the global accuracy for time-depending problems, we consider the metrics computed from actual (at t n+1 ) and previous (at t n ) values, for the same variable used for adaptivity (see also Belhamadia et al. (2004a)).

The mesh adaptation strategy requires to set the interpolation error level errh, the minimum and the maximum edge sizes, hmin and hmax, and the adaptratio parameter, which defines the ratio for a prescribed smoothing of the metrics. The minimum size of triangles and the interpolation error level are set in order to capture the smaller length scale, mainly the boundary layer structure close to the wall, and to solve accurately the large gradient at the solid-liquid interfaces. The error level is adapted for each variable on which the metrics are computed. A default value errh = 0.02 is defined for T and u, and errh = 0.2 for S. The scale for boundary layers in natural convection flow along a vertical wall is given by Eq. (2.25). For high-Prandtl simulations, the thermal boundary layer is thinner than the viscous boundary layer, and conversely for low-Prandtl fluids the viscous boundary layer is thinner than the thermal one. In our simulations, the values of the Rayleigh number, in the range 10 6 to 10 8 , and the Prandtl number of order of 50, lead to a dimensionless thickness of the thermal boundary layer of order of δ θ ∼ 10 -3 . For the same range of Rayleigh numbers and Prandtl number of order of 0.1, we obtain a dimensionless thickness of the viscous boundary layer of δ ν ∼ 10 -2 . The minimum edge-length is therefore set to hmin = 10 -3 . Moreover, we fix the maximum triangle size to avoid generating too large triangles. Mesh adaptivity is performed at each time step and offers a refined discretization of the regularization region where sharp gradients have to be accurately captured. Figure 5.3 shows the adapted mesh during the melting of Gallium PCM in a square cavity. The mesh is remarkably refined around the melting front (Fig 5.3b), localized by the temperature isoline θ = 0. This is also the case in the fluid region around the convection cells, while a coarser mesh is applied in the solid. The typical number of triangles of the generated adaptive mesh is 4, 000. Non-adapted grids offering the same spatial resolution everywhere inside the computational domain would have resulted in N t = 9.94 • 10 10 triangles. Consequently, mesh adaptivity greatly helps in reducing the computational time. The mesh adaptivity capability to capture several interfaces is the striking feature of our method. Its capability to track efficiently two solidification fronts during PCM solidification or the density inversion interface during the water freezing will be discussed in Sec. 5.6.

Initialization

For the melting process, the PCM is initially solid and the temperature is set to a cold temperature θ c , below the temperature of fusion θ f . However, increasing abruptly the boundary temperature to a hot temperature θ h , above the temperature of fusion, in order to initiate the melting process is numerically complicated because of the strong temperature gradient between the wall and the PCM. A first usual approach consists of setting a very thin fluid layer of thickness δx ∼ 0.01 with isothermal temperature in the vicinity of the hot wall (see [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]; Belhamadia et al. (2004a)). A different approach, is applied in the current work. The temperature at the hot wall and the Rayleigh number are increased smoothly over a small time interval to reach the correct values of θ h and Ra. This approach has a physical significance and is more robust when high Rayleigh number simulations are performed.

Concerning the solidification stage, mainly the water freezing case, a 'hot' and a 'cold' restarts are used. The first approach consists of establishing a steady state regime by solving first the steady Eq.

(3.14) and dropping then smoothly the temperature of the cold wall below the temperature of fusion.

The second approach sets directly θ c below the temperature of fusion and considers initial motionless fluid, thus inducing huge gradients near the wall, which requires very small triangles therein.

Melting of a n-octadecane PCM in a square cavity

We start the validation of our algorithm for phase-change systems by simulating the melting of noctadecane PCM within a differentially heated square cavity. Three cases are investigated:

• Benchmark #1: Experimental investigation by [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] of the melting of n-octadecane PCM.

• Benchmark #2: Experimental and numerical investigations by [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice Boltzmann method and performed with CUDA enabled GPU[END_REF] of the melting of a PCM inside a transparent building brick.

• Benchmark #3: Numerical comparison with various numerical methods, presented by [START_REF] Bertrand | Melting driven by natural convection A comparison exercise: first results[END_REF], for the melting of a PCM at high Ra.

The physical parameters used in our simulations are given in Tab. 5.3. The experimental investigation of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] in Benchmark #1 consists of a differentially heated square cavity of dimensions 1.5 cm × 1.5 cm. The left and the right walls are isothermal and the top and the bottom walls are adiabatic. No-slip boundary conditions are prescribed on the velocity over the whole ∂Ω.

Ra

To ensure that we obtain a grid-converged solution, we first perform simulations with different mesh resolutions for Benchmark #1 by keeping the minimum edge length constant and increasing hmax from 10 -2 to 10 -1 . A snapshot of the adapted mesh at t = 78.7 is reported in Fig. 5 For all considered meshes, one can observe the refined grid at the boundary layer region near the walls, due to the high velocity gradients, and along the interface θ = 0. Coarser mesh is applied in the solid region. We note however that a ratio of 10 is set for hmax between Mesh #1 and Mesh #4, resulting into a ratio of 50 for the number of d.o.f. using a P 2 discretization for θ. Simulations are performed using δt = 10 -2 to ensure small time discretization errors. The location of the solid-liquid interface and the time evolution of the Nusselt number related to each of the four adapted meshes are compared in Figs. 5.4e and 5.4f. The maximum difference for the location of the melting front occurs at the top of the cavity since the highest velocity is reached in this region because of the clockwise recirculation of the flow. We also compare the position of the solid-liquid interface along x-direction at y = 1. The coarser grid (Mesh #4) exhibits the highest gap with respect to the most refined Mesh #1, with a relative error of 2.14% for the location of the melting front at the top of the cavity and by overestimating the heat transfer of an order of 3.83%. For Mesh #3, in which the mesh resolution is increased by a factor of 3 compared to the coarse mesh, differences of 0.042% and 0.039% are noticed for both Nu and the interface position. Finally, for Mesh #2, variations lower than 0.01% are obtained. However, since a balance between the accuracy and the computational time is sought, the subsequent simulation are performed using the parameters of Mesh #3, summarized in Tab. 5.4. Solutions using Mesh #3 were performed with only 2876.39 CPU seconds, while 12 CPU hours were needed for Mesh #1.

hmin hmax adaptratio errh nbvx 10 -3 0.1 1.5 0.01 50, 000

Table 5.4: Mesh adaptivity parameters used for Mesh #3. We present in Fig. 5.5 the computed solution for Benchmark #1 case using Mesh #3. Panel (a) displays the temperature distribution of the PCM at dimensionless time t = 78.7. The enclosure is heated from the left and the melting PCM expands from the left to the right (the blue color denotes the solid region). A non-planar shape of the phase-change interface could be observed due to effect of the natural convection in the liquid phase. Further comprehensive descriptions of the influence of the liquid flow on the interface will be discussed in Chapter 6. Panel (b) compares the location of the solid-liquid interface with experimental investigation of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]. For two particular time instants, t = 39.9 and t = 78.7, we could compare our results to available experimental [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and numerical [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]Wang et al., 2010a;[START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] data. In the experimental set-up, the author reported that the top of the PCM was not perfectly insulated and consequently the growth of the upper part of melting front was delayed. The current work agrees well with the experimental results of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] at the bottom part of the melting front. However, our results overestimate the location of the front in the top part of the cavity, which could be related to the experimental heat loss mentioned by the author. Moreover, our results are qualitatively in a better agreement with experimental data than previously published numerical results. This is a direct consequence of the precise tracking of the melting front achieved by the mesh adaptivity performed at each time step. Numerical results for Benchmark #2 are depicted in Fig. 5.6. [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice Boltzmann method and performed with CUDA enabled GPU[END_REF] investigated the melting of octadecane PCM inside a transparent building brick of dimensions 15.2 cm × 3 cm. The numerical simulation is performed using δt = 10 -3 and a total runtime of 19, 231 CPU seconds using 4, 198 triangles. Figure 5.6a illustrates the temperature distribution in the melting PCM, the streamlines, and the phase-change interface obtained by the present finite element simulation. Panel (b) compares the location of the melting front for five particular time instants: t = 0.153, 0.347, 0.465, 0.869, and 1.759, with numerical and experimental data by [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice Boltzmann method and performed with CUDA enabled GPU[END_REF]. Their numerical simulations were based on the thermal lattice Boltzmann method (TLBM), in which the natural convection was solved by LBM and the temperature equation was solved by a finite difference scheme using the enthalpy method. The difficulty here, compared to the first validation case, is the presence of a stronger natural convection flow in the fluid due to the higher value of the Rayleigh number. We notice a very good agreement with both numerical and experimental data of [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice Boltzmann method and performed with CUDA enabled GPU[END_REF].

A last validation case for the n-octadecane PCM is also investigated to test the robustness of the method. [START_REF] Bertrand | Melting driven by natural convection A comparison exercise: first results[END_REF] compiled results provided by five different authors (Lacroix, Le Quéré, Gobin-Vieira, Delannoy and Binnet-Lacroix). Results provided by these authors will be hereafter referred to as (say) 'Lacroix, from [START_REF] Bertrand | Melting driven by natural convection A comparison exercise: first results[END_REF]'. They attempted a first comparison by taking several numerical methods to compute the basic configuration presented in this section. Two investigators among the five failed to predict the process and showed unrealistic behaviors (see Figs. 5.7 and 5.8): Lacroix and Delannoy seem to be insufficiently converged (Fig. 5.7), and Binet-Lacroix overestimates the average Nusselt number by more than 30% (Fig. 5.8). Hence, this collection of results allows us to compare our numerical method and check whether or not realistic results are obtained for complex physical configurations. We inspect the melting front position (Fig. 5.7), the temporal evolution of the liquid fraction L f and the Nusselt number Nu at the left wall (x = 0) (Fig. 5.8), for each of the five results presented by [START_REF] Bertrand | Melting driven by natural convection A comparison exercise: first results[END_REF]. For the liquid fraction, the initial solid state corresponds to L f = 0, while L f = 1 indicates the complete melting of the PCM. The phase-change interface position for four time steps, t = 0.25, t = 1, t = 3 and t = 5 is represented in Fig. 5.7. Our results are for each case in fairly good agreement with those of Gobin and those of Le Quéré. Gobin used a front-tracking method using a coordinate transformation with a finite volume method with a 62 × 42 grid. Le Quéré solved a single domain model using a second order scheme with a finite volume method with a 192 × 192 grid [START_REF] Gobin | Melting from an isothermal vertical wall. Synthesis of numerical comparison exercise[END_REF]). The time evolution of the liquid fraction and the Nusselt number are plotted in Fig. 5.8. A very good agreement is obtained with Gobin and Le Quéré. A relative difference, less than 2% is noticed for the Nusselt number, and a dispersion, smaller than 4%, for the melted fraction.

The high value of the Rayleigh number, Ra = 10 8 , results in a very demanding numerical test. The high velocity, inducing a very narrow thermal boundary layer can lead to unrealistic results and some numerical methods have failed. The interest of the mesh adaptation is clearly evidenced since a typical number of triangles of N t = 7, 000 was used during the simulation, which required 40, 522 CPU seconds.

Melting of Gallium in a rectangular cavity

The melting of tin or Gallium in a rectangular cavity, which corresponds to low-Prandtl fluids, was a controversial case since [START_REF] Dantzig | Modelling liquid-solid phase changes with melt convection[END_REF] raised the question whether the convection in the fluid is mono-cellular or multi-cellular. The experimental result exhibits indeed a mono-cellular structure, while many researchers claim this observation to be incorrect. Prior to [START_REF] Dantzig | Modelling liquid-solid phase changes with melt convection[END_REF] note, both experimental and numerical results supported a single cell solution in the fluid phase. Later, simulations provided solutions with multicellular flow. During the comparison exercise investigated by [START_REF] Bertrand | Melting driven by natural convection A comparison exercise: first results[END_REF], Le Quéré and Couturier-Sadat have predicted a different shape of the front corresponding to low-Prandtl simulation, induced by the multi-cellular structures of the flow in the melted PCM. Le Quéré and Gobin (1999) showed that the multi-cellular structures result from the hydrodynamic instabilities during the conduction regime before the onset of convection. Moreover, the numerical investigation of [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF], testing the influence of the mesh resolution and the numerical method's order, indicates clearly that mono-cellular observation is caused by a problem of convergence of the numerical solution. It can be due to the grid size or inconsistencies in the mathematical model.

Therefore, this test case simulating the melting of Gallium is a relevant exercice to test the consistency of our method. To capture the very small cell during the first step of the melting, [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF] used a 800 × 1, 120 fixed grid in a rectangular domain of dimensions 6.35 cm × 8.89 cm. The authors have reported that the melting solution up to 32 s required about 980 CPU hours on a Compaq Alpha (667 MHz, ev67) processor. A maximum of 4, 820 triangles were necessary with our adaptive method, and only 4 CPU hours were spent to reproduce the numerical result of [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF] up to the final time 280s, on an Intel (2.8 GHz, Core i7) processor.

The time evolution of the flow is presented in Fig. 5.9. The enclosure is heated from the left (vertical) wall and the horizontal walls are adiabatic and no-slip. Dimensionless parameters of the run and properties of Gallium are given in Tabs. 5.1 and 5.3. Temperature field, streamlines and position of the melting front are displayed for several time instants: t = 0.0015, 0.006, 0.01, and 0.019. These values were chosen to visualize the merging of convection cells in the fluid flow and correspond to physical times 20s, 85s, 155s, 280s in [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF]. The number of rolls was considered as a validation criterion by several authors [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF][START_REF] Cerimele | Numerical modelling of liquid/solid phase transitions: Analysis of a gallium melting test[END_REF]Giangi and Stella, 2000). Three cells are observed at t = 0.006 (Fig. 5.9). The number of cells decreases later through a process of roll merging, as it was also reported by Le Quéré and Gobin (1999) and [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF]. Our numerical results are in good agreement with the observations of [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF], [START_REF] Cerimele | Numerical modelling of liquid/solid phase transitions: Analysis of a gallium melting test[END_REF] and Giangi and Stella (2000). The adapted mesh corresponding to t = 0.019 is illustrated in Fig. 5.9e. The mesh is well refined along the isoline θ = 0 denoting the solid-liquid interface, and around the convective cells.

The time evolution of the Nusselt number at the heated vertical wall is plotted in Fig. 5.9f. One can note that the quasi-steady evolution of Nu, observed in Fig. 5.8 for high-Prandtl fluid simulations, is not recovered. This is in agreement with the stability analysis of Le Quéré and Gobin (1999), showing that the Nusselt number becomes unsteady starting from a critical time, for Rayleigh numbers greater than 10 5 and low-Prandtl numbers, while no instability could be found whatever the value of Ra for high-Pr cases. The first slight oscillation at the earlier time steps corresponds to the onset of the convective flow induced by the first instability of the conductive regime. It is followed by a nearly constant evolution of Nu from t = 0.01 to t = 0.02 during the cell-merging process. High oscillations of the heat transfer are then observed, which is a feature of the onset of the second oscillatory instability. A zoom on the periodic evolution of the Nusselt number is displayed in Panel (f). Finally, the amount of heat transfer decreases starting with t = 0.13, when the front reaches the cold (right) vertical wall.

Melting of a cylindrical PCM with inner heated tubes

Cases presented in Secs. 5.3 and 5.2 considered phase change problems evolving in a simple geometry, a rectangular box. A more complex geometry, suggested by [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF], is simulated in this section. It consists of a cylindrical PCM of radius R = 1 with tube inclusions of different arrangements.

The interest in studying this case is not solely the challenge of the complex configuration, but also the possibility to compare our results with those of [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF], obtained using a completely different model based on the Lattice Boltzmann Method. This configuration is also interesting from a practical point of view. [START_REF] Agyenim | A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (LHTESS)[END_REF] pointed out that more than 70% of the PCM containers used for heat storage are using shell-tube systems.

We simulate three cases with one heated tube, four heated tubes, and nine heated tubes with the same total tube area. For the one heated tube case, the radius R i of the inner tube is one-quarter of 56 For a better view of the convection cells, a ratio 2 : 1 was used for the axes. Figure 5.10 shows the temperature field and the position of the solid-liquid interface (black line) related to the three configurations for time instants corresponding to the same liquid fraction L f = 80%. The distribution of the inner tubes in the liquid phase influences directly the fluid motion and the shape of the melting front. The more the number of inner tubes, the stronger the natural convection is in the melted PCM. The shape of the solid-liquid interface displays complex patterns, depending on the space arrangement of the inner tubes. This is linked to the effect of the fluid motion in the presence of obstacles, as it was also noticed in Sec. 4.2. The mesh is adapted following the evolution of the melting interface, even after its separation in several distinct fronts touching the outer boundary (see Figs. 5.10b,c).

Melting of a cylindrical PCM with inner heated tubes

To estimate the efficiency of each configuration, we plot in Fig. 5.11 the time evolution of the liquid fraction L f . By including more heated tubes, the heat transfer is enhanced, inducing a faster melting time. The nine-tube configuration melts 5 times faster than the reference configuration with one tube. A reasonably close match with numerical results of [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF] is confirmed in Fig. 5.11a. We also assess on the effect of tubes arrangements for a fixed number of tubes. Three arrangements are simulated with 9 tubes and reported in Fig. 5.11b: central symmetry, inline, and staggered arrangements. The arrangement with central symmetry appears to provide a better heat transfer compared to the others. The inline and the staggered arrangements exhibit similar trends, with a slightly faster melting rate than the inline arrangement. The difference is explained by [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF] by the contribution of conductive heat transfer induced by the arrangements.

Solid crust formation in a highly distorted mesh

The solid crust formation inside a highly distorted domain, simulated by [START_REF] Nourgaliev | Fullyimplicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change[END_REF] is considered in this section. Our emphasis here is on the ability of our method to tackle also solidification problems, especially in irregular domain (see Fig. 5.12).

The fluid is initially motionless with an initial dimensionless temperature θ 0 = 2, above the temperature of fusion. The temperature of fusion is set to θ f = 1.4, according to [START_REF] Nourgaliev | Fullyimplicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change[END_REF] parameters. It is worth noting that [START_REF] Nourgaliev | Fullyimplicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change[END_REF] used T ref = T f thus θ f = 0. The left side is set at cold temperature θ c = 1.39 in the initial stage while the right wall was kept constant at a hot temperature θ h = 2, so that a steady-state natural circulation should be established. The cold temperature at the left wall is then dropped smoothly to θ c = 1, below the temperature of fusion, starting the formation of a solid crust layer. The top and bottom walls are adiabatic.

The dimensionless parameters of the simulation are: Ra = 1 • 10 6 , Pr = 0.1 and Ste = 4.854. The scaling (4.2), which was also used by [START_REF] Nourgaliev | Fullyimplicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change[END_REF], is used in the current simulation.

The temperature distribution, the streamlines, and the melting front position are reported in Fig. 5.12a and compared with the numerical results of [START_REF] Nourgaliev | Fullyimplicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change[END_REF]. These authors used a second-order discontinuous Galerkin finite element method and a viscosity-based technique to bring the velocity to zero in the solid region. We are qualitatively in good agreement with [START_REF] Nourgaliev | Fullyimplicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change[END_REF]. As expected, the high Ra and Ste dimensionless numbers, combined with the distorted domain, generate a complex vortical flow pattern in the melted PCM. [START_REF] Nourgaliev | Fullyimplicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change[END_REF] presented a mesh convergence analysis and concluded that a 512 × 256 mesh resolution is needed to get a full resolution of the five vortical structures in Fig. 5.12b. Only 2, 769 triangles are used in the present simulation. This case illustrates again the capabilities of our mesh adaptivity algorithm.

Water freezing

We finally consider the difficult case of water freezing in a square cavity. After achieving the convection steady pattern in the cavity (see Fig. Besides the complex fluid flow and the non-linear time evolution of the solid-liquid interface, the code have to handle the non-linear variation of the density presented in Eq. (4.4). Moreover, one can note that in addition to the non-linear formulation of f B (θ) in Eq. (4.5), the coefficient Ra × 1/(βδT ) is very large since the Boussinesq approximation imposes a value of βδT ≪ 1. Since the assumption of constant conductivity is not accurate in the frame of water phase-change, K(θ) is regularized as:

K(T ) = 1 + 1 2 k s k l -1 1 + tanh T f -T R k , (5.4)
with k s = 2.26 and k l = 0.578 and R k = 0.0075. We also define ε 1 = 0 and ε 2 = 0.0075, to ensure a very thin mushy zone (see also [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]).

60

Water freezing

Figure 5.13a superimposes the experimental image from [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF] with our numerical results for the same physical time t ϕ = 2340[s]. The flow pattern in the liquid phase also corresponds very well qualitatively to the experimental image. Since a good agreement with the experiment was sought, the simulation was performed with very small time steps (δt = 10 -5 ≈ 0.014[s]), but still reasonable grids (2, 500 triangles) were generated due to the efficiency of the adaptivity algorithm. The discrepancy between the experimental interface and the current simulation comes principally from the model, which neglects the three-dimensional effects, the supercooling of water and realistic boundary conditions (see Giangi et al. (2000); [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]; [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF]). The source term S(θ) is added to the metrics calculation presented in Sec. 4.2 to track the solidification front. The mesh is thus refined along the line θ m through the function Φ(θ) defined in Eq. (4.7), the solid-liquid interface defined by θ = 0 and in the boundary layer regions. We use both Φ(θ n ) and Φ(θ n+1 ) and S(θ n ) and S(θ n+1 ) in the adaptivity procedure. The final mesh is displayed in Fig. 5.13c, clearly showing that the mesh is refined along the line θ = θ m and the solid-liquid interface (θ = 0). This allows to accurately capture the structure and the extent of the two recirculating zones, features that are difficult to obtain with fixed meshes (see discrepancies described in Giangi et al. (2000); [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]; [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF]).

The time evolution of the solid-liquid interface is depicted in Fig. 5.13d for several time instants from t = 0.05 (t ϕ = 70[s]) to final time t = 1.61 (2340[s]). We observe first a straight shape of the front from t = 0 to t = 0.05. Then, the top and bottom of the freezing front move at different rates from t = 0.3 due to the competing effects of positive and negative buoyancy forces in the liquid water. The upper clockwise circulation has indeed the task to transport the hot liquid to the top wall and back along the extremum of the density variation (localized at the temperature isoline θ = θ m ), and slows down consequently the growth of the solid layer. In the same time, the abnormal recirculation traps the cold liquid water at the bottom part of the cavity, enhancing thus the solidification rate at this region as one can note in Fig. 5.13d.

Concluding remarks

In this chapter, we demonstrated the capability of our numerical method to produce highly accurate solutions of difficult multiphysics problems on several geometrical configurations. The complexity of the problem comes from the non-linear time evolution of the solid-liquid interface that could be highly deformed by the strong convection flow in the fluid phase. Linear and non-linear expressions of the buoyancy force were investigated.

For a linear Boussinesq force f B , melting of pure paraffin (n-octadecane) and metal (Gallium) were simulated in a rectangular cavity heated from the side. The choice of these materials is motivated by their physical properties relatively equal in both solid and liquid phases, and the existence of many numerical and experimental investigations in the literature. Differences between the time evolution of the heat transfer for high and low Pr configurations were investigated. Comparisons with existing benchmarks showed for each case a very good agreement. The unsteadiness of the flow for the Gallium was observed in our simulations by showing multi-cellular structure of the liquid flow and the high oscillating time evolution of the Nusselt number. The performance of the sequential code was also proven by simulating complex geometrical configurations. The melting of a cylindrical PCM including inner heated pipes and the solidification of a liquid PCM inside a distorted cavity showed good agreement with existing numerical data.

The total CPU time and the typical number of triangles corresponding to each PCM cases are summarized in Tab. 5.5. Noticeable is the fact that less than 2 CPU hours are necessary for most of the simulations, even when a sequential algorithm is used. The interest and the power of mesh adaptivity was highlighted by the simulation of the melting of the Gallium, since only 4 CPU hours are necessary to perform this case, while 980 CPU hours were reported by [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF] A non-linear expression of f B was considered for the solidification of pure water. Besides the high Ra and Ste considered in our simulations, the efficiency of the mesh adaptivity was demonstrated by its capability to track simultaneously the interface θ = θ f and the line θ = θ m separating the two recirculating liquid regions. Qualitative comparison with the experimental image exhibited good agreement. The vortical structures in the liquid water were accurately captured and the density inversion interface was precisely solved.

A more comprehensive description of the melting process will be presented in the next chapter. Numerical and analytical tools will be developped to analyse and compare the melting of PCM heated either from the side or from below.

Numerical comparison of basal and lateral melting of a phase change material

We performed extensive validations of our numerical method in Chapters 4 and 5. Very good agreements were obtained for all validation tests against well-known benchmark cases. The robustness of the algorithm was proven by simulating challenging configurations. We now use the code as an investigation tool to analyse the phase-change process during the melting stage. We consider a square cavity of height H filled with n-octadecane PCM and pay a closer attention to the temporal evolution of different physical parameters of the system. Two classes of convective melting systems are of interest in this chapter: (i) lateral melting (LM) or (ii) basal melting (BM). As far as (i) is concerned, the PCM is subject to heating from the left side of the cavity, whereas for BM case the PCM is heated from the bottom generating a top melting boundary. The comparison of both cases is interesting since the dynamics of the melting is known to be fundamentally different for each of the two cases. The investigated configuration, the shape of the interface and the streamlines in the liquid phase are illustrated in Fig. 6.1. LM case (top) exhibits a mono cellular pattern while Rayleigh-Bénard-like convection cells are observed in the BM case (bottom). LM case could be representative for building applications, solar collectors or other thermal energy storage systems. In the frame of building applications, bricks made of PCM melt because of differences between outdoor and indoor temperatures and store the energy in the form of latent heat. [START_REF] Barreneche | In situ thermal and acoustic performance and environmental impact of the introduction of a shape-stabilized PCM layer for building applications[END_REF] showed that a wall made of PCMs reduces the temperature peak of about 20%. With regard to LM configurations, analytical investigation of [START_REF] Bejan | Analysis of melting by natural convection in an enclosure[END_REF] and scaling analysis of [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] permitted to describe the heat transfer during the melting by the mean of N u-Ra correlation. Their approach will be used in the present work. The BM case refers to passive temperature control for electronic devices or for a long list of geophysical problems, such as lava lakes, thermal convection in magma chambers, or ice-melt lakes. In contrast with LM case, notwithstanding linear and weakly nonlinear instability analysis based on a vanishingly small Stefan number assumption [START_REF] Vasil | Dynamic bifurcations and pattern formation in melting-boundary convection[END_REF], exact theoretical expressions of N u and L f are not available during the convective regime because of the important non-linearities in the dynamics. Some comparisons with theory [START_REF] Malkus | The heat transport and spectrum of thermal turbulence[END_REF][START_REF] Grossmann | Scaling in thermal convection: a unifying theory[END_REF] made in the frame of Rayleigh-Bénard convection flow have been however carried out [START_REF] Esfahani | Basal melting driven by turbulent thermal convection[END_REF][START_REF] Madruga | Dynamic of plumes and scaling during the melting of a Phase Change Material heated from below[END_REF][START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF]. In this work, we develop a scale analysis for the so-called 'linear regime', which occurs between the onset of convective and oscillating flows.

In this chapter, we offer a comprehensive comparison between LM and BM configurations. We analyse first the time evolution of the LM process through a scale analysis in Sec. 6.1. Second, the BM case is studied in Sec. 6.2 with theoretical descriptions of the heat transfer involved during the melting. Finally, a comparison between the two cases is presented in Sec. 6.3 where we suggest some practical applications based on our analysis.

Lateral melting of n-octadecane PCM. Case LM.

We consider the physical properties of n-octadecane given in Tab. 5.1 and investigate different heights H of the cavity and different values of δT , to assess on the influence of the Ra number. The numerical configuration is sketched in Fig. 6.1a. Ra numbers ranging from 3.27•10 5 to 3.27•10 6 are investigated. We note that for the range of interest for Ra and Ste numbers, the assumption that (β × δT ) ≪ 0.01 for the Boussinesq approximation is satisfied. A second dimensionless time τ related to the analytical correlation of [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] is introduced:

τ = Ste × Fo = Ste × αt ϕ H 2 = Ste Pr × t, (6.1)
where Fo is the Fourier number.

Analysis of the time evolution of the melting process

We start by describing the time evolution of the melting process for the lowest value of the Stefan and Rayleigh numbers, i. e. Ste = 0.045 and Ra = 3.27 • 10 5 . We are interested in a slow melting of the PCM to capture the transitions between the regimes described by [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF], mainly the onset of the convective regime.

At τ = 0, the material is solid and the initial temperature is set to θ 0 = -0.01 everywhere inside the cavity. Then, the temperature of the left wall is progressively increased to θ h = 1, while the right wall is maintained at the same cold temperature θ c = -0.01. The material starts to melt, with a melting front (identified by the iso-line θ = θ f = 0) propagating from the left to the right side of the domain. Snapshots of the time evolution of the phase-change system are given in Fig. 6.2 for representative We can easily identify three different regimes describing the time evolution of the melting process.

• From τ = 0 to τ = 0.004 (Fig. 6.2a), we note the vertical shape of the melting front, well predicted by the classical conduction model of [START_REF] Stefan | Über die Theorie der Eisbildung, insbesondere über die Eisbildung im Polarmeere[END_REF]. This indicates that, at this stage, heat transfer is dominated solely by conduction.

• Between τ = 0.016 to τ = 0.032 (Fig. 6.2b), the natural convection in the fluid phase starts to alter the shape of the melting front. A mixed conduction and convection regimes rule the heat transfer. Convective flow mainly affects the upper part of the fluid motion, while conduction is still dominating in the lower part. As the volume thermal expansion coefficient β is positive, we expect a clockwise circulation of the liquid inside the convection cell, as noted by [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF]. This also makes the liquid-solid interface to move faster at the top of the cavity, explaining the deformed shape of the melting front, which is a signature of the convection effects (see also [START_REF] Kowalewski | Phase change with convection: modelling and validation[END_REF]).

• After τ = 0.032 (Fig. 6.2c-d), natural convection dominates the heat transfer process and impacts radically the solid-liquid interface shape and motion. The melting front line exhibits four distinct regions characterized by different slopes with respect to the vertical axis. The largest slope is observed at the top of the cavity and is related to the particular shape of the convection cell. Note that top and bottom parts of the interface are normal to the cavity boundaries because of the imposed adiabatic boundary conditions.

• After τ = 0.08 the melting front is nearly touching the right wall of the cavity, firstly at the top (Fig. 6.2e) of the cavity. The melting process continues and the fluid progressively fills the cavity, with a melting front deforming to a vertical line. The simulation of the melting process is stopped at τ = 0.2 (Fig. 6.2f), when it is numerically difficult to separate the melting front from the right wall boundary. At this time instant, the fluid fraction reaches the value of 0.95 and the melting of the PCM is considered to be complete, even though a small region of solid PCM remains at the lower right bottom of the cavity. Note from Fig. 6.2f the existence in the fluid of two recirculating zones instead of a single one observed during previous stages.

Scale analysis of the melting

We further analyse each of the three regimes cited previously and identify the proper scales of the phenomenon. The (dimensionless) location of the interface will be denoted by Γ i . Immediately after the melting starts (see Fig. 6.2a), the melted PCM occupies a thin enclosure of height H and width Γ i . In this configuration, the temperature varies linearly between the two sidewalls and the heat transfer is essentially ruled by conduction. The fluid phase is almost motionless and the horizontal heat fluxes across the incipient melting PCM is balanced by the enthalpy absorbed at the interface. At the solidliquid interface, the energy balance condition which takes into account the released latent heat and the discontinuity of heat fluxes between the solid and the liquid can be expressed by the following Stefan condition:

θ(x = Γ i ) = θ f = 0, (6.2) ∂Γ i ∂t = - Ste Pr ∂θ ∂x .
Since the temperature field during the conduction regime is quasi-steady, Eq. ( 6.2) which is linear could be approximated by

∂Γ i ∂t ≈ - Ste Pr θ f -θ h Γ i ≈ Ste Pr 1 Γ i . (6.3)
The location Γ i of the interface is consequently given by

Γ i = 2 × Ste Pr t = √ 2τ . (6.4)
Moreover, the Nusselt number can be evaluated using the same assumption:

N u = 1 0 ∂θ ∂x x=0 dy = 1 Γ i = (2τ ) -1/2 . (6.5)
To summarize, during the very first stage of the melting, when the heat transfer is led by conduction, the time evolution of the liquid fraction (depending on the location of the interface) and the Nusselt number could be approximated as:

L f ∼ (2τ ) 1/2 , (6.6) N u ∼ (2τ ) -1/2 . (6.7)
While the melting continues to expand to the right side of the domain, a natural convection flow emerges from the top of the cavity since we have a clockwise recirculation of the flow in the liquid phase (see Fig. 6.3). Turning back to the energy Eq. (2.14) in the liquid phase (i. e. C = K = 1), three distinct effects could be identified:

δθ t Inertia , v δθ H Convection , δθ Γ 2 i Conduction . (6.8)
As t increases, the inertia decreases while the convection effect increases since it is proportional to v.

The conduction becomes more and more negligible with the increasing value of Γ i with time. [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] described the convective heat transfer contributing during this regime by defining a Rayleigh number based on y * as Ra y * = Ra × y * 3 , with y * the height of the melted zone altered by the convection flow as shown in Fig. 6.3. At the bottom part of the cavity, the interface remains vertical by the effect of the conductive heat transfer. The total Nusselt number during the mixed conductionconvection regime could be therefore approximated by [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF]:

N u ∼ (2τ ) -1/2 + (2τ ) 3/2 × Ra.
(6.9) Equation (6.9) indicates that the contribution of the conduction ∼ 1/ √ 2τ decreases with the time while the convection contribution is increasing.

Finally, when the natural convection flow is fully developed and dominates the heat transfer along the vertical heated wall, using the correlation of [START_REF] Bejan | Convection heat transfer[END_REF], the dimensionless thickness of thermal boundary layer is δ θ ∼ Ra 1/4 (see also [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF]) and therefore the Nusselt number scales as:

N u ∼ Ra 1/4 . (6.10) Jany and [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] have proposed a more general single correlation, combining the regimes described previously:

N u(τ ) = 1 √ 2τ + c 1 Ra 1/4 - 1 √ 2τ 1 + c 2 Ra 3/4 τ 3/2 n 1/n . (6.11)
The values of the constants were fitted from numerical data: c 1 = 0.27, c 2 = 0.0275, and n = -2. [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] has also suggested from his experimental data the following correlation for N u:

N u =              1 √ 2τ , if τ ≤ τ t , 1 √ 2τ t {1 + C(τ -τ t )} , if τ > τ t ,
c 1 Ra 0.266 , otherwise, (6.12)

with the constant c 1 = 0.234 and the exponent 0.266 fitted from experimental data, and τ t the transition time from conduction to convection as discussed previously.

Predictions of Eqs. (6.12) and (6.11) are compared with our numerical results in Fig. 6.4 showing the time evolution of the Nusselt number at the left wall. Our results perfectly fit with the theoretical prediction of [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] and are also in good agreement with the experimental correlation of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]. The gap between the current simulation and the results of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] could be explained by the experimental heat loss mentioned by the author and the uncertainties of the experimental measurements. The regimes described by the shape of the interface in Sec. 6.1.1 could be followed in the temporal evolution of N u:

1. The pure conduction regime (N u ∼ (2τ ) -1/2 ) for τ 0 to τ ∼ Ra -1/2 = 0.02 (corresponding to Fig. stage of the flow evolution. The signature of this conduction regime is the slow heat transfer characterized by a monotonic decrease of the Nusselt number.

2. The mixed conduction-convection regime (N u ∼ τ -1/2 + Ra τ 3/2 ) for 0.02 ≤ τ ≤ 0.05 (illustrated in Fig. 6.2b).

3. The convection dominated regime (N u ∼ Ra 1/4 ) for τ > Ra -1/2 (corresponding to Figs.

6.2c-e). The plateau at the value of Ra 1/4 corresponds to the pure convective transfer and is observed in Fig. 6.4 for 0.05 ≤ τ ≤ 0.1. Numerical results show a slight decrease of N u in the final stage (τ ≥ 0.1), when the melting front starts to touch the right wall of the cavity (see Figs. 6.2e-f). The correlation model is not valid for this late evolution of the melting process.

Another important basic quantity describing the melting process is the liquid fraction L f . The time evolution of the liquid fraction (Fig. 6.5a) displays three regimes during the melting process. L f initially grows as τ 1/2 , which is a typical law for a conduction-dominated heat transfer. Then, a linear temporal evolution is observed, until the melting front reaches the right wall. This linear regime corresponds to the quasi-steady state observed in the evolution of the Nusselt number (Fig. 6.4).

Using the asymptotic limits of Eq. ( 6.11) for τ → 0 (pure conduction) and τ → ∞ (pure convection), [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] suggested the following correlation law for the time evolution of the liquid fraction:

L f (τ ) = √ 2τ 5 + c 1 Ra 1/4 τ 5 1/5
, (6.13) where c 1 = 0.27 is the same constant as in Eq. (6.11). We compare in Fig. 6.5b our numerical results with the predictions based on Eq. ( 6.13) within the validity domain of the analysis, i. e. before the melting front reaches the right wall of the cavity. A very good agreement is found with theoretical predictions and also with previously published numerical results (Wang et al., 2010a).

Influence of the Rayleigh number

To investigate the influence of the Rayleigh number on the evolution of the melting process, we performed different simulations by multiplying the initial value of Ra = 3. Panels (a) to (c) depict the dynamics of the convective melting flow when half of the initial solid PCM (L f = 0.5) have melted. The top part of the interface moves faster while the bottom one is slowed by the increasing value of Ra. According to the Stefan interface condition in Eq. (6.2), at constant Ste and Pr numbers, the interface velocity is proportional to ∂θ/∂n, which is maximum at the top of the cavity because of the clockwise recirculation of the fluid. This explains the observed trends. Figures 6.7 and 6.8 plot the temporal evolution of the liquid fraction L f (panel a), and the average Nusselt number defined at the hot wall (panel b). The same heat transfer regimes described previously are observed for each case: conduction, mixed conduction-convection and convection. We note that results are plotted with respect to physical time t ϕ instead of τ , because we compare solutions with different values of H. This infers with the definition of the non-dimensional times t and τ , making them not relevant to compare solutions. Figure 6.7a indicates that increasing the Rayleigh number by keeping δT constant induces a slower melting rate. This is the expected behaviour since the size of the PCM is increased by a factor of 2, and the velocity u is decreasing to satisfy the condition Re = 1. We note however a non-monotonic variation of the time necessary to melt a fixed value of PCM. For instance, to achieve L f = 0.5 (50% of the volume is melted), an increase of Ra by a factor of 10 leads to a growth of the time by a factor of 1.7. Nonetheless, when Ra is 5 times larger, the necessary time only increases by a factory of 2. This is most likely due to the non-linear intricacies of the problem and requires further investigation. Furthermore, the Nusselt number reported in Fig. 6.7b shows that the higher the Rayleigh number, the higher is the Nusselt number. This is consistent, since the temperature gradient is integrated along a longer heated wall. .8a shows that by increasing the value of δT , and consequently increasing the Rayleigh number and the Stefan number, the PCM melts faster. We note that the height H of the cavity is kept constant, hence the natural convection flow in the melted PCM is enhanced when the Rayleigh number keep increasing. As a consequence, the convection-dominated regime is reached earlier, as shown in Fig. 6.8b by the shift of the minimum of the N u to lower values of t ϕ . This evolution is also observed for the liquid fraction. As expected, an increase of the Rayleigh number and the Stefan number is followed by an enhancement of the heat transfer during the melting, and consequently an improved efficiency of the PCM.

Basal Melting of n-octadecane PCM. Case BM.

We now pay attention to the BM case. We investigate the melting of pure n-octadecane PCM in a square enclosure subject to heating from the bottom side. When compared to the lateral melting case, the dynamics of the melting is different for such configuration, in which natural convection develops Bénard cells. We recall that the physical parameters and the numerical configuration are reported in Fig. 6.1. We perform two-dimensional numerical simulations, although [START_REF] Gau | Flow visualization during solid-liquid phase change heat transfer II. Melting in a rectangular cavity[END_REF] and [START_REF] Gong | Flow and heat transfer in convection-dominated melting in a rectangular cavity heated from below[END_REF] noticed the existence of three-dimensional convection cells during the very first step of the melting process. These three-dimensional convection cells are usually neglected for relatively moderate Ra numbers, particularly for Ra ≤ 10 8 . In this case, three-dimensional cells survive over a very short time compared to the whole melting process. Therefore the two-dimensional model is realistic. A qualitative description of the dynamics of the natural convection flow and its impact on the melting front is first addressed in Sec. 6.2.1. Then, a scale analysis is conducted to describe the heat transfer that occurs during the melting in Sec. 6.2.2. Finally, a comparison between the LM and BM cases is presented in Sec. 6.3.

Temporal evolution of the melting process

Figure 6.9 displays the structure of the natural convection flow after the primary instability through a sequence of panels for temperature isolines and streamlines in the liquid phase, for Rayleigh numbers ranging from Ra = 3.27 • 10 5 to Ra = 3.27 • 10 6 . An array of lengthening plumes (panels a to c) and counter-rotating convective cells (panels d to f) are located in the liquid phase, in which the number of thermal plumes increases with the Rayleigh number. For Ra = 3.27 • 10 5 , one can observe three equidistant plumes in Fig. 6.9a and five convective cells in Fig. 6.9d, while four and six plumes are observed for Ra = 1.635 • 10 6 in Fig. 6.9b and Ra = 3.27 • 10 6 in Fig. 6.9c. These observations agree well with the numerical results of [START_REF] Gong | Flow and heat transfer in convection-dominated melting in a rectangular cavity heated from below[END_REF] and [START_REF] Madruga | Dynamic of plumes and scaling during the melting of a Phase Change Material heated from below[END_REF] who studied the correlation between the number of thermal plumes and the size of the domain. The shape of the interface is directly linked to the dynamics of these plumes. The mushroom form of the plumes results from the two symmetric counter-rotating convective cells surrounding each of them. We observe an anti-clockwise recirculation of the left convection cell and a clockwise recirculation of the right one. Thus, the liquid is heated to the highest temperature at the bottom and then floats up, reaches the phase change interface and splits into opposite directions. The liquid is cooled as it flows through the phase-change interface. It results a non-planar front with a peak at the center of each couple of counter-rotating convective cells. To describe the temporal evolution of the melting, it is useful to introduce the effective Rayleigh and Nusselt numbers of the fluid layer, based on the height of the melted PCM:

Ra e = Ra × δ3 H , (6.14) N u e = N u × δH , (6.15)
with δH the non-dimensional averaged fluid height. Note that δH could be assimilated to the liquid fraction.

Figure 6.10 depicts in detail the temporal evolution of the melting for Ra = 6.54•10 6 , during which the effective Rayleigh number increases and influences the dynamics of the flow. The temperature field, the location of the interface, and the streamlines are reported in panels (a) to (l).

Before the first instability arises, the melted layer evolves solely by conduction. There is no noticeable fluid flow and the melting front remains straight (panel a). The temperature is distributed linearly in the vertical direction through the incipient liquid layer. The convective regime occurs at Ra e ≈ 7 × 10 3 when the phase-change interface becomes non-planar (panel b). This is in good agreement with the observation of [START_REF] Esfahani | Basal melting driven by turbulent thermal convection[END_REF][START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF]. It is worth noting that, in the framework of the classical Rayleigh-Bénard convection flow, the first instability appears at critical Rayleigh number Ra c ≈ 1707.76, in the limit of vanishing Ste number [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF]. This critical value is however higher for large Ste numbers [START_REF] Vasil | Dynamic bifurcations and pattern formation in melting-boundary convection[END_REF]. The onset of convection is marked by a change in the shape of the interface from straight to a nearly periodic curve.

As the fluid expands upwards, the effective Rayleigh number increases and we observe 10 convective rolls being stretched vertically (panels c to e). For Ra e ≤ 2 • 10 5 , one can note that the number of rolls is time-independent. Such behaviour could be compared to the steady convection regime after the onset of primary instability in the Rayleigh-Bénard system (see [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF]). This regime will be referred to linear regime in the present study.

After the rolls are elongated vertically, they start to oscillate laterally and then merge to create greater rolls, starting with Ra e = 2.9 × 10 5 (panel g). The interface loses periodicity and the structure of the rolls becomes disordered. At this stage, neither the mushroom form of the plumes nor the periodic distribution of the convection cells are no longer observed. The main consequence of this observation is that the melting front is modified. The interface is shaped by the new flow pattern.

Scale analysis

Let us now proceed on a scale analysis for the description of the heat transfer processes during melting. The dynamics of the BM case is usually classified in the literature into five regimes, defined following a power law of the Nu -Ra [START_REF] Esfahani | Basal melting driven by turbulent thermal convection[END_REF][START_REF] Madruga | Dynamic of plumes and scaling during the melting of a Phase Change Material heated from below[END_REF]): a conductive regime, a linear regime, an oscillating regime, a turbulent regime, and finally an ultimate regime. The simulations performed in this chapter cover the first three regimes. During the conductive regime, the heat transfer is fully dominated by conduction. The time evolution of the liquid fraction and the Nusselt number could be approximated by the same scaling obtained in Eq. (6.6). Notwithstanding, after the onset of the convection, the quasi-steady evolution of the heat transfer observed in the LM case (see Fig. 6.4), is no longer observed in the BM case because of non-linear instabilities.

When the convective heat transfer is fully developed in the melted PCM, two distinct heat transfer processes could be identified: a bulk heat transfer for low effective Rayleigh number regime, mainly for Ra e ≤ 10 5 , followed by a boundary layer heat transfer regime, which is predominant when the flow becomes oscillatory or turbulent. We develop in this section a scale analysis for the linear regime, during Figure 6.10: Time evolution of the melting of a PCM heated from below at Ra = 6.54 • 10 6 : temperature field, solid-liquid interface and streamlines. For each time instant, the effective Rayleigh number is indicated.

which bulk heat transfer occurs, to approximate the amount of heat transfer. Beyond this regime, due to the important non-linearities of the dynamics, we rely on theories developed for turbulent Rayleigh-Bénard convection flows [START_REF] Malkus | The heat transport and spectrum of thermal turbulence[END_REF][START_REF] Grossmann | Scaling in thermal convection: a unifying theory[END_REF].

During the linear regime, the convective cells are elongated vertically and the number of rolls is time-independent. In our scale analysis we focus on a single roll convection pattern to assess the amount of heat transfer through it. A schematic overview of the phenomenon is drawn in Fig. 6.11. Let us define L r , the (dimensionless) half-thickness of a cell. An appropriate scaling during this regime is hence

x ∼ L r , y ∼ δH , (6.16) with L r < δH . We consider the dimensionless Navier-Stokes-Boussinesq system (2.12) -(2.14) applied to a single cell (C = K = 1, and S(θ) = A(θ) = 0), which could be rewritten as follows (with scaling (4.6), i. e. Re = 1):

∂u ∂x + ∂v ∂y = 0, (6.17) ∂u ∂t + u ∂u ∂x + v ∂u ∂y = - ∂p ∂x + ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 , (6.18) ∂v ∂t + u ∂v ∂x + v ∂v ∂y = - ∂p ∂y + ∂ 2 v ∂x 2 + ∂ 2 v ∂y 2 + Ra Pr θ, (6.19) ∂θ ∂t + u ∂θ ∂x + v ∂θ ∂y = 1 Pr ∂ 2 θ ∂x 2 + ∂ 2 θ ∂y 2 . (6.20)
The pressure p is first eliminated by deriving Eq. ( 6.18) with respect to y and deriving Eq. ( 6.19) with respect to x. One may refer to the book [START_REF] Bejan | Convection heat transfer[END_REF] for more details, in the frame of the Rayleigh-Bénard convection flow. Substracting one from the other, we get

∂ ∂y ∂u ∂t + u ∂u ∂x + v ∂u ∂y - ∂ ∂x ∂v ∂t + u ∂v ∂x + v ∂v ∂y (6.21) = ∂ ∂y ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 - ∂ ∂x ∂ 2 v ∂x 2 + ∂ 2 v ∂y 2 - Ra Pr ∂θ ∂x .
Since we have L r ≪ δH (consequently ∂ 2 v/∂y 2 ≪ ∂ 2 v/∂x 2 ), three terms are dominant in Eq. (6.21):

∼ ∂ 2 v ∂x∂t inertia ; ∼ ∂ 3 v ∂x 3 f riction ; ∼ Ra Pr ∂θ ∂x buoyancy . (6.22)
In terms of characteristic scales, this reduced momentum balance reads Next, we turn attention back to the energy Eq. (6.20). During the pure conductive regime, the velocity is very small and negligible. The characteristic scale is therefore:

∼ v L r t inertia ; ∼ v L 3
δθ t ∼ δθ PrL 2 r , (6.26) leading to L r ∼ t Pr 1/2 . (6.27)
At this earlier stage, the thermal layer thickness is of the same order as L r , i. e. δ θ ∼ L r . When the convection arises, we can identify three distinct effects

1 t Inertia ; v δH Convection ; 1 L 2 r Pr Conduction .
(6.28) When the melting evolves (i. e. time increases), the inertia decreases in importance. There comes an effective time t e when the energy equation expresses a balance between the convection and the conduction heat transfer (see [START_REF] Bejan | Convection heat transfer[END_REF]): At this effective time t e , during which the inertia becomes negligible and the convection effects increase, the thermal layer thickness is

v δH ∼ 1 δ 2 θ Pr . ( 6
δ θ ∼ t e Pr 1/2 ∼ δH Ra e 1/4 ∼ δH -1/2 Ra -1/4 . (6.32)
Accordingly, before the flow oscillates, the scaling for N u could be written as:

N u ∼ δ1/2 H × Ra 1/4 . (6.33)
Equation (6.33) highlights that the heat transfer through a single convection cell depends on the height of the melted PCM layer during the linear regime. Correlations pertaining to basal melting case in the literature do not distinguish between this contribution of the bulk and the boundary layer heat transfers after the onset of convection and relies directly on empirical correlations for turbulent Rayleigh-Bénard flows. It should be noted that Eq. (6.33) rely on only one cell of the system.

Ra linear regime oscillating regime 3.27 • 10 5 0.286299 -1.635 • 10 6 0.274789 0.269158 3.27 • 10 6 0.279616 0.283658 6.54 • 10 6 0.274043 0.294528 Table 6.1: Exponent of the power laws of the linear and the oscillating regimes.

The exponents of the power laws obtained by our numerical simulations for four Rayleigh numbers are given in Tab. 6.1 during the linear and the oscillating regime. An exponent of 0.28 is observed for the linear regime and decreases with increasing Rayleigh number, which overestimates slightly the predicted scaling exponent of 1/4. This is mostly due to the number of convection cells in the melted PCM and the effect of the non-planar solid-liquid interface, not taken into account in our analysis. However, our numerical results match better with the Grossmann-Lohse theory [START_REF] Grossmann | Scaling in thermal convection: a unifying theory[END_REF], which predicts a scaling exponent of 2/7, within the framework of natural convection flow. Furthermore, during the oscillating regime the exponent increases with the Rayleigh number and tends to the empirical value of 1/3.

Concerning the temporal evolution of the liquid fraction L f , one can introduce the following scaling, from [START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF], obtained by a balance between the total rate of change of the enthalpy in the system and the heat fluxes entering and leaving the domain:

L f (t) = √ 2t (2-3β) + cRa β t 1/(2-3β) , (6.34)
with β the exponent in the N u -Ra power law, c = (2-3β)γ 2+Ste and γ a constant fitted from numerical data. For the conductive regime, by taking γ = β = 0, the well-known behavior L f ∼ t 1/2 is obtained. Replacing β by the exponent value 1/4 for the convective regime, we obtain L f ∼ t 4/5 . The time evolution of the liquid fraction is illustrated in Fig. 6.12a. As expected, the purely diffusive stage displays the scaling L f ∼ t 1/2 . Our simulations exhibit a power-law evolution of L f ∼ t 0.82 when the convection is fully developed in the fluid, in good agreement with L f ∼ t 4/5 . Figure 6.12b shows the temporal evolution of the heat flux, represented by the Nusselt number, and the strength of buoyancy, represented by the Rayleigh number, when the melting evolves. The onset of the convection arises around Ra e = 2 × 10 3 , when a sudden jump in the evolution of N u e is observed. [START_REF] Vasil | Dynamic bifurcations and pattern formation in melting-boundary convection[END_REF] have investigated a weakly non-linear stability analysis and have highlighted a super-exponentional amplitude growth when the Rayleigh number becomes close to the traditional critical value in the limit of vanishing Stefan number. This super-exponential growth is moreover followed by a rapid pattern readjustment. The trend of N u e at the onset of convection is in total agreement with the prediction of [START_REF] Vasil | Dynamic bifurcations and pattern formation in melting-boundary convection[END_REF]. The results of [START_REF] Esfahani | Basal melting driven by turbulent thermal convection[END_REF]; [START_REF] Madruga | Dynamic of plumes and scaling during the melting of a Phase Change Material heated from below[END_REF]; [START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF] exhibit the same trend despite the different boundary conditions (periodic lateral boundary conditions, adiabatic boundary conditions at the top of the cavity, and a low value of Pr for [START_REF] Esfahani | Basal melting driven by turbulent thermal convection[END_REF] and [START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF]). The rapid growth of N u e is followed by a power law with averaged exponent N u ∼ Ra 0.28 , for 10 3 ≤ Ra ≤ 10 5 . Finally, the transition from steady pattern of the convective rolls to oscillating pattern, followed by cell merging, as it is clearly shown in Fig. 6.10(d-f), is also illustrated by a decrease of N u at Ra ∼ 10 5 , followed by high oscillations in the temporal evolution of the heat transfer. The power-law relation N u -Ra is bounded in this stage by the average exponent of 2/7.

Comparison between LM and BM cases

A comparison of the liquid fraction L f , the Nusselt number N u, and the accumulated heat input Q 0 is offered in Fig. 6.13. Blue lines denote the LM case and red lines the BM case. Different values of the Ra number, ranging from Ra = 3.27 • 10 5 to Ra = 6.54 • 10 6 are investigated.

The temporal evolution of L f (left panels) displays a similar trend of both curves up to L f = 0.6 for all cases. For instance, for Ra = 3.27 • 10 5 (panel a), LM and BM cases evolve at the same rate from t = 0 to t = 100, corresponding to 60% of melted PCM. Then, for greater Rayleigh numbers, the gap keeps increasing, albeit slightly. This result agrees with the scaling analysis that predicts the same trend of both systems as long as the number of convective cells is close to one. Since the number of rolls is increasing with the Rayleigh number, an enhancement of the melting rate is thus consistent. Furthermore, differences are noticeable for L f ≥ 0.6, when the interface of the LM case reaches the right (cold) interface. Increasing the Rayleigh number induces a faster advancement of the top of the LM front, which touches the cold wall earlier, slowing hence the evolution of the melting. One can observe during this period a decrease of the slope of the blue curve with respect to the vertical axis. Linear evolution of L f is however noticed throughout the basal melting time since the interface touches the top cold wall very late.

Regarding the time evolution of the heat transfer (middle panels), red and blue curves match perfectly well during the conductive regime (see for e. g. Fig. 6.13 from t = 0 to t = 18 for Ra = 3.27 • 10 5 ). This was expected since the contribution of the convective heat transfer is still negligible at this stage. Differences occur from the emergence of the natural convection flow. The evolution of N u is relatively smooth for the LM case, while the BM case exhibits a sudden increase due to the hydrodynamical instabilities. Significant gap between both curves arises during the oscillating regime: N u displays a highly oscillating value for the BM case when the liquid layer expands upward, while the LM case reaches an asymptotic value. Note that for higher and higher Rayleigh numbers, the Nusselt number increases accordingly.

The evolution of the accumulated heat input Q 0 is consistent with the foregoing remarks. Q 0 is defined as follows:

Q 0 (t) = t 0 Nu(τ ) dτ, (6.35)
The decrease of the heat transfer monitored by the cell merging process in the BM case leads to a lower accumulated heat input. It follows that a more important thermal energy is stored when compared with the LM case.

These observations are of particular importance for PCM design for a specific application. If a fast melting is sought, the BM case is advised. It could be the case for example during quick temperature peak in electronic components. Conversely, when the PCM is expected to melt in a longer period, the LM case is more appropriate.

Concluding remarks

The high accuracy of our numerical method permitted to simulate in this chapter lateral and basal melting cases. Numerical and theoretical comparisons of both cases were performed and have exhibited many differences, mainly on the dynamics of the flow and the heat transfer processes. A main question related to convective melting processes is to predict the temporal evolution of the melting rate, which is correlated to the heat-flux dynamics determined by the flow in the system. Since the dynamics of the flow is completely different for LM and BM systems, we have performed simulations for both cases and compared the evolution of some physical parameters.

We have investigated first the LM case, by simulating the melting of n-octadecane inside a square enclosure heated from the vertical wall. Three principal regimes were identified to describe the dynamics of the melting: conductive regime, mixed regime, and convective regime. The influence of the Rayleigh and the Stefan parameters on the flow was also assessed and indicated clearly that increasing Ra, by increasing either H or δT , enhances the heat transfer. Since stronger convection flow occurs in the liquid region with increasing Rayleigh number, the shape of the phase-change front is altered consequently. The top of the interface moves faster while the bottom part is rather slowed down, resulting in a curved shape of the melting front.

Concerning the simulations of the BM case, simulations for a large range of Rayleigh numbers, Ra = 3.27 • 10 5 to Ra = 6.54 • 10 6 , were carried out for a comprehensive description of the dynamics of the melting. Effective Rayleigh and Nusselt numbers, depending on the fluid height were introduced to describe the flow. The melted fluid layer was shown to be thermally unstable. It quickly develops convective motion of progressively higher intensity as the depth of the melted layer increases. The onset of these instabilities were observable either on the dynamics of the melting flow or the temporal evolution of the Nusselt number. The novelty is the theoretical description of the convective heat transfer during the linear regime, i. e. Ra e ≤ 10 5 , through a scale analysis by separating the contribution of the bulk and the boundary layer heat transfers. For higher Ra e numbers, we used the knowledge acquired on turbulent natural convection system to understand convective melting.

We have finally compared the time evolution of the liquid fraction L f , the heat transfer rate represented by the Nu number, and the accumulated heat input Q 0 . We showed that PCM melts faster when heated from below compared to the lateral melting, since the last case reaches the cold wall earlier, slowing accordingly the melting rate. However, the lateral oscillation of the convective cells in the BM case, yielding to cell merging processes, reduces the heat transfer rate. In the following chapter, we investigate the full melting-solidification cycle for a better understanding of the complete PCM cycle.

Numerical simulation of complete melting-solidification cycle of a phase-change material

This chapter presents an accurate numerical description of the alternate melting and solidification of a PCM. Square and cylindrical PCMs, that were simulated in Chapters 5 and 6 are considered in this chapter to investigate the solidification process. The fundamental operational mode of latent thermal energy storage (LTES) systems based on phase-change materials (PCM) is made of alternate melting and solidification cycles that are not necessarily periodic.

The natural convection occurring in the melting PCM is driven by the temperature difference δT = T h -T f . The dimensionless number that depicts the ratio between the forces creating and those refraining the motion is the Rayleigh number, which appears in the dimensionless form of the Navier-Stokes equations with Boussinesq approximation (Sec. 2, Eq. (2.16)). The higher is its value, the more intense is the heat transfer. Conversely, during the solidification, the phase-change is handled by the discharged temperature T co < T f , where the subscript 'co' stands for 'cooling'. For the square geometry discussed first in Sec. 7.1, this represents the temperature of the left wall, and for the cylindrical geometry presented in Sec. 7.2, it consists of the inner tube temperature. Thus, the relevant temperature difference in the solid phase of the PCM is δT co = T f -T co and the dimensionless temperature in the solid phase should be defined with respect to this δT co . It is then obvious, from Eqs. (2.15) and (2.16), that the Rayleigh number should be defined using the same temperature difference. However, because the Rayleigh number amounts for the motion created by the hot temperature difference, we choose to keep the same definition for the Rayleigh number as for the melting case. This is still relevant for the melted core of the flow, where the persisting motion acts as a boundary condition for the solidification process. Under these conditions, in regard with the solidification process, we introduce a new parameter, r δ = δT co /(T h -T f ), the normalised temperature with respect to T f -T co and the relevant Rayleigh number will be then Ra co = r δ × Ra, where Ra co is the pseudo-Rayleigh number for solidification with a melted boundary. In the following, we will describe the process of solidification using three different values of r δ . A new scaling is moreover introduced:

V ref = α l H ⇒ t = t ϕ ν l H 2 Pr ⇒ Re = 1 Pr . (7.1)
The solidification stage is indeed a slower process compared to the melting, therefore the use of an adapted scaling is more relevant. This leads to a different time scaling for each part of the cycle.

In Sec. 7.1, we perform two operating cases for the solidification process of melted PCM in a differentially heated square cavity. In the first study case the solidification starts after a complete melting of the PCM (liquid fraction of 95%), while in the second case after a partial melting (liquid fraction of 50%). Several practical implications for the two operating modes are drawn. In Sec. 7.2, we investigate the solidification of a cylindrical PCM with an inner obstacle. We provide temporal evolution of solid-liquid interface, liquid fraction, Nusselt number and accumulated heat input.

Melting-solidification cycle of a phase-change material with complete or partial melting in a square cavity

The sketch of the computational domain and boundary conditions are illustrated in Fig. 7.1a, corresponding to the melting of octadecane PCM presented in Sec. 5.2 of Chapter 6. Starting from a melted PCM (Fig. 7.1b), the simulation of the solidification process starts by imposing at the left-wall a constant (cold) temperature θ co as it is shown in Fig. 7.1c. We consider two cases:

-case CM: solidification after a Complete Melting of the material (L f = 0.95, Fig. 6.2f) and -case PM: solidification after a Partial Melting (L f = 0.5, Fig. 6.2d). The solid phase will propagate into the cavity from both left and right sides (Fig. 7.1c), which makes this case computationally challenging. The mesh adaptivity capabilities of our numerical code made possible to accurately track the two solidification fronts identified by the iso-line θ = 0. In the discussion below, the results will be presented using the physical time t ϕ for the same reason mentioned above (note the influence of the size of the domain, since the scaling in Eq. ( 7.1) uses the height of the cavity).

The solidification process starts at physical time t ϕ = 185 min (corresponding to t = 250) for the case CM and at t ϕ = 59 min (t = 78.8) for the case PM.

Solidification after a complete melting (case CM)

The solidification follows after a complete melting corresponding to Fig. 6.2 at t ϕ = 185 min (t = 250). The dynamics of the PCM during the solidification process is reported in Fig. 7.2, with the evolution of temperature distribution and solidification fronts.

At t ϕ = 185 min (Fig. 7.2a), the liquid fraction is L f = 0.95 and the melting/solidification front is close to the right wall of the cavity. Setting smoothly a low temperature θ co = -1 at the left wall, while the right wall is maintained at a constant temperature (θ right = -0.01 ≤ θ f ), triggers the formation of a second solidification front, propagating from the left side of the domain. Figures 7.2b and 7.2c illustrate that the left part of the cavity is solidifying at a faster rate because of the very low temperature imposed at the left wall, inducing a non symmetric evolution of solid-liquid interfaces. The solid part is represented in blue and corresponds to the region of temperature θ ≤ 0. The signature of the conductive heat transfer is characterized by the vertical shape of the left front. Inside the liquid, the initial convection cells facilitate the heat transfer from the boundaries, resulting in a very rapid decrease of the fluid temperature. Temperature gradients being smoothed out during this first stage, the influence of the convection inside the liquid region is considerably reduced. As a result, the velocity inside the liquid is reduced to very low values.

Starting with t ϕ = 430 min (Fig. 7.2d), the shape of both interfaces is almost symmetrical. This is a signature of a conduction dominated process. At t ϕ = 510 min (Fig. 7.2e) the liquid region starts to shrink at the bottom side of the cavity. This process is accelerated and finally the liquid is trapped in a thin pocket and finally disappears completely through the top of the cavity (Fig. 7.2e). The complete solidification ends at t ϕ = 530 min, when the liquid fraction is L f = 0.

The adapted mesh, refined along the two solidification fronts, at t ϕ = 300 min is reported in Fig. 7.2f, illustrating the efficiency of the adaptive mesh tool.

Solidification after a partial melting (case PM)

In this case, the solidification starts from the state corresponding to Fig. 7.3a at t ϕ = 59 min (t = 78.8), when the liquid fraction is L f = 0.5. The temperature of the left wall is progressively lowered at θ co = -1, as in the previous solidification simulation.

The time evolution of the process is illustrated in Figs. 7.3a-e, while the adapted mesh corresponding to t ϕ = 90 min is plotted in Fig. 7.3f. As in the previous case, a second solidification front starts to propagate from the left side of the cavity. The straight shape of the left solid front is always observed while the right solid front is impacted by the convection cell present in the central liquid region (Fig. 7.3b). The stronger convective effect is most likely due to the huge temperature difference that occurs over a smaller space distance (almost half of the volume is occupied by the solid state). This leads to stronger temperature gradients in the liquid region, and consequently to a stronger heat transfer. The two fronts merge to form a pocket of fluid which is connected to the top of the cavity (Fig. 7.3c-e).

It is interesting to note that, as in the previous solidification case, the left part is solidifying at a faster rate, hence the pocket of melted PCM disappears completely from the right at the top side of the cavity (Figs. 7.3c-e).

Analysis of the solidification cycle from two different initial conditions: complete (case CM) and partial melting (Case PM)

The aim of this section is to investigate the temporal evolution of some physical properties of the solidification process, from two different initial conditions: i) completely melted volume (case CM) and ii) partially melted volume (50% of the fluid is melted, case PM). Figure 7.4 represents the temporal evolution of the liquid fraction L f , the Nusselt number Nu, and the accumulated heat input Q 0 (Eq. (6.35)), for the two investigated cases. Note that the average Nusselt number is calculated at the cooled wall, defined similarly to Eq. (2.24), but it can be negative in this case. Simulations for three values r δ = 1, r δ = 5 and r δ = 10 were carried out.

Figure 7.4a illustrates the temporal evolution of the liquid fraction L f for the CM case. Complete melting occurs for t ϕ = 185 min, after which solidification starts, with a continuous decrease of L f untill complete solidification is achieved. For the lowest value of r δ , corresponding to Ra co = 3.27 • 10 5 , the solidification process ends at t ϕ = 530 min. Then, the higher the value of r δ , the faster the discharge process. Final times for solidification are t ϕ = 260 min and t ϕ = 230 min for cases Ra co = 1.62 • 10 6 and Ra co = 3.27 • 10 6 . These values of Ra co correspond to a drop of the cold boundary temperature by a factor of 5 and 10, respectively. The solidification speed, quantified by dL f /dt ϕ is nearly constant during almost the whole process for each case. This uniformity of the process indicates that the natural convection flow vanishes rapidly during the solidification, and conduction remains the only heat transfer mode. Figure 7.4b plots the temporal evolution of L f for the PM case. As previously discussed, 50% of the volume is melted, at time t ϕ = 59 min, when solidification starts. Furthermore, despite that solidification process is started, L f continues to increase slightly at the very beginning of the discharge stage. Then it decreases monotonically towards 0 at t ϕ = 240 min. The heat stored in the melted PCM continues to melt the remaining solid PCM until the convection becomes negligible. It is worth noticing that this behavior is not observed in the complete melting case because of the imposed temperature at the right wall.

Let us now pay attention to the transfers occurring at the left wall, progressively submitted to a lower temperature. This is done through the temporal evolution of the Nusselt number and the accumulated heat input. Panels (c) and (d) of Fig. 7.4 illustrate the Nusselt number for the CM and PM cases. The three investigated Rayleigh numbers are shown, with clear differences between them. This difference corroborates with that already reported for the melting case, over shorter times scales. This indicates that the heat transfer during the solidification process is fundamentally different from the melting one. For the CM case, for Ra co = 3.27 • 10 5 , the Nusselt number first decreases sharply, for t ϕ ≤ 18 min, then it reaches a plateau at Nu = 7 during the complete melting. At t ϕ = 185 min, solidification starts and Nu suddenly decreases over very short times, reaching negative values (Nu ≈ -15). It follows an increase of Nu with time, up to reaching an asymptotic value close to 0 (zero temperature gradients, i. e. uniform temperature at the left wall). The same mechanism is observed over a shorter time interval when Ra co is increased. For the PM case, the Nusselt number also decreases sharply to a negative value when the solidification starts. However, the convection flow remaining in the melted region influences the heat transfer at the very beginning of the solidification process. The hot fluid in the middle of the melted PCM is advected by the natural convection flow to the boundaries and induces a temperature gradient at the left wall, resulting into an oscilating behavior of the Nusselt number before reaching an asymptotic value. This is in agreement with the previous comment about the melting continuing in the right part of the cavity, despite that solidification has started, and the slight increase of the liquid fraction at the very first time steps of the discharging process. Q 0 defined in Eq. (6.35), as it is shown in panels (e) and (f) of Fig. 7.4. Heat is first stored during the melting stage, corresponding to t ϕ ≤ 185 min for CM (Fig. 7.4d) and t ϕ ≤ 59 min for PM (Fig. 7.4e). It is then restored during the solidification stage. The CM case indicates higher value of Q 0 (Q 0 = 1400, for Ra co = 3.27 • 10 5 ) compared to the PM case (Q 0 = 500), meaning that the PCM is more efficient in terms of heat storage. However, the PM case exhibits well balanced characteristic times between the solidification and the melting stages for Ra co = 3.27 • 10 5 . Besides, when the Ra number increases, the stored heat is discharged faster. Moreover, the temperature and the velocity profiles drop sharply during the first step of the cooling process and become almost equal to zero very early in the whole domain. This means that conduction dominates the solidification process, and the convection becomes rapidly negligible. As a consequence, the melting fronts are vertical and have a symmetric position with respect to the center of the cavity.

Melting-solidification cycle of a cylindrical PCM

We are focusing in this section on the solidification of cylindrical PCM of radius R = 1 with an inner cooled tube. This configuration is used for heat storage systems. The simulation of the solidification process starts by imposing progressively at the inner tube a constant cooling temperature θ co = -1.

The outer tube is adiabatic and a no-slip boundary condition is applied for the velocity. The solidification part of the cycle occurs at t = 2.5 when 80% of the initial solid PCM is melted. A solid layer grows around the tube and joins the remaining solid PCM.

Figure 7.5 shows the time evolution of the solidification process for r δ = 1. The simulation resumes from the melting state presented in Fig. 5.10. The low temperature at the inner tube triggers a second solid-liquid interface, propagating from the center to the outer tube (see panels a and b). Both solidification fronts merge at t = 4.34 (panel c), followed by a shrinking process of the liquid (panel d and e) which disappears completely at t = 6.7 (panel d). We note that, as emphasised in Sec. 7.1, the heat transfer is dominated by conductive transfer during the solidification cycle, explaining the circular evolution of the solidification front. The adapted mesh corresponding to t = 3.16 is shown in panel (f), illustrating the refined mesh along the two interfaces θ = 0.

From the scaling analysis (6.4), the radial expansion δ r of the solidification front from the center could be approximated by

δ r = Ste Pr × t. (7.2)
The solidification time could be thus predicted by

t sol = δ 2 r Pr Ste . (7.3)
In our case, by taking into account the remaining solid located at δ r = 0.8, the theoretical solidification time is t sol = 6.4, which is in relatively good agreement with our numerical results. Increasing the Stefan number should hence result in a lower value of the solidification time.

Concluding remarks

The choice of simulating the complete cycle melting-solidification of n-octadecane PCM in this chapter is motivated by the fact that this material is generally used for buildings applications, due to its phase change temperature of 28 o C. [START_REF] Zhu | Dynamic characteristics and energy performance of buildings using phase change materials: a review[END_REF] and [START_REF] Kalnaes | Phase change materials and products for building applications: a state-of-theart review and future research opportunities[END_REF] listed various applications, starting from free cooling, peak load shifting, passive building systems and solar energy storage. For each case, the PCM is assumed to melt during day-time and to solidify during the night-time.

It was noticed that, when the same Rayleigh number is considered to compare the melting and solidification cycles, meaning that the charge and the discharge modes occur at the same δT , the solidification is always slower than the melting. This behavior is linked to the heat transfer mode leading each part of the cycle. Convective heat transfer dominates indeed the melting process, enhancing thus the heat transfer, while conduction is the main heat transfer mode during solidification, resulting into a slower operating process. However, when the discharge temperature is decreased by a factor of 5, i.e both Ra co and r δ × Ste are increased, the solidification and the melting occur over similar times.

A first issue that has been brought up by [START_REF] Ascione | Energy refurbishment of existing buildings through the use of phase change materials: Energy savings and indoor comfort in the cooling season[END_REF] is the difficulty of the PCM systems to completely discharge during night-time. Though, if the PCM does not solidify entirely, the effectiveness of the system may be considerably reduced. In this case, to have a shorter cooling period, it is not advised to melt the PCM completely.

However, for solar energy storage applications, full melting of the PCM is needed to utilize its latent heat storage capacity. Hence, a partial melting is not optimal. For other applications, when a shorter discharge time is desired, the use of external cooling techniques is needed to ensure a colder discharge temperature. We present in this chapter parallel computations of three-dimensional liquid-solid phase-change systems involving natural convection. We use the recent library ffddm that makes available in FreeFem++ a state-of-the-art scalable Schwarz domain decomposition method (DDM). Our motivation to expand our numerical model to 3D configurations is drawn by the lack of publications in the literature presenting accurate 3D simulations of phase-change materials. Also, experimental investigations against which we have validated our numerical method in previous chapters involve three-dimensional effects that we have assumed to be negligible in our comparisons. The latter assumption is however not valid for high Ra numbers.

The main feature of our numerical approach is the use of 3D adaptive mesh, performed using mmg3d library. Mmg3d is a 3D software developed by [START_REF] Dobrzynski | MMG3D: User Guide[END_REF], which allows to remesh an initial mesh made of tetrahedra. The metrics for the mesh adaptation are computed using mshmet library, which provides anisotropic metrics based on solution variations. Since the mesh adaptivity procedure is based on a sequential algorithm, we compute first the metrics with respect to the coarse mesh. Then, finer meshes are generated before the mesh decomposition step in order to reach the desired level of refinement for the subdomains. We use Metis library to split the domain into subdomains. The number of layers of mesh elements in the overlap region between subdomains is set to 2 for all subsequent simulations. For three-dimensional applications, direct solvers used in the framework of 2D problems are not appropriate and iterative methods must be employed since memory requirements for the LU decomposition would rapidly exceed the capacity of available computers. The linear system of equations resulting from the Newton linearization are thus solved using a parallel GMRES Krylov method. Since it is well known that iterative methods can suffer from convergence problems, we adapt the number of subdomain in such a way that each processor could handle 1, 000 tetrahedra. The Optimized Restricted Additive Schwarz (ORAS) preconditioned GMRES solver proved to be extremely efficient since around 30 iterations were necessary to get a residual of 10 -9 for each Newton iteration.

The remainder of this chapter is as follows. In Sec. 8.1, we validate the 3D convection of air in a cubic cavity against the numerical benchmark by [START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF], and compare together the solutions obtained using the parallel and the sequential algorithms. The convection of water is investigated in Sec. 8.2. Then, the melting of n-octadecane PCM inside 3D enclosures is addressed in Sec. 8.3. We finally present in Sec. 8.4 the challenging numerical simulation of water freezing.

Numerical simulation of the natural convection of air in a cubic cavity

Following the same argument of increasing gradually the complexity of computed cases, that was outlined in Chapters 4 and 5, we start by presenting the natural convection of air. This case involves a linear buoyancy force f B . The physical parameters are the same that used in the 2D simulations (Pr = 0.71) and we investigate Rayleigh numbers ranging from Ra = 10 4 to Ra = 10 6 . The walls are rigid and impermeable. The vertical walls at x = 0 and x = 1 are isothermal and have different temperatures θ h = 1 and θ c = 0, respectively. The other walls are considered adiabatic. The fluid is initially at rest and the temperature is linearly distributed from the cold to the hot walls. We solve the steady Eq. (3.14) by increasing smoothly the parameter α in front of the Boussinesq term (which can be assimilated to a Rayleigh continuation step) with a maximum of 6 steps for Ra = 10 6 .

We first compare the current simulation with the numerical data of [START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF]. These authors used a forth order finite difference method, with a vorticity-stream function formulation with different uniform meshes of 120 × 120 × 120 grid nodes. Our results were obtained using uniform grids of 40 × 40 × 40. Since the converged flow pattern and temperature distributions are symmetrical with respect to the center of the cavity for the investigated Rayleigh number, we display in Fig. 8.1 the temperature field at the mid section (y = 0.5), for each of the three Rayleigh numbers Ra = 10 4 (top), Ra = 10 5 (middle), Ra = 10 6 (bottom). On the left we display the numerical results of [START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF] and on the right the present simulation. The comparison with the benchmark solution exhibits a good qualitative agreement. The higher the Rayleigh number, the finer the thermal boundary layer thickness is in the vicinity of the vertical walls. One can also notice the stagnant fluid with stratified temperature in the center of the domain in both numerical solutions.

A more accurate comparison is presented in Tab. 8.1 between the present simulation and numerical data from the literature. and w max of the velocity profiles at mid-domain (x = 0.5 and z = 0.5, respectively) and their locations z max and x max . Comparison with benchmark solutions of [START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF] and numerical data of [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF] for Ra = 10 4 , Ra = 10 5 , Ra = 10 6 .

We compare the maximum velocity values and their corresponding location with benchmark solution of [START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF] and numerical results of [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF]. Our results show a good convergence on a mesh 3 times coarser than the one proposed by [START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF]. A good agreement is found for all values of the Rayleigh number. For Ra = 10 4 the current simulation shows a relative difference of 0.15% for u max and 0.28% for v max . Larger differences of 6.55% is found with the results reported by [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF]. We recall that Moglan (2013) used a second-order finite difference scheme on a 64×64×64 grid. For Ra = 10 5 and 10 6 the results are in good agreement with both references.

Further comparison of the parallel solver with the sequential algorithm is offered in Tab. 8.2 for all cases. We compute L 2 and L ∞ norms for the velocity and the temperature. The number of subdomains varies from 28 to 70 for 1.8 millions of unknowns. The difference between the two algorithms is of order of 10 -6 and we do not observe a large variation of the error when the number of subdomains is increased. The comparison with the sequential algorithm was limited to the 40 × 40 × 40 grid, since the simulations are highly demanding in CPU time. For the highest value of the Rayleigh number, the steady state computation required 57 CPU hours and 3 runs (restarts) with 120 Go of memory for the sequential algorithm. The computational time is considerably reduced using DDM in 3D simulations. This becomes affordable since only 21 CPU minutes were necessary with 70 processors to compute the same case with an error of 6.02152 • 10 -6 on u and 4.

Numerical simulation of the natural convection of water in a cubic cavity

In this section, we simulate the natural convection of water inside a cubic cavity using adaptive 3D meshes. The dimensionless parameters are the same that used in Chapter 4.2 (Ra = 2.518084 • 10 6 and Pr = 6.99). We impose cold dimensionless temperature θ c = 0 at x = 1 (right wall), hot temperature θ h = 1 at x = 0 (left wall), and a homogeneous Neumann boundary condition at the other walls. Dirichlet boundary condition u = 0 is prescribed over the whole boundary ∂Ω. In view of anomalous thermal variation of water density, we adapt the mesh along θ = 0.4 to capture correctly the flow structure. The main limitation of 3D simulations of water convection in the literature is the size of the mesh resolution when fixed grid models are used. As an example, Giangi et al. (2000) and [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF] presented a mesh sensitivity analysis, and concluded that even by using 81 3 grid points, the variation of the velocity was still significant, making the problem computationally expensive.

The temperature distribution and the corresponding adapted mesh for the steady state computation is shown in Fig. 8.2a. The blue region denotes the cold water trapped by the abnormal fluid recirculation and the red region the hot fluid driven by the upper clockwise circulation. A zoom of the mesh at the mid-plane is also offered in Fig. 8.2b. Smaller tetrahedra are clearly observed in the vicinity of the walls, where a spiral movement of the fluid occurs along the walls, and between the two counterrotating circulation patterns. A minimum edge length h min = 3.33 • 10 -3 is used along the maximum density variation, and h max = 0.15 in the stagnant fluid region. The combined mesh adaptation strategy and the efficient parallel algorithm, allow us to reduce considerably the computational time. Only 5012.04 CPU seconds were necessary to compute the steady solution with 500, 000 tetrahedra using 56 processors, while [START_REF] Michalek | Natural convection for anomalous density variation of water: numerical benchmark[END_REF] needed 3.6×10 5 CPU seconds to compute the steady solution with a 81 3 fixed grid.

In many cases, it is convenient to perform two-dimensional numerical simulations, to reduce the computational effort and allow many simulations to be achieved in a realistic time. However, Fig. 8.2b shows that the three-dimensionality of the flow in a cube shape cavity affects the topology of the isosurface θ = 0.4. Giangi et al. (2000) assessed on the effect of three-dimensionality in both convection and freezing of water and they noted that only the flow in the symmetry plane can be considered as two-dimensional. The no-slip velocity and the adiabatic thermal boundary conditions at the side walls enhance the three-dimensional effect near the walls. Furthermore, these authors noted that, even though the temperature profiles at the central cross-section of the cavity match well for both 2D and 3D simulations, important differences between two and three-dimensional velocity profiles are present.

In Figs. 8.2d-f, we plot the velocity profiles along the vertical line passing through the velocity saddle point, where normal and abnormal convection streams collide in the vicinity of the cold wall.

For the e z component of the velocity, the profile along the z-direction at x = 0.93, for both two and three-dimensional simulations is compared with the 3D simulation of [START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF]. Triangular symbols denote the reference solution, obtained with a 3D finite difference code using vorticity-vector potential formulation of the Navier-Stokes and energy equations. 3D (red solid line) profile agrees well with the benchmark solution with a maximum difference of 3%. The 2D (green dashed line) solution exhibits higher discrepancies in the vicinity of the anomalous density variation, where the three-dimensional effect is maximum. However, as noted by Giangi et al. (2000) in their experimental observations, larger differences between 2D and 3D solutions could be observed at x = 0.1 and x = 0.5. Panels (e) and (f) exhibit discrepancies greater than 10% between both solutions, in agreement with the foregoing observation.

Numerical simulations of melting PCM with convection in a 3D configuration

We now proceed to the 3D simulation of solid-liquid phase-change problems, starting with the melting of a PCM in a cube-shape domain. The three-dimensionality of the fluid flow can lead to major differences with the two-dimensional solution in some configurations, as in the example of the melting of Gallium. Nikrityuk (2012) showed that the controversial multi-cellular flow occurs only at the beginning of the melting process in two-dimensional simulations. Three-dimensional simulations exhibit the absence of a multi-cellular structure at any time step. They explained these differences by the presence of walls in the third direction, which suppresses the flow. Also, the occurrence of weak turbulence in the bulk of the melt prevents the formation of the long-living multicellular flow pattern known from 2D simulations. The three-dimensional geometry allows the fluid movement in all directions and the friction effect provided by no-slip walls in the transverse direction influences the flow. Nikrityuk (2012) highlighted the existence of a secondary flow in the transverse direction, as shown in Fig. 8.3, which has its maximal intensity near the side walls.

We expand the 2D simulations presented in Chapter 5 to 3D configurations. Solutions obtained by both configurations are compared to assess on the three-dimensional effects that occur during the melting process. The configuration of the numerical simulation is sketched in Fig. 8 The location of the solid-liquid interface at the mid-plane and the time evolution of the Nusselt number for two and three-dimensional configurations are plotted in Fig. 8.5. The position of the interface in panel (a) displays discrepancies between both melting fronts, but the shape remains in good agreement. These differences are caused by the increase of the velocity in the liquid flow, induced by the three-dimensional effects, as shown by Nikrityuk (2012), who compared the volume-averaged velocity U xyz , defined as follows:

U xyz = Ω u 2 1 + u 2 2 + u 2 3 dV. (8.1)
In our simulations, at t = 78.7, we obtained U xyz = 0.548208 for 3D simulations and U xy = 0.519576 for 2D. However, the time evolution of the heat transfer in panel (b) matches well for 2D and 3D simulations, with a maximum difference lower than 1% at the onset of the convection time, and overlap during all the simulation time. From an engineering point of view, the latter observation means that if only an assessment of the heat transfer during the phase-change process is sought, a 2D simulation is enough since it could be performed in a reasonable computational time. The 3D simulation required 2 days of CPU time with 200 processors, while 45 minutes on a personal computer is enough for the 2D configuration, for the investigated Rayleigh number.

Numerical simulation of water freezing in a cubic cavity

In this section, we present a 3D simulation of the challenging case of ice formation inside a differentially heated cubic cavity. Starting from the steady solution shown in Fig. 8.2, we drop the cold temperature at the right wall from θ c = 0 to θ c = -1, under the temperature of fusion (which correspond to physical temperatures T = 0 o C and T = -10 o C, respectively). Solid crust arises thus from the cold wall and expands towards the left wall.

The temperature distribution and the adapted mesh at final time t ϕ = 2340 s (t = 1.61) are displayed in Figs. 8.6a and 8.6b. When compared with the natural convection case in Sec. 8.1, the mesh is additionally adapted along the solid moving front, localized by the temperature iso-surface θ = 0. For the visual identification of flow structures during the ice formation, panel (c) plots the snapshots of velocity vectors, the iso-surface θ = 0.4 along the anomalous density variation, and the phase-change front at time instant t = 1.61. The effect of the secondary flow, discussed in Sec. 8.3, is also visible on the curved shape of the iso-surface θ = θ m in the transverse ydirection. One can note however, that the solidification front is almost 2D. The buoyancy-induced fluid motion in the abnormal recirculation region is too weak to influence the solid front. Finally, the overlay of the experimental image of [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF] and the current simulation exhibits good agreement for the location of the solid-liquid interface. Differences come mainly from the fact that the undercooling phenomenon during the solidification stage is not taken into account in our model. However, the global shape of the interface agrees qualitatively with the experimental observation.

Concluding remarks

We developed parallel 3D tools using domain decomposition methods to simulate natural convection flows involving melting or solidification boundary conditions. The numerical simulations performed with the recent library ffddm in FreeFem++ showed good agreement with the benchmark solutions and the sequential algorithm for all simulated cases. The qualitative behavior and the topology of the interface for phase-change problems were discussed in detail. Three-dimensional effects in the flow were exhibited by the deformed shape of the interface in the transverse direction. The adaptive 3D mesh procedure proved very efficient to capture several interfaces, as in the example of the water freezing. Grid refinement effort was concentrated only in regions displaying high gradients of the computed variables, and 10 times larger tetrahedra were generated elsewhere.

Conclusion and perspectives

In this thesis, we presented numerical simulations of natural convection with melting boundary. We developed numerical tools for simulating 2D and 3D configurations of solid-liquid phase-change systems involving natural convection, using adaptive finite elements method.

The buoyancy-induced fluid motion in the liquid phase was modeled by the incompressible Navier-Stokes equations with Boussinesq approximation. The phase-change was modeled by an enthalpyporosity model using a temperature-based formulation for the energy conservation equation. A singledomain method, solving the same system of equations in both phases, including a Carman-Kozeny-type penalty term in the momentum equation, was implemented with the FreeFem++ software.

In the numerical approach, the coupled momentum and energy equations were integrated in time using a fully implicit second-order GEAR scheme. The space discretization was based on a Taylor-Hood triangular finite elements, approximating the velocity and the pressure with P 2 and P 1 finite elements, respectively. Since we noticed from the literature review that the published numerical simulations required considerable high computer resources and CPU time, the main feature of our numerical method was the use of a mesh adaptivity algorithm using metrics control to adapt the mesh every time step. We concentrated the grid refinement effort only in regions displaying high gradients of the computed variables, and prescribed coarser meshes in the solid phase and stagnant regions. This reduced considerably the number of degrees of liberty, while keeping high accuracy. The accuracy of the numerical method was tested using manufactured solutions: the expected second-order accuracy in both space and time was obtained.

As a first purpose of the study, we organized the programs as a toolbox for the software FreeFem++. The toolbox was designed to be an easy-to-use tool to handle phase-change problems. All technical issues related to the implementation of the finite element method were hidden, allowing to focus on numerical algorithms and their performance. A sequence of validation tests was performed, and presented by increasing gradually the level of complexity. Using the provided examples, the user could implement easily similar problems involving phase-change phenomena.

We first tested the capability of the Navier-Stokes solver to deal with natural convection problems, without enthalpy and porosity terms in the system of equations. We started by simulating the natural convection of air and water inside differentially heated square enclosures. We validated our results against classical benchmarks of natural convection. Rayleigh numbers up to 10 8 , at the limit of the steadiness for the natural convection of air, were considered. Additionally, we investigated the thermal driven problem including a heated obstacle, to assess for the robustness of the code.

We also computed the natural convection of water, which implies a non-linear variation of the density in the expression of the buoyancy force. The comparison with existing numerical data in the literature showed very good agreement. The mesh adaptivity algorithm proved very efficient, allowing to capture accurately the anomalous variation of the density around T = 4 o C. For all simulations presented in the framework of natural convection problems, the steady state regimes were reached at most after 50 CPU seconds of computation.

We validated extensively the code by simulating the problem of melting and solidification of a PCM.

To compare our simulation with available experimental and numerical results in the literature, we first considered the melting of n-octadecane and Gallium inside vertically heated rectangular enclosures. A mesh convergence test was performed to ensure the grid-convergence of the solution.

Three Rayleigh numbers ranging from 10 5 to 10 8 were considered for the n-octadecane PCM melting. Our simulations were compared with three different benchmarks. A good agreement was found for the comparison with existing experimental and numerical results. Concerning the controversial melting of Gallium in 2D configurations, our simulations showed multi-cellular structures of the melted flow and highly oscillating values of the Nusselt number, in agreement with many observations from the literature.

We also simulated complex geometries, such as highly distorted domain or cylindrical PCM including inner heated pipes. For all considered cases, a variable mesh in the vicinity of the walls and along the phase-change front ensured a good capture of the boundary layer region and the location of the interface. The simulation of the melting of Gallium showed that our simulations run 245 times faster when compared with the simulations of [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF].

We finally presented the challenging case of water freezing inside a differentially heated cavity. A qualitative agreement with the experimental results for the solidification front was obtained. The efficiency of our numerical algorithm allowed us to perform simulations for melting problems, within 45 minutes for the lowest value of Ra to 5 hours for Ra = 10 8 . After a thorough investigation, we could conclude that our numerical method was validated for the natural convection of air and water and for the solid-liquid phase change problems including melting and solidification by providing a good agreement with the benchmarks presented in the literature.

The second purpose of the study was to provide a thorough analysis of the melting and the solidification process and compare with numerical results.

We first compared the behavior of the melting PCM when the latter is subject to heating from the side or from below (lateral or basal melting). For the lateral melting case, we assessed on the influence of the Rayleigh number on the heat transfer, by computing several configurations using different values of H and δT . It was shown that increasing the Rayleigh number by keeping δT constant induces a slower melting rate and a higher heat transfer. However, increasing both Rayleigh and Stefan numbers induced an earlier onset of the convection-dominated regime, improving thus the efficiency of the PCM.

For the basal melting case, we developed a scale analysis for the linear regime and provided a comprehensive description of the heat transfer processes during the melting process. Differences between lateral and basal melting were drawn, mainly for the structure of the flow and the time evolution of the Nusselt number. A mono-cellular flow was observed for the lateral melting case at any investigated Rayleigh numbers, while natural convection developed in the form of Bénard cells for the basal melting case. Moreover, the quasi-steady evolution of Nu for the vertical heating is not observed when the PCM is heated from below, since high oscillating evolutions were exhibited.

Further investigation of the solidification process was done by simulating alternate melting and solidification cycles. Solidification after either complete or partial melting was simulated, with an assessment of the melting rate, the heat transfer and the accumulated heat input. It was observed that the convective heat transfer dominated the melting process, enhancing thus the heat transfer. Conduction was the main heat transfer mode during the solidification, resulting into a slower operating process. However, when the discharge temperature was decreased by a factor of 5, the solidification and the melting occur over similar time intervals. The challenging task for the mesh adaptivity procedure was to track the two moving interfaces during the solidification cycle.

The third purpose of the thesis was to expand the numerical simulation to 3D configurations. Parallel 3D tools for melting or solidification involving natural convection flows, using domain decomposition methods were presented. Numerical simulations were carried out using the recent library ffddm in FreeFem++ . The originality of our numerical approach was the use of 3D adaptive meshes, performed every time step, using mmg3d library. We used Metis library to split the domain into subdomains. The linear system of equations resulting from the Newton linearization was solved using parallel GMRES Krylov methods.

We simulated the natural convection of air and water inside cubic enclosures and found good agreement with classical benchmarks. We also simulated the melting of n-octadecane PCM and the difficult case of water freezing inside a differentially heated cubic cavity. The influence of the three-dimensional effects on the flow, that was neglected in 2D simulations, was highlighted by presenting the 3D shape of the solid-liquid interface. Iso-surfaces of the temperature were also impacted by the secondary flow in the vicinity of the side walls. The adaptive 3D mesh procedure proved very efficient to capture several interfaces, as in the case of the water freezing.

Future work

One of the major limitations of the present physical model is that we neglected the undercooling phenomenon, since we assumed that the temperature of fusion and the temperature of solidification were the same. Though, the liquid could solidify at a temperature significantly lower than the melting temperature, because of a problem of nucleation at the microscopic scale. An enthalpy-based formulation of the enthalpy-porosity model should be implemented to this end, since the undercooling problem can not be implemented in the current temperature-based approach. Additionally, a microsegregation model and a coupling relation between thermal and solute equations could be implemented in the code, in order to be able to simulate solidification or melting of binary mixtures. Such models would allow us to simulate more complex configurations such as dendritic formations during the solidification process.

While the Boussinesq approximation proved to be robust for the considered configuration in this study, the assumption of constant thermo-physical properties limits the range of materials that could be simulated by our code. A variable density code could be a valuable tool to investigate a larger type of PCMs and more industrial configurations.

Concerning the numerical model, the second order GEAR scheme provided accurate solutions for all the considered simulations. Furthermore, a variable time step scheme would increase the robustness of the numerical method. The simulation of the melting PCM requires for example a very small time step during the initialization stage, while larger time steps could be prescribed later in the simulation. Adapting automatically this procedure will increase considerably the efficiency of the code. For the solidification problem, new boundary conditions could be further developed by using some models that take into account realistic boundary conditions.

For the parallel algorithm, other preconditionners and solvers should be investigated. This is not a difficult task with FreeFem++ , since the software is interfaced with popular MUMPS, PETSc, or HPDDM solvers.

Résumé

Nous développons dans ce travail de thèse un outil de simulation numérique pour les matériaux à changement de phase (MCP), en tenant compte du phénomène de convection naturelle dans la phase liquide, pour des configurations en deux et trois dimensions. Les équations de Navier-Stokes incompressible avec le modèle de Boussinesq pour la prise en compte des forces de flottabilité liées aux effets thermiques, couplées avec une formulation de l'équation d'énergie suivant la méthode d'enthalpie, sont résolues par une de méthode d'éléments finis adaptatifs. Une approche mono-domaine, consistant à résoudre les mêmes systèmes d'équations dans les phases solide et liquide, est utilisée. La vitesse est ramenée à zéro dans la phase solide, en introduisant un terme de pénalisation dans l'équation de quantité de mouvement, suivant le modèle de Carman-Kozeny, consistant à freiner la vitesse à travers un milieu poreux. Une discrétisation spatiale des équations utilisant des éléments finis de Taylor-Hood, éléments finis P2 pour la vitesse et éléments finis P1 pour la pression, est appliquée, avec un schéma d'intégration en temps implicite d'ordre deux (GEAR). Le sytème d'équations non-linéaires est résolu par un algorithme de Newton. Les méthodes numériques sont implémentées avec le logiciel libre FreeFem++ (www.freefem.org), disponible pour tout système d'exploitation. Les programmes sont distribués sous forme de logiciel libre, sous la forme d'une forme de toolbox simple d'utilisation, permettant à l'utilisateur de rajouter d'autres configurations numériques pour des problèmes avec changement de phase. Nous présentons dans ce manuscrit des cas de validation du code de calcul, en simulant des cas tests bien connus, présentés par ordre de difficulté croissant : convection naturelle de l'air, fusion d'un MCP, le cycle complet fusion-solidification, chauffage par le bas d'un MCP, et enfin, la solidification de l'eau.

Mots-clefs : Navier-Stokes, Boussinesq, méthode d'enthalpie mono-domaine, Carman-Kozeny, éléments finis, FreeFem++, maillage adaptatif, matériaux à changement de phase (MCP), convection naturelle, algorithme de Newton, fusion, cycle complet, solidification de l'eau.

Abstract

In this thesis we develop a numerical simulation tool for computing two and three-dimensional liquid-solid phase-change systems involving natural convection. It consists of solving the incompressible Navier-Stokes equations with Boussinesq approximation for thermal effects combined with an enthalpy-porosity method for the phase-change modeling, using a finite elements method with mesh adaptivity. A single-domain approach is applied by solving the same set of equations over the whole domain. A Carman-Kozeny-type penalty term is added to the momentum equation to bring to zero the velocity in the solid phase through an artificial mushy region. Model equations are discretized using Galerkin triangular finite elements. Piecewise quadratic (P2) finite-elements are used for the velocity and piecewise linear (P1) for the pressure. The coupled system of equations is integrated in time using a second-order Gear scheme. Non-linearities are treated implicitly and the resulting discrete equations are solved using a Newton algorithm. The numerical method is implemented with the finite elements software FreeFem++ (www.freefem.org), available for all existing operating systems. The programs are written and distributed as an easy-to-use open-source toolbox, allowing the user to code new numerical algorithms for similar problems with phase-change. We present several validations, by simulating classical benchmark cases of increasing difficulty: natural convection of air, melting of a phase-change material, a melting-solidification cycle, a basal melting of a phase-change material, and finally, a water freezing case. Keywords : Navier-Stokes, Boussinesq, enthalpy-porosity method, single-domain method, Carman-Kozeny, finite-elements, Taylor-Hood, FreeFem++, mesh adaptivity, phase-change materials (PCM), natural convection, Newton algorithm, GEAR scheme, open-source toolbox, 2D and 3D simulations, melting, solidification, water freezing.

  for the choice of the numerical method . . . . . . . . . . . . . . . . . 3.2 Finite element algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Mesh adaptivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 A finite-element toolbox for the simulation of phase-change systems with natural convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Program architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Input parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3 Outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Numerical tests of the accuracy of the numerical method . . . . . . . . . . . . . . 3.5.1 Space accuracy: Burggraf stationary flow with thermal effects . . . . . . . . 3.5.2 Time accuracy: manufactured unsteady solution . . . . . . . . . . . . . . . 3.6 Domain decomposition method with FreeFem++: FFDDM . . . . . . . . . . . . . 4 Numerical simulation of natural convection flow in a 2D square cavity 4.1 Natural convection of air . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Differentially heated square cavity . . . . . . . . . . . . . . . . . . . . . 4.1.2 Differentially heated cavity with an inner heated square . . . . . . . . . . . 4.2 Natural convection of water inside a square cavity . . . . . . . . . . . . . . . . . . 4.3 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  .6) where ν l = µ l /ρ is the kinematic viscosity, p = (P + ρ ref gy)/ρ ref includes the hydrostatic pressure ρ ref gy and f B (T ) = gβ(T -T ref ) denotes the buoyancy force.

Figure 3 .

 3 Figure 3.1 gives a schematic overview of the content of the toolbox. All files are provided in a directory called PCM-Toolbox. Many detailed comments are included in the programs, with direct link to the mathematical expressions used in this chapter. The FreeFem++ syntax was intentionally kept at a low level of technicality and supplemented with detailed comments when specific more technical syntax was used.

  Figures

Figure 3

 3 Figure 3.3: Burggraf stationary flow with thermal effects used to test the space accuracy of the numerical scheme. Streamlines (a) and temperature contours (b) of the flow field.

Figure 3 Figure 3

 33 Figure 3.4: Burggraf stationary flow with thermal effects used to test the space accuracy of the numerical scheme. Global error ε (3.23) for the temperature θ and velocity u. Numerical scheme using P 1 finite elements for the temperature (panels a and b) or P 2 for the temperature (panels c and d).

Figure 3 . 6 :

 36 Figure 3.6: Time accuracy of the numerical scheme tested using the time-dependent manufactured solution of Nourgaliev et al. (2016). Discretization using P 2 finite elements for the temperature. Evolution of the global error ε (3.23) for the temperature (a) and the velocity (b) measured at t max = π.

Figure 3 . 7 :

 37 Figure 3.7: Partition of Ω = [0, 1] 3 into 8 subdomains with Metis partitioner.

Figure 4 .

 4 Figure 4.1 offers a comparison of the current simulation with the numerical results of[START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF], who used a fourth-order finite difference method for the spatial discretization and a third-order backward finite difference scheme for the time integration, for Ra ranging from 10 4 to 10 6 . The temperature distribution and the streamlines at the steady state display a good agreement with the benchmark solution. The higher the Ra number, the stronger is the natural convection flow in the cavity. A single convection cell is observed in the center of the cavity for Ra = 10 4 (Fig.4.1d) during which the heat transfer is merely bulk heat transfer. Figures 4.1e and 4.1f however exhibit stronger convection with more rolls, when the boundary layer heat transfer dominates.

Figure 4

 4 Figure 4.1: Natural convection of air for Ra ranging from 10 4 to 10 6 (from left to right) and Pr = 0.71. Temperature fields and streamlines at the steady state. Panels (a) to (c) correspond to the benchmark solution of[START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF]. Panels (d) to (f) display our numerical results.

  Figure 4.2: Natural convection of air in a differentially heated cavity for Ra ranging from 10 4 to 10 6 and Pr = 0.71. (a) Transversal velocity profile along the horizontal symmetry line. (b) Longitudinal velocity profile along the vertical symmetry line. Numerical results obtained using the present Newton method (symbols) with a mesh resolution of nbseg = 80; comparison with the spectral-accurate simulations by Le Queré (1991) (solid lines).

Figure 4

 4 Figure 4.3: Natural convection of air. Temperature isolines at the steady state for (a) Ra = 10 7 , (b) Ra = 4 • 10 7 , (c) and Ra = 10 8 .

Figure 4

 4 Figure 4.4: Natural convection of air. Time evolution of the Nusselt number at x = 0 and x = 1/2 for Ra = 1 • 10 7 .

Figure 4

 4 Figure 4.5: Natural convection of air in a differentially heated cavity with an inner heated square for Ra ranging from 10 4 to 10 6 . Temperature field and streamlines (top) and transversal velocity profile along the horizontal symmetry lines (bottom). Results obtained using the present Newton method (red solid line), with mesh resolution nbseg = 80; comparison with the finite difference code of Moglan (2013).

Figure 4

 4 Figure 4.6: Natural convection of water in a differentially heated cavity with non-dimensional parameters: Ra = 2.518084 • 10 6 and Pr = 6.99. (a) Iso-lines of the temperature at the steady state. (b) Streamline of the steady flow. (c) Illustration of the mesh adaptivity. The mesh is refined along the dimensionless temperature iso-line θ = θ m , due to the density variation. (d) Temperature profile along the horizontal symmetry line. Comparison with the numerical results of Michalek and[START_REF] Michalek | Simulations of the water freezing process -Numerical benchmarks[END_REF].

  6 1.62935 • 10 11 -9.8058 • 10 -14 -3.20216 • 10 -10 10 10 1.62937 • 10 15 1.09832 • 10 -17 -3.23176 • 10 -14 10 20 1.62935 • 10 25 2.7011 • 10 -26 -3.18493 • 10 -24

Figure 5 . 1 :

 51 Figure 5.1: Location of the interface at t = 78.7 for the melting of the PCM. (a) Computations using different values of C CK ranging from 10 6 to 10 20 and comparison with experimental data of[START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]. (b) Velocity profile following the line y = 0.9 with respect to the temperature.

Figure 5

 5 Figure 5.2: Location of the interface at t = 78.7 for the melting of the PCM. Computations using different values of C CK ranging from 10 6 to 10 20 and different values of the smoothing radius R ck : (a) R ck = ε/2, (b) R ck = ε/4.

Figure 5

 5 Figure 5.3: Adapted mesh during simulation of PCM melting. (a) 4, 211 triangles: the mesh is refined around the melting front θ = 0 and in the fluid where velocity gradients exist. (b) zoom showing the variation of the mesh size around the temperature isoline θ = 0.

  .4. -Mesh #1 (Panel a): N t = 48, 247 triangles and 314, 203 degree of freedom (d.o.f), -Mesh #2 (Panel b): N t = 21, 093 triangles and 143, 202 d.o.f, -Mesh #3 (Panel c): N t = 2, 870 triangles and 318, 974 d.o.f, -Mesh #4 (Panel d): N t = 820 triangles and 5, 544 d.o.f.

Figure 5

 5 Figure 5.4: Benchmark #1. Influence of the mesh resolution on the numerical solution at dimensionless time t = 78.7. The mesh is refined around the melting front θ = 0 and in the fluid where velocity gradients exist. Solid lines denote the solid-liquid interface. Number of triangles: (a) 48, 247, (b) 21, 093, (c) 2, 870, (d) 820. The location of the interface (e) and the time evolution of the Nusselt number (f) are compared. 51

Figure 5 . 5 :

 55 Figure 5.5: Benchmark #1. (a) Temperature isoline at t = 78.8. (b) Location of the interface.Comparison with experimental data of[START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and numerical solutions of[START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] andWang et al. (2010a) for two time instants (t = 39.9 and 78.8).

  Wang et al. (2010a) used a FV code with a fixed mesh of resolution [60 × 60], and Danaila et al. (2014) used a FE method with viscosity penalty.

Figure 5 . 6 :

 56 Figure 5.6: Benchmark #2. (a) Temperature distribution, streamlines, and solid-liquid interface at dimensionless time t = 1.759. (b) Comparison of the location of the solid-liquid interface, with both experimental and simulation results of Gong et al. (2015), for five time instants (t = 0.153, 0.347, 0.465, 0.869, and 1.759).

Figure 5 . 7 :

 57 Figure 5.7: Benchmark #3. Location of the solid-liquid interface at dimensionless time (panels a to d) t = 0.25, t = 1, t = 3, and t = 5, compared with five simulations presented by Bertrand et al. (1999).

Figure 5 . 8 :

 58 Figure 5.8: Benchmark #3: Ra = 2 • 10 8 , Pr = 50 and Ste = 0.1. Time evolution of the liquid fraction (a) and the Nusselt number (b).Comparison with five simulations presented by[START_REF] Bertrand | Melting driven by natural convection A comparison exercise: first results[END_REF].

Figure 5

 5 Figure5.9: Melting of Gallium: temperature field, streamlines, and melting front for dimensionless time instants (panels a to d): t = 0.0015, 0.006, 0.01, and 0.019. Adapted mesh corresponding to t = 0.019 (e). Time evolution of the Nusselt number (f).For a better view of the convection cells, a ratio 2 : 1 was used for the axes.

  the outer radius (R i = R/4), for four heated tubes R i = R/8 and for nine heated tubes R i = R/12. A Dirichlet boundary condition is applied to inner tubes θ = θ h , and a Neumann boundary condition ∂θ ∂n = 0 is used at the outers. Homogenous Dirichlet boundary conditions (u = 0) are prescribed on the velocity throughout all ∂Ω. The physical parameters of the simulation are: Ra = 5 • 10 4 , Pr = 0.2 and Ste = 0.02. Only the half of the domain is simulated since the problem is axisymmetric for the investigated Ra number. At t = 0, the PCM is completely solid with isothermal temperature θ 0 lower than the temperature of fusion. Immediately after the pipe placed on the center of the enclosure is heated, the melted layer expands in the radial direction. The mesh is refined initially around inner tubes, and is dynamically adapted at each time step around the melting front and the thermal boundary layer area.

Figure 5 .

 5 Figure 5.10: Temperature fields for the melting of a cylindrical PCM with inner heated tubes. Time instants corresponding to the same liquid fraction L f = 80%. Configurations with (a) one tube (t = 2.5), (b) four tubes (t = 0.99) and (c) nine tubes (t = 0.4). Melting fronts are represented with black lines (only half of the domain is simulated).

Figure 5

 5 Figure 5.11: Time evolution of the liquid fraction for one, four, and nine heated tubes. (a) Comparison with numerical results of Luo et al. (2015) (b) comparison of different arrangements: central symmetry, inline, and staggered.

Figure 5 .

 5 Figure 5.12: Solid crust formation in a distorted mesh. Temperature field and streamlines in our simulation (a) and Nourgaliev et al. (2016) (b).

  4.6), freezing starts by dropping progressively the temperature of the cold (right) wall from T c = 0 • C to T c = -10 • C. The dimensionless parameters describing the problem are: Ra = 2.518084 • 10 6 , Pr = 6.99 and Ste = 0.125.

Figure 5 .

 5 Figure 5.13: Freezing of pure water. Configuration at physical time t ϕ = 2340[s] (t = 1.61): (a) experimental image from Kowalewski and Rebow (1999); the thick red line represents the solid-liquid interface computed with the present method, (b) computed streamlines showing the two recirculating zones in the fluid phase, (c) finite-element mesh refined along the solid-liquid interface (T = 0 o C) and also along the line of maximum water density (T = 4 o C) (d) location of the solidification front for several time instants: from t = 0.05 (corresponding to physical time t ϕ = 70[s]) to t = 1.61 (2340[s])

Figure 6

 6 Figure 6.1: Comparison between lateral (LM) and basal melting (BM). Sketch of the computational domain and boundary conditions for LM (top) and BM (bottom) cases. Streamlines (panels b and d) and solid-liquid interfaces (solid red lines).

Figure 6 . 2 :

 62 Figure 6.2: Temperature iso-lines, streamlines in the fluid phase, and phase-change interface. The solid part is represented in blue and corresponds to the region of temperature θ ≤ θ f = 0. Time instants (panels a to f): τ = 0.004; 0.016; 0.032; 0.063; 0.08; 0.2. Ra = 3.27 • 10 5 , Pr = 56.2 and Ste = 0.045.

  time instants. Panels (a) to (f) depict the flow streamline showing the clockwise recirculation of the fluid, the melting front, and the temperature distribution (the solid phase is denoted by the blue region).

Figure 6

 6 Figure 6.3: Illustration of the mixed regime by Jany and Bejan (1988). Solid red line is the solid-liquid interface, Γ i represents the location of the interface and y * denotes the height of fluid impacted by the emerging convective flow.

Figure 6

 6 Figure 6.4: Complete melting of the PCM. Time evolution of the average Nusselt number defined at the hot (left) wall (cf. Eq. (2.24)) (solid line). Comparison with the experimental results of Okada (1984) (dashed line) and the predictions using the correlation in Eq. (6.11) suggested by Jany and Bejan (1988) (dash-dot line). Ra = 3.27 • 10 5 , Pr = 56.2 and Ste = 0.045.

Figure 6 . 5 :

 65 Figure 6.5: Complete melting of the PCM. (a) Time evolution of the liquid fraction for the complete melting of the PCM. (b) Comparison of our results (solid line) with the numerical results of Wang et al. (2010a) (dashed line) and the predictions using the correlation (6.13) suggested by Jany and Bejan (1988) (dash-dot line).

  27 • 10 5 by a factor of 5 and 10, respectively. The exact values are: Ra = 1.62 • 10 6 and Ra = 3.27 • 10 6 . First, we increase the height H of the cavity by a factor of 3 √ 5 and 3 √ 10 and consider the same δT . Thus the Ste number is kept constant. Second, we increase the temperature difference parameter δT by keeping H constant. It corresponds to an increased value of the Stefan number by a factor of 5 and 10: Ste = 0.223 and Ste = 0.45.

Figure 6

 6 Figure 6.6: PCM melting at L f = 0.5. Illustration of the temperature field, streamlines, and melting front for three Ra numbers: (a) Ra = 3.27 • 10 5 , (b) Ra = 1.62 • 10 6 , and (c) Ra = 3.27 • 10 6 . The Pr and Ste numbers are kept constant: Pr = 56.2 and Ste = 0.045. Snapshots of numerical solutions for Ra = 3.27 • 10 5 , Ra = 1.62 • 10 6 and Ra = 3.27 • 10 6 at constant Ste number are shown in Fig. 6.6. The colours correspond to the temperature distribution, the black lines correspond to streamlines and the white lines correspond to the solid-liquid interface.

Figure 6

 6 Figure 6.7: Complete melting of the PCM. Influence of the value of the Rayleigh number (Ra) on the time evolution of the liquid fraction (a) and the average Nusselt number defined at the hot (left) wall (b). The reference case (Ra = 3.27 • 10 5 ) is represented by red continuous lines. The value of the Ra was increased by a factor of 5 and 10, while the Stefan number Ste is kept constant.

Figure 6

 6 Figure 6.8: Complete melting of the PCM. Influence of the value of the Rayleigh number (Ra) on the time evolution of the liquid fraction (a) and the average Nusselt number defined at the hot (left) wall (b). The reference case (Ra = 3.27 • 10 5 ) is represented by red continuous lines. The value of the Ra and Ste were increased by a factor of 5 and 10.

Figure 6

 6 Figure6.8a shows that by increasing the value of δT , and consequently increasing the Rayleigh number and the Stefan number, the PCM melts faster. We note that the height H of the cavity is kept constant, hence the natural convection flow in the melted PCM is enhanced when the Rayleigh number keep increasing. As a consequence, the convection-dominated regime is reached earlier, as shown in Fig.6.8b by the shift of the minimum of the N u to lower values of t ϕ . This evolution is also observed for the liquid fraction. As expected, an increase of the Rayleigh number and the Stefan number is followed by an enhancement of the heat transfer during the melting, and consequently an improved efficiency of the PCM.

Figure 6 . 9 :

 69 Figure 6.9: Melting of PCM heated from below. Temperature field and solid-liquid interface (top) and streamlines (bottom) in the liquid phase for different Rayleigh numbers (a) Ra = 3.27•10 5 and t = 30, (b) Ra = 1.635•10 6 and t = 15, (c) Ra = 3.27•10 6 and t = 10.

  Figure 6.11: Natural convection flow emerging from melting heated from below.

  23) could be hence simplified by normalising Eq. (6.23) with respect to the friction scale, fluids, the momentum balance is between buoyancy and friction (see e. g.[START_REF] Bejan | Convection heat transfer[END_REF]; Le Quéré and Gobin (1999)) and we obtain

Figure 6

 6 Figure 6.12: BM case for three Rayleigh numbers: Ra = 3.27 • 10 5 , 1.62 • 10 6 , and 3.27 • 10 6 . Time evolution of the liquid fraction (a) and Nusselt number (b).

Figure 6

 6 Figure 6.13: Comparison between LM and BM cases. Rayleigh numbers ranging from Ra = 3.27 • 10 5 to Ra = 6.54 • 10 6 . Liquid fraction(left), Nusselt number (middle), and accumulated heat input (right).

Figure 7

 7 Figure 7.1: Sketch of the computational domain and boundary conditions. General configuration (panel a) with isothermal (θ = cst.) vertical (x = 0 and x = 1) walls and adiabatic (∂θ/∂n = 0) top and bottom walls. Configuration for the melting phase (panel b) with a hot left wall (θ = θ h > 0) and a cold right wall (θ = θ c < 0), followed by a solidification phase (panel c), when the temperature of the left wall is cooled to θ = θ co < 0.

Figure 7 . 2 :

 72 Figure 7.2: Solidification of the PCM: case CM. Temperature iso-lines in the liquid phase. The solid part is represented in blue and corresponds to the region of temperature θ co ≤ θ f = 0. Time instants (panels a to e): t ϕ = 185 min, t ϕ = 231 min, t ϕ = 300 min, t ϕ = 430 min and t ϕ = 510 min. The adapted mesh corresponding to t ϕ = 300 min is plotted in panel (f). Ra co = 3.27 • 10 5 .

Figure 7 . 3 :

 73 Figure 7.3: Solidification of the PCM: case PM. Temperature iso-lines in the liquid phase. The solid part is represented in blue and corresponds to the region of temperature θ co ≤ θ f = 0. Time instants (panels a to e): t ϕ = 59 min, t ϕ = 70 min, t ϕ = 90 min, t ϕ = 131 min and t ϕ = 200 min. The adapted mesh corresponding to t ϕ = 90 min is plotted in panel (f). Ra co = 3.27 • 10 5 .

Figure 7

 7 Figure 7.4: Temporal evolution of the liquid fraction (L f ), the Nusselt number Nu, and the accumulated heat input Q 0 during the entire melting-solidification cycle. Case CM (left) and case PM (right).

Figure 7

 7 Figure 7.5: Solidification of a cylindrical PCM. Temperature iso-lines in the liquid phase. The solid part is represented in blue and corresponds to the region of temperature θ co ≤ θ f = 0. Time instants (panels a to e): t = 2.69, t = 3.16, t = 4.34 , t = 4.94 and t = 6.66. The adapted mesh corresponding to t = 3.16 is plotted in panel (f). Ra co = 5 • 10 4 .

Figure 8

 8 Figure 8.1: Natural convection of air in a 3D differentially heated cavity. Temperature contours at the mid-plane of (y = 0.5). Comparison with the results of Wakashima and Saitoh (2004) (left images).

Figure 8

 8 Figure 8.2: Natural convection of water in a cube-shape cavity. (a) Temperature distribution and adapted mesh. (b) Zoom of the 3D adapted mesh at mid-plane. (c) Iso-surface θ = θ m . Comparison of the profile of the vertical velocity along the z-direction with the numerical benchmark by[START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF] at the mid-plane (y = 0.5) and x = 0.93 (d); x = 0.1 (e); x = 0.5 (f).

  .3. The physical parameters of the run are: Ra = 3.27 • 10 5 , Pr = 56.2, and Ste = 0.045. The PCM is set initially solid throughout the whole domain. As the heating progresses, the natural convection intensifies enough to have a pronounced influence on the shape of the interface. The shape of the phase-change front at t = 78.7 is displayed in Fig. 8.3, showing a nonuniform melting front receding from the top to the bottom of the domain.

Figure 8 . 3 :

 83 Figure 8.3: Melting of a phase-change material. Problem definition, shape of the solid-liquid interface at t = 78.8, and temperature distribution in the symmetry plane (a). Scheme of 3D convection flow near the side wall by Nikrityuk (2012) (b).

Figure 8

 8 Figure 8.4: Melting of a phase-change material. Temperature iso-surfaces (a) and the corresponding adapted mesh (b) at t = 78.8.

Figure 8

 8 Figure 8.5: Comparison between 2D and 3D solutions. Location of the interface at t = 78.7 (a) and the Nusselt number at x = 0 (b).

Figure 8 . 6 :

 86 Figure 8.6: Freezing of pure water in a 3D cubic cavity. (a) Temperature distribution and adapted mesh at t ϕ = 2340[s]. (b) Zoom of the adapted mesh at the mid-plane. (c) Temperature iso-surfaces θ = 0.4 and θ = 0 with velocity vectors in the liquid phase. (d) Superimposition of the interface obtained in the present simulation (red line) on the experimental image of Kowalewski and Rebow (1999).

  Green energy management and phase-change materials . . . . . . . . . . . . . . . 1.2 Theoretical and experimental studies of phase-change systems . . . . . . . . . . . . 1.3 Physical models for phase-change materials . . . . . . . . . . . . . . . . . . . . . 1.4 Numerical approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Present numerical approach . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Contents
	1 Introduction
	1.1

1.6 Purpose of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Thesis plan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Navier-Stokes-Boussinesq and enthalpy-porosity models 2.1 Enthalpy-porosity model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Navier-Stokes equations with Boussinesq approximation . . . . . . . . . . . . . . 2.3 Dimensionless system of equations for the single-domain approach . . . . . . . . . 2.4 Boundary layer approximation and scale analysis . . . . . . . . . . . . . . . . . .

Table 1

 1

	1.1.
	Agriculture Buildings Industry Transportation Other
	4.5	41.8	26	43.8	4.5

.1: Energy consumed by different domains in France in 2017, in million tonnes of oil equivalent. French ministry of ecological transition (Ministère de la transition écologique et solidaire, 2018).

  27 • 10 -17 6.70 • 10 -18 5.61 • 10 -18 3.54 • 10 -18 Ω p 1.15 • 10 -16 5.40 • 10 -14 4.61 • 10 -10 4.37 • 10 -7

Table 3 .

 3 1: Average values of the pressure field and of the divergence of the velocity field, computed for different values of the penalty parameter γ.

Table 4

 4

.1: Physical parameters of air and water at T = 300K used in our simulations. Pr = 0.71 (for air) and Pr = 6.99 (for water).

  617468 • 10 -4 2.418556 • 10 -4 7.203805 • 10 -4 80 1.961756 • 10 -5 5.463075 • 10 -5 3.153492 • 10 -5 9.160633 • 10 -5 120 4.478622 • 10 -6 1.518429 • 10 -5 9.539112 • 10 -6 2.891752 • 10 -5 160 2.186721 • 10 -6 6.181043 • 10 -6 2.802694 • 10 -6 9.789997 • 10 -6 200 1.20101 • 10 -6 3.666731 • 10 -6 2.282512 • 10 -6 5.690341 • 10 -6

Table 4.2: Error associated to the spatial discretization. Comparison with Le Queré (1991).

Table 4

 4

		u max y	v max x	Nu 0	Nu 1/2	Nu max y	Nu min y
	Ra = 10 4	1.61782 0.826772 0.11811 1.96272 2.24487 2.24403	3.53218 0.141732	0.58497 1
	de Vahl Davis (1983)	1.6178 0.823	1.9617 0.119	2.238	2.243 3.528 0.149	0.586 1
	Ra = 10 5	3.47218 0.858268 0.0629921 6.85436 4.52172 4.52176	7.72 0.0787402	0.72797 1
	de Vahl Davis (1983)	3.473 0.855	6.859 0.066	4.519	4.519	7.717 0.085	0.729 1
	Ra = 10 6	6.48369 0.850394 0.0393701 22.0306 8.82521 8.82513	17.5344 0.0393701	0.979543 1
	de Vahl Davis (1983)	6.463 0.850	21.936 0.0379	8.817	8.799	17.925 0.0378	0.989 1
	Le Queré (1991)	6.48344 0.850	22.0559 8.8252 0.038	8.8252	17.5360 0.039	0.97946 1
	Ra = 10 7	4.69349 0.88189	2.19517 16.5237 16.5194 0.023622	39.3455 0.015748	1.36675 1
	Le Queré (1991)	4.69862 0.879	2.21118 16.523 0.021	16.523	39.3947 0.018	1.36635 1
	Ra = 10 8	3.20182 0.929134 0.015748 2.09841 30.2353 30.1604	87.3469 0.00787402	1.92085 1
	Le Queré (1991)	3.21875 0.928	2.22239 30.225 0.012	30.225	87.2355 0.008	1.91907 1

.3: Natural convection of air in a square cavity. Comparison with the results by

[START_REF] De Vahl | Natural convection of air in a square cavity: a benchmark numerical solution[END_REF]

and Le

  CaseNb of iterations Nb of triangles Nb of d.o.f. CPU time (s)Table 4.4: Description of the runs for natural convection cases with Ra = 10 6 .

	Air	unsteady steady	164 6	2, 413 2, 312	16, 133 15, 558	338.527 14.0465
	Air with obstacle	unsteady steady	157 6	3, 489 2, 967	23, 354 19, 933	436.405 15.9369
	Water	unsteady unsteady	1, 164 20	3, 422 2, 960	22, 639 19, 629	3162.86 51.1662

Table 5 .

 5 1: Physical properties of n-octadecane and Gallium PCM.

  Table 5.5: Summary of simulations for the melting PCM: number of triangles and CPU times.

				using a fixed-grid
	algorithm.
	Case		CPU time (s) Number of triangles
		Bench #1	2876.39	2, 900
	N-octadecane	Bench #2	19231	4, 198
		Bench #3	40522.7	7, 000
	Gallium	14621.2	4, 820
	Cylindrical PCM	2060.42	3, 076
	Highly distorted PCM	6657.58	2, 769

8 Numerical simulation of 3D convection and PCM configurations using domain decomposition methods

  Several sections of this chapter are from the paper [G. Sadaka, A. Rakotondrandisa, F. Luddens, C. Lothodé, P-H. Tournier, I. Danaila, Parallel finite-element codes for the simulation of solid-liquid phase-change systems with natural convection, to be submitted, 2019].

  50094 • 10 -6 on θ.Table 8.2: Natural convection of air in a 3D differentially heated cavity. Comparison between sequential and DDM algorithms for uniform grids of 40 × 40 × 40 points.

	Ra nb proc 10 4 28 42 56 70 10 5 28 42 56 70 10 6 28 42 56 70	||u|| 2 1.12496 • 10 -6 1.53698 • 10 -6 1.55576 • 10 -6 1.25622 • 10 -6 1.73254 • 10 -6 2.84973 • 10 -6 3.00832 • 10 -6 3.68118 • 10 -6 6.61804 • 10 -6 1.826 • 10 -5 3.46504 • 10 -6 1.1 • 10 -5 ||u|| ∞ ||T || 2 ||T || ∞ 3.1 • 10 -6 3.09966 • 10 -6 7 • 10 -6 5.1 • 10 -6 3.23352 • 10 -6 8 • 10 -6 5.1 • 10 -6 3.4342 • 10 -6 8 • 10 -6 3.6 • 10 -6 3.56048 • 10 -6 8 • 10 -6 6.1 • 10 -6 2.40467 • 10 -6 7 • 10 -6 7.78 • 10 -6 3.53003 • 10 -6 9 • 10 -6 7.39 • 10 -6 4.17769 • 10 -6 1.1 • 10 -5 9 • 10 -6 4.70846 • 10 -6 1.2 • 10 -5 5.93966 • 10 -6 1.5 • 10 -5 3.98082 • 10 -6 1.2 • 10 -5 7.05144 • 10 -6 1.9247 • 10 -5 5.0044 • 10 -6 2 • 10 -5 6.02152 • 10 -6 1.68 • 10 -5 4.50094 • 10 -6 1.8 • 10 -5

FreeFem++ for different OS can be downloaded from http://www3.freefem.org/.

http://lmrs-num.math.cnrs.fr

In the liquid phase, ∇ • (csU ) = h sl ∇ • U = 0; in the solid phase, s = 0; in the regularization (mushy) region, it is assumed that U = 0.

Remerciements

C'est avec une grande émotion que j'adresse premièrement mes remerciements à mon directeur de thèse, Ionut Danaila. Je tiens à lui exprimer ma gratitude pour la confiance, la patience, ainsi que l'intérêt qu'il m'a témoigné tout au long de ce travail de thèse. L'enthousiasme et le dynamisme qu'il a su communiquer au sein du groupe calcul scientifique ont été pour moi des sources de motivations. J'ai beaucoup appris au-delà de l'aspect scientifique par la rigueur qu'il a su transmettre : "on vise la qualité professionnelle en toute chose"! Merci de m'avoir accordé cette opportunité.

Ma reconnaissance va aussi à Zakaria Belhachmi et Patrick Le Quéré d'avoir accepté d'être les rapporteurs de ma thèse, ainsi qu'à Pascal Frey, Frédéric Hecht, Stéphane Le Masson, Luminita Danaila, et Francky Luddens, qui m'ont fait l'honneur de participer au jury de cette thèse et de porter une attention avisée sur mon travail.