Two-Way PBM–Euler Model for Gas and Liquid Flow in the Ladle

Ladle metallurgy is an important steelmaking technology in high-quality steel production. The blowing of argon at the ladle bottom has been applied in ladle metallurgy for several decades. Until now, the issue of breakage and coalescence among bubbles was still far from being solved. In order to have a deep insight into the complex process of fluid flow in the gas-stirred ladle, the Euler–Euler model and population balance model (PBM) are coupled to investigate the complex fluid flow in the gas-stirred ladle. Here, the Euler–Euler model is applied to predict the two-phase flow, and PBM is applied to predict the bubble and size distribution. The coalescence model, which considers turbulent eddy and bubble wake entrainment, is taken into account to determine the evolution of the bubble size. The numerical results show that if the mathematical model ignores the breakage of bubbles, the mathematical model gives the wrong bubble distribution. For bubble coalescence in the ladle, turbulent eddy coalescence is the main mode, and wake entrainment coalescence is the minor mode. Additionally, the number of the bubble-size group is a key parameter for describing the bubble behavior. The size group number 10 is recommended to predict the bubble-size distribution.


Introduction
Because of the need to meet the demand for high-quality steel, ladle metallurgy has received attention from many researchers [1,2]. Ladle refining is a popular steelmaking technology that has the following features: degassing, desulphurization, adjustment of the chemical composition of alloy elements, inclusion removal, as well as temperature homogenization [1][2][3][4]. Among ladle refining reactors, blowing argon at the ladle bottom is a widely used technology [4]. Due to the buoyancy, the bubbles rise upward and induce the recirculation flow to provide effective mixing and then escape from the free surface [2,3].
With the rapid development of computer hardware and software, some researchers devoted themselves to computational fluid dynamics (CFD). Currently, the relevant mathematical models include the quasi-single-phase model [5], the Euler-Lagrange model, and the Euler-Euler model [5][6][7]. The computational load of the quasi-single-phase model is the lowest among them, but the slip velocity between the liquid phase and the gas phase must be specified. The precondition of the Euler-Lagrange model is that the gas volume fraction is less than 12%. Additionally, the fluid flow is at a steady state, but the solution of the bubble trajectory equation is the function of the time, so the governing equations of the fluid phase have to be the unsteady differential equations. The Euler-Euler model takes into account the interaction between the gas phase and the liquid phase, and one set of partial differential equation equations corresponds to each phase, so the related computational load is heavy.

Assumptions
The Euler-Euler model and PBM are used to describe the two-phase flow in the ladle with bottom blowing. In order to simplify the complex transfer process in the ladle, some assumptions must be applied in the mathematical model. (1) Both the gas and the liquid are the incompressible Newtonian fluid [1][2][3][4][5][6][7][8]; (2) The temperature of the molten steel stays constant [5]; (3) The effect of the top slag on the fluid flow is so weak that it can be ignored [5][6][7]; (4) The bubble is spherical [5,6]; (5) There are no chemical reactions in the molten steel [7,9]; (6) There is only binary breakage of the bubble [3,8].

