Transient 3D hydrodynamic model of a blast furnace main trough

Abstract

A transient 3D CFD model is solved to investigate the features of the complex flow in a blast furnace trough. Hot metal, slag, and air are considered as different phases of the flow. The influence of various taphole stream conditions on the hydrodynamics in the trough is examined, such as different slag-hot metal ratios, taphole diameters, and stream velocities. The special case of a dry trough during its first tapping is also addressed. Attention is devoted to the characterization of the wall shear stress, closely related to mechanical erosion, and to the evolution of the interfaces separating the fluid phases. Intricate flow patterns, characterized by large recirculations and return currents of both slag and hot metal, are observed. The simultaneous tapping of slag and hot metal leads to distinct flow features, deviating from the initial stage, when only hot metal is drained, as well as lower shear stress values, up to 31% less with increasing slag content in the taphole stream. Slag and hot metal levels in the trough evolve rapidly to quasi-steady states. Although the influence of the taphole stream velocity and diameter have a modest impact on the free surface dynamics, variations in the taphole slag fraction lead to significant fluctuations in the depth of slag and hot metal pools in the trough. The height of the slag-iron free surface can increase by up to 40% relative to its initial level, while the variation in the slag pool depth approaches 30%. These changes in the interface positions carry potential implications for refractory wear, as they suggest a substantial shift of the slag-hot metal line during each tapping.

Keywords:

• CFD
• steelmaking
• blast furnace trough
• multiphase flow
• wall shear stress

1. Introduction

The manufacturing of steel is mainly carried out using an intermediate material, known as hot metal or pig iron, produced in large metallurgical furnaces called Blast Furnaces (BF) by reducing iron ore. They are responsible for around 70% of global steel production (Geerdes et al., 2015). In the foreseeable future, the BF is expected to remain the main source of steel production worldwide, as no other steelmaking process can achieve the efficiency of modern BFs (Cameron et al., 2019). To maintain the economy of scale of the process, its productivity and campaign life must satisfy an ever-increasing demand. This places greater requirements on the refractory materials used to line steelmaking vessels, as larger BF throughputs lead to higher rates of mechanical erosion, chemical attack, and also bigger thermal stresses.

Nomenclature

Greek symbols

α	=	Volume fraction [-]
${\gamma }_{tap}$	=	Taphole angle [${}^{\circ }$]
Γ	=	Boundary of Ω
μ	=	Dynamic viscosity [$\mathrm{P}\mathrm{a}\cdot \mathrm{s}$]
${\mu }_T$	=	Eddy viscosity [$\mathrm{P}\mathrm{a}\cdot \mathrm{s}$]
ξ	=	Characteristic function [-]
ρ	=	Density [$\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$]
σ	=	Surface tension [$\mathrm{N}/\mathrm{m}$]
${\tau }_w$	=	Wall shear stress [$\mathrm{P}\mathrm{a}$]
ω	=	Specific dissipation rate [$1/\mathrm{s}$]
Ω	=	Computational domain

Latin symbols

$a_1$	=	Constant of the turbulence model
$C_{\mu }$	=	Constant of the turbulence model
$C^i$	=	ith cell in the grid of the computational domain
$\mathbf{D}\left(\mathbf{v}\right)$	=	Symmetric part of the velocity gradient [$1/\mathrm{s}$]
$D_{tap}$	=	Taphole diameter [$\mathrm{m}\mathrm{m}$]
$F_1,F_2$	=	Blending functions of the SST $K-\omega$ turbulence model
$\mathbf{g}$	=	Gravitational acceleration [$\mathrm{m}/{\mathrm{s}}^2$]
$h_B$	=	Bottom of the trough height [$\mathrm{m}$]
$h_{hm}$	=	Hot metal-slag interface height [$\mathrm{m}$]
$h_{hm}^{\ast }$	=	Hot metal-air interface height, downstream [$\mathrm{m}$]
$h_{jet}$	=	Height of the jet [$\mathrm{m}$]
$h_{skimmer}$	=	Bottom of the skimmer block height [$\mathrm{m}$]
$h_{slag}$	=	Slag-air free surface height [$\mathrm{m}$]
$H_{hm}$	=	Hot metal pool depth (upstream of the skimmer) [$\mathrm{m}$]
$H_{hm}^{\ast }$	=	Hot metal pool depth (downstream of the skimmer) [$\mathrm{m}$]
$H_{slag}$	=	Slag pool depth [$\mathrm{m}$]
$\mathbf{I}$	=	Identity tensor
$I_T$	=	Turbulence intensity [-]
K	=	Turbulent kinetic energy [${\mathrm{m}}^2/{\mathrm{s}}^2$]
$l_{mix}$	=	Mixing length [$\mathrm{m}$]
$\mathcal{L}$	=	Length scale [$\mathrm{m}$]
$\dot{m}$	=	Mass flow rate [$\mathrm{k}\mathrm{g}/\mathrm{s}$]
${\mathcal{M}}_{\mathcal{T}}\left(\varphi \right)$	=	Time average of the variable ϕ over the time scale $\mathcal{T}$
$\mathbf{n}$	=	Outward-pointing unit normal vector [-]
$N_{cell}$	=	Number of cells in the grid of the computational domain
p	=	Pressure [$\mathrm{P}\mathrm{a}$]
Pr	=	Prandtl number
Q	=	Volumetric flow rate [${\mathrm{m}}^3/min$]
Re	=	Reynolds number
S	=	Strain rate magnitude [$1/\mathrm{s}$]
t	=	Time [s]
$t_{end}$	=	Final time value [s]
$\mathcal{T}$	=	Time scale [s]
$\mathbf{v}$	=	Velocity [$\mathrm{m}/\mathrm{s}$]
${\mathbf{v}}^0$	=	Initial velocity of the jet [$\mathrm{m}/\mathrm{s}$]
${\mathbf{v}}_{inlet}$	=	Inlet velocity [$\mathrm{m}/\mathrm{s}$]
$v_{tap}$	=	Taphole stream velocity [$\mathrm{m}/\mathrm{s}$]
$\mathcal{V}$	=	Velocity scale [$\mathrm{m}/\mathrm{s}$]
$\mathbf{x}$	=	Vector of spatial coordinates [m]
${\mathbf{x}}^0$	=	Initial point in the jet trajectory [m]
$y_w$	=	Normal wall distance [m]

Subscripts

${}_{air}$	=	Air
${}_{hm}$	=	Hot metal
${}_{slag}$	=	Slag
${}_{tap}$	=	Taphole
${}_{wall}$	=	Wall

Superscripts

${}^0$	=	Initial value
${}^{imp}$	=	Impact point of the taphole jet

Acronyms

BF	=	Blast Furnace
BFT	=	Blast Furnace Trough
CFD	=	Computational Fluid Dynamics
LES	=	Large Eddy Simulation
PLIC	=	Piecewise Linear Interface Construction
RANS	=	Reynolds-Averaged Navier-Stokes
SST	=	Shear Stress Transport
VOF	=	Volume of Fluid

Therefore, it is of interest to assess and minimize the wear suffered by the refractory materials used in the iron and steel industry. Most of the experimental and numerical studies available in the literature focus on different parts of the BF itself. However, another important vessel used in ironmaking is the Blast Furnace Trough (BFT) or main trough, a trench-like receptacle lined with castable refractories wherein the products of the blast furnace –molten hot metal and slag– are collected after their extraction from the BF and separated by density difference. Avoiding trough breakage and increasing its availability is a major concern for the metallurgical industry (Kumar et al., 2010), as the refractory linings suffer a high degree of degradation. Consequently, they are subject to periodic repairs and operation issues with the consequent costs for the steelmaking companies. To ensure a safe operation and prevent breakouts, the trough campaign life is usually capped after a certain hot metal tonnage is reached. After that, the trough is repaired and the most external refractory layer, which is in contact with the hot fluids, known as working lining, is replaced. Thus, accurately assessing trough refractory wear could increase operation safety, as well as productivity, by reducing the frequency of repairs.

In the BF trough, several coupled phenomena contribute to refractory wear and tear. First, the flow patterns in the trough are complex and intricate. The discharge of the fluids extracted from the BF generates a molten jet stream, which upon its impact onto the BFT generates a highly turbulent mixed gas-liquid multiphase flow with large amounts of entrained air bubbles (Begnis & Brandaleze, 2006). Mechanical erosion, which is closely related to wall shear stress on the jet impingement zone, plays a significant role in refractory wear. Furthermore, as the BF draining process is not continuous and the fluid discharge conditions are not steady, the structure is subject to cyclic thermal loading during its entire campaign life, reaching high temperatures, around 1500${}^{\circ }$C (see e.g. Barral et al., 2022). Lastly, chemical corrosion, enhanced to some extent by the interaction with hot metal and slag, also affects refractory performance. The precise mechanisms that lead to lining degradation are not fully understood yet. Nonetheless, it is generally accepted that the wear of the trough refractories is the consequence of the combined action of the fluid flow, chemical reactions, and heat and mass transfer phenomena (Ranjan et al., 2022).

The main motivation of this work is to investigate the complicated flow that is generated by the impact of the molten jet stream exiting the taphole on the liquid pool held in the trough, as well as the generated shear stress in the refractory linings. These flow patterns have been a matter of research by several authors, as they are believed to have a strong influence on the refractory lining wear. Some publications have relied on physical scale models using oil and water to emulate slag and hot metal, while other researchers have used Computational Fluid Dynamics (CFD) to solve the Navier-Stokes equations with the purpose of characterizing the multiphase flow. Kim and Ozturk (1998) used a physical scale model of the BFT to study the influence of several process variables in slag and hot metal separation, such as the BF discharge rate or the trough geometry. On the other hand, He et al. (2002) used a similar model to study the general flow features. Phase separation, as well as the main characteristics of the flow, were also assessed with a scale model by Luomala et al. (2001). It was concluded that return flows in the trough, generated by the tapping stream, strongly influence the wear of the refractory wall. The effect of longitudinal pads and other modifications in the flow were analyzed. In Begnis and Brandaleze (2006), a physical scale model was developed to investigate the flow patterns. They also inspected the damage on a trough between campaigns, finding that the most severely worn zone correlates well with the higher-turbulence regions resulting from the jet impact.

The hostile environment in the BFT, posing a significant challenge to experimentation, has motivated the surge of mathematical modelling and numerical simulation applied to the main trough, especially in recent years. Rezende et al. (2007) considered a CFD model, involving the multiphase flow of hot metal and air in the trough, and using an algebraic turbulence model. The features of the flow resulting from the jet stream impingement were studied, showing qualitative agreement among the numerical results and those reported by physical models. Moreover, Wang et al. (2017) considered a similar mathematical model to investigate the flow in the trough, including slag as a different fluid phase and buoyancy effects due to thermal expansion of the fluids. Separation efficiency of slag and metal was analyzed, concluding that controlling the height level between the slag runner and the iron dam is essential to achieve separation. In Kou et al. (2019), separation was also addressed with a CFD model. Geometry changes, such as decreasing the trough bottom slope, were shown to have a sizeable impact on the slag-metal separation. Also, Monteiro de Oliveira et al. (2020) combined a CFD model with a physical scale model, which was used to undertake extensive experiments and validate the numerical findings. The main focus was also to study separation efficiency, which was found to be directly related to the residence time in the BFT. Flow modifiers were introduced to contain the region where the mixing is more intense, achieving the typical features of plug flow in a greater region of the trough. Ge et al. (2020) and Yao et al. (2021), followed an analogous approach, modelling the flow in the trough as a single-phase flow in a fixed domain instead. The former considered a mixture of hot metal and slag, while the latter used cast iron properties for the fluid. In both studies, conjugate heat transfer with the refractory lining was incorporated, and Yao et al. (2021) also assessed the thermal stresses in the solids, considering that the refractory behaves as a linear thermoelastic material. Other studies have addressed heat transfer in the BFT considering mathematical models defined on plane cross-sections (e.g. Barral et al., 2022; Li et al., 2023; Vázquez-Fernández et al., 2019), neglecting the fluid flow.

As discussed above, there is a significant body of research regarding the flow conditions in the trough. However, most CFD studies do not include the thick slag layer above the hot metal and, instead, the trough hydrodynamics are modeled considering a two-phase flow (air and hot metal) or a single-phase flow. Those studies that include slag and hot metal as different phases mainly address the issue of slag entrainment with hot metal (e.g. Kou et al., 2019; Wang et al., 2017). The influence of different slag-iron ratios in the taphole stream has not been carefully analyzed. Furthermore, wall shear has not been thoroughly studied in the literature, and large differences between the reported results by different researchers are observed (Ranjan et al., 2022).

In previous work, a thermo-hydrodynamic 3D CFD model was solved to find the position of the critical isotherms in the trough refractory linings (Barral et al., 2019). However, the model only included the latter part of the trough and significantly idealized the fluid flow. For instance, fixed velocity inlet conditions based on hot metal and slag production rates were specified in the middle of the trough. Also, a numerical simulation of the impact of the hot metal jet falling from the blast furnace on the BFT was presented, considering a 2D jet and the Wilcox $K-\omega$ turbulence model. In a subsequent study (Barral et al., 2021), a comparative analysis of the former model with the Shear Stress Transport (SST) $K-\omega$ turbulence model was conducted in two benchmark tests. The SST model was found to be more accurate at simulating the shear stress produced by a jet on a wall. In addition, the mechanical impact of a planar 2D jet on the trough was investigated using this model considering a dry trough.

This study aims to examine the influence of the slag pool on the flow dynamics in the BFT, devoting special attention to the evolution of the positions of the interfaces separating the fluid phases and the wall shear stress, which have not been covered comprehensively in the existing literature. To this end, a realistic geometric model of the trough is considered, and the impact of the 3D jet on the trough is analyzed, dropping some of the idealizations that were introduced in previous works. Different taphole stream conditions are investigated, covering various stages that occur during the BF drainage process. First, the effect of the slag fraction in the taphole stream is investigated. The specific situation during the first tapping in the campaign cycle of the BFT is also addressed, when there is no liquid pool to mitigate the mechanical impact of the hot metal jet. Subsequently, since the taphole diameter and stream velocity vary dynamically during the BF drainage, the influence of these parameters is examined.

Provided that it is not computationally feasible to solve a detailed hydrodynamic model for a full BF tapping, we set a time interval that leads to a solution that can be considered quasi-steady. Under these circumstances, for the given model conditions, the variations of the variables describing the flow, such as the pressure or velocity, can be regarded as fluctuations around a steady mean value. With this purpose, several simplifying assumptions are introduced. Namely, an incompressible and isothermal flow is considered, and the air circulating above the pool of fluids is assumed to be still. Also, surface tension effects are neglected, and the involved fluid phases are regarded as immiscible and continuous. The specific details of the industrial process and the features of the BF trough modelled in this work are described in Section 2. In Section 3, the computational domain, the model equations as well as the boundary and initial conditions that are used to close the mathematical model are presented in detail. Section 4 is devoted to the discussion of the numerical procedure. In Section 5, a grid sensitivity study is presented, and the numerical results for various taphole conditions are examined. Lastly, in Section 6, some concluding remarks and future research lines are addressed.

2. Physical problem

The removal of the accumulated molten slag and hot metal in the BF hearth takes place through the taphole, which is a constructed opening in the hearth wall. This drainage process is frequently known as tapping, often simply referred to as tap or cast. While the taphole is open, the exiting molten stream of hot metal and slag falls onto the BFT, where they undergo separation based on differences in density.

The separation process initiates within the hearth even before the liquids are tapped, given the significant density contrast between slag and hot metal, typically around $2600\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ and $7000\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$, respectively. Since the hot metal and slag are stratified inside the BF hearth, during the beginning of each cast usually only the former is drained, as the metal-slag interface is above the taphole. As iron is removed from the BF, the interface separating it from slag tilts down in the vicinity of the taphole, eventually allowing slag to also be drained (see e.g. Shao & Saxén, 2013b). Thereafter, slag and hot metal are extracted simultaneously until BF gas bursts out from the taphole, point at which it is closed by forcing a plug of refractory clay into the orifice, sealing it, and ending the tapping.

The tapping stream exiting the taphole is driven by the furnace pressure, with velocities usually above 5 m/s. It behaves as an unconfined fluid jet, which travels through air describing a parabolic profile. Ideally, it should have a low arch, with minimal splashing to avoid excessive damage in the trough refractories. When slag and hot metal are tapped simultaneously, the flow regime in the taphole can become more complicated. Some researchers have implicitly assumed a fully dispersed flow of slag and hot metal in the taphole (Iida et al., 2008). Stevenson and He (2005) found that slug flow, i.e. a partially dispersed regime of slag and hot metal, is likely to occur with gas-liquid experiments; characterizing the dynamic behaviour of this type of flow can be challenging, as highlighted in studies such as He et al. (2022). On the other hand, Shao and Saxén (2013a) argue that the substantial density contrast between iron and slag makes it more natural to assume separated flow and that the conclusions drawn from gas-liquid physical models may not be fully extrapolable to the real BF, as they are a rather extreme case for the taphole flow. Moreover, their results indicate that a separated flow of iron and slag is more likely to occur in the taphole of the studied BF. Thus, in this work, we assume that the flow in the taphole remains stratified.

In Figure 1, a sketch of a BF and its trough system is displayed. Also, in Figure 2, a 3D schematic diagram of an empty BFT is shown. Their dimensions and specific features widely depend on the casthouse arrangement and its particular needs. However, the general characteristics are common to most modern troughs. Separation is achieved as the denser hot metal flows underneath a skimmer block, which forces the lighter slag to accumulate on top of the hot metal and drain through a smaller trough that goes sideways, known as the slag runner. The iron dam at the end of the trough allows a residual amount of hot metal and slag to be kept inside the trough between tappings. Most modern BFTs retain a liquid pool or bath between BF tappings, with a depth usually over 20 cm. The fluid pool is not only useful to mitigate the kinetic energy of the falling jet stream, but also leads to better preservation of heat (Geyer & Halifa, 2014). In addition, the liquid bath and reduced slope result in lower hot metal velocities, which increases the residence time of slag and iron in the trough, allowing for better metal-slag separation. Furthermore, refractory wear and tear is greatly reduced at the jet impact point, and most of the refractory degradation is localized at the laterals of the trough (see e.g. Rüther & Lüngen, 1989).

Figure 1. Longitudinal sketch of a BF trough.

Figure 2. Schematic diagram of the main trough.

The model is defined using the data supplied by the company collaborating on the research project PID2019-105615RB-I00. It is a 3-taphole BF, with an estimated hot metal production of 6000 tonnes/day. In Figure 3, a drawing of the troughs operating in the BF is shown. Specifically, a view of the trough from the longitudinal cut highlighted in red in Figure 2 is displayed. Note that the total length of the BFT is about 18 m, whereas its width is 3 m. The part of the empty space in the trough considered as the computational domain in the subsequent sections is highlighted in blue in the aforementioned figure.

Figure 3. Longitudinal view of the BFT from the vertical plane highlighted in Figure 2. Dimensions are indicated in mm.

The physical property values considered for the fluid phases involved in the flow along the BFT (air, slag and hot metal) are collected in Table 1. In this table, we denote as ρ and μ the density and dynamic viscosity, respectively. The reported values for hot metal and slag correspond to plant data estimations at a standard operating temperature of 1500${}^{\circ }$C, supplied by the steelmaking company.

Table 1. Properties of the fluid phases. Slag and hot metal are assumed at a temperature of 1500${}^{\circ }$C.

	$\rho \left[\text{kg}/{\text{m}}^3\right]$	$\mu [\text{Pa}\cdot \text{s}$]
Hot metal	7015	7.15e−3
Slag	2600	0.50
Air	1.21	1.82e−5

3. Mathematical model

3.1. Computational domain

Provided that the scope of this work is to address the hydrodynamics of the fluid phases in the trough, only the empty space within the trough is included in the 3D computational domain Ω, following the drawing in Figure 3. The dimensions indicated there are those of the refractory linings at the time of commissioning a new trough. Therefore, the refractory lining is assumed to be in perfect condition to generate the computational domain. The domain is depicted in Figure 4(a). Note that the lateral and bottom boundaries of the domain coincide with the surface of the refractories, whereas the top boundary is a fictitious boundary, as any of the fluid phases may freely flow through it. In Figure 4(c), the relative position of the computational domain with respect to the actual BFT is shown.

Figure 4. (a) Computational domain (Ω), (b) jet trajectory schematic, and (c) position of the domain with respect to the solid refractory linings.

We consider that the taphole has an inclination angle ${\gamma }_{tap}={10}^{\circ }$ (see Figure 4(b)), which is the drilling angle in the BF corresponding to the trough modelled in this work. It has been observed numerically and experimentally that the hot metal and slag jet closely follows a parabolic profile, as it flows out at a high velocity value. According to the slag and hot metal production data, the velocity at the taphole opening is estimated to be around $v_{tap}=6.29\mathrm{m}/\mathrm{s}$. In addition, an initial taphole diameter of $D_{tap}=62\mathrm{m}\mathrm{m}$ is considered. However, the influence of other taphole diameters and stream velocities will be discussed in Section 5.4.

As it is not the focus of this paper, we do not solve the complete jet starting from the taphole. Instead, the computational domain only includes a portion of air above the liquids (slag and hot metal) flowing through the BFT. A suitable entry point of the jet on its fictitious upper boundary is computed, imposing an inlet boundary condition, as depicted in Figure 4(a). This allows us to reduce the computational cost and facilitate the mesh generation procedure (see Section 4.1).

The numerical results in Barral et al. (2021) show that simulating the complete jet trajectory yields a very similar impingement point compared to that expected from the computations assuming a parabolic profile and neglecting drag. Hence, the drag force plays a small role, and therefore the trajectory of the free jet after leaving the taphole can be reasonably estimated considering only gravitational acceleration, provided that the taphole exiting stream remains as a dense, non-splashy stream until reaching the BF trough. The procedure to determine the jet entry point in the domain, described below, is similar to that previously considered by other authors (e.g. Ge et al., 2020; Yao et al., 2021).

We set the coordinate axes and the origin, $\mathbf{O}$, as in Figure 4, and let $\mathbf{x}\left(t\right)=(x_1,x_2,x_3)\left(t\right)$ be the trajectory of the jet as a function of the time parameter, t. Hence, the known initial position of the jet, corresponding to the position of the taphole centre, is denoted as (1) $\begin{array}{r}\mathbf{x}\left(0\right)={\mathbf{x}}^0=(x_1^0,x_2^0,x_3^0).\end{array}$(1) Specifically, for the selected origin of coordinates, it holds that ${\mathbf{x}}^0=(1.5,2,0)\mathrm{m}$, provided that the total width of the trough is 3 m and that the taphole axis is in the vertical central plane of the trough. In addition, the initial velocity is $\mathbf{v}\left(0\right)={\mathbf{v}}^0=(v_1^0,v_2^0,v_3^0)$, with ${\mathbf{v}}^0={\mathbf{x}}^{\prime }\left(0\right)$. A schematic diagram of the jet trajectory is depicted in Figure 4(b). Since we assume that the trajectory of the jet exiting the taphole is only affected by gravity, it follows that (2) $\begin{array}{rl}& x_1\left(t\right)=x_1^0,\end{array}$(2) (3) $\begin{array}{rl}& x_2\left(t\right)=-9.81\frac{t^2}{2}+v_2^{0}t+x_2^0,\end{array}$(3) (4) $\begin{array}{rl}& x_3\left(t\right)=v_3^{0}t+x_3^0.\end{array}$(4) From (4(4) $\begin{array}{rl}& x_3\left(t\right)=v_3^{0}t+x_3^0.\end{array}$(4) ), it holds that $t=(x_3-x_3^0)/v_3^0$, and substituting this value in (3(3) $\begin{array}{rl}& x_2\left(t\right)=-9.81\frac{t^2}{2}+v_2^{0}t+x_2^0,\end{array}$(3) ) the height of the jet ($x_2$ coordinate) can be described as a function of the distance to the BF, $x_3$. The corresponding function is denoted as $h_{jet}$, which satisfies that (5) $\begin{array}{r}{h}_{jet}\left(x_3\right)=-{\displaystyle \frac{9.81}{2}}{\left({\displaystyle \frac{x_3-x_3^0}{v_3^0}}\right)}^2+v_2^0\left({\displaystyle \frac{x_3-x_3^0}{v_3^0}}\right)+x_2^0.\end{array}$(5) The trough has an inclination of approximately 0.5%. Specifically, with the choice of coordinate origin shown in Figure 4(b), the height of the bottom of the trench through which the hot fluids circulate can be parametrized as a function of the $x_3$ coordinate as follows: (6) $\begin{array}{r}{h}_B\left(x_3\right)=1.07-(x_3-2)\frac{7}{1435},\end{array}$(6) for $x_3\in [2,16.35]\mathrm{m}$ (see Figure 3).

Noting that ${\gamma }_{tap}={10}^{\circ }$ and for an initial velocity modulus equal to $|{\mathbf{v}}^0|=6.29\mathrm{m}/\mathrm{s}$, we find that ${\mathbf{v}}^0=(0,1.09,6.19)$ m/s. Then, using that ${\mathbf{x}}^0=(1.5,2,0)\mathrm{m}$, the solution of $h_{jet}\left(x_3\right)=h_B\left(x_3\right)$ yields the theoretical longitudinal coordinate at which the free jet impacts the bottom surface, denoted as $x_3^{imp}=3.48\mathrm{m}$. Moreover, the corresponding entry point is given by the height of the computational domain, equal to $1.7\mathrm{m}$.

3.2. Turbulence modelling

To model the turbulence in the flow, we use Menter's SST $K-\omega$ model, originally proposed in Menter (1994);Menter et al. (2003). The SST $K-\omega$ turbulence model was conceived as a modification of Wilcox's $K-\omega$ model (Wilcox, 2006). To approximate the eddy viscosity, in addition to the turbulent kinetic energy, K, the model introduces the specific dissipation rate, denoted as ω. It is effectively a hybrid model, which transforms the $K-\omega$ turbulence model into a $K-ϵ$ model away from the walls. This maintains the favorable characteristics of the $K-\omega$ model near the wall while avoiding some of its significant drawbacks. Namely, Wilcox's $K-\omega$ does not perform well under adverse pressure gradients and has been found to be sensitive to the free stream values of the turbulence variables (Versteeg & Malalasekera, 2007).

Most of the CFD studies in the literature concerning the flow in the BFT resort to the standard $K-ϵ$ model (Ge et al., 2020; Monteiro de Oliveira et al., 2020). However, $K-\omega$ models are reported to perform better in wall-bounded flows and can be readily used without wall functions. In particular, the SST model was used to obtain shear stress predictions for the planar hot metal jet in Barral et al. (2021). Nevertheless, to obtain fully accurate predictions for impinging jets, more detailed models should be considered (see e.g. Kaewbumrung & Plengsa-Ard, 2023), although approaches involving RANS turbulence models are popular in the literature (e.g. Mahdavi et al., 2021). In the particular setting considered in this study, the use of the computationally cheaper SST $K-\omega$ was favoured to balance the computational complexity and robustness against the expected accuracy, as in many other studies involving CFD simulations of various situations (e.g. Gao et al., 2016).

The conservation equations that are included to model the evolution of K and ω in the SST turbulence model are stated in Section 3.4. The various constants, source terms and other expressions that are involved in those equations and in the computation of the eddy viscosity are gathered for completeness in the Appendix.

3.3. Volume of fluid

The flow in the BFT has three immiscible fluid phases –air, slag and hot metal– which are separated by density stratification. To account for the multiphase flow, the volume of fluid (VOF) method is used (Hirt & Nichols, 1981). It is a front-capturing method widely used to solve fluid mechanics problems with several involved phases. The interfaces between the phases are not tracked explicitly by moving the mesh or markers. Instead, their position is described implicitly with suitable scalar functions, assuming a single fluid domain with a fixed boundary. This technique has been used in various studies available in the literature featuring the multiphase and turbulent flow in the BFT (e.g. Barral et al., 2021; Monteiro de Oliveira et al., 2020; Rezende et al., 2007).

The domain Ω is separated in open subdomains: ${\mathrm{\Omega }}_{air}\left(t\right)$, ${\mathrm{\Omega }}_{slag}\left(t\right)$ and ${\mathrm{\Omega }}_{hm}\left(t\right)$, one for each different phase (air, slag and hot metal). On each subdomain, we consider a constant density and viscosity, equal to that of the corresponding fluid (see Table 1). Therefore, these properties are defined as piecewise constant functions. For instance, the viscosity reads as follows: (7) $\begin{array}{r}\mu (\mathbf{x},t)=\{\begin{array}{ll}{\mu }_{air},& \mathbf{x}\in {\mathrm{\Omega }}_{air}\left(t\right),\\ {\mu }_{slag},& \mathbf{x}\in {\mathrm{\Omega }}_{slag}\left(t\right),\\ {\mu }_{hm},& \mathbf{x}\in {\mathrm{\Omega }}_{hm}\left(t\right),\end{array}\end{array}$(7) where ${\mu }_{air}$, ${\mu }_{slag}$ and ${\mu }_{hm}$ denote the viscosity of air, slag and hot metal, respectively. The density, $\rho (\mathbf{x},t)$, is defined analogously.

To account for the motion of the interfaces separating the fluid phases, the evolution of a scalar characteristic function for each phase, ξ, is tracked as time advances (see e.g. Tryggvason et al., 2011). Numerically, these functions are approximated by the so-called volume fractions, denoted as ${\alpha }_{slag}$ for the slag phase and ${\alpha }_{hm}$ as well as ${\alpha }_{air}$ for hot metal and air, respectively. For example, considering the slag phase, in a cell-centred finite volume discretization scheme, the corresponding volume fraction is defined as the average of ${\xi }_{slag}$ in each computational cell of the mesh, i.e. if $N_{cell}$ denotes the number of cells in the mesh, then, for $1\le i\le N_{cell}$, we set that (8) $\begin{array}{r}{\alpha }_{slag}^i\left(t\right)={\displaystyle \frac{1}{\text{Vol}\left(C^i\right)}}{\int }_{C^i}{\xi }_{slag}(\mathbf{x},t)\mathrm{d}V,\end{array}$(8) where $C^i$ is the ith cell in the computational grid of the domain and $\text{Vol}\left(C^i\right)$ denotes the volume of such cell.

Thus, a value of ${\alpha }_{slag}^i\left(t\right)=1$ represents that slag fills the cell $C^i$. Values of the volume fraction different from 0 and 1 describe the fraction of the cell volume that is occupied by the slag. Under these circumstances, a mixture of fluid phases is assumed to fill the cell, indicating the transition region where the interfaces are located.

With the VOF method, it is customary to write the model equations in terms of the volume fractions also in the continuous setting. Specifically, from the advection equation for the characteristic function, we write that (9) $\begin{array}{r}{\displaystyle \frac{\mathrm{\partial }{\alpha }_{slag}}{\mathrm{\partial }t}}+\text{div}\left({\alpha }_{slag}\mathbf{v}\right)=0.\end{array}$(9) The VOF method only solves one set of conservation equations for all the fluid phases. It considers mixture properties weighted appropriately using the corresponding volume fractions. For an arbitrary scalar property field ϕ it is considered that (10) $\begin{array}{r}\hat{\varphi }={\alpha }_{slag}{\varphi }_{slag}+{\alpha }_{hm}{\varphi }_{hm}+(1-{\alpha }_{slag}-{\alpha }_{hm}){\varphi }_{air},\end{array}$(10) where we recall that ${\varphi }_{slag}$, ${\varphi }_{hm}$ and ${\varphi }_{air}$ denote the scalar property values corresponding to slag, hot metal, and air, respectively. The volume fraction of air can be directly obtained from the values of the volume fractions of slag and hot metal. Therefore, for the three-phase flow in the BFT, the VOF method adds two new unknowns, ${\alpha }_{slag}$ and ${\alpha }_{hm}$, with the corresponding two advection equations.

Remark 3.1

The importance of surface tension effects can be determined by computing several dimensionless numbers to assess whether the surface tension forces are overcome by other effects in the flow (Tryggvason et al., 2011). Since the flow in the trough is turbulent, especially in the hot metal phase, the value of the Weber number is analyzed, which measures the ratio among the inertia and surface tension forces. It is defined as (11) $\begin{array}{r}We={\displaystyle \frac{\rho \mathcal{L}{\mathcal{V}}^2}{\sigma }},\end{array}$(11) where $\mathcal{L}$ is the characteristic length, $\mathcal{V}$ the characteristic velocity and σ is the surface tension.

Considering a value of $\sigma =1.25\mathrm{N}/\mathrm{m}$ for the interfacial tension between hot metal and air, since ${\rho }_{hm}=7015\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ (see Table 1), it is clear that the surface tension is much less relevant than the inertia forces, both if the length scale and velocity are selected at the taphole or for the characteristic velocity and length scale in the middle of the trough. Hence, in all the results that are presented in this paper, surface tension is neglected. Nonetheless, it should be noted that in the jet impingement area, the flow features are intricate and significant bubble entrainment occurs. In the length scale of these bubbles, surface tension effects are indeed relevant, but a detailed resolution of these features of the flow requires much finer meshes and the use of more complex turbulence models.

3.4. Model equations

In this section, we briefly state the equations composing the considered mathematical model. The unsteady RANS equations are coupled with conservation equations for the volume fractions of slag and hot metal. To address turbulence in the flow, the SST $K-\omega$ model is used. This leads to the following system of equations, solved in the computational domain Ω over the time interval $t\in (0,t_{end}]$, where $t_{end}>0$ is a suitable end time value.

Find $\mathbf{v},p,{\alpha }_{slag},{\alpha }_{hm},K$ and ω such that (12) $\begin{array}{rl}& \text{div}\left(\mathbf{v}\right)=0,\end{array}$(12) (13) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }{\alpha }_{slag}}{\mathrm{\partial }t}}+\text{div}\left({\alpha }_{slag}\mathbf{v}\right)=0,\end{array}$(13) (14) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }{\alpha }_{hm}}{\mathrm{\partial }t}}+\text{div}\left({\alpha }_{hm}\mathbf{v}\right)=0,\end{array}$(14) (15) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }\mathbf{v}\right)}{\mathrm{\partial }t}}+\text{div}(\hat{\rho }\mathbf{v}\otimes \mathbf{v})=-\text{grad}\left(p\right)\\ & +\text{div}\left(\right(\hat{\mu }+{\hat{\mu }}_T\left)2\mathbf{D}\right(\mathbf{v})-{\displaystyle \frac{2}{3}}\hat{\rho }K\mathbf{I})+\hat{\rho }\mathbf{g},\end{array}$(15) (16) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }K\right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }K\mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_K}})\text{grad}\right(K\left)\right]\\ & +G_K-Y_K,\end{array}$(16) (17) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }\omega \right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }\omega \mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_{\omega }}})\text{grad}\right(\omega \left)\right]\\ & +G_{\omega }-Y_{\omega }+D_{\omega }.\end{array}$(17) The various source terms appearing in Equations (16(16) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }K\right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }K\mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_K}})\text{grad}\right(K\left)\right]\\ & +G_K-Y_K,\end{array}$(16) ) and (17(17) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }\omega \right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }\omega \mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_{\omega }}})\text{grad}\right(\omega \left)\right]\\ & +G_{\omega }-Y_{\omega }+D_{\omega }.\end{array}$(17) ), as well as the expression used to compute the eddy viscosity, ${\hat{\mu }}_T$, are briefly stated in the Appendix. Moreover, recall that $\hat{\mu }$ and $\hat{\rho }$ denote mixture properties, computed as indicated in (10(10) $\begin{array}{r}\hat{\varphi }={\alpha }_{slag}{\varphi }_{slag}+{\alpha }_{hm}{\varphi }_{hm}+(1-{\alpha }_{slag}-{\alpha }_{hm}){\varphi }_{air},\end{array}$(10) ). In addition, since we are dealing with the RANS equations, the variables $\mathbf{v}$ and p in (12(12) $\begin{array}{rl}& \text{div}\left(\mathbf{v}\right)=0,\end{array}$(12) )–(17(17) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }\omega \right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }\omega \mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_{\omega }}})\text{grad}\right(\omega \left)\right]\\ & +G_{\omega }-Y_{\omega }+D_{\omega }.\end{array}$(17) ) should be understood as the mean velocity and pressure (Wilcox, 2006).

Remark 3.2

The mathematical model assumes an incompressible flow in the BFT. Although the velocities of hot metal and slag are below the threshold necessary to induce compressibility due to pressure, thermal effects may play a significant role, inducing some degree of buoyancy in the flow due to density variations.

However, experimental measurements at the hot metal pool at the end of the trough indicate that it remains at temperatures very close to 1500${}^{\circ }$C. This suggests that the bulk of the pool of fluids held in the trough during tappings is at a high and mostly constant temperature, except for the slag-air free surface, where significant heat loss occurs due to radiation emission and partial solidification may occur, especially near the skimmer.

The primary flow dynamics in the BF trough are governed by the jet impact in the retained liquid pool, resulting in high velocities that substantially exceed those that would be induced by natural convection in the jet impingement zone. Furthermore, since accurately accounting for thermal effects involves additional modelling difficulties, such as including the solid refractories and accounting for thermal radiation, it is common for similar studies in the literature to neglect them and assume an isothermal flow (see e.g. Monteiro de Oliveira et al., 2020).

3.5. Boundary conditions

The boundary of the computational domain, denoted as Γ, is decomposed as shown in Figure 5. We set the subsequent boundary conditions:

• On ${\mathrm{\Gamma }}_{wall}$, which corresponds to the refractory lining surface, the standard no-slip condition is imposed, i.e. $\mathbf{v}=\mathbf{0}$. Since the SST $K-\omega$ turbulence model is used, this condition is set by blending the linear law valid on the viscous sublayer with the logarithmic wall law as a function of the dimensionless wall distance. A blending function is also used to set a Dirichlet condition on ω (Menter & Esch, 2001). The normal gradients of K and the volume fractions are considered null.

• On ${\mathrm{\Gamma }}_{inlet}$, the fluid velocity corresponding to the intersection of the jet with the fictitious boundary at the top of Ω is set. To this end, we note that this boundary is the intersection between $\overline{\mathrm{\Omega }}$ and the plane $x_2=1.7\mathrm{m}$. Then, the procedure detailed in Section 3.1 enables the computation of the inlet velocity. For instance, for a tapping stream velocity equal to $v_{tap}=6.29\mathrm{m}/\mathrm{s}$, we find that the inlet velocity, denoted as ${\mathbf{v}}_{inlet}$, is such that (18) $\begin{array}{r}{\mathbf{v}}_{inlet}=(0,-2.66,6.19)\text{m/s}.\end{array}$(18) Moreover, the values of the turbulence model variables, K and ω, are also required. For a taphole diameter such that $D_{tap}\approx 0.06\mathrm{m}$, we find that the Reynolds number associated with this characteristic length, a characteristic velocity $\mathcal{V}=6.29\mathrm{m}/\mathrm{s}$, and the physical properties of hot metal is $R{e}_{D_{tap}}\approx 3.7e5$. Following the procedure detailed in Ansys, Inc. (2019), based on the taphole length scale, the mixing length, $l_{mix}$, and the turbulent intensity, $I_T$, are estimated using the following expressions: (19) $\begin{array}{rl}{I}_T& =0.16R{e}_{D_{tap}}^{-1/8}=0.032,\\ l_{mix}& =\frac{0.07D_{tap}}{C_{\mu }^{3/4}}=0.026\text{m},\end{array}$(19) being $C_{\mu }$ the constant defined in Table A1.

  Then, the inlet values of K and ω are estimated using the following relationships: (20) $\begin{array}{rl}K& =\frac{3}{2}(\mathcal{V}{I}_T{)}^2=0.061{\text{m}}^2/{\text{s}}^2,\\ \omega & =\frac{K^{1/2}}{C_{\mu }{l}_{mix}}=105.35{\text{s}}^{-1}.\end{array}$(20) Lastly, adequate ${\alpha }_{slag}$ and ${\alpha }_{hm}$ values are set assuming that the hot metal and slag jet reaches the computational domain without spreading. The values are adjusted depending on the analyzed case to match different mass flow rates through the taphole (see Section 4.3). When slag and hot metal are tapped simultaneously, it is assumed that the flow in the taphole is stratified (see the remark below).

