Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section, you will be able to:

Explain the qualitative difference between laminar and turbulent flow in plain language, with at least two everyday examples.
Derive the Reynolds number $Re = \rho U L / \mu$ from dimensional analysis of the inertial and viscous terms in Navier-Stokes.
Predict whether a given flow is laminar, transitional, or turbulent by computing one number.
Sketch the laminar parabola and the turbulent plug-with-thin-boundary-layer profile and explain the physics behind each shape.
Estimate the Kolmogorov microscale and connect it to the energy cascade.
Implement Reynolds decomposition in PyTorch and differentiate the turbulent kinetic energy with autograd.

The Two Personalities of a Fluid

Light a candle. Look at the column of smoke rising from the flame. Right above the wick the smoke goes up in a clean, glassy, almost-still ribbon. A few centimetres higher, that same ribbon shatters into swirls, billows, and chaotic curls. Two regimes of flow, an inch apart, in the same air.

The lower ribbon is laminar: fluid layers slide past each other like sheets of glass, neighbouring particles stay neighbours, dye dropped in stays in a thin line. The upper curl is turbulent: layers tear, vortices spawn and die, dye is shredded and mixed within seconds. The Navier-Stokes equations govern both pictures — the difference is not in the law, it is in the regime.

One-sentence slogan. A laminar flow is orderly time-marching layers; a turbulent flow is chaotic three-dimensional mixing. The Reynolds number tells you which one you have.

🟦 Laminar — “sheet flow”

Smooth, predictable streamlines
Dye stays in a thin line
Velocity profile is a tall parabola in a pipe
Friction $f = 64/Re$ (exact)
Reversible in time — run the experiment backwards and it works

🟥 Turbulent — “chaotic mixing”

Eddies of every size, never repeating
Dye is shredded and mixed in seconds
Velocity profile is a flat plug with thin boundary layers
Friction $f \sim 0.316 / Re^{1/4}$ (empirical)
Irreversible — kinetic energy cascades down to heat

The honey-versus-water intuition

Pour honey from a spoon. The stream is laminar — slow, viscous, impossible to make it splash. Now pour water from the same height. Water hits the cup and erupts into a chaotic crown of droplets. Same gravity, same height, completely different regime. The only difference is the ratio between how much momentum the fluid carries and how strongly viscosity wants to smooth it out. That ratio has a name. It is the Reynolds number.

Reynolds' Dye Experiment (1883)

In 1883 Osborne Reynolds put a thin streak of coloured dye into water flowing through a glass pipe and slowly increased the flow rate. He photographed three distinct stages of behaviour:

Slow flow — the dye stream stayed perfectly straight all the way down the pipe.
Medium flow — the dye began to oscillate, then spontaneously burst into intermittent puffs of mixed colour.
Fast flow — the dye smeared into a uniform pink mist within the first centimetre of the pipe.

Reynolds discovered that the transition between stages 1 and 3 did not happen at a particular velocity, viscosity, pipe size, or density — but at a particular value of the dimensionless combination $Re = \rho U D / \mu$ . He found $Re \approx 2300$ as the lower transition threshold. The number has stood, with tiny refinements, for 140 years.

What Reynolds really showed

The Navier-Stokes equations contain a parameter (the viscosity $\mu$ ) and the boundary data contain three more ( $\rho, U, L$ ). Reynolds proved by experiment what dimensional analysis suggests: only one dimensionless combination of these four numbers matters. The flow patterns are universal in $Re$ . Two flows with the same $Re$ look the same — even if one is air in a wind tunnel and the other is oil in a pipeline.

Deriving the Reynolds Number From Scratch

Where does $Re = \rho U L / \mu$ come from? Pure dimensional analysis. We start from the incompressible Navier-Stokes momentum equation and ask: which terms balance, and at what scale?

Step 1: write down Navier-Stokes

\displaystyle \rho \!\left(\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v}\!\cdot\!\nabla)\mathbf{v}\right) = -\nabla p + \mu\,\nabla^2 \mathbf{v}

The left-hand side is inertia: how the fluid's momentum changes. The right-hand side has two forces: pressure and viscous diffusion. To pick a regime we only need to compare two of these terms — the convective inertial term $\rho (\mathbf{v}\!\cdot\!\nabla)\mathbf{v}$ and the viscous term $\mu \nabla^2 \mathbf{v}$ .

Step 2: estimate each term by orders of magnitude

Let $U$ be a representative speed and $L$ a representative length. Then:

\displaystyle \underbrace{\rho \,(\mathbf{v}\!\cdot\!\nabla)\mathbf{v}}_{\text{inertia}} \;\sim\; \frac{\rho\, U^2}{L} \qquad \text{vs.} \qquad \underbrace{\mu\, \nabla^2 \mathbf{v}}_{\text{viscous}} \;\sim\; \frac{\mu\, U}{L^2}

The inertial term has velocity twice (one factor of $\mathbf{v}$ , one factor of $\nabla \mathbf{v}$ ), each scaling like $U/L$ . The viscous term has the Laplacian $\nabla^2$ , which contributes $U/L^2$ . The ratio is what defines a regime.

Step 3: form the dimensionless ratio

\displaystyle \boxed{\;Re \;=\; \frac{\text{inertia}}{\text{viscosity}} \;\sim\; \frac{\rho U^2 / L}{\mu U / L^2} \;=\; \frac{\rho\, U\, L}{\mu}\;}

Three multiplications and one division — that is the entire derivation. Notice what each ingredient means in the ratio:

Symbol	Plain meaning	Effect on Re
ρ	How heavy the fluid is	Heavier → more inertia → higher Re
U	How fast the flow moves	Faster → more inertia → higher Re
L	How big the geometry is	Bigger → more room for inertia → higher Re
μ	How viscous (sticky) the fluid is	Stickier → more damping → lower Re

