A Novel Multiphase Methodology Simulating Three Phase Flows in a Steel Ladle

Mixing phenomena in metallurgical steel ladles by bottom gas injection involves three phases namely, liquid molten steel, liquid slag and gaseous argon. In order to numerically solve this three-phase fluid flow system, a new approach is proposed which considers the physical nature of the gas being a dispersed phase in the liquid, while the two liquids namely, molten steel and slag are continuous phases initially separated by a sharp interface. The model was developed with the combination of two algorithms namely, IPSA (inter phase slip algorithm) where the gas bubbles are given a Eulerian approach since are considered as an interpenetrating phase in the two liquids and VOF (volume of fluid) in which the liquid is divided into two separate liquids but depending on the physical properties of each liquid they are assigned a mass fraction of each liquid. This implies that both the liquid phases (steel and slag) and the gas phase (argon) were solved for the mass balance. The Navier–Stokes conservation equations and the gas-phase turbulence in the liquid phases were solved in combination with the standard k-ε turbulence model. The mathematical model was successfully validated against flow patterns obtained experimentally using particle image velocimetry (PIV) and by the calculation of the area of the slag eye formed in a 1/17th water–oil physical model. The model was applied to an industrial ladle to describe in detail the turbulent flow structure of the multiphase system.


Introduction
Ladle steelmaking plays a key role in producing high quality steel grades.The efficiency of metallurgical reactions, such as degassing, deoxidation, and desulphurization that is, the refining capacity of a steel ladle in gas stirred ladles is related to mixing phenomena.Moreover, the homogeneity of temperature or chemical composition is also concerned with mixing phenomena in gas stirred ladles.However, the high temperatures and the opacity of liquid steel make the molten steel processing units such as ladle rather cumbersome for direct experimental measurements and virtually impossible for visual observations for fluid flow patterns.Mathematical models using numerical methods such as computational fluid dynamics (CFD) are, therefore, a valuable tool to carry out studies of steel ladles with a sufficient degree of accuracy and model prediction which makes it possible to understand fluid flow and mass transfer mechanism in a mesoscale level.Accordingly, several mathematical models have been used to investigate mixing phenomena in gas-stirred ladle systems.These models have been classified into three types: Quasi-single or pseudo-single phase models, Eulerian-Lagrangian two-phase models and Eulerian-Eulerian two-phase models.Mazumdar and Guthrie [1] compared the three mathematical modeling approaches and concluded that in terms of the velocity fields all models provide good agreement with the experimental observations.They also reported that the major limitation of current models is the unrealistic prediction of turbulence parameters such as the turbulent kinetic energy.They stated that this limitation is however compensated by a low dependency of heat and mass transfer rates on turbulence inside the ladle.Detailed reviews of modeling-based simulation studies in ladles can be found in [2][3][4].
Most mathematical models require a turbulence model to compute the turbulent kinetic energy (k) and its dissipation rate (ε) in the continuous phase.The k-ε model is currently the most popular two equation turbulence model, for three main reasons: (a) lower computational effort due to its empirical nature, (b) predicts reasonably well a fairly large range of problems and (c) the Prandtl number for ε has a reasonable value which fits the experimental data for the spread of the various entities at locations far from the walls [5].In 1972, Launder and Spalding [6] developed the most widely applied two-equation turbulence model.Modified k-ε turbulence models have been reported more recently [7][8][9] with changes in the empirical constants.Sichen [10] in a review on secondary steelmaking argues that all of the current models have an empirical nature.This empiricism is present in different simplifications and use of empirical constants, for example; (i) free surfaces are usually assumed to be flat and frictionless neglecting spout and wave formation, (ii) a direct validation with real slag/metal systems is still a complicated task, (iii) the vast majority of mathematical models are isothermal, (iv) few works have included the top slag layer.The presence of the top slag layer creates a complex multiphase system (steel/slag/gas).Jonsson and Jönsson [11] were the first to report a numerical 2D model involving a three-phase system to compute the velocity fields.They reported five times higher accuracy (with a deviation of 15%) when the free steel surface was assumed to be flat and slag phase was included in comparison when it was excluded.The interfacial forces on the gas bubbles are the drag, virtual mass, lift and turbulent dispersion forces.The relative importance of each force depends on the size of the bubbles and the degree of turbulence.The drag force is a resistance to the bubble motion relative to the surrounding liquid.The shear lift force is a force perpendicular to the bubble trajectory.The virtual mass force is proportional to half of the mass of the fluid displaced by the bubble immersion.The turbulent dispersion forces describe the effect of turbulent fluctuations in the velocity of the bubbles.Mendez et al. [12] analyzed in detail the effects of the interfacial forces as well as the role of the top slag layer on fluid flow.Their simulations in 2D proved the importance of including all the interfacial forces to properly represent fluid flow.In addition to this, they also indicated that it is mandatory to include the top slag layer in order to obtain a more realistic representation of the industrial process.
Han et al. [13] investigated mixing time in metallurgical ladles in the presence of a top slag layer, employing a Lagrangian frame of reference.They reported an increase in mixing time due to the presence of the top slag layer.Olsen and Cloete [14] investigated a three-phase system composed by top gas/liquid/bubbles, in a Eulerian-Eulerian-Lagrangian frame of reference employing the discrete phase model (DPM) to describe the bubble plume and the Eulerian VOF (volume of fluid) model for tracking the free surface, reporting the critical role of the free surface in a ladle.Based on their simulations, they concluded that an assumption of flat free surface is acceptable only if the interest is in the velocity fields, however, in terms of mixing time, a flat free surface yields lower mixing times.Therefore, if mixing time is of interest, a dynamic treatment for the free surface is required.Recently, Hoang et al. [15] studied mixing time in ladle to observe its influence on mass transfer.They observed that mixing is accomplished in a much shorter time than interphase mass transfer, specifically by two orders of magnitude, which indicated that mass transfer is the rate-limiting step.Li et al. [16] reported a three-phase model focusing on slag eye formation.They reported the advantages of two porous plugs instead of one in terms of a smaller slag eye, for the same gas flow rate.Slag emulsification due to gas injection requires a multiphase modeling approach.Sulasalmi et al. [17] investigated slag emulsification due to argon injection.The mathematical model using the VOF algorithm for a three-phase system (steel/slag/gas) and Large Eddy Simulation (LES) was used to compute turbulence.For example, to track small droplets, a dense computational grid was required in the range of 2 to 4 million computational cells so that small droplets can be tracked.Because of such a high mesh requirement, it is not possible to include the full computational domain but only a small region.A similar approach was followed by Huang et al. [18] who found that when the gas flow rate is relatively low, for example, at 12.3 NL min −1 , the dominant diameter size of entrained droplets was 2-3 mm, and all droplets were smaller than 3 mm.As the gas flow rate further increased, for example, at 16.2 NL min −1 , the droplet size distribution was broadened and the maximum size of droplets exceeded 9 mm.Li et al. [19] reported a water model and a numerical model using the VOF algorithm for a three-phase system (steel/slag/gas) that includes the computation of bubble size distribution during bottom gas injection and its effects on mixing time.The bubble motion considers the momentum transfer due the interfacial forces.Singh et al. [20] reported a 3D mathematical model using VOF for a three-phase system (steel/slag/gas) that provides information on the instantaneous interfacial area and slag eye area as a function of gas flow rate.The model uses a fine grid at the steel-slag interface and at the porous plug region and was validated by comparing with slag eye area measurements.It includes an empirical approach to compute the desulphurization rate using the sulfur partition ratio (equilibrium partition).Their studies indicate the benefits to operate with two plugs to decrease the slag eye area and increase the interfacial area.Haiyan et al. [21] developed a model using VOF and the Pressure Implicit with Splitting of Operator (PISO) algorithm to explain details of the flow structure that promote better mixing efficiency injecting different gas flow rates through the injection elements.
Although the CFD code PHOENICS contains the volume of fluid (VOF) algorithm [22] called scalar equation method (SEM) as well as the Eulerian methodology through the multi-fluid Inter penetrated slip algorithm (IPSA) [23], in order to solve multiphase flow in 3D accurately, both IPSA and SEM algorithms need to be coupled, however, they cannot be activated simultaneously.This limitation forms the motivation of this study.Schwarz [24] used an approach implemented in PHOENICS, where the IPSA technique was used to simulate bubble-liquid two-phase flow, and the algebraic slip model (ASM) (mixture model) was used to simulate the slag-metal interfacial dynamics.This gave a sharp interface (similar to the VOF method) when velocities are relatively low, but allowed mixture of the two liquid phases as entrained droplets when mixing is more intense.The method has been applied in metallurgical systems [25].An even more advanced technique was implemented using Ansys ® CFX ® by Stephens et al. [26] in a HIsmelt reduction vessel for iron production.
In this study, a three-phase mathematical model in 3D integrating both IPSA and VOF algorithms have been developed to describe the fluid flow in a gas-stirred ladle.The mathematical model is validated against experimental measurements in a 1/17th water-oil physical model by comparing the liquid flow patterns and the size of the slag eye.The model was applied to an industrial ladle to describe in detail the turbulent flow characteristics in both liquid and slag phases of the multiphase system.

