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 3646579c..4cea87ff 100644 --- a/docs/src/private_functions.md +++ b/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! ``` \ No newline at end of file diff --git a/examples/billowing.jl b/examples/billowing.jl index abb27b27..e4c2835f 100644 --- a/examples/billowing.jl +++ b/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, 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 7f3b7a42..b6a5b594 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/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 diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 7f21ed28..dc5d59f4 100644 --- a/src/body_aerodynamics.jl +++ b/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 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 ebf021e9..cf3edb46 100644 --- a/src/solver.jl +++ b/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), 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