Navier-Stokes & Numerical Intro

The Navier-Stokes equations are the master equations of fluid motion. They are nothing more than Newton's second law ($F = ma$) and conservation of mass applied to a tiny parcel of fluid — yet they describe everything from blood flow to hurricanes to the airflow over a wing. They are also famously hard: a proof that smooth 3D solutions always exist is one of the Clay Institute's million-dollar Millennium Prize problems.

In this lesson we'll dissect the equations term by term, look at the handful of cases that can be solved by hand, see why everything else needs CFD, and write our own finite-difference solver in Python.

The Continuity Equation

Conservation of mass for an incompressible fluid says velocity has no net divergence — what flows in must flow out:

$$\nabla \cdot \vec{V} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0$$

The Momentum (Navier-Stokes) Equation

For an incompressible, constant-viscosity Newtonian fluid:

$$\underbrace{\rho \frac{\partial \vec{V}}{\partial t}}_{\text{unsteady}} + \underbrace{\rho (\vec{V} \cdot \nabla)\vec{V}}_{\text{convection}} = \underbrace{-\nabla p}_{\text{pressure}} + \underbrace{\mu \nabla^2 \vec{V}}_{\text{diffusion (viscous)}} + \underbrace{\rho \vec{g}}_{\text{body force}}$$

Each term carries a clear physical meaning:

The troublemaker is the convective term: it is nonlinear, which is what makes turbulence, vortices, and chaos possible — and what makes analytical solutions so rare.

Non-Dimensionalization

Scaling the equations by a characteristic velocity $U$ and length $L$ collapses all the parameters into a single dimensionless group, the Reynolds number:

$$Re = \frac{\rho U L}{\mu} = \frac{\text{inertial forces}}{\text{viscous forces}}$$

The non-dimensional momentum equation becomes:

$$\frac{\partial \vec{V}^}{\partial t^} + (\vec{V}^ \cdot \nabla^)\vec{V}^ = -\nabla^ p^ + \frac{1}{Re}\nabla^{2}\vec{V}^*$$

At low $Re$, viscous diffusion dominates (smooth, laminar). At high $Re$, convection dominates (turbulent). Two flows with the same $Re$ and geometry are dynamically similar — the basis of wind-tunnel testing.

Analytical Solutions: Couette and Poiseuille

When the nonlinear convective term vanishes (fully developed flow between parallel plates), Navier-Stokes reduces to a simple ODE with exact solutions.

$$u(y) = U \frac{y}{h}$$

$$u(y) = \frac{1}{2\mu}\frac{dp}{dx}\left(y^2 - h y\right)$$

import numpy as np
import matplotlib.pyplot as plt

h = 0.01            # gap between plates (m)
mu = 1.0e-3         # dynamic viscosity (Pa.s, water)
U = 0.5             # top plate speed (m/s)
dpdx = -200.0       # pressure gradient (Pa/m)

y = np.linspace(0, h, 100)

u_couette    = U * y / h
u_poiseuille = (1/(2*mu)) * dpdx * (y**2 - h*y)
u_combined   = u_couette + u_poiseuille

plt.figure(figsize=(8, 6))
plt.plot(u_couette,    y*1000, label='Couette (shear-driven)', lw=2)
plt.plot(u_poiseuille, y*1000, label='Poiseuille (pressure-driven)', lw=2)
plt.plot(u_combined,   y*1000, '--', label='Combined', lw=2)
plt.xlabel('Velocity u (m/s)'); plt.ylabel('y (mm)')
plt.title('Analytical Flow Between Parallel Plates')
plt.legend(); plt.grid(alpha=0.3)
plt.show()

u_max = u_poiseuille.max()
print(f"Max Poiseuille velocity = {u_max:.4f} m/s at channel centre")

Max Poiseuille velocity = 2.5000 m/s at channel centre

Why We Need CFD

Exact solutions exist only for a few idealized geometries. Add a bend, a wing, a complex boundary, or turbulence and the nonlinear convective term makes the equations intractable by hand. Computational Fluid Dynamics (CFD) solves Navier-Stokes numerically by:

• Discretizing the domain into a grid (mesh) of points or cells.
• Approximating derivatives with algebraic expressions (finite differences, finite volumes, or finite elements).
• Marching in time or iterating until the solution converges.

Finite-Difference Discretization

A finite difference replaces a derivative with the values at neighboring grid points. The three most common forms for a first derivative are:

$$\left.\frac{\partial u}{\partial x}\right|_i \approx

\begin{cases}

\dfrac{u_{i+1} - u_i}{\Delta x} & \text{forward} \\[2ex]

\dfrac{u_i - u_{i-1}}{\Delta x} & \text{backward} \\[2ex]

\dfrac{u_{i+1} - u_{i-1}}{2\Delta x} & \text{central (2nd order)}

\end{cases}$$

The second derivative uses a central stencil:

$$\left.\frac{\partial^2 u}{\partial x^2}\right|_i \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2}$$

Numerical Example 1: 1D Unsteady Diffusion

The 1D diffusion (heat) equation models how a velocity or temperature disturbance spreads by viscosity:

$$\frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial x^2}$$

Using a forward difference in time and central difference in space (FTCS), the explicit update is:

$$u_i^{n+1} = u_i^n + \frac{\nu \Delta t}{\Delta x^2}\left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right)$$