A second viewpoint — time scales

You can derive the same number by comparing two time scales: the convective time $L/U$ over which a particle crosses the geometry, and the diffusive time $\rho L^2 / \mu$ for viscosity to smooth out a velocity bump of size $L$ . Their ratio is again $Re$ . Inertia wins when it acts faster than viscosity can repair the disturbance.

Why a Single Number Decides Everything

Here is the deep idea. Both shear instability (which amplifies tiny disturbances into eddies) and viscous damping (which kills them) compete inside the same Navier-Stokes equation. The Reynolds number is the dimensionless switch that tells you which one wins.

Re range	Regime	What you see
Re < 1	Stokes (creeping) flow	Bacteria swim, paint flows down a wall
1 ≤ Re < ~2300	Laminar pipe flow	Smooth profile, dye stays in line
2300 ≤ Re < ~4000	Transitional	Intermittent puffs of turbulence
~4000 ≤ Re < ~1e5	Fully developed turbulence	Chaotic eddies, mixed dye
Re ≳ 1e5	High-Re turbulence	Cascade of eddies, near-wall boundary layer

The same number for a bee and a Boeing

A bee in still air has $Re \sim 10$ — its world is dominated by viscosity. A Boeing 747 cruising at altitude has $Re \sim 10^9$ — its world is dominated by inertia. The difference is not 10 orders of magnitude of unrelated physics. It is 10 orders of magnitude of this one number.

Interactive Explorer: Drag the Re Slider

The viewer below shows hundreds of tracer particles drifting through a pipe whose mean velocity profile follows the regime you select. Drag the Reynolds-number slider and watch three things change at once:

The particles. Below Re ≈ 2 300 they glide in smooth horizontal layers. Above 4 000 they erupt into cross-stream chaos.
The mean profile. The blue parabola sharpens to a point in laminar flow. The red plug flattens out, with a sudden drop to zero at the walls in turbulent flow.
The friction factor. Watch the readout drop fast in laminar flow (~1/Re) and slow in turbulent (~Re^-1/4) — this is the calling card of the cascade.

Loading interactive Reynolds-number explorer...

3Blue1Brown-style challenge: click the Transitional (Re ≈ 2 800) preset and watch carefully. The particles oscillate between long laminar streaks and sudden bursts of mixing. That is exactly what Reynolds saw in his 1883 dye experiment — intermittency, the signature of the transition regime.

Laminar Parabola vs Turbulent Plug

Two regimes, two utterly different velocity profiles inside a pipe. Both satisfy the no-slip BC at the wall — but the shape between the walls could not be more different.

The laminar parabola (Hagen-Poiseuille)

For steady fully developed flow in a circular pipe of diameter $D$ , the Navier-Stokes equations collapse to a 1D ODE whose solution is a perfect parabola:

\displaystyle u(r) = u_{\max}\!\left(1 - \left(\tfrac{2r}{D}\right)^2\right), \quad u_{\max} = 2\, U_{\text{avg}}

The centreline moves at exactly twice the average — a steep gradient throughout. Friction is large because the velocity has to bend from $u_{\max}$ in the middle all the way to zero at the wall, smoothly.

The turbulent plug (1/7 power law)

For fully developed turbulent pipe flow, the empirical fit is the 1/7-power profile:

\displaystyle u(r) = u_{cl}\,\Bigl(1 - \tfrac{2|r|}{D}\Bigr)^{\!1/7}, \quad u_{cl} \approx 1.224\, U_{\text{avg}}

The centreline only beats the average by ~22 %, because the bulk of the pipe is moving at nearly the same speed. The drop to zero happens in a razor-thin layer right next to the wall — the famous viscous sublayer. The shear in that layer is enormous, which is why turbulent friction is large in absolute terms even though the bulk profile looks flat.

Why turbulence flattens the profile

Eddies are nature's mixers. They scoop slow fluid from near the wall and dump it in the middle, scoop fast fluid from the middle and dump it near the wall. After millions of these eddy exchanges, the bulk velocity is nearly uniform. Only in the very last micrometres next to the wall does viscosity become strong enough to insist on the no-slip condition — that is the boundary sublayer.

Energy Cascade and the Kolmogorov Microscale

Turbulence is not random noise. It is a beautifully organised machine for moving kinetic energy from the largest scales (where the mean flow injects it) to the smallest scales (where viscosity eventually dissipates it as heat). This is the Kolmogorov energy cascade.

large eddies (size L) → medium eddies → small eddies (size η) → heat

Energy flows from left to right at a rate

\varepsilon

per unit mass. Each step is inviscid until the smallest scale, where viscosity finally wins.

The Kolmogorov microscale

At what scale does viscosity catch up? The smallest eddy must have $Re \sim 1$ : its inertia is barely strong enough to overcome the local viscous damping. Setting the Reynolds number of an eddy of size $\eta$ equal to 1 gives the Kolmogorov length:

\displaystyle \eta = \left(\frac{\nu^3}{\varepsilon}\right)^{\!1/4}, \qquad \nu = \mu/\rho

For water in our 5-cm pipe at $U = 0.2$ m/s, $\nu = 10^{-6}$ m²/s and $\varepsilon \sim U^3/L \approx 1.6\times 10^{-1}$ W/kg. Plugging in:

\displaystyle \eta \approx \left(\frac{(10^{-6})^3}{0.16}\right)^{\!1/4} \approx 5\times 10^{-5}\ \text{m} \;=\; 50\ \mu\text{m}

The cascade in one breath

