Give the semi-infinite trailing vortex the bound kernel's linear core, from one shared core definition

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.

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!

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
@@ -173,10 +170,20 @@ 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`.
+
+# 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,
@@ -192,51 +199,71 @@ 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)
+    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)
-    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
+    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
-        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 = 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]
     end
     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
 
 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]

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!, lamb_oseen_core_radius
+using ForwardDiff
 using LinearAlgebra
 using Test
-# using BenchmarkTools
 
 function create_test_filament2()
     x1 = [0.0, 0.0, 0.0]
@@ -13,28 +14,28 @@ 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
+
+"""
+    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 +59,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 +74,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 +83,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 +178,28 @@ end
 
         @test isapprox(normalize(v2), normalize(v1); atol=1e-8)
     end
+
+    @testset "Velocity scales linearly with distance inside core" begin
+        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)
+
+        @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 = 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)
+
+        @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 ≈ velocity_at(1e-5) / 1e-5 rtol = 1e-9
+    end
 end