Turbulence Models -- RANS, LES, DES, y+ discipline

The turbulence problem

Turbulence is the dominant challenge in predictive CFD. The Navier-Stokes equations are deterministic. Direct numerical simulation (DNS) resolves all turbulent scales. DNS predicts everything in principle. The scales range from the integral length scale L down to the Kolmogorov length scale eta. L contains the largest energy-carrying eddies. eta is where viscous dissipation dominates. The grid requirement is:

$$N \in (L/eta)^3 \in Re_L^{9/4}$$

For a low-speed flow over a car with Reynolds number 10^6 based on vehicle length, N ~ 10^{13.5} cells. This number is far beyond any practical solver. Even with HPC, DNS is confined to Re less than 10^4 for canonical flows (channel, cylinder, mixing layer). For all engineering Reynolds numbers, turbulence modelling is mandatory.

Reynolds-Averaged Navier-Stokes (RANS)

RANS is the workhorse of industrial CFD. The concept decomposes each instantaneous variable into a mean and fluctuating component:

$$phi(x, t) = overbar{phi}(x) + phi prime (x, t)$$

Substituting into the Navier-Stokes equations and averaging produces the RANS equations:

$$partial overbar{rho} / partial t + partial / partial x_i (overbar{rho} overbar{U}_i) = 0$$

$$partial overbar{rho} overbar{U}_i / partial t + partial / partial x_j (overbar{rho} overbar{U}_i overbar{U}_j) = - partial overbar{p} / partial x_i + partial / partial x_j \left[ mu \left( partial overbar{U}_i / partial x_j + partial overbar{U}_j / partial x_i - 2/3 overbar{rho} overbar{u prime } prime _i overbar{u prime } prime _j \right) right\right]$$

The new term tau_{t,ij} equals minus overbar{rho} overbar{u prime_i u prime_j} is the Boussinesq Reynolds stress tensor. It is unknown and must be modeled. This is the closure problem of turbulence. The Boussinesq hypothesis from 1877 closes it by introducing turbulent viscosity nu_t:

$$- overbar{rho} overbar{u prime_i u prime_j} = mu_t \left( partial overbar{U}_i / partial x_j + partial overbar{U}_j / partial x_i \right) - 2/3 delta_{ij} rho k$$

k equals 1/2 overbar{u prime_i u prime_i} is the turbulent kinetic energy.

The closure problem reduces from modeling a tensor to modeling a scalar nu_t. But nu_t is unknown. Additional transport equations are needed. Different RANS models solve the closure problem with different numbers of transport equations.

Zero-equation models

Spalart-Allmaras uses a single transport equation for a working variable tilde{nu}. It is related to nu_t. It is not a zero-equation model in the strictest sense. It solves one equation. Still, it is grouped with zero-equation models. No length or time scale transport is required.

Features: extremely robust, widely used in aerospace for external aerodynamics. It is built into Fluent, Star-CCM+, and OpenFOAM as spalartAllmarasFoam. It has limited performance in separation. It performs poorly in strong pressure gradients. This model is default for external aero in many aerospace codes.

One-equation models

Spalart-Allmaras is the only widely-used one-equation model. The Wilcox k-omega 1988 formulation (before SST blending) is technically two-equation.

Two-equation models

k-Epsilon: The standard k-epsilon model is the most widely used model in industrial CFD by sheer numbers. It solves transport equations for turbulent kinetic energy k and its dissipation rate epsilon. The production term is P_k equals mu_t S^2 where S equals square root of 2 S_{ij} S_{ij}.

The standard model has well-known weaknesses:

• It assumes isotropic turbulence (inaccurate for strong rotation or curvature)
• It cannot predict separation from adverse pressure gradients well
• It over-predicts spreading of free shear layers (jets, wakes)
• Wall-bounded flow requires wall functions (cannot integrate to wall)

k-Omega SST: Spalart Shear Stress Transport model blends robustness and near-wall accuracy of k-omega with free-stream independence of k-epsilon. Near the wall (y+ resolved), it uses k-omega formulation. In the free stream (y+ much greater than 300), it blends to k-epsilon.

F_1 equals tanh(arg_1^4) is the blending function. F_1 equals 0 means k-epsilon. F_1 equals 1 means k-omega. d is distance to nearest wall.

The SST model has specific corrections:

• Production limiter: prevents P_k from becoming too large (excessive turbulence generation in adverse pressure gradients)
• Shear stress limitation: limits tau_t to incompressible value 2 rho k near walls

Features of SST:

• Default choice for wall-bounded flows (aerodynamics, turbomachinery, internal flow with separation)
• Performs well for mild-to-moderate separation
• Well-implemented in all platforms (Fluent, Star-CCM+, OpenFOAM)
• Requires y+ about 1 for wall-resolved LES-RANS

Large Eddy Simulation (LES)

LES solves the large energy-carrying eddies directly. It models only the small dissipative subgrid scales (SGS). The fundamental advantage over RANS is that LES resolves anisotropic unsteady structures. The disadvantage is computational cost.

LES filtering is applied to the Navier-Stokes equations:

$$partial overbar{U}_i / partial x_i = 0$$

$$partial overbar{U}_i / partial t + partial / partial x_j (overbar{U}_i overbar{U}_j) = - 1/rho partial overbar{p} / partial x_i + nu partial^2 overbar{U}_i / partial x_j partial x_j - partial tau_{ij}^{SGS} / partial x_j$$

tau_{ij}^{SGS} equals overbar{U_i U_j} minus overbar{U}_i overbar{U}_j is the subgrid-scale stress tensor. This is the missing stress from filtering.

Smagorinsky model (1963)

The simplest SGS model uses the eddy-viscosity concept:

$$tau_{ij}^{SGS} - 1/3 tau_{kk}^{SGS} delta_{ij} = -2 nu_t overbar{S}_{ij}$$