Numerical Method
The three-phase fluid flow model representing a typical gas-stirred steel ladle consists of a combination of two algorithms namely, IPSA for an interpenetrating liquid and a gas phase and the split of the liquid into two liquids separated by a sharp interface (steel and slag) by using the VOF algorithm.Assumptions included in the model are: (a) constant physical properties for each phase, (b) isothermal system that is, thermal gradients are neglected, (c) flat free surface (slag-gas atmosphere), (d) turbulence is considered only in the liquid phases (fluid flow inside the gas bubbles is assumed to be laminar but turbulence affects the motion of the gas phase so turbulence is indirectly present in the average phase conservation equations through the turbulent coefficients of gas dispersion), (e) bubbles are considered to be at constant size for a given gas flow rate according to the correlation proposed by Davidson and Harrison [27] that is, no coalescence or breakage of bubbles takes place and the effect of pressure on the bubble size is also neglected, (f) gas and liquid phases are considered to be interpenetrating.

Eulerian Model Using IPSA Algorithm
The Eulerian approach considers a two-phase interpenetrated fluid flow problem where an interpenetrated discrete gas phase is dispersed into a "single" continuum liquid phase.In every point of the system the sum of the volume fractions is equal to one: where α l and α g are the volume fractions of liquid (the two liquids in our case) and gas phases respectively.Under a Eulerian frame of reference a set of equations for each phase is needed.Mass conservation or continuity equations must be satisfied for both gas (g) and liquid phases (l).In cylindrical coordinates (r, θ, z), are defined as follows: where u, v and w are the equivalent azimuthal, radial and axial velocity components, respectively; ρ represents the density.G Rg and G Rl represent the phase diffusion coefficient.The phase diffusion term models the turbulent dispersion of bubbles by random motion mechanism.Momentum conservation equations for the liquid phase: Momentum conservation equations for the gas phase: Processes 2019, 7, 175 where µ e f f , P, g, represent the effective gas viscosity, pressure and gravity constant, respectively.F rl , F zl , F θl , F rg , F zg , and F θg , represent the components r, z, and θ of the friction forces that transfer momentum between gas and liquid acting either in liquid or gas phase but the same component of both drag forces (for liquid and gas) have the same values but opposite signs for each phase.This momentum exchange, F ki , can be described by the following equation: where u (can be u, v or w) is the velocity component k (r, z, and θ), subscripts "j" and "i" represent gas or liquid and C f is the friction coefficient, obtained from the dispersed flow drag model: where V is the volume of each cell, u slip , the slip velocity and represents the difference in phase velocities, C D is the drag coefficient and D p is the bubble diameter.The drag coefficient correlations used in this work are taken from Kuo and Wallis [28] for tap water.We > 8 (12) where Re (= ρuD p µ ) is the Reynolds number and We (= ) the Weber number: Turbulence is assumed to be exclusively present in the liquids so the standard two equation k-ε turbulence model [6] was used in this simulation, where k is the turbulent kinetic energy and ε is its dissipation rate.However, it is necessary to point out that turbulence is affected by bubbles through bubble dynamics.
where G B is the turbulent generation, P k [29] is an extra turbulence promoted by the presence of bubbles, The five numbers C 1ε , C 2ε , Cµ, σ ε , σ k , are the standard constant values 1.44, 1.92, 0.09, 1.0 and 1.314, respectively [6], µ turb l is the turbulent viscosity.The turbulent viscosity and effective viscosity are taken from the k-ε model: Processes 2019, 7, 175 6 of 17

