From 4965a5fad194ebb0a84e481bdf9ba07e8a206e11 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Thu, 17 Sep 2026 00:42:00 +0200 Subject: [PATCH] Add the Gaunaa et al. 2024 spanwise-flow viscous drag correction as a solver option is_with_viscous_drag_correction (default false) adds spanwise_flow_drag's drag increment and spanwise force to every panel in calc_forces! and calculate_results, from the velocity along y_airf the gamma loop now keeps in LoopResult.v_span_dist. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 7 ++ docs/src/private_functions.md | 1 + src/body_aerodynamics.jl | 78 ++++++--------- src/panel_aerodynamics.jl | 26 ++++- src/settings.jl | 3 + src/solver.jl | 23 ++++- test/runtests.jl | 1 + test/solver/test_viscous_drag_correction.jl | 100 ++++++++++++++++++++ 8 files changed, 185 insertions(+), 54 deletions(-) create mode 100644 test/solver/test_viscous_drag_correction.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 4c00c754..3b837e2a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,13 @@ ## Unreleased +### Added + +- Spanwise-flow viscous drag correction (Gaunaa et al. 2024, + doi:10.1088/1742-6596/2767/2/022068): each section gets a drag increment and a force + along its span from the flow across it, in `solve!`, `solve` and `linearize`. Opt-in + via `is_with_viscous_drag_correction` (default `false`) on the solver settings. + ### Changed - Requires Julia 1.12 or 1.13; 1.10 and 1.11 keep resolving v5.1.1. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index d7807c42..15f636ef 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -52,6 +52,7 @@ effective_alpha panel_inflow dynamic_pressure flow_curvature_cm +spanwise_flow_drag panel_force_directions panel_moment panel_couple_force diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index d04a966d..f0dfad1b 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -748,12 +748,15 @@ end core_radius_fraction, mu, alpha_dist, v_a_dist, chord_dist, x_airf_dist, z_airf_dist, va_dist, va_norm_dist, va_unit_dist, panels::Vector{<:Panel}, is_only_f_and_gamma_output::Bool; - correct_aoa=false, flow_curvature=false) + correct_aoa=false, flow_curvature=false, + is_with_viscous_drag_correction=false, v_span_dist=nothing) 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`. +from `body_aero.omega`. `is_with_viscous_drag_correction` adds +[`spanwise_flow_drag`](@ref) to every section, from the velocity along `y_airf` in +`v_span_dist`. Returns: Dict: Results including forces, coefficients and distributions @@ -777,6 +780,8 @@ function calculate_results( is_only_f_and_gamma_output::Bool; correct_aoa::Bool=false, flow_curvature::Bool=false, + is_with_viscous_drag_correction::Bool=false, + v_span_dist=nothing, ) n_panels = length(panels) @@ -893,8 +898,6 @@ function calculate_results( cross3!(dir_side_ref, dir_lift_ref, va_ref_unit) q_ref = 0.5 * density * va_ref_mag^2 - 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) @@ -906,14 +909,17 @@ function calculate_results( 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_a_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_a_dist[i]), - cl_dist[i], cd_dist[i], cm_dist[i]) - moment_i = loads.moment - @inbounds for k in 1:3 - lift_induced_va[k] = loads.lift * dirs.dir_lift[k] - drag_induced_va[k] = loads.drag * dirs.dir_drag[k] - end + cl_dist[i], cd_dist[i], cm_dist[i]; c_span) + force = loads.force va_panel_mag = va_norm_dist[i] va_panel_mag > 0.0 || throw(ArgumentError( @@ -924,54 +930,32 @@ function calculate_results( panel.va, spanwise_direction) normalize3!(dir_lift_prescribed_va) - lift_prescribed_va = - dot3(lift_induced_va, dir_lift_prescribed_va) + - dot3(drag_induced_va, dir_lift_prescribed_va) - drag_prescribed_va = - (dot3(lift_induced_va, panel.va) + - dot3(drag_induced_va, panel.va)) / va_panel_mag cross3!(temp_vec, dir_lift_prescribed_va, panel.va) inv_vpm = 1.0 / va_panel_mag @inbounds for k in 1:3 temp_vec[k] *= inv_vpm end - side_prescribed_va = - dot3(lift_induced_va, temp_vec) + - dot3(drag_induced_va, temp_vec) - - width = panel.width - @inbounds for k in 1:3 - f_body_3D[k, i] = (lift_induced_va[k] + - drag_induced_va[k]) * width - end + lift_prescribed_va = dot(force, dir_lift_prescribed_va) + drag_prescribed_va = dot(force, panel.va) * inv_vpm + side_prescribed_va = dot(force, temp_vec) - lift_wing_3D_sum += lift_prescribed_va * width * + lift_wing_3D_sum += lift_prescribed_va * dot3(dir_lift_prescribed_va, dir_lift_ref) - drag_wing_3D_sum += drag_prescribed_va * width * - (dot3(panel.va, va_ref_unit) / va_panel_mag) - side_wing_3D_sum += side_prescribed_va * width * + drag_wing_3D_sum += drag_prescribed_va * + (dot3(panel.va, va_ref_unit) * inv_vpm) + side_wing_3D_sum += side_prescribed_va * dot3(temp_vec, dir_side_ref) - inv_qc = 1.0 / (q_panel * panel.chord) - cl_prescribed_va[i] = lift_prescribed_va * inv_qc - cd_prescribed_va[i] = drag_prescribed_va * inv_qc - cs_prescribed_va[i] = side_prescribed_va * inv_qc + 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 - ### Moment ### - # r_vector = panel.aero_center - reference_point - # M_shift = cross(r_vector, f_body_3D[:,i]) - # m_body_3D[:,i] = moment_i * panel.y_airf * width + M_shift - @inbounds for k in 1:3 - dir_lift_prescribed_va[k] = panel.aero_center[k] - - reference_point[k] - drag_induced_va[k] = f_body_3D[k, i] - end - cross3!(temp_vec, - dir_lift_prescribed_va, drag_induced_va) - local_moment_scale = moment_i * width + 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 - m_body_3D[k, i] = local_moment_scale * - panel.y_airf[k] + temp_vec[k] + f_body_3D[k, i] = force[k] + m_body_3D[k, i] = moment[k] end end diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl index 3e19da90..c1e103f6 100644 --- a/src/panel_aerodynamics.jl +++ b/src/panel_aerodynamics.jl @@ -167,18 +167,34 @@ the panel normal at its leading and trailing edge. scale * width * cm * q_dyn * chord """ - panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1) + spanwise_flow_drag(v_a, v_span, chord, density, mu) + +Viscous force increments from spanwise flow (Gaunaa et al. 2024, +doi:10.1088/1742-6596/2767/2/022068), as `(; delta_cd, c_span)`: an addition to the +section drag coefficient and a force coefficient along `y_airf`. Both refer to `v_a`, +the speed normal to the span; `v_span` is the velocity along `y_airf`. +""" +@inline function spanwise_flow_drag(v_a, v_span, chord, density, mu) + f0 = 0.062 * (density * v_a * chord / mu)^(-1 / 7) + skew_factor = (hypot(v_a, v_span) / v_a)^(5 / 7) + return (; delta_cd = f0 * (skew_factor - 1), c_span = f0 * v_span / v_a * skew_factor) +end + +""" + panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1; c_span=0) 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. +`pitching_moment` are the whole panel's, `scale` included. `c_span` adds a force +along `y_airf`, as [`spanwise_flow_drag`](@ref) gives it. """ -@inline function panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1) - (; chord, width) = axes +@inline function panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1; c_span=0) + (; y_airf, 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) + force = (scale * width) .* (lift .* dirs.dir_lift .+ drag .* dirs.dir_drag .+ + (c_span * q_dyn * chord) .* y_airf) return (; lift, drag, moment, force, pitching_moment=scale * width * moment) end diff --git a/src/settings.jl b/src/settings.jl index a5aade8c..ae4b0f7a 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -180,6 +180,8 @@ Solver configuration, used within [`VSMSettings`](@ref). (default `false`) - `flow_curvature`: Add the thin-airfoil pitch-rate moment increment to each section (default `false`) +- `is_with_viscous_drag_correction`: Add the spanwise-flow viscous drag and side + force to each section (default `false`) """ @with_kw mutable struct SolverSettings n_panels::Int64 = 40 @@ -202,6 +204,7 @@ Solver configuration, used within [`VSMSettings`](@ref). 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 + is_with_viscous_drag_correction::Bool=false # spanwise-flow viscous force (Gaunaa 2024) end """ diff --git a/src/solver.jl b/src/solver.jl index c63d44cc..e43f143a 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -85,6 +85,7 @@ end gamma_new::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}) + v_span_dist::MVector{P, T} = zeros(MVector{P, T}) end @with_kw struct BaseResult{P, T} @@ -138,6 +139,8 @@ Main solver structure for the Vortex Step Method.See also: [`solve`](@ref) - `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) +- `is_with_viscous_drag_correction`::Bool = false: Add the spanwise-flow viscous drag + and side force to each section, see: [`spanwise_flow_drag`](@ref) - `reference_point`::MVec3 = [0.0, 0.0, 0.0]: Moment reference point in body frame ## Solution @@ -177,6 +180,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref) is_only_f_and_gamma_output::Bool = false correct_aoa::Bool = false flow_curvature::Bool = false + is_with_viscous_drag_correction::Bool = false reference_point::MVector{3, T} = zeros(MVector{3, T}) # Intermediate results @@ -219,6 +223,7 @@ function Solver(body_aero, settings::VSMSettings) is_only_f_and_gamma_output=ss.calc_only_f_and_gamma, correct_aoa=ss.correct_aoa, flow_curvature=ss.flow_curvature, + is_with_viscous_drag_correction=ss.is_with_viscous_drag_correction, reference_point=reference_point, ) end @@ -405,9 +410,16 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; 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_a_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_a_dist[i]), - cl_dist[i], cd_dist[i], cm_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 @@ -593,7 +605,9 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n body_aero.panels, solver.is_only_f_and_gamma_output; correct_aoa=solver.correct_aoa, - flow_curvature=solver.flow_curvature + flow_curvature=solver.flow_curvature, + is_with_viscous_drag_correction=solver.is_with_viscous_drag_correction, + v_span_dist=solver.lr.v_span_dist, ) # Attach geometric AoA (already computed in calculate_results) to solver.sol if haskey(results, "alpha_geometric") @@ -755,6 +769,10 @@ end 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] end solver.lr.alpha_dist .= atan.(v_normal_dist, v_tangential_dist) @@ -1240,6 +1258,7 @@ function make_dual_shadow(solver::Solver{P, U, Float64}, is_only_f_and_gamma_output = solver.is_only_f_and_gamma_output, correct_aoa = solver.correct_aoa, flow_curvature = solver.flow_curvature, + is_with_viscous_drag_correction = solver.is_with_viscous_drag_correction, 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 aafe3fd8..b77976d4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -53,6 +53,7 @@ function include_selected_tests() 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_viscous_drag_correction.jl") && include("solver/test_viscous_drag_correction.jl") 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") diff --git a/test/solver/test_viscous_drag_correction.jl b/test/solver/test_viscous_drag_correction.jl new file mode 100644 index 00000000..2b6a254a --- /dev/null +++ b/test/solver/test_viscous_drag_correction.jl @@ -0,0 +1,100 @@ +using VortexStepMethod +using VortexStepMethod: spanwise_flow_drag, calc_forces! +using LinearAlgebra +using Test + +@testset "spanwise_flow_drag is the Gaunaa et al. 2024 model" begin + # Gaunaa et al. 2024, J. Phys.: Conf. Ser. 2767 022068, with Re and both + # coefficients built on the speed normal to the span: + # f0 = 0.062 Re^(-1/7), ΔCd = f0 (cos(β)^(-5/7) - 1), C_par = f0 tan(β) cos(β)^(-5/7) + density, mu, chord, speed = 1.225, 1.81e-5, 0.8, 25.0 + for beta in deg2rad.((-40.0, -5.0, 0.0, 10.0, 35.0)) + v_normal = speed * cos(beta) + f0 = 0.062 * (density * v_normal * chord / mu)^(-1 / 7) + correction = spanwise_flow_drag(v_normal, speed * sin(beta), chord, density, mu) + @test correction.delta_cd ≈ f0 * (cos(beta)^(-5 / 7) - 1) atol = 1e-14 + @test correction.c_span ≈ f0 * tan(beta) * cos(beta)^(-5 / 7) atol = 1e-14 + end + @test spanwise_flow_drag(speed, 0.0, chord, density, mu) == (delta_cd=0.0, c_span=0.0) +end + +@testset "Spanwise-flow viscous drag correction in the solver" begin + chord, span, V = 1.0, 8.0, 20.0 + alpha, beta = deg2rad(6.0), deg2rad(12.0) + + 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]) + + solver_off = Solver(body_aero; use_gamma_prev=false) + solver_on = Solver(body_aero; use_gamma_prev=false, + is_with_viscous_drag_correction=true) + va_sideslip = V .* [cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)] + va_straight = V .* [cos(alpha), 0.0, sin(alpha)] + + function force_dist_at(solver, va) + set_va!(body_aero, va) + solve!(solver, body_aero) + return copy(solver.sol.f_body_3D) + end + + @testset "adds the model's drag and spanwise force to each panel" begin + density, mu = solver_on.density, solver_on.mu + for va in (va_sideslip, va_straight) + delta_force = force_dist_at(solver_on, va) .- force_dist_at(solver_off, va) + for (i, panel) in enumerate(body_aero.panels) + v_normal = solver_on.lr.v_a_dist[i] + v_span = solver_on.lr.v_span_dist[i] + cos_beta = v_normal / hypot(v_normal, v_span) + f0 = 0.062 * (density * v_normal * panel.chord / mu)^(-1 / 7) + q_chord_width = 0.5 * density * v_normal^2 * panel.chord * panel.width + span_force = dot(delta_force[:, i], panel.y_airf) + drag_force = norm(delta_force[:, i] .- span_force .* panel.y_airf) + @test span_force ≈ q_chord_width * f0 * (v_span / v_normal) * + cos_beta^(-5 / 7) atol = 1e-12 + @test drag_force ≈ q_chord_width * f0 * (cos_beta^(-5 / 7) - 1) atol = 1e-12 + end + end + end + + @testset "sideslip drives the spanwise flow and raises the drag" begin + delta_force = force_dist_at(solver_on, va_sideslip) .- + force_dist_at(solver_off, va_sideslip) + for (i, panel) in enumerate(body_aero.panels) + @test solver_on.lr.v_span_dist[i] ≈ dot(va_sideslip, panel.y_airf) rtol = 0.05 + @test dot(delta_force[:, i], va_sideslip) > 0 + end + end + + @testset "solve reports the corrected forces" begin + set_va!(body_aero, va_sideslip) + solve!(solver_on, body_aero) + results = solve(solver_on, body_aero) + @test [results["Fx"], results["Fy"], results["Fz"]] ≈ solver_on.sol.force + @test results["F_distribution"] ≈ solver_on.sol.f_body_3D + end + + @testset "linearize reports the corrected forces" begin + y = [va_sideslip; zeros(3)] + results_for(solver) = VortexStepMethod.linearize(solver, body_aero, y; + theta_idxs=nothing, va_idxs=1:3, omega_idxs=4:6)[2] + results_on, results_off = results_for(solver_on), results_for(solver_off) + @test results_on[1:3] ≈ vec(sum(force_dist_at(solver_on, va_sideslip); dims=2)) + @test !(results_on[1:3] ≈ results_off[1:3]) + end + + @testset "calc_forces! stays zero-alloc" begin + set_va!(body_aero, va_sideslip) + solve!(solver_on, body_aero) + calc_forces!(solver_on, body_aero) + @test (@allocated calc_forces!(solver_on, body_aero)) == 0 + end + + @testset "defaults to off" begin + @test Solver(body_aero).is_with_viscous_drag_correction == false + @test VortexStepMethod.SolverSettings().is_with_viscous_drag_correction == false + end +end