Earth Planets Space, 54, 763–770, 2002

Intrinsic eigenvibration frequency in the resonant ultrasound spectroscopy:
Evidence for a coupling vibration between a sample and transducers
Akira Yoneda
Institute for Study of the Earth’s Interior, Okayama University, Misasa, Tottori 682-0193, Japan(Received October 22, 2001; Revised June 4, 2002; Accepted June 6, 2002)

The resonant ultrasound spectroscopy enables us to measure elastic constants of various materials with highaccuracy. One of its curious features is resonance frequency shift caused by modulation of clamp force for a sample,
although its mechanism has not been clarified yet. A coupling vibration model is newly proposed to interpret thisphenomenon and predict its functional form. Further, it is confirmed that the resulted functional forms dependon characteristics of each eigenvibration mode. The extrapolation based on those functional forms enables us tocorrect measured frequency to the intrinsic one with improved accuracy at least by one order of magnitude. Thus
the correction is essential in order to obtain accurate derivatives of frequency with respect to temperature, pressure,composition, and so forth.

1.

Introduction

The resonance method, or resonant ultrasound spectroscopy (RUS), has been extensively used to elucidate elastic properties of various materials (e.g., Sumino et al., 1976;
Migliori and Sarrao, 1997; Ohno et al., 2000; Suzuki et al.,2000). One of the main features of the current RUS is thatsample must contact to a couple of transducers; one is an oscillator, and the other a receiver for sample vibration (Fig. 1).
It has been recognized that an observed resonance frequency f is slightly shifted with magnitude of clamp-force
F to the sample, although its mechanism has not been clarified yet. Therefore intrinsic resonance frequencies, f 0 , havebeen determined through an empirical extrapolation of f towards F = 0 (e.g., Sumino et al., 1976), because the current
RUS requires a finite magnitude of F in order to transmit
vibration energy from transducer to sample, and vice versa.
A coupling vibration between a sample and transducers isinevitably excited on sample-transducer interface. I foundthat the coupling vibration enables us to explain the magnitude of observed frequency shifts. Then I confirmed that thepredicted functional forms explain experimental frequency
shifts observed in a sphere sample of soda glass. In thisstudy, I devised a clamp force controller (Appendix A) to
increase efficiency and reliability.
I classify vibrational modes of an isotropic sphere into
two types according to seismological notation for the earth’s
free oscillation. While spheroidal modes, q Sn , are of radial
motion, troidal modes, q Tn , are of no radial motion. The
former subscript of q is radial order corresponding numberof nodal surface inside the sphere, while the latter angularorder corresponding number of nodal lines on the sphericalsurface. Note that the number of nodal lines on the surfacec The Society of Geomagnetism and Earth, Planetary and Space Sciences
Copy right
(SGEPSS); The Seismological Society of Japan; The Volcanological Society of Japan;
The Geodetic Society of Japan; The Japanese Society for Planetary Sciences.

is n for q Sn modes and n − 1 for q Tn modes.
Note that SI unit is used in this literature except for showing magnitude of clamp force measured in ‘gw’. As 1 gw =0.0098 N  0.01 N, it is not difficult to convert ‘gw’ to ‘N’,or vice versa, without losing perspective idea.

2.

Coupling Vibration Model

We have experienced that the magnitude of  f / f 0 is of
the order of 10−3 against a few gw of F (see Fig. 2 as
an example). Neither deformation of sample nor increaseof elastic wave velocities is able to explain such a large
deviation of resonance frequency (Appendix B).
Only coupling vibration model (Fig. 3) is able to explain
the observed magnitude of frequency shift. The intrinsic angular frequency ω0 (= 2π f 0 ) is defined for the free oscillation state characterized by ks = 0, or, F = 0. Thus we
have
2k
ω02 =
(1)
m
for the decoupling free, state, where m and k are mass andspring constant, respectively, of the simplified system. Thenthe effect of weak spring ks is incorporated into the system,
we have
2k + ksks
ω2 == ω02 +
(2)mm
for a coupling state characterized by ks , which must be afunction of F through an elastic contact between a sample
and transducers. This is the essence of the coupling vibration
model.
The coupling spring ks is a model of the elastic property
at around the contacting region between a sample and transducers; I apply the contacting sphere theory to describe theproperty (see Fig. 4 and Appendix B).
Although actual eigenvibration modes are not identicalwith that simplified configuration shown by Fig. 3, we can

763

764

