Doctorat Paristech

Professeur, CEMEF Michel Michel Bellet

Hervé Combeau

Maître De

Charles-André Gandin

Professeur Jury M Andreas Ludwig

Professeur M Jesper Hattel

Professeur M Henri Nguyen-Thi

Mme Isabelle Poitrault

Professeur, Institut M Hervé Combeau

MINES Jean Lamour

Nancy Examinateur

M Charles-André

Gandin

Hervé

Prof. Jesper Andreas Ludwig

Prof Hattel

Ms. Isabelle Henri Nguyen-Thi

Ms Ms Carole Torrin

Ann-Mari Mulé

Mr. Mr Emmanuel Levrat

Mériem, Yamina, Ghazza Carole

Marie Bedel

Doan Kim

Frédéric Plourde

Multiscale finite element modeling of macrosegregation and grain transport ~~~ Modélisation multi-échelle d'éléments finis de la macroségrégation et du transport des grains

Keywords: solidification, modeling, macrosegregation, steels, multiscale, finite element, transport equations

Special thanks to all the people who I had a chance to meet, to

Research activities and modeling of solidification have rapidly developed since the past halfcentury, not only aiming at improving conventional solidification processes (e.g. ingot casting, continuous casting, rapid solidification), but also addressing advanced materials for which solidification is used as an effective and crucial manufacturing technique (e.g. metalmatrix composites, materials with unique microstructures, superconducting crystals).

Predicting the properties and structures of solidified products needs the knowledge related to many disciplines (such as heat, mass and solute transport, the formation kinetics of new phases, etc.) since a number of physical phenomena simultaneously take place and interact together during casting, from microscopic processes happening in a few seconds at the atomic level (nano scale), over interactions between different phases (between micron and centimeter scales), to the heat extraction and mass transport at the product scale (meter scale) during many hours.

In the context of casting industries, solidification process decisively influences on properties of final products. Mechanical properties of castings are determined by microstructures, e.g. fine-grained structures are favored for fatigue resistance, while coarse-grained structures show a better creep resistance in high temperature applications. The nucleation and evolution of solid structures at the microscopic scale lead to locally change the distribution of alloying elements during phase transformation. It results in an impact on the thermal field due to latent heat release as well as on the behavior of phase flow, e.g. through buoyancy forces and interactive phases drags. Heterogeneities in composition are then manifested at the system scale, that become more or less severe largely depending on the behavior of fluid flow and solid movement, the so-called macrosegregation, which is one of the main defects encountered in casting [Flemings, 1974]. The ability to predict macrosegregation in products obtained by solidification is essential for many industrial applications. This would provide economic benefits and time-saving in production. Macrosegregation in heavy industrial products poses challenges for numerical predictions because it results from various coupled micro-macro phenomena where transport velocity fields become significant and further considerable simulation cost is required.

A schematic of a solidifying ingot in Fig. 1. 1-a) illustrates several simultaneous phenomena occurring at an intermediate state, such as the liquid phase enriched in solute, local segregation channels, dendritic fragments, settling equiaxed grains [Ludwig et al., 2015a]. Fig. 1. 1-b) presents the (left side) macrostructure and (right side) carbon macrosegregation of a steel ingot weighted 65t and containing 0.22wt% nominal carbon. It consists of different structural types and heterogeneous solute distribution: including columnar structure near the surface with a homogenous solute composition; globular equiaxed grains in the center from the bottom up to 75% of the total height, corresponding to a negative segregation zone; and dendritic equiaxed grains on the top of the ingot with a positive segregation [Lesout, 2005].

a) b)

Fig. 1. 1 -Schematic of (a) a solidifying ingot at an intermediate state, issued from [Ludwig et al., 2015a], (b-left) macrostructure and (b-right) carbon macrosegregation of a 65t steel ingot containing 0.22 wt% nominal carbon, issued from [Lesout, 2005].

With an increasing computational power, solidification modeling to predict macrosegregation has been enhanced by properly addressing more and more phenomena at both microscopic and macroscopic levels with respect to their complex dependence. Starting from a simple model that simulates solidification by only solving the heat conservation equation, then other physical issues have been subsequently dealt with, such as melt convection, transport of chemical species, diffusion at different scales, etc. Besides additionally integrating physical aspects, various numerical methods have been also developed to be able to study processes from a very small scale until a practical level of massive castings, in consideration of establishing effective mathematical solutions for a system of non-linear equations used to describe multi-scale multi-physical phenomena.

Reviewing applications to industrial ingots, in 1999, Gu and Beckermann applied the authors' model, a multi-scale multi-component solidification model [Scheneider and Beckermann, 1995], to predict 2D macrosegregation for a heavy steel ingot [Gu and Beckermann, 1999], for which the solid phase was assumed to be fixed. c) display the profiles of macrosegregation in carbon and in sulfur, respectively, along the center of the ingot at the end of solidification. According to the authors, although a similar tendency was generally achieved between numerical solutions and experimental measures, the negative segregation in the lower part of the ingot could not be predicted because of neglecting the sedimentation of solid grains. In addition, the positive segregation in the upper part was underestimated for the case of carbon since the absorption of carbon from insulation materials at the top of the ingot was not included. With a great investigation, recently subsequent simulations applied to large castings demonstrated and confirmed that the negative segregation at the lower zone of ingots is mainly due to sedimentation of equiaxed grains, such as applications to a 3.3t steel ingot [Combeau et al., 2009][Založnik andCombeau, 2009a], a 6.2t steel ingot [Kumar et al., 2012], a 65 steel ingot [Combeau et al., 2012a], a 2.45t steel ingot [Wu et al., 2014] [Ludwig et al., 2015a]. Fig. 1. 3-a) and -b) present the relative segregation profiles along the center line and the relative segregation maps of the 65t steel ingot, respectively. The results are issued from the work of Combeau and co-workers [Combeau et al., 2012a]. Numerical predictions obtained for two case studies were compared to (green curve) experimental measurement, including (blue curve) the case where solid grains are imposed to be fixed and (red curve) the case where the solid motion is accounted for. The latter case reproduced the typical segregation pattern of the ingot (negative segregation at the lower zone and positive segregation at the upper zone). It can be seen in the corresponding experimental characterization (Fig. 1. 1-b)) that the globular grains occupy a large region, compared to other present structures (columnar and dendritic equiaxed structures). In addition, for such ingots of large dimensions, the movements of the solid and liquid phases were shown to be predominantly important. Therefore, accounting for the motion of globular grains due to both liquid convection and sedimentation is essential to predict the macrosegregation. Although significant achievements of solidification modeling have been obtained and numerical simulations can predict essential features of experimental measurements, more investigation would be expected to accurately anticipate properties of final castings in order to control defects and to produce alloys with desired qualities. This is also the domain which is investigated by our present work, focusing on the effect of solid motion. As strong coupling of the solid motion with other phenomena leads to solve a complex set of non-linear equations, it has been limitedly studied, though increasingly, for a detailed understanding about its influences. However, in large industrial casting products such as ingots, the solid motion yet becomes a dominant factor which affects the distribution of structures and solute composition. Moreover, studies were mostly carried out with two-dimensional simulations.

Integrated into a research and development casting project financed by a metallurgical industrial group consisting of ArcelorMittal-Maizières, ArcelorMittal-Industeel, Aperam, Aubert & Duval, and Asco Industries, our work will focus on the prediction of macrosegregation due to the movement of equiaxed solid grains, which is developed to be implemented into a three-dimensional casting software, Thercast® (*), used by the industrial partners. Besides ongoing advancements achieved for this issue from using finite volume methods, our work aims at contributing to this aspect an alternative numerical tool which utilizes a finite element method.

