Measure alpha in the plane square to the bound vortex (Gaunaa CP1)

CHANGELOG.md

@@ -8,6 +8,10 @@
   the right tip, z = x × y up. The separate `KB` definition, whose `Z = Y × X` pointed
   down, is gone; the solver's frame is unchanged. CL, CD and CS are documented as
   wind-axis coefficients, `cfx`, `cfy` and `cfz` as the body-axis ones.
+- The angle of attack is measured in the plane square to the quarter-chord line, from
+  the chord's projection on that plane (Gaunaa, Li & Pirrung 2026, CP1), everywhere it
+  is formed and where it turns the lift direction. Swept panels read their polars at a
+  different angle than before; unswept panels are unchanged.
 
 ## VortexStepMethod v6.0.0 2026-09-23
 

docs/src/private_functions.md

@@ -31,7 +31,6 @@ set_pitch_rate_dist!
 coeffs_at_angles
 nose_down
 bisect_sign_change
-calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
 calculate_stall_angle_list
@@ -41,6 +40,7 @@ calculate_circulation_distribution_elliptical_wing
 _compute_reference_velocity_from_distribution
 smooth_distribution!
 make_dual_shadow
+matrix_row
 ```
 
 ### Panel aerodynamics
@@ -54,6 +54,7 @@ smooth_norm
 panel_span_vector
 panel_chord_weight
 panel_axes
+inflow_angle
 effective_alpha
 panel_inflow
 dynamic_pressure

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() for _ in 1:15]
+- `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache() for _ in 1:12]
 """
 @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() for _ in 1:15]
+    cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:12]
 end
 
 """
@@ -523,20 +523,16 @@ end
 
 """
     update_effective_angle_of_attack!(alpha_corrected, body_aero::BodyAerodynamics, gamma,
-                                      core_radius_fraction, z_airf_dist, x_airf_dist,
-                                      va_vec_dist, va_dist, va_unit_dist)
+                                      core_radius_fraction, va_vec_dist, va_dist,
+                                      va_unit_dist)
 
-Update angle of attack at aerodynamic center for VSM method.
-
-Returns:
-    nothing
+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.
 """
 function update_effective_angle_of_attack!(alpha_corrected,
-    body_aero::BodyAerodynamics, 
+    body_aero::BodyAerodynamics,
     gamma,
     core_radius_fraction,
-    z_airf_dist,
-    x_airf_dist,
     va_vec_dist,
     va_dist,
     va_unit_dist)