• On ${\mathrm{\Gamma }}_{outlet}$, which corresponds to the fictitious boundaries of the computational domain, we specify the so-called $\mathit{p}\mathit{r}\mathit{e}\mathit{s}\mathit{s}\mathit{u}\mathit{r}\mathit{e}\mathit{-}\mathit{o}\mathit{u}\mathit{t}\mathit{l}\mathit{e}\mathit{t}$ boundary condition, imposing a pressure value equal to p = 0. The velocity field required to achieve mass conservation given this pressure value is computed. This can result in fluid entering the computational domain through this boundary. In this scenario, since the upper boundary is located high enough to guarantee that slag and hot metal do not reach it, we assume that inward flow is air, i.e. that ${\alpha }_{air}=1$. Furthermore, if inflow occurs, it is set that $K=0.0001{\text{m}}^2/{\text{s}}^2$ and $\omega =0.1{\text{s}}^{-1}$, which ensure low values of turbulence in the incoming air across the boundary and improve the convergence of the calculation. Otherwise, if the flow is outwards, the normal gradients of all variables (except the pressure) are assumed to be zero.

Figure 5. Boundaries of the computational domain Ω.

Remark 3.3

The refractory lining profile within the taphole is uneven, and substantial roughness can occur due to severe erosion and BF drilling practices. Consequently, the exiting stream becomes highly turbulent, occasionally causing entrained gas to result in significant spreading of the jet before reaching the liquid pool. Moreover, as discussed previously, the flow within the taphole may not be fully stratified. The precise conditions of the flow in the taphole are not fully understood. Although some degree of dispersed flow may be expected, especially if gas entrainment occurs, the stratified taphole flow pattern that is implicitly assumed in this study, with slag sitting on top of the hot metal when both are tapped simultaneously, aligns well with the available literature (see e.g. Shao & Saxén, 2011, 2013a).

In any case, it should be noted that the assumption that the jet enters the computational domain as a cohesive, density-stratified stream when the flow is a combination of slag and hot metal may not be very realistic. An accurate characterization of the jet flow would require emulating the taphole conditions and the flow within the BF, as well as the use of more detailed turbulence models, which are out of scope for our aims in this study due to the high computational costs and modelling uncertainties. If the 3D jet is simulated as exiting the taphole with a RANS model, the shearing with the surrounding air is not enough to induce significant spreading until the impact point is reached (Barral et al., 2021), and the jet closely follows the parabolic profile determined by (5(5) $\begin{array}{r}{h}_{jet}\left(x_3\right)=-{\displaystyle \frac{9.81}{2}}{\left({\displaystyle \frac{x_3-x_3^0}{v_3^0}}\right)}^2+v_2^0\left({\displaystyle \frac{x_3-x_3^0}{v_3^0}}\right)+x_2^0.\end{array}$(5) ).

3.6. Initial conditions

At the beginning of each BF tapping (except for the first one in the BFT campaign), there is a pool of fluids in the trough, which mitigates the impact of the jet on the refractory and also helps slag and hot metal separation. For the initial condition, we assume that the trough is filled with this pool of fluids and, therefore, we set ${\alpha }_{hm}$ and ${\alpha }_{slag}$ with non-zero values at $t=0\mathrm{s}$.

To do so, the positions of the interfaces need to be computed at the initial time. With this purpose, we follow the approach described in Luomala et al. (2001) and Wang et al. (2017). The trough is assumed to be filled with fluids from a previous tapping. Hence, the slag level in the trough should reach the level of the slag runner, denoted as $h_{slag}^0$ in Figure 6. The 0 superscript indicates that it refers to the height of the free surface at the initial time. Moreover, the hot metal level after the skimmer should be approximately at the height of the iron dam, denoted as $h_{hm}^{\ast ,0}$. However, the height of the interface separating slag and hot metal upstream of the skimmer, $h_{hm}^0$, is unknown. If equilibrium conditions hold after the end of each tapping, assuming that the slag and hot metal remain at the constant densities collected in Table 1, the pressure due to the column of slag and hot metal upstream of the skimmer should be equal to that resulting from the column of hot metal downstream (see Figure 6). Thus, it holds that: (21) $\begin{array}{r}{H}_{hm}^0{\rho }_{hm}+H_{slag}^0{\rho }_{slag}=H_{hm}^{\ast ,0}{\rho }_{hm},\end{array}$(21) where $H_{slag}^0$ is the depth of the slag in the trough, $H_{hm}^0$ denotes the height difference between the slag-hot metal interface and the bottom part of the skimmer, and $H_{hm}^{\ast ,0}$ is the height difference between the iron dam and the bottom of the skimmer at t = 0 s.

Figure 6. Schematic diagram of the height values used in the ${\alpha }_{hm}$ and ${\alpha }_{slag}$ initial condition computation.

For our particular trough configuration, it holds that $H_{hm}^{\ast ,0}=0.2\mathrm{m}$ and $h_{slag}^0=1.5\mathrm{m}$ at any time during a standard BF tapping (see the drawing in Figure 3). In addition, it holds that $H_{hm}^0+H_{slag}^0=h_{slag}^0-h_{skimmer}=0.4\text{m},$where $h_{slag}^0$ is the height of the slag trough, as discussed above, and $h_{skimmer}=1.1\mathrm{m}$ is the height of the bottom of the skimmer. Then, from (21(21) $\begin{array}{r}{H}_{hm}^0{\rho }_{hm}+H_{slag}^0{\rho }_{slag}=H_{hm}^{\ast ,0}{\rho }_{hm},\end{array}$(21) ) it follows that (22) $\begin{array}{r}{H}_{hm}^0={\displaystyle \frac{H_{hm}^{\ast ,0}{\rho }_{hm}-0.4{\rho }_{slag}}{{\rho }_{hm}-{\rho }_{slag}}}\approx 8.22\cdot {10}^{-2}\text{m}.\end{array}$(22) Hence, we find that the free surface separating air and slag is at a height $h_{hm}^0=h_{skimmer}+H_{hm}^0=1.182$ m. Using these heights, the values of ${\alpha }_{slag}$ and ${\alpha }_{hm}$ are set at t = 0 s. The height values reported in this section are measured along the $x_2$ coordinate axis, as shown in Figure 4. For clarity, subsequent sections normalize these height values by subtracting the $x_2$ coordinate corresponding to the lowest point of the trench, marked with a red asterisk in Figure 6, where $x_2=1\mathrm{m}$ (see Figure 3).

Finally, we summarize the remaining initial conditions. Since the fluids are initially assumed at rest, the mean velocity is $\mathbf{v}=\mathbf{0}\mathrm{m}/\mathrm{s}$. In addition, the turbulent kinetic energy is $K=0{\mathrm{m}}^2/{\mathrm{s}}^2$ and for the specific dissipation rate we consider a standard value of $\omega =1{\mathrm{s}}^{-1}$.

4. Numerical solution

4.1. Mesh of the domain

A full hexahedral mesh of the computational domain is constructed, allowing the overall number of cells to be reduced while improving the general accuracy of the computation. Hexahedral cells lead to a more accurate interface reconstruction, especially in this specific application, since it is known beforehand that the free surfaces in most of the computational domain have a normal vector whose direction is approximately that of ${\mathbf{e}}_2$.

The above-mentioned consideration does not apply to the jet itself and the region where it impinges the liquid pool, since the kinetic energy of the impact leads to significant air and slag entrainment and mixing with hot metal. Therefore, the interfaces become curved, and the flow is partially buoyancy-driven due to the rising entrained air bubbles, as has been observed experimentally with scale models (Begnis & Brandaleze, 2006; Monteiro de Oliveira et al., 2020). Hence, the grid is refined in those areas, where the most complicated features of the flow are generated, as shown in Figure 7 for the most refined of the constructed meshes, labelled as Mesh 4. Notice that the mesh is also finer near the bottom of the trough and at the expected positions of the slag-air and slag-hot metal interfaces downstream of the jet impingement zone.

Figure 7. Different views of Mesh 4 considered in this paper.

4.2. Discretization and numerical methodology

The software Ansys Fluent, Release 2019 R3, is used to solve the mathematical model. The SIMPLE segregated algorithm is used for the pressure and velocity coupling. The least squares cell-based scheme is utilized for gradient computation, whereas Second-Order Upwind schemes are used for the convective terms in all the conservation equations. Moreover, the PRESTO! scheme is considered to interpolate pressure values at the cell faces. It employs a staggered grid for the pressure values. Time discretization is done using the so-called First-Order implicit scheme. Lastly, to address the volume fraction advection equations, we use the Geo-Reconstruct scheme. It is a variant of the PLIC technique, which produces minimal numerical diffusion in the interface advection (Tryggvason et al., 2011). For the description of these numerical schemes and their implementation in the used software package, we refer to Ansys, Inc. (2019). Some additional information on part of the used discretization schemes can be found in Versteeg and Malalasekera (2007).

The numerical scheme for VOF advection greatly reduces the numerical diffusion of the interface, but it is an explicit scheme. Therefore, small time steps are required to comply with the stability constraints. To reduce the computational effort, the multiphase-specific adaptive time stepping technique is used (Ansys, Inc., 2019). It modifies the time step size based on a time scale calculated each time step by comparing a specified target global Courant number with that computed in the simulation. The global Courant number is defined using a combination of a flux-based and a velocity-based criterion, the latter being the most restrictive. It measures the ratio between the volume of the interface cells and the sum of the outgoing face fluxes. The starting time step size for all simulations is 1e−5 s, with 10 fixed-size time steps being always performed. Then, the step is modified according to the computed time scale, with a maximum change ratio of 1.05 and a minimum of 0.4. A target Global Courant Number of 0.5 is set for the computation of the time scale, and the maximum step size depends on the specific grid (see Table 2).

Table 2. Computational grids, with mesh quality metrics (skewness and orthogonal quality) and maximum specified time-step size used in the adaptive time-stepping method.

		Skewness		Orth. quality
	# cells	Max.	Mean	Min.	Mean	Max. $\mathrm{\Delta }t$ [s]
Mesh 1	200900	0.82	0.38	0.22	0.77	2.5e−4
Mesh 2	412682	0.82	0.38	0.22	0.76	2e−4
Mesh 3	812178	0.82	0.38	0.21	0.77	1.5e−4
Mesh 4	1670712	0.83	0.38	0.21	0.77	1e−4

To assess whether the numerical solution has converged at each time step, we control the locally-scaled residuals corresponding to all the flow variables (Ansys, Inc., 2019), assuming that the time step has converged if all the residuals fall below 1e−5.

4.3. Analyzed cases

Several cases involving different stages of a typical BF tapping are simulated by varying the inlet boundary condition and the initial conditions of the mathematical model. Their main features are highlighted in Table 3. As discussed in Section 2, during the initial portion of the tapping, only hot metal flows out from the taphole. Usually, this stage lasts for around 15–30 minutes. Then, once slag reaches the taphole, both slag and hot metal are drained simultaneously. While the initial tapped amounts are similar (Shao & Saxén, 2012), the slag proportion in the taphole stream may change during the tapping (see e.g. Henriksson et al., 2016). Thus, considering a constant taphole diameter $D_{tap}=62\mathrm{m}\mathrm{m}$ and stream velocity $v_{tap}=6.29\mathrm{m}/\mathrm{s}$, we solve the subsequent cases:

• Run A. The volume fractions of slag and hot metal are set as ${\alpha }_{hm}=1$ and ${\alpha }_{slag}=0$ on ${\mathrm{\Gamma }}_{inlet}$.

• Runs B-D. The volume fractions on ${\mathrm{\Gamma }}_{inlet}$ are adjusted to match the corresponding values in Table 3. The initial condition is set as the state of the computational domain Ω computed in Run A at t = 150 s. The phases are assumed to be stratified in the taphole, with no mixing and separated by a flat interface (see Remark 3.3). Therefore, slag and hot metal are assumed to reach the computational domain maintaining this stratification, with slag on top of hot metal.

• Run A*. The first tapping of the BFT campaign is addressed. In this situation, the trough is empty and there are no fluids to mitigate the jet impact on the refractory. The same boundary conditions as in Run A are considered. Only the initial condition is changed, setting that ${\alpha }_{air}=1$ in Ω instead of that discussed in Section 3.6.

Table 3. Taphole diameter ($D_{tap}$), stream velocity ($v_{tap}$), total volumetric flow rate (Q), volumetric fraction and mass flow rate ($\dot{m}$) in the simulated cases.

				% Volume		${\dot{m}}_{tap}\left[\mathrm{k}\mathrm{g}/\mathrm{s}\right]$
	$D_{tap}\left[\mathrm{m}\mathrm{m}\right]$	$v_{tap}\left[\mathrm{m}/\mathrm{s}\right]$	$Q[{\mathrm{m}}^3/min]$	Slag	Hot metal	Slag	Hot metal
Runs A, A*	62	6.29	1.14	0	100	0	33.2
Run B	62	6.29	1.14	30	70	14.8	93.3
Run C	62	6.29	1.14	50	50	24.7	66.6
Run D	62	6.29	1.14	70	30	34.6	40.0
Run C1	62	6.29	1.14	50	50	24.7	66.6
Run C2	70	6.29	1.45	50	50	31.5	84.9
Run C3	80	6.29	1.90	50	50	41.1	110.9
Run C4	80	4.82	1.45	50	50	31.5	85.0
Run C5	80	5.55	1.67	50	50	36.3	97.8

For each of these cases, we solve a time interval $(0,t_{end}]$, being $t_{end}=200\mathrm{s}$, except for Run A*, for which the end time is set as $t_{end}^{\ast }=100\mathrm{s}$, since we are mainly interested on the first seconds, when there are no fluids.

During normal tapping operation of the BF, the taphole outflow conditions might vary significantly, as the taphole diameter may increase due to severe erosion. Also, the liquid velocity in the taphole, driven by the furnace pressure, might change. Wear is only substantial after the simultaneous tapping of slag begins, with typical wear rates being in the range of 0.11 mm/min (Cameron et al., 2019). Typically, total drainage rates from the BF rise after the slag arrival as the taphole wear increases, up to the time at which the tapping is ended (Henriksson et al., 2016).

To investigate the effect of changes in the taphole diameter and velocity on the flow in the BFT, a constant slag fraction in the taphole flow of $50\mathrm{\%}$ is considered, starting from an initial condition given by the solution of Run C at $t_{end}=200\mathrm{s}$. The following cases are solved:

• Run C1. The same parameters as in Run C are used, solving an additional time interval of 50 s of duration.

• Runs C2, C3. The taphole diameter is assumed to increase due to erosion of the taphole, modelled by enlarging the size of the ${\mathrm{\Gamma }}_{inlet}$ boundary accordingly. As the velocity is kept constant, the jet trajectory is preserved, albeit with significantly higher slag and hot metal mass flow rates.

• Runs C4, C5. For a constant taphole diameter $D_{tap}=80\mathrm{m}\mathrm{m}$, the taphole stream velocity is set to $v_{tap}=4.82\mathrm{m}/\mathrm{s}$ and $v_{tap}=5.55\mathrm{m}/\mathrm{s}$, respectively, to analyze possible variations in the taphole outflow velocity.

For Runs C1–C5, a shorter time interval $(0,t_{end}]$ is solved, considering $t_{end}=50$ s. Note that the conditions of Run C1 are the same as Run C, and its purpose is to serve as a baseline for comparison of the remaining cases. The solution of Runs C2–C5 involves rebuilding the computational mesh described in Section 4.1, as both the diameter and the jet trajectory are being slightly changed. Hence, the initial condition is set by interpolating the final computed state of Run C to the mesh. Furthermore, the boundary condition on ${\mathrm{\Gamma }}_{inlet}$ is updated with the corresponding velocity to the point in the trajectory (see Sections 3.1 and 3.5).

5. Results and discussion

In this section, we discuss the main results that have been obtained after the solution of the mathematical model. First, in Section 5.1, a grid validation study is conducted, considering the initial and boundary conditions corresponding to Run A on four different computational grids. The aim is to assess the dependency of the computed result on the discretization parameter. Then, in Section 5.2, we discuss the numerical findings after solving Runs A-D with the conditions described above. Section 5.3 is devoted to analyze the computed solution for Run A*, where the trough is assumed to be dry as the initial condition. Lastly, in Section 5.4, the results from solving Runs C1–C5 are analyzed.

For the post-processing of the presented numerical results, the open-source software Paraview and the Matplotlib Python library (Hunter, 2007) have been used, as well as the proprietary computing environment MATLAB. The simulations were run in parallel, using 48 cores in a PowerEdge R840 machine, part of the Mariscal computer cluster. The machine is equipped with 4 Intel Xeon Gold 6126 processors and 384 GB of RAM. For each 10 s of time in the simulations, approximately 4 days of computation were required.

5.1. Grid validation

To address the grid dependency of the computed approximations to the solution, we study the differences between the results for the grids collected in Table 2, where the cell count is provided, as well as some mesh quality metrics. We report the mean skewness and orthogonal quality, as well as the maximum skewness and minimum orthogonal quality. Note that with the definition used by the software package, the best possible skewness and orthogonal quality are 0 and 1, respectively. Also, observe in Figure 7 that most of the hexahedral cells in the mesh are very regular, and the only significantly distorted cells are those in the vicinity of the jet trajectory, where some degree of distortion was required to follow the expected jet path with the mesh. In addition, the maximum time step size that we imposed on the adaptive step size controller is shown in Table 2. The conditions for this study are those described in Section 4.3 for Run A, i.e. only hot metal is assumed to be draining from the BF, falling in the trough containing a pool of liquid slag and hot metal. according to the initial condition described in Section 3.6. For this study, the end time is set as $t_{end}^{\ast }=100\mathrm{s}$.

In Figure 8, we compare the mass outflow rate of slag and hot metal for each of the grids during the time interval $(0,100]$ s on ${\mathrm{\Gamma }}_{outlet}$. Note that all meshes predict a swift stabilization of the flow rates, which fluctuate around mean values of $-3\mathrm{k}\mathrm{g}/\mathrm{s}$ and $-120\mathrm{k}\mathrm{g}/\mathrm{s}$ for slag and hot metal, respectively. At the beginning of the simulation, there is no outflow of slag or hot metal. Then, as the jet impacts the pool of fluids within the trough, a small wave is formed, which travels from the impingement point along the trough. It reaches ${\mathrm{\Gamma }}_{outlet}$ in the slag runner at around $t=7\mathrm{s}$ and overflows the iron dam at roughly $t=8.2\mathrm{s}$, leading to hot metal exiting the computational domain through the portion of ${\mathrm{\Gamma }}_{outlet}$ at the end of the trough, as evidenced by the slag and hot metal mass flow rates. The results depicted both in Figure 8(a) for slag and in Figure 8(b) for hot metal show good qualitative agreement among all considered meshes. However, the level of agreement deteriorates as time advances, particularly between Mesh 1 and the finer ones. The results for Mesh 3 and Mesh 4 are very similar.

Figure 8. Outflow rates for the meshes considered in the grid sensitivity study. (a) Slag (b) Hot metal.

In Figure 9, the hot metal-air and slag-hot metal interfaces are depicted at t = 25 s for the four computational grids. The position of these interfaces is assumed to correspond to the contour ${\alpha }_{hm}=0.5$. The interfaces are coloured by pressure values. The legend values are adjusted to the maximum and minimum pressure computed on the visible part of the contour corresponding to Mesh 4. It is observed that Meshes 2, 3, and 4 predict comparable results. Nonetheless, the predicted flow patterns become more intricate as the mesh is refined, with significant air and slag entrainment in the impact zone and complex vortices being formed. These lead to bubble formation, which can only be simulated with the finer meshes. An accurate approximation of these flow characteristics requires including surface tension effects, which have been neglected (see Remark 3.1) and, possibly, the use of more detailed turbulence models, such as LES models to reliably simulate the larger turbulent eddies. In all grids, it is observed that the hot metal in the jet impingement zone is displaced away, and the interface height is significantly lower at the flow stagnation region formed at the beginning of the trough when compared to that computed downstream of the impingement zone. This observation is qualitatively similar to the findings in the scale model studied in Begnis and Brandaleze (2006). On the other hand, the interfaces between slag-iron and slag-air remain clear and horizontal away from the impact zone, as observed in other numerical studies, such as Kou et al. (2019).

Figure 9. Contour ${\alpha }_{hm}=0.5$ at $t=25\mathrm{s}$ for the meshes considered in the grid sensitivity study, coloured by pressure values. The contour is interpreted as the slag-hot metal interface upstream of the skimmer and the hot-metal air free surface downstream of it. (a) Mesh 1 (b) Mesh 2 (c) Mesh 3 (d) Mesh 4.

As mentioned previously, it is of industrial interest to assess the wall shear stress and pressure profiles generated by the hot metal flow in the trough. In particular, the shear stress is closely related to the erosive phenomena, which play a significant role in refractory degradation. In Figure 10, the shear stress and pressure computed with each grid is plotted along the line given by the parametrization $(1.5,h_B(x_3),x_3)$, $x_3\in [3,5]$, with $h_B$ being the bottom of the computational domain, defined in (6(6) $\begin{array}{r}{h}_B\left(x_3\right)=1.07-(x_3-2)\frac{7}{1435},\end{array}$(6) ). It divides the part of ${\mathrm{\Gamma }}_{wall}$ highlighted in blue and labelled as jet impingement area in Figure 5 in two halves.

Figure 10. Wall shear and pressure on the line crossing the jet impingement area. (a) Wall shear, $t=100\mathrm{s}$ (b) Wall shear, averaged over $(75,100]\mathrm{s}$ (c) Pressure, $t=100\mathrm{s}$ (d) Wall shear close to the impingement point, averaged over $(75,100]\mathrm{s}$.

Remark 5.1

The large differences in the instantaneous wall shear values, depicted in Figure 10(a), can be attributed to the observed changes in the volume fractions for the finest meshes, which predict some degree of bubble formation and entrainment. Due to the omission of the surface tension and since the mesh resolution is not fine enough, the size of the bubbles predicted by the model is not fully accurate. As the wall shear stress directly depends on the viscosity of the fluid, fluctuations in the volume fraction due to entrained air reaching the bottom wall have a significant effect on the computed values.

Since we are only concerned about the main characteristics of the flow, we consider averaged variables in the following sense. Let $\varphi (\mathbf{x},t)$ be the computed solution of an arbitrary flow variable. Then, we denote the averaged value of $\varphi (\mathbf{x},t)$ as ${\mathcal{M}}_{\mathcal{T}}\left(\varphi \right)(\mathbf{x},t)$, defined as (23) $\begin{array}{r}{\mathcal{M}}_{\mathcal{T}}\left(\varphi \right)(\mathbf{x},t)=\frac{1}{\mathcal{T}}{\int }_{t-\mathcal{T}}^t\varphi (\mathbf{x},s)\mathrm{d}s,\end{array}$(23) where $\mathcal{T}$ is a suitable time scale.

The interest in these new averaged values is that they allow to filter the large fluctuations observed in the volume fractions. This procedure is applied regardless of whether the original flow value corresponds to a mean value in the RANS equations, such as the velocity or the pressure, as they are still affected by the volume fraction changes. To avoid confusion, we shall refer to these variables as ‘averaged’. The time interval during which the integration (23(23) $\begin{array}{r}{\mathcal{M}}_{\mathcal{T}}\left(\varphi \right)(\mathbf{x},t)=\frac{1}{\mathcal{T}}{\int }_{t-\mathcal{T}}^t\varphi (\mathbf{x},s)\mathrm{d}s,\end{array}$(23) ) is performed is $\mathcal{T}=25\mathrm{s}$ in all cases.

The averaged wall shear stress at $t=100\mathrm{s}$ is depicted in Figure 10(b). The large shear stress oscillations shown previously in Figure 10(a), attributed to entrained air reaching ${\mathrm{\Gamma }}_{wall}$, have been filtered by the averaging. In Figure 10(d), the shear stress is displayed only near the jet impact point. The plot resembles those experimentally measured for submerged jets in various conditions. Specifically, for a normally impinging circular jet, two symmetrical shear stress maximums are observed around a minimum, which corresponds to the impact point (Beltaos & Rajaratnam, 1974; Phares et al., 2000). In the plane hot metal jet numerical experiment considered in Barral et al. (2021), the maximums were not symmetrical due to the oblique impact of the jet on the refractory wall, which is precisely the case in the reported results herein. Notice that, while all meshes seem to agree in the predicted impingement point, Mesh 1 and 2 yield substantially different maximum averaged stresses. On the other hand, pressure values are similar for all meshes, as displayed in Figure 10(c).

Figure 11. Wall shear stress (${\tau }_w$) on the wall at the jet impingement zone, for the different considered grids. Top: at $t=100\mathrm{s}$. Bottom: averaged over $(75,100]\mathrm{s}$.

In Figure 11 the contours of wall shear stress are shown at $t=100\mathrm{s}$ for the different grids. The contours are displayed only on the part of ${\mathrm{\Gamma }}_{wall}$ that corresponds to the jet impingement area (see Figure 5). To better compare the results among the computational grids, the contours appearing in the legend for the top row have been adjusted to a maximum value of 1000 Pa. In the finer meshes, the shear stress values can be locally far higher (up to 1640 Pa in Mesh 4). In agreement with the prior discussion, the observed large differences are a consequence of the finer meshes being able to capture more effects regarding air and slag droplet entrainment. Instead, if the shear stress is averaged over the time interval (75,100] s, a significantly higher degree of similarity is exhibited. As discussed above for the plot over the line passing through the impact point, the averaging filters the unsteady shear stress profiles and yields very similar results for Mesh 3 and Mesh 4.

