*2.2. Numerical Simulation Models*

Figure 2 shows the schematic of the numerical simulation model included in the paper. The numerical simulation model had the following part: the fluid flow model, the bubble and alloy transport model, the bubble and fluid heat exchange model, the alloy melting model, the VOF model, and the alloy concentration diffusion model.

**Figure 2.** Schematic of the simulation model in the argon-stirred ladle.

#### 2.2.1. Fluid Flow and VOF Model

In this paper, the fluid flow was driven by bubble transport. The governing equations are shown in Equations (1)–(7), and a source term was added to simulate the flow driven by the bubble. The drag force was considered in this paper as the bubble-driven flow, and the source term was added. In this paper, three kinds of fluid phases were considered: molten steel phase, molten slag phase, and argon bubble phase. The continuum surface stress and surface tension models were used in the VOF model. The volume fraction of the bubble in the fluid domain was adjusted by the User-Defined Function (UDF) module. The volume fraction of the bubble was updated according to the bubble's position.

$$\frac{1}{\rho\_q} \frac{\partial}{\partial t} (\alpha\_\eta \rho\_\eta) + \nabla \cdot (\alpha\_\eta \rho\_\eta \stackrel{\rightarrow}{\vec{\nu}}) = 0 \tag{1}$$

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

$$\frac{\partial}{\partial t}(\rho \stackrel{\rightarrow}{\boldsymbol{\upsilon}}) + \nabla(\rho \stackrel{\rightarrow}{\boldsymbol{\upsilon}} \stackrel{\rightarrow}{\boldsymbol{\upsilon}}) = -\nabla p + \nabla(\mu(\nabla \stackrel{\rightarrow}{\boldsymbol{\upsilon}} + \nabla \stackrel{\rightarrow}{\boldsymbol{\upsilon}}^{\boldsymbol{\Gamma}}) - \frac{2}{3}\nabla \stackrel{\rightarrow}{\boldsymbol{\upsilon}}I) + \rho \stackrel{\rightarrow}{\boldsymbol{\xi}} + \stackrel{\rightarrow}{F} \tag{3}$$

$$\frac{\partial}{\partial t}(\rho k) + \frac{\partial}{\partial x\_i}(\rho k u\_i) = \frac{\partial}{\partial x\_j}(a\_k \mu\_{eff} \frac{\partial k}{\partial x\_j}) + G\_k + G\_b - \rho \varepsilon \tag{4}$$

$$\frac{\partial}{\partial t}(\rho \varepsilon) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho \varepsilon u\_i) = \frac{\partial}{\partial \mathbf{x}\_j}(\mathbf{a}\_\varepsilon \mu\_{\varepsilon f} \frac{\partial \varepsilon}{\partial \mathbf{x}\_j}) + \mathbb{C}\_{1\varepsilon} \frac{\varepsilon}{k} (\mathbf{G}\_k + \mathbb{C}\_{2\varepsilon} \mathbf{G}\_b) + \mathbb{C}\_{2\varepsilon} \rho \frac{\varepsilon^2}{k} - \mathbb{R}\_{\varepsilon} + \mathbb{S}\_{\varepsilon} \tag{5}$$

$$\stackrel{\rightarrow}{F} = \frac{\alpha\_l}{\Delta V\_{\text{cell}}} \sum\_{i}^{N} (\stackrel{\rightarrow}{F}\_{D,i}) Q \triangle t \tag{6}$$

$$
\stackrel{\rightarrow}{F}\_{D,i} = \frac{18\mu \mathcal{C}\_{D,i} \text{Re}\_{\vec{l}}}{24d\_{b,\vec{l}}^2 \rho\_{b,\vec{i}}} (\stackrel{\rightarrow}{\tilde{u}}\_{l} - \stackrel{\rightarrow}{\tilde{u}}\_{b,\vec{i}}) \tag{7}
$$

where *α*<sup>q</sup> represents the volume fraction of the *q*th phase, *ρ<sup>q</sup>* represents the density in the qth phase, *v* represents the velocity, n represents the phase count in the simulation, *<sup>ρ</sup>* represents the density, *<sup>μ</sup>* represents the molecular viscosity, *<sup>I</sup>* is the unit tensor, <sup>→</sup> *F* is the source term of the bubble-driven flow, *k* represents the turbulent kinetic energy, *Gk* represents the generation of turbulence kinetic energy due to the mean velocity gradient, *Gb* is the generation of turbulence kinetic energy due to buoyancy, *ε* is turbulent kinetic energy dissipation, *α*<sup>l</sup> represents the volume fraction of the liquid phase, Δ*V*cell represents the volume of the cell, *N* represents the bubble count in the cell, *Q* represents the flow rate of the argon, Δ*t* represents the time step, *CD,i* represents the coefficient of drag force, Rei represents the Reynolds number of bubble, *db,i* represents the diameter of the bubble, *ρb,i* represents the density of the bubble, *u*<sup>l</sup> is the velocity of the fluid, and *ub,i* is the velocity of the bubble.

