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

Victor D. Fachinotti, Steven Le Corre, Nicolas Triolet, Manuel Bobadilla, Michel Bellet

► To cite this version:

Victor D. Fachinotti, Steven Le Corre, Nicolas Triolet, Manuel Bobadilla, Michel Bellet. Two-phase thermo-mechanical and macrosegregation modelling of binary alloys solidification with emphasis on the secondary cooling stage of steel slab continuous casting processes. International Journal for Numerical Methods in Engineering, 2006, 67 (10), pp.1341-1384 10.1002/nme.1664. hal-01581686

HAL Id: hal-01581686

https://hal.science/hal-01581686

Submitted on 5 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Two-Phase Thermo-Mechanical and Macrosegregation Modelling
of Binary Alloys Solidification with Emphasis on the Secondary
Cooling Stage of Steel Slab Continuous Casting Processes

Victor D. Fachinotti¹,³, Steven Le Corre¹, Nicolas Triolet², Manuel Bobadilla² and Michel Bellet¹

¹ Ecole des Mines de Paris

Centre de Mise en Forme des Matériaux (CEMEF)

UMR CNRS 7635, BP 207, F-06904 Sophia Antipolis, France

² Arcelor Research, BP 30320, 57283 Maizières-lès-Metz, France

³ Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)

Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)

Parque Tecnológico del Litoral Centro - Ruta Nacional 168, Paraje El Pozo

S3007ABA Santa Fe, Argentina

As an approach towards a better modelling of solidification problems, we introduce a thermo-mechanical and macrosegregation model that considers a solidifying alloy as a binary mixture made of a liquid and a solid phase. Macroscopic conservation laws for mass, momentum and solute are obtained by spatial averaging of the respective microscopic conservation equations. Assuming local thermal equilibrium, a single equation for the conservation of the mixture energy is then written. A single equation can be obtained for the solute as well by invoking a proper microsegregation rule. The numerical implementation in a 2D finite element code is then detailed. Lastly, some examples of simulations of academic tests as well as industrial applications for continuous casting of steel slabs are discussed. They particularly enlighten the ability of the formulation to describe the formation of central macrosegregation during the secondary cooling of slab continuous casting processes.

KEY WORDS: thermo-mechanical, macrosegregation, two-phase, continuous casting

1. Introduction

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

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

A second category of mushy zone models is composed of two-phase models, but assuming a fixed rigid solid phase. Such models have been developed either by using the mixture theory or by applying the spatial averaging method to the conservation equations (mass, momentum, energy, solute) on a representative elementary volume (REV) of the mushy material. They have been proposed to study solidification in columnar or equiaxed regime, accounting in this latter case for the transport and settling of solid grains. The interaction between solid and liquid phases is expressed by the Darcy law. Such models have been initially developed and implemented in finite volume numerical models by Bennon and Incropera [5] and Beckermann et al (Ni and Beckermann [6], Wang and Beckermann [7]). Finite element implementations have been proposed by Ahmad [8] in two dimensions and by Desbiolles et al [9] and Bellet et al [10] in three dimensions.

A third family of models consists of effective two-phase models, in which the deformation of the solid phase is accounted for. These models derive from the soils mechanics theory and apply to the mushy zone when it is "coherent" and saturated with liquid. The "coherency" of a mushy zone is reached when the solid phase is continuous enough to sustain stresses significantly higher than those in the liquid phase. The corresponding solid fraction, $g_{s_{cohe}}$, is typically in the range 0.6-0.7 (the lower values corresponding to a dendritic structure and
the higher ones to a globular structure). At solid fractions higher than $g_{s_{cohe}}$, the solid phase constitutes a continuous solid skeleton which is deformable and in strong interaction with the interstitial fluid. The deformation of this skeleton determines the flow of the liquid phase, which is governed by the Darcy law. The pressure of the interstitial liquid determines in turn the deformation of the solid skeleton. Schematically, the mushy zone then is supposed to behave like a sponge saturated with liquid. This kind of approach has been found quite relevant for slow deformation of semi-solid alloys, especially when liquid segregation can be observed during the deformation. Like the second family, the models proposed in the literature are based on a spatial averaging approach. However, most of them are limited to isothermal conditions and thus cannot be applied to the study of solidification process, as they do not account for mass transfer from the liquid phase to the solid phase that is caused by solidification or fusion. We can quote in this category the following references: Lalli [11], Nguyen et al [12], and Martin et al [13]. The implementation in two-dimensional finite elements numerical models has been developed by Toyoshima and Takahashi [14], Bay et al [15]. In this third family, Mo and coworkers (M’Hamdi et al [16], Nicolli et al [17]) have proposed, in view of studying the continuous casting of aluminium alloys, a model in which mass transfer is accounted for. In their two-dimensional finite element implementation, the conservation of momentum of both the liquid and the solid phase is solved, including Darcy interaction terms, together with the resolution of the conservation of mass, energy and solute for the whole mixture. However, the coupling with solute transport (macrosegregation) is not accounted for: a unique relation between solid fraction and temperature must be prescribed a priori, which is only acceptable for small segregations. In addition, this approach is limited to very low differential velocities between the two phases, the effects of differential advection being neglected in the equation for energy conservation.

In the present paper, we focus on the thermo-mechanical modelling of alloys by the end of solidification, that is to say for cases where the solid fraction $g_s$ exceeds $g_{s_{cohe}}$. The mushy material is then considered as an effective two-phase medium. On one hand, the solid material is considered as an incompressible viscoplastic material, obeying a constitutive equation of power-law type. Invoking homogenization results [18], its macroscopic flow rule is viscoplastic, including compressibility, so that the solid continuum can be seen as a deformable compressible porous medium. On the other hand, the liquid phase is intrinsically Newtonian. At the macroscopic scale, its momentum interaction with the solid skeleton is described by the Darcy law of flow through a porous medium. Then, the present contribution enters the third family of models previously mentioned. In addition, the present formulation includes mass transfer between liquid and solid, and the couplings with heat transfer and solute transport is taken into account. The two-dimensional numerical implementation of the corresponding mass, momentum, energy and solute conservation equations has been carried out in the finite element software R2SOL, developed at CEMEF.

The paper is structured as follows. In Section 2, the coupled conservation equations are presented. In Section 3, their space and time discretisation and the finite element resolution are detailed. In Section 4, the coupling algorithms are discussed. Finally, in Section 5 a validation test of the two-phase mechanical solver is presented, followed by an application to the continuous casting process.
## 2. General macroscopic balance equations

At the microscopic scale, inside each phase, the thermo-mechanical evolution is assumed to be governed by the usual mass, momentum, energy and solute balances. In this work, the balance equations for the mixture, at the (macroscopic) scale of a representative elementary volume (REV), are obtained using the spatial averaging method over a fixed control volume $V_0$. This method is classical and will not be detailed here (see [19, 6, 7, 20] for further details on its basic principles). The notations used in this paper can be summarized as follows. Denoting $\chi_k$ the characteristic function of phase $k$ (= 1 in phase $k$ and 0 elsewhere), for any function $\psi$ we define:

$$ \text{the intrinsic average value:} \qquad \psi_k = \langle \psi \rangle^k = \frac{1}{V_k} \int_{V_0} \psi(\mathbf{x}) \chi_k(\mathbf{x}) \, dV, \qquad (1) $$

$$ \text{the classic average value:} \qquad \langle \psi^k \rangle = \frac{1}{V_0} \int_{V_0} \psi(\mathbf{x}) \chi_k(\mathbf{x}) \, dV, \qquad (2) $$

where $V_k$ denotes the volume occupied by phase $k$ in the REV.

By introducing the volume fraction of phase $k$, defined as

$$ g_k = \frac{V_k}{V_0}, \qquad (3) $$

both average values (1) and (2) become related by

$$ \langle \psi^k \rangle = g_k \psi_k. \qquad (4) $$

Besides, since the solidifying alloy in the mushy state is considered as a saturated two-phase medium, the volume fractions of solid and liquid phases always satisfy the following relationship:

$$ g_s + g_l = 1. \qquad (5) $$

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

$$ \text{Momentum:} \qquad \nabla \cdot (g_k \boldsymbol{\sigma}_k) + \mathbf{M}_k + g_k \rho_k \mathbf{g} = \frac{\partial}{\partial t} (g_k \rho_k \mathbf{v}_k) + \nabla \cdot (g_k \rho_k \mathbf{v}_k \times \mathbf{v}_k), \qquad (6) $$

$$ \text{Mass:} \qquad \frac{\partial}{\partial t} (g_k \rho_k) + \nabla \cdot (g_k \rho_k \mathbf{v}_k) = \Gamma_k, \qquad (7) $$

$$ \text{Energy:} \qquad \frac{\partial}{\partial t} (g_k \rho_k h_k) + \nabla \cdot (g_k \rho_k h_k \mathbf{v}_k) + \nabla \cdot \langle \mathbf{q}^k \rangle = Q_k, \qquad (8) $$

$$ \text{Solute:} \qquad \frac{\partial}{\partial t} (g_k w_k) + \nabla \cdot (g_k w_k \mathbf{v}_k) + \nabla \cdot \langle \mathbf{j}^k \rangle = J_k, \qquad (9) $$

where $\rho$ denotes the density, $\mathbf{v}$ the velocity field, $\boldsymbol{\sigma}$ the Cauchy stress tensor, $\mathbf{g}$ the gravity vector, $h$ the enthalpy per unit mass, $\mathbf{q}$ the heat flow vector, $w$ the solute concentration per unit volume, and $\mathbf{j}$ its flux. The terms $\Gamma$, $\mathbf{M}$, $Q$ and $J$ are associated with the exchanges of mass, momentum, energy and solute, respectively, between the two phases.

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

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

## 2.1. Mass conservation

The local mass balance at the interface between phases ensures that $\Gamma_s + \Gamma_l = 0$ [20]. Then, summing equations (7) for the liquid and solid phases, and assuming that the densities of the two phases remain constant in the solidification interval (but not necessarily equal), we get:

$$ \rho_s \nabla \cdot (g_s \mathbf{v}_s) + \rho_l \nabla \cdot (g_l \mathbf{v}_l) = \frac{\partial g_s}{\partial t} (\rho_l - \rho_s), \quad (10) $$

where the saturation assumption (5) was introduced. Then, dividing by $\rho_l$, equation (10) can be rewritten as

$$ (1 - \Delta\varepsilon^{tr})\nabla \cdot (g_s \mathbf{v}_s) + \nabla \cdot (g_l \mathbf{v}_l) = \frac{\partial g_s}{\partial t} \Delta\varepsilon^{tr}, \quad (11) $$

where

$$ \Delta\varepsilon^{tr} = \frac{\rho_l - \rho_s}{\rho_l} \quad (12) $$

denotes the relative change of volume associated with solidification (often negative for metallic alloys).

It should be noted that the assumption of constant liquid and solid densities affects the solidification interval only. In this interval, $\rho_s$ and $\rho_l$ may be taken equal to the density of the alloy at nominal solidus and liquidus temperature, respectively. Outside this interval, only one phase is present and the mass equation results from equation (7), yielding the usual one-phase mass conservation equation: $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0$.

## 2.2. Momentum conservation

The spatial averaging method used in this work is efficient to obtain in a simple way the macroscopic governing equations of the semi-solid alloy but does not enable to go further in the specifications of the macroscopic model. Reliable constitutive equations would require more sophisticated approaches such as homogenization [21, 18] associated with numerical simulation at the microscopic scale, but this is not within the scope of this work. The full definition of the two-phase model will simply be based on further constitutive assumptions, consistent with previous theoretical works.

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

$$ s = 2\mu_l \dot{\varepsilon}(\mathbf{v}), \quad (13) $$

where $s = \sigma + pI$ is the deviatoric part of the stress tensor $\sigma$, $p = -\text{tr}(\sigma)/3$ being the hydrostatic pressure, $\mu_l$ the viscosity of the liquid, and $\dot{\varepsilon}(\mathbf{v})$ the strain rate tensor, defined as

$$ \dot{\varepsilon}(\mathbf{v}) = \frac{1}{2} (\nabla \mathbf{v} + (\nabla \mathbf{v})^T). \quad (14) $$
As suggested by Ganesan and Poirier [22] and Rappaz et al [20], we adopt the following model for the macroscopic deviatoric stress tensor:

$$ \Sigma^l = \langle \mathbf{s}^l \rangle = 2\mu_l g_l \operatorname{dev}(\dot{\boldsymbol{\varepsilon}}(\mathbf{v}_l)). \quad (15) $$

where dev(*) denotes the deviatoric part of tensor (*).

### 2.2.2. Macroscopic constitutive equation for the solid phase
Experimental studies on the behavior of metallic alloys at high temperature show that the response of the solid phase is well described by constitutive equations of the Norton-Hoff type [23, 24, 25]. Like the liquid, the solid is assumed to be incompressible at the microscopic scale, such that its response can be characterized at this scale by the following constitutive equation:

$$ \mathbf{s} = 2K_s \left(\sqrt{3}\dot{\varepsilon}_{eq}\right)^{m-1} \dot{\varepsilon}(\mathbf{v}), \quad (16) $$

where $K_s$ and $m$ denote the consistency and the strain rate sensitivity, respectively, and $\dot{\varepsilon}_{eq}$ the von Mises equivalent strain rate. The above equation can be rewritten in scalar form in terms of the equivalent stress (in the von-Mises sense) and the equivalent strain rate, as follows:

$$ \sigma_{eq} = 3K_s \left(\sqrt{3}\dot{\varepsilon}_{eq}\right)^{m-1} \dot{\varepsilon}_{eq}. \quad (17) $$

For solid fractions above the coherency fraction $g_{scohe}$, following the theoretical analysis of Geindreau and Auriault [18], the effective stress tensor

$$ \Sigma^s = \langle \mathbf{s}^s \rangle - \langle p^s \rangle \mathbf{I} + g_s p_l \mathbf{I} \quad (18) $$

