Resonance Phenomena in Goupillaud-type Media

by George A. Gazonas and Ani P. Velo

ARL-TR-5384 October 2010

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the Army position unless so designated

by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
Army Research Laboratory

Aberdeen Proving Ground, MD 21005-5066

ARL-TR-5384 October 2010

Resonance Phenomena in Goupillaud-type Media

George A. Gazonas

Weapons and Materials Research Directorate, ARL

Ani P. Velo

Department of Mathematics and Computer Science, University of San Diego

Approved for public release; distribution is unlimited.
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering 

and maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to 

comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To) 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

5c. PROGRAM ELEMENT NUMBER 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

6. AUTHOR(S) 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S) 

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

15. SUBJECT TERMS

16. SECURITY CLASSIFICATION OF: 

19a. NAME OF RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 

17. LIMITATION 

OF ABSTRACT 

18. NUMBER 

OF PAGES 

19b. TELEPHONE NUMBER (Include area code)

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18

October 2010 Final

ARL-TR-5384

Approved for public release; distribution is unlimited.

May 2009-May 2010

Resonance Phenomena in Goupillaud-type Media

George A. Gazonas AH84

Ani P. Velo

U.S. Army Research Laboratory

ATTN: RDRL-WMM-B

Aberdeen Proving Ground, MD 21005-5069

primary author’s email: <george.aristotle.gazonas@us.army.mil>

The discrete resonance frequency spectrum is derived for an m-layered Goupillaud-type elastic medium subjected to a discrete,

time-harmonic forcing function at one end, with the other end fixed. Analytical stress solutions are derived from a global

system of recursion relationships using z-transform methods, where the determinant of the resulting global system matrix |Am|

in the z-space is a palindromic polynomial with real coefficients. The zeros of the palindromic polynomial are distinct and are

proven to lie on the unit circle for 1 ≤ m ≤ 5. An important result is the physical interpretation that all positive angles,

coterminal with the angles corresponding to the zeros of |Am| on the unit circle, represent the resonance frequency spectrum for

the discrete model. The resonance frequency results are then extended to analytically describe the natural frequency spectrum

of a free-fixed m-layered Goupillaud-type medium. The predicted natural frequency spectrum is validated by independently

solving a simplified version of the frequency equation involving a transformation of the spatial variable. The natural frequency

spectrum is shown to depend only on the layer impedance ratios and is inversely proportional to the equal wave travel time for

each layer. A sequence of resonance frequencies for the discrete model is found, which is universal for all multilayered designs

with an odd number of layers, and independent of any design parameters.

wave propagation; z-transform methods; palindromic polynomials with real coefficients; trigonometric sums

UU 76

George A. Gazonas

Unclassified Unclassified Unclassified 410-306-0863

ii
Contents

List of Figures vi

List of Tables viii

Acknowledgments ix

1. Introduction 1

2. Forced Vibrations and Resonance Frequencies: The Discrete Model 3

2.1 Stress Solutions at Non-resonance Frequencies .......................................... 12

2.2 Stress Solutions at Resonance Frequencies ............................................... 13

2.2.1 Case I: cos  ̃ω = −1 = cos θ0 (for m odd only) .......................... 13

2.2.2 Case II: cos  ̃ω = cos θk∗ ........................................................ 15

2.2.3 Case II.a.: sin  ̃ω = sin θk∗ ..................................................... 15

2.2.4 Case II.b.: sin  ̃ω = − sin θk∗ = sin(−θk∗ ) ................................. 16

2.3 The Two-layer Case (m = 2) ........................................................... 17

2.3.1 Stress Solutions at Non-resonance Frequencies ................................. 17

2.3.2 Stress Solutions at Resonance Frequencies ...................................... 18

2.4 The Three-layer Case (m = 3) .......................................................... 21

2.4.1 Stress Solutions at Non-resonance Frequencies ................................. 22

2.4.2 Stress Solutions at Resonance Frequencies ...................................... 24

iii
3. Analytical Formulas for the Natural Frequency Spectrum of a Free-fixed Goupillaudtype Layered Elastic Strip 27

3.1 The Frequency Equation: Explicit Formula for the Natural Frequency Spectrum..... 27

3.2 Heuristic Approach: Implicit General Formula for the Natural Frequency Spectrum . 30

3.3 Interpretation and Illustration of the Natural Frequency Results ........................ 31

3.3.1 Case a. Two-layered Goupillaud-type Elastic Strip with Unequal Layer

Lengths ............................................................................ 31

3.3.2 Case b. Simplified Two-layered Goupillaud-type Elastic Strip with Equal

Layer Lengths ..................................................................... 32

3.3.3 The Predicted Natural Frequency Spectrum ..................................... 32

3.3.4 Case c. Design Modification That Gives a Desired Frequency Spectrum

Within Limitations ................................................................ 33

4. Other Applications of the Frequency Results 34

4.1 One-dimensional Layered Media with a Common Frequency Spectrum ............... 34

4.2 Resonance Frequencies and Optimal Designs ............................................ 35

4.3 Natural Frequencies of a Free-fixed Non-Goupillaud-type Layered Strip with Integer Layer Length Ratios ................................................................... 36

4.4 Representation of the Stress Solutions in Terms of the Chebyshev Polynomials ....... 37

5. Conclusions 38

6. References 41

Appendix A. Stress Formulas for a Heaviside Loading 45

Appendix B. Zeros of the Determinant |Am| 47

iv
Appendix C. Definition of the Matrix Dm 49

Appendix D. Necessary Relations Among the Base Angles and Design Parameters 51

Appendix E. Optimal Base Angles 53

Appendix F. Solutions of the Frequency Equations for a Free-fixed m-layered Goupillaudtype Elastic Strip (2 ≤ m ≤ 5) 55

Appendix G. Material Properties 57

Distribution List 59

v
List of Figures

Figure 1. Lagrangian diagram for an elastic strip made of m-layers of equal wave travel time. 5

Figure 2. Stress time history for a layered strip with impedance ratios α1 = 3, α2 = 2,

α3 = 1.5, α4 = 2.2, α5 = 0.3, α6 = 1.7, α7 = 1.4, α8 = 3.1, α9 = 0.8, and α10 = 4

as applicable; τ = 1 and loading f(n) = cos(nω ̃) applied at ξ = 0 for ω ̃ = π. (a) The

middle of the second layer of a 7-layered strip located at ξ = 3/14. (b) The middle of

the second layer of an 11-layered strip located at ξ = 3/22. .................................. 15

Figure 3. Stress time history for a two-layered Goupillaud-type strip in the middle of the

second layer located at ξ = 3/4. Impedance ratio α1 = 1/3, angle θ1 = 2π/3, τ = 1,

and loading f(n) = sin(nω ̃) at ξ = 0. (a) ω ̃ = π/4. (b) ω ̃ = 0.9 · θ1 = 0.6π. ............. 18

Figure 4. Stress time history for a two-layered Goupillaud-type strip with impedance ratio

α1 = 1/3, τ = 1, subjected to loading f(n) = sin(nω ̃) with ω ̃ = θ1 = 2π/3. (a)

Middle of the first layer located at ξ = 1/4. (b) Middle of the second layer located at

ξ = 3/4. .......................................................................................... 20

Figure 5. The values of the two stress sequences (a) s1(n), (b) s2(n), obtained from the

recursive equations 2 and the analytical equation 38, for the two-layered Goupillaudtype strip with impedance ratio α1 = 1/3, τ = 1, subjected to loading f(n) = sin(nω ̃)

with ω ̃ = θ1 = 2π/3. ............................................................................ 20

Figure 6. Stress time history for a three-layered Goupillaud-type strip with α1 = 3, α2 =

1

2(1−cos(π/4)) − 1, θ1 =

π

4

, τ = 1, subjected to loading f(n) = sin(nω ̃) with ω ̃ =

π

3

. (a)

Middle of the second layer located at ξ = 1/2. (b) Middle of the third layer located at

ξ = 5/6. .......................................................................................... 23

Figure 7. The three-layered Goupillaud-type strip with α1 = 3, α2 =

1

2(1−cos(π/4)) − 1,

θ1 =

π

4

, τ = 1, subjected to loading f(n) = sin(nω ̃) with ω ̃ =

π

3

. (a) The values

of the stress sequence s1(n), obtained from the recursive equations 2 and the analytical

equation 30. (b) The difference of the expressions for s1(n) obtained from the recursive

equations 2 and the analytical equation 30. .................................................... 23

vi
Figure 8. Stress time history for a three-layered Goupillaud-type strip with α1 = 3, α2 =

1

2(1−cos(π/4)) − 1, τ = 1, subjected to loading f(n) = cos(nω ̃) with ω ̃ = π. (a) Middle

of the second layer located at ξ = 1/2. (b) Middle of the third layer located at ξ = 5/6. .. 24

Figure 9. Stress time history for a three-layered Goupillaud-type strip with α1 = 3, α2 =

1

2(1−cos(π/4)) − 1, τ = 1, subjected to loading f(n) = sin(nω ̃) with ω ̃ = θ1 =

π

4

. (a)

Middle of the second layer located at ξ = 1/2. (b) Middle of the third layer located at

ξ = 5/6. .......................................................................................... 26

Figure 10. (a) The values of the stress sequence s1(n), obtained from the recursive equations 2 and the analytical equation 38. (b) The difference of the expressions for s1(n)

obtained from the recursive equations 2 and the analytical equation 38 for a three-layered

Goupillaud-type strip with α1 = 3, α2 =

2

2(1−cos(π/4)) − 1, τ = 1, subjected to loading

f(n) = sin(nω ̃) with ω ̃ = θ1 =

π

4

. ............................................................. 27

Figure 11. Stress time history for a two-layered Goupillaud-type strip with impedance ratio

α1 = 3, θ1 =

π

3

, τ = 1, at the middle of the second layer located at ξ = 3/4 subjected

to loading (a) f(n) = sin(nω ̃) with ω ̃ = θ =

π

3

. (b) f(n) = 1 for n ≥ 0. ................... 36

vii
List of Tables

Table E-1. Optimal base angle values for the m-layer case. ....................................... 53

viii
Acknowledgments

We would like to thank Dr. Steven Segletes for providing the LATEX style files and the U.S. Army

Research Laboratory (ARL) report template for conversion of the original manuscript to ARL

report style. Dr. Segletes has also assisted in solving a variety of LATEX-related format issues in

the Table of Contents, Figures, and Appendices.

ix
INTENTIONALLY LEFT BLANK.

x
1. Introduction

The determination of the natural modes of vibration in layered elastic media is a central problem

of geophysics since earthquake-induced free oscillations of the Earth are often used to constrain

models of the Earth’s interior. Poisson (1) was the first to determine the free radial vibrations of

a homogeneous sphere, and Lamb (2) calculated some of its natural frequencies; the introduction

in Love’s Mathematical Treatise of Elasticity (3) provides an historical account of the early

developments in this field. It was only with the advent of high-speed computers that

planetary-scale vibration problems could be solved with finite, time-dependent sources in

compressible, self-gravitating, inhomogeneous, and/or layered oblate spheroids of revolution

(4, 5). Interestingly, the observed resonance in the free oscillations of the solid Earth, has been

attributed to fluid-structure interaction induced by the atmosphere, rather than the cumulative

effect of microearthquakes (6).

Studies of natural vibrations in elastic media at spatial scales of technological interest focus on

resonance in layered structures (7), anisotropic elastic bodies (8), anisotropic layered crystals

(9, 10), elastic plates (11, 12), periodic media (13), laminated and sandwich plates (14),

composite laminates (15), piezoelectric composites (16), locally resonant acoustic metamaterials

(17, 18), and acoustic wave devices that exploit the phenomenon (19–21). Thus, it is seen that

resonance can enhance the performance of many sensors and devices, yet can damage structures

subjected to sustained, temporally periodic loading (22).

Although there are numerous works that establish the natural modes of vibration of multilayered

media (23, 24), and analogous problems for n-segmented strings (25) and n-strings (26),

analytical solutions for the transient resonance response of multilayered media have been limited

to analyses involving only a few layers. Churchill (27) uses Laplace transform methods to derive

the transient resonance response of a free-fixed elastic bar, and notes the presence of the product

of a temporal term and time-harmonic function in the solution, which becomes unbounded at

large times; similar unbounded behavior is observed in the interaction of a compressible fluid

with an elastic plate (28). Caviglia and Morro provide closed-form, time-domain expressions for

transient waves in a multilayered elastic (possibly anisotropic) medium (29). They derive the

transient resonance response of a single isotropic elastic layer sandwiched between two

half-spaces that is subjected to a temporally periodic sawtooth function, and observe that

resonance makes the elastic layer acoustically transparent. A similar two-dimensional problem is

studied by Kaplunov and Krynkin (30), who examine the influence of layer stiffness and

1
thickness on the resonance vibrations in a symmetrically point-loaded single layer of elastic

material enclosed between two half-spaces. Qiang et al. (31) study the natural frequencies of

anisotropic multilayers, and illustrate beat phenomena in two and three layer systems using an

efficient numerical eigensolution scheme that is based on semi-analytical methods. Exact

expressions for the reflection coefficients for a two-dimensional elastic layer overlying an elastic

half-space are obtained by Fokina and Fokin (32), but the transient response of the system at

resonance is not analyzed. Graff presents a general method (33) to determine the natural

frequencies of composite rods for power ultrasonic applications with a specific solution for a

two-rod system. The method can be extended to multilayered systems, but the general problem

of determining the natural frequencies with specific geometry and material properties can only be

solved using numerical methods. Exact analytical expressions are derivable, however, for the

transient resonance response of a multilayered medium if we specialize the medium to be of

Goupillaud-type (34). The Goupillaud specialization is often used in geophysical applications to

model wave propagation in inhomogeneous media (35). The present treatment uses a global

matrix method that is attributed to Knopoff (36), rather than the Thomsen-Haskell transfer matrix

formalism (37, 38). Since the stresses in the multilayered medium are written in terms of

recursion relations that interlace, only half the number of equations are required relative to those

in classical global matrix methods (36). Finally, Bube and Burridge (39) treat the inverse

problem of finding the coefficients of a first order 2 × 2 hyperbolic system related to reflection

seismology. In order to numerically solve the continuous initial-boundary-value problem, several

difference schemes are applied as discretizations to the corresponding differential equations. The

difference scheme IVp, given in equation (3.1.11), pg. 517 of Bube and Burridge (39), appears

to match our recursive relations for the stress values given in section 2. Our recursive relations

are exact for discretely layered media of Goupillaud-type, and the stress solutions that are derived

herein using z-transform methods, are then used to determine the resonance response for the

discrete model, and subsequently extended to the continuous model.

The remainder of the report is organized as follows. The resonance frequency spectrum for the

discrete model is developed in section 2 by generalizing the global system of recursion relations

derived in Velo and Gazonas (40). Explicit analytical expressions for the stress in a multilayered

elastic medium are obtained using z-transform methods and are quantitatively compared with the

recursion relation solutions at resonance. These results are then extended in section 3 to

analytically describe the natural frequency spectrum of a free-fixed Goupillaud-type strip, also

obtained using a generalization of Graff’s method. Various applications and interpretations of the

frequency results are discussed in section 4. Section 5 summarizes the results, while

miscellaneous references and derivations are included in the appendices, including a proof that

2
the zeros of the palindromic polynomial, representing the determinant of the global system

matrix, lie on the unit circle for 1 ≤ m ≤ 5.

2. Forced Vibrations and Resonance Frequencies: The Discrete Model

A finite strip made of m-layers of homogeneous elastic materials is considered, where the density

and elastic modulus of the material occupying the ith layer are represented by ρi and Ei

respectively, for i = 1, . . . , m. Using the definition of the wave speed c = (E/ρ)

1/2

, we relate to

the ith layer the wave speed ci for i = 1, 2, . . . , m. The characteristic impedance in the ith layer

is given as a product of the material density ρi and wave speed ci, while

αi =

ciρi

ci+1ρi+1

=

√

√

Eiρi

Ei+1ρi+1

represents the impedance ratio between layers i and i + 1 for

i = 1, . . . , m − 1. The m-layered strip is of Goupillaud-type, which means equal wave travel

time for each layer. By replacing the spatial variable x of the wave equation with the new

variable ξ =

