Numerical Simulation of Impinging Gas Jet on a Liquid Bath Using SPH Method

Abstract

As the basis of injection metallurgy, an impinging gas jet on a liquid bath surface was simulated and experimentally verified. The SPH simulation of the gas-liquid two phase flow was developed by using the XSPH method, and the calculation speed was considerably increased with the use of general-purpose computing on graphics processing units. For N₂ gas – water bath system, the simulated shapes and penetration length of the cavity were in good agreement with the experimental results, and the three modes of cavity were reproduced by the calculation. N₂ gas–molten iron bath system was also simulated. The cavity mode was “dimpling mode” for all calculations.

1. Introduction

The impinging gas jet on a liquid bath surface is essential for various metallurgical reactors, especially those used in the basic oxygen furnace (BOF) steelmaking process. Many studies have been reported cold model experiments to address experimental difficulties caused by high temperatures in the practical processes. However, these model experiments cannot reproduce the effects of large surface tension of molten metal and of high temperatures. Because recent developments in computational fluid dynamics (CFD) enable the direct simulation of the gas–liquid two-phase flows, the numerical simulation may replace the difficult and expensive high-temperature experiment. Recently, the Volume of Fluid (VOF) model¹⁾ was applied to an impinging gas jet on a water bath,^2,3,4,5) low-melting point alloy bath,⁴⁾ and molten steel bath.⁵⁾ There was good agreement with the cold model experiments. This model treats a gas–liquid two-phase flow as a homogeneous phase, and the flow is characterized by the liquid fraction φ (0 ≤ φ ≤ 1), where φ = 1 implies a liquid phase, and φ = 0 implies a gas phase. The region between 0 to 1 corresponds to the two-phase mixture. Therefore, the gas–liquid interface should be physically diffused. The suitable value of φ must be given to obtain a sharp interface, which involves physical ambiguity and unavoidable arbitrariness.

Because the smoothed particle hydrodynamics (SPH) method is particle-based, the Navier–Stokes equations for gas and liquid are simultaneously solved. The sharp interface can be obtained easily by the distinction between the liquid and solid particles. However, the application of the SPH method to the impinging gas jet problem is limited because of the following two reasons: 1) The velocity of a gas particle colliding with a liquid particle increases considerably because of the large mass difference. 2) The large number of particles required to fill the space results in a significant increase in calculation time. In this study, the first problem was addressed by using the XSPH method proposed by Monaghan,^6,7,8) and the second by applying a general-purpose computing on graphics processing units (GPGPU) system for the calculations.

2. Simulation Model

2.1. SPH Method for a Gas-Liquid Two Phase Flow

The three-dimensional SPH model⁹⁾ was used for the simulation in this study. The quantic spline function shown in Eq. (1) was used as the kernel function, W(r−r′, h), where r, rʹ, and h are the position vector of the particle, that of a neighboring particle, and the effective radius. q is the dimensionless inter-particle distance defined as |r−r′|/h.   

