From 9af8d69c0904448896875080774d779cf187c486 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:04:48 +0200 Subject: [PATCH] Turn the body about a reference point in set_va!, and apply yaw_rate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit set_va!(body_aero, va, omega; reference_point) gives each panel va - omega × (control_point - reference_point). The point is stored on BodyAerodynamics so reinit!, the omega setter and the ForwardDiff shadow in linearize keep it; it starts at the origin. set_va!(body_aero, settings) passes condition.yaw_rate as omega about body z instead of dropping it. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 9 ++++ docs/src/settings.md | 2 +- src/body_aerodynamics.jl | 46 +++++++++---------- src/solver.jl | 7 ++- .../test_body_aerodynamics.jl | 43 +++++++++++++++-- test/solver/test_forwarddiff.jl | 13 ++++-- test/test_data_utils.jl | 2 + 7 files changed, 85 insertions(+), 37 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4c00c754..d2e36782 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,13 @@ ## Unreleased +### Added + +- `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!` + and `linearize` calls until it is given again. + ### Changed - Requires Julia 1.12 or 1.13; 1.10 and 1.11 keep resolving v5.1.1. @@ -11,6 +18,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]. diff --git a/docs/src/settings.md b/docs/src/settings.md index 3b2fba46..90004454 100644 --- a/docs/src/settings.md +++ b/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 diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index d04a966d..453eef2f 100644 --- a/src/body_aerodynamics.jl +++ b/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 @@ -1066,43 +1070,32 @@ end """ - set_va!(body_aero::BodyAerodynamics, va::VelVector, omega=zeros(MVec3)) + set_va!(body_aero::BodyAerodynamics, va::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 - omega × (control_point - reference_point)`. # Arguments - body_aero::BodyAerodynamics: The [`BodyAerodynamics`](@ref) struct to modify - `va::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::AbstractVector, omega=zeros(MVector{3, T})) where {P, W, T} - n_panels = length(body_aero.panels) - va_distribution = zeros(T, n_panels, 3) +function set_va!(body_aero::BodyAerodynamics{P, W, T}, va::AbstractVector, + 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_distribution .= reshape(va, 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 = -omega × body_aero.panels[j].control_point - va_distribution[j, :] .= omega_va .+ va - end - idx = panel_end + 1 - end - end - - # Update panel velocities + va_distribution = zeros(T, P, 3) for (i, panel) in enumerate(body_aero.panels) - panel.va .= va_distribution[i,:] + panel.va .= va .- omega × (panel.control_point .- body_aero.reference_point) + va_distribution[i, :] .= panel.va end # Update wake elements @@ -1156,6 +1149,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(β) @@ -1183,6 +1178,7 @@ function set_va!(body_aero::BodyAerodynamics, settings::VSMSettings) sin(β), # Y_b (right) sin(α)*cos(β) # Z_b (down) ] - - set_va!(body_aero, va) + omega = [0.0, 0.0, deg2rad(settings.condition.yaw_rate)] + + set_va!(body_aero, va, omega) end diff --git a/src/solver.jl b/src/solver.jl index c63d44cc..57c14146 100644 --- a/src/solver.jl +++ b/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, diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index 037a1b6f..b78efdaa 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/test/body_aerodynamics/test_body_aerodynamics.jl @@ -429,9 +429,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) @@ -442,11 +443,13 @@ end α, β, wind_speed = deg2rad(10.0), deg2rad(5.0), 15.0 expected_va = 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 atol=1e-10 + @test p.va ≈ expected_va .- omega × p.control_point atol=1e-10 end @test body_aero._va ≈ expected_va atol=1e-10 + @test body_aero.omega ≈ omega finally isfile(settings_file) && rm(settings_file; force=true) end @@ -502,6 +505,40 @@ end @test body_aero.omega ≈ new_omega end +""" + test_rigid_body_inflow(body_aero, va, omega, reference_point) + +Test that every panel sees `va` plus the inflow of a body turning at `omega` about +`reference_point`. +""" +function test_rigid_body_inflow(body_aero, va, omega, reference_point) + for panel in body_aero.panels + expected_va = va .- omega × (panel.control_point .- reference_point) + @test panel.va ≈ expected_va 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 = [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, omega; reference_point) + @test body_aero.reference_point ≈ reference_point + test_rigid_body_inflow(body_aero, va, omega, reference_point) + + body_aero.omega = 2 .* omega + test_rigid_body_inflow(body_aero, va, 2 .* omega, reference_point) + + reinit!(body_aero; va, omega) + test_rigid_body_inflow(body_aero, va, omega, reference_point) + + body_aero.reference_point = zeros(3) + test_rigid_body_inflow(body_aero, va, omega, zeros(3)) +end + """ solve_wings(wings) diff --git a/test/solver/test_forwarddiff.jl b/test/solver/test_forwarddiff.jl index c052aae4..8f790a58 100644 --- a/test/solver/test_forwarddiff.jl +++ b/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; 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, omega_op; reference_point) + y_op = [va; 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)) diff --git a/test/test_data_utils.jl b/test/test_data_utils.jl index 82acfb37..268c3bb5 100644 --- a/test/test_data_utils.jl +++ b/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, ), )