Water Model and CFD-PBM Coupled Model of Gas-Liquid-Slag Three-Phase Flow in Ladle Metallurgy

Abstract

In ladle metallurgy, the role of argon gas purging in liquid steel is well recognized. One key aspect of this process is the prediction of gas volume fraction and gas bubble diameter. Calculation of bubble size is important as it is significantly modified by the injection system and gas flow rate, and it also directly determines the buoyancy which influences the physical mixing process and surface reactions. In this study, a porous plug was used for gas injection as the common practice in the steel industry and a one-third scale water model was established. Then a numerical model based on the Eulerian multiphase model was established and the population balance model (PBM) was used to calculate the gas bubble size distribution. The bubble coalescence and breakage were included and the phase interactions were coupled with the PBM to consider the effect of bubble diameter on the fluid flow. A user defined scalar (UDS) transport equation was added to simulate the solute transport in ladle to study the mixing efficiency. The mixing time, wall shear stress and slag entrapment probability were taken into consideration to find a suitable plug position to balance the mixing efficiency and steel cleanness. The results show that the mixing time decreases with increasing the plug radial distance and the maximum wall shear stress appears when plug radial distance is 0.67 R.

1. Introduction

In dispersed gas-liquid flow, the gas bubble diameter distribution plays a significant role in the phase structure and interphase forces, and determines the multiphase hydrodynamic behaviors. In the steel refining process, argon gas stirring has been widely employed to homogenize the chemical composition of alloy elements and temperature, remove inclusions and enhance the rates of refining reactions. The bubble plume is successively formed at the exit of the bottom nozzle. Then they rise upward and entrain the surrounding molten steel into their wakes. During the bubble flotation, the aggregation and breakage happens and causes a wide range bubble size distribution. Still, it is a significant unresolved issue in ladle with gas stirring. Calculation of the bubble size distribution and understanding of bubble aggregation and breakage after bubble come out of the plug are important for reasons as follows: the bubble size is significantly modified by the injection system and flow rates, and it will directly influence the buoyancy and also influence the physical phenomena. It is expected that fine and well-distributed bubbles are good for inclusion removal and homogenizing. With a high gas flow rate, the rising gas bubbles can break through the slag layer, strengthen the interfacial reaction, cause the exposure of molten steel and form a large spout eye to serve the desulphurization. On the other hand, a relatively calm flow without a broken slag layer is needed for inclusion floatation. While surface turbulence and strong interaction between slag and steel are needed to enhance the efficiency of desulphurization, the pick-up of oxygen and nitrogen from atmosphere and the entrapment of slag caused by slag-metal interface fluctuation will possibly happen. Therefore, a comprehensive three-phase model with bubble size prediction is needed to understand the bubbly flow as well as the interfacial behavior of molten steel and slag.

Over the years, many cold physical models^{1,2,3,4,5,6,7)} and numerical simulation^{8,9,10,11,12,13,14)} have been taken to investigate the flow and mixing phenomena in the gas stirred ladle system. Several mathematical models have been developed to study the flow in ladle by both the Eulerian-Lagrange^8,9,14) and Eulerian-Eulerian^15,16,17) approach. Li et al.¹³⁾ developed a transient three-dimensional and three-phase model by Multiphase Volume of Fluid (VOF) method to study the behavior of slag layer. The gas phase was treated as a continuous phase and phase interactions were not included. Cloete et al.¹⁴⁾ developed a mathematical model by employing the Lagrangian discrete phase model (DPM) to describe the bubble plume and the Eulerian VOF model for tracking the free surface of the melt without considering the bubble aggregation and breakage. Liu et al.⁸⁾ developed a similar mathematical model by DPM-VOF coupled method to investigate the mixing time and slag layer behavior. The Population Balance Model (PBM) was first formulated for chemical engineering purposes by Hulburt et al.¹⁸⁾ While the PBM can predict the size distribution of dispersed phases in a wide variety of particulate processes. Recently, it is reported that it is an effective method to simulate the bubble size distribution.^10,17,19) However, in these studies,^9,10,11,12) a flat and free slip surface was always assumed as well as the slag open eye formation and the interfacial behavior were usually neglected. According to the studies of Cloete et al.,¹⁴⁾ it showed that the flat surface assumption is not valid for the investigation of mixing time and interfacial behavior and it is necessary to consider the effect of a dynamic free surface in the ladle metallurgical process research. Furthermore, with bubbles injecting, both the low-Reynolds number and high Reynolds number effects were included in the flow of ladle but were not always considered in previous works.

Detailed knowledge of bubble floatation, mixing, slag behavior and erosion of the ladle lining is necessary to make structural optimization and process optimization to balance the mixing efficiency and the cleanness of steel. It is reported that the mixing time and also the wall shear stress are sensitive to the plug configuration. The present study was undertaken to enhance the understanding of aggregation and breakage after bubble was injected into the liquid through the porous zone as well as the slag layer behavior and mixing phenomena in ladle. Firstly, a one-third scale water model was established. The slag layer was simulated using the bean oil and the gas bubble was formed using the nitrogen gas coming through a porous plug. A high speed camera was used to monitor the bubble floatation and an open software was taken to quantify the bubble size. Next, an Eulerian multiphase flow model coupled with the PBM through the Discrete Method along with the RNG k-ε turbulence model was developed to study the gas-liquid-slag three-phase flow in ladle. Comparisons of gas distribution, bubble size distribution, slag eye size and solute transport with water model experimental data were provided, showing the applicability and accuracy of the model for the fluid flow, slag layer behavior and mixing phenomena inside the ladle.

2. Model Formulation

