Mathematical Modeling of Multiphase Flow in Steel Continuous Casting

Abstract

This paper reviews multiphase flow models for continuous casting of steel from classical models to recent methods. These multiphase flow models are classified into six groups in this paper: quasi-multiphase models, multi-fluid models, moving grid methods, interface tracking methods, particle based models/methods and hybrid models. For each model, the governing equations are summarized and the inherent advantages and disadvantages are discussed. Example applications of each model are presented from previous literature, illustrating typical results and accuracy that can be obtained. The objective of this paper is to guide readers to choose an appropriate multiphase flow model for their application. Argon gas bubble effects on the flow pattern can be modeled with simple mixture models, Eulerian-Eulerian models, Multiple size group models which also track bubble size distributions, and Discrete-Phase Models (DPM). Gas pockets and slug flow, which can occur inside the nozzle are more complex, and hybrid models which combine different models together, such as Eulerian-Eulerian, level-set, Volume of Fluid (VOF), and DPM, appear promising. The shape of the slag/steel interfacial profile, level fluctuations, and slag entrainment can be modeled with free surface methods, such as moving grid, VOF, or other interface tracking methods. Particle transport and entrapment, including inclusions and gas bubbles can be added via DPM models and also require a capture criterion model. Solidification and meniscus phenomena require flow models coupled with heat transfer and solidification.

1. Introduction

In continuous casting of steel, multiphase flow arises in many different ways.¹⁾ Argon gas is often injected into the molten steel, usually through the upper tundish nozzle walls^2,3) or stopper rod tip.⁴⁾ This is beneficial to avoid clogging⁵⁾ and air aspiration.⁶⁾ However, argon bubbles affect the flow pattern in the mold^{7,8,9,10,11,12)} and can lead to the capture of small bubbles by the solidifying shell and defects in the final product.¹³⁾ Introduction of solid phase particles such as inclusions from upstream processes, released nozzle clogs, entrained slag, or re-oxidation via air aspiration, creates another class of multiphase flow behavior.^14,15,16,17) Free surfaces between the flowing molten steel and the surrounding gas, as encountered in free stream pouring,^18,19) or inside the submerged entry nozzle (SEN) where gas pockets may form,^20,21,22) represent a different kind of multiphase flow problem. Interaction between the liquid and powder mold slag layers with the molten steel at the top surface of the mold, and the changing shape of that interface, also requires a multiphase flow analysis.^{23,24,25,26,27)} Stabilizing this interface to avoid excessive surface velocities, level fluctuations and profile variations is of critical importance to avoiding slag entrainment, achieving uniform slag infiltration, and avoiding surface defects, which are important to obtaining quality in the continuous-cast product.^28,29) The solidified shell is another solid phase to consider in the analysis of mold flow as the shape of the solid shell can affect the flow pattern³⁰⁾ and the capture of particles such as bubbles.^31,32,33) The meniscus region is a place where four phases coexist, and is further complicated by the oscillatory motion of the mold and solidification of the liquid steel.^34,35,36)

Numerical simulation is a valuable tool to investigate these multiphase flow phenomena, as the extreme environment of continuous casting makes plant experiments difficult. Water model experiments are accurate for single phase flow, but have accuracy issues in multiphase flow due to the differences of material properties (such as relative density, viscosity, surface tension), contact angle, domain length, and the absence of solidification.^{30,37,38,39,40)} Several experimental efforts in water-air models have been made to overcome the differences above, such as coating walls with paraffin or fluororesin to imitate the poor wettability of liquid steel.^41,42,43) Also, a few valuable studies have used low melting point alloy–argon gas system experiments to reproduce a closer situation to real continuous casting. X-ray^4,44,45) or a transparent silica glass nozzle⁴⁶⁾ have been used to investigate argon bubble behavior in opaque liquid metal flows.

In Computational Fluid Dynamics (CFD), single phase flow is modeled by solving the continuity equation for mass conservation and a set of Navier-Stokes equations for momentum balance. These coupled partial differential equations are solved numerically using a variety of well-established numerical approaches.^47,48,49,50) To handle turbulent flows, relevant to processes such as continuous casting, the most accurate method is Direct Numerical Simulation (DNS).^51,52) DNS resolves the details of the flow using a fine computational grid, and requires no empirical turbulence model, but has a very high computational cost which increases with Reynolds number (decreasing turbulent eddy size) of the flow. The high velocities of continuous casting make it practically difficult to use DNS. Large Eddy Simulation (LES) is an alternative to simulate turbulent flow with good accuracy by using a subgrid-scale model to account for turbulent eddies smaller than the grid size.^53,54,55) Reynolds Averaged Navier-Stokes (RANS) models simulate time-averaged fields by solving additional scalar transport equations for turbulence, which enables coarse grids. They are popularly used for their low computational cost, but are less accurate in transient simulations.^56,57,58,59)

Numerous numerical methods are available to model multiphase flows. However, the variety of choices makes it difficult for users to choose a good model for their application and objectives. The purpose of this study is first to summarize the different multiphase flow models that are currently available, organizing them into groups. Secondly, the advantages and disadvantages of each model are discussed, explaining how the model can be applied, using examples in the context of modeling the continuous casting of steel.

Current multiphase flow models are classified into six groups according to their common features; each group contains several related models/methods, as presented in Fig. 1. These different models are discussed according to this framework together with example applications to the modeling of continuous casting of steel.

Fig. 1.

Overall classification of multiphase flow. (Online version in color.)

2. Quasi-multiphase Models

This early modeling approach approximates multiphase flow as an inhomogeneous mixture of liquid and gas. With the classic Algebraic-Slip Mixture (ASM) model,^60,61) only one continuity equation (Eq. (1)) and set of momentum equations (Eq. (2)) are solved for the liquid-gas mixture, and this model tracks the secondary phase (argon gas) through one extra transport equation for the gas volume fraction (Eq. (3)).   

$$\frac{\partial {\rho }_m}{\partial t}+\nabla \cdot ({\rho }_m{\mathit{u}}_m)=0$$(1)

  

$$\begin{array}{l}\frac{\partial }{\partial t}({\rho }_m{\mathit{u}}_m)+\nabla \cdot ({\rho }_m{\mathit{u}}_m{\mathit{u}}_m)=-\nabla p+\nabla \cdot ({\mu }_m(\nabla {\mathit{u}}_m+\nabla {\mathit{u}}_m{ }^T))\\ +{\rho }_m\mathit{g}+\nabla \cdot ({\alpha }_g{\rho }_g{\mathit{u}}_{dr,g}{\mathit{u}}_{dr,g}+{\alpha }_l{\rho }_l{\mathit{u}}_{dr,l}{\mathit{u}}_{dr,l})\end{array}$$(2)

  

$$\frac{\partial }{\partial t}({\alpha }_g{\rho }_g)+\nabla \cdot ({\alpha }_g{\rho }_g{\mathit{u}}_m)=-\nabla \cdot ({\alpha }_g{\rho }_g{\mathit{u}}_{dr,g})$$(3)

  

$$\begin{array}{l}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}:\mathrm{\hspace{0.17em} }{\mathit{u}}_m=\frac{{\alpha }_g{\rho }_g{\mathit{u}}_g+{\alpha }_l{\rho }_l{\mathit{u}}_l}{{\rho }_m},\mathrm{\hspace{0.17em} }\\ \mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }{\mu }_m={\alpha }_g{\mu }_g+{\alpha }_l{\mu }_l,\mathrm{\hspace{0.17em} }{\rho }_m={\alpha }_g{\rho }_g+{\alpha }_l{\rho }_l\end{array}$$(4)

The material properties (density ρ, viscosity μ) and velocity u for the mixture are mass or volume-weighted averages based on the local liquid and gas phase fractions (Eq. (4)). The averages vary spatially because the volume fraction varies in the mixture. The mixture velocity, the shared pressure field, and the gas volume-fraction field are obtained from the continuity, momentum, and volume-fraction transport equations. This also requires the drift velocity, u_dr,g and the relative velocity between the gas and liquid phases, u_gl, which are found from two further equations: an algebraic equation, Eq. (5), and a slip equation, Eq. (6).   

$${\mathit{u}}_{dr,g}={\mathit{u}}_g-{\mathit{u}}_m=\left(1-\frac{{\alpha }_g{\rho }_g}{{\rho }_m}\right){\mathit{u}}_{gl}$$(5)

  

$$\begin{array}{l}{\mathit{u}}_{gl}={\mathit{u}}_g-{\mathit{u}}_l\\ \hspace{1em}\hspace{0.17em}=\frac{{\rho }_g{d}_g{}^2}{18{\mu }_g{C}_D}\frac{({\rho }_g-{\rho }_m)}{{\rho }_g}\left(\mathit{g}-({\mathit{u}}_m\cdot \nabla ){\mathit{u}}_m-\frac{\partial {\mathit{u}}_m}{\partial t}\right)\end{array}$$(6)

The velocity of the gas phase, u_g can be recovered from the mixture and drift velocities via Eq. (5).

Note that Eq. (3) was rewritten from a classic gas volume fraction transport equation: to solve for the mixture velocity instead of the gas velocity, and to include the accompanying diffusion term that involves drift velocity u_dr. Relative to single-phase flow, this method needs to solve only one additional transport equation (for the volume fraction). Thus, this is one of the most efficient models to handle two phase flows. The main focus of this model is to find a reasonable flow field, with minimum expense to include the effect of the gas phase. This model requires prefixed bubble drag (size and shape) as an initial input, which often oversimplifies bubble behavior, especially when the regime is not bubbly flow, and is not good for studying bubble characteristics.

The Modified-Mixture model⁶²⁾ is another quasi-multiphase model, developed to simplify the ASM approach by substituting the slip equation with a relative velocity fixed to the bubble terminal velocity,^63,64,65) and solving directly for the liquid velocity. In addition, turbulent diffusion is added to the gas-phase transport equation, replacing Eq. (3) with the following:   

$$\frac{\partial {\alpha }_g}{\partial t}+\nabla \cdot ({\alpha }_g({\mathit{u}}_l+{\mathit{u}}_{ter}))=\nabla \cdot (D\nabla {\alpha }_g)$$(7)

where the diffusivity of gas bubbles, D, is calculated by $D=\frac{{\mu }_t}{\rho S{c}_t}$. This method is quite reasonable for bubbly flows where the gas bubbles reach terminal velocity quickly.

The biggest advantage of these mixture models is their low computational cost: they require only simple equations to incorporate the effects of the gas phase. Output is velocity and volume fraction of each phase and the shared pressure field. The limitations of quasi-multiphase models arise from this simplicity: defining gas behavior via a relative velocity is accurate only in bubbly flows, where the gas volume is reasonably small (e.g. < 10%), and the bubble size does not change much from its prefixed value.

In the modeling of continuous casting of steel, quasi-multiphase models have been used effectively to simulate the two-phase flow of argon and liquid steel in the nozzle and mold.^7,62,66,67) Kubo et al.⁶⁷⁾ studied ElectroMagnetic Level Accelerator (EMLA) effects on the flow using ASM for argon bubbles and free surface estimation with a prefixed bubble size. Calculated surface velocity and amount of escaping argon gas from meniscus are compared to measurements from a real caster using an immersed bar and a gas sampling hood. Liu⁷⁾ compared top surface velocity with mold level to nail board measurements of surface velocity from the operating plant and showed reasonable agreement (Fig. 2). Thomas⁶²⁾ and Bessho⁶⁶⁾ used modified mixture models and showed that the flow pattern in the mold is controlled by a balance between the inertia of the steel jet exiting the ports and the buoyancy of regions rich in gas bubbles. Huang et al.⁶⁸⁾ further estimated mold level fluctuations through a turbulent kinetic energy correlation with the modified mixture model. Hwang⁶⁹⁾ and Li⁷⁰⁾ combined ASM with MHD (magnetohydrodynamics) to model two-phase flow in the mold with EMBr (electromagnetic braking) and studied how the removal rate of inclusions tracked by Discrete Phase Model (DPM) is changed by argon injection and a magnetic field.⁷¹⁾ Kubo and coworkers⁷²⁾ compared surface velocity and argon gas floatation of ASM and DPM model (with and without stochastic effects) with real caster measurements, and found good agreement.

Fig. 2.

Validation of ASM surface velocity through nailboard measurement (a) and its gas volume fraction field (b).⁷⁾ (Online version in color.)

3. Multi-fluid Models

3.1. Eulerian-Eulerian (EE) Model

The so-called EE two-fluid model solves one continuity (Eq. (8)) and set of momentum equations (Eq. (9)) for each phase to find a different velocity field for each of the gas and liquid phases, in addition to the volume fraction field and the shared pressure field.⁷³⁾ Eulerian governing equations for the liquid phase are:   

$$\frac{\partial ({\alpha }_l{\rho }_l)}{\partial t}+\nabla \cdot ({\alpha }_l{\rho }_l{\mathit{u}}_l)=0$$(8)

  

$$\begin{array}{l}\frac{\partial ({\alpha }_l{\rho }_l{\mathit{u}}_l)}{\partial t}+\nabla \cdot ({\alpha }_l{\rho }_l{\mathit{u}}_l{\mathit{u}}_l)=-{\alpha }_l\nabla p+\nabla \cdot \\ ({\mu }_l{\alpha }_l(\nabla {\mathit{u}}_l+\nabla {\mathit{u}}_l{ }^T))+{\alpha }_l{\rho }_l\mathit{g}+{\mathit{F}}_{drag}+{\mathit{F}}_{lift}+{\mathit{F}}_{others}\end{array}$$(9)

  

$$\begin{array}{c}{\mathit{F}}_{drag}=\frac{3}{4}\frac{C_D}{d_g}{\alpha }_g{\rho }_l\left|{\mathit{u}}_g-{\mathit{u}}_l\right|({\mathit{u}}_g-{\mathit{u}}_l),\\ {\mathit{F}}_{lift}=-C_L{\alpha }_g{\rho }_l({\mathit{u}}_g-{\mathit{u}}_l)\times (\nabla \times {\mathit{u}}_l)\end{array}$$(10)

A similar set of Eulerian equations are solved for the gas phase, with opposite signs for the momentum-exchange source terms (force per unit volume), F_drag, F_lift and F_others, for the drag, lift and other forces that act between the phases. The choices of drag coefficient C_D and lift coefficient C_L depend on the bubble shape and size and the flow regime. Modeling the bubble drag coefficient is more complex than for solid particles due to its deformability, internal circulation of gas flow and its sensitivity to surfactant concentration distributions, such as associated with water contamination. Lift force can be generated by several different mechanisms, such as velocity gradients, bubble shape, and wall effects. Other momentum interactions (F_others) can be considered, such wall lubrication, and turbulent dispersion. Wall lubrication is added to consider a hydrodynamic force caused by an asymmetric drainage of surrounding liquid at a wall. This force is necessary if an accurate gas volume fraction field near the wall is required. Turbulent dispersion is added to make up for the turbulent fluctuation effects when RANS models are used.⁷⁴⁾ Other forces such as electromagnetic forces, and forces that act on or near interfaces or on particles are discussed in the DPM model section. In general, this model is more accurate than quasi-multiphase models because it allows the gas phase to behave differently than the liquid. For example, in mold flow, some EE bubbles can rise immediately after exiting the port, while other EE bubbles are carried with the jet of molten steel across the mold cavity. This behavior is not possible with a quasi-multiphase model. However, this benefit of EE models comes with the significant extra computational cost of solving another complete set of partial differential equations. Another limitation is that, like quasi-multiphase models, average bubble characteristics (shape, size) must be input to EE models either from measurements of bubble size distribution or from additional models.

Numerous studies have been conducted with EE models to investigate fluid flow in continuous casting. Nozzle design and related issues such as air aspiration, nozzle clogging, and the effects on flow in the mold have been extensively studied by Thomas and coworkers,^6,75) Li and coworkers,^8,76,77) and others.⁷⁸⁾ These EE models also require either a RANS^6,75,77,78) or LES turbulence model for the liquid phase.^8,76) The gas phase usually does not need a turbulence model because of its smaller velocities relative to the liquid phase in the continuous casting process. Liu et al. input the average bubble size measured in water model experiments into their EE model and compared the resulting flow pattern with a photograph (Fig. 3).⁸⁾

Fig. 3.

Validation of EE model using water model experiment with average bubble size from it.⁸⁾ (a) Snapshot of water model experiment, (b) Bubble size distribution from water model experiment, (c) Air volume fraction from EE model with average bubble diameter (2.88 mm).

EE models are readily extended to study particle transport, by adding a transport model for the solid-phase inclusion. Lou and Zhu⁷⁸⁾ coupled their EE model with a Population Balance Theory to simulate the behavior and size distribution of inclusions in gas-stirred ladles. Several researchers have used a DPM model to simulate particles, often together with EE models,⁷⁷⁾ as will be discussed later.

3.2. Interfacial Area Concentration (IAC) Model

A limitation of quasi-multiphase and EE models is that the local bubble size distribution cannot be estimated from the volume fraction field alone. To overcome this limitation, Ishii and Hibiki introduced interfacial area concentration as another bubble parameter.⁷⁴⁾ This model evaluates the evolution of bubble size distribution by solving an additional transport equation for the interfacial area concentration, a_int, (Eq. (11)) in addition to the volume fraction and velocity field(s) from the rest of the multiphase flow model, and an initial bubble size distribution.⁷⁴⁾   

$$\begin{array}{l}\frac{\partial a_{int}}{\partial t}+\nabla \cdot (a_{int}{\mathit{u}}_g)=S_{RC}+S_{WE}+S_{TI}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+S_{RO}+S_{SI}+S_{RV}+S_{LS}+S_{vg}\end{array}$$(11)

Once the a_int field is obtained, average bubble sizes are recovered from the gas volume fraction (α_g), through a Sauter mean diameter, $D_{sm}=\frac{6{\alpha }_g}{a_{int}}$, which is a locally averaged bubble size. If desired, the number of bubbles in a given volume can be recovered by assuming the bubble shape, (typically a sphere). The effects of different coalescence and breakup mechanisms can be modeled as source terms in RHS of Eq. (11). Breakup increases interfacial area, and coalescence decreases it. The advantage of this model is that bubble size distribution is obtained by solving just one additional transport equation, so its computational cost is lower than the MUlti-SIze Group (MUSIG) models or Population Balance Theory, which are discussed later. A disadvantage of IAC is that its accuracy depends on the coalescence and breakup modeling, which require empirically-based calibrations. This calibration has been done only for a limited number of flow conditions in water-air systems. Thus, additional work is required before this model is applied to liquid steel-argon systems such as continuous casting.

3.3. Homogeneous MUSIG Model

The Homogeneous MUlti-SIze Group (MUSIG) model is an extension of the Eulerian-Eulerian model to improve the estimation of evolution of the bubble size distribution.⁷⁹⁾ Using Population Balance Theory in a continuum framework, new modified continuity equations for the gas phase are obtained:   

$$\frac{\partial ({\rho }_g{\alpha }_g{s}_i)}{\partial t}+\nabla \cdot ({\mathit{u}}_g{\rho }_g{\alpha }_g{s}_i)={\check{S}}_i(\mathit{x},t)=B_{i,B}-D_{i,B}+B_{i,C}-D_{i,C}$$(12)

  

$${\alpha }_l+{\alpha }_g{\sum }_{i=1}^M{s}_i=1\mathrm{\hspace{0.17em} },\mathrm{\hspace{0.17em} }{\sum }_{i=1}^M{s}_i=1$$(13)

This model solves a separate continuity equation for the liquid phase and for each predefined bubble size group, i, in the gas phase. By modifying the continuity equations (Eq. (12)) with source terms ${\check{S}}_i(\mathit{x},t)$ (mass per unit volume) for coalescence and breakup, the Homogeneous MUSIG model can track the evolution of bubble size within each size group, in addition to tracking spatial variations in size evolution in the domain. Each bubble size group has its own continuity equation, so it is treated as a different fluid, which exchanges gas volume with the other groups. However, all bubble sizes share a single gas velocity field, and an average bubble size is used for the calculation of momentum interactions between the gas and liquid phases. As a consequence, the different bubble sizes are restricted to behaving in a similar manner, as in the EE model. Another shared limitation is that the initial bubble size distribution, which defines the initial volume fractions of each size group, must be predefined as an initial input. Thus, the model relies on finding information from elsewhere which greatly affects the results.

Several previous studies have investigated multiphase flow in continuous casting with Homogeneous MUSIG models. Shi^9,10) compared results from a Homogeneous MUSIG model using k–ε RANS model for turbulence and 11 discrete bubble groups (0.5 mm to 10.5 mm) with water model PIV measurements and obtained good agreement. As shown in Fig. 4, some of the flow exiting the nozzle port floats is buoyed upward, which is enabled by the large volume fraction of bubbles, especially large bubbles, in that region of the domain. This contrasts with the steel jet that traverses directly across the mold cavity, which has little buoyancy owing to the small gas fraction and mainly smaller bubbles. Liu et al.¹¹⁾ made similar flow-field comparisons, and further compared the Sauter-mean bubble sizes in 16 zones in the upper part of mold with measurements in a water model. As shown in Fig. 5,¹¹⁾ the trends and bubble sizes match well, except near the SEN.

Fig. 4.

Validation of Homogeneous MUSIG model through comparison with water model experiments.^9,10) (a) Flow picture of water model, (b) Simulation result, (c) PIV measurement.

Fig. 5.

Validation of Homogeneous MUSIG model through water model.¹¹⁾ (a) Bubble size distribution in nozzle and mold, (b) Comparison of 16 Sauter-mean diameters from 16 zones water model experiments. (Online version in color.)

3.4. Inhomogeneous MUSIG Model

A limitation of the Homogeneous MUSIG model is that all bubbles share the same velocity field regardless of size. The sophisticated Inhomogeneous model extends the Homogeneous MUSIG model by enabling different bubble size groups to have different velocity fields, by solving more than one set of Navier-Stokes equations for the gas phase, in addition to the modified continuity equations Eqs. (12), (13).⁸⁰⁾   

$$\begin{array}{l}\frac{\partial ({\rho }_g{\alpha }_{gj}{\mathit{u}}_{gj})}{\partial t}+\nabla \cdot ({\rho }_g{\alpha }_{gj}{\mathit{u}}_{gj}{\mathit{u}}_{gj})=-{\alpha }_{gj}\nabla p+\nabla \cdot \\ ({\mu }_g{\alpha }_{gj}(\nabla {\mathit{u}}_{gj}+\nabla {\mathit{u}}_{gj}{ }^T))+{\alpha }_{gj}{\rho }_g\mathit{g}+{\mathit{F}}_D+{\mathit{F}}_L+{\mathit{F}}_M\end{array}$$(14)

The Inhomogeneous MUSIG model has the advantage of allowing different bubble sizes to behave differently. However, the multiple sets of Navier-Stokes equations require greatly increased computational cost, so only two velocity groups have been used in practice.¹²⁾ Moreover, it is a challenge to apply this model with expensive turbulence models such as LES, so only less-accurate RANS models have been used with this method in practice. This model also has the limitation as other EE and homogeneous models in requiring the initial bubble size distribution as input.

A few recent applications of this methodology have been made to flow in the continuous casting nozzle and mold.¹²⁾ Various choices of RANS turbulence model and coefficients in the momentum interactions were tested and adjusted. The resulting Sauter-mean diameter results are compared with water model measurements and Fig. 6 shows overall agreement, similar to the previous simpler models.

Fig. 6.

Comparison of Inhomogeneous MUSIG model with water model experiments.¹²⁾ (a) Comparison of velocity profiles to water model, (b) Comparison of bubble sizes with water model experiments. (Online version in color.)

4. Moving Grid Methods

