stability_derivatives gives roll, pitch and yaw rate derivatives about solver.reference_point

CHANGELOG.md

@@ -5,8 +5,14 @@
 ### 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_vec, 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!`
+  and `linearize` calls until it is given again.
 - Spanwise-flow viscous drag correction (Gaunaa et al. 2024,
   doi:10.1088/1742-6596/2767/2/022068): each section gets a drag increment and a force
   along its span from the flow across it, in `solve!`, `solve` and `linearize`. Opt-in
@@ -21,11 +27,10 @@
 
 ### Fixed
 
+- `set_va!(body_aero, settings)` applies `condition.yaw_rate` as a turn rate about the
+  body z axis; it was read from the settings file and ignored.
 - The `VSMSolution` docstring gives `lift_dist`, `drag_dist` and `panel_moment_dist` in
   the per-unit-span units they hold, [N/m] and [Nm/m], instead of [N] and [Nm].
-
-### Fixed
-
 - With `artificial_damping` on, an iteration whose circulation is already smooth no longer
   re-applies the previous iteration's damping correction.
 

docs/src/settings.md

@@ -32,7 +32,7 @@ condition:
   wind_speed: 10.0                # free-stream velocity magnitude [m/s]
   alpha: 5.0                      # angle of attack [°]
   beta: 0.0                       # sideslip angle [°]
-  yaw_rate: 0.0                   # yaw rate [°/s]
+  yaw_rate: 0.0                   # turn rate about the body z axis [°/s]
 
 wings:
   - name: main_wing               # label the wing carries into plots and output

src/body_aerodynamics.jl

@@ -8,6 +8,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
 - wings::Vector{W}: A vector of wings of type `W <: AbstractWing`; a body can have multiple wings
 - `va::MVec3` = zeros(MVec3):   A vector of the apparent wind speed, see: [`MVec3`](@ref)
 - `omega`::MVec3 = zeros(MVec3): A vector of the turn rates around the kite body axes
+- `reference_point`::MVec3 = zeros(MVec3): The point `omega` turns the body about [m]
 - `gamma_distribution`=zeros(Float64, P): A vector of the circulation
                         of the velocity field; Length: Number of segments. [m²/s]
 - `alpha_uncorrected`=zeros(Float64, P): angles of attack per panel
@@ -35,6 +36,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
     _va::MVector{3, T} = zeros(MVector{3, T})
     has_distributed_va::Bool = false
     omega::MVector{3, T} = zeros(MVector{3, T})
+    reference_point::MVector{3, T} = zeros(MVector{3, T})
     gamma_distribution::MVector{P, T} = zeros(MVector{P, T})
     alpha_uncorrected::MVector{P, T} = zeros(MVector{P, T})
     alpha_corrected::MVector{P, T} = zeros(MVector{P, T})
@@ -155,6 +157,8 @@ function Base.setproperty!(obj::BodyAerodynamics, sym::Symbol, val)
         set_va!(obj, val)
     elseif sym === :omega
         set_va!(obj, obj._va, val)
+    elseif sym === :reference_point
+        set_va!(obj, obj._va, obj.omega; reference_point=val)
     else
         setfield!(obj, sym, val)
     end
@@ -1052,44 +1056,32 @@ end
 
 
 """
-    set_va!(body_aero::BodyAerodynamics, va_vec::VelVector, omega=zeros(MVec3))
+    set_va!(body_aero::BodyAerodynamics, va_vec::VelVector, omega=zeros(MVec3);
+            reference_point=body_aero.reference_point)
 
-Set velocity array and update wake filaments.
+Set a uniform apparent wind and a body turn rate, and update the wake filaments. Each
+panel sees `va_vec - omega × (control_point - reference_point)`.
 
 # Arguments
 - body_aero::BodyAerodynamics: The [`BodyAerodynamics`](@ref) struct to modify
 - `va_vec::VelVector`: Velocity vector of the apparent wind speed       [m/s]
 - `omega::VelVector`: Turn rate vector around x y and z axis            [rad/s]
+- `reference_point`: Point the body turns about, stored on `body_aero`  [m]
 
 `omega` is also projected onto each panel's spanwise axis into
 `pitch_rate_dist`, which the solver reads when `flow_curvature` is enabled.
 """
 function set_va!(body_aero::BodyAerodynamics{P, W, T}, va_vec::AbstractVector,
-                 omega=zeros(MVector{3, T})) where {P, W, T}
-    n_panels = length(body_aero.panels)
-    va_vec_dist = zeros(T, n_panels, 3)
+                 omega=zeros(MVector{3, T});
+                 reference_point=body_aero.reference_point) where {P, W, T}
     body_aero.omega .= omega
