diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c2babbf..f13c68df 100644 --- a/CHANGELOG.md +++ b/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 @@ -20,6 +24,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 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 d7123c47..be9f9156 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 @@ -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 diff --git a/src/solver.jl b/src/solver.jl index d032e697..9f2e5ede 100644 --- a/src/solver.jl +++ b/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, diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index c2c60e05..2fbe4b2c 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/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) diff --git a/test/solver/test_forwarddiff.jl b/test/solver/test_forwarddiff.jl index ee36f4ed..b23fa8b0 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_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)) 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, ), )