I Hamzaoui

S Millet

V Botton

A Benzaoui

D Henry

L Hachani

R Boussaa

K Zaidat

Y Fautrelle

A 2D 1/2 model for natural convection and solidification in a narrow enclosure

Efficient numerical models are derived for problems of natural convection and material solidification in a horizontal differentially heated slender cavity. These 2D 1/2 models are obtained by averaging the equations momentum, heat, and mass conservation along the transverse direction assuming both a constant temperature and a well defined velocity profile in this direction. Based on our former works, the transverse velocity profile is assumed to be either a Poiseuille profile (2D 1/2 P model), or Hartmann-type profiles featuring two boundary layers on the sides of a uniform bulk (2D 1/2

H model). For this 2D

1/2 H model, however, a parameter δ (giving the boundary layer thickness) has to be adjusted: optimal values have been found in a large range of the control parameters and expressed as a reliable fitted function of Gr. The ability of the model to reproduce 3D results in a 2D framework is investigated in a large range of the control parameters (Prandtl number P r and Grashof number Gr); the validity domain of the model in this parameter space is also clarified and rigorously defined. A good precision is obtained for natural convection problems (intensity of the flow, temperature field) as well as for solid-liquid phase change problems (shape, position, and evolution of the front). Yet, a comparison with unpublished experimental data of solidification of pure tin is particularly presented. For this, boundary conditions for the simulation are defined after a post-treatment of time-dependent experimental data. This ensures the implemented boundary conditions to be representative of the experimental process despite a significant and time dependent thermal resistance between the walls of the crucible and the liquid. A very good agreement is observed between the 2D

1/2 H model and the experimental measurements for this pure tin solidification experiment in the AFRODITE set-up.

Introduction

The need of efficient reliable numerical models to perform fast simulations and large parametric studies is prominent and more pronounced when the simulated problems are complex and/or computational time consuming. This is specially the case for the solidification of metallic liquids in convective conditions that can be seen as a coupled problem, with front tracking issues and unsteady and multiscale aspects. Reference experimental investigations of convective effects in solidification are very often made in very slender configurations. This is the case, in particular, for in situ X-ray imaging of solidification [START_REF] Reinhart | Direct simulation of a directional solidification experiment observed in situ and real-time using X-ray imaging[END_REF][START_REF] Reinhart | Influence of natural convection during upward directional solidification: A comparison between in situ X-ray radiography and direct simulation of the grain structure[END_REF][START_REF] Nagira | Impact of melt convection induced by ultrasonic wave on dendrite growth in Sn-Bi alloys[END_REF][START_REF] Shevchenko | The effect of natural and forced melt convection on dendritic solidification in Ga-In alloys[END_REF], which requires very thin samples. This is also the case of the directional solidification benchmark AFRODITE on which the present paper focuses. It was proposed as a 'quasi-2D' experiment by Fautrelle and coworkers [START_REF] Wang | A quasi two-dimensional benchmark experiment for the solidification of a tin-lead binary alloy[END_REF][START_REF] Wang | An investigation of the influence of natural convection on tin solidification using a quasi two-dimensional experimental benchmark[END_REF][START_REF] Hachani | Experimental analysis of the solidification of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Boussaa | 2D and 3D numerical modeling of solidification benchmark of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Hachani | Multiscale statistical analysis of the tin-lead alloy solidification process[END_REF][START_REF] Boussaa | Macrosegregations in Sn-3 wt%Pb alloy solidification: Experimental and 3D numerical simulation investigations[END_REF], where, indeed the solidified samples are 10 cm long, 6 cm high and only 1 cm thick. A very challenging issue when comparing numerical models to this type of experimental investigations is to be able to implement appropriate realistic boundary conditions for the heat transfer problem [START_REF] Haagen | A comprehensive benchmark of fixedgrid methods for the modeling of melting[END_REF][START_REF] Franquet | On the modeling impacts of the energy equation in the simulation of melting[END_REF]. The AFRODITE setup has been designed and continuously improved to overcome this issue, with significant efforts to impose clear boundary conditions and to monitor the temperature field in space and time thanks to a number of thermocouples.

The recent experimental work of Hachani et al. [START_REF] Hachani | Multiscale statistical analysis of the tin-lead alloy solidification process[END_REF] on binary alloy solidification shows that, though regular patterns are expected in the solidified ingot mesostructure, significant statistical variations are also to be expected in their location. This is a strong and new argument in favor of the development of fast efficient models to be used in numerical statistical approaches. AFRODITE has served as a reference experiment case for several numerical model validations [START_REF] Boussaa | 2D and 3D numerical modeling of solidification benchmark of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Boussaa | Macrosegregations in Sn-3 wt%Pb alloy solidification: Experimental and 3D numerical simulation investigations[END_REF][START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF][START_REF] Zheng | Use of a mixed columnarequiaxed solidification model to analyse the formation of as-cast structure and macrosegregation in a Sn-10 wt%Pb benchmark experiment[END_REF][START_REF] Carozzani | Direct simulation of a solidification benchmark experiment[END_REF][START_REF] Gandin | Direct modeling of structures and segregations up to industrial casting scales[END_REF][START_REF] Rad | Validation of a model for the columnar to equiaxed transition with melt convection[END_REF][START_REF] Zheng | Role of fragmentation in as-cast structure: numerical study and experimental validation[END_REF]. Boussaa et al. [START_REF] Boussaa | 2D and 3D numerical modeling of solidification benchmark of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Boussaa | Macrosegregations in Sn-3 wt%Pb alloy solidification: Experimental and 3D numerical simulation investigations[END_REF] used Sn-3wt % Pb alloy solidification results to validate a two-phase volume averaged technique based model, and pointed out that 2D calculations overestimate the velocity field by a 1.5 factor thus leading to inconsistencies with respect to the experimental results. On the contrary, they found that the damping effect of the lateral walls in 3D models leads to better agreement with the experiment as it conserved in particular the prediction of the main mesostructure in the solidified ingot. The same observation is made by Carozzani et al. [START_REF] Carozzani | Direct simulation of a solidification benchmark experiment[END_REF] in their 3D CAFE model validation: since the 2D model did not take into consideration side walls effects, it overestimated convective effects with respect to experiment. This led the authors to consider the 2D model to be limited and only implement their grain-structure prediction model into 3D simulations.

Fully 3D calculations are found to be reliable but still expensive, owing, in particular, to the solid-liquid front tracking requirement in such solidification problems. Botton et al. [START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF] pointed out that considering the temperature field to be 2D is a good approximation, since the cavity is narrow, the liquid metals are good thermal conductors (very low Prandtl number), and the heat flux is along the horizontal longitudinal direction. In contrast, the velocity field cannot be considered as 2D because of its variations along the transverse direction due to the no-slip condition at the side walls. To take into account the three-dimensionality of the velocity field, an alternative to fully 3D simulations can be obtained by integrating the equations of motion along the transverse direction assuming given shapes for the velocity, temperature, and pressure profiles. Note that a trustful representation of the temperature field is one of the main keys for solidification problems modeling. We consider that a very first compulsory step towards the implementation of such an efficient binary alloy solidification model is to validate its ability to: [START_REF] Reinhart | Direct simulation of a directional solidification experiment observed in situ and real-time using X-ray imaging[END_REF] reproduce the temperature field in natural convection configurations; (2) predict the shape and position of the solid-liquid interface in pure metal solidification cases.