(*) Thercast® is a 3-D finite-element software package that simulates ingot casting and continuous casting of steel and other metallic alloys. It is developed jointly by CEMEF and the company Transvalor. (http://www.transvalor.com) Chapter 1 -Introduction

Before considering a literature study of numerical solidification models related to this topic, let us review very briefly basic physical phenomena inducing macrosegregation in alloy solidification, in order to concisely show how the heterogeneous distribution of chemical species is formed at the scale of dendrite arms and then is exhibited at a larger scale.

1.2 Macrosegregation in Alloy Solidification 1.2.1 Segregation at Microscopic Scale A typical equilibrium phase diagram for a simple binary alloy A-B at constant pressure is shown in Fig. 1. 4. Nucleation of solid crystals occurs at a temperature slightly below the liquidus temperature of an alloy at a given composition. Due to the partition of chemical species during the phase transformation, in most cases (with a solute partition coefficient less than unity), the solid phase is formed at a smaller composition than the average value. The rejection of solute from the solid to liquid leads to an increase of the concentration in the liquid phase. There exists a mushy zone consisting simultaneously of the solid and liquid phases over a temperature range.

Under equilibrium and well-mixed conditions that are often assumed for simplicity, the evolution of concentrations in the solid and liquid phases follows the solidus and liquidus lines, respectively, up to their maximum solubility at the eutectic temperature. However, accounting for limited solute diffusion was demonstrated to be important in order to correctly model the formation and development of solidification structures. In fact, mass diffusivity being smaller than thermal diffusivity in metal alloys, solidification at the microscopic scale is principally governed by solute concentration gradients on each side of the solid-liquid interface. Due to a limited solute diffusion in the solid and liquid phases, rejected species from the solid phase are accumulated at the solid-liquid interface and then diffused into the liquid phase. So there exists gradients of composition in the solid and liquid phases as schematized in Fig. 1. 5. Chapter 1 -Introduction

The concentration difference in the liquid phase, ( ), over the one at the interface, (

), is referred to the solutal undercooling. It becomes dominant compared to the thermal counterpart defined by the temperature variation in the context of alloy solidification. Therefore, it governs the nucleation and evolution of solid structures as well as the solute distribution at the scale of dendrite arms. In turn, the morphology of microstructures has a crucial impact on the solute diffusion flux. Solidification microstructures can be divided into columnar (constrained-growth) and equiaxed (unconstrained-growth) structures. In columnar solidification, most of the solute is rejected principally in the lateral direction perpendicular to the heat flow. While in the case of equiaxed crystals that are surrounded by undercooling melt, rejection of solute takes place in all directions on their solid-liquid interface and eventually for an accumulation layer outside the grains. For each structural type, solid grains can grow with more or less dendritic shape, corresponding to different growth and solidification kinetics at the soli-liquid interface. The formed structures and the grain morphology depend on many factors such as phase change conditions, composition, thermodynamic properties of phase transformations, impacts of phase convection, etc. Various subjects related to the nucleation and growth of solidification structures can be consulted in [Kurz and Fisher, 1989] [ Dantzig and Rappaz, 2009].

Advanced solidification models are able to account for the effect of convection on solute gradients in a small layer near the solid-liquid interface as well as for complex shapes of interfacial structures that are associated to instability conditions and curvatures of interfaces.

Studies of these topics can be found in [Flemings, 1974] [ Kurz and Fisher, 1989] [Langer, 1989] [ Pines et al., 1990] [ Wang andBeckermann, 1993][Martorano et al., 2003]. Chapter 1 -Introduction

Segregation at Macroscopic Scale

Due to the movement of the liquid and solid phases during the solidification process, the segregation of chemical components at the scale of solid-liquid interface is manifested on a system scale, the so-called macrosegregation.

Generally, melt flow can be generated by many sources such as external forces, surface tension gradients, residual flow due to filling of the mold, buoyancy forces due to temperature and compositional gradients, drag forces from solid motion, etc. In the present work, the liquid movement is governed by the two later listed mechanisms. The thermal and solutal buoyancy forces can oppose or corporate to each other, depending on heat exchange conditions and on the weight of chemical species rejected in comparison with that of the bulk liquid. Additionally, during solidification these effects can be enhanced or weakened by the solid motion, according to the orientation of these forces relative to the gravity direction.

In consideration of the solid phase, its movement can take place within a small distance, for instance deformation-induced displacements in continuous casting processes. The solid motion also manifests over large domains such as the settling and floating of equiaxed grains in ingot castings. Relating to the later mechanism, its influence on macrosegregation was experimentally observed for different solidified alloys, i.e. undercooled Pb-Sn eutectic alloys [de Groh [START_REF] De Groh | [END_REF] and aluminum alloys (Al-1wt%Cu, Al-10wt%Cu with different amounts of grain refiner) [START_REF] Rerko | [END_REF]. Settling of equiaxed crystals is also known as the principal factor causing the negative segregation cone at the bottom casting in Fe-C alloys solidification [Flemings, 1974]. The origin of the equiaxed crystals has not been entirely clear. Two principal formation mechanisms were proposed including: heterogenous nucleation and detachment and transport of dendrite arms initially formed in a columnar zone.

It can be seen that macrosegregation involves physically different multi-scale phenomena. Consequently, in order to predict macrosegregation, it is necessary to simultaneously model the microscopic growth processes and the macroscopic happenings such as heat transfer, mass transport, phases advection… From the numerical point of view, accounting for the solid motion requires dealing with a complex mathematical system because of its close interaction with other phenomena. Chapter 1 -Introduction

Objectives and Outline

In the aim of modeling macrosegregation with ultimate applications to industrial castings of large sizes or under complex geometries, developing FE solidification model would expand the capacity of numerical predictions, especially being promising for such industrial applications.

With interest in successes achieved from the FV solidification model developed by Založnik and Combeau [Založnik andCombeau, 2010a, 2010b] in simulating macrosegregation, in particular for heavy ingots as presented above, our study's objective is to adapt and implement this model into the FE framework on which the 3D simulating program Thercast® is developed, in order that it is able to model the equiaxed transport during solidification.

Following this chapter in which main physical phenomena related to macrosegregation have been introduced,  a review of solidification models will be presented in Chapter 2, consisting of the multi-scale modeling, volume-averaged method, governing equations used to describe phenomena during solidification and volume-averaged solidification models accounting for solid motion.

 Finite element implementation and adaptation will be detailed in Chapter 3, in which the resolutions for macro-micro equations are established by using the splitting scheme.

 Chapter 4 will present two-dimensional numerical results obtained from the current model that will be compared to references and analyzed for the purpose of verifying and validating our work, taking in first for purely growth process, then with melt convection, then with only transport phenomena and finally solidification cases using the complete growth-transport model.

 Three-dimensional simulations and industrial applications are the subject of Chapter 5.

 Conclusions and perspectives are presented at the end of the report. Chapter 1 -Introduction

Résumé en français

Afin de modéliser la macroségrégation pour des produits industriels de grandes tailles ou de géométries complexes, un modèle de solidification par éléments finis sera développé. Une modélisation par éléments finis peut permettre d'étendre de façon prometteuse la capacité de prédictions numériques du phénomène de macroségrégation dans les applications industrielles.

Au vue de l'intérêt du modèle de volumes finis développé par Založnik et Combeau [Založnik andCombeau, 2010a, 2010b] pour la simulation de la macroségrégation, en particulier pour de larges lingots d'acier, notre étude a pour but d'adapter et de mettre en oeuvre ce modèle dans une formulation de type éléments finis utilisée dans le logiciel industriel Thercast®, afin qu'il soit capable de modéliser la solidification en présence du transport des grains solides.

Suite à ce chapitre où les principaux phénomènes physiques liés à la macroségrégation ont été introduits,

 une revue bibliographique des modèles de solidification sera présentée dans le chapitre 2, composée de quatre parties : Modélisation multi-échelle ; Méthode des prises de moyenne ; Ensemble d'équations gouvernantes ; Modèles de solidification utilisant la technique de prise de moyenne avec la présence du transport des cristaux.

 La mise en oeuvre et l'adaptation numérique seront détaillées dans le chapitre 3, dans lequel la résolution des équations macro-micro est établie en utilisant le schéma de splitting.

 Le chapitre 4 est consacré à présenter et à analyser des résultats numériques obtenus par le modèle actuel dans une configuration bidimensionnelle, en les comparant aux solutions de référence.

 Des simulations tridimensionnelles et des applications industrielles sont menées au chapitre 5.

 Les conclusions et les perspectives seront portées à la fin du manuscrit. Chapter 2 -Literature Reviews Chapter 2

Literature Reviews

2.1 Multi-scale Modeling

In general, continuum theories rely on conservation laws which are expressed in forms of differential equations. Besides long-range macroscopic heterogeneities in the continuum, heterogeneities principally arise from subscale microstructures (e.g. particles, defects, inclusions, etc.) that decisively influence the behavior and performance of systems.

Multi-scale modeling developed aims at formulating a description of material properties or system behavior at the macroscopic level using information from other sub-scale levels that are usually distinguished as following types: molecular level, microscopic level and macroscopic level. Each level addresses phenomena over a specific scale of length and time.

Being developed in the 20 th century, multi-scale modeling methodologies have advanced quickly due to the combination of parallel computing power, experimental capabilities to characterize structure-property relations down to the atomic level and theories that admit multiple length scales. Multi-scale modeling provides a useful tool in predicting macroscopic behavior regarding important smaller-level features. It has been applied and investigated for various disciplines (solid mechanics, fluid mechanics, materials science, physics, mathematics, biological, chemistry, etc.).

In solidification, relevant phenomena occur over a wide range of length and time scales, which we have looked over in Chapter 1. It can be recalled that for the mass distribution, while the local solute segregation due to the growth of solid phase takes place quickly at a solid-liquid interface approximated by some ten microns (the scale of dendrite arms), the liquid motion leads to the transport of mass and heat through the whole domain sized in several meters and also affects the local solidification conditions. During this process, the melt flow interacts with the solid grains on the scale of some hundred microns, e.g. the flow may be dragged to go down or may be pushed up because of the settling of solid grains, and it can be damped inside porous areas. An illustration in Fig. 2. 1 presents ingot castings with several meters size for which the heterogeneity of solute which is manifested at the product scale results from a partition of solute at the scale of dendrite arms during solidification. Analyses of dendrite structures at a smaller scale such as a molecular or atomic level would provide serviceable parameters for numerical modeling as well as give a better understanding of physical behaviors at substantial scales for a global casting. In this context, application of multi-scale modeling to solidification has received a great attention as it allows conveniently linking between microscopic evolutions that can proceed within a few seconds and macroscopic transport phenomena evolving during a few hours. Chapter 2 -Literature Reviews Basic physical phenomena related to the macrosegregation have been reported in Chapter 1.

In this section, we will investigate the involved mathematical and numerical aspects. First let us review the development of multi-scale modeling applied to the alloy solidification process.

Fig. 2. 1 -An illustrative figure presenting casting products (of several meters size) for which heterogeneities in compositions and structures are originated from microscopic processes that are then manifested at the system scale through advective mechanism Two fundamental approaches are proposed. In the first one, macroscopic continuum equations are obtained by applying a mixture theory, which are expressed by mixture quantities [Bennon andIncropera, 1987a, 1987b]. In the second one, continuum equations are derived from local microscopic equations by averaging techniques [Drew, 1983] [Ni andBeckermann, 1991]. In these methods, conservation equations are solved on a single computational domain, thus this precludes tracking phase fronts and discretization on irregular zones that are required when using multi-domain procedures. Between these two methods, while it is difficult to account for microscopic processes in models based on the mixture theory, those using the averaging technique provide a physically natural integration of microscopic processes into macroscopic equations.

Derived from general theories of flow through porous media and of multi-phase systems by [Batchelor, 1970] [Gray, 1975, 1983] [Hassanizadeh and Gray, 1979] and [Drew, 1983], volume averaging technique was first applied to solidification of binary alloys by Beckermann and co-workers, without solid motion [Beckermann and Viskanta, 1988] and then a general model accounting for liquid and solid movements [Ni and Beckermann, 1991], as well as by Ganesan and co-workers [Ganesan andPoirier, 1990] [Poirier et al., 1991]. Simulation applications were performed for equiaxed globular solidification without and with solid movement [Feller andBeckermann, 1993] [Ni andBeckermann, 1993]. After that, the models were extended to solidification of multi-component alloys in considering the stationary solid phase [Schneider and Beckermann, 1995] [Beckermann and Schneider, 1995] [Schneider et al., 1997] [Gu and Beckermann, 1999]. The equiaxed dendritic solidification focusing on microstructure descriptions involving the germination and growth processes controlled by solute diffusion was first investigated by Rappaz and co-workers [Rappaz andThévoz, 1987] [Thévoz et al., 1989]. Later on, with an introduction of a second liquid phase into the volume-averaged two-phase models, the transport of dendritic grains was modeled by considering two distinct liquid phases at both scale in the work of Wang and Beckermann [Wang and Beckermann, 1993, 1996a,1996b,1996c] or by assuming the same macroscopic behavior for two liquid phases as proposed in the work of Combeau and co-workers [Combeau et al., 2009]. Various advanced developments and numerical applications have been subsequently contributed considering different phenomena in the solidification processing in order to improve predictions of casting products' properties, such as accounting for multi-scale solute diffusion at the solid front [Wang and Beckermann, 1993][Martorano et al., 2003] [Bedel, 2014], tracking the grain structure using the CAFE method [START_REF] Gandin | Modelling of Grain Structure Formation during Solidification[END_REF]] [Carozzani et al., 2012], modeling both columnar and equiaxed (globular and dendritic) structures and their transition [Wu and Ludwig, 2009, 2010] [Leriche et al., 2015], taking into account the solid deformation due to external forces as well as interphase interactions [START_REF] Bellet | [END_REF]] [Bellet et al., 2004] [Ludwig et al., 2015b]. More reviews can be read in [Beckermann, 2002] [Pickering, 2013].

As our study is based on the volume averaging technique, we will summarize in the following section its principle concepts and relevant final formulations, further details can be found in the above references.

Volume-Averaged Method

The volume averaging technique is used to establish macroscopic continuum equations on a representative elementary volume (REV) that is much smaller than the system size and large enough to be characteristic of related structures. The REV consists of different phases that may possess various physical properties. These phases are separated from each other by interface boundaries. According to [Ni and Beckermann, 1991], under typical solidification conditions, i.e. around several meters for the size of systems and in several micrometers for characteristic lengths of microstructures, a REV can be found from 10 mm to 1 mm, therefore the REV should contain several crystals.

The average of a microscopic quantity in phase 𝛼 over the REV is defined as:

〈 〉 ∫ (1)
where is the volume of a REV and the phase function is defined as:

( ) { 𝛼 (2)
It should be noticed that the volume average of the phase function is equal to the phase volume fraction, e.g 𝑔 for phase 𝛼. The intrinsic volume average is expressed as the average with respect to the proper volume of phase 𝛼, :

〈 〉 ∫ 〈 〉 𝑔 (3)
The corresponding fluctuation of quantity is calculated as:

̂ ( 〈 〉 ) (4)
The average of the product of two quantities and is expressed as:

〈 〉 〈 〉 〈 〉 〈 ̂ ̂ 〉 (5)
The average of temporal and spatial derivations can be deduced as follows [Ni andBeckermann, 1991] [Rappaz et al., 2003]:

〈 〉 〈 〉 ∫ (6) 〈 〉 〈 〉 ∫ (7)
where is the interfacial area enveloping phase 𝛼; is the interface velocity; is the unit normal vector, directing outward from phase 𝛼.

Conservation Equations

The averaging procedure applied to the equations studied in the present work will be summarized in this section. More details can be consulted in [Ni and Beckermann, 1991], [Beckermann and Viskanta, 1993], [Rappaz et al., 2003], [Založnik and Combeau, 2010a].

For each phase of an alloy in its mushy state (including solid and liquid phases), the local continuum, momentum, energy, and species conservation equations, respectively, are:

( ) (8) ( ) ( ) (9) ( ) ( ) (10) ( ) ( ) (11)
where is the density; the velocity; the pressure; the deviatoric stress tensor; the body force; the specific enthalpy; the heat flux vector; the species concentration; the species flux vector. Note that no solute index is introduced in Eq.( 11), it should be understood that each solute should give rise to the solution of a conservation equation.

At the scale of a REV of the mushy material, by applying the volume-averaged approach to the above local continuum equations and neglecting fluctuation terms, the corresponding volume-averaged conservation equations obtained are as follows:

(𝑔 〈 〉 ) (𝑔 〈 〉 〈 〉 ) (12) (𝑔 〈 〉 〈 〉 ) (𝑔 〈 〉 〈 〉 〈 〉 ) 𝑔 〈 〉 (𝑔 〈 〉 ) 𝑔 〈 〉 (13)
(𝑔 〈 〉 〈 〉 ) (𝑔 〈 〉 〈 〉 〈 〉 ) (𝑔 〈 〉 ) (14)
(𝑔 〈 〉 〈 〉 ) (𝑔 〈 〉 〈 〉 〈 〉 ) (𝑔 〈 〉 ) ( 15
)
where 𝑔 is the volume fraction of phase 𝛼; the mass flux of phase 𝛼 across the phase interface due to phase change; the phase mass transfer rate due to grain nucleation; and the interfacial momentum transfers due to phase change, interfacial stress and grain nucleation respectively; and the interfacial heat transfers due to phase change, diffusion and grain nucleation respectively; and the interfacial solute transfers due to phase change, diffusion and grain nucleation respectively.

In Eqs.( 12)-( 15), the time variation and the change due to macroscopic transport are presented on the LHS. Exchanges between a phase 𝛼 and its connected phases arise from microscopic processes, including the interface movement, diffusion across the interface and nucleation, which are arranged on the RHS.

For those associated with the interface movement, the related fluxes in and out of phase 𝛼 which are defined as the product of the surface integral of quantity at the interface and the velocity of this interface relative to the phase, are given as ( 16), ( 17) and ( 18) for the momentum, solute mass and energy respectively.

∫ ( ) (16)
∫ ( ) (17)
∫ ( ) (18)
where is the density of phase 𝛼; is the velocity of phase 𝛼; is the velocity of the interface; is the unit normal vector at the interface, directing outward from phase 𝛼; is the solute concentration in phase 𝛼;

is the specific enthalpy of phase 𝛼; is the interface area; is the averaging volume.

These integrals can be modeled by the product of the averaged value of at the interface and the phase change rate, which are expressed as ( 19), ( 20) and ( 21) for the fluxes of the momentum, solute mass and energy respectively.

̅ (19) ̅ (20) ̅ (21)
where ̅ ; ̅ ; ̅ are the average values over the interfacial area , in , of velocity, solute concentration and energy, respectively, of phase 𝛼.

Concerning the transfers induced by interfacial stresses or diffusion, they are introduced as ( 22), ( 23) and ( 24) for the quantities of momentum, solute mass and energy respectively.

∫( )

∫ (23) Chapter 2 -Literature Reviews ∫ (24) (22)
where is the deviatoric stress tensor in phase 𝛼; is the pressure in phase 𝛼; is the unit tensor; is the species flux vector in phase 𝛼;

is the heat flux vector in phase 𝛼; is the unit normal vector at the interface, directing outward from phase 𝛼.

These fluxes are often modeled as functions of the difference between the interfacial average value of a quantity and its intrinsic volume average one, leading to equations ( 25), ( 26) and ( 27) for momentum, solute mass and heat transfer respectively.

( ̅ 〈 〉 ) (25) ( ̅ 〈 〉 ) (26) ( ̅ 〈 〉 ) (27)
where is the density of phase 𝛼; is the momentum resistance coefficient of phase 𝛼; is the solute diffusion coefficient in phase 𝛼;

is the characteristic solute diffusion length in phase 𝛼;

is the heat conductivity of phase 𝛼; is the characteristic heat conduction length in phase 𝛼; ⁄ [m -1 ] is the interfacial area concentration (identical for the two phases present in the REV considered).

The contribution of nucleation is evaluated by the corresponding quantities generated or lost due to the formation of a new phase and it can be neglected in comparison with other terms. Chapter 2 -Literature Reviews

Volume-Averaged Solidification Models accounting for Solid Transport

From a simple solidification model that only considers heat diffusion, neglecting melt convection, solid movement, and species redistribution on a macroscopic scale, modeling of alloy solidification accounting for liquid convection began with the work of Flemings and coworkers [START_REF] Flemings | [END_REF][Flemings et al., , 1968a[Flemings et al., , 1968b]]. Inclusion of the solid transport by using the volume averaging technique was first developed by Ni andBeckermann [Ni andBeckermann, 1990, 1991]. Compared to the mixture approach, this formulation has the advantage that microscopic processes are connected in a clear and persistent manner to macroscopic phenomena because the volume-averaged macroscopic equations are directly derived from the local conservation equations.

In the following, we will overview three solidification models taking into account the transport of solid particles and based on the averaging technique as well as the models previously developed at CEMEF. Through this review, we intend to provide a clear overview of the domain's development and to situate the beginning for our present work.

Moreover, before entering into these volume-averaged models, it would be worthy to look over achievements related to the solidification modeling with the presence of solid movement when using the mixture theory. Roberts and Loper were the first ones who developed mathematical models of alloy solidification which take into account the solid motion [START_REF] Roberts | Dynamic Processes in Slurries. Structure and Dynamics of Partially Solidified System[END_REF]. A subsequent mixture model was developed by Vreemann and co-workers, coupling the transport of solid crystals and fluid flow to study the solute redistribution during solidification [Vreemann et al., 2000a]. Numerical applications were performed on direct chill castings of binary alloys [Vreemann et al., 2000b[START_REF] Vreemann | Direct Chill Casting of Aluminum Alloys: Modeling and Eperiments on Industrial Scale Ingots[END_REF], and recently extended for ternary alloys [Krane, 2004]. The authors showed that accounting for the particle motion is important to correctly evaluate the solute distribution, especially for prediction of the negative segregation zone at the ingots' center, which is observed in experimental measurements. Several results from these works were introduced in Chapter 1.

Models of Beckermann and co-workers

In their article published in 1991 [Ni and Beckermann, 1991], the authors presented a twophase model using the averaging method for a general case of the solid and liquid moving at different velocities, in which separate conservation equations are developed for each phase, as presented in Section 2.3, together with the corresponding interfacial balances. While most models assume that the liquid is in equilibrium and locally well mixed, e.g. the lever rule or Scheil models, the non-equilibrium microscopic phenomena, such as nucleation and undercooling which are crucial processes in equiaxed solidification in reality, are accounted for and consistently introduced into the macroscopic equations in this model.

In consideration of solid crystals, their distribution results from the advection and nucleation mechanisms, which is described by a population balance equation ( 28).

( ) ̇ (28)
where is the grain density; is the transport velocity and ̇ is the source term which represents the net generation of crystals and can further account for phenomena such as crystal "birth" (nucleation, dendrite arm breakoff) and crystal "death" (full melting, agglomeration of crystals).

In order to numerically solve a set of aforementioned equations, following the authors' approach, the phase change rate, (in Eq.( 12) in Section 2.3), is first evaluated by using the last known quantities, from which other microscopic interfacial exchange terms are deduced. These quantities are then injected into the conservation equations coupling microscopic processes and macroscopic transport phenomena. After that, a finite volume method, using an implicit upwind scheme, is applied to solve the momentum, solute mass, energy, grain density equations. An iterative procedure is effectuated by re-evaluating the phase change rate and other interfacial quantities, then resolving the micro-macro equations until all fields are converged. Details of the numerical implementation are presented in [Wang andBeckermann, 1994, 1996b].

The two-phase model was then simplified and applied, in a finite volume formulation, to solidification of globular grains for which nucleation takes place with a constant density at the liquidus temperature [Ni andBeckermann, 1993] [Beckermann andNi, 1996]. Numerical simulations of solidification were performed on a 0.05 m x 0.1 m rectangular cavity containing an Al-4wt%Cu alloy, being cooled on the left vertical wall, other walls are adiabatic. Three test cases were carried out, including (Case 1) solidification without solid transport, (Case 2) solidification in the presence of solid transport and with the nucleation rate set up at grains m -3 , (Case 3) also in the presence of solid transport but with a smaller nucleation rate of grains m -3 . Through results obtained, considerable effects caused by the motion and nucleation of solid grains are demonstrated on the evolution of quantities and their interactions during solidification and thus on the formation of macrosegregation and the heterogeneous distribution of grain sizes. As illustrated in Fig. 2. 2, significant differences of the liquid velocity and solid volume fraction, displayed at an intermediate time (30 s), are found from the three calculations. Detailed analyses and results at the end of solidification are referred to the above references. issued from [Beckermann and Ni, 1996]. Liquid velocity vectors and solid volume fraction contours at 30 s for three simulations:

Case 1 without solid motion.

Case 2 with solid motion and nucleation rate at grains m -3 . Case 3 with solid motion and nucleation rate at grains m -3 .

The authors then extended the model to equiaxed dendritic grains [Wang and Beckermann, 1996a], in which three distinguished phases are considered including the solid dendrite, the interdendritic liquid and the extradendritic liquid as schematized in Fig. 2. 3. Interactions between the solid and the total liquid are generally modeled on the entire range of the solid fraction for a single particle as well as for multi particles, as proposed in their development [Wang et al., 1995]. A partitioning of the liquid phase into two separate phases is considered at both micro and macro levels. With assumed simplifications, numerical simulations were performed on an Al-4wt%Cu alloy and on a NH 4 Cl-70wt%H 2 O, reported in [Wang andBeckermann, 1996b, 1996c]. The effect of solid transport can be seen by comparing simulations with and without solid motion, computed on an identical domain containing an Al-4wt%Cu alloy like the previous study. Fig. 2. 4 shows the final macrosegregation maps for three test studies including Case 1 without solid motion, Case 2 with solid motion and the nucleation rate at grains m -3 favoring the formation of dendritic grains, and Case 3 with solid motion and a higher nucleation rate at grains m -3 indicating the case of globular grains. It can be observed that macrosegeration is less pronounced when accounting for solid movement and even becomes less severe with increasing the nucleation rate. Since the solid and liquid motions have the same direction, the macrosegregation due to the relative velocity is thus reduced and more weakened by decreasing the grain size. In other systems where the phase movements are opposite, the segregate solute phenomenon would be amplified. Besides demonstrating the effect of solid motion, the results given by these dendritic studies, compared to those from the model of globular grains, allowed investigating in more detail the formation mechanisms of heterogeneities and complicated impacts related to the floating and sedimentation of dendrites during solidification. Though such a distinct consideration for the liquid phase may provide more realistic descriptions of phase interactions that could be important in certain contexts, it requires multiple supplementary parameters that must be identified through experimental measurements. v dinterdendritic liquid velocity ; v lextradendritic liquid velocity. Figures issued from [Wang and Beckermann, 1996a].

Case 1

Case 2 Case 3

Fig. 2. 4 -Final macrosegregation maps obtained from the three-phase model Results issued from [Beckermann and Ni, 1996b]. Case 1 without solid motion. Case 2 with solid motion and nucleation rate at grains m -3 . Case 3 with solid motion and nucleation rate at grains m -3 . Chapter 2 -Literature Reviews

Models of Wu and co-workers

A simplified version of the previous two-phase volume-averaged model, using a heterogeneous nucleation, was developed only for globular grains by Wu and co-workers [Ludwig andWu, 2002] [Wu et al., 2003] in order to study equiaxed solidification while avoiding uncertainties due to the lack of realistic information required in the complex model proposed by Beckerman and co-workers. A similar resolution approach as the precedent one was proposed to solve the above micro-macro equations, microscopic interfacial exchange terms are first computed in subroutines using the quantities obtained from the latest iteration, and then the full conservation equations are solved and coupled together by a fixed point procedure. Details about the numerical resolution strategy implemented in finite volumes are presented in [START_REF] Ludwig | [END_REF].

In order to study the solid transport, the authors' two test cases (Case 1) without solid motion and (Case 2) with solid motion are introduced here, which are issued from [Wu et al., 2003].

The simulations were performed on a 2D casting cavity, containing an Al-4wt%Cu alloy, being cooled from the bottom and side walls, as schematized in Fig. 2. 5-a), for which the thermo-solutal liquid convection was neglected, only sedimentation-induced flow was taken into account. From the authors' results, the effect of solid movement is visualized and analyzed. The intermediate evolution of variables, at t = 40 s, is displayed in Fig. 2. 5-b) andc), presenting the solid fraction contours and the liquid velocity vectors in Fig. 2. 5-b) and the average composition in Fig. 2. 5-c), those on the LHS for Case 1 and those on the RHS for Case 2. When accounting for the solid motion, liquid flow is enhanced by grain sinking and the flow path is changed with the appearance of vortices, leading to a different solute distribution compared to the fixed solid case. Moreover, the solidification sequence proceeds faster in the bottom than side walls due to solid sedimentation.

After that, the model was extended to include the simultaneous development of columnar and equiaxed structures, assuming simplified morphologies (a globular form attributed to equiaxed grains and a cellular shape to columnar trunks) [Ludwig andWu, 2005] [Wu andLudwig, 2006]. The dendritic morphology was then considered for equiaxed grains [Wu and Ludwig, 2009] and for both equixed and columnar structures [Wu et al., 2010a[Wu et al., , 2010b]].

Recently, the authors investigated their model for ternary alloys [Wu et al., 2013[Wu et al., , 2014]]. Numerical studies performed on small ingots and an application to a large steel ingot considering in a 2D-axissymmetric configuration, both of them verified that the cone-shape negative segregation which is found in solidified ingots results from grain sedimentation, one of results being shown in Fig. 2. 6. Additionally, the authors noticed that the interaction between settling equiaxed grains and either melt flow or columnar tip front can have an important impact to the formation of A-type segregates in ingots. For the simulation applied to a large ingot, the authors restricted their model to a simple consideration of the solute diffusion length and grain morphology, expecting better predictions when these limitations will be overcome. The models implemented in the SOLID software are described in detail by the authors in [Založnik and Combeau, 2010a] and [Combeau et al., 2009] for the two-phase globular and three-phase dendritic models, respectively. The principles of the two models will be summarized in this section.

In these models, the conservation equations, like those presented in Section 2.3, are derived from local continuum equations for each phase using the volume averaging technique and the numerical implementation is carried out based on a finite volume method. The authors use an operator-splitting scheme to split and combine the microscopic and macroscopic phenomena of solidification. According to this scheme, the variation of phases and solute concentrations during solidification is considered to be due to macroscopic transport phenomena and microscopic processes.

With this consideration, the two-stage resolution was developed.

o In the first stage, the evolvement associated with the transport mechanism is calculated by globally solving the averaged equations of phases and solute mass, solely accounting for the transport terms and using a finite volume method with a macro time-step. This is called the transport stage.

o In the second stage, the transport part is no longer taken into consideration, the grain growth is evaluated through the mass balance of interfacial exchanges at the scale of dendrite arms, assuming thermodynamic equilibrium at the solid-liquid interface, locally solving in each control volume with a micro time-step. It is corresponding to the so-called growth stage.

The transport velocity fields are used in the first stage to be computed by solving, outside the splitting scheme, the full average momentum conservation equations or a simplified equation for the solid phase by neglecting the inertial and viscous terms. Two flow regimes are considered according to the solid fraction. If the solid fraction is smaller than the packing value, equiaxed grains are free to move with respect to the liquid. The solid phase is considered as a continuous medium and characterized by its volume fraction and a number of particles contained (or an averaged grain radius). The solid-liquid interaction follows the model of Agarwal and O'Neill [Agarwal and O'Neill, 1988], the transferred momentum due to interfacial stresses is presented by ( 29)-(34). In the other case, i.e. the solid fraction is higher than the packing limit, solid grains are considered to be blocked and form a porous solid bed where the liquid movement is drastically slowed down and modeled by the Carman-Kozeny law for the permeability, as given by ( 35)-( 36).

(*) SOLID is a two-dimentional numerical simulation software dedicated to the casting process. It is used for continuous casting as well for ingot casting and performs an analysis of thermo-metallurgic behavior of the material during the process.

𝑔 ‖〈 〉 〈 〉 ‖(〈 〉 〈 〉 ) (29)
where 𝑔 (30)

{ 𝑔 (𝑔 ) ( 𝑔 ) ( 𝑔 ) ( 𝑔 ) 𝑔 (31) { 𝑔 ( ) [ (𝑔 ) ⁄ ( 𝑔 ) ⁄ ] 𝑔 (32) ( 𝑔 ) (33)
𝑔 ‖〈 〉 〈 〉 ‖ (34)
Packed bed regime :

(𝑔 ) (〈 〉 〈 〉 ) (35)
where

(𝑔 ) ( 𝑔 ) (36)
in which 𝑔 ; 𝑔 are respectively the solid and the liquid fraction; is the density; is the grain diameter; is the Reynolds number; is the dynamic viscosity; is the secondary arm spacing.

Modeling of nucleation can be integrated into the growth stage, together with the modeling of grain growth. New grains are generated when local conditions are favorable, following an instantaneous uniform model with a predefined constant number of nuclei and a predefined constant initial diameter or following an instantaneous non-uniform law in which the number of particles nucleated depends on the undercooling and is classified by grain size. An example for the latter one is illustrated in Fig. 2. 7, issued from Marie Bedel's Ph.D thesis [Bedel, 2014] to which further description is referred. Illustrative figures issued from [Bedel, 2014].

Methodologically initializing one stage with the solution obtained from the other stage creates a consistent connection between these two contributions in order to determine an entire evolution of variables. This approach provides an efficient way when numerically solving a system of non-linear equations that describe strongly coupled multi-scale physical phenomena in solidification. Besides considering the solid motion together with other phenomena (heat transfer, solute transport, liquid movement), the model can treat both globular and dendritic grains, simultaneously incorporating a refine solute diffusion length [Tveito et al., 2011] as well as tracking the columnar-to-equiaxed transition [Leriche et al., 2015]. It was successfully applied to predict the macrosegregation of large industrial ingots, such as a 3.3t steel ingot [Combeau et al., 2009][Založnik andCombeau, 2009a], a 6.2t steel ingot [Kumar et al., 2012], a 65 steel ingot [Combeau et al., 2012]. Fig. 2. 8 presents simulation results for the 3.3t steel ingot at an intermediate time, issued from [Combeau et al., 2009], with three case studies, including Case 1 -fixed solid phase, Case 2 -solid transport by assuming globular grains and Case 3 -solid transport with dendritic grains. It can be seen that the solid motion significantly impacts the solidification sequence and the flow patterns. Off these studies, the macrosegregation predicted by the dendritic model is very close to experimental measurements, as shown in Fig. 2. 9. RHS : segregation ratio map and liquid velocity vectors/streamlines.

Case 1 Case 2 Case 3 Experiment Segregation ratio along center

Fig. 2. 9 -Final macrosegregation ratio maps and segregation ratio along the center line, issued from [Combeau et al., 2009]. Chapter 2 -Literature Reviews 2.4.4 Models at CEMEF Developed in FE formulation, the works related to solidification modeling are implemented in two codes at CEMEF. The first one in Fortran is called R2SOL and used for two-dimensional solidification simulations. The second one is based on the computational platform named CIMLib, a C++ parallel computing library that can address two-dimensional and threedimensional geometries.

The first works carried out in R2SOL consisted in developing the resolution for the Navier-Stokes equation using a velocity-pressure P2+/P1 formulation and an arbitrary Lagrangian-Eulerian method, then implementing the resolution for energy and solute transport equations as well as coupling the resolutions [Gaston, 1997] [Gaston et al., 2000]. After that, subsequent developments investigated various issues and comprised an implementation of a P1+/P1 formulation [Heinrich, 2003], optimizations by studying coupling approaches [Liu, 2005] and by using an adaptive anisotropic remeshing [Gouttebroze, 2005]. The transport of solid phase was not addressed through these developments. However, the deformation of the solid phase and its interaction with the liquid phase gave rise to the development of a solid-liquid concurrent resolution scheme. The approach was limited to the packed bed regime in which the solid phase was treated as a continuum following a viscoplastic compressible constitutive equation, and the liquid flow obeyed a permeability law such as Carman-Kozeny. The concurrent resolution consisted in solving solid velocity, liquid velocity and pressure through a velocity-pressure formulation, using a P1+/P1+/P1 finite element formulation. Details of these works are presented in [Bellet, 2005] [Bellet and Fachinotti, 2006] [Bellet, 2007]. One of the applications was to study the macrosegregation due to solid deformation in continuous casting. Fig. 2. 10 presents simulation results on a thick slab of steel (Fe-0.1%C) in the late stage of solidification [Fachinotti et al., 2006]. Fig. 2. 10-a) shows the relative velocity vector in the mushy zone with two solidus isotherms, displaying the interaction tendency between the solid and liquid phases. Fig. 2. 10-b) depicts the divergence distribution of the intrinsic solid velocity (negative values in blue color indicating contraction zones and positive values in red color for expansion zones). As it can be observed that the intrinsic solid velocity divergence predicted is consistent with the deformation and strain states involved: positive values are found in the stretched mushy core between two face-to-face rolls and negative values in compressed regions between two successive rolls. The segregation map in Fig. 2. 10c) exposes an alternative formation of positive and negative segregations at the center of the slab: negative values (blue color) are predicted for the zones getting out a roll pair and positive values (red color) for the zones entering a roll pair. Figures issued from [Fachinotti et al., 2006].

Related to solidification modeling realized in CIMLib, two recent works including the theses of Rivaux [Rivaux, 2011] and Carozzani [Carozzani, 2012] will be presented. They can be considered as the starting basis for our own developments.

The first contribution aims at studying macrosegregation induced by thermo-solutal liquid convection and solid deformation. This last aspect is specific to the continuous casting process due to deformation of the solidified shell either in the mold or in the machine when passing by the support rolls. The model using the enthalpy-based formulation was first implemented without considering solid deformation, in which the liquid is assumed to be locally well mixed; the microscopic evolution follows the lever rule (or the Scheil approximation). The average solute composition carried along in the liquid flow is determined by solving the summed-up solute equation of the solid and liquid phases, thus there is no explicit description of interfacial exchanges. After that, from the fixed and rigid solid case, in order to account for the solid deformation, the model is extendedly coupled with the resolution of solid mechanics, performed with Thercast®, providing the solid velocity field which is introduced into the Darcy term using the Carman-Kozeny law and in the terms arising from the relative velocity, 〈 〉 〈 〉 , in the solute transport and energy equations. The second work [Carozzani, 2012] is based on the fixed solid assumption. The model was developed to predict grain structures using the Cellular Automaton -Finite Element (CAFE) method, originated from the work of Gandin and Rappaz [Gandin and Rappaz, 1994], in which the mutual interpolations between the FE and CA resolutions allow relating microstructures and macroscopic features. The quantities determined at FE nodes, i.e. enthalpy and solute concentration, are interpolated at cells of the CA grid superimposed onto the FE mesh, at which the nucleation, solute diffusion and grain growth are to be modeled. Once determined, the evolution of variables at CA cells is deduced to obtain the corresponding quantities at FE nodes by an interpolation procedure. We introduce here one of the case applications studied in this thesis to which detailed information and analyses are referred. The case study was investigated to simulate solidification of a parallelepiped cavity of 100 mm (length) x 60 mm (height) x 10 mm (with), containing a Sn-3wt%Pb alloy, as shown in Fig. 2. 11-a). The cooling is controlled by imposing the temperature on two narrow vertical sides, other walls are considered as to be adiabatic. The corresponding experiment benchmark is described in [Hachani et al., 2012].

As presented in Fig. 2. 11-b), the cooling curves predicted by CAFE model (plotted by colored solid curves) retrieve experiment measurements (displayed by black dashed curves), taking at four positions L21, L24, L27, L30 corresponding to four points at 5, 35, 75, 95 mm from the left side at the mid-height of the cavity. Fig. 2. 12 shows a good agreement of the grain structure in the cavity's median plane between numerical predictions given by CAFE simulation (Fig. 2. 12-a) and experimental results (Fig. 2. 12-b). In addition, it is noticed that there are equiaxed grains visualized in the zone marked by a dashed ellipse in Fig. 2. 12-b). The effect of crystal transport during solidification would be accounted for to properly simulate the process. Fig. 2. 12 -Grains structure maps in the cavity's median plane a) predicted by 3D CAFE simulation [Carozzani, 2012] b) obtained from experiment analysis [Hachani et al., 2012] It can be realized that the solid transport in solidification has been almost modeled by using finite volume methods. Till now, modeling with consideration of settlement and floating of solid grains developed on the finite element framework has been limited while the solid movement is known as one of the crucial factors causing macrosegregation and structural heterogeneities, especially in producing large castings. Being inspired by successful applications of the FV model in SOLID software which uses a splitting method to effectively solve highly dependent multi-scale phenomena during solidification, we are interested in adapting and implementing this scheme into the FE context to improve the models at CEMEF in which modeling of the motion of equiaxed grains has not been considered yet. Chapter 2 -Literature Reviews

Summary

