A comparison of the shock response of the Material Point Method

The Lagrangian Material Point Method (MPM) [1, 2] has been implemented into the Eulerian shock physics code CTH [3] at Sandia National Laboratories. Eulerian hydrodynamic methods are useful for large deformation problems, where mesh tangling typically leads to difficulties for Lagrangian finite element methods. However, Eulerian techniques suffer from numerical diffusion due to advection, which can be problematic for many material models requiring the transport of a damage parameter or other state variables that need to remain sharp [4]. The inclusion of the MPM in CTH allows for the accurate simulation of structural response to shock loading in a single framework. This paper presents a comparison of the shock response of the MPM and CPDI to the CTH hydrodynamics code.


Introduction
Eulerian hydrodynamic methods are useful for large deformation problems, where mesh tangling typically leads to difficulties for Finite Element methods. However, Eulerian techniques suffer from diffusion due to numerical advection errors, which can be problematic for many material models requiring the transport of history variables that should remain sharp [4]. The Material Point Method (MPM) alleviates the issue of numerical diffusion in Eulerian methods, and also does not suffer from mesh tangling during large deformation and fragmentation.
The original MPM has been shown to suffer from cell crossing noise due to the lack of smoothness of the shape functions. To alleviate the cell crossing noise Bardenhagen et al. [5] introduced the Generalized Interpolation Material Point (GIMP) method, which allows the particle domains to deform aligned with the background grid, as shown in figure 1. However the GIMP method does not capture the rotation of the particle domains. Sadeghirad et al. [6] introduced the Convected Particle Domain Interpolation (CPDI) method which allows the particle domains to deform as parallelograms and constructs smooth basis functions on the grid, as shown in figure 2. The CPDI method has been shown to simulate large deformation problems accurately and alleviates the cell crossing noise present in the standard MPM.
The MPM and CPDI methods has been implemented into the Eulerian shock physics code CTH to obtain improved response of structural materials. This paper will examine the response of the MPM under shock loading, and compare it to the CPDI method, and the CTH response. The effect of varying the particles per cell (PPC) on the solution is also examined for both the MPM and CPDI methods. This report will be organized as follows. A brief overview of the numerical method will be provided, followed by the problem description. Next the problem results will be presented, and conclusions drawn from the study.

Numerical Method
CTH solves the conservation equations using a finite volume discretization through a Lagrangian step where the mesh is allowed to move, followed by a remap step where the material properties are remapped back to the original grid. Artificial viscosity is used to maintain monotonicity of the solution near discontinuties. Additional details on CTH can be found in [3].
The MPM uses local finite element shape functions to map the material properties from the particles to the computational grid where pressure and stress acceleration are computed. The equation of state is updated on the grid, and the updated state is then interpolated to the particles. The material strength model is evaluated on the particles. For additional details on the MPMs implementation into CTH see [2].

Problem Setup
The one dimensional problem setup for examining the shock response of the methods is shown in figure 3. The problem consists of a 3 cm bar of copper with a velocity discontinuity of 1 km/s at x=1.5 cm. The bar is assumed to be perfectly elastic for most of the cases with a poisson ratio of 0.27. All simulations are conducted in one dimension for simplicity. The CTH solution is determined to give mesh independent results with 500 cells, which is used for all cases examined. One dimensional problem setup for numerical tests.

Results
The effect of varying the PPC on the solution is examined for MPM, as shown in figure 4. The MPM solution for 2 PPC is excessively noisy and the shock speed is incorrectly predicted. Moving to 8 PPC reduces the noise substantially, but the predicted shock speed does not seem to be converging to the CTH solution. The PPC was further increased to 100 and it was found to not improve the solution over the 8 PPC case. For the CPDI method the effect of varying the PPC is shown in figure 5. The 2 and 8 PPC cases lie directly on top of each other. The solution is monotonic except for the velocity solution near x=0. The velocity oscillation is likely due to an initial start up error. The shock speed is correctly predicted. Compared to the CTH solution the CPDI shock profile (∼5-6 cells) is approximately twice the width of the CTH shock profile (∼3-4 cells). Figure 6 compares the MPM, and CPDI solutions using 8 PPC to the CTH solution on the same plot.  The velocity oscillation present in the CPDI solution is not observed to decrease with time. The reason for this can be understood by looking at the background grid data along with the particle velocities as shown in figure 7. When the particle velocities are weighted to the background grid the higher velocity particles cancel out the lower velocity ones resulting in a smooth velocity profile on the background grid. Since the gradients are computed on the background grid data, the gradients are zero and thus the particle velocities do not change. The particle positions are updated using the interpolated smooth grid velocity so the faster particles do not overtake the slower particles in front of them.
Due to the perfectly elastic assumption, pressures of approximately 250 Gigadynes/cm 2 are generated in the solid. In reality almost all materials would undergo a nonlinear deformation process under these conditions. Figure 8 shows the MPM velocity and pressure solutions when an elastic perfectly plastic material model is used for the copper with a yield strength of 70 MPa. By including plasticity the noise in the MPM solution is reduced, and the shock profile is approximately 1 cell wider than the CTH shock profile with the correct shock speed. Excessive noise in the velocity solution is still present which is likely due to initial start up errors as seen for the CPDI solution.

Conclusions
The MPM and CPDI particle methods have been compared to the CTH solution under strong shock loading. The effect of varying PPC on the methods was examined. The MPM method was shown to be noisy and give an incorrect shock speed for the perfectly elastic problem, but the noise was substantially reduced and shock speed correctly predicted by incorporating plasticity. The CPDI method gave a monotonic solution with the correct shock speed. The velocity oscillation exhibited with the CPDI method was most likely due to initial start up errors which was not observed to decrease with time due to the particle weighting to the background grid which results in a smooth grid velocity.