diff --git a/changelog.d/395-the-curved-reference-wing-does-not-conve.md b/changelog.d/395-the-curved-reference-wing-does-not-conve.md new file mode 100644 index 00000000..e7e1fa43 --- /dev/null +++ b/changelog.d/395-the-curved-reference-wing-does-not-conve.md @@ -0,0 +1,7 @@ +### Fixed + +- The `LOOP` solver caps each panel's relaxation at `1 / (1 + π c |z_airf ⋅ AIC_ii|)`, the diagonal step of its self-induced stiffness, so narrow panels (cosine-spaced tips) no longer make it diverge to NaN; the curved reference wing and finely refined wings now converge at the default `relaxation_factor`, with unchanged iteration counts on wings that already converged. A `LOOP` solve that still does not converge keeps halving `relaxation_factor` while it is above 1e-3 and warns once with the factor it stopped at, and each attempt stops at the first non-finite residual. + +### Changed + +- `cl_unrefined_dist`, `cd_unrefined_dist`, `cm_unrefined_dist` and `alpha_unrefined_dist` are area-weighted means over the section's panels instead of plain means, and `chord_unrefined_dist` is the section area over its width, so a coefficient times `chord_unrefined_dist * width_unrefined_dist` is the section's load. `moment_unrefined_dist` is the section's moment, summed over its panels, instead of their mean; `linearize` with `aero_coeffs=false` returns it. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index d8b0c836..53d0b9ad 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -17,6 +17,7 @@ solver_kwargs check_dimensions calculate_AIC_matrices! gamma_loop! +panel_relaxation_factors! build_spanwise_laplacian! local_lift_slope! apply_artificial_viscosity! diff --git a/src/solver.jl b/src/solver.jl index 871e75bb..c507ef10 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -26,12 +26,17 @@ Struct for storing the solution of the [`solve!`](@ref) function. Must contain a - `moment_coeffs`::MVec3: Aerodynamic moment coefficients [CMx, CMy, CMz] [-] - `moment_dist`::Vector{Float64}: Pitching moments around the spanwise vector of each panel. [Nm] - `moment_coeff_dist`::Vector{Float64}: Pitching moment coefficient around the spanwise vector of each panel. [-] -- `moment_unrefined_dist`::MVector{U, Float64}: Averaged moments for unrefined sections [Nm] -- `cl_unrefined_dist`::MVector{U, Float64}: Averaged lift coefficients for unrefined sections [-] -- `cd_unrefined_dist`::MVector{U, Float64}: Averaged drag coefficients for unrefined sections [-] -- `cm_unrefined_dist`::MVector{U, Float64}: Averaged airfoil moment coefficients for unrefined sections [-] +- `moment_unrefined_dist`::MVector{U, Float64}: Pitching moment per unrefined section, + summed over its panels [Nm] +- `cl_unrefined_dist`::MVector{U, Float64}: Area-weighted mean lift coefficient + per unrefined section [-] +- `cd_unrefined_dist`::MVector{U, Float64}: Area-weighted mean drag coefficient + per unrefined section [-] +- `cm_unrefined_dist`::MVector{U, Float64}: Area-weighted mean airfoil moment coefficient + per unrefined section [-] - `moment_coeff_unrefined_dist`::MVector{U, Float64}: Summed `moment_frac`-referenced pitching-moment coefficient per unrefined section [-] -- `alpha_unrefined_dist`::MVector{U, Float64}: Averaged angles of attack for unrefined sections [rad] +- `alpha_unrefined_dist`::MVector{U, Float64}: Area-weighted mean angle of attack + per unrefined section [rad] - `solver_status`::SolverStatus: enum, see [`SolverStatus`](@ref) """ @with_kw mutable struct VSMSolution{P, U, T} @@ -181,7 +186,8 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref) # Intermediate results lr::LoopResult{P, T} = LoopResult{P, T}() br::BaseResult{P, T} = BaseResult{P, T}() - cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:8] + cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = + [LazyBufferCache() for _ in 1:9] cache_base::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache()] cache_lin::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:4] @@ -498,42 +504,45 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; original_section_idx = wing.refined_panel_mapping[local_panel_idx] target_unrefined_idx = section_range[original_section_idx] + area = panel.chord * panel.width + # Accumulate coefficients and moments moment_unrefined_dist[target_unrefined_idx] += moment_dist[panel_idx] - cl_unrefined_dist[target_unrefined_idx] += solver.sol.cl_dist[panel_idx] - cd_unrefined_dist[target_unrefined_idx] += solver.sol.cd_dist[panel_idx] - cm_unrefined_dist[target_unrefined_idx] += solver.sol.cm_dist[panel_idx] + cl_unrefined_dist[target_unrefined_idx] += + solver.sol.cl_dist[panel_idx] * area + cd_unrefined_dist[target_unrefined_idx] += + solver.sol.cd_dist[panel_idx] * area + cm_unrefined_dist[target_unrefined_idx] += + solver.sol.cm_dist[panel_idx] * area moment_coeff_unrefined_dist[target_unrefined_idx] += solver.sol.moment_coeff_dist[panel_idx] - alpha_unrefined_dist[target_unrefined_idx] += solver.sol.alpha_dist[panel_idx] + alpha_unrefined_dist[target_unrefined_idx] += + solver.sol.alpha_dist[panel_idx] * area # Accumulate geometry x_airf_unrefined_dist[target_unrefined_idx] .+= panel.x_airf y_airf_unrefined_dist[target_unrefined_idx] .+= panel.y_airf z_airf_unrefined_dist[target_unrefined_idx] .+= panel.z_airf va_vec_unrefined_dist[target_unrefined_idx] .+= panel.va_vec - chord_unrefined_dist[target_unrefined_idx] += panel.chord + chord_unrefined_dist[target_unrefined_idx] += area width_unrefined_dist[target_unrefined_idx] += panel.width unrefined_count_dist[target_unrefined_idx] += 1 end - # Average coefficients and geometry. width and - # moment_coeff_unrefined_dist stay summed (extensive). for target_unrefined_idx in section_range if unrefined_count_dist[target_unrefined_idx] > 0 count = unrefined_count_dist[target_unrefined_idx] - moment_unrefined_dist[target_unrefined_idx] /= count - cl_unrefined_dist[target_unrefined_idx] /= count - cd_unrefined_dist[target_unrefined_idx] /= count - cm_unrefined_dist[target_unrefined_idx] /= count - alpha_unrefined_dist[target_unrefined_idx] /= count + area = chord_unrefined_dist[target_unrefined_idx] + cl_unrefined_dist[target_unrefined_idx] /= area + cd_unrefined_dist[target_unrefined_idx] /= area + cm_unrefined_dist[target_unrefined_idx] /= area + alpha_unrefined_dist[target_unrefined_idx] /= area x_airf_unrefined_dist[target_unrefined_idx] ./= count y_airf_unrefined_dist[target_unrefined_idx] ./= count z_airf_unrefined_dist[target_unrefined_idx] ./= count va_vec_unrefined_dist[target_unrefined_idx] ./= count - chord_unrefined_dist[target_unrefined_idx] /= count - # width_unrefined_dist is NOT averaged - it is the - # sum of panel widths in the unrefined section + chord_unrefined_dist[target_unrefined_idx] = + area / width_unrefined_dist[target_unrefined_idx] end end end @@ -652,8 +661,9 @@ Converge the circulation distribution and leave it in `solver.lr.gamma_new`, wit turning it into forces. Fills the solver's panel arrays, builds the AIC matrices, starts from `gamma_distribution` (or an elliptical/zero guess when it is `nothing` or `solver.use_gamma_prev` is false) and iterates; a `LOOP` solver that fails to converge -retries once with half the relaxation factor. The circulation half of [`solve!`](@ref), -paired with [`calc_forces!`](@ref). Returns `nothing`. +retries with half the relaxation factor while the factor is above 1e-3, and warns once +with the factor it stopped at. The circulation half of [`solve!`](@ref), paired with +[`calc_forces!`](@ref). Returns `nothing`. """ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_distribution=nothing; log=false) where {P, U, T} @@ -721,11 +731,15 @@ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma solver.lr.gamma_new .= gamma_initial # Run main iteration loop gamma_loop!(solver, body_aero, panels, relaxation_factor; log) - # Try again with reduced relaxation factor if not converged - if solver.solver_type == LOOP && !solver.lr.converged && relaxation_factor > 1e-3 - log && @warn "Running again with half the relaxation_factor = $(relaxation_factor/2)" + initial_relaxation_factor = relaxation_factor + while solver.solver_type == LOOP && !solver.lr.converged && relaxation_factor > 1e-3 + relaxation_factor /= 2 solver.lr.gamma_new .= gamma_initial - gamma_loop!(solver, body_aero, panels, relaxation_factor/2; log) + gamma_loop!(solver, body_aero, panels, relaxation_factor; log) + end + if relaxation_factor < initial_relaxation_factor + @warn "LOOP halved relaxation_factor from $initial_relaxation_factor to " * + "$relaxation_factor, converged: $(solver.lr.converged)" end nothing @@ -883,15 +897,34 @@ function apply_artificial_viscosity!(gamma, panels, alpha_dist, laplacian, visco 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) + panel_relaxation_factors!(panel_relaxation, relaxation_factor, AIC, z_airf_dist, + chord_dist) + +Fill `panel_relaxation` with the LOOP relaxation factor of each panel: `relaxation_factor`, +capped at `1 / stiffness_i`, the diagonal step of the panel's self-induced stiffness +`stiffness_i = 1 + π c_i |z_airf_i ⋅ AIC[i, i, :]|`. +""" +function panel_relaxation_factors!(panel_relaxation, relaxation_factor, AIC, z_airf_dist, + chord_dist) + @inbounds for i in eachindex(panel_relaxation) + normal_self_induction = z_airf_dist[i, 1] * AIC[i, i, 1] + + z_airf_dist[i, 2] * AIC[i, i, 2] + z_airf_dist[i, 3] * AIC[i, i, 3] + stiffness = 1 + π * chord_dist[i] * abs(normal_self_induction) + panel_relaxation[i] = min(relaxation_factor, 1 / stiffness) + end + return panel_relaxation +end + +""" + gamma_loop!(solver::Solver, body_aero::BodyAerodynamics, + panels::AbstractVector{<:Panel}, relaxation_factor; log=true) Main iteration loop for calculating circulation distribution. Both solvers converge on the fixed-point residual `F(gamma) - gamma`, measured -relative to the largest circulation: the LOOP solver takes under-relaxed steps -towards `F(gamma)`, the NONLIN solver a Newton step with a finite-difference +relative to the largest circulation: the LOOP solver steps towards `F(gamma)`, each +panel under-relaxed by [`panel_relaxation_factors!`](@ref), and stops at the first +non-finite residual; the NONLIN solver takes a Newton step with a finite-difference Jacobian, backtracked until it reduces the residual. When `solver.is_with_artificial_viscosity` is set, the LOOP solver replaces the @@ -925,6 +958,7 @@ function gamma_loop!( relative_velocity_crossz = solver.cache[6][va_vec_dist] v_acrossz_dist = solver.cache[7][va_vec_dist] cl_dist = solver.cache[8][solver.lr.gamma_new] + panel_relaxation = solver.cache[9][solver.lr.gamma_new] AIC_x = @view body_aero.AIC[:, :, 1] AIC_y = @view body_aero.AIC[:, :, 2] @@ -1042,6 +1076,8 @@ function gamma_loop!( planform_area += panels[i].width * chord_dist[i] end end + panel_relaxation_factors!(panel_relaxation, relaxation_factor, body_aero.AIC, + z_airf_dist, chord_dist) function f_loop!(gamma_new, gamma) gamma .= gamma_new @@ -1077,7 +1113,7 @@ function gamma_loop!( ) end - @. gamma_new = (1 - relaxation_factor) * gamma + relaxation_factor * gamma_new + @. gamma_new = (1 - panel_relaxation) * gamma + panel_relaxation * gamma_new return nothing end iters = 0 @@ -1090,12 +1126,13 @@ function gamma_loop!( abs_gamma_new .= abs.(solver.lr.gamma_new) reference_error = maximum(abs_gamma_new) reference_error = max(reference_error, solver.tol_reference_error) - # The relaxed step is `relaxation_factor` times the fixed-point residual. - abs_gamma_new .= abs.(solver.lr.gamma_new .- gamma) - residual = maximum(abs_gamma_new) / relaxation_factor + # The relaxed step is `panel_relaxation` times the fixed-point residual. + abs_gamma_new .= abs.(solver.lr.gamma_new .- gamma) ./ panel_relaxation + residual = maximum(abs_gamma_new) normalized_error = residual / reference_error @debug "Iteration: $i, normalized_error: $normalized_error" + isfinite(normalized_error) || break if normalized_error < solver.rtol solver.lr.converged = true diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index 893853c4..10d54532 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/test/body_aerodynamics/test_body_aerodynamics.jl @@ -292,6 +292,19 @@ end end end end + + # CL from the NONLIN solve of the same wing + @testset "VSM LOOP converges on the curved wing at $(alpha_deg)°" for + (alpha_deg, CL) in ((5.7106, 0.422308), (10.0, 0.726321)) + body_aero, _, _, _ = create_geometry(wing_type=:curved) + set_va!(body_aero, 20.0 .* [cosd(alpha_deg), 0.0, sind(alpha_deg)]) + wing = body_aero.wings[1] + solver = Solver(wing.n_panels, wing.n_unrefined_sections; + aerodynamic_model_type=VSM) + sol = solve!(solver, body_aero) + @test sol.solver_status == FEASIBLE + @test sol.force_coeffs[3] ≈ CL rtol=1e-4 + end end @testset "Calculate results against output results" begin diff --git a/test/solver/solver_test_polar.csv b/test/solver/solver_test_polar.csv index 7e043d9c..dd714098 100644 --- a/test/solver/solver_test_polar.csv +++ b/test/solver/solver_test_polar.csv @@ -1,4 +1,5 @@ alpha,cl,cd,cm +-20.0,-0.8,0.04,-0.02 -10.0,0.1,0.02,-0.01 0.0,1.0,0.01,0.0 10.0,1.8,0.025,0.01 diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 25ebf514..506cd66d 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -175,6 +175,18 @@ end end end +@testset "LOOP caps only the panels stiffer than 1 / relaxation_factor" begin + AIC = zeros(2, 2, 3) + AIC[1, 1, 3] = -1 / π # stiffness 1 + π * 1.0 * (1/π) = 2 + AIC[2, 2, 3] = 99 / (2π) # stiffness 1 + π * 2.0 * (99/2π) = 100 + z_airf_dist = [0.0 0.0 1.0; 0.0 0.0 1.0] + chord_dist = [1.0, 2.0] + panel_relaxation = VortexStepMethod.panel_relaxation_factors!(zeros(2), 0.03, AIC, + z_airf_dist, chord_dist) + @test panel_relaxation[1] === 0.03 + @test panel_relaxation[2] ≈ 1 / 100 +end + calc_forces_allocs(solver, body_aero) = (calc_forces!(solver, body_aero); @allocated calc_forces!(solver, body_aero)) diff --git a/test/solver/test_unrefined_dist.jl b/test/solver/test_unrefined_dist.jl index 9bca8a5c..d601cead 100644 --- a/test/solver/test_unrefined_dist.jl +++ b/test/solver/test_unrefined_dist.jl @@ -2,6 +2,15 @@ using VortexStepMethod using LinearAlgebra using Test +""" + area_mean(values, panels, indices) + +Mean of `values[indices]` weighted by the area of `panels[indices]`. +""" +area_mean(values, panels, indices) = + sum(values[i] * panels[i].chord * panels[i].width for i in indices) / + sum(panels[i].chord * panels[i].width for i in indices) + @testset "Unrefined Arrays Tests" begin @testset "Unrefined section array aggregation" begin # Create a simple wing with unrefined sections @@ -28,6 +37,7 @@ using Test va_vec = [10.0, 0.0, 0.0] set_va!(body_aero, va_vec) sol = solve!(solver, body_aero) + @test sol.solver_status == FEASIBLE # Test 1: Unrefined arrays exist and have correct size @test length(sol.cl_unrefined_dist) == wing.n_unrefined_sections @@ -38,50 +48,33 @@ using Test @test !all(sol.cl_unrefined_dist .== 0.0) @test !all(sol.cd_unrefined_dist .== 0.0) - # Test 3: Verify unrefined coefficients are averaged from refined panels + # Test 3: Unrefined coefficients are area-weighted means of their refined panels # refined_panel_mapping maps each refined panel to its unrefined section index + panels = body_aero.panels for unrefined_idx in 1:wing.n_unrefined_sections # Find all refined panels that map to this unrefined section refined_panel_indices = findall(x -> x == unrefined_idx, wing.refined_panel_mapping) if !isempty(refined_panel_indices) - # Calculate expected average from refined panel coefficients - expected_cl = sum(sol.cl_dist[refined_panel_indices]) / length(refined_panel_indices) - expected_cd = sum(sol.cd_dist[refined_panel_indices]) / length(refined_panel_indices) - expected_cm = sum(sol.cm_dist[refined_panel_indices]) / length(refined_panel_indices) - - # Check if unrefined coefficients match expected averages - # Handle NaN values that can occur in INVISCID models - if isnan(expected_cl) - @test isnan(sol.cl_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cl_unrefined_dist[unrefined_idx], expected_cl, rtol=1e-10) - end - if isnan(expected_cd) - @test isnan(sol.cd_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cd_unrefined_dist[unrefined_idx], expected_cd, rtol=1e-10) - end - if isnan(expected_cm) - @test isnan(sol.cm_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cm_unrefined_dist[unrefined_idx], expected_cm, rtol=1e-10) - end + expected_cl = area_mean(sol.cl_dist, panels, refined_panel_indices) + expected_cd = area_mean(sol.cd_dist, panels, refined_panel_indices) + expected_cm = area_mean(sol.cm_dist, panels, refined_panel_indices) + + @test isapprox(sol.cl_unrefined_dist[unrefined_idx], expected_cl, rtol=1e-10) + @test isapprox(sol.cd_unrefined_dist[unrefined_idx], expected_cd, rtol=1e-10) + @test isapprox(sol.cm_unrefined_dist[unrefined_idx], expected_cm, rtol=1e-10) end end # Test 4: Verify physical consistency (lift coefficients should be positive at positive AoA) - # Skip test if values are NaN - if !any(isnan.(sol.cl_unrefined_dist)) - @test all(sol.cl_unrefined_dist .> 0.0) - end + @test all(sol.cl_unrefined_dist .> 0.0) finally rm(settings_file; force=true) end end - @testset "Unrefined force consistency: coeff*width*chord" begin + @testset "Unrefined section loads equal the sums of their panel loads" begin # Verify that for each unrefined section: # coeff_unrefined * width_unrefined * chord_unrefined ≈ # sum(coeff_panel * width_panel * chord_panel) for all panels in section @@ -102,6 +95,7 @@ using Test va_vec = [10.0, 0.0, 0.0] set_va!(body_aero, va_vec) sol = solve!(solver, body_aero) + @test sol.solver_status == FEASIBLE panels = body_aero.panels for unrefined_idx in 1:wing.n_unrefined_sections @@ -133,17 +127,11 @@ using Test panels[i].width for i in pidxs) @test w_u ≈ expected_width rtol=1e-12 - # Force-like quantities should be approximately - # equal (not exact due to averaging) - if !isnan(cl_u) - @test cl_u * c_u * w_u ≈ sum_cl_cw rtol=0.05 - end - if !isnan(cd_u) - @test cd_u * c_u * w_u ≈ sum_cd_cw rtol=0.05 - end - if !isnan(cm_u) - @test cm_u * c_u * w_u ≈ sum_cm_cw rtol=0.05 - end + @test cl_u * c_u * w_u ≈ sum_cl_cw rtol=1e-12 + @test cd_u * c_u * w_u ≈ sum_cd_cw rtol=1e-12 + @test cm_u * c_u * w_u ≈ sum_cm_cw rtol=1e-12 + @test sol.moment_unrefined_dist[unrefined_idx] ≈ + sum(sol.moment_dist[i] for i in pidxs) rtol=1e-12 end finally rm(settings_file; force=true) @@ -176,6 +164,7 @@ using Test va_vec = [10.0, 0.0, 0.0] set_va!(body_aero, va_vec) sol = solve!(solver, body_aero) + @test sol.solver_status == FEASIBLE # Verify arrays have correct size @test length(sol.cl_unrefined_dist) == wing.n_unrefined_sections @@ -183,30 +172,18 @@ using Test @test length(sol.cm_unrefined_dist) == wing.n_unrefined_sections # Verify unrefined coefficients are computed correctly using mapping + panels = body_aero.panels for unrefined_idx in 1:wing.n_unrefined_sections refined_panel_indices = findall(x -> x == unrefined_idx, wing.refined_panel_mapping) if !isempty(refined_panel_indices) - expected_cl = sum(sol.cl_dist[refined_panel_indices]) / length(refined_panel_indices) - expected_cd = sum(sol.cd_dist[refined_panel_indices]) / length(refined_panel_indices) - expected_cm = sum(sol.cm_dist[refined_panel_indices]) / length(refined_panel_indices) - - # Handle NaN for all coefficients - if isnan(expected_cl) - @test isnan(sol.cl_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cl_unrefined_dist[unrefined_idx], expected_cl, rtol=1e-10) - end - if isnan(expected_cd) - @test isnan(sol.cd_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cd_unrefined_dist[unrefined_idx], expected_cd, rtol=1e-10) - end - if isnan(expected_cm) - @test isnan(sol.cm_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cm_unrefined_dist[unrefined_idx], expected_cm, rtol=1e-10) - end + expected_cl = area_mean(sol.cl_dist, panels, refined_panel_indices) + expected_cd = area_mean(sol.cd_dist, panels, refined_panel_indices) + expected_cm = area_mean(sol.cm_dist, panels, refined_panel_indices) + + @test isapprox(sol.cl_unrefined_dist[unrefined_idx], expected_cl, rtol=1e-10) + @test isapprox(sol.cd_unrefined_dist[unrefined_idx], expected_cd, rtol=1e-10) + @test isapprox(sol.cm_unrefined_dist[unrefined_idx], expected_cm, rtol=1e-10) end end