Custom 2D Airfoil Design Code

The direct route — geometry in, pressures out — is an analysis method pressed into service as a design method. It works, but every iteration is guesswork: thickening the nose to soften a suction peak also moves the peak, changes the recovery, shifts transition. The designer is steering through the geometry, one step removed from the quantity that actually matters.

Inverse design removes that indirection. Laminar flow is wanted over the front 40% of the upper surface? Then prescribe a gently accelerating velocity there — a favourable pressure gradient holds the boundary layer laminar almost by definition. The recovery region should avoid separation? Choose a recovery shape that stays safely below the separation criterion. Each aerodynamic requirement maps directly onto a feature of the prescribed distribution, and the geometry that results is, by construction, the section that meets it.

The catch is mathematical: not every velocity distribution corresponds to a physically valid, closed airfoil. Handling that constraint properly — rather than hoping — is where most of the theory in this project lives.

Full viscous flow is governed by the Navier–Stokes equations — accurate and, for a fast design loop, hopelessly expensive. Dropping viscosity gives the Euler equations. But for an airfoil well below stall one can go a step further: outside the paper-thin boundary layer the flow is not only effectively inviscid but also irrotational, and at low Mach numbers incompressible. Under those three assumptions the whole flow field reduces to a single linear equation — Laplace's equation for the velocity potential.

That reduction is what makes the method fast and, more importantly, exact: potential flow around a 2D body can be solved in closed form, with no grid, no discretisation error and no turbulence model in the outer flow. The price is equally clear, and it defines the honest operating envelope of the code: strictly incompressible (trustworthy up to roughly Mach 0.3, stretchable with a compressibility correction, never transonic — no shocks, no wave drag) and inviscid in the outer field, so all viscous physics has to be reintroduced separately through a boundary-layer model.

Viscosity, though, is exactly where drag, transition and stall live. So the architecture is a classic two-layer one: an exact potential-flow solution for the pressure field, coupled to an integral boundary-layer method that runs along the surface afterwards. Chapters 08–10 cover that second layer.

The classical trick for exact 2D potential flow: the flow around a circle is known analytically, and any smooth airfoil-like contour can be related to a circle by a conformal map — an angle-preserving transformation of the complex plane. Map the circle solution through that transformation and you have the exact flow around the airfoil, at any angle of attack, including the circulation that sets the lift via the Kutta condition at the trailing edge.

For inverse design this becomes even more elegant. Points on the circle are identified by their polar angle ν, running from 0° at the trailing edge over the upper surface to 360° back at the trailing edge. The designer prescribes the velocity as a function of ν, and the mapping that turns the circle into the airfoil is derived from that prescription — the surface velocity distribution and the geometry are two faces of the same analytic function.

Working in the circle plane also solves a practical annoyance for free: points equally spaced in ν crowd together near the leading edge in the airfoil plane, exactly where the curvature is highest and resolution is most needed. The 66-point contour the code outputs is dense at the nose and coarse on the flat parts without any explicit grading logic.

Here is the mathematical heart of inverse design: a prescribed velocity distribution is only realisable if it satisfies three integral constraints. The mapping must return to its starting point (two conditions, one per coordinate — otherwise the 'airfoil' is an open spiral), and the flow far away must reduce to a clean uniform stream. Prescribe velocities freely and all three will, in general, be violated.

The code therefore treats a handful of the prescribed quantities not as user inputs but as unknowns — among them the closure contributions controlled by ω and an auxiliary constant K_S — and solves a small system of transcendental equations so that the constraints come out exactly. This is done with an iteration that typically converges in three to five steps, after which the design angles are re-normalised (which is why the α* values shift slightly between the first and the converged pass).

It's a genuinely satisfying piece of numerics to watch: the residual K_S drops from order one to below the rounding threshold in four iterations, and out comes a closed contour — in the example run, the 19.3%-thick section shown in chapter 01, with its zero-lift angle (4.27°) and zero-lift pitching moment (−0.101) delivered as by-products of the same solution.

The potential-flow solution knows nothing about drag — d'Alembert's paradox guarantees it. All viscous physics is confined to the boundary layer, and because that layer is thin, it doesn't need a 2D flow solver of its own: its state at each surface station can be summarised by a few integral thicknesses, and their evolution follows ordinary differential equations driven by the outer velocity distribution.

Concretely, the method marches from the stagnation point along each surface, integrating the momentum-integral equation (von Kármán) together with an energy-integral equation. Two quantities carry the state: the momentum thickness, which measures accumulated momentum loss, and the shape factor built from the energy and momentum thicknesses, which measures how 'healthy' the velocity profile is. Laminar and turbulent stretches use different closure relations for skin friction and dissipation, calibrated against classical boundary-layer data.

Because the boundary layer is fed directly by the designed velocity distribution, cause and effect stay legible: steepen the recovery with μ and you watch the shape factor dive towards its separation limit; extend a favourable gradient with a higher α* and transition walks aft. This tight coupling between the design parameters and the boundary-layer response is what makes the tool educational rather than a black box.

With the boundary layer integrated to the trailing edge, drag follows from momentum accounting. Far downstream, the profile drag of the section equals the momentum deficit carried by the wake — but the boundary-layer solution stops at the trailing edge, where the deficit hasn't finished developing. The classical bridge is the Squire–Young formula: it extrapolates the momentum thickness at the trailing edge to the far wake, using the local trailing-edge velocity and an assumption about how the wake relaxes towards freestream pressure.

Upper and lower surface each contribute their share; the bubble penalties from chapter 09 are added on top. Lift and pitching moment come from the potential-flow circulation, corrected for the boundary layer's decambering effect — the displacement of the outer flow slightly de-loads the aft airfoil, which is why the delivered lift-curve slope sits below the ideal 2π and why the code applies a lift-efficiency correction rather than reporting bare inviscid values.

