# A Discontinuous Galerkin Material Point Method for the solution of impact problems in solid dynamics

Adrien Renaud, Thomas Heuzé, Laurent Stainier

► To cite this version:

Adrien Renaud, Thomas Heuzé, Laurent Stainier. A Discontinuous Galerkin Material Point Method for the solution of impact problems in solid dynamics. Journal of Computational Physics, 2018, 369, pp.80-102. 10.1016/j.jcp.2018.05.001. hal-01913188

HAL Id: hal-01913188

https://hal.science/hal-01913188

Submitted on 8 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>

Distributed under a Creative Commons Attribution 4.0 International License
# A Discontinuous Galerkin Material Point Method for the solution of impact problems in solid dynamics

Adrien Renaud¹, Thomas Heuzé, Laurent Stainier

¹ Research Institute in Civil and Mechanical Engineering (GeM, UMR 6183 CNRS)

Ecole Centrale de Nantes

1 rue de la Noë, Nantes, France

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

## Abstract

An extension of the Material Point Method [1] based on the Discontinuous Galerkin approximation (DG) [2] is presented here. A solid domain is represented by a collection of particles that can move and carry the fields of the problem inside an arbitrary computational grid in order to provide a Lagrangian description of the deformation without mesh tangling issues. The background mesh is then used as a support for the Discontinuous Galerkin approximation that leads to a weak form of conservation laws involving numerical fluxes defined at element faces. Those terms allow the introduction of the characteristic structure of hyperbolic problems within the numerical method by using an approximate Riemann solver [3]. The Discontinuous Galerkin Material Point Method, which can be viewed as a Discontinuous Galerkin Finite Element Method (DGFEM) with modified quadrature rule, aims at meeting advantages of both mesh-free and DG methods. The method is derived within the finite deformation framework for multidimensional problems by using a total Lagrangian formulation. A particular attention is paid to one specific discretization leading to a stability condition that allows to set the CFL number at one. The approach is illustrated and compared to existing or developed analytical solutions on one-dimensional problems and compared to the finite element method on two-dimensional simulations.

**Keywords:** Discontinuous Galerkin approximation; Material Point Method; impact; solid mechanics; hyperelasticity; finite deformation

## 1. Introduction

A wide variety of engineering problems in solid dynamics, such as impact problems or high speed forming techniques are modeled with hyperbolic initial boundary value problems (IBVP). The solutions of this class of mathematical problems involve waves which are important to accurately follow in order to understand physical phenomena occurring within the medium during the deformation. Indeed, knowledge of waves paths provides information about the loading history undergone by each particle, which in turn allows to properly compute the residual strains and stresses, and to optimize the final shape of the body. Furthermore high speed forming processes such as electromagnetic forming [4] often involve finite deformation than can complicate the development of an efficient numerical scheme. A Lagrangian mesh approach may require re-meshing and projection steps in order to deal with mesh entanglement effects. Those issues can be avoided by using an Eulerian method. However interface tracking techniques and convection steps are required in order to follow the boundaries and transport internal variables which is less convenient for solid mechanics.

The numerical simulation of IBVP has been and is still widely performed in solid mechanics with the *Finite Element Method (FEM)* [5] due to, among other assets, its ability to easily deal with nonlinearities such as history-dependent constitutive models [6]. However, despite its well-known advantages a Lagrangian finite element formulation becomes less efficient and accurate for problems involving very large strains due to difficulties mentioned above. *The Material Point Method (MPM)* developed in the 90's [1] can be viewed as an extension of the FEM in which quadrature points can move in an arbitrary computational grid. Those *material points* represent a continuum with a Lagrangian description and are used to store all the fields and internal variables of the problem. The underlying mesh
is only used to compute gradients, apply boundary conditions and solve discrete equations. The method hence uses
many transfer steps of fields between material points and nodes. Since the data are stored at material points, the grid
can be discarded at each time step and re-constructed for computational convenience without additional projections,
thus avoiding mesh entanglement. In addition to problems caused by finite deformations, classical explicit time inte-
grators used in solid dynamics with FEM and MPM introduce high frequency noise in the vicinity of discontinuities.
The removal of these spurious numerical oscillations with additional viscosities is difficult to achieve without loss
of accuracy, and can be troublesome for the wave tracking. Hence, although the material point method has been
successfully applied to a wide range of complex large strain engineering problems [7] this approach seems unable to
accurately deal with discontinuities owing to numerical oscillations, diffusion introduced by mapping steps, and its
low maximum admissible Courant number.

On the other hand, the *Finite Volume Method (FVM)* [8] initially developed for fluid dynamics has proven to be a very efficient tool in order to accurately follow the path of waves. Indeed, the formulation involves a conservative form leading to the same order of accuracy for velocity and gradients, and numerical fluxes allowing to account for the characteristic structure of hyperbolic equations. Moreover, several works on FVM provided high order schemes satisfying the *Total Variation Diminishing (TVD)* stability condition [9] and recent studies extended the approaches to solid mechanics for problems involving history-dependent models [10, 11], or finite deformations with a Lagrangian formulation [12, 13].

The *discontinuous Galerkin (DG)* approximation [2] enables to build numerical schemes that benefit from both finite element and finite volume methods. DG-methods are based on piece-wise polynomial shape functions and their combination with the finite element formulation (DGFEM) led to several development steps to reach *Total Variation Diminishing in the Means (TVDM)* and *Total Variation Bounded (TVB)* high order schemes that ensure convergence to entropy-satisfying solutions [14]. In addition, the same order of accuracy is reached for velocity and gradients within a finite element framework when the weak form is based on a conservation laws system. Although this approach can easily handle mesh-adaption strategies due to the relaxation of fields continuity, the mesh tangling problems do not vanish. Furthermore, the Courant condition that applies is restrictive [14] and prevents to accurately capture discontinuities.

The purpose of this work is to develop a numerical method gathering the strengths of both material point and DG methods and avoiding their respective limits by means of the discontinuous Galerkin approximation. The computational grid of the MPM is used to build a piece-wise linear approximation of the solution of conservation laws, which weak form is written on each element. Boundary integrals arising from integration by parts involve numerical fluxes that connect the elements together and are computed with an approximate Riemann solver [3]. Projection steps between the underlying mesh, used to solve the discrete equations, and material points that store the fields of the problem, provide a Lagrangian description of the deformation in an arbitrary computational grid. The resulting *Discontinuous Galerkin Material Point Method (DGMPM)* hence aims to enable an accurate wave paths tracking in solid media undergoing large strains. A similar approach based on the *Particle-In-Cell (PIC)* method has already been developed in [15] for electromagnetism equations. This discontinuous Galerkin version of PIC, unlike the DGMPM, uses classical Gauss points to integrate the weak form. Moreover, electric and magnetic fields are stored at nodes whereas charge and current densities are stored at particles which prevent the arbitrariness of the grid. In what follows, after a brief review of continuum equations (section 2), the DGMPM formulation is derived within the large strains framework with a total Lagrangian approach for hyperelasticity in section 3. Section 4 is dedicated to the computation of interface fluxes arising in the weak form. First, Riemann problems are written at interfaces with a *quasi-linear* form which spectral analysis is performed. Second, an approximate Riemann solver that can take into account transverse corrections for multi-dimensional problems is derived. Finally, the enforcement of boundary conditions is discussed. Next, after a summary of the numerical procedure in section 5, the stability analysis of the one-dimensional scheme is carried out in section 6. This analysis shows that a reduced integration of the weak form provides an improvement of the stability condition that limits the DGFEM. At last, the method is illustrated in section 7 on one-dimensional cases in order to compare it to analytical solutions, and two-dimensional problems for which comparisons are performed with the Finite Element Method.
## 2. The continuum problem: isothermal hyperelastodynamic

### 2.1. Conservation equations

We consider a three-dimensional isothermal deformation of a homogeneous continuum elastic body from the reference configuration $\Omega_0 \in \mathbb{R}^3$, to the configuration of the same body at a different time $\Omega_t \in \mathbb{R}^3$. The motion of the body is described by the smooth mapping $\varphi(X, t) \equiv x$ from the reference coordinates system $X$ to the spatial coordinates system $x$. The second-order two-point deformation gradient tensor and the material velocity vector are respectively defined as:

$$F_{i\alpha}(X, t) := \frac{\partial \varphi_i(X, t)}{\partial X_{\alpha}} \qquad (1)$$

$$v_i(X, t) := \dot{\varphi}_i(X, t) \qquad (2)$$

where the superposed dot denotes the material time derivative, Greek indices stand for material coordinates and Latin ones for spatial coordinates. Combination of equations (1) and (2) leads to the Lagrangian kinematic equations:

$$\dot{F}(X, t) = \frac{\partial v(X, t)}{\partial X} \quad \forall X \in \Omega_0 \qquad (3)$$

The solid must locally satisfy the conservation of linear momentum written in the reference configuration within the time interval of interest $t \in \tau = [0, T]$:

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

where $\rho_0(X) = \rho(X, t = 0)$ is the reference mass density, $\Pi$ is the first Piola-Kirchhoff stress tensor, $\rho_0 b$ is the body forces vector and $\nabla_X \cdot (\bullet)$ is the divergence vector with respect to reference coordinates. For simplicity, body forces are neglected. The use of the Cartesian coordinates system leads to a homogeneous system of partial differential equations written with equations (3) and (4) [16]:

$$\frac{\partial \mathcal{U}}{\partial t} + \frac{\partial \mathcal{F}_{\alpha}}{\partial X_{\alpha}} = 0 \quad \forall X \in \Omega_0 \qquad (5)$$

where $\mathcal{U}$ is the vector of conserved quantities and $\mathcal{F}_{\alpha}$ is the flux vector in the $\alpha^{th}$ direction:

$$\mathcal{U} = \begin{bmatrix} F \\ \rho_0 v \end{bmatrix} ; \quad \mathcal{F} = - \begin{bmatrix} v \otimes I \\ \Pi \end{bmatrix} ; \quad \mathcal{F}_{\alpha} = \mathcal{F} \cdot e_{\alpha} = - \begin{bmatrix} v \otimes e_{\alpha} \\ \Pi \cdot e_{\alpha} \end{bmatrix} \qquad (6)$$

In what follows, vectors (resp. tensors) denoted by typographic symbols belong to $\mathbb{R}^p$ (resp. $\mathbb{R}^p \times \mathbb{R}^p$) where $p = 2, 3$, whereas calligraphic symbols denote column arrays.

### 2.2. Constitutive model

We are interested here in modeling hyperelastic materials characterized by a local objective stored-energy function $\Psi$. In order to ensure mathematical well-posedness of the problem this function is assumed to be polyconvex [17] namely, $\Psi$ is convex with respect to the deformation gradient, its determinant $J$ and its co-factor tensor $H = JF^{-T}$. Such a stored-energy can then be expressed as a function of the deformation gradient only, leading to the classical definition of hyperelastic materials:

$$\Pi(X, F) = \frac{\partial \Psi(X, F)}{\partial F} \qquad (7)$$

The fourth-order tangent modulus tensor $C$ and the second-order acoustic tensor $A$ are respectively defined as:

$$(C)_{i\alpha j\beta} = \frac{\partial \Pi_{i\alpha}}{\partial F_{j\beta}} = \frac{\partial^2 \Psi(X, F)}{\partial F_{j\beta} \partial F_{i\alpha}} \qquad (8)$$

$$(A)_{ij} = N_{\alpha} C_{i\alpha j\beta} N_{\beta} \qquad (9)$$

where $N_{\alpha}$ are the components of an arbitrary material normal vector $N$.
### 3. Space-time discretization

As for MPM, a continuum body $\Omega_t$ is discretized into a set of $P$ material points in an arbitrary Cartesian grid. The grid is here supposed to be made of $E$ non-overlapping elements, or cells, for convenience. In this grid, the set of faces separating cells that contain particles and having empty neighboring elements defines the boundary of the mesh (see figure 1 for a two-dimensional case). In this section, the formulation of the new approach is developed within the large strains framework.

<figure><img src="image_2.png" /><figcaption>Figure 1: Representation of a continuum body by a set of material points in $\mathbb{R}^2$.</figcaption></figure>

#### 3.1. Weak form of the system of conservation laws

We seek an approximate solution to system (5) in a finite-dimensional function space using discontinuous Galerkin approach. The key idea of DG methods is to allow jump of fields across mesh elements faces by using broken polynomials spaces [2, 18, 19]. In what follows, we will only consider linear polynomials and the corresponding broken spaces:

$$V^1 = \{\mathcal{V} \in H^1(\Omega_e)\} ; \quad V_h^1 = \{\mathcal{V} \in \mathcal{P}^1(\Omega_e)\} \subset V^1 \quad (10)$$

where $H^1(\Omega_e)$ and $\mathcal{P}^1(\Omega_e)$ denote respectively Sobolev space and the linear (Lagrange) polynomials space, in $\Omega_e$. A weak form is written by multiplying equation (5) by a test function $\mathcal{V} \in V_h^1$ and integrating over each element $e$ of the mesh:

$$\int_{\Omega_e} \frac{\partial \mathcal{U}}{\partial t} \cdot \mathcal{V} \, d\Omega + \int_{\Omega_e} \frac{\partial \mathcal{F}_\alpha}{\partial X_\alpha} \cdot \mathcal{V} \, d\Omega = 0 \quad \forall \mathcal{V} \in V_h^1, \forall e \in [1, E] \quad (11)$$

Integration by parts leads to:

$$\int_{\Omega_e} \frac{\partial \mathcal{U}}{\partial t} \cdot \mathcal{V} \, d\Omega - \int_{\Omega_e} \mathcal{F}_\alpha \cdot \frac{\partial \mathcal{V}}{\partial X_\alpha} \, d\Omega + \int_{\Gamma_e} (\mathcal{F} \cdot N) \cdot \mathcal{V} \, d\Gamma = 0 \quad \forall \mathcal{V} \in V_h^1, \forall e \in [1, E] \quad (12)$$

where $N$ is the outward normal vector to the boundary, and $\mathcal{F} \cdot N = \mathcal{F}_N$ denotes a normal flux at cell faces. Those numerical fluxes can be computed from the solution of an approximate Riemann solver [20] (see section 4). In particular, the stationary solution of Riemann problems defined at element boundaries $\Gamma_i$ such that $\cup_i \Gamma_i = \Gamma_e$ gives the well-known Godunov's method.

Specific quantities are introduced and delta Dirac distributions are used as particles characteristic for the mass density [1]:

$$\mathcal{U} = \rho_0 \tilde{\mathcal{U}} ; \quad \mathcal{F}_\alpha = \rho_0 \tilde{\mathcal{F}}_\alpha \quad (13)$$