+    body_aero.reference_point .= reference_point
     set_pitch_rate_dist!(body_aero, omega)
 
-    if all(iszero, omega)
-        va_vec_dist .= reshape(va_vec, 1, 3)
-    else
-        idx = 1
-        for wing in body_aero.wings
-            panel_end = idx + wing.n_panels - 1
-
-            # Calculate velocities for each panel in this wing slice
-            for j in idx:panel_end
-                omega_va_vec = -omega × body_aero.panels[j].control_point
-                va_vec_dist[j, :] .= omega_va_vec .+ va_vec
-            end
-            idx = panel_end + 1
-        end
-    end
-
-    # Update panel velocities
+    va_vec_dist = zeros(T, P, 3)
     for (i, panel) in enumerate(body_aero.panels)
-        panel.va .= va_vec_dist[i,:]
+        panel.va .= va_vec .- omega × (panel.control_point .- body_aero.reference_point)
+        va_vec_dist[i, :] .= panel.va
     end
 
     # Update wake elements
@@ -1147,7 +1139,8 @@ apparent_wind(alpha, beta, wind_speed) =
     set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
 
 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`.
+`beta` [°] and `wind_speed` [m/s] of `settings.condition`, and turn the body at its
+`yaw_rate` [°/s] about Z_b through `body_aero.reference_point`.
 
 # Example
 ```julia
@@ -1160,5 +1153,6 @@ function set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
     condition = settings.condition
     va_vec = apparent_wind(deg2rad(condition.alpha), deg2rad(condition.beta),
         condition.wind_speed)
-    set_va!(body_aero, va_vec)
+    omega = [0.0, 0.0, deg2rad(condition.yaw_rate)]
+    set_va!(body_aero, va_vec, omega)
 end

src/solver.jl

