diff --git a/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md b/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md new file mode 100644 index 00000000..19470cb6 --- /dev/null +++ b/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md @@ -0,0 +1,7 @@ +### Changed + +- `correct_aoa` defaults to `true`: the polars are still read at the 3/4-chord control point, and each section force is turned by the flow at its 1/4-chord aerodynamic centre (Gaunaa, Li & Pirrung, TORQUE 2026). Every VSM result moves, the induced drag most. That flow now leaves out only the section's own bound vortex, not every panel's, which also moves `LLT` results on swept or curved wings. + +### Fixed + +- A vortex segment evaluated on its own axis returns zero rather than NaN when `core_radius_fraction` is near zero, and ForwardDiff sees the velocity's slope there rather than zero. diff --git a/data/TUDELFT_V3_KITE/vsm_settings.yaml b/data/TUDELFT_V3_KITE/vsm_settings.yaml index 5477aed9..9326b2bf 100644 --- a/data/TUDELFT_V3_KITE/vsm_settings.yaml +++ b/data/TUDELFT_V3_KITE/vsm_settings.yaml @@ -53,7 +53,7 @@ solver_settings: # --- Core Aerodynamic Model --- aerodynamic_model_type: VSM # VSM=3D vortex method, LLT=lifting line theory solver_type: LOOP # LOOP=fixed-point iteration, NEWTON=Newton-Raphson - #correct_aoa: false # Defaults to false, meaning the computed AoA at 3/4c is used + #correct_aoa: true # true=forces turned by the 1/4-chord flow # --- Physical Properties --- density: 1.225 # [kg/m³] Air density (ISA sea level) diff --git a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml index d3fb98e8..9c44a7a9 100644 --- a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml +++ b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml @@ -57,7 +57,7 @@ solver_settings: # --- Core Aerodynamic Model --- aerodynamic_model_type: VSM # VSM=3D vortex method, LLT=lifting line theory solver_type: LOOP # LOOP=fixed-point iteration, NONLIN=nonlinear solver - #correct_aoa: false # Defaults to false, meaning the computed AoA at 3/4c is used + correct_aoa: true # true=forces turned by the 1/4-chord flow # --- Physical Properties --- density: 1.225 # [kg/m³] Air density (ISA sea level) @@ -78,5 +78,3 @@ solver_settings: # --- Output Control --- calc_only_f_and_gamma: false # true=compute only forces & circulation - - correct_aoa: false diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 948d5668..b7debcd8 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -28,7 +28,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru influence coefficients, used only for the corrected angle of attack - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the KA body frame [m²] - `c_ref::Float64` = 1.0: Reference chord length (max panel chord) [m] -- `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache()] +- `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache() for _ in 1:3] """ @with_kw mutable struct BodyAerodynamics{P, W<:AbstractWing, T, PN<:Panel{T}} panels::Vector{PN} @@ -49,7 +49,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru AIC_aero_center::Array{T, 3} = zeros(T, P, P, 3) projected_area::T = one(T) c_ref::T = one(T) - cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache()] + cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:3] end """ @@ -429,7 +429,9 @@ end calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model, core_radius_fraction, va_dist, va_unit_dist, target=body_aero.AIC) -Calculate Aerodynamic Influence Coefficient matrices. +Calculate the velocity each panel's unit-strength vortex ring induces at every control +point (`VSM`, less the panel's own 2D bound vortex) or aerodynamic centre (`LLT`, less +the panel's own bound filament). See also: [`BodyAerodynamics`](@ref), [`Model`](@ref) @@ -440,22 +442,19 @@ Returns: nothing va_dist::AbstractVector{T}, va_unit_dist::AbstractMatrix{T}, target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T} - # Determine evaluation point based on model - evaluation_point = model == VSM ? :control_point : :aero_center - evaluation_point_on_bound = model == LLT - - # Allocate work vectors for this function (separate from those used by child functions) - velocity_induced = zeros(MVector{3, T}) - tempvel = zeros(MVector{3, T}) - va_unit = zeros(MVector{3, T}) - U_2D = zeros(MVector{3, T}) + # Slots 1-5 are the filament kernels' scratch. + velocity_induced = body_aero.work_vectors[6] + U_2D = body_aero.work_vectors[7] + tempvel = body_aero.work_vectors[8] + va_unit = body_aero.work_vectors[9] # Python parity: one shared area-weighted wake vector for all panels. - panel_areas = [panel.chord * panel.width for panel in body_aero.panels] - va_vec_dist = zeros(T, length(body_aero.panels), 3) - @inbounds for i in 1:length(body_aero.panels), k in 1:3 - va_vec_dist[i, k] = va_unit_dist[i, k] * va_dist[i] + panel_areas = body_aero.cache[2][va_dist] + for (i, panel) in enumerate(body_aero.panels) + panel_areas[i] = panel.chord * panel.width end + va_vec_dist = body_aero.cache[3][va_unit_dist] + va_vec_dist .= va_unit_dist .* va_dist wake_velocity = _compute_reference_velocity_from_distribution( va_vec_dist, length(body_aero.panels), @@ -472,14 +471,13 @@ Returns: nothing filaments = panel_jring.filaments for icp in eachindex(body_aero.panels) panel_icp = body_aero.panels[icp] - ep = evaluation_point == :control_point ? panel_icp.control_point : - panel_icp.aero_center + ep = model == VSM ? panel_icp.control_point : panel_icp.aero_center calculate_velocity_induced_single_ring_semiinfinite!( velocity_induced, tempvel, filaments, ep, - evaluation_point_on_bound, + model == LLT && icp == jring, va, va_unit, one(T), @@ -487,7 +485,6 @@ Returns: nothing body_aero.work_vectors ) - # Subtract 2D induced velocity for VSM if icp == jring && model == VSM calculate_velocity_induced_bound_2D!(U_2D, panel_jring, ep, body_aero.work_vectors) velocity_induced .-= U_2D @@ -527,7 +524,8 @@ end va_unit_dist) Write each panel's [`inflow_angle`](@ref) at its aerodynamic centre into -`alpha_corrected`, the induced velocity taken from `gamma` on the LLT influence matrix. +`alpha_corrected`, the induced velocity taken from `gamma` without the panel's own bound +filament. """ function update_effective_angle_of_attack!(alpha_corrected, body_aero::BodyAerodynamics, @@ -537,8 +535,7 @@ function update_effective_angle_of_attack!(alpha_corrected, va_dist, va_unit_dist) - # Its own buffer: `AIC` holds the control-point matrix the circulation was solved - # against, so overwriting it here would leave post-solve readers on the LLT one. + # Not `AIC`: that holds the control-point matrix post-solve readers expect. calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_dist, va_unit_dist, body_aero.AIC_aero_center) diff --git a/src/filament.jl b/src/filament.jl index 023addec..d79f612d 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -94,7 +94,8 @@ end Calculate the Biot–Savart velocity induced by a straight vortex segment at `XVP`. Inside the core radius `epsilon` the velocity is evaluated on the core boundary and -scaled linearly with the distance to the axis. +scaled linearly with the distance to the axis. Without a core, it is zero within +1e-12 of the axis, relative to the distance from `x1`. """ @inline function velocity_3D_vortex_segment!( vel, @@ -104,15 +105,18 @@ scaled linearly with the distance to the axis. epsilon, work_vectors ) - r1, r2, r1Xr2, r1Xr0, r1r2norm, r1_proj, r2_proj = work_vectors + r1, r2, r1Xr2, r1Xr0, r1r2norm = work_vectors r0 = filament.r0 nr0 = filament.length r1 .= XVP .- filament.x1 r2 .= XVP .- filament.x2 cross3!(r1Xr0, r1, r0) - nr1Xr0 = norm3(r1Xr0) - if nr1Xr0 / nr0 > epsilon + axis_distance = norm3(r1Xr0) / nr0 + rounding = 1e-12 * norm3(r1) + if epsilon <= rounding && axis_distance <= rounding + vel .= 0.0 + elseif axis_distance > epsilon cross3!(r1Xr2, r1, r2) nr1 = norm3(r1) nr2 = norm3(r2) @@ -124,35 +128,15 @@ scaled linearly with the distance to the axis. @inbounds for k in 1:3 vel[k] = coeff * r1Xr2[k] end - elseif nr1Xr0 / nr0 < 1e-12 * epsilon - vel .= 0.0 else nr0sq = nr0 * nr0 d_r1_r0 = dot3(r1, r0) d_r2_r0 = dot3(r2, r0) - r_rad = r1Xr0 + d_sum = d_r1_r0 / sqrt(d_r1_r0^2 / nr0sq + epsilon^2) - + d_r2_r0 / sqrt(d_r2_r0^2 / nr0sq + epsilon^2) + coeff = -(gamma / (4π)) * d_sum / (epsilon^2 * nr0sq) @inbounds for k in 1:3 - r_rad[k] = r1[k] - d_r1_r0 * r0[k] / nr0sq - end - nr_rad = norm3(r_rad) - @inbounds for k in 1:3 - r1_proj[k] = d_r1_r0 * r0[k] / nr0sq + - epsilon * r_rad[k] / nr_rad - r2_proj[k] = d_r2_r0 * r0[k] / nr0sq + - epsilon * r_rad[k] / nr_rad - end - cross3!(r1Xr2, r1_proj, r2_proj) - nr1_proj = norm3(r1_proj) - nr2_proj = norm3(r2_proj) - d_sum = 0.0 - @inbounds for k in 1:3 - d_sum += r0[k] * (r1_proj[k]/nr1_proj - - r2_proj[k]/nr2_proj) - end - scale = nr1Xr0 / (nr0 * epsilon) - coeff = scale * (gamma / (4π)) / (norm3(r1Xr2)^2) * d_sum - @inbounds for k in 1:3 - vel[k] = coeff * r1Xr2[k] + vel[k] = coeff * r1Xr0[k] end end nothing diff --git a/src/settings.jl b/src/settings.jl index 4062d56d..8e34492a 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -174,8 +174,8 @@ Solver configuration, used within [`VSMSettings`](@ref). - `calc_only_f_and_gamma`: Leave the analysis fields of [`VSMSolution`](@ref) (`lift`, `cl`, their distributions, span and centers of pressure) at their last values (default `false`) -- `correct_aoa`: Perform angle of attack correction (VSM model only) - (default `false`) +- `correct_aoa`: Turn each section force by the flow at its aerodynamic centre + rather than at its control point, VSM model only (default `true`) - `flow_curvature`: Add the thin-airfoil pitch-rate moment increment to each section (default `false`) - `is_with_viscous_drag_correction`: Add the spanwise-flow viscous drag and side @@ -199,7 +199,7 @@ Solver configuration, used within [`VSMSettings`](@ref). core_radius_fraction::Float64 = 0.05 mu::Float64 = 1.81e-5 # dynamic viscosity [N·s/m²] calc_only_f_and_gamma::Bool=false # skip the analysis fields of VSMSolution - correct_aoa::Bool=false # perform aoa correction + correct_aoa::Bool=true # force directions from the aero-centre flow flow_curvature::Bool=false # thin-airfoil pitch-rate moment increment is_with_viscous_drag_correction::Bool=false # spanwise-flow viscous force (Gaunaa 2024) is_with_attached_trailed_force::Bool=false # chordwise bound-vortex force (Gaunaa 2026) diff --git a/src/solver.jl b/src/solver.jl index a8011d04..f91e327f 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -171,6 +171,8 @@ Main solver structure for the Vortex Step Method. See also: [`solve!`](@ref) - mu::Float64 = 1.81e-5: Dynamic viscosity [N·s/m²] - `is_only_f_and_gamma_output`::Bool = false: Whether `solve!` skips the fields of [`VSMSolution`](@ref) that only analysis reads, see [`SolverSettings`](@ref) +- `correct_aoa`::Bool = true: Turn each section force by the flow at its aerodynamic + centre rather than at its control point, VSM model only - `flow_curvature`::Bool = false: Add the thin-airfoil pitch-rate moment increment `-(π/4) q̂` to each section, see: [`flow_curvature_cm`](@ref) - `is_with_viscous_drag_correction`::Bool = false: Add the spanwise-flow viscous drag @@ -210,7 +212,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref) core_radius_fraction::T = T(0.05) mu::T = T(1.81e-5) is_only_f_and_gamma_output::Bool = false - correct_aoa::Bool = false + correct_aoa::Bool = true flow_curvature::Bool = false is_with_viscous_drag_correction::Bool = false is_with_attached_trailed_force::Bool = false @@ -404,7 +406,7 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; aerodynamic_model_type = solver.aerodynamic_model_type # Calculate coefficients for each panel - for (i, panel) in enumerate(panels) # zero bytes + for (i, panel) in enumerate(panels) cl_dist[i] = calculate_cl(panel, alpha_dist[i]) cd_dist[i], cm_dist[i] = calculate_cd_cm(panel, alpha_dist[i]) if solver.flow_curvature @@ -423,8 +425,7 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; drag = solver.sol.drag_dist panel_moment_dist = solver.sol.panel_moment_dist - # Calculate alpha corrections based on model type - if solver.correct_aoa && aerodynamic_model_type == VSM # 64 bytes + if solver.correct_aoa && aerodynamic_model_type == VSM update_effective_angle_of_attack!( alpha_corrected, body_aero, diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index a47650a6..2cd8d749 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/test/body_aerodynamics/test_body_aerodynamics.jl @@ -1,5 +1,7 @@ using VortexStepMethod -using VortexStepMethod: calculate_cl, calculate_cd_cm, calculate_projected_area, calculate_AIC_matrices! +using VortexStepMethod: calculate_cl, calculate_cd_cm, calculate_projected_area, + calculate_AIC_matrices!, velocity_3D_bound_vortex!, velocity_3D_trailing_vortex!, + velocity_3D_trailing_vortex_semiinfinite! using LinearAlgebra using Test using Logging @@ -24,6 +26,18 @@ function inviscid_wing(section_y; n_panels=length(section_y) - 1) return wing end +""" + uniform_AIC!(body_aero, model, core_radius_fraction, va_vec) + +Write into `body_aero.AIC` the influence matrices of `model` under the uniform inflow +`va_vec` [m/s]. +""" +function uniform_AIC!(body_aero, model, core_radius_fraction, va_vec) + n_panels = length(body_aero.panels) + calculate_AIC_matrices!(body_aero, model, core_radius_fraction, + fill(norm(va_vec), n_panels), repeat(normalize(va_vec)', n_panels)) +end + @testset "Induction Matrix Creation" begin # Setup n_panels = 3 @@ -77,17 +91,7 @@ end LLT ) - # Calculate new matrices - va_dist = fill(norm(va_vec), length(body_aero.panels)) - va_unit_dist = repeat(reshape(va_vec ./ norm(va_vec), 1, 3), - length(body_aero.panels)) - calculate_AIC_matrices!( - body_aero, - LLT, - core_radius_fraction, - va_dist, - va_unit_dist - ) + uniform_AIC!(body_aero, LLT, core_radius_fraction, va_vec) AIC_x, AIC_y, AIC_z = @views body_aero.AIC[:, :, 1], body_aero.AIC[:, :, 2], body_aero.AIC[:, :, 3] # Compare matrices @@ -113,17 +117,7 @@ end VSM ) - # Calculate new matrices - va_dist = fill(norm(va_vec), length(body_aero.panels)) - va_unit_dist = repeat(reshape(va_vec ./ norm(va_vec), 1, 3), - length(body_aero.panels)) - calculate_AIC_matrices!( - body_aero, - VSM, - core_radius_fraction, - va_dist, - va_unit_dist - ) + uniform_AIC!(body_aero, VSM, core_radius_fraction, va_vec) AIC_x, AIC_y, AIC_z = body_aero.AIC[:, :, 1], body_aero.AIC[:, :, 2], body_aero.AIC[:, :, 3] # Compare matrices with higher precision for VSM @@ -133,6 +127,64 @@ end end end +@testset "LLT matrix leaves out only each panel's own bound filament" begin + radius = 4.0 + wing = Wing(3) + for theta in range(-π / 4, π / 4, length=4) + section_point = [0.0, radius * sin(theta), radius * cos(theta)] + add_section!(wing, section_point, section_point .+ [1.0, 0.0, 0.0], INVISCID) + end + refine!(wing) + body_aero = BodyAerodynamics([wing]) + va_vec = 10.0 .* [cosd(5), 0.0, sind(5)] + set_va!(body_aero, va_vec) + core_radius_fraction = 0.05 + uniform_AIC!(body_aero, LLT, core_radius_fraction, va_vec) + + work_vectors = body_aero.work_vectors + velocity = zeros(3) + panels = body_aero.panels + for (i, panel_i) in enumerate(panels), (j, panel_j) in enumerate(panels) + filaments = panel_j.filaments + expected = zeros(3) + for filament in filaments[2:3] + velocity_3D_trailing_vortex!(velocity, filament, panel_i.aero_center, 1.0, + norm(va_vec), work_vectors) + expected .+= velocity + end + for filament in filaments[4:5] + velocity_3D_trailing_vortex_semiinfinite!(velocity, filament, normalize(va_vec), + panel_i.aero_center, 1.0, norm(va_vec), work_vectors) + expected .+= velocity + end + if i != j + velocity_3D_bound_vortex!(velocity, filaments[1], panel_i.aero_center, 1.0, + core_radius_fraction, work_vectors) + @test norm(velocity) > 1e-3 + expected .+= velocity + end + @test body_aero.AIC[i, j, :] ≈ expected atol = 1e-12 + end +end + +@testset "LLT matrix is finite on a swept wing without vortex core" begin + wing = Wing(8) + for y in range(-2.0, 2.0, length=9) + leading_edge = [0.2abs(y), y, -0.5abs(y)] + add_section!(wing, leading_edge, leading_edge .+ [1.0, 0.0, 0.0], INVISCID) + end + refine!(wing) + body_aero = BodyAerodynamics([wing]) + va_vec = 10.0 .* [cosd(5), 0.0, sind(5)] + set_va!(body_aero, va_vec) + uniform_AIC!(body_aero, LLT, 1e-10, va_vec) + AIC_small_core = copy(body_aero.AIC) + uniform_AIC!(body_aero, LLT, 1e-20, va_vec) + + @test all(isfinite, body_aero.AIC) + @test body_aero.AIC ≈ AIC_small_core +end + @testset "Wing Geometry Creation" begin @testset "Origin Translation" begin diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl index 397eaf99..e786d247 100644 --- a/test/filament/test_bound_filament.jl +++ b/test/filament/test_bound_filament.jl @@ -1,5 +1,6 @@ using VortexStepMethod: BoundFilament, velocity_3D_bound_vortex!, velocity_3D_trailing_vortex!, reinit!, ALPHA0, NU +using ForwardDiff using LinearAlgebra using Test @@ -26,6 +27,22 @@ function analytical_solution(control_point, gamma) dot(r0, r1/norm(r1) - r2/norm(r2)) end +""" + off_axis_velocity(offset, gamma, core_radius_fraction) + +z velocity [m/s] induced by the unit filament from the origin along x at the point `offset` +[m] off its axis, half a length beyond its end. +""" +function off_axis_velocity(offset, gamma, core_radius_fraction) + T = typeof(offset) + filament = BoundFilament{T}() + reinit!(filament, zeros(T, 3), T[1, 0, 0]) + velocity = zeros(T, 3) + velocity_3D_bound_vortex!(velocity, filament, [1.5, offset, 0.0], gamma, + core_radius_fraction, ntuple(_ -> zeros(T, 3), 10)) + return velocity[3] +end + @testset "BoundFilament Tests" begin gamma = 1.0 core_radius_fraction = 0.01 @@ -260,4 +277,12 @@ end @test v[3] > 0 end end + + @testset "ForwardDiff sees the core's slope on the extended axis" begin + velocity_at(offset) = off_axis_velocity(offset, gamma, core_radius_fraction) + slope = ForwardDiff.derivative(velocity_at, 0.0) + + @test slope != 0 + @test slope ≈ velocity_at(1e-4) / 1e-4 rtol = 1e-9 + end end diff --git a/test/solver/test_polar_vectors.jl b/test/solver/test_polar_vectors.jl index 48439163..2544223a 100644 --- a/test/solver/test_polar_vectors.jl +++ b/test/solver/test_polar_vectors.jl @@ -38,7 +38,7 @@ end sol = solve_three_panel_wing(POLAR_VECTORS, linear_polar(4π, 0.1, -0.02), linear_polar(2π, 0.05, -0.01)) @test sol.solver_status == FEASIBLE - @test sol.cl_dist ≈ [4π, 3π, 2π] .* sol.alpha_dist + @test sol.cl_dist ≈ [4π, 3π, 2π] .* sol.alpha_uncorrected @test sol.cd_dist ≈ [0.1, 0.075, 0.05] @test sol.cm_dist ≈ [-0.02, -0.015, -0.01] @testset "the stronger +y side rolls the wing positive about x" begin diff --git a/test/solver/test_stability.jl b/test/solver/test_stability.jl index 758b68a1..180dcec8 100644 --- a/test/solver/test_stability.jl +++ b/test/solver/test_stability.jl @@ -81,7 +81,7 @@ end @testset "moments about the leading edge: stable trim" begin body_aero, solver = trimmable_wing(0.05) - trims = trim_angle(solver, body_aero, beta, va) + trims = trim_angle(solver, body_aero, beta, va; alpha_tol=1e-7) @test length(trims) == 1 trim = only(trims) trim_coeffs = coeffs_at_angles(solver, body_aero, trim.alpha, beta, va) @@ -93,7 +93,7 @@ end @testset "moments about the trailing edge: unstable trim" begin body_aero, solver = trimmable_wing(-0.05; reference_point=[1.0, 0.0, 0.0]) - trim = only(trim_angle(solver, body_aero, beta, va)) + trim = only(trim_angle(solver, body_aero, beta, va; alpha_tol=1e-7)) trim_coeffs = coeffs_at_angles(solver, body_aero, trim.alpha, beta, va) @test abs(trim_coeffs[5]) < 1e-5 @test trim.dCMy_dalpha > 0 diff --git a/test/verification/test_verification.jl b/test/verification/test_verification.jl index cdf39467..7ca63f41 100644 --- a/test/verification/test_verification.jl +++ b/test/verification/test_verification.jl @@ -49,15 +49,15 @@ function wing_from_coordinates(coordinates, aero_model, aero_data=nothing; end """ - lift_drag_polar(wing, model, alphas; va, relaxation_factor) + lift_drag_polar(wing, model, alphas; va, relaxation_factor, correct_aoa=true) Wing `CL` and `CD` at each angle of attack in `alphas` [deg], solved with the settings of the Python verification cases. """ -function lift_drag_polar(wing, model, alphas; va, relaxation_factor) +function lift_drag_polar(wing, model, alphas; va, relaxation_factor, correct_aoa=true) body_aero = BodyAerodynamics([wing]) solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=model, - relaxation_factor, core_radius_fraction=1e-20) + relaxation_factor, core_radius_fraction=1e-20, correct_aoa) CL = zeros(length(alphas)) CD = zeros(length(alphas)) for (i, alpha) in enumerate(alphas) @@ -105,6 +105,19 @@ max_error(actual, expected) = maximum(abs.(actual .- expected)) @test max_error(CD_llt, CD_vsm) < 1e-2 end + @testset "elliptic AR 4 wing reaches Oswald efficiency 1 (Gaunaa et al. 2026)" begin + aspect_ratio = 4.0 + span = aspect_ratio * π / 4 + coordinates = generate_coordinates_el_wing(1.0, span, 41, "cos") + wing = wing_from_coordinates(coordinates, INVISCID) + oswald_efficiency(CL, CD) = only(CL)^2 / (π * aspect_ratio * only(CD)) + polar(correct_aoa) = lift_drag_polar(wing, VSM, [10.0]; va=10.0, + relaxation_factor=0.05, correct_aoa) + + @test oswald_efficiency(polar(true)...) ≈ 1 atol = 0.01 + @test oswald_efficiency(polar(false)...) < 0.8 + end + @testset "curved Clark Y wing matches RANS" begin alphas = [9.0] coordinates = generate_coordinates_curved_wing(2.18, 6.969, π / 4, 4.673, 60,