Volume of Fluid Model
The liquid phase presented in the Eulerian model is actually formed by two liquid sub-phases namely, steel and slag.The VOF algorithm solves a continuity equation with two marker variables, based on the donor-acceptor technique [30] namely, C 1 and C 3 using the marker equations as mentioned below: The scalar marker variable, C 1 , determines the steel liquid phase present in every cell, while the scalar marker variable, C 3 , determines the fraction of liquid slag present in every cell.If the value of C 1 is zero, slag is present, but if C 1 is one, steel is present in that cell.If C 1 is between 0 and 1 the physical properties in that cell are computed according to the mass fraction of every phase as for example in the case of the density: where subscripts a and e mean steel and slag phases for the overall liquid phase, respectively.Every physical property defined in each part of the system is expressed in a similar way, and consequently the mass and momentum conservation equations are solved.

Coupling between IPSA and VOF Algorithms
This model implies a transient calculation capable to track the interface between the two liquids.In order to reduce smearing at the steel-slag interface, the advanced convection scheme Superbee [31] was employed as a second order scheme for the convection term in the mass conservation equations for the trackers C 1 and C 3 , in marker equations ( 17) and (18), respectively.
The following code snippet presents a detailed sequence of calculation loops indicating the order of solution of every conservation equation of the modified solver in PHOENICS.The coupling of the IPSA and VOF algorithms is indicated through solving initially the VOF model, followed by the solution of the all velocity components of the phases and pressure through IPSA.All scalars are first iteratively solved in a time-continuous manner in the first loop (according to the number of iterations LITC of a slab during a sweep of the solution cycle encompassing the variables with indices from 12 upwards) and these scalars include the turbulent parameters and volume fraction of the phases as well as the VOF model (markers).Once these scalars are solved, all the Navier-Stokes equations are then solved in the outer loop (LSWEEP included in LSTEP) to obtain the three velocity components of the liquids and gas phases and finally, the mass conservation equations are adjusted by correcting the pressure (all phases share the same pressure).