is expressed as a degree $m$ homogeneous function with respect to the strain rate tensor $\langle \dot{\boldsymbol{\varepsilon}} \rangle^s = \dot{\boldsymbol{\varepsilon}}(\mathbf{v}^s)$. Note that such a result is valid as long as the viscosity of the liquid phase remains small compared to the solid phase apparent viscosity $\sigma_{eq}/\dot{\varepsilon}_{eq}$, that is to say for rather small values of the strain rate. Fortunately, this is generally the case for the solidification problems under consideration in this work. This property shows that the solid phase can be modelled as a compressible power-law fluid. We therefore adopt a compressible viscoplastic model [23, 12, 13]. Then, the constitutive equation coincides with equation (17), but with the equivalent stress and strain rate defined as

$$ \Sigma_{eq}^s = \left( A \operatorname{dev}(\Sigma^s) : \operatorname{dev}(\Sigma^s) + B \left( \operatorname{tr}(\Sigma^s) \right)^2 \right)^{1/2}, \quad (19) $$

$$ \langle \dot{\boldsymbol{\varepsilon}} \rangle_{eq}^s = \left( \frac{1}{A} \langle \dot{\boldsymbol{\varepsilon}} \rangle^s : \langle \dot{\boldsymbol{\varepsilon}} \rangle^s + \left( \frac{1}{9B} - \frac{1}{3A} \right) \operatorname{tr}(\langle \dot{\boldsymbol{\varepsilon}} \rangle^s)^2 \right)^{1/2}. \quad (20) $$

Those equations require two rheological functions $A$ and $B$ that depend on the solid volume fraction and for which several models can be found [23, 12, 26, 27]. The constitutive equations of the solid phase at the macroscopic scale then take the form

$$ \Sigma^s = 3K_s \left(\sqrt{3}\langle\dot{\boldsymbol{\varepsilon}}\rangle_{eq}^s\right)^{m-1} \left(\frac{1}{A}\langle\dot{\boldsymbol{\varepsilon}}\rangle^s + \left(\frac{1}{9B} - \frac{1}{3A}\right)\operatorname{tr}(\langle\dot{\boldsymbol{\varepsilon}}\rangle^s)\mathbf{I}\right), \quad (21) $$

or decomposed into deviatoric and spherical parts:

$$ \operatorname{dev}(\Sigma^s) = 3K_s \left(\sqrt{3}\langle\dot{\boldsymbol{\varepsilon}}\rangle_{eq}^s\right)^{m-1} \frac{1}{A} \operatorname{dev}(\langle\dot{\boldsymbol{\varepsilon}}\rangle^s), \quad (22) $$

$$ -\frac{1}{3} \operatorname{tr}(\Sigma^s) = -K_s \left(\sqrt{3}\langle\dot{\boldsymbol{\varepsilon}}\rangle_{eq}^s\right)^{m-1} \frac{1}{3B} \operatorname{tr}(\langle\dot{\boldsymbol{\varepsilon}}\rangle^s). \quad (23) $$
It is worth noticing that the preceding compressible viscoplastic constitutive model tends to the incompressible viscoplastic Norton-Hoff model as solidification completes. The rheological parameters A(g_s) and B(g_s) as defined in [23, 12, 13, 26, 27] satisfy A(1) = 3/2 and B(1) = 0, imposing in that way the incompressibility constraint by means of a penalty factor, the multiplier of tr(⟨ε⟩^s) in equation (23), which becomes infinite. Numerically, this factor is made arbitrarily high but finite, by replacing B(1) by B(g_s_{incomp}), where g_s_{incomp} is a threshold for the solid fraction, above which incompressible behavior is assumed. In R2SOL, we adopt g_s_{incomp} = 0.999.

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

### 2.2.3. Exchange of momentum
According to the work of Ni and Beckermann [6], $M_k$ can be partitioned as:

$$M_k = M_k^d + M_k^p, \qquad (24)$$

the first part being the contribution of deviatoric stresses and the second one the contribution of the isotropic part, and it can easily be shown that

$$M_s^d + M_l^d = 0 \quad \text{and} \quad M_s^p + M_l^p = 0. \qquad (25)$$

The liquid being a Newtonian incompressible fluid with a very low viscosity, we will assume that the pressure equilibrium in the liquid phase is almost instantaneous. Subsequently, the interfacial pressures [6] in both phases ($p_k^*$) equal the intrinsic average value of liquid pressure, i.e. its microscopic value:

$$p_l^* = p_s^* = p_l. \qquad (26)$$

Therefore, $M_k^p$ can be expressed as follows [6]:

$$M_s^p = -M_l^p = -p_l \nabla g_l = p_l \nabla g_s. \qquad (27)$$

On the other hand, depending on the solid fraction, the dissipative terms $M_k^d$ are generally interpreted either as the drag force exerted by the liquid on the isolated solid grains [6] or the filtration force exerted by the liquid flowing through the solid, considered as a rigid porous medium [20, 18]. In both cases, $M_k^d$ can be modelled by a law of the generic type:

$$M_s^d = -M_l^d = g_l C(g_l)(v_l - v_s). \qquad (28)$$

Assuming the macroscopic isotropy of the medium, the interaction tensor C can be reduced to a scalar coefficient C, which is usually defined as

$$C = \frac{g_l \mu_l}{\kappa}, \qquad (29)$$

being $\kappa$ the permeability of the solid matrix in the mushy zone.
The mushy zone is assimilated to an isotropic porous medium whose permeability is defined by the Carman-Kozeny formula

$$ \kappa = \frac{\lambda_2^2 g_l^3}{180(1 - g_l)^2}, \qquad (30) $$

where $\lambda_2$ is the secondary dendrite arms spacing.

### 2.2.4. Mechanical boundary conditions
The boundary conditions of the mechanical problem are:

$$ v_s = v_l = v_{imp} \qquad \text{on } \partial\Omega_u, \qquad (31) $$

$$ \langle \boldsymbol{\sigma}^s \rangle \boldsymbol{n} = \langle \boldsymbol{\sigma}^l \rangle \boldsymbol{n} = \boldsymbol{T}_{imp}, \qquad \text{on } \partial\Omega_s, \qquad (32) $$

where $\langle \boldsymbol{\sigma}^s \rangle = \boldsymbol{\Sigma}^s - g_s p_l \boldsymbol{I}$, $v_{imp}$ and $\boldsymbol{T}_{imp}$ are the imposed surface velocity and traction on the non-overlapped portions $\partial\Omega_u$ and $\partial\Omega_s$ of the boundary $\partial\Omega$ of the analysis domain $\Omega$; $\boldsymbol{n}$ is the normal unit vector pointing outwards $\partial\Omega$.

## 2.3. Energy conservation

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

$$ \boldsymbol{q} = -\lambda \nabla T. \qquad (33) $$

Our approach is based on the local thermal equilibrium assumption, which is valid because the liquid and the solid do exhibit neither too different thermal properties nor strong thermal interfacial barriers [29]. Therefore, at the macroscopic scale, their average temperature is the same:

$$ T_s = T_l = T. \qquad (34) $$

The enthalpies of both phases can then be written as a function of the temperature T as follows:

$$ h_s = \int_{T_0}^{T} c_p(\tau) \, d\tau \quad \text{and} \quad h_l = h_s + L, \qquad (35) $$

where $c_p$ and $L$ are the heat capacity and latent heat per unit mass, respectively. The average mixture enthalpy per unit mass takes the form

$$ \langle h \rangle = f_s h_s + f_l h_l = \int_{T_0}^{T} c_p(\tau) \, d\tau + f_l L, \qquad (36) $$

where $f_l$ and $f_s$ are respectively the liquid and solid mass fraction, satisfying $f_l + f_s = 1$ in a saturated solid-liquid medium. Volume and mass phase fractions are related by

$$ f_k = \frac{\rho_k g_k}{\langle \rho \rangle}, \qquad (37) $$

being

$$ \langle \rho \rangle = g_s \rho_s + g_l \rho_l \qquad (38) $$
the average density of the solid-liquid mixture.

By invoking the thermal equilibrium hypothesis (34), the energy conservation is written as a single equation which is the sum of the energy balance on both phases, where the terms of exchange $Q_s$ and $Q_l$ cancel themselves:

$$ \frac{\partial}{\partial t} \langle \rho h \rangle + \nabla \cdot \langle \rho h \mathbf{v} \rangle + \nabla \cdot \langle \mathbf{q} \rangle = 0, \qquad (39) $$

where

$$ \langle \rho h \rangle = \langle (\rho h)^s \rangle + \langle (\rho h)^l \rangle = g_s \rho_s h_s + g_l \rho_l h_l, \qquad (40) $$

$$ \langle \rho h \mathbf{v} \rangle = \langle (\rho h \mathbf{v})^s \rangle + \langle (\rho h \mathbf{v})^l \rangle = g_s \rho_s h_s \mathbf{v}_s + g_l \rho_l h_l \mathbf{v}_l, \qquad (41) $$

$$ \langle \mathbf{q} \rangle = \langle \mathbf{q}_s \rangle + \langle \mathbf{q}_l \rangle = -(g_s \lambda_s + g_l \lambda_l) \nabla T = -\langle \lambda \rangle \nabla T. \qquad (42) $$

In equation (42), the scalar form of the macroscopic conductivity was obtained one again using the isotropy assumption of the medium at the macroscopic level.

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

$$ \langle \rho h \rangle = \langle \rho \rangle \langle h \rangle, \qquad (43) $$

$$ \langle \rho h \mathbf{v} \rangle = \langle \rho \mathbf{v} \rangle \langle h \rangle + \langle \rho \rangle (\mathbf{v}_l - \mathbf{v}_s) f_s f_l L, \qquad (44) $$

where

$$ \langle \rho \mathbf{v} \rangle = g_s \rho_s \mathbf{v}_s + g_l \rho_l \mathbf{v}_l. \qquad (45) $$

After elementary calculations, we obtain the following advection-diffusion equation for the average mixture enthalpy:

$$ \langle \rho \rangle \frac{\partial \langle h \rangle}{\partial t} + \langle \rho \mathbf{v} \rangle \cdot \nabla \langle h \rangle - \nabla \cdot (\tilde{\kappa} \nabla \langle h \rangle) + S_T = 0, \qquad (46) $$

being

$$ \tilde{\kappa} = \langle \lambda \rangle \frac{\partial T}{\partial \langle h \rangle}, \qquad (47) $$

$$ S_T = L \nabla \cdot (\langle \rho \rangle (\mathbf{v}_l - \mathbf{v}_s) f_s f_l). \qquad (48) $$

As shown by the above equation, the source term $S_T$ is originated by the relative movement of both phases in the mushy zone.

### 2.3.1. Thermal initial and boundary conditions

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

$$ T = T_0 \quad \text{at } t = 0 \text{ in } \Omega, \qquad (49) $$

and the following boundary conditions:

$$ T = T_w \quad \text{at } t > 0 \text{ on } \partial \Omega_T, \qquad (50) $$

$$ (-\langle \lambda \rangle \nabla T) \cdot \mathbf{n} = q_w \quad \text{at } t > 0 \text{ on } \partial \Omega_q, \qquad (51) $$

$$ (-\langle \lambda \rangle \nabla T) \cdot \mathbf{n} = h(T - T_{\text{ext}}) \quad \text{at } t > 0 \text{ on } \partial \Omega_c, \qquad (52) $$

prescribing the temperature $T_w$ on $\partial \Omega_T$, the heat flux $q_w$ through $\partial \Omega_q$, and the heat exchange through $\partial \Omega_c$ due to convection to the environment at temperature $T_{\text{ext}}$ with $h$ as the convection coefficient; $\partial \Omega_T$, $\partial \Omega_q$, and $\partial \Omega_c$ are non-overlapping portions of the boundary $\partial \Omega$ of $\Omega$, being $\mathbf{n}$ the unit normal vector pointing outwards to $\partial \Omega$.
## 2.4. Solute conservation

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

$$j = -D\nabla w, \qquad (53)$$

being $D$ the solute diffusion coefficient. By summing the solute conservation equations averaged for each phase, and considering that $J_s = -J_l$, we obtain:

$$\frac{\partial \langle w \rangle}{\partial t} + \nabla \cdot \langle w \mathbf{v} \rangle + \nabla \cdot \langle \mathbf{j} \rangle = 0, \qquad (54)$$

where

$$\langle w \rangle = \langle w^s \rangle + \langle w^l \rangle = g_s w_s + g_l w_l, \qquad (55)$$

$$\langle w \mathbf{v} \rangle = \langle (w \mathbf{v})^s \rangle + \langle (w \mathbf{v})^l \rangle = g_s w_s \mathbf{v}_s + g_l w_l \mathbf{v}_l, \qquad (56)$$

$$\langle \mathbf{j} \rangle = \langle \mathbf{j}^s \rangle + \langle \mathbf{j}^l \rangle = -\langle (D \nabla w)^s \rangle - \langle (D \nabla w)^l \rangle = -g_s D_s \nabla w_s - g_l D_l \nabla w_l. \qquad (57)$$

The diffusion coefficients $D_k$ are very small (particularly $D_s$), so that they are usually neglected. However, in this work we neglect the diffusion in the solid phase, but maintain the diffusion in the liquid phase since it has a stabilizing effect on the numerical solution. Then:

$$\langle \mathbf{j} \rangle \approx -g_l D_l \nabla w_l. \qquad (58)$$

### 2.4.1. Solute exchange at the solid-liquid interface: the microsegregation model
The solute exchange between the solid and the liquid phase is governed by interfacial phenomena occurring at the scale of the dendrite arms. This is a non-equilibrium process involving a very complex physics. However, several simplifying hypotheses are allowed for macrosegregation studies. Besides the uniformity of temperature within the REV (already assumed for the energy balance), let us assume that both phases are at thermodynamic equilibrium at the interface, according to the alloy phase diagram. Then, at the interface, it satisfies:

$$w_s^* = k(T)w_l^* \quad \text{and} \quad w_l^* = \psi(T), \qquad (59)$$

