From 1c49ba235196bbbb20c212c58f24224dd766fdf5 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 23 Sep 2026 10:57:20 +0200 Subject: [PATCH 1/2] Blend POLY cl and cd into the flat-plate post-stall form and bound cm Past |alpha| = 20 deg a POLY section swapped its polynomial for the flat-plate form without blending, so cl dropped by 0.85 and cd by 0.10 across one step, and the post-stall forms had the wrong parity (cl > 0 and cd < 0 for alpha < -20 deg). cl and cd now blend over 20-25 deg with a cubic smoothstep into 2 sin|sin| cos and 2 |sin|^3, and cm holds its +-20 deg value. Inside +-20 deg nothing changes. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 9 +++++++++ docs/src/private_functions.md | 1 + docs/src/private_types.md | 2 ++ src/VortexStepMethod.jl | 4 +++- src/panel.jl | 34 ++++++++++++++++++++++--------- test/panel/test_panel.jl | 38 +++++++++++++++++++++++++++++++++++ 6 files changed, 77 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4578e6ad..6fc81498 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # Changelog +## Unreleased + +### Fixed + +- A `POLY` section's cl and cd no longer jump at |α| = 20°: they blend from the + α-polynomial into the flat-plate post-stall form over 20°–25°. Past 20°, cl now has + the sign of α and cd stays positive, and cm holds its ±20° value instead of running + its quadratic unbounded. Results inside ±20° are unchanged. + ## VortexStepMethod v6.0.0 2026-09-23 ### Added diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index f7482533..e34ebf04 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -124,6 +124,7 @@ reads_from same_knots polar_knots window_alpha +blend_post_stall assemble_polar_matrix load_polar_data load_matrix_polar_data diff --git a/docs/src/private_types.md b/docs/src/private_types.md index 0941c427..3b6d810a 100644 --- a/docs/src/private_types.md +++ b/docs/src/private_types.md @@ -19,6 +19,8 @@ SemiInfiniteFilament ```@docs LEI_AIRFOIL_BREUKELS SCAN_KNOT_MAX +POLY_FIT_LIMIT +POLY_BLEND_WIDTH ``` ### Airfoil aerodynamics (AirfoilAero) diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 6f473f1c..8fcf6075 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -254,7 +254,9 @@ Enumeration of the implemented aerodynamic models. See also: [`AeroData`](@ref) # Elements - `POLY`: α-polynomial coefficients for cl/cd/cm (e.g. Breukels LEI coeffs, generated - by the `AirfoilAero` package). Core only evaluates the polynomial. + by the `AirfoilAero` package), fitted inside ±20°. Past that, cl and cd blend into the + flat-plate forms 2 sin α |sin α| cos α and 2 |sin α|³ by ±25°, and cm holds its ±20° + value. - `POLAR_VECTORS`: Polar vectors as function of alpha (lookup tables with interpolation). A panel's table may be rewritten at run time by [`set_polar!`](@ref VortexStepMethod.set_polar!), which is how a live polar source regenerates it from the diff --git a/src/panel.jl b/src/panel.jl index 3f0b413a..5895fcb6 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -488,6 +488,25 @@ function calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity) return flow.alpha, flow.v_eff end +"Edge of the `POLY` α-polynomial fit [deg]." +const POLY_FIT_LIMIT = 20.0 +"Width over which a `POLY` cl and cd blend into the flat-plate post-stall form [deg]." +const POLY_BLEND_WIDTH = 5.0 + +""" + blend_post_stall(coeffs, flat_plate, alpha) + +`POLY` coefficient at `alpha` [rad]: the α-polynomial `coeffs` (α in degrees) inside +±`POLY_FIT_LIMIT`, the flat-plate value `flat_plate` past ±(`POLY_FIT_LIMIT` + +`POLY_BLEND_WIDTH`), and a cubic smoothstep between the two. +""" +function blend_post_stall(coeffs, flat_plate, alpha) + alpha_deg = rad2deg(alpha) + progress = clamp((abs(alpha_deg) - POLY_FIT_LIMIT) / POLY_BLEND_WIDTH, 0, 1) + weight = progress^2 * (3 - 2progress) + return (1 - weight) * evalpoly(alpha_deg, coeffs) + weight * flat_plate +end + """ calculate_cl(panel::Panel, alpha) calculate_cl(panel::Panel, alpha, delta) @@ -506,11 +525,8 @@ function calculate_cl(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} isnan(alpha) && return R(NaN) alpha = window_alpha(panel, alpha) if panel.aero_model == POLY - cl = evalpoly(rad2deg(alpha), panel.cl_coeffs) - if abs(alpha) > (π/9) - cl = 2 * cos(alpha) * sin(alpha)^2 - end - return R(cl) + flat_plate = 2 * cos(alpha) * sin(alpha) * abs(sin(alpha)) + return R(blend_post_stall(panel.cl_coeffs, flat_plate, alpha)) elseif panel.aero_model == INVISCID return R(2π * alpha) end @@ -536,10 +552,7 @@ function calculate_cd(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} isnan(alpha) && return R(NaN) alpha = window_alpha(panel, alpha) if panel.aero_model == POLY - if abs(alpha) > (π/9) # Outside ±20 degrees - return R(2 * sin(alpha)^3) - end - return R(evalpoly(rad2deg(alpha), panel.cd_coeffs)) + return R(blend_post_stall(panel.cd_coeffs, 2 * abs(sin(alpha))^3, alpha)) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cd_interp = panel.cd_interp cd_interp === nothing && @@ -566,7 +579,8 @@ function calculate_cm(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} isnan(alpha) && return R(NaN) alpha = window_alpha(panel, alpha) if panel.aero_model == POLY - return R(evalpoly(rad2deg(alpha), panel.cm_coeffs)) + return R(evalpoly(clamp(rad2deg(alpha), -POLY_FIT_LIMIT, POLY_FIT_LIMIT), + panel.cm_coeffs)) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cm_interp = panel.cm_interp cm_interp === nothing && diff --git a/test/panel/test_panel.jl b/test/panel/test_panel.jl index 63fc8231..b3acc4b6 100644 --- a/test/panel/test_panel.jl +++ b/test/panel/test_panel.jl @@ -1,4 +1,7 @@ using VortexStepMethod: Panel, Section, calculate_relative_alpha_and_relative_velocity, calculate_cl, calculate_cd_cm, reinit!, INVISCID, POLAR_VECTORS, MVec3 +using VortexStepMethod: Wing, add_section!, refine!, BodyAerodynamics, POLY, LINEAR, + calculate_cd, calculate_cm +using VortexStepMethod.AirfoilAero: lei_poly_coeffs using VortexStepMethod: panel_axes, panel_span_vector using Interpolations: linear_interpolation, Line using LinearAlgebra @@ -216,3 +219,38 @@ end @test isapprox(flipped.x_airf, axes.x_airf; atol=1e-12) end end + +@testset "POLY polar is continuous past its ±20° fit range" begin + coeffs = lei_poly_coeffs(0.1, 0.08) + wing = Wing(2; spanwise_distribution=LINEAR) + add_section!(wing, [0.0, 1.0, 0.0], [1.0, 1.0, 0.0], POLY, coeffs) + add_section!(wing, [0.0, -1.0, 0.0], [1.0, -1.0, 0.0], POLY, coeffs) + refine!(wing) + panel = BodyAerodynamics([wing]).panels[1] + polar(alpha) = (calculate_cl(panel, alpha), calculate_cd(panel, alpha), + calculate_cm(panel, alpha)) + + @testset "the polynomial is untouched inside ±20°" begin + for alpha in deg2rad.(-20:0.5:20) + @test polar(alpha) == evalpoly.(rad2deg(alpha), coeffs) + end + end + @testset "no coefficient jumps between samples 0.01° apart" begin + alphas = deg2rad.(-60:0.01:60) + values = reduce(hcat, collect.(polar.(alphas))) + @test maximum(abs, diff(values; dims=2)) < 0.01 + end + @testset "post-stall lift is odd and drag even and positive" begin + for alpha in deg2rad.(25:5:85) + @test calculate_cl(panel, -alpha) ≈ -calculate_cl(panel, alpha) + @test calculate_cd(panel, -alpha) ≈ calculate_cd(panel, alpha) + @test calculate_cd(panel, alpha) > 0 + end + end + @testset "cm holds its ±20° value past the fit range" begin + for alpha_deg in 20:5:90, side in (-1, 1) + @test calculate_cm(panel, deg2rad(side * alpha_deg)) ≈ + evalpoly(side * 20.0, coeffs[3]) + end + end +end From 758806d631fa9e56bad27f39bf13e0bf6b6f45d5 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 23 Sep 2026 11:07:09 +0200 Subject: [PATCH 2/2] Skip the POLY flat-plate form inside the fit range and test the blend is C1 blend_post_stall now returns the polynomial early inside +-20 deg, so the flat-plate form is only evaluated past it; calculate_cl at 5 deg is back to main's time. The panel test uses create_panel instead of a refined Wing, and checks the second difference of cl and cd over -60..60 deg so a C0 blend with a slope kink fails. Co-Authored-By: Claude Opus 5.5 --- docs/src/private_functions.md | 2 ++ src/panel.jl | 18 ++++++++++++------ test/panel/test_panel.jl | 18 ++++++++++-------- 3 files changed, 24 insertions(+), 14 deletions(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index e34ebf04..5e889e49 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -125,6 +125,8 @@ same_knots polar_knots window_alpha blend_post_stall +flat_plate_cl +flat_plate_cd assemble_polar_matrix load_polar_data load_matrix_polar_data diff --git a/src/panel.jl b/src/panel.jl index 5895fcb6..e90794cf 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -493,18 +493,25 @@ const POLY_FIT_LIMIT = 20.0 "Width over which a `POLY` cl and cd blend into the flat-plate post-stall form [deg]." const POLY_BLEND_WIDTH = 5.0 +"Flat-plate post-stall lift coefficient at `alpha` [rad]." +flat_plate_cl(alpha) = 2 * cos(alpha) * sin(alpha) * abs(sin(alpha)) +"Flat-plate post-stall drag coefficient at `alpha` [rad]." +flat_plate_cd(alpha) = 2 * abs(sin(alpha))^3 + """ blend_post_stall(coeffs, flat_plate, alpha) `POLY` coefficient at `alpha` [rad]: the α-polynomial `coeffs` (α in degrees) inside -±`POLY_FIT_LIMIT`, the flat-plate value `flat_plate` past ±(`POLY_FIT_LIMIT` + +±`POLY_FIT_LIMIT`, the flat-plate form `flat_plate(alpha)` past ±(`POLY_FIT_LIMIT` + `POLY_BLEND_WIDTH`), and a cubic smoothstep between the two. """ function blend_post_stall(coeffs, flat_plate, alpha) alpha_deg = rad2deg(alpha) - progress = clamp((abs(alpha_deg) - POLY_FIT_LIMIT) / POLY_BLEND_WIDTH, 0, 1) + polynomial = evalpoly(alpha_deg, coeffs) + abs(alpha_deg) <= POLY_FIT_LIMIT && return polynomial + progress = min((abs(alpha_deg) - POLY_FIT_LIMIT) / POLY_BLEND_WIDTH, 1) weight = progress^2 * (3 - 2progress) - return (1 - weight) * evalpoly(alpha_deg, coeffs) + weight * flat_plate + return (1 - weight) * polynomial + weight * flat_plate(alpha) end """ @@ -525,8 +532,7 @@ function calculate_cl(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} isnan(alpha) && return R(NaN) alpha = window_alpha(panel, alpha) if panel.aero_model == POLY - flat_plate = 2 * cos(alpha) * sin(alpha) * abs(sin(alpha)) - return R(blend_post_stall(panel.cl_coeffs, flat_plate, alpha)) + return R(blend_post_stall(panel.cl_coeffs, flat_plate_cl, alpha)) elseif panel.aero_model == INVISCID return R(2π * alpha) end @@ -552,7 +558,7 @@ function calculate_cd(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} isnan(alpha) && return R(NaN) alpha = window_alpha(panel, alpha) if panel.aero_model == POLY - return R(blend_post_stall(panel.cd_coeffs, 2 * abs(sin(alpha))^3, alpha)) + return R(blend_post_stall(panel.cd_coeffs, flat_plate_cd, alpha)) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cd_interp = panel.cd_interp cd_interp === nothing && diff --git a/test/panel/test_panel.jl b/test/panel/test_panel.jl index b3acc4b6..3f5c237d 100644 --- a/test/panel/test_panel.jl +++ b/test/panel/test_panel.jl @@ -1,6 +1,5 @@ using VortexStepMethod: Panel, Section, calculate_relative_alpha_and_relative_velocity, calculate_cl, calculate_cd_cm, reinit!, INVISCID, POLAR_VECTORS, MVec3 -using VortexStepMethod: Wing, add_section!, refine!, BodyAerodynamics, POLY, LINEAR, - calculate_cd, calculate_cm +using VortexStepMethod: POLY, calculate_cd, calculate_cm using VortexStepMethod.AirfoilAero: lei_poly_coeffs using VortexStepMethod: panel_axes, panel_span_vector using Interpolations: linear_interpolation, Line @@ -222,13 +221,11 @@ end @testset "POLY polar is continuous past its ±20° fit range" begin coeffs = lei_poly_coeffs(0.1, 0.08) - wing = Wing(2; spanwise_distribution=LINEAR) - add_section!(wing, [0.0, 1.0, 0.0], [1.0, 1.0, 0.0], POLY, coeffs) - add_section!(wing, [0.0, -1.0, 0.0], [1.0, -1.0, 0.0], POLY, coeffs) - refine!(wing) - panel = BodyAerodynamics([wing]).panels[1] + panel = create_panel(Section([0.0, 1.0, 0.0], [1.0, 1.0, 0.0], POLY, coeffs), + Section([0.0, -1.0, 0.0], [1.0, -1.0, 0.0], POLY, coeffs)) polar(alpha) = (calculate_cl(panel, alpha), calculate_cd(panel, alpha), calculate_cm(panel, alpha)) + alphas = deg2rad.(-60:0.01:60) @testset "the polynomial is untouched inside ±20°" begin for alpha in deg2rad.(-20:0.5:20) @@ -236,10 +233,15 @@ end end end @testset "no coefficient jumps between samples 0.01° apart" begin - alphas = deg2rad.(-60:0.01:60) values = reduce(hcat, collect.(polar.(alphas))) @test maximum(abs, diff(values; dims=2)) < 0.01 end + @testset "cl and cd slopes do not kink at the blend edges" begin + for coefficient in (calculate_cl, calculate_cd) + values = coefficient.(Ref(panel), alphas) + @test maximum(abs, diff(diff(values))) < 1e-4 + end + end @testset "post-stall lift is odd and drag even and positive" begin for alpha in deg2rad.(25:5:85) @test calculate_cl(panel, -alpha) ≈ -calculate_cl(panel, alpha)