Torsional Vibrational Modes of an Elastic Sphere Last updated: January 21, 2008
Torsional Vibrational Modes of an Elastic Sphere
Torsional modes (also known as Love waves and called "toroidal" in seismology) have divergence free displacement fields. The lower frequency torsional modes of an isotropic elastic sphere are calculated exactly. Additional details are added to an earlier derivation [Zhen 2000].

1. Introduction

   The classic problem of the vibrational modes of an isotropic continuum elastic sphere has been extensively studied [Lamb 1882] [Love 1944] [Auld 1973] [Nishiguchi and Sakuma 1981] [Bastrukov 1994] [Graff 1994] [Bastrukov et al. 1998] [Ye 2000 pdf ::], and an accurate list of the exact frequencies has recently been reported [Murray 2002a]. This problem and its generalizations have applications in geophysics and seismology [L. Hanyk et al. 1998 pdf ::]. Molecular dynamical simulations of a Lennard Jones like crystalline solid [Murray 2002b] provide numerical evidence to support the correctness of the exact continuum analysis.
   Consider a freely suspended isotropic elastic sphere not subject to body forces. The sphere has radius R, uniform density ρ, shear modulus μ and Poisson ratio ν. These parameters fully describe the problem. If the sphere is disturbed, it will begin to vibrate. A given eigenfunction has frequency ω (in radians per second). The speed of transverse waves (shear waves) is

Cs = sqrt(μ/ρ) [eq.1.1]

while the speed of longitudinal (compression) waves is

Cl = sqrt((λ+2μ)/ρ) [eq.1.2]

   It is convenient to express normal mode frequencies in terms of the dimensionless variable

η = ω R / Cs [eq.1.3]

and to note that the wavevector, ks (in inverse metres) of shear (transverse) waves is

ks = ω / Cs = η / R [eq.1.4]

while the wavevector of longitudinal (compression) waves is:

kl = ω / Cl [eq.1.5]

   For a normal mode with frequency ω, the displacement field u(r,t) in the sphere is

u(r,t) = u(r) exp(i ω t). [eq.1.6]

The vector field u(r) can with generality be expressed as the sum of the gradient of a scalar field and the curl of a divergence free vector field.

u(r) = div Φ + curl A [eq.1.7]

2. Pulsating, Spheroidal and Torsional Modes:

   It is at this point that the distinction is made between two classes of normal modes, "spheroidal" (also know as Rayleigh waves) and "torsional" (also known as Love waves). The simplest and easiest to visualize kind of vibrations are where u = grad( Φ(r) ), where r = ||r||, and it is possible to get a closed form expression for the frequencies of these modes, which are called the "pulsating modes" [Love 1944][equation (24) in Ye 2000]

  tan(ξ)
-------
   ξ		  =      1
----------
1 - η2/4		 [eq.2.1]
where ξ = (Cs/Cl) η.
   It is completely natural to guess that other normal modes can be found of the more general form u = grad( Φ(r) ). However, it turns out that there are no other modes of this form apart from the "pulsating" ones! The problem is that the boundary conditions for stress cannot be satisfied at all points on the surface of the sphere. However, it is possible to satisfy the boundary conditions for modes (called the "spheroidal" modes) of the form:

u(r,θ,φ) = Al grad [ j(l, kl r ) P(l,cos(θ)) ]
+ Bl curl curl [ r j(l, ks r ) P(l,cos(θ)) ]
[2.2]
for l=0, 1, 2, 3, ... . Here j(l,x) are spherical Bessel functions of the first kind and P(l,x) are the Legendre polynomials. The ratio of constants Bl/Al must be chosen so as to match the boundary conditions. The "pulsating" modes are the l=0 case of the "spheroidal" modes.
   The other class of modes (called the "torsional" modes) corresponds to choosing the displacement field, (given in spherical coordinates by (ur,uθ,uφ) ) to have the form:

ur = 0
uθ = 0[2.3]
uφ(r,θ,φ) = - j(l,ks r) d/dθ ( P(l,cos(θ) )

where spherical coordinates are related to Cartesian coordinates by x = r sin(θ) cos(φ), y = r sin(θ) sin(φ) and z = r cos(θ). The wavevector ks = η/R, where the dimensionless frequency η takes on an eigenvalue for each mode.
   It is possible to back up one step further. Suppose that u = curl(r Q(r) ) where vector field r = (x,y,z) and scalar field Q = j(l,ks r) P(l,cos(θ)), (or to derive the result in even greater generality suppose Q(r,θ) = R(ks r)g(θ) for arbitrary functions R() and g(). Steps later in this paper show that R(x) must be the spherical Bessel functions of the first kind and g(θ) must be related to the Legendre polynomials.) Applying the formula for curl in spherical coordinates to r Q, we get the above expression for u.

3. Solving the Navier-Stokes Equation:

   The differential equation for the displacement field [after Ye 2000 equation (3) ] is:

(λ+2μ) grad(div u) - μ curl(curl u) = ρ (d2 u/dt2 ) [3.1]

where λ is the first Lamé constant. Note that div u = 0 for the above torsional form. Then, assuming that u is proportional to exp(i ω t),

curl(curl u) = (ks)2 u [3.2]

Where ks is the wavevector for shear waves. Note that (ρ ω2 / μ) = (ω/Cs)2 = (ks)2
   The full expression for (curl u) = B = (Br,Bθ,Bφ) in spherical coordinates is:
Br =     1
---------
r sin(θ) 		[  d
---
dθ		 (sin(θ) uφ) -   d
---
dφ		 uθ ]

i.e.
Br =     1
---------
r sin(θ) 		[  d
---
dθ		 (sin(θ) uφ) ]

Bθ= (1/r) [ (1/sin(θ)) (d/dφ) ur - (d/dr)(r uφ) ]
i.e.
Bθ= (-1/r) [ (d/dr)(r uφ) ]

Bφ = (1/r) [ (d/dr)(r uθ) - (d/dθ) ur ]

where the terms that are zero are highlighted in red. Let (curl curl u) = (curl B) = D = (Dr,Dθ,Dφ). Once again using the formula for curl in spherical coordinates
Dr =     1
---------
r sin(θ) 		[  d
----
dθ		(sin(θ) Bφ) -   d
---
dφ		 Bθ ]

Dθ= (1/r) [(1/sin(θ)) (d/dφ) Br - (d/dr)(r Bφ) ]

Dφ =  1
---
 r		 [  d
---
dr		 ( r Bθ )  -   d
---
dθ		 ( Br ) ]
substituting...

Dr = (1/(r sin(θ))) [ (d/dφ) ((1/r) [ (d/dr)(r uφ) ]) ]

Dθ= (1/r) [(1/sin(θ)) (d/dφ){ (1/(r sin(θ))) [(d/dθ)(sin(θ) uφ) ] } ]

Dφ =  1
---
 r		 [  d
---
dr		 ( r (-1/r) [ (d/dr)(r uφ) ] )  -   d
---
dθ		 (     1
---------
r sin(θ) 		[  d
----
dθ		(sin(θ) uφ) ] ) ]
Both Dr and Dθ vanish because there is no φ dependence anywhere. We are finally left only with

Dr = 0
Dθ=0
Dφ =  -   1
---
 r		 [  d
---
dr		 ( r (1/r) [  d
---
dr		 (r uφ) ] ) ]
 -    1
-----
 r*r		 [  d
---
dθ		 (    1
-------
sin(θ) 		[  d
---
dθ		(sin(θ) uφ) ] ) ]
Now, since D = curl curl u = (ks)2 u, we simply equate the φ component of each.
 -   1
---
 r		 [  d
---
dr		 (  d
---
dr		 (r uφ) ) ]
 -    1
----
r2		 [  d
---
dθ		 (    1
-------
sin(θ) 		[  d
---
dθ		(sin(θ) uφ) ] ) ]
= (ks)2 uφ
The next step is to substitute for uφ after inserting the assumed form:
 -   1
---
 r		 [  d
---
dr		 (  d
---
dr		 (r ( (j(l,ks r)  d
---
dθ		 P(l,cos(θ)) ) ) ) ]
 -    1
----
r2		 [  d
---
dθ		 (    1
-------
sin(θ) 		[  d
---
dθ		 (sin(θ) ( (j(l,ks r)  d
---
dθ		 P(l,cos(θ)) ) )] ) ]
= (ks)2 ( (j(l,ks r)  d
---
dθ		 P(l,cos(θ)) )

Now divide every term by j(l,ks r)(d/dθ) P(l,cos(θ)) to get

-       1
-------------
 r j(l,ks r)		 [ d2
-----
dr2		 (r j(l,ks r) ) ]
-  (1/[(d/dθ)P(l,cos(θ))]) (1/(r2)) [  d
---
dθ		 [   1
-------
sin(θ)		 [  d
---
dθ		 ( sin(θ)  d
---
dθ		 P(l,cos(θ) ) ] ] ]
= (ks)2

Finally, multiply each term by r2.

-      r
-----------
 j(l,ks r)		 [  d2
----
dr2		 (r j(l,ks r) ) ]
-  (1/[(d/dθ)P(l,cos(θ))]) [  d
---
dθ		 [   1
-------
sin(θ)		 [  d
---
dθ		 ( sin(θ)  d
---
dθ		 P(l,cos(θ) ) ] ] ]
= (ks r)2
And next move one term to the other side.
-      1
---------------------
(d/dθ)P(l,cos(θ))		 [  d
---
dθ		 [   1
-------
sin(θ)		 [  d
---
dθ		 ( sin(θ)  d
---
dθ		 P(l,cos(θ) ) ] ] ]
=      r
-----------
 j(l,ks r)		 [  d2
----
dr2		 (r j(l,ks r) ) ] + (ks r)2
   The two sides are equal to one another. But the green side depends on θ and not on r. The red side depends on r but not on θ. The only possible resolution of this is that neither side depends on either r or θ. They are both equal to the same constant. A propitious choice for this turns out to be l(l+1). This turns the previous single equation into two separate equations, the first of which is:
0 =  r [  d2
----
dr2		 (r j(l,ks r) ) ] + [ (ks r)2 - l(l+1) ] j(l,ks r)
   This is the Spherical Bessel differential equation. So it is certainly satisfied by the j(l,k r) which are Spherical Bessel functions of the first kind.
   The other equation we get is:
 d
---
dθ		 [   1
-------
sin(θ)		 [  d
---
dθ		 ( sin(θ)  d
---
dθ		 P(l,cos(θ) ) ] ] = - l(l+1) (d/dθ)P(l,cos(θ))
Take the indefinite integral of both sides over θ, and multiply by sin(θ),
 d
---
dθ		 ( sin(θ)  d
---
dθ		 P(l,cos(θ) ) = - l(l+1) sin(θ) P(l,cos(θ))
   This is the differential equation for Legendre polynomials, so it is certainly satisfied by P(l,cos(θ)).

4. Boundary Conditions

   The next step is to consider the boundary conditions at the surface of the sphere. There is no external force applied to the surface of the sphere, and therefore some of the components of the stress tensor must vanish there. It remains then only to find the values of ks for which the boundary condition is satisfied. These correspond to the frequencies of the modes of vibration.
   The stress tensor σ is given, in Cartesian coordinates, by

λ div u + 2 μ dux/dxμ(dux/dy+duy/dx)μ(dux/dz+duz/dx)
μ(dux/dy+duy/dx)λ div u + 2 μ duy/dyμ(duy/dz+duz/dy)
μ(dux/dz+duz/dx)μ(duy/dz+duz/dy)λ div u + 2 μ duz/dz

where λ is the first Lamé constant and μ is shear modulus.
   Consider some point on the surface. Suppose that the X-axis is normal to the surface at that point, and Y and Z are tangent to the surface. Since no external force acts on the surface, σXX = 0, σXY = 0 and σXZ = 0. I will use an X Y Z axis system that is tilted up by angle δ. Normal cartesian coordinates are x, y and z.

x = r sin(θ) cos(φ)
y = r sin(θ) sin(φ)
z = r cos(θ)

X = cos(δ) x + sin(δ) z
Y = y
Z =-sin(δ) x + cos(δ) z

x = cos(δ) X - sin(δ) Z
y = Y
z = sin(δ) X + cos(δ) Z

so that

r sin(θ) cos(φ) = cos(δ) X - sin(δ) Z
r sin(θ) sin(φ) = Y
r cos(θ) = sin(δ) X + cos(δ) Z

r = sqrt(x2+y2+z2)
θ = arctan(sqrt(x2+y2)/z)
φ = arctan(y/x)

   I wrote a short and easy to follow C++ computer program (wetdog4.cpp C++ source code ) to numerically calculate u and all of its derivatives at a given value of X, Y, Z.
   Then, I chose an arbitrary value of δ and evaluate the strain components at points along the X axis (Y=Z=0). It always seems to be the case that the diagonal strain components are zero. The XZ component of strain was also always zero. The XY strain varies with X. I then numerically searched for values of X at which the XY strain is zero. This has to be repeated for each value of l. I tried l from 1 to 5. For numerical reasons, I only searched up as far as η = 12.

Table I.
Torsional mode frequencies calculated using wetdog4.cpp
l η
(1st root)		 η
(2nd root)		 η
(3rd root)
15.7649.096-
22.5027.13710.514
33.8658.44411.910
45.0949.713-
56.26610.951-

   The values of η turn out to be exactly the same as the frequencies for the torsional modes calculated earlier [Murray 2002 link to article ]. I verified this to three digits of precision for quite a number of different roots and different values of l. The YZ strain component is not zero, but that is not required by the boundary conditions.

5. Duality of Forms for Displacement Field

   There is a duality relationship between the zero divergence part of the spheroidal modes and the torsional modes. Suppose that a zero divergence displacement field u satisfies

(λ+2μ) grad(div u) - μ curl(curl u) = ρ (d2 u/dt2 ) [3.1]

so that, also having assumed frequency ω = Cs ks,

μ curl(curl u) = (ks)2 u [3.2]

Take the curl of both sides of this equation.

μ curl(curl(curl u)) = (ks)2 curl u

   Thus, if zero divergence u is a solution of the equation [3.1] and [3.2], then u2 = curl u also satisfies equation [3.2]. But since div (curl( A ) ) = 0 for any vector field A, u2 satisfies equation [3.1] as well. It would also be good to check the boundary conditions for stress at the surface of the sphere. This turns out to be a problem. The stress of the dual form is not the same, and so the boundary conditions are not met. To see that this is a duality, recall that

curl (curl ( A ) ) = - laplacian A

for any divergence free vector field A. And also note from equation [3.2] that the curl of the curl of a torsional mode just results in the same torsional mode multiplied by a constant. From this, we see that just as the curl of a torsional mode of order l gives the spheroidal zero divergence form of order l, in the same way the curl of a spheroidal zero divergence form of order l gives the torsional mode of order l.
References:
H. Lamb, Proc. London Math. Soc. 13, 187 (1882).)
A. E. H. Love, A Treatise on the Mathematical Theory of Elasticity, (Dover, New York, 1944).
B. A. Auld, Acoustic Fields and Waves in Solids, (John Wiley & Sons, New York, 1973).
N. Nishiguchi and T. Sakuma, Sol. Stat. Comm. 38, 1073 (1981). [Modal series expansion]
Sergey I Bastrukov, Phys. Rev. E 49, 3166 (1994).
F. K. Graff, Wave Motion in Elastic Solid, (Ohio University Press, Ohio, 1994).
Bastrukov et al., Physica A 250, 435 (1998).
Zhen Ye, "On the Low Frequency Elastic Response of a Spherical Particle" Chinese Journal of Physics Vol. 38, pages 103-110, (2000). link to pdf ::
Daniel B. Murray,"Vibrational Frequencies of an Elastic Sphere" 2002 ( link to article )
Daniel B. Murray,"Molecular Dynamical Estimates of Vibrational Frequencies of a Crystalline Sphere" ( link to article )
Daniel B. Murray, "Velocity Fields of Vibrational Modes of an Elastic Sphere" 2002 ( link to article )
L. Hanyk, C. Matyska and D. A. Yuen "Initial value approach for viscoelastic responses of the earth's mantle" 1998. pdf ::
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