Initial Conditions
Initially, all components of the velocity (u, v, w) for both phases are zero.Turbulence parameters are also zero for a static fluid.The entire domain is occupied by the two liquid phases (steel and a layer of slag above the steel) and no gas is originally present in the calculation (α l = 1, α l = 0).In the slag region; C 1 = 0 and C 3 = 1; in the steel region C 1 = 1 and C 3 = 0. Gas injection starts as the computation begins.

Boundary Conditions
The no-slip boundary condition is applied to the bottom and lateral walls of the ladle.This condition means that the components parallel and perpendicular to each wall are at zero speed so that all components of the liquid and gas phases are zero [6].A standard wall function is applicable since the near-wall flows are not expected to be subjected to severe pressure gradients [6,32].At the free surface zero axial gradients (or zero fluxes) of all variables are set except for the mass of gas, which is allowed to leave the system when bubbles reach the free surface.At the symmetry planes located at 0 • and 180 • in the azimuthal direction, as well as at the symmetry axis, all fluxes are zero.Finally, at the gas inlet (for single plug gas injection) or inlets (for dual plug gas injection), a mass flow rate of gas is set at the nozzle positions located at the bottom wall of the ladle.Table 1 gives a full description of the boundary conditions, based on the geometry shown in Figure 1a.

Solution Technique
The mathematical model was implemented in the CFD software PHOENICS (Version 2017, Concentration, Heat and Momentum Limited (CHAM), London, England).For validation purposes the experimental data from a scale water model were employed.Once the validation was completed, fluid flow phenomena in an industrial ladle furnace process was explored in more detail.The numerical solution involved the discretization of a 200 tons ladle into around 33,000 nodes in a polar cylindrical coordinate system with 22, 34 and 44 elements in angular radial and axial directions

Solution Technique
The mathematical model was implemented in the CFD software PHOENICS (Version 2017, Concentration, Heat and Momentum Limited (CHAM), London, England).For validation purposes the experimental data from a scale water model were employed.Once the validation was completed, fluid flow phenomena in an industrial ladle furnace process was explored in more detail.The numerical solution involved the discretization of a 200 tons ladle into around 33,000 nodes in a polar cylindrical coordinate system with 22, 34 and 44 elements in angular radial and axial directions forming a non-uniform grid obtained after a thorough sensitivity analysis of the grid.A time step size of 0.5 s was used and convergence was achieved when imbalances of all conservation equations were below 1%. Figure 1b shows the grid of the ladle.As mentioned earlier, the commercial software PHOENICS used in this study is unable to couple IPSA and the VOF methods, due to inherent limitation of the solver.As such the VOF algorithm described in Section 2.2 was programmed through a separate subroutine, written in Fortran, and linked to the solver, in order to solve the marker equations associated to the VOF algorithm and then coupled to both the IPSA and VOF models.

