Turbulence Runs Three Rulers at Once — The Kolmogorov Cascade and the -5/3 Law

1941, a notebook in Moscow#

Two months before the war reached his door, Kolmogorov mailed three short notes to a journal. Together they ran fewer than sixteen pages. There was no wind tunnel, no measurement data, no figure. Only dimensional analysis and a pair of statistical assumptions.

One line from that notebook still pins CFD down eighty years later. At sufficiently high Reynolds number, the way energy is distributed across scales is fixed by a single length — the Kolmogorov microscale $\eta$. This post walks through how that line arises, why it sets your grid budget and DNS (Direct Numerical Simulation) cost, and finishes by fitting a -5/3 slope to a synthetic spectrum.

Raise the time-speed slider and the gap between breakups shrinks, while the order from large eddy to small stays put.

Large eddies (cyan) live a long turnover, $\tau \sim L^{2/3}$. Medium eddies (orange) split faster, and small eddies (pink) hand off to gray (viscous dissipation) almost instantly. While each eddy is alive, the energy flux it passes down equals every other scale's flux — that is the definition of the inertial subrange.

Three rulers stuffed into one flow#

Lay out a chunk of turbulence and three lengths show up.

• Integral scale $L$ — the size of the largest eddies. Set by the flow's geometry (pipe diameter, chimney width). This is where the energy sits.
• Taylor microscale $\lambda$ — defined through the curvature of the two-point velocity autocorrelation at the origin. A mediator between large and small.
• Kolmogorov microscale $\eta$ — the smallest ruler, where viscosity grinds kinetic energy into heat. Far above the molecular mean free path, but far below $L$.

The three rulers coexist in the same flow. As the Reynolds number climbs, the gap between $L$ and $\eta$ widens. That gap is the inertial subrange.

Dimensional analysis bridges the gap#

Kolmogorov's hypotheses fit on two lines.

1. At small enough scales, statistics depend only on viscosity $\nu$ and the energy dissipation rate $\varepsilon$.
2. In the inertial subrange viscosity drops out, and statistics depend only on $\varepsilon$.

Hypothesis 1 builds the smallest scale. $\nu$ has units $\mathrm{m^2/s}$ and $\varepsilon$ has units $\mathrm{m^2/s^3}$. To extract a length:

$\eta$ is length, $u_\eta$ velocity, $\tau_\eta$ time. Units alone close the system.

Hypothesis 2 cuts deeper. The inertial-range spectrum $E(k)$ must depend only on $\varepsilon$ and wavenumber $k$. $E$ has units $\mathrm{m^3/s^2}$. Matching dimensions:

$C$ is the Kolmogorov constant, measured around $C \approx 1.5$. The -5/3 law falls out from assumptions and units. Not a single data point was used.

Scale separation — the gap Re opens#

The large-eddy estimate $\varepsilon \sim u'^3 / L$ (energy injection equals dissipation in steady state) gives:

with $\mathrm{Re} = u'L/\nu$. A tenfold rise in Reynolds shrinks $\eta/L$ to $10^{-3/4} \approx 0.18$ of its previous value. For a grid to resolve every scale, cell count scales as $\mathrm{Re}^{9/4}$. A laboratory wind tunnel at $\mathrm{Re} = 10^4$ has $L/\eta \approx 1{,}000$, and a 3D DNS needs $\sim 10^9$ cells. An aircraft wing at $\mathrm{Re} = 10^7$ would demand $10^{15}$ cells — beyond any current supercomputer.

The fit window's width in decades is fixed by its edges. From $10/L$ to $0.1/\eta$ the width is $\log_{10}(L/\eta)-2$. At $\mathrm{Re}=10^{4}$, $L/\eta\approx 1{,}000$, so the window is one decade. Each tenfold rise in Reynolds multiplies $\eta/L$ by $0.18$. Apply that twice, at $\mathrm{Re}=10^{6}$, and $\eta/L$ is about $0.032$ of the old value, so $L/\eta$ opens to about $3\times 10^{4}$ and the window to about 2.5 decades.

With the snippet's $u'=1$ and $L=1$, $\varepsilon=u'^{3}/L=1$. Before $f_L$ and $f_\eta$, the inertial piece is $E(k)=1.5\,k^{-5/3}$. The commented fit, −1.661, sits 0.006 off the expected −1.667. Ten percent noise still leaves three digits if the window stays shut.

How that gap first opens inside a pipe continues in the critical Reynolds where a water jet breaks into layers. As the nominal Reynolds grows, cells between $L$ and $\eta$ scale as $\mathrm{Re}^{9/4}$, and at $10^{4}$ a 3D grid is already $10^{9}$.

The road energy travels#

A large eddy spins. Viscosity is too weak to drain it directly. Instead the large eddy breaks into smaller eddies, those into smaller still. The breakup time is the eddy turnover time $\tau_\ell \sim \ell/u_\ell$. That breakup binds every scale through a single number — the dissipation rate $\varepsilon$.

