Aina Rakotondrandisa
email: aina.rakotondrandisa@etu.univ-rouen.fr

Georges Sadaka
email: georges.sadaka@univ-rouen.fr

Ionut Danaila
email: ionut.danaila@univ-rouen.fr

A finite-element toolbox for the simulation

Keywords: melting, solidification, PCM, Newton method, finite element, mesh adaptivity, Navier-Stokes, Boussinesq, FreeFem++ Programming language: FreeFem++ (free software, www3.freefem.org) Computer: PC, Mac, Super-computer. Operating system: Windows, Mac OS, Linux. Classification:

We present and distribute a new numerical system using classical finite elements with mesh adaptivity for computing two-dimensional liquid-solid phase-change systems involving natural convection. The programs are written as a toolbox for FreeFem++ (www3.freefem.org), a free finite-element software available for all existing operating systems. The code implements a single domain approach. The same set of equations is solved in both liquid and solid phases: the incompressible Navier-Stokes equations with Boussinesq approximation for thermal effects. This model describes naturally the evolution of the liquid flow which is dominated by convection effects. To make it valid also in the solid phase, a Carman-Kozeny-type penalty term is added to the momentum equations. The penalty term brings progressively (through an artificial mushy region) the velocity to zero into the solid. The energy equation is also modified to be valid in both phases using an enthalpy (temperature-transform) model introducing a regularized latent-heat term. 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. For the temperature both P2 or P1 discretizations are possible. 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. An efficient mesh adaptivity algorithm using metrics control is used to adapt the mesh every time step. This allows us to accurately capture multiple solid-liquid interfaces present in the domain, the boundary-layer structure at the walls and the unsteady convection cells in the liquid. We present several validations of the toolbox, by simulating benchmark cases of increasing difficulty: natural convection of air, natural convection of water, melting of a phase-change material, a melting-solidification cycle, and, finally, a water freezing case. Other similar cases could be easily simulated with this toolbox, since the code structure is extremely versatile and the syntax very close to the mathematical formulation of the model.

Introduction

Solid-liquid phase-change problems are encountered in numerous practical applications, ranging from metal casting and thermal energy storage (phase-change materials) to food freezing and cryosurgery. Melting and solidification are also fundamental processes in geophysical problems, such as Earth's mantle formation, lava lakes or magma chambers. In many of these problems, temperature gradients induce buoyancy forces in the liquid (melted) phase generating a significant convective flow. Consequently, the appropriate mathematical description of the liquid phase is the usual model for the natural convection flow: the incompressible Navier-Stokes system of equations with Boussinesq approximation for thermal (buoyancy) effects. 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 below a short description of the two main approaches suggested in the literature to deal with this 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]. 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 (see the solidification of a phase-change material presented in this paper) 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 post-processing 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 3 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 temperature-based 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 | An overview of numerical methods for solving phase change problems[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 | An overview of numerical methods for solving phase change problems[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 fixed-grid methods for the modeling of melting[END_REF].

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 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 [START_REF] Gartling | Finite element analysis of convective heat transfer problems with change of phase[END_REF]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 were theoretically 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]. A new formulation of the phase-change problem using as variables the pseudo-stress, strain rate and velocity for the Navier-Stokes-Brinkman equations was recently suggested by [START_REF] Alvarez | New mixed finite element methods for natural convection with phase-change in porous media[END_REF].

A third technique used to ensure a zero velocity field in the solid phase is the so-called enthalpy-porosity 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 is a considerable increase of the computational time, even for two dimensional cases. Finite element (FE) methods offer the possibility to dynamically refine the mesh only in the regions of the domain where sharp phenomena take place (e. g. solid-liquid interface, recirculation zones). 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 suggested for classical two-phase Stefan problem and phase-change systems with convection. In [START_REF] Belhamadia | Anisotropic mesh adaptation for the solution of the Stefan problem[END_REF] a hierarchical error estimator was used for 2D Stefan problems (without convection). A different mesh adaptivity method, based on the definition of edge length from a solution dependent metric, was used in [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF] to deal with the same Stefan problems in 3D and also with 2D phase-change systems with convection. The time-dependent mesh adaptivity strategy based on metrics control suggested by [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] proved very effective in refining the mesh near walls and interfaces, when simulating melting and solidification phase-change systems (see also [START_REF] Rakotondrandisa | Numerical modelling of a melting-solidification cycle of a phase-change material with complete or partial melting[END_REF]. [START_REF] Zimmerman | Monolithic simulation of convection-coupled phase-change: Verification and reproducibility[END_REF] implemented the variable viscosity model suggested by [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] in a different finite-element framework using FEniCS and an AMR (Adaptive Mesh Refinement) technique based on a dual-weighted residual method. In a very recent contribution, [START_REF] Belhamadia | A two-dimensional adaptive remeshing method for solving melting and solidification problems with convection[END_REF] derived a time-dependent adaptive remeshing method for phase-change problem with convection based on an error estimator applicable to second or higher order variables.

However, to the best of our knowledge, no adaptive finite-element programs exist in the CPC Program Library for phase-change problems. The purpose of this paper is therefore to distribute a finite-element toolbox to solve two-dimensional solid-liquid phase-change problems.

The present toolbox is based on a single-domain enthalpy-porosity model for solidliquid 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. The choice of using Taylor-Hood finite elements for the fluid flow and P 2 for the temperature was also made in [START_REF] Woodfield | Stability and finite element approximation of phase change models for natural convection in porous media[END_REF]; [START_REF] Belhamadia | A two-dimensional adaptive remeshing method for solving melting and solidification problems with convection[END_REF]. It offers second order accuracy and (inf-sup) stability of the space discretization of the fluid flow. Other possible choices for the finite-element discretization could be tested for the implementation of the algorithm presented in this paper. For instance, the mini-element (see [START_REF] Arnold | A stable finite element for the Stokes equations[END_REF][START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF], offering a global linear convergence, could be an interesting alternative to replace the Taylor-Hood element and thus optimize the computational time. This possibility could be explored with the present toolbox, since the mini-element exists in FreeFem++ .

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 feature of the toolbox offers the possibility to deal with complicated phase-change cases, involving multiple solid-liquid interfaces. There are two 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]).

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 code new numerical algorithms for similar problems with phase-change.

