chapter. Turbulence Modelling 第7章 湍流建模[cfd-0008]

The first approach for the approximate treatment of turbulent flows was presented by Reynolds in 1895. The methodology is based on the decomposition of the flow variables into a mean and a fluctuating part. The governing equations (7.1) are then solved for the mean values, which are the most interesting for engineering applications. Thus, considering first incompressible flows, the velocity components and the pressure in Eq. (7.1) are substituted by [13]

where the mean value is denoted by an overbar and the turbulent fluctuations by a prime. The mean values are obtained by an averaging procedure. There are three different forms of the Reynolds averaging:

1. Time averaging -- appropriate for stationary turbulence (statistically steady turbulence)

As a consequence, the mean value \(\bar{v}_i\) does not vary in time, but only in space. The situation is sketched in Fig. 7.2. In practice, \(T\to\infty\) means that the time interval \(T\) should be large as compared to the typical time-scale of the turbulent fluctuations.

2. Spatial averaging -- appropriate for homogeneous turbulence

with \(\Omega\) being a control volume. In this case, \(\bar{v}_i\) is uniform in space, but it is allowed to vary in time.

3. Ensemble averaging -- appropriate for general turbulence

Here, the mean value \(\bar{v}_i\) still remains a function of time and of space coordinates.

For all three approaches, the average of the fluctuating part is zero, i.e., \(\overline{v'_i}=0\). However, it can be easily seen that \(\overline{v'_iv'_i}\neq 0\). The same is true for \(\overline{v'_iv'_j}\), if both turbulent velocity components are correlated.

In cases where the turbulent flow is both stationary and homogeneous, all three averaging forms are equivalent. This is called the ergodic hypothesis.

In cases where the density is not constant, it is advisable to apply the density (mass) weighted or Favre decomposition [14], [15] to certain quantities in Eq. (7.1) instead of Reynolds averaging. Otherwise, the averaged governing equations would become considerably more complicated due to additional correlations involving density fluctuations. The most convenient way is to employ Reynolds averaging for density and pressure, and Favre averaging for other variables such as velocity, internal energy, enthalpy and temperature. Favre averaged quantities, for example the velocity components, are obtained from the relation [14], [15]

where \(\bar{\rho}\) denotes the Reynolds-averaged density. Hence, the Favre decomposition reads

where \(\tilde{v}_i\) represents the mean value and \(v''_i\) the fluctuating part of the velocity \(v_i\). Again, the average of the fluctuating part is zero, i.e., \(\widetilde{v''_i}=0\). Furthermore, the average of the product of two fluctuating quantities is not zero, if the quantities are correlated. Hence, for example, \(\widetilde{v''_iv''_i}\neq 0\) and in general \(\widetilde{v''_iv''_j}\neq 0\).

The following relationships can be derived for a mix between Favre and Reynolds averaging

These relations will be utilised in later subsections.

If we apply either the time averaging Eq. (7.8) or the ensemble averaging Eq. (7.10) to the incompressible Navier-Stokes equations (7.6), we obtain the following relations for the mass and momentum conservation

These are known as the Reynolds-Averaged Navier-Stokes equations (RANS). The equations (7.14) are formally identical to the Navier-Stokes equations (7.6) with the exception of the additional term

which constitutes the so-called Reynolds-stress tensor. It represents the transfer of momentum due to turbulent fluctuations. The laminar viscous stresses are evaluated according to Eqs. (7.2) and (7.3) using Reynolds-averaged velocity components, i.e.,

The Reynolds-stress tensor consists in 3D of the nine components

However, since \(v'_i\) and \(v'_j\) in the correlations can be interchanged, the Reynolds-stress tensor contains only six independent components. The sum of the normal stresses divided by density defines the turbulent kinetic energy, i.e.,

As we can see, the fundamental problem of turbulence modelling based on the Reynolds-averaged Navier-Stokes equations is to find six additional relations in order to close the equations (7.14). We shall introduce the basic methodologies in the Subsections 7.1.5--7.1.7.

In turbulence modelling, it is quite common to assume that Morkovin's hypothesis [16] is valid. It states that the turbulent structure of a boundary layer is not notably influenced by density fluctuations if \(\rho'\ll\bar{\rho}\). This is generally true for wall-bounded flows up to a Mach number of about five. However, in the case of hypersonic flows or for compressible free shear layers, density fluctuations have to be taken into account. The same holds also for flows with combustion or with significant heat transfer.

Application of the Reynolds averaging (Eq. (7.8) or (7.10)) to density and pressure, and of the Favre averaging Eq. (7.11) to the remaining flow variables in the compressible Navier-Stokes equations (7.1) yields [17]

These are the Favre- and Reynolds-Averaged Navier-Stokes equations. Similarly to the Reynolds averaging, the viscous stress tensor in the momentum (and energy) equation is extended by the Favre-averaged Reynolds-stress tensor, i.e.,

Its form is similar to Eq. (7.17) with Favre instead of Reynolds averaging. The components of the laminar (molecular) viscous stress tensor \(\bar{\tau}_{ij}\) are evaluated by Eq. (7.2) using Favre-averaged velocity components.

If we employ the definition of the Favre-averaged turbulent kinetic energy, i.e.,

we can express the total energy in Eq. (7.19) as

The total enthalpy is defined as

The individual parts of the Favre- and Reynolds-averaged Navier-Stokes equations (7.19) have the following physical meaning [17]:

• \(\dfrac{\partial}{\partial x_j}\left(k\dfrac{\partial\tilde{T}}{\partial x_j}\right)\) - molecular diffusion of heat
• \(\dfrac{\partial}{\partial x_j}\left(\bar{\rho}\widetilde{v''_jh''}\right)\) - turbulent transport of heat
• \(\dfrac{\partial}{\partial x_j}\left(\widetilde{\tau_{ij}v''_i}\right)\) - molecular diffusion of \(\tilde{K}\)
• \(\dfrac{\partial}{\partial x_j}\left(\bar{\rho}\widetilde{v''_jK}\right)\) - turbulent transport of \(\tilde{K}\)
• \(\dfrac{\partial}{\partial x_j}\left(\tilde{v}_i\bar{\tau}_{ij}\right)\) - work done by the molecular stresses
• \(\dfrac{\partial}{\partial x_j}\left(\tilde{v}_i\tau^F_{ij}\right)\) - work done by the Favre-averaged Reynolds stresses

The molecular diffusion and turbulent transport of \(\tilde{K}\) are very often neglected. This is a valid approximation for transonic and supersonic flows. In order to close the Favre- and Reynolds-averaged equations (7.19), we also have to supply six components of the Favre-averaged Reynolds-stress tensor (Eq. (7.20)) and three components of the turbulent heat-flux vector. We shall discuss the three basic approaches in the next subsections.

One of the most significant contributions to turbulence modelling was presented in 1877 by Boussinesq [11], [12]. His idea is based on the observation that the momentum transfer in a turbulent flow is dominated by the mixing caused by large energetic turbulent eddies. The Boussinesq hypothesis assumes that the turbulent shear stress depends linearly on the mean rate of strain, as in a laminar flow. The proportionality factor is the eddy viscosity. The Boussinesq hypothesis for Reynolds averaged incompressible flow (Eq. (7.14)) can be written as

where \(\overline{S}_{ij}\) denotes the Reynolds-averaged strain-rate tensor (Eq. (7.3), cf. also Eq. (7.16)), \(K\) is the turbulent kinetic energy (\(K=(1/2)\overline{v'_iv'_i}\)), and \(\mu_T\) stands for the eddy viscosity. Unlike the molecular viscosity \(\mu\), the eddy viscosity \(\mu_T\) represents no physical characteristic of the fluid, but it is a function of the local flow conditions. Additionally, \(\mu_T\) is also strongly affected by flow history effects.

In the case of the compressible Favre- and Reynolds-averaged Navier-Stokes equations (7.19), the Boussinesq eddy-viscosity hypothesis reads

where \(\tilde{S}_{ij}\) and \(\tilde{K}\) are the Favre-averaged strain rate and turbulent kinetic energy, respectively. Note the similarity to Eq. (7.2). The term \((2/3)\rho K\delta_{ij}\) in Eqs. (7.24) and (7.25) is required in order to obtain the proper trace of \(\tau^R_{ij}\) or \(\tau^F_{ij}\). This means that we must have

in the case of \(\overline{S}_{ii}=0\) (continuity equation) or \(\tilde{S}_{ii}=0\), in order to fulfil the relations Eq. (7.18) or (7.21) for the turbulent kinetic energy. However, the term \((2/3)\rho K\delta_{ij}\) is often neglected, particularly in connection with simpler turbulence models (like algebraic ones).

The approximation, which is commonly used for the modelling of the turbulent heat-flux vector, is based on the classical Reynolds analogy [18]. Hence, we may write

with the turbulent thermal conductivity coefficient \(k_T\) being defined as

In Equation (7.27), \(c_p\) denotes the specific heat coefficient at constant pressure and \(Pr_T\) is the turbulent Prandtl number. The turbulent Prandtl number is in general assumed to be constant over the flow field (\(Pr_T=0.9\) for air).

By applying the eddy-viscosity approach to the Reynolds- (and Favre-) averaged form of the governing equations (2.19) or Eq. (7.1), the dynamic viscosity coefficient \(\mu\) in the viscous stress tensor Eq. (2.15) or Eq. (7.2) is simply replaced by the sum of a laminar and a turbulent component, i.e.,

The laminar viscosity \(\mu_L\) is computed, for example, with the aid of the Sutherland formula (2.30). Furthermore, according to the Reynolds analogy given by Eq. (7.26), the thermal conductivity coefficient \(k\) in Eq. (2.24) or Eq. (7.2) is evaluated as

The eddy-viscosity concept of Boussinesq is, at least from the engineering point of view, very attractive since it requires ``only'' the determination of \(\mu_T\) (the turbulent kinetic energy \(K\) needed for the term \((2/3)\rho K\delta_{ij}\) in Eq. (7.24) or (7.25) is either obtained as a by-product of the turbulence model or is simply omitted). Once we know the eddy viscosity \(\mu_T\), we can easily extend the Navier-Stokes equations (2.19) or (7.1) to the simulation of turbulent flows by introducing averaged flow variables and by adding \(\mu_T\) to the laminar viscosity. Therefore, Boussinesq's approach became the basis for a large variety of first-order turbulence closures. However, there are applications for which the Boussinesq hypothesis is no longer valid (see, e.g., [17] p. 214 or [19] p. 111):

• flows with sudden change of mean strain rate,
• flows with significant streamline curvature,
• flows with rotation and stratification,
• secondary flows in ducts and in turbomachinery,
• flows with boundary layer separation and reattachment.

The limitations of the eddy-viscosity approach are caused by the assumption of equilibrium between the turbulence and the mean strain field, as well as by the independence on system rotation. The results can be notably improved by using appropriate correction terms in the turbulence models [20], [21]. Further increased accuracy of predictions can be achieved through the application of non-linear eddy-viscosity models which are described next.

In order to remove the restrictions imposed by the assumption of equilibrium between the turbulence and the mean strain rate, Lumley [22], [23] proposed to extend the linear Boussinesq approach by higher-order products of strain and rotation tensors. This can viewed as a Taylor series expansion. Following the idea of Lumley, numerous non-linear eddy-viscosity models were proposed, see, for example, Refs. [24]--[29].

In the following, we shall present one recent approach proposed by Shih et al. [25]. It includes up to third-order terms in the general eddy-viscosity formulation and is particularly suited to swirling flows. As already pointed out in [19], p. 194, cubic terms are essential for high accuracy. The Reynolds stresses \(\tau^R_{ij}\) can be expressed as [25], [30] (cf. Eq. (7.24))

with \(\overline{S}_{ij}\), \(\overline{\Omega}_{ij}\) according to Eqs. (7.3) and (7.4). Furthermore,

The values of the turbulent kinetic energy \(K\) and the dissipation rate \(\varepsilon\) are obtained from low-Reynolds \(K\)-\(\varepsilon\) turbulence model (cf. Subsection 7.2.2). The factors \(C_1\) to \(C_5\) in Eq. (7.30) are provided in Refs. [25], [30].

In comparison to the linear eddy-viscosity approach, the non-linear models are computationally only slightly more expensive, but they offer a substantially improved prediction capabilities for complex turbulent flows.

It is possible to derive exact equations for the Reynolds stresses by taking the time average (second-order moment)

where \(\mathcal{N}(v_i)\) denotes the Navier-Stokes operator, i.e.,

Using the average Eq. (7.32) together with Eq. (7.33), we obtain the following Reynolds-stress transport equation [31]

for incompressible flow. The formulation for compressible flows can be found in Ref. [17], p. 179, or in [32], [33]. The production of the turbulent kinetic energy \(P_{ij}\), the pressure-strain term \(\Pi_{ij}\), the dissipation-rate term \(\varepsilon_{ij}\), and the third-order diffusion term \(C_{ijk}\) in Eq. (7.34) are defined as

In Eq. (7.35), \(S'_{ij}\) denotes the fluctuating part of the strain-rate tensor. The first part of \(C_{ijk}\), the triple velocity term, represents transport driven by fluctuating convection, the two other parts are the pressure transport terms (pressure-velocity correlations), respectively.

As we can see, the exact Reynolds-stress equation contains new unknown higher-order correlations (e.g., \(\overline{v'_iv'_jv'_k}\)). Therefore, Equation (7.34) can be closed only by using empirical models. This is caused by the non-linear nature of the Navier-Stokes equations. The second-order closures -- the Reynolds-stress models -- provide the necessary framework for solving Eq. (7.34). Examples of implementations can be found, e.g., in Refs. [34]--[36].

The Spalart-Allmaras one-equation turbulence model [38] employs transport equation for an eddy-viscosity variable \(\tilde{\nu}\). It was developed based on empiricism, dimensional analysis and Galilean invariance. It was calibrated using results for 2-D mixing layers, wakes and flat-plate boundary layers. The Spalart-Allmaras model also allows for reasonably accurate predictions of turbulent flows with adverse pressure gradients. Furthermore, it is capable of smooth transition from laminar to turbulent flow at user specified locations. The Spalart-Allmaras model has several favourable numerical features. It is "local" which means that the equation at one point does not depend on the solution at other points. Therefore, it can be readily implemented on structured multi-block or on unstructured grids. It is also robust, converges fast to steady-state and requires only moderate grid resolution in the near-wall region.

The Spalart-Allmaras turbulence model can be written in tensor notation as follows [38]

The terms on the right-hand side represent eddy-viscosity production, conservative diffusion, non-conservative diffusion, near-wall turbulence destruction, transition damping of production, and transition source of turbulence. Furthermore, \(\nu_L = \mu_L/\rho\) denotes the laminar kinematic viscosity and \(d\) is the distance to the closest wall (for the computation of wall distances see Ref. [39]). The turbulent eddy viscosity in Eq. (7.28) and (7.29) is obtained from

The production term is evaluated with the following formulae

In Equation (7.38), \(S\) stands for the magnitude of the mean rotation rate, i.e.,

where \(\Omega_{ij}\) is given by Eq. (7.4). Note that \(\tilde{S}\) differs from its original definition in [38]. The modification was suggested by Spalart in order to prevent \(\tilde{S}\) from reaching zero (cf. Ref. [40], p. 155).

The terms controlling the destruction of the eddy viscosity read

Functions used for modelling the laminar-turbulent transition are given by

where \(\omega_t\) represents the vorticity at the wall at the trip point (position has to be specified by the user), \(\|\Delta\vec{v}\|_2\) denotes the 2-norm of the difference between the velocity at the trip point and the current field point, \(d_t\) is the distance to the nearest trip point, and \(\Delta x_t\) stands for the spacing along the wall at the trip point.

Finally, the various constants in Eqs. (7.36)-(7.40) are defined as

In order to account for non-equilibrium effects on the production term, it was recently proposed to express \(C_{b1}\) as a function of the strain rate [41].

As pointed out in [38], it is convenient to substitute the non-conservative diffusion term in Eq. (7.36), i.e.,

by the following expression

In this way, difficulties with the discretisation of the term \(\left(\partial\tilde{\nu}/\partial x_j\right)^2\) are circumvented.

The Spalart-Allmaras turbulence model (Eq. (7.36)) reads after the transformation into the finite-volume framework as follows

where \(\Omega\) represents the control volume, \(\partial\Omega\) its surface, and \(dS\) is a surface element of \(\Omega\). The convective flux is defined as

with \(V\) being the contravariant velocity (see Eq. (2.22)). The convective flux is in general discretised using a first-order upwind scheme. The viscous flux is given by

where \(n_x\), \(n_y\), and \(n_z\) are the components of the unit normal vector. The normal viscous stresses read

Finally, the source term in Eq. (7.43) becomes

Alternatively, the non-conservative diffusion can be formulated as suggested by Eq. (7.42). The model constants are provided in Eq. (7.41).

The initial value of \(\tilde{\nu}\) is usually taken as \(\tilde{\nu} = 0.1\,\nu_L\). The same value is also specified at inflow boundaries. At outflow boundaries, \(\tilde{\nu}\) is simply extrapolated from the interior of the computational domain. At solid walls, it is appropriate to set \(\tilde{\nu} = 0\) and hence \(\mu_T = 0\).

The K-\(\varepsilon\) turbulence model is the most widely employed two-equation eddy-viscosity model. It is based on the solution of equations for the turbulent kinetic energy \(K\) and the turbulent dissipation rate \(\varepsilon\). The historic roots of the K-\(\varepsilon\) model reach to the work of Chou [42]. During the 1970's, various formulations of the model were proposed. The most important contributions were due to Jones and Launder [43], [44], Launder and Sharma [45] as well as due to Launder and Spalding [46].

The K-\(\varepsilon\) turbulence model requires addition of the so-called damping functions in order to stay valid through the viscous sublayer to the wall. The aim of the damping functions is to assure proper limiting behaviour of \(K\) and \(\varepsilon\) at the wall, i.e.,

where \(y\) represents the coordinate normal to the wall. Further, it can be shown that the Reynolds shear stress behaves like (see, e.g., [17] pp. 138-139)

The K-\(\varepsilon\) models with damping functions are also denoted as low Reynolds number models. The most widely used formulations of the damping functions were proposed by Jones and Launder [43], Launder and Sharma [45], Lam and Bremhorst [47], and by Chien [48]. The reader may find a comparison of seven different low Reynolds number K-\(\varepsilon\) models in Ref. [49].

The K-\(\varepsilon\) turbulence model is more difficult to solve numerically than the previously discussed Spalart-Allmaras model (Subsection 7.2.1). Particularly, the damping functions lead to turbulence equations with stiff source terms. This, and the necessary high grid resolution nearby walls (in order to resolve the viscous sublayer), requires the utilisation of at least point-implicit or better full-implicit time-stepping schemes. Reference [50] contains useful hints on the explicit time discretisation of the K-\(\varepsilon\) equations. Examples of implementations of the K-\(\varepsilon\) model on structured as well as on unstructured grids can be found, e.g., in [51]-[59]. Finally, it is important to note that the accuracy of the K-\(\varepsilon\) model degrades for flows with adverse pressure gradient [49], [17].

A low Reynolds number K-\(\varepsilon\) model can be written as

The terms on the right-hand side represent conservative diffusion, eddy-viscosity production and dissipation, respectively. Furthermore, \(\phi_{\varepsilon}\) denotes the so-called explicit wall term. The Favre-averaged turbulent stresses \(\tau_{ij}^F\) are given by Eq. (7.25) and the strain-rate tensor \(S_{ij}\) follows from Eq. (7.3). The turbulent eddy viscosity in Eq. (7.28) and (7.29) results from

The turbulent kinetic energy is also employed for the evaluation of the eddy viscosity according to Eq. (7.24) or Eq. (7.25). The quantity \(\varepsilon^{*}\) is related to the turbulent dissipation rate \(\varepsilon\) by

The term \(\varepsilon_w\) is the value of the dissipation rate at the wall. The definition in Eq. (7.52) greatly simplifies the application of wall-boundary conditions (see further below).

The constants, the near-wall damping functions as well as the wall term differ between the various K-\(\varepsilon\) models. Here, we choose the Launder-Sharma model because it gives good results for a wide range of applications [49]. For the Launder-Sharma model, the constants and the turbulent Prandtl number are given by [45]

Furthermore, the near-wall damping functions read

with \(Re_T = \rho K^2/(\varepsilon^{*}\mu_L)\) being the turbulent Reynolds number.

Finally, the explicit wall term \(\phi_{\varepsilon}\) and the value \(\varepsilon_w\) are defined as

where \(v_s\) stands for the velocity parallel to the wall, and \(y_n\) represents the coordinate normal to the wall. In order to avoid an explicit knowledge of the wall distance and orientation, it is common to compute the wall term and \(\varepsilon_w\) from the following Cartesian tensor form [60], [57]

instead of Eq. (7.55).

Written in time-dependent integral form for a control volume \(\Omega\) with a surface element \(dS\), the low Reynolds number K-\(\varepsilon\) turbulence model reads

The vector of the conservative variables takes the form

The vector of the convective fluxes is defined

where \(V\) denotes the contravariant velocity (see Eq. (2.22)). The vector of the viscous fluxes is given by

with the normal turbulent viscous stresses

In Eq. (7.60), \(n_x\), \(n_y\), \(n_z\) represent the components of the outward-facing unit normal vector of the surface \(\partial\Omega\). The source term is evaluated from

where \(P\) denotes the production term of the turbulent kinetic energy. It is defined as

with the Favre-averaged turbulent stresses \(\tau_{ij}^F\) given by Eq. (7.25). The constants, the near-wall damping functions as well as the wall term follow for the Launder-Sharma model from the definitions in Eqs. (7.52)-(7.56). The turbulent eddy viscosity \(\mu_T\) is obtained from Eq. (7.51).

The simplest approach is to initialise \(K\) and \(\varepsilon^{*}\) with their freestream values. A better alternative consists of prescribing profiles for \(K\) and \(\varepsilon^{*}\) near solid walls. The profiles can be obtained from analogy to turbulent flat-plate boundary layer [51]. However, this requires the knowledge of wall distances which may not be readily available like it is the case on unstructured grids.

The proper boundary conditions at solid walls are \(K = 0\) and \(\varepsilon^{*} = 0\), provided the transformation in Eq. (7.52) is utilised. This also implies \(\mu_T = 0\) at walls. At inflow boundaries, \(K\) and \(\varepsilon^{*}\) can be computed from relations for the turbulent intensity and length scale, i.e.,

where we assumed \(\varepsilon_{\infty}^{*} = \varepsilon_{\infty}\). In turbomachinery, \(\left(l_T\right)_{\infty}\) is chosen between \(10^{-3}\) and \(10^{-2}\) times the mean radial blade spacing [61]. The values of \(K\) and \(\varepsilon^{*}\) are extrapolated from the interior at outflow boundaries.

As we already noted, the low Reynolds number models require very fine grids at walls. The standard condition is that the first node (or cell centroid) should be located at the distance \(y^{+} \le 1\) from the wall. In order to reduce the the stiffness of the turbulence equations and to save a number of grid points/cells, coarser grids with \(10 \le y^{+} \le 100\) are often employed. In such a case, the K-\(\varepsilon\) model Eq. (7.50) or Eq. (7.57) is applied without the damping functions (\(f_{\mu} = f_{\varepsilon1} = f_{\varepsilon2} = 1\); \(\varepsilon_w = 0\)) and the wall term (\(\phi_{\varepsilon} = 0\)). We speak here of a high Reynolds number turbulence model. Apparently, the distance between the first node (cell centroid) and the wall has to be bridged by the so-called wall functions. The wall functions deliver the values of \(K\) and \(\varepsilon^{*}\) at the node (cell centroid) adjacent to the wall. The turbulence equations are not solved at the wall itself and at the first layer of nodes (cells). Various formulation of the wall functions are used, in general based on the logarithmic wall-law. One example is the function of Spalding [62], which models the viscous sublayer, the transition region as well as the logarithmic layer. Implementations of high Reynolds number models were described, e.g., in [63]-[67] or [58].

The application of the wall functions leads (provided the grid is not too coarse) to reasonably accurate results for attached boundary layers. It also allows the utilisation of purely explicit time-stepping schemes. However, the use of wall functions becomes highly questionable for separated flows.

The K-\(\omega\) Shear Stress Transport (SST) turbulence model of Menter [68], [69] merges the K-\(\omega\) model of Wilcox [70], [17] with a high Reynolds number K-\(\varepsilon\) model (transformed into the K-\(\omega\) formulation). The SST model seeks to combine the positive features of both models. Therefore, the K-\(\omega\) approach is employed in the sublayer of the boundary layer. The reason is that the K-\(\omega\) model needs no damping function. This leads, for similar accuracy, to significantly higher numerical stability in comparison to the K-\(\varepsilon\) model. Furthermore, the K-\(\omega\) model is also utilised in logarithmic part of the boundary layer, where it is superior to the K-\(\varepsilon\) approach in adverse pressure flows and in compressible flows. On the other hand, the K-\(\varepsilon\) model is employed in the wake region of the boundary layer because the K-\(\omega\) model is strongly sensitive to the freestream value of \(\omega\) [71]. The K-\(\varepsilon\) approach is also used in free shear layers since it represents a fair compromise in accuracy for wakes, jets, and mixing layers.

One distinct feature of the SST turbulence model is the modified turbulent eddy-viscosity function. The purpose is to improve the accuracy of prediction of flows with strong adverse pressure gradients and of pressure-induced boundary layer separation. The modification accounts for the transport of the turbulent shear stress. It is based on the observation of Bradshaw that the principal shear stress is proportional to the turbulent kinetic energy.

A certain disadvantage of the SST model is that distances to the nearest wall have to be known explicitly. This requires special provisions on multiblock structured or on unstructured grids. See, e.g., Ref. [39] for the computation of wall distances. Examples for applications of the SST turbulence model can be found in [72]-[74].

The transport equations for the turbulent kinetic energy and the specific dissipation of turbulence read in differential form [68]

The terms on the right-hand side of Eq. (7.65) represent conservative diffusion, eddy-viscosity production and dissipation, respectively. Furthermore, the last term in the \(\omega\)-equation describes the cross diffusion. The Favre-averaged turbulent stresses \(\tau_{ij}^F\) are given by Eq. (7.25) and the strain-rate tensor \(S_{ij}\) follows from Eq. (7.3). The turbulent eddy viscosity in Eq. (7.28) and (7.29) is obtained from [68]

This definition of the turbulent viscosity guarantees that in an adverse pressure gradient boundary layer, where the production of \(K\) is larger than its dissipation \(\omega\) (hence \(a_1\omega < \|\text{curl}\,\vec{v}\|_2\)), Bradshaw's assumption, i.e., \(\tau = a_1\rho K\) (shear stress proportional to turbulent kinetic energy) is satisfied.

The function \(f_1\) in Eq. (7.65), which blends the model coefficients of the K-\(\omega\) model in boundary layers with the transformed K-\(\varepsilon\) model in free-shear layers and freestream zones, is defined as

where \(d\) stands for the distance to the nearest wall and \(CD_{K\omega}\) is the positive part of the cross-diffusion term in Eq. (7.65), i.e.,

The auxiliary function \(f_2\) in Eq. (7.66) is given by

The model constants are as follows

Finally, the coefficients of the SST turbulence model \(\beta\), \(C_{\omega}\), \(\sigma_K\), and \(\sigma_{\omega}\) are obtained by blending the coefficients of the K-\(\omega\) model, denoted as \(\phi_1\), with those of the transformed K-\(\varepsilon\) model (\(\phi_2\)). The corresponding relation reads

The coefficients of the inner model (K-\(\omega\)) are given by

The coefficients of the outer model (K-\(\varepsilon\)) are defined as

The integral formulation of the SST turbulence model corresponds, in principle, to that of the K-\(\varepsilon\) model from Subsection 7.2.2. Therefore, it is not repeated here.

The boundary conditions for the kinetic turbulent energy and the specific dissipation at solid walls are

with \(d_1\) being the distance of the first node (cell centroid) from the wall. The grid has to be refined such that \(y^{+} < 3\).

For the inflow boundaries, the following freestream values are recommended

where \(L\) denotes the length of the computational domain, \(1 \le C_1 \le 10\) and \(2 \le C_2 \le 5\), respectively. The values of \(K\) and \(\omega\) are extrapolated from the interior at outflow boundaries.

LES is based on a spatial filtering operation, which decomposes any flow variable \(U\) into a filtered (large-scale, resolved) part \(\overline{U}\) and into a sub-filter (unresolved) part \(U'\), i.e.,

The filtered variable at the location \(\vec{r}_0\) in space is defined as

where \(\Omega\) denotes the entire flow domain, \(G\) represents the filter function, and \(\vec{r}\) is the position vector, respectively. The filter function determines the structure and size of the small scales. The filter function depends on the difference \(\vec{r}_0 - \vec{r}\) and on the filter width \(\Delta = \left(\Delta_1\,\Delta_2\,\Delta_3\right)^{1/3}\), with \(\Delta_i\) being the filter width in the \(i\)-th spatial coordinate. The following filter functions are the mostly used ones (see Fig. 7.3):

• the tophat filter:
• The sharp Fourier cut-off filter:
• The Gaussian filter:

The tophat and the Gaussian filter smooth the large-scale fluctuations as well as the small scales below the filter width. The cut-off filter affects only the scales below the cut-off wave-number. In practice, the Gaussian filter is always employed in conjunction with a sharp Fourier cut-off. Filters suitable for grids with varying cell sizes were proposed in Refs. [114], [115].

The spatial filtering, defined by Eq. (7.76) and Eq. (7.77), has to be applied to the Navier-Stokes equations in order to remove the small turbulent scales. The filter width \(\Delta\) as well as the filter function are considered as free parameters. In fact, the governing equations are usually not explicitly filtered. Instead, the grid as well as the discretisation errors are assumed to define the filter \(G\). For the discussion of explicit filtering see Refs. [116], [117].

Because of the differing treatment, we shall distinguish in the following between compressible (7.1) and incompressible (7.6) formulation of the Navier-Stokes equations.

For an incompressible flow of a Newtonian fluid, the filtered governing equations (7.6) take the form

where \(\nu\) denotes the kinematic viscosity coefficient. The equations (7.81) describe the temporal and spatial evolution of the large, energy-carrying scales of motion. The non-linearity of the convective term leads to the appearance of the so-called subgrid-scale stress (SGS) tensor

which describes the effects of the unresolved scales. The SGS tensor has to be modelled (see Subsection 7.3.3) in order to close the equations.

The SGS tensor can be decomposed into three parts [118], namely

The individual parts have the following physical meaning:

is the so-called Leonard stress term and represents the interactions between large-scale eddies which produce small-scale turbulence. This term only can be evaluated explicitly from the filtered velocity field \(v_i\). Further, the cross-stress term

describes interactions between large- and small-scale eddies. Finally,

is the so-called SGS Reynolds-stress tensor. It reflects interactions between the small-scale structures. The above decomposition (7.83) is no longer used mainly because \(L_{ij}\) and \(C_{ij}\) are not invariant with respect to Galilean transformation¹.

If LES is to be applied to compressible flows, we have to apply Favre averaging (Subsection 7.1.2) together with the spatial filtering to the Equations (7.1). Otherwise, the filtered Navier-Stokes equations would contain products between density and other variables like velocity or temperature. Thus, the velocity components, the energy and the temperature in Eq. (7.1) is decomposed as

The filtered variable at the location \(\vec{r}_0\) in space is given by

where the overbar denotes the filtering in Eq. (7.77). The Favre-filtered Navier-Stokes equations (7.1) read [111], [113]

with the terms

and

In the above equations (7.89)-(7.90), \(e\) denotes internal energy per unit mass, \(\tilde{S}_{ij}\) is the Favre-filtered strain-rate tensor, and \(\tau_{ij}^{SF} = \overline{\rho}\left(\widetilde{v_i v_j} - \tilde{v}_i\tilde{v}_j\right)\) represents the Favre-averaged subgrid-scale stress. Furthermore, \(\mu\), \(\mu_B\), and \(k\) stand for the molecular viscosity, the bulk viscosity, and for the thermal conductivity, respectively. Finally, \(\tilde{\mu}\), \(\tilde{\mu}_B\), and \(\tilde{k}\) are the corresponding values at the filtered temperature \(\tilde{T}\).

The right-hand side of Eq. (7.89) contains terms which have to be modelled. In the momentum equation, the SGS stresses \(\tau_{ij}^{SF}\) are approximated, but the second term, i.e., \(\left(\overline{\sigma}_{ij} - \hat{\sigma}_{ij}\right)\) is usually neglected. In the energy equation, term \(\mathcal{A}\) can be expressed through the SGS stresses [119], term \(\mathcal{B}\) can be neglected, and terms \(\mathcal{C}\), \(\mathcal{D}\) can be modelled as proposed in [120].

The main task of a subgrid-scale model is to simulate energy transfer between the large and the subgrid scales. On the average, the energy is transported from the large scales to the small ones (turbulent cascade process). Therefore, a subgrid-scale model has to provide means of adequate energy dissipation. However, in some instances the energy also flows from small to large scales - a process called backscatter. Thus, the model should account for this effect as well. Backscatter models are discussed, e.g., in [121].

Various subgrid-scale models were proposed in the past and the research still continues. The models can be divided into two basic classes. The first one consists of approaches which model the SGS tensor \(\tau_{ij}^S\) explicitly. A necessary condition is then that the numerical dissipation caused by the spatial discretisation scheme must be much lower than the subgrid-scale dissipation. The majority of explicit SGS models is based on the eddy-viscosity concept, which is explained next. Furthermore, we shall present the Smagorinsky model, which forms the basis of all subgrid-scale models. We shall also briefly discuss the basics of the so-called dynamic subgrid-scale models. A comparison of six different explicit subgrid-scale models was presented recently in [122].

The second class of subgrid-scale models consists of approaches where the SGS stresses are modeled implicitly by an appropriate discretisation of the convective fluxes (thus \(\tau_{ij}^S\) is omitted). These models are referred to as Monotonically Integrated LES (MILES). The methodology was first presented by Boris et al. [123] and recently advocated by Grinstein and Fureby [124], [125]. The condition in this case is that the numerical dissipation correctly models the subgrid-scale dissipation. This is not quite easy to achieve, as the investigations in Ref. [126] and [127] prove. Nevertheless, MILES was applied with some success to a variety of flow problems [128]-[134].

These explicit models are able to represent the global dissipative effects of the small scales, but they cannot reproduce the local details of the energy exchange. In the case of incompressible flows, the eddy-viscosity models relate the SGS stresses to the large-scale strain-rate \(\overline{S}_{ij}\) as follows

The strain-rate \(\overline{S}_{ij}\) is obtained from Eq. (7.3) by using filtered velocity components. The eddy viscosity \(\nu_T\) is in general evaluated from algebraic relations in order to save numerical costs. The isotropic part of the SGS stresses (\(\tau_{kk}^S\)) can either be added to the filtered pressure [135], modelled [136] or neglected.

Relation similar to Eq. (7.91) applies also in the case of compressible Navier-Stokes equations. The components of the Favre-averaged SGS stress tensor are approximated as

The components of the strain-rate tensor \(\tilde{S}_{ij}\) are given in Eq. (7.90).

The Smagorinsky model [75] is based on the equilibrium hypothesis which implies that the small scales dissipate entirely and instantaneously all the energy they receive from the large scales. The algebraic model assumes the form

where \(|\overline{S}| = \left(2\overline{S}_{ij}\overline{S}_{ij}\right)^{1/2}\) is the magnitude of the strain-rate tensor and \(C_s\) denotes the Smagorinsky constant. The theoretical value found by Lilly [137] is \(C_s \approx 0.18\). However, the Smagorinsky constant depends on the type of the flow. For example, in shear flows \(C_s\) has to be reduced to approximately 0.1. The filter width \(\Delta\) in Eq. (7.93) is usually chosen to be twice the average grid size, i.e., \(\Delta = 2\left(\Delta x_1\,\Delta x_2\,\Delta x_3\right)^{1/3}\).

In order to account for the reduced growth of the small scales near walls, the value of the eddy viscosity \(\nu_T\) has to be reduced. Thus, the Smagorinsky model Eq. (7.93) is modified according to Van Driest damping as

where \(y^{+}\) represents the dimensionless wall distance (for the computation of wall distances see, e.g., Ref. [39]).

The Smagorinsky model is numerically cheap and easy to implement. However, it has several serious disadvantages:

• it is too dissipative in laminar regions with mean shear;
• it requires special provisions near walls and at laminar-turbulent transition;
• the parameter \(C_s\) is not uniquely defined;
• the process of energy backscatter is not modelled.

Because of these shortcomings, various other approaches were proposed (see, e.g., [111]). Very popular are the following dynamic models.

The dynamic SGS models employ the same relation as the Smagorinsky model (Eq. (7.93)) for the evaluation of the eddy viscosity \(\nu_T\) in Eq. (7.91) or (7.92). The difference is that the Smagorinsky constant (adjusted a priori) is replaced by a parameter, which evolves dynamically in space and in time. Hence,

The parameter \(C_d\) is computed based on the energy content of the smallest scale of the turbulence. For this purpose, Germano et al. [138] proposed to employ a second filter - the so-called test filter \(\hat{\Delta}\). The width of the test filter has to be larger than that of the filter \(\Delta\) applied to the governing equations (usually \(\hat{\Delta} = 2\Delta\)). The application of the test filter to the filtered equations leads to the so-called subtest-scale stresses \(\tau_{ij}^{ST}\)

The subtest-scale stresses are related to the SGS stresses \(\tau_{ij}^S\) (Eq. (7.83)) via the Germano identity [139]

where \(\hat{L}_{ij}\) denotes the Leonard stresses associated with the test filter. It represents the contribution to the Reynolds stresses by the scales whose length is intermediate between the filter width \(\Delta\) and the test filter width \(\hat{\Delta}\).

If we express the subtest-scale and SGS stresses in Eq. (7.97) using the eddy-viscosity approach Eq. (7.91) together with Eq. (7.95), we obtain

with

The notation \([\,]^{\wedge}\) means that the whole term enclosed in the square brackets is test-filtered. The parameter \(C_d\) can be derived from Eq. (7.98) by using the least-squares minimisation of Lilly [140]. This leads to

The above formulation (7.100) is mathematically inconsistent since the parameter \(C_d\) was taken outside the test filter in Eq. (7.98). In practice, the numerator and denominator in Eq. (7.100) are therefore ensemble-averaged in the homogeneous directions, i.e.,

Improved dynamic SGS models were proposed, e.g., by Ghosal et al. [141], Carati et al. [142], Piomelli and Liu [143], and Held [90].

The computational costs of LES of wall-bounded flows at high Reynolds numbers (\(Re > 10^6\)) are still too high for engineering purposes. The reason is the excessively large number of grid points (cells) required to resolve the wall layer appropriately. In order to reduce the costs, it is possible to model the wall layer by specifying a correlation between the velocity in the outer flow and the stress at the wall. This approach is quite similar to using wall functions in RANS simulations. The basic assumption is that there is only a weak interaction between the near-wall and the outer region, which is supported by the investigations in [144] and [145].

Earlier implementations of the wall models were based on the assumption that the dynamics of the wall layer are universal and hence they can be approximated by a generalised law-of-the-wall. Basically, the models utilised the logarithmic law (see [77] and [146]-[148]). Balaras et al. [149] proposed recently a new zonal approach. Within the two-layer model, the filtered Navier-Stokes equations (7.81) are solved up to the first grid point above the wall. From this point to the wall 2-D boundary layer equation are solved on a refined embedded grid. The solution on the embedded grid is then used to prescribe the wall shear stress as a boundary condition for the LES. The zonal approach of Balaras et al. [149] allows it to place the first point in a region \(20 < y^{+} < 100\), which leads to significantly reduced grid size and hence computational time. The methodology was applied with success to turbulent flows in a plane channel, square duct and rotating channel. Later on, it was also employed for the LES of separated flows with encouraging results [150]-[153].

Even though the above wall models help to reduce the number of grid points (cells) considerably, LES still remains too costly for complex engineering configurations. For this reason, Spalart recently suggested another approach, the so-called Detached Eddy Simulation (DES), which is aimed at the simulation of high Reynolds-number massively separated flows [154], [155]. The methodology represents a hybrid between the RANS and LES. The idea is to employ highly stretched grids together with a RANS turbulence model (mostly the Spalart-Allmaras model from Subsection 7.2.1 or Menter's SST model from Subsection 7.2.3) to resolve the attached boundary layer(s), and to use LES outside the wall region together with an isotropic grid to capture the detached 3-D eddies. Thus, DES tries to combine the strengths of both methods in a single framework.

The implementation of DES for the Spalart-Allmaras turbulence model consists of replacing the wall distance \(d\) in the formulae (7.36) and (7.38)-(7.40) by the DES length scale

The length scale is dependent on the largest dimension of the control volume, i.e., \(\Delta = \max(\Delta x, \Delta y, \Delta z)\). The constant \(C_{DES}\) depends to some extent on the type of the flow. For a homogeneous turbulence, the value \(C_{DES} = 0.65\) was found optimal [156]. On the other hand, \(C_{DES} = 0.1\) was recommended for transonic and supersonic jets [157]. The definition of the length scale \(l\) in Eq. (7.102) makes sure that within the boundary layer, where \(d < C_{DES}\Delta\) and hence \(l = d\), the original RANS model is recovered. On the other hand, outside the boundary layer \(l = C_{DES}\Delta\) and the Spalart-Allmaras model serves as a one-equation SGS model for the LES (cf. Eq. (7.91) and (7.24)). When integrating the governing equations in time, the global time step has to be adjusted such as to resolve the time scales of the detached eddies. This usually means that the time step would by far exceed the stability margin of an explicit scheme for the boundary layer region. It is therefore more efficient to employ a time-accurate implicit scheme, such as the dual time-stepping approach described in Section 6.3. More details regarding the DES methodology and examples of simulations can be found in the above references or in [158]-[161].