Thi-Thuy-My Nguyen
email: nguyenttmy@gmail.com

Charles-André Gandin

Hervé Combeau
email: cherve.combeau@univ-lorraine.frdmiha.zaloznik@univ-lorraine.fremichel.bellet@mines-paristech.fr

Miha Založnik

Michel Bellet

Finite element multi-scale modelling of chemical segregation in steel solidification taking into account the transport of equiaxed grains

Keywords: Solidification, Modeling, Grain structure, Steel, Macrosegregation, Finite Element Method Solidification, Modeling, Grain structure, Steel, Macrosegregation, Finite Element Method

The transport of solid crystals in the liquid pool during solidification of large ingots is known to have a significant effect on their final grain structure and macrosegregation. Numerical modeling of the associated physics is challenging since complex and strong interactions between heat and mass transfer at the microscopic and macroscopic scales must be taken into account. The paper presents a finite element multi-scale solidification model coupling nucleation, growth and solute diffusion at the microscopic scale, represented by a single unique grain, while also including transport of the liquid and solid phases at the macroscopic scale of the ingots. The numerical resolution is based on a splitting method which sequentially describes the evolution and interaction of quantities into a transport and a growth stage. This splitting method reduces the nonlinear complexity of the set of equations and is, for the first time, implemented using the finite element method. This is possible due to the introduction of an artificial diffusion in all conservation equations solved by the finite element method.

Simulations with and without grain transport are compared to demonstrate the impact of solid phase transport on the solidification process as well as the formation of macrosegregation in a binary alloy (Sn-5wt%Pb). The model is also applied to the solidification of the binary alloy Fe-0.36wt%C in a domain representative of a 3.3-ton steel ingot.

Introduction

The casting industry commonly faces difficulties in the production of products free from macrosegregation [1]. Microsegregation naturally takes place during solidification as a redistribution of the chemical species at the solid-liquid interface occurs, the result of thermodynamic equilibrium. Various processes take place which enhance segregation far from its interfacial origin including long range diffusion, liquid flow due to shrinkage, melt convection, solid deformation and transport of solid grains/fragments in the casting. The later induce macrosegregation from both solid and liquid transport over long distances. While these phenomena are difficult to avoid during conventional casting of metallic alloys, their magnitude needs to be controlled. The local average composition defined in a small representative volume must not deviate from the nominal composition of the alloy by more than a few percent, otherwise the properties could vary significantly and subsequent thermomechanical heat treatments may not be able to restore the desired properties everywhere in the product. This is true for various classes of metallic alloys, including large steel products [2]. Recently, efforts have been made to provide detailed experimental characterizations of macrosegregation in large steel ingots [3,4,5].

Numerical modeling of solidification accounting for the transport of the equiaxed grains at the process scale remains limited. Work was first done by Ni and Beckermann who proposed a volume-averaged model that consistently coupled microscopic phenomena with macroscopic transport [6,7]. Other solidification models, based on similar principles, have also been developed [8][9][10][11][12][13][14][15][16][17]. In volume-averaged models the transport of solid equiaxed grains is described by grain population balances along with mass and solute mass conservation equations for the solid phase. These equations consist of contributions from advection at the macroscopic (process) scale and grain growth governed by diffusion and phase change phenomena at the microscopic (grain) scale. These are strongly coupled with the transport of heat, mass, chemical species and momentum in the liquid and solid phases. Identifying the complexity in solving the set of coupled equations, Založnik and Combeau [8] proposed an operator splitting scheme as a flexible method for integration of the macroscopic transport terms and the local growth terms. This method was successfully implemented to simulate large ingot casting [18]. Modeling and simulation of steel ingots is particularly demanding due to the size of the castings and the complexity of the multiphase flow, however, models and applications have been improving in recent years [19][20][21][22][23]. These models are based on the finite volume method (FVM), while the finite element method (FEM) has not yet been considered for volume-averaged multiscale modelling of solidification with transport of equiaxed grains. An implementation using FEM may be attractive for multiple reasons. First, FEM generally offers more flexibility and versatility than FVM in describing the boundaries of the domain to be analyzed, and defining the boundary conditions which prevail there. This is particularly true when considering structured FVM where the "staircase" effect along the boundaries is detrimental. FEM also offers opportunities to more simply couple with stress/strain structural analyses, eg. to model the occurrence of thermomechanical defects in solidified regions, as such analyses are generally also conducted using FEM [24,25].

In this paper, a numerical FEM solidification model is presented, accounting for microscopic phenomena as well as for the motion of solid and liquid phases. First, the macroscopic conservation equations and the constitutive relations describing interfacial interactions and exchanges are summarized. The resolution method and the numerical implementation for the set of non-linear equations is then detailed. This implies introducing an artificial diffusion which deals with the discontinuities at the packing front by FEM. Numerical simulation using the current model is then performed to validate the numerical implementation. Results demonstrate an efficient FEM resolution scheme implemented for the purely convective transport problem, which is difficult to solve numerically by FEM in the absence of diffusive effects.

Two-phase model of solidification

This section summarizes the governing equations taken into consideration in this case when modeling solidification in the presence of solid transport. The full nomenclature and a complete set of notations can be found in Appendix A. Eqs. ( 1)-( 4) represent the conservation of total mass, momentum, solute mass and energy, respectively, for each phase 𝛼 (𝛼 = 𝑠 for the solid or 𝛼 = 𝑙 for the liquid phase) [8]. These macroscopic conservation equations are obtained by averaging the microscopic equations over a representative elementary volume (REV). A detailed derivation of the governing equations can be found elsewhere [6,26].

Mass

𝜕 𝜕𝑡 (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 ) + ∇ • (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 ) = 𝛤 𝛼 + 𝛷 𝛼

Momentum 𝜕 𝜕𝑡 (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 ) + ∇ • (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 ⊗ 〈𝐯 𝛼 〉 𝛼 )

= -𝑔 𝛼 ∇〈𝑝 𝛼 〉 𝛼 + ∇ • (𝑔 𝛼 〈𝛕 𝛼 〉 𝛼 ) + 𝑔 𝛼 〈𝐛 𝛼 〉 𝛼 + 𝐌 𝛤,𝛼 + 𝐌 𝑑,𝛼 + 𝐌 𝛷,𝛼 Species 𝜕 𝜕𝑡 (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝑤 𝛼 〉 𝛼 ) + ∇ • (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝑤 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 )

= -∇ • (𝑔 𝛼 〈𝐣 𝛼 〉 𝛼 ) + 𝐽 𝛤,𝛼 + 𝐽 𝑗,𝛼 + 𝐽 𝛷,𝛼 Energy 𝜕 𝜕𝑡 (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈ℎ 𝛼 〉 𝛼 ) + ∇ • (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈ℎ 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 )

= -∇ • (𝑔 𝛼 〈𝐪 𝛼 〉 𝛼 ) + 𝑄 𝛤,𝛼 + 𝑄 𝑞,𝛼 + 𝑄 𝛷,𝛼

In the above equations ( 1)-( 4), the notation 〈 𝛼 〉 𝛼 indicates the intrinsic volume average in phase 𝛼, 𝑔 denotes the volume fraction, 𝜌 the density, 𝐯 the velocity, 𝑝 the pressure, 𝛕 the deviatoric part of the stress tensor, 𝐛 the body force per unit volume, 𝑤 the solute mass concentration, 𝐣 the solute flux vector, ℎ the specific enthalpy, 𝐪 the heat flux vector. The right-hand side of the above equations gathers the exchange terms rising from different microscopic processes: 𝛤 denotes the mass exchange rate due to phase change, 𝛷 the mass exchange rate due to grain nucleation, 𝐌 the vector for interfacial momentum exchange, 𝐽 the solute exchange rate, and 𝑄 the heat exchange rate. The contributions of nucleation (terms with superscript 𝛷) can be neglected compared to other terms in theses equations. The microscopic exchange contributions are modeled as follows, more detail can be found in

Appendix B.

due to phase change due to interfacial stresses or diffusion 𝐌 𝛤,𝛼 = 𝐯 ̅ 𝛼 * 𝛤 𝛼 𝐌 𝑑,𝛼 = 𝜌 𝛼 𝑅 𝛼 𝑆 v (𝐯 ̅ 𝛼 * -〈𝐯 𝛼 〉 𝛼 ) 𝐽 𝛤,𝛼 = 𝑤 ̅ 𝛼 * 𝛤 𝛼 𝐽 𝑗,𝛼 = 𝜌 𝛼 𝐷 𝛼 𝛿 𝛼 𝑆 v (𝑤 ̅ 𝛼 * -〈𝑤 𝛼 〉 𝛼 ) 𝑄 𝛤,𝛼 = ℎ ̅ 𝛼 * 𝛤 𝛼 𝑄 𝑞,𝛼 = ҡ 𝛼 𝑙 𝛼 𝑆 v (𝑇 ̅ 𝛼 * -〈𝑇 𝛼 〉 𝛼 )

Here 𝐯 ̅ * , 𝑤 ̅ * , ℎ ̅ * , 𝑇 ̅ * are the average values over the interface, R is the momentum resistance coefficient, D is the solute diffusion coefficient,  is the characteristic solute diffusion length, ҡ is the heat conductivity, l is the characteristic heat conduction length, 𝑆 v = 𝐴 𝑉 𝑒 ⁄ is the interfacial area concentration in the REV considered, of volume 𝑉 𝑒 .

The model is closed by balances of mass, momentum, solute, and heat at the solid-liquid interface:

∑ 𝛤 𝛼 + 𝛷 𝛼 𝛼=𝑠,𝑙 = 0 ∑ 𝐌 𝛤,𝛼 + 𝐌 𝑑,𝛼 + 𝐌 𝛷,𝛼 𝛼=𝑠,𝑙 = 0

∑ 𝐽 𝛤,𝛼 + 𝐽 𝑗,𝛼 + 𝐽 𝛷,𝛼 𝛼=𝑠,𝑙 = 0 ∑ 𝑄 𝛤,𝛼 + 𝑄 𝑞,𝛼 + 𝑄 𝛷,𝛼 𝛼=𝑠,𝑙 = 0

The population of grains is described by an average density per unit volume, N. The population balance equation writes:

𝜕𝑁 𝜕𝑡 + ∇ • (𝑁〈𝐯 𝑠 〉 𝑠 ) = 𝑁 ̇ (7)
where the nucleation rate is defined as:

𝑁 ̇= {
𝑁 0 is the nucleation density, 𝛿𝑡 the time step, and 𝑇 𝑛𝑢𝑐𝑙 the nucleation temperature. In Eq. (8) nucleation occurs under two conditions: The first nucleation event occurs when the local temperature drops below the nucleation temperature for the first time. Further nucleation events follow if the local grain density drops to zero (due to grain transport or remelting) and the local temperature is below the nucleation temperature. The latter criterion is a heuristic treatment used to prevent unphysical results [11].

In order to resolve the coupled microscopic and macroscopic phenomena, the microscopic exchange terms can be considered as source terms in the governing equations. The equations can then be integrated directly in a coupled way. However, due to the stiffness of the microscopic terms, such a solution scheme requires very small time steps. This makes the computation time realistically too long to simulate industrial sized castings. In the present work, an operator splitting method is applied to solve the entire system of equations [8].

According to this method, the phase fractions and solute compositions are considered to evolve in two subsequent stages: the macroscopic transport stage and the microscopic growth stage. When using this assumption, each stage can be numerically integrated by operating on different scales of time and space. Therefore, such a splitting technique provides an efficient way to solve multi-scale problems. This splitting method is used for the solution of the mass conservation equation for the solid, the conservation equation for the grain population density and the solute conservation equations for both phases. The method is summarized as follows:

first, in the macroscopic transport stage, only variation due to the macroscopic transport is integrated, which is determined by solving Eqs.( 9)- (11) on the global finite element mesh using a macro time step. This gives an intermediate quantity with index tr:

𝜕𝑔 𝑡𝑟 𝑠 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝑠 〈𝐯 𝑠 〉 𝑠 ) = 0 (9)
𝜕𝑁 𝑡𝑟 𝜕𝑡 + ∇ • (𝑁 𝑡𝑟 〈𝐯 𝑠 〉 𝑠 ) = 0 (10) 𝜕(𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 ) 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 〈𝐯 𝛼 〉 𝛼 ) = 0 (11) Second, in the microscopic growth stage, the contribution of microscopic processes, nucleation and growth, are integrated through Eqs.( 12)- (14). These equations are solved, locally, at each node of the finite element mesh, leading to quantity with index gr:

𝜕𝑔 𝑔𝑟 𝑠 𝜕𝑡 = 𝛤 𝑠 + 𝛷 𝑠 𝜌 𝑠 (12)
𝜕𝑁 𝑔𝑟 𝜕𝑡 = 𝑁 ̇ (13)
𝜕(𝑔 𝑔𝑟 𝛼 〈𝑤 𝛼 〉 𝑔𝑟 𝛼 ) 𝜕𝑡 = 𝐽 𝛤,𝛼 + 𝐽 𝑗,𝛼 + 𝐽 𝛷,𝛼 𝜌 𝛼

A smaller time step (micro time step) must be used for the microscopic growth stage. The solution of the macroscopic transport stage is the initial condition for the integration of the microscopic growth stage. The sequence of both integration steps thus gives the solution over a macroscopic time step. The whole modeling algorithm, using this splitting approach, is shown schematically in Figure 1.

Each iteration 𝜈 begins with an implicit finite element resolution for energy and transport conservation equations in the transport stage by using a macro time step. This gives the solutions denoted with superscript 𝜈 + 1 and subscript 𝑡𝑟 when they are associated with the transport stage. Then the nucleation-and-growth stage is solved locally (i.e., at each node), with variables initialized by values obtained from the transport stage at 𝜈 + 1. This local resolution proceeds through micro time steps, assuming that the average quantities for the solid plus liquid mixture, composition 〈𝑤〉 = 𝑔 𝑙 〈𝑤 𝑙 〉 𝑙 + 𝑔 𝑠 〈𝑤 𝑠 〉 𝑠 and enthalpy 〈ℎ〉 = 𝑔 𝑙 〈ℎ 𝑙 〉 𝑙 + 𝑔 𝑠 〈ℎ 𝑠 〉 𝑠 , no longer evolve during the macro time step. Finally, the momentum equations are solved with a semi-implicit solver on the macro time step to compute the new estimation of velocity fields of the liquid and the solid phases at iteration 𝜈 + 1. A complete evolution of the different quantities over the time step is evaluated by the final results obtained from the growth stage, as these solutions already include the change from the transport stage. The splitting scheme is only used to solve the evolution of phase fractions, grain density and solute concentrations since it involves very different scales of time and space. The resolution of the energy and momentum conservation equations do not require operator splitting because the constitutive relations coupling the micro-and macroscopic scales are simpler.

Regarding the transport stage, Eqs. ( 9)-( 11) are of pure convective nature, and notoriously difficult to solve numerically in the absence of diffusive effects. Moreover, another numerical difficulty arises from the discontinuity of transport velocities due to the solid packing phenomenon. Indeed, when forming a packed solid layer, solid grains suddenly change from a moving state to a fixed state. Solving these issues in the framework of FEM for solidification

simulations remains an open issue and will be addressed in the following sections.

Treatment of discontinuities at the packing front by the FEM

The solution of convection-dominated equations by FEM can encounter problems of unphysical oscillations, especially in zones with steep gradients. Different resolution techniques have been developed to overcome such issues: e.g. discretizing by upwind methods [27][28][29], by stabilized methods such as Streamline Upwind Petrov-Galerkin [30],

Galerkin Least-Squares [31,32], Residual Free Bubbles [33], or by using a corrected flux approach [34]. A specific complexity in the resolution of Eqs. ( 9)-( 11) arises from discontinuities due to the solid packing phenomenon taking place at the interface between the liquid pool and the packed bed of solid grains. On one side, small grains move freely in the liquid phase, the solid fraction being smaller than a characteristic packing fraction. On the other side, grains accumulate and grow to form a steady and fixed packed zone, the solid fraction being higher than the packing limit. Discontinuities are then related to abrupt changes in the velocity and fraction of the solid phase. Therefore, an adaptive artificial diffusion is introduced to stabilize the finite element resolutions without unreasonably smearing results.

This added diffusion detects and reduces discontinuities at locations where the solid phase is being packed. The added diffusion is then not present everywhere, it is restricted to critical zones of packing. The diffusion coefficient, 𝐷 𝑀 , as expressed below, consists of the gradient of solid velocities, i.e. a combination of the divergence of the average solid velocity ∇ • 〈𝐯 𝑠 〉 and the divergence of the intrinsic average solid velocity ∇ • (〈𝐯 𝑠 〉 𝑠 ), which contain information about the variation of the related quantities: velocity and fraction of solid phase.