@@ -551,31 +547,11 @@ function update_effective_angle_of_attack!(alpha_corrected,
         mul!(view(induced_velocity, :, k), view(body_aero.AIC_aero_center, :, :, k), gamma)
     end
 
-    # In-place relative velocity calculation
-    relative_velocity = body_aero.cache[2][va_vec_dist]
-    relative_velocity .= va_vec_dist .+ induced_velocity
-
-    # Preallocate and compute dot products manually
-    n = size(relative_velocity, 1)
-    v_normal     = body_aero.cache[3][relative_velocity]
-    v_tangential = body_aero.cache[4][relative_velocity]
-
-    @inbounds for i in 1:n
-        vn = 0.0
-        vt = 0.0
-        for j in 1:3
-            vn += z_airf_dist[i, j] * relative_velocity[i, j]
-            vt += x_airf_dist[i, j] * relative_velocity[i, j]
-        end
-        v_normal[i] = vn
-        v_tangential[i] = vt
-    end
-
-    # Direct angle calculation without temporary arrays
-    @inbounds for i in 1:n
-        alpha_corrected[i] = atan(v_normal[i], v_tangential[i])
+    @inbounds for (i, panel) in enumerate(body_aero.panels)
+        v_rel = matrix_row(va_vec_dist, i) .+ matrix_row(induced_velocity, i)
+        alpha_corrected[i] = inflow_angle(v_rel, SVector{3}(panel.y_airf),
+                                          SVector{3}(panel.z_airf))
     end
-
     nothing
 end
 
@@ -832,21 +808,21 @@ function calculate_results(
 )
 
     n_panels = length(panels)
-    if length(body_aero.cache) < 15
-        append!(body_aero.cache, [LazyBufferCache() for _ in 1:(15 - length(body_aero.cache))])
+    if length(body_aero.cache) < 12
+        append!(body_aero.cache, [LazyBufferCache() for _ in 1:(12 - length(body_aero.cache))])
     end
 
-    cl_dist = body_aero.cache[5][alpha_dist]
-    cd_dist = body_aero.cache[6][alpha_dist]
-    cm_dist = body_aero.cache[7][alpha_dist]
-    panel_width_dist = body_aero.cache[8][alpha_dist]
-    alpha_corrected = body_aero.cache[9][alpha_dist]
-    cl_prescribed_va = body_aero.cache[10][alpha_dist]
-    cd_prescribed_va = body_aero.cache[11][alpha_dist]
-    cs_prescribed_va = body_aero.cache[12][alpha_dist]
-    f_body_3D = body_aero.cache[13][alpha_dist, (3, length(alpha_dist))]
-    m_body_3D = body_aero.cache[14][alpha_dist, (3, length(alpha_dist))]
-    alpha_geometric = body_aero.cache[15][alpha_dist]
+    cl_dist = body_aero.cache[2][alpha_dist]
+    cd_dist = body_aero.cache[3][alpha_dist]
+    cm_dist = body_aero.cache[4][alpha_dist]
+    panel_width_dist = body_aero.cache[5][alpha_dist]
+    alpha_corrected = body_aero.cache[6][alpha_dist]
+    cl_prescribed_va = body_aero.cache[7][alpha_dist]
+    cd_prescribed_va = body_aero.cache[8][alpha_dist]
+    cs_prescribed_va = body_aero.cache[9][alpha_dist]
+    f_body_3D = body_aero.cache[10][alpha_dist, (3, length(alpha_dist))]
+    m_body_3D = body_aero.cache[11][alpha_dist, (3, length(alpha_dist))]
+    alpha_geometric = body_aero.cache[12][alpha_dist]
 
     fill!(f_body_3D, 0.0)
     fill!(m_body_3D, 0.0)
@@ -861,19 +837,8 @@ function calculate_results(
                 body_aero.pitch_rate_dist[i], chord_dist[i], v_rel_dist[i])
         end
         panel_width_dist[i] = panel.width
-        va = va_dist[i]
-        x_norm = norm3(panel.x_airf)
-        z_norm = norm3(panel.z_airf)
-        if va == 0.0 || x_norm == 0.0 || z_norm == 0.0
-            alpha_geometric[i] = NaN
-        else
-            inv_va = 1.0 / va
-            v_tangential = -dot3(panel.x_airf, panel.va_vec) *
-                           inv_va / x_norm
-            v_normal = -dot3(panel.z_airf, panel.va_vec) *
-                       inv_va / z_norm
-            alpha_geometric[i] = atan(-v_normal, -v_tangential)
-        end
+        alpha_geometric[i] = iszero(va_dist[i]) ? NaN : inflow_angle(panel.va_vec,
+            SVector{3}(panel.y_airf), SVector{3}(panel.z_airf))
     end
 
     # Calculate alpha corrections based on model type
@@ -883,8 +848,6 @@ function calculate_results(
             body_aero,
             gamma_new,
             core_radius_fraction,
-            z_airf_dist,
-            x_airf_dist,
             va_vec_dist,
             va_dist,
             va_unit_dist

src/panel.jl

@@ -479,16 +479,6 @@ function calculate_relative_alpha_and_relative_velocity(
     return flow.alpha, flow.v_eff
 end
 
-"""
-    calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity)
-
-Calculate relative angle of attack and relative velocity of the panel.
-"""
-function calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity)
-    flow = panel_inflow(panel_axes(panel), panel.va_vec, panel.va_vec, induced_velocity)
-    return flow.alpha, flow.v_eff
-end
-
 """
     calculate_cl(panel::Panel, alpha)
     calculate_cl(panel::Panel, alpha, delta)

src/panel_aerodynamics.jl

@@ -95,6 +95,18 @@ already sampled at three-quarter chord.
     ifelse(v_rel > 0, -0.25π * pitch_rate * chord / (2 * smooth_norm(v_rel)),
            zero(pitch_rate * chord))
 
+"""
+    inflow_angle(velocity, y_airf, z_airf)
+
+Angle of attack of `velocity` in the plane square to `y_airf`: from the chord's
+projection on that plane, `y_airf × z_airf`, towards `z_airf`.
+"""
+@inline inflow_angle(velocity, y_airf, z_airf) =
+    atan(dot(velocity, z_airf), dot(velocity, cross(y_airf, z_airf)))
+
+"""Row `i` of an `n×3` matrix as a 3-vector."""
+@inline matrix_row(matrix, i) = SVector(matrix[i, 1], matrix[i, 2], matrix[i, 3])
+
 """
     effective_alpha(alpha, deficiency)
 
@@ -115,9 +127,9 @@ leaves it zero. `deficiency` feeds [`effective_alpha`](@ref).
 """
 @inline function panel_inflow(axes, va_vec_1, va_vec_2, v_ind, dva_vec_1=nothing,
                               dva_vec_2=nothing, deficiency=0)
-    (; x_airf, y_airf, z_airf, chord) = axes
+    (; y_airf, z_airf, chord) = axes
     v_eff = 0.5 .* (va_vec_1 .+ va_vec_2) .+ v_ind
-    alpha = atan(dot(v_eff, z_airf), dot(v_eff, x_airf))
+    alpha = inflow_angle(v_eff, y_airf, z_airf)
     v_span = cross(v_eff, y_airf)
     pitch_rate = isnothing(dva_vec_1) ? zero(chord) :
         section_pitch_rate(0.5 .* (dva_vec_1 .+ dva_vec_2), z_airf, chord)
@@ -143,8 +155,8 @@ Lift and drag unit vectors, `(; dir_lift, dir_drag)`. `alpha_dir` is the angle t
 turns the force, `spanwise` the wing's spanwise direction.
 """
 @inline function panel_force_directions(axes, alpha_dir, spanwise)
-    (; x_airf, y_airf, z_airf) = axes
-    dir_inflow = cos(alpha_dir) .* x_airf .+ sin(alpha_dir) .* z_airf
+    (; y_airf, z_airf) = axes
+    dir_inflow = cos(alpha_dir) .* cross(y_airf, z_airf) .+ sin(alpha_dir) .* z_airf
     lift_cross = cross(dir_inflow, y_airf)
     dir_lift = lift_cross ./ smooth_norm(lift_cross)
     drag_cross = cross(spanwise, dir_lift)

src/solver.jl

@@ -178,7 +178,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref)
     # Intermediate results
     lr::LoopResult{P, T} = LoopResult{P, T}()
     br::BaseResult{P, T} = BaseResult{P, T}()
-    cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:10]
+    cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:8]
     cache_base::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}  = [LazyBufferCache()]
     cache_lin::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:4]
 
