First created: Mar. 18, 2004     Updated: Wed. Apr. 7, 2004
Elastic wave scattering from a solid spherical inclusion
Support information and on-line reference for "cross" paper

The calculation required to generate these curves assumes that the incident wave is a spherical Bessel function of the first kind, rather than a Hankel function. That means that the incident wave is not a travelling wave.

Letters A, B and C represent the amplitudes of the various parts of the wave. The amplitude of the incident wave is 1, which is plotted as a grey horizontal line on each graph. These amplitudes are complex numbers. The absolute value is what is plotted.

Note that the horizontal axis for these figures shown below is kloa, which is a C++ variable name chosen to represent the dimensionless quantity klo a, where klo is the longitudinal wavevector in the matrix and a is the radius of the sphere.

These three figures are based on Figure 5 in Kafesaki and Economou PRB52 13317 1995. This is the situation of an steel sphere inside an epoxy matrix. The l=0 curve doesn't have any sharp structure. The l=1 curve does not seem to show the peak in Kafesaki, but I think this is because of the way the waves are normalized. I think there is a factor of k missing. A quick experiment (not plotted here) seemed to show this is right. The l=2 curve has the desired peak at kloa = 3.3 in agreement with Kafesaki.


These two figures are based on Figure 6 in Kafesaki and Economou PRB52 13317 1995. This is the situation of an epoxy sphere inside a steel matrix. The l=0 curve shows quite a sharp feature around kloa = 1.8 in apparent agreement with Kafesaki, although Kafesaki left the peak unlabelled. The l=1 curve has three sharp features in agreement with those peaks labelled in Kafesaki. The l=2 curve also seems to match since it has a peak at kloa = 1.7, as shown in Kafesaki. However it has a second smaller feature at kloa = 1.3 which is not visible in Kafesaki most likely because of the large l=1 peak in that region.

J. Sessarego, J. Sageloli, R. Guillermin, and H. Überall
"Scattering by an elastic sphere embedded in an elastic isotropic medium"
J. Acoust. Soc. Am. 104, pages 2836-2844 (1998)

   In their formalism they have ordinary ψ and bold ψ, which happily does work correctly on my web browser. ψ is a scalar quantity while ψ is a vector field. φ is also a vector field. They did not use Φ as some other papers have done. Literal quotes from their paper are in blue color. u is the displacement field. They say that the incident field is a plane compressional wave.
u = ∇ φ + ∇ × ψ
(1)
where φ and ψ satisfy the equations
∇2φ + kL2φ = 0       ∇2ψ + kT2ψ = 0
(2)
However, the right side of equation (2) is incorrect since ψ is not in boldface. It should instead read:
∇2ψ + kT2ψ = 0
In a problem with spherical symmetry the vector potential ψ [sic - should be ψ] has only one component and can be written in the form: ψ = (0, 0, ψ) where ψ is a scalar.
    In these conditions the scalar potentials φ and ψ satisfy the two equations:
∇2φ + kL2φ = 0       ∇2ψ + kT2ψ = 0
(3)
    The problem with this is that they are using the wrong form for the zero divergence portion of the displacement field. What they should be using is  ∇ × ∇ × (r χ)  where χ is a scalar field. It is traditional to use this form together with  ∇ × (r ψ)  when constructing a general solution of the equation of motion in spherical coordinates. Their choice of   ∇ × (0, 0, ψ)   is a correct solution of the equation of motion, assuming that ψ satisfies ∇2ψ + kT2ψ = 0, but if ψ is a function of r and θ only then   ∇ × (0, 0, ψ)   has displacement purely along the φ direction. So it is completely unsuitable for mixing with the zero curl term   ∇ Φ   so as to match boundary conditions.
March 27, 2004:
  This is about numerical verification of decomposition formula for scalar plane wave:
ei k r cos(θ) =  Σ
l = 0 		il (2 l + 1) jl(k r) Pl(cos(θ))
To verify this formula, and also test my numerical method, I used a program plane1.cpp which uses Monte Carlo to estimate the coefficient for each order l and compare it to the expected value of il (2 l + 1). I consider a spherical region of radius Rmax. I pick some value of k. Then I evaluate the inner product of the plane wave with a given spherical wave and divide by the inner product of the spherical wave with itself. The results worked out quite well.

   Next I numerically verified the expansion of a longitudinal plane wave in p-wave partial spherical waves. The form of the incident wave is as in the cross paper, [eqLP2] Eq. (21).
(0, 0, A exp(i kL r cos(θ))) = ∇ [ (  1
---
kL 		)  Σ
l = 0 		A i(l-1) (2 l + 1) jl (kL r) Pl (cos(θ)) ]
In other words, the partial p-wave amplitudes are
Bl = A i(l-1) (2 l + 1)
as in [eqLP6] equation (26) in the cross paper.
   The C++ program is PLongPwB.cpp. It worked quite well for l up to 7. I also was able to reproduce the expected factor of i(l-1) (2 l + 1).
C++ source code: PlongPwB.cpp

For a transverse oncoming wave, the case of s-wave amplitudes is checked by:
C++ source code: PTranSwH.cpp

   The case of decomposition of a transverse plane wave into spherical partial waves is handled for t-waves by C++ program PTranTwI.cpp. I am using x-polarization.
C++ source code: PTranTwI.cpp

C++ Monte Carlo programs for checking partial wave power (P in watts) formulae:
Longitudinal incident wave:
fluxp0.cpp - p-waves [eqEF18] (80)
fluxs0.cpp - s-waves [eqEF19] (81)
Transverse incident wave:
fluxp1.cpp - p-waves [eqEF20] (82)
fluxs1.cpp - s-waves [eqEF21] (83)
fluxt1.cpp - t-waves [eqEF22] (84)
Daniel Murray
Associate Professor
Math, Stats & Physics Unit
University of British Columbia - Okanagan
Kelowna, BC, Canada
daniel "dot" murray "at" ubc "dot" ca

For a list of related articles click on the link.
Hosted by www.Geocities.ws

1