A. YONEDA: INTRINSIC EIGENVIBRATION FREQUENCY IN RUS

still expect the similar relation under the reasonable assumpx = −g(ν, Y, R)F 2/3(3)
tions (Appendix C).
According to the contacting sphere theory, the distance where the factor g is a function of Poisson ratio ν, Youngbetween two spheres, x, varies with F as
modulus Y , and curvature radius R. The area of contact, A,
is similarly expressed as
A = h(ν, Y, R)F 2/3 .

(4)

As long as we do not disturb the sample-transducer couplingstate, we can assume constant values for g and h during
modulation of F.
First let consider the case of 0 S0 mode characterized as
pure normal displacement of sample surface. In this case,the modeled spring constant, ks , is given as
ks = −dF3= g −1 F 1/3 .
dx
2

(5)

Thus we have the relation:
ks ∝ F 1/3 .

(6)

Combining (2) and (6), we have a result:ω2 = ω2 − ω02 ∝ F 1/3 ,or  f 2 = f 2 − f 02 ∝ F 1/3 .
(7)

The characteristics of T modes are completely different
withthose of S modes; T modes are free from normal surface
Fig. 1. Block diagram for a RUS experiment. Thick lines correspond tocoaxial cable for sinusoidal wave transmission. Thin lines are GPIB cable displacement. Therefore, some modification is required forfor digitized data. A set of longitudinal wave transducers (M121 from T modes, because the weak spring constant ks corresponds
Panametrics Co.) and shear wave transducers (V156 from Panametrics
to tangential displacement between sample and transducer. It
Co.) are used. Transformer is used to cut low frequency noise.

Fig. 2. An example of frequency shift caused by clamp-force modulation; they are 0 T2 peaks of a soda glass sample. The clamp force values are shown ineach subplot. We can see that even a few gw clamp force variation causes ∼1 kHz peak shift.

A. YONEDA: INTRINSIC EIGENVIBRATION FREQUENCY IN RUS

3.

765

Experimental

A soda glass sphere (diameter, d: 2.890 mm; mass, m:
3.18 × 10−5 kg; density, ρ: 2.52 × 103 kg/m3 ) is used toevaluate the above model, because it is homogeneous andelastically isotropic. In the resonance frequency measurements, clamp force for sample is automatically controlled by
a combination of an electric actuator, a load cell, and a voltage regulator (Appendix A). The sample temperature was
kept at 25 ± 0.1◦ C by placing the main part of the system ina temperature-controlled chamber in order to remove a bias
caused by temperature fluctuation.
The resonance frequencies measured at various clamp
forces are curve-fitted by assuming the proposed functional
forms; note that the present correction method does not require any additional information such as physical propertiesof sample, etc.

4.
Fig. 3. Illustration for a coupling vibration model. In this figure, an eigen
vibration mode is simplified to a harmonic oscillator; the stronger spring
constant k corresponds to intrinsic elasticity of sample, and the weaker
one ks to extrinsic property at around contact region between sample andtransducer. Note that total sample mass is 2m in this configuration.

Fig. 4. Illustration for the contacting sphere theory. R is the radius of eachsphere. The distance x between the center of the two spheres decreases
with F. Note that a contacting circle of diameter 2a is formed at theinterface of the two sphere owing to their elastic deformation. Therefore,
the area of contacting circle, A, is πa 2 .

may be reasonable to assume ks proportional to the contact
area, A, between sample and transducer. The contacting
sphere theory shows A is in proportion with F 2/3 . Therefore
we may expectω2 ∝ ks ∝ F 2/3(8)

for T modes.
The S modes except 0 S0 are constituted with both normal
and tangential surface displacements. Therefore it is difficult
to predict the functional form theoretically. It should be
revisited after experimental observation.

Results and Discussion

0 S0 mode

Figure 5 shows the observed resonance frequency variations observed in two runs separately conducted; we can recognize ∼0.5 kHz shift of resonance frequency even at a few
gw of F. The resulted trends are consistent with the relation(7). Please note that the each trend in Fig. 5 was measured ina single run keeping contacting condition; otherwise, we canrecognize that  f may differ even at the same F.
The extrapolated free oscillation frequencies converge to1611.3–1611.4 kHz. In the following discussion, the intrinsic eigenfrequency of 0 S0 is fixed at 1611.4 kHz. As a sample
is clamped vertically by a couple of transducers, the sample
weight may be a possible source of error in f 0 . However its
small weight is negligible, because it causes only ∼0.02 kHz
deviation.
As shown in Appendices B and C, numerical estimation
gives theoretical value of 2.2 × 1011 Hz2 /N1/3 for ω2 /F 1/3 .
The observed gradients in Fig. 5 are calculated to be 2.9 ×
1011 Hz2 /N1/3 and 1.8 × 1011 Hz2 /N1/3 , respectively, foreach run. This agreement in magnitude between theory andexperiments is additional support for the coupling vibration
model.
T modesResonance frequency variations in T modes (0 T4 , 1 T1 ) are
shown in Fig. 6; those two modes are selected for comparison, because they have similar resonance frequencies in spiteof quite different vibrational character. The main difference
is in radial order q; q = 0 for 0 T4 and q = 1 for 1 T1 . Although the amplitude information is not presented in Fig. 6,
the signal amplitude of 0 T4 is much larger by ∼20 times thanthat of 1 T1 .
It is impossible to judge which is better model, the 1/3
model or the 2/3 model, just from a first look of Fig. 6.
However we can recognize that the slope of f 2 variationfor 1 T1 is much smaller than that for 0 T4 . Note that 1 T1 ischaracteristic mode having a nodal surface inside, and thus
relatively small surface motion. The extrapolated value for1 T1 is 2190.2 kHz for the 1/3 model and 2190.3 kHz for the
2/3 model; the difference is as small as 0.1 kHz.
Note that every T mode frequency is determined only
by sample diameter, d, and shear wave velocity, vs (e.g.,
Schreiber et al., 1973). Therefore the resonance frequency of766

A. YONEDA: INTRINSIC EIGENVIBRATION FREQUENCY IN RUS

12

2.599

x 10

f 2, Hz2

2.59852.598
2.59752.597
2.5965
2.596
0

0.2

0.4

0.6

0.8
F1/3 , gw1/3

1

1.2

1.4

1.6

0.5

1

1.5

2
F, gw

2.5

3

3.5

4

(a)
6

1.6122

x 10

f, Hz

1.612
1.61181.6116
1.6114
1.61120(b)Fig. 5. Frequency shifts for 0 S0 mode in two measurements. Upper frame: f 2 − F 1/3 plot; Lower frame: f − F plot.

12

4.81

12

1T1

x 10

3.755
3.75
f 2, Hz2

f 2, Hz2

4.805
4.8
4.795

3.745
3.74

4.79
0

0.5

1
1.5
F 1/3 , gw1/3

(a1)

3.735
0

2

0.5
(b1)

12

4.805

0T4

x 10

1
1.5
F 1/3 , gw1/3

2

2
F 2/3 , gw 2/3

3

12

x 10

3.756

x 10

3.752

4.8

f 2, Hz2

f 2, Hz2

3.754

4.795

3.75
3.748
3.746
3.744

4.79
0

1
(a2)

2
F 2/3 , gw 2/3

3

3.742
0

1
(b2)

Fig. 6. Frequency shifts for 1 T1 and 0 T4 . (a1) f 2 − F 1/3 plot for 1 T1 ; (a2) f 2 − F 2/3 plot for 1 T1 ; (b1) f 2 − F 1/3 plot for 0 T4 ; (b2) f 2 − F 2/3 plot for
0 T4 .

A. YONEDA: INTRINSIC EIGENVIBRATION FREQUENCY IN RUS

767

Table 1. The results for corrected frequencies of T modes. Note that normalized ratio is based on f 0 of 1 T1 in each column.

Modes

f 0 byf 0 byNormalized ratio in

Normalized ratio in

Theoreticalthe 1/3 modelthe 2/3 modelthe 1/3 modelthe 2/3 model

ratio

0 T2

948.9 kHz950.3 kHz

0.43330.43390.43400 T31466.91468.20.67000.67030.6706

0 T4

1932.91934.80.88250.88340.88401 T12190.22190.3

1.0

1.0

1.0

Table 2. The results for corrected frequencies of selected S modes. The theoretically calculated frequencies are also shown in the last column forcomparison.

Modes

f 0 byf 0 byAverage of the 1/3

f 0 by theoreticalthe 1/3 modelthe 2/3 modeland the 2/3 models

calculation

0 S21001.7 kHz

1003.0 kHz

1002.4 kHz

1002.1 kHz

1 S1

1282.31282.81282.6

1282.8

1 S2

1823.2

1824.1

1823.7

1824.30 S31483.51485.61484.61484.5

1 T1 can constrain vs most reliably, because its frequency shift

5.

is much smaller than that of other T modes. Table 1 summarizes f 0 ratio with comparing the theoretical one. We can distinguish which is better model from the f 0 ratio among the Tmodes. In every case shown in Table 1, the 2/3 model gives
better internal consistency; differences with the theoreticalratios are only a few parts of 10−4 . Therefore I concluded
that the 2/3 model is the right model for T modes. Thus theintrinsic eigenfrequency of 1 T1 is fixed at 2190.3 kHz.
Now we can obtain longitudinal wave velocity v p =5811.7 m/s, shear velocity vs = 3450.4 m/s, Poisson ratio v = 0.228, and Young modulus Y = 73.7 GPa from
f 0 = 1611.4 kHz for 0 S0 , f 0 = 2190.3 kHz for 1 T1 ,d = 2.890 mm, and ρ = 2.52 × 103 kg/m3 ; these values
are used in the following discussion and Appendix B.
S modes except 0 S0The S modes except 0 S0 are a combination of normal
and tangential surface displacement. Therefore it is quite
probable that those f − F trends are composite of the 1/3
and the 2/3 models. Table 2 compares experimental resultswith calculated frequencies based on v p and vs determinedin this work.
Comparing those experimental values with calculatedones, we can recognize that the average of the 1/3 and the
2/3 models is consistent with the calculated frequencies. The
differences are of only a few parts of 10−4 . Therefore I propose the averaged value for general S modes as a practicaland reasonable solution.
The above averaging method may be a kind of compromise; in order to develop more sophisticated correction
method, we have to know vibration pattern at the interface between a sample and transducers. Unfortunately it isquite difficult to generate reproducible vibration pattern inthe present experimental technique.

The coupling vibration model clearly explains the frequency shift caused by clamp force modulation. The phenomenon is quite similar to the bond effect recognized in ultrasound velocity measurement (e.g., Yoneda, 1990; Spetzler
et al., 1993), because both the problems are originated fromvibration energy transfer.
If you just want to obtain elastic constants with accuracyof 10−3 , the resonance frequency correction may be unnecessary. However if you need better accuracy, for instanceto obtain precise derivatives of elastic constants on temperature, pressure, or composition, etc., it is very important toapply the present correction method in order to eliminate theambiguity in resonance frequency.
Last of all, I would like to stress the following three notes
in order to enhance accuracy and efficiency:

Concluding Remarks(1) 0 S0 and T modes are special modes realized only forelastically isotropic sphere sample. If the sample shapeis not sphere or physical property of sample in not
isotropic, eigenviblation modes are basically combination of normal and tangential surface displacement.
Therefore the average of 1/3 and 2/3 models should berecommended in general in order to conduct resonancefrequency correction.
(2) Sample and transducers should be clean enough. If thesurface is contaminated by greasy substance, you may
observe substantial frequency shifts even at negative
clamp force. It is not surprising, because “grease” cankeep substantial contact between a sample and transducers against tensional force.
(3) As far as my experience, larger amplitude of eigenfrequency peaks tends to have larger frequency shift; an
example is amplitude contrast between 1 T1 and 0 T4 . It

768

A. YONEDA: INTRINSIC EIGENVIBRATION FREQUENCY IN RUS

is reasonable, because harder coupling between sample
and transducer may give simultaneously larger amplitude and larger frequency shift. Therefore smaller amplitude modes may be more important in order to constrain sample properties in spite of relative difficulty indetermining the resonant frequency.
Acknowledgments. The author is indebted to Mr. Maeda for designing and assembling the electronic circuit for clamp-force controller. He also thanks Prof. I. Suzuki at Okayama Univ. and Prof.
I. Ohno at Ehime Univ. for helpful comments to the present model.
This work was partially supported by a Grant-in Aid for Scientific
Research No. 12440121 from the Japan Society for the Promotionof Science.

Appendix A. Automatic Clamp-Force Controller

Appendix B. Contacting Sphere Theory and its
Characteristics
The contacting sphere theory itself is found in many text
book (e.g., Timoshenko and Goodier, 1970). It is wellknown that the theory is applicable as well for sphere-plane
contact by assuming infinitely large curvature radius to theplane. Note that sample and transducer are assigned by subsuffixes, 1 and 2, respectively. Elastic properties of the sample are ν1 = 0.228 and Y1 = 73.8 GPa. Its curvature radius,
R1 , is identical with its radius r0 ; thus R1 = r0 = d/2 =1.445 mm. Note that the transducers’ surface is covered bya thin tungsten-carbide plate (ν2 = 0.22; Y2 = 534 GPa;
R2 = ∞).
By substituting those values, the contacting radius a isgiven as
3 3π(k 1 + k 2 )R1
a=(B.1)
F 1/3 ≈ 2.5 × 10−5 F 1/3
4

So far the clamp force has been controlled by manualadjusting of weight on a balance equipped in the resonancemeasurement system. This is not only inconvenient but also
where F is force, and k1 and k2 are defined as
inaccurate due to frictional bias in a balance. Note that theuncertainty of an ordinary commercial balance is as large as
1 − νi2
ki =,
i = 1, 2.
(B.2)
0.2 gw, which may cause trouble in the present frequency
π Yicorrection.
Therefore I devised an automatic clamp force controller Then d x, the normal displacement between sample andfree from frictional error. Figure A1 shows the concept of transducer, is given as
the instrument. Figure A2 shows the stability of F, or output
2
2
of load-cell. We can see that F is controlled automatically
3 9π (k 1 + k 2 )dx =F 2/3 ≈ 4.4 × 10−7 F 2/3 .
(B.3)
within 0.1 gw, or ∼0.001 N, both for long-term variation and
16R1
short period fluctuation.
Then we have ks as∂F
3F 1/3ks ==≈ 3.4 × 106 F 1/3 . (B.4)∂x
2 × (4.4 × 10−7 )In Fig. B1, I present variations of a, A (= πa 2 ), d x, and ks
with F in the present experimental configuration. These parameters are useful to evaluate effect of other possible mechanisms for resonance frequency shift. First, length change
x/x is estimated to be ∼10−5 at 5 gw for a soda glass sample of ∼3 mm diameter. Second, mean pressure at contact
area is estimated to be 106 Pa at 5 gw as well. We can assume that v/v ∼ p/C, where v is elastic wave velocity,
p pressure increment, and C a representative elastic constant ∼100 GPa. Then substitute p = 106 Pa, we havev/v ∼ 10−5 . Note that relative resonance frequency shift f / f is given by x/x + v/v. Therefore even the combination of the two mechanisms cannot explain the observed
frequency shift of order of  f / f ∼ 10−3 .

Appendix C. Coupling Vibration for a General
Eigenvibration Mode

Fig. A1. Schematic diagram of the clamp-force controller. The voice coil isa clamp-force stabilizing actuator; it is operated by a feedback DC current
supplier minimizing the difference between the input target voltage andthe load cell output. The spring is used to minimize weight bias causedby coaxial cable for transducer.

Actual eigen vibration of sample is much more complicated than the simplified model shown in Fig. 3. Therefore I
show here that the coupling vibrational model is still reasonable to any eigenvibration mode. Let suppose the following
conditions for eigenvibration modes and the coupling vibration.
(I) The modeled spring ks yields only small perturbation
for an eigenvibration mode.
(II) It does not have its own mass like an ideal spring, orcannot store any kinetic energy.

A. YONEDA: INTRINSIC EIGENVIBRATION FREQUENCY IN RUS-6769

sodaglass4

x 10
10

amplitude, V

9
8
7
6
5

0

50

100

150

200

250

300

50

100

150
time, min

200

250

300

3.2
3

load, gw

2.8
2.6
2.4
2.2
2
1.8
1.6
0

Fig. A2. An example of resonance signals by using the clamp-force controller. We can see that the clamp force is controlled within 0.1 gw. Upper frame:
Resonance signals. Lower frame: Load, or clamp force F, during the measurement. We can see that amplitude of resonance peaks decreases withdecreasing F.

6

0.008

5
dx, mm

a, mm

-5

0.01

0.006
0.004

x 10

4
3
2

0.002

1

0
0

1

2

(a)

3
F, gw

4

0
0

5

1

2

3
F, gw

4

5

1

2

3
F, gw

4

5

(c)

-4

3

x 10

1400
1200
1000

2

ks, N/ mm

A, mm2

2.5

1.5
1

800
600
400

0.5

200

0
0

1
(b)

2

3
F, gw

4

5

0
0
(d)

Fig. B1. Useful parameters in the coupling model as a function of clamp force F. (a) radius of contacting circle, a; (b) area of contacting circle, A; (c)displacement of sphere, d x; (d) spring constant, ks .

770

A. YONEDA: INTRINSIC EIGENVIBRATION FREQUENCY IN RUS

According to the eigenvibration analysis (e.g., Ohno, 1.445 mm) into the following equation
1976; Visscher et al., 1991; Yoneda, 2000) intrinsic angu
1 2
2
lar frequency must satisfy
ρu r d V = ks δ(u)2 = ks u r (r0 )2
(C.6)
ω
2


1
1
E k = ω02
ρu i u i d V =Ci jkl εi j εkl d V = E S (C.1)
we have ω2 /ks = 6.5 × 104 . Combining (B.4) and (C.6),
2
2
we have
where E K and E S are kinetic energy and strain energy, respectively, ρ density, u displacement, C elastic stiffness conω2 ≈ 2.2 × 1011 F 1/3 .
(C.7)
stants, ε strain. The assumption (I) allows us to modify (C.1)
Note that the above estimation is only possible for 0 S0 , beascause it is the unique mode of uniform surface displacement,
111
ω2
ρu i u i d V =Ci jkl εi j εkl d V + ks1 δ(u 1 )2
u r (r0 ).
2
2
2
1+ ks2 δ(u 2 )2(C.2) References2
Migliori, A. and J. Sarrao, Resonant ultrasound spectroscopy, John Wiley &Sons, New York, 1997.

where δ(u) is an effective displacement for spring ks , and Ohno, I., Free vibration of a rectangular parallelepiped crystal and its applisub-suffixes, 1 and 2, specifies contact states against the
cation to determination of elastic constants of orthorhombic crystals, J.
Phys. Earth, 24, 355–379, 1976.
two transducers. Because δ(u 1 )/δ(u 2 ) is considered to bea constant in each eigenvibration mode, we can combine the Ohno, I., M. Abe, M. Kimura, Y. Hanayama, H. Oda, and I. Suzuki, Elasticity measurement of silica glass under gas pressure, American Mineraladditional terms in (C.2) asogist, 85, 288–291, 2000.
11
ks1 δ(u 1 )2 + ks2 δ(u 2 )2 = ks δ(u)2 .
2
2

(C.3)

Note that there is no additional term for kinetic energy because of assumption (II). Combining (C.1), (C.2), and (C.3),
we can deriveω2 = ω2 − ω02 ∝ ks

(C.4)

for any eigenvibration mode. Therefore we can apply the
coupling vibration model in general.
Last of all, let estimate ω2 /F 1/3 for 0 S0 . In the 0 S0mode, the surface displacements u is simply expressed in thespherical coordinates:
ur = −

h J1.5 (hr ),√ruθ = uφ = 0(C.5)

where h (= ω/v p ) is wave number. Then substituting (C.5)
and sample’s parameters (ρ = 2.52×103 kg/m3 ; r0 = d/2 =Schreiber, E., O. L. Anderson, and N. Soga, Elastic constants and theirmeasurement, McGraw Hill, New York, 1973.
Spetzler, H. A., G. Chen, S. Whitehead, and I. C. Getting, A new ultrasonicinterferometry for the determination of equation of state parameter ofsub-milimeter single crystal, Pure Appl. Geophys., 141, 341–377, 1993.
Sumino, Y., I. Ohno, T. Goto, and M. Kumazawa, Measurement of elastic
constants and internal frictions on single-crystal MgO by rectangularparallelepiped resonance, J. Phys. Earth, 24, 263–273, 1976.
Suzuki, I., I. Ohno, and O. L. Anderson, Harmonic and anharmonic properties of spinel MgAl2 O4 , American Mineralogist, 85, 304–311, 2000.
Timoshenko, S. P. and J. N. Goodier, Theory of Elasticity third edition,
p. 412, McGraw Hill, New York, 1970.
Visscher, W. M., A. Migliori, T. M. Bell, and R. A. Reinert, On the normal
modes of free vibration of inhomogeneous and anisotropic elastic objects,
J. Acoust. Soc. Am., 90, 2154–2162, 1991.
Yoneda, A., Pressure derivatives of elastic constants of single crystal MgOand MgAl2 O4 , J. Phys. Earth, 38, 19–55, 1990.
Yoneda, A., The xyzr algorithm specialized for eigenvibration problem ofbored and laminated objects, J. Sound and Vibration, 236, 431–441,
2000.

A. Yoneda (e-mail: yoneda@misasa.okayama-u.ac.jp)