Victor D Fachinotti

Steven Le Corre

Nicolas Triolet

Manuel Bobadilla

Michel Bellet

Two-phase thermo-mechanical and macrosegregation modelling of binary alloys solidification with emphasis on the secondary cooling stage of steel slab continuous casting processes

Keywords: thermo-mechanical, macrosegregation, two-phase, continuous casting

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

Introduction

Macrosegregation, i.e. the lack of homogeneity of the solute concentration at the whole scale of a solidified product, is a central problem since it strongly influences the further workability of the cast products and their mechanical properties. Macrosegregation is the result of slow interdendritic flow of molten liquid and transport of alloying elements at the product scale. In most numerical models, only the natural convection induced by thermal and solutal gradients is taken into account. The influence of the solid kinetics on the fluid flow is rarely modelled. However, as summarized by Flemings [START_REF] Flemings | Our understanding of macrosegregation: past and present[END_REF], the macrosegregation of chemical species may depend strongly on the deformation of the solid skeleton in the mushy zone. This is especially the case in continuous casting: as bulging occurs between the supporting rolls, there exists a large mushy zone which deforms together with the solid shell.

In the literature, the deformation of a semi-solid mushy metallic alloy has been most frequently studied at low solid fraction, typically lower than 0.6, in view of modelling processes such as the injection of metallic slurries, also known as "thixoforming" processes. In this context, the semi-solid alloy is considered as an equivalent one-phase continuous medium. At low solid fraction, less than 0.2, the material is generally considered as Newtonian. Following Einstein [START_REF] Einstein | Eine neue bestimmung der moleküldimensionen[END_REF], many authors have proposed models in which the viscosity depends on the solid fraction and on microstructural parameters such as the mean diameter of the solid particles or of their specific surface. At higher solid fractions, typically between 0.2 and 0.6, the previous models have been complemented by a shear rate dependency of the viscosity, which is generated by the interactions of solid particles and by their viscous deformation. The shear-thinning property of metallic alloys in the mushy state, that is the decrease of the apparent viscosity with shear rate, or more generally with the equivalent strain rate, is hence introduced. Alloys in such semi-solid conditions are characterized by values of strain rate sensitivities between 0.2 (typical value in the solid state at high temperature) and 1 (limit value in the Newtonian case). The rheological models derive from the suspension theory, like the model of Kumar et al. [START_REF] Kumar | Constitutive modelling and characterisation of the flow behavior of semi-solid metal alloy slurries[END_REF], in which the apparent viscosity results from two contributions: the hydrodynamic one and a second contribution due to the interactions between the particles in suspension. This model is then complemented by an equation for the kinetics of agglomeration-disagglomeration between the particles. Besides the fact that such models are not suited for high solid fractions, their drawback is the large number of parameters which makes their practical identification extremely difficult. Zavaliangos and Lawley [START_REF] Zavaliangos | Numerical simulation of thixoforming[END_REF] have applied this model to the thixoforming process.

A second category of mushy zone models is composed of two-phase models, but assuming a fixed rigid solid phase. Such models have been developed either by using the mixture theory or by applying the spatial averaging method to the conservation equations (mass, momentum, energy, solute) on a representative elementary volume (REV) of the mushy material. They have been proposed to study solidification in columnar or equiaxed regime, accounting in this latter case for the transport and settling of solid grains. The interaction between solid and liquid phases is expressed by the Darcy law. Such models have been initially developed and implemented in finite volume numerical models by Bennon and Incropera [START_REF] Bennon | A continuum model for momentum, heat and species transport in binary solid-liquid phase change systems. I. Model formulation[END_REF] and Beckermann et al (Ni and Beckermann [6], Wang and Beckermann [START_REF] Wang | Equiaxed dendritic solidification with convection: Part 1. multiscale/multiphase modeling[END_REF]). Finite element implementations have been proposed by Ahmad [START_REF] Ahmad | Numerical simulation of transport processes in multicomponent systems related to solidification problems[END_REF] in two dimensions and by Desbiolles et al [START_REF] Desbiolles | Micro-/macrosegregation modeling in casting: a fully coupled 3D model[END_REF] and Bellet et al [START_REF] Bellet | A 3D-FEM model solving thermomechanics and macrosegregation in binary alloys solidification[END_REF] in three dimensions.

