Evaluation of the Material Point Method within CTH to 

Model 2-Dimensional Plate Impact Problems

by Denzell Bolling

ARL-MR-0880 September 2014

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-5069

ARL-MR-0880 September 2014

Evaluation of the Material Point Method within CTH to 

Model 2-Dimensional Plate Impact Problems

Denzell Bolling

Weapons and Materials Research Directorate

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 2014

2. REPORT TYPE

Final

3. DATES COVERED (From - To)

21 May 2014–15 Aug 2014

4. TITLE AND SUBTITLE

Evaluation of the Material Point Method within CTH to Model 2-Dimensional

Plate Impact Problems

5a. CONTRACT NUMBER

W911NF-10-2-0076

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Denzell Bolling

5d. PROJECT NUMBER

SEAP-CQL

5e. TASK NUMBER

5f. WORK UNIT NUMBER

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

US Army Research Laboratory

ATTN: RDRL-WMP-C

Aberdeen Proving Ground, MD 21005-5069

8. PERFORMING ORGANIZATION

 REPORT NUMBER

ARL-MR-0880

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

Participated in the College Qualified Leaders Program at Howard University.

14. ABSTRACT

The material point method (MPM) is a mixed Eulerian and Lagrangian computational method that allows for the discretization 

of solid bodies into material points carrying the full material state of the original body. This method is being evaluated as a new 

implementation within CTH, a code originally developed by Sandia National Laboratories, as a way of eventually dealing with 

fracture and other difficult solid mechanics scenarios. Applying CTH with MPM to a 2-dimensional plate impact problem 

yields multiple issues that must be examined. Issues with void compression upon initial impact, inhomogeneous particle states 

within a single cell, and inconsistencies within results obtained at different mesh resolutions are all discussed.

15. SUBJECT TERMS

impact, particles, CTH, markers, material point method, particle in cell

16. SECURITY CLASSIFICATION OF: 

17. LIMITATION

OF ABSTRACT

UU

18. NUMBER

OF PAGES

18

19a. NAME OF RESPONSIBLE PERSON

Richard Becker

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

410-278-7980

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
iii

Contents

List of Figures iv

Acknowledgments v

Student Bio vi

1. Introduction/Background 1

2. Experiment/Calculations 2

2.1 Mesh ................................................................................................................................2

2.2 Diatom.............................................................................................................................2

2.3 Particles...........................................................................................................................2

2.4 Equation of State .............................................................................................................3

2.5 Tracers.............................................................................................................................3

3. Results and Discussion 4

3.1 Plate Contact....................................................................................................................4

3.2 Lowered Pressure and Mesh Resolution .........................................................................4

3.3 Tracer Location ...............................................................................................................5

3.4 Raised Pressure................................................................................................................6

3.5 Uintah: A Particle Code ..................................................................................................6

4. Summary and Conclusions 8

5. References 9

Distribution List 10
iv

List of Figures

Fig. 1 A continuum pressure contour plot as the impact wave travels after plate impact .............3

Fig. 2 Pressure impact image from CTH, where the red indicates particle areas of high 

pressure, and the green indicates areas of zero pressure (scale units in gigapascals)................4

Fig. 3 (Left) Coarse mesh and (right) fine mesh continuum comparisons ....................................5

Fig. 4 Mesh resolution comparison shows disagreement ..............................................................5

Fig. 5 Tracer 4 pressure comparison based on location.................................................................6

Fig. 6 Fine mesh continuum comparison after the change in tracer location ................................6

Fig. 7 Pressure contour from Uintah shows ringing as the impact wave travels outward.............7

Fig. 8 Individual particle pressure within consecutive cells after plate impact for Uintah and 

CTH with MPM.........................................................................................................................8
v

Acknowledgments

I wish to acknowledge the mentorship of Dr Richard Becker, Stephan Bilyk, Dr Shane 

Schumacher, Dr Kevin Rugirello, Dr Andrew Tonge, Dr Jeffrey Lloyd, Dr Mary Jane Graham, 

and Dr Gbadebo Owolabi.
vi

Student Bio

I am currently attending Howard University in Washington, DC, as a second-year graduate 

student within mechanical engineering. I also attended Howard University for my undergraduate 

degree in mechanical engineering. My past research experience includes an internship at Sandia 

