diff --git a/CHANGELOG.md b/CHANGELOG.md index 50f48fdf..d593e4a7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,28 @@ # Changelog +## Unreleased + +### Added +- `flow_curvature` solver setting (default `false`). When enabled, each section + gets the thin-airfoil pitch-rate moment increment `Δcm = -(π/4) q̂` with + `q̂ = q c / (2 v_rel)` and `q = ω ⋅ y_airf`. A section rotating about its own + spanwise axis sees an incidence that varies linearly along the chord, which is + equivalent to parabolic camber and produces a quarter-chord moment that a + single control point cannot represent. The lift response to `q` was already + exact because the inflow is sampled at the three-quarter-chord point, so only + the moment was missing. +- `pitch_rate_dist` field on `BodyAerodynamics`: each panel's rotation rate about + its own spanwise axis, which is what `flow_curvature` reads. `set_va!(body_aero, + va, omega)` fills it by projecting `omega` onto every panel's `y_airf`, so + panels at different dihedral see different rates from one body rate. The + distributed `set_va!(body_aero, va_distribution; pitch_rate_dist)` takes it + directly, so twist and flapping rates of a deforming wing — which no single + body rate can express — reach the moment. Omitting the keyword zeroes it rather + than reusing a stale `omega`. +- `section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord)` builds + one entry of that distribution from a section's edge velocities, and reduces to + `ω ⋅ y_airf` for rigid motion. + ## VortexStepMethod v4.0.0 2026-08-03 ### Added diff --git a/docs/src/functions.md b/docs/src/functions.md index 72de8257..b5cc0dd4 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -82,6 +82,7 @@ CurrentModule = VortexStepMethod ## Setting the inflow conditions and solving ```@docs set_va! +section_pitch_rate solve solve! reinit!(body_aero::BodyAerodynamics{P, W, T}) where {P, W, T} diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index bec6cb73..c6ba50f9 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -17,6 +17,8 @@ calculate_cl calculate_cd calculate_cm calculate_cd_cm +flow_curvature_cm +set_pitch_rate_dist! calculate_relative_alpha_and_velocity calculate_relative_alpha_and_relative_velocity update_effective_angle_of_attack! diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index c8b8de59..7f707116 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -30,7 +30,7 @@ export ObjWing, Section, Wing, refine!, reinit! export BodyAerodynamics export Solver, VSMSolution, linearize, solve, solve!, solve_base!, calc_forces! export calculate_results -export add_section!, set_va! +export add_section!, set_va!, section_pitch_rate export calculate_projected_area, calculate_span export MVec3 diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 044d4df7..96e4c60b 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -15,6 +15,9 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru - `stall_angle_list`=zeros(Float64, P): stall angle per panel - `alpha_dist::MVector{P, Float64}` = zeros(Float64, P) - `v_a_dist::MVector{P, Float64}` = zeros(Float64, P) +- `pitch_rate_dist::MVector{P, Float64}` = zeros(Float64, P): rotation rate of each + panel about its own spanwise axis, positive nose-up [rad/s]; set by + [set_va!](@ref) and read when the solver has `flow_curvature` enabled - `work_vectors`::NTuple{10, MVec3} = ntuple(_ -> zeros(MVec3), 10) - `AIC::Array{Float64, 3}` = zeros(3, P, P) - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the kite body reference frame [m²] @@ -34,6 +37,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru stall_angle_list::MVector{P, T} = zeros(MVector{P, T}) alpha_dist::MVector{P, T} = zeros(MVector{P, T}) v_a_dist::MVector{P, T} = zeros(MVector{P, T}) + pitch_rate_dist::MVector{P, T} = zeros(MVector{P, T}) work_vectors::NTuple{10, MVector{3, T}} = ntuple(_ -> zeros(MVector{3, T}), 10) AIC::Array{T, 3} = zeros(T, 3, P, P) projected_area::T = one(T) @@ -714,7 +718,53 @@ function compute_panel_center_of_pressures( end """ - calculate_results(body_aero::BodyAerodynamics, gamma_new, + flow_curvature_cm(pitch_rate, chord, v_rel) + +Quarter-chord moment increment of a section rotating about its own spanwise axis, +from thin airfoil theory. The rotation makes the local incidence vary linearly +along the chord, which is equivalent to parabolic camber and yields +`Δcm = -(π/4) q̂` with `q̂ = q c / (2 v_rel)` and `q` positive nose-up. +Independent of the pivot location; the lift response to `q` needs no correction +because the inflow is already sampled at the three-quarter-chord control point. +""" +@inline function flow_curvature_cm(pitch_rate, chord, v_rel) + v_rel > 0 || return zero(chord) + return -0.25π * pitch_rate * chord / (2v_rel) +end + +""" + section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord) + +Rate at which a section rotates about its own spanwise axis, from the velocities +of its leading and trailing edge. Positive nose-up, matching +[`flow_curvature_cm`](@ref). Use this to build a `pitch_rate_dist` for +[`set_va!`](@ref) from a deforming structure, where twist and flapping rates +differ per section and no single body rate describes them. +""" +@inline function section_pitch_rate(velocity_leading, velocity_trailing, + z_airf, chord) + chord > 0 || return zero(chord) + normal_rate = dot3(velocity_trailing, z_airf) - + dot3(velocity_leading, z_airf) + return -normal_rate / chord +end + +""" + set_pitch_rate_dist!(body_aero, omega) + +Fill `body_aero.pitch_rate_dist` from a rigid-body turn rate by projecting it +onto each panel's own spanwise axis. Panels with different dihedral see +different rates from the same `omega`. +""" +function set_pitch_rate_dist!(body_aero::BodyAerodynamics, omega) + for (i, panel) in enumerate(body_aero.panels) + body_aero.pitch_rate_dist[i] = dot3(omega, panel.y_airf) + end + return nothing +end + +""" + calculate_results(body_aero::BodyAerodynamics, gamma_new, density, core_radius_fraction, mu, alpha_dist, v_a_dist, @@ -726,6 +776,9 @@ end Calculate final aerodynamic results. Reference point is in the kite body (KB) frame. +`flow_curvature` adds [`flow_curvature_cm`](@ref) to every section moment, read +from `body_aero.omega`. + Returns: Dict: Results including forces, coefficients and distributions """ @@ -747,6 +800,7 @@ function calculate_results( panels::Vector{<:Panel}, is_only_f_and_gamma_output::Bool; correct_aoa::Bool=false, + flow_curvature::Bool=false, ) n_panels = length(panels) @@ -775,6 +829,10 @@ function calculate_results( cl_array[i] = calculate_cl(panel, alpha_dist[i]) cd_array[i], cm_array[i] = calculate_cd_cm( panel, alpha_dist[i]) + if flow_curvature + cm_array[i] += flow_curvature_cm( + body_aero.pitch_rate_dist[i], chord_array[i], v_a_dist[i]) + end panel_width_array[i] = panel.width va_norm = va_norm_array[i] x_norm = norm3(panel.x_airf) @@ -1058,11 +1116,15 @@ Set velocity array and update wake filaments. - 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] + +`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) body_aero.omega .= omega + set_pitch_rate_dist!(body_aero, omega) if all(iszero, omega) va_distribution .= reshape(va, 1, 3) @@ -1092,9 +1154,28 @@ function set_va!(body_aero::BodyAerodynamics{P, W, T}, va::AbstractVector, omega return nothing end -function set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix) +""" + set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix; + pitch_rate_dist=nothing) + +Set a per-panel inflow distribution. `pitch_rate_dist` gives each panel's rotation +rate about its own spanwise axis [rad/s], positive nose-up; build it with +[`section_pitch_rate`](@ref) when the structure deforms, since twist and flapping +rates differ per section and no single body rate describes them. It is reset to +zero when omitted, because this method takes no `omega` and a stale one would +silently feed the `flow_curvature` moment. +""" +function set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix; + pitch_rate_dist=nothing) size(va_distribution, 1) != length(body_aero.panels) && throw(ArgumentError("Number of rows in va distribution should be equal to number of panels.")) + if isnothing(pitch_rate_dist) + body_aero.pitch_rate_dist .= 0 + else + length(pitch_rate_dist) != length(body_aero.panels) && + throw(ArgumentError("Length of pitch rate distribution should be equal to number of panels.")) + body_aero.pitch_rate_dist .= pitch_rate_dist + end for (i, panel) in enumerate(body_aero.panels) panel.va .= va_distribution[i, :] diff --git a/src/settings.jl b/src/settings.jl index 83b05c90..ebc65ed5 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -89,6 +89,8 @@ Solver configuration, used within [`VSMSettings`](@ref). (default `false`) - `correct_aoa`: Perform angle of attack correction (default `false`) +- `flow_curvature`: Add the thin-airfoil pitch-rate moment increment to each + section (default `false`) """ @with_kw mutable struct SolverSettings n_panels::Int64 = 40 @@ -110,6 +112,7 @@ Solver configuration, used within [`VSMSettings`](@ref). mu::Float64 = 1.81e-5 # dynamic viscosity [N·s/m²] calc_only_f_and_gamma::Bool=false # whether to only output f and gamma correct_aoa::Bool=false # perform aoa correction + flow_curvature::Bool=false # thin-airfoil pitch-rate moment increment end """ diff --git a/src/solver.jl b/src/solver.jl index b3e36d03..ba31e1fd 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -135,6 +135,8 @@ Main solver structure for the Vortex Step Method.See also: [solve](@ref) - `core_radius_fraction`::Float64 = 1e-20: - mu::Float64 = 1.81e-5: Dynamic viscosity [N·s/m²] - `is_only_f_and_gamma_output`::Bool = false: Whether to only output f and gamma +- `flow_curvature`::Bool = false: Add the thin-airfoil pitch-rate moment + increment `-(π/4) q̂` to each section, see: [flow_curvature_cm](@ref) - `reference_point`::MVec3 = [0.0, 0.0, 0.0]: Moment reference point in body frame ## Solution @@ -173,6 +175,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [solve!](@ref) mu::T = T(1.81e-5) is_only_f_and_gamma_output::Bool = false correct_aoa::Bool = false + flow_curvature::Bool = false reference_point::MVector{3, T} = zeros(MVector{3, T}) # Intermediate results @@ -214,6 +217,7 @@ function Solver(body_aero, settings::VSMSettings) mu=ss.mu, is_only_f_and_gamma_output=ss.calc_only_f_and_gamma, correct_aoa=ss.correct_aoa, + flow_curvature=ss.flow_curvature, reference_point=reference_point, ) end @@ -288,6 +292,11 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; for (i, panel) in enumerate(panels) # zero bytes cl_dist[i] = calculate_cl(panel, alpha_dist[i]) cd_dist[i], cm_dist[i] = calculate_cd_cm(panel, alpha_dist[i]) + if solver.flow_curvature + cm_dist[i] += flow_curvature_cm( + body_aero.pitch_rate_dist[i], solver.sol._chord_dist[i], + v_a_dist[i]) + end width_dist[i] = panel.width # Geometric AoA using panel-local axes and prescribed @@ -599,7 +608,8 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n solver.br.va_unit_dist, body_aero.panels, solver.is_only_f_and_gamma_output; - correct_aoa=solver.correct_aoa + correct_aoa=solver.correct_aoa, + flow_curvature=solver.flow_curvature ) # Attach geometric AoA (already computed in calculate_results) to solver.sol if haskey(results, "alpha_geometric") @@ -1218,6 +1228,7 @@ function make_dual_shadow(solver::Solver{P, U, Float64}, mu = TD(solver.mu), is_only_f_and_gamma_output = solver.is_only_f_and_gamma_output, correct_aoa = solver.correct_aoa, + flow_curvature = solver.flow_curvature, reference_point = MVector{3, TD}(solver.reference_point), ) return body_aero_d, solver_d diff --git a/test/runtests.jl b/test/runtests.jl index 6635b965..d613c2b5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -52,6 +52,7 @@ function include_selected_tests() should_run_test("ram_geometry/test_kite_geometry.jl") && include("ram_geometry/test_kite_geometry.jl") should_run_test("settings/test_settings.jl") && include("settings/test_settings.jl") should_run_test("solver/test_solver.jl") && include("solver/test_solver.jl") + should_run_test("solver/test_flow_curvature.jl") && include("solver/test_flow_curvature.jl") should_run_test("solver/test_forwarddiff.jl") && include("solver/test_forwarddiff.jl") should_run_test("solver/test_backend_comparison.jl") && include("solver/test_backend_comparison.jl") should_run_test("solver/test_unrefined_dist.jl") && include("solver/test_unrefined_dist.jl") diff --git a/test/solver/test_flow_curvature.jl b/test/solver/test_flow_curvature.jl new file mode 100644 index 00000000..7a5c6b08 --- /dev/null +++ b/test/solver/test_flow_curvature.jl @@ -0,0 +1,108 @@ +using VortexStepMethod +using VortexStepMethod: flow_curvature_cm +using LinearAlgebra +using Test + +@testset "section_pitch_rate reduces to the rigid-body rate" begin + # a rigid section rotating at q about y_airf must give back exactly q, which + # is what lets one expression serve both rigid and deforming wings + x_airf, y_airf, z_airf = [1.0, 0, 0], [0, 1.0, 0], [0, 0, 1.0] + chord, q = 1.7, 0.8 + omega = q * y_airf + leading = [0.0, 0.0, 0.0] + trailing = chord * x_airf + @test section_pitch_rate(cross(omega, leading), cross(omega, trailing), + z_airf, chord) ≈ q + @test section_pitch_rate([0.0, 0, 0], [0.0, 0, 0], z_airf, chord) == 0.0 + @test section_pitch_rate([0.0, 0, 0], [0.0, 0, 1.0], z_airf, 0.0) == 0.0 +end + +@testset "Flow curvature pitch-rate moment" begin + chord, span, V, q = 1.0, 8.0, 20.0, 0.5 + + wing = Wing(12) + for y in range(span / 2, -span / 2, length=5) + add_section!(wing, [0.0, y, 0.0], [chord, y, 0.0], INVISCID) + end + refine!(wing) + body_aero = BodyAerodynamics([wing]) + set_va!(body_aero, [V, 0.0, 0.0], [0.0, q, 0.0]) + + panel = body_aero.panels[6] + q_local = dot(body_aero.omega, panel.y_airf) + + @testset "increment matches thin airfoil theory" begin + @test body_aero.pitch_rate_dist[6] ≈ q_local + @test flow_curvature_cm(q_local, chord, V) ≈ + -0.25π * q_local * chord / (2V) + @test flow_curvature_cm(0.0, chord, V) == 0.0 + @test flow_curvature_cm(q_local, chord, 0.0) == 0.0 + end + + @testset "positive rate about y_airf raises aft incidence" begin + # nose-up rotation loads the aft chord more, which is what makes the + # negative increment a damping rather than a driving moment + aft = panel.aero_center + 0.1chord * panel.x_airf + dv = cross(body_aero.omega, panel.aero_center) - + cross(body_aero.omega, aft) + @test sign(dot(dv, panel.z_airf)) == sign(q_local) + end + + solver_off = Solver(body_aero; flow_curvature=false) + solver_on = Solver(body_aero; flow_curvature=true) + + function moment_at(solver, omega) + set_va!(body_aero, [V, 0.0, 0.0], omega) + solve!(solver, body_aero) + return copy(solver.sol.moment) + end + + @testset "opposes the rotation" begin + # closed form of the summed increment, pi*rho*V*S*c^2/16 per rad/s, + # exact in the limit where the induced velocity is small against V + expected = π * 1.225 * V * span * chord * chord^2 / 16 * q + for omega in ([0.0, q, 0.0], [0.0, -q, 0.0]) + dM = moment_at(solver_on, omega) .- moment_at(solver_off, omega) + @test dot(dM, omega) < 0 + @test norm(dM) ≈ expected rtol = 0.05 + end + end + + @testset "no effect without rotation" begin + dM = moment_at(solver_on, zeros(3)) .- moment_at(solver_off, zeros(3)) + @test all(iszero, dM) + end + + @testset "defaults to off" begin + @test Solver(body_aero).flow_curvature == false + @test VortexStepMethod.SolverSettings().flow_curvature == false + end + + @testset "distributed rates drive a deformation mode" begin + n = length(body_aero.panels) + va_dist = repeat([V 0.0 0.0], n) + + set_va!(body_aero, va_dist) + @test all(iszero, body_aero.pitch_rate_dist) + solve!(solver_on, body_aero) + base = copy(solver_on.sol.cm_dist) + + # antisymmetric twist rate: no rigid-body omega can express this + rates = [panel.aero_center[2] > 0 ? 1.0 : -1.0 for panel in body_aero.panels] + set_va!(body_aero, va_dist; pitch_rate_dist=rates) + @test body_aero.pitch_rate_dist ≈ rates + solve!(solver_on, body_aero) + + for i in 1:n + @test solver_on.sol.cm_dist[i] - base[i] ≈ + flow_curvature_cm(rates[i], solver_on.sol._chord_dist[i], + solver_on.lr.v_a_dist[i]) + end + # the two half-wings must be driven in opposite senses + @test sign(solver_on.sol.cm_dist[1] - base[1]) == + -sign(solver_on.sol.cm_dist[n] - base[n]) + + @test_throws ArgumentError set_va!(body_aero, va_dist; + pitch_rate_dist=rates[1:end-1]) + end +end