R x

0

ds

c(s)

, the problem can be converted to a Goupillaud-type strip with the same length

L as the transit time τ through the strip, equal layer lengths of τ

m

, and an equal wave travel time of

τ

m

for each layer in either direction, which further implies (equal) wave speed of unity for each

layer. Here c = c(s) is the piecewise constant wave speed function taking values c1, c2, . . . , cm in

each layer, respectively. In sections 2 and 3, it is shown that the stress values and impedance

ratios remain unchanged as a result of this simplification. Therefore, although the spatial variable

transformation mentioned previously simplifies the problem, and is used to analyze a

Goupillaud-type strip with length τ and unit wave speed in each layer, the results can be

applied/interpreted for a Goupillaud-type strip of unequal layer lengths.

In prior work (40), the strip was subjected to zero initial conditions, and a Heaviside step in stress

loading p = constant at time t = 0 at one end, and held fixed at the other end. This is illustrated

in figure 1 with f(n) = p for n ≥ 0. The present work generalizes this system by replacing the

Heaviside step in stress loading with a discrete loading function f(n), n ≥ 0. After transforming

the spatial variable x of the wave equation to the new variable ξ =

R x

0

ds

c(s)

, the initial-boundary

3
value problem becomes







∂

2u(ξ,t)

∂t2 =

∂

2u(ξ,t)

∂ξ2 for (i−1)τm < ξ < iτm

, i = 1, . . . , m,

σ(0, t) = E1

c1

∂u

∂ξ (0, t) = E ̃

1

∂u

∂ξ (0, t) = f(n) for (n − 1) 2τm < t < n2τm

, n ≥ 1,

u(τ, t) = 0,

u(ξ, 0) = ∂u

∂t (ξ, 0) = 0.

(1)

The functions u and σ represent the displacement and stress, respectively, while f(n) represents

the discrete applied loading. Here, Ei and ρi represent the actual material properties for the

physical case, while ρ ̃i and E ̃

i represent the material properties of the ith layer for the simplified

case. We will use this notation whenever we want to distinguish the physical case from the

simplified case. Having ρ ̃i = E ̃

i =

Ei

ci

=

√

Eiρi, it implies that the impedance ratio

αi =

√

√

Eiρi

Ei+1ρi+1

=

√

E ̃ √ iρ ̃i

E ̃

i+1ρ ̃i+1

=

E ̃

i

E ̃

i+1

, as well as the stress values E1

∂u

∂x =

E1

c1

∂u

∂ξ = E ̃

1

∂u

∂ξ , remain

the same for i = 1, 2, . . . , m − 1. As seen in figure 1, the previous coordinates of the boundary

and layer interfaces x0 = 0, xj, and xm = L are replaced by the new coordinates ξ0 = 0, ξj =

jτ

m

,

and ξm = τ for j = 1, 2, . . . , m. The time variable t is represented on the vertical axis and τ

m

represents the (equal) wave travel time for each layer of the m-layered strip in either direction. In

figure 1, the stress jump-discontinuities are propagated along the dashed characteristic lines of the

wave equation, while the inner vertical solid lines represent the layer interfaces. Due to the zero

initial conditions, the stress values below the first characteristic line segment with equation t = ξ

for 0 < ξ < τ are zero. Above it, however, the intersection of characteristics with each other and

the boundaries splits the region into square subregions. Due to the continuity of stress and

displacement at each layer interface, the stress value is constant in each square subregion and is

represented by si(n) for i = 1, 2, . . . , m and n ≥ 1. In a given layer i, the stress values alternate

(interlace) between si−1(n) and si(n) for n = 1, 2, . . ., i = 1, 2, . . . , m, and s0(n) = f(n). When

positioned in the middle of the ith or (i + 1)th layer, the time interval at which si(n) is reached is

4n+2i−3

2m < t < 4n+2i−12m

, as depicted later in figures 2, 3, 4, 6, 8, 9, and 11. Here n ≥ 1,

i = 0, 1, 2, . . . , m, and s0(n) = f(n).

4
 



 

 

 

 

  

	
 

	
 

	 

	
 



	
 

	
 



 





 

!

 

 

" #  









$ 

#	
 



#	
 



#

	
 



#	
 



#

	
 



#

	 















	 

  

"	 

  

# " 	 


  

" 	 


  

 

      

  









	
 



#	
 



	
 



Figure 1. Lagrangian diagram for an elastic strip made of m-layers of equal wave

travel time.

From the superposition of the transmitted and reflected waves at the layer interfaces and

boundary, a system of coupled first order difference equations for the stress terms was derived in

Velo and Gazonas (40); the system was then solved using z-transform methods. Herein, we

generalize this system of equations by incorporating an arbitrary discrete loading function f(n),

n ≥ 0 which appears in the first line of equation 2:







s1(n + 1) = −s1(n) + 2α1

1+α1

s2(n) + 2

1+α1

f(n + 1),

s2(n + 1) = −s2(n) + 2α2

1+α2

s3(n) + 2

1+α2

s1(n + 1),

. . .

si(n + 1) = −si(n) + 2αi

1+αi

si+1(n) + 2

1+αi

si−1(n + 1),

. . .

sm−1(n + 1) = −sm−1(n) + 2αm−1

1+αm−1

sm(n) + 2

1+αm−1

sm−2(n + 1),

sm(n + 1) = −sm(n) + 2sm−1(n + 1),

(2)

5
subject to zero-stress initial conditions si(0) = 0 for 1 ≤ i ≤ m. Here n ≥ 0 is a suitable time

index and si(n) represent the stress terms for 1 ≤ i ≤ m (see figure 1). The discrete function

f(n), n ≥ 0, represents the stress loading applied at ξ = 0. Notice that by applying the limiting

condition αi → 0 to the recursive relation si(n + 1) = −si(n) + 2αi

1+αi

si+1(n) + 2

1+αi

si−1(n + 1)

for i = m, we recover the last line of equation 2, corresponding to the fixed end. A similar

system of difference equations may be obtained for other initial and boundary conditions. For

example, a free-free multilayered system can be studied by setting the last line in equation 2 to

sm(n) = sm(0) = 0, and a multilayer on a half-space can be obtained using

sm(n + 1) = sm−1(n + 1), and the rest of the equations for 1 ≤ i < m remain the same.

Equation 2 can then be solved using z-transform methods. Indeed, due to the zero-stress initial

conditions si(0) = 0 for 1 ≤ i ≤ m, the z-transform1 of si(n + 1) becomes

Z(si(n + 1)) = zSi(z) − zsi(0) = zSi(z), (3)

where Si(z) is the z-transform of si(n), while

Z(f(n + 1)) = zF(z) − zf(0) = z(F(z) − f(0)), (4)

where F(z) is the z-transform of f(n). As a result, the z-transform of equation 2 is expressed as







zS1(z) = −S1(z) + 2α1

1+α1

S2(z) + 2

1+α1

z(F(z) − f(0)),

zS2(z) = −S2(z) + 2α2

1+α2

S3(z) + 2

1+α2

zS1(z),

. . .

zSi(z) = −Si(z) + 2αi

1+αi

Si+1(z) + 2

1+αi

zSi−1(z),

. . .

zSm−1(z) = −Sm−1(z) + 2αm−1

1+αm−1

Sm(z) + 2

1+αm−1

zSm−2(z),

zSm(z) = −Sm(z) + 2zSm−1(z).

(5)

After reorganizing the terms in equation 5, the linear system is written in matrix-vector form as

Am ̃xm = b ̃m, (6)

1The definition for the z-transform of g(n), n ≥ 0, is Z(g(n)) = G(z) = P∞

n=0

g(n)z

−n for |z| > R in the complex

plane (see Jury (41)).

6
where Am is a tridiagonal matrix given as

Am =

















































z + 1 −η1α1 0 · · · 0 0 0 0

−η2z z + 1 −η2α2 0 · · · 0 0 0

0 −η3z z + 1 −η3α3 0 · · · 0 0

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

0 0 0 · · · 0 −ηm−1z z + 1 −ηm−1αm−1

0 0 0 · · · 0 0 −ηmz z + 1

















































m×m

, (7)

while

 ̃xm =

































S1(z)

S2(z)

.

.

.

Sm(z)

































m×1

, b ̃m =

































η1z(F(z) − f(0))

0

.

.

.

0

































m×1

, ηi =

2

1+αi

for i = 1, . . . , m − 1, ηm = 2. (8)

Due to the sparseness of b ̃m, the solution of the linear system in equation 6 is

 ̃xm = Am

−1b ̃m =

η1z(F(z)−f(0))

|Am|





























(−1)1+1|A1,1|

(−1)1+2|A1,2|

.

.

.

(−1)1+m|A1,m|





























=

η1z(F(z)−f(0))

|Am|





























|A1,1|

−|A1,2|

.

.

.

(−1)1+m|A1,m|





























. (9)

Here |Am| is the determinant of Am, A1,i for i = 1, . . . , m are minors of Am, and |A1,i| are their

corresponding determinants. The following recursive relation (see e.g., Muir (42)) holds

|Am| = (z + 1)|Am−1| − η1η2α1z|Am−2|, (10)

7
where Am corresponds to the m-layered strip with impedance ratios α1, α2, . . . , αm−1; and Am−1

corresponds to the remaining (m − 1)-layered strip with impedance ratios α2, α3, . . . , αm−1,

obtained after the removal of the first layer from the m-layered strip; Am−2 corresponds to the

remaining (m − 2)-layered strip with impedance ratios α3, . . . , αm−1, obtained after the removal

of the first two layers from the m-layered strip. One can derive by induction that the determinant

|Am| is a palindromic polynomial with real coefficients, i.e., the coefficients in front of z

m−j

and

z

j

are real and equal to each-other for j = 0, . . . , m and m ≥ 1:

|Am| =







z

m + am,1zm−1 + am,2zm−2 + · · · + am, m

2 −1z

m

2 +1 + am, m

2

z

m

2 +

+am, m

2 −1z

m

2 −1 + · · · + am,2z

2 + am,1z + 1, for m-even,

(z + 1)[z

m−1 + bm,1zm−2 + · · · + bm, m−1

2 −1

z

m−1

2 +1 + bm, m−1

2

z

m−1

2 +

+bm, m−1

2 −1

z

m−1

2 −1 + · · · + bm,2z

2 + bm,1z + 1], for m-odd.

(11)

The complex roots of a polynomial with real coefficients occur in conjugate pairs, and if the

coefficients are palindromic, then the roots occur in inverse pairs. Since the inverse pairing is not

necessarily the same as the conjugate pairing, if a root’s inverse is not its conjugate, then it must

be one of a set of four related roots that satisfy a palindromic quartic with real coefficients. As a

result, each palindromic polynomial of even degree in equation 11 is expected to factor into

quadratic and quartic palindromes with real coefficients. In appendix B, the quadratic and quartic

factors of |Am| for 1 ≤ m ≤ 5 are shown to satisfy the necessary and sufficient conditions for all

their roots to be distinct (complex conjugate), and on the unit circle. As a result, all the zeros of

|Am| for 1 ≤ m ≤ 5 are on the unit circle, and their corresponding angles are provided in

appendix D.

In this report, we only consider m-layered designs for which all the zeros of |Am| are distinct,

have no multiples, and are on the unit circle. Based on the results from appendix B, this class of

designs is not empty. Under this assumption, zk = e

Iθk = cos θk + I sin θk, and

z

−1

k =  ̄zk = e

−Iθk = cos θk − I sin θk, where z ̄k is the complex conjugate of zk and 1 ≤ k ≤ bm

2

c.

From the fact that

zk + z

−1

k = zk +  ̄zk = 2 cos θk, (12)

a newly factored representation of |Am| is obtained, with linear factor (z + 1) for m-odd, and only

8
palindromic quadratics with real coefficients involving the cosines of the angles θ1, θ2, . . . , θb m

2

c,

|Am| =







b m

2Q

c

k=1

[z

2 − 2z cos θk + 1] for m even,

(z + 1)

b m

2Q

c

k=1

[z

2 − 2z cos θk + 1] for m odd.

(13)

Based on the factorization of the determinant |Am| given in equation 13, for any m-layered

design, if there are m distinct roots on the unit circle, then there are m distinct angles. The angles

{±θk}

m

2

k=1 correspond to the roots z = e

±Iθk for m even, while the angles θ0 = π, {±θk}

b m

2

c

k=1

correspond to the roots z = e

Iθ0 = −1 and z = e±Iθk , respectively, for m odd. The b

m

2

c

essential angles 0 < θk < π for k = 1, . . . , b

m

2

c were referred to as the base angles in Velo and

Gazonas (40). The equations for the angles θ0 = π and {θk}

b m

2

c

k=1 in terms of the design

parameters (combinations of the impedance ratios) can be derived from the condition:

|Dm| = 0. (14)

The definition for the matrix Dm as well as the angle-design parameter equations for up to five

layers are included in appendices C–D.

In the rest of this section, we prove that the resonance frequency spectrum for the discrete model

consists of all the positive angles coterminal with θ0 = π or {±θk}

b m

2

c

k=1 , as shown below







ω ̃l0 = θ0 + 2l0π = (2l0 + 1)π (for m odd only),

ω ̃l1,k = θk + 2l1π, l0, l1, l2 = 0, 1, 2 . . . ,

ω ̃l2,k = −θk + 2(l2 + 1)π, k = 1, . . . , b

m

2

c.

(15)

Equivalently, the resonance frequency ω ̃ is expected to satisfy one of the following relations:







cos  ̃ω = −1 (for m odd only),

cos  ̃ω = cos θk for k = 1, . . . , b

m

2

c.

(16)

Since the degree of |Am| is m and the degree of its minors |A1,i| is m − 1, for i = 1, 2, . . . , m,

the substitution of equation 13 into equation 9 results in the following expansion of the

components  ̃xm(i) of  ̃xm into partial fractions:

9
 ̃xm(i) = z(F(z) − f(0)) h

η1(−1)i|A1,i|

|Am|

i

= z(F(z) − f(0)) "

bi,0

z+1 +

b m

2P

c

k=1

a

∗

i,kz+b

∗

i,k

z

2−2z cos θk+1 #

= (F(z) − f(0)) "

bi,0z

z+1 +

b m

2P

c

k=1

a

∗

i,kz

2+b∗

i,kz

z

2−2z cos θk+1 #

,

(17)

where bi,0 = 0 for m even and i = 1, 2, . . . , m. Applying the inverse z-transform, one can derive

that

Z

−1

"

bi,0z

z+1 +

b m

2P

c

k=1

a

∗

i,kz

2+b∗

i,kz

z

2−2z cos θk+1 #

= Z

−1

"

bi,0· z

z+1 +

b m

2P

c

k=1

3

ai,k· z(z−cos θk)

z

2−2z cos θk+1 +

bi,k·z sin θk

z

2−2z cos θk+1 ́

#

= bi,0(−1)n +

b m

2P

c

k=1

ai,k cos (nθk) + bi,k sin (nθk),

(18)

for the following choices of the coefficients and i = 1, 2, . . . , m,

ai,k = a

∗

i,k, bi,k =

a

∗

i,k cos θk + b

∗

i,k

sin θk

for k = 1, . . . , b

m

2

c. (19)

From the fact that 0 < θk < π it follows that sin θk 6= 0 for all k = 1, . . . , b

m

2

c. As a result, after

applying the inverse z-transform to equation 17 and based on equation 18, a general

representation for the stress terms in equation 2 is obtained,

si(n) = Z

−1

( ̃xm(i)) = f(n) ∗

"

bi,0(−1)n +

b m

2P

c

k=1

ai,k cos (nθk) + bi,k sin (nθk)

#

−f(0) ·

"

bi,0(−1)n +

b m

2P

c

k=1

ai,k cos (nθk) + bi,k sin (nθk)

#

,

(20)

where operation ∗ means convolution, i = 1, 2, . . . m, and n ≥ 0. The stress representation

previously derived in Velo and Gazonas (40) can now be recovered for the special choice of

loading f(n) = p for n ≥ 0 (see appendix A for further details). The stress representation in

equation 20 involves trigonometric sums which relate to the Chebyshev polynomials of the first

and second kind (see subsection 4.4).

10
In the work that follows, we use equations 17–18 to obtain stress solutions when a harmonic

