Bubbly Mold Flow in Continuous Casting: Comparison of Numerical Flow Simulations with Water Model Measurements

The impact on the mold flow of inert gas bubbles originating from a gas injection at the top of the submerged entry nozzle is investigated. Results from a physical model using water and air are compared with corresponding numerical flow simulation results. In numerical models, the bubble velocities are determined by calculating the force equilibrium between buoyancy, drag force, and other forces acting on the bubbles. For the observed bubble size in the physical model, the so‐determined rising velocity of the bubbles is significantly too high in comparison to the water model experiment. Various effects can influence the rising velocity of bubbles. One of them is that the presence of turbulence obviously reduces the rising velocity. The influence of turbulence models and of a turbulence‐induced bubble drag modification is analyzed in numerical flow simulations and compared to water model results.


Introduction
The injection of gas at the stopper rod tip in continuous casting of steel slabs is a common practice to prevent clogging and to control the produced steel quality. The injected gas forms bubbles and is transported by the liquid steel through the submerged entry nozzle in the mould region of the strand ( Figure 1). Most of the gas bubbles rise to the mold surface and leave the liquid steel there.
In a resting liquid, single gas bubbles start rising due to their buoyancy and reach quickly the terminal rising velocity where the buoyancy force and drag force are in equilibrium. Basically, the terminal rising velocity is higher for larger bubbles. Drag forces and resulting terminal rising velocities for bubbles in liquid have been intensively investigated for air bubbles in water, [1,2] but also for bubbles in liquid metals, [1] as shown in Figure 2a. The unstable shape of bubbles in contrast to rigid spheres [3] causes a saturation-like effect as the terminal rising velocity is only slightly increasing for bubbles larger than 2 mm. As a consequence, bubbles in the diameter range of 2-30 mm have a terminal rising velocity between 0.2 and 0.3 m s À1 , which is drastically lower than for rigid spheres with the same size and density. Nevertheless, drag laws for rigid particles are used in numerical simulations for gas bubbles. [4] The difference in surface tension of liquid metal and water obviously has a low influence on the terminal rising velocity. The bubble size in the mold flow in the real casting situation depends on the gas injection process and on the transport from the injection position [usually at the top of the submerged entry nozzle (SEN)] to the mold. The bubble's size and motion behavior in the real casting situation is hard to determine. Laboratory experiments of the flow situation are either downscaled in combination with liquid metals [5] or unscaled as well as downscaled (typically by a factor of 3) with water as the liquid. The model should represent the flow structure in the SEN and the mold. Similarity criteria can help to check if the same flow situation can be reached in the laboratory setup by choosing appropriate process parameters. The main differences between the real situation and the physical model are due to the geometrical scale (if not equal), the difference in surface tension, and the difference in absolute pressure levels. [6] Due to these differences, the bubble sizes may be different in the model in comparison to the real situation. As a consequence, it is difficult to ensure that a laboratory experiment is similar to the real casting situation.
In numerical simulations of the mold flow, the injection process is mainly not considered and therefore, a bubble size and density distribution at the inflow boundary condition in the SEN cross-section is required. In many cases a bubble size is chosen in combination with drag law to calculate the bubble motion. The bubble size is determined by either measuring the bubble size in an unscaled [7] or downsized water model; [8] observing or measuring the bubble distribution in a water model and choosing the bubble size in the simulation model that fits best to DOI: 10.1002/srin.202000415 The impact on the mold flow of inert gas bubbles originating from a gas injection at the top of the submerged entry nozzle is investigated. Results from a physical model using water and air are compared with corresponding numerical flow simulation results. In numerical models, the bubble velocities are determined by calculating the force equilibrium between buoyancy, drag force, and other forces acting on the bubbles. For the observed bubble size in the physical model, the sodetermined rising velocity of the bubbles is significantly too high in comparison to the water model experiment. Various effects can influence the rising velocity of bubbles. One of them is that the presence of turbulence obviously reduces the rising velocity. The influence of turbulence models and of a turbulence-induced bubble drag modification is analyzed in numerical flow simulations and compared to water model results.
the experimental situation; [9] choosing the bubble size according to the solidified bubbles found in the cast metal; choosing the bubble size in a way that "a reasonable effect" by the gas bubbles on the flow can be obtained; [10] or choosing the bubble size without any experimental validation. [4] Alternatively, in some simulations the bubble drift velocity is specified without calculating the drift velocity during the simulation. [7] The flow rate of the injected gas under casting conditions is usually specified in a volumetric flow rate at standard conditions (room temperature und atmospheric pressure). Due to the liquid steel temperature, the gas should expand approximately by a factor of 6 if it has the same temperature as the liquid steel (1800 K/300 K ¼ 6). The hydrostatic pressure in the liquid steel inside of the SEN in the height of the mold level is approximately equal to the atmospheric pressure; 30 cm below the mould level (a typical penetration depth for gas bubbles with some millimeter diameter in mold flows) the expansion factor is %5, while it is around 12 at 700 mm distance above the mould level (typical position for gas injection inside the top region of the SEN). Some authors use an expansion factor of 6, [7] some use a lower factor, e.g., 4.1, [11] and some authors do not correlate their model injection rates to real casting injection rates.
The motion of bubbles in turbulent flows seems to be more complex than in a resting liquid: As reported in ref. [12], turbulence is observed to reduce the bubble rising velocity, [13,14] to have no influence on the bubble drift velocity, [15,16] or to even increase the rising velocity. [17] The mechanisms of these   [1] calculated with drag law [2] and drag modification [13,14] for different Kolmogorov length scales λ and for particles (rigid spheres); [3] b) bubble drift velocity under turbulent flow conditions a function of the turbulent dissipation rate ϵ for different bubble sizes calculated with drag-law [2] and drag modification. [13,14] www.advancedsciencenews.com observed phenomena are only partially understood. Therefore, mathematical models to describe these effects are often based on empirical relations. A reduced bubble drift velocity due to turbulence in a mold-like flow was observed by Kwon et al. [18] and the model of Brucato [14] (see also Figure 2a) was applied to consider the effect in the numerical flow simulation.
Experiments show that the presence of bubbles in turbulent flows also increases the turbulence. [19] Sato's model [20] adds an extra eddy viscosity term proportional to the bubble volume fraction and drift velocity and is often applied to bubbly mold flow simulations. [9,[21][22][23] Although the impact of the bubble size on the bubble drift velocity is weak over a large diameter range, some investigations consider different size classes of bubbles and their exchange due to coalescence and breakup. [21,22] 2. Mathematical Modeling

