Vector calculus in fluid mechanics

Fluid mechanics is where vector calculus is most at home. A flowing fluid is described by a velocity field that varies in space and time, and every one of the operators of Lessons 8 and 9 has a direct physical meaning in it: divergence measures expansion, curl measures spin, the gradient of pressure pushes the fluid, and the Laplacian describes viscous diffusion. This lesson shows how the governing equations of fluid flow are built from these operators. It introduces the material derivative, which follows a moving fluid particle; derives the continuity equation from the divergence theorem; describes two-dimensional flows with a stream function and irrotational flows with a velocity potential; combines simple potential flows to model flow around a cylinder, with and without spin; and derives Bernoulli's equation from Euler's equation of motion. Groundwater flow toward a well, governed by the same mathematics, completes the picture.

Velocity fields and the material derivative

In the Eulerian description used throughout fluid mechanics, the velocity \(\mathbf{v}(x, y, z, t)\) is a vector field: at each point and time it gives the velocity of whatever fluid is passing that point. Streamlines are curves everywhere tangent to the velocity; in steady flow they are the paths that particles follow.

A fluid particle moves through the field, so the rate at which a quantity changes for that particle, the material derivative, has two parts: the change at a fixed point, and the change due to moving to a new point. By the chain rule (Lesson 3):

Applied to the velocity itself, it gives the acceleration of a fluid particle, \(D\mathbf{v}/Dt = \partial\mathbf{v}/\partial t + (\mathbf{v} \cdot \nabla)\mathbf{v}\). The second, convective term means that fluid accelerates even in a perfectly steady flow whenever it moves into a region of different velocity.

Conservation of mass: the continuity equation

Consider any fixed region E with boundary surface S. The mass inside can change only through flow across S, so the rate of increase of mass equals the inward mass flux. Using the divergence theorem (Lesson 8):

Since this holds for every region, the integrands are equal, giving the continuity equation:

For liquids, and for gases at speeds well below the speed of sound, the density is nearly constant and the equation reduces to \(\nabla \cdot \mathbf{v} = 0\): an incompressible flow has zero divergence, so whatever flows into any small region must flow out.

The stream function

For two-dimensional incompressible flow, the condition \(\partial u/\partial x + \partial v/\partial y = 0\) is satisfied automatically if the velocity comes from a stream function \(\psi(x, y)\):

The stream function is constant along streamlines, so its contour map is a picture of the flow. Better still, the difference in ψ between two streamlines equals the volume flow rate (per unit depth) passing between them. Where streamlines crowd together the flow is fast; where they spread apart it is slow.

Vorticity and the velocity potential

The vorticity \(\boldsymbol{\omega} = \nabla \times \mathbf{v}\) is twice the local angular velocity of fluid elements (Lesson 9). Flow outside thin boundary layers and wakes is often nearly irrotational, \(\boldsymbol{\omega} = \mathbf{0}\), and then the velocity is the gradient of a velocity potential, \(\mathbf{v} = \nabla\phi\). If the flow is also incompressible, \(\nabla \cdot \nabla\phi = 0\):

Ideal flow therefore reduces to Laplace's equation, which is linear, so solutions can be added. A few building blocks, written in polar coordinates, cover a remarkable range of flows:

Flow past a cylinder

Adding a uniform flow and a doublet with \(K = Ua^2\) gives

On the circle r = a, ψ = 0, so that circle is a streamline and can be replaced by a solid cylinder: this is ideal flow around a cylinder of radius a. The velocity on the surface is purely tangential, with speed \(2U|\sin\theta|\): zero at the front and rear stagnation points, and twice the free-stream speed at the top and bottom.

Bernoulli's equation (below) then gives the pressure on the surface. With the pressure coefficient \(C_p = (p - p_\infty)/(\tfrac{1}{2}\rho U^2)\), the result is \(C_p = 1 - 4\sin^2\theta\): +1 at the stagnation points and −3 at the top and bottom. The pressure distribution is symmetric front to back, so it produces no drag at all, a famous result called d'Alembert's paradox. Real cylinders do feel drag, because the viscous boundary layer separates from the rear surface and leaves a low-pressure wake; ideal flow is accurate over the front of the cylinder and fails behind it. This is a good reminder that every model has a range of validity.

