Most developers treat Computational Fluid Dynamics (CFD) as a black box. You hand a mesh to commercial software and wait for the colored plots to come back. The usual advice is “just use a library,” which is how people end up with slow, unstable simulations they cannot debug. If you want to understand how fluid motion is actually computed, build a Navier-Stokes Solver in Python yourself.
I have spent 14 years wrestling with complex systems, and one lesson keeps repeating: when a simulation breaks, restarting the software is not a fix, it is giving up. Translating the partial differential equations into discretized code is what lets you see the physics instead of guessing at it. So that is what the rest of this post does, in NumPy.
The physics: momentum and continuity
The Navier-Stokes equations are Newton’s second law applied to fluids. For incompressible flow, two things matter: momentum, which balances pressure against viscosity, and continuity, which keeps mass from appearing out of nowhere. Pressure and velocity are coupled, and that coupling is the hardest part of any Navier-Stokes Solver in Python.
The derivation is written out on Physics StackExchange, but the part that bites you in code is the Pressure-Poisson equation. It has to be solved at every timestep, or the velocity field stops being divergence free.
Discretization on a grid
On a computer this turns into a finite difference scheme over a uniform grid. Each term gets its own treatment:
- Time: forward difference, explicit Euler.
- Advection: backward or upwind difference, which is what keeps it stable.
- Diffusion: central difference.
As I wrote in the post on Fast Python performance, raw Python loops are the bottleneck here. Everything below is vectorized with NumPy, otherwise the run crawls.
Implementing the solver logic
The loop always runs in the same order: build the source term, solve for pressure with Jacobi iteration, then update the velocity field. Vectorized, the source term calculation comes out like this.
# bbioon_calc_source_term
# rho: density, dt: time step, dx/dy: grid spacing
b[1:-1, 1:-1] = (rho * (
1 / dt * ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) +
(v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
(v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx)) -
((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2
))With the source term in hand, the pressure gets iterated. People cut the Jacobi iteration count to save runtime, which is a mistake. Too few iterations and the pressure field never balances the velocity, so the fluid blows up numerically.
# bbioon_pressure_poisson
for _ in range(nit):
pn = p.copy()
p[1:-1, 1:-1] = (
(pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +
(pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2 -
b[1:-1, 1:-1] * dx**2 * dy**2
) / (2 * (dx**2 + dy**2))
# Gauge pressure boundary conditions
p[:, -1] = 0; p[:, 0] = 0; p[-1, :] = 0; p[0, :] = 0Simulating airflow around a wing
A wing is really just a mask. You mark a set of grid points as solid and force the velocity there to zero, which is the no-slip condition. Run this Navier-Stokes Solver in Python and the pressure field does what the textbook says it should: high underneath the wing, low above it, Bernoulli’s principle showing up on its own.
The solver is laminar, though, so it has no turbulence model. Push the Reynolds number too high and it starts to oscillate, then crashes. Worth knowing before you trust anything it prints.
If this kind of work is eating your dev hours, I can take it on. I have been working with WordPress and high-performance computation since the 4.x days.
Pragmatic takeaway
Writing a solver by hand pays off well past the coursework. You learn to handle ordering problems in your data and to manage large arrays with NumPy discretization. Mostly, the commercial tools stop looking like magic, and a failing simulation becomes something you can debug instead of restart.
