diff --git a/src/filament.jl b/src/filament.jl index e2c5a772..c1a5d5d5 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -84,17 +84,19 @@ function velocity_3D_bound_vortex!( @debug "inside core radius" @debug "distance from control point to filament: $(nr1Xr0 / nr0)" - # Project onto core radius - cross3!(r2Xr0, r2, r0) nr0sq = nr0 * nr0 - nr2Xr0 = norm3(r2Xr0) 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!(r1_projXr2_proj, r1_proj, r2_proj) @@ -289,13 +291,13 @@ function velocity_3D_trailing_vortex_semiinfinite!( else r1_proj = work_vectors[4] cross_tmp = work_vectors[5] - # temp = r1/nr1 - Vf + nVfsq = nVf * nVf @inbounds for k in 1:3 - cross_tmp[k] = r1[k]/nr1 - Vf[k] + 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] + + r1_proj[k] = d_r1_Vf * Vf[k] / nVfsq + epsilon * cross_tmp[k] / n_tmp end cross3!(cross_tmp, r1_proj, Vf) diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl index 511af4c3..11996c51 100644 --- a/test/filament/test_bound_filament.jl +++ b/test/filament/test_bound_filament.jl @@ -74,7 +74,7 @@ end filament = BoundFilament{Float64}() reinit!(filament, [0.0, 0.0, 0.0], [1e6, 0.0, 0.0]) control_point = [5e5, 1.0, 0.0] - + velocity_3D_bound_vortex!( induced_velocity, filament, @@ -85,8 +85,8 @@ end ) @test !any(isnan.(induced_velocity)) @test isapprox(induced_velocity[1], 0.0, atol=1e-8) - @test abs(induced_velocity[2]) < 1e-8 - @test isapprox(induced_velocity[3], 0.0) + @test isapprox(induced_velocity[2], 0.0, atol=1e-8) + @test abs(induced_velocity[3]) > 1e-10 end @testset "Point Far from Filament" begin @@ -133,35 +133,35 @@ end @testset "Around Core Radius" begin filament = create_test_filament() delta = 1e-5 - + points = [ [0.5, core_radius_fraction - delta, 0.0], [0.5, core_radius_fraction, 0.0], [0.5, core_radius_fraction + delta, 0.0] ] - + velocities = [zeros(3) for p in points] [ velocity_3D_bound_vortex!(velocities[i], filament, p, gamma, core_radius_fraction, work_vectors) for (i, p) in enumerate(points) ] - + # Check for NaN and finite values @test all(!any(isnan.(v)) for v in velocities) @test all(all(isfinite.(v)) for v in velocities) - + # Check magnitude is maximum at core radius @test norm(velocities[2]) > norm(velocities[1]) @test norm(velocities[2]) > norm(velocities[3]) - + # Check continuity around core radius @test isapprox(velocities[1], velocities[2], rtol=1e-2) - + # Check non-zero velocities @test !all(isapprox.(velocities[1], zeros(3), atol=1e-10)) @test !all(isapprox.(velocities[2], zeros(3), atol=1e-10)) @test !all(isapprox.(velocities[3], zeros(3), atol=1e-10)) - + # Check symmetry v_neg = zeros(3) velocity_3D_bound_vortex!( @@ -174,4 +174,80 @@ end ) @test isapprox(velocities[2], -v_neg) end + + @testset "Velocity is azimuthal (perpendicular to axis and radius)" 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) + end + end + end + + @testset "Solid-body rotation inside core" begin + filament = create_test_filament() + + v_half = zeros(3) + v_quarter = zeros(3) + velocity_3D_bound_vortex!( + v_half, filament, + [0.5, 0.5 * core_radius_fraction, 0.0], + gamma, core_radius_fraction, work_vectors) + velocity_3D_bound_vortex!( + v_quarter, filament, + [0.5, 0.25 * core_radius_fraction, 0.0], + gamma, core_radius_fraction, work_vectors) + + @test isapprox(norm(v_quarter), 0.5 * norm(v_half); rtol=1e-3) + @test isapprox(normalize(v_quarter), normalize(v_half); atol=1e-8) + end + + @testset "Branch 1 / Branch 3 agree at core boundary" begin + filament = create_test_filament() + delta = 0.05 * core_radius_fraction + v_inside = zeros(3) + v_outside = zeros(3) + velocity_3D_bound_vortex!( + v_inside, filament, + [0.5, core_radius_fraction - delta, 0.0], + gamma, core_radius_fraction, work_vectors) + velocity_3D_bound_vortex!( + v_outside, filament, + [0.5, core_radius_fraction + delta, 0.0], + gamma, core_radius_fraction, work_vectors) + + @test isapprox(norm(v_inside), norm(v_outside); rtol=0.1) + @test isapprox(normalize(v_inside), normalize(v_outside); atol=1e-2) + end + + @testset "Matches Rankine vortex profile inside core" begin + L = 1e6 + r_c = 0.01 * L + filament_long = BoundFilament{Float64}() + reinit!(filament_long, [0.0, 0.0, 0.0], [L, 0.0, 0.0]) + + for r in (1.0, 10.0, 100.0, 1000.0) + v = zeros(3) + velocity_3D_bound_vortex!( + v, filament_long, [L/2, r, 0.0], + gamma, 0.01, work_vectors) + + expected_mag = gamma * r / (2π * r_c^2) + + @test isapprox(v[1], 0.0; atol=expected_mag * 1e-4) + @test isapprox(v[2], 0.0; atol=expected_mag * 1e-4) + @test isapprox(abs(v[3]), expected_mag; rtol=1e-3) + @test v[3] > 0 + end + end end \ No newline at end of file diff --git a/test/filament/test_semi_infinite_filament.jl b/test/filament/test_semi_infinite_filament.jl index ab870fb5..43b2149f 100644 --- a/test/filament/test_semi_infinite_filament.jl +++ b/test/filament/test_semi_infinite_filament.jl @@ -118,7 +118,7 @@ end filament = create_test_filament2() vel_pos = zeros(3) vel_neg = zeros(3) - + velocity_3D_trailing_vortex_semiinfinite!( vel_pos, filament, @@ -137,7 +137,42 @@ end filament.vel_mag, work_vectors ) - + @test isapprox(vel_pos, -vel_neg) end + + @testset "Velocity is azimuthal (perpendicular to axis and radius)" begin + filament = create_test_filament2() + direction = filament.direction + + for d in (1e-4, 1e-3, 5e-3, 1e-2, 1e-1) + for phi in (0.0, π/4, π/2, π, -π/3) + p = [0.5, d * cos(phi), d * sin(phi)] + v = zeros(3) + velocity_3D_trailing_vortex_semiinfinite!( + v, filament, direction, p, gamma, + filament.vel_mag, work_vectors) + + r_radial = [0.0, p[2], p[3]] + @test isapprox(dot(v, direction), 0.0; atol=1e-10) + @test isapprox(dot(v, r_radial) / norm(v), 0.0; atol=1e-6) + end + end + end + + @testset "Constant azimuthal direction inside core" begin + filament = create_test_filament2() + v_a = filament.vel_mag + + d_inside = 1e-4 + v1 = zeros(3); v2 = zeros(3) + velocity_3D_trailing_vortex_semiinfinite!( + v1, filament, filament.direction, + [0.5, d_inside, 0.0], gamma, v_a, work_vectors) + velocity_3D_trailing_vortex_semiinfinite!( + v2, filament, filament.direction, + [0.5, 2 * d_inside, 0.0], gamma, v_a, work_vectors) + + @test isapprox(normalize(v2), normalize(v1); atol=1e-8) + end end