From 8cfb9eed944aa8edd603b6c4bda7b08a3d13f022 Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 24 Aug 2026 13:47:34 +0200 Subject: [PATCH] Write the panel aerodynamics once, as traceable functions SymbolicAWEModels re-expressed this package's per-panel force symbolically so it could carry it in an ODE right-hand side: the airfoil axes, the chord blend, the angle of attack, the lift and drag directions, the force and the pitching moment, all transcribed by hand into ModelingToolkit equations. Two spellings of one derivation, and the recent billow-axis fix is exactly the kind of edit that silently changes only one of them. Pull that algebra into src/panel_aerodynamics.jl as pure functions of the section geometry and the flow: panel_axes, panel_inflow, panel_force_directions, panel_loads, and the small helpers around them. They are generic in the number type, free of in-place buffers, and branch-free apart from two ifelse guards, so the same definitions serve a Float64 solver, a ForwardDiff.Dual linearization and a symbolic trace. update_panel_properties!, init_pos!, calc_forces! and calculate_results now call them instead of spelling the algebra out three times, which is one copy fewer inside this package too. Panel geometry is unchanged bit for bit. Forces, moments and coefficients agree to within 1 ulp, the products having been reassociated: on a POLAR_VECTORS wing 47 of 55 reported quantities are bitwise identical and the largest relative difference is 2.1e-16. calc_forces! stays zero-allocation, and the geometry pass got slightly cheaper (update_panel_properties! 64 -> 0 bytes, reinit! 4208 -> 4144) because the chord no longer subtracts MVectors. section_pitch_rate gains a three-argument form taking the trailing minus leading edge apparent wind directly, which is the shape a symbolic consumer has; the four-argument form is unchanged. effective_alpha and panel_inflow's deficiency argument carry an unsteady lag into the angle the polars are read at. Nothing in the solver has unsteady state to feed them, but they belong beside the rest of the panel physics rather than in whoever grows it first. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 26 +++++ Project.toml | 2 +- docs/src/private_functions.md | 22 +++- src/VortexStepMethod.jl | 1 + src/body_aerodynamics.jl | 67 ++----------- src/panel.jl | 38 +++---- src/panel_aerodynamics.jl | 184 ++++++++++++++++++++++++++++++++++ src/solver.jl | 82 +++------------ src/wing_geometry.jl | 89 ++++++---------- 9 files changed, 309 insertions(+), 202 deletions(-) create mode 100644 src/panel_aerodynamics.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 86785ef9..c95a39db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,31 @@ # Changelog +## Unreleased + +### Added +- Shared panel aerodynamics (`src/panel_aerodynamics.jl`): the per-panel physics + written once as pure, branch-free, number-type-generic functions of the section + geometry and the flow — `panel_axes`, `panel_inflow`, `panel_force_directions`, + `panel_loads` and the small helpers around them. `update_panel_properties!`, + `init_pos!`, `calc_forces!` and `calculate_results` now call them instead of + spelling the algebra out three times, and `SymbolicAWEModels` traces the same + functions with symbolic arguments to build its equations, so the two packages + can no longer drift apart. Panel geometry is unchanged bit for bit; forces, + moments and coefficients agree to within 1 ulp, the products having been + reassociated. `calc_forces!` stays zero-allocation. +- `effective_alpha` and the `deficiency` argument of `panel_inflow` carry an + unsteady lag (a Wagner indicial deficiency) into the angle the polars are read + at, while the geometric angle still turns the force. Unused by the solver, + which has no unsteady state; it is shared so a symbolic consumer that does have + one reads the same definition. + +### Changed +- `section_pitch_rate` gained a three-argument form taking the trailing minus + leading edge apparent wind directly. The four-argument form is unchanged. +- Norms inside the shared panel aerodynamics are floored (`smooth_norm`) rather + than guarded by branches, so the expressions are differentiable and traceable. + The floor is `1e-12` m, which moves no reported quantity. + ## VortexStepMethod v4.1.2 2026-08-22 ### Fixed diff --git a/Project.toml b/Project.toml index 1a4da1cf..84001964 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "VortexStepMethod" uuid = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" authors = ["1-Bart-1 ", "Oriol Cayon and contributors"] -version = "4.1.2" +version = "4.2.0" [workspace] projects = ["examples", "docs", "test"] diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index d4683042..b48ed3e8 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -17,7 +17,6 @@ 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 @@ -31,6 +30,27 @@ smooth_distribution! make_dual_shadow ``` +### Panel aerodynamics +The per-panel aerodynamics, written once as pure, branch-free functions of the +section geometry and the flow. `SymbolicAWEModels` traces the same functions with +symbolic arguments to build its equations, so both packages evaluate one +definition of the physics. +```@docs +SMOOTH_FLOOR +smooth_norm +panel_span_vector +panel_chord_weight +panel_axes +effective_alpha +panel_inflow +dynamic_pressure +flow_curvature_cm +panel_force_directions +panel_moment +panel_couple_force +panel_loads +``` + ### Induced velocities ```@docs velocity_3D_bound_vortex! diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 79998e86..087a346a 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -414,6 +414,7 @@ function help(url) end # Include core functionality +include("panel_aerodynamics.jl") include("settings.jl") include("section_aero.jl") include("wing_geometry.jl") diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 9c3c1dc1..a1f3cdb9 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -727,38 +727,6 @@ function compute_panel_center_of_pressures( return panel_cp_locations end -""" - 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) @@ -927,43 +895,26 @@ function calculate_results( cross3!(dir_side_ref, dir_lift_ref, va_ref_unit) q_ref = 0.5 * density * va_ref_mag^2 - induced_va_airfoil = body_aero.work_vectors[4] - dir_lift_induced_va = body_aero.work_vectors[5] - dir_drag_induced_va = body_aero.work_vectors[6] lift_induced_va = body_aero.work_vectors[7] drag_induced_va = body_aero.work_vectors[8] dir_lift_prescribed_va = body_aero.work_vectors[9] temp_vec = body_aero.work_vectors[10] + spanwise_unit = SVector{3}(spanwise_direction) # Main calculation loop for (i, panel) in enumerate(panels) panel_area = panel.chord * panel.width area_all_panels += panel_area - alpha_corrected_i = alpha_corrected[i] - c_alpha = cos(alpha_corrected_i) - s_alpha = sin(alpha_corrected_i) - @inbounds for k in 1:3 - induced_va_airfoil[k] = c_alpha * panel.x_airf[k] + - s_alpha * panel.z_airf[k] - end - normalize3!(induced_va_airfoil) - - cross3!(dir_lift_induced_va, - induced_va_airfoil, panel.y_airf) - normalize3!(dir_lift_induced_va) - cross3!(dir_drag_induced_va, - spanwise_direction, dir_lift_induced_va) - normalize3!(dir_drag_induced_va) - - q_lift = 0.5 * density * v_a_dist[i]^2 - lift_i = cl_array[i] * q_lift * chord_array[i] - drag_i = cd_array[i] * q_lift * chord_array[i] - moment_i = cm_array[i] * q_lift * chord_array[i]^2 - + axes = panel_axes(panel) + dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit) + loads = panel_loads(axes, dirs, + dynamic_pressure(density, density, v_a_dist[i]), + cl_array[i], cd_array[i], cm_array[i]) + moment_i = loads.moment @inbounds for k in 1:3 - lift_induced_va[k] = lift_i * dir_lift_induced_va[k] - drag_induced_va[k] = drag_i * dir_drag_induced_va[k] + lift_induced_va[k] = loads.lift * dirs.dir_lift[k] + drag_induced_va[k] = loads.drag * dirs.dir_drag[k] end va_panel_mag = va_norm_array[i] diff --git a/src/panel.jl b/src/panel.jl index 782ef330..83edc848 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -100,16 +100,15 @@ function init_pos!( panel.LE_point_1 .= section_1.LE_point panel.TE_point_2 .= section_2.TE_point panel.LE_point_2 .= section_2.LE_point - panel.chord = ( - norm(panel.TE_point_1 - panel.LE_point_1) + - norm(panel.TE_point_2 - panel.LE_point_2) - ) / 2 + panel.chord = panel_chord( + SVector{3}(section_1.LE_point), SVector{3}(section_1.TE_point), + SVector{3}(section_2.LE_point), SVector{3}(section_2.TE_point)) panel.corner_points[:, 1] = panel.LE_point_1 panel.corner_points[:, 2] = panel.TE_point_1 panel.corner_points[:, 3] = panel.TE_point_2 panel.corner_points[:, 4] = panel.LE_point_2 vec .= bound_point_2 .- bound_point_1 - panel.width = norm(vec) + panel.width = smooth_norm(vec) reinit!(panel.filaments[1], bound_point_2, bound_point_1, vec) reinit!(panel.filaments[2], bound_point_1, panel.TE_point_1, vec) reinit!(panel.filaments[3], panel.TE_point_2, bound_point_2, vec) @@ -268,6 +267,18 @@ function reinit!( end +""" + panel_axes(panel::Panel) + +The panel's stored airfoil frame and size in the shape [`panel_axes`](@ref) +returns, so a panel built by the geometry pass feeds the same functions as one built +straight from section points. +""" +@inline panel_axes(panel::Panel) = + (x_airf = SVector{3}(panel.x_airf), y_airf = SVector{3}(panel.y_airf), + z_airf = SVector{3}(panel.z_airf), chord = panel.chord, + width = panel.width) + """ calculate_relative_alpha_and_relative_velocity(panel::Panel, induced_velocity::Vector{Float64}) @@ -286,14 +297,8 @@ function calculate_relative_alpha_and_relative_velocity( panel::Panel{T}, induced_velocity::AbstractVector{T} ) where T - # Calculate relative velocity and angle of attack - # Constants throughout iterations: panel.va, panel.x_airf, panel.y_airf - relative_velocity = panel.va .+ induced_velocity - v_normal = dot(panel.z_airf, relative_velocity) - v_tangential = dot(panel.x_airf, relative_velocity) - alpha = atan(v_normal, v_tangential) - - return alpha, relative_velocity + flow = panel_inflow(panel_axes(panel), panel.va, panel.va, induced_velocity) + return flow.alpha, flow.v_eff end """ @@ -302,11 +307,8 @@ end Calculate relative angle of attack and relative velocity of the panel. """ function calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity) - relative_velocity = panel.va + induced_velocity - v_normal = dot(panel.z_airf, relative_velocity) - v_tangential = dot(panel.x_airf, relative_velocity) - alpha = atan(v_normal, v_tangential) - return alpha, relative_velocity + flow = panel_inflow(panel_axes(panel), panel.va, panel.va, induced_velocity) + return flow.alpha, flow.v_eff end """ diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl new file mode 100644 index 00000000..3e19da90 --- /dev/null +++ b/src/panel_aerodynamics.jl @@ -0,0 +1,184 @@ +"""Length [`smooth_norm`](@ref) adds in quadrature to keep a norm positive at zero.""" +const SMOOTH_FLOOR = 1e-12 + +""" + smooth_norm(vec, floor=SMOOTH_FLOOR) + +Norm floored in quadrature: positive and differentiable at the origin, and +branch-free, which is what makes the functions below traceable. +""" +@inline smooth_norm(vec, floor=SMOOTH_FLOOR) = sqrt(sum(abs2, vec) + floor^2) +@inline smooth_norm(x::Number, floor=SMOOTH_FLOOR) = sqrt(x^2 + floor^2) + +""" + panel_span_vector(le_1, te_1, le_2, te_2) + +Quarter-chord vector from section 2 to section 1: the panel's spanwise axis, its +length the span width. +""" +@inline panel_span_vector(le_1, te_1, le_2, te_2) = + (0.75 .* le_1 .+ 0.25 .* te_1) .- (0.75 .* le_2 .+ 0.25 .* te_2) + +@inline panel_chord(le_1, te_1, le_2, te_2) = + 0.5 * (smooth_norm(te_1 .- le_1) + smooth_norm(te_2 .- le_2)) + +""" + panel_chord_weight(width_prev, width_own, width_next) + +Section 1's share of the edge blend setting a panel's chord direction, weighted by +panel spacing so a panel leans towards a narrower neighbour. `nothing` for a +missing neighbour at a tip; a lone panel gets `0.5`. +""" +@inline function panel_chord_weight(width_prev, width_own, width_next) + share = if isnothing(width_prev) && isnothing(width_next) + one(width_own) / 2 + elseif isnothing(width_prev) + width_own / (width_own + width_next) + elseif isnothing(width_next) + width_prev / (width_prev + width_own) + else + 0.25 * (width_prev / (width_prev + width_own) + + width_own / (width_own + width_next) + 1) + end + return 1 - share +end + +""" + panel_axes(le_1, te_1, le_2, te_2, chord_weight=0.5, orient=1) + +Airfoil frame and size of the panel between two sections, as +`(; x_airf, y_airf, z_airf, chord, width)`. + +`chord_weight` ([`panel_chord_weight`](@ref)) enters as an offset from the +midpoint rather than as `w·p₁ + (1-w)·p₂`: the two are equal, but the offset form +leaves a constant term to fold, which a symbolic consumer builds several times +faster. `orient` is `±1`, flipping `y_airf`/`z_airf` so the frame does not depend +on section ordering. +""" +@inline function panel_axes(le_1, te_1, le_2, te_2, chord_weight=0.5, orient=1) + lean = chord_weight - 0.5 + chord_vec = (0.5 .* (te_1 .+ te_2) .- 0.5 .* (le_1 .+ le_2)) .+ + lean .* ((te_1 .- te_2) .- (le_1 .- le_2)) + span_vec = panel_span_vector(le_1, te_1, le_2, te_2) + width = smooth_norm(span_vec) + x_airf = chord_vec ./ smooth_norm(chord_vec) + y_airf = orient .* (span_vec ./ width) + z_cross = cross(x_airf, le_1 .- le_2) + z_airf = orient .* (z_cross ./ smooth_norm(z_cross)) + return (; x_airf, y_airf, z_airf, + chord=panel_chord(le_1, te_1, le_2, te_2), width) +end + +""" + section_pitch_rate(delta_va, z_airf, chord) + section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord) + +Rate a section rotates about its own spanwise axis, positive nose-up. `delta_va` +is the trailing minus leading edge apparent wind; apparent wind is +`wind - velocity`, so that is the leading minus trailing edge velocity, hence the +reversed order in the four-argument form. Chordwise wind variation enters here +too. Builds a `pitch_rate_dist` for [`set_va!`](@ref) on a deforming structure, +where no single body rate describes every section. +""" +@inline section_pitch_rate(delta_va, z_airf, chord) = + ifelse(chord > 0, dot(delta_va, z_airf) / smooth_norm(chord), zero(chord)) +@inline section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord) = + section_pitch_rate(velocity_leading .- velocity_trailing, z_airf, chord) + +""" + flow_curvature_cm(pitch_rate, chord, v_rel) + +Thin-airfoil quarter-chord moment increment of a section pitching about its own +spanwise axis: `Δcm = -(π/4)·q̂` with `q̂ = q·c/(2·v_rel)`, `q` positive nose-up. +Pivot-independent, and lift needs no matching correction because the inflow is +already sampled at three-quarter chord. +""" +@inline flow_curvature_cm(pitch_rate, chord, v_rel) = + ifelse(v_rel > 0, -0.25π * pitch_rate * chord / (2 * smooth_norm(v_rel)), + zero(pitch_rate * chord)) + +""" + effective_alpha(alpha, deficiency) + +Angle the polars are read at: the geometric inflow angle less an unsteady lag. The +geometric angle still turns the force, so a lag shifts the coefficients only. +""" +@inline effective_alpha(alpha, deficiency) = alpha - deficiency + +""" + panel_inflow(axes, va_1, va_2, v_ind, dva_1=nothing, dva_2=nothing, + deficiency=0) + +Flow a panel sees, as `(; v_eff, alpha, alpha_eff, v_span, pitch_rate)`, where +`v_span` is the effective velocity across the span. `axes` is a +[`panel_axes`](@ref) result. `dva_1`/`dva_2` are the sections' trailing minus +leading edge apparent wind, giving the [`section_pitch_rate`](@ref); `nothing` +leaves it zero. `deficiency` feeds [`effective_alpha`](@ref). +""" +@inline function panel_inflow(axes, va_1, va_2, v_ind, dva_1=nothing, + dva_2=nothing, deficiency=0) + (; x_airf, y_airf, z_airf, chord) = axes + v_eff = 0.5 .* (va_1 .+ va_2) .+ v_ind + alpha = atan(dot(v_eff, z_airf), dot(v_eff, x_airf)) + v_span = cross(v_eff, y_airf) + pitch_rate = isnothing(dva_1) ? zero(chord) : + section_pitch_rate(0.5 .* (dva_1 .+ dva_2), z_airf, chord) + return (; v_eff, alpha, alpha_eff=effective_alpha(alpha, deficiency), + v_span, pitch_rate) +end + +""" + dynamic_pressure(rho_1, rho_2, v_span) + +Panel dynamic pressure from its two section densities and the +[`panel_inflow`](@ref) span-wise velocity, as that vector or its magnitude. +""" +@inline dynamic_pressure(rho_1, rho_2, v_span::AbstractVector) = + 0.25 * (rho_1 + rho_2) * dot(v_span, v_span) +@inline dynamic_pressure(rho_1, rho_2, v_rel::Number) = + 0.25 * (rho_1 + rho_2) * v_rel^2 + +""" + panel_force_directions(axes, alpha_dir, spanwise) + +Lift and drag unit vectors, `(; dir_lift, dir_drag)`. `alpha_dir` is the angle that +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 + lift_cross = cross(dir_inflow, y_airf) + dir_lift = lift_cross ./ smooth_norm(lift_cross) + drag_cross = cross(spanwise, dir_lift) + return (; dir_lift, dir_drag=drag_cross ./ smooth_norm(drag_cross)) +end + +"""Panel pitching moment per unit span about `y_airf`, positive nose-up.""" +@inline panel_moment(cm, q_dyn, chord) = cm * q_dyn * chord^2 + +""" + panel_couple_force(cm, q_dyn, chord, width, scale=1) + +[`panel_moment`](@ref) as a force couple: the magnitude two opposed forces one +chord apart need to produce it. A particle model places the moment this way, along +the panel normal at its leading and trailing edge. +""" +@inline panel_couple_force(cm, q_dyn, chord, width, scale=1) = + scale * width * cm * q_dyn * chord + +""" + panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1) + +Panel load from its polar coefficients, [`panel_axes`](@ref) and +[`panel_force_directions`](@ref), as `(; lift, drag, moment, force, +pitching_moment)`. The first three are per unit span; `force` and +`pitching_moment` are the whole panel's, `scale` included. +""" +@inline function panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1) + (; chord, width) = axes + lift = cl * q_dyn * chord + drag = cd * q_dyn * chord + moment = panel_moment(cm, q_dyn, chord) + force = (scale * width) .* (lift .* dirs.dir_lift .+ drag .* dirs.dir_drag) + return (; lift, drag, moment, force, pitching_moment=scale * width * moment) +end diff --git a/src/solver.jl b/src/solver.jl index 02227f0b..48607f3e 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -331,16 +331,10 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; end - # create an alias for the three vertical output vectors lift = solver.sol.lift_dist drag = solver.sol.drag_dist panel_moment_dist = solver.sol.panel_moment_dist - # Compute using fused broadcasting (no intermediate allocations) - @. lift = cl_dist * 0.5 * density * v_a_dist^2 * solver.sol._chord_dist - @. drag = cd_dist * 0.5 * density * v_a_dist^2 * solver.sol._chord_dist - @. panel_moment_dist = cm_dist * 0.5 * density * v_a_dist^2 * solver.sol._chord_dist^2 - # Calculate alpha corrections based on model type if solver.correct_aoa && aerodynamic_model_type == VSM # 64 bytes update_effective_angle_of_attack!( @@ -369,76 +363,32 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; projected_area = body_aero.projected_area c_ref = body_aero.c_ref - wv = body_aero.work_vectors - dir_iva = wv[1] - dir_lift = wv[2] - dir_drag = wv[3] - lift_va = wv[4] - drag_va = wv[5] - r_vec = wv[6] - f_tmp = wv[7] - cross_tmp = wv[8] + spanwise_unit = SVector{3, T}(spanwise_direction) for (i, panel) in enumerate(panels) - - ### Lift and Drag ### panel_area = panel.chord * panel.width area_all_panels += panel_area panel_areas[i] = panel_area - # Calculate induced velocity direction - alpha_corrected_i = alpha_corrected[i] - c_alpha = cos(alpha_corrected_i) - s_alpha = sin(alpha_corrected_i) - @inbounds for k in 1:3 - dir_iva[k] = c_alpha * panel.x_airf[k] + - s_alpha * panel.z_airf[k] - end - normalize3!(dir_iva) - - # Calculate lift and drag directions - cross3!(dir_lift, dir_iva, panel.y_airf) - normalize3!(dir_lift) - cross3!(dir_drag, spanwise_direction, dir_lift) - normalize3!(dir_drag) - - # Calculate force vectors - li = lift[i] - di = drag[i] - @inbounds for k in 1:3 - lift_va[k] = li * dir_lift[k] - drag_va[k] = di * dir_drag[k] - end - - # Body frame forces - width = panel.width - @inbounds for k in 1:3 - solver.sol.f_body_3D[k, i] = (lift_va[k] + - drag_va[k]) * width - end - - # Calculate the moments - m_scale = panel_moment_dist[i] * width + axes = panel_axes(panel) + dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit) + loads = panel_loads(axes, dirs, + dynamic_pressure(density, density, v_a_dist[i]), + cl_dist[i], cd_dist[i], cm_dist[i]) + lift[i] = loads.lift + drag[i] = loads.drag + panel_moment_dist[i] = loads.moment + + force = loads.force + arm_vec = SVector{3, T}(panel.aero_center) - SVector{3, T}(reference_point) + moment = loads.pitching_moment .* axes.y_airf .+ cross(arm_vec, force) @inbounds for k in 1:3 - r_vec[k] = panel.aero_center[k] - - reference_point[k] - f_tmp[k] = solver.sol.f_body_3D[k, i] - end - cross3!(cross_tmp, r_vec, f_tmp) - @inbounds for k in 1:3 - solver.sol.m_body_3D[k, i] = m_scale * - panel.y_airf[k] + cross_tmp[k] + solver.sol.f_body_3D[k, i] = force[k] + solver.sol.m_body_3D[k, i] = moment[k] end - # Moment distribution (moment on each panel) arm = (moment_frac - 0.25) * panel.chord - ftotal_dot_z = 0.0 - @inbounds for k in 1:3 - ftotal_dot_z += (lift_va[k] + drag_va[k]) * - panel.z_airf[k] - end - moment_dist[i] = (ftotal_dot_z * arm + - panel_moment_dist[i]) * width + moment_dist[i] = dot(force, axes.z_airf) * arm + loads.pitching_moment end # Python parity: normalize with area-weighted reference velocity for distributed inflow. diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index f64b01b8..1025ed41 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -147,6 +147,7 @@ Structure to hold calculated panel properties. - `x_airf`::Matrix{Float64}: Vector of unit vectors tangential to chord line - `y_airf`::Matrix{Float64}: Vector of unit vectors in spanwise direction - `z_airf`::Matrix{Float64}: Vector of unit vectors pointing up (cross of x_airf and y_airf) +- `widths`::Vector{Float64}: Span width of each panel """ @with_kw mutable struct PanelProperties{P, T} aero_centers::Matrix{T} = zeros(T, P, 3) @@ -156,6 +157,7 @@ Structure to hold calculated panel properties. x_airf::Matrix{T} = zeros(T, P, 3) y_airf::Matrix{T} = zeros(T, P, 3) z_airf::Matrix{T} = zeros(T, P, 3) + widths::Vector{T} = zeros(T, P) coords::Matrix{T} = zeros(T, 2(P+1), 3) end @@ -180,76 +182,47 @@ function update_panel_properties!(panel_props::PanelProperties{P,T}, section_lis x_airf = panel_props.x_airf y_airf = panel_props.y_airf z_airf = panel_props.z_airf - vec = zeros(MVector{3, T}) - vec2 = zeros(MVector{3, T}) + widths = panel_props.widths @debug "Shape of coordinates: $(size(coords))" - + for i in 1:n_panels coords[2i-1, :] .= section_list[i].LE_point coords[2i, :] .= section_list[i].TE_point coords[2i+1, :] .= section_list[i+1].LE_point coords[2i+2, :] .= section_list[i+1].TE_point end - + @debug "Coordinates: $coords" - + + point(row) = SVector{3, T}(coords[row, 1], coords[row, 2], coords[row, 3]) + section_points(i) = (point(2i-1), point(2i), point(2i+1), point(2i+2)) + for i in 1:n_panels - # Define panel points - @views begin - LE_1 = coords[2i-1, :] # LE_1 - LE_2 = coords[2i+1, :] # LE_2 - TE_2 = coords[2i+2, :] # TE_2 - TE_1 = coords[2i, :] # TE_1 - end - - # Calculate control point position - @views @. vec = coords[2i-1, :] * 0.75 + coords[2i, :] * 0.25 - - (coords[2i+1, :] * 0.75 + coords[2i+2, :] * 0.25) - di = norm(vec) - - ncp = if i == 1 - @views @. vec = coords[2i+1, :] * 0.75 + coords[2i+2, :] * 0.25 - - (coords[2i+3, :] * 0.75 + coords[2i+4, :] * 0.25) - diplus = norm(vec) - di / (di + diplus) + widths[i] = smooth_norm(panel_span_vector(section_points(i)...)) + end + + for i in 1:n_panels + le_1, te_1, le_2, te_2 = section_points(i) + weight = if n_panels == 1 + panel_chord_weight(nothing, widths[1], nothing) + elseif i == 1 + panel_chord_weight(nothing, widths[1], widths[2]) elseif i == n_panels - @views @. vec = coords[2i-3, :] * 0.75 + coords[2i-2, :] * 0.25 - - (coords[2i-1, :] * 0.75 + coords[2i, :] * 0.25) - dimin = norm(vec) - dimin / (dimin + di) + panel_chord_weight(widths[i-1], widths[i], nothing) else - @views @. vec = coords[2i-3, :] * 0.75 + coords[2i-2, :] * 0.25 - - (coords[2i-1, :] * 0.75 + coords[2i, :] * 0.25) - dimin = norm(vec) - @views @. vec = coords[2i+1, :] * 0.75 + coords[2i+2, :] * 0.25 - - (coords[2i+3, :] * 0.75 + coords[2i+4, :] * 0.25) - diplus = norm(vec) - 0.25 * (dimin / (dimin + di) + di / (di + diplus) + 1) - end - ncp = 1 - ncp - - # Calculate points - @. begin - aero_centers[i, :] = (LE_2 * (1 - ncp) + LE_1 * ncp) * 0.75 + - (TE_2 * (1 - ncp) + TE_1 * ncp) * 0.25 - control_points[i, :] = (LE_2 * (1 - ncp) + LE_1 * ncp) * 0.25 + - (TE_2 * (1 - ncp) + TE_1 * ncp) * 0.75 - - bound_points_1[i, :] = LE_1 * 0.75 + TE_1 * 0.25 - bound_points_2[i, :] = LE_2 * 0.75 + TE_2 * 0.25 - end - - # Calculate reference frame vectors - @views begin - @. vec = (control_points[i, :] - aero_centers[i, :]) - @. vec2 = (LE_1 - LE_2) - vec .= vec × vec2 - z_airf[i, :] .= normalize(vec) - @. vec = control_points[i, :] .- aero_centers[i, :] - x_airf[i, :] .= normalize(vec) - @. vec = bound_points_1[i, :] - bound_points_2[i, :] - y_airf[i, :] .= normalize(vec) + panel_chord_weight(widths[i-1], widths[i], widths[i+1]) end + axes = panel_axes(le_1, te_1, le_2, te_2, weight) + + le_blend = weight .* le_1 .+ (1 - weight) .* le_2 + te_blend = weight .* te_1 .+ (1 - weight) .* te_2 + aero_centers[i, :] .= 0.75 .* le_blend .+ 0.25 .* te_blend + control_points[i, :] .= 0.25 .* le_blend .+ 0.75 .* te_blend + bound_points_1[i, :] .= 0.75 .* le_1 .+ 0.25 .* te_1 + bound_points_2[i, :] .= 0.75 .* le_2 .+ 0.25 .* te_2 + x_airf[i, :] .= axes.x_airf + y_airf[i, :] .= axes.y_airf + z_airf[i, :] .= axes.z_airf end return nothing end