$$\begin{array}{l}W\left(\mathit{r}-{\mathit{r}}^{\prime },h\right)=\\ {\displaystyle \frac{3}{359\pi h^3}}\{\begin{array}{r}\hfill \left[{\left(3-q\right)}^5-6{\left(2-q\right)}^5+15{\left(1-q\right)}^5\right]\hspace{1em}for\mathrm{\hspace{0.17em} } 0\le q\le 1\\ \hfill \left[{\left(3-q\right)}^5-6{\left(2-q\right)}^5\right]\hspace{1em}for\mathrm{\hspace{0.17em} } 1<q\le 2\\ \hfill {\left(3-q\right)}^5\hspace{1em}for\mathrm{\hspace{0.17em} } 2<q\le 3\\ \hfill 0\hspace{1em}for\mathrm{\hspace{0.17em} } 3<q\hspace{1em}\hspace{0.5em}\mathrm{\hspace{0.17em} }\end{array}\end{array}$$(1)

The governing equation is the Navier–Stokes equation:   

$$\frac{D\mathit{u}}{Dt}=-\frac{1}{\rho }\nabla P+\frac{\mu }{\rho }\mathrm{\Delta }\mathit{u}+\frac{\mathit{f}}{\rho }$$(2)

where u is the velocity, ρ is the density, P is the pressure, μ is the viscosity, and f is an external force. The surface tension was calculated by the pairwise potential model⁹⁾ and substituted into f.

Usually, the density of the i-th particle ρ_i is calculated using Eq. (3) in SPH calculations; however, for a gas–liquid two-phase flow, calculation errors arise because the density of a liquid particle is thousand times higher than that of a gas particle. In this study, Eq. (4)¹⁰⁾ was used instead of Eq. (3) to prevent errors in the calculation, where m_i and r_i are the mass and position vectors of the i-th particle.   

$${\rho }_i={\sum }_j{m}_{j}W\left({\mathit{r}}_i-{\mathit{r}}_j,h\right)$$(3)

  

$${\rho }_i=m_i{\sum }_{j}W\left({\mathit{r}}_i-{\mathit{r}}_j,h\right)$$(4)

The density difference problem in the calculation of the pressure term in Eq. (2) was addressed by using Eq. (5),¹¹⁾ where V_i is the volume of i-th particle and Π_ij is an artificial viscosity introduced by Monagahan⁶⁾ for numerical stability.   

$$\begin{array}{l}-{\displaystyle \frac{1}{\rho }}\nabla P_i=\\ -\left\{{\displaystyle \frac{1}{m_{\text{i}}}}{\sum }_j\left(P_i{V}_i{}^2+P_j{V}_j{}^2\right)\nabla W\left({\mathit{r}}_i-{\mathit{r}}_j,h\right)+m_{\text{i}}{\sum }_j{\mathrm{\Pi }}_{ij}\nabla W\left({\mathit{r}}_i-{\mathit{r}}_j,h\right)\right\}\end{array}$$(5)

The artificial viscosity in Eq. (5) was calculated by using Eqs. (6), (7), (8), (9).   

$${\mathrm{\Pi }}_{ij}=\{\begin{array}{r}\hfill {\displaystyle \frac{-{\alpha }_1\overline{c_{ij}}{\omega }_{ij}+{\alpha }_2{\omega }_{ij}{}^2}{\overline{{\rho }_{ij}}}} \hspace{1em}{u}_{ij}\cdot r_{ij}<0\\ \hfill 0\hspace{1em}{u}_{ij}\cdot r_{ij}\ge 0\end{array}$$(6)

  

$${\omega }_{ij}=\frac{h{\mathit{u}}_{ij}\cdot {\mathit{r}}_{ij}}{{\mathit{r}}_{ij}{}^2+{\beta }^2}$$(7)

  

$$\overline{c_{ij}}=\frac{c_i+c_j}{2}$$(8)

  

$$\overline{{\rho }_{ij}}=\frac{{\rho }_i+{\rho }_j}{2}$$(9)

In the equations, c_i, r_ij, and u_ij are the sonic velocity of i-th particle, distance between the i-th and j-th particles, and their velocity difference, respectively. The parameters α₁, α₂ and β were set to be 0.01, 0.01 and 0.1h, respectively.

The same problem in the calculation of the viscous term in Eq. (2) was also overcome by using Eq. (10),^10,12) where μ_i is the viscosity of the i-th particle.   

$$\frac{\mu }{\rho }\mathrm{\Delta }\mathit{u}=\frac{1}{m_{\text{i}}}{\sum }_j\frac{2{\mu }_i{\mu }_j}{{\mu }_i+{\mu }_j}\left(V_i{}^2+V_j{}^2\right)\frac{{\mathit{r}}_{ij}\cdot \nabla W\left({\mathit{r}}_i-{\mathit{r}}_j,h\right)}{{\mathit{r}}_{ij}{}^2+{\eta }^2}{\mathit{u}}_{ij}$$(10)

The right-hand side of this equation can be calculated without using the particle density, which maintains the symmetry of the force. η was set to be 0.01h in the calculation,

2.2. XSPH Method

The density difference problem, especially near the gas–liquid interface, causes the improbable behavior of a gas particle, which results in the serious instability of the gas–liquid interface. To address for this instability, the XSPH method was proposed by Monaghan.^6,7,8) The velocity was weight-averaged for particles in the vicinity of the interface, and the revised velocity of the i-th particle, $\widehat{{\mathit{u}}_i}$, was calculated using Eq. (11), where ε is a constant and was set to be 0.5 in this study.   

$$\widehat{{\mathit{u}}_i}={\mathit{u}}_i+\epsilon {\sum }_j\frac{2V_i{V}_j}{V_i+V_j}\nabla W\left({\mathit{r}}_i-{\mathit{r}}_j,h\right)$$(11)

2.3. Utilization of GPGPU

The GPGPU was utilized for the simulation and acceleration of computer processing in this study. After transferring the data stored -on the computer- in the CPU to the GPU, the calculations were parallel processed in the GPU for every particle at one time. The calculation space was divided into cubes with sides equal in length to the effective length, h, to minimize the neighboring particle search step. Because the SPH method employs the explicit method, the effect of parallel processing is remarkable. As a result, when the GPU (Tesla P100) was added to the CPU (Intel Xenon E5-2620 v4) system, the calculation speed increased to more than three hundred times its initial value, for the 3D dam breaking simulation with one million particles.

3. Results and Discussion

3.1. Impingement of the N2 Gas Jet on the Water Bath

The cylindrical vessel (i.d. = 100 mm, height = 80 mm) made of acrylic resin, was used for the cold model experiment. It was filled with ultra-pure water was filled to a depth of 50 mm, and N₂ gas was blown through the knife edged stainless steel nozzle (ID = 3 mm) to the center of the water surface, at a height of 10 mm, as shown in Fig. 1. N₂ gas flow rate was precisely controlled by the thermal mass flow controller. The formed cavity was observed by using the video system, and its time-averaged depth was obtained from a series of snap shots.

Fig. 1.

Comparison of calculated snapshots with experimental results. (Online version in color.)

The calculation space was the cylinder (i.d. = 100 mm, height = 80 mm) same as that used in the experiment. The walls were constructed by immobile particles; however, the slit (10 mm) was set on top of the side wall as the gas outlet. Initially, the liquid particles were filled up to a depth of 50 mm, and the remaining space was filled with gas particles. Then, the gas particles were continuously injected from the circular nozzle with a 3-mm diameter. Gas particles were released via the slit to maintain a constant pressure of the space.

Table 1 summarizes the simulation conditions.¹³⁾

Table 1. Calculation conditions for the simulation.

	N2 gas	water	Fe
Temperature (K)	298	298	1873
Density (kg/m3)13)	1.15	997	7020
Viscosity (mPa·s)13)	0.0018	0.089	4.8
Surface tension (mN/m)13)	–	73	1800
Initial inter-particle distance (m)		2.0×10−3
The number of particles	288815	301000	301000
Time step (s)		2.0×10−6

Figure 1 compares the typical snap shots of the experiments and calculations.

For the convenience of visualization, gas particles are eliminated and only liquid particles are shown in the figure. The injection velocity u₀ is the gas velocity at the nozzle exit, which was calculated from the gas flow rate in the experiment and was given to the gas particles in the simulation. The maximum value of u₀ was 20 m/s to maintain the stability of the calculation.

When u₀ = 5–10 m/s, a slight cavity and no droplet formation were observed, which corresponds to the “dimpling mode” according to Molloy’s classification¹⁴⁾ of cavity. When u₀ = 15 m/s, a shallow cavity and droplet formation were observed; however, the calculation did not show the droplet formation. The current hypothesis for this outcome is that the size of the liquid particles were comparable to actual droplet size, and this was solved by the reduction of particle size. This cavity mode appears to correspond to the “splashing mode.” The shapes of the cavity show a good agreement up to u₀ = 15 m/s.

When u₀ = 18 m/s, the cavity vertically penetrated the bath, which corresponds to “penetration mode.” Violent splashing was observed both in the experiment and simulation. The penetrated cavity oscillated in the experiment; however, the cavity tended to stay as a plume in the bath in the simulation.

The cavity depth was defined as the lowest position of the gas particle from the steady surface of the bath. Figure 2 shows the variation in the calculated cavity depth with time. Initially, the cavity penetrated deep into the bath, after which it began to oscillate. The cavity depth increased, and the oscillation period decreased with an increase in injection velocity, which agrees with the experimental observation. When u₀ = 18 m/s, the cavity submerged into the bath as a plume after 0.4 s. The precise measurement of the oscillation period in the experiment was not achieved because of the poor temporal resolution of the video system. The observations using a high-speed camera system will be the future work.

Fig. 2.

Variation of calculated cavity depth with time for N₂ gas–water system. (Online version in color.)

The time averaged cavity depth was calculated by averaging the calculation results between 0.1 s to 0.5 s, which is plotted against the injection velocity as open circles in Fig. 3. In the same figure, the solid circles indicate the observed cavity depth, and the dotted line is the reported empirical equation given by Eq. (12),¹⁵⁾ where, L, H, M, ρ_l, and g are the cavity depth, nozzle height, momentum of gas jet, density of the bath liquid, and acceleration of gravity, respectively.   

$$\frac{L}{H}{\left(1+\frac{L}{H}\right)}^2=\frac{{15.2}^2}{2\pi }\frac{M}{{\rho }_{l}g{H}^3}$$(12)

Fig. 3.

Comparison of calculated average cavity depth with experimental results. (Online version in color.)

The experimental and simulation results show a good agreement up to u₀ = 15 m/s; however, the simulation overestimates at higher u₀ values because of the submergence of the plume. Because the circulation flow of the bath may affect this phenomenon, a larger bath design would address this problem. It is thought that the difference between the calculation and experimental result for u₀ = 6–9 m/s is because of the interference of the reflected surface wave with the cavity. This interference can be overcome by the same treatment.

3.2. Impingement of N₂ Gas Jet on the Iron Bath

The same simulation program was applied to the calculation for the molten iron bath system. The calculation conditions are listed in Table 1. The behavior of the cavity position and the snapshots are shown in Fig. 4. The cavity depth increased once and returned to a constant level. It oscillated very slightly because of a large surface tension, and increased with an increase in u₀. The cavity mode was “dimpling mode” for all calculations as shown in the snapshots. The decrease in particle size and time step enabled the calculation of considerably higher gas velocity jets; however, the enhancement of the calculation performance to compensate for the increasing calculation cost is necessary.

Fig. 4.

Variation of the calculated cavity depth with time for N₂ gas–molten iron system. (Online version in color.)

4. Conclusion

The GPGPU was utilized to accelerate the SPH simulation. The calculation speed increased more than three hundred times, compared with that of the CPU alone. The SPH simulation for the gas–liquid two-phase flow was developed by using the XSPH method and applied to the impinging gas jet on a liquid bath surface. The obtained results are as follows:

(1) The simulated shapes and penetration length of the cavity were in good agreement with the experiment up to u₀ = 15 m/s.

(2) The three modes of cavity proposed by Molloy were reproduced by the simulation of an N₂ gas–water bath system.

(3) N₂ gas–molten iron bath was simulated. The cavity mode was “dimpling mode” for all calculations.

Acknowledgement

The authors are grateful to the JSPS Grant-in-Aid for scientific research programs (Grant No. 17K06881) and the Waseda University Grant for special research projects (Grant No. 2017K-180) for financial support.

References

• 1)   C. W.  Hirt and  B. D.  Nichols: J. Comput. Phys., 39 (1981), 201.
• 2)   M.  Ando,  D.  Komagata,  K.  Takahashi,  T.  Ishii and  K.  Ueno: Tetsu-to-Hagané, 101 (2015), 82 (in Japanese).
• 3)   O.  Nakamura,  T.  Yamamoto,  K.  Ueno and  K.  Takatani: Tetsu-to-Hagané, 101 (2015), 88 (in Japanese).
• 4)   S.  Sato,  M.  Ando,  J.  Okada,  Y.  Ueda and  M.  Iguchi: Tetsu-to-Hagané, 105 (2019), 10 (in Japanese).
• 5)   Q.  Li,  M.  Li,  S.  Kuang and  Z.  Zou: Metall. Mater. Trans. B, 46 (2015), 1494.
• 6)   J. J.  Monaghan: J. Comput. Phys., 82 (1989), 1.
• 7)   J. J.  Monaghan: Annu. Rev. Astron. Astrophys., 30 (1992), 543.
• 8)   J. J.  Monaghan: Mon. Not. R. Astron. Soc., 335 (2002), 843.
• 9)   M.  Nakano and  K.  Ito: ISIJ Int., 56 (2016), 1297.
• 10)   X. Y.  Hu and  N. A.  Adams: J. Comput. Phys., 213 (2006), 844.
• 11)   Y.  Hashimoto,  S.  Yamada,  Y.  Yamaguchi,  K.  Sakai,  T.  Nakajima,  H.  Fujimura and  T.  Kuroda: Jpn. J. Multiph. Flow, 25 (2012), 451 (in Japanese).
• 12)   K.  Szewc: Task Q., 18 (2014), 67.
• 13)   Y.  Kawai and  Y.  Shiraishi: Handbook of Physico-chemical Properties at High Temperatures, ISIJ, Tokyo, (1988).
• 14)   N. A.  Molloy: J. Iron Steel Inst., 56 (1970), 208.
• 15)   T.  Tanaka and  K.  Okane: Tetsu-to-Hagané, 74 (1988), 1593 (in Japanese).