Energy enters turbulence at the geometric scale $L$ . It exits as heat at the Kolmogorov scale $\eta$ . The ratio $L/\eta \sim Re^{3/4}$ tells you how many orders of magnitude of eddies live between them. For $Re = 10^4$ , that is about three decades — for an aircraft at $Re = 10^9$ , almost seven. That is why direct numerical simulation of high-Re turbulence is impossibly expensive.

Worked Example: Water in a 5-cm Pipe

Try this by hand first. Then expand the solution below to compare.

Problem

Room-temperature water ( $\rho = 1000$ kg/m³, $\mu = 10^{-3}$ Pa·s) flows through a straight pipe of diameter $D = 0.05$ m and length $L = 10$ m. Consider two operating points:

Slow valve: $U_{\text{avg}} = 0.02$ m/s
Fast valve: $U_{\text{avg}} = 0.20$ m/s

For each, find (a) the Reynolds number and the regime, (b) the Darcy friction factor, (c) the pressure drop along the 10 m of pipe, and (d) the pumping power per unit volume of fluid delivered.

Click to reveal the full hand solution

Step 1 — Reynolds numbers

Apply $Re = \rho U D / \mu$ :

Re_{\text{slow}} = \frac{1000 \cdot 0.02 \cdot 0.05}{10^{-3}} = 1000

Re_{\text{fast}} = \frac{1000 \cdot 0.20 \cdot 0.05}{10^{-3}} = 10\,000

One is comfortably laminar (Re = 1 000 < 2 300), the other is comfortably turbulent (Re = 10 000 > 4 000).

Step 2 — Friction factors

Laminar uses $f = 64/Re$ :

f_{\text{slow}} = \frac{64}{1000} = 0.064

Turbulent uses Blasius $f = 0.316/Re^{1/4}$ :

f_{\text{fast}} = \frac{0.316}{10\,000^{1/4}} = \frac{0.316}{10} = 0.0316

Step 3 — Pressure drops over 10 m

Use $\Delta p = f \cdot (L/D) \cdot \tfrac{1}{2} \rho U^2$ :

\Delta p_{\text{slow}} = 0.064 \cdot \tfrac{10}{0.05} \cdot \tfrac{1}{2} \cdot 1000 \cdot (0.02)^2 = 2.56\ \text{Pa}

\Delta p_{\text{fast}} = 0.0316 \cdot \tfrac{10}{0.05} \cdot \tfrac{1}{2} \cdot 1000 \cdot (0.20)^2 = 126.4\ \text{Pa}

Speeding up by 10× cost us 50× more pressure drop. That is the curse of turbulence — the friction factor drops only slowly while the dynamic pressure (½ρU²) grows quadratically.

Step 4 — Pumping power per volume

Volumetric flow rate $Q = (\pi D^2/4)\, U$ :

Q_{\text{slow}} \approx 3.93\times 10^{-5}\ \text{m}^3/\text{s} = 0.039\ \text{L/s}

Q_{\text{fast}} \approx 3.93\times 10^{-4}\ \text{m}^3/\text{s} = 0.393\ \text{L/s}

Pumping power $P = Q \cdot \Delta p$ :

P_{\text{slow}} = 3.93\times 10^{-5} \cdot 2.56 \approx 0.10\ \text{mW}

P_{\text{fast}} = 3.93\times 10^{-4} \cdot 126.4 \approx 49.6\ \text{mW}

Five hundred times more pumping power for ten times more flow. If you ever wondered why oil pipelines run slowly, this is the reason.

Sanity check against the explorer

In the interactive above, set Pure laminar and High-Re presets and compare the friction-factor readout. The numbers (0.080 vs 0.030 here, 0.064 vs 0.0316 in the hand calculation) match the regime predictions to within rounding.

Python: Reynolds, Friction, Profiles

The cleanest way to internalise the regime split is to type the four formulas into Python, run them on the engineering scenario from the worked example, and verify every number by hand. Each line below is annotated — click any line in the code panel to jump to its explanation.

Reynolds number, friction factor, and both pipe profiles from scratch

🐍reynolds_pipe.py

Explanation(43)

Code(99)

1NumPy for vectorised math

We will evaluate u(y) at 200 points and want every arithmetic operator to act elementwise on the whole array. NumPy gives us that for free, which is why the profile formulas remain one-liners.

EXAMPLE

y = np.linspace(-D/2, D/2, 200) — one call, 200 floats.

2Matplotlib for the comparison plot

The end of the script draws the laminar parabola and the turbulent plug on the same axes so the qualitative jump between regimes is visible at a glance.

5reynolds_number — the universal regime detector

One function, three uses: pipe flow (L = diameter), aircraft wings (L = chord), tiny spheres (L = radius). The Reynolds number is dimensionless and tells you, at a glance, whether viscous or inertial physics dominates.

6Inputs: rho, U, L, mu

ρ is the fluid density (kg/m³). U is a representative speed (the mean pipe velocity, the free-stream velocity, ...). L is a representative length (pipe diameter, body chord, ...). μ is the dynamic viscosity (Pa·s). These four numbers carry every piece of physics the equation needs to decide regime.

EXAMPLE

Water at 20°C: ρ=1000 kg/m³, μ=1e-3 Pa·s. Air: ρ=1.2, μ=1.8e-5.

7Docstring: the formula itself

Re = ρ U L / μ. This is the only equation in the entire field of fluid mechanics that gets quoted in every chapter, every textbook, every research paper. Memorise the symbols, not just the value.

15return rho * U * L / mu

One line, four multiplications, one division. The return value has units of (kg/m³)·(m/s)·(m) / (Pa·s) = (kg/(m·s)) / (kg/(m·s)) — dimensionless, exactly as designed.