To compare the velocity predictions among the considered grids, in Figure 12, the contours of the velocity magnitude are displayed on a cross-section of the trough corresponding to $x_3=7.35\mathrm{m}$ and determined by the normal vector $\mathbf{n}={\mathbf{e}}_3$. Only the part of the cross-section occupied by slag and hot metal is shown. The contours at the top of the figure are instantaneous ones at $t=100\mathrm{s}$, whereas those at the bottom are averaged values over $(75,100]\mathrm{s}$. In both cases, the range of the contours appearing in the legend is adjusted to the maximum and minimum value of velocity magnitude found for Mesh 4. There are substantial differences between all meshes, including the more refined ones, although the qualitative characteristics of the velocity are similar between Mesh 3 and 4. Also, as this section is located 4 m away from the jet impingement point, we observe that the computed instantaneous values show much less significant fluctuations.

Figure 12. Plane cross-section of the computational domain ($x_3=7.35\mathrm{m}$), only considering the part occupied by slag and hot metal. Top: velocity magnitude [m/s] contours at $t=100\mathrm{s}$. Bottom: velocity magnitude [m/s], averaged over $(75,100]\mathrm{s}$.

The presented results evidence a non-negligible degree of mesh dependency, even among the finer meshes (Mesh 3 and 4), especially if the instantaneous shear stress results or the velocities near the impingement point are compared. Far away from the impingement point, the velocity results are similar. If the computed fields are averaged in time, the fluctuations due to entrained bubbles are eliminated and small differences are observed for all the analyzed flow variables between Mesh 3 and Mesh 4, including the wall shear in the jet impingement area. Therefore, we conclude that the latter mesh is sufficiently fine to guarantee a reasonable degree of accuracy for the prediction of the general features of the flow. In the subsequent results, we solve the mathematical model for all the cases gathered in Table 3 using Mesh 4 (Runs A-D, C1) or a mesh of similar size and features (Runs C2-C5).

5.2. Effects of taphole stream slag-hot metal ratio

In this section, the simulations concerning variable volumetric fraction of slag and hot metal in the taphole flow are discussed, i.e. Runs A-D (see Section 4.3). The main characteristics of the flow are studied, as well as the evolution of the free surfaces separating the fluid phases, and the shear stress due to the impinging jet for the taphole conditions described in Table 3.

Concerning the evolution of the positions of the air-slag and slag-hot metal free surfaces, they are computed as follows. The contour ${\alpha }_{slag}=0.5$ indicates the position of both free surfaces. Since we are interested in assessing the position of the part of these interfaces that is not affected by the strong circulations formed in the jet impingement point, after computing the coordinates of the points defining the contour, we only consider those points that satisfy $x_3>5\mathrm{m}$. Then, we assume that a point in such contour corresponds to the slag-hot metal interface if ${\alpha }_{hm}>{\alpha }_{air}$ for the corresponding interpolation of the volume fraction values. If the previous condition is not satisfied, the point is assumed to belong to the air-slag interface. Separating the points in this manner, we repeat the procedure for several time values on Runs A-D constructing the boxplots shown in Figure 13 by computing the corresponding statistics involving the $x_2$ coordinates of all points in the contour. Atypical values have not been included in the boxplots for clarity purposes. The free surface heights reported in the plots correspond to normalized values obtained by adjusting the $x_2$ coordinate value of the contour by subtracting the $x_2$ value corresponding to the lowest point in the computational domain (marked with a red asterisk in Figure 6). As discussed in Section 3.6, the $x_2$ coordinate for this point is $x_2=1\mathrm{m}$.

Figure 13. Evolution of the interface height for Runs A-D. (a) Slag-air, Run A (b) Slag-hot metal, Run A (c) Slag-air, Run B (d) Slag-air, Run C (e) Slag-air, Run D (f) Slag-hot metal, Run B (g) Slag-hot metal, Run C (f) Slag-hot metal, Run D.

We observe that in around only 50 s, both the slag and hot metal interfaces reach a quasi-steady position in Run A, with very low dispersion among all points composing the corresponding contours. During these first 50 seconds, the interface height rises, especially the one separating slag and hot metal. Notice that the maximum initial thickness of hot metal is approximately 18 cm and, at the end of Run A, it has grown to almost 26 cm, which is a relative increase of around 40%. In Runs B-D, the arrival of slag in the trough results in an increase in the depth of the slag pool and a decrease in the depth of the hot metal layer. Moreover, we find that the dispersion in the computed heights of both interfaces is substantially larger, indicating that they are less plane. This is possibly due to slag also flowing downstream of the BFT, as in Run A the slag is mostly static in the whole computational domain. In Table 4, the median final heights of the interfaces are reported for each case. The initial heights, at t = 0 s in Run A, are $h_{hm}=0.182$ m and $h_{slag}=0.5$ m, considering the normalization discussed previously. At the end of Run D, the median height of the slag-hot metal interface is approximately equal to the initial value in Run A, while the slag level is 8 cm higher. Thus, since the slag pool depth is approximately equal to the difference between $h_{slag}$ and $h_{hm}$, a relative increase roughly below 30% is predicted from the initial condition to the end of Run D.

Table 4. Final median height of the slag-air ($h_{slag}$) and hot metal-slag ($h_{hm}$) interfaces, for Runs A-D.

	Run A	Run B	Run C	Run D
$h_{slag}$ [m]	0.535	0.554	0.568	0.584
$h_{hm}$ [m]	0.257	0.226	0.201	0.174

In Figure 14, the ${\alpha }_{slag}=0.5$ contour is displayed at t = 150 s for Runs A, C, and D. The global position of the slag-hot metal and slag-air interfaces is shown. While the qualitative features observed in the figure are similar for all cases, the flow generated by the jet impact is apparently less complex when slag is being tapped. Also, note that the jet impact produces significant air entrainment in the pool. The air bubbles predicted by the simulations resurface at a distance around 1-2 m away from the impact point, similar to that observed in He et al. (2002). Furthermore, for all the taphole slag-iron ratios, the jet reaches the bottom wall maintaining a coherent structure, without substantial spreading.

Figure 14. Contour ${\alpha }_{slag}=0.5$ at t = 150 s (Runs A, C, D) and ${\alpha }_{hm}=0.5$ at t = 100 s (Run A*), coloured by pressure values. (a) Run A (b) Run C (c) Run D (d) Run A*.

In Figure 15, the median height as a function of the longitudinal coordinate $x_3$ of all the interfaces is displayed for different time values for Runs A-D. To obtain the height values, the points in the contour ${\alpha }_{slag}=0.5$ are examined, identifying whether the point corresponds to the iron-air free surface or iron-slag interface as discussed above. Also, those points in the contour ${\alpha }_{hm}=0.5$ after the start of the skimmer ($x_3>14.85\mathrm{m}$) are represented, as the corresponding $x_2$ values yield the height of the iron-air interface downstream of the skimmer (see Figure 3). Then, the selected points in the contours are separated into subintervals according to their $x_3$ coordinate values. The median value in each subinterval is plotted in the figure. The results in Figure 15(a), which shows the interface location evolution for Run A, are consistent with the previously discussed findings. From the initial condition, note how the large kinetic energy unleashed by the jet impact first empties most of the hot metal upstream of the impingement point, as the slag-hot metal interface, shown at the left of the plot, is substantially higher at t = 1 s than in t = 50 s. The remaining interfaces show a very stable position, especially far away from the jet impingement area, with an almost completely flat profile for $x_3>5$ m. Also, we remark that the large inclination of the iron-air free surface at the right of the plot is due to the slope of the iron runner, which is 3.5%. In Runs B, C, and D, depicted in Figure 15(b–d), the position of the hot metal-air interface, downstream of the skimmer, is much more stable in comparison. This is a consequence of the initial condition that was selected for these simulations, as the initial values for all flow variables were those computed at t = 150 s in Run A. These numerical results are in qualitative agreement with the approximate positions of the interfaces separating oil, water, and air in the physical scale model analyzed in Begnis and Brandaleze (2006) (which emulate slag, hot metal, and air in the real setting, respectively), where flat profiles in the quiet region away from the impingement zone were also reported.

Figure 15. Median height of the interfaces at various time values for Runs A-D. (a) Run A (b) Run B (c) Run C (d) Run D.

The position of the interfaces is compared at the end time of the simulation for each run, i.e. at time t = 200 s, in Figure 16. The results show again that during a standard tapping, while the variations in height of the slag-air free surface are rather small, the hot metal-slag free surface suffers non-negligible changes during the process. Specifically, we observe that the heights computed for most of the hot metal-slag interface are those already reported at the end time in the boxplots in Figure 13. The changes in the slag and hot metal levels reported in this work have not been observed in the analyzed literature. Measuring these variables in the industrial operation setting is challenging, given that the hot metal-slag interface remains concealed beneath a thick slag layer within the trough. These factors can significantly affect BFT operation. Towards the end of the tapping process, when the slag ratio in the taphole may increase, there is a heightened risk of the slag descending below the skimmer and contaminating the hot metal. This poses a safety concern (Kumar et al., 2010). Moreover, if the volume fraction of slag in the taphole exceeds 70% (Run D), we anticipate the interface position descending even further. Additionally, the cyclic alternation during each tapping between slag and hot metal within the same wall region can contribute to an escalation in refractory wear rates.

Figure 16. Median height of the interfaces for Runs A-D, at $t=200\mathrm{s}$.

Concerning the flow characteristics in the trough, recirculation patterns and secondary eddies that are strongly dependent on the taphole conditions are observed. In Figure 17, the velocity vectors on the plane $x_1=1.5\mathrm{m}$ are shown for time $t=150\mathrm{s}$ in the zone where the iron and slag jet impinges the wall, for Runs A-D. The vectors are only displayed on those points where ${\alpha }_{hm}>0$ or ${\alpha }_{slag}>0$. To better highlight the differences between the simulated situations, the maximum velocity magnitude contour that is displayed in the legend does not correspond to the actual maximum, which is around 7 m/s for all runs. The contour ${\alpha }_{slag}=0.5$ is also depicted as a white curve. All the displayed fields correspond to averages over $(125,150]\mathrm{s}$. When only hot metal is being tapped (Run A), two slag recirculations in the stagnation region upstream of the impingement point are generated, which merge in a single eddy during simultaneous tapping of slag and hot metal at the BFT (Runs B-D). In all cases, a backward-advancing current is observed downstream of the jet impact zone. This current is more pronounced when only hot metal is being drained from the BF and becomes weaker as the slag-iron ratio in the taphole increases. This phenomenon has been observed in several studies (e.g. Luomala et al., 2001). In Run A, we also observe a backward circulation of iron near the slag. Notably, while the primary flow characteristics remain consistent between Runs B-D, the arrival of slag from the BF at the trough induces significant differences between these cases and Run A.

Figure 17. Velocity vectors on the plane $x_1=1.5\mathrm{m}$, with white contour lines corresponding to ${\alpha }_{slag}=0.5$ and velocity magnitude contours. All the results correspond to averages over $(125,150]\mathrm{s}$. (a) Run A (b) Run B (c) Run C (d) Run D

Similar remarks apply to Figure 18, where the same contours and vectors are displayed at t = 150 s but in the final part of the BFT, corresponding to the part where the slag runner goes sideways, the skimmer and the iron dam. Here, the contour ${\alpha }_{hm}=0.5$ is shown in addition to the contour ${\alpha }_{slag}=0.5$ as a white curve to also depict the iron-air free surface downstream of the skimmer. Immediately preceding the iron dam, we observe the formation of two hot metal circulations around the primary current, akin to those identified in the physical scale model investigated in Begnis and Brandaleze (2006). The flow underneath the skimmer is substantially stronger in Run A, as more hot metal is being drained from the BF. Before the skimmer, it is observed that the hot metal flow is plug-like, without any complicated recirculating patterns. In contrast, a large circulation of slag forms as it reaches the skimmer wall at the end of the BFT and, in Run A, slag flows backward in the free surface with the air.

Figure 18. Skimmer and slag runner zone: velocity vectors on the plane $x_1=1.5\mathrm{m}$, with contour lines corresponding to ${\alpha }_{slag}=0.5$ and ${\alpha }_{hm}=0.5$, coloured by velocity magnitude contours. The results correspond to averages over $(125,150]\mathrm{s}$. (a) Run A (b) Run B (c) Run C (d) Run D.

To further analyze the flow features, in Figure 19, velocity vectors are depicted on the plane $x_2=1.5\mathrm{m}$, with velocity magnitude contours. As in the previously discussed cases, we consider a maximum velocity contour of 1.5 m/s, and the results correspond to the averaged variables over the interval $(125,150]\mathrm{s}$. Since the displayed results are close to the slag-air interface, except near the jet entry region, the predicted velocities are low, especially in Run A. Slag is flowing backward along the sidewalls, and in Run D, the flow seems slightly more developed, with symmetrical circulations both downstream and upstream of the jet entry region.

