diff --git a/CHANGELOG.md b/CHANGELOG.md index 4ca2653a..84eb614b 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 15f636ef..86ae4c1a 100644 --- a/docs/src/private_functions.md +++ b/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! diff --git a/src/filament.jl b/src/filament.jl index c675a860..b4993542 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, 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 @@ -132,7 +70,7 @@ Calculate induced velocity by a trailing vortex filament. - `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,60 +81,67 @@ 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 @@ -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 @@ -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) diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl index 11996c51..397eaf99 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 @@ -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 \ No newline at end of file +end