The discrete population balance model in the present work is integrated with the Eulerian multiphase model. In bubbly flows, interfacial momentum transfer exhibits a dominant effect in the multiphase momentum equations. The total interfacial forces considered in the present study can be categorized into three main terms: drag, lift and virtual mass force. The RNG k-ε turbulence model^20,21) was used in present model. It is similar in form to the standard k-ε model but the accuracy for rapidly strained flows and swirling flows are improved. And also it provides an analytically-derived differential formula for effective viscosity that accounts for low-Reynolds number effects. The turbulence interaction was included using the Sato model.²²⁾ To calculate the non-uniform bubble size distribution, the Population Balance Model by the Discrete Method^23,24,25) was employed. The Luo aggregation kernel and Luo breakage krenel²⁶⁾ were both included. The effect of bubble diameter on flow field was included by transferring the bubble Sauter mean diameter data between the PBM and the momentum equation. Figure 1 shows the gas stirring phenomenon in ladle and the procedure details of solution methodology of this model.

Fig. 1.

Solution methodology of this model.

2.1. Eulerian Multiphase Model

In the current work, the mathematical model is based on the Eulerian multiphase approach. In this approach, pressure and turbulence fields are shared among phases whereas phasic equations are solved for continuity and momentum equations. The continuity equation for phase q can be written as:   

$$\frac{\partial }{\partial t}({\alpha }_q{\rho }_q)+\nabla \cdot ({\alpha }_q{\rho }_q{\overrightarrow{u}}_q)=0$$(1)

where α_q, ρ_q and ${\overrightarrow{u}}_q$ are the volume fraction, density and velocity of the q^th phase respectively. The momentum balance for the q^th phase is:   

$$\begin{array}{l}\frac{\partial }{\partial t}({\alpha }_q{\rho }_q{\overrightarrow{u}}_q)+\nabla \cdot ({\alpha }_q{\rho }_q{\overrightarrow{u}}_q{\overrightarrow{u}}_q)=\\ -{\alpha }_q\nabla p+\nabla \cdot {\overline{\tau }}_q+{\alpha }_q{\rho }_q\overrightarrow{g}+{\overrightarrow{F}}_{drag,q}+{\overrightarrow{F}}_{lift,q}+{\overrightarrow{F}}_{VM,q}\end{array}$$(2)

where ${\overline{\tau }}_q$ , p and $\overrightarrow{g}$ are the q^th phase stress-strain tensor, pressure and gravity acceleration respectively. ${\overrightarrow{F}}_{drag,q}$ is the drag force, ${\overrightarrow{F}}_{lift,q}$ is the lift force and ${\overrightarrow{F}}_{VM,q}$ is the virtual mass force acting on the q^th phase. The drag force between phases can be described as:   

$${\overrightarrow{F}}_{drag,q}=\frac{3f{\mu }_q}{L}{A}_i({\overrightarrow{u}}_p-{\overrightarrow{u}}_q)$$(3)

where A_i is the interfacial area, L and μ_q are the diameter and viscosity of q^th phase respectively, f is the drag function. In this work, the Schiller-Naumann model²⁷⁾ coupled with the population balance model was used. The diameter and interfacial area can be calculated directly from the population balance variables so the influence of bubble diameter on the fluid flow is included.

The lift force is mainly due to velocity gradient in the primary phase flow field and the virtual mass force is caused by the acceleration of a phase relative to another. They are defined as:   

$${\overrightarrow{F}}_{lift,q}=C_l{\rho }_q{\alpha }_p({\overrightarrow{u}}_q-{\overrightarrow{u}}_p)\times (\nabla \times {\overrightarrow{u}}_q)$$(4)

  

$${\overrightarrow{F}}_{VM,q}=-0.5{\alpha }_p{\rho }_q\left(\frac{d_q{\overrightarrow{u}}_q}{dt}-\frac{d_p{\overrightarrow{u}}_p}{dt}\right)$$(5)

where C_l is the lift force coefficient. The term d_q/dt denotes the substantial derivative of the q^th phase.

The RNG-based k-ε turbulence model is derived from the instantaneous Navier-Stokes equations, using a mathematical technique called “renormalization group” (RNG) methods. It has greater potential to give accurate predictions for complex flows. More details can be found in Yakhot et al.²¹⁾ The effect of swirl and also the low-Reynolds number effects are included in the RNG k-ε turbulence model by modeling the effective viscosity. In the high-Reynolds number limit, the turbulent viscosity for mixture, μ_t,m is computed using:   

$${\mu }_{t,m}={\rho }_m{C}_{\mu }\frac{k^2}{\epsilon }+{\mu }_{t,p}$$(6)

where C_μ=0.0845. The extra term μ_t,p on the right side is included when the influence of the dispersed phase on the multiphase turbulence equations is considered. In the present work, the model proposed by Sato et al.²²⁾ has been used to take the turbulence induced by the movement of bubbles into account. The expression is:   

$${\mu }_{t,p}=C_{\mu ,p}\rho {\alpha }_{p}L\left|{\overrightarrow{u}}_p-{\overrightarrow{u}}_q\right|$$(7)

where L is the diameter of the bubbles and C_μ,p=0.6.

2.2. Tomiyama Lift Force Model

