CFD with Python

From scratch

A computational fluid dynamics solver built from scratch in Python: a direct numerical simulation of 2D incompressible flow past a cylinder, to watch vortex shedding (the Von Kármán vortex street) appear.

• It solves the 2D incompressible Navier–Stokes equations at $Re = 200$, in a channel $20D$ long and $10D$ wide with the cylinder at one fifth of the length.
• Chorin's projection method on a staggered grid, with the cylinder added by an immersed boundary method.
• Numba compiles the grid loops to parallel machine code.
• The run saves velocity and vorticity frames, a video of the simulation, and a dataset (.npz) for later surrogate modelling with a PINN.
• The repo has the explanation (Explanation.ipynb) and the code (CFD.ipynb) side by side.

The cell below runs the same method live in your browser. Drag the cylinder to move it.

The equations

With the free-stream speed $U$, the cylinder diameter $D$ and the density $\rho$ all set to 1, the equations are already non-dimensional, and the only parameter left is the Reynolds number $Re = UD/\nu = 1/\nu$:

$$\nabla \cdot \mathbf{u} = 0, \qquad \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\,\mathbf{u} = -\nabla p + \frac{1}{Re}\,\nabla^{2}\mathbf{u} + \mathbf{f}$$

The first says the fluid is incompressible (no velocity divergence). The second is Newton's second law for each bit of fluid: it is pushed by pressure differences, slowed by viscosity, and pushed by $\mathbf{f}$, the force that stands in for the cylinder (see below).

A staggered grid

The $20 \times 10$ domain is cut into $384 \times 192$ cells ($h \approx 0.052$). Pressure $p$ lives at cell centres, the horizontal velocity $u$ on the left and right faces, and the vertical velocity $v$ on the top and bottom faces. This staggered (MAC) layout puts each velocity exactly between the two pressures that push it, so the pressure gradient and the divergence of a cell are plain differences of neighbours:

$$\left.\frac{\partial p}{\partial x}\right|_{i+\frac12,j} = \frac{p_{i+1,j} - p_{i,j}}{\Delta x}, \qquad (\nabla\cdot\mathbf{u})_{i,j} = \frac{u_{i+\frac12,j} - u_{i-\frac12,j}}{\Delta x} + \frac{v_{i,j+\frac12} - v_{i,j-\frac12}}{\Delta y}$$

Diffusion uses the 5-point stencil $\nabla^2 u \approx \dfrac{u_{\text{right}} + u_{\text{left}} - 2u}{\Delta x^2} + \dfrac{u_{\text{top}} + u_{\text{bottom}} - 2u}{\Delta y^2}$, and advection ($u\,\partial u/\partial x + v\,\partial u/\partial y$) uses first-order upwind differences: each derivative looks at the neighbour the flow is coming from. Where a velocity component isn't stored (for example $v$ at a $u$ face), it is the average of the four nearest values.

One time step: Chorin's projection

Stepping the momentum equation forward on its own would not keep the flow divergence-free, and it gives no way to update the pressure. Chorin's projection method splits each step in three:

1. Predict a velocity $\mathbf{u}^{\ast}$ from advection, diffusion and the cylinder's force, ignoring pressure: $\mathbf{u}^{\ast} = \mathbf{u}^{n} + \Delta t\left(-(\mathbf{u}^{n}\cdot\nabla)\mathbf{u}^{n} + \tfrac{1}{Re}\nabla^{2}\mathbf{u}^{n}\right)$.
2. Solve for the pressure that removes the divergence of $\mathbf{u}^{\ast}$. Asking for $\nabla\cdot\mathbf{u}^{n+1} = 0$ in step 3 gives a Poisson equation: $\nabla^{2}p^{\,n+1} = \dfrac{\nabla\cdot\mathbf{u}^{\ast}}{\Delta t}$.
3. Correct the velocity with the pressure gradient: $\mathbf{u}^{n+1} = \mathbf{u}^{\ast} - \Delta t\,\nabla p^{\,n+1}$, which is now divergence-free.

The Poisson equation becomes a linear system $A\mathbf{p} = \mathbf{b}$ with one row per cell (73,728 of them). $A$ never changes, so it is built once as a sparse matrix and factorised once with SciPy; each step then only solves for a new $\mathbf{b}$. Stored as a dense NumPy array, $A$ would need 40.5 GiB.

The walls have $\partial p/\partial n = 0$ (each wall cell copies its pressure into a ghost cell, which just adds the missing neighbour's coefficient to the diagonal), and one cell is pinned to $p = 0$, because pressure is only defined up to a constant.

The cylinder: an immersed boundary

The grid is a plain rectangle and never changes. The cylinder is 200 points on its surface, $(X_k, Y_k) = (C_x + R\cos\theta_k,\; C_y + R\sin\theta_k)$, and each step:

1. Interpolate the velocity at every surface point from the four nearest grid values (bilinear weights).
2. Compute the force that would stop the fluid there (the cylinder is still): $f_k = \alpha\,\dfrac{0 - u_k}{\Delta t}$, with a relaxation factor $\alpha = 0.8$ so interpolation errors don't produce excessive forces.
3. Spread each force back onto the same four grid points with the same weights, and add it to $\mathbf{u}^{\ast}$.

Because the cylinder is only a list of points, moving it means moving the points. That is why you can drag it around in the live cell.

Boundaries and stability

The fluid comes in from the left at speed 1 and a small angle of $2.5°$, so the flow isn't perfectly symmetric and the wake can start to oscillate. The walls are frictionless (no flow through them), and the outlet is open: velocities don't change across it.

The explicit time step has two limits, from diffusion and from advection:

$$\Delta t \le \min\!\left(\frac{Re\,h^{2}}{4},\; \frac{h}{2\max(u, v)}\right) = \min(0.1352,\; 0.0026) = 0.0026$$

(taking $\max(u, v) = 10$ as a worst case). The notebook uses $\Delta t = 0.001$ to stay safely inside it and to capture the shedding smoothly.

Results: lift, drag and the Strouhal number

The force the cylinder exerts on the fluid, summed over its surface points, gives the drag and lift, and their coefficients $C_d = 2f_{\text{drag}}$ and $C_l = 2f_{\text{lift}}$ (since $\tfrac12\rho U^2 = \tfrac12$). As vortices peel off one side and then the other, $C_l$ oscillates. Its frequency $f$, found with a fast Fourier transform, gives the Strouhal number $St = fD/U$.

The empirical relation $St = 0.2665 - 1.018/\sqrt{Re}$ predicts about $0.194$ at $Re = 200$. The notebook measures $0.159$, which matches $Re \approx 80$. The reason is the first-order upwind scheme: its truncation error acts like extra, numerical viscosity, so the flow behaves as if the Reynolds number were lower:

$$\frac{1}{Re_{\text{eff}}} = \frac{1}{Re} + \nu_{\text{numerical}}$$

The live cell shows the same effect. Its grid is coarser (256 × 128, $h \approx 0.078$), so it has more numerical viscosity and measures about $St \approx 0.145$ at $Re = 200$ (watch St in the Variables pane). Raise the Reynolds number and it creeps towards the textbook value.

A second method: lattice-Boltzmann

The cell below is a different solver for the same flow, written for this site to compare (it is not in the notebook). Instead of velocities and pressure, the lattice-Boltzmann method tracks how many fluid particles move in each of 9 directions at every grid point (D2Q9). Each step they stream to the neighbouring point and collide, relaxing towards equilibrium at a rate set by the viscosity. The cylinder is cells that bounce particles back. Pressure is never solved for, which makes each step very cheap.

References

From the notebook's sources:

1. C. H. K. Williamson and G. L. Brown, A series in $1/\sqrt{Re}$ to represent the Strouhal–Reynolds number relationship of the cylinder wake, Journal of Fluids and Structures, 1998.
2. A. Placzek, J.-F. Sigrist and A. Hamdouni, Numerical simulation of vortex shedding past a circular cylinder in a cross-flow at low Reynolds number with finite volume technique, ASME Pressure Vessels and Piping Conference, 2007.
3. J. Kim, D. Kim and H. Choi, An immersed-boundary finite-volume method for simulations of flow in complex geometries, Journal of Computational Physics, 2001.
4. F. H. Harlow and J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Physics of Fluids, 1965.
5. Flow Science, What are artificial and numerical viscosities?