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/docs/src/private_functions.md b/docs/src/private_functions.md index 2e96c6d9..fe6963d6 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -8,6 +8,9 @@ CurrentModule = VortexStepMethod ```@docs calculate_AIC_matrices! gamma_loop! +build_spanwise_laplacian! +local_lift_slope! +apply_artificial_viscosity! frozen_wake! calc_forces! calculate_cl diff --git a/examples/billowing.jl b/examples/billowing.jl index fa42e8c0..219b7e54 100644 --- a/examples/billowing.jl +++ b/examples/billowing.jl @@ -8,8 +8,8 @@ using MakieControlPlots 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) @@ -119,7 +119,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) @@ -142,8 +142,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/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index 2116ae37..4055eb2b 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -600,19 +600,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 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..0e9f0915 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,97 @@ 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 + +""" + 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}, 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 +1036,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 +1079,19 @@ function gamma_loop!( cl_dist, chord_array, ) + # Gate the linear solve on any(mu > 0): fires exactly in post-stall. + if use_viscosity + 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 - @. 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 +1254,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..d5abf838 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -91,3 +91,95 @@ 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 + +# 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) + + 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 + +@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