$$\rho_0(\mathbf{X}) = \sum_{p=1}^{P} m_p \delta(\mathbf{X}_p - \mathbf{X}) \quad (14)$$
where $m_p$ is the mass of the $p^{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, and equation (12) thus reads:

$$ \int_{\Omega_e} \rho_0 \frac{\partial \bar{\mathcal{U}}}{\partial t} \cdot \mathcal{V} \, d\Omega - \int_{\Omega_e} \rho_0 \bar{\mathcal{F}}_\alpha \cdot \frac{\partial \mathcal{V}}{\partial X_\alpha} \, d\Omega + \int_{\Gamma_e} \mathcal{F}_N \cdot \mathcal{V} \, d\Gamma = 0 \quad \forall \mathcal{V} \in \mathcal{V}_h^1, \forall e \in [1, E] \quad (15) $$

which, with Dirac integration properties transforms into:

$$ \sum_{p=1}^{p} m_p \left[ \left. \frac{\partial \bar{\mathcal{U}}}{\partial t} \cdot \mathcal{V} - \bar{\mathcal{F}}_{\alpha} \cdot \frac{\partial \mathcal{V}}{\partial X_{\alpha}} \right|_{X=X_p} + \int_{\Gamma_e} \mathcal{F}_N \cdot \mathcal{V} \, d\Gamma = 0 \quad \forall \mathcal{V} \in \mathcal{V}_h^1, \forall e \in [1, E] \quad (16) $$

## 3.2. Semi-discrete and discrete equations

The particles are viewed as interpolation points so that fields of the weak form (16) are evaluated from nodal values:

$$ \bar{\mathcal{U}}(X_p) = \bar{\mathcal{U}}_p = \sum_{i} S_i(X_p) \bar{\mathcal{U}}_i = \sum_{i} S_{ip} \bar{\mathcal{U}}_i \quad i = 1, \dots, N_{\text{nodes}} \quad (17) $$

where $S_{ip}$ denotes the shape function attached to the $i^{th}$ node evaluated at $X = X_p$, and $\bar{\mathcal{U}}_i$ the specific vector of conserved quantities at node $i$. In equation (17) and in the following of the paper, index $p$ stands for a quantity defined at material points whereas $(i, j)$ stand for nodal quantities. Combination of equations (16)-(17) and arbitrariness of the test field yield:

$$ \sum_{p=1}^{p} \left[ S_{ip} m_p S_{jp} \frac{\partial \bar{\mathcal{U}}_j}{\partial t} - \frac{\partial S_{ip}}{\partial X_\alpha} m_p S_{jp} \bar{\mathcal{F}}_{\alpha,j} \right] + \int_{\Gamma_e} S_i(X) \mathcal{F}_N \, d\Gamma = 0 \quad \forall e \in [1, E] \quad (18) $$

$$ \Leftrightarrow M_{ij} \frac{\partial \bar{\mathcal{U}}_j}{\partial t} - K_{ij}^{\alpha} \bar{\mathcal{F}}_{\alpha,j} + \hat{\mathcal{F}}_i = 0 \quad (19) $$

which is the semi-discrete system that must be solved on the grid. Note that the consistent mass matrix $M_{ij}$ can be singular when only one material point lies in an element (reduced integration) [21]. However the use of a diagonally lumped mass matrix permits to avoid the singularity problem, and in what follows the lumped matrix will be written: $M_i^L = \sum_j M_{ij}$. Furthermore this configuration is equivalent to the DGFEM [22] with a mass matrix integrated with only one Gauss point.

Finally, the time interval $\tau$ is discretized in $N_T$ sub-intervals of size $\Delta t^n$ such that $\sum_{n=1}^{N_T} \Delta t^n = T$ and time integration is performed with an explicit forward Euler algorithm leading to discrete equations:

$$ M_i^L \frac{\bar{\mathcal{U}}_i^{n+1} - \bar{\mathcal{U}}_i^n}{\Delta t^n} = K_{ij}^{\alpha} \bar{\mathcal{F}}_{\alpha,j}^{n} - \hat{\mathcal{F}}_i^n \quad (20) $$

where the superscripts denotes the time step at which a field is evaluated.

## 4. Computation of interface fluxes

Interface fluxes involved in boundary integrals of equation (18) can be used to introduce the characteristic structure of hyperbolic problems within the numerical scheme, providing that they satisfy some mathematical requirements, namely they must be E-fluxes [23]. One possibility is to solve a Riemann problem normal to each boundary of a given cell and to compute the Godunov flux of the solution. An approximate Riemann solver that computes a solution which can be used to calculate those fluxes is presented in this section.
### 4.1. The Riemann problem

For the multi-dimensional case the problem at a given face consists in finding the solution of the following auxiliary initial value problem (IVP):

$$ 
\begin{gathered}
\frac{\partial \tilde{\mathcal{U}}}{\partial t} + \frac{\partial \mathcal{F}_N}{\partial \xi} = \mathbf{0}, \quad \xi \in [-\infty, \infty[, \quad t > 0 \\
\tilde{\mathcal{U}}(\xi, 0) = \begin{cases} \tilde{\mathcal{U}}_L & \text{if } \xi < 0 \\ \tilde{\mathcal{U}}_R & \text{if } \xi > 0 \end{cases}
\end{gathered}
\qquad (21) $$

where the top equation is the projection of system (5) along the normal vector to an interface with the auxiliary variable $\xi = X \cdot N$. Initial conditions $(\tilde{\mathcal{U}}_L, \tilde{\mathcal{U}}_R)$ arise from the DG approximation which introduces multi-valued fields on interfaces and can be viewed as a duplication of nodes (figure 2). Due to the expressions of flux and conserved quantities vectors (6), the calculation of Godunov fluxes based on the solution $\tilde{\mathcal{U}}$ to problem (21) requires the computation of constitutive equations. However, this procedure can be costly for nonlinear problems, especially for history-dependent models. Hence, the quasi-linear form is written introducing an auxiliary vector $\mathcal{Q}$ [24]:

$$ 
\begin{gathered}
\frac{\partial \tilde{\mathcal{U}}}{\partial \mathcal{Q}} \frac{\partial \mathcal{Q}}{\partial t} + \frac{\partial \mathcal{F}_N}{\partial \mathcal{Q}} \frac{\partial \mathcal{Q}}{\partial \xi} = \mathbf{0}, \quad \xi \in [-\infty, \infty[, \quad t > 0 \\
\mathcal{Q}(\xi, 0) = \begin{cases} \mathcal{Q}_L & \text{if } \xi < 0 \\ \mathcal{Q}_R & \text{if } \xi > 0 \end{cases}
\end{gathered}
\qquad (22) $$

A wise choice of the auxiliary vector, namely $\mathcal{Q}$ is a rearrangement of the flux components, permits to compute the interface flux from the solution $\mathcal{Q}^*$ of problem (22) with a simple assembly:

$$ \mathcal{F}_N^* = \frac{\partial \mathcal{F}_N}{\partial \mathcal{Q}} \mathcal{Q}^* \qquad (23) $$

In what follows, the auxiliary vector will be taken as $\mathcal{Q}^T = [\Pi^T, \nu^T]$. One single Riemann problem is solved at each interface center with initial conditions obtained by averaging the nodal values on each side as depicted in figure 2, so that the computational cost is reduced.

<figure><img src="image_3.png" /><figcaption>Figure 2: Duplication of nodes at an interface and building of initial conditions of the Riemann problem (2D).</figcaption></figure>

IVP (22) admits self-similar solutions, constant along each ray traced from the origin of the $(\xi, t)$ plane. A linearized Riemann solver allows to determine explicitly at a moderate cost an approximate value of the stationary solution [3] providing that the characteristic structure of the solution is known.

### 4.2. Spectral analysis

The characteristic structure of the solution of problem (22) is given by the spectral analysis of the Jacobian matrix defined for the general three-dimensional case as:

$$ \mathcal{J}_N^Q(\mathcal{Q}) = \left(\frac{\partial \tilde{\mathcal{U}}}{\partial \mathcal{Q}}\right)^{-1} \frac{\partial \mathcal{F}_N}{\partial \mathcal{Q}} = - \left[ \begin{array}{cc} \mathbf{0}_4 & \frac{\partial \Pi}{\partial F} \cdot N \\ \frac{1}{\rho_0} \mathbf{I} \otimes N & \mathbf{0}_2 \end{array} \right] \qquad (24) $$
where $\mathbf{0}_2$ and $\mathbf{0}_4$ are null second-order and fourth-order tensors respectively. In what follows, the dependency of the Jacobian matrix on $\mathcal{Q}$ is omitted to simplify the notation. The problem (21), and hence (22), is known to be hyperbolic if and only if [25]:

$$ \frac{\partial^2 \Psi}{\partial F_{i\alpha} \partial F_{j\beta}} N_{\alpha} N_{\beta} m_i m_j > 0 \quad \forall \mathbf{N}, \mathbf{m} \in \mathbb{R}^3; \mathbf{N}, \mathbf{m} \neq \mathbf{0} \qquad (25) $$

which is known as the *Legendre-Hadamard or ellipticity condition*. The prerequisite that the stored-energy function is polyconvex ensures the satisfaction of the Legendre-Hadamard condition (25) [17] hence, hyperbolicity is guaranteed. Moreover, the Jacobian matrix (24) has in this case 6 non zero and 6 null eigen-values with associated independent eigen-vectors [25]. The *M* eigen-values of the Jacobian matrix (24) give rise to *M* characteristic fields given by the resolution of the following eigen-system:

$$ \mathcal{J}_N^\mathcal{Q} \mathcal{R}_\mathcal{Q}^K = \lambda_K \mathcal{R}_\mathcal{Q}^K ; \quad \mathcal{R}_\mathcal{Q}^K = \begin{bmatrix} \mathbf{\Pi}^K \\ \mathbf{v}^K \end{bmatrix} \qquad (26) $$

$\mathcal{R}_\mathcal{Q}^K$ being the right eigen-vector associated to the $K^{th}$ eigen-value $\lambda_K$ of $\mathcal{J}_N^\mathcal{Q}$. Each of these characteristic fields corresponds to one wave emerging from the origin of the $(\xi, t)$ plane. Equation (26) leads to the two independent systems:

$$ -C_{i\alpha j\beta} N_{\beta} v_j^K = \lambda_K \Pi_{i\alpha}^K \qquad (27) $$

$$ -\frac{\Pi_{i\alpha}^K}{\rho_0} N_{\alpha} = \lambda_K v_i^K \qquad (28) $$

and by introducing equation (28) in (27), one gets the acoustic tensor eigen-system:

$$ N_{\alpha} C_{i\alpha j\beta} N_{\beta} Y_{j}^{p} = \omega_{p} Y_{i}^{p} \quad p = 1, \dots, 3 \qquad (29) $$

where $(\omega_p, Y^p)$ are the acoustic tensors eigen-values $(\omega = \rho_0 \lambda^2)$ and associated eigen-vectors. Since Legendre-Hadamard condition (25) is satisfied, the acoustic tensor admits three positive eigen-values leading to six characteristic speeds $\lambda_{2p, 2p-1} = \pm \sqrt{\omega_p / \rho_0}$, $(p = 1, 2, 3)$. We then see that one eigen-vector $Y^p$ of the acoustic tensor is associated to two eigen-values $\lambda_K$ and two velocity vectors $v^K$, thus it comes: $Y^p = v^{2p-1} = v^{2p}$. Finally, the stress components of the Jacobian eigen-vectors $\Pi^K$ are determined with equation (27). Then, we have for the non zero eigen-values the following characteristic fields:

$$ \left\{ \pm \sqrt{\frac{\omega_p}{\rho_0}} \quad ; \quad \begin{bmatrix} -C : (\mathbf{Y}_p \otimes \mathbf{N}) \\ \pm \sqrt{\frac{\omega_p}{\rho_0}} \mathbf{Y}_p \end{bmatrix} \right\} \qquad (30) $$

The last step consists in building a complete basis of the kernel of $\mathcal{J}_N^\mathcal{Q}$ by finding six independent vectors solving equation (28) for the null eigen-value:

$$ \Pi_{i\alpha}^K N_{\alpha} = 0 \qquad (31) $$

### 4.3. Approximate-state Riemann solver

The solution of problem (22) consists of *M* discontinuous waves emanating from the origin of the $(\xi, t)$ plane, from which the stationary solution reads:

$$ \mathcal{Q}^* = \mathcal{Q}_L + \sum_{\substack{m=1 \\ \lambda_m < 0}}^{M} \mathcal{W}_m = \mathcal{Q}_R - \sum_{\substack{m=1 \\ \lambda_m > 0}}^{M} \mathcal{W}_m \qquad (32) $$

where $\mathcal{W}_m = \delta_m \mathcal{R}_\mathcal{Q}^m$ are the waves formed of the wave strength coefficients $\delta_m$ weighting the right eigen-vectors $\mathcal{R}_\mathcal{Q}^m$ associated to the *M* eigen-values $\lambda_m$ of the Jacobian matrix (24). By rearranging equation (32), one can compute the wave strength coefficients by solving:

$$ \mathcal{Q}_R - \mathcal{Q}_L = \sum_{m=1}^{M} \delta_m \mathcal{R}_\mathcal{Q}^m \qquad (33) $$

and then calculate the stationary solution $\mathcal{Q}^*$ and associated normal fluxes $\mathcal{F}_N$ with assembly (23). Those fluxes are assumed to be constant on the face and the boundary integrals of the semi-discrete equation (18) can be computed.
### 4.4. Transverse Riemann solver

The method derived above for the computation of normal fluxes can be viewed as the donor-cell Upwind (DCU) method [8] meaning that only contributions from the upwind cells sharing an edge (in two dimensions) with the current one are considered. For multidimensional problems waves can travel in several directions such that contributions coming from corner cells must be taken into account in order to improve accuracy and stability of the numerical scheme. The corner transport upwind (CTU) method [26] consists in considering contributions coming from upwind cells sharing only a node (in two dimensions) with the considered grid cell and propagating in bias. This approach allows to improve the Courant condition especially for solid mechanics problems for which strain components are coupled through Poisson's effect. One can define at each cell interface left-going and right-going fluctuations defined as:

$$ \mathcal{A}^{-}(\Delta \mathcal{U}) = \mathcal{F}(\mathcal{U}^{*}) - \mathcal{F}(\mathcal{U}_{\xi^{-}}) \quad ; \quad \mathcal{A}^{+}(\Delta \mathcal{U}) = \mathcal{F}(\mathcal{U}_{\xi^{+}}) - \mathcal{F}(\mathcal{U}^{*}) \quad (34) $$

Let's consider the patch of grid cells shown in figure 3, and focus on the edge denoted (i) which local normal vector

<figure><img src="image_4.png" /><figcaption>Figure 3: Normal and transverse fluctuations defined from edge i.</figcaption></figure>

$N^i$ is shown. The Riemann problem defined at this edge gives rise to normal fluctuations $\mathcal{A}_i^-(\Delta \mathcal{U})$ and $\mathcal{A}_i^+(\Delta \mathcal{U})$ contributing to cells L and R respectively. These terms lead to the computation of transverse fluctuations giving contribution to neighboring cells across edges (j) and (k) for cell L, and across edges (m) and (l) for cell R. Transverse fluctuations are computed by projecting normal fluctuations on the characteristic basis associated to the Riemann problem (21) defined on the adjacent edge, hence the name transverse Riemann solver. Spectral analysis of corresponding Jacobian matrix is carried out in [27]. The negative normal fluctuation is, for instance, decomposed on the characteristic basis associated to edge (j) as:

$$ \mathcal{A}_i^-(\Delta \mathcal{U}) = \sum_{m=1}^{M} \beta_m \mathcal{R}_j^m = \sum_{m=1}^{M} \mathcal{W}_j^m \quad (35) $$

where $\mathcal{R}_j^m$ is based on the normal vector $N^j$ but also on different tangent moduli between grid cells L and its neighbor T. The calculation of transverse fluctuations is performed by solving equation (35) for coefficients $\beta_m$, and by multiplying each wave $\mathcal{W}_j^m$ by the associated characteristic speed. Since only waves with positive celerities with respect to the orientation defined by outward normal vector to considered edge will contribute to the transverse fluctuation, only the positive operator $\mathcal{B}^+$ is used:

$$ \mathcal{B}_j^+ \mathcal{A}_i^-(\Delta \mathcal{U}) = \sum_{m=1}^{M} \lambda_m \beta_m \mathcal{R}_j^m \quad (36) $$

An additional numerical flux defined at edges is hence built from these transverse fluctuations:

$$ \mathcal{F}_j^{\text{tran}} = \frac{\Delta t}{2\Delta X_j} \mathcal{B}^+ \mathcal{A}^-(\Delta \mathcal{U}) \quad (37) $$
which contributes to the grid cell T (ΔX<sub>j</sub> being the length of edge (j)). With duplicated nodes introduced by the DG approximation this contribution must be counted negatively for nodes belonging to cell L (outgoing fluctuation) and positively for nodes belonging to cell T (incoming fluctuation).

### 4.5. Boundary conditions enforcement

As in MPM [21, 28], boundary conditions (BCs) are treated at nodes that compose the domain boundary (see figure 1). Introducing ghost nodes on boundary interfaces allows to widen the use of the Approximate Riemann solver. Similarly to finite volumes [8], a vector Q<sub>G</sub> is enforced on those nodes so that the stationary solutions of Riemann problems (22) are consistent with boundary conditions. This vector is determined for one ghost node by (i) setting the free quantities (i.e: velocity (resp. stress) components for Neumann (resp. Dirichlet) BCs) equal to that of the associated internal node and (ii) solving system (32) for enforced quantities on ghost nodes.

Notice that the use of conservation laws (21) requires to inverse constitutive equations in order to enforce Neumann boundary conditions, which is not convenient. Hence, the use of the quasi-linear form (22) also allows to easily handle boundary conditions. Moreover, the transverse Riemann solver needs to introduce corner ghost nodes equivalent to finite volume corner cells [8].

## 5. Numerical procedure

Let us assume that the vector of specific conserved quantities $\tilde{\mathcal{U}}^n$ is known at every material points that discretize a continuum body $\Omega$ in a grid made of $N_{nodes}$ nodes at a given time $t^n$. We can now derive the computational steps to obtain the DGMPM solution, which requires to (i) solve the hyperbolic system (5) in $\Omega$ with associated boundary conditions and (ii) update the vector $\tilde{\mathcal{U}}^n$ on material points.

The scheme has been established with a total Lagrangian formulation therefore, the weak form (16) is written on the reference configuration. Hence, the lumped mass matrix $M_i^L$ and the pseudo-stiffness matrices $K_{ij}^\alpha$ are computed once and for all at the beginning of the calculation:

$$M_i^L = \sum_j \left( \sum_p S_{ip} m_p S_{jp} \right) ; \quad K_{ij}^\alpha = \sum_p \frac{\partial S_{ip}}{\partial X_\alpha} m_p S_{jp} \quad i, j = 1, \dots, N_{nodes} \quad (38)$$

Then, the procedure follows these steps:

The known at material points

(a) Convective phase: the discrete equation (20) and Riemann problems (22) require a projection of fields onto the grid:

$$M_i^L \tilde{\mathcal{U}}_i = \sum_{p=1}^{p} S_{ip} m_p \tilde{\mathcal{U}}_p \quad \text{and} \quad M_i^L \mathcal{Q}_i = \sum_{p=1}^{p} S_{ip} m_p \mathcal{Q}_p \quad (i = 1, \dots, N_{nodes}) \quad (39)$$

to be solved for each $\tilde{\mathcal{U}}_i$ and $\mathcal{Q}_i$ respectively. This weighted least squares interpolation is similar the one introduced in FLIP (Fluid Implicit Particle method) to avoid a diffusion-like error [1, 29, 30].

(b) Compute time step: in the general case the tangent modulus depends on the deformation gradient as well as the waves speeds. It is therefore needed to compute the time step so that the fastest wave can at most cross the smallest cell of the mesh according to Courant condition.

(c) The specific flux vectors $\tilde{\mathcal{F}}_{\alpha,i}$ involved in the discrete equation (20) are computed from $\mathcal{Q}_i$ knowing $\rho_0$, thus avoiding the computation of constitutive equations.

(d) Enforce boundary conditions on mesh nodes (see section 4.5).

(e) Computation of interface fluxes:

1. Build $\mathcal{Q}_{\xi^*}$ states based on $\mathcal{Q}_i^n$ where i denotes the nodes belonging to the face on side L or R.
2- Compute the stationary solution by means of the approximate Riemann solver.

3- Calculate the corresponding normal flux.

(f) Advance solution in time by solving the discrete equation (20) at each node.

(g) Map the updated solution to material points with a classical interpolation:

$$ \bar{u}_p^{n+1} = \sum_{i=1}^{N} s_{ip} \bar{u}_i^{n+1} \qquad (40) $$

At this point the mesh has virtually moved, but since fields have been transferred back to particles, the underlying grid can be discarded and rebuilt for the next time step for computational convenience, thus involving the convective phase (a) at the next time step. Notice that adaptive grid algorithms can be applied in the reference configuration, so that to improve wave front tracking in the current one.

(h) Material point kinematics and constitutive model: The new solution $\bar{u}_p^{n+1}$ allows to increment the deformation $\varphi(X, t)$ and to update stress components which will be used in the auxiliary vector for the next time step, through hyperelastic constitutive equations:

$$ \varphi_p^{n+1} = X_p + \Delta t \nu_p^{n+1} \qquad (41) $$

$$ \Pi_p^{n+1} = \frac{\partial \Psi}{\partial F}(F_p^{n+1}) \qquad (42) $$

This summary highlights significant differences with respect to the MPM. First, while the DGMPM uses a classical interpolation in order to map back the updated solution to material points (step f), FLIP method and MPM require an additional time integration on material points. Indeed, in the original schemes the changes in grid values (i.e: accelerations and strain rates) are mapped back to particles rather than the new value in order to avoid numerical dissipation [30]. Furthermore, the use of conservation laws (5) instead of momentum equation in the weak form implies that both velocity and gradients are solved at nodes making this new approach close to finite volume methods. Finally, the deformation gradient is no longer calculated with shape functions gradients so that the task of choosing between *Update Stress First* (USF) and *Update Stress Last* (USL) algorithms [31, 32] vanishes. The USL formulation consists in integrating the constitutive equations at material points with the updated nodal velocity field given by the resolution of the MPM discrete equation thanks to shape functions gradients. Within the USF algorithm, the material points stress tensors are computed with the nodal velocity field resulting from the convective phase at the beginning of each time step. Since the stress tensors at particles are used in the discrete equation through the internal forces vector, these two formulations modify the accuracy of the MPM. In that sense the DGMPM scheme is simpler though additional cost is linked to Riemann solutions. Fortunately, the use of discontinuous Galerkin approximation makes this numerical method highly parallelizable [14].

# 6. Stability analysis for one-dimensional linear problems

A stability analysis of the one-dimensional DGMPM scheme is performed on a scalar linear hyperbolic equation. Following [33], the discretized equation is written in a finite difference sense in order to perform a von Neumann linear stability analysis.

## 6.1. *Model equation - Space discretization*

We consider the scalar linear advection equation for an arbitrary quantity $q = \rho \bar{q}$ moving at the constant speed $a \in \mathbb{R}^{+*}$ in a homogeneous one-dimensional medium of length $l$:

$$ \frac{\partial q}{\partial t} + a \frac{\partial q}{\partial X} = 0 \qquad (43) $$
The specific flux function is $\bar{\mathcal{F}} = a\bar{q}$ and equation (43) is discretized with the discontinuous Galerkin material point method. Thus, it is assumed that the medium has been divided with P material points arbitrarily distributed in E two-nodes elements of constant length $\Delta X$ (figure 4). The grid mesh is such that at least one particle lies in every cell during the computation in order to ensure that there is no hole in the bar. Moreover, periodic boundary conditions are considered to simplify the analysis.

<figure><img src="image_5.png" /><figcaption>Figure 4: One-dimensional mesh made of E elements of constant length $\Delta X = \frac{1}{E}$.</figcaption></figure>

Since fields are carried by particles, we seek for the scheme equation that gives the solution $\bar{Q}_{\alpha}^{n+1}$ at material point $\alpha$ and time step $n + 1$, with respect to the solutions of other particles at previous time step. In this section, Latin and Greek indices denote respectively nodes and material points, $S_{i\alpha}$ is thus the shape function of node $i$ evaluated at the position of material points $\alpha$. Linear shape functions defined in element $k$ yield:

$$S_{2k-1}(X) = \frac{X_{2k} - X}{\Delta X} \qquad S_{2k}(X) = \frac{X - X_{2k-1}}{\Delta X} \qquad X \in [X_{2k-1}, X_{2k}] \qquad (44)$$

Furthermore, the cell containing a given particle $\beta$ will be denoted by $c(\beta)$ and this notation will be also used for nodes, for instance: $c(2k - 1) = c(2k) = k$.

## 6.2. DGMPM scheme equation

The method followed to write the scheme equation is to trace backward the numerical procedure described in section 5 in order to get an expression for material point $\alpha$:

$$\bar{Q}_{\alpha}^{n+1} = f(\bar{Q}_{\mu}^{n}) \quad \mu = 1, \dots, P \qquad (45)$$

where $f$ is a function representing the DGMPM scheme. Quantities at time $t^{n+1}$ are obtained by interpolating nodal solutions of the discrete equation (20) thanks to mapping (17):

$$\bar{Q}_{\alpha}^{n+1} = \sum_{i=1}^{E} (S_{2i-1,\alpha}\bar{q}_{2i-1}^{n+1} + S_{2i,\alpha}\bar{q}_{2i}^{n+1}) \qquad (46)$$

Interfaces fluxes are $\mathcal{F}_N = (aq^*)N$, where $q^*$ is the stationary solution of Riemann problem at cell interface and $N = \pm 1$ the outward unit normal. Nodal solutions at time $n + 1$ in equation (46) are the result of the discrete system (20):

$$\bar{q}_{i}^{n+1} = \bar{q}_{i}^{n} + \frac{\Delta t}{M_{i}^{L}} a \left( \sum_{j=1}^{2E} K_{i,j} \bar{q}_{j}^{n} - q_{i}^{*} N_{i} \right) \quad (\text{no sum on } i) \qquad (47)$$

Provided linear shape functions, the lumped mass and the pseudo-stiffness matrices are:

$$M_{i}^{L} = \sum_{\beta=1}^{P} S_{i\beta} m_{\beta} \qquad (48)$$

$$K_{2i-1,j} = \sum_{\lambda=1}^{P} \frac{\partial S_{2i-1,\lambda}}{\partial X} m_{\lambda} S_{j\lambda} = - \sum_{\lambda=1}^{P} \frac{m_{\lambda} S_{j\lambda}}{\Delta X} \qquad (49)$$

$$K_{2i,j} = \sum_{\lambda=1}^{P} \frac{\partial S_{2i,\lambda}}{\partial X} m_{\lambda} S_{j\lambda} = \sum_{\lambda=1}^{P} \frac{m_{\lambda} S_{j\lambda}}{\Delta X} \qquad (50)$$
where $m_\alpha$ is the mass of $\alpha^{th}$ material point. Since the pseudo-stiffness matrix component $K_{ij}$ is zero if $i$ and $j$ are not connected by an element, due to the discontinuities of shape functions, the sum over $j$ in equation (47) reduces to: $\sum_{j=1}^{2E} K_{2i-1,j} \bar{q}_j^n = K_{2i-1,2i-1} \bar{q}_{2i-1}^n + K_{2i-1,2i} \bar{q}_{2i}^n$. We can thus write:

$$ \frac{\sum_j K_{2i-1,j} \bar{q}_j^n}{M_{2i-1}^L} = - \frac{\sum_\lambda (S_{2i-1,\lambda} \bar{q}_{2i-1}^n + S_{2i,\lambda} \bar{q}_{2i}^n)}{\Delta X \sum_\beta S_{2i-1,\beta}} ; \quad \frac{\sum_j K_{2i,j} \bar{q}_j^n}{M_{2i}^L} = \frac{\sum_\lambda (S_{2i-1,\lambda} \bar{q}_{2i-1}^n + S_{2i,\lambda} \bar{q}_{2i}^n)}{\Delta X \sum_\beta S_{2i,\beta}} \quad (51) $$

Since waves travel from left to right ($a > 0$), the stationary solution is equal to $q_L$ on an interface, that is:

$$ q_{2i-1}^* = \rho \bar{q}_{2i-2}^n ; \quad q_{2i}^* = \rho \bar{q}_{2i}^n \quad (52) $$

$$ N_{2i-1} = -1 ; \quad N_{2i} = 1 \quad (53) $$

Moreover, in a homogeneous medium the mass density can be replaced by $\rho = N^{c(\alpha)}m_\alpha/\Delta X$ where $N^{c(\alpha)}$ is the number of particles in the cell that contains the $\alpha^{th}$ material point. Then, we have:

$$ \frac{q_{2i-1}^*}{M_{2i-1}^L} = \frac{N^{c(\alpha)}}{\Delta X \sum_{\beta} S_{2i-1,\beta}} \bar{q}_{2i-2}^n ; \quad \frac{q_{2i}^*}{M_{2i}^L} = \frac{N^{c(\alpha)}}{\Delta X \sum_{\beta} S_{2i,\beta}} \bar{q}_{2i}^n \quad (54) $$

Equations (51) and (54) are introduced into discrete system (47) so that equation (46) is rewritten as:

$$ \bar{Q}_{\alpha}^{n+1} = \sum_{i=1}^{2E} S_{i\alpha} \bar{q}_i^n + a \frac{\Delta t}{\Delta X} \sum_{i=1}^{E} \sum_{\lambda=1}^{P} \left( S_{2i-1,\lambda} \bar{q}_{2i-1}^n + S_{2i,\lambda} \bar{q}_{2i}^n \right) \left[ \frac{S_{2i,\alpha}}{\sum_{\beta} S_{2i,\beta}} - \frac{S_{2i-1,\alpha}}{\sum_{\beta} S_{2i-1,\beta}} \right] \\ + a \frac{\Delta t}{\Delta X} N^{c(\alpha)} \sum_{i=1}^{E} \left[ \frac{S_{2i-1,\alpha}}{\sum_{\beta} S_{2i-1,\beta}} \bar{q}_{2i-2}^n - \frac{S_{2i,\alpha}}{\sum_{\beta} S_{2i,\beta}} \bar{q}_{2i}^n \right] \quad (55) $$

In equation (55) the solutions at nodes result from the convection step (39):

$$ \bar{q}_i^n = \frac{\sum_{\mu} S_{i\mu} m_{\mu} \bar{Q}_{\mu}^n}{\sum_{\gamma} S_{i\gamma} m_{\gamma}} = \frac{\sum_{\mu} S_{i\mu} \bar{Q}_{\mu}^n}{\sum_{\gamma} S_{i\gamma}} \quad (56) $$

Thus, introduction of mapping (56) in equation (55) and permutation of sums over $\mu$ and $i$ lead after some simplifications to the scheme equation:

$$ \bar{Q}_{\alpha}^{n+1} = \sum_{\mu=1}^{P} \bar{Q}_{\mu}^{n} \left\{ \sum_{i=1}^{2E} S_{i\alpha} \frac{S_{i\mu}}{\sum_{\gamma} S_{i\gamma}} + a \frac{\Delta t}{\Delta X} \sum_{i=1}^{E} (S_{2i-1,\mu} + S_{2i,\mu}) \left[ \frac{S_{2i,\alpha}}{\sum_{\beta} S_{2i,\beta}} - \frac{S_{2i-1,\alpha}}{\sum_{\beta} S_{2i-1,\beta}} \right] \right. \\ \left. + a \frac{\Delta t}{\Delta X} N^{c(\alpha)} \sum_{i=1}^{E} \left[ \frac{S_{2i-1,\alpha}}{\sum_{\beta} S_{2i-1,\beta}} \frac{S_{2i-2,\mu}}{\sum_{\gamma} S_{2i-2,\gamma}} - \frac{S_{2i,\alpha}}{\sum_{\beta} S_{2i,\beta}} \frac{S_{2i,\mu}}{\sum_{\gamma} S_{2i,\gamma}} \right] \right\} \quad (57) $$

This recurrence relation can be further simplified by viewing that the term in parentheses is non zero and equal to one, due to partition of unity, only for element $i = c(\mu)$. On the other hand, $S_{2i-1,\alpha}S_{2i-2,\mu}$ is non zero only if cell $i$ contains material point $\alpha$ and is such that $i = c(\mu) + 1$. Hence, the scheme equation reads:

$$ \bar{Q}_{\alpha}^{n+1} = \sum_{\mu=1}^{P} \bar{Q}_{\mu}^{n} \left\{ \sum_{i=1}^{2E} S_{i\alpha} \frac{S_{i\mu}}{\sum_{\gamma} S_{i\gamma}} + a \frac{\Delta t}{\Delta X} \left[ \frac{S_{2c(\mu),\alpha}}{\sum_{\beta} S_{2c(\mu),\beta}} - \frac{S_{2c(\mu)-1,\alpha}}{\sum_{\beta} S_{2c(\mu)-1,\beta}} \right] \right. \\ \left. + a \frac{\Delta t}{\Delta X} N^{c(\alpha)} \left[ \frac{S_{2c(\mu)+1,\alpha}}{\sum_{\beta} S_{2c(\mu)+1,\beta}} \frac{S_{2c(\mu),\mu}}{\sum_{\gamma} S_{2c(\mu),\gamma}} - \frac{S_{2c(\mu),\alpha} S_{2c(\mu),\mu}}{(\sum_{\beta} S_{2c(\mu),\beta})^2} \right] \right\} \quad (58) $$

which can be written for simplicity:

$$ \bar{Q}_{\alpha}^{n+1} = \sum_{\mu=1}^{P} \bar{Q}_{\mu}^{n} D_{\alpha\mu} \quad (59) $$
The numerical solution $\tilde{Q}_{\alpha}^{n+1}$ depends on the Courant number $a\Delta t/\Delta X$ and on the position of material points within the grid mesh through the shape functions in the expression of coefficients $D_{\alpha\mu}$. Furthermore, the first (resp. second) brackets in equation (58) involve shape functions that are non zero if material point $\mu$ and $\alpha$ lie in the same cell (resp. previous cell). This property means that the numerical domain of dependence of the DGMPM for scalar linear advection with $a > 0$ covers two cells regardless of the number of material points.

### 6.3. The von Neumann linear stability analysis

The computational domain is now repeated periodically by mapping it to the domain $[-l, 0]$ so that the solution at material point $\beta$ and time step $n$ is expanded into a discrete Fourier basis of $2E + 1$ harmonics over the domain $X \in [-l, l]$:

$$ \tilde{Q}_{\beta}^{n} = \sum_{j=-E}^{E} A_{j}^{n} e^{I\beta k_{j} \Delta X} \qquad (60) $$

with $A_j^n$, the magnitude of the jth harmonic at time step $n$, $I = \sqrt{-1}$, and $k_j$ a wave number. Introduction of this expansion in equation (59) yields:

$$ A_j^{n+1} e^{I \alpha k_j \Delta X} = \sum_{\mu=1}^{P} A_j^n D_{\alpha\mu} e^{I \mu k_j \Delta X} \quad \forall j = -E, \dots, E \qquad (61) $$

The amplification factor between two time steps at a given point is defined as:

$$ \frac{A_j^{n+1}}{A_j^n} = \sum_{\mu=1}^{P} e^{I(\mu-\alpha)k_j \Delta X} D_{\alpha\mu} \quad \forall j = -E, \dots, E \qquad (62) $$

A necessary condition to ensure the stability of a numerical scheme is that the amplification factor must be lower than one in modulus: $|A^{n+1}/A^n| < 1$. This upper bound allows to ensure that eventual errors do not increase during the computation. For expression (62), this leads to:

$$ \left| \sum_{\mu=1}^{P} e^{I(\mu-\alpha)k_j \Delta X} D_{\alpha\mu} \right| \le \sum_{\mu=1}^{P} \left| e^{I(\mu-\alpha)k_j \Delta X} D_{\alpha\mu} \right| = \sum_{\mu=1}^{P} |D_{\alpha\mu}| \qquad (63) $$

where the triangle inequality, and the unit modulus of the complex number $e^{I(\mu-\alpha)k_j \Delta X}$ have been used. Hence, the DGMPM scheme is stable for a given discretization if the Courant number is set so that the following condition is satisfied:

$$ \sum_{\mu=1}^{P} |D_{\alpha\mu}| \le 1 \qquad (64) $$

According to the scheme equation (58), such a stability condition can be very hard to deal with analytically for general discretizations. However, the particular configuration for which one single particle lies in each cell can be studied analytically. In this case, condition (64) reads:

$$ \sum_{\mu=1}^{P} |D_{\alpha\mu}| = \left|1 - a \frac{\Delta t}{\Delta X}\right| + \left|a \frac{\Delta t}{\Delta X}\right| \le 1 \qquad (65) $$

which is satisfied for all $a \frac{\Delta t}{\Delta X} \le 1$ regardless of the positions of particles. Hence, this space discretization leads to a stable scheme for the optimal Courant condition while the classical DGFEM scheme developed in [22] is restricted to condition $\Delta t/\Delta X = O(\sqrt{\Delta X})$. Recall that this DGMPM space discretization can be viewed as a DGFEM with reduced integration of the weak form (see section 3). Hence, combination of reduced integration with projections between mesh nodes and integration points allows to improve the restrictive stability condition of original DGFEM [22, 34]. Note that by considering space-time DG formulations [35, 36], it is possible to relax constraints of pure space DGFEM and obtain a critical CFL at 1. This optimal CFL condition for DGMPM is however no longer valid for discretizations
involving more than one material point per cell. Indeed, figure 5 shows the evolution of the amplification factor with
respect to the Courant number for several (regular) discretizations. Since the scheme equation (58) depends on the
particle considered, those curves are plotted for each material point contained in a given cell. In figures 5(e) and
5(f) we can see that the amplification factor satisfies the stability condition (64) for CFL numbers lower than one.
Two material points in each cell are considered in figure 5(e), with either a symmetrical or an offset configuration,

<figure><img src="image_6.png" /><figcaption>Figure 5: Evolution of amplification factor with Courant number for several space discretizations</figcaption></figure>

respectively shown in figures 5(a) and 5(b). It also appears that the more restrictive condition on CFL number is
given by the upwind (first) material point (CFL = 0.43 for symmetrical configuration and CFL = 0.34 for shifted
one). Moreover, the material points positions influence the stability condition as shown by the shifted configuration
for which the result is more restrictive. Similar observations can be made with figure 5(f) in which three particles
are considered. In that case, Courant number decreases to CFL = 0.29 for configuration 5(c) and CFL = 0.26 for
configuration 5(d).

The von Neumann linear stability analysis performed in this section provides an explicit (nonlinear) condition (64) that must be satisfied in order to ensure the stability of the DGMPM scheme. First, this relation allows to fully exploit the ability of the method to rebuild the grid mesh arbitrarily. Indeed, after such a procedure, the number and positions of material points in grid cells can change and one must be able to adapt properly the CFL number in order to ensure stability. An advantage on the original MPM is hence highlighted since no stability condition exists for the latter scheme. Second, it has been shown that some DGMPM discretizations provide an improvement of the restrictive CFL number that applies to the DGFEM scheme. In particular, the optimal condition allowing to capture discontinuities (CFL = 1) can be reached when one particle lies in each cell of the grid mesh. This property turns out to be a strength for the DGMPM since it aims at follow waves accurately as shown with the following numerical results.

7. Numerical results

In this section, comparisons between DGMPM, MPM and FEM are performed. The first simulation considers
one-dimensional linear elasticity and is purely used for analysis since MPM has little interest for small strains. Next,
finite deformations are considered for one-dimensional and two-dimensional problems.

7.1. Compression impact on a linear elastic bar

We consider a bar made of a linear elastic material with density $ρ_0 = 7800kg.m^{-3}$, Young's modulus $E = 2×10^{11}Pa$ and arbitrary cross section $S$, suddenly subjected to a compression stress on its left end $T^d = 1×10^9 Pa$ while the right one is fixed. We assume that the impact is such that the linearized geometrical framework holds and mass density is constant $ρ(x,t) = ρ_0(X)$. The analytical solution of this problem is known [37] and consists of an elastic disturbance propagating rightward in the bar at constant speed $c = \sqrt{E/ρ}$ before reflection on the right end. The length of the bar is $l = 6m$ and it is discretized with 150 regular grid cells containing one single centered material point. A comparison is performed between the DGMPM, the MPM and analytical solution. The MPM has been used for this comparison with the USL and USF implementation in order to see which procedure gives the most suitable results for impact problems. The Cauchy stress profile along the bar given by the three numerical schemes are plotted in figure 6 at two different time steps. Figure 6(a) and 6(b) show respectively the incident and reflected waves. The DGMPM solution

<figure><img src="image_7.png" /><figcaption>Figure 6: Cauchy stress profile along the bar at different times with single-particle per cell discretization. Comparison between DGMPM ($CFL = 1$), MPM-USF ($CFL = 1$), MPM-USL ($CFL = 0.5$), and analytical solution [37].</figcaption></figure>

fits perfectly the analytical profile due to the condition $CFL = 1$ while the MPM-USF, which suffers from diffusion and oscillations, can not capture a discontinuity though the same Courant number is used. On the other hand, the USL formulation of the material point method is unstable for $CFL ≥ 0.7$, although this result is not proved analytically. The scheme is then unable to capture discontinuities. However, the oscillations resulting from the MPM-USL procedure decrease with time and the correct stress level is obtained faster than for the USF solution. Furthermore, numerical diffusion is quite similar for the two MPM algorithms. Hence, the USL implementation of the material point method seems to be better suited to impact problems than the USF.

It is of interest to see the influence of the number of material points lying in each element. Thus, we consider now the same problem as before, with 150 material points regularly spaced in 75 cells as space discretization. The particles are placed symmetrically with respect to cell centers (figure 7). In that case, the optimal stability condition no longer holds. A comparison is performed between the same numerical methods as before and the analytical solution. In

<figure><img src="image_8.png" /><figcaption>Figure 7: A bar discretized with two-particles per cell discretization.</figcaption></figure>

addition, a second order in time DGMPM version is used by means of a second order Runge-Kutta time discretization
(RK2-DGMPM). This particular implementation leads for this space discretization to a Courant number set to one, though this result has not been proved analytically. On the other hand, the CFL numbers used for the material point methods are unchanged whereas the first order in time DGMPM one is decreased (CFL = 0.5). The decrease of the

<figure><img src="image_9.png" /><figcaption>Figure 8: Cauchy stress profile along the bar at different times with two-particles discretization. Comparison between DGMPM with forward Euler time discretization (CFL = 0.5), RK2-DGMPM (CFL = 1), MPM-USF (CFL = 1), MPM-USL (CFL = 0.5) and analytical solution.</figcaption></figure>

CFL number introduces more diffusion than in the MPM solutions but no spurious oscillations appear. On the other hand, MPM solutions are quite similar to that of the previous discretization, though the USF algorithm leads to more oscillations. Furthermore, figure 8 shows that using a second order time integrator within the DGMPM scheme allows to obtain the optimal stability condition and hence, to fit the analytical solution. However, this property no longer holds when material points are moved in each cell so that the symmetry with respect to cells centers disappears. Thus, a condition equivalent to (64) is required for second order time discretization. This will be the object of future works.

## 7.2. Impact problems on a hyperelastic one-dimensional medium

In this section we consider a one-dimensional medium made of a hyperelastic Saint-Venant-Kirchhoff (SVK) material in order to see how the discontinuous Galerkin material point method behaves for nonlinear hyperbolic problems. The following deformation gradient tensor occurs in this medium:

$$F = F e_1 \otimes e_1 + e_2 \otimes e_2 + e_3 \otimes e_3 \qquad (66)$$

The Saint-Venant-Kirchhoff constitutive equation leads to the following PK1 tensor (see Appendix A):

$$\Pi = \Pi e_1 \otimes e_1 + \Pi_2 (e_2 \otimes e_2 + e_3 \otimes e_3) \qquad (67)$$

where $\Pi = \frac{2\mu+\lambda}{2}F(F^2 - 1)$, $\Pi_2 = \frac{\lambda}{2}(F^2 - 1)$ and $(\mu, \lambda)$ are the Lamé's coefficients. The SVK model is known to not be polyconvex however, the problem treated remains hyperbolic if condition (A.10), which depends on the loads, holds. Hence, we will only consider loading conditions that satisfy hyperbolicity condition. Moreover, the flux function resulting from this constitutive model is concave ($-\frac{\partial^2\Pi}{\partial F^2} < 0$) and leads to monotonically increasing characteristic speeds with the deformation gradient. Since the evolution of wave celerities with respect to the deformation gradient governs the wave pattern (i.e: either a rarefaction or a shock wave), the structure of the solution will differ from other constitutive laws with convex flux function such as the neo-Hookean model. As a consequence, the Saint-Venant-Kirchhoff model can lead to non-physical solutions. Nevertheless, the use of SVK model in spite of its limitations is here justified since an analytical solution can be found (see Appendix A), which is not possible for a neo-Hookean material.
### 7.2.1. Compression load

The one-dimensional medium is first submitted to a compression load on its left end $T^d = 2 \times 10^{10} Pa$ so that the solution of this problem consists of a right-going rarefaction wave as shown in Appendix A. The domain is discretized with 150 material points lying in 150 uniform grid cells and a comparison between DGMPM and MPM is performed. Figure 9 shows the stress profiles and the analytical solution along a horizontal line at different time steps. Slight oscillations appear in the DGMPM solution at the beginning of the computation (figure 9(a)). However, these oscillations are much lower than the MPM ones. We then see that the DGMPM stress profile is close to the analytical solution but does not fit it exactly. Indeed, the approximate Riemann solver does not account for the

<figure><img src="image_10.png" /><figcaption>Figure 9: First Piola-Kirchhoff stress profile along the bar at different times solution of the compression impact problem. Comparison between DGMPM ($CFL = 1$), MPM ($CFL = 0.5$) and analytical solution (Appendix A).</figcaption></figure>

information contained in the rarefaction fan and cannot provide the exact solution. Figure 10 shows the error made on Riemann problems at cell interfaces, characterized by the distance between the intersection of integral curves on the one hand, and the intersection of approximated curves on the other hand. The 1-integral curve (resp. 2-integral curve) corresponds to all the states $\mathcal{U}$ connected through a shock wave (resp. rarefaction wave) to left (resp. right) initial condition at the left end of the bar. The stationary state corresponds to the intersection of those curves in the phase plane ($F, \rho_0 v$) and we see that the approximate and analytical stationary solutions are different. One can imagine to use an exact Riemann solver based on the analytical solution in order to reduce this error. However, such an implementation requires to solve a nonlinear problem at each cell interface and is very costly. Moreover, an analytical solution is available for SVK material but it is not the case for other constitutive models, which prevents the use of an exact solver.

### 7.2.2. Tensile load

We now consider a tensile load of same magnitude that gives rise to a shock for which the celerity of the fastest wave is given by the deformation gradient on the upwind side. The MPM solution of this problem shows the same limitations as before, it is hence not presented. The Figure 11 shows results obtained for the same discretization as before with the DGMPM compared to the analytical solution. We see that slight oscillations appear in the DGMPM solution after the shock. Note that those oscillations (figure 11(a)) decrease with time (figure 11(b)). However, these oscillations can be removed with a Courant number smaller than one. Moreover, the numerical scheme do not perfectly capture the discontinuity, though a good behavior is shown, due to the wave celerity used for the time step evaluation which is an upper bound of the shock speed $\omega_s$ (A.13). Nevertheless, the loading conditions considered so far are deliberately extreme for the purpose of visualization. A lower tensile load will lead to a numerical solution closer to analytical results as can be seen in figure 12 which shows the comparison for a load magnitude $T^d = 2 \times 10^{10} Pa$.

## 7.3. Impact problem on a hyperelastic two-dimensional medium

Let's now move to multi-dimensional simulations by considering an infinite medium in directions $e_2$ and $e_3$, and of width $l$ in direction $e_1$ (see figure 13), undergoing a tensile load. The domain is modeled by a finite medium, with
<figure><img src="image_11.png" /><figcaption>Figure 10: Analytical and approximate stationary solutions corresponding to the problem of compression impact on the left end of a SVK hyperelastic one-dimensional medium in the phase plane $(F, \rho_0 v)$. Initial conditions $u_L$ and $u_R$ are such that the analytical stationary solution corresponds to the applied load.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: First Piola-Kirchhoff stress profile along a horizontal line of the one-dimensional medium at different times. Comparison between DGMPM CFL = 1 and analytical solution.</figcaption></figure>

zero shear stress and transverse velocity components prescribed on the top, bottom and right boundaries (figure 13). In addition zero out-of-plane velocity and strain components are assumed so that plane waves can be simulated. The solid is discretized with 61 × 11 material points such that only one particle lies in each grid cell. The simulation has been performed with the DGMPM-DCU, since BCs reproduce a plane wave, which is compared to the analytical solution. Figures 14(a) and 14(b) show the comparison between numerical and analytical results for longitudinal and transversal stresses $\Pi_{11}$ and $\Pi_{22}$, with a tensile load set to $T^d = 4 \times 10^9$ Pa. First, top figures in subplots show the longitudinal stress isovalues within the medium. The results have been post-treated with the software *Paraview* [38] which filters allow to represent material points with spheres deliberately big for the purposes of visualization. We see in those stress maps that the one-dimensional dependency of the solution is not disturbed within the numerical scheme. Next, bottom graphs in subplots 14(a) and 14(b) show analytical and numerical stress components along the bottom boundary of the medium. As for finite volume and finite element methods, two-dimensional simulations require a decrease of Courant number with the DGMPM. Since no two-dimensional stability analysis is available so far, the CFL number was set to 0.7 by testing various values. Hence, the shock wave is described less accurately than for one-dimensional computations. However, no oscillations appear in the numerical solutions and the analytical
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12: First Piola-Kirchhoff stress profile along the bar at different times solution of the tensile impact problem for a surface force magnitude $T^d = 4 \times 10^8$ Pa. Comparison between DGMPM CFL = 1 and analytical solution.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Two-dimensional medium modeling.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Two-dimensional simulation of an impact problem on a hyperelastic medium. Comparison between DGMPM-DCU (CFL = 0.7) and analytical solutions.</figcaption></figure>

stress levels are obtained for both components.

### 7.4. Impact problem on a hyperelastic plate

Previous simulations showed that the DGMPM provides good results for one-dimensional problems compared to an analytical solution. In this section, an extension to more general two-dimensional problems is made, for which no analytical solution is available, thus comparison with a finite element solution (bilinear approximation) is performed. The two-dimensional plate is made of the compressible hyperelastic neo-Hookean material, which is known to be a polyconvex model, and considered here in the plane strain case. Geometry, boundary and loading conditions, and
material parameters are given in figure 15(a). The solid is discretized such that DGMPM material points are equivalent

<figure><img src="image_17.png" /><figcaption>Figure 15: Impact on a two-dimensional plate made of a neo-Hookean hyperelastic material. Parameters : $E = 2 \times 10^{11}$ Pa, $\nu = 0.3$, $\rho = 7800$ kg $\cdot$ m$^{-3}$, $\nu^d = 1000$ m $\cdot$ s$^{-1}$.</figcaption></figure>

to finite element nodes, thus the plate is represented with $l \times h \equiv 33 \times 25$ material points and $l \times h \equiv 32 \times 24$ Q1 elements for FEM (see figure 15(b) for a coarse example). Once again, only one material point lies in each cell of the arbitrary computational grid. The finite element computation is performed with the software Abaqus [39], using an explicit time discretization with no artificial viscosity added. These numerical results are compared to those obtained from DGMPM computations using both CTU and DCU methods. The Courant number is set to unity in DGMPM-CTU and to 0.7 in DGMPM-DCU leading to *average* time steps $\Delta t_{CTU} = 1.57 \times 10^{-5}$ s and $\Delta t_{DCU} = 1.13 \times 10^{-5}$ s, whereas the *constant* time step used in the FEM simulation is $\Delta t_{FEM} = 1.43 \times 10^{-5}$ s. Figure 16 shows numerical results in terms of the Cauchy stress tensor $\sigma = \frac{1}{J} \Pi F^T$ exported from Abaqus to Paraview with the code developed in [40]. Numerical Cauchy stress components ($\sigma_{11}, \sigma_{22}$) along bottom boundary are gathered in the left plots and longitudinal stress $\sigma_{11}$ isovalues in the domain resulting from CTU and FEM computations are compared in right colored maps for various time steps. Note that the deformed shape of finite element method is used as a background for the CTU stress map in order to compare shapes at different time steps. At the beginning of the computation (figure 16(a)), stress profiles are quite similar despite slight oscillations in the FEM solutions and more numerical diffusion in the DCU solution than in the CTU one. Then, interactions between shear and tension-compression waves occurring in the left part of the plate give rise to a decrease of the stress level on the wave front which is managed differently by the two methods. Thus, additional oscillations are introduced in the FEM solution leading to a more diffuse wave profile (right colored map in figure 16(b)). On the other hand, the deformed shapes of the plate computed with the CTU and the FEM remain close. After the reflection of waves on the fixed boundary, significant differences can be seen in numerical solutions (figure 16(c)). We also see that final shapes differ, which can be explained by the velocity fields plotted in figure 17.

Those curves show the horizontal component of velocity along bottom boundary at different times. First, figures 17(a) and 17(b) show that spurious oscillations also appear in FEM velocity solution before reflection of the wave on fixed boundary. Those oscillations remain in the solution until the end of the computation (figure 17(d)). On the other hand, the prescribed velocity boundary condition at the right end of the plate is not enforced exactly in the DGMPM solutions as shown in figure 17(d). Indeed, the enforcement of boundary conditions in MPM is still a challenging question [28] which does not spare the DGMPM. However, this point has not been investigated yet. Nonetheless, boundary conditions do not disturb much waves profiles as seen with all the simulations presented above. Hence, the final DGMPM stresses plotted in figure 16(d), as well as previous ones, do not exhibit overshoots. The absence of such over-stresses in hyperelastic regime is expected to allow a better management of history-dependent constitutive models including threshold such as elastoplasticity. In those models, non physical overshoots can lead to a wrong assessment of residual stresses and strains.
<figure><img src="image_18.png" /><figcaption>Figure 16: Cauchy stress tensor axial components along bottom boundary of the plate and longitudinal component in the whole domain. Comparison between FEM, DGMPM-CTU and DGMPM-DCU solutions. Velocity is prescribed on the bottom part of the left boundary.</figcaption></figure>

## 8. Conclusion

A new discontinuous Galerkin material point method was presented for the numerical simulation of impacts on hyperelastic media undergoing finite deformations. This approach is based on the weak form of a system of conservation laws written in a solid domain represented by particles in an arbitrary computational grid, as in the original material point method. The background mesh is used as a support for the discontinuous Galerkin approximation leading to a weak form element-by-element that involves numerical fluxes at cell faces. These are computed from an approximate Riemann solver leading to the donor-cell-upwind method. In addition, transverse corrections can be added for multidimensional problems so as to improve accuracy of the scheme by adapting the corner-transport-upwind method. The space discretization is based on particles mass density and an explicit time discretization is performed leading to the discrete equations. The stability analysis of this numerical scheme has been carried out in section 6 leading to an explicit condition on CFL number with respect to the discretization used. Moreover, it was shown that when one single particle lies in each cell of the grid mesh for one-dimensional problems, the optimal Courant number is allowed. This result highlights that the DGMPM stability condition is improved compared to that of the DGFEM. The enhanced stability combined with the use of the characteristic structure of the solution lead to a scheme that provides the ability to accurately follow waves in media undergoing finite deformations under total Lagrangian approach. Comparisons with an analytical one-dimensional solution in sections 7.2 and 7.3, and with the finite element method for two-dimensional problems in section 7.4 confirmed the good behavior of the method.

The approach can also be adapted to arbitrary high order approximations. However, well-suited polynomials com-
<figure><img src="image_19.png" /><figcaption>Figure 17: Horizontal velocity component along bottom boundary of the plate. Comparison between FEM, DGMPM-CTU and DGMPM-DCU solutions. Velocity is prescribed on the bottom part of the left boundary.</figcaption></figure>

bined with reconstruction based quadrature such as Moving Least Squares or Spline interpolation will be required in order to avoid numerical instabilities due to the nonphysical negative nodal masses [41, 42]. Moreover, the arbitrariness of the grid can be fully exploited by employing adaptive mesh techniques on the reference configuration in order to track accurately waves in the current one for problems involving complex geometries.

## Appendix A. Analytical solution of the one-dimensional problem

We focus here on the analytical structure of the solution of the hyperbolic system written in the conservative form (5) for one-dimensional strain states considering a hyperelastic Saint-Venant-Kirchhoff material (SVK). In that case, the deformation gradient and the PK1 tensor are respectively:

$$F = F e_1 \otimes e_1 + e_2 \otimes e_2 + e_3 \otimes e_3 \qquad (A.1)$$

$$\Pi = \Pi_{11} e_1 \otimes e_1 + \Pi_{22} (e_2 \otimes e_2 + e_3 \otimes e_3) \qquad (A.2)$$

which corresponds to a *plane wave* solution. The stress tensor is given by the following SVK constitutive equations:

$$\Pi = F\Sigma = F[2\mu E + \lambda \text{trace}(E)I] \qquad (A.3)$$

where $\lambda, \mu$ are the Lamé's coefficients, $\Sigma$ is the second Piola-Kirchhoff stress tensor and $E$ the Green-Lagrange strain tensor defined as:

$$E = \frac{1}{2}(F^T F - I) \qquad (A.4)$$
These definitions for the one-dimensional strain problem lead to the stress components:

$$ \Pi_{11} = \frac{2\mu + \lambda}{2} F(F^2 - 1) \quad ; \quad \Pi_{22} = \frac{\lambda}{2}(F^2 - 1) = \frac{\lambda}{F(2\mu + \lambda)} \Pi_{11} \qquad (A.5) $$

Thus, only the axial stress will be considered since the other component can be deduced from it. For simplicity in what follows, the stress will read:

$$ \Pi = \frac{\lambda + 2\mu}{2} F(F^2 - 1) \qquad (A.6) $$

The Jacobian matrix is then:

$$ \mathcal{J}_N(\tilde{\mathcal{U}}) = \frac{\partial \mathcal{F}}{\partial \mathcal{U}} = - \begin{bmatrix} 0 & \frac{1}{\rho_0} \\ \frac{\lambda+2\mu}{2}(3F^2-1) & 0 \end{bmatrix} \quad \text{with} \quad \mathcal{U} = \begin{bmatrix} F \\ \rho_0 v \end{bmatrix}, \quad \mathcal{F} = \begin{bmatrix} -v \\ -\Pi \end{bmatrix} \qquad (A.7) $$

with eigen-values and associated right eigen-vectors:

$$ \omega^1 = - \sqrt{\frac{\lambda + 2\mu}{2\rho_0}(3F^2 - 1)} , \quad \mathcal{R}^1 = \begin{bmatrix} 1 \\ -\rho_0 \omega^1 \end{bmatrix} \quad ; \quad \omega^2 = \sqrt{\frac{\lambda + 2\mu}{2\rho_0}(3F^2 - 1)} , \quad \mathcal{R}^2 = \begin{bmatrix} 1 \\ -\rho_0 \omega^2 \end{bmatrix} \qquad (A.8) $$

The two characteristic fields of the system are genuinely nonlinear since the following condition is satisfied for each field:

$$ \frac{\partial \omega^i(\mathcal{U})}{\partial \mathcal{U}} \cdot \mathcal{R}^i \neq 0 \quad \forall \mathcal{U} \in \mathbb{R}^2 \qquad (A.9) $$

Hyperbolicity of the system is ensured by real eigen-values, leading to the following condition on the deformation gradient:

$$ \omega^i \in \mathbb{R} \Leftrightarrow F \geq 1/\sqrt{3} \qquad (A.10) $$

strict hyperbolicity being ensured with the strict inequality. Moreover, the dependency on deformation gradient of characteristic speed can lead to two types of wave in the solution. Indeed, considering the Riemann problem:

$$ \frac{\partial \mathcal{U}}{\partial t} + \frac{\partial \mathcal{F}}{\partial X} = 0, \quad X \in ]-\infty, \infty[, \quad t > 0 \\ \mathcal{U}(X, 0) = \begin{cases} \mathcal{U}_L & \text{if } X < 0 \\ \mathcal{U}_R & \text{if } X > 0 \end{cases} \qquad (A.11) $$

If $\mathcal{U}(X, 0)$ is such that $F_L > F_R$, then $\omega_L^1 > \omega_R^1$ and the first characteristic field corresponds to a centered rarefaction wave meaning that the characteristics move away from each other. Furthermore, $\omega_L^2 > \omega_R^1$ implies that the second characteristic fields is related to a shock wave in which the characteristics collide and generate a shock (see figure A.18). The same reasoning leads to an opposite configuration for the case $F_L < F_R$ (see figure A.19).

<figure><img src="image_20.png" /><figcaption>Figure A.18: Wave pattern generated by Riemann problem initial condition $F_L > F_R$: 1-rarefaction, 2-shock</figcaption></figure>

The two structures of the solution (1-rarefaction, 2-shock or 1-shock, 2-rarefaction) are similar to the wave pattern that result from the *Dam-Break problem* of shallow water equations. Hence, the analytical solution of the hyperelastic bar impact problem for the SVK material is close to the dam-break solution developed in [8].
<figure><img src="image_21.png" /><figcaption>Figure A.19: Wave pattern generated by Riemann problem initial condition $F_L < F_R$: 1-shock, 2-rarefaction</figcaption></figure>

## Shock waves

Two states $\mathcal{U}$ lying on each side of a shock wave are connected by the so-called Rankine-Hugoniot condition:

$$s_i[\mathcal{U}] = [\mathcal{F}] \qquad (A.12)$$

where $[\bullet]$ denotes a jump across the $i^{th}$ discontinuous wave traveling at speed $s_i$. The jump condition must be satisfied between a known state $\mathcal{U}$ and an unknown state on the other side of the shock $\mathcal{U}_*$:

$$s_i(F - F_*) = -(v - v_*) \qquad (A.13)$$

$$s_i(\rho_0 v - \rho_0 v_*) = -\frac{\lambda + 2\mu}{2} [F^3 - F - (F_*^3 - F_*)] \qquad (A.14)$$

Substitution of the shock speed in the second equation with its expression given by the first one leads to a second order polynomial in $v_*$:

$$(v - v_*)^2 = \frac{\lambda + 2\mu}{2\rho_0} (F - F_*) [F^3 - F - (F_*^3 - F_*)] \qquad (A.15)$$

It is easy to verify that the Lax entropy condition which prevents from an expansion shock ($\omega(\mathcal{U}_L) > s_i > \omega(\mathcal{U}_R)$), ensures that the discriminant of this polynomial is positive and the solution $v_*$ then reads:

$$v_* = v \pm \sqrt{\frac{\lambda + 2\mu}{2\rho_0}(F - F_*)[F^3 - F - (F_*^3 - F_*)]} \qquad (A.16)$$

corresponding to the two families of shock. Finally, by assuming a weak shock ($\mathcal{U} \approx \mathcal{U}_*$) one gets the conditions to be satisfied:

$$\rho_0 v_* = \rho_0 v_L + \sqrt{\rho_0 \frac{\lambda + 2\mu}{2} (F_L - F_*) [F_L^3 - F_L - (F_*^3 - F_*)]} \quad \text{across a 1-shock} \qquad (A.17)$$

$$\rho_0 v_* = \rho_0 v_R - \sqrt{\rho_0 \frac{\lambda + 2\mu}{2} (F_R - F_*) [F_R^3 - F_R - (F_*^3 - F_*)]} \quad \text{across a 2-shock} \qquad (A.18)$$

## Centered rarefaction waves

A centered rarefaction wave is a self-similar solution defined in genuinely nonlinear fields such that:

$$\mathcal{U}(X, t) = \tilde{\mathcal{U}}(\xi = X/t) \qquad (A.19)$$

where $\tilde{\mathcal{U}}(\xi = X/t)$ traces out a curve which tangent vector at each point is collinear to a given eigen-vector $\mathcal{R}^i$. Such a smooth curve is an integral curve defined by $\tilde{\mathcal{U}}'(\xi) = \alpha(\xi)\mathcal{R}^i(\tilde{\mathcal{U}}(\xi))$. With the parameter $\xi = X/t$ defined before, one shows that the solution across a rarefaction wave must satisfy:

$$\tilde{\mathcal{U}}'(\xi) = \frac{\mathcal{R}^i(\tilde{\mathcal{U}}(\xi))}{\nabla \omega^i(\tilde{\mathcal{U}}) \mathcal{R}^i} \quad \text{for} \quad \omega^i(\mathcal{U}_L) \le \xi \le \omega^i(\mathcal{U}_R) \qquad (A.20)$$
where the denominator is nonzero in a genuinely nonlinear field due to condition (A.9). The integration between $\xi = \omega^i(\mathcal{U}_L)$ and $\xi = \omega^i(\mathcal{U}_R)$ of this system of ordinary differential equations, after some algebra, leads to the expression of integral curves. These expression can be used to compute a state $\mathcal{U}_*$ connected to a known state $\mathcal{U}$ through a centered rarefaction wave:

$$ \rho_0 v_* = \rho_0 v_L + \sqrt{\rho_0 \frac{\lambda + 2\mu}{24}} \left[ \sqrt{3} \left( F_* \sqrt{3F_*^2 - 1} - F_L \sqrt{3F_L^2 - 1} \right) - \ln \frac{\sqrt{3}F_* + \sqrt{3F_*^2 - 1}}{\sqrt{3}F_L + \sqrt{3F_L^2 - 1}} \right] \quad \text{for a 1-rarefaction} \quad (A.21) $$

$$ \rho_0 v_* = \rho_0 v_R - \sqrt{\rho_0 \frac{\lambda + 2\mu}{24}} \left[ \sqrt{3} \left( F_* \sqrt{3F_*^2 - 1} - F_R \sqrt{3F_R^2 - 1} \right) - \ln \frac{\sqrt{3}F_* + \sqrt{3F_*^2 - 1}}{\sqrt{3}F_R + \sqrt{3F_R^2 - 1}} \right] \quad \text{for a 2-rarefaction} \quad (A.22) $$

Finally, solutions inside the rarefaction fan are given by $\xi = \pm\omega(F(\xi))$ and therefore:

$$ F(\xi) = \sqrt{\frac{2\rho_0}{3(\lambda + 2\mu)}} \xi^2 + 1 \quad (A.23) $$

and $\rho_0 v(\xi)$ is computed with either equation (A.21) or (A.22) by taking $F_* = F(\xi)$.

It is then possible to compute the solution for a 1-shock,2-rarefaction characteristic fields (resp. 1-rarefaction,2-shock) by solving the nonlinear system composed of equations (A.17) and (A.22) (resp. equations (A.21) and (A.18)). The solution inside the rarefaction fan being given by equation (A.23).
[1] D. Sulsky, Z. Chen, H. L. Schreyer, A particle method for history-dependent materials, Computer methods in applied mechanics and engineering 118 (1-2) (1994) 179-196.

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

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

[4] V. Psyk, D. Risch, B. Kinsey, A. Tekkaya, M. Kleiner, Electromagnetic forminga review, Journal of Materials Processing Technology 211 (5) (2011) 787-829.

[5] T. Belytschko, W. K. Liu, B. Moran, K. Elkhodary, Nonlinear finite elements for continua and structures, John wiley & sons, 2013.

[6] J. C. Simo, T. J. Hughes, Computational inelasticity, Vol. 7, Springer Science & Business Media, 2006.

[7] Z. Więkowski, The material point method in large strain engineering problems, Computer methods in applied mechanics and engineering 193 (39) (2004) 4417-4438.

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

[9] A. Harten, High resolution schemes for hyperbolic conservation laws, Journal of computational physics 49 (3) (1983) 357-393.

[10] N. Favrie, S. Gavrilyuk, Dynamics of shock waves in elastic-plastic solids, in: ESAIM: Proceedings, Vol. 33, EDP Sciences, 2011, pp. 50-67.

[11] T. Heuzé, Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids, Journal of Computational Physics 346 (2017) 369-388.

[12] C. H. Lee, A. J. Gil, J. Bonet, Development of a cell centred upwind finite volume algorithm for a new conservation law formulation in structural dynamics, Computers & Structures 118 (2013) 13-38.

[13] J. Haider, C. H. Lee, A. J. Gil, J. Bonet, A first-order hyperbolic framework for large strain computational solid dynamics: An upwind cell centred total lagrangian scheme, International Journal for Numerical Methods in Engineering 109 (3) (2017) 407-456.

[14] B. Cockburn, Discontinuous galerkin methods for convection-dominated problems, in: High-order methods for computational physics, Springer, 1999, pp. 69-224.

[15] G. Jacobs, J. S. Hesthaven, High-order nodal discontinuous galerkin particle-in-cell method on unstructured grids, Journal of Computational Physics 214 (1) (2006) 96-121.

[16] B. J. Plohr, D. H. Sharp, A conservative eulerian formulation of the equations for elastic flow, Advances in Applied Mathematics 9 (4) (1988) 481-499.

[17] J. M. Ball, Convexity conditions and existence theorems in nonlinear elasticity, Archive for rational mechanics and Analysis 63 (4) (1976) 337-403.

[18] D. A. Di Pietro, A. Ern, Mathematical aspects of discontinuous Galerkin methods, Vol. 69, Springer Science & Business Media, 2011.

[19] X. Feng, O. Karakashian, Y. Xing, Recent Developments in Discontinuous Galerkin Finite Element Methods for Partial Differential Equations: 2012 John H Barrett Memorial Lectures, Vol. 157, Springer Science & Business Media, 2013.

[20] J. A. Trangenstein, Numerical solution of hyperbolic partial differential equations, Cambridge University Press, 2009.

[21] E. Love, D. Sulsky, An unconditionally stable, energymomentum consistent implementation of the material-point method, Computer Methods in Applied Mechanics and Engineering 195 (3336) (2006) 3903-3925.

[22] G. Chavent, G. Salzano, A finite-element method for the 1-d water flooding problem with gravity, Journal of Computational Physics 45 (3) (1982) 307-344.

[23] S. Osher, Riemann solvers, the entropy condition, and difference, SIAM Journal on Numerical Analysis 21 (2) (1984) 217-235.

[24] J. Trangenstein, P. Colella, A higher-order godunov method for modeling finite deformation in elastic-plastic solids, Communications on Pure and Applied Mathematics 44 (1991) 41-100.

[25] G. Kluth, B. Després, Perfect plasticity and hyperelastic models for isotropic materials, Continuum Mechanics and Thermodynamics 20 (3) (2008) 173-192.

[26] P. Colella, Multidimensional upwind methods for hyperbolic conservation laws, Journal of Computational Physics 87 (1) (1990) 171-200.

[27] G. Kluth, Analyse mathmatique et numrique de systmes hyperlastiques et introduction de la plasticité, Ph.D. thesis (2008).

[28] M. Steffen, P. Wallstedt, J. Guilkey, R. Kirby, M. Berzins, Examination and analysis of implementation choices within the material point method (MPM), Computer Modeling in Engineering and Sciences 32 (2) (2008) 107-127.

[29] D. Burgess, D. Sulsky, J. Brackbill, Mass matrix formulation of the flip particle-in-cell method, Journal of Computational Physics 103 (1) (1992) 1-15.

[30] J. Brackbill, H. Ruppel, Flip: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions, Journal of Computational Physics 65 (2) (1986) 314-343.

[31] O. Buzzi, D. M. Pedroso, A. Giacomini, et al., Caveats on the implementation of the generalized material point method, Computer Modeling in Engineering and Sciences 1 (1) (2008) 1-21.

[32] X. Zhang, Z. Chen, Y. Liu, The Material Point Method: A Continuum-based Particle Method for Extreme Loading Cases, Academic Press, 2016.

[33] C. Hirsch, Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics, Butterworth-Heinemann, 2007.

[34] W. H. Reed, T. Hill, Triangular mesh methods for the neutron transport equation, Los Alamos Report LA-UR-73-479.

[35] G. R. Richter, An explicit finite element method for the wave equation, Applied Numerical Mathematics 16 (1) (1994) 65-80.

[36] R. Abedi, B. Petracovici, R. Haber, A spacetime discontinuous galerkin method for linearized elastodynamics with element-wise momentum balance, Computer Methods in Applied Mechanics and Engineering 195 (25) (2006) 3247-3273, discontinuous Galerkin Methods.

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

[38] U. Ayachit, The paraview guide: a parallel visualization application (2015).

[39] M. Smith, ABAQUS/Standard User's Manual, Version 6.9, Simulia, 2009.

[40] Q. Liu, J. Li, J. Liu, Paraview visualization of abaqus output on the mechanical deformation of complex microstructures, Comput. Geosci. 99 (C) (2017) 135-144. doi:10.1016/j.cageo.2016.11.008.
URL https://doi.org/10.1016/j.cageo.2016.11.008

[41] R. Tielen, High-order material point method, Master's thesis, Delft University of Technology (2016).

[42] R. Tielen, E. Wobbes, M. Möller, L. Beuth, A high order material point method, Procedia Engineering 175 (2017) 265-272.