Analytical solution to acoustic scattering by a bubble

Sorin Mitran (mitran@unc.edu, 2/29/24)

Clear["Global`*"]

Plane waves

p[x_,t_,k_,n_,c_,P_] = P Exp[I k (x.n - c t)]

Examples

Plane  wave  in the x-direction in 2D

Manipulate[ContourPlot[Re[p[{x,y},t,1,{1,1}/Sqrt[2],1,1]],{x,-3π,3π},{y,-3π,3π},PlotPoints->50,Contours->Table[v,{v,-1,1,0.2}]],{t,0,2π}]

Ultrasound excitation of intravenously injected microbubbles

Import["/Users/mitran/courses/MATH529L/MBvideo.avi"]

Spherical vector system

(r,θ,ϕ)
0⩽θ⩽π, -π<ϕ⩽π
x = r sin(θ) cos(ϕ), y = r sin(θ) sin(ϕ), z = r cos(θ)

FromSphericalCoordinates[{r,θ,ϕ}]

Y[l_,m_,θ_,ϕ_]=Exp[I m ϕ] LegendreP[l,m,Cos[θ]]

grad = (OverHat[r]) (∂_r)+(OverHat[(θ/r)]) (∂_θ)+(OverHat[(ϕ/(r sin(θ)))]) (∂_ϕ)

gY[l_,m_,r_,θ_,ϕ_]=Grad[Y[l,m,θ,ϕ],{r,θ,ϕ},"Spherical"]

Y0 = r x grad(Ylm), (OverHat[r])×(OverHat[θ])=(OverHat[ϕ]), (OverHat[r])×(OverHat[ϕ])=- (OverHat[θ])

Y0[l_,m_,θ_,ϕ_]={0, -r Part[gY[l,m,r,θ,ϕ],3], -r Part[gY[l,m,r,θ,ϕ],3]}

Y0[l_,θ_]=Y0[l,0,θ,0]

Yp[l_,m_,θ_,ϕ_]={(l+1) Y[l,m,θ,ϕ],-r Part[gY[l,m,r,θ,ϕ],2],-r Part[gY[l,m,r,θ,ϕ],3]}

Yp[l_,θ_]=Yp[l,0,θ,0]

Ym[l_,m_,θ_,ϕ_]={l Y[l,m,θ,ϕ],r Part[gY[l,m,r,θ,ϕ],2],r Part[gY[l,m,r,θ,ϕ],3]}

Ym[l_,θ_]=Ym[l,0,θ,0]

Canonical waves

P[l_,m_,r_,θ_,ϕ_,k_]=SphericalBesselJ[l+1,k r ] Yp[l,m,θ,ϕ]-SphericalBesselJ[l-1,k r] Ym[l,m,θ,ϕ]

P[l_,r_,θ_,k_]=SphericalBesselJ[l+1,k r ] Yp[l,θ]-SphericalBesselJ[l-1,k r] Ym[l,θ]

uxz[urθϕ_]:={urθϕ[[1]] Sin[urθϕ[[2]]],urθϕ[[1]] Cos[urθϕ[[2]]]}

uxz[P[1,1.,π/4.,1.]]

ContourPlot[ If[(x^2)+(z^2)≤1.,uxz[P[1,Sqrt[(x^2)+(z^2)],ArcTan[x/z],1]][[1]],0.],{x,0.001,1.},{z,-1.,1.},AspectRatio->Automatic,PlotPoints->200]

ContourPlot[ If[(x^2)+(z^2)≤1.,uxz[P[2,Sqrt[(x^2)+(z^2)],ArcTan[x/z],1]][[1]],0.],{x,0.001,1.},{z,-1.,1.},AspectRatio->Automatic,PlotPoints->200]

ContourPlot[ If[(x^2)+(z^2)≤1.,uxz[P[4,Sqrt[(x^2)+(z^2)],ArcTan[x/z],1]][[1]],0.],{x,0.001,1.},{z,-1.,1.},AspectRatio->Automatic,PlotPoints->200]