𝐷 𝑀 = (ℎ 𝑒 𝐯 𝑠 ) 2 (𝛼 |∇ • 〈𝐯 𝑠 〉| + 𝛽|∇ • 〈𝐯 𝑠 〉 𝑠 |) (15)
In this expression, ℎ 𝑒 𝐯 𝑠 denotes the characteristic size of element 𝑒 in the direction of the solid velocity 𝐯 𝑠 , as proposed in [24]: (16) where 𝜑 𝑖 is the interpolation function associated with node 𝑖, 𝐯 𝑐𝑒𝑛𝑡𝑒𝑟 𝑠 is the solid velocity at the center of the element, and 𝑁𝑛 is the number of nodes per element. The coefficients 𝛼 and 𝛽 allow control over the amount of diffusion in a direct way, thus offering more flexibility than an implicit diffusion introduced by the upwind method. An alternative option would be to use the divergence of the intrinsic solid velocity, ∇ • 〈𝐯 𝑠 〉 𝑠 , and the solid fraction gradient, ∇𝑔 𝑠 . However, the former expression is preferred because it does not exist without solid transport. Hence the model is still valid in cases without solid motion. The formulation of the coefficient 𝐷 𝑀 is inspired by the work of Cook and Cabot who developed an artificial nonlinear diffusion using the entropy gradient to treat issues associated with discontinuities of temperature and mass fraction in supersonic reacting flows [35]. For the aforementioned issues experienced when simulating solidification there has, as of yet, been no relevant investigations.

ℎ 𝑒 𝐯 𝑠 = 2‖𝐯 𝑐𝑒𝑛𝑡𝑒𝑟 𝑠 ‖ ∑ |𝐯 𝑐𝑒𝑛𝑡𝑒𝑟 𝑠 • ∇𝜑 𝑖 | 𝑁𝑛 𝑖
It should be noted that it is necessary to use the same artificial diffusion coefficient for all relevant transport equations in order to ensure consistency between the transport of related quantities and the conservation of mass and energy. Consequently, the transport equations for solid phase, grain density and solute, derived from Eqs. ( 9)- (11), are modeled by Eqs. ( 17)- (19). These equations are solved in the finite element framework using an implicit scheme in time and a weighted residual approach with P1 linear elements (triangles in 2D, tetrahedra in 3D) and SUPG stabilization method.

𝜕𝑔 𝑡𝑟 𝑠 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝑠 〈𝐯 𝑠 〉 𝑠 ) -∇ • (𝐷 𝑀 𝛻𝑔 𝑡𝑟 𝑠 ) = 0 ( 17
)
𝜕𝑁 𝑡𝑟 𝜕𝑡 + ∇ • (𝑁 𝑡𝑟 〈𝐯 𝑠 〉 𝑠 ) -∇ • (𝐷 𝑀 𝛻𝑁 𝑡𝑟 ) = 0 (18) 𝜕(𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 ) 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 〈𝐯 𝛼 〉 𝛼 ) -∇ • (𝐷 𝑀 𝛻(𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 )) = 0 (19)

Formulation of Energy Conservation

Regarding heat transfer, the equation for energy conservation, Eq.( 4), for the solid (𝛼 = 𝑠)

and for the liquid (𝛼 = 𝑙) phases, assuming thermal equilibrium between both phases in the REV, and introducing the added diffusion, the following mixture energy equation can be established:

𝜌 [ 𝜕〈ℎ〉 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝑠 〈ℎ 𝑠 〉 𝑠 〈𝐯 𝑠 〉 𝑠 + 𝑔 𝑡𝑟 𝑙 〈ℎ 𝑙 〉 𝑙 〈𝐯 𝑙 〉 𝑙 ) -∇ • (𝐷 𝑀 ∇(𝑔 𝑡𝑟 𝑠 〈ℎ 𝑠 〉 𝑠 + 𝑔 𝑡𝑟 𝑙 〈ℎ 𝑙 〉 𝑙 ))] -∇ • (〈ҡ〉∇𝑇) = 0 (20)
where the average enthalpy and thermal conductivity are defined by 〈ℎ〉 = 𝑔 𝑠 〈ℎ 𝑠 〉 𝑠 + 𝑔 𝑙 〈ℎ 𝑙 〉 𝑙 and 〈ҡ〉 = 𝑔 𝑠 〈ҡ 𝑠 〉 𝑠 + 𝑔 𝑙 〈ҡ 𝑙 〉 𝑙

Additional it is assumed that the densities of phases are constant and equal and that the heat diffusion follows the Fourier law. Using this method, the phase fractions in the advection terms should be taken as those calculated in the transport stage, so that mass conservation is maintained [8]. In the present work, the enthalpy formulation of the energy equation is used.

Like the preceding transport equations, the energy equation is solved in the framework of the finite element formulation with an implicit scheme for time integration as well as a weighted residual approach with P1 linear elements and SUPG stabilization method.

Formulation of Momentum Conservation

Several assumptions are made when solving the momentum equations, including:

i.

The phase densities are constant and equal, except for the buoyancy forces for which the liquid and solid densities are modeled as follows, respectively.

𝜌 𝐵 𝑙 = 𝜌 𝑟𝑒𝑓 [1 -𝛽 𝑇 (𝑇 -𝑇 𝑟𝑒𝑓 ) -𝛽 𝑤 (〈𝑤 𝑙 〉 𝑙 -𝑤 𝑟𝑒𝑓 )] (22)
𝜌 𝐵 𝑠 = 𝜌 𝑟𝑒𝑓 1 -𝛽 𝑠ℎ𝑟 (23) ii. The pressure is assumed to be identical in the solid and liquid phases.

〈𝑝 𝑙 〉 𝑙 = 〈𝑝 𝑠 〉 𝑠 = 𝑝 (24)
iii. The liquid behaves as a Newtonian fluid with a constant viscosity. Neglecting the interfacial momentum transfer due to phase change, the divergence of the average deviatoric stress tensor is modeled as:

∇ • 〈𝛕 𝑙 〉 = 𝜇 𝑙 ∇ • (∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) 𝑇 ) (25)
iv.

The momentum transfer due to nucleation and growth is considered negligible relative to other terms, the momentum balance at the solid-liquid interface is thus described as:

𝐌 𝑑,𝑙 + 𝐌 𝑑,𝑠 = 𝟎 (26)
Using these assumptions, the liquid momentum equation is derived as Eq. ( 27) from the general formulation for phase 𝛼, Eq. (2).

𝜌 [ 𝜕 𝜕𝑡 (𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇ • (𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ⊗ 〈𝐯 𝑙 〉 𝑙 )] = -𝑔 𝑙 ∇𝑝 + 𝜇 𝑙 ∇ • (∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) 𝑇 ) + 𝑔 𝑙 𝜌 𝐵 𝑙 𝐠 -𝐌 𝑟𝑒𝑔𝑖𝑚𝑒 𝑑,𝑠 (27)
where the momentum transfer due to the drag force is described as

𝐌 𝑟𝑒𝑔𝑖𝑚𝑒 𝑑,𝑠 = ℳ 𝑟𝑒𝑔𝑖𝑚𝑒 (𝑔 𝑙 ) 2 (〈𝐯 𝑙 〉 𝑙 -〈𝐯 𝑠 〉 𝑠 ) (28)
in which 3 in the slurry regime (𝑔 𝑠 < 𝑔 𝑐 𝑠 )

ℳ 𝑟𝑒𝑔𝑖𝑚𝑒 = { ℳ 𝑠𝑙𝑢𝑟𝑟𝑦 = 3𝑔 𝑠 𝜇 𝑙 𝐶 𝑑 𝑅𝑒 4(𝑑 𝑔 ) 2 (𝑔 𝑙 )
ℳ 𝑝𝑎𝑐𝑘𝑒𝑑 = 180(1 -𝑔 𝑙 ) 2 𝜇 𝑙 𝜆 2 2 (𝑔 𝑙 ) 3 in the packed -bed regime (𝑔 𝑠 ≥ 𝑔 𝑐 𝑠 ) (29)
where the solid-liquid interaction follows the model of Agarwal and O'Neill [36].

In the finite element context, solving the above equation, Eq. ( 27), and distinguishing the two regimes on a single discretized domain may cause numerical difficulties because of the assembly of dissimilar and irregular terms (i.e. slurry regime contributions and packed-bed regime contributions). In order to overcome this issue, a smoothing procedure was developed, similar to the one introduced by Plotkowski and Krane [15]. It consists in using a transition function, 𝛼 𝑡 , to switch the momentum transfer, 𝐌 𝑟𝑒𝑔𝑖𝑚𝑒 𝑑,𝑠

, between the two regimes. The following expression is proposed for this transition function, which provides a compromise between having a sufficiently sharp change of regimes while avoiding an abrupt switch between them:

𝛼 𝑡 = 1 -0.5(1 -𝑡𝑎𝑛ℎ[𝛼 0 (𝑔 𝑐 𝑠 -𝑔 𝑠 )]) (30)
where the value 𝛼 0 = 100 is chosen, as plotted in Figure 2. The liquid momentum equation can then be expressed as

𝜌 [ 𝜕 𝜕𝑡 (𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇ • (𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ⊗ 〈𝐯 𝑙 〉 𝑙 )] (31)
= -𝑔 𝑙 ∇𝑝 + 𝜇 𝑙 ∇ • (∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) 𝑇 ) + 𝑔 𝑙 𝜌 𝐵 𝑙 𝐠 -[𝛼 𝑡 ℳ 𝑠𝑙𝑢𝑟𝑟𝑦 + (1 -𝛼 𝑡 )ℳ 𝑝𝑎𝑐𝑘𝑒𝑑 ](𝑔 𝑙 ) 2 (〈𝐯 𝑙 〉 𝑙 -〈𝐯 𝑠 〉 𝑠 )
Which is solved by a semi-implicit time integration and a weighted residual method, precisely by the P1/P1 velocity-pressure formulation stabilized with the Variational Multi Scale method [37]. In the equation for the conservation of solid phase momentum, the inertial and viscous terms are neglected [8]. Therefore, the solid momentum equation, which is only considered in the slurry regime, reduces to:

-𝑔 𝑠 ∇𝑝 + 𝑔 𝑠 𝜌 𝐵 𝑠 𝐠 + ℳ 𝑠𝑙𝑢𝑟𝑟𝑦 (𝑔 𝑙 ) 2 (〈𝐯 𝑙 〉 𝑙 -〈𝐯 𝑠 〉 𝑠 ) = 0 (32) From this equation, the solid velocity can be calculated locally, at each node of the finite element mesh. In order to avoid a sudden change of the velocity when solid packing occurs, a gradual transition is applied, by using the same switching function 𝛼 𝑡 introduced in Eq. (30).

Therefore the solid velocity is calculated at each node by:

〈𝐯 𝑠 〉 𝑠 = 𝛼 𝑡 (〈𝐯 𝑙 〉 𝑙 + 𝑔 𝑠 (𝜌 𝐵 𝑠 𝐠 -∇𝑝) ℳ 𝑠𝑙𝑢𝑟𝑟𝑦 (𝑔 𝑙 ) 2 ) ( 33
)
Figure 3. Adjustment strategy for the velocity of convected grains in the vicinity of the packed bed.

Furthermore, it is necessary to ensure that moving solid grains do not penetrate the preexisting packed bed, to avoid the solid fraction at the packing front exceeding the packing fraction. The velocity of the grains is therefore adjusted so that they land smoothly on the packed bed, i.e. reaching 〈𝐯 𝑠 〉 𝑠 = 0 at the packing limit. The algorithm for this adjustment is presented in Figure 3. It consists first of the calculation of the distance that the grains would travel in the direction of the solid fraction gradient, moving with their current velocity during the time step. This settling distance is 〈𝑣 𝑠 〉 𝑝𝑟𝑜𝑗 𝑠 ∆𝑡, where 〈𝑣 𝑠 〉 𝑝𝑟𝑜𝑗 𝑠 is the projected solid velocity and ∆𝑡 the time step. Then the settling distance is compared to the distance between the grains and the packing limit, 𝑑 𝑃-𝐿𝑖𝑚𝑖𝑡 , defined by the packing solid fraction, 𝑔 𝑐 𝑠 . If the settling distance is higher, the velocity is corrected by the factor 𝛼 as shown in Figure 3. This correction prevents the grains from traveling beyond the packing limit.

Verification of the adaptive artificial diffusion method for particle packing

The model described above is sophisticated as it involves artificial diffusive terms in conservation equations, a smoothing procedure for the transition from a slurry regime to a packed bed regime in the momentum conservation, as well as an adjustment strategy for the convection of solid in the vicinity of the packed bed. Several studies are then needed to evaluate the numerical parameters of the model. The first of which is presented below for a 1D sedimentation configuration where a simple analytical solution can be derived.

Test Case Description

In order to test these new parameters, the model is first applied to simulate a one-dimensional (1D) sedimentation process, considering only transport and neglecting nucleation and growth processes. It consists in the settling of a predefined number of globular grains with equal and constant size at uniform and constant velocity. A schematic of the test is presented in Figure 4 and the 1D analytical solution is derived in Appendix D. The computational domain is two-dimensional (2D), with a width and length of 1 mm × 100 mm. There is no heat exchange through the boundaries of the domain. The velocities at the top and bottom faces as well as the normal velocities along the vertical walls are set to zero and a perfect slip condition is applied to the tangential velocities on the vertical walls, i.e. returning to the 1D configuration described in Appendix D. Initial conditions are given in Figure 4. The binary alloy Sn-5wt%Pb is considered, its thermophysical properties can be found in Appendix C [8]. The enthalpy is evaluated according to the solid fraction and the temperature. The simulation parameters are given in Table 1.

The present test case is defined to benefit from the simple analytical solution shown in Figure 5 and Figure 6 at several times with dashed lines. The downward velocity of the settling grains is arbitrarily imposed to be a constant value in the unpacked region, equal to 1 mm s -1 and directed toward the -y-axis. The corresponding upward liquid velocity in the unpacked region was computed and found to be constant, equal to 0.111 mm s -1 , in agreement with the total mass balance. When the fraction of solid reaches 0.3 in the packed bed the velocity of the phases falls to zero. As there is no solidification, the fraction of solid cannot increase further. Conservation of the initial mass of solid thus defines the height of the packed bed. It reaches 20 mm once settling is complete. Also considering the adiabatic boundary conditions for heat transfer as well as the absence of phase change by solidification/remelting, the temperature is expected to remain constant and uniform throughout the simulation domain, equal to its initial value 498 K (224°C) shown in Figure 5 and

Analysis of Simulation Results

The initial solid velocity is set to 1 mm s -1 for the vertical component in the downward direction. The liquid velocity is then determined by solving the momentum equation, Eq. (31).

The results are presented in Figure 5 and Figure 6. Those obtained from the numerical simulation are shown as solid lines while those calculated by the analytical solution are dashed lines. The sedimentation process is illustrated in Figure 5-a as profiles of solid volume fraction at t = 0 s, 10 s (before grains reach the bottom of the domain), 30 s (accumulation below the packing limit), and 60 s (end of packing). From the initial state, the solid grains fall downward while the liquid moves upward in such a way that the continuity equation is respected, shown in Figure 5-b and -c. Since the solid fraction in the two-phase region is 0.1, the liquid velocity is ten times smaller than the solid velocity. The first grains reach the bottom after 20 s and begin to accumulate until the packing fraction (chosen as 0.3) is reached. During the process, the change of the solid phase from the moving to the packed state has an impact on the liquid movement, which is revealed by the peaks of liquid velocity at the packing front. The sedimentation is complete after 60 s, although the solid and liquid velocities remain non-zero within a small layer where the gravitational force is balanced with the diffusive effect. This phenomenon is maintained due to the persistent gradient of the solid fraction at the transition interface between the solid packed bed and the solid-free region.

The analytical and numerical results show similar trends, however, there exist differences between the two solutions. The differences found within the transition zones are due to diffusive effects. The simulation results contain an inevitable numerical diffusion and in the current case, at the boundary of the packed bed, an additional artificial diffusion. Furthermore, it can be observed that the solid fraction in the packed zone exceeds the predefined packing value of 0.3. This over accumulation is related to purely numerical issues when packing the solid phase and it will be discussed further in the next section.

As grain motion also involves the transport of solute and heat, it is important to verify the consistency of all transported quantities. The profiles of the average composition in about 30% solid at 0.364 wt% Pb and 70% liquid at 5 wt% Pb. Furthermore, as expected for a pure transport phenomenon, the temperature does not change during this process. Only a slight deviation from the initial temperature, smaller than 1 K, can be seen at 60 s, as shown in Figure 5-e. A good overall conservation of all quantities is verified, the maximal relative errors for the global solute mass and energy being about 10 -6 and 10 -5 , respectively.

Effects of Artificial Diffusion

Two simulations were performed in order to study the impact of the artificial diffusion term, by separately assessing the effect of the constant parameters introduced with variable 𝐷 𝑀 . In the first case, the coefficients are 𝛼 = 1 and 𝛽 = 0 while in the second one, 𝛼 = 0 and 𝛽 = 1. These problems are overcome in Case 2, the results of which are shown in Figure 6-d, -e and -f. This reveals that the term related to the variation of the intrinsic velocity has a more important role in solving numerical singularities than the average velocity. An over-packing relative to the predefined fraction is also stated with a higher excessive quantity than when 𝛼 = 20; 𝛽 = 1. This observation further shows the role of artificial diffusion in treating numerical issues related to the packing of solid, since the higher the added diffusion the better the simulation respects the predefined packing critical solid fraction. However, using high values of 𝛼 and 𝛽 lead to an excess of diffusion and produce unphysical solutions.

This set of simulations confirms that numerical concerns of finite element resolution in the presence of sharp discontinuities of the transport velocity field in an absence of diffusion in hyperbolic equations can be solved by adding a supplementary diffusive component.

Modeling Applications -Macrosegregation Simulations

The following section presents test cases, first a two-dimensional (2D) and then three-dimensional (3D) applications for the solidification of a small Sn-Pb ingot. The objective is to further test the model while maintaining small geometries. Finally, the model being validated, a simulation for a 3D steel ingot on the scale of an industrial process will be given in order to discuss the application to a real casting geometry and to identify future possible improvements of the model.

Two-dimensional Test Case

Description

The studied case, presented in Figure 7, is configured according to the Hebditch-Hunt benchmark study [38]. A 100 mm x 60 mm cavity contains the same Sn-5wt%Pb alloy as in the previous sections, initially in the liquid state at 499.15 K (226 °C). Cooling takes place from the left side, an environment at 25 K, with a heat transfer coefficient of 300 W m -2 K -1 .

The rest of the boundary is assumed adiabatic. The simulation is performed with the complete model, which accounts for nucleation and growth processes, and transport phenomena, shown schematically in Figure 1. Nucleation happens at sites where the liquid is cooled below the liquidus temperature and where there are no existing grains. A homogenous grain density of [8] (right column, FV): (a) solid fraction , 𝑔 𝑠 , and intrinsic solid velocity vectors, 〈𝐯 𝑠 〉 𝑠 , (b) temperature, 𝑇, and intrinsic liquid velocity vectors, 〈𝐯 𝑙 〉 𝑙 . Black curves are isolines of solid fraction (0.1; 0.2; 0.3). The simulation case is defined in Figure 7 and Table 2. stationary mushy zone has a negative segregation which results from accumulation of solutedepleted grains. However, the average composition is not uniform. At the very bottom of the cavity, there is an accumulation of solute which results from intergranular melt flow localized along the bottom wall and oriented in the direction of the temperature gradient. In the left part of the packed layer the negative segregation is amplified by the flow of intergranular liquid, which is oriented against the temperature gradient in this region. As the process advances, the channel continues to extend along the bottom side.

Results and discussion

When comparing FE and FV methods, the map produced at the end of solidification shows a similar tendency of segregation formation, including the negative segregation located in the left zone of the cavity and the positive channel formed at the bottom. However, the FE solution produces a less marked negative segregation, and a larger positive channel along the bottom wall. Additionally, in the upper zone of the cavity, the FE solution shows positive segregation near the upper-left corner and a slightly negative segregation nearby. In this region, the FV computation produces positively segregated channels with higher solute content, distributed horizontally and alternating with negatively segregated zones. Differences between these two results can be partly explained by the influence of numerical factors.

Because of the extremely high nonlinearity of the problem, any differences in the numerical solution methods can lead to noticeable differences between the solutions. The artificial diffusion used in the FE method inevitably leads to smoothing of the macroscopic fields and thus to a smaller degree of segregation than in the FV solution. In addition, many other factors, including discretization schemes, iteration procedures, etc., can be the cause of the differences between the FE and the FV results. The reader interested in such aspects can refer to literature in the context of simulations with a fixed solid phase [39,40].

When considering the distribution of solid grains it is useful to look at the distribution of grain density in Figure 10 at different instants. At 10 s, 20 s, and 200 s, it can be first observed that a large number of grains are gathered in the packed layer. It can also be seen that the transport of crystals by liquid advection induces several zones of significant grain density in the slurry region. As solid grains are transported, such a heterogeneous distribution is expected.

Although the results obtained from the FE and FV methods are still very similar at t = 10 s, larger differences between the two solutions are perceived later on. Compared to the FV results, a higher grain density in the right region is predicted by the FE simulation. It is known that the number of grains is governed and influenced by different coupled processes, consisting of nucleation, transport, re-melting and re-nucleation mechanisms. Although the physical parameters and the numerical procedure for nucleation are identical in both simulations, the larger quantity of grains in the FE solution might be caused by a higher frequency of nucleation events at nodes that were emptied of grains because of transport or remelting. Nevertheless, further investigations should be carried out, as the sources of those differences still remain unconfirmed. Despite these differences, clear similarities between both results can be observed, including a high grain density in the lower-right corner and a low grain density in the upper-left corner at the end of solidification. 3.

Table 3. Simulation parameters for the 3D test case presented in Figure 11. Results show very similar trends. However, it can be observed that the distributions of variables are non-uniform within the thickness of the 3D geometry. In the solid fraction maps (Figure 12-a), across the horizontal planes, it can be seen that the solid zone advances further along the specimen walls (plane P3) than in the interior. Thus the iso-surface of solid fraction exhibits a concave shape in the central zone of the cavity. This phenomenon can be explained by considering the distributions of other quantities since all relevant variables are closely related. First, as seen on the maps of velocities (Figure 12-e and f), there is no motion of neither solid nor liquid along the walls due to the no-slip boundary conditions. Because the grains remain attached to the cavity walls and do not settle, the solid fraction is higher than in the interior of the cavity (plane P6). The solid and liquid phases are not constrained in the immediate vicinity of the external surfaces. With a high quantity of mobile solid grains, the settling velocity is thus stronger there than that in the central zone as observed on the solid velocity map (Figure 12-f). This solid motion then enhances the downward movement of the solute-enriched liquid phase. Comparison of the maximum velocities at 10 s is possible, showing (2D, solid phase) 40.4 mm s -1 versus (3D, solid phase)

Results and discussion

41.0 mm s -1 and (2D, liquid phase) 39.6 mm s -1 versus (3D, solid phase) 31.5 mm s -1 . In both cases, the location of the highest velocities for the solid and liquid phases are very close. The solid velocity being directly computed from eq. ( 33) by neglecting the inertial and viscous terms, i.e. only accounting for the solid-liquid interaction through the transferred momentum due to interfacial stresses, maximum values are almost equal in the 2D and 3D simulations.

However, a lower value is found in the 3D simulation for the maximum velocity of the liquid phase. This is due to the interaction of the liquid flow with the two largest cavity walls (plane P3 and its symmetric), not accounted for in the 2D approximation. Thus, the transport of the liquid phase by the solid phase, while being obviously present, plays a less important role in comparison to the 2D approximation.

It is also interesting to observe the distribution of the flow in the horizontal cross sections.

The maximum velocities for both phases are not observed at the symmetry plane of the cavity.

Instead they are localized at about 1/3 of the half-cavity thickness from the cavity surfaces.

This uneven distribution is enhanced when successively considering the cross sections from top (height 57 mm) to bottom (height 17 mm). 3.

Three-dimensional Simulation of an Industrial Steel Ingot

Description

The numerical model is now applied to predict macrosegregation in an industrially cast steel ingot as produced by the company Aubert & Duval. Physical features and alloying components of the casting are detailed in [18], while its schematic is illustrated in Figure 14a. The bottom and lower region of the steel ingot are cooled via a mold while the upper region is surrounded by an insulating refractory material. The top surface is covered by a layer of exothermal powder. The mold outer surface thermally exchanges with the environment by convection and radiation.

For the current study, the steel is considered as a binary alloy composed of an iron base and a 0.36wt% nominal carbon content, which plays a dominant role in determining buoyancy force compared to other chemical elements [18]. A simulation is performed on one quarter of a cylinder (0.3 m radius and 1.8 m height), an approximation of the octagonal cross-section of the real ingot, which is bounded by two symmetric planes (P2 and P3 in Figure 14-b). The The simulation is carried out with a uniform mesh size of 20 mm and a constant time step of 0.01 s. 4.

Results and discussion

Table 4. Simulation parameters for the ingot case presented in Figure 14.

As it can be observed in these figures, the flow descends along the cooled side and ascends along the centerline, resulting in a global circulation loop in the shape of an elongated torus.

The maximum velocity is around 150 mm/s along the cooled wall. The persistence of this circulation loop is remarkable and its flow direction is the reverse of the direction observed when solid transport is not taken into account. If the solid is assumed to be fixed to the mold, the flow is driven only by natural convection induced by density differences in the liquid. The density of the liquid depends on the temperature and on the chemical composition; the concentration of carbon is most important when compared to all solutes in the multicomponent steel [18]. In the mushy zone the liquid density decreases with decreasing temperature because the liquid is enriched in carbon as solidification progresses. This creates a lateral liquid density gradient from the cooled walls towards the core of the casting, which drives the thermosolutal natural convection. If the solid is fixed, lighter liquid ascends along the walls and heavier liquid descends in the core. When solid grains move, the flow is completely modified. Solid crystals, heavier than the liquid, sediment along the cooled walls and entrain the liquid, thus inducing a downward flow. This phenomenon leads to an overall flow loop descending at the solidification front and ascending in the ingot core. This clearly shows that besides the thermo-solutal effect the solid transport plays a significant role in the formation of the natural convective flow during solidification. Moreover, the motion of solid and liquid are strongly coupled; the moving phases can entrain one another via drag forces.

The predicted flow structure is clearly three-dimensional and is not axisymmetric, although a four-fold symmetry is implicitly assumed by the choice of the computational domain. (about 150 mms -1 ) along the cooled wall. Despite a higher solid fraction near the wall (which can be seen in the transverse cross-sections) there is no permanent solid layer attached to this cooled wall. Additionally, in the velocity maps in Figure 15 it can be seen that after sinking to the bottom along the cooled side, the mobile solid phase is transported towards the center zone by liquid flow. In this way, solid grains coming from the outer solidified region first accumulate at the center area and then extend to the side wall, resulting in a packed solid built up from the bottom side. Since solute-depleted grains settle and occupy the lower zone, the liquid enriched in solute is ejected upwards. This gives rise to the formation of a negative segregation cone in the lower zone of the ingot, as shown in Figure 17-a. This is a typical phenomenon experimentally found in steel ingots. cylindrical FVM considering a spherical globular solid grains and (green curve) 2D cylindrical FVM simulation accounting for a dendritic morphology of the solid phase [18].

The simulation case is defined in Figure 14 and Table 4.

Conclusions

In this study, a finite element solidification model which takes into account the transport of equiaxed grains is presented. This model consists of  the resolution of a set of highly nonlinear and strongly coupled equations over multiple scales in time and space, including those of energy, phase movement, phase transport, grain density transport, solute transport, nucleation and solid growth,  the coupling of the equations based on the operator splitting algorithm, previously developed by Založnik and Combeau [8], demonstrated as an effective way for the numerical resolution of the evolution of solidification structures in the growth stage and the transport stage.

Different issues make the finite element implementation challenging:

 solution of pure transport equations by the finite element method,  sharp discontinuities in the velocity fields due to the packing of solid grains.

Propositions were introduced to overcome these difficulties:

 addition of an adaptive artificial diffusion to the transport equations,  implementation of a specific treatment to deal with the packing issue, consisting in using a transition function and adjusting the solid velocity.

A careful investigation was conducted to ensure the consistency between related quantities during the process as well as to guarantee the conservation of mass and energy. It progressively consisted of:

 a 1D pure transport simulation of sedimentation to evaluate its effects and to propose appropriate values for the adaptive artificial diffusion,  application of the complete transport-growth model to simulate macrosegregation in a 2D configuration [11],

 extension of the above simulation in 3D, resulting in a heterogeneous distribution of variables in the third direction which could not be captured by 2D simulations,  3D simulation of the solidification of a 3.3 ton Fe -0.36 wt% C steel ingot, representative of a real ingot [18], showing macrosegregation prediction in reasonable agreement with experimental measurements.

To our knowledge, the present development is original in the context of the finite element method. It should be noted that the computational time reached 35 days for the simulation of the binary Fe-0.36wt%C alloy. Improvements are expected when using a combination of adaptive techniques for the macroscopic time step and the FE mesh. The present finite element model could then become a promising tool to simulate solidification, especially for industrial applications such as ingots of complex geometries and large size. It also has potential for coupling with segregation due to thermomechanical deformation while accounting for the grain structure formed during casting.

𝐛

Appendix B -Solute Diffusion Lengths and Area Concentration

The solute diffusion lengths are taken from the work of Tveito and co-workers [41], as the following formulations.

Solute diffusion length in the liquid phase

𝛿 𝑙 = 𝑤 𝑙 * -〈𝑤 𝑙 〉 𝑙 - 𝜕𝑤 𝑙 𝜕𝐧 | * (C. 1) = 𝑑 𝑑 𝑅 𝑔 - 𝑓(𝑅 𝑔 , ∆) + 𝑔(𝑅 𝑓 , 𝑅 𝑔 , ∆) 𝑑[𝑅 𝑔 + 𝑑 -(𝑅 𝑔 + ∆ + 𝑑)𝑒 -∆ 𝑑 ⁄ ] -𝑓(𝑅 𝑔 , ∆) + 𝑔(𝑅 𝑓 , 𝑅 𝑔 , ∆)(𝑒 -∆ 𝑑 ⁄ -1) (C. 2)
where

𝑑 = 𝐷 𝑙 𝑣 𝑓(𝑅 𝑔 , ∆) = (𝑅 𝑔 + ∆) 2 -𝑅 𝑔 2 2 (C. 3) 𝑔(𝑅 𝑓 , 𝑅 𝑔 , ∆) = 𝑅 𝑓 3 -(𝑅 𝑔 + ∆) 3 3(𝑅 𝑔 + ∆) ∆= 𝑚𝑖𝑛 (𝑅 𝑓 -𝑅 𝑔 ; 2𝑅 𝑔 𝑆ℎ 𝑐𝑜𝑛𝑣 ) (C. 4)
𝑆ℎ 𝑐𝑜𝑛𝑣 = 2 3𝑔 𝑙 𝑆𝑐 1 3 ⁄ 𝑅𝑒 𝑛(𝑅𝑒) 𝑆𝑐 = 𝜇 𝑙 𝜌 𝑙 𝐷 𝑙 The area concentration is calculated as:

𝑆 v = 4𝜋(𝑅 𝑔 ) 2 𝑁 (C. 8)

Appendix D -Analytical solution for the 1D Test Case

The 1D Test Case consists in pure sedimentation of a column of preexisting globular grains with fixed size in a uniform temperature domain. Considering constant and equal densities of the solid and liquid phases, as well as no phase change and no nucleation, the average total mass conservation simplifies to 𝑔 𝑠 〈𝒗 𝑠 〉 𝑠 + 𝑔 𝑙 〈𝒗 𝑙 〉 𝑙 = 0. For the sake of simplicity, a constant settling value of the solid velocity is imposed, set to 〈𝑣 𝑠 〉 0 𝒔 = -1 𝑚𝑚 𝑠 -1 . The 1D domain height and the initial conditions are defined in Figure 4: a continuous and uniform 60 mm mushy zone region is initially present between heights 20 mm and 80 mm, with a uniform average grain density per unit volume, 𝑁 0 = 10 9 𝑔𝑟𝑎𝑖𝑛𝑠 𝑚 3 , and volume fraction of solid, 𝑔 0 𝒔 = 0.1. One can easily derive the value for the liquid velocity in the mushy zone,

〈𝑣 𝑙 〉 𝑙 = -𝑔 0 𝒔 〈𝑣 𝑠 〉 0 𝒔 /(1 -𝑔 0 𝒔 ) = 0.11 𝑚𝑚 𝑠 -1 .
Similarly, the radius of the grains, 𝑅 0 , is simply given by using the definition of the fraction of solid, 𝑔 0 𝒔 = 𝑁 0 (4/3)𝜋𝑅 0 3 , leading to the value 𝑅 0 = 0.288 𝑚𝑚. Considering the fixed settling velocity and the packing limit at which the grain stop, 𝑔 𝑐 𝑠 = 0.3, the time evolution of the distribution of the mushy zone is simply derived by considering that the total fraction of the solid phase is unchanged over the entire domain, while not exceeding 𝑔 𝑐 𝑠 in the packed bed. Values are reported in Table D1.

The temperature is fixed to 498 K (224.856 °C), i.e. below the liquidus temperature of the 𝑆𝑛 -5 𝑤𝑡% 𝑃𝑏 alloy, that is 498.72 K (225.57 °C) according to the thermophysical properties listed in Table C1 of Appendix C [8]. The average solute mass composition is defined by 〈𝑤〉 = 𝑔 𝑠 〈𝑤 𝑠 〉 𝑠 + 𝑔 𝑙 〈𝑤 𝑙 〉 𝑙 . At any time, as the system is closed with respect to mass transfer, integration over the entire domain must retrieve the nominal composition of the alloy, 𝑤 0 = 5 wt% Pb. The initial composition profile assumes no macrosegregation. This means that the average composition is equal to 𝑤 0 at any position along the domain.

However, assuming complete mixing in both liquid and solid phases, the lever rule holds and one can derive the equilibrium intrinsic composition of the liquid and solid phases, 〈𝑤 𝑙 〉 𝑙 = 5.556 wt% Pb and 〈𝑤 𝑠 〉 𝑠 = 0.364 wt% Pb, respectively. Knowing the distribution of solid and liquid and their initial and intrinsic compositions, one can directly compute the average compositions by tracking the change of phases due to sedimentation. Computed values are reported in Table D1.

Table D1. Time evolution of the distribution of the solid along with 1D simulation domain (dashed lines in Figure 5 and Figure 6). given in Table 1. 0.3). The simulation case is defined in Figure 7 and Table 2. simulations. Black curves are isolines of solid fraction (0.1; 0.2; 0.3). The simulation case is defined in Figure 7 and Table 2. 3. (d) 300 s. Iso-surfaces are displayed for 3, 4, 7, 9, and 10 wt%Pb. Iso-lines on the plane at height 20 mm from the bottom are from 3 to 4.2 wt%Pb with equiinterval of 0.1 wt%Pb). The simulation case is defined in Figure 11 and Table 3. 4.

List of figure captions

Figure 15 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy.

