Numerical simulation of the shear stress produced by the hot metal jet on the blast furnace runner

Computers and Mathematics with Applications 102 (2021) 146–159 Contents lists available at ScienceDirect Computers and Mathematics with Applications www.elsevier.com/locate/camwa Numerical simulation of the shear stress produced by the hot metal jet on the blast furnace runner P. Barral a,b, B. Nicolásc, P. Quintelaa,b,∗ aDepartamento de Matemática Aplicada and Instituto de Matemáticas (IMAT), Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain bTechnological Institute for Industrial Mathematics (ITMATI), Rúa Constantino Candeira s/n, 15782 Santiago de Compostela, Spain cDepartament Matemàtiques i Informàtica, Universitat de Barcelona, Gran Via de les Corts Catalanes 585, 08007 Barcelona, Spain A R T I C L E I N F O A B S T R A C T Keywords: Blast furnace Jet impact Numerical simulation Multiphase flow Turbulent flow Shear stress During steel casting process a jet of molten metal runs out of the blast furnace hearth and strikes the runner. The continuous impact of hot fluids causes significant damage to its surface, which is made of refractory concrete. In particular, the initial impact on the dry runner is expected to be critical. This work deals with the analysis of the mechanical impact on the runner through the numerical simulation of the process. We propose an incompressible turbulent isothermal Navier-Stokes model, where turbulence is modelled considering two 𝑘 −𝜔models (standard and SST). The interface dynamics is described by applying the Volume of Fluid (VOF) method, while the surface tension vector is provided by the Continuum Surface model (CSF). Their numerical results are performed in 2D. A comparative analysis of the most suitable transient turbulent multiphase model is presented by simulating benchmark physical experiments. The shear stress arising from the impact of the jet on the runner is also analyzed. An improvement of the classical analytical expression given in [1]is proposed. Both, the chosen turbulence model, and the formulas to compute the shear stress are validated using two benchmark laboratory tests and three numerical experiments. Numerical results are given for the impact of the jet on the dry runner of the blast furnace. 1. Introduction A blast furnace is a furnace used for smelting ores to produce industrial metals. This paper focuses on blast furnaces operating by a steelmaking company, similar to that considered in [2], where the smelted ores, iron and carbon (coke), are used to generate pig iron, also known as hot metal, which is the raw material of steel. At the top of the blast furnace there is the throat, see Fig. 1, through which the ores and fluxes are introduced. Hot air enriched with oxygen is injected under great pressure through the tuyeres, placed in the bosh. This air allows to hold the loads while the smelting process is performed. When the carbon ore contacts the air, several exothermic chemical reactions take place increasing the temperature inside the furnace, that can reach up to 1500 ◦C. When the iron ore is smelted, the bottom of the blast furnace, called hearth, is drilled and metal fluids run out through the taphole. Its inclination is around 10 degrees upwards (see Fig. 1). So the hot fluids run out of the blast furnace like a jet describing a parabolic path until reaching the runner, also known as trough (see Fig. 2). After approx- *Corresponding author at: Departamento de Matemática Aplicada and Instituto de Matemáticas (IMAT), Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain. E-mail address: [email protected] (P. Quintela). imately 60 minutes of casting, when the furnace gas starts going out, the taphole is plugged with clay. During the stop, the metal level in the furnace rises again. Once the necessary level is reached, a new casting cycle starts. In this way, the level of fluids in the hearth is kept as constant as possible for process safety [3]. So, the pressure inside the hearth and the velocity of fluids when running out through the taphole can be considered to be constant. Inside the furnace, in addition to the hot metal, there is slag, formed by fluxes, gangues of minerals and ashes of coke. The density of the hot metal is much greater than that of the slag, therefore the metal is located in the lower part of the furnace and is the first to go out when the taphole is open. Afterwards, a mixture of metal and slag bubbles goes out. Since the diameter of the taphole is very small, we can assume that it is a mixture. Later, they will be separated again on the casting runner because of their different densities. It should be noted that the runner is slightly inclined to encourage fluids to move towards the end of the runner. At the end of the runner, the slag at the top is deflected through the slag runner, while the hot metal goes down underneath the skimmer https://doi.org/10.1016/j.camwa.2021.10.013 Received 27 January 2021; Received in revised form 28 September 2021; Accepted 9 October 2021 Available online 22 October 2021 0898-1221/©2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). P. Barral, B. Nicolás and P. Quintela Computers and Mathematics with Applications 102 (2021) 146–159 Table 1 Nomenclature. Magnitude Description Magnitude Description 𝑎1Adjustment coefficient for SST turbulence model 𝛼,𝛼∗Adjustment coefficients for SST turbulence model 𝐛𝑝(m/s2) Body force per unit of mass acting on phase 𝑝𝛼 𝑝Volume fraction of 𝑝-phase 𝐛(m∕s2) Effective body force per unit of mass 𝛽,𝛽∗Adjustment functions for SST turbulence model 𝐶𝑙Cell 𝛿(adim) Dirac delta function 𝑑(m) Taphole diameter Δ𝑡(s) Time step 𝑑𝑛(mm) Nozzle diameter ΓSolid-fluid interface 𝐃(𝐯)(1∕s) Symmetrical part of the velocity 𝐯gradient Γaxis Boundary corresponding to the axis of the jet 𝐷𝜔(kg∕(m3s2)) Cross-diffusion term of 𝜔Γin,Γout,Γwall Inlet, outlet and wall boundaries 𝐹1,𝐹 2Blending functions for SST turbulence model Γ𝑘(kg∕(ms)) Effective diffusivity term of 𝑘 𝐹𝑟 (adim) Froude number Γ𝜔(kg∕(ms)) Effective diffusivity term of 𝜔 𝐠,𝑔 (m∕s2) Gravitational acceleration vector and modulus 𝜅(1∕m2) Interface curvature between two phases 𝐺𝑘(kg∕(ms3)) Generation term of 𝑘𝜆(adim) Relation between 𝑟and 𝐻 𝐺𝜔(kg∕(m3s2)) Generation term of 𝜔𝜇 𝑒𝑓𝑓 (Pa s) Effective fluid dynamic viscosity 𝐻(mm) Height from the nozzle to the wall 𝜇𝑓(Pa s) Generical fluid dynamic viscosity 𝐈(adim) Identity matrix 𝜇ℎ(Pa s) Hot metal dynamic viscosity 𝑘(m2∕s2) Turbulence kinetic energy 𝜇𝑇(Pa s) Turbulent fluid dynamic viscosity 𝑘𝑖𝑛𝑙𝑒𝑡 (m2∕s2) Inlet turbulence kinetic energy 𝜇𝑇𝑖𝑛𝑙𝑒𝑡 (Pa s) Inlet turbulent fluid dynamic viscosity 𝑘𝑑(m2s∕kg) Erosion kinetics coefficient per unit of mass 𝜋,Π,Π∗(Pa) Flow pressure, mean pressure and reduced pressure 𝑘𝑒𝑟 (s∕m) Erosion kinetics coefficient 𝜌(kg∕m3)Effectivemassdensity 𝐿w(mm) Size of the wall segment in normal jets simulation 𝜌𝑓(kg∕m3) Generical fluid mass density 𝑚 (kg∕(m2s)) Eroded mass flux 𝜌ℎ(kg∕m3) Hot metal mass density 𝐦𝜎(kg∕(m2s2)) Surface tension vector 𝜌𝑝(kg∕m3)𝑝-Phase mass density 𝑀(adim) Mach number 𝜌𝑠(kg∕m3) Generical solid porous media mass density 𝐌𝑝(N∕m3) Force of phase 𝑝due to interaction with other phases 𝜎(N∕m) Surface tension coefficient between two phases  𝐧(adim) Unit vector normal to the interface between two phases 𝜎ℎ,𝑎 (N∕m) Surface tension coefficient between hot metal and air 𝑁𝑐Cells number 𝜎𝑘,𝜎 𝜔Prandtl numbers for 𝑘and 𝜔 𝑁𝑝Number of immiscible phases 𝜎𝜔,2Adjustment coefficient for SST turbulence model 𝑟(mm) Radial coordinate for normal jets simulations 𝜏(Pa) Shear stress 𝑅𝑒 (adim) Reynolds number 𝜏𝑐(Pa) Critical shear stress 𝑆(1∕s) Strain rate magnitude 𝜏𝑚𝑎𝑥,𝜏 𝑚𝑖𝑛 (Pa) Maximum and minimum shear stress 𝑡(s) Time 𝝉𝑅(Pa) Reynolds stress tensor 𝐓𝑝(Pa) Stress tensor for phase 𝑝𝜔(1∕s) Specific turbulence dissipation rate 𝐓(Pa) Effective stress tensor 𝜔𝑖𝑛𝑙𝑒𝑡 (1∕s) Inlet turbulence dissipation rate 𝑣0(m∕s) Exit velocity from the nozzle ΩComputational domain for jet impact simulation 𝑣𝑠𝑜𝑢𝑛𝑑 (m∕s) Speed of sound Ω𝑒Computational domain for experiments simulation 𝑣Γ(m∕s) Solid-fluid interface modification velocity Ω𝑠𝑡𝑎𝑔 Stagnation zone 𝐯,𝐕,𝐯′(m∕s) Flow velocity, mean velocity and velocity fluctuation 𝐯𝑝(m∕s) Velocity of phase 𝑝 𝐯′𝑖𝑛𝑙𝑒𝑡 (m∕s) Hot metal velocity inlet fluctuation 𝐕𝑖𝑛𝑙𝑒𝑡 ,𝑣 𝑖𝑛𝑙𝑒𝑡 (m∕s) Hot metal velocity inlet, vector and modulus 𝑊𝑒(adim) Weber number 𝑥(m) Horizontal coordinate for hot metal simulation 𝑦(m) Vertical coordinate for hot metal simulation 𝑌𝑘(kg∕(ms3)) Dissipation term of 𝑘 𝑌𝜔(kg∕(m3s2)) Dissipation term of 𝜔 𝑧(m) Vertical coordinate for normal jets simulations Fig. 1. Schematic view of a blast furnace. and reaches the hot metal runner (see Fig. 2). Notice that after the skimmer there is a small vertical wall, so that hot metal level in this zone has to be high enough to go beyond it and to access to the hot metal runner. When a casting cycle finishes, the level of accumulated fluid is not enough for this to happen and a layer of fluid remains over the runner between two casting cycles. So, this hot fluid layer makes impossible to access the runner bottom surface until its useful life ends and, therefore, it is not feasible to take wear measurements while it is operating. For this reason, the numerical simulation of the jet impact is an essential tool. Useful life of a casting runner is limited by several wear phenomena that take place during casting process. Wear is mainly due to chemical, thermal and mechanical phenomena. Chemical wear is due to the composition of the slag, which contains highly corrosive substances. On the other hand, extremely high temperatures are reached during casting process and they can cause the appearance of fluid phases inside the runner refractory concrete. Furthermore, when the hot metal jet first hits the runner, a thermal shock could be produced, although in practice this is mitigated by the use of heaters for raising the temperature of the refractory concrete. Finally, mechanical wear has its origin in the impact produced on the runner by the molten metal jet when running out of the blast furnace under pressure. Our interest is focused on the latter phenomenon, the mechanical wear. A more detailed description of the entire casting process can be found in [3]. Although there is abundant literature on numerical simulation of furnace behaviour (see for example [4–7]), there are not many papers on numerical simulation of the runner. So, the authors have devoted considerable effort to understanding its thermomechanical behaviour: first performing a stationary two-dimensional thermal analysis in the middle section of the runner in [8]; then, a three-dimensional ther147 P. Barral, B. Nicolás and P. Quintela Computers and Mathematics with Applications 102 (2021) 146–159 Fig. 2. Schematic longitudinal view of the runner. Impact zone is surrounded and computational domain is marked with dash-line. mohydrodynamic model was solved to find the position of critical isotherms within solid refractory layers in [2]. In the latter and in the proceedings [9], a first numerical simulation of the impact of the hot metal jet falling from the blast furnace on the runner was presented. A simpler turbulence model, the Wilcox 𝑘-𝜔model, was used there [10]. A first objective of this work is to improve the jet hydrodynamic simulation. Since it is not possible to obtain experimental measures in blast furnace jets to validate numerical results, in this paper the authors present a comparative analysis of the Wilcox 𝑘-𝜔method with the Shear Stress Transport (SST) 𝑘-𝜔turbulence model on two benchmark laboratory tests registered in [1]. The analysis shows that the SST 𝑘-𝜔 turbulence model is better at simulating shear stress produced by a jet on a wall. In the previous papers [2,9]two scenarios were compared: the jet impact on dry trough versus on a narrow pool of hot metal. Numerical simulations shown that the maximum shear stresses produced by mechanical effects was obtained at the impact instant in the first scenario. Therefore the second objective of this work is to accurately compute the mechanical impact of the jet on the dry runner. More precisely, our interest is focused on the calculation of the exerted shear stress since its magnitude is considered in the bibliography the essential agent in the mechanical erosion, as can be seen in several references, from the earliest works by [11]to recent ones like the PhD thesis of [12]and her subsequent paper [13]. So, it is interesting to find an analytical formula to approximate them. In this sense, it is remarkable [1]work. They proposed a normalized shear stress equation adjustment from the experimental study of submerged air jets impinging normally over a rigid wall. Their equation has been widely applied, however, it presents some discrepancies with the laboratory experimental data as we move away from the impact point of centerline of the jet. In this paper we also propose a new adjustment equation to overcome these discrepancies. The proposed new formula is validated not only with the laboratory data from the benchmark cases mentioned above, but also with three additional numerical experiments. Numerical results with this method are presented and, in addition, a dynamic adaptive meshing is carried out that allows to improve the tracking of the jet and the calculation times. Section 2is devoted to the description of the jet impact problem. In order to carry out an accurate study of the impact, some preliminary analysis of the bibliography about erosion and turbulence models are performed in Section 3. Subsection 3.1 consists of an analysis of the mechanical erosion mechanisms, which allows to identify shear stress as the main magnitude to focus on. In Subsection 3.2, three benchmark laboratory tests are introduced and the classical formula of [1]is analyzed. Subsection 3.3 is devoted to propose an improvement of the later. In Subsection 3.4, a discussion of the most suitable turbulence model for this problem is included. Section 3is completed with the validation of the proposed formula for the calculation of the shear stress on three numerical tests. In Section 4, we come back to the jet impact problem arising in the casting process taking into account the conclusions drawn from the analysis of the previous section. Details about the involved models related to incompressible fluids, turbulence, and immiscible and multiphase fluids behaviour are given in its first subsections. The complete mathematical model associated to the jet impact problem is presented in Subsection 4.4. It corresponds to an incompressible turbulent multiphase flow in transient regime. Its numerical results can be found in Section 5. In particular, results are included on the shear stress with respect to the jet exit distance in the first moments of impact as well as the maximum pressure points. Finally, some conclusions about the main results of this paper are presented in Section 6. 2. Jet impact problem In this section we introduce the jet impact problem for its subsequent modelling and numerical simulation. Here, the explanation about how fluids run out of the blast furnace is presented, as well as their characterization in basis of the dimensionless numbers theory. To facilitate the reading of this document, a list of the notation used is included in Table 1. For this study, computational domain corresponding to dash-line marked box in Fig. 2is considered, focusing our attention on impact zone, surrounded by a circle in the same figure. In this paper we focus on the first few seconds of a new casting, since afterwards, as it was proved in [9]the liquid pool on the runner cushions the impact of the jet. Notice that during this time interval just hot metal runs out the blast furnace, so two phases are involved: hot metal and air. Since hot metal density is almost three times slag density, let us remark that in this way, we are considering the worst case from the impact point of view. As it was announced in the previous section, it is extremely difficult to measure the velocity at which the fluids run out the furnace, therefore it is assumed to be constant. To compute its value we have considered the following data: •Hot metal and slag production amount at each casting cycle: 450𝐸03 kg and 100𝐸03 kg, respectively. • Taphole diameter: 𝑑=0.06 m. • Material properties: 𝜌ℎ= 7015 kg∕m3,𝜇 ℎ=7.15𝐸−03Pas, where 𝜌ℎand 𝜇ℎare the hot metal density and its dynamic viscosity, respectively. •Casting cycle duration: 5400 s. All this information leads to the run out velocity modulus 𝑣𝑖𝑛𝑙𝑒𝑡 = 6.29 m∕s, named as velocity inlet, since this is the velocity at which hot metal enters the computational domain, shown in Fig. 2. Its value is considered constant and estimated using the mass flow computed from the previous data during the casting cycle, assuming that only hot metal comes out for the first 20 minutes and then a mixture of hot metal and slag comes out (see [14]). Dimensionless numbers study must be made in order to characterize the hot metal jet behaviour. Since we only know the data relative to hot metal in the taphole zone, the characterization is realized at the beginning of the fluid trajectory. Therefore, the characteristic velocity and length of the fluid considered are 𝑣𝑖𝑛𝑙𝑒𝑡 and the taphole diameter, 𝑑, respectively. More details about dimensionless numbers characterising fluids can be found, for example, in [15]. • Mach number, 𝑀, is the relation between the fluid characteristic velocity and the speed of sound, 𝑣𝑠𝑜𝑢𝑛𝑑 . Therefore, considering the air near the hot metal at a temperature of 30 ◦Cits value is: 𝑀=𝑣𝑖𝑛𝑙𝑒𝑡 𝑣𝑠𝑜𝑢𝑛𝑑 =0.018. 148 P. Barral, B. Nicolás and P. Quintela Computers and Mathematics with Applications 102 (2021) 146–159 Taking into account that the speed of sound in molten liquids is higher (of 4200 m∕sfor liquid Fe, see [16]), the Mach number is quite low for both phases, so the fluid can be considered incompressible with constant density. •Froude number, 𝐹𝑟, gives the relation between inertia and gravity terms: 𝐹𝑟=𝑣2 𝑖𝑛𝑙𝑒𝑡 𝑔𝑑 =67.21, where 𝑔(m∕s2) is the gravitational acceleration. The value of 𝐹𝑟is bigger than 1, so inertia forces of the hot metal jet overcome those of the gravity near the taphole. As the fluid leaves the taphole, the characteristic length increases and Froude number decreases. Therefore, gravity forces overcome inertia ones. So, the fluid describes a parabolic path. • Reynolds number, 𝑅𝑒, gives the relation between inertia and viscosity terms: 𝑅𝑒 =𝜌ℎ𝑣𝑖𝑛𝑙𝑒𝑡𝑑 𝜇ℎ =3.7𝐸+05. The value of 𝑅𝑒 is quite high and consequently, a turbulent flow must be considered. •To study surface tension importance in turbulent flows, Weber number, 𝑊𝑒, must be analysed. It gives the relation between inertia and surface tension terms: 𝑊𝑒=𝜌ℎ𝑑𝑣2 𝑖𝑛𝑙𝑒𝑡 𝜎ℎ,𝑎 =1.23𝐸+04, where 𝜎ℎ,𝑎 (N∕m) is the surface tension coefficient between hot metal and air, whose value is 𝜎ℎ,𝑎 =1.25 N∕m. Since Weber number is quite high, hot metal jet inertia overcomes its surface tension with the air. After this dimensionless analysis, we can conclude that the hot metal jet running out of the blast furnace consists of a transient problem for an incompressible, turbulent and multiphase flow, where surface tension forces do not play an important role. However, surface tension will be included in the mathematical model in order to achieve an accurate numerical simulation of the interaction between the jet and the thin pool of fluid reposing over the runner. Before introducing the mathematical model associated to jet behaviour, a previous analysis is included in the next section. It allows to answer two main questions: the relevant physical magnitude to characterize the runner wear and the most suitable choice of turbulence model. Once these doubts are clarified, in Section 4we proceed to present the mathematical model for the hot metal jet numerical simulation. 3. Previous analysis based on benchmark tests In this section we analyze the main erosion mechanisms, the magnitudes involved, and the more appropriated hydrodynamic models to carry out an accurate simulation of the actual jet impact. In Subsection 3.1, shear stress is identified as the main mechanical parameter responsible of running wear. In Subsection 3.2, three benchmark laboratory tests, introduced in [1]and [17], are included. An improvement of the classical formula given in [1]is proposed in Subsection 3.3. In Subsection 3.4, a discussion of the most suitable turbulence model in order to compute shear stress produced by a turbulent jet impacting over a rigid wall is included. For that, we analyze two turbulence models to simulate the physical experiments performed by [1]. The results show that Shear Stress Transport (SST) 𝑘 −𝜔turbulence model is the best choice. This section is completed with the validation of the proposed formula for the calculation of the shear stress on three numerical tests in Subsection 3.5. Fig. 3. Upper curve: Profile of the shear stress exerted by a jet impacting normally on the solid-fluid interface Γ. Lower curve: Profile of the theoretical modification of solid-fluid interface, Γ, according to classical erosion law, (1) or (2). 3.1. Mechanical erosion mechanisms In this section some erosion laws collected in the literature are analyzed in order to identify the mechanical parameters responsible of runner wear. As the runner is made of refractory concrete, we focus on the erosion on a porous solid (runner) produced by a fluid (hot metal). In this field, a large number of empirical and semi-empirical erosion laws and models have been developed, but their application is restricted to very specific scenarios. An erosion law that is applicable to our problem is the shear stress excess law, also known as classical erosion law, proposed by [11]. This law relates the eroded mass flux, 𝑚 (kg∕(m2s)), with the shear stress excess, (𝜏−𝜏𝑐)(Pa), exerted on the solid-fluid interface, Γ: 𝑚 ={𝑘𝑒𝑟(𝜏−𝜏𝑐),if 𝜏>𝜏 𝑐, 0,otherwise, on Γ,(1) being 𝜏𝑐(Pa) the critical shear stress, from which the erosion takes place, and 𝑘𝑒𝑟 (s∕m) the erosion kinetics coefficient. Both, 𝜏𝑐and 𝑘𝑒𝑟, only depend on the porous material. Eroded mass flux can be written in terms of the solid-fluid interface modification velocity, 𝑣Γ(m∕s): 𝑚 =𝑣Γ𝜌𝑠on Γ, being 𝜌𝑠(kg∕m3) the porous solid density. So, classical erosion law can be also written as 𝑣Γ={𝑘𝑑(𝜏−𝜏𝑐),if 𝜏>𝜏 𝑐, 0,otherwise, on Γ,(2) with 𝑘𝑑=𝑘𝑒𝑟∕𝜌𝑠. It is well known that the centerline of a jet has null velocity at the impact point, as well as null shear stress. As we move away from this point, these magnitudes increase until they reach a maximum, 𝜏𝑚𝑎𝑥 (Pa), and then they decrease again. So, if we look at the center plane of a jet with normal impact over a wall, we can see two symmetrical maxima and a minimum in between, see Fig. 3. The area between the two maxima is called stagnation zone, Ω𝑠𝑡𝑎𝑔 . Notice that the application of classical erosion laws, (1)or (2), translates into a peak of non-eroded material, as shown in Fig. 3, what in practice is not observed. To overcome this difficulty, [13] recently introduced the following modification to the classical erosion law: 𝑣Γ=⎧ ⎪ ⎨ ⎪ ⎩ 𝑘𝑑(𝜏𝑚𝑎𝑥 −𝜏𝑐),if 𝜏𝑚𝑎𝑥 >𝜏 𝑐in Ω𝑠𝑡𝑎𝑔, 𝑘𝑑(𝜏−𝜏𝑐),if 𝜏>𝜏 𝑐out of Ω𝑠𝑡𝑎𝑔, 0,otherwise, on Γ. Since most of the bibliography erosion laws point the shear stress as the responsible for the mechanical erosion process, this one is the magnitude selected to focus on and its accurate computation is the objective of the following subsections. 149 P. Barral, B. Nicolás and P. Quintela Computers and Mathematics with Applications 102 (2021) 146–159 Fig. 4. Scheme of a normal circular jet impacting over a wall. Table 2 Data of the experimental setup for experiments performed by [1]and [17]. Experiment 𝐻(mm) 𝑑𝑛(mm) 𝑣0(m∕s) BR-RUN No. 4 496.062 23.444 50.630 BR-RUN No. 5 422.275 6.426 89.365 BL 457.2 25.4 106.75 3.2. Benchmark laboratory tests In order to choose the most suitable turbulence model for studying shear stress exerted by a jet over a rigid wall, we analyse the results of simulating three laboratory experiments: two performed by [1], more precisely those named BR RUN No. 4and BR RUN No. 5, and one by [17], identified in the following by the BL acronym. These experiments consist on a circular turbulent jet composed by air, being 𝑑𝑛(mm) the nozzle diameter, from which the fluid goes down at a velocity 𝑣0(m∕s) and impacts normally over a rigid wall placed at a distance 𝐻(mm) from the nozzle, see Fig. 4. Data of the experimental setup are shown in Table 2. In [1], experimental measurements are presented as well as their adjustment equation for the shear stress exerted over the wall, 𝜏(Pa), which are normalized by the maximum shear stress, 𝜏𝑚𝑎𝑥 (Pa). A formula to approximate the maximum shear stress is also deduced. The adjustment equation they got is 𝜏 𝜏𝑚𝑎𝑥 =0.18(1−𝑒−114𝜆2 𝜆)−9.43𝜆𝑒−114𝜆2,(3) where 𝜆is an adimensional parameter defined as 𝜆 =𝑟∕𝐻, being 𝑟(mm) the distance to the jet center over the wall (see Fig. 4). Besides, maximum shear stress was approximated by: 𝜏𝑚𝑎𝑥 =0.16 𝜌𝑓𝑣2 0 (𝐻 𝑑𝑛)2,(4) where 𝜌𝑓=1.225 (Kg∕m3) is the fluid mass density. In their work, maximum shear stress was experimentally found at 𝜆𝑚𝑎𝑥 ≈0.14. For their shear stress adjustment equation, they also used the data obtained by [17]. Experimental data and adjustment curve, given by equation (3), are shown in Fig. 11 of [1]. There it is quite visible how the adjustment curve fits well for low values of 𝜆, but it does not for high values. Although equation (3)was criticized by other authors, see for example [18], the agreement among shear measurements for impinging jets Fig. 5. Experimental data obtained by [1]and by [17], together with the adjustment equation (3)and the improved adjustment equation (6). that are fully developed is acceptable, even compared with other experiments, especially for small 𝜆values. Notice that adjustment in equation (3)was introduced in the seventies, so it is totally normal to find some discrepancies, like the mentioned differences for high values of 𝜆or the fact that the maximum value for 𝜏∕𝜏𝑚𝑎𝑥 is 1.007 instead of the expected value 1, and that it occurs at 𝜆 =0.1371 instead of 0.14. Nowadays we have access to much more precise adjustment tools without requiring too much complication. These reasons have led the authors to look for a new formula to adjust the experimental data. According to classical erosion law (1), if exerted shear stress for high values of 𝜆is greater than critical shear stress of one material, this area should be eroded. So, it is important to propose a new adjustment equation to accurately approximate shear stress and therefore the subsequent wear. In order to improve their shear stress formula, experimental data were asked to Dr. N. Rajaratnam, but, unfortunately, the data collection was not available and he proposed us to recalculate the data points coordinates from Fig. 11 of their article. To achieve this goal, we used the free software Engauge Digitizer [19]. Using this software we localized the experimental points of that graphic for the three experiments, collected in Tables 3–5. Fig. 5shows the data of Tables 3–5and the goodness of the adjustment given by equation (3)(identified with the BR acronym) for the three benchmark cases considered. Besides, we can see its agreement with Fig. 11 of [1]. 3.3. Improved formula to calculate the normalized shear stress This subsection has the objective of proposing a new adjustment equation by fitting the three benchmark experimental tests. Following the mathematical profile of the adjustment equation proposed by [1], the following non linear curve to fit the data of Tables 3–5is proposed: 𝜏 𝜏𝑚𝑎𝑥 =𝑎(1−𝑒−𝑏𝜆2 𝜆)−𝑐𝜆𝑒−𝑑𝜆2,(5) where 𝑎, 𝑏, 𝑐and 𝑑are adjustment constants. Their computed values using MATLAB®are included in Table 6for each experiment along with their average and information about the goodness of fit. To propose the new adjustment shear stress equation, identified by the BNQ acronym, we use the average of the three computed values for each constant. BNQ normalized shear stress equation adjustment is the following one: 𝜏 𝜏𝑚𝑎𝑥 =0.2502(1−𝑒−54.49𝜆2 𝜆)−1.344𝜆𝑒−3.318𝜆2.(6) 150 P. Barral, B. Nicolás and P. Quintela Computers and Mathematics with Applications 102 (2021) 146–159 Table 3 Data obtained from Fig. 11 of [1]for the BR-RUN No. 4experiment. 𝜆𝜏∕𝜏𝑚𝑎𝑥 𝜆𝜏∕𝜏𝑚𝑎𝑥 𝜆𝜏∕𝜏𝑚𝑎𝑥 𝜆𝜏∕𝜏𝑚𝑎𝑥 0,01339 0,12089 0,09445 0,92072 0,16912 0,96832 0,27384 0,60710 0,03810 0,42174 0,10851 0,97103 0,18776 0,90817 0,29960 0,53967 0,05280 0,58277 0,11980 0,99640 0,20872 0,86532 0,35063 0,41244 0,06662 0,71514 0,13871 1,00688 0,22931 0,78237 0,38590 0,33022 0,08373 0,86102 0,15392 0,98855 0,25415 0,70536 0,41107 0,30475 Table 4 Data obtained from Fig. 11 of [1]for the BR-RUN No. 5experiment. 𝜆𝜏∕𝜏𝑚𝑎𝑥 𝜆𝜏∕𝜏𝑚𝑎𝑥 𝜆𝜏∕𝜏𝑚𝑎𝑥 𝜆𝜏∕𝜏𝑚𝑎𝑥 0,02758 0,29144 0,10342 0,9307 0,19057 0,91785 0,32527 0,50469 0,03484 0,4044 0,12225 0,9679 0,22259 0,81448 0,37532 0,38885 0,04835 0,4795 0,14730 0,9825 0,24554 0,73547 0,43049 0,30189 0,08428 0,83243 0,16438 0,9719 0,28425 0,61143 0,47513 0,25449 Table 5 Data obtained from Fig. 11 of [1]for the BL experiment. 𝜆𝜏∕𝜏𝑚𝑎𝑥 𝜆𝜏∕𝜏𝑚𝑎𝑥 𝜆𝜏∕𝜏𝑚𝑎𝑥 𝜆𝜏∕𝜏𝑚𝑎𝑥 0,03108 0,39467 0,08834 0,90419 0,16614 0,93477 0,30568 0,48368 0,04329 0,51169 0,09888 0,94288 0,18905 0,87198 0,33364 0,39155 0,05440 0,60194 0,10963 0,99304 0,22087 0,75428 0,39077 0,28084 0,06476 0,70455 0,13445 1,00858 0,24998 0,67271 0,44540 0,21771 0,07783 0,77199 0,15537 0,97717 0,27847 0,56248 0,50000 0,16603 Table 6 Values of constants obtained from fitting Eq. (5). SSE is the sum square of residuals and 𝑅2measures the success of the fit, being 1its ideal value. Experiment a b c d SSE 𝑅2 BR-RUN No. 4 0.2586 53.54 1.485 3.663 0.0048 0.996 BR-RUN No. 5 0.2525 51.24 1.197 3.196 0.0041 0.996 BL 0.2395 58.69 1.351 3.094 0.0056 0.996 Average 0.2502 54.49 1.344 3.318 The plot of the experimental data points together with BR adjustment equation (3)and BNQ adjustment equation (6)is shown in Fig. 5. We can observe that the new one fits very well, even for high values of 𝜆. 3.4. Turbulence model choice In order to carry out the numerical simulation of both experiments, in what follows, a steady problem is solved. Taking into account the axial symmetry of a normal circular jet, the simulations are performed using a two-dimensional domain, Ω𝑒=[0, 𝐿w] ×[0, 𝐻], being 𝐿wthe length of the wall segment considered for the simulation. This length is chosen long enough so that the numerical results are not affected by the condition imposed on the right boundary of the computational domain; a value of 𝐿w= 250 mm is assumed for both experiments. The software used for the simulations is ANSYS Fluent®, version 15.0. The boundaries of the geometry shown in Fig. 6are the following: •Γaxis (𝑟 =0mm and 𝑧 ∈(0, 𝐻)) corresponds to the jet axis. •Γin (𝑟 ∈(0, 𝑑𝑛∕2) and 𝑧 =𝐻), corresponds to the boundary through which the air gets into the domain at the velocity 𝑣0. •Γout (𝑟 ∈(𝑑𝑛∕2, 𝐿w)and 𝑧 =𝐻) ∪(𝑟 =𝐿wand 𝑧 ∈(0, 𝐻)), is a fictitious boundary considered to be open air. A pressure outlet condition is imposed on it. •Γwall (𝑟 ∈(0, 𝐿w)and 𝑧 =0mm) is the wall segment where the air jet impacts and shear stress is developed. Since our intention is to analyze shear stress on the wall, it is necessary to choose a turbulence model for an accurate treatment of boundary layer effects. In the following, two different turbulence models are compared: Standard 𝑘 −𝜔, proposed by [10], and one of its variants, Shear Stress Transport (SST) 𝑘 −𝜔model, proposed by [20]. The reaFig. 6. Computational domain Ω𝑒considered for the simulation of the experiments performed in [1]. sons for this election are two: the first one is that Standard 𝑘 −𝜔model1 was developed for wall-bounded shear flows, and the second one arises from the PhD thesis of [12], who concluded that, as time progresses, Standard 𝑘 −𝜔model provides better results than Standard 𝑘 −𝜖model (one of the most common turbulence models, developed for free shear flows); similar conclusions are reported in [21]. The main difference between the Standard 𝑘 −𝜔and the SST 𝑘 −𝜔turbulence models is that the last one combines the former near the walls with the Standard 𝑘 −𝜀 turbulence model far from them, and its use is recommended in predicting boundary layers under strong adverse pressure gradients (see [22]). As will be seen in this work, the profiles of the shear stress exerted by the pressurized air and hot metal jets are very similar, except for the lack of symmetry with respect to the centre of the jet, due to the discharge angle of the hot metal. Therefore, and given the unavailability of 1Abbreviation “Std.” is used in graphics and tables to designate Standard 𝑘 −𝜔turbulence model to save space. 151 P. Barral, B. Nicolás and P. Quintela Computers and Mathematics with Applications 102 (2021) 146–159 Table 7 Characteristics of the meshes used for the simulation of the laboratory experiments presented in [1]. Mesh 1 Mesh 2 Mesh 3 Mapped Yes Yes No Inflation No No Yes Axis 300 250 350 divisions Wall 150 + 100 + 200 + divisions Bias Factor 2Bias Factor 2Bias Factor 2 RUNNo.454545 No. elements 43,808 43,955 24,700 24,255 18,083 18,716 Table 8 Maximum shear stress values and the corresponding value of 𝜆. The first column corresponds to the estimated values by Beltaos and Rajaratnam equation (4), and the following ones to those obtained by the numerical simulations. BR-RUN No. 4 BR Eq. (4) Mesh 1 Mesh 2 Mesh 3 Std. SST Std. SST Std. SST 𝜏𝑚𝑎𝑥 (Pa) 1.122 2.916 1.166 2.567 1.051 2.215 1.042 𝜆≈0.14 0.059 0.131 0.065 0.128 0.078 0.133 BR-RUN No. 5 BR Eq. (4) Mesh 1 Mesh 2 Mesh 3 Std. SST Std. SST Std. SST 𝜏𝑚𝑎𝑥 (Pa) 0.363 0.496 0.344 0.468 0.301 0.260 0.387 𝜆≈0.14 0.093 0.157 0.101 0.146 0.087 0.140 Fig. 7. Results of the BR-RUN No. 4simulations, using the three meshes shown in Table 7and the Standard (Std.) and SST 𝑘 −𝜔turbulence models. Left: Shear stresses exerted on the wall. Right: Normalized shear stresses. experimental data for the hot metal, these benchmarks allow a proper evaluation of which of the methods will best capture the shear stress in the real case. In the numerical simulation, three different meshes are used, with the characteristics shown in Table 7. Term Inflation refers to a type of mesh refinement that generates layers of parallel elements on a surface (or line, in 2D). It is usually used near the walls when there is a boundary layer, as the case of the turbulent boundary layer. It is highly recommendable to use inflations when generating no mapped meshes, like in Mesh 3, where a refinement of 10 layers with a growth rate of 1.05 among them in 5mm of total thickness is considered. Bias Factor is the ratio of the largest edge to the smallest one of the mesh. The laboratory benchmark tests have been numerically reproduced using commercial software ANSYS Fluent®, version 15.0 for the three meshes. The results obtained are shown in Figs. 7and 8. In spite of not having the shear stress experimental measures, we compare the results of the simulations with the shear stress calculated by equation (3) and their maximum value, calculated from equation (4). All shear stress curves have a minimum at 𝑟 =0mm. This minimum corresponds to the impact point of the jet centerline, where the velocity is null, as well as shear stress, what is in good agreement with real observations. As we move away from this point, shear stress increases until reaching a maximum value and then decreases again as we have previously announced. It is quite clear that the results obtained using the SST 𝑘 −𝜔turbulence model are closer to those predicted by equation (4). Furthermore, the improved fitting curve introduced in (6)is closer to the results of the numerical simulations with the SST 𝑘 −𝜔turbulence model than those obtained from (3), even for 𝜆values far from the centre of the jet. In addition, in the Figs. 7and 8, we can see how the Standard 𝑘 −𝜔turbulence model overestimates the shear stress, and it is more dependent on the mesh than the SST 𝑘 −𝜔turbulence model. In Table 8maximum shear stress values and the corresponding value of 𝜆for each experiment, turbulence model and mesh are summarized. We can see how Mesh 3 provides good results for both experiments, in spite of being the most coarse of the three ones. The reason is the Inflation application, that allows to solve accurately the boundary layer. In conclusion, to model the behaviour of the hot metal jet, the SST 𝑘 −𝜔turbulence model seems to be the most appropriated. In Section 4, a detailed description of its associated mathematical model is presented. 3.5. Numerical validation of the shear stress calculation in other scenarios To make sure that the new adjustment equation for normalized shear stress is verified in other scenarios of jet impact, three new numerical experiments are proposed. They are similar to those performed in the laboratory experiments described above, although the jet nozzle height, its diameter, and the output speed are modified. The values of 𝐻, 𝑑𝑛 152 P. Barral, B. Nicolás and P. Quintela Computers and Mathematics with Applications 102 (2021) 146–159 Fig. 8. Results of the BR-RUN No. 5simulations, using the three meshes shown in Table 7and the Standard (Std.) and SST 𝑘 −𝜔turbulence models. Left: Shear stresses exerted on the wall. Right: Normalized shear stresses. Table 9 Details of the numerical experiments data performed to validate the new adjustment equation (5). Experiment 𝐻(mm) 𝑑𝑛(mm) 𝑣0(m∕s) No. elements Exp. 1 500 20 50 18,026 Exp. 2 400 10 100 19,456 Exp. 3 350 20 60 22,931 and 𝑣0for each numerical experiment are shown in Table 9. The magnitude of these values incorporates substantial changes with respect to laboratory data, although the jets are still considered to be air. Notice that in these cases, physical experiments were not performed, but their simulation using the commercial software ANSYS Fluent®version 15.0. For the numerical simulations, we have used the SST 𝑘 −𝜔turbulence model. As done in the previous subsection, we consider the bidimensional domain described in Fig. 6for the values given in Table 9, and perform each numerical experiment in steady regime using a mesh based on Mesh 3 specifications (see Table 7), such that the number of elements corresponds to the last column in Table 9for each experiment. Normalized shear stress obtained for the three trials as well as the adjustment curves proposed for BR and for us are shown in Fig. 9, for each performed simulation. The new normalized shear stress equation adjustment (6)fits very well the normalized numerical shear stress profile for the three numerical experiments. It also fits better than the classical formula (3)given by [1]. 4. Mathematical model In this section, we come back over the blast furnace jet problem introduced in Section 2. A complete mathematical model for an incompressible multiphase turbulent flow in transient regime is explained. After the analysis presented in Subsection 3.4, SST 𝑘 −𝜔turbulence model is considered for the turbulence problem. The multiphase problem is solved by applying the Volume of Fluid (VOF) method and the Continuum Surface Force (CSF) model. For sake of simplicity, we consider the bidimensional problem posed on the longitudinal central section of the runner, inside the area marked with dash-line in Fig. 2. Computational domain Ωis presented in Fig. 10, and their boundaries are: •Γin corresponds to the end of the taphole, through which fluids run out from the blast furnace and enter the computational domain. It is placed 1.27 m above the bottom runner surface and has a diameter of 0.06 m. •Γout is the part in contact with the air. It is a fictitious boundary, which is considered far enough from the jet in order to not to affect its numerical simulation. •Γwall corresponds to the bottom runner surface, defined by the points: (0, 0.8), (0.615, 0.8), (1.615, 0.07) and (15.965, 0), and its final vertical wall. 4.1. Incompressible model An incompressible flow behaviour in laminar regime under the force of gravity is described by Navier-Stokes equations for mass and momentum balances: ⎧ ⎪ ⎨ ⎪ ⎩ div𝐯=0, 𝜌𝑓 𝜕𝐯 𝜕𝑡 +𝜌𝑓div(𝐯⊗𝐯)+grad𝜋−2𝜇𝑓div(𝐃(𝐯)) = 𝜌𝑓𝐠,(7) where 𝐯(m∕s) is the flow velocity, 𝐃(𝐯)(s−1) is the symmetrical part of its gradient, 𝜋(Pa) is the flow pressure, 𝜇𝑓(Pa s) is the fluid dynamic viscosity, 𝑡(s) is the time and 𝐠(m∕s2) is the gravitational acceleration vector. Operator ⊗denotes tensor product of two vectors. For more details, see for example [23]. 4.2. Turbulence model Turbulent flows are developed in very different size scales and a direct numerical simulation would be extremely expensive from a computational point of view. So, magnitudes of interest (for example, the velocity 𝐯) are usually decomposed into their mean value (denoted by 𝐕or ⟨𝐯⟩) and their fluctuation (𝐯′): 𝐯=𝐕+𝐯′. Decomposed magnitudes are replaced in Navier-Stokes equations, (7), and the resulting equations are averaged, obtaining the following ones: ⎧ ⎪ ⎨ ⎪ ⎩ div𝐕=0, 𝜌𝑓 𝜕𝐕 𝜕𝑡 +𝜌𝑓div(𝐕⊗𝐕)+𝜌𝑓div(⟨𝐯′⊗𝐯′⟩)+gradΠ−2𝜇𝑓div(𝐃(𝐕)) =𝜌𝑓𝐠, (8) where operator Π(Pa) denotes the mean pressure [23]. It must be observed that averaged equations are analogous to those for a laminar incompressible flow, (7), in terms of averaged velocity and gradient. There is only a new term, proportional to the tensor product of velocity fluctuations, named Reynolds stress tensor (𝝉𝑅): 𝝉𝑅=−𝜌𝑓⟨𝐯′⊗𝐯′⟩. This term introduces six new unknowns, while the number of equations remains the same. To overcome this situation, it is necessary to introduce the turbulence model. 153 P. Barral, B. Nicolás and P. Quintela Computers and Mathematics with Applications 102 (2021) 146–159 Fig. 9. Normalized shear stress for numerical experiments using the SST 𝑘−𝜔turbulence models. Top left: Exp. 1. Top right: Exp. 2. Bottom: Exp. 3. Fig. 10. Computational domain Ωfor jet impact problem. The Shear Stress Transport (SST) 𝑘 −𝜔turbulence model is within the first order Reynolds Averaged Navier Stokes (RANS) models. These models use Boussinesq hypothesis to approximate the Reynolds stress tensor. According to this hypothesis, this tensor is assumed to be similar to viscous stress tensor: 𝝉𝑅=2𝜇𝑇𝐃(𝐕)+ 1 3tr(𝝉𝑅)𝐈, being 𝜇𝑇(Pa s) the turbulent dynamic viscosity and 𝐈the identity matrix. This expression corresponds to the decomposition of a tensor into its deviatoric and spherical parts, being the deviatoric part proportional to 𝐃(𝐕), with a proportionality factor of 2𝜇𝑇. The spherical part is related to the turbulence kinetic energy 𝑘(m2∕s2), that is defined as 𝑘=1 2⟨|𝐯′|2⟩=− 1 2𝜌𝑓 tr(𝝉𝑅). Taking into account these assumptions, averaged Navier-Stokes equations presented in (8) can be written as: ⎧ ⎪ ⎨ ⎪ ⎩ div𝐕=0, 𝜌𝑓 𝜕𝐕 𝜕𝑡 +𝜌𝑓div(𝐕⊗𝐕)+gradΠ∗−2𝜇𝑒𝑓𝑓 div(𝐃(𝐕)) = 𝜌𝑓𝐠,(9) where Π∗=Π −tr(𝝉𝑅)𝐈∕3 is the reduced mean pressure, and 𝜇𝑒𝑓 𝑓 = 𝜇𝑓+𝜇𝑇is the effective dynamic viscosity. In addition, SST 𝑘 −𝜔turbulence model consists of two transport equations that have the objective of finding the turbulent dynamic viscosity, 𝜇𝑇[20]. One equation is proposed for 𝑘and another one for the specific turbulence dissipation rate, 𝜔(1∕s), that is proportional to 𝜌𝑓𝑘 𝜇𝑇 through the relation: 𝜇𝑇=1 𝐴 𝜌𝑓𝑘 𝜔, where 𝐴is a function computed as: 𝐴=max[1 𝛼∗,𝑆𝐹2 𝑎1𝜔], being 𝛼∗and 𝑎1model constants experimentally obtained, 𝑆(1∕s) the strain rate magnitude, and 𝐹2a blending function that depends on variables 𝑘and 𝜔, fluid properties, and the vertical distance to the wall. Notice that this function allows to switch from Standard 𝑘 −𝜔 turbulence model formulation to Standard 𝑘 −𝜀one, and vice versa, depending on the distance to the wall [20]. Transport equations for 𝑘and 𝜔are the following: ⎧ ⎪ ⎨ ⎪ ⎩ 𝜌𝑓 𝜕𝑘 𝜕𝑡 +𝜌𝑓𝐕⋅grad𝑘=div(Γ𝑘grad𝑘)+𝐺𝑘−𝑌𝑘, 𝜌𝑓 𝜕𝜔 𝜕𝑡 +𝜌𝑓𝐕⋅grad𝜔=div(Γ𝜔grad𝜔)+𝐺𝜔−𝑌𝜔+𝐷𝜔, (10) where Γ𝑘(kg∕ms) and Γ𝜔(kg∕ms) are the effective diffusivity terms of 𝑘and 𝜔, 𝐺𝑘(kg∕ms3) and 𝐺𝜔(kg∕m3s2) are the generation terms of 𝑘and 𝜔, 𝑌𝑘(kg∕ms3) and 𝑌𝜔(kg∕m3s2) are the dissipation terms of 𝑘 and 𝜔due to turbulence, and 𝐷𝜔(kg∕m3s2) is the cross-diffusion term, that depends on 𝑘and 𝜔gradients. These terms have the following expressions: 154