The lift force is known to be responsible for the segregation of small and larger bubbles in bubbly flows²⁸⁾ as that in the gas-string ladle. In the gas-stirring ladle for steel refining, the liquid is in batch and the gas is passed continuously from the bottom plug with deformation. The mean axial liquid velocity decreases from the bubble plume to the wall and induces bubble rotation. The Magnus force arising out of this rotation along with the axial slip velocity helps the bubble move toward the wall. So the direction of lift force acting on the dispersed gas phase is opposite to the direction of the lift force in Eq. (4) where C_l is a positive value and although considerable research efforts have been made, there are not generally accepted values or correlations to be used.^29,30,31) According to the correlations proposed by Tomiyama et al.,^32,33) the sign of C_l depends on the bubble size and the model proposed by them is applicable to the lift force on larger-scale deformation bubbles in the ellipsoidal and spherical cap regimes. In this work, the bubble size distribution is predicted by the population balance model and the influence of bubble size distribution on the lift force is also taken into consideration. The Tomiyama lift force model lightly modified by Frank et al.³⁴⁾ is given by:   

$$C_l=\hspace{0.17em}\{\begin{array}{l}min\left[0.288tanh(0.121\mathrm{R}{\mathrm{e}}_p),f\left(E{o}^{\prime }\right)\right]\hfill \\ f\left(E{o}^{\prime }\right)\hfill \\ -0.27\hfill \end{array}\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\begin{array}{c}E{o}^{\prime }\le 4\\ 4<E{o}^{\prime }\le 10\\ 10<E{o}^{\prime }\end{array}$$(8)

where   

$$\begin{array}{c}f\left(E{o}^{\prime }\right)=0.00105E{o^{\prime }}^3-0.0159E{o^{\prime }}^2\\ -0.0204E{o}^{\prime }+0.474\end{array}$$(9)

Re_p is the bubble Reynolds number and $E{o}^{\prime }$ is a modified Eötvös number based on the long axis of the deformable bubble L_h:   

$$L_h=L{\left(1+0.163E{o}^{0.757}\right)}^{1/3}$$(10)

where Eo=g(ρ_q−ρ_p)L²/σ and $E{o}^{\prime }$ =g(ρ_q−ρ_p) $L_h^2$ /σ. In this model, for small air bubbles in water (L<4 mm), 0<C_l ≤ 0.288 and C_l becomes negative when L becomes larger. As the lift force is only significant for phases that separate quickly, only the lift force between water and gas bubble is included.

2.3. Discrete Population Balance Model

The Population Balance Model (PBM) was employed and the discrete method was used by discretizing the bubble population into a finite number of size intervals. Assuming the gas phase to be made of spherical bubbles of diameter L, the dispersed phase volume fraction can be expressed as:   

$${\alpha }_p(\overrightarrow{x},t)={\int }_0^{+\infty }n\left(\overrightarrow{x},\mathrm{\Phi },t\right)\frac{\pi }{6}{L}^{3}dL$$(11)

The bubble state vector is characterized by a set of external coordinates ( $\overrightarrow{x}$ ) which denote the spatial position of the bubble and (Φ) is referred to as an internal coordinate. n( $\overrightarrow{x}$ ,Φ,t) is the number density function. Assuming that Φ is the bubble volume, the transport equation for n( $\overrightarrow{x}$ ,Φ,t) can be given as:   

$$\begin{array}{l}\frac{\partial }{\partial t}\left[n\left(V,t\right)\right]+\nabla \cdot \left[\overrightarrow{u}n\left(V,t\right)\right]=\\ B_{ag}\left(V,t\right)+B_{br}\left(V,t\right)-D_{ag}\left(V,t\right)-D_{br}\left(V,t\right)\end{array}$$(12)

2.3.1. Luo Aggregation Kernel

The birth and death terms due to aggregation B_ag and D_ag in Eq. (12) can be written as:   

$$B_{ag}=\frac{1}{2}{\int }_0^V{\mathrm{\Omega }}_{ag}\left(V-V^{\prime },V^{\prime }\right)n\left(V-V^{\prime },t\right)n\left(V^{\prime },t\right)d{V}^{\prime }$$(13)

  

$$D_{ag}={\int }_0^{\infty }{\mathrm{\Omega }}_{ag}\left(V,V^{\prime }\right)n\left(V,t\right)n\left(V^{\prime },t\right)d{V}^{\prime }$$(14)

where Ω_ag is the aggregation kernel with the unite of m³/s, the Luo model²⁶⁾ defined the aggregation kernel as the rate of bubble volume formation as a result of binary collision of bubbles with volumes V_i and V_j:   

$${\mathrm{\Omega }}_{ag}\left(V_i,V_j\right)={\omega }_{ag}\left(V_i,V_j\right)P_{ag}\left(V_i,V_j\right)$$(15)

where ω_ag is the frequency of collision and is defined as:   

$${\omega }_{ag}\left(V_i,V_j\right)=\frac{\pi }{4}{\left(L_i+L_j\right)}^2{n}_i{n}_j{\overline{u}}_{ij}$$(16)

where ${\overline{u}}_{ij}$ is the characteristic velocity of collision of two bubbles with diameter L_i and L_j. The expression is as follows:   

$${\overline{u}}_{ij}=1.43{\epsilon }^{1/3}{\left(L_i^{2/3}+L_j^{2/3}\right)}^{1/2}$$(17)

where ε is the turbulence dissipation.

P_ag is the probability that the collision results in coalescence. It can be written as:   

$$P_{ag}=exp\left\{-c_1\frac{{\left[0.75\left(1+L_i^2/L_j^2\right)\left(1+L_i^3/L_j^3\right)\right]}^2}{{\left({\rho }_p/{\rho }_q+0.5\right)}^{1/2}{\left(1+L_i/L_j\right)}^3}W{e}^{1/2}\right\}$$(18)

where c₁ is a constant of order unity, L_i/L_j is the bubble size ratio and ρ_p/ρ_q is the ratio of bubble density and liquid density. We is the Weber number defined as:   

$$We=\frac{{\rho }_q{L}_i{\left({\overline{u}}_{ij}\right)}^2}{\sigma }$$(19)

where σ is the interfacial tension.

2.3.2. Luo and Lehr Breakage Kernel

One of the main hypotheses is the binary nature of the phenomenon. Two bubbles coalesce to form a bigger bubble and a single bubble breaks up to form two smaller bubbles. The birth rate of bubbles of volume V due to breakage is:   

$$B_{br}={\int }_{{\mathrm{\Omega }}_V}\lambda g^{\prime }\left(V^{\prime }\right)\beta \left(V|V^{\prime }\right)n(V^{\prime },t)d{V}^{\prime }$$(20)

where g′(V′) is the fraction of bubbles of volume V′ breaking per unit time (breakage frequency) and β(V|V′) is the probability density function (PDF, daughter size distribution) of bubbles breaking from volume V′ to V. From the expression, g′(V′)n(V′)dV′ bubbles of volume V′ break per unit time, producing λg′(V′)n(V′)dV′ bubbles where λ is the number of child bubbles produced per parent bubble (two for binary breakage).

The death rate of volume V due to breakage is given by:   

$$D_{br}=g^{\prime }\left(V\right)n\left(V,t\right)$$(21)

In present work, the Luo and Lehr breakage model^26,35) is used encompassing both breakage frequency and the PDF of breaking bubbles. The general breakage rate per unit volume is usually written as:   

