Writing · Flexcompute Foil · Part II
The Geometry
Problem
A solver that matches XFOIL to four decimal places on clean airfoils can completely fall apart with a five-degree control surface deflection. The fix wasn’t in the solver. It was in how I prepared the geometry.
Harry Smith · Part I
Part I covered the viscous solver — the boundary layer coupling, the Newton iteration, the eight bugs that kept it from matching XFOIL. I got the solver to four-decimal-place parity on clean airfoils. Then I deflected a flap five degrees and the whole thing fell over.
Zero out of seventeen alpha points converged. Not “a few points diverged near stall.” Zero. Every single operating point. Every flap deflection I tested — positive, negative, five degrees, twenty degrees. The solver was producing garbage.
The first thing I did was the diagnostic that saved me weeks of looking in the wrong place. I took XFOIL’s own panelled geometry — after it had run its own GDES FLAP and PANE commands — saved those 160 coordinates, loaded them into Rustfoil, and ran the polar. Seventeen out of seventeen converged. The solver was fine. The geometry was the problem.
What GDES FLAP Actually Does
If you’ve used XFOIL, you know the flap command. Type GDES, type FLAP, enter a hinge location and a deflection angle, type X to execute, go back to the main menu, type PANE to repanel. Takes about five seconds. Looks like “rotate the points aft of the hinge.”
It isn’t. That’s what I tried first, and it’s why nothing converged. What XFOIL actually does is roughly four hundred lines of Fortran that I ended up reading three times before I understood it. Here’s what the naive “split and rotate” looks like — drag the deflection slider and watch what happens at the hinge:
Flap Deflection
NACA 0012scroll to zoom · drag to pan
See the problem? One side folds back on itself. The other side opens a gap. Toggle naive spline and zoom to the hinge to see what happens next — the red curve forces C2 continuity through the corner and oscillates, creating a tiny x-reversal that kills the solver. Toggle GDES fix to see what XFOIL’s procedure actually produces: fold trimmed, break points placed, clean geometry. XFOIL handles both of these with a careful procedure:
First, GETXYF uses Newton inversion on the buffer spline (SINVRT) to find the exact arc-length positions where x = hinge on the upper and lower surfaces. Then it determines which surface folds (disappearing segment) and which opens (gap). For a trailing-edge-down deflection with the hinge on the camber line, the bottom surface folds and the top opens.
Then SSS — 150 lines of Newton iteration on the spline — finds the exact break arc-lengths. On the fold side, it finds two arc-length values S1 and S2 that bracket the segment which “disappears” when the flap rotates (the included angle between the two hinge-to-break vectors equals the deflection angle). On the gap side, it finds the single arc-length where the hinge-to-surface vector is perpendicular to the surface — the point where the surface has to break to let a gap open.
After finding the breaks, XFOIL places explicit corner points with adjacent helper nodes at 33% spacing from the break to the next existing node. It rotates the flap, deletes the disappeared segment on the fold side, adds circular-arc fill points on the gap side, runs SCHECK to remove splinter segments (micro-panels that occur when a break lands nearly on top of an existing node), and re-splines the whole buffer.
Only then does PANE run PANGEN to distribute the final 160 panel nodes by curvature. The buffer it operates on has already been cleaned up by the FLAP procedure. This is the step I was skipping entirely.
The Spline Problem
Here’s what I tried first: split the coordinates at the hinge x-position, rotate the aft points, detect any self-intersection on the fold side and trim it, then repanel through the cubic spline. Seemed reasonable. A flap is just a rotation, right?
The problem is that a flap deflection creates a slope discontinuity at the hinge. The surface on the fixed side approaches the break from one direction; the rotated surface leaves in another. A standard cubic spline enforces C2 continuity at every interior knot — it forces the second derivative to match on both sides of every point, including the break. When the data has a genuine slope discontinuity, the spline can’t represent it. It oscillates through the corner (essentially a Gibbs phenomenon for splines).
When I dumped the panelled coordinates near the hinge, the x-values were going backwards. Not by a lot — a reversal of about 0.002 chord — but enough to create a tiny self-intersecting loop in the surface. The inviscid solver would produce a pressure spike there, and the BL equations would just diverge trying to march through it.
I only knew to look for this because I’d written my own cubic splining routines before and been bitten by the same issue in a different context. If you’ve ever tried to fit a cubic spline through data with a genuine kink and watched the interpolated values overshoot wildly on either side, you know the feeling. The spline is doing exactly what it was told to do — enforce C2 continuity — but the data isn’t C2, and the resulting oscillation is catastrophic for a boundary layer solver that needs monotonic surface coordinates.
The key diagnostic was simple and I wish I had run it first. Feed XFOIL’s own geometry to Rustfoil’s solver. Seventeen out of seventeen converged for every deflection. The solver was fine. The geometry pipeline was the entire problem.
The Fix, Procedurally
I tried several things before I got it right (naturally). First I added explicit break points at the hinge in the input geometry, hoping the spline’s curvature-based node bunching would handle the rest. It didn’t. Then I implemented corner detection in the spline — scanning for isolated tangent-angle spikes and breaking the tridiagonal solve there, matching XFOIL’s SEGSPL. That helped for moderate deflections but failed at small angles (corner below the detection threshold) and at large angles (fold geometry was wrong).
What finally worked was what should have been obvious from the start: understand the actual procedure and recreate it. I read XFOIL’s FLAP subroutine until I understood what every section was doing and why, then rebuilt it in Rust. SSS with its Newton iteration. SINVRT for spline inversion. The index management for deleting disappeared points and inserting new ones. SCHECK for splinter removal. The fold-versus-gap logic with its 33%-spacing helper points.
One pragmatic divergence: I skip XFOIL’s arc-fill on the gap side. XFOIL adds circular-arc points to smoothly bridge the opening between the fixed and rotated surfaces. These arc-fill points create a temporary x-reversal in the buffer (they trace a circular arc that curves backward), and XFOIL’s SEGSPL handles this because it allows derivative discontinuities at segment joints. My repaneller uses a standard cubic spline that can’t handle the reversal without oscillating. So I leave the gap empty and let the repaneller’s curvature-based node placement bridge it naturally. Works fine.
The result is a 450-line Rust module (flap.rs) that matches XFOIL’s FLAP procedure closely enough for the solver to converge at every deflection we tested.
Convergence
NACA 0012, hinge at 75% chord, Re = 106, alpha from −4° to 12° in 1° steps, 17 operating points per deflection.
| δ | XFOIL 6.99 | Rustfoil |
|---|---|---|
| −25° | 7/17 | 15/17 |
| −20° | 14/17 | 17/17 |
| −15° | 6/17 | 17/17 |
| −10° | 16/17 | 17/17 |
| −5° | 17/17 | 17/17 |
| +5° | 16/17 | 17/17 |
| +10° | 16/17 | 17/17 |
| +15° | 16/17 | 16/17 |
| +20° | 16/17 | 17/17 |
| +25° | 15/17 | 15/17 |
Before the GDES port, every cell in the Rustfoil column was 0/17.
The Polars
Eleven flap deflections, −25° to +25° in 5° steps, alpha from −6° to +16° in 0.5° steps. XFOIL 6.99 is dashed, Rustfoil is solid. Every comparison uses the same Reynolds number, the same N-crit, the same alpha sequence.
Flap Polars Across Eleven Deflections
CL vs α for the NACA 0012 with flap deflections from −25° to +25° in 5° steps, hinge at 75% chord, Re = 10⁶. XFOIL is dashed, Rustfoil is solid. For moderate deflections (±10°) the curves lie on top of each other. At the extremes you can see where XFOIL gives up and Rustfoil keeps going — but also where the two start to diverge near stall.
Drag Polars
CL vs CD tells you whether the solver is getting the viscous drag right, not just the lift. The drag bucket shape and the drag rise at high CL are both sensitive to boundary layer transition and the pressure recovery aft of the hinge. For ±10° the agreement is very tight. For larger deflections the drag numbers start to diverge near CL_max, which is where the geometry difference (I skip XFOIL’s arc-fill; it doesn’t) shows up most.
Error Envelopes
The honest view. ΔCL and ΔCD between Rustfoil and XFOIL at every alpha where both converge. For |δ| ≤ 10° the errors are in the fourth decimal place — effectively zero. For |δ| > 15° near stall, differences up to ΔCL ≈ 0.2 appear. These are real and I don’t paper over them. The two solvers are solving slightly different geometries at the hinge, so at the edge of validity they give slightly different answers.
The “Better Than XFOIL” Caveat
Let’s be honest: the convergence table looks flattering. At −20° I get 17/17 where XFOIL gets 14/17. At −15° I get 17/17 where XFOIL gets 6/17. It’d be easy to claim I’ve “beaten” XFOIL. I want to be careful about that.
For moderate deflections — anything up to about ±10° — the story is clean. Where both solvers converge, the CL values agree to the fourth decimal place and CD to the fifth. The extra converged points are at the edges of the alpha range where XFOIL’s Newton solver gives up a couple of iterations early. That’s a genuine win: same answer, slightly wider convergence envelope.
For larger deflections, it gets murkier. At ±15° and beyond, the CL differences at individual alpha points can reach 0.2. CD differences can reach 0.01. These aren’t rounding errors. They’re real, and they happen because the two solvers are solving slightly different problems.
The geometry difference is in the gap side of the hinge. XFOIL fills it with a circular arc; I leave it open and let the repaneller bridge it. For small deflections the arc is tiny and the difference is negligible. For large deflections the arc spans a meaningful fraction of the local panel spacing, and the two solvers are genuinely looking at different surface shapes near the hinge. Near stall — where the boundary layer is on the edge of separation anyway — that’s enough to tip the balance.
The honest summary: matching or exceeding XFOIL’s convergence at every tested deflection, with excellent accuracy where both converge for moderate deflections, and a known geometry approximation that introduces measurable differences at extreme conditions. I could close the gap by implementing arc-fill with a segmented spline that supports derivative discontinuities — XFOIL’s SEGSPL — but that’s a separate piece of work.
Why This Actually Matters
I have spent months doing on-and-off analyses with XFOIL on aerofoils with hinge deflections. Control surface optimisation, decambering effects, figuring out what a flap schedule looks like across an operating envelope. This is bread-and-butter work if you are designing anything with moveable surfaces — which is essentially every aeroplane.
The XFOIL workflow for this is painful. Load a foil, GDES, deflect one flap angle, go back, repanel, OPER, set up viscous, run an alpha sweep, save the polar, then do the whole dance again for the next deflection angle. For a matrix sweep — say, seven deflection angles times twenty alpha points times four Reynolds numbers — that’s 560 operating points and you’re either scripting the session file by hand or clicking through the terminal for an hour. I’ve done this. Many times. It’s not fun.
With the GDES port working, the web UI now has a one-click matrix sweep. You pick any two parameters — alpha and flap deflection, alpha and Reynolds number, flap deflection and hinge position, whatever combination you need — set the ranges, hit Generate Sweep, and it runs the full matrix. Every combination. The results land in a local database and show up immediately in the Plot Builder, where you can colour by deflection, group by Reynolds number, overlay L/D curves, and actually see what the design space looks like.
The Python API is the same story. A flap study that would have been an afternoon of XFOIL terminal scripting is now five lines:
foil = flexfoil.naca("2412")
for defl in [-10, -5, 0, 5, 10, 15, 20]:
f = foil.with_flap(hinge_x=0.75, deflection=defl)
polar = f.polar(alpha=(-4, 14, 1), Re=1e6)
print(f"δ={defl:+d}° CL_max={max(polar.cl):.3f}")That runs 119 viscous solves in about two seconds. Each one goes through the full XFOIL FLAP procedure, repanels, runs the coupled inviscid/viscous solver, and stores the result. You can sweep flap deflection against alpha, or against hinge position, or against Reynolds number. You can do it in a Jupyter notebook or in the browser. You can click on any point in the resulting polar plot and see the pressure distribution and flow field for that specific condition.
None of this is possible if the flap geometry breaks the solver. That 0/17 convergence wasn’t a theoretical problem — it was the thing standing between “XFOIL in a browser” and “XFOIL in a browser that you can actually use for control surface design.”
The Fix Is Always Understand What XFOIL Is Doing
Same lesson as Part I. The solver bugs were fixed by reading the Fortran more carefully. The geometry failure was fixed by reading the FLAP subroutine, understanding what it does and why, and recreating it in Rust. Not translating it — the line-by-line translation is the easy part, and it’s what the AI agents are good at. Understanding why SSS exists, why a simple “split at hinge x” fails, what problem the arc-fill solves — that’s what actually fixes things.
SSS is 150 lines of Newton iteration on a spline to find two arc-length break points. Doesn’t look like much. But when I tried to replace it with something simpler and watched the solver diverge on every flap case, I understood why it exists. The break points aren’t at hinge x. They’re at the precise arc-length positions where the geometry can be cleanly split, rotated, and reassembled without introducing oscillations that the BL solver can’t tolerate.
There’s a version of this story where I spent a week being clever with corner-detecting splines and spike-ratio thresholds and dense buffers, and none of it worked reliably. There’s a much shorter version where I read the Fortran, understood it, recreated it, and it worked. I did both, in that order.
When in doubt, read the Fortran. When not in doubt, also read the Fortran.