Governing Equations 2.2.1. Eulerian Multiphase Hydrodynamic Equations
The Euler-Euler model has a separate set of continuity and momentum equations for the gas phase and liquid phase, respectively [5][6][7][8][9]. The coupling of the phases comes from the shared interfacial momentum exchange source terms and the pressure field. The related physical parameters in the numerical simulation are listed in Table 1. The continuity equation and the momentum conservation model can be expressed as follows: where α k (α g + α l = 1), ρ k and u k are the volume fraction, the density and velocity of gas (k = g), and liquid (k = l), respectively. p is the pressure, g is the gravity acceleration, and R k is the interfacial momentum exchange source term between the gas and the liquid. µ eff,l is the effective viscosity of liquid phase, which is determined with the standard k − ε turbulence model. The effective viscosity of gas phase, µ eff,g , is calculated using µ eff,l . The interfacial force source terms are expressed as: where R l ( R g ) denotes the interfacial momentum exchange source terms from the gas (liquid) phase to the liquid (gas) phase. F D , F L , F VM , and F TD are the drag force, the lift force, the virtual mass force, and the turbulence dispersion force, respectively. Figure 1 shows the complex physical phenomena.  (1) The drag force, D F , which is the dominant interphase force, provides the resistance to the flow due to the motion of the bubbles relative to the molten steel, acts in the opposite direction to motion of the bubbles, and can be calculated with [21,22,23,24,25]:  (1) The drag force, F D , which is the dominant interphase force, provides the resistance to the flow due to the motion of the bubbles relative to the molten steel, acts in the opposite direction to motion of the bubbles, and can be calculated with [21][22][23][24][25]: where the Sauter mean diameter of bubbles d g is calculated from PBM. The drag coefficient C D , which is given by Tomiyama [23], Re g min 1 + 0.15Re g 0.687 , 3 , 8Eo g 3(Eo g + 4) (5) is related to the bubble Reynolds number (Re g = ρ l u g − u l d g /µ l ), the bubble Eotvos number (Eo g = g ρ l − ρ g d g 2 /σ), and the bubble Weber number (We g = ρ l u g − u l 2 d g /σ).
µ l is the viscosity of molten steel. σ is the surface tension coefficient. g is the magnitude of gravity.
(2) The interaction between the bubbles and the liquid shear field results in the lift force, F L , which is perpendicular to the relative motion of the bubbles and the liquid. The lift force is the function of the vector product of slip velocity and the curl of liquid velocity: where the lift coefficient C L is set to 0.1 [6].
(3) According to the Drew model [24,26], the virtual mass force, F VM , arises from the inertia of molten steel relative to the acceleration of bubbles, and is defined as follows: where the virtual mass coefficient C VM is set to 0.5., d g /dt and d l /dt operators denote the substantial derivatives in the gas phase and the molten steel.
(4) The turbulent dispersion force, F TD , comes from turbulent fluctuations in the molten steel velocity due to bubble-eddy interactions, which is given by Burns equations [27], where C D is the drag force coefficient, and µ t,l is the turbulent viscosity of liquid phase.