$${\mathrm{\Omega }}_{br}\left(V,V^{\prime }\right)=g^{\prime }\left(V^{\prime }\right)\beta \left(V|V^{\prime }\right)\left[m^{-3}{s}^{-1}\right]$$(22)

where the parent bubble has a volume V′ and the daughter particle has a volume V.

The occurrence of breakage is determined by the energy level of the arriving eddy and only eddies of length scale smaller than the bubble diameter can induce bubble oscillations. The general form is the integral over the size of eddies l hitting the bubble with diameter L. The integral is taken over the dimensionless eddy size ξ=l/L. The general form is as follows:   

$$\begin{array}{l}{\mathrm{\Omega }}_{br}\left(V,V^{\prime }\right)=\\ 0.923\alpha {\left(\frac{\epsilon }{L^2}\right)}^{1/3}{\int }_{l_{min}/L}^1\frac{{\left(1+\xi \right)}^2}{{\xi }^{11/3}}exp\left(-b{\xi }^{-11/3}\right)d\xi \end{array}$$(23)

where b is modeled as:   

$$b=12\frac{\left[f_v{}^{2/3}+{(1-f_v)}^{2/3}-1\right]\sigma }{2.047\rho {\epsilon }^{2/3}{L}^{5/3}}$$(24)

where f_v is the breakage volume fraction defined as the ratio of one of the daughter bubble volumes and the parent volume. For binary breakage, the breakage kernel must be symmetrical with respect to f_v=0.5.

2.3.3. Numerical Method for PBM

As mentioned above, the approach used in present work is the discrete method developed by Hounslow et al.,²³⁾ Litster et al.²⁴⁾ and Kumar et al.²⁵⁾ Aggregation and breakage involve the migration of bubbles from a given class to the neighboring classes. The PBE is written in terms of volume fraction of bubble size i:   

$$\begin{array}{l}\frac{\partial ({\rho }_p{\alpha }_i)}{\partial t}+\nabla \cdot \left({\rho }_p{\overrightarrow{u}}_i{\alpha }_i\right)=\\ {\rho }_p{V}_i\left(B_{ag,i}-D_{ag,i}+B_{br,i}-D_{br,i}\right)\end{array}$$(25)

where α_i is the volume fraction of the bubble size i and it is defined as:   

$${\alpha }_i=N_i{V}_i\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }i=0,1,\cdot \cdot \cdot ,N-1$$(26)

where V_i is the volume of the bubble size i and N_i is defined as:   

$$N_i(t)={\int }_{V_i}^{V_{i+1}}n\left(V,t\right)dV$$(27)

The breakage formulation for the discrete method in this work is based on the Hagesather method.³⁶⁾ And it is assumed that there is no breakage for the smallest bubble class.

Finally, for coupling with the fluid dynamics, the Sauter mean diameter of the bubbles is calculated in every time step and returned to the functions including bubble diameter. The Sauter mean diameter can be calculated as follows:   

$$d_{32}=\frac{\sum N_i{L}_i^3}{\sum N_i{L}_i^2}$$(28)

Thus a complete two-way coupling model between CFD and PBM is assured.

3. Experimental Details

In order to measure the bubble distribution inside the ladle and validate the mathematical model, a one-third scale water model was established. The entire experimental system is shown in Fig. 2. The N₂ (25°C, 1 atm) was chosen to simulate argon gas in the real ladle, and the gas was injected into the liquid from a nozzle made of the porous mullite. The mullite used in the experiment is the same as the real ladle. And the radial position of nozzle could be changed with an adjustable bottom as shown in Fig. 3(a). Water and oil were used to simulate the molten steel and the top slag layer respectively. The momentum of the water was obtained from the injection and floatation of gas bubbles. The details about the geometric parameter and the material properties are shown in Table 1.

Fig. 2.

The experimental system.

Fig. 3.

The (a) bottom and plug of water model (b) mesh and boundary conditions.

Table 1. Parameters of both experimental and numerical simulation.

Parameters	Experimental and numerical model
Bottom diameter	617 mm
Slope angle	2.44 degree
Liquid depth	700 mm
Slag layer thickness	50 mm
Porous plug diameter	43.4 mm
Plug radial position	0 R–0.73 R
Water density	1000 kg·m−3
Water viscosity	0.001 kg·m−1·s−1
Oil density	900 kg·m−3
Oil viscosity	0.006 kg·m−1·s−1
Gas density	1.138 kg·m−3 (25°C)
Gas viscosity	1.663×10−5 kg·m−1·s−1
Gas flow rate	70–130 NL/h
Interfacial tension	0.072 N/m (water/air)

