Volume-Averaged Modeling of Multiphase Flow Phenomena during Alloy Solidification

The most recent developments and applications in volume-averaged modeling of solidification processes have been reviewed. Since the last reviews of this topic by Beckermann and co-workers [Applied Mech. Rev. 1993, p. 1; Annual Rev. Heat Transfer 1995, p. 115], major progress in this area has included (i) the development of a mixed columnar-equiaxed solidification model; (ii) further consideration of moving crystals and crystal dendritic morphology; and (iii) the model applications to analyze the formation mechanisms of macrosegregation, as-cast structure, shrinkage cavity and porosity in different casting processes. The capacity of computer hardware is still a limiting factor. However, many calculation examples, as verified by the laboratory casting experiments, or even by the casting processes at a small industrial scale, show great application potential. Following the trend of developments in computer hardware (projection according to Moore’s law), a full 3D calculation of casting at the industry scale with the multiphase volume-averaged solidification models will become practically feasible in the foreseeable future.


Introduction
One scientific challenge in developing a solidification model which is applicable for engineering casting processes lies in the requirement to bridge the length scales, i.e., to incorporate microstructure-scale phenomena in large-scale processes relevant to engineering or vice versa [1][2][3][4]. Numerous numerical models have been proposed targeting the phenomena at a certain level of length scale [5][6][7][8], but no single model has been able to span the entire range of length scales from the process scale of industrial interest (~m) down to the atomic scale (~Å) [9]. Even with an optimistic projection of the advancement of the computational capacity, a microstructure-scale model that accounts for the evolution of the dendritic crystal morphology via phase field or level set methods would not be feasible for simulating a 2D casting (~m) until the year 2055 [10]. In the view of engineering applications, a model targeting large-scale phenomena is favorable. Another issue of solidification is the involved multiphase flow: bulk melt flow incorporated with moving equiaxed crystals, interdendritic flow in the mushy region, formation of open cavity or enclosed gas pores, etc. With the increasing computational capacity, a scientific discipline, called computational multi-fluid dynamics (CMFD) [11][12][13], since it appeared in the late last century, was introduced into the field of solidification by Beckermann and co-workers [3,4]. The frequently used CMFD approaches can be summarized in three categories: Eulerian-Lagrange; free-surface tracking; and, Euler-Euler.

Governing Equations and Volume Average Operation
The transport equations for mass, momentum, energy, species etc. (Table 1), valid at the microscopic scale, can be written in a general form where, ψ is the volume-weighted quantity. This quantity can be directly solved in the Eulerian framework with fixed grid. However, physical quantities (ξ) are usually defined as mass-weighted ones. They are related by the density ρ, i.e., ψ = ρξ. Note that in momentum conservation equation the pressure gradient force is treated as a part of source term.
In order to describe the multiphase transport phenomena during solidification, a representative volume is shown in Figure 1. For simplicity, here only a two-phase system is referred to. Equation (1) is only valid locally in each phase region. For the motion of the solid phase, we further assume the solid phase as a pseudo-fluid. In principle, the issue of solidification would not be solvable without tracking the solid-liquid interface information. As mentioned before, solving the details of dendritic morphology by tracking the moving liquid-solid interface with phase field or level set methods is computationally costly. Therefore, there should be greater focus on the global transport phenomena by treating the multiple phases as interpenetrating continua. A computational cell with the similar size of the representative volume of Figure 1 can be referred to. Tracking the morphological details of each crystal is not necessary, but their effect (e.g., as described as functions of the characteristic morphological lengths) on the global transport phenomena must be taken into account. Note that in the multiphase transport system, the term "phase", as primarily defined by the material scientists, is used differently in the field of CMFD and here in this paper. For example, in order to describe the hydrodynamic behavior of mixed columnar-equiaxed solidification, the primary crystals of equiaxed and columnar structures are treated as two different phases, i.e., equiaxed phase and columnar phase, because the former can move freely while the latter is considered not movable, although both of them belong to the same phase (crystal structure) in the material physics point of view.
By performing the volume averaging operation V 0 of Equation (1) over the representative volume V 0 of Figure 1, we obtain ∂ψ ∂t Note that the volume (superficial) average of a variable ψ k of the phase k over the representative volume V 0 is ψ k V 0 = 1 V 0 V k ψ k dV; and the intrinsic average of the variable ψ k of the phase k is: ψ k dV; they are related by ψ k V 0 = f k · ψ k V k . Therefore, a valid transport equation at the macroscopic scale can be derived (referring to phase k), with ψ = ρξ, dropping the intrinsic operator symbol V k for simplicity, it can be written in: The intrinsic quantity ξ k V k , written as ξ k hereafter, represents the mean ξ of phase k in the representative volume. Details about derivation of the macroscopic conservation equations were made by Beckermann et al. [3,40]. Meanings of each term in Equation (2) are summarized in Table 2.
where, ψ is the volume-weighted quantity. This quantity can be directly solved in the Eulerian framework with fixed grid. However, physical quantities ( ξ ) are usually defined as mass-weighted ones. They are related by the density ρ, i.e., ρξ ψ = . Note that in momentum conservation equation the pressure gradient force is treated as a part of source term.  −3 ] c [-] ( ) ' n [kg −3 ] 0 N In order to describe the multiphase transport phenomena during solidification, a representative volume is shown in Figure 1. For simplicity, here only a two-phase system is referred to. Equation (1) is only valid locally in each phase region. For the motion of the solid phase, we further assume the solid phase as a pseudo-fluid. In principle, the issue of solidification would not be solvable without tracking the solid-liquid interface information. As mentioned before, solving the details of dendritic morphology by tracking the moving liquid-solid interface with phase field or level set methods is computationally costly. Therefore, there should be greater focus on the global transport phenomena by treating the multiple phases as interpenetrating continua. A computational cell with the similar size of the representative volume of Figure 1 can be referred to. Tracking the morphological details of each crystal is not necessary, but their effect (e.g., as described as functions of the characteristic morphological lengths) on the global transport phenomena must be taken into account. Note that in the multiphase transport system, the term "phase", as primarily defined by the material scientists, is used differently in the field of CMFD and here in this paper. For example, in order to describe the hydrodynamic behavior of mixed columnar-equiaxed solidification, the primary crystals of equiaxed and columnar structures are treated as two different phases, i.e., equiaxed phase and columnar phase, Figure 1. A representative volume element, V 0 , contains two phases, k and j. They move at different velocities u k and u j . The volume of both phases are V k and V j , and their volume fraction are f k (= V k /V 0 ) and f j , and the sum of them makes 1.0. The total k-j interface area being enclosed in the volume element is A k , while the k-phase border area being enveloped by the volume element (red) is ∆A k .
Note: u int k is the velocity of the k-j interface, which is the sum of two parts: the moving velocity of the k-phase u k and the growth velocity of the k-j interface due to solidification/melting v int k . The j-k interface area concentration is calculated as S k (= A k /V 0 ). Table 2 are as follows.

Further explanations to
i.
To calculate the diffusive fluxes J k , effective properties µ eff k , λ eff k , D eff k are to be defined. They take usually the values of their physical counterparts. For solid phase k, µ eff k is modeled as a function of f k by considering the solid particle impingement and their packing behavior. ii.
The dispersive terms J d k were not considered in most multiphase solidification models, or treated as an additional part of the effective diffusive fluxes. When a highly turbulent flow is involved, they must be modeled explicitly. iii.
The superscript * indicates the interface quantities ξ * k , which are different from the volume averaged quantity ξ k . It is the ξ * k and its difference to ξ k that govern the solidification/melting induced exchange terms ψ M k ; it is the difference (ξ * k − ξ k ) that determines diffusive transfer terms ψ D k . If the k phase is solid, u * k = u k . iv.
The surface average quantity is noted with an over bar ξ * k , and ξ * Latent heat (∆h) is treated by the global definition of enthalpy, h k = h k, ref + T k T k, re f c p,k dT.
Liquid and solid at the same temperature have different enthalpies, and their difference makes the ∆h. vi. The pressure p is assumed to be shared by all phases. The pressure gradient force ∇p is treated as a part of source term in momentum conservation equations. vii.
Some symbols such as v int k will be described and discussed in late sections.
The striking feature of this volume average approach is that the volume-averaged quantities ξ k can be directly solved at the macroscopic scale, while the local solidification phenomena, which happen at the microscopic (interfacial) scale, can be considered in different terms J k , J d k , ψ M k , ψ D k , . ψ k . Table 2. Conservation equations of macroscopic transport system (multiphase) [3,40].

