Turbulence Modeling

Turbulence is the dominant flow regime in engineering applications. When you drive a car at highway speeds, the Reynolds number based on vehicle length exceeds $10^6$ — firmly in the turbulent regime. Yet turbulence remains one of physics' greatest unsolved problems. CFD must model turbulence because directly resolving it is computationally prohibitive.

What Is Turbulence?

Characteristics

Turbulent flow exhibits:

• Irregularity — Chaotic, seemingly random fluctuations
• Three-dimensionality — Even in nominally 2D geometries
• Diffusivity — Enhanced mixing of momentum, heat, and species
• Dissipation — Kinetic energy cascades to small scales and dissipates
• Wide range of scales — From domain size to Kolmogorov microscale

The Scale Separation Problem

The Kolmogorov length scale (smallest turbulent eddies):

$$\eta = \left(\frac{\nu^3}{\epsilon}\right)^{1/4}$$

The ratio of largest to smallest scales:

$$\frac{L}{\eta} \sim Re^{3/4}$$

For $Re = 10^6$ (typical car):

• $L/\eta \sim 30,000$
• 3D mesh requirement: $(30,000)^3 \sim 10^{13}$ cells
• Time steps: Proportional to smallest scales

Reynolds Averaging

Decomposition

Every turbulent quantity can be split:

$$u = \bar{u} + u'$$

Where:

• $\bar{u}$ = time-averaged (mean) velocity
• $u'$ = fluctuating component

Reynolds-Averaged Navier-Stokes (RANS)

Applying Reynolds averaging to the Navier-Stokes equations:

$$\frac{\partial \bar{u}_i}{\partial t} + \bar{u}_j \frac{\partial \bar{u}_i}{\partial x_j} = -\frac{1}{\rho}\frac{\partial \bar{p}}{\partial x_i} + \nu \frac{\partial^2 \bar{u}_i}{\partial x_j^2} - \frac{\partial \overline{u'_i u'_j}}{\partial x_j}$$

The last term is the Reynolds stress tensor:

$$\tau_{ij}^{turb} = -\rho \overline{u'_i u'_j}$$

The Closure Problem

We have 4 equations (continuity + 3 momentum) but introduced 6 new unknowns (symmetric Reynolds stress tensor). This is the turbulence closure problem — we need additional equations or assumptions.

The Boussinesq Hypothesis

Most RANS models assume the Reynolds stresses are proportional to mean strain:

$$-\overline{u'_i u'_j} = \nu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) - \frac{2}{3}k\delta_{ij}$$

Where:

• $\nu_t$ = turbulent (eddy) viscosity
• $k = \frac{1}{2}\overline{u'_i u'_i}$ = turbulent kinetic energy

Two-Equation Models

Why Two Equations?

Dimensional analysis suggests:

$$\nu_t \propto \text{velocity scale} \times \text{length scale}$$

Two equations provide two independent scales to determine $\nu_t$.

The k-epsilon Model

• $k$ = turbulent kinetic energy
• $\epsilon$ = dissipation rate of $k$

$$\nu_t = C_\mu \frac{k^2}{\epsilon}$$

$$\frac{\partial k}{\partial t} + \bar{u}_j \frac{\partial k}{\partial x_j} = \frac{\partial}{\partial x_j}\left[\left(\nu + \frac{\nu_t}{\sigma_k}\right)\frac{\partial k}{\partial x_j}\right] + P_k - \epsilon$$

$$\frac{\partial \epsilon}{\partial t} + \bar{u}_j \frac{\partial \epsilon}{\partial x_j} = \frac{\partial}{\partial x_j}\left[\left(\nu + \frac{\nu_t}{\sigma_\epsilon}\right)\frac{\partial \epsilon}{\partial x_j}\right] + C_{\epsilon 1}\frac{\epsilon}{k}P_k - C_{\epsilon 2}\frac{\epsilon^2}{k}$$

$$P_k = \nu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) \frac{\partial \bar{u}_i}{\partial x_j}$$

k-epsilon Variants

Limitations of k-epsilon

• Poor near-wall performance — Requires wall functions
• Overpredicts spreading rates in round jets
• Fails for strongly separated flows
• Insensitive to adverse pressure gradients

The k-omega Family

k-omega Model (Wilcox)

Uses specific dissipation rate $\omega = \epsilon / (C_\mu k)$:

$$\nu_t = \frac{k}{\omega}$$

• Better near-wall behavior
• No wall functions needed (integrates to wall)
• Good for boundary layers with pressure gradients

• Sensitive to freestream $\omega$ values

k-omega SST (Shear Stress Transport)

Menter's SST model blends k-epsilon in the freestream with k-omega near walls:

$$F_1 = \tanh\left[\left(\min\left[\max\left(\frac{\sqrt{k}}{\beta^* \omega y}, \frac{500\nu}{y^2 \omega}\right), \frac{4\rho \sigma_{\omega 2} k}{CD_{k\omega}y^2}\right]\right)^4\right]$$

$$\phi = F_1 \phi_1 + (1 - F_1)\phi_2$$

Where $\phi_1$ = k-omega constants, $\phi_2$ = transformed k-epsilon constants.

$$\nu_t = \frac{a_1 k}{\max(a_1 \omega, S F_2)}$$

This prevents overprediction of shear stress in adverse pressure gradient flows.

Model Comparison

When to Use What

Model Hierarchy

One-Equation Models

Spalart-Allmaras

Solves one transport equation for modified turbulent viscosity $\tilde{\nu}$:

$$\frac{D\tilde{\nu}}{Dt} = c_{b1}\tilde{S}\tilde{\nu} - c_{w1}f_w\left(\frac{\tilde{\nu}}{d}\right)^2 + \frac{1}{\sigma}\left[\nabla \cdot ((\nu + \tilde{\nu})\nabla\tilde{\nu}) + c_{b2}(\nabla\tilde{\nu})^2\right]$$

• Simple (one equation)
• Good for attached boundary layers
• Widely validated for aerospace

• Not general-purpose
• Poor for free shear flows

Reynolds Stress Models (RSM)

Instead of assuming Boussinesq, solve transport equations for all six Reynolds stresses:

$$\frac{\partial \overline{u'_i u'_j}}{\partial t} + ... = P_{ij} + D_{ij} + \Pi_{ij} - \epsilon_{ij}$$

• Captures anisotropy
• Better for swirling, rotating flows
• More physics

• 7 additional equations
• More expensive
• Harder to converge
• Requires careful setup

Large Eddy Simulation (LES)

Philosophy

Resolve large, energy-containing eddies directly; model only small, universal scales.

Filtering

Apply spatial filter to Navier-Stokes:

$$\bar{u}_i(\mathbf{x}, t) = \int G(\mathbf{x} - \mathbf{x}', \Delta) u_i(\mathbf{x}', t) d\mathbf{x}'$$

Subgrid-Scale Models

The unresolved scales appear as subgrid-scale stress:

$$\tau_{ij}^{sgs} = \overline{u_i u_j} - \bar{u}_i \bar{u}_j$$

$$\nu_{sgs} = (C_s \Delta)^2 |\bar{S}|$$

LES Requirements

• Turbulent spectra
• Instantaneous flow features
• Better accuracy for separated flows

Wall Treatment

The Near-Wall Challenge

Turbulent boundary layers have distinct regions:

• Viscous sublayer (y+ < 5): Viscosity dominates
• Buffer layer (5 < y+ < 30): Transition
• Log layer (y+ > 30): Universal log law

$$y^+ = \frac{y u_\tau}{\nu}, \quad u_\tau = \sqrt{\tau_w / \rho}$$

Wall Functions

$$u^+ = \frac{1}{\kappa}\ln(y^+) + B$$

• First cell centroid at y+ = 30-300
• Faster computation
• Less accurate for separated flows

Low-Reynolds-Number Models

• First cell at y+ < 1
• Many cells in boundary layer
• More accurate but expensive

y+ Guidelines

Practical Tips

Model Selection Flowchart

• Is flow attached?

- No → Use k-omega SST

• Strong rotation/curvature?

• Need unsteady details?

• Validation data available?

Common Mistakes

Turbulence Intensity

At inlet, specify turbulence based on upstream conditions:

$$k = \frac{3}{2}(U \cdot I)^2$$

Key Takeaways

• Turbulence is chaotic with wide scale ranges — DNS impractical for engineering
• Reynolds averaging introduces Reynolds stresses — the closure problem
• Boussinesq hypothesis models Reynolds stresses with eddy viscosity
• k-epsilon is robust for free shear flows; poor for separation
• k-omega SST is the workhorse for external aerodynamics and separated flows
• Wall treatment (y+) is critical for accurate results
• LES resolves large eddies directly — expensive but more accurate
• Model selection depends on flow physics and available resources

What's Next

Turbulence models introduce approximations, and so does every other aspect of CFD. The next lesson covers Verification & Validation — how to quantify numerical errors and build confidence in your results.