The Arbitrary Lagrangian Eulerian (ALE) method features movement of the nodes or cells which make up the mesh or grid of the computational domain.⁸¹⁾ The domain lies somewhere between pure Lagrangian (nodes and cells move with the fluid) and pure Eulerian (fluid moves through the mesh/grid, which remains fixed in space and time with respect to the other points in the grid, typically in the laboratory frame of reference). Thus, in ALE, the grid points move with respect to each other in order to track some aspect of the fluid, such as keeping the gas/liquid interface at the boundary between cells in multiphase flow. On that boundary, the grid/cell velocities are set to the interface velocity. Physically, this means that those cell surfaces are attached to the interface between phases and the grid deforms with time. Adaptive gridding or rezoning process(es) can be applied to recover quality of the deformed mesh, if the interface cell boundaries move too much.

Muzaferija and Peric⁸²⁾ applied this ALE concept to derive a finite-volume based method to simulate free-surface motion, applied to inviscid fluids with no surface tension, but with no discussion of rezoning. Liu²³⁾ extended this method to a practical Moving Grid method, which includes viscous effects, surface tension and rezoning of the grid within each time step. The location of the interface is found by solving exact boundary conditions at the interface (Eqs. (15), (16)) together with a single continuity equation and set of single-phase momentum equations within each phase region above and below the interface (Eqs. (8) and (9) with α_l = 1 or α_g = 1 and all F = 0) through iteration within each time step.   

$$\mathrm{K}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{\hspace{0.17em} }\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{\hspace{0.17em} }\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}:\mathrm{\hspace{0.17em} }({\mathit{u}}_l-{\mathit{u}}_g)\cdot \mathit{n}=0$$(15)

  

$$\begin{array}{l}\mathrm{D}\mathrm{y}\mathrm{n}\mathrm{a}\mathrm{m}\mathrm{i}\mathrm{c}\mathrm{\hspace{0.17em} }\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{\hspace{0.17em} }\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}:\mathrm{\hspace{0.17em} }\\ \mathit{n}(p_l-p_g)=-\mathit{n}\cdot ({\tau }_l-{\tau }_g)+\mathit{n}({\kappa }_1-{\kappa }_2)\sigma \end{array}$$(16)

Once the interface location is obtained, the mesh coordinates are adjusted to spread the deformation of the interface cells over more cells. Liu²³⁾ regridded a large portion of the domain using the energy-minimization-based rezoning method in ANSYS-FLUENT.⁸³⁾ This enables relatively large deformations of the interface relative to the cell size. The biggest advantage of this method is sharp and accurate capture of the exact interface shape, including surface tension effects. This is due to the exact boundary conditions, enabled by alignment of the interface with the cell boundaries. However, the moving grid is subject to convergence problems, and has significant computational cost. Furthermore, this methodology is limited to simple interfaces, such as the relatively calm free surface at the top of the mold, as it cannot accommodate interfaces that have excessive deformation, break apart or intersect.

Takatani et al.⁸⁴⁾ dropped the RHS terms in Eq. (16), so that vertical deformation of the grid is calculated by a simple pressure method. This was found to simulate transient flows in the continuous casting mold as good as the moving grid method.⁸⁴⁾ Further simplification treating the free surface as a flat wall and converting the pressure variations across the surface to potential energy gives a simple analytical method to easily predict free surface level variations as an uncoupled post processing step.^7,23,85) This pressure method estimates the surface height through energy conversion from pressure to potential energy.   

$$h_{steel}=\frac{p-p_{avg}}{\mathrm{\Delta }\rho g},\mathrm{\hspace{0.17em} }\mathrm{\Delta }\rho ={\rho }_{steel}-{\rho }_{ref}$$(17)

This pressure method shows good accuracy when the change of surface level is small so that surface tension effect is negligible because curvature of the top surface is small. It has been used successfully by many researchers to predict profile variations and variations in surface level in continuous casting.^7,23,85,86) Displacement of the slag layer that floats above the liquid steel surface can be considered by adjusting the reference density used in Eq. (17) from ρ_ref = 0 to a density difference (ρ_ref = c ⁎ ρ_slag), where c is an empirical constant between 0 and 1. Average surface profiles have matched well with plant measurements when c=1.⁸⁵⁾ However, the displacement effect of the slag layer is often negligible compared to simple lifting of the slag layer by the steel, especially during rapid level fluctuations, when comparison with measurements shows that c is nearly 0.⁸⁶⁾ Owing to the instant response of the calculated level to pressure changes, and its neglect of transverse flow effects, this method cannot handle complex transient behaviors and instabilities, such as wave crashing or breakup.

Theodorakakos²⁶⁾ and Panaras²⁵⁾ applied a 2D steady-state version of the moving grid method to capture the interface wave between water and oil, validated it with measurements of interface shape in a water-oil model of a continuous-casting slab mold, and suggested a minimum oil layer thickness to cover the free surface. Liu²³⁾ implemented a 3D transient version of this moving grid method with a Detached Eddy Simulation (DES) turbulence model, and applied it to simulate transient level fluctuations of the liquid steel pool in 3D. They investigated the effect of dithering on transient behavior of the free surface and mold level fluctuations (Fig. 7).

Fig. 7.

Effect of dithering on mold level fluctuation using Moving grid model.²³⁾ (a) Single-phase flow, (b) Two-phase flow (bubble diameter = 5 mm). (Online version in color.)

The Spines method, developed by Engleman for the commercial finite-element program, FIDAP,⁸⁷⁾ is a similar type of moving grid method, but uses a different procedure for rezoning. Specifically, lines of nodes are allowed to move along pathways called spines, according to constraints which are mesh- and problem-dependent. The general procedure for all of these moving grid methods is as follows: 1) define an initial interface shape and a pressure field as initial guesses, 2) solve momentum equations with the dynamic boundary condition (Eq. (16)) based on the current shape of interface and assumed pressure field, 3) solve the continuity equation and calculate mass fluxes passing through the interface, 4) move the interface according to the amount and direction of mass fluxes, 5) reconstruct rest of the mesh according to the rezoning method, 6) iterate on 1–5 until Eq. (15) is satisfied (i.e., no mass flux across the interface). This method captures the surface tension effect more accurately than other methods, but is more susceptible to convergence problems.

Rietow⁸⁸⁾ used the Spines method with finite-element analysis of turbulent flow to capture the interface between the liquid steel and slag to simulate the nail board experiment in continuous casting. This work showed how to use the lump height difference around each nail to provide an accurate estimate of surface velocity, which was validated in subsequent work.⁷⁾ Surface tension effects were found to be important for the short distance across the nail diameter.

5. Interface Tracking Methods

An important class of multiphase models uses a fixed grid for solving the continuity equation and a single set of momentum equations in the 3D domain volume, Eqs. (1), (2), combined with some other method to approximate the location of the 2D interface between phases within that domain.

5.1. Volume of Fluid (VOF) Model

The VOF model, developed by Hirt⁸⁹⁾ is a quasi-multiphase model that solves Eqs. (1), (2), (3), (4), finding the gas fraction field by solving a gas volume fraction transport equation, similar to Eq. (3), except that the terms with drift velocity in Eqs. (2) and (3) drop because there is only one velocity field. This enables the interface between gas and liquid regions to be estimated from the cells with volume fractions between 0 and 1, generally found between cells on both sides which have 0 (purely liquid) and 1 (purely gas). Surface tension forces, which affect the momentum (Eq. (2)) require a method to estimate the curvature of the calculated interface shape⁹⁰⁾   

$$S_{surf}=\int \sigma \kappa \mathit{n}\delta (\mathit{x}-{\mathit{x}}^{\prime })ds$$(18)

The big advantage of VOF is that complex interface shapes can be tracked from a simple volume fraction transport equation, without the convergence problems associated with moving grids. However, it is difficult to resolve the interface, unless a very refined mesh is used, especially for small gas volumes, such as bubbles, which require at least 8 cells across their diameter.^91,92) In addition, numerical diffusion of the gas volume fraction tends to cause smearing of the interface, which makes calculation of surface tension difficult. Higher-order solution methods and methods to redistribute the gas volume, such as geometric reconstruction, are useful to sharpen the interface.^93,94,95) To avoid spurious interface velocities, the sharp surface force (SSF) method⁹⁶⁾ has been applied successfully to bubble motion.^91,92)

Another limitation of VOF is imposed from the treatment of the advection terms in the volume fraction transport equation to avoid mass loss. This restricts the time step according to the Courant number limit,⁹⁷⁾ $\frac{u\mathrm{\Delta }t}{\mathrm{\Delta }x}+\frac{v\mathrm{\Delta }t}{\mathrm{\Delta }y}+\frac{w\mathrm{\Delta }t}{\mathrm{\Delta }z}<1$, which requires the time step to decrease with the cell size. Thus, the computational cost becomes prohibitive for modeling a real continuous caster where argon bubbles can be smaller than 3 mm, which would require cell size ~0.1 mm and time step ~10⁻⁶ sec. Moreover, bubble interactions (coalescence and break up) cannot be captured physically in VOF, as mesh resolution on the molecular scale of the interfaces would be required. Even with well-resolved bubbles, coalescence always occurs when two bubbles collide, for example.⁹⁸⁾ Thus, VOF models applied to continuous casting can simulate gas pockets and steel-slag interface motion, but not small bubbles.

Several VOF models have been successfully applied to the continuous casting of steel. Ilegbusi and Szekely⁹⁹⁾ investigated free surface modeling on steel flow in a mold, comparing results with and without deformation of the free surface meniscus using VOF. Wang et al.¹⁰⁰⁾ studied start-up of the continuous casting process, tracking the rising steel/air interface using VOF and simulating dummy bar withdrawal with time-dependent generation of new cells and a fixed grid. Some²²⁾ VOF models have simulated gas pocket formation near the top of an SEN water model. Others^{101,102,103,104)} have used VOF to investigate the shape and stability of the top free surface in a casting mold (with/without gas injection) and compared their results with water model experiments. Zhao et al. studied slag entrainment by benchmarking an water-oil experiment using VOF and simulating the formation of slag eyes, and detachment mechanisms of slag droplets depending on the casting speed.²⁹⁾

Similar investigations have been done to track motion of the molten metal/slag interface in electroslag remelting (ESR),¹⁰⁵⁾ in horizontal belt casting,¹⁰⁶⁾ and in continuous casting of steel billets.^107,108) The latter investigated the effect of nozzle misalignment on the flow pattern and surface motion in a billet mold¹⁰⁷⁾ and on shell thickness, using a solidification model.¹⁰⁸⁾ Singh et al.²⁷⁾ captured not only the water-air free surface motion, but also large air bubbles in a 1/3 scale water model using a 2D VOF model.

VOF models have also been applied to investigate the detailed transient behavior of the meniscus region during mold oscillation. Ojeda et al.¹⁰⁹⁾ initiated computational modeling of this complex problem, including the solidifying steel shell, liquid steel, liquid slag and solid slag. Ramirez-Lopez et al.^34,35) extended this approach to include solidification of the liquid steel, and evaluated shell solidification, oscillation mark formation, and slag infiltration (Fig. 8). Jonayat et al.³⁶⁾ further included temperature-dependent slag properties (viscosity, thermal conductivity and specific heat), and predicted mold temperature variations and slag consumption during mold oscillation that matched with lab and plant measurements. Recent work has extended this model to include shell solidification and using a very refined grid were able to predict solidified hook formation,^110,111) oscillation-mark formation under various conditions¹¹²⁾ and evolution of solidification at the meniscus from cast-start to steady state.¹¹³⁾

Fig. 8.

Comparison of transient behavior of meniscus region using VOF from Ramirez-Lopez et al. (2010)³⁴⁾ (right) and Ojeda et al. (2006)¹⁰⁹⁾ (left) at different periods of cycle.

5.2. Level Set (LS) Method

The Level Set method is similar to VOF, except that it estimates the interface location using a transport equation of a scalar function ϕ, that represents distance from the interface,¹¹⁴⁾   

$$\frac{\partial \varphi }{\partial t}+{\mathit{u}}_f·\nabla \varphi =0$$(19)