@@ -1216,10 +1216,9 @@ function make_dual_shadow(solver::Solver{P, U, Float64},
     length(body_aero.wings) == 1 || throw(ArgumentError(
         "make_dual_shadow currently supports body_aero with one wing"))
     wing_d = _wing_with_eltype(body_aero.wings[1], TD)
-    body_aero_d = BodyAerodynamics([wing_d];
-        va = MVector{3, TD}(body_aero._va),
-        omega = MVector{3, TD}(body_aero.omega),
-    )
+    body_aero_d = BodyAerodynamics([wing_d])
+    set_va!(body_aero_d, MVector{3, TD}(body_aero._va), MVector{3, TD}(body_aero.omega);
+            reference_point=body_aero.reference_point)
     solver_d = Solver(body_aero_d;
         solver_type = solver.solver_type,
         aerodynamic_model_type = solver.aerodynamic_model_type,

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_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...)
+    set_va!(body_aero, va_vec, body_aero.omega; reference_point=solver.reference_point)
+    jac, results, converged = linearize(solver, body_aero, [va_vec; 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
 
 """
@@ -55,11 +62,12 @@ end
 
 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.
+`body_aero.omega` about `solver.reference_point`. 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)
+    set_va!(body_aero, apparent_wind(alpha, beta, wind_speed), body_aero.omega;
+        reference_point=solver.reference_point)
     sol = solve!(solver, body_aero; throw_on_fail=true)
     return [sol.force_coeffs; sol.moment_coeffs]
 end

test/body_aerodynamics/test_body_aerodynamics.jl

@@ -431,9 +431,10 @@ end
     @test length(results_NEW["cd_distribution"]) == length(body_aero.panels)
 end
 
-@testset "set_va! with VSMSettings" begin
+@testset "set_va! with VSMSettings applies the yaw rate about body z" begin
     settings_file = create_temp_wing_settings("body_aerodynamics", "test_wing.yaml";
-                                              alpha=10.0, beta=5.0, wind_speed=15.0)
+                                              alpha=10.0, beta=5.0, wind_speed=15.0,
+                                              yaw_rate=30.0)
     try
         settings   = VSMSettings(settings_file)
         wing       = Wing(settings)
@@ -444,11 +445,13 @@ end
 
         α, β, wind_speed = deg2rad(10.0), deg2rad(5.0), 15.0
         expected_va_vec = wind_speed .* [cos(α)*cos(β), sin(β), sin(α)*cos(β)]
+        omega = [0.0, 0.0, deg2rad(30.0)]
 
         for p in body_aero.panels
-            @test p.va ≈ expected_va_vec atol=1e-10
+            @test p.va ≈ expected_va_vec .- omega × p.control_point atol=1e-10
         end
         @test body_aero._va ≈ expected_va_vec atol=1e-10
+        @test body_aero.omega ≈ omega
     finally
         isfile(settings_file) && rm(settings_file; force=true)
     end
@@ -504,6 +507,40 @@ end
     @test body_aero.omega ≈ new_omega
 end
 
+"""
+    test_rigid_body_inflow(body_aero, va_vec, omega, reference_point)
+
+Test that every panel sees `va_vec` plus the inflow of a body turning at `omega` about
+`reference_point`.
+"""
+function test_rigid_body_inflow(body_aero, va_vec, omega, reference_point)
+    for panel in body_aero.panels
+        expected_va_vec = va_vec .- omega × (panel.control_point .- reference_point)
+        @test panel.va ≈ expected_va_vec atol=1e-12
+    end
+end
+
+@testset "set_va! rotates the body about reference_point" begin
+    body_aero = BodyAerodynamics([inviscid_wing([0.0, 1.0, 2.0]),
+                                  inviscid_wing([10.0, 11.0, 12.0])])
+    va_vec = [10.0, 0.0, 1.0]
+    omega = [0.1, 0.2, 1.0]
+    reference_point = [0.25, 6.0, -0.5]
+
+    set_va!(body_aero, va_vec, omega; reference_point)
+    @test body_aero.reference_point ≈ reference_point
+    test_rigid_body_inflow(body_aero, va_vec, omega, reference_point)
+
+    body_aero.omega = 2 .* omega
+    test_rigid_body_inflow(body_aero, va_vec, 2 .* omega, reference_point)
+
+    reinit!(body_aero; va=va_vec, omega)
+    test_rigid_body_inflow(body_aero, va_vec, omega, reference_point)
+
+    body_aero.reference_point = zeros(3)
+    test_rigid_body_inflow(body_aero, va_vec, omega, zeros(3))
+end
+
 """
     solve_wings(wings)
 

test/solver/test_forwarddiff.jl

@@ -21,19 +21,24 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs,
     omega = [0.0, 0.0, 0.0]
     y0 = [va_vec; omega]
 
-    @testset "AutoForwardDiff matches AutoFiniteDiff (LOOP, INVISCID)" begin
-        solver = Solver(body_aero;
+    turns = ((omega, zeros(3)), ([0.0, 0.0, 0.2], [0.5, 4.0, 0.0]))
+    @testset "ForwardDiff matches FiniteDiff about $reference_point (LOOP, INVISCID)" for
+            (omega_op, reference_point) in turns
+        pivot_body = BodyAerodynamics([wing])
+        set_va!(pivot_body, va_vec, omega_op; reference_point)
+        y_op = [va_vec; omega_op]
+        solver = Solver(pivot_body;
             use_gamma_prev=false,
             type_initial_gamma_distribution=ELLIPTIC)
 
         jac_fwd, _, fwd_converged = VortexStepMethod.linearize(
-            solver, body_aero, y0;
+            solver, pivot_body, y_op;
             theta_idxs=nothing, va_idxs=1:3, omega_idxs=4:6,
             aero_coeffs=true, backend=AutoForwardDiff())
         @test fwd_converged
 
         jac_fd, _, fd_converged = VortexStepMethod.linearize(
-            solver, body_aero, y0;
+            solver, pivot_body, y_op;
             theta_idxs=nothing, va_idxs=1:3, omega_idxs=4:6,
             aero_coeffs=true,
             backend=AutoFiniteDiff(absstep=1e-5, relstep=1e-5))

test/solver/test_stability.jl

@@ -1,5 +1,5 @@
 using VortexStepMethod
-using VortexStepMethod: coeffs_at_angles
+using VortexStepMethod: apparent_wind, coeffs_at_angles
 using Test
 
 """
@@ -43,6 +43,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_vec = apparent_wind(alpha, beta, wind_speed)
+    omega = [0.1, -0.05, 0.08]
+    set_va!(body_aero, va_vec, 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_vec, 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
 
@@ -68,6 +99,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 "a NONLIN solver with backend=nothing finds the same trim" begin
         body_aero = trimmable_wing_aero(0.05)
         trim_loop = only(trim_angle(Solver(body_aero), body_aero, beta, wind_speed))

test/test_data_utils.jl

@@ -135,6 +135,7 @@ function create_temp_wing_settings(module_name, wing_file;
     alpha=10.0,
     beta=5.0,
     wind_speed=15.0,
+    yaw_rate=0.0,
 )
     wing_file_path = isabspath(wing_file) ? wing_file : test_data_path(module_name, wing_file)
     wing_file_path = replace(normpath(wing_file_path), '\\' => '/')
@@ -158,6 +159,7 @@ function create_temp_wing_settings(module_name, wing_file;
             "alpha" => alpha,
             "beta" => beta,
             "wind_speed" => wind_speed,
+            "yaw_rate" => yaw_rate,
         ),
     )