Solution Procedures/Strategy
Diverse methods solving Eulerian multiphase flow are available [11,12,[42][43][44][45]. Most recent multiphase problems can be solved with the commercial solvers, e.g., the CFD software, ANSYS FLUENT (ANSYS Inc., Canonsburg, PA, USA) [42]. As shown in Figure 2, the solver provides a platform to solve the global transport equations (Equation (2) and Table 2). In the meantime it provides flexibility (open program interface) in defining the exchange and source terms for the transport equations, and even allows modification of the solution procedure. The multiphase transport equations are highly coupled and non-linear; they must be solved in an iterative manner. For each time step, up to 60 iterations are sometimes necessary to decrease the normalized residuals of f k , u k , c k , p and n k to a value below the convergence limit of 10 −4 , and h k below that of 10 −7 . In each iteration some intermediate (auxiliary) quantities, e.g., morphological parameters describing the shape and average diameter of crystals, are updated, and the exchange terms and the source terms are calculated. The formulation of the solver is fully implicit, and no stability criterion has to be met. However, due to the complexity of the multiphase coupling, the time step ∆t should be limited to ensure a sufficient computational accuracy. The optimal time step can only be obtained empirically by trial simulations in such multiphase system, or by using dynamic time step control. To perform the volume average operation, the representative volume must be large enough to include the characteristic size of the interfacial structures (primary or secondary dendrite arm spacing) ( Figure 1) [3,40], but the computational cell (volume element) as used for solving the global transport equation system is independent from this constraint. In some cases, the finer the cell, the better the resolution that will be obtained. Note that the calculated results for the transport quantities, e.g., f k , u k , h k , c k represent the volume averaged quantities. It does not matter how fine a cell is (even finer than secondary arm spacing of the dendrites); those quantities cannot be interpreted as field quantities of the interdendritic region. Note that the calculated results for the transport quantities, e.g., k f , k u  , k h , k c represent the volume averaged quantities. It does not matter how fine a cell is (even finer than secondary arm spacing of the dendrites); those quantities cannot be interpreted as field quantities of the interdendritic region.

Models at a Glance
A series of volume-averaged solidification models have been derived, as summarized in Table  3 [46,47]. This paper is not going to describe details of all those models, but will select some modeling examples to explore their functionalities/capabilities and limitations.

Models at a Glance
A series of volume-averaged solidification models have been derived, as summarized in Table 3 [46,47]. This paper is not going to describe details of all those models, but will select some modeling examples to explore their functionalities/capabilities and limitations. [ [75][76][77][84][85][86][87][88] * This reference list cannot include all contributions. Some pioneer works, on which the recent solidification models were further developed, can be read in the previous reviews of this topic [3,4].

Two-Phase Model for Equiaxed Solidification (Non-Dendritic)
Many industrial alloys solidify with equiaxed structures. In order to promote homogenous mechanical properties, the alloys are often modified by adding a grain refiner to achieve fine near-globular crystal morphology instead of the distinct dendrites. As a numerical model, the grain morphology can be simplified as spherical. The growth kinetic is dominated by diffusion around the globular crystal. A two-phase solidification model for globular equiaxed solidification was first proposed by Ni and Beckermann [40,41], then expanded by Ludwig and Wu for refining some closure laws [55][56][57], recently by Wu et al. [64,65] and Bedel et al. [89] for considering transport of both inoculants (refiner) and globular crystals.

Model Description
The general description and model assumptions are as follows: -Two phases are the primary liquid melt (phase-) and solid grains (phase-s). Referring to Figure 1 and Table 2. -A continuous heterogeneous nucleation law is applied for the grain nucleation [55,56]. -A 1D (spherical coordinate) steady state analytical solution is used to approximate the diffusion-governed grain growth ( Figure 3, Equation (3a)). The solid phase is treated as a pseudo-fluid with an artificial viscosity, µ s [3,90], which increases with the solid volume fraction from the physical value of viscosity of the bulk melt to the 'infinite value' at a so-called packing limit, f c s . It takes 0.637 for ideally-dense packed spheres of equal diameter. To mimic the coalescence of dendritic crystals, a packing limit with much lower solid volume fraction (e.g., 0.2) can be used [91,92].
-Viscous interaction between settling grains and the melt is modeled with a drag law, and the flow through the packed grains is modeled as the flow through porous media with Darcy law. -A very large (artificial) volume heat exchange rate is applied between the two phases to mimic the thermal equilibrium condition. -A Boussinesq approximation is made for the thermo-solutal convection and grain sedimentation. c .
-The solid phase is treated as a pseudo-fluid with an artificial viscosity, s μ [3,90], which increases with the solid volume fraction from the physical value of viscosity of the bulk melt to the 'infinite value' at a so-called packing limit, c s f . It takes 0.637 for ideally-dense packed spheres of equal diameter. To mimic the coalescence of dendritic crystals, a packing limit with much lower solid volume fraction (e.g., 0.2) can be used [91,92]. -Viscous interaction between settling grains and the melt is modeled with a drag law, and the flow through the packed grains is modeled as the flow through porous media with Darcy law.
-A very large (artificial) volume heat exchange rate is applied between the two phases to mimic the thermal equilibrium condition.
-A Boussinesq approximation is made for the thermo-solutal convection and grain sedimentation. As shown in Table 2 [70,71]. Rs is the radius of the equiaxed grain, as calculated by . As shown in Figure 3, the average growth velocity int s v (the same as the local growth velocity int s v  ) can be treated as a Stefan problem according to energy, species balance, and thermodynamics at the solid-liquid interface [40,41].
(3) Figure 3. Schematic of the spherical crystal growth of a binary alloy and solute diffusion profile.
As shown in Table 2, the solidification rate (M s = v int s ρ s S s ) is the key parameter for determining other solidification-induced exchange terms ψ M s . The solid-liquid interface area concentration is S s = n s · 4πR 2 s · Φ imp , where Φ imp is a surface impingement factor. Φ imp can be simply modeled as f , or treated as an function of f [70,71]. R s is the radius of the equiaxed grain, as calculated by R s = 3 3 f s /4πn s . As shown in Figure 3, the average growth velocity v int s (the same as the local growth velocity v int s ) can be treated as a Stefan problem according to energy, species balance, and thermodynamics at the solid-liquid interface [40,41].
There are 4 equations for 4 unknowns v int s , T * , c * , c * s . The quantities c , c s , T , T s are obtained from solving the global transport equations. Actually, most recent models assume a thermal equilibrium Metals 2019, 9, 229 9 of 43 condition in the representative volume element (T ≈ T s ). The species diffusion is the only governing factor for the crystal growth. Therefore, Equation (3) is further simplified as: With the known v int s , namely M s , other solidification/remelting induced exchange terms U M s , Q M s , C M s can be calculated. The diffusion lengths are estimated as l j = R s , l j s = R s /2, assuming a steady state diffusion field during spherical crystal growth (Figure 3). If a parabolic diffusion profile inside the solid is assumed, l q s (= l j s ) can be estimated as R s /5 [40]; and if the inter-grain impingement of the diffusion field is considered, l q (= l j ) should be calculated as R s (1 − R s /R f ) [70,71]. This diffusion length must be modified if the effect of turbulent flow is to be considered [41]. Referring to Table 2, the average interfacial velocity u * s takes the liquid velocity u for the case of solidification, the solid velocity u s for the case of remelting; the average interfacial enthalpy h * s takes the liquid enthalpy h for the case of solidification, the solid enthalpy h s for the case of remelting (h − h s = ∆h, the latent heat); the average interfacial concentration c * s takes the thermodynamic equilibrium concentration c * s of the solid phase for the case of solidification, the volume-averaged solid concentration c s for the case of remelting. For most metal alloys, the second term of right-hand-side in Equation (3a), which is used to consider back diffusion in solid, is negligible. One known exceptional case is that the carbon back diffusion in steel is important. The law for heterogeneous nucleation has been verified for most technical alloys. The inoculants (embryos), with their initial number density (n max ) existing in the parent melt, will be activated with undercooling as equiaxed nuclei. The remaining inoculants and the as-activated equiaxed crystals are quantified by their number densities: n em and n s . The transport of n em is calculated according to the liquid velocity u , differently from n s , which is transported according to u [64,65].
The source term N nu determines the nucleation rate, which follows a Gaussian distribution function of undercooling.
n max , ∆T N and ∆T σ are nucleation parameters, to be determined experimentally [79,93]. The solid phase is treated as a pseudo-fluid. Therefore, the calculation of the corresponding momentum conservation equation requires the definition of a viscosity of the solid phase, µ s , which is caused by collisions between individual equiaxed grains. Ishii and Zuber [90] found for the viscosity of a solid/liquid mixture: Assuming that the mixing rule is valid here [3], µ mix = f µ + f s µ s , the viscosity of the solid phase is derived as It increases exponentially with f s until it reaches infinite value at the packing limit ( f c s ). The drag force (viscous stress transfer) between the liquid and solid phases is calculated by where K s is drag coefficient. Depending on the special configuration of the fluid/solid mixture, people have suggested different models for K s [94][95][96][97][98][99]. To calculate K s during globular equiaxed solidification, two scenarios are distinguished: the approach for the submerged-objects (Kozeny-Carman) for the case of small f s ; the approach for the flow-through-porous-medium (Blake-Kozeny) for the case of large f s .
where d s (= 2 × R s ) is the diameter of the equiaxed grain. K is the Darcy permeability, as calculated by K = K 0 f 3 / f 2 s . In order to obtain a smooth transition of K s at f c s , we choose the empirical parameter K 0 = d 2 s /180. The above Equation (7) is valid for laminar flow, for which the Reynolds number, R e = d s ρ f u − u s /(µ · f s ) [99], is less or equal to 10. For the turbulent flow, other formulations must be considered [41,94,100].

Application of Globular Equiaxed Solidification Model
The first simulation with this model was made by Ni and Beckermann [41] to explore the crystal sedimentation. Classical knowledge about the crystal sedimentation is superficial, as plotted by Ohno in Figure 4a [101]. No quantitative correlation between the sedimentation induced solidification results (e.g., grain size distribution, macrosegregation) and the process variables (e.g., temperature, constitutional undercooling) could be derived, because the two-phase flow and its interplay with the local crystal growth kinetics and the rearrangement of the grain distribution were missing. Therefore, such exploratory simulations as done by Ni and Beckermann have become crucially important to achieve just a mere basic understanding of this correlation.  This globular equiaxed solidification model was later extended by different authors to investigate the formation of globular equiaxed crystals during channel cooling of A356 aluminum alloy (rheo-casting) [102,103], decomposition and macrosegregation during monotectic solidification [104][105][106], the crystal sedimentation induced macrosegregation in steel ingots [92,107,108], and premature solidification phenomenon during mold filling of shape casting [109]. Evaluation efforts were also made by comparing the simulations with laboratory casting experiments [110,111]. To some extent, the model has explained the grain sedimentation and its influence on the grain size distribution, but the use of this model for prediction and control of the as-cast structure in real casting processes is still . Using a two-phase globular equiaxed solidification model to study the crystal sedimentation phenomenon [55,56]. Numerical simulation was performed in a 2D (0.18 × 0.16 m 2 ) domain for an Al-4.0wt.%Cu alloy. (a) Schematic of the crystal sedimentation phenomenon [101]. Simulation results (80 s) for (b) the velocity of liquid ( u ,max = 2.37 mm/s) and isolines of fraction solid, (c) the velocity of equiaxed grains ( u e,max = 2.38 mm/s) and isolines of fraction solid, (d) the mixture concentration c mix (dark for positive segregation, light for negative segregation) and its isolines, (e) isolines of grain size d s in µm. Nucleation parameters: n max = 10 14 m −3 , ∆T N = 10 K, ∆T σ = 4 K.
Simulation of a 2D domain (Al-4.0wt.%Cu alloy), as cooled from side and bottom, was performed [55,56]. The metallic die at a constant temperature of 290 K was filled instantaneously with a melt at 925 K. The heat transfer coefficient at the casting-die interface was 750 W/(m 2 ·K). The simulation results to some extent have "reproduced" the scenario depicted in Figure 4a. Only the results at 80 s are shown in Figure 4b-e; nucleation and solidification start along the die wall. Some of the grains sink and settle at the bottom near the corners, while some are brought by the melt flow to the center region. Both liquid and solid phases are coupled through the momentum exchange terms. The melt is drawn by the sinking grains, forming two vortices: one clockwise in the right half and the other anticlockwise in the left half of the domain. Although thermo-solutal convection contributes to the flow, the crystal sedimentation and sedimentation induced melt convection dominate. The equiaxed grains accumulate when the local f s exceeds the packing limit f c s . The grain accumulation is responsible for the formation of negative segregation (A and D zones). Two mechanisms are responsible for positive segregation: feeding the packed zones with segregated melt and squeezing it out of the segregated melt by settling grains. The positive segregation zones are not stationary during solidification, they move with the flow. As shown in Figure 4e, the grain size distribution is relatively uniform (about 55 µm) along the side wall. Relatively large grains appear at the bottom of the casting center. The grains as nucleated near the die wall sink and grow. As they reach the lower central regions, they have grown to a grain size of 88 µm.
This globular equiaxed solidification model was later extended by different authors to investigate the formation of globular equiaxed crystals during channel cooling of A356 aluminum alloy (rheo-casting) [102,103], decomposition and macrosegregation during monotectic solidification [104][105][106], the crystal sedimentation induced macrosegregation in steel ingots [92,107,108], and premature solidification phenomenon during mold filling of shape casting [109]. Evaluation efforts were also made by comparing the simulations with laboratory casting experiments [110,111]. To some extent, the model has explained the grain sedimentation and its influence on the grain size distribution, but the use of this model for prediction and control of the as-cast structure in real casting processes is still quite limited. Firstly, the globular morphology is overestimated; secondly, most cases of alloy solidification have a mixed columnar-equiaxed structure.

Mixed Columnar-Equiaxed Solidification (Non Dendritic)
The most general case of as-cast structure for industry castings is a mixed columnar-equiaxed microstructure. During solidification crystals grow in columnar and/or equiaxed manner. The equiaxed crystals can move, while columnar are stationary (or move with a predefined velocity). Due to this difference, both equiaxed and columnar crystals are necessarily treated as separated phases, although in thermodynamic point of view, they belong to the same phase.

Model Description
This model was proposed by the current authors [64,65]. As schematically shown in Figure 5, 3 phases are considered. The general description and model assumptions are as follows.
-Phase definition: primary liquid (phase-), equiaxed as the first secondary phase (phase-e), columnar as the second secondary phase (phase-c). -Both primary liquid and equiaxed phases are moving phases, for which the corresponding momentum equations are solved. The columnar phase is assumed to stick to the mold wall (stationary), and solidify from the wall towards the bulk melt. -Similar to the Section 3.1.1, the origin of the equiaxed crystals is modeled by a continuous heterogeneous nucleation law [55,56]. - The columnar grains are assumed to originate from the mold wall. The columnar tip position is explicitly tracked. Hunt's blocking mechanism [112] is applied for predicting the columnar-to-equiaxed transition (CET). -Ideal morphologies for both solid phases are assumed: spheres for equiaxed (globular) grains, and cylinders for columnar (cellular) dendrite trunks. The radial growth is controlled by diffusion according to the analytically solved concentration profiles around the corresponding crystal. A simplified approach is suggested to treat some hydrodynamic behaviors of equiaxed dendrites [113]. -Both the liquid melt and the solid phases have different (volume averaged) concentrations, c , c e (or c c ); the thermodynamic equilibrium condition applies at the liquid-solid interface with interfacial concentrations, c * , c * e (or c * c ). -A very large (artificial) volume heat exchange rate is applied among the three phases to mimic the thermal equilibrium condition. -Thermo-solutal convection and grain sedimentation are usually modeled with the Boussinesq approach. Feeding flow due to different densities between liquid and solid must be treated differently. For example, during continuous casting, an open calculation domain is considered, and hot melt is allowed to be conducted from an inlet into the domain [60]. In ingot/shape castings, an additional gas/slag phase must be added to feed the shrinkage cavity [114].

Model Description
This model was proposed by the current authors [64,65]. As schematically shown in Figure 5, 3 phases are considered. The general description and model assumptions are as follows.
-Phase definition: primary liquid (phase- ), equiaxed as the first secondary phase (phase-e), columnar as the second secondary phase (phase-c).
-Both primary liquid and equiaxed phases are moving phases, for which the corresponding momentum equations are solved. The columnar phase is assumed to stick to the mold wall (stationary), and solidify from the wall towards the bulk melt.
-Similar to the Section 3.1.1, the origin of the equiaxed crystals is modeled by a continuous heterogeneous nucleation law [55,56].
-The columnar grains are assumed to originate from the mold wall. The columnar tip position is explicitly tracked. Hunt's blocking mechanism [112] is applied for predicting the columnar-toequiaxed transition (CET).
-Ideal morphologies for both solid phases are assumed: spheres for equiaxed (globular) grains, and cylinders for columnar (cellular) dendrite trunks. The radial growth is controlled by diffusion according to the analytically solved concentration profiles around the corresponding crystal. A simplified approach is suggested to treat some hydrodynamic behaviors of equiaxed dendrites [113]. -A very large (artificial) volume heat exchange rate is applied among the three phases to mimic the thermal equilibrium condition. -Thermo-solutal convection and grain sedimentation are usually modeled with the Boussinesq approach. Feeding flow due to different densities between liquid and solid must be treated differently. For example, during continuous casting, an open calculation domain is considered, and hot melt is allowed to be conducted from an inlet into the domain [60]. In ingot/shape castings, The way to treat the equiaxed solidification is same as that described in Section 3.1.1. Solidification of the columnar and the treatment of its interaction with the equiaxed are described below.

Columnar Tip Front Tracking Algorithm
The columnar tips advance with v tip (tip growth velocity). The columnar front is assumed to grow in the preferred direction that is closest to the temperature gradient ( Figure 6). An algorithm to track the columnar tip front, which applies to both structured and unstructured meshes and to both 2D and 3D, is proposed.
-Each control volume is marked with a status marker, i c . When i c = 1, it indicates that the control volume contains the columnar tip front; when i c = 2, it means that the volume element has been passed by the columnar tip front; otherwise when i c = 0, it is still in the bulk region where no columnar tip has yet reached it. All elements are initialized with i c = 0, except for the boundary (wall) elements, which are marked as initial condition with i c = 1. -To each control volume, a reference length, l ref , is assigned. l ref is used to determine the length for the columnar front to grow out of a volume element. l ref is defined by seeing the control volume as an equivalence sphere in 3D (or circle in 2D) with a diameter of l ref . The volume of the sphere must be equal to that of the corresponding control volume: 4π The growth velocity of the primary dendrite columnar tip, v tip , and the tip radius, R tip , are determined by the LGK (Lipton-Glicksman-Kurz) model [72]. Note that in the late works of the authors, the KGT model [88] was used to track the columnar tip front [84,85]. The average length of the columnar grains in the volume element as marked as columnar tip (i c = 1) is calculated by the integral l = t v tip dt, starting from the moment when the front enters the control volume element. The neighboring control volume elements, which are still marked with i c = 0, will be converted into the columnar tip elements and their status is changed to i c = 1. In the meantime the status of the first volume element, just passed by the columnar tip front, is set to i c = 2. -A mass transfer between the liquid to the columnar phases (solidification/melting) occurs only for those elements of i c = 0.

Columnar Tip Front Tracking Algorithm
The columnar tips advance with vtip (tip growth velocity). The columnar front is assumed to grow in the preferred direction that is closest to the temperature gradient ( Figure 6). An algorithm to track the columnar tip front, which applies to both structured and unstructured meshes and to both 2D and 3D, is proposed.
-Each control volume is marked with a status marker, c i . When -The growth velocity of the primary dendrite columnar tip, tip v , and the tip radius, tip R , are determined by the LGK (Lipton-Glicksman-Kurz) model [72]. Note that in the late works of the authors, the KGT model [88] was used to track the columnar tip front [84,85]. The average length of the columnar grains in the volume element as marked as columnar tip ( , starting from the moment when the front enters the control volume element.
-As soon as l exceeds ref l , the columnar tip front grows out of the considered volume element.
The neighboring control volume elements, which are still marked with 0 c = i , will be converted into the columnar tip elements and their status is changed to In the meantime the status Figure 6. Schematic of the advancing columnar tip front, and the algorithm of the columnar front tracking during mixed columnar-equiaxed solidification [65,84]. Reproduced from [84], with permission from Elsevier, 2019.

Competition Between Columnar and Equiaxed
Crystals start to grow in the chill zone (casting surface) as an equiaxed morphology. As the solidification progresses, a temperature gradient is gradually established. The growth of the crystals becomes unidirectionally against the heat flux direction, and then columnar structure develops. This is called as equiaxed-to-columnar transition, i.e., ECT. This chill zone is here ignored, and it is assumed that the columnar grains originate directly from the mold wall. In the late stage of solidification, new equiaxed crystals nucleate and grow in the liquid ahead of columnar zone. The growth of a columnar structure will be blocked by the equiaxed crystals leading to columnar-to-equiaxed transition, i.e., CET. When the motion of equiaxed crystals and their interactions with the growing columnar phase are considered, the aforementioned ECT and CET become more complex. The following algorithm is suggested.
-Both columnar and equiaxed phases can coexist in the same volume element. When f c (local volume fraction of the columnar phase) is more than a critical value, f free c (e.g., 0.2), an infinite drag force coefficient is applied between the two solid phases in the corresponding momentum equations, and thus the equiaxed grains are "captured" by the columnar phase. When f c < f free c , no drag force is applied, and thus, the motion of the equiaxed grains is not affected by the columnar dendrite tips. The validation of the critical value f free c = 0.2 is discussed elsewhere [65]. - The columnar tip blocking mechanism of Hunt [112] is implemented for the CET. In those volume elements, that contain the columnar primary dendrite tips, the tip growth velocity, v tip , is set to zero as soon as the local volume fraction of the equiaxed grains, f e , exceeds the critical threshold of f e,CET = 0.49. A recent study suggested a value of 0.2 should be taken for f e,CET [115,116]. The impingement of the solute concentration fields between equiaxed and columnar growths, also known as the 'soft blocking mechanism' for CET as proposed by Martorano et al. [117], is automatically included in the model. - Equiaxed crystals coagulate with each other. When f e (volume fraction of equiaxed) becomes larger than the packing limit (0.637), the equiaxed phase builds a rigid network. As the equiaxed phase hits the mold wall (stationary), the velocity of the equiaxed phase is set zero, and then the equiaxed phase is considered as rigid and stationary as well. This rigid and stationary equiaxed phase region is tracked (numerically marked). New columnar dendrite tip front will initiate/grow from the outer contour of the marked rigid-equiaxed crystal region, triggering the ECT again [63]. - The columnar solidification competes with the equiaxed solidification. When the growth of the columnar tip front "overtakes" the growth of equiaxed phase, ECT is triggered. Otherwise, the growth of the columnar is suppressed and CET is triggered. The following algorithm is implemented to judge if the growth of columnar can "overtake" the growth of equiaxed phase. If the columnar tip front can grow out of the considered volume element before f e reaches a blocking limit ( f e,CET ), the columnar primary tip front can continue to grow into the neighboring cell or cells. In contrast, if the columnar tip front cannot grow out of the considered volume element before f e reaches f e,CET , CET occurs.

Diffusion Governed Growth of Columnar Phase
Referring to Table 2, the solidification rate of a solid phase is calculated by M s = v int s ρ s S s . Here different subscript, M c , is used to describe the interphase mass transfer (instead of the symbol M s , as here more than one solid phases are considered). The subscript s will be replaced by corresponding c or e. To calculate M c , 2 states of the volume elements have to be distinguished: (i) volume elements that have been passed by the columnar tips (i c = 2) and (ii) volume elements which contain the columnar tips (i c = 1).
For the elements of i c = 2, where the columnar trunks have passed the volume element: while for the elements of i c = 1, which include columnar primary dendrite tips: The radial growth velocity of the columnar trunk is calculated as where l j ,c and l j s,c are species diffusion lengths in the interdendritic liquid and in the solid columnar dendrites, and they are calculated as l j ,c = R c · ln(R f /R c ) and l j s,c = R c /2 (with an assumption of parabolic diffusion profile inside the solid dendrite [40], l j s,c can be taken as R c /5). Depending on different arrangements (aligned and staggered) of the columnar trunks, the radius of the columnar trunk (R c ) can be calculated: array. R f is the maximum radius that the columnar trunk can reach at f c = 1. The growth area S trunk c is 2πR c /λ 2 1 · Φ c,imp for the aligned array, 2/ √ 3 · 2πR c /λ 2 1 · Φ c,imp for the staggered array. The growth area impingement factor Φ c,imp is min[ f c / (1− π/4), 1] for the aligned array, min f c / (1− π/ √ 12 , 1 for the staggered array. The growth velocity of the primary dendrite tip (v tip ) and the tip radius (R tip ) are known from the LGK model [72]. The growth area (projection in the growth direction) of the columnar dendrite tip is calculated as S c, tip = n c · πR 2 tip · f , where n c is the number density of the columnar trunks ( f c / πR 2 c · l ), and f is here considered as the growth area impingement factor in the volume element considering columnar dendrite tips. The surface area of columnar trunks in the volume element containing columnar primary dendrite tips is calculated as S trunk c,tip = n c · 2πR c l · f , where l is the length of the columnar trunk in the considered volume element.

Permeability in the Columnar Mushy Zone
The drag force between interdendritic liquid and the columnar dendrites is calculated by The Blake-Kozeny approach is used to calculate the drag coefficient K c = − f 2 µ /K, where the isotropic permeability K is described as a function of the fraction liquid and the primary dendrite arm spacing [118]: Literature data on K in solidifying dendrites show a wide scatter [119][120][121][122][123]. Permeability is a volume-averaged quantity, which was originally defined to describe the flow through porous medium of "homogenous structure". The solidifying crystal structure (morphology) is not homogenous. In addition to the fraction liquid (or solid) and primary arm spacing, other morphological parameters, such as structural anisotropy, tortuosity, etc., influence the interdendritic flow direction as well. This paper is not intended to represent comprehensive review on this topic, but the importance of choosing an appropriate permeability law should not be underestimated.

Application of Mixed Columnar-Equiaxed Solidification Model
The first exploratory simulation using this model to explain mixed columnar-equiaxed solidification was made [64,65]. As shown in Figure 7, a model ingot with an under-scaled size was considered. The overall solidification sequence is governed by heat transfer, but it is strongly influenced by the melt convection and the grain sedimentation. The solidification, i.e., the growth of columnar phase, starts when the surface temperature drops below liquidus. In the meantime, the equiaxed grains nucleate and grow. Near the wall, they sink and induce melt convection. The sinking crystals drag the melt downwards along the wall, and induce a rising melt flow in the middle of the casting. Additionally, the thermal-solutal buoyancy also contributes to the overall melt convection. With the progress of solidification, the volume fraction of the columnar phase (f c ) in the mold wall region increases. In the vicinity of the columnar tip front the equiaxed grains continue to nucleate, grow and sink. The sedimentation of the equiaxed grains plays an important role in the final as-cast structure. The sinking grains pile up in the lower region of the ingot. This pile-up of the equiaxed grains can block the growth of columnar tip front in the lower region, causing CET in an enclosed conic shape. It is this kind of convection and grain sedimentation phenomena that are responsible for the formation of macrosegregation as typically observed in steel ingots, i.e., a positive segregation at the top region and a negative segregation cone, corresponding to the equiaxed zone, in the bottom. This model has been recently extended to simulate the formation of the as-cast structure and macrosegregation in industry ingots [113,114,124]. Major extensions to the original model are to treat the dendritic morphology of the equiaxed crystals and the formation of shrinkage cavity. The previous assumption of the ideal globular morphology of the equiaxed grain leads to overestimation of the macrosegregation [64,65,69,93].
As inspired by the "grain envelope approach" [73][74][75][76][77][78][79][80][81][82], a simplified dendrite envelope approach [113,124] is suggested to treat (i) the liquid-equiaxed (dendritic) drag force, (ii) the blocking mechanism of the columnar tip growth by the dendritic equiaxed grains, and (iii) the 'solid viscosity' and the packing limit. Ignoring the solidification shrinkage by the previous model would cause error estimation of the as-cast structure and macrosegregation in the top region, and would also neglect the segregation, as caused by the shrinkage-induced feeding flow in the deep mush zone and near the end of the solidification region. To solve this problem, an additional phase (the phase number 4), i.e., gas phase (or covering liquid slag phase), is introduced into the model. The gas phase has no mass and species exchange with the metal phases, but it is coupled with them by taking into account the momentum and energy exchanges. This four phase mixed columnar-equiaxed solidification model was described elsewhere [114].
drag the melt downwards along the wall, and induce a rising melt flow in the middle of the casting. Additionally, the thermal-solutal buoyancy also contributes to the overall melt convection. With the progress of solidification, the volume fraction of the columnar phase (fc) in the mold wall region increases. In the vicinity of the columnar tip front the equiaxed grains continue to nucleate, grow and sink. The sedimentation of the equiaxed grains plays an important role in the final as-cast structure. The sinking grains pile up in the lower region of the ingot. This pile-up of the equiaxed grains can block the growth of columnar tip front in the lower region, causing CET in an enclosed conic shape. It is this kind of convection and grain sedimentation phenomena that are responsible for the formation of macrosegregation as typically observed in steel ingots, i.e., a positive segregation at the top region and a negative segregation cone, corresponding to the equiaxed zone, in the bottom. This model has been recently extended to simulate the formation of the as-cast structure and macrosegregation in industry ingots [113,114,124]. Major extensions to the original model are to treat the dendritic morphology of the equiaxed crystals and the formation of shrinkage cavity. The previous assumption of the ideal globular morphology of the equiaxed grain leads to overestimation of the macrosegregation [64,65,69,93].
As inspired by the "grain envelope approach" [73][74][75][76][77][78][79][80][81][82], a simplified dendrite envelope approach [113,124] is suggested to treat (i) the liquid-equiaxed (dendritic) drag force, (ii) the blocking mechanism of the columnar tip growth by the dendritic equiaxed grains, and (iii) the 'solid viscosity' and the packing limit. Ignoring the solidification shrinkage by the previous model would cause error estimation of the as-cast structure and macrosegregation in the top region, and would also neglect the segregation, as caused by the shrinkage-induced feeding flow in the deep mush zone and near the end of the solidification region. To solve this problem, an additional phase (the phase number 4), i.e., gas phase (or covering liquid slag phase), is introduced into the model. The gas phase has no mass and species exchange with the metal phases, but it is coupled with them by taking into account the   [125]. Experimental analysis of this ingot was recently published by Duan et al. [126]. Some most simulation relevant process parameters including casting geometry, alloy composition (Fe-0.51 wt.% C-0.23 wt.%Si-0.6 wt.%Mn-0.13 wt.%Ni-0.006 wt.%S-0.009 wt.%P), pouring temperature (1833 K), mold configuration, thermal boundary condition as derived from temperature measurements were reported; and macrosegregation map across a half vertical section as analyzed using infrared carbon-sulphur analyzer with good resolution (30 × 30 mm 2 ) were presented.
The solidification sequence is shown in Figure 8. The cooling and solidification starts from the mold wall. A columnar structure develops there and grows towards casting center. Equiaxed crystals nucleate and grow in front of columnar tips, and those equiaxed grains sink and settle at bottom of the ingot. The melt is dragged downwards by the sinking grains, which in turn induces a rising melt flow in the ingot center region. Both u and u e fields are instable and non-axisymmetric. As the melt flow and the motion of the equiaxed crystals are fully coupled with energy, species and mass transport phenomena, both the melt flow and crystal motion will influence the solidification sequence. Sedimentation of crystals in the bottom region will cause f e to increase significantly. As f e in front of the columnar tip front is high enough, CET (columnar-to-equiaxed transition) occurs. In the upper part of the ingot, the columnar tips can continue to grow. The flow (bulk and interdendritic) and crystal sedimentation are key mechanisms for the formation of macrosegregation in such ingots.
The simulation-experiment agreement is very promising, either regarding the segregation map (Figure 9) or the segregation profiles along different lines ( Figure 10). This approach also provides reasonable shrinkage cavity information. The macrosegregation map in the ingot is by no means axisymmetric. It is worth mentioning that the experimentally-determined segregation map was mirrored from half section of measured result. The modeling result, Figure 9a Section C, shows a strong tendency of A-segregates in the upper part of the ingot. momentum and energy exchanges. This four phase mixed columnar-equiaxed solidification model was described elsewhere [114].   [125]. Experimental analysis of this ingot was recently published by Duan et al. [126]. Some most simulation relevant process parameters including casting geometry, alloy composition (Fe-0.51 wt.% C-0.23 wt.%Si-0.6 wt.%Mn-0.13 wt.%Ni-0.006 wt.%S-0.009 wt.%P), pouring temperature (1833 K), mold configuration, thermal boundary condition as derived from temperature measurements were reported; and macrosegregation map across a half vertical section as analyzed using infrared carbon-sulphur analyzer with good resolution (30 × 30 mm 2 ) were presented.
The solidification sequence is shown in Figure 8. The cooling and solidification starts from the mold wall. A columnar structure develops there and grows towards casting center. Equiaxed crystals nucleate and grow in front of columnar tips, and those equiaxed grains sink and settle at bottom of the ingot. The melt is dragged downwards by the sinking grains, which in turn induces a rising melt melt flow and the motion of the equiaxed crystals are fully coupled with energy, species and mass transport phenomena, both the melt flow and crystal motion will influence the solidification sequence. Sedimentation of crystals in the bottom region will cause fe to increase significantly. As fe in front of the columnar tip front is high enough, CET (columnar-to-equiaxed transition) occurs. In the upper part of the ingot, the columnar tips can continue to grow. The flow (bulk and interdendritic) and crystal sedimentation are key mechanisms for the formation of macrosegregation in such ingots. The simulation-experiment agreement is very promising, either regarding the segregation map ( Figure 9) or the segregation profiles along different lines ( Figure 10). This approach also provides reasonable shrinkage cavity information. The macrosegregation map in the ingot is by no means axisymmetric. It is worth mentioning that the experimentally-determined segregation map was mirrored from half section of measured result. The modeling result, Figure 9a Section C, shows a strong tendency of A-segregates in the upper part of the ingot.
This mixed columnar-equiaxed solidification model was mostly applied to study macrosegregation in steel ingots [69,124,125,[127][128][129], but it was also used analyze solidification during direct chill casting of copper alloys [130,131] and continuous casting of steel [132]. In the meantime, this model is being continuously improved by considering the multicomponent alloy [69,70] and the fragmentation rates [133]. Laboratory evaluations of the model were made in addition to the aforementioned procedures [128,[134][135][136].

Incorporation of Dendritic Morphology
When incorporating a dendritic morphology into the multiphase volume-averaged solidification models, it is necessary to consider different phase regions. For example, the liquid melt, which is enclosed in the deep dendrite (interdendritic melt), has a far different solute concentration than the This mixed columnar-equiaxed solidification model was mostly applied to study macrosegregation in steel ingots [69,124,125,[127][128][129], but it was also used analyze solidification during direct chill casting of copper alloys [130,131] and continuous casting of steel [132]. In the meantime, this model is being continuously improved by considering the multicomponent alloy [69,70] and the fragmentation rates [133]. Laboratory evaluations of the model were made in addition to the aforementioned procedures [128,[134][135][136].

Incorporation of Dendritic Morphology
When incorporating a dendritic morphology into the multiphase volume-averaged solidification models, it is necessary to consider different phase regions. For example, the liquid melt, which is enclosed in the deep dendrite (interdendritic melt), has a far different solute concentration than the extra-dendritic melt. Although the inter-and extra-dendritic melts are in the same state of liquid, they have different thermodynamic properties and govern different crystal growth kinetics. In addition to its variety of concentrations, the diffusion length associated with the interdendritic melt is different from that with the extra-dendritic melt. Their hydrodynamic behavior may even differ: the interdendritic melt is more prone to be entrapped in the dendrite network moving together with or sticking to the crystal dendrite, while the extra-dendritic melt may be free to flow in other direction. It is necessary to consider the inter-and extra-dendritic melts as separate "phases".
As shown in Figures 11 and 12, a smooth hypothetical contour connecting primary (or primary, secondary and tertiary) dendrite tips (dashed line), called as grain envelope, is constructed to separate the inter-and extra-dendritic melts. Definition of the grain envelope seems subjective, but it influences the modeling accuracy. During metal alloy solidification, the Lewis number (ratio between thermal and mass diffusivities) is in the order of magnitude 3. A consideration of complete 'thermal mixing' at the scale of one grain (or even inside the representative volume of Figure 1) is widely accepted; hence, the growth of dendrite is dominantly controlled by solute diffusion, either inside or at the boundary of the grain envelope. A key assumption, necessarily made here, is that the interdendritic melt has a homogenous solute distribution which is significantly different from that of the extra-dendritic melt. The evolution of the dendritic grain (i.e., the expansion of grain envelope) is governed by the growth kinetics of the dendrite tips, which is dependent on the constitutional undercooling of the extra-dendritic melt and the diffusion length scale at the dendrite tip; while the solidification of the interdendritic melt is governed by the growth kinetics of the solid-liquid interface, which occurs at the diffusion length scale of the inter-secondary dendrite arm spacing. If a reasonable grain envelope could be defined (the volume averaged solute concentration c d can represent the local solute concentration of the interdendritic melt), and the diffusion or mixing length scale of the interdendritic melt is much smaller than the extra-dendritic melt, an assumption of complete 'solute-mixing' of the interdendritic melt (c d = c * = c env ) would be valid [73,74,[76][77][78][79][80]. In that case, solidification of the interdendritic melt (i.e., evolution of the solid phase inside envelope) could be derived indirectly from the mass/solute conservations. Otherwise the solidification of the interdendritic melt must be calculated explicitly according to the diffusion-governed growth kinetics at the solid-liquid interface [81,82,84,85]. interdendritic melt is much smaller than the extra-dendritic melt, an assumption of complete 'solutemixing' of the interdendritic melt ( d c = *  c = env c ) would be valid [73,74,[76][77][78][79][80]. In that case, solidification of the interdendritic melt (i.e., evolution of the solid phase inside envelope) could be derived indirectly from the mass/solute conservations. Otherwise the solidification of the interdendritic melt must be calculated explicitly according to the diffusion-governed growth kinetics at the solid-liquid interface [81,82,84,85]. One further modeling step is that the grain envelope is necessarily simplified as a volume equivalent sphere with a radius of R e . From the known growth velocity of the dendrite tips ( Φ is called as Figure 11. Concept of an equiaxed grain envelope and schematic of solute distribution in different phase regions. An equiaxed dendrite is presented by a "natural grain envelope" separating melts in the interand extra-dendritic regions, and the envelope contour connects the primary (or primary and secondary) dendrite tips (dashed line). This natural grain envelope is further modeled as a volume-equivalent sphere. The growth velocity of the volume equivalent sphere (v e env,M ) can be derived from the dendrite tip velocity (v e tip ) if the shape of the envelope is preserved/known. The grain envelope encloses two different phase regions: solid dendrites and interdendritic melt. Solute distribution profiles in different phase regions are shown with dotted lines, whose intrinsic averages (c s , c d , c ) are shown with solid lines. c * s , c * are thermodynamic equilibrium concentrations, valid at the solid-liquid interfaces; c env is the averaged concentration of the grain envelope. Reproduced from [84], with permission from Elsevier, 2019. Φ that helps estimation of the area of the "natural" grain envelope from the area of the volume equivalent sphere. If the envelope is shape-preserving during crystal growth, both parameters are constant, and they can be derived from the geometry of the envelope shape.    The shape of the columnar trunk envelope is simplified as a step-wise cylinder with cross sectional area equivalent to the tree-trunk envelope (dashed line), which connects the secondary and tertiary dendrite tips. The contour of the columnar envelope near the primary dendrite tip (dashed line connecting the primary and secondary dendrite tips) is simplified as a paraboloid. Reproduced from [84], with permission from Elsevier, 2019.
One further modeling step is that the grain envelope is necessarily simplified as a volume equivalent sphere with a radius of R e . From the known growth velocity of the dendrite tips (v e tip ), the radial growth velocity of the volume equivalent sphere (v e env,M = Φ e env · v e tip ) can be derived by conserving the volume of the 'natural' grain envelope with the volume equivalent sphere (4πR e3 /3), where Φ e env is the shape factor. As shown in Figure 11, the solute transfer by diffusion from the grain into the extra-dendritic melt is carried out at the surface of the 'natural' grain envelope (S e env,D : area concentration), which is identical to the interface between d and phases. Different from the surface area of the volume equivalent sphere (S e env,M = 3R e2 /R 3 f ), S e env,D = S e env,M /Φ e sph , where Φ e sph is called as sphericity of the grain envelope. It is crucial to have two morphological parameters, i.e., the shape factor Φ e env and the sphericity Φ e sph , to describe the grain envelope. It is Φ e env that connects the radial growth velocity of the volume equivalent sphere with the dendrite tip growth velocity; it is Φ e sph that helps estimation of the area of the "natural" grain envelope from the area of the volume equivalent sphere. If the envelope is shape-preserving during crystal growth, both parameters are constant, and they can be derived from the geometry of the envelope shape.
Examples of grain envelopes for both equiaxed and columnar dendrite structures are collected in Table 4. Note that two superscripts, e and c, are used to distinguish the equiaxed and columnar dendrites. For columnar dendrite (Figure 12), an area equivalent circle is used to model the cross section of columnar grain envelope. Again, two morphological parameters, i.e., shape factor Φ c env and the circularity Φ c circ , are used to describe the columnar grain envelope. Table 4. Morphological parameters for envelopes of selected dendrite structures [84]. Reproduced from [84], with permission from Elsevier, 2019.

Equiaxed † Columnar (Cross Section)
Sphere Φ e env = 1; Φ e sph = 1 Figure 12. The shape of the columnar trunk envelope is simplified as a step-wise cylinder with cross sectional area equivalent to the tree-trunk envelope (dashed line), which connects the secondary and tertiary dendrite tips. The contour of the columnar envelope near the primary dendrite tip (dashed line connecting the primary and secondary dendrite tips) is simplified as a paraboloid. Reproduced from [84], with permission from Elsevier, 2019.   Figure 12. The shape of the columnar trunk envelope is simplified as a step-wise cylinder with cross sectional area equivalent to the tree-trunk envelope (dashed line), which connects the secondary and tertiary dendrite tips. The contour of the columnar envelope near the primary dendrite tip (dashed line connecting the primary and secondary dendrite tips) is simplified as a paraboloid. Reproduced from [84], with permission from Elsevier, 2019.  Octahedral Figure 12. The shape of the columnar trunk envelope is simplified as a step-wise cylinder with cross sectional area equivalent to the tree-trunk envelope (dashed line), which connects the secondary and tertiary dendrite tips. The contour of the columnar envelope near the primary dendrite tip (dashed line connecting the primary and secondary dendrite tips) is simplified as a paraboloid. Reproduced from [84], with permission from Elsevier, 2019.  10 c circ π = Φ † equiaxed images taken from reference [6]. † † 6 orthogonal square pyramids (OSP6) [142,143] with pyramid angle 18.43°. † † † 4 orthogonal square wedges (OSW4) with wedge angle 60°. Figure 12. The shape of the columnar trunk envelope is simplified as a step-wise cylinder with cross sectional area equivalent to the tree-trunk envelope (dashed line), which connects the secondary and tertiary dendrite tips. The contour of the columnar envelope near the primary dendrite tip (dashed line connecting the primary and secondary dendrite tips) is simplified as a paraboloid. Reproduced from [84], with permission from Elsevier, 2019.  10 c circ π = Φ † equiaxed images taken from reference [6]. † † 6 orthogonal square pyramids (OSP6) [142,143] with pyramid angle 18.43°. † † † 4 orthogonal square wedges (OSW4) with wedge angle 60°.

Square rod
√ 3 Figure 12. The shape of the columnar trunk envelope is simplified as a step-wise cylinder with cross sectional area equivalent to the tree-trunk envelope (dashed line), which connects the secondary and tertiary dendrite tips. The contour of the columnar envelope near the primary dendrite tip (dashed line connecting the primary and secondary dendrite tips) is simplified as a paraboloid. Reproduced from [84], with permission from Elsevier, 2019.  10 c circ π = Φ † equiaxed images taken from reference [6]. † † 6 orthogonal square pyramids (OSP6) [142,143] with pyramid angle 18.43°. † † † 4 orthogonal square wedges (OSW4) with wedge angle 60°. Figure 12. The shape of the columnar trunk envelope is simplified as a step-wise cylinder with cross sectional area equivalent to the tree-trunk envelope (dashed line), which connects the secondary and tertiary dendrite tips. The contour of the columnar envelope near the primary dendrite tip (dashed line connecting the primary and secondary dendrite tips) is simplified as a paraboloid. Reproduced from [84], with permission from Elsevier, 2019.  10 c circ π = Φ † equiaxed images taken from reference [6]. † † 6 orthogonal square pyramids (OSP6) [142,143] with pyramid angle 18.43°. † † † 4 orthogonal square wedges (OSW4) with wedge angle 60°. † equiaxed images taken from reference [6]; † † 6 orthogonal square pyramids (OSP6) [142,143] with pyramid angle 18.43 • ; † † † 4 orthogonal square wedges (OSW4) with wedge angle 60 • .

Phase Definition
To describe the multiphase transport phenomena of the mixed columnar-equiaxed solidification ( Figure 5), three hydrodynamic phases are defined: e-, c-and -phases. They are quantified with their volume fractions f e , f c , f . They move with corresponding velocity: u e , u c , and u . The velocity of the c-phase, u c , is predefined: zero for the shape/ingot casting; the withdrawing speed for continuous casting; or as input field value from a thermal mechanical model (deforming dendrites). In order to consider the dendritic morphology, two distinct phase regions exist in each equiaxed or columnar structure: the solid dendrite and interdendritic melt. We assume that both solid dendrite and interdendritic melt inside one envelope share the same velocity. For the same mixed columnar-equiaxed solidification system, five phase regions are defined: the interdendritic melt in the equiaxed grain, the solid dendrites in the equiaxed grain, the interdendritic melt in the columnar dendrite trunk, the solid dendrites in the columnar dendrite trunk, and the extra-dendritic melt. We name the five phase regions as "thermodynamic" phases, and

• Equiaxed Growth
An equiaxed grain, immediately after nucleation, starts to grow with globular morphology. The growth of a globular grain (v e glob ) was described in detail in Section 3.1 by Equation (3a) [64,65]. A different symbol v int s was used for v e glob there. Following the growth-mode transition from globular to dendritic [81], the grain becomes dendritic. Then, the growth velocity of the primary dendrite tips is determined according to the LGK model [72,144].
here v e env, M is the growth velocity of the volume equivalent sphere, Ω is constitutional supersaturation c * − c / c * − c * s , Iv −1 (Ω) is the inverse Ivantsov function, which can be approximated as a Ω 1−Ω b with a = 0.4567 and b = 1.195 [74]. A simple approach is used to determine the globular-to-dendritic transition (GDT) by making a direct comparison between two growth velocities, v e glob and v e env, M . Therefore, the general formulation for the growth velocity of the equiaxed grain is v e env = max v e glob , v e env, M .
In another study, the stability criterion for growing spheres as proposed by Mullins and Sekerka was used to determine the GDT [145][146][147]. With the known v e env and S e env, M , the volume averaged mass transfer rate from -phase to e-phase can be calculated by M e = v e env · S e env, M · ρ e . • Columnar Growth An algorithm for tracking the columnar primary tip front was described in Section 3.2.1.1 ( Figure 6). Columnar trunks start to grow with cellular morphology. This morphology is simplified as a stepwise cylinder. The diffusion-governed growth of a cylinder is detailed in Section 3.2.1.3 [64,65]. Following the cellular-to-dendritic transition (CDT), the morphology of the trunk becomes dendritic. To model the dendritic growth of columnar structure, two zones are distinguished: the zone containing only columnar trunks, and the zone containing columnar primary dendrite tips (Figure 12). Near the primary dendrite tip, the envelope borders the tips of primary and secondary dendrite. In the zone far away from the primary dendrite tip, the envelope borders the tips of secondary and tertiary dendrites. The longitudinal section of the envelope near the primary dendrite tip is assumed to be a paraboloid. The shape of the dendrite trunk away from the primary dendrite tip can be described as being either a cylinder, square rod or four orthogonal square wedges connecting the secondary and tertiary dendrite tips ( Table 4).
The radial growth velocity of the volume equivalent cylinder (v c env, M ) is calculated from the growth velocity of the secondary dendrite tip: v c env, M = Φ c env · v c tip , where Φ c env is a shape factor. The diffusion area of the columnar envelope (S c env,D ) is estimated according to the surface area of the volume-equivalent cylinder (S c env,M ) by considering a circularity factor Φ c circ : S c env,D = S c env,M /Φ c circ . For cellular trunk, both Φ c env and Φ c circ are equal to one. Assuming a constant primary dendrite arm spacing, λ 1 , the volume-equivalent cylinder has an average diameter of d c (= 2 R c ). They are correlated to the fraction of columnar phase by the surface concentration S c env, M of the volume equivalent cylinder is calculated with where Φ c Imp , approximated as f , is an impingement factor. Equations (13) and (14) applies for both cellular and dendritic growth. Equation (14) is derived by assuming that the array of the columnar dendrite trunks is aligned. For the staggered array of the columnar trunks, the πd c /λ 2 1 is replaced by 2/ √ 3 · πd c /λ 2 1 . For the volume elements containing primary dendrite tips, a fictitious grain envelope (dashed line) enclosing the primary and secondary dendrite tips is shown in Figure 12. This envelope is simplified as a volume equivalent paraboloid, described elsewhere [84]. In a volume element that contains primary dendrite tips, a paraboloid with a diameter of d c and a length of l is used to represent the contour of the primary dendrite tip. The length of the paraboloid is explicitly tracked with the method described in Section 3.2.1.1. The same shape factor Φ c env and circularity Φ c circ are employed to simplify the shape of the parabolic envelope. d c and S c env, M of the volume-equivalent paraboloid in the volume element containing primary dendrite tips are calculated as Note that these formulations differ from Equations (13) and (14) for columnar trunk growth.

• Solidification of Melt in the Interdendritic Region
Solidification of the interdendritic melt-it does not matter if it is in both equiaxed and columnar growth-is driven by supersaturation and governed by solute diffusion in the interdendritic melt. The solidification rate of interdendritic melt is determined by the d-s interface growth velocity v e sd and the d-s interface area concentration S e s . The driving force for v e sd is c * − c e d / c * − c * s , but it is governed by diffusion at the length scale l e d , which is related to the secondary arm spacing λ 2 [75,81] by l e d = β 2 · (λ 2 − d 2 )/2, where β 2 is a constant on the order of unity and d 2 is the diameter of the secondary dendrite arms. It is assumed that d 2 is correlated to λ 2 by λ 2 − d 2 = λ 2 · α e d , thus: the d-s interface area inside the grain envelope is also proportional to the secondary arm spacing (∝ 2/λ 2 ). Considering an impingement factor Φ s Imp (= α e d ) for the growing surface, the d-s interfacial surface concentration, in reference to the total volume, is calculated as S e s = 2·Φ s Imp λ 2 · f e . For globular growth, the above equation does not apply; S e s must be equal to S e env,M . Note that Equation (17) is equally valid for columnar growth, although a superscript "e" indicating the equiaxed growth is referred to.

Conservation Equations and Solution Strategy
Sixteen independent transport quantities, n, f e , f c , f e s , f c s , p, u , u e , c , c e , c c , c e s , c c s , h , h e , h c , are to be computed. Here a velocity vector ( u ) is counted as one transport quantity, although its components (u , v , w ) are necessarily solved individually for each momentum conservation equation. Other intermediate variables can be derived from the above transport quantities.
Temperatures, T , T e , T c , are derived directly from the enthalpies, h , h e , h c . All hydrodynamic phases share a single pressure field p.
The above transport quantities are solved on the base of the hydrodynamic phases. The formulations of the conservation equations for the hydrodynamic phases are quite similar to those in the previous three-phase mixed columnar-equiaxed (non-dendritic) solidification model (Section 3.2) [64,65]. Here, special attention must be paid for the mass and species conservations.
The mass conservation equations are ∂ ∂t ( f e ρ e ) + ∇ · f e ρ e u e = M e − M ec , ∂ ∂t ( f e s ρ e ) + ∇ · f e s ρ e u e = M e ds , ∂ ∂t The species conservation equations are ∂ ∂t ( f e ρ e c e ) + ∇ · f e ρ e u e c e = C e − C ec , ∂ ∂t ( f e s ρ e c e s ) + ∇ · f e s ρ e u e c e s = C e ds , Neglecting crystal fragmentation and attachment phenomena, the mass transfer rate between the columnar and equiaxed phases M ec is equal to zero, and C ec ≡ 0. In order to close the above conservation equations, the corresponding source terms, M e , M c , M e ds , M c ds , C e , C c , C e ds , and C c ds , must be defined according to the aforementioned growth kinetics.
The mass and species transfers for equiaxed solidification are summarized in Tables 5 and 6. For dendritic growth, species transfer between the extra-dendritic melt and the equiaxed grain C e includes the transfer into the grain envelope due to the expansion of the envelope, C M e , and the diffusive flux from the inter-to the extra-dendritic melts, C D e . Back diffusion in solid dendrites is ignored here. Table 5. Mass transfer rates for equiaxed solidification.

Mass Transfer Rate
For Globular Growth For Dendritic Growth Table 6. Species transfer rates for equiaxed solidification.

Species Transfer Rate For Globular Growth For Dendritic Growth
The mass and species transfers for columnar solidification are summarized in Tables 7 and 8. For dendritic growth, the species transfer between the extra-dendritic melt and the columnar dendrite trunks C c includes species transfer into the dendrite trunk envelope due to expansion of the envelope, C M c , and the diffusive flux from the inter-to the extra-dendritic melts, C D c . Back diffusion in solid dendrites is ignored.
In Tables 7 and 8, the contribution to the mass transfer due to growth of the columnar dendrite tips, M c tip , is calculated as As shown in Tables 6 and 8, the diffusion length around the envelope l (l e or l c ) is estimated by the envelope growth velocity, D /v env , where v env represents v e env or v c env . Different methods for estimating l have been suggested [73][74][75]117,142,143]. In the present model, the D /v env formulation is preferred for its reasonable approximation and numerical simplicity [81]. The physical bounds of l dictate that it should never be larger than half of the inter-grain spacing. Two extreme cases where l = D /v env may lead to an unrealistic estimation should be avoided. One is the infinitely small envelope growth velocity; the other is the late stage of solidification when growing grains/trunks impinge upon one another. Additionally, we assume that l should not be smaller than the diffusion length of the interdendritic melt, l e d or l c d , therefore, the following corrections to l are made, where y 1 is the interdendritic spacing, which can be estimated as 2l e d for equiaxed grain, 2l c d for columnar dendrite; y 2 is the inter-grain spacing, which can be estimated as (λ 1 − d c ) for columnar growth zone and ( 3 √ 6/πn e − d e ) for equiaxed growth zone.

Equiaxed Solidification in Isolated "Representative Spherical Cell"
Equiaxed solidification in an isolated "representative spherical cell" without convection and solid movement was studied by Rappaz and Thevoz (R-T model) [73,74]. No columnar structure is involved, and solidification shrinkage is ignored (ρ = ρ s = ρ e ).
Add Equations (23) and (24), with f e c e ≡ f e d c e c + f e s c e s , we have Consider the mass and species conservations of the solid dendrite (add Equations (21) and (26)), and C e ds = c * s · M e ds , we get ∂ ∂t ( f e s c e s ) = c * s ∂ f e s ∂t . Equations (30) becomes If a "complete mixing" of interdendritic melt is considered, i.e., c e d = c e env = c * , the same species conservation equation as R-T model (Equation (6) in reference [73]) is derived: This initial model casts light on many important features of equiaxed solidification. This model was used to calculate the phase evolution during equiaxed solidification of alloy Al-5.0wt.%Si, as shown in Figure 13. Recalescence is predicted and explained. This recalescence is not due to nucleation undercooling, but is associated with the necessary undercooling to drive the equiaxed solidification. Further explanation to the modeling results of Figure 13 can be read in original papers [73,74].
shown in Figure 13. Recalescence is predicted and explained. This recalescence is not due to nucleation undercooling, but is associated with the necessary undercooling to drive the equiaxed solidification. Further explanation to the modeling results of Figure 13 can be read in original papers [73,74].
"The extension of the interdendritic melt, where complete mixing is assumed, to the spherical envelope of the grain ( 1 e env = Φ and 1 e sph = Φ ) is certainly not correct" [73]. The assumption of an ideal spherical envelope leads to an overestimation of the interdendritic melt region. Underestimation of the interdendritic melt concentration ( *  c ) and the associated tip growth undercooling (to satisfy the solute balance Equation (30b)). Due to this assumption, it is difficult to quantitatively reproduce some experimental results, especially for alloys without adding sufficient grain refiner (no inoculation and large grains). Wu et al. performed a parameter study on the influence of grain morphology parameters [82]. Findings show that the predicted total eutectic is insensitive, but the extra-and inter-dendritic eutectics are quite sensitive to the assumption of a grain envelope. By defining a suitable shape of the grain envelope, a satisfied quantitative simulation-experiment agreement can be achieved (Figure 14). "The extension of the interdendritic melt, where complete mixing is assumed, to the spherical envelope of the grain (Φ e env = 1 and Φ e sph = 1) is certainly not correct" [73]. The assumption of an ideal spherical envelope leads to an overestimation of the interdendritic melt region. Underestimation of the interdendritic melt concentration (c * ) and the associated tip growth undercooling (to satisfy the solute balance Equation (30b)). Due to this assumption, it is difficult to quantitatively reproduce some experimental results, especially for alloys without adding sufficient grain refiner (no inoculation and large grains). Wu et al. performed a parameter study on the influence of grain morphology parameters [82]. Findings show that the predicted total eutectic is insensitive, but the extra-and inter-dendritic eutectics are quite sensitive to the assumption of a grain envelope. By defining a suitable shape of the grain envelope, a satisfied quantitative simulation-experiment agreement can be achieved ( Figure 14). This initial model casts light on many important features of equiaxed solidification. This model was used to calculate the phase evolution during equiaxed solidification of alloy Al-5.0wt.%Si, as shown in Figure 13. Recalescence is predicted and explained. This recalescence is not due to nucleation undercooling, but is associated with the necessary undercooling to drive the equiaxed solidification. Further explanation to the modeling results of Figure 13 can be read in original papers [73,74].
"The extension of the interdendritic melt, where complete mixing is assumed, to the spherical envelope of the grain ( 1 e env = Φ and 1 e sph = Φ ) is certainly not correct" [73]. The assumption of an ideal spherical envelope leads to an overestimation of the interdendritic melt region. Underestimation of the interdendritic melt concentration ( *  c ) and the associated tip growth undercooling (to satisfy the solute balance Equation (30b)). Due to this assumption, it is difficult to quantitatively reproduce some experimental results, especially for alloys without adding sufficient grain refiner (no inoculation and large grains). Wu et al. performed a parameter study on the influence of grain morphology parameters [82]. Findings show that the predicted total eutectic is insensitive, but the extra-and inter-dendritic eutectics are quite sensitive to the assumption of a grain envelope. By defining a suitable shape of the grain envelope, a satisfied quantitative simulation-experiment agreement can be achieved (Figure 14). The dendritic solidification model was firstly applied to study the columnar-to-equiaxed transition by Wang and Beckermann (W-B model) [75][76][77]. Again, the flow and grain movement were not considered. Some other assumptions/considerations are as follows.
-As shown in Figure 15, 1D solidification is considered. The columnar tip front was tracked

Columnar-to-Equiaxed Transition
The dendritic solidification model was firstly applied to study the columnar-to-equiaxed transition by Wang and Beckermann (W-B model) [75][76][77]. Again, the flow and grain movement were not considered. Some other assumptions/considerations are as follows.
-As shown in Figure 15, 1D solidification is considered. The columnar tip front was tracked simply by the isotherm Tcf. This Tcf isotherm separates the calculation domain into different regions: either equiaxed growth in right side, or columnar growth in the left side. It is assumed that when the equiaxed crystals are small, they are swallowed by the columnar front and transferred into columnar structure. Therefore, the whole model system considers only three "thermodynamic phases": extra-dendritic melt ( ), solid dendrites (s) and interdendritic melt (d). This assumption leads to the consequence that the growth kinetic laws for columnar and equiaxed dendrite tips become identical, and can thus be combined into one The dendrite envelope is spherical (Φ e env = 1 and Φ e sph = 1). -Primarily, to determine the dendrite tip growth velocity, different far field liquid concentrations for equiaxed and columnar are used: the nominal composition (c 0 ) for the growth of columnar tip front, but the intrinsic extra-dendritic liquid concentration (c ) for the growth of equiaxed [77]. This part is later modified by considering the same far field liquid concentration (c ) for both modes of growth; hence, the solute field impingement between the equiaxed growth and the growth of the primary columnar tip front is considered [117]. -A two time-step scheme is used: a larger time step for thermal field, a small time step for the grain nucleation and crystal growth kinetics. Note that a similar scheme was also adopted by Rappaz and Thevoz [73,74] for their solute diffusion model.  Figure 15. Schematic of the physical problem and illustration of the multiphase approach for modeling columnar-to-equiaxed transition [77]. Reproduced from [77], with permission from Springer Nature, 2019.
With this model, some important correlations of CET, as established from classical theory and experiments [112,[148][149][150][151][152], were confirmed. For example, (1) the increase of nuclei number density promote remarkably the CET; (2) a large heat extraction from the chill mold favors the formation of columnar region; (3) the length of the columnar region increases with decreasing the nominal alloy composition; (4) lower pouring temperature (superheat) favors the formation of an equiaxed zone. A quantitative comparison of the numerically predicted CET position with some experiments was made [77,151,152]. The results rendered strong agreement, given reasonable nucleation parameters. Figure 15. Schematic of the physical problem and illustration of the multiphase approach for modeling columnar-to-equiaxed transition [77]. Reproduced from [77], with permission from Springer Nature, 2019.
With this model, some important correlations of CET, as established from classical theory and experiments [112,[148][149][150][151][152], were confirmed. For example, (1) the increase of nuclei number density promote remarkably the CET; (2) a large heat extraction from the chill mold favors the formation of columnar region; (3) the length of the columnar region increases with decreasing the nominal alloy composition; (4) lower pouring temperature (superheat) favors the formation of an equiaxed zone. A quantitative comparison of the numerically predicted CET position with some experiments was made [77,151,152]. The results rendered strong agreement, given reasonable nucleation parameters.
Drawbacks of the early study include (1) the treatment of the grain envelope as ideal sphere, although the original authors have outlined the possibility of consideration of non-spherical grain envelope [75]; and (2) exclusion of the effect of flow and grain movement. Further studies on CET with the five-phase mixed columnar-equiaxed dendritic solidification model (Section 3.3.2) by including non-spherical grain envelopes and/or flow/grain transport phenomena can be found in the references [85,153].
The equiaxed phase (hydrodynamic) volume fraction and its velocity ( f e , u e ) are tracked, as shown in Figure 16. The advancement of the columnar primary dendrite tip front (c-front) resembles the growth of the columnar structure region. The temperature development within the ingot is indicated by the liquidus isotherm (922 K), as shown in Figure 16a-c. The melt between the liquidus and the c-front is undercooled. It seems that the superheat of the entire ingot is released in less than 100 s. Growth of e-and c-phases starts to compete at the bottom of the mold. Equiaxed sedimentation and induced melt flow are highly effective during entire solidification process. At very beginning (70 s, Figure 16a), there is only thermal buoyancy flow. The liquid close to the side wall sinks while that at the ingot center rises. Further cooling (85 s) increases both f e and f c at the wall. The downward liquid stream along the side wall brings the solid crystals towards the mold bottom. The larger crystals settle at the bottom, whereas the smaller ones are carried up by the liquid due to drag force. The columnar zone extends progressively, especially at the bottom corner region. The advancing c-front "captures" the equiaxed crystals, and the captured equiaxed phase remains in the columnar mushy zone. After 103 s, the columnar structure reached the top, while equiaxed crystals are showering in the bulk. This process causes stacking of equiaxed solid ahead of the growing c-front at the bottom. Further growth of equiaxed crystals and their sedimentation increase f e ahead of the c-front in the ingot bottom region, leading to the columnar-to-equiaxed transition (CET) by the mechanism of "hard blocking". It seems that no equiaxed crystals nucleate at the ingot top any more (Figure 16e-h) in the late stage. Nevertheless, c-front continues to advance and converges towards the centerline at the ingot top.
The simulation is compared with the experiment (Figure 17). The numerically predicted phase distribution agrees to some extent with the experiment: (i) the upper ingot half is fully columnar, (ii) the center region of lower half is primarily equiaxed, (iii) the upper part of the equiaxed region in the ingot center is extended upward (iv) a mixed columnar/equiaxed structure is found in the peripherals of the lower ingot half, and (v) the mixed zone close to the mold wall is extended upward. The measured and predicted macrosegregation patterns show some similarity: (i) the measured concentration, cmix, falls in a range between 3.67-4.41 wt.% Cu compared to a predicted range of 3.739-4.392 wt.% Cu; (ii) some negative macrosegregation is found in the upper region of the ingot, (iii) the equiaxed region in the ingot core exhibits severe negative macrosegregation, (iv) the mixed columnar-equiaxed zone at the top boundary of the CET line exhibits positive macrosegregation, (v) the bottom boundary of CET contains dispersed regions of positive macrosegregation, (vi) the mixed columnar-equiaxed structure between CET line and mold wall is positively segregated, and (vii) several discrete locations of positive-negative macrosegregation exist in the upper part of the ingot. distribution agrees to some extent with the experiment: (i) the upper ingot half is fully columnar, (ii) the center region of lower half is primarily equiaxed, (iii) the upper part of the equiaxed region in the ingot center is extended upward (iv) a mixed columnar/equiaxed structure is found in the peripherals of the lower ingot half, and (v) the mixed zone close to the mold wall is extended upward. The measured and predicted macrosegregation patterns show some similarity: (i) the measured concentration, cmix, falls in a range between 3.67-4.41 wt.% Cu compared to a predicted range of 3.739-4.392 wt.% Cu; (ii) some negative macrosegregation is found in the upper region of the ingot, (iii) the equiaxed region in the ingot core exhibits severe negative macrosegregation, (iv) the mixed columnar-equiaxed zone at the top boundary of the CET line exhibits positive macrosegregation, (v) the bottom boundary of CET contains dispersed regions of positive macrosegregation, (vi) the mixed columnar-equiaxed structure between CET line and mold wall is positively segregated, and (vii) several discrete locations of positive-negative macrosegregation exist in the upper part of the ingot.  It is fundamental to state one important point: the origin of equiaxed crystals plays a key role in this kind of casting. Previous simulations (Figures 16-17) have ignored mold filling. The origin of the equiaxed crystal is modeled as heterogeneous nucleation during solidification. Nevertheless, a good simulation-experiment agreement is obtained, because a very large superheat (very high pouring temperature, 1073 K) is applied. When the casting is poured at a low superheat (pouring temperature, 973 K), the as-cast structure cannot be 'reproduced' numerically without considering the 'premature solidification' during mold filling [86,154]. This so-called premature solidification includes (i) a large amount of equiaxed nuclei originated during mold filling, and thus, (ii) an amount of solid that has formed during mold filling. A recent study [109] has confirmed the so-called 'big bang' theory of crystal origin, as proposed by Chalmers [155]. It posits that large amounts of crystals originate from the undercooled melt at the mold wall during pouring, and those crystals survive the low superheat, serving as nuclei for further growth during the subsequent solidification. A simulation by including the premature solidification during mold filling can successfully "reproduce" the experiment of a low pouring temperature as well [109].

Other Activities and Future Aspects
Some further issues of volume-averaged modeling of solidification, which were not presented in the previous examples, are briefly discussed below.
-Consideration of multi-component alloys. The volume averaging approach provides flexibility for handling the diffusion-governed crystal growth kinetics, which is fully coupled with thermodynamics, as described by Equations (3), (11), (17). This coupling is easily realized for binary alloys, It is fundamental to state one important point: the origin of equiaxed crystals plays a key role in this kind of casting. Previous simulations (Figures 16 and 17) have ignored mold filling. The origin of the equiaxed crystal is modeled as heterogeneous nucleation during solidification. Nevertheless, a good simulation-experiment agreement is obtained, because a very large superheat (very high pouring temperature, 1073 K) is applied. When the casting is poured at a low superheat (pouring temperature, 973 K), the as-cast structure cannot be 'reproduced' numerically without considering the 'premature solidification' during mold filling [86,154]. This so-called premature solidification includes (i) a large amount of equiaxed nuclei originated during mold filling, and thus, (ii) an amount of solid that has formed during mold filling. A recent study [109] has confirmed the so-called 'big bang' theory of crystal origin, as proposed by Chalmers [155]. It posits that large amounts of crystals originate from the undercooled melt at the mold wall during pouring, and those crystals survive the low superheat, serving as nuclei for further growth during the subsequent solidification. A simulation by including the premature solidification during mold filling can successfully "reproduce" the experiment of a low pouring temperature as well [109].

Other Activities and Future Aspects
Some further issues of volume-averaged modeling of solidification, which were not presented in the previous examples, are briefly discussed below. - Consideration of multi-component alloys. The volume averaging approach provides flexibility for handling the diffusion-governed crystal growth kinetics, which is fully coupled with thermodynamics, as described by Equations (3), (11), (17). This coupling is easily realized for binary alloys, but it needs a further modeling step for handling multicomponent alloys, as the solidification path (i.e., the dependency of phase evolution on the local thermal or concentration condition) is a priori unknown. Previously, a predefined solidification path to handle the ternary [156,157] or multicomponent alloy systems has been implemented [158][159][160], i.e., the f s − T or c * i − T functions are given or pre-calculated from level rule, Gulliver-Scheil or some thermodynamic models. The drawback is that the ignorance of the effect of diffusion-governed crystal growth kinetics on the solidification path would exclude its application for some alloys whose solidification path is dependent on the diffusion kinetics. Therefore, the direct incorporation of the thermodynamics with the diffusion-governed growth kinetics, i.e., coupling the microsegregation with macrosegregation calculation, is proposed by many authors [161][162][163][164][165][166]. Application of this new approach for studying the solidification of Fe-C-Mn ternary alloy [70,71] has shown the significant role of the diffusion-governed growth kinetics in the calculation of the solidification path, which in turn influences the subsequence formation of macrosegregation. -Fragmentation vs. nucleation. The early volume-averaged solidification models have considered the origin of the equiaxed grains as a continuous undercooling-dependent heterogeneous nucleation [2,55,56] or a simultaneous nucleation event [76,77,91], although more and more evidences of laboratory experiments [167][168][169][170][171][172] and industry practices [173] have confirmed the origin of equiaxed crystals by fragmentation. Hence, it becomes clearer that the solute-driven remelting of the secondary or high-order of dendrite arms is the dominant mechanism for fragmentation; and the melt flow plays an important role in the formation of fragments: (i) the flow influences the transport of the fragments and (ii) it promotes or retards the remelting of the dendrite stems through solute transportation in the interdendritic mushy zone. Therefore, Campanella et al. [174], based on Flemings' theory for the local remelting of mushy zone [175], derived an onset criterion for remelting-induced fragmentation. Unfortunately, this criterion cannot provide further quantitative information about the production rate of the number density of fragments and their initial size. Lesoult reported a very valuable experimental study for steel based on a hypothesis that stirring of molten steel ahead of the columnar solidification front would result in seeding of the liquid with dendritic fragments [173]. Presumably, those fragments were eroded from the columnar tip front. He derived from experimental data that the fragment flux, quantified by the number of crystals, eroded per unit of time and area of the columnar solidification front (cm −2 ·s −1 ), correlated with the tangential velocity of the liquid along the solidification front. This concept was later implemented into a volume-averaged solidification model by Leriche et al. to calculate the fragmentation during a mixed columnar-equiaxed solidification of steel [140]. As the size and shape of the newly eroded fragments were not available, they were assumed to be spheres with diameter of 1 µm. Inspired by this work [140] and the work of Campanella et al. [174], a local remelting-based formulation for the fragmentation was suggested and implemented in a three-phase mixed columnar-equiaxed solidification model by the current authors [133], and it is applied to calculate the as-cast structure of a Sn-10 wt.% Pb benchmark, where the equiaxed crystals come only from crystal fragmentation [176,177]. The calculation results were verified by the benchmark experiments [178,179]. -Multiphase hydrodynamics. Knowledge of multiphase hydrodynamics came mostly from the field of CFD, which were gradually introduced into the solidification research. Studying the hydrodynamic interaction between the melt flow and growing crystals or a developing mushy zone, although in its infant stage, has made significant progress. For example, there have been attempts to couple the experimental and numerical methods to determine the permeability of mushy zone. By performing direct flow calculation through the 3D dendrite network, reconstructed from series scanning of metallography on postmortem analysis of the as-cast structure [180,181] or high energy X-ray microtomography of solidifying dendrites [182][183][184] or phase field simulations [119], the permeability of the mushy zone can be derived. Very recent studies based on this method have shown that the detailed structure of the mushy zone, such as the intermetallic precipitates, influences the permeability significantly [185,186]. Another interesting study was recently made on the packing behavior of settling equiaxed dendrites [187]. A so-called packing limit (volume fraction at packing point), usually taking the value (0.637) for random close packing of equal spheres, is often used to calculate the packing equiaxed grain (characterized by the shape of grain envelope). It is found that the real value (0.2~0.6) during solidification of technical alloys can be much smaller than that value of ideal spheres. It is highly dependent on the grain morphology (shape factor) and the hydrodynamic conditions (Stokes number, or settling velocity). Additionally, to model the packing behavior of equiaxed crystals by the volume averaging approach, a special numerical scheme would be necessary to implement [188,189]. In addition to defining the packing limit, the effect of the moving direction of the crystals, either advent towards the rigid solid (already packed) region or away from it, must be taken into account as well. As the conventional volume average method does not distinguish the above effect, any crystal/crystals being brought into a volume element of the packing front (the grain volume fraction at this volume element just reaches the packing limit) are considered to be immediately "captured" by the packing front. This immediate capturing might not occur in the scenario when the flow or the motion of crystals is outward the packing front. -Incorporation of thermo-mechanics. In classical volume-averaged solidification models, the solid phases are usually treated as rigid objects, either stationary such as columnar structure and rigidly-packed crystals, or movable such as suspended equiaxed crystals. The crystals or dendrite network of the mushy zone is considered as non-deformable. Most solidification processes, however, are subject to mechanical deformation, which comes either from the imposed deformation such as the shell bulging between supporting rolls during continuous casting [58][59][60][61]190] or passively from the thermo-mechanical deformation (solidification shrinkage or thermal expansion/shrinkage of solidified region leads to deformation). The first volume-averaged solidification model incorporating thermo-mechanical deformation was introduced by Bellet et al. [191,192]. The solid phase is treated as a continuous solid skeleton, which is deformable and in strong interaction with the interdendritic flow. The deforming dendritic skeleton in turn influences the interdendritic flow, as governed by Darcy's law. This model targets the calculation of deformation-induced macrosegregation. Unfortunately, the pioneer model of Bellet et al. was applied only within a narrow range, close to the end solidification, that is to say for the cases where the lowest solid volume fraction in the calculation domain must exceed the coherent point (solid volume fraction of 0.6-0.7). This work was recently extended by Ludwig et al. to cover the entire range of solidification from pure liquid melt to complete solidification [193,194]. Note that some other thermal mechanical models (mostly finite element based) were used to analyze the surface cracks [195][196][197], hot tearing [198][199][200] and other related phenomena of different castings. In those models, however, the nature of multiphase flow during solidification is ignored or simplified; hence, the formation of structural and compositional heterogeneities (mixed columnar-equiaxed structure, macrosegregation, porosity) cannot be accounted. Further discussions on those models can be read elsewhere [201].
Future research should pay more attention to following aspects: -Enhancement of model accuracy. Pioneer work proposing the volume-averaged model is often exploratory, and the modeling results are qualitative. The aim of developing such a model is, however, to target at the engineering processes at the macroscopic scale. The modeling results must meet the accuracy criterion of engineering applications. On the one hand, some model assumptions should be looser, so that the model can represent the reality as closely as possible.
On the other hand, plenty of closure laws need to be further refined. The aforementioned studies, such as on the hydrodynamic interactions between fluid and solid phases during solidification [181][182][183][184][185][186][187][188][189][190], are to be reinforced through rigorous scientific work. -Dedicated laboratory benchmark experiments for model evaluation. Comparison of the modeling results with the engineering castings, as shown in Figures 8-10 and Figure 17, is helpful, but model validation against dedicated laboratory benchmark experiments [178,179] is more desirable. The reason is that the influencing factors for the as-cast structure of engineering castings are too complex. It is extremely difficult to correlate the as-modeled solidification result of an engineering casting to each individual model assumption, closure law or parameter. 'Curve-fitting' of the modeling result with the experiment could become a pitfall. Due to the complexity of engineering casting process, if one calculated quantity of the solidification process, e.g., a cooling curve or a segregation profile, which agrees coincidently with the experiment, it does not mean that the numerical model is fully verified. - Modeling and process parameters, and materials properties. With the enhanced modeling capacity through micro-macro coupling, more input data will be required. In additional to the classical mean of thermal physical properties and heat transfer boundary conditions, etc., one needs physical properties for the solute diffusion and growth kinetics. Furthermore, parameters describing the origin of crystals (nucleation or fragmentation) are needed as well.
-Turbulence in multiphase flow. As described in Table 2 of Section 2.1, the dispersive flux terms ( J d k ) were not considered in most multiphase solidification models. How significant are these terms in modeling of solidification, and how should they be treated in the solidification model? Further studies are needed.
The last point to be discussed about multiphase solidification models is the high calculation cost. The full 3D calculation of the steel ingot of 36 ton (Figures 8-10) with the four-phase solidification model took~4 weeks in parallel on 12 cores (2.9 GHz). One reason is the large tonnage (36 ton) of the casting and the long solidification time (7.6 h); another reason is the sophisticated, non-linear coupling of the multiple equation system. 21 transport equations were solved simultaneously. It seems that this calculation time is too costly with contemporary computer hardware. Following the study of Voller et al., [8,9], with the projection of Moore's law, using the volume-averaged models to capture the average effects of small scale phenomena in larger scale processes will continue to be a key area of research in computational physics. In comparison to other physics-based solidification models, such phase field or cellular automata, using the volume-averaged models for engineering castings, is most feasible. Maybe this will be the only physics-based solidification model which is applicable to casting processes at an engineering scale in the foreseeable future.

Summary
An introduction of the computational multi-fluid dynamics and volume average technique into the field of solidification research provides a modeling tool for the solidification processes by bridging the global multiphase transport phenomena with the solidification thermodynamics and crystal growth kinetics. The main advantage of volume-averaged modeling method, vs. phase field, cellular automata or other microstructure-oriented modeling methods, is its computational efficiency and applicability to casting processes relevance to engineering; without losing its necessary physics-based background. Thanks to the pioneer work of Beckermann and other researchers [3,4], the framework of multiphase volume-averaged modeling of solidification has been outlined. In the last two decades, the following major progress has been made: (i) the development of solidification models which consider the formation of mixed columnar-equiaxed structure; (ii) further consideration of moving equiaxed crystals and crystal dendritic morphology were refined; (iii) application of the models to analyze the formation mechanisms of macrosegregation, as-cast structure, shrinkage cavity/porosity. This review has presented several modeling examples. It is clearly demonstrated that the volume-averaged models are not only used to perform exploratory simulations, showing the probability for some solidification phenomena (achieving/strengthening solidification knowledge), but are used to calculate practical castings (e.g., Figures 9, 10 and 17) and the modeling results agree quantitatively with the as-cast results. In the foreseeable future, these applications will become more prospective with increasing computational capacity.
The following research activities in the area of the volume-averaged modeling should continue: model extension for multi-component alloy system, consideration of fragmentation as a new source of equiaxed crystals, improvement of the closure laws (multiphase hydrodynamics), and incorporation of the effect of thermo-mechanics. Further research aspects are: designing/performing dedicated laboratory benchmark experiments for model evaluation; providing more reliable modeling, processing parameters and materials properties; including turbulence models in the multiphase flow.