Historically speaking, the use of slender experimental configurations was highlighted by Hele-Shaw [START_REF] Hele-Shaw | On the motion of viscous fluid between two parallel plates[END_REF] to get rid of inertia effects. By integrating the inertialess equations of motion along the transverse direction, he showed that the flow acts as a potential flow. Schlichting et al. [START_REF] Schlichting | Boundary layer theory[END_REF] give a more extensive description of this approach and show that inertia terms can be added to the formulation. This approach has since been largely used in experimental, theoretical and numerical studies in a wide range of fields. The possibility to apply such a method to solidification problems is indicated by Ockendon et al. [START_REF] Ockendon | Viscous Flow[END_REF]. Katz et al. in further studies, [START_REF] Katz | Simulation of directional solidification, thermochemical convection, and chimney formation in a Hele-Shaw cell[END_REF] used the Hele-Shaw approach to model the directional solidification experiment of aqueous ammonium chloride and pointed out its ability to reproduce chimney formation observed in the experiment. In instability studies, Gondret et al. [START_REF] Gondret | Shear instability of two-fluid parallel flow in a Hele-Shaw cell[END_REF] and Ruyer-Quil [START_REF] Ruyer-Quil | Inertial corrections to the Darcy law in a Hele-Shaw cell[END_REF] used it to study the shear instability of two immiscible fluids, a viscous liquid and gas. The model was also used in flows presenting interfaces: Bensimon et al. [START_REF] Bensimon | Viscous flows in two dimensions[END_REF] and Homsy et al. [START_REF] Homsy | Viscous fingering in porous media[END_REF] studied the viscous fingering, Maxworthy [START_REF] Maxworthy | The nonlinear growth of a gravitationally unstable interface in a Hele-Shaw[END_REF] studied the nonlinear growth of an unstable interface in a large Hele-shaw cell and Howison [START_REF] Howison | Cusp development in Hele-Shaw flow with a free surface[END_REF] studied the cusp development in Hele-Shaw flow with a free surface; in complex fluid flow: Savage et al. [START_REF] Savage | The motion of a finite mass of granular material down a rough incline[END_REF] studied the motion of a granular material over a rough incline; and in continuous casting by Deam et al. [START_REF] Deam | Instability in filling thin walled continuous castings[END_REF]. It has also been used in coupled problems such as chemical reaction-diffusion phenomena coupled with hydrodynamic flow problems: for example Huang et al. [START_REF] Huang | Pattern formation and evolution near autocatalytic reaction fronts in a narrow vertical slab[END_REF] and Bockmann et al. [START_REF] Böckmann | Growth rates of the buoyancy-driven instability of an autocatalytic reaction front in a narrow cell[END_REF] studied the formation and growth of autocatalytic reaction front. Bockmann et al. [START_REF] Böckmann | Growth rates of the buoyancy-driven instability of an autocatalytic reaction front in a narrow cell[END_REF] found that the agreement between the numerical and the experimental results, for the autocatalytic reaction front growth, is far better with a 2D 1/2 model (appellation used from now on for models including inertia terms in the averaged equations) than with the Hele-Shaw model (i.e. an inertialess model), even though some discrepancies persist. Similar methods are also widely used in free surface hydraulics [START_REF] Chow | Open-channel Hydraulics[END_REF] under the denomination of shallow-water equations. In most of the previously mentioned works, the parabolic Poiseuille profile is used when integrating the equations of motion along the transverse direction, which is well justified in pressure driven duct flows at steady state. Botton et al. [START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF] pointed out that this assumption is questionable in buoyancy driven flows in closed cavity. They showed that the use of a transverse profile for the velocity featuring a uniform bulk and two boundary layers (called Hartmann type profile) rather than the Poiseuille profile when averaging the 3D equations led to a better consideration of the side walls effect and an improvement of the simulation results, i.e. better agreement with the full 3D results. In their Hartmann type profile formulation appears a parameter δ, assimilated to a boundary layer thickness. They suggested that this parameter can be tuned to an optimal value, but only two practical cases were considered. Here we investigate the existence and the evolution of this optimal value of δ in natural convection problems over a large range of the control parameters covering many practical situations. The validity domain of the model in the (Grashof, Prandtl) parameter space is clarified and rigorously defined based on a strict criteria. We also present quantitative comparisons between the results obtained with the 3D, 2D and 2D 1/2 models in these natural convection situations.

As our main interest is in solidification problems, two solid-liquid phase change situations in differentially heated cavities are then considered and discussed. In both of them, the solid-liquid front tracking is based on a multi-domain approach as presented by Wolf et al. [START_REF] Wolf | Solidification of a pure metal at a vertical wall in the presence of liquid superheat[END_REF] and more recently used by Avnaim et al. [START_REF] Avnaim | Comparison of threedimensional multidomain and single-domain models for the horizontal solidification problem[END_REF][START_REF] Avnaim | Numerical and experimental study of the traveling magnetic field effect on the horizontal solidification in a rectangular cavity. Part 2: Acting forces ratio and solidification parameters[END_REF]. In the first situations, the hot and cold walls temperatures are above and below the melting temperature of the considered pure material, respectively. The initial liquid domain at rest will evolve towards a steady solution corresponding to an equilibrium between the solid and the liquid phases. The different solid-liquid front positions and shapes, obtained with the 3D, 2D and 2D 1/2 models are reported and compared. The use of the 2D 1/2 models is shown to lead to a better estimate of the 3D front shape and position than that obtained with the 2D model. This is attributed to the fact that the 2D 1/2 models consider the damping effects at the side walls. The other situation is an experimental solidification process performed within the 'AFRODITE' setup: the hot and cold walls temperatures are decreased at a same rate until they are both below the melting temperature and, starting from a pure liquid tin sample, a tin ingot is eventually obtained. The original, yet unpublished results of the experiment are reported. A post processing of the experimental thermocouples data is used to define the effective time dependent boundary conditions for the 2D 1/2 model. Results from the 2D 1/2 model are presented and compared to the experimental results, illustrating the ability of the 2D 1/2 model to simulate such a solidification case.

Geometrical and thermal configuration

The geometrical configuration considered in this study is presented in Fig. 1. It is a differentially heated parallelepipedic cavity. The dimensions are chosen in order to correspond to the crucible of the benchmark experiment proposed by Fautrelle and co-workers [START_REF] Wang | A quasi two-dimensional benchmark experiment for the solidification of a tin-lead binary alloy[END_REF][START_REF] Wang | An investigation of the influence of natural convection on tin solidification using a quasi two-dimensional experimental benchmark[END_REF][START_REF] Hachani | Experimental analysis of the solidification of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Boussaa | 2D and 3D numerical modeling of solidification benchmark of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Hachani | Multiscale statistical analysis of the tin-lead alloy solidification process[END_REF][START_REF] Boussaa | Macrosegregations in Sn-3 wt%Pb alloy solidification: Experimental and 3D numerical simulation investigations[END_REF]. This cavity is confined in the transverse x 3 direction. Its aspect ratios are: A w = W/H = 1/6 and A l = L/H = 10/6, where L is the length in the x 1 direction (L = 10 cm), H is the height in the x 2 direction (H = 6 cm) and W is the width in the x 3 direction (W = 1 cm). The top, bottom and side walls are adiabatic, while the two end walls at x 1 = 0 cm and x 1 = 10 cm are at constant temperatures T * h and T * c , respectively. Low Prandtl number materials will be considered in the study. As a first step, the natural convection initiated in the cavity will be studied. Solid-liquid phase change cases will then be examined. The control parameters for the study of natural convection are the Prandtl number:

P r = ν α , (1)
and the Grashof number:

Gr = gβ(∆T * /L)H 4 ν 2 , ( 2
)
where ν, α, and β are the kinematic viscosity, the thermal diffusivity and the coefficient of thermal expansion, respectively, and ∆T * = T * h -T * c is the temperature difference between the hot and cold walls of the cavity. The 3D dimensionless governing equations, obtained using ν/H, H, H 2 /ν, ρ(ν/H) 2 , and ∆T * as characteristic quantities for velocity, length, time, pressure, and temperature, respectively and for a Newtonian incompressible fluid in the Boussinesq approximation, are:

∂u i ∂t + u.∇u i = - ∂p ∂x i + A l Gr θ δ i,2 + ∇ 2 u i for i = 1, 2 and 3, (3)
∂θ ∂t + u.∇θ = 1 P r ∇ 2 θ, (4)
∇.u = 0, (5)
where u = (u 1 , u 2 , u 3 ), t, and p are the dimensionless velocity, time, and pressure, respectively, δ i,2 is the Kronecker symbol, while θ is the dimensionless temperature:

θ = [T * -1/2(T * h + T * c )
]/∆T * . The thermal boundary conditions for θ are: θ = 0.5 at the hot wall, θ = -0.5 at the cold wall and ∂θ/∂η = 0 at all other walls.

3 The 2D 1/2 models for natural convection

Limitations of the 2D model and formulation of the 2D 1/2 models

The inability of the 2D model to predict the velocity and temperature fields within the cavity (already reported by Botton et al. [START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF]) is clearly illustrated in Fig. 2 in the case P r = 0.01 and Gr = 10 7 . In this figure, the velocity fields are given as coloured magnitude levels and vectors, while the temperature fields are depicted by black isotherms. The same colorbar scale is used for Fig. 2(a-d). Figure 2(a) reports the result obtained with the 3D model and averaged in the transverse direction x 3 over the small width W , while Fig. 2(b) shows the result obtained with the 2D model. Note first that the 3D solution is steady (as it will be in all the cases studied in this section 3), while, for these parameters, the 2D solution is unsteady. This indicates that the 2D flow is so intense that it even becomes unsteady in this case. To highlight these large velocities, a representative snapshot has been chosen during the oscillatory evolution for the 2D plot. The contrast between the results shown in Fig. 2(a) and Fig. 2(b) is striking. The 2D model gives a highly convective flow with two recirculation zones of significant size, one in the bottom left corner near the hot wall and the other in the top right corner near the cold wall; a more elongated convective circulation with a lower velocity magnitude is obtained with the 3D model. These differences in the velocity fields also lead to differences in the temperature fields; we can indeed see that the 2D isotherms are highly deformed compared to the 3D isotherms, especially in the regions where the velocity is maximum in the 2D results. These differences confirm that the 2D model, although interesting in terms of simulation time, is not reliable for the study of natural convection within such a confined cavity. The side walls indeed affect the obtained flow structure via the no-slip condition they induce, and this effect modifies both the velocity and temperature fields.

The so called 2D 1/2 model is proposed [START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF] in order to take into consideration the dissipative effects caused by the side walls without losing the main advantages of 2D modeling in terms of simulation time and cost and then keeping the possibility to make large range parametric studies. The 2D 1/2 models are obtained by averaging the equations of momentum, heat and mass conservation (Eqs. [START_REF] Nagira | Impact of melt convection induced by ultrasonic wave on dendrite growth in Sn-Bi alloys[END_REF][START_REF] Shevchenko | The effect of natural and forced melt convection on dendritic solidification in Ga-In alloys[END_REF][START_REF] Wang | A quasi two-dimensional benchmark experiment for the solidification of a tin-lead binary alloy[END_REF]) along the P models are steady solutions while the 2D solution is unsteady. A representative snapshot has then be chosen for the 2D plot.

