From 030da6168591d5dd9d521df9611397cb1e882667 Mon Sep 17 00:00:00 2001 From: Bart Date: Sat, 27 Jun 2026 10:16:45 +0200 Subject: [PATCH 1/3] Add Li/Gaunaa spanwise artificial viscosity for post-stall stabilization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Port the opt-in implicit artificial-viscosity regularization (Li, Gaunaa, Pirrung & Lønbæk, TORQUE 2026) from the Python Vortex-Step-Method to the LOOP solver. Post-stall (negative lift-slope) circulation distributions otherwise develop non-physical sawtooth oscillations and never converge; the implicit scheme (I - diag(mu) L) gamma = F(gamma) stays stable at relaxation factors of order one. Adds the discrete spanwise Laplacian with second-order tip closures (Eq. 15), per-panel lift-slope evaluation, and gating so attached flow is a no-op. The expensive linear solve is gated on any(mu > 0) rather than on precomputed stall angles, since the Julia polar is an interpolation object. Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/private_functions.md | 2 + src/settings.jl | 6 ++ src/solver.jl | 109 ++++++++++++++++++++++++++++++++-- test/solver/test_solver.jl | 57 ++++++++++++++++++ 4 files changed, 170 insertions(+), 4 deletions(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 7cebb81d..435f6277 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -8,6 +8,8 @@ CurrentModule = VortexStepMethod ```@docs calculate_AIC_matrices! gamma_loop! +build_spanwise_laplacian! +local_lift_slope! frozen_wake! calc_forces! calculate_cl diff --git a/src/settings.jl b/src/settings.jl index 062c6213..6e8429b6 100644 --- a/src/settings.jl +++ b/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 diff --git a/src/solver.jl b/src/solver.jl index 5ad73422..268d2b95 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -126,6 +126,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 @@ -159,6 +165,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 @@ -199,6 +209,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, @@ -347,8 +359,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 @@ -804,12 +816,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}, @@ -939,6 +1002,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!( @@ -967,10 +1045,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) @@ -1133,6 +1232,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), diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index d74769af..9d2c5b05 100644 --- a/test/solver/test_solver.jl +++ b/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 From fbaf1f5a26832cf83fce9dd932e42bf416de4ab6 Mon Sep 17 00:00:00 2001 From: Bart Date: Sat, 27 Jun 2026 18:31:28 +0200 Subject: [PATCH 2/3] Fix plotting --- data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml | 2 +- examples/billowing.jl | 10 +++++----- ext/VortexStepMethodControlPlotsExt.jl | 6 +++--- ext/VortexStepMethodMakieExt.jl | 6 +++--- src/body_aerodynamics.jl | 2 +- 5 files changed, 13 insertions(+), 13 deletions(-) diff --git a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml index 13d757b7..4a2308cb 100644 --- a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml +++ b/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) diff --git a/examples/billowing.jl b/examples/billowing.jl index 1f6b886f..f8fc056b 100644 --- a/examples/billowing.jl +++ b/examples/billowing.jl @@ -6,8 +6,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) @@ -117,7 +117,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) @@ -140,8 +140,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, diff --git a/ext/VortexStepMethodControlPlotsExt.jl b/ext/VortexStepMethodControlPlotsExt.jl index dba0be96..e0f5d177 100644 --- a/ext/VortexStepMethodControlPlotsExt.jl +++ b/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 diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index fce74bb1..8017851b 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -607,19 +607,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 diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 39cd3696..b384162c 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -785,7 +785,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 From 9680bf70603e949fe990b85193251fa0adf8299d Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 13 Jul 2026 13:55:30 +0200 Subject: [PATCH 3/3] Extract apply_artificial_viscosity! and test it via real solver code The previous post-stall test reimplemented the whole scheme inline and never touched the solver, so the actual artificial-viscosity path in gamma_loop! had 0% patch coverage. - Extract the viscosity step from gamma_loop! into apply_artificial_viscosity!, so it is directly testable and no longer duplicated inline. - Solve in place with ldiv!(gamma, lu!(M), rhs) instead of M \ rhs, cutting the post-stall allocation from ~3.7 kB to the LU pivot vector; the attached hot path stays zero-alloc. - Replace the reimplementing test with a unit test of apply_artificial_viscosity! on a real post-stall wing (asserts fire/no-op, smoothing, and zero allocation on the attached path) and a solve! integration test (attached no-op is bit-identical to viscosity off; high-AoA solve stays finite). The firing test regenerates the NeuralFoil polar over alpha_range -5:40 deg because cl clamps flat past the tabulated range, so the default 15 deg polar never reaches negative lift slope and the feature never activates. Co-Authored-By: Claude Opus 4.8 --- docs/src/private_functions.md | 1 + src/solver.jl | 56 ++++++++++++------ test/solver/test_solver.jl | 107 ++++++++++++++++++++++------------ 3 files changed, 111 insertions(+), 53 deletions(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 435f6277..f7c96b22 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -10,6 +10,7 @@ calculate_AIC_matrices! gamma_loop! build_spanwise_laplacian! local_lift_slope! +apply_artificial_viscosity! frozen_wake! calc_forces! calculate_cl diff --git a/src/solver.jl b/src/solver.jl index 268d2b95..0e9f0915 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -860,6 +860,40 @@ function local_lift_slope!(slopes, panels, alpha_dist, delta=deg2rad(0.5)) return slopes end +""" + apply_artificial_viscosity!(gamma, panels, alpha_dist, laplacian, viscosity_matrix, + lift_slope, mu_array, gamma_target, planform_area, factor) + +Apply one implicit Li/Gaunaa artificial-viscosity step to `gamma` in place and return +`true` when it fired. The per-panel viscosity is +`mu_i = max(0, -factor * planform_area * Cl'_i / width_i^2)`, with the local lift slope +`Cl'_i` from [`local_lift_slope!`](@ref) evaluated at `alpha_dist`. When any panel is +post-stall (`mu_i > 0`), `gamma` is replaced by the solution of +`(I - diag(mu) L) gamma = gamma`, with `L` the spanwise Laplacian in `laplacian` +(see [`build_spanwise_laplacian!`](@ref)); otherwise `gamma` is left unchanged. The +remaining arguments are preallocated work buffers reused across iterations. +""" +function apply_artificial_viscosity!(gamma, panels, alpha_dist, laplacian, viscosity_matrix, + lift_slope, mu_array, gamma_target, planform_area, factor) + n_panels = length(panels) + local_lift_slope!(lift_slope, panels, alpha_dist) + any_stalled = false + @inbounds for i in 1:n_panels + m = -factor * planform_area * lift_slope[i] / panels[i].width^2 + mu_array[i] = max(zero(eltype(mu_array)), m) + mu_array[i] > 0 && (any_stalled = true) + end + any_stalled || return false + gamma_target .= gamma + one_t, zero_t = one(eltype(viscosity_matrix)), zero(eltype(viscosity_matrix)) + @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 + ldiv!(gamma, lu!(viscosity_matrix), gamma_target) + return true +end + """ gamma_loop!(solver::Solver, AIC_x::Matrix{Float64}, AIC_y::Matrix{Float64}, AIC_z::Matrix{Float64}, @@ -1047,23 +1081,11 @@ function gamma_loop!( ) # 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 + apply_artificial_viscosity!( + gamma_new, panels, solver.lr.alpha_dist, laplacian, + viscosity_matrix, lift_slope, mu_array, gamma_target, + planform_area, solver.artificial_viscosity_factor, + ) end # Update gamma with relaxation and damping diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 9d2c5b05..d5abf838 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -109,42 +109,77 @@ end @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 +# A deep-alpha polar is required: the viscosity only fires where the local lift +# slope is negative, and `cl` clamps flat past the tabulated alpha range. The +# default 15-deg range never stalls, so we extend it to 40 deg. +poststall_wing = ram_air_matrix_wing(; n_panels=20, n_sections=4, + alpha_range=deg2rad.(-5:2:40), delta_range=deg2rad.(-3:3:3)) +refine!(poststall_wing) + +roughness(v) = sum(abs, @views v[1:end-2] .- 2 .* v[2:end-1] .+ v[3:end]) + +@testset "apply_artificial_viscosity! smooths post-stall, no-op attached" begin + body_aero = BodyAerodynamics([poststall_wing]) + panels = body_aero.panels + n = length(panels) - base_converged, _ = iterate_loop(false, 0.1) - av_converged, gamma_av = iterate_loop(true, 0.5) + laplacian = zeros(n, n) + VortexStepMethod.build_spanwise_laplacian!(laplacian, n) + viscosity_matrix = zeros(n, n) + lift_slope = zeros(n) + mu_array = zeros(n) + gamma_target = zeros(n) + planform_area = sum(p.width * p.chord for p in panels) + + spiky() = [isodd(i) ? 1.0 : -1.0 for i in 1:n] + attached = fill(deg2rad(2.0), n) + post_stall = fill(deg2rad(22.0), n) + + # Attached flow: the local slope is positive everywhere, so mu stays zero, + # the solve is skipped, and gamma is returned untouched. + gamma_attached = spiky() + fired_attached = VortexStepMethod.apply_artificial_viscosity!(gamma_attached, + panels, attached, laplacian, viscosity_matrix, lift_slope, mu_array, + gamma_target, planform_area, 0.035) + @test !fired_attached + @test gamma_attached == spiky() + + # Post-stall: the solve fires and smooths the sawtooth. + gamma_stalled = spiky() + rough_before = roughness(gamma_stalled) + fired_stalled = VortexStepMethod.apply_artificial_viscosity!(gamma_stalled, + panels, post_stall, laplacian, viscosity_matrix, lift_slope, mu_array, + gamma_target, planform_area, 0.035) + @test fired_stalled + @test roughness(gamma_stalled) < rough_before + + # The attached (hot) path must not allocate. + gamma_alloc = spiky() + VortexStepMethod.apply_artificial_viscosity!(gamma_alloc, panels, attached, + laplacian, viscosity_matrix, lift_slope, mu_array, gamma_target, + planform_area, 0.035) + allocs = @allocated VortexStepMethod.apply_artificial_viscosity!(gamma_alloc, + panels, attached, laplacian, viscosity_matrix, lift_slope, mu_array, + gamma_target, planform_area, 0.035) + @test allocs == 0 +end - @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 +@testset "solve! artificial viscosity: attached no-op, post-stall finite" begin + body_aero = BodyAerodynamics([poststall_wing]) + solver_off = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM, + is_with_artificial_viscosity=false) + solver_on = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM, + is_with_artificial_viscosity=true) + + # Attached flow: viscosity never fires, so results are bit-identical. + set_va!(body_aero, [10.0, 0.0, 0.0]) + gamma_off = copy(solve!(solver_off, body_aero).gamma_distribution) + set_va!(body_aero, [10.0, 0.0, 0.0]) + gamma_on = copy(solve!(solver_on, body_aero).gamma_distribution) + @test gamma_on == gamma_off + + # High angle of attack: the viscosity path runs and stays finite. + set_va!(body_aero, [10.0 * cosd(25), 0.0, 10.0 * sind(25)]) + sol = solve!(solver_on, body_aero) + @test all(isfinite, sol.gamma_distribution) end