Fluid Dynamics for CFD Engineers
Navier-Stokes to Turbulence Modeling

Category: Fundamental Theory | Updated: 2026-03-25 | NovaSolver Contributors

Whether you're designing a cooling channel for a battery pack, predicting drag on a car, or analyzing a supersonic nozzle, the same governing equations are at work. Understanding the physics behind those equations — and the approximations your CFD software makes to solve them — is the key to getting reliable results.

1. Fluid Properties

🧑‍🎓

Before we get into the equations, what are the key fluid properties I need to define in a CFD solver? I always see density, viscosity, and sometimes thermal conductivity — but I'm not sure which ones actually matter for which problems.

🎓

Simply put: density drives inertia (how the fluid resists acceleration), and dynamic viscosity drives internal friction (resistance to shearing). For isothermal incompressible flow — water through a pipe, say — those two are enough. Add thermal conductivity and specific heat once you care about temperature. Add equation of state (link density to pressure/temperature) once the flow is compressible. So the physics of your problem dictates exactly which properties you need.

1.1 Key Fluid Properties

Property	Symbol	Units (SI)	Water (20°C)	Air (20°C, 1 atm)
Density	ρ	kg/m³	998	1.204
Dynamic viscosity	μ	Pa·s = kg/(m·s)	1.002×10⁻³	1.825×10⁻⁵
Kinematic viscosity	ν = μ/ρ	m²/s	1.004×10⁻⁶	1.516×10⁻⁵
Thermal conductivity	k	W/(m·K)	0.598	0.0257
Specific heat (const. pressure)	c_p	J/(kg·K)	4182	1005
Prandtl number	Pr = μc_p/k	—	7.01	0.713
Speed of sound	a	m/s	~1482	343

Newtonian vs. non-Newtonian fluids: Most engineering fluids (water, air, oils) are Newtonian — dynamic viscosity is constant with shear rate. Non-Newtonian fluids (blood, paint, polymer melts, concrete slurry) have viscosity that depends on shear rate. Modeling these requires specialized material models (power-law, Bingham plastic, Carreau model) beyond standard CFD.

2. Continuity Equation — Conservation of Mass

🧑‍🎓

I've seen "divergence-free velocity field" mentioned for incompressible flow. What does that actually mean physically?

🎓

Imagine a small control volume in the fluid. "Divergence-free" means exactly as much fluid enters as leaves — no accumulation, no source. That's just mass conservation for an incompressible fluid. If you compress a gas (compressible flow), density can change and the equation has an extra term. In practice, if the Mach number is below about 0.3, you can safely treat the flow as incompressible — which makes the math much simpler.

The general continuity equation (conservation of mass) in differential form:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0$$

For incompressible flow (ρ = constant):

$$\nabla \cdot \mathbf{u} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0$$

This constraint is what makes incompressible CFD tricky — pressure and velocity are coupled through the divergence-free condition, not through an explicit pressure equation. The SIMPLE and projection algorithms exist precisely to handle this coupling.

3. The Navier-Stokes Equations

🧑‍🎓

I know Navier-Stokes equations are "the equations of fluid motion," but I've never really understood each term. Can you break it down in plain English?

🎓

Think of it as Newton's second law — F = ma — written for a fluid parcel. Left side is acceleration (two terms: unsteady change and convective change as the parcel moves to a different location). Right side is the sum of forces: pressure gradient (pushes fluid from high to low pressure), viscous stresses (internal friction), and body forces like gravity or buoyancy. That's literally all there is. The mathematical complexity comes from the nonlinear convective term, which is the root cause of turbulence.

3.1 Incompressible Navier-Stokes (Newtonian Fluid)

$$\rho\underbrace{\left(\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}\right)}_{\text{inertia}} = \underbrace{-\nabla p}_{\text{pressure}} + \underbrace{\mu \nabla^2 \mathbf{u}}_{\text{viscous}} + \underbrace{\rho \mathbf{g}}_{\text{body force}}$$

Expanded in index notation (Einstein summation convention):

$$\rho\left(\frac{\partial u_i}{\partial t} + u_j\frac{\partial u_i}{\partial x_j}\right) = -\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j}\left[\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right] + \rho g_i$$

