Commit 2fac86d8 authored by ollyl's avatar ollyl

Updated refraction_correction

parent 7f3e6a94
function shear!(fvec,xi,xm,Ma,h)
a = sqrt(xi[1]^2+(1-Ma^2)*(xi[2]^2+xi[3]^2))
b = norm(xm-xi)
fvec[1] = (1/a)*xi[1]-(1/b)*(1-Ma^2)*(xm[1]-xi[1]) - Ma
fvec[2] = (1/a)*xi[2]-(1/b)*(xm[2]-xi[2])
fvec[3] = h - xi[3]
ξ(Ma,a,b) = sqrt((1-Ma*a*b)^2-a^2)
function f!(F,x,Δxm,Δz2,Δz3,Ma,Δy)
if (-0.8<=x[1]<=0.8) #&& (0.2<=abs(x[2])<=1)
Δx20 = Δz2 * x[1]*x[2]/ξ(Ma,x[1],x[2])
Δx21 = Δz2 * Ma*(1-Ma*x[1]*x[2])/ξ(Ma,x[1],x[2])
Δx3 = Δz3 * x[1]*x[2]/sqrt(1-x[1]^2)
F[1] = Δx20 + Δx21 + Δx3 - Δxm
F[2] = x[2] - (Δxm - Δx21)/sqrt((Δxm-Δx21)^2+Δy^2)
return nothing
else
return nothing
end
end
"""
shearlayercorrection(E::Environment)
r,pc_pm = refraction_correction(E::Environment[,h1=1.5,h2=z0-h1])
Compute shear layer correction using Amiet's derivation for a zero-thickness shear layer. Returns time-delay and amplitude correction. Relies on `shear=true`, `Ma` and `h` set in `Environment()`. The coordinate system is assumed to be at the microphone array center and `h` is the distance from the array center to the shear layer (this is contrary to the references).
Compute propagation time correction and amplitude correction for 2D planar jet.
References:
- R. K. Amiet, “Refraction of sound by a shear layer,” Journal of Sound and Vibration, vol. 58, no. 4, pp. 467–482, 1978.
- C. Bahr, et. al “Shear Layer Correction Validation Using A Non-Intrusive Acoustic Point Source,” presented at the 16th AIAA/CEAS Aeroacoustics Conference, Reston, Virigina, 2012.
- Glegg, S., & Devenport, W. (2017). Aeroacoustics of low mach number flows: Fundamentals, analysis, and measurement. Chap. 10.
"""
function shearlayercorrection(E::Environment)
@unpack N,M,micgeom_s,Rxy,Ma,h,c,z0,ampcorr = E
ta = Array{Float64,2}(undef,N,M)
pcpm2 = Array{Float64,2}(undef,N,M)
xn = zeros(3,M)
for n = 1:N
xn .= Rxy[:,n] .- micgeom_s
for m = 1:M
res = nlsolve((fvec,x)->shear!(fvec,x,xn[:,m],-Ma,h),[.1,.1,h])
xi = res.zero
r1 = norm(xi)
d = -xi[1]*Ma*c/r1
c1 = d + sqrt(d^2+c^2-(Ma*c)^2)
r2v = xn[:,m]-xi
r2 = norm(r2v)
ta[n,m] = r1/c + r2/c1
if ampcorr
theta = acos(-sign(Ma)*xi[1])
hH = (z0-h)/z0
zeta = sqrt((1 - abs(Ma)*cos(theta))^2 - cos(theta)^2)
pcpm2[n,m]=1/4/zeta^2*hH^2*(1+(1/hH-1)*zeta.^3/sin(theta)^3)*(1+(1/hH-1)*zeta/sin(theta))*(zeta+sin(theta)*(1-abs(Ma)*cos(theta))^2)^2
end
function refraction_correction(E,h1=1.5,h2=E.z0-h1)
r = zeros(E.N,E.M)
pc_pm = similar(r)
for n = 1:E.N
for m = 1:E.M
Δx,Δy,Δz = E.Rxy[:,n]-E.micgeom[:,m] # vector from mic to grid point
r0 = sqrt(Δx^2+Δy^2+Δz^2)
a0 = sqrt(Δx^2+Δy^2)/r0
b0 = Δx/sqrt(Δx^2+Δy^2)
Δz2,Δz3 = h1,h2
res = nlsolve((F,x)->f!(F,x,Δx,Δz2,Δz3,E.Ma,Δy),[a0,b0])
a,b = res.zero
Δx20 = Δz2*a*b/ξ(E.Ma,a,b)
Δx3 = Δz3*a*b/sqrt(1-a^2)
r[n,m] = Δz3/sqrt(1-a^2) + Δz2*(1-E.Ma*a*b)/ξ(E.Ma,a,b)
# Amplitude correction:
pc_pi = Δz2^2/(Δz2+Δz3)^2
pi_pt = (sqrt(1-a^2)*(1-E.Ma*a*b)^2+ξ(E.Ma,a,b))^2/(4*ξ(E.Ma,a,b)^2)
K(R,A,D) = sqrt(a^2*(D*sqrt(1-a^2)*(1-b^2)/ξ(E.Ma,a,b) + A*(b-D/R) )^2 + a^2*(1-b^2)*(A-D*sqrt(1-a^2)*b/ξ(E.Ma,a,b))^2 + (1-a^2)*(R-D*b)^2)
d = 1-E.Ma*a*b
r2 = (1-E.Ma*a*b)*Δz2/ξ(E.Ma,a,b)
r3 = Δz3/sqrt(1-a^2)
ρ2 = r2/d^2
R2 = ρ2
R3 = r3+ρ2
A2 = sqrt(1-a^2)*ρ2/ξ(E.Ma,a,b)
A3 = r3+sqrt(1-a^2)*ρ2/ξ(E.Ma,a,b)
D2 = ρ2*E.Ma*a
D3 = ρ2*E.Ma*a
pt_pm = R3*K(R3,A3,D3)/(R2*K(R2,A2,D2))
pc_pm[n,m] = pc_pi * pi_pt * pt_pm
end
end
return ta,pcpm2
return r,pc_pm
end
"""
......@@ -64,20 +84,19 @@ function steeringvectors!(E::Environment)
if E.shear
w = 2pi*fn
ta,pcpm2 = AeroAcoustics.shearlayercorrection(E)
dx,pc_pm = refraction_correction(E,h)
ta = dx./c
for j in 1:Nf
vi[:,:,j] .= 1 ./(ta.*c.*D0.*sum(1 ./(ta.*c).^2,dims=2)).*exp.(-im.*w[j].*ta)
#vi[:,:,j] .= (1 ./M).*(ta.*c./D0).*exp.(-im.*w[j].*ta)
end
if E.ampcorr
for j in 1:Nf
vi[:,:,j] .*= sqrt.(pcpm2)
vi[:,:,j] .*= pc_pm
end
end
else
kw = 2pi*fn/c
for j in 1:Nf
#vi[:,:,j] .= (1 ./M).*(D./D0).*exp.(-im.*kw[j].*(D.-D0))
vi[:,:,j] .= 1 ./(D.*D0.*sum(1 ./D.^2,dims=2)).*exp.(-im.*kw[j].*(D.-D0))
end
end
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment