# A Discontinuous Galerkin Material Point Method for the simulation of impact problems in solid mechanics

Adrien Renaud, Thomas Heuzé

► To cite this version:

Adrien Renaud, Thomas Heuzé. A Discontinuous Galerkin Material Point Method for the simulation of impact problems in solid mechanics. 6th International Conference on Computational Methods in Structural Dynamics and Earthquake Engineering Methods in Structural Dynamics and Earthquake Engineering (COMPDYN 2017), Jun 2017, Rhodes Island, Greece. 10.7712/120117.5678.17188 . hal-01913197

HAL Id: hal-01913197

https://hal.science/hal-01913197

Submitted on 6 Nov 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A DISCONTINUOUS GALERKIN MATERIAL POINT METHOD (DGMPM) FOR THE SIMULATION OF IMPACT PROBLEMS IN SOLID MECHANICS

A. Renaud¹, T. Heuzé¹

¹ Research Institute in Civil and Mechanical Engineering (GeM, UMR 6183 CNRS)
Ecole Centrale de Nantes
1 rue de la Noë, Nantes

e-mail: {adrien.renaud,thomas.heuze}@ec-nantes.fr

**Keywords:** Material Point Method, Discontinuous Galerkin, dynamic impact, elastic solid, hyperelastic solid.

**Abstract.** The material point method is extended in this work to the Discontinuous Galerkin approximation framework for the simulation of impacts on elastic and hyperelastic solids. The formulation is based on the weak form of conservation laws on each cell of an eulerian grid in which volume integrals are discretized on a set of material points lying in that cell, and on the computation of Godunov fluxes at cells faces. The resulting method is first derived within the small strains framework and illustrated on a one-dimensional and a two-dimensional problem of impact on an elastic media. Then a one-dimensional hyperelastic problem of a solid undergoing large strains is presented, and a comparison is performed with an analytical solution.
# 1 INTRODUCTION

The numerical simulation of hyperbolic initial boundary value problems (IBVP) in solid mechanics has been and is again mainly performed with the finite element method (FEM). Despite its well-known advantages, this method becomes less efficient and accurate for problems involving large strains because of mesh tangling and hence, the required re-meshing and projection steps. Meshless methods have been developed to circumvent these shortcomings. Among them, the material point method [1] (MPM) uses a collection of particles that can move in an eulerian grid to represent a continuum with a lagrangian description. The material points carry all the fields of the problem while the grid is only used to compute gradients and solve the discrete equations. Fields travel from material points to nodes and *vice versa*, and since history is known at every single particle, the grid can be discarded and reconstructed for convenience.

Hyperbolic IBVP including extreme loading conditions such as impacts require the ability to accurately track the wave front of shock waves to provide a better understanding of physical phenomena occurring in dynamic forming processes such as electromagnetic forming. However, classical time integrators used with the FEM introduce high frequency noise in the vicinity of discontinuities which is hard to remove without loss of accuracy. Furthermore, the MPM space discretization leads to a decrease of the Courant number that prevents the capture of discontinuities. The discontinuous galerkin finite element method (DGFEM) has been developed for the neutron transport equation [2] and is now used in fluid and solid mechanics [3]. This extension of classical finite element method is based on the discontinuous Galerkin (DG) approximation in which shape functions are discontinuous across elements boundaries. The continuity of fields is not enforced on the whole mesh but only element-by-element and it allows to capture sharp solutions like shocks. Moreover, the characteristic structure of hyperbolic problems can be introduced within the numerical scheme as what has been first done with the finite volume method [4]. Hence, the DG approximation is well-suited to accurately capture wave fronts but is still limited by mesh tangling in a lagrangian approach where large deformations occur.

The purpose of this work is to develop a numerical method which allows to accurately follow wave fronts in a solid that is subjected to large deformations in a lagrangian framework. The eulerian grid of MPM can be used as a support for the DG approximation and the DGMPM introduced here aims at meeting advantages of both MPM and DG methods while freeing of their respective limits. We first derive the DGMPM formulation within the large strain framework with a total lagrangian formulation which is then linearized to apply to small strains framework. Then, we focus on the computation of interface fluxes resulting from the weak form. The method is finally illustrated on one-dimensional and two-dimensional problems, and comparisons are performed with the material point method, the finite element method and the analytical solution when it is known.

