Dynamics of Rising Bubbles in a Quiescent Slag Bath with Varying Thermo-Physical Properties

The motion of bubbles in a liquid slag bath with temperature gradients is investigated by means of 3D fluid dynamic computations. The goal of the work is to describe the dynamics of the rising bubbles, taking into account the temperature dependency of the thermo-physical properties of the slag. Attention is paid to the modeling approach used for the slag properties and how this affects the simulation of the bubble motion. In particular, the usage of constant values is compared to the usage of temperature-dependent data, taken from models available in the literature and from in-house experimental measurements. Although the present study focuses on temperature gradients, the consideration of varying thermo-physical properties is greatly relevant for the fluid dynamic modeling of reactive slag baths, since the same effect is given by heterogeneous species and solid fraction distributions. CFD is applied to evaluate the bubble dynamics in terms of the rising path, terminal bubble shape, and velocity, the gas–liquid interface area, and the appearance of break-up phenomena. It is shown that the presence of a thermal gradient strongly acts on the gas–liquid interaction when the temperature-dependent properties are considered. Furthermore, the use of literature models and experimental data produces different results, demonstrating the importance of correctly modeling the slag’s thermo-physical properties.


I. INTRODUCTION
THE design and control of several metallurgical processing operations involve the knowledge and application of multiphase fluid mechanics, almost always linked to heat and mass transfer phenomena. It is precisely to boost these transfer processes that bubbling flows are realized in various reactors, with the aim of increasing the interfacial area between the gas and liquid phases. Examples include desulfuration and decarbonization processes in the ladle steel metallurgy, or pyrometallurgical operations, such as the primary smelting in Top Submerged Lance (TSL) furnaces and secondary matte and slag conversion with bottom-, side-, and top-blown vessels. The motion of gas bubbles in liquid media is a key aspect of these processes, directly affecting design, operation, and productivity. This deep level of knowledge is nowadays required for a full understanding of metallurgical operations and the viability of a circular economy system in the metals market. Thermo-fluid dynamics, together with thermodynamic and kinetic studies, process simulation, and plant optimization, is indeed part of the technical tools needed to optimize metallurgical processes and evaluate their resource efficiency. [1,2] Because of its relevance in energy and process engineering, the fluid dynamics of rising bubbles have been widely studied for decades. [3,4] Early investigations can be found in the work of Haberman and Morton, [5] who experimentally investigated the dynamics of rising bubbles in various liquids, ranging from water to oils or syrup. For each system studied, a drag correlation was proposed as a function of the Reynolds number. The results, categorized in terms of the bubble terminal velocity, path, and shape, offer a broad and detailed overview of bubble dynamics. An important piece of research was carried out by Bhaga and Weber, [6,7] who experimentally determined the dynamics of bubbles rising in viscous liquids, studying the shapes, terminal velocities, and wake development. At present, this pioneering work represents the main valuable reference source for such flows and a well-established benchmark test case for the validation of numerical solutions. The bubble analysis is based on images taken with a camera moving at the same bubble rising velocity. The properties of the liquid were varied by using different concentrations of aqueous sugar solution, covering a range of Morton numbers from 7:4 Á 10 À4 to 850. Thanks to this, one of the first complete bubble shape regime maps was proposed.
The increase in computing power and the development of modeling techniques allowed numerical investigations of flows around bubbles to be carried out, starting under the assumption of a spherical shape and moving on to the modeling of ellipsoidal bubbles. [8][9][10][11] One of the first numerical studies is by Moore, [12,13] who proposed a drag coefficient for spheroidal bubbles as a function of the Reynolds number and bubble aspect ratio, using a mathematical analysis of the boundary layer around the bubble and the wake region. Numerical simulations of complex bubble shapes became possible with the development of front-tracking methods, such as the Volume of Fluid (VOF) and Level-Set (LS) methods, [14,15] which can track the interface between the gas and the liquid. Several works can be found in the literature where Computational Fluid Dynamics (CFD) is applied to analyze bubbles' shape and motion. [16][17][18][19][20] Most of the investigations present a combined experimental and numerical study, proving the reliability of the modeling approaches. Others are purely computational analysis, as in the work of Hysing et al., [21] who proposed a benchmark for the two-dimensional modeling of bubble dynamics, with the purpose of establishing a validation test case as a reference source, considering the absence of analytical solutions and the difficulties involved in measuring such flows experimentally. Ohta and Sussman [22] carried out a CFD investigation based on the Coupled Level-Set Volume of Fluid (CLSVOF) method on the motion of single-skirted bubbles rising in a viscous liquid. They performed a sensitivity analysis on the skirt thickness by varying the density and viscosity ratio of the two fluids, finding out that the skirt thickness increased in line with the density ratio, whereas the effect of viscosity ratio also depended on the density ratio. In a recent comprehensive study, Sharaf et al. [23] carried out an extensive experimental and numerical investigation of single bubbles rising in quiescent liquids. Using a glycerol-water solution, more than 300 different configurations were tested, leading to the compilation of a phase diagram in the Galilei-Eo¨t-vo¨s plane, where all the different bubble regimes are reported. Moreover, Direct Numerical Simulations (DNS) were conducted, which agreed well with the empirical observations. In agreement with Ohta and Sussman, Sharaf et al. found out that the viscosity and density ratio are two important parameters for the characterization of bubble dynamics, especially for transitions from one regime to another. The same group presented further numerical investigations focusing on path instabilities and break-up phenomena. [24,25] Their results also showed the importance of a 3D approach in modeling bubble dynamics, since except for the axi-symmetric one, all the regimes present an asymmetric and inherently 3D nature. Bubble dynamics have also been investigated in molten metal systems, [26][27][28][29][30][31][32][33][34] which require sophisticated apparatus, such as X-ray radiography and ultrasound velocimetry, to measure the flow and provide data to validate numerical investigations, as in the previous work of the author. [35] In the present study, the authors focus on the motion of bubbles in the fayalitic slag, typical of smelting processes in the non-ferrous metallurgy. In addition to the metals' high density and surface tension, liquid slags are heavily viscous fluids, with values in the range of 0.1 to 1 PaÁs. Due to the high temperatures, necessary to liquify the oxides, and the opacity of the liquid, experimental observations are challenging, and such research has not been described in the literature to the best of our knowledge. A few numerical investigations have been performed with CFD analysis. Cold models have been applied to study the hydrodynamics of gas injection in iso-thermal liquid slag, [36][37][38] while simulations of hot and reactive systems were conducted by Huda et al. [39,40] They studied the process of zinc fuming using CFD for a Tuyere-Blown and a TSL furnace, solving the heat transfer and the multiphase chemical interactions. Although the application of an Euler-Euler framework does not allow the interface between the phases to be reconstructed, or offer a detailed insight into the bubbly flow, the physical behaviors of the system were reproduced, and the fuming rate trends obtained were similar to experimental campaigns. The physical properties of the slag were kept constant, at the process operating temperature. Nevertheless, metallurgical furnaces such as the TSL smelter present a temperature distribution which is not necessarily homogeneous. Together with the presence of a certain solid fraction and a heterogeneous species distribution, this causes local variations in the thermo-physical properties in the bath, influencing the dynamics of the bubbles which move in it.
In other technical fields beyond metallurgical applications, the dynamics of bubbles rising in a bath with varying properties have been studied, albeit not extensively. Sahu [41] wrote a review article on the dynamics of bubbles rising in viscosity-stratified liquids. He went through different application cases, considering the variation of viscosity in space, and examining in detail the variation due to temperature gradient and the non-Newtonian nature of the fluid. He concluded by underlining the lack of related works in the literature, despite the fundamental importance of this knowledge for describing processes with such features. Regarding viscosity, Premlata et al. [42] applied a 2D axi-symmetric CFD code to study the dynamics of rising bubbles in viscosity-stratified liquid. They compared the terminal bubble shape obtained from simulations with constant and linearly increasing viscosity. The bubble shape and, therefore, dynamics were found to be very different in the two configurations at high Eo¨tvo¨s and Galilei numbers. This was attributed to the fact that, while the bubble rises, liquid with a different viscosity migrates in the wake region, influencing the bubble dynamics.
Minor effects were observed at lower Eo¨tvo¨s and Galilei numbers. Tripathi et al. [43] carried out a theoretical and numerical analysis of the bubble rise, taking into account how the surface tension depends on the temperature. They studied the case of a bubble rising in a cylindrical vessel with a vertical variation of surface tension, considering its linear and parabolic dependence on the temperature. They found out that the motion of a bubble rising in a self-rewetting liquid (second case) can be reversed. The work was extended by Balla et al., [44] who used a 3D model. The migration of bubbles towards the wall was observed with small Bond numbers and was not visible with a 2D approach. The study of convective flows generated by surface tension gradients, also known as the Marangoni effect, is present in the References 45 and 46. Hardy et al. [47] experimentally investigated the motion of bubbles in a vertical temperature gradient. They injected air bubbles into viscous silicone oil subject to a temperature gradient. For certain gradients, an oscillatory behavior by the bubbles was observed, with their diameter and vertical position varying with an almost sinusoidal trend. They found that it was possible to stop the bubble motion because of a balance between the buoyancy and Marangoni forces.
The goal of this investigation is to evaluate the correct modeling of the thermo-physical properties of a fayalitic slag, by means of a CFD study of rising bubbles. Given the thermodynamic complexity of a slag, its physical properties depend on its composition, temperature, and pressure. Metallurgical operations, such as a TSL smelting unit, operate at a constant pressure, which excludes one variable from the problem. Assuming that there are no chemical interactions with the gas phase, no feed stream, and, therefore, no change in composition, the slag properties depend only on the temperature. This is the case when process gas (fuel, air) is injected into a fully converted TSL slag bath. This is a test case used by the authors to investigate the hydrodynamics of a multiphase flow in a TSL furnace, taking into account a real slag and including the submerged fuel combustion, necessary to keep the bath liquid. [48,49] Therefore, particular attention is paid here to the temperature dependency of the physical properties, comparing the implementation of models available in the literature and data experimentally measured in-house. The analysis of bubble dynamics in such flows generates proposed guidelines for the physical modeling of metallurgical processes.