## 2.2.2. Bubble and Alloy Transport Model

The discrete phase model (DPM) was applied in the bubble and inclusion transport modeling. The bubble and alloy transport equations are shown in Equations (8)–(10). The discrete random walk model was applied in the bubble transport simulation. The randomness of the inclusion transport will increase if the *CL* is increased. The bubble will be removed when the bubble reaches the steel–slag interface. The steel–slag interface is when the volume fraction of slag is 0.5. The bubble diameter will increase if the bubble floats. The bubble diameter variation mechanism is shown in Equation (11).

$$\frac{d\overrightarrow{u}\_p}{dt} = \frac{\overrightarrow{u}\_I - \overrightarrow{u}\_p}{\tau\_r} + \frac{\overrightarrow{g}(\rho\_p - \rho\_\mathbf{f})}{\rho\_p} \tag{8}$$

$$T\_E = -T\_L \ln(r) \tag{9}$$

$$T\_L = \mathcal{C}\_L \frac{k}{\varepsilon} \tag{10}$$

$$d\_{b,i} = \sqrt[3]{\frac{p\_0 + \rho\_1 g h}{p\_0 + \rho\_1 g (H - y\_i)}} d\_{\text{bottom}} \tag{11}$$

where <sup>→</sup> *up* and <sup>→</sup> *ul* are, respectively, the particle vector and the fluid vector; *ρ<sup>p</sup>* is the particle density; *ρ*<sup>f</sup> is the fluid density; *T*<sup>E</sup> is the eddy lifetime; *r* is the random Gaussian number; *TL* is the integral time scale; *CL* (integral time-scale constant) is a constant.

#### 2.2.3. Energy Conservation and Bubble Heat Exchange Model

The governing equation of the temperature is shown in Equations (12) and (13). In this paper, a heat transport model of molten steel and the bubble was applied. The equations of the bubble heat transport are shown in Equations (14)–(19). The bubble temperature is lower than the molten steel temperature, and the heat will transfer from the molten steel to the bubble. The heat transfer energy is based on the energy conservation law. The bubble will stop heat transfer if the bubble temperature is the same as that of the molten steel. The effect of bubble initial injection temperature on the molten steel temperature simulation is discussed in Section 3.2.2.

$$\frac{\partial}{\partial t}(\rho E) + \nabla(\stackrel{\rightarrow}{\upsilon}(\rho E + p)) = \nabla(k\_{eff}\nabla T - \sum\_{j} h\_{j}\stackrel{\rightarrow}{l}\_{j} + (\overline{\mathbb{F}}\_{eff} \cdot \stackrel{\rightarrow}{\upsilon})) \tag{12}$$

$$E = h - \frac{p}{\rho} + \frac{v^2}{2} \tag{13}$$

$$E\_T = C\_{b,i} \cdot m\_{b,i} \cdot (T\_{\text{cell}} - T\_{b,init}) \tag{14}$$

$$E\_{step} = Nu \cdot k\_m \cdot \pi \cdot d\_{b,i} \cdot \left(T\_{cell} - T\_{b,pre\,temp} \right) \cdot \triangle t \tag{15}$$

$$N\mu = 2 + 0.6 \cdot \text{Re}^{1/2} \cdot \text{Pr}^{1/3} \tag{16}$$

$$T\_{\rm b,cur.temp} = T\_{\rm b,pr.temp} + \frac{E\_{step}}{(\mathbb{C}\_{b,i} \cdot m\_{b,i})} \tag{17}$$

$$\text{Pr} = \frac{\mu\_f \cdot \mathbb{C}\_p}{k\_m} \tag{18}$$

$$\text{Re}\_{\mathbb{P}} = \frac{\rho\_{b,i} d\_{b,i} (\stackrel{\rightarrow}{\tilde{u}\_I} - \stackrel{\rightarrow}{\tilde{u}\_{b,i}})}{\mu\_l} \tag{19}$$