The interface function ϕ describes the shape of the interface through a smooth and continuous curve ϕ(x,t) = 0. Like gas fraction in VOF, the fluid properties are calculated in Level-Set models through ϕ-weighted averages. In contrast with VOF, the Level Set method involves additional reinitialization processes: smoothing of the interface and sharpening of the ϕ field. Through these processes, the Level Set method achieves more realistic material property gradients across the interface. The biggest advantage of the Level Set method is an accurate curvature calculation, which is important for surface-tension calculation. Also, the LS method is free from the Courant number limit for time step size, as it does not use donor-acceptor style advection methods. However, ϕ is not conservative through these reinitialization processes, which can lead to mass conservation problems. This method requires a fine mesh with high order numerical schemes to avoid numerical smearing of the interface. In addition, LS has the same weakness as VOF regarding inability to resolve small bubbles and bubble interactions.

The LS method has been applied, together with an energy equation solution, to investigate the behavior of bubbles rising into an inclined solidification front.¹¹⁵⁾ A porosity-dependent body force is activated when the bubble enters the mushy zone, where it may become entrapped.

5.3. Front Tracking (FT) Method

The Front Tracking method combines a standard fixed (Eulerian) mesh of the flow domain with a moving (Lagrangian) mesh to track the moving interface between phases.¹¹⁶⁾ Tryggvason¹¹⁷⁾ suggested the generation of a Lagrangian surface mesh by connecting massless markers (v_marker = u_m) on the interface, which are tracked every time step in the Eulerian mesh. The main advantage of this Front Tracking method is a sharp interface compared to other interface tracking methods, and the accompanying benefits. A disadvantage of this method is that the interaction between the fluid properties in the Eulerian domain and the Lagrangian interface needs to be handled carefully. A complex algorithm is required to reconstruct the Lagrangian surface mesh dynamically from the markers. Also, gas regions can never merge owing to the separation between each interface. Thus, additional modeling is required to handle interactions such as coalescence and breakup. As with the LS method, mass loss can happen in the process of advection of the moving interface and its reconstruction. Sufficient resolution of the Eulerian mesh is required for accuracy, and the time step size is restricted by the Courant number limit. No application to continuous casting has yet been made with this method.

6. Particle Based Models/Methods

6.1. Discrete Phase Model (DPM)

The DPM model treats the secondary phase as a discrete particle in a Lagrangian manner.^118,119) It inherently requires another model for the velocity field dominated by the primary phase, typically involving an Eulerian approach. In DPM, each particle is tracked individually as a point mass moving through the flow domain by integrating two simple ordinary differential equations: Newton’s equation of motion for the velocity (Eq. (20)) and the particle trajectory equation (Eq. (21)) for the particle position.¹²⁰⁾   

$${\rho }_p\frac{d{\mathit{v}}_i}{dt}=\sum \mathit{F}={\mathit{F}}_D+{\mathit{F}}_L+{\mathit{F}}_B+{\mathit{F}}_V+{\mathit{F}}_P$$(20)

  

$$\frac{d{\mathit{r}}_i}{dt}={\mathit{v}}_i$$(21)

The drag (F_D), lift (F_L) and buoyancy (F_B) forces are usually dominant, but the virtual mass (F_V) and pressure gradient (F_P) forces are important in transient simulations when the particles experience very curved trajectories with sharp accelerations.^121,122) These forces require separate models to account for the particle shape and the boundary layer around each particle, which are not resolved directly.

With a RANS model, the discrete random walk (DRW) method¹²³⁾ is a popular way to account for turbulent fluctuations, by adding a Gaussian-distributed random velocity component on a time averaged liquid velocity U_l calculated from RANS models.   

$${\mathit{u}}_l={\mathit{U}}_l+\zeta \sqrt{2k/3}$$(22)

Here the random number ζ is normally-distributed between −1 and +1, and is changed after a time interval chosen by the smaller of the eddy crossing time and the eddy life time.^31,83) Because isotropic turbulence is assumed, the DRW method generates equal fluctuation components in every direction. Near boundaries, the extra velocity component normal to the wall therefore becomes excessive, as discussed later.

The DPM model can be combined with several different models to find the required velocity field that is dominated by the primary molten steel phase. When the particle fraction is small, (α_g < 10⁻⁶),^124,125) such as the modeling of inclusions in the continuous-casting process, DPM may be used in a one-way coupled approach, after solving a single continuity equation and a set of momentum equations (Eqs. (8) and (9) with α_l = 1 and ignoring the momentum exchange source terms F_drag F_lift and F_others in Eq. (9)). For larger gas fractions, a more accurate one-way coupled method is to solve for the velocity field using a full EE model.¹²⁶⁾ Alternatively, the momentum exchange source terms in a single Eulerian model can be found by a two-way coupling approach, summing over all of the particles to find the net forces acting on the continuum phase per unit volume $\left(-\frac{\sum \mathit{F}}{V_{cell}}\right)$. The latter method is limited to dilute phases because the liquid volume fraction is always assumed to equal 1, which neglects the particle volume, and introduces inaccuracies for high particle volume fractions.

The intuitive DPM model is easy to apply to simulate bubble behavior with realistic distributions of size and shape, as the forces may be adjusted as appropriate for each bubble. Also, the behavior of all the bubbles can be easily and clearly visualized during post processing. As it avoids solving partial differential equations for the secondary phase, computational cost is lower, unless the number of particles tracked is very large.

A DPM approach is often applied to simulate particles (oxide inclusions, bubbles and slag entrainment) in the liquid steel in continuous casting systems. In the studies of flow pattern in the mold, Pfeiler et al.¹²⁷⁾ studied how bubbles affect the flow pattern in the mold through a comparison of one-way and two-way coupled simulations. They showed that the two-way coupling causes more spreading of bubbles and inclusions, and gives more realistic results by including the effects of large bubbles on the flow pattern. The DPM method has been used to study the effect of EMBr on two phase flow in the mold:¹²⁸⁾ bubbles or inclusions are tracked by DPM, and the magnetohydrodynamics (MHD) equations^86,129) are solved together with Eqs. (8), (9) for the liquid steel phase.^{67,86,129,130,131,132,133,134)} Kubo et al. investigated ElectroMagnetic Level Stabilizer (EMLS)¹³⁴⁾ and ElectroMagnetic Level Accelerator (EMLA) effects⁶⁷⁾ including argon bubbling with a prefixed Rosin-Rammler DPM bubble-size distribution.

Yu¹³¹⁾ found that EMBr stabilizes flow pattern, but it does not help bubbles to escape to the top surface. Toh¹³⁰⁾ compared numerical simulation results to measurements from a mercury-argon model system. The numerical result is validated by comparing a frequency of bubbles floating up to the top surface. Cho applied LES turbulence model and EMBr for the liquid steel with the DPM gas bubbles and obtained good agreement to nailboard measurements from plants.¹³²⁾ Lopez et al.²⁴⁾ used DPM with VOF: DPM was used for tracking bubbles, and VOF was used for capturing the liquid slag-liquid steel interface in 2D simulations. The simulation result was compared with measurements from a low melting point alloy (58%Bi–42%Sn)-argon-oil system.

The transport of inclusions and bubbles in the mold has been studied extensively using DPM.^{14,31,33,135,136,137,138,139)} The capture of particles into the solidification front is an important extension of these particle simulations. The simple capture criterion that particles become entrapped when they touch the solidification front is reasonable when the particle diameter is smaller than the primary dendrite arm spacing (PDAS). This was demonstrated by comparison of the results from Large Eddy Simulations of the turbulent flow together with thousands of DPM particles with measurements in both continuous casting water models and as-cast slabs.³³⁾ For large particles, however, capture is difficult due to cross flow. Thus, a more advanced capture criterion was developed for large particles, that involved a force balance on the particle balancing on three dendrite tips. This criterion involves 3 forces acting near the solidification front, in addition to the forces in Eq. (20). These forces are the lubrication force, Van-der-Walls interfacial force, and the surface-tension gradient force, which arises from the strong temperature and concentration gradients found near the solidification front, and are also called Marangoni forces.^138,140) Particles are captured only when the reaction force at the dendrite tips remains positive, so that the particle is unable to rotate away and escape back into the flow. Thomas et al.³²⁾ and Jin et al.³¹⁾ applied this methodology to model turbulent flow, particle transport and capture in a slab caster. Jin et al. firstly obtained the flow pattern using a two-way coupled DPM calculaiton, and released inclusion particles with the capture criteria on the solidified shell (Fig. 9). The sizes and locations of captured particles reasonably match experimentdal data from plant measurements. Thomas et al. compared the bubble capture results of different turbulence models (LES and RANS) and showed that LES is more accurate.³²⁾ This is because RANS models slightly overestimate particle capture due to excessive normal velocities towards the wall with the DRW method that arises from the assuming isotropic turbulence near the wall. Miki et al.¹³⁸⁾ investigated the effect of flow velocity across the solidification front on inclusion entrapment by modeling a physical laboratory experiment with solidifying steel. Several works couple a simple solidification model to calculate solid fraction for a particle entrapment model and compare the predictions with commercial steel samples.^{136,137,138,139)}

Fig. 9.

Results of bubble capture on shell faces with DPM.³¹⁾ (a) Wide face inner radius, (b) Wide face outer radius and (c) Narrow face using advanced capture criterion. (Online version in color.)

The DPM is also used to track inclusion particles to predict nozzle clogging, which is another important detrimental phenomenon in continuous casting. Zhang et al.¹⁴⁾ predicted clogging locations in the SEN, and investigated the effects of one-sided clogged nozzle on flow and particle capture. Long et al.¹³⁵⁾ studied clogging using semi-analytical solutions of pipe flow. Mohammadi-Ghaleni et al. simulated particle deposition with DPM, matched clogging in a real casting nozzle and revealed how the thickest clogging forms below the closed part of the slide-gate, where a large recirculation zone develops by flow separation.¹⁴¹⁾ These models adopt the simple criterion that particles are captured when they touch the nozzle wall.

Zhang¹⁶⁾ and Xu¹⁵⁾ studied attchment of inclusions to a single bubble by simulating flow around a sphere with small solid inclusions, which are tracked using DPM. They showed that smaller bubbles with larger inclusions increases the removal rate. Many works have implemented DPM for simulating inclusion behavior in upstream processes, such as tundish flow.^{142,143,144,145,146)}

6.2. Dense Discrete Phase Model (DDPM)

The Dense Discrete Phase Model was developed to improve the DPM when particle fractions are high.¹⁴⁷⁾ Like DPM, the DDPM governing equations uses an Eulerian model, Eqs. (8), (9) for the continuum phase. In DDPM, however, the DPM particles are used to calculate a liquid volume fraction field, found by summing the volume contribution of every particle according to its location. In addition, DDPM imposes a volume fraction limit based on packing structure which is designed to prevent the overlapping of particles in any computational cell. A limitation of DDPM is that particle sizes larger than the computational cell size are not allowed, in order to avoid negative liquid volume fractions. Instead, special treatments for large particles such as smoothing functions or cutoff volume fractions function can be used to distribute the particle volume to neighboring cells.

Furthermore, it is possible to add bubble interactions such as coalescence and breakup to DPM or DDPM through additional modeling. However, this increases the computational cost in proportion to O(n²) when interactions are considered between n particles. By applying a method to identify likely collision partners, typically found in neighbor cells, the computational cost can be dropped to O(n).¹⁴⁸⁾ Recently, Zhang et al.^149,150) included coalescence and breakup modeling based on Weber number criteria with DDPM (Fig. 10) to simulate bubbles in modeling turbulent flow in a continuous slab caster using a standard k − ε model. Lai et al. extended this model by considering bubble/inclusion interaction in a continuous casting mold with entrapment of particles into the solidified shell.¹⁵¹⁾

Fig. 10.

Simulation of bubble breakup at port by DDPM.¹⁴⁹⁾ (Online version in color.)

6.3. Discrete Element Method (DEM)