where the thermodependent partition coefficient $k$ as well the equation of the liquidus line $\psi(T)$ can be determined from the phase diagram. Usually, both $k$ and the slope of the liquidus line ($m_l$) vary only slightly with temperature, letting us assume them constant, as shown in Figure 1.

Further, we assume the perfect diffusion of solute in the liquid phase at the microscopic scale, so that the solute concentration in the liquid inside the REV is uniform and particularly $w_l^* = w_l$. We obtain then:

$$w_s^* = kw_l \qquad (60)$$

and

$$w_l = \frac{T - T_m}{m_l}, \qquad (61)$$

where $T_m$ is the melting temperature of the pure substance.

Regarding solute diffusion in the solid phase, there are two classic limiting cases:
<figure><img src="image_2.png" /><figcaption>Figure 1. Local thermodynamic equilibrium at the interface, governed by the phase diagram (on the right), where both solidus and liquidus lines are linearized.</figcaption></figure>

* perfect diffusion in solid (lever rule), with the concentration in solid inside the REV uniform as well, and particularly $w_s^* = w_s$, so we can write:

$$ \langle w \rangle = (g_l + kg_s)w_l; \qquad (62) $$

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

$$ \langle w \rangle = g_l w_l + \int_0^{g_s} w_s^* dg_s = g_l w_l + \int_0^{g_s} kw_l dg_s, \qquad (63) $$

or, in differential form, by

$$ \frac{\partial \langle w \rangle}{\partial t} = \frac{\partial}{\partial t}(g_l w_l) + kw_l \frac{\partial g_s}{\partial t}. \qquad (64) $$

Besides these limiting cases, to account for a more realistic diffusion rate in the solid phase, one could solve the unsteady solute diffusion equation with moving interface in a one-dimensional domain representing the characteristic length of the microstructure and/or that of the diffusion mechanism we consider. We refer to the works of Poirier et al [30], Mo [31], Combeau et al [32] and more recently the one-dimensional finite-difference model of Thuinet et al [33], who solve the diffusion equation in the solid and the liquid phases for alloys with multiple constituents and including the $\delta/\gamma$ transformation in the case of steels.

### 2.4.2. Final form of the solute conservation equation
We can reduce the microsegregation model to a general microsegregation function relating the average concentration in the solid-liquid mixture to the liquid fraction $g_l$ and the average concentration in the liquid phase $w_l$:

$$ \langle w \rangle = \varphi(g_l, w_l), \qquad (65) $$

which is given by equations (62) and (63) for the particular cases of lever rule and Scheil law, respectively.

Now, we can write the solute balance equation in the form of an advection-diffusion equation for the average concentration, as follows:

$$ \frac{\partial \langle w \rangle}{\partial t} + v^w \cdot \nabla \langle w \rangle - \nabla \cdot (\tilde{D} \nabla \langle w \rangle) + S_w = 0, \qquad (66) $$
being

$$ \mathbf{v}^w = \mathbf{v}_s + g_l(\mathbf{v}_l - \mathbf{v}_s) \frac{\partial w_l}{\partial \langle w \rangle}, \qquad (67) $$

$$ \tilde{D} = g_l D_l \frac{\partial w_l}{\partial \langle w \rangle}, \qquad (68) $$

$$ S_w = \langle w \rangle \nabla \cdot \mathbf{v}_s + w_l \nabla \cdot (g_l(\mathbf{v}_l - \mathbf{v}_s)). \qquad (69) $$

In the same way as $S_T$ for the energy balance, the source term $S_w$ is originated by the relative movement and compressibility of both phases in the mushy zone.

### 2.4.3. Initial and boundary conditions for the solute problem
Equation (66) is subject to the initial condition

$$ \langle w \rangle = w_0 \quad \text{at } t = 0 \text{ in } \Omega, \qquad (70) $$

being $w_0$ the initial solute distribution, usually equal to the nominal solute content of the alloy.

Further, we assume that there is no solute flux through the domain boundary, such that the following natural boundary condition applies:

$$ \nabla \langle w \rangle \cdot \mathbf{n} = 0 \quad \text{at } t > 0 \text{ on } \partial\Omega. \qquad (71) $$

#### Final macroscopic balance equation set

$$ \text{Solid momentum:} \qquad \nabla \cdot \mathbf{\Sigma}^s - g_s \nabla p_l + g_l C(\mathbf{v}_l - \mathbf{v}_s) + g_s \rho_s \mathbf{g} = \mathbf{0} $$

$$ \text{Liquid momentum:} \qquad \nabla \cdot \mathbf{\Sigma}^l - g_l \nabla p_l - g_l C(\mathbf{v}_l - \mathbf{v}_s) + g_l \rho_l \mathbf{g} \\ = \rho_l \frac{\partial}{\partial t} (g_l \mathbf{v}_l) + \rho_l \nabla \cdot (g_l \mathbf{v}_l \times \mathbf{v}_l), $$

$$ \text{Mixture mass:} \qquad (1 - \Delta\epsilon^{tr})\nabla \cdot (g_s \mathbf{v}_s) + \nabla \cdot (g_l \mathbf{v}_l) = \frac{\partial g_s}{\partial t} \Delta\epsilon^{tr}, $$

$$ \text{Mixture energy:} \qquad \langle \rho \rangle \frac{\partial \langle h \rangle}{\partial t} + \langle \rho \mathbf{v} \rangle \cdot \nabla \langle h \rangle - \nabla \cdot (\tilde{\kappa} \nabla \langle h \rangle) + S_T = 0, $$

$$ \text{Mixture solute:} \qquad \frac{\partial \langle w \rangle}{\partial t} + \mathbf{v}^w \cdot \nabla \langle w \rangle - \nabla \cdot (\tilde{D} \nabla \langle w \rangle) + S_w = 0. $$

## 3. Numerical implementation

### 3.1. Mechanical model

This formulation was implemented in the two-dimensional code R2SOL, a finite element software using linear triangles and based on the $P1^+/P1$ mixed formulation [34, 35, 36]. These elements are shown in Figure 2.

As it is implicit in all the preceding formulation, the present model adopts the intrinsic average velocities of the solid and liquid phases as primal variables. This is remarked since in the isothermal version by Le Corre et al [37], which constitutes the starting point for the present
<figure><img src="image_3.png" /><figcaption>Figure 2. Triangular P1$^{+}$/P1 finite element.</figcaption></figure>

work, a different choice for the velocity in the liquid phase was made, adopting as primary unknown not the intrinsic average $v_l$ but the classic average $\langle v^l \rangle = g_l v_l$. This previous choice was motivated by the convenience of using the model already developed in CEMEF [36] for melt flow in a rigid porous medium, which has $\langle v^l \rangle$. However, for the solid phase, the use of the intrinsic value $v_s$ was compulsory in order to avoid the division by $g_l$ in the interaction term (28), source of numerical problems in the vicinity of the solidus interface.

The advantages of the new model are evident when dealing with the Dirichlet boundary condition (31), since these conditions are posed directly over the primal variables. It is not the case with the previous model [37], where this condition becomes dependent on the liquid fraction $g_l$.

Another advantage, in the authors' opinion, lies on the consistency of using the same interpolation functions for both intrinsic velocities $v_s$ and $v_l$, contrary to the work of Le Corre et al [37], where variables of different nature were identically interpolated.

Then, let $\mathcal{V}$ be the space of “kinematically admissible” velocity fields and $\mathcal{V}_0$ the space of “zero kinematically admissible” velocity fields. The virtual power principle states that the solution $(v_s, v_l, p_l) \in \mathcal{V} \times L^2(\Omega)$ of the problem must satisfy:

$$ \int_{\Omega} \boldsymbol{\Sigma}^s : \nabla v_s^* \, dV - \int_{\Omega} p_l \nabla \cdot (g_s v_s^*) \, dV - \int_{\Omega} g_l C(v_l - v_s) \cdot v_s^* \, dV \\ - \int_{\Omega} g_s \rho_s \boldsymbol{g} \cdot v_s^* \, dV = \int_{\partial \Omega_s} \langle \boldsymbol{\sigma}^s \rangle \boldsymbol{n} \cdot v_s^* \, dS, \qquad (72) $$

$$ \int_{\Omega} \boldsymbol{\Sigma}^l : \nabla v_l^* \, dV - \int_{\Omega} p_l \nabla \cdot (g_l v_l^*) \, dV + \int_{\Omega} g_l C(v_l - v_s) \cdot v_l^* \, dV \\ - \int_{\Omega} g_l \rho_l \boldsymbol{g} \cdot v_l^* \, dV + \int_{\Omega} \rho_l \gamma(g_l, v_l) \cdot v_l^* \, dV = \int_{\partial \Omega_s} \langle \boldsymbol{\sigma}^l \rangle \boldsymbol{n} \cdot v_l^* \, dS, \qquad (73) $$

$$ \int_{\Omega} \left( (1 - \Delta\epsilon^{tr})\nabla \cdot (g_s v_s) + \nabla \cdot (g_l v_l) - \frac{\partial g_s}{\partial t} \Delta\epsilon^{tr} \right) p_l^* \, dV = 0, \qquad (74) $$

for all $(v_s^*, v_l^*, p_l^*) \in \mathcal{V}_0 \times L^2(\Omega)$, with

$$ \gamma(g_l, v_l) = \frac{\partial}{\partial t}(g_l v_l) + \nabla \cdot (g_l v_l \times v_l) = g_l \frac{\partial_l v_l}{\partial t} + v_l \left( \frac{\partial g_l}{\partial t} + \nabla \cdot (g_l v_l) \right), \qquad (75) $$

where we introduce the “k-phase”-material derivative $(k = s, l)$ of any variable $\psi$ defined in the phase $k$ as

$$ \frac{d_k \psi}{dt} = \frac{\partial \psi}{\partial t} + v_k \cdot \nabla \psi. \qquad (76) $$
### 3.1.1. Finite element formulation of the mechanical problem
The above mixed solid velocity-liquid velocity-pressure formulation requires interpolation functions that satisfy the Babuška-Brezzi stability conditions [38]. The present finite element formulation will be called $(P1^+)^2/P1$, since it constitutes the extension to the case of two velocity fields of the $P1^+/P1$ formulation previously used for the “one-velocity” models [34, 35, 36]†.

On each finite element, the velocity fields $v_s$ and $v_l$ are approximated as follows:

$$v_s \approx \sum_{i=1}^{3} N_i(x) V_{s_i} + N_b(x) V_s^b = v_s^{P1} + b^s, \quad (77)$$

$$v_l \approx \sum_{i=1}^{3} N_i(x) V_{l_i} + N_b(x) V_l^b = v_l^{P1} + b^l, \quad (78)$$

where $V_{s_i}$ and $V_{l_i}$ ($i = 1, 2, 3$) are the unknown velocities of the solid and liquid phase at the vertex node $i$, $V_s^b$ and $V_l^b$ are the additional solid and liquid velocity unknowns associated to the bubble node located at the center of each element, $N_i$ is the standard linear (P1) shape function associated to the vertex node $i$, and $N_b$ is the bubble interpolation function, chosen here to be linear on each sub-triangle defined by the vertices and the center of the element. According to the definition of a bubble function, it is constructed such that $N_b = 0$ on the edges of the element and $N_b = 1$ at the bubble node.

The interstitial pressure $p_l$ is approximated by a standard linear interpolation:

$$p_l \approx \sum_{i=1}^{3} N_i(x) P_i, \quad (79)$$

where $P_i$ is the unknown pressure associated to node $i$.

Test functions $v_s^*$, $v_s^*$ and $p_l^*$ are approximated in the same way (this is the so-called Galerkin formulation [44], where shape and test functions coincide).

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

* Momentum in the solid phase:

<footnote>†“One-velocity” stands for models having only one velocity variable, such as the “one-phase” [34, 35, 39] and the “pseudo-two-phase” [40, 41, 42, 10, 43, 36] models. For the “one-phase” model, there are only two phases, the solid-like and the liquid-like (with the interface defined by the coherence temperature), and the velocity variable coincides with the velocity of one or the other phase. On the other hand, although the “pseudo-two-phase” considers the mushy zone as a mixture of solid and liquid phases, the solid phase is assumed to be rigid; in such a way, the liquid velocity remains the only velocity variable.</footnote>
– P1 contribution:

$$ \mathbf{R}_{v_s} = \mathbf{R}_{v_s v_s} + \mathbf{R}_{v_s v_l} + \mathbf{R}_{v_s b_s} + \mathbf{R}_{v_s b_l} + \mathbf{R}_{v_s p} - \mathbf{F}_{v_s} \quad (80) $$

$$ \mathbf{R}_{v_s v_s} = \int_{\Omega} \mathbf{B}^T \mathbf{\Sigma}^s (\mathbf{v}_s^{P1}) \, dV + \left( \int_{\Omega} g_l C \mathbf{N}^T \mathbf{N} \, dV \right) \mathbf{V}_s \quad (80a) $$

$$ \mathbf{R}_{v_s v_l} = \left( - \int_{\Omega} g_l C N^T N \, dV \right) \mathbf{V}_l \quad (80b) $$

$$ \mathbf{R}_{v_s b_s} = \int_{\Omega} \mathbf{B}^T \mathbf{\Sigma}^s (\mathbf{b}_s) \, dV + \left( \int_{\Omega} g_l C \mathbf{N}^T \mathbf{N}_b \, dV \right) \mathbf{V}_s^b \quad (80c) $$

$$ \mathbf{R}_{v_s b_l} = \left( - \int_{\Omega} g_l C \mathbf{N}^T \mathbf{N}_b \, dV \right) \mathbf{V}_l^b \quad (80d) $$

$$ \mathbf{R}_{v_s p} = \left( - \int_{\Omega} (g_s \mathbf{B}^T \mathbf{I} + \mathbf{N}^T \nabla g_s) \mathbf{N}_p \, dV \right) \mathbf{P} \quad (80e) $$

