NACA 0012 CFD: From Commercial CFD to a C++ Solver From First Principles

One aerofoil, a NACA 0012 at 10° incidence and a chord Reynolds number of 6×10⁶, and two goes at it. The first was in ANSYS Fluent, validated against NASA surface-pressure data, and it answered whether the result could be trusted. The second was a solver written from scratch, and it exists because the first study kept raising a question it had no way of answering: what was happening between the mesh and the answer.

That precursor is the short section immediately below. The rest of the page is the solver. It implements the numerical machinery a commercial code keeps out of sight: the finite-volume discretisation, the pressure-velocity coupling that makes incompressible flow solvable at all, the correction that keeps a skewed mesh honest, and a two-equation turbulence model with its near-wall treatment. There is no numerical library underneath it. The linear systems are solved by a Gauss-Seidel sweep with over-relaxation written for this project, and no discretisation, no boundary condition and no turbulence closure comes from a library.

The discipline that makes it worth showing is not that it runs. It is that every stage was checked against a problem whose answer was already known, before being trusted on a problem whose answer was not. At the end the two solutions meet: the same section, the same incidence, the same Reynolds number, one answer from a commercial code and one from a few thousand lines of C++.

Before any of the code below existed, this aerofoil was solved in ANSYS 2026 R1: two-dimensional steady RANS, 10° incidence, chord Reynolds number 6×10⁶, on a mesh of roughly 27,000 cells with ten inflation layers and a sphere of influence around the section. The case was chosen because high-quality NASA experimental data exists for it, so the argument could be carried to the end instead of stopping at plausibility.

It worked. The computed surface pressure distribution overlaps the NASA data across the full chord, and the integrated lift came out at CL ≈ 1.06 against a measured 1.07 to 1.08, about 1.4 % low, having been predicted at 1.097 by a thin-aerofoil hand calculation done before the solver was opened. Matching the whole Cp distribution is the stronger claim: an integrated coefficient can agree through two errors cancelling, a point-by-point match across the chord cannot.

So the study answered its own question: the result is trustworthy for lift and not for drag, and it can say which is which. What it could not answer is what happened between the mesh and the answer. Choosing a second-order scheme, a coupled solver or a y⁺ target was still selecting from a menu. Everything below is what came of writing the thing each of those options switches on.

The solver is cell-centred and collocated: velocity components and pressure all live at cell centres. Every transport equation is integrated over a control volume and converted by the divergence theorem into a sum of fluxes through faces, so the discrete equations conserve mass and momentum by construction rather than by accident.

The aerofoil is generated analytically from the NACA 0012 thickness equation with the closed-trailing-edge coefficient, distributed with cosine spacing so that points cluster where the curvature is: the leading edge and the trailing edge.

The first O-grid ran. It produced a converged Laplace solution and had no inverted cells, which is exactly the situation in which a mesh gets used without further thought. Because the generator also computed quality metrics, it was measured instead, and the measurement was bad enough to send the generator back for a rewrite.

Verification asks a question that has nothing to do with aerodynamics: given these equations, does the code solve them correctly, and at the rate the discretisation promises? It is answered against problems with closed-form solutions, where any discrepancy is the code's fault rather than the physics model's.

Once the operators were generalised to the body-fitted grid, the geometric machinery was tested on its own terms: face ownership and connectivity, cell areas, the closure of face vectors around every cell, mass conservation under a uniform flow imposed exactly, and the Green-Gauss gradient of a linear pressure field with a known answer.

SIMPLE was not written directly onto the aerofoil grid. It was built first on the lid-driven cavity, a square box with a sliding top, on a staggered Cartesian arrangement where the pressure-velocity coupling is simplest and the failure modes are best understood. The cavity at Re = 100 has been a reference case since Ghia, Ghia and Shin published tabulated centreline velocities in 1982.

Before turbulence, the collocated solver was exercised on the aerofoil at Re = 100, where the flow is steady and laminar and the equations being solved are unambiguous. This is where Rhie-Chow interpolation and the non-orthogonal correction were shaken out.

An earlier attempt to run the momentum equations on their own, with pressure frozen, is worth recording because of how it failed. Without a pressure correction enforcing continuity there is nothing to stop the velocity field growing without bound, and the run diverged to non-finite values in every one of the 15,800 cells. That is the correct behaviour for an incompressible momentum predictor with no coupling, and it is the most direct demonstration of why SIMPLE exists.

Moving to a chord Reynolds number of six million changes the problem entirely. The boundary layer becomes a thin turbulent sheet, and the solver has to gain a turbulence model and a mesh capable of resolving that sheet without wall functions.

Before running the case of interest, the turbulent solver was run at zero incidence. A symmetric section at zero angle of attack must produce exactly zero lift, so any lift the solver reports is pure numerical and modelling error, with no physics hiding it. It is the cheapest severe test available for this configuration.

The turbulent case was solved twice on the same mesh with the same turbulence model, changing only how convected quantities are evaluated at faces. The first-order upwind solution converged to every tolerance and looks entirely reasonable. It is also wrong in a way that only a comparison exposes.

The pressure drag falls by a factor of nearly six. That number is the numerical diffusion of the first-order scheme made visible: smearing gradients across cells thickens the effective boundary layer, weakens the suction peak and leaves a pressure imbalance that is a discretisation artefact rather than aerodynamics.

The peak turbulent viscosity ratio tells the same story from the model's side. The first-order solution carries almost twice the peak eddy viscosity of the second-order one. Numerical diffusion thickened the shear layer, the turbulence model saw the thickened profile and responded by producing more turbulence, and the two errors reinforced each other.

The NASA Turbulence Modeling Resource publishes this exact case, a NACA 0012 at 10° and Re = 6×10⁶, solved with the SST model by three independent codes on a common 897 × 257 grid. Because the three agree closely with each other, they form a usable reference band.