diff --git a/CHANGELOG.md b/CHANGELOG.md index dd9c7964..cf3905c0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,10 @@ the right tip, z = x × y up. The separate `KB` definition, whose `Z = Y × X` pointed down, is gone; the solver's frame is unchanged. CL, CD and CS are documented as wind-axis coefficients, `cfx`, `cfy` and `cfz` as the body-axis ones. +- The angle of attack is measured in the plane square to the quarter-chord line, from + the chord's projection on that plane (Gaunaa, Li & Pirrung 2026, CP1), everywhere it + is formed and where it turns the lift direction. Swept panels read their polars at a + different angle than before; unswept panels are unchanged. ## VortexStepMethod v6.0.0 2026-09-23 diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index f7482533..a498dff0 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -31,7 +31,6 @@ set_pitch_rate_dist! coeffs_at_angles nose_down bisect_sign_change -calculate_relative_alpha_and_velocity calculate_relative_alpha_and_relative_velocity update_effective_angle_of_attack! calculate_stall_angle_list @@ -41,6 +40,7 @@ calculate_circulation_distribution_elliptical_wing _compute_reference_velocity_from_distribution smooth_distribution! make_dual_shadow +matrix_row ``` ### Panel aerodynamics @@ -54,6 +54,7 @@ smooth_norm panel_span_vector panel_chord_weight panel_axes +inflow_angle effective_alpha panel_inflow dynamic_pressure diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index bd390826..e37caa70 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -28,7 +28,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru influence coefficients, used only for the corrected angle of attack - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the KA body frame [m²] - `c_ref::Float64` = 1.0: Reference chord length (max panel chord) [m] -- `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache() for _ in 1:15] +- `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache() for _ in 1:12] """ @with_kw mutable struct BodyAerodynamics{P, W<:AbstractWing, T, PN<:Panel{T}} panels::Vector{PN} @@ -49,7 +49,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru AIC_aero_center::Array{T, 3} = zeros(T, P, P, 3) projected_area::T = one(T) c_ref::T = one(T) - cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:15] + cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:12] end """ @@ -523,20 +523,16 @@ end """ update_effective_angle_of_attack!(alpha_corrected, body_aero::BodyAerodynamics, gamma, - core_radius_fraction, z_airf_dist, x_airf_dist, - va_vec_dist, va_dist, va_unit_dist) + core_radius_fraction, va_vec_dist, va_dist, + va_unit_dist) -Update angle of attack at aerodynamic center for VSM method. - -Returns: - nothing +Write each panel's [`inflow_angle`](@ref) at its aerodynamic centre into +`alpha_corrected`, the induced velocity taken from `gamma` on the LLT influence matrix. """ function update_effective_angle_of_attack!(alpha_corrected, - body_aero::BodyAerodynamics, + body_aero::BodyAerodynamics, gamma, core_radius_fraction, - z_airf_dist, - x_airf_dist, va_vec_dist, va_dist, va_unit_dist) @@ -551,31 +547,11 @@ function update_effective_angle_of_attack!(alpha_corrected, mul!(view(induced_velocity, :, k), view(body_aero.AIC_aero_center, :, :, k), gamma) end - # In-place relative velocity calculation - relative_velocity = body_aero.cache[2][va_vec_dist] - relative_velocity .= va_vec_dist .+ induced_velocity - - # Preallocate and compute dot products manually - n = size(relative_velocity, 1) - v_normal = body_aero.cache[3][relative_velocity] - v_tangential = body_aero.cache[4][relative_velocity] - - @inbounds for i in 1:n - vn = 0.0 - vt = 0.0 - for j in 1:3 - vn += z_airf_dist[i, j] * relative_velocity[i, j] - vt += x_airf_dist[i, j] * relative_velocity[i, j] - end - v_normal[i] = vn - v_tangential[i] = vt - end - - # Direct angle calculation without temporary arrays - @inbounds for i in 1:n - alpha_corrected[i] = atan(v_normal[i], v_tangential[i]) + @inbounds for (i, panel) in enumerate(body_aero.panels) + v_rel = matrix_row(va_vec_dist, i) .+ matrix_row(induced_velocity, i) + alpha_corrected[i] = inflow_angle(v_rel, SVector{3}(panel.y_airf), + SVector{3}(panel.z_airf)) end - nothing end @@ -832,21 +808,21 @@ function calculate_results( ) n_panels = length(panels) - if length(body_aero.cache) < 15 - append!(body_aero.cache, [LazyBufferCache() for _ in 1:(15 - length(body_aero.cache))]) + if length(body_aero.cache) < 12 + append!(body_aero.cache, [LazyBufferCache() for _ in 1:(12 - length(body_aero.cache))]) end - cl_dist = body_aero.cache[5][alpha_dist] - cd_dist = body_aero.cache[6][alpha_dist] - cm_dist = body_aero.cache[7][alpha_dist] - panel_width_dist = body_aero.cache[8][alpha_dist] - alpha_corrected = body_aero.cache[9][alpha_dist] - cl_prescribed_va = body_aero.cache[10][alpha_dist] - cd_prescribed_va = body_aero.cache[11][alpha_dist] - cs_prescribed_va = body_aero.cache[12][alpha_dist] - f_body_3D = body_aero.cache[13][alpha_dist, (3, length(alpha_dist))] - m_body_3D = body_aero.cache[14][alpha_dist, (3, length(alpha_dist))] - alpha_geometric = body_aero.cache[15][alpha_dist] + cl_dist = body_aero.cache[2][alpha_dist] + cd_dist = body_aero.cache[3][alpha_dist] + cm_dist = body_aero.cache[4][alpha_dist] + panel_width_dist = body_aero.cache[5][alpha_dist] + alpha_corrected = body_aero.cache[6][alpha_dist] + cl_prescribed_va = body_aero.cache[7][alpha_dist] + cd_prescribed_va = body_aero.cache[8][alpha_dist] + cs_prescribed_va = body_aero.cache[9][alpha_dist] + f_body_3D = body_aero.cache[10][alpha_dist, (3, length(alpha_dist))] + m_body_3D = body_aero.cache[11][alpha_dist, (3, length(alpha_dist))] + alpha_geometric = body_aero.cache[12][alpha_dist] fill!(f_body_3D, 0.0) fill!(m_body_3D, 0.0) @@ -861,19 +837,8 @@ function calculate_results( body_aero.pitch_rate_dist[i], chord_dist[i], v_rel_dist[i]) end panel_width_dist[i] = panel.width - va = va_dist[i] - x_norm = norm3(panel.x_airf) - z_norm = norm3(panel.z_airf) - if va == 0.0 || x_norm == 0.0 || z_norm == 0.0 - alpha_geometric[i] = NaN - else - inv_va = 1.0 / va - v_tangential = -dot3(panel.x_airf, panel.va_vec) * - inv_va / x_norm - v_normal = -dot3(panel.z_airf, panel.va_vec) * - inv_va / z_norm - alpha_geometric[i] = atan(-v_normal, -v_tangential) - end + alpha_geometric[i] = iszero(va_dist[i]) ? NaN : inflow_angle(panel.va_vec, + SVector{3}(panel.y_airf), SVector{3}(panel.z_airf)) end # Calculate alpha corrections based on model type @@ -883,8 +848,6 @@ function calculate_results( body_aero, gamma_new, core_radius_fraction, - z_airf_dist, - x_airf_dist, va_vec_dist, va_dist, va_unit_dist diff --git a/src/panel.jl b/src/panel.jl index 3a211857..531ee1be 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -479,16 +479,6 @@ function calculate_relative_alpha_and_relative_velocity( return flow.alpha, flow.v_eff end -""" - calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity) - -Calculate relative angle of attack and relative velocity of the panel. -""" -function calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity) - flow = panel_inflow(panel_axes(panel), panel.va_vec, panel.va_vec, induced_velocity) - return flow.alpha, flow.v_eff -end - """ calculate_cl(panel::Panel, alpha) calculate_cl(panel::Panel, alpha, delta) diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl index 46cfba89..8aabc7e8 100644 --- a/src/panel_aerodynamics.jl +++ b/src/panel_aerodynamics.jl @@ -95,6 +95,18 @@ already sampled at three-quarter chord. ifelse(v_rel > 0, -0.25π * pitch_rate * chord / (2 * smooth_norm(v_rel)), zero(pitch_rate * chord)) +""" + inflow_angle(velocity, y_airf, z_airf) + +Angle of attack of `velocity` in the plane square to `y_airf`: from the chord's +projection on that plane, `y_airf × z_airf`, towards `z_airf`. +""" +@inline inflow_angle(velocity, y_airf, z_airf) = + atan(dot(velocity, z_airf), dot(velocity, cross(y_airf, z_airf))) + +"""Row `i` of an `n×3` matrix as a 3-vector.""" +@inline matrix_row(matrix, i) = SVector(matrix[i, 1], matrix[i, 2], matrix[i, 3]) + """ effective_alpha(alpha, deficiency) @@ -115,9 +127,9 @@ leaves it zero. `deficiency` feeds [`effective_alpha`](@ref). """ @inline function panel_inflow(axes, va_vec_1, va_vec_2, v_ind, dva_vec_1=nothing, dva_vec_2=nothing, deficiency=0) - (; x_airf, y_airf, z_airf, chord) = axes + (; y_airf, z_airf, chord) = axes v_eff = 0.5 .* (va_vec_1 .+ va_vec_2) .+ v_ind - alpha = atan(dot(v_eff, z_airf), dot(v_eff, x_airf)) + alpha = inflow_angle(v_eff, y_airf, z_airf) v_span = cross(v_eff, y_airf) pitch_rate = isnothing(dva_vec_1) ? zero(chord) : section_pitch_rate(0.5 .* (dva_vec_1 .+ dva_vec_2), z_airf, chord) @@ -143,8 +155,8 @@ Lift and drag unit vectors, `(; dir_lift, dir_drag)`. `alpha_dir` is the angle t turns the force, `spanwise` the wing's spanwise direction. """ @inline function panel_force_directions(axes, alpha_dir, spanwise) - (; x_airf, y_airf, z_airf) = axes - dir_inflow = cos(alpha_dir) .* x_airf .+ sin(alpha_dir) .* z_airf + (; y_airf, z_airf) = axes + dir_inflow = cos(alpha_dir) .* cross(y_airf, z_airf) .+ sin(alpha_dir) .* z_airf lift_cross = cross(dir_inflow, y_airf) dir_lift = lift_cross ./ smooth_norm(lift_cross) drag_cross = cross(spanwise, dir_lift) diff --git a/src/solver.jl b/src/solver.jl index d6f07f34..508de43f 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -178,7 +178,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref) # Intermediate results lr::LoopResult{P, T} = LoopResult{P, T}() br::BaseResult{P, T} = BaseResult{P, T}() - cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:10] + cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:8] cache_base::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache()] cache_lin::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:4] @@ -373,35 +373,9 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; end width_dist[i] = panel.width - # Geometric AoA using panel-local axes and prescribed - # freestream — scalar ops to avoid allocations - begin - va1 = solver.sol.va_vec_dist[i,1] - va2 = solver.sol.va_vec_dist[i,2] - va3 = solver.sol.va_vec_dist[i,3] - va = sqrt(va1^2 + va2^2 + va3^2) - x1 = solver.sol._x_airf_dist[i,1] - x2 = solver.sol._x_airf_dist[i,2] - x3 = solver.sol._x_airf_dist[i,3] - x_norm = sqrt(x1^2 + x2^2 + x3^2) - z1 = solver.sol._z_airf_dist[i,1] - z2 = solver.sol._z_airf_dist[i,2] - z3 = solver.sol._z_airf_dist[i,3] - z_norm = sqrt(z1^2 + z2^2 + z3^2) - if va == 0 || x_norm == 0 || z_norm == 0 - alpha_geometric_dist[i] = NaN - else - inv_va = -1.0 / va - vu1 = va1 * inv_va - vu2 = va2 * inv_va - vu3 = va3 * inv_va - v_tangential = (x1*vu1+x2*vu2+x3*vu3) / x_norm - v_normal = (z1*vu1+z2*vu2+z3*vu3) / z_norm - alpha_geometric_dist[i] = atan( - -v_normal, -v_tangential) - end - end - + va_vec = matrix_row(solver.sol.va_vec_dist, i) + alpha_geometric_dist[i] = iszero(va_vec) ? NaN : inflow_angle(va_vec, + matrix_row(solver.sol._y_airf_dist, i), matrix_row(solver.sol._z_airf_dist, i)) end lift = solver.sol.lift_dist @@ -415,8 +389,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; body_aero, gamma_new, solver.core_radius_fraction, - solver.sol._z_airf_dist, - solver.sol._x_airf_dist, solver.sol.va_vec_dist, solver.br.va_dist, solver.br.va_unit_dist @@ -775,9 +747,6 @@ end relative_velocity_crossz, v_acrossz_dist, z_airf_dist, - x_airf_dist, - v_normal_dist, - v_tangential_dist, va_magw_dist, cl_dist, chord_dist, @@ -806,20 +775,11 @@ end end @inbounds for i in 1:n_panels - v_normal_dist[i] = - z_airf_dist[i,1]*relative_velocity_dist[i,1] + - z_airf_dist[i,2]*relative_velocity_dist[i,2] + - z_airf_dist[i,3]*relative_velocity_dist[i,3] - v_tangential_dist[i] = - x_airf_dist[i,1]*relative_velocity_dist[i,1] + - x_airf_dist[i,2]*relative_velocity_dist[i,2] + - x_airf_dist[i,3]*relative_velocity_dist[i,3] - solver.lr.v_span_dist[i] = - y_airf_dist[i,1]*relative_velocity_dist[i,1] + - y_airf_dist[i,2]*relative_velocity_dist[i,2] + - y_airf_dist[i,3]*relative_velocity_dist[i,3] + v_rel = matrix_row(relative_velocity_dist, i) + y_airf = matrix_row(y_airf_dist, i) + solver.lr.alpha_dist[i] = inflow_angle(v_rel, y_airf, matrix_row(z_airf_dist, i)) + solver.lr.v_span_dist[i] = dot(y_airf, v_rel) end - solver.lr.alpha_dist .= atan.(v_normal_dist, v_tangential_dist) @inbounds for i in 1:n_panels solver.lr.v_rel_dist[i] = smooth_sqrt( @@ -945,7 +905,6 @@ function gamma_loop!( ) where {P, U, T} va_vec_dist = solver.sol.va_vec_dist chord_dist = solver.sol._chord_dist - x_airf_dist = solver.sol._x_airf_dist y_airf_dist = solver.sol._y_airf_dist z_airf_dist = solver.sol._z_airf_dist solver.lr.converged = false @@ -961,8 +920,6 @@ function gamma_loop!( relative_velocity_crossz = solver.cache[6][va_vec_dist] v_acrossz_dist = solver.cache[7][va_vec_dist] cl_dist = solver.cache[8][solver.lr.gamma_new] - v_normal_dist = solver.cache[9][solver.lr.gamma_new] - v_tangential_dist = solver.cache[10][solver.lr.gamma_new] AIC_x = @view body_aero.AIC[:, :, 1] AIC_y = @view body_aero.AIC[:, :, 2] @@ -991,9 +948,7 @@ function gamma_loop!( velocity_view_x, velocity_view_y, velocity_view_z, va_vec_dist, induced_velocity_all, relative_velocity_dist, y_airf_dist, relative_velocity_crossz, v_acrossz_dist, - z_airf_dist, x_airf_dist, - v_normal_dist, v_tangential_dist, - va_magw_dist, cl_dist, chord_dist, + z_airf_dist, va_magw_dist, cl_dist, chord_dist, ) @inbounds for i in 1:n_panels residual[i] -= gamma_iter[i] @@ -1013,9 +968,7 @@ function gamma_loop!( velocity_view_x, velocity_view_y, velocity_view_z, va_vec_dist, induced_velocity_all, relative_velocity_dist, y_airf_dist, relative_velocity_crossz, v_acrossz_dist, - z_airf_dist, x_airf_dist, - v_normal_dist, v_tangential_dist, - va_magw_dist, cl_dist, chord_dist, + z_airf_dist, va_magw_dist, cl_dist, chord_dist, ) inv_step = 1.0 / step for i in 1:n_panels @@ -1042,9 +995,7 @@ function gamma_loop!( velocity_view_x, velocity_view_y, velocity_view_z, va_vec_dist, induced_velocity_all, relative_velocity_dist, y_airf_dist, relative_velocity_crossz, v_acrossz_dist, - z_airf_dist, x_airf_dist, - v_normal_dist, v_tangential_dist, - va_magw_dist, cl_dist, chord_dist, + z_airf_dist, va_magw_dist, cl_dist, chord_dist, ) @inbounds for i in 1:n_panels residual_perturbed[i] -= gamma_perturbed[i] @@ -1108,9 +1059,6 @@ function gamma_loop!( relative_velocity_crossz, v_acrossz_dist, z_airf_dist, - x_airf_dist, - v_normal_dist, - v_tangential_dist, va_magw_dist, cl_dist, chord_dist, diff --git a/test/panel/test_panel.jl b/test/panel/test_panel.jl index 63fc8231..05c8b42a 100644 --- a/test/panel/test_panel.jl +++ b/test/panel/test_panel.jl @@ -1,8 +1,11 @@ using VortexStepMethod: Panel, Section, calculate_relative_alpha_and_relative_velocity, calculate_cl, calculate_cd_cm, reinit!, INVISCID, POLAR_VECTORS, MVec3 -using VortexStepMethod: panel_axes, panel_span_vector +using VortexStepMethod: panel_axes, panel_span_vector, panel_inflow, panel_force_directions using Interpolations: linear_interpolation, Line using LinearAlgebra using Test +if !@isdefined(test_data_path) + include("../test_data_utils.jl") +end function create_panel(section1::Section, section2::Section) # Calculate panel geometry @@ -215,4 +218,21 @@ end @test isapprox(flipped.z_airf, -axes.z_airf; atol=1e-12) @test isapprox(flipped.x_airf, axes.x_airf; atol=1e-12) end + + @testset "alpha is measured in the plane square to the bound vortex" begin + for v_eff in ([1.0, 0.0, 0.0], [1.0, 0.3, 0.1]), orient in (1, -1) + oriented = panel_axes(le_1, te_1, le_2, te_2, 0.5, orient) + (; x_airf, y_airf, z_airf) = oriented + alpha = panel_inflow(oriented, v_eff, v_eff, zeros(3)).alpha + @test alpha ≈ crossflow_alpha(v_eff, x_airf, y_airf, z_airf) + along_span = v_eff .+ 2.0 .* y_airf + @test panel_inflow(oriented, along_span, along_span, zeros(3)).alpha ≈ alpha + dir_lift = panel_force_directions(oriented, alpha, y_airf).dir_lift + v_crossflow = v_eff .- dot(v_eff, y_airf) .* y_airf + @test abs(dot(dir_lift, v_crossflow)) < 1e-12 + @test abs(dot(dir_lift, y_airf)) < 1e-12 + end + @test rad2deg(panel_inflow(axes, [1.0, 0, 0], [1.0, 0, 0], zeros(3)).alpha) ≈ + -20.097 atol = 1e-3 + end end diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 089b3353..25ebf514 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -328,3 +328,33 @@ end @test VortexStepMethod.finite_full(ForwardDiff.Dual(1.0, 2.0)) @test !VortexStepMethod.finite_full(ForwardDiff.Dual(1.0, Inf)) end + +@testset "alpha is the crossflow angle on a swept, tapered wing" begin + polar = (deg2rad.([-20.0, 0.0, 20.0]), [-1.8, 0.2, 2.2], fill(0.02, 3), fill(-0.05, 3)) + wing = Wing(8; spanwise_distribution=LINEAR) + add_section!(wing, [1.5, 4.0, 0.3], [2.2, 4.0, 0.1], POLAR_VECTORS, polar) + add_section!(wing, [0.0, 0.0, 0.0], [2.0, 0.0, 0.0], POLAR_VECTORS, polar) + add_section!(wing, [1.5, -4.0, 0.3], [2.2, -4.0, 0.1], POLAR_VECTORS, polar) + refine!(wing) + body_aero = BodyAerodynamics([wing]) + va = [15.0, 1.0, 1.5] + set_va!(body_aero, va) + + for correct_aoa in (false, true) + solver = Solver(length(body_aero.panels), 3; correct_aoa) + sol = solve!(solver, body_aero) + @test sol.solver_status == FEASIBLE + gamma = solver.lr.gamma_new + for (i, panel) in enumerate(body_aero.panels) + frame = (panel.x_airf, panel.y_airf, panel.z_airf) + induced = body_aero.AIC[i, :, :]' * gamma + @test solver.lr.alpha_dist[i] ≈ + crossflow_alpha(va .+ induced, frame...) atol = 1e-6 + @test sol.alpha_geometric_dist[i] ≈ crossflow_alpha(va, frame...) atol = 1e-10 + correct_aoa || continue + induced_aero_center = body_aero.AIC_aero_center[i, :, :]' * gamma + @test sol.alpha_dist[i] ≈ + crossflow_alpha(va .+ induced_aero_center, frame...) atol = 1e-8 + end + end +end diff --git a/test/test_data_utils.jl b/test/test_data_utils.jl index 00cdf7f6..4ef182b8 100644 --- a/test/test_data_utils.jl +++ b/test/test_data_utils.jl @@ -213,3 +213,14 @@ function get_complete_settings_file(module_name) return create_temp_wing_settings(module_name, basename(get_standard_wing_file(module_name))) end end + +""" + crossflow_alpha(velocity, x_airf, y_airf, z_airf) + +Angle of `velocity` in the plane square to `y_airf`, from the chord `x_airf` projected +onto that plane towards `z_airf`: Gaunaa et al. 2026 Eqs. (1) and (3), formed directly. +""" +function crossflow_alpha(velocity, x_airf, y_airf, z_airf) + chord_in_plane = normalize(x_airf .- dot(x_airf, y_airf) .* y_airf) + return atan(dot(velocity, z_airf), dot(velocity, chord_in_plane)) +end