EXAMPLE

Re for water flowing at U=0.2 m/s in a D=0.05 m pipe = 1000·0.2·0.05/1e-3 = 10000.

18friction_factor — encapsulates two regimes

Engineering pipe-flow design needs one number — the Darcy friction factor f — to predict pressure drop. The catch is that f changes formula across the transition. We bundle both formulas behind one function call.

19Docstring: dp = f · (L/D) · ½ ρ U²

The Darcy-Weisbach pressure-drop equation. The (½ ρ U²) part is the dynamic pressure of the flow; (L/D) is geometry; f is the dimensionless 'how rough is this regime?' coefficient. All three multiply.

22Docstring: laminar branch — f = 64/Re

Derived analytically from the Hagen-Poiseuille solution by integrating the parabolic profile. It is *exact*. The fact that f drops as 1/Re means doubling the speed only modestly raises pressure drop while in the laminar regime.

EXAMPLE

Re=1000 → f=0.064. Re=2000 → f=0.032. Pure 1/Re scaling.

23Docstring: turbulent branch — Blasius f = 0.316/Re^(1/4)

Empirical for smooth pipes in 4000 ≤ Re ≤ 10⁵, fitted to Nikuradse's 1933 experiments. Note the much weaker scaling: f drops only as Re^(-1/4). That is why turbulent pressure drop scales nearly like U² — bad news for pumping costs.

EXAMPLE

Re=10000 → f=0.0316. Re=20000 → f=0.0266. Barely halves over a 2× speed jump.

26if Re < 2300: laminar branch

Re = 2300 is the canonical lower transition threshold for a circular pipe. Below it, *any* tiny disturbance decays exponentially: the laminar solution is asymptotically stable.

27return 64.0 / Re

Closed-form laminar friction factor. The number 64 comes from integrating the Hagen-Poiseuille parabola — see the worked example below for the full derivation.

28if Re >= 4000: turbulent branch

Above Re ≈ 4000 the flow is fully developed turbulent and the Blasius correlation is accurate. Between 2300 and 4000 the flow is intermittent — sometimes laminar, sometimes turbulent — and no single formula works.

29return 0.316 / Re**0.25

Blasius's 1913 fit. Re**0.25 is a slowly varying function — across a factor of 16 in Reynolds number it changes by only a factor of 2. That weak dependence is the calling card of turbulence.

30return NaN for transitional

An honest answer: we do not have a clean closed form here. In engineering practice you either avoid the regime, use the Moody diagram, or fall back to one of the laminar/turbulent formulas as a conservative bound.

33laminar_profile — the parabola

Closed-form Hagen-Poiseuille velocity. y is measured from the pipe centreline (not from the wall), so the formula uses (2y/D)² and is symmetric. We will pass y from −D/2 to D/2 below.

41u_max = 2 · U_avg

Integrating the parabola over the pipe cross-section and dividing by area gives the *average* velocity. The arithmetic comes out to exactly u_max / 2. So if your flow-meter reads 0.02 m/s, the centreline is actually moving at 0.04 m/s.

EXAMPLE

U_avg = 0.02 m/s → u_max = 0.04 m/s.

42u(y) = u_max · (1 − (2y/D)²)

The parabola. Vectorised: y is an array of 200 floats, (2*y/D)**2 broadcasts elementwise, the subtraction broadcasts, the multiplication broadcasts. NumPy returns u as a 200-element array.

45turbulent_profile — the 1/7 power law

Empirical mean velocity profile for fully developed pipe turbulence, named for the exponent 1/7. It captures the *flat plug shape*: the velocity is nearly constant across the core and crashes to zero only in a thin near-wall layer.

53u_cl = 60/49 · U_avg ≈ 1.224 · U_avg

Integrating the 1/7 power law over the pipe cross-section gives U_avg / u_cl = 49/60. So the centreline velocity is only ~22% higher than the average — utterly unlike the laminar case where it is double.

EXAMPLE

U_avg = 0.20 m/s → u_cl ≈ 0.2449 m/s.

54u(y) = u_cl · (1 − 2|y|/D)^(1/7)

np.abs(y) handles the symmetry across the centreline. Raising to the 1/7 power makes the profile rise *very fast* near the wall — that is exactly the thin boundary-layer signature of turbulence.

EXAMPLE

At y = 0.49·D/2: (1 − 0.49)^(1/7) ≈ 0.905 — still 90% of u_cl very close to the wall.

58Engineering scenario: water in a 5-cm pipe

Real numbers: room-temperature water (ρ=1000, μ=1e-3) flowing through a household-scale pipe (D=5 cm). We scan three speeds and watch the regime flip.

59rho = 1000 kg/m³

Density of water at 20 °C. Air is ~800× lighter, mercury is ~13× heavier; for any liquid you handle daily, 1000 is your default.

60mu = 1e-3 Pa·s

Dynamic viscosity of water at 20 °C. Honey is ~10⁴ times this, olive oil is ~100×. Tiny number — but multiplied by an enormous velocity gradient in a boundary layer, it produces real forces.

61D = 0.05 m

Five-centimetre pipe — kitchen-faucet scale, big enough that turbulence is easy to trigger but small enough that laminar flow is reachable at very low speeds.

62L = 10 m

Pipe length over which we will compute total pressure drop. Doubling L doubles dp in either regime — the L/D factor is linear in the Darcy formula.

64speeds = [0.02, 0.05, 0.20]

Three speeds chosen to land squarely in laminar (0.02), transitional (0.05 → Re=2500), and turbulent (0.20 → Re=10000). Watch the friction factor change formulas as the loop steps through these.

66Header for the results table

Right-aligned column widths so the numbers line up nicely. We will print one row per speed.

68for U in speeds: — scan three regimes

Three iterations, each producing one row. The loop body is the entire Re → regime → f → dp computation chain.

69Re = reynolds_number(...)

First iteration: Re = 1000·0.02·0.05/1e-3 = 1000. Second iteration: Re = 2500. Third iteration: Re = 10000. The numbers walk us right across the transition.

70f = friction_factor(Re)

Iteration 1: Re=1000 → f = 64/1000 = 0.064. Iteration 2: Re=2500 (transitional) → f = NaN. Iteration 3: Re=10000 → f = 0.316/10000^0.25 = 0.0316.

71Regime selection — three-branch if/elif/else

We name the regime alongside the number. This is the part of the output that is most informative for engineers — Re=2500 in a NaN row is a red flag that the design lives in the unstable middle.

74dp/L pressure-drop intensity

We compute pressure drop per metre, so the comparison is apples-to-apples regardless of pipe length. The 0.5 · ρ · U² is dynamic pressure, (1/D) comes from the L/D factor with L absorbed into 'per-metre'.

EXAMPLE

Laminar: dp/L = 0.064·(1/0.05)·0.5·1000·0.02² = 0.256 Pa/m. Turbulent: 0.0316·(1/0.05)·0.5·1000·0.20² = 12.64 Pa/m — fifty times more!

75f-string printout

Each iteration prints one tidy row. The visual: as you walk Re from 1000 → 10000, the regime column flips, the friction column drops, but the pressure-drop column rises *much* faster than the velocity — the curse of turbulence.

78Plot the two canonical profiles

Now we make the plot. We pick a laminar case (U_avg = 0.02) and a turbulent case (U_avg = 0.20), evaluate both profiles at 200 points across the pipe, and overlay them.

79U_avg = 0.20

Pick a representative average velocity for the turbulent case. Could be any number — the profile *shapes* are universal, only the magnitudes differ.

80y = np.linspace(-D/2, D/2, 200)

200 sample points across the diameter, including the walls. y = 0 is the centreline, y = ±D/2 are the walls. The arrays u_lam and u_tur will inherit this shape.

81u_lam = laminar_profile(y, D, U_avg=0.02)

Vectorised call: NumPy applies the parabola elementwise. The peak (at y=0) is 2·0.02 = 0.04 m/s. The endpoints (at y = ±D/2) are exactly zero — no-slip again.

82u_tur = turbulent_profile(y, D, U_avg=0.20)

Same shape array, totally different shape *as a function*. The peak (at y=0) is 0.20 · 60/49 ≈ 0.2449 m/s. The flat-plug character will be visible immediately in the plot.

84plt.figure(figsize=(6, 4))

Slightly wider than tall — gives the velocity axis some room and matches the orientation of the pipe (long horizontal, narrow vertical).

85plt.plot(u_lam, y*1000, ...)

Note: u on the horizontal axis, y on the vertical axis. This is the fluid-mechanics convention — the plot visually aligns with the pipe geometry. Multiplying y by 1000 converts metres to millimetres for readability.

86plt.plot(u_tur, y*1000, ...)

Overlay the turbulent profile in red. The visual: a tall, pointy blue parabola squeezed between two flat regions, vs a wide red plug that drops vertically near the walls. Two completely different organisations of the same equations.

89axhline at ±D/2 — the walls

Two grey horizontal lines mark the pipe walls. Both curves *must* touch u=0 there — the no-slip condition is still in force, even in turbulence.

56 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4
5def reynolds_number(rho, U, L, mu):
6    """
7    Dimensionless ratio of inertial forces to viscous forces.
8
9        Re = rho * U * L / mu
10
11    Inputs (SI units):
12        rho : fluid density       [kg / m^3]
13        U   : characteristic speed [m / s]
14        L   : characteristic length [m]   (pipe diameter, chord, ...)
15        mu  : dynamic viscosity   [Pa.s]
16
17    For pipe flow the convention is L = pipe diameter D.
18    """
19    return rho * U * L / mu
20
21
22def friction_factor(Re):
23    """
24    Darcy friction factor f used in the pipe pressure-drop formula
25        dp = f * (L / D) * (1/2) * rho * U^2.
26
27    Laminar  (Re < 2300):       f = 64 / Re                (Hagen-Poiseuille)
28    Turbulent (4000 <= Re):     f = 0.316 / Re^(1/4)       (Blasius, smooth pipes)
29    Transitional region:        no clean closed form -- return NaN.
30    """
31    if Re < 2300:
32        return 64.0 / Re
33    if Re >= 4000:
34        return 0.316 / Re**0.25
35    return float("nan")
36
37
38def laminar_profile(y, D, U_avg):
39    """
40    Parabolic Hagen-Poiseuille profile in a pipe of diameter D.
41
42        u(y) = u_max * (1 - (2y/D)^2),    u_max = 2 * U_avg.
43
44    y is measured from the pipe centreline (y in [-D/2, D/2]).
45    """
46    u_max = 2.0 * U_avg
47    return u_max * (1.0 - (2.0 * y / D) ** 2)
48
49
50def turbulent_profile(y, D, U_avg):
51    """
52    1/7-power-law approximation to the turbulent mean velocity profile.
53
54        u(y) = u_cl * (1 - 2|y|/D)^(1/7),   u_cl = (60/49) * U_avg ~= 1.224 U_avg.
55
56    Reasonable for fully developed pipe turbulence with Re < 1e5.
57    """
58    u_cl = U_avg * 60.0 / 49.0
59    return u_cl * (1.0 - 2.0 * np.abs(y) / D) ** (1.0 / 7.0)
60
61
62# ---- Engineering setup: water flowing in a 5-cm pipe -----------------
63rho = 1000.0        # water density       [kg/m^3]
64mu  = 1.0e-3        # water viscosity     [Pa.s]
65D   = 0.05          # pipe diameter       [m]
66L   = 10.0          # pipe length         [m]
67
68speeds = [0.02, 0.05, 0.20]   # average velocities to scan, [m/s]
69
70print(f"{'U [m/s]':>9} {'Re':>8} {'regime':>14} {'f':>8} {'dp/L [Pa/m]':>14}")
71print("-" * 60)
72for U in speeds:
73    Re = reynolds_number(rho, U, D, mu)
74    f  = friction_factor(Re)
75    if Re < 2300:      regime = "laminar"
76    elif Re < 4000:    regime = "transitional"
77    else:              regime = "turbulent"
78    dpdx = f * (1.0 / D) * 0.5 * rho * U * U if not np.isnan(f) else np.nan
79    print(f"{U:9.3f} {Re:8.0f} {regime:>14} {f:8.4f} {dpdx:14.3f}")
80
81# ---- Plot the two canonical profiles side by side --------------------
82U_avg = 0.20
83y = np.linspace(-D / 2.0, D / 2.0, 200)
84u_lam = laminar_profile(y, D, U_avg=0.02)         # slow flow -> laminar
85u_tur = turbulent_profile(y, D, U_avg=0.20)       # fast flow -> turbulent
86
87plt.figure(figsize=(6, 4))
88plt.plot(u_lam, y * 1000, label="laminar (Re=1000)",   color="#2563eb", lw=2.5)
89plt.plot(u_tur, y * 1000, label="turbulent (Re=10000)", color="#dc2626", lw=2.5)
90plt.axhline(0,            color="k", lw=0.7)
91plt.axhline( D * 500,     color="gray", lw=1.2)
92plt.axhline(-D * 500,     color="gray", lw=1.2)
93plt.xlabel("u(y) [m/s]")
94plt.ylabel("y [mm]   (centreline at 0)")
95plt.title("Laminar parabola vs turbulent plug")
96plt.legend()
97plt.grid(alpha=0.3)
98plt.tight_layout()
99plt.show()

Expected output

Running the script with the three pipe speeds prints:

  U [m/s]       Re         regime        f    dp/L [Pa/m]
------------------------------------------------------------
    0.020     1000        laminar   0.0640          0.256
    0.050     2500   transitional      nan            nan
    0.200    10000      turbulent   0.0316         12.640

The middle row deliberately prints NaN — a polite reminder that the transitional regime resists closed-form treatment. The plot then shows the tall blue parabola squeezed inside the broad red plug, both touching $u = 0$ at the walls.

PyTorch: Reynolds Decomposition With Autograd

Now the modern twist. Real turbulence analysis is built on a single change of variables — Reynolds decomposition: split every velocity into a time-mean and a fluctuation, $u(t) = U + u'(t)$ . The averages of the fluctuations form a tensor (the Reynolds stress) that drives every closure model from k-ε to k-ω to large-eddy simulation.

We can do every step of this analysis with PyTorch tensors, and then — because the whole chain is differentiable — use autograd to compute the sensitivity of the turbulent kinetic energy $k = \tfrac{1}{2}\langle u'^2 + v'^2 \rangle$ with respect to every measured sample. That sensitivity map is how researchers find coherent structures in real data.

Reynolds decomposition, Reynolds stress, TKE, and autograd-based sensitivity

🐍reynolds_stress.py

Explanation(31)

Code(62)

1import torch — autograd is the whole point

We will compute the turbulent kinetic energy k as a function of the velocity samples and then ask PyTorch for ∂k/∂u and ∂k/∂v automatically. No autograd, no story.

3torch.manual_seed(0)

Determinism. The synthetic noise we add to the fake turbulent signal is sampled with a fixed seed so the printed numbers reproduce exactly.

7N = 20000 samples

20 000 samples mimics a 20-second recording from a hot-wire anemometer running at 1 kHz — typical for laboratory turbulence measurements.

8t = linspace(0, 20, N)

Uniform time grid from 0 to 20 seconds. We do not actually use t for arithmetic — only for the sinusoidal drift term below — but it is helpful to picture the signal as a time series.

10U_mean = 0.40, V_mean = 0.0

Pretend we are inside a turbulent pipe flow. The mean streamwise velocity is 0.4 m/s; the mean cross-stream velocity is 0 by symmetry (the flow does not drift sideways on average).

13u_fluct = 0.08·sin(...) + 0.05·randn

Synthetic streamwise fluctuation. The sin term gives a slow, large eddy at 0.3 Hz; the randn term sprinkles in broadband noise representing all the small eddies. The sum is a passable cartoon of real u'(t).

14v_fluct = -0.6 · u_fluct + 0.03·randn

Crucially, v' is *correlated* with u'. The negative sign mimics the near-wall reality: when a high-speed fluid parcel hits the wall (u' > 0), it tends to move toward the wall (v' < 0). This correlation is exactly what makes ⟨u'v'⟩ nonzero — the heart of the Reynolds stress.

EXAMPLE

If u'=+0.1, then on average v' ≈ -0.06. Product u'v' ≈ -0.006 < 0.

16u = U_mean + u_fluct

Reynolds decomposition assembled. u(t) is the total instantaneous velocity — the thing you would actually measure with an anemometer.

17v = V_mean + v_fluct

Same for the cross-stream component. We now have two synthetic time series ready to be analysed.

20u.requires_grad_(True)

Tell autograd to track every operation that touches u from here on. Without this, the gradient of k with respect to u in Step 3 would fail.