The Discrete Element Method is a Lagrangian based particle tracking method similar to DPM but without requiring a continuum-phase velocity field. DEM was originally developed to model highly packed solid particles such as granular materials (collection of solid particles) or fluidized beds.¹⁵²⁾ In the case of thick fluids, or low-density particles, a velocity field of a continuum phase could also be solved, and used to augment the forces on the particles.

This model involves tracking every particle every time step by solving the DPM equations (Eqs. (20) and (21)), together with collisions between all of the particles, which involves O(n²) extra computations. To incorporate the collision effect, a collision force F_coll is modeled in addition to other forces in Eq. (20).   

$${\mathit{F}}_{coll}={\mathit{F}}_n+{\mathit{F}}_t$$(23)

The normal collision force F_n depends on mechanical properties of the particle: hard particle collisions generate elastic collisions where kinetic energy is conserved, while soft particles experience damped, inelastic collisions.   

$$\mathrm{H}\mathrm{a}\mathrm{r}\mathrm{d}\mathrm{\hspace{0.17em} }\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{l}\mathrm{e})\mathrm{\hspace{0.17em} }\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{\hspace{0.17em} }\mathrm{s}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{\hspace{0.17em} }\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}:\mathrm{\hspace{0.17em} }{\mathit{F}}_n=(k_n\mathrm{\Delta }n){\mathit{e}}_{12}$$(24)

  

$$\begin{array}{l}\mathrm{S}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{\hspace{0.17em} }\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{l}\mathrm{e})\mathrm{\hspace{0.17em} }\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{\hspace{0.17em} }\mathrm{s}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}-\mathrm{d}\mathrm{a}\mathrm{s}\mathrm{h}\mathrm{p}\mathrm{o}\mathrm{t}\mathrm{\hspace{0.17em} }\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}:\mathrm{\hspace{0.17em} }\\ {\mathit{F}}_n=(k_n\mathrm{\Delta }n+\gamma ({\mathit{v}}_{\mathit{r}\mathit{e}\mathit{l}}\cdot {\mathit{e}}_{12})){\mathit{e}}_{12}\end{array}$$(25)

The tangential force, F_t, is due to friction (Coulomb’s law):   

$${\mathit{F}}_{\mathit{t}}={\mu }_f\left|{\mathit{F}}_{\mathit{n}}\right|$$(26)

This method is suited to simulate the delivery and behavior of solid or granulated mold powder particles on the surface of a continuous casting mold. However, no study could be found yet using this model in the modeling of continuous casting.

6.4. Smooth Particle Hydrodynamics (SPH)

Smooth Particle Hydrodynamics (SPH) is a hybrid model that combines a continuum approach with a discrete approach.^153,154) Like DPM, a set of equations is solved for every particle (Eq. (21)), but this model solves Lagrangian versions of the continuity equation (Eq. (27), Navier-Stokes equations (Eq. (28)) and an equation of state (Eq. (29)) to calculate pressure from the interaction forces among the particles.   

$$\frac{d{\rho }_i}{dt}=-{\rho }_i{\sum }_{j=1}^N{V}_j({\mathit{v}}_i-{\mathit{v}}_j)\cdot {\nabla }_{i}w({\mathit{r}}_i-{\mathit{r}}_j,h)$$(27)

  

$$\begin{array}{l}{\displaystyle \frac{d{\mathit{v}}_i}{dt}}=-{\displaystyle \frac{1}{{\rho }_i}}{\sum }_{j=1}^N{V}_j(p_j+{\rho }_i{\Pi }_{ij})\cdot \\ {\nabla }_{i}w({\mathit{r}}_i-{\mathit{r}}_j,h)+{\sum }_{j=1}^N{\displaystyle \frac{2V_j{\nu }_j{\mathit{r}}_{ij}\cdot {\nabla }_{i}w({\mathit{r}}_i-{\mathit{r}}_j,h)}{\left({\left|{\mathit{r}}_{ij}\right|}^2+{\eta }^2\right)}}{\mathit{v}}_{ij}+{\mathit{g}}_i\end{array}$$(28)

  

$$p_i=c_s{}^2({\rho }_i-{\rho }_0)+p_0$$(29)

where pressure is found with the help of an artificial viscosity function, Eq. (30), and density is interpolated with a kernel function (Eq. (31)),   

$$\begin{array}{l}{\Pi }_{ij}= {\displaystyle \frac{-\alpha c_{s}h{\mathit{v}}_{ij}\cdot {\mathit{r}}_{ij}}{{\left|{\mathit{r}}_{ij}\right|}^2+{\eta }^2}}\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }({\mathit{v}}_{ij}{\mathit{r}}_{ij}<\mathbf{0})\\ \hspace{1em}\hspace{0.5em}=0\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }({\mathit{v}}_{ij}{\mathit{r}}_{ij}>\mathbf{0})\end{array}$$(30)

  

$$\begin{array}{l}w(r,h)=\frac{1}{3\pi h^2}(q^3-6q+6),\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }(0\le q<1)\\ =\frac{1}{3\pi h^2}{(2-q)}^3,\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }(1\le q<2)=0,\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }(2\le q)\end{array}$$(31)

  

$$q=\frac{\left|\mathit{r}-{\mathit{r}}_b\right|}{h}$$(32)

The steps of this model proceed as follows: 1) guess initial velocity and location of each particle (Eq. (21)); 2) calculate density from the fluid particle distribution (Eq. (27)) using the kernel function (Eq. (31)); 3) calculate particle pressures from the density equation of state (Eq. (29)); 4) solve for particle velocities with Eq. (28) using the new pressures; 5) calculate the new positions of each particle (Eq. (21)); 6) Repeat 2–5.

The great advantage of the SPH model is that it can handle large deformation of gas/liquid interfaces in complex free surface flows, such as splashing inlet streams. Since the fluid is treated as set of particles and tracked in a Lagrangian framework, a computational mesh is not needed. The interface is a consequence of the particle distribution, without no modeling of the interface. Also, it captures continuum-level variations such as pressure and density using state functions, which integrate the contributions of all of the particles in the local region. Furthermore, SPH requires a relatively small number of particles compared to other particle-based methods as this model was developed to capture macro-scale phenomena. The disadvantage is that cost increases greatly with the number of particles, and details of the free surface, such as bubble size, cannot be resolved better than the particle size.

Recently, Natsui et al.¹⁵⁵⁾ investigated transient behavior of interfaces with in molten slag/molten matte including suspended slag droplets in copper smelting processes. Parametric studies on slag properties and droplet size concluded that the interface tension is the primary factor to decide interface shape, droplet size and settling time. Prakash et al applied SPH models to free surface flows in die casting and wheel casting of aluminum ingots.^18,19) Oxide generation by exposure to the air is modeled at the free surfaces calculated by SPH.¹⁸⁾ Solidification can be added into SPH through a temperature-dependent viscosity model.¹⁹⁾ This method has good potential for capturing free-surface flows with reactions in other systems such as transfer operations in steel continuous casting.

6.5. Lattice Boltzmann Method (LBM)

The Lattice Boltzmann method is a linealized version of the Boltzmann equation that models the flow of liquid and gas via two sets of fluid particles on a predefined lattice.¹⁵⁶⁾ It solves for the particle distribution function f, which has seven dimensions (x (space), c (velocity), t (time)).   

$$f_b^{}(\mathit{x}+{\mathit{v}}_b\mathrm{\Delta }t,t+\mathrm{\Delta }t)-f_b^{}(\mathit{x},t)=-\frac{f_b^{ }-f_{b,eq}^{ }}{\tau }+F_b$$(33)

The space and velocity fields are discretized by choosing a lattice that is denoted by DaQb, where a is the number of dimensions, and b is the number of possible directions the particle can move. Time typically marches explicitly with a time step Δt. Macroscopic variables such as velocity v and density ρ are recovered by summing the particle mass and momentum at each node. To treat multiphase problems, a special body force modeling¹⁵⁶⁾ is required which recovers interface curvature and surface tension accurately even at high density ratios.^157,158,159)   

$$F_b^{ }(\mathit{x})=-Gc{}_s^2\psi (\mathit{x}){\sum }_b\omega ({\left|{\mathit{v}}_b\right|}^2)\psi \mathrm{\hspace{0.17em} }(\mathit{x}+\mathit{v}{}_b)v_b$$(34)

The advantage of LBM is that it can simulate meso-scale phenomena such as flow in boundary layers, and can be easily parallelized. A limitation of LBM is its inherent anisotropy due to the choice of lattice. Also, this method suffers from numerical instability in multiphase flow systems with high-density ratios such as molten steel and argon gas. Furthermore, LBM has the same problem as the VOF and Level set methods regarding bubble interactions.⁹⁵⁾

Wang et al.¹⁷⁾ studied the behavior of different shapes of inclusion clusters floating in stagnant liquid using LBM. The net force acting on each solid particle is calculated by summing the forces acting on interface nodes. They found a correlation for terminal velocity as a function of the density and number of particles in the cluster. Zhang et al.¹⁶⁰⁾ showed that LBM can capture transient turbulent asymmetric single-phase flow behavior in a continuous casting mold. Inamuro et al.¹⁶¹⁾ used LBM to simulate oscillations and binary collisions of liquid droplets, and rising bubbles in a square duct and matched with theoretical predictions.

6.6. Dissipative Particle Dynamics (DPD)

Dissipative Particle Dynamics (DPD) is another particle-based method used to simulate multiphase flow. It is a coarse-grained version of Molecular Dynamics (MD), and is used to simulate meso-scale phenomena.¹⁶²⁾ Each particle in DPD contains a large group of molecules, which contrasts with MD where each particle is a single atom or molecule. Although DPD does not resolve molecular level phenomena, it can treat meso-scale phenomena. Like DEM, SPH, and MD, DPD is a Lagrangian method that solves the particle trajectory equation (Eq. (20)), together with a together with the Newton’s equation that includes interaction forces among neighboring groups of particles according to their spacing (Eq. (35)).   

$$\frac{d{\mathit{v}}_i}{dt}={\sum }_{j\ne i}{\mathit{F}}_{ij}^{\mathit{D}}+{\mathit{F}}_{ij}^{\mathit{C}}+{\mathit{F}}_{ij}^{\mathit{R}}$$(35)

where the force components include: 1) dissipative force F^D: works to reduce relative velocities between particles and generates a viscous effect, 2) conservative force F^C: from an inter-particle potential such as Lennard Jones in MD, 3) random force F^R: provides extra degrees of freedom to account for particle internal behavior. Details of these forces are given elsewhere.¹⁶³⁾

The potential advantage of DPD is its ability to capture meso-scale/molecular-based phenomena in flows with any number of phases, large deformations, and it does not have the inherent anisotropic motion problems of LBM. Its disadvantage is that DPD would require a large number of particles to even approach the phenomena captured by MD simulations, so would have great computational cost. Perhaps it could be applied to mechanistic modeling of thin films such as arising during bubble coalescence and breakup.

7. Hybrid Models

Considering that every method presented above has both advantages and disadvantages, it seems logical to improve existing models by combining them together with other models. A few such hybrid models have been developed. For example, VOF has been combined with the Level Set method (CLSVOF)¹⁶⁴⁾ to better calculate curvature. VOF has also been combined, in separate works, with Eulerian-Eulerian models, with MUSIG models, or with DPM models to better handle bubbles that are smaller than the grid size. In these hybrid models, the VOF resolves mesh scale bubbles, while the other model treats sub-grid scale bubbles. Recently, another approach that combines EE with DPM has been introduced (EEDPM).

7.1. Coupled Level Set and Volume of Fluid Model (CLSVOF)

Combining VOF with Level Set methods aims to exploit the better ability of VOF for mass conservation, with the better ability of Level set methods to track interface shape in order to calculate curvature more accurately. This is important to calculate interfacial tension for applications such as modeling the detailed behavior of individual bubbles. Recently Sussmann’s group combined these two methods to simulate free surfaces.¹⁶⁴⁾ However, CLSVOF still has several other limitations of the VOF and the Level Set methods, such as difficulty handling complex free-surface shapes involving large deformations, such as encountered during bubble coalescence and breakup, as mentioned previously. CLSVOF methods have not yet been applied to continuous casting.