In this chapter, we first reviewed the multi-scale modeling and its application to the solidification context, using the volume averaging technique whose principles were summarized. In the following, the development of solidification models was presented, focusing on those using the volume averaging technique with the presence of solid transport.

Lastly, models at CEMEF were summarized, on which our work will be developed consisting in adapting the FV growth-transport model in SOLID software into the FE framework. The realization steps can be outlined as follows.

 First considering a stationary system, in order to describe the formation and evolution of equiaxed grains, we will implement the nucleation-growth stage from the FV model which is able to account for an incomplete mixing solute diffusion within both solid and liquid phases.

 The next operation is to take into account the melt convection but still no solid motion. We will additionally solve the liquid momentum equation and the transport equations due to the liquid movement driven by thermo-solutal buoyancy.

 Following this work, the solid motion will be integrated into the model through two principal steps: first the FE resolution for purely advective equations will be established without grain growth. Then the implementation for the complete model will consist in combining solid transport with growth stage.

 2D numerical simulations will be performed and compared to reference results in each step to verify our implementations and our resolutions.

 After that, based on this model for two phases, the development will be extended to a three-phase model which allows accounting for the morphology of grains.

 Finally, 3D simulations will be studied and industrial applications will be carried out for the 3.3t steel ingot produced by Aubert & Duval (Fig. 2. 8). Chapter 2 -Literature Reviews

Résumé en français

Dans ce chapitre, nous avons tout d'abord présenté le principe de la modélisation multiéchelle et de son application à l'étude de la solidification en utilisant la méthode de la prise de moyenne volumique dont les principes ont été résumés. La partie suivante a porté sur les modèles de solidification qui utilisent la technique de la prise de moyenne et prennent en compte le transport des grains équiaxes. Enfin, nous avons fait un résumé des modèles préexistants au CEMEF, auxquels sera intégré notre travail dont les étapes de réalisation sont les suivants.

 Dans un premier temps, sans mouvement des phases, les modélisations de la germination et de la croissance des grains équiaxes seront mises en oeuvre avec un modèle de croissance qui permet de prendre en compte une diffusion incomplète du soluté dans les phases de liquide et de solide.

 A l'étape suivante, la convection thermo-solutale sera prise en compte mais sans mouvement du solide. Pour ce faire, le modèle de croissance sera combiné à la résolution par éléments finis des équations de conservation de la quantité de mouvement et de transport du soluté associées à l'advection du liquide.

 Subséquemment, l'intégration du transport des grains solides dans le modèle sera effectuée en deux étapes principales: premièrement, la résolution par éléments finis des équations du transport pur sera mise en place en négligeant la croissance des grains; deuxièmement, l'implémentation du modèle complet, constitué du transport et de la croissance, sera menée.

 Au cours du fil de ces travaux, des calculs numériques correspondants seront établis et réalisés pour vérifier nos implémentations.

 Après ces étapes, à partir du modèle à deux phases, le développement sera étendu à un modèle à trois phases qui permet de prendre en compte la morphologie des grains dendritiques.

 La partie finale s'intéressera aux simulations en 3D et aux applications industrielles. Chapter 3 -Finite Element Implementation and Adaptation Chapter 3

Finite Element Implementation and Adaptation

In this chapter, we will present our model adaptation and its numerical implementation at CEMEF. We first focus on the adaptation to the FE context of the model for the macroscopic phenomena including the transport of phases, grain density and solute, neglecting terms related to the contribution of microscopic processes; then we present the FE solvers for the conservation equations of energy and the one of solid and liquid momentums. After that, the microscopic evolution models, nucleation and growth stages, are described. Then the total coupling between these stages will be addressed and discussed. Finally, from the present model with two phases, we implemented its extension to the three-phase model issued from the thesis work of Marie Bedel [Bedel, 2014], which allows accounting for the morphology of dendritic grains.

Two-Phase Solidification Model

Before studying each stage of the transport-growth model, we introduce first several points related to the use of the splitting scheme, applied to the solidification context, developed by Založnik and Combeau [Založnik and Combeau, 2010a].

1. Using this scheme, the evolution of phase fraction, grain density and solute concentrations is solved by splitting the transport and growth stages since it results from the progress of phenomena at very different scales of time and space (from a few seconds for nucleation and growth processes occurring at dendrite arms, to a few hours for those of transport taking place on systems sized in meters). The resolutions for the energy and momentum conservation equations are effectuated outside the splitting. An illustration of the resolution algorithm is schematized in Fig. 3. 1. Chapter 3 -Finite Element Implementation and Adaptation Fig. 3. 1 -Schematic of the resolution algorithm using the splitting method.

2. An iterative procedure of fixed point type is utilized to solve the entire set of nonlinear equations. The notation ( ) is hereafter referred to a previous iteration relative to the actual one considered ( ).

3. Additionally, when applying this splitting method, the use of phase fractions associated to velocity fields should be noticed. As introduced in Section 2.4.3, in the transport stage where only the advective terms are accounted for, the transport of phase 𝛼 (𝛼 represents either the solid or liquid phase, i.e. 𝛼 or ) is thus described by Eq.( 37) from which a continuity equation ( 38) can be deduced by summing equations written for the solid and liquid phases, with assumption of constant and identical density for each phase. It is expressed as a function of phase fractions obtained from the transport stage, 𝑔 . Moreover, to determine which value of phase fractions should be used in the advective terms of the transport equation for a general quantity in phase 𝛼, 〈 〉 , it is attributed by an undetermined value 𝑔 as written in Eqs.(39).

On the one side, these equations can be further performed as Eqs.( 40). On the other side, they must satisfy Eqs.( 41) since the quantity must be conserved during the process. Substituting Eqs.( 37) and ( 41) into Eqs.( 40), we can deduce that Eqs.( 40) are fulfilled if 𝑔 𝑔 . Consequently, the phase fractions computed from the transport stage 𝑔 are used in the continuity equation as well as in advective terms of all equations to ensure the conservation of quantities and the consistent scheme implemented.

(𝑔 ) (𝑔 〈 〉 ) (37)

(𝑔 〈 〉 ) (𝑔 〈 〉 ) (38)
(𝑔 〈 〉 ) (𝑔 〈 〉 〈 〉 ) (39) ⇔ 𝑔 〈 〉 〈 〉 [ 𝑔 (𝑔 〈 〉 )] ((𝑔 𝑔 )〈 〉 〈 〉 ) (40) 〈 〉 〈 〉 〈 〉 〈 〉 (41)
where the subscripts and indicate the transport stage and an undetermined value, respectively.

Modeling of Macroscopic Transport Phenomena with FEM

Pure Transport Equations and Artificial Diffusion

To model the transport of particles, like solid grains, Lagrangian-Lagrangian and Eulerian-Lagrangian methods are a useful analysis tool for studying in detail interaction processes, e.g. particle-particle or particle-liquid interactions. In Lagrangian-Lagrangian approaches, the fluid dynamics is simulated by the Smoothed Particle Hydrodynamics (SPH) method using a set of moving makers in a Lagrangian framework. Each marker has an associated kernel function, W(r, h), defined over a support domain, S(h), where r is the distance from the SPH marker and h is a characteristic length that defines the kernel smoothness, as presented in Fig. 3. 2-a. The optimal value for h is such that every fluid particle has about 30-50 neighbors within smoothing volume. Each marker carries fluid properties like pressure, density, etc. For the solid phase, its behavior results from the sum of forces and moments acting on it, which are described by the Newton-Euler equations for each single solid body. The dynamics of the two phases, fluid and solid, are modeled through Boundary Condition Enforcing (BCE) markers placed on and close to the solid surface, as schematized in Fig. 3 and SHP neighbors (white circles). Figures issued from [Pazouki and Negrut, 2015].

Further details related to this technique can be found in [Abel, 2011] [Bauer and Springel, 2014] [Pazouki and Negrut, 2015].

Another approach to simulate the particles transport is using Eulerian-Lagrangian methods in which the integration of transport equations is performed for each individual particle moving in a continuous fluid. The motion of particles is modeled by accounting for various governing forces that depend on study cases. Different applications were investigated by using this method: studying the behavior and trajectory of growing solid grains transported in the bulk liquid [START_REF] Stomp | Hydrodynamics of moving structures during a solidification process[END_REF]; simulating the motion of bubbles (light particles) and inclusions (heavy particles) in a rotating fluid flow during the Rotating Continuous Casting process [Fidel et al., 2011]; predicting the entrapment of slag inclusion and bubbles in the continuous casting od steel slabs in [Zhang et al., 2004] [Yuan and Thomas, 2006] [Pfeiler et al., 2007] [Thomas et al., 2014].

However Eulerian-Eulerian approaches are commonly preferred in industrial applications because on the one hand they can provide essential information, e.g. the distribution and velocity of particles and on the other hand they are more advantageous for dealing with high Reynolds numbers and a large number of particles. Our development aims mainly at studying the effects of liquid and solid movements on the macrosegregation formation during solidification, which is applied first to average-sized samples and ultimately to large industrial castings. With this objective in mind, we use the Eulerian-Eulerian method in which the solid is considered as a continuum composed by a large number of particles; it is characterized by a fraction, a grain density (or a mean radius) and its solute content. The solid phase moves in a continuous fluid and solid-liquid interactions are introduced through drag forces in the momentum equations.

In the transport stage, neglecting the contribution of microscopic phenomena (which will be treated separately, in a splitting manner), the transport of quantities including the solid phase, grain density and solutes are described by Eqs.( 42), ( 43) and ( 44) respectively.

𝑔

(𝑔 〈 〉 ) (42)

( 〈 〉 ) (43)
(𝑔 〈 〉 ) (𝑔 〈 〉 〈 〉 ) (44)
where the subscript indicates the transport stage, 𝑔 is the solid fraction, is the grain density, 〈 〉 is the solute concentration in phase 𝛼, 〈 〉 is the velocity of phase 𝛼.

Due to the lack of diffusive effects, FEM resolutions for purely advective equations or hyperbolic equations as those above manifest numerical problems as instabilities and possibly non-unique weak solutions.

Concerning convection-dominated problems, numerical instabilities can be effectively decreased by using upwind methods since the discretization of advective terms with upwind difference schemes introduces by itself an artificial diffusion that smears oscillations, e.g. 1D studies introduced in [START_REF][END_REF], 2D problems specified in [Henrich et al., 1977] and a more general upwind technique applied to both 1D and 2D cases presented in [Hughes, 1978]. However resulting solutions are excessively over-diffused. For this reason, FE stabilizations have been developed to overcome difficulties encountered in solving convection-dominated equations. A great pioneering investigation for this issue was devoted by Hughes and co-workers whose works have been popularly applied in various domains. For incompressible flows, streamline upwind methods were proposed to eliminate the crosswind diffusion by constructing an artificial diffusion operator acting only in the flow direction [START_REF] Hughes | [END_REF]. Instead of modifying the diffusion flux like in these approaches, which can cause incoherence problems, the Streamline Upwind/Petrov-Galerkin (SUPG) method was developed [START_REF] Hughes | [END_REF], in which the introduction of streamline upwind diffusion can be considered as a modification of the interpolation functions of the standard Galerkin method. This approach ensures that the residual of the associated weak formulation vanish with the exact solution. Further numerical application examples for these methods were reported in [Brooks and Hughes, 1982]. FE stabilization techniques dedicated to incompressible flows can be consulted in [Tezduyar andHughes, 1982] [Hughes et al., 1987].

An alternative approach developed to solve hyperbolic equations is the flux corrected transport scheme which provides a compromise solution for problems of spurious oscillations and excessive numerical diffusion [START_REF] Kuzmin | [END_REF]. Unlike methods derived from the Galerkin method, this technique involves modifying the algebraic matrices of system in order to have desired effects, i.e. manipulating on matrices to be solved, first adding diffusion to eliminate negative off-diagonal coefficients of mass matrices, then removing excessive diffusion, but it is applicable to the case of divergence-free velocities.

Firstly, through the above studies, it is worthy to notice that the presence of diffusion, whether being generated due to computational procedures or an artificial introduction, it seems essential when numerically solving convection-dominated equations. In our case, we solve the transport equations ( 42)-( 44) in which there is no diffusion, numerical resolutions can produce unphysical instabilities. Working on the existing FE solver which was developed for a convection-diffusion-reaction equation and already integrated the SUPG stabilization scheme, we propose to introduce an artificial diffusion coefficient into hyperbolic equations to be solved. Secondly, while the total average velocity, including both solid and liquid velocities, is divergence-free, as described by Eq.( 38), the transport velocities belonging to each phase are non-divergence-free fields Eqs.(37). In addition, complexities come from discontinuities due to the solid packing phenomenon taking place at a transition zone between two regions: in one region small grains move freely in the liquid phase (the solid fraction is smaller than a packing fraction); in another region large grains settle, accumulate to create a packed zone and become stationary (the solid fraction is higher than the packing limit). It can be seen that discontinuities are related to quick changes in the velocity and fraction of the solid phase while it forms a packed bed.

Accordingly, we developed in our context an adaptive artificial diffusion in order to stabilize FE resolutions while not smearing results unreasonably. The coefficient , as expressed in (45), consists of the gradient of solid velocities (the average solid velocity divergence 〈 〉 and the intrinsic solid velocity divergence 〈 〉 ) which contain information about the variation of the related quantities (both velocity and fraction of the solid phase). It allows detecting and reducing discontinuities where the solid phase is being packed. Thus the artificial diffusion is not solved everywhere. It is only restricted to critical zones where the solid velocity is changing abruptly. An alternative option could be obtained by using the intrinsic solid velocity divergence 〈 〉 and the solid fraction gradient 𝑔 . However, the former expression provides a further advantage as it vanishes in the absence of solid transport; the model is strictly valid for the case without solid motion. The formulation of this coefficient is based on the principle from the work of Cook and Cabot [Cook and Cabot, 2003] in which the authors developed an artificial non-linear diffusion using the entropy gradient to treat issues associated with discontinuities of the temperature and mass fraction in supersonic reacting flows. The constant factors 𝛼 and in this coefficient (45), allows controlling the amount of diffusion in a direct way. Moreover, this approach provides more flexibility compared to the implicit diffusion introduced by upwind methods.

( ) (𝛼 | 〈 〉| | 〈 〉 |) (45)
where denotes the characteristic size of each mesh element in the direction of the solid velocity , the formulation is issued from [Fachinotti et al., 2006].

‖ ‖ ∑ | | (46)
with is the interpolation function associated with node is the solid velocity at the center of the element It should be noticed that using a single artificial diffusion coefficient in the different transport equations is necessary, in order to ensure the consistent transport of related quantities and the conservation of mass and energy.

Consequently, by adding the artificial diffusion coefficient , the transport equations for solid phase, grain density and solute as issued from ( 42)-( 44) are modeled by Eqs.( 47)-( 49).

𝑔 〈 〉 𝑔 𝑔 〈 〉 ( 𝑔 ) (47)
〈 〉 〈 〉 ( ) (48)
(𝑔 〈 〉 ) 〈 〉 (𝑔 〈 〉 ) 𝑔 〈 〉 〈 〉 ( (𝑔 〈 〉 )) (49)
Assuming constant densities for the solid and liquid phases, it is possible to transform Eq.( 49) in the form:

𝑔 [ 〈 〉 〈 〉 〈 〉 ] ( 𝑔 〈 〉 ) 〈 〉 ( 𝑔 ) ( 〈 〉 𝑔 ) (50)
The third term in Eq.( 50) is identified by using Eq.( 47). Its introduction allows reducing the need of calculating the divergence of potentially discontinuous velocity fields, which is known as being not favorable for the FEM. The solute transport equation is thus to be implemented with the formulation of Eq.( 50).

Since Eqs.( 47), ( 48) and ( 50) are of similar general form, as given by Eq.( 51), we summarize here the discretization implemented for Eq.( 51) and we can obtain the discretization of the above equations by replacing the variable and the factors 𝛼 𝛼 𝛼 with the corresponding quantities, as summarized in Tab.3. 1. Details of FE discretization and resolution methods can be consulted in [Rappaz et al., 2003].

Tab.3. 1 -Corresponding quantities to those of the generic equation ( 51)

Eq. 𝛼 𝛼 𝛼 (47) 𝑔 〈 〉 (48) 〈 〉 (50) 〈 〉 𝑔 𝑔 ( 𝑔 ) ( 〈 〉 𝑔 )

Let

, where is the number of space dimensions, the studied domain which is bounded by the border ,

𝛼 ( ) (𝛼 ) 𝛼 (51)
with the boundary conditions : on and on (52

)
where is the variable; the transport velocity; 𝛼 the diffusion coefficient; 𝛼 the reaction term; the source term; and are the boundary's subsets on which the Dirichlet and Neumann conditions apply respectively.

The corresponding generalized solution is supposed to satisfy the strong form of Eq. ( 51) for sufficiently smooth data and to exist even if the divergence theorem is not applicable and the underlying conservation law holds only in an integral sense. The solution must belong to the space of functions which are continuous with continuous partial derivatives of first and second order. In order to broaden the class of admissible functions, it is worthwhile to consider an integral or weak form of the conservation law. An approach to obtain the derivation of weak forms for a given partial differential equation is the method of weighted residuals. Consider the residual of Eq.( 51):

𝛼 ( ) (𝛼 ) 𝛼 (53)
The above weak-form equation is solved by the FEM, which consists in approximating the continuous spaces and , of infinite dimensions, by the corresponding discrete spaces and .

(

) 62
The domain is decomposed into small elements and the discretized domain is defined as:

⋃ (63)
where denotes the domain of element and the set of elements .

The variable in the studied domain is approximated by the following sum expressed over Nn nodes.

∑

(64) where is the interpolation functions of the variable associated with node and is the value of at node . By using an implicit temporal discretization and applying the divergence theorem [Rappaz et al., 2003] to the diffusion term, Eq.( 58) is expressed as:

∫ 𝛼 ∫ 𝛼 ( ) (65)
∫ 𝛼 ( ) ∫ 𝛼 ( ) ∫ 𝛼 ∫
Here the repeated index stands for a summation and denotes the nodal value of at the beginning of the time step when discretizing the time derivative . It should be noticed that for linear elements and on boundaries. Chapter 3 -Finite Element Implementation and Adaptation

The discretization of the problem leads to Eq.( 66) being expressed with matrix formulation.

[ ]{ } { } { } (66)
where { } is the vector of variable to be determined

[ ] ∑ ∫ 𝛼 (67) [ ] ∑ ∫ 𝛼 (68) [ ] ∑ [ ∫ 𝛼 ∫ 𝛼 ] (69) [ ] ∑ ∫ 𝛼 (70) { } ∑ ∫ (71) { } ∑ ∫ (72)
We have seen the FE implementation for the equations in the transport stage accompanied with particular numerical treatments. The fixed point procedure presented at the beginning of this chapter is now updated as shown in Fig. 3. 3, by introducing the transport solutions determined for an actual iteration , including the solid fraction 𝑔 , the grain density and the intrinsic compositions 〈 〉 . Next we will study the FE solver for the energy equation, recalling that this equation is solved outside the splitting scheme, according to the developed method in [Založnik and Combeau, 2010a]. Chapter 3 -Finite Element Implementation and Adaptation 3.1.1.2 Energy Equation 3.1.1.2-a) Description of the Implemented Equation

The averaged energy equation for phase 𝛼 (recalling that 𝛼 represents either the liquid or the solid phase in our problem) is given as the following equation, previously introduced as Eq.( 14) in Chapter 2:

(𝑔 〈 〉 〈 〉 ) (𝑔 〈 〉 〈 〉 〈 〉 ) (𝑔 〈 〉 ) (73)
where 〈 〉 is the intrinsic average of the heat flux vector in phase 𝛼; and are the interfacial heat transfers due to phase change, diffusion and nucleation respectively.

Two considerations are taken into account for Eq.( 73).

i. Firstly, as demonstrated above, when using the splitting scheme, the phase fractions in advection terms must be the values obtained from the transport stage. Eq.( 73) is thus expressed as ( 74).

(𝑔 〈 〉 〈 〉 ) (𝑔 〈 〉 〈 〉 〈 〉 ) (𝑔 〈 〉 ) (74)
ii. Secondly, for numerical stabilization, the origin transport equations ( 42)-( 44) are transformed to Eqs.( 47)-( 49) by introducing the same artificial diffusion coefficient. This operation implies that the related quantities (phase fraction, grain density and solute concentrations) would be modified in a same manner by effects of numerical diffusion relative to solutions without diffusion. With respect to the treatment of the transport equations, the energy transported in each phase must be also modified by the same artificial diffusion in order that the transports of variables are consistent. Therefore the modified energy equation for each phase is given by:

(𝑔 〈 〉 〈 〉 ) (𝑔 〈 〉 〈 〉 〈 〉 ) ( (𝑔 〈 〉 〈 〉 )) (𝑔 〈 〉 ) (75)
Furthermore, other assumptions when solving the energy equation are listed below:

 The densities of the solid and liquid phases, respectively 〈 〉 and 〈 〉 , are constant and equal, denoted :

〈 〉 〈 〉 (76)
 The heat diffusion is modeled by the Fourier law:

〈 〉 〈 〉 〈 〉 (77)
 The temperature of both phases is assumed to be identical due to a large thermal diffusivity.

〈 〉 〈 〉 (78)

 Heat exchanges at the solid/liquid interface are assumed to be in equilibrium.

(79)

By using the assumptions (77)-( 79) and adding up Eqs.(75) for the solid (𝛼 ) and liquid (𝛼 ) phases, we obtain the mixture energy equation:

[ 〈 〉 (𝑔 〈 〉 〈 〉 𝑔 〈 〉 〈 〉 ) ( (𝑔 〈 〉 𝑔 〈 〉 ))] (〈 〉 ) (80)
where

〈 〉 𝑔 〈 〉 𝑔 〈 〉 (81) 〈 〉 𝑔 〈 〉 𝑔 〈 〉 (82)
For convenience of numerical implementation, the advection term is expressed using the average values of enthalpy and velocity. The corresponding correction is put on the RHS as a source term, noticing that since the average velocity is a divergence-free field, the term associated to the convection on the LHS 〈 〉 〈 〉 is equal to (〈 〉〈 〉) on the RHS. In addition, in the added diffusion term, the mixture enthalpy using the transport phase fractions is approximated by the unknown average enthalpy. These transformations lead to Eq.( 83). Chapter 3 -Finite Element Implementation and Adaptation

[ 〈 〉 〈 〉 〈 〉 ( 〈 〉)] (〈 〉 ) [ (𝑔 〈 〉 〈 〉 𝑔 〈 〉 〈 〉 ) (〈 〉〈 〉)] (83)
where

〈 〉 𝑔 〈 〉 𝑔 〈 〉 (84)
Compared to Eq.( 51) presented above, Eq.(83) presents differences. It contains two unknowns in the LHS terms, the average enthalpy 〈 〉 and the temperature , which are strongly and non-linearly related together as well as to other variables in the solidification context. An iterative resolution is operated in order to determine a solution for these variables and we present here the approach's principle.

For simplicity, the RHS of Eq.( 83) is denoted as (85). In the iterative procedure, this source term is computed using the latest known values.

[ (𝑔 〈 〉 〈 〉 𝑔 〈 〉 〈 〉 ) (〈 〉〈 〉)] (85)
The strong formulation is defined as to finding 〈 〉 over a bounded domain and its boundary such that:

[ 〈 〉 〈 〉 〈 〉 ( 〈 〉)] (〈 〉 ) (86)
with the boundary conditions

(Dirichlet condition) 〈 〉 (Neumann condition) 〈 〉 ( ) (Fourier condition)
where is the unit outward normal to .

Using the interpolation and test functions in the corresponding functional spaces, which are presented in Section 3.1.1 with the SUPG method, the weak formation consists in finding

〈 〉 such that ∫ { [ 〈 〉 〈 〉 〈 〉 ( 〈 〉)] (〈 〉 ) } (87)
By applying the boundary conditions (considering the Diriclet and Fourrier conditions) to Eq.( 87) and using the divergence theorem, we obtain Eq.( 88) after several transformations.

∫ 〈 〉 ∫ 〈 〉 〈 〉 ∫ 〈 〉 (88) ∫〈 〉 ∫ ( ) ∫ 3.1.1.2-b) Discretization
The studied domain is approximated by the discretized domain , the enthalpy 〈 〉 by The approximated values of enthalpy, 〈 〉 , and temperature, , are evaluated by the following sums expressed over the nodes.