# 2 THE HYPERELASTIC DYNAMIC CONTINUUM PROBLEM

Let $\Omega_0 \in \mathbb{R}^3$ denote the reference position of a continuum body, with reference coordinates system $\mathbf{X}$, and $\Omega_t \in \mathbb{R}^3$ the position of the same body at a different time with the current coordinates system $\mathbf{x}$. The motion of the body is described by the smooth function $\varphi(\mathbf{X}, t)$ which gives the position $\mathbf{x}$ at time $t$ of a particle initially located at $\mathbf{X}$. The mapping between current and reference configurations is $\mathbf{x} \equiv \varphi(\mathbf{X}, t)$. The deformation gradient tensor and the
velocity vector are defined as :

$$F(\mathbf{X}, t) := \frac{\partial \varphi(\mathbf{X}, t)}{\partial \mathbf{X}} \qquad (1)$$

$$v(\mathbf{X}, t) := \dot{\varphi}(\mathbf{X}, t) \qquad (2)$$

where the superposed dot denotes the material time derivative. Combination of equations (1) and (2) yields :

$$\dot{F}(\mathbf{X}, t) = \frac{\partial v(\mathbf{X}, t)}{\partial \mathbf{X}} \qquad (3)$$

The solid must satisfy locally the conservation of linear momentum written in the reference configuration within the time interval of interest $t \in ]0, T]$, and with neglected body forces :

$$\frac{\partial \rho_0 v}{\partial t} - \nabla_X \cdot \Pi = 0 \quad \forall \mathbf{X} \in \Omega_0 \qquad (4)$$

where $\rho_0(\mathbf{X}) = \rho(\mathbf{x}, t = 0)$ is the reference mass density, $\Pi$ is the first Piola-Kirchhoff stress tensor and $\nabla_X \cdot (\bullet)$ is the divergence vector with respect to reference coordinates. The Plohr and Sharp formulation [5] consists in viewing equation (3) as a conservation law in order to write the following homogeneous system of partial differential equations :

$$\frac{\partial q}{\partial t} + \sum_{l=1}^{3} \frac{\partial \mathcal{F}_l}{\partial X_l} = 0 \qquad (5)$$

where $q$ is the vector of conserved quantities and $\mathcal{F}_l$ the $l^{th}$ flux vector :

$$q = \begin{cases} F \cdot e_1 \\ F \cdot e_2 \\ F \cdot e_3 \\ \rho_0 v \end{cases} \quad ; \quad \mathcal{F}_l = - \begin{cases} (v \otimes e_1) \cdot e_l \\ (v \otimes e_2) \cdot e_l \\ (v \otimes e_3) \cdot e_l \\ \Pi \cdot e_l \end{cases} = \mathcal{F} \cdot e_l \qquad (6)$$

In the particular case of linear elasticity within the small strains framework, the system is built by replacing equation (3) with the geometrical compatibility equations $\dot{\varepsilon} = \frac{1}{2}(\nabla v + \nabla v^T)$, that can also be combined with the elastic constitutive equation $\sigma = \mathbb{H} : \varepsilon$, where $\mathbb{H}$ is the fourth-order elasticity tensor. The vectors expressions are :

$$q = \begin{cases} \sigma \cdot e_1 \\ \sigma \cdot e_2 \\ \sigma \cdot e_3 \\ \rho v \end{cases} \quad ; \quad \mathcal{F}_l = - \begin{cases} \mu (v_1 e_l + v_i e_i) + \lambda v_l e_1 \\ \mu (v_2 e_l + v_i e_i) + \lambda v_l e_2 \\ \mu (v_3 e_l + v_i e_i) + \lambda v_l e_3 \\ \sigma \cdot e_l \end{cases} \qquad (7)$$

with $\sigma$ being the Cauchy stress tensor, $\varepsilon$ the linearized strain tensor, and $(\lambda, \mu)$ are Lamé's coefficients.

### 3 THE DISCONTINUOUS GALERKIN MATERIAL POINT METHOD

A continuum body $\Omega_t$ is described by a set of $P$ material points in an eulerian grid composed of $N$ quadrangular 4-nodes elements (or cells) such that the volume covered by the grid is $V = \cup_{i=1}^N V_i$ with $V_i \cap V_j = \emptyset \forall i, j$ (fig. 1).
<figure><img src="image_2.png" /><figcaption>Figure 1: Representation of a continuum body by a set of material points in $\mathbb{R}^3$.</figcaption></figure>

## 3.1 Weak form of equation of motion

The discontinuous galerkin approach relies on approximation functions that are discontinuous at cells interfaces and hence, the weak form of equation (5) is written on each element of the grid [3]:

$$ \int_{\mathcal{V}_e} \frac{\partial \boldsymbol{q}}{\partial t} \cdot \boldsymbol{w} \, dV + \int_{\mathcal{V}_e} \sum_l \frac{\partial \mathcal{F}_l}{\partial X_l} \cdot \boldsymbol{w} \, dV = 0 \quad \forall \boldsymbol{w}, \forall \mathcal{V}_e \quad (8) $$

Integration by parts leads to:

$$ \int_{\mathcal{V}_e} \frac{\partial \boldsymbol{q}}{\partial t} \cdot \boldsymbol{w} \, dV - \int_{\mathcal{V}_e} \sum_l \mathcal{F}_l \cdot \frac{\partial \boldsymbol{w}}{\partial X_l} \, dV + \int_{\Gamma_e} (\mathcal{F} \cdot \boldsymbol{n}) \cdot \boldsymbol{w} \, d\Gamma = 0 \quad \forall \boldsymbol{w}, \forall \mathcal{V}_e \quad (9) $$

where $\boldsymbol{n}$ is the outward normal vector to the boundary and $\mathcal{F} \cdot \boldsymbol{n} = \mathcal{F}_n$ denotes a normal flux at cell faces (which computation is developed in section 4). We introduce specific quantities that will be injected in the weak form (9):

$$ \boldsymbol{q}^s = \rho_0 \boldsymbol{q} \quad ; \quad \mathcal{F}_k^s = \rho_0 \mathcal{F}_k \quad (10) $$

Moreover, a particular representation of mass density based on material points mass is used:

$$ \rho(\boldsymbol{x}, t) = \sum_{\alpha=1}^{P} m_{\alpha} \delta(\boldsymbol{x}_{\alpha} - \boldsymbol{x}) \implies \rho_0(\boldsymbol{X}) = \rho(\boldsymbol{x}, 0) = \sum_{\alpha=1}^{P} m_{\alpha} \delta(\boldsymbol{X}_{\alpha} - \boldsymbol{X}) \quad (11) $$

where $\delta$ is the delta dirac distribution and $m_\alpha$ is the mass of the $\alpha^{th}$ material point. Such a discretization of the reference mass density combined with the writing of conservation laws (5) leads to a total lagrangian formulation. Equation (9) thus reads:

$$ \int_{\mathcal{V}_e} \rho_0 \frac{\partial \boldsymbol{q}^s}{\partial t} \cdot \boldsymbol{w} \, dV - \int_{\mathcal{V}_e} \rho_0 \sum_l \mathcal{F}_l^s \cdot \frac{\partial \boldsymbol{w}}{\partial X_l} \, dV + \int_{\Gamma_e} \mathcal{F}_n \cdot \boldsymbol{w} \, d\Gamma = 0 \quad \forall \boldsymbol{w}, \forall \mathcal{V}_e \quad (12) $$

which, with delta dirac integration properties transforms to:

$$ \sum_{\alpha=1}^{P} \left( m_{\alpha} \frac{\partial \boldsymbol{q}^s(\boldsymbol{X}_{\alpha})}{\partial t} \cdot \boldsymbol{w}(\boldsymbol{X}_{\alpha}) - \sum_{l} \mathcal{F}_{l}^{s}(\boldsymbol{X}_{\alpha}) \cdot \frac{\partial \boldsymbol{w}(\boldsymbol{X}_{\alpha})}{\partial X_{l}} \right) + \int_{\Gamma_e} \mathcal{F}_n \cdot \boldsymbol{w} \, d\Gamma = 0 \quad \forall \boldsymbol{w}, \forall \mathcal{V}_e \quad (13) $$
## 3.2 Nodes to material points mapping

As in MPM we need to transfer fields from nodes to material points and we use shape functions which support reduces to only one cell :

$$ \mathbf{q}^s(\mathbf{X}_\alpha) = \sum_{i=1}^{N} S_{i\alpha} \mathbf{q}_i^s \qquad (14) $$

where $S_{i\alpha}$ denotes the shape functions attached to the $i^{th}$ node evaluated at $\mathbf{X} = \mathbf{X}_\alpha$, and $\mathbf{q}_i^s$ the specific vector of conserved quantities at node $i$. Once introduced in (13), it comes :

$$ w_i \cdot \sum_{\alpha=1}^{P} S_{i\alpha} m_\alpha S_{j\alpha} \frac{\partial \mathbf{q}_j^s}{\partial t} - w_i \cdot \sum_l \frac{\partial S_{i\alpha}}{\partial X_l} m_\alpha S_{i\alpha} \mathbf{F}_{l,j}^s + w_i \cdot \int_{\Gamma_e} S_i(\mathbf{X}) \mathbf{F}_n d\Gamma = 0 \quad \forall w_i \quad (15) $$

$$ \Leftrightarrow M_{ij} \frac{\partial \mathbf{q}_j^s}{\partial t} - \sum_{l=1}^{3} K_{ij}^l \mathbf{F}_{l,j}^s + \mathbf{F}_j^* = \mathbf{0} \qquad (16) $$

which is the semi-discrete system that must be solved on the grid. Time discretization is performed with an explicit algorithm (forward Euler, Runge-Kutta 2), leading to discrete equations.

**Remark :** In particular configurations of material points, the consistent mass matrix $M_{ij}$ can be singular [6] and it is avoided by lumping this matrix : $M_i^L = \sum_j M_{ij}$.

## 3.3 Material points to nodes mapping

In order to solve system (16), fields have to be defined at nodes and an interpolation of conserved quantity is performed :

$$ \mathbf{q}_i = M_i^L \mathbf{q}_i^s = \sum_{\alpha=1}^{P} S_{i\alpha} m_\alpha \mathbf{q}_\alpha^s \qquad (17) $$

to be solved for each $\mathbf{q}_i^s$. This mapping is a mass weighted interpolation that stabilized the MPM.

**Remark :** Every fields are defined at nodes in DGMPM and constitutive equations are integrated in fluxes contained in volume and surface integrals whereas in MPM, velocity is updated at nodes and constitutive equations are integrated at material points.

# 4 INTERFACE FLUXES COMPUTATION

Fields continuity is relaxed across elements boundary and interface fluxes appear in the weak form. They can be used to introduce the characteristic structure of hyperbolic problems in the numerical scheme. One possibility is to solve a Riemann problem normal to each edge of a given cell and to compute the Godunov flux corresponding to the stationary solution ($x/t = 0$) of the Riemann problem [7]. For the two-dimensional case (fig 2), the problem at a given edge consists in finding the solution of the auxiliary initial value problem (IVP):
$$ \frac{\partial \tilde{q}}{\partial t} + \frac{\partial \mathcal{F}_n}{\partial \xi} = 0, \quad \xi = \mathbf{X} \cdot \mathbf{n} \in ]-\infty, \infty[, \quad t > 0 \qquad (18) $$

$$ \tilde{q}(\xi, 0) = \begin{cases} \tilde{q}_L & \text{if } \xi < 0 \\ \tilde{q}_R & \text{if } \xi > 0 \end{cases} \qquad (19) $$

where L and R denote states belonging to the Left cell and the Right cell to the edge respectively. These states are constructed by averaging the values defined at nodes for each side (remind that the DG approximation gives one nodal field value for each cell connected to a given node). The

<figure><img src="image_3.png" /><figcaption>Figure 2: Edge separating two cells indexed L (left) and R (right). A local coordinate $\xi$ and a normal vector $\mathbf{n}$ are defined.</figcaption></figure>

normal flux vector can be decomposed with the consideration that $\mathbf{n} = n_1\mathbf{e}_1 + n_2\mathbf{e}_2$:

$$ \mathcal{F}_n = \mathcal{F} \cdot \mathbf{n} = n_1 \mathcal{F}_1 + n_2 \mathcal{F}_2 \qquad (20) $$

The conservation laws system (18) combined with (20) yields:

$$ \frac{\partial \tilde{q}}{\partial t} + \sum_p n_p \frac{\partial \mathcal{F}_p}{\partial \xi} = 0 $$

$$ \frac{\partial \tilde{q}}{\partial t} + \sum_p n_p \frac{\partial \mathcal{F}_p}{\partial \tilde{q}} \cdot \frac{\partial \tilde{q}}{\partial \xi} = 0 $$

$$ \frac{\partial \tilde{q}}{\partial t} + \sum_p n_p \mathbf{A}^p \cdot \frac{\partial \tilde{q}}{\partial \xi} = 0 $$

The characteristic analysis of this system written with vectors (6) can be conducted with the study of the jacobian matrix $J = \sum_p n_p \mathbf{A}^p$. In two dimensions:

$$ \mathbf{A}^p = \begin{bmatrix} \mathbf{0} & \mathbf{0} & -H_{ip1j} \mathbf{e}_i \otimes \mathbf{e}_j \\ \mathbf{0} & \mathbf{0} & -H_{ip2j} \mathbf{e}_i \otimes \mathbf{e}_j \\ -\frac{\delta_{1p}}{\rho_0} \mathbf{1} & -\frac{\delta_{2p}}{\rho_0} \mathbf{1} & \mathbf{0} \end{bmatrix} \qquad (21) $$

and $H_{klmn}$ are the fourth-order material elasticity tensor components. The eigenvalues and eigenvectors of the Jacobian matrix are used to obtain the solution of the Riemann problem for linear (elastic) problems, or to build an approximate Riemann Solver [7] for hyperelastic nonlinear problems. The solution is used to compute Godunov fluxes at cells interfaces $\mathcal{F}_n(\tilde{q}(x/t = 0))$.
# 5 NUMERICAL RESULTS

## 5.1 One dimensional compression wave in a bar (small strain)

Let us consider an elastic bar of length $l = 6m$ made of a linear elastic material with Young's modulus $E = 2.10^{11}$ Pa, density $\rho = 7800$ kg.m$^{-3}$ and unit cross section area $S$. The left end of the bar is suddenly loaded with a surface force $T^d = 1.10^9$ Pa. The bar is discretized with 150 regular grid cells (each containing one single centered material point) and the time integration is carried out with an explicit Euler scheme. This problem has been solved with DGMPM and a comparison is performed with MPM and the analytical solution [8] for different time steps (fig 3).

<figure><img src="image_4.png" /><figcaption>Figure 3: Axial stress in a bar at different times. Comparison of DGMPM, MPM and analytical solution.</figcaption></figure>

The solution consists of a discontinuous elastic disturbance propagating rightward at speed $c = \sqrt{\frac{E}{\rho}}$ and is here accurately captured by the DGMPM solution. It is known that a numerical scheme is able to capture a discontinuity providing the Courant number be one without loss of stability ($\alpha = c \frac{\Delta t}{\Delta x}$). The MPM is unstable with $\alpha = 1$ and hence, it is unable to track discontinuities. The stability analysis of DGMPM will not be developed here, but it can be shown that when a single material point is located in each cell, the Courant number can be equal to one and a wave front can be represented. This property is lost with more points in each cell but it can be improved by using a Runge-Kutta 2 time integrator.
## 5.2 Two dimensional partial impact on a plate (small strain)

We are now interested in two-dimensional plane strains cases and the problem treated concerns a linear elastic rectangular plate undergoing an impact on a part of its left boundary (fig 4). The plate dimensions are $l \times h = 4 \times 3 \, m^2$ discretized with 33 x 25 material points, and the width of impact area is $a = 1 \, m$. A velocity Dirichlet boundary condition is enforced on material points directly in the MPM while ghost nodes are used in the DGMPM grid [4]. The problem has been solved with the finite element code Cast3M [9] (explicit time integration), the material point method and its discontinuous Galerkin extension and the comparison is shown in fig 5 for one point per cell.

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Elastic plate undergoing a compressive impact on a part of its left boundary. Parameters: $E = 2.10^{11} Pa$, $\nu = 0.3$, $\rho = 7800 kg \cdot m^{-3}$, $v^d = 5m.s^{-1}$.</figcaption></figure>

As mentioned above about FEM and MPM, the time integration used involves oscillations in the solution and the Courant condition ($\alpha < 1$) introduces numerical diffusion. This two
<figure><img src="image_7.png" /><figcaption>Figure 5: Longitudinal stress profile resulting from MPM (top left), FEM (top right) and DGMPM (bottom left) and superposed plot along the bottom boundary of the plate at different times for DGMPM and MPM and FEM (bottom right). $\alpha_{MPM} = 0.5$, $\alpha_{FEM} = 0.9$, $\alpha_{DGMPM} = 1$</figcaption></figure>
dimensional problem involves two families of waves : tensile/compression waves traveling at celerity $c_p = \sqrt{(\lambda + 2\mu)/\rho}$ and shear waves traveling at celerity $c_s = \sqrt{\mu/\rho}$. These waves can make the DGMPM solution scheme unstable when the Courant number equals one. To avoid it, a transverse Riemann solver based on the Corner Transport Upwind method [4] has been used and sharp solutions can be captured more accurately. These results confirm that the extension of MPM to DG framework provides a better wave front tracking.

## 5.3 One dimensional compression wave in a hyperelastic bar (large strain)

The last numerical test is an extension of the first one to the large strains framework. The bar is now made of a hyperelastic Saint-Venant-Kirchhoff material for which the constitutive equation is :

$$ \Pi = 2\mu F \cdot E + \lambda \operatorname{trace}(E)F \quad (22) $$

where $E = \frac{1}{2}(F^T \cdot F - 1)$ is the Green-Lagrange strain tensor and $\Pi$ is the first Piola-Kirchhoff stress tensor. For a one-dimensional bar, this relation reduces to : $\Pi_{xx} = \frac{E}{2}F_{xx}(F_{xx}^2 - 1)$. The bar is suddenly loaded on its left side ($T^d = 8.10^8 Pa$). In figure 6, the upper subplot shows the material points positions at different times updated with the velocity solution. The bottom one shows the Piola-Kirchhoff stress from MPM, DGMPM (both in total Lagrangian formulation) and an analytical solution.

<figure><img src="image_8.png" /><figcaption>Figure 6: Saint-Venant-Kirchhoff hyperelastic bar undergoing large deformations (positions - Piola stress) plot at different times. Material parameters are the same as in the linear elastic case.</figcaption></figure>

The DGMPM is no longer exactly on the analytical solution of this nonlinear problem and the reason might be the use of an approximate Riemann solver instead of the nonlinear one (which is more expensive). Once again classical MPM fails to represent the discontinuity and introduces oscillations in the results.

# 6 CONCLUSION

In this work, the Material Point Method has been extended to the Discontinuous Galerkin approximation framework in order to improve wave front tracking. This extension is based on the use of shape functions defined element-by-element and discontinuous across cells interfaces.
The weak form of the problem within the large strains framework is written for the equation of motion (4) and for the lagrangian kinematic equation (3). The DGMPM hyperbolic system is built with a total lagrangian formulation (13) and contains numerical fluxes at cells interfaces which are computed with the solution of a Riemann problem. The main difference with the original scheme is that the constitutive equations are contained in nodal fluxes and the characteristic structure of hyperbolic problems can be introduced in the numerical solution, providing a better description of discontinuities.

## REFERENCES

[1] D. Sulsky, S.-J. Zhou, and H. L. Schreyer, “Application of a particle-in-cell method to solid mechanics,” *Computer Physics Communications*, vol. 87, no. 12, pp. 236–252, 1995.

[2] W. H. Reed and T. Hill, “Triangular mesh methods for the neutron transport equation,” *Los Alamos Report LA-UR-73-479*, 1973.

[3] J. S. Hesthaven and T. Warburton, *Nodal discontinuous Galerkin methods: algorithms, analysis, and applications*. Springer Science & Business Media, 2007.

[4] R. J. Leveque, *Finite volume methods for hyperbolic problems*. Cambridge university press, 2002.

[5] B. J. Plohr and D. H. Sharp, “A conservative eulerian formulation of the equations for elastic flow,” *Advances in Applied Mathematics*, vol. 9, no. 4, pp. 481–499, 1988.

[6] E. Love and D. Sulsky, “An unconditionally stable, energymomentum consistent implementation of the material-point method,” *Computer Methods in Applied Mechanics and Engineering*, vol. 195, no. 3336, pp. 3903–3925, 2006.

[7] E. F. Toro, *Riemann solvers and numerical methods for fluid dynamics: a practical introduction*. Springer Science & Business Media, 2013.

[8] L. Wang, *Foundations of stress waves*. Elsevier, 2011.

[9] Cast3M. User's manual, 2016.