Turbulence Modeling
In many simulations, even very recent ones, the standard k-ε turbulence model is used for the turbulence modeling of the liquid phase or of the mixture of liquid and gas. [7,9,21,22] Recently some groups have been using large eddy simulations for multiphase mold flow simulations, [8,23,[25][26][27] while hybrid turbulence models such as the scale adaptive simulation (SAS) model are very rare. [23] For the results presented in this publication, two turbulence models are used: the realizable k-ε [28,29] (RKE) and the SAS model. [29][30][31] Both models are used to calculate the transient single phase flow in the mold water model described in the "Experimental Section" for the comparison with single phase flow measurements. Therefore, the geometrical shape of the mold is equivalent to the water model (dimensions see Table 1) and, for example, solidification (similarly to ref. [32]) is not considered. Figure 3 shows the grid used for both simulations with a cell size of 5 mm in the interior of the SEN and with boundary layers of 0.3 mm at the SEN interior walls, at the wide and narrow faces of the mold, and at the top free surface of the mold, which is modeled as a planar wall with zero shear stress. The RKE model is used in combination with the "enhanced wall treatment" model, [29,33,34] which allows the correct treatment of the turbulent boundary layer for resolved and nonresolved boundary layers. Therefore, the dimensionless wall distance of wall-adjacent cells (the so-called y þ value) can be lower and higher than 1, in contrast to the classical wall treatment, where y þ must be higher than 1 to ensure that the wall-adjacent cell lies in the logarithmic boundary layer. The SAS model is based on k-ω turbulence equations, so there is no special wall treatment independent of the boundary layer resolution. The y þ values reached in the numerical simulations are between 3 and 10 at the inner walls of the SEN and lower than 8 at the mold walls.
For both simulations a zero velocity field is used as an initial condition and 100 s are calculated. The time-step size for the RKE model is 0.01 and 0.001 s for the SAS model. The time-step chosen for the SAS model ensures that the fluid needs 2.5 time-steps to pass a 5 mm finite volume cell at a typical flow velocity of 2 m s À1 . This is necessary to properly consider the propagation of vortices resolved by the grid without numerical diffusion. In contrast, the time-step size for the RKE model can be 10 times larger as the higher turbulent viscosity used by this model prevents the self-induced formation of vortices even if small time-steps are used. As a consequence, the computational effort is much higher for the SAS than for the RKE model, but turbulence can be better modeled with the SAS model as will be shown in the results.

Gas Bubbles
For the simulation of the bubbles a Eulerian approach is used to calculate the bubble distribution by solving a transport equation for the bubble volume fraction field where α is the gas volume fraction, u x , u y , and u z are the liquid velocity components, v is the relative velocity between the liquid and gas bubble, μ t is the turbulent viscosity, Sc ¼ 0.7 is the turbulent Schmidt number, and S is a source term used for the boundary condition in grid cells adjacent to the top surface, where bubbles escape through the boundary. Gravity is assumed to act in the negative z-direction. The gas volume fraction at the  inlet boundary inside the SEN is calculated from the relation of the gas and liquid volumetric flow rates, assuming equal velocities for both phases and a homogeneous bubble distribution. The terminal rising velocity of a bubble in a liquid without turbulence where d is the bubble diameter, g is the gravitational acceleration, Δρ is the density difference between liquid and gas, ρ is the liquid density (for water/air and argon/steel Δρ ρ % 1Þ, C D ðReÞ is the drag coefficient as a function of the particle Reynolds number Re, and ν is the kinematic viscosity of the liquid. This equation is an implicit equation as v can be found at the left hand side of the equation as well as on the right hand side as Re is a function of v.
Brucato et al. [13] performed measurements of particle sedimentation with different particle sizes and densities in a liquid, where they could show that the sedimentation velocity decreases if the turbulence increases, induced by a cylinder rotating about its vertical axis. They proposed a drag coefficient modification that depends on the turbulent length scale, assuming a homogeneous turbulence over their experimental setup. The experimental results fit well to their proposed drag reduction formula. Lane et al. [14] performed simulations of a stirred tank containing a liquid with gas bubbles with inhomogeneous turbulence and compared the simulation results to corresponding measurements of the bubble distribution. They found out that Brucato's relation only gives results closer to their measurements if they significantly reduce the numerical factor in Brucato's drag reduction equation from 8.76 Â 10 À4 to 6.5 Â 10 À6 where C D,0 is the drag coefficient without any turbulence, λ is the Kolmogorov length scale, and ϵ is the turbulent dissipation rate. The observation that the numerical coefficient in Equation (3) has to be adopted to each flow situation implies that there are probably additional flow or turbulence quantities influencing the drag reduction not considered in the equation. However, every turbulence model (except direct numerical simulations) is somehow based on empirical relations and using Equation (3) with an adopted coefficient and the local length scale leads to significantly better results, as will be shown in the results section. Inserting Equation (3) in Equation (2) and neglecting the change of the drag coefficient due to a change of the particle Reynolds number due to the change of the drift velocity, the bubble rising velocity in a turbulent liquid becomes where v 0 is the rising velocity in a liquid without turbulence. A numerical evaluation of the exact and simplified equation using Tomiyama's drag law [2] showed a good agreement of the simplified and exact solution over a wider range of typical bubble sizes and turbulent dissipation rates (Figure 2b). For the simulation results, v 0 ¼ 0.23 m s À1 and d ¼ 5 mm are used.
The bubble size is a typical size observed in the water model experiments ( Figure 5), and the terminal rising velocity v 0 for this size is determined using the diagram in Figure 2a. As the terminal rising velocity of bubbles with this size varies only weakly with the bubble diameter (refer to Figure 2a), the consideration of different bubble diameters and of a bubble size distribution is expected to have a negligible influence and therefore, a uniform bubble size is chosen. Similar values for bubble size and rising velocity were also derived from experiments in ref. [7]. In a later phase of the project, the numerical model will be adopted to simulate the real situation with liquid steel, and if necessary, different bubble sizes will be considered.
In the simulations with the SAS turbulence model the dissipation rate is calculated by ϵ ¼ 0.09 kω according to the definition of ω. [29] In a bubbly liquid flow, gas bubbles displace liquid. If the same liquid flow rate flows through the SEN without and with 10% volume fraction, for example, and if the gas bubbles are homogeneously distributed in the SEN, the liquid velocity has to increase by 10% due to the gas volume displacement. This effect is neglected in most simulations using Lagrangian particle trajectory calculations (see e.g., ref. [18,25]; 17 out of 31 analyzed publications on numerical simulation of mold flows with gas bubbles use Lagrangian particle trajectory calculations). In the presented simulation results, a constant mixture density is used and therefore, the gas volume displacement is neglected.

Experimental Section
At voestalpine a 1:1 scaled water model of a tundish and strand is available as sketched in Figure 4. Table 1 lists the geometrical and process parameters. A stopper rod was used to control the flow through the SEN. Air can be injected at the stopper rod tip to simulate the argon injection. Figure 5 shows a close-up of the SEN port region. In this picture it can be seen that the bubble size is relatively uniform around 5-6 mm. As mentioned previously, the bubble size in the water model may differ from the bubble size in the real situation, but here only the situation in the water model was considered in the numerical simulation as well as a validation for the numerical simulation.
A video camera was positioned in front of the mold water model (Figure 4). The bubble distribution can be qualitatively seen in the video. A long time average shows the time-averaged distribution; see Results.
Two kind of particle image velocimetry (PIV) measurements were performed. The first used the classical PIV method for the single phase flow (no gas injection) with a particulate seeding, a pulsed laser light sheet, and a PIV camera to measure the 2D velocity field in the right half of the mold center plane ( Figure 4). The multiphase PIV method turned out to be inapplicable for typical gas injection rates because the bubble density was so high that they obstructed the camera view on the center plane. [24] For the second type of PIV measurement, two PIV cameras were located above the top surface of the mold (Figure 4). The bubbles on the surface were illuminated by a linear flash light synchronized with the PIV cameras and used as tracer particles for the PIV algorithm to measure the flow field at the top surface in the multiphase situation. For typical gas injection rates the bubble density was high enough so that the tracer density was sufficiently high and bubbles below the surface were well covered by the bubbles on the surface.

Results
In the following comparisons of measured und simulated results, time averaging is used for time series where u i is a time series of i ¼ 1 .. n values like a velocity component or the camera image intensity at a certain position, andū is the corresponding time average value at this position. Figure 6 shows the time-averaged flow fields in the interior of the SEN at the center plane parallel to the narrow faces in the region of the SEN outlet ports by velocity vectors and contours of the velocity magnitude and the turbulent kinetic energy: For the PIV measurement, the turbulent kinetic energy is estimated by evaluating the two velocity fluctuation components in the measuring plane and multiplying it by 3/2 to account for the third (not measured) component perpendicular to the measuring plane [35] k ¼

Single Phase Flow
where u 0 x , u 0 y , u 0 z are the velocity fluctuation components (actual value minus time average) in a local coordinate system where the x,y-coordinates are aligned with the measuring plane and the z-coordinate is perpendicular to the measuring plane. For the RKE result, an instantaneous value of the turbulent kinetic energy as calculated by the turbulence model is used. For the SAS model, the turbulent kinetic energy is calculated from the resolved fluctuations of the velocity fields analogous to the PIV measurement according to Equation (6). In the considered region, the turbulent kinetic energy calculated by the SAS model equations only represents the nonresolved turbulence and is much lower than the turbulent kinetic energy calculated from the velocity fluctuations. The PIV measurements show a strongly asymmetrically fluctuating flow pattern, but a symmetrical double vortex structure in the time-averaged flow field (Figure 6a) with high turbulent energy in the vortex centers. Only the SAS turbulence model is able to capture this flow structure, which significantly influences the free jets leaving the SEN ports and, therefore, the mold flow structure. The downstream velocity magnitude in the upper region of the displayed area is lower in the SAS simulations than in the measured PIV flow field. In contrast, the velocities are similar in the upstream regions near the side walls of the SEN. The magnitude of the maximum turbulent kinetic energy is only about 70% in the SAS simulation result in comparison to the measured PIV result. This difference seems acceptable due to the insufficiencies in the determination of the turbulent kinetic energy and due to the fact that the turbulent kinetic energy of the SAS model in the bottom area of the SEN results mainly from the flow velocity fluctuations induced by the self-instability of the flow in this region and is nearly independent of the inlet turbulence boundary condition located in the SEN cross-section.
In Figure 7 the time-averaged flow field in the right half of the vertical mold center plane parallel to the wide faces is outlined by velocity vectors and velocity magnitude contours obtained from a PIV measurement (300 measurements during 75 s) and from numerical simulations using the RKE and the SAS turbulence model. The PIV method has a limited velocity range, depending on the time distance between two subsequent camera images, which determines the displacement of the tracer particles between two images, used for the local velocity calculation. In regions of high velocities, the displacement gets too high, and the PIV correlation algorithm fails. On the other hand, if displacements become too low in regions of low velocities, the calculation of the velocities becomes inaccurate. The image time distance was chosen so that the mold velocities could be properly resolved. Therefore, the velocities in the jet core are too high and the displayed values are interpolated values from the surrounding velocities. As a consequence, the velocity magnitude can only be compared to the numerical simulations outside of the jet core region.
The following criteria are used to compare the accuracy of the numerical models in comparison to the measured data. 1) A dashed line starting at the upper port edge marks the center of the free jet defined by the maximum jet velocity of the PIV data. This line is also drawn for the numerical results at the same geometrical position. The jet center resulting from the SAS model velocity field is much closer to this line than the RKE. 2) Near the impingement area of the jet at the narrow face, the velocities decay too much in the RKE results, while the SAS result is closer to the measured data. 3) The maximum While simple turbulence models such as the RKE are very popular as the computational effort is much lower than for the SAS model, the results obtained by the SAS model are much closer to the PIV measurement in comparison to the RKE model. The proper reproduction of the free jet characteristics and the top surface velocity is essential for the multiphase situation as they determine the propagation of the gas bubbles in the mold. Figure 8 shows instantaneous snapshots from the gas bubble distribution in the water model and from corresponding transient multiphase simulations using the RKE and the SAS turbulence model with and without the Brucato drag modification for a gas injection rate of 40 l min À1 , corresponding to 6.7 l min À1 argon injection rate at standard conditions in liquid steel if the thermal expansion of the gas is considered by a factor of 6, as described previously. The color camera image is converted to a grayscale image and then colored with a color map like the one used for the numerical results (dark ¼ blue, intermediate brightness ¼ green, red ¼ maximum brightness). In the camera image, bubbles appear as bright spots, while the dark background can be seen in areas without bubbles. Basically, brighter areas of the image qualitatively indicate the presence of bubbles in these areas somewhere over the whole mold thickness. As the lighting situation is not homogeneous and the camera image brightness is determined by bubbles over the whole mold thickness (not only in the center plane), the image brightness is not quantitatively related to the volume fraction in the center plane, but can be qualitatively compared to get an impression of the bubble distribution in the mold.