〈 〉 with 〈 〉 { ( ) | ( ) | ( )}
〈 〉 ∑ 〈 〉 ∑ (89)
By using an implicit temporal discretization and choosing at nodes, we have for any node index (here the repeated index stands for a summation):

∫ 〈 〉 〈 〉 (90) ∫ 〈 〉 (〈 〉 ) ∫ (〈 〉 ) ∫〈 〉 ( ) ∫ (( ) )
∫ Chapter 3 -Finite Element Implementation and Adaptation

Eq.( 90) can be expressed as a vector equation using matrix formulation as:

[ ]{ } [ ]{ } { } { } { } (91)
where the different matrixes and vectors result from an assembly procedure processed on all finite elements:

{ } and { } are the vectors of the unknowns including enthalpy and temperature, respectively. In the present work, we solve for the set of nodal enthalpies which minimizes the residual vector ( 99). An alternative solution would consist in choosing nodal temperatures as principal unknown [START_REF] Saad | [END_REF].

[ ] ∑ ∫ (92) [ ] ∑ ∫ 〈 〉 (93) [ ] ∑ ∫ (94) [ ] ∑ [ ∫〈 〉 ∫ ] (95) { } ∑ ∫ 〈 〉 (96) { } ∑ ∫ (97) { } ∑ ∫ (
{ } [ ]{ } [ ]{ } { } { } { } (99)
As mentioned at the beginning of this chapter, due to a close relationship of all relevant variables, an iterative scheme -a fixed point approach -is implemented to solve the set of related equations, which is also effectuated at the same time to iteratively solve the non-linear energy equation. Therefore, the same notations are defined to call iterations, including for an actual iteration and for a previous one.

In each actual iteration , we determine the enthalpy, { } , or the variation of the enthalpy, { } ( { } { } ), which satisfies Eq.( 100), which derives from a first order development of the vector { } :

{ } { } [ ] ( { } { } ) { } [ ] { } (100) ⇒ [ ] { } { } (101) Knowing that [ ] [ ̃] (102)
where the tensor [ ̃] is expressed with its components :

[ ̃] [ ] [ ]
where there is no sum over (103) we rewrite Eq.( 101) as For simplicity, the terms on the RHS of Eq.( 105) can be denoted as follows:

[ ̃] { } { } (104) ⇒ [ ̃] { } [[ ]{ } [ ]{ } { } { } { }] (
{ } [ ]{ } (106) { } [ ]{ } (107) { } [ ]{ } (108) { } [ ]{ } (109) ⇒ [ ̃] { } { } { } { } { } { } { } { } (110)
As { } { } do not change in iterations over a timestep, it is not necessary to indicate their number of iterations. Eq.( 110) is then rewritten as ( 111), which is also arranged in the same order as in the implemented code. Although the terms related to boundary conditions which are separately treated in the code, here they are presented here together with the others.

[ ̃] { } { } { } { } { } { } { } { } (111)
where ] are presented as (92);( 93);( 94) and ( 103) respectively, for which, the velocity and the conductivity are updated in each iteration by using the latest values.

[ ] [ ] [ ] [ ̃
( { } { } ) ∑ ∫ ( 〈 〉 〈 〉 ) (112) { } ∑ ∫ 〈 〉 (113) ( { } { } ) ∑ [ ∫ 〈 〉 〈 〉 ∫ 〈 〉 ∫ ] (114)
{ } { } are identical to (97) and ( 98).

The enthalpy obtained by solving the energy equation is put in the context of the fixed point scheme as presented in Fig. 3. 4. Fig. 3. 4 -Schematic of the resolution algorithm using the splitting method, introduced the energy solution.

Momentum Equations

Depending on the solid fraction, the behaviors of the solid and liquid phases are modeled in two distinguished regimes:

 If the solid fraction is smaller than a critical packing value (𝑔 𝑔 ), solid grains are considered to move freely in the melt. This defines the so-called slurry regime.  In the other case (𝑔 𝑔 ), solid grains form a porous packed bed and they become stationary. The velocity of the solid phase is then considered null :

. For this, we have the so-called packed bed regime (or porous regime).

Actually, the reason for choosing such a simple condition to signify the transition from the slurry regime to the packed bed regime may be questionable. In the present work, we keep the same criterion as that in [Založnik and Combeau, 2010a]. In the application examples, we will see that despite of its simplicity, this criterion seems enough to represent the progressive packing of solid grains after transport and settlement.

Several assumptions are supposed when solving the momentum equations as below: i. The phase densities are constant and equal, except for the buoyancy forces for which the liquid and solid densities are modeled as ( 115) and ( 116) respectively.

[ ( ) (〈 〉 )] (115) (116)
where is the density at a reference temperature, , and at a reference concentration, ; is the shrinkage coefficient.

ii. The pressures in the solid and liquid phases are assumed to be equal.

〈 〉 〈 〉 (117)
iii. The liquid behaves as a Newtonian fluid with a constant viscosity. By neglecting the interfacial momentum transfer due to phase change, the divergence of the average deviatoric stress tensor is modeled as:

〈 〉 ( (𝑔 〈 〉 ) (𝑔 〈 〉 ) ) (118)
iv.

In addition, the momentum transfer due to nucleation is also negligible relative to others, the equilibrium at the solid/liquid interface is thus described as:

(119)

According to [Založnik and Combeau, 2010a], the authors solve two separate equations to calculate the liquid velocity, i.e.:

 In the slurry regime, the equation to be solved is obtained by summing up the momentum equations of the liquid and solid phases in order that the momentum transfer due to interfacial interactions is not explicitly required to be computed.

 In the packed bed regime, only the momentum equation for the liquid is solved ( ) with the use of a Darcy term to model the interaction with the static solid.

In the FE context, calculating the velocity of each phase by solving two different equations on a discretized domain can provoke numerical problems because of assembling dissimilar and irregular terms (i.e. one equation for slurry-regime elements and another for their packedregime neighbors). In order to overcome this issue, we further analyze the situation as follows.

3.1.1.3-a) Conservation of the momentum for the liquid phase

Using the aforementioned assumptions, the liquid momentum equation which is presented by Eq.( 13) in Chapter 2 will be deduced as (120) for the slurry regime, (𝑔 𝑔 ).

[

(𝑔 〈 〉 ) (𝑔 〈 〉 〈 〉 )] (120) 𝑔 ( (𝑔 〈 〉 ) (𝑔 〈 〉 ) ) 𝑔
Additionally, when using the splitting method, the liquid fraction associated to the velocity on the LHS of Eq.( 120) must be the value coming from the transport stage with respect to the continuity equation as demonstrated at the beginning of this chapter. Meanwhile, the liquid fraction obtained from the growth stage is used for terms on the RHS of Eq.( 120) since these quantities are related to a global evolution including both transport and growth processes.

Accounting for these aspects, Eq.( 120) is deduced as:

[ 〈 〉 (〈 〉 〈 〉 )] (121) 𝑔 ( 〈 〉 〈 〉 ) 𝑔
where 〈 〉 denotes the product of the liquid fraction computed from the transport stage 𝑔 and the intrinsic liquid velocity 〈 〉 .

Similarly, the liquid momentum equation for the packed bed regime (𝑔 𝑔 ) is given as:

[ 〈 〉 (〈 〉 〈 〉 )] (122) 𝑔 ( 〈 〉 〈 〉 ) 𝑔
Noting that Eqs.( 121) and ( 122) only differ from each other in the momentum transfer term due to interfacial stresses, we can thus model the movement of the liquid phase with only one equation including both regimes by using a function, 𝛼 , switching between different regimes:

𝛼 { (𝑔 𝑔 ) (𝑔 𝑔 ) (123)
Then, it can be envisaged to solve for a 𝛼 -weighted sum of Eqs.( 121) and ( 122) on the whole domain including slurry regime and packed bed regime, as expressed by (124):

𝛼 ( ) ( 𝛼 ) ( ) (124)
For numerical implementation, since a sudden change between the two regimes using the above switching function may cause discontinuities, we use a less brutal switching function, defined as:

𝛼 ( [𝛼 (𝑔 𝑔 )]) (125)
where we choose 𝛼 , as plotted in Fig. 3. 5. On the other hand, in order to effectuate the proposed combination (124) of Eq.( 121) and Eq.( 122) in a consistent way, we perform the momentum transfers of drag forces in the general form (126) applicable to both regimes.

(𝑔 ) (〈 〉 〈 〉 ) (126)
where the subscript 𝑔 denotes either the slurry regime or the packed bed regime.

The expression of for each regime is deduced by using the formulations of momentum transfers due to interfacial stresses presented in the article [Založnik and Combeau, 2010a], which were introduced in Chapter 2 (Eqs.( 29)-( 36)).

For the slurry regime :

𝑔 ‖〈 〉 〈 〉 ‖ (𝑔 ) 𝑔 ( ) (𝑔 ) (127)
For the packed bed regime :

( 𝑔 ) (𝑔 ) (128)
with , and defined as specified in Chapter 2.

Then the general liquid momentum equation including both flow regimes is obtained by substituting ( 121),( 122) and ( 126) into ( 124):

[ 〈 〉 (〈 〉 〈 〉 )] (129)
𝑔 ( 〈 〉 〈 〉 ) 𝑔 where [𝛼 ( 𝛼 ) ](𝑔 ) (〈 〉 〈 〉 ) (130)
and are evaluated as ( 127) and ( 128) respectively.

Moreover, for the liquid momentum equation, the existing FE solver (referred to [Hachem et al., 2010]) was developed by choosing the average velocity, 〈 〉 𝑔 〈 〉 , as an unknown to be determined and the entire equation was divided by the liquid fraction for a numerical reason. Its implemented formulation is given as (131).

𝑔 [ 〈 〉 𝑔 (〈 〉 〈 〉)] (131)
( 〈 〉 〈 〉 ) 𝑔 (〈 〉 𝑔 〈 〉 )
where is presented in Chapter 2 (Eq.( 36)).

With respect to the previously developed solver, two considerations are proposed for the present work which consists in implementing Eq.( 129), as follows:

 the average transport velocity, 〈 〉 𝑔 〈 〉 , is selected as the unknown.

 since the viscosity and drag terms also contain the average velocity, so it is as well approximated by the average transport unknown, i.e.:

( 〈 〉 〈 〉 ) ( 〈 〉 〈 〉 ) (132) [𝛼 ( 𝛼 ) ]𝑔 (〈 〉 𝑔 〈 〉 ) (133)
Therefore, the final formulation of the liquid momentum equation implemented is given as:

𝑔 [ 〈 〉 (〈 〉 〈 〉 )] (134)
( 〈 〉 〈 〉 ) 𝑔 𝑔
where

𝑔 [𝛼 ( 𝛼 ) ](〈 〉 𝑔 〈 〉 ) (135)
Solving the system of equations including Eq. ( 134) and the continuity equation Eq.( 38), we determine the liquid velocity, 〈 〉 , and the pressure, .

3.1.1.3-b) Conservation of the momentum for the solid phase

Since the solid motion is only existent in the slurry regime, we only solve the momentum equation for the solid phase in this regime. Neglecting the inertial and viscous terms, Eq.( 13), as introduced in Chapter 2, reduces here to:

𝑔 𝑔 (136)
Substituting ( 126) into (136), we can deduce the solid velocity in the slurry regime as:

〈 〉 〈 〉 𝑔 ( 𝑔 ) (𝑔 ) (137)
To avoid a sudden change of the solid velocity when the packing process occurs, we also apply a gradual switch from the moving to packed state of solid grains by using the transition function, 𝛼 , as designated in (125). Therefore the solid velocity is expressed as:

〈 〉 𝛼 (〈 〉 𝑔 ( 𝑔 ) (𝑔 ) ) (138)
Furthermore, to ensure that moving solid grains cannot penetrate into the area where the solid has already packed, the velocity of the grains going toward the zone near the packing region is adjusted so that these grains can land smoothly in the packed bed, i.e. reaching 〈 〉 at the packing limit. In order to do so, we check the distance which grains can realize with their actual velocity in the direction of the gradient of the solid fraction, 〈 〉

, and compare it to that between them and the packing limit, . If the former is higher, it indicates that grains go so fast and they will enter into the packed zone, their velocity should be calibrated so that they just arrive at the limit. This adjustment is illustrated by the scheme in Fig. 3. 6. Lastly, the fact that the same modified solid velocity ( 138) is used to model the momentum transfer induced by interfacial stresses in the slurry regime for the momentum equations of both phases, ensures the consistency in modeling solid-liquid interactions.

The solutions of the momentum equations are arranged in the present fixed point scheme, as illustrated in Fig. 3. 7.

Before adjustment

After adjustment Fig. 3. 6 -Schematic of the adjustment for the solid velocity. Fig. 3. 7 -Schematic of the resolution algorithm using the splitting method, introduced the solutions of the momentum equations.

These above developments correspond to adaptations realized for the transport stage from the FV model. The implementation of the FE resolutions has been detailed, in which solving purely advective equations and numerical treatments for discontinuities are delicate issues, which require careful investigation. In the next section, the relevant equations used for the description of solidification evolution during the nucleation and growth stages will be presented.

First nucleation occurs if local favorable conditions are satisfied. In the present work, instantaneous nucleation at the liquidus temperature is considered, in which a predefined number of grains are distributed at places where the temperature is smaller than the liquidus temperature. Nucleation is limited to positions where it did not happen yet. In addition, a redistribution of grains will be effectuated at sites that are empty of grains due to transport. Under these assumptions, the rate of grain generation due to nucleation ̇ can be mathematically expressed as ( 140). This nucleation model can be improved by using more realistic laws, e.g. those in which grain densities are classified by size and generated depending upon the undercooling [Bedel, 2014]. The following set of equations is considered:

̇ (139) ̇ { ( ) (∫ ̇ ) (140) 𝑔 (141) (𝑔 〈 〉 ) (142)
(𝑔 〈 〉 ) ( 143
)
where the subscript indicates the values obtained from nucleation; ̇ is the rate of grain generation due to nucleation; is a micro time-step; is the temperature at which nucleation occurs;

is the rate of solid formation associated with nucleation; is the composition at the solid-liquid interface in the solid side;

is the interfacial solute transfer of phase 𝛼 due to grain nucleation.

Next, the growth process is considered at each node as follows, in which the variables continue evolving from their states at the end of nucleation.

𝑔

(144)

𝑔 (145) (𝑔 〈 〉 ) 〈 〉 〈 〉 (146) (𝑔 〈 〉 ) 〈 〉 〈 〉 (147)
where the subscript 𝑔 indicates the values obtained from the growth stage. Relevant quantities were introduced in Chapter 2. For recalling, is the mass flux of phase 𝛼 across the phase interface due to phase change, is the interfacial solute transfer of phase 𝛼 due to phase change, is the solute diffusion coefficient in phase 𝛼, ⁄ ( ) [m - 1 ] is the interfacial area concentration (identical for the two phases present in the REV considered);

( 𝑔 ⁄ ) ⁄ [m] is the grain radius; is the characteristic solute diffusion length in phase 𝛼 (details in Annex);

is the solute partition coefficient.

Solving Eqs.( 144)-( 147) allows determining the quantities of interest consisting of the solid fraction and the phase concentrations. However these equations also contain other unknowns: the mass flux, , and the interfacial concentration, .

The phase change rate during primary solidification is principally controlled by species exchanges at the interface. The variation of the enthalpy and the interfacial solute balance provide relations connecting variables and .

〈 〉 𝑔 (148) (149) (150)
where : the specific heat and : the latent heat of fusion, which are assumed constant.

Note that the symbol is used to refer to variations within a micro time-step.

Additional steps of manipulation to establish the system of equations to be solved are summarized as follows.

The enthalpy variation on a micro-time step can be analyzed using Eq.( 148) as follows:

〈 〉 𝑔 ( 𝑔 𝑔 𝑔 ) ( ) ( 𝑔 𝑔 𝑔 ) (151) ⇒ 〈 〉 〈 〉 ( ) ( 𝑔 𝑔 𝑔 𝑔 ) (152) ⇒ ( ) 𝑔 𝑔 〈 〉 〈 〉 (153)
where is the number of micro-time steps per a macro time-step ; is the temperature at the previous micro time-step. In Eq.( 151) the interfacial liquid composition is assumed to follow equilibrium and is thus directly linked to the temperature by the phase diagram.

Note also that no transport terms remain in these equations. The main driving force for solidification is due to the enthalpy variation and the mass flux determined by manipulation of Eqs. ( 149) and (150) as follows, assuming a thermal equilibrium at the solid-liquid interface (i.e. the interfacial solid composition has been linked to that of the liquid phase as relation ).

〈 〉 〈 〉 (154) ⇔ ( ) 〈 〉 〈 〉 (155)
The two nonlinearly related variables and are determined by solving the system of equations ( 153) and ( 155) with an iterative approach, for which in the first term of Eq.( 155) is updated by the latest value. When doing this repetitive calculation process, as the phase compositions have not been determined yet, they are estimated by the last known values. If the solid formation rate is negative, , this indicates the occurrence of remelting for which the resolution procedure is re-effectuated considering that the interface and average solute concentrations on the solid side in the solute balances ( 146) and ( 154) are equal and unchanged relative to the previous iteration, i.e.

〈 〉 instead of considering . During the growth stage, the local grain density given by the solution of Eq.( 48) from the transport stage does not change, except in case of total remelting. The iterative procedure is repeated until these variables have converged.

The interests obtained from the nucleation and growth stage, including the solid fraction, the grain density, the intrinsic compositions and the temperature are integrated into the entire iterative scheme as shown in Fig. 3. 8. Fig. 3. 8 -Schematic of the resolution algorithm using the splitting method, introduced the solutions obtained from the nucleation and growth stage.

Coupling between Microscopic and Macroscopic Scales

Through the previous sections, we have studied the implemented model for the macroscopic and microscopic stages. Using a splitting method developed by Založnik and Combeau [Založnik and Combeau, 2010a], the coupling between these stages is conducted as schematized in Fig. 3. 9. Each iteration begins with an implicit FE resolution using a macro time-step for the energy and transport conservation equations in the transport stage neglecting the nucleation and growth terms. This gives the solutions with the subscript and further adding the subscript for those associated to the transport stage. Then the nucleation and growth stages are solved locally (at nodes), with variables initialized by values obtained from the transport stage at . The solution proceeds through micro time-steps, assuming that the average composition and the total enthalpy no longer evolve on a macro time-step. Finally, the momentum equations are solved with a semi-implicit solver on the macro timestep to compute the velocity fields at .

Fig. 3. 9 -Schematic of the resolution algorithm using the splitting method, with the growth and transport coupling.

Because of the close interaction between these physical phenomena, a strong connection is effectuated for these stages through the initialization procedure: final values obtained from one stage are used to initialize the other stage. For example considering a quantity , the variation of which results from transport, nucleation and growth processes we can express that its time derivative results from three contributions: Chapter 3 -Finite Element Implementation and Adaptation (156

)
where are the operators representing the advection, nucleation and growth respectively.

The splitting method is organized as follows on a macro time-step [ ] Initial value: Final value:

⏟ [ ] (157)
Initial value: Final value:

Initial value: Final value: (158)
where the subscripts , and 𝑔 represent values obtained from the transport, nucleation and growth stage respectively; the superscript indicates the time.

For clarity, nucleation is considered as a separate stage. While the variation of the grain density due to nucleation is important, this process only results in a slight change in amount of other variables. Therefore, nucleation can be joined with the growth stage to form only one stage representing for microscopic processes and the splitting scheme thus applies to two stages. In addition, a full evolution of quantities in each time step is evaluated by final results obtained from the growth stage since these solutions also include the change due to the transport stage. So in the case where there is no ambiguity between values of the transport and growth stages, the subscript 𝑔 is simplified and without subscript implies a value computed from the growth stage.

The resolution algorithm of the complete growth-transport model in presence of relevant equations and the flowchart of the nucleation-growth stage are presented in Fig. 3. 10 and

Extension to a Three-Phase Model

Through the previous study, we have worked on the two-phase model in which equiaxed grains are assumed to be fully globular. This is known as a two-phase model since only a spherical solid shell is surrounded by a liquid, as schematized in Fig. 3. 12-a). However, in most solidification cases, the morphology of the equiaxed grains is more complicated, like dendrites as the simplified illustration in Fig. 3. 12-b). A dendritic grain consists of the solid dendrite and its interdendritic liquid phase. This grain is separated from the extradendritic liquid phase by an envelope whose surface passes over the active tips of the primary and secondary dendrite arms [Combeau et al., 2009]. In order to account for the grain morphology, the present model is extended to the three-phase model which is derived from the work of Combeau and co-workers [Combeau et al., 2009] and Bedel [Bedel, 2014]. Details of the three-phase model are presented in these references and we summarize here the equations to be used and the main points which are additionally developed from the two-phase model.

As introduced in the literature review of solidification models, according to these references the distinction into three phases is considered at the microscopic scale and only two phases including the solid and the total liquid are to be present at the macroscopic scale. Therefore the interfacial stresses between the solid phase and the total liquid are modeled in the same way as in the two-phase model by using the relations ( 29)-( 36). However, it should be noticed that the characteristic length used in the three-phase model is the radius of a solid sphere that has the same volume as the solid dendrite within a grain. Additionally, in the two-phase model, the change of regimes occurs once the solid fraction has reached the packing limit. In the three-phase model, the interactions between the moving solid grains and the total liquid take place only when the envelop fraction, being higher than the solid fraction, attains the packing value, then the solid grains are blocked.

As existing the two phases at the macroscopic scale, the transport thus due to the movement of the two phases, consisting of the total liquid and solid phases, is almost unchanged with respective to that in the preceding two-phase model, except for adding a transport equation to be solved. In fact considering the two-phase model, the solid phase can be figured out from the solid fraction and the grain density. Therefore, its motion corresponds to the transport of the solid fraction and the grain density. While for the dendritic case, the grains are additionally characterized by the internal solid fraction or the envelope fraction. Consequently, modeling of the transport of dendritic grains requires the simultaneous description of transport of this extra quantity. The general resolution scheme is presented in Fig. 3. 13 where the parts highlighted in dark yellow are added for the three-phase model from the two-phase one. The detailed algorithm consisting of equations to be solved is shown in Fig. 3. 15.

Concerning the microscopic scale, the solute diffusion is modeled by idealizing a well-mixed solute in the interdendritic liquid, for which the 1D profile of concentrations can be presented as Fig. 3. 14-a) in consideration of the real configuration where a gradient of solute composition exists as schematized in Fig. 3. 14-b).

Similarly to the two-phase model, the evolution of variables in the growth stage is computed by integrating locally the conservation equations of mass and solute contents without the transport terms. Nevertheless the corresponding quantities in the liquid phase are separately evaluated for those in the interdendritic and extradendritic liquids, given as Eqs.( 159)-( 164).

(𝑔 ) (159)

(𝑔 ) (160) (𝑔 ) (161)
(𝑔 〈 〉 ) 〈 〉 (162) (𝑔 〈 〉 ) 〈 〉 〈 〉 (163) (𝑔 〈 〉 ) 〈 〉 ( 164
)
where the superscripts , and indicate the interdendritic liquid, the extradendritic liquid and the envelope respectively ; are the variation rates of the solid, the total liquid, the envelope in the interdendritic side and the envelope in the extradendritic side respectively; are the solute diffusion lengths in the solid dendrite, in the interdendritic liquid associated to the solid side, in the interdendritic liquid associated to the extradendritic liquid side and in the extradendritic liquid respectively; Moreover, the balances of total mass and solute mass flux exchanged at the interfaces between the solid and the interdendritic liquid and between the interdendritic and extradendritic liquids are described by Eqs.( 165) and ( 168).

(

) (166) ( ) 〈 〉 〈 〉 (167) 〈 〉 〈 〉 165
From the above equations ( 159)-( 168), by assuming a well-mixed solute diffusion in the interdendritic liquid (i.e. 〈 〉 ), as schematized in Fig. 3. 14-a), and neglecting the time variation of the concentration at the grain-extradendritic liquid interface (since this does not result in important errors as demonstrated in [Tveito et al., 2013]), we obtain the set of equations ( 169)-( 171).