@@ -373,35 +373,9 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
         end
         width_dist[i] = panel.width
 
-        # Geometric AoA using panel-local axes and prescribed
-        # freestream — scalar ops to avoid allocations
-        begin
-            va1 = solver.sol.va_vec_dist[i,1]
-            va2 = solver.sol.va_vec_dist[i,2]
-            va3 = solver.sol.va_vec_dist[i,3]
-            va = sqrt(va1^2 + va2^2 + va3^2)
-            x1 = solver.sol._x_airf_dist[i,1]
-            x2 = solver.sol._x_airf_dist[i,2]
-            x3 = solver.sol._x_airf_dist[i,3]
-            x_norm = sqrt(x1^2 + x2^2 + x3^2)
-            z1 = solver.sol._z_airf_dist[i,1]
-            z2 = solver.sol._z_airf_dist[i,2]
-            z3 = solver.sol._z_airf_dist[i,3]
-            z_norm = sqrt(z1^2 + z2^2 + z3^2)
-            if va == 0 || x_norm == 0 || z_norm == 0
-                alpha_geometric_dist[i] = NaN
-            else
-                inv_va = -1.0 / va
-                vu1 = va1 * inv_va
-                vu2 = va2 * inv_va
-                vu3 = va3 * inv_va
-                v_tangential = (x1*vu1+x2*vu2+x3*vu3) / x_norm
-                v_normal = (z1*vu1+z2*vu2+z3*vu3) / z_norm
-                alpha_geometric_dist[i] = atan(
-                    -v_normal, -v_tangential)
-            end
-        end
-
+        va_vec = matrix_row(solver.sol.va_vec_dist, i)
+        alpha_geometric_dist[i] = iszero(va_vec) ? NaN : inflow_angle(va_vec,
+            matrix_row(solver.sol._y_airf_dist, i), matrix_row(solver.sol._z_airf_dist, i))
     end
 
     lift = solver.sol.lift_dist
@@ -415,8 +389,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
             body_aero,
             gamma_new,
             solver.core_radius_fraction,
-            solver.sol._z_airf_dist,
-            solver.sol._x_airf_dist,
             solver.sol.va_vec_dist,
             solver.br.va_dist,
             solver.br.va_unit_dist
@@ -775,9 +747,6 @@ end
     relative_velocity_crossz,
     v_acrossz_dist,
     z_airf_dist,
-    x_airf_dist,
-    v_normal_dist,
-    v_tangential_dist,
     va_magw_dist,
     cl_dist,
     chord_dist,
@@ -806,20 +775,11 @@ end
     end
 
     @inbounds for i in 1:n_panels