II. BUBBLE DYNAMICS AND SHAPE REGIMES
The dynamics of a bubble rising in a quiescent liquid can be described with non-dimensional numbers representing the relationships of the hydrodynamic forces acting on it. The Reynolds number (ratio of inertial and viscous forces), the Eo¨tvo¨s number (ratio of gravitational and interfacial forces) , and the Morton number (ratio of viscous and interfacial forces) where q, l , and r are the density, viscosity, and surface tension of the liquid, respectively, and D eq and u t are the equivalent diameter and the terminal velocity of the bubble. Mathematical correlations of the three numbers are available in the literature with which one can estimate the terminal velocity or the shape of the bubbles, knowing the equivalent diameter and the physical properties. The same information has been translated in graph, [4] which is well known as the bubble-shaped regime map ( Figure 1). Depending on Re, Eo, and Mo, the shapes of bubbles have been classified as follows: Spherical, when viscous and interfacial forces dominate over inertia and the ratio of the minor and major axis is higher than 0.9; Ellipsoidal, for oblate bubble shapes where the interface is convex overall, viewed from the inside; Spherical-cap or Ellipsoidal-cap, when shapes similar to segment cuts of spheroids are formed. If a skirt is formed at the trailing edge, the bubble can also be defined as dimpled or skirted; Random wobbling, which appears at lower Mo where the non axi-symmetric shape leads to oscillatory motions of the bubble. Examining the TSL process, furnaces generally operate at temperatures between 1150-1350 C, in order to establish a desired slag viscosity value. Nevertheless, the temperature is not homogeneously distributed in the bath during time. [50] The combustion region around the lance and the oxidation zones are heat sources in the furnace, whereas colder spots can be found at the peripheral parts of the bath, as well as in the feed stream entry point. In addition, the solid fraction in the bath, generated from the feedstock injection and from the formation of magnetite, is not necessarily homogeneous. Though this is not the subject of the present work, in a real reactive system, the composition of the slag is not homogeneous in the bath either. These phenomena produce varying local physical properties in the slag and, therefore, affect the motion of the bubbles.
Considering the material properties of the slag examined here (see next section), the specified temperature range relates to one portion of the bubble regime map, shown in red in Figure 1. The area of interest for the present work concerns, therefore, dimpled and ellipsoidal-cap bubbles.

