A Quantitative Approach to Comparing High Velocity Impact Experiments and Simulations Using XCT Data

Available online at www.sciencedirect.com

ScienceDirect Procedia Engineering 103 (2015) 610 – 617

The 13th Hypervelocity Impact Symposium

A Quantitative Approach to Comparing High Velocity Impact Experiments and Simulations Using XCT Data Andrew L. Tongea*, Brian Leavya, Jerry LaSalviaa, K.T. Rameshb, Rebecca Brannonc b

a US Army Research Lab, ATTN: RDRL-WMP-C, BLDG 390 RM 129, APG, MD 21005-5066, [email protected] Hopkins Extreme Materials Institute, The Johns Hopkins University, Hackerman Hall, Room 200, 3400 North Charles Street, Baltimore, MD 21218-2608, [email protected] c University of Utah Mechanical Engineering, 50 S. Central Campus Drive, 2254 MEB, Salt Lake City, UT 84112, [email protected]

Abstract While computational models of impact events have the potential to accelerate the design cycle, one's confidence in a material model should be related to the extent of validation work that has been performed for that model. Quantities of interest used for validation are often either scalar volume-averaged quantities (such as the average density, or the force applied to a boundary) or field quantities (such as the strain field obtained from digital image correlation, or density maps computed from X-ray computed tomography (XCT)). Volume averaged quantities are easy to compare quantitatively since they are either a single value or a simple time series. However, these averaged quantities do not capture differences in the failure process within a material and can be blunt instruments for validation efforts. Field quantities provide spatial information, but are difficult to reduce to a scalar that quantifies the goodness of a particular model with respect to another model. This work describes an approach to using XCT data to quantify how well a particular simulation agrees with simulation data while accounting for the statistical nature of failure in brittle materials. © 2015 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license © 2015 Published by Elsevier Ltd. Selection and/or peer-review under responsibility of the Hypervelocity Impact Society. (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the Curators of the University of Missouri On behalf of the Missouri University of Science and Technology Keywords: Ceramic, Boron Carbide, Balistic Impact, Validation

1. Introduction Computational modeling of impact events has great potential to reduce the time and cost required to bring a new product from the concept stage to the production stage. However, to realize this potential, computational models must undergo validation testing to build confidence that they reproduce observed physical behavior. This process is particularly difficult for those cases where stochastic behavior leads to sampling and variability issues. As an example, an important physical feature of ceramic materials used for impact applications is the variability in their response. Without accounting for the inherent variability in ceramic materials, it is doubtful that a simulation can predict the distribution of possible responses necessary for uncertainty quantification and risk assessment activities [1]. In the presence of variability, a particular simulation should not be expected to match a particular experiment at each spatial location; both the simulations and the experiments should produce a distribution of responses. The high computational cost associated with determining the converged distribution of responses for a particular quantity of interest using an approach similar to that of [2] is justified for some applications, but for other applications the analyst only needs an estimate of the “expected response.” Even with the lower fidelity requirements for the latter case, one needs an approach for defining terms such as “more similar to experiments” [3] when faced with a choice between many different models (or the challenge of selecting material parameters for a particular model). Defining these terms necessitates an objective, quantitative, and mathematically rigorous basis for quantifying similarity between experiments and simulations.

1877-7058 © 2015 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license