The paper is organized as follows. Section 2 introduces the governing equations and the regularization technique used in the enthaphy-porosity model. Section 3 presents the Newton algorithm for the Navier-Stokes-Boussinesq system of equations. The mesh adaptivity technique is also discussed in this section. A description of the programs is given in Section 4. The accuracy of the numerical method is tested in Section 5 using manufactured solutions. Finally, Section 6 is devoted to extensive numerical validations of the method. The robustness of the algorithm is demonstrated by comparing our results with well defined classical benchmarks. The capabilities of the toolbox to deal with complex geometries is also illustrated. The main features of the software and possible extensions are summarised in Section 7.

Enthalpy-porosity model

We consider a solid-liquid system placed in a two-dimensional domain Ω. In the following, subscripts s and l will refer to the solid and the liquid phases, respectively. We present in this section the single-domain model, using the same system of equations in both liquid and solid phases.

Energy equation

The phase change process is modelled using an enthalpy method with a temperaturebased 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 with the classical energy equation:

∂(ρh) ∂t ϕ + ∇ • (ρhU ) -∇ • (k∇T ) = 0, (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)
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 this paper 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]). Equation (1) can be further simplified under the following assumptions: (i) the density difference between solid and liquid phases is negligible, i. e. ρ l = ρ s = ρ is constant; (ii) the mushy region is narrow and the velocity inside this zone is zero. Consequently, the final temperature transforming model is obtained by combining (2) and (1) and neglecting the convection term ∇ • (csU )2 :

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

Navier-Stokes equations with Boussinesq approximation

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

∂U ∂t ϕ + (U • ∇)U + 1 ρ ∇p - µ l ρ ∇ 2 U -f B (T )e y = A(T )U , (4)
where p denotes the pressure, µ l the dynamic viscosity of the liquid (assumed to be constant) and f B (T ) the Boussinesq force. The Carman-Kozeny penalty term A(T )u is artificially introduced in (4) to extend this equation into 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 equation (4) to A(T )u = 0, equivalent to u = 0. Exact expressions for f B and A will be given in the next section.

Finally, the conservation of mass in the liquid phase is expressed by the continuity equation:

∇ • U = 0.
(5)

Final system of equations for the single-domain approach

The present numerical code solves a dimensionless form of the previous equations. After choosing a reference length H (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 H , u = U V ref , θ = T -T ref δT , t = V ref H t ϕ . ( 6
)
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. Temperature difference δT defines a temperature scale, that will be set differently for melting and solidification cases.

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

∇ • u = 0, (7) ∂u ∂t + (u • ∇)u + ∇p - 1 Re ∇ 2 u -f B (θ) e y -A(θ)u = 0, (8)
∂ (Cθ) ∂t + ∇ • (Cθu) -∇ • K ReP r ∇θ + ∂ (CS) ∂t = 0. ( 9
)
The linearised Boussinesq buoyancy force (f B ), the Reynolds (Re) and Prandtl (P r) numbers are defined as:

f B (θ) = Ra P rRe 2 θ, Re = ρV ref H µ l = V ref H ν l , P r = ν l α l , ( 10
)
with ν the kinematic viscosity and α = k/(ρc) the thermal diffusivity. In the expression of f B , the Rayleigh number of the flow is defined as:

Ra = gβH 3 δT ν l α l , (11)
with β the thermal expansion coefficient and g the gravitational acceleration.

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 . ( 12
)
The non-dimensional function S = s/s l in the energy equation ( 9) takes a similar nondimensional form:

S(θ) = s s l =    h sl /c l δT = 1 Ste , θ ≥ θ f , 0, θ < θ f , (13)
with Ste the Stefan number. Discontinuous step-functions defined in ( 12) and ( 13) 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 s , θ s , R s ) = f l + f s -f l 2 1 + tanh a s θ s -θ R s , (14)
where f l , f s are the imposed values in the liquid and solid phases, a s a smoothing parameter, θ s the central value (around which we regularize) and R s the smoothing radius. For example, we use for the non-dimensional source term in (9) the following regularisation over the artificial mushy region θ ∈ [-ε, ε]:

S(θ) = 1 Ste - 1 2 Ste 1 + tanh θ r -θ R s , ( 15
)
where θ r is the central value around which we regularize (typically θ r = θ f = 0) and R r the smoothing radius (typically R s = ε).

Finally, the penalty term A(θ)u in momentum equation ( 8) follows from the Carman-Kozeny model (Voller et al., 1987;[START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF][START_REF] Kheirabadi | The effect of the mushy-zone constant on simulated phase change heat transfer[END_REF]:

A(θ) = - C CK (1 -L f (θ)) 2 L f (θ) 3 + b , ( 16
)
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 (15). 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.

Numerical method

To solve the system of equations ( 7)-( 9) 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]). Among the numerous numerical tools offered by FreeFem++ , the use of the powerful mesh adaptivity function proved mandatory in this study 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 artificial mushy region, offering an optimal resolution of the large gradients of all regularized functions (S, K, L f ), while the mesh was coarsened (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 pressure with P 1 finite elements (piecewise linear). The temperature and the enthalpy are discretized using either P 1 or P 2 finite elements. For the time integration, we use a second-order Gear (BDF2) scheme (see also [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF]):

dφ dt 3φ n+1 -4φ n + φ n-1 2δt , (17)
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 equations ( 7)-( 9) are treated implicitly (i. e. taken at time t n+1 ). The resulting non-linear equations are solved using a Newton algorithm.

Finite element algorithm

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. Theory and numerical 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]. To simplify the presentation, we consider in the following that C = 1 and K = 1. For phase-change materials like paraffins this is a physically valid assumption. Nevertheless, we keep in the equations below the variable K, since for the water freezing case it is necessary to take into account the dependence K(θ). The derivation of supplementary terms associated to K(θ) is straightforward.

We use the second-order finite-difference scheme ( 17) and obtain the following implicit semi-discretization in time of the single-domain model ( 7)-( 9):

∇ • u n+1 = 0, (18) 3 2 u n+1 δt + (u n+1 • ∇)u n+1 + ∇p n+1 - 1 Re ∇ 2 u n+1 -A(θ n+1 )u n+1 -f B (θ n+1 ) e y = 2 u n δt - u n-1 2δt , (19)
3 2

θ n+1 + S(θ n+1 ) δt + ∇ • u n+1 θ n+1 -∇ • K ReP r ∇θ n+1 = 2 θ n + S(θ n ) δt - θ n-1 + S(θ n-1 ) 2δt . ( 20
)
This system of non-linear equations is solved at time t n+1 = (n + 1)δt, using two previous states: t n and t n-1 . We denote by f B (θ) the Boussinesq force that can be non-linear in the general case (e. g. natural convection or freezing of water).