As shown in Fig. 2, the mirror and a camera were used to capture the slag open eye and a high speed camera was used to capture the bubble distribution when the flow reached a steady state. Figure 4 shows the bubble aggregation phenomenon and the bubble size analysis process of the open software ImageJ, which can distinguish the chosen area in picture between the color gamut values. Firstly, the scale should be set, with the help of the rule set in ladle, to tell the real size of picture. Then the bubbles can be marked according to the color difference. Finally, the area of every bubble can be reported and the mean bubble diameter can be calculated. So that the bubble diameter predicted by the numerical model can be compared to the experimental bubble diameter.

Fig. 4.

(a) Bubble aggregation phenomenon (b) bubble size analysis process.

The saturated brine was added to the liquid from the slag eye when flow reaches to be steady and two monitor points (as shown in Fig. 2) were placed in both sides of the plug in the lower zone of the ladle called the dead zone. Two conductivity meters connected to computers were used to measure the curves of the conductivity change caused by the solute transport. Then the measured conductivity quantities were changed to non-dimensional quantities as well as the numerical calculation quantities so that the non-dimensional quantities of both experimental measurement and numerical calculation can be compared.

4. Numerical Details

To verify the computational model, the model size for the computation was set the same to the water model as shown in Table 1. All computations were performed using the commercial CFD software FLUENT. The grid densities of the 3 D mesh were determined as follows: the maximum mesh size of 12 mm was used, and a stretching ratio of 1.1 was used near the wall and the slag layer resulting in an approximate grid size of 4 mm and the mesh size of the inlet was also set to 4 mm. One of the meshes and the boundary conditions are shown in Fig. 3(b). All the meshes for different plug positions were set to the same grid density to avoid the influence of mesh size. The number of cells was between 250000 to 300000 for the various cases. All the cases were carried out by transient simulations and the double-precision solver was used.

As the flow reaches a steady state, a user defined scalar (UDS) was added to simulate the solute transport. It is assumed that the UDS only transports in water and the equation in the multiphase flow can be written as:   

$$\frac{\partial {\alpha }_q{\rho }_q\chi }{\partial t}+\nabla \cdot \left({\alpha }_q{\rho }_q{\overrightarrow{u}}_q\chi -{\alpha }_q{\mathrm{\Gamma }}_q\nabla \chi \right)=0$$(29)

where χ is the value of UDS and Γ_q is the diffusion coefficient fixed by experiment in the present work. The UDS was added at the position where the slag eye formed as the experiment did. Meanwhile the UDS variable was monitored by 5 monitor points and the mixing time was determined as the time when the gaps between the non-dimensional quantities were less that 3%.

In the population balance model, the number of bubble bins was specified 6 for balancing the computational cost and the accuracy for predicting the bubble size distribution. The range of bubble size was determined by the analysis of the experimental results. And it was assumed that the bubble diameter in the initial state was always the smallest. The conservation equations were discretized using the control volume technique. The phase coupled SIMPLE scheme was used for the pressure-velocity coupling of the three-phase flow and the hybrid-upwind discretization scheme was used for the convective terms. A physical time scale of 0.01 s was adopted for all the unsteady simulations. All simulations were stopped when flows became steady and the total time was 100 s.

5. Results and Discussion

5.1. Model Verification

The computation parameters for model verification are shown in Table 2. It is assumed that the smallest bubble diameter in the computational domain and the initial bubble diameter is 1 mm. In the experiment, the bubbles floating up to the top surface were captured by a high speed camera and the slag open eye was also captured by a camera through a mirror. Figure 5 shows the bubble volume fraction and size distribution and the slag layer of both experimental and numerical result. The numerical result was predicted using the Tomiyama lift force model. It is found that as bubble floating upward the top surface, most bubbles do not float straightly with the influence of lift force found in this work. And the influence on larger bubbles is more significant resulting in the diameter distribution in the upper part where the bubble diameter is smaller in the center of bubble plume.

Table 2. Computation parameters for model verification.

Parameters	Value
Gas flow rate	90 NL/h
Slag layer thickness	50 mm
Liquid depth	700 mm
Range of bubble diameter	1.0–3.1 mm
Initial Bubble size	1.0 mm
Number of bins	6

Fig. 5.

Comparison between calculation and experiment.

The curves in Fig. 6(a) show the bubble diameter distribution of six different positions on height direction and the points show the bubble diameter measured through the picture captured by high speed camera. It is found that, in upper zone, the bubble diameter in the center of bubble plume is smaller as the larger lift force on larger bubbles. As well as shown in Fig. 6(b), the curves predicted by computation fit well with the profile of experimental bubble plume. Figure 7 shows both the change of gas volume fraction and bubble Sauter mean diameter with the height of ladle. The values of the curves are on the location of the centerline and the experimental results are obtained from analyzing the bubble mean diameter at the positions along with the centerline as shown. Compare the bubble size distribution of computation result to experiment result, it is found that bubble diameter predicted fit well with the experiment. Moreover, it is found that bubble aggregation rate is much higher and gas volume fraction decreases quickly in the zone where the height is less than 0.1 m. It may cause the inclusion removal rate decrement because smaller bubble is good for inclusion removal.

Fig. 6.

Comparison of bubble diameter distribution between experiment and calculation.

Fig. 7.

The change of gas volume fraction and bubble diameter with height during bubble floatation.

5.2. Effect of Gas Flow Rate