where *h* is sensible enthalpy, *ET* is the total energy of the bubble, *Cb.i* is the heat capacity of the bubble, *mb.i* is the bubble mass, *T*cell is the cell's temperature where the bubble is injected into the fluid domain, *Tb.init* is the initial bubble temperature, *Estep* is the energy exchange in each timestep, *Nu* is the Nusselt number, *k*m is the fluid thermal conductivity, *Tb, pre.temp* is the bubble temperature in the previous timestep, and *Cp* is the heat capacity of the fluid.

#### 2.2.4. Alloy Melting and Alloy Concentration Diffusion Model

In this paper, a solid crust will form when the alloy is added to the high-temperature molten steel. At this time, the alloy began melting inside the solid crust. The alloy concentration was released until the solid crust melted at a specific time. The alloy concentration numerical simulation governing equation is shown in Equations (20)–(24). The source term was added in a cell where the alloy melted.

$$t\_{\rm melft} = \frac{\rho\_A C\_{P,A} d\_A}{2\pi h\_\odot} \frac{T\_s - T\_0}{T\_M - T\_s} \tag{20}$$

$$Nu = \frac{d\_A h\_\odot}{k\_m} \tag{21}$$

$$Nu = 2 + \left(0.4 \text{Re}\_A{}^{1/2} + 0.06 \text{Re}\_A{}^{2/3}\right) \cdot \text{Pr}^{2/5} \tag{22}$$

$$\frac{\partial \rho \mathbb{C}}{\partial t} + \frac{\partial}{\partial x\_i} (\rho u\_i \mathbb{C} - \Gamma\_k \frac{\partial \mathbb{C}}{\partial x\_i}) = S\_\mathbb{C} \tag{23}$$

$$S\_{\subset} = \frac{m\_{alloy}}{m\_{\text{cell}}} \tag{24}$$

where *tmelt* is the alloy solid crust melting time, *ρ<sup>A</sup>* is the alloy density, *CP,A* is the heat capacity of the alloy, *dA* is the diameter of the alloy, *hc* is the coefficient of the heat transfer, *Ts* is the solidification temperature of molten steel, *TM* is the temperature of molten steel, T0 is the initial alloy temperature, *C* is the mass concentration of alloy, *Sc* is the source term of the alloy melting, *m*alloy is the mass of the melted alloy, and *m*cell is the mass of the cell where the alloy is melting.

#### *2.3. Boundary Condition and Solution Strategy*

Table 2 shows the relevant parameters that are required in the numerical simulation. The initial bubble diameter and the initial bubble velocity were calculated by Equations (23) and (24). The numerical simulation was calculated by ANSYS Fluent software (ANSYS, Inc., Canonsburg, PA, USA). The inlet boundary condition was set as follows: the bubbles were injected from the plugs, and the inject position was defined by the User-Defined Function (UDF) that was built in the ANSYS Fluent software (ANSYS, Inc., Canonsburg, PA, USA). The boundary name, position, and condition are shown in Figure 3. The bubble flow rate was 5 L·min−1. The bubble injection velocity was 2.63 m·s−<sup>1</sup> in this paper. The initial bubble diameter was 3.924 mm. The position of the alloys thrown was above the slag eye by 2.5 m; based on the law of free fall, the velocity of the alloy contacting with the molten steel was 7 m·s−1. Figure <sup>3</sup> shows the hexahedral mesh of the fluid zone in the prototype and hydraulic model. The argon bubble was injected from the bottom of the ladle, and the injection position is shown in Figure 3. In the prototype, argon was injected at two places, and the angle between them was 135◦. The solution initialization used standard initialization and computing from all zones. The solution step was 0.001 s, the injection rate of the bubble was one bubble per solution step, and the *Nper* was 1000. The free surface was used on the top of the slag (in prototype modeling) or the top of the air (in the hydraulic modeling). The rest boundary was taken as the wall boundary. The SIMPLE scheme was used in the pressure-velocity coupling. The residual in the governing equations was 0.001. The heat flux in the ladle wall, the molten slag surf, and the ladle bottom was 0.006 W2·mm<sup>−</sup>2.

$$d\_{b,init} = \left(0.35 \cdot \left(\frac{Q^2}{9.8}\right)^{0.2}\right) \tag{25}$$

$$v\_{b,init} = \frac{Q}{\left(N\_{per} \cdot \frac{1}{6} \cdot \pi \cdot d\_{b,init}^3\right)}\tag{26}$$

where *vb,init* is the initial velocity when the bubble is injected into the fluid zone and *Nper* is the bubble injection count in 1 s.


**Table 2.** Relevant parameters in the numerical simulation.

**Figure 3.** Hexahedral mesh of the fluid zone, boundary name, boundary condition, and the schematic of the argon injection position: (**a**) prototype; (**b**) hydraulic model.

#### **3. Results and Discussion**

*3.1. Model Validation*

3.1.1. Effect of Discrete Random walk Model on the Simulation Results

Figure 4 shows the simulation results of the water model simulation and the PIV measurements results of the water model. Figure 4a shows the slice position of the numerical simulation results. The PIV measurements area is identical to the slice position. Figure 4b shows the velocity contour of the slice position, the velocity is higher in the center, and the bubbles are transported in the center, which accelerates the fluid velocity when the bubble random transport model is deactivated. The bubble is upward in a straight line, and the velocity contour is symmetrical. Figure 4c shows the velocity contour when the bubble random transport model is activated, and the *CL* is 0.15. The area of the velocity at a higher speed is in a reverse cone shape. The maximum velocity is lower than when the random walk model is deactivated. The bubbles are transported in a random way; the degree of dispersion of the bubble position is higher, which results in a lower maximum velocity. Figure 4d shows the PIV measurement result, and the measurement area is plotted in Figure 4c in a red scatter line. Figure 4d shows that the velocity is higher at the top of the ladle model, the velocity in a higher-speed area is in a reverse cone shape, and the velocity contour distribution of the numerical simulation is similar to the PIV measurement results.

**Figure 4.** (**a**) The slice position of the simulation results of the water model, (**b**) the velocity magnitude contour in the slice when the random walk model was deactivated, (**c**) the velocity magnitude contour in the slice when the random walk model was activated and the *CL* was 0.15, and (**d**) the time-averaged results of the PIV measurements in the water model.

Figure 5 shows the quantitative data of the velocity magnitude of the numerical simulation and PIV measurement results. Figure 5a–c show the simulation results of the bubble random walk model being deactivated and the PIV measurement results. The results show that the velocity magnitude distribution is sharp, and the velocity peak is narrow and higher compared to the PIV results. The flow of the fluid in the water model is driven by the bubbles. When the bubble random walk model is deactivated, the flow of the bubbles in the fluid presents a vertical upward movement, so the velocity peaks appear high and narrow. Figure 5d–f show the simulation results when the bubble random walk model is activated, the *CL* is 0.15, and the PIV measurements results. The velocity magnitude peak is smaller and broader than when the bubble random walk model is deactivated.

**Figure 5.** The quantitative velocity magnitude data of numerical simulation of water model and PIV measurements results, when the random walk model was deactivated (**a**) near the water surface, (**b**) under the water surface by 100 mm, and (**c**) under the water surface by 200 mm; when the random walk model was activated and the *CL* was 0.15 (**d**) near the water surface, (**e**) under the water surface by 100 mm, and (**f**) under the water surface by 200 mm.

#### 3.1.2. Effect of Time Scale Factor on the Simulation Results

Figure 6 shows the velocity magnitude of the numerical simulation results and the PIV measurements. The numerical simulation result under the case of *CL* is 0.3 or 0.45. Figure 7 shows the variance in the velocity difference between the numerical simulation and the PIV measurement results. The calculation range is the velocity data from −100 mm to 100 mm.

In Figure 7, *DL* represents the variance in the velocity magnitude difference between the numerical simulation and the PIV measurements, where a *CL* of 0 illustrates the case where the bubble random transport model is deactivated. The results show that the variance is decreased when the *CL* increases from 0 to 0.3, then the variance is increased when the *CL* increases from 0.3 to 0.45. The *CL* has a significant influence on the surf velocity. The randomness of the bubble transport is higher at the top of the ladle, then the variation in *CL* has a significant impact on the surf velocity. The velocity difference is less when the *CL* is 0.3.

**Figure 6.** The quantitative velocity magnitude data of numerical simulation of water model and PIV measurements results, when the random walk model was activated and *CL* was 0.3 (**a**) near the water surface, (**b**) under the water surface by 100 mm, and (**c**) under the water surface by 200 mm; when the random walk model was activated and the *CL* was 0.45 (**d**) near the water surface, (**e**) under the water surface by 100 mm, and (**f**) under the water surface by 200 mm.

**Figure 7.** The variance in the velocity difference between PIV measurements and numerical simulation at −100 mm to 100 mm in Figures 5 and 6.