The space discretization of variables over the domain Ω ⊂ R 2 uses a finite-element method based on a weak formulation of the system ( 18)-( 20). 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 = L 2 0 (Ω) = q ∈ L 2 (Ω) Ω q = 0 . ( 21
)
Following the generalization of the Stokes problem [START_REF] Temam | Navier-Stokes equations. Theory and numerical 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 ( 18)-( 20) 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 (22) 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 = 2 δt (u n , v) - 1 2δt u n-1 , v , ∀ v ∈ V (23) 3 2δt θ n+1 + S(θ n+1 ), φ + u n+1 • ∇θ n+1 , φ + K ReP r ∇θ n+1 , ∇φ = 2 δt (θ n + S(θ n ), φ) - 1 2δt θ n-1 + S(θ n-1 ), φ , ∀ φ ∈ V , ( 24
)
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 Ω ∂juj • ∂jvi, b : V × Q → R, b(u, q) = - Ω ∇ • u q = - 2 i=1 Ω ∂iui • q, c : V × V × V → R, c(w; z, v) = Ω [(w • ∇) z] • v = 2 i,j=1 Ω wj(∂jzi)vi.
Note that we introduced in Eq. ( 22) a penalty term on the pressure. Its role is analysed below at different levels of the implementation of the numerical method. At the discrete level of model equations, we infer from the incompressibility constraint ( 18) that the pressure is defined up to an additive constant. Requesting that p n+1 ∈ Q, i. e. the pressure field is of zero average in Ω, removes this uncertainty. For the discretized equations, the zero average constraint is easy to apply in spectral methods by, for instance, cancelling the first Fourier mode, corresponding to the average value [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]. This task becomes more difficult when local approximation numerical methods (as finite differences or finite elements) are used. In these cases, a usual technique is to explicitly prescribe the value of the pressure at a given discrete point. This implies to explicitly modify the matrix of the final linear system. Equation ( 22) has exactly this role. It is theoretically justified by the fact that the range space of ∇ • u is contained in Q = L 2 0 (Ω), which can be identified isometrically with L 2 (Ω)/R [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]. At a strictly algebraic level, the penalty term in ( 22) is introduced to avoid a zero lower diagonal block in the final global matrix. Even though the non-penalized problem is well-posed for finite-elements satisfying the inf-sup condition (as the Taylor-Hood element used here), the zero entries on this (pressure) diagonal block would require pivoting strategies when a direct solver is used for the solution of the linear algebraic system. Using ( 22), this block is filled with the matrix (-γM ), where M corresponds to the mass matrix of the P 1 discretization used for the pressure. Efficient direct methods using LU-type solvers become then available to solve the linear system (UMFPACK solver in our case and [START_REF] Aldbaissy | A full discretisation of the time-dependent Boussinesq (buoyancy) model with nonlinear viscosity[END_REF], SuperLU in [START_REF] Woodfield | Stability and finite element approximation of phase change models for natural convection in porous media[END_REF]).