The nonlinear convective term $(\mathbf{u}\cdot\nabla)\mathbf{u}$ — velocity times velocity-gradient — is responsible for the enormous mathematical difficulty. It's this term that drives the cascade of energy from large eddies to small ones in turbulent flow.

3.2 Dimensionless Navier-Stokes

Non-dimensionalizing with reference velocity $U_0$, length $L$, pressure $\rho U_0^2$:

$$\frac{\partial \mathbf{u}^*}{\partial t^*} + (\mathbf{u}^*\cdot\nabla^*)\mathbf{u}^* = -\nabla^* p^* + \frac{1}{Re}\nabla^{*2}\mathbf{u}^*$$

where $Re = \rho U_0 L / \mu$ is the Reynolds number. This shows that two geometrically similar flows at the same Re are physically equivalent — the foundation of wind tunnel testing and scale-model CFD validation.

3.3 Stokes Flow (Re → 0)

When inertia is negligible (microfluidics, lubrication), the convective term drops out:

$$\nabla p = \mu \nabla^2 \mathbf{u}, \qquad \nabla \cdot \mathbf{u} = 0$$

These linear equations are much easier to solve analytically. Stokes flow applies in: bearing lubrication films (very thin gaps, high viscosity), microfluidic channels, slow groundwater seepage.

4. Energy Equation

The transport equation for total energy (or temperature, for low-Mach incompressible flows with Boussinesq approximation):

$$\rho c_p\left(\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T\right) = \nabla\cdot(k\nabla T) + \dot{q} + \Phi$$

where $\dot{q}$ is volumetric heat generation (W/m³) and $\Phi$ is the viscous dissipation function (important in high-speed or high-viscosity flows):