(𝑔 〈 〉 ) 〈 〉 (169) (𝑔 〈 〉 ) (𝑔 〈 〉 ) (𝑔 〈 〉 ) 𝑔 𝑔 〈 〉 (170) ( ) 〈 〉 𝑔 𝑔 〈 〉 (171)
By regrouping the terms associated to the variables to be determined, and , on the LHS, Eq.( 171) is rewritten as Eq.( 172), assuming an interfacial equilibrium condition (i.e.

).

( 𝑔 𝑔 ) 〈 〉 𝑔 𝑔 〈 〉 (172)
Using the same approach as for the two-phase model, the system of equations to be solved, consisting of Eqs.( 169), ( 170) and ( 172), is closed by the enthalpy variation Eq.( 153). The solution of interests including the solid formation rate, , the interfacial liquid concentrations, and the intrinsic concentrations, 〈 〉 .

The solution algorithm is presented in Fig. 3. 16 where the terms in the yellow zones are the ones modified when comparing to those of the two-phase model. Applying an iterative procedure as presented in the two-phase model, first two equations ( 172) and ( 153) are solved to compute the solid formation rate, , the interfacial concentration, . Then the intrinsic compositions in the solid dendrite and in the total liquid are determined by solving Eqs.( 169) and ( 170).

Since the evolvement inside of grains is microscopically governed by solute exchanges at dendrite arms, the solid formation rate is thus calculated by solving mass balance equations at the solid -interdendritic liquid interface. While as being controlled by the thermal and solutal undercoolings, the grain growth is evaluated by the envelop variation rate or the envelop fraction 𝑔 (as denoted in the algorithm Fig. 3. 16) and modeled by the development of dendrite tips that can be assumed in a hemispherical or paraboloidal form. In our current work, the implemented formulations used to calculate the tip growth velocity with two different approximations are issued from the works of Kurz and co-workers [START_REF] Kurz | [END_REF]], which are summarized in the following.

For both cases, considering metal alloys, the undercoolings due to the effect of kinetics and curvature of dendrite tips can be neglected, so the total undercooling consists of the solutal and thermal parts, as given ( 173).

(173)

From the definition of the solutal and thermal supersaturations -the quantities being used to evaluate the gradients of solute and heat which are driving forces of the grain growth -, the corresponding undercoolings are expressed by the relations ( 174) and ( 175), as introduced in [START_REF] Kurz | [END_REF]]:

〈 〉 ( ( ( ) ) ) (174) (175)
where and are the solutal and thermal supersaturations respectively.

In a simple case, the shape of a dendrite tip is assumed to be hemispherical, for which the supersaturations are approximated by the Pléclet numbers (176), defined as the ratio of tip radius to the solute (or thermal) diffusion length. Under this assumption, the tip growth velocity can be computed as (177), issued from the Marie Bedel's thesis [Bedel, 2014].

(176)

[ ( 〈 〉 )] ( ) (177)
where the subscript indicates the solutal or thermal parts.

A closer morphology of a dendrite tip to its real form looks a likely paraboloid of revolution, as proposed by Papapetrou. With this consideration, the supersaturation being the resolution of the diffusion equation is found as a function of the Péclet number from the Ivantsov's work the so-called Ivantsov function [START_REF] Kurz | [END_REF].

( ) ( ) (178)
where the definition of the function of ( ) is referred to [START_REF] Kurz | [END_REF]]. Chapter 3 -Finite Element Implementation and Adaptation

The undercooling equation ( 173) accompanied with the constitutive relations ( 174), ( 175) and ( 178) allows determining the product of the tip growth velocity, , and the tip radius, . Then, the tip growth velocity, , is deduced once the radius tip being determined according to the relation developed in [START_REF] Kurz | [END_REF]], as follows:

[ ( ) ] (179)
where

( ) (180)
( ( ) ) (181)
Finally, for both cases, knowing the tip growth velocity, we can calculate the envelope fraction, 𝑔 , as presented in the algorithm Fig. 3. 16. Using this model, the internal solid fraction, being defined as the ratio of the solid fraction and the envelope fraction, which is computed from the growth stage, provides information on the grain morphology resulting from a growth competition between the solid inside the grain and the grain envelope. If the internal solid fraction is close to one, grains are considered as fully globular (or fully spherical). If it is smaller than one, the grain growth adopts the dendritic morphology.  Development of an adaptive artificial diffusion to deal with numerical issues encountered when using FE solvers for purely transport equations which are disturbed by discontinuities of divergence-non-free velocity fields.

 Implementation of the computational module for the solid velocity.

 Realization of necessary modifications in solving the energy equation and the momentum equation of the liquid phase.

 Implementation of the microscopic model accounting for a controlled solute diffusion in the liquid and solid phases.

 Achievement of the complete two-phase model.

 Implementation of the three-phase model for dendritic grains.

In the next chapter, numerical simulations will be performed and the results will be compared with references from the FVM simulation and the literature in order to verify our implementations. Chapter 3 -Finite Element Implementation and Adaptation

Résumé en francais

Dans ce chapitre, nous avons présenté un modèle de solidification dérivé du travail de Založnik et Combeau [Založnik et Combeau, 2010a] dans lequel les auteurs ont développé une méthode de splitting permettant de résoudre la complexité du couplage fort entre des phénomènes de transport et de germination-croissance. Outre une brève description du modèle existant, nous sommes rentrés dans les détails des adaptations effectuées pour mettre en oeuvre le modèle en utilisant la méthode des éléments finis. Les principes du modèle sont ainsi établis et une implémentation numérique est proposée. Nos contributions apportées aux modèles du CEMEF peuvent être résumées ainsi :

 Développement d'un terme de diffusion artificielle pour traiter des problèmes numériques lors de la résolution des équations hyperboliques par la méthode des éléments finis dans lesquelles les matières sont transportées par des champs vitesses dont d'une part le comportement présente des discontinuités et d'autre part la divergence est non nulle.

 Mise en oeuvre du module de calcul de la vitesse du solide.

 Réalisation des modifications nécessaires concernant la résolution de l'équation de l'énergie et de celle de la quantité de mouvement du liquide.

 Mise en oeuvre du modèle microscopique prenant en compte la diffusion du soluté finie dans les phases liquide et solide.

 Implémentation du modèle à deux phases pour des grains globulitiques.

 Implémentation du modèle à trois phases des grains dendritiques.

Ce modèle sera appliqué aux simulations numériques présentées dans le chapitre suivant dont les résultats seront comparés à ceux de cas de référence obtenus par la méthode des volumes finis et à ceux dans littérature afin de vérifier nos implémentations. For verification purposes the model is first tested for a case with pure heat diffusion, i.e. solving the energy equation and the nucleation and growth stage in the resolution algorithm for the two-phase model in Fig. 3. 10 (presented in Chapter 3). This case tests its capacity to account for the nucleation and grain growth processes at the microscopic scale. Both liquid convection and solid transport are neglected. Solidification occurs in a one-dimensional manner. This is achieved by defining a 1 mm thick and 100 mm wide domain cooled from one of its narrow sides with an exterior temperature at 100°C and a heat transfer coefficient of 500 W m -2 K -1 . The initial temperature of the melt is set at 750°C. Six sensors are defined along the sample, two at both ends and four between them with a regular spacing of 20 mm, as seen in Fig. 4. 1. The thermophysical properties of the material are taken from the reference [Ni and Beckermann, 1993] (reported in Annex -A.1) and are representative of a binary Al-4wt%Cu alloy, except for the solute diffusivity coefficients and the grain density. These two parameters are control parameters. They are used to switch between a case with a very high solute diffusion, so-called "infinite" diffusion, and the other case for a lower value of the diffusion coefficient, so-called "finite" diffusion. Results of these test cases obtained from the present model will be compared to those given by classical models including the lever rule and Gulliver-Scheil approximations [START_REF] Kurz | [END_REF]].

This test is calculated on a 2D mesh with the mesh size of 1 mm, the solution being made 1D through the control of boundary conditions. Simulation parameters are listed in Annex -A.1. First, for the simulation with high solute diffusion in the solid phase, the cooling curves and solid fraction evolutions obtained at the six sensors from the "infinite" solute diffusion case (green curves) are compared to those predicted by the lever rule (orange curves) as well as to those given by the "finite" diffusion case (red curves) in Fig. 4. 2.

Compared cases in Fig. 4 The same methodology is carried out for the test without solute diffusion in the solid phase, i.e. the solutions acquired from the "infinite" solute diffusion case (green curves) are evaluated by comparing to those given by the Gulliver-Scheil approximation (orange curves) and those obtained from the "finite" diffusion case (red curves) in Fig. 4 It can be seen in these figures that the green curves resulting from the present model exactly retrieve the lever rule for the case of high solid solute diffusion and the Gulliver-Scheil approximation in the study without solute back-diffusion. The cooling curves present a quick decrease of the temperature in the liquid zone due to the heat extraction and then a slow cooling which corresponds to the formation and evolution of the solid phase with the latent heat release. These simulations also illustrate the effect of the cooling rate as observed that the solid fraction increases quickly near the cold wall due to a high cooling rate and decreases gradually away from this side. With "finite" diffusion, it can be noticed that the first red curve on the cold side presents a recalescence close to the liquidus temperature, which is not computed when the liquid composition is well mixed, i.e. for the orange and green curves. This phenomenon, which characterizes an equiaxed solidification, can be captured by the present model since the nucleation and grain growth kinetics are taken into account. It depends on both the values of the grain density and the liquid diffusion coefficient and nucleation undercooling. When an embryo forms initially, the latent heat released is too small to prevent further cooling. When the solid grain growth is sufficient, i.e. for a higher undercooling, the latent heat released cannot be instantly compensated by heat extraction, which causes the temperature to increase. For the points farther away from the cold wall, cooling rate is decreased so recalescence disappears.

Yet the red curves are found to depart very little from the green curves. With high values of the Lewis number (the ratio between heat diffusivity and mass diffusivity) in metallic alloys, solidification at the microscopic scale is dominated by solute diffusion. A limited solute diffusion in the liquid phase results in a slow phase transformation at the beginning of solidification, illustrated by a difference in the solid fraction evolutions between the "infinite" and "finite" diffusion cases. This difference decreases as solidification progresses. During the formation of solid grains, the solute rejected from the solid to the liquid is accumulated at the grain boundary. When the process advances, the solutal gradient progressively increases around growing grain surfaces, leading to accelerate solute mixing.

Additionally, without solute diffusion in the solid phase, a formation of eutectic structure is present toward the end of solidification when the eutectic temperature is reached. This effect is reproduced by the present model is well correlated with the Gulliver-Scheil approximation, as shown in Fig. 4. 3. The amount of eutectic is proportional to the final vertical jump of the solid fraction curves to reach unity. Furthermore, as discussed elsewhere [START_REF] Thévoz | Modeling of equiaxed microstructure formation in casting[END_REF]], the choice of the time-step depends upon situations. It should be small enough in order to efficiently capture microscopic evolutions, e.g. the recalescence phenomenon, and must be sufficiently large to avoid certain numerical issues, as bouncing effects or false eutectic plateaus, that can be experienced when too small time-steps are employed in modeling of microstructures. Compared to the use of a single time-step scheme, the two-time-step technique is less time consuming because a reasonable time-step for macro transport phenomena can be retained.

With a good agreement between the reference solutions and those achieved from the present model and through the above analyses related to the effects of a partially mixed solute diffusion, it is demonstrated that the present model can be used as an alternative tool for classical approximations and additionally gives more flexibility in studying microscopic phenomena. Consequently, the verification of our implementation for the growth stage is qualified, from which the movement of the liquid phase driven by the thermo-solutal buoyancy is taken into account and the numerical simulations which are performed for this step will be presented in the next section.

Thermo-solutal Liquid Convection during Solidification

Test Cases Description

In order to verify our implementation when integrating natural liquid convection into the current model, two 2D solidification simulations have been carried out, which take into account buoyancy driven liquid convection, assuming a fixed solid phase.

+ Case 1 : The computational domain sized by 0.1 m x 0.06 m is cooled from the right and left sides, supposing a heat transfer coefficient of 400 W m -2 K -1 . An environmental temperature is 25°C and an initial liquid temperature is 285.49°C. Due to its symmetry, only haft of the cavity is simulated. Fig. 4. 4-a) shows the simulated geometry, the boundary and initial conditions. The sample material is a Pb-18wt.%Sn alloy, thermophysical properties of which are given in Annex -A.2. The temperature and concentration gradients act oppositely on the buoyancy forces: the liquid enriched in Sn becomes lighter and induces an upward flow along the solidification front, while the cooling flow driven by thermal effect goes downward. This benchmark was proposed as a reference test case for the verification and validation of different numerical solidification codes, presented in [START_REF] Bellet | [END_REF].

+ Case 2 : This case corresponds to an identical geometry 0.1 m x 0.06 m cooled from the left side only, a heat exchange coefficient of 300 W m -2 K -1 and the melt at an initial temperature of 260°C, as illustrated in Fig. 4. 4-b). The simulation is performed for the whole piece, considering a Sn-5wt.%Pb alloy that is characterized by cooperative thermal and solutal buoyancy forces. Related material properties are presented in Annex -A.2.

Although there is no solid motion, transports due to the solid movement are kept in the model (resolution algorithm in Fig. 3. 10) in order to verify the application of the numerical procedure to a general case, in which the solid velocity and the packing solid fraction are imposed equal to zero. Simulation parameters of Case 1 and Case 2 are reported in Annex -A.2.

Mesh form

Simulation Parameters Mesh size

[mm] Macro time step As we have just seen in the 1D simulations, the present model when using high values of solute diffusion coefficients, m 2 s -1 for the liquid phase and m 2 s -1 for the solid phase, produces results almost identical to those given by the lever rule. This tendency is also found in these 2D simulations using the FE formulation. In this section, we only present the comparison between the solutions obtained from the FE and FV formulations (FV results are obtained by using SOLID software). Similarities and differences will be pointed out. This comparison is intended as a based assessment for those upcoming in the subsequent parts including solid transport and three-phase model since the FE results will be evaluated by comparing with those of the FVM. Fig. 4. 5 depicts the final macrosegregation maps obtained from the present FE model with high values of the solute diffusion coefficients and those predicted with the lever rule approximation using the FVM. It can be seen that a global agreement is achieved between these results. A high content of Sn is predicted at the top of the cavity close to the symmetry plane for Case 1 since during solidification the light solute-rich liquid rises up and two flows coincide at the center of the cavity due to its symmetrical configuration. As a result, the solute is continuously accumulated in the upper zone as well as in the middle area and a solutedepleted distribution is formed in the lower part of the domain. Because of the flow dominated by the heavy solutal element, a different macrosegregation pattern is obtained in Case 2: a strong positive segregation occupies the bottom and a large zone on the right hand side while a solute-poor region is located in the upper part of the domain.

There exist differences between the results obtained with FEM and FVM that can be attributed to the influence of numerical factors, i.e. discretization schemes, interpolation functions, iteration procedures, etc. Similar differences have already been pointed out and discussed in [Combeau et al., 2012b] where several numerical models have been benchmarked. We reported in Fig. 4. 6 the results obtained from different methods, issued from this reference. The same macrosegregation pattern is globally predicted. The solutions differ from each other in local variations (the number of mesosegregation channels, their intensity and inclination, etc). One of the main reasons is supposed to be related to solving the momentum conservation equation for which the permeability is differently discretized and interpolated upon modeling schemes. It has a strong influence on the prediction of mesosegregated channels. It is worth noting that for the Pb-18wt.%Sn alloy the segregated channels simulated with an unstructured FEM mesh extend over a larger zone compared to the one obtained by using a structured FVM mesh. The opposite tendency is observed for the Sn-5wt%Pb alloy. This phenomenon could be explained by the fact that the unstructured FEM mesh is more sensible to perturbation due to liquid movement. While the counteracting thermal-solutal advection of Pb-18wt%Sn alloy results in instabilities, the cooperating convection of Sn-5wt%Pb might be less unstable. Furthermore, the concentration profiles at the three height levels of the cavity (10 mm, 30 mm, 50 mm) obtained from the present model are in good agreement with those predicted in the reference [Combeau et al., 2012b], plotted in Fig. 4. 7.

FEM (present approach) FVM (SOLID software)

Case 1

Case 2 Fig. 4. 5 -Average composition maps obtained (left) from the present FE model using high values of solute diffusion coefficients compared with those (right, issued from [Combeau et al., 2012] and [Založnik and Combeau, 2010b]) from the lever rule with FVM.

FVM (FLUENT software) FVM (THETIS software) FVM (OpenFOAM software -upwind) FVM (OpenFOAM software -QUICK) Fig. 4. 6 -Average composition maps obtained from different numerical solutions for the 2D benchmark of binary columnar solidification (using a Pb-18wt%Sn alloy), described in [START_REF] Bellet | [END_REF]. The results are issued from [Combeau et al., 2012b]. As studied in the 1D simulations, the limited solute diffusion is stated having effects on the evolution of solidification. Therefore, its influence will be also analyzed in the 2D context and compared with the case of infinite solute diffusion, as presented hereafter. This topic was investigated for a case study on a Sn-5wt%Pb alloy in the FV context by Založnik and Combeau [Založnik and Combeau, 2009b].

Limited Solute Diffusion at the Microscopic Scale

Two simulations for the two alloys are performed, which are identical to the above cases, except for the solute diffusion in the liquid phase and the grain density decreasing to m 2 s -1 and grains m -3 , respectively (instead of m 2 s -1 and grains m -3 in the previous section). Fig. 4. 8 presents the final macrosegregation maps obtained from the present model for the two alloys. These maps are globally similar to those of the infinite solute diffusion presented in the previous section. However, we observe that mesosegregates are not formed in the case of low solute diffusion. This tendency was also found in the work of Založnik and Combeau [Založnik and Combeau, 2009b], which further allows verifying our implementation. Based on analyses in this reference, we will present in this part a study related to the impact of the limited solute diffusion on the formation of mesosegregation channels.

Test case with a Pb-18wt%Sn alloy

Test case with a Sn-5wt%Pb alloy In fact, for the case of fast solute diffusion in the liquid at the microscopic scale, studied on a binary alloy, the variation of the temperature and of the intrinsic composition in the liquid phase are in proportion given by the relations ( 182) and ( 183). So we can see that a perturbation of the liquid concentration 〈 〉 results in an immediate adaptation of the solidification kinetics following the relation (184) deduced from ( 182) and ( 183). In the zones where 〈 〉 is positive, the enrichment in solute decelerates solidification. The longer the liquid rests, the more the mushy zone is perturbed by the fluid flow so that the segregated channels have convenient conditions to form and become later the favorite paths for the fluid flow. Consequently, they continue to develop and finally form mesosegregated channels.

𝑔

(182)

〈 〉 (183) 𝑔 〈 〉 ⁄ (184)
Conversely, in the low solute diffusion case, the adaptation of liquid fraction to a perturbation of the solute concentration, 〈 〉 , can be deduced from relation (185). By neglecting a back diffusion, we have the relation ( 186) which shows that the response of 𝑔 to the variation of composition in the liquid is not immediate, like in the previous case of high diffusion. Under such conditions, there is less perturbation by the liquid flow, the mushy zone is more stable, so the finite diffusion in the liquid phase induces certain constraints on the phase transformation and thus restrains the formation of mesosegregated channels.

∑ ( ( 〈 〉 ) ) (185) 𝑔 ∫ ( ) 〈 〉 (186)
Additionally, for the symmetric simulation with the Pb-18wt.%Sn alloy, the positive segregation is still found at the top of the cavity but does not expand on the center line as much as it does when the solute is assumed to be well diffused. This difference is also suggested due to the more steadiness of the semi-solidified zone in the case where the diffusion is controlled, so the less quantity of the solute enriched liquid penetrates into the mushy zone.

Through this section, it is shown that the present model reproduces the reference solutions of 2D simulations accounting for the thermo-solutal liquid convection. The comparisons were reported between the FE and FV solutions with a global agreement besides differences due to numerical aspects. Moreover, apart from its advantages illustrated in the above 1D simulations, this model, via analyses on the 2D test applications, provides a useful tool to study the influence of the finite solute diffusion on segregation, which is the case in the reality while most models consider a perfect solute diffusion in the liquid phase with also a perfect solute diffusion in the solid phase or no diffusion in the solid phase (lever rule and Gulliver-Scheil approximations). Unlike the infinite solute mixing, the partial diffusion likely delays the transient solidification process and limits the formation of mesosegregates.

These studies showed the ability of the present model to simulate solidification considering the fixed solid as well as allowed verifying our implementation for this step. In the next section, the integration of the solid motion will be examined and an application to simulate the sedimentation of solid grains will be carried out.

Purely Convective Transport during Sedimentation

Test Case Description

Simulation of sedimentation first neglects the nucleation and growth stages, i.e. with a given number of grains being assumed to be globular with constant size. The gravitational force is the only body force that acts on both phases. The velocities at the top and bottom sides as well as the normal velocities at the vertical walls are set to zero and a perfect slip condition is applied to the tangential velocities at the vertical walls, as presented in Tab.4. 2 and Fig. 4. 9.

Fig. 4. 9 displays the sample geometry and the initial state of variables, companied with thermal and mechanical boundary conditions in Tab.4. 2. At the beginning, the grains are distributed with a fixed density, 10 9 grains m -3 , within the full width of a subdomain from height 20 mm to 80 mm in a 10 mm × 100 mm domain. The volume fraction of solid in this subdomain is set to 0.1, thus fixing the grain radius to 288 µm. The intrinsic solute compositions and the temperature are related to the solid fraction through the lever rule. These quantities are disposed in order that there is a homogeneous average composition at 5 wt.% Pb on the entire sample: i.e. supposing the initial temperature to be uniform at 498.057 K, the intrinsic composition in the liquid phase is 5.515 wt.% Pb where there exists the solid and 5 wt.% Pb elsewhere; thus the intrinsic composition of the solid is 0.328 wt% Pb where the solid is present and zero elsewhere. Lastly, the enthalpy is evaluated according to the solid fraction and the temperature. The simulation is operated on an unstructured triangular mesh with a size of 0.5 mm, a constant time step 0.01 s and 𝛼 and for the two constant parameters of the diffusion artificial coefficient defined in (45) (in Chapter 3 -Section 3.1.1). The simulation parameters and the thermophysical properties of the binary alloy Sn-5wt%Pb (the same alloy used in Section 4.2) are reported in Annex -A.3.

Tab.4. 2 -Thermal and mechanical boundary conditions

Analyses and Coherency Verification

The sedimentation process is illustrated through the maps of solid volume fraction distribution and its profiles at t = 0 s, 2 s (before the grains have reached the bottom of the domain), 5 s and 10 s (accumulation below the packing limit), and 15 s (end of packing), as shown in Fig. 4. 10.

From the initial state, the solid grains fall with a velocity of around 6 mm s -1 and the smaller liquid velocity has the opposite direction so that the continuity equation is respected, as displayed in Fig. 4. 11. The first grains arrive at the sample bottom after 3 s and begin to accumulate until the characteristic packing fraction is reached, taken as 0.3. When the process advances, we can observe a deviation of the liquid velocity from the vertical direction due to the interactions with the settling grains, especially an emergence of a large vortex is observed at 5 s. The sedimentation is complete after 15 s, although the solid and liquid velocities still rest within a small layer where the gravitational force is stably balanced with the diffusive effect, following Eq.( 187) which is derived from Eq.( 47) for the steady regime. This phenomenon is maintained due to the persistent gradient of the solid fraction at the transition surface.

