Bound and trailing vortices share one segment kernel, so the trailing core induces azimuthal velocity

CHANGELOG.md

@@ -20,9 +20,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].
-
-### Fixed
-
+- 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.
 - With `artificial_damping` on, an iteration whose circulation is already smooth no longer
   re-applies the previous iteration's damping correction.
 

docs/src/private_functions.md

@@ -63,6 +63,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!

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,87 +53,24 @@ 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, va, work_vectors)
+    velocity_3D_trailing_vortex!(vel, filament::BoundFilament,
+                                 XVP, gamma, va, 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
 - `gamma`: Vortex strength
 - `va`: 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,71 +81,78 @@ as implemented in KiteAeroDyn".
     va,
     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 / va)
+    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 / va)
-
     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
         @inbounds for k in 1:3
             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
@@ -263,8 +208,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 / va)

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
 
@@ -175,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)
+        va = 1e-4
+        trailing_core_radius = sqrt(4 * ALPHA0 * NU * 0.5 / va)
+        vortices = (
+            (velocity_3D_bound_vortex!, core_radius_fraction, core_radius_fraction),
+            (velocity_3D_trailing_vortex!, va, 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
@@ -250,4 +260,4 @@ end
             @test v[3] > 0
         end
     end
-end
+end