diff --git a/CHANGELOG.md b/CHANGELOG.md index e08187bc..fb6baad7 100644 --- a/CHANGELOG.md +++ b/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. 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 8a64ce89..d3cf4e2f 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 @@ -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 diff --git a/src/solver.jl b/src/solver.jl index 84cb5297..bf7cf596 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/src/stability.jl b/src/stability.jl index 9ca803cf..c2da2694 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_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 diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index 28121e1d..5943aaba 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 fae465f6..b5014e65 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/solver/test_stability.jl b/test/solver/test_stability.jl index 3e26529f..ac172631 100644 --- a/test/solver/test_stability.jl +++ b/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)) 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, ), )