A last comment on the penalty term on the pressure in ( 22) concerns its stabilizing effects. Since the Taylor-Hood finite elements used in our simulations satisfy the inf-sup condition, this technique is assimilated to a stable penalty, or, in other words, just a computational trick to obtain a good solution (see Remark 4.6, page 86 in [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF] (see Example 6.3.5, page 372 in [START_REF] Boffi | Mixed finite element methods and applications[END_REF]. The value of the penalty parameter γ could be arbitrarily small, without affecting the overall convergence rate. It is interesting to note that this technique could be used to stabilize finite element formulations that do not satisfy the inf-sup condition (see Remark 4.3, page 67 in [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]. In this case, it is assimilated to a brute-force penalty (see also Example 6.3.5 in [START_REF] Boffi | Mixed finite element methods and applications[END_REF] and the value of the penalty parameter γ has to be carefully chosen to avoid altering the convergence of the method: it depends on the grid size and the theoretical convergence rate of the used finite elements [START_REF] Boffi | Mixed finite element methods and applications[END_REF], and also, when applied to Stokes or Navier-Stokes equations, on the viscosity [START_REF] Oden | Penalty-finite element methods for the analysis of stokesian flows[END_REF].

As a final remark, we should mention that the large body of literature on penalty methods for Stokes or Navier-Stokes equations originates from the slightly compressible fluid model introduced by [START_REF] Temam | Navier-Stokes equations. Theory and numerical analysis[END_REF]. In this model, Eq. ( 22) is used to eliminate the pressure, by assuming that p n+1 = (-1/γ)∇ • u n+1 , and finally obtain a momentum equation ( 23) with a regularization term (1/γ) ∇ • v, ∇ • u n+1 replacing the pressure term. The equivalence between the two formulations was studied by [START_REF] Bercovier | Perturbation of mixed variational problems. application to mixed finite element methods[END_REF] (see also [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]. In our calculations, we use in Eq. ( 22) very low values of the penalty parameter (γ = 10 -7 ). This ensures very low values of the average on Ω for the pressure and also for the divergence of the velocity field, e. g. Ω p of order of 10 -11 and Ω ∇ • u n+1 of order of 10 -18 for the case of the natural convection of air.

The system of non-linear equations ( 22)-( 24) 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: b (u k+1 , q) -γ(p k+1 , q) = 0, (25) 3 2δt (u k+1 , v) + c (u k+1 ; u k , v) + c (u k ; u k+1 , v) + 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 ) - dfB dθ (θ k ) θ k+1 ey, v = 1 δt 2u n - 1 2 u n-1 , v +c (u k ; u k , v) - dA dθ (θ k ) θ k u k , v - dfB dθ (θ k ) θ k -fB(θ k ) ey, v , (26) 3 2δt θ k+1 + dS dθ (θ k ) θ k+1 , φ + (u k • ∇θ k+1 , φ) + (u k+1 • ∇θ k , φ) + K ReP r ∇θ k+1 , ∇φ = 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 ), φ . ( 27
)
Note that the last term of Eq. ( 26) cancels in the case of a linear Boussinesq force f B (see Eq. ( 10)); 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 equations ( 25)-( 27) 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. The Newton loop (following k) has to be 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 (25)-(27) to get w k+1 stop when w k+1 -w k < ξ N actualize (u n+1 , p n+1 , θ n+1 ) = w k+1 . ( 28
)
The FreeFem++ 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 ( 25)-( 27). 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.

Mesh adaptivity

We use the standard function (adaptmesh) which is a very convenient tool offered by FreeFem++ to efficiently adapt 2D meshes by metrics control [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]; [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), ( 29
)
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 trianguler 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. ( 30
)
Note that the FreeFem++ function buildmesh efficiently implements the Delaunay algorithm for generating triangular meshes. The mesh generator offered by FreeFem++ is generally sufficient for 2D applications, even when domains with complicated shapes are considered, but meshes may be imported in FreeFem++ from other software, e. g. Gmsh.

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 algorithm defined in [START_REF] Frey | Maillages. Hermès[END_REF]).

The possibility to use the standard FreeFem++ adaptmesh function was one of the main advantages that decided us in choosing this software to implement our algorithm to solve phase-change problems. No supplementary script layers or modifications of the standard function were necessary and the adaptivity part of the FreeFem++ scripts contains only a single instruction line. Recently, Zimmerman and Kowalski (2018) used the FEniCs finite element software offering a dual-weighted residual method for mesh adaptivity purposes. Using this method with an AMR (Adaptive Mesh Refinement) technique for phase change problems displayed a major drawback: the mesh was refined during the advancement of the liquid-solid front, but never coarsened behind. The metrics control used by the FreeFem++ adaptmesh function offers the ability to simultaneously refine the mesh in regions with strong gradients and coarsen the mesh in regions with low gradients. This function proved very efficient and versatile for phase-change problems by offering the possibility to take into account several metrics computed from different variables monitoring the time-evolution of the system. For natural convection systems, the mesh was adapted using the values of the two velocity components and the temperature. For phase-change systems, to accurately track the solid-liquid interface we added the enthalphy source term as a variable in the adaptivity criterion. For water systems (convection or freezing), we also added an extra function tracking the anomalous change of density around 4 o C. To monitor the time change of these variables used to compute the metrics intersection, we considered for each variable the metrics computed at actual (t n+1 ) and previous (t n ) time instants. This technique also reduces the impact of the interpolation on the global accuracy for time-depending problems (see also [START_REF] Belhamadia | Anisotropic mesh adaptation for the solution of the Stefan problem[END_REF]). 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 §6. It is important to emphasize that the standard mesh adaptivity and interpolation tools offered by FreeFem++ are very efficient. The CPU time needed to adapt the mesh, using up to six variables for the metrics evaluation, do not exceed 4% of the CPU time requested to advance the solution by a time step. The exact values of CPU times for the melting cases presented in this paper are reported in Tab. 3. Given this computational efficiency of the adaptivity tool, we could afford to refine the mesh every time step during the computation, even when complex domain shapes or physics were considered.

Description of the programs

In this section, we first describe the architecture of the programs and the organisation of provided files. Then we focus on the list of input parameters and the structure of output files.

Program architecture

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

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:

• test of the convergence of the method using manufactured solutions,

• 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 stationary 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 Figure 2. The obtained solutions are saved in the folder OUTPUT/Data. Depending on the output format selected by the user, data files are generated in specific folders for being visualized with: Tecplot, Paraview, Gnuplot or Medit. We also provide in the folder Figures ready-made layouts for these visualisation softwares. The user can thus obtain the figures from this paper using newly generated data. More details about the output structure are given below.

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 (6). 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, (31)
(2) :

V (2) ref = α H =⇒ t (2) ref = t (1) ref Pr =⇒ Re = 1/Pr, (32)
(3) :

V (3) ref = ν l H Ra Pr =⇒ t (3) ref = t (1) ref Pr Ra =⇒ Re = Ra Pr , (33)
• 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 ( 11) and ( 10),

• T hot , T cold : are dimensionless temperatures according to (6),

• 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 ( 14) used to regularize discontinuous functions are set by default as follows:

• 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]:

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 - - ρ(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 Figure 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 (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 Figure 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 (18). Default value = 10 -7 ,

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

• newtonMax: limits the maximum number of iterations in the Newton algorithm (28). 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 Figure 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 Figure 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. [START_REF] Roache | Verification and Validation in Computational Science and Engineering[END_REF]) 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 equations ( 7)-( 9) for a stationary case (Burggraf flow) and a time-dependent one [START_REF] Nourgaliev | Fully-implicit 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 , (34)
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 P r = 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:

u 1 (x, y) = σg (x)h (y), ( 35
)
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),
with σ > 0 a scaling parameter and functions

g(x) = x 5 5 - x 4 2 + x 3 3 , ( 36
)
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, ( 37
)
which ensures the continuity of the velocity at corners (u(0, 1) = u(1, 1) = 0), since non-slip walls are imposed for the other borders: u(x, 0) = u(0, y) = u(1, y) = 0. 20

The forcing terms that have to be added to the momentum and energy (temperature) equation are derived by injecting the exact solution (35) into the system ( 7)-( 9): 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 3a and3b illustrate the streamlines and the temperature field, respectively. Figure 4 plots the discretization error ε as a function of the grid size h = δx = δy. The grid was not adapted for these computations. Errors are measured in natural norms: H 1 -norm for the velocity components u 1 , u 2 and temperature θ, and L 2 -norm for the pressure p. Figures 4a and4b correspond to the simulation using P 1 finite elements for the temperature, while in Figs. 4c and4d P 2 finite elements were used for the temperature. Convergence rates computed by a lest-squares fit of the slopes of the curves are displayed for reference. The optimal order of accuracy (O(h 2 )) corresponding to the Taylor-Hood discretization of the flow field is obtained (Figs. 4a and4c), independently of the temperature discretization. Expected orders of accuracy for the temperature are confirmed in Figs. 4b (first order) and 4d (second order) for the two implementations. In the subsequent simulations, we use the P 2 discretization of the temperature, offering an optimal coupling between momentum and energy equations (note in Figs. 4 the lower values of the errors obtained with this discretization). Similar discretizations using Taylor-Hood finite elements for the fluid flow and P 2 for the temperature were used in recent contributions by [START_REF] Woodfield | Stability and finite element approximation of phase change models for natural convection in porous media[END_REF]; [START_REF] Belhamadia | A two-dimensional adaptive remeshing method for solving melting and solidification problems with convection[END_REF]. The global error ε (34) was computed using natural norms: L 2 -norm for the pressure p and H 1 -norm for velocity components u 1 , u 2 and temperature θ. Taylor-Hood P 2 /P 1 finite elements are used for the velocity and pressure. Discretization of the temperature using P 1 (plots a and b) and P 2 (plots c and d) finite elements.

f u1 = 0, ( 38
)
f u2 = σ 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 P rRe 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 ReP r (a (x)b(y) + a(x)b (y)) .

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 | Fully-implicit 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), ( 39
)
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 1. The corresponding forcing source

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

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 u1 = α u cos(t) cos(a) sin(b) -U c γ 1 sin(a) sin(b) + U c γ 2 cos(a) cos(b) (40) -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 u2 = -α u cos(t)
Guided by the results obtained in §5.1 for the space accuracy, we fixed the grid size to h = δx = 0.01 to ensure small spatial discretization errors and used P 2 finite elements for the temperature. 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 (34) was computed. The time convergence is displayed in Fig. 5a for the temperature and in Fig. 5b for the total velocity u = u 2 1 + u 2 2 . The expected second order convergence in time is obtained for all variables of the problem. Figure 5: Time accuracy of the numerical scheme tested using the time-dependent manufactured solution of [START_REF] Nourgaliev | Fully-implicit orthogonal reconstructed discontinuous galerkin method for fluid dynamics with phase change[END_REF]. Global error ε (34) computed at tmax = π for the temperature θ (a) and total velocity u = u 2 1 + u 2 2 (b) using L 2 and infinity norms. Discretization using P 2 finite elements for the temperature.

Numerical simulations of natural convection and phase-change problems

In this section we test the robustness of the distributed toolbox. We consider well defined benchmarks used to validate numerical codes for natural convection and phasechange problems. The difficulty of the computed cases is increased progressively by considering the following physical systems: (i) natural convection of air ( §6.1), (ii) melting of a phase-change material ( §6.2), (iii) alternate melting and solidification of a phasechange system ( §6.3) and (iv) the convection and the freezing of pure water ( §6.4). This approach allows us to test the programs by adding progressively non-linearities in the Newton algorithm. For each test case, we compare our results with experimental and previously published numerical data.

Natural convection of air

We start by testing the Newton algorithm ( 25)-( 27) for the case of natural convection, i. e. C = K = 1, A(θ) = S(θ) = 0. We consider the classical problem of the thermally driven square cavity [0, 1] × [0, 1], filled with air. The Boussinesq term f B (θ) is then linear and takes the form (10). Top and bottom walls are adiabatic, while the temperature is fixed on the left (hot) wall and the right (cold) wall. Natural convection flows are computed for three values of the Rayleigh number: Ra = 10 4 , 10 5 , 10 6 . The Prandtl number is set to P r = 0.71. It was shown in Le [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high Rayleigh number[END_REF] that the flow in this configuration becomes unsteady for Ra = 10 8.5 . Therefore, steady states could be computed for the chosen values of the Rayleigh number.

We provide programs for both steady (time-independent) and time-dependent cases. The steady case is performed using a continuation following the boundary value of the temperature: a small value for θ c (or θ h ) is set initially and is then smoothly increased until reaching the wanted value. The time-dependent case is computed until a steady state with a single convection cell is reached. To make the simulation more complicated, we present two configurations: i) the classical differentially heated square cavity and ii) a differentially heated cavity with an inner heated obstacle. All computations for the natural convection cases are performed using the stationary solver and mesh adaptivity. The initial mesh was generated using nbseg = 80 segments on each side, i. e. h min = 1/80 = 0.0125.

Classical differentially heated square cavity

