Merge pull request #353 from OpenSourceAWE/agent/329-set-va-rotates-the-body-about-the-origin

set_va! turns the body about a stored reference_point, and set_va!(body_aero, settings) applies yaw_rate

Author: 1-Bart-1

CHANGELOG.md

@@ -4,6 +4,10 @@
 
 ### Added
 
+- `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
@@ -23,6 +27,8 @@
 
 ### 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].
 - Inside its vortex core, `velocity_3D_trailing_vortex!` induces an azimuthal velocity

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
@@ -1083,44 +1087,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
@@ -1175,6 +1167,8 @@ 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`
 
 The velocity vector is constructed as:
 - X_b (forward): wind_speed * cos(α) * cos(β)  
@@ -1202,6 +1196,7 @@ function set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
         sin(β),         # Y_b (right)
         sin(α)*cos(β)   # Z_b (down)
     ]
-    
-    set_va!(body_aero, va_vec)
+    omega = [0.0, 0.0, deg2rad(settings.condition.yaw_rate)]
+
+    set_va!(body_aero, va_vec, omega)
 end

src/solver.jl

@@ -1215,10 +1215,9 @@ function make_dual_shadow(solver::Solver{P, U, Float64},
                           body_aero::BodyAerodynamics{P, W, Float64},
                           ::Type{TD}) where {P, U, W, TD}
     wings_d = [_wing_with_eltype(wing, TD) for wing in body_aero.wings]
-    body_aero_d = BodyAerodynamics(wings_d;
-        va = MVector{3, TD}(body_aero._va),
-        omega = MVector{3, TD}(body_aero.omega),
-    )
+    body_aero_d = BodyAerodynamics(wings_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,

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