Experimental Work
Ladle eye: A water-engine oil 1/17th physical model was employed to measure ladle eye.The geometric characteristics are listed in Table 2.The ladle eye area was measured using a photographic technique.The area was computed using the image processing tools included in MATLAB.Formation of ladle eye is a highly dynamic process due to the discontinuous release of gas bubbles.Ladle eye was reported based on the average value of 15 photographs.Gas injection was carried out with one single nozzle located at 2/3 R. The oil thickness was constant equivalent to 4% of the height of water.The gas flow rate was 13.8 × 10 −5 m 3 /s (8.3 NL min −1 ).
Velocity patterns: Flow velocities and turbulent kinetic energy contours were obtained experimentally using particle image velocimetry (PIV, Dantec Dynamics, Denmark) with one single nozzle located at 2/3 R. PIV equipment was manufactured by Dantec Dynamics®and operated in the single frame with a time interval between pulses of 100 µs and using 30% of the total power of the laser of 10 mJ.Images were captured with a High Speed camera SpeedSense M320 (Skovlunde, Denmark) with a resolution of 1920 pixels × 1200 pixels.A total of 500 pair of photos were captured and statistically processed with the software DynamicStudio version 4.0 (Skovlunde, Denmark) to get the final results on flow patterns and turbulent intensity.Additional details can be found in [33].

Model Validation
The mathematical model was validated using experimental data on ladle eye and velocity fields.Figure 2 compares the slag eye area measured and that predicted by the model with a single plug gas injection.The experimental ladle eye was around 30%.The mathematical model also estimates a ladle eye area of around 30%.As seen, the agreement is quite good.
Processes 2019, 7, x FOR PEER REVIEW 10 of 17 Figure 2 compares the slag eye area measured and that predicted by the model with a single plug gas injection.The experimental ladle eye was around 30%.The mathematical model also estimates a ladle eye area of around 30%.As seen, the agreement is quite good.Figure 3 compares measured and computed turbulent kinetic energy contours in the longitudinal plane containing the nozzle for a single plug gas injection with a nozzle located 2/3 R and for a 1/17th physical model.Agreement is reasonably good since the maximum turbulent kinetic energy (k) is found in the plume and the order of magnitude is similar in both measured and computed contours as expected from the predictions and measurements.As seen in Figure 3, the mathematical model overestimates turbulence at the nozzle and at the free surface.The reason for the high values of k near the nozzle has been pointed out in Schwarz and Taylor [34].Bubbles are formed at the nozzle, and the diameter of the bubbles is significantly larger than the nozzle.This means that the initial velocity of bubbles at the nozzle is comparatively smaller than the gas velocity through the nozzle.This deviation can be corrected using an "effective nozzle diameter".This issue Figure 3 compares measured and computed turbulent kinetic energy contours in the longitudinal plane containing the nozzle for a single plug gas injection with a nozzle located 2/3 R and for a 1/17th physical model.Agreement is reasonably good since the maximum turbulent kinetic energy (k) is found in the plume and the order of magnitude is similar in both measured and computed contours as expected from the predictions and measurements.As seen in Figure 3, the mathematical model overestimates turbulence at the nozzle and at the free surface.The reason for the high values of k near the nozzle has been pointed out in Schwarz and Taylor [34].Bubbles are formed at the nozzle, and the diameter of the bubbles is significantly larger than the nozzle.This means that the initial velocity of bubbles at the nozzle is comparatively smaller than the gas velocity through the nozzle.This deviation can be corrected using an "effective nozzle diameter".This issue has been neglected in this study as at the nozzle, the gas holdup is very high so the local value of k computed does not seem to have any physical significance.At the free surface, the model does not allow the deformation of the free surface and therefore there is an overestimation that is explained by the flat free surface assumption.In spite of the deviation, it can be concluded that the mathematical model offers a good prediction of the experimental data with a maximum deviation of around 10% and if the over predicted zones of high k at the nozzle and at the free surface were eliminated, the predicted and measured k maps would have agreed much better.The Euler-Euler part of the model (IPSA) presented in this study has been extensively validated against experimental results of liquid velocity, gas holdup, and trace evolution in mixing experiments as indicated in [35].The combination of the IPSA and VOF algorithms allowing the presence of two liquids is validated in this study through the slag eye area.

Flow Patterns of an Industrial Ladle
The mathematical model developed was applied to the industrial size ladle.Table 3 shows important physical and geometrical information of the 200 tons ladle from a steel plant in México used to run the calculations.The model was employed to get a full characterization of flow patterns, in terms of velocity fields, turbulence parameters and so on.Two simulations were performed.The first case includes one nozzle located at 2/3 R, a gas flow rate of 5.5 × 10 −3 m 3 /s and a slag thickness of 4% H.The second simulation includes two nozzles located at 2/3 R, an angle of 180° between the nozzles, gas flow rates of 5.5 × 10 −3 m 3 /s in each nozzle and the same slag thickness.

