From 34b10a23af8d55efc3aa52e167f3febb48ef6205 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Thu, 17 Sep 2026 01:31:01 +0200 Subject: [PATCH 1/2] Add roll, pitch and yaw rate derivatives to stability_derivatives MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The omega columns of linearize, scaled to p̂ = pb/2V, q̂ = q c_ref/2V and r̂ = rb/2V, with the body turning about solver.reference_point. trim_angle's bracketing solves turn about the same point. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 6 +++-- src/stability.jl | 30 +++++++++++++++--------- test/solver/test_stability.jl | 44 +++++++++++++++++++++++++++++++++++ 3 files changed, 67 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cbf09ee6..57b581ea 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,8 +5,10 @@ ### 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. + with respect to angle of attack, sideslip and the nondimensional roll, pitch and yaw + rates p̂ = pb/2V, q̂ = q c_ref/2V, r̂ = rb/2V, turning about `solver.reference_point`, + and `trim_angle` the angles of attack at which `CMy` changes sign, with the slope that + says whether each trim is stable. - `set_va!(body_aero, va, omega; reference_point)` turns the body about `reference_point` [m] instead of the origin. The point is stored on `BodyAerodynamics`, starts at the origin, and is kept by later `set_va!`, `reinit!` diff --git a/src/stability.jl b/src/stability.jl index 8046e935..aabe3b8a 100644 --- a/src/stability.jl +++ b/src/stability.jl @@ -2,25 +2,32 @@ 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. +`alpha` [rad], sideslip `beta` [rad], `wind_speed` [m/s] and rotation rate `body_aero.omega` +[rad/s], and their derivatives: `dalpha` and `dbeta` [1/rad], and `dp`, `dq`, `dr` with +respect to the body rates about x, y and z as p̂ = pb/2V, q̂ = q c_ref/2V and r̂ = rb/2V, +with b the wing span, c_ref `body_aero.c_ref` and V `wind_speed`. Moments are about, and +the body turns about, `solver.reference_point`. `kwargs` go to [`linearize`](@ref), which +leaves `body_aero` at this inflow. -Returns `(coeffs, dalpha, dbeta, converged)`. +Returns `(coeffs, dalpha, dbeta, dp, dq, dr, converged)`. """ function stability_derivatives(solver::Solver, body_aero::BodyAerodynamics, alpha, beta, wind_speed; kwargs...) va = apparent_wind(alpha, beta, wind_speed) - jac, results, converged = linearize(solver, body_aero, va; - theta_idxs=nothing, va_idxs=1:3, aero_coeffs=true, kwargs...) + set_va!(body_aero, va, body_aero.omega; reference_point=solver.reference_point) + jac, results, converged = linearize(solver, body_aero, [va; body_aero.omega]; + theta_idxs=nothing, va_idxs=1:3, omega_idxs=4:6, 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) + span = body_aero.wings[1].span + return (coeffs=results[1:6], dalpha=coeff_jac[:, 1:3] * dva_dalpha, + dbeta=coeff_jac[:, 1:3] * dva_dbeta, + dp=coeff_jac[:, 4] * 2wind_speed / span, + dq=coeff_jac[:, 5] * 2wind_speed / body_aero.c_ref, + dr=coeff_jac[:, 6] * 2wind_speed / span, converged) end """ @@ -54,10 +61,11 @@ end pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed) `CMy` 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`. +`wind_speed` [m/s], at the rotation rate `body_aero.omega` about `solver.reference_point`. """ function pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed) - set_va!(body_aero, apparent_wind(alpha, beta, wind_speed), body_aero.omega) + set_va!(body_aero, apparent_wind(alpha, beta, wind_speed), body_aero.omega; + reference_point=solver.reference_point) return solve!(solver, body_aero).moment_coeffs[2] end diff --git a/test/solver/test_stability.jl b/test/solver/test_stability.jl index 07634e72..fbe5898e 100644 --- a/test/solver/test_stability.jl +++ b/test/solver/test_stability.jl @@ -48,6 +48,37 @@ end @test derivatives.dbeta ≈ dbeta rtol = 1e-4 atol = 1e-6 end +@testset "rate derivatives match central differences of solve! about reference_point" begin + body_aero = trimmable_wing_aero(0.05) + reference_point = [0.25, 0.5, 0.1] + solver = Solver(body_aero; reference_point, use_gamma_prev=false) + alpha, beta, wind_speed, step = deg2rad(4.0), deg2rad(3.0), 20.0, 1e-4 + va = apparent_wind(alpha, beta, wind_speed) + omega = [0.1, -0.05, 0.08] + set_va!(body_aero, va, omega) + + derivatives = stability_derivatives(solver, body_aero, alpha, beta, wind_speed) + @test derivatives.converged + @test body_aero.reference_point == reference_point + + function coeffs_at_rate(rate) + set_va!(body_aero, va, rate; reference_point) + sol = solve!(solver, body_aero) + return [sol.force_coeffs; sol.moment_coeffs] + end + rate_scales = 2wind_speed ./ [body_aero.wings[1].span, body_aero.c_ref, + body_aero.wings[1].span] + rate_derivatives = map(1:3) do axis + rate_step = step .* (1:3 .== axis) + (coeffs_at_rate(omega + rate_step) - coeffs_at_rate(omega - rate_step)) / + 2step * rate_scales[axis] + end + @test all(!iszero, rate_derivatives) + @test derivatives.dp ≈ rate_derivatives[1] rtol = 1e-4 atol = 1e-6 + @test derivatives.dq ≈ rate_derivatives[2] rtol = 1e-4 atol = 1e-6 + @test derivatives.dr ≈ rate_derivatives[3] rtol = 1e-4 atol = 1e-6 +end + @testset "trim_angle finds where CMy changes sign" begin beta, wind_speed = 0.0, 20.0 @@ -71,6 +102,19 @@ end @test trim.dCMy_dalpha > 0 end + @testset "pitching body: trim turning about the reference point" begin + body_aero = trimmable_wing_aero(-0.05) + reference_point = [1.0, 0.0, 0.0] + solver = Solver(body_aero; reference_point) + omega = [0.0, 0.5, 0.0] + set_va!(body_aero, apparent_wind(0.0, beta, wind_speed), omega) + trim = only(trim_angle(solver, body_aero, beta, wind_speed)) + set_va!(body_aero, apparent_wind(trim.alpha, beta, wind_speed), omega; + reference_point) + cmy = solve!(solver, body_aero).moment_coeffs[2] + @test abs(cmy) < 1e-5 * abs(trim.dCMy_dalpha) + end + @testset "no sign change in alpha_range: no trim" begin body_aero = trimmable_wing_aero(0.05) solver = Solver(body_aero) From 953b2d1a9dce02dc7e44c89e128b9a9b52753fcf Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 16:01:26 +0200 Subject: [PATCH 2/2] Say where stability_derivatives leaves the pivot and which V scales the rates MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The docstrings of stability_derivatives, trim_angle and coeffs_at_angles now say they store solver.reference_point on body_aero, and that V in p̂, q̂, r̂ is va rather than the panel inflow speed the coefficients are normalised by. The rate test steps body_aero.omega and reuses coeffs_at_angles instead of its own closure. Co-Authored-By: Claude Opus 5 --- src/stability.jl | 17 ++++++++++------- test/solver/test_stability.jl | 23 +++++++++-------------- 2 files changed, 19 insertions(+), 21 deletions(-) diff --git a/src/stability.jl b/src/stability.jl index b239a0d3..7c0d8658 100644 --- a/src/stability.jl +++ b/src/stability.jl @@ -4,10 +4,11 @@ Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` at angle of attack `alpha` [rad], sideslip `beta` [rad], apparent wind speed `va` [m/s] and rotation rate `body_aero.omega` [rad/s], and their derivatives: `dalpha` and `dbeta` [1/rad], and `dp`, -`dq`, `dr` with respect to the body rates about x, y and z as p̂ = pb/2V, q̂ = q c_ref/2V -and r̂ = rb/2V, with b the wing span, c_ref `body_aero.c_ref` and V `va`. Moments are -about, and the body turns about, `solver.reference_point`. `kwargs` go to -[`linearize`](@ref), which leaves `body_aero` at this inflow. +`dq`, `dr` with respect to p̂ = pb/2V, q̂ = q c_ref/2V and r̂ = rb/2V, where p, q, r are the +rates about body x, y, z, b the wing span, c_ref `body_aero.c_ref` and V is `va`, not the +area-weighted panel inflow speed the coefficients are normalised by. Moments are about, +and the body turns about, `solver.reference_point`. `kwargs` go to [`linearize`](@ref). +Leaves `body_aero` at this inflow with `body_aero.reference_point = solver.reference_point`. Returns `(coeffs, dalpha, dbeta, dp, dq, dr, converged)`. """ @@ -39,7 +40,8 @@ sign between neighbouring entries of `alpha_range`, bisected to `alpha_tol` [rad sideslip `beta` [rad] and apparent wind speed `va` [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. +[`SolveFailure`](@ref) if a solve misses the solver's tolerances. Leaves +`body_aero.reference_point = solver.reference_point`. """ function trim_angle(solver::Solver, body_aero::BodyAerodynamics, beta, va; alpha_range=deg2rad.(-5:2:15), alpha_tol=1e-5, backend=AutoForwardDiff()) @@ -62,8 +64,9 @@ end Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` solved at angle of attack `alpha` [rad], sideslip `beta` [rad] and apparent wind speed `va` [m/s], at the -rotation rate `body_aero.omega` about `solver.reference_point`. Throws a -[`SolveFailure`](@ref) if the solve misses the solver's tolerances. +rotation rate `body_aero.omega`, which it turns about and stores as +`body_aero.reference_point`. Throws a [`SolveFailure`](@ref) if the solve misses the +solver's tolerances. """ function coeffs_at_angles(solver, body_aero, alpha, beta, va) set_va!(body_aero, apparent_wind(alpha, beta, va), body_aero.omega; diff --git a/test/solver/test_stability.jl b/test/solver/test_stability.jl index 8df2f2e8..758b68a1 100644 --- a/test/solver/test_stability.jl +++ b/test/solver/test_stability.jl @@ -56,29 +56,24 @@ end reference_point = [0.25, 0.5, 0.1] body_aero, solver = trimmable_wing(0.05; reference_point, use_gamma_prev=false) alpha, beta, va, step = deg2rad(4.0), deg2rad(3.0), 20.0, 1e-4 - va_vec = apparent_wind(alpha, beta, va) omega = [0.1, -0.05, 0.08] - set_va!(body_aero, va_vec, omega) + body_aero.omega = omega derivatives = stability_derivatives(solver, body_aero, alpha, beta, va) @test derivatives.converged @test body_aero.reference_point == reference_point - function coeffs_at_rate(rate) - set_va!(body_aero, va_vec, rate; reference_point) - sol = solve!(solver, body_aero) - return [sol.force_coeffs; sol.moment_coeffs] - end rate_scales = 2va ./ [body_aero.wings[1].span, body_aero.c_ref, body_aero.wings[1].span] - rate_derivatives = map(1:3) do axis + for (axis, derivative) in enumerate((derivatives.dp, derivatives.dq, derivatives.dr)) rate_step = step .* (1:3 .== axis) - (coeffs_at_rate(omega + rate_step) - coeffs_at_rate(omega - rate_step)) / - 2step * rate_scales[axis] + body_aero.omega = omega + rate_step + coeffs_plus = coeffs_at_angles(solver, body_aero, alpha, beta, va) + body_aero.omega = omega - rate_step + coeffs_minus = coeffs_at_angles(solver, body_aero, alpha, beta, va) + central_difference = (coeffs_plus - coeffs_minus) / 2step * rate_scales[axis] + @test !iszero(central_difference) + @test derivative ≈ central_difference rtol = 1e-4 atol = 1e-6 end - @test all(!iszero, rate_derivatives) - @test derivatives.dp ≈ rate_derivatives[1] rtol = 1e-4 atol = 1e-6 - @test derivatives.dq ≈ rate_derivatives[2] rtol = 1e-4 atol = 1e-6 - @test derivatives.dr ≈ rate_derivatives[3] rtol = 1e-4 atol = 1e-6 end @testset "trim_angle finds where CMy changes sign" begin