forcing function of the form f(n) = sin(nω ̃) with f(0) = 0 and n ≥ 0 is applied, unless

mentioned otherwise. A similar approach can be used for a harmonic forcing function of the

form f(n) = cos(nω ̃) for n ≥ 0. As discussed in Proakis and Manolakis (43), the nature of time

(continuous or discrete) is expected to affect the nature of the frequency. The frequency for the

analog (continuous) time signal is related to the frequency of the digital (discrete) time signal

through a linear transformation involving the time step or the so-called sampling interval. There

is a linear transformation in frequency

ω ̃ =

2τ

m

ω, (21)

from the continuous forcing function with frequency ω to the discrete forcing function with

frequency ω ̃, and time step 2τ

m

:

sin(tω) = sin μ

t

2τ

m

·

2τ

m

ω

¶

→ sin μd

t

2τ

m

e · 2τ

m

ω

¶

= sin μn ·

2τ

m

ω

¶

= sin(nω ̃). (22)

As seen in figure 1, t ≥ 0 represents the time variable, n = d

t

2τ

m

e = 0, 1, 2, . . . represents a

time-related index, τ represents the wave transit time through the layered strip, and m represents

the number of layers.

For the choice of f(n) = sin(nω ̃), n ≥ 0, unless mentioned otherwise, we analytically prove the

equivalent equations 15 and 16, and validate them through various numerical experiments.

Substituting f(0) = 0 into equation 17, and

F(z) = Z(f(n)) = Z(sin(nω ̃)) = z sin( ̃ω)

(z

2 − 2z cos  ̃ω + 1), (23)

after using equation 18, it can be shown that

si(n) = Z

−1

( ̃xm(i))

= Z

−1

"

z sin( ̃ω)

(z

2−2z cos  ̃ω+1) ·

Ã

bi,0z

z+1 +

b m

2P

c

k=1

a

∗

i,kz

2+b∗

i,kz

z

2−2z cos θk+1 ! #

si(n) = bi,0Z

−1

h

z

2

sin( ̃ω)

(z

2−2z cos  ̃ω+1)(z+1) i

+

b m

2P

c

k=1

Z

−1

h

z sin( ̃ω)·(a

∗

i,kz

2+b∗

i,kz)

(z

2−2z cos θk+1)(z2−2z cos  ̃ω+1) i

,

(24)

or equivalently,

si(n) = bi,0Z

−1

h

z

2

sin( ̃ω)

(z

2−2z cos  ̃ω+1)(z+1) i

+

b m

2P

c

k=1

Z

−1

h

z sin( ̃ω)·(a

∗

i,kz

2+b∗

i,kz)

(z

2−2z cos θk+1)(z2−2z cos  ̃ω+1) i

. (25)

11
In the subsections that follow, we calculate the inverse z-transform of both expressions in

equation 25 in order to derive the stress formulas at resonance and non-resonance frequencies.

2.1 Stress Solutions at Non-resonance Frequencies

Here we show that the values of the driving frequency ω ̃, which satisfy cos  ̃ω 6= −1 = cos θ0 (for

m odd) and cos  ̃ω 6= cos θk for all k = 1, . . . , b

m

2

c, represent non-resonance frequency values. In

the set of real numbers, this set of frequency values is the complement of the set of values

identified in equation 16. We begin by finding the inverse z-transform of the two expressions

given in equations 27 and 29.

When cos  ̃ω 6= −1, the partial fraction expansion

z

2

sin( ̃ω)

(z

2−2z cos  ̃ω+1)(z+1) =

sin  ̃ω

2(1+cos  ̃ω)

·

h

z

2+z

z

2−2z cos  ̃ω+1 −

z

z+1 i

=

1

2(1+cos  ̃ω)

·

h

z(z−cos  ̃ω) sin  ̃ω

z

2−2z cos  ̃ω+1 +

z sin  ̃ω(1+cos  ̃ω)

z

2−2z cos  ̃ω+1 − sin  ̃ω

z

z+1 i

,

(26)

implies that

Z

−1

h

z

2

sin( ̃ω)

(z

2−2z cos  ̃ω+1)(z+1) i

=

sin  ̃ω

2(1+cos  ̃ω)

cos(nω ̃) + 1

2

sin(nω ̃) −

sin  ̃ω

2(1+cos  ̃ω)

(−1)n. (27)

When cos  ̃ω 6= cos θ, the partial fraction expansion,

z sin  ̃ω(az2 + bz)

(z

2 − 2z cos  ̃ω + 1)(z2 − 2z cos θ + 1) =

Az(z − cos  ̃ω)

z

2 − 2z cos  ̃ω + 1

+

Bz sin  ̃ω

z

2 − 2z cos  ̃ω + 1

+

Cz(z − cos θ)

z

2 − 2z cos θ + 1

+

Dz sin θ

z

2 − 2z cos θ + 1

, (28)

implies that

Z

−1

·

z sin  ̃ω(az2 + bz)

(z

2 − 2z cos  ̃ω + 1)(z2 − 2z cos θ + 1)  ̧

= A cos(nω ̃) + B sin(nω ̃)

+C cos(nθ) + D sin(nθ), (29)

where A = −C =

a sin  ̃ω

2(cos  ̃ω−cos θ)

, B =

a cos  ̃ω+b

2(cos  ̃ω−cos θ)

, and D = −

(a cos θ+b) sin  ̃ω

2 sin θ(cos  ̃ω−cos θ)

. Substituting the

results from equations 27 and 29 into equation 25, the following stress solutions for

12
cos  ̃ω 6= −1 = cos θ0 and cos  ̃ω 6= cos θk for all k = 1, . . . , b

m

2

c are obtained:

si(n) = −

bi,0 sin  ̃ω

2(1+cos  ̃ω)

(−1)n +

"

bi,0 sin  ̃ω

2(1+cos  ̃ω) +

b m

2P

c

k=1

Ai,k #cos(nω ̃)

+

"

bi,0

2 +

b m

2P

c

k=1

Bi,k #sin(nω ̃) +

b m

2P

c

k=1

Ci,k cos(nθk) +

b m

2P

c

k=1

Di,k sin(nθk),

(30)

where

Ai,k = −Ci,k =

a

∗

i,k sin  ̃ω

2(cos  ̃ω−cos θk)

, Bi,k =

a

∗

i,k cos  ̃ω+b

∗

i,k

2(cos  ̃ω−cos θk)

, Di,k = −

(a

∗

i,k cos θk+b

∗

i,k) sin  ̃ω

2 sin θk(cos  ̃ω−cos θk)

. (31)

As before, bi,0 = 0 when the number of layers m is even. For a given value of ω ̃, the stress

solutions in equation 30, including all the coefficients and angles θk for k = 1, . . . , b

m

2

c, can be

also expressed in terms of the impedance ratios only (see equations 6–9, 14, and appendices D

and E).

Given that 0 < θk < π implies that sin θk 6= 0 for all k = 1, . . . , b

m

2

c. This, together with the fact

that cos  ̃ω 6= −1 and cos  ̃ω 6= cos θk for all k = 1, . . . , b

m

2

c, implies that the coefficients in the

stress solutions are always bounded (see equations 30–31). Therefore, for a given m-layered

design and driving frequency ω ̃, the stress solutions in equations 30 given by trigonometric sums

of sines and cosines are bounded. However, as cos  ̃ω approaches −1 = cos θ0 (for m odd) or any

of the cos θk values, k = 1, . . . , b

m

2

c, the denominators of some coefficients of the stress solutions

in equations 30–31 get close to zero, allowing these solutions to take larger values and approach

resonance.

2.2 Stress Solutions at Resonance Frequencies

The stress solutions in equation 30 suggest that the driving frequencies ω ̃ identified in equation 16

are resonance frequencies. Here we prove equation 16 and its equivalent form equation 15, by

deriving the stress solutions at the expected resonance frequencies, demonstrating that these

solutions grow without bound over time.

2.2.1 Case I: cos  ̃ω = −1 = cos θ0 (for m odd only)

The fact that cos( ̃ω) = −1 or equivalently that ω ̃ takes values ω ̃l0 = (2l0 + 1)π for l0 = 0, 1, 2 . . .,

implies that f(n) = sin(nω ̃) = sin(nπ) = 0 for all n ≥ 0. As a result, due to the zero input in

loading we are unable to detect resonance, as confirmed by equation 20. However, by choosing

f(n) = cos(nω ̃) = cos(nπ) we overcome this obstacle. Based on equations 17–20, we expect

13
any unbounded stress terms to come from

F(z) ·





bi,0z

z + 1

+

b m

2 X

c

k=1

a

∗

i,kz

2 + b∗

i,kz

z

2 − 2z cos θk + 1



 , (32)

after the inverse z-transform is applied. Here F(z) = Z(f(n)) = Z(cos((nπ)) = z

z+1 , which

after substituting into equation 32, gives

z

z+1 ·

"

bi,0z

z+1 +

b m

2P

c

k=1

a

∗

i,kz

2+b∗

i,kz

z

2−2z cos θk+1 #

=

bi,0z

2

(z+1)2 +

b m

2P

c

k=1

z·(a

∗

i,kz

2+b∗

i,kz)

(z+1)(z

2−2z cos θk+1)

=

bi,0z

2

(z+1)2 +

b m

2P

c

k=1

a

∗

i,kz

2+b∗

i,kz

z

2−2z cos θk+1 −

b m

2P

c

k=1

a

∗

i,kz

2+b∗

i,kz

(z+1)(z

2−2z cos θk+1) .

(33)

Furthermore, for the last sum above, we have that

b m

2P

c

k=1

a

∗

i,kz

2+b∗

i,kz

(z+1)(z

2−2z cos θk+1)

=

b m

2P

c

k=1

a

∗

i,k

sin θk

·

z

2

sin θk

(z+1)(z

2−2z cos θk+1) +

b m

2P

c

k=1

b

∗

i,k ·

z

(z+1)(z

2−2z cos θk+1)

=

b m

2P

c

k=1

a

∗

i,k

sin θk

·

z

2

sin θk

(z+1)(z

2−2z cos θk+1) +

b m

2P

c

k=1

b

∗

i,k

2 cos θk

·

h

z

z+1 +

−z

2+(1+2 cos θk)z

z

2−2z cos θk+1 i

.

(34)

Applying the inverse z-transform to the sums in equations 33–34 and using equations 18 and 27,

we conclude that the only term in equation 33 that generates unbounded stress values is bi,0z

2

(z+1)2 for

bi,0 6= 0, due to the fact that

Z

−1

·

bi,0z

2

(z + 1)2

 ̧

= bi,0 · Wn(−1) = bi,0(−1)n(n + 1). (35)

Here Wn(y) represents the Chebyshev polynomial of the second kind evaluated at y = −1. The

presence of resonance for the choice of frequency values ω ̃l0 = (2l0 + 1)π, l0 = 0, 1, 2, . . . ,

proves that they are part of the resonance frequency spectrum identified in equation 15. Notice

that the resonance frequency spectrum ω ̃l0 = (2l0 + 1)π for l0 = 0, 1, 2, . . . is universal for all

the designs with an odd number of layers as it is independent of any design parameters (e.g.,

impedance ratios). This is illustrated in figure 2 for the 7- and 11-layer case and later on in

figure 8 for the 3-layer case.

14
0 1 2 3 4 5 6 7 8

−6

−4

−2

0

2

4

6

8

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

0 1 2 3 4 5 6 7 8

−10

−8

−6

−4

−2

0

2

4

6

8

10

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

(a) (b)

Figure 2. Stress time history for a layered strip with impedance ratios α1 = 3,

α2 = 2, α3 = 1.5, α4 = 2.2, α5 = 0.3, α6 = 1.7, α7 = 1.4, α8 = 3.1,

α9 = 0.8, and α10 = 4 as applicable; τ = 1 and loading

f(n) = cos(nω ̃) applied at ξ = 0 for ω ̃ = π. (a) The middle of the

second layer of a 7-layered strip located at ξ = 3/14. (b) The middle

of the second layer of an 11-layered strip located at ξ = 3/22.

2.2.2 Case II: cos  ̃ω = cos θk∗

The following expansion

z sin θ(az2+bz)

(z

2−2z cos θ+1)2 =

Az ̃ [−2z+(1+z

2

) cos θ

(z

2−2z cos θ+1)2 ] + Bz ̃ (z

2−1) sin θ

(z

2−2z cos θ+1)2 +

Cz ̃ sin θ

z

2−2z cos θ+1 (36)

implies that

Z

−1

·

z sin θ(az2 + bz)

(z

2 − 2z cos θ + 1)2

 ̧

= A ̃ · n cos(nθ) + B ̃ · n sin(nθ) + C ̃ sin(nθ), (37)

where A ̃ = −

(a cos θ+b)

2 sin θ

, B ̃ =

a

2

, and C ̃ =

a+b cos θ

2 sin2 θ

.

2.2.3 Case II.a.: sin  ̃ω = sin θk∗

Applying equations 27, 29, and 37 to equation 25, the following stress solutions for the frequency

spectrum ω ̃l1,k∗ = θk

∗ + 2l1π, l1 = 0, 1, 2, . . . are obtained,

15
si(n) = −

bi,0 sin θk∗

2(1+cos θk∗ )

(−1)n +

"

bi,0 sin θk∗

2(1+cos θk∗ ) + nA ̃

i,k∗ +

b m

2P

c

k=1,k6=k

∗

Ai,k #cos(nθk

∗ )

+

"

bi,0

2 + nB ̃

i,k∗ + C ̃i,k∗ +

b m

2P

c

k=1,k6=k

∗

Bi,k #sin(nθk

∗ )

+

b m

2P

c

k=1,k6=k

∗

Ci,k cos(nθk) +

b m

2P

c

k=1,k6=k

∗

Di,k sin(nθk).

(38)

Here Ai,k, Bi,k, Ci,k, Di,k are given from equation 31 for ω ̃ = θk

∗ + 2l1π, l1 = 0, 1, 2, . . ., while

A ̃

i,k∗ = −

(a

∗

i,k∗ cos θk

∗ + b

∗

i,k∗ )

2 sin θk

∗

, B ̃

i,k∗ =

a

∗

i,k∗

2

, C ̃

i,k∗ =

a

∗

i,k∗ + b

∗

i,k∗ cos θk

∗

2 sin2θk

∗

. (39)

As before, bi,0 = 0 for m even and sin θk 6= 0 for all k = 1, . . . , b

m

2

c. The stress solutions in

equation 38, including all the coefficients and angles θk for k = 1, . . . , b

m

2

c, can be also expressed

in terms of the impedance ratios only (see equations 6–9, equation 14, and appendices D and E).

All the terms in the stress solution in equation 38 are bounded, except for those that are multiplied

by the time-related linear factor n given below

nA ̃

i,k∗ cos(nθk

∗ ) + nB ̃

i,k∗ sin(nθk

∗ ) = −

n

2

·

(a

∗

i,k∗ cos θk

∗ + b

∗

i,k∗ )

sin θk

∗

cos(nθk

∗ ) − a

∗

i,k∗ sin(nθk

∗ )

 ̧

= −

n

2

·

a

∗

i,k∗ cos [(n + 1)θk

∗ ] + b

∗

i,k∗ cos [(nθk

∗ )]

sin θk

∗

 ̧

, (40)

where a

∗

i,k∗ 6= 0 or b

∗

i,k∗ 6= 0 for i = 1, 2, . . . , m. This implies that the stress amplitude is

expected to grow without bound over time, proving resonance for the frequency spectrum

ω ̃l1,k∗ = θk

∗ + 2l1π, l1 = 0, 1, 2, . . ., included in equation 15.

2.2.4 Case II.b.: sin  ̃ω = − sin θk∗ = sin(−θk∗ )

A similar stress representation to equation 38 can be obtained at the frequency values

ω ̃l2,k∗ = −θk

∗ + 2(l2 + 1)π, l2 = 0, 1, 2, . . ., after a sign adjustment for some of the coefficients.

The presence of resonance proves that this frequency spectrum is also part of the resonance

frequency spectrum identified in equation 15.

In summary, we related to an m-layered Goupillaud-type strip with all the zeros of |Am| on the

unit circle, its corresponding set of angles θ0 = π (for m odd only) and {±θk}

b m

2

c