National Laboratories studying the expansion of oil reserve reservoirs. I also worked on 2 

undergraduate research experiences at the Georgia Institute of Technology studying alternative 

energy. In the future, I plan on finishing my master’s degree and applying for PhD programs 

within mechanical engineering.
1

1. Introduction/Background

CTH is a continuum hydrocode originally developed by Sandia National Laboratories. The 

software system can model multidimensional, multimaterial, large deformation, and strong shock 

wave physics.

1 The existing code is used widely across the US Army Research Laboratory; 

increasing the code’s capabilities would be a great benefit. For example, in cases that involve 

material fracture, continuum Eulerian codes such as CTH do not track fractures cleanly and can 

smear cracks and change their extent or growth characteristics. These adverse numerical artifacts

can introduce error into computational results.

2

The recent implementation of the material point method (MPM) into CTH may ameliorate some 

of these difficulties. This method involves the discretization of solid bodies into Lagrangian 

material points. These points (also called markers) carry the full state of the material, including 

mass, density, velocity, stress, and strain. Interactions among adjacent cells of the numerical 

background grid enforce compatibility and momentum balance, and the constitutive relations are 

evaluated on the markers. The markers receive their strain rate from the continuum grid during 

the Lagrangian part of the step, and the stress and material history are updated for the markers. 

The material points remain fixed to the material as the grid is returned to its original spatial 

location during a remap step. The remapped markers provide the stress to the grid for the 

momentum calculation to complete the cycle. The transient coupling between the particles and 

the grid enables large deformations of the Lagrangian particles without mesh tangling. Simply 

stated here, the inclusion of MPM has the potential to assist in accurately solving all kinds of 

large deformation mechanics problems where the internal state of material is history dependent.

3

The incorporation of MPM within CTH creates a new capability to be evaluated within this 

report. One of the central differences between continuum CTH and CTH with MPM is the 

superposition of 2 separate pressure terms within the momentum equation (Eq. 1) used in the

analysis.

4

∆u = ∆t

ms

∫ −∇ ∙ (PI + Q)dV + ∆t

ms

∫ ∇ ∙ (σs + PI)dV + ∆t

ms

∫ ρsb dV . (1)

This equation can be broken down into 3 terms. The first term is the continuum term; it is the 

traditional way of computing accelerations induced by the pressure gradient on the background 

grid. The second term is the particle term; it holds the full stress tensor, σs, at the particle level.

The third term is a body force term, which does not play a role in the following computations. 

The additive contributions of the positive pressure at the particle level and the negative pressure 

at the continuum level play a significant role in how CTH with MPM computes acceleration of 

the grid. In a uniform field, the pressure contribution removed from the particle term is exactly 

replaced by the continuum pressure from the background grid. For inhomogeneous
2

deformation fields, the result is more complex, as will be suggested in this work. A thorough 

understanding of the CTH solution algorithm with markers will lead to greater confidence for 

using the markers in a wide range of Army problems.

2. Experiment/Calculations

The evaluation of CTH with MPM is conducted on 2 main levels. The first level is a direct 

response comparison of CTH with and without particles for a zero-obliquity impact problem 

between 2 similar materials. The second level of the evaluation investigates mesh sensitivity for 

the particle models using 2 different mesh resolutions, one fine and another coarse. A 

comparison between these 2 resolutions ensures that the code is converging to a physically 

reasonable solution based on different levels of mesh refinement. Comparisons between the 

different paired simulations briefly described here demonstrate how accurately the new 

capabilities of CTH operate. The input deck built for these tests is composed of multiple 

components, including the mesh, diatom, markers, tracers, and the equation of state needed to 

compute the stress state of the materials.

2.1 Mesh

The mesh is a user-supplied spatial discretization of the material geometries for the problem of 

interest. For a Lagrangian domain, the mesh moves with the material, and for a Eulerian domain, 

the material advects through a stationary mesh. Mesh creation for the current evaluation uses a 

Eulerian domain with a 2-dimensional (2-D) rectangular geometry. Its mesh extends from –0.01 

to 0.01 m along the flyer plate propagation direction. There are 400 elements of 5 × 10-5 m

dimension for the fine mesh and 100 elements of 2 × 10-4

m dimension for the coarse mesh. The 

mesh extends 0.002 m from the origin in the y-direction. There are 40 elements of 5 × 10-5 m 

height for the fine mesh and 10 elements of 2 × 10-4 m height for the coarse mesh. Symmetry 

conditions were applied for the 4 boundaries of the computational domain.

2.2 Diatom

The diatom block of the input deck focuses on construction of the bodies of material that make 

up the problem. Here the initial model is a symmetric 2-D plate impact experiment between 2

copper plates. The plates are infinite in the y-direction and measure 0.005 m in the x-direction. 

They travel across a gap of 0.0003 m toward each other at 500 m/s each until impact. 

2.3 Particles

The continuum solution for this problem is evaluated using the standard CTH code, and it is to 

be compared with the version of CTH that uses the MPM (CTH with MPM) particle solution. 

The particles are distributed 3 per cell in the x and y directions and completely represent the 

copper plates described in the diatom for the particle runs. The particles are capable of
3

transporting the full state of the material at any specific particle location, but the equation of state 

(EOS) is evaluated on the background grid. The particles can also analyze multiple material 

pressures and temperatures with the MMP4 option enabled in this portion of the input deck.

2.4 Equation of State 

The EOS used to analyze the problem is the Mie-Gruneisen EOS. The strength of the copper is

taken from the Johnson-Cook material parameter database with a Poisson’s ratio of 0.3, density 

of 8950 kg/m3

, and a sound speed of 3560 m/s.

2.5 Tracers

Eight tracers are originally placed along the lower symmetry boundary of the left copper plate in 

all of the simulations. These tracers point to the locations along the body on which resultant data 

will be saved, recorded, and later plotted. Recorded variables on the tracers include time, stress, 

position, pressure, temperature, and density. Each of these values is recorded every 10 ns during the 

simulation. 

In addition to using an appropriate EOS and placing tracers, it was also necessary to compute the 

analytical equilibrium pressure for this plate impact. This calculation establishes a standard 

baseline around which computational agreement can be met. This pressure was evaluated using 

Eq. 2. It multiplies the density, ρ, of the copper plate with the longitudinal wave speed, CL, and 

velocity, V, of the plate.

5 An equilibrium pressure of 2.1 × 104 MPa was calculated.

P = ρ ∗ CL ∗ V . (2)

Not only should the numerically calculated peak pressure value result in an answer on the order 

of the analytic solution (see Eq. 2), but it should also exhibit a waveform like the one shown in 

Fig. 1. In Fig. 1 a continuum impact of 2 copper plates has taken place, and the resulting wave 

has a relatively smooth high-pressure response as it propagates through the plates.

Fig. 1 A continuum pressure contour plot as the impact wave travels after plate impact
4

3. Results and Discussion

3.1 Plate Contact

It is immediately evident that CTH with MPM is vulnerable to one of the inherent problems with 

current particle methods: void/particle mixed cell issues. Inspecting the pressure contours from 

CTH for the initial test shows that the 2 copper plates are not entirely coming into contact 

(Fig. 2). However, particles away from the impact face experience a pressure increase indicative 

of impact. The mixed cells at the interface of the 2 plates contain both particles and void space. 

An investigation into this part of the result will be discussed in more depth for a less complicated 

problem later in this report. 

Fig. 2 Pressure impact image from CTH, where the red indicates particle areas of high pressure,

and the green indicates areas of zero pressure (scale units in gigapascals)

Upon further investigation, it was also discovered that the particles at the impact face were at 

zero pressure despite the impact having already taken place, as shown by the green markers in 

Fig. 2. This particular issue may be linked to the superposition of the particle and background 

grid pressures discussed previously. Determining why particles are indicating zero pressure at 

the impact face requires further investigation, but the issue associated with the void gap between 

the plates is mitigated by physically closing the gap within the problem setup. The plates are

placed directly next to one another and move toward each other for impact.

3.2 Lowered Pressure and Mesh Resolution

With the void gap issues having been noted and avoided, it is now possible to continue the 

original comparison between the particle method of CTH with MPM and continuum CTH. Even 

though 8 tracers were placed along the left plate within the model, only tracer 4 will be used for 

data comparison within this report for simplicity. This tracer is located near the middle of the left 

plate’s horizontal width, and it is indicated by the yellow circle within the key of the following

figures. The trends discussed for this tracer are assumed to be indicative of the activity at the 

other tracers.
5

The left graph in Fig. 3 depicts the pressure impact response at tracer 4 for the coarse mesh 

compared with the continuum model. For this calculation the response of the particle method is 

well below that of the continuum solution with an approximate 17% difference in peak pressure 

value. It is also important to note that the continuum solution’s equilibrium pressure value is 

approximately equal to the analytical solution, previously calculated as 2.08 × 104 MPa. The fine 

mesh pressure result comparison shown in the right graph of Fig. 3 has better agreement with the 

continuum solution with only a half percent difference in peak pressure; however, its rise time is 

much slower, leading to a pressure difference of 15% after the initial rise. Direct comparisons 

between the 2 mesh resolutions for CTH with MPM are shown in Fig. 4. This result indicates 

that the solutions from these 2 simulations are not identical. There is approximately a 14% 

difference in peak pressure values between the 2 resolution models. Additionally, neither of 

these particle results exhibits the more step-like increase, plateau, and step decrease in pressure 

obtained using continuum CTH. 

Fig. 3 (Left) Coarse mesh and (right) fine mesh continuum comparisons

Fig. 4 Mesh resolution comparison shows disagreement

3.3 Tracer Location

The results from tracer 4 at both mesh resolutions are generally below the continuum result. In 

order to investigate any influence from the boundary conditions, the group of 8 tracers was 

shifted vertically up. Now instead of taking resultant data from the lower symmetry boundary of 

the left plate, we gather it from the middle line of the left plate’s geometry. This change should 
6

not greatly affect the resultant data from the tracers, but as Fig. 5 shows, the shift causes a 

noticeable increase in impact pressure by approximately 7%. This indicates that there are

boundary condition effects on marker tracers placed on the lower symmetry boundary. This is 

not observed using continuum CTH. 

Fig. 5 Tracer 4 pressure comparison based on location

3.4 Raised Pressure 

The change in pressure output for tracer 4 as a result of its change in position requires a 

reevaluation of the original comparison between CTH with MPM and continuum CTH. This time 

we focus on the results of the fine mesh since it more closely matches the continuum results.

Figure 6 shows that the increase in pressure response obtained by tracer relocation has placed the 

particle method pressure slightly above the continuum result by approximately 4%. This near

reversal of the original pressure discrepancy can be qualified further only by investigating the 

behavior of the particle method specifically, without the coupling to the continuum grid through 

Eq. 1.

Fig. 6 Fine mesh continuum comparison after the 

change in tracer location

3.5 Uintah: A Particle Code

To gain additional insight into the performance of the CTH with MPM particle method, a second 

code is used with the current plate impact problem. Uintah is a pure particle code developed at 

the University of Utah, and it differs from CTH in that it does not factor in a continuum grid 

pressure in its momentum equation.6 The original symmetric impact across a gap problem is built 
7

within Uintah’s environment, and the results are analyzed to elucidate what is causing the 

pressure discrepancy with CTH with MPM.

First, the issue of the mixed cells halting visible plate contact was examined. This behavior can 

be explained by the method by which particle codes move bodies like the copper plates within 

the mesh. Uintah allows the plates to initially move toward each other by compressing the void 

space in between them. This method of resolving the void works well until mixed cells 

containing both particles and void are encountered just before “physical” impact. Since the 

particle code applies the same deformation rate onto copper particles and void within any cell, 

the individual particles within these mixed cells are compressed just like void space. The stress 

on these particles increases until the average stress on their cells equals the equilibrium pressure 

signaling to the code that impact has now occurred. Essentially, the code believes that the plates 

have come in contact even though there is still void space between them.

The gap issue was again sidestepped by setting the 2 plates directly next to each other for the

next Uintah simulation. Fig. 7 shows that after the impact of the 2 plates, the resulting waveform 

is not quite smooth. There is noticeable ringing occurring as the wave propagates through the 

plates. A more in-depth investigation into this code’s particle behavior is necessary. 

Fig. 7 Pressure contour from Uintah shows ringing as the impact wave travels outward 

To begin, 27 tracers are placed at consecutive particle locations in the wake of the impact wave

at 0.3 μs after impact. The pressure response of these tracers shows that after the shock passes, 2 

particles in a given cell achieve a similar stress level while the final particle falls well above or 

below in stress. It seems as though the large shock caused by the plate impact raises the stress so

quickly with respect to the mesh size that some particles comprising the plate bodies within a cell

are not able to reach the final equilibrium pressure state before being pushed into an adjoining 

cell. However, in this new cell the existing 2 particles have already achieved their final state 

either at a much higher or lower pressure. The result is a set of high and low pressures that 

combines to create the correct average pressure. As this effect propagates outward from the 

impact face, it results in the ringing seen in the pressure plot from Uintah in Fig. 7.

Placing tracers in the same locations for the same time step in CTH with MPM’s input deck 

yielded similar pressure ringing results as obtained with Uintah. The 2 sets of data were then 

plotted for graphical comparison. The scatter plot of Fig. 8 shows that ringing in CTH with MPM 

is still apparent, but it has been lessened to some degree. The smoothing of the CTH with MPM

results may be a result of the superposition of the pressure values taken from the individual 
8

particles and the background grid through Eq. 1. This mixing of particle and continuum pressure 

is not done in Uintah. Although the comparison of Uintah as a particle method brought about an 

unexpected result, further inspection of the particle code’s output brought more insight into what 

may be happening within the particle implementation of CTH.

Fig. 8 Individual particle pressure within consecutive cells after plate impact for Uintah 

and CTH with MPM

4. Summary and Conclusions

Several issues have been identified within the new implementation of the material point method 

for the historically continuum-based hydrocode CTH. The current MPM implementation within 

CTH exhibited the same lack of mixed cell void closure between impacted bodies as Uintah. 

Boundary condition effects on particle output were discovered upon changing tracer locations 

within the model. When the resolution of the mesh was changed from coarse to fine, the pressure 

response from CTH with MPM approached the CTH continuum solution. This mesh sensitivity 

may indicate that a high level of simulation resolution is necessary with this particle 

implementation in order to achieve the stress field resolution given by the traditional continuum 

method. Also, particles appeared to exhibit zero pressure values at the impact face as well as 

noticeably varied values within a single cell after impact. A number of these concerns revolve 

around the addition of variables evaluated at the particle and background grid levels within CTH 

with MPM’s solution algorithm. Further research into this implementation will be required to

gain a fundamental understanding of how this averaging method couples to particle models. 

Continued development of CTH with MPM will hopefully solve the problems encountered 

during the application of this 2-D plate impact problem and eventually much more challenging 

mechanics scenarios.
9

5. References

1. McGlaun JM, Thompson SL, Elrick MG. CTH: a three-dimensional shock wave physics 

code. International Journal of Impact Engineering. 1990;10:351–360.

2. Crawford DA, Brundage AL, Harstad EN, Ruggirello KP, Schmitt RG, Schumacher SC, 

Simmons JS. CTH user’s manual and input instructions. Ver. 10.3. Albuquerque (NM): 

Sandia National Laboratories; 2013.

3. Chen Z, Brannon R. An evaluation of the material point method. Albuquerque (NM): Sandia 

National Laboratories; 2002.

4. Schumacher SC, Rugirello KP. Marker methods. Albuquerque (NM): Sandia National 

Laboratories; 2014. 

5. Ahrens TJ. High pressure shock compression of solids. New York (NY): Springer-Verlag;

1993. p. 75–114.

6. Guilkey J, Harman T, Luijens J, Schmidt J, Thornock J, Davidson de St. Germain J, 

Shankar S, Peterson J, Brownlee C, Charles R. Uintah user guide. Ver. 1.5. Salt Lake City 

(UT): University of Utah; 1997.
10

1 DEFENSE TECHNICAL

(PDF) INFORMATION CTR

DTIC OCA

2 DIRECTOR

(PDF) US ARMY RESEARCH LAB

RDRL CIO LL

IMAL HRA MAIL & RECORDS MGMT

1 GOVT PRINTG OFC

(PDF) A MALHOTRA

1 HOWARD UNIV

 (PDF) D BOLLING

1 DIR USARL

(PDF) RDRL WMP A

S BILYK

RDRL WMP C

R BECKER

T BJERKE

RDRL WMP E

D HORNBAKER