From c6ad0251604b742097afebff2e86ca5d439786f0 Mon Sep 17 00:00:00 2001 From: arturgower Date: Wed, 10 Mar 2021 18:53:07 +0000 Subject: [PATCH] wrote scattering from effective plate. Still needs testing. --- Project.toml | 4 +- docs/src/manual/reflection.md | 55 ++++++------------- src/EffectiveWaves.jl | 2 +- .../planewaves/boundary_condition.jl | 26 ++++----- ...jl => material_scattering_coefficients.jl} | 55 +++++++++++++++++++ src/acoustics/export.jl | 2 +- src/acoustics/low_frequency.jl | 21 +++---- src/effective_wave/wavemode.jl | 2 +- 8 files changed, 98 insertions(+), 69 deletions(-) rename src/acoustics/effective_wave/planewaves/{reflection_coefficient.jl => material_scattering_coefficients.jl} (72%) diff --git a/Project.toml b/Project.toml index d200393c..6911ba16 100644 --- a/Project.toml +++ b/Project.toml @@ -19,7 +19,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [compat] -ApproxFun = "0.9, 0.10, ^0.11" +ApproxFun = "0.9, 0.10, ^0.11,^0.12" GSL = "0.6, 0.7, 0.8, 0.9, 1.0" MDBM = "0.1" MultipleScattering = "^0.1.5" @@ -28,7 +28,7 @@ Optim = "0.22, 1" Polynomials = "0.8" RecipesBase = "1.0" Reexport = "0.2" -SpecialFunctions = "0.7, 0.8, 0.9, 10.0" +SpecialFunctions = "0.3, 0.7, 0.8, 0.9, 0.10, 1.0" StaticArrays = "0.8,0.9,0.10,0.11,0.12" julia = "1" diff --git a/docs/src/manual/reflection.md b/docs/src/manual/reflection.md index 249dfdfe..281f091b 100644 --- a/docs/src/manual/reflection.md +++ b/docs/src/manual/reflection.md @@ -119,30 +119,32 @@ s1 = Specie( ); radius2 = 0.002 s2 = Specie( - Acoustic(spatial_dim; ρ=0.2, c=4.1), radius2; + Acoustic(spatial_dim; ρ=0.00001, c=0.1), radius2; + #Acoustic(spatial_dim; ρ=0.2, c=4.1), radius2; volume_fraction=0.15 ); species = [s1,s2] +species = [s2] # Choose the frequency ω = 1e-2 k = ω / real(medium.c) -# Check that we are indeed in a low frequency limit -k * radius2 < 1e-4 - +# For the limit of low frequencies we can define eff_medium = effective_medium(medium, species) -normal = [0.0,0.0,-1.0] # an outward normal to the surface - -# Define the material region -material = Material(Halfspace(normal),species) +# Define a plate +normal = [0.0,0.0,-1.0] # an outward normal to both surfaces of the plate +width = 1.0 # plate width +plate1 = Plate(normal,width) # define a plane wave source travelling at a 45 degree angle in relation to the material source = PlaneSource(medium, [cos(pi/4.0),0.0,sin(pi/4.0)]) -R = reflection_coefficient(ω,source, eff_medium, material.shape) +amps = planewave_amplitudes(ω, source, eff_medium, plate1); +Ramp = amps[1] +Tamp = amps[2] # output @@ -155,41 +157,18 @@ Using only one plane wave mode we can calculate both reflection and transmission ```julia 2 k_effs = wavenumbers(ω, medium, species; tol = 1e-6, num_wavenumbers = 1, basis_order = 1) - -# Define a plate -normal = [0.0,0.0,-1.0] # an outward normal to both surfaces of the plate -width = 1.0 # plate width - -# Define the material region -material = Material(Plate(normal,width),species) - -basis_order = 1; -kws = Dict(:basis_order => basis_order) - -MM = eigensystem(ω, source, material; kws...) - -# calculate eigenvectors -using LinearAlgebra - k_eff = k_effs[1] -MM_svd = svd(MM(k_eff)) -MM_svd.S[end] -v1 = MM_svd.V[:,end] - -k_eff = -k_effs[1] -MM_svd = svd(MM(k_eff)) -MM_svd.S[end] -v2 = MM_svd.V[:,end] - +abs(k_eff - ω / eff_medium.c) < 1e-12 - -vecs = eigenvectors(ω, k_effs[1], source, material; kws...) - +plate = Plate(normal,width) +material = Material(plate,species) # Calculate the wavemode for the first wavenumber -wave1 = WaveMode(ω, k_effs[1], source, material; tol = 1e-6, basis_order = 1) +# the WaveMode function calculates the types of waves and solves the needed boundary conditions +wavemodes = WaveMode(ω, k_eff, source, material; tol = 1e-6, basis_order = 1); +material_scattering_coefficients(wavemodes, source, material) ``` Currently implementing... formulas from [Gower & Kristensson 2020](https://arxiv.org/pdf/2010.00934.pdf). diff --git a/src/EffectiveWaves.jl b/src/EffectiveWaves.jl index b6b1d158..596f3246 100644 --- a/src/EffectiveWaves.jl +++ b/src/EffectiveWaves.jl @@ -20,7 +20,7 @@ export dispersion_equation, dispersion_complex, eigensystem # supplies a matrix export scattering_amplitudes_average -export reflection_coefficient, reflection_coefficients +export reflection_coefficient, reflection_coefficients, planewave_amplitudes export effective_medium # List of shorthand for some materials diff --git a/src/acoustics/effective_wave/planewaves/boundary_condition.jl b/src/acoustics/effective_wave/planewaves/boundary_condition.jl index 49dbdbdb..0c88e9f4 100644 --- a/src/acoustics/effective_wave/planewaves/boundary_condition.jl +++ b/src/acoustics/effective_wave/planewaves/boundary_condition.jl @@ -43,23 +43,23 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C scale_number_density = one(T) - one(T) / material.numberofparticles - # First we calculate the innerward pointing normal + # First we calculate the outward pointing normal n = material.shape.normal; - n = n .* sign(real(dot(n,psource.direction))); + n = - n .* sign(real(dot(n,psource.direction))); # calculate the distances from the origin to the two faces of the plate - Z0 = dot(material.shape.normal, material.shape.origin) - Z1 = Z0 - material.shape.width / T(2) - Z2 = Z0 + material.shape.width / T(2) + Z0 = dot(-n, material.shape.origin) + Z1 = Z0 - material.shape.width / 2 + Z2 = Z0 + material.shape.width / 2 # next the components of the wavenumbers in the direction of the inward normal - direction1 = transmission_direction(k_eff, (ω / psource.medium.c) * psource.direction, -n) - direction2 = transmission_direction(- k_eff, (ω / psource.medium.c) * psource.direction, -n) + direction1 = transmission_direction(k_eff, (ω / psource.medium.c) * psource.direction, n) + direction2 = transmission_direction(- k_eff, (ω / psource.medium.c) * psource.direction, n) k = (ω / psource.medium.c) - kz = k * dot(conj(n), psource.direction) - k1_effz = k_eff * dot(conj(n), direction1) - k2_effz = k_eff * dot(conj(n), direction2) + kz = k * dot(-conj(n), psource.direction) + k1_effz = k_eff * dot(-conj(n), direction1) + k2_effz = k_eff * dot(-conj(n), direction2) rθφ = cartesian_to_radial_coordinates(psource.direction); @@ -67,14 +67,14 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C lm_to_n = lm_to_spherical_harmonic_index I1(F::Array{Complex{T}}, k_effz::Complex{T}) = sum( - - F[lm_to_n(dl,dm),j,1] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(1)^dl * + - F[lm_to_n(dl,dm),j,1] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(-1)^dl * exp(im * (k_effz - kz)*(Z1 + outer_radius(material.species[j]))) * scale_number_density * number_density(material.species[j]) / (k*kz * (kz - k_effz)) for dl = 0:basis_order for dm = -dl:dl, j in eachindex(material.species)) I2(F::Array{Complex{T}}, k_effz::Complex{T}) = sum( - - F[lm_to_n(dl,dm),j,1] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(1)^dm * - exp(im * (k_effz - kz)*(Z2 - outer_radius(material.species[j]))) * + - F[lm_to_n(dl,dm),j,1] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(-1)^dm * + exp(im * (k_effz + kz)*(Z2 - outer_radius(material.species[j]))) * scale_number_density * number_density(material.species[j]) / (k*kz * (kz + k_effz)) for dl = 0:basis_order for dm = -dl:dl, j in eachindex(material.species)) diff --git a/src/acoustics/effective_wave/planewaves/reflection_coefficient.jl b/src/acoustics/effective_wave/planewaves/material_scattering_coefficients.jl similarity index 72% rename from src/acoustics/effective_wave/planewaves/reflection_coefficient.jl rename to src/acoustics/effective_wave/planewaves/material_scattering_coefficients.jl index f3016adf..c05a20c4 100644 --- a/src/acoustics/effective_wave/planewaves/reflection_coefficient.jl +++ b/src/acoustics/effective_wave/planewaves/material_scattering_coefficients.jl @@ -30,6 +30,61 @@ function reflection_coefficient(ω::T, wave_eff::EffectivePlaneWaveMode{T}, psou return R end +function material_scattering_coefficients(wavemodes::Vector{E}, psource::PlaneSource{T,3,1,Acoustic{T,3}}, material::Material{3,Plate{T,3}}) where {T<:AbstractFloat,Dim, E<:EffectivePlaneWaveMode{T,Dim}} + + # Unpacking parameters + k = real(wavemodes[1].ω / psource.medium.c) + + species = material.species + S = length(species) + rs = outer_radius.(species) + + ho = maximum(w.basis_order for w in wavemodes) + ls, ms = spherical_harmonics_indices(ho) + + # make the normal outward pointing + plate = material.shape + n = plate.normal / norm(plate.normal) + if real(dot(n,psource.direction)) > 0 + n = -n + end + + rθφ = cartesian_to_radial_coordinates(psource.direction) + Ys = spherical_harmonics(ho, rθφ[2], rθφ[3]); + + direction_ref = psource.direction - 2 * dot(n,psource.direction) * n + rθφ = cartesian_to_radial_coordinates(direction_ref) + Yrefs = spherical_harmonics(ho, rθφ[2], rθφ[3]); + + Z0 = dot(-n,plate.origin - psource.position) + Z1 = Z0 - plate.width / 2 + Z2 = Z0 + plate.width / 2 + + kcos_in = k * dot(- conj(n), psource.direction) + + RTs = map(wavemodes) do w + kcos_eff = w.wavenumber * dot(- conj(n), w.direction) + + Rp = sum( + number_density(species[i[2]]) * w.eigenvectors[i] * 2pi * (1.0im)^(ls[i[1]]-1) * + Yrefs[i[1]] * (exp(im*(kcos_eff + kcos_in)*(Z2 - rs[i[2]])) - exp(im*(kcos_eff + kcos_in)*(Z1 + rs[i[2]]))) / + ((kcos_in + kcos_eff) * k * kcos_in) + for i in CartesianIndices(w.eigenvectors)) + + Tp = sum( + number_density(species[i[2]]) * w.eigenvectors[i] * 2pi * (-1.0)^ls[i[1]] * (1.0im)^(ls[i[1]]+1) * + Ys[i[1]] * (exp(im*(kcos_eff - kcos_in)*(Z2 - rs[i[2]])) - exp(im*(kcos_eff - kcos_in)*(Z1 + rs[i[2]]))) / + ((kcos_in - kcos_eff) * k * kcos_in) + for i in CartesianIndices(w.eigenvectors)) + + [Rp, Tp] + end + + Ramp, Tamp = sum(RTs) + [0.0,1.0] + + return [Ramp, Tamp] +end + "The average reflection coefficient" function wienerhopf_reflection_coefficient(ω::T, psource::PlaneSource{T,2,1,Acoustic{T,2}}, material::Material{2,Halfspace{T,2}}; tol::T = T(1e-7), diff --git a/src/acoustics/export.jl b/src/acoustics/export.jl index 2fcae28b..9e8b6643 100644 --- a/src/acoustics/export.jl +++ b/src/acoustics/export.jl @@ -11,7 +11,7 @@ include("low_frequency.jl") include("effective_wave/asymptotics.jl") include("effective_wave/planewaves/eigensystem.jl") -include("effective_wave/planewaves/reflection_coefficient.jl") +include("effective_wave/planewaves/material_scattering_coefficients.jl") include("effective_wave/planewaves/boundary_condition.jl") include("effective_wave/planewaves/wavemode-wienerhopf.jl") diff --git a/src/acoustics/low_frequency.jl b/src/acoustics/low_frequency.jl index 6d174b3d..077bb725 100644 --- a/src/acoustics/low_frequency.jl +++ b/src/acoustics/low_frequency.jl @@ -54,11 +54,11 @@ function reflection_coefficient(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,D end """ - planewave_coefficients(psource::PlaneSource, reflect_medium::Acoustic, plate::Plate) + planewave_amplitudes(psource::PlaneSource, reflect_medium::Acoustic, plate::Plate) -Calculates the coefficients, or amplityudes, of the plane waves scattering from and transmitted inside the plate, for the plane wave source `psource` and the `plate`. The function returns `[R, T, P1, P2]` where `R` is the reflection coefficient, `T` is the coefficient of the transmitted wave, and `P1` (`P2`) are the amplitudes of the wave travling forward (backward) inside the plate. +Calculates the coefficients, or amplitudes, of the plane waves scattering from and transmitted inside the plate, for the plane wave source `psource` and the `plate`. The function returns `[R, T, P1, P2]` where `R` is the reflection coefficient, `T` is the coefficient of the transmitted wave, and `P1` (`P2`) are the amplitudes of the wave travling forward (backward) inside the plate. """ -function planewave_coefficients(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,Dim}}, plate_medium::Acoustic{T,Dim}, plate::Plate{T,Dim}) where {T<:AbstractFloat,Dim} +function planewave_amplitudes(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,Dim}}, plate_medium::Acoustic{T,Dim}, plate::Plate{T,Dim}) where {T<:AbstractFloat,Dim} k0 = ω / psource.medium.c k1 = if abs(plate_medium.c) == zero(T) @@ -85,17 +85,12 @@ function planewave_coefficients(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,D Z1 = Z0 - plate.width / T(2) Z2 = Z0 + plate.width / T(2) - UR = (C0 - C1)*(C0 + C1) * exp(2im*k0*Z1*cos(θ0)) * (exp(2im*k1*Z1*cos(θ1)) - exp(2im*k1*Z2*cos(θ1))) / - ((C0 + C1)^2 * exp(2im*k1*Z1*cos(θ1)) - (C0 - C1)^2 * exp(2im*k1*Z2*cos(θ1))) + denom = ((C0 + C1)^2 * exp(2im*k1*Z1*cos(θ1)) - (C0 - C1)^2 * exp(2im*k1*Z2*cos(θ1))); - UP1 = 2C0*(C0 + C1) * exp(im*Z1*(k0*cos(θ0) + k1*cos(θ1))) / - ((C0 + C1)^2*exp(2im*k1*Z1*cos(θ1)) - (C0 - C1)^2*exp(2im*k1*Z2*cos(θ1))) - - UP2 = 2*C0*(C0 - C1)*exp(im*(k0*Z1*cos(θ0) + k1*(Z1 + 2*Z2)*cos(θ1))) / - (-((C0 + C1)^2*exp((2*im)*k1*Z1*cos(θ1))) + (C0 - C1)^2*exp(2im*k1*Z2*cos(θ1))) - - UT = 4*C0*C1*exp(im*(k0*(Z1 - Z2)*cos(θ0) + k1*(Z1 + Z2)*cos(θ1))) / - ((C0 + C1)^2*exp(2im*k1*Z1*cos(θ1)) - (C0 - C1)^2*exp(2im*k1*Z2*cos(θ1))) + UR = (C0 - C1)*(C0 + C1) * exp(2im*k0*Z1*cos(θ0)) * (exp(2im*k1*Z1*cos(θ1)) - exp(2im*k1*Z2*cos(θ1))) / denom + UP1 = 2C0*(C0 + C1) * exp(im*Z1*(k0*cos(θ0) + k1*cos(θ1))) / denom + UP2 = - 2*C0*(C0 - C1)*exp(im*(k0*Z1*cos(θ0) + k1*(Z1 + 2*Z2)*cos(θ1))) / denom + UT = 4*C0*C1*exp(im*(k0*(Z1 - Z2)*cos(θ0) + k1*(Z1 + Z2)*cos(θ1))) / denom return [UR, UT, UP1, UP2] end diff --git a/src/effective_wave/wavemode.jl b/src/effective_wave/wavemode.jl index 62a2c8ef..248e9126 100644 --- a/src/effective_wave/wavemode.jl +++ b/src/effective_wave/wavemode.jl @@ -52,7 +52,7 @@ function WaveMode(ω::T, wavenumber::Complex{T}, psource::PlaneSource{T,Dim,1}, direction2 = transmission_direction(- wavenumber, (ω / psource.medium.c) * psource.direction, material.shape.normal) eigvectors2 = eigenvectors(ω, - wavenumber, psource, material; direction_eff = direction2, kws...) - α = solve_boundary_condition(ω, k_eff, eigvectors1, eigvectors2, psource, material; kws...) + α = solve_boundary_condition(ω, wavenumber, eigvectors1, eigvectors2, psource, material; kws...) # apply normalisation eigvectors1[:,:,1] = eigvectors1[:,:,1] .* α[1]