Two-Phase Flow with Gas Injection
For the simulation results, the gas bubble volume fraction in the center plane is displayed. Therefore, the water model and simulation results can only be compared qualitatively. While the SAS model resolves the inhomogeneous bubble distribution as observed in the water model, the RKE model produces a significantly smoother gas distribution. Without the Brucato drag modification, the gas bubbles rise too fast near the SEN and therefore, the bubble distribution remains much too concentrated around the SEN in comparison to the water model observation. Figure 8c,e shows the bubble drift velocity resulting from the Brucato drag modification. In regions of higher turbulence, the drift velocity decreases. For the simulation results without the Brucato drag modification, the bubble drift velocity is not shown as it is constant over the whole domain.
As a snapshot of a randomly fluctuating process is arbitrary and therefore difficult to compare, time-averaged results are expected to be more meaningful. Thus, a 120 s long timeaveraged camera image is compared to time-averaged gas distributions from simulations, again with the RKE and SAS model and each with and without the Brucato drag modification, respectively. Again, the numerical simulation results with the Brucato drag modification fit much better to the water-model camera picture than those without drag modification. Figure 9 shows results similar to Figure 8 for different gas injection rates (32,40,60, and 100 L min À1 ). Only time-averaged simulation results calculated with the SAS turbulence model and the Brucato drag modification are shown in comparison to the time-averaged camera pictures taken at the mold water model. The gas distribution changes with the injection rates, which seems qualitatively well reproduced by the simulations. The slightly asymmetric left-right color distribution in the time-averaged camera images is probably due to an asymmetric lighting situation.
For a quantitative comparison of measured and calculated flow structures under the influence of gas injection, the top surface velocity fields are considered. Figure 10 shows a comparison of measured and calculated time-averaged flow velocities on the top surface of the mold for different gas flow rates between 0 and 100 L min À1 . The colors in the contour plots represent the    velocity component parallel to the wide face of the mold in the direction of the mold center u x · signumðÀxÞ, where x is the coordinate axis parallel to the wide face of the mold. Starting with low gas injection rates (upper rows), the PIV measurements (not available without any gas as there is no tracer for the PIV method in this case) show a typical double roll flow pattern with flow velocities directed to the mold center (SEN), as indicated by the dark blue color. With increasing gas injection rates (lower rows), the magnitude of the flow velocities directed to the mold center becomes lower, until the flow pattern changes to a single roll pattern with large areas where the flow velocities are directed from the mold center toward the narrow faces, as indicated by the red color. While this transition between double and single roll flow pattern is qualitatively reproduced by the SAS turbulence model simulations (with Brucato drag modification), the flow pattern calculated with the RKE turbulence model (also with Brucato drag modification) remains a double roll flow pattern even for the highest gas flow rates.
Although there is a good agreement of the mold top surface velocity without gas injection between PIV measurement und SAS simulation as demonstrated previously, there is a nonnegligible difference for even small gas injection rates. According to the PIV measurements in Figure 10, the top surface maximum velocity for 16 L min À1 gas injection is only 0.24 m s À1 , while the maximum calculated by the SAS turbulence is still 0.3 m s À1 like in the single-phase flow situation. The gas injection seems to reduce the top surfaces' velocities even for small gas rates. For higher gas flow rates, the numerical simulation results continue to reach a similar flow pattern and velocity magnitude somehow delayed (with respect to the gas injection rate) in comparison to the water model results. One reason could be that the bubbles induce additional turbulence to the liquid flow as observed, for example, in ref. [24]. The model of Sato [20] will be used in future simulations to consider this effect.
In Figure 11 the total gas volume in the mold ∫ αdV and the kinetic energy of the flow ∫ 1 2 ρu 2 dV are shown as a function of the gas injection rate for the RKE and the SAS turbulence model. The kinetic energy of the flow field calculated by the RKE is consequently higher than the kinetic energy of the flow field calculated by the SAS model, resulting, for example, in overestimated top surface velocities, as shown previously in comparison with measurements. Obviously damping effects in the free jets exiting the SEN ports are better resolved by the SAS model. The choice of the turbulence model also influences the gas holdup for higher injection rates.

Conclusion
Water model results are compared with corresponding numerical flow simulations using the RKE and SAS turbulence models for the single-phase situation without gas injection. For the considered SEN and mold geometry, essential characteristics of the turbulent mold flow such as the SEN jet flow pattern and top surface velocities can only be captured by the SAS turbulence model, which resolves coarser turbulence structures and is therefore about ten times more time consuming than the RKE model. For the two-phase situation with argon injection, the Brucato model is used to consider the reduced bubble drift velocity in regions with higher turbulence, yielding gas bubble distributions in the mold that are much closer to the water model than without the Brucato model (turbulence independent bubble drift velocity).
For higher gas injection rates, the top surface flow structure changes from double toward single roll, which can be observed in both PIV measurements and SAS results. In contrast, the RKE results cannot satisfyingly reproduce the transition from double to single roll, as the flow structure remains double-roll-like even for very high gas injection rates. Therefore, k-ε turbulence models such as the RKE do not seem to be appropriate for both single-and two-phase mold flow simulations.
In future investigations, we will try to bring numerical simulation results and water model measurements closer together by further enhanced measurements and numerical simulation models. The validation has to be performed over the whole relevant process parameter range for gas injection rates, casting speed, SEN immersion depths, and mold dimensions to guarantee the robustness of the numerical simulation model. A grid independence test similar to that in ref. [36] is also planned to ensure that the used grid is fine enough and the results do not depend on the grid cell size. Finally, the numerical simulation model will be transferred to the real process situation (liquid steel and argon injection, including solidification effect).