Default correct_aoa to true, with the aero-centre flow leaving out only the section's own bound vortex

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.

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)

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

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,22 +471,20 @@ 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),
                 core_radius_fraction,
                 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)
 

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

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)

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,

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