ux=Table[{r Sin[θ],r Cos[θ],uxz[P[6,r,θ,1.]][[1]]},{r,0.02,1.,0.02},{θ,0.001,π/2-0.001,π/36.}];

Limit[P[l,r,θ,k],θ->0]

SV[l_,m_,r_,θ_,ϕ_,k_]=l SphericalBesselJ[l+1,k r] Yp[l,m,θ,ϕ]+(l+1)SphericalBesselJ[l-1,k r] Ym[l,m,θ,ϕ]

SV[l_,r_,θ_,k_]=l SphericalBesselJ[l+1,k r] Yp[l,θ]+(l+1)SphericalBesselJ[l-1,k r] Ym[l,θ]

SH[l_,m_,r_,θ_,ϕ_,k_]= SphericalBesselJ[l+1,k r] Y0[l,m,θ,ϕ]

SH[l_,r_,θ_,k_]= SphericalBesselJ[l+1,k r] Y0[l,θ]

Boundary conditions

Define the (l^th) component of incident, scattered, refracted

Incident plane wave

U0[l_,r_,θ_,ωα2_]=j[l+1,ωα2 r] Yp[l,θ] - j[l-1,ωα2 r] Ym[l,θ];
ψ0p[l_,r_,ωα2_]=j[l+1,ωα2 r]; ψ0m[l_,r_,ωα2_]=- j[l-1,ωα2 r];

Refracted wave

U1[l_,r_,θ_,ωα1_,ωβ1_]=(al1 j[l+1,ωα1 r] + l bl1 j[l+1,ωβ1 r] ) Yp[l,θ] - (- al1 j[l-1,ωα1  r] + (l+1)  bl1 j[l-1,ωβ1  r] )Ym[l,θ];
ψ1p[l_,r_,ωα1_,ωβ1_]=al1 j[l+1,ωα1 r] + l bl1 j[l+1,ωβ1 r] ; ψ1m[l_,r_,ωα1_,ωβ1_]=- al1 j[l-1,ωα1  r] + (l+1)  bl1 j[l-1,ωβ1  r];

Scattered wave

U2[l_,r_,θ_,ωα2_]=al2  h[l+1,ωα2 r] Yp[l,θ] - al2   h[l-1,ωα2 r] Ym[l,θ];
ψ2p[l_,r_,ωα2_]=al2  h[l+1,ωα2 r]; ψ2m[l_,r_,ωα2_]=- al2   h[l-1,ωα2 r];

Continuity of radial displacement

Define the projection of a displacement component along radial direction

Pr[l_,ψp_,ψm_]=(l+1) ψp + l ψm

eq1 = Pr[l,ψ0p[l,r,ωα2]+ψ2p[l,r,ωα2],ψ0m[l,r,ωα2]+ψ2m[l,r,ωα2]] -
         Pr[l,ψ1p[l,r,ωα1,ωβ1],ψ1m[l,r,ωα1,ωβ1]]

eq1S = Simplify[eq1 /. {j->SphericalBesselJ, h->SphericalHankelH2}]

Continuity of radial stress

Define the divergence of a displacement component

div[l_,ψp_,ψm_]:= (l+1)(D[ψp,r]+((l+2)/r) ψp)+l (D[ψm,r]-((l-1)/r) ψp)

Define the radial derivative of the radial displacement component

dr[l_,ψp_,ψm_]:= (l+1) D[ψp,r] +l D[ψm,r]

eq2 =Simplify[ λ2 div[l,ψ0p[l,r,ωα2]+ψ2p[l,r,ωα2],ψ0m[l,r,ωα2]+ψ2m[l,r,ωα2] ] - λ1 div[l,ψ1p[l,r,ωα1,ωβ1],ψ1m[l,r,ωα1,ωβ1]] + 2 μ1 dr[l,ψ1p[l,r,ωα1,ωβ1],ψ1m[l,r,ωα1,ωβ1]]]

eq2S = Simplify[eq2 /. {j->SphericalBesselJ, h->SphericalHankelH2}]

