diff --git a/CHANGELOG.md b/CHANGELOG.md index 3d1ec3f7..b7c2129c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -68,6 +68,11 @@ ### Fixed +- On a body whose wings span different directions, such as a wing and a vertical fin, + `solve!`, `solve` and `linearize` take each panel's lift, drag and side directions + from its own wing's `spanwise_direction`, not the first wing's. `solve` computes + `wing_span` and `aspect_ratio_projected` from the extent of all wings along the + first wing's span, through the new `calculate_span(wings, spanwise_direction)`. - `ELLIPTIC` initial circulation works on a body with more than one wing, where it threw an `ArgumentError`: each wing gets an ellipse over its own span, along its own `spanwise_direction` and centred on its own mid-span, also for a single wing off y = 0. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index b6e76f60..6174bb8e 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -60,6 +60,7 @@ dynamic_pressure flow_curvature_cm spanwise_flow_drag panel_force_directions +prescribed_va_directions panel_moment panel_couple_force panel_loads @@ -91,6 +92,7 @@ calculate_filaments_for_plotting ```@docs unrefined_deform! unrefined_section_range +panel_range deform! compute_refined_panel_mapping! compute_refined_section_interpolation! diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 0d628876..2ce94f20 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -259,11 +259,20 @@ Indices of the unrefined sections of wing `wing_idx` in a distribution that runs unrefined sections of all wings in order, such as `moment_unrefined_dist`. """ function unrefined_section_range(body_aero::BodyAerodynamics, wing_idx) - offset = 0 - for i in 1:wing_idx-1 - offset += body_aero.wings[i].n_unrefined_sections - end - return offset .+ (1:body_aero.wings[wing_idx].n_unrefined_sections) + wings = body_aero.wings + offset = sum(wing.n_unrefined_sections for wing in view(wings, 1:wing_idx-1); init=0) + return offset .+ (1:wings[wing_idx].n_unrefined_sections) +end + +""" + panel_range(body_aero::BodyAerodynamics, wing_idx) + +Indices of the panels of wing `wing_idx` in `body_aero.panels`. +""" +function panel_range(body_aero::BodyAerodynamics, wing_idx) + wings = body_aero.wings + offset = sum(wing.n_panels for wing in view(wings, 1:wing_idx-1); init=0) + return offset .+ (1:wings[wing_idx].n_panels) end """ @@ -500,16 +509,14 @@ points measured along that wing's `spanwise_direction` from the wing's mid-span. """ function calculate_circulation_distribution_elliptical_wing(gamma_i, body_aero::BodyAerodynamics, gamma_0=1.0) - panel_offset = 0 - for wing in body_aero.wings + for (wing_idx, wing) in enumerate(body_aero.wings) lo, hi = spanwise_extent(wing) axis = normalize(wing.spanwise_direction) - for i in panel_offset .+ (1:wing.n_panels) + for i in panel_range(body_aero, wing_idx) span_position = dot(body_aero.panels[i].control_point, axis) - (lo + hi) / 2 # Clamped: a control point can lie outside the span of the unrefined sections gamma_i[i] = gamma_0 * sqrt(max(0.0, 1 - (2span_position / (hi - lo))^2)) end - panel_offset += wing.n_panels end return nothing end @@ -772,6 +779,17 @@ function set_pitch_rate_dist!(body_aero::BodyAerodynamics, omega) return nothing end +""" + prescribed_va_directions(va, spanwise) + +Lift and side unit vectors, `(; dir_lift, dir_side)`, of inflow `va` on a wing along +`spanwise`: lift normal to both, side normal to lift and `va`. +""" +@inline function prescribed_va_directions(va, spanwise) + dir_lift = normalize(cross(va, spanwise)) + return (; dir_lift, dir_side=cross(dir_lift, va) / norm(va)) +end + """ calculate_results(body_aero::BodyAerodynamics, gamma_new, reference_point, density, core_radius_fraction, mu, alpha_dist, v_rel_dist, chord_dist, @@ -880,8 +898,7 @@ function calculate_results( drag_wing_3D_sum = 0.0 side_wing_3D_sum = 0.0 - # Get wing properties and reference velocity - spanwise_direction = body_aero.wings[1].spanwise_direction + reference_spanwise = SVector{3}(body_aero.wings[1].spanwise_direction) va_ref_vec = MVec3(0.0, 0.0, 0.0) weighted_speed_sq = 0.0 total_area = 0.0 @@ -914,77 +931,62 @@ function calculate_results( @inbounds for k in 1:3 va_ref_unit[k] = va_ref_vec[k] * inv_va_ref end - dir_lift_ref = body_aero.work_vectors[2] - cross3!(dir_lift_ref, va_ref_vec, spanwise_direction) - dir_lift_ref_norm = norm3(dir_lift_ref) - dir_lift_ref_norm > 0.0 || throw(ArgumentError( + reference_dirs = prescribed_va_directions(SVector{3}(va_ref_vec), reference_spanwise) + all(isfinite, reference_dirs.dir_lift) || throw(ArgumentError( "Reference lift direction is undefined because " * "reference flow is parallel to spanwise direction.")) - @inbounds for k in 1:3 - dir_lift_ref[k] /= dir_lift_ref_norm - end - dir_side_ref = body_aero.work_vectors[3] - cross3!(dir_side_ref, dir_lift_ref, va_ref_unit) q_ref = 0.5 * density * va_ref^2 - 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 - - axes = panel_axes(panel) - dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit) - c_span = 0.0 - if is_with_viscous_drag_correction - viscous = spanwise_flow_drag(v_rel_dist[i], v_span_dist[i], panel.chord, - density, mu) - cd_dist[i] += viscous.delta_cd - c_span = viscous.c_span - end - loads = panel_loads(axes, dirs, - dynamic_pressure(density, density, v_rel_dist[i]), - cl_dist[i], cd_dist[i], cm_dist[i]; c_span) - force = loads.force - - va_panel = va_dist[i] - va_panel > 0.0 || throw(ArgumentError( - "Panel $i has non-positive apparent " * - "velocity magnitude.")) - q_panel = 0.5 * density * va_panel^2 - cross3!(dir_lift_prescribed_va, - panel.va_vec, spanwise_direction) - normalize3!(dir_lift_prescribed_va) - - cross3!(temp_vec, dir_lift_prescribed_va, panel.va_vec) - inv_va_panel = 1.0 / va_panel - @inbounds for k in 1:3 - temp_vec[k] *= inv_va_panel - end - lift_prescribed_va = dot(force, dir_lift_prescribed_va) - drag_prescribed_va = dot(force, panel.va_vec) * inv_va_panel - side_prescribed_va = dot(force, temp_vec) - - lift_wing_3D_sum += lift_prescribed_va * - dot3(dir_lift_prescribed_va, dir_lift_ref) - drag_wing_3D_sum += drag_prescribed_va * - (dot3(panel.va_vec, va_ref_unit) * inv_va_panel) - side_wing_3D_sum += side_prescribed_va * - dot3(temp_vec, dir_side_ref) - - inv_q_area = 1.0 / (q_panel * panel_area) - cl_prescribed_va[i] = lift_prescribed_va * inv_q_area - cd_prescribed_va[i] = drag_prescribed_va * inv_q_area - cs_prescribed_va[i] = side_prescribed_va * inv_q_area - - arm = SVector{3}(panel.aero_center) - SVector{3}(reference_point) - moment = loads.pitching_moment .* axes.y_airf .+ cross(arm, force) - @inbounds for k in 1:3 - f_body_3D[k, i] = force[k] - m_body_3D[k, i] = moment[k] + for (wing_idx, wing) in enumerate(body_aero.wings) + spanwise_unit = SVector{3}(wing.spanwise_direction) + for i in panel_range(body_aero, wing_idx) + panel = panels[i] + panel_area = panel.chord * panel.width + area_all_panels += panel_area + + axes = panel_axes(panel) + dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit) + c_span = 0.0 + if is_with_viscous_drag_correction + viscous = spanwise_flow_drag(v_rel_dist[i], v_span_dist[i], panel.chord, + density, mu) + cd_dist[i] += viscous.delta_cd + c_span = viscous.c_span + end + loads = panel_loads(axes, dirs, + dynamic_pressure(density, density, v_rel_dist[i]), + cl_dist[i], cd_dist[i], cm_dist[i]; c_span) + force = loads.force + + va_panel = va_dist[i] + va_panel > 0.0 || throw(ArgumentError( + "Panel $i has non-positive apparent " * + "velocity magnitude.")) + q_panel = 0.5 * density * va_panel^2 + panel_va = SVector{3}(panel.va_vec) + inv_va_panel = 1.0 / va_panel + drag_prescribed_va = dot(force, panel_va) * inv_va_panel + wing_dirs = prescribed_va_directions(panel_va, spanwise_unit) + body_dirs = prescribed_va_directions(panel_va, reference_spanwise) + + lift_wing_3D_sum += dot(force, body_dirs.dir_lift) * + dot(body_dirs.dir_lift, reference_dirs.dir_lift) + drag_wing_3D_sum += drag_prescribed_va * + (dot(panel_va, va_ref_unit) * inv_va_panel) + side_wing_3D_sum += dot(force, body_dirs.dir_side) * + dot(body_dirs.dir_side, reference_dirs.dir_side) + + inv_q_area = 1.0 / (q_panel * panel_area) + cl_prescribed_va[i] = dot(force, wing_dirs.dir_lift) * inv_q_area + cd_prescribed_va[i] = drag_prescribed_va * inv_q_area + cs_prescribed_va[i] = dot(force, wing_dirs.dir_side) * inv_q_area + + arm = SVector{3}(panel.aero_center) - SVector{3}(reference_point) + moment = loads.pitching_moment .* axes.y_airf .+ cross(arm, force) + @inbounds for k in 1:3 + f_body_3D[k, i] = force[k] + m_body_3D[k, i] = moment[k] + end end end @@ -997,7 +999,7 @@ function calculate_results( # Calculate wing geometry properties projected_area = body_aero.projected_area - wing_span = body_aero.wings[1].span + wing_span = calculate_span(body_aero.wings, reference_spanwise) aspect_ratio_projected = wing_span^2 / projected_area # Calculate Reynolds number diff --git a/src/solver.jl b/src/solver.jl index 69dfb1b1..7cc8810c 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -429,46 +429,44 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; area_all_panels = zero(T) panel_areas = solver.sol.panel_area_dist - # Get wing properties - spanwise_direction = body_aero.wings[1].spanwise_direction - - # Calculate wing geometry properties projected_area = body_aero.projected_area c_ref = body_aero.c_ref - - spanwise_unit = SVector{3, T}(spanwise_direction) - for (i, panel) in enumerate(panels) - panel_area = panel.chord * panel.width - area_all_panels += panel_area - panel_areas[i] = panel_area - - axes = panel_axes(panel) - dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit) - c_span = zero(T) - if solver.is_with_viscous_drag_correction - viscous = spanwise_flow_drag(v_rel_dist[i], solver.lr.v_span_dist[i], - panel.chord, density, solver.mu) - cd_dist[i] += viscous.delta_cd - c_span = viscous.c_span - end - loads = panel_loads(axes, dirs, - dynamic_pressure(density, density, v_rel_dist[i]), - cl_dist[i], cd_dist[i], cm_dist[i]; c_span) - 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 - solver.sol.f_body_3D[k, i] = force[k] - solver.sol.m_body_3D[k, i] = moment[k] - end + for (wing_idx, wing) in enumerate(body_aero.wings) + spanwise_unit = SVector{3, T}(wing.spanwise_direction) + for i in panel_range(body_aero, wing_idx) + panel = panels[i] + panel_area = panel.chord * panel.width + area_all_panels += panel_area + panel_areas[i] = panel_area + + axes = panel_axes(panel) + dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit) + c_span = zero(T) + if solver.is_with_viscous_drag_correction + viscous = spanwise_flow_drag(v_rel_dist[i], solver.lr.v_span_dist[i], + panel.chord, density, solver.mu) + cd_dist[i] += viscous.delta_cd + c_span = viscous.c_span + end + loads = panel_loads(axes, dirs, + dynamic_pressure(density, density, v_rel_dist[i]), + cl_dist[i], cd_dist[i], cm_dist[i]; c_span) + 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 + solver.sol.f_body_3D[k, i] = force[k] + solver.sol.m_body_3D[k, i] = moment[k] + end - arm = (moment_frac - 0.25) * panel.chord - moment_dist[i] = dot(force, axes.z_airf) * arm + loads.pitching_moment + arm = (moment_frac - 0.25) * panel.chord + moment_dist[i] = dot(force, axes.z_airf) * arm + loads.pitching_moment + end end # Python parity: normalize with area-weighted reference velocity for distributed inflow. @@ -515,11 +513,11 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; width_unrefined_dist .= 0.0 fill!(unrefined_count_dist, 0) - panel_idx = 1 for (wing_idx, wing) in enumerate(body_aero.wings) if wing.n_unrefined_sections > 0 section_range = unrefined_section_range(body_aero, wing_idx) - for local_panel_idx in 1:wing.n_panels + wing_panels = panel_range(body_aero, wing_idx) + for (local_panel_idx, panel_idx) in enumerate(wing_panels) panel = body_aero.panels[panel_idx] original_section_idx = wing.refined_panel_mapping[local_panel_idx] target_unrefined_idx = section_range[original_section_idx] @@ -541,7 +539,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; width_unrefined_dist[target_unrefined_idx] += panel.width unrefined_count_dist[target_unrefined_idx] += 1 - panel_idx += 1 end # Average coefficients and geometry. width and @@ -563,9 +560,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; # sum of panel widths in the unrefined section end end - else - # Skip panels for wings with no unrefined sections - panel_idx += wing.n_panels end end end diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index c6b3e05d..c849c001 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -1686,25 +1686,31 @@ end """ spanwise_extent(wing::AbstractWing) + spanwise_extent(wings, spanwise_direction) Lowest and highest projection of the unrefined sections' LE and TE points on `wing`'s -`spanwise_direction`, as `(lo, hi)` [m]. +`spanwise_direction`, or of all `wings` on `spanwise_direction`, as `(lo, hi)` [m]. """ -function spanwise_extent(wing::AbstractWing) - axis = normalize(wing.spanwise_direction) - return extrema(dot(point, axis) for section in wing.unrefined_sections +function spanwise_extent(wings, spanwise_direction) + axis = normalize(spanwise_direction) + return extrema(dot(point, axis) for wing in wings + for section in wing.unrefined_sections for point in (section.LE_point, section.TE_point)) end +spanwise_extent(wing::AbstractWing) = spanwise_extent((wing,), wing.spanwise_direction) """ calculate_span(wing::AbstractWing) + calculate_span(wings, spanwise_direction) -Wing span along `spanwise_direction` [m]. +Extent [m] of the unrefined sections of `wing` along its spanwise direction, or of all +`wings` together along `spanwise_direction`. """ -function calculate_span(wing::AbstractWing) - lo, hi = spanwise_extent(wing) +function calculate_span(wings, spanwise_direction) + lo, hi = spanwise_extent(wings, spanwise_direction) return hi - lo end +calculate_span(wing::AbstractWing) = calculate_span((wing,), wing.spanwise_direction) # Project point onto plane @inline function project_onto_plane!(point_proj, point, normal) diff --git a/test/runtests.jl b/test/runtests.jl index 3077e22b..76981ba9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,6 +57,7 @@ function include_selected_tests() should_run_test("solver/test_moment_units.jl") && include("solver/test_moment_units.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_wing_directions.jl") && include("solver/test_wing_directions.jl") should_run_test("solver/test_unrefined_dist.jl") && include("solver/test_unrefined_dist.jl") should_run_test("solver/test_stability.jl") && include("solver/test_stability.jl") should_run_test("verification/test_verification.jl") && include("verification/test_verification.jl") diff --git a/test/solver/test_wing_directions.jl b/test/solver/test_wing_directions.jl new file mode 100644 index 00000000..e1ba4e61 --- /dev/null +++ b/test/solver/test_wing_directions.jl @@ -0,0 +1,74 @@ +using VortexStepMethod +using LinearAlgebra +using Test + +""" + rectangular_wing(rotation, offset) + +A rectangular wing spanning `y`, with every point turned by `rotation` and moved by +`offset`, and its spanwise direction turned with it. +""" +function rectangular_wing(rotation, offset) + chord, span = 1.5, 6.0 + alpha_range = deg2rad.([-10.0, 0.0, 10.0]) + polar = (alpha_range, [-0.6, 0.4, 1.4], fill(0.02, 3), fill(-0.08, 3)) + wing = Wing(10; spanwise_direction=rotation * [0.0, 1.0, 0.0]) + for y in (span / 2, -span / 2) + add_section!(wing, rotation * [0.0, y, 0.0] + offset, + rotation * [chord, y, 0.0] + offset, POLAR_VECTORS, polar) + end + refine!(wing) + return wing +end + +@testset "each wing's panels take their own wing's spanwise direction" begin + # 90° about x leaves the inflow along x unchanged; the second wing is far enough + # away that the two wings do not induce on each other + rotation = [1.0 0.0 0.0; 0.0 0.0 -1.0; 0.0 1.0 0.0] + va = [20.0, 0.0, 0.0] + solo_aero = BodyAerodynamics([rectangular_wing(I(3), zeros(3))]) + pair_aero = BodyAerodynamics([rectangular_wing(I(3), zeros(3)), + rectangular_wing(rotation, [0.0, 1000.0, 0.0])]) + set_va!(solo_aero, va) + set_va!(pair_aero, va) + n_panels = length(solo_aero.panels) + rotated = n_panels .+ (1:n_panels) + + @testset "solve!" begin + solo, pair = map((solo_aero, pair_aero)) do body_aero + solver = Solver(length(body_aero.panels), 2length(body_aero.wings)) + solve!(solver, body_aero) + end + @test pair.solver_status == FEASIBLE + @test pair.f_body_3D[:, 1:n_panels] ≈ solo.f_body_3D rtol = 1e-5 + @test pair.f_body_3D[:, rotated] ≈ rotation * solo.f_body_3D rtol = 1e-5 + @test pair.lift_dist[rotated] ≈ solo.lift_dist rtol = 1e-5 + @test pair.drag_dist[rotated] ≈ solo.drag_dist rtol = 1e-5 + end + + @testset "solve" begin + solo, pair = map((solo_aero, pair_aero)) do body_aero + solver = Solver(length(body_aero.panels), 2length(body_aero.wings)) + solve(solver, body_aero) + end + rotated_forces = pair["F_distribution"][:, rotated] + @test rotated_forces ≈ rotation * solo["F_distribution"] rtol = 1e-5 + for key in ("cl_distribution", "cd_distribution", "cs_distribution") + @test pair[key][rotated] ≈ solo[key] rtol = 1e-5 atol = 1e-8 + end + # along x the inflow makes z the body's lift and y its side direction + @test pair["lift"] ≈ pair["Fz"] rtol = 1e-10 + @test pair["side"] ≈ pair["Fy"] rtol = 1e-10 + @test pair["drag"] ≈ pair["Fx"] rtol = 1e-10 + end +end + +@testset "the projected aspect ratio spans every wing" begin + body_aero = BodyAerodynamics([rectangular_wing(I(3), zeros(3)), + rectangular_wing(I(3), [0.0, 6.0, 0.0])]) + set_va!(body_aero, [20.0, 0.0, 0.0]) + solver = Solver(length(body_aero.panels), 2length(body_aero.wings)) + results = solve(solver, body_aero) + @test results["wing_span"] ≈ 12.0 + @test results["aspect_ratio_projected"] ≈ 12.0^2 / 18.0 +end