$$\Phi = \mu\left[2\left(\frac{\partial u}{\partial x}\right)^2 + 2\left(\frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2\right]$$

🧑‍🎓

I'm doing a simulation of air cooling for electronics. Do I need to solve the full energy equation, or is there a simpler approach?

🎓

For low-speed air cooling (Mach << 1), the Boussinesq approximation works well — you solve the incompressible Navier-Stokes plus the energy equation, but with one modification: density changes only appear in the gravity term as buoyancy: $\rho \mathbf{g} \approx \rho_0[1 - \beta(T-T_0)]\mathbf{g}$. This captures natural convection without solving the full compressible equations. Viscous dissipation is also negligible for air cooling. All major CFD solvers (Fluent, OpenFOAM, STAR-CCM+) implement this by default when you enable the energy equation with the Boussinesq option.

5. Reynolds Number and Flow Regimes

$$Re = \frac{\rho U L}{\mu} = \frac{UL}{\nu} = \frac{\text{inertial forces}}{\text{viscous forces}}$$

The Reynolds number is the single most important dimensionless parameter in fluid dynamics. It determines whether flow is laminar or turbulent:

Re Range	Flow Regime	Typical Application
Re < 2,300	Laminar (pipe flow)	Microfluidics, viscous oil flows, human blood vessels
2,300 < Re < 4,000	Transitional	Critical design region — both regimes possible, unstable
Re > 4,000	Turbulent (pipe flow)	Most industrial pipe flows (water, air in ducts)
Re > 5×10⁵	Turbulent (flat plate)	Aircraft boundary layers, automotive external flow

🧑‍🎓

So for flow around a car at highway speed, what's the Reynolds number roughly?

🎓

Let's calculate: 120 km/h = 33 m/s, car length L ≈ 4.5 m, air kinematic viscosity ν ≈ 1.5×10⁻⁵ m²/s. So Re = 33 × 4.5 / 1.5×10⁻⁵ ≈ 10⁷. That's highly turbulent — you absolutely need a turbulence model. Compare that to a blood vessel: U ≈ 0.1 m/s, D ≈ 3 mm = 0.003 m, blood ν ≈ 3×10⁻⁶ m²/s, so Re ≈ 100 — nearly laminar, no turbulence model needed. The same equation, completely different physics.

6. Turbulence Fundamentals and RANS Models

🧑‍🎓

In Fluent I see options like k-epsilon, k-omega SST, RSM, LES... It's overwhelming. Can you explain what turbulence modeling actually is and why we need it?

🎓

Here's the core problem: turbulent flow has eddies at scales from meters down to micrometers. Resolving all of them directly (DNS — Direct Numerical Simulation) requires a mesh with billions of cells and a supercomputer running for weeks. So instead, we time-average the equations (RANS — Reynolds-Averaged Navier-Stokes) and model the effect of turbulence on the mean flow. That modeling introduces additional unknowns — the Reynolds stresses — that we need extra equations to close. Those extra equations are the turbulence models. Think of it as parametrizing the chaos rather than resolving it.

6.1 Reynolds Decomposition and RANS

Every variable is split into mean and fluctuating parts: $\mathbf{u} = \overline{\mathbf{u}} + \mathbf{u}'$. After time-averaging the N-S equations:

$$\rho\left(\frac{\partial \overline{u}_i}{\partial t} + \overline{u}_j\frac{\partial \overline{u}_i}{\partial x_j}\right) = -\frac{\partial \overline{p}}{\partial x_i} + \frac{\partial}{\partial x_j}\left[\mu\frac{\partial \overline{u}_i}{\partial x_j} - \rho\overline{u_i' u_j'}\right]$$

The term $-\rho\overline{u_i' u_j'}$ is the Reynolds stress tensor — 6 additional unknowns that require closure equations. This is the "turbulence closure problem."

6.2 The Boussinesq Turbulent Viscosity Hypothesis

The most common closure: model Reynolds stresses analogously to viscous stresses:

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

where $\mu_t$ is the turbulent (eddy) viscosity and $k = \frac{1}{2}\overline{u_i'u_i'}$ is the turbulent kinetic energy. The task of the turbulence model is to compute $\mu_t$.

6.3 Standard k-ε Model

Two transport equations for turbulent kinetic energy $k$ and dissipation rate $\varepsilon$:

$$\frac{\partial(\rho k)}{\partial t} + \nabla\cdot(\rho \mathbf{u} k) = \nabla\cdot\!\left[\left(\mu+\frac{\mu_t}{\sigma_k}\right)\nabla k\right] + P_k - \rho\varepsilon$$ $$\frac{\partial(\rho\varepsilon)}{\partial t} + \nabla\cdot(\rho \mathbf{u} \varepsilon) = \nabla\cdot\!\left[\left(\mu+\frac{\mu_t}{\sigma_\varepsilon}\right)\nabla\varepsilon\right] + C_{\varepsilon 1}\frac{\varepsilon}{k}P_k - C_{\varepsilon 2}\frac{\rho\varepsilon^2}{k}$$

The turbulent viscosity: $\mu_t = C_\mu \rho k^2/\varepsilon$

Standard constants: $C_\mu = 0.09$, $C_{\varepsilon 1} = 1.44$, $C_{\varepsilon 2} = 1.92$, $\sigma_k = 1.0$, $\sigma_\varepsilon = 1.3$.

Strength: Robust, computationally cheap, widely validated for free shear flows (jets, wakes). Weakness: Overestimates turbulent viscosity in adverse pressure gradient flows (causes boundary layers to separate too late), poor near-wall behavior.

6.4 k-ω SST Model (Menter 1994)

The Shear Stress Transport model blends k-ε (in the freestream) with k-ω (near walls) using a blending function $F_1$:

$$\frac{\partial(\rho k)}{\partial t} + \nabla\cdot(\rho\mathbf{u}k) = \nabla\cdot[(\mu+\sigma_k\mu_t)\nabla k] + \tilde{P}_k - \beta^*\rho k\omega$$ $$\frac{\partial(\rho\omega)}{\partial t} + \nabla\cdot(\rho\mathbf{u}\omega) = \nabla\cdot[(\mu+\sigma_\omega\mu_t)\nabla\omega] + \alpha\frac{\rho\omega}{k}\tilde{P}_k - \beta\rho\omega^2 + 2(1-F_1)\frac{\rho\sigma_{\omega 2}}{\omega}\nabla k\cdot\nabla\omega$$

The key feature: turbulent viscosity is limited using the SST modification:

$$\mu_t = \frac{\rho a_1 k}{\max(a_1\omega,\, SF_2)}$$

where $S$ is the strain rate magnitude and $F_2$ is another blending function. This prevents over-prediction of eddy viscosity in adverse pressure gradients. k-ω SST is the industry default for external aerodynamics, turbomachinery, and any separated flow.

6.5 LES and Hybrid Methods

Method	What it resolves	Cost (relative to DNS)	Best use case
DNS	All scales	1 (reference)	Academic research, Re < 10,000
LES	Large eddies; models small	~1/100 of DNS	Noise prediction, combustion, LES-required physics
DES/IDDES	RANS near walls + LES far field	~1/1000 of DNS	Massively separated flows, bluff bodies
RANS (k-ω SST)	None (all modeled)	~1/10⁶ of DNS	Industrial design, parametric studies

LES filters the N-S equations: $\tilde{u}_i = \int G(\mathbf{x},\mathbf{x}')u_i(\mathbf{x}')d\mathbf{x}'$. The sub-grid scale stresses require a closure model (Smagorinsky: $\tau_{ij}^{SGS} = -2C_s^2\bar{\Delta}^2|\tilde{S}|\tilde{S}_{ij}$).

🧑‍🎓

So for a typical industrial problem — say, predicting drag on a wing mirror — which model would you recommend?

🎓

For a wing mirror specifically — a bluff body with massive separation — k-ω SST gives reasonable drag and pressure distribution, but it tends to under-predict the noise because it doesn't capture the unsteady vortex shedding. In production automotive aero, companies typically use DDES (Delayed Detached Eddy Simulation) for the detailed mirror study, but start with RANS SST for screening dozens of geometries quickly. The rule of thumb: RANS for design iteration, LES/DES for final validation of important configurations.

7. Compressible Flow and Shock Waves

When the Mach number exceeds ~0.3, density variations become significant. The compressible N-S equations add an equation of state and couple velocity, pressure, density, and temperature:

$$\frac{\partial \rho}{\partial t} + \nabla\cdot(\rho\mathbf{u}) = 0$$ $$\frac{\partial(\rho\mathbf{u})}{\partial t} + \nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u}) = -\nabla p + \nabla\cdot\boldsymbol{\tau}$$ $$\frac{\partial(\rho E)}{\partial t} + \nabla\cdot[(\rho E + p)\mathbf{u}] = \nabla\cdot(\boldsymbol{\tau}\cdot\mathbf{u} - \mathbf{q})$$ $$p = \rho R T \quad \text{(ideal gas)}$$

