(CFD) is often seen as a black box of complex commercial software. However, implementing a solver “from scratch” is one of the most powerful ways to learn the physics of fluid motion. I started this as a personal project, and as a part of a course on Biophysics, I took it as an opportunity to finally understand how these beautiful simulations work.
This guide is designed for data scientists and engineers who want to move beyond high-level libraries and understand the underlying mechanics of numerical simulations by translating partial differential equations into discretized Python code. We will also explore fundamental programming concepts like vectorized operations with NumPy and stochastic convergence, which are essential skills for everyone interested in broader scientific computing and machine learning architectures.
We will walk through the derivation and Python implementation of a simple incompressible Navier-Stokes (NS) solver. And then, we will then apply this solver to simulate airflow around a bird’s wing profile.
The Physics: Incompressible Navier-Stokes
The fundamental equations of CFD are the Navier-Stokes equations, which describe how velocity and pressure evolve in a fluid. For steady flight (like a bird gliding), we assume the air is incompressible (constant density) and laminar. They can be understood as just Newton’s motion law, but for an infinitesimal element of fluid, with the forces that affect it. That’s mostly pressure and viscosity, but depending on the context you could add in gravity, mechanical stresses, or even electromagnetism if you’re feeling like it. The author can attest to this being very much not recommend for a first project.
The equations in vector form are:
\[
\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{v} \\
\nabla \cdot \mathbf{v} = 0
\]
Where:
- v: Velocity field (u,v)
- p: Pressure
- ρ: Fluid density
- ν: Kinematic viscosity
The first equation (Momentum) balances inertia against pressure gradients and viscous diffusion. The second equation (Continuity) enforces that the fluid density remains constant.
The Pressure Coupling Problem
A major challenge in CFD is that pressure and velocity are coupled: the pressure field must adjust constantly to ensure the fluid remains incompressible.
To solve this, we derive a Pressure-Poisson equation by taking the divergence of the momentum equation. In a discretized solver, we solve this Poisson equation at every single timestep to update the pressure, ensuring the velocity field remains divergence-free.
Discretization: From Math to Grid
To solve these equations on a computer, we use Finite Difference schemes on a uniform grid.
- Time: Forward difference (Explicit Euler).
- Advection (Nonlinear terms): Backward/Upwind difference (for stability).
- Diffusion & Pressure: Central difference.
For example, the update formula for the u (x-velocity) component looks like this in finite-difference form
\[
u_{i,j}^n \frac{u_{i,j}^n - u_{i-1,j}^n}{\Delta x}
\]
In code, the advection term u∂x∂u uses a backward difference:
\[
u_{i,j}^n \frac{u_{i,j}^n - u_{i-1,j}^n}{\Delta x}
\]
The Python Implementation
The implementation proceeds in four distinct steps using NumPy arrays.
1. Initialization
We define the grid size (nx, ny), time step (dt), and physical parameters (rho, nu). We initialize velocity fields (u,v) and pressure (p) to zeros or a uniform flow.
2. The Wing Geometry (Immersed Boundary)
To simulate a wing on a Cartesian grid, we need to mark which grid points lie inside the solid wing.
- We load a wing mesh (e.g., from an STL file).
- We create a Boolean mask array where
Trueindicates a point inside the wing. - During the simulation, we force velocity to zero at these masked points (no-slip/no-penetration condition).
3. The Main Solver Loop
The core loop repeats until the solution reaches a steady state. The steps are:
- Build the Source Term (b): Calculate the divergence of the velocity terms.
- Solve Pressure: Solve the Poisson equation for p using Jacobi iteration.
- Update Velocity: Use the new pressure to update u and v.
- Apply Boundary Conditions: Enforce inlet velocity and zero velocities inside the wing.
The Code
Here is how the core mathematical updates look in Python (vectorized for performance).
Step A: Building the Pressure Source Term This represents the Right-Hand Side (RHS) of the Poisson equation based on current velocities.
# b is the source term
# u and v are current velocity arrays
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
))Step B: Solving for Pressure (Jacobi Iteration) We iterate to smooth out the pressure field until it balances the source term.
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))
# Boundary conditions: p=0 at edges (gauge pressure)
p[:, -1] = 0; p[:, 0] = 0; p[-1, :] = 0; p[0, :] = 0Step C: Updating Velocity Finally, we update the velocity using the explicit discretized momentum equations.
un = u.copy()
vn = v.copy()
# Update u (x-velocity)
u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
un[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, 0:-2]) +
nu * (dt / dx**2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
dt / dy**2 * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])))
# Update v (y-velocity)
v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
un[1:-1, 1:-1] * dt / dx * (vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
dt / (2 * rho * dy) * (p[2:, 1:-1] - p[0:-2, 1:-1]) +
nu * (dt / dx**2 * (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
dt / dy**2 * (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))Results: Does it Fly?
We ran this solver on a rigid wing profile with a constant far-field inflow.
Qualitative Observations The results align with physical expectations. The simulations show high pressure beneath the wing and low pressure above it, which is exactly the mechanism that generates lift. Velocity vectors show the airflow accelerating over the top surface (Bernoulli’s principle).
Forces: Lift vs. Drag By integrating the pressure field over the wing surface, we can calculate lift.
- The solver demonstrates that pressure forces dominate viscous friction forces by a factor of nearly 1000x in air.
- As the angle of attack increases (from 0∘ to −20∘), the lift-to-drag ratio rises, matching trends seen in wind tunnels and professional CFD packages like OpenFOAM.
Limitations & Next Steps
While making this solver was great for learning, the tool itself has its limitations:
- Resolution: 3D simulations on a Cartesian grid are computationally expensive and require coarse grids, making quantitative results less reliable.
- Turbulence: The solver is laminar; it lacks a turbulence model (like k−ϵ) required for high-speed or complex flows.
- Diffusion: Upwind differencing schemes are stable but numerically diffusive, potentially “smearing” out fine flow details.
Where to go from here? This project serves as a starting point. Future improvements could include implementing higher-order advection schemes (like WENO), adding turbulence modeling, or moving to Finite Volume methods (like OpenFOAM) for better mesh handling around complex geometries. There are lots of clever techniques to get around the plethora of scenarios that you may want to implement. This is just a first step on the direction of really understanding CFD!
