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