(http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the Curators of the University of Missouri On behalf of the Missouri University of Science and Technology

doi:10.1016/j.proeng.2015.04.079

Andrew L. Tonge et al. / Procedia Engineering 103 (2015) 610 – 617

611

Advanced X-ray computed tomography (XCT) has made it possible to reconstruct three-dimensional density profiles within a specimen [4]. Since simulations of impact events can also produce post-impact density maps within a specimen, comparing these two outputs is a natural choice for model validation efforts. XCT imaging of brittle failure has long been recognized for its model-validation potential (cf. [3,5,6]), but it is typically used in only a qualitative manner, where plots of measured and computed damage fields might be displayed side-by-side without any quantitative measure of similarity. Perhaps the greatest focus to date on 2D comparisons of binary fields (which, for our purposes, would be damaged and undamaged zones) has been confined to image-analysis literature (cf. [7]) where the goal has been to define a mapping that will transform a test image into a reference image for the purposes of segmentation or tagging. There has been relatively little work on quantitative 3D comparisons except in cases of target fragmentation, in which quantitative comparisons of fragment size statistics are relatively straightforward (cf. [8,9] and references therein). As discussed in the verification, validation, and uncertainty quantification literature (cf. [10]), one's confidence in model predictions should not be any greater than one’s confidence in the associated validation efforts. In this work, we develop a quantitative measure of the similarity between the results of impact simulations and experiments. We propose a comparison metric that can quantitatively rank the results from any two simulations for their similarity to experiments. The comparison metric is demonstrated using impact velocities of 300 m/s, but the technique applies to the postmortem analysis of any impact event, including hypervelocity impacts. 1.1. Example material model In this work, we apply the Tonge-Ramesh material model [11] to boron carbide, an advanced structural ceramic. The material model accounts for the interaction and growth of a distribution of flaws through a micromechanics-based damage model. Once the damage level reaches a critical value, granular flow is activated, and the material is treated as an elastic viscoplastic material with a Drucker-Prager flow surface and associative flow behavior. These constitutive features introduce porosity as the material is sheared. The dilatation produced through granular flow is used to compute the distension in a P−α type porosity model [12] to allow inelastic re-compression as needed. Similar to work by Hild et al. [13], the material model explicitly introduces variability into the constitutive response by using a local Poisson process to create local realizations of the flaw distribution. The material parameters used for these simulations are the same as the ones used in [14, ch. 5]. They are listed in table 1 for completeness. These material parameters were taken from the literature when possible. The shape of the flaw size distribution is based on experimental observations of the distribution of carbon flakes by J. Hogan. The flaw density was calibrated based on the median strength for dynamic uniaxial stress compression at a strain rate of 500 1/s assuming a homogeneous stress state. 1.2. Experimental Method For comparison purposes, a commercial hot-pressed boron-carbide variant manufactured by Coorstek Vista with cylindrical geometry was impacted with a cemented carbide sphere at 281 m/s and examined by X-ray Computed Tomography (XCT). The density, average grain size, Young’s modulus, Poisson’s ratio, fracture toughness, and bend strength for this boron-carbide variant were measured to be 2.52 g/cm3, 15 μm, 460 GPa, 0.17, 2.9 MPa∙m1/2, and 450 MPa, respectively. The primary minor phase was determined to be graphite. The cylinder was nominally 38.1 mm in diameter by 19.1 mm thick. The surface roughness of the impact surface was Ra value of 0.18 μm. In addition, the cylinder was encased on its sides with a cold-mount epoxy 6.35 mm thick to aide in its intact recovery. The cemented carbide sphere (tungsten carbide with 6 wt.% cobalt, Grade 25, 0.000635 mm roundness tolerance, Class C-2) was supplied by Machining Technologies Inc. The sphere was launched on a sabot using a gas-gun with a 12.7 mm bore. Further details of the boron6.34 mm

38.1 mm

19.1 mm

(a) Default particle domains

(b) Conforming particle domains

Figure 1: The default cuboidal particle domains produce stair-stepped boundaries that can contribute to mesh bias that propagates from the contact discretization. Conforming boundaries based on a hexahedral mesh can help mitigate this issue.

612

Andrew L. Tonge et al. / Procedia Engineering 103 (2015) 610 – 617

carbide variant, ballistic test set-up and procedures are described elsewhere [19,20]. Following the ballistic test, the boroncarbide cylinder was examined by XCT using a system built by North Star Imaging. The X-ray scans were conducted at 450 kV. The 2D and 3D spatial resolutions were 125 μm and 114 μm, respectively. 2. Simulation approach and setup In this work we use the Convected Particle Domain Interpolation (CPDI) technique [21] to solve the weak form of the balance of linear momentum. CPDI uses average quantities over particle domains to approximate the integrals used in the weak form to compute nodal forces. In the Uintah [22] implementation of CPDI, particle domains are parallelepipeds. For our baseline simulations, we use initially cuboidal particle domains to discretize the geometry. This leads to a poor representation of the curved cylindrical and spherical surfaces as discussed in [23]. Figure 1 shows a comparison between the default cuboidal particle geometry (figure 1a) and the same geometry discretized using particle domains that more closely conform to the boundaries of the geometry (figure 1b). The conforming particle domains do a better job of representing the spherical impactor than the cuboidal domains; however, generating the conforming particle domains adds an additional processing step of generating a hexahedral mesh, from which the particle domains are computed. In all simulations, we use a background grid spacing of 0.25 mm and particle sizes that are close to 0.125 mm on a side. The particles are not all the same size, because they were generated from a mapped hexahedral mesh. We have performed simulations using a mesh that is two times coarser (0.25 mm particles) and observed that the thickness of the radial features is tied to the mesh size. As the purpose of this work is presenting an approach for comparing simulations and experiments, we leave a discussion of the convergence or non-convergence of the particular numerical technique and material model for a separate discussion. To account for the statistical nature of the material response, we performed two simulations using different random seeds. These two simulations were then compared with each other and with the experimental results. Table 1: Material model parameters for boron carbide Material Parameter

Value

Reference

Density (ρ)

2520 kg/m3

Theoretical Density [15]

Specific Heat Capacity (Cv)

962 J/(kg K)

Dandekar [16]

Bulk Sound Speed (C0)

9.6 km/s

Calculated

Us-Up Slope (S)

0.914

Dandekar [16]

Gruneisen Parameter (Γ0)

1.28

Dandekar [16]

Shear Modulus (G0)

197 GPa

Paliwal and Ramesh [17]

Bulk Modulus (K0)

232 GPa

Dandekar [16]

Minimum Flaw Size (smax)

1.0 μm

Maximum Flaw Size (smin)

25 μm

Distribution Exponent (α)

2.6

Flaw Density (η)

22×10121/m3

Fit to dynamic strength

Fracture Toughness (KIC)

2.5 MPa m1/2

Paliwal and Ramesh [17]

Rayleigh Wave Speed (Cr)

8.0 km/s

Calculated

Maximum Crack Velocity (vm)

0.2Cr

Experiments

Crack Growth Exponent (γc)

1.0

Crack Face Friction Coefficient (μ)

0.8

Crack Orientation (ϕ)

60o

Most Damaging

Granular Slope (A)

0.8

Fit to Chocron et al. [18]

Cohesive Strength (B)

3 MPa

Fit to Chocron et al. [18]

Relaxation time (τ)

7 ns

Activation Damage (Dc)

0.125

Maximum Damage (Dmax)

1.0

Reference Crush Pressure (P0)

100 MPa

Reference Distension (JGP0)

2.0

Consolidation pressure (Pc)

10 GPa

Andrew L. Tonge et al. / Procedia Engineering 103 (2015) 610 – 617

613

Figure 2: Illustration of the image alignment process using triangles. Images are first shifted so that the centroids coincide with the origin. The second image is then rotated to minimize the error measure.

3. A quantitative measure of the difference between two scalar fields For this work, we begin by selecting the L2 norm as a quantitative measure of the relative difference between two scalar fields, Asim and Aexp:

³ ( A A ) dV ³ ( A ) dV 2

L2 =

sim

exp 2

exp

Since this is a point-wise comparison, both the simulation and experimental data must be mapped to a common grid for the purposes of computing the integral in the above equation. When comparing two spatially varying fields, a coordinate system’s origin and axes orientations constitute a set of arbitrary degrees of freedom, constrained only by applicable symmetry groups. For example, when comparing two cubic domains of a statistically isotropic medium subjected to nonshear loading (so that a cuboid remains a cuboid), there are eight possible choices for the “top” ( ) face and four possible choices for an orientation of that face for a total of 32 possible orientations for comparing two cubes. With a cylinder there are two choices for the top and bottom of the target, but once those are chosen, the in-plane rotation is arbitrary. The coordinate system should be chosen from the coordinate systems that respect the symmetries (rotational for this problem) in

Figure 3: XCT data shown in a rendered three-dimensional view using a color scale to highlight details next to a gray-scale slice used for comparison to the simulation results.

614

Andrew L. Tonge et al. / Procedia Engineering 103 (2015) 610 – 617

the loading to minimize the difference between the experimental observations and the simulations. With the goal of minimizing the number of tuning parameters that are available for adjusting the level of agreement between the simulations and experiments, we define a consistent approach for choosing the reference grid (or coordinate system). We are interested in comparing simulations and experiments where a spherical projectile impacts a cylindrical target. We select the direction of the z-coordinate axis so the projectile traveled in the negative z direction prior to impact. We then select the center of mass of the image intensity (I(x)) as the origin of the coordinate system. The center of mass is given by: x = ³ xI(x)dA / ³ I(x)dA We choose the rotation angle about the z-axis that minimizes L2 where the quantity of interest A(x) is the image intensity. This process is shown in figure 2 by aligning two equilateral triangles that have been initially offset and rotated. 4. Quantity for comparison Simulations and experiments should be compared using equivalent quantities that are directly measured in both the experiment and the simulation. XCT assigns a gray scale value to a voxel based on the average density within the volume; therefore, in principle, one should directly compare the measured density from the XCT scan with the density field at the end of the simulation. However, we do not have a calibration table that relates a gray-scale value to a density, and therefore, we choose to compare quantities that are proportional to the volume fraction of void space that has been introduced into the material. For the experiments, we assume that dark regions within the boron-carbide cylinder correspond to areas of reduced density, which are the result of voids formed by crack growth and sliding (including deformation of the comminuted material). Figure 3 shows a rendering of the three-dimensional XCT data and a two-dimensional slice taken from near the impact surface. The color scale, shading, and transparency in the three-dimensional image were chosen to highlight the features in the XCT image while the gray-scale slice does not have similar image processing applied. Within the Tonge-Ramesh material model, granular flow (relative sliding and rolling of sub-grid scale material fragments) produces dilatation, which is then tracked using an internal variable denoted JGP. This quantity is the ratio of the density of the solid material supporting the applied load (ρe) to the average density in the representative volume element including the void space (ρ). It is analogous to the distension in [11]; however, the two quantities differ by the volume change ratio in the plastically incompressible matrix material, which is typically small compared to JGP. When there are no residual stresses and no thermal expansion, ρe=ρ0 (where ρ0 is the initial density), and JGP=ρ0/ρ. This indicates that it is reasonable to compare JGP to the XCT results. After identifying the quantities for comparison, one must define a method for converting these two different measurements into a common system of units. For simplicity we perform the comparison in a binary fashion. Physically, we are attempting to identify regions that ether have or have not undergone dilatation as a result of the impact event. For the experimental results, we first define a region of interest that separates the boron-carbide cylinder from the rest of the image. This selection is done manually to reduce the complexity of the process, but it could be done using existing image segmentation algorithms. After identifying a region of interest, we apply a threshold to the gray-scale image. The threshold represents an additional input parameter, and its effect on the results will be discussed in section 6. To obtain a similar

Figure 4: Image processing procedure for both simulation images (top) and the experimental images (bottom).

Andrew L. Tonge et al. / Procedia Engineering 103 (2015) 610 – 617

binary representation from the simulation results, we first use the visualization tool VisIt [24], to create a typical contour plot of the field variable JGP. For simplicity, we use a gray-scale color scheme with white indicating regions with JGP greater than or equal to 1.2 and black indicating JGP equal to 1.0. After exporting the visualization image, we crop the image to isolate the boron-carbide target, fill the white background color with black, and process the image using a very small threshold so that regions that experienced bulking due to granular flow are white. These steps are shown in figure 4 for both the simulations (top) and the XCT data (bottom). Once the experimental and simulation images are converted to binary images, they can be aligned using the procedure discussed in section 3. In the next section, we discuss the results of a comparison between two simulations of nominally identical conditions where the only difference is a change in the random seed used to assign the initial flaw distribution within the material. 5. Comparison between two simulations We begin by comparing two simulation results, because comparing simulation results to each other requires less interpretation than comparing simulation results to experiments. The upper-left image in figure 5 shows the granular flow region for the first simulation in light blue, the second simulation in yellow, and overlapping granular flow regions in red after aligning the centroids prior to performing any rotations. The upper-right image shows the same comparison after rotating the results from the second simulation 9 degrees to minimize L2. There are regions that agree and regions that disagree because the two simulations represent different specimens prepared from the same material. The simulations agree on the approximate size of the central granular flow region and some of the radial features agree, but the radial features do not all line up exactly. This lack of agreement along with a rotation angle that does not correspond to the symmetry of the background mesh provides additional confidence that the observed features represent physical stochastic processes. If a sufficiently large number of simulations were performed to compute the distribution of rotation angles (e.g. [2]), then one should expect a uniform distribution of rotation angles. Any deviation from the physically expected distribution is an undesirable expression of a bias introduced by the discretization. The radial features emanating from the central damage region represent a loss of radial symmetry in the problem. In brittle materials, radial symmetry is lost when cracks begin to grow. This occurs when sufficient energy can be dissipated by a few discrete cracks instead of a uniform damage front. The frequency of radial cracking is related to both the rate sensitivity of the failure process and to the variability of the material [25]. 6. Comparison between simulations and experiments Comparing simulations to the experimental results introduces an additional parameter, which is the threshold intensity

Figure 5: Image comparison results for two simulations that use a different random seed but are otherwise identical.

615

616

Andrew L. Tonge et al. / Procedia Engineering 103 (2015) 610 – 617

used to distinguish the intact, undamaged boron carbide from the damaged regions (in the simulations). To illustrate the effect of this parameter, we perform the comparison using two different values for the threshold. The lower threshold value of 0.259 (on the left) retains most of the linear crack-like features in the image, but also includes some regions that seem to be noise. The higher threshold value of 0.294 (on the right) loses some of the linear features, but contains fewer noise regions. Figure 6 shows the comparison between the simulation and the experiment using these two threshold values. With the lower threshold value, the error measure L2=1.31×10-3, while the higher threshold results in an error measure of L2=1.28×10-3. Although there is clearly a visible effect of the threshold value on the region that is marked as damaged, the quantitative error is only moderate. Both threshold values provide enough information to align the experimental and simulation images and perform the comparison. In practice, experimental measures of damage also involve the use of some threshold, which sometimes is the result of an instrument limit such as resolution, and sometimes is the result of operator choices during the imaging process. Further, there is also substantial variation from experiment to experiment, so that eventually one may need to compare within a suite of experiments to determine which observed features are robust with respect to the material and the loading. The radial damaged regions in the experiments are thinner than similar regions in the simulations. Additionally, the central region inside of the radial features is smaller in the experiments than it is in the simulations. This may indicate that the material model used for these simulations does not localize as readily as the real material or that the computational framework used for the computations does not support the localization. A comparison between even simplified ballistic impact problems is sufficiently complicated that such ballistic impact makes for a much better validation problem than a calibration problem. Ideally, models should be calibrated to experiments that are able to isolate as many mechanisms as possible, while validation experiments strive to exercise the interactions between the different mechanisms. 7. Summary In this work, we have presented an approach for comparing simulations with experimental of impact events on brittle targets. By quantifying the difference between simulations and experiments using the L2 error norm and the alignment procedure discussed in section 3, we have defined a scalar that represents the quality of the agreement. This scalar can be used to support claims that simulation set A is more similar to experiments than simulation set B, which will facilitate the development of high fidelity material models for the failure of brittle materials. The approach consists of first processing images from either the simulations or the experiments so that they are binary representations of where a density change associated with material failure has occurred. These images are then aligned first by aligning the centroids of the damaged regions and then by applying a rotation about the centroid to minimize the error measure.

Figure 6: Comparison between the simulation results and the experimental results for two different threshold values.

Andrew L. Tonge et al. / Procedia Engineering 103 (2015) 610 – 617

617

Acknowledgements The efforts of A.L. Tonge, K.T. Ramesh, and R.M. Brannon were partially supported by the Army Research Laboratory under Cooperative Agreement Number W911NF-12-2-0022. In addition to support from CRA W911NF-12-2-0022, A.L. Tonge acknowledges support through the Oak Ridge Institute for Science and Education postdoctoral fellowship program at the Army Research Laboratory. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory or the U.S. government. The U.S. Government is authorized to reproduce and distribute reprints for Governments purposes notwithstanding any copyright notation herein. References [1] Strack, O.E. and Leavy, R.B. and Brannon, R.M. 2014. Aleatory uncertainty and scale effects in computational damage models for failure and fragmentation, International Journal for Numerical Methods in Engineering, doi: 10.1002/nme.4699 [2] Bishop, J. E. and Strack, O. E. 2011. A statistical method for verifying mesh convergence in monte carlo simulations with application to fragmentation, International Journal for Numerical Methods in Engineering, 88, p. 279–306. [3] Brannon, R. M., Wells, J. M., and Strack, O. E. 2007. Validating theories for brittle damage, Metallurgical and Materials Transactions A: Physical Metallurgy and Materials Science, 38, p. 2861–2868. [4] Wells, J. M. and Brannon, R. M. 2007. Advances in x-ray computed tomography diagnostics of ballistic impact damage. Metallurgical and Materials Transactions A: Physical Metallurgy and Materials Science, 38, p. 2943–2949. [5] Leavy, R. B., Clayton, J. D., Strack, O. E., Brannon, R. M., and Strassburger, E. 2013. Edge on impact simulations and experiments. In Hyper Velocity Impact Society Symposium Procedings, 58, p. 445–452. [6] O’Masta, M. R., Deshpande, V.S., and Wadley, H.N.G. 2014. Mechanisms of projectile penetration in dyneema encapsulated aluminum structures. International Journal of Impact Engineering, 74, p. 16-35. [7] Lee, E. and Gunzburger, M. 2010. An optimal control formulation of an image registration problem. Journal of Mathematical Imaging and Vision. 36, p. 69–80. [8] Meyer, H. W., Jr. and Brannon, R. M. 2012. A model for statistical variation of fracture properties in a continuum mechanics code. Intternational Journal of Impact Engineering, 42, p. 48-58. [9] Zhou, F., Molinari, J., and Ramesh, K. T. 2006. Characteristic fragment size distributions in dynamic fragmentation. Applied Physics Letters, 88. [10] Schwer, L.E. 2007. An overview of the PTC 60/V&V 10: Guide for verification and validation in computational solid mechanics: Transmitted by L. E. Schwer, chair PTC 60/V&V 10, Engineering with Computers, 23, p. 245–252. [11] Tonge, A. L. and Ramesh, K. T. 2015. Multi-scale defect interactions in high rate brittle material failure, Part 1: Model formulation and application to AlON. Journal of the Mechanics and Physics of Solids, In Process. [12] Carroll, M. and Holt, A.C. 1972. Suggested modification of the P-α model for porous materials. Journal of Applied Physics, 43, p. 759–761. [13] Hild, F., Denoual, C., Forquin P., and Brajer, X. 2003. On the probabilistic–deterministic transition involved in a fragmentation process of brittle materials, Computers & Structures: Advanced Computational Models and Techniques in Dynamics, 81, p. 1241–1253. [14] Tonge, A.L. 2014. A unified framework which uses multi-scale microstructural information for modeling dynamic failure in brittle materials. PhD Thesis, The Johns Hopkins University. [15] Thevenot, F. 1990. Boron carbide-A comprehensive review, Journal of the European Ceramic Society, 6, p. 205–225. [16] Dandekar, D. P. Shock response of boron carbide: Technical Report ARL-TR-2456, Army Research Laboratory, April 2001. [17] Paliwal, B. and Ramesh, K. T. 2007. Effect of crack growth dynamics on the rate-sensitive behavior of hot-pressed boron carbide. Scripta Materialia, 57, p. 481–484. [18] Chocron, S., Anderson, C. E. Jr., Dannemann, K. A., Nicholls, A. E., and King, N. L. 2011. Intact and predamaged boron carbide strength under moderate confinement pressures, Journal of the American Ceramics Society, 95(2011) doi:10.1111/j.1551-2916.04931.x [19] Normandia, M. J. 2004. Impact response and analysis of several silicon carbides, International Journal Applied Ceramic Technology, 1, p. 226-234. [20] LaSalvia, J. C., Normandia, M. J., Miller, H. T., and MacKenzie, D. E. 2009 “Sphere impact induced damage in ceramics: II. Armor-grade B4C and WC,” Advances in Ceramic Armor: A Collection of Papers Presented at the 29th International Conference on Advanced Ceramics and Composites, January 23-28, 2005, Cocoa Beach, Florida, Ceramic Engineering and Science Proceedings. [21] Sadeghirad, A., Brannon, R.M., and Burghardt, J. 2000. A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations. International Journal for Numerical Methods in Engineering, 86, p. 1435–1456. [22] de St. Germain, J., McCorquodale, J., Parker, S., and Johnson, C. 2000. “Uintah: A massively parallel problem solving environment,” Ninth IEEE International Symposium on High Performance and Distributed Computing, p. 33–41. [23] Kamojjala, K., Brannon, R., Sadeghirad, A., and Guilkey, J. 2013. Verification tests in solid mechanics, Engineering with Computers, doi: 10.1007/s00366-013-0342-x. [24] Childs, H., Brugger, E., Whitlock, B., Meredith, J., Ahern, S., Pugmire, D, Biagas K., Miller, M., Harrison C., Weber, G. H., Krishnan H., Fogal, T., Sanderson, A., Garth, C., Bethel, E. W., Camp, D., Rübel, O., Durant, M., Favre, J. M., and Navrátil, P. 2012. “Visit: An end-user tool for visualizing and analyzing very large data,” In High Performance Visualization–Enabling Extreme-Scale Scientific Insight, p. 357–372. [25] Leavy, R.B. and Brannon, R.M. and Strack, O.E. 2010. The use of sphere indentation experiments to characterize ceramic damage models. International Journal of Applied Ceramic Technology, 7, p. 606-615.