A third family of models consists of effective two-phase models, in which the deformation of the solid phase is accounted for. These models derive from the soils mechanics theory and apply to the mushy zone when it is "coherent" and saturated with liquid. The "coherency" of a mushy zone is reached when the solid phase is continuous enough to sustain stresses significantly higher than those in the liquid phase. The corresponding solid fraction, g s cohe , is typically in the range 0.6-0.7 (the lower values corresponding to a dendritic structure and the higher ones to a globular structure). At solid fractions higher than g s cohe , the solid phase constitutes a continuous solid skeleton which is deformable and in strong interaction with the interstitial fluid. The deformation of this skeleton determines the flow of the liquid phase, which is governed by the Darcy law. The pressure of the interstitial liquid determines in turn the deformation of the solid skeleton. Schematically, the mushy zone then is supposed to behave like a sponge saturated with liquid. This kind of approach has been found quite relevant for slow deformation of semi-solid alloys, especially when liquid segregation can be observed during the deformation. Like the second family, the models proposed in the literature are based on a spatial averaging approach. However, most of them are limited to isothermal conditions and thus cannot be applied to the study of solidification process, as they do not account for mass transfer from the liquid phase to the solid phase that is caused by solidification or fusion. We can quote in this category the following references: Lalli [START_REF] Lalli | A model for deformation and segregation of solid-liquid mixtures[END_REF], Nguyen et al [START_REF] Nguyen | Theoretical and experimental study of the isothermal mechanical behaviour of alloys in the semi-solid state[END_REF], and Martin et al [START_REF] Martin | Viscoplastic behaviour of porous metallic materials saturated with liquid. part i: Constitutive equations[END_REF]. The implementation in two-dimensional finite elements numerical models has been developed by Toyoshima and Takahashi [START_REF] Toyoshima | A numerical simulation of forming processes for semi-solid materials[END_REF], Bay et al [START_REF] Bay | Numerical modelling of extrusion for a semi-solid alloy. application to rheological parameters identification in Sn-Pb alloys[END_REF]. In this third family, Mo and coworkers (M'Hamdi et al [START_REF] Hamdi | Two-phase modeling directed toward hot tearing formation in aluminum direct chill casting[END_REF], Nicolli et al [START_REF] Nicolli | Modeling of macrosegregation caused by volumetric deformation in a coherent mushy zone[END_REF]) have proposed, in view of studying the continuous casting of aluminium alloys, a model in which mass transfer is accounted for. In their two-dimensional finite element implementation, the conservation of momentum of both the liquid and the solid phase is solved, including Darcy interaction terms, together with the resolution of the conservation of mass, energy and solute for the whole mixture. However, the coupling with solute transport (macrosegregation) is not accounted for: a unique relation between solid fraction and temperature must be prescribed a priori, which is only acceptable for small segregations. In addition, this approach is limited to very low differential velocities between the two phases, the effects of differential advection being neglected in the equation for energy conservation.

In the present paper, we focus on the thermo-mechanical modelling of alloys by the end of solidification, that is to say for cases where the solid fraction g s exceeds g s cohe . The mushy material is then considered as an effective two-phase medium. On one hand, the solid material is considered as an incompressible viscoplastic material, obeying a constitutive equation of power-law type. Invoking homogenization results [START_REF] Geindreau | Investigation of the viscoplastic behaviour of alloys in the semi-solid state by homogenization[END_REF], its macroscopic flow rule is viscoplastic, including compressibility, so that the solid continuum can be seen as a deformable compressible porous medium. On the other hand, the liquid phase is intrinsically Newtonian. At the macroscopic scale, its momentum interaction with the solid skeleton is described by the Darcy law of flow through a porous medium. Then, the present contribution enters the third family of models previously mentioned. In addition, the present formulation includes mass transfer between liquid and solid, and the couplings with heat transfer and solute transport is taken into account. The two-dimensional numerical implementation of the corresponding mass, momentum, energy and solute conservation equations has been carried out in the finite element software R2SOL, developed at CEMEF. The paper is structured as follows. In Section 2, the coupled conservation equations are presented. In Section 3, their space and time discretisation and the finite element resolution are detailed. In Section 4, the coupling algorithms are discussed. Finally, in Section 5 a validation test of the two-phase mechanical solver is presented, followed by an application to the continuous casting process.

General macroscopic balance equations

At the microscopic scale, inside each phase, the thermo-mechanical evolution is assumed to be governed by the usual mass, momentum, energy and solute balances. In this work, the balance equations for the mixture, at the (macroscopic) scale of a representative elementary volume (REV), are obtained using the spatial averaging method over a fixed control volume V 0 . This method is classical and will not be detailed here (see [START_REF] Hassanizadeh | General conservation equations for multi-phase systems: 1. averaging procedure[END_REF][START_REF] Ni | A volume-averaged two-phase model for transport phenomena during solidification[END_REF][START_REF] Wang | Equiaxed dendritic solidification with convection: Part 1. multiscale/multiphase modeling[END_REF][START_REF] Rappaz | Numerical modelling in materials science and engineering[END_REF] for further details on its basic principles). The notations used in this paper can be summarized as follows. Denoting χ k the characteristic function of phase k (= 1 in phase k and 0 elsewhere), for any function ψ we define: the intrinsic average value:

ψ k = ψ k = 1 V k V 0 ψ(x)χ k (x) dV, (1)
the classic average value:

ψ k = 1 V 0 V 0 ψ(x)χ k (x) dV, (2)
where V k denotes the volume occupied by phase k in the REV. By introducing the volume fraction of phase k, defined as

g k = V k V 0 , ( 3
)
both average values (1) and ( 2) become related by

ψ k = g k ψ k . ( 4
)
Besides, since the solidifying alloy in the mushy state is considered as a saturated twophase medium, the volume fractions of solid and liquid phases always satisfy the following relationship:

g s + g l = 1. (5)
Now, applying the spatial averaging process to microscopic balance equations in each phase k (k = s, l), and using the previous notations, one obtains the following set of macroscopic equations:

Momentum: ∇ • (g k σ k ) + M k + g k ρ k g = ∂ ∂t (g k ρ k v k ) + ∇ • (g k ρ k v k × v k ), (6)
Mass:

∂ ∂t (g k ρ k ) + ∇ • (g k ρ k v k ) = Γ k , (7)
Energy:

∂ ∂t (g k ρ k h k ) + ∇ • (g k ρ k h k v k ) + ∇ • q k = Q k , ( 8
) Solute: ∂ ∂t (g k w k ) + ∇ • (g k w k v k ) + ∇ • j k = J k , ( 9
)
where ρ denotes the density, v the velocity field, σ the Cauchy stress tensor, g the gravity vector, h the enthalpy per unit mass, q the heat flow vector, w the solute concentration per unit volume, and j its flux. The terms Γ , M , Q and J are associated with the exchanges of mass, momentum, energy and solute, respectively, between the two phases.

Let us notice that the present model being applied to coherent mushy zones only, the velocities of the solid phase remain very slow (typically 0.01 m/s in continuous casting for instance). Therefore the inertial terms can be neglected in the momentum equation for the solid phase.

The following subsections will now detail the additional assumptions and constitutive models adopted in this paper for each variable appearing in those macroscopic balance equations.

Mass conservation

The local mass balance at the interface between phases ensures that Γ s + Γ l = 0 [START_REF] Rappaz | Numerical modelling in materials science and engineering[END_REF]. Then, summing equations [START_REF] Wang | Equiaxed dendritic solidification with convection: Part 1. multiscale/multiphase modeling[END_REF] for the liquid and solid phases, and assuming that the densities of the two phases remain constant in the solidification interval (but not necessarily equal), we get:

ρ s ∇ • (g s v s ) + ρ l ∇ • (g l v l ) = ∂g s ∂t (ρ l -ρ s ), (10)
where the saturation assumption (5) was introduced. Then, dividing by ρ l , equation ( 10) can be rewritten as

(1 -Δε tr )∇ • (g s v s ) + ∇ • (g l v l ) = ∂g s ∂t Δε tr , ( 11
)
where

Δε tr = ρ l -ρ s ρ l ( 12
)
denotes the relative change of volume associated with solidification (often negative for metallic alloys).

It should ne noted that the assumption of constant liquid and solid densities affects the solidification interval only. In this interval, ρ s and ρ l may be taken equal to the density of the alloy at nominal solidus and liquidus temperature, respectively. Outside this interval, only one phase is present and the mass equation results from equation [START_REF] Wang | Equiaxed dendritic solidification with convection: Part 1. multiscale/multiphase modeling[END_REF], yielding the usual one-phase mass conservation equation: ∂ρ ∂t + ∇ • (ρv) = 0.

Momentum conservation

The spatial averaging method used in this work is efficient to obtain in a simple way the macroscopic governing equations of the semi-solid alloy but does not enable to go further in the specifications of the macroscopic model. Reliable constitutive equations would require more sophisticated approaches such as homogenization [START_REF] Auriault | Etude du comportement macroscopique d'un milieu poreux saturé déformable[END_REF][START_REF] Geindreau | Investigation of the viscoplastic behaviour of alloys in the semi-solid state by homogenization[END_REF] associated with numerical simulation at the microscopic scale, but this is not within the scope of this work. The full definition of the two-phase model will simply be based on further constitutive assumptions, consistent with previous theoretical works.

Macroscopic constitutive equation for the liquid phase

At the microscopic scale, we assume that the liquid metal behaves as an incompressible Newtonian fluid. The incompressibility hypothesis is valid as long as the temperature range of the solidification interval remains narrow enough. In such a case, we can write:

s = 2μ l ε(v), (13)
where s = σ + pI is the deviatoric part of the stress tensor σ, p =tr (σ)/3 being the hydrostatic pressure, μ l the viscosity of the liquid, and ε(v) the strain rate tensor, defined as

ε(v) = 1 2 ∇v + (∇v) T . ( 14
)
As suggested by Ganesan and Poirier [START_REF] Ganesan | Conservation of mass and momentum for the flow of interdendritic liquid during solidification[END_REF] and Rappaz et al [START_REF] Rappaz | Numerical modelling in materials science and engineering[END_REF], we adopt the following model for the macroscopic deviatoric stress tensor:

Σ l = s l = 2μ l g l dev( ε(v l )). ( 15
)
where dev( * ) denotes the deviatoric part of tensor ( * ).

Macroscopic constitutive equation for the solid phase

Experimental studies on the behavior of metallic alloys at high temperature show that the response of the solid phase is well described by constitutive equations of the Norton-Hoff type [START_REF] Abouaf | Finite element simulation of hot isostatic pressing of metal powders[END_REF][START_REF] Geindreau | Constitutive behaviour of metal powder during hot forming. Part I: Experimental investigation with lead powder as a simulation material[END_REF][START_REF] Geindreau | Constitutive behaviour of metal powder during hot forming. Part II: Unified viscoplastic modeling[END_REF]. Like the liquid, the solid is assumed to be incompressible at the microscopic scale, such that its response can be characterized at this scale by the following constitutive equation:

s = 2K s √ 3 εeq m-1 ε(v), (16)
where K s and m denote the consistency and the strain rate sensitivity, respectively, and εeq the von Mises equivalent strain rate. The above equation can be rewritten in scalar form in terms of the equivalent stress (in the von-Mises sense) and the equivalent strain rate, as follows:

σ eq = 3K s √ 3 εeq m-1 εeq . ( 17
)
For solid fractions above the coherency fraction g s cohe , following the theoretical analysis of Geindreau and Auriault [START_REF] Geindreau | Investigation of the viscoplastic behaviour of alloys in the semi-solid state by homogenization[END_REF], the effective stress tensor

Σ s = s s -p s I + g s p l I ( 18
)
is expressed as a degree m homogeneous function with respect to the strain rate tensor ε s = ε(v s ). Note that such a result is valid as long as the viscosity of the liquid phase remains small compared to the solid phase apparent viscosity σ eq / εeq , that is to say for rather small values of the strain rate. Fortunately, this is generally the case for the solidification problems under consideration in this work. This property shows that the solid phase can be modelled as a compressible power-law fluid. We therefore adopt a compressible viscoplastic model [START_REF] Abouaf | Finite element simulation of hot isostatic pressing of metal powders[END_REF][START_REF] Nguyen | Theoretical and experimental study of the isothermal mechanical behaviour of alloys in the semi-solid state[END_REF][START_REF] Martin | Viscoplastic behaviour of porous metallic materials saturated with liquid. part i: Constitutive equations[END_REF]. Then, the constitutive equation coincides with equation ( 17), but with the equivalent stress and strain rate defined as

Σ s eq = A dev(Σ s ) : dev(Σ s ) + B ( tr (Σ s )) 2 1/2 , ( 19
)
ε s eq = 1 A ε s : ε s + 1 9B - 1 3A tr ( ε s ) 2 1/2 . ( 20
)
Those equations require two rheological functions A and B that depend on the solid volume fraction and for which several models can be found [START_REF] Abouaf | Finite element simulation of hot isostatic pressing of metal powders[END_REF][START_REF] Nguyen | Theoretical and experimental study of the isothermal mechanical behaviour of alloys in the semi-solid state[END_REF][START_REF] Martin | Viscoplastic behaviour of porous metallic materials saturated with liquid. part ii: Experimental identification on a Sn-Pb model alloy[END_REF][START_REF] Martin | Fracture behaviour in tension of viscoplatic porous metallic materials saturated with liquid[END_REF]. The constitutive equations of the solid phase at the macroscopic scale then take the form

Σ s = 3K s √ 3 ε s eq m-1 1 A ε s + 1 9B - 1 3A tr ( ε s )I , ( 21
)
or decomposed into deviatoric and spherical parts:

dev(Σ s ) = 3K s √ 3 ε s eq m-1 1 A dev( ε s ), ( 22
) - 1 3 tr (Σ s ) = -K s √ 3 ε s eq m-1 1 3B tr ( ε s ). ( 23
)
It is worth noticing that the preceding compressible viscoplastic constitutive model tends to the incompressible viscoplastic Norton-Hoff model as solidification completes. The rheological parameters A(g s ) and B(g s ) as defined in [START_REF] Abouaf | Finite element simulation of hot isostatic pressing of metal powders[END_REF][START_REF] Nguyen | Theoretical and experimental study of the isothermal mechanical behaviour of alloys in the semi-solid state[END_REF][START_REF] Martin | Viscoplastic behaviour of porous metallic materials saturated with liquid. part i: Constitutive equations[END_REF][START_REF] Martin | Viscoplastic behaviour of porous metallic materials saturated with liquid. part ii: Experimental identification on a Sn-Pb model alloy[END_REF][START_REF] Martin | Fracture behaviour in tension of viscoplatic porous metallic materials saturated with liquid[END_REF] satisfy A(1) = 3/2 and B(1) = 0, imposing in that way the incompressibility constraint by means of a penalty factor, the multiplier of tr ( ε s ) in equation [START_REF] Abouaf | Finite element simulation of hot isostatic pressing of metal powders[END_REF], which becomes infinite. Numerically, this factor is made arbitrarily high but finite, by replacing B(1) by B(g s incomp ), where g s incomp is a threshold for the solid fraction, above which incompressible behavior is assumed. In R2SOL, we adopt g s incomp = 0.999.

This simple model permits a first introduction to the compressibility of the solid skeleton in the mechanical behavior of a coherent mushy zone. The interested reader can refer to the recent work of Ludwig et al [START_REF] Ludwig | Rheological behavior of partially solidified al-cu alloys : experimental and numerical study[END_REF], who suggest complementing the previous model by introducing a state variable representing the cohesion of the skeleton. An evolution law of this variable is also proposed in order to represent the resistance of the inter-granular liquid films and solid bridges to tension. These modifications permit a better response of the model for small deformations, such as those encountered in continuous casting.

Exchange of momentum

According to the work of Ni and Beckermann [START_REF] Ni | A volume-averaged two-phase model for transport phenomena during solidification[END_REF], M k can be partitioned as:

M k = M d k + M p k , (24)
the first part being the contribution of deviatoric stresses and the second one the contribution of the isotropic part, and it can easily be shown that

M d s + M d l = 0 and M p s + M p l = 0. ( 25
)
The liquid being a Newtonian incompressible fluid with a very low viscosity, we will assume that the pressure equilibrium in the liquid phase is almost instantaneous. Subsequently, the interfacial pressures [START_REF] Ni | A volume-averaged two-phase model for transport phenomena during solidification[END_REF] in both phases (p * k ) equal the intrinsic average value of liquid pressure, i.e. its microscopic value:

p * l = p * s = p l . ( 26
)
Therefore, M p k can be expressed as follows [START_REF] Ni | A volume-averaged two-phase model for transport phenomena during solidification[END_REF]:

M p s = -M p l = -p l ∇g l = p l ∇g s . ( 27
)
On the other hand, depending on the solid fraction, the dissipative terms M d k are generally interpreted either as the drag force exerted by the liquid on the isolated solid grains [START_REF] Ni | A volume-averaged two-phase model for transport phenomena during solidification[END_REF] or the filtration force exerted by the liquid flowing through the solid, considered as a rigid porous medium [START_REF] Rappaz | Numerical modelling in materials science and engineering[END_REF][START_REF] Geindreau | Investigation of the viscoplastic behaviour of alloys in the semi-solid state by homogenization[END_REF]. In both cases, M d k can be modelled by a law of the generic type:

M d s = -M d l = g l C(g l )(v l -v s ). ( 28
)
Assuming the macroscopic isotropy of the medium, the interaction tensor C can be reduced to a scalar coefficient C, which is usually defined as

C = g l μ l κ , ( 29
)
being κ the permeability of the solid matrix in the mushy zone.

The mushy zone is assimilated to an isotropic porous medium whose permeability is defined by the Carman-Kozeny formula

κ = λ 2 2 g 3 l 180(1 -g l ) 2 , ( 30
)
where λ 2 is the secondary dendrite arms spacing.

Mechanical boundary conditions

The boundary conditions of the mechanical problem are:

v s = v l = v imp on ∂Ω u , ( 31
)
σ s n = σ l n = T imp , on ∂Ω s , ( 32
)
where σ s = Σ sg s p l I, v imp and T imp are the imposed surface velocity and traction on the non-overlapped portions ∂Ω u and ∂Ω s of the boundary ∂Ω of the analysis domain Ω; n is the normal unit vector pointing outwards ∂Ω.

Energy conservation

At the microscopic level, the Fourier law is used. Being λ the thermal conductivity, we have:

q = -λ∇T. ( 33
)
Our approach is based on the local thermal equilibrium assumption, which is valid because the liquid and the solid do exhibit neither too different thermal properties nor strong thermal interfacial barriers [START_REF] Auriault | Macroscopic modelling of heat transfer in composites with interfacial thermal barrier[END_REF]. Therefore, at the macroscopic scale, their average temperature is the same:

T s = T l = T. ( 34
)
The enthalpies of both phases can then be written as a function of the temperature T as follows:

h s = T T 0 c p (τ ) dτ and h l = h s + L, ( 35
)
where c p and L are the heat capacity and latent heat per unit mass, respectively. The average mixture enthalpy per unit mass takes the form

h = f s h s + f l h l = T T 0 c p (τ ) dτ + f l L, ( 36
)
where f l and f s are respectively the liquid and solid mass fraction, satisfying f l + f s = 1 in a saturated solid-liquid medium. Volume and mass phase fractions are related by

f k = ρ k g k ρ , ( 37
) being ρ = g s ρ s + g l ρ l ( 38
)
the average density of the solid-liquid mixture. By invoking the thermal equilibrium hypothesis [START_REF] Perchat | Mixed formulation and iterative solver for the parallelisation of a forging simulation software by mesh partitioning[END_REF], the energy conservation is written as a single equation which is the sum of the energy balance on both phases, where the terms of exchange Q s and Q l cancel themselves:

∂ ∂t ρh + ∇ • ρhv + ∇ • q = 0, (39)
where

ρh = (ρh) s + (ρh) l = g s ρ s h s + g l ρ l h l , ( 40
) ρhv = (ρhv) s + (ρhv) l = g s ρ s h s v s + g l ρ l h l v l , ( 41
)
q = q s + q l = -(g s λ s + g l λ l )∇T = -λ ∇T. ( 42
)
In equation ( 42), the scalar form of the macroscopic conductivity was obtained one again using the isotropy assumption of the medium at the macroscopic level.

By operating, we can rewrite equations ( 40) and ( 41) in the form

ρh = ρ h , ( 43
)
ρhv = ρv h + ρ (v l -v s )f s f l L, ( 44
)
where

ρv = g s ρ s v s + g l ρ l v l . ( 45
)
After elementary calculations, we obtain the following advection-diffusion equation for the average mixture enthalpy:

ρ ∂ h ∂t + ρv • ∇ h -∇ • (κ∇ h ) + S T = 0, ( 46
) being κ = λ ∂T ∂ h , ( 47
)
S T = L∇ • ( ρ (v l -v s )f s f l ) . ( 48
)
As shown by the above equation, the source term S T is originated by the relative movement of both phases in the mushy zone.

Thermal initial and boundary conditions

The heat balance equation ( 46) is subject to the initial condition:

T = T 0 at t = 0 in Ω, ( 49
)
and the following boundary conditions :

T = T w at t > 0 on ∂Ω T , (50)
(-λ ∇T ) • n = q w at t > 0 on ∂Ω q , ( 51
) (-λ ∇T ) • n = h(T -T ext ) at t > 0 on ∂Ω c , (52)
prescribing the temperature T w on ∂Ω T , the heat flux q w through ∂Ω q , and the heat exchange through ∂Ω c due to convection to the environment at temperature T ext with h as the convection coefficient; ∂Ω T , ∂Ω q , and ∂Ω c are non-overlapping portions of the boundary ∂Ω of Ω, being n the unit normal vector pointing outwards to ∂Ω.

Solute conservation

The solute flux vector j is determined at the microscopic scale by the first law of Fick, which can be written for the isotropic case as follows:

j = -D∇w, ( 53
)
being D the solute diffusion coefficient. By summing the solute conservation equations averaged for each phase, and considering that J s = -J l , we obtain:

∂ w ∂t + ∇ • wv + ∇ • j = 0, (54)
where

w = w s + w l = g s w s + g l w l , ( 55
)
wv = (wv) s + (wv) l = g s w s v s + g l w l v l , ( 56
)
j = j s + j l = -(D∇w) s -(D∇w) l = -g s D s ∇w s -g l D l ∇w l . ( 57
)
The diffusion coefficients D k are very small (particularly D s ), so that they are usually neglected. However, in this work we neglect the diffusion in the solid phase, but maintain the diffusion in the liquid phase since it has a stabilizing effect on the numerical solution. Then:

j ≈ -g l D l ∇w l . ( 58
)

Solute exchange at the solid-liquid interface: the microsegregation model

The solute exchange between the solid and the liquid phase is governed by interfacial phenomena occurring at the scale of the dendrite arms. This is a non-equilibrium process involving a very complex physics. However, several simplifying hypotheses are allowed for macrosegregation studies.

Besides the uniformity of temperature within the REV (already assumed for the energy balance), let us assume that both phases are at thermodynamic equilibrium at the interface, according to the alloy phase diagram. Then, at the interface, it satisfies:

w * s = k(T )w * l and w * l = ψ(T ), (59)
where the thermodependent partition coefficient k as well the equation of the liquidus line ψ(T ) can be determined from the phase diagram. Usually, both k and the slope of the liquidus line (m l ) vary only slightly with temperature, letting us assume them constant, as shown in Figure 1. Further, we assume the perfect diffusion of solute in the liquid phase at the microscopic scale, so that the solute concentration in the liquid inside the REV is uniform and particularly w * l = w l . We obtain then:

w * s = kw l ( 60
)
and

w l = T -T m m l , ( 61
)
where T m is the melting temperature of the pure substance.

Regarding solute diffusion in the solid phase, there are two classic limiting cases: • perfect diffusion in solid (lever rule), with the concentration in solid inside the REV uniform as well, and particularly w * s = w s , so we can write: w = (g l + kg s )w l ;

(62)

• zero diffusion in solid (Scheil law), being the average concentration given by

w = g l w l + g s 0 w * s dg s = g l w l + g s 0 kw l dg s , (63)
or, in differential form, by

∂ w ∂t = ∂ ∂t (g l w l ) + kw l ∂g s ∂t . ( 64
)
Besides these limiting cases, to account for a more realistic diffusion rate in the solid phase, one could solve the unsteady solute diffusion equation with moving interface in a onedimensional domain representing the characteristic length of the microstructure and/or that of the diffusion mechanism we consider. We refer to the works of Poirier et al [START_REF] Poirier | The energy and solute conservation equations for dendritic solidification[END_REF], Mo [START_REF] Mo | An internal variable description of solidification suitable for macrosegregation modeling[END_REF], Combeau et al [START_REF] Combeau | Modeling of microsegregation in macrosegregation computations[END_REF] and more recently the one-dimensional finite-difference model of Thuinet et al [START_REF] Thuinet | Computer simulation of microsegregation in the case of columnar growth involving a peritectic transformation for multicomponent steels[END_REF], who solve the diffusion equation in the solid and the liquid phases for alloys with multiple constituents and including the δ/γ transformation in the case of steels.

Final form of the solute conservation equation

We can reduce the microsegregation model to a general microsegregation function relating the average concentration in the solidliquid mixture to the liquid fraction g l and the average concentration in the liquid phase w l :

w = ϕ(g l , w l ), (65)
which is given by equations ( 62) and (63) for the particular cases of lever rule and Scheil law, respectively. Now, we can write the solute balance equation in the form of an advection-diffusion equation for the average concentration, as follows:

∂ w ∂t + v w • ∇ w -∇ • D∇ w + S w = 0, (66)
being

v w = v s + g l (v l -v s ) ∂w l ∂ w , ( 67
) D = g l D l ∂w l ∂ w , ( 68
)
S w = w ∇ • v s + w l ∇ • (g l (v l -v s )) . ( 69
)
In the same way as S T for the energy balance, the source term S w is originated by the relative movement and compressibility of both phases in the mushy zone. 66) is subject to the initial condition