The temperature is imposed at the left (hot) wall as θ h = 0.5 and at the right (cold) wall as θ c = -0.5. Top and bottom walls are adiabatic. The initial condition models a cavity filled with motionless air (u = 0), with a linear distribution of the temperature. Both steady and time-dependent codes converge to the same flow state with a single convection cell. this final state, horizontal u(y) and vertical v(x) velocity profiles were extracted at mid-domain (y = 0.5 and x = 0.5, respectively) and plotted in Figure 6. Our results are in very good agreement with reference numerical results obtained by Le [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high Rayleigh number[END_REF] with a spectral code. Table 2 offers a quantitative assessment of the accuracy of the present Newton method. The values of u max and its location Y are compared to reference values from Le [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high Rayleigh number[END_REF]. The Newton method gives results identical to reference values, with a difference less than 0.01%. Run u max at x=0.5 (error) Y (error) Reference values spectral 0.0648344 0.850 Newton (Steady) nbseg = 80 0.0648297 (0.007 %) 0.850394 (0.05 %) Newton (Unsteady) nbseg = 80 0.0648296 (0.007 %) 0.850532 (0.06 %)

Table 2: Natural convection of air in a differentially heated cavity. Maximum value umax of the horizontal velocity profile at mid-domain (x = 0.5) and location Y of this maximum. Comparison to reference values by Le [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high Rayleigh number[END_REF].

Differentially heated cavity with an inner heated obstacle

We consider the same differentially heated cavity as previously and add a centred square obstacle. The boundaries of the inner square are non-slip isothermal walls, maintained at a dimensionless hot temperature θ h = 0.8. The solution computed for Ra = 10 6 and P r = 0.71 is compared with the results obtained by [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF], who used a 6th order finite-difference method with an immersed boundary method to model the obstacle. The temperature distribution in the cavity is shown in Figure 7a. The vertical velocity profile v(x) along the horizontal symmetry line (y = 0.5) is displayed in Figure 7b and shows a very good agreement with the numerical results reported by [START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF]. Results obtained using the present Newton method (red solid line), with an initial mesh resolution nbseg = 80. Comparison with the results obtained using a 6th order finite-difference method by Moglan (2013).

Melting of a phase-change material (PCM)

We continue our validation tests by considering the full system (25)-( 27) for the case of the melting of a phase-change material. Two new non-linearities are now present in the system: the Carman-Kozeny penalty term A(θ) and the enthalpy source term S(θ). The function S is regularized using (15). We also consider that the material properties in the liquid and solid are the same, i. e. C = K = 1. This is a frequent assumption (Wang et al., 2010a;[START_REF] Ma | Solid velocity correction schemes for a temperature transforming model for convection phase change[END_REF]. Five cases were computed (the exact values of the defining parameters are summarized in Table 3):

• PCM-Case #1 simulates the experimental study of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]. It consists of a differentially heated square cavity, filled with octadecane paraffin.

• PCM-Case #2 is extracted from the collective publication by [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF], in which the results of different numerical approaches were compared for the simulation of the melting of a PCM.

• PCM-Case #3 simulates the melting of a cylindrical PCM with heated inner tubes, as in [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF].

• PCM-Case #4 simulates the melting of Gallium in a rectangular cavity heated by the side-wall, as in [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF].

• PCM-Case #5 reproduces the simulation of [START_REF] Nourgaliev | Fully-implicit orthogonal reconstructed discontinuous galerkin method for fluid dynamics with phase change[END_REF] To guide the user of the toolbox, we also provide in Table 3 the values of the CPU time necessary to run each case (with default parameters) on a personal computer and the typical number of triangles of the generated adaptive mesh.

Supplemental Material for animations depicting the dynamics of the cases presented in this section are provided at http://lmrs-num.math.cnrs.fr/2019CPCP1.html. 6.2.1. PCM-Case #1: Melting of an octadecane PCM in a square cavity [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] studied experimentally the melting of an octadecane PCM in a square cavity of height H = 1.5 cm. His results were often used to validate numerical methods [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]Wang et al., 2010a;[START_REF] Ma | Solid velocity correction schemes for a temperature transforming model for convection phase change[END_REF][START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]. The material is initially solid (θ 0 = -0.01) and melts progressively starting from the left boundary, maintained at a hot temperature θ h = 1. The right boundary is also isothermal, with cold temperature θ c = -0.01. Horizontal boundaries are adiabatic. The other parameters of this case are reported in Table 3.

The computation starts from a refined mesh near the hot boundary. Mesh adaptivity is applied at each time step using metrics computed from three variables: the two fluid velocities and the enthalpy source term S. To reduce the impact of the interpolation on the global accuracy, we use two successive fields (S n ) and (S n+1 ) in the adaptivity procedure. This allows us to refine the mesh in the fluid part of the domain and inside the artificial mushy. Figure 8a gives an illustration of the adapted mesh at dimensionless time t = 78.7. In Figure 8b we compare the position of the solid-liquid interface at t = 39.9 and t = 78.7 with the experimental data of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and previously published numerical results (Wang et al., 2010a;[START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]. The obtained shape and position of the liquid-solid interface is closer to experimental results than numerical results reported in Wang et al. (2010a). This is a direct consequence of the mesh adaptivity capabilities of our method. [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and numerical results of [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] and Wang et al. (2010a) for two time instants, t = 39.9 and t = 78.7.

This comparison also allowed us to finely tune the value of the constants used in the model ( 16). 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 [START_REF] Kheirabadi | The effect of the mushy-zone constant on simulated phase change heat transfer[END_REF][START_REF] Mathura | Influence of mushy zone constant on thermohydraulics of a PCM[END_REF]. This choice of the value of this constant is a still open problem. 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 varying in the range [10 6 , 10 8 ]. Imposing a too large value (C CK = 10 10 ) results in artificially slowing the propagation of the melting front. Consequently, we set for all subsequent simulations C CK = 10 6 .

PCM-Case #2: Melting of an octadecane PCM with high Rayleigh number

This case considers the same problem of the melting of a PCM, but with a very high value of the Rayleigh number Ra = 10 8 (see Table 3). This case is very challenging since the natural convection becomes important in the fluid flow, and enhances considerably the heat transfer. [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). These results 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 have 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 behaviours in Figures 9a and9b: Lacroix and Delannoy seem to be insufficiently converged as shown by Figure 9a, and Binet-Lacroix overestimates the average Nusselt number by more than 30% (Figure 9b). Hence, this collection of results allows us to validate our numerical method and check if realistic results are obtained for complex physical configurations. For comparison purpose, we extract from simulations the position of the melting front and the Nusselt number Nu at the left wall (x = 0) for each of the five methods presented by [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF]. The Nusselt number Nu is defined as follows:

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

(41) The position of the melting front for three time instants, t = 1, 3 and 5 is reported in Figure 9a. Our results are for each case in fairly good agreement with those of Gobin and Le Quéré. Details of their numerical method are presented in [START_REF] Gobin | Melting from an isothermal vertical wall. synthesis of numerical comparison exercise[END_REF]. Gobin used a front-tracking method based on a coordinate transformation with a finite volume method in a 62 × 42 grid and Le Quéré solved a single domain method using a second-order scheme with a finite volume method and a 192 × 192 grid. The time evolution of the Nusselt number is presented in Figure 9b. A very good agreement is obtained with the results of Gobin and Le Quéré. A relative difference, less than 2%, is observed for the Nusselt number.

The high Rayleigh number Ra = 10 8 is a very demanding numerical test. The high velocity, inducing a very narrow thermal boundary layer, can lead to unrealistic results if under-resolved. This explains why some numerical methods failed in [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF]. The interest of the mesh adaptation is clearly demonstrated for this case, since we initially used a coarse 40 × 40 grid.

PCM-Case #3: Melting of cylindrical PCM with inner heated tubes

Previous cases #1 and #2 considered phase change problems evolving in a simple geometry, a square cavity. 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 configurations with one, four and nine heated tubes. The size of the tubes is adjusted to have the same total tube area for all configurations. The radius R i of the inner tube is R i = R/4 for the case with one tube, R i = R/8 for the four heated tubes case and R i = R/12 for the case with nine tubes. A Dirichlet boundary condition (θ = θ h ) is applied to the boundary of inner tubes. A Neumann boundary condition (∂θ/∂n = 0) is used for the outer boundary. For the velocity, all boundaries are considered as non-slip walls (u = 0). Only half of the domain is simulated since all configurations are symmetric with respect to the vertical axis (see Figure 10). The mesh is refined initially around the inner tubes, and is dynamically adapted at each time step around the melting front and the thermal boundary layer area. The same metrics presented in Section 6.2.1 are used for the mesh adaptivity.

Figure 10 shows the temperature field and the position of the solid-liquid interface (black line) for 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. The mesh is nicely adapted following the evolution of the melting interface, even after its separation in several distinct fronts touching the outer boundary (see Figures 10b,c).

To estimate the efficiency of each configuration, we plot in Figure 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 configurations melts 5 times faster than the reference configuration with one tube. Note also from Figure 11 a good agreement between our results and those reported by [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF] for the evolution of the liquid fraction.

PCM-Case #4: Melting of Gallium in a rectangular cavity

The melting of the Gallium in a rectangular cavity 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. Experimental results exhibited indeed a mono-cellular structure, while many researchers claimed that this observation was incorrect. Prior to Dantzig (1989), both experimental and numerical studies reported a single convection cell in the fluid phase. Later, simulations provided solutions with multi-cellular flow. [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF] concluded that the mono-cellular observation was caused by a problem of convergence of the numerical solution, due to coarse grids or inconsistencies in the mathematical model. Therefore, this test case simulating the melting of the Gallium is a relevant exercise to test the accuracy of our method. The parameters of this case are reported in Table 3. To capture the very small convection cells 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. With our adaptive method, a maximum number of 4820 triangles is necessary 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]. The grid size is thus reduced with our method by a factor of 100. The time evolution of the flow is presented in the Figure 12. Temperature field, streamlines and position of the melting front are plotted for several time instants: t = 0.0015, 0.006, 0.01, and 0.019. These values were chosen to visualise 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 (Figure 12). The number of cells decreases later through a process of roll merging, as it was also reported by [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). 6.2.5. PCM-Case #5: Solid crust formation in a highly distorted mesh [START_REF] Nourgaliev | Fully-implicit orthogonal reconstructed discontinuous galerkin method for fluid dynamics with phase change[END_REF] used a discontinuous Galerkin finite element method to simulate the solid crust formation inside a highly distorted domain (Figure 13). The fluid is initially motionless with an initial dimensionless temperature θ 0 = 2. The temperature of fusion is set to θ f = 1.4, according to Nourgaliev et al. ( 2016) parameters (see Table 3 for the values of all parameters). The left side boundary is maintained at a cold temperature θ c = 1.39 in the initial stage. The right wall is isothermal, with hot temperature θ h = 2. A nearly steady-state natural circulation is induced in the early time evolution of the flow. Then, the cold temperature at the left wall is decreased to θ c = 1, below the temperature of solidification. At this point, the formation of a solid crust layer starts at the left boundary. Figure 13 shows the temperature field and the streamlines for the time instant t = 30. Our results are in a very good agreement with those of [START_REF] Nourgaliev | Fully-implicit orthogonal reconstructed discontinuous galerkin method for fluid dynamics with phase change[END_REF].

Melting-solidification cycle of a PCM

We address in this section the challenging problem of simulating a complete meltingsolidification cycle of a PCM. The simulation starts from the final state obtained in §6.2.1. At t = 78.7, the PCM is partially melted, with the liquid fraction L f = 0.5 (see Figure 14a). The temperature of the left boundary is then dropped to a cold temperature θ c = -0.01 < θ f , identical to that of the right wall. The solidification starts and the solid phase propagates into the cavity from both left and right sides. Panels (a) to (e) in Figure 14 depict the time evolution of the solidification process. The solid phase is represented in blue and corresponds to the region of temperature θ ≤ θ f = 0, while the solid-liquid interfaces (i. e. contour lines of θ = 0) are represented by bold solid lines. The solidification is completely achieved at t = 4600. Note that the mesh adaptivity of our toolbox is able to accurately track the two solidification fronts. The adapted mesh shown in Figure 14f illustrates that a finer mesh is well generated along the interface (θ = 0), while a coarser mesh is obtained in the regions of low gradients.

Natural convection of water and water freezing

We consider in this section the natural convection and the solidification of water. With these cases, we test the ability of our numerical system to deal with additional non-linear terms. Since pure water exhibits non-linear density variation for T < 10.2 • C, with a maximum at T m = 4.0293 • C, the Boussinesq force becomes non-linear. We used the following density-temperature relationship suggested by [START_REF] Gebhart | A new density relation for pure and saline water[END_REF]:

ρ(T ) = ρ m (1 -w |T -T m | q ) , (42)
with ρ m = 999.972 [kg/m 3 ], w = 9.2793 • 10 -6 [( • C) -q ], and q = 1.894816. Hence, the bouyancy term f B = g(ρ ref -ρ)/ρ ref in ( 7) is not any more linear and becomes after scaling:

f B (θ) = Ra P r Re 2 1 βδT ρ(θ f ) -ρ(θ) ρ(θ f ) , ( 43
)
where β = (1/ρ m ) (dρ/dT ) is the thermal expansion coefficient taking the value β = 6.91 • 10 -5 [(K) -1 ] ( [START_REF] Scanlon | A numerical analysis of buoyancy-driven melting and freezing[END_REF]. We simulate a differentially heated square cavity 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 non-dimensional parameters describing the problem are (see [START_REF] Michalek | Simulations of the water freezing process -numerical benchmarks[END_REF] for physical details): Ra = 2.518084 • 10 6 , Pr = 6.99 and Ste = 0.125.

Natural convection of water

The initial temperature is linearly distributed in the square cavity, with a hot temperature T h = 10 • C at the left wall and a cold temperature T c = T f = 0 • C at the right wall. The temperature field and the streamlines of the steady state are presented in Figure 15a. The isoline θ = θ m , corresponding to the line of maximum density is represented by a dashed line. Due to the anomalous thermal variation of water density, 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 ).

A more precise comparison with previously published results is shown in Figure 15b. The obtained temperature profile θ(x) along the horizontal symmetry line of the cavity (y = 0.5) is in good agreement with the numerical results of [START_REF] Michalek | Simulations of the water freezing process -numerical benchmarks[END_REF]. Their results were obtained with finite-volume and finite-difference codes (FLUENT and FRECONV3V). 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 1,807 vertices (3,430 triangles).

Four test cases were presented, by adding progressively non-linearities in the system of equations: (i) natural convection of air in a differentially heated cavity, (ii) melting of a PCM, (iii) melting-solidification cycle of a PCM, (iv) natural convection and freezing of water. The computations for case (ii) were rendered more challenging by considering complex geometries (highly distorted mesh, cylindrical PCM with inner heated tubes) and computationally demanding cases (high Rayleigh numbers). The efficiency of the adaptivity method by metric control was investigated by tracking simultaneously several interfaces (two melting fronts during the solidification cycle and density inversion interface for water flows).

For each test case, we provided a separate folder containing all the necessary files (parameters, restart files) necessary to run them directly. We described in the text body of the paper the expected results and their validation. A very good agreement with experimental data or numerical results was obtained for all considered test cases, proving the capability of our method to tackle a large range of problems. Ready-made scripts and layouts are provided with the toolbox to allow the user to generate the figures presented in this paper with newly generated data after running the programs. Validation data sets from experiments or previous publications are included in these layouts. Movies depicting the dynamics of some cases simulated in this paper are provided as Supplemental Material at http://lmrs-num.math.cnrs.fr/2019CPCP1.html.

Since FreeFem++ is a free software, the method could be easily implemented and tested by anyone interested in simulating phase-change problems. All technical issues related to the implementation of the finite element method are hidden, allowing to focus on numerical algorithms and their performance. This offers the possibility to address other computational challenges related to different physical or mathematical models in this field.

The extension of the method for 2D and 3D cases, using domain decomposition methods adapted to parallel computing, will be reported in a forthcoming contribution.

Figure 1 :

 1 Figure 1: Folder tree structure of the PCM-Toolbox-2D to solve solid-liquid phase-change problems.

Figure 2 :

 2 Figures

Figure 3 :

 3 Figure 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 4 :

 4 Figure4: Space accuracy of the numerical scheme tested using the Burggraf manufactured solution. The global error ε (34) was computed using natural norms: L 2 -norm for the pressure p and H 1 -norm for velocity components u 1 , u 2 and temperature θ. Taylor-Hood P 2 /P 1 finite elements are used for the velocity and pressure. Discretization of the temperature using P 1 (plots a and b) and P 2 (plots c and d) finite elements.

  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 P rRe 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, y, t) cos(a) cos(b) + 2K ReP r T c cos(a) sin(b),

Figure 6 :

 6 Figure 6: Natural convection of air in a differentially heated cavity for values of Ra ranging from 10 4 to 10 6 . (a) Vertical velocity profile v(x) along the horizontal symmetry line (y = 0.5). (b) Longitudinal velocity profile u(y) along the vertical symmetry line (x = 0.5). Results obtained using the present Newton method (symbols), with an initial mesh resolution nbseg = 80. Comparison with the spectral simulations by Le Quéré (1991) (solid lines).

Figure 7 :

 7 Figure 7: Natural convection of air in a differentially heated cavity with an inner heated square for Ra = 10 6 . Temperature field (a) and vertical velocity profile along the horizontal symmetry line (b).Results obtained using the present Newton method (red solid line), with an initial mesh resolution nbseg = 80. Comparison with the results obtained using a 6th order finite-difference method by[START_REF] Moglan | Modeling and numerical simulation of flow and heat phenomena in a telecommunication outdoor cabinet[END_REF].

Figure 8 :

 8 Figure 8: PCM-Case #1. Melting of the PCM. (a) Adapted mesh at time instant t = 78.7. (b) Position of the solid-liquid interface. Comparison with experimental data of[START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and numerical results 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 t = 78.7.

Figure 9 :

 9 Figure 9: PCM-Case #2. Melting of the PCM with high Rayleigh number (Ra = 10 8 ). Comparison with five sets of results presented in Bertrand et al. (1999). (a) Location of the solid-liquid interface at dimensionless time instants t = 1, 3 and 5. (b) Temporal evolution of the Nusselt number.

Figure 10 Figure 11 :

 1011 Figure 10: PCM-Case #3. 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 localized with black lines (only half of the domain is simulated).

Figure 12 :

 12 Figure12: PCM-Case #4. Melting of Gallium: temperature field, streamlines, and melting front for dimensionless time instants (panels a to d): 0.0015, 0.006, 0.01, and 0.019. For a better view of the convection cells, a ratio 2:1 was used for the axis dimensions.

Figure 13 :

 13 Figure 13: PCM-Case #5. Solid crust formation in a distorted mesh. Temperature field and streamlines at dimensionless time t = 30: our simulation (a) and Nourgaliev et al. (2016) (b).

  using highly distorted meshes to compute a natural convection flow with solid crust formation.

		case #1	case #2 case #3 case #4 case #5
	Ra	3.27 • 10 5	10 8	5 • 10 4	7 • 10 5		10 6
	P r	56.2	50	0.2	0.0216		0.1
	Ste	0.045	0.1	0.02	0.046	4.854
	δt	0.1	10 -3	10 -3	10 -5	10 -3
	V ref	ν l H	ν l H	ν l H	ν l H	ν l H	Ra P r
	CPU time	01:09:05 18:17:56 00:49:21 08:32:58 02:49:20
	Number of triangles	2, 900	7, 000	2, 076	3, 600	2, 769
	Number of iterations	798	5, 000	491	2, 001	3, 027
	ratio CPU-time adapt/iteration	2.5%	2.5%	2%	2%		4%

Table 3 :

 3 Parameters for the cases simulating the melting of a phase-change material.

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

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

Acknowledgements

This project was co-financed by the European Union with the European regional development funds and by the Normandy Regional Council via the M2NUM (ERDF, HN0002081) and M2SiNUM (ERDF, 18P03390/18E01750/18P02733) projects. Part of this work was performed using computing resources of CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie, Normandy, France).

A comprehensive physical description of this case is given in our recent paper [START_REF] Rakotondrandisa | Numerical modelling of a melting-solidification cycle of a phase-change material with complete or partial melting[END_REF]. The case of the solidification process starting after a complete melting of the PCM is also presented in this paper. We can conclude from this section that the Newton method is able to deal with either melting or solidification process. The simulations are in good agreement with existing experimental and numerical results for the melting and show consistent behaviour for the solidification.

Water freezing

We finally consider the difficult case of water freezing in a square cavity. The initial state for this computation is the convection steady pattern in the cavity presented in Figure 15. The freezing starts by dropping smoothly the temperature of the cold right wall from

Figure 16(a) superimposes the experimental image from Kowalewski and Rebow (1999) with our numerical results for the same physical time t ϕ = 2340 [s]. The flow pattern in the liquid phase, shown in Figure 16(b) also corresponds very well qualitatively to the experimental image. The simulation was performed with small time steps (δt = 10 -2 ≈ 15[s]) and reasonable-size grids (less than 3000 nodes) due to the efficiency of the adaptivity algorithm.

The two recirculating zones being separated by the line T = T m we used the mesh adaptivity capability of the method to refine the grid along this line. The metrics used for adaptivity were computed from the two components of the velocity and a P 1 function φ(T ) "tracking" the value T m , defined by the general regularization expression:

with R φ = 0.02. To reduce the impact of the interpolation on the global accuracy (see also [START_REF] Belhamadia | Anisotropic mesh adaptation for the solution of the Stefan problem[END_REF]), we used both φ(T n ) and φ(T n+1 ) in the adaptivity procedure. The final mesh is displayed in Figure 16c, clearly showing that the mesh is refined along the lines T = T m and T = 0 • C. This allows to accurately capture the structure and the extent of the two recirculating zones, features that are difficult to obtain with fixed meshes (discrepancies in numerical results are 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]). As a consequence, the temperature contours lines in Figure 16d are smooth and clearly define the two interfaces in the system: the liquid-solid interface (θ = 0) and the density inversion interface (θ = 0.4) separating the two recirculating liquid regions. [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]; 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 • C) and also along the line of maximum water density (T = 4 • C) (d) temperature iso-lines.

Summary and conclusions

We provide with this paper an adaptive finite-element toolbox for solving twodimensional phase-change problems with convection. The programs were written using FreeFem++ , a free software offering a programming syntax close to the mathematical formulation. A single domain numerical approach was first derived. The details of the finite-element formulation were then presented. The key ingredients of the implemented method are: (i) a second order accuracy in space and time; (ii) the use of an adaptive finite element method with a well chosen regularization of the functions representing the variation of thermodynamic properties at the solid-liquid interface, and (iii) a fully implicit discretization with a Newton algorithm for solving the non-linear system of equations.