The default solve now gives the same results as the Python VSM v2.3.0

CHANGELOG.md

@@ -15,9 +15,24 @@
   is not drawn, and still names the saved file and window.
 - The panel plots (`plot!(ax, panel)`, `plot!(ax, body)`, `plot(panel)`, `plot(body)`)
   take `border_color`, default `:black`, for the panel outlines.
+- `relaxation_factor_limit` gives the largest stable relaxation factor of the `LOOP` solver,
+  and `VSMSolution.f_attached_trailed` the attached trailed part of each panel's force.
 
 ### Changed
 
+- Default results now agree with the Python VSM v2.3.0, most visibly on swept and curved
+  wings: the circulation is `0.5 |V_eff| c Cl`; the panel chord and chordwise axis are
+  square to the bound vortex; `correct_aoa` and `is_with_attached_trailed_force` default to
+  `true`; the flow turning the forces includes a curved bound vortex correction; drag acts
+  along that flow; each wake ring follows its own panel's inflow; trailing vortex cores grow
+  with the distance from their line; `q_ref` and the reference directions leave out body
+  rates; attached trailed forces are shared per section boundary; `cd` beyond a
+  `POLAR_VECTORS` table is held, not extrapolated. V3 kite at 8°: `cl` 1.051 → 1.099, `cd`
+  0.491 → 0.458.
+- `max_iterations` defaults to 5000 and `relaxation_factor` to `nothing`, meaning 0.8 of
+  `relaxation_factor_limit`, which is lower than the old 0.03 on narrow or high-slope panels.
+- `*_unrefined_dist` coefficients and angles are area-weighted means of their panels.
+- `velocity_3D_trailing_vortex_semiinfinite!` returns zero on its own axis, not `NaN`.
 - BREAKING: `solve` and `calculate_results` are removed; use `solve!`, which returns the
   solver's `VSMSolution`, and `solve!(solver, body_aero, nothing)` to start from a fresh
   circulation as `solve` did. Keys that became other fields: `Fx`/`Fy`/`Fz` → `force`,

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                # false turns the forces by the flow at 3/4c instead of 1/4c
 
     # --- 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                # false turns the forces by the flow at 3/4c instead of 1/4c
 
   # --- Physical Properties ---
   density: 1.225 # [kg/m³] Air density (ISA sea level)
@@ -79,4 +79,3 @@ solver_settings:
   # --- Output Control ---
   calc_only_f_and_gamma: false # true=compute only forces & circulation
 
-  correct_aoa: false

docs/src/private_functions.md

@@ -33,9 +33,15 @@ nose_down
 bisect_sign_change
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
+quarter_chord_AIC!
+curved_bound_correction!
+unrolled_ring_difference
+unit_bound_velocity
+unit_trailing_velocity
+thin_airfoil_strips
+relaxation_factor_limit
 induced_velocity_at
-attached_trailed_loads
-panel_force_moment
+attached_trailed_loads!
 calculate_stall_angle_list
 wing_span_flip
 spanwise_extent
@@ -64,6 +70,8 @@ dynamic_pressure
 flow_curvature_cm
 spanwise_flow_drag
 panel_force_directions
+panel_chord
+perpendicular_part
 prescribed_va_directions
 analysis_fields!
 find_center_of_pressure!

src/filament.jl

@@ -62,7 +62,8 @@ end
                                  XVP, gamma, va, work_vectors)
 
 Calculate induced velocity by a trailing vortex filament, with a Lamb–Oseen core
-radius grown over the axial distance of `XVP` from the filament start.
+radius grown over the distance of `XVP` from the filament's line. A point on that line
+gets no velocity.
 
 # Arguments
 - `XVP`: Control point coordinates
@@ -82,9 +83,15 @@ as implemented in KiteAeroDyn".
     work_vectors
 )
     r1 = work_vectors[1]
+    r1Xr0 = work_vectors[4]
     r1 .= XVP .- filament.x1
-    axial_distance = abs(dot3(r1, filament.r0)) / filament.length
-    epsilon = sqrt(4 * ALPHA0 * NU * axial_distance / va)
+    cross3!(r1Xr0, r1, filament.r0)
+    distance = norm3(r1Xr0) / filament.length
+    if distance <= 1e-10 * max(norm3(r1), filament.length)
+        vel .= 0.0
+        return nothing
+    end
+    epsilon = sqrt(4 * ALPHA0 * NU * distance / va)
     velocity_3D_vortex_segment!(vel, filament, XVP, gamma, epsilon, work_vectors)
 end
 
