Numerical Simulation of Bubble Size Distribution in Single Snorkel Furnace (SSF) with Population Balance Model (PBM)

Abstract

Single Snorkel Furnace (SSF) vacuum refining furnace is a novel external refining equipment for high clean steel production. RH is a molten steel refining technology developed by Rheinstahl-Heraeus company. Compared with the traditional RH furnace, the SSF furnace has the advantages of a simple structure, high refining efficiency, and low production cost. However, because the upward flow and the downward flow are in a single snorkel, the flow phenomenon is more complex than that in the RH device. Therefore, the gas–liquid two-phase flow law in SSF furnaces plays an important role in improving equipment efficiency and accurate control. In addition, the evolution and movement behavior of bubbles have an important influence on the two-phase flow. In this study, the Population Balance Model (PBM) model is employed to study the bubble properties, taking into account the effect of bubble coalescence and breakup on the flow field. The simulation results with this model are consistent with the experimental values, and the comparison with the results of the model without the PBM is revealed to be closer with less error. The results show that with the PBM model the flow field is more homogeneously distributed, the flow velocity is more stable, and the area distribution of the upward flow and downward flow in the snorkel is more symmetrical. In the case of this study, as the fluid level rises, the bubble diameter will increase due to the decrease in hydrostatic pressure.

1. Introduction

In 2021, China’s crude steel output [1] fell by 2.8% to 1033, million tons, while steel production increased by 0.9% to 1336 million tons. Low-carbon steel is one of the widely adopted types of steel that has had the carbon and impurity content reduced through metallurgy processes. Its structure and mechanical properties are superior to those of regular steel, making it a popular choice for national defense, aviation, and space applications, as well as other industries where stringent steel quality is required. Applications of low-carbon steel are expected to develop further, and its demand will soar in the near future. Decarburization is a crucial process in the production of low-carbon steel where a traditional Rheinstahl–Heraeus (RH) refining furnace, which was first developed by Rheinstahl and Heraeus company in Germany, is used. Lately, benefiting from its superior refining efficiency compared to RH furnace, the single-snorkel vacuum refining furnace (SSF) has attracted a lot of interest in both industrial and research communities. Figure 1 shows the comparison of the traditional RH and the SSF furnaces. The main structure of both furnaces is similar, consisting of a ladle and a vacuum chamber. The major difference between the two is the arrangement of snorkels. Figure 1a shows that there are two snorkels in the RH furnace, one up-leg snorkel (left) and one down-leg snorkel (right). The fluid enters the vacuum chamber along the up-leg snorkel as the argon gas is injected from the inlet, and sequentially circulates back to the ladle through the down-leg snorkel. On the other hand, there is only one snorkel in the SSF furnace (see in Figure 1b). Argon gas is injected through the eccentric inlet on the bottom of the lade and circulates within the snorkel. The liquid steel is drawn into the ladle due to the pressure difference between the vacuum chamber and ambient air. Powered by buoyancy forces, argon gas goes up to the free surface of the liquid steel from one side of the ladle (right) and then flows back from the opposite side (left), causing the liquid steel to circulate and achieving degassing and decarburization purposes. The SSF furnace has a relatively simple structure and simple operation. More importantly, it makes it easier to increase the area of contact between argon gas and steel; reducing the issue of spattering in the decarburization process and posing a high refining efficiency compared to the conventional RH refining furnace.

1.1. Previous Works in Flow Structures and Decarbonization Modelling

With the advantages of a simple structure and outstanding refining efficiency, the SSF furnace is gradually replacing the conventional RH furnace. Chen et al. [3] compared the SSF furnace’s circulation efficiency to that of the RH furnace during the modeling process. He concluded that the SSF furnace’s circulation efficiency is 15% higher and the associated mixing time is 20% shorter in comparison to the RH furnace. Through analyzing the circulation flow characteristics, Dai et al. [4] concluded that steel in the RH furnace is totally mixed between the vacuum chamber and the ladle, while it is not completely exchanged in the SSRF furnace. Yang et al. [5] and Qin et al. [6] performed a series of parametric studies to optimize the circulation efficiency of an 80 t single-snorkel vacuum refining furnace (SSF). The circulation efficiency was found dependent on the eccentric inlet position, the diameter, and the depth of the snorkel. Zhang et al. [7] and Ouyang et al. [8] traced the salt tracer in the water model and evaluated the effect of tracer concentration on mixing time to figure out the flow behavior in the SSF furnace.

