Introduction

The downward injection of process gas through a top lance is the main feature of the Top-Submerged-Lance (TSL) smelting reactor. When the lance is submerged in the liquid slag bath, a violent bubbly flow is generated by the injected gas, creating a large gas–liquid interface, which plays a key role for the kinetics of the smelting process. Clearly, a full understanding of this complex multiphase flow is needed to describe the entire process. Indeed, several studies concerning TSL injection have been published over the years, showing the effort that has been put into gaining a better understanding of these multiphase flows. It should be noted that the high temperatures in the furnaces and the presence of chemically aggressive compounds in the slag impede any direct observation of the slag–gas multiphase interaction in the real process. To gather detailed information on the flow, it is therefore necessary to work on lab-scale models with fluids which are easy to manage, in non-reactive and non-toxic conditions. A complementary approach is the usage of Computational Fluid Dynamics (CFD), which allows multiphase flows to be studied under conditions where detailed experimental measurements are difficult, if not impossible. After thorough verification and validation, which needs to be based on well-defined experimental activity, the numerical model can provide additional data for the analysis of the flow and also be used as a predictive tool. TSL gas injection has been widely investigated for air–water systems and with other viscous transparent liquids, such as water–glycerol solutions.[1,2,3,4,5,6,7] Abranches de Antunes studied top gas injection into a cylindrical vessel, focusing on the bubble formation frequency and its possible dependence on the scale of the reactor.[8] The investigations were carried out with water and both air and helium, to quantify the momentum transfer by injecting gases of different density. He found out that the diameter of the lance barely impacts the bubble frequency, whereas lowering the submerged depth leads to a linear increase in the frequency. Higher bubble detachment frequencies were found in the case of helium, showing the importance of fluid properties in the hydrodynamic description of this system. The frequency of bubble detachment was also analyzed by Neven,[9] who in addition to the water–air system, also worked with glycerin and helium. No significant influence on the frequency was observed while using liquids of different viscosity. Furthermore, the process of bubble detachment was studied in detail and the interaction of the bubbles with the vessel wall was observed, showing how the reactor size affected the dynamics. Iguchi et al. carried out an experimental work where helium was blown into Wood’s metal, which is liquid at around \(100\,^{\circ }{\hbox {C}}.\) They used an electromagnetic sensor (Vives probe) to measure the velocity and turbulence fields in the liquid bath.[10] An electroresistivity probe was used to measure the void fraction, and thus to estimate the bubble size and frequency. The authors suggest that correlations derived from the water–air model can be applied to a Wood’s metal-helium model to predict axial velocity, turbulence distributions, and bubble characteristics.

Some numerical works on liquid metal multiphase flows are also present in the literature. Schwarz and Fröhlich[11,12] investigated the effect of an electromagnetic field on rising bubbles in liquid metal, using the Immersed Boundary Method (IBM) to describe the gas–metal interface. The results of the model are in good agreement with experimental data obtained by means of ultrasound Doppler velocimetry.[13] The authors present a comprehensive study of local bubble movement phenomena, a valid study which complements the experiments.

To date, there has been almost no fundamental research on the hydrodynamic modeling of top-submerged gas injection into slag or metal baths. Wang et al.[14] recently studied top gas injection into viscous liquids using CFD, based on a combined LES-VOF approach. After a comparison with experimental PIV data obtained in a lab-scaled bath of paraffin oil, the model was directly applied to an industrial scale vessel with liquid slag. The hydrodynamics of the gas injection were analyzed in terms of vortex structures, splashing, and surface sloshing. A semi-empirical equation was derived to correlate the gas penetration depth with the modified Froude number. The authors of the present study performed previous numerical investigations on TSL gas injection for non-ferrous applications with the aim of evaluating the effect of viscous and interfacial forces on the multiphase flow.[15] The results showed that important phenomena of the smelting process, such as the gas–liquid interface area, liquid mixing and splashing, strongly depend on the properties of the liquid. This means that the usage of lab vessels which are scaled down based on Froude number similarity may not be appropriate for studying TSL flows in non-ferrous systems, because such scaling does not consider the effect of viscosity and surface tension. These considerations also encourage research on various fluids whose fluid dynamic properties are closer to those of smelting slags.

In another recent publication, the authors also report the application of X-ray radiography to visualize the flow of top-submerged gas injection into liquid metal.[16] This measurement technique has been successfully applied to investigate liquid metal multiphase flows.[17,18,19,20,21,22] The experimental activity of Akashi et al. reveals the complexity of the multiphase flow of TSL Ar injection in the eutectic metal alloy GaInSn. The post-processing applied to the X-ray images allows the flow regime to be characterized and provides the properties of individual bubbles such as size, shape, and trajectory. The bubble detachment frequency is measured as a function of the gas flow rate and submersion depth. A method is proposed for determining the local void fraction based on the Beer–Lambert law, allowing the bubble volume to be estimated. The results are compared with published data for water baths.[23]

The present work is an extension of these previous activities, focusing on the detailed analysis of the flow structure occurring at a high gas flow rate injected at three different vertical lance positions. A CFD model based on the VOF approach is employed and compared with the experimental data, proving the reliability of the numerical approach when dealing with highly turbulent liquid metal bubbly flows and providing a valuable insight into multiphase interaction.

The Experimental Setup

The experiments were carried out in the X-ray laboratory at HZDR. The experimental configuration is almost identical to the measuring arrangement developed by Keplinger et al.[21,22] to investigate the behavior of bubble plumes generated by bottom gas injection. The ternary alloy GaInSn at the eutectic composition was used as a model fluid which is liquid at room temperature. Plevachuk et al.[24] reported a compilation of the main thermophysical properties over a wide temperature range. At ambient temperature and at its eutectic point, GaInSn has a viscosity \(\mu = 2.16\times 10^{-3}\,{\hbox {kg}}/{\text {ms}},\) a density \(\rho = 6360\,{\hbox {kg}}/{\text {m}}^{3},\) and a surface tension \(\sigma = 0.533\,{\hbox {N}}/{\text {m}}.\) Figure 1 shows a schematic drawing of the experimental setup, the vessel and the positions of the gas injection nozzle. The liquid metal was filled to a height of \(H = 144\,{\hbox {mm}}\) into a narrow vessel made of acrylic glass, with a horizontal rectangular cross section of \(144\times 12\,{\hbox {mm}}^{2}\) (see Figure 1(b)). The high attenuation of X-ray radiation in the liquid metal and the dynamics of the bubble motion restrict the dimension of the fluid container along the X-ray beam. Images with a good intensity contrast and low signal-to-noise ratio can only be achieved by using a flat geometry of this kind. Previous measurements[18,21,22] have shown that a gap of \(D = 12\,{\hbox {mm}}\) is a good compromise. On one hand, the image contrast obtained in this configuration is sufficient to clearly visualize bubbles and gas–liquid interfaces; on the other hand, the moderate geometric restriction still allows for the development of sufficient dynamics of both the bubble plume and the free surface of the melt. The acrylic glass material of the walls does not cause any significant attenuation of the X-ray beam intensity. The inert gas Ar (density \(\rho _{\text {g}} = 1.78\,{\hbox {kg}}/{\text {m}}^{3}\) and viscosity \(\mu _{\text {g}} = 1.25\times 10^{-5}\,{\hbox {kg}}/{\text {ms}}\)) is injected through a straight lance made of stainless steel, with an outer diameter of \(d_{\text {out}} = 5.5\,{\hbox {mm}}\) and an inner diameter of \(d_{\text {in}} =5.0\,{\hbox {mm}}.\) A gas flow control system (MKS Instruments) is used to control the gas flow rate of \(6850\,{\hbox {cm}}^{3}/{\hbox {min}}\) (\(Q_{\text {gas}} = 0.11\,{\hbox {l}}/{\text {s}}\)). The system is normed to nitrogen at standard conditions (1 bar, \(0\,^{\circ }{\hbox {C}}\)). Thus, the desired gas flow is adjusted by applying a correction factor. The lance was aligned vertically and immersed in the melt from the top in the center of the container’s horizontal cross section. Three injection positions were chosen with different submergence depths L: \(L = 1/4 H\) (top), 1/2H (middle), and 3/4H (bottom). The measurement time was 40 seconds, starting with the initiation of the bubble injection. The applicability of X-ray radiography for quantitative measurement of bubbles rising in a stagnant liquid was already validated in previous studies.[21]

Fig. 1
figure 1

Schematic drawings of (a) the experimental setup and (b) the vessel and positions of the gas injection nozzle

In all the experiments presented here, the size of the observation window was approximately \(106\times 156\,{\hbox {mm}}^{2}\) (\(1730\times 2560\,{\hbox {px}}^{2}\)), yielding a spatial resolution of each pixel of \(l_{\text {px}} = 0.061\,{\hbox {mm}}.\) The difference between the very high attenuation of X-ray radiation of the liquid metal and the extremely weak attenuation of Ar leads to a clear contrast in the resulting images. The local image intensity correlates with the fraction of gas at a specific position in the \(y{-}z\) plane. The quantitative analysis of bubble parameters is performed by means of digital image processing. The fact that large parts of the volume in the fluid vessel, which are above the free surface, are almost completely filled with gas (with the exception of a few small liquid metal droplets) would lead to a strong overexposure of the recordings. To prevent this, this zone was covered with a thin lead foil and the corresponding image area is not taken into account in the analysis. A detailed description of the image processing algorithm can be found in Akashi et al.[16].

The Numerical Setup

The VOF Model

To simulate the multiphase flow, CFD is used with the Volume Of Fluid (VOF) method to track the interface between Ar and the liquid metal alloy. The VOF approach belongs to the so-called “one-fluid” models, in which only one set of momentum equations is solved:

$$\nabla \cdot {\vec {u}} = 0,$$
(1)
$$\begin{aligned}&\frac{\partial \rho (\alpha ) {\vec {u}}}{\partial t} + \nabla \cdot \rho (\alpha ) {\vec {u}} {\vec {u}} = - {\vec {\nabla p}} + \nabla \cdot \mu (\alpha ) [\nabla {\vec {u}} \\&\quad +(\nabla {\vec {u}})^{\text {T}}] + {\vec {f_{\text {g}}}} + {\vec {f_{\sigma }}}, \end{aligned}$$
(2)

where \({\vec {u}}\) is the velocity vector, p is the pressure, \({\vec {f_{\text {g}}}}\) is the gravitational force, and \( {\vec {f_{\sigma }}}\) is the interfacial force due to the surface tension. In the present work, a RANS approach is used to model turbulence. The \({\textit{k}}\epsilon \)-Realizable model is chosen, in which two additional transport equations are solved: one for the turbulent kinetic energy k and one for its rate of dissipation \(\epsilon \)

$$\begin{aligned}&\frac{\partial \rho (\alpha ) k}{\partial {t}} + \nabla \cdot \rho (\alpha ) k {\vec {u}} = \nabla \cdot \left[ \left( \mu (\alpha )+\frac{\mu _{t}(\alpha )}{\sigma _{\text {k}}}\right) {\vec {\nabla k}}\right] \\&\quad + G_{\text {k}} + G_{\text {b}} - \rho (\alpha )\epsilon , \end{aligned}$$
(3)
$$\begin{aligned}&\frac{\partial \rho (\alpha ) \epsilon }{\partial t} + \nabla \cdot \rho (\alpha ) \epsilon {\vec {u}} = \nabla \cdot \left[ \left( \mu (\alpha )+\frac{\mu _{t}(\alpha )}{\sigma _\epsilon }\right) {\vec {\nabla \epsilon }}\right] \\&\quad + \rho (\alpha )C_{1}S\epsilon - \rho (\alpha )C_{2}\frac{\epsilon ^{2}}{k + \sqrt{\nu \epsilon }} \\&\quad + C_{1 \epsilon }\frac{\epsilon }{k}C_{3\epsilon }G_{\text {b}}, \end{aligned}$$
(4)

where

$$ C_{1}=\max \left[ 0.43, \frac{\eta }{\eta +5}\right] ,\quad \eta =S\frac{k}{\epsilon },\quad S=\sqrt{2S_{ij}S_{ij}} . $$
(5)

\(G_{\text {k}}\) is the generation of turbulent kinetic energy due to the mean velocity gradients, \(G_{\text {b}}\) is the generation of turbulent kinetic energy due to buoyancy, \(\sigma _{\text {k}}\) and \(\sigma _{\epsilon }\) are the turbulent Prandtl numbers, and \(C_{1},\)\(C_{2},\) and \(C_{3}\) are constants. The local fluid properties of density and viscosity are defined as a function of the local fraction of gas and liquid:

$$ \rho (\alpha )= {} \rho _{\text {l}} \alpha + \rho _{\text {g}} (1-\alpha ),$$
(6)
$$ \mu (\alpha )= {} \mu _{\text {l}} \alpha + \mu _{\text {g}} (1-\alpha ), $$
(7)

where \(\alpha \) is the volume fraction, \(\rho _{\text {l}}\) is the density of the liquid phase, \(\rho _{\text {g}}\) is the density of the gas phase, \(\mu _{\text {l}}\) is the dynamic viscosity of the liquid phase, and \(\mu _{\text {g}}\) is the dynamic viscosity of the gas phase.

The volume fraction \(\alpha \) is a discrete function which varies from 0 to 1 and is defined as follows,

$$ \alpha =\left\{ \begin{array}{ll} \,0 &{} \quad \text {if the phase gas is in the cell},\\ \,0< \alpha < 1 &{} \quad \text {at the interface cell},\\ \,1 &{} \quad \text {if the phase liquid is in the cell} \end{array}\right. $$
(8)

and it is advected, according to the advection equation in the conservative form

$$ \frac{\partial \alpha }{\partial t} + \nabla \cdot {\vec {u}} \alpha = 0.$$
(9)

The commercial package ANSYS \({\hbox {Fluent}}{\textregistered }\) (v.17.2) was used to perform the calculations and the numerical setup described in Table I was applied.

Table I Numerical Setup

Figure 2 shows the computational domain of the top configuration and the boundary conditions applied, which faithfully reproduce the experimental setup. The inlet is directly located at the lance tip, avoiding cells in the inner side of the nozzle in order to reduce the grid size. All the three top, middle, and bottom configurations were calculated with the CFD model.

Fig. 2
figure 2

Overview of the computational domain. Common boundary conditions are applied in the calculations, the inlet is placed directly at the tip of the lance, and the outlet is on the upper wall of the vessel. The three dashed lines represent the different lance positions analyzed in this study, 1/4 H,  2/4 H, and 3/4 H (Color figure online)

Post-processing the Data

To make the upcoming sections easier to understand, some steps used to post-process the data obtained from the simulations are explained here.

The averaged chordal void fraction extracted from the X-ray imaging[16] is compared with the time-averaged distribution of the Ar volume fraction, calculated with the CFD model. Indeed, the average volume fraction of Ar can be interpreted as the probability of finding Ar in a specific location during the time of observation. The distribution is then further averaged in the direction of the depth (12 mm) to obtain one single image representing the 3D vessel. A slice averaging procedure over 10 steps is applied, as shown in Figure 3. As described in Reference 16, the distributions derived from the X-ray images are computed after a process of bubble detection and filtering, while the CFD average maps include the overall void fraction. The CFD solutions will therefore show higher values of void fraction in some regions, due to the presence of small bubbles not captured by X-ray radiography. Nevertheless, a comparison of the bubble plume can still be properly accomplished.

Fig. 3
figure 3

Example of slice averaging procedure for the void fraction. For each saved time step (100 Hz), the field is sectioned with 10 slices equally distributed over the vessel parallel to the X-ray beam (a), and an average is computed (b). The resulting distribution contains the information about the variation on that dimension, which makes it comparable with the images provided by X-ray visualization (Color figure online)

The procedure to calculate the bubble frequency from the X-ray imaging is described in detail in Reference 16. Briefly, the image intensity is monitored via a thin observation window (\(100\,{\hbox {mm}}\times 1.5\,{\hbox {mm}}\)), placed at 10 mm above the lance tip. The values are spatially integrated for each time step to obtain a time-dependent signal S(t) as

$$ S(t) = \int _{z=10}^{z=11.5}\int _{y=0}^{y=100} I(y,z)\,{\text {d}}y\,{\text {d}}z.$$
(10)

In the CFD calculation, the void fraction is tracked over a monitoring surface, placed at the exact location of the pixel window used for the X-ray images. Both time signals are then post-processed with a Fast Fourier Transform and Gaussian fit, to determine the main frequency, as shown in Figure 4. In order to remove the initial transient development of the flow from the calculated data, the temporal range of observation is 2 to 30 seconds. The signal was recorded with a frequency of 125 Hz in the X-ray images and 100 Hz in the CFD simulations.

Fig. 4
figure 4

Example of bubble frequency analysis. The power-frequency spectrum obtained with FFT is filtered with a Gaussian function (Color figure online)

Grid Convergence Study

A grid refinement study was accomplished to ensure that the results were grid-independent. To this end, simulations were conducted on three different block-structured grids for the top configuration. The spatial discretization was varied, with major focuses on the injection area, in the bubble plume region and in proximity of the free surface. A short summary is reported in Table II. The simulations were run for 10 seconds (real time) and the grid convergence was established by comparing the pressure drop of the Ar gas between the inlet and outlet and the bubble frequency calculated as described earlier. Figure 5 shows the results of the study: the bubble frequency spectrum was evaluated for the three grids, showing a clear convergence towards the experimentally determined value reported in the plot. A further refinement level would have been prohibitive with the available computing resources. Grid 3 is therefore used in this work, with a minimum cell length of 0.2 mm in the zones of interest and a size of round 1.9 M cells. A time step of \(1\times 10^{-5}\) seconds is employed to keep the Courant number below 1. In Figure 6, a close-up view of the computational grid shows the mesh refinement at the gas inlet boundary. The simulations were run on 240 CPUs, allocated at the HPC Cluster of the Center for Information Services and High Performance Computing (ZIH) at TU Dresden, with an order of computing days of O(10) to complete 30 seconds of real time for the following simulations.

Table II Grid Convergence Study: Numerical Details of the Grids and Corresponding Pressure Drops Are Listed Here
Fig. 5
figure 5

Grid convergence study: bubble frequency analysis (Color figure online)

Fig. 6
figure 6

Details of the computational grid 3, used for simulations in this work. A view of the injection region shows the level of refinement, which permits a resolution of approx. 50–60 pt. per bubble diameter (Color figure online)

Main Properties of the Two-Phase Flow

Qualitative Comparison of Bubble Shape and Ascent Behavior

Figure 7 offers a snapshot of an early stage of the injected bubble for the middle configuration obtained from both the experiment (Figure 7(a)) and the simulation (Figure 7(b)). Corresponding movie sequences are available as Electronic Supplementary Materials. The snapshots show a clear similarity for the bubble shape comparing the numerical results with the experimental data. After injection, the Ar bubble adopts a circular shape on the top section, similarly to cut segments of spheroids, and extends up to 50 mm in the axial direction. The terms spherical cap and ellipsoidal cap are used for such bubble shapes.[25] When they are indented at the bottom, as in the present results, they are often defined as dimpled. These terms have been used for bubbles in three-dimensional geometries. In the present study, however, the vessel has a quasi-2D geometry and the bubble shape could be significantly influenced by the walls. A slight difference of shapes is observed in vicinity of the lance. In the snapshots obtained experimentally, the top part of the bubble hardly attaches to the lance. The contact angle between the gas and the lance is approximately 180\(^{\circ },\) meaning a condition of no wetting. In the numerical snapshot, however, the top part of the bubble clearly attaches to the lance. Indeed, the wettability of the lance and side walls was not explicitly modeled in the CFD calculations, resulting in this slight difference in the bubble shapes. Despite this discrepancy, the integral results shown in the next sections do not deviate significantly.

Fig. 7
figure 7

Snapshot of initial bubbles arising in the liquid metal for the middle configuration from (a) X-ray and (b) CFD results (Color figure online)

Figure 8 contains sequences of bubble snapshots obtained for the middle configuration, from both the experiment (Figures 8(a)–(c)) and simulation (Figure 8(d)–(f)). Both image sequences exhibit similar behavior for the rising bubbles. After detachment from the outlet of the lance, the bubble is strongly deformed. While the center of the bubble is accelerated towards the free surface by the buoyancy force, the surrounding fluid decelerates the bubble at its lateral interface. As a result, an annular rim forms at the lower end of the bubble, the radial extension of which becomes thinner and thinner over time. In some cases, the wake of the bubble already sucks in the next bubble. This leads to an intensive bubble–bubble interaction with the annular rim bulging inwards as shown in Figure 8. The process of coalescence, when two bubbles approach each other, triggers break-up events. If the vertical elongation increases and the lateral expansion of the bubble becomes too small in this area, the shear forces cause smaller bubbles to detach. The newly formed smaller bubbles usually remain in the fluid much longer than the large bubbles, which reach the free surface very quickly. Such intense behavior is hardly seen in lower gas flow rate conditions, as observed by Akashi et al.[16] A detailed analysis of the break-up phenomenon is presented in the next section of the article.

Fig. 8
figure 8

Sequence of images for the middle configuration from (a)–(c) X-ray and (d)–(f) CFD results (Color figure online)

Average Distribution of Void Fraction

A comparison of time-averaged distributions of void fraction obtained with X-ray imaging and CFD modeling is shown in Figures 9, 10, and 11, for the three vertical lance positions. The VOF model is able to predict the average distributions well and thus the position and extent of the bubbles plume formed after Ar injection. In the top configuration, the distribution is clearly symmetrical in both results. Because of the short rising path to the free surface, the bubbles ascend vertically around the lance. A lateral displacement of the plume is observed when the lance is immersed at a deeper position, as for middle and bottom configurations. The numerical model provides a consistent prediction of this asymmetry as well. Indeed, this type of behavior is also well known for multiphase flows in gas-agitated baths for other types of gas injection.[26,27]

Fig. 9
figure 9

Distribution of void fraction for the top configuration. (a) X-ray radiography and (b) CFD (Color figure online)

Fig. 10
figure 10

Distribution of void fraction for the middle configuration. (a) X-ray radiography and (b) CFD (Color figure online)

Fig. 11
figure 11

Distribution of void fraction for the bottom configuration. (a) X-ray radiography and (b) CFD (Color figure online)

Schwarz derived a mechanism to describe the formation of sloshing waves in a liquid bath where gas is injected from the bottom, based on a theoretical analysis of the system.[28] He postulated that the wave motion results from a resonance between the natural motion of the plume and the natural sloshing motion of the free surface. Considering a 2D ideal tank, when the plume is displaced, the surface will raise on the same side. Because of gravity, the bath will swing to the opposite side, carrying the plume with it. If the wave–plume interaction reaches a resonance condition, a self-sustained oscillation can be established. This depends on various factors such as the gas flow rate, the height of the bath (immersion depth for a TSL system), the vessel geometry, and the fluid properties. Schwarz also extended the investigation to three-dimensional geometries with axi-symmetry where rotational motion is possible, for example, cylindrical vessels. He considered the rotational wave as the sum of two non-rotating waves at right angles and 90\(^{\circ }\) out of phase. The rotary sloshing of the wave induces the rotation of the plume, which generates a swirl motion in the liquid bath. The authors believe that the mechanism proposed for bottom injection can also be extended to top injection, but a simple transfer to the specific geometry of the fluid container considered here is certainly not possible. In Figures 10 and 11, the results of the current work show that the time-averaged shape of the plume is asymmetrical, which is equivalent to the fact that the bubble propagation has a preferential direction during the time of observation (28 seconds). Hence, the middle and bottom configurations do not provide a sloshing motion; instead, a permanently asymmetric deflection of free surface is observed. Any rotary sloshing is obviously prevented by the geometry of the fluid container. In fact, as the bubbles rise to the side, corresponding recirculation vortexes are formed, which in turn ensure that the plume is held on the side. The vortex structure of the velocity field within the fluid is also stabilized by the strong confinement of the side walls. To better understand this particular phenomenon, further investigations using three-dimensional cylindrical vessels are needed.

Bubble Frequency

Figure 12 presents the results of the bubble frequency analysis. The bubble frequency was obtained as described in Section III–B. The main frequencies are evaluated for the three top, middle, and bottom configurations for both the experimental data and the simulation results. Excellent agreement is found between experimentally obtained and simulation-based frequencies, further strengthening the validation of the VOF model for Ar-metal systems. The results show a linear increase in the bubble frequency with the lance submersion depth. This behavior is in contrast with the observations by Gosset et al. and Abranches de Antunes,[8,23] who noticed a slight decrease in the frequency when lowering the lance into the bath. They studied the injection of helium/air in a water bath using a 3D cylindrical vessel. Two aspects have to be kept in mind when examining this difference: a different fluid and a different vessel geometry are used in the present work. First, the high liquid density and surface tension in the Ar-GaInSn system can indeed have a strong influence on the bubble detachment mechanism. Obiso et al. pointed out the importance of viscosity and surface tension in the hydrodynamics of a TSL injection and showed that the liquid properties strongly influence the multiphase flow.[15] A deviation from Gosset’s work was also already observed by Akashi et al.[16]. When a dimensionless correlation was applied to link the bubble frequency f to the injected gas velocity u (hence the flow rate), different coefficients were obtained. The dimensionless correlation proposed by Gosset et al., see Eq. [11], between the Strouhal, Weber, and Bond numbers indeed does not take into account either viscosity or surface tension, as

$$ St = C_{1} \left( \frac{We}{Bo}\right) ^{C_{2}}, $$
(11)

where \(C_{1}, C_{2}\) are constants, \(St=fd/u,\)\(We=\rho d u^{2}/\sigma, \) and \(Bo=\Delta \rho g d^{2}/\sigma .\) Nevertheless, the usage of a quasi-2D vessel in the work of Akashi et al.[16] could have a dominant effect on the flow, as already observed in the previous section of this work. Future investigations on a 3D cylinder will also provide further important insights into this question.

Fig. 12
figure 12

Bubble frequency analysis: comparison between X-ray and CFD results, obtained for top, middle, and bottom configurations (Color figure online)

Bubble Size Distributions

This section compares experimental data and numerical simulations regarding the estimation of equivalent bubble diameters. Experimental values of the equivalent bubble diameter were obtained from the three-dimensional void fraction data by post-processing X-ray images, as described by Akashi et al.[16] The corresponding numerical data for the bubble volume are based on the evaluation of the depth-averaged images (presented in Figure 3), which also include the void distribution along y direction, parallel to the X-ray beam. Knowing the distribution of the depth-dependent void fraction \(\alpha (y,z),\) the volume of the ith bubble is computed according to Eq. [12],

$$ V_{{\text {B}}_{i}}=D \int _{A_{i}} \alpha (y,z)\,{\text {d}}A, $$
(12)

where D is the vessel depth and A is the bubble cross section. The equivalent diameter \(d_{{\text {V}}_{{\text {B}}_{i}}}\) is therefore calculated as follows,

$$ d_{{\text {V}}_{{\text {B}}_{i}}} = \root 3 \of {\frac{6\,V_{{\text {B}}_{i}}}{\pi }}. $$
(13)

Bubble size distributions are obtained by applying the post-processing step to each frame of the observation time of 28 seconds. The histograms in Figure 13 illustrate the average number of bubbles per frame \(\bar{N}\) as a function of the equivalent diameter \(d_{\text {V}},\) comparing the experiments and simulations for the three configurations studied.

Fig. 13
figure 13

Comparison of bubble size distributions between X-ray imaging (a, c, e) and CFD simulations (b, d, f) for different configurations of the lance submersion depth: top (a, b), middle (c, d), and bottom (e, f) (Color figure online)

Figure 13 contains the results for the bubble size distributions obtained from the X-ray data (a, c, e) and from the CFD (b, d, f), respectively. Comparable tendencies are observed for both outcomes. First, the gross number of bubbles increases with the lance submersion depth. A plausible explanation for this growth in the bubble number may be the increasing turbulence in the bath as the lance is immersed deeper into the metal. This surely leads to higher probabilities of bubble break-up and gas entrainment from the free surface. The same trend was observed by Akashi et al.[16] when increasing the gas flow rate. Second, two peaks appear in all distributions at large and smaller equivalent diameters, indicating the presence of primary and secondary bubbles, where the primary bubbles are considered to be those immediately after detachment, while the secondary bubbles appear later due to bubble–bubble interaction. The manifestation of this bipolar nature of the bubble size distributions also increases with the lance immersion depth, which confirms the above statement. To calculate the average bubble diameter of the primary bubbles, the same procedure as in Reference 16 was applied to all distributions. Threshold values of 9 mm for the top case and 14 mm for the middle and bottom cases were determined in each configuration as the minimum between the peaks associated with the primary and secondary bubbles in the histograms of the numerical results. The result shown in Figure 14 confirms that the numerical model agrees well with the experiments and shows a slight increase in the bubble diameter corresponding with the lance immersion depth.

Fig. 14
figure 14

Average large bubble diameter (Color figure online)

A discrepancy to the experimental data can be seen both in the total number of bubbles per image, which is overestimated by CFD, and in the position of the peak for the secondary bubbles, which is shifted towards larger bubbles. The scale of the bubble sizes captured by the VOF model depends on the cell size of the computational grid. This acts de facto as a space filter restricting the interface resolution to larger scales and filtering out unresolved small structures.[29] As a consequence, the homogenization of the smaller scales might explain the clearer and higher peak of secondary bubbles in the CFD results around 4 mm. On the other hand, X-ray imaging is also limited at the lower size scale. Some of the smaller bubbles might not be caught and some of them might be lost in the several processing steps required to obtain a binary image (Gaussian blurring, binarization via adaptive threshold, morphological transformation).[16] These considerations make it difficult to compare the results in the lower part of the histograms in detail, without questioning the results discussed above.

Two-Phase Fluid Dynamics of the Bath

Velocity Field

For a better understanding of the fluid flow, some numerical results related to fluid velocities are presented in this section. Figure 15 contains streamlines of the flow field for the three different submergence depths. Average distributions of the void fraction for each case are also shown to explain the mechanism of the bubbly driven flow. Figure 16 contains velocity contours in z direction (Figures 16(a)–(c)) and in x direction (Figures 16(d)–(f)). The domain with an average void fraction above 0.5 was blanked to focus the attention on the region occupied by the liquid metal. For the top configuration, the formation of two symmetric recirculating vortexes is visualized (Figure 15(a)). The bubbles detached from the exit of the lance rise straight, accelerating upward the surrounding liquids. At the same time, the depression created at the wake of the bubbles sucks in liquid metal from the bottom of the vessel. The combination of these two phenomena leads to the formation of recirculating vortexes. The contour plots of the velocities in x direction and z direction show that there is a higher velocity region above the lance exit, where the bubbles have lost most of their initial momentum. The magnitudes of the velocities in z direction are higher than the velocities in x direction. In fact, the liquid flow is mainly driven by the buoyancy forces acting on the bubbles. The negative axial velocities at the walls and the high radial velocities at the free surface show the presence of the recirculating cell. About half of the liquid volume at the bottom of the vessel remains almost unaffected in this case. For the middle and bottom configurations, two recirculating vortexes are asymmetrically formed due to the asymmetric bubble path explained in the last Section IV–B (Figures 15(b) and (c)). The distributions of the axial velocities at the walls and the radial velocities at the bottom of the vessel show that the recirculation involves the liquid throughout the vessel (Figures 16(b), (c), (e), (f)).

Fig. 15
figure 15

Streamlines and time-averaged void fraction for top (a), middle (b), and bottom (c) configurations (Color figure online)

Fig. 16
figure 16

(ac) Contour plots of time-averaged z velocity for the top, middle, and bottom configurations. (df) Contour plots of time-averaged x velocity for the top, middle, and bottom configurations (Color figure online)

The agitation of the slag bath is certainly the central desired feature of the TSL furnace, because of the enhanced reaction rates. However, when the lance is in a deep vertical position, the mixing can involve the lower liquid layers as well, where the slag–matte interface is established. In copper smelting, for example, a slight agitation of this interface could be favorable for the further oxidation of iron sulphide,[30] hence increasing the matte grade:

$$\begin{aligned}&{\text {FeS}}\,({\text {matte}}) + 6{\text {FeO}}_{1.5}\,({\text {slag}}) \\&\quad \longleftrightarrow 7{\text {FeO}}\,({\text {slag}}) + {\text {SO}}_{2}\,({\text {offgas}}). \end{aligned}$$
(14)

A strong recirculation of the lower zone of the reactor is, however, not desired, to avoid the mixing of the two phases. As a consequence, optimally positioning the lance can maximize the performance of the process. Nevertheless, one must consider the complexity of the process, namely other variables (gas flow, temperature, solid fraction) and phenomena (splashing, lifetime of lance, and refractory lining).

Bubble Break-up

The event of break-up at the bottom of the bubbles was briefly described in the qualitative comparison of the numerical results with the experimental data. Additional clarifications of the flow characteristics are given here by means of CFD modeling.

Several bubble break-up mechanisms have been proposed in the literature: surface instability, the shearing-off process, viscous shear forces, and turbulent fluctuations or collisions.[31,32] In liquids with low viscosity, the interfacial shear force is usually negligible. In these cases, shearing-off phenomena appear and the gas velocity profile inside the bubble is thought to be the mechanism for the break-up. The gas inside the slug bubble moves globally upward at the bubble terminal velocity, except at the boundary near the air–water interface, where downward velocities develop because of internal recirculation, dragging the liquid as well. Moving with the liquid film, the gas may therefore penetrate the liquid slug at the rim of the bubble, causing break-up. Figure 17 features a sequence of images during break-up. The iso-surface outlines the location of the bubble boundaries, and a slice located at the middle plane gives detailed information about the velocity field, with visualizations of velocity vectors and a contour map of axial velocity. As described before, negative axial velocities are observed at the peripheral boundaries of the bubble, as a consequence of internal recirculation vortexes. The interface is hence stretched at the bottom and the resulting elongated gas rim detaches in smaller structures, also aided by the upcoming bubble with higher kinetic energy.

Fig. 17
figure 17

Sequence of images during break-up at the trailing edge of the bubble. The distribution map of vertical velocity shows that the gas inside the bubble moves downward at the lateral boundaries (Color figure online)

This brings us to the presence of collisional break-up, which can happen in the case of high turbulent kinetic energy or inertial forces of the hitting bubble, as described by Liao et al.[31] Figure 18 shows the axial component of velocity mapped onto the gas–liquid interface for the last two snapshots of the previous sequence. It can be observed how the trailing bubble, sucked in the wake of the previous one, rises at a higher speed, approaching and colliding with the following skirt and thus contributing to the detachment of the smaller structure. The separated bubble then rises independently to the free surface.

Fig. 18
figure 18

Axial velocity mapped on the bubble surface for the third and fourth snapshot from Fig. 17. The higher velocity of the trailing bubble leads to a collisional break-up event at the bottom rim of the leading bubble (Color figure online)

Hence, the majority of the break-up events may be caused by a combination of the two mechanisms described. As seen in Section IV–D, the presence of secondary bubbles increases with the lance immersion depth. Indeed, lowering the lance into the vessel leads to stronger turbulence and to higher rising velocities for the gas phase.

Analysis of the Free Surface Motion

As already discussed in Section IV–B, the injection of gas into a liquid bath might excite oscillatory modes of the free surface, leading to sloshing waves in 2D domains and to rotational waves in 3D domains. A statistical analysis of the free surface motion was performed to provide a better understanding of this aspect. Figure 19 shows a comparison of the average position of the free surface during the observation time of 28 seconds, between top and middle configurations, extracted as an iso-surface at a void fraction of 0.5 and averaged along the depth. As a consequence of the average gas plume positioning, the free surface presents a perfectly symmetric profile for the top configuration, and a clearly asymmetric step between the two sides of the vessel for the middle configuration. From this information, it can already be assumed that asymmetric average profiles of the free surface (as in Figure 19(b)) cannot result from sloshing behavior. The well-balanced bath height of top configuration might instead represent the average position of a sloshing wave.

Fig. 19
figure 19

Location of the average-free surface for two configurations: (a) top and (b) middle. The dashed line was used to track the volume fraction along time, permitting a frequency analysis of the surface oscillation (Color figure online)

To quantify this statement, the metal void fraction value was tracked during the observation time over a vertical line of 2 cm, placed across the interface at the side of the vessel, close to the wall. These data therefore contain the information about any sloshing behavior of the bath, since it presents a value of 1 in the case of a high wave and 0 for a low wave. The time signal obtained was hence post-processed using the Fast Fourier Transform to get the sloshing frequency of the wave. The result is shown in Figure 20 and confirms the presence of a sloshing wave for the top configuration with a frequency of 2 Hz. No oscillation is observed for the middle configuration. Corresponding movie sequences focusing on the free surface are available as Electronic Supplementary Material, offering the reader visual clarification.

Fig. 20
figure 20

Frequency analysis of the free surface motion. The red curve (top) shows clear sloshing behavior with a frequency of 2 Hz, whereas the flat spectrum of the blue curve (middle) indicates a stable (and asymmetric) wave (Color figure online)

Summary and Conclusions

The hydrodynamics of top-submerged gas injection in liquid metal were investigated by means of X-ray radiography and CFD simulations. The study was performed in a quasi-2D vessel, where Ar gas was blown into a bath of GaInSn, a metal alloy that is liquid at ambient conditions. The CFD model is based on the VOF approach, where the interface between fluids is tracked. An in-depth validation of the modeling approach was possible thanks to close cooperation during the experimental and simulation analysis. The numerical investigation allowed an extensive understanding to be gained of the phenomena observed in the flow, given the full access to local flow variables.

The key outcomes are highlighted as follows:

  • The results of this study demonstrate a high degree of reliability of the numerical modeling of turbulent metallurgical flows with a strong multiphase interaction. This is of fundamental importance in order to close the gap between the study of laboratory scale air–water systems and industrial furnaces, in which slags with very different material properties compared to water are present.

  • Validation by comparing the numerical results with the experimental data produces excellent agreement. X-ray radiography provides an attractive and efficient tool for code validation, since it can image a bubbly flow in a dense and opaque liquid such as a metal alloy and yields CFD-grade data. The CFD model was able to predict time-averaged distributions of the void fraction for three different lance submersion depths, and also predict the presence of asymmetries in the bubble plume. Perfect agreement was also observed for the bubble frequency, an important control parameter for TSL smelting furnaces. Similar trends as observed for the bubble size distributions in the experiments were also predicted by CFD. The numerically calculated equivalent diameters of the primary bubbles were consistent with the experimental data.

  • The analysis of the bubble frequency showed a good agreement between the X-ray data and CFD simulation results. The frequency is found to increase linearly with the lance submersion depth, in contrast to studies available in the literature reporting on gas–water systems. It cannot be definitively clarified at this point whether the observed linear increase is due to the liquid metal used in this study, since the previous investigations on gas–water systems were carried out on cylindrical tanks, and it remains to be clarified to what extent the side walls of the fluid container used here also influences the results.

  • The time-averaged void fraction distributions show an asymmetric bubble plume for the middle and bottom configurations, whereas a symmetric plume appears in top configuration. This behavior is further confirmed and explained with a statistical analysis of the free surface motion. An oscillatory sloshing wave with a frequency of 2 Hz is formed in the top configuration, in which the motion of the bubbles, exposed to a minor buoyancy resistance, move upwards entering into a two-dimensional synchronism with the bath. When the lance is lowered, however, the resistance offered by the liquid increases and the bubble plume, not in synchronism with the surface motion, becomes asymmetric. In a 3D cylindrical geometry it would tend to get a rotational component in its motion. For bottom upward injection, a similar phenomenon has also been observed and a theory explaining it has been developed by Schwarz.[28] The presence of the side walls might therefore create a confinement effect, which blocks the rotational motion of the plume and bath surface, causing a stable asymmetric wave.

  • Besides integral results, the model was also able to reproduce local phenomena, for example, the bubble break-up. When the bubble starts rising, the trailing edge rim tends to bend inward. The high surface tension of GaInSn, together with high relative velocities between the gas and liquid phases, leads to the detachment of smaller structures at the trailing edge. These mechanisms, widely discussed in the literature on bubble dynamics, can be visualized using CFD models.

  • With the aim of obtaining a better understanding of the flow in a TSL smelting furnace, further investigations are desirable in the future. The essential limitations of this study concern the variability of the material parameters and the geometry of the fluid container. Unlike conventional physical modeling in water, the metal alloy GaInSn allows for flow measurements in a fluid of high density and surface tension, but the viscosity of this melt is still several orders of magnitude below that of typical molten slags. Therefore, further studies should also consider the modification of the alloy viscosity. The quasi-two-dimensional geometry of the fluid container represents a significant restriction with regard to the development of three-dimensional flow structures. Phenomena such as the rotary sloshing or the 3D dynamics of the bubble plume can obviously not be investigated in this setup. Therefore, further experiments in cylindrical tanks are planned to approach a geometrically scaled setup in which the hydrodynamics of the TSL smelter are reproduced.