{'fraction_non_alphanumeric': 0.029915240152900115, 'fraction_numerical': 0.014542130629882001, 'mean_word_length': 5.023123123123123, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Turbulence modeling and simulation advances in CFD during the past 50 years

This paper is a short retrospective review of the predictive methods of turbulent flows in Computational Fluid Dynamics over the last 50 years since the first development of computers. The different schools of turbulence modeling are presented with the aim to guide both users and researchers involved in numerical simulation of turbulent flows.


Introduction
Most of the turbulence numerical predictions are based on fully or partially-averaged statistical equations. Due to the turbulence closure problem, the equations to be solved more or less require modeling. Although one can conceptually distinguish the physical modeling corresponding to the constitutive laws of the turbulence from the numerical modeling corresponding to the algorithmic technique of the computational problem to solve, the two concepts are closely linked and progress together. The development of powerful computers allows solving more complex models. This retrospective landscape tries to point out the capabilities of the different models and the main cornerstones of their evolution in CFD research in relation to developments in computational capabilities. Evidently this review cannot be exhaustive at all, and references are given for the interested reader. One of the main properties of turbulence in fluids is an efficient mixing, giving rise to many studies in heat and mass transfer and countless applications in industrial and environmental flows.

Description and approaches in turbulent flows
The turbulent field is usually supposed to be an unsteady solution of the Navier-Stokes equations. Direct numerical simulation (DNS) of turbulent flows requires huge computer power and even nowadays is not possible for complex real flows [1]. Statistical modeling such as the Reynoldsaveraged Navier-Stokes (RANS) equations has long been the main practical way to get useful flow predictions in spite of some weaknesses. This is to this day often the main method used for the simulation of turbulent flows encountered in engineering and industrial applications [2][3][4]. Large eddy simulations (LES) have been developed in which the fine-grained eddies are modeled while the filtered flow field is simulated. This modeling was initially developed to simulate atmospheric and geophysical flows and is now more and more used to get particular details and some insights on flow structures [5,6]. The Scale resolving simulation (SRS) methods including DNS and LES are beginning to be applied in aerospace industries [7,8]. In the past two decades, more recent hybrid RANS/LES models that combine in various ways the RANS and the LES methods have been proposed for simulating turbulent flows of practical interest allowing a second life to RANS closures. These models take into account the advantages of these two methodologies [9][10][11]. In parallel with these methods, spectral turbulence models are applied in the spectral space that are mainly used to study laboratory flows from a fundamental point of view with emphasis on the physical aspects of turbulence [12,13] based on two-point statistics. According to these authors, each of these methods has its own respective advantages and limitations and specific field of application so that they should be considered as complementary tools in computational fluid dynamics. In that sense, the most appropriate tool does not depend merely on the intrinsic performances of the method itself but more precisely on both the required computational resources, the nature of the flow and the question to address.

Principles of RANS
RANS are based on statistical averaging of the instantaneous Navier-Stokes equations giving rise to an open hierarchy of equations of moments. Numerous varieties of closure models have been developed depending on the level of closure and the approximations of unknown moments [2][3][4][5]. The mean value is theoretically defined as ensemble averaging and approximated by time or space averaging in experiments. Each variable φ in incompressible turbulence, is then decomposed into a mean part in a statistical sense 〈φ〉 and a fluctuation part φ . The mean value 〈φ〉 is then computed in practice by means of time averaging 〈φ(x)〉 = 1 In (1), the period of time T is assumed to be long in comparison with the characteristic turbulent time scale τ = k/ , where k denotes the turbulent kinetic energy and , its dissipation rate, i.e., T τ. Physically, this method relies on the fact that each mean flow variable describing the flow properties is associated to the integration of the entire energy density spectrum E (κ) over the wave number range κ ∈ [0, ∞[. In practice, the integration is being performed over all the small to the large wave numbers so that individual eddy scales cannot be distinguished. The transport equations (Reynolds equations) of the mean statistical velocity in incompressible flows, reads where τ ij = 〈u i u i 〉 − 〈u i 〉〈u j 〉 is the Reynolds stress tensor. In (2), ν, ρ, p, u i , denote the molecular viscosity, density, pressure, velocity component, respectively. The equations are then solved numerically in conjunction with modeled equations of moments up to n's order (n rarely exceeds n = 2).

Principles of LES
LES is based on the filtering process applied to the Navier-Stokes equations in order to distinguish large-scale fluctuations from the more universal fine-scale turbulence [5,6]. Mean values are obtained from averaging the unsteady solution in a post-treatment. Large eddies and Lagrangian tracers can be used to study the detailed behaviour of the flow. In this method, each variable φ is then decomposed into a filtered partφ and a small-scale fluctuation part φ > , but as a result of importance, in generalφ =φ and 〈φ > 〉 = 0. Contrarily to the RANS method, each filtered variable is associated with the integration of the energy density spectrum E (κ) over the wave number range κ ∈ [0, κ c ] where κ c is the cutoff wave number defined by means of the local grid size ∆ as κ c = π/∆, while the subgrid-scale (SGS) variable is computed from an integration over the range [κ c , ∞[. In practice, the largest wave number that can be simulated is given by the smallest grid step size of the mesh spacing. The filtering operator [14] applied to the Navier-Stokes equations allows to dissociate the simulated large eddies from the smaller eddies which are modeled. These calculated large eddies being a part of the fluctuation are random and the filtering process does not retain the nice operative properties of the statistical mean value [15][16][17].
The analytical definition of the filter in physical space is generally defined as which reduces to a simple convolution in homogeneous turbulence when the grid spacing ∆ becomes constant in space. Formally, the filtered motion equation then reads [16] ∂ū i ∂t where (τ ij ) sfs = u i u j −ū iū j is the subfilter-scale (SFS) stress tensor and where D/D t = ∂/∂t +ū j ∂/∂x j .

Principles of hybrid RANS/LES
Hybrid RANS/LES reconciles the two previous methods by introducing a mechanism allowing to switch from one method to the other depending of the zones of interest in the flow [9][10][11].

Mathematical framework of RANS and LES
As a result, the RANS and LES momentum equations (2) and (4) take exactly the same mathematical form if the commutation terms appearing in the function E arising from the filtering process in the material derivative of any variable are neglected [16][17][18]. In that sense, the difference between the RANS and LES approaches relies in practice on the closure of the equations and not on the basic equations themselves. In fact, a key difference is that SGS models for LES usually depend on the grid spacing, but RANS models do not. The numerical solution of (2) in RANS or (4) in LES requires to determine the turbulent stress τ ij or (τ ij ) sfs all along the calculation. This issue to address is known as the turbulence closure problem. Two main basic levels of modeling are currently used in computational fluid dynamics and transport processes. The first approach is based on the eddy viscosity models (EVM) known as first-order models while the second refers to second-moment closure (SMC) known as secondorder models [2][3][4]. Although RANS, LES and hybrid RANS/LES methods have been developed independently of each other, it is possible to see some connections between these different schools of modeling [10]. This is a topic that is receiving more and more attention these days.
From a numerical point of view, it is worth mentioning that RANS closures need robust numerical methods to solve stiff turbulence transport equations while LES require high-precision methods with good conservative properties to perform correctly.

The pre-computer era
The first attempts towards turbulent flow prediction based on the statistical closure of the Reynolds equations were made by Taylor in 1915 [19] who introduced a turbulence eddy viscosity represented as the product of a characteristic length and a characteristic velocity. But the wellknown Prandtl mixing length model introduced by Prandtl in 1925 [20] proved to be particularly fruitful. It reads ν t = l 2 m (2S ij S ij ) where S ij is the mean rate of strain and l m the mixing length to be prescribed in each particular flow. In 1932, Taylor proposed that the shear stresses are created by vorticity transfer rather than momentum [21]. A new theory of Prandtl in 1942 [22] later represented eddy viscosity more simply by the product of the width of the mixing zone and the difference in mean velocities across the shear layer. Also, Reichardt [23] supposed that the shear stress was proportional to the gradient of momentum flow rate rather than the gradient of mean velocity. All these methods can be viewed as algebraic eddy viscosity methods based on intuitive physics. These are relatively simple and when joined with similarity hypotheses can lead to analytical solutions of basic 2D flows. Such solutions flourished in this early period, before around 1965. Many examples of such analytical solutions can be found in [24]. Another example of simple representation is the well-known logarithmic boundary layer in which the shear stress is uniform and can be related to the friction velocity at the wall.

Extension of algebraic eddy viscosity and mixing length approaches to real flows
With the emergence of scientific computers, the numerical prediction of turbulent flows was developed in several research groups as a complement to experimental studies. The 1968 Stanford conference [25] mainly devoted to 2D boundary layers was a strong reference point in which many models based on eddy viscosity concept, as described before, were applied to various 2D flows and their predictions compared. For low Reynolds number turbulence, the Van Driest correction is used for approaching a wall [26]. It was the time of emerging CFD computational codes like the well-known Patankar and Spalding method [27] judiciously conceived to perform numerical solutions of parabolic expanding flows like jets, wakes and boundary layers. Its efficiency came from the use of the normalized stream function as a variable and the tridiagonal algorithm as a solver in a finite volume framework.

One-equation RANS models
The algebraic EVM tuned for each different type of flow suffered from a lack of generality, the mixing length had to be prescribed empirically in each case. So the route was open to introduce transport modeling. The first step is to model the transport equation for the mean kinetic energy of turbulence k, allowing to account for turbulent diffusion of energy, while the turbulence viscosity was obtained from the Prandtl-Kolmogorov hypothesis [28] ν t = C µ l k. But as the characteristic length scale l was still prescribed empirically, this model was soon abandoned for the benefit of two-equation models. However, another one-equation model was singled out, the Bradshaw model [29] which used a transport equation for the shear stress itself and was dedicated essentially to boundary layer problems.

Two-equation RANS modeling
The main idea underpinning the two-equation models is the need to get rid of empirical scale specification and to provide a general means to get the characteristic length scale of turbulence via an additional transport equation. For computer solvers, the Patankar and Spalding procedure [27] could be easily extended to this kind of models in the case of plane or axisymmetric parabolic flows with the boundary layer type approximations. Using the reduced stream function ω as the variable, the equations to solve take the general form with a marching procedure in X . In the case of 2D or 3D recirculating flows, the equations to solve become elliptic and they usually read where Φ is any scalar or tensorial fluid property. At that time, for calculating recirculating flows, an elliptic solver was initiated in the group of Spalding [30] based on finite volume discretization with staggered grids and the SIMPLE algorithm (semi-implicit method for pressure-linked equations) to solve pressure in incompressible flows and widely known as TEACH code (teaching elliptic axisymmetric characteristics heuristically). The technique has then been further developed by Patankar [31] and a non-staggered grid version also exists with colocated arrangement. Most of the applications in RANS modeling are solved using the finite volume method because of its conservation properties and robustness for solving turbulence transport equations with dominating source terms, even if the method is of first-order precision only. However, some research works have used finite difference expansions for higher-order methods or finite element methods for complex geometry applications.

The k-modeling
This is probably the most widespread and used turbulence model in practice. In its standard original form coming from Jones and Launder [32], the modeled transport equation of the turbulent kinetic energy at high turbulence Reynolds number reads ∂k ∂t where the terms appearing in the right-hand side of this equation are identified as the processes of production P , dissipation-rate and turbulent diffusion J k . The modeled companion transport equation of the dissipation rate reads where J is the turbulent diffusion, c 1 and c 2 are constant coefficients. The turbulent Reynolds stress tensor is then computed as Still based on the turbulence viscosity concept, this model can be numerically solved efficiently using generally finite volume techniques. The characteristic length scale of turbulence defined as l = k 3/2 / is obtained from the modeled transport equation of the energy dissipation rate which had been studied in the pioneering research work at Los Alamos laboratory [33]. The kmodel was also the starting point for some variant models (see hereafter) such as the k-ω model developed by Wilcox [34] using the characteristic frequency ω = /(c µ k) as well as the well-known shear-stress transport (SST) model developed by Menter [35] that is a combination of k-and k-ω models.

Some improvements to the k-model
Numerical predictions of various flows had exhibited some lack of universality like in the wellknown round jet anomaly and the spreading rate of a wake. To improve the round jet prediction without altering plane jets, several proposals have been made by adding extra terms in the dissipation-rate equation [36] or by sensitizing the numerical coefficients to the second and third invariants of anisotropy defined as A 2 = a ij a j i and A 3 = a ij a j k a ki , respectively (Einstein's summation convention is used), where a ij = (τ ij − 2/3kδ ij )/k. In the search for more universality, or more practical use, many variants have been introduced.

Additional term to improve detached flows
Another important correction to single out was needed in detached flows in which the length scale predicted from the model was overestimated, implying discrepancies in associated heat transfer calculations. Some noteworthy modifications to the basic k-model were made for this purpose by introducing additional source terms in the dissipation rate -equation [37].

Other two-equation closures
The wall boundary conditions for the dissipation-rate corresponds to a constant value in the immediate vicinity of the wall in the case where the detailed sublayer is fully resolved. However, for coarse grid resolution, this procedure is no longer possible. A simple approach is then to consider an equilibrium hypothesis between the production and the dissipation-rate (P ≈ ) when approaching the log law near the wall. These practical difficulties had led researchers to develop alternate models using other quantities. Usually, the transport equation for the turbulent energy k is considered in addition to a transport equation for a variable z defined in the general case by z = k m l n , the usual k-model is recovered for (m = 3/2, n = −1) but other combinations may be easier to solve numerically. The length scale itself (m = 0, n = 1) or combined as (m = 1, n = 1) has been proposed in the two-equation k-kl Rotta's model [38,39] and later by Rodi and Spalding [40], or in the two-equation k-l model by Smith [41], or (m = 1, n = −2) in [42], a characteristic frequency (m = 1/2, n = −1) in the k-ω model with ω = /(c µ k) by Wilcox [4,34] used in many industrial flow predictions among others [4]. The series of models [35] proposed by Menter combines the advantages of k-ω and k-models, the SST model widely used in industries accounts in a simplified manner, for the influence of shear stress transport. Other turbulence models have been devised replacing the length scale in the turbulence equations by more complex quantities, such as for instance the one developed by Lin and Wolfshtein [43] using a tensorial volume of turbulence.

Non-linear k-models
Non-linear k-models introduce additional terms in the constitutive relation of the Reynolds stresses in which the gradient term with the eddy viscosity coefficient becomes the first term of an extended development. These models were mainly developed for rotating and corner flows and can be derived from the general formulation as where Ω ij denotes the vorticity tensor and a i are coefficients, eventually of functions of other invariants. In its most elaborate form developed in the Manchester group [44], the model, now free from the pure gradient hypothesis, allows to consider more complex geometries including, for instance, the effect of streamline curvature. This class of models, also studied in [45], in the form k-l or k-ε, is valuable to predict secondary flows in non-circular ducts.

Wall treatment
Special wall treatments were necessary to reduce the number of discretization grid points at the time of low capacity computers. The original wall functions treatment detailed in [27] is based on the hypothesis of constant flux with log-law of the wall. More elaborate treatments have been introduced subsequently, in particular Craft et al. [46] used an analytical method based on integration of simplified mean flow and energy equations in the control volumes adjacent to the wall, more general flows can thus be tackled and extended to roughness and heat transfer [47]. Low Reynolds number k-models, like for instance [48], allow accounting for the influence of molecular viscosity and in particular to describe the details of the near-wall region including the viscous sublayer and the buffer layer joining the logarithmic zone. Useful comparisons of performances of various two-equation models are given in [49] and an up-to-date review of EVM is given by Hanjalic and Launder [50].

Advanced RANS modeling
Several weaknesses appeared in two-equation modeling and in particular in axisymmetric flows, in wake flows, in adverse pressure gradient boundary layers, separated flows, complex geometries, rotating flows, and others, when the physics of energy transfer between the different components of the velocity plays a crucial role in the determination of the flow. It soon appeared that a full account of the effect of these complexities required a more advanced description of the turbulence field. In particular the development of SMC that considers the transport equations of the Reynolds stresses and their dissipation rate by solving the statistical equations of all the components of the Reynolds stress tensor allowed gaining more generality. An overview of the numerous advanced closure methods can be gained through several references [2,3,[51][52][53][54][55][56][57][58]. The basis of development of these models makes extensive use of tensor calculus, anisotropy developments such as the Lumley's invariant modeling [59,60] together with fundamental phenomenology of turbulence. Two major families of models have been developed, Reynolds stress models (RSM) using a set of evolution equations and algebraic stress models (ASM) also called explicit algebraic Reynolds stress models (EARSM) that use an algebraic set of equations for the stress components coupled with a two-equation transport model like a k--type model, free of the eddy viscosity hypothesis that is a simplification of RSM. However, despite its high degree of sophistication, SMC did not always guarantee systematic improvements over high-end RANS models and numerical difficulties were sometimes involved in the cases where the numerical procedure was not appropriate. What is gained in universality is sometimes lost in precision for a specific application. An important aspect in modeling turbulent stresses is the realizability constraints ensuring that the modeled stresses are indeed moments of a probability law. There are two main lines of study, the invariant theory approach [61] and the stochastic analysis based on Langevin equations [62].

Reynolds stress models (RSM)
These models are known as second moment transport modeling closures, sometimes called Differential stress models (DSM). Main physical hypotheses needed to close the model are related to the pressure-strain correlations, diffusion terms and the dissipation of the stress components. The modeled transport equation of the turbulent stress τ ij can be written in a synthetic compact form as where in this equation, the different terms P ij , Π ij , J ij , ij appearing in the right-hand side are the production, redistribution, diffusion and dissipation rate, respectively. In this equation, the redistribution term corresponds to the fluctuating pressure-strain correlation and plays a major role in the correct prediction of the flow anisotropy. This term is decomposed into a slow and a rapid contribution that characterize the return to isotropy. The first ideas in this respect were proposed by the Los Alamos group in New Mexico [33,63]. Hanjalic and Launder [64] as well as Launder et al. [65] proposed the well-known set of hypotheses that is considered as the pioneering works in the formulation of the RSM models. Even to this day, their respective work serves as a reference prototype for subsequent developments. The fundamentals of the methodology can be found in Launder's synthetic overview [51]. Indeed, many proposals can be found in the literature, but broadly speaking, most of these models keep the same basic terms as in [64,65] and simply extend the developments by using higher-order approximations. The pressure-strain correlation modeling is a determinant for a successful closure, it includes three contributions: a linear one (rapid term), a non-linear one (slow term) and a wall reflection term (with rapid and slow term counterparts). The emergence of some large computer centers in the 1970s allowed handling the numerical solution of turbulent flows using second-order closures, for instance on the well-known CDC7600 computer (10 up to 36 max Mflops and 512 K 60-bit words). More advanced closures have subsequently been developed, including some refined features. Low Reynolds number versions of RSM models have been developed and [66] can be considered as a prototype. The use of low Reynolds number correction functions were often useful for approaching walls [67,68]. More advanced forms of low Reynolds number models are more complex, see for instance [69]. Besides, quadratic terms were first introduced in modeling the pressure-strain correlation [70][71][72]. A thorough analysis based on invariant theory is developed in [73]. Then, cubic terms give rise to the two-component limit (TCL) model [56] which is compatible with the tendency of turbulence to become 2D near a wall. The no wall-reflection redistribution topographical terms allow discarding explicitly the wall distance in the approximations using non-topographic wall detectors [74][75][76][77]. Extensive applications using these models can be found in [78]. In practice, the numerical solution of the stress tensor equations brings in a new difficulty because, in the absence of eddy viscosity, the momentum equation loses its diffusive dominant form and becomes stiff. The problem can been solved by numerical stabilization practices such as a fourth shifted grid for shear stress discretization and introducing apparent viscosities in the discretized equations [79,80]. This procedure can be extended to non-shifted grids [81]. In order to illustrate the capabilities of these type of methods, Figure 1 shows some application results of different closures made in the Leschziner Imperial College group to the supersonic fin plate junction using the RANS-SST and the RANS-MCL (modified component limit) which is a modification for compressible flow [82] of the TCL closure employing a cubic pressure-strain model and entirely topology free. As a result of interest, these results show that only the secondorder closure is able to reproduce multiple separation/reattachments ahead of the fin. Figure 2 applies successfully the same MCL closure to the complex flow around a 3D afterbody with issuing square jet which shows the advantages of the MCL closure over simpler models as discussed in [83]. More generally, second-order closures are beneficial for dealing with complexities due to geometry or interaction with other phenomena. Applications are numerous and widely represented in the scientific literature. Among them, some specific examples can be found in several references [68,[84][85][86][87][88][89][90][91].

Algebraic stress models (ASM)
First introduced intuitively by Rodi [92] and subsequently developed by Gatski and Speziale [93], Durbin and Pettersson-Reif [94], these ASM models rely on two-equation models. They were intialy developed with the aim to reduce computational costs in comparison with the one required by RSM models that need to solve seven coupled equations. In its simplest form, the ASM can be deduced from stress transport equations using the so-called Rodi hypothesis [92,95] originally written as approximately equivalent in homogeneous turbulence to assuming that the flow anisotropy remains constant along the streamline da ij /dt = 0. Using (8) and (11), the ASM model then consequently reads But despite the simplification of ASM in regards with RSM, it has been found in practice that some problems in the numerical solution of the equations were still acute and even worsen in some particular cases. In the framework of invariant modeling, it is worth noting the use of Rodi's hypothesis by Pope [96] to derive a non-linear viscosity model. The ASM model of Wallin and Johansson [97] is developed in k-ω form including near-wall treatment ensuring realizability of the stress components and particularly dedicated to compressible boundary layers.

Renormalization group theory
Turbulence models were also developed in the framework of the renormalization group (RNG) theory by Yakhot [98] with scale expansions for the Reynolds stress and the source of dissipation terms as a complementary tool in turbulence modeling [99]. These methods allowed suggesting new additional terms in existing models and also to give analytical expressions for the coefficients. The k--RNG model is a well-known example.

Compressible turbulent flows modeling
Several averaging procedures exist [100,101], but there are mainly two approaches, the densityweighted Favre-averaged equations [102] and Reynolds-averaged equations. In strongly compressible flows extra terms are to be included in the model equations [100,103,104]. However, in some cases, when density fluctuations are weak and can be neglected, modeling is greatly simplified. An important change, with respect to CFD, is the type of numerical methods to be used due to the fact that pressure acquires a thermodynamic meaning and is closely tied with the energy equation. In this case, the continuity equation is solved for density and pressure is obtained from an equation of state, but this is not specific to turbulent flows.

Scalar transport modeling
Associated turbulent heat and/or mass transfer are treated using the same principles. In addition to the turbulent scalar fluxes, the scalar variance and its dissipation rate can be modeled separately in the case of passive scalars θ. In its basic form [51,70,105], the model equations to solve closely resemble their dynamic counterpart. More advanced closures [106,107] have also been developed subsequently. The generalized gradient diffusion hypothesis (GGDH) first introduced by Daly and Harlow [33] is often applied in a first approximation to compute the heat flux where c τ θ is a numerical coefficient. DNS still remains an useful tool to validate turbulence models in this field [108,109]. The extension to the case of active scalar transport is also a practical field for many applications, a typical example being the modeling of turbulence subjected to buoyancy [110]. When dealing with turbulent fluxes of a transported scalar, realizability constraints have also to be considered, like in dynamical problems. For instance the Langevin equation approach [111] may give straightforward results.

Elliptic relaxation models
The concept of elliptic relaxation was proposed by Durbin [112] to model non-homogeneous turbulence by means of an additional elliptic equation that is coupled with the k-model. This method was then extended to RSM models with the modeling of the pressure-strain correlation term [113]. Complex geometries can be cleanly treated at the price of higher computational effort. A simplified approach consisting of the elliptic-blending second-order closure has been then developed to alleviate the computational requirements [114,115]. This method entails the solution of a single one-elliptic equation for a blending function which controls transition between near-wall and far-from-the-wall pressure strain model (see [58]).

Multiple scale models
The fact that usual RANS closures are devised as single-scale closures is justified by the Kolmogorov theory of universal cascade determined by the dissipation-rate only. Departures from this hypothesis are however expected when the turbulence is out of equilibrium. A first attempt to distinguish fine dissipative scales from energetic scales was proposed by Schiestel [116] and later developed in [117] and [118]. The analytic work made in [118,119] explains the link with spectral closures. The final equations are then obtained by partial integration of the spectral spectrum with m = 1, n slices (n = 2 or 3 in practice) in the wave number ranges [κ m−1 , κ m ] and have the same overall structure as standard models. For instance, the transport equation for the partial turbulent stress τ (m) ij reads [119] ∂τ (m) where in this equation, the different terms ij appearing on the righthand side are the production, in and out transfer fluxes, redistribution, diffusion and dissipation rate, respectively. This equation can be solved in conjunction with transport equations for the transfer fluxes F (m) , using the same numerical procedures as for usual statistical models. From the present section to subsequent sections including statistical multiple scale models, statistical spectral models, LES and hybrid subfilter models, an account of the differing eddy scales is possible as opposed to single-scale closures, and these approaches can be viewed as various multiscale and multiresolution methods [120].

What about URANS?
According to the acronym, URANS are unsteady solutions of RANS models, thus they produce time varying statistical fields. Interesting results can be obtained in flows such as separated layers, wakes, rotating flows, free convection [121] in which some periodic behaviour may happen. But, in some cases the unsteadiness is irregular and non-periodic, looking like a macro-simulation of turbulent large eddies, in this case, although the resulting calculation looks physically realistic, it is an unsuited application of URANS in contradiction to its appellation. It can be viewed as a large-scale eddy simulation without any reference to the grid, and thus the physical interpretation becomes difficult. For instance, true URANS allowed simulating self-sustained oscillations of a turbulent plane jet issuing into a rectangular cavity [122] as well as the vortex shedding in solid rocket motors [123]. The conceptual problem in URANS is that the separation of scales is not always clearly possible. For instance vortex shedding in pure URANS may break down into true turbulence (LES).

Spectral turbulence models
In conjunction with the development of turbulence models in the physical space for users involved in engineering or industrial applications, spectral turbulence models have also been devised to study turbulent flows in the laboratory from a fundamental point of view with emphasis on the physical aspects of turbulence. All the methods of closure considered so far were dealing with one-point statistics (except multiscale models, sometimes referred to as 1.5-point closures). Two-point closures have been developed in a more theoretical framework mainly in Fourier space giving rise to spectral modeling. Due to the increased complexities, most of these models are originally limited to homogeneous and isotropic turbulence [12], though these limitations have been removed now. These spectral theories initially developed by Jeandel et al. [13] may be used to deduce simpler one-point closures from spectral integration and spherical averaging showing a hierarchy in levels of description. Extension to anisotropic turbulence was considered in particular by Cambon et al. [124].

Early spectral theories
These theories limited to homogeneous isotropic turbulence have been introduced long ago in order to approximate the spectral flux due to the inertial cascade and thus calculate the mean energy spectrum compatible with the Kolmogorov spectrum in equilibrium flows. These models are fully described in [125,126] and rely on the equation of the two-point correlation tensor φ ij = 〈 u i (x, κ) u j (x, κ )〉δ(κ + κ ), where u i denotes the Fourier transform of the fluctuating velocity u i .

The EDQNM model
Among these spectral models, one of the most popular is the eddy damped quasi-normal Markovianized (EDQNM) model [127] which focusses on the non-linear inertial terms of closure. Extension to non-homogeneous turbulence has been also considered [128]. This approach embodies many important properties of the inertial cascade interactions that justifies its interest in fundamental laboratory studies. But its analytical complexities prevent more extended practical applications.

Models with hidden parameters
Stimulated by the peculiar properties of rotating turbulence, two-point closures were further studied and developed by introducing two anisotropy tensors for the polarization anisotropy and for the directional anisotropy in wave vector space [129][130][131] in order to figure out the complex tensorial properties of the fluctuating turbulent field. A similar concept was also developed in physical space [132]. Then, integrated in Fourier space to get a one-point closure, these models known as structure-based models, introduce hidden parameters linked to the two types of anisotropies mentioned earlier.

Later computer developments: direct and large eddy simulations
In this section we leave the purely statistical approach considered previously to perform numerical simulation of a realization of the fluctuating turbulent flow either completely (DNS) or partially (LES).

Computational resources in fine grid simulations
Fine grid numerical simulations of turbulence require high computational resources. They can be a priori estimated in the following way. For a DNS, the calculation must solve the smallest flow eddies down to the Kolmogorov scale η K . But also the dimensions of the computational domain must be large enough to comprise the largest turbulence scales. Taking into account these constraints, in the case of homogeneous turbulence in a box domain, the necessary number of grid points is found to be N 1 N 1 N 3 = 64R 9/4 t , where R t is the turbulence Reynolds number (usually of order Re/10), the computational time being proportional to T ∝ R 11/4 t . In the case of LES, the number of necessary grid points is obviously reduced in the ratio (η K /∆) 3 where ∆ is the grid step. When considering real shear flows like boundary layers these estimates need to be revised [133,134]. These requirements can be checked against the evolution of computer power [10]. From a practical point of view, the numerical methods suited for such calculation are somehow different from the ones used in RANS modeling. RANS modeling with complex transport equations and stiff source terms achieve stability with relatively low-order numerical schemes. On the contrary, for fine grid simulations using either the pure Navier-Stokes equations (DNS) or viscosity hypothesis (LES), the use of higher-order methods is necessary in order to get precision, avoid numerical viscosity and dissipation.

Developments in LES
The very first LES calculation goes back to the 1970s with the pioneering works of Orszag and Patterson [135] for homogeneous flows and Deardorff [136] for channel flow followed by the work of Moin et al. [137] as a first milestone in LES development. Since that time, substantial progress has been made by several groups summarised in [138,139]. Some thoughts about the conceptual foundations in LES are discussed in [140].

The original Smagorinsky model and the dynamic Smagorinsky model
The aim of LES is to spare computer resources while simulating the non-universal large scales as much as possible in order that only the most universal smaller scales corresponding to the end of the energy spectrum are modeled. Physically, this means that the cutoff wave number κ c is placed in the region of the Kolmogorov law given by E (κ) = C κ 2/3 κ −5/3 where C κ stands for the Kolmogorov constant. So, simpler models are usually sufficient for a good account of subfilter turbulence. The first to be considered and largely used is the well-known Smagorinsky model [141,142] which looks like the mixing length hypothesis but the length scale being given by the grid step. The SGS turbulent stress tensor is then computed using the Boussinesq hypothesis as where the modeling of the subgrid turbulent eddy viscosity inspired from the mixing length hypothesis reads In (17), C s is the Smagorinsky coefficient that takes on a constant value [141] in the standard Smagorinsky model. Germano et al. [143] developed the dynamic Smagorinsky version of the model where the coefficient C s is evaluated locally and dynamically in time and space by introducing a superfilter used to estimate the Smagorinsky "constant" directly from the simulated flow, allowing a far better universality. This approach was improved by Lilly [144] and often used for highly resolved LES, like for instance in [145].

The structure-function models
The "structure-function" models introduced by Lesieur's team find their basic foundation in the spectral space [6,139]. These models indeed extend in physical space the model of spectral viscosity [146] in Fourier space by the same research group. They have been developed in several variants [6] with mainly the "selective structure function" model in which the viscosity is reduced in near 2D turbulence and the "filtered structure function" model in which the large scales are filtered out before computing viscosity to overcome the too dissipative behaviour of the pure structure function model. This model was applied to simulate a large variety of flows [6] such as for instance rotating channel flows [147] and mixing layers [148].

The Bardina model
The scale-similarity Bardina model [149] bases its formalism using twice filtering on the idea that local interactions near the cutoff are dominant. In practice, the model proved to be not very dissipative and it is mainly useful in combination with the Smagorinsky model. Indeed, while viscosity-based models were efficient to account for the energetic dissipative effects in scales interactions, the structural aspects of subfilter turbulence is better represented by scale-similarity models.

Subfilter-scale transport models
Looking for more advanced turbulence description, more complex LES models have been developed using transport equations. Yoshizawa and Horiuti [150] proposed a subfilter model using a transport equation for the SGS turbulent energy k sgs where the length scale of turbulence is given by the grid step size ∆. In fact, Deardorff's work [151] involved with geophysical flows, had the merit of considering early the transport equation of the SGS stress tensor (τ ij ) sgs in its full formulation, still with the use of ∆. These approaches however implicitly assume approximate spectral equilibrium between the production and dissipation rate because is only deduced from ∆, and has thus opened new routes in LES subgrid-scale modeling. All these types of models have their advantages and drawbacks that can be explained on the basis of their statistical properties [152].

The role of DNS
DNS has long been considered as a tool for analyzing detailed laboratory turbulent flows as a substitute to expensive experiments. In this sense, it contributes to fundamental knowledge of turbulent flows. Its use in practical real life situations such as engineering and environment is more difficult, considering the additional geometric and physical complexities often present and not always necessary with respect to the answers that are sought. DNS development began with the pioneering works of the Stanford team [153,154] with the investigation of the fully developed turbulent channel flow. The statistics of the turbulent energy, dissipation rate and correlations of the fluctuating velocities were worked out to determine the flow characteristics. DNS data are often used for the validation of turbulence models to this day. This research field is growing fast despite the difficulties and limitations imposed by the need for increased power of the super-computers as well as the development of computational techniques including vectorisation and parallelisation. But even with modern super-computers, the applicability of DNS still remains limited to flows with relatively low or moderate Reynolds numbers [155][156][157]. Besides investigating the detailed turbulence field in fundamental laboratory turbulent flows, another important application of DNS is in providing reference tests and benchmarks for evaluating simpler models. In many cases, DNS can be viewed as a complement or a substitute to experimental invertigations, the hot wire being replaced by the discretization point. DNS are now a tool in turbulence research allowing to devise novel numerical experiments [134] that are not possible in the real laboratory.

Developments in hybrid RANS/LES simulations
At the end of the 20th century, the concept of hybrid RANS/LES began to take shape. A detailed up-to-date review of hybrid models can be found in [9][10][11] and examples for flow applications in practice in [158,159].

Two types of models: hybrid zonal and non-zonal models
Zonal methods split the computational domain into several subdomains in which different models are applied with the hard problem of control of boundaries. The question that is raised is then to match the different flow regions by means of artificial turbulent fluctuations [160]. Non-zonal methods may embody an automatic RANS-LES switch parameter or use a progressive change in the model so that seamless coupling is achieved. An attempt to unify these respective formalisms which encompass RANS and LES can be found in [161]. The important issue of noncommutation errors mentioned above and coming from variable filtering in simulations may be aggravated in hybrid models specially at the internal boundaries in zonal approaches. This problem has been analyzed by Hamba [162] on the basis of DNS comparisons showing that these non-commutation errors increase near the interface. Another aspect of this problem is the log-law mismatch in near-wall flows for which Hamba [163] justifies additional filtering from approximation of commutation terms.

Early hybrid modeling
An early contribution to hybrid models is the VLES proposed by Speziale [164], which combines RANS and DNS by damping the turbulent stresses in regions where the grid step is finer so that the calculation runs between RANS and DNS depending on the grid spacing. In VLES the unresolved scales region embodies the energy containing eddies and unsteady closures are necessary.

The detached eddy simulation (DES)
One of the most popular hybrid models, widely applied in practical flow calculation and especially in aeronautics, is the detached eddy simulation (DES) developed by Spalart [53,165], Spalart et al. [166]. This approach makes use of the Spalart-Allmaras RANS model using one-equation transport of turbulent viscosity in the wall region in which the turbulence length scale is given by the wall distance d w and is replaced by the grid step far from the wall where ∆ = max(∆ 1 , ∆ 2 , ∆ 3 ) and C DES is a constant coefficient. Since the same model is used in both zones, their junction is continuous and hence the pure DES approach is no longer zonal. This method has been extended to two-equation models using the SST k-ω model of Menter [35] with some adaptations introduced in the sink term of the transport equation for the subgrid turbulence energy. Another extension named delayed detached eddy simulation (DDES) [167] uses a parameter to delay the LES function in boundary layers, including the molecular and turbulent viscosity information into the switching mechanism. Then, the improved delayed detached eddy simulation (IDDES) brings improved wall-modeling capabilities with also an SST-IDDES variant. For purposes of illustration, Figure 3 shows an interesting application of the transonic flow over an axisymmetric bump with shock-induced separation described in [168].
In order to get the precise description of the shock region a three zone RANS-DNS-IDDES calculation has been performed using DNS in the shock region, RANS upstream and IDDES downstream. Figure 4 displays an instantaneous view of the fine turbulent eddies in the shock vicinity.

The partially integrated transport modeling (PITM)
Inspired from the multiple scale concept developed in spectral space, the partially integrated transport modeling (PITM) method developed by Schiestel and Dejoan [169] for the subfilterscale k sfs -model and by Chaouat and Schiestel [16,[170][171][172] for the (τ ij ) sfs -model can be applied to almost any existing RANS model using the dissipation-rate transport equation. The present modeling approach of the PITM method finds its basic foundations in the spectral space by considering the transport equation of the the two-point fluctuating velocity correlations in the physical space [161,173]. This method was initially formulated in the case of anisotropic homogeneous flows and then extended to non-homogeneous flow considering the concept of the tangent homogeneous space [17,161]. Using the variational calculus, it has been demonstrated that the properties of the model established in homogeneous turbulence can be extended mutatis mutandis to the case of non-homogeneous turbulence [174]. Taking its Fourier transform and averaging over spherical shells, the resulting spectral equation reads [119,161] ∂ϕ ij (X , κ) where the different terms appearing in the right-hand side of this equation are respectively the production, transfer, redistribution, diffusion and dissipation contributions, acting in the spectral space associated with the scalar wave number κ (modulus of the wave vector). Spectral splitting with partial integration over [κ c , ∞[ gives rise to the equation for subfilter-scale stresses in the physical space after some algebra. The dissipation-rate transport equation derived using the multiscale technique [119,169] looks like (9) but the c 2sfs coefficient is variable and acts as a dynamical parameter to control the relative amount of subfilter energy. More precisely, this coefficient is now a function of the grid-step size of the mesh ∆ ratio to the turbulence length-scale l = k 3/2 / , so that c 2sfs = c 2sfs (η c ), where η c = κ c l = πl /∆ [170][171][172]. From a physical point of view, the dissipation rate interpreted here as a flux of energy that is transferred from the large scale to the small scale remains the same as the one returned by (9) because it is independent of the cutoff-wave number ∂ /∂κ c = 0 in the equilibrium range.

Two-equation subfilter model
The two-equation subfilter method embodies a subfilter-scale EVM that has been first developed for applications in standard turbulent flows and engineering flows with ease of calculation. It is simple to use in the framework of two-equation models and allows combining the advantages of both RANS and LES in a practical manner. This model has simulated fairly well, for instance, the turbulent pulsed channel flows [169,175] and the mixing of turbulent flow streams involving differing scales [175].

Subfilter stress transport model
This subfilter-scale stress model is more advanced in terms of modeling than EVM because it relies on the transport equations for the subfilter stress (τ ij ) sfs and the dissipation-rate . The where the terms appearing in the right-hand side of this equation are identified as the subfilter production, redistribution, diffusion and dissipation, respectively. The transport equation for the dissipation-rate is still given by (20) but the diffusion term assumes now a tensorial diffusivity hypothesis. The tensorial dissipation rate is approached by ij = 2/3 δ ij . This model allows describing more accurately the physical mechanisms of the turbulence processes. In particular, it encompasses the pressure-strain correlation term that redistributes the energy among the stress components and the anisotropy of the dissipation in stress components for reproducing the flow anisotropy [170][171][172]. This model was applied to a large variety of both internal and external flows with success, accompanied by a drastic reduction of the grid-points and computational time in comparison with standard LES models, thus showing promising perspectives. Various applications were tackled, such as flow in a plane channel with appreciable fluid injection through a permeable wall corresponding to the propellant burning in solid rocket motors [170], rotating channel flows encountered in turbomachinery [176], flow over periodic hills with separation and reattachment of the boundary layer at several Reynolds numbers as shown in Figure 5 [177,178] corresponding to the experiment reference [179], airfoil flows [180], flow in small axisymmetric contraction [181]. The PITM method has been recently extended to the turbulent transfer of a passive scalar including transport of variance and its dissipation rate [182].

The partially-averaged Navier-Stokes (PANS)
The partially-averaged Navier-Stokes (PANS) model was introduced by Girimaji [183] and is based on the transport equations for the SGS turbulent energy k sgs and its dissipation rate. The first applications can be found in [184]. In this method, contrarily to the PITM method, the ratio of subgrid modeled energy to the total energy f k = k sgs /k is rather imposed at a constant arbitrary value. This hypothesis once inserted into the turbulence model, the resulting transport equations formally look almost similar to the PITM ones in spite of a totally different approach. Some applications handled by PANS have been summarized in the recent paper [185] such as for instance, the turbulent flows past a square and circular cylinder [186], the flow around a rudimentary landing gear [187], swirling confined flows [188] as well as the flows over periodic hills [189].

The scale adaptive simulation (SAS)
The scale adaptive simulation (SAS) was derived by Menter et al. [190,191] from the concept of the k-kl Rotta model. It adds another scale L νK using the second derivative of the velocity field known as the Von Kármán length scale defined as to the traditional input of the velocity gradient tensor, K being the von Kármán constant. Contrarily to previous models, there is no explicit dependency on the grid spatial resolution [190,191].
There are essentially two variants of the approach, the KSKL-SAS (K-square-root-kL-SAS) two equation model and the SST-SAS (shear-stress transport-SAS) model [190]. A typical example of application, among many others, given in [191] and shown in Figure 6 allows to predict the intensive mixing caused by the turbulence generated in the unstable regime in a combustion chamber. More recently, Menter et al. [192] demonstrated that the stress-blended eddy simulation (SBES) approach based on the blending of existing RANS and LES models using linearly-weighted stress components is optimal for applications with a mix of boundary layers and free shear flows. The case of the mixing layer displayed in Figure 7 shows how easily the SBES model develops finegrained turbulence eddies that are impossible to obtain with more standard approaches.

Numerical methods for the simulation of turbulent flows
Different types of numerical methods in fluids have been applied for the modeling and simulation of turbulent flows that are essentially the finite volume method, the finite difference method, and the finite element method [193]. Finite volumes and finite differences are the most popular in practice while finite element methods using weight functions are more oriented towards mathematical properties of the numerical method. If the finite differences method is essentially used in conjunction with structured grids in Cartesian or curvilinear coordinates [194], the finite volume method is now more and more used with unstructured grids, as obviously is the finite element method. The numerical schemes differ from RANS to LES both in time and space. In RANS, low-order upwind schemes are often applied considering that the flow is often steady. Mean while in LES, high-order centered schemes are emphasised to accurately capture the unsteady regime of the flow and the insight into the evolving turbulent structures [195,196]. Besides explicit highorder schemes in space, mention has to be made of the Hermitian and compact schemes [197] that allow precision from implicit additional relations. Among the advanced methods developed for the Navier-Stokes equations [196][197][198][199][200][201][202][203][204], various specific techniques have been adapted. In all methodologies including transport equations of the turbulent stresses, a special numerical treatment is necessary because of the mathematical complexity of solving these equations, which are strongly coupled leading to a lack of robustness of the numerical scheme, both in cases of structured meshes [205][206][207] and unstructured meshes [208,209]. The use of spectral numerical methods (not to be confused with spectral closures) known for their high precision [199,200] is useful for DNS in relatively simple geometries, but their extensions to more complex models and geometries is difficult. Spectral methods and spectral element methods are thus specialized for   [192].) applications in LES and DNS flows in simple geometries. As regards time integration, various explicit and implicit time discretization schemes are used, multilevel methods may have advantages at the price of more complexity. Good precision of time advancement schemes is crucial in DNS and LES calculations for which small spurious errors may compromise the long duration simulation. Generally, it is advisable to avoid too dissipative schemes which smoothen out some frequencies and introduce errors. Intensive use of vectorization and parallelization programming allowed reducing the computational time required for the simulation. In particular, CFD codes are now often optimized with computational techniques such as the message passing interface (MPI).

Validation of turbulence models
In turbulent flow numerical predictions, the use of appropriate numerical methods is mandatory, but the inherent uncertainties of physical modeling (full statistics or subfilter statistics) has also to be considered and appreciated. In former RANS calculations, the experimental reference was the only reference to check the performance of a model, but as mentioned above, it is now complemented by DNS. The requirements for successful simulation are thus twofold. Practical examples are innumerable, benchmarks of well documented turbulent flows remain an invaluable reference. Isotropic homogeneous turbulence and the turbulent plane channel flow have been very standard test cases. Among many others, we may cite here the flow over periodic hills illustrated in Figure 5 studied in [179] and which led to extended testing of models [210]. On the pure numerical point of view LES-type methods can be evaluated on "a priori" tests using data from a previous DNS and "a posteriori" tests using the results of the LES calculation [211].

Concluding remarks and future prospects
The last 50 years of CFD turbulence development in the research area have seen a huge variety of turbulent closures and associated numerical techniques roughly evolving from less statistical to more simulational, following the constant increase in scientific computing power. LES has now arrived at maturity. But in spite of the huge development of this computing power, many engineering and environmental turbulent flows remain still out of scope of DNS full simulation and even fine LES. The older and simpler RANS methods remain useful tools and can still be recommended as the starting point (and sometimes also the finishing point) for engineering simulations [159]. Higher-level closures allow to account for superimposed phenomena and complex flows [3,212] and even nowadays, RANS model developments are still needed. For this reason, a large array of methods with different levels of description is necessary. Indeed, the choice among the different methodologies available strongly depends on the physical problem considered and the type of answers which are expected. In particular the hybrid range of techniques, considering its flexibility, can encompass many practical situations [10,11,213,214] and are a good alternative, with the accuracy of LES and the speed of RANS. Some problems like intermittency would need more investigation. Many of the problems discussed are encountered in aerospace industry in which CFD plays an increasingly crucial role [7] and in which various types of methods have their place for different uses [8,58,215,216]. The question is also open whether modern data-driven techniques may be able to find optimal models in a user-defined sense [217,218].

Conflicts of interest
The authors declare no competing financial interest.

Dedication
The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.