III. THERMO-PHYSICAL PROPERTIES OF THE SLAG
The slag examined in this study originates from a secondary slag converter from a copper production line, and so is practically free of copper. The slag is used in the pilot-scale TSL plant, operating at TU Bergakademie Freiberg, where the hydrodynamics of gas injection with submerged combustion are studied by combining experimental investigations and CFD modeling. [48,49] For the purpose of determining the physical properties, the chemical composition is simplified to the major components and used as a reference for the calculations, see Table I.
Under the assumption that chemical interactions do not take place between the liquid and gas phases, at a given constant pressure, the composition is fixed and thermo-physical properties only depend on the temperature.
As already discussed, here the authors investigate the dynamics of bubbles in a slag bath with temperature gradients by means of CFD simulations. One crucial aspect of the study is the comparison between the usage of models available in the literature and data measured in-house for the implementation of the physical properties. An overview is presented here.

A. Models
Given the composition, a thermodynamic calculation was performed in FactSage Ó to determine the solidus and liquidus points and the solid fraction-temperature curve in the mushy region. FToxid and FTmisc were used, allowing all phases and possible two-and three-phase immiscibility. The physical property models considered in this analysis are selected from the literature on slag. The temperature-dependent viscosity is calculated applying the Riboud model, [51,52] based on the Arrhenius type function where the parameters A and B are calculated from the mole fractions of the components. To take into account the presence of the solid fraction in the mushy region, a second-order power series approach (extension of the Einstein equation) was adopted here. [53] The correction factor, the so-called effective viscosity, is defined as where ½l is the intrinsic viscosity (½l ¼ 2:5 for spheres), B is the Huggins coefficient (B ¼ 6:2 for spheres), and / is the solid fraction, calculated in Fact-Sage Ó . This will be referred to, here, as the Einstein-Riboud model. The model proposed by Lange is employed for the temperature-dependent density of the slag. [54] The surface tension, thermal conductivity, and specific heat are modeled as proposed by Mills in References 55 and 56.

B. Experimental Measurements
The slag temperature-dependent viscosity, density, and surface tension were also measured experimentally to provide more accurate data. The experiments were conducted in the laboratories of TU Bergakademie Freiberg. Given the inability to perform experimental measures for specific heat and thermal conductivity, and assuming that these have a secondary effect on the bubble dynamics, the models from literature are always applied.