21v.requires_grad_(True)

Same for v. The two tensors are now nodes in a computation graph PyTorch is silently building.

24U = u.mean() — the Reynolds operator

Replace 'time average' with 'sample mean'. By construction U should land very close to U_mean = 0.4. The slight wobble is statistical noise from the 20 000 samples.

EXAMPLE

U.item() ≈ 0.4001 — agreement to four decimal places.

25V = v.mean()

Mean cross-stream velocity. Should land very close to 0.

26u_prime = u − U

Define the fluctuating part. By construction ⟨u'⟩ = 0 to machine precision. This is the Reynolds decomposition in one line: every velocity is (mean) + (fluctuation around the mean).

27v_prime = v − V

Same for v. Now we have u' and v' — the actual stars of turbulence theory.

30uu = ⟨u'²⟩

Mean square of the fluctuation, also called the *streamwise variance*. This is half of the streamwise contribution to TKE. Large uu means strong eddies in the flow direction.

31vv = ⟨v'²⟩

Cross-stream variance. In real pipe flow it is roughly 60–70% of uu near the wall — turbulence is anisotropic.

32uv = ⟨u'v'⟩ — the Reynolds stress

The single most important number in engineering turbulence. It quantifies how much streamwise momentum is being transported in the cross-stream direction by the eddies. If the correlation we built in is correct, uv should be negative.

EXAMPLE

Expected: roughly -0.6 · uu = -0.6 · 0.0124 ≈ -0.0074.

35k = ½ (uu + vv)

Turbulent kinetic energy per unit mass. In 3D it would include ww too, but for this 2D illustration uu + vv is the whole budget. k is what RANS closures (k-ε, k-ω) ultimately predict.

37print U_mean

First sanity check. Should agree with 0.40 to several decimals. If it does not, your synthetic data has a bug.

38print uu

Streamwise variance. Roughly 0.08² (from the sin) + 0.05² (from the noise) ≈ 0.0089. Will be very close in the output.

39print vv

Cross-stream variance. From v' = -0.6·u' + 0.03·noise we expect 0.36·uu + 0.0009 ≈ 0.0041.

40print Reynolds <u'v'>

The headline turbulence quantity. Sign matters: negative ⟨u'v'⟩ in this near-wall analogy means momentum is being carried toward the wall, which is the mechanism that flattens the turbulent mean profile in the explorer above.

41print TKE k

Total kinetic energy of the fluctuations, per unit mass. About 0.0065 J/kg in this synthetic case. In a real flow, k(y) drives every closure-based turbulence model.

45torch.autograd.grad(k, [u, v])

The autograd magic. PyTorch returns two tensors of shape (N,) — sens_u and sens_v — each telling you ∂k/∂u_i and ∂k/∂v_i sample-by-sample. This is a per-sample sensitivity map: how much would TKE change if I nudged the i-th measurement by ε?

48torch.topk(sens_u.abs(), 5)

Find the five samples whose fluctuations contribute the most to k. These correspond to the largest peaks of u' in the time series — the strongest eddies in our 20-second window.

49print top-5 sample indices

Sanity check: the top-5 should cluster around the maxima of the sin term, which dominates the eddy energy. This is how researchers identify *coherent structures* in real data.

50print top-5 sensitivities

The actual ∂k/∂u_i values at those samples. Tiny per-sample (≈ 1e-5) because N is large and each sample is one of N contributions to the mean. Multiplied by N, they recover u'_i / N times an N-fold sum.

54dU/dy = 20 /s — typical near-wall shear

Mean shear rate of the underlying flow at this point. In a pipe with U_avg = 0.4 m/s and D = 0.05 m, the wall shear is comfortably in the tens of inverse seconds.

55P = -⟨u'v'⟩ · dU/dy

The *production* term of the turbulent kinetic energy equation. It tells you how fast TKE is being created from the mean flow's kinetic energy by the action of Reynolds stress on the mean shear. This is how turbulence stays alive — feeding on the mean flow's gradient.

EXAMPLE

P > 0 (because uv < 0 and dU/dy > 0) means the flow is energetically *making* turbulence right now.

56print production rate P

If this number is positive and large, your flow is in an energy-injecting state — typical for boundary layers in adverse-pressure-gradient regions. Negative would mean the mean flow is sucking energy out of the turbulence (relaminarisation).

31 lines without explanation