@@ -189,10 +196,13 @@ function reinit!(filament::SemiInfiniteFilament{T}, x1::AbstractVector,
 end
 
 """
-    velocity_3D_trailing_vortex_semiinfinite(filament::SemiInfiniteFilament, 
-                                             Vf, XVP, GAMMA, va, work_vectors)
+    velocity_3D_trailing_vortex_semiinfinite!(vel, filament::SemiInfiniteFilament,
+                                              Vf, XVP, GAMMA, va, work_vectors)
 
-Calculate induced velocity by a semi-infinite trailing vortex filament.
+Calculate induced velocity by a semi-infinite trailing vortex filament at `XVP`, for an
+inflow speed `va` > 0 [m/s]. Inside the core radius, grown over the distance from the
+filament's line, the velocity is that on the core boundary, with no upper limit as
+`XVP` nears the line. A point on the line gets no velocity.
 """
 function velocity_3D_trailing_vortex_semiinfinite!(
     vel,
@@ -208,22 +218,24 @@ function velocity_3D_trailing_vortex_semiinfinite!(
     GAMMA = -GAMMA * filament.filament_direction
     r1 .= XVP .- filament.x1
 
-    # Core radius, grown with the axial distance of `XVP` along `Vf`.
     d_r1_Vf = dot3(r1, Vf)
     nVf = norm3(Vf)
-    epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_Vf) * nVf / va)
-
     cross3!(r1XVf, r1, Vf)
-
     nr1XVf = norm3(r1XVf)
     nr1 = norm3(r1)
-    if nr1XVf / nVf > epsilon
+    distance = nr1XVf / nVf
+    if distance < 1e-16
+        vel .= 0.0
+        return nothing
+    end
+    # Core radius, grown with the distance of `XVP` from the filament's line.
+    epsilon = sqrt(4 * ALPHA0 * NU * distance / va)
+
+    if distance > epsilon
         K = GAMMA / (4π) / (nr1XVf^2) * (1 + d_r1_Vf / nr1)
         @inbounds for k in 1:3
             vel[k] = K * r1XVf[k]
         end
-    elseif nr1XVf / nVf < 1e-12 * epsilon
-        vel .= 0.0
     else
         r1_proj = work_vectors[4]
         cross_tmp = work_vectors[5]

src/panel.jl

@@ -201,7 +201,7 @@ function build_interps(section_1::Section, section_2::Section, remove_nan)
             cm = Vector{Float64}((aero_1[4] .+ aero_2[4]) ./ 2)
             knots = polar_knots(alphas)
             cl_i = linear_interpolation(knots, cl; extrapolation_bc=extrap_flat)
-            cd_i = linear_interpolation(knots, cd; extrapolation_bc=extrap_line)
+            cd_i = linear_interpolation(knots, cd; extrapolation_bc=extrap_flat)
             cm_i = linear_interpolation(knots, cm; extrapolation_bc=extrap_flat)
         else
             (aero_1 isa Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} &&

src/panel_aerodynamics.jl

@@ -19,8 +19,21 @@ length the span width.
 @inline panel_span_vector(le_1, te_1, le_2, te_2) =
     (0.75 .* le_1 .+ 0.25 .* te_1) .- (0.75 .* le_2 .+ 0.25 .* te_2)
 
-@inline panel_chord(le_1, te_1, le_2, te_2) =
-    0.5 * (smooth_norm(te_1 .- le_1) + smooth_norm(te_2 .- le_2))
+"""
+    panel_chord(le_1, te_1, le_2, te_2)
+
+Mean length of the two section chords, each with its component along the quarter-chord
+line removed.
+"""
+@inline function panel_chord(le_1, te_1, le_2, te_2)
+    span_vec = panel_span_vector(le_1, te_1, le_2, te_2)
+    span_unit = span_vec ./ smooth_norm(span_vec)
+    return 0.5 * (smooth_norm(perpendicular_part(te_1 .- le_1, span_unit)) +
+                  smooth_norm(perpendicular_part(te_2 .- le_2, span_unit)))
+end
+
+"""Part of `vec` perpendicular to the unit vector `unit`."""
+@inline perpendicular_part(vec, unit) = vec .- dot(vec, unit) .* unit
 
 """
     panel_chord_weight(width_prev, width_own, width_next)
@@ -47,9 +60,9 @@ end
     panel_axes(le_1, te_1, le_2, te_2, chord_weight=0.5, orient=1)
 
 Airfoil frame and size of the panel between two sections, as
-`(; x_airf, y_airf, z_airf, chord, width)`. `x_airf` is chordwise, `y_airf` runs
-along the quarter-chord line the bound vortex sits on — not square to `x_airf` on
-a swept panel — and `z_airf` is `x_airf × y_airf` normalised. `chord_weight`
+`(; x_airf, y_airf, z_airf, chord, width)`, an orthonormal frame: `y_airf` runs along
+the quarter-chord line the bound vortex sits on, `x_airf` is chordwise and square to
+it, and `z_airf` is `x_airf × y_airf`. `chord_weight`
 ([`panel_chord_weight`](@ref)) is section 1's share of the chord-direction blend;
 `orient` is `±1` and flips `y_airf` and `z_airf` together.
 """
@@ -59,8 +72,10 @@ a swept panel — and `z_airf` is `x_airf × y_airf` normalised. `chord_weight`
         lean .* ((te_1 .- te_2) .- (le_1 .- le_2))
     span_vec = panel_span_vector(le_1, te_1, le_2, te_2)
     width = smooth_norm(span_vec)
-    x_airf = chord_vec ./ smooth_norm(chord_vec)
-    y_airf = orient .* (span_vec ./ width)
+    span_unit = span_vec ./ width
+    chord_square = perpendicular_part(chord_vec, span_unit)
+    x_airf = chord_square ./ smooth_norm(chord_square)
+    y_airf = orient .* span_unit
     z_cross = cross(x_airf, span_vec)
     z_airf = orient .* (z_cross ./ smooth_norm(z_cross))
     return (; x_airf, y_airf, z_airf,
@@ -149,18 +164,18 @@ Panel dynamic pressure from its two section densities and the
     0.25 * (rho_1 + rho_2) * v_rel^2
 
 """
-    panel_force_directions(axes, alpha_dir, spanwise)
+    panel_force_directions(axes, alpha_dir)
 
-Lift and drag unit vectors, `(; dir_lift, dir_drag)`. `alpha_dir` is the angle that
-turns the force, `spanwise` the wing's spanwise direction.
+Lift and drag unit vectors, `(; dir_lift, dir_drag)`. `alpha_dir` is the angle of the
+flow that turns the force: drag runs along it in the section plane and lift square to
+it and to `y_airf`.
 """
-@inline function panel_force_directions(axes, alpha_dir, spanwise)
+@inline function panel_force_directions(axes, alpha_dir)
     (; 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)
-    return (; dir_lift, dir_drag=drag_cross ./ smooth_norm(drag_cross))
+    return (; dir_lift=lift_cross ./ smooth_norm(lift_cross),
+            dir_drag=dir_inflow ./ smooth_norm(dir_inflow))
 end
 
 """Panel pitching moment per unit span about `y_airf`, positive nose-up."""

src/settings.jl

@@ -153,12 +153,12 @@ Solver configuration, used within [`VSMSettings`](@ref).
     [`LLT`](@ref Model) (default `VSM`)
 - `solver_type`: `"LOOP"` or `"NONLIN"` (default `"LOOP"`)
 - `density`: Air density (kg/m^3) (default `1.225`)
-- `max_iterations`: Maximum solver iterations (default `1500`)
+- `max_iterations`: Maximum solver iterations (default `5000`)
 - `rtol`: Relative tolerance on the fixed-point residual (default `1e-5`)
 - `tol_reference_error`: Reference error tolerance
     (default `0.001`)
-- `relaxation_factor`: Convergence relaxation factor
-    (default `0.03`)
+- `relaxation_factor`: Convergence relaxation factor, `nothing` for 0.8 of
+    [`relaxation_factor_limit`](@ref) (default `nothing`)
 - `is_with_artificial_viscosity`: Enable Li/Gaunaa post-stall
     artificial viscosity (default `false`)
 - `artificial_viscosity_factor`: Viscosity scaling coefficient k
@@ -174,35 +174,35 @@ 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 the panel forces by the flow at the quarter chord rather than at the
+    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
     force to each section (default `false`)
 - `is_with_attached_trailed_force`: Add the force on the chordwise trailed vortex
-    segments bound to each panel (default `false`)
+    segments bound to each panel (default `true`)
 """
 @with_kw mutable struct SolverSettings
     n_panels::Int64 = 40
     aerodynamic_model_type::Model = VSM
     solver_type::String = "LOOP"    # type of solver
     density::Float64 = 1.225                # air density  [kg/m³] 
-    max_iterations::Int64 = 1500
+    max_iterations::Int64 = 5000
     rtol::Float64 = 1e-5                    # relative residual tolerance [-]
     tol_reference_error::Float64 = 0.001
-    relaxation_factor::Float64 = 0.03       # relaxation factor for convergence
+    relaxation_factor::Union{Nothing, Float64} = nothing # relaxation factor for convergence
     is_with_artificial_viscosity::Bool = false  # Li/Gaunaa post-stall artificial viscosity
     artificial_viscosity_factor::Float64 = 0.035 # viscosity scaling coefficient k
     type_initial_gamma_distribution::InitialGammaDistribution = ZEROS # see: [`InitialGammaDistribution`](@ref)
     use_gamma_prev::Bool = true             # if false, always reinitialize gamma from type_initial_gamma_distribution
     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                  # turn the forces by the quarter-chord 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)
+    is_with_attached_trailed_force::Bool=true # chordwise bound-vortex force (Gaunaa 2026)
 end
 
 """