Stability requires the diffusion number $\nu \Delta t / \Delta x^2 \le 0.5$.

import numpy as np
import matplotlib.pyplot as plt

# Domain and physical parameters
L  = 2.0            # length (m)
nx = 81             # grid points
nu = 0.1            # diffusion coefficient (m^2/s)
dx = L / (nx - 1)

# Stability: choose dt so that nu*dt/dx^2 <= 0.5
sigma = 0.4
dt = sigma * dx**2 / nu
nt = 400

x = np.linspace(0, L, nx)

# Initial condition: a square "hat" pulse
u = np.ones(nx)
u[(x >= 0.5) & (x <= 1.0)] = 2.0

plt.figure(figsize=(9, 6))
snapshots = [0, 40, 120, 250, nt]
u_run = u.copy()
for n in range(nt + 1):
    if n in snapshots:
        plt.plot(x, u_run, label=f't = {n*dt:.3f} s')
    un = u_run.copy()
    u_run[1:-1] = un[1:-1] + sigma * (un[2:] - 2*un[1:-1] + un[:-2])
    # Dirichlet boundaries held at 1.0

plt.xlabel('x (m)'); plt.ylabel('u')
plt.title('1D Unsteady Diffusion (FTCS) — Pulse Spreading Out')
plt.legend(); plt.grid(alpha=0.3)
plt.show()

print(f"dx = {dx:.4f} m, dt = {dt:.5f} s, diffusion number = {sigma}")
print(f"Peak value: start = 2.000, end = {u_run.max():.4f}")

dx = 0.0250 m, dt = 0.00250 s, diffusion number = 0.4
Peak value: start = 2.000, end = 1.4399

The sharp pulse smooths and spreads — exactly how viscosity diffuses momentum. The peak drops as the disturbance is smeared across the domain.

Numerical Example 2: 1D Burgers' Equation

Burgers' equation adds the nonlinear convective term to diffusion, making it a 1D toy model of the full Navier-Stokes momentum balance:

$$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial x^2}$$

We discretize convection with a backward (upwind) difference and diffusion with a central difference.

import numpy as np
import matplotlib.pyplot as plt

nx = 201
L = 2.0
dx = L / (nx - 1)
nu = 0.02
dt = 0.4 * dx**2 / nu       # diffusion-limited stable step
nt = 600

x = np.linspace(0, L, nx)

# Initial condition: smooth sine bump that steepens into a shock-like front
u = 1.0 + 0.5 * np.sin(np.pi * x)

plt.figure(figsize=(9, 6))
snapshots = {0, 150, 300, 450, nt}
for n in range(nt + 1):
    if n in snapshots:
        plt.plot(x, u.copy(), label=f't = {n*dt:.3f} s')
    un = u.copy()
    # Upwind convection + central diffusion (interior points)
    u[1:-1] = (un[1:-1]
               - un[1:-1] * dt/dx * (un[1:-1] - un[:-2])
               + nu * dt/dx**2 * (un[2:] - 2*un[1:-1] + un[:-2]))
    # Periodic-like fixed ends
    u[0] = un[0]
    u[-1] = un[-1]

plt.xlabel('x (m)'); plt.ylabel('u')
plt.title("1D Burgers' Equation — Nonlinear Steepening vs Diffusion")
plt.legend(); plt.grid(alpha=0.3)
plt.show()

print(f"dt = {dt:.5f} s,  CFL (u*dt/dx) ~ {u.max()*dt/dx:.3f}")
print(f"Final max gradient = {np.max(np.abs(np.gradient(u, dx))):.3f}")

dt = 0.00200 s,  CFL (u*dt/dx) ~ 0.200
Final max gradient = 5.931

Watch the wave steepen: faster fluid catches up to slower fluid (the nonlinear convective term), while viscosity ($\nu$) fights to keep the front smooth. This tug-of-war is the essence of real fluid dynamics in one dimension.

Common Pitfalls

• Violating the stability limit. For explicit diffusion, $\nu\Delta t/\Delta x^2 \le 0.5$; for convection the CFL number $u\Delta t/\Delta x \le 1$. Exceed these and the solution blows up.
• Central differencing convection on coarse grids. It produces unphysical oscillations; upwind differencing is more robust (at the cost of numerical diffusion).
• Forgetting boundary conditions. Every numerical scheme needs explicit treatment of the domain edges.
• Confusing dynamic and kinematic viscosity. $\nu = \mu/\rho$ appears in the non-dimensional and diffusion forms.
• Expecting the discrete answer to be exact. Finite differences carry truncation error that shrinks as $\Delta x, \Delta t \to 0$.

Key Takeaways

• The Navier-Stokes equations are $F=ma$ plus mass conservation for a fluid, with unsteady, convective, pressure, viscous, and body-force terms.
• The convective term is nonlinear, which is why turbulence exists and why most flows have no analytical solution.
• Non-dimensionalization reduces the parameters to the Reynolds number, governing the balance of inertia and viscosity.
• Couette and Poiseuille flows are exact solutions; everything more complex needs CFD.
• Finite differences turn derivatives into algebra on a grid; explicit schemes (FTCS, upwind) are simple but must respect stability limits.

Next, we'll bring the whole course together in a capstone project: a complete Python pipe-network solver.