From 40e03aedf4e6d6a2bc11156f786739081d854261 Mon Sep 17 00:00:00 2001 From: pivaps Date: Thu, 8 Aug 2024 17:47:52 +0100 Subject: [PATCH] Fixing reflection transmission from plate two media and adding test to compare with one medium case --- .../plane_waves/boundary-condition.jl | 25 +++---- .../plane_waves/reflection-transmission.jl | 12 ++-- .../acoustics-3D/Two-media/planar-symmetry.jl | 67 ++++++++++++++++--- 3 files changed, 75 insertions(+), 29 deletions(-) diff --git a/src/effective_waves/plane_waves/boundary-condition.jl b/src/effective_waves/plane_waves/boundary-condition.jl index 9493c83..19de702 100644 --- a/src/effective_waves/plane_waves/boundary-condition.jl +++ b/src/effective_waves/plane_waves/boundary-condition.jl @@ -252,7 +252,7 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C k_effz_m = -k_effz_p # Needed coefficients - Δ = 2 * k * k0 * ρ * ρ0 * cos(k0*Z) - 1im * ((k0 * ρ)^2 + (k * ρ0)^2 * sin(k0 * Z)) + Δ = 2 * k * k0 * ρ * ρ0 * cos(k0*Z) - 1im * ((k0 * ρ)^2 + (k * ρ0)^2) * sin(k0 * Z) D_p = k0 * ρ * (k0 * ρ + k * ρ0) / Δ D_m = k0 * ρ * (k0 * ρ - k * ρ0) / Δ D_1 = ((k0 * ρ)^2 - (k * ρ0)^2) / 2Δ @@ -264,7 +264,7 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C Ys = spherical_harmonics(basis_order, rθφ[2], rθφ[3]) lm_to_n = lm_to_spherical_harmonic_index - Pr(x::Complex{T}, y::Complex{T}, r::Float64) = exp(2im * (x + y) / Z) * sin((x + y) * (Z / 2 - r)) / (x + y) + Pr(x::Complex{T}, y::Complex{T}, r::T) = exp(1im * (x + y) * Z / 2) * sin((x + y) * (Z / 2 - r)) / (x + y) Bp(F::Array{Complex{T}}, p_or_m::Int) = (4π / (k0^2)) * sum( 1im^(-T(l)) * Ys[lm_to_n(l, m)] * Pr(p_or_m * k_eff, -k0, rs[j]) * F[lm_to_n(l, m), j, p] * nf[j] @@ -274,22 +274,23 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C 1im^(T(l)) * Ys[lm_to_n(l, m)] * Pr(p_or_m * k_eff, k0, rs[j]) * F[lm_to_n(l, m), j, p] * nf[j] for p = 1:size(F, 3), l = 0:basis_order for m = -l:l, j in eachindex(species)) - Iu(F::Array{Complex{T}}, p_or_m::Int) = 4π * sum( - 1im^(-T(dl)) * Ys[lm_to_n(dl, dm)] * exp(1im * (k0 + p_or_m * k_eff) * rs[j]) * (k_eff + p_or_m * k0) * F[lm_to_n(dl, dm), j, p] * nf[j] + Iu(F::Array{Complex{T}}, p_or_m::Int) = (2π / (k0^2)) * sum( + 1im^(-T(dl)) * Ys[lm_to_n(dl, dm)] * exp(1im * (p_or_m * k_eff - k0) * rs[j]) * F[lm_to_n(dl, dm), j, p] * nf[j] / (k_eff - p_or_m * k0) for p = 1:size(F, 3), dl = 0:basis_order for dm = -dl:dl, j in eachindex(species)) - Il(F::Array{Complex{T}}, p_or_m::Int) = 4π * sum( - 1im^(T(dl)) * Ys[lm_to_n(dl, dm)] * exp(1im * (k0 + p_or_m * k_eff) * (Z - rs[j])) * (k_eff - p_or_m * k0) * F[lm_to_n(dl, dm), j, p] * nf[j] + Il(F::Array{Complex{T}}, p_or_m::Int) = (2π / (k0^2)) * sum( + 1im^(T(dl)) * Ys[lm_to_n(dl, dm)] * exp(1im * (p_or_m * k_eff + k0) * (Z - rs[j])) * F[lm_to_n(dl, dm), j, p] * nf[j] / (k_eff + p_or_m * k0) for p = 1:size(F, 3), dl = 0:basis_order for dm = -dl:dl, j in eachindex(species)) - MI11 = D_2 * Bp(F_p, 1) * exp(1im * k0 * Z) + D_1 * Bm(F_p, 1) * exp(-1im * k0 * Z) - Iu(F_p, 1) - MI12 = D_2 * Bp(F_m, -1) * exp(1im * k0 * Z) + D_1 * Bm(F_m, -1) * exp(-1im * k0 * Z) - Iu(F_m, -1) - MI21 = (D_1 * Bp(F_p, 1) + D_2 * Bm(F_p, 1)) * exp(1im * k0 * Z) + Il(F_p, 1) - MI22 = (D_1 * Bp(F_m, -1) + D_2 * Bm(F_m, -1)) * exp(1im * k0 * Z) + Il(F_m, -1) + MI11 = D_2 * Bp(F_p, 1) * exp(1im * k0 * Z) + D_1 * Bm(F_p, 1) * exp(-1im * k0 * Z) + 1im * Iu(F_p, 1) + MI12 = D_2 * Bp(F_m, -1) * exp(1im * k0 * Z) + D_1 * Bm(F_m, -1) * exp(-1im * k0 * Z) + 1im * Iu(F_m, -1) + MI21 = (D_1 * Bp(F_p, 1) + D_2 * Bm(F_p, 1)) * exp(1im * k0 * Z) - 1im * Il(F_p, 1) + MI22 = (D_1 * Bp(F_m, -1) + D_2 * Bm(F_m, -1)) * exp(1im * k0 * Z) - 1im * Il(F_m, -1) - MIs = [MI11 MI12; MI21 MI22] + MIs = [MI11 MI12; + MI21 MI22] - forcing = G * γ0 * [D_p * exp(-1im * k0 * Z), D_m * exp(1im * k0 * Z)] + forcing = -G * γ0 * [D_p * exp(-1im * k0 * Z), D_m * exp(1im * k0 * Z)] α = MIs \ forcing diff --git a/src/effective_waves/plane_waves/reflection-transmission.jl b/src/effective_waves/plane_waves/reflection-transmission.jl index 40a3c46..bd79c00 100644 --- a/src/effective_waves/plane_waves/reflection-transmission.jl +++ b/src/effective_waves/plane_waves/reflection-transmission.jl @@ -167,20 +167,20 @@ function reflection_transmission_coefficients(wavemodes::Vector{E}, psource::Pla k0z = sqrt(k0^2 - k^2 + kz^2) # Needed coefficients - Δ = 2 * k * k0 * ρ * ρ0 * cos(k0 * Z) - 1im * ((k0 * ρ)^2 + (k * ρ0)^2 * sin(k0 * Z)) + Δ = 2 * k * k0 * ρ * ρ0 * cos(k0 * Z) - 1im * ((k0 * ρ)^2 + (k * ρ0)^2) * sin(k0 * Z) D_p = k0 * ρ * (k0 * ρ + k * ρ0) / Δ D_m = k0 * ρ * (k0 * ρ - k * ρ0) / Δ D_0 = 2 * k * k0 * ρ * ρ0 / Δ D_1 = ((k0 * ρ)^2 - (k * ρ0)^2) / 2Δ - Pr(x::Complex{T}, y::Complex{T}, r::T) = exp(2im * (x + y) / Z) * sin((x + y) * (Z / 2 - r)) / (x + y) + Pr(x::Complex{T}, y::Complex{T}, r::T) = exp(1im * (x + y) * Z / 2) * sin((x + y) * (Z / 2 - r)) / (x + y) - Bp(Fp::Array{Complex{T}}, Fm::Array{Complex{T}}) = (4π / (k0 * k0z)) * sum( + Bp(Fp::Array{Complex{T}}, Fm::Array{Complex{T}}) = (4π / (k0^2)) * sum( 1im^(-T(l)) * Ys[lm_to_n(l, m)] * (Pr(k_eff, -k0, rs[j]) * Fp[lm_to_n(l, m), j, p] + Pr(-k_eff, -k0, rs[j]) * Fm[lm_to_n(l, m), j, p]) * nf[j] for p = 1:size(Fp, 3), l = 0:basis_order for m = -l:l, j in eachindex(species)) - Bm(Fp::Array{Complex{T}}, Fm::Array{Complex{T}}) = (4π / (k0 * k0z)) * sum( + Bm(Fp::Array{Complex{T}}, Fm::Array{Complex{T}}) = (4π / (k0^2)) * sum( 1im^(T(l)) * Ys[lm_to_n(l, m)] * (Pr(k_eff, k0, rs[j]) * Fp[lm_to_n(l, m), j, p] + Pr(-k_eff, k0, rs[j]) * Fm[lm_to_n(l, m), j, p]) * nf[j] for p = 1:size(Fp, 3), l = 0:basis_order for m = -l:l, j in eachindex(species)) @@ -193,10 +193,6 @@ function reflection_transmission_coefficients(wavemodes::Vector{E}, psource::Pla D_p * Bp(wavemodes[1].eigenvectors, wavemodes[2].eigenvectors) + D_0 * G) * exp(-1im * k * Z) - # direction_ref = psource.direction - 2 * dot(n,psource.direction) * n - # reflected_wave = PlaneSource(psource.medium; direction = direction_ref, amplitude = Ramp) - # transmitted_wave = PlaneSource(psource.medium; direction = psource.direction, amplitude = Tamp) - return [Ramp, Tamp] else # Unpacking parameters diff --git a/test/acoustics-3D/Two-media/planar-symmetry.jl b/test/acoustics-3D/Two-media/planar-symmetry.jl index f8d2663..5097f07 100644 --- a/test/acoustics-3D/Two-media/planar-symmetry.jl +++ b/test/acoustics-3D/Two-media/planar-symmetry.jl @@ -2,7 +2,7 @@ using EffectiveWaves, Test, LinearAlgebra @testset "Compare one and two media halfspace" begin - # Low frequency test + # Two media test spatial_dim = 3 # Choose the frequency ωs = 0.001:0.5:3.001 @@ -15,7 +15,7 @@ using EffectiveWaves, Test, LinearAlgebra radius1 = 1.0 s1 = Specie( Acoustic(spatial_dim; ρ=10.0, c=10.0), radius1; - volume_fraction=0.05 + volume_fraction=0.2 ) species = [s1] @@ -43,8 +43,56 @@ using EffectiveWaves, Test, LinearAlgebra Reffs1 = [reflection_coefficient(w, source1, material) for w in wavemodes1] Reffs2 = [reflection_coefficient(w, source2, material) for w in wavemodes2] - @test abs(Reffs2[1] - Reffs1[1]) / abs(Reffs1[1]) < 0.05 - @test norm(Reffs2 - Reffs1) / norm(Reffs1) < 0.1 + @test abs(Reffs2[1] - Reffs1[1]) / abs(Reffs1[1]) < 0.01 + @test norm(Reffs2 - Reffs1) / norm(Reffs1) < 0.01 +end + +@testset "Compare one and two media plate" begin + + # Two media test + spatial_dim = 3 + # Choose the frequency + ωs = 0.001:0.5:3.001 + basis_order = 2 + + medium = Acoustic(spatial_dim; ρ=1.001, c=1.001) + medium0 = Acoustic(spatial_dim; ρ=1.0, c=1.0) + + # Choose the species + radius1 = 1.0 + s1 = Specie( + Acoustic(spatial_dim; ρ=0.01, c=0.01), radius1; + volume_fraction=0.2 + ) + species = [s1] + + # Define the halfspace + r = maximum(outer_radius.(species)) + normal = [0.0, 0.0, -1.0] # an outward normal to both surfaces of the plate + width = 200.0; + plate = Plate(normal, width, [0.0, 0.0, width / 2]) + + # define a plane wave sources with slightly different material + source1 = PlaneSource(medium0, [0.0, 0.0, 1.0]) + source2 = PlaneSource(medium, [0.0, 0.0, 1.0]) + + # Calculate the effective wavenumber and wavemode numerically from the general methods + kp_arr = [wavenumbers(ω, medium0, species; tol=1e-6, num_wavenumbers=2, basis_order=basis_order) for ω in ωs] + k_effs = [kps[1] for kps in kp_arr] + + # Define material + material = Material(medium0, plate, species) + + # Compute both wavemodes, for one and two slightly different media + wavemodes1 = [WaveMode(ωs[i], k_effs[i], source1, material; tol=1e-6, basis_order=basis_order) for i in eachindex(ωs)] + wavemodes2 = [WaveMode(ωs[i], k_effs[i], source2, material; tol=1e-6, basis_order=basis_order) for i in eachindex(ωs)] + + # Compute reflection coefficients for both cases + RTs1 = [reflection_transmission_coefficients(w, source1, material) for w in wavemodes1] + RTs2 = [reflection_transmission_coefficients(w, source2, material) for w in wavemodes2] + + @test norm(RTs2[1] - RTs1[1]) / norm(RTs1[1]) < 0.001 + @test sum(norm.(RTs2 - RTs1)) / sum(norm.(RTs1)) < 0.01 end @testset "Compare low frequency halfspace two media" begin @@ -52,17 +100,17 @@ end # Low frequency test spatial_dim = 3 # Choose the frequency - ωs = 0.001:0.02:0.101 + ωs = 0.0001:0.0002:0.0011 basis_order = 2 - medium = Acoustic(spatial_dim; ρ = 2.0, c = 2.5) + medium = Acoustic(spatial_dim; ρ = 3.0, c = 3.0) medium0 = Acoustic(spatial_dim; ρ = 1.0, c = 1.0) # Choose the species radius1 = 1.0 s1 = Specie( Acoustic(spatial_dim; ρ=10.0, c=10.0), radius1; - volume_fraction=0.05 + volume_fraction=0.2 ); species = [s1] @@ -90,7 +138,8 @@ end Reffs = [reflection_coefficient(w, source, material) for w in wavemodes] - @test abs(Reffs[1] - Rlows[1]) / abs(Rlows[1]) < 0.0001 - @test norm(Reffs - Rlows) / norm(Rlows) < 0.005 + @test abs(Reffs[1] - Rlows[1]) / abs(Rlows[1]) < 0.00002 + @test norm(Reffs - Rlows) / norm(Rlows) < 0.0002 end +