-        v_normal_dist[i] =
-            z_airf_dist[i,1]*relative_velocity_dist[i,1] +
-            z_airf_dist[i,2]*relative_velocity_dist[i,2] +
-            z_airf_dist[i,3]*relative_velocity_dist[i,3]
-        v_tangential_dist[i] =
-            x_airf_dist[i,1]*relative_velocity_dist[i,1] +
-            x_airf_dist[i,2]*relative_velocity_dist[i,2] +
-            x_airf_dist[i,3]*relative_velocity_dist[i,3]
-        solver.lr.v_span_dist[i] =
-            y_airf_dist[i,1]*relative_velocity_dist[i,1] +
-            y_airf_dist[i,2]*relative_velocity_dist[i,2] +
-            y_airf_dist[i,3]*relative_velocity_dist[i,3]
+        v_rel = matrix_row(relative_velocity_dist, i)
+        y_airf = matrix_row(y_airf_dist, i)
+        solver.lr.alpha_dist[i] = inflow_angle(v_rel, y_airf, matrix_row(z_airf_dist, i))
+        solver.lr.v_span_dist[i] = dot(y_airf, v_rel)
     end
-    solver.lr.alpha_dist .= atan.(v_normal_dist, v_tangential_dist)
 
     @inbounds for i in 1:n_panels
         solver.lr.v_rel_dist[i] = smooth_sqrt(
@@ -945,7 +905,6 @@ function gamma_loop!(
 ) where {P, U, T}
     va_vec_dist = solver.sol.va_vec_dist
     chord_dist = solver.sol._chord_dist
-    x_airf_dist = solver.sol._x_airf_dist
     y_airf_dist = solver.sol._y_airf_dist
     z_airf_dist = solver.sol._z_airf_dist
     solver.lr.converged   = false
@@ -961,8 +920,6 @@ function gamma_loop!(
     relative_velocity_crossz = solver.cache[6][va_vec_dist]
     v_acrossz_dist           = solver.cache[7][va_vec_dist]
     cl_dist                  = solver.cache[8][solver.lr.gamma_new]
-    v_normal_dist            = solver.cache[9][solver.lr.gamma_new]
-    v_tangential_dist        = solver.cache[10][solver.lr.gamma_new]
 
     AIC_x = @view body_aero.AIC[:, :, 1]
     AIC_y = @view body_aero.AIC[:, :, 2]
@@ -991,9 +948,7 @@ function gamma_loop!(
             velocity_view_x, velocity_view_y, velocity_view_z,
             va_vec_dist, induced_velocity_all, relative_velocity_dist,
             y_airf_dist, relative_velocity_crossz, v_acrossz_dist,
-            z_airf_dist, x_airf_dist,
-            v_normal_dist, v_tangential_dist,
-            va_magw_dist, cl_dist, chord_dist,
+            z_airf_dist, va_magw_dist, cl_dist, chord_dist,
         )
         @inbounds for i in 1:n_panels
             residual[i] -= gamma_iter[i]
@@ -1013,9 +968,7 @@ function gamma_loop!(
                     velocity_view_x, velocity_view_y, velocity_view_z,
                     va_vec_dist, induced_velocity_all, relative_velocity_dist,
                     y_airf_dist, relative_velocity_crossz, v_acrossz_dist,
-                    z_airf_dist, x_airf_dist,
-                    v_normal_dist, v_tangential_dist,
-                    va_magw_dist, cl_dist, chord_dist,
+                    z_airf_dist, va_magw_dist, cl_dist, chord_dist,
                 )
                 inv_step = 1.0 / step
                 for i in 1:n_panels
@@ -1042,9 +995,7 @@ function gamma_loop!(
                     velocity_view_x, velocity_view_y, velocity_view_z,
                     va_vec_dist, induced_velocity_all, relative_velocity_dist,
                     y_airf_dist, relative_velocity_crossz, v_acrossz_dist,
-                    z_airf_dist, x_airf_dist,
-                    v_normal_dist, v_tangential_dist,
-                    va_magw_dist, cl_dist, chord_dist,
+                    z_airf_dist, va_magw_dist, cl_dist, chord_dist,
                 )
                 @inbounds for i in 1:n_panels
                     residual_perturbed[i] -= gamma_perturbed[i]
@@ -1108,9 +1059,6 @@ function gamma_loop!(
                 relative_velocity_crossz,
                 v_acrossz_dist,
                 z_airf_dist,
-                x_airf_dist,
-                v_normal_dist,
-                v_tangential_dist,
                 va_magw_dist,
                 cl_dist,
                 chord_dist,