$$ \mathbf{F}_{v_s} = \int_{\Omega} g_s \rho_s \mathbf{N}^T \mathbf{g} \, dV + \int_{\partial \Omega_{\sigma}} \mathbf{N}^T \mathbf{T}_{imp} \, dS \quad (80f) $$

– Bubble contribution:

$$ \mathbf{R}_{b_s} = \mathbf{R}_{b_s v_s} + \mathbf{R}_{b_s v_l} + \mathbf{R}_{b_s b_s} + \mathbf{R}_{b_s b_l} + \mathbf{R}_{b_s p} - \mathbf{F}_{b_s} \quad (81) $$

$$ \mathbf{R}_{b_s v_s} = \int_{\Omega} \mathbf{B}_b^T \mathbf{\Sigma}^s (\mathbf{v}_s^{P1}) \, dV + \left( \int_{\Omega} g_l C \mathbf{N}_b^T \mathbf{N} \, dV \right) \mathbf{V}_s \quad (81a) $$

$$ \mathbf{R}_{b_s v_l} = \left( - \int_{\Omega} g_l C \mathbf{N}_b^T \mathbf{N} \, dV \right) \mathbf{V}_l \quad (81b) $$

$$ \mathbf{R}_{b_s b_s} = \int_{\Omega} \mathbf{B}_b^T \mathbf{\Sigma}^s (\mathbf{b}_s) \, dV + \left( \int_{\Omega} g_l C \mathbf{N}_b^T \mathbf{N}_b \, dV \right) \mathbf{V}_s^b \quad (81c) $$

$$ \mathbf{R}_{b_s b_l} = \left( - \int_{\Omega} g_l C \mathbf{N}_b^T \mathbf{N}_b \, dV \right) \mathbf{V}_l^b \quad (81d) $$

$$ \mathbf{R}_{b_s p} = \left( - \int_{\Omega} (g_s \mathbf{B}_b^T \mathbf{I} + \mathbf{N}_b \nabla g_s) \mathbf{N}_p \, dV \right) \mathbf{P} \quad (81e) $$

$$ \mathbf{F}_{b_s} = \int_{\Omega} g_s \rho_s \mathbf{N}_b \mathbf{g} \, dV \quad (81f) $$

• Momentum in the liquid phase:
– P1 contribution:

$$ \mathbf{R}_{v_l} = \mathbf{R}_{v_l v_s} + \mathbf{R}_{v_l v_l} + \mathbf{R}_{v_l b_s} + \mathbf{R}_{v_l b_l} + \mathbf{R}_{v_l p} - \mathbf{F}_{v_l} \quad (82) $$

$$ \mathbf{R}_{v_l v_s} = \left( - \int_{\Omega} g_l C \mathbf{N}^T \mathbf{N} \, \mathrm{d}V \right) \mathbf{V}_s \quad (82a) $$

$$ \mathbf{R}_{v_l v_l} = \int_{\Omega} \mathbf{B}^T \boldsymbol{\Sigma}^l (\mathbf{v}_l^{P1}) \, \mathrm{d}V + \left( \int_{\Omega} g_l C \mathbf{N}^T \mathbf{N} \, \mathrm{d}V \right) \mathbf{V}_l + \underbrace{\int_{\Omega} \mathbf{N}^T \rho_l \gamma(g_l, \mathbf{v}_l^{P1}) \, \mathrm{d}V}_{\mathbf{R}_{v_l v_l}^{\gamma}} \quad (82b) $$

$$ \mathbf{R}_{v_l b_s} = \left( - \int_{\Omega} g_l C \mathbf{N}^T \mathbf{N}_b \, \mathrm{d}V \right) \mathbf{V}_s^b \quad (82c) $$

$$ \mathbf{R}_{v_l b_l} = \int_{\Omega} \mathbf{B}^T \boldsymbol{\Sigma}^l (\mathbf{b}_l) \, \mathrm{d}V + \left( \int_{\Omega} g_l C \mathbf{N}^T \mathbf{N}_b \, \mathrm{d}V \right) \mathbf{V}_l^b + \underbrace{\left( \int_{\Omega} \mathbf{N}^T \rho_l \gamma(g_l, \mathbf{b}_l) \, \mathrm{d}V \right)}_{\mathbf{R}_{v_l b_l}^{\gamma}} \quad (82d) $$

$$ \mathbf{R}_{v_l p} = \left( - \int_{\Omega} (g_l \mathbf{B}^T \mathbf{I} + \mathbf{N} \nabla g_l) \mathbf{N}_p \, \mathrm{d}V \right) \mathbf{P} \quad (82e) $$

$$ \mathbf{F}_{v_l} = \int_{\Omega} g_l \rho_l \mathbf{N}^T \mathbf{g} \, \mathrm{d}V + \int_{\partial \Omega_{\sigma}} \mathbf{N}^T \mathbf{T}_{imp} \, \mathrm{d}S \quad (82f) $$

– Bubble contribution:

$$ \mathbf{R}_{b_l} = \mathbf{R}_{b_l v_s} + \mathbf{R}_{b_l v_l} + \mathbf{R}_{b_l b_s} + \mathbf{R}_{b_l b_l} + \mathbf{R}_{b_l p} - \mathbf{F}_{b_l} \quad (83) $$

$$ \mathbf{R}_{b_l v_s} = \left( - \int_{\Omega} g_l C \mathbf{N}_b^T \mathbf{N} \, \mathrm{d}V \right) \mathbf{V}_s \quad (83a) $$

$$ \mathbf{R}_{b_l v_l} = \int_{\Omega} \mathbf{B}_b^T \boldsymbol{\Sigma}^l (\mathbf{v}_l^{P1}) \, \mathrm{d}V + \left( \int_{\Omega} g_l C \mathbf{N}_b^T \mathbf{N} \, \mathrm{d}V \right) \mathbf{V}_l + \underbrace{\left( \int_{\Omega} \mathbf{N}_b^T \rho_l \gamma(g_l, \mathbf{v}_l^{P1}) \, \mathrm{d}V \right)}_{\mathbf{R}_{b_l v_l}^{\gamma}} \quad (83b) $$

$$ \mathbf{R}_{b_l b_s} = \left( - \int_{\Omega} g_l C \mathbf{N}_b^T \mathbf{N}_b \, \mathrm{d}V \right) \mathbf{V}_s^b \quad (83c) $$
$$ \mathbf{R}_{b_l b_l} = \int_{\Omega} \mathbf{B}_b^T \boldsymbol{\Sigma}^l(\mathbf{b}_l) \, \mathrm{d}V + \left( \int_{\Omega} g_l C \mathbf{N}_b^T \mathbf{N}_b \, \mathrm{d}V \right) \mathbf{V}_l^b + \underbrace{\left( \int_{\Omega} \mathbf{N}_b^T \rho_l \boldsymbol{\gamma}(g_l, \mathbf{b}_l) \, \mathrm{d}V \right)}_{\mathbf{R}_{b_l b_l}^{\gamma}} \quad (83d) $$

$$ \mathbf{R}_{b_l p} = \left( - \int_{\Omega} (g_l \mathbf{B}_b^T \mathbf{I} + \mathbf{N}_b \nabla g_l) \mathbf{N}_p \, \mathrm{d}V \right) \mathbf{P} \quad (83e) $$

$$ \mathbf{F}_{b_l} = \int_{\Omega} g_l \rho_l \mathbf{N}_b \mathbf{g} \, \mathrm{d}S \quad (83f) $$

• Mass in the solid-liquid mixture:

$$ \mathbf{R}_p = \mathbf{R}_{pv_s} + \mathbf{R}_{pv_l} + \mathbf{R}_{pb_s} + \mathbf{R}_{pb_l} + \mathbf{R}_{pp} - \mathbf{F}_p \quad (84) $$

$$ \mathbf{R}_{pv_s} = \left( \int_{\Omega} (\Delta\varepsilon^{tr} - 1) \mathbf{N}_p (g_s \mathbf{B}^T \mathbf{I} + \mathbf{N} \nabla g_s)^T \, \mathrm{d}V \right) \mathbf{V}_s \quad (84a) $$

$$ \mathbf{R}_{pv_l} = \left( - \int_{\Omega} \mathbf{N}_p (g_l \mathbf{B}^T \mathbf{I} + \mathbf{N} \nabla g_l)^T \, \mathrm{d}V \right) \mathbf{V}_l \quad (84b) $$

$$ \mathbf{R}_{pb_s} = \left( \int_{\Omega} (\Delta\varepsilon^{tr} - 1) \mathbf{N}_p (g_s \mathbf{B}_b^T \mathbf{I} + \mathbf{N}_b \nabla g_s)^T \, \mathrm{d}V \right) \mathbf{V}_s^b \quad (84c) $$

$$ \mathbf{R}_{pb_l} = \left( - \int_{\Omega} \mathbf{N}_p (g_l \mathbf{B}_b^T \mathbf{I} + \mathbf{N}_b \nabla g_l)^T \, \mathrm{d}V \right) \mathbf{V}_l^b \quad (84d) $$

$$ \mathbf{R}_{pp} = \mathbf{0} \quad (84e) $$

$$ \mathbf{F}_p = - \int_{\Omega} \Delta\varepsilon^{tr} \frac{\partial g_s}{\partial t} \mathbf{N}_p \mathbf{g} \, \mathrm{d}V \quad (84f) $$

In the above expressions, $N$, $N_b$ and $N_p$ are respectively the matrices of linear velocity, bubble velocity and pressure shape functions, $B$, $B_b$ and $B_p$ the corresponding gradient matrices, and $V_s$, $V_l$ and $P$ are the vectors of nodal solid velocity, liquid velocity and pressure unknowns; further, the tensors such as $\Sigma^s$, $\Sigma^l$ and $\mathbf{I}$ are mapped into vectors as usual in the finite element practice (see [44] for details).

### 3.1.2. Treatment of inertia and advection effects
For $g_s > g_{scohe}$, the porosity of the solid matrix is relatively low, slowing the interstitial melt flow. In that case, it is licit to neglect the contributions to advection and inertia from the bubble degree of freedom:

$$ \mathbf{R}_{b_l v_l}^{\gamma} = \mathbf{R}_{b_l b_l}^{\gamma} = \mathbf{0}. \quad (85) $$

Consequently, considering the previous hypothesis, the $\gamma$-term contribution reduces to

$$ \mathbf{R}_{v_l v_l}^{\gamma} = \int_{\Omega} \mathbf{N}^T \rho_l g_l \frac{\partial_l \mathbf{v}_l}{\partial t} \, \mathrm{d}V + \int_{\Omega} \mathbf{N}^T \rho_l \mathbf{v}_l \left( \frac{\mathrm{d}g_l}{\mathrm{d}t} + \nabla \cdot (g_l \mathbf{v}_l) \right) \, \mathrm{d}V, \quad (86) $$

Time discretization is achieved via the fully implicit (then, unconditionally stable) backward Euler method. Let $t = t_0 + \Delta t$ be the current instant, $t_0$ the previous instant, and $\Delta t$ the time step size. Then, the time derivative of any variable $\psi$ at time $t$ is approximated as follows

$$ \frac{\partial \psi}{\partial t}(\mathbf{x}, t) \approx \frac{\psi(\mathbf{x}, t) - \psi(\mathbf{x}, t_0)}{\Delta t}, \quad (87) $$
being $x$ the spatial coordinates with respect to a fixed reference frame.

In the computation of time derivatives, the relative movement between the material and spatial mesh must be taken into account. This is achieved by using the Arbitrary Lagrangian-Eulerian (ALE) method we proposed in [45]. First, let us introduce the concept of mesh derivative (i.e., the derivative along the nodal trajectories [46]), related to the k-phase-material derivative by the expression:

$$ \frac{d_k \psi}{dt} = \frac{\partial_{\text{mesh}} \psi}{\partial t} + (\mathbf{v}_k - \mathbf{v}_{\text{mesh}}) \cdot \nabla \psi. \qquad (88) $$

where $\mathbf{v}_{\text{mesh}}$ is the mesh velocity field. In this work, we assume the mesh velocity equal to the intrinsic solid velocity, i.e.

$$ \mathbf{v}_{\text{mesh}} \equiv \mathbf{v}_s. \qquad (89) $$

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

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

$$ \mathbf{R}(\mathbf{V}^{(k)}) \approx \mathbf{R}(\mathbf{V}^{(k-1)}) + \mathbf{K} \Delta \mathbf{V} = 0, \qquad (90) $$

where $\mathbf{V}$ is the vector of generic unknowns and $\mathbf{K} = d\mathbf{R}/d\mathbf{V}$ is the tangent matrix.

In this work, convergence is further improved by using a line-search algorithm [47]. Actually, let us assume that $\Delta\mathbf{V}$ defines not the increment but the search direction, in such a way that the unknown variable is updated to the k-th iteration as follows:

$$ \mathbf{V}^{(k)} = \mathbf{V}^{(k-1)} + \alpha \Delta \mathbf{V}, \qquad (91) $$

where the parameter $\alpha$ is determined by requiring the orthogonality between the residual vector and the search direction, i.e.

$$ \mathbf{R}(\mathbf{V}^{(k-1)} + \alpha \Delta \mathbf{V}) \cdot \Delta \mathbf{V} = 0. \qquad (92) $$

This scalar equation is solved using the simple regula-falsi method detailed in [44].

Particularized to the nonlinear system of equations (80) to (84), the Newton-Raphson method implies to solve at each iteration the following linear equation system:

$$ \begin{bmatrix} \mathbf{K}_{v_s v_s} & \mathbf{K}_{v_s v_l} & \mathbf{K}_{v_s b_s} & \mathbf{K}_{v_s b_l} & \mathbf{K}_{v_s p} \\ \mathbf{K}_{v_l v_s} & \mathbf{K}_{v_l v_l} & \mathbf{K}_{v_l b_s} & \mathbf{K}_{v_l b_l} & \mathbf{K}_{v_l p} \\ \mathbf{K}_{b_s v_s} & \mathbf{K}_{b_s v_l} & \mathbf{K}_{b_s b_s} & \mathbf{K}_{b_s b_l} & \mathbf{K}_{b_s p} \\ \mathbf{K}_{b_l v_s} & \mathbf{K}_{b_l v_l} & \mathbf{K}_{b_l b_s} & \mathbf{K}_{b_l b_l} & \mathbf{K}_{b_l p} \\ \mathbf{K}_{p v_s} & \mathbf{K}_{p v_l} & \mathbf{K}_{p b_s} & \mathbf{K}_{p b_l} & \mathbf{K}_{p p} \end{bmatrix} \begin{bmatrix} \Delta \mathbf{V}_s \\ \Delta \mathbf{V}_l \\ \Delta \mathbf{V}_s^b \\ \Delta \mathbf{V}_l^b \\ \Delta \mathbf{P} \end{bmatrix} = - \begin{bmatrix} \mathbf{R}_{v_s} \\ \mathbf{R}_{v_l} \\ \mathbf{R}_{b_s} \\ \mathbf{R}_{b_l} \\ \mathbf{R}_p \end{bmatrix}, \qquad (93) $$
where

$$
\begin{align}
\boldsymbol{K}_{v_s v_s} &= \frac{\mathrm{d}\boldsymbol{R}_{v_s}}{\mathrm{d}\boldsymbol{V}_s}, & \boldsymbol{K}_{v_s v_l} &= \frac{\mathrm{d}\boldsymbol{R}_{v_s}}{\mathrm{d}\boldsymbol{V}_l}, & \boldsymbol{K}_{v_s b_s} &= \frac{\mathrm{d}\boldsymbol{R}_{v_s}}{\mathrm{d}\boldsymbol{V}_s^b}, \nonumber \\
\boldsymbol{K}_{v_s b_l} &= \frac{\mathrm{d}\boldsymbol{R}_{v_s}}{\mathrm{d}\boldsymbol{V}_l^b}, & \boldsymbol{K}_{v_s p} &= \frac{\mathrm{d}\boldsymbol{R}_{v_s}}{\mathrm{d}\boldsymbol{P}}, & \text{etc.} \tag{94}
\end{align}
$$

For a space of dimension *d* (*d* = 2 for 2D and axisymmetric cases), there are *n*<sub>*n*</sub> = 2*d* + 1 unknowns associated to the nodes and *n*<sub>*e*</sub> = 2*d* (bubble) unknowns associated to the elements of the spatial mesh, that is one more vectorial unknown per node and per element with respect to the one-velocity problem. This corresponds either to the additional liquid velocity field regarding the solid mechanics model [34, 35] or the solid velocity considering the fluid dynamics work [36]. However, by rewriting equation (93) in the form

$$
\begin{bmatrix} \mathbf{K}_{vv} & \mathbf{K}_{vb} & \mathbf{K}_{vp} \\ \mathbf{K}_{bv} & \mathbf{K}_{bb} & \mathbf{K}_{bp} \\ \mathbf{K}_{pv} & \mathbf{K}_{pb} & \mathbf{K}_{pp} \end{bmatrix} \begin{bmatrix} \Delta V \\ \Delta V^b \\ \Delta P \end{bmatrix} = - \begin{bmatrix} \mathbf{R}_v \\ \mathbf{R}_b \\ \mathbf{R}_p \end{bmatrix}, \quad (95)
$$

we realize that, despite the supplementary unknowns, it remains identical in form to that of one-velocity problems.

Furthermore, since the bubble unknowns are internal to the element, they can be eliminated at the element level by the same procedure, known as static condensation, used for one-velocity problems [34, 35]. From equation (95) applied to each element, the iterative bubble increment results

$$
\Delta \mathbf{V}^b = -\mathbf{K}_{bb}^{-1} (\mathbf{R}_b + \mathbf{K}_{bv} \Delta \mathbf{V} + \mathbf{K}_{bp} \Delta \mathbf{P}). \quad (96)
$$

Then, the system (95) reduces to

$$
\begin{bmatrix} \mathbf{K}_{vv} - \mathbf{K}_{vb}\mathbf{K}_{bb}^{-1}\mathbf{K}_{bv} & \mathbf{K}_{vp} - \mathbf{K}_{vb}\mathbf{K}_{bb}^{-1}\mathbf{K}_{bp} \\ \mathbf{K}_{pv} - \mathbf{K}_{pb}\mathbf{K}_{bb}^{-1}\mathbf{K}_{bv} & \mathbf{K}_{pp} - \mathbf{K}_{pb}\mathbf{K}_{bb}^{-1}\mathbf{K}_{bp} \end{bmatrix} \begin{bmatrix} \Delta V \\ \Delta P \end{bmatrix} = \begin{bmatrix} \mathbf{K}_{vb}\mathbf{K}_{bb}^{-1}\mathbf{R}_b - \mathbf{R}_v \\ \mathbf{K}_{pb}\mathbf{K}_{bb}^{-1}\mathbf{R}_b - \mathbf{R}_p \end{bmatrix}. \quad (97)
$$

The resulting linear and symmetric system of equations is solved using the PETSc library [48], which allows parallel processing.

3.1.4. Transport of the solution in case of moving meshes
By solving the mechanical problem
defined by the equation system (97), we obtain an estimate to the liquid-phase-material
derivative of the velocity in the liquid phase:

$$
\frac{d_l \boldsymbol{v}_l}{dt} \approx \frac{\boldsymbol{v}_l(\boldsymbol{x} + \boldsymbol{v}_l \Delta t, t) - \boldsymbol{v}_l(\boldsymbol{x}, t_0)}{\Delta t}, \quad (98)
$$

from which an estimate to the mesh derivative is computed:

$$
\frac{\partial_{\text{mesh}} \boldsymbol{v}_l}{\partial t} = \frac{d_l \boldsymbol{v}_l}{d t} - (\boldsymbol{v}_l - \boldsymbol{v}_{\text{mesh}}) \cdot \nabla \boldsymbol{v}_l. \quad (99)
$$

Finally, the time-discretized form of the mesh derivative,

$$
\frac{\partial_{\text{mesh}} \boldsymbol{v}_l}{\partial t} \approx \frac{\boldsymbol{v}_l(\boldsymbol{x} + \boldsymbol{v}_{\text{mesh}} \Delta t, t) - \boldsymbol{v}_l(\boldsymbol{x}, t_0)}{\Delta t}, \quad (100)
$$
determines the updated value of the liquid velocity:

$$v_l(\mathbf{x} + v_{\text{mesh}}\Delta t, t) = v_l(\mathbf{x}, t_0) + \frac{\partial_{\text{mesh}}v_l}{\partial t}\Delta t. \qquad (101)$$

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

## 3.2. Thermal model

### 3.2.1. Treatment of the mesh movement
As done for the mechanical problem, an ALE method allowing the description of the relative movement between the solid and the liquid regions and the mesh must be used. However, as demonstrated by Gouttebroze et al [43], the performance of the Streamline-Upwind Petrov/Galekin (SUPG [49]) for solving the heat balance equation is superior to that of the ALE method previously described for the mechanical problem. In order to apply the SUPG formulation, we introduce the mesh derivative of the average enthalpy

$$\frac{\partial_{\text{mesh}}\langle h \rangle}{\partial t} = \frac{\partial\langle h \rangle}{\partial t} + v_{\text{mesh}} \cdot \nabla\langle h \rangle, \qquad (102)$$

into the heat balance equation (46), obtaining

$$\langle \rho \rangle \frac{\partial_{\text{mesh}} \langle h \rangle}{\partial t} + \langle \rho \rangle \tilde{v} \cdot \nabla \langle h \rangle - \nabla \cdot (\tilde{\kappa} \nabla \langle h \rangle) + S_T = 0, \qquad (103)$$

where

$$\langle \rho \rangle \tilde{v} = \langle \rho v \rangle - \langle \rho \rangle v_{\text{mesh}}. \qquad (104)$$

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

$$\frac{\partial_{\text{mesh}}\langle h \rangle}{\partial t} \approx \frac{\langle h \rangle (\mathbf{x} + v_{\text{mesh}}\Delta t, t) - \langle h \rangle (\mathbf{x}, t_0)}{\Delta t}, \qquad (105)$$

Then, the solution of equation (103) directly yields the updated value of the average enthalpy at the new nodal positions $\mathbf{x} + v_{\text{mesh}}\Delta t$, and, unlike mechanical analysis, no transport step is needed.

### 3.2.2. Finite element formulation of the thermal problem
The solution of the thermal problem consists in finding $\langle h \rangle \in L^2(\Omega)$ satisfying:

$$\int_{\Omega} \langle \rho \rangle \frac{\partial_{\text{mesh}} \langle h \rangle}{\partial t} \omega \, dV + \int_{\Omega} \langle \rho \rangle \tilde{v} \cdot \nabla \langle h \rangle \omega \, dV + \int_{\Omega} \tilde{\kappa} \nabla \langle h \rangle \cdot \nabla \omega \, dV \\ + \int_{\Omega} S_T \omega \, dV - \int_{\partial \Omega_q} q_w \omega \, dS - \int_{\partial \Omega_c} h(T - T_{\text{ext}}) \omega \, dS = 0, \qquad (106)$$

for all $\omega \in L^2(\Omega)$ such that $\omega = 0$ on $\partial \Omega_T$.
As usual, the solution is approximated inside each linear triangular finite element as

$$ \langle h \rangle \approx \sum_{i=1}^{3} N_i(\boldsymbol{x})H_i, \qquad (107) $$

where $H_i$ is the unknown value of the average enthalpy at node $i$ located at the vertex $i$ of each element, $N_i$ is the standard linear shape function associated to the node $i$, the same as used for the mechanical problem.

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

$$ \omega = N_i + \alpha(\tilde{\boldsymbol{v}} \cdot \nabla N_i) = \tilde{N}_i, \qquad (108) $$

where $\tilde{\boldsymbol{v}}$ is evaluated at the center of the element and $\alpha$ is the upwind coefficient, defined as [49]

$$ \alpha = \frac{\Delta x}{\sqrt{15}\|\tilde{\boldsymbol{v}}\|} \left( \frac{1}{\tanh Pe} - \frac{1}{Pe} \right). \qquad (109) $$

where $\| \cdot \|$ denotes the $L_2$-norm of $(\cdot)$,

$$ Pe = \frac{\|\tilde{\boldsymbol{v}}\| \Delta x}{2\tilde{\kappa}/\langle \rho \rangle} \qquad (110) $$

is the element Peclet number, and $\Delta x$ is a characteristic size of the element in the advection direction, defined as [50]

$$ \Delta x = \frac{2\|\tilde{\boldsymbol{v}}\|}{\sum_{i=1}^{3} |\tilde{\boldsymbol{v}} \cdot \nabla N_i|}. \qquad (111) $$

Introducing the trial and weighting functions defined by equations (107) and (108) into the weak form (106) of the heat balance, we obtain the discrete heat equation:

$$ \mathbf{R}_T = M \frac{\partial_{\text{mesh}} \mathbf{H}}{\partial t} + (\mathbf{C} + \mathbf{K}) \mathbf{H} + \mathbf{S}(\mathbf{H}) + \mathbf{c}(\mathbf{H}) - \mathbf{F} = \mathbf{0}, \qquad (112) $$
the components of the above arrays being given by:

$$M_{ij} = \int_{\Omega} \langle \rho \rangle \tilde{N}_i N_j \, dV, \qquad (113)$$

$$C_{ij} = \int_{\Omega} \tilde{N}_i (\langle \rho \rangle \tilde{v} \cdot \nabla N_j) \, dV, \qquad (114)$$

$$K_{ij} = \int_{\Omega} \tilde{\kappa} \nabla \tilde{N}_i \cdot \nabla N_j \, dV, \qquad (115)$$

$$S_i = \int_{\Omega} S_T(\langle h \rangle) \tilde{N}_i \, dV, \qquad (116)$$

$$c_i = \int_{\partial \Omega_c} h T(\langle h \rangle) N_i \, dS, \qquad (117)$$

$$F_i = \int_{\partial \Omega_q} q_w N_i \, dS + \int_{\partial \Omega_c} h T_{ext} N_i \, dS. \qquad (118)$$

In the above, the existence of a function $T = T(\langle h \rangle)$, reciprocal to that given by the enthalpy definition (36), is assumed. Further, for a P1 element, $\nabla \tilde{N}_i \equiv \nabla N_i$.

The nonlinear and non-symmetric heat equation is solved using the Newton-Raphson technique together with line-searching, in a way analogous to that described for the mechanical problem in Section 3.1.3. Once again, at each Newton-Raphson iteration, linear solvers from the PETSc library [48] are used.

### 3.3. Macrosegregation model

Since an analogy can be made between energy and solute balance equations, (46) and (66) respectively, the latter is solved in an identical way to that described in the previous section.

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

$$\frac{\partial_{mesh} \langle w \rangle}{\partial t} \approx \frac{\langle w \rangle (\boldsymbol{x} + \boldsymbol{v}_{mesh} \Delta t, t) - \langle w \rangle (\boldsymbol{x}, t_0)}{\Delta t}, \qquad (119)$$

into the solute balance equation (66), we have

$$\frac{\partial_{mesh} \langle w \rangle}{\partial t} + \tilde{v}^w \cdot \nabla \langle w \rangle - \nabla \cdot (\tilde{D} \nabla \langle w \rangle) + S_w = 0, \qquad (120)$$

with

$$\tilde{v}^w = v^w - v_{mesh}. \qquad (121)$$

It is worth noting here that taking $v_{mesh} = v_s$, the solid intrinsic velocity, the time derivative in equation (120) appears as the rate of variation of the average concentration when following the displacement of the solid phase. Equation (120) is then very close to the equation obtained by Lesoult et al [51], who developed a macrosegregation equation during the deformation of a coherent mushy zone using a Lagrangian approach.

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