k=1 . In general,

16
the same set of angles could correspond to more than one layered design. When the strip was

subjected to a harmonic forcing f(n) = sin(nω ̃) or f(n) = cos(nω ̃) at one end and held fixed at

the other end, we showed that the positive angles coterminal with θ0 = π (for m odd only) and

{±θk}

b m

2

c

k=1 represented values of the driving frequency ω ̃, which produced resonance (see equation

15 or 16). The complement of this set of values in the set of real numbers was shown to represent

the set of non-resonance frequencies.

In section 2.3 we demonstrate this approach and validate some of the results mainly for the twoand three-layer case. According to appendix B and equations D-1 and D-2 in appendix D, a

Goupillaud-type strip with m = 2, 3 has all the zeros of |Am| on the unit circle, hence a

resonance spectrum given by equation 15. For simplicity, we continue to choose

f(n) = sin(nω ̃), n ≥ 0, unless mentioned otherwise.

2.3 The Two-layer Case (m = 2)

For the two-layer case, it follows from equations 21–22 that ω ̃ = τω and

f(n) = sin(n · τω) = sin(nω ̃). As a result, equation 17 becomes

 ̃x2 =

z sin( ̃ω)

(z

2 − 2z cos  ̃ω + 1) ·









2z

2+2z

χ1(z

2−2z cos θ1+1)

4z

2

χ1(z

2−2z cos θ1+1)









, (41)

where χ1 = 1 + α1 > 1 and cos θ1 =

χ1−2

χ1

6= ±1 from equation D-1.

2.3.1 Stress Solutions at Non-resonance Frequencies

When cos  ̃ω 6= cos θ1, we obtain the following bounded stress solutions







s1(n)

s2(n)





 = Z

−1

( ̃x2)

=

1

χ1(cos  ̃ω−cos θ1)







sin  ̃ω cos(nω ̃) + (1 + cos  ̃ω) sin(nω ̃) − sin  ̃ω cos(nθ1) −

(1+cos θ1) sin  ̃ω

sin θ1

sin(nθ1)

2 sin  ̃ω cos(nω ̃) + 2 cos  ̃ω sin(nω ̃) − 2 sin  ̃ω cos(nθ1) −

2 cos θ1 sin  ̃ω

sin θ1

sin(nθ1)





 ,

(42)

by applying equations 25 and 30 to equation 41. As the driving frequency ω ̃ approaches any of

the values from the resonance frequency spectrum in equation 15, the denominator

(cos  ̃ω − cos θ1) of the coefficients of the stress solutions gets close to zero, allowing the stress

amplitude to take larger values.

17
Illustrations of the bounded stress solutions at non-resonance frequencies appear in figure 3. As

depicted in the plots in figure 3 and discussed previously in this section, when positioned in the

middle of the second layer, the time interval at which si(n) is reached is 4n+2i−3

4 < t < 4n+2i−14

for n ≥ 1 and i = 1, 2 (see figure 1). Beat phenomena are clearly visible in figure 3(b), as the

driving frequency approaches a resonance frequency.

0 10 20 30 40 50

−2.5

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

2.5

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

0 10 20 30 40 50 60

−15

−10

−5

0

5

10

15

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

(a) (b)

Figure 3. Stress time history for a two-layered Goupillaud-type strip in the middle

of the second layer located at ξ = 3/4. Impedance ratio α1 = 1/3,

angle θ1 = 2π/3, τ = 1, and loading f(n) = sin(nω ̃) at ξ = 0. (a)

ω ̃ = π/4. (b) ω ̃ = 0.9 · θ1 = 0.6π.

2.3.2 Stress Solutions at Resonance Frequencies

When cos  ̃ω = cos θ1 and sin  ̃ω = sin θ1, the stress solutions become







s1(n)

s2(n)





 =

1

χ1









−

(1+cos θ1)

sin θ1

n cos(nθ1) + hn +

1

1−cos θ1

i

sin(nθ1)

−

2 cos θ1

sin θ1

n cos(nθ1) + h2n +

2

sin2 θ1

i

sin(nθ1)









, (43)

18
by applying equations 25 and 38 to equation 41 for θk

∗ = θ1. Further simplifications of equation

43 implies







s1(n)

s2(n)





 = −

n

χ1









cos[(n+ 1

2 )θ1]

sin θ1

2

2 cos[(n+1)θ1]

sin θ1









+

sin(nθ1)

χ1







1

(1−cos θ1)

2

sin2 θ1







(44)

= −











cosh

(n+ 1

2 ) cos−1

3

χ1−2

χ1

 ́i

√χ1

cosh(n+1) cos−1

3

χ1−2

χ1

 ́i

√

χ1−1











· n +







1

2

χ1

2(χ1−1)





 · sin ·

n cos−1

μ

χ1 − 2

χ1

¶ ̧.

Given that 0 < θ1 < π, χ1 = (1 + α1) and cos θ1 =

χ1−1

χ1

, it follows that sin θ1 =

2

√

χ1−1

χ1

and

sin θ1

2 =

1

√χ1

. Therefore, all the stress solutions can be expressed in terms of the design

parameter χ1 or equivalently in terms of the impedance ratio α1 only, as shown in the last equality

of equation 44.

As seen in figure 1, the stress values alternate between f(n) and s1(n) in the first layer, and

between s1(n) and s2(n) in the second layer, for n = 1, 2, . . .. The terms of both stress

sequences s1(n) and s2(n) become unbounded over time due to the (time-related) linear factor n

(see equation 44). This implies that the stress amplitude will grow without bound in each layer,

confirming resonance. Based on these results, we expect to detect resonance in both, the first and

second layer, as illustrated in figure 4. The good agreement between the analytical results

generated from equation 38, or equivalently equation 44, and the stress values generated from the

recursive equations 2 for s1(n) and s2(n), n ≥ 1, is displayed in figure 5. The fact that the linear

factor n in equation 44 multiplies a cosine function explains the sinusoidal fluctuations of the

stress solutions in figures 4–5.

19
0 10 20 30 40 50

−30

−20

−10

0

10

20

30

40

50

Time

Stress Values in the First Layer

f(n)

s

1

(n)

0 10 20 30 40 50

−100

−50

0

50

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

(a) (b)

Figure 4. Stress time history for a two-layered Goupillaud-type strip with

impedance ratio α1 = 1/3, τ = 1, subjected to loading f(n) = sin(nω ̃)

with ω ̃ = θ1 = 2π/3. (a) Middle of the first layer located at ξ = 1/4.

(b) Middle of the second layer located at ξ = 3/4.

0 10 20 30 40 50

−30

−20

−10

0

10

20

30

40

50

n

s

1(n)

s

1

(n) Recursive

s

1

(n) Analytical

0 10 20 30 40 50

−100

−50

0

50

n

s

2(n)

s

2

(n) Recursive

s

2

(n) Analytical

(a) (b)

Figure 5. The values of the two stress sequences (a) s1(n), (b) s2(n), obtained

from the recursive equations 2 and the analytical equation 38, for the

two-layered Goupillaud-type strip with impedance ratio α1 = 1/3,

τ = 1, subjected to loading f(n) = sin(nω ̃) with ω ̃ = θ1 = 2π/3.

20
2.4 The Three-layer Case (m = 3)

For the three-layer case, it follows from equations 21–22 that ω ̃ =

2τ

3

ω and

f(n) = sin(n ·

2τ

3

ω) = sin(nω ̃). As a result, the first equation of 17 becomes

 ̃x3 =

z

2

sin( ̃ω)

(z

2 − 2z cos  ̃ω + 1) ·



















2(1+α2)z

2+4(1−α2)z+2(1+α2)

χ2(z+1)(z

2−2z cos θ1+1)

4z

χ2(z

2−2z cos θ1+1)

8z

2

χ2(z+1)(z

2−2z cos θ1+1)



















, (45)

where χ2 = (1 + α1)(1 + α2) > 1, cos θ0 = −1, and cos θ1 =

χ2−2

χ2

6= ±1 from equation D-2 in

appendix D. Applying in equation 45 partial fraction expansion as shown and equation 25,

results in

















s1(n)

s2(n)

s3(n)

















= Z

−1

( ̃x3) = Z

−1{

1

χ2

·

z

2

sin( ̃ω)

(z

2−2z cos  ̃ω+1) ·



















2

1+cos θ1

·

h

2α2

z+1 +

[(1+α2) cos θ1+(1−α2)](z+1)

z

2−2z cos θ1+1 i

4z

z

2−2z cos θ1+1

4

1+cos θ1

·

h

1

z+1 +

(1+2 cos θ1)z−1

z

2−2z cos θ1+1 i



















}.

(46)

21
2.4.1 Stress Solutions at Non-resonance Frequencies

When cos  ̃ω 6= −1 = cos θ0 and cos  ̃ω 6= cos θ1, bounded stress solutions are obtained by applying

equation 30 to equation 46 and further simplifying the results,

















s1(n)

s2(n)

s3(n)

















=

1

χ2

{



















−2α2 sin  ̃ω

(1+cos  ̃ω)(1+cos θ1)

0

−

2 sin  ̃ω

(1+cos  ̃ω)(1+cos θ1)



















(−1)n +



















[(1+α2) cos  ̃ω+(1−α2)] sin  ̃ω

(cos  ̃ω−cos θ1)(1+cos  ̃ω)

2 sin  ̃ω

(cos  ̃ω−cos θ1)

2(2 cos  ̃ω+1) sin  ̃ω

(cos  ̃ω−cos θ1)(1+cos  ̃ω)



















cos(nω ̃)

+



















(1+α2) cos  ̃ω+(1−α2)

(cos  ̃ω−cos θ1)

2 cos  ̃ω

(cos  ̃ω−cos θ1)

2(2 cos  ̃ω−1)

(cos  ̃ω−cos θ1)



















sin(nω ̃) −



















[(1+α2) cos θ1+(1−α2)] sin  ̃ω

(cos  ̃ω−cos θ1)(1+cos θ1)

2 sin  ̃ω

(cos  ̃ω−cos θ1)

2(2 cos θ1+1) sin  ̃ω

(cos  ̃ω−cos θ1)(1+cos θ1)



















cos(nθ1)

−



















[(1+α2) cos θ1+(1−α2)] sin  ̃ω

(cos  ̃ω−cos θ1) sin θ1

2 cos θ1 sin  ̃ω

(cos  ̃ω−cos θ1) sin θ1

2(2 cos θ1−1) sin  ̃ω

(cos  ̃ω−cos θ1) sin θ1



















sin(nθ1) }.

(47)

As the driving frequency ω ̃ approaches any of the values from the resonance frequency spectrum

in equation 15, the denominators (cos  ̃ω + 1) or (cos  ̃ω − cos θ1) of the coefficients of the stress

solutions in equation 47 get close to zero, allowing the stress amplitude to take larger values.

Illustrations of the bounded stress solutions at non-resonance frequencies are given in figure 6.

As depicted in figure 6 and discussed previously in this section, when positioned in the middle of

the second or third layer, the time interval at which si(n) is reached is 4n+2i−3

6 < t < 4n+2i−16

, as

applicable, for n ≥ 1 and i = 1, 2, 3 (see figure 1). The good agreement between the analytical

results generated from equation 30, or equivalently equation 47, and the stress values generated

from the recursive equations 2 for s1(n), n ≥ 1, as well as the error introduced by the

computational approximation of quantities (such as sine and cosine values), are displayed in

figure 7.

22
0 10 20 30 40 50

−3

−2

−1

0

1

2

3

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

0 10 20 30 40 50

−4

−3

−2

−1

0

1

2

3

4

Time

Stress Values in the Third Layer

s

2

(n)

s

3

(n)

(a) (b)

Figure 6. Stress time history for a three-layered Goupillaud-type strip with

α1 = 3, α2 =

1

2(1−cos(π/4)) − 1, θ1 =

π

4

, τ = 1, subjected to loading

f(n) = sin(nω ̃) with ω ̃ =

π

3

. (a) Middle of the second layer located at

ξ = 1/2. (b) Middle of the third layer located at ξ = 5/6.

0 10 20 30 40 50 60 70 80 90 100

−2.5

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

2.5

n

s

1(n)

s

1

(n) Recursive

s

1

(n) Analytical

0 20 40 60 80 100

−2

−1

0

1

2

3

x 10−14

n

(Recursive − Analytical) for s

1(n)

(a) (b)

Figure 7. The three-layered Goupillaud-type strip with α1 = 3,

α2 =

1

2(1−cos(π/4)) − 1, θ1 =

π

4

, τ = 1, subjected to loading

f(n) = sin(nω ̃) with ω ̃ =

π

3

. (a) The values of the stress sequence

s1(n), obtained from the recursive equations 2 and the analytical

equation 30. (b) The difference of the expressions for s1(n) obtained

from the recursive equations 2 and the analytical equation 30.

23
2.4.2 Stress Solutions at Resonance Frequencies

2.4.2.1 Case I. cos  ̃ω = −1 = cos θ0

The stress formulas in the z-space for the three-layer case in equation 45, indicate that the

sequences s1(n) and s3(n) increase without bound because cos( ̃ω) = −1 = cos θ0 corresponds to

a zero of the denominator z = e

Iθ0 = −1. Alternatively, we do not expect unbounded growth for

the sequence s2(n), since z = e

Iθ0 = −1 is not a zero of the denominator. Figure 8 illustrates

this behavior for each of the stress sequences. Similar behavior was observed in the numerical

experiments for the seven- and eleven-layer case (see illustration in figure 2) where the stress

sequences with odd subindices increase without bound while the stress sequences with even

subindices are bounded.

0 5 10 15 20

−8

−6

−4

−2

0

2

4

6

8

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

0 5 10 15 20

−10

−5

0

5

10

15

Time

Stress Values in the Third Layer

s

2

(n)

s

3

(n)

(a) (b)

Figure 8. Stress time history for a three-layered Goupillaud-type strip with

α1 = 3, α2 =

1

2(1−cos(π/4)) − 1, τ = 1, subjected to loading

f(n) = cos(nω ̃) with ω ̃ = π. (a) Middle of the second layer located at

ξ = 1/2. (b) Middle of the third layer located at ξ = 5/6.

2.4.2.2 Case II.a.

When cos  ̃ω = cos θ1 and sin  ̃ω = sin θ1, and by applying equation 38 to equation 46 the stress

solutions become,

24
















s1(n)

s2(n)

s3(n)

















= −

2 sin θ1

χ2(1+cos θ1)

2

















α2

0

1

















(−1)n +

2

χ2

















α2 sin θ1

(1+cos θ1)

2 −

[(1+α2) cos θ1+(1−α2)]

2 sin θ1

n

−

cos θ1

sin θ1

n

sin θ1

(1+cos θ1)

2 −

(2 cos θ1−1)

sin θ1

n

















cos(nθ1)

+

2

χ2

















(1−α2) cos θ1+(1+α2)

2 sin2 θ1

+

[(1+α2) cos θ1+(1−α2)]

2(1+cos θ1)

n

n +

1

sin2 θ1

2−cos θ1

sin2 θ1

+

(2 cos θ1+1)

1+cos θ1

n

















sin(nθ1).

(48)

Since the focus here is the study of resonance, the unbounded part of each of the stress solutions

containing the time-related factor n can be simplified as follows:

2n

χ2

















−

[(1+α2) cos θ1+(1−α2)]

2 sin θ1

cos(nθ1) + [(1+α2) cos θ1+(1−α2)]

2(1+cos θ1)

sin(nθ1)

−

cos θ1

sin θ1

cos(nθ1) + sin(nθ1)

−

(2 cos θ1−1)

sin θ1

cos(nθ1) + (2 cos θ1+1)

1+cos θ1

sin(nθ1)

















=

−

2n

χ2





















[(1+α2) cos θ1+(1−α2)]

2 cos θ1

2

sin θ1

cos £¡n +

1

2

¢

θ1

¤

cos[(n+1)θ1]

sin θ1

cos[(n+ 3

2 )θ1]

cos

θ1

2

sin θ1





















= −n





















2α1

√χ2

(1+α1)(χ2−1) cos h¡

n +

1

2

¢

cos−1 χ2−2

χ2

i

cosh(n+1) cos−1 χ2−2

χ2

i

√

χ2−1

√χ2

(χ2−1) · cos h¡