Key dimensionless parameters for compressible flow:

$$Ma = \frac{U}{a} = \frac{U}{\sqrt{\gamma RT}}, \qquad \gamma = \frac{c_p}{c_v}$$

Mach regime	Ma range	Characteristics
Subsonic	Ma < 1	Smooth flow, pressure waves travel upstream
Transonic	0.8 < Ma < 1.2	Mixed sub/supersonic zones, shock waves form
Supersonic	1 < Ma < 5	Oblique shocks, expansion fans, no upstream influence
Hypersonic	Ma > 5	High-temperature dissociation effects, intense heating

7.1 Normal Shock Relations

Across a normal shock, Rankine-Hugoniot conditions apply:

$$Ma_2^2 = \frac{(\gamma-1)Ma_1^2 + 2}{2\gamma Ma_1^2 - (\gamma-1)}, \qquad \frac{p_2}{p_1} = 1 + \frac{2\gamma}{\gamma+1}(Ma_1^2-1)$$

For a normal shock at Ma = 2 in air (γ = 1.4): Ma₂ = 0.577, p₂/p₁ = 4.5. Total pressure is lost across a shock — a key consideration in supersonic inlet design.

8. Numerical CFD: FVM, SIMPLE, and Convection Schemes

🧑‍🎓

I know FEM for structures, but CFD mostly uses FVM — Finite Volume Method. Why the different approach?

🎓

In practice, FVM is preferred in CFD because it directly discretizes conservation laws — mass, momentum, energy — in integral form over control volumes. This gives exact conservation by construction on any mesh. FEM can achieve this too, but historically FVM was simpler to implement conservatively for convection-dominated problems and unstructured meshes. Ansys Fluent, OpenFOAM, and STAR-CCM+ are all FVM-based.

8.1 Finite Volume Discretization

Integrating the continuity equation over a control volume $V$ with surface $S$:

$$\frac{d}{dt}\int_V \rho\, dV + \oint_S \rho\mathbf{u}\cdot d\mathbf{S} = 0$$

Discretized over cell $P$ with neighbor $N$ sharing face $f$:

