7.2 First-Order Closures 一阶封闭[cfd-7-2]

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.