7.2. NP2 Phase Model

The NP2 phase model is a hybrid of the VOF and inhomogeneous MUSIG models.¹⁶⁵⁾ As implied by its name, the NP2 model solves N + 2 continuity equations and N + 1 sets of momentum equations: one set for the continuous phase (VOF liquid and gas) and N fields for the discrete phase (inhomogeneous MUSIG model bubble-size phases). The latter inhomogeneous MUSIG model equations are to model the sub-grid-scale bubbles (bubble diameter ≤ grid size) while large gas regions and free surfaces (interface >> grid size) are captured by the VOF equations. Mass and momentum exchange source terms are needed to exchange gas between the VOF and MUSIG phases, and when the MUSIG bubbles become too large, they are converted to VOF gas. The NP2 Phase Model has advantages of both models in capturing the interface shape of large bubbles and free surfaces using an appropriate relatively-coarse grid for the VOF method, together with tracking the evolution of small sub-grid bubbles including coalescence and breakup. The main limitation of the NP2 model is its computational cost solving so many sets of continuity and momentum equations, and converging all of the mass and momentum exchanges. To relax this high computational cost, Tomiyama et al.¹⁶⁵⁾ assumed that the discrete (bubble) phases are all inviscid. Ojima et al.¹⁶⁶⁾ investigated the effects of hydrophilic particles on slurry flows through a vertical bubble column and found reasonable agreement with measurements of volume fraction fields in a water-air experiment.

7.3. Multifluid VOF Model

The Multifluid VOF model is a hybrid of the EE model and the VOF model.¹⁶⁷⁾ Actually, this model is very similar to the EE model except the continuity equations of both phases (Eq. (8)) has an additional term involving the divergence of an artificial compression velocity, $\nabla \cdot ({\mathit{u}}_{ck}{\alpha }_k(1-{\alpha }_k))$. This artificial velocity for phase k, ${\mathit{u}}_{ck}=C_{\alpha }\left|{\mathit{u}}_k\right|\frac{\nabla {\alpha }_k}{\left|\nabla {\alpha }_k\right|}$ is only activated (C_α = 1) when the cell is located on the interface between phases, and compresses the spread in the volume fraction field in the normal direction to the interface, to make the interface sharp. This interface tracking happens when the bubble size is resolvable by the mesh, which is determined with the aid of another model such as the IAC model or population balance theory. When the bubbles are too small, the artificial compression is deactivated (C_α = 0), so the phase k velocity field is modeled with the EE model.

The Multifluid VOF model does not use the reconstruction schemes found in VOF models, and it has good mass conservation. It generates a lesser accurate interface shape, but performs faster. Because of this idea, Multifluid VOF model does not lose its accuracy for the sub-grid bubbles. Limitation of this model is that it still needs a low Courant number for the interface tracking, and an artificial drag force is required to let both fluid velocities become equal at the interface.

7.4. DPM + VOF Model

Pirker et al.¹⁶⁸⁾ combined VOF with DPM for simulation of steel flow with argon bubbles in the nozzle of a continuous caster. They transferred gas volume from DPM to VOF when the volume fraction of DPM gas particles exceeds a certain threshold (bubble concentration > 0.9 × density of gas in a cell). This method enables bubbles to accumulate to become VOF gas regions, giving rise to large gas pockets and small DPM bubbles flowing together, as shown in Fig. 11. A limitation of this model is that it allows single cells to transform into VOF gas regions, which cannot be properly resolved: VOF requires at least 8 cells in each direction, meaning over 100 cells total to properly resolve a bubble shape.^91,92) Furthermore, small gas regions are more prone to smearing of volume fraction. Another serious limitation that this method shares with other VOF models is that the breakup of VOF gas regions into DPM bubbles is not considered.

Fig. 11.

Visualization of small discrete bubbles and continuous gas bubbles near port region using DPM + VOF model with their vorticity.¹⁶⁸⁾

7.5. EEDPM Model

As already discussed, the velocity field from an EE model can be supplied to a DPM model to simulate particle transport. Recently, Yang et al.^20,169) extended that approach by combining an EE model of gas-liquid flow with a DPM of gas bubble transport, that is fully stepwise-coupled between time steps. The velocity, pressure, and gas-volume-fraction fields from the EE-model is supplied to the DPM model of every bubble in the system. The spatial and time-varying DPM bubble size distribution is supplied to the EE model to enable better estimate of the interfacial areas for momentum exchange between the gas and liquid phases in the EE model. The bubbles are tracked individually as point masses in DPM, and the evolution of bubble size distributions is found with the combination of additional models of coalescence, breakup and volumetric expansion of the DPM bubbles. New DPM bubbles are created by modeling the shearing-off process that occurs at the boundaries of the gas pockets, found from the EE model. DPM bubbles which enter a gas pocket are consumed and removed from the DPM model. The continuity equation for the mass balances of the EE model is solved independently from the DPM equations, as there is no gas exchange between models, which avoids convergence problems. This hybrid model was recently applied to bubbly pipe flow experiments,^169,170) and to a lab-scale stopper-rod system using liquid Galinstan and argon gas.^21,171,172) The gas volume fraction distributions, flow pattern and bubble size distributions all compare well with the measurements.

8. Conclusions

This paper reviews the many different computational methods and models that can be used to simulate multiphase turbulent flows, as applied to the continuous casting of steel. More than 20 different models are summarized and discussed here, as classified into six different groups. Each model has advantages and disadvantages, so researchers should choose a model carefully, based on the governing phenomena important to their objectives, and according to available computational resources. Multiphase flow modeling is becoming better at tackling complex real-world phenomena, owing to the increase in model sophistication combined with more powerful computational hardware, such as parallel solution techniques. However, more research is needed to introduce fundamental phenomena to enable models to accurately predict complex multiphase flow behavior in processes such as continuous casting of steel. Several promising methods and models have not yet been applied to continuous casting problems. Meshless particle-based methods such as SPH or DPD have the potential to handle free surfaces with gas interactions such as open-stream pouring, or bubble coalescence and breakup, better than the continuum-based methods. Hybrid models are another promising approach to solve complex multiphase problems in continuous casting. The attraction of hybrid models is that they combine the advantages of two existing or new methods to overcome their disadvantages. New developments with hybrid multiphase models, including those reviewed in this paper, are expected in the near future.

Acknowledgements

The authors thank the Continuous Casting Consortium at the University of Illinois at Urbana-Champaign, and the National Science Foundation (Grant CMMI 15-63553) for support to enable this work.

Nomenclature

Symbols

a_int: interfacial area concentration

B_i,B: bubble birth rate in size group i due to breakup of larger bubbles

B_i,C: bubble birth rate in size group i due to coalescence of smaller bubbles

C_α: interface compression coefficient

C_D: drag coefficient

c_s: speed of sound

D_i,B: bubble death rate in size group i due to breakup (becomes smaller bubbles)

D_i,C: bubble death rate in size group i due to coalescence (becomes larger bubbles)

D: turbulence diffusion coefficient

D: number of space dimensions

d: diameter of bubble

d_c: damping constant

e₁₂: unit normal vector between particle 1 and 2

f: particle distribution function

F: force term for LBM

F: momentum source (force per unit volume)

G: strength of the interaction

g: gravity

h: smooth length parameter

K_ij: momentum transfer coefficient from phase i to j

k_n: normal stiffness

m: mass

M: total number of bubble size groups

N: total number of discrete particles

n: normal vector to interface

n: total number of phases

Δn: deformation of a particle in normal direction

p: pressure

Q: number of discrete velocities

r: discrete particle position vector

r_ij: relative position (r_i − r_j)

r_b: position vector for a neighbor particle

s_i: size fraction of bubble size i = $\frac{\text{a}\mathrm{\hspace{0.17em} }\text{number}\mathrm{\hspace{0.17em} }\text{of}\mathrm{\hspace{0.17em} }\text{bubble}\mathrm{\hspace{0.17em} }\text{in}\mathrm{\hspace{0.17em} }\text{size}\mathrm{\hspace{0.17em} }\text{group}\mathrm{\hspace{0.17em} }i}{\text{total}\mathrm{\hspace{0.17em} }\text{number}\mathrm{\hspace{0.17em} }\text{of}\mathrm{\hspace{0.17em} }\text{bubbles}}$

${\check{S}}_i$: continuity equation source term for size group i

S_RC: coalescence through random collision driven by turbulence eddies

S_WE: coalescence through collision due to acceleration by wake

S_TI: breakup due to turbulent eddies

S_RO: shearing off around the base rim of the cap bubbles

S_SI: breakup of large cap bubbles due to surface instability

S_RV: collision due to the difference in the bubble rise velocity

S_LS: breakup due to the laminar shear in viscous fluid

S_vg: collision due to the velocity gradient

S_Mg: momentum transfer to velocity group j due to coalescence and breakup of bubbles

S_ct: turbulent Schmidt number

t: time

u: fluid velocity field

U: time averaged fluid velocity

u_ck: compression velocity for phase k

V: volume

v: discrete particle velocity

v_ij: relative velocity (v_i − v_j)

w: kernel function

x: position

α: volume fraction

γ: damping constant

κ: curvature

μ: viscosity

μ_f: friction coefficient

Π_ij: artificial viscosity between particle i and j

ρ: density

σ: surface tension

ϕ: interface function

ν: kinematic viscosity

η: SPH model constant (η = 0.1 h)

τ: relaxation time

ψ: density dependent interaction potential

ω: weight function

Subscripts

b: bin number (for LBM)

B: buoyancy

D: drag

dr,k: drift velocity of phase k

eq: equilibrium

i: particle ID

j: pair particle of particle i

k: arbitrary phase (liquid or gas)

q: the other phase

G: gravity

g: gas

gj: gas phase velocity group j

L: lift

l: liquid

m: mixture

M: momentum transfer from different size groups by breakup or coalescence

p: particle

P: pressure gradient force

rel: relevant

t: turbulent

TD: turbulence dispersion

tan: tangential

ter: terminal

V: virtual mass

WL: wall lubrication

0: reference

Superscript

n: new

o: old

T: transpose

Math

∇: gradient

∇_i: gradient in particle i ‘s velocity direction

∇·: divergence

||: absolute value

$\frac{\partial }{\partial t}$: partial derivative with respect to time

$\frac{d}{dt}$: exact derivative with respect to time

$\frac{D}{Dt}$: material derivative

×: cross product

Σ: summation

References