Adding circulation: the Magnus effect

Adding a free vortex of circulation Γ (for example by spinning the cylinder) breaks the top-bottom symmetry. With clockwise circulation and flow from left to right, the flow speeds up over the top and slows underneath, the stagnation points move down to where \(\sin\theta = -\Gamma/(4\pi Ua)\), and the lower pressure on top produces an upward force perpendicular to the flow. The Kutta-Joukowski theorem gives this lift per unit length:

The same theorem, applied to an aerofoil whose shape creates circulation, is the foundation of wing theory. Spinning balls curve for the same reason, and rotor ships have used spinning cylinders as sails.

import numpy as np

rho, U, a = 1.2, 10.0, 0.1                     # air, 10 m/s, cylinder radius 0.1 m
theta = np.linspace(0, 2 * np.pi, 20001)

for Gamma in [0.0, 2.0, 6.0]:
    # tangential surface speed with clockwise circulation Gamma
    vt = -2 * U * np.sin(theta) - Gamma / (2 * np.pi * a)
    p = 0.5 * rho * (U**2 - vt**2)             # Bernoulli, relative to free stream
    # force per unit length = -integral of p n ds, with n = (cos, sin) and ds = a dtheta
    drag = 0.0 - np.trapezoid(p * np.cos(theta) * a, theta)
    lift = 0.0 - np.trapezoid(p * np.sin(theta) * a, theta)
    s = Gamma / (4 * np.pi * U * a)
    print(f"Gamma = {Gamma:.1f} m2/s: drag {round(drag, 3) + 0.0:.3f} N/m, lift {round(lift, 2) + 0.0:.2f} N/m,"
          f" rho U Gamma = {rho * U * Gamma:.2f} N/m, stagnation at sin(theta) = {round(-s, 3) + 0.0:.3f}")

Gamma = 0.0 m2/s: drag 0.000 N/m, lift 0.00 N/m, rho U Gamma = 0.00 N/m, stagnation at sin(theta) = 0.000
Gamma = 2.0 m2/s: drag 0.000 N/m, lift 24.00 N/m, rho U Gamma = 24.00 N/m, stagnation at sin(theta) = -0.159
Gamma = 6.0 m2/s: drag 0.000 N/m, lift 72.00 N/m, rho U Gamma = 72.00 N/m, stagnation at sin(theta) = -0.477

The numerical integration confirms the theory: zero drag for every circulation, and lift equal to ρUΓ. Once Γ exceeds 4πUa, the stagnation points merge below the cylinder and leave its surface, and a region of fluid circulates with the cylinder.

Euler's equation and Bernoulli's equation

Newton's second law for a fluid particle, with pressure and gravity forces and no viscosity, is Euler's equation:

The pressure term is a force per unit volume pushing fluid from high to low pressure. Using the identity \((\mathbf{v} \cdot \nabla)\mathbf{v} = \nabla(\tfrac{1}{2}|\mathbf{v}|^2) - \mathbf{v} \times \boldsymbol{\omega}\) and writing \(\mathbf{g} = -\nabla(gz)\), steady incompressible flow gives

The right side is perpendicular to the velocity, so along any streamline the bracketed quantity is constant: this is Bernoulli's equation. In an irrotational flow the right side is zero everywhere, and the same constant applies throughout the flow. Adding viscosity contributes a term \(\mu\nabla^2\mathbf{v}\) to Euler's equation, giving the Navier-Stokes equations, which computational fluid dynamics programs solve numerically.

Groundwater flow toward a well

Water seeping through soil obeys Darcy's law, \(\mathbf{v} = -K\nabla h\), where h is the hydraulic head and K the hydraulic conductivity. Combined with incompressibility, \(\nabla \cdot \mathbf{v} = 0\), this gives \(\nabla^2h = 0\): the same Laplace equation as ideal flow. For a well pumping Q from a confined aquifer of thickness b, the flow is radial, and the cylindrical Laplacian (Lesson 9) gives \(h = A\ln r + B\). Fitting the flow rate \(Q = 2\pi rb\,K\,dh/dr\) gives the Thiem equation for the head difference between two observation wells: