NACA 0012 Airfoil CFD Solver

A 2D incompressible Navier–Stokes solver for laminar flow past a NACA 0012 airfoil, implemented from scratch in Python using a fractional-step (projection) method on a staggered Cartesian grid.

Overview

This solver was built as a computational physics project to study the core numerical methods behind incompressible CFD. The simulation produces physically reasonable velocity fields, streamlines, and vorticity distributions consistent with classical laminar airfoil aerodynamics.

Key features:

• Staggered grid (MAC layout) for pressure–velocity coupling
• Fractional-step (projection) method for time integration
• NACA 4-digit airfoil geometry generation with cosine spacing
• Ray-casting airfoil masking for no-slip boundary enforcement
• Configurable angle of attack, Reynolds number, and grid resolution

Physics

The solver integrates the incompressible Navier–Stokes equations:

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

Numerical scheme:

1. Predictor step — advance momentum without pressure gradient to get intermediate velocity $\mathbf{u}^*$
2. Pressure solve — solve Poisson equation $\nabla^2 p = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}^*$ via Gauss–Seidel iteration
3. Projection step — correct velocity: $\mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t , \nabla p$

Convective terms use centered finite differences; diffusion uses a second-order Laplacian.

Results

Velocity Magnitude	Streamlines	Vorticity
Flow accelerates around the leading edge	Streamlines remain attached	Vorticity concentrated at surface

The simulation at Re ≈ 1000 and 6° angle of attack shows:

• Smooth flow acceleration around the leading edge and along both surfaces
• Attached streamlines with no far-field recirculation
• Vorticity layers of correct sign concentrated near the airfoil surface

Known limitation: The Poisson solver does not fully converge to machine precision within the iteration budget, leaving residual divergence. Dominant flow features are captured correctly, but strict incompressibility is not enforced. See Discussion in the paper for details.

Getting Started

Requirements

numpy
matplotlib

Install with:

pip install numpy matplotlib

Running the simulation

Open airfoil_flow.ipynb in Jupyter and run all cells in order. The notebook is structured as follows:

Cell	Contents
1	Imports
2	Configuration parameters
3	Mesh generation and boundary conditions
4	Airfoil geometry and masking
5	Navier–Stokes solver functions
6	Simulation setup and execution
7–9	Visualization (velocity, streamlines, vorticity)

Note: The default grid (400×200) and 500 time steps take several minutes to run due to the Python loop in the convection and masking routines. Reduce Nx, Ny, or nt for a faster test run.

Key parameters

Defined in the configuration cell:

Parameter	Default	Description
Nx, Ny	400, 200	Grid resolution
dt	0.0005	Time step
nt	500	Number of time steps
nu	0.001	Kinematic viscosity (Re ≈ 1000)
alpha	-0.1047 rad	Angle of attack (−6°)
inlet_velocity	1.0	Freestream velocity

Project Structure

.
├── airfoil_flow.ipynb   # Main simulation notebook
├── README.md
└── Reiter_CP2025_Final_Project_CFD.pdf  # Written report

Potential Improvements

• Replace Gauss–Seidel with SOR or a multigrid Poisson solver for better divergence control
• Vectorize the convection and ray-casting loops for significant speedup
• Use a body-fitted grid to eliminate stair-stepping errors at the airfoil surface
• Add angle-of-attack sweep and lift/drag coefficient computation
• Extend to unsteady simulation and higher Reynolds numbers (vortex shedding)

References

• Chorin, A. J. (1968). Numerical solution of the Navier-Stokes equations. Mathematics of Computation, 22(104), 745–762.
• Abbott, I. H., & Von Doenhoff, A. E. (1959). Theory of Wing Sections. Dover.
• NACA 4-digit airfoil formulation: airfoiltools.com