$$nu_t = (C_s Delta)^2 |overbar{S}|$$

Delta is the grid filter width. Delta is typically (Delta_x Delta_y Delta_z)^{1/3}. C_s is in range 0.1 to 0.2 and is the Smagorinsky constant.

Problem: C_s is not universal. The model dissipates too much energy near walls. It does not vanish as y+ goes to 0. Wall-modeled LES requires a damping function: f_nu equals 1 minus exp(-(y+/A+)^3) with A+ about 26.

WALE model (2000, Nicoud and Ducros)

WALE improves on Smagorinsky. It gives the correct near-wall scaling nu_t ~ (y+)^3. It needs no damping functions:

$$nu_t = (C_w Delta)^2 overline{(S^d_{ij} S^d_{ij})^{3/2} / ((S_{ij} S_{ij})^{5/2} + (S^d_{ij} S^d_{ij})^{5/4})}$$

S^d_{ij} is the symmetric part of the square of the velocity gradient tensor. WALE is now the recommended SGS model for most LES applications. It has better near-wall behavior than Smagorinsky. It has fewer tuning parameters than dynamic models.

Dynamic models (Germano et al., 1991; Lilly, 1992)

The dynamic Smagorinsky model computes C_s locally from the resolved flow field. It uses a two-level filter:

1. Apply a test filter at scale hat{Delta} equals 2 Delta 2. Compare resolved stresses at Delta and hat{Delta} 3. Apply the Germano identity to solve for C_s 4. Average C_s in homogeneous directions (or minimize error in the least-squares sense)

Dynamic models are the most accurate in principle. They adapt to the local flow. They suffer from practical problems:

• Need for directional averaging (not practical for fully 3D unstructured meshes)
• Negative eddy-viscosity regions (cause instability without clipping)
• Higher cost per time step

Practical LES guide

LES is fundamentally transient. It demands much finer grids than RANS. It needs much smaller time steps (CFL less than 1). It is only practical for moderate-Re geometries. The computational cost is 20 to 100 times RANS in wall-clock time. Justification is required.

See the OpenFOAM Turbulence note for exact turbulenceProperties and momentumTransport configuration. See the FVM Overview note for scheme considerations (unique convection schemes are needed for LES).

Detached Eddy Simulation (DES / DDES)

DES is a hybrid RANS-LES approach. RANS handles the attached boundary layer. LES handles separated regions. This combines RANS cost-efficiency for boundary layers with LES fidelity for unsteady separation.

The modification is to the turbulence model length scale:

$$d_{DES} = d_{RANS} - F_{DES} max(C_{DES} Delta, psi d_{RANS})$$

d_{RANS} is the RANS distance to nearest wall. Delta is the local grid spacing. C_{DES} is about 0.618. F_{DES} is a switching function.

Delayed DES (DDES, Menter et al. 2006) adds a shielding term. The term prevents premature switching to LES in regions where RANS should persist. For example, near-wall regions on grids too coarse for LES. This prevents Grid-Induced Separation (GIS). GIS is a well-known failure mode of DES where the coarse RANS grid causes premature transition to LES.

OpenFOAM supports kOmegaSST with DES and DDES options in the momentumTransport dictionary. Fluent and Star-CCM+ implement delayed DES. Fluent uses the SST-SAS variant. Star-CCM+ uses DES and IDDES.

y+ discipline -- wall treatment and mesh requirements

The dimensionless wall distance y+ is the single most important parameter in wall-bounded simulation practice:

$$y+ = y sqrt(tau_w / rho) / nu = y u_tau / nu$$

y is first cell center distance from wall. tau_w is wall shear stress. u_tau equals sqrt(tau_w / rho) is friction velocity. The value of y+ determines which wall treatment is valid:

Wall-resolved LES and DNS require y+ about 1 everywhere. Wall-function RANS requires y+ in range 30 to 300 at every wall cell center.

The danger zone is particularly insidious. The mesh places first cells in the buffer layer. Neither the law-of-the-wall (log-layer) model nor the viscous-sublayer-resolving model is valid. The result is systematically incorrect drag and heat-transfer predictions.

nutUSpaldingWallFunction (OpenFOAM v14 and newer) is a continuous wall function. It bridges the viscous sublayer to the log layer. It eliminates the explicit y+ requirement. It is the safest default for snappyHexMesh cases where first-cell height cannot be guaranteed everywhere.

Inlet turbulence estimation

Two-equation RANS models require inlet boundary values for k and epsilon (or omega). Specifying them incorrectly can bias the entire solution. This is particularly true for separation prediction:

$$k = 3/2 (I U_{ref})^2, \quad epsilon = C_mu^{3/4} k^{3/2} / L, \quad omega = k^{1/2} / (C_mu^{1/4} L)$$

I is turbulence intensity. The value is square root of (2k/3) over U_{ref}. L is turbulence length scale. The value is about 0.07 times characteristic dimension for internal flow. For external flow, it is about 0.07 L_{body}.

Typical external freestream intensity values:

These are starting estimates. Experimental data from PIV, hot-wire, or LDV always beats a calculation.

Cross-references

• FVM Overview -- turbulence transport equations in the general transport framework
• Mesh Quality -- y+ and boundary-layer cell height
• TVD Limiters -- bounded schemes for turbulence equations
• OpenFOAM Turbulence -- kEpsilon, kOmegaSST, LES settings for v14
• OpenFOAM Case Setup -- turbulenceProperties to momentumTransport
• Fluent vs OpenFOAM -- turbulence model comparison across platforms
• Pope, S. B. (2000). Turbulent Flows. Cambridge University Press.
• Menter, F. R. (1994). Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications. AIAA Journal, 32(8), 1598-1605.