Calculated velocities of (a-c) the liquid phase and (d-f) the solid phase at process times (a, d) 10 s, (b, e) 100 s and (c, f) 200 s. In the vertical longitudinal crosssection, vectors indicate the velocity direction, while their color reflects the velocity magnitude. In the four horizontal transverse sections, the maps present the distribution of tangential velocities. The pink surface represents the packing surface at the characteristic solid fraction 0.4. The simulation case is defined in Figure 14 and Table 4.

Figure 16 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy.

Calculated solid fractions (a-c) and average solute composition (d-f) at process times (a, d) 10 s, (b, e) 100 s and (c, f) 200 s. The pink surface represents the packing surface at solid fraction 0.4. The simulation case is defined in Figure 14 and Table 4. (orange curve) 2D cylindrical FVM considering a spherical globular solid grains and (green curve) 2D cylindrical FVM simulation accounting for a dendritic morphology of the solid phase [18]. The simulation case is defined in Figure 14 and Table 4.

Introduction

The casting industry commonly faces difficulties in the production of products free from macrosegregation [1]. Microsegregation naturally takes place during solidification as a redistribution of the chemical species at the solid-liquid interface occurs, the result of thermodynamic equilibrium. Various processes take place which enhance segregation far from its interfacial origin including long range diffusion, liquid flow due to shrinkage, melt convection, solid deformation and transport of solid grains/fragments in the casting. The later induce macrosegregation from both solid and liquid transport over long distances. While these phenomena are difficult to avoid during conventional casting of metallic alloys, their magnitude needs to be controlled. The local average composition defined in a small representative volume must not deviate from the nominal composition of the alloy by more than a few percent, otherwise the properties could vary significantly and subsequent thermomechanical heat treatments may not be able to restore the desired properties everywhere in the product. This is true for various classes of metallic alloys, including large steel products [2]. Recently, efforts have been made to provide detailed experimental characterizations of macrosegregation in large steel ingots [3,4,5].

Numerical modeling of solidification accounting for the transport of the equiaxed grains at the process scale remains limited. Work was first done by Ni and Beckermann who proposed a volume-averaged model that consistently coupled microscopic phenomena with macroscopic transport [6,7]. Other solidification models, based on similar principles, have also been developed [8][9][10][11][12][13][14][15][16][17]. In volume-averaged models the transport of solid equiaxed grains is described by grain population balances along with mass and solute mass conservation equations for the solid phase. These equations consist of contributions from advection at the macroscopic (process) scale and grain growth governed by diffusion and phase change phenomena at the microscopic (grain) scale. These are strongly coupled with the transport of heat, mass, chemical species and momentum in the liquid and solid phases. Identifying the complexity in solving the set of coupled equations, Založnik and Combeau [8] proposed an operator splitting scheme as a flexible method for integration of the macroscopic transport terms and the local growth terms. This method was successfully implemented to simulate large ingot casting [18]. Modeling and simulation of steel ingots is particularly demanding due to the size of the castings and the complexity of the multiphase flow, however, models and applications have been improving in recent years [19][20][21][22][23]. These models are based on the finite volume method (FVM), while the finite element method (FEM) has not yet been considered for volume-averaged multiscale modelling of solidification with transport of equiaxed grains. An implementation using FEM may be attractive for multiple reasons. First, FEM generally offers more flexibility and versatility than FVM in describing the boundaries of the domain to be analyzed, and defining the boundary conditions which prevail there. This is particularly true when considering structured FVM where the "staircase" effect along the boundaries is detrimental. FEM also offers opportunities to more simply couple with stress/strain structural analyses, eg. to model the occurrence of thermomechanical defects in solidified regions, as such analyses are generally also conducted using FEM [24,25].

phenomena as well as for the motion of solid and liquid phases. First, the macroscopic conservation equations and the constitutive relations describing interfacial interactions and exchanges are summarized. The resolution method and the numerical implementation for the set of non-linear equations is then detailed. This implies introducing an artificial diffusion which deals with the discontinuities at the packing front by FEM. Numerical simulation using the current model is then performed to validate the numerical implementation. Results demonstrate an efficient FEM resolution scheme implemented for the purely convective transport problem, which is difficult to solve numerically by FEM in the absence of diffusive effects.

Two-phase model of solidification

This section summarizes the governing equations taken into consideration in this case when modeling solidification in the presence of solid transport. The full nomenclature and a complete set of notations can be found in Appendix A. Eqs. ( 1)-( 4) represent the conservation of total mass, momentum, solute mass and energy, respectively, for each phase 𝛼 (𝛼 = 𝑠 for the solid or 𝛼 = 𝑙 for the liquid phase) [8]. These macroscopic conservation equations are obtained by averaging the microscopic equations over a representative elementary volume (REV). A detailed derivation of the governing equations can be found elsewhere [6,26].

Mass 𝜕 𝜕𝑡 (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 ) + ∇ • (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 ) = 𝛤 𝛼 + 𝛷 𝛼 (1)
Momentum 𝜕 𝜕𝑡 (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 ) + ∇ • (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 ⊗ 〈𝐯 𝛼 〉 𝛼 )

= -𝑔 𝛼 ∇〈𝑝 𝛼 〉 𝛼 + ∇ • (𝑔 𝛼 〈𝛕 𝛼 〉 𝛼 ) + 𝑔 𝛼 〈𝐛 𝛼 〉 𝛼 + 𝐌 𝛤,𝛼 + 𝐌 𝑑,𝛼 + 𝐌 𝛷,𝛼 Species 𝜕 𝜕𝑡 (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝑤 𝛼 〉 𝛼 ) + ∇ • (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈𝑤 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 )

= -∇ • (𝑔 𝛼 〈𝐣 𝛼 〉 𝛼 ) + 𝐽 𝛤,𝛼 + 𝐽 𝑗,𝛼 + 𝐽 𝛷,𝛼

Energy 𝜕 𝜕𝑡 (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈ℎ 𝛼 〉 𝛼 ) + ∇ • (𝑔 𝛼 〈𝜌 𝛼 〉 𝛼 〈ℎ 𝛼 〉 𝛼 〈𝐯 𝛼 〉 𝛼 ) (4)
= -∇ • (𝑔 𝛼 〈𝐪 𝛼 〉 𝛼 ) + 𝑄 𝛤,𝛼 + 𝑄 𝑞,𝛼 + 𝑄 𝛷,𝛼

In the above equations ( 1)-( 4

The population of grains is described by an average density per unit volume, N. The population balance equation writes:

𝜕𝑁 𝜕𝑡 + ∇ • (𝑁〈𝐯 𝑠 〉 𝑠 ) = 𝑁 ̇ (7)
where the nucleation rate is defined as:

𝑁 ̇= { 𝑁 0 𝛿𝑡 if (𝑇 ≤ 𝑇 𝑛𝑢𝑐𝑙 ) and (𝑁 = 0 or ∫ 𝑁 ̇𝑑𝑡 = 0 𝑡 0 ) 0 otherwise ( 8
)
𝑁 0 is the nucleation density, 𝛿𝑡 the time step, and 𝑇 𝑛𝑢𝑐𝑙 the nucleation temperature. In Eq. ( 8) nucleation occurs under two conditions: The first nucleation event occurs when the local temperature drops below the nucleation temperature for the first time. Further nucleation events follow if the local grain density drops to zero (due to grain transport or remelting) and the local temperature is below the nucleation temperature. The latter criterion is a heuristic treatment used to prevent unphysical results [11].

In order to resolve the coupled microscopic and macroscopic phenomena, the microscopic exchange terms can be considered as source terms in the governing equations. The equations can then be integrated directly in a coupled way. However, due to the stiffness of the microscopic terms, such a solution scheme requires very small time steps. This makes the computation time realistically too long to simulate industrial sized castings. In the present work, an operator splitting method is applied to solve the entire system of equations [8].

According to this method, the phase fractions and solute compositions are considered to evolve in two subsequent stages: the macroscopic transport stage and the microscopic growth stage. When using this assumption, each stage can be numerically integrated by operating on different scales of time and space. Therefore, such a splitting technique provides an efficient way to solve multi-scale problems. This splitting method is used for the solution of the mass conservation equation for the solid, the conservation equation for the grain population density and the solute conservation equations for both phases. The method is summarized as follows:

first, in the macroscopic transport stage, only variation due to the macroscopic transport is integrated, which is determined by solving Eqs.( 9)-( 11) on the global finite element mesh using a macro time step. This gives an intermediate quantity with index tr:

𝜕𝑔 𝑡𝑟 𝑠 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝑠 〈𝐯 𝑠 〉 𝑠 ) = 0 (9)
𝜕𝑁 𝑡𝑟 𝜕𝑡 + ∇ • (𝑁 𝑡𝑟 〈𝐯 𝑠 〉 𝑠 ) = 0 (10)

𝜕(𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 ) 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 〈𝐯 𝛼 〉 𝛼 ) = 0 (11)
Second, in the microscopic growth stage, the contribution of microscopic processes, nucleation and growth, are integrated through Eqs.( 12)-( 14). These equations are solved, locally, at each node of the finite element mesh, leading to quantity with index gr:

𝜕𝑔 𝑔𝑟 𝑠 𝜕𝑡 = 𝛤 𝑠 + 𝛷 𝑠 𝜌 𝑠 ( 12
)
𝜕𝑁 𝑔𝑟 𝜕𝑡 = 𝑁 ̇ (13) 𝜕(𝑔 𝑔𝑟 𝛼 〈𝑤 𝛼 〉 𝑔𝑟 𝛼 ) 𝜕𝑡 = 𝐽 𝛤,𝛼 + 𝐽 𝑗,𝛼 + 𝐽 𝛷,𝛼 𝜌 𝛼 (14)
A smaller time step (micro time step) must be used for the microscopic growth stage. The solution of the macroscopic transport stage is the initial condition for the integration of the microscopic growth stage. The sequence of both integration steps thus gives the solution over a macroscopic time step. The whole modeling algorithm, using this splitting approach, is shown schematically in Regarding the transport stage, Eqs. ( 9)-( 11) are of pure convective nature, and notoriously difficult to solve numerically in the absence of diffusive effects. Moreover, another numerical difficulty arises from the discontinuity of transport velocities due to the solid packing phenomenon. Indeed, when forming a packed solid layer, solid grains suddenly change from a moving state to a fixed state. Solving these issues in the framework of FEM for solidification

simulations remains an open issue and will be addressed in the following sections.

Treatment of discontinuities at the packing front by the FEM

The solution of convection-dominated equations by FEM can encounter problems of unphysical oscillations, especially in zones with steep gradients. Different resolution techniques have been developed to overcome such issues: e.g. discretizing by upwind methods [27][28][29], by stabilized methods such as Streamline Upwind Petrov-Galerkin [30],

Galerkin Least-Squares [31,32], Residual Free Bubbles [33], or by using a corrected flux approach [34]. A specific complexity in the resolution of Eqs. ( 9)-( 11) arises from discontinuities due to the solid packing phenomenon taking place at the interface between the liquid pool and the packed bed of solid grains. On one side, small grains move freely in the liquid phase, the solid fraction being smaller than a characteristic packing fraction. On the other side, grains accumulate and grow to form a steady and fixed packed zone, the solid fraction being higher than the packing limit. Discontinuities are then related to abrupt changes in the velocity and fraction of the solid phase. Therefore, an adaptive artificial diffusion is introduced to stabilize the finite element resolutions without unreasonably smearing results.

This added diffusion detects and reduces discontinuities at locations where the solid phase is being packed. The added diffusion is then not present everywhere, it is restricted to critical zones of packing. The diffusion coefficient, 𝐷 𝑀 , as expressed below, consists of the gradient of solid velocities, i.e. a combination of the divergence of the average solid velocity ∇ • 〈𝐯 𝑠 〉 and the divergence of the intrinsic average solid velocity ∇ • (〈𝐯 𝑠 〉 𝑠 ), which contain information about the variation of the related quantities: velocity and fraction of solid phase.

𝐷 𝑀 = (ℎ 𝑒 𝐯 𝑠 ) 2 (𝛼 |∇ • 〈𝐯 𝑠 〉| + 𝛽|∇ • 〈𝐯 𝑠 〉 𝑠 |) (15)
In this expression, ℎ 𝑒 𝐯 𝑠 denotes the characteristic size of element 𝑒 in the direction of the solid velocity 𝐯 𝑠 , as proposed in [24]: (16) where 𝜑 𝑖 is the interpolation function associated with node 𝑖, 𝐯 𝑐𝑒𝑛𝑡𝑒𝑟 𝑠 is the solid velocity at the center of the element, and 𝑁𝑛 is the number of nodes per element. The coefficients 𝛼 and 𝛽 allow control over the amount of diffusion in a direct way, thus offering more flexibility than an implicit diffusion introduced by the upwind method. An alternative option would be to use the divergence of the intrinsic solid velocity, ∇ • 〈𝐯 𝑠 〉 𝑠 , and the solid fraction gradient, ∇𝑔 𝑠 . However, the former expression is preferred because it does not exist without solid transport. Hence the model is still valid in cases without solid motion. The formulation of the coefficient 𝐷 𝑀 is inspired by the work of Cook and Cabot who developed an artificial nonlinear diffusion using the entropy gradient to treat issues associated with discontinuities of temperature and mass fraction in supersonic reacting flows [35]. For the aforementioned issues experienced when simulating solidification there has, as of yet, been no relevant investigations.

ℎ 𝑒 𝐯 𝑠 = 2‖𝐯 𝑐𝑒𝑛𝑡𝑒𝑟 𝑠 ‖ ∑ |𝐯 𝑐𝑒𝑛𝑡𝑒𝑟 𝑠 • ∇𝜑 𝑖 | 𝑁𝑛 𝑖
It should be noted that it is necessary to use the same artificial diffusion coefficient for all relevant transport equations in order to ensure consistency between the transport of related quantities and the conservation of mass and energy. Consequently, the transport equations for solid phase, grain density and solute, derived from Eqs. ( 9)- (11), are modeled by Eqs. ( 17)- (19). These equations are solved in the finite element framework using an implicit scheme in time and a weighted residual approach with P1 linear elements (triangles in 2D, tetrahedra in 3D) and SUPG stabilization method.

𝜕𝑔 𝑡𝑟 𝑠 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝑠 〈𝐯 𝑠 〉 𝑠 ) -∇ • (𝐷 𝑀 𝛻𝑔 𝑡𝑟 𝑠 ) = 0 ( 17
)
𝜕𝑁 𝑡𝑟 𝜕𝑡 + ∇ • (𝑁 𝑡𝑟 〈𝐯 𝑠 〉 𝑠 ) -∇ • (𝐷 𝑀 𝛻𝑁 𝑡𝑟 ) = 0 (18) 𝜕(𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 ) 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 〈𝐯 𝛼 〉 𝛼 ) -∇ • (𝐷 𝑀 𝛻(𝑔 𝑡𝑟 𝛼 〈𝑤 𝛼 〉 𝑡𝑟 𝛼 )) = 0 (19)

Formulation of Energy Conservation

Regarding heat transfer, the equation for energy conservation, Eq.( 4), for the solid (𝛼 = 𝑠)

and for the liquid (𝛼 = 𝑙) phases, assuming thermal equilibrium between both phases in the REV, and introducing the added diffusion, the following mixture energy equation can be established:

𝜌 [ 𝜕〈ℎ〉 𝜕𝑡 + ∇ • (𝑔 𝑡𝑟 𝑠 〈ℎ 𝑠 〉 𝑠 〈𝐯 𝑠 〉 𝑠 + 𝑔 𝑡𝑟 𝑙 〈ℎ 𝑙 〉 𝑙 〈𝐯 𝑙 〉 𝑙 ) -∇ • (𝐷 𝑀 ∇(𝑔 𝑡𝑟 𝑠 〈ℎ 𝑠 〉 𝑠 + 𝑔 𝑡𝑟 𝑙 〈ℎ 𝑙 〉 𝑙 ))] -∇ • (〈ҡ〉∇𝑇) = 0 ( 20
)
where the average enthalpy and thermal conductivity are defined by 〈ℎ〉 = 𝑔 𝑠 〈ℎ 𝑠 〉 𝑠 + 𝑔 𝑙 〈ℎ 𝑙 〉 𝑙 and 〈ҡ〉 = 𝑔 𝑠 〈ҡ 𝑠 〉 𝑠 + 𝑔 𝑙 〈ҡ 𝑙 〉 𝑙

Additional it is assumed that the densities of phases are constant and equal and that the heat diffusion follows the Fourier law. Using this method, the phase fractions in the advection terms should be taken as those calculated in the transport stage, so that mass conservation is maintained [8]. In the present work, the enthalpy formulation of the energy equation is used.

Like the preceding transport equations, the energy equation is solved in the framework of the finite element formulation with an implicit scheme for time integration as well as a weighted residual approach with P1 linear elements and SUPG stabilization method.

Formulation of Momentum Conservation

Several assumptions are made when solving the momentum equations, including:

i.

The phase densities are constant and equal, except for the buoyancy forces for which the liquid and solid densities are modeled as follows, respectively.

𝜌 𝐵 𝑙 = 𝜌 𝑟𝑒𝑓 [1 -𝛽 𝑇 (𝑇 -𝑇 𝑟𝑒𝑓 ) -𝛽 𝑤 (〈𝑤 𝑙 〉 𝑙 -𝑤 𝑟𝑒𝑓 )] (22)
𝜌 𝐵 𝑠 = 𝜌 𝑟𝑒𝑓 1 -𝛽 𝑠ℎ𝑟 (23) ii. The pressure is assumed to be identical in the solid and liquid phases.