n +

3

2

¢

cos−1 χ2−2

χ2

i





















(49)

Given that 0 < θ1 < π, χ2 = (1 + α1)(1 + α2), and cos θ1 =

χ2−2

χ2

, it follows that

sin θ1 =

2

√

χ2−1

χ2

. Therefore, the stress solutions can be expressed in terms of the impedance

ratios α1 and α2 only. This is demonstrated in the last equality (equation 49).

25
As seen in figure 1, the stress values alternate between f(n) and s1(n) in the first layer, between

s1(n) and s2(n) in the second layer, and between s2(n) and s3(n) in the third layer, for

n = 1, 2, . . .. All three stress sequences/solutions s1(n), s2(n), and s3(n) become unbounded

over time due to the (time-related) linear factor n. This implies that the stress amplitude will

grow without bound in all three layers (figure 9). The good agreement between the analytical

results generated from equation 38, or equivalently equation 48, and the stress values generated

from the recursive equations 2 for s1(n), n ≥ 1, as well as the error introduced by the

computational approximation of quantities (such as sine and cosine values) are displayed in figure

10. The fact that the linear factor n multiplies a cosine function as shown in equation 49,

explains the sinusoidal fluctuations of the stress solutions in figures 9–10.

0 10 20 30 40 50

−30

−20

−10

0

10

20

30

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

0 10 20 30 40 50

−30

−20

−10

0

10

20

30

40

Time

Stress Values in the Third Layer

s

2

(n)

s

3

(n)

(a) (b)

Figure 9. Stress time history for a three-layered Goupillaud-type strip with

α1 = 3, α2 =

1

2(1−cos(π/4)) − 1, τ = 1, subjected to loading

f(n) = sin(nω ̃) with ω ̃ = θ1 =

π

4

. (a) Middle of the second layer

located at ξ = 1/2. (b) Middle of the third layer located at ξ = 5/6.

26
0 20 40 60 80 100

−30

−20

−10

0

10

20

30

40

n

s

1(n)

s

1

(n) Recursive

s

1

(n) Analytical

0 20 40 60 80 100

−3

−2

−1

0

1

2

3

4

x 10−13

n

(Recursive − Analytical) for s

1(n)

(a) (b)

Figure 10. (a) The values of the stress sequence s1(n), obtained from the

recursive equations 2 and the analytical equation 38. (b) The

difference of the expressions for s1(n) obtained from the recursive

equations 2 and the analytical equation 38 for a three-layered

Goupillaud-type strip with α1 = 3, α2 =

2

2(1−cos(π/4)) − 1, τ = 1,

subjected to loading f(n) = sin(nω ̃) with ω ̃ = θ1 =

π

4

.

3. Analytical Formulas for the Natural Frequency Spectrum of a

Free-fixed Goupillaud-type Layered Elastic Strip

In this section, we derive the natural frequencies of a free-fixed m-layered Goupillaud-type strip

in two different ways:

1. By simplifying and then solving the frequency equation for 2 ≤ m ≤ 5.

2. By extending the results of the resonance frequency spectrum from the discrete model to

the continuous model for any m ≥ 2. In the continuous model, the strip is subjected to a

continuous harmonic forcing function at one end and held fixed at the other end.

3.1 The Frequency Equation: Explicit Formula for the Natural Frequency Spectrum

The frequency equation for a free-fixed m-layered Goupillaud-type strip is derived in this section,

by generalizing the procedure described in Graff (33). The frequency equation is simplified

through a transformation of the spatial variable, to derive the natural frequency spectrum for

2 ≤ m ≤ 5.

27
The governing equations for the displacement u(x, t) in each layer of the strip are:

∂

2ui(x,t)

∂t2 = c

2

i

∂

2ui(x,t)

∂x2 for xi−1 < x < xi

, i = 1, . . . , m, (50)

where 0 = x0 < x1 < · · · xi < · · · < xm = L, L represents the length of the strip, ci represents

the wave speed in the ith layer, x represents the spatial variable, and t ≥ 0 represents the time

variable. Seeking harmonic solutions, we may represent u(x, t) = U(x)ψ(t) and the spatial

component(s) of the displacement for each layer as

Ui(x) = Ai sin βix + Bi cos βix for xi−1 < x < xi, i = 1, . . . , m, (51)

where βi = ω/ci for i = 1, . . . , m, and ω represents the frequency of the harmonic motion.

Similarly, ψ(t) can be represented in the form ψ(t) = A∗sin ωt + B∗cos ωt or in exponential

form as ψ(t) = A∗∗e

Iωt + B∗∗e−Iωt, where I =

√

−1. As a result, after applying the free-fixed

boundary conditions at x = 0 and x = L, respectively,

dU1(0)

dx = Um(L) = 0, (52)

and the continuity of stress and displacement at each layer interface, results in the following linear

system of equations:







A1 = 0,

Ai sin 3

ω

ci

xi

 ́

+ Bi cos 3

ω

ci

xi

 ́

= Ai+1 sin 3

ω

ci+1

xi

 ́

+ Bi+1 cos 3

ω

ci+1

xi

 ́

for i = 1, . . . , m − 1,

αi

h

Ai cos 3

ω

ci

xi

 ́

− Bi sin 3

ω

ci

xi

 ́i = Ai+1 cos 3

ω

ci+1

xi

 ́

− Bi+1 sin 3

ω

ci+1

xi

 ́

for i = 1, . . . , m − 1,

Am sin 3

ω

cm

xm

 ́

+ Bm cos 3

ω

cm

xm

 ́

= 0,

(53)

The relations αi =

ciρi

ci+1ρi+1

=

Ei

ci

·

ci+1

Ei+1

were applied equation 53, using the fact that c · ρ = E/c.

The condition that the determinant of the system matrix for equation 53 must vanish generates the

frequency equation for ω. Solving the frequency equation and obtaining a general formula could

be a challenging task. However, if the spatial variable x of the wave equation is replaced with the

new variable ξ =

R x

0

ds

c(s)

, the frequency equation can be simplified. This transformation converts

the problem to a Goupillaud-type layered strip with the same length L as the wave transit time τ

through the strip (L = τ ) (see Velo and Gazonas (40) and section 2). Under this transformation,

the layers have equal lengths of τ

m

and wave speed of unity. In addition, the previous coordinates

x0 = 0, xj, and xm = L are replaced by the new coordinates ξ0 = 0, ξj =

jτ

m

, and ξm = τ for

28
j = 1, 2, . . . , m, while the impedance ratio αi =

E ̃

i

E ̃

i+1

=

√

√

Eiρi

Ei+1ρi+1

remains the same for

i = 1, 2, . . . , m − 1. As before, Ei, ρi represent the actual material properties for the physical

case, while E ̃

i

, ρ ̃i represent the material properties of the ith layer for the simplified case, with

E ̃

i =  ̃ρi =

√

Eiρi. As seen next, both designs of the strip share the same natural frequency

spectrum under free-fixed boundary conditions.

As a result, equation 53 can be transformed in a new matrix-vector form as follows:

∆m ̃v =  ̃0, (54)

where

∆m =





















































1 0 0 0 . . . 0

sin γ1

α1 cos γ1

0

0

0

.

.

. ∆ ̃ m

0





















































2m×2m

,  ̃v =













































A1

B1

.

.

.

Ai

Bi

.

.

.

Am

Bm













































2m×1

, (55)

and

∆ ̃ m =

=



















































cos γ1 − sin γ1 − cos γ1 0 0 0 0

−α1 sin γ1 − cos γ1 sin γ1 0

0 sin γ2 cos γ2 − sin γ2 − cos γ2

0 α2 cos γ2 −α2 sin γ2 − cos γ2 sin γ2

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

0 . . . 0 0

0 . . . 0 sin γm−1 cos γm−1 − sin γm−1 − cos γm−1

0 . . . 0 αm−1 cos γm−1 −αm−1 sin γm−1 − cos γm−1 sin γm−1

0 . . . 0 0 0 sin γm cos γm



















































.

(56)

29
Here γi = ξiω =

iτω

m

for i = 1, . . . , m. Seeking a nonzero solution for equation 54, the

determinant of the system matrix ∆m must vanish:

|∆m| = 0. (57)

Equation 57 gives the transformed/simplified frequency equation for ω. Equation 57 can be

successfully solved, using symbolic algebra software, for a strip with up to five layers (see

appendix F). Solving the frequency equation 57 for the general m-layer case, poses a

computational challenge.

Based on these derivations, a comparison of equation 14 with equation 57, equation D-1 with

equation F-1, equation D-2 with equation F-3, equation D-4 with equation F-5, and equation D-6

with equation F-9 reveals that the two matrices Dm and ∆m, defined in equation C-1 and

equations 55–56, have the same dimensions and cosine roots of their determinants subject to the

transformation

cos θ = cos

2τ

m

ω. (58)

Also notice that the pairs of the frequency equations F-1, F-3 and F-5, F-9, have a similar pattern

to equations D-1, D-2 and D-4, D-6, respectively. Equations 58 can be justified from equations

21–22 and the fact that the natural frequencies of a free-fixed elastic strip correspond to the

resonance frequencies when a (continuous) forcing, which varies harmonically with time is

applied at the free end of the strip.

3.2 Heuristic Approach: Implicit General Formula for the Natural Frequency Spectrum

Using the linear frequency transformation in equation 21, the resonance frequency spectrum for

the continuous model (where the strip is subjected to a continuous harmonic forcing at one end

and held fixed at the other end) can be recovered from the resonance frequency results in equation

15 for the discrete model. This results in the following implicit formula for the natural frequency

spectrum of a free-fixed m-layered Goupillaud-type strip:







ωl0 =

m

2τ

· [θ0 + 2l0π] = 1

2(

τ

m )

· (2l0 + 1)π, (for m odd only),

ωl1,k =

m

2τ

· [θk + 2l1π] = 1

2(

τ

m )

· [θk + 2l1π], l0, l1, l2 = 0, 1, 2, . . . ,

ωl2,k =

m

2τ

· [−θk + 2(l2 + 1)π] = 1

2(

τ

m )

· [−θk + 2(l2 + 1)π], k = 1, . . . , b

m

2

c.

(59)

30
Substitution of the angle values from equations D-1, D-3, D-5, and D-7 into equation 59, results

in the natural frequency spectrum for the free-fixed Goupillaud-type layered strip with up to five

layers, equations F-2, F-4, and F-6, derived independently in appendix F from the frequency

equation 57. Although a rigorous proof for the equation 59 is not provided, our theoretical and

numerical results validate it as the general natural frequency spectrum of a free-fixed m-layered

Goupillaud-type strip with the zeros of the corresponding matrix |Am| on the unit circle (see

matrix Am defined in the z-space in equation 7). The zeros of |Am| on the unit circle generate

the respective angles θ0 = π (for m odd only) and {±θk}

b m

2

c

k=1 needed to recover the natural

frequency spectrum.

3.3 Interpretation and Illustration of the Natural Frequency Results

The explicit natural frequency equations (F-2, F-4, and F-6) derived from the frequency equation

57, as well as the general frequency spectrum derived from the implicit equation 59, suggest that

for a free-fixed m-layered Goupillaud-type strip, the natural frequency spectrum does not depend

on the actual length of the strip. According to these formulas, the natural frequency spectrum

depends only on certain combinations of the impedance ratios α1, α2, . . . , αm−1 and (it is

inversely proportional to) the (equal) wave travel time τ

m

for each layer.

The following illustration validates these expectations and demonstrates the simplicity of

predicting the natural frequency spectrum when using the general equation 59, or solving the

simplified frequency equation 57, versus solving directly the frequency equation generated from

equation 53. Our numerical experiments with up to five layers have shown that the natural

frequency spectrum can be predicted using either equation 59 and appendix D, or appendix F,

even when a symbolic algebra software such as Maple cannot solve the frequency equation

obtained from equation 53.

3.3.1 Case a. Two-layered Goupillaud-type Elastic Strip with Unequal Layer Lengths

• Layer lengths: L1 =

√

2, L2 = 3√2

• Length of the strip: L = L1 + L2 = 4√2

• Wave speed per layer: c1 = 1, c2 = 3

• Elastic modulus per layer: E1 = 2 + √3, E2 = 3(2 −

√

3)

• Wave travel time per layer: L1

c1

=

L2

c2

=

τ

m =

√

2

• Transit time through the strip: τ = 2√2

• Impedance ratio: α1 =

E1

c1

·

c2

E2

=

2+√3

2−

√

3

31
• Frequency equation from equation 53:

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

1 0 0 0

sin ¡√2 ω

¢

cos ¡√2 ω

¢

− sin 3 √

2

3

ω

 ́

− cos 3 √

2

3

ω

 ́

α1 cos ¡√2 ω

¢

−α1 sin ¡√2 ω

¢

− cos 3 √

2

3

ω

 ́

sin 3 √

2

3

ω

 ́

0 0 sin 3

4

√

2

3

ω

 ́

cos 3

4

√

2

3

ω

 ́

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

= 0. (60)

3.3.2 Case b. Simplified Two-layered Goupillaud-type Elastic Strip with Equal Layer

Lengths

This free-fixed strip has the same natural frequency spectrum as the one in Case a.

• Layer lengths: L ̃

1 =

L1

c1

= L ̃

2 =

L2

c2

=

τ

2 =

√

2

• Length of the strip: L ̃ = L ̃

1 + L ̃2 = τ = 2√

2

• Wave speed per layer: c ̃1 =  ̃c2 = 1

• Elastic modulus per layer: E ̃

1 = 2 + √

3, E ̃

2 = 2 −

√

3

• Wave travel time per layer: L ̃1

c ̃1

=

L ̃2

c ̃2

=

τ

2 =

√

2

• Transit time through the strip: τ = 2√2

• Impedance ratio: α ̃1 =

E ̃1

E ̃2

=

2+√3

2−

√

3

= α1

• Frequency equation from equation 57:

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

1 0 0 0

sin ¡√

2 ω

¢

cos ¡√2 ω

¢

− sin ¡√2 ω

¢

− cos ¡√2 ω

¢

α1 cos ¡√2 ω

¢

−α1 sin ¡√2 ω

¢

− cos ¡√2 ω

¢

sin ¡√2 ω

¢

0 0 sin ¡2

√

2 ω

¢

cos ¡2

√

2 ω

¢

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

 ̄

= 0. (61)

3.3.3 The Predicted Natural Frequency Spectrum

The common natural frequency spectrum for the free-fixed strip, for Cases a. and b.:







ωl1 =

1

2

√

2

¡

π

6 + 2l1π

¢

=

π

12√2

+ √π

2

l1,

l1, l2 = 0, 1, 2, . . . ,

ωl2 =

1

2

√

2

¡

−

π

6 + 2(l2 + 1)π

¢

=

11π

12√2

+ √π

2

l2,

(62)

32
can be derived easily and accurately from equation 59 and equation D-1 for

θ1 = cos−1

3

α1−1

α1+1 ́

=

π

6

, or equation F-2. The Maple software recovers most of these results by

finding the zeros of the determinants previously given in Cases a. and b. As expected, the

frequency equation in the (simplified) Case b. is much easier to solve in Maple than the original

frequency equation in Case a.

3.3.4 Case c. Design Modification That Gives a Desired Frequency Spectrum Within

Limitations

One way to modify the natural frequency spectrum equation 62 is to modify any of the two

layered strips given in Cases a. and b., with χ1 = (1 + α1) and θ1 =

π

6

, by adding a new layer in

front of the existing two layers. The length of the new layer must be chosen so that the strip

remains of Goupillaud-type. The new three-layered strip is characterized by

χ2 = (1 + α

∗

1

)(1 + α

∗

2

) with α

∗

1

to be determined and α

∗

2 = α1. For a desirable value of the

cosine of the base angle θ

∗

1 given by V = cos θ

∗

1 =

χ2−2

χ2

, one can derive the necessary formula for

the unknown impedance ratio α

∗

1

:

α

∗

1 = −1 +

2

χ1(1 − V )

, where χ1 − 2

χ1

= cos θ1 < V = cos θ

∗

1 < 1, (63)

using the fact that χ2 = (1 + α

∗

1

)χ1. For instance, for the choice of θ

∗

1 =

π

9 we have that

cos θ1 = cos ¡

π

6

¢

< V = cos θ

∗

1 = cos ¡

π

9

¢

, which results in the natural frequency spectrum:







ωl0 =

1

2

√

2

(π + 2l0π) = π

2

√

2

+ √π

2

l0,

ωl1 =

1

2

√

2

¡

π

9 + 2l1π

¢

=

π

18√2

+ √π

2

l1, l0, l1, l2 = 0, 1, 2, . . . ,

ωl2 =

1

2

√

2

¡

−

π

9 + 2(l2 + 1)π

¢

=

17π

18√2

+ √π

2

l2,

(64)

for the free-fixed three-layered strip with impedance ratios α

∗

1 = −1 + 2−

√

3

2(1−cos(

π

9 )), α

∗

2 =

2+√3

2−

√

3

and the (equal) wave travel time per layer τ

∗

3 =

τ

2 =

√

2.

33
4. Other Applications of the Frequency Results

The properties of the materials used in this section appear in appendix G. Using the definition of

the wave speed c =

p

E/ρ, the characteristic impedance ρc in each layer is calculated in terms of

the material properties as √Eρ. Once the length of one layer is chosen, the length of each

remaining layer is determined according to the material properties to provide equal wave travel

time: Li

ci

=

Li+1

ci+1

or Li+1 =

qEi+1

ρi+1

·

ρi

Ei

Li. Here Lj, cj, Ej and ρj, for j = i, i + 1, represent the

length, wave speed, elastic modulus, and material density for the ith and (i + 1)th layer,

respectively, i = 1, . . . , m − 1.

4.1 One-dimensional Layered Media with a Common Frequency Spectrum

According to equations 14 and 15, the resonance frequency spectrum for the discrete model

depends only on certain combinations of the impedance ratios α1, α2, . . . , αm−1. Such

combinations of impedance ratios when 2 ≤ m ≤ 5 are represented by the design parameters

χm−1 for m = 2, 3, and χm−1, Γm−1 for m = 4, 5. There are (theoretically) infinitely many

m-layered designs of a Goupillaud-type elastic strip, which have the same values of the design

parameters χm−1, Γm−1 for 2 ≤ m ≤ 5, hence the same resonance frequency spectrum for the

discrete model.

As for the natural frequency spectrum of a free-fixed m-layered Goupillaud-type strip based on

equation 59, it depends not only on combinations of the impedance ratios α1, α2, . . . , αm−1 but

also on the (equal) wave travel time τ

m

for each layer. Therefore, to obtain designs with a

common natural frequency spectrum, besides having the same values of the relevant design

parameters, the layer lengths have to be chosen so that they all have the same wave travel time.

Another way to find designs with a common natural frequency spectrum was previously

illustrated in subsection 3.3 using a transformation of the spatial variable.

For the three-layer case where χ2 = (1 + α1)(1 + α2), the tungsten-lead-aluminum configuration

with α1 ≈ 5.828 and α2 ≈ 0.952 has the same value of χ2 ≈ 13.33 as the lead-aluminum-copper

configuration with α1 ≈ 0.952 and α2 ≈ 5.828. Another material design with χ2 ≈ 13.33 is the

aluminum-steel-brass configuration with α1 ≈ 0.331 and α2 ≈ 9.016. In addition, in order for

these material designs to have a common natural frequency spectrum for the free-fixed strip, only

the length of one layer in one of the designs may be chosen at random. The lengths of the

remaining layers in all the three designs are then determined to provide equal wave travel time.

34
For the four-layer case, we have that χ3 = (1 + α1)(1 + α2)(1 + α3) and Γ3 = α1α3 − 1. The

aluminum-(glass fiber reinforced material)-nickel-bronze configuration with α1 ≈ 2.27,

α2 ≈ 0.698 and α3 = 1/α1 ≈ .44, has almost the same value of χ3 = 8 and Γ3 = 0, as the

homogeneous strip (all “four layers” are made of the same elastic material). As a result, these

designs share the same resonance frequency spectrum when subjected to a discrete harmonic

forcing at one end and held fixed at the other end. The common frequency spectrum can be

derived from equation 15 after substituting χ3 = 8 and Γ3 = 0 into equation D-5. In addition, a

common wave travel time for each of the four layers in both designs would ensure the same

natural frequency spectrum under free-fixed boundary conditions.

Furthermore, according to equation D-5 and table E-1, these designs with parameters χ3 = 8 and

Γ3 = 0, and base angle values θ1,opt = π/4 and θ2,opt = 3π/4 obtained for j = 2, are expected to

be optimal, i.e., provide a stress amplitude of double the loading when a Heaviside step in stress

loading is applied at one end of the strip while the other end is held fixed.

4.2 Resonance Frequencies and Optimal Designs

For the m-layered strip subjected to a Heaviside loading (40), it was shown that when the optimal

values of the base angles θk,opt for k = 1, . . . , b

m

2

c, given in appendix E, were substituted into the

respective angle-design parameter equations in appendix D, optimal m-layered designs were

generated for 2 ≤ m ≤ 5. These optimal designs provide the smallest stress amplitude of double

the loading in all layers, for all times.

However, when these optimal angle values are substituted into equations 15 and 59, they generate

the resonance frequency spectrum for the discrete and continuous model, respectively. This

means that for the optimal designs identified in appendices D–E, the displacement and stress

amplitude will grow without bound over time when a harmonic forcing with a resonance

frequency is applied instead of the Heaviside loading. As a result, for these designs, depending

on the type of loading, one can get either the best or the worst outcome when it comes to

controlling the stress amplitude (see the illustration in figure 11). As previously seen in figure 1,

the stress values in a two-layered strip are given by the sequences f(n), s1(n), and s2(n), n ≥ 0.

For the example considered in figure 11(a), we have that |f(n)| ≤ 1 for all n ≥ 0, therefore the

stress amplitude, which exceeds unity, is determined by the other two stress sequences.

Similarly, for the example considered in figure 11(b) we have that f(n) = 1 for all n ≥ 0,

therefore the stress amplitude, which exceeds unity, is also determined by the other two stress

sequences. Notice that as the stress amplitude in figure 11(a) grows without bound over time for

a harmonic forcing condition, the stress amplitude in figure 11(b) does not exceed double the

35
loading for a Heaviside loading condition. Examples of such material designs were given in the

previous subsection for the four-layer case.

0 10 20 30 40 50

−30

−20

−10

0

10

20

30

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

0 10 20 30 40 50

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Time

Stress Values in the Second Layer

s

1

(n)

s

2

(n)

(a) (b)

Figure 11. Stress time history for a two-layered Goupillaud-type strip with

impedance ratio α1 = 3, θ1 =

π

3

, τ = 1, at the middle of the second

layer located at ξ = 3/4 subjected to loading (a) f(n) = sin(nω ̃) with

ω ̃ = θ =

π

3

. (b) f(n) = 1 for n ≥ 0.

4.3 Natural Frequencies of a Free-fixed Non-Goupillaud-type Layered Strip with Integer

Layer Length Ratios

The natural frequency results for a free-fixed Goupillaud-type layered strip with equal layer

lengths given in equations F-1 through F-9, may be extended to a free-fixed non-Goupillaud-type

layered strip with integer layer length ratios. We illustrate this with the two examples included

below.

A three-layered Goupillaud-type strip with equal layer lengths (α2 = 1 or α1 = 1), represents a

two-layered non-Goupillaud-type strip with 1:2 or 2:1 wave travel time ratios respectively. The

case with 1:2 ratio requires that α2 = 1 and equation F-3 becomes

cos

2τ

3

ω = −1 or cos

2τ

3

ω =

α1

1 + α1

, for α1 > 0. (65)

Similar results to equation 65 can be obtained for the case with 2:1 ratio which requires that

α1 = 1.

The frequency results for a free-fixed four-layered Goupillaud-type strip with equal layer lengths

given in equation F-5, can be extended to a free-fixed three-layered non-Goupillaud-type strip

with wave travel time ratios 1:2:1 by choosing α2 = 1. With these conditions, equation F-5

36
becomes

(1 + α1)(1 + α3) cos2

τ

2

ω − Γ3 cos

τ

2

ω + (Γ3 − (1 + α1)(1 + α3) + 2) = 0, (66)

where Γ3 = (α1α3 − 1) is chosen so that the (cosine) solutions vary in the interval [−1, 1]. This

can be achieved when Γ3 = 0, which implies that (1 + α1)(1 + α3) > 2. As a result, equation 66

reduces to:

0 < cos2

3τ

2

ω

 ́

=

(1 + α1)(1 + α3) − 2

(1 + α1)(1 + α3)

< 1. (67)

4.4 Representation of the Stress Solutions in Terms of the Chebyshev Polynomials

The formulas related to the stress solutions in equations 30, 35, and 38 consist of sums with terms

of the form cos(nφ), sin(nφ), (−1)n, n and (−1)n(n + 1). These terms relate to the Chebyshev

polynomials of the first and second kind Tn(y) and Wn(y), respectively (see Schaum (44)), as

shown below:

cos(nφ) = Tn(cos φ), (−1)n = Tn(−1), n = Wn(1) − 1, (−1)n(n + 1) = Wn(−1). (68)

In addition, from the definition of the Chebyshev polynomials when cos φ 6= 0, one can also

derive that

sin(nφ) = [ Wn(cos φ) − Tn(cos φ) ] · tan φ. (69)

As a result, the stress formulas given by equations 30, 35, and 38 can be written in terms of

Chebyshev polynomials. We demonstrate this for the general stress representation in equation

20, which can be written in terms of the Chebyshev polynomials as shown below:

si(n) = f(n) ∗ [ bi,0Tn(−1) +

b m

2P

c

k=1

(ai,k − bi,k tan θk) · Tn(cos θk) + bi,k tan θk · Wn(cos θk) ]−

−f(0) · [ bi,0Tn(−1) +

b m

2P

c

k=1

(ai,k − bi,k tan θk) · Tn(cos θk) + bi,k tan θk · Wn(cos θk) ],

(70)

for i = 1, 2, . . . m and n ≥ 0. This representation is valid when cos θk 6= 0 for all 1 ≤ k ≤ bm

2

c.

37
5. Conclusions

The resonance frequency spectrum in equation 15 of the discrete model was derived by

developing an analytical representation for the stress solutions when a discrete sinusoidal forcing

function was applied at one end of the strip while the other end was held fixed. The approach

was based on the z-transform methods. Two key elements of this approach were that the

determinant of the resulting global matrix |Am| in the z-space was a palindromic polynomial with

real coefficients and that its zeros were distinct and on the unit circle. An important connection

was made in equation 15 between the angles (θ0 = π and {±θk}

b m

2

c

k=1 ), representing the distinct

zeros of |Am| and the natural/resonance frequency spectrum. We showed that the resonance

frequency spectrum in equation 15 for the discrete model consists of all the positive angles

coterminal with θ0 = π, for m, odd only, and {±θk}

b m

2

c

k=1 , for any m. It is important to notice that

the positive angles coterminal with θ0 = π represent resonance frequencies for any strip with an

odd number of layers and do not depend on the material properties. According to equation 15,

the resonance frequency spectrum for the discrete model depends only on (certain combinations

of) the impedance ratios and not on any other parameters such as the length of the strip, etc. The

design parameters representing such combinations for an m-layered strip were

χm−1 =

mQ−1

i=1

(1 + αi) for 2 ≤ m ≤ 5 and Γm−1 for m = 4, 5. This means that as long as

m-layered designs with different impedance ratios have such design parameters in common, they

are guaranteed to have the same resonance frequency spectrum for the discrete model.

The natural frequency spectrum of a free-fixed Goupillaud-type strip with up to five layers was

first derived by simplifying and then solving the frequency equation. Independently of that, the

linear transformation in frequency in equation 21 enabled us to heuristically describe the natural

frequency spectrum of a free-fixed m-layered Goupillaud-type strip from the discrete model in

equation 15. Both approaches provided an easy way to analytically describe the natural

frequency spectrum versus solving the frequency equation directly. The frequency spectrum in

equations 59 and appendix D, or appendix F, suggests that the natural frequency spectrum of a

free-fixed m-layered Goupillaud-type strip depends only on (combinations of) the impedance

ratios and (it is inversely proportional to) the equal wave travel time for each layer τ

m

. Designs

with these quantities in common share the same natural frequency spectrum for the free-fixed

boundary conditions (see subsections 3.3.1 and 3.3.2) and subsection 4.1 for illustrations of two

different ways of obtaining such designs). A possible design modification that provides a desired

natural frequency spectrum is also demonstrated in subsection 3.3.4.

38
When the resonance frequency results are combined with the stress optimization results obtained

in Velo and Gazonas (40), we conclude that for a given optimal design of a Goupillaud-type

layered strip with one fixed end, depending on the type of loading at the other end, one can get

either the best or worst outcome when it comes to the stress amplitude. Although this work

focused on a Goupillaud-type layered elastic strip, the frequency results were shown to extend to

a non-Goupillaud-type layered strip. The stress formulas involving trigonometric sums were also

interpreted using the Chebyshev polynomials of the first and second kind.

In this work, we have proven that the palindromic polynomial with real coefficients |Am| has all

its zeros on the unit circle for 1 ≤ m ≤ 5, see appendix B. Whether this is true for any

m-layered Goupillaud-type strip or only for certain classes of designs remains to be investigated.

Proving rigorously/theoretically the conversion of the resonance frequency spectrum from the

discrete model in equation 15 to the continuous model in equation 59 also remains to be shown.

Finally, the experimental validation of the frequency results proposed here would be a natural and

valuable extension of this work.

39
INTENTIONALLY LEFT BLANK.

40
6. References

1. Poisson, S. D. Mémoire sur l’équilibre et le Mouvement des Corps élastiques, Paris. Paris,

Mém. de l’Acad. 1829, 8 357–570.

2. Lamb, H. On the Vibrations of an Elastic Sphere. P. Lond. Math. Soc. 1882, 13 189–212.

3. Love, A. E. H. A Treatise on the Mathematical Theory of Elasticity; 4

th ed.; Dover

Publications: New York, 2002 (originally Cambridge University Press: 1927).

4. Ben-Menahem, A.; Singh, S. J. Seismic Waves and Sources; Springer-Verlag: New York,

1981.

5. Eringen, A. C.; Suhubi, E. S. Elastodynamics, Vol. II, Linear Theory; Academic Press:

New York, 1975.

6. Nishida, K.; Kobayashi, N.; Fukao, Y. Resonant Oscillations Between Solid Earth and the

Atmosphere. Science 2000, 287 2244–2246.

7. Brekhovskikh, L. M. Waves in Layered Media; Academic Press: New York, 1960.

8. Visscher, W. M.; Migliori, A.; Bell, T. M.; Reinert, R. A. On the Normal Modes of Free

Vibration of Inhomogeneous and Anisotropic Elastic Objects. J. Acous. Soc. Am. 1991, 90

(4) 2154–2162.

9. Kosevich, A. M.; Minaev, P. A.; Syrkin, E. S. Interaction of Acoustic Waves with an

Interface in Highly Anisotropic Layered Crystals. Low Temp. Phys. 2004, 30 (6) 500–504.

10. Ouchani, N.; Nougaoui, A.; Aynaou, H.; Bria, D.; El Boudouti, E. H.; Velasco, V. R.;

Daoudi, A. Acoustic Waves in (110) Layered Structures. Surf. Sci. 2008, 602 2107–2113.

11. Fiorito, R.; Madigosky, W.; Uberall, H. Resonance Theory of Acoustic Waves Interacting

with an Elastic Plate. J. Acous. Soc. Am. 1979, 66 (6) 1857–1866.

12. Jackins, P. D.; Gaunaurd, G. C. Resonance Acoustic Scattering from Stacks of Bonded

Elastic Plates. J. Acous. Soc. Am. 1986, 80 (6) 1762–1776.

13. Mead, D. J. The Forced Vibration of One-Dimensional Multi-Coupled Periodic Structures:

An Application to Finite Element Analysis. J. Sound Vibrat. 2009, 319 282–304.

41
14. Rao, M. K.; Scherbatiuk, K.; Desai, Y. M.; Shah, A. H. Natural Vibrations of Laminated

and Sandwich Plates. J. Eng. Mech. 2004, 130 (11) 1268–1278.

15. Chandrashekara, S.; Santosh, U. An Efficient Solution to the Free Vibration of Thick

Angle-ply Laminae. Int. J. Solids Struct. 1991, 27 (8) 999–1010.

16. Monsivais, G.; Rodriguez-Ramos, R.; Esquivel-Sirvent, R.;

Fernandez-Alvarez, L. Stark-Ladder Resonances in Piezoelectric Composites. Phys. Rev. B

2003, 68 174109.

17. Ambati, M.; Fang, N.; Sun, C.; Zhang, X. Surface Resonant States and Superlensing in

Acoustic Metamaterials. Phys. Rev. B 2007, 75 195447.

18. Wang, Z. G.; Lee, S. H.; Kim, C. K.; Park, C. M.; Nahm, K.; Nikitov, S. A. Effective

Medium Theory of the One-Dimensional Resonance Phononic Crystal. J. Phys.: Condens.

Mat. 2008, 20 055209.

19. Lange, K.; Rapp, B. E.; Rapp, M. Surface Acoustic Wave Biosensors: A Review. Anal.

Bioanal. Chem. 2008, 391 1509–1519.

20. Leisure, R. G.; Willis, F. A. Resonant Ultrasound Spectroscopy. J. Phys.: Condens. Mat.

1997, 9 6001–6029.

21. Yang, J.; Yang, Z. Analytical and Numerical Modeling of Resonant Piezoelectric Devices in

China-A Review. Sci. China Ser. G 2008, 51 (12) 1775–1807.

22. Abdulrehem, M. M.; Ott, E. Low Dimensional Description of Pedestrian-Induced

Oscillation of the Millenium Bridge. Chaos 2009, 19 013129.

23. Broer, W.; Hoenders, B. J. Natural Modes and Resonances in a Dispersive Stratified

N-Layer Medium. J. Phys. A: Math. Theor. 2009, 42 (24) 245207.

24. Fedorchenko, A. I.; Stachiv, I.; Wang, A. B.; Wang, W. H. Fundamental Frequencies of

Mechanical Systems with N-Piecewise Constant Properties. J. Sound Vibrat. 2008, 317

490–495.

25. Hsueh, W. Analytical Solution of Harmonic Travelling Waves in N-segment Strings. Proc.

R. Soc. Lond. A 2000, 456 2115–2126.

26. Gaudet, S.; Gauthier, C.; LeBlanc, V. G. On the Vibrations of an N-String. J. Sound Vibrat.

2000, 238(1) 147–169.

42
27. Churchill, R. R. Operational Mathematics; McGraw-Hill: New York, 1972.

28. Dhia, A. S. B. B.; Mercier, J. F. Resonances of an Elastic Plate in a Compressible Confined

Fluid. Q. J. Mech. Appl. Math. 2007, 60 (4) 397–421.

29. Caviglia, G.; Morro, A. A Closed-Form Solution for Reflection and Transmission of

Transient Waves in Multilayers. J. Acous. Soc. Am. 2004, 116 (2) 643–654.

30. Kaplunov, J.; Krynkin, A. Resonance Vibrations of an Elastic Interfacial Layer. J. Sound

Vibrat. 2006, 294 663–677.

31. Qiang, G.; Wanxie, Z.; Howson, W. P. A Precise Method for Solving Wave Propagation

Problems in Layered Anisotropic Media. Wave Motion 2004, 40 191–207.

32. Fokina, M. S.; Fokin, V. N. Resonances of Acoustic Waves Interacting with an Elastic

Seabed. J. Comput. Acoust. 2001, 9 (3) 1079–1093.

33. Graff, K. F. Wave Motion in Elastic Solids; Dover Publications: New York, 1975.

34. Goupillaud, P. An Approach to Inverse Filtering of Near-Surface Layer Effects from Seismic

Records. Geophysics 1961, 36 754–760.

35. Wang, N.; Lin, J. X.; Sadayuki, U. Goupillaud Inverse Problem with Arbitrary Input. J.

Acous. Soc. Am. 1997, 101 (6) 3255–3260.

36. Knopoff, L. A Matrix Method for Elastic Wave Problems. Bull. Seismol. Soc. Am. 1964, 54

431–438.

37. Thomson, W. T. Transmission of Elastic Waves Through a Stratified Solid Medium. J. Appl.

Phys. 1950, 21 89–93.

38. Haskell, N. A. The Dispersion of Surface Waves on Multilayered Media. Bull. Seismol. Soc.

Am. 1953, 43 17–34.

39. Bube, K.; Burridge, R. The One Dimensional Inverse Problem of Reflection Seismology.

Siam Rev. 1983, 25 (4) 497–559.

40. Velo, A. P.; Gazonas, G. A.; Ameya, T. z-Transform Methods for the Optimal Design of

One-Dimensional Layered Elastic Media. J. SIAM Appl. Math. 2009, 70 (3) 762–788.

41. Jury, E. I. Theory and Application of the z-Transform Method; John Wiley & Sons, Inc.:

1964.

43
42. Muir, T. A Treatise on the Theory of Determinants; Dover Publications: New York, 1960.

43. Proakis, J. G.; Manolakis, D. G. Introduction to Digital Signal Processing; Macmillan

Publishing Company: 1988.

44. Spiegel, M. R.; Liu, J. Mathematical Handbook of Formulas and Tables, Schaum’s Outline

Series; 2nd Edition ed.; McGraw-Hill: 1999.

45. Ashby, M. F. Materials Selection in Mechanical Design; 2nd Edition

ed.; Butterworth-Heinemann: Oxford, UK, 2000.

44
Appendix A. Stress Formulas for a Heaviside Loading

Some of the definitions and formulas in appendices A–E refer to the results from Velo and

Gazonas (40). These results were obtained when studying stress propagation in a

Goupillaud-type layered strip, subjected to zero initial conditions, a Heaviside step in stress

loading p = constant at time t = 0 at one end and held fixed at the other end. As before, the

design parameters are given by χm−1 =

mQ−1

i=1

(1 + αi) for 2 ≤ m ≤ 5, Γ3 = α1α3 − 1, and

Γ4 = α1α3α4 + α1α2α4 + α1α4 + α2α4 + α1α3 − 1, where χm−1 > 1 for m ≥ 2 and Γm−1 > −1

for m = 4, 5. Here we recover from the general formulation in equations 5 and 9, the stress

formulas derived in Velo and Gazonas (40) when the strip is subjected to a Heaviside step in

stress loading p at time t = 0 on one end and held fixed on the other end. This case is illustrated

in figure 1 for the choice of f(n) = p for n ≥ 0. Indeed, in this case we have

F(z) = Z(f(n)) = z

z−1

p and f(0) = p, which imply that (F(z) − f(0)) = p

z−1

. As a result,

equation 9 becomes

~xm =

η1z(F(z) − f(0))

|Am|





























|A1,1|

−|A1,2|

.

.

.

(−1)1+m|A1,m|





























=

η1pz

(z − 1)|Am|





























|A1,1|

−|A1,2|

.

.

.

(−1)1+m|A1,m|





























, (A-1)

simplifies to,

~xm(i) = hpz

η1(−1)i|A1,i|

(z−1)|Am|

i

= pz "

ai,0

z−1 +

bi,0

z+1 +

b m

2P

c

k=1

a

∗

i,kz+b

∗

i,k

z

2−2z cos θk+1 #

. (A-2)

After applying the inverse z-transform,

si(n) = Z

−1

(~xm(i)) = Z

−1

"

ai,0z

z−1 +

bi,0z

z+1 +

b m

2P

c

k=1

a

∗

i,kz

2+b∗

i,kz

z

2−2z cos θk+1 #

· p

= Z

−1

"

ai,0z

z−1 +

bi,0· z

z+1 +

b m

2P

c

k=1

3

ai,k· z(z−cos θk)

z

2−2z cos θk+1 +

bi,k·z sin θk

z

2−2z cos θk+1  ́

#

· p,

(A-3)

45
the stress formulas follow,

si(n) = "ai,0 + bi,0(−1)n +

b m

2P

c

k=1

ai,k cos (nθk) + bi,k sin (nθk)

#

· p, (A-4)

for the choices of the coefficients given below

ai,k = a

∗

i,k, bi,k =

a

∗

i,k cos θk + b

∗

i,k

sin θk

, for k = 1, . . . , b

m

2

c. (A-5)

As before, sin θk 6= 0 for all k = 1, . . . , b

m

2

c and bi,0 = 0 for m even. In addition, the

trigonometric sums in the stress equations A-4, allow us to express them in terms of the

Chebyshev polynomials of the first and second kind Tn(y) and Wn(y), by applying equation 68

and equation 69 as follows

si(n) = "ai,0 + bi,0Tn(−1) +

b m

2P

c

k=1

(ai,k − bi,k tan θk) · Tn(cos θk) + bi,k tan θk · Wn(cos θk)

#

· p, (A-6)

for i = 1, 2, . . . m and n ≥ 0. This representation is valid when cos θk 6= 0 for all 1 ≤ k ≤ bm

2

c.

The coefficients ai,j , bi,j , as well as the angles θk depend only on the impedance ratios αi,

i = 1, 2, . . . , m, j = 0, 1, . . . , b

m

2

c and k = 1, . . . , b

m

2

c.

46
Appendix B. Zeros of the Determinant |Am|

The expressions for the determinant |Am| using the recursive equation 10 for 1 ≤ m ≤ 5 are:

|A1| = z + 1, |A2| = (z + 1)2 − η2η1α1z = z

2 −

2(χ1−2)

χ1

z + 1,

|A3| = (z + 1) · [ z

2 −

2(χ2−2)

χ2

z + 1 ], |A4| = z

4 −

4Γ3

χ3

z

3 +

2(4Γ3−χ3+8)

χ3

z

2 −

4Γ3

χ3

z + 1,

|A5| = (z + 1) · [ z

4 −

4Γ4

χ4

z

3 +

2(4Γ4−χ4+8)

χ4

z

2 −

4Γ4

χ4

z + 1 ].

(B-1)

The determinant |Am| is a palindromic polynomial with real coefficients, as expected. All the

zeros of |Am| for 1 ≤ m ≤ 5 are on the unit circle and distinct based on the following criteria.

1. The zeros of a quadratic factor with real coefficients of the form z

2 + μz + 1 are distinct

(complex conjugate) on the unit circle iff μ

2 − 4 < 0 or equivalently iff |μ| < 2. When m = 2, 3

we have that |μ| < 2 ⇔

 ̄

 ̄

 ̄

χm−1−2

χm−1

 ̄

 ̄

 ̄

< 1, which is true because χm−1 > 1.

2. The zeros of a quartic factor with real coefficients of the form z

4 + μz3 + νz2 + μz + 1 are

distinct (complex conjugate) on the unit circle 2iff [ μ

2 − 4ν + 8 > 0 and |

−μ±

√

μ2−4ν+8

2

| < 2 ].

When m = 4, 5 we have the condition:

μ

2 − 4ν + 8 > 0 ⇔ [ (χm−1 − Γm−1)2 − 4χm−1 ] = [ (Υm−1 − 2)2 − 4(1 + Γm−1) ] > 0, where

Υm−1 = χm−1 − Γm−1.

The inequality follows from the following relations:

(Υ3 − 2)2 − 4(1 + Γ3) > (α1 + α3)

2 − 4α1α3 = (α1 − α3)2 ≥ 0 for m = 4, and

(Υ4 − 2)2 − 4(1 + Γ4) > (α1 + α2 + α1α2 + α3 + α4 + α3α4)

2 − 4(1 + Γ4) >

(−α1 − α2 − α1α2 + α3 + α4 + α3α4)

2 ≥ 0 for m = 5.

The condition

 ̄

 ̄

 ̄

 ̄

−μ±

√

μ2−4ν+8

2

 ̄

 ̄

 ̄

 ̄

< 2 ⇔

 ̄

 ̄

 ̄

 ̄

Γm−1±

√

(χm−1−Γm−1)

2−4χm−1

χm−1

 ̄

 ̄

 ̄

 ̄

< 1, which is true because

χm−1 > 1 > 0 and Γm−1 + 1 > 0 for m = 4, 5.

2This criterion is derived after dividing the quartic polynomial by z2

and then applying the transformation y =

z +

1

z

.

47
INTENTIONALLY LEFT BLANK.

48
Appendix C. Definition of the Matrix Dm

The following equation is the definition of matrix Dm:

Dm =



































F B1 0 . . . . . . 0

C2 F B2 0 . . . . . . 0

0 C3 F B3 0 . . . . . . 0

.

.

.

.

.

.

.

.

.

0 . . . . . . 0 Ci F Bi 0 . . . . . . 0

.

.

.

.

.

.

.

.

. 0

0 . . . . . . 0 Cm−1 F Bm−1

0 . . . . . . 0 Cm F



































2m×2m

(C-1)

where F = I + RT, Bi =

4α

2

i

(1+αi)

2 I, Ci =

4

(1+αi)

2 RT, Cm = 4RT, 1 ≤ i ≤ m − 1.

Here θ represents any of the angles θj, j = 0, . . . , b

m

2

c, as applicable; I represents the 2 × 2

identity matrix; and RT represents the transpose of the rotation matrix R =

"

cos θ − sin θ

sin θ cos θ

#

.

49
INTENTIONALLY LEFT BLANK.

50
Appendix D. Necessary Relations Among the Base Angles and Design

Parameters

The angles identified here for a strip with up to five layers relate to the zeros of |Am|, already

shown to be on the unit circle in appendix B. The 2b

m

2

c order equation for finding the non-trivial

zeros of |Am| in equation 13 is reduced to an b

m

2

c order equation for the cosine of the base

angles. The inequalities satisfied by the design parameters χm−1 for m = 2, 3 and χm−1, Γm−1

for m = 4, 5 in appendix B, validate the existence of these angles.

The two-layer case. The condition |D2| = 0 equivalent to hcos θ −

χ1−2

χ1

i2

= 0, implies that

−1 < cos θ =

χ1 − 2

χ1

< 1 =⇒ 0 < θ1 = cos−1

μ

χ1 − 2

χ1

¶

< π, (D-1)

The three-layer case. The condition |D3| = 0 is equivalent to

·

cos θ = −1 or 3cos θ −

χ2−2

χ2

 ́2

= 0  ̧, which based on the fact that

χ2 = (1 + α1)(1 + α2) > 1, leads to the relations

cos θ = −1 or − 1 < cos θ =

χ2 − 2

χ2

< 1, (D-2)

and the corresponding angle solutions

θ0 = cos−1(−1) = π, 0 < θ1 = cos−1

μ

χ2 − 2

χ2

¶

< π. (D-3)

The four-layer case. The condition |D4| = 0 ⇔ [χ3 cos2θ − 2Γ3 cos θ + (2Γ3 − χ3 + 4)]2 = 0,

leads to the quadratic equation:

cos2θ −

2Γ3

χ3

cos θ +

2Γ3 − χ3 + 4

χ3

= 0, (D-4)

and implies the following base angle solutions for any given set of the design parameters χ3 and

Γ3,

0 < θ1,2 = cos−1

Ã

Γ3 ±

p

(χ3 − Γ3)

2 − 4χ3

χ3

!

< π. (D-5)

51
The five-layer case.

|D5| = 0 ⇔ [ cos θ = −1 or [ χ4 cos2θ − 2Γ4 cos θ + (2Γ4 − χ4 + 4) ]2 = 0 ], leads to the

equations:

cos θ = −1 or cos2θ −

2Γ4

χ4

cos θ +

2Γ4 − χ4 + 4

χ4

= 0, (D-6)

and the corresponding angle solutions:

θ0 = cos−1(−1) = π, 0 < θ1,2 = cos−1

Ã

Γ4 ±

p

(χ4 − Γ4)

2 − 4χ4

χ4

!

< π, (D-7)

for any given set of the design parameters χ4 and Γ4.

52
Appendix E. Optimal Base Angles

For an m-layered Goupillaud-type elastic strip subjected to a Heaviside loading at one end and

held fixed at the other end, the optimal designs have the smallest stress amplitude of double the

loading. These optimal designs are characterized by the following optimal values of their base

angles for j = 2, 3, 4, . . .:

Table E-1. Optimal base angle values for the m-layer case.

m 2 3 4 5

Optimal Base Angles θ1,opt =