• 1)   B. G.  Thomas and  L.  Zhang: ISIJ Int., 41 (2001), 1181.
• 2)   R.  Liu and  B. G.  Thomas: Metall. Mater. Trans. B, 46 (2015), 388.
• 3)   H.  Bai and  B. G.  Thomas: Metall. Mater. Trans. B, 32 (2001), 1143.
• 4)   K.  Timmel,  N.  Shevchenko,  M.  Roder,  M.  Anderhuber,  P.  Gardin,  S.  Eckert and  G.  Gerbeth: Metall. Mater. Trans. B, 46 (2015), 700.
• 5)   K.  Rackers and  B.  Thomas: 78th Steelmaking Conf., Iron and Steel Society, Warrendale, PA, (1995), 723.
• 6)   H.  Bai and  B. G.  Thomas: Metall. Mater. Trans. B, 32 (2001), 707.
• 7)   R.  Liu,  B. G.  Thomas,  J.  Sengupta,  S. D.  Chung and  M.  Trinh: ISIJ Int., 54 (2014), 2314.
• 8)   Z.  Liu,  B.  Li,  M.  Jiang and  F.  Tsukihashi: ISIJ Int., 53 (2013), 484.
• 9)   Q.  Yuan,  T.  Shi,  S. P.  Vanka and  B. G.  Thomas: Computational Modeling of Materials, Minerals and Metals, The Minerals, Metals, and Materials Society, Warrendale, PA, (2001), 491.
• 10)   T.  Shi: Master thesis, University of Illinois at Urbana Champaign, (2001), http://ccc.illinois.edu/PDF%20Files/Reports/LTV_report.pdf, (accessed 2017-05-27).
• 11)   Z. Q.  Liu,  F. S.  Qi,  B. K.  Li and  S. C. P.  Cheung: Int. J. Multiph. Flow, 79 (2016), 190.
• 12)   Z.  Liu,  F.  Qi,  B.  Li and  M.  Jiang: Metall. Mater. Trans. B, 46 (2015), 933.
• 13)   L.  Zhang,  S.  Yang,  K.  Cai,  J.  Li,  X.  Wan and  B. G.  Thomas: Metall. Mater. Trans. B, 38 (2007), 63.
• 14)   L.  Zhang,  Y.  Wang and  X.  Zuo: Metall. Mater. Trans. B, 39 (2008), 534.
• 15)   Y.  Xu,  M.  Ersson and  P. G.  Jönsson: ISIJ Int., 56 (2016), 1982.
• 16)   L.  Zhang,  J.  Aoki and  B. G.  Thomas: Metall. Mater. Trans. B, 37 (2006), 361.
• 17)   Y.  Wang,  H.  Li,  L.  Guo,  X.  Yu and  B.  Zhou: Steel Res. Int., 86 (2015), 732.
• 18)   P. W.  Cleary,  M.  Prakash,  J.  Ha,  N.  Stokes and  C.  Scott: Prog. Comput. Fluid Dyn., 7 (2007), 70.
• 19)   M.  Prakash,  P. W.  Cleary,  J.  Grandfield,  P.  Rohan and  V.  Nguyen: Prog. Comput. Fluid Dyn., 7 (2007), 101.
• 20)   H.  Yang and  B. G.  Thomas: The 7th Int. Conf. on Modelling and Simulation of Metallurgical Processes in Steelmaking (STEELSIM 2017), AIST, Warrendale, PA, (2017), 173.
• 21)   K.  Timmel,  N.  Shevchenko,  M.  Röder,  M.  Anderhuber,  P.  Gardin,  S.  Eckert and  G.  Gerbeth: Metall. Mater. Trans. B, 46 (2015), 700.
• 22)   A.  Ramos-Banderas,  R. D.  Morales,  R.  Sánchez-Pérez,  L.  García-Demedices and  G.  Solorio-Diaz: Int. J. Multiph. Flow, 31 (2005), 643.
• 23)   R.  Liu: Ph.D. thesis, University of Illinois at Urbana Champaign, (2014), https://www.ideals.illinois.edu/handle/2142/72988, (accessed 2017-05-27).
• 24)   P. E. R.  Lopez,  P. N.  Jalali,  J.  Björkvall,  U.  Sjöström and  C.  Nilsson: ISIJ Int., 54 (2014), 342.
• 25)   G. A.  Panaras,  A.  Theodorakakos and  G.  Berggeles: Metall. Mater. Trans. B, 29 (1998), 1117.
• 26)   A.  Theodorakakos and  G.  Bergeles: Metall. Mater. Trans. B, 29 (1998), 1321.
• 27)   V.  Singh,  S. K.  Dash,  J. S.  Sunitha,  S. K.  Ajmani and  A. K.  Das: ISIJ Int., 46 (2006), 210.
• 28)   L. C.  Hibbeler and  B. G.  Thomas: Iron Steel Technol., 10 (2013), 121.
• 29)   P.  Zhao,  Q.  Li,  S. B.  Kuang and  Z.  Zou: High Temp. Mater. Process., 36 (2017), 551.
• 30)   R.  Chaudhary,  B.  Rietow and  B. G.  Thomas: Materials Science and Technology Conf. and Exhibition 2009, MS&T, AIST, Warrendale, PA, (2009), 1090.
• 31)   K.  Jin,  B. G.  Thomas and  X.  Ruan: Metall. Mater. Trans. B, 47 (2016), 548.
• 32)   B. G.  Thomas,  Q.  Yuan,  S.  Mahmood,  R.  Liu and  R.  Chaudhary: Metall. Mater. Trans. B, 45 (2014), 22.
• 33)   Q.  Yuan,  B. G.  Thomas and  S. P.  Vanka: Metall. Mater. Trans. B, 35 (2004), 703.
• 34)   P. E.  Ramirez-Lopez,  P. D.  Lee and  K. C.  Mills: ISIJ Int., 50 (2010), 425.
• 35)   P. E.  Ramirez-Lopez,  P. D.  Lee,  K. C.  Mills and  B.  Santillana: ISIJ Int., 50 (2010), 1797.
• 36)   A. S. M.  Jonayat and  B. G.  Thomas: Metall. Mater. Trans. B, 45 (2014), 1842.
• 37)   S.-M.  Cho,  B. G.  Thomas and  S.-H.  Kim: ISIJ Int., 58 (2018), 1443.
• 38)   M.  Iguchi,  J.  Yoshida,  T.  Shimizu and  Y.  Mizuno: ISIJ Int., 40 (2000), 685.
• 39)   N.  Kasai and  M.  Iguchi: Tetsu-to-Hagané, 91 (2005), 847 (in Japanese).
• 40)   N.  Kasai and  M.  Iguchi: Tetsu-to-Hagané, 91 (2005), 546 (in Japanese).
• 41)   K.  Ishiguro and  M.  Iguchi: ISIJ Int., 43 (2003), 663.
• 42)   K.  Fukushi,  K.  Nogi and  M.  Iguchi: ISIJ Int., 46 (2006), 620.
• 43)   G. G.  Lee,  B. G.  Thomas and  S. H.  Kim: Met. Mater. Int., 16 (2010), 501.
• 44)   Z.  Wang,  K.  Mukai and  D.  Izu: ISIJ Int., 39 (1999), 154.
• 45)   K.  Timmel,  S.  Eckert,  G.  Gerbeth,  F.  Stefani and  T.  Wondrak: ISIJ Int., 50 (2010), 1134.
• 46)   T.  Kato,  M.  Hara,  A.  Muto,  S.  Hiraki and  M.  Kawamoto: ISIJ Int., 47 (2007), 840.
• 47)   S. P.  Vanka: Comput. Methods Appl. Mech. Eng., 55 (1986), 321.
• 48)   F. H.  Harlow and  J. E.  Welch: Phys. Fluids, 8 (1965), 2182.
• 49)   S. V.  Patankar and  D. B.  Spalding: Int. J. Heat Mass Transf., 15 (1972), 1787.
• 50)   J.  Kim and  P.  Moin: J. Comput. Phys., 59 (1985), 308.
• 51)   S. A.  Orszag and  G.  Patterson, Jr.: Phys. Rev. Lett., 28 (1972), 76.
• 52)   P.  Moin and  K.  Mahesh: Annu. Rev. Fluid Mech., 30 (1998), 539.
• 53)   J.  Smagorinsky: Mon. Weather Rev., 91 (1963), 99.
• 54)   B.  Galperin: Large Eddy Simulation of Complex Engineering and Geophysical Flows, Cambridge University Press, Cambridge, UK, (1993).
• 55)   K.  Akselvoll and  P.  Moin: J. Fluid Mech., 315 (1996), 387.
• 56)   W. P.  Jones and  B. E.  Launder: Int. J. Heat Mass Transf., 15 (1972), 301.
• 57)   D. C.  Wilcox: AIAA J., 26 (1988), 1299.
• 58)   F. R.  Menter: 23rd Fluid Dynamics, Plasmadynamics, and Lasers Conf., NASA, Ames Research Center, CA, (1993), 2906.
• 59)   W.  Rodi: Turbulence Models and Their Application in Hydraulics, CRC Press, Boca Raton, FL, (1993).
• 60)   M.  Ishii,  T. C.  Chawla and  N.  Zuber: AIChE J., 22 (1976), 283.
• 61)   M.  Manninen,  V.  Taivassalo and  S.  Kallio: On the Mixture Model for Multiphase Flow, Technical Research Centre of Finland, Espoo, (1996).
• 62)   B. G.  Thomas,  X.  Huang and  R. C.  Sussman: Metall. Mater. Trans. B, 25 (1994), 527.
• 63)   D.  Geng,  H.  Lei and  J.  He: Metall. Mater. Trans. B, 41 (2010), 234.
• 64)   S.  Soo: Fluid Dynamics of Multiphase Systems, Blaisdell Pub., Waltham, MA, (1967).
• 65)   J.  Szekely: Fluid Flow Phenomena in Metals Processing, Academic Press, New York, (1979).
• 66)   N.  Bessho,  R.  Yoda,  H.  Yamasaki,  T.  Fujii,  T.  Nozaki and  S.  Takatori: ISIJ Int., 31 (1991), 40.
• 67)   N.  Kubo,  J.  Kubota,  M.  Suzuki and  T.  Ishii: ISIJ Int., 47 (2007), 988.
• 68)   X.  Huang and  B.  Thomas: Can. Metall. Q., 37 (1998), 197.
• 69)   Y.  Hwang,  P.  Cha,  H.  Nam,  K.  Moon and  J.  Yoon: ISIJ Int., 37 (1997), 659.
• 70)   B.  Li,  T.  Okane and  T.  Umeda: Metall. Mater. Trans. B, 31 (2000), 1491.
• 71)   B.  Li and  F.  Tsukihashi: ISIJ Int., 43 (2003), 923.
• 72)   N.  Kubo,  T.  Ishii,  J.  Kubota and  N.  Aramaki: ISIJ Int., 42 (2002), 1251.
• 73)   F. H.  Harlow and  A. A.  Amsden: J. Comput. Phys., 17 (1975), 19.
• 74)   M.  Ishii and  T.  Hibiki: Thermo-fluid Dynamics of Two-phase Flow, Springer Science+Business Media, New York, (2006).
• 75)   H.  Bai and  B. G.  Thomas: Metall. Mater. Trans. B, 32 (2001), 269.
• 76)   Z.  Liu and  B.  Li: Metall. Mater. Trans. B, 48 (2017), 1833.
• 77)   Z.  Liu,  F.  Qi,  B.  Li and  M.  Jiang: J. Iron Steel Res. Int., 21 (2014), 1081.
• 78)   W.  Lou and  M.  Zhu: Metall. Mater. Trans. B, 44 (2013), 762.
• 79)   S.  Lo: Application of Population Balance to CFD Modelling of Bubbly Flow via the MUSIG model, AEA Technology, AEAT-1096, (1996).
• 80)   E.  Krepper,  D.  Lucas,  T.  Frank,  H.-M.  Prasser and  P. J.  Zwart: Nucl. Eng. Des., 238 (2008), 1690.
• 81)   C. W.  Hirt,  A. A.  Amsden and  J. L.  Cook: J. Comput. Phys., 14 (1974), 227.
• 82)   S.  Muzaferija and  M.  Peri’c: Numer. Heat Transf. B., 32 (1997), 369.
• 83)   ANSYS  Fluent Theory Manual, ANSYS Inc., Pittsburgh, PA, (2012).
• 84)   K.  Takatani,  Y.  Tanizawa,  H.  Mizukami and  K.  Nishimura: ISIJ Int., 41 (2001), 1252.
• 85)   Q.  Yuan,  B. G.  Thomas and  S.  Vanka: Metall. Mater. Trans. B, 35 (2004), 685.
• 86)   S.-M.  Cho,  B. G.  Thomas and  S.-H.  Kim: Metall. Mater. Trans. B, 47 (2016), 3080.
• 87)  FIDAP 7.0: Fluid Dynamics Analysis Package; Theory Manual, Fluid Dynamics International, Boulder, CO, (1993).
• 88)   B.  Rietow and  B. G.  Thomas: Proc. Iron and Steel Technology Conf. (AISTech 2008), AIST, Warrendale, PA, (2008).
• 89)   C. W.  Hirt and  B. D.  Nichols: J. Comput. Phys., 39 (1981), 201.
• 90)   J. U.  Brackbill,  D. B.  Kothe and  C.  Zemach: J. Comput. Phys., 100 (1992), 335.
• 91)   K.  Jin,  P.  Kumar,  S. P.  Vanka and  B. G.  Thomas: Phys. Fluids, 28 (2016), 093301.
• 92)   P.  Kumar and  S.  Vanka: Int. J. Multiph. Flow, 77 (2015), 32.
• 93)   E.  Aulisa,  S.  Manservisi,  R.  Scardovelli and  S.  Zaleski: J. Comput. Phys., 225 (2007), 2301.
• 94)   K.  Yokoi: J. Comput. Phys., 226 (2007), 1985.
• 95)   M.  van Sint Annaland,  N. G.  Deen and  J. A. M.  Kuipers: Chem. Eng. Sci., 60 (2005), 2999.
• 96)   Z.  Wang and  A. Y.  Tong: Int. J. Therm. Sci., 47 (2008), 221.
• 97)   R. H.  Pletcher,  J. C.  Tannehill and  D.  Anderson: Computational Fluid Mechanics and Heat Transfer, CRC Press, Boca Raton, FL, (2012).
• 98)   M.  van Sint Annaland,  W.  Dijkhuizen,  N. G.  Deen and  J. A. M.  Kuipers: AIChE J., 52 (2006), 99.
• 99)   O. J.  Ilegbusi and  J.  Szekely: ISIJ Int., 34 (1994), 943.
• 100)   Y.  Wang and  L.  Zhang: ISIJ Int., 50 (2010), 1777.
• 101)   P.  Ramírez-López,  L. G.  Demedices,  O.  Dávila,  R.  Sánchez-Pérez and  R. D.  Morales: Metall. Mater. Trans. B, 36 (2005), 787.
• 102)   I.  Calderon-Ramos and  R. D.  Morales: Metall. Mater. Trans. B, 47 (2016), 1866.
• 103)   E.  Torres-Alonso,  R. D.  Morales,  L. G.  Demedices,  A.  Nájera,  J.  Palafox-Ramos and  P.  Ramirez-Lopez: ISIJ Int., 47 (2007), 679.
• 104)   S.  Garcia-Hernandez,  J. d. J.  Barreto,  R. D.  Morales and  H.  Arcos-Gutierrez: ISIJ Int., 53 (2013), 809.
• 105)   Q.  Wang,  Z.  He,  B.  Li and  F.  Tsukihashi: Metall. Mater. Trans. B, 45 (2014), 2425.
• 106)   S.  Ge,  M.  Isac and  R.  Guthrie: Metall. Mater. Trans. B, 46 (2015), 2264.
• 107)   I.  Calderon-Ramos,  J. d. J.  Barreto and  S.  Garcia-Hernandez: ISIJ Int., 53 (2013), 802.
• 108)   A.  Najera-Bastida,  R. D.  Morales,  S.  Garcia-Hernandez,  E.  Torres-Alonso and  A.  Espino-Zarate: ISIJ Int., 50 (2010), 830.
• 109)   C.  Ojeda,  J.  Sengupta,  B. G.  Thomas,  J.  Barco and  J. L.  Arana: Proc. Iron and Steel Technology Conf. (AISTech 2006), AIST, Warrendale, PA, (2006), 1017.
• 110)   X.  Yan: Master thesis, University of Illinois at Urbana Champaign, (2017), https://www.ideals.illinois.edu/handle/2142/97791, (accessed 2018-04-05).
• 111)   X.  Yan,  A. S. M.  Jonayat and  B. G.  Thomas: 2016 TMS Annual Meeting & Exhibition, TMS, Warrendale, PA, (2016), 181.
• 112)   C.  Blaes: Master thesis, Colorado School of Mines, (2017), https://mountainscholar.org/handle/11124/171023, (accessed 2018-04-05).
• 113)   J.  Yang,  Z.  Cai and  M.  Zhu: ISIJ Int., 58 (2018), 299.
• 114)   S.  Osher and  J. A.  Sethian: J. Comput. Phys., 79 (1988), 12.
• 115)   Z.  Lei,  C.  Chen,  X.  Jin,  Y.-b.  Zhong and  Z.  Ren: ISIJ Int., 53 (2013), 830.
• 116)   S. O.  Unverdi and  G.  Tryggvason: J. Comput. Phys., 100 (1992), 25.
• 117)   G.  Tryggvason,  B.  Bunner,  A.  Esmaeeli,  D.  Juric,  N.  Al-Rawahi,  W.  Tauber,  J.  Han,  S.  Nas and  Y. J.  Jan: J. Comput. Phys., 169 (2001), 708.
• 118)   M. R.  Maxey: Phys. Fluids, 30 (1987), 1915.
• 119)   J. J.  Riley: Phys. Fluids, 17 (1974), 292.
• 120)   G. H.  Yeoh and  J.  Tu: Computational Techniques for Multiphase Flows, Elsevier, Amsterdam, (2009).
• 121)   R.  Clift,  J. R.  Grace and  M. E.  Weber: Bubbles, Drops, and Particles, Dover Publication, Minesota, NY, (2005).
• 122)   C. T.  Crowe,  J. D.  Schwarzkopf,  M.  Sommerfeld and  Y.  Tsuji: Multiphase Flows with Droplets and Particles, CRC Press, Boca Raton, FL, (2011).
• 123)   A.  Gosman and  E.  Loannides: J. Energy, 7 (1983), 482.
• 124)   S.  Elghobashi: Appl. Sci. Res., 52 (1994), 309.
• 125)   S.  Balachandar and  J. K.  Eaton: Annu. Rev. Fluid Mech., 42 (2010), 111.
• 126)   Z.  Liu,  B.  Li,  M.  Jiang and  F.  Tsukihashi: ISIJ Int., 54 (2014), 1314.
• 127)   C.  Pfeiler,  M.  Wu and  A.  Ludwig: Mater. Sci. Eng. A, 413–414 (2005), 115.
• 128)   J.  Szekely and  S.  Asai: Tetsu-to-Hagané, 61 (1975), 2012 (in Japanese).
• 129)   K.  Jin,  S. P.  Vanka and  B. G.  Thomas: Metall. Mater. Trans. B, 49 (2018), 1360.
• 130)   T.  Toh,  H.  Hasegawa and  H.  Harada: ISIJ Int., 41 (2001), 1245.
• 131)   H.  Yu and  M.  Zhu: ISIJ Int., 48 (2008), 584.
• 132)   S. M.  Cho,  B. G.  Thomas and  S. H.  Kim: Metall. Mater. Trans. B, 47 (2016), 3080.
• 133)   B. G.  Thomas: Steel Res. Int., 89 (2018), 170032.
• 134)   N.  Kubo,  T.  Ishii,  J.  Kubota and  T.  Ikagawa: ISIJ Int., 44 (2004), 556.
• 135)   M.  Long,  X.  Zuo,  L.  Zhang and  D.  Chen: ISIJ Int., 50 (2010), 712.
• 136)   Y.  Yin,  J.  Zhang,  S.  Lei and  Q.  Dong: ISIJ Int., 57 (2017), 2165.
• 137)   S.  Lei,  J.  Zhang,  X.  Zhao and  K.  He: ISIJ Int., 54 (2014), 94.
• 138)   Y.  Miki,  H.  Ohno,  Y.  Kishimoto and  S.  Tanaka: Tetsu-to-Hagané, 97 (2011), 423 (in Japanese).
• 139)   L. F.  Zhang and  Y. F.  Wang: JOM, 64 (2012), 1063.
• 140)   S.-M.  Lee,  S.-J.  Kim,  Y.-B.  Kang and  H.-G.  Lee: ISIJ Int., 52 (2012), 1730.
• 141)   M.  Mohammadi-Ghaleni,  M.  Asle Zaeem,  J. D.  Smith and  R.  O’Malley: Metall. Mater. Trans. B, 47 (2016), 3384.
• 142)   S.  Chang,  L.  Zhong and  Z.  Zou: ISIJ Int., 55 (2015), 837.
• 143)   K.  Chattopadhyay,  M.  Hasan,  M.  Isac and  R. I. L.  Guthrie: Metall. Mater. Trans. B, 41 (2010), 225.
• 144)   D.  Chen,  X.  Xie,  M.  Long,  M.  Zhang,  L.  Zhang and  Q.  Liao: Metall. Mater. Trans. B, 45 (2014), 392.
• 145)   S. N.  Lekakh and  D. G. C.  Robertson: ISIJ Int., 53 (2013), 622.
• 146)   S.  López-Ramirez,  J.  Palafox-Ramos,  R. D.  Morales,  J. d. J.  Barreto and  D.  Zacharias: Metall. Mater. Trans. B, 32 (2001), 615.
• 147)   N. A.  Patankar and  D. D.  Joseph: Int. J. Multiph. Flow, 27 (2001), 1685.
• 148)   J.  Capecelatro and  O.  Desjardins: J. Comput. Phys., 238 (2013), 1.
• 149)   T.  Zhang,  Z.  Luo,  H.  Zhou,  B.  Ni and  Z.  Zou: ISIJ Int., 56 (2016), 116.
• 150)   T.  Zhang,  Z. G.  Luo,  C. L.  Liu,  H.  Zhou and  Z. S.  Zou: Powder Technol., 273 (2015), 154.
• 151)   Q.  Lai,  Z.  Luo,  Q.  Hou,  T.  Zhang,  X.  Wang and  Z.  Zou: ISIJ Int., 58 (2018), 2062.
• 152)   P. A.  Cundall and  O. D. L.  Strack: Géotechnique, 29 (1979), 47.
• 153)   L. B.  Lucy: Astron. J., 82 (1977), 1013.
• 154)   Z.  Chen,  Z.  Zong,  M. B.  Liu,  L.  Zou,  H. T.  Li and  C.  Shu: J. Comput. Phys., 283 (2015), 169.
• 155)   S.  Natsui,  R.  Nashimoto,  T.  Kumagai,  T.  Kikuchi and  R. O.  Suzuki: Metall. Mater. Trans. B, 48 (2017), 1792.
• 156)   X.  Shan and  H.  Chen: Phys. Rev. E, 47 (1993), 1815.
• 157)   M. R.  Swift,  W. R.  Osborn and  J. M.  Yeomans: Phys. Rev. Lett., 75 (1995), 830.
• 158)   X.  He,  X.  Shan and  G. D.  Doolen: Phys. Rev. E, 57 (1998), R13.
• 159)   D.  Lycett-Brown and  K. H.  Luo: Phys. Rev. E, 91 (2015), 023305.
• 160)   X.  Zhang and  H.  Shen: J. Hydrodyn., 18 (2006), 431.
• 161)   T.  Inamuro,  T.  Ogata,  S.  Tajima and  N.  Konishi: J. Comput. Phys., 198 (2004), 628.
• 162)   R. D.  Groot and  P. B.  Warren: J. Chem. Phys., 107 (1997), 4423.
• 163)   A.  Tiwari and  J.  Abraham: Phys. Rev. E, 74 (2006), 056701.
• 164)   M.  Sussman and  E. G.  Puckett: J. Comput. Phys., 162 (2000), 301.
• 165)   A.  Tomiyama,  N.  Shimada and  H.  Asano: 4th ASME JSME Joint Fluids Engineering Conf., ASME, New York, (2003), 501.
• 166)   S.  Ojima,  K.  Hayashi and  A.  Tomiyama: Int. J. Multiph. Flow, 58 (2014), 154.
• 167)   K. E.  Wardle and  H. G.  Weller: Int. J. Chem. Eng., 2013 (2013), 1.
• 168)   S.  Pirker,  D.  Kahrimanovic and  S.  Schneiderbauer: Metall. Mater. Trans. B, 46 (2015), 953.
• 169)   H.  Yang,  S. P.  Vanka and  B. G.  Thomas: Proc. ASME 2017 Int. Mechanical Engineering Cong. and Exposition (IMECE), ASME, New York, (2017), 70337.
• 170)   H.  Yang,  S. P.  Vanka and  B. G.  Thomas: J. Fluids Eng., 140 (2018), 101202.
• 171)   H.  Yang,  S. P.  Vanka and  B. G.  Thomas: JOM, 70 (2018), 2148.
• 172)   H.  Yang,  S. P.  Vanka and  B. G.  Thomas: TMS Annual Meeting & Exhibition, TMS, Warrendale, PA, (2018), 119.