# **Modeling and Simulation of Metallurgical Processes in Ironmaking and Steelmaking**

Edited by Thomas Echterhof, Ko-Ichiro Ohno and Ville-Valtteri Visuri Printed Edition of the Special Issue Published in *Metals*

www.mdpi.com/journal/metals

## **Modeling and Simulation of Metallurgical Processes in Ironmaking and Steelmaking**

## **Modeling and Simulation of Metallurgical Processes in Ironmaking and Steelmaking**

Editors

**Thomas Echterhof Ko-Ichiro Ohno Ville-Valtteri Visuri**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Thomas Echterhof RWTH Aachen University Germany

Ko-Ichiro Ohno Kyushu University Japan

Ville-Valtteri Visuri University of Oulu Finland

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Metals* (ISSN 2075-4701) (available at: https://www.mdpi.com/journal/metals/special issues/modeling simulation ironmaking steelmaking).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-5153-1 (Hbk) ISBN 978-3-0365-5154-8 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**

## **About the Editors**

#### **Thomas Echterhof**

Thomas Echterhof is an academic director and deputy head of the Department for Industrial Furnaces and Heat Engineering at RWTH Aachen University, Germany. He received his Ph.D. (Dr.-Ing.) in engineering science in 2010. He is heading the research group on mass and energy balances with a special emphasis on the sustainability of electric arc furnace steelmaking including circular economy aspects. His main research interests include the optimization of energy and resource efficiency and the environmental impact of energy and resource-intensive processes such as EAF steelmaking. He is also involved in teaching and is giving a lecture on electric arc furnace technology at RWTH Aachen University. He has authored or co-authored more than 70 scientific papers in journals and conference proceedings.

#### **Ko-Ichiro Ohno**

Ko-ichiro Ohno is a professor in the Department of Materials in the Faculty of Engineering at Kyushu University. Before his current position, he worked at Kyushu University, as an assistant professor from 2006 to 2012. At the same time, stayed at RWTH Aachen University, Germany, as a postdoctoral researcher from 2010 to 2011. He worked at the Kyushu university as an associate professor from 2013 to 2021. Ohno´s field of specialization is metallurgical process engineering based on kinetics, transport phenomenology, and thermodynamics. His research fields are kinetics analysis and thermodynamics of the ironmaking process by "in situ" observation technique. The future target is to develop a novel technology in the iron-making process for the reduction of CO2 emissions through our study and education. He has authored or co-authored over 80 peer-reviewed journal and conference papers. Ohno holds a doctoral degree in engineering as well as master's and bachelor's degrees from Tohoku University.

#### **Ville-Valtteri Visuri**

Ville-Valtteri Visuri is an associate professor and deputy head of the Process Metallurgy Research Unit at the University of Oulu. Before his current position, he worked as an R&D manager at Outokumpu Stainless Oy and as a senior research fellow and postdoctoral researcher at the University of Oulu. He has also had two research visits to RWTH Aachen University. Visuri's research is focused on the chemical and physical fundamentals of steelmaking processes, particularly the modeling of unit processes. He has authored or co-authored 59 peer-reviewed journal and conference papers. Alongside his main occupation, he serves as the chairman of the board in the Division of Metallurgical Engineers at the Finnish Association of Mining and Metallurgical Engineers. Visuri holds the title of Docent in Process Metallurgy, a doctoral degree in process engineering as well as master's and bachelor's degrees in Industrial Engineering and Management, all from the University of Oulu.

## **Preface to "Modeling and Simulation of Metallurgical Processes in Ironmaking and Steelmaking"**

This Special Issue focuses on the modeling and simulation of ironmaking and steelmaking processes, providing a selection of peer-reviewed papers on topical research questions. The idea for co-editing a Special Issue was borne out of the bilateral collaborations between RWTH Aachen University, Kyushu University and the University of Oulu, which have included both experimental and modeling aspects.

The manuscript submission was opened amidst the COVID-19 pandemic in August 2020 and closed in March 2022 with still no end to the pandemic in sight. For this reason, all aspects of the Special Issue had to be organized without meeting in person and face-to-face marketing of the Special Issue was limited. Nevertheless, we managed to receive good visibility for our Special Issue and were thus able to attract a fair number of manuscript submissions for peer review. Eventually, a total of 14 manuscripts were accepted for publication based on positive recommendations from external subject matter experts. These papers and an accompanying editorial were published from 2021 to 2022 in volumes 11 and 12 of *Metals* and are reprinted in this Special Issue.

We are grateful to both the peer reviewers and authors, without whom this Special Issue could not have materialized. Mr. Toliver Guo from MDPI is acknowledged for taking care of all the practicalities related to the publication of the Special Issue both in digital and printed form. Finally, we hope that the Special Issue will make interesting reading and give our readers new ideas and tools for driving scientific and practical progress in the iron and steel industry.

> **Thomas Echterhof, Ko-Ichiro Ohno, and Ville-Valtteri Visuri** *Editors*

## *Editorial* **Modeling and Simulation of Metallurgical Processes in Ironmaking and Steelmaking**

**Thomas Echterhof 1,\*, Ko-Ichiro Ohno 2,\* and Ville-Valtteri Visuri 3,\***

#### **1. Introduction and Scope**

The UN's 2030 Sustainable Development Goals, the Paris Agreement, and the European Green Deal, among other goals, all aim to improve the sustainability of industrial production and reduce CO2 emissions. The European Union, for example, aims to reach carbon neutrality and a circular economy by 2050. This goal cannot be achieved without a significant reduction in the CO2 emissions created by the steel industry. To reach this goal, further process optimizations regarding energy and resource efficiency, as well as the development of new processes or process routes, are needed.

The parameters necessary for the analysis and optimization of the existing and new metallurgical processes in ironmaking and steelmaking often cannot be measured directly because of the harsh conditions inside the furnaces and metallurgical vessels. Typically, the direct information sources in ironmaking and steelmaking are off-gas analysis and spot measurements, which suffer from the delay associated with the analysis of the sample. Due to the harsh environment, opportunities to determine the flow conditions in the vessels by measurements are even more limited. While new methods for the direct and continuous measurement of some of these parameters are under development, for many processes, they are currently unavailable. Furthermore, plant trials that would be necessary to evaluate the impact of different optimization strategies may be impossible because of the prohibitive cost or safety concerns.

Modelling and simulation have established themselves as an invaluable source of information regarding otherwise unknown process parameters, and as an alternative to plant trials with a lower associated cost, risk, and duration. Models are also applicable for model-based control of metallurgical processes.

This Special Issue aimed to cover recent advances in the modelling and simulation of unit processes in ironmaking and steelmaking. To this end, fourteen articles have been published in the present Special Issue of *Metals*. The articles give voice to a total of 67 authors, representing affiliations from ten countries (Figure 1a) with the majority of contributions originating from universities (Figure 1b). The subjects include reviews on the fundamentals of modelling and simulation of metallurgical processes, as well as the fields of iron reduction/ironmaking, steelmaking via the primary and secondary route, and continuous casting.

**Citation:** Echterhof, T.; Ohno, K.-I.; Visuri, V.-V. Modeling and Simulation of Metallurgical Processes in Ironmaking and Steelmaking. *Metals* **2022**, *12*, 1185. https://doi.org/ 10.3390/met12071185

Received: 6 July 2022 Accepted: 6 July 2022 Published: 12 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** Authors of the contributions of the Special Issue: (**a**) by country of affiliation; (**b**) by affiliation type.

#### **2. Contributions**

The fundamentals of modelling and simulation are covered by several articles. Wang et al. [1] review the application of turbulence modelling in metallurgical applications. The article focuses on the validation of models by experiments and the applicability of models to industrial cases. Bubble-induced turbulence, supersonic jet transport, and electromagnetic suppression of turbulence are three specific flow problems discussed in the article. Haas et al. [2] present "A Review of Bubble Dynamics in Liquid Metals". The article presents a detailed discussion of the state of knowledge of bubble dynamics in liquid metals, measurement methods, and bubble formation mechanism at nozzles and purging plugs, and discusses modelling uncertainties using the steel ladle as an example. It is concluded that, while a lot of valuable information has been extracted from both aqueous and liquid metal systems, there are still plenty of uncertainties regarding the behaviour of bubbles in liquid steel. The uncertainties highlighted in the paper include the effect of impurities present in steel, the interaction of gas bubbles in bubble swarms, and the deformation of bubbles due to injectors or walls.

Three articles explore the reduction processes in ironmaking. Quatravaux et al. [3] investigated the Midrex NGTM process and adapted a blast furnace operating diagram to describe the direct reduction of iron oxide in the Midrex shaft furnace. The authors employed the developed graphical tool for a sensitivity study of the reduction in the shaft furnace. Wang et al. [4] present a combination of numerical and experimental investigations of the charging of carbon composite briquettes in a blast furnace. The authors studied the replacement of ore by the carbon composite briquettes and report on carbon conversion, coke rate, blast furnace operation, and productivity, etc. Liu et al. [5] studied the reduction process within a pre-reduction rotary kiln employing numerical simulation methods. The numerical model was validated based on measurement data and then used to investigate process details such as temperature gradients or the progress of the reduction process within the rotary kiln.

The field of steelmaking is covered by eight articles. The first two are more related to the primary steelmaking route, including the blast furnace and basic oxygen furnace while the following five concern steelmaking in the electric arc furnace (EAF). The final article addresses the overarching topic of steelmaking slags. Linnestad et al. [6] report on the development of a process model describing the Composition Adjustment by Sealed argon bubbling with Oxygen Blowing (CAS-OB) process. In this study, a simplified mathematical model for the heat-up stage of the CAS-OB process was coupled with adaptive estimation of state variables using a Kalman filter. The authors present the model's capabilities, as well as model results from the prediction of a data series with more than 1000 heats based on the model implementation at two CAS-OB stations at SSAB Europe Oy in Raahe, Finland. The predicted temperature profiles are displayed for the operator to assist in their decision-making, thus helping to reduce the number of rejected heats. Gao et al. [7] investigated the melting of scrap in hot metal and determined mass transfer coefficients during scrap melting based on laboratory-scale experiments. Subsequently, the authors determined the influence of the variables molten pool stirring rate, bath temperature, and scrap type on the mass transfer coefficient. Tomažiˇc et al. [8] studied the optimization of the energy consumption in an EAF employing data-driven modelling approaches. The models were developed and validated based on industrial data from an EAF steel plant and were used to determine the optimal duration of the transformer profile during melting to reduce energy consumption, thereby increasing EAF efficiency. Jawahery et al. [9] present a first-principles-based process model of the EAF steelmaking designed for online real-time optimization of the process. The model focuses on the integration of auxiliary process data to predict energy efficiency and heat transfer limitations in the EAF. The model was validated using steel temperature and weight measurements achieving reasonable agreement. Al Nasser et al. [10] developed a simplified arc impingement model to study the direct-current electric arc in computational fluid dynamics simulations. The model was used to investigate the influence of the factors arc gap, the density of the gas and total electric current on the behaviour of the arc, and the overall process, including arc impingement depth, velocity magnitude, and arc stability. Reimann et al. [11] investigated statistical modelling approaches based on operational data from five industrial-size EAFs for the prediction of the electric energy demand. Results of the modelling approaches applied to the different industrial EAFs are presented, compared and differences are discussed. Schubert et al. [12] describe the development of a fast modelling approach that was subsequently used to simulate and predict the scrap preheating of a continuously charged EAF. Modelling applications, assumptions, possible enhancements, and limitations, as well as initial simulation results, are presented and discussed. Finally, Safavi Nick et al. [13] studied the modelling and simulation of slag heat recovery to design an optimized heat recovery system. The simulation focused on the heat exchanger and fluid medium used to transfer the heat. This study is a good example of how modern computational fluid dynamics (CFD)-based simulation tools can be used as a design tool for heat recovery processes.

The final field – covered by one article – is continuous casting. Guthrie and Isac [14] present a historical review of casting methods used for sheet steel production. The review shows the development from conventional continuous casting to thin-slab casting as well as different concepts of direct strip casting, and discusses the pros and cons of current methods.

#### **3. Conclusions and Outlook**

A variety of connected topics have been compiled in the present Special Issue of *Metals*, providing an overview of recent developments in different aspects of modelling and simulation of metallurgical processes in ironmaking and steelmaking. Despite being a longterm research interest, modelling and simulation could become even more important with the upcoming transformation to carbon-neutral ironmaking and steelmaking processes.

As Guest Editors of this Special Issue, we are very pleased with the contributions received and hope that the presented articles will be useful to researchers and operators of metallurgical plants working to optimize plant operation. We would like to thank all the authors for their contributions, as well as all reviewers for their efforts in maintaining the high quality of published papers. We would also like to give special thanks to all staff at the *Metals* Editorial Office, especially to Toliver Guo, Assistant Editor, who managed and facilitated the publication process of this Special Issue.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Review* **Mathematical Methodology and Metallurgical Application of Turbulence Modelling: A Review**

**Yannan Wang, Lingling Cao \*, Zhongfu Cheng, Bart Blanpain and Muxing Guo**

Department of Materials Engineering, Katholieke Universiteit (KU) Leuven, 3000 Leuven, Belgium<EMAIL_ADDRESS>(Y.W<EMAIL_ADDRESS>(Z.C<EMAIL_ADDRESS>(B.B<EMAIL_ADDRESS>(M.G.)

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** This paper focusses on three main numerical methods, i.e., the Reynolds-Averaged Navier-Stokes (RANS), Large Eddy Simulation (LES), and Direct Numerical Simulation (DNS) methods. The formulation and variation of different RANS methods are evaluated. The advantage and disadvantage of RANS models to characterize turbulent flows are discussed. The progress of LES with different subgrid scale models is presented. Special attention is paid to the inflow boundary condition for LES modelling. Application and limitation of the DNS model are described. Different experimental techniques for model validation are given. The consistency between physical experimentation/modelling and industrial cases is discussed. An emphasis is placed on the model validation through physical experimentation. Subsequently, the application of a turbulence model for three specific flow problems commonly encountered in metallurgical process, i.e., bubble-induced turbulence, supersonic jet transport, and electromagnetic suppression of turbulence, is discussed. Some future perspectives for the simulation of turbulent flow are formulated.

**Keywords:** turbulence modelling; RANS/LES/DNS; inflow condition; model validation; model application

#### **1. Introduction**

Since turbulence plays a significant role in flow transport phenomena, considerable efforts have been devoted to understanding flow turbulence in both engineering and academic communities [1–4]. Nevertheless, turbulence is still not fully understood due to its complexities in nature, e.g., instantaneous and intermittency characteristics, strong nonlinearity, and a wide range of temporal and spatial scales [4,5].

Turbulent flow can be numerically resolved with different levels of accuracy. Many numerical approaches for solving turbulence have been proposed, such as the Reynolds-Averaged Navier–Stokes (RANS) [6–11], the Large Eddy Simulation (LES) [12–14], and Direct Numerical Simulation (DNS) approaches [15]. Among these numerical methods, the RANS approach, specifically the Eddy Viscosity Model (EVM), is widely used for calculating turbulent flows thanks to its relatively high accuracy in predicting the mean flow features and its more limited computational demands. However, this approach suffers from several weaknesses, e.g., compromised accuracy and uncertainties due to assumptions in the model construction and insufficient incorporation of the fluid physics. In the LES approach, the whole eddy range is separated into two parts, namely, the large-scale eddy and subgrid-scale (SGS) eddy. The former can be directly resolved, while the latter is computed using the SGS model. As the computing power rapidly increases, this approach is extensively used to study turbulence physics and to resolve low-to-medium Reynolds number flows. In order to obtain meaningful results, additional attention is paid to inflow boundary conditions and mesh density for LES modelling. Compared to the RANS and LES, the DNS approach is the most accurate numerical method because it directly resolves

**Citation:** Wang, Y.; Cao, L.; Cheng, Z.; Blanpain, B.; Guo, M. Mathematical Methodology and Metallurgical Application of Turbulence Modelling: A Review. *Metals* **2021**, *11*, 1297. https:// doi.org/10.3390/met11081297

Academic Editors: Thomas Echterhof, Ko-Ichiro Ohno and Ville-Valtteri Visuri

Received: 10 June 2021 Accepted: 21 July 2021 Published: 17 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

all the turbulent eddies without using any models. However, DNS requires extensive computational demands to solve turbulent engineering flows, which is difficult to be satisfied by the current computing power. In addition, the large amount of data generated by DNS should be carefully analyzed.

Due to the uncertainties of the RANS method, many model variants have been proposed to solve specific turbulent flows [6–10], which can easily confuse the users of Computational Fluids Dynamics (CFD) codes to choose an appropriate one for their own cases. A comprehensive and well-organized description of the model formulation and development is very helpful in judging the applicability of different models, although there are some works with respect to model description [1,3]. Apart from SGS models of the LES method, the inflow boundary condition is an indispensable part of the simulation. Considering a diversity of methods generating inflow information [16,17], an elaborate analysis and discussion on the methods is necessary and conducive for the LES model user. In order to clearly understand the potential of the DNS method, the capabilities and current limitations of DNS method need to be clarified. In addition to the numerical solutions for turbulent flows, model validation is needed to warrant the accuracy of the simulation. Depending on different studies, different physical experimentations assisted with measuring techniques are performed to validate the numerical model. There are a number of studies focusing on the fundamentals of different measuring techniques, data interpretation, and applications of the techniques [18,19]. However, the consistency between physical experimentation/modelling and industrial applications needs to be discussed. The limitations and development of the measuring techniques are critical for obtaining reliable data and should be reviewed.

In this work, we review the formulation and development of three main numerical approaches (i.e., the RANS, LES, and DNS) for turbulence modelling. The advantages and disadvantages of the approaches are systematically discussed. Different methods of the inflow boundary condition are described for LES modelling. For the purpose of validating numerical models, different physical experimentation methods are presented. The consistency between the physical modelling and industrial applications is discussed. Limitations and progress of the experimental validation techniques are shown. Three turbulencerelated flow problems commonly encountered in metallurgical fields (i.e., bubble-induced turbulence (BIT), supersonic jet transport, and electromagnetic damping of turbulence) are discussed to demonstrate how to customize a conventional turbulence model for solving a specific flow problem. Finally, perspectives for modelling turbulent flows are proposed. With this review, we intend to help the current and potential CFD users to understand the modelling techniques for turbulence flows better and to expand the insight into the physics of turbulence.

#### **2. Turbulence Simulation Approaches**

#### *2.1. RANS*

In the RANS approach, instantaneous solution variables in the governing equations are decomposed into the mean and fluctuating components, as expressed in Equation (1).

$$f = f + f'\tag{1}$$

Substituting this variable expression into the instantaneous continuity and Navier-Stokes (N-S) equations yields the ensemble- or time-averaged forms for single-phase Newtonian flow, as shown in Equations (2) and (3). Henceforth, repeated-suffix summation convention is used in the formulae.

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x\_i} (\rho \overline{u}\_i) = 0 \tag{2}$$

$$\frac{\partial}{\partial t}(\rho \overline{u}\_i) + \frac{\partial}{\partial \mathbf{x}\_j}(\rho \overline{u}\_i \overline{u}\_j) = -\frac{\partial \overline{p}}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left(\mu \frac{\partial \overline{u}\_i}{\partial \mathbf{x}\_j} - \rho \overline{u\_i' u\_j'}\right) \tag{3}$$

The fluctuating quantities are included in the Reynolds stress tensor (−*ρu*- *i u*- *j* ) with six components. In order to close the equation set, the Reynolds stress tensor needs to be appropriately solved. One of the solutions for this closure problem employs the Boussinesq assumption [20], which relates the Reynolds stresses to the mean velocity gradients. The advantage of the Boussinesq assumption is the relatively low computational cost due to its simplicity. This works well for the engineering flows, which are dominated by only one turbulent shear stress such as the jet flow, wall boundary layer flow, and mixing layers flow. However, this approach is insensitive to the streamline curvature, rotation and body forces, and it exhibits a poor performance in the flows with a strong anisotropy or stress transport effect [21,22]. It also has difficulty in predicting transitional flows.

It is worth noting that the turbulent viscosity used in solving the Reynolds stress terms is a function of the space and flow features, rather than a physical parameter such as the fluid viscosity, which is dependent on the molecular structure of the fluid. Obviously, the turbulent viscosity needs to be solved before computing the Reynolds stress terms. In this paper, two-equation models, which include two additional transport equations, are reviewed. Usually the turbulence kinetic energy (*k*) is adopted as one equation and the turbulent kinetic energy dissipation rate (*ε*) or the specific dissipation rate (*ω*) as another one. The modified versions of the *ε*/*ω* model will also be presented here. Due to length restrictions, the zero- and one-equation models are not included in the article, but they may be found elsewhere [23–25].

#### 2.1.1. The *k*-*ε* Model

The standard *k*-*ε* (SKE) model was originally proposed by Launder and Spalding [6]. The model has been widely applied for resolving turbulent flows without a severe pressure gradient or strong swirling effect (e.g., plane jet, mixing layer, and boundary layer flows) because of its relatively high robustness, low computational cost, and reasonable accuracy. Equations (4) and (5) show the general form of *k* and *ε*. By solving these two transport equations, the turbulent viscosity can be calculated as expressed in Equation (6).

$$\frac{\partial}{\partial t}(\rho k) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho k \overline{u}\_l) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_l}{\sigma\_k} \right) \frac{\partial k}{\partial \mathbf{x}\_j} \right] - \rho \overline{u\_i' u\_j'} \frac{\partial \overline{u}\_l}{\partial \mathbf{x}\_j} - \rho \varepsilon + S\_{\mathbf{k}} \tag{4}$$

$$\frac{\partial}{\partial t}(\rho \varepsilon) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho \varepsilon \overline{\boldsymbol{u}}\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_l}{\sigma\_\varepsilon} \right) \frac{\partial \varepsilon}{\partial \mathbf{x}\_j} \right] - \mathbf{C}\_{1\varepsilon} \frac{\varepsilon}{k} \rho \overline{\boldsymbol{u}\_i' \boldsymbol{u}\_j'} \frac{\partial \overline{\boldsymbol{u}\_l}}{\partial \mathbf{x}\_j} - \mathbf{C}\_{2\varepsilon} \rho \frac{\varepsilon^2}{k} + \mathbf{S}\_\varepsilon \tag{5}$$

$$
\mu\_t = \mathbb{C}\_{\mu} \rho \frac{k^2}{\varepsilon} \tag{6}
$$

The terms from left to right in Equations (6) and (7) are the respective local time derivative, convection, diffusion, production, sink, and source terms. The SKE model is derived from a fully turbulent flow with high Reynolds numbers. The viscous effect is ignored in the model. However, this cannot be applied in the vicinity of the wall, where the viscous force dominates the flow characteristics. In order to deal with this problem, either a wall function is adopted with the SKE model, or a low Reynolds number model is used. The former confuses the users' judgement whether the weakness of this method lies in the basic SKE model itself or in the wall function. The latter requires additional functions to modify the standard transport equations. With respect to the low Reynolds number models [26–32], an example proposed by Lam and Bremhorst (LB model) [31] is presented in Equations (7) and (8).

$$\frac{\partial}{\partial t}(\rho k) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho k \overline{\mathbf{u}}\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_k} \right) \frac{\partial k}{\partial \mathbf{x}\_j} \right] + \mu\_t \left( \frac{\partial \overline{\mathbf{u}}\_i}{\partial \mathbf{x}\_j} + \frac{\partial \overline{\mathbf{u}}\_j}{\partial \mathbf{x}\_i} \right) \frac{\partial \overline{\mathbf{u}}\_i}{\partial \mathbf{x}\_j} - \rho \varepsilon \tag{7}$$

$$\frac{\partial}{\partial t}(\rho \varepsilon) + \frac{\partial}{\partial x\_i}(\rho \varepsilon \overline{u}\_i) = \frac{\partial}{\partial x\_j} \left[ \left( \mu + \frac{\mu\_l}{\sigma\_\varepsilon} \right) \frac{\partial \varepsilon}{\partial x\_j} \right] + \mathbb{C}\_{1\varepsilon} f\_1 \frac{\varepsilon}{k} \mu\_l \left( \frac{\partial \overline{u}\_i}{\partial x\_j} + \frac{\partial \overline{u}\_j}{\partial x\_i} \right) \frac{\partial \overline{u}\_i}{\partial x\_j} - \mathbb{C}\_{2\varepsilon} f\_2 \rho \frac{\varepsilon^2}{k} \tag{8}$$

where:

$$
\mu\_t = \mathbb{C}\_{\mu} f\_{\mu} \rho \frac{k^2}{\varepsilon} \tag{9}
$$

Compared with the SKE model, different formulations of functions *f*1, *f*2 and *fμ* are developed in the LB model to describe the near-wall behavior better. It has been confirmed that the function *fμ* has a predominant influence on the model performance, and functions *f*1, *f*2 play a secondary role in the performance [11].

There are also other modified versions of the SKE model, amongst which the widely used Renormalization Group (RNG) *k-ε* model [7] and the Realizable *k*-*ε* (RKE) model [8] are introduced in the following section. The main differences between the RNG *k-ε* model and the SKE model are the modifications of the turbulent viscosity and *ε* sink term. A differential equation is analytically derived for effective viscosity *μeff* to account for the low Reynolds number effect. This feature can improve the predictive ability of the RNG *k-ε* model for low Reynolds number flows or near-wall flows. Additionally, a new *ε* destruction term is used to account for the rapid strain by modifying the constant of this term. The RKE model is modified mainly with regard to the turbulent viscosity and the *ε* equation. By defining a variable *Cμ* [8,33,34] instead of a constant value in turbulent viscosity formulation, the RKE model satisfies the realizability constraints, i.e., positive values for the normal stresses and the Schwartz inequality for the shear stresses. In order to increase the robustness of the model, a new *ε* equation is employed based on a dynamic equation for fluctuating vorticity. The new *ε* equation describes turbulent vortex stretching and turbulent dissipation more appropriately compared to the *ε* equation in the SKE model. With the modified *ε* equation, the well-known round-jet anomaly that is a poor prediction of the spreading rate of round/axisymmetric jet may be solved [8].

#### 2.1.2. The *k*-*ω* Model

The *k*-*ω* model is widely used for turbulence modelling [35]. Different versions of this model have been developed in the last decades [9,10,36–40]. In this paper, the most well-known *k*-*ω* model proposed by Wilcox [9] is reviewed. The transport equations of this model are presented in Equations (10) and (11), where the calculation of *β*1 refers to the work of Wilcox [9].

$$\frac{\partial}{\partial t}(\rho k) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho k \overline{\mathbf{u}\_i}) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_l}{\sigma\_k} \right) \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_j} \right] - \rho \overline{\mathbf{u}\_i^{\prime} \mathbf{u}\_j^{\prime}} \frac{\partial \overline{\mathbf{u}\_i}}{\partial \mathbf{x}\_j} - \rho \boldsymbol{\beta}^\* k \boldsymbol{\omega} + \mathbf{S}\_k \tag{10}$$

$$\frac{\partial}{\partial t}(\rho\omega) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho\omega\overline{\mathbf{u}\_i}) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_l}{\sigma\_\omega} \right) \frac{\partial \omega}{\partial \mathbf{x}\_j} \right] - \alpha \frac{\omega}{k} \rho \overline{\mathbf{u}\_i^\prime \mathbf{u}\_j^\prime} \frac{\partial \overline{\mathbf{u}\_i}}{\partial \mathbf{x}\_j} - \rho \beta\_1 \omega^2 + S\_\omega \tag{11}$$

Unlike the *k-ε* model, the *k*-*ω* model can be integrated through the viscous sublayer without any damping function to account for the low Reynolds number effect with high numerical stability. Therefore, it is well applied in aerodynamic flows [10,35]. However, the *k-ω* model is highly sensitive to the empirical value of *ω* at the free edge of the turbulent shear layer, which can lead to a large prediction error. In order to solve this problem, a modified model was proposed with combination of the original *k*-*ω* model and the SKE model by adding a blending function [10]. This new model is termed the Baseline *k*-*ω* model, which applies the original *k*-*ω* model in the near-wall region and switches to the SKE model in the outer region. The Baseline *k*-*ω* model has a similar performance to the original *k*-*ω* model in boundary layer flows, but the former one avoids the strong freestream dependence. However, both *k-ω* models fail to predict the onset and amount of separation in adverse pressure gradient flows. Based on the Baseline *k*-*ω* model, further modification to eddy viscosity is proposed to account for the transport effects of the principal turbulent shear stress, leading to a significant improvement in predicting the adverse pressure gradient flows [10]. However, the introduced blending function depends on empiricism (e.g., the distance to wall), limiting its application to flows in complex geometries.

2.1.3. Advanced Eddy Viscosity Model

All the aforementioned methods belong to the category of the Eddy Viscosity Model (EVM). Some other advanced EVMs were developed [41–45]. In order to account for the strong anisotropy in the near-wall region, Durbin [41] adopted a new turbulent viscosity term defined in Equation (12), which is considered to be more appropriate than that defined in Equation (6) in a near-wall region. A separate transport equation for a wallnormal turbulent stress *υ*2 was proposed and solved with the aid of the elliptic relaxation concept. This model is termed the *υ*2-*f* model, where *f* represents the elliptic relaxation function. Subsequently, several modified versions were proposed with respect to the velocity scale [44], the characteristic length [42], the function *f* [45], and the variable *υ*2 [43]. This model category performs well for pressure-induced separating flow, buoyancy impairing turbulent flow, and backstep flow [44–46].

$$
\mu\_{\ell} = \mathbb{C}\_{\mu} \rho \overline{\upsilon^2} \tau \tag{12}
$$

where:

$$\pi = \max\left(\frac{k}{\varepsilon}, \mathbb{C}\_{\tau} \left(\frac{\mu}{\rho \varepsilon}\right)^{0.5}\right) \tag{13}$$

In addition to the modified versions of the linear EVM, the idea of non-linear EVM has been substantially used [47–53]. Even though these modified models demonstrated certain improvements over linear EVMs in predicting flows with a strong streamline curvature or turbulent stresses in the near-wall sublayer, they are still inferior to the more advanced model, e.g., the RSM model as seen in Section 2.1.4.

#### 2.1.4. Reynolds Stress Model

In order to overcome the limitations of the EVM, Second-Moment Closure (SMC) models abandoning the Boussinesq assumption have been developed. The SMC model directly solves the transport equation for each of the Reynolds stress terms. Since the SMC approach considers the effects of streamline curvature, rotation, and rapid change of strain rate in a more rigorous manner, it is long expected to replace the currently widely applied two-equation models. The SMC model class consists of the Algebraic Stress Model (ASM) and the Differential Stress Model (DSM). ASM is derived from differential stress transport equations by invoking the weak-equilibrium assumption [54–56]. It ignores the transport terms of the anisotropy by assuming that the transport of the Reynolds stress is proportional to that of turbulent kinetic energy [57,58]. In general, the ASM is considered an intermediate tool between the LEVM and the DSM. Due to the space limitation, only the DSM is presented in this paper. A symbolic representation of the stress transport equation is expressed in Equation (14). In addition, a scale-determining equation, i.e., the *ε* equation, is needed to complete the DSM.

$$L\_{i\bar{j}} + C\_{i\bar{j}} = D\_{T,\bar{i}\bar{j}} + D\_{L,\bar{i}\bar{j}} + P\_{i\bar{j}} + G\_{i\bar{j}} + \phi\_{i\bar{j}} + \varepsilon\_{i\bar{j}} + s\_{i\bar{j}} \tag{14}$$

The terms from left to right represent the local time derivative of Reynolds stress, convection, turbulent diffusion, molecular diffusion, stress production, buoyancy production, pressure strain, dissipation and source, respectively. The *sij* term is user-defined for a specific stress transport source. If there is no source, this term becomes zero. It is required to model *DT*,*ij*, *Gij*, *φij* and *εij* to close the equation, while it is not necessary to model *Lij*, *Cij*, *DL*,*ij* and *Pij*.

The *DT*,*ij* term includes the velocity transport and the pressure transport. The velocity triple moments can be measured, whereas the pressure transport is intractable. Usually, the pressure transport is considered to be negligible [59]. Therefore, the model is mainly designed for the velocity triple moments. The most popular model is the generalized gradient-diffusion model proposed by Daly and Harlow (DH) [60]. The DH model has a symmetry problem in the indices, leading to dependence on the coordinate frame. Subsequently, some variants of this model were developed by Hanjalic and Launder [61], Shir [62], Mellor and Herring [63]. More complex models were also put forward by Nagano and Tagawa [64] and Magnaudet [65]. Due to uncertainties in the modelling equations, the complex models may not necessarily outperform the simplified models. Given more computing resources consumed by the complex models and a poor convergence, application of these models in engineering has been doubted [3,66].

Compared to the *εij* and *Gij* terms, it is necessary to pay extra attention to model *φij*. Usually, the *φij* term is decomposed into three parts, namely, the slow pressurestrain term *φij*,1, rapid pressure-strain term *φij*,2, and wall-reflection term *φij*,*w*. Not all of the models introduced below include the third term. Rotta [67] proposed a linear model for *φij*,1, which considers that *φij*,1 is proportional to the stress anisotropy tensor. However, this linear model is unable to satisfy the realizability constraints. A general quadratic model [68] is proposed to solve this problem. Linear [61,69–71] and nonlinear models [72–74] were proposed to model the *φij*,2 term. Even though the nonlinear model is considered to be theoretically advanced, the complexity of the formulation prohibits its application for engineering computation. The turbulence anisotropy is enhanced due to the damping effect of the normal stress by the wall. This damping affects both the pressure-strain terms. In order to account for the damping effect, a commonly used model [62,75] is presented to model the wall-reflection term *φij*,*w*. However, this model involves a variable, i.e., the normal distance to the wall. This is believed to be a major weakness. For purpose of overcoming this weakness, the elliptic relaxation concept and elliptic-blending method were proposed to account for the near-wall inhomogeneity, which is described in Section 2.1.3, and more information on that can be found in [76,77].

The DSM is the most elaborate model in the RANS approach, which has an indisputable superiority over the rudimentary two-equation models in predicting complex flows, e.g., highly swirling and rotating flow, separating flow, and secondary flow. However, its application is limited by (1) a high degree of uncertainty in modelling the high-order correlation terms (e.g., pressure-strain and dissipative correlation) due to an insufficient knowledge of physics; (2) a high demand for computational resource (approximately 50–150% more computing time than a two-equation EVM [21]). Fortunately, due to the use of more advanced models (e.g., the DNS), the accuracy and robustness of the DSM have been improved. A rapid development in computer science (e.g., parallel processing and improved performance) satisfies the high computational need for the use of DSM model. The DSM has received more attention recently because of the failure of the EVM in predicting complex turbulent flows.

#### *2.2. LES Approach*

#### 2.2.1. Formulation and Subgrid-Scale Model

Turbulent flow features a wide range of eddy scales from the Kolmogorov length scale to the size comparable to the characteristic length of the mean flow. The large eddies contain most of the turbulent energy and are mainly responsible for the momentum and energy transfer. They are strongly affected by boundary conditions. The small eddies tend to be more isotropic and homogeneous, and their dissipation process is linked to fluid viscosity. For this reason, it is very difficult for the RANS approach to model all the eddies in a single model. The LES approach separates the large eddies from the small ones by employing a spatial filtering method [78] for the instantaneous governing equations. After that, the large eddies are directly resolved by the filtered equation, and the small ones (i.e., the Subgrid-Scale (SGS) eddies) are modelled by the SGS model. The filtered variable (donated by an overbar) is defined by Equation (15). The resulting continuity and momentum equations are expressed in Equations (16) and (17), showing similar forms but a different physical meaning for those in the RANS approach.

$$\overline{f}(\mathbf{x}) = \int\_{D} f\left(\mathbf{x'}\right) G\left(\mathbf{x}, \mathbf{x'}; \overline{\boldsymbol{\Delta}}\right) d\mathbf{x'} \tag{15}$$

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x\_j} (\rho \overline{u}\_i) = 0 \tag{16}$$

$$\frac{\partial}{\partial t}(\rho \overline{u}\_i) + \frac{\partial}{\partial \mathbf{x}\_j}(\rho \overline{u}\_i \overline{u}\_j) = -\frac{\partial \overline{p}}{\partial \mathbf{x}\_i} + \frac{\partial \sigma\_{ij}}{\partial \mathbf{x}\_j} - \frac{\partial \mathbf{r}\_{ij}}{\partial \mathbf{x}\_j} \tag{17}$$

where *σij* refers to the stress tensor due to molecular viscosity, and *τij* represents the SGS stress term. Most of current SGS models adopt the Boussinesq assumption (termed the Eddy-Viscosity model), which relates the SGS stress to the large-scale strain-rate tensor, as shown in Equation (18). Based on the definition of the eddy viscosity, various SGS models [12–14,79,80] have been proposed. Smagorinsky [12] developed the first SGS model by assuming a local energy equilibrium between the large scale and the subgrid scale. The eddy viscosity in this model is defined in Equation (20). This model becomes very popular to date due to its simplicity, numerical robustness, and stability. However, it has several drawbacks: (1) the model constant varies with different flows; (2) the model cannot predict the inverse energy transfer (i.e., backscatter) due to its purely dissipative nature; (3) the model has difficulty in reproducing the correct mean quantities (e.g., SGS dissipation) as the grid scale approaches the integral scale; (4) the model does not yield a zero eddy viscosity in near-wall regions.

$$
\pi\_{ij} - \frac{1}{3}\pi\_{kk}\delta\_{ij} = -2\mu\_t\overline{S}\_{ij} \tag{18}
$$

where:

$$\overline{S}\_{ij} = \frac{1}{2} \left( \frac{\partial \overline{u}\_i}{\partial x\_j} + \frac{\partial \overline{u}\_j}{\partial x\_i} \right) \tag{19}$$

$$
\mu\_t = \rho \left(\mathbb{C}\_s \Lambda\right)^2 \left|\overline{S}\right|\tag{20}
$$

In Equation (20), where:

$$|\overline{\mathcal{S}}| = \sqrt{2\overline{\mathcal{S}}\_{ij}\overline{\mathcal{S}}\_{ij}}\tag{21}$$

$$
\Delta = \left(\Delta x \Delta y \Delta z\right)^{1/3} \tag{22}
$$

where Δ represents the local grid scale. In order to solve the model constant problem, a Dynamic Smagorinsky–Lilly (DSL) model was proposed [13,14], where the model constant is dynamically calculated by using the resolved eddies with the scale size between the grid filter and test filter. The main advantage of this model is that it is not necessary to prescribe and/or tune the model constant. However, the DSL model is subjected to a numerical instability and a variable model constant. Germano [13] proposed an averaging method to overcome this weakness. A good performance was achieved in a channel flow simulation [81]. Another variant of Smagorinsky–Lilly model is the Dynamic Kinetic Energy (DKE) model [82–84]. Unlike the algebraic form in Smagorinsky–Lilly and DSL models, the DKE model solves an additional transport equation for the SGS turbulent kinetic energy instead of adopting the local equilibrium assumption. This model can better account for the energy transfer from the large-scale eddy at the cost of computational expenses. Some other variants of the Smagorinsky–Lilly model were formulated to solve the low Reynolds number effect in the near-wall region, one of which was based on the square of the velocity gradient tensor named the Wall-Adapting Local Eddy-viscosity (WALE) model [79]. Compared to the original Smagorinsky–Lilly model, the WALE model can produce a zero eddy viscosity in the vicinity of the wall or in a pure shear flow. Hence, this model does not need a damping function. In addition to the WALE model, a hybrid model [85] was proposed by combining the Smagorinsky–Lilly model with a damping function [86] to improve the predictive capability for wall-bounded flows. This hybrid model demonstrated a good performance in plane channel flow with different Reynolds numbers. However, this model involves a variable, i.e., the wall normal

distance. The determination of this wall normal distance requires an empirical approach for specific flow.

In the framework of the eddy-viscosity SGS model, there are several alternatives to the Smagorinsky-type SGS model, such as Vreman's model [87], the QR model [88,89], the σ-model [90], and the S3PQR model [91]. Compared to the Smagorinsky-type SGS model, Vreman's model can predict zero eddy viscosity in near-wall regions or in transitional flows without explicit filtering, averaging or clipping procedures. However, it was found that the model coefficient in Vreman's model is far from universal. To solve this problem, two procedures were proposed to dynamically determine the model coefficient, i.e., the one based on the global equilibrium between the subgrid-scale dissipation and the viscous dissipation [92,93] and the other one based on the Germano identity [94]. It was reported that the latter is better suited for transient flows [94]. The QR model, which is a minimum-dissipation eddy-viscosity model, gives the minimum eddy dissipation required to dissipate the energy of sub-filter scales. The advantages of this model lie in appropriately switching off for laminar and transitional flows, the low computational complexity, and consistency with the exact sub-filter tensor on isotropic grids. The disadvantage of this model is the insufficient eddy dissipation, which can be corrected by increasing the model constant. Moreover, the QR model requires an approximation of the filter width to be consistent with the exact sub-filter tensor on anisotropic grids. It was noted that the accuracy of the model result for anisotropic grids is highly dependent on the used filter width approximation. By modifying the Poincaré inequality used in the QR model, the dependence can be removed, leading to the construction of an anisotropic minimum-dissipation model that generalizes the desirable properties of the QR model to anisotropic grids [95]. For the purpose of meeting a set of properties based on the practical/physical considerations, the σ-model based on the singular values of the velocity gradient tensor was developed [90]. Owing to its unique properties, ease of implementation, and low computational cost, the σmodel is considered to be suitable for complex flow configurations. Subsequently, through comparison between static and dynamic σ-models, it was found that the local dynamic procedure is not suited for the σ-model, and a global dynamic procedure is suggested [96]. Trias et al. [91] built a general framework for LES eddy-viscosity models, which is based on the 5D phase space of invariants. By imposing appropriate restrictions in this space, a new eddy-viscosity model, i.e., the S3PQR model, was developed. In addition to meeting a set of desirable properties such as positiveness, locality, Galilean invariance, proper near-wall behavior, and automatic switch-off for laminar, 2D and axisymmetric flows, this new model is well-conditioned and has a low computational cost, with no intrinsic limitations for statistically inhomogeneous flows. Despite of all the merits of this model, special attention should be given to the calculation of the characteristic length scale and the determination of the model constant before engineering applications.

Alternatives to the eddy-viscosity SGS model are the similarity model [97–100], the velocity estimation model [101,102], the Approximate Deconvolution model (ADM) [103–105], and the regularization model [106–108]. The similarity model class adopts the idea that an accurate approximation for a SGS model can be reconstructed from the information contained in the resolved field. Therefore, the similarity models [97–100] approximate the SGS stress tensor by a stress tensor calculated from the resolved scales. Due to this nature, the similarity model can naturally account for the inverse energy transfer (i.e., backscatter). This is different from the eddy viscosity model, which only considers the global SGS dissipation (i.e., the net energy flux from the resolved scales to the subgrid scales). It is worth mentioning that, due to the importance of accurate prediction of the inverse energy transfer, a dynamic two-component SGS model was proposed to include the non-local and local interactions between the resolved scales and subgrid scales [100]. The model correctly predicted the breakdown of the net transfer into forward and inverse contributions in a priori tests. In some cases, however, the similarity model underestimated the SGS dissipation. An extra dissipative term was added to solve this issue. This new model formulated is also referred to as the mixed model. Furthermore, the similarity model

and the mixed model need additional computational resources due to the implementation of the second filtering. Special attention should be paid to choose an appropriate filter type and size. Following the same idea, Domaradzki et al. [101] improved the SGS stress approximation by replacing the unknown unfiltered variables by approximately deconvolved filtered variables. Subsequently, this SGS model based on the estimation of the unfiltered velocity, which was originally formulated in spectral space, was extended to the physical space [102]. It was found that both versions of this velocity estimation model perform better than or are comparable to classical eddy viscosity models for most physical quantities. This model can account for backscatter without any adverse effects on the numerical stability. Several questions for improving the model need to be addressed, such as the modelling of nonequilibrium and high Reynolds number turbulence in three Cartesian directions. Stolz et al. [103–105] proposed a formulation of the ADM for LES, in which an approximation of the unfiltered solution is obtained from the filtered solution by a series expansion involving repeated filtering. Given a good approximation of the unfiltered solution at a time instant, the nonlinear flux terms of the filtered N-S equations can be computed directly, avoiding the explicit computation of the SGS closures. The effect of the non–represented scales is modelled by a relaxation regularization involving a second filtering and a dynamically estimated relaxation parameter. The ADM is evaluated for incompressible wall-bounded flow [104] and compressible flows [103,105]. The results showed that the ADM can have a significant improvement over the standard and dynamic Smagorinsky models, while at less computational cost compared with that of the dynamic models or the velocity estimation model. The high-Reynolds-number supersonic flow [109] and transitional flow [110] were investigated by the ADM. Agreement was observed between the ADM and experiments or DNS. For the former flow, a rescaling and recycling method was used to have a better control on the desired inflow data. Recently, the ADM was extended for a two-phase flow simulation [111]. By comparing the macroscopic flow characteristics, the ADM showed a better performance than the conventional LES model. However, further investigations should be performed on the relaxation term model for a two-phase simulation and microscopic characteristics of the dispersed phase in a 3D simulation. Another important class of the SGS model is the regularization model, which combines a regularization principle with an explicit filter and its inversion. The regularization model includes many versions, such as the Leray model [106], the Leray-α model [107,112], the Clark-α model [108], the Navier–Stokes-α (NS-α) model [113], etc. For the last one, several variants are proposed, including the NS-α deconvolution model [114–117] and the reduced order NS-α model [118,119]. It was found that the NS-α deconvolution model can significantly improve the prediction accuracy by carefully choosing the filtering radius and by correctly selecting the approximate deconvolution order [117]. Given the difficulties in efficient and stable simulation of the NS-α model for incompressible flows on coarse grids, the reduced order NS-α model is introduced by using deconvolution as an approximation to the filter inverse, reducing the fourth-order NS-α formulation to a second-order model. In spite of the success of the reduced order NS-α model, future work needs to be conducted on locally and dynamically choosing α and numerical testing on different benchmark flows, to name but a few. In addition, comparative studies have been performed between different regularization models [120–122], in which the capability of the regularization model has been demonstrated for a specific flow.

The LES is considered a compromise between the RANS and the DNS. It is more accurate than the RANS and it needs less computational resources than the DNS. However, the LES model has not reached the maturity stage as a numerical tool for the design or the parametric study of complex engineering flows, due to not only a high computational requirement, but also many unresolved issues such as ill-defined boundary conditions, wall-resolved flow, turbulent flow with chemical reactions, and compressible flow. Nevertheless, the LES model has been successfully applied in transitional flow [123–126], separated flow [127,128], and bubbly flow [129,130]. Figure 1 shows the calculation results in a separated boundary layer transition on a flat plate with a semi-circular leading edge of

radius of 5 mm under elevated free-stream turbulence. A periodic boundary condition was adopted in the spanwise direction. Free-slip and no-slip conditions were used at lateral boundaries and the plate surface, respectively. The simulation agrees well with experimental data on mean and fluctuating streamwise velocities for an Enhanced-Turbulence-Level (ETL) case, demonstrating a good performance of the LES model for the transitional flow. Figure 2 shows a comparison between a modified *k*-*ε* model and the LES model in a bubbly flow. Compared with the experimental data, the LES model is superior in predicting the turbulence.

**Figure 1.** *Cont*.

**Figure 1.** (**a**) Mean streamwise velocity, (**b**) rms streamwise velocity at different streamwise stations (ETL-case) (LES: solid lines, experimental data: symbols, x/l: normalized streamwise, y/l: wall-normal direction, U/U0: normalized mean streamwise velocity, u'/U0: normalized rms streamwise velocity) [126].

**Figure 2.** Comparison of model predictions with experimental data for (**a**) axial and (**b**) lateral liquid velocity fluctuations at a height of 0.25 m in a square cross-sectioned bubble column with a superficial gas velocity of 4.9 × 10−3 m/s [129].

#### 2.2.2. Inlet Boundary Condition

The fluid behavior in the domain is largely determined by the inflow condition [131]. The treatment of the inflow condition is of significant importance for LES modelling. Currently, there are two main categories for generating the inflow data, namely, precursor simulation and a synthetic method. The former involves a separate simulation, where the

periodic boundary condition or the recycling method can be used. The flow data are stored at each time step in this simulation and then introduced to the inlet boundary for modelling the flow of interest. The main advantage of this method is to obtain more realistic inflow conditions, which represent the required flow characteristics (e.g., velocity profile, turbulence intensity, shear stress, power spectrum, and turbulent structures). The inlet boundary, however, needs to be placed in an equilibrium region for scaling arguments in the precursor simulation, which may even not exist in some flows. This method may lead to a spurious periodicity for the time series [132]. In addition, running a separate simulation requires high computational costs especially for a high-Reynolds number flow. This restrains its application to complex engineering flows. The limitation can be reduced by an internal mapping method. This method integrates the precursor simulation into the main domain, mapping the data downstream of the inlet back to the inlet boundary [133,134].

The synthetic method as an alternative to the precursor simulation is expected to construct the inflow conditions for practical flows. The simplest way is to impose a white-noise random component on the inlet velocity. The magnitude of this random component is determined by the turbulent intensity. Since the turbulence-like component is rapidly dissipated due to the lack of spatial and temporal correlation, this white-noise method is inappropriate to generate the inflow data [135]. In order to impose realistic inflow data on the inlet boundary, other advanced synthetic methods have been developed. These advanced methods consist of the Fourier technique [136,137], principal orthogonal decomposition (POD) method [138,139], digital filtering technique [140,141], and synthetic eddy method (SEM) [142–144]. Several comparative studies have been performed on different synthetic methods [143–146]. Jarrin et al. [143] used the SEM in the hybrid RANS/LES simulations for turbulent flows from simple channel and square duct flows to the flow over an airfoil trailing edge. Compared to other synthetic methods (i.e., Batten's method (Fourier method) [136] and random method), the SEM can substantially reduce the inlet section, leading to a large decrease in the CPU time. Figure 3 shows a better performance of the SEM in the inlet velocity vector compared to that of the others. The SEM realistically reproduced the magnitude and length scale of the fluctuations. The fluctuations in the Batten's model decorrelated in space in the near-wall region due to the decomposition in Fourier's mode. Recently, Skillen et al. [144] improved the SEM of Jarrin et al. (Original SEM) [142] with respect to the normalization algorithm and the eddy placement. The former leads to an improvement over the original SEM model. The latter saves a cost of around 1–2 orders of magnitude. Figures 4 and 5 show the turbulent shear stresses and skin-friction coefficients from the original SEM [142], the improved SEM of Skillen et al. [144], and the precursor LES of Kaltenbach et al. [147]. In comparison with the original SEM, the improved SEM shows a better agreement with the precursor LES results both for the turbulent shear stress and the skin-friction coefficient. Although there is a rapid development of synthetic methods, the available synthetic methods are limited to construct the inflow condition with all required turbulence characteristics mentioned above. Further development is needed. It is unnecessary to claim which method is the best; however, the most appropriate method can be selected by considering the accuracy and the computational cost.

**Figure 3.** Velocity vectors of the LES inlet condition for hybrid simulation of channel flow: (**a**) precursor LES, (**b**) SEM, (**c**) Batten's method, (**d**) random method (y/*δ*: normalized y distance, z/*δ*: normalized z distance, *δ*: initial boundary layer thickness) [143].

**Figure 4.** Profiles of turbulent shear stress (the number at the top of the figure indicates profile locations (the distance from the start of the incline, normalized by the initial boundary layer thickness)) [144].

**Figure 5.** Wall shear stress for the inclined wall (*Cf*: wall shear stress, *x*/*δ*: normalized distance from the start of the incline) [144].

#### *2.3. DNS Approach*

The DNS model numerically solves the three-dimensional, time-dependent N-S equations without using a turbulence model. The DNS method captures all the turbulence scales present in the given flow by directly using a very fine mesh and very small time step. The application of the DNS approach is hindered by the requirement of extremely high computational resource. It has been estimated that the number of grid points is proportional to *Re*9/4 in a DNS case [148]. Since the eddy scale in the near-wall region is much smaller than that in the outer domain, a refinement of the mesh is needed in the near-wall region to fully resolve the turbulence, which further increases the number of grid points. Given the large computational domain, complex geometry, and high Reynolds number in practical engineering, the application of DNS approach is currently impossible for most practical flows. However, with the continuous development of parallel computing technique [149], hybrid CPU + GPU computing architecture [150,151], and advanced numerical algorithm [152], a remarkably high computing performance has been reached for DNS. Meanwhile, new challenges have arisen regarding the high-performance computing. Discussion on this point is excluded due to the limited space in this work; however, readers can be directed to relevant work for more information [153,154].

To demonstrate how DNS resolves flow turbulence, an example case is presented for bubble-induced turbulence (BIT), which is one of the important research topics in gas– liquid flow. In the study by Feng et al. [155], a precursor DNS simulation was performed on a homogeneous single-phase turbulent flow. The results of this precursor simulation were used as the inflow condition to calculate the turbulent field around a fully resolved bubble. Figure 6a,b shows the mesh used in the BIT study and the turbulent eddy generated on the highly deformable bubble surface. The results showed that the bubble created new vortices in the wake region, leading to turbulence enhancement. The magnitude of the turbulence enhancement increased with the liquid turbulent intensity and the relative velocity [155]. Apart from the direct study of turbulence, DNS plays an indispensable role in evaluating and developing turbulence models [156–161] and in providing complementary information for experimental study [162–164].

**Figure 6.** (**a**) Mesh used in the BIT study, (**b**) turbulent eddy generated on the highly deformable bubble surface (Q-criterion with a value of 100) [155].

#### **3. Model Validation Technique**

In order to guarantee the fidelity of the established model, model validation is a must before practical simulation. This can be done numerically or experimentally. In the numerical way, an advanced approach such as the DNS can be used to verify a less advanced model such as the RANS. The merit of this validation technique is to reproduce the real conditions of the practical flow, leading to more comprehensive and realistic results. Since the DNS is not fully ready for the simulation of complex engineering flows due to the computing power problem, this technique is less popular compared to the experimental technique. Table 1 lists several commonly used experimental techniques. Depending on the purpose, different experimental methods are employed for the validation, such as an optical fiber probe to measure solid phase holdup and solid velocity in a gas–solid system [165], LDA (Laser-Doppler Anemometry) and PIV (Particle Imaging Velocimetry) to measure phase velocity in a single- or two-phase system [166,167], an electroresistivity probe to monitor mixing behavior in metallurgical processes (converter and ladle) [168,169], and video recording to capture cavity shape and dimension (converter) [170].

**Table 1.** Experimental validation techniques.

Note: DIA: Digital Image Analysis.

Experimental validation is very useful for understanding model development and turbulence physics, while it also suffers from many limitations. Improvement is needed to represent more realistic situations. In the metallurgical field, it is very difficult to measure flow velocity and flow pattern due to the aggressive and complex operating environment (e.g., high temperature, opaque vessel, and multiphase coexistence). Physical modelling experimentation is adopted based on the similarity principle to obtain flow information. This validation method, however, has several shortcomings: (1) A full-scale physical model is prohibitive to set up in laboratory study due to the difficulty in building a large industry-scale vessel, and the difficulty in mimicking the industrial operational conditions. In general, a scaled-down low temperature physical model is usually employed. This can only reproduce part of the flow dynamics since not all of necessary flow dimensionless numbers can be simultaneously satisfied [202]; (2) Most physical modelling experiments are conducted at room temperature, at which it is impossible to study the heat transfer and melt solidification behavior. In addition, the effect of temperature on the gas phase is not considered in such experiments; (3) The reliability of physical modelling depends on a selection of materials used in the simulation of the real system, resulting in an additional experimental error. For instance, water is frequently used to mimic liquid steel because the kinematic viscosity of water is very close to that of the liquid steel at 1600 ◦C. Other properties of water, however, such as the density and surface tension, are very different from liquid steel, making the similarity criteria difficult to be entirely fulfilled. Compared to water, low melting point alloys, such as Bi-Sn and Ga-In-Sn alloys, have a closer resemblance to the physical properties of liquid steel [203,204]. The Bi-Sn alloy system can be operated in a temperature range of 200–400 ◦C, compared to the GaInSn alloy system at room temperature. The effect of temperature on liquid viscosity can be investigated. Due to similar electrical conductivity between the low-melting point alloy and liquid steel, the alloy can be used to model the electromagnetic stirring or breaking in the continuous casting process. With the aid of measuring techniques such as Ultrasound Doppler Velocimetry (UDV) [205] and Contactless Inductive Flow Tomography

(CIFT) [206,207], a significant progress has been achieved with this physical modelling for simulating the real flow behavior in the continuous casting process.

There are other errors induced by the measurement techniques. For PIV and LDA measurement, a correct phase discrimination can reduce errors, improving the measurement quality. Different phase discrimination methods have been proposed. Kulick et al. [208] exploited the large difference in the amplitude of the Doppler burst pedestals obtained from the solid particle and the tracer and then took the ratio of the Doppler signal amplitude to the pedestal amplitude as the discriminator. Regarding the PIV technique, Kiger and Pan employed a two-dimensional median filter to correctly identify and separate the dispersed particles from the two-phase image [209]. Khalitov and Longmire [210] adopted a two-parameter (size and brightness) algorithm to separate the tracer from the solid phase. This algorithm has been proved to be applicable in gas–solid and liquid–liquid systems [210,211]. In addition, the Turbulence Kinetic Energy (TKE) dissipation rate is very difficult to be measured by the PIV or LDA technique since it strongly depends on the spatial resolution [212]. Tanaka and Eaton [213] performed High-Resolution PIV (HR-PIV) measurements with a spatial resolution smaller than the Kolmogorov scale *η*k. A modified method was used for the phase identification algorithm to eliminate the common error in the HR-PIV measurement. It was concluded that the measurement error of the TKE dissipation rate can be reduced to a few percent if a proper spatial resolution is employed (in the range of *η*k/10 to *η*k/2). The HR-PIV is also used to measure the carrier-phase velocity and turbulence structures near the particle surface [214,215]. The attenuations of the TKE and its dissipation rate were experimentally obtained, leading to a reasonable prediction of the macroscopic turbulence modification. Even though the HR-PIV technique works well for diluted flow, further development of the HR-PIV technique is needed for investigating the dense flow or the flow region with a high volume fraction of the dispersed phase. By using the HR-PIV technique, Wang et al. [130,216] found a large deviation of the velocity magnitude in the bubbly zone of a top-submerged gas injection flow. This deviation was partly attributed to the less reliable HR-PIV experimental data caused by the refraction of gas bubbles. For brevity's sake, other techniques, for instance, the measurement of the pressure drop and fluctuation and holdup of solid phase are excluded from this work. Interested readers can find relevant studies elsewhere [217,218].

#### **4. Applications**

Most of the flows involved in industrial production are turbulent. Therefore, the choice of a proper turbulence model used to accurately represent the flow characteristics is essential to CFD applications. In order to draw attention to this point, the applications of turbulence models in solving three commonly encountered turbulence-related flow issues are discussed. It is necessary to mention that this work aims to demonstrate how to customize a conventional turbulence model for a specific flow problem instead of comprehensively reviewing the single phase or multiphase phase flow characteristics.

#### *4.1. Bubble-Induced Turbulence (BIT)*

In an engineering bubbly flow such as the flow in a bubble column or the flow in a metallurgical ladle, the BIT has to be taken into consideration to fully reproduce the flow characteristics. There are two common methods to account for the BIT: (1) the Effective Viscosity Modified Method (EVMM), where BIT is added in the effective viscosity term [191,219]; (2) the Transport Equation Modified Method (TEMM), where BIT is modelled with a source term in the turbulent transport equation [220–225]. Both methods can be used with the RANS-based approach (*k*-*ε*/*ω*/SST (Shear-Stress Transport) [224] and RSM, [226–228]) termed RANS-BIT, and LES approach, termed LES-BIT [129,191,229]. Many RANS-BIT comparative studies were performed on the EVMM and TEMM [224,225]. The results showed that the EVMM is not suitable for the turbulence prediction since it depends on an algebraic model (e.g., the model of Sato et al. [219]). LES-BIT studies were performed with the EVMM, and it was observed that including the BIT slightly changed the

numerical results [129,191]. Since different experimental validations suggest different formulations of the source term, no universal formulation of the source term has been found. Compared to the *k*-*ε*/*ω*/SST-based BIT modelling, the RSM-based BIT modelling [225,228] considers the effect of the anisotropic property of turbulence. This improves the overall model performance in simulating both the void fraction distribution and turbulent kinetic energy. This model is more suitable for modelling the complex multiphase flow with prominent anisotropy. Nevertheless, it introduces more unknown parameters, making this model less readily applicable. Reliable experiments with individual Reynolds stress measurement are needed for the model validation. Niceno et al. [229] compared the EVMM and TEMM of the LES-based approach and concluded that both methods predict the turbulent kinetic energy well qualitatively, but the TEMM is superior for the quantitative prediction.

Recently, the DNS study became available for bubbly flows; several DNS studies [161,230,231] on disperse bubbly channel flows were performed, and the obtained TKE budgets were used to account for the BIT. With the aid of the DNS data, Ma et al. developed a model for the BIT and incorporated it as a source/sink term in the *k*-*ω* SST model transport equation. This model was adopted by Liao et al. to evaluate its performance on bubbly flows in containers and vertical pipes [232]. The results showed improvements regarding the radial gas volume fraction and velocity profile in high-volume fraction cases were achieved. Later, Ma et al. [231] extended the BIT model to a second moment level. During the development of this full SMC for BIT in the Euler-Euler framework, particular attention was given to the treatment of the pressure-strain term for bubbly flows and the form of the interfacial term to account for BIT. For the latter, an effective BIT source term was proposed, which largely simplified the modelling work. To understand the anisotropic behavior of the bubbly flow, an anisotropy-invariant analysis was conducted, based on which the BIT closure was improved. This new SMC with the proposed BIT model was compared with the experimental data of Akbar et al. [233]. A good agreement was achieved in predicting the gas void fraction, phase streamwise velocity, and liquid phase Reynolds stresses (see Figure 7). As the computational power increases rapidly, the DNS study will play a more important role in uncovering the turbulence physics of bubbly flows. Therefore, more DNS data are expected.

**Figure 7.** Comparison between the model of Ma et al. [231] and the experimental data of Akbar et al. [233]: (**a**) gas void fraction, (**b**) streamwise velocities of the liquid and gas phases, (**c**) Reynolds normal stresses of the liquid phase, (**d**) Reynolds shear stress of the liquid phase.

#### *4.2. Supersonic Jet Transport*

Supersonic gas jets are applied in Basic Oxygen Furnace (BOF) steelmaking processes [234,235]. Due to its vital role in the refining efficiency [236] and the service lifetime of the furnace lining [237], it is imperative to understand the behavior of the supersonic

jet flow. During supersonic jet transport, it interacts with the surrounding medium to produce a turbulent mixing region, as shown in Figure 8. An accurate prediction of the growth rate of the mixing region is a challenge in modelling the jet transport. The standard *k*-*ε* turbulence model suffers from a typical weakness, i.e., an overestimated growth rate of the turbulent mixing region around the supersonic gas jet [238]. This leads to a large error in predicting the volume of cavity formed by the gas jet impinging on the liquid surface [239]. Several modifications have been proposed to solve this problem [240–242]. Sarkar et al. [240] considered that the compressible dissipation term should be included to account for the effect of the compressibility on the supersonic jet flow. Based on the asymptotic analysis and DNS data, an algebraic model was proposed for the compressible dissipation by using the turbulent Mach number. Subsequently, Sarkar [241] found that the reduced turbulence generation is the main cause of preventing the turbulent mixing layer from the growth rather than the effect of the dilatation term (pressure dilatation and compressible dissipation) on the growth. A gradient Mach number was induced to describe the inhibiting effect of compressibility on the turbulence growth rate. Heinz [242] used the gradient Mach number to evaluate the effect of compressibility on the turbulence distribution. An expression of the model constant *C*μ as a function of the gradient Mach number was obtained. However, the *k*-*ε* turbulence model with the compressibility correction failed to predict the supersonic jet flow with a temperature gradient, such as the potential flow core length and the cavity shape and dimension [243]. This is because the temperature gradient between the gas jet and the ambient is ignored by the original turbulence model. To sensitize the turbulence model with temperature fluctuation, efforts [244–246] have been made to modify the turbulent viscosity term or turbulent heat flux term. Abdol-Hamid et al. [247] corrected the turbulent viscosity term with the temperature gradient for the case of a hot gas jet entering into low temperature ambient. The results of modelling agreed well with experimental data for subsonic and supersonic jet flows. However, this model did not give reasonable results for the BOF supersonic jet flow [243], under which a room temperature gas jet enters into a high temperature ambient (see Figure 9a). Alam et al. [243] modified Abdol-Hamid's model for the gas jet of BOF. The predicted distributions of axial velocity and dynamic pressure along the central axis of the jet closely agreed with the experimental data (see Figure 9b,c). It is worth mentioning that the data compared in Figure 9 were extracted along the central axis of the jet (see the dashed line in Figure 8). Wang et al. [248] adopted the model of Alam et al. to study the multiple supersonic oxygen jets in the BOF process, where a better prediction of the shape of cavity caused by the gas jet impinging on the liquid surface was achieved. To fit the empirical data available for a cold jet to hot environment [249], Lebon et al. [250] derived another expression of *C*μ as a function of the enthalpy ratio of ambient gas to gas jet. It was found that the modified model is adequate to model a compressible jet to hot environment. In addition to the two-equation turbulence model, the LES model with the compressibility and temperature corrections is also widely used to study the jet noises [251]. Wang et al. [252] and Bodony et al. [253] critically reviewed the roles of the SGS model and inflow boundary condition in predicting jet noise. The open issues and future directions were also included in the papers. Therefore, the LES model for jet noise prediction is excluded from this article to avoid repetition.

**Figure 8.** Illustration of the supersonic jet transport.

**Figure 9.** (**a**) and (**b**) Axial velocity distribution, (**c**) dynamic pressure distribution along the central axis of the jet (*x*: distance along the central axis, *de*: nozzle diameter) [243].

#### *4.3. Electromagnetic Suppression of Turbulence*

A static magnetic field has been commonly used to stabilize the turbulent flow in the continuous casting process, aiming to improve the product quality. The application of the magnetic field not only suppresses the mean flow, but also dampens the flow turbulence. For the former, including the electromagnetic force/Lorentz force in the momentum equations, as reported by many studies [254–258], solves the suppression of the mean flow. The latter, also called Joule dissipation, can be resolved by adding the electromagnetic damping effect in a conventional turbulence model. In order to tailor the turbulence model towards magnetohydrodynamic (MHD) flows, many fundamental investigations have been performed [259–263]. Ji and Gardner [259] modified a standard low Re *k*-*ε* model [27] to account for the Joule dissipation. Additional source/sink terms were added in the *k* and *ε* equations, as well as a damping factor for turbulent viscosity. Both the new terms and the damping factor contain an exponential expression *e*−CN, where C is an empirical constant determined from experimental data, and N is the interaction parameter defined as the ratio of the time scale of large eddies (*L*/*U*) to the characteristic magnetic braking time (*ρ*⁄σB2). This interaction parameter represents the strength of turbulence damping due to the magnetic field. However, as reported by Kenjereš and Hanjali´c [260], the use of the time scale of large eddies makes the model dependent on the bulk flow properties, restricting the model to homogeneous magnetic fields and simple geometries. Replacing the bulk time scale with the local turbulent time scale (*k*/*ε*) can overcome the above-mentioned deficiency. As noted in the work of Ji and Gardner [259], the modified turbulence model failed to accurately predict the turbulent kinetic energy when the strength of the magnetic field increased. This is because the model is formulated based on the assumption of isotropic turbulence, while MHD flows exhibit strong anisotropy caused by the damping of turbulent fluctuations nonparallel to the magnetic field. The anisotropy of the Joule dissipation has to be considered in the turbulence model. Widlund et al. [261] proposed a scalar dimensionality anisotropy parameter carrying the information with regard to the magnitude and anisotropy of the Joule dissipation tensor. A scalar transport equation for the anisotropy parameter was also proposed. With this scalar transport equation, it is easy to extend a conventional turbulence model towards MHD engineering applications. It should be noted that the scalar transport equation is based on phenomenological reasoning, and special attention should be paid to the hydrodynamic part of the proposed equation. Miao et al. [264] coupled the scalar transport equation with the RANS SST *k-ω* turbulence model [10] to solve the MHD flow in a continuous casting mold. The Joule dissipation was taken into account by adding source/sink terms in the *k* and *ω* equations, and the anisotropic behavior of the Joule dissipation was represented by the anisotropy parameter contained in the source/sink terms. The computational domain of the continuous casting mold is shown in Figure 10, and the measured region and data extraction line for liquid metal (low melting point alloy Ga68In20Sn12) velocity are depicted in the figure. The averaged velocities obtained from the three-equation turbulence model (i.e., *k*-*ω*-anisotropy parameter) qualitatively and quantitatively agreed with the experimental data (see Figures 11 and 12). The turbulence models described above are modified in the framework of the RANS approach. Due to the time-averaging treatment in the RANS approach, the change of turbulence structure and transient flow behavior are difficult to be captured by the RANS-based turbulence model. Given that engineering flows are usually very complex and have a high Re number, the RANS-based turbulence models are still considered practical candidates for solving MHD flows.

**Figure 10.** The computational domain of the continuous casting mold and the measured region and data extraction line (figure adapted from [264]).

**Figure 11.** The time-averaged horizontal velocity in the midplane parallel to the mold wide face with electrically conducting mold walls, B0 = 310 mT (within the dotted lines): (**a**) experimental results, (**b**) numerical results from the SST *k*-*ω* model + Joule dissipation, (**c**) numerical results from the conventional SST *k*-*ω* model [264] (The measured region in Figure 10).

**Figure 12.** The time-averaged vertical velocity along a vertical line in the midplane at *x* = 0.045 m (B0 = 310 mT with electrically conducting mold) [264] (The data extraction line in Figure 10).

In addition to the RANS-based turbulence model, the LES model, which is able to capture turbulence structures, has been modified or developed to solve the MHD turbulent flow. Shimomura [262] incorporated the magnetic damping effect in the form of a locally determined damping factor for SGS eddy viscosity. Compared with the original Smagorinsky model [12], the new model showed better performance in both turbulent and laminar states. The new model also successfully predicted the anisotropic laminarization caused by the applied magnetic field, for which the RANS-based model is incompetent. Kobayashi [263] developed a new SGS model based on the coherent structure extracted by the second invariant of a velocity gradient tensor in grid scale flow field. Compared with the conventional Smagorinsky model [12], the new SGS model does not need an explicit walldamping function or change the model parameter depending on the flow. Compared with the dynamic Smagorinsky model [13], the new SGS model is numerically stable due to the fact that the model parameter is always positive. Considering the advantages of the new SGS model, Singh et al. [265,266] used it for the MHD turbulent flow in a continuous caster. The numerical results closely matched the measured data. Although the LES models can

capture transient flow behavior and turbulent structures, the computational cost required is much higher than that of the RANS-based turbulence model. With the fast increase in computing power, more LES simulations are expected for engineering flows in the near future. The DNS studies for MHD flows as presented by Chaudhary et al. [267,268] are also necessary since important information regarding turbulent quantities can be obtained for formulating the RANS or LES model.

#### **5. Conclusions**

Three main turbulence approaches (i.e., the RANS, LES, and DNS) have been reviewed in this paper. The formulations and variations of the RANS approach have been described, evaluated, and discussed. The eddy viscosity models still remain the most widely used methods for calculating simple engineering flows. The DSM overcomes the limitations of the eddy viscosity model and can be used to predict complex anisotropic flows. However, since a high degree of uncertainty is introduced in the DSM, special attention should be paid to its application. Although the LES is not completely ready for the calculation of high Re number engineering flow in the current stage, it can be used for the studies on fundamental turbulence physics and the low Re number flow in simple geometry. The LES is highly recommended to quantify the turbulent quantities and transient flow behavior. Since the inflow conditions are very important for LES modelling, a further improvement of the method is needed. With the improvement of the model formulation and accurate specification of inlet boundary conditions, the LES shows great potential to realistically solve complex turbulent flows. Due to the high computing power needed, the DNS is often adopted to understand the turbulence physics and to evaluate less advanced turbulence models, rather than being applied for a real case study. It contributes to the development of the turbulence model. Depending on the actual situation, the corresponding experimental validation should be adopted to ensure the fidelity of the used turbulence model before a formal simulation. To demonstrate the application of turbulence modelling, the conventional turbulence models have been customized to resolve three important turbulence-related flow issues, namely, BIT, supersonic jet transport, and electromagnetic damping. Success has been achieved by the customized turbulence model. The general outlook for turbulent flow simulation is listed as follows:

caused by magnetic fields. The mechanisms of turbulence modulation with respect to these aspects remain open for fundamental investigation.

**Author Contributions:** Conceptualization, Y.W.; writing, Y.W., and L.C.; review, and editing, Z.C., B.B. and M.G.; supervision, M.G. and B.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Belgian IWT project, grant number 140514, and China Scholarship Council, grant number<PHONE_NUMBER>05.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The financial support from Belgian IWT project (No.140514) and China Scholarship Council (CSC) is highly acknowledged. The authors also would like to give their thanks to the anonymous reviewers whose comments helped to improve the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**

#### **References**

### *Review* **A Review of Bubble Dynamics in Liquid Metals**

**Tim Haas \*, Christian Schubert, Moritz Eickhoff and Herbert Pfeifer**

Department for Industrial Furnaces and Heat Engineering, RWTH Aachen University, Kopernikusstraße 10, D-52074 Aachen, Germany<EMAIL_ADDRESS>(C.S<EMAIL_ADDRESS>(M.E<EMAIL_ADDRESS>(H.P.)

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Gas bubbles are of major importance in most metallurgical processes. They promote chemical reactions, homogenize the melt, or float inclusions. Thus, their dynamics are of crucial interest for the optimization of metallurgical processes. In this work, the state of knowledge of bubble dynamics at the bubble scale in liquid metals is reviewed. Measurement methods, with emphasis on liquid metals, are presented, and difficulties and shortcomings are analyzed. The bubble formation mechanism at nozzles and purging plugs is discussed. The uncertainty regarding the prediction of the bubble size distribution in real processes is demonstrated using the example of the steel casting ladle. Finally, the state of knowledge on bubble deformation and interfacial forces is summarized and the scalability of existing correlations to liquid metals is critically discussed. It is shown that the dynamics of bubbles, especially in liquid metals, are far from understood. While the drag force can be predicted reasonably well, there are large uncertainties regarding the bubble size distribution, deformation, and lift force. In particular, the influence of contaminants, which cannot yet be quantified in real processes, complicates the discussion and the comparability of experimental measurements. Further open questions are discussed and possible solutions are proposed.

**Keywords:** liquid metals; bubble generation; bubble size distribution; porous plugs; bubble deformation; drag force; lift force

#### **1. Introduction**

Gas injection and stirring is an integral part of most metallurgical processes. In some cases, the gas even determines metallurgical tasks. Examples are the removal of hydrogen in the degassing of aluminum or the removal of non-metallic inclusions from steel in the ladle or tundish. Thus, metallurgical bubble column reactors are of major importance for the process industry. Their main advantages are a simple construction, low maintenance costs, good mass, and heat transfer [1] and applicability in case mechanic stirring is prohibited by high temperature or reactive flows.

During the buoyancy-driven rise of bubbles, momentum is exchanged between the fluid and bubbles, which promotes a stirring effect. Characteristic of bubbles are their non-rigid phase boundaries. Thus, bubbles react to pressure gradients with a deformation. This leads to fascinating phenomena such as shape and path oscillation. These were already described by Leonardo DaVinci, which is why it is sometimes referred to as Leonardo's paradox [2]. On the other hand, this additional degree of freedom also leads to a very high complexity of the interaction between bubbles and fluid. Therefore, despite many years of extensive research, a comprehensive description of all phenomena is not yet possible. However, fundamental understanding of the bubble properties, shapes, rising velocities, or interfacial forces is crucial for the understanding, modelling, and optimization of bubble column reactors [3]. In metallurgy, the optimization of processes is often carried out using computational fluid dynamics (CFD). However, gaps in the knowledge of the behavior of gas bubbles in metals lead to considerable model uncertainties.

Bubble reactors are a multi-scale problem. It can be roughly divided into four different scales, though not all phenomena can be sharply assigned to one scale. On the microscale,

**Citation:** Haas, T.; Schubert, C.; Eickhoff, M.; Pfeifer, H. A Review of Bubble Dynamics in Liquid Metals. *Metals* **2021**, *11*, 664. https://doi.org/ 10.3390/met11040664

Academic Editors: Fernando Castro and Gunter Gerbeth

Received: 19 March 2021 Accepted: 14 April 2021 Published: 19 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the mass and heat transfer between fluid and bubble takes place. The bubble scale determines the formation mechanism of bubbles, bubble deformation, and interfacial forces. The swarm scale deals with the interaction of bubbles, as well as the characterization of the bubble column. At the macro or reactor scale, the stirring effect generated by the bubbles is considered. In order to limit the length of this paper, this review focuses on the bubble scale, largely omitting the other scales.

Different dimensionless numbers are used to describe bubble dynamics or bubble generation. The most important dimensionless numbers are the Reynolds number (Re), Eötvös number (Eo), Morton number (Mo), Weber number (We), and Froude number (Fr):

$$\begin{aligned} \text{Re} &= \frac{\text{Inertial force}}{\text{Viscous force}} = \frac{\text{u}\_{\text{C}}\text{L}\_{\text{C}}\text{o}\_{1}}{\mu\_{\text{I}}}\\ \text{Eo} &= \frac{\text{Buoyancy force}}{\text{Surface tension force}} = \frac{\text{g}(\rho\_{\text{I}} - \rho\_{\text{b}})\text{L}\_{\text{C}}^{2}}{\sigma} \\ \text{Mo} &= \frac{\text{g}\mu\_{\text{I}}^{4}(\rho\_{\text{I}} - \rho\_{\text{b}})}{\rho\_{\text{I}}^{2}\sigma^{3}} \\ \text{We} &= \frac{\text{Inertial force}}{\text{Surface tension force}} = \frac{\rho\_{\text{I}}\text{u}\_{\text{C}}^{2}\text{L}\_{\text{C}}}{\sigma} \\ \text{Fr} &= \frac{\text{Inertial force}}{\text{Gravity force}} = \frac{\text{u}\_{\text{C}}^{2}}{\text{g}\text{L}\_{\text{C}}} \end{aligned} \tag{1}$$

where uC and LC are the characteristic velocity and length, respectively. For this, different values are employed, though for bubbly flows, the bubble size, expressed as the diameter of a sphere with the same volume and the relative velocity between bubble and liquid, is commonly utilized. If other quantities are used, this is explicitly mentioned in the text.

In the following, different analytical methods and measurement techniques are discussed, and specific advantages, disadvantages, and restrictions are analyzed. Afterwards, different bubble scale phenomena are discussed. This review covers the bubble formation mechanism at nozzles and purging plugs, bubble deformation, and finally the interfacial forces between bubbles and fluid. Due to the small number of studies in liquid metals, conclusions drawn from aqueous systems are summarized first. In the second step, their transferability to liquid metals is discussed. This allows guidelines to be derived as to which models should currently be used for CFD, and the model uncertainty can be analyzed in more detail. Finally, open questions are identified and possible solution strategies are proposed.

#### **2. Measurement Methods**

Fundamental knowledge about the dynamics of bubbles can be achieved theoretically by analytical approximations, by experimental measurements, or by means of direct numerical simulation (DNS). The advantages and drawbacks of the different approaches are outlined briefly below.

#### *2.1. Theoretical Investigations*

Theoretical work is based on the analytical solution of the flow equations with certain boundary conditions at the phase boundary. However, this is only possible for special cases where certain simplifications can be made. Usually this is the assumption of spherical bubbles and inviscid (Re→∞) or creeping flows (Re→0). Due to these limitations, theoretical work is usually of limited use for industrial applications. However, these theoretical results serve as benchmarks for numerical models, the foundation for semi-empirical models, or reveal functional relationships.

#### *2.2. Experimental Measurements*

In experimental measurements, bubble characteristics are measured directly or indirectly and phenomenological conclusions can be drawn. However, there are numerous difficulties with measurements. First, the experimental boundary conditions are almost impossible to control. In some cases, influencing factors are unknown or cannot be determined. One example is the influence of contamination. It has been known for a long time that contamination has a significant influence on the bubble deformation, the rising

velocity, and the lateral movement of the bubbles. Therefore, distilled water or deionized water has often been used for measurements. However, some studies [4,5] show that with these apparently pure systems, the influence of contamination can still be significant. In liquid metals, contamination by oxidation can hardly be avoided. Another example is the influence of injection, which often affects the bubble shape or movement some distance above the injector. As these effects cannot be fully quantified yet, this results in a large measurement uncertainty and similar studies are often not directly comparable. This applies in particular to liquid metals where the experimental setups are usually very small.

For the actual measurements, different methods can be found in the literature, which can be categorized into intrusive and non-intrusive methods. Intrusive methods, such as conductivity probes or fiber optic probes, interact to some extent with the flow and are usually restricted to single or few bubbles. Thus, the results may be biased by the applied measurement technique. On the other hand, intrusive methods can be applied in a broader range of setups and with different liquids. Non-intrusive methods do not affect the flow but have more constrains regarding their applicability. Amongst the non-intrusive methods, imaging methods, either based on a multi-stage image processing [6–10] or on machine learning [11,12], are the most frequently used for non-opaque systems because multiple objects can be observed simultaneously and the equipment is comparably cheap and flexible.

A major difficulty is that a large number of different phenomena must be considered for a comprehensive analysis of rising gas bubbles. These takes place on different temporal and spatial scales and require different measurement techniques:

*Bubble path and velocity*: For a statistically significant path, recordings on a decimeter or better a meter scale are necessary. The frame rate should be at least 200 frames per second (FPS), whereas the internal memory of high-speed cameras limits the duration of the recording. To reconstruct the three-dimensional rise, which can be zigzagging or spiraling, at least two cameras are required.

*Bubble shape:* In order to analyze the bubble oscillation, recording must be made in the millimeter range with high frame rates. The inaccuracy of the reconstruction of the threedimensional shape of non-spherical bubbles from two-dimensional projections depends strongly on the number of cameras. Fu and Liu [13] showed that the bubble volume error of a single bubble is about 25% for a one-camera system and decreases to about 10% for a two-camera system.

*Velocity distribution of the fluid*: For a comprehensive understanding, not only the bubble but also the flow field around the bubble must be measured. For this purpose, the fluid must be seeded with some tracer, which in turn is likely to influence the flow. Moreover, the wake can stretch from a millimeter to a few centimeters.

*Bubble swarms*: A particular challenge for postprocessing are bubble swarms with high void fractions. Bubbles may overlap on the two-dimensional projection and are visible as a cluster of bubbles. Their segmentation is usually achieved by multi-staged image processing techniques [6–10] or machine learning [11,12]. The accuracy of the different approaches are discussed in [14]. However, no generally applicable workflow was proposed yet. Thus, experimental measurements are often restricted to relative dilute bubble swarms.

A simultaneous investigation of all phenomena is therefore only possible with many cameras, additional equipment, and a sophisticated triggering and could not be realized so far. Different approaches for non-opaque systems are discussed by Bröder and Sommerfeld [7].

In liquid metals, where the fluid is opaque and high temperatures and corrosivity prevent most measurement techniques, measurements are particularly cumbersome. Most early studies employed indirect measurement methods, for example by analyzing the pressure fluctuation [15,16] or the noise [17,18] when a bubble was released or by measuring the time between the generation of the bubble and the breaking of the bath surface [19]. However, these measurements are characterized by a high measurement uncertainty and only

allow a macroscopic analysis of most phenomena. Furthermore, their applicability is limited to a few bubbles or low gas flow rates. Therefore, subsequent studies (e.g., [16,20,21]) used resistivity measurements. However, this technique also has disadvantages. First, the method is intrusive and provides only local information [22]. Second, the sensors have a short lifetime, especially at high temperatures [21]. Third, it is difficult to establish the functional relationship between the measured chord (or pierced) length and the actual bubble size [23]. On the other hand, this method enables the investigation of significantly higher flow rates. In addition, multi-needle sensors, like those developed by Iguchi et al. [24], allow a reconstruction of the bubble shape. X-ray measurements (e.g., [25–27]) are nonintrusive and allow an image-based analysis like described above. However, given the high absorption of the liquid metals, measurements are restricted to quasi two-dimensional experiments where the setup's thickness is not more than 12 mm [27]. Neutron radiography (e.g., [28,29]) employs a similar technique but allows thicker experimental setups up to 25 mm [30]. On the other hand, the resulting images are noisier and have less contrast than X-ray images, making subsequent image processing more difficult. Sarma et al. [29] showed that neutron radiography is capable of measuring the motion of tracers in the fluid, which would allow simultaneous measurements of fluid velocity and bubble velocity. This method, called NeuPIV, could provide interesting insights into bubble dynamics in liquid metals in the future. Other non-intrusive measurement techniques are ultrasound Doppler velocimetry (e.g., [31,32]) and inductive methods (e.g., [33,34]). However, they can only be used to measure the bubble velocity and the evaluation of the measurement signals can be very difficult with a higher number of bubbles. A problem with many measurement techniques is that sensors must be placed in close vicinity of the experimental setup, which prevents high temperature measurements. Therefore, measurements are often limited to Galinstan (GaInSn) or Mercury (Hg), as it is already liquid at room temperature.

Due to the described measurement difficulties in liquid metal, the majority of studies have been performed in aqueous systems. However, even for these systems, a comprehensive understanding of the complex bubble–fluid interaction by experimental measurements alone is almost impossible.

#### *2.3. Direct Numerical Simulation*

Due to the constant increase in the computing capacity of modern high-performance clusters, the possibility of an analysis by means of DNS is opening up. Here the phase boundary between fluid and gas bubble is directly modelled by means of a suitable interface tracking method. Thus, many phenomena can be analyzed far better than with experimental measurements, as quantities are directly accessible. Furthermore, the fluid properties can be easily varied, so that fluids can be examined, which cannot by experiments. However, DNS is a mathematical model that is affected by modelling and numerical errors. In particular, the interface tracking method and the treatment of surface tension are still being continuously developed. Therefore, existing results from DNS calculations must be critically examined for the influence of these models. They are therefore largely omitted in this review. Furthermore, the analysis is still limited by the computing capacity. According to Kolmogorov's length scales, the number of necessary cells increases with Re9/4. This implies that computations are limited to relatively small computational areas with low turbulence (low Re) and only a few bubbles.

#### **3. Bubble Formation Mechanisms**

In most metallurgical processes, gas bubbles are generated by single nozzles (e.g., regular nozzles, open lances, Laval nozzles, impellers) or purging plugs. Gas bubbles generated by metallurgical reactions may be assigned to the micro scale and are omitted in this review. The bubble size distribution is an important parameter of the metallurgical process. It directly affects the removal efficiency of non-metallic inclusions [35] and degassing reactions [36] and can influence the fluid flow [37,38]. Generally, small bubbles and a narrow size distribution are desirable. In addition, the bubble size distribution

is an important boundary condition of numerical models. Despite that importance, the bubble size distribution in actual processes is still largely unknown and comprehensive theoretical models have not been derived yet. A straightforward approach would be to assume bubbles in the model as well as in the melt in pressure equilibrium. A scaling of the bubble size from a physical model (M) to the real process (R) could then be maintained by [39]:

$$\mathbf{d}\_{\mathbf{b},\mathbf{R}} = \mathbf{d}\_{\mathbf{b},\mathbf{M}} \sqrt{\frac{\sigma\_{\mathbf{R}} (\rho\_{\mathbf{l}} - \rho\_{\mathbf{b}})\_{\mathbf{M}}}{\sigma\_{\mathbf{M}} (\rho\_{\mathbf{l}} - \rho\_{\mathbf{b}})\_{\mathbf{R}}}} \tag{2}$$

However, this approach implies that the bubble size is independent of the injector, which contradicts most experimental results for relevant flow rates.

Since the generation of gas bubbles by injectors plays an important role in a variety of processes, there is a great deal of research on the generation of gas bubbles in water or aqueous systems. These were discussed in the review by Kulkarni and Joshi [40]. The majority of these studies deal with the generation of gas bubbles at nozzles. Studies on the generation of gas bubbles on porous plugs or slot plugs are scarce. Very few studies, summarized in Table 1, investigated the generation of gas bubbles in molten metals, since measurements are cumbersome. These studies have common ground in the fact that they investigate the generation of bubbles at single nozzles. To the best of the author's knowledge, there are no studies on the generation of gas bubbles at purging plugs in liquid metals.

#### *3.1. Bubble Generation at Single Nozzles*

Before discussing the results in liquid metals, some important observations and correlations from water models are introduced, as these are often used to check the scalability of water experiment correlations to liquid metals or serve as foundation for correlations for liquid metals.

In an early and widely recognized study, Leibson et al. [41] correlated the bubble size with the orifice Reynolds number. They found two regimes. At laminar flow conditions, the bubble size increases with increasing orifice Reynolds number Reor, where the nozzle diameter and the exit velocity are the characteristic scales, and depend on the orifice inner diameter dni. In contrast to that, at turbulent outflow conditions, the bubble size is considerably smaller and almost independent of the Reynolds number and the orifice diameter.

$$\mathrm{d}\_{\mathrm{b}} = 0.029 \mathrm{d}\_{\mathrm{ni}}^{0.5} \mathrm{Re}\_{\mathrm{or}}^{0.33} \mathrm{Re}\_{\mathrm{or}} < 2100 \mathrm{d}\_{\mathrm{b},32} = 0.0071 \mathrm{Re}\_{\mathrm{or}}^{-0.05} \mathrm{Re}\_{\mathrm{or}} > 10000 \tag{3}$$

However, all comparisons with measurements in liquid metals [17,42] showed poor agreement with Equation (3), so that it can be concluded that a scaling by the Reynolds number alone is insufficient.

Numerous studies on the generation of gas bubbles at nozzles have shown that a distinction can be made between the constant volume regime and the constant frequency regime [40]. In the former, an increase of the flow rate leads to an increase in bubble frequency, but the bubble size remains constant. In that regime, the bubble size is given by Tate's law [43]:

$$\mathbf{d}\_{\mathbf{b}} = \left(\frac{6\sigma \mathbf{d}\_{\mathrm{nl}}}{\varrho\_{\mathrm{l}} \mathbf{g}}\right)^{1/3} \tag{4}$$

With a further increase beyond a critical flow rate, the bubble detachment frequency remains constant but the bubble size increases. For that regime, Davidson and Amick [44] derived a correlation for intermediate gas Qg flow rates up to 250 cm3/s:

$$\mathbf{d}\_{\rm b} = 0.0054 \left( 10^7 \mathbf{Q}\_{\rm g} \mathbf{d}\_{\rm ri}^{0.5} \right)^{0.289} \tag{5}$$

Mersmann [45] derived a correlation which covers both regimes:

$$\mathbf{d}\_{\mathbf{b}} = 0.01 \cdot \left[ \frac{3 \cdot 10^6 \,\mathrm{ord}\_{\mathrm{ri}}}{\rho\_{\mathrm{l}} \mathrm{g}} + \left( \left( \frac{3 \cdot 10^6 \,\mathrm{ord}\_{\mathrm{ri}}}{\rho\_{\mathrm{l}} \mathrm{g}} \right)^2 + \mathrm{K} \frac{10^6 \mathrm{Q}\_{\mathrm{g}} \mathrm{d}\_{\mathrm{ri}}}{\mathrm{g}} \right)^{0.33} \right]^{0.33} \tag{6}$$

where K is an empirical constant.

Investigations on gas bubble generation in liquid metals, as well as their boundary conditions and the most important results, are summarized in Table 1. Hg [15,20,42], silver (Ag) [15], iron (Fe) [18,46], copper (Cu), lead (Pb), tin (Sn) [17], Wood's metal [47], Li17Pb83 [16], GaInSn [22], or aluminum (Al) [48] were investigated. Thus, a broad range of physical properties is covered.

The earliest studies used indirect measurement methods. Here, either the pressure pulse in the gas supply was measured, which results from the rupture of the gas bubble [15,16], or the resulting vibration [17,18]. The bubble volume was calculated from the detachment frequency and the flow rate. Due to this indirect measurement technique, however, the gas flow rate was limited to rather small values, so that the rupture signals of individual gas bubbles could be distinguished.

Sano and Mori [15] compared the bubble size in liquid silver and mercury with theoretical correlations (Equations (4)–(6)) derived in aqueous systems. They reported good agreement between their measurements and theory. However, they pointed out that there is a difference in the wettability of the injector. Therefore, the inner diameter of the injector dni in Equations (4)–(6) has to be replaced by the outer diameter dno in non-wetting systems. Irons and Guthrie [18] made similar observations, but found less satisfying agreement with the correlation of Davidson and Amick (Equation (5)) or Mersmann (Equation (6)). This was attributed to the fact that the experimental setup, in particular the chamber volume of the injector, differed between the experiments, which had a large influence on the resulting bubble size. Interestingly, both measurements confirm the observations from the water experiments that the bubble size becomes independent of the fluid's material properties for sufficiently large flow rates. Thus, this observation seems to be valid for systems with significantly higher surface tensions and density differences like molten metals. Based on the correlation of Mersmann (Equation (6)), Mori et al. [46] derived a first correlation of the bubble size in liquid metals (note that this correlation was first published in 1977 in Japanese by Sano and Mori and is therefore often called Sano and Mori equation):

$$\mathbf{d}\_{\mathbf{b}} = 0.01 \cdot \left[ \left( 10^6 \frac{6 \sigma \cdot \mathbf{d}\_{\mathrm{no}}}{\rho\_{\mathrm{l}} \mathrm{g}} \right)^2 + 0.0242 \left( 10^7 \mathrm{Q}\_{\mathrm{g}} \mathrm{d}\_{\mathrm{no}}^{0.5} \right)^{1.734} \right]^{1/6} \tag{7}$$

Note that this correlation was originally derived using cm and dyn. The use of SI units can lead to deviations of the order of ±5%. This equation was derived for liquid iron and validated with values from previous measurements for Hg, Ag, and H2O. It should be noted, however, that this correlation was derived for very small flow rates Qg ≤ 70 cm3/s. This is primarily due to the indirect measurement technique. This has the additional disadvantage that no conclusions can be drawn whether the bubbles disintegrate during injection or on their rising path. Therefore, subsequent studies used resistivity measurements.

The validity of Equation (7) has been critically reviewed in several subsequent studies, where different results were found. Mori et al. [20] found good agreement, even for higher flow rates. Irons and Guthrie [18] found larger gas bubbles than predicted by Equation (7), but this was attributed to differences in chamber volume, as mentioned above. A similar result was found by Iguchi et al. [21]. The reasons for that can only be speculated, since some details describing the experimental setup are not given in their paper. Tschuchiya et al. [16] also found different values. Zhang et al. [49] compared the values of Equation (7) with measured values of bubbles generated with a porous plug in a water model. For that, the diameter in Equation (7) was replaced by the pore size. However, poor agreement was found. The gas bubbles generated by the porous plug were significantly smaller than predicted by Equation (7).

Moreover, there are also studies that contradict the observation that the initial bubble size is independent of the fluid properties. Tsuchiya et al. [16] argues that the fitting parameters of Equation (7) depend on the fluid properties and the nozzle orientation. However, it should be noted that the number of measurements was so small that these parameters might be overfitted. Andreini et al. [17] found relatively poor agreement with water model correlations, in particular with Equation (3). They argued that the bubble size should not be scaled with the orifice Reynolds number but with the orifice Froude and Weber number, where the characteristic length and velocity are the nozzle diameter and the exit velocity of the gas. They fitted their data to:

$$\mathbf{d}\_{\rm b} = \mathbf{d}\_{\rm ni} \text{Fr}\_{\rm or}^{0.224} \text{We}\_{\rm or}^{-0.109} \tag{8}$$

However, Irons and Guthrie [18] point out that the poor agreement with the water model correlations may also be due to different nozzles. Unfortunately, these parameters are not given in [14], so this cannot be further verified. Iguchi et al. [26] measured the bubble size by X-ray imaging. The results of their study suggest that the bubble size depends not only on gas flow rate but also on gas and fluid properties. However, the latter point is contradicted by the measurements of Xie et al. [47].

Resistivity measurements for higher flow rates suggest that the detachment bubble size can be reasonably predicted by Equation (5) [42]. However, the actual bubble size is strongly determined by bubble breakup and coalescence [42,47]. Therefore, Sano and Mori [42] correlated the actual bubble size with a proportional factor that takes into account bubble breakup and coalescence and fitted their data to:

$$\mathrm{d}\_{\mathrm{b},32} = 0.00091 \left( 10^6 \frac{\sigma}{\rho\_{\mathrm{l}}} \right)^{0.5} \left( \frac{200 \,\mathrm{Q}\_{\mathrm{g}}}{\pi \mathrm{d}\_{\mathrm{ni}}^2} \right)^{0.44} \tag{9}$$

It should be noted that this correlation was derived for the Sauter mean diameter. However, as shown in Figure 1, the deviation is small in case the bubble distortion is small.

**Figure 1.** Arithmetic mean and Sauter mean diameter as a function of bubble aspect ratio.

Xie et al. [47] investigated bubble swarms in Wood's metal and reported that the bubble size distribution followed a log-normal distribution. For the mean bubble diameter, they correlated:

$$\mathbf{d}\_{\mathbf{b}} = 0.000146 \left[ 10^{4} \frac{\mathbf{Q}\_{\text{g}}}{\mathbf{z} + \mathbf{H}\_{0}} \right]^{0.1} \tag{10}$$

where H0 is the distance between the nozzle and the mathematical origin of the plume in m (correlation given in [47]).

In their study, no significant influence of the gas properties or the nozzle size was found. Using a similar experimental setup, Iguchi et al. [50] confirmed the applicability of Equation (10). However, a major criticism of Equation (10) is that it was only derived for one fluid and certainly overfits this system. Fluid property dependent bubble coalescence or breakup, which was found to be important, cannot be considered. Thus, Equation (10) is unlikely to generalize to other fluids.

The studies carried out so far on the formation of gas bubbles in liquid metals show that the theory established for water models for the formation at single nozzles at low gas flow rates also applies to liquid metals as long as the different wettability of the nozzles is considered. At very low flow rates, bubbles of identical size are formed, depending on the fluid properties and the injector. At slightly higher flow rates, however, the bubble size depends primarily on the flow rate and not on the injector or the fluid properties. Because of the difficult measurement conditions, especially at high flow rates, this observation refers to very low flow rates, and for more relevant high flow rates there are very few measurements [42,47]. At these high flow rates, the bubble size distribution seems to be largely dominated by bubble coalescence and breakup. However, there are too few studies in this area to evaluate whether existing models are applicable to liquid metals. Although both studies derived an empirical equation to determine the bubble size, it can be assumed that these correlations significantly overfit the measured data and that their scaling to other fluid systems or to higher flow rates is not valid. At very high flow rates, single nozzles go into jetting mode and the bubble size distribution results from breakup of the jet and sometimes subsequent coalescence [51].

**Table 1.** Studies on bubble generation in liquid metals (Pres Pul = indirect pressure pulse method, Ac = indirect acoustic method, Res = Electric resistivity, dno = Outer diameter of the nozzle, dni = Inner diameter of the nozzle, BSD = Bubble size distribution).

#### *3.2. Bubble Generation from Porous Plugs*

So far, bubble formation in liquid metals at single nozzles has been discussed. However, transferring conclusions or even correlations to industrial vessels in which purging plugs are employed is difficult. It is known from water model experiments that from single nozzles and purging plugs, completely different bubble size distributions are generated, especially at high purging rates. Therefore, bubble formation by purging plugs is discussed in the following. The results presented are all based on observations in aqueous systems or theoretical considerations. Thus, the difficulty of scaling those results to other material systems arises. To the best of the authors' knowledge, measurements in liquid metals that could verify the scalability of existing models are not available.

Early studies [52] observed that at low flow rates, more pores become activated with increasing flow rate while the bubble size remains approximately constant. When all pores are activated, a high void fraction occurs in the direct vicinity of the porous plug. At a certain flow rate, the bubbles start to coagulate to form larger bubbles. A descriptive model for low flow rates was developed by Loimer et al. [53]. Here, a wetted sieve tray is considered first and the theory is then modified to represent the more complex case of a wetted and later a non-wetted porous plug. The model assumes a chamber below the sieve tray in which a certain isostatic pressure prevails, which is dependent on the gas flow rate. As soon as this pressure exceeds the capillary pressure in one of the orifices of the sieve tray, the first orifice gets activated and gas flows through it. If the volume flow is sufficient to lower the pressure in the chamber, the volume flow stops and the pressure increases again. If the volume flow is not sufficient to lower the pressure, successive orifices are activated. In mathematical terms, the minimum flow rate is given by:

$$\mathbf{Q}\_{\rm min} = 2\pi \sqrt{\frac{\sigma \mathbf{r}\_{\rm or}^3}{\rho\_{\rm g}}} \tag{11}$$

The number of activated orifices is given by:

$$\mathbf{m} = \frac{\mathbf{Q}\_{\mathbf{g}}}{\mathbf{Q}\_{\text{min}}} \tag{12}$$

For a porous plug, the assumption of an isostatic pressure cannot be made. Instead, it is assumed that the pressure around an activated pore decreases with 1/R [53]. Furthermore, the capillary pressure in a pore is different from those in a simple orifice. With these considerations, Equation (11) becomes:

$$\mathbf{Q}\_{\rm min} = \frac{4\pi\sigma\kappa}{\mu\_{\rm K}}\tag{13}$$

where κ is the permeability of the porous plug. For non-wetted porous plugs, the bubble formation mechanism is different as shown in Figure 2. Bubbles can migrate vertically over the porous plug [54] and coagulate on the plug surface already [55], which results in the release of larger bubbles. Equation (13) has to be extended by a factor ymax/ror, where ymax is the maximum radius of the contact line, which depends on the wetting angle θ. The contact line is the place at which the plug, the bubble, and the liquid are in three-phase contact. The relation between ymax and θ is summarized in [53].

**Figure 2.** Bubble generation from a wetted (**a**) and non-wetted (**b**) sieve tray (Reproduced from [53], with permission from Elsevier, 2021).

Applying their theories to some typical plugs, Loimer et al. [53] concluded that, similar to the bubble formation at single nozzles, the bubbles at porous plugs are first formed in a quasi-steady regime at which the bubble size is given by Tate's law (Equation (4)). In case the superficial velocity, that is to say the quotient of flow rate and active plug area, exceeds a critical value, the bubble size starts to increase with the flow rate.

However, a rather mediocre agreement between experiments and theory was found. This was explained by the anisotopic permeability of the medium, which can cause pressure gradients in the chamber and in the porous material, respectively. Nevertheless, this model provides a descriptive idea of bubble formation at low flow rates. For higher flow rates, however, bubble coalescence in the vicinity of the porous plug gains importance. According to the theory of Koide et al. [56], all bubbles are formed with a similar size. However, in the vicinity of the plug, the void fraction is high and increases with the flow rate. At the same time, bubble velocity is low since bubbles have not reached their terminal rising velocity. This promotes coalescence of bubbles, which becomes more significant with increasing flow rates.

Since different physical mechanisms are dominant, it is useful to divide the bubble formation into different regimes depending on the gas flow rate. A first distinction [57] divided bubble formation into four different regimes: Quiescent column of discrete bubbles, pulsating column, onset of coalescence, and 'blanketing' at which the entire plug is covered by bubbles and coalescence is dominant. Images of the different regimes are shown in Figure 3.

The transition between the regimes was attributed to the superficial gas velocity. The transition from the quiescent to the pulsating column regime was found at 16 cm/s [57], 15 cm/s [59], or 14 cm/s [58]. The onset of coalescence was found between 30 and 38 cm/s [57], above 25 cm/s [59], or 40 cm/s [58]. 'Blanketing' was found above 67 cm/s [57]. Furthermore, it was shown that the superficial velocities at which transitions take place depend on the pore size, while the ambient pressure or the plug diameter had little effect [60]. For comparison, the superficial gas velocity in liquid steel during soft bubbling is approximately 90 cm/s. However, these studies investigated the same fluid system. Therefore, the transferability of these results to other fluids should be made rather critical. Furthermore, the differences between the various results, especially for the onset of coalescence, suggest that other factors besides the superficial gas velocity need to be considered.