Ordered and Highly Scalable Granular Media 

for Shock Mitigation 

by Robert Doney and Surajit Sen 

ARL-TR-3612 September 2005 

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-3612 September 2005 

Ordered and Highly Scalable Granular Media 

for Shock Mitigation 

Robert Doney 

Weapons and Materials Research Directorate, ARL 

Surajit Sen

State University of New York, Buffalo, NY 

Approved for public release; distribution is unlimited. 
ii 

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) 

September 2005 

2. REPORT TYPE 

Final 

3. DATES COVERED (From - To) 

June 2004–April 2005 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

Ordered and Highly Scalable Granular Media for Shock Mitigation 

5c. PROGRAM ELEMENT NUMBER 

5d. PROJECT NUMBER 

AH80 

5e. TASK NUMBER 

6. AUTHOR(S) 

Robert Doney and Surajit Sen*

5f. WORK UNIT NUMBER 

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

U.S. Army Research Laboratory 

ATTN: AMSRD-ARL-WM-TA 

Aberdeen Proving Ground, MD 21005-5066 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-TR-3612 

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 

Approved for public release; distribution is unlimited. 

13. SUPPLEMENTARY NOTES 

*

State University of New York, Buffalo, NY 

14. ABSTRACT 

We provide an update to our current efforts in studying shock mitigation using one-dimensional arrays of adjacent metal spheres that 

may progressively shrink in size. These tapered chains are characterized by the number of spheres, N, the amount of tapering, q, and 

restitutive losses, ω. Spheres are assumed to interact via the nonlinear Hertz potential, V α δ5/2, where δ is the amount of overlap of 

adjacent grains and the constant of proportionality varies with the material. To gauge the ability of such chains to absorb transient 

pulses, we look at the normalized kinetic energy, KEN, as a function of N, q, ω. A hard-sphere approximation and numerical analysis 

is therefore performed for 3 ≤ N ≤ 20 and 0 ≤ q, ω ≤ 0.1 for two tapered chain architectures. Intended experimental studies are 

outlined as those that have already been performed and focus on solitary wave propagation in monodisperse (q = 0) chains. Results 

are quite encouraging and independent of system size; consequently, we propose a tapered chain armor panel design consisting of

many tapered chains working in tandem. The specific and gravimmetric absorbed energies are calculated for an example prototype. 

Intended work and further discussion are provided in the conclusion. 

15. SUBJECT TERMS

shock, mitigation, granular media, ballistic shock, Hertz, tapered chain, solitary wave 

16. SECURITY CLASSIFICATION OF: 

19a. NAME OF RESPONSIBLE PERSON

Robert Doney 

a. REPORT 

UNCLASSIFIED 

b. ABSTRACT 

UNCLASSIFIED

c. THIS PAGE 

UNCLASSIFIED

17. LIMITATION 

OF ABSTRACT 

UL 

18. NUMBER 

OF PAGES 

68 

19b. TELEPHONE NUMBER (Include area code) 

410-278-7309 

 Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18
iii 

Contents 

List of Figures v

List of Tables vi

1. Introduction 1

2. Granular Media and TCs 1

3. Hard-Sphere Approximation 3

3.1 The Simple Tapered Chain..............................................................................................3 

3.1.1 Lossless STC Hard-Sphere Approximation ..........................................................3 

3.1.2 Lossy STC Hard-Sphere Approximation ..............................................................4 

3.1.3 KE Parameter Space for STC Hard-Spheres.........................................................5 

3.2 The Decorated Tapered Chain.........................................................................................5 

4. Numerical Solution to the Equations of Motion 11

4.1 STC................................................................................................................................14 

4.1.1 STC KE and Force Parameter Spaces: KE N,q, ( ω) , F N,q, ( ω) ......................16 

4.1.2 STC Mathematical Model ...................................................................................19 

4.2 DTC...............................................................................................................................20 

5. Experimental Efforts 22

5.1 STC................................................................................................................................22 

5.2 DTC...............................................................................................................................23 

6. Precompressed Chains 27

7. TC Armor Panels 28

8. Concluding Remarks 31

9. References 32
iv 

Appendix A. Simple Tapered Chain Code 35

Appendix B. Decorated Tapered Chain Modifications 45

Appendix C. Practical Extration and Report Language (PERL) Script for Parametric 

Studies 49

Appendix D. MATLAB Code for Generating Numerical Kinetic Energy Surfaces 53

Appendix E. Normalization 57

Distribution List 58 
v 

List of Figures 

Figure 1. The simple tapered chain: N = 10, q = 8%, L = 70.7 mm, and ri = 5 mm. .....................2 

Figure 2. KE (N,q,EL) parameter space for the STC hard-sphere approximation. N varies 

from 3 to 25 and q from 0% to 20%. .........................................................................................6 

Figure 3. KE(N,q) parameter space for the STC hard-sphere approximation where the initial 

velocity is supplied to the smaller end of the chain...................................................................7 

Figure 4. The decorated tapered chain: N = 13, q = 8%, f = 0.7, L = 80.7 mm, and ri = 5 

mm. ............................................................................................................................................7 

Figure 5. KE as a function of N,q and f for DTC hard-spheres. Note that the colorscale is 

calibrated for each subplot.......................................................................................................12 

Figure 6. KE as a function of time for the smallest particle in a chain with 

N = 20, ω = 0.05, and q = [ ] 0 : 0.02 : 0.1 . Each plot is evaluated for a different tapering 

with initial KE as 0.0838 J.......................................................................................................15 

Figure 7. Numerical solution of KE N( ,q) parameter space for constant ω. Note that it 

would be very difficult and costly to produce a similar plot empirically. Each node on 

the surface represents a different TC and experiment. ............................................................17 

Figure 8. Numerical solution of F N( ) ,q parameter space for constant ω...................................18 

Figure 9. Difference plot of the lossless KE parameter spaces for the hard-sphere 

approximation (figure 2a) and numerical solution (7a). Note that the azimuthal view has 

been rotated 180° for clarity. ...................................................................................................19 

Figure 10. The simulated and modeled STC. ...............................................................................21 

Figure 11. DTC numerical solution of KE N( ,q, f) and sample DTCs for various f (inset). 

The scale is the same as that in figure 7...................................................................................22 

Figure 12. STC experimental setup by Job et al. (12). .................................................................23 

Figure 13. STC experimental setup. The left half of the figure represents a monodisperse 

chain while the right a TC. Lower plots denote later time......................................................24 

Figure 14. Preliminary DTC housing chamber.............................................................................25 

Figure 15. DTC constituents.........................................................................................................26 

Figure 16. DTC interstitial grains. ................................................................................................27 

Figure 17. KE as a function of tapering and restitution for constant N = 20 and ∆ = 0.03% . ....28 

Figure 18. STC armor plate panel.................................................................................................29 

Figure 19. SAE and VAE as a function of M, the total STC armor panel mass, and uf , the 

incident flyer plate velocity. ....................................................................................................30 
vi 

List of Tables 

Table 1. Material properties (28). .................................................................................................13 
1 

1. Introduction 

Our recent work (1, 2) in the area of granular physics indicates that certain (ordered with respect 

to size) arrangements of smooth, metal spheres adjacent to each other in one-dimensional (1-D) 

chains demonstrate noticeable shock absorption. This has been observed and measured 

analytically (1–5) and through numerical integration of the equations of motion (1, 6–8). 

Experimental efforts (9–12) have focused primarily on solitary wave propagation and 

demonstrating reduced impact force using tapered chains (TCs) (4). The system is scalable, as 

there does not appear to be any restrictions on particle size. There is, of course, limits imposed 

by quantum physics and manufacture of the spheres. Up to this point in our efforts, only single 

chains of varying length have been considered. However, increasingly encouraging results have 

prompted us to consider the collective disruptive effects of many chains as a defeat mechanism 

against ballistic shock and other undesirable transient pulses. 

The goal of this report is to summarize both published and unpublished results through 

normalized energy calculations and include an analytical description of a potential armor panel. 

The contents are organized as follows. A small overview of granular media and the ordered 

systems of interest are first discussed. This is followed by a hard-sphere approximation and the 

numerical solution to the equations of motion. In addition, our mathematical model provides 

some predictive capability for a certain range of systems. Past experimental arrangements and 

intended investigations are reviewed. In conclusion, we introduce a potential armor panel 

configuration with an estimate of the specific and volumetric absorbed energy and summarize 

the salient features of the report. 

2. Granular Media and TCs 

Granular media (13–17)—not to be confused with microscopic grains in a metal—are a curious 

thing. In general, its members consist of discrete particles that can range in size from 

micrometers to meters and number in quantity from several to the uncountable. The most wellknown constituents of this group are sands and powders and are utilized across many disciplines. 

In this communication, components of granular media may be referred to as grains, particles, 

beads, or spheres. These media are particularly intriguing in that among other things it has 

properties of liquids and solids and calls have been made to add it as a fundamental state of 

matter (13). Within the U.S. Army Research Laboratory, granular media have been investigated 

for understanding wave propagation and its mathematical similarity to nonlinear rods (18). 

While granular media have been fundamental in saving lives for decades against ballistics 

(sandbags), new features are being observed in the growing literature which suggest a more 
2 

technological role may be at hand. Particularly appealing and relatively simple to address is the 

behavior of 1-D systems where smooth elastic spheres intersect at a point in their initial state. 

Even though we may constrain the basis of variables in this system to just a few parameters, the 

complexity of its dynamics is formidable. We call these systems tapered chains (TCs). 

We define TCs (figure 1) as 1-D granular arrays of elastic spheres that touch at a single point in 

their initial state and grow to a disk under compression in the plane perpendicular to the figure. 

The chains can be characterized by the number of grains, N, the successive decrease in size of 

the grains or tapering, q, and restitutive losses, ω. Restitution (19–20) represents the expansion 

phase in a collision where energy is converted from elastic potential to kinetic and energy losses 

can be incorporated. 

Figure 1. The simple tapered chain: N = 10, q = 8%, L = 70.7 mm, and ri = 5 mm. 

It is intriguing and suggestive that, conceptually, granular media can be envisioned as the inverse 

of a porous material. In a porous material, there are gaseous (air) voids within a solid matrix. 

For dry granular media, solid voids (grains) exist within a gaseous matrix (air). Packing, 

however, limits the analogy; for grains, each entity is in contact with at least one other, whereas 

voids are not necessarily in contact with one another in porous materials. This is an interesting 

comparison to draw because the considerable amount of work required to collapse the voids in a 

porous material makes it a favorable technology against ballistic shock (21). Are there then 

opportunities in granular media for the armor designer to exploit? 

The interaction potential between adjacent elastic spheres was first identified by Hertz (22) circa 

1881. Modern and eloquent derivations were later performed by both Landau and Lifshitz (23) 

and Love (24), with a simplified order of magnitude approach by Leroy (25). More recently, 

Nesterenko (26) has written extensively on the dynamics of granular chains. If one can show 

significant reduction in output energy while ensuring that the structure of granular columns can 

be maintained, it is quite possible to use these TCs against propagating shocks. For now, we 

assume that one remains within the elastic regime. Future work will begin to relax this 

restriction.
3 

3. Hard-Sphere Approximation 

These approximations differ from the numerically simulated systems in two major ways. The 

first is that the chain is not bounded by fixed rigid walls. As a result, energy will not continue to 

be transmitted up and down the chain. The second is that the potential becomes infinite and as a 

consequence, the energy packet is only 1 grain in width. The system therefore propagates energy 

as independent collisions. This is congruent to the independent collision model proposed by Wu 

(4). In the numerical results, width of the energy pulse is a function of the tapering—and when 

all grains are the same size it is about 5 grains wide. By generating an iterative form of the 

conservation equations, one can arrive at an expression for the normalized kinetic energy (KE), 

KEN = KEout KEin . This ratio will be the primary variable determining the absorptive quality of 

TCs. 

3.1 The Simple Tapered Chain 

The simple tapered chain (STC) is displayed in figure 1. To generate an initial disturbance, an 

input velocity, vi

 , is applied to the rightmost and largest grain with radius, ri. It propagates to 

the left, encountering an initially stationary grain of radius,ri+1. The radius of the i +1 particle 

may be reduced by q% from ri

. This tapering q will be constant along the entire length of the 

chain. During the transmission of the impulse along the chain, there may be energy losses and 

we consider two cases described in the following subsections. 

3.1.1 Lossless STC Hard-Sphere Approximation 

Ignoring any energy loss during a collision, we perform an STC hard-sphere approximation for a 

TC. Masses and radii are expressed as 

ri+1 = ri − ri

q = (1− q)ri = εri , 

mi = ρVi = 4

3

π ri

3

ρ = ηri

3

, (1) 

and 

mi+1 = ηri+1

3 = ηε

3 ri

3

, (2) 

where ε =1− q. Evaluating the conservation of momentum with a single prime denoting post 

collision values and the initial condition that the i +1 particle is stationary before a collision 

( vi+1 = 0), all η cancel and we obtain 

mivi + mi+1vi+1 = mivi

′ + mi+1vi+1′

ri

3

vi = ri

3

vi

′ + ε

3

ri

3

vi+1

′

vi = vi

′ + ε

3

vi+1

′

, (3)
4 

where equations 1 and 2 have been used. Following the same procedure for the conservation of 

energy while ignoring the factor of one-half yields 

vi

2 = v ′

i

2 + ε

3

v ′

i+1

2 . (4) 

Letting A = ε

3

v ′ i+1, we can rewrite equation 3 in terms of v ′ i and substitute the resulting 

expression into equation 4, 

vi

2 = v( i − A)

2

+ Av ′ i+1

= vi

2 − 2Avi + A2 + Av ′ i+1

2vi = A + v ′ i+1

, 

v ′

i+1

vi

= 2

1+ ε

3 . (5) 

Note that for one collision, 

2 2

out i i i i 1 11 1 3

in i i i i

KE KE m v v

KE KE m v v

ε + ++ +

′ ⎛⎞ ⎛⎞ ′ ′ == = ⎜⎟ ⎜⎟ ⎝⎠ ⎝⎠

 (6) 

( )

3

1

2 3

4

1

i

i

KE

KE

ε

ε

+

′ = + . 

For N particles there will be N-1 collisions, each of which has the ratio in equation 6. Therefore, 

the normalized KE, KEN , for the lossless STC hard-sphere approximation is given as 

( )

( )

–1

3

2 3

4 1– =

1+ 1–

N

N

q KE

q

⎧ ⎫ ⎪ ⎪ ⎨ ⎬ ⎪⎡ ⎤ ⎪ ⎩ ⎭ ⎣ ⎦

. (7) 

3.1.2 Lossy STC Hard-Sphere Approximation 

The same approximation can be performed with some amount of energy loss, E  ̃ 

L , included such 

that the system is still conservative and may represent a better approximation. Consequently, the 

momentum equation 3 is unchanged, but equation 4 becomes 

vi

2 = v ′

i

2 + ε

3

v ′

i+1

2 + E  ̃ 

L . (8) 

One then obtains a more complicated expression replacing equation 5: 

v ′ i+1

vi

=

2 − E  ̃ 

L

ε

3

vi

v ′ i+1

1+ ε

3 . (9) 

We can make the substitution, E  ̃ 

L ∝vi

v ′

i+1or E  ̃ L = ELvi

v ′

i+1, where EL is the constant of 

proportionality. This adjustment yields
5 

( )

3

1

3 3

2 –

1

i L

i

v E

v

ε

ε ε

+

′ = + . 

The corresponding result for the normalized KE for N particles is 

( )

() ()

1 2 3

2 3 3

2 1

1 11

N

L

N

q E

KE

q q

−

⎧ ⎫ ⎡ ⎤ − − ⎪ ⎣ ⎦ ⎪ = ⎨ ⎬ ⎪ − +− ⎡ ⎤ ⎪ ⎩ ⎭ ⎣ ⎦

. (10) 

In the limit EL = 0, equation 10 reduces to the lossless case, equation 7, as one would expect. 

Note that results are independent of initial velocity and size of the grains. 

3.1.3 KE Parameter Space for STC Hard-Spheres 

Figure 2 highlights the behavior of equation 10 for 0 ≤ q ≤ 0.1, 3 ≤ N ≤ 20 and selected EL . The 

tapering q resembles a sigmoid or half-gaussian and is stretched to infinity for small N. In the 

limits that q = 0 or N = 1 (not shown), KE is unity. For the lossless, monodisperse chain (panel 

(a), q = 0), that energy is completely transferred regardless of the number of spheres. 

Interestingly, if the initial velocity is supplied to the smaller end of the TC, one also observes 

shock absorption similar to the system in figure 1, albeit with less efficiency. This is 

accomplished by adjusting the definition of tapering. In this case, subsequent particles are 

growing in size, i.e., ri+1 = (1 + q)ri. Equations 7 and 10 are modified accordingly, and the result 

for a lossless system is illustrated in figure 3. Apparently, both configurations mitigate a 

propagating pulse. 

3.2 The Decorated Tapered Chain 

The second and last TC geometry that we consider is the decorated tapered chain (DTC) 

(figure 4). This can be assembled from the STC by introducing a single-sized interstitial grain of 

radius reven between every member of an STC. We constrain the system to an odd number of 

particles such that the inserted grains are not at the boundary. Additionally, we presume that 

these interstitial grains will be equal to or smaller than the smallest member (rN ) of the tapered 

part of the chain; therefore, reven = frN , where 0 < f ≤1.0—although the flexibility is already 

built in for f to be any size. 
6 

Figure 2. KE (N,q,EL) parameter space for the STC hard-sphere approximation. N varies from 3 to 25 and q from 

0% to 20%. 
7 

Figure 3. KE(N,q) parameter space for the STC hard-sphere approximation where the initial 

velocity is supplied to the smaller end of the chain. 

Figure 4. The decorated tapered chain: N = 13, q = 8%, f = 0.7, L = 80.7 mm, and ri = 5 mm. 

It is immediately clear that size mismatch between neighboring grains is a function of position 

along the DTC. This is in stark contrast to the STC, where successive grains are always smaller 

(or larger) by the same amount. It is possible then to have DTC chains that appear to resemble 

monodisperse chains for only a portion of the chain. 
8 

Again, our primary interest is in deriving an expression for the normalized KE, 

2

1 11

out N N N

in

KE KE m v

KE KE m v

⎛ ⎞ ′ = = ⎜ ⎟ ⎝ ⎠

2

2 1 2

11 1 1

NN i i

N ii

mv vv v

mv v v v

+ +

− +

⎡⎛ ⎞ ⎛ ⎞⎛ ⎞ ′ ′′ ⎛ ⎞′ ⎤ = ⎢⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎥ ′ ′′ ⎣⎝ ⎠ ⎝ ⎠⎝ ⎠ ⎝ ⎠⎦

" " , (11) 

which can again be obtained by the conservation of momentum and energy where collisions are 

treated as independent. Since these are not trivial iterative solutions as in the STC hard-sphere 

approximation, we must solve for successive collisions and determine the overall pattern. We 

will eventually look for forms of v ′ i+1 vi and then generalize for N particles or N −1 collisions. 

First, the relationship among masses and radii must be evaluated. Given N(odd) particles in a 

DTC, every other grain will reduce in size by q% such that ri+2 = ε ri

, where ε = (1− q). 

Assembling the radii, we have 

ir (12) 

i N 1 r fr + = (13) 

2 (1 ) i ii i i r r qr q r r ε + = − =− = (14) 

i N 3 r fr + = (15) 

2

42 2 2 (1 ) ii i i i r r qr q r r ε ++ + + = − =− = (16) 

i N 5 r fr + = (17) 

3

64 4 4 (1 ) ii i i i r r qr q r r ε ++ + + = − =− = (18) 

#

N N 1 r fr − = (19) 

( ) 1 2

22 2 (1 ) NNN N N i r r qr q r r ε − = − =− = −− − (20) 

Equations 12–14 imply N = 3, equations 12–16 imply N = 5, etc. The major equations for radii 

are therefore 

( )

( )( ) ( )

( )

1 2

–1 2

1, 3, , 1

N

N i

N

ii N i

r r

r fr

ε

ε

−

++ −

=

= "

. (21) 
9 

Recall for masses that mi = ρVi = 4

3

πri

3

ρ = ηri

3

. Note that since η is just a constant and will 

cancel once the conservation equations are put into use, we will ignore it from now on. This 

expression for mi combined with equations for ri

 provide the relations 

( )

0 3

3 3

1 1

3 3

2

3 3

3 3

6 3

4

3 3

–1 1

3 12 3

i i

ii i

i i

ii i

i i

NN i

N

N i

m r

m r Ar

m r

m r Ar

m r

m r Ar

m r

ε

ε

ε

ε

+ +

+

+ +

+

−

−

∝

∝ =

∝

∝ =

∝

∝ =

∝

#

 (22) 

where A = f 3

ε

3( ) N−1 2

. We may now use these to evaluate the conservation equations. Beginning 

with momentum and assuming that each subsequent particle in the chain begins at rest, we solve 

for the first five collisions. Primes and double-primes indicate post-collision states. A primed 

quantity denotes the first post-collision state of a sphere which serves as input to the next 

collision. To keep track of its velocity after the second collision, it is denoted by a double-prime 

and will eventually be eliminated. With ε = (1− q), we obtain 

11 1 m v m v m v v v Av ii ii i i i i i + + + = + → =+ ′ ′ ′′ (23) 

3 m v m v m v Av Av v ii ii i i i i i 11 11 2 2 1 1 2 ε ++ ++ + + + + + ′ ′′ ′ ′ ′′ ′ = + → =+ (24) 

3 3 m v m v m v v v Av ii ii ii i i i 22 22 33 2 2 3 ε ε ++ ++ ++ + + + ′ ′′ ′ ′ ′′ ′ = + → =+ (25) 

6 m v m v m v Av Av v ii ii ii i i i 33 33 44 3 3 4 ε ++ ++ ++ + + + ′ ′′ ′ ′ ′′ ′ = + → =+ (26) 

6 6 m v m v m v v v Av ii ii ii i i i 44 44 55 4 4 5 ε ε ++ ++ ++ + + + ′ ′′ ′ ′ ′′ ′ = + → =+ (27) 

#

From the pattern in the previous formulas, equation 23 can be rewritten as ε

0

v ′ i = ε

0

v ′ i + Av ′

i+1. 

Evaluating the energy conservation equations yields the same form as equations 23–27 except 

velocities are squared: 

22 2

ii i 1 v v Av + = + ′ ′ , (28) 

2 2 32 Av Av v iii 11 2 ε + + + ′ = + ′′ ′ , (29) 

32 3 2 2

i ii 2 23 ε ε v v Av + + + ′ = + ′′ ′ , (30) 
10 

2 2 62 Av Av v iii 33 4 ε + + + ′ = + ′′ ′ , (31) 

62 6 2 2

i ii 4 45 ε ε v v Av + + + ′ = + ′′ ′ . (32) 

#

We can combine equations 23–27 and 28–32 to eliminate the double-primed terms and form the 

velocity ratios: v ′ i+1

vi

 , v ′ i+2

v ′ i+1

, etc. Beginning with equation 23, we isolate v ′ i and square it to 

obtain v ′ i

2 = vi2 − 2Avi

v ′ i+1 + A2v ′ i+1

2 . Next, substitute this into equation 28 and rearrange to obtain 

v ′ i+1

vi

. This is then repeated for equation pairs 24 and 29, 25 and 30, etc., to obtain the following 

ratios: 

0

1

0

i 2

i

v

v A

ε

ε

+

′ = +

, (33) 

2

3

1

i 2

i

v A

v A ε

+

+

′ = ′ +

, (34) 

3

3

3

2

i 2

i

v

v A

ε

ε

+

+

′ = ′ +

, (35) 

4

6

3

i 2

i

v A

v A ε

+

+

′ = ′ +

, (36) 

#

where we insert a term of ε

0

 in equation 33. With our goal being relation 11, we combine 

equations 33–36: 

12 1

1 1 21

Ni i N N

ii N N

v vv v v

v vv v v

++ −

+ −−

′ ′′ ′ ′ ⎛ ⎞⎛ ⎞ ⎛ ⎞⎛ ⎞ = ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ′ ′′ ⎝ ⎠⎝ ⎠ ⎝ ⎠⎝ ⎠

" (37) 

036

03 36 69

3

5

7

2 22 22 2

N

N

N

AAA

A AA AA A

εεε

εε εε εε

=

=

=

⎛⎞ ⎛⎞ ⎛⎞ ⎛⎞ ⎛⎞ ⎛⎞ = ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎝⎠ ⎝⎠

+ ++ ++ + 	 ⎝⎠ ⎝⎠ ⎝⎠ 


	


	


 (38) 

#
11 

After some observation, the ratio can be put into closed form to obtain 

( )

( )

( ) 1 2 3 1

3 1 3

1 1

2 2 N j

N

j j

j

v A

v A A

ε

ε ε

− −

−

=

′ ⎛ ⎞⎛ ⎞ = ⎜ ⎟⎜ ⎟ ∏ ⎝ ⎠ + ⎝ ⎠ +

 (39) 

( ) ( )

( )

( ) 1 2 3 1

1 1 2

3 1 3

1

1 2

N j

N N

j j

j

A

A A

ε

ε ε

− −

− −

−

=

⎛ ⎞⎛ ⎞ = ⎜ ⎟⎜ ⎟ ∏ ⎝ ⎠ + ⎝ ⎠ +

. (40) 

Turning to the mass ratios, it appears that most terms cancel. 

1 2

1 11

Ni i N

ii N

m mm m

m mm m

+ +

+ −

⎛ ⎞⎛ ⎞ ⎛ ⎞ = ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎝ ⎠⎝ ⎠ ⎝ ⎠

" (41) 

3 6

0 3

A A

A A

ε ε

ε ε

⎛⎞ ⎛⎞ ⎛⎞ ⎛⎞ = ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎝⎠ ⎝⎠ ⎝⎠", (42) 

leading to the simple expression 

3 12 ( )

1

mN N

m

ε − = . (43) 

We can now identify the normalized KE, equation 11, by combining equations 40 and 43 to form 

( )( ) ( )

( )

( )

2

1 2 3 1

1 3 2

3 1 3

1

4

N j N

N j j

j

KE A

A A

ε

ε

ε ε

− − −

−

=

⎧ ⎫ ⎪ ⎪ = ⎨ ⎬ ⎡ ⎤ + + ⎡ ⎤ ⎩ ⎭ ⎪ ⎣ ⎦⎣ ⎦ ⎪ ∏ , (44) 

with A = f 3

ε

3( ) N−1 2

 and ε = (1− q). This is plotted in figure 5, where it is clear that as f 

decreases—so that inertial mismatches increase—the normalized KE decays rapidly. Note that 

the colorscale is calibrated for each subplot. 

4. Numerical Solution to the Equations of Motion 

A more accurate method of evaluating the dynamics of TCs is to formulate and solve the 

differential equations of motion (EOM). Given that N is typically small, ∑F mz i ii =  can be 

solved, allowing us to easily keep track of velocities, forces, energies, and positions. The 

problem will consist of TCs barely touching in their initial configuration and affixed between 

rigid walls or spheres of infinite radius. 
12 

Figure 5. KE as a function of N,q and f for DTC hard-spheres. Note that the colorscale is calibrated for each 

subplot. 

The always-repulsive potential between adjacent grains i and i +1 is due to geometric effects 

and derived in Landau and Lifshitz (23). It can be written as 

( ) 1 52 52

,1 ,1 ,1 ,1

1

2

5

i i

ii ii ii ii

i i

R R V a

DR R

δ δ δ +

+ + ++

+

= ≡ +

, (45) 

where Rj : , 1, , j ii N = + { " } are the radii and 

2 2

1

1

3 1 1

4

i i

i i

D

E E

σ σ +

+

⎛ ⎞ − − = + ⎜ ⎟ ⎝ ⎠

. (46)
13 

Here, δi,i+1 = Ri + Ri+1 − zi+1 − z ( i) represents the overlap between successive grains, where zj is 

their position. Additionally, i i, 1 a + has been defined for material properties: Ej , the Young’s 

modulus and σ j , the Poisson’s ratio; and radii, Rj . Note also that j can refer to either particle i

or i + 1. If δi,i+1 ≤ 0, then V = 0 since adjacent grains i and i +1 begin to lose contact. For our 

particular study, all materials are the same in any single TC, so it reduces to 

2 3 1

2

D

E

⎛ ⎞ −σ = ⎜ ⎟ ⎝ ⎠ . (47) 

Under the Hertzian potential, equation 45, a squeezed sphere initially appears soft. As 

compression increases, its stillness increases substantially. In this sense, a Hertz potential can be 

thought of as a nonlinear spring. 

At the boundaries, the enclosing walls are represented by spheres of infinite radius: R0 and RN +1. 

The corresponding potentials then for grains against the boundary scale as 

V( ) δ wall,1 ∝ R1 and V( ) δ N,wall ∝ RN . 

The force on grain i, therefore, is the sum of influences from its neighboring grains: 

( ) 32 32

1, 1, , 1 , 1

5

2 i i i i i i ii ii

dV

mz a a

d

δ δ

δ == − −− ++  . (48) 

These equations are solved numerically using the Velocity-Verlet algorithm (27). The original 

source code (8) which applies to the STC is listed in appendix A. Modifications to this for the 

DTC are listed in appendix B. 

Results have been obtained for a large selection of chains consisting entirely of Ti6Al4V or SiC 

spheres. Arbitrarily, we have chosen to use Ti6Al4V when restitutive losses are ignored (ω = 0) 

and SiC otherwise. The following material properties were assumed (reference D in 

equation 47). 

Table 1. Material properties (28). 

Material ρ

(mg/mm3

) 

D

(mm2

/kN) 

Occurrence 

SiC 3.2 0.003266 ω ≠ 0

Ti6Al4V 4.42 0.01206 ω = 0

Note that the fundamental unit of force in our simulations is the kilonewton and an initial 

velocity of 0.01 mm/μs is applied to the largest end of each chain. 

Simulations were performed for 3 ≤ N ≤ 20 and 0 ≤ q,ω ≤ 0.1 over a system time of 1 ms where 

the timestep was set to 10 ps, corresponding to 108

 steps in the integration loop.
14 

Restitution, the mechanism for introducing energy losses, is tuned by an asymmetry between the 

(contact) forces of loading and unloading during a collision: 

1 unloading

loading

F

F = −ω . (49) 

Therefore, the unloading or expansive force in a collision is some fraction of the loading force 

and perfectly elastic collisions correspond to ω = 0. 

In order to handle the hundreds of TC calculations, a Practical Extraction and Report Language 

(PERL) script (appendix C) was created to automate the process. It consists of nested loops that 

iterate through N and q inserting their new values into a template containing the C code which 

solves the EOM of the TCs. It then copies this to a directory labeled by some appropriate 

mnemonic. This was done so that when visualizing the data in MATLAB, it would be easy to 

change directories as part of a loop. An example directory structure is, 

/w0/N20/N20q10/taperchain1.cpp. Data files are created for every grain in each TC 

and their number is encapsulated in the filename. Appendix D contains the MATLAB code to 

generate KE surfaces. 

Recall that we are interested in forming the ratio KEN = KEout KEin , the normalized KE. Here, 

KEout is the first* kinetic energy peak felt by the last (smallest) grain and KEin is a constant 

determined from the initial velocity and density of the first grain. Conversely, Fout is the first 

minimum felt by the last grain since the direction of a negative force is into the wall or force 

sensor. The algorithm to pick out the first turning point for each of these is straightforward. 

Since KE t( ) is a column vector, one can iterate through each element until the first occurrence of 

( ) i +1 < i is true. In that case, i represents the extremum. For force, we compare against 

( ) i +1 > i. The only data files then that are of interest for KEN are those that represent the 

dynamics of the first (KEin ) and last particle (KEout) in each chain. 

4.1 STC 

Figure 6 plots a small cross section of the results for several chains withN = 20, ω = 0.05, 

andq = { } 0, 0.02, 0.04, 0.06, 0.08, 0.1 . These plots show the kinetic energy spectrum of the 

tail particle as a function of time for the first half of the simulation (500 μs). From equation 49, 

each collision only transmits 95% of the force. Note the decreasing scale of KE in each 

underlying panel. One might expect KE to increase since particle velocity is increasing. 

However, KE scales as v

2

 and m ∝ r

3 ∝(1− q)

3

, so the latter term dominates, and consequently 

KE decreases. 

In the top panel of figure 6, we see at about 95 μs that the tail particle in a monodisperse chain 

first receives the energy bundle and undergoes translational motion. Its contact force with the

 *

 Normalization is discussed further in appendix E. 
15 

Figure 6. KE as a function of time for the smallest particle in a chain with N = 20, ω = 0.05, and 

q = [ ] 0 : 0.02 : 0.1 . Each plot is evaluated for a different tapering with initial KE as 0.0838 J.
16 

boundary increases as it slows down eventually coming to rest with all energy stored as elastic 

potential at about 100 μs. This is then converted back to kinetic as the tail particles begin to 

rebound. 

By selecting the first peak in each of these subplots as well as the many other chains not shown 

here, we can form the normalized KE surfaces and observe how well STCs perform as energy 

absorbers. 

Pulses in monodisperse chains (N ≥15, q = 0) represent a special class of waves known as a 

solitary waves which are remarkably stable and have been studied extensively (2–7, 9, 11–12). 

For TCs ( q ≠ 0), the solitary wave cannot quite fully form because subsequent grains within the 

wave envelope move at higher, disparate velocities with increasingq—leading to greater 

dispersion and destructive wave broadening. This translational symmetry breaking was also 

reported by and is in agreement with Nakagawa et al. (4). 

4.1.1 STC KE and Force Parameter Spaces: KE N,q, ( ω) , F N,q, ( ω)

Figure 7 highlights the numerical results for KE(N, q, constant ω). Note that it would be very 

difficult and costly to produce a similar plot empirically. Each node on the surface represents a 

different TC and experiment. There appears to be a sigmoidal and exponential dependence on q

andN, respectively, in agreement with the hard-sphere approximation. These KE surfaces 

represent STC chains that thermalize more than half the incident energy introduced into the 

system. For example, the least effective shock mitigating geometry that we've simulated here—

KE(ω = 0, N = 3, q = 0)—reduces the output KE by about 60%. That amount, after restitution is 

accounted for in cases where ω ≠ 0, is distributed among the other grains as kinetic and potential 

energy. Note again that results are independent of initial grain size. For monodisperse chains 

( q = 0) in figure 7a, there is asymptotic behavior as N increases which is a result of energy 

propagating as solitary waves, when N is of sufficient number (about 15), energy can propagate 

with negligible loss. Additionally, there is a rapid drop in KE(3 ≤ N ≤ 5, q = 0, ω = 0) which is 

likely due to the width of the incident pulse being longer than the TC. This implies that wave 

reflection begins before the impulse is fully applied. 

A similar surface can be plotted representing the normalized force for various TCs. This is 

visible in figure 8 and similar to the KE plots. Again, the functionality of N here is essentially a 

one to two phase exponential, but the sigmoidal nature of q is less obvious and can be 

approximated linearly throughout the parameter space. In general, the magnitude of F appears 

to be twice that of KE . 

The numerical results can be compared to the hard-sphere approximation by scaling the latter to 

the largest numerical solution—KE (N = 3, q = 0)—and plotting their difference (figure 9). 

Energy loss is not accounted for in this comparison. The distribution is due to differences in the 

interaction potential and width of the energy bundle. Therefore, it appears that for smaller 

particles (i.e., larger q) compressive effects and width of the energy packet becomes increasingly 

important. This is compounded further by increasing the number of particles, N.
17 

Figure 7. Numerical solution of KE N( ) ,q parameter space for constant ω. Note that it would be very difficult 

and costly to produce a similar plot empirically. Each node on the surface represents a different TC and 

experiment.
18 

Figure 8. Numerical solution of F N( ) ,q parameter space for constant ω. 
19 

Figure 9. Difference plot of the lossless KE parameter spaces for the hard-sphere 

approximation (figure 2a) and numerical solution (7a). Note that the azimuthal view 

has been rotated 180° for clarity. 

4.1.2 STC Mathematical Model 

It would be convenient to have a predictive capability for at least a selection of individual TCs as 

a function of N, q, and ω. We therefore propose a mathematical model for the 

KE = KE( ) ω,q,N = 20 parameter space of the form 

( ) , C E Bq D KE q Ae e ω ω = , (50) 

which corresponds to a two-dimensional Weibull distribution. One can obtain the sigmoidal and 

exponential behavior of q and ω, respectively, when the exponents are restricted to C >1 

and E ≤1. For simplicity, we set E =1 and C = 3 2 since these values provide the general 

functionality. The coefficients B and D were evaluated using 4th-order polynomial fits and the 

scaling coefficient A is essentially the point KE N( = 20,ω = q = 0) determined in the simulation. 

It turns out that a 2nd-order fit was not sufficiently robust and a higher-order fit yielded marginal 

gains at the cost of mathematical encumbrance. 

This model currently lacks the rigor sufficient for planar behavior in the limit of small N. It is 

unclear at this point if N can be completely decoupled from ω and q and written as an additional 

exponential term. In all likelihood, the coefficients B and D would be written as functions of N. 
20 

Ultimately, it would be useful to expand equation 50 and have a model that completely describes 

the normalized KE as a function of N, q, ω, and ∆ where the latter is a constant external loading 

or precompression of the chain. 

In evaluating the fits, we find that KE(ω, q, N = 20) is described quite well by 

( ) ( )

( )

5 4 4 3 3 2 32

54 33 2

 1.5055 10 4.016 10 3.98110 0.0147 0.05435

 4.144 10 1.955 10 0.03962 0.02887 8.341

, 0.35544

 

q q qq q

q qqq

KE q e

e ω

ω − −−

− −

⎡ ⎤ − ⋅ + ⋅ −⋅ + − ⎣ ⎦

⎡ ⎤ ⋅ +⋅ − + − ⎣ ⎦

=⋅ ⋅

⋅

 (51) 

It is compelling to rewrite equation 51 in the more suggestive form 

( ) 11 2 9 2 7 2 5 2 3 2

, q qqqq KE q e e e e e α βγδε ω ∝ " (52) 

( ) ( ) () () n n n nn 3 21 1 q q q qq e e e ee

α β γ δε + ++ −

∝ " (53) 

where the powers of q have been adjusted to be written as a series in n, the interaction 

potential—equivalent to 2.5 for Hertzian spheres. The occurrence of an apparent series in n is 

quite striking. Missing from this of course is a (n − 2) or q

1 2 term which requires a different fit 

for B in equation 50. Note also the mixed term qω in equation 51. This is indicative of a manybody effect; one cannot have restitution in a single particle chain since, by definition, it requires 

at least two to evaluate the loading and unloading. 

Figure 10 compares the numerical results and the model for the case of N = 20. The results are 

practically indistinguishable except for minor timestep errors in the simulation which appear as 

noise for large q. At these values for N and q, the tail particle has a large velocity and for 

timesteps not small enough, the KE plots lose their smoothness. The algorithm therefore selects 

a peak close to what the real extremum would be in the limit that the timestep approaches zero. 

As such, the model is more accurate. 

4.2 DTC 

In this section, we address the more capable DTC. KE time spectra is excluded since an analysis 

of the complicated motion is outside the scope of this report. Instead, let us focus immediately 

on the KE parameter space and compare with the STC. 

Figure 11 illustrates the intriguing nature of the DTC KE N( ,q, f) surfaces where the 

dissimilarity from the DTC hard-sphere approximation is quite clear. Immediately obvious is a 

ripple across the surface for large f which translates towards the origin as f decreases. The effect 

vanishes at about f = 0.6 and the KE surface resembles that of a much-improved STC for

f ≤ 0.5. We propose qualitatively that this is a consequence of the inertial mismatch changing 

as a function of position in the DTC—in contrast to the STC where it is constant along the chain 

and the wave effect is not observed. 
21 

Figure 10. The simulated and modeled STC. 

In the side panel in figure 11, several DTCs are shown. In the top case, f = 1, so that as the 

incident pulse propagates to the left, the chain becomes increasingly monodisperse. That degree 

of monodispersity is a function of N and q since wave amplitude and velocity is sensitive to 

them. If N is large enough, then whatever remaining energy is left could propagate nearly 

lossless like a solitary wave. As f gets smaller (next panel down), any apparent monodispersity 

begins to disappear. As one moves to the right in the chain, the interstitial grains become small 

compared to the neighboring grains. As f becomes smaller still, the interstitial grain is small 

compared to all grains such that reducing f has little effect—little enough so that the wave 

behavior on the KE surface is not observed. Thus there appears to be a transition so that 

when f ≤ 0.5, the chain will never appear monodisperse and inertial mismatches can always be 

considered large.
22 

Figure 11. DTC numerical solution of KE N( ,q, f) and sample DTCs for various f (inset). The scale is the same 

as that in figure 7. 

5. Experimental Efforts 

Up to this point, most experimental work reported in the open literature has focused on 

monodisperse chains by themselves or as a precursor to TCs. Here we only discuss the setup 

used in such studies as well as the prototype for our intended investigation of the DTC. 

5.1 STC 

There are two common ways of performing STC experiments. The first (11) is based on a 

Newton's cradle but more epic in extent (figure 12). In this particular case, there is a 

monodisperse sequence used to set up a well-defined solitary wave which propagates to the left. 

A small striking sphere on the far right is applied to generate the initial pulse. 
23 

Figure 12. STC experimental setup by Job et al. (12).

The second common STC configuration (4) consists of a V-shaped track which can be slanted to 

ensure proper (centerline) contact of the tapered spheres (figure 13). Here, the left column 

highlights a monodisperse chain and the right column shows the same experiment but for a TC—

subsequent rows represent later time. In the middle row, the striker (gold sphere) has hit the 

chain; in the bottom row, the tail particle in the monodisperse chain has perforated the paper 

target (black hole in upper left of picture), while for the TC on the right, it merely dents and 

bounces off the paper. This is a remarkable demonstration of the impulse decimation capability 

inherent to TCs. 

5.2 DTC 

Figures 14–16 highlight the constituents and current housing of a DTC configuration. 

Within the assembly, interstitial grains are radially constrained along the axis of the DTC by thin 

discs with small holes bored in them (figure 16). The apparatus can be tested in any orientation; 

however, when vertical, asymmetric loading due to gravity will become a factor. Uniform 

precompression enhances the shock absorption capability of TCs and will be discussed in the 

next section. In figure 15, we see that precompression of the DTC prototype can be controlled 

by the large screw at the top of the assembly. 
24 

Source: Picture courtesy of Prof. M. Nakagawa of the Colorado School of Mines and Dr. Juan Agui of the National Aeronautics 

and Space Administration–Glenn Research Center. 

Figure 13. STC experimental setup. The left half of the figure represents a monodisperse chain while the right a 

TC. Lower plots denote later time. 
25 

Figure 14. Preliminary DTC housing chamber. 
26 

Figure 15. DTC constituents. 
27 

Figure 16. DTC interstitial grains. 

6. Precompressed Chains 

Apparently, when TCs are precompressed, shock absorption increases substantially (3). We 

measure precompression, ∆, as a percentage of particle overlap in the initial state. The EOM in 

equation 48 are modified accordingly and become 

{ [ ] ( ) [ ] ( ) } 32 32

1, 1 1 , 1 1 1

5

2

mz a R R z z a R R z z i i i i i i i i ii i i i i = + −∆− − − + −∆− − −− − + + +  . (54) 

Figure 17 highlights the impressive results for KE N( = 20, ∆ = 0.03%) as tapering and 

restitution are varied. Precompression in a DTC could therefore provide the best shock 

absorption geometry that we have currently considered.
28 

Figure 17. KE as a function of tapering and restitution for 

constant N = 20 and ∆ = 0.03% . 

7. TC Armor Panels 

We have seen that Hertzian TCs act as shock absorbers. In particular, recall that at the end of 

section 3.1.3 we noted that TCs operate regardless of a pulse's direction through the chain. Thus, 

to maximize efficiency, we can place neighboring STCs in alternating directions. 

For them to become realizable in armor applications, however, one should evaluate them based 

on performance of the whole system. The specific absorbed energy (SAE)—which has units of 

J/g—is one such metric and commonly used by ballisticians. It is quantified by equation 55, 

while the volumetric absorbed energy (VAE) is described by equation 56 and has units of J/cm3

. 

( ) 1

1 ,,

out

in

in out in in

KE KE

KE KE KE KE N q KE

MM M

ω

⎛ ⎞ ⎜ ⎟ − − ⎝ ⎠ − = = . (55) 

 1 ,, ( ) KE N q KEin

V

− ω

. (56)
29 

Here, KE N( ) ,q,ω is the normalized KE that has been evaluated in the numerical simulations. 

Figure 18 illustrates what a TC armor prototype might resemble along with an incident flyer 

plate. 

Figure 18. STC armor plate panel. 

In equations 55 and 56, KEin represents the energy of the incoming flyer plate. For convenience, 

we assume that all materials are Ti6Al4V ( ρ = 4.42 mg/mm3

), that uf = 100 m/s, tf = 24.5 mm, 

and that l = w = 230 mm. Since m = ρV = ρlwtf , the KE of the flyer plate, KEin , is 28.7 kJ. 

Also, in both cases, we set KE N( ) ,q,ω = KE(20,0.1,0.06) ≈ 0.05, which is obtained from the 

numerical results. Thus, the sizes of the end spheres in such a chain are d20 =10 mm and 

d1 =1.35 mm. 

In order to evaluate M in equation 55, several assumptions are made and then later relaxed to 

determine a range for M and better approximate what type of performance a TC armor panel 

might achieve. Figure 18 displays the relevant variables and configuration for a preliminary 

analysis, where t = 95.8 mm is chosen based on the chain parameters used in picking 

KE N( ) ,q,ω . 

Each STC would require some type of mechanical support to maintain the overall geometry. 

One solution is to house them in hollow cones as is displayed in the diagram. Conical springs, or 

Belville washers, are another possibility. We can absorb this added material into d1 as an 

effective diameter which we set to  ̃ d 1 =1.5 mm. The additional dimensions, w,l, are arbitrarily 

assigned to be a function of the number of TCs: l = w = 20 d20 +  ̃ d ( 1)= 230mm, so that each side 

is spanned by 40 TCs. Therefore, the total number of TCs in the armored panel is 1600. Finally, 
30 

δt

1,δt2 represent thicknesses for enclosing plates so that the TCs aren't exposed. This amount is 

arbitrary but chosen to be small, δt

1 = δt2 = 2.5 mm, so that the total thickness of the panel is 

100.8 mm, or about 4 in. 

Now, M can be written as the sum, M = NSTCmSTC + 2mδt + msupport + mfill, where the first term is 

the product of the mass of one STC and the total number of them, the second is the contribution 

from the enclosing plates, the third is due to the conical housing of STCs, and the last is the filler 

material which we ignore as being air. Using conservative values, we let mδt = ( ) δt lwρ = 0.58

kg, msupport =1.0 kg so that M =13.63 kg. 

Since equation 55 scales with uf , the velocity of the flyer plate, a wide range of SAE can be 

determined. If uf =100 m/s, the SAE=2. Let us assume then that 6.63 ≤ M ≤ 20.63 and 

100 ≤ uf ≤ 300 m/s and plot the SAE surface which is visible in figure 19a. 

The VAE can be calculated by dividing this quantity by the density of the material—since we are 

assuming the whole panel is Ti6Al4V (figure 19b). Another measure, however, would be to 

divide by the volume of the entire panel since that space is already claimed and cannot be 

occupied by anything else. In dividing by the density, we are only counting that portion of the 

armor panel contained by Ti6Al4V; it doesn't include the air fills. Note that the SAE and VAE 

have been determined primarily analytically. 

Figure 19. SAE and VAE as a function of M, the total STC armor panel mass, and uf , the incident flyer plate 

velocity. 
31 

8. Concluding Remarks 

This report has focused on TCs and their inherent ability to decimate propagaing impulses. Due 

to a lack of dependence on system size, it appears that TCs can provide meter to submicron scale 

shock absorption. TCs take well-defined pulses, such as shocks, and turn them into noise. The 

signal is spread out over time and space and the process is known as thermalizing an impulse (3). 

Essentially, the shock is stretched from a few microseconds to several hundred with its extent 

increasing from several grains to all members in the TC. This is in opposition to monodisperse 

chains which act as highly efficient shock transmitters due to the resilient nature of solitary 

waves. In fact, Pöschel and Brilliantov have calculated (29) the optimal transmission of KE for 

TCs. 

It turns out that regardless of whether an impulse travels forwards or backwards along a chain, 

both directions mitigate a propagating pulse. At later times, this likely competes with 

coordinated motion of adjacent grains which may create unpredictable energy or force spikes for 

certain configurations, as figure 7 (lower panel) demonstrates. Antiparallel orientations of TCs 

can then be arranged to exploit the collective effects of many TCs as an armor panel. The 

capabliity of such a panel is analytically derived and yields values for the SAE on the same order 

100

-102 J/g as other technologies. However, such a determination is still rather suggestive 

without a full scale prototype. TCs may be housed in hollow cones or Belville washers (conical 

springs) for added absorption. Additionally, there is evidence (30) that nitinol (a nickel-titanium 

alloy) may be a better material to choose for the spheres and housing. 

Hard-sphere approximations were analytically obtained for two chain geometries, the STC and 

DTC. This approximation simplifies the problem to one where the impulse delivered is a single 

grain in width and elastic effects are ignored. For the STC, this approximation is comparable to 

the numerical results. However, for the DTC, the differences are quite large and consequently 

elastic effects play a major role in the system. This is likely due to the increased velocity of 

members in the chain and therefore the number of collisions. 

The authors see ample opportunity for study along several lines to move this technology forward 

and rate it as a useful technology in blast mitigation system. Quasistatic compression to 

plasticity and failure, ALEGRA simulations, and measuring the effects of many TCs working in 

cooperation are areas we are beginning to investigate. 
32 

9. References 

1. Doney, R.; Sen, S. The Tapered Chain: Hard Sphere Approximation and Simulation. Phys. 

Rev. E., submitted for publication, July 2005. 

2. Sokolow, A.; Pfannes, J. M. M.; Doney, R.; Nakagawa, M.; Agui, J. H.; Sen, S. Impulse 

Absorption by Small Tapered Granular Chains. Science, submitted for publication, 2005. 

3. Sen, S.; Manciu, M. Thermalizing an Impulse. Phys. A 2001, 299, 551–558. 

4. Nakagawa, M.; Agui, J.; Wu, D. T.; Extramiana, D. V. Impulse Dispersion in a Tapered 

Granular Chain. Gran. Mat. 2003, 4, 167–174. 

5. Wu, D. Conservation Principles in Solitary Impulse Propagation Through Granular Chains. 

Phys. A 2002, 315, 194–202. 

6. Manciu, M.; Sen, S.; Hurd, A. J. The Propagation and Backscattering of Soliton-Like Pulses 

in a Chain of Quartz Beads and Related Problems (I). Propagation. Phys. A 1999, 274, 588–

606. 

7. Manciu, M.; Sen, S.; Hurd, A. J. The Propagation and Backscattering of Soliton-Like Pulses 

in a Chain of Quartz Beads and Related Problems (II). Backscattering. Phys. A 1999, 274, 

607–618. 

8. Pfannes, J. Energy Propagation in Granular Chains. M.S. Thesis, State University of New 

York, Buffalo, NY, May 2003. 

9. Coste, C.; Falcon, E.; Fauve, S. Solitary Waves in a Chain of Beads Under Hertz Contact. 

Phys. Rev. E 1997, 56 (5), 6104–6117. 

10. Rossmanith, H. P.; Shukla, A. Photoelastic Investigation of Dynamic Load Transfer in 

Granular Media. Acta Mech. 1982, 42, 211–225. 

11. Sen, S.; Manciu, M.; Wright, J. D. Solitonlike Pulses in Perturbed and Driven Hertzian 

Chains and Their Possible Applications in Detecting Buried Impurities. Phys. Rev. E 1998, 

57 (2), 2386–2397. 

12. Job, S.; Melo, F.; Sokolow, A.; Sen, S. How Hertzian Solitary Waves Interact with 

Boundaries in a 1D Granular Medium. Phys. Rev. Let. 2005, 94 (178002), 1–4. 

13. Duran, J. Sands, Powders, and Grains; Springer-Verlag: New York, NY, 2000. 

14. Mehta, A., Ed. Granular Matter; Springer-Verlag: New York, NY, 1994. 
33 

15. Kakalios, J. Resource Letter GP-1: Granular Physics or Nonlinear Dynamics in a Sandbox. 

Am. J. Phys. 2005, 73 (1), 8–22. 

16. Jaeger, H. M.; Nagel, R. Physics of the Granular State. Science 1992, 255 (5051), 1523–

1531. 

17. Jaeger, H. M.; Nagel, R.; Behringer, R. P. Granular Solids, Liquids, and Gases. Rev. Mod. 

Phys. 1996, 68 (4), 1259–1273. 

18. Wright, T. Weak Shocks and Steady Waves in a Nonlinear Elastic Rod or Granular 

Material; BRL-TR-02584; U.S. Army Ballistics Research Laboratory: Aberdeen Proving 

Grounds, MD, 1984. 

19. Walton, O. R.; Braun, R. L. Viscosity, Granular-Temperature, and Stress Calculations for 

Shearing Assemblies of Inelastic, Frictional Disks. J. Rheol. 1986, 30, 949–980. 

20. Stronge, W. J. Impact Mechanics; Cambridge University Press: Cambridge, MA, 2000. 

21. Frey, R.; Gniazdowski, N.; Li, T. H.; Tarzian, F. Ballistic Shock from Reactive Armor and 

Explosive Launchers; U.S. Army Research Laboratory: Aberdeen Proving Ground, MD, to 

be submitted for publication, 2005. 

22. Hertz, H. Über Die Berührung Fester Elastischer Körper (On the Behavior of Solid Elastic 

Bodies). J. Reine Angew. Math. 1882, 92, 156–171. 

23. Landau, L. D.; Lifshitz, E. M. Theory of Elasticity; Pergamon Press: Oxford, 1970; 30–36. 

24. Love, A. E. H. A Treatise on the Mathematical Theory of Elasticity; Dover: New York, NY, 

1955. 

25. Leroy, B. Collision Between Two Balls Accompanied by Deformation: A Qualitative 

Approach to Hertz’s Theory. Am. J. Phys. 1984, 53 (4), 346–349. 

26. Nesterenko, V. F. Dynamics of Heterogeneous Materials; Springer-Verlag: New York, NY, 

2001. 

27. Allen, M.; Tildesley, D. Computer Simulation of Liquids; Oxford: New York, NY, 1987. 

28. MatWeb Home Page. http://www.matweb.com (accessed August 2004). 

29. Pöschel, T.; Brilliantov, N. V. Extremal Collision Sequences of Particles on a Line: 

Optimal Transmission of Kinetic Energy. Phys. Rev. E 2001, 63 (021505), 1–9. 

30. Jackson, C. M.; Wagner, H. J.; Wasilewski, R. J. 55-Nitinol—The Alloy Tith a Memory: Its 

Physical Metallurgy, Properties, and Applications; NASA-SP-5110; National Aeronautics 

and Space Administration: Washington, DC, 1972. 
34 

INTENTIONALLY LEFT BLANK. 
35 

Appendix A. Simple Tapered Chain Code 1

 

 This appendix appears in its original form, without editorial change. 

1

 Pfannes, J. Energy Propagation in Granular Chains. M.S. Thesis, State University of New York, Buffalo, NY, May 2003. 
36 

/* PROGRAM taperchain.cpp 

 This program consideres an one dimensional chain of spheres that 

 shrink succesively in radius ("tapered chain"). Initially the spheres 

 are barely in contact, i. e. they just touch each other and are not 

 compressed (zero loading). 

 The chain ends at both edges at fixed walls. 

 The program calculates the interaction of the system once disturbed 

 by an instantaneous (delta) impulse exerted on one end of the chain. 

 Restitution both between the spheres and between the edge spheres 

 and the corresponding wall can be introduced. 

 The program does not consider gravity. 

 The EOM are solved with the Velocity Verlet algorithm. 

 scale of problem: mm-mg-musec (mimimu) 

 in this scale unit of force: 1000 N 

 in this scale unit of energy: 1 J 

 */ 

#include <cmath> 

#include <iostream> 

#include <fstream> 

#include <cstdlib> 

#include <string> 

#include <sstream> 

using namespace std; 

/************************** ALTERABLE PARAMETER: ****************************/ 

const int nptles=20; // total number of particles 

const double rho=3.2 /* SiC (mg/mm^3) */, D=0.00326603139013 /* (mm^2/N) */; 

const double rlarge = 5.0; // (radius of large ptle (mm)) 

const double q = 0.0; // (tapering factor (%)) 

const double xn = 2.5; // (exponent in potential) 

const double dt = 0.00001; // (timestepwidth (musec)) 

const unsigned int nsteps = 100000000; // (# steps integration loop) 

const int idiagp = 20000; // (stepwidth diagnostics) 

const int idump = 20000; // (stepwidth dump) 

const double v1in = 0.0; // (initial v small ptle (mm/musec)) 

const double vnin = -0.01; // initial v large ptle (mm/musec)) 

const double epsilon = 1.0; // ((1 - restitution factor) all ptles) 

/****************************************************************************/ 

ofstream readme("taperchain.readme"); // global scope fcts 

ofstream EnergyImpulse("taperchain.EneImp"); 

void radii (double rlocal[]) { 

 rlocal[nptles-1] = rlarge; 

 if (q == 0) // avoid roundoff errors w/out tapering 

 for (int i = 0; i < nptles-1; i++) 

 rlocal[i] = rlarge; 

 else 

 for (int i = 2; i < nptles+1; i++) 
37 

 rlocal[nptles-i] = rlocal[nptles-i+1] * (1 - q*0.01); 

} 

void masses (double r[], double masslocal[]) { 

 const double pi = 4 * atan(1.0); 

 const double masslarge = (4.0/3.0) * pi * rlarge*rlarge*rlarge * rho; 

 masslocal[nptles-1] = masslarge; 

 if (q == 0) // avoid roundoff errors w/out tapering 

 for (int i = 0; i < nptles-1; i++) 

 masslocal[i] = masslarge; 

 else 

 for (int i = 0; i < nptles-1; i++) 

 masslocal[i] = r[i]*r[i]*r[i] * masslarge / (rlarge*rlarge*rlarge); 

} 

void strenghtfac (double r[], double alocal[]) { 

 alocal[0] = (2.0 / (5.0 * D)) * (sqrt(r[0])); 

 alocal[nptles] = (2.0 / (5.0 * D)) * (sqrt(r[nptles-1])); 

 if (q == 0) // avoid roundoff errors w/out tapering 

 for (int i = 1; i < nptles; i++) 

 alocal[i] = (2.0 / (5.0*D)) * (sqrt(0.5*rlarge)); 

 else 

 for (int i = 1; i < nptles; i++) 

 alocal[i] = (2.0 / (5.0 * D)) * (sqrt((r[i]*r[i-1])/(r[i]+r[i-1]))); 

} 

// initialpos prints absolute initial positions, not for calculations 

void initialpos (double r[], double xInitiallocal[]) { 

 if (q == 0) // avoid roundoff errors w/out tapering 

 for (int i = 0; i < nptles; i++) 

 xInitiallocal[i] = (2.0*(i+1) - 1) * rlarge; 

 else { 

 xInitiallocal[0] = r[0]; 

 for (int i = 1; i < nptles; i++) 

 xInitiallocal[i] = xInitiallocal[i-1] + r[i-1] + r[i]; 

 } 

} 

// absolutpos prints absolute positions to ptle files, not for calculations 

void absolutpos (double r[], double x[], 

 double xInitial[], double xAbsolutlocal[]) { 

 for (int i = 0; i < nptles; i++) 

 xAbsolutlocal[i] = xInitial[i] + x[i]; 

} 

void computeAccelerations (double x[], double a[], double r[], 

 double acc[], double overbefore[], 

 double mass[], double& pot) { 

 pot = 0.0; // every call calculates new pot contributions 

 for (int i = 0; i < nptles; i++) // zeroing all acc in every call 

 acc[i] = 0.0; // (= every timestep) 

 /******* potential/force between neighboring ptles *****************/ 

 for (int i = 0; i < nptles-1; i++) { 

 if (x[i] > x[i+1]) { // only when overlap 
38 

 double over = x[i] - x[i+1]; 

 double overnm1 = pow(over, (xn - 1.0)); 

 pot += over * overnm1 * a[i+1]; 

 double forceBetw = a[i+1] * xn * overnm1; 

 double forceFactor; 

 if (overbefore[i+1] < over) // when compressing 

 forceFactor = 1.0; 

 else forceFactor = epsilon; // when decompressing 

 forceBetw *= forceFactor; 

 // dim acc: force 

 acc[i] -= forceBetw; // sign(-): towards smaller x 

 acc[i+1] += forceBetw; // sign(+): towards larger x 

 overbefore[i+1] = over; // update for next timestep 

 } 

 else overbefore[i+1] = 0.0; // reset when no overlap 

 } 

 /** potential/force between fixed wall (small, x=0) <-> small ptle **/ 

 if (x[0] < 0) { 

 double over = - x[0]; 

 double overnm1 = pow(over, (xn - 1.0)); 

 pot += over * overnm1 * a[0]; 

 double forceSmall = a[0] * xn * overnm1; 

 double forceFactor; 

 if (overbefore[0] < over) 

 forceFactor = 1.0; 

 else forceFactor = epsilon; 

 forceSmall *= forceFactor; 

 acc[0] += forceSmall; 

 

 overbefore[0] = over; 

 } 

 else overbefore[0] = 0.0; 

 

 /*** potential/force between fixed wall (large) <-> large ptle ******/ 

 if (x[nptles-1] > 0) { 

 double over = x[nptles-1]; 

 double overnm1 = pow(over, (xn - 1.0)); 

 pot += over * overnm1 * a[nptles]; 

 double forceLarge = a[nptles] * xn * overnm1; 

 double forceFactor; 

 if (overbefore[nptles] < over) 

 forceFactor = 1.0; 

 else forceFactor = epsilon; 

 forceLarge *= forceFactor; 

 acc[nptles-1] -= forceLarge; 

 overbefore[nptles] = over; 
39 

 } 

 else overbefore[nptles] = 0.0; 

 /***** real dim of acc: division by mass **********/ 

 for (int i = 0; i < nptles; i++) 

 acc[i] /= mass[i]; 

} 

void velocityVerletStep (double x[], double v[], double acc[], 

 double a[], double r[], double overbefore[], 

 double mass[], double& pot) { 

 for (int j = 0; j < nptles; j++) { 

 x[j] += v[j] * dt + 0.5 * acc[j] * dt*dt; 

 v[j] += 0.5 * acc[j] * dt; 

 } 

 computeAccelerations (x, a, r, acc, overbefore, mass, pot); 

 for (int j = 0; j < nptles; j++) 

 v[j] += 0.5 * acc[j] * dt; 

} 

void ptleHeader (ofstream* print, int k) { 

 (* print) << "# ptle " << k+1 << ": time (musec)" << '\t' << "x (mm)" 

 << '\t' << "v (mm/musec)" << '\t' << "a (mm/musec^2)" 

 << '\t' << "kin. E. (J)" << '\t' << "f (kN)" << '\t' 

 << "impulse (mg*mm/musec)" << '\t' << "xRelative (mm)" << '\n'; 

} 

void dumpData (double t, double mass[], double v[], double acc[], 

 double r[], double x[], // scope absolutpos 

 double xInitial[], double xAbsolut[], 

 ofstream* print, int k) { 

 double keDumpPt[nptles]; // new arrays for dumping data 

 double vDumpPt[nptles]; // since arrays pass by argument 

 double accDumpPt[nptles]; // dump data manipulated 

 double xDumpPt[nptles]; 

 keDumpPt[k] = 0.5 * mass[k] * v[k]*v[k]; 

 if (keDumpPt[k] < 1.0e-20) keDumpPt[k] = 0.0; // set small values to zero 

 

 vDumpPt[k] = v[k]; 

 if (vDumpPt[k] < 1.0e-20 && vDumpPt[k] > -1.0e-20) vDumpPt[k] = 0.0; 

 

 accDumpPt[k] = acc[k]; 

 if (accDumpPt[k] < 1.0e-20 && accDumpPt[k] > -1.0e-20) accDumpPt[k] = 0.0; 

 

 xDumpPt[k] = x[k]; 

 if (xDumpPt[k] < 1.0e-20 && xDumpPt[k] > -1.0e-20) xDumpPt[k] = 0.0; 

 absolutpos (r, x, xInitial, xAbsolut); // calculate absolute pos. 

 (* print) << t << '\t' << xAbsolut[k] << '\t' << vDumpPt[k] << '\t' 

 << accDumpPt[k] << '\t' << keDumpPt[k] << '\t' 

 << accDumpPt[k]*mass[k] << '\t' 
40 

 << mass[k]*vDumpPt[k] << '\t' << xDumpPt[k] << '\n'; 

} 

void dumpEnergyImpulse (double t, double kelocal, double telocal, double pot, 

 double ptotallocal, double mass[], double v[]) { 

 kelocal = 0.0; 

 ptotallocal = 0.0; 

 double absptotallocal = 0.0; // scope only within function 

 for (int j = 0; j < nptles; j++) { 

 kelocal += mass[j] * v[j]*v[j]; 

 ptotallocal += mass[j] * v[j]; 

 absptotallocal += mass[j] * abs(v[j]); 

 } 

 kelocal *= 0.5; 

 telocal = kelocal + pot; 

 double potDump = pot; 

 if (kelocal < 1.0e-20) kelocal = 0.0; // set very small values to zero 

 if (ptotallocal < 1.0e-20 && ptotallocal > -1.0e-20) ptotallocal = 0.0; 

 if (telocal < 1.0e-20) telocal = 0.0; 

 if (potDump < 1.0e-20) potDump = 0.0; 

 if (absptotallocal < 1.0e-20) absptotallocal = 0.0; 

 EnergyImpulse.precision(16); 

 EnergyImpulse << t << '\t' << kelocal << '\t' << potDump << '\t' 

 << telocal << '\t' << absptotallocal << '\t' 

 << ptotallocal << '\n'; 

} 

void readmeInfo (double ke, double pot, double te, double ptotal, 

 double ke1in, double kenin, double p1in, double pnin, 

 double r[], double xInitial[], double mass[], 

 double a[]) { 

 readme << '\t' << ":-) *** TAPERCHAIN within walls *** (-:" 

 << '\n' << '\n'; 

 readme << "parameter of this run: " << '\n' << '\n'; 

 readme << "total number of particles: " << '\t' << '\t' << nptles << '\n'; 

 readme << "density of particles (mg/mm^3): " << '\t' << rho << '\n'; 

 readme << "quantity D of particles (mm^2/N): " << '\t' << D << '\n'; 

 readme << "radius of large ptle (mm): " << '\t' << '\t' << rlarge << '\n'; 

 readme << "tapering factor (%): " << '\t' << '\t' << '\t' << q << '\n'; 

 readme << "exponent in potential: " << '\t' << '\t' << '\t' << xn << '\n'; 

 readme << "timestepwidth (musec): " << '\t' << '\t' << '\t' << dt << '\n'; 

 readme << "# steps integration loop: " << '\t' << '\t' << nsteps << '\n'; 

 readme << "stepwidth diagnostics: " << '\t' << '\t' << '\t' << "every " 

 << idiagp << " timesteps" << '\n'; 

 readme << "stepwidth dump: " << '\t' << '\t' << '\t' << "every " 

 << idump << " timesteps" << '\n'; 

 readme << "initial v small ptle (mm/musec): " << '\t' << v1in << '\n'; 

 readme << "initial v large ptle (mm/musec): " << '\t' << vnin << '\n'; 

 readme << "restitution factor for all ptles: " << '\t' << 1-epsilon << '\n' 

 << '\n'; 

 readme << "total length of run (musec): " << '\t' << dt * nsteps << '\n'; 

 readme << "total rows recorded for .EneImp file: " << '\t' << '\t' 
41 

 << nsteps/idiagp +1 << '\n'; 

 readme << "total rows recorded for particle files: " << '\t' 

 << nsteps/idump +1 << '\n' << '\n'; 

 readme << "Initial system info (t=0): " << '\n' << '\n'; 

 readme << "kin. E. (J)" << '\t' << "pot. E. (J)" << '\t' 

 << "tot. E. (J)" << '\t' << "total impulse (mg*mm/musec)" << '\n'; 

 readme << ke << '\t' << pot << '\t' << te << '\t' << ptotal 

 << '\n' << '\n'; 

 readme << "kin. E. of small particle (J): " << '\t' << '\t' << ke1in 

 << '\n'; 

 readme << "kin. E. of large particle (J): " << '\t' << '\t' << kenin 

 << '\n'; 

 readme << "impulse of small particle (mg*mm/musec): " << '\t' 

 << p1in << '\n'; 

 readme << "impulse of large particle (mg*mm/musec): " << '\t' 

 << pnin << '\n' << '\n'; 

 readme << "particle radii (mm): " << '\n'; 

 for (int i=0; i < nptles; i++) { 

 readme << r[i] << '\t'; 

 } 

 readme << '\n' << '\n'; 

 readme << "initial particle positions (mm): " << '\n'; 

 for (int i=0; i < nptles; i++) { 

 readme << xInitial[i] << '\t'; 

 } 

 readme << '\n' << '\n'; 

 readme << "total length of one dimensional alignment (mm): " << '\t' 

 << xInitial[nptles-1] + r[nptles-1] << '\n' << '\n'; 

 readme << "particle masses (mg): " << '\n'; 

 for (int i=0; i < nptles; i++) { 

 readme << mass[i] << '\t'; 

 } 

 readme << '\n' << '\n'; 

 readme << "particle interaction strenghts (0.0316*N/mm^(3/2)): " << '\n'; 

 for (int i=0; i < nptles+1; i++) { 

 readme << a[i] << '\t'; 

 } 

 readme << '\n'; 

} 

int main ( ) { 

 double r[nptles], x[nptles], xAbsolut[nptles], xInitial[nptles], 

 v[nptles], acc[nptles], mass[nptles], a[nptles+1], 

 overbefore[nptles+1], pot; 

 /******************** functions *****************************/ 

 radii (r); 

 masses (r, mass); 

 strenghtfac (r, a); 

 initialpos (r, xInitial); 

 /********************** output files ****************************/ 

 /****** ptle-files *********/ 

 ofstream print [nptles]; 

 for (int k = 0; k < nptles; k++) { 
42 

 string filename; 

 ostringstream buffer; 

 buffer << "taperchain_" << k+1 << ".dat"; 

 filename = buffer.str(); 

 print[k].open(filename.c_str()); // convert string to char 

 } 

 for (int k = 0; k < nptles; k++) // header for particles 

 ptleHeader(& print[k], k); 

 /***************************/ 

 // header for .EneImp file 

 EnergyImpulse << "# time" << '\t' << "kin. E. (J)" << '\t' 

 << '\t' << "pot. E. (J)" << '\t' << '\t' 

 << "total E. (J)" << '\t' << '\t' 

 << "|(total imp.)| (mg*mm/musec)" << '\t' 

 << "total imp. (mg*mm/musec)" << '\n'; 

 /*****************************************************************/ 

 for (int i = 0; i < nptles; i++) { // zeroing 

 x[i] = 0.0; // relative particle positions for calculation 

 v[i] = 0.0; 

 acc[i] = 0.0; 

 overbefore[i] = 0.0; 

 } 

 overbefore[nptles] = 0.0; 

 double t = 0.0; 

 v[0] = v1in; // mind special input data 

 v[nptles-1] = vnin; 

 /**************** initial energy info *********************************/ 

 double ke1in = 0.5 * mass[0] * v[0]*v[0]; // initial kE of edge ptles 

 double kenin = 0.5 * mass[nptles-1] * v[nptles-1]*v[nptles-1]; 

 double p1in = mass[0] * v[0]; // initial impulse of edge particles 

 double pnin = mass[nptles-1] * v[nptles-1]; 

 computeAccelerations (x, a, r, acc, overbefore, mass, pot); 

 // call for initial potential energy 

 

 double ke = 0.0; // checking for initial system energy 

 double ptotal = 0.0; 

 for (int i = 0; i < nptles; i++) { 

 acc[i] = 0.0; // reset acc for verlet in case of initial overlap 

 ke += mass[i] * v[i]*v[i]; 

 ptotal += mass[i] * v[i]; 

 } 

 ke *= 0.5; 

 double te = pot + ke; // initial total system energy 

 readmeInfo (ke, pot, te, ptotal, ke1in, kenin, p1in, pnin, 

 r, xInitial, mass, a); 

 /******************* begin of timestep loop *********************/ 

 for (unsigned int i = 0; i < nsteps+1; i++) { 

 
43 

 t = i * dt; 

 velocityVerletStep (x, v, acc, a, r, overbefore, mass, pot); 

 

 /*********** check system energy, plot data *******************/ 

 if ((i % idiagp) == 0) { 

 dumpEnergyImpulse(t, ke, te, pot, ptotal, mass, v); 

 } 

 /********* check particle energy, plot data *******************/ 

 if ((i % idump) == 0) { 

 for (int k = 0; k < nptles; k++) { 

 dumpData(t, mass, v, acc, // scope dumpData 

 r, x, xInitial, xAbsolut, // scope absolutpos 

 & print[k], k); 

 } 

 } 

 }/******************** end of timestep loop *********************/ 

 readme.close(); 

 EnergyImpulse.close(); 

 

 for (int k = 0; k < nptles; k++) // closing particle files 

 print[k].close(); 

} 
44 

INTENTIONALLY LEFT BLANK. 
45 

Appendix B. Decorated Tapered Chain Modifications 

 

 This appendix appears in its original form, without editorial change. 
46 

/************************** ALTERABLE PARAMETER: ****************************/ 

int nptles = 3; // total number of particles 

const double f = 0.9; // fractional size of interstitial grain w.r.t last 

grain 

const double rho=4.42; // TiAlV (mg/mm^3) 

const double D = 0.01206; // TiAlV (mm^2/N) 

const double rlarge = 5.0; // (radius of large ptle (mm)) 

const double q = 0.0; // (tapering factor (%)) 

const double xn = 2.5; // (exponent in potential) 

const double dt = 0.00001; // (timestepwidth (musec)) 

const unsigned int nsteps = 100000000; // (# steps integration loop) 

const int idiagp = 20000; // (stepwidth diagnostics) 

const int idump = 20000; // (stepwidth dump) 

const double v1in = 0.0; // (initial v small ptle (mm/musec)) 

const double vnin = -0.01; // initial v large ptle (mm/musec)) 

const double epsilon = 1.0; // ((1 - restitution factor) all ptles) 

/****************************************************************************/ 

ofstream readme("taperchain.readme"); // global scope fcts 

ofstream EnergyImpulse("taperchain.EneImp"); 

// Generate radii and masses for DTC 

void spheres (double rlocal[], double masslocal[]) { 

 rlocal[nptles-1] = rlarge; // shifts everything to index starting at zero 

 double tapering = 1 - q*0.01; 

 const double pi = 4 * atan(1.0); 

 const double masslarge = (4.0/3.0) * pi * pow(rlarge,3) * rho; 

 masslocal[nptles-1] = masslarge; 

 if (q==0 && f==1.0) // Monodisperse in DTC 

 for (int i=0; i<nptles-1; i++) { 

 rlocal[i] = rlarge; 

 masslocal[i] = masslarge; 

 } 

 else if (q==0) { // Quasi-Monodisperse: Avoid 

roundoff errors without tapering 

 for (int i=1; i<nptles-1; i=i+2) { // Interstitial grains 

 rlocal[i] = rlarge*f; 

 masslocal[i] = (4.0/3.0) * pi * pow(rlocal[i],3) * rho; 

 } 

 for (int i=0; i<nptles-1; i=i+2) { // Non-interstitial 

grains 

 rlocal[i] = rlarge; 

 masslocal[i] = masslarge; 

 } 

 } 

 // non-monodisperse chains 

 else { 

 for (int i=1; i<nptles-1; i=i+2) { // Find radii of 

interstitial grains 

 rlocal[i] = f*pow(tapering,(nptles-1)/2)*rlarge; 

 masslocal[i] = (4.0/3.0) * pi * pow(rlocal[i],3) * rho; 

 } 

 for (int i=nptles-1; i>=0; i=i-2) { //Find radii of 

non-interstitital grains 

 rlocal[i-2] = rlocal[i] * tapering; 

 masslocal[i-2] = (4.0/3.0) * pi * pow(rlocal[i-2],3) * rho; 

 } 

 } 

} 

void strenghtfac (double r[], double alocal[]) { 

 alocal[0] = (2.0 / (5.0 * D)) * (sqrt(r[0])); 

 alocal[nptles] = (2.0 / (5.0 * D)) * (sqrt(r[nptles-1])); 
47 

 if (q == 0 && f == 1.0) // avoid roundoff errors w/out 

tapering 

 for (int i = 1; i < nptles; i++) 

 alocal[i] = (2.0 / (5.0*D)) * (sqrt(0.5*rlarge)); 

 else 

 for (int i = 1; i < nptles; i++) 

 alocal[i] = (2.0 / (5.0 * D)) * (sqrt((r[i]*r[i-1])/(r[i]+r[i-1]))); 

} 

// initialpos prints absolute initial positions, not for calculations 

void initialpos (double r[], double xInitiallocal[]) { 

 if (q == 0 && f == 1.0 ) // avoid roundoff errors w/out 

tapering 

 for (int i = 0; i < nptles; i++) 

 xInitiallocal[i] = (2.0*(i+1) - 1) * rlarge; 

 else { 

 xInitiallocal[0] = r[0]; 

 for (int i = 1; i < nptles; i++) 

 xInitiallocal[i] = xInitiallocal[i-1] + r[i-1] + r[i]; 

 } 

}
48 

INTENTIONALLY LEFT BLANK. 
49 

Appendix C. Practical Extration and Report Language (PERL) Script 

for Parametric Studies 

 

 This appendix appears in its original form, without editorial change. 
50 

#!/usr/bin/env perl 

# This program automates the considerable task of setting up parametric 

studies on the tapered 

# spherical elastic 1-D grain problem. It takes an input file: 

taperchain29.cpp and searches 

# through it replacing the values of epsilon, q, and N in for loops and 

spitting out a file 

# in the appropriate directory. The directories are created on the fly. 

# This version uses the updated directories and is looking at initial velocity 

on the small 

# grain. 

$w = 0.1; 

$FILE_NAME = "taperchain28_w01.cpp"; 

$SOURCE_DIR = "/home/rldoney"; 

#$SOURCE_DIR = "/Users/bob/Work/Classes/Spring 2004/Dr. Sen study/runs"; 

$FILE_IN = "$SOURCE_DIR/$FILE_NAME"; 

$OUT_DIR = "/nfs/scratch/rldoney/TiAlV/D.SimpleTapered"; 

#$OUT_DIR = "$SOURCE_DIR/TiAlV/D.SimpleTapered/D.Vin_large/"; 

$EPSILON_CHK = "const double epsilon ="; # Set pattern to match line with 

epsilon 

$Q_CHK = "const double q ="; # Set pattern to match 

line with q 

$N_CHK = "int nptles="; # Set pattern to match line 

with N 

system("clear"); # clear the screen; 

#for($w = 0.0; $w<=0.02; $w+=0.01) { # Restitution 

$epsilon = 1.0 - $w; # taperchain.cpp program uses epsilon instead of w 

directly 

print"\n\n:::::::: w=$w \t epsilon = $epsilon :::::::::\n"; 

system("date '+DATE: %m/%d/%y%nTIME: %H:%M:%S'"); 

print"\n"; 

chdir("$OUT_DIR") or die "Cant open $OUT_DIR"; 

mkdir("w$w"); 

chdir("w$w") or die "Cant open $OUT_DIR/w$w"; 

 

for($N=3; $N <=20; $N++) { 

 mkdir("N$N"); 

 chdir("N$N") or die "Cant open N$N"; 

 for($q=0; $q <=10; $q++) { 

 mkdir("N$N\q$q"); 

 chdir("N$N\q$q") or die "Cant open N$N\q$q"; 

 $CURRENT_DIR = `pwd`; 

 chop($CURRENT_DIR); # remove trailing \n 

 print"Current Dir: $CURRENT_DIR\n"; 

 $FILE_OUT= "$CURRENT_DIR/$FILE_NAME"; 

 open(FROM, "$FILE_IN") or die "Cant open $FILE_IN: $!"; 

 open(TO, ">$FILE_OUT") or die "Cant open $CURRENT_DIR/$FILE_OUT: 

$!"; 

 print" "; 

 system("date '+TIME: %H:%M:%S'"); 

 

 while(<FROM>) { # Read in the file line by line into $_ 

 # Replacements (Regular expression matching) 

 if($w != 0) { 

 s/$EPSILON_CHK \d+\.\d+/$EPSILON_CHK $epsilon/; # change 

epsilon 

 } 

 
51 

 s/$N_CHK\d+/$N_CHK $N/; # change N 

 s/$Q_CHK \d+\.\d+/$Q_CHK $q\.0/; # change q 

 

 print TO $_; # Write the current line with any 

changes to TO 

 close TO; 

 } 

 close FROM; 

 # Escape to the shell, compile the file, and run it 

 # It was unexpected, but we need the ' ' because of the spaces in the 

path 

 system ("g++ '$FILE_OUT'"); 

 system ("./a.out"); 

 chdir("../"); 

 } 

 chdir("../"); 

 print("\n"); 

} 

system("date '+DATE: %m/%d/%y%nTIME: %H:%M:%S'"); 
52 

INTENTIONALLY LEFT BLANK. 
53 

Appendix D. MATLAB Code for Generating Numerical Kinetic 

Energy Surfaces 

 

 This appendix appears in its original form, without editorial change. 
54 

clc; format compact; format long 

% STC directory 

cd '/Volumes/Xternal/Individual Study/TiAlV/D.SimpleTapered/D.Vin_large' 

pwd 

wi = 0; 

qmin = 0; 

qmax = 10; 

qtot = qmax-qmin + 1; 

fileprefix = 'taperchain_'; 

for w=0.0 : .02 : 0.1 

 wi = wi + 1; 

 cd (['w',num2str(w)]); 

 Ni = 0; 

 for N=3:1:20 

 Ni = Ni + 1; % Since N skips, we need an index for N 

 filesuffix = N; % For simple tapered chain 

 cd (['N',num2str(N)]) 

 foldername = (['N',num2str(N)]); % Build name of N# directory 

 foldernameq = ([foldername,'q']); % Build name of N#q 

directory 

 for q=qmin:qmax 

 %qi = q; % If we ignore the monodisperse case 

 qi = q + 1; % Indices can't be zero, so create qi as an index 

 foldername = ([foldernameq,num2str(q)]); % Foldername is determined 

by current q 

 cd([foldername]) 

 file = ([fileprefix,num2str(filesuffix)]); 

 filename = ([file,'.dat']); % Add the .dat suffix to 

current filename 

 pwd 

 % Largest Grain 

 KEinL(Ni,qi,:) = dlmread([filename],'\t','E2..E12'); % 

E(J) - KE 

 %FinL(Ni,qi,:) = dlmread([filename],'\t','F2..F1202'); % 

F(kN) - Force 

 

 % Smallest Grain 

 %t(Ni,qi,:) = dlmread('taperchain_1.dat','\t','A2..A5002'); 

 % t(us) - Time 

 KEoutS(Ni,qi,:) = dlmread('taperchain_1.dat','\t','E2..E1502'); 

 % E(J) - KE 

 %FoutS(Ni,qi,:) = dlmread('taperchain_1.dat','\t','F2..F1202'); 

 % F(kN) - Force 

 

 % INPUTS 

 % ======= 

 % Maximum input kinetic energy is 1st element (initial velocity) 

 % in the largest grain of every chain 

 KEmax_in_L(Ni,qi,:) = KEinL(Ni,qi,1); 

 

 % Maximum input Force is somewhere early on for 1st element. 

 %% ANDed portion is to make sure a small peak doesn't occur earlier 

 %for i=1:1200 

 % if (FinL(Ni,qi,i+1) < FinL(Ni,qi,i)) && (FinL(Ni,qi,i) > 0.001) 

 % Fmax_in_L(Ni,qi) = FinL(Ni,qi,i); 

 % break 

 % end 

 %end 

 

 % OUTPUTS 

 % ======== 
55 

 % New technique of normalizing: wrt first peak hitting last grain. 

 % peak occurs when comparing the difference between the ith and 

 % (i+1)st elements. when it is less than zero we have just passed a 

 % peak 

 %% ANDed portion is to make sure a small peak doesn't occur earlier 

 for i=1:1200 % Only need a small part of t 

 

 if (KEoutS(Ni,qi,i+1) < KEoutS(Ni,qi,i)) && (KEoutS(Ni,qi,i) > 

0.001) 

 KEmaxF_out_S(Ni,qi) = KEoutS(Ni,qi,i); 

 break 

 end 

 

 % For force, since the acceleration is the other way, the sign 

 % < goes to >. And we want only want the magnitude of the force 

 %FmaxF_out_S(Ni,qi) = abs(min(FoutS(Ni,qi,:))); 

 %if (FoutS(Ni,qi,i+1) > FoutS(Ni,qi,i)) 

 % FmaxF_out_S(Ni,qi) = abs(FoutS(Ni,qi,i)); 

 % break 

 %end 

 i = i+1; 

 end 

 cd .. 

 clear foldername 

 end 

 cd .. 

 end 

 %Fnorm_in_L = FmaxF_out_S ./ Fmax_in_L; % Normalize Fout/Fin per 

each specific grain 

 KEnorm_in_L = KEmaxF_out_S ./ KEmax_in_L; % Normalize Fout/Fin per 

each specific grain 

 cd '/Volumes/Xternal/Individual Study/TiAlV/D.SimpleTapered/D.Vin_large' 

 

 Ni = linspace(3,20,18); 

 qi2 = linspace(qmin,qmax*0.01,qtot); 

 % Figure labeling 

 if wi == 1 

 letter = 'a'; 

 elseif wi == 2 

 letter = 'b'; 

 elseif wi == 3 

 letter = 'c'; 

 elseif wi == 4 

 letter = 'd'; 

 elseif wi == 5 

 letter = 'e'; 

 else 

 letter = 'f'; 

 end 

 % For some reason, these 3-d plots switch the x and y plotting. The 

 % variables are correct, but the ordering was unexpected 

 %plot( squeeze( Ft(1,2,:)) , squeeze( Fin(1,2,:) ) ) % N = 3 , q = 

1 

 subplot(3,2,wi) 

 %surf(qi2,Ni,Fnorm_in_L) 

 surf(qi2,Ni,KEnorm_in_L) 

 xlabel('\bf q') 

 ylabel('\bf N_i') 

 %if (wi == 1 | wi == 3 | wi == 5) 

 zlabel('\bf \it KE_N') 

 %end 

 text(0,18,0.4,['\it \omega = ',num2str(w)]) 
56 

 % subplot labeling based on index of wi 

 text(0,18,0.5,['\bf (',letter,')'],'FontSize',14,'FontName','Times') 

 axis([0.0 0.1 3 21 0 0.5]) 

 view(135,10) 

 caxis([0 0.5]) % Adjust the limits of the color scheme for Z, must 

come before colorbar 

 %colorbar 

end 

cd '/Applications/Physics/Matlab7' 
57 

Appendix E. Normalization 

Normalization has posed some challenges in trying to properly assess the absorption quality of a 

simple tapered chain (STC). Proper normalization schemes will help the architect better 

determine which STC is best for which application. Recall that the goal is to measure the energy 

at the last grain versus the energy put into the system by the first. In general, the functional form 

will stay the same regardless of the strategy and adjustments in the normalization will simply 

scale the kinetic energy (KE) surface. In some cases, the output force of each chain is based on 

the maximum value felt by any possible chain under consideration (i.e., monodisperse and no 

energy loss).1

 This gives a measure of how one chain is better than another. 

In this communication, we have chosen to form the ratio based on the output KE and force felt 

by each specific chain. This serves to grade the individual effectiveness of any chain without 

reference to another. In choosing a peak value with this method, one could identify either when 

the impulse first hits the last grain or look for the absolute maximum peak whenever it may 

occur. We have chosen to use the former for several reasons. First, it ignores the complexity of 

nonlinear reverberations which can lead to large peaks at unpredictable times. Second, we argue 

that this is just as realistic as selecting the maximum value anywhere in the time spectrum. 

It turns out that in most cases, the absolute maximum is the first peak. There are special cases 

where the maximum may occur at later times and this needs to be investigated further. In 

figure 6 (in the body of the report) for q = 0.1, for example, we see the striking occurrence of the 

secondary pulse about 225 μs being much stronger than the initial arrival at 35 μs. This is one of 

those instances that disagrees with the way we choose to normalize our KE surfaces. We have 

investigated this particular case further without including extraneous plots and report the 

following observations. The effect exists for N =15 − 20 for constant ω. When N = 20 is held 

fixed and restitution is increased, the peak KE once again occurs for the first arrival of the pulse. 

As q increases, so do the number of collisions and the requisite energy loss (sinceω ≠ 0). 

Therefore, it is less likely to find a global peak later in the simulation. The situation is further 

complicated by the interplay between q and ω.

 1

 Pfannes, J. Energy Propagation in Granular Chains. M.S. Thesis, State University of New York, Buffalo, NY, May 2003. 
NO. OF 

COPIES ORGANIZATION 

58 

 1 DEFENSE TECHNICAL 

 (PDF INFORMATION CTR 

 ONLY) DTIC OCA 

 8725 JOHN J KINGMAN RD 

 STE 0944 

 FORT BELVOIR VA 22060-6218 

 1 US ARMY RSRCH DEV & 

 ENGRG CMD 

 SYSTEMS OF SYSTEMS 

 INTEGRATION 

 AMSRD SS T 

 6000 6TH ST STE 100 

 FORT BELVOIR VA 22060-5608 

 1 INST FOR ADVNCD TCHNLGY 

 THE UNIV OF TEXAS 

 AT AUSTIN 

 3925 W BRAKER LN STE 400 

 AUSTIN TX 78759-5316 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 IMNE ALC IMS 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 DIRECTOR 

 US ARMY RESEARCH LAB 

 AMSRD ARL CI OK TL 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 DIRECTOR 

 US ARMY RESEARCH LAB 

 AMSRD ARL CS IS T 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

ABERDEEN PROVING GROUND

 1 DIR USARL 

 AMSRD ARL CI OK TP (BLDG 4600) 
NO. OF NO. OF 

COPIES ORGANIZATION COPIES ORGANIZATION

59 

 5 INST FOR DEFENSE ANALYSES 

 H BERTRAND 

 G BOEZER 

 J HEAGY 

 I KOHLBERG 

 M RIGDON 

 4850 MARK CENTER DR 

 ALEXANDRIA VA 22311-1882 

 7 STATE UNIVERSITY OF NEW YORK 

 DEPT OF PHYSICS 

 R GONSALVES 

 B POWELL 

 S SEN (5 CPS) 

 239 FRONCZIAK HALL 

 BUFFALO NY 14260-1500 

 1 STATE UNIVERSITY OF NEW YORK 

 DEPT OF MATHEMATICS 

 B PITMAN 

 244 MATHEMATICS BLDG 

 BUFFALO NY 14260-2900 

ABERDEEN PROVING GROUND

 37 DIR USARL 

 AMSRD ARL WM 

 T WRIGHT 

 AMSRD ARL WM MD 

 B CHEESEMAN 

 AMSRD ARL WM T 

 B BURNS 

 M ZOLTOSKI 

 AMSRD ARL WM TA 

 A BARD 

 P BARTKOWSKI 

 R DONEY (10 CPS) 

 G FILBEY 

 B FREY 

 V HERNANDEZ 

 E HORWATH 

 T JONES 

 M KEELE 

 P KINGMAN 

 D KLEPONIS 

 C KRAUTHAUSER 

 B LEAVY 

 I LISKO 

 A MIHALCIN 

 J RUNYEON 

 S SCHOENFELD 

 AMSRD ARL WM TB 

 R SKAGGS 

 

 AMSRD ARL WM TC 

 M FERMEN-COKER 

 AMSRD ARL WM TE 

 P BERNING 

 C HUMMER 

 J POWELL 

 AMSRD ARL SL BE 

 B BRUCHEY 

 T BJERKE 
60

INTENTIONALLY LEFT BLANK.