x 3 direction. For that, we assume both a constant temperature in this direction, which is realistic for low values of the Prandtl number corresponding to liquid metals and metallic alloys in such a confined cavity (this assumption will be checked later in subsection 3.3, and a well defined transverse profile for the velocity expressed as a function f (x 3 ). This approach reduces the problem from a 3D formulation to a 2D modified model referred to as 2D 1/2 model, where all the variables and dependencies in x 3 have disappeared. This 2D 1/2 formulation is also known as shallow-water formulation, especially in the field of free surface hydraulics.

The dimensionless governing equations in the 2D 1/2 modeling are:

∂ ūi ∂t + ū.∇ ūi = - ∂ p ∂x i + A l Gr θ δ i,2 + ∇ 2 ūi + S i for i = 1, and 2, (6)
∂ θ ∂t + ū.∇ θ = 1 P r ∇ 2 θ, (7)
∇.ū = 0,

where ū, p and θ are the mean values of velocity, pressure, and temperature along the x 3 direction, and are the functions of x 1 and x 2 only. The term S i in the momentum conservation equation (Eq. ( 6)) is the only difference between the classical 2D model and the 2D 1/2 model and is expressed as:

S i =   1 - 1 A w Aw 0 f 2 (x 3 )dx 3   ū.∇ ūi + 1 A w ∂f (x 3 ) ∂x 3 Aw 0 ūi . (9)
Two different 2D 1/2 models have been tested, based on different functions f (x 3 ) for the velocity profile along x 3 . The first model (referred to as the 2D

1/2 P model) is based on the parabolic Poiseuille profile (see Fig. 4) given by:

f P (x 3 ) = 6 A w 2 (A w -x 3 )x 3 . (10)
The corresponding expression for S i noted S P,i is:

S P,i = - 1 5 ū.∇ ūi - 12 A w 2 ūi . (11)
The second model (referred to as the 2D

1/2
H model) is based on a Hartmann-type profile presenting two boundary-layer-like regions along the side walls and a uniform bulk velocity:

f H (x 3 ) = F   1 - cosh( x 3 -Aw 2 δ ) cosh( Aw 2δ )   where F = 1 1 -2δ Aw tanh( Aw 2δ ) . ( 12
)
δ is a free parameter which corresponds to the thickness of the boundary layer normalized by H and can be adjusted (see Fig. 4 and discussion in the next section). A second expression of S i noted S H,i is then obtained:

S H,i = 1 -F 2 3 2 - 1 2 tanh 2 ( A w 2δ ) - 3δ A w tanh( A w 2δ ) ū.∇ ūi -F 2 δA w tanh( A w 2δ ) ūi . ( 13
)
The 3D, 2D, 2D

1/2
H and 2D

1/2 P models are implemented in COMSOL Multiphysics 5.0 software. The equations are discretized using first order elements for velocity, pressure and temperature, while both streamline and crosswind diffusions are used as stabilization methods. The Delaunay-based mesh generator is used to generate the different meshes: Delaunay triangulation for the 2D models and Delaunay tetrahedralization for the 3D model. The time stepping method used in the implicit time-dependent solver is a free backward differentiation formula (free BDF). For the 2D models, a fully coupled approach is used: at each time step, the non-linear problem is solved by a Newton method with an updated Jacobian matrix, each linear Newton correction being obtained by the parallel sparse direct solver (PARDISO). In contrast, a segregated solver is used for the 3D model: a Newton method is still applied at each time step, but for the different variables separately and with minimal Jacobian update, and the Newton corrections are obtained by an iterative solver (GMRES) for the velocities and pressure and by a direct solver (PARDISO) for the temperature. Several mesh refinements were used for the different models, and the minimum refinement ensuring the mesh independence of the obtained solution was used.

The results obtained with the 2D

1/2 H model and the 2D 1/2 P model for the case of Fig. 2 (P r = 0.01 and Gr = 10 7 ) are shown in Fig. 2(c) and Fig. 2(d), respectively. An improvement in the obtained velocity and temperature fields, which are now closer to the 3D fields, is seen for both the 2D 1/2 models. The highest velocity magnitudes observed in the purely 2D model is reduced and the recirculation zones have disappeared. The additional energy dissipation induced by the no-slip condition at the side walls has considerably improved the results. A further comparison of the results shows that the 2D

1/2
H model gives a more accurate prediction of the 3D velocity and temperature field than the 2D 1/2 P model. Indeed, a closer look at the isotherms shows that there are still differences between the 2D 1/2 P and 3D temperature fields, while an almost perfect prediction of the 3D temperature field is obtained with the 2D

1/2 H model. As reported by Botton et al. [START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF], the 2D 1/2 model featuring a Hartmann-type profile to describe the shape of the transverse velocity field is best suited to predict convection for low Prandtl number fluids in such a confined cavity, compared to the Poiseuille profile based 2D

1/2 P model. It is important to note, however, that the Hartmann-type profile used for this case was adjusted by choosing an optimal value for the parameter δ, as explained in the following section.

Adjustment of the parameter δ in the 2D

1/2

H model

For low Prandtl number fluids at any value of the Grashof number, the parameter δ can be optimized in order to get the best comparison between the 2D

1/2

H results and the 3D results. The idea is to obtain a map in the form: δ optimal = g(P r, Gr). For given Pr and Gr, we first investigate the influence of the parameter δ on the precision of the 2D

1/2 H model compared to the 3D model. The comparison is performed through the maximum relative error between the 2D and 3D temperature fields. This maximum relative error is estimated on the values of the dimensionless temperature θ at a set of comparison points. In practice, this set of comparison points consists in a 50 × 30 regular grid, defined in the plane of the 2D simulations on one side, and in the vertical middle plane of the 3D simulations on the other side. The grid used was checked to be delicate enough for an accurate comparison. This maximum relative error is then given by:

∆ max (2D) = M ax{|θ 3D (i) -θ 2D (i)|} for i = 1, 1500 (14)
where i is the index of the grid points. Such error can be estimated as well for the purely 2D model, the 2D

1/2 P model and the 2D

1/2 H model. The variation of ∆ max with δ for the 2D 1/2
H (δ) model is shown in Fig. 3 for P r = 0.01 and Gr = 10 6 . It is compared with the error obtained with the 2D 1/2 P and purely 2D models. Figure 3 shows that there is only one value of δ minimizing the maximum relative error in temperature ∆ max . This value of δ is denoted δ optimal . For this case, P r = 0.01 and Gr = 10 6 , δ optimal = 0.235×A w . The 2D 1/2 H (δ optimal ) model gives a temperature field with a maximal relative error of 0.206% compared to the 3D temperature field, while the 2D 1/2 P and the purely 2D models give a maximal relative error in temperature of 2.6% and 14.2%, respectively.

The change in the value of the parameter δ modifies the shape of the assumed transverse velocity profile in the 2D 1/2 H formulation. Different shapes of this profile (as illustrated in Fig. 4) imply a different contribution of the wall effect to the intensity of the convective flow within the cavity. As in the case presented in Fig. 3, for any Gr and P r values, there is only one value of δ, δ optimal , optimizing the wall Botton et al. [START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF] only reported two values of δ optimal , δ optimal = 0.1×A w for the case Gr = 1.296×10 7 and P r = 0.0047, and δ optimal = 0.15 × A w for the case Gr = 4.34 × 10 6 and P r = 0.0129. They also proposed the use of the latter value of δ optimal for other (Gr, P r) numbers as a good compromise to approach 3D results. In the present work, a large range of Grashof and Prandtl numbers is investigated. The Prandtl number is varied from P r = 10 -3 to P r = 10, which includes a wide variety of fluids, from liquid metals to gases and aqueous liquids, while the Grashof number is taken in the range 10 3 ≤ Gr ≤ 5 × 10 7 . In order to get the map δ optimal = g(P r, Gr), the value of δ optimal is first determined as described in Fig. 3 for each pair (P r, Gr).

Figure 4 gives the shape of the profiles assumed to model the different transverse velocity profiles along x 3 . The parabolic Poiseuille profile given by Eq. ( 10), as well as different Hartmann-type profiles expressed by Eq. ( 12) and obtained for several δ values, are plotted. We see that by varying the value of δ, we can go from really flat profiles (as for δ/A w = 0.052) to profiles close to the Poiseuille profile (as for δ/A w = 0.5). A closer look to the profiles in Fig. 4(b) reveals that slight differences with the Poiseuille profile still remain for δ = 0.5 and that these differences can be further decreased for larger values of δ.

The variations of the chosen profiles affect the source term denoted S i in Eq. ( 6) and expressed by Eq. ( 13) for the 2D

1/2
H models and Eq. ( 11) for the 2D 1/2 P model. The expression of S i can be written in the more general form as:

S i = a ū.∇ ūi + b ūi , (15)
where a and b are constants in the case of the 2D

1/2 P model and are the functions of δ in the 2D

1/2 H model. Table 1 gives a comparison between the values of a and b in Eq. [START_REF] Carozzani | Direct simulation of a solidification benchmark experiment[END_REF] for high δ values and those for the Poiseuille type model. The slight differences shown in Fig. 4(b) are reflected on the values of a and b which directly influence the momentum conservation equation (Eq. ( 6)). This justifies the fact that δ optimal is found to be greater than A w /2 in some cases, while the maximum relative error ∆ max thus obtained is still improved compared to that found with the 2D P r = 0.4 and Gr = 10 5 , we obtain δ optimal = 0.57 × A w with an error ∆ max of 0.72%, while the 2D 1/2 P model error is 1.39%. Other values of δ optimal were found to reach up to 0.69 × A w and ∆ max is still improved compared to the 2D These results indicate that δ can be interpreted as a boundary layer thickness (as defined by Botton et al. [START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF]) only for the small values of δ featuring flat profiles. It is otherwise nothing more than an adjustment parameter allowing to fit Hartmann-type profiles to the transverse variation of the velocity in our 3D results.

Domain of validity of the 2D 1/2 models

The main interest of the 2D 1/2 H,P models is to predict the 3D temperature field assuming that the temperature is uniform along the transverse x 3 direction. For low values of the Prandtl number, moderate Grashof numbers and within a cavity confined in this x 3 direction, the temperature is indeed expected to be independent of x 3 . However, for larger Prandtl and Grashof numbers, it is not so clear. We then explore the (Gr, P r) parameter space in order to define the parameter range in which the temperature variation along x 3 remains small enough for our 2D 1/2 approach to be valid.

A first qualitative insight is shown in Fig. 5. The transverse temperature contours obtained with the 3D model in the horizontal section plane at x 2 = 0.2 (see Fig. 1) are given for a Prandtl number equal to 0.01 and Grashof numbers in the range from 10 3 to 5 × 10 7 . For each case, we also give C, the maximum curvature of the dimensionless temperature contours along the x 3 direction. C is defined as the maximum difference of the dimensionless temperature between the vertical middle plane of the 3D cavity and its side wall and can be expressed as: As shown in Fig. 5, the value of C increases with the increase of Gr. Low values of C (typically 1% or less) indicate that the isotherms along the x 3 direction are almost straight lines. It is the case for Gr = 10 3 , Gr = 10 5 , Gr = 10 6 and even Gr = 5.5 × 10 6 , where we have C equal to 0.00033%, 0.033%, 0.26%, and 1.2%, respectively. Larger values of the Grashof number such as Gr = 1.5 × 10 7 , Gr = 2.5 × 10 7 , and Gr = 5 × 10 7 give a maximum curvature C of 2.9%, 4.1%, and 6.1%, respectively, indicating larger curvatures of the isotherms and then a less uniform temperature field along x 3 .

C = M ax{|θ center (i) -θ wall (i)|} = M ax{ (T * center (i) -T * wall (i) ∆T * }. (16
The change of the heat transfer process within the cavity between the two differentially heated end walls with the variation of the Grashof number is also seen in Fig. 5. The heat transfer evolves from mainly diffusive for Gr = 10 3 (straight equidistant isotherms), to highly convective for Gr ≥ 10 7 (more curved isotherms tightened close to the end walls).

From the results shown in Fig. 5, C = 2.5% is arbitrarily chosen to define the limit below which the use of the 2D 1/2 H,P models is appropriate. We use the fact that the qualitative view of the case at Gr = 10 7 (which gives C = 2.1%) shows that the curvature of the isotherms still remains moderate, even though the heat transfer is already highly convective. On the other side, we have to note that C is a maximum value of the curvature obtained from all the comparison grid points in the vertical middle plane and those at the side walls. C is then a maximum localized value, much greater than the mean curvature C mean obtained by averaging the differences of temperature at the grid points. For example for Gr = 10 7 , we have C = 2.1% and C mean = 0.29% only.

Using the just defined tolerance level (C = 2.5%) and for given values of P r, we can calculate the critical Grashof number below which the use of the 2D 1/2 H,P models is appropriate. As an example, from Fig. 5, we see that for P r = 0.01 the critical value of the Grashof number is between Gr = 10 7 and Gr = 1.5 × 10 7 . Such calculations will be done for chosen values of P r within the range [0.001, 10], namely 0.001, 0.004, 0.01, 0.04, 0.1, 0.4, 1, 4, and 10. In each case, the Grashof number is increased by steps from Gr = 10 3 to values of Gr for which the curvature C is about the tolerance level C = 2.5%.

For each studied value of the Prandtl number, Figure 6(a) presents the variation of C (the maximum curvature of the 3D isotherms along x 3 ) with the increase of the Grashof number until the tolerance level is exceeded. Figure 6(a) shows that, as expected, the increase of the Grashof number increases the curvature C for all the Prandtl numbers. It also shows that the tolerance level is more quickly exceeded for higher Prandtl numbers. This is justified by the fact that in the energy conservation equation (Eq. ( 4)), the convective term responsible for the isotherms deformations will be more efficient for a lower diffusive term contribution, which is associated with a stronger P r value. parameter space for the validity of the 2D 1/2 models. These points are extracted from (a) and correspond, for each P r value, to the last point below the 2.5% tolerance level. Red line: approximated limit of validity in the (Gr, P r) plane.

The highest Grashof numbers reached with a maximum curvature C just below the tolerance level (defined as the critical Grashof numbers) are associated with their corresponding Prandtl numbers and presented in Fig. 6(b) as black dots, giving an estimation of the validity limit for the 2D 1/2 approach. This limit can be fairly well approximated by a straight line (red line in Fig. 6(b)) obtained by a linear fit to all the critical points. The condition of validity thus obtained is:

P r × Gr 0.8749 < 10 4.074 . (17)

Results and discussion

In subsection 3.1, a qualitative comparison between the results obtained with the different models (3D, 2D, 2D

1/2

H and 2D

1/2 P ) has been presented for the case where P r = 0.01 and Gr = 10 7 (Fig. 2). In this section, detailed values of the maximum relative error in temperature, ∆ max , obtained by the different models when compared to the 3D model, are presented, first for the studied Grashof numbers at P r = 0.01 (Fig. 7), and then for all the studied cases within the range of variation of (P r, Gr) (Fig. 8). Note that in all these cases, three cases, which belong to the extreme border of the parameter range, will give an unsteady solution with the 2D model, namely (P r = 0.01 and Gr = 10 7 , the case of Fig. 2), (P r = 0.004 and Gr = 2.5 × 10 7 ) and (P r = 0.001 and Gr = 5 × 10 7 ). In these cases, a representative snapshot is used for the comparisons as in Fig. 2. The variation of ∆ max with Gr for the different 2D and 2D 1/2 models is shown in Fig. 7 For all the considered Grashof numbers from 10 3 to 10 7 , the relative error obtained with the 2D 1/2 H model is lower than that obtained using the 2D

1/2 P model, itself lower than the 2D error. Although for Gr = 10 3 the 2D error is only 0.25% and is quite acceptable, it exceeds 10% for Gr = 10 5 to reach an error of 18.5% for Gr = 10 7 . The error is reduced with the 2D H models, respectively. Each case (Gr, P r) is associated with the value ∆ max given by the considered model, which gives a set of points defined as (Gr, P r, ∆ max ), or as ∆ max function of Gr and P r, a form which can be easily handled to obtain the given contours. It should be noted that a logarithmic scale is used for the contours of ∆ max , and the numbers on the isolines (represented as black thin lines) are the corresponding powers of 10 (for example, the line labeled '-3' corresponds to an isocontour of ∆ max = 10 -3 ). Figure 8(a) gives the contours of ∆ max obtained by the 2D model. ∆ max is higher than 10% for a large part of the studied (Gr, P r) cases and exceeds 25% near the validity limit. Once more, this illustrates the lack of reliability of the 2D model to study natural convection in this kind of confined cavity.

As shown in Fig. 8(b), the 2D 1/2 P model allows a better prediction of the 3D temperature field than the pure 2D modeling. ∆ max , however, still exceeds 7% when Gr > 5.5 × 10 6 for both P r = 0.004 and P r = 0.01, and reaches a maximum value of 10.5% for Gr = 2.5 × 10 7 and P r = 0.004. The 2D 1/2 H model gives a still better prediction of the 3D results than the 2D 1/2 P model and, as shown in Fig. 8(c), this is valid for all the considered (Gr, P r) cases. The highest value, ∆ max = 2.2%, is reached for this model for Gr = 2 × 10 6 and P r = 0.04. The comparison between Fig. 8(b) and Fig. 8(c) demonstrates that an adequate choice of the transverse velocity profile allows to take into consideration the right amount of dissipative effects induced by the side walls and then to have a better prediction of the 3D temperature field. optimal /A w , fitted curve Figure 9: δ optimal obtained for all the studied (Gr, P r) cases as a function of the Grashof number (red points). The existence of more than one point for the same Grashof number is due to the different Prandtl numbers studied; the fact that these points are all near a unique master curve shows the small influence of the Prandtl number. The black line corresponding to Eq. ( 18) gives a good fit to all these points.

The additional cost when choosing the Hartmann type profile is that one more parameter has to be adjusted, namely δ. It would then be interesting to find a general correlation giving the optimal value of δ for any values of (Gr, P r). All the values of δ optimal obtained from all the cases considered in Fig. 8(c) are plotted in Fig. 9 as a function of the Grashof number (red points). We can notice that the dependence of δ optimal on the Prandtl number is so small that all the data appear to be fairly well fitted by a single curve δ optimal = f (Gr) (black line in Fig. 9) given by δ optimal varies from 0.69 × A w for the lowest studied value of the Grashof number (Gr = 10 3 ), which corresponds to a transverse velocity profile quite close to the Poiseuille profile, to 0.052 × A w for Gr = 5×10 7 and P r = 0.001 , which corresponds to a very flattened Hartmann-type transverse velocity profile (see the profiles in Fig. 4).

δ optimal /A w = a erf(b (Log(Gr) -3) + c) -exp(Log(Gr) -3 -d) + e, (18)
Equation 18 then allows a direct estimation of δ optimal as a function of the fluid and the working conditions and therefore a direct use of the 2D

1/2
H model. This is particularly useful for parametric studies where good predictions can then be obtained without having to make 3D simulations, highly demanding in terms of computational time and resources.

Test in solid-liquid phase change configurations

The efficiency of the 2D 1/2 models for the study of natural convection within a transversally confined cavity was highlighted in the previous section. The comparison of these models with the 3D model was mainly made in terms of temperature field rather than velocity field. This is justified by our interest in solidification processes: in this context it is the temperature field which indeed controls the solid-liquid phase change. In this section we investigate the possibility to use the 2D 1/2 models for the prediction of the solid-liquid front shape and position within the cavity during a solidification process. In a first step, we intend to determine the steady state position of the interface obtained when the cold (hot) wall is set at a constant temperature below (above) the melting temperature of the studied material. After that, we investigate the ability to simulate a full tin solidification process performed in the 'AFRODITE' experimental setup [START_REF] Wang | A quasi two-dimensional benchmark experiment for the solidification of a tin-lead binary alloy[END_REF][START_REF] Wang | An investigation of the influence of natural convection on tin solidification using a quasi two-dimensional experimental benchmark[END_REF][START_REF] Hachani | Experimental analysis of the solidification of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Boussaa | 2D and 3D numerical modeling of solidification benchmark of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Hachani | Multiscale statistical analysis of the tin-lead alloy solidification process[END_REF][START_REF] Boussaa | Macrosegregations in Sn-3 wt%Pb alloy solidification: Experimental and 3D numerical simulation investigations[END_REF]. In both cases, solidification of pure tin is concerned. Its melting temperature is T * m = 505 K (see its main physical properties in Table 2).

Steady state situation

In this first situation, where the time evolution of the solidification will lead to a steady state equilibrium between the two phases, the results obtained with the different 2D models (2D, 2D

H and 2D

1/2 P ) will be compared to those obtained with the 3D model, mainly in terms of the position and shape of the solidliquid front at equilibrium. Two cases will be considered: in the first case (referred to as 'case (1)'), the fixed hot wall temperature (above T m ) is T h = 512 K and the fixed cold wall temperature (below T m ) is T c = 502 K, whereas in the second case (referred to as 'case (2)'), these temperatures are T h = 507.5 K and T c = 503.5 K, respectively. In both cases, the diffusive solution will be used as an initial field for the temperature, giving thus a perfectly vertical initial solid-liquid front. P numerical models for the phase change study are derived from the previous models for natural convection and are again implemented in COMSOL Multiphysics 5.2 software. The Comsol 'Deformed Geometry' physics is now used: it allows mesh deformation and solid-liquid front tracking. This corresponds to the multi-domain approach formerly proposed be Wolf et al. [START_REF] Wolf | Solidification of a pure metal at a vertical wall in the presence of liquid superheat[END_REF] and used by Avnaim et al. [START_REF] Avnaim | Comparison of threedimensional multidomain and single-domain models for the horizontal solidification problem[END_REF][START_REF] Avnaim | Numerical and experimental study of the traveling magnetic field effect on the horizontal solidification in a rectangular cavity. Part 2: Acting forces ratio and solidification parameters[END_REF]. The defined domain, corresponding to the chosen geometry, is divided into two regions (one for each phase). The advantage of this method is that the physical properties of the liquid and solid phases are defined separately in each region. It also allows the use of two different mesh refinements. This is very useful as, in the solid phase region, only the heat conduction equation is solved:

Symbols

∂T * s ∂t * = k s ρ Cp s ∇ * 2 T * s , (19)
and no high mesh refinement is required, whereas in the liquid region, the energy equation is coupled with the momentum and mass conservation equations, which implies a higher mesh refinement. Following the results of Assael et al. [START_REF] Assael | Reference data for the density and viscosity of liquid copper and liquid tin[END_REF] suggesting a linear variation of the density for liquid tin on a large range of temperature near T m , the Boussinesq approximation can still be applied in the liquid domain to describe the density variation in the buoyant force. The density is otherwise considered as a constant.

The equations in the liquid phase for the 2D 1/2 models are:

∂ ū * i ∂t * + ū * .∇ * ū * i = - 1 ρ ∂ p * ∂x * i + g β (T * l -T * m ) δ i,2 + ν ∇ * 2 ū * i + S * i for i = 1 and 2, ( 20
)
∂ T * l ∂t * + ū * .∇ * T * l = k l ρ Cp l ∇ * 2 T * l , (21)
∇ * . ū * = 0. (22)
The term S * i in Eq. ( 20) refers to the source term of the 2D 1/2 models. By removing it we obtain the 2D conservation equations and by increasing the variation of the index i to 3 and removing the bars over the variables, we obtain the 3D equations.

The solid and liquid regions share a moving boundary which is the phase change front. The energy balance at this front is expressed by the Stefan condition [START_REF] Wolf | Solidification of a pure metal at a vertical wall in the presence of liquid superheat[END_REF]:

ρ ∆H f v f .n = (φ l -φ s ).n, (23)
where ∆H f is the latent heat of fusion, v f is the solid-liquid front velocity and n is the normal vector at the front. φ l and φ s are the heat fluxes coming from the liquid and solid regions, respectively.

The dimensional term S * i can be deduced from S i already defined in Eqs. [START_REF] Carozzani | Direct simulation of a solidification benchmark experiment[END_REF][START_REF] Haagen | A comprehensive benchmark of fixedgrid methods for the modeling of melting[END_REF][START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF]. We can write:

S * i = ν 2 H 3 S i = a ū * .∇ * ū * i + ν H 2 b ū * i .
Using A w = W/H and δ = δ * /H, we can express a and b as a function of W and δ * and finally get for S * i :

S * H,i = 1 -F 2 3 2 - 1 2 tanh 2 ( W 2δ * ) - 3δ * W tanh( W 2δ * ) ū * .∇ * ū * i -F 2 δ * W tanh( W 2δ * ) ν ū * i (24) with F = [1 - 2δ * W tanh( W 2δ * )] -1 and S * P,i = - 1 5 ū * .∇ * ū * i - 12 W 2 ν ū * i . ( 25
)
As shown in the previous section, the choice of δ in the 2D

1/2

H model has to be optimized in order to get reliable results. The ratio δ * optimal /W is given by Eq. ( 18) and is only a function of the Grashof number Gr which needs then to be known.

The calculation of Gr (Eq. ( 2)) requires the knowledge of ∆T * /L l , i.e. the temperature difference between the hot and the cold boundaries of the liquid zone over the distance between these boundaries, L l . In our phase change situation, the hot boundary is at T * h and the cold boundary is the solid-liquid front at T * m . Only L l has then to be estimated. A simple way to obtain an estimation of L l in this first test case is to calculate the steady state position of the front in the case of a purely diffusive heat transfer. The heat balance at steady state (v f = 0) then gives:

φ l d = φ s ⇔ k l T * h -T * m L l = k s T * m -T * c L -L l ⇔ L l = L [k l (T * h -T * m )/(k s (T * m -T * c ))] 1 + k l (T * h -T * m )/(k s (T * m -T * c )) . ( 26
)
Due to the presence of an additional convective heat transfer within the liquid region, L l calculated from Eq. ( 26) underestimates the size of the zone occupied by molten tin. The estimation of L l can be improved by using a heat flux in the liquid region, φ lc , taking into consideration the contribution of the convective heat transfer. In Eq. ( 5.63) in the chapter 'Internal natural convection', Bejan [START_REF] Bejan | Convection Heat Transfer[END_REF] reports the Berkovsky-Polevikov correlation giving the Nusselt number in rectangular cavities with differentially heated side walls:

N u H = 0.18 ( P r 0.2 + P r Ra H ) 0.29 ( L H ) -0.13 [START_REF] Maxworthy | The nonlinear growth of a gravitationally unstable interface in a Hele-Shaw[END_REF] where N u H and Ra H are the Nusselt and Rayleigh numbers based on the height H, respectively. We will use this correlation at the final steady state of the solidification where phase change and movement of the front are no more involved. The main approximation will be that one of the endwalls is the solid-liquid front which is not a vertical flat wall. Valuable estimations of L l will nevertheless be obtained.

With our notations, N u H can be expressed as

N u H = φ lc H k l (T * h -T * m ) , (28)
and Ra H does not depend on L l . The heat balance is then written:

φ lc = φ s ⇔ k l T * h -T * m H N u H = k s T * m -T * c L -L l . ( 29
)
Using Eq. ( 27) with L = L l , we obtain an equation for L l of the form

L -0.13 l (L -L l ) + C 1 = 0, (30)
where C 1 is a constant. The solution of Eq. ( 30) then allows to have a convecto-diffusive estimation of L l . Gr can then be calculated, and δ * is estimated with Eq. ( 18). The results obtained with the purely diffusive and the convecto-diffusive estimations, for case (1) and case [START_REF] Reinhart | Influence of natural convection during upward directional solidification: A comparison between in situ X-ray radiography and direct simulation of the grain structure[END_REF], are summarized in Table 3. Figure 10 illustrates the steady states reached by the different models in case [START_REF] Reinhart | Direct simulation of a directional solidification experiment observed in situ and real-time using X-ray imaging[END_REF]. For each model, both the temperature field (isotherms in the solid and liquid regions) and the velocity field (colored levels of magnitude and vectors in the liquid region) are shown. The comparison between Fig. 10(a) and Fig. 10(b), once more, illustrates the poor ability of the 2D model to predict the 3D steady state. The velocity magnitude within the melt is too high in the 2D model, inducing an irrelevant flow structure with marked recirculation zones in the upper-right and lower-left corners. As a consequence, the isotherms are too much deformed by the flow, the convective heat transfer between the hot wall and the solid-liquid front is too high, and the position of the front is not well estimated with a significantly smaller solidified volume than in the 3D model. P model, by imposing a dissipative effect to take into account the friction at the side walls, reduces the intensity of the convective structure within the melt and improves the 2D predictions for the front position. Differences, however, still remain when compared to the 3D model (intensity of the flow a little too strong and deformation of the isotherms). H model where δ * is obtained from the convecto-diffusive approach estimation. The better account of the walls induced dissipation in this case results in a still slightly reduced velocity magnitude and a better prediction of the solidification front.

A better comparison of the different solid-liquid front shapes and positions for both case (1) and case (2) is given in Fig. 11. More quantitative comparisons are also given in Table 3. Figure 11 clearly shows the improved estimation of the 3D front when changing the 2D models from a purely 2D model, to the 2D

1/2 P model, the 2D 1/2
H model with diffusive estimation of δ * , and finally the 2D

1/2
H model with convectodiffusive estimation of δ * . For each model, we measure the difference between the front position given by the considered model and the front position given by the 3D model (at different heights) and normalize it by the 3D position chosen as a reference. In case (1), the 2D model gives a front with a mean error of almost 12% and a maximum error higher than 18%. The 2D 1/2 P model reduces these errors to 4.84% and 6.06%, respectively. As shown in Fig. 11, in this range of error, the difference with the 3D front position is still not negligible. The 2D

1/2
H model with a basic diffusive approach to estimate δ * reduces significantly these errors to only 0.79% for the mean error and to 1.01% for the maximum error. Finally, the 2D

1/2
H model with a convecto-diffusive approach to estimate δ * is the most accurate in the prediction of the 3D solid-liquid front in comparison to the other models: the matching between the two fronts is better (see Fig. 11) and the mean and maximum errors are only 0.11% and 0.22%, respectively. Similar comparisons have been obtained in case (2) (see Fig. 11 and Table 3), which confirms the reliability of the 2D Table 3 also shows the different estimations of the liquid length L l used for the 2D

1/2 H model. As expected L l is underestimated with the diffusive approximation. The values obtained by the two methods for Gr and δ * are different but reasonably close. It should be noted that the condition of validity of the 2D 1/2 model given by Eq. ( 17) is fulfilled for both the diffusive and the convecto-diffusive approaches and in both cases (1) and ( 2

). 2D 2D 1/2 P 2D 1/2 H (diffusive) 2D 1/2 H (convecto-diffusive) L l - - (1) 0.0538 [m] (2) 0.0455 [m] (1) 0.0751 [m] (2) 0.0584 [m] Gr - - (1
) 5.818 × 10 6

(2) 2.461 × 10 6

(1) 4.174 × 10 6

(2) 2.041 × 10 6

δ * -- (1) 0.1106 × W (2) 0.1474 × W

(1) 0.1219 × W (2) 0.1596 × W Mean relative error of the front position compared to 3D (1) 11.80% (2) 25.80% (1) 4.84% (2) 6.00% (1) 0.79% (2) 0.74% (1) 0.11% (2) 0.14% Maximum relative error of the front position compared to 3D (1) 18.61% (2) 32.39%

(1) 6.06% (2) 6.82%

(1) 1.01% (2) 1.06% (1) 0.22% (2) 0.34%

Solidification process

The 2D

1/2

H model has been proven, therefore, to be the most appropriate to predict solid-liquid front shape and position. In this section the 2D 1/2 H model is used to simulate a tin solidification process from pure liquid to a full solid ingot. The results of the simulation are directly confronted to yet unpublished experimental results obtained with the experimental benchmark AFRODITE [START_REF] Wang | A quasi two-dimensional benchmark experiment for the solidification of a tin-lead binary alloy[END_REF][START_REF] Wang | An investigation of the influence of natural convection on tin solidification using a quasi two-dimensional experimental benchmark[END_REF][START_REF] Hachani | Experimental analysis of the solidification of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Hachani | Multiscale statistical analysis of the tin-lead alloy solidification process[END_REF].

Figure 12: Sketch of the experimental benchmark with the crucible (light grey) in the middle, along with the location and the spacing between the fifty thermocouples (orange dots). The two lateral heat (dark grey) each contain nine additional thermocouples with a specific arrangement (blue dots).

The setup, developed in the SIMAP-EPM laboratory, was prominent in numerical models validation [START_REF] Boussaa | 2D and 3D numerical modeling of solidification benchmark of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Boussaa | Macrosegregations in Sn-3 wt%Pb alloy solidification: Experimental and 3D numerical simulation investigations[END_REF][START_REF] Botton | A 2D 1/2 model for low Prandtl number convection in an enclosure[END_REF][START_REF] Zheng | Use of a mixed columnarequiaxed solidification model to analyse the formation of as-cast structure and macrosegregation in a Sn-10 wt%Pb benchmark experiment[END_REF][START_REF] Carozzani | Direct simulation of a solidification benchmark experiment[END_REF][START_REF] Gandin | Direct modeling of structures and segregations up to industrial casting scales[END_REF][START_REF] Rad | Validation of a model for the columnar to equiaxed transition with melt convection[END_REF][START_REF] Zheng | Role of fragmentation in as-cast structure: numerical study and experimental validation[END_REF]. It consists of mainly parts :

(1) The parallelepiped crucible of L = 10 cm length, W = 1 cm width and H = 6 cm height, represented in light grey in Fig. 12 and made of 0.5 mm thick stainless steel sheet and containing the sample of pure tin.

(2) Two lateral heat exchangers represented in dark grey in Fig. 12 made of copper, each one with a heating element part and a controlled temperature water box. They both have the same cross section (1 cm × 6 cm) than the crucible to ensure a more homogeneous heat flux release.

(3) Sixty-eight K-type thermocouples with a recording frequency of 1 Hz for the temperature measurement. Fifty of them are laser-welded on the largest surface of the steel sheet in the crucible and arranged in a 5 rows × 10 columns regular lattice with 1 cm spacing between them in both the horizontal and the vertical direction. Their location is represented with orange dots labelled from L 1 to L 50 in Fig. 12. These thermocouples will give an approximate view of the temperature field evolution inside the sample during the solidification process. The eighteen other thermocouples are equally divided between the two heat exchangers, labelled from FR1 to FR9 and from FL1 to FL9 for the right-hand side and the left-hand side exchangers, respectively and are represented with blue dots in Fig. 12. They allow the measurement of the heat flux extracted from the sample. The temperature difference between the two heat exchangers (monitored between FL4 and FR4) is held at 20 K and the cooling rate is set to -0.03 K/s. (4) A Kirchhoff box insures thermal insulation by minimizing radiative heat loss via the lateral surfaces.

(5) The whole setup is put inside a controlled atmosphere to reduce heat transfer by air convection.

For more details see the publications of Hachani et al. [START_REF] Hachani | Experimental analysis of the solidification of Sn-3 wt%Pb alloy under natural convection[END_REF][START_REF] Hachani | Multiscale statistical analysis of the tin-lead alloy solidification process[END_REF]. The experimental process includes different steps: the initial pure tin sample is first enclosed in the crucible. It is then slowly heated, using both the right and left heat exchangers, its temperature (at about 550 K) exceeds the melting temperature T * m of tin: the sample is therefore fully melted. A temperature difference of 20 K is then between the two heat exchangers and held constant until the end of the process. The two heat exchangers are first kept at constant temperatures (both above T * m ) for more than half an hour to let natural convection establish. After that, the same cooling rate of -0.03 K/s is applied to both exchangers, so that the temperatures of the hot and cold walls decrease with time. Due to the applied temperature difference, the cold wall temperature is the first to go below T * m , prior to the hot wall temperature. Solid zones will then progressively grow from the end walls until the sample is fully solidified. A typical evolution of the temperature at a thermocouple is shown in Fig. 13 (green line corresponding to the thermocouple L 26 ). After the initial plateau corresponding to the settling of the convection, the applied cooling rate induces a regular decrease of the temperature down to the melting temperature T * m . A kind of plateau is then observed at T * m , before a further decrease of the temperature. This plateau at the melting temperature, more or less pronounced depending on the thermocouple position (see Fig. 16 and the corresponding comments), was used to precisely calibrate the thermocouples in the sample. An important point for the reliability of our numerical simulation is to determine appropriate thermal boundary conditions. This is not an easy task due to a varying thermal contact resistance between the exchangers and the melt. This is achieved using the heat flux conservation between the left/right exchanger and the hot/cold wall along x 1 direction, respectively. For the left-hand side we write:

φ lef t exchanger = φ hot wall ⇔ k copper dT * dx * 1 0 - = k tin dT * dx * 1 0 + (31)
The derivative of T * with respect to x * 1 is estimated at 0 -with a second order backward differentiation 31)) during the pure tin solidification process. The nominal cooling rate is -0.03 K/s, T * m is the melting temperature. The curves of T * h (t) and T * c (t) are used as boundary conditions in our numerical simulations.

using each series of three thermocouples at the same height in the left heat exchanger:

dT * dx * 1 0 - = 1 3 [ 3 (T F L1 + T F L3 + T F L5 ) -4 (T F L2 + T F L4 + T F L6) + (T F L7 + T F L8 + T F L9 ) 2 e * 1 ], (32)
where e * 1 = 1.5 cm is the horizontal distance between the thermocouples of the heat exchangers. In contrast, dT * /dx * 1 is estimated at 0 + with a first order forward difference between T * h , the temperature at the hot wall, and the temperature at the three first thermocouples in the sample, at the same heights than those on the exchanger:

dT * dx * 1 0 + = 1 3 [ 3 T * h -(T L 1 + T L 21 + T L 41 ) e * 2 ], (33)
where e * 2 = 0.5 cm is the distance between the hot wall and the first column of thermocouples. At any time, these two estimations of the derivative of T * on either side of the hot wall at x * 1 = 0 are used in Eq. ( 31) to deduce T * h , the temperature at the hot wall; the same method is used at the cold wall for T * c . The obtained estimations of the temperature at the hot and cold walls, starting from the time where the differential heating is imposed to induce natural convection, are plotted in Fig. 13 as a function of time. These temperatures are used as thermal boundary conditions in our numerical model. The numerical simulation is done with a 2D

1/2 H model quite similar to the one used for the steady state situation where the two phases co-exist during the whole simulation time. For the solidification process considered here, only the liquid phase is present within the cavity at the beginning, but, in order to use the deformed geometry and the dynamic mesh method, a solid region has to be defined and be present in the studied geometry domain all along the simulation time. To fulfil this constraint, two very thin layers of solid regions, on either side of the cavity where the solidification is expected to start, were added from the beginning. Their thickness of 0.1 mm represents only one thousandth of the domain length. ). The thickness of the layers has been reduced successively from 3 mm to 0.1 mm, i.e. until it no longer affects the solidification time at all. An automatic remeshing technique is used: remeshing occurs when the mesh has undergone a distortion larger than a given limit, allowing the increase of the number of elements in the solid regions as they grow in size, while the number of elements in the liquid zone is progressively reduced allowing to keep an optimal density of elements.

