A Tour of Fluid Flow Regimes, Checked Against Exact Solutions#

Fluid dynamics has few exact solutions, and each one describes a single flow regime: slow viscous flow in a channel, a thin boundary layer, inviscid flow around a body, vortex motion, shock waves, and decaying 2D flow. This tutorial visits each regime with physicskit.fluids and compares the numerical output with its exact solution. Those comparisons tell you which solver fits which regime, and how far each one can be trusted.

Note

Units. The functions in physicskit.fluids take and return values in whatever consistent unit system you pass in. SI is used for the dimensional examples below, and the idealized problems (the cylinder, vortices, Sod, Taylor-Green) are non-dimensional, as they are in the literature.

Which regime? The Reynolds number#

The Reynolds number \(Re = UL/\nu\) compares inertia with viscosity, and it largely decides which of the following sections applies. Water (\(\nu = 10^{-6}\) m²/s) flowing at 1 m/s through a 1 cm pipe:

import numpy as np
import physicskit.fluids as fl

print(round(fl.reynolds_number(velocity=1.0, length=0.01, nu=1.0e-6)))
# 10000

That is well above the pipe-flow transition near \(Re \approx 2300\), so this flow would be turbulent. Laminar, viscous-dominated solutions need much smaller \(Re\).

Viscous flow: Poiseuille and Blasius#

Between two plates a distance \(h\) apart, a constant pressure gradient drives the parabolic Poiseuille profile \(u(y) = -\frac{1}{2\mu}\frac{dp}{dx}\,y(h - y)\), an exact solution of the full Navier-Stokes equations. Integrating it gives the flow rate per unit depth, \(Q = -\frac{dp}{dx}\frac{h^3}{12\mu}\):

h, mu, dpdx = 1.0e-3, 1.0e-3, -100.0          # 1 mm gap, water, 100 Pa/m
y = np.linspace(0.0, h, 201)
u = fl.poiseuille_flow_velocity(y, dpdx=dpdx, mu=mu, h=h)
Q = fl.poiseuille_flow_rate(dpdx=dpdx, mu=mu, h=h)
print(f"u_max = {u.max():.4e} m/s, Q = {Q:.4e} m^2/s, "
      f"integral of u = {np.trapezoid(u, y):.4e}, -dpdx h^3/(12 mu) = {-dpdx * h**3 / (12 * mu):.4e}")
# u_max = 1.2500e-02 m/s, Q = 8.3333e-06 m^2/s, integral of u = 8.3331e-06, -dpdx h^3/(12 mu) = 8.3333e-06

