(CFD) is usually seen as a black box of complex business software. Nonetheless, implementing a solver “from scratch” is one of the crucial powerful ways to learn the physics of fluid motion. I began this as a private project, and as a component of a course on Biophysics, I took it as a possibility to finally understand how these beautiful simulations work.
This guide is designed for data scientists and engineers who wish to move beyond high-level libraries and understand the underlying mechanics of numerical simulations by translating partial differential equations into discretized Python code. We may even explore fundamental programming concepts like vectorized operations with NumPy and stochastic convergence, that are essential skills for everybody keen on broader scientific computing and machine learning architectures.
We are going to walk through the derivation and Python implementation of an easy incompressible Navier-Stokes (NS) solver. After which, we’ll then apply this solver to simulate airflow around a bird’s wing profile.
The Physics: Incompressible Navier-Stokes
The basic equations of CFD are the Navier-Stokes equations, which describe how velocity and pressure evolve in a fluid. For regular flight (like a bird gliding), we assume the air is incompressible (constant density) and laminar. They might 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 may add in gravity, mechanical stresses, and even electromagnetism in case you’re feeling prefer it. The writer can attest to this being very much not recommend for a primary 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 primary equation (Momentum) balances inertia against pressure gradients and viscous diffusion. The second equation (Continuity) enforces that the fluid density stays constant.
The Pressure Coupling Problem
A significant challenge in CFD is that pressure and velocity are coupled: the pressure field must adjust continuously to make sure the fluid stays incompressible.
To resolve 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 each timestep to update the pressure, ensuring the speed field stays divergence-free.
Discretization: From Math to Grid
To resolve these equations on a pc, 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 instance, 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 4 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 want to mark which grid points lie the solid wing.
- We load a wing mesh (e.g., from an STL file).
- We create a Boolean mask array where
Trueindicates a degree contained in the wing. - In the course of the simulation, we force velocity to zero at these masked points (no-slip/no-penetration condition).
3. The Predominant Solver Loop
The core loop repeats until the answer reaches a gentle state. The steps are:
- Construct the Source Term (b): Calculate the divergence of the speed terms.
- Solve Pressure: Solve the Poisson equation for p using Jacobi iteration.
- Update Velocity: Use the brand new pressure to update u and v.
- Apply Boundary Conditions: Implement inlet velocity and 0 velocities contained in the wing.
The Code
Here is how the core mathematical updates look in Python (vectorized for performance).
Step A: Constructing 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 speed using the express 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 relentless far-field inflow.
Qualitative Observations The outcomes align with physical expectations. The simulations show high pressure beneath the wing and low pressure above it, which is strictly the mechanism that generates lift. Velocity vectors show the airflow accelerating excessive surface (Bernoulli’s principle).
Forces: Lift vs. Drag By integrating the pressure field over the wing surface, we are able to calculate lift.
- The solver demonstrates that pressure forces dominate viscous friction forces by an element of nearly 1000x in air.
- Because the angle of attack increases (from 0∘ to −20∘), the lift-to-drag ratio rises, matching trends seen in wind tunnels and skilled 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 effective flow details.
Where to go from here? This project serves as a place to begin. Future improvements could include implementing higher-order advection schemes (like WENO), adding turbulence modeling, or moving to Finite Volume methods (like OpenFOAM) for higher mesh handling around complex geometries. There are plenty of clever techniques to get across the plethora of scenarios that it’s possible you’ll wish to implement. That is just a primary step on the direction of really understanding CFD!