The value of the parameter δ * used in the 2D 1/2 H model is estimated from the experimental temperature data (Fig. 13) at the beginning of the solidification process when T * c reaches the melting temperature T * m . The effective temperature difference between the hot and cold walls is then equal to 3.6 K, which gives a Grashof number Gr = 2.06 × 10 6 and yields δ * = 0.1590. Figure 15 presents the results obtained with the 2D 1/2 H model (on the left) along with the corresponding experimental results (on the right). The solid phase is represented in brown and the liquid phase in blue. It should be noted that the experimental results are given within the zone covered by the lattice of 50 thermocouples. This zone is also represented with the black rectangle on the numerical results. The time evolution of the solidification process is illustrated in Fig. 15 showing a fair agreement between the numerical and experimental results. As T * c is the first temperature in the sample to go below T * m , the solidification starts from the right-hand side. A second front appears later during the process at the left-hand side of the cavity. Both fronts progress, reducing the liquid region, and eventually meet at a position x * 1 closer to the left side of the cavity. In the experiment x * H model) solid-liquid front positions at different times during the solidification process. The black rectangle inside the cavity domain on the numerical plots represents the zone where the experimental data (given by the lattice of fifty thermocouples) is available.

t * = 1201 s, while the numerical simulations give x * 1 2.9 cm and t * = 1283 s. We can also note that, at t = 4750 and 5000 s, the shape of the fronts in the experiment is not as flat as in the simulations. This effect cannot be due to convection at these times where the two fronts have already appeared and the remaining liquid zone is almost isothermal, but could be a consequence of imperfect thermal boundary conditions. Moreover, the shape of the experimental front, determined from the array of thermocouples with only five elements along the height of the sample, is also obtained with a limited precision. This fair agreement between the numerical results, obtained with the 2D 1/2 H model, and the experimental results is confirmed in Fig. 16, which gives the time evolution of the temperatures monitored at different thermocouples (L 21 , L 22 , L 23 , L 27 , L 28 , and L 29 , all shown in Fig. 12) and the corresponding temperatures calculated with the 2D and the melting temperature T * m are also represented. The temperatures of all the thermocouples decrease slowly until reaching the melting temperature T * m . They stay constant for a while at T * m , before continuing to decrease. The different steps of the solidification process are also visible in Fig. 16. T * c (bold blue curve) is the first temperature signal to go below T * m , indicating the front apparition at the right-hand side. This front will move to the left inside the sample: as it will pass at a thermocouple position, the temperature of that thermocouple will leave the pseudo-plateau value at T * m to decrease linearly. For the chosen thermocouples in Fig. 16, as expected, it will first affect T * L 29 (thermocouple closer to the cold wall), and then T * L 28 and T * L 27 . The second front on the left-hand side then appears when T * h (bold red curve) decreases below T * m , and, similarly, the displacement of this front to the right will successively affect T * L 21 (thermocouple closer to the hot wall), T * L 22 and T * L 23 . This pseudo-plateau behaviour at T * m , almost not marked for T * L 29 , is particularly visible for the thermocouples L 21 to L 28 that are still in the liquid zone when the second front appears. Indeed, the numerical and experimental results have shown that, past this time, the whole liquid region between the two fronts is quasi-isothermal at a temperature close above T * m . The temperature at the thermocouples in this remaining liquid region will then stay quasi-constant, close to T * m , until one of the fronts reaches the thermocouple position, inducing then a linear decrease of the temperature characteristic of the solid zone. In the array of thermocouples, T * L 23 at x * 1 = 2.5 cm is the last temperature to go below T * m , giving a long pseudo-plateau at T * m and indicating that the two fronts meet at x * 1 close above this value. Note that the numerical temperature data are only available up to the full solidification of the sample at t * = 5170 s: this is due to the fact that the numerical model must keep a thin liquid zone and has then to stop just before the junction of the two fronts.

Conclusion

We have presented 2D 1/2 models for the simulation of natural convection and pure material solidification in a horizontal differentially heated slender cavity. The configuration corresponds to a reference solidification benchmark, namely the AFRODITE experimental setup. The 2D 1/2 models are obtained by averaging the equations of momentum, heat and mass conservation in the transverse x 3 direction. The integration is made assuming a constant temperature along x 3 , which seems realistic in such a confined cavity and for low values of the Prandtl number corresponding to liquid metals and metallic alloys. A well defined velocity profile has also to be assumed. Two velocity profiles have been considered, the usual parabolic Poiseuille profile (2D 1/2 P model) and the Hartmann-type profile featuring two boundary layers of adjustable size δ on the sides of a uniform bulk (2D 1/2 H model). The first assumption stating the uniformity of the temperature along the transverse x 3 direction has been checked to be valid in a large range of the control parameters P r (the Prandtl number) and Gr (the Grashof number). Qualitative but also quantitative comparisons between the 2D, 2D H and 3D calculation results obtained in the study of natural convection have been presented. It has been found that the pure 2D model poorly reproduces the physics of the flow, leading to maximum relative deviations in temperature with respect to the 3D results greater than 10% for a large part of the studied (Gr, P r) cases and even exceeding 25% in some cases. The 2D 1/2 P model improves the results, but deviations in temperature higher than 10% can be obtained in some cases. The 2D 1/2 H model, with the adequate choice of δ, is the more efficient 2D model. The highest temperature deviation in the studied parameter region does not exceed 2.2%, proving its ability to be a reliable substitute to the 3D calculations in the prediction of the temperature field for natural convection. An empirical relationship giving δ optimal as a function of the control parameters (Gr in practice) is derived for a direct use in practical cases. The use of 2D, 2D H and 3D models for pure material solid-liquid phase change problems is also presented. The solid-liquid front tracking method is based on a dynamic mesh method. A first typical situation where the hot and cold walls temperatures are above and below the melting temperature, respectively, evolves towards an equilibrium state between the solid and liquid phases. For this situation two different cases with different boundary conditions were considered. For the first case, it has been found that the equilibrium position and shape of the interface is hardly reproduced by the 2D model, with a smaller solid zone and a maximum relative error of 18.61% on the position. This is attributed to an overestimate of convective effects since this model does not account for viscous dissipation at the side walls. In contrast, the prediction of the 2D 1/2 P model is fair (error reduced to 6.06%) and that of the 2D 1/2 H model is very good, with a relative error smaller than 1%. These results were confirmed in the second case. Finally, a full experimental solidification process performed in the 'AFRODITE setup' and its numerical 2D

1/2 H counterpart are reported. A thorough treatment of the experimental temperature data has been made to define appropriate unsteady thermal boundary conditions for the numerical model and deduce the optimal value of δ. A fair agreement is found between the numerical and experimental results in terms of the time evolution of the solid-liquid fronts position, despite discrepancies in the fronts shape. The 2D

1/2 H model, with the right estimation of the dissipative effects induced by the side walls, is thus found to be very efficient for the study of natural convection and solid-liquid phase change problems. It has the advantage of keeping a 2D formulation, which makes it far less computer time consuming than 3D models. We see this work as a validation step on the way to very efficient modeling for binary alloys solidification in slender configurations. The obtained model could indeed be light enough to be embedded in meta-models for a statistical approach. This could be convenient to simulate the recent statistical investigation of solidification by Hachani et al. [START_REF] Hachani | Multiscale statistical analysis of the tin-lead alloy solidification process[END_REF] in the AFRODITE setup, but also to model numerous very thin setups used for X-ray in-situ observations of solidification processes.

Figure 1 :

 1 Figure 1: Geometrical configuration of the study: Differentially heated parallelepipedic cavity with left and right end walls held at constant temperatures T * h and T * c , respectively. The other walls are adiabatic. Characteristic planes are shown: the vertical middle plane (green plane) and a horizontal section plane (hatched plane).