Population Balance Model
In order to predict the bubble-size distribution in the gas-string ladle, the population balance model (PBM) is employed, and the classes method (CM) [26,[28][29][30] is implemented for solving PBM in OpenFOAM. The bubble distribution is represented by bubble classes, and the coalescence rate and breakage rate are transformed into the birth rate and death rate in each class. PBE for the ith bubble class can be given as follows: where B C,i and B B,i are not only the birth rates caused by the breakage and coalescence of bubbles, but also the death rates caused by the breakage and coalescence of bubbles. Bubble number density, n i , is the number of bubbles in group i per unit volume. PBE has the same form as the generic convection-diffusion equation for multiphase flow.
Here, the bubble number density and bubble volume fraction of the dispersed phase have the following relationship: where f i is the volume fraction of size group i in the total dispersed phase fraction, f i = α i /α g. They satisfy ∑ The source terms on the coalescence and breakage of the bubble can be modeled as where a(v − v , v ) represents the coalescence rate between the bubble of size v and the bubble of size v , b(v ) is the breakage rate of a bubble of size v , m(v ) is the number of daughter bubbles generated from the breakage of a bubble of size v , and p(v , v ) is the probability density function for a bubble of size v, generated by the breakage of a bubble of size v .
(1) Coalescence kernel functions In the previous literature, many researchers only considered the turbulent eddy mechanism [8]. However, the coalescence of bubbles resulting from wake entrainment is significant for the formation of large bubbles. Thus, there are the turbulent eddy mechanism and the bubble wake entrainment in the current paper.
The coalescence rate a(v i , v j ) between a bubble of size v i and a bubble of size v j is usually expressed as where The bubble collisions frequency resulting from turbulent eddy is calculated with where u ij = (u 2 i + u 2 j ) 1/2 is the characteristic velocity of the collision between bubbles of size d i and d j ,. If the bubble size d i is equal to the eddy size, the mean turbulent velocity u i of the bubble is equal to the mean turbulent velocity u λ = √ 2(ε l λ) 1/3 of the eddy. The bubble coalescence efficiency model is based on a phenomenological analysis. The coalescence between two bubbles depends on the ratio of the contact time τ ij and the coalescence time t ij for the drainage of the liquid film between them to reach a critical rupture thickness [8,27,31]. The coalescence probability function of the bubble of size v i and the bubble of size v j is written as where The bubble collision frequency resulting from wake entrainment can be calculated with where the parameter K is set to be 6.0., u i,W = 0.71(gd i ) 1/2 is the rising velocity of the leading bubble. The big bubbles have the effective wake region for the bubble to coalesce, so the parameter Θ is valid for the bubbles larger than d c and can be expressed as The bubble coalescence efficiency is given by Hibiki and Ishii [32], (2) Breakup kernel functions In the turbulent flow, the bubble breakage rate is determined by the interaction among bubbles with a size close to the turbulent eddy. Therefore, the turbulent eddy, whose size is much smaller than the bubble size, does not have enough energy to break the bubble. On the other hand, the turbulent eddy whose size is greater than the bubble size can merely carry bubbles and does not cause bubble breakage. The breakage rate is described as a product of collision frequency and breakage efficiency. In the current paper, the breakup model is derived from theories of isotropic turbulence [26,31]. The breakup rate function proposed by Lehr [35] can be applied to describe the bubble breakup. In this case, a bubble of volume v j breaks into two bubbles, one of volume v i and the other of volume v j − v i .
where T = ( σ ρ l ) 2/5 1 According to the definition of breakage kernel, b(v ) and p(v, v ) can be obtained with According to the bubble sizes and volume fractions for various bubble classes, the Sauter mean diameter is defined as [36]

Numerical Details
The numerical simulation for the gas-stirred ladle is carried out using the free, opensource software OpenFOAM (version 6.0.0). The computational domain includes the two-phase flow region and free board above the liquid in the ladle. The free-board height is one-third of the bath height. All the computational domains are covered with nonuniform hexahedral grids, as shown in Figure 2. The convergence criterion is that the root mean-square normalized residual is less than 10 −5 . The time step is set to 0.001 s. There are five types of boundary conditions. The veloc total gas flow rate at the nozzle of ladle bottom, the pressu slip boundary condition is applied to the solid wall, and the gas phase and liquid phase, as shown in Figure 2. The det Table 1.
Sano and Mori's empirical equation is used to calculate nozzle exit, where dn is the nozzle diameter, and Vg is the gas flowrate.

Validation of Grid Independence
In order to ensure that the numerical results are not th resolved grids, we carried out the grid sensitivity experimen an inlet flow rate of 9000 mL/s. Based on the five grid syste mum velocity and maximum gas volume fraction of molte There are five types of boundary conditions. The velocity inlet is determined by the total gas flow rate at the nozzle of ladle bottom, the pressure outlet is specified, the no-slip boundary condition is applied to the solid wall, and the pressure outlet is used for the gas phase and liquid phase, as shown in Figure 2. The detailed information is listed in Table 1.
Sano and Mori's empirical equation is used to calculate the initial bubble size at the nozzle exit, where d n is the nozzle diameter, and V g is the gas flowrate.

Validation of Grid Independence
In order to ensure that the numerical results are not the spurious artifacts of poorly resolved grids, we carried out the grid sensitivity experiments for the industrial ladle with an inlet flow rate of 9000 mL/s. Based on the five grid systems, Figure 3 gives the maximum velocity and maximum gas volume fraction of molten steel at 1 m below the free surface.
The variation of the velocity and the variation of the gas volume fraction is 0.179 m/s and 0.002 between the grid numbers 224,636 and 275,640. When the grid is refined from 275,640 to 423,972, the variation of the velocity between other adjacent grids is 0.012 m/s, 0.002 m/s, and 0.006 m/s, and the gas volume fraction differences between other adjacent grids are 0.006, 0.002, and 0.005. There are no perceptible differences between the velocity and the gas volume fraction in the cases of 338,778 grids and 423,972 grids. Therefore, the 338,778 grids are fine enough for the coupling of the flow field and PBM.

Validation of Bubble-Size Distribution
The bubble-size distribution and the flow field reported in the literature [33,34] are applied to validate the numerical result. Table 1 gives the related geometrical parameters and physical parameters.
In Anagbo's experiment [34], the water model was a one-sixth scale model of a 150-ton steelmaking ladle. The experimental apparatus consisted of a cylindrical plexiglass tank containing deionized water with a depth of 0.4 mm and a porous plug located at the center of the bottom. The inlet flow rate was 600 mL/s. Additionally, the axial bubble diameter was measured using the double-contact electro-resistivity probe. Figure 4 gives the difference between the experimental data and the numerical results: (1) The bubble behaviors consist of the bubble breakage and the coalescence among bubbles, and the coalescence among bubbles has two modes: turbulent eddy coalescence and wake entrainment. Model E only considers the coalescence and ignores the breakage, so the bubble diameter predicted with model E is the greatest among the models, and the bubble diameter increases with the increase in the axial distance. Such a phenomenon disagrees with the experimental data. The bubble size predicted with model B is greater than that by model C, so the effect of the turbulent eddy on the coalescence among bubbles is greater than wake entrainment. (5) Model D considers turbulent eddy coalescence and wake entrainment coalescence, so the bubble size predicted with model D is greater than that by model B (or model C). The bubble maximum diameter predicted with model (D) is 0.03 m, and the related experimental data are 0.034 m. In this way, the relative error is 11.8%. Several reasons lead to such a difference: (1) There is a spherical bubble assumption in the mathematical model. (2) It is difficult for the probe to identify spherical bubbles and non-spherical bubbles.   Figure 5 shows that the numerical results agree well with experimental data [33]. In the water model of a geometrical scale of 1/10, a flush-mounted orifice with a 4 mm inner  Figure 5 shows that the numerical results agree well with experimental data [33]. In the water model of a geometrical scale of 1/10, a flush-mounted orifice with a 4 mm inner diameter is placed at the center of the ladle bottom, and the measurements are made using a combined electrical probe technique and the laser Doppler anemometer (LDA). Figure 5a indicates that the predicted maximum gas volume fraction is 0.44, the related experimental value is 0.41, and the relative error is 7%. Figure 5b indicates the maximum difference of liquid axial velocity between the experimental data and the predicted result is 0.04 m/s, and the maximum relative error is 12%.  Figure 5 shows that the numerical results agree well with experimental data [33]. In the water model of a geometrical scale of 1/10, a flush-mounted orifice with a 4 mm inner diameter is placed at the center of the ladle bottom, and the measurements are made using a combined electrical probe technique and the laser Doppler anemometer (LDA). Figure  5a indicates that the predicted maximum gas volume fraction is 0.44, the related experimental value is 0.41, and the relative error is 7%. Figure 5b indicates the maximum difference of liquid axial velocity between the experimental data and the predicted result is 0.04 m/s, and the maximum relative error is 12%.

Bubble-Size Distribution for Different Bubble-Size Classes
In order to balance the computational cost and the computational accuracy, the Sauter mean diameter of argon bubbles at the ladle axis is applied in the four bubble-size

Bubble-Size Distribution for Different Bubble-Size Classes
In order to balance the computational cost and the computational accuracy, the Sauter mean diameter of argon bubbles at the ladle axis is applied in the four bubblesize groups. Based on the Sauter mean diameter, the bubble size is divided into G groups, v i+1 = θv i (i = 1, 2, 3, · · · G−1).
The Sauter mean diameter of argon bubbles has the following features at the ladle axis, as shown in Figure 6: (1) The diameter of argon bubbles at the ladle axis follows a similar trend for the four bubble-size groups. The bubble diameter decreases in the initial rising stages and then remains unchanged gradually. Such a distribution feature was observed by Anagbo [34]. (2) In the four cases, the diameter of argon bubbles in the case of G = 5 is the smallest because the bigger geometric discretization factor θ leads to the greater error.
(3) The distribution of the diameter of argon bubbles is almost similar in the cases of G = 10, 15, 20. The maximum relative error of the Sauter mean diameter of argon bubbles between M = 10 and M = 20 is only 8.8%. Thus, the bubble size is divided into ten groups in the following numerical calculation. In the previous paper [3,11], the geometric discretization factor on bubble size was not considered, and the bubble was directly divided into fixed groups in the numerical calculation. Such a simplification may lead to a bad result, such a simplification may lead to a bad result that the predicted bubble size is less than the actual bubble size.

Flow Field in the Industrial Ladle
The flow field of molten steel and the argon volume fraction in the ladle are shown in Figure 7. Argon bubbles leave the ladle bottom as a strong jet flow and spread outward until they reach the free surface. As argon bubbles float upward, the gas volume fraction decreases gradually, but the plume region increases gradually. Because of the buoyancy of argon bubbles, the molten steel also flows upwards. Consequently, there are vortexes on both sides of the gas jet flow, and these vortexes are near the free surface. Some researchers [37,38] used the Euler-Euler approach to describe the fluid flow in the gasstirred ladle and indicated that the width of the plume remained almost unchanged along the vertical direction. Such a result is different from the experimental result. The reason is that their mathematical model ignored the effect of the turbulent dissipation force on the fluid flow.

Flow Field in the Industrial Ladle
The flow field of molten steel and the argon volume fraction in the ladle are shown in Figure 7. Argon bubbles leave the ladle bottom as a strong jet flow and spread outward until they reach the free surface. As argon bubbles float upward, the gas volume fraction decreases gradually, but the plume region increases gradually. Because of the buoyancy of argon bubbles, the molten steel also flows upwards. Consequently, there are vortexes on both sides of the gas jet flow, and these vortexes are near the free surface. Some researchers [37,38] used the Euler-Euler approach to describe the fluid flow in the gasstirred ladle and indicated that the width of the plume remained almost unchanged along the vertical direction. Such a result is different from the experimental result. The reason is that their mathematical model ignored the effect of the turbulent dissipation force on the fluid flow.

Disperse Phase in the Industrial Ladle
It can be observed in Figure 8 that the bubble diameter and the gas volume fraction have some similar features in the plume region of the industrial ladle: (1) The bubble diameter and the gas volume fraction are symmetrically distributed on both sides of the ladle centerline. (2) The bubble diameter and the gas volume fraction along the plume centerline are greater than that in other regions. (3) The bubble diameter and the gas volume fraction gradually decrease with the increase in the axial distance from the ladle bottom.

Disperse Phase in the Industrial Ladle
It can be observed in Figure 8

Disperse Phase in the Industrial Ladle
It can be observed in Figure

The Effect of Bubble Size on Flow Field in the Industrial Ladle
In the past references [5][6][7], there were many research works based on the cons bubble size or the turbulent eddy coalescence among bubbles. Figure 9 gives the dif ences among the mathematical models with a constant bubble size, the mathemat model with a bubble turbulent eddy coalescence (model B), and the mathematical mo with a bubble turbulent eddy coalescence and wake entrainment coalescence (model Figure 9a shows that, in the case of the constant bubble size, there is an abrupt increas the axial velocity of molten steel. But an abrupt increase in the axial velocity of mo steel is not obvious in the cases of model B and model D. Such a difference comes fr the effect of bubble breakage and bubble coalescence on the flow of molten steel. Fig  9b also shows that, with the increase in the distance from the ladle bottom, the gas volu fraction in the case of the constant bubble size decreases faster than that in the cas model B and model D. Such phenomena indicate that the interaction between the f flow and the bubble breakage/coalescence can prompt the uniformity of fluid velocity bubble distribution. Figures 4 and 9c show the importance of the wake entrainment lescence of bubbles. In other words, the bubble behaviors can be described accurately o when the bubble breakage and the bubble coalescence (turbulent eddy and wake entr ment) are considered in the mathematical model.

The Effect of Bubble Size on Flow Field in the Industrial Ladle
In the past references [5][6][7], there were many research works based on the constant bubble size or the turbulent eddy coalescence among bubbles. Figure 9 gives the differences among the mathematical models with a constant bubble size, the mathematical model with a bubble turbulent eddy coalescence (model B), and the mathematical model with a bubble turbulent eddy coalescence and wake entrainment coalescence (model D). Figure 9a shows that, in the case of the constant bubble size, there is an abrupt increase in the axial velocity of molten steel. But an abrupt increase in the axial velocity of molten steel is not obvious in the cases of model B and model D. Such a difference comes from the effect of bubble breakage and bubble coalescence on the flow of molten steel. Figure 9b also shows that, with the increase in the distance from the ladle bottom, the gas volume fraction in the case of the constant bubble size decreases faster than that in the case of model B and model D. Such phenomena indicate that the interaction between the fluid flow and the bubble breakage/coalescence can prompt the uniformity of fluid velocity and bubble distribution. Figures 4 and 9c show the importance of the wake entrainment coalescence of bubbles. In other words, the bubble behaviors can be described accurately only when the bubble breakage and the bubble coalescence (turbulent eddy and wake entrainment) are considered in the mathematical model.

Conclusions
In order to solve the issue of breakage and coalescence among the bubbles in ladle metallurgy, this work proposed a strategy that the Euler-Euler model and population balance model are coupled to investigate the complex fluid flow in the gas-stirred ladle. Research shows that there is coalescence and the breakage of bubbles in the molten steel in the ladle. Bubble breakage is a necessary factor. If the mathematical model ignores this factor, the mathematical model gives the wrong bubble distribution. Furthermore, the bubble coalescence in the ladle consists of turbulent eddy coalescence and wake entrainment coalescence. The turbulent eddy coalescence is the main mode, and the wake entrainment coalescence is the minor mode. The bubble breakage, turbulent eddy coalescence, and wake entrainment coalescence are introduced to describe bubble behavior in order to ensure that PBM can give the transfer behavior in the ladle accurately. Meanwhile, it is clarified that the number of the bubble-size group is the key parameter to describe the bubble behavior, and the size group number 10 is recommended to predict the bubble-size distribution in the ladle.

Conclusions
In order to solve the issue of breakage and coalescence among the bubbles in ladle metallurgy, this work proposed a strategy that the Euler-Euler model and population balance model are coupled to investigate the complex fluid flow in the gas-stirred ladle. Research shows that there is coalescence and the breakage of bubbles in the molten steel in the ladle. Bubble breakage is a necessary factor. If the mathematical model ignores this factor, the mathematical model gives the wrong bubble distribution. Furthermore, the bubble coalescence in the ladle consists of turbulent eddy coalescence and wake entrainment coalescence. The turbulent eddy coalescence is the main mode, and the wake entrainment coalescence is the minor mode. The bubble breakage, turbulent eddy coalescence, and wake entrainment coalescence are introduced to describe bubble behavior in order to ensure that PBM can give the transfer behavior in the ladle accurately. Meanwhile, it is clarified that the number of the bubble-size group is the key parameter to describe the bubble behavior, and the size group number 10 is recommended to predict the bubble-size distribution in the ladle.