Figure 19. Velocity vectors on the plane $x_2=1.5\mathrm{m}$, with velocity magnitude contours. All the results correspond to averages over $(125,150]\mathrm{s}$. (a) Run A (b) Run C (c) Run D.

The computed velocities clearly show that the complex flow in the trough is heavily influenced by the proportions of slag and hot metal in the taphole, despite the assumption of total constant volumetric flow rates in Runs A-D. The simultaneous tapping of slag and hot metal leads to different flow patterns compared to the initial stage of the tapping, when only hot metal is being drained. In all examined scenarios, intricate flow patterns characterized by large recirculations and return currents of both slag and hot metal are observed.

In Figure 20, the averaged $v_3$ contours are shown at times t = 150 s (top row) and t = 200 s (bottom row) on the cross-section $x_3=11\mathrm{m}$. Little time dependency is observed, and similar velocity profiles are found in slag and hot metal over time. In Run A, in accordance with the previously discussed results, only hot metal flows downstream of the trough at a significant speed, and two backward-advancing slag currents form in the walls. In Runs B-D, the slag flow develops, although small backward slag currents still appear in Run B.

Figure 20. Longitudinal component of the velocity ($v_3$) [m/s] on plane cross-section at $x_3=11\mathrm{m}$, only considering the part occupied by slag and hot metal. Top: averaged over $(125,150]\mathrm{s}$. Bottom: averaged over $(175,200]\mathrm{s}$.

Figure 21 displays the contours of averaged shear stress and pressure at $t=200\mathrm{s}$, illustrating the influence of the taphole slag fraction. As the fraction of slag in the taphole stream increases, the maximum shear stress values decrease, and the area with high stress values becomes more localized. While the pressure contours exhibit similar qualitative trends, the maximum pressure values, particularly at the impingement point, decrease significantly as the taphole slag-iron ratio increases. Table 5 summarizes the actual maximum values for each case shown in the figure. Note that, for Run D (70% of slag in the taphole), the shear stress is 31% lower than considering a hot metal jet (Run A). In Figure 22, the averaged shear stress and pressure are compared at times $t=100\mathrm{s}$ and $t=200\mathrm{s}$ on the line dividing the jet impingement area. We find that both the shear stress and pressure remain fairly stable during the simulation, with small variations over time in all cases. In Run D, the contour shape showing a clear local minimum of the shear stress in the impingement point vanishes, whereas in Run C, it is only barely visible.

Figure 21. Averaged wall shear stress (top) and pressure (bottom) on the wall at the jet impingement zone at t = 200 s for Runs A-D.

Figure 22. Averaged wall shear (left) and pressure (right) on a line crossing the jet impingement area, for Runs A-D.

Table 5. Maximum averaged wall shear stress and pressure for Runs A-D.

	Run A	Run B	Run C	Run D
Max. shear stress [Pa]	948.7	931.8	708.0	658.4
Max. pressure [Pa]	$6.64\mathrm{e}4$	$5.15\mathrm{e}4$	$3.55\mathrm{e}4$	$2.62\mathrm{e}4$

The notably lower shear stress and pressure values may be attributed to the high viscosity of slag, which prevents the jet from reaching the surface at higher speeds, causing deceleration due to the substantial shear forces between the slag in the jet and that in the fluid pool. The velocity contours depicted in Figure 17 also illustrate this effect, where a significantly wider spread of the high velocity zone is observed when the jet is composed of both slag and hot metal.

5.3. Dry BF trough

In this section, the results obtained after solving Run A* are discussed. As described in Section 4.3, the aim is to investigate the shear stress and pressure on the jet impingement region during the first tapping of a dry BFT, i.e. with no residual pool of liquids. In this situation, higher values of the shear stress are expected, especially during the first stages of the simulation.

First, we note that, as depicted in Figure 14(d), a small hot metal layer is formed at the impingement point. Thus, in contrast with the prior cases, no significant air entrainment is predicted, and the phenomena associated with trapped slag droplets in the hot metal have also disappeared. Consequently, there is no need to filter the results to eliminate fluctuations in the volume fractions using the procedure detailed in Remark 5.1.

In Figure 23, the shear stress and pressure are depicted along the line crossing the jet impingement area for different time values. Again, as in Run A, a minimum of the shear stress curve in the stagnation point due to the jet impact on ${\mathrm{\Gamma }}_{wall}$ is observed, with non-symmetrical maximums around it. The time dependency is essentially negligible and the computed shear stress values are only slightly above those discussed in Section 5.2. This suggests that the role played by the pool of fluids in the impact mitigation is small. Nevertheless, when analyzing the shear stress and pressure contours across the entire impingement zone, displayed in Figure 24, we observe significantly higher shear stress values. The zone displaying high wall shear stress extends further, and the maximum computed shear stress is around 20% larger in Run A* compared to Run A. This maximum value is very similar to that found in Barral et al. (2021) for a planar hot metal jet on the dry trough. In both cases, no significant time dependency of the results was observed beyond the initial impact stage. However, it is worth noting that the computed pressure values in Run A* are notably higher, reaching $1.9\mathrm{e}4$ Pa, in contrast to the values observed in Runs A-D. In particular, directly comparing the pressure curves in Figure 22 with those in Figure 23, it is observed that the pressure in Run A* is around three times higher than that obtained in Run A.

Figure 23. Averaged wall shear (left) and pressure (right) on a line crossing the jet impingement area, for Run A*.

Figure 24. Wall shear stress and pressure on the jet impingement area for Run A*, $t=25\mathrm{s}$. (a) Shear stress (b) Pressure.

5.4. Effects of taphole stream velocity and diameter

Subsequently, we analyze the results for Runs C1–C5 (see Section 4.3). These cases involve variations in the taphole stream velocity and diameter while keeping a constant total slag volumetric fraction of $50\mathrm{\%}$, as specified in Table 3.

To properly capture the jet entry in the computational domain, all of these cases require the reconstruction of the mesh of the domain, except for Run C1, which serves as a baseline for comparison by maintaining the same parameter values as Run C. For the mesh reconstruction, the same sizes and features of Mesh 4 are preserved (see Section 4.1). In particular, to build the mesh for Runs C4-C5 and to ensure that the mesh is sufficiently fine in the vicinity of the jet trajectories, the calculations outlined in Section 3.1 are repeated to obtain their entry points and paths in the computational domain.

In Figure 25(a), the median height of the free surfaces is displayed as a function of the $x_3$ coordinate at time t = 50 s. Since the initial condition for Runs C1–C5 is the end value computed in Run C, the initial state of these free surfaces matches that shown in Figure 16. Hence, the time value t = 50 s would correspond to time t = 250 s in Section 5.2. Small differences are observed among the heights displayed for Runs C1, C2, and C3 compared to the initial free surface positions. Consequently, given that the only parameter that varies in these runs is $D_{tap}$, the numerical results suggest that the mass flow rate has minimal influence on the free surface evolution. Moreover, in Figure 25(b), the free surface median heights for Runs C4 and C5 are displayed. The height variations downstream of the impingement point are again minimal. However, due to the lower $v_{tap}$ values and the resulting different impact points, the behaviour of the free surfaces at the beginning of the BFT differs, resulting in a notable shrinkage of the hot metal pool in the stagnation region upstream of the jet impingement point.

Figure 25. Median height of the interfaces for Runs C1–C5, at t = 50 s. (a) Runs C1–C3 (b) Runs C4, C5.

In Figure 26, the velocity vectors on the plane $x_1=1.5\mathrm{m}$ are displayed near the impingement region. To facilitate comparison across cases, as in Figure 17, the maximum velocity magnitude contour in the legend is not the actual maximum, which is approximately equal to the taphole stream velocity. The displayed fields are averages over the time interval $(25,50]$ s. For Runs C1–C3, similar flow features are observed as in cases B-D, with complicated recirculation patterns induced by the jet impingement. The increased diameter in Runs C2 and C3 results in higher mass flow rates, leading to more pronounced circulations and higher velocities in the liquid pool. Nevertheless, the primary flow characteristics remain consistent. In Run C3, the region immediately following the jet impact point, where only a thin film of hot metal trails along the bottom wall of the trough, extends further downstream as a consequence of the larger mass flow rate. On the other hand, in Runs C4 and C5, most of the slag pool is within the backward-advancing recirculation generated by the impinging jet. The thickness of the hot metal layer that remains in the stagnation region upstream of the jet is very small, as already observed in Figure 25(b).

Figure 26. Velocity vectors on the plane $x_1=1.5\mathrm{m}$, with white contour lines corresponding to ${\alpha }_{slag}=0.5$ and velocity magnitude contours. All the results correspond to averages over $(25,50]$ s. (a) Run C1 (b) Run C2 (c) Run C3 (d) Run C4 (e) Run C5.

Away from the impingement zone, in the cases analyzed in this section, the flow of slag and iron closely resembles the results observed in Runs B-D, showing features akin to those typically observed in plug flow. The $v_3$ component of the velocity on the plane cross-section $x_3=11$ m is displayed in Figure 27. The flow behaviour is consistent across all runs, with velocity magnitude variations that can be attributed to the taphole mass flow rate values. In particular, Run C2 and C4, characterized by identical mass flow rates but differing taphole diameters and stream velocities, exhibit closely matched velocity values. In the latter part of the trough, the circulation patterns are identical to those presented for Run C in Figure 18(c).

Figure 27. Longitudinal component of the velocity ($v_3$) [m/s] averaged over $(25,50]$ s on the plane cross-section at $x_3=11$ m, only considering the part occupied by slag and hot metal.

Figure 28 shows the averaged wall shear stress contours on the jet impingement region at t = 50 s for the various analyzed cases. The wall shear is slightly larger in Run C3 than in Runs C1 and C2, with the high-wall shear zone being more spread. However, the differences in the observed values are relatively small, suggesting that the taphole diameter and, thus, the mass flow rate alone has a limited impact on wall shear. Conversely, in Runs C4 and C5, where the taphole stream velocity is reduced, significantly lower shear stress values are observed. Furthermore, the contour shapes exhibit distinct patterns, possibly due to the different impingement angle of the jet (see Figure 26). The maximum shear values depicted in the contours are gathered in Table 6, along with the maximum pressure values for each case. The pressure contour distribution closely resembles that presented in Figure 21 for Runs A-D, showing a pressure maximum in the jet impingement point.

Figure 28. Averaged wall shear stress on the wall on the jet impingement zone at t = 50 s, for Runs C1–C5.

Table 6. Maximum averaged wall shear stress and pressure for Runs C1–C5.

	Run C1	Run C2	Run C3	Run C4	Run C5
Max. shear stress [Pa]	725.5	779.5	744.1	469.2	550.5
Max. pressure [Pa]	$3.85\mathrm{e}4$	$4.77\mathrm{e}4$	$5.63\mathrm{e}4$	$5.98\mathrm{e}4$	$6.17\mathrm{e}4$

In Figure 29, we compare the averaged shear stress and pressure profiles for Runs C1–C5 at $t=50\mathrm{s}$, along the line dividing the jet impingement area. Due to the variations in the jet trajectory in Runs C4 and C5, the curves corresponding to these cases appear shifted relative to the corresponding impingement point. As in the previous results, the wall shear stress profiles show significant differences depending on the stream velocity. Nonetheless, the shapes of the profiles remain highly similar in all cases. Conversely, the pressure values exhibit substantial sensitivity to the taphole diameter but show remarkable similarity for the various analyzed taphole stream velocities. Furthermore, the wall shear is nearly identical in Runs C1–C3, as already discussed for Figure 28.

Figure 29. Averaged wall shear (left) and pressure (right) on a line crossing the jet impingement area at $t=50\mathrm{s}$, for Runs C1–C5.

The results presented in this section indicate that the taphole diameter alone exerts only a small influence on the primary flow characteristics within the trough. Specifically, keeping a constant velocity with larger taphole diameter, which results in higher mass flow rates, has a negligible impact on wall shear stress, free surface evolution, and flow characteristics in the BFT. On the other hand, variations in taphole stream velocity result in notable differences in shear stress and significant alterations in flow structures, but the positions of the free surfaces remain similar. This contrasts with the findings reported in Section 5.2, where it was observed that changing the taphole slag-hot metal ratio had a large influence on the wall shear stress and the evolution of the free surfaces.

As opposed to other previous studies dealing with the numerical simulation of the trough, such as Kou et al. (2019); Wang et al. (2017), we have found no significant slag entrainment with hot metal in any of the analyzed cases. This discrepancy can be attributed to the geometrical particularities of our trough configuration and jet impact point, since Kou et al. (2019) observed the former to have a substantial effect in the computed amount of entrained slag, finding less slag entrainment with deeper troughs. In particular, Wang et al. (2017) found a 100% separation efficiency of hot metal from slag when the iron dam and slag runner height difference was above 0.2 m. For the trough dimensions considered herein, this difference is precisely 0.2 m, as shown in Figure 3. Also, the different numerical approaches may play a substantial role in the analysis of separation. Since most of the entrainment is due to small droplets formed due to the intense shearing in the jet entry region (Ranjan et al., 2022), accurate simulation of this phenomenon may require very fine meshes and carefully considering surface tension effects. In any case, the settling velocity of the small droplets in the quiet zone that is generated downstream of the jet impact region can be estimated through Stokes' law (Cameron et al., 2019). Therefore, it is of interest to design BFTs that maximize the residence time of slag and hot metal, increasing the size of the quiet zone.

