LOOP caps each panel's relaxation at its own stiffness, so narrow tip panels no longer diverge; unrefined section loads are area-weighted

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.

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!

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

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

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

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))