The gas flow rate is a very important parameter to homogenize the melt and temperature. Figure 8 shows the slag open eye both predicted by computation and captured by experiment with four different gas flow rates: 70, 90, 110 and 130 NL/h. The result shows that when the gas flow rate is less than 90 NL/h, the slag eye is not steady. And then the slag eye becomes steady when gas flow rate is larger. Moreover, the slag eye increases with the increment of gas flow rate.

Fig. 8.

The slag open eye of different gas flow rate: (a) 70 NL/h (b) 90 NL/h (c) 110 NL/h (d) 130 NL/h.

The Weber number can be used to study the drop break-up rate in turbulent dispersed flows.^37,38) The liquid mixed with bubbles flows to the top surface and breaks up the top slag layer. Then bubbles break into the air while water flows back down. This will improve the interfacial reaction rate of the melt but meanwhile will cause the slag entrapment. Large slag drops will quickly float up back to the slag but small drops will get deep down into the ladle which is bad to the steel quality. The slag entrapment mechanism is the same to the drop break-up so the present work also uses the Weber number to study the slag entrapment. As shown in Eqs. (17) and (19), the Weber number is in direct proportion to ε^2/3. So the turbulence dissipation rate is used to study the slag entrapment in this work. Figure 9 shows the turbulence dissipation rate distribution and the liquid velocity in the top slag layer. From the velocity plotted on the slag layer, it can be found that water flows down back after arrived at the top surface, which results in the velocity at the center towards the wall and the slag flows to the center. The high turbulence dissipation rate at the center position of the slag layer will lead to forming the slag open eye and also re-oxidation. And the turbulence dissipation rate, which can predict slag break-up probability, increases with the increment of gas flow rate. A suitable gas flow rate is needed to balance the desulphurization efficiency and the cleanness of the steel during ladle operation and the gas critical flow rate (70 NL/h in present water model) is considered as an important operating parameter.

Fig. 9.

The turbulence dissipation rate and velocity in the slag layer of gas flow rate (a) 70 NL/h (b) 90 NL/h (c) 110 NL/h (d) 130 NL/h.

Figure 10 shows the bubble diameter change in the vertical direction with different gas flow rates. It is found that gas bubble diameter distribution in upper zone increases with increasing the gas flow rate. But the diameter distribution in lower zone shows a different change law that bubble Sauter mean diameter in the zone of y<0.4 m decreases when gas flow rate increases from 110 NL/h to 130 NL/h. As gas flow rate increases, both the momentum and number density of gas bubble increase, and when the gas flow rate increases from 70 NL/h to 110 NL/h, the increment of bubble collision and aggregation rate because of the bubble number density increment plays a more significant role. On the contrary, when the flow rate increases from 110 NL/h to 130 NL/h, the decrement of bubble collision and aggregation rate because of the bubble momentum increment plays a more significant role.

Fig. 10.

The bubble diameter distribution of different gas flow rates.

5.3. Effect of Plug Radial Position

Mixing is very important in the steel-making refining processes, and the intensity of mixing is usually defined by mixing time. Many studies^{1,2,3,4,5,6,7,8,9,10,11,12,13,14)} have been taken to understand the mixing characteristics in gas-stirring systems and a great attention has been paid to study the effect of nozzle position and number on the mixing in the gas- stirred systems. The novelty and originality of this work are as follows: 1) Both the mixing efficiency and the erosion of ladle wall are taken into consideration to find a suitable plug radial position, while previous works always ignore the ladle wall erosion. 2) The effect of bubble diameter distribution and the slag layer on the flow field are both taken into consideration. 3) The study of plug radial position in this study is more detailed and a flow pattern change is found when the position changes from 0.67 R to 0.73 R. In this work, the plug radial positions including 0.5 R, 0.56 R, 0.62 R, 0.67 R and 0.73 R are studied because the real ladle studied in this work is 0.67 R and we try to ensure the equal interval.

Figure 11 shows the simulation result of liquid flow field, gas distribution and slag transporting with centric and eccentric blowing. The result shows that the liquid flow in ladle with eccentric blowing is different from that with centric blowing. With centric blowing, a small circulation flow on the upper zone near the wall forms and will entrap the slag into the circulation flow zone. The white line shows the liquid flow formed due to the liquid flowing to the two sides of the slag eye. It will also cause a circulation zone under the slag layer. Figure 12 shows the numerical result of liquid flow field with four plug radial positions of 0 R, 0.56 R, 0.67 R and 0.73 R. It can be found that with the increment of plug radial distance, the circulation flow near the wall becomes weaker. As the radial position reaches at 0.73 R, the flow pattern changes. Figure 13 shows the turbulence dissipation rate of the same four plug positions to study the slag layer break-up probability. And the result shows that with increasing the plug radial distance, the turbulence dissipation rate increases which shows the higher break-up probability.

Fig. 11.

The flow field of liquid and the stream lines of slag with (a) centric (b) eccentric blowing.

Fig. 12.

Numerical result of liquid flow field with different plug radial positions of (a) 0 R (b) 0.56 R (c) 0.67 R (d) 0.73 R.

Fig. 13.

Turbulence dissipation rate with different plug radial positions of (a) 0 R (b) 0.56 R (c) 0.67 R (d) 0.73 R.

Figure 14(a) shows the sample points and the coordinate system of simulation with the slag layer shown by an iso-surface. The origin is located at the center of the bottom and the coordinates of the sample points are: 1(0.29, 0.015, 0), 2(−0.29, 0.015, 0), 3(0, 0.35, 0), 4(0, 0.65, −0.29), 5(−0.29, 0.65, 0). Figure 14(b) shows the two monitor points located at the same location to the numerical model. The two points are located in the widely acknowledged dead zone and can avoid the influence of plug position change. Also an examination from sample points shows that the mixing time is mainly retarded by the sample points of Point 1 and Point 2.

Fig. 14.

The monitor points in (a) numerical model (b) water model.

Figure 15 shows the experimental and numerical non-dimensional variable change of Point 1 and Point 2 in the case using the diffusion coefficient Γ_q=2. The UDS is added at the time that flow reaches a steady state (t=100 s) and the time that shown in the figure is recorded after that. The non-dimensional variable is described as χ/χ_∞ where χ_∞ is a steady value of the quantity obtained from both experiment and numerical calculation. The points are the measured average values of five experiments. The result shows that the UDS transport equation fixed can predict an accurate solute transport in liquid. Figure 16 shows the UDS change of the five sample points with time passing and, in this study, the mixing time is defined as the time of attaining a 97% degree of homogenization. As shown in Fig. 16, the mixing time with the plug radial position of 0.67 R, which was the same to the real ladle studied in this work, is 99.8 s. It is also found that the sample points of Point 1 and Point 2 are the hardest mixing places.

Fig. 15.

The computational and experimental non-dimensional variable of (a) point 1 (b) point 2.

Fig. 16.

The UDS change with physical time of 0.67 R radial plug position.

In order to find a suitable plug radial position, the same method was taken to calculate the mixing time of plug radial positions of 0.5 R, 0.56 R, 0.62 R, 0.67 R and 0.73 R. Meanwhile, the wall shear stress is also studied to predict the erosion of ladle lining. Figure 17 shows the wall shear stress distribution of different plug radial positions. Large wall shear stress will not only influence the ladle life-span but also the steel quality. The maximum wall shear stress was considered as an important parameter to predict the erosion and increases with the increasing the plug radial distance from 0.56 R to 0.67 R. Especially, the maximum wall shear stress of plug radial position at 0.73 R is lower than 0.67 R and higher to 0.62 R because of the liquid field change that the near wall circulation flow of 0.73 R becomes weak. Figure 18 shows both the mixing time and maximum wall shear stress of different plug radial positions. It is found that the mixing time decreases with increasing the plug radial distance and the maximum wall shear stress of radial position at 0.67 R is the largest.

Fig. 17.

Wall shear with different plug radial positions of (a) 0.56 R (b) 0.62 R (c) 0.67 R (d) 0.73 R.

Fig. 18.

The mixing time and maximum wall shear stress of different plug radial positions.

6. Conclusions

In this study, a relatively comprehensive Eulerian three-phase model coupled with the PBM was established and predictions of the gas volume fraction, bubble size distribution, slag open eye, slag break-up probability, wall shear stress and solute transporting were realized. The RNG k-ε turbulence model was taken to simulate the complex fluid flow with both high-Reynolds number and low-Reynolds number and with a rapidly strained flow under the slag layer. The comparison of bubble size distribution, slag open eye size and solute transporting between experiment and numerical calculation was made, and the influence of gas flow rate and plug radial position was studied. The conclusions can be drawn as:

(1) The phase interactions of Tomiyama lift force and Schiller-Naumann drag force are included in the numerical model and coupled with the PBM to take the bubble diameter effect into account. The result shows that bubbles aggregate quickly in the lower zone of ladle. After that, they tend to float separately because the velocity gradient in flow field helps the bubble move from high velocity field to low velocity field and the effect is more significant on larger bubbles.

(2) The slag open eye predicted by the numerical model agreed well with experiment result. And the critical gas flow rate to form a steady open eye was found (70 NL/h in the present water model condition).

(3) The slag break-up probability was studied using the Weber number and it was found that slag break-up probability increases with the increment of the gas flow rate and the radial distance.

(4) The mixing time and wall shear stress of different plug radial positions were studied and it was found that the mixing time decreases with increasing the plug radial distance. And the plug radial position at 0.67 R has a maximum wall shear stress as the position at 0.73 R forms a different flow field and decreases the wall shear stress.

Acknowledgement

Authors are grateful to the National Natural Science Foundation of China for support of this research, Grant No. 51210007.

Nomenclature

A_i: interfacial area, m²

B_ag, B_br: birth term due to aggregation and breakage

c, C: model constant

D_ag, D_br: death term due to aggregation and breakage

d₃₂: Sauter mean diameter, m

Eo: Eötvös number

Eo′: modified Eötvös number

${\overrightarrow{F}}_{drag}$ : drag force, N/m³

${\overrightarrow{F}}_{lift}$ : lift force, N/m³

${\overrightarrow{F}}_{VM}$ : virtual mass force, N/m³

f: drag function

f_v: breakage volume fraction=0.5

g, $\overrightarrow{g}$ : gravitational constant/acceleration, m/s²

g′: breakage frequency

k: turbulence kinetic energy, m² s⁻²

l: eddy size, m

L: bubble diameter, m

N_i: number of bubble size i

n: number density, m⁻³

p: pressure, Pa

P_ag: collision probability

Re: Reynolds number

t: time, s

$\overrightarrow{u}$ : velocity, m s⁻¹

V: volume of bubble(s), m³

We: Weber number

Greek letters

α: volume fraction

ρ: density, kg m⁻³

$\overline{\tau }$ : stress-strain tensor

μ: viscosity of q^th phase, kg m⁻¹ s⁻¹

μ_t: turbulent viscosity, kg m⁻¹ s⁻¹

ε: turbulence dissipation rate, m² s⁻³

σ: interfacial tension, N m⁻¹

Ω_ag: aggregation kernel, m³ s⁻¹

Ω_br: breakage rate per unit volume, m⁻³ s⁻¹

ω_ag: frequency of collision

λ: number of child bubbles per parent bubble