The -5/3 slope and the traps at either end#

Plot $E(k)$ log-log and a straight band at slope -5/3 appears in the middle. Different physics rule the two ends.

• Around $k \sim 1/L$ — the energy-containing range. The geometry-imposed largest eddies live here. $E(k)$ peaks.
• Around $k \sim 1/\eta$ — the viscous dissipation range. $E(k)$ falls exponentially. Viscosity drinks the energy.

Move the Reynolds-number slider up and the -5/3 band gets longer.

At $\mathrm{Re} = 10^3$ the -5/3 band is barely visible. At $\mathrm{Re} = 10^6$ a nearly four-decade straight line opens. Experimentally, if a clean -5/3 slope holds over even a single decade, the flow counts as fully developed turbulence.

Python — fitting -5/3 to a synthetic spectrum#

How does one extract a -5/3 slope from real data? Linear regression on log axes. Let's try on synthetic data.

import numpy as np
 
C_K = 1.5
c_L, p0 = 6.78, 2
beta_e, c_eta = 5.2, 0.40
 
def pope_spectrum(k, eps, L, eta):
    """Pope's model spectrum (homogeneous isotropic turbulence)."""
    kL = k * L
    kEta = k * eta
    f_L = (kL / np.sqrt(kL**2 + c_L)) ** (5/3 + p0)
    f_eta = np.exp(-beta_e * ((kEta**4 + c_eta**4)**0.25 - c_eta))
    return C_K * eps**(2/3) * k**(-5/3) * f_L * f_eta
 
# Setup
L_int = 1.0      # integral length
u_rms = 1.0      # rms velocity
Re = 1e5
nu = u_rms * L_int / Re
eps = u_rms**3 / L_int
eta = (nu**3 / eps) ** 0.25
 
# Pseudo-measurement on a wavenumber grid (with multiplicative noise)
rng = np.random.default_rng(0)
k = np.logspace(-1, np.log10(1.0 / eta) - 0.5, 80)
E_true = pope_spectrum(k, eps, L_int, eta)
E_meas = E_true * np.exp(0.10 * rng.standard_normal(k.size))   # ~10% lognormal noise
 
# Fit -5/3 only on the inertial range  (10/L < k < 0.1/eta)
mask = (k > 10 / L_int) & (k < 0.1 / eta)
slope, _ = np.polyfit(np.log(k[mask]), np.log(E_meas[mask]), 1)
print(f"L/eta = {L_int / eta:.1f}")
print(f"fitted slope = {slope:+.3f}  (expected -1.667)")
# L/eta = 17783.3
# fitted slope = -1.661  (expected -1.667)

Even with 80 noisy points the fitted slope is -1.66, three digits in agreement with the theoretical -5/3. Crank the noise up to 30% and the slope still lands between -1.59 and -1.74 — the inertial range is sturdy.

One question. Why mask out $k < 10/L$ and $k > 0.1/\eta$? Because $f_L$ and $f_\eta$ bend the curve at both ends. In real measurements, a fit window that is too wide drags the slope toward -1.4. When a colleague reports "my data doesn't show -5/3", half the time the culprit is the fit window.

Surviving where the grid can't reach#

Once the cost of DNS is laid out, the workarounds become legible.

• RANS (Reynolds-Averaged Navier–Stokes) — compute only the time-averaged field and close the fluctuations with a model. $\eta$ is not resolved, so grid cost scales like $\mathrm{Re}^0$. The industry workhorse, weak on anisotropic Reynolds stresses.
• LES (Large Eddy Simulation) — resolve the large eddies and close subgrid eddies through an SGS (subgrid-scale) model. Grid cost ranges from $\mathrm{Re}^{0.5}$ to $\mathrm{Re}^1$ depending on whether walls are modeled. The middle ground.
• DNS — resolve down to $\eta$ and integrate down to $\tau_\eta$. Cost scales as $\mathrm{Re}^{3}$ (3D plus time). An academic tool.

All three rest on the same -5/3 line. RANS replaces the entire inertial range with a model, LES resolves only part of it, DNS resolves all of it. Where you cut is the modeling decision.

When the fit window opens to both ends and the slope sags to −1.4#

$\eta/L = \mathrm{Re}^{-3/4}$ sets how wide the window can be. Leave the mask outside $10/L < k < 0.1/\eta$ and $f_L$ plus $f_\eta$ bend the ends, so the slope that should have followed only $\varepsilon$ sags toward −1.4. The −1.661 in the code comment is what that closed window keeps.

Even $\mathrm{Re} = 10^4$ is already a heavy 3D grid, and $\mathrm{Re} = 10^7$ pushes it to $10^{15}$. Cut below that and hand the rest to a subgrid model, and the place where the coefficient comes out negative is written in where the dynamic coefficient went negative under Germano's identity. Trusting the cascade flux and shipping the clipped coefficient unchanged are different decisions.