Add stability_derivatives (angle of attack, sideslip) and trim_angle, built on linearize

CHANGELOG.md

@@ -4,6 +4,11 @@
 
 ### Added
 
+- `stability_derivatives` gives the force and moment coefficients and their derivatives
+  with respect to angle of attack and sideslip, and `trim_angle` the angles of attack at
+  which `CMy` changes sign, with the slope that says whether each trim is stable.
+- `apparent_wind(alpha, beta, wind_speed)` gives the body-frame inflow vector at an angle
+  of attack and sideslip, as `set_va!(body_aero, settings)` sets it.
 - `Solver(settings)` and `Solver(n_panels, n_unrefined_sections)` build a solver without
   a `BodyAerodynamics`; keyword arguments override the settings.
 - `set_va!(body_aero, va_vec, omega; reference_point)` turns the body about

README.md

@@ -131,7 +131,7 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+va_vec = apparent_wind(alpha, 0.0, va)
 set_va!(body_aero, va_vec)
 ```
 It is possible to import the wing geometry using an `.obj` file as shown in the example `ram_air_kite.jl`. During the import the polars are calculated automatically using XFoil. This approach is valid for rigid wings and ram-air kites, but not for leading edge inflatable kites.

docs/src/examples.md

@@ -77,7 +77,7 @@ multiple wings.
 ###### Set inflow conditions
 
 ```julia
-julia> va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+julia> va_vec = apparent_wind(alpha, 0.0, va)
 julia> set_va!(body_aero, va_vec, [0, 0, 0.1])
 ```
 

docs/src/functions.md

@@ -96,12 +96,15 @@ CurrentModule = VortexStepMethod
 ## Setting the inflow conditions and solving
 ```@docs
 set_va!
+apparent_wind
 section_pitch_rate
 solve
 solve!
 solve_base!
 reinit!(body_aero::BodyAerodynamics{P, W, T}) where {P, W, T}
 linearize
+stability_derivatives
+trim_angle
 calculate_results
 ```
 

docs/src/index.md

@@ -137,7 +137,7 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+va_vec = apparent_wind(alpha, 0.0, va)
 set_va!(body_aero, va_vec)
 ```
 

docs/src/private_functions.md

@@ -28,6 +28,9 @@ calculate_cd
 calculate_cm
 calculate_cd_cm
 set_pitch_rate_dist!
+coeffs_at_angles
+nose_down
+bisect_sign_change
 calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!

examples/bench.jl

@@ -37,7 +37,7 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+va_vec = apparent_wind(alpha, 0.0, va)
 set_va!(body_aero, va_vec)
 
 # Step 4: Initialize solvers for both LLT and VSM methods