Figure 2 : 2 H 2 P

 222 Figure 2: Velocity fields (coloured magnitude levels and vectors) and temperature fields (isotherms) for P r = 0.01 and Gr = 10 7 obtained with different models: (a) 3D model (averaged over the width W ), (b) 2D model (snapshot), (c) 2D 1/2 H model with a Hartmann-type transverse velocity profile and (d) 2D 1/2 P model with a Poiseuille transverse velocity profile. For these parameters, the solutions obtained with the 3D, 2D 1/2 H and 2D 1/2

Figure 3 : 2 H 2 P

 322 Figure 3: Maximum relative errors in temperature ∆ max of the different models compared to the 3D temperature field in the central middle plane for P r = 0.01 and Gr = 10 6 . Blue curve: variation of ∆ max with the parameter δ for the 2D 1/2 H model; orange line: ∆ max for the 2D 1/2 P model; black line: ∆ max for the 2D model.

1 / 2 P 1 Figure 4 : 2 P

 12142 Figure 4: The profiles used in the 2D 1/2 approach to model the velocity variation along the transverse x 3 direction: the parabolic Poiseuille profile used in the 2D 1/2 P model and the Hartmann-type profiles (depending on the value of δ) used in the 2D 1/2 H models. (a) General view of the profiles; (b) zoom for larger values of δ.

