Formation of Gaps at the Specimen-Bar Interfaces in 

Numerical Simulations of Compression Hopkinson Bar 

Tests on Soft, Nearly Incompressible Materials

by Mike Scheidler and Martin N. Raftenberg

ARL-TR-5301 September 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-5301 September 2010

Formation of Gaps at the Specimen-Bar Interfaces in

Numerical Simulations of Compression Hopkinson Bar

Tests on Soft, Nearly Incompressible Materials

Mike Scheidler and Martin N. Raftenberg

Weapons and Materials Research Directorate, ARL

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 2010

2. REPORT TYPE

Final

3. DATES COVERED (From - To)

July 2008–December 2009

4. TITLE AND SUBTITLE

Formation of Gaps at the Specimen-Bar Interfaces in Numerical Simulations of 

Compression Hopkinson Bar Tests on Soft, Nearly Incompressible Materials

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Mike Scheidler and Martin N. Raftenberg

5d. PROJECT NUMBER

105

5e. TASK NUMBER

5f. WORK UNIT NUMBER

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

U.S. Army Research Laboratory

ATTN: RDRL-WMP-B

Aberdeen Proving Ground, MD 21005-5066

8. PERFORMING ORGANIZATION

 REPORT NUMBER

ARL-TR-5301

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

14. ABSTRACT

Numerical simulations of compression Hopkinson bar tests on soft, nearly incompressible materials were performed in an effort

to understand inertial effects in these tests. The specimen was modeled as either a linear elastic or a nonlinear elastic solid; the 

nonlinear model was calibrated to ballistic gelatin data. The simulations revealed a previously unreported phenomenon, 

namely, for sufficiently high strain rates in the specimen and sufficiently short rise times to the final strain rate, small gaps 

formed at both the specimen-incident bar and the specimen-transmission bar interfaces. The size of these gaps ranged from 

sub-micron to as large as 43 m. Gaps formed at small specimen strains but persisted out to large strains, closing and reopening multiple times. At some instants the gaps extended over most of the face of the specimen. Unloading of the specimen 

and the pressure bars accompanied these gaps. In an effort to rule out numerical artifacts, the sensitivity of these results to 

mesh size and contact algorithm parameters was studied.

15. SUBJECT TERMS

Hopkinson bar, gaps, nearly incompressible, Mooney-Rivlin, inertial effects

16. SECURITY CLASSIFICATION OF: 

17. LIMITATION

OF ABSTRACT

UU

18. NUMBER

OF PAGES

122

19a. NAME OF RESPONSIBLE PERSON

Mike Scheidler

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

410-278-5436

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
iii

Contents

List of Figures v

List of Tables viii

Acknowledgments ix

1. Introduction 1

2. Problem Description 3

2.1 Problem Geometry ..........................................................................................................3

2.2 Loading and Boundary Conditions..................................................................................5

2.3 Material Properties..........................................................................................................6

2.4 Wave Speeds and Travel Times......................................................................................8

2.5 Measurement of Gap Size .............................................................................................10

2.6 Stress, Stretch and Strain...............................................................................................11

3. The Mooney-Rivlin Constitutive Model for the Specimen 14

3.1 Incompressible and Compressible Versions..................................................................14

3.2 Calibration of the Mooney-Rivlin Model for Ballistic Gelatin.....................................15

4. Computational Parameters 16

4.1 Geometry and Meshing .................................................................................................16

4.2 Boundary Conditions.....................................................................................................17

4.3 Miscellaneous................................................................................................................18

5. Wave Propagation in the Pressure Bars 18

5.1 General Remarks...........................................................................................................18

5.2 Considerations Governing the Bar Lengths ..................................................................19

5.3 Velocity Histories in the Incident Bar...........................................................................20

6. Simulations with the Moony-Rivlin Model and a 1-s Rise Time 25

6.1 Results on the Centerline...............................................................................................25
iv

6.1.1 General Remarks...............................................................................................26

6.1.2 The Initial Loading of the Specimen: Stage I (140–165 s)............................32

6.1.3 Gap Opening and Closure (Stages II–VI) .........................................................34

6.2 Results at Other Radial Locations.................................................................................35

6.3 Pressure Contours..........................................................................................................38

6.4 Effects of Pressure Bar Radius on Gap Formation........................................................47

7. Simulations with the Moony-Rivlin Model and a 25-s Rise Time 49

7.1 Results on the Centerline...............................................................................................49

7.1.1 General Remarks...............................................................................................49

7.1.2 The Initial Loading of the Specimen: Stage I (143–184 s)............................54

7.1.3 Gap Opening and Closure (Stages II–IV) .........................................................55

7.2 Gap Size vs. Radius at Selected Times .........................................................................56

7.3 Force History in the Transmission Bar..........................................................................57

8. Simulations with a Linear Elastic Model for the Specimen 61

8.1 The 1-s Initial Rise Time ............................................................................................61

8.2 The 25-s Initial Rise Time ..........................................................................................67

9. Discussion and Concluding Remarks 72

10. References 75

Appendix A. The Mooney-Rivlin Constitutive Model 79

Appendix B. Sensitivity Studies 85

Appendix C. Additional Computational Studies 99

List of Symbols, Abbreviations, and Acronyms 104

Distribution List 106
v

List of Figures

Figure 1. Geometry and coordinate system for the SHPB simulations (not drawn to scale). ........4

Figure 2. Meshes for the Incident Bar (IB), Transmission Bar (TB), and Specimen (S). The 

baseline mesh for the specimen is shown in case (a)...............................................................17

Figure 3. The scaled axial velocity vs. non-dimensional time at three locations on the 

centerline of the incident bar: the loading end with a 1-s rise time, 1/3 of the way from 

the loading end, and 2/3 of the way from the loading end (blue curve). .................................21

Figure 4. The scaled axial velocity vs. non-dimensional time at three locations on the 

centerline of the incident bar: the loading end with a 25-s rise time, 1/3 of the way from 

the loading end, and 2/3 of the way from the loading end (blue curve). .................................22

Figure 5. The scaled axial velocity history on the centerline of the incident bar at the S-IB 

interface for a 1-s (initial) rise time......................................................................................23

Figure 6. The scaled axial velocity history on the centerline of the incident bar at the S-IB 

interface for a 25-s (initial) rise time.....................................................................................24

Figure 7. The histories of the axial displacement of the incident bar and the specimen on the 

centerline at the S-IB interface. Also shown is the history of the corresponding gap size 

(blue)........................................................................................................................................26

Figure 8. The histories of the axial displacement of the transmission bar and the specimen on 

the centerline at the S-TB interface. Also shown is the history of the corresponding gap 

size (blue).................................................................................................................................27

Figure 9. A comparison of the histories of the gap sizes on the centerline at the S-IB and 

S-TB interfaces. Also shown is the history of mean axial stretch in the specimen (blue)......28

Figure 10. A comparison of the histories of the gap sizes on the centerline at the S-IB and 

S-TB interfaces in the vicinity of the time at which they first form. Also shown is the 

history of mean axial stretch in the specimen (blue). ..............................................................29

Figure 11. Histories of the axial stress (negative in compression) and the gap size on the 

centerline at the S-IB interface. The stress is measured at the centroids of the specimen 

and incident bar elements adjacent to the centerline. ..............................................................31

Figure 12. Histories of the axial stress (negative in compression) and the gap size on the 

centerline at the S-TB interface. The stress is measured at the centroids of the specimen 

and transmission bar elements adjacent to the centerline. .......................................................32

Figure 13. Deformed mesh plots in the vicinity of the specimen at t  236 s and an axial 

strain of 21%: (a) full specimen mesh; (b) enlargement of the mesh near the centerline, 

where the gap size is about 37 m...........................................................................................36

Figure 14. Gap size vs. deformed specimen radius at the S-IB interface at selected times..........37

Figure 15. Gap size vs. deformed specimen radius at the S-IB interface at selected times..........37

Figure 16. Gap size vs. deformed specimen radius at the S-TB interface at selected times.........38

Figure 17. Pressure contours in the vicinity of the specimen at selected times............................40
vi

Figure 18. Pressure contours in the vicinity of the specimen at selected times............................41

Figure 19. Pressure contours in the vicinity of the specimen at selected times............................42

Figure 20. Pressure contours in the vicinity of the specimen at selected times............................43

Figure 21. Pressure contours in the vicinity of the specimen at selected times............................45

Figure 22. Pressure contours in the vicinity of the specimen at selected times............................46

Figure 23. Pressure contours in the vicinity of the specimen at selected times............................47

Figure 24: Histories of the gap size on the centerline at the S-IB interface (top) and S-TB 

interface (bottom) for two different bar radii...........................................................................48

Figure 25. Histories of the gap size on the centerline at the S-IB interface for two loading 

waves with different initial rise times......................................................................................50

Figure 26. Histories of the gap size on the centerline at the S-TB interface for two loading 

waves with different initial rise times......................................................................................51

Figure 27. Histories of the axial stress (negative in compression) and the gap size on the 

centerline at the S-IB interface. The stress is measured at the centroids of the specimen 

and incident bar elements adjacent to the centerline. ..............................................................52

Figure 28. Histories of the axial stress (negative in compression) and the gap size on the 

centerline at the S-TB interface. The stress is measured at the centroids of the specimen 

and incident bar elements adjacent to the centerline. ..............................................................53

Figure 29. Gap size vs. deformed specimen radius at the S-IB interface at selected times..........57

Figure 30. History of the total force (positive in compression) at the stress gage location in 

the transmission bar (2.52 mm from the specimen interface) compared with the history of 

the gap size on the centerline at the S-TB interface.................................................................58

Figure 31. A comparison of the histories of the gap sizes on the centerline at the S-IB and

S-TB interfaces for a 1-s rise time and   0.49999 (linear elastic model). Also shown 

is the history of mean axial stretch in the specimen (blue)......................................................62

Figure 32. A comparison of the histories of the gap sizes on the centerline at the S-IB 

interface for a 1-s rise time and selected values of Poisson’s ratio (linear elastic model). ...63

Figure 33. A comparison of the histories of the gap sizes on the centerline at the S-TB 

interface for a 1-s rise time and selected values of Poisson’s ratio (linear elastic model). ...64

Figure 34. The maximum gap size at the centerline for each interface as a function of 

Poisson’s ratio (linear elastic model, 1-s rise time)...............................................................65

Figure 35. Deformed mesh in the vicinity of the specimen at t  276 s for Poisson’s ratio 

 (linear elastic model, 1-s rise time). ....................................................................66

Figure 36. A comparison of the histories of the gap sizes on the centerline at the S-IB and 

S-TB interfaces for a 25-s rise time and   0.49999 (linear elastic model). Also shown 

is the history of mean axial stretch in the specimen (blue)......................................................67

Figure 37. A comparison of the histories of the gap sizes on the centerline at the S-IB 

interface for a 25-s rise time and selected values of Poisson’s ratio (linear elastic 

model). .....................................................................................................................................68
vii

Figure 38. A comparison of the histories of the gap sizes on the centerline at the S-TB 

interface for a 25-s rise time and selected values of Poisson’s ratio (linear elastic 

model). .....................................................................................................................................69

Figure 39. The maximum gap size at the centerline for each interface as a function of 

Poisson’s ratio (linear elastic model, 25-s rise time).............................................................70

Figure 40. Deformed mesh in the vicinity of the specimen at t  286 s for Poisson’s ratio 

 (linear elastic model, 25-s rise time). ..................................................................71

Figure A-1. Stress components as a function of axial stretch z for a uniaxial stress test on a 

single element using the compressible Mooney-Rivlin model with   0.495 (colored 

curves). Comparison with axial stress zz for the incompressible model (dashed curve).......84

Figure B-1. Histories of the gap sizes on the centerline at the S-IB interface for three 

specimen meshes and a 1-s initial rise time for the loading wave.........................................86

Figure B-2. Histories of the gap sizes on the centerline at the S-TB interface for three 

specimen meshes and a 1-s initial rise time for the loading wave.........................................87

Figure B-3. Histories of the gap sizes on the centerline at the S-IB interface for three 

specimen meshes and a 25-s initial rise time for the loading wave.......................................88

Figure B-4. Histories of the gap sizes on the centerline at the S-TB interface for three 

specimen meshes and a 25-s initial rise time for the loading wave.......................................89

Figure B-5. Histories of the gap sizes on the centerline at the S-IB interface for three values 

of the contact parameter SFACT and a 1-s initial rise time for the loading wave. ...............91

Figure B-6. Histories of the gap sizes on the centerline at the S-TB interface for three values 

of the contact parameter SFACT and a 1-s initial rise time for the loading wave. ...............92

Figure B-7. Histories of the gap sizes on the centerline at the S-IB interface for two values of 

the contact parameter VDC and a 1-s initial rise time for the loading wave.........................93

Figure B-8. Histories of the gap sizes on the centerline at the S-TB interface for two values 

of the contact parameter VDC and a 1-s initial rise time for the loading wave. ...................94

Figure B-9. Histories of the gap sizes on the centerline at the S-IB interface for three values 

of the contact parameter SFACT and a 25-s initial rise time for the loading wave. .............95

Figure B-10. Histories of the gap sizes on the centerline at the S-TB interface for three 

values of the contact parameter SFACT and a 25-s initial rise time for the loading wave. ..96

Figure B-11. Histories of the gap sizes on the centerline at the S-IB interface for two values 

of the contact parameter VDC and a 25-s initial rise time for the loading wave. .................97

Figure B-12. Histories of the gap sizes on the centerline at the S-TB interface for two values 

of the contact parameter VDC and a 25-s initial rise time for the loading wave. .................98

Figure C-1. Pressure contours in the vicinity of an annular specimen at 7% axial strain. Any 

positive pressure (compression) is shaded red; other colors indicate negative (tensile) 

pressure in units of GPa. The loading wave had a 1-s initial rise time. Observe the gap 

along most of the S-IB interface. ...........................................................................................102
viii

List of Tables

Table 1. Linear elastic constants for materials used in the numerical simulations.........................7

Table 2. Geometric parameters and wave travel times. ..................................................................9

Table 3. Parameters for the compressible Mooney-Rivlin Model................................................15
ix

Acknowledgments

The authors express their appreciation to Bryan Love for performing independent numerical 

simulations to verify some of the results presented here.
x

INTENTIONALLY LEFT BLANK.
1

1. Introduction

The split Hopkinson pressure bar (SHPB), also known as the Kolsky bar, is an experimental 

technique that is widely used for studying the strain rate sensitivity of inelastic materials in a 

state of compressive uniaxial stress. The standard relations for analyzing data from SHPB tests 

rely on the assumption of uniform stress, strain and strain rate throughout the specimen. These 

conditions, as well as a nearly constant nominal strain rate, can often be achieved after an initial 

“ringing-up” period. The review articles by Gray (1) and Gama et al. (2) are recommended for 

historical background, discussions of experimental procedures and data analysis techniques, and 

additional references to the literature. 

The most common version of the SHPB test assumes that the two pressure bars undergo small 

strain, linear elastic deformations, that the specimen remains in contact with the pressure bars, 

and that the faces of the bars in contact with the specimen remain planar. These conditions 

require that the specimen be softer than the pressure bars. However, difficulties arise in the 

analysis of the experimental data when the specimen is extremely soft relative to the bars. These 

difficulties as well as techniques that have been developed to overcome them are discussed in 

Gray (1) and Gama et al. (2) and also in Gray and Blumenthal (3), Song and Chen (4), Moy et al. 

(5), and Song et al. (6). One source of difficulty is that the signal in the transmission bar may be 

too weak to provide accurate measurements of stress in the specimen; this issue is not addressed 

here. Another difficulty when the imposed axial strain rate is sufficiently high is that the stress 

and strain within a soft specimen may be non-uniform and the stress state non-uniaxial 

throughout the test. Such conditions either invalidate the test or at least require appropriate 

“inertial” corrections to the test data. 

The original objective of the present work was to perform numerical simulations on very soft 

materials in order to test the validity of some approximate, analytical, inertial corrections for 

SHPB test data. These relations, which were developed by Scheidler (7), utilized the simplifying 

approximation of constant volume and thus are limited to nearly incompressible materials, that 

is, to materials that are substantially stiffer in dilatation than in shear. Examples of soft, nearly 

incompressible materials include many types of rubber as well as many biological materials due 

to their high water content. Ballistic gelatin, which is 80–90% water by mass, also falls into this 

category; it is widely used as a tissue surrogate in impact and penetration tests.

We performed numerous axisymmetric finite-element simulations of SHPB tests on soft, nearly 

incompressible materials using the commercial code LS-DYNA (8). The initial simulations used 

a nonlinear elastic model for the specimen, namely, LS-DYNA’s compressible version of the 

Mooney-Rivlin model (8, 9). The model was calibrated to give rough agreement with the 

uniaxial compression data on ballistic gelatin obtained by Moy et al. (5) and with the pressurevolume data on ballistic gelatin in Aihaiti and Hemley (10). We used a relatively thin specimen 
2

(1.45-mm thickness, 12.7-mm diameter) and a nominal strain rate of 2500/s in all simulations. 

These were chosen to agree with the test conditions in Moy et al. (5). For these initial 

simulations we used solid (as opposed to annular) specimens.

Much to our surprise, we observed that small gaps formed along both the specimen-incident bar 

(S-IB) interface and the specimen-transmission bar (S-TB) interface. The size of these gaps 

ranged from sub-micron to as large as 43 m, depending on the rise time of the loading wave in 

the incident bar, the time after the arrival of the wave at the specimen, the particular interface, 

and the radial location on that interface. The gaps closed and re-opened multiple times, and did 

not necessarily form over the entire face of the specimen. However, for the S-IB interface and 

the shortest rise time, a gap existed over most of the face of the specimen for axial strains 

ranging from 4% to 19%, and out to half of the specimen radius for strains up to 28%.1 The 

faces of the pressure bars remained very nearly planar, as expected; the gaps formed as the 

specimen moved axially inward toward its center. Large pressure spikes were observed on the 

axis of the bars and the specimen as these gaps closed up and the specimen “slapped” the bars. 

To the best of our knowledge, this phenomenon has not been reported in either the experimental 

or the computational literature on SHPB tests,2 although certain features of previously reported 

experimental data appear to be consistent with the opening and closing of gaps.

The assumption of contact between the specimen and the bars is implicit in any analysis of the 

SHPB test, including the inertial corrections that we had intended to validate. Consequently, the 

focus of our study shifted to the gap phenomenon. In particular, we sought to determine whether 

the gaps were a numerical artifact. The dependence of gap formation on the following 

conditions was examined:

1. Constitutive Model for the Specimen:

a) Model type:

i) nonlinear elastic (compressible Mooney-Rivlin)

ii) linear elastic

b) Model parameters (linear elastic case only)

2. Loading Condition:

a) Linear ramp (1 or 25 s) to constant velocity at the far end of the incident bar

b) Impact of a striker bar on the far end of a long incident bar

c) Direct impact of a striker bar on the specimen

 

1 Gaps at the S-IB interface are clearly visible in figure 13 (solid specimen) and figure C-1 (annular specimen).

2 However, a summary of some of the results in this report is given the proceedings paper (11).
3

3. Geometry:

a) Pressure bar diameter 

b) Specimen geometry:

i) solid

ii) annular

4. Boundary Condition on the Lateral Surfaces:

a) Stress-free 

b) Radially constrained

5. Computational Parameters:

a) Specimen mesh size

b) Contact algorithm parameters

6. The code used for the simulations

The organization of the report can be inferred from the table of contents. Sections 2, 3, 4, 5.1, 

5.2, and appendix A contain preliminary and background material. The results of a selected set 

of numerical simulations are presented in detail in sections 5.3, 6, 7, 8, and appendix B. 

Appendix C consists of less detailed summaries of simulations for which some of the conditions 

differed from those in the main body of the report. These conditions are italicized in the list 

above. In particular, all of the numerical simulations discussed in the main body of this report 

were performed with LS-DYNA. In appendix C-3, we summarize the results of some analogous 

simulations performed by Bryan Love at the U.S. Army Research Laboratory (ARL) using the 

finite element code PRESTO (from Sandia National Laboratories) with loading condition 

described in 2.b above.

2. Problem Description

2.1 Problem Geometry

The SHPB experiment was modeled using the three parts shown in figure 1, which is not to 

scale. As indicated in the figure, a cylindrical coordinate system was used with the origin 

located on the axis of symmetry or centerline of the specimen and the pressure bars at the S-IB

interface. The reference or material coordinates (R, Z ) are the coordinates in the undeformed 

state at time t  0; the corresponding coordinates in the deformed state are denoted by (r, z). For 

axisymmetric deformations, the Z-axis, the z-axis, the set of points with r  0, and the set of 

points with R  0 all describe the centerline. The positive direction of the Z-axis is along the 

direction of propagation of the initial loading wave, that is, from incident bar to specimen to 
4

transmission bar. All coordinate values and problem dimensions are given in millimeters, 

whereas element sizes (section 4.1) and gap sizes are specified in microns.

Figure 1. Geometry and coordinate system for the SHPB simulations (not drawn to scale).

A solid, disk-shaped specimen (S) was used for most of the simulations described in this report. 

A brief summary of simulations involving annular, washer-shaped specimens is given in 

appendices C-2 and C-3. Unless specified otherwise, a solid specimen is assumed in subsequent 

discussions. The specimen has initial thickness (or length) LS  1.45 mm and an initial diameter 

of 12.7 mm, giving a length-to-diameter ratio of 0.114. Relatively thin specimens are commonly 

used to decrease the time it takes to “ring-up” to a uniform state.

The incident bar (IB) and the transmission bar (TB) have the same diameter. Two bar diameters 

were considered, 25.6 mm and 19.0 mm. Most of the simulations were done with the larger bar 

diameter. Discussion of simulations with the smaller bar diameter is confined to section 6.4 and 

appendix C-2. The length of the incident bar used in the simulations is LIB  768 mm. This 

gives a length-to-diameter ratio of 30.0 for the larger diameter bar and 40.4 for the smaller 

diameter bar. The length of the transmission bar used in the simulations is LTB  256 mm, which 

is 1/3 the length of the incident bar and gives length-to-diameter ratios of 10.0 and 13.5. 

The dimensions of the specimen and the smaller diameter bar were chosen to agree with those 

used in the experimental study on gelatin by Moy et al. (5). However, we used the larger 

diameter bar in most of our simulations. This was motivated by our original goal of studying 

inertial effects at large strains. We wished to obtain nominal axial strains of up to 75% in 

compression, which corresponds to an axial stretch of 0.25, that is, a specimen compressed to 1⁄4 

of its original thickness. This results in a radial stretch of 2.0, assuming incompressibility and a 

homogeneous deformation. Hence the diameter of the specimen would increase by a factor of 2, 

that is, from 12.7 to 25.4 mm. And because the deformed specimen must not extend beyond the 

pressure bars, we chose a diameter of 25.6 mm for the bars. 

RS

LTB

LS

Specimen

(S) R

RB RB

LIB

non-reflecting 

boundary

prescribed 

velocity

Z

Incident Bar

(IB)

Transmission Bar

(TB)
5

We will often refer to the specimen and bar radii rather than the diameters. The initial radius of 

the specimen is RS  6.35 mm, and the two pressure bar radii considered are RB  12.8 mm and 

RB  9.5 mm.

2.2 Loading and Boundary Conditions

In the simplest version of the SHPB test, a striker bar impacts the incident bar and generates a 

compressive stress pulse which propagates down the length of the incident bar and subsequently 

into the specimen and the transmission bar. A direct impact of a striker bar on the incident bar 

would have been easy to simulate, but it can also be approximated reasonably well by a 

prescribed step in the axial velocity at the end of the incident bar.

For the purpose of reducing inertial effects in soft specimens, it is beneficial to insert a pulse 

shaper between the striker and incident bars (4–6). This is a thin disk of another soft material 

(generally not the same as the specimen) which serves to smooth the incident pulse and increase 

its rise time. To avoid the complexity of modeling the pulse shaping material and, in particular, 

the computational issues associated with the extremely large deformations of the pulse shaper, 

we chose instead to approximate the effect the pulse shaper on the incident pulse. As in the 

computational study by Song et al. (6), we prescribed an axial velocity, vz (t), on the impact end 

of the incident bar. This prescribed velocity has an initial value of zero and rises linearly with 

time until attaining a plateau value v0 at time tR:

0

0

v ; 0

v ( )

v ;

R

R

R

z

t

t t

t t

t t





  

 



 

(1)

We set v0 to 1.8125 m/s in all calculations. This value was chosen to produce a nominal axial 

strain rate of approximately 2500/s in the specimen, in agreement with the tests in Moy et al. (5). 

We considered two values for the initial rise time tR, 1 and 25 s. 

In an SHPB test, the lateral surfaces of the pressure bars (R  RB) and the specimen (R  RS) are 

unconstrained and hence stress-free. The same condition was used in the numerical simulations 

with two exceptions. For these two exceptional cases (discussed in appendix C-1.1), the lateral 

surfaces were constrained in order to study the effects of eliminating the unloading waves from 

these surfaces and reducing the radial acceleration in the specimen.

Friction between the faces of the specimen and the bars can result in bulging of the specimen on 

compression. In SHPB tests the specimen faces are lubricated to eliminate (or at least 

substantially reduce) the friction. Consequently, in the numerical simulations the S-IB and S-TB 

interfaces were modeled with frictionless contact (see section 4.2).

At the far end of the transmission bar, a non-reflecting boundary condition was used; refer to 

sections 4.2 and 5.2 for further discussion.
6

2.3 Material Properties

For soft specimens, pressure bar materials with low mechanical impedance are often used to 

improve the signal to noise ratio (1–6). For example, the tests reported in Moy et al. (5) and 

Song et al. (6) utilized 7075-T6 aluminum bars. For the imposed velocity history discussed 

above, the stress in an aluminum bar remains well below yield (see section 5.1). We used an 

isotropic, linear elastic constitutive model for the incident and transmission bars, with material 

properties appropriate for aluminum. These properties (obtained from Song et al. [6]) are a 

density of 2.7 g/cm3, a Young’s modulus of 68 GPa, and a Poisson’s ratio of 0.33.

As mentioned in the Introduction, we chose ballistic gelatin as a representative soft, nearly 

incompressible material, and for this we used a compressible version of the classical MooneyRivlin constitutive model for incompressible materials. This isotropic, nonlinear elastic 

constitutive model and its calibration are discussed briefly in section 3 and in more detail in 

appendix A, and the results of simulations with this model are given in sections 6 and 7 and 

appendices B and C. At this point we simply note that the density of the specimen was chosen to 

be that of water, 1 g/cm3. This same density was also used in subsequent simulations with a 

linear elastic model for the specimen (see section 8).3 These simulations were performed to 

determine whether or not the formation of gaps was attributable to the nonlinearity and/or near 

incompressibility of the constitutive model. It turned out that gaps formed for all but one of the 

six sets of material parameters used in the linear elastic model for the specimen.

Table 1 lists the density and linear elastic properties used for the pressure bars and the specimen 

in the various simulations. For the ballistic gelatin specimen, these parameters represent the 

small strain, linear elastic approximation to the compressible Mooney-Rivlin model (section 3) 

and so do not accurately reflect the large strain response. Nevertheless, the table allows one to 

quickly assess the relative stiffnesses of the materials. E, G, K, and L denote the Young’s 

modulus, shear modulus, bulk modulus, and longitudinal modulus; they govern the response to 

uniaxial stress, simple shear, pure dilatation, and uniaxial strain, respectively.  and  denote 

Poisson’s ratio and the initial density. Since all the materials are isotropic, any two of the 

material parameters E, G, K, L, and  determine the other three. For example,

 

2 1

3 1 2

K G 









 and 

1 4

(1 )(1 2 ) 3

L E K G 

 



  

 

. (2)

 

3 However, the linear elastic constants used in these simulations are not necessarily those appropriate for the small strain 

response of ballistic gelatin.
7

Table 1. Linear elastic constants for materials used in the numerical simulations.

Material

Density

(g/cm3)







Elastic Moduli

(GPa)

Dimensionless Elastic 

Constants

Wave Speeds

(mm/s)

Impedance Ratios

(specimen/bar)

E G K L  G/K cL cE cG

Based 

on cL

Based 

on cE

Aluminum 

Bars 2.70 68.0 25.6 66.7 101. 0.33 0.383 6.11 5.02 3.08

 

Ballistic 

Gelatin 

Specimen 

1.00

4

2.40 10

5

8.00 10

4.00 4.00 0.49999

5

2.00 10

2.00 0.0155 0.00894 0.121 0.00114

Linear 

Elastic 

Specimen

1.00

4

4.30 10

4

1.43 10

7.17 7.17 0.49999

5

2.00 10

2.68

0.0207

0.0120

0.162

0.00153

0.717 0.717 0.49990

4

2.00 10

0.847 0.0513

0.143 0.144 0.49950 0.00100 0.379 0.0230

0.0717 0.0719 0.49900 0.00200 0.268 0.0163

4

1.44 10

0.00717 0.00736 0.49000 0.0201 0.0858 0.00520

4

1.45 10

0.00358 0.00378 0.48000 0.0405 0.0615 0.0121 0.00373
8

These and other relations between the moduli and Poisson’s ratio may be found in most 

textbooks on elasticity or elastic wave propagation. It is clear from table 1 that the aluminum 

pressure bars are substantially stiffer than the specimen materials. This is particularly true in 

shear. The shear modulus of the specimen materials is 4–5 orders of magnitude smaller than the 

shear modulus of the bars.

The dimensionless parameter G/K, the ratio of the shear to the bulk modulus, is a measure of the 

incompressibility of the material. It is a strictly decreasing function of Poisson’s ratio  and 

approaches zero as  and approaches 1⁄2 from below, as follows from equation 21. Materials for 

which G/K



0.01 are typically referred to as nearly incompressible. From equation 21, we find 

that this condition is equivalent to 



0.495.4 Note that for the specimen, all cases but the last 

two in table 1 satisfy this condition. For nearly incompressible materials, volumetric and shear 

strains of the same magnitude will generate a pressure that is several orders of magnitude larger 

than the shear stress; conversely, pressure and shear stresses of the same magnitude will generate 

a volumetric strain that is several orders of magnitude smaller than the shear strain.

2.4 Wave Speeds and Travel Times

Table 1 also lists three different linear elastic wave speeds for each material: the longitudinal 

wave speed5 cL, the shear wave speed cG, and the “bar wave speed” cE (see section 5.1) The 

subscript indicates which elastic modulus is used to compute the wave speed from the relation 

/ Mc M  

, where M denotes one of the moduli L, G, or E. The table does not list the 

Rayleigh wave speeds, which govern the propagation of surface waves. For all the specimen 

materials considered here, the Rayleigh wave speed is about 95% of the corresponding shear 

wave speed; for the aluminum bars, it is about 93% of the shear wave speed.6

The mechanical impedances corresponding to the longitudinal modulus and the Young’s 

modulus are cL and cErespectively. The ratio of specimen impedance to bar impedance 

governs the transmission and reflection of the loading wave. Table 1 lists the values of this ratio. 

In the one-dimensional theory of wave propagation the different cross-sectional areas of the 

specimen and the bar would also be taken into account, and in such analyses the impedances 

based on the “bar speed” (i.e., cE) are typically used. However, at early times after the arrival 

of the loading wave, the specimen deformation is far from one-dimensional and the longitudinal 

wave speed cL is more relevant, particularly in the interior of the specimen and the bars.

Table 2 summarizes the geometric parameters of the problem and also gives the corresponding 

wave travel times based on the wave speeds in table 1. The travel times (in s) are listed in two 

groups: the time for a wave to travel the length of the bar or the specimen, and the time for a 

 

4 Some authors use the less restrictive condition  

0.49 as the criterion for near incompressibility; this is equivalent to the 

condition G/K



0.02.

5 This is also referred to as the dilatational wave speed.

6 See table 7.5.1 in Eringen and Suhubi (12). 
9

wave to travel the radius of the bar or the specimen. Of course, these travel times depend on the 

type of wave under consideration. Travel times for the lengths are given for the longitudinal, 

shear, and bar waves. Travel times for the radius are given for the longitudinal and shear waves 

only, since the bar wave speed is not relevant in this case.

Table 2. Geometric parameters and wave travel times.

Material

Length

(mm)



Travel Times

(s) Radius

(mm)

Travel Times

(s)

Based

on cL

Based

on cE

Based

on cG

Based

on cL

Based

on cG

Incident Bar

768.0 125.7 153.0 249.6 12.8 2.10 4.16

768.0 125.7 153.0 249.6 9.5 1.56 3.09

Transmission Bar

256.0 41.9 51.0 83.2 12.8 2.10 4.16

256.0 41.9 51.0 83.2 9.5 1.56 3.09

Ballistic Gelatin Specimen 1.45 0.72 93.60 162.1 6.35 3.17 710.0

Linear Elastic Specimen



1.45

0.54

69.93

121.1

6.35

2.37 530.4

 1.71 121.1 7.50 530.4

 3.83 121.1 16.76 530.3

 5.41 121.1 23.69 530.2

 16.91 120.7 74.02 528.6

 23.59 120.3 103.3 526.6

The wave speeds listed in tables 1 and 2 for the Mooney-Rivlin model for the specimen are only 

approximate since they are based on the linear elastic approximation to a nonlinear constitutive 

model. The wave speeds, and hence the impedance ratios and travel times, will change with the 

deformation. Also, since the specimen undergoes large deformations (for both the nonlinear and 

linear elastic models), and since the travel times are based on the undeformed lengths, these 

times are only approximate for the deformed specimen. Nevertheless, the tabulated values 
10

provide an order of magnitude estimate for the wave speeds and travel times in the specimen. 

These estimates for the wave travel times are useful for determining whether certain features of 

the specimen response (to be discussed later) correlate with the arrival of certain waves.

2.5 Measurement of Gap Size

Many of the figures in this report give (for one or both interfaces) the gap size as a function of 

time on the centerline or the gap size as a function of radius at fixed times. In this section, we 

specify the convention used in measuring these gaps. We begin by restricting attention to the 

centerline.

Material points initially on the centerline must remain on the centerline, that is, r  0 if and only 

if R  0. Hence, the radial displacement of points on the centerline is zero. The axial 

displacements of the specimen and the incident bar at their interface (Z  0) on the centerline are 

denoted by uS-IB and uIB, respectively; similarly, zS-IB and zIB denote the corresponding deformed 

axial coordinates. Since the specimen and the bar are initially in contact, they will continue to be 

in contact at some later time t if and only if their axial displacements coincide at time t or, 

equivalently, their deformed axial coordinates coincide at time t. If the specimen and the 

incident bar separate, that is, if a gap forms between them, then the size of this gap, S-IB, is 

given by the difference in their displacements or, equivalently, by the difference in their 

deformed axial coordinates:

S-IB  uS-IB  uIB  zS-IB  zIB . (3)

Similarly, the axial displacements of the transmission bar and the specimen at their interface 

(Z = LS = 1.45) on the centerline are denoted by uTB and uS-TB, respectively; and zTB and zS-TB 

denote the corresponding deformed axial coordinates. If the specimen and the transmission bar 

separate, then the size of the gap between them, S-TB, is given by the difference in their 

displacements or, equivalently, by the difference in their deformed axial coordinates:

S-TB  uTB  uS-TB  zTB  zS-TB . (4)

As will be seen in subsequent figures, when gaps do form it is the result of the specimen moving 

away from the bars, so that for the coordinate system in figure 1, uS-IB  uIB for a gap at the 

incident bar, and uS-TB  uTB for a gap at the transmission bar. The definitions of the gap sizes 

given above were chosen so that the sign of the gap would be positive in these cases; of course, a 

value of zero for the gap implies that the specimen and bar are in contact.

If, as above, we restrict attention to the centerline, then the gap sizes are a function of time only. 

The majority of our plots of gap size are for this case. However, we have also included a few 

plots of gap size versus radial location at selected times. For this case it makes a difference 

whether the axial displacements are referred to points on the specimen and bar with the same 

initial radius R or to points on the specimen and bar with the same deformed radius r. The latter 

case, that is, the axial separation of the specimen and bar along a line of constant radius in the 
11

deformed state, is of more interest and is used here. With this caveat, the definitions above still 

apply.7

2.6 Stress, Stretch and Strain

When the specimen remains in contact with the bars, the deformed thickness (or length) of the 

specimen, S, is determined either from the current positions of its faces or from its initial length 

LS and the current displacements of its faces:

S  zTB – zIB  LS  uTB – uIB . (5)

These variables may be measured at any radial location since the faces of the bars remain very 

nearly planar when the specimen is relatively soft, as in the cases considered here. At radial 

locations where there are substantial gaps between the specimen and the bars, equation 5 will 

overestimate the deformed thickness.

The term stretch denotes the local ratio of deformed to undeformed length, which, of course, 

varies with direction. The stretch in all directions is unity in the undeformed state. A stretch 

greater than 1 corresponds to extension in the given direction; a stretch between 0 and 1 

corresponds to compression. The stretches in the radial, hoop and axial directions are denoted by 

r,  and z, respectively. They are the diagonal components of the deformation gradient tensor 

F and are given by the relations

, , r z

r r z

R R Z

   

 

  

 

. (6)

We are primarily interested in the stretches and strains in the specimen and, in particular, in the 

axial stretch and corresponding axial strain. The nominal (or engineering) strain ez and the 

logarithmic (or true) strain z in the axial direction are given by

ez  1z , z  –ln z . (7)

Both of these strain measures are used in the experimental literature for reporting SHBP test 

data. The radial and hoop strains are defined similarly. Note that these strain measures have 

been taken positive if the particular coordinate direction is in compression. In the general 

continuum mechanics literature, the sign convention for strain is “positive in extension” and 

“negative in compression”. However, in the SHPB literature the opposite convention is typically 

used since the specimen is under axial compression; we have followed that convention here.

The average or mean value of the axial stretch in the specimen is the ratio of the deformed to 

undeformed thickness:

 

7 The deformed axial coordinates are computed at element nodes. The specimen and bar nodes off the centerline do not 

necessarily line up. For a specimen node with deformed radius r, we used linear interpolation between the two incident bar nodes 

whose deformed radius bounded r in order to determine zIB(r); similarly for the transmission bar. 
12 S TB IB

S S

1 z

u u

L L





  

l

, (8)

where the relation on the right follows from equation 5. The mean value of the nominal axial 

strain8 is given by

S IB TB

S S

1 1 z z

u u

e

L L



    

l



. (9)

Note that we use the same notation for local and mean values; it should be clear from the context 

which meaning is intended.

We use a superposed dot to denote the material time derivative. It follows from equation 9 that 

the mean value of the nominal axial strain rate is given by

IB TB IB v v v

z z

S S

e

L L





    , (10)

where vIB and vTB denote the axial component of velocity of the incident and transmission bars at 

the specimen interfaces. The approximation on the right follows from the fact that for soft 

specimens, vIB will generally be much larger than vTB.

9 This approximate relation yields an even 

cruder approximation for the “strain acceleration”: 

IB S v /z z

e L    

. For the imposed velocity

history at the far end of the incident bar, as described in equation 1, vIB eventually oscillates 

around a mean value of roughly twice v0 (see section 5.3).

Note that it is the mean values of axial strain and axial strain rate that are actually measured10 in 

an SHPB test. Similarly, the values of the axial stretch given in our figures are the mean values 

computed from equation 8; consequently, the corresponding nominal axial strains cited in 

subsequent discussions are the mean values determined from equation 9. These mean values 

may differ substantially from the local values if the specimen deformation is highly non-uniform. 

This will certainly be the case during the “ring-up” period when wave propagation dominates. 

Also, at radial locations where there are substantial gaps between the specimen and the bars, the 

relation 8 will overestimate the axial stretch, and consequently, the relation 9 will underestimate 

the axial strain. Likewise, the relations for the mean value of the nominal axial strain rate in 

equation 10 will not be valid when a gap forms.

The Cauchy stress tensor is denoted by . Its components, which give force per unit deformed 

area, are often referred to as true stresses. We follow the sign convention in the general 

 

8 A mean value of the logarithmic axial strain is provided by equation 72 with z given by equation 8.

9 The displacement histories at the S-IB and S-TB interfaces in figures 7 and 8 are consistent with this assertion.

10 More precisely, they are inferred from strain gage measurements in the pressure bars.
13

continuum mechanics literature, that is, the stress components are taken positive in tension and 

negative in compression,11 whereas the pressure p is taken to be positive in compression:

 

1 1

tr

3 3 rr zz p           

. (11)

Here, “tr” denotes the trace, and rr, , and zz are the radial, hoop and axial components 

of Since axisymmetric deformations are imposed and the constitutive models for all materials 

are isotropic,

σθz  σzθ  σrθ  σθr   

However, the shear stress component σrz σzr need not be zero in general; in particular, rz will 

be nonzero in the specimen whenever there is bulging.

Since only isotropic elastic constitutive relations are used here, the principal axes of stress and 

strain coincide. The  (i.e., hoop) direction is always a principal axis of stress and strain for 

axisymmetric deformations, but the r and z (i.e., radial and axial) directions will be principal 

axes if and only if the shear stress rz is zero. In this case, the radial, hoop and axial stretches are

the principal stretches, and the Jacobian12 of the deformation, J, is given by

J  rz. (13)

There are two situations where the condition rz  0, the relation 13, and the relations

σrr σθθ  er  eθ  

r

z

J

 



 

 

can be guaranteed to hold exactly. One of these occurs when the deformation is uniform, which 

is the desired state in an SHPB test. The other occurs at any point on the centerline of a solid 

specimen, regardless of whether or not the deformation is uniform.13 Furthermore, for the nearly 

incompressible specimen materials considered here, J will generally be close to one even for 

large deformations, in which case equation 143 yields the approximation

1

r

z

 



 

(15)

 

11 This is also the sign convention used in LS-DYNA. However, in the SHPB literature it is more common to take stress 

positive in compression, since the specimen is under axial compression during the useful part of the test. 

12 This is the determinant of the deformation gradient F and represents the local ratio of deformed to undeformed volume.

13 That rr  and rz  on the centerline follow from the radial and axial components, respectively, of the momentum 

balance equation, after multiplying by r and taking the limit as r approaches zero. That r  on the centerline follows from 

equations 6, since r/R approaches 

  r R as R and r approach zero. Also, both   r Z

and 

  z R

must be zero on the 

centerline to prohibit singularities. 
14

on the centerline. When the deformation in the specimen is non-uniform, the relations 13-15 will 

hold approximately at points near the centerline, and they may even hold approximately at points 

substantially off the centerline.

3. The Mooney-Rivlin Constitutive Model for the Specimen

3.1 Incompressible and Compressible Versions

The classical Mooney-Rivlin model is a nonlinear elastic constitutive model for isotropic, 

incompressible solids. It is one of the earliest models of this type, and was originally applied to 

the response of rubber under large deformations; see the original papers by Mooney (13) and 

Rivlin (14) and the book by Treloar (15). Because of its simplicity, the model is widely used in 

theoretical studies and is often applied to nearly incompressible materials other than rubber. The 

Mooney-Rivlin model is discussed in many textbooks on nonlinear elasticity and continuum 

mechanics; see (16–19). A description of the model is provided in appendix A-2.

The incompressibility assumption is, of course, an idealization. For quasi-static problems with at 

least one free surface, the assumption of incompressibility typically yields reasonable results if 

the bulk modulus is at least two orders of magnitude larger than the shear modulus. However, 

the incompressibility constraint is inconsistent with the propagation of longitudinal waves, so for 

dynamic problems involving nearly incompressible materials, some degree of compressibility 

must be added to the constitutive model.

The simulations described in sections 6 and 7 and in appendices B and C used LS-DYNA’s 

compressible version of the Mooney-Rivlin model for the specimen; this is Material Type 27 in 

the manuals (8, 9). The model has three material constants (aside from the density): the elastic 

moduli14 A1 and A2, and the Poisson’s ratio . The two moduli characterize the response to 

volume-preserving deformations (just as in the incompressible version), whereas the nonlinear 

pressure-volume relation involves all three constants. We refer to this as the “compressible 

Mooney-Rivlin model”; it is discussed in more detail in appendix A-3. Verification of the model 

implementation is discussed in appendix A-4. For the purposes of the present discussion, we 

simply note that the shear modulus G in the small strain, linear elastic approximation to the 

model is related to the moduli A1 and A2 by

G  2(A1  A2) . (16)

Since the linear elastic bulk modulus K is given in terms of the shear modulus and Poisson’s 

ratio by equation 21,  can be adjusted so as to yield the bulk modulus for the material in 

question, once A1 and A2 (and hence G) have been determined. For the nearly incompressible 

materials of interest here, this means choosing  very close to 1⁄2.

 

14 The LS-DYNA manuals (8, 9) use the symbols A and B for the constants A1 and A2, respectively.
15

3.2 Calibration of the Mooney-Rivlin Model for Ballistic Gelatin

Ballistic gelatin is a prime example of a soft, nearly incompressible material. Gelatin powder is 

an aggregate of large protein molecules (collagen) of various sizes, extracted from animal tissues 

(20). In this report, the term gelatin refers to the solution of this powder in water. The common 

formulations of ballistic gelatin are 20% and 10% powder by mass (80% and 90% water by 

mass). Both formulations are used as tissue surrogates in impact and penetration tests.

LS-DYNA’s compressible Mooney-Rivlin model was calibrated to yield a response typical of 

20% ballistic gelatin. The parameters A1 and A2 listed in table 3 were chosen to give rough 

agreement with quasi-static uniaxial compression data on ballistic gelatin obtained in an 

experimental study by Moy et al. (5). More precisely, these values give a stress-strain curve that 

lies slightly above their curve for a strain rate of 1/s.15 This represents a rough extrapolation of 

their quasi-static data (strain rates of 0.001/s, 0.01/s, 1/s) to strain rates on the order of 2500/s.16 

By equation 16, the corresponding small strain shear modulus is G = 80 kPa. Substantially lower 

values for these moduli would be expected for 10% gelatin.

Table 3. Parameters for the compressible Mooney-Rivlin Model.

 (g/cm3) A1 (kPa) A2 (kPa) 

1.00 12.0 28.0 0.49999

Since gelatin is mostly water, one would expect its density and (initial) bulk modulus to be close 

to that of water, namely, 1 g/cm3and 2.3 GPa, respectively. The room temperature densities for 

20% gelatin reported in Aihaiti and Hemley (10) and Winter and Shifler (20) are about 

1.00 g/cm3and 1.06 g/cm3, respectively.17 We used the value 1.00 g/cm3. However, the 

estimated bulk modulus for 20% ballistic gelatin in Aihaiti and Hemley (10) is more than twice 

that of water.18 We chose a Poisson’s ratio of   0.49999, which together with the shear 

modulus of 80 kPa (calculated above) yields a bulk modulus of K = 4.0 GPa. Note that this is 

more than five orders of magnitude larger than the shear modulus.

 

15 The corresponding stress-stretch curve is plotted in figure A-1 in our appendix A-4. As discussed their, these curves are 

largely insensitive the value of  provided that it is close to 1⁄2.

16 This extrapolated curve lies well below their stress-strain curve at a strain rate of 2500/s, obtained in an SHPB test on an 

annular specimen. We believe that there are large inertial effects in that test and, in particular, that the stress state is not uniaxial 

and so not comparable with the quasi-static data. The quasi-static data referred to here was provided by Tusit Weerasooriya of 

ARL. This data is not included in (5), except for the upper stress-strain curve in figure 6(a) in that report, which is described 

there as corresponding to a strain rate of 1/s but which appears to be the 0.001/s curve in the data provided to us.

17 The value 1.00 g/cm3

is based on a slight extrapolation of the data given in Aihaiti and Hemley (10); the authors attribute 

this data to Dana Dattlebaum at Los Alamos National Laboratory (LANL).

18 The value reported in Aihaiti and Hemley (10) varied from about 4.6 to 5.2 GPa at room temperature, depending on the 

functional form of the equations used to fit the data. However, it appears that even the lowest of these values may be too high. 

After are computational study was completed, we were informed that the data in Aihaiti and Hemley (10) is inconsistent with the 

shock Hugoniot data from LANL, and that the latter implies an initial bulk modulus much closer to that of water (21).
16

We emphasize that our intention in this study was neither to develop a constitutive model for 

ballistic gelatin nor to simulate a particular SHPB experiment on gelatin. Ballistic gelatin is 

known to exhibit viscoelastic behavior. The rate-dependence of the large strain compression 

data in (5) and elsewhere is consistent with this, as is the small strain rheometer data in Juliano et 

al. (22). An accurate model for ballistic gelatin would require the addition of a viscoelastic 

component to the constitutive model. Consequently, numerical simulations capable of 

duplicating high-rate SHPB experiments would likely require a properly calibrated nonlinear 

viscoelastic model. The use of a purely elastic model was motivated by our original objective of 

quantifying inertial effects in SHPB tests. Since there are no strain rate effects in an elastic 

model, any differences observed between quasi-static and high-rate uniaxial compression tests 

must be due to inertial effects in the latter.

4. Computational Parameters

4.1 Geometry and Meshing

The meshes for the specimen and pressure bars consisted of 4-node quadrilateral elements used 

in 2-D axisymmetric mode. Single-point integration was used throughout. Three different 

meshes were used for the specimen, as shown in figure 2. Since the axial strain in the specimen 

is approximately twice the radial strain, it was beneficial to use specimen elements whose initial 

length in the axial direction was twice that in the radial direction, so that the element became 

more nearly square under severe axial compression. The default or baseline mesh is shown in 

figure 2a and referred to as the “25×50 m mesh”: all specimen elements had a radial edge 

length of 25 m and an axial edge length of 50 m. This specimen mesh was used in all 

simulations with the exception of the mesh sensitivity study discussed in appendix B-1, where 

the results for the baseline mesh are compared with the results for a coarser mesh (50×100 m in 

figure 2b) and a finer mesh (12.3×12.5 m in figure 2c)
17

Figure 2. Meshes for the Incident Bar (IB), Transmission Bar (TB), and Specimen (S). The 

baseline mesh for the specimen is shown in case (a).

All incident bar and transmission bar elements were squares with a 100-m edge length. 

Differences between bar and specimen element sizes in the axial direction could result in a 

computational impedance mismatch if the specimen and bars had similar mechanical 

impedances. This was not an issue here since the specimen and bar impedances were already 

quite dissimilar (see table 1). Also, since the faces of the stiffer bars are expected to remain 

nearly flat, elements with larger radial edge length could be tolerated. The element size in the 

bars does need to be small enough to capture the essential features of the wave propagation (see 

section 5). The 100-m edge length seemed to be sufficient for this.

4.2 Boundary Conditions

The interaction of the specimen and the bars at the S-IB interface and the S-TB interface was 

governed by the LS-DYNA segment-based penalty algorithm CONTACT_2D with the option

AUTOMATIC_SURFACE_TO_SURFACE (8, 9). This algorithm requires three input

parameters: SFACT, a scale factor for the penalty force stiffness; VDC, a viscous damping 

coefficient (in percent of the level for critical damping); and a Coulomb friction coefficient, 

which was set to zero for frictionless contact (see section 2.2). The LS-DYNA default values 
18

SFAC = 1 and VDC = 10 were used for all simulations with the exception of the contact 

algorithm parameter sensitivity study in appendix B-2.

As discussed in section 2.2, the boundary condition at the impact end of the incident bar is a 

prescribed axial velocity history. No radial constraint was placed on the nodes there.

At the far end of the transmission bar, we used the LS-DYNA non-reflecting boundary condition 

BOUNDARY_NON_REFLECTING_2D (8, 9) to simulate a semi-infinite transmission bar. 

This algorithm imposes normal and shearing stresses given by

normal normal vL     c , 

shear tangential vG     c . (17)

Here cL and cG are the longitudinal and shear wave speeds, respectively, of the transmitting 

material (in our case aluminum), and vnormal and vtangential are the particle velocities in the normal 

and tangential direction, respectively. We conducted some preliminary tests on this boundary 

condition and found that it performed extremely well in one-dimensional uniaxial strain, and 

reasonably well for waves of the type generated in our simulations. Once the loading wave 

travels through the specimen and reaches the transmission bar, it takes about 42 s for a 

longitudinal wave and 51 s for the bar wave to travel the length of the transmission bar (see 

table 2), and then another 42–51 s for any slight reflections from the supposedly non-reflecting 

boundary to arrive back at the specimen. Thus any influence that slight reflections from the far 

end of the transmission bar might have on the interactions at the S-TB interface could not 

possibly be felt earlier than 84–102 s after the loading wave first reached that interface.

4.3 Miscellaneous

The default values, 1.5 and 0.06, were used for the quadratic and linear artificial viscosity 

coefficients, respectively. The time step factor of safety, TSSFAC, was set to 0.1, meaning that 

the time step computed using the Courant condition is then multiplied by 0.1.

LS-DYNA’s Orthotropic Elastic model (Material Model 2 in LS-DYNA [8] and Hallquist [9]) 

was used for the pressure bars. The material constants for this orthotropic, linear elastic model 

were chosen in such a way as to yield an isotropic, linear elastic model. Values of the material 

constants appropriate for aluminum were used in all the simulations and are given in table 1.

5. Wave Propagation in the Pressure Bars

5.1 General Remarks

The loading wave in the incident bar generates an axial stress on the order of cEv0, where v0 is 

the plateau of the imposed velocity at the impact end, as described in equation 1. For the value 

v0  1.8125 m/s used in our simulations this gives an axial stress on the order of 25 MPa for 
19

aluminum, which is well below the dynamic yield stress.19 This justifies the use of a linear 

elastic model for the pressure bars in our simulations. Linear elastic deformation of the bars is 

also a requirement in actual SHPB tests, since the analysis of test data is based on linear elastic 

wave propagation theory (1, 2).

A good reference for linear elastic wave propagation in bars is Chapters 2 and 8 in Graff (24); 

see also Eringen and Suhubi (12). The one-dimensional theory of wave propagation in long bars 

predicts that a longitudinal elastic wave travels undistorted at the “bar speed” cE, that is, at the 

speed governed by the Young’s modulus E. The more accurate results of the three-dimensional 

theory are much more complex, particularly near the loading end where the wave becomes 

highly distorted with large oscillations. These oscillations decay somewhat as the wave 

progresses down the bar, but they do not go away entirely. Near the centerline of the bar the 

wave initially propagates at the longitudinal elastic wave speed cL. But release waves from the 

stress-free lateral surface interact with the longitudinal wave, so that by the time the wave has 

traveled 10–20-bar diameters only a trace of the initial longitudinal wave remains. Most of the 

energy is contained in a wave moving at the slower bar speed cE. Due to wave dispersion, the 

rise time of the loading pulse increases as it moves down the bar and the axial stress, strain and 

particle velocity oscillate about the value predicted by the one-dimensional theory.

According to table 2, for the material properties and incident bar length considered here, the 

longitudinal and bar waves should take about 126 and 153 s, respectively, to arrive at the 

specimen. The average of these arrival times is about 140 s. As will be seen in section 5.3, the 

velocity histories in the incident bar in the simulations are consistent with the preceding 

observations.

5.2 Considerations Governing the Bar Lengths

In an SHPB experiment the length of the striker bar determines the duration of the loading pulse, 

which in turn determines the duration of the test. The length of the incident bar is determined by 

the requirement that the incident and reflected pulses measured by a strain gage mounted 

midway down the bar must not overlap. As a result, incident bar lengths of several meters are 

commonly used.

In order to reduce the size of the output files and particularly the run times of the simulations, we 

did not wish to model the full length of a real incident bar. Indeed, there was no need to do so, 

as the prescribed velocity history (equation 1) is not a pulse but instead has essentially infinite 

duration, and there is no strain gage to consider. Nevertheless, two considerations govern the 

minimum length required for the incident bar in our numerical simulations. First, a bar of length 

at least 20 diameters is desirable in order to simulate the effects of dispersion on the loading 

wave. Second, in the simulations the total strain that can be imposed on the specimen is 

proportional to the length of incident bar. To see this, note that the compressive incident wave 

 

19 The value of the Hugoniot elastic limit for 7075-T6 aluminum listed in Steinberg (23) is 420 MPa.
20

reflects from the lower impedance specimen as a tensile wave, travels back down the bar, reflects 

from the loading end as tensile wave,20 and eventually arrives back at the specimen, after which 

time the simulation is no longer useful. The time for this round trip is thus twice the time to 

travel one length of the bar, and this time window limits the total strain that can be imposed on 

the specimen.

We chose an incident bar length of LIB = 768 mm, which gives a length-to-diameter ratio of 30.0 

for the larger diameter bar and a time window of roughly 2



140 s = 280 s, which suffices for 

large strains in the specimen. Since the incident wave arrives at the specimen at about 

t  140 s, specimen loading would occur from about t  140 s to t  420 s. However, when 

gaps formed (as they did in most of our simulations), they developed shortly after the arrival of 

the incident wave. And since this became the focus of our study, we terminated most of the 

simulations at 340 s. Thus most of our time history plots are restricted to the 200-s time 

window from t  140 s to t  340 s.

As discussed at the end of section 4.2, slight reflections from the supposedly non-reflecting 

boundary at the far end of the transmission bar could not possibly be felt at the S-TB interface 

earlier than 84–102 s after the loading wave first reaches that interface. Using t  140 s for 

the latter arrival time,21 we see that the effects of any deficiencies in the non-reflecting boundary 

condition cannot be felt earlier than about t  230 s. As will be seen below, gaps formed at the 

S-TB interface well before this time, although we cannot rule out the possibility of some 

influence at later times.

5.3 Velocity Histories in the Incident Bar

Figures 3–6 are plots of the time history of the axial component of the particle velocity, vz, at 

four locations on the centerline of the incident bar (r  R  0). The velocity is scaled by the 

plateau value (v0 = 1.8125 m/s) of the input velocity. All plots are for the bar with the larger 

radius (RB = 12.8 mm), although the results for the smaller radius should be similar.

Figures 3 and 4 give velocity histories at the loading end (Z  LIB), 1/3 of the way from the 

loading end (Z  0.666 LIB), and 2/3 of the way from the loading end (Z  0.333 LIB). On 

these two plots the time is scaled by the transit time of a longitudinal wave along the length of 

the bar, which is LIB /

Al

L

c  125.7 s, where 

Al

L

c

is the longitudinal wave speed of aluminum. 

Figure 3 is for an initial rise time tR of 1 s, which is too short for the linear ramp to be 

distinguishable on the time scale of this plot. Figure 4 is for the 25-s initial rise time; in this 

case the linear ramp is clearly distinguishable. For both cases, the velocity at the two interior 

locations oscillates about the plateau value v0. 

 

20 The fixed velocity boundary condition at the loading end produces the same effect on the stress wave as a rigid boundary 

would, that is, a tensile wave is reflected as a tensile wave. Indeed, relative to an observer moving with the imposed velocity v0, 

that boundary would appear rigid.

21 Note that the travel time through the specimen is less than 1 s for the Mooney-Rivlin model as well as for the linear 

elastic model with   0.49999 (see table 2).
21

Figure 3. The scaled axial velocity vs. non-dimensional time at three locations on the centerline of the 

incident bar: the loading end with a 1-s rise time, 1/3 of the way from the loading end, and 2/3 

of the way from the loading end (blue curve).
22

Figure 4. The scaled axial velocity vs. non-dimensional time at three locations on the centerline of the 

incident bar: the loading end with a 25-s rise time, 1/3 of the way from the loading end, and 2/3 

of the way from the loading end (blue curve).
23

Figure 5. The scaled axial velocity history on the centerline of the incident bar at the S-IB interface for a 

1-s (initial) rise time.

Figures 5 and 6 give the velocity history on the centerline of the incident bar at the S-IB interface 

for the 1- and 25-s initial rise times, respectively.22 In these two figures the time is given in s 

for later comparison with other plots. Note that the particle velocity ramps up to nearly double 

its initial value plateau v0, as would be expected from reflection of the wave at a free surface. 

This is consistent with the small specimen/bar impedance ratio (see table 1). This doubling of 

the velocity on reflection is also evident at later times for the two interior locations in figures 

3 and 4. The Mooney-Rivlin model was used for the specimen in these simulations, but in view 

of the low impedance ratios for all specimen models considered here, only small differences in 

the particle velocity histories in the incident bar after reflection would be expected for the other 

specimen models.

 

22 The incident bar velocity vz

in these figures was denoted by vIB in equation 10. 
24

Figure 6. The scaled axial velocity history on the centerline of the incident bar at the S-IB interface for a

25-s (initial) rise time.

A comparison of figures 3 and 5 (1-s initial rise time) with figures 4 and 6 (25-s initial rise 

time) shows that the computational “pulse shaping” has substantially smoothed the loading wave 

in the latter case. Also of interest is the increase in the final rise times of these waves due to 

dispersion in the bar. This final rise time could be defined as the difference between the time at 

the first peak in velocity at the S-IB interface and the arrival time of the wave at the S-IB 

interface. However, precise arrival times are difficult to determine from the plots in figures 

5 and 6. A close inspection of figure 5, for example, reveals barely discernable peaks in the

oscillatory toe of the velocity history starting at about 135 s, although there is no substantial 

increase in velocity until about 140 s. Thus, for the 1-s initial rise time, most of the final rise 

in particle velocity at the S-IB interface occurs from 140 to 164 s, that is, over a 24-s time 

interval. For the 25-s initial rise time, most of the final rise in particle velocity at the S-IB 

interface occurs from 143 s to 184 s, that is, over a 41-s time interval.
25

Of the three curves in figures 3 and 4, the blue curves correspond to a point in the incident bar 

closest to the specimen, namely, 1/3 of the bar length back from the S-IB interface. In the 

following discussion, we consider the velocity history at this interior location at early times, that 

is, prior to the arrival of the reflected wave at a non-dimensional time of about 1.6. This 

represents the loading wave shortly before it arrives at the S-IB interface. For a given initial rise 

time, we compare this loading wave with the velocity history at the S-IB interface over a time 

interval that includes the first two major peaks in velocity at these two locations. Clearly, these 

velocity histories are qualitatively similar.23 In particular, for the 25-s rise time, a close 

inspection of the plots reveals that both curves have three inflection points prior to the first peak 

in velocity. For the 1-s rise time, the situation is more complicated due to the small oscillations 

superimposed on the main curves; but if these are neglected, both curves would appear to have 

only one inflection point between the toe and the first peak. These inflection points in the 

velocity history correspond to extreme values in the acceleration of the bar, which turn out to be 

important for interpreting the stress histories in the specimen. Furthermore, the fact that these 

features are also exhibited prior to the arrival of the wave at the S-IB interface implies that they 

are not a consequence of interaction with the specimen, but instead are imposed upon the 

specimen.

Finally, we note that the axial velocity and axial acceleration of the specimen at the S-IB 

interface will coincide with the axial velocity and acceleration of the incident bar as long as the 

specimen and bar are in contact, but when a gap forms at this interface the velocities and 

accelerations of the bar and the specimen will generally be different.

6. Simulations with the Moony-Rivlin Model and a 1-s Rise Time

All the simulations in this section are for the loading wave with a 1-s initial rise time and the 

compressible Mooney-Rivlin model for the specimen. Conventions pertaining to the gaps are 

discussed in section 2.5. The results for the larger radius bar (the default radius RB = 12.8 mm in 

this report) are presented in subsections 6.1–6.3. These results are compared with those for the 

smaller radius bar in subsection 6.4.

6.1 Results on the Centerline

All results in this subsection pertain to points on the centerline at the S-IB or the S-TB interface. 

We begin with some remarks about the figures and a few general observations, and then proceed 

to discuss the details of gap formation on the centerline.

 

23 Of course, they differ quantitatively by a factor of nearly 2 due to a near doubling of the particle velocity on reflection of 

the wave at the S-IB interface. Also, as will be demonstrated in sections 6 and 7, the large velocity spikes at later times in figures 

5 and 6 coincide with the closing of gaps on the centerline at the S-IB interface.
26

6.1.1 General Remarks

Figure 7 shows the time histories of the axial displacement of the incident bar (uIB) and the 

specimen (uS-IB) at the S-IB interface, as well as the history of the corresponding gap size 

(S-IB  uS-IB  uIB). Similarly, figure 8 shows the time histories of the axial displacement of the 

transmission bar (uTB) and the specimen (uS-TB) at the S-TB interface, as well as the history of the 

corresponding gap size (S-TB  uTB  uS-TB). In all of the figures in this report, the displacements 

and gaps are measured in microns. The histories of the gaps at the two interfaces are compared 

in figure 9. Figure 10 enlarges the plot in figure 9 in the vicinity of the time at which the gaps 

initially form. From figures 7 to 10, it is clear that gaps first form at both interfaces at about 

165 s. These gaps open and close several times, and they do so more frequently at the S-TB 

interface.

Figure 7. The histories of the axial displacement of the incident bar and the specimen on the centerline at 

the S-IB interface. Also shown is the history of the corresponding gap size (blue).
27

Figure 8. The histories of the axial displacement of the transmission bar and the specimen on the centerline 

at the S-TB interface. Also shown is the history of the corresponding gap size (blue).
28

Figure 9. A comparison of the histories of the gap sizes on the centerline at the S-IB and S-TB interfaces. 

Also shown is the history of mean axial stretch in the specimen (blue).
29

Figure 10. A comparison of the histories of the gap sizes on the centerline at the S-IB and S-TB interfaces 

in the vicinity of the time at which they first form. Also shown is the history of mean axial 

stretch in the specimen (blue).

The largest gaps observed were 43 m at the S-IB interface and 38 m at the S-TB interface. 

Since these maximum gap sizes are smaller than the initial 50-m axial dimension of the 

specimen elements, the issues of mesh dependence and dependence on the contact algorithm 

naturally arise. These issues are addressed in appendix B. The results described there support 

the conclusion that the gap phenomenon is not a numerical artifact introduced by the contact 

algorithm or stemming from an insufficiently fine mesh. Additional support for the latter

conclusion will be given in section 8 (linear elastic model for the specimen), where gap sizes are 

found to exceed the mesh size.

Since the loading wave in the incident bar arrives at the S-IB interface at (or slightly before) 

t  140 s, whereas the gaps first form about 25 s later, gap formation does not coincide with 

arrival of the loading wave. Nor does gap formation coincide with the first wave reflections 

from either interface. Indeed, since the travel time of an axial longitudinal wave across the 

specimen is about 0.7 s (see table 2), there is time for about 18 round trips of this wave before 
30

the gaps form. In any case, an initially compressive longitudinal wave in the specimen should 

reflect from the transmission bar (and subsequently from the incident bar) as a compressive wave 

due to the higher bar impedance, so tensile stresses would not be expected to develop from these 

wave reflections. On the other hand, it is clear from table 3 that longitudinal unloading waves 

from the lateral (stress-free) surface of the specimen should arrive at the centerline at both 

interfaces about 3–4 s after the loading wave in the bar first reaches the specimen, but this does 

not correlate with gap formation either. Finally, it is clear from table 2 that 25 s is insufficient 

time for a shear wave to propagate more than a fraction of the thickness or the radius of the 

specimen; consequently, these waves cannot influence the formation of gaps. The possible 

influence of wave propagation in the bars on the formation of gaps is discussed in section 6.3.

The history of the mean axial stretch z (see equation 8) is plotted on the right axis in figures 

9 and 10. The slope of the axial stretch curve in figure 9 is nearly constant for t



200 s, and 

we estimate that the nominal strain rate for these times is between 2470/s and 2480/s. This is 

within 1% of the desired 2500/s strain rate.

Figure 11 plots the histories of the axial stress zz in the specimen and the incident bar at points 

near the centerline and the specimen-bar interface; the corresponding gap size is also shown. 

Figure 12 is an analogous plot at the transmission bar interface. The stress is actually measured 

at the centroids of the specimen and bar elements adjacent to the centerline and the interface. 

Since the coordinates of these centroids do not coincide, slight differences in the values are not 

unexpected. Also, recall (section 2.6) that the axial stress is taken to be negative in compression. 

When a gap opens between the faces of the specimen and the bar, the axial stress should drop to 

zero on both faces, regardless of the size of the gap. The plots in figures 11 and 12 are consistent 

with this expectation, keeping in mind that these stresses are actually measured at a slight 

distance away from the faces.
31

Figure 11. Histories of the axial stress (negative in compression) and the gap size on the centerline at

the S-IB interface. The stress is measured at the centroids of the specimen and incident bar 

elements adjacent to the centerline.
32

Figure 12. Histories of the axial stress (negative in compression) and the gap size on the centerline at 

the S-TB interface. The stress is measured at the centroids of the specimen and 

transmission bar elements adjacent to the centerline.

A comparison of the velocity history at the S-IB interface in figure 5 with the displacement, gap, 

stress and stretch histories in figures 7–12 leads to the following additional observations. For 

convenience of discussion we have divided the time histories into several stages.

6.1.2 The Initial Loading of the Specimen: Stage I (140–165 s)

This first stage covers the initial loading of the specimen prior to gap formation. It begins with 

the first significant increase in the axial velocity of the incident bar at the S-IB interface at about 

140 s (see figure 5)24 and ends with the onset of the gaps at 165 s. From the axial stretch 

histories in figures 9 and 10, we see that the axial stretch decreases from 1 to about 0.97 over this 

25-s time interval; hence, the nominal axial strain increases from 0% to about 3%.

The first observable increase in axial stress at the S-IB interface occurs at about 143 s (see 

figure 11). The stress increases in magnitude to a peak value of 4.0 MPa at around 158 s, 

which is 7 s prior to the formation of a gap at 165 s. After this peak, the stress decreases more 

rapidly than it had risen and drops to zero at about 164 s. Superimposed on this behavior are 

 

24 However, as observed at the end of section 5.3, the loading wave may have arrived as early as 135 s.
33

smaller oscillations with a period of 1–2 s; similar oscillations can be seen in the velocity 

history at the S-IB interface (figure 5). The history of the axial stress at the S-TB interface is 

qualitatively and quantitatively similar, with the times delayed by about 1 s (see figure 12). In 

particular, the axial stresses at the two interfaces are approximately equal throughout this stage, 

from which we may infer that zz is approximately uniform through the thickness of the 

specimen. In the experimental literature, this condition is referred to as dynamic equilibrium.

25

At t  158 s (the instant of the stress peak in figure 11), the axial stretch z is about 0.99, as can 

be inferred from figure 9. This stretch corresponds to a compressive nominal strain of 0.01, that 

is, 1% axial strain. A glance at the stress-stretch curve for uniaxial stress in figure A-1 reveals 

that when z = 0.99, the axial stress zz is a very small fraction of 1 MPa. Indeed, from equation 

A-13 and the values of A1 and A2 in table 3, we find that zz  2.4 kPa for a state of uniaxial 

stress when z = 0.99.26 This is three orders of magnitude smaller than the peak value of 

zz  4.0 MPa in figure 11, which (as previously noted) occurs at z = 0.99. Consequently, the 

stress state at the instant of this peak (t  158 s) is nowhere close to the desired state of uniaxial 

stress. In fact, as established in the next paragraph, at this instant the stress state on the 

centerline is very nearly hydrostatic.

Recall from section 2.6 that at points on the centerline, the radial, hoop and axial directions are 

principal axes of stress and strain, the relations 14 hold exactly, and equation 15 holds 

approximately (for small volume changes, i.e., 

J 1

). As discussed in appendices A-2 and A-3, 

the stress difference zz  rr on the centerline is given exactly by equation A-14 when J  1 and 

approximately by that equation when 

J 1

. When z = 0.99, equation A-14 yields the value

2.4 kPa.27 Summarizing, we have zz  rr  2.4 kPa and   rr exactly, whereas 

zz  4.0 MPa. Hence, the differences between the principal stresses rr, , and zz are either 

zero or at least three orders of magnitude smaller than the stresses themselves. It follows that all 

three principal stresses are nearly equal to each other and hence, by equation 11, nearly equal to 

–p. The volumetric strain 1J at this peak pressure is approximately 0.001, that is, 0.1%.28

The preceding quantitative estimates apply only at the stress peak on the centerline. However, at 

neighboring times and at radial locations near the centerline we arrive at the same qualitative 

conclusion, namely, that the stress state is nearly hydrostatic for most of this stage (certainly 

 

25 It should be kept in mind that the stresses in figures 11 and 12 are measured on the centerline. On the other hand, only the

mean axial stress on an interface can be inferred from experimental data, so it is the equality of the mean axial stresses at the two 

interfaces that is implied by the term dynamic equilibrium. 

26 As discussed in appendix A-3, equation A-13 holds approximately for uniaxial strain when 

J 1

(i.e., small volume 

changes). Also, note that the axial strain is small enough that the exact relation in equation A-13 and the linear elastic 

approximation give essentially the same result.

27 Since the right-hand side of equation A-14 equals the right-hand side of equation A-13, this value follows from the 

calculation in the preceding paragraph.

28 This follows from setting p equal to zz and using equation A-21 with K = 4.0 GPa (see section 3.2), which yields 

1J = 4.0 MPa/4.0 GPa.
34

during the time interval 147–163 s in which |zz| exceeds 0.4 MPa), and the axial stress exceeds 

that for a uniaxial stress state by several orders of magnitude.

From figure 5, we estimate that the peak slope of the velocity-time curve during Stage I also 

occurs around 158 s.29 Thus, the peak in the axial acceleration of the incident bar at the S-IB 

interface coincides with the peak in the axial stress on the centerline, and as that acceleration 

decreases towards zero so does the axial stress. The axial acceleration of the incident bar 

induces an axial acceleration in the specimen, which in turn induces a radial acceleration in the 

specimen. Therefore, it seems likely that during this stage the value of the axial stress and the 

nearly hydrostatic stress state in the neighborhood of the stress peak are inertial effects 

associated with the accelerations imposed on the specimen by the incident bar.30 Furthermore, 

since the specimen is approximately in dynamic equilibrium (at least along the centerline), these 

would appear to be primarily radial inertia effects. Song et al. (6) also observed a correlation 

between axial acceleration and axial stress in the specimen at early times in their SHPB tests and 

numerical simulations on gel specimens and likewise concluded that this was a radial inertia 

effect. This conclusion is also consistent with various theoretical analyses of inertial effects 

applied to soft materials (26–33). We will return to these issues in section 7, where we discuss 

the case with a 25-s initial rise time.

6.1.3 Gap Opening and Closure (Stages II–VI)

Stage II (165–280 s): This stage begins with the opening of gaps at both interfaces and ends 

with the first closing of the gap at the S-IB interface. Over this 115-s time interval, the nominal 

axial strain increases from 3% to 32% (see figure 9).

From figure 10, we see that a gap begins to form at the S-TB interface at about 165.0 s and at 

the S-IB interface about 0.1 s later. Prior to this time, S-IB is negative, which means that the 

specimen and incident bar have interpenetrated. The maximum amount of this interpenetration is 

0.06 m, which is 0.12% of the initial axial edge length of a specimen element. The penaltybased contact algorithm allowed this degree of interpenetration; it is also evident in several 

subsequent figures. This interpenetration begins to decrease at about 164.8 s, which is 0.2 s 

prior to the instant at which the gap forms at the S-TB interface. The axial stretch at the onset of 

the gaps is about 0.971. The corresponding nominal axial strain is 0.029, that is, very close to

3%.

From figure 5, we see that the peak value of the axial velocity of the incident bar occurs at about 

164 s. The axial acceleration of the bar is zero at this instant and negative for some time 

thereafter. The axial stress at both interfaces also reduces to zero at about 164 s (figures 11 and 

12). Hence, both the axial stress and the axial acceleration at the S-IB interface drop to zero at 

about the same time, and a gap forms at that interface about 1 s later, that is, about 1 s after 

 

29 Here, we disregard the smaller oscillations superimposed on the main curve.

30 In this regard, see the comments in the next to last paragraph in section 5.3.
35

the incident bar begins to decelerate. That a gap forms just after the axial stress drops to zero is 

not surprising. The important observation here is that this occurs when the axial acceleration of 

the S-IB interface passes through zero. We will return to this point in section 7.

With the exception of some minor oscillations, the incident bar continues to decelerate until 

about 173 s and then re-accelerates until 180 s (figure 5). However, the S-IB gap does not

close when the incident bar re-accelerates. The bar velocity goes through many oscillations 

about its mean value of about 2v0 before the S-IB gap first closes at 280 s. The gap at the S-TB 

interface is also closed at 280 s, although it has opened and closed several times prior to this. 

Nevertheless, a gap exists at the S-TB interface for most of this stage. Observe the minor spikes 

in stress at the S-TB interface when the gap closes there. Note that the magnitude of the S-TB 

gap never exceeds 4 m in this stage, whereas the S-IB gap reaches a maximum of 38 m at 

about 226 s.

Stage III (280–286 s): The gaps at both interfaces close at 280 s and remain closed, for the 

most part, over this 6-s time interval. The S-IB gap does re-open momentarily at the beginning 

of this stage, and sub-micron gaps re-open momentarily at both interfaces near the end of the 

stage. The first large spike in the velocity at the S-IB interface (figure 5) occurs at the beginning 

of this stage, and large spikes in the stress (figures 11 and 12) occur throughout this stage.

Stage IV (286–324 s): The gaps re-open and remain open at both interfaces over this 30-s 

time interval. The axial stresses near the interfaces again drop to nearly zero. The gaps are now 

roughly the same size at both interfaces. The S-IB gap reaches a peak value of 43 m at 308 s; 

the S-TB gap reaches a peak value of 38 m about 4 s later. These are the maximum gap sizes 

achieved in this simulation. During this stage, the nominal axial strain increases from about 33% 

to 43%.

Stage V (324–333 s): The gaps at both interfaces close and remain closed, for the most part, 

over this 9-s time interval. This gap closure results in the second large spike in the velocity at 

the S-IB interface (figure 5) and the second group of large spikes in the stress (figures 11 and 12) 

as the specimen “slaps” the bars.

Stage VI (333–340 s): The gaps re-open and remain open at both interfaces. The axial stresses 

near the interface again drop to nearly zero. At 340 s, the axial stretch is 0.535 and the nominal 

axial strain is 46.5%. Thus, gaps persist to large strains at both interfaces.

6.2 Results at Other Radial Locations

Figure 13 shows deformed mesh plots in the vicinity of the specimen at 236 s (Stage II). The 

centerline lies along the bottom of each plot. Thus, figure 13a shows the entire (deformed) 

specimen mesh in the 2-D axisymmetric simulation; note the bulging at the free surface of the 

specimen near the S-IB interface. A close examination of figure 13a reveals that the gap at the 

S-IB interface extends outward from the centerline for about 3⁄4 of the specimen radius, and that a 

gap does not form at any radial location along the S-TB interface at this instant. The S-IB gap is 
36

clearly visible in figure 13b and is seen to be approximately uniform in size for much of its 

length.31 From figure 9, we conclude that the gap at the S-IB interface is about 37 m at the 

centerline; this is roughly the axial dimension of the deformed specimen elements. The axial 

stretch at this instant is about 0.79, corresponding to a nominal axial strain of 21%. Figures 

8 and 9 also confirm that there is no gap on the centerline at the S-TB interface at this time.

Figure 13. Deformed mesh plots in the vicinity of the specimen at t  236 s and an axial strain of

21%: (a) full specimen mesh; (b) enlargement of the mesh near the centerline, where the gap size 

is about 37 m.

Figures 14 and 15 plot the gap size at the S-IB interface vs. the deformed specimen radius r at 

selected times. Figure 16 gives the same information at the S-TB interface.

 

31 The latter observation is consistent with the results in figure 14.
37

Figure 14. Gap size vs. deformed specimen radius at the S-IB interface at selected times.

Figure 15. Gap size vs. deformed specimen radius at the S-IB interface at selected times.
38

Figure 16. Gap size vs. deformed specimen radius at the S-TB interface at selected times.

The selected times in figure 14 span Stage II and also include the instant t = 160 s, which is 

prior to any gap formation. Note that at each instant after 160 s, the largest gap occurs at a nonzero radial coordinate, that is, off the centerline. The largest gap size during this stage is nearly 

42 m, which occurs at 230 s and 2 mm from the centerline. From 170 to 230 s, the size of 

the gap increases with time and the gap extends along most of the specimen radius, that is, across 

most of the face of the specimen. Over this time interval, the nominal axial strain increases from 

4% to 19%. From 230 to 280 s, both the radial extent of the gap and the maximum size of the 

gap decrease with time. Over this time interval, the nominal axial strain increases from 19% to 

32%. We estimate that the gap extends over roughly half the radius of the specimen at t = 265 s 

and a nominal axial strain of 28%.

The selected times in figures 15 and 16 span Stages III–VI. From these figures, we conclude that 

the radial locations where the S-IB and S-TB gaps in Stage IV exceed a few microns are 

confined to the vicinity of the centerline. The largest gap sizes during this time interval occur at 

310 s and are nearly equal: 43 m at the S-IB interface and 41 m at the S-TB interface.

6.3 Pressure Contours

Figures 17–23 display pressure contours in the vicinity of the specimen at specific times between 

152 and 166.4 s. Also indicated at the top of each contour plot is the mean axial stretch z at 

the selected time. Each contour plot shows only a small portion of the length of the bars but does 

span the bars from the centerline (bottom) to the stress-free lateral surface (top). The incident 
39

bar is on the left and the transmission bar is on the right. Recall that the radius of the bars is 

slightly more than twice the initial radius of the specimen. Hence, the large, rectangular, white 

region between the bars in these contour plots is just empty space above the specimen. The 

specimen-bar interfaces are just barely distinguishable as the thin, light lines below this region. 

The length scale on these figures can be inferred from the thickness of the specimen, which is 

initially 1.45 mm and does not change substantially over the time interval considered here.

Since there is no load on the face of a bar where it overhangs the specimen, the stress on that part 

of the surface is zero, and since the bar faces remain very nearly planar, we may conclude that 

zz, rz, and z are essentially zero on the overhanging faces. But the stress tensor  need not be 

identically zero at points on the overhanging surface, so the pressure need not be zero there. 

Indeed, the radial and hoop stresses in the bars could be nonzero at some points on the 

overhanging surface, which would result in a nonzero pressure there; cf. equation 11.32 

Similarly, at those places and times at which the specimen and the bars are in contact, the 

pressure need not be continuous across the interface. Since the contact is frictionless, the only 

nonzero component of the traction vector at the specimen-bar interface is the axial stress zz, so it 

is this stress component that must be continuous across the interface when there is contact. 

Likewise, it is the axial stress (but not necessarily the pressure) that must be zero on the 

specimen and bar faces at points where a gap forms.

Recall that pressure is taken to be positive in compression (equation 11). The light green areas in 

the contour plots are regions of positive pressure, whereas the dark green areas are regions of 

negative pressure. The dark (and typically jagged) curves bounding these two green regions are 

the zero-pressure level curves. Clearly, the pressure field in the specimen and the bars is an 

extremely complex function of position and time. Reflection of waves from the interfaces and 

the free surfaces and oscillations induced in the specimen and bar by the passage of these waves 

make it difficult to sort out cause and effect here, particularly in view of the fact that the pressure 

is only part of the full stress state.

Figures 17–20 display pressure contours at specific times between 152 and 164.8 s, which is 

prior to the formation of a gap. Negative (tensile) pressures develop in the bars prior to 

developing in the specimen. Even at the early time of 152 s (figure 17), there are small regions 

of negative pressure in portions of the bars that overhang the specimen. These tensile regions are 

most likely generated by unloading waves from the stress-free portion of the bar faces. By 

163.4 s (figure 18), negative pressure has developed in the incident bar in a region bordering 

the centerline, about 6–8 specimen lengths back from the S-IB interface. From 164.2 to164.6 s 

(figures 19 and 20), the pressure in the incident bar oscillates (with axial distance from the 

specimen) between positive and negative values.

 

32 Likewise, since the low amplitude waves in the bars should induce little bulging of the stress-free lateral surface of the 

bars, rr, r, and zr are essentially zero on this surface, although the pressure need not be zero there.
40

Figure 17. Pressure contours in the vicinity of the specimen at selected times.
41

Figure 18. Pressure contours in the vicinity of the specimen at selected times.
42

Figure 19. Pressure contours in the vicinity of the specimen at selected times.
43

Figure 20. Pressure contours in the vicinity of the specimen at selected times.
44

Negative pressure does not develop in the specimen until 164.4 s (figure 20), although it is 

limited to very small regions. By 164.8 s, most of the specimen and the regions of the bars 

adjacent to it are in a state of negative pressure. Gap formation begins after 164.8 s (see section 

6.1.3). From 165.0 s onward, (figures 21–23), the regions of the bars adjacent to the specimen 

are primarily in a state of negative pressure. However, the sign of the pressure in the specimen 

oscillates (with time) from primarily negative (164.8–165.2 s) to primarily positive 

(165.4–165.8 s) and back again to primarily negative (166.0–166.4 s). In view of this 

oscillation in the sign of the pressure, it is possible that the axial stress also oscillates between 

tensile and compressive values at points in the interior of the specimen,33 although this cannot be 

confirmed from the pressure contours.

 

33 Recall that tensile axial stresses cannot be supported at the specimen and bar faces.
45

Figure 21. Pressure contours in the vicinity of the specimen at selected times.
46

Figure 22. Pressure contours in the vicinity of the specimen at selected times.
47

Figure 23. Pressure contours in the vicinity of the specimen at selected times.

Finally, note that at 152 and 162 s (figure 17), the pressure on the centerline in (and near) the 

specimen is around 1 MPa. This is consistent with the conclusion that the stress state on the 

centerline at these instants is nearly hydrostatic (see section 6.1.2). Indeed, from figures 

11 and 12, we see that zz is about 1 MPa at these times.

6.4 Effects of Pressure Bar Radius on Gap Formation

The results in section 6.3 suggest that the large overhang of the bars (particularly, the incident 

bar) may contribute to the formation of gaps in the specimen. In order to explore this hypothesis, 

we reduced the bar radius RB from 12.8 to 9.5 mm (see section 2.1). The histories of the gap size 

on the centerline for both cases are compared in figure 24. It is clear that the reduction in the bar 

radius substantially reduced the size of the gaps at both interfaces at later times.
48

Figure 24. Histories of the gap size on the centerline at the S-IB interface (top) and S-TB interface 

(bottom) for two different bar radii.
49

Somewhat surprisingly, however, figure 24 reveals that this reduction in the bar radius increased 

both the size and the duration of the initial gaps at both interfaces. Also, for the smaller diameter 

bar, the gaps at both interfaces form about 2 s earlier. These observations suggest a stronger 

interaction with the stress-free lateral surfaces of the bars, which are now closer to the specimen.

The effect of the bar radius on the formation of gaps for annular specimens is discussed in 

appendix C-2, but otherwise all results that follow are for the larger radius RB = 12.8 mm.

7. Simulations with the Moony-Rivlin Model and a 25-s Rise Time

All the simulations in this section are for the loading wave with a 25-s initial rise time and the 

compressible Mooney-Rivlin model for the specimen. For this case, the axial velocity histories 

at four locations on the centerline of the incident bar are given in figures 4 and 6. As observed in 

section 5.3, the velocity history at the S-IB interface for this case (figure 6) is substantially 

smoother than the corresponding velocity history for the 1-s initial rise time (figure 5). Recall

that the final rise time for the 1-s case is 24 s. For the 25-s initial rise time, the final rise 

time has increased to about 41 s. The purpose of this section is to examine the effects of this 

“pulse shaping” on the formation of gaps.

7.1 Results on the Centerline

All results in this subsection pertain to points on the centerline at the S-IB or S-TB interface. We 

begin with some remarks about the figures and a few observations, and then proceed to discuss 

the details of gap formation on the centerline.

7.1.1 General Remarks

Figure 25 compares the gap size histories on the centerline at the S-IB interface for the 1- and 

25-s initial rise times; figure 26 does the same at the S-TB interface. Clearly, the pulse shaping 

has substantially reduced the size and the duration of the gaps at both interfaces.
50

Figure 25. Histories of the gap size on the centerline at the S-IB interface for two loading waves with 

different initial rise times.
51

Figure 26. Histories of the gap size on the centerline at the S-TB interface for two loading waves with 

different initial rise times.

The gap size histories for the 25-s initial rise time are more easily discernable as the blue curves 

in figure 27 (S-IB interface) and figure 28 (S-TB interface). The largest gap at the S-IB interface 

is now only 2.45 m, while the largest gap at the S-TB interface is now only 1.27 m. The 

sensitivity of the gap size histories to the mesh size and the contact algorithm parameters is 

examined in appendix B. The results described there support the conclusion that the gap 

phenomenon is not a numerical artifact introduced by the contact algorithm or stemming from an 

insufficiently fine mesh.
52

Figure 27. Histories of the axial stress (negative in compression) and the gap size on the centerline at the S-IB 

interface. The stress is measured at the centroids of the specimen and incident bar elements 

adjacent to the centerline.
53

Figure 28. Histories of the axial stress (negative in compression) and the gap size on the centerline at the 

S-TB interface. The stress is measured at the centroids of the specimen and incident bar 

elements adjacent to the centerline.

Figure 27 also plots the histories of the axial stress zz in the specimen and the incident bar at 

points near the centerline and the specimen-bar interface. Figure 28 contains an analogous plot 

at the transmission bar interface. Recall that the stress is actually measured at the centroids of 

the specimen and bar elements adjacent to the centerline and the interface. Since the coordinates 

of these centroids do not coincide, slight differences in the values are possible. Also, recall that 

the axial stress is taken to be negative in compression. As expected, the axial stress is zero or 

nearly zero34 when a gap exists, even when the gap is submicron in size. Sharp stress spikes 

occur when the gap closes and the specimen “slaps” the bar.

We did not plot the mean axial stretch history for the case considered in this section. All stretch 

(and corresponding strain) values cited below were estimated from the axial stretch history for a 

 

34 Since the element centroids do not lie on the interface, the axial stress may have a small but nonzero value at the centroid 

even if it is zero on the face of the element due to the formation of a gap.
54

linear elastic specimen with a 25-s initial rise time;35 see figure 36 in section 8.2. The slope of 

the axial stretch curve in that figure is constant for t



190 s, and we estimate that the nominal 

strain rate for these times is about 2490/s.36 This is within 0.5% of the desired 2500/s strain rate.

A comparison of the velocity history at the S-IB interface in figure 6 with the gap and stress 

histories in figures 27 and 28 leads to the following additional observations. For convenience of 

discussion we have divided the time histories into four stages. The first two stages are analogous 

to those for the 1-s initial rise time (see sections 6.1.2 and 6.1.3).

7.1.2 The Initial Loading of the Specimen: Stage I (143–184 s)

This first stage covers the initial loading of the specimen prior to gap formation. It begins with 

the first significant increase in the axial velocity of the incident bar at the S-IB interface at about 

143 s (see figure 6) and ends with the onset of the first gap at 184 s. The axial stretch 

decreases from 1 to about 0.95 over this time 41-s interval; hence, the nominal axial strain 

increases from 0% to about 5%.

The first observable increase in axial stress occurs at about 149 s at the S-IB interface and 

150 s at the S-TB interface. The stress in the specimen increases to a sharp peak of 

1.6 and 1.5 MPa at the S-IB and S-TB interfaces, respectively, at about 164 s. Thereafter, 

the stress at both interfaces decreases to a broad local minimum centered at about 173 s, then 

increases again to a local maximum around 177–178 s, and finally decreases to zero at 184 s. 

The axial stress at this second (minor) peak is 1.1 and 1.0 MPa at the S-IB and S-TB 

interfaces, respectively. Smaller oscillations are superimposed on this general behavior.

The axial stretch at 164 s, the time of the first (major) peak in stress, is about 0.99. Since this is 

the same value of the stretch at the stress peak for the 1-s initial rise time (section 6.1.2), we 

conclude that at this instant, zz  rr 2.4 kPa and   rr. However, as noted above, 

zz  1.5 MPa. Hence, the differences between the principal stresses rr, , and zz are either 

zero or nearly three orders of magnitude smaller than the stresses themselves. It follows that all 

three principal stresses are nearly equal to each other and, hence, to –p. That is, the stress state is 

nearly hydrostatic.

While the preceding quantitative estimates apply only at the stress peak, similar calculations at 

neighboring times lead to the same qualitative conclusion: the stress state is nearly hydrostatic 

for most of this stage, as opposed to the desired state of uniaxial stress.

A close examination of the axial velocity history of the incident bar at the specimen interface 

(figure 6) reveals that the extrema in the axial stress in the specimen at 164, 173, and 177–178 s 

coincide (at least approximately) with inflection points in the velocity vs. time curve, that is, with 

 

35 Recall that the mean axial stretch is computed from equation 8. Because of the low impedance of the specimen for both

models, the displacements and hence the mean stretches for these two cases should be close.

36 We get essentially the same value using the approximation for the nominal strain rate in equation 10, together with 

vIB = 1.99 v0, which is the mean value for the velocity of the incident bar at later times, estimated from figure 6.
55

extrema in the axial acceleration vs. time curve. Similarly, as the acceleration decreases towards 

zero after 178 s, so does the axial stress. Thus, just as for the case of the 1-s initial rise time 

(section 6.1.2), it seems likely that during this stage, the value of the axial stress and the nearly 

hydrostatic stress state in the neighborhood of the stress peak are inertial effects associated with 

the accelerations imposed on the specimen by the incident bar.37 Furthermore, since the axial 

stresses at either face are nearly equal, we may conclude that the axial stress is approximately 

uniform through the thickness of the specimen (at least on the centerline), so these would appear 

to be primarily radial inertia effects.

The preceding conclusion is reinforced by comparing the peak stresses and accelerations for the 

two different rise times. The peak axial stress at the S-IB interface for the 1-s case (4.0 MPa) 

is about 21⁄2 times the peak stress (1.6 MPa) for the 25-s case. From the velocity histories in 

figures 5 and 6, we can obtain rough estimates for the axial accelerations of the incident bar at 

the time of these peak stresses. We find that the acceleration for the 1-s case is between 21⁄2 to 

3 times that for the 25-s case. Thus the peak axial stress in this early stage increases in (rough) 

proportion to the peak axial acceleration. Song et al. (6) reached a similar conclusion from their 

numerical simulations of SHPB tests on gel rubber specimens. This conclusion is also consistent 

with various theoretical analyses of inertial effects for soft materials (26–33).

7.1.3 Gap Opening and Closure (Stages II–IV)

Stage II (184–218 s): This stage begins with the opening of the first gap, which occurs at the 

S-TB interface, and ends with the first closing of the gap at the S-IB interface. During this 34-s 

time interval, the axial strain increases from 5% to about 14%. Note that the duration of this 

stage is about a third of the duration of Stage II for the case of a 1-s initial rise time, and gap 

formation occurs about 20 s later than it did for the 1-s case.

The axial stress at both interfaces drops to zero at about 184 s. The gap at the S-TB interface 

forms at this time, whereas the gap at the S-IB interface forms about 2 s later. Analogous to the 

case for the 1-s initial rise time, the contact algorithm allowed a 0.10-m interpenetration at the 

S-IB interface prior to gap formation. This interpenetration begins to decrease at about 184 s. 

From figure 6, we see that the peak in the axial velocity of the incident bar at the S-IB interface 

also occurs at 184 s. The axial acceleration of the bar is zero at this instant and then becomes 

negative for some time thereafter.

This stage ends with the first closing of the gap at the S-IB interface at about 218 s. The gap at 

the S-TB interface is also closed at this instant, although it has opened and closed several times 

previously. Nevertheless, it is open for most of this stage. The largest gap size at the S-IB 

interface (but not at the S-TB interface) is attained during this stage: about 2.45 m at 200 s.

 

37 In this regard, see the comments in the next to last paragraph in section 5.3.
56

Stage III (218–302 s): This stage begins with the first closing of the gap at the S-IB interface, 

which produces large spikes in the axial stress at both interfaces (figures 27 and 28) and a large 

spike in the axial velocity at the S-IB interface (figure 6). Subsequently, the gaps at both 

interfaces open and close many times, again with large spikes in stress and velocity 

corresponding to gap closure. This stage ends with the permanent closing of the gaps at both 

interfaces at about 302 s. During this 84 s time interval, the nominal axial strain increases 

from 14% to 34%. The largest gap at the S-TB interface occurs during this stage: 1.27 m at 

t  247 s.

Stage IV (302–340 s): The gaps at both interfaces close and remain closed during this final 

stage. The axial stress oscillates (without spikes) about a steadily increasing mean value. There 

is no analog of this stage for the case of a 1-s initial rise time, at least for times less than 340 s.

During Stage IV, the axial stretch z decreases monotonically from 0.66 to 0.57, and hence the 

nominal axial strain increases from 34% to 43%. From equation A-3 and the values of A1 and A2

in table 3, we find that for this range of z, zz would increase monotonically from about 

117 kPa to 175 kPa if the stress state were uniaxial (see also figure A-1). However, from figures 

27 and 28 we see that in fact zz oscillates between 0 and 1.1 MPa during this stage, with an 

average value around 0.5 MPa = 500 kPa. Hence, the specimen is not in a state of uniaxial stress 

during this stage. Indeed, using equation A-14 and arguing as in section 6.1.2, we conclude that 

for most Stage IV, the stress state is closer to hydrostatic than it is to uniaxial.

7.2 Gap Size vs. Radius at Selected Times

Figure 29 plots the gap size vs. deformed specimen radius at the S-IB interface at selected times. 

There is no gap at any radial location for times less than or equal to 185 s. Also, note the 

0.10 m interpenetration allowed by the contact algorithm at the S-IB interface prior to gap 

formation.
57

Figure 29. Gap size vs. deformed specimen radius at the S-IB interface at selected times.

At 190 s, the gap extends along most of the face of the specimen. Thereafter, the radial extent 

of the gap decreases with time. By 215 s, the gap extends only 3 mm out from the centerline, 

which is slightly less than half the radius. The maximum gap size increases with time up to 

200 s and then decreases with time until about 220 s, when the gap appears to close 

temporarily. Also, note that at each instant from 190–215 s, the maximum size of the gap 

occurs off the centerline. From 230 s onward, gap opening and closing exhibits a more 

complicated radial and temporal behavior. The results in figure 29 are consistent with the those 

on the centerline in figures 25 and 27.

7.3 Force History in the Transmission Bar

Figure 30 is somewhat similar to figure 28; in particular, the history of the gap size on the 

centerline at the S-TB interface is repeated in figure 30. But instead of plotting the axial stress in 

the transmission bar at a point very near the S-TB interface, figure 30 plots the total axial force 

Fz acting on the cross-section of the transmission bar at an axial location 2.52 mm from the 

specimen-bar interface. This is about 25 bar elements back from the interface. Also, note that Fz

has been taken positive in compression, whereas zz in figure 28 is positive in tension.38

 

38 When comparing figures 28 and 30, keep in mind that the colors of the stress-force curves and the gap size curves have 

been switched, as have the left-right locations of the vertical axes. 
58

Figure 30. History of the total force (positive in compression) at the stress gage location in the 

transmission bar (2.52 mm from the specimen interface) compared with the history of the 

gap size on the centerline at the S-TB interface.

On taking the difference in the sign conventions into account, we see (not unexpectedly) that the 

force history in figure 30 is qualitatively similar to the stress history in figure 28. In particular, 

the force is positive from 150 to 186 s, corresponding to the Stage I loading discussed in 

section 7.1.2, and the timing of the local maximum and minimum values of Fz andzz coincide 

during this stage. A gap opens and closes repeatedly at the centerline starting at 184 s; 

subsequently, the total force oscillates with repeated dips below zero (i.e., a tensile force). This 

continues until about 260 s, after which the force remains positive and oscillates about a 

steadily increasing mean value.

When comparing the force and gap histories in figure 30, it should be kept in mind that the total 

force is the cumulative effect of the stress at all radial locations in the transmission bar, whereas 

the gap has been measured only on the centerline. Nevertheless, from 184 to 302 s (Stages II 

and III in section 7.1.3), the local maximum and minimum values of the total force coincide, for 

the most part, with the closure and opening, respectively, of a gap on the centerline. In 

particular, the times at which the total force becomes tensile coincide with gaps, for the most 

part. It seems reasonable to conclude that when a gap exists on the centerline at the S-TB 
59

interface, it extends over a significant portion of the specimen face, although we do not have any 

plots (analogous to figure 29) to verify this.

The force in the transmission bar in figure 30 is measured at an axial location 2.52 mm from the 

specimen-bar interface. This location was selected because it is typical for the placement of a 

quartz stress gage in the bar in SHPB tests (Casem et al. [25]). Hence, we refer to this axial 

location as the gage location, even though a stress gage was not included in the simulation. The 

signal from a quartz gage is proportional to the total force applied to it, regardless of whether or 

not the axial stress in the transmission bar is radially uniform (25). With appropriate 

calibrations, the force measured by a quartz gage should also coincide with the force inferred 

from semi-conductor strain gage measurements midway down the bar; see Moy et al. (5) for the 

case of a ballistic gelatin specimen.

In view of the proximity of the gage location to the S-TB interface, we expect that the force on 

the transmission bar at the specimen interface may be approximated by the force at the gage 

location, with some qualifications (see below). Hence, keeping in mind the sign conventions and 

letting A0 denote the original cross-sectional area of the specimen, we see that Fz /A0 provides 

an estimate for the average or mean value of the axial component of the nominal stress in the 

specimen.39 Similarly, 

zz

, the mean value of the true axial stress zz in the specimen, is 

estimated by dividing Fz by the deformed cross-sectional area of the specimen, which in turn is 

approximated by A0 /z ,

40 so that 

  zz z Fz /A0. These estimates for the mean value of the 

nominal and true axial stress in the specimen are consistent with those used in reporting data in 

the SHPB literature.

Let us apply this estimate for 

zz

to the first peak in Fz at 164-165 s, which is prior to gap 

formation. The force at this instant is 25 N, z is about 0.99 (see section 7.1.2), and the 

undeformed specimen radius is 6.35 mm, which gives 

  zz 0.20 MPa = 200 kPa. This is 

almost two orders of magnitude above the value of 2.4 kPa that would correspond to a state of 

uniaxial stress at this axial stretch. On the other hand, the axial stress zz on the centerline at the 

S-TB interface at this instant is 1.5 MPa (section 7.1.2), which is 7–8 times the estimated mean 

stress. Hence, if this estimate for the mean stress is valid (at least approximately), then the axial 

stress must be radially non-uniform, with higher values near the centerline. This qualitative 

behavior is consistent with theoretical analyses of radial inertia effects for soft materials (7, 26, 

27, 32, 33) and with the results of other numerical simulations of SHPB tests (6, 34).

For early times, the force history in figure 30 is qualitatively similar to the mean nominal stress 

histories on solid specimens in figure 3 of Moy et al. (5) and figure 5 of Song et al. (6). The 

latter two figures report experimental SHPB data on a ballistic gelatin specimen and a gel rubber 

specimen, respectively. All three figures show a large “inertial spike” at small strains and early 

 

39 This is also referred to as the engineering stress. It is the zz-component of the 1st Piola-Kirchhoff stress tensor.

40 This approximate relation is exact when the axial stretch in the specimen is uniform and there is no volume change. 
60

times (165 s in figure 30). In Moy et al. (5) and Song et al. (6) it is demonstrated that this 

inertial spike can be substantially reduced by hollowing out the specimen, the rational being that 

radial inertia contributions to the axial stress are largest in the center of the specimen (cf. the 

preceding paragraph). Unlike figure 30, the experimental data in Moy et al. (5) does not exhibit 

the large oscillations or any evidence of gap formation after the inertial spike. This is not 

surprising in view of the fact that the pulse shaping in those tests resulted in a much smoother 

loading wave than the one used in our simulation. On the other hand, the entire force history in 

figure 30 is qualitatively similar to nominal stress history for the solid specimen in figure 5 of 

Song et al. (6), for which no pulse shaping was used. These similarities include the oscillations 

and the dips to negative (i.e., tensile) forces following the inertial spike. Thus the data in Song et 

al. (6) may be exhibiting indirect evidence of gap formation at the S-TB interface, although the 

possibility of gaps was not discussed in that paper.41

Based on the discussion above and on the previous discussion of the histories of the stress, 

velocity, and gap size on the centerline for both the 25- and 1-s initial rise times, the following 

conclusions seem plausible. Not only is the inertial spike a radial inertia effect, but so also is the 

reduction from this large compressive stress to zero axial stress, even though the specimen is 

subject to compressive axial strain. This reduction to zero axial stress is driven by the rapid drop 

in the axial acceleration of the incident bar to (and below) zero. Once the axial stress at the 

interface drops to zero,42 a gap can form there. For cases with more gradual deceleration via 

better pulse shaping, as in Moy et al. (5), there is a less severe drop in stress that never reaches 

zero before the stress begins to increase again.

Finally, it should be kept in mind that when a gap exists over a substantial portion of the face of 

the specimen, the contact area between the specimen and the bar is reduced, so the local value of 

the axial stress in the specimen is magnified at those points in contact with the bar, whereas it is 

zero at points on the specimen face where a gap has formed. Consequently, in this case the mean 

axial stress computed as above does not provide a meaningful estimate of the highly non-uniform 

stress state in the specimen. Furthermore, differences between the force at the gage location and 

the force at the S-TB interface could arise due to axial acceleration of the transmission bar near 

the interface (25). This acceleration is often regarded as negligible (based the assumption that 

the transmitted pulses are weak), but it must be significant at least part of the time in the present 

 

41 The stress history at the S-TB interface in Song et al. (6) was inferred from strain gages mounted on the transmission bar. 

Similarly, the stress history at the S-IB interface in a SHPB test can be inferred from strain gages mounted on the incident bar. 

Large stress oscillations and tensile stresses are also seen in some incident bar strain gage data on soft materials; cf. figure 6 

(lead) in Gray (1) and figure 2 (estane-based polymer binder) in Gray and Blumenthal (3). This could be indicative of gap 

formation at the S-IB interface. However, the estimate for the stress at the S-IB interface based on 2-wave calculations from 

strain gage data is subject to larger uncertainty than the estimate for the stress at the S-TB interface; cf. Gray (1), Gray and 

Blumenthal (3), and Chen et al. (35). Quartz gages were used to measure force histories in the incident and transmission bars in 

SHPB tests on RTV 630 silicone rubber in Chen et al. (35). Large stress oscillations and tensile stresses are seen in the incident 

bar data here as well (figures 6 and 8); and some of the transmission bar data shows less frequent and smaller amplitude dips to 

tensile stress. These results could be indicative of gap formation at both interfaces. However, the authors point out that the 

specimen may have failed in some of these tests, and this possibility complicates the interpretation of the data.

42 Note that tensile axial stresses cannot be transmitted across the specimen-bar interface. 
61

case. Indeed, the tensile forces at the gage location at various times between 186 and 260 s 

cannot possibly represent the force at the S-TB interface, since tensile forces cannot be

transmitted across the interface: the force on exerted on the transmission bar by the specimen 

(and hence, on the specimen by the transmission bar) must be compressive (positive) or zero.

8. Simulations with a Linear Elastic Model for the Specimen

For all of the results presented up to this point, the compressible Mooney-Rivlin model was used 

for the specimen. Recall that this is a nonlinear elastic constitutive model that was calibrated to 

yield a nearly incompressible material (see section 3 and appendix A). The question naturally 

arises as to whether the nonlinearity in the material response and the near incompressibility of 

the material are necessary for the formation of gaps or for the persistence of the gaps once they 

form. The simulations in this section address these issues by replacing the Mooney-Rivlin model 

for the specimen with an isotropic, linear elastic model and by varying the Poisson’s ratio for 

that model.43

We used LS-DYNA’s Orthotropic Elastic model (Material Model 2 in LS-DYNA [8] and 

Hallquist [9]) for the specimen, with the material constants chosen in such a way as to yield an 

isotropic, linear elastic model. For all of the simulations in this section, the Young’s modulus E

was fixed at 0.43 MPa, which is slightly less than twice the initial value of the Young’s modulus 

used for the Mooney-Rivlin model (see table 1). The density of the specimen remained the same 

(1 g/cm3). Six different values of the Poisson’s ratio  were used in the simulations, ranging 

from 0.49999 to 0.48. The largest value,   0.49999, is the same value used for the 

compressible Mooney-Rivlin model. The six values for Poisson’s ratio as well as the 

corresponding values of the other elastic constants determined from E and  are listed in table 1. 

Note that as Poisson’s ratio decreases from 0.49999 to 0.48, the bulk modulus K and the 

longitudinal modulus L (which are nearly identical) decrease by a factor of 2000, whereas the 

changes in the shear modulus G are insignificant; the ratio of the shear to bulk modulus, G/K, 

which is a measure of incompressibility (see section 2.3), increases from 0.0002 to 0.04.

8.1 The 1-s Initial Rise Time

Figures 31–35 present results for the loading wave with a 1-s initial rise time. Figure 31 

compares the histories of the gap size on the centerline at the two interfaces for the case 

  0.49999. Also plotted in that figure (right axis) is the history of the axial stretch z; see 

equation 8. Clearly, nonlinearity in the constitutive model is not required for the formation of a 

gap or for the persistence of the gap out to large strains. The gap forms at about 165 s (3% 

nominal strain). The gap at the S-IB interface remains open for a substantial time interval, 

 

43 The elastic constants used in these simulations are not necessarily those appropriate for the small strain response of 

ballistic gelatin.
62

eventually closing at about 263 s (27% nominal strain). The S-TB gap is also closed at this

time, but has opened and closed several times prior to that, and the size of the gap at the S-TB 

interface during this first stage is substantially less than at the S-IB interface. All of these 

features are qualitatively similar to those observed for the nonlinear model (see section 6.1).

Figure 31. A comparison of the histories of the gap sizes on the centerline at the S-IB and S-TB interfaces 

for a 1-s rise time and   0.49999 (linear elastic model). Also shown is the history of mean 

axial stretch in the specimen (blue).
63

Figure 32. A comparison of the histories of the gap sizes on the centerline at the S-IB interface for 

a 1-s rise time and selected values of Poisson’s ratio (linear elastic model).

Figure 32 compares the histories of the gap size on the centerline at the S-IB interface for the six 

values of Poisson’s ratio. Figure 33 does the same at the S-TB interface. Gaps form at both 

interfaces for 



 that is, for G/K



0.02. But they do not form at either interface for the 

lowest value of Poisson’s ration,   0.48 that is, for G/K = 0.04. Thus it appears that the 

specimen must be nearly incompressible in order for gaps to form.
64

Figure 33. A comparison of the histories of the gap sizes on the centerline at the S-TB interface for 

a 1-s rise time and selected values of Poisson’s ratio (linear elastic model).

The three values of Poisson’s ratio that resulted in the largest gap sizes at both interfaces were 

0.499, 0.4995, and 0.4999. Note that the corresponding gap histories in figures 32 and 33 do not 

extend all the way to 360 s. We terminated these histories prior to that time because the 

specimen elements in direct contact with the bar exhibited hourglassing. The gap size curves in 

the figures were terminated well before this occurred.

It is interesting to note that the maximum value (over the duration of the simulation) of the size 

of the gap at the centerline does not vary monotonically with Poisson’s ratio. Figure 34 plots this 

maximum gap size at each interface as a function of Poisson’s ratio . At the S-IB interface, the 

maximum gap size increases with  up to a maximum of 140 m for  and decreases 

with thereafter. At the S-TB interface, the maximum gap size increases with  up to a 

maximum of 180 m for  and decreases with thereafter. For the 

maximum gap size at both interfaces is essentially the same (about 140 m). For smaller values 

of the largest gap occurs at the S-TB interface; for larger values of the largest gap occurs at 

the S-IB interface.
65

Figure 34. The maximum gap size at the centerline for each interface as a function of Poisson’s ratio 

(linear elastic model, 1-s rise time).

The gap of 180 m on the centerline at the S-TB interface for the case is the largest 

gap size observed in all of the simulations in this report. This occurred at 276 s, which is about 

110 s after the loading wave arrived at the specimen. Figure 35 plots the deformed mesh in the 

vicinity of the specimen at this instant for this value of Poisson’s ratio. The centerline lies along 

the bottom of the plot, so the figure shows the entire deformed specimen mesh in the 2-D 

axisymmetric simulation. It is clear that at this instant the gaps at either interface are largest on 

the centerline. Also, the gap extends further outward radially at the S-TB interface. Along most 

of the radial extent of the gap at the S-TB interface and along half of the radial extent of the gap 

at the S-IB interface, the gap size exceeds the 100-m length of bar elements and hence is 

several times larger than the 50-m initial length of the undeformed specimen elements, which in 

turn is much larger than the deformed element length at this instant. This observation supports 

the assessment that the gap phenomenon cannot be attributed to numerical artifact associated 
66

with an insufficiently fine mesh. Finally, note that according to figure 31, the estimated axial 

stretch at this instant is about 0.69; the corresponding nominal strain is 31%. However, these 

estimates are based on equations 8 and 9, that is, on the distance between the two bars. It is clear 

from figure 35 that for radial locations near the centerline, the axial strain in the specimen is 

substantially higher than this estimated value.

Figure 35. Deformed mesh in the vicinity of the specimen at t  276 s for Poisson’s ratio 

 (linear elastic model, 1-s rise time).
67

8.2 The 25-s Initial Rise Time

Figures 36–40 present results for the loading wave with a 25-s initial rise time. Figure 36 

compares the histories of the gap size on the centerline at the two interfaces for the case 

  0.49999. Also plotted in that figure (right axis) is the history of the axial stretch z. The gap 

forms at about 185 s (5% nominal strain). The pulse shaping has delayed the onset of gap 

formation, has substantially reduced the gap size at both interfaces (the largest value is now 

about 2.6 m), and has also substantially reduced the duration over which the gap first remains 

open at the S-IB interface. These features are similar to those observed for the nonlinear model 

(see section 7).

Figure 36. A comparison of the histories of the gap sizes on the centerline at the S-IB and S-TB 

interfaces for a 25-s rise time and   0.49999 (linear elastic model). Also shown is 

the history of mean axial stretch in the specimen (blue).
68

Figure 37. A comparison of the histories of the gap sizes on the centerline at the S-IB interface 

for a 25-s rise time and selected values of Poisson’s ratio (linear elastic model).

Figures 37 compares the histories of the gap size on the centerline at the S-IB interface for the 

six values of Poisson’s ratio. Figure 38 does the same at the S-TB interface. Just as for the 1-s 

initial rise time, gaps form at both interfaces for 



 but they do not form at either interface 

for the lowest value of Poisson’s ration,   0.48. For simulations that exhibited hourglassing, 

the gap history curves were terminated well before hourglassing occurred.
69

Figure 38. A comparison of the histories of the gap sizes on the centerline at the S-TB interface 

for a 25-s rise time and selected values of Poisson’s ratio (linear elastic model).

Just as for the 1-s initial rise time, the maximum value (over the duration of the simulation) of 

the size of the gap at the centerline does not vary monotonically with Poisson’s ratio. Figure 39

plots this maximum gap size at each interface as a function of Poisson’s ratio . At both 

interfaces the maximum gap size increases with  up to , and decreases with 

 thereafter. The maximum gap sizes for each value of Poisson’s ratio are roughly the same at 

both interfaces, with the largest differences occurring for .
70

Figure 39. The maximum gap size at the centerline for each interface as a function of Poisson’s 

ratio (linear elastic model, 25-s rise time).

Figure 40 plots the deformed mesh in the vicinity of the specimen at 286 s for the case 

. For this value of and this instant, the gap size at the S-TB interface achieves its 

largest value of 62 m, and the gap at the S-IB interface is 71 m, which is nearly equal to its 

largest value of 72 m (see figures 38 and 39). From figure 36, we see that the estimated axial 

stretch at this time is about 0.70, which implies a nominal strain of 30%. The centerline lies 

along the bottom of the plot, so the figure shows the entire (deformed) specimen mesh in the 2-D 

axisymmetric simulation. Both gaps extend to about half the radius of the specimen. Along 

most of the radial extent of the gap at either interface, the size of the gap exceeds the initial

length of the specimen elements and is about twice the deformed length of the elements.
71

Figure 40. Deformed mesh in the vicinity of the specimen at t  286 s for Poisson’s ratio 

 (linear elastic model, 25-s rise time).
72

9. Discussion and Concluding Remarks

The SHPB test is designed to impose a state of compressive uniaxial stress on the specimen, at 

least after an initial ringing-up period. Even during the ring-up, one expects that the specimen is 

in a (generally non-uniform) state of compression, since the initial compressive longitudinal 

wave in the specimen reflects from the higher impedance pressure bars as a compressive wave. 

Thus the implicit assumption that the specimen and pressure bars remain in contact would appear 

to be reasonable.

Nevertheless, in this report we have presented conclusive numerical evidence that gaps may open 

between the specimen and the bars under certain conditions. These gaps formed at small strains 

but existed out to large strains, closing and re-opening multiple times. In some cases the gaps 

persisted for over 100 s and extended over much of a face of the specimen for most of that 

time. The results in the appendix support the conclusion that this gap phenomenon is not an

artifact of an insufficiently fine mesh (appendix B-1) or of the contact algorithm (appendix B-2) 

or of the particular code used for the simulations (appendix C-3). If a gap does exist over a 

substantial portion of the specimen face in a SHPB test, then that test is invalid for the purposes 

of material property characterization since the estimates for the axial stress, strain, and strain rate 

are no longer representative of the non-uniform conditions in the specimen.

The results of this study indicate that for soft specimens, gaps are most likely to form under a 

combination of two conditions, one pertaining to the specimen properties and the other to the 

loading conditions. For a given loading condition, gaps seem more likely to form in specimens 

which are nearly incompressible. This condition was satisfied in the two simulations with the 

nonlinear elastic (Mooney-Rivlin) model for the specimen, and gaps formed in both cases 

(sections 6 and 7). When a linear elastic model was used for the specimen, gaps formed for five 

of the six sets of elastic moduli considered (section 8). The case where no gaps formed 

corresponded to the lowest Poisson’s ratio (0.48) and hence to the highest ratio of shear to bulk 

modulus (0.04).

The loading condition that seems to promote gap formation is axial deceleration of the specimen. 

As the time rate of change of the velocity imposed by the incident bar on the specimen decreases 

from peak positive values (acceleration) through zero to negatives values (deceleration), the axial 

stress at the specimen-bar interfaces drops to zero, allowing the specimen to separate from the 

bars. It appears that the driving mechanism behind this inertial effect is the corresponding radial 

deceleration induced by the axial deceleration. This conclusion is not inconsistent with existing 

literature on radial inertial effects for soft specimens (26–33). Furthermore, we found that 

constraining the lateral surface of the specimen (to eliminate the radial motion) suppresses gap 

formation, at least in direct impact tests (appendix C-1.1).
73

On the other hand, the simple qualitative conclusions above belie the extremely complex stress 

states in the specimen and pressure bars (particularly the incident bar) at times in the vicinity of 

gap formation, as illustrated by the pressure contour plots in section 6.3.44 Starting about 1.6 s 

prior to the opening of gaps at both interfaces and at least 1 s prior to the onset of negative 

pressure in the specimen, the pressure in the incident bar oscillates (with axial distance from the 

specimen) between positive and negative values, eventually becoming negative in the entire 

vicinity of the specimen immediately before the gaps form (section 6.3).

To the best of our knowledge, the gap phenomenon documented here has not been reported in 

either the experimental or the computational literature. However, certain features of previously 

reported experimental data appear to be consistent with the opening and closing of gaps, namely, 

axial stress histories (inferred from strain or stress gage measurements in the bars) which 

oscillate between positive and negative values (see section 7.3). Since our initial concern was to 

determine whether or not these gaps were a numerical artifact, we focused on examining a few 

cases in detail rather than simulating a broad range of loading conditions and constitutive models 

for the specimen. Consequently, our study is by no means exhaustive. We conclude by listing 

some of the limitations of this work, which in turn suggest areas for further study:

1) We certainly expect that gaps will no longer form if the strain rate in the specimen is 

sufficiently small, but only one axial strain rate was considered here, a nominal rate of 

2500/s (after the initial ring-up). The imposed strain rate can be reduced either by 

decreasing the plateau velocity imposed at the far end of the incident bar or by increasing 

the thickness of the specimen.

2) We considered only a linear ramp to a plateau velocity in the main body of the report (and 

more severe loadings in the appendices C-1 and C-3). It would be of interest to consider 

smoother loading histories at the far end of the incident bar, since experimental evidence 

indicates that gaps do not form with better pulse shaping.

3) We considered only one specimen diameter (12.7 mm). Theories on radial inertial effects 

in SHPB tests predict that the inertial stresses are proportional to the square of the 

specimen diameter. Thus if radial inertia is a mechanism for gap formation, we would 

expect that gaps would no longer form if the specimen diameter is sufficiently small.

4) All the simulations in this report involved specimens that were “soft” relative to the bars 

(particularly in shear), so we cannot draw any conclusions about gap formation in stiffer 

specimens.

 

44 Recall that these pressure contours were for the loading wave with a 1-s initial rise time.
74

5) Only elastic (linear or nonlinear) constitutive models were considered for the specimen, 

since absence of strain rate effects in the model allowed for easier interpretation of the 

results. It is not clear how the incorporation of viscoelasticity in the model would affect 

the formation of gaps, although we have no reason to believe that this would suppress gap 

formation altogether.

6) For the nonlinear elastic specimen model (the compressible Mooney-Rivlin model), we 

considered only one set of material parameters. It would be instructive to fix the bulk 

modulus and increase the shear modulus until gap formation is suppressed. Note that in the 

simulations with a linear elastic specimen, we essentially fixed the shear modulus and 

decreased the bulk modulus until gap formation was suppressed.

7) As an approximation to a well lubricated specimen, the specimen-bar interfaces were 

treated as frictionless in the simulations. It is not clear what effect (if any) a small amount 

of friction would have on the formation of gaps.
75

10. References

1. Gray III, G. T. Classic Split-Hopkinson Pressure Bar Testing, ASM Handbook Vol. 8, 

Mechanical Testing and Evaluation; American Society for Metals: Materials Park, Ohio, 

462–476, 2000.

2. Gama, B. G.; Lopatnikov, S. L.; Gillespie Jr., J. W. Hopkinson bar experimental 

technique: a critical review. Applied Mechanics Reviews 2004, 57, 223–250.

3. Gray III, G. T.; Blumenthal, W. R. Split-Hopkinson Pressure Bar Testing of Soft 

Materials, ASM Handbook Vol. 8, Mechanical Testing and Evaluation; American Society 

for Metals: Materials Park, Ohio, 488–496, 2000.

4. Song, B.; Chen, W. Split Hopkinson pressure bar techniques for characterizing soft 

materials. Latin American Journal of Solids and Structures 2005, 2, 113–152.

5. Moy, P.; Weerasooriya, T.; Juliano, T. F.; VanLandingham, M. R.; Chen, W. Dynamic 

Response of an Alternative Tissue Simulant, Physically Associating Gels (PAG); ARL-RP136; U.S. Army Research Laboratory: Aberdeen Proving Ground, MD, September 2006.

6. Song, B.; Ge, Y.; Chen, W. W.; Weerasooriya, T. Radial Inertia Effects in Kolsky Bar 

Testing of Extra-soft Specimens. Experimental Mechanics 2007, 47, 659–670.

7. Scheidler, M. Inertial Effects in the Characterization of Incompressible Soft Materials in 

Hopkinson Bar Tests. presented at the 44th Annual Meeting of the Society of Engineering 

Science, Texas A&M University, October 21–24, 2007.

8. Livermore Software Technology Corporation, LS-DYNA Keyword User’s Manual, Volume 

II, Version 971, May 2007.

9. Hallquist, J. O. LS-DYNA Theory Manual, March 2006.

10. Aihaiti, M.; Hemley, R. J. Equation of state of ballistic gelatin; Report No. 50533-EG.1;

U.S. Army Research Office, June 2008.

11. Raftenberg, M. N.; Scheidler, M. Gap Formation in Simulations of SHPB Tests on Elastic 

Materials Soft in Shear, Shock Compression of Condensed Matter–2009. AIP Conference 

Proceedings, American Institute of Physics; Melville, New York, 2009; Vol. 1195: pp 

715–718.

12. Eringen, A. C.; Suhubi, E. S. Elastodynamics vol. II: Linear Theory; Academic Press: 

New York, 1975.
76

13. Mooney, M. A Theory of Large Elastic Deformation. Journal of Applied Physics 1940, 

11, 582–592.

14. Rivlin, R. S. Torsion of a Rubber Cylinder. Journal of Applied Physics 1947, 18, 444–

449.

15. Treloar, L.R.G. The Physics of Rubber Elasticity, 3rd Ed.; Clarendon Press: Oxford, 1975.

16. Batra, R. C. Elements of Continuum Mechanics, AIAA, Reston, VA, 2005.

17. Ogden, R. Non-Linear Elastic Deformations; Dover Publications: Mineola, NY, 1997.

18. Truesdell, C.; Noll, W. The Non-Linear Field Theories of Mechanics, 2nd Ed.; Springer:

New York, 1992.

19. Holzapfel, G. A. Nonlinear Solid Mechanics; John Wiley & Sons: New York, 2000.

20. Winter, J.; Shifler, D. The Material Properties of Gelatin Gels, BRL Contractor Report 

No. 217, March, 1975.

21. Rapacki, E. Private communication, ARL, 2009.

22. Juliano, T. F.; Forster, A. M.; Drzal, P. L.; Weerasooriya, T.; Moy, P.; VanLandingham, M. 

R. Multiscale Mechanical Characterization of Biomimetic Physically Associating Gels. J. 

Mater. Res. 2006, 21, 2084–2092.

23. Steinberg, D. J. Equation of State and Strength Properties of Selected Materials; UCRLMA-106439; Lawrence Livermore National Laboratory Report, Feburary 1996.

24. Graff, K. F. Wave Motion in Elastic Solids; Dover Publications: New York, 1991.

25. Casem, D.; Weerasooriya, T.; Moy, P. Inertial Effects of Quartz Force Transducers 

Embedded in a Split Hopkinson Pressure Bar. Experimental Mechanics 2005, 45, 368–

376.

26. Forrestal, M. J.; Wright, T. W.; Chen, W. The Effect of Radial Inertia on Brittle Samples 

During the Split Hopkinson Pressure Bar Test. International Journal of Impact 

Engineering 2007, 34, 405–411.

27. Dharan, C.K.H.; Hauser, F. E. Determination of Stress-Strain Characteristics at Very High 

Strain Rates. Experimental Mechanics 1970, 10, 370–376.

28. Gorham, D.A.P.; Pope, P. H.; Cox, O. Sources of Error in Very High Strain Rate 

Compression Tests. Institute of Physics Conference Series No. 70, 1984, 151–158.

29. Walley, S. M.; Field, J. E.; Pope, P. H; Safford, N. A. The Rapid Deformation Behaviour 

of Various Polymers. J. Phys. III France 1, 1991, 1889-1925.
77

30. Maliowski, J. Z.; Klepaczko, J. R.; Kowalewski, Z. L. Miniaturized Compression Test at 

Very High Strain Rates by Direct Impact. Experimental Mechanics 2007, 47, 451–463.

31. Jia, D.; Ramesh, K. T. A Rigorous Assessment of the Benefits of Miniaturization in the 

Kolsky Bar System. Experimental Mechanics 2004, 44, 445–454.

32. Scheidler, M.; Kraft, R. Inertial Effects in Compression Hopkinson Bar Tests on Soft 

Materials. to appear in the Proceedings of the 1st Annual ARL Ballistic Technologies 

Workshop, 2010.

33. Warren, T. L.; Forrestal, M. J. Comments on the Effect of Radial Inertia in the Kolsky Bar 

Test for an Incompressible Material, Experimental Mechanics, published online 29 

December 2009, DOI 10.1007/s11340-009-9322-x.

34. Chung, D.-T. The Effect of Specimen Shape on Dynamic Flow Stress, Shock Compression 

of Condensed Matter–1995, AIP Conference Proceedings, American Institute of Physics: 

Woodbury, New York, 1996; Vol. 370: pp 483–486.

35. Chen, W.; Lu, F.; Frew, D.; Forrestal, M. Dynamic Compression Testing of Soft 

Materials. Journal of Applied Mechanics 2002, 69, 214–223.

36. Scheidler, M.; Raftenberg, M.; Love, B.; Kraft, R. Inertial Effects in Numerical 

Simulations of Hopkinson Bar Tests on Nearly Incompressible Soft Materials, presented at

The Joint ASCE-ASME-SES Conference on Mechanics and Materials, Virginia Tech, June 

24–28, 2009.

37. Love, B. Unpublished results, ARL, 2008.
78

INTENTIONALLY LEFT BLANK.
79

Appendix A. The Mooney-Rivlin Constitutive Model

A-1 Theoretical Background

This subsection contains a brief discussion of the concepts and notation needed to describe the 

two nonlinear elastic models considered below; for additional background the reader may refer 

to the books (16–19).

Let 1, 2, 3 denote the principal stretches, that is, the ratios of the deformed to undeformed 

length along the principal axes of strain. Each stretch i is unity in the undeformed state, with 

i  1 in tension, and 0  i  1 in compression. The principal nominal (engineering) strains are 

iwhen taken positive in tension, or 1i when taken positive in compression. The left 

Cauchy-Green deformation tensor B is defined in terms of the deformation gradient F by 

B = FFT, where the superscript T denotes the transpose. The principal axes of B are the 

principal axes of strain in the deformed state. The principal values (or eigenvalues) of B are 

1

2

, 2

2

, 3

2

, the squares of the principal stretches. The three principal invariants of B are

222

1 1 2 3 I     tr B  , (A-1)

2 2 2 2 2 2 2 2

2 1 2 2 3 3 1

1

( tr ) tr

2

I      [ ] B B       , (A-2)

2 2 2 2

3 1 2 3 I J    det B  . (A-3)

Here “tr” denotes the trace, “det” denotes the determinant, and

J  det F  123 (A-4)

is the Jacobian of the deformation, that is, the local ratio of deformed to undeformed volume. 

Thus, J  1 at points where there is no volume change. A useful alternative expression for I2 is

1

2 3 3 222

1 2 3

111 I I I tr





 

       B . (A-5)

The principal stresses (i.e., the eigenvalues) of the Cauchy stress tensor  are denoted by 

, 2,  3, and the pressure p is given by

 1 2 3 

1 1

tr

3 3

p           . (A-6)

For any (not necessarily linear) isotropic elastic material, the principal axes of stress coincide 

with the principal axes of B.
80

A-2 The Mooney-Rivlin Model for Incompressible Elastic Materials

The classical Mooney-Rivlin model is a nonlinear elastic constitutive model for isotropic, 

incompressible solids; cf. the original papers by Mooney (13) and Rivlin (14), as well as the 

books (15–19). Since the material is assumed to be incompressible, there is no change in 

volume, so equations A-3 and A-4 reduce to

I3  J  123  1. (A-7)

The strain energy function W for the Mooney-Rivlin model has the simple form

W A I A I     1 1 2 2  3 3   , (A-8)

The coefficients A1 > 0 and A2 ≥ 0 are constants with the units of stress (i.e., elastic moduli); the 

inequalities are required for physically realistic behavior. The Cauchy stress tensor 

corresponding to this strain energy function is

1

p A A 2 2 1 2

      I B B , (A-9)

where I denotes the identity tensor. The shear modulus G in the small strain, linear elastic 

approximation to the constitutive relation A-9 is given by G  2(A1  A2). The special case of 

equation A-9 with A2 = 0 (and hence 2A1 = G is called a neo-Hookean material.

The coefficient 

p

in equation A-9 denotes an indeterminate scalar related to the pressure. The 

indeterminacy of 

p

is a consequence of the incompressibility constraint: 

p

is not determined 

by B, but the value of 

p

at each place and time must be such that the momentum balance 

equations and the boundary conditions are satisfied. Although 

p

in equation A-9 is often 

referred to as a “pressure”, it is generally not equal to the pressure p as defined in equation A-6, 

since the other terms in equation A-9 contribute to the trace of .

The relation A-9 implies that the principal stresses are given by

2

1 2 2

1

2 2 i i

i

p A A 



     (i  1,2,3). (A-10)

The difference of any two principal stresses is independent of 

p

:

 

2 2

1 2 2 2

1 1 2 2 i j i j

i j

A A  

 

 

            

. (A-11)

Uniaxial Stress: The general relation A-11 may be used to obtain the stress-stretch relation in a 

uniaxial stress test. Consider a cylindrical coordinate system with the direction of applied stress 

parallel to the z-axis. Then zz is the only non-zero stress component. If zz  0, this is a simple 

tension test; if zz  0, it is a simple compression test. Since the material is isotropic, the 

principal axes of strain lie along the coordinate axes in this case and the strains (and hence 

stretches) in the radial and hoop directions are equal. Thus, in the relations above, we may set 
81

1  r    2 and 3  z, as well as 1  rr    2  0 and 3  zz. Then by equations 

A-6 and A-7, we see that p  zz/3 and

1

r

z

 



 

. (A-12)

Then on setting i  3 and j  1 in equation A-11 and using the relations above, we obtain the 

following expression for the axial stress zz in terms of the axial stretch z:

2

1 2 2

1 1 2 2 3 ( 1)

   

               z z z

z z

zz A A G   

 

. (A-13)

The expression on the right is the linearly elastic approximation, which is valid for small strains 

only, that is, z close to 1. For an incompressible material, 3G  E, where E is the Young’s 

modulus, so we recover the well-known result that in the linear elastic approximation, the axial 

stress is the Young’s modulus times the axial strain.

Stress Difference on the Centerline: Consider an axisymmetric (but possibly non-uniform) 

deformation of a solid cylindrical specimen. As discussed at the end of section 2.6, at points on 

the centerline (i.e., the z-axis), the radial, hoop and axial directions are principal axes of stress 

and strain, and the relations 14 in section 2.6 hold exactly. Thus, the principal stress difference 

zz  rr is given by equation A-11 with the indices i and j set to z and r, respectively; and since 

J  1, equation A-12 holds. Therefore,

2

1 2 2

1 1 2 2 3 ( 1)

z z z

z z

zz rr A A G   

 

   

                

, (A-14)

where the linearly elastic approximation on the right holds for small strains.

Note that the right-hand side of this equation is the same as that of equation A-13. In particular, 

for a uniaxial stress test, equation A-14 reduces to equation A-13 since rr   in that case. 

However, for the more general case considered here, rr need not be zero, and equation A-14 is 

only valid on the centerline.

A-3 The Compressible Version of the Mooney-Rivlin Model in LS-DYNA

The simulations described in sections 6 and 7 used LS-DYNA’s compressible version of the 

Mooney-Rivlin model for the specimen; this is Material Type 27 in (8, 9). The model has three 

material constants (aside from the density): the elastic moduli45 A1 and A2 and the Poisson’s ratio

 The strain energy (per unit undeformed volume) is given by

       

2 2

1 1 2 2 3 3 W A I A I C I D I 3 3 1 1 

        , (A-15)

 

45 The LS-DYNA manuals (8, 9) use the symbols A and B for the constants A1 and A2, respectively.
82

which differs from the incompressible model by the addition of the terms with coefficients C and 

D. The constant C is determined from A1 and A2 by

1

2

2

A

C A   . (A-16)

The constant D is determined from A1, A2, and  by

1 2 (5 2) (11 5)

2(1 2 )

  





A A

D

 



. (A-17)

We refer to this as the “compressible Mooney-Rivlin model”.

The Cauchy stress tensor corresponding to this strain energy function is46

1 1

2

2 2 Ap A J J

      I B B , (A-18)

where

 

3 2

5 2

4 2

4

C A

p D J J I

J J

    . (A-19)

When J  1, equation 18 yields the relations A.10 and A.11 for the principal stresses, which were 

obtained previously for the incompressible model. If J is close to 1, then these relations hold as 

approximations. Similarly, when J  1, we recover the relation A-13 for the axial stress zz in a 

uniaxial stress test and the relation A-14 for the stress difference zz  rr on the centerline for a 

general axisymmetric deformation. If J is close to 1, then47

1/   r z

, and we find that A-13 

and A-14 hold as approximations. In this regard, note that 

3G E 

for a nearly incompressible 

material; indeed, for the ballistic gelatin specimen they agree to at least three significant figures 

(see table 1).

Note that the scalar 

p

in equations A-18–A-19 is generally not equal to the pressure p, since the 

other terms in equation A-18 contribute to the trace of . Indeed, from equations A-6, A-18, 

A-5, and A-19, we obtain

1 2

1 2

2 2

3 3

A A p p I I

J J

   . (A-20)

 

46 None of the relations in this paragraph appear in the LS-DYNA User’s Manual (8) or in the Theory Manual by Hallquist

(9). However, equation A-18 is equivalent to the relation for the second Piola-Kirchhoff stress tensor given in (9), and equations 

A-19–A-20 imply the relation for the pressure in a state of pure dilation given in (9).

47 See equation 15 in section 2.6.
83

Thus the pressure p is a complicated function of both the volumetric strain and the shear strain. 

In particular, p does not necessarily reduce to zero when there is no volume change (i.e., when 

J = 1). However, p does equal zero in the undeformed state, as expected. This follows from 

equations A-20, A-19, A-16, and the fact that J = 1 and I1 = I2 = 3 when the material is 

undeformed.

For nearly incompressible materials ( close to 1⁄2), A1 and A2 (and hence C) are small compared 

to D, so for small volume changes (J close to but not equal to 1), equations A-19 and A-20 yield 

the simple approximations

p p K J K D     (1 ), 8 , (A-21)

where the constant K is the initial (or linear elastic) bulk modulus. It is easily verified that the 

approximation for K in equation A-21 is consistent with equations 21, 16, and A-17 when  is 

close to 1⁄2.

The LS-DYNA manuals do not provide any explanation or references for the unconventional 

form of the two compressibility terms in the strain energy function in equation A-15, that is, the 

terms with coefficients C and D, although the term involving C appears to have been added so 

that the pressure reduces to zero in the undeformed state. However, the relation for the pressure 

obtained from this strain energy function (equations A-19–A-20) is inappropriate for large 

volumetric strains, since the 

3

J J  term in equation A-19 has an absolute maximum at

1

0.58

3

J   . (A-22)

If  is sufficiently close to 1⁄2 (e.g., 



0.4999), the

3

J J  term dominates except at very large 

compressions (i.e., J  0.5).48 This physically unrealistic behavior was not an issue for the 

SHPB simulations reported here. For the parameters used in our calibration of the model (see 

table 3), the value of the pressure at the maximum described above is about 770 MPa, whereas 

the peak pressures in our simulations were around 4 MPa, with the exception of the sharp spikes 

on the centerline which did not exceed 10 MPa.49 In other words, the volume changes were 

small enough that the approximation A-21 was valid.

A-4 Verification of the Compressible Mooney-Rivlin Model

Prior to our simulations of SHPB tests, we performed a partial verification of the implementation 

of the compressible version of the Mooney-Rivlin model in LS-DYNA. We used the fact that 

for a quasi-static uniaxial stress test, the stress predicted by the compressible model should 

approach that predicted by the incompressible model as Poisson’s ratio  approaches 1⁄2 with 

A1 and A2 fixed. The constants A1 and A2 were chosen as in table 3. We used a single 8-node 

 

48 This unrealistic behavior is caused by the last term in equation A.15; it could be eliminated by replacing I3 there by its 

square root, that is, by J, but then the term with coefficient C would also need to be modified to guarantee that the pressure is 

zero in the undeformed state.

49 Refer to figures 11, 12, 17–23, 27, and 28, and recall that in the axial stress history plots, the stress state at the stress peaks 

prior to the formation of gaps is nearly hydrostatic.
84

hexagon element that was initially cubic with a 100 m edge length. In order to diminish any 

inertial effects, the initial density was reduced from 1 to 0.001 g/cm

3

for these simulations only.

There were no discernable differences in the results for nominal axial strain rates of 100/s, 

1000/s, and 10,000/s.

For various values of Poisson’s ratio  between 0.495 and 0.49999, the stress components were 

computed for axial stretches z ranging from 0.2 (large compression) to 4.0 (large extension). 

These results were compared with the theoretical relation A-13 for the incompressible model. 

The results for   0.495 are shown in figure A-1. The stress-stretch curves for the axial stress 

zz for the compressible model (red) and incompressible model (dashed line) are 

indistinguishable except for a slight difference when z is below 0.35. This difference decreased 

for higher values of Poisson’s ratio; in particular, for the value   0.49999 used in the SHPB 

simulations, the stress-stretch curves were indistinguishable for all values of z. As expected, the 

other stress components50 (the blue line in the figure) were indistinguishable from zero on this 

scale.

Figure A-1. Stress components as a function of axial stretch z for a uniaxial stress 

test on a single element using the compressible Mooney-Rivlin model 

with   0.495 (colored curves). Comparison with axial stress zz for the 

incompressible model (dashed curve).

Finally, note that for small volume changes (

J 1

), the curve in figure A-1 also represents the

stress difference zz  rr on the centerline for a general axisymmetric deformation.51

 

50 These components are taken relative to a Cartesian coordinate system aligned with the (initially cubic) specimen.

51 See equation A-14 and the discussion following equation A-19.
85

Appendix B. Sensitivity Studies

All of the results presented in the main body of this report were obtained with the baseline 

25×50-m mesh in the specimen; see section 4.1 and figure 2a. All of the results presented in 

the main body were obtained with the LS-DYNA default values of the contact algorithm 

parameters (SFAC = 1, VDC = 10); see section 4.2. Here we examine the sensitivity of the gap 

size on the centerline to the size of the specimen mesh (appendix B-1) and the contact algorithm 

parameters (appendix B-2).

For the mesh sensitivity study the default values of the contact algorithm parameters are used. 

Likewise, for the contact algorithm sensitivity study, the baseline mesh is used. For all of the 

results in this appendix we used the compressible Mooney-Rivlin constitutive model for the 

specimen and the larger (12.8-mm) bar radius. Both the 1- and 25-s initial rise times for the 

loading wave are considered.

B-1 Mesh Sensitivity

Refer to section 4.1 for a discussion of the specimen and bar meshes. Figures B-1–B-4 compare 

the histories of the gap size on the centerline for three specimen meshes: the baseline 

25×50 m mesh (figure 2a), the coarser 50×100m mesh (figure 2b), and the finer 

12.3×12.5 m mesh (figure 2c). The mesh size in the bars is 100×100 m for all simulations. 

Figures B-1 and B-2 are for the 1-s initial rise time for the loading wave; figures B-3 and B-4 

are for the 25-s initial rise time. Figures B-1 and B-3 give the gap size at the S-IB interface; 

figures B-2 and B-4 give the gap size at the S-TB interface.
86

Figure B-1. Histories of the gap sizes on the centerline at the S-IB interface for three specimen meshes 

and a 1-s initial rise time for the loading wave.
87

Figure B-2. Histories of the gap sizes on the centerline at the S-TB interface for three specimen 

meshes and a 1-s initial rise time for the loading wave.
88

Figure B-3. Histories of the gap sizes on the centerline at the S-IB interface for three specimen meshes

and a 25-s initial rise time for the loading wave.
89

Figure B-4. Histories of the gap sizes on the centerline at the S-TB interface for three specimen meshes 

and a 25-s initial rise time for the loading wave.

First, we consider the 1-s initial rise time for the loading wave. A detailed discussion of the 

results for this case using the baseline mesh was given in section 6.1, where the gap formation 

was broken down into six stages. At both interfaces, there is mesh convergence in the onset of 

the gaps (at 165 s). For the gap at the S-IB interface (figure B-1), we see that there is 

substantial mesh convergence in both gap size and duration for Stage II (165–280 s, during 

which there is a persistent gap) and Stage III (280–286 s, during which the gap remains closed). 

For the gap at the S-TB interface (figure B-2), we see that there is substantial mesh convergence 

in the duration of Stages II and III, but no convergence in the size of the gap during Stage II. 

However, the results for all three meshes are qualitatively similar during Stage II, that is, the gap 

opens and closes several times. For Stage IV (286–324 s, during which both gaps re-open), 

there is a degree of convergence in the duration of the gaps at either interface but not in gap size. 

Convergence has not been achieved at later times.
90

Next, we consider the 25-s initial rise time for the loading wave. A discussion of the results for 

this case using the baseline mesh was given in section 7. At both interfaces, there is mesh 

convergence in the onset of the gaps. For the gap at the S-IB interface (figure B-3), we see that 

there is substantial mesh convergence in both gap size and duration during the Stage II 

(184–218 s, during which there is a persistent gap); in particular, the timing of the peak value 

(2.2–2.4 m) agrees within 1 s for all three meshes. For other times at the S-IB interface and 

all times at the S-TB interface (figure B-4), the gap sizes remain under about 1.2 m; the local 

peaks in gap size generally occur at similar times; and the gaps persist for similar durations for 

the three meshes. However, gap sizes generally display less mesh convergence than gap 

durations.

The results above support the conclusion that the gap phenomenon is not a numerical artifact 

introduced by an insufficiently fine mesh in the specimen.

B-2 Sensitivity to the Contact Algorithm Parameters

Refer to section 4.2 for a brief discussion of the contact algorithm. Figures B-5, B-6, B-9, and 

B-10 compare the histories of the gap size on the centerline for three values of the scale factor 

SFACT: 0.1, 1 and 10, which correspond to a less stiff, the default, and a stiffer penalty force, 

respectively. Figures B-7, B-8, B-11, and B-12 compare the histories of the gap size on the 

centerline for two values of the viscous damping coefficient VDC: 10 and 20, which correspond 

to the default damping and twice the default damping, respectively. Figures B-5–B-8 are for the 

1-s initial rise time for the loading wave; figures B-9–B-12 are for the 25-s initial rise time. 

The results at the S-IB interface are given in figures B-5, B-7, B-9, and B-11; those for the S-TB 

interface are given in figures B-6, B-8, B-10, and B-12.

First, we consider the 1-s initial rise time for the loading wave (figures B-5–B-8). A detailed 

discussion of the results for this case using the default values for the contact parameters was 

given in section 6.1, where the gap formation was broken down into six stages. The durations of 

Stage II (165–280 s), Stage III (280–286 s) and Stage IV (286–324 s) are unaffected by the 

changes in SFACT and VDC. The size of the gap at the S-IB interface during Stage II 

(165–280 s) is also unaffected by the changes in SFACT and VDC. At other times, the size of 

the gap at both interfaces varies somewhat with the values of these parameters, but the results are 

qualitatively similar.
91

Figure B-5. Histories of the gap sizes on the centerline at the S-IB interface for three values of the 

contact parameter SFACT and a 1-s initial rise time for the loading wave.
92

Figure B-6. Histories of the gap sizes on the centerline at the S-TB interface for three values of the 

contact parameter SFACT and a 1-s initial rise time for the loading wave.
93

Figure B-7. Histories of the gap sizes on the centerline at the S-IB interface for two values of the 

contact parameter VDC and a 1-s initial rise time for the loading wave.
94

Figure B-8. Histories of the gap sizes on the centerline at the S-TB interface for two values of the 

contact parameter VDC and a 1-s initial rise time for the loading wave.

Next, we consider the 25-s initial rise time for the loading wave (figures B-9–B-12). A 

discussion of the results for this case using the default values for the contact parameters was 

given in section 7. At the S-IB interface, we see that doubling VDC has no effect on the duration 

or size of the gap during the Stage II (184–218 s). The changes in SFACT also have no effect 

on the duration and only a slight effect on the size of the gap during this stage; in particular, the 

peak value occurs at about 200 s for all three cases. For other times at the S-IB interface and all 

times at the S-TB interface, doubling VDC has little effect on the duration of the gaps and only a 

slight effect on their size. The changes in SFACT have more effect on both duration and size, 

but the results are qualitatively similar.

The results above support the conclusion that the gap phenomenon is not a numerical artifact

introduced by the contact algorithm.
95

Figure B-9. Histories of the gap sizes on the centerline at the S-IB interface for three values of the 

contact parameter SFACT and a 25-s initial rise time for the loading wave.
96

Figure B-10. Histories of the gap sizes on the centerline at the S-TB interface for three values of the 

contact parameter SFACT and a 25-s initial rise time for the loading wave.
97

Figure B-11. Histories of the gap sizes on the centerline at the S-IB interface for two values of the 

contact parameter VDC and a 25-s initial rise time for the loading wave.
98

Figure B-12. Histories of the gap sizes on the centerline at the S-TB interface for two values of the 

contact parameter VDC and a 25-s initial rise time for the loading wave.
99

Appendix C. Additional Computational Studies

This section contains summaries of additional numerical simulations in which gaps formed at the 

specimen-bar interfaces. These simulations involve conditions that differ in one or more ways 

from those considered previously, and in appendix C-3 they involve the use of a different code.

C-1 Direct Impact on the Specimen

Recall that in a conventional SHPB test, the striker bar impacts the incident bar and generates a 

compressive stress pulse which propagates along the incident bar and subsequently into the 

specimen and the transmission bar; and for soft specimens, a pulse shaper is often inserted 

between the striker and incident bars to produce a smoother loading wave. A variation of the 

SHPB test, known as the direct impact Hopkinson bar or direct impact compression test, involves 

the elimination of the incident bar, so that the striker bar directly impacts the specimen (27–30). 

The transmission bar is retained and serves the same purpose as in a conventional SHPB test. 

The main advantage of this technique is that it permits higher strain rates to be imposed on the 

specimen. There are some disadvantages as well; see the papers cited above and also the 

discussions in Gray (1) and Jia and Ramesh (31).

We performed several simulations in which a striker bar directly impacts the specimen. The 

specimen geometry, bar diameters, bar properties, and mesh sizes were the default values used 

previously (sections 2.1, 2.3, and 4.1). The compressible Mooney-Rivlin model (section 3) was 

used for the specimen. The initial velocity of the striker bar is vSB = 3.625 m/s, which yields a 

nominal axial strain rate of approximately vSB/LS = 2500/s, as before. Note that vSB is twice the 

plateau velocity v0 imposed at the far end of the incident bar in the previous SHPB simulations; 

on the other hand, the particle velocity in the incident bar nearly doubles as the loading wave 

reflects from the S-IB interface (section 5.3). Thus, in either case an axial velocity of roughly 

3.625 m/s is ultimately imposed on the face of the specimen. However, in the SHPB simulations 

the imposed velocity ramps up smoothly from zero (figures 5 and 6), whereas the direct impact 

imposes an initial step in velocity (and hence an acceleration impulse) on the specimen. In the 

subsequent discussion, the SHPB simulation mentioned for purposes of comparison will be the 

case with a 1-s initial rise time (section 6), as this case imposed the most severe accelerations 

on the specimen.

In the SHPB simulation, the velocity at the S-IB interface exhibits large oscillations as a result of 

wave dispersion in the incident bar (figures 3 and 5). This dispersion effect will be absent for 

direct impact. Nevertheless, there will be some slight deceleration of the striker bar as well as 

(possibly large) oscillations in velocity throughout the specimen at early times due to multiple 

reflections from the specimen-bar interfaces. Unfortunately, we do not have any plots of the 

velocity history for this case.
100

In the SHPB simulation, gaps formed at the S-IB and S-TB interfaces about 25 s after the 

arrival of the loading wave at the S-IB interface. In the direct impact simulation, gaps formed 

much sooner at both the specimen-striker bar and S-TB interfaces. In fact, by 6 s after impact, 

the entire face of the specimen had separated from the transmission bar. This complete 

separation continued for the short duration of the simulation, which was stopped at about 13 s 

after impact, corresponding to a compressive axial strain of about 3%. Prior to gap formation, 

the peak compressive axial stress in the specimen was about 8 times that for the SHPB 

simulation. Also, tensile axial stresses on the centerline were first observed in the transmission 

bar, analogous to the SHPB simulation (section 6.3).

Our goal in doing these direct impact simulations was not to study the direct impact compression 

test per se, but rather to examine the effect of a different loading history on the formation of 

gaps. However, the results described above call into question the use of direct impact on soft, 

nearly incompressible specimens. Indeed, whereas inertial effects at early times in SHPB tests 

can be reduced (though not eliminated) with better pulse shaping, for a direct impact test there is 

no way to reduce the acceleration impulse imparted to the specimen.

C-1.1 Lateral Constraints

The lateral surfaces of the specimen and the bars are unconstrained and hence stress-free in both 

the SHPB and direct impact compression tests. We performed two additional direct impact 

simulations to examine the effect of these free surfaces on gap formation. The conditions in both 

simulations were as described above for the direct impact test, with the exceptions noted below.

In the first simulation, the lateral surfaces of the bars and the specimen were constrained to 

prevent any radial motion, whereas the axial component of the motion at these surfaces was 

unconstrained as before. No gaps formed in this case. Axial stress contours at 10 s revealed a 

uniform compressive stress state throughout the specimen and the bars (in the neighborhood of 

the specimen), with the exception of a small region of tensile stress adjacent to the portion of the 

faces of the bars overhanging the specimen. These results are not surprising. Indeed, under 

these conditions, the specimen essentially undergoes a compressive, uniaxial strain, just as in a 

normal plate impact test. Since the radial motion throughout the specimen is eliminated for the 

most part, so are the associated radial inertia effects.

In the second simulation, the lateral constraint on the specimen was removed, but the constraint 

on the bars was retained. Gaps formed at both interfaces. The results were similar to those for 

the unconstrained direct impact simulation described previously.

These two simulations support the conclusion that radial inertia effects are responsible for gap 

formation in the direct impact tests as well.
101

C-2 Annular Specimens

All of the results presented up to this point have been for solid (disc shaped) specimens. Annular 

(washer shaped) specimens have also been used in SHPB tests on soft materials in an effort to 

reduce inertial effects (5, 6). The motivation for hollowing out the specimen comes from 

approximate theoretical estimates for the inertially induced radial stress in solid specimens, 

which predict that the radial stress is largest on the centerline and decreases parabolically to zero 

at the lateral (stress-free) surface (6, 7, 26, 27, 32, 33). Furthermore, numerical simulations 

reveal that in situations where gaps do not form but the axial acceleration is sufficiently large, 

the pressure at the center of a solid specimen may undergo large oscillations (32, 34). Therefore, 

removing the center portion of the specimen might reduce these inertial effects, resulting in a 

stress state that is closer to uniaxial. Experimental results indicate that this works for early times 

(5, 6).52

We performed several numerical simulations on an annular specimen with the same initial 

thickness (LS = 1.45 mm) and initial outer radius (RS = 6.35 mm) as the solid specimens 

discussed previously. The inner radius of the annular specimen was taken to be 2.638 mm. As 

before, we considered two values for the radius of the pressure bars, RB = 12.8 mm and 

RB = 9.5 mm. The dimensions of the annular specimen and the smaller bar radius were chosen to 

agree with those used in the experimental study on ballistic gelatin by Moy et al. (5).53 The 

compressible Mooney-Rivlin model was used for all annular specimen simulations.

For the larger diameter pressure bars and a 25-s initial rise time for the loading wave, the total 

force acting on a cross-section of the incident and transmission bars was measured at an axial 

location of 2.52 mm from the specimen-bar interfaces (cf. section 7.3). The force history in the 

transmission bar was qualitatively similar to that for the solid specimen (figure 30), although the 

negative forces (indicating gap opening) persisted for a shorter time. The force history in the 

incident bar oscillated between positive and negative from 157 s to the end of the simulation at 

340 s (43% nominal axial strain), indicating repeated gap opening and closure.

For the smaller diameter pressure bars and a 25-s initial rise time, no gaps were observed 

except momentarily (at about 7% axial strain) at the S-IB interface over a small region midway 

between the inner and outer surfaces.

For the smaller diameter pressure bars and a 1-s initial rise time for the loading wave, a 

substantial gap was observed at the S-IB interface but not the S-TB interface. At t  182 s and 

an axial strain of 7%, this gap extends over 3⁄4 of the S-IB interface, as seen in figure C-1. This 

figure also reveals bulging of the annular specimen at the outer (top) and inner (bottom) stressfree lateral surfaces, near the S-IB interface. The centerline of the specimen and bars lies outside 

 

52 However, theoretical and computational results demonstrate that the presence of a stress-free inner surface leads to 

additional inertial effects at later times (32, 36). This complication will not be discussed further here.

53 Actually, the inner radius of the annular specimen used here is 1% smaller than the value used in Moy et al. (5). 

The reason for this slight difference is that it allowed us to use the same mesh size for the solid and annular 

specimens.
102

of (below) the figure. Figure C-1 also shows contours of the pressure. In order to reveal more 

clearly the regions of negative (tensile) pressure, any point with zero or positive pressure 

(compression) is shaded red; any color other than red indicates a tensile state. The figure reveals 

a region of tensile pressure in the specimen bordering the gap at the S-IB interface.

Figure C-1. Pressure contours in the vicinity of an annular specimen at 7% axial strain. Any 

positive pressure (compression) is shaded red; other colors indicate negative (tensile) 

pressure in units of GPa. The loading wave had a 1-s initial rise time. Observe the 

gap along most of the S-IB interface.
103

C-3 Simulations with the PRESTO Code

With the exception of the results to be discussed in this section, all of the numerical simulations 

in this report used the commercial finite-element code LS-DYNA (8, 9) in the 2-D axisymmetric

mode. In an effort to verify that the formation of gaps at the specimen-bar interfaces is not some 

artifact of this particular code, Bryan Love (37) performed analogous numerical simulations with 

PRESTO, a finite-element code from Sandia National Laboratories. All of the PRESTO 

simulations were fully three-dimensional, which necessitated larger element sizes than those 

used in the LS-DYNA simulations. The two codes also have different contact algorithms.

A compressible version of the Mooney-Rivlin model was used for the specimen, with the same 

calibration used in the LS-DYNA simulations (see section 3).54 Both solid and annular 

specimens were considered; these had the same dimensions as in the LS-DYNA simulations.

The diameter of the aluminum pressure bars used in these simulations was 1 inch (25.4 mm), 

which is essentially the larger of the two diameters considered previously. The pressure bars 

were much longer (2438 mm) than in the LS-DYNA simulations, corresponding to the bar 

lengths in an actual Hopkinson bar setup (25). The loading condition in the PRESTO 

simulations involved an aluminum striker bar (of the same diameter) impacting the far end of the 

incident bar, as in a real SHPB test. However, no pulse shaper was used between the striker and 

incident bars. A few direct impact simulations were also performed.

The results of these PRESTO simulations confirmed the findings of the previous LS-DYNA 

simulations with regard to the formation of gaps at the specimen-bar interfaces. The fact that the 

PRESTO simulations were fully three-dimensional rules out the possibility that gap formation in 

the LS-DYNA simulations was in some way associated with issues at the centerline in the 2-D 

axisymmetric mode.

 

54 However, the PRESTO version of the compressible Mooney-Rivlin model does not suffer from the unrealistic behavior of 

the pressure at very large volumetric compression and high Poisson’s ratio that is exhibited by the LS-DYNA version

(see appendix A-3).
104

List of Symbols, Abbreviations, and Acronyms

A1, A2 constants in the Mooney-Rivlin model

cL, cG, cE longitudinal, shear, and bar wave speeds

E Young’s modulus 

ez nominal axial strain in the specimen

Fz total force on a cross-section of the transmission bar

G shear modulus

IB incident bar

J Jacobian of the deformation

K bulk modulus

L longitudinal modulus

LIB length of the incident bar

LS initial length (thickness) of the specimen

LTB length of the transmission bar

p pressure

r radial coordinate of the deformed material 

R radial coordinate of the undeformed material 

RB radius of the incident and transmission bars

RS initial radius of the specimen

S specimen

SFACT scale factor for the penalty force stiffness in the contact algorithm

SHPB split Hopkinson pressure bar

S-IB specimen-incident bar 

S-TB specimen-transmission bar 

t time

TB transmission bar

tR rise time of the axial velocity prescribed at the end of the incident bar
105

uIB, uS-IB axial displacements of the incident bar and specimen at the S-IB interface

uTB, uS-TB axial displacements of the transmission bar and specimen at the S-TB interface

v0 plateau value of the axial velocity prescribed at the end of the incident bar

VDC viscous damping coefficient in the contact algorithm

vIB, vTB axial velocity of the incident and transmission bars at the specimen interface

vz axial velocity of the incident bar

z axial coordinate of the deformed material 

Z axial coordinate of the undeformed material 

zIB, zS-IB deformed axial coordinates of the incident bar and specimen at the S-IB interface

zTB, zS-TB deformed axial coordinates of the transmission bar and specimen at the S-TB 

interface

S-IB gap size at the specimen-incident bar interface

S-TB gap size at the specimen-transmission bar interface

 Poisson’s ratio 

S deformed length (thickness) of the specimen

z, r,  axial, radial, and hoop stretches in the specimen

 Cauchy (true) stress tensor

 initial density 

zz, rr,  axial , radial, and hoop components of the Cauchy stress
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

106

1 DEFENSE TECHNICAL

(PDF INFORMATION CTR

ONLY) 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

3 INSTITUTE FOR ADVANCED 

TECHNOLOGY

 S BLESS

 S SATAPATHY

 S LEVINSON

 3925 WEST BRAKER ST

SUITE 400

AUSTIN TX 78759-5316

3 SOUTHWEST RESEARCH INSTITUTE

 C A ANDERSON

 J WALKER

 K DANNEMANN

6220 CULEBRA RD

DRAWER 28510

SAN ANTONIO TX 78228-0510

4 SOUTHWEST RESEARCH INSTITUTE

 G JOHNSON

 T HOLMQUIST

 S BEISSEL

 C GERLACH

5353 WAYZATA BLVD

MINNEAPOLIS MN 55416-1340

1 INTERNATIONAL RESEARCH 

ASSOCIATES

 D ORPHAL

4450 BLACK AVE STE E

PLEASANTON CA 94566-6145

4 JOHNS HOPKINS UNIVERSITY

MECHANICAL ENGINEERING

 K T RAMESH

 A DOUGLAS

 K HEMKER

 T NGUYEN

223 LATROBE HALL

3400 N CHARLES ST

BALTIMORE MD 21218-2682

1 JOHNS HOPKINS UNIVERSITY

CIVIL ENGINEERING

 L GRAHAM-BRADY

210 LATROBE HALL

3400 N CHARLES ST

BALTIMORE MD 21218-2682

1 VIRGINIA POLYTECHNIC INST

COLLEGE OF ENGRG

 R BATRA

BLACKSBURG, VA 24061-0219

1 LIVERMORE SOFTWARE

TECHNOLOGY CORP

 M JENSEN

7374 LAS POSITAS RD

LIVERMORE CA 94550

1 UNIVERSITY OF ALABAMA AT 

BIRMINGHAM

MECHANICAL ENGINEERING

 D LITTLEFIELD

HOEN ENGINEERING BLDG

1530 3RD AVE S

BIRMINGHAM AL 35294-4440
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

107

1 PURDUE UNIVERSITY

AERONAUTICS & ASTRONAUTICS

 W CHEN

701 W STADIUM AVE

W. LAFAYETTE IN 47907-2045

1 UNIVERSITY OF RHODE ISLAND

MECHANICAL ENGINEERING

 A SHUKLA

206 WALES HALL

KINGSTON RI 02881

1 UNIVERSITY OF MISSISSIPPI

MECHANICAL ENGINEERING

 A RAJENDRAN

201 CARRIER HALL

UNIVERSITY MS 38677

2 UNIVERSITY OF CALIFORNIA 

SAN DIEGO

MECHANICAL AND AEROSPACE 

ENGINEERING

 M MEYER

 V NESTERENKO

9500 GILMAN DR

LA JOLLA CA 92093-0411

1 UNIVERSITY OF CALIFORNIA 

SAN DIEGO

DEPT APPLIED MECH & ENG SCI

 S NEMAT-NASSER

9500 GILMAN DR

LA JOLLA CA 92093-0416

1 UNIVERSITY OF UTAH

MECHANICAL ENGINEERING

 R BRANNON

50 S CENTRAL CAMPUS DR

SALT LAKE CITY UT 84112

1 BROWN UNIVERSITY

ENGINEERING DIV

R CLIFTON

182 HOPE ST

PROVIDENCE, RI 02912

1 UNIVERSITY OF TEXAS-AUSTIN

AEROSPACE ENGINEERING AND 

ENGINEERING MECHANICS

 K RAVI-CHANDAR

1 UNIVERSITY STATION, C0600

AUSTIN TX 78712-0235

2 UNIVERSITY OF NEBRASKALINCOLN

DEPT OF ENGINEERING MECHANICS

 R FENG

 M NEGAHBAN

LINCOLN NE 68588-0526

1 STANFORD UNIVERSITY 

MECHANICS AND COMPUTATION 

GROUP

ATTN A LEW

WILLIAM F DURAND BLDG

STANFORD CA 94305

1 UNIV OF DELAWARE

CTR FOR COMPOSITE MTRLS

ATTN B GAMA

201 SPENCER LABORATORY

NEWARK DE 19716

1 UNIV OF DELAWARE

CTR OF COMPOSITE MTRLS

ATTN S LOPATNIKOV

215 CONPOSITES MANUFAC SCI LAB

NEWARK DE 19716

1 UNIV OF DELAWARE CTR FOR 

COMPOSITE MATERIALS

ATTN J GILLESPIE

201 SPENCER LAB

NEWARK DE 19716

1 UNIVERSITY OF FLORIDA

MECHANICAL AND AEROSPACE 

ENGINEERING

ATTN G SUBHASH

GAINESVILLE FL 32611

1 UNIVERSITY OF WATERLOO

DEPT OF MECHANICAL 

ENGINEERING

ATTN C P SALISBURY

WATERLOO ONTARIO CANADA 

2 AIR FORCE RESEARCH 

LABORATORY

RWMW

 J L JORDAN

 L CHHABILDAS

101 WEST EGLIN BLVD

EGLIN AFB FL 32542
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

108

5 SANDIA NATIONAL LABORATORIES

 M R BAER MS 0836

 S SILLING MS 1322

 E S HERTEL MS 1185

 S SCHUMACHER MS 0836

 P A TAYLOR MS 1160

P O BOX 5800

ALBUQUERQUE NM 97185

2 SANDIA NATIONAL LABORATORIES

 T J VOGLER MS 9402

 B SONG MS 9404

P O BOX 969

LIVERMORE CA 94551-0969

8 LOS ALAMOS NATIONAL 

LABORATORY

 E N BROWN

 G GRAY III

 E MAS

 D DATTLEBAUM

 T MASON

 P RAE

 J WALTER

 B E CLEMENTS

P O BOX 1663

LOS ALAMOS, NM 87545

2 LAWRENCE LIVERMORE NATIONAL 

LABORATORY

 D GOTO

 M MURPHY

P O BOX 808

LIVERMORE, CA 94551-0808

1 APPLIED RESEARCH ASSOCIATES

SOUTHWEST DIVISION

 D E GRADY

4300 SAN MATEO BLVD NE

ALBUQUERQUE NM 87110

1 T W WRIGHT

4906 WILMSLOW RD

BALTIMORE MD 21210

1 M J FORRESTAL

1805 NEWTON PL NE

ALBUQUERQUE NM 87106

ARMY RESEARCH OFFICE:

1 RDRL ROE M

 D STEPP

1 RDRL ROE N

 B LAMATTINA

ABERDEEN PROVING GROUND

93 DIR USARL

RDRL CIM G (BLDG 4600)

RDRL CIH C

 D GROVE

 J CAZAMIAS

 P CHUNG

 J KNAP

RDRL SL

 R COATES

RDRL SLB W

 W MERMAGEN

RDRL WM

 P PLOSTINS

 J MCCAULEY

 J SMITH

RDRL WML

 J NEWILL

 M ZOLTOSKI

RDRL WML H 

 M FERMEN-COKER 

 L MAGNESS

 D SCHEFFLER

 S SCHRAML

RDRL WMM

 J BEATTY

 R DOWDING

RDRL WMM B

 B CHEESEMAN

 C FOUNTZOULIS

 G GAZONAS

 D HOPKINS

 P MOY

 B POWERS

 C RANDOW 

 M VANLANDINGHAM

 R WILDMAN

RDRL WMM C

 R JENSEN

RDRL WMM D

 E CHIN

RDRL WMM G

 J LENHART
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

109

RDRL WMP

 P BAKER 

 S SCHOENFELD 

RDRL WMP A

 B RINGERS

RDRL WMP B

 C ADAMS

 R BECKER

 S BILYK 

 J BRADLEY

 D CASEM 

 J CLAYTON

 A DAGRO 

 D DANDEKAR 

 A DWIVEDI

 M GREENFIELD

 C A GUNNARSSON 

RDRL WMP B

 C HOPPEL

 J HOUSKAMP

 Y HUANG 

 R KRAFT 

 B LEAVY

 B LOVE 

 M LYNCH

 J NIEMCZURA

 M RAFTENBERG (5 COPIES) 

 M SCHEIDLER (10 COPIES) 

 T WEERASOORIYA 

 C WILLIAMS

 K ZIEGLER

RDRL WMP C

 T BJERKE

 K KIMSEY

 S SEGLETES 

 W WALTERS

RDRL WMP D

 R DONEY

 J RUNYEON

 B SCOTT

RDRL WMP E

 P BARTKOWSKI

 M BURKINS

 T JONES

 E HORWATH

 C KRAUTHAUSER

 M LOVE

RDRL WMP F 

 R BITTING

 E FIORAVENTE

 D FOX

 A FRYDMAN

 N GNIAZDOWSKI 

 R GUPTA

RDRL WMP G 

 R BANTON

 N ELDREDGE

 S KUKUCK
110

INTENTIONALLY LEFT BLANK.