@@ -76,13 +76,8 @@ vsm_solver = Solver(
 va = 15.0
 aoa = 15.0
 side_slip = 0.0
-yaw_rate = 0.0
 aoa_rad = deg2rad(aoa)
-va_vec = [
-    cos(aoa_rad) * cos(side_slip),
-    sin(side_slip),
-    sin(aoa_rad)
-] * va
+va_vec = apparent_wind(aoa_rad, side_slip, va)
 set_va!(body_aero, va_vec)
 
 # Solving

examples/billowing.jl

@@ -89,7 +89,7 @@ sideslip_deg = settings.condition.beta
 
 α0 = deg2rad(angle_of_attack_deg)
 β0 = deg2rad(sideslip_deg)
-va_vec = wind_speed .* [cos(α0) * cos(β0), sin(β0), sin(α0) * cos(β0)]
+va_vec = apparent_wind(α0, β0, wind_speed)
 set_va!(body_aero_flat, va_vec)
 set_va!(body_aero_bill, va_vec)
 

examples/linearize_check.jl

@@ -35,11 +35,7 @@ va         = 15.0
 aoa_deg    = 10.0
 aoa_rad    = deg2rad(aoa_deg)
 side_slip  = 0.0
-va_vec_b_0 = [
-    cos(aoa_rad) * cos(side_slip),
-    sin(side_slip),
-    sin(aoa_rad),
-] * va
+va_vec_b_0 = apparent_wind(aoa_rad, side_slip, va)
 omega_b_0  = zeros(3)
 theta_0    = zeros(n_unrefined)
 

examples/obj_to_yaml_kite.jl

@@ -74,7 +74,7 @@ VortexStepMethod.reinit!(body_aero)
 
 solver = Solver(wing.n_panels, wing.n_unrefined_sections;
     aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP)
-set_va!(body_aero, [cos(deg2rad(8)) * va, 0.0, sin(deg2rad(8)) * va])
+set_va!(body_aero, apparent_wind(deg2rad(8), 0.0, va))
 results = VortexStepMethod.solve(solver, body_aero; log=true)
 
 if PLOT

examples/rectangular_wing.jl

@@ -41,7 +41,7 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+va_vec = apparent_wind(alpha, 0.0, va)
 set_va!(body_aero, va_vec, [0, 0, 0.1])
 
 # Step 4: Initialize solvers for both LLT and VSM methods

examples/stall_model.jl

@@ -61,11 +61,7 @@ aoa = 17.0
 side_slip = 0.0
 yaw_rate = 0.0
 aoa_rad = deg2rad(aoa)
-va_vec = [
-    cos(aoa_rad) * cos(side_slip),
-    sin(side_slip),
-    sin(aoa_rad)
-] * va
+va_vec = apparent_wind(aoa_rad, side_slip, va)
 set_va!(body_aero, va_vec)
 
 # Plotting geometry

src/VortexStepMethod.jl

@@ -32,9 +32,10 @@ export slice_args, preview_args
 export ObjWing, Section, Wing, refine!, reinit!
 export BodyAerodynamics
 export Solver, VSMSolution, linearize, solve, solve!, solve_base!, calc_forces!
+export stability_derivatives, trim_angle
 export SolveFailure
 export calculate_results
-export add_section!, set_va!, section_pitch_rate
+export add_section!, apparent_wind, set_va!, section_pitch_rate
 export calculate_projected_area, calculate_span
 export MVec3
 
@@ -437,6 +438,7 @@ include("panel.jl")
 include("body_aerodynamics.jl")
 include("wake.jl")
 include("solver.jl")
+include("stability.jl")
 
 include("plotting_helpers.jl")
 

src/body_aerodynamics.jl

@@ -1158,26 +1158,20 @@ function set_va!(body_aero::BodyAerodynamics, va_vec_dist::AbstractMatrix;
 end
 
 """
-    set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
-
-Set velocity array from VSM settings configuration.
+    apparent_wind(alpha, beta, wind_speed)
 
-This convenience method extracts flight conditions from VSMSettings and 
-constructs the velocity vector in the body reference frame based on:
-- Wind speed from settings.condition.wind_speed
-- Angle of attack from settings.condition.alpha (converted from degrees)
-- Sideslip angle from settings.condition.beta (converted from degrees)
-- Yaw rate from settings.condition.yaw_rate (converted from °/s), applied as `omega`
-  about Z_b and turning the body about `body_aero.reference_point`
+Apparent wind vector in the body frame [m/s] at angle of attack `alpha` [rad], sideslip
+`beta` [rad] and `wind_speed` [m/s].
+"""
+apparent_wind(alpha, beta, wind_speed) =
+    wind_speed .* [cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)]
 
-The velocity vector is constructed as:
-- X_b (forward): wind_speed * cos(α) * cos(β)  
-- Y_b (right): wind_speed * sin(β)
-- Z_b (down): wind_speed * sin(α) * cos(β)
+"""
+    set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
 
-# Arguments
-- `body_aero::BodyAerodynamics`: The aerodynamic body to modify
-- `settings::VSMSettings`: Settings object containing flight conditions
+Set the uniform inflow of `body_aero` to the [`apparent_wind`](@ref) at the `alpha` and
+`beta` [°] and `wind_speed` [m/s] of `settings.condition`, turning the body about
+`body_aero.reference_point` at its `yaw_rate` [°/s] about Z_b.
 
 # Example
 ```julia
@@ -1187,16 +1181,8 @@ set_va!(body_aero, settings)
 ```
 """
 function set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
-    α = deg2rad(settings.condition.alpha)
-    β = deg2rad(settings.condition.beta)
-    wind_speed = settings.condition.wind_speed
-
-    va_vec = wind_speed * [
-        cos(α)*cos(β),  # X_b (forward)
-        sin(β),         # Y_b (right)
-        sin(α)*cos(β)   # Z_b (down)
-    ]
-    omega = [0.0, 0.0, deg2rad(settings.condition.yaw_rate)]
-
-    set_va!(body_aero, va_vec, omega)
+    condition = settings.condition
+    va_vec = apparent_wind(deg2rad(condition.alpha), deg2rad(condition.beta),
+        condition.wind_speed)
+    set_va!(body_aero, va_vec, [0.0, 0.0, deg2rad(condition.yaw_rate)])
 end

src/stability.jl

@@ -0,0 +1,92 @@
+"""
+    stability_derivatives(solver, body_aero, alpha, beta, wind_speed; kwargs...)
+
+Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` at angle of attack
+`alpha` [rad], sideslip `beta` [rad] and `wind_speed` [m/s], and their derivatives with
+respect to `alpha` and `beta` [1/rad], at the rotation rate `body_aero.omega` and with
+moments about `solver.reference_point`. `kwargs` go to [`linearize`](@ref), which leaves
+`body_aero` at this inflow.
+
+Returns `(coeffs, dalpha, dbeta, converged)`.
+"""
+function stability_derivatives(solver::Solver, body_aero::BodyAerodynamics, alpha, beta,
+        wind_speed; kwargs...)
+    va_vec = apparent_wind(alpha, beta, wind_speed)
+    jac, results, converged = linearize(solver, body_aero, va_vec;
+        theta_idxs=nothing, va_idxs=1:3, aero_coeffs=true, kwargs...)
+    dva_dalpha = ForwardDiff.derivative(
+        angle -> apparent_wind(angle, beta, wind_speed), alpha)
+    dva_dbeta = ForwardDiff.derivative(
+        angle -> apparent_wind(alpha, angle, wind_speed), beta)
+    coeff_jac = jac[1:6, :]
+    return (coeffs=results[1:6], dalpha=coeff_jac * dva_dalpha,
+        dbeta=coeff_jac * dva_dbeta, converged)
+end
+
+"""
+    trim_angle(solver, body_aero, beta, wind_speed; alpha_range=deg2rad.(-5:2:15),
+               alpha_tol=1e-5, backend=AutoForwardDiff())
+
+Angles of attack [rad] at which `CMy` of `body_aero` about `solver.reference_point` changes
+sign between neighbouring entries of `alpha_range`, bisected to `alpha_tol` [rad], at
+sideslip `beta` [rad] and `wind_speed` [m/s]. Returns one `(alpha, dCMy_dalpha)` per trim,
+the slope [1/rad] from [`stability_derivatives`](@ref) with `backend`; a trim is statically
+stable where `dCMy_dalpha < 0`. Throws a [`SolveFailure`](@ref) if a solve misses the
+solver's tolerances.
+"""
+function trim_angle(solver::Solver, body_aero::BodyAerodynamics, beta, wind_speed;
+        alpha_range=deg2rad.(-5:2:15), alpha_tol=1e-5, backend=AutoForwardDiff())
+    is_nose_down = alpha -> nose_down(solver, body_aero, alpha, beta, wind_speed)
+    nose_down_range = is_nose_down.(alpha_range)
+    trims = @NamedTuple{alpha::Float64, dCMy_dalpha::Float64}[]
+    for i in 1:length(alpha_range)-1
+        nose_down_range[i] == nose_down_range[i+1] && continue
+        alpha = bisect_sign_change(is_nose_down, alpha_range[i], alpha_range[i+1],
+            alpha_tol)
+        derivatives = stability_derivatives(solver, body_aero, alpha, beta, wind_speed;
+            backend, throw_on_fail=true)
+        push!(trims, (alpha=alpha, dCMy_dalpha=derivatives.dalpha[5]))
+    end
+    return trims
+end
+
+"""
+    coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
+
+Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` solved at angle of
+attack `alpha` [rad], sideslip `beta` [rad] and `wind_speed` [m/s], at the rotation rate
+`body_aero.omega`. Throws a [`SolveFailure`](@ref) if the solve misses the solver's
+tolerances.
+"""
+function coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
+    set_va!(body_aero, apparent_wind(alpha, beta, wind_speed), body_aero.omega)
+    sol = solve!(solver, body_aero; throw_on_fail=true)
+    return [sol.force_coeffs; sol.moment_coeffs]
+end
+
+"""
+    nose_down(solver, body_aero, alpha, beta, wind_speed)
+
+Whether `CMy` from [`coeffs_at_angles`](@ref) is negative.
+"""
+nose_down(solver, body_aero, alpha, beta, wind_speed) =
+    coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)[5] < 0
+
+"""
+    bisect_sign_change(predicate, low, high, tol)
+
+Bisect `[low, high]`, across which the boolean `predicate` flips, to a width of `tol` and
+return the midpoint.
+"""
+function bisect_sign_change(predicate, low, high, tol)
+    predicate_low = predicate(low)
+    while high - low > tol
+        middle = (low + high) / 2
+        if predicate(middle) == predicate_low
+            low = middle
+        else
+            high = middle
+        end
+    end
+    return (low + high) / 2
+end

test/runtests.jl

@@ -58,6 +58,7 @@ function include_selected_tests()
     should_run_test("solver/test_forwarddiff.jl") && include("solver/test_forwarddiff.jl")
     should_run_test("solver/test_backend_comparison.jl") && include("solver/test_backend_comparison.jl")
     should_run_test("solver/test_unrefined_dist.jl") && include("solver/test_unrefined_dist.jl")
+    should_run_test("solver/test_stability.jl") && include("solver/test_stability.jl")
     should_run_test("verification/test_verification.jl") && include("verification/test_verification.jl")
     should_run_test("VortexStepMethod/test_VortexStepMethod.jl") && include("VortexStepMethod/test_VortexStepMethod.jl")
     should_run_test("wake/test_wake.jl") && include("wake/test_wake.jl")