Flow Patterns of an Industrial Ladle
The mathematical model developed was applied to the industrial size ladle.Table 3 shows important physical and geometrical information of the 200 tons ladle from a steel plant in México used to run the calculations.The model was employed to get a full characterization of flow patterns, in terms of velocity fields, turbulence parameters and so on.Two simulations were performed.The first case includes one nozzle located at 2/3 R, a gas flow rate of 5.5 × 10 −3 m 3 /s and a slag thickness of 4% H.The second simulation includes two nozzles located at 2/3 R, an angle of 180 • between the nozzles, gas flow rates of 5.5 × 10 −3 m 3 /s in each nozzle and the same slag thickness.Figure 4 shows the two-liquid phase distribution in the ladle for both simulations.Only half of the ladle is shown due to symmetry.The two liquid phases are immiscible in each other and ideally a line should describe the interface.However, all mathematical models report some degree of error to represent an interface, called numerical diffusion.A finer grid around the interface, and higher-order convective schemes were used to mitigate the numerical diffusion.The results in this figure show that numerical diffusion is low, which is desired and is an indication of the accuracy of this model.
Figure 5a shows the liquid velocity fields for the single plug gas injection (i.e., one nozzle) case.The gas injected at the bottom forms bubbles that rise to the top surface.The two-phase region is called the plume.The motion of the bubbles drags liquid from the bottom to the top surface.On the top surface, gas escapes from the liquid opening the slag layer, forming the "slag eye".Liquid molten steel collides with the top slag layer and is pushed backwards, forming a recirculation loop on each side of the plume.The position of the central part of each recirculation loop changes depending on the gas flow rate and slag thickness.Figure 6b shows the predicted velocity fields when two nozzles at 2/3 R are considered for the simulation.Here it can be seen that, due to the same effects explained above concerning the dragging of the liquid by the ascending currents of bubbles and to the collision of the liquid steel with the top slag layer, four recirculation loops are formed, each one on each side of the two plumes.Regarding the velocity fields in the slag layer, in Figure 5a it is observed in the case of one nozzle, that liquid steel flows away in opposite direction from the ascending plume when the top slag layer is reached and the slag appears to behave in a similar way flowing away from the spout to the opposite wall.Figure 5b shows a more complex velocity field in the slag layer when two nozzles are considered in the simulation, including the presence, in the slag layer region located between the plumes, of two recirculation loops.
a line should describe the interface.However, all mathematical models report some degree of error to represent an interface, called numerical diffusion.A finer grid around the interface, and higherorder convective schemes were used to mitigate the numerical diffusion.The results in this figure show that numerical diffusion is low, which is desired and is an indication of the accuracy of this model.Figure 5a shows the liquid velocity fields for the single plug gas injection (i.e., one nozzle) case.The gas injected at the bottom forms bubbles that rise to the top surface.The two-phase region is called the plume.The motion of the bubbles drags liquid from the bottom to the top surface.On the top surface, gas escapes from the liquid opening the slag layer, forming the "slag eye".Liquid molten steel collides with the top slag layer and is pushed backwards, forming a recirculation loop on each side of the plume.The position of the central part of each recirculation loop changes depending on the gas flow rate and slag thickness.Figure 6b shows the predicted velocity fields when two nozzles at 2/3 R are considered for the simulation.Here it can be seen that, due to the same effects explained above concerning the dragging of the liquid by the ascending currents of bubbles and to the collision of the liquid steel with the top slag layer, four recirculation loops are formed, each one on each side of the two plumes.Regarding the velocity fields in the slag layer, in Figure 5a it is observed in the case of one nozzle, that liquid steel flows away in opposite direction from the ascending plume when the top slag layer is reached and the slag appears to behave in a similar way flowing away from the spout to the opposite wall.Figure 5b shows a more complex velocity field in the slag layer when two   Figure 6 shows more details of the velocity patterns in the slag layer.In Figure 6a, with single plug gas injection located at 1/2 R, it can be observed that when the slag flows away from the spout it reaches the reactor wall, turns to the left and follows through the peripheral region of the slag layer located near the wall, returning to the spout region.Figure 6b, with dual plug gas injection, shows the prediction of more complex flow patterns in the slag layer due to the interactions of the slag currents leaving each one of the two plumes regions and colliding with each other in the symmetry axis zone, and colliding with the reactor walls.
Figure 7 shows the volume gas fraction in the plumes or gas hold up for the two cases taken as examples of the model output (single gas injection with a nozzle at 1/2 R and dual injection with nozzles located at 2/3 R).The plume region is clearly determined from these plots, which forms a Figure 6 shows more details of the velocity patterns in the slag layer.In Figure 6a, with single plug gas injection located at 1/2 R, it can be observed that when the slag flows away from the spout it reaches the reactor wall, turns to the left and follows through the peripheral region of the slag layer located near the wall, returning to the spout region.Figure 6b, with dual plug gas injection, shows the prediction of more complex flow patterns in the slag layer due to the interactions of the slag currents leaving each one of the two plumes regions and colliding with each other in the symmetry axis zone, and colliding with the reactor walls.
Figure 7 shows the volume gas fraction in the plumes or gas hold up for the two cases taken as examples of the model output (single gas injection with a nozzle at 1/2 R and dual injection with nozzles located at 2/3 R).The plume region is clearly determined from these plots, which forms a plume in the shape of a cylinder that shows a lateral dispersion at the free surface.The shape and size of the plume is an important parameter to estimate the degree of mixing of liquid steel.As the diameter of this plume increases, it is expected to get higher mixing as noted by Conejo et al. [36].It is noticed that the gas fraction is higher in the central axis of each plume.This is because there are no diffusion barriers for the gas phase it would be expected a lower gas concentration in the central plume when the free surface is reached.The high value of gas hold up in this figure requires further analysis.Figure 8; Figure 9 shows characteristic slip velocities and the forces dominating the flow dynamics such as the ratio of inertial over viscous forces (Reynolds number) and inertial over surface tension forces (Weber number) in the plume regions for both simulations (single and dual plug gas injections) with the aim of providing with further description of the fluid flow behavior of the ladle.Weber numbers are below 10 with both the injections, while average Reynolds numbers are around 1500 and relative velocities of around 0.15 m/s are found in this regions.
analysis.Figure 8; Figure 9 shows characteristic slip velocities and the forces dominating the flow dynamics such as the ratio of inertial over viscous forces (Reynolds number) and inertial over surface tension forces (Weber number) in the plume regions for both simulations (single and dual plug gas injections) with the aim of providing with further description of the fluid flow behavior of the ladle.Weber numbers are below 10 with both the injections, while average Reynolds numbers are around 1500 and relative velocities of around 0.15 m/s are found in this regions.dynamics such as the ratio of inertial over viscous forces (Reynolds number) and inertial over surface tension forces (Weber number) in the plume regions for both simulations (single and dual plug gas injections) with the aim of providing with further description of the fluid flow behavior of the ladle.Weber numbers are below 10 with both the injections, while average Reynolds numbers are around 1500 and relative velocities of around 0.15 m/s are found in this regions.
(a) (b)   The results depicted in Figures 8 and 9 suggest that the numerical procedure proposed in this study is able to provide detailed information on the dominating forces and relative motion of the phases.The results depicted in Figures 8 and 9 suggest that the numerical procedure proposed in this study is able to provide detailed information on the dominating forces and relative motion of the phases.