β: probability density function(PDF)

ξ: dimensionless eddy size=l/L

χ: user defined scalar (UDS)

Γ: UDS diffusion coefficient

Subscripts

q, p: continuous and dispersed phase

i, j: bubble class

m: mixture

References

• 1)   J. W.  Han,  S. H.  Heo,  D. H.  Kam,  B. D.  You,  J. J.  Pak and  H. S.  Song: ISIJ Int., 41 (2001), 1165.
• 2)   Y.  Xie and  F.  Oeters: Steel Res., 63 (1992), 93.
• 3)   K.  Yonezawa and  K.  Schwerdtfeger: Metall. Meter. Trans. B, 30B (1999), 411.
• 4)   D.  Mazumdar and  J. W.  Evans: Metall. Meter. Trans. B, 35B (2004), 400.
• 5)   M.  Peranandhanthan and  D.  Mazumdar: ISIJ Int., 50 (2010), 1622.
• 6)   K.  Krishnapisharody and  G. A.  Irons: ISIJ Int., 48 (2008), 1807.
• 7)   S. P.  Patil,  D.  Satish,  M.  Peranandhanathan and  D.  Mazumdar: ISIJ Int., 50 (2010), 1117.
• 8)   H.  Liu,  Z.  Qi and  M.  Xu: Steel Res. Int., 82 (2011), 440.
• 9)   D.  Guo and  G. A.  Irons: Metall. Meter. Trans. B, 31B (2000), 1457.
• 10)   T.  Wang,  J.  Wang and  Y.  Jin: Ind. Eng. Chem. Res., 44 (2005), 7540.
• 11)   S.  Ganguly and  S.  Chakraborty: ISIJ Int., 44 (2004), 537.
• 12)   M.  Madan,  D.  Satish and  D.  Mazumdar: ISIJ Int., 45 (2005), 677.
• 13)   B.  Li,  H.  Yin,  C.  Zhou and  F.  Tsukihashi: ISIJ Int., 48 (2008), 1704.
• 14)   S. W. P.  Cloete,  J. J.  Eksteen and  S. M.  Bradshaw: Pro. Comput. Fluid Dy., 9 (2009), 345.
• 15)   Z.  Liu,  B.  Li,  M.  Jiang and  F.  Tsukihashi: ISIJ Int., 54 (2014), 1314.
• 16)   Z.  Liu,  B.  Li,  M.  Jiang and  F.  Tsukihashi: ISIJ Int., 53 (2013), 484.
• 17)   Z.  Liu,  L.  Li,  F.  Qi,  B.  Li,  M.  Jiang and  F.  Tsukihashi: Metall. Meter. Trans. B, 46 (2015), 406. DOI: 10.1007/s11663-014-0192-y.
• 18)   H.  Hulburt and  S.  Katz: Chem. Eng. Sci., 19 (1964), 555.
• 19)   A.  Mukhopadhyay,  E. W.  Grald,  K.  Dhanasekharan,  S.  Sarkar and  J.  Sanyal: Steel Res. Int., 76 (2005), 22.
• 20)   V.  Yakhot and  S. A.  Orszag: J. Sci. Comput., 1 (1986), 3.
• 21)   V.  Yakhot,  S. A.  Orszag,  S.  Thangam,  T. B.  Gatski and  C. G.  Speziale: Phys. Fluids A, 4 (1992), 1510.
• 22)   Y.  Sato and  K.  Sekoguchi: Int. J. Multiphase Flow, 2 (1979), 79.
• 23)   M. J.  Hounslow,  R. L.  Ryall and  V. R.  Marshall: AIChE J., 34 (1988), 1821.
• 24)   J. D.  Litster,  D. J.  Smit and  M. J.  Hounslow: AIChE J., 41 (1995), 591.
• 25)   S.  Kumar and  D.  Ramkrishna: Chem. Eng. Sci., 51 (1996), 1311.
• 26)   H.  Luo and  H. F.  Svendsen: AIChE J., 42 (1996), 1225.
• 27)   L.  Schiller and  Z.  Naumann: Z. Ver. Deutsch. Ing., 77 (1935), 318.
• 28)   M. E.  Díaz,  F. J.  Montes and  M. A.  Galán: Chem. Eng. Sci., 64 (2009), 930.
• 29)   D. A.  Drew and  R. T.  Lahey: Int. J. Multiphase Flow, 13 (1987), 113.
• 30)   I.  Zun: Int. J. Multiphase Flow, 6 (1980), 583.
• 31)   M. T.  Dhotre and  B. L.  Smith: Chem. Eng. Sci., 62 (2007), 6615.
• 32)   A.  Tomiyama and  N.  Shimada: J. Press. Vesse. Technol., 123 (2001), 510.
• 33)   A.  Tomiyama,  H.  Tamai,  I.  Zun and  S.  Hosokawa: Chem. Eng. Sci., 57 (2002), 1849.
• 34)   Th.  Frank,  P. J.  Zwart,  E.  Krepper,  H.-M.  Prasser and  D.  Lucas: Nucl. Eng. Des., 238 (2008), 647.
• 35)   F.  Lehr,  M.  Millies and  D.  Mewes: AIChE J., 8 (2002), 2426.
• 36)   L.  Hagesaether,  H. A.  Jakobsen and  H. F.  Svendsen: Chem. Eng. Sci., 57 (2002), 3251.
• 37)   S.  Galinat,  O.  Masbernat,  P.  Guiraud,  C.  Dalmazzone and  C.  Noïk: Chem. Eng. Sci., 60 (2005), 6511.
• 38)   P.  Angeli and  G. F.  Hewitt: Chem. Eng. Sci., 55 (2000), 3133.