〈𝑝 𝑙 〉 𝑙 = 〈𝑝 𝑠 〉 𝑠 = 𝑝 (24)
iii. The liquid behaves as a Newtonian fluid with a constant viscosity. Neglecting the interfacial momentum transfer due to phase change, the divergence of the average deviatoric stress tensor is modeled as:

∇ • 〈𝛕 𝑙 〉 = 𝜇 𝑙 ∇ • (∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) 𝑇 ) (25)
iv.

The momentum transfer due to nucleation and growth is considered negligible relative to other terms, the momentum balance at the solid-liquid interface is thus described as:

𝐌 𝑑,𝑙 + 𝐌 𝑑,𝑠 = 𝟎 (26)
Using these assumptions, the liquid momentum equation is derived as Eq. ( 27) from the general formulation for phase 𝛼, Eq. (2).

𝜌 [ 𝜕 𝜕𝑡 (𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇ • (𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ⊗ 〈𝐯 𝑙 〉 𝑙 )] = -𝑔 𝑙 ∇𝑝 + 𝜇 𝑙 ∇ • (∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) 𝑇 ) + 𝑔 𝑙 𝜌 𝐵 𝑙 𝐠 -𝐌 𝑟𝑒𝑔𝑖𝑚𝑒 𝑑,𝑠 (27)
where the momentum transfer due to the drag force is described as

𝐌 𝑟𝑒𝑔𝑖𝑚𝑒 𝑑,𝑠 = ℳ 𝑟𝑒𝑔𝑖𝑚𝑒 (𝑔 𝑙 ) 2 (〈𝐯 𝑙 〉 𝑙 -〈𝐯 𝑠 〉 𝑠 ) (28)
in which 3 in the slurry regime (𝑔 𝑠 < 𝑔 𝑐 𝑠 )

ℳ 𝑟𝑒𝑔𝑖𝑚𝑒 = { ℳ 𝑠𝑙𝑢𝑟𝑟𝑦 = 3𝑔 𝑠 𝜇 𝑙 𝐶 𝑑 𝑅𝑒 4(𝑑 𝑔 ) 2 (𝑔 𝑙 )
ℳ 𝑝𝑎𝑐𝑘𝑒𝑑 = 180(1 -𝑔 𝑙 ) 2 𝜇 𝑙 𝜆 2 2 (𝑔 𝑙 ) 3 in the packed -bed regime (𝑔 𝑠 ≥ 𝑔 𝑐 𝑠 ) (29)
where the solid-liquid interaction follows the model of Agarwal and O'Neill [36].

In the finite element context, solving the above equation, Eq. ( 27), and distinguishing the two regimes on a single discretized domain may cause numerical difficulties because of the assembly of dissimilar and irregular terms (i.e. slurry regime contributions and packed-bed regime contributions). In order to overcome this issue, a smoothing procedure was developed, similar to the one introduced by Plotkowski and Krane [15]. It consists in using a transition function, 𝛼 𝑡 , to switch the momentum transfer, 𝐌 𝑟𝑒𝑔𝑖𝑚𝑒 𝑑,𝑠

, between the two regimes. The following expression is proposed for this transition function, which provides a compromise between having a sufficiently sharp change of regimes while avoiding an abrupt switch between them:

𝛼 𝑡 = 1 -0.5(1 -𝑡𝑎𝑛ℎ[𝛼 0 (𝑔 𝑐 𝑠 -𝑔 𝑠 )]) (30)
where the value 𝛼 0 = 100 is chosen, as plotted in Figure 2. The liquid momentum equation can then be expressed as

𝜌 [ 𝜕 𝜕𝑡 (𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇ • (𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ⊗ 〈𝐯 𝑙 〉 𝑙 )] (31)
= -𝑔 𝑙 ∇𝑝 + 𝜇 𝑙 ∇ • (∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) + ∇(𝑔 𝑙 〈𝐯 𝑙 〉 𝑙 ) 𝑇 ) + 𝑔 𝑙 𝜌 𝐵 𝑙 𝐠 -[𝛼 𝑡 ℳ 𝑠𝑙𝑢𝑟𝑟𝑦 + (1 -𝛼 𝑡 )ℳ 𝑝𝑎𝑐𝑘𝑒𝑑 ](𝑔 𝑙 ) 2 (〈𝐯 𝑙 〉 𝑙 -〈𝐯 𝑠 〉 𝑠 )
Which is solved by a semi-implicit time integration and a weighted residual method, precisely by the P1/P1 velocity-pressure formulation stabilized with the Variational Multi Scale method [37]. In the equation for the conservation of solid phase momentum, the inertial and viscous terms are neglected [8]. Therefore, the solid momentum equation, which is only considered in the slurry regime, reduces to:

-𝑔 𝑠 ∇𝑝 + 𝑔 𝑠 𝜌 𝐵 𝑠 𝐠 + ℳ 𝑠𝑙𝑢𝑟𝑟𝑦 (𝑔 𝑙 ) 2 (〈𝐯 𝑙 〉 𝑙 -〈𝐯 𝑠 〉 𝑠 ) = 0 (32)
From this equation, the solid velocity can be calculated locally, at each node of the finite element mesh. In order to avoid a sudden change of the velocity when solid packing occurs, a gradual transition is applied, by using the same switching function 𝛼 𝑡 introduced in Eq. ( 30).

Therefore the solid velocity is calculated at each node by: Furthermore, it is necessary to ensure that moving solid grains do not penetrate the preexisting packed bed, to avoid the solid fraction at the packing front exceeding the packing fraction. The velocity of the grains is therefore adjusted so that they land smoothly on the packed bed, i.e. reaching 〈𝐯 𝑠 〉 𝑠 = 0 at the packing limit. The algorithm for this adjustment is presented in Figure 3. It consists first of the calculation of the distance that the grains would travel in the direction of the solid fraction gradient, moving with their current velocity during the time step. This settling distance is 〈𝑣 𝑠 〉 𝑝𝑟𝑜𝑗 𝑠 ∆𝑡, where 〈𝑣 𝑠 〉 𝑝𝑟𝑜𝑗 𝑠 is the projected solid velocity and ∆𝑡 the time step. Then the settling distance is compared to the distance between the grains and the packing limit, 𝑑 𝑃-𝐿𝑖𝑚𝑖𝑡 , defined by the packing solid fraction, 𝑔 𝑐 𝑠 . If the settling distance is higher, the velocity is corrected by the factor 𝛼 as shown in Figure 3. This correction prevents the grains from traveling beyond the packing limit.

〈𝐯 𝑠 〉 𝑠 = 𝛼 𝑡 (〈𝐯 𝑙 〉 𝑙 + 𝑔 𝑠 (𝜌 𝐵 𝑠 𝐠 -∇𝑝) ℳ 𝑠𝑙𝑢𝑟𝑟𝑦 (𝑔 𝑙 ) 2 ) (33)

Verification of the adaptive artificial diffusion method for particle packing

The model described above is sophisticated as it involves artificial diffusive terms in conservation equations, a smoothing procedure for the transition from a slurry regime to a packed bed regime in the momentum conservation, as well as an adjustment strategy for the convection of solid in the vicinity of the packed bed. Several studies are then needed to evaluate the numerical parameters of the model. The first of which is presented below for a 1D sedimentation configuration where a simple analytical solution can be derived.

Test Case Description

In order to test these new parameters, the model is first applied to simulate a one-dimensional (1D) sedimentation process, considering only transport and neglecting nucleation and growth processes. It consists in the settling of a predefined number of globular grains with equal and constant size at uniform and constant velocity. A schematic of the test is presented in Figure 4 and the 1D analytical solution is derived in Appendix D. The computational domain is two-dimensional (2D), with a width and length of 1 mm × 100 mm. There is no heat exchange through the boundaries of the domain. The velocities at the top and bottom faces as well as the normal velocities along the vertical walls are set to zero and a perfect slip condition is applied to the tangential velocities on the vertical walls, i.e. returning to the 1D configuration described in Appendix D. Initial conditions are given in Figure 4. The binary alloy Sn-5wt%Pb is considered, its thermophysical properties can be found in Appendix C [8]. The enthalpy is evaluated according to the solid fraction and the temperature. The simulation parameters are given in Table 1.

The present test case is defined to benefit from the simple analytical solution shown in Figure 5 and Figure 6 at several times with dashed lines. The downward velocity of the settling grains is arbitrarily imposed to be a constant value in the unpacked region, equal to 1 mm s -1 and directed toward the -y-axis. The corresponding upward liquid velocity in the unpacked region was computed and found to be constant, equal to 0.111 mm s -1 , in agreement with the total mass balance. When the fraction of solid reaches 0.3 in the packed bed the velocity of the phases falls to zero. As there is no solidification, the fraction of solid cannot increase further. Conservation of the initial mass of solid thus defines the height of the packed bed. It reaches 20 mm once settling is complete. Also considering the adiabatic boundary conditions for heat transfer as well as the absence of phase change by solidification/remelting, the temperature is expected to remain constant and uniform throughout the simulation domain, equal to its initial value 498 K (224°C) shown in Figure 5 and

Analysis of Simulation Results

The initial solid velocity is set to 1 mm s -1 for the vertical component in the downward direction. The liquid velocity is then determined by solving the momentum equation, Eq. ( 31).

The results are presented in Figure 5 and Figure 6. Those obtained from the numerical simulation are shown as solid lines while those calculated by the analytical solution are dashed lines. The sedimentation process is illustrated in Figure 5-a as profiles of solid volume fraction at t = 0 s, 10 s (before grains reach the bottom of the domain), 30 s (accumulation below the packing limit), and 60 s (end of packing). From the initial state, the solid grains fall downward while the liquid moves upward in such a way that the continuity equation is respected, shown in Figure 5-b and-c. Since the solid fraction in the two-phase region is 0.1, the liquid velocity is ten times smaller than the solid velocity. The first grains reach the bottom after 20 s and begin to accumulate until the packing fraction (chosen as 0.3) is reached. During the process, the change of the solid phase from the moving to the packed state has an impact on the liquid movement, which is revealed by the peaks of liquid velocity at the packing front. The sedimentation is complete after 60 s, although the solid and liquid velocities remain non-zero within a small layer where the gravitational force is balanced with the diffusive effect. This phenomenon is maintained due to the persistent gradient of the solid fraction at the transition interface between the solid packed bed and the solid-free region.

The analytical and numerical results show similar trends, however, there exist differences between the two solutions. The differences found within the transition zones are due to diffusive effects. The simulation results contain an inevitable numerical diffusion and in the current case, at the boundary of the packed bed, an additional artificial diffusion. Furthermore, it can be observed that the solid fraction in the packed zone exceeds the predefined packing value of 0.3. This over accumulation is related to purely numerical issues when packing the solid phase and it will be discussed further in the next section.

As grain motion also involves the transport of solute and heat, it is important to verify the consistency of all transported quantities. The profiles of the average composition in about 30% solid at 0.364 wt% Pb and 70% liquid at 5 wt% Pb. Furthermore, as expected for a pure transport phenomenon, the temperature does not change during this process. Only a slight deviation from the initial temperature, smaller than 1 K, can be seen at 60 s, as shown in Figure 5-e. A good overall conservation of all quantities is verified, the maximal relative errors for the global solute mass and energy being about 10 -6 and 10 -5 , respectively.

Effects of Artificial Diffusion

Two simulations were performed in order to study the impact of the artificial diffusion term, by separately assessing the effect of the constant parameters introduced with variable 𝐷 𝑀 . In the first case, the coefficients are 𝛼 = 1 and 𝛽 = 0 while in the second one, 𝛼 = 0 and 𝛽 = 1. These problems are overcome in Case 2, the results of which are shown in Figure 6-d, -e and -f. This reveals that the term related to the variation of the intrinsic velocity has a more important role in solving numerical singularities than the average velocity. An over-packing relative to the predefined fraction is also stated with a higher excessive quantity than when 𝛼 = 20; 𝛽 = 1. This observation further shows the role of artificial diffusion in treating numerical issues related to the packing of solid, since the higher the added diffusion the better the simulation respects the predefined packing critical solid fraction. However, using high values of 𝛼 and 𝛽 lead to an excess of diffusion and produce unphysical solutions.

This set of simulations confirms that numerical concerns of finite element resolution in the

Modeling Applications -Macrosegregation Simulations

The following section presents test cases, first a two-dimensional (2D) and then three-dimensional (3D) applications for the solidification of a small Sn-Pb ingot. The objective is to further test the model while maintaining small geometries. Finally, the model being validated, a simulation for a 3D steel ingot on the scale of an industrial process will be given in order to discuss the application to a real casting geometry and to identify future possible improvements of the model.

Two-dimensional Test Case

Description

The studied case, presented in Figure 7, is configured according to the Hebditch-Hunt benchmark study [38]. A 100 mm x 60 mm cavity contains the same Sn-5wt%Pb alloy as in the previous sections, initially in the liquid state at 499.15 K (226 °C). Cooling takes place from the left side, an environment at 25 K, with a heat transfer coefficient of 300 W m -2 K -1 .

The rest of the boundary is assumed adiabatic. The simulation is performed with the complete model, which accounts for nucleation and growth processes, and transport phenomena, shown schematically in Figure 1. Nucleation happens at sites where the liquid is cooled below the liquidus temperature and where there are no existing grains. A homogenous grain density of . Black curves are isolines of solid fraction (0.1; 0.2; 0.3). The simulation case is defined in Figure 7 and Table 2. When comparing FE and FV methods, the map produced at the end of solidification shows a similar tendency of segregation formation, including the negative segregation located in the left zone of the cavity and the positive channel formed at the bottom. However, the FE solution produces a less marked negative segregation, and a larger positive channel along the bottom wall. Additionally, in the upper zone of the cavity, the FE solution shows positive segregation near the upper-left corner and a slightly negative segregation nearby. In this region, the FV computation produces positively segregated channels with higher solute content, distributed horizontally and alternating with negatively segregated zones. Differences between these two results can be partly explained by the influence of numerical factors.

Results and discussion

Because of the extremely high nonlinearity of the problem, any differences in the numerical solution methods can lead to noticeable differences between the solutions. The artificial diffusion used in the FE method inevitably leads to smoothing of the macroscopic fields and thus to a smaller degree of segregation than in the FV solution. In addition, many other factors, including discretization schemes, iteration procedures, etc., can be the cause of the differences between the FE and the FV results. The reader interested in such aspects can refer to literature in the context of simulations with a fixed solid phase [39,40].

When considering the distribution of solid grains it is useful to look at the distribution of grain density in Figure 10 at different instants. At 10 s, 20 s, and 200 s, it can be first observed that a large number of grains are gathered in the packed layer. It can also be seen that the transport of crystals by liquid advection induces several zones of significant grain density in the slurry region. As solid grains are transported, such a heterogeneous distribution is expected.

Although the results obtained from the FE and FV methods are still very similar at t = 10 s, larger differences between the two solutions are perceived later on. Compared to the FV results, a higher grain density in the right region is predicted by the FE simulation. It is known that the number of grains is governed and influenced by different coupled processes, consisting of nucleation, transport, re-melting and re-nucleation mechanisms. Although the physical parameters and the numerical procedure for nucleation are identical in both simulations, the larger quantity of grains in the FE solution might be caused by a higher frequency of nucleation events at nodes that were emptied of grains because of transport or remelting. Nevertheless, further investigations should be carried out, as the sources of those differences still remain unconfirmed. Despite these differences, clear similarities between both results can be observed, including a high grain density in the lower-right corner and a low grain density in the upper-left corner at the end of solidification. 3). The simulation case is defined in Figure 7 and Table 2. 3.

Table 3. Simulation parameters for the 3D test case presented in Figure 11. Results show very similar trends. However, it can be observed that the distributions of variables are non-uniform within the thickness of the 3D geometry. In the solid fraction maps (Figure 12-a), across the horizontal planes, it can be seen that the solid zone advances further along the specimen walls (plane P3) than in the interior. Thus the iso-surface of solid fraction exhibits a concave shape in the central zone of the cavity. This phenomenon can be explained by considering the distributions of other quantities since all relevant variables are closely related. First, as seen on the maps of velocities (Figure 12-e and f), there is no motion of neither solid nor liquid along the walls due to the no-slip boundary conditions. Because the grains remain attached to the cavity walls and do not settle, the solid fraction is higher than in the interior of the cavity (plane P6). The solid and liquid phases are not constrained in the immediate vicinity of the external surfaces. With a high quantity of mobile solid grains, the settling velocity is thus stronger there than that in the central zone as observed on the solid velocity map (Figure 12-f). This solid motion then enhances the downward movement of the solute-enriched liquid phase. Comparison of the maximum velocities at 10 s is possible, showing (2D, solid phase) 40.4 mm s -1 versus (3D, solid phase) 41.0 mm s -1 and (2D, liquid phase) 39.6 mm s -1 versus (3D, solid phase) 31.5 mm s -1 . In both cases, the location of the highest velocities for the solid and liquid phases are very close. The solid velocity being directly computed from eq. ( 33) by neglecting the inertial and viscous terms, i.e. only accounting for the solid-liquid interaction through the transferred momentum due to interfacial stresses, maximum values are almost equal in the 2D and 3D simulations.

Results and discussion

However, a lower value is found in the 3D simulation for the maximum velocity of the liquid phase. This is due to the interaction of the liquid flow with the two largest cavity walls (plane P3 and its symmetric), not accounted for in the 2D approximation. Thus, the transport of the liquid phase by the solid phase, while being obviously present, plays a less important role in comparison to the 2D approximation.

It is also interesting to observe the distribution of the flow in the horizontal cross sections.

The maximum velocities for both phases are not observed at the symmetry plane of the cavity.

Instead they are localized at about 1/3 of the half-cavity thickness from the cavity surfaces.

This uneven distribution is enhanced when successively considering the cross sections from top (height 57 mm) to bottom (height 17 mm). 3.

Three-dimensional Simulation of an Industrial Steel Ingot

Description

The numerical model is now applied to predict macrosegregation in an industrially cast steel ingot as produced by the company Aubert & Duval. Physical features and alloying components of the casting are detailed in [18], while its schematic is illustrated in Figure 14a. The bottom and lower region of the steel ingot are cooled via a mold while the upper region is surrounded by an insulating refractory material. The top surface is covered by a layer of exothermal powder. The mold outer surface thermally exchanges with the environment by convection and radiation.

For the current study, the steel is considered as a binary alloy composed of an iron base and a 0.36wt% nominal carbon content, which plays a dominant role in determining buoyancy force compared to other chemical elements [18]. A simulation is performed on one quarter of a cylinder (0.3 m radius and 1.8 m height), an approximation of the octagonal cross-section of the real ingot, which is bounded by two symmetric planes (P2 and P3 in Figure 14-b The simulation is carried out with a uniform mesh size of 20 mm and a constant time step of 0.01 s. 4.

Results and discussion

Table 4. Simulation parameters for the ingot case presented in Figure 14.

As it can be observed in these figures, the flow descends along the cooled side and ascends along the centerline, resulting in a global circulation loop in the shape of an elongated torus.

The maximum velocity is around 150 mm/s along the cooled wall. The persistence of this circulation loop is remarkable and its flow direction is the reverse of the direction observed when solid transport is not taken into account. If the solid is assumed to be fixed to the mold, the flow is driven only by natural convection induced by density differences in the liquid. The density of the liquid depends on the temperature and on the chemical composition; the concentration of carbon is most important when compared to all solutes in the multicomponent steel [18]. In the mushy zone the liquid density decreases with decreasing temperature because the liquid is enriched in carbon as solidification progresses. This creates a lateral liquid density gradient from the cooled walls towards the core of the casting, which drives the thermosolutal natural convection. If the solid is fixed, lighter liquid ascends along the walls and heavier liquid descends in the core. When solid grains move, the flow is completely modified. Solid crystals, heavier than the liquid, sediment along the cooled walls and entrain the liquid, thus inducing a downward flow. This phenomenon leads to an overall flow loop descending at the solidification front and ascending in the ingot core. This clearly shows that besides the thermo-solutal effect the solid transport plays a significant role in the formation of the natural convective flow during solidification. Moreover, the motion of solid and liquid are strongly coupled; the moving phases can entrain one another via drag forces.

The predicted flow structure is clearly three-dimensional and is not axisymmetric, although a four-fold symmetry is implicitly assumed by the choice of the computational domain. (about 150 mms -1 ) along the cooled wall. Despite a higher solid fraction near the wall (which can be seen in the transverse cross-sections) there is no permanent solid layer attached to this cooled wall. Additionally, in the velocity maps in Figure 15 it can be seen that after sinking to the bottom along the cooled side, the mobile solid phase is transported towards the center zone by liquid flow. In this way, solid grains coming from the outer solidified region first accumulate at the center area and then extend to the side wall, resulting in a packed solid built up from the bottom side. Since solute-depleted grains settle and occupy the lower zone, the liquid enriched in solute is ejected upwards. This gives rise to the formation of a negative segregation cone in the lower zone of the ingot, as shown in Figure 17-a. This is a typical phenomenon experimentally found in steel ingots. cylindrical FVM considering a spherical globular solid grains and (green curve) 2D cylindrical FVM simulation accounting for a dendritic morphology of the solid phase [18].

The simulation case is defined in Figure 14 and Table 4.

Conclusions

In this study, a finite element solidification model which takes into account the transport of equiaxed grains is presented. This model consists of  the resolution of a set of highly nonlinear and strongly coupled equations over multiple scales in time and space, including those of energy, phase movement, phase transport, grain density transport, solute transport, nucleation and solid growth,  the coupling of the equations based on the operator splitting algorithm, previously developed by Založnik and Combeau [8], demonstrated as an effective way for the numerical resolution of the evolution of solidification structures in the growth stage and the transport stage.

Different issues make the finite element implementation challenging: 4.

𝐛
Figure 16 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy.

Calculated solid fractions (a-c) and average solute composition (d-f) at process times (a, d) 10 s, (b, e) 100 s and (c, f) 200 s. The pink surface represents the packing surface at solid fraction 0.4. The simulation case is defined in Figure 14 and Table 4. (orange curve) 2D cylindrical FVM considering a spherical globular solid grains and (green curve) 2D cylindrical FVM simulation accounting for a dendritic morphology of the solid phase [18]. The simulation case is defined in Figure 14 and Table 4 Packing solid fraction, .

[-]

Figure 1 .

 1 Figure 1. Schematic of the resolution algorithm using the splitting method.

  167

Figure 2 .

 2 Figure 2. Transition function 𝛼 𝑡 vs. solid fraction supposing a packing solid fraction 𝑔 𝑐 𝑠 = 0.3.

Figure 6 .

 6 Due to the formation of a packed bed of grains with solid fraction 0.3 in a liquid with intrinsic composition 5 wt% Pb, the average composition reaches 3.609 wt% Pb. Consequently, total mass conservation leads to an average composition above the packed bed in the initially twophase region equal to 5.464 wt% Pb.

Figure 4 .

 4 Figure 4. Schematic of the 1D sedimentation test showing the sample geometry and initial distribution of solid grains. Additional conditions and simulation parameters are given in Table1.

Figure 5 -

 5 d present consistent evolution during the sedimentation process. When solid grains settle (Figure5-b), the solute-rich liquid moves upward (Figure5-c) and fills the region left by the grains, leading to an increase of the average composition in the upper zone. The average composition of 3.609 wt% Pb in the bottom zone corresponds to the final state where there is

Figure 6 -

 6 Figure 6-a, -b and -c present respectively the vertical profiles of solid fraction, average composition and temperature along the sample height for Case 1. In this case, the solid accumulation cannot be simulated since numerical problems occur when solid grains reach the bottom boundary, consequently the temperature does not remain constant and uniform throughout the domain.

Figure 5 .Figure 6 .

 56 Figure 5. Vertical profiles at times 0 s, 10 s, 30 s and 60 s for (a) solid fraction, (b) solid velocity, (c) liquid velocity, (d) average composition and (e) temperature with parameter 𝛼 = 20; 𝛽 = 1. The analytical solution (dashed curves) is derived in Appendix D.

10 9

 9 grains m -3 is generated. The transport of heat, mass, and solute is due to the motion of the solid and liquid phases, which are controlled by both thermo-solutal convection and sedimentation. It is assumed that there is no phase movement on the sides of the cavity (sticking contact with the boundary). The two-dimensional computation is carried out on a non-structured triangular mesh with a mean mesh size of 1 mm (including 15 143 elements and 7 730 nodes) and 10 micro time steps per a constant macro time step of 0.01 s.

Figure 7 .

 7 Figure 7. Schematics of the 2D cavity test for Sn-5wt%Pb alloy solidification showing the geometry and initial values. Simulation parameters are given inTable 2.

Figure 8 Figure 8 -

 88 Figure 8 presents the results at t = 10 s, including (a) solid fraction map and solid velocity vectors; (b) temperature map and liquid velocity vectors. The three black isolines represent the solid fractions 0.1 (upper isoline), 0.2 (intermediate) and 0.3 (lower). It can be observed inFigure8-a that after nucleation along the left cooled wall, solid grains settle to the bottom under the combined effect of gravity and downward solutal convection. Some of the grains that are still small are directly transported towards the right wall, they then continue to be carried by the liquid and move upward to about mid-height of the specimen. Due to the transport and settling of solid grains, a layer of packed grains begins building up along the

Figure 9

 9 Figure9 displays the time evolution of the distribution of the average solute composition. As for Figure8, the FE results are presented in the left column and the corresponding FV solutions in the right column. After 10 s of solidification, negative segregation forms at the bottom, while a large area of the cavity still remains at the initial composition of 5 wt%. The rejection of Pb from the solid phase during solidification enriches the liquid, increasing its density. However, the bottom area is occupied by the solid phase, composed of solutedepleted grains, even denser than the solute-rich liquid. Accordingly, a negative segregation layer progressively builds up from the bottom side, observed from 10 to 20 s. Additionally, it can be observed that the transition between the slurry zone, which is at a solid fraction of ~ 0.1, and the packed layer at a solid fraction of around 0.3, is rather thin. Such a narrow layer also indicates a prompt transition between the two flow regimes. While the solid velocity in the stationary packed bed is zero, an inter-granular liquid flow through the permeable packed bed persists. This flow creates a semicircular anti-clockwise circulation that brings solute from the upper to the lower regions of the packed bed in the left part of the domain, and from the lower to the upper regions in the right part. The resulting macrosegregation can be seen at time 200 s in Figure9-c: the average composition map shows that a large area in the

Figure 9 .Figure 10 .

 910 Figure 9. Simulations for the 2D cavity test for Sn-5wt%Pb alloy solidification showing maps of average Pb composition (wt%) at (a) 10 s, (b) 20 s, (c) 200 s, and (d) the end of solidification for the finite element (FE) and the finite volume (FV) simulations. Blackcurves are isolines of solid fraction (0.1; 0.2; 0.3). The simulation case is defined in Figure7and Table2.

  developed to simulate macrosegregation in three dimensions. The case considered is an extension of the previous 2D case, giving the cavity a thickness of 10 mm in the third dimension. Due to symmetry, the computational domain occupies one half of the thickness of the specimen, as shown in Figure11-a. The two largest surfaces are the median plane with symmetry conditions (numbered 6 in the figure, further denoted P6), and the front wall with a no-slip condition (plane numbered 3 in the figure, further denoted P3). Heat is extracted from the left wall, others being assumed adiabatic. The calculation is performed on a nonstructured mesh with a uniform mesh size of 1 mm, (294 935 elements and 57 626 nodes) and using a constant macro time step of 0.01 s and five micro increments per macro time step.

Figure 11 .

 11 Figure 11. Schematics of the 3D cavity test for Sn-5wt%Pb alloy solidification showing the geometry and initial values. Simulation parameters are given in Table3.

Figure 12 Figure 10 -

 1210 Figure12 presents the different fields along three horizontal planes and the vertical median plane at time 10 s. Comparison is possible with 2D simulations presented in Figure8-a

Figure 12 -Figure 12 -Figure 12 -

 121212 Figure 12-d is coherent with the temperature map given in Figure 12-c and the position of the liquidus isosurface. The solid fraction thus remains very low in a large undercooled zone, as shown when comparing Figure 12-a with the liquidus positionor the nucleation front in Figure 12-d -. This is due to the very low driving force for growth at low undercooling.Consequently, liquid flow is possible in this undercooled region, and is stronger in the vicinity of the symmetry plane. Transport of Pb solute that segregates between the grains is also preferentially taking place in the vicinity of the symmetry plane, explaining the isocontour drawn in Figure12-b, concave along the vertical growth front and convex in the bottom right region of the cavity. Close to the cavity walls, the sedimentation of the grains is the main cause for the liquid flow, leading to liquid velocity higher than at the center of the cavity.

Figure 12 .

 12 Figure 12. Simulations for the 3D cavity test for Sn-5wt%Pb alloy solidification showing snapshots at 10 s in the vertical symmetry plane and in 3 horizontal transversal planes at heights 17 mm, 37 mm, and 57 mm from the bottom. Variables drawn are (a) solid fraction, (b) average composition, (c) temperature, (d) grain density, (e) liquid velocity, and (f) solid velocity. Black iso-lines in the planar representations are (a) 𝑔 𝑠 = {0.02 ; 0.03}, (b) 〈𝑤〉 = {5.04 ; 5.06 ; 5.08} wt%Pb, (c) 𝑇 = {494.15 ; 495.15 ; 497.15 } K ({221 ; 222 ; 224 } °C) and (d) 𝑁 = {2 ; 3} 10 8 grains m -3 . Iso-surfaces in the four top views are defined by (a) 𝑔 𝑠 = 0.01, (b) 〈𝑤〉 = 5.02 wt%Pb , (c) 𝑇 = 496.15 K (223 °C) , and (d) 𝑁 = 10 7 grains m -3 . Velocity vectors in the bottom views are only displayed in the symmetry plane. The simulation case is defined in Figure 11 and Table 3.

Figure 13

 13 Figure13 displays the evolution upon solidification of the iso-surfaces of the average composition. The formation of a segregated channel at the bottom of the cavity is revealed. It forms first around the central zone of the cavity where the movements of the solid and liquid phases are not limited, unlike those restricted on the external surfaces. Inside the cavity, there always exists, until the end of the solidification process, some liquid flow that circulates through the packed bed and transports the solute from the upper to the lower regions, although its intensity is significantly reduced in comparison with the slurry zone. Consequently, the free recirculation of inter-granular flow enriched in the heavy element Pb leads to the formation of a positively segregated channel with a high solute content in the interior domain.When the process advances, this channel becomes a preferential path for the liquid phase and continues to expand along the width of the cavity (following the direction of inter-granular flow) as well as to develop in the thickness of the cavity with a content decreasing progressively when approaching the lateral surface. The various distribution of solute composition, in turn, induces different rates of solidification in the cavity, resulting in subsequent heterogeneities of quantities in all three dimensions.

Figure 13 .

 13 Figure 13. Simulations for the 3D cavity test for Sn-5wt%Pb alloy solidification showing snapshots of the average composition at times (a) 150 s, (b) 200 s, (c) 250 s, and (d) 300 s. Iso-surfaces are displayed for 3, 4, 7, 9, and 10 wt%Pb. Iso-lines on the plane at height 20 mm from the bottom are from 3 to 4.2 wt%Pb with equi-interval of 0.1 wt%Pb). The simulation case is defined in Figure 11 and Table 3.

  Figure14-b) is assumed to be adiabatic as the thermal cooling is restrained by use of the exothermal powder layer. At the beginning of the simulation, the ingot is assumed to be already filled by the liquid alloy at 1776.15 K (1503 °C). The exterior temperature is modeled to be gradually changed during the process: beginning at 900 K (626.85 °C) until 2000 s, then imposed to be 300 K (26.85 °C) when the cooling time is over 4000 s, and decreased linearly with time during the intermediate period. Grain nucleation is modeled following an instantaneous nucleation law with an initial density of 10 9 grains m -3 and a nucleation undercooling of 10 -3 K. Crystals are assumed to be blocked when the solid fraction reaches a packing value of 0.4. In this investigation, only equiaxed spherical crystals are considered.

Figure 15

 15 Figure15 presents the velocities of liquid and solid phases in a vertical cross-section, for which the vectors indicate the velocity direction and the color reflects the velocity magnitude.The three upper sub-figures (a-c) are for the liquid phase and the lower (d-f) for the solid phase. The pink surface displays the packing limit interface at a solid fraction of 0.4, below which solid grains are blocked and piled up from the bottom to this interface. Additionally, the tangential component of velocities is also illustrated in four transversal cross-sections at 0.4, 0.8, 1.2 and 1.6 m from the ingot's bottom.

Figure 14 .

 14 Figure 14. Schematics of the 3.3 ton steel ingot cast by Aubert & Duval [18] presenting (a) dimensions and materials and (b) simulated settings. Simulation parameters are given in Table4.

  The 3D structure of the flow is indicated in the horizontal slices of Figure 15, showing the tangential velocity component for both phases. It is around one order of magnitude smaller than that of the vertical velocities. The 3D structure is even more clearly observable on the shape of the packing front (pink surface in Figure 15) and in the distribution of macrosegregation, shown in Figures 16-17. It is possible that a certain degree of destabilization is induced by the numerics due to the use of a relatively coarse mesh size (20 mm).

Figure 16

 16 Figure16 presents the distribution of solid fraction (a-c) and of average composition (d-f) at different instants (10, 100 and 200 s). It can be seen that in the beginning of the process (at time 10 s) the whole population of grains are transported and sediment at quite high speed

Figure 17 -

 17 Figure17-b) globular grain morphology with a 2D-FVM[18] is shown in Figure17-b. It reveals larger deviation from the measurements compared to the present 3D-FEM simulation.However, when compared with (green curve in Figure17-b) dendritic morphology of solid crystals performed with a 2D-FVM implementation[18], clear improvement is seen and the prediction of segregation comes closer to experimental data.

Figure 15 .

 15 Figure 15. 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy. Calculated velocities of (a-c) the liquid phase and (d-f) the solid phase at process times (a, d) 10 s, (b, e) 100 s and (c, f) 200 s. In the vertical longitudinal cross-section, vectors indicate the velocity direction, while their color reflects the velocity magnitude. In the four horizontal transverse sections, the maps present the distribution of tangential velocities. The pink surface represents the packing surface at the characteristic solid fraction 0.4. The simulation case is defined in Figure 14 and Table 4.