Zero shear stress in sphere to match μ1=0 in fluid

The shear stress is defined as (Subscript[σ, rθ])=μ((1/r) ((∂ (Subscript[u, r]))/(∂ θ))+ ((∂ (Subscript[u, θ]))/(∂ r))-((Subscript[u, θ])/( r)))

σrθ[u_]:=Simplify[(1/r)D[u[[1]],θ]+D[u[[2]],r]-(1/r) u[[2]]]

A component along Ym[l,θ] gives contribution

θrθm=σrθ[jm[r] Ym[l,θ]]

A component along Yp[l,θ] gives contribution

θrθp=σrθ[jp[r] Yp[l,θ]]

These are colinear in θ, with ratio

θrθp/θrθm

Based on above, define effect of shear stress operator on coefficients in U1 expansion

ss[l_,ψp_,ψm_]:=(l+2) ψp -r D[ψp,r] + (l-1) ψm + r D[ψm,r]

eq3 = Simplify[ss[l,ψ1p[l,r,ωα1,ωβ1],ψ1m[l,r,ωα1,ωβ1]]]

eq3S = Simplify[eq3 /. j->SphericalBesselJ]

Form linear system

A=IdentityMatrix[3]; b={0,0,0}; AS=IdentityMatrix[3]; bS={0,0,0};

A[[1,1]]=Simplify[Coefficient[eq1,al1]]; A[[1,1]] /. l->L
A[[1,2]]=Simplify[Coefficient[eq1,bl1]]; A[[1,2]] /. l->L
A[[1,3]]=Simplify[Coefficient[eq1,al2]]; A[[1,3]]/. l->L
b[[1]]=Simplify[-(eq1 - al1 A[[1,1]] - bl1 A[[1,2]] - al2 A[[1,3]])]; b[[1]] /. l->L

AS[[1,1]]=Simplify[Coefficient[eq1S,al1]]; AS[[1,1]] /. {SphericalBesselJ->j, SphericalHankelH2->h}
AS[[1,2]]=Simplify[Coefficient[eq1S,bl1]]; AS[[1,2]] /. {SphericalBesselJ->j, SphericalHankelH2->h}
AS[[1,3]]=Simplify[Coefficient[eq1S,al2]]; AS[[1,2]] /. {SphericalBesselJ->j, SphericalHankelH2->h}
bS[[1]]=Simplify[-(eq1S - al1 AS[[1,1]] - bl1 AS[[1,2]] - al2 AS[[1,3]])]; bS[[1]] /. {SphericalBesselJ->j, SphericalHankelH2->h}

A[[2,1]]=Simplify[Coefficient[eq2,al1]]; A[[2,1]] /. l->L
A[[2,2]]=Simplify[Coefficient[eq2,bl1]]; A[[2,2]] /. l->L
A[[2,3]]=Simplify[Coefficient[eq2,al2]]; A[[2,3]] /. l->L
b[[2]]=Simplify[-(eq2 - al1 A[[2,1]] - bl1 A[[2,2]] - al2 A[[2,3]])]; b[[2]] /. l->L

AS[[2,1]]=Simplify[Coefficient[eq2S,al1]]; AS[[2,1]] /. {SphericalBesselJ->j, SphericalHankelH2->h}
AS[[2,2]]=Simplify[Coefficient[eq2S,bl1]]; AS[[2,2]] /. {SphericalBesselJ->j, SphericalHankelH2->h}
AS[[2,3]]=Simplify[Coefficient[eq2S,al2]]; AS[[2,2]] /. {SphericalBesselJ->j, SphericalHankelH2->h}
bS[[2]]=Simplify[-(eq2S - al1 AS[[2,1]] - bl1 AS[[2,2]] - al2 AS[[2,3]])]; bS[[2]] /. {SphericalBesselJ->j, SphericalHankelH2->h}

A[[3,1]]=Simplify[Coefficient[eq3,al1]]; A[[3,1]] /. l->L
A[[3,2]]=Simplify[Coefficient[eq3,bl1]]; A[[3,2]] /. l->L
A[[3,3]]=Simplify[Coefficient[eq3,al2]]; A[[3,3]] /. l->L
b[[3]]=Simplify[-(eq3 - al1 A[[3,1]] - bl1 A[[3,2]] - al2 A[[3,3]])]; b[[3]] /. l->L

AS[[3,1]]=Simplify[Coefficient[eq3S,al1]]; AS[[3,1]] /. {SphericalBesselJ->j, SphericalHankelH2->h}
AS[[3,2]]=Simplify[Coefficient[eq3S,bl1]]; AS[[3,2]] /. {SphericalBesselJ->j, SphericalHankelH2->h}
AS[[3,3]]=Simplify[Coefficient[eq3S,al2]]; AS[[3,3]] /. {SphericalBesselJ->j, SphericalHankelH2->h}
bS[[3]]=Simplify[-(eq3S - al1 AS[[3,1]] - bl1 AS[[3,2]] - al2 AS[[3,3]])]; bS[[3]] /. {SphericalBesselJ->j, SphericalHankelH2->h}

Solve linear system

The analytical solution can be found, but is unwieldy

sol=Simplify[Solve[A.{al1,bl1,al2}==b,{al1,bl1,al2}] ];

Plane P-wave in ϕ=0 plane (m=0)

U0[r_,θ_,k_,n_]:=Sum[Re[N[P[l,r,θ,k] Exp[-I (π/2)(l+1)]]],{l,0,n}]

Table[U0[2/Sin[θ],θ,1,50],{θ,π/18,17π/18,π/18}] // TableForm

Uzx[Usp_,θ_]:={Usp[[1]] Sin[θ] + Usp[[2]] Cos[θ],Usp[[1]]  Cos[θ] -  Usp[[2]] Sin[θ]}

Table[Uzx[U0[2/Sin[θ],θ,1,50],θ],{θ,π/18,17π/18,π/18}] // TableForm

zUz[Usp_,r_,θ_]:={N[r Cos[θ]],Part[Uzx[Usp,θ],2]}

ListPlot[Table[Table[zUz[U0[z,0.000001,k,50],z,0.000001],{z,0.01,6.28,0.01}],{k,1,4}] ,Joined->True,PlotLegends->Automatic]

δθ=π/1440; x=0.5; ph=ListPlot[Table[zUz[U0[x/Sin[θ],θ,4,50],x/Sin[θ],θ],{θ,ArcTan[x/(2π)],π/2,δθ}],Joined->True,PlotRange->All,PlotLegends->Automatic,PlotMarkers->Automatic]

δθ=π/1440; x=1.; p1=ListPlot[Table[zUz[U0[x/Sin[θ],θ,4,50],x/Sin[θ],θ],{θ,ArcTan[x/(2π)],π/2,δθ}],Joined->True,PlotRange->All,PlotLegends->Automatic,PlotMarkers->Automatic]

δθ=π/1440; x=2.; p2=ListPlot[Table[zUz[U0[x/Sin[θ],θ,4,50],x/Sin[θ],θ],{θ,ArcTan[x/(2π)],π/2,δθ}],Joined->True,PlotRange->All,PlotLegends->Automatic,PlotMarkers->Automatic]

δθ=π/1440; x=4.; p4=ListPlot[Table[zUz[U0[x/Sin[θ],θ,4,50],x/Sin[θ],θ],{θ,ArcTan[x/(2π)],π/2,δθ}],Joined->True,PlotRange->All,PlotLegends->Automatic,PlotMarkers->Automatic]

Verify plane wave behavior

Show[{ph,p1,p2,p4}]

Investigate number of summation terms required to obtain plane wave

n=50

δθ=π/1440; x=2.; ListPlot[Table[zUz[U0[x/Sin[θ],θ,10,50],x/Sin[θ],θ],{θ,ArcTan[x/(2π)],π/2,δθ}],Joined->True,PlotRange->All,PlotLegends->Automatic,PlotMarkers->Automatic]

n=100

δθ=π/1440; x=2.; ListPlot[Table[zUz[U0[x/Sin[θ],θ,10,100],x/Sin[θ],θ],{θ,ArcTan[x/(2π)],π/2,δθ}],Joined->True,PlotRange->All,PlotLegends->Automatic,PlotMarkers->Automatic]

Roughly 10 times wavenumber needed to satisfy Nyquist criterion up to z=2π.

Displacement field inside the sphere (in ϕ=0 plane, i.e. axisymmetric)

ASn[l_,r_,ωα1_,ωα2_,ωβ1_,λ1_,λ2_,μ1_]=AS;

bSn[l_,r_,ωα1_,ωα2_,ωβ1_,λ1_,λ2_,μ1_]=bS;

coef[l_,r_,T_,ρ1_,Vp1_,Vs1_,ρ2_,Vp2_]:= Module[{ω,α1,α2,β1,ωα1,ωα2,ωβ1,λ1,λ2,μ1,sol,al1,bl1,al2,AS,bS},
ω=2.π/T;α1=1./Vp1; α2=1./Vp2; β1=1./Vs1; ωα1=ω α1; ωα2 = ω α2; ωβ1=ω β1;
μ1= ρ1 (Vs1^2); λ1=ρ1 (Vp1^2)-2μ1; λ2=ρ2 (Vp2^2);
AS=ASn[l,r,ωα1,ωα2,ωβ1,λ1,λ2,μ1];
bS=bSn[l,r,ωα1,ωα2,ωβ1,λ1,λ2,μ1];
sol= Solve[ AS.{al1,bl1,al2} == bS,{al1,bl1,al2}];
Return[ {al1,bl1,al2} /. sol] ]

bval[l_,r_,T_,ρ1_,Vp1_,Vs1_,ρ2_,Vp2_]:= Module[{ω,α1,α2,β1,ωα1,ωα2,ωβ1,λ1,λ2,μ1,sol,al1,bl1,al2,AS,bS},
ω=2.π/T;α1=1./Vp1; α2=1./Vp2; β1=1./Vs1; ωα1=ω α1; ωα2 = ω α2; ωβ1=ω β1;
μ1= ρ1 (Vs1^2); λ1=ρ1 (Vp1^2)-2μ1; λ2=ρ2 (Vp2^2);
bS=bSn[l,r,ωα1,ωα2,ωβ1,λ1,λ2,μ1];
Return[bS]];

Aval[l_,r_,T_,ρ1_,Vp1_,Vs1_,ρ2_,Vp2_]:= Module[{ω,α1,α2,β1,ωα1,ωα2,ωβ1,λ1,λ2,μ1,sol,al1,bl1,al2,AS,bS},
ω=2.π/T;α1=1./Vp1; α2=1./Vp2; β1=1./Vs1; ωα1=ω α1; ωα2 = ω α2; ωβ1=ω β1;
μ1= ρ1 (Vs1^2); λ1=ρ1 (Vp1^2)-2μ1; λ2=ρ2 (Vp2^2);
AS=ASn[l,r,ωα1,ωα2,ωβ1,λ1,λ2,μ1];
Return[ AS] ]

A1=Aval[4,0.01,(10.^(-3)),1700.,3000.,1500.,1000.,1500.] ; A1// MatrixForm

bval[4,0.01,(10.^(-3)),1700.,3000.,1500.,1000.,1500.] // MatrixForm

coef[4,0.01,(10.^(-3)),1700.,3000.,1500.,1000.,1500.]

BVPs for the Poisson equation in spherical coordinates

Dirichlet BVP problem(∇^2)u= f(r,θ,ϕ) for 0< r<1, u(r=1,θ,ϕ)=g(θ,ϕ)

Simple analytical solutions

Homogeneous, constant BC

f(r,θ,ϕ)=0, g(θ,ϕ)=c
Solution: u=c

Inhomogeneous, constant BC

f(r,θ,ϕ)=r, g(θ,ϕ)=c
Solution: u=u(r)

(∇^2)u=(((∂^2)u)/(∂(r^2)))+(2/r)((∂u)/(∂r))=r,u(r=1)=c,u(0) finite

Eq=u''[r]+2 u'[r]/r == r;
sol=DSolve[{Eq,u[1]==c},u[r],r][[1,1]]

usol[r_]=Expand[u[r] /. sol /. (C[1])->0]

Above can be expanded on eigenbasis of Sturm-Liouville problem along r.

Homogeneous, variable BC on surface of sphere

(∇^2)u= f(r,θ,ϕ) for 0< r<1, u(r=1,θ,ϕ)=g(θ,ϕ) 
f(r,θ,ϕ)=0, g(θ,ϕ)=sin(cos(θ))
Solution: u=u(r,θ)

(∇^2)u=(((∂^2)u)/(∂(r^2)))+(2/r)((∂u)/(∂r))+(1/(r^2))(((∂^2)u)/(∂(θ^2)))+((cot θ)/(r^2))((∂u)/(∂θ))=0,u(r=1)=sin(cos(θ)),u(0) finite

u(r,θ)=(∑_(l=0)^∞)(Subscript[A, l]) (r^l)(Subscript[P, l])(cos θ)

u(1,θ)=sin(cos(θ)) = (∑_(l=0)^∞)(Subscript[A, l]) (Subscript[P, l])(cos θ)

(b,c)=(∫_0^π)b(θ)c(θ)sin θdθ

Problems

Problem 1: Homogeneous, variable azimuthal BC on surface of sphere

Solve BVP below, plot solution through contour lines, plot gradient of solution, plot norm of gradient, evaluate norm of gradient over entire sphere.

(∇^2)u= f(r,θ,ϕ) for 0< r<1, u(r=1,θ,ϕ)=g(θ,ϕ) 
f(r,θ,ϕ)=0, g(θ,ϕ)=A sin(cos(k θ)),k=13
Solution: u=u(r,θ)

(∇^2)u=(((∂^2)u)/(∂(r^2)))+(2/r)((∂u)/(∂r))+(1/(r^2))(((∂^2)u)/(∂(θ^2)))+((cot θ)/(r^2))((∂u)/(∂θ))=0,u(r=1)=sin(cos(θ)),u(0) finite

u(r,θ)=(∑_(l=0)^∞)(Subscript[A, l]) (r^l)(Subscript[P, l])(cos θ)

u(1,θ)=sin(cos(θ)) = (∑_(l=0)^∞)(Subscript[A, l]) (Subscript[P, l])(cos θ)

(b,c)=(∫_0^π)b(θ)c(θ)sin θdθ

g[θ_]=Sin[Cos[3θ]];
PolarPlot[{1,1+4 (10^(-2)) g[θ]},{θ,0,2π}]

Problem 2: Homogeneous, general variable BC on surface of sphere

Solve BVP below, plot solution through contour lines, plot gradient of solution, plot norm of gradient, evaluate norm of gradient over entire sphere.

(∇^2)u= f(r,θ,ϕ) for 0< r<1, u(r=1,θ,ϕ)=g(θ,ϕ) 
f(r,θ,ϕ)=0, g(θ,ϕ)=A sin(cos(k θ)) cos(sin(j ϕ));
Solution: u=u(r,θ,ϕ)

(∇^2)u=(((∂^2)u)/(∂(r^2)))+(2/r)((∂u)/(∂r))+(1/(r^2))(((∂^2)u)/(∂(θ^2)))+((cot θ)/(r^2))((∂u)/(∂θ))=0,u(r=1)=sin(cos(θ)),u(0) finite

u(r,θ)=(∑_(l=0)^∞)(Subscript[A, l]) (r^l)(Subscript[P, l])(cos θ)

u(1,θ)=sin(cos(θ)) = (∑_(l=0)^∞)(Subscript[A, l]) (Subscript[P, l])(cos θ)

(b,c)=(∫_0^π)b(θ)c(θ)sin θdθ

g[θ_]=Sin[Cos[13θ]];
PolarPlot[{1,1+4 (10^(-2)) g[θ]},{θ,0,2π}]