Experimental validation in the industrial setting of the BF trough poses significant challenges, particularly in directly measuring flow velocities. In the absence of direct measurements, we have compared our numerical findings with previous studies, particularly those involving physical scale models of the BFT (e.g. Begnis & Brandaleze, 2006; He et al., 2002; Monteiro de Oliveira et al., 2020). While the reported results are qualitatively well-aligned with the existing literature, further work should be aimed at validating the numerical findings, especially given the numerous simplifying assumptions within this study.

Many of these simplifications, such as considering continuous, Newtonian, and isothermal phases, or assuming a trough in perfect condition, have been used by other researchers in the field (see e.g. Monteiro de Oliveira et al., 2020). However, it remains uncertain whether these multiphase flow models precisely capture the industrial BFT's complexity, especially regarding velocity and turbulence parameters, as noted by Ranjan et al. (2022). In particular, further work should be also directed at assessing the extent to which thermal effects influence the hydrodynamics in the BFT (see Section 3.4).

Despite these limitations, the proposed mathematical model and the numerical methodology offer insights into understanding the flow characteristics within the BFT under various taphole conditions. The introduced methodology can serve as a practical tool for evaluating different trough geometrical configurations, which can minimize the region exposed to both slag and hot metal, reduce shear stress, and enhance slag-metal separation. Our findings can also be adapted for more complex studies. For instance, instead of simulating only the BFT hydrodynamics, the computed shear stress and pressure can be used in mechanical models (e.g, as demonstrated by Yao et al., 2021), and details such as the free surface positions can inform models aiming to assess the heat transfer in the BFT.

6. Conclusion

In this study, a transient 3D hydrodynamic model is solved to gain insights into the complex multiphase flow in a BF trough. Incorporating slag as an additional phase enabled us to investigate the flow behaviour under various taphole stream conditions, such as different slag fractions, taphole diameters, and stream velocities. The numerical simulations provide valuable insights into the phenomena occurring in the BFT that have not been addressed in the literature, such as the free surface evolution and the wall shear stress due to the jet impingement, where direct experimentation is difficult due to the harsh operating conditions.

In all the examined scenarios, intricate flow patterns characterized by large recirculations and return currents of both slag and hot metal are observed. The main features of the flow align qualitatively with previous studies, including both experimental work with physical scale models and numerical investigations by other authors. The results show that the complex flow in the trough is heavily influenced by the proportions of slag and hot metal in the taphole. Some of the most significant findings are summarized below:

• The simultaneous tapping of slag and hot metal leads to distinct flow patterns, deviating from the initial stage, when only hot metal is drained. When slag is also being tapped, the flow in the trough becomes increasingly similar to plug flow in the laminar region away from the impingement zone. In the surroundings of the impingement region, the kinetic energy of the impact leads to significant air and slag entrainment and mixing with hot metal in all cases. While the taphole diameter has only a minor impact on the primary flow characteristics, varying the taphole stream velocity results in differences in the jet impingement region due to the changes in its trajectory. However, the flow features in the laminar region away from the impact point remain similar.

• The region subjected to high wall shear stress is located on the bottom of the BFT, surrounding the impingement point. The computed shear stress values exhibit a pattern similar to studies of submerged impinging jets on planar surfaces, characterized by a minimum shear stress value directly at the jet impact point. Higher slag content or reduced velocity of the taphole stream results in decreased shear stress and a more confined high-shear zone near the jet impact point. For a 0.7 slag fraction, the shear stress in the jet impingement region is 31% lower than considering a hot metal jet. Taphole diameter changes, resulting in higher mass flow rates, have minimal influence on shear stress values. Jet impact in the dry trough leads to substantially increased shear stress and pressure values.

• Slag and hot metal levels in the trough evolve rapidly to quasi-steady states across all analyzed taphole flow conditions. The influence of taphole stream velocity and diameter on free surface dynamics is found to be limited. However, the depth of the slag and hot metal pools display large variations in response to changes in the taphole slag fraction. Specifically, a 40% relative increase between the height of the free surface separating slag and iron in the trough with respect to the initial value is found. Considering the slag pool depth, this value is slightly below 30%. These dynamic changes in the interface positions can have implications for refractory wear, as the location of the slag-hot metal line in the BFT is predicted to shift significantly during each tapping.

The contributions in this study offer valuable information on the various and complicated phenomena in the BFT, direct characterization of refractory wear remains an open challenge. The proposed methodology can serve as a practical tool for evaluating different trough geometrical configurations, aiming to minimize the region exposed to both slag and hot metal, reduce shear stress, or enhance slag-metal separation. The CFD model can also be coupled with other models in further work including the solid refractories, allowing for a comprehensive analysis of thermal and shear stress effects on the trough refractory linings. However, further work should be aimed at validating the numerical findings, especially given the numerous simplifying assumptions within this paper. In particular, the extent to which thermal effects influence the hydrodynamics in the BFT should be addressed, as well as the impact of more complicated taphole flow regimes.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was partially supported by ERDF and Xunta de Galicia funds under the ED431C 2017/60 and ED431C 2021/15 grants, by the Ministerio de Ciencia, Innovación y Universidades through the Plan Nacional de I+D+i (MTM2015-68275-R) and the grant BES-2016-077228, and by the Agencia Estatal de Investigación through project PID2019-105615RB-I00/ AEI /10.13039/501100011033.

Appendix. Turbulence model

In the SST $K-\omega$ model, Equations (16(16) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }K\right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }K\mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_K}})\text{grad}\right(K\left)\right]\\ & +G_K-Y_K,\end{array}$(16) ) and (17(17) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }\omega \right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }\omega \mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_{\omega }}})\text{grad}\right(\omega \left)\right]\\ & +G_{\omega }-Y_{\omega }+D_{\omega }.\end{array}$(17) ) are introduced to model the evolution of the turbulent kinetic energy, K, and the specific dissipation rate, ω, respectively (Menter et al., 2003). From the values of K and ω the turbulent viscosity, ${\hat{\mu }}_T$, is computed as (A1) $\begin{array}{r}{\hat{\mu }}_T={\displaystyle \frac{\hat{\rho }K}{max(\omega ,{\displaystyle \frac{S{F}_2}{a_1}})}},\end{array}$(A1) where S denotes the mean strain rate magnitude, defined as $S=\sqrt{2\mathbf{D}\left(\mathbf{v}\right):\mathbf{D}\left(\mathbf{v}\right)}$, and $a_1$ is a constant (see Table A1). Moreover, $F_2$ is the subsequent blending function: (A2) $\begin{array}{r}{F}_2=\tanh (max{({\displaystyle \frac{\sqrt{K}}{0.045\omega y_w}},{\displaystyle \frac{500\hat{\mu }}{\hat{\rho }{y}_w^2\omega }})}^2),\end{array}$(A2) being $y_w$ the distance to the nearest wall. Furthermore, we recall that due to the use of the VOF method to account for the different phases of the flow, $\hat{\rho }$ and $\hat{\mu }$ are the mixture properties corresponding to the density and dynamic viscosity, as introduced in (10(10) $\begin{array}{r}\hat{\varphi }={\alpha }_{slag}{\varphi }_{slag}+{\alpha }_{hm}{\varphi }_{hm}+(1-{\alpha }_{slag}-{\alpha }_{hm}){\varphi }_{air},\end{array}$(10) ).

Table A1. Turbulence model constants.

$a_1$	${\sigma }_{\omega 1}$	${\sigma }_{\omega 2}$	${\sigma }_{K1}$	${\sigma }_{K2}$	${\text{β}}_1$	${\text{β}}_2$	${\text{β}}_{\mathrm{\infty }}^{\ast }$	${\alpha }_1$	${\alpha }_2$	E	κ	$C_{\mu }$
0.31	2.0	1.168	1.176	1.0	0.075	0.0828	0.09	5/9	0.44	9.793	0.41	0.09

Source: Ansys, Inc. (2019).

An additional blending function, denoted as $F_1$, is also required. It is defined as (A3) $\begin{array}{r}{F}_1=\tanh (min{({\chi }_1,{\displaystyle \frac{4\hat{\rho }K}{{\sigma }_{\omega 2}{D}_{\omega }^{+}{y}_w^2}})}^4),\end{array}$(A3) where (A4) $\begin{array}{r}{\chi }_1=max({\displaystyle \frac{\sqrt{K}}{0.09\omega y_w}},{\displaystyle \frac{500\hat{\mu }}{\hat{\rho }\omega y_w^2}}),D_{\omega }^{+}=max(D_{\omega }^{\ast },1e-10),\end{array}$(A4) with ${\sigma }_{\omega 2}$ being the constant defined in Table A1. Moreover, $D_{\omega }^{\ast }$ is such that (A5) $\begin{array}{r}{D}_{\omega }^{\ast }=2\hat{\rho }{\displaystyle \frac{1}{{\sigma }_{\omega 2}}}{\displaystyle \frac{\text{grad}\left(K\right)\cdot \text{grad}\left(\omega \right)}{\omega }}.\end{array}$(A5)

In the conservation Equations (16(16) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }K\right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }K\mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_K}})\text{grad}\right(K\left)\right]\\ & +G_K-Y_K,\end{array}$(16) ) and (17(17) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }\omega \right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }\omega \mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_{\omega }}})\text{grad}\right(\omega \left)\right]\\ & +G_{\omega }-Y_{\omega }+D_{\omega }.\end{array}$(17) ), the turbulent Prandtl numbers ${\sigma }_K$ and ${\sigma }_{\omega }$ are defined as (A6) $\begin{array}{r}{\sigma }_K={\displaystyle \frac{1}{\frac{F_1}{{\sigma }_{K1}}+\frac{1-F_1}{{\sigma }_{K2}}}},{\sigma }_{\omega }={\displaystyle \frac{1}{\frac{F_1}{{\sigma }_{\omega 1}}+\frac{1-F_1}{{\sigma }_{\omega 2}}}}.\end{array}$(A6) In the previous equation, ${\sigma }_{K1}$, ${\sigma }_{K2}$, and ${\sigma }_{\omega 1}$ are model constants that have been calibrated with experimental measurements (see Table A1). The production terms of K and ω in (16(16) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }K\right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }K\mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_K}})\text{grad}\right(K\left)\right]\\ & +G_K-Y_K,\end{array}$(16) )–(17(17) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }\omega \right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }\omega \mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_{\omega }}})\text{grad}\right(\omega \left)\right]\\ & +G_{\omega }-Y_{\omega }+D_{\omega }.\end{array}$(17) ), denoted respectively as $G_K$ and $G_{\omega }$, are (A7) $\begin{array}{r}{G}_K=min({\hat{\mu }}_T{S}^2,10\hat{\rho }{\text{β}}_{\mathrm{\infty }}^{\ast }K\omega ),G_{\omega }={\displaystyle \frac{\hat{\rho }{\alpha }^{\ast }}{{\hat{\mu }}_T}}{G}_K,\end{array}$(A7) where ${\text{β}}_{\mathrm{\infty }}^{\ast }$ is a constant and ${\alpha }^{\ast }$ is such that (A8) $\begin{array}{r}{\alpha }^{\ast }={\alpha }_1{F}_1+{\alpha }_2(1-F_1).\end{array}$(A8) The limiter in the turbulent kinetic energy production term is introduced to avoid excessive build-up of turbulence in the stagnation regions (Ansys, Inc., 2019; Versteeg & Malalasekera, 2007). Additionally, in (16(16) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }K\right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }K\mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_K}})\text{grad}\right(K\left)\right]\\ & +G_K-Y_K,\end{array}$(16) )–(17(17) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }\omega \right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }\omega \mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_{\omega }}})\text{grad}\right(\omega \left)\right]\\ & +G_{\omega }-Y_{\omega }+D_{\omega }.\end{array}$(17) ), the turbulence dissipation terms $Y_K$ and $Y_{\omega }$ are defined as (A9) $\begin{array}{r}{Y}_K=\hat{\rho }{\text{β}}_{\mathrm{\infty }}^{\ast }K\omega ,Y_{\omega }=\hat{\rho }{\text{β}}^{\ast }{\omega }^2.\end{array}$(A9) Moreover, ${\text{β}}^{\ast }$ is defined using the blending function $F_1$ as (A10) $\begin{array}{r}{\text{β}}^{\ast }=F_1{\text{β}}_1+(1-F_1){\text{β}}_2,\end{array}$(A10) where ${\text{β}}_1$ and ${\text{β}}_2$ are constants taking the values collected in Table A1.

Lastly, the cross diffusion term, $D_{\omega }$, appearing in (17(17) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\left(\hat{\rho }\omega \right)}{\mathrm{\partial }t}}+\text{div}\left(\hat{\rho }\omega \mathbf{v}\right)=\text{div}\left[(\hat{\mu }+{\displaystyle \frac{{\hat{\mu }}_T}{{\sigma }_{\omega }}})\text{grad}\right(\omega \left)\right]\\ & +G_{\omega }-Y_{\omega }+D_{\omega }.\end{array}$(17) ), is such that (A11) $\begin{array}{r}{D}_{\omega }=(1-F_1)D_{\omega }^{\ast }.\end{array}$(A11)