Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
1-Bart-1 merged 4 commits into from
Sep 18, 2026
Merged

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

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
30b0756
Share one straight-segment kernel between bound and trailing vortices
1-Bort-1 Sep 16, 2026
465d199
Test both vortex functions for azimuthal velocity in one testset
1-Bort-1 Sep 17, 2026
c02a2ab
Merge origin/main: carry the va rename into the shared vortex kernel
1-Bort-1 Sep 17, 2026
3dc74fc
Fold the second Unreleased "Fixed" changelog heading into the first
1-Bort-1 Sep 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
1 change: 1 addition & 0 deletions docs/src/private_functions.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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!
Expand Down
152 changes: 48 additions & 104 deletions src/filament.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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, 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
- `gamma`: Vortex strength
- `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!(
Expand All @@ -143,71 +81,78 @@ 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
@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
Expand Down Expand Up @@ -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)
Expand Down
25 changes: 23 additions & 2 deletions test/filament/test_bound_filament.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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

Expand Down Expand Up @@ -250,4 +251,24 @@ end
@test v[3] > 0
end
end
end

@testset "Trailing vortex velocity is azimuthal inside and outside the core" begin

@1-Bort-1 1-Bort-1 Sep 16, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

MINOR: This testset nearly copies 'Velocity is azimuthal (perpendicular to axis and radius)' at line 179, changing only the function called and the core radius. Since both functions now share one kernel, one loop over (function, core radius) pairs would state the invariant once instead of keeping two copies in step (§2, §7).

@1-Bart-1 1-Bart-1 Sep 17, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fix this

@1-Bort-1 1-Bort-1 Sep 17, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed in 465d199: one testset now loops over both vortex functions, each at its own core radius.

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
Loading