From 30b0756d0cebfd449cb28efbc73f9b41186f9a52 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:39:17 +0200 Subject: [PATCH 1/3] Share one straight-segment kernel between bound and trailing vortices MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `velocity_3D_bound_vortex!` and `velocity_3D_trailing_vortex!` were two copies of the same Biot–Savart segment formula that differ only in how the core radius is chosen. #241 corrected the in-core projection in one copy only, so inside its core the trailing vortex still placed the field point along r1 x r0 (azimuthal) and induced a radial velocity. Both now compute their core radius and call `velocity_3D_vortex_segment!`, which projects radially. Outside the core the trailing velocity is bit-identical to before; inside it only the direction changes. The two `@debug` calls in the in-core branch are dropped: in the shared kernel they would have made the trailing in-core path ~2.5x slower. Fixes #333 Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 3 + docs/src/private_functions.md | 1 + src/filament.jl | 152 +++++++++------------------ test/filament/test_bound_filament.jl | 25 ++++- 4 files changed, 75 insertions(+), 106 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4c00c754..74c48de7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,9 @@ - The `VSMSolution` docstring gives `lift_dist`, `drag_dist` and `panel_moment_dist` in the per-unit-span units they hold, [N/m] and [Nm/m], instead of [N] and [Nm]. +- Inside its vortex core, `velocity_3D_trailing_vortex!` induces an azimuthal velocity + instead of a radial one. Only points within the millimetre-scale Oseen core of a + panel's chordwise trailing segment were affected. ## VortexStepMethod v5.1.1 2026-09-12 diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index d7807c42..c9cce3d7 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -62,6 +62,7 @@ panel_loads ```@docs velocity_3D_bound_vortex! velocity_3D_trailing_vortex! +velocity_3D_vortex_segment! velocity_3D_trailing_vortex_semiinfinite! calculate_velocity_induced_bound_2D! calculate_velocity_induced_single_ring_semiinfinite! diff --git a/src/filament.jl b/src/filament.jl index 6ea8bceb..e001b2b3 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -39,10 +39,11 @@ function reinit!(filament::BoundFilament{T}, x1, x2, vec=zeros(MVector{3, T})) w end """ - velocity_3D_bound_vortex(vel, filament::BoundFilament, XVP, - gamma, core_radius_fraction, work_vectors) + velocity_3D_bound_vortex!(vel, filament::BoundFilament, XVP, + gamma, core_radius_fraction, work_vectors) -Calculate induced velocity by a bound vortex filament at a point in space. +Calculate induced velocity by a bound vortex filament at a point in space, with a +core radius of `core_radius_fraction` times the filament length. """ function velocity_3D_bound_vortex!( vel, @@ -52,79 +53,16 @@ function velocity_3D_bound_vortex!( core_radius_fraction, work_vectors ) - r1, r2, r1Xr2, r1Xr0, r2Xr0, r1r2norm, r1_proj, r2_proj, - r1_projXr2_proj, vel_ind_proj = work_vectors - r0 = filament.r0 - nr0 = filament.length - r1 .= XVP .- filament.x1 - r2 .= XVP .- filament.x2 - - epsilon = core_radius_fraction * nr0 - - cross3!(r1Xr0, r1, r0) - - # Check point location relative to filament - nr1Xr0 = norm3(r1Xr0) - if nr1Xr0 / nr0 > epsilon - cross3!(r1Xr2, r1, r2) - nr1 = norm3(r1) - nr2 = norm3(r2) - @inbounds for k in 1:3 - r1r2norm[k] = r1[k]/nr1 - r2[k]/nr2 - end - nr1Xr2 = norm3(r1Xr2) - coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, r1r2norm) - @inbounds for k in 1:3 - vel[k] = coeff * r1Xr2[k] - end - elseif nr1Xr0 / nr0 < 1e-12 * epsilon - vel .= 0.0 - else - @debug "inside core radius" - @debug "distance from control point to filament: $(nr1Xr0 / nr0)" - - nr0sq = nr0 * nr0 - d_r1_r0 = dot3(r1, r0) - d_r2_r0 = dot3(r2, r0) - r_rad = r1Xr0 - @inbounds for k in 1:3 - r_rad[k] = r1[k] - d_r1_r0 * r0[k] / nr0sq - end - nr_rad = norm3(r_rad) - @inbounds for k in 1:3 - r1_proj[k] = d_r1_r0 * r0[k] / nr0sq + - epsilon * r_rad[k] / nr_rad - r2_proj[k] = d_r2_r0 * r0[k] / nr0sq + - epsilon * r_rad[k] / nr_rad - end - cross3!(r1_projXr2_proj, r1_proj, r2_proj) - - nr1pXr2p = norm3(r1_projXr2_proj) - nr1_proj = norm3(r1_proj) - nr2_proj = norm3(r2_proj) - d_sum = 0.0 - @inbounds for k in 1:3 - d_sum += r0[k] * (r1_proj[k]/nr1_proj - - r2_proj[k]/nr2_proj) - end - coeff = (gamma / (4π)) / (nr1pXr2p^2) * d_sum - @inbounds for k in 1:3 - vel_ind_proj[k] = coeff * r1_projXr2_proj[k] - end - - scale = nr1Xr0 / (nr0 * epsilon) - @inbounds for k in 1:3 - vel[k] = scale * vel_ind_proj[k] - end - end - nothing + epsilon = core_radius_fraction * filament.length + velocity_3D_vortex_segment!(vel, filament, XVP, gamma, epsilon, work_vectors) end """ - velocity_3D_trailing_vortex(vel, filament::BoundFilament, - XVP, gamma, v_a, work_vectors) + velocity_3D_trailing_vortex!(vel, filament::BoundFilament, + XVP, gamma, v_a, work_vectors) -Calculate induced velocity by a trailing vortex filament. +Calculate induced velocity by a trailing vortex filament, with a Lamb–Oseen core +radius grown over the axial distance of `XVP` from the filament start. # Arguments - `XVP`: Control point coordinates @@ -132,7 +70,7 @@ Calculate induced velocity by a trailing vortex filament. - `v_a`: Inflow velocity magnitude - work_vectors: preallocated array of intermediate variables -Reference: Rick Damiani et al. "A vortex step method for nonlinear airfoil polar data +Reference: Rick Damiani et al. "A vortex step method for nonlinear airfoil polar data as implemented in KiteAeroDyn". """ @inline function velocity_3D_trailing_vortex!( @@ -143,60 +81,67 @@ as implemented in KiteAeroDyn". v_a, work_vectors ) - r1 = work_vectors[2] - r2 = work_vectors[3] - r_perp = work_vectors[4] - r1Xr2 = work_vectors[5] - r1Xr0 = work_vectors[6] - r2Xr0 = work_vectors[7] - normr1r2 = work_vectors[8] + r1 = work_vectors[1] + r1 .= XVP .- filament.x1 + axial_distance = abs(dot3(r1, filament.r0)) / filament.length + epsilon = sqrt(4 * ALPHA0 * NU * axial_distance / v_a) + velocity_3D_vortex_segment!(vel, filament, XVP, gamma, epsilon, work_vectors) +end +""" + velocity_3D_vortex_segment!(vel, filament::BoundFilament, XVP, + gamma, epsilon, work_vectors) + +Calculate the Biot–Savart velocity induced by a straight vortex segment at `XVP`. +Inside the core radius `epsilon` the velocity is evaluated on the core boundary and +scaled linearly with the distance to the axis. +""" +@inline function velocity_3D_vortex_segment!( + vel, + filament::BoundFilament, + XVP, + gamma, + epsilon, + work_vectors +) + r1, r2, r1Xr2, r1Xr0, r1r2norm, r1_proj, r2_proj = work_vectors r0 = filament.r0 nr0 = filament.length r1 .= XVP .- filament.x1 r2 .= XVP .- filament.x2 - nr0sq = nr0 * nr0 - d_r1_r0 = dot3(r1, r0) - - # Cut-off radius. The perpendicular component has length |r1.r0|/|r0|, so the - # vector itself is only needed inside the core. - epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_r0) / nr0 / v_a) - cross3!(r1Xr0, r1, r0) - - # Check point location relative to filament nr1Xr0 = norm3(r1Xr0) if nr1Xr0 / nr0 > epsilon cross3!(r1Xr2, r1, r2) nr1 = norm3(r1) nr2 = norm3(r2) @inbounds for k in 1:3 - normr1r2[k] = r1[k]/nr1 - r2[k]/nr2 + r1r2norm[k] = r1[k]/nr1 - r2[k]/nr2 end nr1Xr2 = norm3(r1Xr2) - coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, normr1r2) + coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, r1r2norm) @inbounds for k in 1:3 vel[k] = coeff * r1Xr2[k] end elseif nr1Xr0 / nr0 < 1e-12 * epsilon vel .= 0.0 else - # Project onto core radius — reuse r_perp, normr1r2 - r1_proj = r_perp - r2_proj = normr1r2 - cross3!(r2Xr0, r2, r0) - nr2Xr0 = norm3(r2Xr0) + nr0sq = nr0 * nr0 + d_r1_r0 = dot3(r1, r0) d_r2_r0 = dot3(r2, r0) + r_rad = r1Xr0 + @inbounds for k in 1:3 + r_rad[k] = r1[k] - d_r1_r0 * r0[k] / nr0sq + end + nr_rad = norm3(r_rad) @inbounds for k in 1:3 r1_proj[k] = d_r1_r0 * r0[k] / nr0sq + - epsilon * r1Xr0[k] / nr1Xr0 + epsilon * r_rad[k] / nr_rad r2_proj[k] = d_r2_r0 * r0[k] / nr0sq + - epsilon * r2Xr0[k] / nr2Xr0 + epsilon * r_rad[k] / nr_rad end - cross3!(r1Xr2, r1_proj, r2_proj) - nr1Xr2_val = norm3(r1Xr2) nr1_proj = norm3(r1_proj) nr2_proj = norm3(r2_proj) d_sum = 0.0 @@ -204,10 +149,10 @@ as implemented in KiteAeroDyn". d_sum += r0[k] * (r1_proj[k]/nr1_proj - r2_proj[k]/nr2_proj) end - coeff = (gamma / (4π)) / (nr1Xr2_val^2) * d_sum scale = nr1Xr0 / (nr0 * epsilon) + coeff = scale * (gamma / (4π)) / (norm3(r1Xr2)^2) * d_sum @inbounds for k in 1:3 - vel[k] = scale * coeff * r1Xr2[k] + vel[k] = coeff * r1Xr2[k] end end nothing @@ -262,8 +207,7 @@ function velocity_3D_trailing_vortex_semiinfinite!( GAMMA = -GAMMA * filament.filament_direction r1 .= XVP .- filament.x1 - # Core radius. `r_perp` is `(r1.Vf) Vf`, so its length is `|r1.Vf| |Vf|` and - # the vector itself is only needed inside the core. + # Core radius, grown with the axial distance of `XVP` along `Vf`. d_r1_Vf = dot3(r1, Vf) nVf = norm3(Vf) epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_Vf) * nVf / v_a) diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl index 11996c51..3ec65aa9 100644 --- a/test/filament/test_bound_filament.jl +++ b/test/filament/test_bound_filament.jl @@ -1,4 +1,5 @@ -using VortexStepMethod: BoundFilament, velocity_3D_bound_vortex!, reinit! +using VortexStepMethod: BoundFilament, velocity_3D_bound_vortex!, + velocity_3D_trailing_vortex!, reinit!, ALPHA0, NU using LinearAlgebra using Test @@ -250,4 +251,24 @@ end @test v[3] > 0 end end -end \ No newline at end of file + + @testset "Trailing vortex velocity is azimuthal inside and outside the core" begin + filament = create_test_filament() + r0 = [1.0, 0.0, 0.0] + v_a = 1e-4 + core_radius = sqrt(4 * ALPHA0 * NU * 0.5 / v_a) + + for d in (0.25, 0.5, 0.99, 1.01, 2.0) .* core_radius + for phi in (0.0, π/4, π/2, π, -π/3) + p = [0.5, d * cos(phi), d * sin(phi)] + v = zeros(3) + velocity_3D_trailing_vortex!(v, filament, p, gamma, v_a, work_vectors) + + r_radial = [0.0, p[2], p[3]] + @test norm(v) > 1e-3 + @test isapprox(dot(v, r0), 0.0; atol=1e-10) + @test isapprox(dot(v, r_radial), 0.0; atol=1e-10) + end + end + end +end From 465d1995cb8d82dd95abd1c8db110174a19164e5 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Thu, 17 Sep 2026 09:53:32 +0200 Subject: [PATCH 2/3] Test both vortex functions for azimuthal velocity in one testset The bound and trailing azimuthal testsets differed only in the function called and its core radius; loop over both instead. Co-Authored-By: Claude Opus 5 --- test/filament/test_bound_filament.jl | 57 +++++++++++----------------- 1 file changed, 23 insertions(+), 34 deletions(-) diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl index 3ec65aa9..c1ffded7 100644 --- a/test/filament/test_bound_filament.jl +++ b/test/filament/test_bound_filament.jl @@ -176,21 +176,30 @@ end @test isapprox(velocities[2], -v_neg) end - @testset "Velocity is azimuthal (perpendicular to axis and radius)" begin + @testset "Velocity is azimuthal inside and outside the core" begin filament = create_test_filament() - r0 = [1.0, 0.0, 0.0] - - for d in (0.25, 0.5, 0.99, 1.0, 1.01, 2.0) .* core_radius_fraction - for phi in (0.0, π/4, π/2, π, -π/3) - p = [0.5, d * cos(phi), d * sin(phi)] - v = zeros(3) - velocity_3D_bound_vortex!( - v, filament, p, gamma, - core_radius_fraction, work_vectors) - - r_radial = [0.0, p[2], p[3]] - @test isapprox(dot(v, r0), 0.0; atol=1e-10) - @test isapprox(dot(v, r_radial), 0.0; atol=1e-8) + v_a = 1e-4 + trailing_core_radius = sqrt(4 * ALPHA0 * NU * 0.5 / v_a) + vortices = ( + (velocity_3D_bound_vortex!, core_radius_fraction, core_radius_fraction), + (velocity_3D_trailing_vortex!, v_a, trailing_core_radius), + ) + + for (velocity_3D_vortex!, core_parameter, core_radius) in vortices + @testset "$velocity_3D_vortex!" begin + for distance in (0.25, 0.5, 0.99, 1.0, 1.01, 2.0) .* core_radius + for phi in (0.0, π/4, π/2, π, -π/3) + radial = [0.0, distance * cos(phi), distance * sin(phi)] + point = [0.5, 0.0, 0.0] + radial + velocity = zeros(3) + velocity_3D_vortex!(velocity, filament, point, gamma, + core_parameter, work_vectors) + + @test norm(velocity) > 1e-3 + @test isapprox(dot(velocity, filament.r0), 0.0; atol=1e-10) + @test isapprox(dot(velocity, radial), 0.0; atol=1e-10) + end + end end end end @@ -251,24 +260,4 @@ end @test v[3] > 0 end end - - @testset "Trailing vortex velocity is azimuthal inside and outside the core" begin - filament = create_test_filament() - r0 = [1.0, 0.0, 0.0] - v_a = 1e-4 - core_radius = sqrt(4 * ALPHA0 * NU * 0.5 / v_a) - - for d in (0.25, 0.5, 0.99, 1.01, 2.0) .* core_radius - for phi in (0.0, π/4, π/2, π, -π/3) - p = [0.5, d * cos(phi), d * sin(phi)] - v = zeros(3) - velocity_3D_trailing_vortex!(v, filament, p, gamma, v_a, work_vectors) - - r_radial = [0.0, p[2], p[3]] - @test norm(v) > 1e-3 - @test isapprox(dot(v, r0), 0.0; atol=1e-10) - @test isapprox(dot(v, r_radial), 0.0; atol=1e-10) - end - end - end end From 3dc74fc7b902a1dc5c35117376405f0451c850ca Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Thu, 17 Sep 2026 10:50:18 +0200 Subject: [PATCH 3/3] Fold the second Unreleased "Fixed" changelog heading into the first Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 3 --- 1 file changed, 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5a83f2c3..84eb614b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,9 +23,6 @@ - Inside its vortex core, `velocity_3D_trailing_vortex!` induces an azimuthal velocity instead of a radial one. Only points within the millimetre-scale Oseen core of a panel's chordwise trailing segment were affected. - -### Fixed - - With `artificial_damping` on, an iteration whose circulation is already smooth no longer re-applies the previous iteration's damping correction.