$$\langle w \rangle \approx \sum_{i=1}^{3} N_i(\boldsymbol{x}) W_i, \qquad (122)$$
being $W_i$ the unknown value at node $i$. Further, a SUPG weighting function is used, identical to that defined by equation (108) with

$$ \alpha = \frac{\Delta x}{\sqrt{15} \|\tilde{v}^w\|} \left( \frac{1}{\tanh Pe} - \frac{1}{Pe} \right), \qquad (123) $$

$$ Pe = \frac{\|\tilde{v}^w\| \Delta x}{2\tilde{D}}, \qquad (124) $$

$$ \Delta x = \frac{2 \|\tilde{v}^w\|}{\sum_{i=1}^{3} |\tilde{v}^w \cdot \nabla N_i|}. \qquad (125) $$

Now, the discrete solute balance equation can be expressed as

$$ \mathbf{R}_w = \mathbf{M}^w \frac{\partial_{\text{mesh}} \mathbf{W}}{\partial t} + (\mathbf{C}^w + \mathbf{K}^w) \mathbf{W} + \mathbf{S}^w(\mathbf{W}) = \mathbf{0}, \qquad (126) $$

where

$$ M_{ij}^w = \int_{\Omega} \tilde{N}_i N_j \, dV, \qquad (127) $$

$$ C_{ij}^w = \int_{\Omega} \tilde{N}_i (\tilde{v}^w \cdot \nabla N_j) \, dV, \qquad (128) $$

$$ K_{ij}^w = \int_{\Omega} \tilde{D} \nabla N_i \cdot \nabla N_j \, dV, \qquad (129) $$

$$ S_i^w = \int_{\Omega} S_w(\langle w \rangle) \tilde{N}_i \, dV. \qquad (130) $$

Finally, equation (126) is solved here again using the Newton-Raphson method with line-searching, using at each iteration the linear solvers from the PETSc library.

## 4. Coupled solution of thermal, mechanical, and macrosegregation problems

### 4.1. Intermediate thermo-macrosegregation coupling

Thermal and macrosegregation solutions interact due to the simultaneous influence of both enthalpy and solute concentration on the local liquid or solid fraction. Assuming $\langle h \rangle$ and $\langle w \rangle$ to be known from the solution of thermal and macrosegregation problems, respectively, we invoke the enthalpy dependency on $g_l$ and $T$ defined by equation (36), the relationship (61) between $w_l$ and $T$ derived from the equilibrium phase diagram, and the generic microsegregation rule (65) relating relationship $\langle w \rangle$ to $w_l$ and $g_l$, in order to pose a system of scalar equations from which $w_l$, $g_l$ and $T$ can be computed. The so-defined “microsegregation box” can be expressed as follows:

$$ \mathbf{m}(w_l, T, g_l; \langle h \rangle, \langle w \rangle) = \begin{bmatrix} \langle h \rangle - \int_{T_0}^T c_p(\tau) \, d\tau - g_l L \\ T - T_m - m_l w_l \\ \langle w \rangle - \varphi(g_l, w_l) \end{bmatrix} = \mathbf{0}. \qquad (131) $$

Henceforth, by writing $(\cdot; \bullet)$, it is implied that $(\cdot)$ are variables and $(\bullet)$ are parameters.
This system has a closed analytical solution just in a few particular cases, for instance, when the microsegregation is described by the lever rule (62) and the heat capacity $c_p$ is assumed to be constant. Otherwise, it is solved numerically using a regula-falsi method, spanning between $g_l = 0$ at the solidus temperature $T_{sol}$ (or $g_l = g_{l_{eut}}$ at the beginning of the eutectic transformation at temperature $T_{eut}$) and $g_l = 1$ at the liquidus temperature $T_{liq}$.

Above system also determines the derivatives $\partial T/\partial \langle h \rangle$ and $\partial w_l/\partial \langle w \rangle$. As before, these derivatives have analytical expressions for a few cases, and are in general computed numerically.

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

## 4.2. Global coupling scheme

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

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

* **New time instant $t \leftarrow t + \Delta t$:**

- **Initialization:**

$$k = 0,$$

$$ (\mathbf{v}_s, \mathbf{v}_l, p_l, \langle h \rangle, \langle w \rangle; T, w_l, g_l)_{t+\Delta t}^{(0)} = (\mathbf{v}_s, \mathbf{v}_l, p_l, \langle h \rangle, \langle w \rangle; T, w_l, g_l)_t $$

- **New coupling iteration $k \leftarrow k + 1$:**

1. Obtain $\langle h \rangle^{(k)}$ by solving the thermal problem

$$ \mathbf{R}_T(\langle h \rangle^{(k)}; \langle w \rangle^{(k-1)}, \mathbf{v}_s^{(k-1)}, \mathbf{v}_l^{(k-1)}) = \mathbf{0}. $$

2. Obtain $\langle w \rangle^{(k)}$ by solving the macrosegregation problem

$$ \mathbf{R}_w(\langle w \rangle^{(k)}; \langle h \rangle^{(k)}, \mathbf{v}_s^{(k-1)}, \mathbf{v}_l^{(k-1)}) = \mathbf{0}. $$

3. Obtain $w_l^{(k)}$, $T^{(k)}$, and $g_l^{(k)}$ by calling the “microsegregation box”:

$$ \mathbf{m}(w_l^{(k)}, T^{(k)}, g_l^{(k)}; \langle h \rangle^{(k)}, \langle w^{(k)} \rangle) = \mathbf{0}. $$

4. Obtain $\mathbf{v}_s^{(k)}$, $\mathbf{v}_l^{(k)}$, and $p_l^{(k)}$ by solving the mechanical problem

$$ \mathbf{R}_M(\mathbf{v}_s^{(k)}, \mathbf{v}_l^{(k)}, p_l^{(k)}; w_l^{(k)}, T^{(k)}, g_l^{(k)}) = \mathbf{0}. $$

5. Check convergence:

* If $\|(*)^{(k)} - (*)^{(k-1)}\| > \text{tolerance}$, go to the next coupling iteration.
* Otherwise, update the variables

$$ (\mathbf{v}_s, \mathbf{v}_l, p_l, \langle h \rangle, \langle w \rangle; T, w_l, g_l)_{t+\Delta t} = (\mathbf{v}_s, \mathbf{v}_l, p_l, \langle h \rangle, \langle w \rangle; T, w_l, g_l)_{t+\Delta t}^{(k)} $$

Then, go to the next time step.

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

In addition, it is worth noting that the “microsegregation box” is also called inside the thermal and macrosegregation solvers at each Newton-Raphson iteration, not only for computing $g_l$, $T$, and $w_l$, but also the derivatives $\partial T / \partial \langle h \rangle$ and $\partial w_l / \partial \langle w \rangle$, as mentioned in Section 4.1.

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

## 5. Application tests

### 5.1. Validation: Axial compression of a semi-solid specimen

Let us consider a simple compression test on a porous specimen saturated with liquid, illustrated in Figure 3. This problem has a closed-form solution discussed in the Appendix, and serves then for the validation of the finite element implementation of the mechanical model previously presented. The following application aims to be representative of the mechanical conditions encountered in the mushy core of a continuously cast steel slab at the end of the secondary cooling stage when passing along the support rolls. The values of the considered geometrical and material parameters are given in Table I.

<table><caption>Table I. Dimensions and mechanical parameters used in the validation test.</caption><tr><td>$l = 0.25$ m</td><td>$K_s = 10^7$ Pa s</td><td>$\mu_l = 10^{-3}$ Pa s</td></tr><tr><td>$h = 0.1$ m</td><td>$A = 20$</td><td>$\lambda_2 = 100$ to $200$ μm</td></tr><tr><td>$\|\mathbf{v}_{imp}\| = 1$ mm/s</td><td>$B = 0.5$</td><td></td></tr><tr><td>$g_s = 0.8$</td><td></td><td></td></tr></table>

Figure 4 shows the influence of the permeability

$$ \kappa = \frac{\lambda_2^2 g_l^3}{180(1 - g_l)^2} $$
<figure><img src="image_4.png" /><figcaption>Figure 3. Two-phase simple compression test.</figcaption></figure>

as it varies due to a variation of the secondary dendrite arm spacing $\lambda_2$. As observed, the lower this spacing (and the permeability as well), the slower the interdendritic flow (characterized Figure 4 by the relative velocity $v_l - v_s$) and the higher the interstitial pressure $p_l$.

The results of the two-phase finite element model match exactly the closed form solution, as illustrated in Figure 5 for $v_l - v_s$ and $p_l$.

<figure><img src="image_5.png" /><figcaption>Figure 4. Analytical solutions for the relative velocity $v_l - v_s$ and the interstitial pressure $p_l$ showing the influence of the secondary dendrite arm spacing $\lambda_2$.</figcaption></figure>

## 5.2. Application to liquid segregation

The second test aims to simulate a situation frequently encountered in continuous casting problems. It consists in imposing a deformation on a reduced portion of the surface of a partially solidified alloy (Figure 6). It thus simulates the action of a roll. Symmetry conditions are imposed on all the boundaries except on the upper face, which is a free surface. We
<figure><img src="image_6.png" /><figcaption>Figure 5. Analytical vs. two-phase finite element solution for the relative velocity $v_l - v_s$ and the interstitial pressure $p_l$ ($\lambda_2 = 150$ µm).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Redistribution of the liquid inside a semi-solid domain with an initial solid fraction distribution. The vectors on the right represent the average relative velocity $\langle v_l \rangle - g_l v_s$.</figcaption></figure>

perform here again an isothermal calculation, but imposing an initial distribution of liquid fraction along the y-axis described in Figure 6. Here, the solid phase rheology is still somewhat arbitrary, but it is now non-linear, with a strain rate sensitivity $m = 0.2$. The interaction coefficient is modelled by a classical Carman-Kozeny law and reaches a very high value in the solidified zone.

Black arrows plotted on the deformed domain on the right Figure 6 show the relative average velocity vectors $\langle v_l \rangle - g_l v_s$. It can be seen that the applied pressure leads to a redistribution of the liquid with respect to the solid. In the zone where the solid is under pressure (left side of the deformed sample), it undergoes densification and the liquid is rejected to the right side of the sample.
5.3. Application to the late stage of solidification in steel continuous casting

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

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

This problem has already been addressed by several authors. Following the pioneer work of Miyazawa and Schwerdtfeger [52], and the analytical developments of Lesoult and Sella [53], Kajitani et al. [54] proposed the more advanced model. However, they did not take into account the deformation of the solid phase in the mushy zone in their approach, and made strong hypotheses on its velocity field, which was deduced directly and homothetically from the velocity field of the solid crust. The present two-phase formulation aims to improve the modelling of macrosegregation by considering not only the relative but also the concurrent movement between the solid and the liquid phases.

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

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

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

Up to this two-phase section, a classical “one-phase” analysis [36] is performed, which is used as an initialization for the two-phase analysis. In Figure 7 we observe thermal results from the “one-phase” analysis, particularly concerning the liquid fraction distribution along the strand. The enthalpy distribution at the initial two-phase section (located 15.2 m-far from the meniscus) serves as initial condition for the two-phase thermal analysis. Since the solid fraction must be high enough to let us assume the solid phase to be coherent, we eliminate the higher values of enthalpy at the core of the slab.
<table><caption>Table II. Material data for continuous casting test case I.</caption><tbody><tr><td><strong>Alloy</strong></td><td>Fe-0.1%C</td></tr><tr><td>ρ<sub>s</sub> = ρ<sub>l</sub></td><td>7060. kg/m<sup>3</sup></td></tr><tr><td>c<sub>p</sub></td><td>675. J/(kg K)</td></tr><tr><td>λ<sub>s</sub> = λ<sub>l</sub></td><td>30. W/(m K)</td></tr><tr><td>L</td><td>2.6×10<sup>5</sup> J/kg</td></tr><tr><td>T<sub>m</sub></td><td>1538.°C</td></tr><tr><td>m<sub>l</sub></td><td>-80. K/wt%C</td></tr><tr><td>k</td><td>0.18</td></tr><tr><td>λ<sub>2</sub></td><td>1.×10<sup>-4</sup> m</td></tr><tr><td>Microsegregation law</td><td>Lever rule</td></tr><tr><td><strong>Solid phase</strong></td><td>Compressible viscoplastic [12]</td></tr><tr><td>K(T), m(T)</td><td>from Kozlowski et al [55] (K(T) intentionally decreased)</td></tr><tr><td>A(g<sub>s</sub>), B(g<sub>s</sub>)</td><td>from Nguyen et al [12]</td></tr><tr><td>g<sub>scohe</sub></td><td>0.65</td></tr><tr><td><strong>Liquid phase</strong></td><td>Newtonian</td></tr><tr><td>μ<sub>l</sub></td><td>0.001 Pa s</td></tr></tbody></table>

Further, the current version of the “one-phase” solver [36] is not able up to date to perform a macrosegregation analysis under continuous casting conditions, so we will assume a uniform solute concentration $w_0$ (equal to the nominal carbon content) in the initial two-phase section as initial condition. Although it implies an important simplification, this hypothesis allows us to distinguish the effect of solid deformation from the other sources of macrosegregation that are relevant in the earlier stages of the continuous casting process (for instance, macrosegregation by natural convection driven by thermal and solutal gradients).

In order to use the global unsteady-state approach [35, 39], a buffer zone is defined, consisting of a thick slice of the slab that ends at the initial two-phase section, as shown in Figure 8. Unlike in the one-phase analysis, where the buffer zone is located above the meniscus in the mold region, the buffer zone now lies in the curved portion of the caster. As before, the buffer length must be greater than the maximum displacement per time step of any particle in the initial two-phase section, otherwise it could be impossible to find the particle antecedent in the advection calculations.

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

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

Figure 9 plots the region with positive and negative values of the divergence of the intrinsic
<figure><img src="image_8.png" /><figcaption>Figure 7. Test case I. One-phase analysis of steel continuous casting process. On the left, distribution of liquid fraction along the slab. On the top, the enthalpy distribution at the initial two-phase section serving as initial thermal condition for the two-phase analysis.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Test case I. Two-phase analysis of steel continuous casting process: finite element mesh for the buffer zone and mechanical boundary conditions.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9. Test case I. Final distribution of the negative (blue or dark grey) and positive (red or light grey) values of the divergence of the intrinsic velocity in the solid phase, $\nabla \cdot v_s$.</figcaption></figure>

velocity in the solid phase ($\nabla \cdot v_s$), indicating respectively expansion and contraction of the solid skeleton in the mushy zone. As it can be seen there, the solid phase is subject to alternating expansion and contraction at the core of the slab. Expansion usually develops in the intervals between two successive rolls, and is associated to the local influx of liquid melt caused by the bulging of the slab. Conversely, contraction normally develops under the rolls, due the local expulsion of liquid.

Figure 10 evidences the complexity of the melt flow in the mushy region. Between two successive rolls, the bulging of the solidified shell has two main effects: the development of tensions in the mushy core of the slab, and concurrently, of compressions in the periphery of this region. The latter effect can also be explained by assimilating the longitudinal section of the solidified shell as a beam supported by the rolls, whose fibers adjacent to the solidification front suffers compression. This combined central expansion-peripheral compression causes the relative velocity $v_l - v_s$ point to the core, implying the liquid enrichment of this region.

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

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

The axial shift between the zones of positive $\nabla \cdot v_s$ and those of positive segregations is probably due to a delay which corresponds to the transport of solute, taking into account the difference between the solid and the liquid velocity and the general casting velocity. Indeed,
<figure><img src="image_11.png" /><figcaption>Figure 10. Test case I. Detail of the distribution of relative velocity vectors ($v_l - v_s$) in the mushy zone.</figcaption></figure>

the solute transport is not instantaneous.

Further, major macrosegregation defects are shown to appear at the center of the slab, with a slighter correlate of opposite sign in its periphery. In real continuous casting processes, only positive central macrosegregation is expected. The spurious alternate sign of central macrosegregation observed here may be explained by the excessive bulging.

Let us also notice the existence of some oscillations of the numerical solution for the solute distribution, particularly in the periphery of the mushy core. In fact, as shown by Figure 13, the mesh density decreases from the core, where the mesh is adequately refined as proved by the smoothness of the solution, to the outside, where it is evidently insufficient. This problem should be overcome when the new adaptive mesh refinement tools developed by Liu [36] in R2SOL for one-phase analysis be operational for the current two-phase analysis. This would allow to keep an adequate refinement of the continuously evolving mushy region at each time step.

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

### 5.3.2. Test case II
The second application intends to assume more realistic bulging conditions in order to get closer to the real physics of a continuously cast steel slab, although we must
<figure><img src="image_12.png" /><figcaption>Figure 11. Test case I. Effect of the contraction (blue or dark grey) and the expansion (red or light grey) of the solid skeleton on the relative velocity $v_l - v_s$ in the mushy zone.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12. Test case I. Regions of positive (red or light grey) and negative (blue or dark grey) carbon segregation $\langle w \rangle - w_0$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13. Test case I. Influence of the finite element mesh on the smoothness of the macrosegregation solution.</figcaption></figure>

keep in mind the limitations imposed by the adopted hypotheses, particularly the coherency of the solid phase and the lack of hardening effect in the constitutive equation for the solid phase.

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

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

Regarding macrosegregation, the nominal carbon content $w_0 = 0.18\%$ is set as initial condition inside the buffer zone. The final carbon segregation $\langle w \rangle - w_0$ is depicted in Figure 16. We observe there the ability of the current model to capture the phenomenon of central
<figure><img src="image_15.png" /><figcaption>Figure 14. Test case II. One-phase analysis of steel continuous casting process. On the left, distribution of liquid fraction along the slab. On the top, the enthalpy distribution at the initial two-phase section serving as initial thermal condition for the two-phase analysis.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. Test case II. Final distribution of liquid fraction in the analyzed portion of the slab.</figcaption></figure>
<table><caption>Table III. Material data for continuous casting test case II.</caption><tbody><tr><td><strong>Alloy</strong></td><td>Fe-0.18%C</td></tr><tr><td>ρ<sub>s</sub> = ρ<sub>l</sub></td><td>7300. kg/m<sup>3</sup></td></tr><tr><td>c<sub>p</sub></td><td>675. J/(kg K)</td></tr><tr><td>λ<sub>s</sub> = λ<sub>l</sub></td><td>35. W/(m K)</td></tr><tr><td>L</td><td>2.6×10<sup>5</sup> J/kg</td></tr><tr><td>T<sub>m</sub></td><td>1538.°C</td></tr><tr><td>m<sub>l</sub></td><td>-80. K/wt%C</td></tr><tr><td>k</td><td>0.18</td></tr><tr><td>λ<sub>2</sub></td><td>1.×10<sup>-4</sup> m</td></tr><tr><td>Microsegregation law</td><td>Lever rule</td></tr><tr><td><strong>Solid phase</strong></td><td>Compressible viscoplastic [12]</td></tr><tr><td>K(T), m(T)</td><td>from Kozlowski et al [55]</td></tr><tr><td>A(g<sub>s</sub>), B(g<sub>s</sub>)</td><td>from Nguyen et al [12]</td></tr><tr><td>g<sub>s</sub><sub>cohe</sub></td><td>0.65</td></tr><tr><td><strong>Liquid phase</strong></td><td>Newtonian</td></tr><tr><td>μ<sub>l</sub></td><td>0.001 Pa s</td></tr></tbody></table>

(positive) segregation, characteristic of continuously cast steel slabs [56].

Even if the solidification is not completed so that the solute concentrations can still evolve, the computed level of macrosegregation seems to be significantly lower than that observed in reality. Nevertheless, we feel that this could be strongly improved by the use of more realistic constitutive models and parameters, especially for the solid phase compressible behavior and the mechanical interaction term. Furthermore, as already mentioned, a finer mesh would be necessary to confirm such trends.

## 6. Conclusions and perspectives

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

However, further research is needed before obtaining a robust computational tool with the accuracy required in real casting applications. First, although the viscoplastic compressible model [23, 12, 13] seems to be adequate to characterize the macroscopic mechanical behavior of the coherent solid phase in the mushy zone, this is contested by the lack of experimental (and
<figure><img src="image_17.png" /><figcaption>Figure 16. Test case II. Final distribution of carbon segregation ($(w) - w_0$) in the analyzed portion of the slab.</figcaption></figure>

literature) data. Further, this model does not account for strain hardening, a mechanism that becomes important once the material is completely solidified. It is not only isotropic hardening that must be modelled, as done in “one-phase” models, but also kinematic hardening in order to account for alternate loading.

Concerning the viscoplastic incompressibility constraint in the fully solid region, the current two-phase model enforces it by means of a penalty factor for the volume change rate tr($\dot{\varepsilon}_s$) into the constitutive equation for the solid phase, while in the “one-phase” model [35, 39] it is implied in the mass conservation equation. This difference blocks the use of the available “one-phase” mechanical model [35, 39], which is already capable of modelling isotropic hardening, for the fully solid region in two-phase analysis. The development of a library of constitutive models, suitable for one-phase and two-phase mechanical analysis, should be encouraged as future work.

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

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

# APPENDIX

## Analytical solution of a simple compression test on a porous specimen

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

$$ \Sigma^s = K_s (\beta \langle \dot{\varepsilon} \rangle^s + \alpha \operatorname{tr}(\langle \dot{\varepsilon} \rangle^s) \mathbf{I}), \quad (132) $$

where

$$ \beta = \frac{3}{A}, \quad \alpha = \frac{1}{3B} - \frac{1}{A}. \quad (133) $$

The liquid behavior is described by equation (15) [22, 20], linear as well, rewritten here as

$$ \Sigma^l = g_l \mu_l (b \langle \dot{\varepsilon} \rangle^l + a \operatorname{tr}(\langle \dot{\varepsilon} \rangle^l) \mathbf{I}), \quad (134) $$

where

$$ b = 2, \quad a = -\frac{2}{3}. \quad (135) $$

Now, let us consider the simple compression test illustrated in Figure 3, adopting the following assumptions:

* plane strain conditions;

* the liquid fraction is uniform throughout the specimen;

* gravity and inertia effects are neglected;

* no mass transfer is considered between the two phases;

* the normal stress on each phase is supposed to be null at $x = \pm l$.

We start with the following kinematic hypothesis:

$$ \mathbf{v}_s(x) = \begin{bmatrix} u(x) \\ -\dot{\varepsilon}y \end{bmatrix}, \quad \mathbf{v}_l(x) = \begin{bmatrix} v(x) \\ -\dot{\varepsilon}y \end{bmatrix}, \quad (136) $$

where $\dot{\varepsilon} = \|\mathbf{v}_{imp}\|/h$ denotes the nominal compression strain rate. The strain rate tensors in both phases have then the following expression:

$$ \langle \dot{\varepsilon} \rangle^s = \dot{\varepsilon}(\mathbf{v}_s) = \begin{bmatrix} \frac{dv}{dx} & 0 \\ 0 & -\dot{\varepsilon} \end{bmatrix}, \quad \langle \dot{\varepsilon} \rangle^l = \dot{\varepsilon}(\mathbf{v}_l) = \begin{bmatrix} \frac{dv}{dx} & 0 \\ 0 & -\dot{\varepsilon} \end{bmatrix}, \quad (137) $$
Applying momentum and mass conservation equations, we get:

$$ \nabla \cdot \boldsymbol{\Sigma}^s - g_s \nabla p_l - \boldsymbol{M}_l^d = \mathbf{0}, \qquad (138) $$

$$ \nabla \cdot \boldsymbol{\Sigma}^l - g_l \nabla p_l + \boldsymbol{M}_l^d = \mathbf{0}, \qquad (139) $$

$$ \nabla \cdot (g_s \boldsymbol{v}_s) + \nabla \cdot (g_l \boldsymbol{v}_l) = 0, \qquad (140) $$

with

$$ \boldsymbol{M}_l^d = -\frac{g_l^2 \mu_l}{\kappa} (\boldsymbol{v}_l - \boldsymbol{v}_s). \qquad (141) $$

Replacing the constitutive equations (132) and (134) as well as kinematic hypothesis (136) and (137) into the above conservation equations, we obtain

$$ K_s(\alpha + \beta) \frac{d^2 u}{dx^2} - g_s \frac{\partial p_l}{\partial x} + \frac{g_l^2 \mu_l}{\kappa} (v - u) = 0, \qquad (142) $$

$$ \mu_l g_l (a + b) \frac{d^2 v}{dx^2} - g_l \frac{\partial p_l}{\partial x} - \frac{g_l^2 \mu_l}{\kappa} (v - u) = 0, \qquad (143) $$

$$ \frac{\partial p_l}{\partial y} = 0 \quad \Rightarrow \quad p_l = p_l(x), \qquad (144) $$

$$ g_s \frac{du}{dx} + g_l \frac{dv}{dx} - \dot{\varepsilon} = 0. \qquad (145) $$

Summing the two momentum conservation equations (142) and (143) and using then the derivative of the mass balance equation (145) with respect to x, we have

$$ \frac{\partial p_l}{\partial x} = (K_s(\alpha + \beta) - g_s \mu_l (a + b)) \frac{d^2 u}{dx^2}, \qquad (146) $$

which is integrated to give

$$ p_l - p_0 = (K_s(\alpha + \beta) - g_s \mu_l (a + b)) \frac{du}{dx}. \qquad (147) $$

The integration constant $p_0$ is determined by the free surface conditions at $x = l$:

$$ \langle \sigma_{xx}^s \rangle = \Sigma_{xx}^s - g_s p_l = 0, \qquad (148) $$

$$ \langle \sigma_{xx}^l \rangle = \Sigma_{xx}^l - g_l p_l = 0. \qquad (149) $$

By introducing the constitutive laws (132) and (134), the mass balance equation (145), and equation (147) into the previous boundary conditions, after a few operation we get

$$ \tilde{\mu} \left. \frac{du}{dx} \right|_{x=l} - \alpha K_s \dot{\varepsilon} = g_s p_0, \qquad (150) $$

$$ - \tilde{\mu} \left. \frac{du}{dx} \right|_{x=l} + \mu_l \dot{\varepsilon} (ag_s + b) = g_l p_0, \qquad (151) $$

with

$$ \tilde{\mu} = g_l K_s (\alpha + \beta) + g_s^2 \mu_l (a + b). \qquad (152) $$
Now, the sum of equations (150) and (151) yields

$$p_0 = (\mu_l(ag_s + b) - \alpha K_s) \dot{\varepsilon}. \qquad (153)$$

On the other hand, injecting equation (146) together with the integral form of the mass balance equation (145) given by

$$g_s u(x) + g_l v(x) - \dot{\varepsilon}x = 0, \qquad (154)$$

into the solid momentum equation (142), we arrive to

$$\tilde{\mu} \frac{d^2 u}{dx^2} + \frac{g_l \mu_l}{\kappa} (\dot{\varepsilon}x - u) = 0, \qquad (155)$$

which, subject to the boundary condition $u(0) = 0$, has the following solution

$$u(x) = \dot{\varepsilon}x + a_0 \sinh(rx), \qquad (156)$$

where

$$r = \sqrt{\frac{g_l \mu_l}{\kappa \tilde{\mu}}}. \qquad (157)$$

We determine the integration constant $a_0$ by replacing this expression for $u$ into equation (150) and using the expression (153) for $p_0$. This leads to

$$a_0 = \frac{g_l(bg_s\mu_l - \beta K_s)}{\tilde{\mu}r \cosh(rL)} \dot{\varepsilon}. \qquad (158)$$

Finally, from the integral mass balance (154), we get

$$v(x) = \dot{\varepsilon}x - \frac{g_s}{g_l} a_0 \sinh(rx). \qquad (159)$$

Summarizing, the full solution of this problem is given by