The main purpose of refining furnaces is decarburization, degassing, and desulfurization, so the research on decarburization efficiency has been an important research direction in recent years. Duan et al. investigated the impact of several model parameters (e.g., gas flow rate, ladle slag thickness, vacuum chamber height, and pressure) on the overall decarburization process [9,10]. Chen et al. [11] employed the Volume-Of-Fluid (VOF) coupling with Discrete Particle Model (DPM) to study the multiphase flow structure in the RH furnace. Their results predicted showed that the decarburization effect in the furnace is weakened in the low-pressure environment. In a separate study on the metallurgical process inside the RH furnace, Lei et al. [12] revealed that the carbon content in the molten steel steadily decreases as the vacuum level decreases. Geng et al. [13] found that decarburization on the surface of the Ar bubble accounted for most of the decarburization in the RH furnace model, and the amount of decarburization was about twice as much as that inside the vacuum chamber. Rui and Duan et al. [14,15] constructed an SSF furnace mathematical model and investigated the variation of desulfurization rate with time under the influence of solvent sulfur content, steel stirring energy, and other factors using the Euler–Euler method. Dou et al. [16] investigated the concentration distribution of carbon and oxygen as well as the mixing procedure in a furnace. Their results revealed that the argon flow rate is crucial in enhancing the mixing efficiency; lowering the mixing time and carbon mass distribution.

Chen et al. [17] observed that CO₂ would improve the decarburization efficiency in the early stages compared with the conventional Ar gas injection. Considering the economic and environmental impact, the carbon content is only 2.6 ppm higher in the RH furnace after the injection of CO₂ than that of Ar. Therefore, CO₂ can be considered instead of argon in the decarburization process. In order to investigate flow structures and non-equilibrium decarburization behavior within the furnace, Chen et al. [18,19] considered the stirring effect of CO gas on the flow field in the furnace. The results show that the effect of CO gas on the flow field is not negligible. The above studies revealed that bubble dynamics and their interfacial effects have a significant impact on the mixing flow structures and decarburization process.

1.2. Significance of Bubble Dynamics and Its Impact on Decarburization

Theoretically speaking, the decarburization reaction in the refining furnace mainly occurs at three interfacial locations [20,21,22]: the surface of CO bubbles generated inside the steel; the surface of Ar bubbles; and the free surface of the steel in the vacuum chamber. Among these locations, previous studies in RH furnace have shown the majority of decarburization reactions take place at the interface of Ar bubbles. The decarburization efficiency is significantly affected by the interfacial area of Ar bubbles. Moreover, interfacial momentum exchanges are also sensitive to the size of bubbles. Obviously, compared with small diameter bubbles, large diameter bubbles are driven by higher buoyance force inducing higher flow velocity and better mixing efficiency in the furnace. Large bubbles could also deviate from the spherical shape and bubbly flow will transform into cap-bubbly flow or possibly slug flow [23]. Therefore, studying the bubble size distribution is beneficial for gaining a better understanding of the flow field, mixing efficiency, and the sequential decarbonization process.

Nevertheless, the overall bubble size in the system is governed by complex bubble dynamic processes: bubble coalescence and breakage processes. In the coalescence process, bubbles with different sizes travel at different speeds which might collide and merge together forming bigger diameter bubbles. On the other hand, during the breakage process, larger bubbles could break and separate into two or multiple small bubbles via the shear induced surface instability. In other fields, Wu et al. [24] has used experiments to investigate specifically the velocity and size distribution of bubbles in continuous casting molds. Nevertheless, previous studies generally assumed average bubble size to simulate the dispersed gas–liquid flow, rather than the local bubble size distribution. Qi and Dai et al. [2,25,26] analyzed the features of bubble size, density, and velocity within the SSRF furnace. Their studies confirmed the significant effect of gas expansion on the circulation velocity within the vacuum chamber. Aiming to resolve local bubble size distribution, Chen et al. [27] adopted the Population Balance Model (PBM) to consider the effects of bubble breakup, coalescence, and expansion in the SSRF furnace. Unfortunately, the study did not focus on the bubble dynamic in depth. Limited studies on the bubble dynamics have been carried out and its effect on the decarburization process remains elusive. Aiming to understand the complex bubble dynamics and its role in the decarburization process, this article presents a mathematical framework to resolve the variation of bubble size and its associated interfacial effects using the PBM technique. The numerical prediction is first verified for the mixing time in the furnace, and then it is a detailed explanation of the phenomenon of bubble coalescence and breakage.

2. Mathematical Modeling

2.1. Euler–Euler Two-Fluid Model

The mathematical model developed is based on the two-fluid Euler–Euler method. Regarding both liquid and gas phases as continuous terms, the continuity equations and momentum equations for liquid and gas phases are described as:

$$\frac{\partial \left({\rho }_l{\alpha }_l\right)}{\partial t}+\nabla \cdot \left({\alpha }_l{\rho }_l{\stackrel{\rightharpoonup }{u}}_l\right)=0,$$

(1)

$$\frac{\partial \left({\rho }_g{\alpha }_g\right)}{\partial t}+\nabla \cdot \left({\alpha }_g{\rho }_g{\stackrel{\rightharpoonup }{u}}_g\right)=0,$$

(2)

In the above equation, the subscripts $l$ and $g$ represent the liquid and gas phases, respectively, and $\rho$, α, t, and u represent the density, phase fraction, time, and velocity, respectively.

Momentum equation:

$$\frac{\partial }{\partial t}\left({\alpha }_l{\rho }_l{\stackrel{\rightharpoonup }{u}}_l\right)+\nabla \cdot \left({\alpha }_l{\rho }_l{\stackrel{\rightharpoonup }{u}}_l{\stackrel{\rightharpoonup }{u}}_l\right)=-{\alpha }_l\nabla p+\nabla \cdot {\stackrel{=}{\tau }}_l+{\alpha }_l{\rho }_l\stackrel{\rightharpoonup }{g}+{\stackrel{\rightharpoonup }{F}}_{\lg },$$

(3)

$$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g{\stackrel{\rightharpoonup }{u}}_g\right)+\nabla \cdot \left({\alpha }_g{\rho }_g{\stackrel{\rightharpoonup }{u}}_g{\stackrel{\rightharpoonup }{u}}_g\right)=-{\alpha }_g\nabla p+\nabla \cdot {\stackrel{=}{\tau }}_g+{\alpha }_g{\rho }_g\stackrel{\rightharpoonup }{g}+{\stackrel{\rightharpoonup }{F}}_{gl},$$

(4)

where $\stackrel{=}{\tau }$ denotes the stress-strain tensor, whose expression is given by:

$${\stackrel{=}{\tau }}_l={\alpha }_l{\mu }_l\left(\nabla {\stackrel{\rightharpoonup }{u}}_l+{\left(\nabla {\stackrel{\rightharpoonup }{u}}_l\right)}^T-\frac{2}{3}I\left(\nabla \cdot {\stackrel{\rightharpoonup }{u}}_l\right)\right),$$

(5)

${\mu }_l$ denotes the effective viscosity, and the relationship between ${\mu }_l$ and ${\mu }_g$ can be expressed as:

$${\mu }_g=\frac{{\rho }_g}{{\rho }_l}{\mu }_l,$$

(6)

2.2. Momentum Transfer

Multi-phase flow in the interface momentum transfer its main effect, the momentum Equation (3) ${\stackrel{\rightharpoonup }{F}}_{\mathit{lg}}$ denotes the inter-phase forces, the main forces between the two phases include drag, lift, and virtual mass force, the main relationship is expressed as:

$$F_{\lg }=-F_{gl}=F_{\lg }{}^D+F_{\lg }{}^L+F_{\lg }{}^{VM},$$

(7)

The right-hand side of the equation is represented as drag(D), lift(L), and virtual mass (VM) forces, and the forces are described below.

2.2.1. Drag Force

Drag force between the gas-liquid phases due to relative motion and as a consequence momentum transfer:

$$F_{\lg }{}^D=\frac{1}{8}{C}_D{\rho }_l{A}_i\left|{\stackrel{\rightharpoonup }{u}}_g-{\stackrel{\rightharpoonup }{u}}_l\right|\left({\stackrel{\rightharpoonup }{u}}_g-{\stackrel{\rightharpoonup }{u}}_l\right),$$

(8)

where $C_D$ denotes the drag coefficient and represents the area of the interface between the two phases $A_i$. For this drag coefficient calculation, using the Schiller and Naumann Model, the expression of the drag coefficient $C_D$ is:

$$C_D=\left\{\begin{array}{cc}24\left(1+0.15{\mathrm{Re}}^{0.678}\right)/\mathrm{Re}& \mathrm{Re}\le 1000\\ 0.44& \mathrm{Re}>1000\end{array}\right.,$$

(9)

2.2.2. Lift Force

After the gas bubble enters the liquid from the inlet, the gas is subjected to a lateral force in the liquid due to the radial velocity gradient of the liquid phase, i.e., the lift force, the expression of which is:

$$F_{\lg }{}^L=-C_L{\rho }_l{\alpha }_g\left({\stackrel{\rightharpoonup }{u}}_l-{\stackrel{\rightharpoonup }{u}}_g\right)\times \left(\nabla \times {\stackrel{\rightharpoonup }{u}}_l\right),$$

(10)

where $C_L$ denotes the lift force coefficient. The lift force model we currently employ is the Tomiyama Lift Force Model, which is mainly applicable to the lift force of larger scale deformable bubbles in ellipsoidal and spherical cap regions, and its expression is:

$$C_D=\left\{\begin{array}{cc}\min \left[0.288\tanh \left(0.121{\mathrm{Re}}_p\right),f\left(E{o}^{\prime }\right)\right]& E{o}^{\prime }\le 4\\ f\left(E{o}^{\prime }\right)& 4<E{o}^{\prime }\le 10\\ -0.27& 10<E{o}^{\prime }\end{array}\right.,$$

(11)

where

$$f\left(E{o}^{\prime }\right)=0.00105E{o}^{\prime }{}^3-0.0159E{o}^{\prime }{}^2-0.0204E{o}^{\prime }+0.474,$$

(12)

The model relies on the Eötvös number, whose defining equation is:

$$E{o}^{\prime }=\frac{g\left({\rho }_l-{\rho }_g\right)d_h^2}{\sigma },$$

(13)

$$d_h=d_b{\left(1+0.163E{o}^{0.757}\right)}^{1/3},$$

(14)

$$Eo=\frac{g\left({\rho }_l-{\rho }_g\right)d_b^2}{\sigma },$$

(15)

where $\sigma$ is the surface tension, $d_b$ is the bubble diameter, and $d_h$ is the long axis of the deformable bubble.

2.2.3. Virtual Mass Force

In the study of two-phase flow, the gas phase enters the liquid phase and accelerates relative to the liquid phase, and the inertia of the accelerated bubble liquid phase mass exerts a “virtual mass force” relative to it, expressed as:

$$F_{\lg }{}^{VM}=C_{VM}{\alpha }_g{\rho }_l\left(\frac{d_l{\stackrel{\rightharpoonup }{u}}_l}{dt}-\frac{d_g{\stackrel{\rightharpoonup }{u}}_g}{dt}\right),$$

(16)

where $C_{\mathit{VM}}$ is the virtual mass force coefficient with a constant value of 0.5.

2.3. User-Defined Scalar (UDS) Transport Equations

In this paper, the mixing time of the model is used to verify the correctness of the model by adding a tracer to the free surface of the ladle after the flow field in the furnace is stabilized, in which the UDS diffusion equation is utilized to model the solute transport in the furnace with the equation:

$$\frac{\partial {\alpha }_l{\rho }_l{\varphi }_l}{\partial t}+\nabla \cdot \left({\alpha }_l{\rho }_l{\stackrel{\rightharpoonup }{u}}_l{\varphi }_l-{\alpha }_l{\mathsf{\Gamma }}_l\nabla {\varphi }_l\right)=0,$$

(17)

where ${\varphi }_l$ is the scalar of the UDS equation and ${\mathsf{\Gamma }}_l$ is the set diffusion coefficient. In this paper, in order to calculate the mixing time, the normalized concentration equation in Qi [2] is used:

$$C_n=\frac{C_t}{C_{\infty }},$$

(18)

where $C_n$ is the normalized concentration, $C_t$ is the concentration at the time of t, and $C_{\infty }$ is the steady concentration; and the mixing time was determined as the time when the gaps between the normalized concentration were less than 5%.

2.4. Population Balance Model (PBM) for Bubble Dynamics

In this multiphase flow system, the size of the bubble transport may change with the occurrence of transport and chemical reactions, the evolution of bubble size depends on the coalescence, breakage, and other phenomena, so a balance equation is needed to describe the overall change in bubbles. As mentioned earlier, this paper adopted the PBM model to resolve the local bubble size distribution and its evolution due bubble interactions. The PBM technique uses a discrete method to divide the bubble population into 10 bubble size groups and calculate the bubble size distribution using the number density equation. The connection between the PBM model and the Euler–Euler Two-Fluid model is the connection between the number density (n) and the volume fraction:

$$\alpha \left(\overrightarrow{x},t\right)={\displaystyle \underset{{\mathsf{\Omega }}_{\varphi }}{\int }n}V\left(\varphi \right)d{V}_{\varphi },$$

(19)

where “external coordinates” ($\overrightarrow{x}$) denotes the spatial position of the bubble, “internal coordinates” ($\varphi$) denotes bubble size, and $\mathrm{V}\left(\varphi \right)$ is the volume of a bubble in state.

The transport equation for the number density function is given as:

$$\frac{\partial }{\partial t}\left[n\left(V,t\right)\right]+\nabla \cdot \left[\stackrel{\rightharpoonup }{u}n\left(V,t\right)\right]=B_C+B_B-D_C-D_B,$$

(20)

where n denotes the number density of bubbles, $B_C{,B}_B$ denote the birth of bubbles due to coalescence and breakup, respectively, and $D_C,D_B$ denote the death of bubbles due to coalescence and breakup, respectively.

2.4.1. Breakage

The birth and death of bubbles due to breakage are mainly obtained from the following equation:

$$B_B={\displaystyle \underset{{\mathsf{\Omega }}_v}{\int }\rho g\left(V^{\prime }\right)\beta \left(\left.V\right|V^{\prime }\right)}n\left(V^{\prime }\right)d{V}^{\prime },$$

(21)

$$D_B=g\left(V\right)n\left(V\right),$$

(22)

$\mathrm{g}\left(V^{\prime }\right)$ is the breakup frequency, $V^{\prime }$ denotes the volume. The ratio of fragmented bubbles per unit of time, and $\mathsf{\beta }\left(V∕V^{\prime }\right)$ denotes the probability density function of bubbles fragmented from volume $V^{\prime }$ to volume $\mathrm{V}$.

The breakage rates of bubbles due to turbulence induced breakage is modeled according to the Luo and Lehr Breakage Kernels [28] method.

$${\mathsf{\Omega }}_{br}\left(V,V^{\prime }\right)=K{\displaystyle \underset{{\xi }_{\min }}{\overset{1}{\int }}\frac{{\left(1+\xi \right)}^2}{{\xi }^n}\exp \left(-b{\xi }^m\right)d\xi },$$

(23)

where:

$$K=0.9238{\epsilon }^{1/3}{d}^{-2/3}\alpha ,$$

(24)

$$b=12\left[f^{2/3}+{\left(1-f\right)}^{2/3}-1\right]\sigma {\rho }^{-1}{\epsilon }^{-2/3}{d}^{-5/3}{\beta }^{-1},$$

(25)

where ${\Omega }_{\mathit{br}}$ is the breakup rate, $\xi =\lambda /d_j$ is the size ratio between an eddy and a bubble in the inertial sub-range and consequently ${\xi }_{\min }={\lambda }_{\min }/d_j$, the values of n and m are $\frac{3}{11}$ and $-\frac{11}{3}$, respectively, and the expression of $f$ is $f=V∕V^{\prime }$.

2.4.2. Coalescence

The birth and death of bubbles due to coalescence are described by the following equation:

$$B_C=\frac{1}{2}{\displaystyle \underset{0}{\overset{V}{\int }}a\left(V-V^{\prime },V^{\prime }\right)}n\left(V-V^{\prime }\right)n\left(V^{\prime }\right)d{V}^{\prime },$$

(26)

$$D_C={\displaystyle \underset{0}{\overset{\infty }{\int }}a\left(V,V^{\prime }\right)}n\left(V\right)n\left(V^{\prime }\right)d{V}^{\prime },$$

(27)

where the $\frac{1}{2}$ factor in Equation (26) is added to avoid performing two collision interpretations.

The coalescence rate of bubbles due to turbulence induced coalescence is modeled according to the Prince and Blanch Coalescence Kernel [29] method. It is clear that collisions may result from the random motion of bubbles due to turbulence. The coalescence model assumes that the coalescence of two bubbles occurs in three steps: 1. Bubbles collide with each other, trapping a small amount of liquid between them; 2. The liquid film separating the bubbles gradually disappears until it reaches a critical thickness; 3. The film ruptures and the bubbles combine to form one large bubble. This coalescence kernel is modeled by a collision rate of two bubbles and a collision efficiency related to the time required for coalescence:

$$a\left(V_i,V_j\right)=\left({\theta }_{ij}^T+{\theta }_{ij}^B+{\theta }_{ij}^S\right){\eta }_{ij},$$

(28)

The collision efficiency is modeled by comparing the time required for coalescence and the actual contact time during the collision:

$${\eta }_{ij}=e^{-t_{ij}/{\tau }_{ij}},$$

(29)

where:

$$t_{ij}={\left(\frac{{\rho }_l{r}_{ij}^3}{16\sigma }\right)}^{1/2}\ln \left(\frac{h_0}{h_f}\right),$$

(30)

and:

$${\tau }_{ij}=\frac{r_{ij}^{2/3}}{{\epsilon }_l^{1/3}},$$

(31)

where ${\rho }_l$ is the density of the liquid, $\sigma$ is the surface tension, $h_0$ is the initial film thickness, $h_f$ is the critical film thickness when rupture occurs, ${\epsilon }_l$ is the turbulent eddy dissipation of the liquid, and $r_{ij}={\left(\frac{1}{2}\left(\frac{1}{r_i}+\frac{1}{r_j}\right)\right)}^{-1}$ is the equivalent radius, $r_i$ and $r_j$ are the radii of bubbles i and j, respectively.

The turbulent contributions to collision frequency are modeled as:

$${\theta }_{ij}^T=\frac{\pi }{4}{\mathrm{F}}_{CT}{\left(d_i+d_j\right)}^2{\left[{\left(\sqrt{2}{\epsilon }_l^{1/3}{d}_i^{1/3}\right)}^2+{\left(\sqrt{2}{\epsilon }_l^{1/3}{d}_j^{1/3}\right)}^2\right]}^{1/2},$$

(32)

The buoyancy contribution to the collision frequency is calculated as:

$${\theta }_{ij}^B=\frac{\pi }{4}{\mathrm{F}}_{CB}{\left(d_i+d_j\right)}^2\left|\sqrt{\frac{2.14\sigma }{{\rho }_l{d}_j}+0.505g{d}_j}-\sqrt{\frac{2.14\sigma }{{\rho }_l{d}_i}+0.505g{d}_i}\right|,$$

(33)

where $d_i$ and $d_j$ are the diameters of bubbles i and j, respectively.

The shear contribution to collision frequency (${\theta }_{ij}^S$) is currently neglected.

2.5. Numerical Details

In this study, the bubble size distribution in an SSF mold was simulated using the commercial Computational Fluid Dynamics (CFD) software program ANSYS-FLUENT 20.2 [30]. The geometry was the same as Qi’s model [2] where a vacuum vessel with a diameter of 0.475 m and a height of 0.5 m, a submerged pipe with a diameter of 0.425 m and a height of 0.412 m, and a label with a diameter of 0.77 m and 0.86 m and the height of 0.5 m. Appropriate meshing techniques were employed to mesh the computational domain with structured grids (hexahedral mesh), resulting in 351,160 grids in total. Local grid refinements were also employed to better resolve the flow structures within the snorkel and vacuum chamber. In the current study, a mass flow boundary was applied to the gas inlet. A constant pressure (−4000 Pa) outlet condition was applied to the up of the vacuum chamber domain. The population balance equations were calculated by dividing the bubbles into ten bins. The minimum diameter of the bubble in group 9 was set to 1 mm and the ratio exponent was 1.2. Therefore, the bubble diameters of groups 0 to 9 are 12.1 mm, 9.19 mm, 6.96 mm, 5.28 mm, 4 mm, 3.03 mm, 2.30 mm, 1.74 mm, 1.32 mm, and 1 mm. The velocity–pressure linkage was conducted through the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) procedure. All the cases were carried out by transient simulations, and the double-precision solver was used. For all flow conditions, reliable convergence criterion based on the RMS (root mean square) residual of 1.0 × 10⁻⁵ was adopted for the termination of numerical calculations.

3. Results and Discussion

3.1. Validation of Numerical Model

In this study, the normalized concentration of NaCl with time is simulated for different inlet flows both with and without the PBM model. The concentration of two points in the ladle is detected, and the mixing time is calculated when $c_n=\frac{c_t}{c_{\infty }}\le \pm 5\%$. The simulated predicted values with and without the PBM technique are compared with the Qi [2] experimental values in Figure 2. In general, both numerical predictions have successfully captured the overall trend of the mixing time at various argon flow rates (i.e., 360 NL/h, 720 NL/h, 900 NL/h, and 1080 NL/h). For lower flow rates (i.e., from 360 to 900 NL/h), the mixing rate is almost linearly dependent on the flow rate. This indicates that a higher flow rate induces more effective circulating flow and reduces the time required for mixing in the system. Nevertheless, with a higher flow rate of 1080 NL/h, the drop in mixing rate becomes plateaus; showing additional gas injection has a relatively insignificant effect on enhancing the circulation. It may be due to the fact that as the gas flow rate increases, the bubble diameter inside the plume becomes larger and the bubble column moves significantly toward the wall. In addition, when the gas flow rate increases from 900 NL/h to 1080 NL/h, the bubble diameter does not increase further, and the stirring for the furnace does not improve. Comparing both results, one can easily notice that the PBM captures the embedded bubble dynamics and produces excellent agreement with the experimental data. Without bubble dynamic considerations, the numerical result tends to over-predict the mixing rate and the plateau point at a lower flow rate (i.e., 900 NL/h). Noticeable errors still persisted for the PBM simulations. As a result, the PBM could underestimate the bubble size, and the gas phase flow in the liquid phase primarily takes the form of bubbly flow. In contrast, when the argon flow rate is high enough in practice, the bubble deforms and, after its gas content hits a certain critical level, slug flow emerges. The experimental mixing time will be shorter than the predicted value because the slug flow accelerates the flow in comparison to the bubble flow. The results of the model with the PBM are more accurate than the results of the mixing time without the PBM. The mixing time trend of the PBM model is the same as the experimental mixing time trend at argon flow rates. Additionally, there is a maximum of 5% inaccuracy relative to the experimental value.

3.2. Flow Field

Figure 3 shows the liquid superficial velocity vector distribution at the center plane of the SSF at the flow rate of 900 NL/h. As depicted, following the entry into the furnace, the fluid flows upward along the right side of the snorkel and then begins to flow downward along the left side of the snorkel after reaching the free surface of the vacuum chamber, resulting in a significant circulation in the ladle, snorkel, and vacuum chamber. Comparing the two figures (a) and (b) in Figure 3, one can observe that the predicted liquid phase plume width by the PBM technique is significantly wider. Due to the vacuum status of the furnace, rising bubbles are subject to a vertical pressure gradient causing bubble expansion. The expanding bubble size increases the likelihood of bubble coalescence, forming larger bubbles.

In Figure 3, one can also observed that the superficial velocity at the center of the plume with PBM model consideration is relatively slower. Figure 4 shows a closer examination of the superficial liquid velocity distribution at the center phase above the gas injection inlet (i.e., X = 0.127 m). As depicted, it can be concluded that considering the bubble breakup and coalescence, it causes a more homogeneous velocity field in the X = 0.127 m plane. It may be due to the larger diameter of the bubble predicted by the PBM technique, which pushes the bubble to the center due to the lift effect. The fluid in the center will be carried upward by the bubbles, resulting in a more uniform flow field predicted by the PBM technique.

The determination of the circulation flow rate is related to the area of the upward flow and downward flow. Figure 5 shows the comparison of the effect of different flow rates on the area of the upward and downward flow at the Y = 0.3 m plane. Comparing predictions with or without the PBM model, it is clear that the predicted upward flow and downward flow areas are asymmetric without bubble dynamic consideration. Meanwhile, with the PBM model, the upward flow and downward flow areas are more evenly distributed. This is mainly because bubbles are subjected to turbulence at the plume boundary. When the bubble diameter is smaller than 3 mm, the bubble is shifted toward the center by a stronger lifting force on the bubble (see also in Equation (11)), and it causes the boundary of the upward flow to move toward the axis. Meanwhile, the results of Figure 4 are well interpreted from Figure 5. The same effect can be also observed in other flow rates as depicted in Figure 5.

Figure 6 illustrates the distribution of bubble diameters at the Y = 0.3 m plane for different flow rates. It can be seen from the figure that there is an annular region on the right side of the plane with a larger bubble diameter, which is mainly due to the gas inlet in the negative Y direction in this region. Due to a large number of bubbles in this region, bubbles are subject to rigorous interactions, causing more coalescence and forming larger bubbles. Meanwhile, in the downward region (i.e., left-hand side), smaller diameter bubbles are prevalent as bubble number density is low in the region, reducing the change of coalescence. In general, as more bubbles are injected at a higher flow rate, the overall bubble diameter increases with the increase in flow rate.

Figure 7 shows the comparison of the free surface level height with respect to different flow rates and bubble dynamic considerations. As depicted, it can be observed that there is some effect between the anticipated free surface height and the effect of bubbles on breakup and coalescence. The experimental data of Qi [2], which both exhibit a trend of the low left side and high right side, are in agreement with the free surface fluctuations anticipated by this model. It can be found that at each gas flow rate, the smaller fluctuations in the liquid level are predicted by the PBM technique. This is mainly because of a small diameter bubble predicted by the PBM technique at the edge of the plume, with some of the bubbles moving toward the center. The bubbles predicted by the technique without PBM are concentrated and do not move toward the center. As a result, a large number of bubbles escape when the liquid level is reached, causing greater fluctuations.

3.3. Local Bubble Size Distribution

Several previous articles have demonstrated that the turbulent eddies collision mechanism is prevalent in fully turbulent flow [31] and have concluded that the bubble coalescence process can be broken into three steps: the formation of a liquid film between two approaching bubbles by entraining liquid from the continuous phase; the continuation of the bubbles’ approach while the liquid film continues to drain thin; and the breakage of the liquid film when it is sufficiently thin for the bubbles to coalesce. Since the bubble plume inside the furnace is very wide, there is a high chance of bubble coalescence. At the edges of the bubble plume, the bubbles are subject to turbulent vortices, causing the rupture of the bubbles due to turbulent impact at the plume edges.

Figure 8 shows the local bubble size distribution at the center plane (i.e., Z = 0) of the SSF under different flow rates. From the figure, it can be seen that once a bubble enters the furnace interior (i.e., inlet diameter of 3 mm), its diameter gradually increases as it rises due to expansion. As mentioned before, at the plume’s central region, which is the dominant position to coalesce, the high chance of collision tends to increase the bubble size in the region.

Figure 9 shows the distribution of number density and mean bubble diameter at four different cross-section lines. Notice that line 4 represents the intersection line at the free liquid surface (i.e., plane Z = 0). The bubble characteristics are analyzed from two perspectives in the following. As shown in Figure 9a,b, both the number density and diameter of bubbles reach the maximum value at the location of X = 0.05–0.1 m. It is because a large amount of gas enters from the inlet (i.e., X = 0.127 m) and moves vertically along the Y-direction, causing a large number of bubbles to accumulate in this region and move toward the center of the ladle. In Figure 9b, the bubble diameter has a peak at X = 0.125 m, mainly because of the spread of the bubble plume and bubble expansion. At a flow rate of 360 NL/h, the flow rate is relatively low for the gas to carry liquid to reach here. Therefore, the peak of bubble diameter here is less prominent compared to other higher flow rates. Line 3 is located at the snorkel (see in Figure 9c). The bubbles’ number density and diameter both reach the maximum value at X = 0.05–0.1 m. Similar to other locations, the peak bubble diameter increases with the increasing flow rate. Figure 9d shows bubble number density and diameter at the free liquid surface. A large number of bubbles escape at X = 0.05–0.1 m and bubble breakup occurs. In general, the average bubble diameter increases as it rises with the liquid vertically and reaches the maximum size at the free surface.

Figure 10 depicts the interfacial area at different gas flow rates, from which it can be seen that the main bubble size in the furnace is between 4 and 5 mm. It is around a 33% to 67% increment compared to the bubble diameter at the inlet (i.e., 3 mm). At a flow rate of 360 NL/h, the peak of the interfacial area in the furnace is 2.5 m², which corresponds to a bubble size of 4 mm. However, the peak of the interfacial area at other flow rates corresponds to a bubble size of about 5 mm, which indicates that the coalescence rate in the furnace is dominant in comparison to the breakup process. As the flow rate increases, the available interfacial area for the decarburization reaction to occur also increases, and theoretically, the decarburization that occurs on the surface of Ar bubbles is more efficient. In the next step, we will further explore the relationship between the decarburization behavior of the Ar bubble surface and the interfacial area.

4. Conclusions

In this paper, a three-dimensional mathematical model of two-fluid multiphase flow with coupled CFD-PBM method is established to study the bubble size distribution in the SSF furnace, which is characterized by the ability to study the effects of bubble breakup and coalescence on the flow field. From the results, the following conclusions can be drawn: (1) The mixing times of the with and without PBM model are calculated and compared with the experimental values, and both results are found to be in good agreement with the experimental results, but the results of the with PBM model are closer to the experimental values. (2) The flow field of the PBM model is more homogeneous, and the distribution of the upward flow and downward flow in the same cross-section is more reasonable. (3) The bubble diameter in the furnace increases with the increase of the height of the liquid surface, and the majority of the bubble diameter in the furnace is larger than 3 mm.