Add Li/Gaunaa spanwise artificial viscosity for post-stall stabilization

data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml

@@ -37,7 +37,7 @@
 # Define the flight state for the aerodynamic analysis
 condition:
   wind_speed: 10.0 # [m/s] Free stream velocity magnitude
-  alpha: 5.0 # [°] Angle of attack (pitch angle relative to flow)
+  alpha: 0.0 # [°] Angle of attack (pitch angle relative to flow)
   beta: 0.0 # [°] Sideslip angle (yaw angle relative to flow)
   yaw_rate: 0.0 # [°/s] Yaw rate (for dynamic analysis, 0 for static)
 

docs/src/private_functions.md

@@ -8,4 +8,6 @@ calculate_AIC_matrices!
 update_panel_properties!
 calculate_inertia_tensor
 unrefined_deform!
+build_spanwise_laplacian!
+local_lift_slope!
 ```

examples/billowing.jl

@@ -2,8 +2,8 @@ using LinearAlgebra
 using VortexStepMethod
 
 PLOT = true
-SAVE_ALL = false
-USE_TEX = false
+SAVE_ALL = true
+USE_TEX = true
 OUTPUT_DIR = joinpath(dirname(@__DIR__), "output")
 
 # Data paths (all within this repo)
@@ -113,7 +113,7 @@ solver_bill = make_solver(body_aero_bill)
 
 # --- Set flight conditions ---
 wind_speed = condition_cfg["wind_speed"]
-angle_of_attack_deg = condition_cfg["alpha"]
+angle_of_attack_deg = 10.0
 sideslip_deg = condition_cfg["beta"]
 
 α0 = deg2rad(angle_of_attack_deg)
@@ -136,8 +136,8 @@ println("Billowed:  CL=$(round(results_bill["cl"]; digits=4)), " *
 if PLOT
     # Plot geometry (flat wing)
     plot_geometry(
-        body_aero_flat,
-        "Flat wing geometry";
+        body_aero_bill,
+        "Billowing wing geometry";
         save_path=OUTPUT_DIR,
         is_save=false || SAVE_ALL,
         is_show=true,

ext/VortexStepMethodControlPlotsExt.jl

@@ -426,7 +426,7 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la
     for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list)
         axs[1, 0].plot(
             y_coordinates_i,
-            result_i["alpha_geometric"],
+            rad2deg.(result_i["alpha_geometric"]),
             label=label_i
         )
     end
@@ -439,7 +439,7 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la
     for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list)
         axs[1, 1].plot(
             y_coordinates_i,
-            result_i["alpha_at_ac"],
+            rad2deg.(result_i["alpha_at_ac"]),
             label=label_i
         )
     end
@@ -452,7 +452,7 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la
     for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list)
         axs[1, 2].plot(
             y_coordinates_i,
-            result_i["alpha_uncorrected"],
+            rad2deg.(result_i["alpha_uncorrected"]),
             label=label_i
         )
     end

ext/VortexStepMethodMakieExt.jl

@@ -605,19 +605,19 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la
 
     # Plot alpha geometric
     for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list)
-        lines!(ax_alpha_geo, Vector(y_coords), Vector(results["alpha_geometric"]),
+        lines!(ax_alpha_geo, Vector(y_coords), rad2deg.(Vector(results["alpha_geometric"])),
             label=label)
     end
 
     # Plot alpha at ac
     for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list)
-        lines!(ax_alpha_ac, Vector(y_coords), Vector(results["alpha_at_ac"]),
+        lines!(ax_alpha_ac, Vector(y_coords), rad2deg.(Vector(results["alpha_at_ac"])),
             label=label)
     end
 
     # Plot alpha uncorrected
     for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list)
-        lines!(ax_alpha_unc, Vector(y_coords), Vector(results["alpha_uncorrected"]),
+        lines!(ax_alpha_unc, Vector(y_coords), rad2deg.(Vector(results["alpha_uncorrected"])),
             label=label)
     end
 

src/body_aerodynamics.jl

@@ -782,7 +782,7 @@ function calculate_results(
                            inv_va_norm / x_norm
             v_normal = -dot3(panel.z_airf, panel.va) *
                        inv_va_norm / z_norm
-            alpha_geometric[i] = pi + atan(v_normal, v_tangential)
+            alpha_geometric[i] = atan(-v_normal, -v_tangential)
         end
     end
 

src/settings.jl

@@ -70,6 +70,10 @@ Solver configuration, used within [`VSMSettings`](@ref).
 - `artificial_damping`: Enable artificial damping
     (default `false`)
 - `k2`, `k4`: Artificial damping parameters
+- `is_with_artificial_viscosity`: Enable Li/Gaunaa post-stall
+    artificial viscosity (default `false`)
+- `artificial_viscosity_factor`: Viscosity scaling coefficient k
+    (default `0.035`)
 - `type_initial_gamma_distribution`:
     [`ELLIPTIC`](@ref InitialGammaDistribution) or `ZEROS`
     (default `ELLIPTIC`)
@@ -95,6 +99,8 @@ Solver configuration, used within [`VSMSettings`](@ref).
     artificial_damping::Bool = false        # whether to apply artificial damping
     k2::Float64 = 0.1                       # artificial damping parameter
     k4::Float64 = 0.0                       # artificial damping parameter
+    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 = ELLIPTIC # see: [InitialGammaDistribution](@ref)
     use_gamma_prev::Bool = true             # if false, always reinitialize gamma from type_initial_gamma_distribution
     core_radius_fraction::Float64 = 1e-20

src/solver.jl

@@ -123,6 +123,12 @@ Main solver structure for the Vortex Step Method.See also: [solve](@ref)
 - `is_with_artificial_damping`::Bool = false: Whether to apply artificial damping
 - `artificial_damping`::NamedTuple{(:k2, :k4), Tuple{Float64, Float64}} = (k2=0.1, k4=0.0): Artificial damping parameters
 
+## Artificial viscosity settings
+- `is_with_artificial_viscosity`::Bool = false: Enable the Li/Gaunaa spanwise artificial
+    viscosity (TORQUE 2026) for post-stall stabilization in the LOOP solver
+- `artificial_viscosity_factor`::Float64 = 0.035: Coefficient k in the viscosity scaling
+    (the conservative envelope from the paper)
+
 ## Additional settings
 - `type_initial_gamma_distribution`::InitialGammaDistribution = ELLIPTIC: see: [InitialGammaDistribution](@ref)
 - `use_gamma_prev`::Bool = true: reuse provided previous gamma as initial guess when available
@@ -156,6 +162,10 @@ sol::VSMSolution = VSMSolution(): The result of calling [solve!](@ref)
     is_with_artificial_damping::Bool = false
     artificial_damping::NamedTuple{(:k2, :k4), Tuple{Float64, Float64}} =(k2=0.1, k4=0.0)
 
+    # Li/Gaunaa spanwise artificial viscosity (TORQUE 2026) settings
+    is_with_artificial_viscosity::Bool = false
+    artificial_viscosity_factor::T = T(0.035)
+
     # Additional settings
     type_initial_gamma_distribution::InitialGammaDistribution = ZEROS
     use_gamma_prev::Bool = true
@@ -196,6 +206,8 @@ function Solver(body_aero, settings::VSMSettings)
         relaxation_factor=ss.relaxation_factor,
         is_with_artificial_damping=ss.artificial_damping,
         artificial_damping=(k2=ss.k2, k4=ss.k4),
+        is_with_artificial_viscosity=ss.is_with_artificial_viscosity,
+        artificial_viscosity_factor=ss.artificial_viscosity_factor,
         type_initial_gamma_distribution=ss.type_initial_gamma_distribution,
         use_gamma_prev=ss.use_gamma_prev,
         core_radius_fraction=ss.core_radius_fraction,
@@ -302,8 +314,8 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
                 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] = pi + atan(
-                    v_normal, v_tangential)
+                alpha_geometric_dist[i] = atan(
+                    -v_normal, -v_tangential)
             end
         end
 
@@ -759,12 +771,63 @@ end
     return nothing
 end
 
+"""
+    build_spanwise_laplacian!(laplacian, n_panels)
+
+Fill the `n_panels × n_panels` matrix `laplacian` with the discrete spanwise
+Laplacian used by the Li/Gaunaa post-stall artificial-viscosity regularization.
+Interior rows use the three-point stencil `gamma[i-1] - 2 gamma[i] + gamma[i+1]`;
+the tip rows use the second-order closures of Li, Gaunaa, Pirrung & Lønbæk
+(TORQUE 2026, Eq. 15) that enforce `gamma -> 0` at the wing tips. Panels are
+assumed ordered consecutively along the span and approximately uniformly spaced.
+Returns the matrix unchanged (all zeros) when `n_panels < 3`.
+"""
+function build_spanwise_laplacian!(laplacian, n_panels)
+    laplacian .= 0
+    n_panels < 3 && return laplacian
+    @inbounds for i in 2:n_panels-1
+        laplacian[i, i-1] = 1.0
+        laplacian[i, i]   = -2.0
+        laplacian[i, i+1] = 1.0
+    end
+    laplacian[1, 1] = -4.0
+    laplacian[1, 2] = 4.0 / 3.0
+    laplacian[n_panels, n_panels]   = -4.0
+    laplacian[n_panels, n_panels-1] = 4.0 / 3.0
+    return laplacian
+end
+
+"""
+    local_lift_slope!(slopes, panels, alpha_dist, delta=deg2rad(0.5))
+
+Per-panel local lift-curve slope `dCl/dalpha`, evaluated from each panel's own
+2-D polar at its current effective angle of attack by central differences. The
+slope turns negative in post-stall, which is what activates the
+artificial-viscosity regularization in [`gamma_loop!`](@ref).
+"""
+function local_lift_slope!(slopes, panels, alpha_dist, delta=deg2rad(0.5))
+    inv_2delta = 1.0 / (2delta)
+    @inbounds for i in eachindex(panels)
+        cl_plus  = calculate_cl(panels[i], alpha_dist[i] + delta)
+        cl_minus = calculate_cl(panels[i], alpha_dist[i] - delta)
+        slopes[i] = (cl_plus - cl_minus) * inv_2delta
+    end
+    return slopes
+end
+
 """
     gamma_loop!(solver::Solver, AIC_x::Matrix{Float64},
               AIC_y::Matrix{Float64}, AIC_z::Matrix{Float64},
               panels::AbstractVector{<:Panel}, relaxation_factor::Float64; log=true)
 
 Main iteration loop for calculating circulation distribution.
+
+When `solver.is_with_artificial_viscosity` is set, the LOOP solver replaces the
+explicit target `F(gamma)` with the implicit Li/Gaunaa solution
+`(I - diag(mu) L) gamma = F(gamma)` before relaxation, stabilizing post-stall
+distributions. The Python reference gates this on precomputed per-panel stall
+angles; here the polar is an interpolation object, so the expensive linear solve
+is gated on `any(mu > 0)` instead, which is behaviourally identical.
 """
 function gamma_loop!(
     solver::Solver{P, U, T},
@@ -894,6 +957,21 @@ function gamma_loop!(
     end
 
     if solver.solver_type == LOOP
+        # Work buffers allocated once: the geometry is frozen during iteration.
+        use_viscosity = solver.is_with_artificial_viscosity
+        laplacian        = use_viscosity ? zeros(T, n_panels, n_panels) : zeros(T, 0, 0)
+        viscosity_matrix = use_viscosity ? zeros(T, n_panels, n_panels) : zeros(T, 0, 0)
+        lift_slope       = use_viscosity ? zeros(T, n_panels) : zeros(T, 0)
+        mu_array         = use_viscosity ? zeros(T, n_panels) : zeros(T, 0)
+        gamma_target     = use_viscosity ? zeros(T, n_panels) : zeros(T, 0)
+        planform_area    = zero(T)
+        if use_viscosity
+            build_spanwise_laplacian!(laplacian, n_panels)
+            @inbounds for i in 1:n_panels
+                planform_area += panels[i].width * chord_array[i]
+            end
+        end
+
         function f_loop!(gamma_new, gamma, damp)
             gamma .= gamma_new
             update_gamma_candidate!(
@@ -922,10 +1000,31 @@ function gamma_loop!(
                 cl_dist,
                 chord_array,
             )
+            # Gate the linear solve on any(mu > 0): fires exactly in post-stall.
+            if use_viscosity
+                local_lift_slope!(lift_slope, panels, solver.lr.alpha_dist)
+                any_stalled = false
+                @inbounds for i in 1:n_panels
+                    m = -solver.artificial_viscosity_factor * planform_area *
+                        lift_slope[i] / panels[i].width^2
+                    mu_array[i] = max(zero(T), m)
+                    mu_array[i] > 0 && (any_stalled = true)
+                end
+                if any_stalled
+                    gamma_target .= gamma_new
+                    @inbounds for col in 1:n_panels, row in 1:n_panels
+                        viscosity_matrix[row, col] =
+                            (row == col ? one(T) : zero(T)) -
+                            mu_array[row] * laplacian[row, col]
+                    end
+                    gamma_new .= viscosity_matrix \ gamma_target
+                end
+            end
+
             # Update gamma with relaxation and damping
-            @. gamma_new = (1 - relaxation_factor) * gamma + 
+            @. gamma_new = (1 - relaxation_factor) * gamma +
                     relaxation_factor * gamma_new + damp
-            
+
             # Apply damping if needed
             if solver.is_with_artificial_damping
                 smooth_circulation!(damp, gamma, 0.1, 0.5)
@@ -1087,6 +1186,8 @@ function make_dual_shadow(solver::Solver{P, U, Float64},
         atol = TD(solver.atol),
         is_with_artificial_damping = solver.is_with_artificial_damping,
         artificial_damping = solver.artificial_damping,
+        is_with_artificial_viscosity = solver.is_with_artificial_viscosity,
+        artificial_viscosity_factor = TD(solver.artificial_viscosity_factor),
         type_initial_gamma_distribution = solver.type_initial_gamma_distribution,
         use_gamma_prev = false,
         core_radius_fraction = TD(solver.core_radius_fraction),

test/solver/test_solver.jl

@@ -91,3 +91,60 @@ calc_forces_allocs(solver, body_aero) =
         rm(settings_file; force=true)
     end
 end
+
+@testset "Spanwise Laplacian tip closures" begin
+    # Interior three-point stencil plus the Eq. 15 tip closures.
+    n = 5
+    laplacian = zeros(n, n)
+    VortexStepMethod.build_spanwise_laplacian!(laplacian, n)
+
+    @test laplacian[3, :] == [0.0, 1.0, -2.0, 1.0, 0.0]
+    @test laplacian[1, :] == [-4.0, 4.0/3.0, 0.0, 0.0, 0.0]
+    @test laplacian[end, :] == [0.0, 0.0, 0.0, 4.0/3.0, -4.0]
+    @test sum(laplacian[3, :]) == 0.0
+
+    # Degenerate spans stay all-zero (no regularization possible).
+    small = zeros(2, 2)
+    VortexStepMethod.build_spanwise_laplacian!(small, 2)
+    @test all(small .== 0.0)
+end
+
+@testset "Artificial viscosity stabilizes post-stall" begin
+    # Paper case b: AR=10, Cl'=-pi; plain iteration diverges, implicit converges.
+    AR, clp, N, V, c = 10.0, -pi, 50, 1.0, 2.0
+    b = AR * c
+    dz = b / N
+    alpha_g = deg2rad(5.0)
+    idx = collect(0:N-1)
+    induction = [1.0 / (pi * dz) / (4.0 * (j - i)^2 - 1.0) for i in idx, j in idx]
+
+    residual(gamma) =
+        0.5 .* V .* c .* (clp .* ((alpha_g .+ (induction * gamma) ./ V) .- alpha_g) .+ 1.0)
+
+    laplacian = zeros(N, N)
+    VortexStepMethod.build_spanwise_laplacian!(laplacian, N)
+    mu = -0.035 * N^2 / AR * clp                  # Eq. (16), > 0 since clp < 0
+    system = Matrix(1.0I, N, N) .- mu .* laplacian
+
+    function iterate_loop(use_viscosity, omega; max_iter=20000)
+        gamma = zeros(N)
+        for _ in 1:max_iter
+            prev = gamma
+            nxt = use_viscosity ? (system \ residual(prev)) : residual(prev)
+            gamma = (1 - omega) .* prev .+ omega .* nxt
+            ref = max(maximum(abs.(gamma)), 1e-4)
+            maximum(abs.(gamma .- prev)) / ref < 1e-6 && return true, gamma
+        end
+        return false, gamma
+    end
+
+    base_converged, _ = iterate_loop(false, 0.1)
+    av_converged, gamma_av = iterate_loop(true, 0.5)
+
+    @test !base_converged    # base fixed point diverges (sawtooth)
+    @test av_converged
+    peak = maximum(abs.(gamma_av))
+    sawtooth = sum(abs.(gamma_av[1:end-2] .- 2 .* gamma_av[2:end-1] .+
+        gamma_av[3:end])) / (length(gamma_av) - 2) / peak
+    @test sawtooth < 0.02
+end