Mathematical Methodology and Metallurgical Application of Turbulence Modelling: A Review

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.


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][2][3][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][7][8][9][10][11], the Large Eddy Simulation (LES) [12][13][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 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][7][8][9][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.

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).
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.

∂ρ ∂t
Metals 2021, 11, 1297 3 of 32 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][24][25].

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).
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][27][28][29][30][31][32], an example proposed by Lam and Bremhorst (LB model) [31] is presented in Equations (7) and (8). where: 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 µ e f f 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].

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][37][38][39][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].
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.  [41][42][43][44][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][45][46]. where: In addition to the modified versions of the linear EVM, the idea of non-linear EVM has been substantially used [47][48][49][50][51][52][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.

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][55][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.
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 s ij term is user-defined for a specific stress transport source. If there is no source, this term becomes zero. It is required to model D T,ij , G ij , φ ij and ε ij to close the equation, while it is not necessary to model L ij , C ij , D L,ij and P ij .
The D T,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. Sub-sequently, 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 G ij 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][70][71] and nonlinear models [72][73][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.

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.
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][13][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. where: In Equation (20), where: 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][83][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][98][99][100], the velocity estimation model [101,102], the Approximate Deconvolution model (ADM) [103][104][105], and the regularization model [106][107][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][98][99][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][104][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][115][116][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][121][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][124][125][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.
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][124][125][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.

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

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

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][143][144]. Several comparative studies have been performed on different synthetic methods [143][144][145][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.  [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.

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 gasliquid 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][157][158][159][160][161] and in providing complementary information for experimental study [162][163][164].

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 gasliquid 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][157][158][159][160][161] and in providing complementary information for experimental study [162][163][164].

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].  [169,[184][185][186][187][188] PIV [189][190][191][192][193]
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].

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.

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][221][222][223][224][225]. Both methods can be used with the RANS-based approach (k-ε/ω/SST (Shear-Stress Transport) [224] and RSM, [226][227][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. [226][227][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. .

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][241][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][245][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.

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][255][256][257][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][260][261][262][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 (ρ⁄σB 2 ). This interaction parameter represents the strength of turbulence damping due to the magnetic field. However, as reported by Kenjereš and Hanjalić [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.

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][255][256][257][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][260][261][262][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 (ρ⁄σB 2 ). This interaction parameter represents the strength of turbulence damping due to the magnetic field. However, as reported by Kenjereš and Hanjalić [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.

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][255][256][257][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][260][261][262][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 (ρ⁄σB 2 ). This interaction parameter represents the strength of turbulence damping due to the magnetic field. However, as reported by Kenjereš and Hanjalić [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. 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.    [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 wall-damping 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   [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 wall-damping 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  [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. In addition, detailed experiment data, especially the turbulent quantities, are expected for validating numerical simulation of MHD flow.

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: 1.
Due to insufficient knowledge of the turbulence physics, there is a high degree of uncertainty in modelling the higher-order correlations in a DSM simulation. It is necessary to perform fundamental studies on turbulence physics, and the DNS is an important method for those studies. Besides, advanced numerical schemes are expected to minimize the diffusive errors when solving the non-linear higherorder correlations.

2.
An inflow condition method, which can generate the required realistic turbulence characteristics at a reasonable computational cost, is needed for an accurate LES simulation.

3.
It is of significant importance to experimentally obtain accurate turbulence data such as the Reynolds stress term, turbulence kinetic energy, and its dissipation rate, especially for a densely dispersed flow. Advanced non-intrusive experimental techniques are needed to measure the flow field information and phase distribution. The obtained data can be used, for instance, to develop a reliable model considering the effect of the dispersed phase. 4.
The BIT is important for the turbulent flow. Due to a small size of the bubble, the induced turbulence is characterized with small spatial and temporal scales. An approach (experimental technique and/or numerical simulation) is required to obtain insight into the induced turbulence.

5.
Comprehensive and realistic modelling is indispensable for engineering applications, such as the inclusion of bubble coalescence and breakage, the consideration of temperature and pressure affecting the bubble behavior, the involvement of compressibility and the temperature gradient in supersonic jet flow, and the turbulence damping caused by magnetic fields. The mechanisms of turbulence modulation with respect to these aspects remain open for fundamental investigation. 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.