Scattering by spheres

References
[1] J. J. Faran, J. Ac. Soc. Am., 23(4):405-418, 1951

SetDirectory["/home/mitran/papers/2017/LithotripsyValidation/notebooks"];

Analytical evaluation

Definitions

Elastic displacement in sphere

Represent displacement through scalar, vector potentials (Overscript[u, →])=−∇Ψ+∇×(Overscript[A, →]), and introduce Bessel series expansion of the potentials

j[n_,x_]=SphericalBesselJ[n,x]; y[n_,x_]=SphericalBesselY[n,x];
jp[n_,x_] = Simplify[D[SphericalBesselJ[n,x],x]]; jpp[n_,x_]=Simplify[D[SphericalBesselJ[n,x],{x,2}]];
yp[n_,x_] = Simplify[D[SphericalBesselY[n,x],x]];
P[n_,θ_]=LegendreP[n,Cos[θ]]; Q[n_,θ_]=LegendreQ[n,Cos[θ]];
Pp[n_,θ_]=Simplify[D[LegendreP[n,Cos[θ]],θ]];
Ppp[n_,θ_]=Simplify[D[LegendreP[n,Cos[θ]],{θ,2}]];

Plot[Evaluate[Table[P[n,θ],{n,0,6}]],{θ,0,2π},PlotLegends->Table[n,{n,0,6}]]

Limiting behavior of spherical Legendre polynomial derivatives

Limit[Pp[n,θ],θ->0]

Limit[Pp[n,θ],θ->π]

Ppp[n,θ]

Limit[Ppp[n,θ],θ->0]

Introduce regularized forms

Ppr[n_,θ_]:= If[Abs[θ-Round[θ/π]π]<(10^(-8)),0,Pp[n,θ]];
Pppr[n_,θ_]:= If[Abs[θ-Round[θ/π]π]<(10^(-8)),-(1/2) n (1+n) (((-1))^(Round[θ/π])),Ppp[n,θ]];

Plot[Evaluate[Table[Q[n,θ],{n,0,6}]],{θ,0,2π}]

Plot[Evaluate[Table[Pp[n,θ],{n,0,6}]],{θ,0,2π}]

Plot[Evaluate[Table[Ppr[n,θ],{n,0,6}]],{θ,0,2π}]

Q[0,0]

Plot[Evaluate[Table[Ppp[n,θ],{n,0,6}]],{θ,0,2π}]

Plot[Evaluate[Table[Pppr[n,θ],{n,0,6}]],{θ,0,2π}]

Ψ[r_,θ_,k1_,L_]:=Sum[a[n] j[n,k1 r] P[n,θ],{n,0,L}];
Aϕ[r_,θ_,k2_,L_]:=Sum[b[n] j[n,k2 r] Ppr[n,θ],{n,0,L}];

Verify above as solutions of the appropriate Helmholtz equation in axisymmetric spherical coordinates

LaplaceOperator[f_]:=Simplify[(1/(r^2))D[(r^(2 )) D[f,r],r]+(1/((r^2)Sin[θ]))D[Sin[θ] D[f,θ],θ]];

FullSimplify[LaplaceOperator[j[n,k1 r] P[n,θ]]+(k1^2) j[n,k1 r] P[n,θ],Assumptions->{θ∈Reals,k1>0,n∈Integers}]

FullSimplify[LaplaceOperator[j[n,k2 r] Q[n,θ]]+(k2^2) j[n,k2 r] Q[n,θ],Assumptions->{θ∈Reals,k2>0,n∈Integers}]

Displacement components

ur[r_,θ_,k1_,k2_,L_]:=-Sum[(a[n]k1 jp[n,k1 r]+((n(n+1) b[n])/r)j[n,k2 r] ) P[n, θ],{n,0,L}];
uθ[r_,θ_,k1_,k2_,L_]:=-Sum[(((a[n])/r)j[n,k1 r]+b[n] k2 jp[n,k2 r]+((b[n])/r) j[n,k2 r] ) Ppr[n, θ],{n,0,L}];

uz[r_,θ_,k1_,k2_,L_]:=ur[r,θ,k1,k2,L]Cos[θ]-uθ[r,θ,k1,k2,L] Sin[θ];

Displacement derivatives

durdr[r_,θ_,k1_,k2_,L_]:=Sum[(-a[n](k1^2)jpp[n,k1 r]-b[n] ((n(n+1))/r)k2 jp[n,k2 r]+ b[n] ((n(n+1))/(r^2))j[n,k2 r]) P[n ,θ],{n,0,L}];
durdθ[r_,θ_,k1_,k2_,L_]:=-Sum[(a[n]k1 jp[n,k1 r]+b[n](( n(n+1))/r)j[n,k2 r] )Ppr[n, θ],{n,0,L}];
duθdr[r_,θ_,k1_,k2_,L_]:=Sum[(((a[n])/(r^2))j[n,k1 r]-((a[n])/r)k1 jp[n,k1 r]-b[n] (k2^2) jpp[n,k2 r]-((b[n])/r) k2 jp[n,k2 r]+((b[n])/(r^2))j[n,k2 r] ) Ppr[n ,θ],{n,0,L}];
duθdθ[r_,θ_,k1_,k2_,L_]:=-Sum[((( a[n])/r)j[n,k1 r]+b[n] k2 jp[n,k2 r]+((b[n])/r) j[n,k2 r] ) Pppr[n ,θ],{n,0,L}];

Displacement divergence

Δ[r_,θ_,k1_,L_]:=(k1^2) Sum[a[n]j[n,k1 r]  P[n, θ],{n,0,L}];

Stress components

σrr[r_,θ_,k1_,k2_,λ_,μ_,L_]:= λ Δ[r,θ,k1,L] + 2 μ durdr[r,θ,k1,k2,L];
σθθ[r_,θ_,k1_,k2_,λ_,μ_,L_]:= λ Δ[r,θ,k1,L] + ((2 μ)/r) (duθdθ[r,θ,k1,k2,L]+ur[r,θ,k1,k2,L]);
σϕϕ[r_,θ_,k1_,k2_,λ_,μ_,L_]:= λ Δ[r,θ,k1,L] + ((2 μ)/r)(ur[r,θ,k1,k2,L]+Cot[θ]uθ[r,θ,k1,k2,L]);
σrθ[r_,θ_,k1_,k2_,λ_,μ_,L_]:=μ (((durdθ[r,θ,k1,k2,L])/r)+duθdr[r,θ,k1,k2,L]-((uθ[r,θ,k1,k2,L])/r));

Principal stresses

σT[σrr_,σrθ_,σθθ_,σϕϕ_]:=Max[(1/2)(σrr+σθθ)+(Sqrt[((σrθ)^2)+(1/4)(((σrr-σθθ))^2)]),σϕϕ];
σC[σrr_,σrθ_,σθθ_,σϕϕ_]:=-Min[(1/2)(σrr+σθθ)-(Sqrt[((σrθ)^2)+(1/4)(((σrr-σθθ))^2)]),σϕϕ];

Pressure and displacement in fluid

Incident plane wave

pi[P0_,k_,r_,θ_]=P0 Exp[-I k r Cos[θ]];

Bessel expansion of plane wave

pi[P0_,k_,r_,θ_,L_]:= P0 Sum[ (2n+1) (((-I))^n) j[n,k r] P[n, θ],{n,0,L}];

Scattered wave

ps[r_,θ_,k3_,L_]:=Sum[c[n] (j[n,k3 r]-I y[n,k3 r]) P[n, θ],{n,0,L}];

Sphere boundary condition

M [n_,x1_,x2_,x3_,a_,ω_,ρ1_,ρ3_,c2_,ν_]= { {x1 jp[n,x1] , n(n+1) j[n,x2] , ((x3 )/((ω^2) ρ3))(jp[n,x3] - I yp[n,x3]) },
{2 ρ1 (c2^2)(x1^2) ((ν/(1-2ν)) j[n,x1] - jpp[n,x1]),-2 ρ1 (c2^2) n(n+1) ( x2 jp[n,x2]-j[n,x2]),(a^2)(j[n,x3]-I y[n,x3])},
{2(j[n,x1]-x1 jp[n,x1]),-(x2^2) jpp[n,x2]+(2-n(n+1)) j[n,x2],0}} ;

f[n_,x3_,a_,ω_,ρ3_,P0_]={-((P0 x3)/((ω^2) ρ3)) (2n+1) (((-I))^n) jp[n,x3],-P0 (2n+1) (((-I))^n) (a^2) j[n,x3],0};

f[n,x3,a,ω,ρ3,P0]

sol[n_,x1_,x2_,x3_,a_,ω_,ρ1_,ρ3_,c2_,ν_,P0_]=Simplify[LinearSolve[ M[n,x1,x2,x3,a,ω,ρ1,ρ3,c2,ν],f[n,x3,a,ω,ρ3,P0]]];

Save["SphereBCsolution.m",sol];

tanα[n_,x_]= -x jp[n,x]/j[n,x];
tanβ[n_,x_]= -x yp[n,x]/y[n,x];
tanδ[n_, x_]= -j[n,x]/y[n,x];
ζScatteringAngle[n_,x1_,ν_]:=Module[{tanα1,tanα2,tanζ1,x2},
x2=x1 Sqrt[2((1-ν)/(1-2ν))];
tanα1=If[x1==0,-n,tanα[n,x1]];
tanα2=If[x2==0,-n,tanα[n,x2]];
tanζ1=If[x1==0,If[n==1,-π/4,0],((-(x2^2))/2) ((((tanα1)/(tanα1+1))-(((n^2)+n)/(tanα2+(n^2)+n-1-(x2^2)/2)))/(((2tanα1+(n^2)+n-(x2^2)/2)/(tanα1+1))-((((n^2)+n)(tanα2+1))/(tanα2+(n^2)+n-1-(x2^2)/2))))];
Return[ArcTan[tanζ1]180/π];
];

Scattering angle, Eq (30) of Farran. Limiting behavior

ltanα[n_]=Limit[tanα[n,x],x->0]

Table[ltanα[n],{n,0,100}]

expr=Simplify[((((tanα1)/(tanα1+1))-(((n^2)+n)/(tanα2+(n^2)+n-1-(x2^2)/2)))/(((2tanα1+(n^2)+n-(x2^2)/2)/(tanα1+1))-((((n^2)+n)(tanα2+1))/(tanα2+(n^2)+n-1-(x2^2)/2)))) /. {x2->0,tanα1->(1/2)-((Gamma[(3/2)+n])/(Gamma[(1/2)+n])),tanα2->(1/2)-((Gamma[(3/2)+n])/(Gamma[(1/2)+n]))}]

Factor[(n^2)+n-2]

NζScatteringAngle[n_,x1_,ν_]:=If[Abs[x1]<=(10^(-6)),If[n==1,-45,0],ζScatteringAngle[n,x1,ν]];

Compute scattering angles, Table II of Farran (1951)

angles=Table[Chop[NζScatteringAngle[n,x1,1./3.]],{x1,0.,5.,0.2},{n,0,9}] ;TableForm[Round[100 angles]/100.,TableHeadings->{Table[x1,{x1,0.,5.,0.2}],Table[n,{n,0,9}]}]

angles=Table[Chop[NζScatteringAngle[n,x1,1./3.]],{x1,0.,0.2,0.01},{n,0,9}] ;TableForm[Round[100 angles]/100.,TableHeadings->{Table[x1,{x1,0.,0.2,0.01}],Table[n,{n,0,9}]}]

Exact scattering amplitude

cExact[n_,x1_,x2_,x3_,ρ1_,ρ3_,P0_]:=Module[{tanα1,tanα2,tanα3,tanβ3,tanδ3,tanζ1,tanΦ,tanη,η,c},
tanα1=If[x1==0,-n,tanα[n,x1]];
tanα2=If[x2==0,-n,tanα[n,x2]];
tanα3=If[x3==0,-n,tanα[n,x3]];
tanβ3=tanβ[n,x3];
tanδ3=tanδ[n,x3];
tanζ1=If[x1==0,If[n==1,-π/4,0],((-(x2^2))/2) ((((tanα1)/(tanα1+1))-(((n^2)+n)/(tanα2+(n^2)+n-1-(x2^2)/2)))/(((2tanα1+(n^2)+n-(x2^2)/2)/(tanα1+1))-((((n^2)+n)(tanα2+1))/(tanα2+(n^2)+n-1-(x2^2)/2))))];
tanΦ=((-ρ3)/(ρ1)) tanζ1;
tanη=tanδ3 ((tanΦ+tanα3)/(tanΦ+tanβ3));
η=ArcTan[tanη];
c=-P0 (2n+1) Re[ (((-I))^(n+1)) Sin[η] Exp[I η]];
Return[c]];

Sample calculations

15:3 Begostone in water sample calculations

Material properties

Properties for solid medium (cylinder)

ρ1=1995.; c1=4159.; c2=2319.; μ = ρ1 (c2^2); λ = ρ1 (c1^2) - 2 μ; ν=(λ/(2(λ+μ))); R=3.5 (10^(-3));

Properties for fluid medium

ρ3=1000.; c3=1500.;

Incident plane wave at mid frequency (1 MHz)

Compute series expansion terms

Set period, pulsation, wave numbers, amplitude

Clear[a,b,c];
T=1.0 (10^(-6)); ω=((2 π)/T); k1=(ω/c1); k2=(ω/c2);k3=(ω/c3); P0=40 (10^6); x1 = k1 R; x2=k2 R; x3 = k3 R;

Compare analytical and numerical solution to boundary condition matching. Analytical solution

sol[10,x1,x2,x3,R,ω,ρ1,ρ3,c2,ν,P0]

Numerical solution

LinearSolve[ M[10,x1,x2,x3,R,ω,ρ1,ρ3,c2,ν],f[10,x3,R,ω,ρ3,P0]]

Even though matrix is ill-conditioned, numerical solution closely matches analytical solution.

Compute expansion coefficients using analytical solution to boundary condition system. Compare analytical and numerical c coefficients of the scattered wave.

coef=Monitor[Table[sol[n,x1,x2,x3,R,ω,ρ1,ρ3,c2,ν,P0],{n,0,50}],n];

a[n_]:=a[n]=Re[coef[[n+1,1]]];
b[n_]:=b[n]=Re[coef[[n+1,2]]];
c[n_]:=c[n]=Re[coef[[n+1,3]]];

Table[{n,a[n],b[n],c[n],cExact[n,x1,x2,x3,ρ1,ρ3,P0]},{n,0,50}] // TableForm

Overall 2-norm of the error is small

Norm[Table[c[n]-cExact[n,x1,x2,x3,ρ1,ρ3,P0],{n,0,50}]]/51

Construct a graphical representation of the expansion coefficients. The a coefficients correspond to the scalar potential Ψ expansion, the b coefficients correspond to the vector potential A expansion.

acoef=Table[{n,Log10[Abs[a[n]]]},{n,0,50}];
bcoef=Table[{n,Log10[Abs[b[n]]]},{n,0,50}];
ccoef=Table[{n,Log10[Abs[c[n]]]},{n,0,50}];

ListPlot[{acoef,bcoef,ccoef},Joined->True,PlotMarkers->Automatic,Axes->False,Frame->True,FrameLabel->{"n","lg(coef)"},PlotRange->All,PlotLegends->{"lg((Subscript[a, n]))","lg((Subscript[b, n]))","lg((Subscript[c, n]))"}]

Displacement expansions

Above table and plot shows oscillatory expansion coefficients (Subscript[a, n]),(Subscript[b, n]). The displacement series use linear combinations of these coefficients

urcoef=Table[Log10[Abs[n(n+1) ((b[n])/R)j[n,x2] + a[n]k1 jp[n,x1]]],{n,0,50}] ;
uθcoef=Table[Log10[Abs[((a[n])/R)j[n,x1] +((b[n])/R)j[n,x2]+ b[n]k2 jp[n,x2]]],{n,0,50}] ;

ListPlot[{urcoef,uθcoef},Joined->True,Axes->False,Frame->True,FrameLabel->{"n","lg(coef)"},PlotRange->All,PlotLegends->{"lg((Subscript[u, r]) (n^th) coef)","lg((Subscript[u, θ]) (n^th) coef)"}]

Though the (Subscript[a, n]),(Subscript[b, n]) coefficients diverge with n, the linear combinations arising in the displacement expansions decay with n

urC=ur[R,0.,k1,k2,50];
Plot[Evaluate[Table[ur[rR*R,0.,k1,k2,L]/urC,{L,5,40,5}]],{rR,0.,1.},Axes->False,Frame->True,FrameLabel->{"r/R","(Subscript[u, r]) (m)"},PlotLegends->Table[L,{L,5,40,5}]]

Plot[Evaluate[Table[Log10[Abs[(ur[rR*R,0.,k1,k2,L]-ur[rR*R,0.,k1,k2,50])/ur[rR*R,0.,k1,k2,50]]],{L,5,45,5}]],{rR,0.,1.},Axes->False,Frame->True,FrameLabel->{"r/R","lg rel. err. (Subscript[u, r]) "},PlotLegends->Table[L,{L,5,45,5}]]

From the above error plot, series expansions with 40 terms have relative error less than (10^(-20)), those with 30 terms have relative error less than (10^(-10)).

Compare displacements obtained with various number of terms in series expansion

Δr=0.0001; Δθ=π/180.; rmin=0.00001;

urSphR=Monitor[Table[{θ,ur[R,θ,k1,k2,30]},{θ,0,2π,Δθ}],θ];
ListPlot[urSphR,Joined->True,Axes->False,Frame->True,FrameLabel->{"θ","(Subscript[u, r])(R)"},GridLines->Automatic]

uθSphR=Monitor[Table[{θ,uθ[R,θ,k1,k2,30]},{θ,0,2π,Δθ}],θ];
ListPlot[uθSphR,Joined->True,Axes->False,Frame->True,FrameLabel->{"θ","(Subscript[u, r])(R)"},GridLines->Automatic]

urSph=Monitor[Table[{r Cos[θ],r Sin[θ],ur[r,θ,k1,k2,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];
uθSph=Monitor[Table[{r Cos[θ],r Sin[θ],uθ[r,θ,k1,k2,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];

urSphPlot=ListContourPlot[Flatten[urSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(Subscript[u, r]) (m)"]

uθSphPlot=ListContourPlot[Flatten[uθSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(Subscript[u, θ]) (m)"]

uzSph=Monitor[Table[{r Cos[θ],r Sin[θ],uz[r,θ,k1,k2,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];

uzSphPlot=ListContourPlot[Flatten[uzSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(Subscript[u, z]) (m)"]

durdrSph=Monitor[Table[{r Cos[θ],r Sin[θ],durdr[r,θ,k1,k2,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];
durdθSph=Monitor[Table[{r Cos[θ],r Sin[θ],durdθ[r,θ,k1,k2,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];
duθdrSph=Monitor[Table[{r Cos[θ],r Sin[θ],duθdr[r,θ,k1,k2,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];
duθdθSph=Monitor[Table[{r Cos[θ],r Sin[θ],duθdθ[r,θ,k1,k2,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];

durdrSphPlot=ListContourPlot[Flatten[durdrSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(∂_r)(Subscript[u, r]) (m)",PlotRange->All]

durdθSphPlot=ListContourPlot[Flatten[durdθSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(∂_θ)(Subscript[u, r]) (m)",PlotRange->All]

duθdrSphPlot=ListContourPlot[Flatten[duθdrSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(∂_r)(Subscript[u, θ]) (m)",PlotRange->All]

duθdθSphPlot=ListContourPlot[Flatten[duθdθSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(∂_θ)(Subscript[u, θ]) (m)",PlotRange->All]

σrrSph=Monitor[Table[{r Cos[θ],r Sin[θ],σrr[r,θ,k1,k2,λ,μ,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];
σrθSph=Monitor[Table[{r Cos[θ],r Sin[θ],σrθ[r,θ,k1,k2,λ,μ,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];
σθθSph=Monitor[Table[{r Cos[θ],r Sin[θ],σθθ[r,θ,k1,k2,λ,μ,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];
σϕϕSph=Monitor[Table[{r Cos[θ],r Sin[θ],σϕϕ[r,θ,k1,k2,λ,μ,30]},{r,rmin,R,Δr},{θ,0,2π,Δθ}],r];

Stress expansions

σrr[0.01,0.,k1,k2,λ,μ,45]

σrrR=Table[σrr[R ,0.,k1,k2,λ,μ,5n],{n,1,8}];
σrrCL=Monitor[Table[{rR,σrr[R rR,0.,k1,k2,λ,μ,5n]/σrrR[[n]]},{n,1,8},{rR,0.000001,1.,0.001}],{n,rR}];

ListPlot[σrrCL[[2;;]],Joined->True,PlotLegends->{"n=10","n=15","n=20","n=25","n=30","n=35","n=40"}]

σrrSphPlot=ListContourPlot[Flatten[σrrSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(Subscript[σ, rr]) (N/(m^2))"]

σrθSphPlot=ListContourPlot[Flatten[σrθSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(Subscript[σ, rθ]) (N/(m^2))",PlotRange->All]

σθθSphPlot=ListContourPlot[Flatten[σθθSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(Subscript[σ, θθ]) (N/(m^2))",PlotRange->All]

σTSph=Table[{σrrSph[[i,j,1]],σrrSph[[i,j,2]],
σT[σrrSph[[i,j,3]],σrθSph[[i,j,3]],σθθSph[[i,j,3]],σϕϕSph[[i,j,3]]]},{i,1,Dimensions[σrrSph][[1]]},{j,1,Dimensions[σrrSph][[2]]}];
σCSph=Table[{σrrSph[[i,j,1]],σrrSph[[i,j,2]],
σC[σrrSph[[i,j,3]],σrθSph[[i,j,3]],σθθSph[[i,j,3]],σϕϕSph[[i,j,3]]]},{i,1,Dimensions[σrrSph][[1]]},{j,1,Dimensions[σrrSph][[2]]}];

σTSphPlot=ListContourPlot[Flatten[σTSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(Subscript[σ, T]) (N/(m^2))",PlotRange->All]

σCSphPlot=ListContourPlot[Flatten[σCSph,1],ColorFunction->"LightTemperatureMap",PlotLegends->Automatic,PlotLabel->"(Subscript[σ, C]) (N/(m^2))",PlotRange->All]

SetDirectory["/home/mitran/papers/2017/LithotripsyValidation/notebooks"];
DumpSave["FarranSpheres.mx","Global`"];