Five solvers run on this page. Each one is small enough to read and large enough to show the physics: an incompressible flow with a dye source, a cylinder wake shedding vortices at a frequency the page measures, a boundary layer solved two independent ways that agree, a potential flow whose lift is checked by a momentum integral, and the energy spectrum of the first solver's field. The equations beside each figure are the ones it runs.
For an incompressible fluid with kinematic viscosity \nu,
\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}, \qquad \nabla \cdot \mathbf{u} = 0The four terms on the left and right are local acceleration, advection, pressure gradient, viscous diffusion and body force. The continuity condition is not an evolution equation; it is a constraint the pressure enforces, which is what the projection step does.
Stam's scheme (1999) advances a provisional velocity \mathbf{u}^* with forcing, implicit diffusion and semi-Lagrangian advection, then removes its divergence with a pressure solve:
\nabla^2 p = \nabla \cdot \mathbf{u}^*, \qquad \mathbf{u} = \mathbf{u}^* - \nabla pThis is the Helmholtz–Hodge decomposition: any field is a divergence-free part plus a gradient, and subtracting the gradient leaves the part that satisfies continuity. Implicit diffusion and back-traced advection make the scheme stable for any time step, at the cost of numerical dissipation, which is why the solver's Reynolds number is nominal and the divergence residual is shown.
Re = \frac{\rho U L}{\mu} = \frac{U L}{\nu}The ratio of inertial to viscous stresses, and the one parameter that fixes the flow regime once geometry is set. For pipe flow the classical transition sits near Re \approx 2300, with fully turbulent flow above about 4000; for the cylinder in the next figure the wake goes unsteady near Re \approx 47. Both figures take Re from a slider and set \nu from it.
\frac{D\boldsymbol{\omega}}{Dt} = (\boldsymbol{\omega} \cdot \nabla)\mathbf{u} + \nu \nabla^2 \boldsymbol{\omega}, \qquad \boldsymbol{\omega} = \nabla \times \mathbf{u}The curl of the momentum equation. The first term on the right is vortex stretching and tilting; in two dimensions the vorticity is normal to the plane and that term vanishes, so vorticity is only advected and diffused. That is why the vorticity view of the solver is the clearest one, and why two-dimensional turbulence cascades differently from three-dimensional, as the spectrum figure below shows.
\Gamma = \oint_C \mathbf{u} \cdot d\mathbf{l}The line integral of velocity around a closed contour; by Stokes' theorem, the vorticity flux through it.
\omega_z = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}The colour in the figure above: blue for clockwise, mint for counter-clockwise, shed alternately from the two sides of the cylinder.
St = \frac{f D}{U}Measured live from the sign changes of the transverse velocity at the probe. For a circular cylinder in the laminar shedding range it sits near 0.2.
Near a flat plate the boundary-layer equations admit a similarity solution. With \eta = y\sqrt{U/\nu x} and stream function \psi = \sqrt{\nu U x}\, f(\eta), so that u/U = f'(\eta),
f''' + \tfrac{1}{2} f f'' = 0, \qquad f(0) = f'(0) = 0, \qquad f'(\eta \to \infty) = 1A two-point boundary value problem: the figure shoots on f''(0) until the far condition holds. From the solution,
\delta_{99} = \eta_{99}\sqrt{\frac{\nu x}{U}}, \qquad C_f = \frac{2 f''(0)}{\sqrt{Re_x}}, \qquad \delta^* = \int_0^\infty \Big(1 - \frac{u}{U}\Big)\,dy, \qquad \theta = \int_0^\infty \frac{u}{U}\Big(1 - \frac{u}{U}\Big)\,dywith H = \delta^*/\theta the shape factor. Every one of these appears in the readouts as a quadrature of the computed f, not as a typed constant; the familiar values, \eta_{99} \approx 4.91, C_f\sqrt{Re_x} \approx 0.664, H \approx 2.59, are what the quadratures return.
The boundary-layer equations themselves, without the similarity assumption, are
u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \nu\frac{\partial^2 u}{\partial y^2}, \qquad \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0Parabolic in x, so they can be marched downstream. At each station the figure solves the y-direction implicitly, a tridiagonal system, then integrates continuity for v. Started from a uniform profile the march forgets its start within a few percent of the plate and lands on the Blasius curve; the maximum deviation at the trailing edge is a readout.
Past transition the layer is turbulent, and the empirical power-law fits replace the similarity solution:
\frac{u}{U} \approx \Big(\frac{y}{\delta}\Big)^{1/7}, \qquad \delta \approx 0.37\, x\, Re_x^{-1/5}, \qquad C_f \approx 0.059\, Re_x^{-1/5}, \qquad H \approx 1.3A turbulent layer grows as x^{4/5} rather than x^{1/2} and carries several times the wall shear.
Irrotational, incompressible flow has a velocity potential with \nabla^2\phi = 0, and Laplace's equation is linear: solutions add. In the complex plane the elements on the page are
W(z) = U e^{-i\alpha} z + \sum_k \frac{m_k}{2\pi}\ln(z - z_k) + \frac{i\Gamma}{2\pi}\ln z + \frac{U a^2 e^{i\alpha}}{z}, \qquad u - iv = \frac{dW}{dz}uniform stream, sources and sinks, a clockwise vortex and a doublet. Stream plus doublet is a cylinder of radius a; adding the vortex moves its stagnation points to \sin\theta_s = -\Gamma/4\pi U a, and they leave the surface when |\Gamma| > 4\pi U a.
The lift per unit span on any closed body in a uniform stream with circulation \Gamma around it is
L' = \rho U \Gammaand the drag is zero. The figure does not take that on trust: it integrates pressure from Bernoulli, p = p_\infty + \tfrac{1}{2}\rho(U^2 - |\mathbf{V}|^2), together with the momentum flux around a circle enclosing the body,
\mathbf{F} = -\oint_C \big( p\,\hat{\mathbf{n}} + \rho\,\mathbf{V}(\mathbf{V} \cdot \hat{\mathbf{n}}) \big)\, dsand reports both the lift and the drag that come out. The lift matches \rho U \Gamma to the quadrature's precision and the drag is zero to the same precision, which is d'Alembert's paradox measured rather than stated. The pressure coefficient in the colour field is C_p = 1 - |\mathbf{V}|^2/U^2, so the cylinder's minimum of −3 at the shoulders with no circulation is a check on the field.
In three dimensions, energy injected at a large scale is passed to smaller eddies until viscosity dissipates it at the Kolmogorov scale. In the inertial range between, Kolmogorov's 1941 argument gives the spectrum from the dissipation rate \varepsilon alone:
E(k) = C_K\, \varepsilon^{2/3} k^{-5/3}, \qquad \eta = \Big(\frac{\nu^3}{\varepsilon}\Big)^{1/4}, \qquad \tau_\eta = \Big(\frac{\nu}{\varepsilon}\Big)^{1/2}, \qquad u_\eta = (\nu\varepsilon)^{1/4}In two dimensions vortex stretching is absent and enstrophy is conserved as well as energy, so there are two cascades (Kraichnan 1967): energy goes to larger scales as k^{-5/3} and enstrophy to smaller scales as k^{-3}. The solver above is two-dimensional and forced at the largest scale, so the mid range should sit near −3, steeper still where numerical dissipation bites. The figure's fitted slope is the measurement.
Engineering solvers rarely resolve the cascade. Splitting u_i = \bar{u}_i + u_i' and averaging gives the Reynolds-averaged 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 \partial x_j} - \frac{\partial}{\partial x_j}\overline{u_i' u_j'}with the Reynolds stress \overline{u_i' u_j'} left to a model. The two-equation closures carry the turbulent kinetic energy k and a dissipation variable:
\frac{Dk}{Dt} = P_k - \varepsilon + \frac{\partial}{\partial x_j}\Big[\Big(\nu + \frac{\nu_t}{\sigma_k}\Big)\frac{\partial k}{\partial x_j}\Big], \qquad \frac{D\varepsilon}{Dt} = C_{\varepsilon 1}\frac{\varepsilon}{k}P_k - C_{\varepsilon 2}\frac{\varepsilon^2}{k} + \frac{\partial}{\partial x_j}\Big[\Big(\nu + \frac{\nu_t}{\sigma_\varepsilon}\Big)\frac{\partial \varepsilon}{\partial x_j}\Big]with eddy viscosity \nu_t = C_\mu k^2/\varepsilon in the k–\varepsilon model, or in the k–\omega form better behaved at walls,
\frac{Dk}{Dt} = P_k - \beta^* k\omega + \frac{\partial}{\partial x_j}\Big[(\nu + \sigma_k \nu_t)\frac{\partial k}{\partial x_j}\Big], \qquad \frac{D\omega}{Dt} = \alpha\frac{\omega}{k}P_k - \beta\omega^2 + \frac{\partial}{\partial x_j}\Big[(\nu + \sigma_\omega \nu_t)\frac{\partial \omega}{\partial x_j}\Big], \qquad \nu_t = \frac{k}{\omega}None of these run on the page; they are here because they are what the direct solvers above are usually replaced by at engineering Reynolds numbers.
| Number | Definition | Ratio | On this page |
|---|---|---|---|
| Reynolds | Re = \dfrac{\rho U L}{\mu} = \dfrac{U L}{\nu} | inertia to viscous stress | the slider on the solver and the cylinder wake |
| Strouhal | St = \dfrac{f L}{U} | shedding frequency to flow rate | measured at the probe in the cylinder wake |
| Mach | Ma = \dfrac{U}{c} | flow speed to sound speed | below 0.3 the incompressible equations on this page hold |
| Froude | Fr = \dfrac{U}{\sqrt{g L}} | inertia to gravity | free-surface flows, not modelled here |
| Weber | We = \dfrac{\rho U^2 L}{\sigma} | inertia to surface tension | droplet break-up, not modelled here |
| Prandtl | Pr = \dfrac{\nu}{\alpha} = \dfrac{c_p \mu}{k} | momentum to heat diffusivity | coupled heat transfer, not modelled here |
1// One projection step, as it runs on this page (fluidCore.ts).2// u, v are cell-centred velocities; the channel is periodic in x.3privateproject(u, v, p, div, iters){4for(let j =1; j < ny -1; j++){5const row = j * nx;6for(let i =0; i < nx; i++){7const l = row +(i ===0? nx -1: i -1);8const r = row +(i === nx -1?0: i +1);9// divergence of the provisional field, h = 110 div[row + i]=-0.5*(u[r]- u[l]+ v[row + i + nx]- v[row + i - nx]);11 p[row + i]=0;12}13}14this.bounds(0, div);this.bounds(0, p);15// Gauss–Seidel on ∇²p = ∇·u*16this.linSolve(0, p, div,1,4, iters);17// subtract the gradient: u ← u* − ∇p18for(let j =1; j < ny -1; j++){19const row = j * nx;20for(let i =0; i < nx; i++){21const l = row +(i ===0? nx -1: i -1);22const r = row +(i === nx -1?0: i +1);23 u[row + i]-=0.5*(p[r]- p[l]);24 v[row + i]-=0.5*(p[row + i + nx]- p[row + i - nx]);25}26}27this.bounds(1, u);this.bounds(1, v);28}
The full solver is under three hundred lines: forcing, two implicit diffusion solves, two semi-Lagrangian advections, this projection, then the dye. The lattice-Boltzmann wake is a separate three-hundred-line D2Q9 kernel. Both run at the frame rate on a laptop, which is the point of building them in the page rather than showing pictures of them.
Semi-Lagrangian advection with linear interpolation smears small scales; the effective Reynolds number is well below the nominal one and the spectrum falls off steeper than −3. Higher-order interpolation, vorticity confinement or a MacCormack correction recover some of it; a spectral or finite-volume scheme recovers more at a cost the page cannot pay at frame rate.
Grid cells set the smallest eddy the solvers can hold. At 256 × 96 and 240 × 96 there is roughly a decade of scales between the forcing and the grid, which is enough to see a cascade and not enough to measure an exponent to two figures.
Both direct solvers are planar. Real wakes and shear layers go three-dimensional above modest Reynolds numbers, and the cascade changes direction. The Kármán street's Strouhal number stays close to its two-dimensional value, which is why the measurement in the wake figure is meaningful; the spectrum's slope does not.
The boundary-layer march stops where the pressure gradient turns adverse and the layer separates; the marching equations become singular there (Goldstein 1948). Predicting separation and the transition to turbulence upstream of it is what most of practical aerodynamics is about, and neither is on this page.