$$\frac{\rho_P^{n+1} - \rho_P^n}{\Delta t} V_P + \sum_f (\rho\mathbf{u})_f \cdot \mathbf{S}_f = 0$$

8.2 Pressure-Velocity Coupling: SIMPLE Algorithm

For incompressible flow, the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm:

Guess pressure field $p^*$
Solve momentum equations → velocity $\mathbf{u}^*$ (doesn't satisfy continuity)
Compute pressure correction $p' = p - p^*$ from continuity
Correct velocity: $\mathbf{u} = \mathbf{u}^* + f(p')$
Update scalars (T, k, ε, ...)
Check convergence; repeat from step 2

8.3 Convection Scheme Selection

Scheme	Order	Bounded?	Recommended for
Upwind (UDS)	1st	Yes	Initial convergence, highly skewed meshes
Linear (CDS)	2nd	No	Low-Re smooth flows
QUICK	3rd	Partial	Structured hexahedral meshes
Linear-upwind (LUDS)	2nd	Partial	OpenFOAM default, general purpose
MUSCL	3rd	Yes (with limiter)	Compressible flows, shocks

The false diffusion trap: First-order upwind smears sharp gradients (like a temperature front in a jet) due to numerical diffusion. Always switch to second-order after initial convergence in RANS simulations.

9. Practical CFD Pitfalls

🧑‍🎓

I keep getting "turbulent kinetic energy residuals diverging" when I start my simulation. What does that usually mean?

🎓

That's almost always bad initial conditions or bad inlet turbulence boundary conditions. The k and ω (or ε) equations are stiff — if you initialize k=0 exactly, or set it to an unrealistically high value, the production-destruction balance in the turbulence equations goes haywire. Start with turbulence intensity of 5% and a turbulent length scale of 10% of the hydraulic diameter — that's a safe, physically reasonable initial condition. Also, always start with first-order upwind and low relaxation factors (0.3 for pressure, 0.7 for velocity), then switch to second-order once residuals drop two orders of magnitude.

9.1 Near-Wall Treatment and y⁺

The dimensionless wall distance:

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

Wall treatment	Required y⁺	Turbulence model compatibility
Low-Re (resolve sublayer)	y⁺ ≈ 1	k-ω SST (recommended), k-ε with two-layer
Standard wall function	30 < y⁺ < 300	k-ε standard, realizable k-ε
Enhanced wall function	1 < y⁺ < 300	Fluent enhanced, OpenFOAM nutWallFunction

9.2 Common CFD Mistakes Checklist

Inlet turbulence conditions unspecified — solver uses defaults that may not match your setup
Domain too small — blockage ratio > 5% in external flows distorts results; outlet should be 15+ diameters downstream
Only first-order schemes — gives smooth convergence but wrong answer due to numerical diffusion
Residuals as only convergence criterion — always also monitor a physical quantity (drag coefficient, outlet temperature) that should plateau
Ignoring mesh quality — max aspect ratio, skewness, and orthogonality matter; 90°+ skewness cells cause divergence
Wrong turbulence model for separated flow — k-ε overestimates recovery in separated regions; use SST or DES instead

Key Takeaways

Navier-Stokes = Newton's 2nd law for fluids: inertia = pressure + viscous + body forces
Reynolds number determines laminar vs. turbulent; most industrial flows are turbulent
RANS models (k-ω SST recommended) time-average and model turbulence effects
FVM discretizes conservation laws in integral form — conservative by construction
y⁺ ≈ 1 for SST; 30–300 for wall functions
Always monitor physical quantities, not just residuals, to confirm convergence

Related Interactive Tools

Reynolds Number Calculator — Compute Re and visualize flow regime transitions
y⁺ First Cell Height Calculator — Given Re and turbulence model, get Δs
Flat Plate Boundary Layer — Visualize δ(x) growth, laminar vs. turbulent
Pipe Flow Pressure Loss — Darcy-Weisbach + Moody chart

Cross-Topics

Heat Transfer Numerical Methods CFD Analysis OpenFOAM Guide Ansys Fluent Guide Coupled Analysis

Written by NovaSolver Contributors (Anonymous Engineers & AI) | CAE Technical Encyclopedia

Fluid Dynamics for CFD EngineersNavier-Stokes to Turbulence Modeling