From 5b60042880b0084cadff72c83ca4f951de066690 Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 24 Aug 2026 15:02:45 +0200 Subject: [PATCH] Reconstruct NeuralFoil's surface pressure from edge velocity in arc length NeuralFoil samples at the cell centres of a uniform chordwise grid, so it reports nothing over the first and last 1/2N of chord -- the nose, which is almost all of the axial force, and the trailing edge. The old reconstruction extrapolated Cp linearly off the end of each surface, independently. Through a stagnation point the inviscid surface speed is linear in arc length and Cp = 1 - ue^2 is quadratic, so extrapolating Cp runs a straight line through the curve's own turning point. On an SK100 mid-span section at 5 deg that put Cp = -1.60 at the leading edge, where it has to approach +1, and left the two surfaces at 0.057 and 0.117 at the trailing edge where a sharp edge carries one pressure. Chord fraction is also the wrong coordinate there: the unsampled nose is 4.1% of arc against 1.6% of chord. Interpolate ue instead, in arc length, over both surfaces at once. Signing the lower surface negative makes them one continuous curve whose zero is the stagnation point, so the nose is interpolated between the innermost station on each side rather than extrapolated off the end of one, and Cp passes through exactly 1 at stagnation without that being imposed -- including the fact that stagnation sits on the lower surface at positive incidence. The stagnation point goes in as its own knot, placed where the two innermost stations interpolate linearly to zero, which is the stagnation-point-flow result rather than a fit. The trailing edge takes the mean of the two surfaces' extrapolated speeds, the Kutta condition. Between stations the interpolation is a shape-preserving monotone cubic, the family NeuralFoil's own training pipeline resamples XFoil boundary layers with. Measured on that section, integrating the surface traction: max Cp goes from 0.689 to 1.000, leading-edge Cp from -1.60 to +0.77 (XFoil 0.41), the two trailing-edge values from 0.057/0.117 to a matched 0.093, and the integrated drag from -0.0199, the wrong sign, to +0.0269 against NeuralFoil's own reported 0.0276. The lift integral loses about a point of accuracy (Cl_int 0.934 -> 0.919 against a reported 0.949), and on smooth sections whose true drag is small the integral now overshoots rather than inverting. Both are far smaller than what they replace. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 18 +++ docs/src/private_functions.md | 4 + .../airfoil_solvers/neuralfoil_solver.jl | 110 ++++++++++++++++-- src/airfoil_aero/neuralfoil.jl | 10 +- 4 files changed, 128 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 038f8be6..e0536e5a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,23 @@ # Changelog +## Unreleased + +### Fixed +- The NeuralFoil surface pressure is now reconstructed from the predicted edge + velocities in arc length, not by interpolating `Cp` in chord fraction. + NeuralFoil samples at the cell centres of a uniform grid, so it reports nothing + over the first and last `1/2N` of chord; the old code extrapolated `Cp` off the + end of each surface independently. Through a stagnation point `ue` is linear in + arc length while `Cp` is quadratic, so that extrapolated a parabola through its + own turning point: on an SK100 section it put `Cp = -1.60` at the leading edge, + where it has to approach `+1`, and the two surfaces disagreed at the trailing + edge in violation of the Kutta condition. Signing the lower surface negative + makes both surfaces one continuous curve through stagnation, so the nose is + interpolated rather than extrapolated, and the trailing edge takes one speed for + both surfaces. Integrating that section's surface traction now recovers 98% of + NeuralFoil's own reported drag where it previously came out with the wrong sign. + `neuralfoil_section` returns `ue_upper`/`ue_lower` alongside the pressures. + ## VortexStepMethod v4.2.0 2026-08-25 ### Added diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index b48ed3e8..a623fed2 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -175,6 +175,10 @@ write_aero_matrix write_node_table flat_plate_cf neuralfoil_contour_solution +contour_arc +arc_at_chord +trailing_edge_speed +velocity_knots fill_node_nans! ``` diff --git a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl index 9a625600..ae3f35a4 100644 --- a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl +++ b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl @@ -24,9 +24,10 @@ end analyze_sweep(solver::NeuralFoilSolver, def, alpha_range, Re) -> Vector{SectionSolution} Evaluate all angles (radians) in one vectorized NeuralFoil call on the deformed -Kulfan parameters, then assemble a single per-node `cp` on the deformed contour nodes -(`def.x`, `def.y`) from NeuralFoil's separate upper/lower surface pressures. `cf` is a -flat-plate closure ([`flat_plate_cf`](@ref)), NeuralFoil not exposing skin friction. +Kulfan parameters, then assemble a per-node `cp` on the deformed contour nodes +(`def.x`, `def.y`) from the predicted edge velocities +([`neuralfoil_contour_solution`](@ref)). `cf` is a flat-plate closure +([`flat_plate_cf`](@ref)), NeuralFoil not exposing skin friction. """ function analyze_sweep(solver::NeuralFoilSolver, def::DeformedSection, alpha_range, Re) res = neuralfoil_section(def.kulfan, rad2deg.(collect(alpha_range)), Re; @@ -34,22 +35,109 @@ function analyze_sweep(solver::NeuralFoilSolver, def::DeformedSection, alpha_ran n_crit=solver.n_crit, xtr_upper=solver.xtr_upper, xtr_lower=solver.xtr_lower) x, y = collect(float.(def.x)), collect(float.(def.y)) le = argmin(x) + arc = contour_arc(x, y, le) cf = [flat_plate_cf(clamp(xk, 0.0, 1.0), Re) for xk in x] - return [neuralfoil_contour_solution(alpha_range[i], res, i, x, y, le, cf) + return [neuralfoil_contour_solution(alpha_range[i], res, i, x, y, le, arc, cf) for i in eachindex(alpha_range)] end """ - neuralfoil_contour_solution(alpha, res, i, x, y, le, cf) -> SectionSolution + contour_arc(x, y, le) -> Vector + +Signed arc length of every contour node from the leading edge node `le`: positive +along the upper surface toward the trailing edge, negative along the lower. The +natural coordinate across a blunt nose, where a small chordwise step is a long one +along the skin. +""" +function contour_arc(x, y, le) + arc = zeros(Float64, length(x)) + for k in (le - 1):-1:1 + arc[k] = arc[k + 1] + hypot(x[k] - x[k + 1], y[k] - y[k + 1]) + end + for k in (le + 1):length(x) + arc[k] = arc[k - 1] - hypot(x[k] - x[k - 1], y[k] - y[k - 1]) + end + return arc +end + +""" + arc_at_chord(x, arc, indices, fractions) -> Vector + +Signed arc length at each chord fraction, read off the contour nodes `indices` — the +one surface the fractions belong to. +""" +function arc_at_chord(x, arc, indices, fractions) + nodes = collect(indices) + order = sortperm(x[nodes]) + xs, as = x[nodes][order], arc[nodes][order] + return [begin + j = clamp(searchsortedfirst(xs, f), 2, length(xs)) + gap = max(xs[j] - xs[j - 1], eps()) + as[j - 1] + (f - xs[j - 1]) / gap * (as[j] - as[j - 1]) + end for f in fractions] +end + +""" + trailing_edge_speed(res, i, upper_arc, lower_arc, upper_te, lower_te) -> Float64 + +One edge speed for both surfaces at the trailing edge: each surface's last two +stations extrapolated to its own trailing-edge arc length, the two magnitudes +averaged, so a sharp edge carries a single pressure (Kutta). +""" +function trailing_edge_speed(res, i, upper_arc, lower_arc, upper_te, lower_te) + reach(a1, a2, v1, v2, target) = + v2 + (target - a2) / (a2 - a1 + eps()) * (v2 - v1) + n = length(upper_arc) + n < 2 && return abs(res.ue_upper[end, i]) + upper = reach(upper_arc[n - 1], upper_arc[n], + res.ue_upper[n - 1, i], res.ue_upper[n, i], upper_te) + lower = reach(lower_arc[n - 1], lower_arc[n], + res.ue_lower[n - 1, i], res.ue_lower[n, i], lower_te) + return (abs(upper) + abs(lower)) / 2 +end + +""" + velocity_knots(res, i, x, arc, le) -> (knots, values) + +The whole contour's edge velocity as one curve of (signed arc length, signed +`ue/u∞`), strictly increasing in arc length. + +NeuralFoil reports nothing over the first and last `1/2N` of chord. Signing the lower +surface negative joins both surfaces into one curve through the stagnation point, so +the unsampled nose is interpolated between the innermost station on each side rather +than extrapolated off the end of one. The stagnation point enters as its own knot, +where `ue` interpolated linearly between those two stations reaches zero — the +stagnation-point-flow result, `ue` being linear in arc length there. +""" +function velocity_knots(res, i, x, arc, le) + upper_arc = arc_at_chord(x, arc, 1:le, res.x) + lower_arc = arc_at_chord(x, arc, le:length(x), res.x) + speed = trailing_edge_speed(res, i, upper_arc, lower_arc, arc[1], arc[end]) + lower_first, upper_first = abs(res.ue_lower[1, i]), abs(res.ue_upper[1, i]) + stagnation = lower_arc[1] + lower_first / max(lower_first + upper_first, eps()) * + (upper_arc[1] - lower_arc[1]) + knots = [arc[end]; reverse(lower_arc); stagnation; upper_arc; arc[1]] + values = [-speed; reverse(-abs.(res.ue_lower[:, i])); 0.0; + abs.(res.ue_upper[:, i]); speed] + order = sortperm(knots) + knots, values = knots[order], values[order] + keep = [1; [k for k in 2:length(knots) if knots[k] > knots[k - 1] + 1e-9]] + return knots[keep], values[keep] +end + +""" + neuralfoil_contour_solution(alpha, res, i, x, y, le, arc, cf) -> SectionSolution Assemble a full-contour [`SectionSolution`](@ref) for case `i` of a NeuralFoil sweep -`res`: interpolate `res.cp_upper`/`res.cp_lower` (at `res.x`) onto the contour nodes -`(x, y)` split at the leading edge `le`, carrying the precomputed `cf`. +`res`: build the one-curve edge velocity ([`velocity_knots`](@ref)), interpolate it +with a shape-preserving monotone cubic in arc length, and square it into +`Cp = 1 - ue²` at every contour node. `ue` is the interpolated quantity because it is +linear in arc length through a stagnation point, where `Cp` is quadratic. """ -function neuralfoil_contour_solution(alpha, res, i, x, y, le, cf) - up = linear_interpolation(res.x, res.cp_upper[:, i]; extrapolation_bc=Line()) - lo = linear_interpolation(res.x, res.cp_lower[:, i]; extrapolation_bc=Line()) - cp = [(k <= le ? up : lo)(clamp(x[k], 0.0, 1.0)) for k in eachindex(x)] +function neuralfoil_contour_solution(alpha, res, i, x, y, le, arc, cf) + knots, values = velocity_knots(res, i, x, arc, le) + speed = interpolate(knots, values, FritschButlandMonotonicInterpolation()) + cp = [1 - speed(clamp(a, first(knots), last(knots)))^2 for a in arc] return SectionSolution(alpha, res.cl[i], res.cd[i], res.cm[i], res.confidence[i], x, y, cp, cf) end diff --git a/src/airfoil_aero/neuralfoil.jl b/src/airfoil_aero/neuralfoil.jl index aea9d062..08a2c7ec 100644 --- a/src/airfoil_aero/neuralfoil.jl +++ b/src/airfoil_aero/neuralfoil.jl @@ -365,8 +365,10 @@ Full NeuralFoil evaluation returning integrated coefficients and the surface pressure distribution reconstructed from the predicted edge-velocity ratios (`Cp = 1 - (ue/vinf)^2`) at NeuralFoil's `N` fixed station x/c. -Returns `(; alpha, cl, cd, cm, confidence, x, cp_upper, cp_lower)`, with `x` of -length `N` and `cp_upper`/`cp_lower` sized `N × n_alpha`. +Returns `(; alpha, cl, cd, cm, confidence, x, cp_upper, cp_lower, ue_upper, +ue_lower)`, with `x` of length `N` and the four matrices sized `N × n_alpha`. A +reconstruction should interpolate `ue` and square afterwards: it is linear in arc +length through a stagnation point, where `Cp` is quadratic. """ function neuralfoil_section(params::KulfanParameters, alpha, Re; model_size::String="large", weights_dir=nothing, @@ -384,7 +386,9 @@ function neuralfoil_section(params::KulfanParameters, alpha, Re; confidence = Vector{Float64}(sigmoid.(y[1, :])), x = compute_optimal_x_points(N), cp_upper = Matrix{Float64}(1 .- upper_ue .^ 2), - cp_lower = Matrix{Float64}(1 .- lower_ue .^ 2)) + cp_lower = Matrix{Float64}(1 .- lower_ue .^ 2), + ue_upper = Matrix{Float64}(upper_ue), + ue_lower = Matrix{Float64}(lower_ue)) end """