1import torch
2
3torch.manual_seed(0)
4
5# ---- Step 1. Generate a synthetic turbulent velocity signal --------
6# We pretend to be standing at one point in a turbulent flow and
7# recording velocity components (u, v) over time.
8
9N = 20000                     # number of samples (e.g. 20 seconds at 1 kHz)
10t = torch.linspace(0, 20.0, N)
11
12U_mean, V_mean = 0.40, 0.0    # mean velocities in two directions [m/s]
13
14# Synthetic fluctuations: low-frequency drift + broadband noise.
15u_fluct = 0.08 * torch.sin(2 * torch.pi * 0.3 * t) + 0.05 * torch.randn(N)
16v_fluct = -0.6 * u_fluct + 0.03 * torch.randn(N)    # correlated v' with u'
17
18u = U_mean + u_fluct                  # streamwise velocity
19v = V_mean + v_fluct                  # cross-stream velocity
20
21# Mark u and v as parameters of the analysis -- we will differentiate
22# Reynolds-stress and TKE expressions through them.
23u.requires_grad_(True)
24v.requires_grad_(True)
25
26# ---- Step 2. Reynolds decomposition: u = U + u' --------------------
27U = u.mean()
28V = v.mean()
29u_prime = u - U
30v_prime = v - V
31
32# Reynolds stress tensor entries.
33uu = (u_prime ** 2).mean()
34vv = (v_prime ** 2).mean()
35uv = (u_prime * v_prime).mean()
36
37# Turbulent kinetic energy per unit mass.
38k = 0.5 * (uu + vv)
39
40print(f"U_mean        = {U.item():.4f}")
41print(f"<u'^2>        = {uu.item():.4f}")
42print(f"<v'^2>        = {vv.item():.4f}")
43print(f"Reynolds <u'v'> = {uv.item():.4f}")
44print(f"TKE   k         = {k.item():.4f}")
45
46# ---- Step 3. Use autograd to extract sensitivities -----------------
47# d(TKE) / d(u_i) for every sample i. PyTorch returns a tensor of
48# shape (N,) whose entries tell you how much each instantaneous sample
49# contributes to k.
50sens_u, sens_v = torch.autograd.grad(k, [u, v])
51
52# Pick out the five most "energetic" samples for u'.
53top = torch.topk(sens_u.abs(), k=5)
54print("Top-5 sample indices contributing to k:", top.indices.tolist())
55print("Top-5 sensitivities:", [round(v.item(), 5) for v in top.values])
56
57# ---- Step 4. Production term P = -<u'v'> * dU/dy -------------------
58# In a real flow dU/dy is the mean shear; here we just plug in a number
59# to illustrate the calculation.
60dU_dy = torch.tensor(20.0)            # 1/s, a typical near-wall shear
61P = -uv * dU_dy                       # turbulence production rate (Pa/(kg/m^3))
62print(f"\nturbulence production P = {P.item():.4f}")

Why doing this in PyTorch (instead of NumPy)

Two reasons. First, autograd: as soon as we mark $u$ and $v$ as tracked tensors, PyTorch can give us $\partial k / \partial u_i$ for every sample with no extra code. That sensitivity map is the same kind of object an LES-modelling researcher fits with deep nets. Second, GPU scaling: turbulence datasets routinely run to gigabytes, and the same one-liners scale to GPU tensors without modification.

Synthetic vs measured data

The data above is a cartoon — a single sinusoid plus Gaussian noise, with an artificially imposed cross-correlation between $u'$ and $v'$ . Real near-wall turbulence has a much richer spectrum, intermittent bursts, and anisotropy that depends on wall distance. The pipeline is the same; only the numbers shift.

Why This Matters Everywhere

🛢️ Oil & gas pipelines

Operators sit on the edge of the turbulent regime on purpose: low Re is too slow to deliver volume, but the pumping cost scales nearly as $U^3$ once Re crosses 4 000. Every percentage point of friction factor is millions of dollars per year.

✈️ Aircraft drag

The wing of a passenger jet has Re ≈ 10⁷ — turbulent boundary layer everywhere. Roughly half of all fuel burn is paid to overcome turbulent skin friction. Engineers spend careers on laminar-flow control: pushing the transition just a few percent downstream saves billions of dollars in fuel.

❤️ Cardiovascular flow

Blood in the aorta has Re ≈ 4 000 — right in the transition zone. That is why a healthy heart produces a smooth swoosh and a stenosed valve makes a noisy murmur: turbulent fluctuations generate the audible sound a stethoscope picks up.

⛳ Golf-ball dimples

Dimples deliberately trip the boundary layer to turbulence at a lower Reynolds number. Turbulent boundary layers separate later, shrinking the wake and cutting pressure drag — a dimpled ball flies twice as far as a smooth one.

🌪 Atmosphere & weather

The atmospheric boundary layer is the largest turbulent flow most humans interact with. Re for the troposphere is ≳ 10¹². Weather prediction is, at its core, a turbulence-closure problem on a rotating sphere.

🔬 Microfluidics

Lab-on-a-chip channels have Re ≪ 1 — fully Stokes flow. Mixing is a serious problem because eddies do not exist. Engineers add herringbone grooves to force chaotic advection by geometry rather than instability.

Summary

Laminar and turbulent flow are not two different theories — they are two regimes of the same Navier-Stokes equation, separated by a single dimensionless number.

Key Takeaways

Laminar flow is orderly sheets; neighbouring fluid particles stay neighbours; dye stays in a line. Friction factor $f = 64/Re$ is exact.
Turbulent flow is chaotic three-dimensional mixing; dye is shredded in seconds; the mean profile is a flat plug with a thin viscous sublayer. Friction factor $f \approx 0.316/Re^{1/4}$ is empirical.
The Reynolds number $Re = \rho U L / \mu$ is the dimensionless ratio of inertial to viscous forces. It comes out of pure dimensional analysis of the Navier-Stokes equation.
The transition in a pipe occurs near $Re \approx 2300$ (instability begins) and $Re \approx 4000$ (fully turbulent). Between them the flow is intermittent.
Turbulence is a kinetic-energy cascade from the geometric scale $L$ down to the Kolmogorov microscale $\eta = (\nu^3/\varepsilon)^{1/4}$ , where viscosity finally converts motion into heat.
Reynolds decomposition $u = U + u'$ is the foundation of engineering turbulence modelling. The covariance $\langle u'v' \rangle$ is the Reynolds stress, and $k = \tfrac{1}{2}\langle u'^2 + v'^2 \rangle$ is the turbulent kinetic energy.

Laminar vs Turbulent Flow in One Sentence:

“Two regimes of one equation, separated by one number — and that number is named Reynolds.”

Coming Next: In the next section we confront the deepest mystery in fluid mechanics: the Millennium Prize problem. Do smooth solutions to the 3D Navier-Stokes equations always exist for all time? Nobody knows — and there is a one-million-dollar prize for whoever proves it.