Figure 5 :

 5 Figure 5: Qualitative view of the transverse temperature contours obtained with the 3D model in the horizontal section plane at x 2 = 0.2 (see Fig. 1) for P r = 0.01 and different Grashof numbers. For eachGrashof number the corresponding maximum curvature C of the isotherms along the x 3 direction is given.

Figure 6 :

 6 Figure 6: (a) Variation of C (the maximum curvature of the 3D isotherms along the x 3 direction) with the Grashof number Gr for different Prandtl numbers from P r = 0.001 to P r = 10. The black bold line represents the chosen level of tolerance for the curvature C. (b) Limit points in the (P r, Gr)parameter space for the validity of the 2D 1/2 models. These points are extracted from (a) and correspond, for each P r value, to the last point below the 2.5% tolerance level. Red line: approximated limit of validity in the (Gr, P r) plane.

2 HFigure 7 : 2 H

 272 Figure 7: Maximum relative error in temperature ∆ max obtained with the 2D (black curve), 2D 1/2 P (red curve) and 2D 1/2 H (green curve) models as a function of the Grashof number for P r = 0.01.

1 / 2 P 2 H

 122 model. It is only 9.5 × 10 -6 for Gr = 10 3 , does not exceed 1% until Gr > 4 × 10 5 , but still reaches 7.9% for Gr = 10 7 . The 2D 1/model still reduces this error by almost one order of magnitude: the error is 1.2 × 10 -6 for Gr = 10 3 and does not exceed 1.05% for the highest considered Grashof number at P r = 0.01, i.e. Gr = 10 7 .Instead of plotting other graphs (∆ max vs Gr), for the other values of the Prandtl number, all the values of ∆ max obtained for the different (Gr, P r) cases studied are used to plot iso-contours of ∆ max in the (Gr, P r) parameter plane. The results obtained for each type of modelization are shown in Fig.8.

Figure 8 : 2 P

 82 Figure 8: Contours of ∆ max given by the studied models in the (Gr, P r) plane, restricted to the region of validity of the 2D 1/2 approach. The black dots in (a), (b) and (c) represent the studied (Gr, P r) cases, the bold line represents the limit of validity of the 2D 1/2 model, and the black thin lines correspond to ∆ max values given as powers of 10. (a) Contours of ∆ max given by the 2D model, the error exceeds 20% near the limit of validity; (b) contours of ∆ max given by the 2D 1/2 P model: a maximum error of 10.5% is reached for Gr = 2.5 × 10 7 and P r = 0.004; (c) contours of ∆ max given by the 2D 1/2 H model: a maximum error of 2.2% only is reached over all the region of validity for Gr = 2 × 10 6 and P r = 0.04; (d) colorbar giving the values of ∆ max .

Figure 8 (

 8 Figure 8 (a), (b) and (c) give the contours of the maximum relative error in temperature ∆ max for the 2D, 2D 1/2 P and 2D 1/2

  with a = -0.2662, b = 1.256, c = -3.158, d = 7.018, and e = 0.04084.

Figure 10 : 2 H

 102 Figure 10: Velocity fields (colored magnitude level and vectors) in the liquid phase region and temperature fields (black isotherms (1 K step)) obtained with the different models for a steady state phase change situation. The solid-liquid front position and shape correspond to the isotherm T * = T * m . (a) 3D model (averaged over the third dimension x 3 ), (b) 2D model, (c) 2D 1/2 H model and (d): 2D 1/2 P

Figure 10

 10 Figure 10(d) shows that the 2D 1/2

Figure 10 (

 10 Figure 10(c) corresponds to the 2D 1/2

Figure 11 :

 11 Figure 11: Solid-liquid front shapes and positions obtained by the different models for the two steady state phase change situations, case (1) and case (2).

Figure 13 :

 13 Figure 13: Typical evolution of the temperature at a thermocouple (L 26 ) and estimated experimental temperatures T * h (t) and T * c (t) at the hot and cold walls (see Eq. (31)) during the pure tin solidification process. The nominal cooling rate is -0.03 K/s, T * m is the melting temperature. The curves of T * h (t) and T * c (t) are used as boundary conditions in our numerical simulations.

Figure 14 ( 2 HFigure 14 : 2 H

 142142 Figure 14(a) shows the initial mesh used in the 2D 1/2 H model; one can barely see the two layers at the hot

1 2. 5 Figure 15 :

 1515 Figure 15: Experimental (right) and numerical (left and obtained with the 2D 1/2

1 / 2 HFigure 16 :

 1216 Figure 16: Numerical and experimental time evolution of the temperature at mid-height in the cavity at some thermocouples positions during the solidification process, along with the hot and cold walls temperatures (T * h and T * c , respectively). The melting temperature T * m is also represented.

Table 1 :

 1 Variation of the variables a and b in S

		1/2 P model.
	a constants a a a	a a a 2D 1/2
	a	-0.1893	-0.1951	-0.1972	-0.2
	b	-460.01	-444.64	-439.15	-432
	comparison with the constant values in the 2D	1/2 P model.	1/2 H model in

δ/A w = 0.5 δ/A w = 0.75 δ/A w = 1 Poiseuille i (Eq. (

15

)) with the increase of δ in the 2D

Table 2 :

 2 Physical properties used to model tin in the vicinity of its melting temperature T * m . The density ρ is considered to have the same value in the two phases in order to ensure the overall mass conservation.

		Value	Description
	ν	2.5 × 10 -7 m 2 /s	kinematic viscosity
	ρ	7070 kg/m 3	density
	β	2.2 × 10 -5 1/K	coefficient of thermal expansion
	k s	60 W/(m×K)	solid phase thermal conductivity
	k l	30 W/(m×K)	liquid phase thermal conductivity
	Cp s	228 J/(kg×K)	solid phase heat capacity at constant pressure
	Cp l	242 J/(kg×K)	liquid phase heat capacity at constant pressure
	∆H f	59.6 kJ/kg	latent heat of fusion
	T * m	505 K	melting temperature
	The 3D, 2D, 2D 1/2 H and 2D 1/2

Table 3 :

 3 Comparison of the estimated liquid length L l , Grashof number Gr and parameter δ * between the 2D H (convective) approaches. Mean and maximum error in solid-liquid front position for the different 2D and 2D 1/2 models in comparison to the 3D front position. (1) and (2)refer to case (1) and case (2), respectively.

	1/2 H (diffusive) and the 2D 1/2

Acknowledgement

This collaborative work was supported by the PHC Maghreb Partnership Program No. 36951NG. A grant of the PROFAS B+ program given to I. H. is also gratefully acknowledged.