At high \(Re\), viscosity only matters in a thin boundary layer next to a wall (Prandtl, 1904). Over a flat plate, Blasius (1908) reduced the layer to the ODE \(2f''' + f f'' = 0\) with \(f(0) = f'(0) = 0\), \(f'(\infty) = 1\). blasius_solve() solves it by shooting. Its wall-shear constant \(f''(0) = 0.33206\) gives the skin-friction law \(C_f = 0.664/\sqrt{Re_x}\) and the thickness \(\delta_{99} \approx 4.91\,x/\sqrt{Re_x}\):

blasius = fl.blasius_solve()
print(round(float(blasius["fpp"][0]), 5), round(float(blasius["fp"][-1]), 6))
# 0.33206 1.0

x = np.array([0.1, 1.0])                        # air at 1 m/s, 10 cm and 1 m from the edge
print(fl.blasius_boundary_layer_thickness(x, U_inf=1.0, nu=1.5e-5).round(5))
# [0.00601 0.01902]

print(fl.blasius_skin_friction_coefficient(np.array([1.0e5])).round(6))
# [0.0021]

The layer is 6 mm thick after 10 cm and grows like \(\sqrt{x}\), to 19 mm after one metre.

Inviscid flow: the cylinder, d’Alembert and Kutta-Joukowski#

Outside the boundary layer the flow is nearly inviscid and irrotational, so it can be built by superposing elementary potential flows (PotentialFlow). A uniform stream plus a doublet gives flow around a cylinder of radius \(a\). On its surface the speed is \(2U\sin\theta\), so Bernoulli gives \(C_p = 1 - 4\sin^2\theta\):

flow = fl.flow_past_cylinder(U_inf=1.0, radius=1.0)
theta = np.array([0.0, np.pi / 6, np.pi / 2, np.pi])
R = 1.0 + 1e-9                                  # just outside the surface
print(flow.pressure_coefficient(R * np.cos(theta), R * np.sin(theta)).round(4))
print((1 - 4 * np.sin(theta) ** 2).round(4))
# [ 1.  0. -3.  1.]
# [ 1.  0. -3.  1.]

The stagnation points (\(C_p = 1\)) and the shoulder (\(C_p = -3\)) match exactly. Integrating the surface pressure gives the force on the cylinder, and adding a bound vortex of circulation \(\Gamma\) breaks the top-bottom symmetry:

theta = np.linspace(0.0, 2 * np.pi, 2001)[:-1]
dtheta = theta[1] - theta[0]
for Gamma in [0.0, -4.0]:
    flow = fl.flow_past_cylinder(U_inf=1.0, radius=1.0, circulation=Gamma)
    cp = flow.pressure_coefficient(R * np.cos(theta), R * np.sin(theta))
    drag = -0.5 * np.sum(cp * np.cos(theta)) * dtheta     # per unit span, rho = U = a = 1
    lift = -0.5 * np.sum(cp * np.sin(theta)) * dtheta
    print(f"Gamma={Gamma:+.1f}: drag = {drag:+.6f}, lift = {lift:+.6f}, "
          f"Kutta-Joukowski = {fl.kutta_joukowski_lift(rho=1.0, U_inf=1.0, circulation=Gamma):+.1f}")

# Gamma=+0.0: drag = +0.000000, lift = -0.000000, Kutta-Joukowski = -0.0
# Gamma=-4.0: drag = +0.000000, lift = +4.000000, Kutta-Joukowski = +4.0

The drag is zero in both cases. This is d’Alembert’s paradox (1752): ideal flow exerts no drag on a body, and real drag comes from the boundary layer and its separation. The lift is \(L = -\rho U\Gamma\) exactly (Kutta, 1902; Joukowski, 1906), which is how circulation produces the lift on a wing.

Vortex dynamics: point vortices#

In 2D, a flow’s vorticity can be concentrated in point vortices, which move only through the velocities they induce on each other (Helmholtz, 1858; Kirchhoff, 1876). PointVortexSystem integrates their motion. Two equal vortices orbit their midpoint at \(\Omega = \Gamma/(\pi d^2)\), and two opposite vortices travel together in a straight line at \(U = \Gamma/(2\pi d)\):

pair = fl.PointVortexSystem(positions=[[-0.5, 0.0], [0.5, 0.0]], circulations=[1.0, 1.0])
t, traj = pair.trajectory(dt=0.01, n_steps=1000)
sep = traj[:, 1] - traj[:, 0]
angle = np.unwrap(np.arctan2(sep[:, 1], sep[:, 0]))
print(f"Omega = {angle[-1] / t[-1]:.5f}, Gamma/(pi d^2) = {1 / np.pi:.5f}")
# Omega = 0.31831, Gamma/(pi d^2) = 0.31831

dipole = fl.PointVortexSystem(positions=[[0.0, -0.5], [0.0, 0.5]], circulations=[1.0, -1.0])
t, traj = dipole.trajectory(dt=0.01, n_steps=1000)
print(f"U = {(traj[-1, 0, 0] - traj[0, 0, 0]) / t[-1]:+.5f}, Gamma/(2 pi d) = {1 / (2 * np.pi):.5f}")
# U = -0.15915, Gamma/(2 pi d) = 0.15915

Both match to all five digits. The dipole is the 2D analogue of a smoke ring, and a staggered double row of such vortices is the von Kármán street (von_karman_vortex_street()).

Compressible flow: shocks#

Once the Mach number exceeds 1, the flow can steepen into shocks, discontinuities across which mass, momentum and energy are conserved (the Rankine-Hugoniot conditions). For a normal shock at \(M_1 = 2\) in air (\(\gamma = 1.4\)) the jumps are standard textbook table entries:

print({k: round(v, 4) for k, v in fl.normal_shock_relations(2.0).items()})
# {'p2_p1': 4.5, 'rho2_rho1': 2.6667, 'T2_T1': 1.6875, 'M2': 0.5774}

The Sod shock tube (Sod, J. Comput. Phys. 27, 1, 1978) tests a time-dependent solver. At \(t = 0\) a diaphragm separates gas at \((\rho, u, p) = (1, 0, 1)\) from gas at \((0.125, 0, 0.1)\), and removing it produces a rarefaction, a contact discontinuity and a shock. The exact Riemann solution has, between the contact and the shock, \(p^* = 0.30313\), \(u^* = 0.92745\) and \(\rho = 0.26557\), with the shock at \(x = 0.8504\) at \(t = 0.2\):

sod = fl.sod_shock_tube(nx=800)
x, rho, u, p = sod["x"], sod["rho"], sod["u"], sod["p"]
plateau = (x > 0.72) & (x < 0.78)
print(f"p* = {p[plateau].mean():.4f}, u* = {u[plateau].mean():.4f}, rho = {rho[plateau].mean():.4f}")
# p* = 0.3031, u* = 0.9276, rho = 0.2701

shock = x[np.argmax(np.abs(np.gradient(p)) * (x > 0.75))]
print(f"shock at x = {shock:.3f}")
# shock at x = 0.852

Pressure and velocity are right to four digits and the shock is in the right place. The density is 2% high because the first-order Lax-Friedrichs scheme smears the contact discontinuity over many cells, and the smeared contact reaches into this window. That is the known weakness of the scheme, and the reason higher-order Godunov-type schemes exist.

Viscous decay: Taylor-Green vortices#

Finally, the full 2D incompressible Navier-Stokes equations. NavierStokes2D is a pseudo-spectral solver in vorticity-streamfunction form on a periodic box. The Taylor-Green vortex \(\omega = 2\sin x\sin y\) (Taylor and Green, 1937) is an exact solution: its nonlinear term vanishes identically, so it keeps its shape and decays as \(e^{-2\nu t}\):

nu = 0.05
solver = fl.NavierStokes2D(n=64, length=2 * np.pi, nu=nu)
omega0 = 2 * np.sin(solver.X) * np.sin(solver.Y)
out = solver.simulate(omega0, dt=0.01, steps=200)
print(f"{np.abs(out['omega']).max():.6f} {2 * np.exp(-2 * nu * 2.0):.6f}")
# 1.637462 1.637462

The peak vorticity after \(t = 2\) matches to six digits, the spectral accuracy expected for a smooth solution. Start the same solver from a random vorticity field instead and it evolves into decaying 2D turbulence; energy_spectrum() measures the resulting kinetic-energy spectrum.

Where to go next#