Viscosity
The viscosity behavior of a slag heavily depends on the partial pressure of oxygen as this affects the ratio of ferric to ferrous iron. [57] To recreate a gas atmosphere that resembles this partial pressure, this ratio was measured from the original process sample via Mo¨ßbauer spectroscopy. The spectrometer was built out of modules from the WissEl GmbH with a movable 100 mCi Co 57 -gamma radiation source within a Rh-matrix and sample holder. The result of this measurement was that 96.2 pct (molar) of the total iron inside the original sample was ferrous iron. The gas composition for the viscosity measurement was calculated via a FactSage calculation using the FactPS, FToxid, and FTmisc databases. All possible solution species were selected, including the sample composition and a mixture of CO and CO 2 , where the mass of the gas is 20 times the mass of the slag. The increased amount of gas ensures that the gas composition changes negligibly when the resulting thermodynamic equilibrium is calculated. [58] The resulting gas mixture contained 31 vol pct of CO and 61 vol pct of CO 2 in order to achieve a FeO/Fe 2 O 3 ratio similar to the original process slag within the viscosity measurement.
For the viscosity determination, the measurement chamber was continuously purged with the calculated gas mixture. The viscometer is of Searle-type (Anton Paar MCR 302 coupled with HTM Reetz LORA resistively heated furnace). Both the crucible and the rotating bob that come into contact with the slag are made from a platinum-rhodium (80/20) alloy. The material warrants sturdiness and immiscibility with the sample. [59] The temperature was measured with type B thermocouple at the bottom of the crucible that is cross-checked against a second type B thermocouple above the crucible. The wide-gap system was calibrated using the DGG standard glass I (lime soda silicate glass) as well as silicon oil viscosity standards. For the viscosity measurement, the slag was heated to 1550 C and stirred for 1 hour to achieve homogeneity. At each temperature point, the slag viscosity was measured for 30 min at a shear rate of 24,8 1/s to reach the equilibrium viscosity. Between the measuring points, the cooling rate was 5 K/min. Data are shown in Figure 2.

Density and Surface Tension
The Maximum Bubble Pressure (MBP) technique was employed for the measurement of both the density and surface tension, using molybdenum crucible and capillary tube. A detailed description of the experimental setup can be found in the work of Korobeinikov et al. [60] At each individual temperature 7 different immersion depths of the capillary tube were tested in order to determine the density from the linear interpolation of the data, following the Young-Laplace equation derived based on the meniscus formed at the injection point being spherical in shape. To account for the gravitational effects on the meniscus shape, which was   in fact non-spherical, the Schro¨dinger correction was applied, where P r ¼ P max À qgh. The capillary depth was set in the range from 10 to 22 mm and the temperature varied in the range from 1100°C to 1400 C. Data are shown in Figure 3.

A. Numerical Method and Test Cases
The commercial software ANSYS Fluent Ó (v.19.2) was used to perform the calculations, with the setup shown in Table II. For the sake of simplicity, the full mathematical formulation is not reported here and can be found elsewhere. [61] The CFD model is based on using the VOF as a multiphase model to track the interface between the gas bubble and the liquid. Since the slag bath is not iso-thermal, the energy equation is solved, and the temperature-dependent physical properties, discussed in the previous section, are implemented as User-Defined Functions (UDFs). When using the VOF method to calculate bubble dynamics, all forces acting at the interface result directly from the integration of the thermo-fluid dynamic equations and the specification of sub-models is not needed. In order to capture any asymmetries in the shape and motion of the bubbles, the simulations are fully 3D. The computational domain is shown in Figure 4: an initially spherical bubble of diameter D = 0.025 m rises up in a quiescent slag bath inside a cylindrical vessel, whose diameter is equal to 5D, to prevent the wall from affecting the dynamics of the bubble. [62] A wall boundary condition is applied at the bottom and sides of the cylinder, with the pressure outlet on the top. As described in Sect. 4.3, a grid resolution of 0.5 mm is necessary to ensure that a grid-independent solution is found. Nevertheless, the application of this discretization level over the entire domain would generate a prohibitive grid size with the actual available resources. To circumvent the issue, a structured grid with the desired resolution is used in the central vane of the vessel, where the bubble moves during the ascent. The external domain extending as far as the wall is meshed with polyhedrons ( Figure 5). This produces an acceptable mesh size of 9:6 Á 10 6 cells. A time step of 1 Á 10 À4 s is applied and the Courant-Firedrichs-Lewy (CFL) number is always kept below 0.1. The simulations were run on 320 CPUs, and for each case, an order of computing days of O(1) was accomplished to complete the simulation of 1 s of real time.
Slag metallurgical reactors are considered as reference to determine the thermal configurations. For example, in TSL smelting reactors, the bath shows local temperature gradients even though the average reference temperature is kept constant to ensure a certain slag viscosity. In order to safely operate the furnaces, wall temperatures must be low enough-usually close to slag solidus point-to form a freeze lining that protects the refractories. [63][64][65][66][67] The temperature gradients will then lie between the maximum bath temperature in the combustion zone and the slag freezing temperature at the wall. Investigating a TSL zinc-fuming process, Huda et al. observed three thermal zones in the slag bath: the submerged combustion zone (1680 K to 1780 K), the central bath zone (1500 K to 1600 K) and the bottom bath zone (1380 K to 1430 K). [50]   To evaluate the impact of the temperature variation of the bath on the bubble dynamics, three benchmark test cases are analyzed, applying radial, horizontal, and vertical temperature gradients, where the thermal effects are decomposed, in order to have a detailed insight and a fundamental understanding of the thermo-physical interactions. A bath reference temperature of 1500 K is chosen, with a variation range of ± 100 K. Radial and horizontal gradients develop over the reference length D, the vertical gradient over the axial rising path. Figure 6 shows the initial condition of the described setups. The temperature of the gas phase in the bubble is set as constant and equal to 1600 K. In a gas injection system with submerged combustion, the gas phase is hotter than the liquid bath, as a consequence of the heat of reaction. The gas inside the bubble is air and the chemical interaction with the liquid phase is not considered.
According to the dependencies described in the previous section, the physical properties of the slag vary in the direction of the applied temperature gradient. Hence, for each configuration test case, the simulation is run with three setups. First, the properties of the slag are set constant and equal to the values from the literature models at the reference temperature of 1500 K. Second, the temperature-dependent properties from the literature models are applied. Third, the properties of density, viscosity, and surface tension measured in-house are adopted.

B. Validation Assessment
The use of the VOF method for the modeling of bubbly flows in complex liquids has been largely validated by the author in Reference 35. To also prove its applicability to fluid systems more similar to the fayalitic slag under investigation, three different ellipsoidal-cap bubbles were selected from the References 6 and 7. This allows the CFD model to be validated in the red area of the bubble regime map, shown in Figure 1, which covers the Re-Eo-Mo range of the slag physical properties object of the present analysis. The studied configurations are summarized in Table III, which reports the physical properties of the liquid phases used in the experiments.   7-Qualitative comparison of bubble shapes between CFD results (b) and experimental observations (a). [7] The configuration of cases B,C, and D is reported in Table III. (a) is reprinted with permission from Ref. [7]. regard to both the leading top curvature and the cavity formed at the rear. A quantitative comparison was also conducted in the evaluation of the Reynolds number at the bubble terminal velocity and the bubble aspect ratio, reported in Figure 8, which shows a perfect agreement with the reference.

C. Grid Independence Study
The test case B from the previous section was used to perform a grid refinement study and establish grid-independent solutions. Three levels of mesh refinements were calculated varying the spatial step of the block-structured volume from 0.5 to 1 and 1.5 mm, as summarized in Table IV. The bubble shape and position are compared in Figure 9 at the time of 1.5 s. It is evident that the coarse resolution is not enough to correctly resolve the bubble shape. As a matter of fact, an overestimation of the aspect ratio of 3.09 pct means a higher drag force acting on the bubble and therefore a lower velocity. As shown in Figure 10, both middle and fine grids predict the aspect ratio of the bubble with an error rate below 1 pct. The errors in the plot are relative to the reference test case. The error rate for the Reynolds number stops at around 6 pct with the middle level of refinement and remains constant with further resolution. It follows that grid convergence would be accomplished with the middle-sized grid. Nevertheless, the finer grid is used in the simulations in the present work, since break-up phenomena may occur under certain circumstances and the resolution of detached smaller structures must be ensured.

V. RESULTS AND DISCUSSION
The analysis of the results is based on the observation of several characteristics of bubble dynamics, such as the rising path, bubble shape, terminal velocity, and break-up, and the evolution of the gas-liquid interface area. Obviously, all the aspects are interconnected and help the reader to get an overview on the bubble-slag interaction.
A. Single Rising Bubble

Rising path
The paths of the rising bubbles are represented in Figures 11(a) through (c), where the three modeling approaches for the slag properties are compared on baths with a radial, horizontal, and vertical temperature gradient. The center of mass (CoM) of the gas phase was monitored over time and plotted over the three spatial coordinates. The first phenomenon which becomes apparent is the straight rise of the bubble when constant values are used for the physical properties of the slag (green curves). Indeed, the presence of a temperature gradient has no effect and the motion of the bubble that is perfectly axi-symmetric. The vertical gradient (Figure 11(c)) has no effect on the direction of the path, while deviations from the center line in the order of 1 mm are observed for the radial temperature gradient (Figure 11(a)) when the experimentally measured properties are used. On the other hand, the presence of a transverse gradient generates a larger deviation from the axis for both cases using literature models and experimental data, with a maximum of 18 mm. At first glance, this behavior might seem to be related to the Marangoni effect, which causes a bubble to migrate towards areas with a lower surface tension. A surface tension gradient along the bubble surface generates thermocapillary stresses on it, modifying the dynamics of the bubble itself. [45] However, this is not the case, because the variation in the surface tension with the temperature is minimal and presents an opposite derivative between the Mills models and the experimental data ( Figure 3). This should lead to an opposite path deviation from the axis, which is not observed. The departure from the center line has to be attributed to the strong viscosity variation, which extends almost an order of magnitude over the considered temperature range. As observed for the case with a radial gradient, the usage of experimental data produces additional oscillations in the order of 1 mm from a certain point of the rise.   Figures 12, 13, and 14 show the bubble shapes for all the three gradient configurations, as well as for the three modeling approaches of the slag properties. All snapshots are extracted at t = 0.75 s and do not compare the axial position of the bubbles, but only their shape. Corresponding movie sequences are available as electronic supplementary material and help the reader to picture the subject matter.

Bubble shape
The explanation of the unstable path oscillations, which are observed when using the experimentally measured data, relies on the presence of bubble break-up phenomena, which are strongly evident for the configurations of the radial and horizontal thermal gradients. A slight detachment of smaller structures also appears with the vertical thermal gradient, although practically without affecting the bubble dynamics. The smaller bubbles detached from the trailing edge remain in the wake area of the main bubble, causing the fluctuation in the center of mass of the gas.
In addition to the break-up events, the main body of the bubbles exhibits some differences in shape depending on the slag modeling approaches. To quantify these discrepancies, the aspect ratio was calculated as follows: considering only the main bubble and not the smaller structures in its wake. High values of A R correspond to flattened and oblate shapes, whereas smaller values represent prolate bubbles. Table V summarizes the outcomes. A general decrease in the aspect ratio is observed when using the temperature-dependent slag properties and, as expected, stronger deviations are related to the radial and horizontal gradient configurations, with an A R reduction of 47 and 32 pct with the usage of constant properties and experimental data. A strong contribution is certainly due to the break-up events, which modify the volume of the main bubble and hence its exposure to the temperature range. This drastic change depends entirely on the physical properties of the liquid phase and reveals the importance of their correct modeling in a CFD calculation based on front-tracking methods.

Bubble terminal velocity
The bubble rising velocities for all cases are plotted in Figure 15. As a consequence of the shape change, the bubbles rise at different velocities depending on the gradient and configuration of physical properties. In the radial thermal gradient bath (Figure15(a)), an increase of around 20 and 45 pct is observed when changing from constant properties modeling to the literature and experimental data models, respectively. An increase of around 23 pct is also present in the horizontal gradient configuration (Figure15(b)), until a sudden deceleration in the last 0.1 s due to the increased detachment of smaller bubbles in the wake. The change of temperature along the axis almost never affects the rising velocity when temperature-dependent slag physical properties are applied (Figure15(c)).  The relationship between the aspect ratio and rising velocity is evident and has been widely discussed in the References 68 through 70. Several analytical and empirical correlations have been proposed, where the terminal velocity u t is a function of the aspect ratio of the bubble, as well as of the liquid properties and the bubble type (spherical, ellipsoidal). Tomiyama et al. [69] derived a theoretical u t correlation for distorted spheroidal bubbles rising in a quiescent liquid bath. The model is based on the momentum jump condition at the gas-liquid interface, under the assumption of a surface tension-dominated flow. The terminal bubble velocity is, therefore, calculated as: where E ¼ 1=A R , d eq is the bubble equivalent diameter, and c is a geometric parameter equal to 2 under the assumption of spherical-cap shape. When the correlation is applied to the cases in this study, the liquid properties r, q l , and Dq are considered as averaged in the temperature range, and the bubbles are simplified to spherical-caps. The correlation graph reported in Figure 16 compares the terminal velocities at t ¼ 0:75 s obtained from CFD and from the Tomiyama correlation. Despite the assumptions of hemispherical bubble shape, surface tension-dominated flow and averaged physical properties, the results correlate well, also showing a certain consistency of Tomiyama's theory in the case of slightly different systems and demonstrating the reliability of the current computational approach.   The correct calculation of the bubble rising velocity is fundamental in terms of reactor modeling, since it determines the bubbling frequency, a key control parameter of metallurgical furnaces based on bubbling flows.

Gas-liquid interface
The gas-liquid interface represents another crucial parameter of reacting bubbling flows because it controls heat and mass transfer from one phase to the other. The usage of a front-tracking multiphase model, such as the VOF, permits the geometrical reconstruction of the interface and, therefore, enables it to be tracked over time. Figure 17 shows how the interface evolves over time in all the cases studied. In the first part of the rise, the interface is higher when constant properties are used, which agrees with the A R data discussed above. At a constant volume, the surface of an oblate ellipsoidal-cap shape is indeed greater than a prolate one. Nevertheless, Figures 17(a) and (b) reveal a sudden increase in the size of the interface at around 0.35 seconds with the experimentally measured physical properties. This has to be attributed to the appearance of break-up events, which produce smaller bubbles and, therefore, increase the surface-to-volume ratio. For both cases, the usage of experimental data leads to a surface that is estimated to be 40 to 50 pct greater compared with the usage of literature-based physical properties. It is relevant to note that, although temperature-dependent variables were used in both cases, the results can significantly vary regarding the interfacial surface. In the case of a vertical temperature gradient, shown in Figure 17(c), the usage of literature and experimental data gives comparable results, slightly below these calculated with constant physical properties.

Bubble break-up phenomena
As already mentioned, break-up events at the bottom skirt of the main bubble appear in some of the configurations, mostly when applying the temperature-dependent physical properties from experiments. The formation of satellite structures has a direct influence on the bubble dynamics, increasing the gas-liquid interface and bringing instabilities to the motion. According to the literature, four main mechanisms of break-up can occur depending on the flow features: turbulent fluctuation and collision, viscous shear stress, shearing-off process, and interfacial instability. [71] Given the high values of slag viscosity and surface tension, it is evident that the observed break-up falls within the viscous shear stress and shearing-off typologies, which are equivalent for bubble diameters larger than 5 mm. [72] The mechanism provides that small bubble structures shear from the bottom skirt of the cap-shaped bubble for viscosity ratios much smaller than 1 (l g =l l ). The detachment is due to an unbalance between the viscous shear and the surface tension at the interface. This process can be precisely observed in the studied cases. Figures 18, 19, and 20 help illustrate this explanation. The local Capillary number defined as is computed and mapped on the gas-liquid interface.
Here l l is the liquid viscosity, d the bubble diameter, _ c the shear rate, and r the surface tension. The deformation of a bubble because of viscous shear stress mainly depends on the Capillary number, which represents the ratio of viscous and surface tension forces. It has been shown that deformations increase with Ca. [71] For each thermal gradient configuration, the local value of Ca is below 10 for the cases of constant and literature-based properties, as shown from Figures 18, 19, and 20(a), (b), and break-up appearance is negligible. However, when the value of Ca increases, detachment starts at the bottom rim. Moreover, as one can see from Figures 18, 19, and 20(c), the departure of secondary bubbles takes place in the areas with highest Ca, where the shear rate is maximum. The reason why this occurs when using the experimentally measured data is, therefore, to be attributed to the respectively higher and lower values of viscosity and surface tension over all the temperature range, in comparison to the other two approaches. Eventually, this results in a flow condition at which the shear stresses and, therefore, Ca are higher, leading to break-up. After shearing-off from the bottom skirt, the satellite bubbles find themselves in the wake of the leading bubble. This is an area of relative depression, where the adjacent liquid is dragged up at the same velocity of the gas, as can be observed from Figure 21. Therefore, they will not experience any resistance and will follow the main bubble upwards due to inertia.

B. Two Rising Bubbles
Additionally, the simulation of two consecutive rising bubbles was performed under the configurations already discussed. In fact, in all metallurgical applications where bubbling flows are involved, the bubbles evolve in a sort of bubble chain because of the continuous gas injection. In some cases, for example in the TSL furnace, the interactions are so strong that two or more bubbles may collide, break or coalesce. [35] In the presence of thermal gradients, depending on which kind of bubble   interaction takes place, the gas is exposed to different physical properties distributions, which might lead to stronger deviations than those observed with one single rising bubble. Figures 22, 23, and 24 show the status of the bubbles at t = 0.75 s, for all the configurations analyzed. Corresponding movie sequences are available as electronic supplementary material and help the reader to visualize the ascent and the behavior of the two bubbles. In each case, after release, the second bubble is sucked into the wake of the leading one and hence accelerated. This behavior was observed in previous works focused on TSL gas injection. [34,35] Because of the acceleration, the two bubbles collide and coalesce, continuing the ascent as a single bubble with double the volume. Only when constant slag properties are used does the first coalesced bubble break into four main smaller structures, which rise up separately. The differences observed in the bubble shape, break-up, and therefore dynamics are evident and considerably higher than during the ascent of the single bubble, confirming the points made above. For each thermal gradient configuration, the bubble shapes change considerably along with the modeling approach of physical properties, which entails the establishment of different terminal velocities. This, together with the collision of the first two bubbles, enhances the break-up events at the trailing edges, which are present in all cases.

C. Effect of Bubble Diameter
The initial diameter of the released bubble is a key parameter in the study of bubble dynamics, since this and the fluid properties play a role in the determination of the terminal bubble shape and velocity. In the present    study, the variation of the diameter has an additional consequence, as it subjects the bubble to different ranges of temperature and therefore liquid properties. The simulations of two smaller diameters, namely D = 0.020 m and D = 0.015 m, were performed under the configuration of the radial thermal gradient, with the three modeling approaches for the slag's physical properties. Figure 25(a) shows the development of the gas-liquid interface over time for the three tested bubble diameters. The green curves represent the bubbles with D = 0.025 m, as already discussed in the previous section, while the red and blue curves are related to the smaller diameters. It is seen that break-up events do not appear when the initial diameter is reduced. Indeed, at D = 0.025 m smaller bubbles detach from the trailing edge when the experimentally measured slag properties are applied and the gas-liquid interface suddenly increases. This is not observed at D = 0.020 m and D = 0.015 m. Nevertheless, the development of the bubble shape follows the same trend at each diameter. Regardless of break-up, in a slag bath with constant properties the bubble assumes an oblate shape, whereas the usage of temperature-dependent properties leads to prolate bub-ble shapes. As already discussed above, a change in the bubble shape has an influence on the rising velocity, which Figure 25(b) clearly illustrates. At each diameter, the bubble calculated with experimentally measured properties rises with a higher velocity than the ones with literature models and constant properties. It follows that, despite a reduction in the effects, the usage of different modeling approaches for the slag properties leads to different bubble dynamics even for smaller initial diameters.

D. Effect of Gas in the Bubble
In all previous calculations, the gaseous phase consists of air at 1600 K. To evaluate the effect of the type of gas on the bubble dynamics, the configuration with a radial thermal gradient and experimentally measured slag properties was also performed with CO 2 and SO 2 as the gas phase, initially at the same temperature of 1600 K. These are gases that can typically be found in sulfide smelting processes and mostly differ, for the purpose of the present analysis, in terms of density (at room temperature  q air ¼ 1:225 kg/m 3 , q CO 2 ¼ 1:7878 kg/m 3 , q SO 2 ¼ 2:77 kg/m 3 ). Figures 26 and 27 show the results of the analysis regarding the bubble shape, gas-liquid interface, and terminal velocity. As expected, the dynamics of the bubble evolve in an identical manner, whatever the gas contained within. The gravitational force, responsible for buoyancy, is the only term where the density of the gas plays a role. Expressed as it is clear that considering the high densities of the slag, the consequences of having different gases inside the bubble does not affect its dynamics. Tomiyama's equation (Eq. [9]) shows the same for the determination of the terminal velocity.

VI. CONCLUSIONS
The dynamics of rising bubbles in a slag bath were investigated by means of CFD simulation, considering the heterogeneity of the liquid thermo-physical properties. The present work focused on the presence of temperature gradients as the cause for the variation in the slag properties. The main goal was to understand how different models of the temperature-dependent properties affect the solution of the bubble dynamics and, therefore, offer guidance on the approach used to simulate metallurgical multiphase systems involving bubbling flows. Three thermal gradient configurations were simulated, each of which involved three different modeling approaches for the slag physical properties: constant values, temperature-dependent data from the literature, and temperature-dependent data from measurements. The analysis was based on the observation of the bubble characteristics during its ascent, and the major outcomes can be briefly summarized as follows: The analysis of the bubble rising path shows that, deviations from the center line occur in the presence of transversal thermal gradients, when using temperature-dependent slag properties. On one hand, the application of data from the literature and experiments produces comparable results; on the other, the usage of constant values leads to a wrong solution. In the perspective of a reactor simulation, this would compromise the evaluation of the mixing inside the bath and therefore the mass and heat transfers.
Given the properties of the liquid, the bubble assumes a certain shape, which directly determines its dynamics. The shape and velocity are strictly linked, as also observed in the present study. The calculated bubbles are all of the ellipsoidal-cap type, but while the usage of constant properties produces an oblate profile, prolate shapes are obtained with temperature-dependent data. Furthermore, break-up events at the bubble skirt appear when the experimentally measured slag data are used. These are to be attributed to the increased unbalance between viscous stresses and surface tension, resulting from the higher viscosity and lower surface tension, in comparison with the literature-based models. All these aspects play a role in the determination of the bubble terminal velocity, which shows differences depending on the different modeling approaches. A correct calculation of the bubble rising velocity is crucial in the perspective of reactor modeling, because it directly affects the prediction of the bubbling frequency, often a key process control parameter.
The quantitative evaluation of the gas-liquid interface shows major discrepancies between the models. The oblate bubble shape obtained with the constant slag properties creates a higher surface-to-volume ratio than a prolate shape formed when considering the temperature dependency. In addition, as already discussed above, break-up phenomena take place when the measured slag data are used. This increases by up to 50 pct the gas-liquid interface in comparison with the case where temperature-dependent slag data from the literature are used. This marked difference would strongly compromise the calculation of hot and reactive bubbling flows, if interface heat and mass transfers were included.
Larger differences are observed when simulating two consecutive bubbles, since interaction between the two bubbles additionally takes place. The subsequent bubble is accelerated due to the depression created in the wake of the leading one. This causes the two to coalesce. The results are quite different between all the slag models, which confers more importance to the choice of an adequate model when simulating a rector process, since almost in each case bubble flows are established in the form of bubble chains.
A parameter analysis shows that varying the bubble diameter, and hence exposing the bubble to a different temperature range, still leads to differences between the slag modeling approaches. The presence of air, CO 2 or SO 2 as gas phase does not affect the bubble dynamics at all.
The analysis is limited to the temperature dependency of the slag's physical properties, but such considerations may also be extended to species and solid fraction gradients, which are present in pyrometallurgical slag baths. This investigation has shown that the choice of a specific slag modeling approach leads to markedly different results in terms of the bubble dynamics and may considerably affect the reliability of a rector model, where gas-liquid interactions occur. A correct evaluation of the interface area, for example, is essential for the determination of mass and heat transfer phenomena and, therefore, the kinetics of bath processes. Eventually, for the simulation of pyrometallurgical flows with interface-tracking methods, the author recommends the usage of slag property models which include temperature dependency, and, depending on the resources and aim of the project, to base them on experimental measurements.

ACKNOWLEDGMENTS
The authors would like to thank the Faculty of Mathematics and Computer Science of TU Bergakademie Freiberg and the Center for Information Services and High Performance Computing (ZIH) at TU Dresden for the allocation of the computing time. The German Federal Ministry of Education and Research (BMBF) has funded this research within the framework of the Center for Innovation Competence Virtuhcon (Virtuhcon II, Projects Number: 03Z22FN11, 03Z22FN12).

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativec ommons.org/licenses/by/4.0/.

FUNDING
Open Access funding provided by Projekt DEAL.