Conclusions
A comprehensive mathematical model has been developed with an improved methodology to study multiphase flow applying two algorithms, IPSA and VOF using the CFD code PHOENICS.This unconventional approach to represent a two-immiscible liquid phase system stirred by gas injection through single and dual plugs located at the bottom of the ladle was developed by a hybrid combination of a Eulerian representation of a two-phase (dispersed gas bubbles into a continuous single liquid) and a VOF method to separate the liquid into two immiscible liquids (steel and slag).A comparison of experimental data (both centric and off-centric) from a 1/17th water-oil stirred ladle model, involving the values of slag eye area and turbulent kinetic energy, show a reasonable agreement and it is considered that the model has been successfully validated.The application of the developed model in an industrial size ladle indicates the potential to reach a full representation of the flow structure in both liquids, steel and slag.The model provides information on the geometry of the plume, velocity fields and turbulence parameters.

Figure 1 .
Figure 1.(a) Schematic view of the gas stirred ladle with nozzle (both single and dual plug gas injections), (b) Cylindrical grid coordinates of the ladle (axisymmetric section) used in the simulations.

Figure 1 .
Figure 1.(a) Schematic view of the gas stirred ladle with nozzle (both single and dual plug gas injections), (b) Cylindrical grid coordinates of the ladle (axisymmetric section) used in the simulations.

Figure 2 .
Figure 2. (a) Physical model with a scale factor of 1/17th to observe steel-slag behavior using waterengine oil respectively (b) Calculated ladle eye with mathematical model.

