Journal of Sound and <ibration (2000) 000(0), 000}000
doi:10.1006/jsvi.2000.2963, available online at http://www.idealibrary.com on

THE xyzr ALGORITHM SPECIALIZED FOR
EIGENVIBRATION PROBLEM OF BORED AND
LAMINATED OBJECTSAKIRA YONEDA
Institute for Study of the Earth1s Interior, Okayama ;niversity, Misasa, ¹ottori 682-0193, Japan(Received 27 July 1999, and in ,nal form 8 February 2000)
We can determine elastic properties of an object from its eigenvibration frequencies, aslong as a numerical calculation algorithm is available for the eigenvibration frequencies. The
xyz algorithm developed by Visscher et al. [1] is an extremely versatile method applicableeven for anisotropic, inhomogeneous, and irregular-shaped objects, with the exception ofbored and laminated systems. I examined the characteristics and performance of the original
xyz algorithm, and succeeded in developing the xyzr algorithm for bored and laminatedsphere and cylinder. The accuracy and convergency of the xyzr algorithm are examined bycomparison with analytical solutions; the method is accurate enough for us to develop a newexperimental technique to obtain the elastic properties of soft materials accommodated ina solid container.
 2000 Academic Press

1. INTRODUCTION

The eigenvibration problem has been analyzed by means of the Rayleigh}Ritz method
based on Hamilton's principle. Holland [2] solved the eigenvibration problem for a cubewith cubic crystallographic symmetry by using trigonometric functions as the basefunctions. Demarest [3] improved convergence of the same problem by introducing
Legendre polynomials as base functions. Ohno [4] expanded Demarest's solution toa parallelepiped rectangle with orthorhombic crystallographic symmetry. Senoo andNishimura [5] obtained a numerical solution for a cylinder by combining the Rayleigh}Ritz
method and the "nite element method by using Legendre polynomials for base functions as
well. However, further progress has not yet been achieved in this analyzing scheme owing to
di$culty in manipulating Legendre polynomials.
The xyz algorithm by Visscher et al. [1] provided a breakthrough in this "eld.
They showed that simple power functions in Cartesian co-ordinates can be used forthe base functions (Appendix A). The xyz algorithm enables an easy analysis of
eigenvibration, even for objects with internal inhomogenenity, elastic anisotropy, andirregular shape.
The versatility of the xyz algorithm inspired me to devise a new experimental technique
to investigate the rheological properties of liquid and soft materials by measuring
eigenvibration frequencies of a container "lled with the target material. However, I foundthat the xyz algorithm is not so accurate for an object with inner and outer interfaces suchas a spherical shell or tube (see section 3).
Through examining the characteristics of the xyz algorithm, I developed the xyzr
algorithm as a solution for a spherical shell or a tube. Here, I discuss "rst the characteristics
0022-460X/00/000000#00 $35.00/0
RAVINDRA BG yamuna 2000 Academic Press
JSV 20002963
DATE ....................
DISC USED YES ! NO !

2

A. YONEDA

Figure 1. An example showing how a sinusoidal function is "tted well by power functions. Note that odd powers+x, x, x, x, were used in the "tting because the sine function (thick line) is an odd function. Thin lines are
the resultant polynomial functions after "tting; N is the maximum power used in the "tting. For instance,N"5 speci"es that the set of +x, x, x, are used in the "tting. We can see that the "ttable range increases withincreasing N.

of the xyz algorithm by comparison with analytical solutions, and then introduce the xyzr
algorithm and show its performance.
In the following discussion, please note the di!erences among solid, bored, and laminatedobjects with spherical or cylindrical symmetry. Namely, solid sphere and cylinder do nothave an empty region inside, whereas, bored objects are identical with a shell and a tube andare empty inside. Laminated objects are two layered without an empty region inside.

2. CHARACTERISTICS OF THE xyz METHOD

It is quite curious why the xyz method works so widely and accurately in various
eigenvibration problems. A displacement function for standing waves, u (x, y, z), is generallyGgivenu (x, y, z)"e SR sin (lx#p ) sin (my#p ) sin (nz#p ) e ,
G
V
W
X G

(1)

where u is the angular frequency, and e a directional unit vector in the ith direction
G(i"1, 2, 3) in the Cartesian co-ordinate. Any eigenvibration mode can be expressed asa superposition of the above equations. Fundamental modes of eigenvibration are speci"edby wave numbers standing in a vibrating body, or numbers of nodal surfaces. Figure 1
shows how "nite power functions approximate well a sinusoidal function, or a standing
wave excited in a body. We can judge that a mode of w "0)5 (w : wave number over sample
LL
dimension) may be calculated accurately even at N"3; identically, for w "1)0, at N"5L
and for w "1)5, at N"7 respectively. A simpli"ed relation is that N is required to be 5}6
Ltimes that of w .
L
I con"rmed the above discussion by examining eigenvibrations for an isotropic sphere,
whose eigenvibration modes can be calculated analytically. In the following discussion,please note that physical properties of an object are always "xed at p"0)25 (p the Poisson
ratio), or l "(3l (l the longitudinal wave velocity; l the shear wave velocity), andN
Q N
Q
analytical eigenfrequencies are always normalized by l /d (d the diameter of sphere orQ  cylinder).
In spherical co-ordinates, the analytical solutions of displacement are expressed asa product of the spherical Bessel function, j ( j (x)"J (x)/x; J being the "rst kind of Bessel
L
L
JSV 20002963

3

THE xyzr ALGORITHM

Figure 2. An example showing how the "rst kind of spherical Bessel function is "tted well by power functions.
Thick lines are j (r) and j (r), in the displacement functions of S and ¹ respectively. Procedures are similar to


O 
 L
that in Figure 1. Note that N"3, 4, and 5 specify the set of power functions used, or +r, r,, +r, r, r,, and+r, r, r, respectively.

function), the Legendre polynomial function, P, and the trigonometric functions.
Displacement functions of some fundamental modes are as follows (Appendix B)
hJ (hr)S : ; "! , ; "0,
? 
P
F
(r
¹:
 

; "0,
(

< "0,
< "0,
P
F
J (kr) +(n#1!m)PK (cos h)!(n#1) cos hPK (cos h),
L>
L< " 
e K(.
(
sin h
(r

(2)

(3)

The "tness of power functions to the spherical Bessel functions is shown in Figure 2. The"rst and second extrema of j (h) satisfy the boundary condition of S and 
S respectively. Therefore, accurate frequencies are expected to be given at N"3 and 5, 
respectively, for those non-angular modes. This prediction is consistent with theresults shown in Table 1. On the other hand, for the ¹ modes, it takes 43 to L
approximate their radial term. Considering that it requires an additional n3 for the
angular term, the eigenfrequencies of ¹ and ¹ are expected to be given accurately  
at N"6 and N"11, respectively, as shown in Table 1. Although only the four modes aredemonstrated here, I have checked that the present discussion is exactly valid for othermodes.
Before proceeding the numerical evaluation for bored and laminated systems, please notethat the inner}outer radius ratio is always "xed at r /r "0)5 as a representative case,G 
where the superscripts i and o refer to the inner and outer regions respectively. As
shown for the case of k"0 in Table 2, the xyz algorithm is not as accurate for a boredsphere, or a spherical shell, as it is for a solid sphere (see Table 1). In this case, a second
kind of Bessel function, Y, is introduced in the solutions as explained in Appendix B.
In general, y ("Y(x)/x) is hard to be "tted using only positive powers; Figure 3 shows
the situation and the necessity of introducing inverse powers to yield better "tting of y.
Therefore, I have devised a new algorithm for the bored body as described in the next
section.
JSV 20002963

4

A. YONEDA

TABLE 1Frequencies of the xyz algorithm for some selected modes in a homogeneous
isotropic solid sphere. Normalized analytical frequencies are shown in theparentheses below mode indices. The numerical frequencies are normalized
by the analytical one of each modeS (1)4133)

N"3: 1)0009
N"5: 1)0000¹
 (0)79614)

N"3: 1)0578
N"5: 1)0002
N"7: 1)0000

S (3)3403)

N"5: 1)0007
N"7: 1)0000¹
 (1)8346)

N"5: 1)0160
N"7: 1)0003
N"9: 1)0000

S (0.84030)

N"3: 1)0037
N"5: 1)0000¹
 (2)7120)

N"9: 1)0065
N"11: 1)0001
N"13: 1)0000

TABLE 2
Normalized numerical frequencies by the xyz and xyzr algorithms for someselected modes in a bored sphere or a shell with rG /r"0)5. Note that the k"0
corresponds to the xyz algorithm. Normalized analytical frequencies are shownin the parentheses below mode indices

N

k

S
 
(1)1228)S (062139)¹ (0)77520)

7

0
3
5
7
0
3
5
7
0
3
5
7
0
5
7
9
0
5
7
9

1)01061)00001)00001)00071)00231)00001)00001)00001)00041)00001)00001)00001)00011)00001)00001)00001)00001)00001)00001)00001)01941)01051)00371)01971)01151)00421)00061)00001)00701)00121)00011)00001)00321)00031)00001)00001)00121)00001)00001)00001)00221)00141)00081)06121)00091)00011)00011)00041)00031)00001)00001)00001)00011)00001)00001)00001)00001)00001)00001)0000

9

11

13

15

3. THE xyzr ALGORITHM FOR A BORED BODY AND A LAMINATED BODY

Assume a spherical shell with inner and outer radii of r and r respectively. I modi"edGbase functions as
xJ yKzL
u"
e,
H
G
rI
JSV 20002963(4)

THE xyzr ALGORITHM

5Figure 3. A demonstration showing why inverse powers are required in the "tting of the second kind of
spherical Bessel function. Thick lines are j (r) and y (r). The "tted interval is 2(r(5. The dash-dotted linesare the resultant polynomial functions in terms of the set of +r\, r\, r, r,, while the dotted lines are+r, r, r, r,. We can see that the set of +r, r, r, r, is not suitable to "t y (r) accurately, while the set of +r\, r\,r, r, gives reasonable "tting both for j and y functions.

Figure 4. Schematic illustration of various base functions. The numbers correspond to those de"ned inequations (6)}(8).

Figure 5. Comparison of convergency between the xyz (dotted lines with open symbols) and the xyzr (solid lines
with solid symbols) algorithms for some selected modes. The data in Table 2 are used. Note that k is chosen to be&N/2 for the xyzr algorithm: - - 䊐 - - *䊏* S ; - - -䉭-- - - - 2䉱2 S ; - -䊊- - - 2䊉2 ¹ .
 
 
 
JSV 20002963

6

A. YONEDA

where r"(x#y#z) and k is an adjustable parameter to yield the best result. Thisxyzr base function is su$cient for a spherical shell, precluding the co-ordinate origin ofr"0. Table 2 shows the accuracy of the xyzr algorithm in comparison with that of the xyz
algorithm for the special case (k"0) of the xyzr algorithm. We can recognize a signi"cantimprovement using the xyzr algorithm. In my experience, the best performance is usually
given at k&N/2; this suggests that the symmetric set of power functions, +r\I,
r\I\, 2 , rI\, rI,, is usually most advantageous for "tting j(r) and y(r) simultaneously.
In any case, results are not so sensitive to small variations in k.
The xyzr algorithm in cylindrical co-ordinates are the same as those in sphericalco-ordinates, except for the de"nition of r in equation (4). Obviously, r"(x#y) is theright choice in the cylindrical co-ordinate; the conversion rule is simply given asr"(x#y#z)Pr"(x#y).

(5)

For a laminated body, we can expect that the xyz and the xyzr algorithms are useful inthe inner and the outer regions respectively. When a laminated sphere is composed ofmaterials with di!erent elastic constants, some derivatives of displacement, or strains, must
be discontinuous at the interface, because there the corresponding stresses that are required
to be continuous. Thus, we have to prepare base functions in order to be able to adjust forstrain discontinuity at the interface. This requirement is satis"ed by preparing, at least, twotypes of base functions. I devised the following three base functions, u , u and u :H
H
Hu"xJyKzLe ,G
H

0Gr(r ,G

(6a)xJyKzLrI
G e , r GrGr ,
u"
G
G

H
rI

(6b)

xJyKzLrHu"e , 0Gr(r ,
G
G
H
rH
G

(7a)

xJyKzLrI
G e , r GrGr ,
u"
G
G

H
rI

(7b)

u"xJyKzLe ,
G
H

(8)

and
0GrGr ,


where displacement continuity at r"r is realized by introducing rI or rH into equations (6b)GGG
and (7).
Assuming that u is the necessary base function in this system, we have three possiblecombinations among them; u}u, u}u, and u}u}u. A preliminary assessment
was conducted in the spherical co-ordinate. Although the u}u}u combination seems
to be most e!ective, it was found to be too redundant to decrease internal independency
among them, or to be more likely to fail in solving eigenvalue problem at lower N. u}u
and u}u combinations usually provided a similar performance.
The con"guration of extremely low acoustic impedance is of prime interest for developing
the liquid}solid hybrid resonance experiments [6]. Therefore, I examined the case ofo /o "0)1, l /l "0)1, l /l "0)1 as an extreme case. As shown in Table 3, the u}u
G M
NG NM
QG QM
combination yields much better results than those by the u}u combination, for the caseJSV 20002963

7

THE xyzr ALGORITHM

TABLE 3Comparison of performance between the two base function sets for a laminated sphere;o /o "0)1, l /o "0)1, l /s "0)1, and r /r "0)5. The numbers after the mode indices are
G M
NG NM
QG QM
G M
the normalized analytical frequencies

Mode index
S (0)49493)
 ¹ (0)36680)
 ¹ (0)44477)
 ¹ (0)57878)
 u}u base functions
N"15, k"8, j"15u}u base functions
N"15, k"71)00001)00001)00001)00001)00001)00001)0001
1)0031TABLE 4
Convergent frequencies and their partial derivatives for several typical modes in a laminatedcylinder, whose height is equal to its diameter. Mode indices are an analogy of those of a sphere(Appendix B). The calculation conditions are N"15, k"11, j"17; N"17, k"11, j"21;
N"19, k"15, j"21; the values of k and j are not as signixcant as N
Mode index

Convergent
frequency, f* log f* log lNG

* log f* log lQG

* log f* log lNM

* log f* log l
QM

S
 
S
 
¹
 
¹
 

0)270630)340530)244310)44675

0)46
0)41
0
0

0)60
0)48
1)0
1)0

0)21e-30)34e-1000)86e-30)74e-10)52e-30)55e-3

with large wave number w owing to relatively low acoustic velocities in the inner region.
GL
Therefore, the u}u combination is the best recommended algorithm for a laminatedsystem in the spherical co-ordinate.
Thus, the performance of the u}u combination is examined as well in the cylindricalco-ordinate; the conditions except for shape are similar to those of Table 3. Although thegross convergency is worse than that in the spherical co-ordinate, some modes arecon"rmed to be convergent by N"15. Table 4 shows convergent numerical frequencies,
and its partial derivatives over velocities, for some typical modes free from angulardependency. Note that we have "ve digits of the same frequency over N"15, 17, 19 forthese modes. The convergent numerical frequencies of ¹ and ¹ are exactly the same as  those analytical solutions. Among these modes, ¹ converges as fast as by N"7, and 
S and ¹ by N"11. Besides the typical modes shown in Table 4, we have an additional  
12 modes matching the present convergent criterion over N"15, 17, 19. From themagnitude of partial derivatives, we can conclude that S modes are most sensitive to l ,O NGwhile ¹ modes are insensitive to l . These modes are really important to constrain l and
O 
NG
NG
l in the liquid}solid hybrid resonance technique, as reported in the preliminary work [6].
QG

ACKNOWLEDGMENTS

The author thanks Drs N. Suda, M. Kumazawa, H. Kumagai, and K. Ishihara, at
Nagoya University, for discussions. He is also grateful to Profs. I Suzuki, H. Oda, andJSV 20002963

8

A. YONEDA

M. Walter, at Okayama University, and Prof. I. Ohno, at Ehime University, for helpfulcomments.

REFERENCES
1. W. H. VISSCHER, A. MIGLIORI, T. M. BELL and R. A. REINERT 1991 Journal of the AcousticalSociety of America 90, 2154}2162. On the normal modes of free vibration of inhomogeneous andanisotropic elastic objects.
2. R. HOLLAND 1968 Journal of the Acoustical Society of America 43, 988}997. Resonant properties of
piezoelectric deramic rectangular parallelepipeds.
3. H. H. DEMAREST JR 1971 Journal of the Acoustical Society of America 49, 768}775. Cube-resonancemethod to determine the elastic constants of solids.
4. I. OHNO 1976 Journal of Physics of the Earth 24, 355}379. Free vibration of a rectangularparallelepiped crystal and its application to determination of elastic constants of orthorhombic
crystals.
5. M. SENOO and T. NISHIMURA 1984 Bulletin of Japanese Society of Mechanical Engineers 27,2339}2346. Measurement of elastic constants of polycrystals by the resonance method ina cylindrical specimen.
6. A. YONEDA 1998 Reviews of High Pressure in Science and ¹echnology 7, 79}82, Liquid}solid hybrid
resonance method for investigating acoustic properties of magmas.

APPENDIX A: BRIEF REVIEW OF THE xyz ALGORITHM DEVELOPEDBY VISSCHER E¹ A¸.

In the xyz algorithm, the displacement in the ith direction u (x, y, z) is expanded as
Gu" a u ,HG HG
Gu "xJ yK zL e ,HG
G(A1, A2)

where e (i"1, 2, 3) is a directional unit vector in the Cartesian co-ordinate. The integers l,
Gm, and n are limited by an integer N as
l#m#n)N.

(A3)

The Lagrangian ¸ is written as



¸"



1
*u *u
1
G I d<,ouu u ! C
G G 2 GHIJ *x *x
2
H J

(A4)

where o is the density and u the angular frequency. Substituting equation (A1) into
equation (A4), we have¸" (ua2Ea!a2Ca),(A5)where the (j , j ) component of the matrices E and C are given as
G I

E "d ou u d<,
HGHI
GI
HG HI

*u *u
G I d<.
C " C
HG HI
GHIJ *x *x
H J

(A6, A7)

As we express the displacement as equation (A1), the variation of displacement is given by
aPa#da.
JSV 20002963

(A8)

9

THE xyzr ALGORITHM

Thus, the variation of ¸ isd¸" (a2#da2) (uE!C) (a#da)! a2(uE!C)a" a2 (uE!C)da# da2 (uE!C)a.



(A9)

Note that both E and C are asymmetric matrices, the stationary condition of d¸"0 issatis"ed withuEa"Ca.

(A10)As E is not a diagonal matrix owing to the non-orthogonality of the present base functions,or power functions, equation (A10) is a generalized eigenvalue problem. Nowadays, we
can solve equation (A10) on an ordinary workstation by means of a ready-madesubroutine found in various numerical computation packages.
Note that the Legendre base functions provide us with a diagonal matrix of E owing totheir orthogonality. Thus equation (A10) is simpli"ed as an ordinary eigenvalue problem;
it used to be critical to save memory and CPU time even on the most powerful computer20}30 years ago.

APPENDIX B: GENERALIZED SOLUTIONS OF EIGENVIBRATIONFOR A SPHERE AND A CYLINDER

Sphere:Note that ; corresponds to P-waves, and < and = are S-waves. The potential for P-waves
is written asJ(hr) uU"e SRPK (cos h) L>
eK .
L
(r
For S-waves, h should be replaced by k. Therefore, angular frequency is given asu"l h"l k for h and k satisfying the boundary conditions
N
Q
*U
;" ,
P *r

*U
;" ,
F r*h

*U
;"
,
( r sin h*

< "0,
P

*U
<"
,
F sin h**U< "! ,
(
*h

 





*U
1 *U
1
*
= "!
sin h
#
,
P
r sin h *h*hsin h * *(rU)=",
F r*r*h;"P*(rU)=",
r sin h*r *hnJ(hr)!hrJ
(hr)
>L
>L
PK (cos h) e K(e SR,Lr
JSV 20002963

10

A. YONEDA

J(hr) +(n#1!m)PK (cos h)!(n#1) cos hPK (cos h),
L>
L; " >Le K(e SR ,
F
sin h
r
; "i
(

mJ

>L
r

(hr) PK (cos h)Le K( e SR ,
sin h

< "0,
P
< "i
F

mJ

>L
r

(kr) PK (cos h)Le K( e SR,sin hJ(kr) +(n#1!m)PK (cos h)!(n#1) cos hPK (cos h),
L>
L< " >Le K(e SR ,
(
r
sin h
n(n#1) J
(kr)
>L
="
PK (cos h)e K(e SR,PLr+(n#1)J(kr)!krJ
(kr),
>L
>L
="
F
r
+(n#1!m)PK (cos h)!(n#1) cos hPK (cos h),L>
L;e K(e SR ,
sin h
= "im
(



(n#1)J>L(kr)!krJ(kr) PK (cos h)>LLe K( e SR .
r
sin h

We can classify two types of vibration modes for an isotropic sphere as seismologicalnotation for the earth's free oscillation. <s have no radial component, and are calledtoroidal modes ¹ . The combinations of ;s and =s have the radial component, and are
O L
called spheroidal modes S . The former subscript q is radial order corresponding to the
O L
number of nodal surfaces inside the sphere, while the latter angular order corresponds tothe number of nodal lines on the spherical surface. Note that the number of nodal lines on
the surface is n for S modes and n!1 for ¹ modes.
O L
O L
Note that there are two types of Bessel functions; the "rst kind Bessel function J(r) is "niteat r"0, while the second kind of Y(r) is in"nite at r"0. Thus the above expression,
written with the "rst kind of Bessel function, J, is only valid for a solid sphere, which
comprises r"0 in the system. On the other hand, the second kind Bessel function, Y,
must be introduced for a bored sphere or a shell, because it is free from the in"nity at
r"0. Thus, J should be replaced with AJ #BY , such as
L
L
L
; "i
(

m(AJ
(hr)#BY(hr)) PK (cos h)
>L
>L
L
e K( e SR ,
r
sin h

where A and B are constants adjusted to satisfy the inner and the outer boundary
conditions. The situation is the same in the cylindrical co-ordinate as well.
JSV 20002963

11

THE xyzr ALGORITHM

Cylinder:The notations are similar as in the case of the spherical co-ordinate. However, the potential
is given as
U"e SR J (hr)e LXe KF.
KThe angular frequency is given as u"l (h#n"l (k#n for h, k and n satisfying
QN
the boundary conditions. Although the following solutions are written for J, it should be
replaced with AJ#BY in the case of a tube, as discussed in the spherical co-ordinate:
*U
;" ,
P *r

*U
;" ,
F r*h

*U
;" ,
X *z

*U
,
<"
P r*h

*U
< "! ,
F
*r

< "0,
X

*U
="
,
P *r*z*U=",
F r*h*z



;"
P



mJ (hr)!hrJ
(hr)
K
K>
e KFe LXe SR,r; "nJ (hr) e KF(ie LX)e SR,
X
K








*U *U
*U
= "!
# #
,
X
*r r*r r*hmJ (hr);" K(ie KF) e LXe SR,FrmJ (kr)(ie KF) e LXe SR,
<" K
P
r




< "!
F

mJ (kr)!krJ
(kr)
K
K>
e KFe LXe SR,
r

< "0,
X

= "n
P

mJ (kr)!krJ
(kr)
K
K>
e KFe LXe SR,rmnJ (kr)K="(ie KF) e LXe SR,Fr= "kJ (kr) e KF (ie LX) e SR.
X
K

JSV 20002963