None of this is exotic — these are the standard tools of classical airfoil analysis — but implementing them yourself, and watching the drag polar respond to a single slider nudge in the recovery shape, teaches more section aerodynamics than any number of canned CFD runs.

The numerical core alone is a batch program: edit a text input file, run, read a few hundred lines of tabulated output. Usable, but it throttles exactly what inverse design is good at — rapid, intuition-building iteration. So the second half of the project was wrapping the solver in a GUI that closes the loop.

The interface exposes the aft-design parameters as sliders — λ, μ, ω for each surface, plus the two ramp controls — and one button: run. On each run the tool writes the parameters into the input file, executes the solver, parses the full output, and refreshes four linked plots: velocity distributions, lift curve with transition/separation overlay, moment curve, drag polar — plus a live preview of the section itself. Standard file handling (new/open/save/save-as, modified-state tracking), persistent settings and a dark mode round it off.

The whole cycle takes about a second, and that speed changes how the tool is used. Sweeping μ from −0.5 to +1 in ten steps and watching the polar's low-drag bucket widen, shift and finally collapse into separation teaches the recovery trade-off viscerally — the kind of understanding that batch runs with minutes of turnaround never produce.

The first complete version was Python: the numerics prototyped where they're quickest to get right, the interface in Tkinter with embedded Matplotlib figures, pandas doing the output wrangling. Python was the right choice for working out the method — every intermediate quantity is one print statement away, and refactoring the parametrisation cost minutes instead of hours.

For the final version the code moved to C++. Partly for raw speed — parameter sweeps and future optimisation loops want the solver call to be milliseconds, not the better part of a second — and partly for robustness: a compiled, statically-typed numerical core with no interpreter and no dependency environment is simply a nicer thing to ship and to keep working years later. The Python version survives as the reference implementation the C++ results were validated against, which proved its worth more than once when a subtle indexing bug in the port produced polars that were almost, but not quite, right.

The two-stage approach is one of the main process lessons of the project: prototype the physics in the language where thinking is cheap, port to the language where running is cheap, and keep the prototype alive as the ground truth.

It's worth being precise about what physics is actually in the loop. Captured: the exact inviscid pressure field at any incidence, including the correct nonlinear behaviour near the leading edge that thin-airfoil theory misses; laminar and turbulent boundary-layer growth driven by that field; natural transition as a function of Reynolds number and pressure-gradient history; laminar separation bubbles (flagged, with empirical drag penalty); turbulent trailing-edge separation onset; profile drag from the wake momentum deficit; lift and moment including the boundary layer's decambering correction; and the dependence of all of the above on Reynolds number.

That list covers the phenomena that dominate section behaviour for essentially the whole low-speed envelope — sailplanes, propellers, UAVs, wind-turbine blades, general aviation. It's the same modelling level that generations of proven laminar airfoil families were designed with, which was exactly the ambition: not a toy, but the real classical toolchain, built end to end.

Every simplification in chapter 03 has a bill attached, and part of the project was learning to state it precisely. Compressibility: the flow model is strictly incompressible. Below Mach ≈ 0.3 the error is negligible; towards Mach 0.5–0.6 a Prandtl–Glauert or Kármán–Tsien correction could stretch the pressure predictions, but the moment the local flow goes supersonic the model is silently, fundamentally wrong — no shocks, no wave drag, no drag divergence. Transonic design is out of scope by construction, not by tuning.

Viscous coupling: the boundary layer listens to the potential flow, but the potential flow never hears back (beyond the global decambering correction). This one-way coupling is accurate while the layer stays thin and attached — and progressively fiction once significant separation develops. Consequently the code cannot predict maximum lift or stall shape; it can only warn where separation begins. Bubbles are flagged and penalised, not resolved, so at low Reynolds numbers (below roughly 10⁵, where bubbles dominate) the drag numbers become qualitative. Transition is an empirical correlation, valid for smooth surfaces in low-turbulence flow — roughness or a dirty environment shifts it in reality but not in the model, except by forcing transition manually.

And structurally: everything is 2D, steady and single-element. No sweep, no finite-wing effects, no flaps in the model, no unsteady aerodynamics. None of these limits is a flaw — each is the direct price of the one-second design loop — but the difference between a dangerous tool and a useful one is that the boundaries are known, stated, and respected.

The project delivered exactly the intuition it was started for. Concepts that stay abstract in lectures — why a favourable gradient delays transition, why an aggressive recovery buys thickness but risks the bucket, why moment and camber trade against each other, why low-Reynolds sections need their recovery softened — became things felt through the sliders, because every one of them is a visible cause-and-effect in the four plots.

It also sharpened engineering judgement about models in general: an exact solution of the wrong equations (pure potential flow) is useless for drag, while a stack of well-chosen approximations — integral boundary layer, empirical transition, wake extrapolation — predicts real polars remarkably well inside a known envelope. Fidelity is not one number; it's an itemised list.

With a solver call this cheap, the natural next step is to stop moving the sliders by hand: wrap the design loop in an optimiser and search the λ–μ–ω–α* space directly against polar objectives — widest bucket at a target cl, best endurance figure of merit, mildest moment. A compressibility correction for the high-subsonic end and a proper bubble drag model for the low-Reynolds end are the two most valuable physics upgrades on the list.

The code has already earned its keep beyond its own page: it became the workhorse behind the Horten H3 airfoil modernisation project, where the design loop ran against a historic flying wing's stability constraints. That's the quiet test of a self-built tool — whether it gets reached for when the next real problem shows up.