π

j

θ1,opt =

π

2j−1

θ1,opt =

π

2j

, θ1,opt =

π

2j+1 ,

θ2,opt = π −

π

2j

θ2,opt =

(2j−1)π

2j+1

53
INTENTIONALLY LEFT BLANK.

54
Appendix F. Solutions of the Frequency Equations for a Free-fixed

m-layered Goupillaud-type Elastic Strip (2 ≤ m ≤ 5)

F-1 Natural Frequency Spectrum for the Free-fixed Two-layer Case

After a few mathematical manipulations involving trigonometric identities, the frequency

equation 57 for m = 2 becomes

cos τω =

α1 − 1

α1 + 1

=

χ1 − 2

χ1

, where χ1 = (1 + α1) > 1 or α1 > 0. (F-1)

This suggests the following natural frequency spectrum

ωl1 =

1

τ

·

·

cos−1

μ

χ1 − 2

χ1

¶

+ 2l1π

 ̧

, ωl2 =

1

τ

·

·

− cos−1

μ

χ1 − 2

χ1

¶

+ 2(l2 + 1)π

 ̧

(F-2)

for l1, l2 = 0, 1, 2, . . .. Here α1 represents the impedance ratio between the first and the second

layer, while χ1 = 1 + α1. When both layers are occupied by the same material (α1 = 1),

equation F-1 recovers a known literature result for the natural frequency spectrum of a

homogeneous strip of length τ given in Graff (33). In this case the frequency equation F-1

becomes cos τω = 0, which implies the natural frequency spectrum ωl =

(2l+1)

τ

·

π

2

, l = 0, 1, 2 . . ..

Notice that equation F-1 implies that cos τω 6= ±1.

F-2 Natural Frequency Spectrum for the Free-fixed Three-layer Case

Using the fact that |∆m| = |∆ ̃ m|, the frequency equation 57 for m = 3 becomes equivalent to

solving the following equations:

cos

2τ

3

ω = −1 or cos

2τ

3

ω =

χ2 − 2

χ2

, where χ2 = (1 + α1)(1 + α2) > 1. (F-3)

From equation F-3, the frequency spectrum for the three-layer is







ωl0 = (2l0 + 1) 3π

2τ

,

ωl1 =

3

2τ

·

h

cos−1

3

χ2−2

χ2

 ́

+ 2l1π

i

, ωl2 =

3

2τ

·

h

− cos−1

3

χ2−2

χ2

 ́

+ 2(l2 + 1)π

i

,

(F-4)

55
where l0, l1, l2 = 0, 1, 2, . . ..

F-3 Natural Frequency Spectrum for the Free-fixed Four-layer Case

Use of the double angle formula for the cosine for the four-layer case (m = 4), results in

|∆4| =

1

4

·

£

χ3 cos2 τ

2

ω − 2Γ3 cos τ

2

ω + (2Γ3 − χ3 + 4)¤= 0, or equivalently

cos2

τ

2

ω −

2Γ3

χ3

cos

τ

2

ω +

2Γ3 − χ3 + 4

χ3

= 0. (F-5)

For any given values of the design parameters χ3 and Γ3, the frequency spectrum is:







ωl1 =

2

τ

·

·

cos−1

μ

Γ3±

√

(χ3−Γ3)

2−4χ3

χ3

¶

+ 2l1π

 ̧

,

l1, l2 = 0, 1, 2 . . . ,

ωl2 =

2

τ

·

·

− cos−1

μ

Γ3±

√

(χ3−Γ3)

2−4χ3

χ3

¶

+ 2(l2 + 1)π

 ̧

.

(F-6)

When the first two layers are occupied by one material (α1 = 1) and the other two are occupied

by another material (α2 6= 1, α3 = 1), the four-layer case reduces to the two-layer case and

equation F-5 reduces to equation F-1. From equation F-5 if follows that

cos2

τ

2

ω =

α2

α2 + 1

, (F-7)

which after applying the double angle formula for cosine becomes consistent with equation F-1,

cos τω =

α2 − 1

α2 + 1

. (F-8)

F-4 Natural Frequency Spectrum for the Free-fixed Five-layer Case

Similar to the previous cases, using symbolic algebra software and trigonometric identities we

obtain for m = 5 that |∆5| = cos τω

5

·

£

χ4 cos4 τ

5

ω − (χ4 + Γ4) cos2 τ

5

ω + (Γ4 + 1)¤= 0. After

using the double angle trigonometric formulas we have that

cos

2τ

5

ω = −1 or cos2

2τ

5

ω −

2Γ4

χ4

cos

2τ

5

ω +

2Γ4 − χ4 + 4

χ4

= 0. (F-9)

The frequency spectrum can then be obtained for any given values of χ4 and Γ4.

56
Appendix G. Material Properties

The elastic modulus E and density ρ of the materials used in subsection 4.1 (see Ashby (45)) are

as follows:

• Aluminum alloy: E = 70 GPa and ρ = 2, 500 kg/m3,

• Brass: E = 20 GPa and ρ = 984.19 kg/m3,

• Bronze: E = 50 GPa and ρ = 7, 189.35 kg/m3,

• Copper alloy: E = 10 GPa and ρ = 515.19 kg/m3,

• Glass fiber reinforced material: E = 30 GPa and ρ = 1, 130 kg/m3,

• Lead: E = 14 GPa and ρ = 11, 340 kg/m3,

• Nickel alloy: E = 200 GPa and ρ = 348.18 kg/m3,

• Tungsten alloy: E = 275 GPa and ρ = 19, 610 kg/m3,

• Steel: E = 200 GPa and ρ = 8, 000 kg/m3.

57
INTENTIONALLY LEFT BLANK.

58
NO. OF

COPIES ORGANIZATION

1

(PDF

ONLY)

DEFENSE TECHNICAL

INFORMATION CTR

DTIC OCA

8725 JOHN J KINGMAN RD

STE 0944

FORT BELVOIR VA 22060-6218

1 DIRECTOR

US ARMY RESEARCH LAB

IMNE ALC HRR

2800 POWDER MILL RD

ADELPHI MD 20783-1197

1 DIRECTOR

US ARMY RESEARCH LAB

RDRL CIM L

2800 POWDER MILL RD

ADELPHI MD 20783-1197

1 DIRECTOR

US ARMY RESEARCH LAB

RDRL CIM P

2800 POWDER MILL RD

ADELPHI MD 20783-1197

1 DIRECTOR

US ARMY RESEARCH LAB

RDRL D

2800 POWDER MILL RD

ADELPHI MD 20783-1197

ABERDEEN PROVING GROUND

1 DIR USARL

RDRL CIM G (BLDG 4600)

59
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

2 NSF

S MCKNIGHT

G PAULINO

4201 WILSON BLVD, STE 545

ARLINGTON, VA, 22230-0002

2 DARPA

L CHRISTODOULOU

W COBLENZ

3701 N FAIRFAX DR

ARLINGTON VA 22203-1714

1 DIRECTOR

US ARMY ARDEC

AMSRD AAR AEE W

E BAKER

BLDG 3022

PICATINNY ARSENAL NJ

07806-5000

2 US ARMY TARDEC

AMSTRA TR R MS 263

K BISHNOI

D TEMPLETON MS 263

WARREN MI 48397-5000

1 COMMANDER

US ARMY RSRCH OFC

RDRL ROE N

B LAMATTINA

PO BOX 12211

RESEARCH TRIANGLE PARK NC

27709-2211

1 COMMANDER

US ARMY RSRCH OFC

RDRL ROI M

J LAVERY

PO BOX 12211

RESEARCH TRIANGLE PARK NC

27709-2211

1 COMMANDER

US ARMY RSRCH OFC

RDRL ROE M

D STEPP

PO BOX 12211

RESEARCH TRIANGLE PARK NC

27709-2211

5 NAVAL RESEARCH LAB

E R FRANCHI CODE 7100

M H ORR CODE 7120

J A BUCARO CODE 7130

J S PERKINS CODE 7140

S A CHIN BING CODE 7180

4555 OVERLOOK AVE SW

WASHINGTON DC 20375

1 DTRA

M GILTRUD

8725 JOHN J KINGMAN RD

FORT BELVOIR VA 22060

1 ERDC

US ARMY CORPS OF ENGINEERS

USACEGSL

P PAPADOS

7701 TELEGRAPH RD

ALEXANDRIA VA 22315

1 AFOSR/NL

875 NORTH RANDOLPH ST

SUITE 325, RM 3112

F FAHROO

ARLINGTON VA 22203

1 CLEMSON UNIV

DEPT MECH ENGINEERS

M GRUJICIC

241 ENGRG INNOVATION BLDG

CLEMSON SC 29634-0921

1 UNIV OF CALIFORNIA

CTR OF EXCELLENCE FOR

ADV MATLS

S NEMAT NASSER

SAN DIEGO CA 92093-0416

5 DIRECTOR

LANL

P MAUDLIN

R GRAY

W R THISSELL

A ZUREK

F ADDESSIO

PO BOX 1663

LOS ALAMOS NM 87545

60
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

7 DIRECTOR

SANDIA NATL LABS

J BISHOP MS 0346

E S HERTEL JR MS 0382

W REINHART MS 1181

T VOGLER MS 1181

L CHHABILDAS MS 1811

M FURNISH MS 1168

M KIPP MS 0378

PO BOX 5800

ALBUQUERQUE NM 87185-0307

1 DIRECTOR

LLNL

M J MURPHY

PO BOX 808

LIVERMORE CA 94550

3 CALTECH

M ORTIZ MS 105 50

G RAVICHANDRAN

T J AHRENS MS 252 21

1201 E CALIFORNIA BLVD

PASADENA CA 91125

5 SOUTHWEST RSRCH INST

C ANDERSON

K DANNEMANN

T HOLMQUIST

G JOHNSON

J WALKER

PO DRAWER 28510

SAN ANTONIO TX 78284

1 TEXAS A&M UNIV

DEPT OF MATHEMATICS

J WALTON

COLLEGE STATION TX 77843

1 UNIVERSITY OF MISSISSIPPI

DEPT OF MECH ENGRG

A M RAJENDRAN

201-B CARRIER HALL

UNIVERSITY, MS 38677

2 SRI INTERNATIONAL

D CURRAN

D SHOCKEY

333 RAVENSWOOD AVE

MENLO PARK CA 94025

1 VIRGINIA POLYTECHNIC INST

COLLEGE OF ENGRG

R BATRA

BLACKSBURG VA 24061-0219

7 UNIV OF NEBRASKA

DEPT OF ENGRG MECH

F BOBARU

Y DZENIS

G GOGOS

M NEGAHBAN

R FENG

J TURNER

Z ZHANG

LINCOLN NE 68588

1 JOHNS HOPKINS UNIV

DEPT OF MECH ENGRG

K T RAMESH

LATROBE 122

BALTIMORE MD 21218

1 WORCESTER POLYTECHNIC INST

MATHEMATICAL SCI

K LURIE

WORCESTER MA 01609

4 UNIV OF UTAH

DEPT OF MATH

A CHERKAEV

E CHERKAEV

E S FOLIAS

R BRANNON

SALT LAKE CITY UT 84112

1 PENN STATE UNIV

DEPT OF ENGRG SCI & MECH

F COSTANZO

UNIVERSITY PARK PA 168023

4 UNIV OF DELAWARE

DEPT OF MECH ENGRG

T BUCHANAN

T W CHOU

A KARLSSON

M SANTARE

126 SPENCER LAB

NEWARK DE 19716

1 UNIV OF DELAWARE

CTR FOR COMPST MATRLS

J GILLESPIE

NEWARK DE 19716

61
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

1 COMPUTATIONAL MECH

CONSULTANTS

J A ZUKAS

PO BOX 11314

BALTIMORE MD 21239-0314

1 LOUISIANA STATE UNIV

R LIPTON

304 LOCKETT HALL

BATON ROUGE LA 70803-4918

1 INST OF ADVANCED TECH

UNIV OF TX AUSTIN

S BLESS

3925 W BRAKER LN STE 400

AUSTIN TX 78759-5316

1 APPLIED RSCH ASSOCIATES

D E GRADY

4300 SAN MATEO BLVD NE

STE A220

ALBUQUERQUE NM 87110

1 INTERNATIONAL RSRCH

ASSOC INC

D L ORPHAL

4450 BLACK AVE

PLEASANTON CA 94566

3 ORNL

ENVIRONMENTAL SCI DIV

W DOLL

T GAMEY

L BEARD

PO BOX 2008

OAK RIDGE TN 37831

2 WASHINGTON ST UNIV

INST OF SHOCK PHYSICS

Y M GUPTA

J ASAY

PULLMAN WA 99164-2814

1 NORTHWESTERN UNIV

DEPT OF CIVIL & ENVIRON ENGRG

Z BAZANT

2145 SHERIDAN RD A135

EVANSTON IL 60208-3109

1 UNIV OF DAYTON

RSRCH INST

N S BRAR

300 COLLEGE PARK

MS SPC 1911

DAYTON OH 45469

2 TEXAS A&M UNIV

DEPT OF GEOPHYSICS MS 3115

F CHESTER

T GANGI

COLLEGE STATION TX 778431

1 UNIV OF SAN DIEGO

DEPT OF MATH & CMPTR SCI

A VELO

5998 ALCALA PARK

SAN DIEGO CA 92110

1 NATIONAL INST OF

STANDARDS & TECHLGY

BLDG & FIRE RSRCH LAB

J MAIN

100 BUREAU DR MS 8611

GAITHERSBURG MD 20899-8611

1 MIT

DEPT ARNTCS ASTRNTCS

R RADOVITZKY

77 MASSACHUSETTS AVE

CAMBRIDGE MA 02139

1 UNIV OF DELAWARE

DEPT ELECTRICAL & CMPTR ENGRG

D WEILE

NEWARK DE 19716

1 T W WRIGHT

4906 WILMSLOW RD

BALTIMORE MD 21210

1 UNIV OF TEXAS-PAN AMERICAN

COLLEGE OF ENGRG

& COMPUTER SCI

D H ALLEN

1201 WEST UNIVERSITY DR

EDINBURG, TX 78539-2999

3 RDRL D

J CHANG

R SKAGGS

V WEISS

BLDG 205

2800 POWDER MILL RD

ADELPHI MD 20783-1197

ABERDEEN PROVING GROUND

76 DIR USARL

62
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

RDRL WM

B FORCH

S KARNA

J MCCAULEY

P PLOSTINS

M ZOLTOSKI

RDRL WML

D LYONS

J NEWILL

RDRL WML B

I BATYREV

B RICE

N WEINGARTEN

RDRL WML D

P CONROY

M NUSCA

RDRL WML G

M BERMAN

W DRYSDALE

RDRL WML H

D SCHEFFLER

S SCHRAML

B SCHUSTER

RDRL WMM

J BEATTY

R DOWDING

J ZABINSKI

RDRL WMM A

M MAHER

J TZENG

E WETZEL

RDRL WMM B

T BOGETTI

B CHEESEMAN

C FOUNTZOULAS

A GAZONAS

G GAZONAS

D HOPKINS

P MOY

B POWERS

C RANDOW

T SANO

M VANLANDINGHAM

R WILDMAN

C F YEN

RDRL WMM C

J LA SCALA

RDRL WMM D

E CHIN

K CHO

RDRL WMM E

M COLE

T JESSEN

J LASALVIA

M PEPI

J SANDS

RDRL WMM F

L KECSKES

H MAUPIN

RDRL WML G

J ANDZELM

A RAWLETT

RDRL WMP

P BAKER

S SCHOENFELD

RDRL WMP B

R BECKER

S BILYK

D CASEM

J CLAYTON

M GREENFIELD

C HOPPEL

R KRAFT

B LEAVY

B LOVE

M RAFTENBERG

M SCHEIDLER

T WEERASOOYIA

RDRL WMP C

T BJERKE

S SEGLETES

RDRL WMP D

R DONEY

J RUNYEON

B SCOTT

RDRL WMP E

M BURKINS

RDRL WMP F

M CHOWDHURY

A FRYDMAN

N GNIAZDOWSKI

R GUPTA

RDRL WMP G

N ELDREDGE

D KOOKER

S KUKUCK

G R PEHRSON

63
INTENTIONALLY LEFT BLANK.

64