$$\begin{aligned} v_{s_x}(x) &= \dot{\varepsilon}x + a_0 \sinh(rx), \\ v_{l_x}(x) &= \dot{\varepsilon}x - \frac{g_s}{g_l} a_0 \sinh(rx), \\ p_l(x) &= p_0 + (K_s(\alpha + \beta) - g_s \mu_l(a+b)) (\dot{\varepsilon} + a_0 r \cosh(rx)), \end{aligned}$$

with

$$\begin{aligned} r &= \sqrt{\frac{g_l \mu_l}{\kappa \tilde{\mu}}}, \\ \tilde{\mu} &= g_l K_s (\alpha + \beta) + g_s^2 \mu_l (a + b), \\ a_0 &= \frac{g_l(bg_s \mu_l - \beta K_s)}{\tilde{\mu}r \cosh(rL)} \dot{\varepsilon}, \\ p_0 &= (\mu_l(ag_s + b) - \alpha K_s) \dot{\varepsilon}. \end{aligned}$$

REFERENCES
1. M. C. Flemings. Our understanding of macrosegregation: past and present. *ISIJ International*, 40(9):833-841, 2000.

2. A. Einstein. Eine neue bestimmung der moleküldimensionen. *Annalen der Physik*, 19:289-306, 1906. See also corrective note: *Annalen der Physik*, 34:591-592, 1911.

3. P. Kumar, C. L. Martin, and S. B. Brown. Constitutive modelling and characterisation of the flow behavior of semi-solid metal alloy slurries. *Acta Metall. Mater.*, 42:3595-3614, 1994.

4. A. Zavaliangos and A. Lawley. Numerical simulation of thixoforming. *J. Materials Engrg. Performance*, 4:40-47, 1995.

5. W. D. Bennon and F. P. Incropera. A continuum model for momentum, heat and species transport in binary solid-liquid phase change systems. *I. Model formulation. Int. J. Heat Mass Transfer*, 30(10):2161-2170, 1987.

6. J. Ni and C. Beckermann. A volume-averaged two-phase model for transport phenomena during solidification. *Metall. Trans. B*, 22:349-361, 1991.

7. C. Y. Wang and C. Beckermann. Equiaxed dendritic solidification with convection: Part 1. multiscale/multiphase modeling. *Metall. and Mat. Trans. A*, 27:2754-2764, 1996.

8. N. Ahmad. *Numerical simulation of transport processes in multicomponent systems related to solidification problems*. PhD thesis, Ecole Polytechnique Fédérale de Lausanne, Switzerland, 1995.

9. J.-L. Desbiolles, Ph. Thevoz, and M. Rappaz. Micro-/macrosegregation modeling in casting: a fully coupled 3D model. In D. M. Stefanescu, J. A. Warren, M. R. Jolly, and M. J. M. Krane, editors, *Proc. MCWASP X, 10th Int. Conf. on Modeling of Casting, Welding and Advanced Solidification Processes*, pages 242-252, Destin, Florida, USA, May 2003.

10. M. Bellet, V. D. Fachinotti, S. Gouttebroze, W. Liu, and H. Combeau. A 3D-FEM model solving thermomechanics and macrosegregation in binary alloys solidification. In M. Rappaz, C. Beckermann, and R. Trivedi, editors, *Proc. Symposium on Solidification Processes and Microstructures, in Honor of Wilfried Kurz, TMS Annual Meeting*, pages 41-46, Charlotte, NC, USA, March 2004.

11. L. A. Lalli. A model for deformation and segregation of solid-liquid mixtures. *Metall. Trans. A*, 16:1393-1403, 1985.

12. T. G. Nguyen, D. Favier, and M. Suery. Theoretical and experimental study of the isothermal mechanical behaviour of alloys in the semi-solid state. *Int. J. Plasticity*, 10:663-693, 1994.

13. C. L. Martin, D. Favier, and M. Suéry. Viscoplastic behaviour of porous metallic materials saturated with liquid. part i: Constitutive equations. *Int. J. Plasticity*, 13(3):215-235, 1997.

14. S. Toyoshima and Y. Takahashi. A numerical simulation of forming processes for semi-solid materials. *ISIJ Int.*, 31:577-582, 1991.

15. F. Bay, J. Barlier, Y. Chastel, M. Bobadilla, and G. Lovato. Numerical modelling of extrusion for a semi-solid alloy. application to rheological parameters identification in Sn-Pb alloys. In B. G. Thomas and C. Beckermann, editors, *Proc. MCWASP VIII, 8th Int. Conf. on Modeling of Casting, Welding and Advanced Solidification Processes*, pages 941-948, San Diego, California, USA, June 1998.

16. M. M'Hamdi, A. Mo, and C. L. Martin. Two-phase modeling directed toward hot tearing formation in aluminum direct chill casting. *Metall. and Mat. Trans. A*, 33:2081-2093, 2002.

17. L. C. Nicolli, A. Mo, and M. M'Hamdi. Modeling of macrosegregation caused by volumetric deformation in a coherent mushy zone. *Metall. Mater. Trans. A*, 36:433-442, 2005.

18. C. Geindreau and J.-L. Auriault. Investigation of the viscoplastic behaviour of alloys in the semi-solid state by homogenization. *Mechanics of Materials*, 31:535-551, 1999.

19. M. Hassanizadeh and W. G. Gray. General conservation equations for multi-phase systems: 1. averaging procedure. *Adv. in Water Resources*, 2:131-144, 1979.

20. M. Rappaz, M. Bellet, and M. Deville. *Numerical modelling in materials science and engineering*. Springer Verlag, New York, 2003.

21. J.-L. Auriault and E. Sanchez-Palencia. Etude du comportement macroscopique d'un milieu poreux saturé déformable. *J. de Mécanique*, 16:575-603, 1977.

22. S. Ganesan and D. R. Poirier. Conservation of mass and momentum for the flow of interdendritic liquid during solidification. *Metall. Trans. B*, 21:173-181, 1990.

23. M. Abouaf, J.-L. Chenot, G. Raisson, and P. Bauduin. Finite element simulation of hot isostatic pressing of metal powders. *Int. J. Numer. Methods Engrg.*, 25:191-212, 1988.

24. C. Geindreau, D. Bouvard, and P. Doremus. Constitutive behaviour of metal powder during hot forming. Part I: Experimental investigation with lead powder as a simulation material. *Eur. J. Mech. A/Solids*, 18:581-596, 1999.

25. C. Geindreau, D. Bouvard, and P. Doremus. Constitutive behaviour of metal powder during hot forming. Part II: Unified viscoplastic modeling. *Eur. J. Mech. A/Solids*, 18:597-615, 1999.

26. C. L. Martin, D. Favier, and M. Suéry. Viscoplastic behaviour of porous metallic materials saturated with liquid. part ii: Experimental identification on a Sn-Pb model alloy. *Int. J. Plasticity*, 13(3):237-259, 1997.

27. C. L. Martin, D. Favier, and M. Suéry. Fracture behaviour in tension of viscoplastic porous metallic
materials saturated with liquid. *Int. J. Plasticity*, 15:981-1008, 1999.

28. O. Ludwig, B. Commet, J. M. Drezet, C. L. Martin, and M. Suéry. Rheological behavior of partially solidified al-cu alloys : experimental and numerical study. In D. M. Stefanescu, J. A. Warren, M. R. Jolly, and M. J. M. Krane, editors, *Proc. MCWASP X, 10th Int. Conf. on Modeling of Casting, Welding and Advanced Solidification Processes*, pages 183-190, Destin, Florida, USA, May 2003.

29. J.-L. Auriault and H. I. Ene. Macroscopic modelling of heat transfer in composites with interfacial thermal barrier. *Int. J. Heat Mass Transfer*, 37(18):2885-2892, 1994.

30. D. R. Poirier, P. J. Nandapurkar, and S. Ganesan. The energy and solute conservation equations for dendritic solidification. *Metall. Trans. B*, 22:889-900, 1991.

31. A. Mo. An internal variable description of solidification suitable for macrosegregation modeling. *Metall. and Mat. Trans. B*, 25:597-605, 1994.

32. H. Combeau, J.-M. Drezet, A. Mo, and M. Rappaz. Modeling of microsegregation in macrosegregation computations. *Metall. and Mat. Trans. A*, 27:2314-2327, 1996.

33. L. Thuinet, G. Lesoult, and H. Combeau. Computer simulation of microsegregation in the case of columnar growth involving a peritectic transformation for multicomponent steels. In D. M. Stefanescu, J. A. Warren, M. R. Jolly, and M. J. M. Krane, editors, *Proc. MCWASP X, 10th Int. Conf. on Modeling of Casting, Welding and Advanced Solidification Processes*, pages 237-244, Destin, Florida, USA, May 2003.

34. E. Perchat, L. Fourment, and T. Coupez. Mixed formulation and iterative solver for the parallelisation of a forging simulation software by mesh partitioning. In B. H. V. Topping, editor, *3rd Euro Conf. on Parallel and Distributed Computing for Computational Mechanics*, pages 67-72, Weimar, Germany, 1999.

35. A. Heinrich. *Modélisation thermomécanique de la coulée continue d'acier en deux dimensions*. PhD thesis, Ecole des Mines de Paris, France, 2003.

36. W. Liu. *Finite element modelling of macrosegregation and thermomechanical phenomena in solidification processes*. PhD thesis, Ecole des Mines de Paris, France, 2005.

37. S. Le Corre and M. Bellet. Two-phase modeling of metals solidification: a numerical approach for the thermo-mechanical problem. In S. Ghosh, J. M. Castro, and J. K. Lee, editors, *Proc. NUMIFORM'2004, 8th Int. Conf. on Numerical Methods in Industrial Forming Processes*, pages 1185-1190, Columbus, Ohio, USA, June 2004.

38. D. N. Arnold. Mixed finite element methods for elliptic problems. *Comput. Methods Appl. Mech. Engrg.*, 82:281-300, 1990.

39. M. Bellet and A. Heinrich. A two-dimensional finite element thermomechanical approach to a global stress-strain analysis of steel continuous casting. *ISIJ Int.*, 44(10):1686-1695, 2004.

40. N. Ahmad, H. Combeau, J.-L. Desbiolles, T. Jalanti, G. Lesoult, M. Rappaz, and C. Stomp. Numerical simulation of macrosegregation: a comparison between finite volume method and finite element method predictions and a confrontation with experiments. *Metall. and Mat. Trans. A*, 29:617-630, 1997.

41. V. D. Fachinotti and M. Bellet. Three-dimensional finite element modelling of thermomechanics and macrosegregation in binary alloys solidification. In *Proc. ESAFORM 6th Int. Conf. On Material Forming*, pages 511-514, Salerno, Italy, April 2003. University of Salerno.

42. V. D. Fachinotti and M. Bellet. Modélisation tridimensionnelle par éléments finis de la macroségrégation au cours de la solidification d'un alliage binaire. In *Proc. 16ème Congrès Français de Mécanique*, Nice, France, September 2003. Université de Nice-Sophia Antipolis-Ecole des Mines de Paris.

43. S. Gouttebroze, V. D. Fachinotti, M. Bellet, and H. Combeau. 3D-FEM modeling of macrosegregation in solidification of binary alloys. *Sent to Int. J. Forming Proc.*, 2005.

44. O. C. Zienkiewicz and R. L. Taylor. *The Finite Element Method*, volume 1: The Basis. Butterworth-Heinemann, 5th. edition, 2000.

45. M. Bellet and V. D. Fachinotti. ALE method for solidification modeling. *Comput. Methods Appl. Mech. Engrg.*, 193:4355-4381, 2004.

46. M. Bellet and J.-L. Chenot. The ALE method for the numerical simulation of material forming processes. In S.-F. Shen and P. R. Dawson, editors, *Proc. NUMIFORM'95, 5th Int. Conf. on Numerical Methods in Industrial Forming Processes*, pages 39-48, Ithaka, NY, USA, June 1995.

47. H. Matthies and G. Strang. The solution of nonlinear finite element equations. *Int. J. Numer. Methods Engrg.*, 14(11):1613-1626, 1979.

48. S. Balay, K. Buschelman, V. Eijkhout, W. Gropp, D. Kaushik, M. Knepley, L. Curfman McInnes, B. Smith, and H. Zhang. *PETSc Users Manual. Revision 2.2.1*. Mathematics and Computer Science Division, Argonne National Laboratory, http://www.mcs.anl.gov/petsc, 2004.

49. A. N. Brooks and T. J. R. Hughes. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. *Comput. Methods Appl. Mech. Engrg.*, 32:199-259, 1982.

50. T. Tezduyar and S. Sathe. Stabilization parameters in SUPG and PSPG formulations. *J. Computat. Appl. Mech.*, 4(1):71-88, 2003.

51. G. Lesoult, C.-A. Gandin, and N. T. Niane. Segregation during solidification with spongy deformation of
the mushy zone. *Acta Mater.*, 51:5263–5283, 2003.

52. K. Miyazawa and K. Schwerdtfeger. Macrosegregation in continuously cast steel slabs: preliminary theoretical investigation on the effect of steady state bulging. *Arch. Eisenhüttenwes*, 52:415–422, 1981.

53. G. Lesoult and S. Sella. Analysis and prevention of centreline segregation during continuous casting of steel related to deformation of the solid phase. In *Proc. 6th Int. Iron and Steel Congress*, volume 1, pages 681–688, Nagoya, 1990.

54. T. Kajitani, J.-M. Drezet, and M. Rappaz. Numerical simulation of deformation induced segregation in continuous casting of steel. *Metall. and Mat. Trans. A*, 32:1479–1491, 2001.

55. P. F. Kozlowski, B. G. Thomas, J. A. Azzi, and H. Wang. Simple constitutive equations for steel at high temperature. *Metall. and Mat. Trans. A*, 23:903–918, 1992.

56. M. Reza Aboutalebi, M. Hasan, and R. I. L. Guthrie. Coupled turbulent flow, heat, and solute transport in continuous casting processes. *Metall. and Mat. Trans. B*, 26:731–744, 1995.