Figure 16 .Figure 17 .

 1617 Figure 16. 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy. Calculated solid fractions (a-c) and average solute composition (d-f) at process times (a, d) 10 s, (b, e) 100 s and (c, f) 200 s. The pink surface represents the packing surface at solid fraction 0.4. The simulation case is defined in Figure 14 and Table 4.

(C. 5 ) 5 (

 55 𝑅𝑒 =𝑔 𝑙 2𝑅 𝑔 ‖〈𝐯 𝑙 〉 𝑙 -〈𝐯 𝑠 〉 𝑠 ‖ 𝜈 𝑛(𝑅𝑒) = 2𝑅𝑒 0.28 + 4.65 3(𝑅𝑒 0.28 + 4.65)(C. 6)Solute diffusion length in the solid phase 𝛿 𝑠 = 𝑅 𝑔

Figure 1

 1 Figure 1Schematic of the resolution algorithm using the splitting method.

Figure 2

 2 Figure 2 Transition function 𝛼 𝑡 vs. solid fraction supposing a packing solid fraction 𝑔 𝑐 𝑠 = 0.3.

Figure 3

 3 Figure 3 Adjustment strategy for the velocity of convected grains in the vicinity of the packed bed.

Figure 4

 4 Figure 4 Schematics of the 1D sedimentation test showing the sample geometry and initial distribution of solid grains. Additional conditions and simulation parameters are

Figure 5

 5 Figure 5 Vertical profiles at times 0 s, 10 s, 30 s and 60 s for (a) solid fraction, (b) solid velocity, (c) liquid velocity, (d) average composition and (e) temperature with parameter α=20; β=1. The analytical solution (dashed curves) is derived in Appendix D.

Figure 6

 6 Figure 6 Vertical profiles at times 0 s, 10 s, 30 s and 60 s for (a, d) the solid fraction, (b, e) the average composition and (c, f) the temperature with (a-c) 𝛼 = 1; 𝛽 = 0 and (d-f) 𝛼 = 0; 𝛽 = 1. The analytical solution (dashed curves) is derived in Appendix D.

Figure 7

 7 Figure 7 Schematics of the 2D cavity test for Sn-5wt%Pb alloy solidification showing the geometry and initial values. Simulation parameters are given in Table 2.

Figure 8

 8 Figure 8 Simulations for the 2D cavity test for Sn-5wt%Pb alloy solidification showing snapshots at time 10 s with the present finite element model (left column, FE) and with a reference finite volume model [8] (right column, FV): (a) solid fraction , 𝑔 𝑠 , and intrinsic solid velocity vectors, 〈𝐯 𝑠 〉 𝑠 , (b) temperature, 𝑇, and intrinsic liquid velocity vectors, 〈𝐯 𝑙 〉 𝑙 . Black curves are isolines of solid fraction (0.1; 0.2;

Figure 9

 9 Figure 9 Simulations for the 2D cavity test for Sn-5wt%Pb alloy solidification showing maps of average Pb composition (wt%) at (a) 10 s, (b) 20 s, (c) 200 s, and (d) the end of solidification for the finite element (FE) and the finite volume (FV)

Figure 10

 10 Figure 10 Simulations for the 2D cavity test for Sn-5wt%Pb alloy solidification showing maps of grain density (m -3 ) at (a) 10 s, (b) 20 s, (c) 200 s, and (d) end of solidification for (FE) the finite element simulation and (FV) the finite volume simulation. The simulation case is defined in Figure 7 and Table 2.

Figure 11

 11 Figure 11 Schematics of the 3D cavity test for Sn-5wt%Pb alloy solidification showing the geometry and initial values. Simulation parameters are given in Table3.

Figure 12

 12 Figure 12 Simulations for the 3D cavity test for Sn-5wt%Pb alloy solidification showing snapshots at 10 s in the vertical symmetry plane and in 3 horizontal transversal planes at heights 17 mm, 37 mm, and 57 mm from the bottom. Variables drawn are (a) solid fraction, (b) average composition, (c) temperature, (d) grain density, (e) liquid velocity, and (f) solid velocity. Black iso-lines in the planar representations are (a) 𝑔 𝑠 = {0.02 ; 0.03}, (b) 〈𝑤〉 = {5.04 ; 5.06 ; 5.08} wt%Pb, (c) 𝑇 = {494.15 ; 495.15 ; 497.15 } K ({221 ; 222 ; 224 } °C) and (d) 𝑁 = {2 ; 3} 10 8 grains m -3 . Iso-surfaces in the four top views are defined by (a) 𝑔 𝑠 = 0.01, (b) 〈𝑤〉 = 5.02 wt%Pb , (c) 𝑇 = 496.15 K (223 °C) , and (d) 𝑁 = 10 7 grains m -3 . Velocity vectors in the bottom views are only displayed in the symmetry plane. The simulation case is defined in Figure 11 and Table3.

Figure 13

 13 Figure 13 Simulations for the 3D cavity test for Sn-5wt%Pb alloy solidification showing snapshots of the average composition at times (a) 150 s, (b) 200 s, (c) 250 s, and

Figure 14

 14 Figure 14 Schematics of the 3.3 ton steel ingot cast by Aubert & Duval [18] presenting (a) dimensions and materials and (b) simulated settings. Simulation parameters are given in Table4.

Figure 17 3D

 17 Figure 17 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy with (a) the final macrosegregation map and (b) segregation profiles along the center line with (red points) measurements, (blue curve) present 3D FEM simulation and

∑

  ), the notation 〈 𝛼 〉 𝛼 indicates the intrinsic volume average in phase 𝛼, 𝑔 denotes the volume fraction, 𝜌 the density, 𝐯 the velocity, 𝑝 the pressure, 𝛕 the deviatoric part of the stress tensor, 𝐛 the body force per unit volume, 𝑤 the solute mass concentration, 𝐣 the solute flux vector, ℎ the specific enthalpy, 𝐪 the heat flux vector. The right-hand side of the above equations gathers the exchange terms rising from different microscopic processes: 𝛤 denotes the mass exchange rate due to phase change, 𝛷 the mass exchange rate due to grain nucleation, 𝐌 the vector for interfacial momentum exchange, 𝐽 the solute exchange rate, and 𝑄 the heat exchange rate. The contributions of nucleation (terms with superscript 𝛷) can be neglected compared to other terms in theses equations. The microscopic exchange contributions are modeled as follows, more detail can be found in Appendix B. due to phase change due to interfacial stresses or diffusion 𝐌 𝛤,𝛼 = 𝐯 ̅ 𝛼 * 𝛤 𝛼 𝐌 𝑑,𝛼 = 𝜌 𝛼 𝑅 𝛼 𝑆 v (𝐯 ̅ 𝛼 * -〈𝐯 𝛼 〉 𝛼 ) 𝐽 𝛤,𝛼 = 𝑤 ̅ 𝛼 * 𝛤 𝛼 𝐽 𝑗,𝛼 = 𝜌 𝛼 𝐷 𝛼 𝛿 𝛼 𝑆 v (𝑤 ̅ 𝛼 * -〈𝑤 𝛼 〉 𝛼 ) 𝑄 𝛤,𝛼 = ℎ ̅ 𝛼 * 𝛤 𝛼 𝑄 𝑞,𝛼 = ҡ 𝛼 𝑙 𝛼 𝑆 v (𝑇 ̅ 𝛼 * -〈𝑇 𝛼 〉 𝛼 ) Here 𝐯 ̅ * , 𝑤 ̅ * , ℎ ̅ * , 𝑇 ̅ * are the average values over the interface, R is the momentum resistance coefficient, D is the solute diffusion coefficient,  is the characteristic solute diffusion length, ҡ is the heat conductivity, l is the characteristic heat conduction length, 𝑆 v = 𝐴 𝑉 𝑒 ⁄ is the interfacial area concentration in the REV considered, of volume 𝑉 𝑒 . The model is closed by balances of mass, momentum, solute, and heat at the solid-liquid interface: 𝐽 𝛤,𝛼 + 𝐽 𝑗,𝛼 + 𝐽 𝛷,𝛼 𝛼=𝑠,𝑙 = 0 ∑ 𝑄 𝛤,𝛼 + 𝑄 𝑞,𝛼 + 𝑄 𝛷,𝛼 𝛼=𝑠,𝑙 = 0

Figure 1 .

 1 Each iteration 𝜈 begins with an implicit finite element resolution for energy and transport conservation equations in the transport stage by using a macro time step. This gives the solutions denoted with superscript 𝜈 + 1 and subscript 𝑡𝑟 when they are associated with the transport stage. Then the nucleation-and-growth stage is solved locally (i.e., at each node), with variables initialized by values obtained from the transport stage at 𝜈 + 1. This local resolution proceeds through micro time steps, assuming that the average quantities for the solid plus liquid mixture, composition 〈𝑤〉 = 𝑔 𝑙 〈𝑤 𝑙 〉 𝑙 + 𝑔 𝑠 〈𝑤 𝑠 〉 𝑠 and enthalpy 〈ℎ〉 = 𝑔 𝑙 〈ℎ 𝑙 〉 𝑙 + 𝑔 𝑠 〈ℎ 𝑠 〉 𝑠 , no longer evolve during the macro time step. Finally, the momentum equations are solved with a semi-implicit solver on the macro time step to compute the new estimation of velocity fields of the liquid and the solid phases at iteration 𝜈 + 1. A complete evolution of the different quantities over the time step is evaluated by the final results obtained from the growth stage, as these solutions already include the change from the transport stage. The splitting scheme is only used to solve the evolution of phase fractions, grain density and solute concentrations since it involves very different scales of time and space. The resolution of the energy and momentum conservation equations do not require operator splitting because the constitutive relations coupling the micro-and macroscopic scales are simpler.

Figure 1 .

 1 Figure 1. Schematic of the resolution algorithm using the splitting method.

  167

Figure 2 .

 2 Figure 2. Transition function 𝛼 𝑡 vs. solid fraction supposing a packing solid fraction 𝑔 𝑐 𝑠 = 0.3.

Figure 3 .

 3 Figure 3. Adjustment strategy for the velocity of convected grains in the vicinity of the packed bed.

Figure 6 .

 6 Due to the formation of a packed bed of grains with solid fraction 0.3 in a liquid with intrinsic composition 5 wt% Pb, the average composition reaches 3.609 wt% Pb. Consequently, total mass conservation leads to an average composition above the packed bed in the initially twophase region equal to 5.464 wt% Pb.

Figure 4 .

 4 Figure 4. Schematic of the 1D sedimentation test showing the sample geometry and initial distribution of solid grains. Additional conditions and simulation parameters are given in Table1.

Figure 5 -

 5 d present consistent evolution during the sedimentation process. When solid grains settle (Figure5-b), the solute-rich liquid moves upward (Figure5-c) and fills the region left by the grains, leading to an increase of the average composition in the upper zone. The average composition of 3.609 wt% Pb in the bottom zone corresponds to the final state where there is

Figure 6 -

 6 Figure 6-a, -b and -c present respectively the vertical profiles of solid fraction, average composition and temperature along the sample height for Case 1. In this case, the solid accumulation cannot be simulated since numerical problems occur when solid grains reach the bottom boundary, consequently the temperature does not remain constant and uniform throughout the domain.

Figure 5 .Figure 6 .

 56 Figure 5. Vertical profiles at times 0 s, 10 s, 30 s and 60 s for (a) solid fraction, (b) solid velocity, (c) liquid velocity, (d) average composition and (e) temperature with parameter 𝛼 = 20; 𝛽 = 1. The analytical solution (dashed curves) is derived in Appendix D.

10 9

 9 grains m -3 is generated. The transport of heat, mass, and solute is due to the motion of the solid and liquid phases, which are controlled by both thermo-solutal convection and sedimentation. It is assumed that there is no phase movement on the sides of the cavity (sticking contact with the boundary). The two-dimensional computation is carried out on a non-structured triangular mesh with a mean mesh size of 1 mm (including 15 143 elements and 7 730 nodes) and 10 micro time steps per a constant macro time step of 0.01 s.

Figure 7 .

 7 Figure 7. Schematics of the 2D cavity test for Sn-5wt%Pb alloy solidification showing the geometry and initial values. Simulation parameters are given inTable 2.

Figure 8 .

 8 Figure 8. Simulations for the 2D cavity test for Sn-5wt%Pb alloy solidification showing snapshots at time 10 s with the present finite element model (left column, FE) and with a reference finite volume model [8] (right column, FV): (a) solid fraction , 𝑔 𝑠 , and intrinsic solid velocity vectors, 〈𝐯 𝑠 〉 𝑠 , (b) temperature, 𝑇, and intrinsic liquid velocity vectors, 〈𝐯 𝑙 〉 𝑙. Black curves are isolines of solid fraction (0.1; 0.2; 0.3). The simulation case is defined in Figure7and Table2.

Figure 8 Figure 8 -

 88 Figure 8 presents the results at t = 10 s, including (a) solid fraction map and solid velocity vectors; (b) temperature map and liquid velocity vectors. The three black isolines represent the solid fractions 0.1 (upper isoline), 0.2 (intermediate) and 0.3 (lower). It can be observed inFigure8-a that after nucleation along the left cooled wall, solid grains settle to the bottom under the combined effect of gravity and downward solutal convection. Some of the grains that are still small are directly transported towards the right wall, they then continue to be carried by the liquid and move upward to about mid-height of the specimen. Due to the transport and settling of solid grains, a layer of packed grains begins building up along the

Figure 9

 9 Figure9 displays the time evolution of the distribution of the average solute composition. As for Figure8, the FE results are presented in the left column and the corresponding FV solutions in the right column. After 10 s of solidification, negative segregation forms at the bottom, while a large area of the cavity still remains at the initial composition of 5 wt%. The rejection of Pb from the solid phase during solidification enriches the liquid, increasing its density. However, the bottom area is occupied by the solid phase, composed of solutedepleted grains, even denser than the solute-rich liquid. Accordingly, a negative segregation layer progressively builds up from the bottom side, observed from 10 to 20 s. Additionally, it can be observed that the transition between the slurry zone, which is at a solid fraction of ~ 0.1, and the packed layer at a solid fraction of around 0.3, is rather thin. Such a narrow layer also indicates a prompt transition between the two flow regimes. While the solid velocity in the stationary packed bed is zero, an inter-granular liquid flow through the permeable packed bed persists. This flow creates a semicircular anti-clockwise circulation that brings solute from the upper to the lower regions of the packed bed in the left part of the domain, and from the lower to the upper regions in the right part. The resulting macrosegregation can be seen at time 200 s in Figure9-c: the average composition map shows that a large area in the stationary mushy zone has a negative segregation which results from accumulation of solutedepleted grains. However, the average composition is not uniform. At the very bottom of the cavity, there is an accumulation of solute which results from intergranular melt flow localized along the bottom wall and oriented in the direction of the temperature gradient. In the left part of the packed layer the negative segregation is amplified by the flow of intergranular liquid, which is oriented against the temperature gradient in this region. As the process advances, the channel continues to extend along the bottom side.

Figure 9 .

 9 Figure 9. Simulations for the 2D cavity test for Sn-5wt%Pb alloy solidification showing maps of average Pb composition (wt%) at (a) 10 s, (b) 20 s, (c) 200 s, and (d) the end of solidification for the finite element (FE) and the finite volume (FV) simulations. Blackcurves are isolines of solid fraction (0.1; 0.2; 0.3). The simulation case is defined in Figure7and Table2.

Figure 10 .

 10 Figure 10. Simulations for the 2D cavity test for Sn-5wt%Pb alloy solidification showing maps of grain density (m -3 ) at (a) 10 s, (b) 20 s, (c) 200 s, and (d) end of solidification for (FE) the finite element simulation and (FV) the finite volume simulation. The simulation case is defined in Figure 7 and Table 2.

  developed to simulate macrosegregation in three dimensions. The case considered is an extension of the previous 2D case, giving the cavity a thickness of 10 mm in the third dimension. Due to symmetry, the computational domain occupies one half of the thickness of the specimen, as shown in Figure11-a. The two largest surfaces are the median plane with symmetry conditions (numbered 6 in the figure, further denoted P6), and the front wall with a no-slip condition (plane numbered 3 in the figure, further denoted P3). Heat is extracted from the left wall, others being assumed adiabatic. The calculation is performed on a nonstructured mesh with a uniform mesh size of 1 mm, (294 935 elements and 57 626 nodes) and using a constant macro time step of 0.01 s and five micro increments per macro time step.

Figure 11 .

 11 Figure 11. Schematics of the 3D cavity test for Sn-5wt%Pb alloy solidification showing the geometry and initial values. Simulation parameters are given in Table3.

Figure 12 Figure 10 -

 1210 Figure12 presents the different fields along three horizontal planes and the vertical median plane at time 10 s. Comparison is possible with 2D simulations presented in Figure8-a

Figure 12 -Figure 12 -Figure 12 -

 121212 Figure 12-d is coherent with the temperature map given in Figure 12-c and the position of the liquidus isosurface. The solid fraction thus remains very low in a large undercooled zone, as shown when comparing Figure 12-a with the liquidus positionor the nucleation front in Figure 12-d -. This is due to the very low driving force for growth at low undercooling.Consequently, liquid flow is possible in this undercooled region, and is stronger in the vicinity of the symmetry plane. Transport of Pb solute that segregates between the grains is also preferentially taking place in the vicinity of the symmetry plane, explaining the isocontour drawn in Figure12-b, concave along the vertical growth front and convex in the bottom right region of the cavity. Close to the cavity walls, the sedimentation of the grains is the main cause for the liquid flow, leading to liquid velocity higher than at the center of the cavity.

Figure 12 .

 12 Figure 12. Simulations for the 3D cavity test for Sn-5wt%Pb alloy solidification showing snapshots at 10 s in the vertical symmetry plane and in 3 horizontal transversal planes at heights 17 mm, 37 mm, and 57 mm from the bottom. Variables drawn are (a) solid fraction, (b) average composition, (c) temperature, (d) grain density, (e) liquid velocity, and (f) solid velocity. Black iso-lines in the planar representations are (a) 𝑔 𝑠 = {0.02 ; 0.03}, (b) 〈𝑤〉 = {5.04 ; 5.06 ; 5.08} wt%Pb, (c) 𝑇 = {494.15 ; 495.15 ; 497.15 } K ({221 ; 222 ; 224 } °C) and (d) 𝑁 = {2 ; 3} 10 8 grains m -3 . Iso-surfaces in the four top views are defined by (a) 𝑔 𝑠 = 0.01, (b) 〈𝑤〉 = 5.02 wt%Pb , (c) 𝑇 = 496.15 K (223 °C) , and (d) 𝑁 = 10 7 grains m -3 . Velocity vectors in the bottom views are only displayed in the symmetry plane. The simulation case is defined in Figure 11 and Table 3.

Figure 13

 13 Figure13 displays the evolution upon solidification of the iso-surfaces of the average composition. The formation of a segregated channel at the bottom of the cavity is revealed. It forms first around the central zone of the cavity where the movements of the solid and liquid phases are not limited, unlike those restricted on the external surfaces. Inside the cavity, there always exists, until the end of the solidification process, some liquid flow that circulates through the packed bed and transports the solute from the upper to the lower regions, although its intensity is significantly reduced in comparison with the slurry zone. Consequently, the free recirculation of inter-granular flow enriched in the heavy element Pb leads to the formation of a positively segregated channel with a high solute content in the interior domain.When the process advances, this channel becomes a preferential path for the liquid phase and continues to expand along the width of the cavity (following the direction of inter-granular flow) as well as to develop in the thickness of the cavity with a content decreasing progressively when approaching the lateral surface. The various distribution of solute composition, in turn, induces different rates of solidification in the cavity, resulting in subsequent heterogeneities of quantities in all three dimensions.

Figure 13 .

 13 Figure 13. Simulations for the 3D cavity test for Sn-5wt%Pb alloy solidification showing snapshots of the average composition at times (a) 150 s, (b) 200 s, (c) 250 s, and (d) 300 s. Iso-surfaces are displayed for 3, 4, 7, 9, and 10 wt%Pb. Iso-lines on the plane at height 20 mm from the bottom are from 3 to 4.2 wt%Pb with equi-interval of 0.1 wt%Pb). The simulation case is defined in Figure 11 and Table 3.