Figure 2 .
Figure 2. (a) Physical model with a scale factor of 1/17th to observe steel-slag behavior using water-engine oil respectively (b) Calculated ladle eye with mathematical model.

Figure 3 .
Figure 3. Model validation using turbulent kinetic energy data obtained by (a) particle image velocimetry (PIV) and (b) numerical simulation.The computational fluid dynamics (CFD) model results show a good prediction of experimental results.

Figure 3 .
Figure 3. Model validation using turbulent kinetic energy data obtained by (a) particle image velocimetry (PIV) and (b) numerical simulation.The computational fluid dynamics (CFD) model results show a good prediction of experimental results.

Figure 4 .
Figure 4. Steel-slag distribution (Steel, red phase; blue, slag phase) according to the mathematical model predictions considering the two cases in this study (a) single plug gas injection with nozzle at 1/3 R and (b) dual plug gas injection with nozzles at 2/3 R, observing a small numerical diffusion due to the use of a finer grid in the said interface, improving the degree of precision.

Figure 4 .Figure 5 .
Figure 4. Steel-slag distribution (Steel, red phase; blue, slag phase) according to the mathematical model predictions considering the two cases in this study (a) single plug gas injection with nozzle at 1/3 R and (b) dual plug gas injection with nozzles at 2/3 R, observing a small numerical diffusion due to the use of a finer grid in the said interface, improving the degree of precision.Processes 2019, 7, x FOR PEER REVIEW 13 of 17

Figure 5 .
Figure 5. Velocity fields are shown in the plane of symmetry of the ladle furnace (a) single plug gas injection, the liquid metal shows two recirculation on the sides of the single nozzle, generated by the drag of the bubbles on the liquid steel, (b) dual plug gas injection creates two plumes with two circulating loops on each side of each plume.The flow patterns in the slag layer are different for both injections (single or dual).

Figure 5 .Figure 6 .
Figure 5. Velocity fields are shown in the plane of symmetry of the ladle furnace (a) single plug gas injection, the liquid metal shows two recirculation on the sides of the single nozzle, generated by the drag of the bubbles on the liquid steel, (b) dual plug gas injection creates two plumes with two circulating loops on each side of each plume.The flow patterns in the slag layer are different for both injections (single or dual).

Figure 6 .
Figure 6.Details of 3D velocity patterns on the top slag layer (a) Single plug gas injection (i.e., one nozzle) located at 1/2 R; (b) dual plug as injection (two nozzles) located at 2/3 R. In both the cases, a complex slag flow pattern is observed.

Figure 7 .
Figure 7. Gas plume that forms in the mathematical simulation of the ladle furnace adopted a cylindrical shape and shows a greater fraction of gas in the central axis of the plume and has a single plume for (a) single plug gas injection and two plumes for (b) dual plug gas injection.

Figure 8 .
Figure 8. Single plug gas injection (one nozzle) simulation: Gas-bubble characteristics in the system: (a) Weber number, (b) Reynolds number and (c) Slip velocity in m/s.

Figure 7 .
Figure 7. Gas plume that forms in the mathematical simulation of the ladle furnace adopted a cylindrical shape and shows a greater fraction of gas in the central axis of the plume and has a single plume for (a) single plug gas injection and two plumes for (b) dual plug gas injection.

Figure 7 .
Figure 7. Gas plume that forms in the mathematical simulation of the ladle furnace adopted a cylindrical shape and shows a greater fraction of gas in the central axis of the plume and has a single plume for (a) single plug gas injection and two plumes for (b) dual plug gas injection.

Figure 8 .
Figure 8. Single plug gas injection (one nozzle) simulation: Gas-bubble characteristics in the system: (a) Weber number, (b) Reynolds number and (c) Slip velocity in m/s.

Figure 8 . 17 Figure 9 .
Figure 8. Single plug gas injection (one nozzle) simulation: Gas-bubble characteristics in the system: (a) Weber number, (b) Reynolds number and (c) Slip velocity in m/s.Processes 2019, 7, x FOR PEER REVIEW 15 of 17

Figure 9 .
Figure 9. Dual plug gas injection (two nozzles) simulation: Gas-bubble characteristics in the system: (a) Weber number, (b) Reynolds number and (c) Slip velocity in m/s.

Table 1 .
Boundary conditions of the mathematical model.

Table 2 .
Geometrical and physical variables in the physical model.

Table 3 .
Geometrical and physical variables needed for the simulations.

Table 3 .
Geometrical and physical variables needed for the simulations.