Initial and boundary conditions for the solute problem Equation (

w = w 0 at t = 0 in Ω, ( 70
)
being w 0 the initial solute distribution, usually equal to the nominal solute content of the alloy. Further, we assume that there is no solute flux through the domain boundary, such that the following natural boundary condition applies:

∇ w • n = 0 at t > 0 on ∂Ω. ( 71
)
Final macroscopic balance equation set

Solid momentum: ∇ • Σ s -g s ∇p l + g l C(v l -v s ) + g s ρ s g = 0 Liquid momentum: ∇ • Σ l -g l ∇p l -g l C(v l -v s ) + g l ρ l g = ρ l ∂ ∂t (g l v l ) + ρ l ∇ • (g l v l × v l ), Mixture mass: (1 -Δε tr )∇ • (g s v s ) + ∇ • (g l v l ) = ∂g s ∂t Δε tr , Mixture energy: ρ ∂ h ∂t + ρv • ∇ h -∇ • (κ∇ h ) + S T = 0, Mixture solute: ∂ w ∂t + v w • ∇ w -∇ • D∇ w + S w = 0.

Numerical implementation

Mechanical model

This formulation was implemented in the two-dimensional code R2SOL, a finite element software using linear triangles and based on the P 1 + /P 1 mixed formulation [START_REF] Perchat | Mixed formulation and iterative solver for the parallelisation of a forging simulation software by mesh partitioning[END_REF][START_REF] Heinrich | Modélisation thermomécanique de la coulée continue d'acier en deux dimensions[END_REF][START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF]. These elements are shown in Figure 2.

As it is implicit in all the preceding formulation, the present model adopts the intrinsic average velocities of the solid and liquid phases as primal variables. This is remarked since in the isothermal version by Le Corre et al [START_REF] Le Corre | Two-phase modeling of metals solidification: a numerical approach for the thermo-mechanical problem[END_REF], which constitutes the starting point for the present work, a different choice for the velocity in the liquid phase was made, adopting as primary unknown not the intrinsic average v l but the classic average v l = g l v l . This previous choice was motivated by the convenience of using the model already developed in CEMEF [START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF] for melt flow in a rigid porous medium, which has v l . However, for the solid phase, the use of the intrinsic value v s was compulsory in order to avoid the division by g l in the interaction term [START_REF] Ludwig | Rheological behavior of partially solidified al-cu alloys : experimental and numerical study[END_REF], source of numerical problems in the vicinity of the solidus interface.

The advantages of the new model are evident when dealing with the Dirichlet boundary condition [START_REF] Mo | An internal variable description of solidification suitable for macrosegregation modeling[END_REF], since these conditions are posed directly over the primal variables. It is not the case with the previous model [START_REF] Le Corre | Two-phase modeling of metals solidification: a numerical approach for the thermo-mechanical problem[END_REF], where this condition becomes dependent on the liquid fraction g l .

Another advantage, in the authors' opinion, lies on the consistency of using the same interpolation functions for both intrinsic velocities v s and v l , contrary to the work of Le Corre et al [START_REF] Le Corre | Two-phase modeling of metals solidification: a numerical approach for the thermo-mechanical problem[END_REF], where variables of different nature were identically interpolated.

Then, let V be the space of "kinematically admissible" velocity fields and V 0 the space of "zero kinematically admissible" velocity fields. The virtual power principle states that the solution (v s , v l , p l ) ∈ V × L 2 (Ω) of the problem must satisfy:

Ω Σ s : ∇v * s dV - Ω p l ∇ • (g s v * s ) dV - Ω g l C(v l -v s ) • v * s dV - Ω g s ρ s g • v * s dV = ∂Ω s σ s n • v * s dS, ( 72
) Ω Σ l : ∇v * l dV - Ω p l ∇ • (g l v * l ) dV + Ω g l C(v l -v s ) • v * l dV - Ω g l ρ l g • v * l dV + Ω ρ l γ(g l , v l ) • v * l dV = ∂Ω s σ l n • v * l dS, ( 73
) Ω (1 -Δε tr )∇ • (g s v s ) + ∇ • (g l v l ) - ∂g s ∂t Δε tr p * l dV = 0, ( 74
) for all (v * s , v * l , p * l ) ∈ V 0 × L 2 (Ω), with γ(g l , v l ) = ∂ ∂t (g l v l ) + ∇ • (g l v l × v l ) = g l d l v l dt + v l ∂g l ∂t + ∇ • (g l v l ) , ( 75
)
where we introduce the "k-phase"-material derivative (k = s, l) of any variable ψ defined in the phase k as

d k ψ dt = ∂ψ ∂t + v k • ∇ψ. ( 76
)
3.1.1. Finite element formulation of the mechanical problem The above mixed solid velocityliquid velocity-pressure formulation requires interpolation functions that satisfy the Babuška-Brezzi stability conditions [START_REF] Arnold | Mixed finite element methods for elliptic problems[END_REF]. The present finite element formulation will be called (P 1 + ) 2 /P 1, since it constitutes the extension to the case of two velocity fields of the P 1 + /P 1 formulation previously used for the "one-velocity" models [34, 35, 36] † .

On each finite element, the velocity fields v s and v l are approximated as follows:

v s ≈ 3 i=1 N i (x)V s i + N b (x)V b s = v P 1 s + b s , ( 77
)
v l ≈ 3 i=1 N i (x)V l i + N b (x)V b l = v P 1 l + b l , ( 78
)
where V s i and V l i (i = 1, 2, 3) are the unknown velocities of the solid and liquid phase at the vertex node i, V b s and V b l are the additional solid and liquid velocity unknowns associated to the bubble node located at the center of each element, N i is the standard linear (P 1) shape function associated to the vertex node i, and N b is the bubble interpolation function, chosen here to be linear on each sub-triangle defined by the vertices and the center of the element. According to the definition of a bubble function, it is constructed such that N b = 0 on the edges of the element and N b = 1 at the bubble node.

The interstitial pressure p l is approximated by a standard linear interpolation:

p l ≈ 3 i=1 N i (x)P i , ( 79
)
where P i is the unknown pressure associated to node i.

Test functions v * s , v * s and p * l are approximated in the same way (this is the so-called Galerkin formulation [START_REF] Zienkiewicz | The Finite Element Method[END_REF], where shape and test functions coincide).

By replacing the trial and test functions by their finite element approximations given by equations (77) to (79) into equations (72) to (74), we obtain the spatially discretized form:

• Momentum in the solid phase: † "One-velocity" stands for models having only one velocity variable, such as the "one-phase" [START_REF] Perchat | Mixed formulation and iterative solver for the parallelisation of a forging simulation software by mesh partitioning[END_REF][START_REF] Heinrich | Modélisation thermomécanique de la coulée continue d'acier en deux dimensions[END_REF][START_REF] Bellet | A two-dimensional finite element thermomechanical approach to a global stress-strain analysis of steel continuous casting[END_REF] and the "pseudo-two-phase" [START_REF] Ahmad | Numerical simulation of macrosegregation: a comparison between finite volume method and finite element method predictions and a confrontation with experiments[END_REF][START_REF] Fachinotti | Three-dimensional finite element modelling of thermomechanics and macrosegregation in binary alloys solidification[END_REF][START_REF] Fachinotti | Modélisation tridimensionnelle par éléments finis de la macroségrégation au cours de la solidification d'un alliage binaire[END_REF][START_REF] Bellet | A 3D-FEM model solving thermomechanics and macrosegregation in binary alloys solidification[END_REF][START_REF] Gouttebroze | 3D-FEM modeling of macrosegregation in solidification of binary alloys[END_REF][START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF] models. For the "one-phase" model, there are only two phases, the solid-like and the liquid-like (with the interface defined by the coherence temperature), and the velocity variable coincides with the velocity of one or the other phase. On the other hand, although the "pseudo-twophase" considers the mushy zone as a mixture of solid and liquid phases, the solid phase is assumed to be rigid; in such a way, the liquid velocity remains the only velocity variable.

-P 1 contribution:

R v s =R v s v s + R v s v l + R v s b s + R v s b l + R v s p -F v s (80) R v s v s = Ω B T Σ s (v P 1 s ) dV + Ω g l CN T N dV V s (80a) R v s v l = - Ω g l CN T N dV V l (80b) R v s b s = Ω B T Σ s (b s ) dV + Ω g l CN T N b dV V b s (80c) R v s b l = - Ω g l CN T N b dV V b l (80d) R v s p = - Ω g s B T I + N T ∇g s N p dV P ( 80e
)
F v s = Ω g s ρ s N T g dV + ∂Ω σ N T T imp dS (80f)
-Bubble contribution:

R b s = R b s v s + R b s v l + R b s b s + R b s b l + R b s p -F b s (81) R b s v s = Ω B T b Σ s (v P 1 s ) dV + Ω g l CN T b N dV V s (81a) R b s v l = - Ω g l CN T b N dV V l (81b) R b s b s = Ω B T b Σ s (b s ) dV + Ω g l CN T b N b dV V b s (81c) R b s b l = - Ω g l CN T b N b dV V b l (81d) R b s p = - Ω g s B T b I + N b ∇g s N p dV P (81e) F b s = Ω g s ρ s N b g dV (81f)
• Momentum in the liquid phase:

-P 1 contribution:

R v l = R v l v s + R v l v l + R v l b s + R v l b l + R v l p -F v l (82) R v l v s = - Ω g l CN T N dV V s (82a) R v l v l = Ω B T Σ l (v P 1 l ) dV + Ω g l CN T N dV V l + Ω N T ρ l γ(g l , v P 1 l ) dV R γ v l v l (82b) R v l b s = - Ω g l CN T N b dV V b s (82c) R v l b l = Ω B T Σ l (b l ) dV + Ω g l CN T N b dV V b l + Ω N T ρ l γ(g l , b l ) dV R γ v l b l (82d) R v l p = - Ω g l B T I + N ∇g l N p dV P ( 82e
)
F v l = Ω g l ρ l N T g dV + ∂Ω σ N T T imp dS (82f) -Bubble contribution: R b l =R b l v s + R b l v l + R b l b s + R b l b l + R b l p -F b l (83) R b l v s = - Ω g l CN T b N dV V s (83a) R b l v l = Ω B T b Σ l (v P 1 l ) dV + Ω g l CN T b N dV V l + Ω N T b ρ l γ(g l , v P 1 l ) dV R γ b l v l (83b) R b l b s = - Ω g l CN T b N b dV V b s (83c) R b l b l = Ω B T b Σ l (b l ) dV + Ω g l CN T b N b dV V b l + Ω N T b ρ l γ(g l , b l ) dV R γ b l b l (83d) R b l p = - Ω g l B T b I + N b ∇g l N p dV P (83e) F b l = Ω g l ρ l N b g dS (83f)
• Mass in the solid-liquid mixture:

R p =R pv s + R pv l + R pb s + R pb l + R pp -F p (84) R pv s = Ω Δε tr -1 N p g s B T I + N ∇g s T dV V s (84a) R pv l = - Ω N p g l B T I + N ∇g l T dV V l (84b) R pb s = Ω Δε tr -1 N p g s B T b I + N b ∇g s T dV V b s (84c) R pb l = - Ω N p g l B T b I + N b ∇g l T dV V b l (84d) R pp = 0 (84e) F p = - Ω Δε tr ∂g s ∂t N p g dV (84f)
In the above expressions, N , N b and N p are respectively the matrices of linear velocity, bubble velocity and pressure shape functions, B, B b and B p the corresponding gradient matrices, and V s , V l and P are the vectors of nodal solid velocity, liquid velocity and pressure unknowns; further, the tensors such as Σ s , Σ l and I are mapped into vectors as usual in the finite element practice (see [START_REF] Zienkiewicz | The Finite Element Method[END_REF] for details).

Treatment of inertia and advection effects

For g s > g s cohe , the porosity of the solid matrix is relatively low, slowing the interstitial melt flow. In that case, it is licit to neglect the contributions to advection and inertia from the bubble degree of freedom:

R γ b l v l = R γ b l b l = 0. ( 85
)
Consequently, considering the previous hypothesis, the γ-term contribution reduces to

R γ v l v l = Ω N T ρ l g l ∂ l v l ∂t dV + Ω N T ρ l v l dg l dt + ∇ • (g l v l ) dV, ( 86
)
Time discretization is achieved via the fully implicit (then, unconditionally stable) backward Euler method. Let t = t 0 + Δt be the current instant, t 0 the previous instant, and Δt the time step size. Then, the time derivative of any variable ψ at time t is approximated as follows

∂ψ ∂t (x, t) ≈ ψ(x, t) -ψ(x, t 0 ) Δt , ( 87
)
being x the spatial coordinates with respect to a fixed reference frame.

In the computation of time derivatives, the relative movement between the material and spatial mesh must be taken into account. This is achieved by using the Arbitrary Lagrangian-Eulerian (ALE) method we proposed in [START_REF] Bellet | ALE method for solidification modeling[END_REF]. First, let us introduce the concept of mesh derivative (i.e., the derivative along the nodal trajectories [START_REF] Bellet | The ALE method for the numerical simulation of material forming processes[END_REF]), related to the k-phase-material derivative by the expression:

d k ψ dt = ∂ mesh ψ ∂t + (v k -v mesh ) • ∇ψ. ( 88
)
where v mesh is the mesh velocity field. In this work, we assume the mesh velocity equal to the intrinsic solid velocity, i.e.

v mesh ≡ v s . ( 89
)
The solid phase assumed to be coherent all along the analysis, this hypothesis makes the mesh become rarely degenerated, with the consequent saving of the computational time associated to re-meshing.

Solution of the mechanical problem

The discretized mechanical problem stated by the system of equations ( 80) to ( 84) is nonlinear in general, and will be solved iteratively. To this end, we make use of the Newton-Raphson technique, widely recognized as a very efficient nonlinear solver thanks to its quadratical convergence rate [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. This method is based on the first order Taylor expansion of the residual vector, say R, which can be approximated at k-th iteration as follows:

R(V (k) ) ≈ R(V (k-1) ) + KΔV = 0, ( 90
)
where V is the vector of generic unknowns and K = dR/dV is the tangent matrix.

In this work, convergence is further improved by using a line-search algorithm [START_REF] Matthies | The solution of nonlinear finite element equations[END_REF]. Actually, let us assume that ΔV defines not the increment but the search direction, in such a way that the unknown variable is updated to the k-th iteration as follows:

V (k) = V (k-1) + αΔV , ( 91
)
where the parameter α is determined by requiring the orthogonality between the residual vector and the search direction, i.e.

R(V

(k-1) + αΔV ) • ΔV = 0. ( 92
)
This scalar equation is solved using the simple regula-falsi method detailed in [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. Particularized to the nonlinear system of equations (80) to (84), the Newton-Raphson method implies to solve at each iteration the following linear equation system:

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ K v s v s K v s v l K v s b s K v s b l K v s p K v l v s K v l v l K v l b s K v l b l K v l p K b s v s K b s v l K b s b s K b s b l K b s p K b l v s K b l v l K b l b s K b l b l K b l p K pv s K pv l K pb s K pb l K pp ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ΔV s ΔV l ΔV b s ΔV b l ΔP ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = - ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ R v s R v l R b s R b l R p ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 93
)
where

K v s v s = dR v s dV s , K v s v l = dR v s dV l , K v s b s = dR v s dV b s , K v s b l = dR v s dV b l , K v s p = dR v s dP , etc. (94)
For a space of dimension d (d = 2 for 2D and axisymmetric cases), there are n n = 2d + 1 unknowns associated to the nodes and n e = 2d (bubble) unknowns associated to the elements of the spatial mesh, that is one more vectorial unknown per node and per element with respect to the one-velocity problem. This corresponds either to the additional liquid velocity field regarding the solid mechanics model [START_REF] Perchat | Mixed formulation and iterative solver for the parallelisation of a forging simulation software by mesh partitioning[END_REF][START_REF] Heinrich | Modélisation thermomécanique de la coulée continue d'acier en deux dimensions[END_REF] or the solid velocity considering the fluid dynamics work [START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF]. However, by rewriting equation (93) in the form

⎡ ⎣ K vv K vb K vp K bv K bb K bp K pv K pb K pp ⎤ ⎦ ⎡ ⎣ ΔV ΔV b ΔP ⎤ ⎦ = - ⎡ ⎣ R v R b R p ⎤ ⎦ , ( 95
)
we realize that, despite the supplementary unknowns, it remains identical in form to that of one-velocity problems. Furthermore, since the bubble unknowns are internal to the element, they can be eliminated at the element level by the same procedure, known as static condensation, used for one-velocity problems [START_REF] Perchat | Mixed formulation and iterative solver for the parallelisation of a forging simulation software by mesh partitioning[END_REF][START_REF] Heinrich | Modélisation thermomécanique de la coulée continue d'acier en deux dimensions[END_REF]. From equation (95) applied to each element, the iterative bubble increment results

ΔV b = -K -1 bb (R b + K bv ΔV + K bp ΔP ) . ( 96
)
Then, the system (95) reduces to

K vv -K vb K -1 bb K bv K vp -K vb K -1 bb K bp K pv -K pb K -1 bb K bv K pp -K pb K -1 bb K bp ΔV ΔP = K vb K -1 bb R b -R v K pb K -1 bb R b -R p . ( 97
)
The resulting linear and symmetric system of equations is solved using the PETSc library [START_REF] Balay | PETSc Users Manual[END_REF], which allows parallel processing.

Transport of the solution in case of moving meshes

By solving the mechanical problem defined by the equation system (97), we obtain an estimate to the liquid-phase-material derivative of the velocity in the liquid phase:

d l v l dt ≈ v l (x + v l Δt, t) -v l (x, t 0 ) Δt , ( 98
)
from which an estimate to the mesh derivative is computed:

∂ mesh v l ∂t = d l v l dt -(v l -v mesh ) • ∇v l . ( 99
)
Finally, the time-discretized form of the mesh derivative,

∂ mesh v l ∂t ≈ v l (x + v mesh Δt, t) -v l (x, t 0 ) Δt , ( 100
)
determines the updated value of the liquid velocity:

v l (x + v mesh Δt, t) = v l (x, t 0 ) + ∂ mesh v l ∂t Δt. ( 101
)
The solutions for the velocity in solid phase as well as the interstitial pressure can be updated in the same way, although it is not strictly necessary. In fact, we do not need to know the previous values of the solid velocity due to the hypotheses of negligible inertia effects in the solid phase (see Section 3.1.2), nor those of the pressure, since no initial conditions appear for this variable.

Thermal model 3.2.1. Treatment of the mesh movement

As done for the mechanical problem, an ALE method allowing the description of the relative movement between the solid and the liquid regions and the mesh must be used. However, as demonstrated by Gouttebroze et al [START_REF] Gouttebroze | 3D-FEM modeling of macrosegregation in solidification of binary alloys[END_REF], the performance of the Streamline-Upwind Petrov/Galekin (SUPG [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF]) for solving the heat balance equation is superior to that of the ALE method previously described for the mechanical problem. In order to apply the SUPG formulation, we introduce the mesh derivative of the average enthalpy

∂ mesh h ∂t = ∂ h ∂t + v mesh • ∇ h , (102)
into the heat balance equation ( 46), obtaining

ρ ∂ mesh h ∂t + ρ ṽ • ∇ h -∇ • (κ∇ h ) + S T = 0, (103)
where

ρ ṽ = ρv -ρ v mesh . ( 104
)
As before, equation ( 103) is discretized in time using the backward-Euler finite-difference scheme, being the mesh derivative approximated as follows

∂ mesh h ∂t ≈ h (x + v mesh Δt, t) -h (x, t 0 ) Δt , (105)
Then, the solution of equation ( 103) directly yields the updated value of the average enthalpy at the new nodal positions x + v mesh Δt, and, unlike mechanical analysis, no transport step is needed.

Finite element formulation of the thermal problem

The solution of the thermal problem consists in finding h ∈ L 2 (Ω) satisfying:

Ω ρ ∂ mesh h ∂t ω dV + Ω ρ ṽ • ∇ h ω dV + Ω κ∇ h • ∇ω dV + Ω S T ω dV - ∂Ω q q w ω dS - ∂Ω c h(T -T ext )ω dS = 0, ( 106
)
for all ω ∈ L 2 (Ω) such that ω = 0 on ∂Ω T .

As usual, the solution is approximated inside each linear triangular finite element as

h ≈ 3 i=1 N i (x)H i , ( 107
)
where H i is the unknown value of the average enthalpy at node i located at the vertex i of each element, N i is the standard linear shape function associated to the node i, the same as used for the mechanical problem.

However, unlike the mechanical problem, the presence of the advection term (second term in the l.h.s. of the above equation) may induce numerical instabilities when the standard Galerkin formulation is used. These defects take the form of local node-to-node oscillations (or "wiggles") in presence of a high temperature gradient as downstream boundary condition. In that case, the advection term must be stabilized. In this work, we make use of the widely accepted SUPG formulation, originally proposed by Brooks and Hughes [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF], where the weighting function (different from the shape function, as meant by "Petrov/Galerkin") is defined in such a way to add a stabilizing diffusivity in the streamline direction. For a typical element, we have:

ω = N i + α(ṽ • ∇N i ) = Ñi , ( 108
)
where ṽ is evaluated at the center of the element and α is the upwind coefficient, defined as [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] α = Δx √ 15 ṽ

1 tanhP e - 1 P e . ( 109
)
where * denotes the L 2 -norm of ( * ),

P e = ṽ Δx 2κ/ ρ (110)
is the element Peclet number, and Δx is a characteristic size of the element in the advection direction, defined as [START_REF] Tezduyar | Stabilization parameters in SUPG and PSPG formulations[END_REF] Δx = 2 ṽ

3 i=1 |ṽ • ∇N i | . ( 111
)
Introducing the trial and weighting functions defined by equations ( 107) and (108) into the weak form (106) of the heat balance, we obtain the discrete heat equation:

R T = M ∂ mesh H ∂t + (C + K)H + S(H) + c(H) -F = 0, ( 112
)
the components of the above arrays being given by:

M ij = Ω ρ Ñi N j dV, (113)
C ij = Ω Ñi ( ρ ṽ • ∇N j ) dV, (114)
K ij = Ω κ∇ Ñi • ∇N j dV, ( 115
)
S i = Ω S T ( h ) Ñi dV, ( 116
)
c i = ∂Ω c hT ( h )N i dS, ( 117
)
F i = ∂Ω q q w N i dS + ∂Ω c hT ext N i dS. ( 118
)
In the above, the existence of a function T = T ( h ), reciprocal to that given by the enthalpy definition [START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF], is assumed. Further, for a P 1 element, ∇ Ñi ≡ ∇N i . The nonlinear and non-symmetric heat equation is solved using the Newton-Raphson technique together with line-searching, in a way analogous to that described for the mechanical problem in Section 3.1.3. Once again, at each Newton-Raphson iteration, linear solvers from the PETSc library [START_REF] Balay | PETSc Users Manual[END_REF] are used.

Macrosegregation model

Since an analogy can be made between energy and solute balance equations, [START_REF] Bellet | The ALE method for the numerical simulation of material forming processes[END_REF] and (66) respectively, the latter is solved in an identical way to that described in the previous section.

First, by introducing the mesh derivative of the average concentration, in the discrete form

∂ mesh w ∂t ≈ w (x + v mesh Δt, t) -w (x, t 0 ) Δt , ( 119
)
into the solute balance equation (66), we have

∂ mesh w ∂t + ṽw • ∇ w -∇ • D∇ w + S w = 0, ( 120
) with ṽw = v w -v mesh . ( 121
)
It is worth noting here that taking v mesh = v s , the solid intrinsic velocity, the time derivative in equation ( 120) appears as the rate of variation of the average concentration when following the displacement of the solid phase. Equation ( 120) is then very close to the equation obtained by Lesoult et al [START_REF] Lesoult | Segregation during solidification with spongy deformation of the mushy zone[END_REF], who developed a macrosegregation equation during the deformation of a coherent mushy zone using a Lagrangian approach.

Then, proceeding as before, the trial function is defined as a linear combination of P 1 shape functions:

w ≈ 3 i=1 N i (x)W i , ( 122
)
being W i the unknown value at node i. Further, a SUPG weighting function is used, identical to that defined by equation ( 108) with

α = Δx √ 15 ṽw 1 tanhP e - 1 P e , ( 123
) P e = ṽw Δx 2 D , ( 124
)
Δx = 2 ṽw 3 i=1 |ṽ w • ∇N i | . ( 125
)
Now, the discrete solute balance equation can be expressed as

R w = M w ∂ mesh W ∂t + (C w + K w )W + S w (W ) = 0, ( 126
)
where

M w ij = Ω Ñi N j dV, ( 127
)
C w ij = Ω Ñi (ṽ w • ∇N j ) dV, (128)
K w ij = Ω D∇N i • ∇N j dV, ( 129
)
S w i = Ω S w ( w ) Ñi dV. ( 130
)
Finally, equation ( 126) is solved here again using the Newton-Raphson method with linesearching, using at each iteration the linear solvers from the PETSc library.

Coupled solution of thermal, mechanical, and macrosegregation problems

Intermediate thermo-macrosegregation coupling

Thermal and macrosegregation solutions interact due to the simultaneous influence of both enthalpy and solute concentration on the local liquid or solid fraction. Assuming h and w to be known from the solution of thermal and macrosegregation problems, respectively, we invoke the enthalpy dependency on g l and T defined by equation [START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF], the relationship (61) between w l and T derived from the equilibrium phase diagram, and the generic microsegregation rule (65) relating relationship w to w l and g l , in order to pose a system of scalar equations from which w l , g l and T can be computed. The so-defined "microsegregation box" can be expressed as follows:

m(w l , T, g l ; h , w ) = ⎡ ⎣ h - T T 0 c p (τ ) dτ -g l L T -T m -m l w l w -ϕ(g l , w l ) ⎤ ⎦ = 0. ( 131
)
Henceforth, by writing ( * ; •), it is implied that ( * ) are variables and (•) are parameters.

This system has a closed analytical solution just in a few particular cases, for instance, when the microsegregation is described by the lever rule (62) and the heat capacity c p is assumed to be constant. Otherwise, it is solved numerically using a regula-falsi method, spanning between g l = 0 at the solidus temperature T sol (or g l = g l eut at the beginning of the eutectic transformation at temperature T eut ) and g l = 1 at the liquidus temperature T liq .

Above system also determines the derivatives ∂T /∂ h and ∂w l /∂ w . As before, these derivatives have analytical expressions for a few cases, and are in general computed numerically.

Let us remark that the system (131), containing three scalar unknowns, has to be solved locally for each sampling point of the mesh at each Newton-Raphson iteration of the thermal and macrosegregation analysis.

Global coupling scheme

Besides the interaction between them, thermal and macrosegregation models are also coupled to the mechanical model. On the one hand, mechanical properties are in general dependent on temperature, and temperature and solute gradients induce natural convection. On the other hand, the velocities computed from the mechanical analysis define the advection velocities ṽ and ṽw in the energy and solute balance equations, respectively.

The coupled solution of all the balance equations is achieved using a simple staggered scheme, where these equations are solved successively at each time step, as described as follows:

• New time instant t ← t + Δt:

-Initialization: 1) >tolerance, go to the next coupling iteration. * Otherwise, update the variables

k = 0, (v s , v l , p l , h , w ; T, w l , g l ) (0) t+Δt = (v s , v l , p l , h , w ; T, w l , g l ) t -New coupling iteration k ← k + 1: 1. Obtain h (k) by solving the thermal problem R T ( h (k) ; w (k-1) , v (k-1) s , v (k-1) l ) = 0.

Obtain w (k) by solving the macrosegregation problem

R w ( w (k) ; h (k) , v (k-1) s , v (k-1) l ) = 0.

Obtain w

m(w (k) l , T (k) , g (k) l ; h (k) , w (k) ) = 0.

Obtain v

R M (v (k) s , v (k) l , p (k) l ; w (k) l , T (k) , g (k) l ) = 0. 5. Check convergence: * If ( * ) (k) -( * ) (k-
(v s , v l , p l , h , w ; T, w l , g l ) t+Δt = (v s , v l , p l , h , w ; T, w l , g l ) (k) t+Δt
Then, go to the next time step.

In the convergence checking, ( * ) represents any or some of the involved variables. However, in this work the time step is assumed to be small enough, and only one coupling iteration per time step is performed, making that checking superfluous.

In addition, it is worth noting that the "microsegregation box" is also called inside the thermal and macrosegregation solvers at each Newton-Raphson iteration, not only for computing g l , T , and w l , but also the derivatives ∂T /∂ h and ∂w l /∂ w , as mentioned in Section 4.1.

Finally, let us remark the different sizes of the problems to solve. Being N the total number of nodes of the mesh, the solution of the thermal and macrosegregation problems requires to solve systems of size N , while the system defined by mechanical problem contains 5N unknowns for 2D and axisymmetric analysis. On the other hand, the system to solve within the "microsegregation box" has only 3 scalar unknowns, and has to be solved at each sampling or Gauss point of the mesh. Usually, there is only one sampling point per element, making the computational cost associated to this operation marginal compared to the cost of solving the global thermal, macrosegregation and mechanical systems for typical meshes containing thousands of nodes.

Application tests

Validation: Axial compression of a semi-solid specimen

Let us consider a simple compression test on a porous specimen saturated with liquid, illustrated in Figure 3. This problem has a closed-form solution discussed in the Appendix, and serves then for the validation of the finite element implementation of the mechanical model previously presented. The following application aims to be representative of the mechanical conditions encountered in the mushy core of a continuously cast steel slab at the end of the secondary cooling stage when passing along the support rolls. The values of the considered geometrical and material parameters are given in Table I. as it varies due to a variation of the secondary dendrite arm spacing λ 2 . As observed, the lower this spacing (and the permeability as well), the slower the interdendritic flow (characterized Figure 4 by the relative velocity v lv s ) and the higher the interstitial pressure p l .

The results of the two-phase finite element model match exactly the closed form solution, as illustrated in Figure 5 for v lv s and p l .

Application to liquid segregation

The second test aims to simulate a situation frequently encountered in continuous casting problems. It consists in imposing a deformation on a reduced portion of the surface of a partially solidified alloy (Figure 6). It thus simulates the action of a roll. Symmetry conditions are imposed on all the boundaries except on the upper face, which is a free surface. We perform here again an isothermal calculation, but imposing an initial distribution of liquid fraction along the y-axis described in Figure 6. Here, the solid phase rheology is still somewhat arbitrary, but it is now non-linear, with a strain rate sensitivity m = 0.2. The interaction coefficient is modelled by a classical Carman-Kozeny law and reaches a very high value in the solidified zone. Black arrows plotted on the deformed domain on the right Figure 6 show the relative average velocity vectors v lg l v s . It can be seen that the applied pressure leads to a redistribution of the liquid with respect to the solid. In the zone where the solid is under pressure (left side of the deformed sample), it undergoes densification and the liquid is rejected to the right side of the sample.

Application to the late stage of solidification in steel continuous casting

Test case I

Macrosegregation has always been one of the major defects we have to face up in the development of steel continuous casting processes. This defect is influenced by the bulging of the strand between the support rolls, which produces tension in the core of the strand that is not completely solidified, inducing the inflow of the molten liquid rich in segregated solutes. Casting engineers have developed several techniques to control this phenomenon, acting either on heat transfer (intense local cooling of the shell) or by mechanical means ("soft" reduction of the strand width).

Nowadays, the understanding and control of macrosegregation is particularly pertinent for high-strength steels, where the high content of alloying elements promotes the formation of this defect, as well as that of hot-cracking, also associated with the development of tensile stresses at the end of solidification.

This problem has already been addressed by several authors. Following the pioneer work of Miyazawa and Schwerdtfeger [START_REF] Miyazawa | Macrosegregation in continuously cast steel slabs: preliminary theoretical investigation on the effect of steady state bulging[END_REF], and the analytical developments of Lesoult and Sella [START_REF] Lesoult | Analysis and prevention of centreline segregation during continuous casting of steel related to deformation of the solid phase[END_REF], Kajitani et al. [START_REF] Kajitani | Numerical simulation of deformation induced segregation in continuous casting of steel[END_REF] proposed the more advanced model. However, they did not take into account the deformation of the solid phase in the mushy zone in their approach, and made strong hypotheses on its velocity field, which was deduced directly and homothetically from the velocity field of the solid crust. The present two-phase formulation aims to improve the modelling of macrosegregation by considering not only the relative but also the concurrent movement between the solid and the liquid phases.

Test case I consists of a preliminary study of continuous casting of steel slabs in the late stage of solidification. A typical casting velocity V cast = 1.25 m/min is imposed, and a uniform time step Δt = 1 s is adopted all along the computation.

Steel is modelled as a binary iron-carbon alloy, whose main properties are listed in Table II. However, in this preliminary study, the flow stress of the material has been artificially decreased. For this purpose, the viscoplastic consistency K has been divided by a factor 10. This has been done intentionally, in order to magnify the bulging of the slab between the supporting rolls, which is supposed to be the basic phenomenon giving rise to mushy zone deformation and macrosegregation. This will permit a better visualization of the interstitial melt flow, providing a concrete illustration of the potential and interest of the two-phase thermomechanical approach. In a second step, the nominal values of K will be used (see next section).

Further concessions must be made due to the constraints inherent to the model. For instance, the hypothesis of coherent solid phase has required to start the two-phase analysis at a transverse section, say the initial two-phase section, far enough from the meniscus. Besides, the current constitutive model for the solid phase in two-phase analysis does not include hardening. Finally, solidification shrinkage and thermal strains have not been accounted for, by taking ρ s = ρ l equal and constant.

Up to this two-phase section, a classical "one-phase" analysis [START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF] is performed, which is used as an initialization for the two-phase analysis. In Figure 7 we observe thermal results from the "one-phase" analysis, particularly concerning the liquid fraction distribution along the strand. The enthalpy distribution at the initial two-phase section (located 15.2 m-far from the meniscus) serves as initial condition for the two-phase thermal analysis. Since the solid fraction must be high enough to let us assume the solid phase to be coherent, we eliminate the higher values of enthalpy at the core of the slab. Further, the current version of the "one-phase" solver [START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF] is not able up to date to perform a macrosegregation analysis under continuous casting conditions, so we will assume a uniform solute concentration w 0 (equal to the nominal carbon content) in the initial two-phase section as initial condition. Although it implies an important simplification, this hypothesis allows us to distinguish the effect of solid deformation from the other sources of macrosegregation that are relevant in the earlier stages of the continuous casting process (for instance, macrosegregation by natural convection driven by thermal and solutal gradients).

In order to use the global unsteady-state approach [START_REF] Heinrich | Modélisation thermomécanique de la coulée continue d'acier en deux dimensions[END_REF][START_REF] Bellet | A two-dimensional finite element thermomechanical approach to a global stress-strain analysis of steel continuous casting[END_REF], a buffer zone is defined, consisting of a thick slice of the slab that ends at the initial two-phase section, as shown in Figure 8. Unlike in the one-phase analysis, where the buffer zone is located above the meniscus in the mold region, the buffer zone now lies in the curved portion of the caster. As before, the buffer length must be greater than the maximum displacement per time step of any particle in the initial two-phase section, otherwise it could be impossible to find the particle antecedent in the advection calculations.

In the buffer zone, the enthalpy and solute distributions are assumed not to vary in the longitudinal (i.e. casting) direction, remaining equal to those of the initial two-phase section. Concerning the mechanical analysis, the metallostatic pressure is applied to the upper face of the buffer zone (H in Figure 8 denotes the metallostatic head), while at the current lower face a fictitious extracting tool is modelled by imposing an uniform velocity v imp in the axis direction, with v imp = V cast S/S 0 , being S 0 the area of the meniscus section and S the current lower face area.

The two-phase analysis stops at a casting distance of 2.20 m, after exceeding the maximum number of finite elements admitted by the current version of R2SOL (50000).

Figure 9 plots the region with positive and negative values of the divergence of the intrinsic velocity in the solid phase (∇•v s ), indicating respectively expansion and contraction of the solid skeleton in the mushy zone. As it can be seen there, the solid phase is subject to alternating expansion and contraction at the core of the slab. Expansion usually develops in the intervals between two successive rolls, and is associated to the local influx of liquid melt caused by the bulging of the slab. Conversely, contraction normally develops under the rolls, due the local expulsion of liquid. Figure 10 evidences the complexity of the melt flow in the mushy region. Between two successive rolls, the bulging of the solidified shell has two main effects: the development of tensions in the mushy core of the slab, and concurrently, of compressions in the periphery of this region. The latter effect can also be explained by assimilating the longitudinal section of the solidified shell as a beam supported by the rolls, whose fibers adjacent to the solidification front suffers compression. This combined central expansion-peripheral compression causes the relative velocity v lv s point to the core, implying the liquid enrichment of this region.

On the other hand, when passing through the rolls, the mushy core is compressed, while its periphery (in contact with the solid shell) is under tension. Then, the liquid in the region under the rolls is expelled from the mushy core. This behavior is in complete agreement with the distribution of positive and negative ∇ • v s , as evidenced in Figure 11, where Figures 9 and 10 are superimposed.

Therefore, each time the slab passes through a roll pair, there is a redistribution of interstitial liquid within the mushy zone. This liquid, rich in segregated solute (carbon), logically contributes to the formation of macrosegregations, as shown by Figure 12.

The axial shift between the zones of positive ∇ • v s and those of positive segregations is probably due to a delay which corresponds to the transport of solute, taking into account the difference between the solid and the liquid velocity and the general casting velocity. Indeed, the solute transport is not instantaneous. Further, major macrosegregation defects are shown to appear at the center of the slab, with a slighter correlate of opposite sign in its periphery. In real continuous casting processes, only positive central macrosegregation is expected. The spurious alternate sign of central macrosegregation observed here may be explained by the excessive Let us also notice the existence of some oscillations of the numerical solution for the solute distribution, particularly in the periphery of the mushy core. In fact, as shown by Figure 13, the mesh density decreases from the core, where the mesh is adequately refined as proved by the smoothness of the solution, to the outside, where it is evidently insufficient. This problem should be overcome when the new adaptive mesh refinement tools developed by Liu [START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF] in R2SOL for one-phase analysis be operational for the current two-phase analysis. This would allow to keep an adequate refinement of the continuously evolving mushy region at each time step.

Solidus isotherms

To conclude with this test, we recall that this first two-phase application to continuous casting just aims to give a qualitative representation of the phenomena associated to the concurrent solid deformation and liquid flow. Nevertheless, this application has served to highlight the great potentiality of the two-phase modelling.

Test case II

The second application intends to assume more realistic bulging conditions in order to get closer to the real physics of a continuously cast steel slab, although we must keep in mind the limitations imposed by the adopted hypotheses, particularly the coherency of the solid phase and the lack of hardening effect in the constitutive equation for the solid phase.

The material to be considered is a Fe-0.18wt%C steel. Material parameters are listed in Table III. Slabs with a thickness of 222 mm are cast in a vertical-curved machine whose curvature radius before unbending is 12.5 m. A lower casting velocity, equal to 0.86 m/min, promotes the formation of a wider mushy zone earlier in the caster. Subsequently, the initial two-phase section is located now at a metallurgical length of 11 m, and there is no need to cut the initial enthalpy curve, as shown in Figure 14.

By using a slightly coarser mesh and increasing the maximum number of finite elements allowed by R2SOL to 80000, the analysis reached a casting distance of 3.25 m, as shown in Figure 15. This figure also plots on the right top the initial liquid fraction, highlighting an important problem concerning the coupling of one-phase and two-phase models. The one-phase model used for initialize the two-phase analysis [START_REF] Heinrich | Modélisation thermomécanique de la coulée continue d'acier en deux dimensions[END_REF] is mainly devoted to the thermo-mechanical analysis in the solidified shell, the mushy and liquid zones being treated in an approximative way, with rather few elements. Due to the huge computational time and memory requirements, it is nowadays unfeasible to perform a one-phase analysis up to the initial two-phase section with a mushy zone as refined as the two-phase analysis requires. This problem is open and should be the matter of future developments.

Regarding macrosegregation, the nominal carbon content w 0 = 0.18% is set as initial condition inside the buffer zone. The final carbon segregation ww 0 is depicted in Figure 16. We observe there the ability of the current model to capture the phenomenon of central Even if the solidification is not completed so that the solute concentrations can still evolve, the computed level of macrosegregation seems to be significantly lower than that observed in reality. Nevertheless, we feel that this could be strongly improved by the use of more realistic constitutive models and parameters, especially for the solid phase compressible behavior and the mechanical interaction term. Furthermore, as already mentioned, a finer mesh would be necessary to confirm such trends.

Conclusions and perspectives

This work contributes to the understanding of the complex phenomena observed in a solidifying medium, particularly in the late stages of solidification once the solid phase has already developed a coherent structure. This assumption was thought to be adequate to characterize the secondary cooling region of continuous casting processes, when the cast strand is subject to the alternate effects of rolling and bulging, which induces alternate compression and expansion states in the mushy core of the strand. The solid skeleton in the mushy region is seen as a viscoplastic compressible body and can be assimilated to a sponge that absorbs or expels the interstitial liquid, rich in segregated solutes. This mechanism is the main source of macrosegregation in the secondary cooling zone, and its responsibility on central macrosegregation in continuously cast steel slab has been demonstrated in this work.

However, further research is needed before obtaining a robust computational tool with the accuracy required in real casting applications. First, although the viscoplastic compressible model [START_REF] Abouaf | Finite element simulation of hot isostatic pressing of metal powders[END_REF][START_REF] Nguyen | Theoretical and experimental study of the isothermal mechanical behaviour of alloys in the semi-solid state[END_REF][START_REF] Martin | Viscoplastic behaviour of porous metallic materials saturated with liquid. part i: Constitutive equations[END_REF] seems to be adequate to characterize the macroscopic mechanical behavior of the coherent solid phase in the mushy zone, this is contested by the lack of experimental (and literature) data. Further, this model does not account for strain hardening, a mechanism that becomes important once the material is completely solidified. It is not only isotropic hardening that must be modelled, as done in "one-phase" models, but also kinematic hardening in order to account for alternate loading.

Concerning the viscoplastic incompressibility constraint in the fully solid region, the current two-phase model enforces it by means of a penalty factor for the volume change rate tr ( εs ) into the constitutive equation for the solid phase, while in the "one-phase" model [START_REF] Heinrich | Modélisation thermomécanique de la coulée continue d'acier en deux dimensions[END_REF][START_REF] Bellet | A two-dimensional finite element thermomechanical approach to a global stress-strain analysis of steel continuous casting[END_REF] it is implied in the mass conservation equation. This difference blocks the use of the available "onephase" mechanical model [START_REF] Heinrich | Modélisation thermomécanique de la coulée continue d'acier en deux dimensions[END_REF][START_REF] Bellet | A two-dimensional finite element thermomechanical approach to a global stress-strain analysis of steel continuous casting[END_REF], which is already capable of modelling isotropic hardening, for the fully solid region in two-phase analysis. The development of a library of constitutive models, suitable for one-phase and two-phase mechanical analysis, should be encouraged as future work.

A last remark concerns mesh discretization: a quite fine mesh is needed for an accurate description of the mushy zone. In the current version of R2SOL, the mesh is refined by regions defined a priori in the buffer zone. Then, as the solidification progresses and the width of the mushy zone decreases, there is an excessively fine mesh outside this zone. This suggests that the computational cost could be reduced by using adaptive remeshing, i.e., by reducing the mesh size only where it is necessary. This is a work in progress [START_REF] Liu | Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes[END_REF].

Acknowledgements

The financial support of ARCELOR Research, ASCOMETAL (from LUCCHINI group) and the French Ministère de l'Economie, des Finances et de l'Industrie, in the frame of the OSC-Continuous Casting project, is acknowledged. Victor D. Fachinotti is also granted by the Argentine Council for Scientific and Technical Research (CONICET).

APPENDIX

Analytical solution of a simple compression test on a porous specimen

Let us consider the compression test on a porous specimen saturated with liquid numerically solved in Section 5.1. For convenience, let us assume a linear constitutive law for the solid phase behavior, i.e. m = 1 in equation ( 21), yielding

Σ s = K s (β ε s + α tr ( ε s )I) , ( 132
)
where

β = 3 A , α = 1 3B - 1 A . ( 133
)
The liquid behavior is described by equation ( 15) [START_REF] Ganesan | Conservation of mass and momentum for the flow of interdendritic liquid during solidification[END_REF][START_REF] Rappaz | Numerical modelling in materials science and engineering[END_REF], linear as well, rewritten here as Now, let us consider the simple compression test illustrated in Figure 3, adopting the following assumptions:

Σ l =
• plane strain conditions;

• the liquid fraction is uniform throughout the specimen; • gravity and inertia effects are neglected; • no mass transfer is considered between the two phases;

• the normal stress on each phase is supposed to be null at x = ±l.

We start with the following kinematic hypothesis:

v s (x) = u(x) -εy , v l (x) = v(x) -εy , ( 136
)
where ε = v imp /h denotes the nominal compression strain rate. The strain rate tensors in both phases have then the following expression:

ε s = ε(v s ) = du dx 0 0 -ε , ε l = ε(v l ) = dv dx 0 0 -ε , ( 137
)
Applying momentum and mass conservation equations, we get:

∇ • Σ s -g s ∇p l -M d l = 0, (138) ∇ • Σ l -g l ∇p l + M d l = 0, (139) ∇ • (g s v s ) + ∇ • (g l v l ) = 0, (140)
with

M d l = - g 2 l μ l κ (v l -v s ). ( 141
)
Replacing the constitutive equations ( 132) and (134) as well as kinematic hypothesis ( 136) and (137) into the above conservation equations, we obtain

K s (α + β) d 2 u dx 2 -g s ∂p l ∂x + g 2 l μ l κ (v -u) = 0, ( 142
)
μ l g l (a + b) d 2 v dx 2 -g l ∂p l ∂x - g 2 l μ l κ (v -u) = 0, ( 143
)
∂p l ∂y = 0 ⇒ p l = p l (x), (144)

Figure 1 .

 1 Figure 1. Local thermodynamic equilibrium at the interface, governed by the phase diagram (on the right), where both solidus and liquidus lines are linearized.

Figure 2 .

 2 Figure 2. Triangular P 1 + /P 1 finite element.

l

  , T(k) , and g

  calling the "microsegregation box":

  by solving the mechanical problem

l = 0.25 m K s = 10 7 8 Figure 4 = λ 2 2 g 3 l180( 1 -g l ) 2 Figure 3 .

 7843123 Figure4shows the influence of the permeabilityκ = λ 2 2 g 3 l 180(1g l ) 2

Figure 4 .

 4 Figure 4. Analytical solutions for the relative velocity v lv s and the interstitial pressure p l showing the influence of the secondary dendrite arm spacing λ 2 .

Figure 5 .Figure 6 .

 56 Figure 5. Analytical vs. two-phase finite element solution for the relative velocity v lv s and the interstitial pressure p l (λ 2 = 150 μm).

Figure 7 .Figure 8 .

 78 Figure 7. Test case I. One-phase analysis of steel continuous casting process. On the left, distribution of liquid fraction along the slab. On the top, the enthalpy distribution at the initial two-phase section serving as initial thermal condition for the two-phase analysis.

Figure 9 .

 9 Figure 9. Test case I. Final distribution of the negative (blue or dark grey) and positive (red or light grey) values of the divergence of the intrinsic velocity in the solid phase, ∇ • v s .

Figure 10 .

 10 Figure 10. Test case I. Detail of the distribution of relative velocity vectors (v lv s ) in the mushy zone.

Figure 11 .Figure 12 .Figure 13 .

 111213 Figure 11. Test case I. Effect of the contraction (blue or dark grey) and the expansion (red or light grey) of the solid skeleton on the relative velocity v lv in the mushy zone.

Figure 14 .Figure 15 .

 1415 Figure 14. Test case II. One-phase analysis of steel continuous casting process. On the left, distribution of liquid fraction along the slab. On the top, the enthalpy distribution at the initial two-phase section serving as initial thermal condition for the two-phase analysis.

Figure 16 .

 16 Figure 16. Test case II. Final distribution of carbon segregation ( w -w 0 ) in the analyzed portion of the slab.

  g l μ l b ε l + a tr ( ε l )I ,

  momentum conservation equations (142) and (143) and using then the derivative of the mass balance equation (145) with respect to x, we have∂p l ∂x = (K s (α + β)g s μ l (a + b)) d 2 u dx 2 , (146)which is integrated to givep lp 0 = (K s (α + β)g s μ l (a + b)) du dx . (147)The integration constant p 0 is determined by the free surface conditions at x = l:σ s xx = Σ s xxg s p l = 0, (148)σ l xx = Σ l xxg l p l = 0. (149)By introducing the constitutive laws (132) and (134), the mass balance equation (145), and equation (147) into the previous boundary conditions, after a few operation we getμ du dx x=l -αK s ε = g s p 0 , (150) μ du dx x=l + μ l ε(ag s + b) = g l p 0 , (151) with μ = g l K s (α + β) + g 2 s μ l (a + b). (152)Now, the sum of equations (150) and (151) yieldsp 0 = (μ l (ag s + b) -αK s ) ε. (153)On the other hand, injecting equation (146) together with the integral form of the mass balance equation (145) given byg s u(x) + g l v(x)εx = 0,(154)into the solid momentum equation (142), we arrive toμ d 2 u dx 2 + g l μ l κ ( εxu) = 0,(155)which, subject to the boundary condition u(0) = 0, has the following solutionu(x) = εx + a 0 sinh(rx), (156)wherer = g l μ l κμ . (157)We determine the integration constant a 0 by replacing this expression for u into equation (150) and using the expression (153) for p 0 . This leads toa 0 = g l (bg s μ l -βK s ) μr cosh(rL) ε.(158)Finally, from the integral mass balance (154), we getv(x) = εx -g s g l a 0 sinh(rx). (159)Summarizing, the full solution of this problem is given byv s x (x) = εx + a 0 sinh(rx), v l x (x) = εx -g s g l a 0 sinh(rx), p l (x) = p 0 + (K s (α + β)g s μ l (a + b)) ( ε + a 0 r cosh(rx)) ,with r = g l μ l κμ , μ = g l K s (α + β) + g 2 s μ l (a + b), a 0 = g l (bg s μ l -βK s ) μr cosh(rL) ε, p 0 = (μ l (ag s + b) -αK s ) ε.

Table I .

 I Dimensions and mechanical parameters used in the validation test.

Table II .

 II Material data for continuous casting test case I.

	Alloy	Fe-0.1%C
	ρ s = ρ l	7060. kg/m 3
	c p	675. J/(kg K)
	λ s = λ l	30. W/(m K)
	L	2.6×10 5 J/kg
	T m	1538. o C
	m l	-80. K/wt%C
	k	0.18
	λ 2	1.×10 -4 m
	Microsegregation law Lever rule
	Solid phase	Compressible viscoplastic [12]
	K(T ), m(T )	from Kozlowski et al [55]
		(K(T ) intentionally decreased)
	A(g s ), B(g s )	from Nguyen et al [12]
	g s cohe	0.65
	Liquid phase	Newtonian
	μ l	0.001 Pa s

Table III .

 III Material data for continuous casting test case II.

	Alloy	Fe-0.18%C
	ρ s = ρ l	7300. kg/m 3
	c p	675. J/(kg K)
	λ s = λ l	35. W/(m K)
	L	2.6×10 5 J/kg
	T m	1538. o C
	m l	-80. K/wt%C
	k	0.18
	λ 2	1.×10 -4 m
	Microsegregation law Lever rule
	Solid phase	Compressible viscoplastic [12]
	K(T ), m(T )	from Kozlowski et al [55]
	A(g s ), B(g s )	from Nguyen et al [12]
	g s cohe	0.65
	Liquid phase	Newtonian
	μ l	0.001 Pa s

(positive) segregation, characteristic of continuously cast steel slabs

[START_REF] Aboutalebi | Coupled turbulent flow, heat, and solute transport in continuous casting processes[END_REF]

.