Figure 14 -

 14 Figure14-b) is assumed to be adiabatic as the thermal cooling is restrained by use of the exothermal powder layer. At the beginning of the simulation, the ingot is assumed to be already filled by the liquid alloy at 1776.15 K (1503 °C). The exterior temperature is modeled to be gradually changed during the process: beginning at 900 K (626.85 °C) until 2000 s, then imposed to be 300 K (26.85 °C) when the cooling time is over 4000 s, and decreased linearly with time during the intermediate period. Grain nucleation is modeled following an instantaneous nucleation law with an initial density of 10 9 grains m -3 and a nucleation undercooling of 10 -3 K. Crystals are assumed to be blocked when the solid fraction reaches a packing value of 0.4. In this investigation, only equiaxed spherical crystals are considered.

Figure 15

 15 Figure15 presents the velocities of liquid and solid phases in a vertical cross-section, for which the vectors indicate the velocity direction and the color reflects the velocity magnitude.The three upper sub-figures (a-c) are for the liquid phase and the lower (d-f) for the solid phase. The pink surface displays the packing limit interface at a solid fraction of 0.4, below which solid grains are blocked and piled up from the bottom to this interface. Additionally, the tangential component of velocities is also illustrated in four transversal cross-sections at 0.4, 0.8, 1.2 and 1.6 m from the ingot's bottom.

Figure 14 .

 14 Figure 14. Schematics of the 3.3 ton steel ingot cast by Aubert & Duval [18] presenting (a) dimensions and materials and (b) simulated settings. Simulation parameters are given in Table4.

  The 3D structure of the flow is indicated in the horizontal slices of Figure 15, showing the tangential velocity component for both phases. It is around one order of magnitude smaller than that of the vertical velocities. The 3D structure is even more clearly observable on the shape of the packing front (pink surface in Figure 15) and in the distribution of macrosegregation, shown in Figures 16-17. It is possible that a certain degree of destabilization is induced by the numerics due to the use of a relatively coarse mesh size (20 mm).

Figure 16

 16 Figure16 presents the distribution of solid fraction (a-c) and of average composition (d-f) at different instants (10, 100 and 200 s). It can be seen that in the beginning of the process (at time 10 s) the whole population of grains are transported and sediment at quite high speed

Figure 17 -Figure 17 -

 1717 Figure17-b) globular grain morphology with a 2D-FVM[18] is shown in Figure17-b. It reveals larger deviation from the measurements compared to the present 3D-FEM simulation.However, when compared with (green curve in Figure17-b) dendritic morphology of solid crystals performed with a 2D-FVM implementation[18], clear improvement is seen and the prediction of segregation comes closer to experimental data.

Figure 15 .

 15 Figure 15. 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy. Calculated velocities of (a-c) the liquid phase and (d-f) the solid phase at process times (a, d) 10 s, (b, e) 100 s and (c, f) 200 s. In the vertical longitudinal cross-section, vectors indicate the velocity direction, while their color reflects the velocity magnitude. In the four horizontal transverse sections, the maps present the distribution of tangential velocities. The pink surface represents the packing surface at the characteristic solid fraction 0.4. The simulation case is defined in Figure 14 and Table 4.

  681

Figure 16 .

 16 Figure 16. 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy. Calculated solid fractions (a-c) and average solute composition (d-f) at process times (a, d) 10 s, (b, e) 100 s and (c, f) 200 s. The pink surface represents the packing surface at solid fraction 0.4. The simulation case is defined in Figure 14 and Table 4.

  682

Figure 17 .

 17 Figure 17. 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy with (a) the final macrosegregation map and (b) segregation profiles along the center line with (red points) measurements, (blue curve) present 3D FEM simulation and (orange curve) 2Dcylindrical FVM considering a spherical globular solid grains and (green curve) 2D cylindrical FVM simulation accounting for a dendritic morphology of the solid phase[18].The simulation case is defined in Figure14and Table4.

Figure 12

 12 Figure 12 Simulations for the 3D cavity test for Sn-5wt%Pb alloy solidification showing snapshots at 10 s in the vertical symmetry plane and in 3 horizontal transversal planes at heights 17 mm, 37 mm, and 57 mm from the bottom. Variables drawn are (a) solid fraction, (b) average composition, (c) temperature, (d) grain density, (e) liquid velocity, and (f) solid velocity. Black iso-lines in the planar representations are (a) 𝑔 𝑠 = {0.02 ; 0.03}, (b) 〈𝑤〉 = {5.04 ; 5.06 ; 5.08} wt%Pb, (c) 𝑇 = {494.15 ; 495.15 ; 497.15 } K ({221 ; 222 ; 224 } °C) and (d) 𝑁 = {2 ; 3} 10 8 grains m -3 . Iso-surfaces in the four top views are defined by (a) 𝑔 𝑠 = 0.01, (b) 〈𝑤〉 = 5.02 wt%Pb , (c) 𝑇 = 496.15 K (223 °C) , and (d) 𝑁 = 10 7 grains m -3 . Velocity vectors in the bottom views are only displayed in the symmetry plane. The simulation case is defined in Figure 11 and Table3.

Figure 13

 13 Figure 13 Simulations for the 3D cavity test for Sn-5wt%Pb alloy solidification showing snapshots of the average composition at times (a) 150 s, (b) 200 s, (c) 250 s, and(d) 300 s. Iso-surfaces are displayed for 3, 4, 7, 9, and 10 wt%Pb. Iso-lines on the plane at height 20 mm from the bottom are from 3 to 4.2 wt%Pb with equiinterval of 0.1 wt%Pb). The simulation case is defined in Figure11and Table3.

Figure 14

 14 Figure 14 Schematics of the 3.3 ton steel ingot cast by Aubert & Duval [18] presenting (a) dimensions and materials and (b) simulated settings. Simulation parameters are given in Table4.

Figure 15 3DFigure 14

 1514 Figure 15 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy. Calculated velocities of (a-c) the liquid phase and (d-f) the solid phase at process times (a, d) 10 s, (b, e) 100 s and (c, f) 200 s. In the vertical longitudinal crosssection, vectors indicate the velocity direction, while their color reflects the velocity magnitude. In the four horizontal transverse sections, the maps present the distribution of tangential velocities. The pink surface represents the packing surface at the characteristic solid fraction 0.4. The simulation case is defined in Figure 14 and Table4.

Figure 17 3D

 17 Figure 17 3D FEM simulation of the solidification of a binary Fe-0.36wt%C alloy with (a) the final macrosegregation map and (b) segregation profiles along the center line with (red points) measurements, (blue curve) present 3D FEM simulation and

[

  wt.%)-1 ] Solute diffusion coefficient at micro scale, in liquid,-[m 2 s -1 ]Solute diffusion coefficient at micro scale, in solid,

Table 1 .Table 1 .

 11 Simulation parameters for the 1D test case presented in

Figure 4.

Table 2 .Table 2 .

 22 Simulation parameters for the 2D test case presented in

Figure 7. Figure 8. Simulations for the 2D cavity test for Sn-5wt%Pb alloy solidification showing snapshots at time 10 s with the present finite element model (left column, FE) and with a reference finite volume model

Table 1 .Table 1 .

 11 Simulation parameters for the 1D test case presented in Figure 4.

Table 2 .Table 2 .

 22 Simulation parameters for the 2D test case presented in Figure 7.

Table 4

 4 .Density of the liquid phase,

	Mesh size Macro time step	.	.5 [mm] [s]
	Macro/Micro time ratio		[-]
			[-]
			[-]
	Table 1

Table D1

 D1

	Time	Interval		Solid fraction	Average composition
	t [s]	[mm]		[-]	[wt% Pb]
		[ , ]			5
		[ , ]	.	5
		[ ,	]		5
		[ , ]			5
		[ , ]	.1	4.5364
		[ , ]	.	5
		[ , ]		5.556
		[ ,	]		5
		[ , ]		.3	3.6092
		[ , ]		.	4.5364
		[ , ]	.	5
		[ , ]		5.556
		[ ,	]		5
		[ , ]		.3	3.6093
		[ , ]		5.556
		[ ,	]		5
				Table C1

Acknowledgements

The authors gratefully acknowledged the financial support to this study from the following industrial partners: ArcelorMittal, Aubert & Duval, AscoIndustries and Aperam. The finite volume simulations were performed with software SOLID developed at Institut Jean Lamour, Université de Lorraine, Nancy, France. The authors thank Laurent Heyvaert for his help with the simulations with SOLID as well as Jacob Kennedy for his careful reading of the manuscript.

Acknowledgements

The authors gratefully acknowledged the financial support to this study from the following industrial partners: ArcelorMittal, Aubert & Duval, AscoIndustries and Aperam. The finite volume simulations were performed with software SOLID developed at Institut Jean Lamour, Université de Lorraine, Nancy, France. The authors thank Laurent Heyvaert for his help with the simulations with SOLID as well as Jacob Kennedy for his careful reading of the manuscript.

Appendix C -Thermophysical Data

Table C1. Thermophysical data of Sn -5 wt.% Pb alloy [8].

811

List of table captions

Table 1

Simulation parameters for the 1D test case presented in Figure 4.

Table 2

Simulation parameters for the 2D test case presented in Figure 7.

Table 3

Simulation parameters for the 3D test case presented in Figure 11.

Table 4

Simulation parameters for the ingot case presented in Figure 14.

Table C1 Thermophysical data of Sn -5 wt.% Pb alloy [8].

Table D1 Time evolution of the distribution of the solid along with 1D simulation domain (dashed lines in Figure 5 and Figure 6).

 solution of pure transport equations by the finite element method,  sharp discontinuities in the velocity fields due to the packing of solid grains.

Propositions were introduced to overcome these difficulties:

 addition of an adaptive artificial diffusion to the transport equations,  implementation of a specific treatment to deal with the packing issue, consisting in using a transition function and adjusting the solid velocity.

A careful investigation was conducted to ensure the consistency between related quantities during the process as well as to guarantee the conservation of mass and energy. It progressively consisted of:

 a 1D pure transport simulation of sedimentation to evaluate its effects and to propose appropriate values for the adaptive artificial diffusion,  application of the complete transport-growth model to simulate macrosegregation in a 2D configuration [11],

 extension of the above simulation in 3D, resulting in a heterogeneous distribution of variables in the third direction which could not be captured by 2D simulations,  3D simulation of the solidification of a 3.

Appendix B -Solute Diffusion Lengths and Area Concentration

The solute diffusion lengths are taken from the work of Tveito and co-workers [41], as the following formulations.

Solute diffusion length in the liquid phase

where

𝑆ℎ 𝑐𝑜𝑛𝑣 = 2 3𝑔 𝑙 𝑆𝑐 1 3 ⁄ 𝑅𝑒 𝑛(𝑅𝑒) 𝑆𝑐 = 𝜇 𝑙 𝜌 𝑙 𝐷 𝑙 The area concentration is calculated as:

Appendix C -Thermophysical Data 810 Table C1. Thermophysical data of Sn -5 wt.% Pb alloy [8].

811

Appendix D -Analytical solution for the 1D Test Case

The 1D Test Case consists in pure sedimentation of a column of preexisting globular grains with fixed size in a uniform temperature domain. Considering constant and equal densities of the solid and liquid phases, as well as no phase change and no nucleation, the average total mass conservation simplifies to 𝑔 𝑠 〈𝒗 𝑠 〉 𝑠 + 𝑔 𝑙 〈𝒗 𝑙 〉 𝑙 = 0. For the sake of simplicity, a constant settling value of the solid velocity is imposed, set to 〈𝑣 𝑠 〉 0 𝒔 = -1 𝑚𝑚 𝑠 -1 . The 1D domain height and the initial conditions are defined in Figure 4: a continuous and uniform 60 mm mushy zone region is initially present between heights 20 mm and 80 mm, with a uniform average grain density per unit volume, 𝑁 0 = 10 9 𝑔𝑟𝑎𝑖𝑛𝑠 𝑚 3 , and volume fraction of solid, 𝑔 0 𝒔 = 0.1. One can easily derive the value for the liquid velocity in the mushy zone,

Similarly, the radius of the grains, 𝑅 0 , is simply given by using the definition of the fraction of solid, 𝑔 0 𝒔 = 𝑁 0 (4/3)𝜋𝑅 0 3 , leading to the value 𝑅 0 = 0.288 𝑚𝑚. Considering the fixed settling velocity and the packing limit at which the grain stop, 𝑔 𝑐 𝑠 = 0.3, the time evolution of the distribution of the mushy zone is simply derived by considering that the total fraction of the solid phase is unchanged over the entire domain, while not exceeding 𝑔 𝑐 𝑠 in the packed bed. Values are reported in Table D1.

The temperature is fixed to 498 K (224.856 °C), i.e. below the liquidus temperature of the 𝑆𝑛 -5 𝑤𝑡% 𝑃𝑏 alloy, that is 498.72 K (225.57 °C) according to the thermophysical properties listed in Table C1 of Appendix C [8]. The average solute mass composition is defined by 〈𝑤〉 = 𝑔 𝑠 〈𝑤 𝑠 〉 𝑠 + 𝑔 𝑙 〈𝑤 𝑙 〉 𝑙 . At any time, as the system is closed with respect to mass transfer, integration over the entire domain must retrieve the nominal composition of the alloy, 𝑤 0 = 5 wt% Pb. The initial composition profile assumes no macrosegregation. This means that the average composition is equal to 𝑤 0 at any position along the domain.

However, assuming complete mixing in both liquid and solid phases, the lever rule holds and one can derive the equilibrium intrinsic composition of the liquid and solid phases, 〈𝑤 𝑙 〉 𝑙 = 5.556 wt% Pb and 〈𝑤 𝑠 〉 𝑠 = 0.364 wt% Pb, respectively. Knowing the distribution of solid and liquid and their initial and intrinsic compositions, one can directly compute the average compositions by tracking the change of phases due to sedimentation. Computed values are reported in Table D1.

List of table captions

Table 1

Simulation parameters for the 1D test case presented in Figure 4.

Table 2

Simulation parameters for the 2D test case presented in Figure 7.

Table 3

Simulation parameters for the 3D test case presented in Figure 11.

Table 4

Simulation parameters for the ingot case presented in Figure 14.

Table C1 Thermophysical data of Sn -5 wt.% Pb alloy [8].

Table D1 Time evolution of the distribution of the solid along with 1D simulation domain (dashed lines in Figure 5 and Figure 6). given in Table 1. 0.3). The simulation case is defined in Figure 7 and Table 2. simulations. Black curves are isolines of solid fraction (0.1; 0.2; 0.3). The simulation case is defined in Figure 7 and Table 2. 3.

List of figure captions