(𝑔 〈 〉 ) ( 𝑔 ) (187)
As the grain motion involves the transport of other quantities, it is important to verify the consistency of these transports. Compared to the evolution of the solid fraction, an identical behavior of the grain density as well as the same profiles of the enthalpy are shown in Fig. 4. 12 and Fig. 4. 13 (top-left), which demonstrate a compatible transport exerted on these quantities during the process. In addition, the profiles of the average composition in Fig. 4. 13 (top-right) also present its accordance with the others. When solid grains sediment, the soluterich liquid moves upward 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.6 wt.% Pb in the negative segregation zone at the bottom corresponds to the final state of about 30% of solid at 0.362 wt.%Pb and 70% of liquid at 5 wt.% Pb. Furthermore, as expected for a pure transport simulation, the temperature almost does not change during this process although several oscillations of its value can be seen at the transition zones, but kept smaller than 1°C. This is displayed in Fig. 4. 13 (bottom-left). A good overall conservation of all quantities is also verified, as shown in Fig. 4. 13 (bottom-right), the maximum relative errors of solute mass and energy being about 10 -6 and 10 -5 , respectively. Compared to the evolution of the solid fraction given from the FVM (SOLID software), shown in Fig. 4. 14, the sedimentation simulated by the present FEM model has globally similar tendency from the beginning until the completion after 15 s. However, the results obtained from our FEM model are likely more diffused and progress faster than those predicted by SOLID regarding the intermediate solutions at 5 s and 10 s. Moreover, the solid layer should be accumulated at the packing fraction 0.3 but it is observed that the maximum solid fraction slightly exceeds this predefined value, reaching 0.31 in Fig. 4. 10. Although the treatments for the solid velocity, as presented in the previous chapter, are carried out, which comprise an addition of the artificial diffusion to purely advective equations as well as an adjustment to avoid the penetration of moving grains into the packed region. This excess of the solid fraction demonstrates that the solid accumulation is not well controlled at the packing limit, or in other words, singularities are not sufficiently reduced. As a consequence, it would be necessary to evaluate the sensibility of the model to the amount of artificial diffusion, which will be conducted in the following.

Solid fraction, [-] 0 s 2 s 5 s 10 s 15 s

Height -Solid fraction

Simulation Parameters Mesh size (structured mesh)

[mm] Time step

[s] Fig. 4. 14 -Results obtained from the FVM (SOLID software) at 0 s, 2 s, 5 s, 10 s and 15 s with (maps) the solid fraction distributions and (graph) its profiles along the sample's height.

Effects of Artificial Diffusion

At first, two simulations were performed in order to study the impact of the artificial coefficient by separately evaluating the effect of its constant parameters, recalling the formulation expressed as

( ) (𝛼 | 〈 〉| | 〈 〉 |).
Case 1: 𝛼 ;

Case 2: 𝛼 ; Fig. 4. 15-top shows the profiles along the sample's height of the solid fraction and the temperature for Case 1. The simulation is stopped after 2.7s because of numerical problems occurring that are also manifested by a large variation of the temperature before this time.

Although also exposing remarkable oscillations of the temperature but with a smaller order than Case 1, Case 2 can overcome problems encountered in Case 1. This demonstrates that the term related to the variation of the intrinsic velocity has a more important role in solving numerical singularities than its counterpart. Additionally, it can be seen that the solid packing not only exceeds the maximum fraction, which has already been stated for the case where 𝛼 ;

, but also has a heterogeneous distribution that further reveals instabilities in numerical resolutions.

For this reason, a next computation is carried out by increasing the artificial diffusion with 𝛼 ;

. Fig. 4. 16-left shows the profiles of the solid fraction distributed along the height of the sample at different times. It can be seen that the solid layer is accumulated at the fraction which is almost exactly the same as the predefined packing value. However, the transition zone becomes more extended due to diffusion compared with Fig. 4. 14. This set of simulations confirms that numerical concerns due to discontinuities of the transport velocity field and an absence of diffusion in hyperbolic equations can be solved by adding a supplementary diffusive component. Nevertheless, taking a high value of , for instance, risks generating too much diffusion. Therefore, it is important to regulate the amount of additional diffusion to avoid unreasonably diffusing quantities.

After testing on these parameters, we found that and 𝛼 being approximated from 20 to 70 can provide a compromise for issues related to diffusion and instabilities. As plotted in Fig. 4. 16-right the profiles of the solid fraction where 𝛼 ; , which become more homogeneous with a smaller excess, around 2% relative to the packing value reference, compared to those in Fig. 4. 15, for the case 𝛼 ; , excessing more than 10%. Following the integration of solid motion for the transport stage, the combination of transport stage and growth stage is realized to achieve a complete two-phase model. Simulations in the presence of solid transport due to sedimentation and natural liquid convection are performed using the same Sn-5wt%Pb alloy. The packing solid fraction is also assigned at the constant of 0.3 as the previous simulation. Studying the effect of this parameter on solidification can be found elsewhere [Ilegbusi and Mat, 1998][Arnberg et al., 1993][Bedel, 2014]. Based on experimental observations and comparisons, low values of the packing fraction, from 0 to 0.3, are recommended by the work of Krane [Krane, 2004]. The initial and boundary conditions are figured out as Fig. 4. 17 and Tab.4. 3; simulation parameters are listed in Annex -A.4.

Height

Mesh form

Simulation Parameters Mesh size -c, -d, -e, -f present the distributions of the solid fraction, the grain density, the temperature, the average composition, the liquid velocity and the solid velocity, respectively at t = 10s. Considering these subfigures, we can see that cooling starts on the LHS. As soon as solid grains form, they settle to the bottom due to the gravitational force and the downward thermo-solutal convection. Some of them are still small and are transported toward the RHS. They are also carried by the liquid and move upward to a certain level.

In details, Fig. 4. 18-a displays the solid fraction at t = 10 s, in which the three black isolines represent the solid fraction at 0.1 (upper isoline), 0.2 (intermediate isoline) and 0.3 (lower isoline). Due to the settling and transport of solid grains, the solid layer has reached the packing state and then being accumulated along the bottom side where the highest solid fraction is found. The coolest zone is yet seen on the LHS of the cavity, as illustrated by the temperature map in Fig. 4. 18-c, showing that the solid fraction distribution is indirectly related to the temperature. Additionally, it can be seen in Fig. 4. 18-b that besides a large number of grains gathered in the lower regions by sedimentation, the transport of crystals due to the liquid advection, as shown in Fig. 4. 18-e and -f, results in a rather high grain density in several zones of the slurry region, yet not reflected by the sold fraction map.

Furthermore, the average composition in Fig. 4. 18-d shows that while a large part of the cavity still remains at the initial content, we can find a negative segregation at bottom and solute-rich liquid just above this zone. During the solidification, the rejection of solute from the solid to the liquid enriches the liquid and makes it heavier, inducing a downward liquid convection. However the bottom area is occupied by the solid phase composed of solutedepleted grains, even heavier than the solute-rich liquid, as plotted in Annex -A.4.

The evolution of the quantities at t = 20 s is shown in Fig. 4. 19. The solid layer has progressively advanced from the bottom side and negative segregation has extended in the same manner. At this moment, we can see in Fig. 4. 19-e and -f that the liquid and solid velocities become disturbed; the counter-clockwise rotation of the liquid convection observed at the early state has been broken due to effects of the sedimentation. This disturbance thus influences the transport of grains and results in a heterogeneous distribution, besides the causes considered in the previous study at t = 10 s. In addition, it can be seen that there exist instabilities of the temperature at the packing surface. This could be due to numerical effects as noticed and discussed in the 1D simulation section, that are related to discontinuities of the velocity field when the solid packing takes place.

Fig. 4. 20 presents the state of the quantities at t = 200 s, which are distributed in two distinct areas: slurry and packed areas. While the solid fraction is still small in the upper part, much smaller than 0.1, it is fairly high in the lower one, reaching around the solid fraction of 0.8 near to the cooled wall. The layer of the solid fraction from 0.1 to 0.3 remains rather thin in comparison with the domain's size. This indicates a quite prompt transition from one to another region. In the packed stationary bed, the solid velocity is zero but there exists an intergranular downward flow (displayed by the white streamlines) which brings the solute from the upper to the lower region. Although the large part in the stationary mushy zone has a negative segregation due to the accumulation of solute-depleted grains, at the very bottom of the cavity, the intensity of the flow through the mushy zone is sufficient to lead to remelting and finally forms a positive segregation with the form of a horizontal channel. It is also noticed that the counterclockwise liquid velocity is less intense, its maximum magnitude approaching 8 mm/s while it is around 30 mm/s for the moving solid case, because it is only driven by the buoyancy force without reinforcement due to the settling of the solid grains. In addition by comparing the results at 200 s in Fig. 4. 20 and Fig. 4. 24, it can be seen that solidification progresses faster when accounting for the solid transport than in other case. This can be explained by the fact that solid grains that settle in the vicinity of the cold mold are replaced by a hotter liquid phase. Consequently, the global heat exchange between the metal and the mold takes place more efficiently and results in faster solidification than the case of no solid movement where the heat is extracted by conduction through the solid layer fixed at the region near the cooled wall; the isolines of the solid fraction and of the isotherms are fairly vertical and superimposed. The final distributions of solute and grains for the fixed solid case are exposed in Fig. 4. 25 displayed with the same scales used in the solid transport simulation for a comparison. The large difference between these two cases demonstrates a significant effect of the solid transport on macrosegregation. Furthermore, although a higher solute content is found at the bottom area for both cases, its formation mechanism is quite different. In the case without solid movement, it is primarily governed by natural convection in the zone with relatively high liquid fraction. The solid zone extends from left to right and the rejected solute is progressively accumulated in the bottom-right region due to higher density of the soluteenriched liquid. A high solute concentration is therefore found at the bottom-right corner. With grain transport, a part of the initial negative segregation in the lower pattern is replaced by the positive one which is due to the circulation of the liquid phase in the mushy zone and the remelting, leading to an expense of the solute content at the top of the cavity. Additionally, the heterogeneous distribution of grains is totally confined to the solid movement.

at the end of the solidification

Average composition 〈 〉 [wt.%]

Grain density [grains m -3 ] Fig. 4. 25 -Case with fixed solid phase -Distributions of (left) the average composition (min value: 3.25 ; max value: 18.88) and (right) the grain density.

Effects of Artificial Diffusion

For this part considering the constant parameter as in the previous study, simulations are performed with two different values of 𝛼 and compared to the previous one where 𝛼 as well as to those obtained using the FVM from SOLID software.

Fig. 4. 27 presents the distribution of the solid fraction and the grain density at t = 10 s.

Overall it can be seen that the smaller the artificial diffusion coefficient, the less the distribution of solid grains expands and the closer the distance between isopleths of the solid fraction.

Intermediate results at t = 100 s in Fig. 4. 28 illustrate that the diffusion coefficient strongly influences on the inclination angle of the slurry-porous interface: the larger the diffusion is applied, the less the interfacial slope is observed. It can be explained by the fact that besides a small quantity of the settling solid after nucleation on the cooled wall, the solid zone in the left part is built up by solid grains that almost grow during a long travel: starting from nucleation on the left side, being advected to the right wall and carried along by the slurry liquid from right to left, then they settle when coming back to the left region. Thereby this phenomenon takes place consecutively during the process. In the case of high diffusion, grains can diffuse more easily on the packed bed surface before accumulation. This effect prevents grains from a quick packing after sedimentation, thus the piling up is spread more horizontally than in the situation with small diffusion.

Altogether, it can be observed in Fig. 4. 29 that the final maps of macrosegregation and grain density obtained from using these three coefficients 𝛼 are similar. As anticipated, those given by 𝛼 are smoother than the others and the simulation results with 𝛼 and 20 are closer to those given by FVM. However, while the energy and the mass are well conserved in the case where 𝛼 , as seen in Fig. 4. 30, these quantities reveal a slight deviation from the conservation curves with 𝛼 . An additional test with 𝛼 cannot be carried out until the end of solidification because problems of coherency occur that are provoked by numerical instabilities. Moreover, the most remarkable difference between the FEM and FVM solutions can be realized on the distribution of the grain density, as seen in Fig. 4. 29. Compared to the FVM results, a larger grain density situated on the RHS of the cavity is found by the FEM simulation. As the number of grains is governed by the nucleation and then influenced by transport, remelting and re-nucleation mechanism, a large quantity of grains generated may be due to a high frequency of crystals emptiness because of the transport or remelting phenomenon, thus the re-nucleation would happen more often. However, yet it has been not known clearly which reason leads to this state. Fig. 4. 30 -Time evolutions of (red curve) the energy in a conserved case, (green curve) the energy and (blue curve) the average mass calculated from simulations.

Through these simulations, it has been shown that introducing an artificial diffusion is necessary and 𝛼 can be considered as a lower limit of the proposed formulation. However this value is only partially concluded, it may be not favorable for all situations and can be tested in order to estimate a suitable quantity. Therefore a further development of a more general artificial diffusion coefficient would be expected to optimize numerical resolutions.

Dendritic Solidification Modeling

After studying a two-phase model intended for globular solidification, this section is contributed to the extended model for three phases which allows taking into account the morphology of dendritic grains. The simulations performed are principally used to verify our implementation by comparing the results to those obtained from the FVM (SOLID software).

Purely Diffusive Solidification

The first verification is realized for the implementation of the growth stage by simulating a pure thermal solidification where a sample of Al-5wt.%Si is initially set at 621.5°C and cooled with a heat extraction rate equal to 42300 W kg -1 . The test cases come from those figured out in the Ph.D. thesis of Marie Bedel [Bedel, 2014] for three final grain sizes: 0.1 mm, 1 mm and 10 mm, to which the nucleated grain densities correspond. With each grain size, three simulations are carried out for different growth models: one for the two-phase model, two for the three-phase model where dendritic tips are assumed to be either hemispherical or paraboloidal. Thermophysical data and simulation parameters are presented in Annex -A.5. Fig. 4. 31 shows the temporal evolution of the temperature for these simulations, the graphs on the left display the results obtained from our implementation and those on the right present the solutions from the thesis of Marie Bedel. The fractions of the solid phase, the envelope and the internal solid are plotted versus time in Fig. 4. 32-Fig. 4. 34 for the three cases. It can be seen that the present implementation produces the results which retrieve those of the reference, which allows validating the step.

Analyses for these results are detailed in the study of Marie Bedel, here we summarize the main remarks of her study related to the three-phase model.

 Overall, the three-phase model predicts a faster recalescence with respect to the twophase model because the solid phase increases more rapidly due to a higher gradient of solute.

 Additionally, when the grain density is small (i.e. when the final grain size is large), solid grains develop rather with a dendritic form. Therefore the three-phase model would allow a better prediction of the growth evolution than the two-phase model in which grains are supposed to have a full globular morphology.

 For a three-phase configuration, the modeling of the morphology and the recalescence is significantly affected by the assumption made for the velocity of dendritic tips. With a given undercooling, the tips velocity computed using the hemispherical model is higher than calculated with the paraboloidal model. On the one hand it results in a faster expansion of envelopes and the grains become thus more dendritic; on the other hand it induces a quick development of the solid fraction leading to a quick recalescence.

FEM results

FVM results

Final grain size : 0.1 mm

Final grain size : 1 mm

Final grain size : 10 mm Furthermore, the 1D test case which was established in Section 4.1 for the two-phase model is also performed here to verify our implementation with respect to the FEM solution of the pure heat diffusion without the solid and liquid movement. Results obtained from the present threephase model (assuming either hemispherical tips or paraboloidal tips) are compared to those given by the lever rule and Gulliver-Scheil approximations, for which the time evolutions of the temperature and of the solid fraction are presented in Through these results, it can be seen that the three-phase model, either with hemispherical tips or with paraboloidal tips, reproduces the lever rule and Gulliver-Scheil solutions. This study additionally confirms and validates our implementation.

Complete Growth-Transport Solidification

After studying the growth stage, the complete transport-growth dendritic model is applied to simulate a 2D case that was figured out in Marie Bedel's thesis [Bedel, 2014]. A 76 mm square cavity contains a Al-22wt.%Cu alloy in the liquid state initially at 700 °C, and is cooled from the left side in an environment at 500 °C with a heat transfer coefficient of 14000 W m -2 K -1 , as schematized in Fig. 4. 37. Thermophysical properties of the material and simulation parameters are presented in Annex -A.6. Two simulations of the three-phase model which assume hemispherical tips (Case 1.a) and paraboloidal tips (Case 1.b) are performed and compared to the two-phase case (Case 2). The computation is operated on a non-structured triangular mesh size of 1 mm and a macro-time step of 0.01 s.

Mesh form

Simulation Parameters Mesh size

[mm] . 38 presents the distribution of variables at t = 40 s for the three simulations, including from top to bottom the temperature, the solid fraction, the internal solid fraction and the grains density. Overall, while the solutions given by the two dendritic models are different, those obtained assuming paraboloidal tips are very close to those predicted by the two-phase globular model. It can be seen that the internal solid fraction calculated from the paraboloidaltip model has a high value close to one for nearly the whole domain, as shown in Fig. 4. 38 (the third row). As a result, the evolution of variables simulated by using the paraboloidal-tip model is similar to that described by the two-phase globular model.

Macro time step [s] Macro/Micro time [-] ( ) [-] ( ) [-]
Since the calculation given from the hemispherical-tip model predicts a rapid growth of envelopes, the grains with hemispherical tips possess a morphology more dendritic than paraboloidal-tip dendrites, as deduced from the internal solid fraction maps in Fig. 4. 38 (the third row). On the other hand, it also leads to a quick evolvement of the solid phase and thus a large quantity of solid is predicted to form on the cooled wall, being observed in Fig. 4. 38 (the second row). As consequence, the solid layer which attaches to the cooled wall makes the heat exchange between the hot liquid and the environment less efficient; the cooling takes place more slowly compared to the other cases, as presented in Fig. 4. 38 (the first row).

Regarding the grain density maps in Fig. 4. 38 (the fourth row), since the density of the solid phase is smaller than the liquid density, solid grains thus have a tendency to float and accumulate in the upper zone, unlike the previous simulations in which solid grains are heavier than the melt bulk. In addition, during solidification the liquid phase has continuously increased its density due to the effect of thermo-solutal buoyancy. Therefore finally a higher grain density is found at the top of the cavity and a positive segregation is located in the bottom-right area where the last solidification takes place, as shown in Fig. 4. 39.

Moreover, since hemispherical grains grow faster, the packed bed zone is formed earlier, leading to a weaker segregation, compared to those obtained from the paraboloidal-tip model and the two-phase model. This is observed in the maps of macrosegregation of Fig. 4. 39 and being confirmed when plotting the time evolution of segregation indexes, as displayed in Fig. 4. 40, the lowest is given by the hemispherical-tip model. However, the cooling and solidification are faster at the upper zone in the FE simulation. It is also noticed that the grain density given from the FEM is higher than that from the FVM. As a high grain density leads to a small undercooling and thus reduces the tip growth velocity, the development of envelopes will be limited. For this reason, although the solid fraction maps are quite similar, the internal solid fraction is rather different, which is considerably higher with the FE solution. Differences of the grain density predicted by the FEM and FVM formulations were also remarked in Section 4.4 for the two-phase globular model, for which the exact origin has not been known yet.

Summary

This chapter presented different simulations in order to verify our implementation for the twophase globular model and the three-phase dendritic model in the FEM context.

First, the present model was validated for a case of pure heat conduction, by comparing to the lever rule and Gulliver-Scheil approximations. It was demonstrated that the current model can retrieve the solutions of these classical models by using appropriate parameters consisting of the solute diffusion coefficients and grain density. In addition, this model can predict the recalescence phenomenon thanks to accounting for grain growth kinetics, which cannot be captured by the simplified models.

Then, validation test cases accounting for the liquid natural convection were realized on the Pb-18wt%Sn and Sn-5wt%Pb benchmarks. The results presented a good agreement between our current FE formulation and other different approaches. One of the advantages of the present model in solidification modeling with the controlled solute diffusion was also considered and analyzed, showing its impact on the formation of meso-segregated channels.

The following part was investigated to study the model with the integration of the solid transport. The pure solid transport was applied to simulate the sedimentation phenomenon, in which a good coherency between the variables during the transport process was verified. The FE solutions were further compared to those using the FVM, presenting a close agreement between these results, which allowed verifying our implementation. Moreover, the effect of the newly developed parameters in the FE model was studied and analyzed, illustrating that numerical issues related to FE solutions for hyperbolic equations can be overcome by using a suitable adaptive artificial diffusion.

After that, the complete model coupling the solid transport with other phenomena was performed for the reference test case which was realized by Založnik and Combeau, using the same model with the FVM [Založnik and Combeau, 2010b]. Finally, the application of the extended model -the three-phase dendritic model -was studied in the similar manner for the two-phase model, realizing different simulations from the pure thermal conduction to the full transport-growth coupling. In general, a correspondence between solutions obtained by the FEM and FVM formulations was achieved. Differences between these results as well as the sensibility to the artificial diffusion coefficient of the FE model were pointed out and analyzed. Through our supplementary test cases and studies in literature, the solutions were perceived to partially depend on numerical factors, such as time step, mesh size, discretization schemes, artificial diffusion … However, there exist other aspects that have been not thoroughly identified, for which we suggested several possibilities which may induce these differences.

Résumé en français

Dans ce chapitre, nous avons présenté des différents cas test de simulation dans le but de vérifier notre implémentation pour le modèle globulitique à deux phases et pour le modèle dendritique à trois phases, dans le contexte de la formulation des éléments finis.

Dans un premier temps, le modèle actuel a été validé sur un cas de conduction thermique pure, par une comparaison avec la loi des leviers et l'approximation de Gulliver-Scheil. En utilisant des paramètres appropriés incluant les coefficients de diffusion du soluté et la densité de grains, le présent modèle permet de retrouver les solutions des modèles classiques. De plus, ce modèle peut prédire le phénomène de recalescence grâce à la prise en compte des cinétiques de croissance de grains, ce qui n'est pas capturé par les modèles simplifiés.

Ensuite, d'autres cas test de validation en présence de la convection naturelle du liquide ont été réalisés sur les benchmarks de Pb-18wt%Sn et de Sn-5wt%Pb, présentés dans [START_REF] Bellet | [END_REF] et [Založnik and Combeau, 2010b]. Les résultats ont montré un bon accord avec ceux obtenus par les autres approches. L'un des avantages du présent modèle, la capacité de prendre en compte une diffusion du soluté limitée, a été considéré et analysé. Les cas de simulation correspondants ont illustré l'impact de ce phénomène sur la formation des canaux de ségrégation.

La partie suivante a été consacrée à étudier le modèle lors de l'intégration du transport du solide. Tout d'abord, le transport du solide pur a été appliqué à la simulation du phénomène de sédimentation, dans laquelle une bonne cohérence entre les quantités pendant le transport a été vérifiée. Ces solutions obtenues par la méthode des éléments finis ont aussi présenté une concordance avec celles que donne la méthode des volumes finis. Via ce cas test, l'effet des nouveaux paramètres développés dans la formation des éléments finis a été également étudié. Cela a démontré que des problèmes numériques concernant les solutions des éléments finis pour des équations hyperboliques pouvaient être surmontés en utilisant une diffusion artificielle convenablement adaptée.

Dans une dernière étape, le modèle complet dans lequel le transport du solide est mis en couplage avec les autres phénomènes a été appliqué au cas de référence qui a été réalisé par le travail de Založnik et Combeau en utilisant le même modèle avec la méthode des volumes finis [Založnik and Combeau, 2010b]. Puis, l'application du modèle étendule modèle dendritique à trois phasesa été menée de même manière pour le modèle à deux phases, dans laquelle différentes simulations ont été mises en place, du cas thermique pur à la simulation avec un couplage complet transport-croissance. En général, une cohérence entre les solutions données par les méthodes des éléments finis et des volumes finis a été obtenue. Des différences entre les résultats ainsi que la sensibilité du modèle des éléments finis au coefficient de diffusion artificielle ont été remarquées et étudiées. Via nos cas test supplémentaires et des études dans la littérature, les solutions se sont avérées dépendre partiellement des facteurs numériques, par exemple le pas de temps, le maillage, les schémas des discrétisation, la diffusion artificielle… Toutefois, il y aurait d'autres aspects dont leurs influences ne sont pas encore clairement identifiées, pour lesquels nous avons proposé des possibilités qui pourraient induire ces différences.

Chapter 5

Tests 3D et Applications Industrielles Suite aux études pour une configuration 2D du chapitre précédent, le modèle est tout d'abord appliqué à des cas 3D du benchmark de Hebditch-Hunt [Hebditch and Hunt, 1974]. Des applications industrielles sont ensuite menées sur un lingot d'acier 3t3 coulé par Aubert & Duval [Combeau et al., 2009] pour des géométries 2D et 3D. Ces simulations montrent que la prise en compte des configurations réelles de l'écoulement via les parois a un impact sur la prédiction de l'évolution de la solidification. Des phénomènes se déroulant dans l'espace sont bien pris en compte avec des simulations 3D, par exemple le ralentissement des cristaux solides par le frottement le long de parois, la formation 3D des canaux de ségrégation (par exemple le canal enrichi en soluté en bas de la cavité dans cette étude) ou la distribution non uniforme des quantités à travers l'épaisseur, ce qui ne peut pas être capturé par des calculs 2D dans lesquels les variables ne se développent que dans un plan, c'est-à-dire pour un volume infini dans la direction normale à la simulation. En outre, nous remarquons que pour les deux cas 3D étudiés la solution dans le plan de symétrie est proche de celle prédite par le calcul 2D puisque pour ce plan les champs de vitesse sont contraints à rester dans ce plan. Dans nos simulations, nous reprenons l'approximation de [Combeau et al., 2009] L'échange thermique à la surface du lingot est modélisé en appliquant sur la paroi latérale et la surface inférieure une condition aux limites de type Fourier avec deux coefficients de transfert thermique différents: W m -2 K -1 sur la partie supérieure représentant un échange limité via le réfractaire et W m -2 K -1 sur la partie inférieure figurant un échange plus efficace via le corps de la lingotière. Il est supposé que le lingot est déjà rempli par l'alliage liquide à sa température de liquidus de 1503 °C au début de la simulation. La température extérieure évolue au cours du temps : elle est fixée à 626.85 °C (900 K) jusqu'à 2000 s et à 26.85 °C (300 K) à partir de 4000 s ; dans la période intermédiaire elle diminue linéairement en fonction du temps, comme résumé dans le Tab.5. 7.

Cette évolution de la température extérieure sert à modéliser le changement de la température de l'environnement (un moule en fonte) avec lequel le lingot est mis au contact au cours du refroidissement. Une loi instantanée est utilisée pour la germination avec une densité de grains égale à et une surfusion égale à . La fraction d'empilement limite à laquelle la phase solide est bloquée est de 0.4. Seuls des grains équiaxes globulaires sont modélisés dans cette étude. Macro/Micro temps [-] ( )

[-] ( )

[-] Fig. 5. 18 -Schéma 2D du lingot en configuration cartésienne plane. Conditions aux limites en thermique au Tab.5. 7. Conditions aux limites en mécanique (repérées pour les numéros 1 à 4) à l'annexe -A.8. En présence du mouvement du solide, l'écoulement du liquide est également gouverné par des grains solides dans la zone où ceux-ci sont en suspension. Ainsi nous voyons que la vitesse du liquide est toujours dirigée vers le bas le long du front de solidification sous l'effet de la sédimentation des grains, entraînant l'écoulement à remonter au centre du lingot. Dans le cas où la phase solide est fixe (analysée dans les références citées dans notre étude bibliographique Section 2.4 dont un des résultats est montré à la Fig. 2. 8, Case 1), l'écoulement est contrôlé par la convection naturelle thermosolutale. Il descend le long de la paroi refroidie au début de la solidification du fait de l'effet thermique, mais lorsque la solidification est suffisamment avancée, le déplacement du liquide est dominé par l'effet solutal ce qui induit un écoulement ascendant le long du front de solidification et descendant au centre du lingot.

Tab.5. 7 -Conditions aux limites thermiques

[ ) [s] [K] [ ) [s] [K] [s] [K] ; [W m -2 K -1 ]
[m s -1 ] [m s -1 ] [m s -1 ] [m s -1 ] [m s -1 ] [m s -1 ]
Il est constaté que les vitesses de déplacement sont élevées (elles peuvent atteindre un maximum de 250 mm s -1 ), il n'y pas une couche de grains germés constamment attachée à la paroi verticale, comme observé sur les cartes de la fraction de solide au cours du temps (à 10, 20, 50, 100 et 200 s) dans la Fig. 5. 20 dont la courbe blanche représente l'endroit où la fraction solide atteint la valeur limite (0.4) à partir de laquelle les grains sont bloqués. Les champs de vitesse (Fig. 5. 19) et la distribution temporelle de la fraction de solide (Fig. 5. 20) montrent que des grains solides sédimentent près de la paroi refroidie où ils sont initialement formés, et sont ensuite transportés par le liquide vers le centre. Les grains se rencontrent et s'accumulent d'abord au centre, puis s'étendent du centre à la paroi, conduisant à la formation d'une couche solide progressant de bas en haut dans le lingot.

Ces résultats, bien que qualitatifs (analyses en 2D cartésien), confirment l'hypothèse du présent modèle qui suppose que le temps caractéristique de la croissance peut être négligé par rapport à celui du mécanisme de transport. D'autre part, on peut visualiser à la Fig. 5. 21 les distributions des variables à 250 s (de gauche à droite : la fraction de solide, la densité de grains, la composition moyenne, la température et le rayon moyen des grains). On constate que la densité de grains se trouve élevée en bas au centre du lingot et graduellement diminuée vers la paroi. Pendant le processus, le déplacement des grains induit une distribution non homogène de la quantité de grains et influence la cinétique de croissance des grains (e.g. le transport des grains loin du front de solidification entraîne une augmentation de la surfusion dans cette région, la vitesse de croissance est alors accélérée). De plus, les grains solides, pauvres en soluté, occupent la région inférieure du lingot. Le liquide enrichi en soluté est ainsi globalement transféré vers la région supérieure. Cela induit une ségrégation négative dans la zone inférieure du lingot qui est compensée par une ségrégation positive dans sa partie supérieure avec une augmentation progressive de la concentration moyenne du centre à la paroi. Fig. 5. 21 -Calcul 2D plan cartésien. Distribution des variables à 250 s, de gauche à droite : fraction de solide, densité de grains, composition moyenne en carbone, température et rayon moyen des grains. L'iso-ligne blanche affiche la limite de la zone packée à la fraction de solide de 0.4.

Les cartes finales de densité de grains, de macroségrégation et la courbe de ségrégation au centre du lingot sont présentées sur la Fig. 5. 22. A part une grande densité de grains qui se situe en bas du lingot due au dépôt des grains solides, nous voyons qu'une densité de grains également élevée se trouve dans la partie supérieure face à la zone du réfractaire. Cela pourrait être expliqué par le fait que le refroidissement via le réfractaire est limité. Ainsi, le liquide est maintenu à une température élevée et homogène tout en étant continuellement enrichi en soluté dans cette partie supérieure du lingot. Lorsque le liquidus de l'alliage est atteint dans cette région, le front de la zone packée se trouve juste en dessous. L'extraction de la chaleur entraîne alors une germination des grains qui ne pouvant pas sédimenter puisque leur déplacement est empêché par le solide déjà bloqué. D'un autre côté, la sédimentation peut provoquer une raréfaction des grains dans la région en haut du lingot. Aux endroits où des grains sont tous transportés alors que la condition de germination (présentée à la relation (140) dans le Chapitre 3) est encore remplie -la température étant inférieure à la température de liquidus -la germination est à nouveau activée. Cette région est alors composée des cristaux qui germent du fait du départ local de grains. Par conséquent, plus la sédimentation est dominante, plus la fréquence de germination devient importante. L'accumulation continuelle de grains entraîne une densité de grains élevée dans cette région, ce qui n'est pas représentatif de la réalité.

A la fin de solidification

Densité de grains

Composition en carbone Courbe de ségrégation au centre du lingot [grains m -3 ] 〈 〉 [wt.%] Fig. 5. 22 -Calcul 2D plan cartésien. Cartes finales (gauche) de densité de grains et (centre) de macroségrégation.

A droite, courbe de ségrégation au centre du lingot.

En regardant les résultats prédits à la fin de la solidification, nous voyons que les points essentiels des mesures expérimentales sont reproduits par notre simulation numérique pour ce qui concerne la ségrégation, négative en bas et positive en haut du lingot. Cependant, l'intensité de la macroségrégation prédite par notre modèle 2D semble trop élevée en comparaison des observations (cf Fig. 5. 22). Par la suite, nous allons étudier la solidification du lingot dans un espace 3D, toujours en considérant des grains globulaires. Par la suite, des résultats 3D seront analysés et comparés aux résultats obtenus dans le cas 2D précédent. Il est à noter que les configurations de solidification sont différentes pour ces deux simulations: le calcul 2D effectué en cartésien plan, le calcul 3D est réalisé avec une géométrie cylindrique. La Fig. 5. 24 présente les champs de vitesse visualisés sur une coupe verticale, pour lesquels le vecteur indique la direction et la couleur des vecteurs reflète la norme des vitesses. Trois figures (a-c) sont présentées pour la vitesse du liquide et trois autres (d-f) pour la phase solide. Le plan rose représente la surface limite de la zone packée (à la fraction de solide de 0.4) dans laquelle les grains sont bloqués. De plus, la valeur de la composante orthoradiale des vitesses dans un référentiel ( ) est également calculée et affichée sur la coupe verticale et pour quatre sections transversales à 0.4, 0.8, 1.2 et 1.6 m de hauteur depuis le bas du lingot. Il convient d'observer que les vitesses sont moins élevées et semblent plus régulières que celles obtenues par le calcul 2D cartésien (Fig. 5. 19). Des vortex de vitesse sont clairement visibles dans la simulation 2D alors que dans le cas 3D cylindrique, l'écoulement a tendance à descendre le long de la paroi et à remonter au centre du lingot, formant ainsi une cellule globale. Comme attendu et constaté sur les cartes de température de la Fig. 5. 25, le refroidissement est plus rapide avec la configuration 3D cylindrique qu'avec le domaine 2D cartésien. Nous parlerons dans la suite de la différence de l'avancement du refroidissement entre ces deux configurations via l'observation de la fraction de solide. La force motrice générée sous l'effet thermique est moins perturbée dans le cas 3D cylindrique, conduisant à un écoulement plus stable. Plus tard à 200 s, la moindre stabilité des vitesses observée sur les Fig. 5. 24-c) et -f) montre un développement complexe des phénomènes dans la partie haute du lingot en contact avec le réfractaire où la solidification aura lieu plus tard en présence d'un liquide plus enrichi en soluté. D'une part, le maillage utilisé n'est pas idéalement symétrique, et l'écoulement pas parfaitement structuré, donnant lieu à d'inévitables composantes de la vitesse. D'autre part, il se peut également que des instabilités se développant donnent naissance à des composantes orthoradiales non nulles de la vitesse. La Fig. 5. 26 présente les distributions de la fraction de solide au cours de la solidification à 10, 20, 50, 100 et 200s. Au début du refroidissement, jusqu'à 20 s, l'ensemble des grains solides évolue librement dans le lingot cylindrique pour le calcul 3D. Dans le calcul 2D (Fig. 5. 20), à 20 s, une bande packée de solide commence à se former au pied du lingot (comme le révèle la courbure blanche de la fraction de blocage). La formation rapide de la zone packée du calcul 2D résulte d'une forte vitesse de sédimentation des grains, ce que nous avons analysé sur les cartes de vitesse en configuration cartésienne à la A l'issue de cette étude, il est possible de constater que les simulations ont bien reproduit l'essentiel du phénomène de transport des grains quel que soit le calcul, cartésien 2D, ou axisymétrique 3D sur un quart de lingot. Le rôle de la géométrie est également démontré par cette comparaison.

Faute de temps, l'application du modèle des grains dendritiques au lingot n'a pas pu être présentée dans ce manuscrit. Ce modèle à trois phases est cependant nécessaire afin de prendre en compte la morphologie des grains équiaxes puisque comme montré par les mesures (Fig. 2. 8), la partie en bas du lingot est occupée par des grains globulaires alors que des grains dendritiques se trouvent dans la partie supérieure du lingot.

En outre, les écarts entre la solution numérique et les résultats expérimentaux peuvent provenir de plusieurs autres raisons : facteurs numériques (maillage et pas de temps), conditions aux limites estimées sans moule et réfractaire, … ainsi que des phénomènes physiques comme le retrait, l'existence à la fois des structures colonnaires et des grains équiaxes globulaires, dendritiques, les transitions des structures, la présence d'autres éléments d'alliage … autant de facteurs qui ne sont pas pris en compte dans le cadre de notre travail, et qui devraient donc faire l'objet d'études ultérieures.

Résumé

A travers ce chapitre, nous avons appliqué le modèle pour simuler la solidification d'une pièce 3D basée sur les configurations du benchmark de Hebditch-Hunt. Le premier calcul était effectué avec une petite épaisseur pour tester le fonctionnement de la méthode implémentée et vérifier la validité de résultats numériques. Comme attendu, nous avons trouvé une distribution uniforme des variables en traversant l'épaisseur pour une pièce mince en imposant des conditions symétriques aux bords. La simulation du benchmark réel est ensuite réalisée pour la moitié de la cavité. Les résultats obtenus montrent l'influence importante de la prise en compte des parois physiques, ce qui permet de capturer l'effet du frottement entre des parois et des grains mobiles ainsi que ses conséquences sur la macroségrégation et l'avancement de la solidification.

D'autre part, l'application du modèle à un lingot industriel a été réalisée sur un domaine plan rectangulaire 2D dans un repère cartésien et sur un cylindre 3D. L'impact du transport des grains sur le comportement de l'écoulement a été mis en évidence. Globalement les solutions numériques retrouvent la même tendance de ségrégation identifiée par les mesures expérimentales, dans lesquelles les résultats du calcul 3D cylindrique donnent effectivement une meilleure prédiction. Ainsi les configurations géométriques influencent largement le déroulement de la solidification.

Les résultats présentés montrent la faisabilité et la capacité du modèle mis en place pour la simulation de la macroségrégation en 2D et 3D, pour des domaines de petite taille et de grandes dimensions. En particulier, le fonctionnement du modèle actuel en 3D apporte une nouvelle possibilité de simuler et d'étudier la macroségrégation due au transport des grains équiaxes, ce qui joue un rôle important dans les distributions solutale et aussi structurale des produits solidifiés. Au cours des différentes simulations nous avons essayé de donner des premières analyses pour ce type de calculs 3D. Ce travail ouvre une perspective au renforcement de la prédiction des défauts dus au transport des grains solides puisque l'effet de ce phénomène est souvent limité à des simulations 2D alors que dans la situation réelle les pièces solidifiées sont souvent de forme 3D anisotrope.

Cependant ces simulations ont pris un temps de calcul considérable, résumé dans les Tab.5. 8 et Tab.5. 9. Il est possible de voir que la taille de maille est un paramètre qui influe de manière importante sur le temps de calcul. Sur le calcul du benchmark de Hebditch-Hunt 3D, lorsque la taille de maille est diminuée par 2 (de 2 mm à 1 mm), les nombres d'éléments et de noeuds sont environ 7.5 fois plus grands (de 38 619 éléments à 294 935 éléments ; de 8 552 noeuds à 57 626 noeuds), le temps de calcul augmente approximativement de 3 fois (de 9 jours à 29 jours bien qu'un plus grand nombre de coeurs soient utilisés pour le calcul avec la taille de maille de 1 mm).

Plusieurs perspectives pour optimiser le temps de calcul peuvent être évoquées.

 La première possibilité est de raffiner le maillage aux seuls endroits où il existe une variation importante des variables, le maillage étant raisonnablement grossier dans le reste du domaine. Cependant, pour mettre en pratique le remaillage dans le contexte du transport des grains, d'autres études approfondies seront nécessaires. le problème est d'identifier les zones considérées comme critiques, de définir les critères de remaillage correspondants ainsi que d'étudier la compatibilité du remaillage développé à la résolution globale du problème. D'une part, l'application du remaillage influencerait la qualité de la macroségregation prédite, en particulier lorsque cette procédure est fréquemment effectuée. D'autre part, cette réalisation aurait un impact sur la distribution de la diffusion artificielle qui dépend de la taille de maille, le raffinage du maillage aux zones d'empilement peut entraîner une diffusion insuffisante pour éviter des problèmes numériques. L'étude pour cette voie pourra poursuivre le travail de thèse de Gouttebroze [Gouttebroze, 2005] dans lequelle une stratégie de remaillage a été mise en place en supposant une phase solide fixe.

 La deuxième possibilité est l'utilisation d'un pas de temps adaptatif. La résolution de la mécanique des fluides est l'une des étapes prenant du temps de calcul puisque la détermination simultanée des champs de vitesses et de pression serait coûteuse pour obtenir une convergence de ces variables en même temps. Une technique de pas de temps adaptés correspondante au solveur de la mécanique des fluides de CIMLib a été étudiée dans la thèse de Jannoun [Jannoun, 2015]. Par ailleurs, les résolutions de l'ensemble du problème sont obtenues en appliquant un schéma itératif de point fixe, ce qui demande également du temps de calcul.

 Dans l'étude présente, la simulation 3D pour le lingot est faite sur un quart du domaine. Les phénomènes pourraient cependant être approximés sur un sixième ou un huitième du domaine. La réalisation des calculs avec des portions angulaires différentes est possible lorsque l'on travaille avec le logiciel Thercast®. Cela permettrait de considérablement réduire le temps de calcul d'une part et il serait intéressant d'étudier l'effet angulaire afin de choisir une géomètre optimale d'autre part.

Tab.5. 8 -Temps de calcul sur la cavité Sn -5 wt.% Pb Tab.5. 9 -Temps de calcul sur le lingot d'acier 3t3

Chapter 6

Conclusions and Perspectives

Conclusions

During solidification, the solid phase is formed which can be stationary or mobile. The stationary state is often referred to columnar structures which are attached to cooled walls and fixed there. For the later state, the solid movement can be due to the liquid convection and sedimentation. In solidification modeling, while significant developments related to the fixed solid case have been made, less investigation has been considered for the mobile solid case since accounting for the solid motion requires more complex mathematical models as well as more powerful numerical approaches to solve a strongly nonlinear problem. This issue has been investigated for mathematical and computational aspects by several groups as presented in the literature review (Section 2.4), based on the finite volume method. Moreover, as involving high computing costs, applications to predict macrosegregation to industrial large sized castings, for which the phases motion plays an important role, have been limited and often two-dimensional simulations are carried out. In reality, three-dimensional devices associated with complex components are often cast. Finite element methods, in which modeling can be performed on non-structured meshes, provide a useful tool to treat irregular and complex geometries.

In this situation, the present thesis was proposed which aims at modeling solidification accounting for the transport of equiaxed grains in the finite element context. Following the model in 2D axisymmetric structured finite volume developed by Založnik and Combeau [Založnik and Combeau, 2010a], we adapted and implemented this FV model into the FE framework. At the microscopic scale, the model is able to deal with the nucleation and growth of grains including either globular or dendritic structures with a limited solute diffusion in both liquid and solid phases. The local microscopic processes are integrated into macroscopic transport equations by a volume averaging technique. The entire issue is solved using a splitting scheme that provides an effective approach for numerical implementation to deal with a set of non-linear equations describing highly dependent multiscale phenomena.

Through our work, three main tasks accomplished can be summarized:

1. The two-phase model of globular grains was solved by splitting two stagesgrowth and transport stagesand was adapted and implemented in the FE code as three principal steps:

 Implementation of the growth model

The growth model consists in solving a set of equations established to describe the growth kinetics which is controlled by solute diffusion limited in both solid and liquid phases at the microscopic scale. In this stage, two phases are considered: the solid phase which is constituted of equiaxed grains being assumed to grow under a fully globular form, and the surrounding liquid phase.

The numerical implementation was effectuated to locally solve the system at each node of the discretized domain. One-dimensional numerical studies for pure heat conduction illustrated that the solutions obtained by the present model can retrieve those of the standard models (lever rule and Gulliver-Scheil approximations). In addition, thanks to accounting for the growth kinetics, it can predict the recalescence phenomenon which cannot be captured by standard solidification paths.

 Resolution for the transport phenomena, including following sub-works:

 Integration of the thermo-solutal liquid convection with a fixed solid Considering the fixed solid phase, the implemented growth model was then combined with the resolution of the liquid momentum equation. Simulations to predict macrosegregation due to the natural convection were performed, including benchmarks [START_REF] Bellet | Call for contributions to a numerical benchmark problem for 2D columnar solidification of binary alloys[END_REF]Combeau, 2010b]. The results obtained from the present model when simulating with high solute diffusion gave good agreement to those in the literature which use the lever rule. Furthermore, the capability of accounting for a limited solute diffusion in the model was also tested and the finite element solutions showed the same tendencies predicted by using a finite volume method, e.g. a delayed solidification or a favorable formation of meso-segregated channels.

 Development of the FE resolution for transport equations, introducing the solid motion into the model with the presence of the liquid movement while neglecting the growth process Relating to the macroscopic scale, the average transport velocities are intrinsically non-divergence-free fields and present discontinuities when the solid phase builds up as a packed bed. In the actual linear finite element context, evaluating the convection flux on elements of a discretized domain can provoke important irregularities that induce numerical instabilities. To deal with these issues, we introduced an adaptive artificial diffusion into the purely convective equations to be solved, which allows detecting zones subjected to sharp gradients and helps treating reasonably such discontinuities. A special care was made to ensure a consistency between the resolution of these transport equations and other relevant equations when introducing this artificial diffusion. A one-dimensional simulation of the sedimentation phenomenon demonstrated that without this treatment the solution cannot be obtained.

Accompanying this work, a computational module of the solid velocity was implemented. The solid phase, gathering a number of grains having the same characteristics (size, solute content), is considered as a continuous medium. Calculating its velocity consists in solving the solid averaged momentum equation and implementing necessary adaptations, e.g. developing a technique to avoid the penetration into the packed bed when solid grains approach the packed layer.

 Accomplishment of the complete growth-transport two-phase model

Combining the two implemented stages, one for the microscopic growth processes and another for the macroscopic transport phenomena, the complete two-phase model was achieved and then applied to simulate the reference case which was performed by using a finite volume method [Založnik and Combeau, 2010b]. A global agreement was obtained between solutions of these two formulations. Besides analyzing differences, a study of the effect of new numerical parameters developed in the finite element model was also carried out. The study showed that adding a suitable diffusion is important; it should be large enough to deal with numerical problems but well controlled to avoid physically unrealistic solutions.

2. In order to account for the dendritic morphology of solid grains, a three-phase model, being extended from the two-phase model, was then implemented based on the works of Combeau and co-workers [Combeau et al., 2009] and [Bedel, 2014] which use a finite volume method. Numerical implementation was also considered for two principal stages. At the microscopic scale, the growth model is modified in order to account for the evolution of three distinct phases, including the dendritic solid, the interdendritic liquid and the extradendritic liquid. A grain envelope is constituted of the dendritic solid and the interdendritic liquid and its morphology is characterized by the solid fraction contained in it. At the macroscopic scale, two liquid phases are considered to have a same hydrodynamic behavior in order that the like two-phase model is applicable. Different simulations were carried out for zero-, one-, and twodimensional configurations. The solutions obtained from the present work were in reasonable agreement with the finite volume results.

3. Finally, three-dimensional simulations were performed on the Hebditch-Hunt benchmark. Differences relative to the two-dimensional results as well as non-uniform distributions of variables in the third spatial direction were visualized and analyzed. Such numerical evaluations would provide more practical information for experimental characterizations, especially for that realize on anisotropic-form objects or on complex-shape systems with nonsymmetrical heat transfers.

Industrial applications were then carried out to predict the macrosegreation for a 3.3t steel ingot. The typical macrosegregation profile of steel ingots was predicted by the present model: a negative segregation in the lower zone and a positive segregation in the upper zone. The formation of these segregations is predominantly controlled by grains sedimentation and fluid convection. In addition, two studies, one for a twodimensional symmetrical ingot and another for a three-dimensional axisymmetrical ingot, showed a considerable effect of geometries on the solidification sequence.

In our work, the main challenges consisted in the implementation of the FE resolution for purely advective equations and in numerically treating discontinuities due to the solid packing. It required a careful investigation under consideration of less investigation involved using the FE method. Up to our knowledge, the present development is original in the context of the FE method.

A good coherency achieved when comparing the results obtained from the present work to the references verifies our implementation and illustrates the capacity of the model. In a near future, these models will be integrated in the commercial version of Thercast® software, thus providing industry with extended capacities of numerical simulations.

Perspectives

A solidification model accounting for grain transport has been implemented. Although various issues have been worked out in the present thesis, there are still points that would be subjected to subsequent improvements. Furthermore, diverse themes would be interesting for future investigations. Several propositions can be made:

 Computational procedures would require improvement since applications to large industrial product is time-consuming. Using time and mesh adaptation techniques could be tested. Anisotropic remeshing application can be based on the work by Gouttebroze [Gouttebroze, 2015] for solidification modeling with a fixed solid phase.

A technique of adaptive time-steps for the fluid mechanics solver in CIMLib has been developed in the Jannoun's thesis [Jannoun, 2015].

 Homogenous nucleation at the liquidus temperature is considered in the present work. However, in practice, nucleation of equiaxed grains commonly takes place heterogeneously in the bulk undercooled melt. Advanced nucleation models would need being implemented in order to have a better description of physical phenomena. Investigation for this topic can be found in [Avrami, 1940] [Oldfield, 1966] [Rappaz and Thévoz, 1987] [Wu et al., 2003] [Dagner et al., 2008].

 Studying an elasto-viscoplastic behavior of the solid phase, e.g. for the packed bed which is formed by the sedimentation of equiaxed grains, and its resulting effect like induced shrinkage as well as residual stresses and corresponding deformations would provide further insight into phenomena during the solidification process. Relevant studies about this matter can be consulted in [Bellet, 2007] [Ludwig et al., 2015b].

 The presence of multi-chemical species and the formation of multiple phases are always part of alloy solidification. To take into account these phenomena, the micromodel should be extended in order to describe the microscopic exchange of different chemical species and then combined with the macroscopic resolution in which the transport of relevant solute elements requires more equations to be integrated (at least two transport equations for each alloying element).

 Application of the development to real casting processes (through Thercast® implementation), including application during mold filling and application in the context of continuous casting.

 In addition, according to the work by Leriche and co-workers [Leriche et al., 2015], the present model could be applicable to model a microstructural transition (columnarto-equiaxed transition), that would provide a better prediction of solidification structures, allowing for a further control to produce castings with expected qualifications.

  Fig.1. 2-b) and Fig.1. 2-c) display the profiles of macrosegregation in carbon and in sulfur, respectively, along the center of the ingot at the end of solidification. According to the authors, although a similar tendency was generally achieved between numerical solutions and experimental measures, the negative segregation in the lower part of the ingot could not be predicted because of neglecting the sedimentation of solid grains. In addition, the positive segregation in the upper part was underestimated for the case of carbon since the absorption of carbon from insulation materials at the top of the ingot was not included.

  Fig.1. 2 -Macrosegregation modeling without solid motion, results issued from [Gu and Beckermann, 1999]. (a) Schematic of the 2D simulated ingot. Comparison between numerical and experimental segregation profiles (b) in carbon and (c) in sulfur along the center line of the ingot.

  3 -Segregation simulation applied to a 65t steel ingot[Combeau et al., 2012a] with a) final relative segregation profiles along the center line of the ingot: (green curve) experiment results, (blue curve) without solid motion, (red curve) with solid motion; and b) final relative segregation maps (left) without solid motion and (right) with solid motion.

Fig. 1

 1 Fig.1. 4 -Phase diagram for a binary alloy (A-B).

Fig. 1

 1 Fig.1. 5 -1D illustrative schematic of solute profiles in a mushy zone at temperature T (referred to Fig.1. 4)

Fig. 2

 2 Fig.2. 2 -Results obtained from the two-phase modelissued from[Beckermann and Ni, 1996]. Liquid velocity vectors and solid volume fraction contours at 30 s for three simulations:Case 1 without solid motion.Case 2 with solid motion and nucleation rate at grains m -3 . Case 3 with solid motion and nucleation rate at grains m -3 .

  Fig.2. 3a) Schematic of an averaging volume containing several dendritic grains. Three phases including solid dendrite, interdendritic liquid and extradendritic liquid. b) Partitioning of the liquid flow : v ftotal liquid velocity ; v dinterdendritic liquid velocity ; v lextradendritic liquid velocity. Figures issued from[Wang and Beckermann, 1996a].

  Fig.2. 5 -Simulations issued from [Wu et al., 2003]. a) Schematic of the studied cavity. b) Contours of the solid fraction and liquid velocity vectors. c) Solute distribution. for b) and c) : Case 1 without solid transport on LHS & Case 2 with solid transport on RHS

3

 3 Models of Combeau and co-workers (SOLID software (*))

Fig. 2

 2 Fig.2. 7 -Grain nucleation with an instantaneous non-uniform law. Grain density depending on undercooling and arranged in classes with different grain sizes,Illustrative figures issued from[Bedel, 2014].

Fig

  Fig.2. 8 -Ingot geometry with the map of grain structures identified from experiments and simulation results at 1800 s, issued from[Combeau et al., 2009], including three case studies: Case 1 fixed solid phase. Case 2 globular grain transport.  Case 3 dendritic grain transport. For each case, LHS : solid fraction map and solid velocity vectors.RHS : segregation ratio map and liquid velocity vectors/streamlines.

Fig. 2 .

 2 Fig.2. 10 -Schematic of a continuous casting slab & the liquid fraction distribution along the slab. Results illustrated on a thick slice of the slab (marked by the red rectangle) including: a) Relative velocity vectors in the mushy zone displayed with (red) solidus isotherms. b) Divergence distribution of the intrinsic solid velocity (negative values in blue color for contraction zones ; positive values in red color for expansion zones). c) Segregation map (negative values in blue color ; positive values in red color).Figures issued from[Fachinotti et al., 2006].

  Fig.2. 11a) Configurations of the simulated cavity. b) Cooling curves (black dashed curves) measured by experiment and (colored solid curves) predicted by CAFE simulation, recorded at four locations L21, L24, L27, L30 being situated at 5, 35, 75, 95 mm from the left side at the mid-height of the cavity, issued from [Carozzani, 2012].

Fig. 3 .

 3 Fig.3. 2a) Illustration of a support domain, S, with black dot SPH makers.b) Schematics of BCE makers (black circles) distributed on a solid object (gray area) and SHP neighbors (white circles). Figures issued from[Pazouki and Negrut, 2015].

Fig. 3 . 3 -

 33 Fig.3. 3 -Schematic of the resolution algorithm using the splitting method, introduced the solutions obtained from the transport stage.

  98) Chapter 3 -Finite Element Implementation and Adaptation 3.1.1.2-c) Resolution Method

  105) Chapter 3 -Finite Element Implementation and Adaptation

Fig. 3 .

 3 Fig.3. 5 -Transition function, 𝛼 ( [𝛼 (𝑔 𝑔 )]), vs solid fraction where the value of the packing solid fraction 𝑔 is equal to 0.3.

Fig. 3 .

 3 Fig.3. 10 -Transport-growth resolution algorithm for the two-phase model.

Fig. 3 .

 3 Fig.3. 12a) Schematic of the two-phase model (spherical solid shells and the surrounding liquid). b) Schematic of the three-phase model (the solid dendrite, the interdendritic liquid and the extradendritic liquid).

  between the solid phase and the interdendritic liquid phase.

Fig. 3 .Fig. 3 .

 33 Fig.3. 13 -Schematic of the resolution algorithm using the splitting method, for the three-phase model.

Chapter 3 -

 3 Fig.3. 16 -Resolution algorithm at the growth stage for the three-phase model.

  Fig.4. 1 -Schematic of the simulation case.

Fig. 4

 4 Fig.4. 2 -Time evolutions of (left) the temperature and (right) the solid fraction obtained from (orange curves) the lever rule and from (green and red curves) the present two-phase model considering a high solute diffusion in the solid phase.

Fig.4. 3 -

 3 Fig.4. 3 -Time evolutions of (left) the temperature and (right) the solid fraction obtained from (orange curves) the Gulliver-Scheil approximation and from (green and red curves) the present two-phase model considering no solute diffusion in the solid phase.

  Fig.4. 4 -Schematics of system dimensions and boundary conditions of 2D simulation cases (accompanied with Tab.4. 1). a) Case 1 performed on a half of a 0.1 m x 0.06 m sample using a Pb-18wt.%Sn. b) Case 2 performed on an entire 0.1 m x 0.06 m sample with a Sn-5wt.%Pb.

  Fig.4. 7 -Case 1 -Concentration profiles at the three height levels of the cavity, including (top) 10 mm, (center) 30 mm and (bottom) 50 mm, obtained (left) from the present model and (right) from various numerical solutions issued from [Combeau et al., 2012b].

Fig. 4

 4 Fig.4. 8 -Final average composition maps obtained by the present two-phase model, considering a limited solute diffusion.

Fig

  Fig.4. 9 -Geometry and initial states of variables set up for the sedimentation simulation.

Fig

  Fig.4. 15 -Vertical profiles of (left) the solid fraction and (right) the temperature at 0 s, 2 s, 5 s, 10 s and 15 s.

  Fig.4. 17 -Recalling schematic of system dimensions (identical to Fig.4. 4-b) Same thermal boundary conditions were described in Section 4.2.1. Mechanical boundary conditions are presented in Tab.4. 3.

  Fig.4. 18 -Case with solid transport -Distributions of variables including a) Solid fraction b) Gain density c) Temperature d) Average composition e) Liquid velocity f) Solid velocity The three black isolines represent the solid fraction at 0.1 (upper isoline), 0.2 (intermediate isoline) and 0.3 (lower isoline). Arrow showing the direction of velocities, color depicting their magnitude.

  Fig.4. 19 -Case with solid transport -Distributions of variables including a) Solid fraction b) Grain density c) Temperature d) Average composition e) Liquid velocity f) Solid velocity The three black isolines represent the solid fraction at 0.1 (upper isoline), 0.2 (intermediate isoline) and 0.3 (lower isoline). Arrow showing the direction of velocities, color depicting their magnitude.

  Fig.4. 22 -Case with fixed solid phase -Distributions of variables including a) Solid fraction b) Grain density c) Temperature d) Average composition e) Liquid velocity f) Solid velocity The three black isolines represent the solid fraction at 0.1, 0.2 and 0.3. Arrow showing the direction of the velocities, color depicting their magnitude.

Fig. 4 .FinalFinalFinal

 4 Fig.4. 31 -Time evolution of the temperature (left) the present FEM results and (right) the FVM results issued from [Bedel, 2014] (green curves) two-phase model, (orange curves) three-phase model with hemispherical tips, (purple curves) three-phase model with paraboloidal tips. The maximal difference between the FE and FV results is smaller than 1°C.

  Fig.4. 36 -Time evolutions of (left) the temperature and (right) the solid fraction obtained from (green curves) the present three-phase model with paraboloidal tips and (orange curves) the lever rule and Gulliver-Scheil models

Fig. 4 .

 4 Fig.4. 37 -Schematic of system dimensions Thermal and mechanical boundary conditions are presented in Tab.4. 4.

  Fig.4. 38 -FEM results obtained using (left) three-phase model with hemispherical tips, (center) three-phase model with paraboloidal tips, (right) two-phase model from top to bottom: temperature, solid fraction, internal solid fraction and grain density

  Fig.4. 41 -FVM (SOLID software) results obtained using the three-phase model with hemispherical tips, a) Temperature b) Solid fraction c) Internal solid fraction d) Grain density

5. 1 3 -

 13 Fig.5. 1 -Schéma du cas 3D pour une pièce mince de 1 mm d'épaisseur Les conditions aux limites thermiques et mécaniques sont repérées pour les numéros 1 à 6 accessibles aux Tab.5. 1 et Tab.5. 2.

  des simulations Dans cette partie, le modèle est appliqué pour simuler la solidification d'un lingot d'acier de 3.3 tonnes coulé par Aubert & Duval. Nous avons présenté les travaux dédiés à ce lingot établis par Combeau et al.[Combeau et al., 2009] dans notre revue de la littérature (Section 2.4.3). Les points principaux concernant les analyses sur ce lingot sont rappelés et résumés ici. La configuration expérimentale du lingot est représentée sur la Fig.5. 17-a). L'acier liquide remplit le lingot à partir du pied, en source. La partie inférieure du lingot est en contact avec le moule et celle supérieure est en contact avec un réfractaire. Les parois externes du moule échangent thermiquement avec l'air extérieur par convection naturelle et rayonnement. La surface supérieure est couverte par une couche de poudre exothermique. Le temps de solidification total est de l'ordre de 90 minutes. La Fig.5. 17-b) et la courbe rouge dans la Fig.5. 17-c) présentent respectivement la carte de macroségrégation interpolée à partir des mesures expérimentales (114 points) et le profil de ségrégation au centre du lingot dont une ségrégation négative se situe en bas et une ségrégation positive élevée se trouve dans la zone supérieure en contact du réfractaire. Les courbes bleue et verte sur la Fig.5. 17-c) sont les résultats numériques du travail de Combeau et al. en supposant la croissance de grains globulitiques d'une part et dendritiques d'autre part.

  Fig.5. 17 -Configuration du lingot et résultats issus de l'article [Combeau et al., 2009] : a) Configuration du lingot d'acier de 3.3 tonnes. b) Carte de macroségrégation expérimentale interpolée dans une demisection longitudinale du lingot. c) Courbes de ségrégation au centre du lingot.

  5.2.2 Etude sur le lingot en configuration plane cartésienne 2DLa Fig.5. 19 présente les vitesses du liquide et du solide à 10 s, 100 s et 200 s avec les isolignes noires de la fraction de solide à 0.1, 0.2, 0.3 et 0.4. Les résultats à 10 s montrent l'écoulement du liquide descendant proche de la paroi refroidie du fait de l'effet thermique, renforcé par la sédimentation des grains solides. Le liquide remonte au centre du lingot et emporte une partie du solide vers le haut. On peut constater la présence de petits vortex générés par les interactions entre les deux sens de l'écoulement. Plus tard, à 100 s et 200 s, les vortex se développent de manière plus nette, les mouvements du solide et du liquide adoptent une forme de cellules. Ceci révèle des interactions complexes des phénomènes. En outre, les évolutions similaires des vitesses du liquide et du solide reflètent une dépendance bien couplée du mouvement entre les deux phases -un entraînement important du mouvement des grains solides par le liquide et vice-versa.

Fig. 5 .

 5 Fig.5. 19 -Calcul 2D plan cartésien. Vitesses du liquide et du solide aux instants 10 s, 100 s et 200 s. Les isolignes noires de la fraction de solide correspondent à 0.1, 0.2, 0.3, 0.4 et 1.

Fraction

  Fig.5. 20 -Calcul 2D plan cartésien. Distribution de la fraction de solide aux différents instants 10, 20, 50, 100 et 200 s.L'iso-ligne blanche affiche la limite de la zone packée à la fraction de solide de 0.4.

5. 2 . 3

 23 Fig.5. 23 -Schéma 3D d'un quart de lingot à étudier.Les conditions aux limites en thermique sont présentées dans le Tab.5. 7 (identiques à celles du cas 2D) et celles en mécanique (repérées pour les numéros 1 à 5) à l'annexe A.9.

  Fig.5. 19, par rapport à celles en 3D cylindrique à la Fig.5. 24. Par contre jusqu'à 200 s, le front d'empilement du solide dans le cas 3D (plan rose, Fig.5. 26) dépasse la moitié de la hauteur lorsqu'au même moment ce front n'arrive pas encore à la moitié du lingot dans le calcul 2D cartésien (courbe blanche, Fig.5. 20). La solidification plus avancée dans un lingot cylindrique 3D est justifiée par le fait que l'extraction de chaleur surfacique rapportée au volume est plus élevée en 3D qu'en 2D plan.

  Fig.5. 24 -Calcul 3D, résultats à 10 s, 100 s et 200 s (les vecteurs indiquent la direction des vitesses, la couleur des vecteurs indique la norme des vitesses, les cartes présentent la valeur des composantes de la vitesse) (a-c) résultats pour la phase liquide et (d-f) résultats pour la phase solide Le plan rose affiche la limite de la zone packée (à la fraction solide de 0.4).

FigFraction

  Fig.5. 25 -Cartes de température à 10 s, 100 s et 200 s (a-c) résultats du calcul 3D cylindrique et (d-f) résultats du calcul 2D cartésien.

  It is developed by Institut Jean Lamour in Nancy and Sciences Computers Consultants. (http://www.scconsultants.com) Chapter 2 -Literature Reviews

	Slurry regime :

  . 3.

			Compared cases in Fig.4. 3
	two top sub-figures		two bottom sub-figures
	"infinite" case (present two-phase model)	"finite" case (present two-phase model)
		& Gulliver-Scheil		& "infinite" case (present two-phase model)
	where
	Cases	[m 2 s -1 ]	[grain m -3 ]	[m 2 s -1 ]	Plotted curves
	"infinite"					green curves
	"finite"					red curves

  : l'acier est considéré comme un alliage binaire Fe-0.36wt.%C dont les propriétés sont données à l'annexe A.8. Ce choix est justifié car, selon Combeau et al., l'élément carbone joue un rôle significativement dominant par rapport aux autres éléments contenus dans la force motrice résultant des gradients thermique et solutaux. Par ailleurs, le lingot réel est de section octogonale. Nous avons pour notre part considéré deux types de modélisation : une en approche 2D plane (l'approche 2D axisymétrique n'étant pas disponible) et une en approche 3D à section transverse circulaire, donc plus proche du problème réel. La surface en haut du lingot est supposée adiabatique. L'assemblage des éléments moulants à corps multiples, dans lequel chaque composant est caractérisé par un matériau différent, ne peut pas être considéré. Dans le travail actuel, nous avons simulé un domaine qui ne représente que le lingot.

Recalescence

Acknowledgments Chapter 3 -Finite Element Implementation and Adaptation Fig.3. 15 -Transport-growth resolution algorithm for the three-phase model. Chapter 4 -Numerical Simulation and Validation Fig.4. 35 -Time evolutions of (left) the temperature and (right) the solid fraction obtained from (green curves) the present three-phase model with hemispherical tips and (orange curves) the lever rule and Gulliver-Scheil models.

where ( ) is the Sobolev space which is given as

with ( ) the Lebesgue space :

A zero residual remains unchanged when being multiplied with a test function and integrated over the domain of interest. By multiplying by the test function (using the SUPG method), the weak formulation is defined as to find such that Eq. ( 58) is satisfied.

with (59) where is a positive stabilization parameter. An optimal formulation of has not been generalized for all situations. By default in the existing FE codes, the parameter is proposed as (60), which is generally considered as a guaranteed choice for most cases with unstructured meshes. Further discussion and test comparisons for different SUPG stabilization parameters can be found in [Akin andTezduyar, 2004] [John andKnobloch, 2007].

where is the characteristic mesh size of an element in the direction of the velocity , as defined in (46);

is the velocity at the center of the element.

If the set space has infinite dimension and is composed of independent functions, Eq. ( 58) and Eq.( 51) are equivalent. If it's not the case, the weak solution of Eq.( 58) is only an approximate solution to the problem.

Comparison FEM and FVM solutions for the solid transport case

On the other hand, regarding the FEM solution for the solid transport case in Fig. 4. 21 compared to those given from the FVM in Fig. 4. 26, a similar tendency is stated: a negative segregation in the left zone resulting from the accumulation of solute-depleted grains and a positive segregate channel at the bottom of the cavity due to the circulation of the liquid phase in the mushy zone and the remelting. Concerning the grain density maps, a similar number of grains are distributed in the bottom-right and upper-left corners for both cases.

The differences between these two results can be explained by the influence of numerical factors. Besides those mentioned in the previous section for the fixed solid case, introducing an artificial diffusion in FEM resolutions induces a less intense segregation and a more diffused distribution of grain density, when comparing to the FVM results.

at the end of the solidification

Average composition 〈 〉 [wt.%]

Grain density [grains m -3 ] Fig. 4. 26 -Case with solid transport -FVM results obtained from the SOLID software. Distributions of (left) the average composition and (right) the grain density.

In order to choose an appropriate artificial diffusion, meaning that it must be sufficient to stabilize numerical issues and reasonably reduce discontinuities but not cause so much diffusion, a sensitivity analysis of the model to the artificial diffusion coefficient was presented for the 1D calculations in Section 4.3.3. This aspect will be further studied with the actual 2D simulation case in the next section.

Annexes Models for the Solute Diffusion Lengths

For the two-phase globular model, the solute diffusion length in the liquid phase is issued from the work of Tveito and co-workers [Tveito et al., 2011]. It is expressed as the following formulation.

(195)

where is the liquid composition at the liquid-solid interface in the liquid side; 〈 〉 is the average composition in the liquid phase; is the liquid composition; is the unit vector normal to the liquid-solid interface;

is the solute diffusion in the liquid phase; is the growth velocity of the solid phase; is the grain radius; is the final grain radius; is the stagnant film thickness; is the kinematic viscosity; 〈 〉 is the intrinsic liquid velocity; 〈 〉 is the intrinsic solid velocity.

Related to the area concentration:

In the three-phase dendritic model, the exchange surface of solute is the interface between the grain and the extradendritic liquid phase, instead of the one between the solid and the liquid in the two-phase model. The grain radius used in Error! Reference source not found. is replaced by the approximated envelop radius which is equal to the radius of a spherical grain having the same volume as follows.

Data and Simulation Parameters

A.1 Mono-dimensional Solidification (Section 4.1)

Thermophysical data of an Al -4 wt.% Cu alloy

Simulation Parameters

Mesh size Modélisation multi-échelle d'éléments finis de la macroségrégation et du transport des grains RESUME : Ce travail de thèse a pour but de modéliser la macroségrégation des produits obtenus par solidification en prenant en compte le transport des grains équiaxes. Le modèle de solidification à deux phases (solide et liquide) est traité par une méthode d'éléments finis, consistant à résoudre les équations de conservation moyennées de l'énergie, de la quantité de mouvement et de la masse, dans lesquelles les évolutions multiéchelles de la masse des phases et des solutés sont modélisées en utilisant une approche de splitting. D'après cette technique, la variation des quantités est considérée comme résultant de la contribution de deux étapes : l'étape de croissance et l'étape de transport. L'implémentation numérique du modèle a été réalisée avec trois opérations principales : tout d'abord implémenter le modèle de croissance des grains, ensuite intégrer des phénomènes de transport résultant de la convection thermo-solutale du liquide et du mouvement du solide, enfin mettre en oeuvre le modèle complet en combinant les étapes de croissance et de transport. Lors de ces opérations, une investigation attentive a été consacrée à l'établissement de la résolution par éléments finis pour les équations de transport avec champs discontinus de vitesse à divergence non nulle, afin de surmonter des problèmes numériques en respectant la qualité des solutions physiques. Parallèlement à ces travaux, différents tests de simulation 2D ont été effectués à chaque étape d'implémentation. De bons accords ont été globalement obtenus entre les solutions données par le modèle présent et celles de référence dans la littérature. Finalement, des applications industrielles et des simulations 3D ont été menées, pour lesquelles les résultats numériques reproduisent les configurations caractéristiques des mesures expérimentales : un profil typique de macroségrégation des lingots composé d'une ségrégation négative occupant de la zone inférieure et d'une ségrégation positive en zone supérieure. Ce profil est la signature caractéristique de la sédimentation des cristaux et de la convection thermo-solutale. En outre, un modèle à trois phases étendu à partir du modèle à deux phases précédenten distinguant la phase liquide interdendritique -s'est avéré capable de décrire la morphologie des grains dendritiques.

Mots clés : solidification, modélisation, macroségrégation, aciers, multi-échelle, éléments finis, transport

Multiscale finite element modeling of macrosegregation and grain transport ABSTRACT :

The present work aims at modeling macrosegregation of castings, accounting for the transport of equiaxed grains. A two-phase (solid and liquid) finite element solidification model is presented, consisting in solving a system of volume-averaged conservation equations of energy, momentum, solute, in which the multiscale evolutions of phase and solute mass are modeled by using a splitting method. According to this approach, the variation of quantities is considered as due to the contribution of two stages: the growth stage and the transport stage. The numerical implementation was realized with three principal steps: first implementing growth processes, then integrating transport phenomena including the thermo-solutal liquid convection and the solid movement, lastly combing the growth and transport stages to achieve a complete growth-transport model. Of these steps, solving the transport equations with discontinuous and non-divergence-free velocity fields by using finite element method required an attentive investigation in order to overcome numerical issues while respecting for physical solutions. Parallel to these works, various two-dimensional simulation tests were carried out in each implementation step. Agreements were globally found between results obtained from the present model and those of reference from the literature. Finally, industrial applications and three-dimensional simulations were performed, which show that computational solutions can predict essential features of experimental measurements. In particular, a typical macrosegregation profile of steel ingots, containing a negative segregation in the lower zone and a positive segregation in the upper zone, which is predominantly characterized by crystals sedimentation and fluid circulation was retrieved. Moreover, a three-phase model considering two different liquid phases, extended from the above-mentioned two-phase model, which enables to describe the morphology of dendritic solid crystals was implemented.