From 0008f2ac160c30aaba996a7ae5d39992dfb24491 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 30 Sep 2026 15:17:52 +0200 Subject: [PATCH 1/2] Give the semi-infinite trailing vortex the bound kernel's linear core velocity_3D_trailing_vortex_semiinfinite! now uses the same zero-core guard and closed-form core as velocity_3D_vortex_segment!: zero at its start point rather than NaN, and inside the core the velocity scales linearly with the distance to the axis, so ForwardDiff sees its slope there. Co-Authored-By: Claude Opus 5.5 --- ...elocity-3d-trailing-vortex-semiinfinite.md | 7 ++ src/filament.jl | 49 +++++------ test/filament/test_semi_infinite_filament.jl | 81 +++++++++++++------ 3 files changed, 83 insertions(+), 54 deletions(-) create mode 100644 changelog.d/399-velocity-3d-trailing-vortex-semiinfinite.md diff --git a/changelog.d/399-velocity-3d-trailing-vortex-semiinfinite.md b/changelog.d/399-velocity-3d-trailing-vortex-semiinfinite.md new file mode 100644 index 00000000..6942d9cc --- /dev/null +++ b/changelog.d/399-velocity-3d-trailing-vortex-semiinfinite.md @@ -0,0 +1,7 @@ +### Changed + +- Inside the core of a semi-infinite trailing vortex the induced velocity scales linearly with the distance to the axis, as it does for bound and finite trailing filaments, rather than holding the core-boundary value. + +### Fixed + +- A semi-infinite trailing vortex evaluated at its start point returns zero rather than NaN, and ForwardDiff sees the velocity's slope on its axis rather than zero. diff --git a/src/filament.jl b/src/filament.jl index d79f612d..4bae7324 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -173,10 +173,13 @@ function reinit!(filament::SemiInfiniteFilament{T}, x1::AbstractVector, end """ - velocity_3D_trailing_vortex_semiinfinite(filament::SemiInfiniteFilament, - Vf, XVP, GAMMA, va, work_vectors) + velocity_3D_trailing_vortex_semiinfinite!(vel, filament::SemiInfiniteFilament, + Vf, XVP, GAMMA, va, work_vectors) -Calculate induced velocity by a semi-infinite trailing vortex filament. +Calculate the velocity induced at `XVP` by a semi-infinite trailing vortex filament along +`Vf`, with a Lamb–Oseen core radius grown over the axial distance of `XVP` from `x1`. +Inside the core the velocity is scaled linearly with the distance to the axis. Without a +core, it is zero within 1e-12 of the axis, relative to the distance from `x1`. """ function velocity_3D_trailing_vortex_semiinfinite!( vel, @@ -192,51 +195,37 @@ function velocity_3D_trailing_vortex_semiinfinite!( GAMMA = -GAMMA * filament.filament_direction r1 .= XVP .- filament.x1 - # 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 / va) cross3!(r1XVf, r1, Vf) - nr1XVf = norm3(r1XVf) + axis_distance = nr1XVf / nVf nr1 = norm3(r1) - if nr1XVf / nVf > epsilon - K = GAMMA / (4π) / (nr1XVf^2) * (1 + d_r1_Vf / nr1) - @inbounds for k in 1:3 - vel[k] = K * r1XVf[k] - end - elseif nr1XVf / nVf < 1e-12 * epsilon + rounding = 1e-12 * nr1 + if epsilon <= rounding && axis_distance <= rounding vel .= 0.0 + return nothing + elseif axis_distance > epsilon + K = GAMMA / (4π) / (nr1XVf^2) * (1 + d_r1_Vf / nr1) else - r1_proj = work_vectors[4] - cross_tmp = work_vectors[5] - nVfsq = nVf * nVf - @inbounds for k in 1:3 - cross_tmp[k] = r1[k] - d_r1_Vf * Vf[k] / nVfsq - end - n_tmp = norm3(cross_tmp) - @inbounds for k in 1:3 - r1_proj[k] = d_r1_Vf * Vf[k] / nVfsq + - epsilon * cross_tmp[k] / n_tmp - end - cross3!(cross_tmp, r1_proj, Vf) - K = GAMMA / (4π) / (norm3(cross_tmp)^2) * - (1 + dot3(r1_proj, Vf) / norm3(r1_proj)) - @inbounds for k in 1:3 - vel[k] = K * cross_tmp[k] - end + K = GAMMA / (4π) / (nVf * epsilon)^2 * + (1 + d_r1_Vf / sqrt(d_r1_Vf^2 / nVf^2 + epsilon^2)) + end + @inbounds for k in 1:3 + vel[k] = K * r1XVf[k] end nothing end - """ cross3!(result::AbstractVector{T}, a::AbstractVector{T}, b::AbstractVector{T}) where T Compute cross product of 3D vectors in-place. """ -@inline function cross3!(result::AbstractVector{T}, a::AbstractVector{T}, b::AbstractVector{T}) where T +@inline function cross3!(result::AbstractVector{T}, a::AbstractVector{T}, + b::AbstractVector{T}) where T x = a[2]*b[3] - a[3]*b[2] y = a[3]*b[1] - a[1]*b[3] z = a[1]*b[2] - a[2]*b[1] diff --git a/test/filament/test_semi_infinite_filament.jl b/test/filament/test_semi_infinite_filament.jl index 11014cae..e593746d 100644 --- a/test/filament/test_semi_infinite_filament.jl +++ b/test/filament/test_semi_infinite_filament.jl @@ -1,7 +1,8 @@ -using VortexStepMethod: SemiInfiniteFilament, velocity_3D_trailing_vortex_semiinfinite!, reinit! +using VortexStepMethod: SemiInfiniteFilament, velocity_3D_trailing_vortex_semiinfinite!, + reinit!, ALPHA0, NU +using ForwardDiff using LinearAlgebra using Test -# using BenchmarkTools function create_test_filament2() x1 = [0.0, 0.0, 0.0] @@ -13,28 +14,36 @@ function create_test_filament2() return filament end -function analytical_solution(control_point, gamma, x1, direction, filament_direction, va) +function analytical_solution(control_point, gamma, x1, direction, filament_direction) gamma = -gamma # Sign convention difference r1 = control_point - x1 r1_cross_direction = cross(r1, direction) - r_perp = dot(r1, direction) * direction - - alpha0 = 1.25643 - nu = 1.48e-5 - epsilon = sqrt(4 * alpha0 * nu * norm(r_perp) / va) - - if norm(r1_cross_direction) > epsilon - K = (gamma / (4π * norm(r1_cross_direction)^2)) * - (1 + dot(r1, direction) / norm(r1)) - return K * r1_cross_direction * filament_direction - else - r1_proj = dot(r1, direction) * direction + epsilon * - (r1/norm(r1) - direction) / norm(r1/norm(r1) - direction) - r1_cross_direction_proj = cross(r1_proj, direction) - K_proj = (gamma / (4π * norm(r1_cross_direction_proj)^2)) * - (1 + dot(r1_proj, direction) / norm(r1_proj)) - return K_proj * r1_cross_direction_proj * filament_direction - end + K = (gamma / (4π * norm(r1_cross_direction)^2)) * (1 + dot(r1, direction) / norm(r1)) + return K * r1_cross_direction * filament_direction +end + +""" + core_radius(axial_distance, va) + +Lamb–Oseen core radius [m] of a trailing vortex `axial_distance` [m] downstream of its +start. +""" +core_radius(axial_distance, va) = sqrt(4 * ALPHA0 * NU * axial_distance / va) + +""" + off_axis_velocity(offset, gamma) + +z velocity [m/s] induced by the unit-speed trailing filament from the origin along x at +the point `offset` [m] off its axis, half a metre downstream. +""" +function off_axis_velocity(offset, gamma) + T = typeof(offset) + filament = SemiInfiniteFilament{T}() + reinit!(filament, zeros(T, 3), T[1, 0, 0], one(T), 1) + velocity = zeros(T, 3) + velocity_3D_trailing_vortex_semiinfinite!(velocity, filament, filament.direction, + [0.5, offset, 0.0], gamma, filament.va, ntuple(_ -> zeros(T, 3), 10)) + return velocity[3] end @testset "SemiInfiniteFilament Tests" begin @@ -58,7 +67,7 @@ end analytical = analytical_solution( control_point, gamma, filament.x1, filament.direction, - filament.filament_direction, filament.va + filament.filament_direction ) @test isapprox(induced_velocity, analytical, rtol=1e-6) @@ -73,7 +82,6 @@ end ] induced_velocity = zeros(3) - # Filament start point is singular in the current implementation. velocity_3D_trailing_vortex_semiinfinite!( induced_velocity, filament, @@ -83,7 +91,7 @@ end filament.va, work_vectors ) - @test all(isnan.(induced_velocity)) + @test induced_velocity == zeros(3) for point in test_points velocity_3D_trailing_vortex_semiinfinite!( @@ -178,4 +186,29 @@ end @test isapprox(normalize(v2), normalize(v1); atol=1e-8) end + + @testset "Velocity scales linearly with distance inside core" begin + epsilon = core_radius(0.5, 1.0) + v_half = off_axis_velocity(0.5 * epsilon, gamma) + v_quarter = off_axis_velocity(0.25 * epsilon, gamma) + + @test v_quarter ≈ 0.5 * v_half rtol = 1e-12 + @test off_axis_velocity(-0.25 * epsilon, gamma) ≈ -v_quarter rtol = 1e-12 + end + + @testset "Velocity is continuous at the core boundary" begin + epsilon = core_radius(0.5, 1.0) + v_inside = off_axis_velocity(epsilon * (1 - 1e-9), gamma) + v_outside = off_axis_velocity(epsilon * (1 + 1e-9), gamma) + + @test v_inside ≈ v_outside rtol = 1e-6 + end + + @testset "ForwardDiff sees the core's slope on the axis" begin + velocity_at(offset) = off_axis_velocity(offset, gamma) + slope = ForwardDiff.derivative(velocity_at, 0.0) + + @test slope != 0 + @test slope ≈ velocity_at(1e-5) / 1e-5 rtol = 1e-9 + end end From 32243fa6f1359a9592bb0b6bfa8171eb3fce2e22 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 30 Sep 2026 16:30:10 +0200 Subject: [PATCH 2/2] Share the core radius, zero-core guard and core terms between the filament kernels MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Both velocity_3D_vortex_segment! and velocity_3D_trailing_vortex_semiinfinite! now assemble their inside-core velocity from core_end_term and core_coefficient and guard the zero-core case with on_axis_without_core. lamb_oseen_core_radius replaces the two copies of the Lamb–Oseen formula and the test's own copy. The semi-infinite docstring lists its arguments and units. Co-Authored-By: Claude Opus 5.5 --- docs/src/private_functions.md | 4 ++ src/filament.jl | 64 ++++++++++++++++---- test/filament/test_semi_infinite_filament.jl | 15 +---- 3 files changed, 58 insertions(+), 25 deletions(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index c70033dc..5bbd818f 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -79,6 +79,10 @@ velocity_3D_bound_vortex! velocity_3D_trailing_vortex! velocity_3D_vortex_segment! velocity_3D_trailing_vortex_semiinfinite! +lamb_oseen_core_radius +on_axis_without_core +core_end_term +core_coefficient calculate_velocity_induced_bound_2D! calculate_velocity_induced_single_ring_semiinfinite! cross3! diff --git a/src/filament.jl b/src/filament.jl index 4bae7324..2288f3d2 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -84,7 +84,7 @@ as implemented in KiteAeroDyn". r1 = work_vectors[1] r1 .= XVP .- filament.x1 axial_distance = abs(dot3(r1, filament.r0)) / filament.length - epsilon = sqrt(4 * ALPHA0 * NU * axial_distance / va) + epsilon = lamb_oseen_core_radius(axial_distance, va) velocity_3D_vortex_segment!(vel, filament, XVP, gamma, epsilon, work_vectors) end @@ -113,8 +113,7 @@ scaled linearly with the distance to the axis. Without a core, it is zero within cross3!(r1Xr0, r1, r0) axis_distance = norm3(r1Xr0) / nr0 - rounding = 1e-12 * norm3(r1) - if epsilon <= rounding && axis_distance <= rounding + if on_axis_without_core(axis_distance, epsilon, norm3(r1)) vel .= 0.0 elseif axis_distance > epsilon cross3!(r1Xr2, r1, r2) @@ -130,11 +129,9 @@ scaled linearly with the distance to the axis. Without a core, it is zero within end else nr0sq = nr0 * nr0 - d_r1_r0 = dot3(r1, r0) - d_r2_r0 = dot3(r2, r0) - d_sum = d_r1_r0 / sqrt(d_r1_r0^2 / nr0sq + epsilon^2) - - d_r2_r0 / sqrt(d_r2_r0^2 / nr0sq + epsilon^2) - coeff = -(gamma / (4π)) * d_sum / (epsilon^2 * nr0sq) + end_terms = core_end_term(dot3(r1, r0), nr0sq, epsilon) - + core_end_term(dot3(r2, r0), nr0sq, epsilon) + coeff = -core_coefficient(gamma, end_terms, nr0sq, epsilon) @inbounds for k in 1:3 vel[k] = coeff * r1Xr0[k] end @@ -180,6 +177,13 @@ Calculate the velocity induced at `XVP` by a semi-infinite trailing vortex filam `Vf`, with a Lamb–Oseen core radius grown over the axial distance of `XVP` from `x1`. Inside the core the velocity is scaled linearly with the distance to the axis. Without a core, it is zero within 1e-12 of the axis, relative to the distance from `x1`. + +# Arguments +- `Vf`: unit direction of the filament [-] +- `XVP`: evaluation point [m] +- `GAMMA`: vortex strength [m²/s] +- `va`: apparent wind speed [m/s] +- `work_vectors`: preallocated 3-vectors for intermediate results """ function velocity_3D_trailing_vortex_semiinfinite!( vel, @@ -197,21 +201,21 @@ function velocity_3D_trailing_vortex_semiinfinite!( d_r1_Vf = dot3(r1, Vf) nVf = norm3(Vf) - epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_Vf) * nVf / va) + epsilon = lamb_oseen_core_radius(abs(d_r1_Vf) / nVf, va) cross3!(r1XVf, r1, Vf) nr1XVf = norm3(r1XVf) axis_distance = nr1XVf / nVf nr1 = norm3(r1) - rounding = 1e-12 * nr1 - if epsilon <= rounding && axis_distance <= rounding + if on_axis_without_core(axis_distance, epsilon, nr1) vel .= 0.0 return nothing elseif axis_distance > epsilon K = GAMMA / (4π) / (nr1XVf^2) * (1 + d_r1_Vf / nr1) else - K = GAMMA / (4π) / (nVf * epsilon)^2 * - (1 + d_r1_Vf / sqrt(d_r1_Vf^2 / nVf^2 + epsilon^2)) + nVfsq = nVf * nVf + K = core_coefficient(GAMMA, 1 + core_end_term(d_r1_Vf, nVfsq, epsilon), nVfsq, + epsilon) end @inbounds for k in 1:3 vel[k] = K * r1XVf[k] @@ -219,6 +223,40 @@ function velocity_3D_trailing_vortex_semiinfinite!( nothing end +""" + lamb_oseen_core_radius(axial_distance, va) + +Lamb–Oseen core radius [m] of a trailing vortex `axial_distance` [m] downstream of its +start, in an apparent wind of `va` [m/s]. +""" +lamb_oseen_core_radius(axial_distance, va) = sqrt(4 * ALPHA0 * NU * axial_distance / va) + +""" + on_axis_without_core(axis_distance, epsilon, nr1) + +True when both the distance to the filament axis and the core radius `epsilon` are +within 1e-12 of `nr1`, the distance from the filament start. +""" +@inline on_axis_without_core(axis_distance, epsilon, nr1) = + max(axis_distance, epsilon) <= 1e-12 * nr1 + +""" + core_end_term(d, nr0sq, epsilon) + +Contribution of one filament end to the core velocity, for `d` the dot product of the +end-to-point vector with the axis vector of squared length `nr0sq`. +""" +@inline core_end_term(d, nr0sq, epsilon) = d / sqrt(d^2 / nr0sq + epsilon^2) + +""" + core_coefficient(gamma, end_terms, nr0sq, epsilon) + +Factor on the cross product of the point vector and the axis vector that gives the +velocity inside a core of radius `epsilon`, linear in the distance to the axis. +""" +@inline core_coefficient(gamma, end_terms, nr0sq, epsilon) = + (gamma / (4π)) * end_terms / (epsilon^2 * nr0sq) + """ cross3!(result::AbstractVector{T}, a::AbstractVector{T}, b::AbstractVector{T}) where T diff --git a/test/filament/test_semi_infinite_filament.jl b/test/filament/test_semi_infinite_filament.jl index e593746d..0c6e4fb2 100644 --- a/test/filament/test_semi_infinite_filament.jl +++ b/test/filament/test_semi_infinite_filament.jl @@ -1,5 +1,5 @@ using VortexStepMethod: SemiInfiniteFilament, velocity_3D_trailing_vortex_semiinfinite!, - reinit!, ALPHA0, NU + reinit!, lamb_oseen_core_radius using ForwardDiff using LinearAlgebra using Test @@ -22,14 +22,6 @@ function analytical_solution(control_point, gamma, x1, direction, filament_direc return K * r1_cross_direction * filament_direction end -""" - core_radius(axial_distance, va) - -Lamb–Oseen core radius [m] of a trailing vortex `axial_distance` [m] downstream of its -start. -""" -core_radius(axial_distance, va) = sqrt(4 * ALPHA0 * NU * axial_distance / va) - """ off_axis_velocity(offset, gamma) @@ -188,7 +180,7 @@ end end @testset "Velocity scales linearly with distance inside core" begin - epsilon = core_radius(0.5, 1.0) + epsilon = lamb_oseen_core_radius(0.5, 1.0) v_half = off_axis_velocity(0.5 * epsilon, gamma) v_quarter = off_axis_velocity(0.25 * epsilon, gamma) @@ -197,7 +189,7 @@ end end @testset "Velocity is continuous at the core boundary" begin - epsilon = core_radius(0.5, 1.0) + epsilon = lamb_oseen_core_radius(0.5, 1.0) v_inside = off_axis_velocity(epsilon * (1 - 1e-9), gamma) v_outside = off_axis_velocity(epsilon * (1 + 1e-9), gamma) @@ -208,7 +200,6 @@ end velocity_at(offset) = off_axis_velocity(offset, gamma) slope = ForwardDiff.derivative(velocity_at, 0.0) - @test slope != 0 @test slope ≈ velocity_at(1e-5) / 1e-5 rtol = 1e-9 end end