**2. Methodology: Numerical Model Development**

A Eulerian multi-phase mathematical model simulating the physical model described in this study was developed under the following assumptions: (i) steady state, (ii) a symmetry plane is considered for the dual gas injection system and thus only half of the ladle is solved; (iii) Newtonian and incompressible fluids for both liquids and gas phases; (iv) isothermal flow and (v) constant bubble diameter of 0.01 m. The latter is certainly an oversimplification, by neglecting bubble disintegration and coalescence phenomena under real dynamic bubble size distributions. However, using the Eulerian model most regions of the ladle are fairly well predicted, as presented in Section 3.1, both in turbulence and velocity magnitude, except for the water-oil-bubble open eye regions. Governing equations for the 3D Eulerian-Eulerian multi-phase algorithm include mass conservation equation, Navier–Stokes equations and the *k*-ε realizable turbulence model for the water phase.

The volume of the *q*-phase, *Vq*, is given by the volume integral:

$$V\_{\eta} = \int\_{V} \kappa\_{\eta} \, dV \tag{1}$$

where ∝*<sup>q</sup>* is the volume fraction of phase *q*, and the sum of volume fractions must be equal to one according to:

$$\sum\_{q=1}^{n} \kappa\_q = 1\tag{2}$$

The continuity equation for each *q*-phase is:

$$\nabla \cdot \left( \ltimes\_q \rho\_q \overrightarrow{\upsilon}\_q \right) = 0 \tag{3}$$

where ρ*<sup>q</sup>* and <sup>→</sup> *v q* are the density and velocity vector of the *q*-phase, respectively.

The momentum conservation equation for the *q*-phase is:

$$\nabla \cdot \left( \propto\_q \rho\_q \overrightarrow{\boldsymbol{\upsilon}}\_q \overrightarrow{\boldsymbol{\upsilon}}\_q \right) = -\propto\_q \nabla P + \nabla \cdot \left( \propto\_q \mu\_{ef,q} \left( \overrightarrow{\boldsymbol{\upsilon}} \overrightarrow{\boldsymbol{\upsilon}}\_q + \left( \overrightarrow{\boldsymbol{\upsilon}} \overrightarrow{\boldsymbol{\upsilon}}\_q \right)^T \right) \right) + \propto\_q \rho\_q \overrightarrow{\boldsymbol{\chi}}\_q + \overrightarrow{F}\_T \tag{4}$$

For the water phase:

$$
\mu\_{cf,l} = \mu\_{lam,l} + \mu\_{t,l} \tag{5}
$$

For other phases:

$$
\mu\_{ef,q} = \mu\_{lam,q} \tag{6}
$$

where *P*, μ*e f*,*q*, → *g* are the pressure, effective viscosity of the *q*-phase and the gravity acceleration, respectively. The effective viscosity for the water phase is the sum of the molecular viscosity (μ*lam*,*l*) and the turbulent viscosity (μ*t*,*l*) defined by the turbulence model employed. For the other phases μ*lam*,*<sup>q</sup>* is only the laminar viscosity of every fluid, and the subscripts *l* and *lam* stand for water and laminar, respectively. The term <sup>→</sup> *FT* is the contribution of all interphase forces. The only force considered is the drag force, and the virtual mass, lift and turbulent dispersion forces are neglected.

In this study, a sensitivity analysis was performed to choose the best turbulence model that predicts correctly the turbulence measured in the ladle. From all the models tested, the best option is the realizable *k*-ε turbulence model [21]. This model solves two additional conservation equations applicable only to the water liquid phase, one of these equations being the conservation of turbulent kinetic energy, *k*:

$$\nabla \cdot \left( a\_l \rho\_l k \overrightarrow{\boldsymbol{\sigma}}\_l \right) = \nabla \cdot \left( a\_l \left( \mu\_{l \text{m}, l} + \frac{\mu\_{l, l}}{\sigma\_k} \right) \nabla k \right) + a\_l \mathbf{G}\_{k, l} + a\_l \mathbf{G}\_b - a\_l \rho\_l \boldsymbol{\varepsilon} + a\_l \rho\_l \pi\_{k, l} \tag{7}$$

The conservation equation for the dissipation of the turbulent kinetic energy, ε:

$$\nabla \cdot \left( \alpha\_l \eta\_l \varepsilon \stackrel{\scriptstyle \rightarrow}{\upsilon}\_l \right) = \nabla \cdot \left( \alpha\_l \left( \mu\_{lm,l} + \frac{\mu\_{l,l}}{\sigma\_\varepsilon} \right) \nabla \varepsilon \right) + \alpha\_l \mathbb{C}\_1 \left( \frac{\varepsilon}{k} \mathbb{G}\_b \mathbb{C}\_3 \right) - \alpha\_l \eta\_l \mathbb{C}\_2 \frac{\varepsilon^2}{k + \sqrt{\frac{\varepsilon \mu\_{lm,l}}{\rho\_l}}} + \alpha\_l \eta\_l \pi\_{\varepsilon,l} \tag{8}$$

*C*1, *C*2, *C*3, σ*<sup>k</sup>* and σε are constants of the model, with values of 1.44, 1.92, 1.3, 1.0 and 1.2 respectively. *Gk* is the production of turbulent kinetic energy due to the mean velocity gradients of the water phase and *Gb* is the additional turbulent kinetic energy produced by the bubbles. π*k*,*<sup>l</sup>* and πε,*<sup>l</sup>* are the turbulent interaction terms defined by Troshko-Hassan [22] as:

$$\pi\_{k,l} = \mathbb{C}\_{kc} \sum\_{p=1}^{m} K\_{\mathcal{J}} \left| \overrightarrow{\boldsymbol{\upsilon}}\_{\mathcal{J}} - \overrightarrow{\boldsymbol{\upsilon}}\_{l} \right|^{2} \tag{9}$$

$$
\pi\_{\varepsilon,l} = \mathbb{C}\_{td} \frac{\mathfrak{NC}\_A \Big| \overrightarrow{\boldsymbol{\upsilon}}\_{\mathcal{S}} - \overrightarrow{\boldsymbol{\upsilon}}\_I}{2d\_{\mathcal{S}}} \pi\_{k,l} \tag{10}
$$

where *Cke* and *Ctd* are model constants with values of 0.75 and 0.45, respectively. *Kgl* is the covariance of the velocities of the continuous phase *l* and the disperse phase *g*. *dg* is the bubble diameter; the velocity difference between the gas and liquid phase, <sup>→</sup> *<sup>v</sup> <sup>g</sup>* <sup>−</sup> <sup>→</sup> *v l* , is defined as the relative velocity, <sup>→</sup> *v rel*; whereas *CA* is the drag coefficient determined through a sensitivity analysis (not mentioned in this study) performed that revealed the best option was the well-known symmetric model:

$$C\_A = \left\{ \begin{array}{c} \frac{24 \left(1 + 0.15 Re^{0.687} \right)}{Re} \text{ } Re \le 1000 \\ 0.44 \text{ } Re > 1000 \end{array} \right. \tag{11}$$

where *Re* is the Reynolds number.

The interphase forces, <sup>→</sup> *FT* are limited to the drag force, <sup>→</sup> *FA* as follows:

$$
\stackrel{\rightharpoonup}{F}\_T = \stackrel{\rightharpoonup}{F}\_A = \frac{3\alpha\_\mathcal{S}\alpha\_l\rho\_l\mathcal{C}\_A}{4d\_\mathcal{\mathcal{S}}} (\stackrel{\rightharpoonup}{\upsilon\_\mathcal{S}} - \stackrel{\rightharpoonup}{\upsilon\_l})\tag{12}
$$

Finally, to compute mixing time in the water phase *l*, a single conservation equation for a tracer chemical species *i* is solved in transient mode with initial conditions of zero concentration of solute everywhere except for the pulsed tracer injected at the free surface above the plume of high flow, as follows:

$$\frac{\partial}{\partial t}(\rho\_l w\_{i,l}) + \nabla \cdot (\rho\_l \overrightarrow{u\_l} w\_{i,l}) = \nabla \cdot \left(\rho\_l D\_{i,l} + \frac{\mu\_{l,l}}{S c\_l}\right) \nabla w\_{i,l} \tag{13}$$

The geometry of the water (physical) model that was used in PIV experiments is made in 3D using ANSYS Design Modeler 19.0 (Ansys Inc., Canonsburg, PA, USA) for the numerical model. The diameter and height of the cylindrical ladle are 0.185 m and 0.214 m (0.17 m liquid level), respectively, similar to that of the physical model [15] that corresponds to a 1:17 scale ratio of an industrial-size ladle of 200-ton capacity. Figure 1 shows an illustration of the mesh built with approximately 350,000 elements with average orthogonality, skewness and aspect ratios of 0.98, 0.1 and 1.99, respectively. The choice of the mesh elements was based on a sensitivity analysis using approximately 200,000 and 500,000 mesh elements.

**Figure 1.** Computational domain of the ladle used for the numerical simulations presented in this study. The two vertical sections indicating the gas injection inlets and the top horizontal section indicating the slag layer are comparatively denser than the remaining mesh domain. The inset shows the mesh density for the slag zone which is different from the melt zone. The top portion of the mesh (half along the symmetry plane) is shown separately.

Non-slip boundary conditions at the bottom and lateral walls have been used, whereas the standard wall functions have been used to connect the turbulent core of the fluid with the laminar flow near the walls. The gas injection inlets at the nozzle positions and an open boundary to the atmosphere at the top of the system allow the outflow of the gas phase. The complete set of boundary conditions can be found in Table 1.


**Table 1.** Boundary conditions used in the numerical model presented in this study.

The system of partial differential equations was numerically solved and the solution in pseudo-transient mode was considered to be converged when the residuals of all conservation equations were below 1 <sup>×</sup> <sup>10</sup><sup>−</sup>3. Approximately 1200 iterations were required to get the convergence in around 36 h of CPU time in a computer with 8 MB in RAM with an Intel Core® i7-3770 processor of 3.4 GHZ. Table 2 lists all the numerical simulations performed in this study which is based on a full factorial experimental design at two levels with the three variables, namely, gas flow rate, dual gas injection ratio and (slag) oil thickness, as mentioned in Jardón-Pérez et al. [15].

**Table 2.** Design of experiment with high and low values of the three variables, namely, gas flow rate, dual gas injection ratio and (slag) oil thickness for the eight case studies presented in this study.

