# Implicit Material Point Method for non-associated plasticity of geomaterials

Louis Guillet, Vincent Acary, Franck Bourrier, Oliver Goury

► To cite this version:

Louis Guillet, Vincent Acary, Franck Bourrier, Oliver Goury. Implicit Material Point Method for non-associated plasticity of geomaterials. 2025. hal-05070887

HAL Id: hal-05070887

https://hal.science/hal-05070887v1

Preprint submitted on 16 May 2025

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
Implicit Material Point Method for non-associated plasticity of
geomaterials

Louis Guillet¹, Vincent Acary¹, Franck Bourrier², and Oliver Goury¹

¹Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, Institute of Engineering, LJK, 38000, Grenoble, {louis.gillet,vincent.acary,olivier.goury}@inria.fr

²Univ. Grenoble Alpes, INRAE, CNRS, IRD, Grenoble INP, IGE, 38000 Grenoble, France, franck.bourrier@inrae.fr

**Abstract**

This paper presents an implicit solver for non-associated plasticity problems based on the semi-smooth Newton method, using a formulation of a second order cone complementarity problem. The presented method is used to solve the incremental problem derived from the Material Point Method (MPM). Geomaterials are often described by non-associated elasto-plastic models, that are usually simulated with the classic return mapping algorithm. The goal of this work is to propose an alternative approach, developing a monolithic solver capable of solving precisely the plastic constitutive equations and the Dirichlet boundary conditions, expressed as equality constraints, while avoiding, in particular, the difficult computation of the consistent tangent operator. This kind of method can be adapted to include unilateral contact, material hardening and Coulomb friction in the monolithic solver. Inspired by the Implicit Standard Material (ISM) framework, the constitutive equations are written as differential inclusion and a semi-smooth Newton's method is used to solve the constitutive laws within an implicit time scheme. The capacities of the method are illustrated on classical problems in the field of geomechanics, namely the footing and slope stability problems, and compared to classical results from the literature, in particular analytical solutions based on limit analysis, showing good agreement.

**keywords** : Non-smooth solver, Implicit material point method, non conforming Dirichlet boundary condition, non-associated plasticity

**Acknowledgement**

We would like to acknowledge Inria and the exploratory action GRANIER, which have provided funding for the work conducted in the scope of this paper.
# Contents

<table><tr><td>1 Introduction</td><td>3</td></tr><tr><td>2 Modelling the dynamics of elastoplastic associated and non-associated materials</td><td>5</td></tr><tr><td>3 Numerical Method : Material Point Method</td><td>11</td></tr><tr><td>4 Solving of the incremental problem</td><td>15</td></tr><tr><td>4.1 Solving of the complementarity conditions</td><td>16</td></tr><tr><td>4.2 Dirichlet Boundary Conditions</td><td>17</td></tr><tr><td>4.3 Non smooth Newton's method</td><td>19</td></tr><tr><td>5 Results and discussion</td><td>20</td></tr><tr><td>5.1 Footing</td><td>20</td></tr><tr><td>5.2 Slope stability</td><td>32</td></tr><tr><td>6 Conclusion</td><td>38</td></tr><tr><td>A Yield criterion in plane strain plasticity</td><td>43</td></tr><tr><td>B Retrieving a SOC using Voigt's notation</td><td>44</td></tr><tr><td>C Sub-gradient of natural map</td><td>45</td></tr></table>
# 1 Introduction

The objective of this work is to develop a robust, stable and implicit numerical framework, in continuity of the work of [1], to simulate non-associated elasto-plastic material undergoing large deformations. This framework, based on the writing of variational inequalities, is aimed to be compatible with different discretisations, namely the Material Point Method (MPM) and the Finite Element Method (FEM).

The MPM is a particle based method dedicated to simulate continuous media. It has been firstly introduced in [2] and is a combination of two existing methods, the Particle in Cell (PIC) developed in [3] and the Fluid Interpolated Particles (FLIP) developed in [4]. It relies on two discretizations, an Eulerian grid and Lagrangian material points (also referred as particles in the literature). The computation is performed on the grid nodes and the information is stored on the particles. This allows to get the benefit of Lagrangian methods while avoiding any mesh distortion. This method is well suited for large deformation problems and can be associated with numerous material constitutive models. In addition to its suitability for parallel computing, the simplicity of its formulation, its adaptation for both small and large strains and its capacity to handle complex constitutive models have pushed many researchers to use this method in the last two decades. It has been used in the field of geosciences, for modelling soils deformations [5], snow avalanches, [6, 7, 8], glacier creeping [9, 10], or in the field of computer graphics mainly for granular media simulation [11, 12]. Compared with the Finite Element Method (FEM), this method has been the subject of fewer theoretical studies. However, several studies [13, 14] have shown the good energy conservation properties of the implicit version of this method for elasto-plastic materials subjected to large deformations, and others [15, 16] have focused on error estimation in MPM methods. Numerous studies [17, 12] have relied on numerical experiments to focus on the stability or convergence of different interpolation schemes within the scope of MPM.

Most studies based on the MPM use explicit time integrators. It allows the method to handle very large amount of material points (more than hundreds of million), thereby facilitating the investigation of large scale phenomena. However, explicit time integrators can suffer from lack of stability and do not accurately solve at each time step the governing equations (since there is decoupling between the conservation equation and the constitutive relation). Other studies use an implicit time integrator to gain stability and therefore handle larger time steps (see for example [18, 19, 20]). In the previous references, the elasto-plastic behaviour of the material is solved with the classical Return Mapping Algorithm (RMA) [21, 22]. Even if it suffers from several limitations, namely difficulties in the computation of the so-called consistent tangent operator and in the convergence for any initial guess, the quadratic local convergence allows to solve the constitutive laws with machine accuracy.

In the last two decades, some studies [23, 24] used tools from convex cone programming to solve variational inequalities for associated materials. These tools often combine a FEM discretisation with commercial solvers for convex optimization based on the Interior Point Method. More recently, this kind of approach has been used with a MPM, or with a MPM-FEM coupling [25]. In this type of approach, the non-associated character of the material, if considered, is treated in an uncoupled manner.

This paper presents an approach for non-associated elasto-plasticity based on non-smooth optimisation. Plasticity is presented as a variational inequality, as a prolongation of the work described in [1]. This reformulation of the plastic law is not new [26, 27, 28] and associated plasticity is written within the scope
of the Generalised standard materials (GSM). With a change of variable, non-associated plasticity can also
be derived from a variational formulation with the introduction of the Implicit Standard material (ISM)
described in [29]. In the latter formulation, plasticity is written in the same way as Coulomb Friction, with
the advantage of benefiting from the numerous studies on implicit solvers done for this problem [30, 31, 11].
Choosing a solver often comes down to striking a balance between robustness and convergence rate [31].
For the problem of frictional contact, the structure of the discrete system often imposes the use of a first
order method such as fix point or a Gauss-Siedel type of solver [32]. In our framework, the elasticity of the
material induces a discrete problem with full rank, direct second order method [33] as Newton's method are
used. It is worth noticing that the presented method does not differ much from the classic RMA. It presents
wider formulation where RMA can appear as a particular case of the semi-smooth Newton's method [22]. In
addition, this approach does not require the explicit computation of the consistent tangent operator. The
main contributions of this work are :

* the formulation of the elasto-plastic incremental problem for Material Point Method as a variational principle

• the formulation of the Newton’s Method for non-associated elasto-plasticity

* the comparison with the FEM for classical case studies in the field of geomechanics, namely footing and slope stability problems
## 2 Modelling the dynamics of elastoplastic associated and non-associated materials

The constitutive model described in this work is a linear elasto-plastic model with perfect plasticity. The notation used in this work are inspired from [34]. We introduce the following material quantities: $\rho$ the density field, $v$ the velocity field, $u$ the displacement field, $\sigma$ the internal stress tensor field, $f_v$ the external volumetric forces.

The material belongs to a set $\Omega \in \mathbb{R}^d$ where $d$ corresponds to the dimension. The material domain $\Omega$ is moving through time as the material deforms. We have the following conservation equations:

$$ \frac{D\rho}{Dt} + \rho\nabla \cdot v = 0, \qquad (1) $$

$$ \rho \frac{Dv}{Dt} = \nabla \cdot \sigma + f_v. \qquad (2) $$

Under the small strain assumption, the Euler-Almansi tensor $e$, classicaly used when choosing an Eulerian framework, can be approximated by:

$$ e \sim \varepsilon = \frac{1}{2}(\nabla u + (\nabla u)^T). \qquad (3) $$

The classical assumption of additive decomposition of the strain tensor into an elastic par $\varepsilon^E$ and a plastic part $\varepsilon^P$ is also made:

$$ \varepsilon = \varepsilon^E + \varepsilon^P. \qquad (4) $$

The definition of an elastic potential energy $\psi$ gives a relation between the elastic part of the strain tensor and the stress tensor $\sigma$:

$$ \sigma = \frac{\partial\psi(\varepsilon^E)}{\partial\varepsilon^E}. \qquad (5) $$

In the case of linear elasticity, the elastic potential energy is a quadratic form and the relation between $\varepsilon^E$ and $\sigma$ is linear:

$$ \sigma = \mathcal{D} : \varepsilon^E. \qquad (6) $$

### Associated Plasticity : Generalised Standard Material

In this section an elasto-plastic model following the framework of Generalised Standard Material (GSM) is presented [27]. As a first step, we do not include hardening in this model, although it can also be done [1]. The intrinsic dissipation per unit of volume $d$ is

$$ d = \sigma : \dot{\varepsilon} - \dot{\psi}(\varepsilon^E). \qquad (7) $$

The chain rule and the additive decomposition allow to rewrite the intrinsic dissipation as $d = \sigma : \dot{\varepsilon}^P$. The Clausius-Duhem inequality, equivalent to the second principle of thermodynamics, ensures that the intrinsic dissipation remains always non-negative.

The plasticity condition is described by a plasticity criterion or a yield criterion $f$ ensuring that $f$ remains non negative. The set

$$ \mathcal{C} := \{\sigma \text{ s.t. } f(\sigma) \le 0\}. \qquad (8) $$
is called the admissible set. The only assumption made on $f$ is that $\mathcal{C}$ is a convex set including the origin. In the case of associated plasticity, following the work of [27, 28, 35], the plastic deformation rate is described by the normal inclusion

$$ \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} \in \mathcal{N}_{\mathcal{C}}(\boldsymbol{\sigma}). \qquad (9) $$

In the latter equation, $\mathcal{N}_{\mathcal{C}}(\boldsymbol{\sigma})$ represents the normal cone to the set $\mathcal{C}$ at the point $\boldsymbol{\sigma}$ and is defined by

$$ \mathcal{N}_{\mathcal{C}}(\boldsymbol{\sigma}) = \{y \text{ s.t. } \forall \hat{\boldsymbol{\sigma}} \in \mathcal{C}, y : (\boldsymbol{\sigma} - \hat{\boldsymbol{\sigma}}) \geq 0\}. \qquad (10) $$

The relation between the plastic strain increment $\dot{\boldsymbol{\epsilon}}^{\mathrm{P}}$ and $\boldsymbol{\sigma}$ is called the plastic flow rule. Note that the inclusion in (9) induces the maximum dissipation principle as

$$ d = \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} : \boldsymbol{\sigma} \geq \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} : \hat{\boldsymbol{\sigma}} \quad \forall \hat{\boldsymbol{\sigma}} \in \mathcal{C}. \qquad (11) $$

For a convex set defined by a yield criterion as in (8), the definition of the normal cone yields the standard flow rule for plasticity as

$$ \left\{ \begin{array}{l} \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} = \nabla f(\boldsymbol{\sigma}) \lambda \\ 0 \geq f(\boldsymbol{\sigma}) \perp \lambda \geq 0, \end{array} \right. \qquad (12) $$

where $\lambda$ are the plastic multipliers.

The normality condition written in (9) ensures the existence of two conjugate convex potentials $V$ and $W$ such that for all $\hat{\boldsymbol{\sigma}}$ and $\hat{\boldsymbol{\epsilon}}^{\mathrm{P}}$ admissible

$$ V(\hat{\boldsymbol{\sigma}}) + W(\hat{\boldsymbol{\epsilon}}^{\mathrm{P}}) \geq \hat{\boldsymbol{\sigma}} : \hat{\boldsymbol{\epsilon}}^{\mathrm{P}}, \qquad (13) $$

and, at the solution, we have

$$ V(\boldsymbol{\sigma}) + W(\dot{\boldsymbol{\epsilon}}^{\mathrm{P}}) = \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}}^{\mathrm{P}}, \qquad (14) $$

which is equivalent to

$$ \boldsymbol{\sigma} \in \partial_{\hat{\boldsymbol{\epsilon}}^{\mathrm{P}}} W(\dot{\boldsymbol{\epsilon}}^{\mathrm{P}}) \text{ and } \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} \in \partial_{\boldsymbol{\sigma}} V(\boldsymbol{\sigma}). \qquad (15) $$

This formulation, combined with other assumptions on the free energy $\psi$, ensures the well-posedness of the problem, see [27] for more details.

## Non-associated plasticity : Implicit Standard Materials

The normality rule described in (9) is characteristic of associated materials and does not allow to accurately model the volume changes that strongly condition the response of geomaterials (sands, soils,...). It is thus common in the field of geomechanics to introduce a relaxation of the normality rule. To this aim, the concept of Implicit Standard Material (ISM) was first introduced in the work of [29]. It defines a broader theoretical framework for non-associated materials by replacing the sum of the two potential $V + W$ by a function $b$ of both variables $\dot{\boldsymbol{\epsilon}}^{\mathrm{P}}$ and $\boldsymbol{\sigma}$ that satisfies the following properties :

* for all $\boldsymbol{\sigma}$ admissible fixed, the function $b(\boldsymbol{\sigma}, \bullet)$ is convex,

* for all $\dot{\boldsymbol{\epsilon}}^{\mathrm{P}}$ admissible fixed, the function $b(\bullet, \dot{\boldsymbol{\epsilon}}^{\mathrm{P}})$ is convex,

* for all $(\hat{\boldsymbol{\sigma}}, \hat{\boldsymbol{\epsilon}}^{\mathrm{P}})$ admissible, $b(\hat{\boldsymbol{\sigma}}, \hat{\boldsymbol{\epsilon}}) \geq \hat{\boldsymbol{\sigma}} : \hat{\boldsymbol{\epsilon}}$,
* at the optimum, we have the equivalence : $b(\sigma, \dot{\epsilon}^P) = \sigma : \dot{\epsilon}^P$ if and only if $\dot{\epsilon}^P \in \partial_\sigma b(\sigma, \dot{\epsilon}^P)$ and $\sigma \in \partial_{\dot{\epsilon}^P} b(\sigma, \dot{\epsilon}^P)$.

Such function b is called a bi-potential. The GSM model is a particular case of the ISM model with $b = V + W$.

## Yield criterion : Drucker-Prager model

In this research work, we choose to use the Drucker-Prager yield criterion as it is extensively used to model both associated and non-associated materials, and in particular geomaterials. The criterion, originally developed in [36], is expressed as a function of the first invariant $I_1$ of the stress tensor and of the second invariant $J_2$ of the deviatoric part of the stress tensor called s. These invariants are computed as follows :

$$I_1(\boldsymbol{\sigma}) = \operatorname{tr}(\boldsymbol{\sigma}) = \sum_i \sigma_{ii}, \quad \boldsymbol{s} = \boldsymbol{\sigma} - \frac{I_1(\boldsymbol{\sigma})}{3} \boldsymbol{I}, \quad J_2(\boldsymbol{\sigma}) = \frac{1}{2} \boldsymbol{s} : \boldsymbol{s}. \qquad (16)$$

When the dimension d is equal to 3, the invariant $J_2(\boldsymbol{\sigma})$ can be rewritten as a function of the components of $\boldsymbol{\sigma}$ as :

$$J_2(\boldsymbol{\sigma}) = \frac{1}{6} \sum_{i,j} (\sigma_{ii} - \sigma_{jj})^2 + (\sigma_{12}^2 + \sigma_{13}^2 + \sigma_{23}^2). \qquad (17)$$

Considering two parameters $\eta$ and $k$, depending only on the material properties, the yield criterion of the Drucker-Prager model is

$$f(\boldsymbol{\sigma}) = \sqrt{J_2(\boldsymbol{\sigma}) + \eta I_1(\boldsymbol{\sigma}) - k}, \qquad (18)$$

where the parameters $\eta$ and $k$ are expressed depending on the internal frictional angle and the cohesion allowing thus to relate Drucker-Prager and Mohr-Coulomb criteria, see [37, 36, 35] for more details. The admissible set $\mathcal{C} = \{\boldsymbol{\sigma} \text{ s.t. } f(\boldsymbol{\sigma}) \le 0\}$ for the stress tensor is a translated cone due to cohesion. As illustrated in Figure 1, we have

$$\boldsymbol{\sigma} \in \mathcal{C} \iff \frac{k}{\eta} \boldsymbol{e} + \boldsymbol{\sigma} \in \boldsymbol{K} \qquad (19)$$

where $\boldsymbol{K}$ is the second cone defined as $\boldsymbol{K} = \{\boldsymbol{\sigma} \text{ s.t. } -\eta I_1(\boldsymbol{\sigma}) \ge \sqrt{J_2(\boldsymbol{\sigma})}\}$ and $\boldsymbol{e}$ writes as

$$\boldsymbol{e} = -\frac{1}{d} \boldsymbol{I}. \qquad (20)$$

The admissible set for the plastic strain rate can also be expressed in terms of invariant. Let us focus first on the associated plasticity : the normality rule written in (9) is equivalent to

$$\dot{\boldsymbol{\epsilon}}^P \in \mathcal{N}_K \left( \boldsymbol{\sigma} + \frac{k}{\eta} \boldsymbol{e} \right), \qquad (21)$$

which is equivalent to the cone complementarity problem

$$\boldsymbol{K}^* \ni -\dot{\boldsymbol{\epsilon}}^P \perp \left( \boldsymbol{\sigma} + \frac{k}{\eta} \boldsymbol{e} \right) \in \boldsymbol{K}. \qquad (22)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic of the yield criterion and of the associated plastic flow rule. In blue is represented the convex set $C$ of admissible stress and in orange the corresponding second order cone $K$</figcaption></figure>

The complementarity written in (22) ensures that $-\dot{\varepsilon}^P \in K^*$ the dual cone of $K$. This is equivalent in terms of invariant as

$$ \sqrt{J_2(\dot{\varepsilon}^P)} - \frac{1}{\eta} I_1(\dot{\varepsilon}^P) \leq 0. \qquad (23) $$

**Remark 1.** There is no effect of the cohesion on the definition of the admissible set for the plastic strain rate. This is due to the fact that the cohesion within the term $k$ will just translate the admissible set, keeping the normality rule unchanged
Dilatancy control : non-associated Drucker-Prager model

<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic of the yield criterion with the non-associated flow rule and the corresponding change of variables.</figcaption></figure>

In the case of non-associated plasticity, the inclusion described in (9) is no more valid. The new admissible set $K_\epsilon$ for the plastic strain rate will depend on an additional coefficient $\eta_\theta$ accounting for dilatancy

$$ \sqrt{J_2(\dot{\epsilon}^P)} - \frac{k_d}{\eta_\theta} I_1(\dot{\epsilon}^P) \le 0, \qquad (24) $$

where $k_d$ is a normalisation coefficient depending on the dimension. The coefficient $\eta_\theta$ and the admissible set $K_\epsilon$ are illustrated in Figure 2. In practice, the coefficient $\eta_\theta$ will depend on $\theta$ in the same manner as $\eta$ depends on $\phi$, $\theta$ being the dilatancy angle defined as in [29] or [38].

When the set $C$ is a cone, as it is the case for the Drucker Prager model, the normal cone inclusion is a cone complementarity problems, which is a well studied problem in optimization with numerous numerical solvers. The objective is to retrieve a normal cone inclusion in our plastic model. Following the work of [29], this is possible to define new variables as

$$ \tilde{\epsilon}^P = \dot{\epsilon}^P + (\eta - \eta_\theta) \sqrt{J_2(\dot{\epsilon}^P)} I. \qquad (25) $$

With this new variable, we obtain the following normal cone inclusion

$$ \tilde{\epsilon}^P \in \mathcal{N}_C(\sigma). \qquad (26) $$

The change of variable is illustrated with the schematic represented in Figure 2. By introducing $\mu = \eta - \eta_\theta$, the change of variable in (25) has a similar formulation as the complementarity problem for dry
Coulomb friction presented in [31] with $\mu$ being a friction coefficient. The following system summarises the constitutive model and the governing equations

$$ \left\{ \begin{array}{l} \rho \frac{Dv}{Dt} = \nabla \cdot \boldsymbol{\sigma} + f_v \\ \boldsymbol{\sigma} = \mathcal{D} : \nabla_s(u) - \dot{\boldsymbol{\varepsilon}}^{\mathrm{P}} \\ \begin{aligned} \tilde{\boldsymbol{\varepsilon}}^{\mathrm{P}} &= \dot{\boldsymbol{\varepsilon}}^{\mathrm{P}} + \mu \sqrt{J_2(\dot{\boldsymbol{\varepsilon}}^{\mathrm{P}})} \mathbf{I} \\ \tilde{\boldsymbol{\varepsilon}}^{\mathrm{P}} &\in \mathcal{N}_e(\boldsymbol{\sigma}) \end{aligned} \end{array} \right. \qquad (27) $$

The formulation of the constitutive equations remains general with respect to dimension. In Appendix A, the relation between $\eta, k, \phi$ in the case of plane strain assumption are detailed.
## 3 Numerical Method : Material Point Method

The Material Point Method (MPM) is a particle-based method dedicated to simulate the response of con- tinuum media to dynamical loadings developed in [2]. It relies on two spatial discretisations, a Eulerian grid and material points (also referred as particles in the literature). The material information such as stress, strain and velocity are stored in the materials points and the equilibrium equation are solved at the grid nodes. The solving scheme of the MPM can be separated into four main steps:

1. Interpolation of the quantities required to solve equilibrium equations from the material points to the grid nodes

2. Dynamic phase : solving of the equilibrium equations at the grid nodes and of the constitutive equation on the material points

3. Interpolation from the grid nodes to the material points

4. Advection of the material points with the newly interpolated velocity

5. Update of the density of the material points according to the continuity equation

The main scheme of the MPM described above can be interpreted in different ways. One convenient way is to consider a standard formulation of FEM with advected quadrature points (represented by the material points). The main advantage of this approach is to keep linear operations for advection but it induces a quadrature error in addition to the discretisation error already existing in the FEM [15, 16]. Note that, for classical explicit time integrators with MPM [39], the equilibrium equations and the constitutive laws in step 2 are uncoupled and a choice has to be made on the order of solving of these equations, leading to the well known Update Stress First (USF) or Update Stress Last (USL) schemes and their alternatives.

In order to obtain detailed expressions for step 1 to 5, it is necessary to develop the spatial discretisation of the weak form the equilibrium equation and the constitutive equation.

### Momentum Conservation Equation

The MPM formulation is classicaly obtained from the weak form of the momentum conservation equation:

$$ \forall w \in V, \int_{\Omega} \rho \frac{Dv}{Dt} \cdot w + \int_{\Omega} \boldsymbol{\sigma} : D(w) = \int_{\Omega} \rho b \cdot w + \int_{\partial\Omega} w \sigma_c \cdot n. \quad (28) $$

The last term in (28) corresponds to the forces applied at the boundaries of the domain and will be treated as a boundary condition and developed in section 4.2. The set V is the appropriate functional space. In equation (28), the operator $\frac{D\bullet}{Dt}$ and $D(\bullet)$ correspond respectively to the material derivative and the symmetric gradient. They are defined by

$$ \frac{D\bullet}{Dt} = \frac{\partial\bullet}{\partial t} + v\nabla(\bullet) \quad D(\bullet) = \nabla_s(\bullet) = \frac{1}{2}(\nabla(\bullet) + \nabla(\bullet)^T). \quad (29) $$
We now consider that the domain $\Omega$ is divided in space in $n_p$ subsets $(\Omega_p)_{p \in [1, n_p]}$, where $n_p$ is the number of material points. Note that these subdomains are dependent on time. The weak form described in (28) becomes

$$ \forall w \in V, \quad \sum_{p=1}^{n_p} \int_{\Omega_p} \rho \frac{Dv}{Dt} \cdot w + \sum_{p=1}^{n_p} \int_{\Omega_p} \boldsymbol{\sigma} : D(w) = \sum_{p=1}^{n_p} \int_{\Omega_p} \rho b \cdot w \quad (30) $$

As classically done in FEM, we suppose that the test functions and the velocity field belong to a finite dimensional space $V^h = \text{span}((N_i)_{i \in [1, n_g]}) \subset V$ where the functions $N_i$ are called interpolation functions centred at each grid node $(x_i)_{i \in [1, n_g]}$ and $n_g$ is the number of grid nodes.

$$ \forall w \in V^h \quad \sum_{p=1}^{n_p} \int_{\Omega_p} \rho \sum_{i=1}^{n_g} \dot{v}_i N_i \cdot \sum_{j=1}^{n_g} w_j N_j + \sum_{p=1}^{n_p} \int_{\Omega_p} \boldsymbol{\sigma} : \sum_{j=1}^{n_g} w_j \nabla N_j = \sum_{p=1}^{n_p} \int_{\Omega_p} \rho b \cdot \sum_{j=1}^{n_g} w_j N_j. \quad (31) $$

**Remark 2.** In (31), the components of the material derivative of the velocity become $\dot{v}_i$, the partial time derivative. This is only possible because the velocity at the nodes is interpolated from the velocity of advected particles.

As (31) is valid for any test function $w \in H$, we have for all grid node $j$

$$ \sum_{i=1}^{n_g} \dot{v}_i \sum_{p=1}^{n_p} \int_{\Omega_p} \rho N_i \cdot N_j + \sum_{p=1}^{n_p} \int_{\Omega_p} \boldsymbol{\sigma} : \nabla N_j = \sum_{p=1}^{n_p} \int_{\Omega_p} \rho b \cdot N_j. \quad (32) $$

In this work, as described in the original paper of [2], we choose to approximate the integrals in (32) by a zero order quadrature (the material quantities are considered as constant over each material domain $\Omega_p$). If $x_p$ denotes the position, depending on time, of the material point $p$, we make the approximation :

$$ \int_{\Omega_p} \phi \approx V_p \phi(x_p), \quad (33) $$

with $V_p = |\Omega_p|$, the volume of the material domain. Equation (32) becomes

$$ \sum_{p=1}^{n_p} \sum_{i=1}^{n_g} m_p N_i(x_p) \dot{v}_i N_j(x_p) = - \sum_{p=1}^{n_p} V_p \boldsymbol{\sigma}_p \nabla N_j(x_p) + \sum_{p=1}^{n_p} m_p b_p N_j(x_p), \quad (34) $$

where $m_p$ is the (constant) mass of each material domain $\Omega_p$. Hence (34) can be written as :

$$ \forall j \in [1, n_g] \quad \sum_{i=1}^{n_g} m_{ij} \dot{v}_j = f_j^{\text{int}} + f_j^{\text{ext}}, \quad (35) $$

where $M = (m_{ij})_{i,j \in [1, n_g]}$ is the grid mass matrix, $f_j^{\text{int}}$ and $f_j^{\text{ext}}$ represent respectively the internal and
external forces at the grid node j. They are defined as

$$f_j^{int} = - \sum_{p=1}^{n_p} V_p \sigma_p \nabla N_j(x_p), \qquad (36)$$

$$f_j^{ext} = \sum_{p=1}^{n_p} m_p b_p N_j(x_p), \qquad (37)$$

$$m_{ij} = \sum_{p=1}^{n_p} m_p N_i(x_p) N_j(x_p). \qquad (38)$$

Note that, in (35), the matrix $M := (m_{ij})$ is the consistent mass matrix on the grid nodes. In this work, as many others, we will use the so called lumped mass matrix (or diagonal mass matrix) where each diagonal component is defined by

$$m_i = \sum_{j=1}^{n_g} m_{ij} = \sum_{p=1}^{n_p} m_p N_i(x_p) \qquad (39)$$

$$\forall i \in [1, n_g], \quad m_i \dot{v}_i = f_i^{ext} + f_i^{int}. \qquad (40)$$

The use of the lumped mass matrix increases the energy dissipation during the different steps of the MPM [13]. However, the consistent mass matrix can be singular, leading to numerical issues. One solution is to use a weighting between the two matrices, which decreases energy losses but requires a better handling of the boundary conditions [40].

## Constitutive equations

The constitutive equations are solved at the material points which is similar to the option chosen in the FEM to solve these equations at quadrature points [1]. We discretise (6) over the material points

$$S : \sigma_p = D \left( \sum_{i=1}^{n_g} v_i N_i(x_p) \right) - \dot{\varepsilon}_p^P \qquad (41)$$

where $S$ is the inverse of the fourth order elasticity tensor and $D$ is the symmetric part of the gradient. Equation (41) can be rewritten as

$$S : \sigma_p = \frac{1}{2} \sum_{i=1}^{n_g} v_i^T \nabla N_i(x_p) + \nabla N_i(x_p)^T v_i - \dot{\varepsilon}_p^P. \qquad (42)$$

The plastic criterion and the plastic flow rule described in Section 2 are also expressed at all material points $p$ for both associated and non-associated plasticity

$$\dot{\varepsilon}_p^P \in \mathcal{N}_c(\sigma_p). \qquad (43)$$
## Description of one MPM time step

In this section, we present in details the MPM algorithm (steps 1 to 5) for a time step $[t^n, t^{n+1}]$. Let us denote $v_p^n, \sigma_p^n$ the quantities of the material points at $t^n, h := t^{n+1} - t^n$ the time step size. The interpolation scheme used in this work is the pure FLIP scheme. Even though the FLIP scheme is not stable with an explicit time integrator (where a combination between PIC and FLIP or APIC [12] is preferred), this scheme shows stability with an implicit time integrator [14].

As exposed in Section 3, the conservation and constitutive equations are discretized using the MPM following a scheme divided into five steps:

1. We compute the interpolated velocities $v_i^L$ and masses $m_i$ at each grid nodes from the material point velocities $v_p^n$ and masses $m_p$

$$m_i = \sum_p m_p N_i(x_p), \qquad (44)$$

$$v_i^L = \frac{1}{m_i} \sum_p v_p N_i(x_p^n) \qquad (45)$$

The interpolated velocities $v_i^L$ are also called the Lagrangian velocities, hence the notation.

2. The conservation and constitutive equations are solved by computing a solution $v_i^{n+1}, \sigma_p, \dot{\varepsilon}_p^P$ of the complementarity problem (40, 42, 43). Detailed explanations are given on the numerical solution procedure in Section 4. The nodal acceleration is then defined at each grid node $i$ by

$$\dot{v}_i^{n+1} = \frac{v_i^{n+1} - v_i^L}{h}. \qquad (46)$$

3. The new velocities at each material point are computed following the FLIP scheme, that interpolates the velocity increment at each material point $p$ with

$$v_p^{n+1} = v_p^n + \sum_i N_i(x_p^n)(v_i^{n+1} - v_i^L). \qquad (47)$$

4. The material points are advected by

$$x_p^{n+1} = x_p + h v_p^{n+1}. \qquad (48)$$

5. The material points density is updated with the new strain increments $\dot{\varepsilon}_p^{n+1}$ by

$$\rho_p^{n+1} = \frac{\rho_p^n}{1 + h \operatorname{tr}(\dot{\varepsilon}_p^{n+1})}. \qquad (49)$$
## 4 Solving of the incremental problem

This section is dedicated to the numerical procedure solving the incremental problem, in the form of a complementarity problem, that we have to solve in step 2 of the method. In order to write the incremental problem in a matrix form, we introduce the standard Voigt's notation for a two dimensional problem.

$$ \sigma = \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{pmatrix} \qquad \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} = \begin{pmatrix} \dot{\boldsymbol{\epsilon}}_{11}^{\mathrm{P}} \\ \dot{\boldsymbol{\epsilon}}_{22}^{\mathrm{P}} \\ 2\dot{\boldsymbol{\epsilon}}_{12}^{\mathrm{P}} \end{pmatrix} \qquad (50) $$

such that the inner product remains unchanged

$$ \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} = \boldsymbol{\sigma} \cdot \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} \qquad (51) $$

We define the operator $S$ as the inverse of the elasticity operator where

$$ \mathcal{D}^{-1} : \boldsymbol{\sigma} \xrightarrow{\text{Voigt's Notation}} S\boldsymbol{\sigma}. \qquad (52) $$

The operator $B$ is defined in a similar manner :

$$ \forall j \in \mathcal{G}, \quad \sum_{p} V_p \boldsymbol{\sigma}_p \nabla N_j(x_p) \xrightarrow{\text{Voigt's Notation}} (B^T \boldsymbol{\sigma})_j. \qquad (53) $$

The incremental problem, dedicated to calculate $v, \sigma, \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} \in \mathbb{R}^{2n_g} \times \mathbb{R}^{3n_p} \times \mathbb{R}^{3n_p}$ reads as

$$ \left\{ \begin{array}{l} M \left( \frac{v - v^L}{h} \right) + B^T \boldsymbol{\sigma} = F \\ S \left( \frac{\boldsymbol{\sigma} - \boldsymbol{\sigma}^n}{h} \right) = Bv - \dot{\boldsymbol{\epsilon}}^{\mathrm{P}} \\ \mathcal{K}^* \ni \tilde{\boldsymbol{\epsilon}}^{\mathrm{P}} \perp \left( \boldsymbol{\sigma} + \frac{k}{\eta} e \right) \in \mathcal{K} \end{array} \right. \qquad (54) $$

where $F \in \mathbb{R}^{2n_g}$ is the concatenation of the vector of external force at each material point, $e \in \mathbb{R}^{3n_p}$ is the concatenation of the Voigt's notation of $e$ defined in (22), the matrix $M$ is the lumped mass matrix on the grid nodes defined by (39). The set $\mathcal{K}$ is a second order cone of $\mathbb{R}^{3n_p}$ constructed by the Cartesian product of the second order cones at the different material points material points in Voigt's notation:

$$ \sigma + \frac{k}{\eta} e \in \mathcal{K} \iff \forall p \in [1, n_p], \quad \boldsymbol{\sigma}_p + \frac{k}{\eta} e \in \mathcal{K}. \qquad (55) $$

The modified plastic strain $\tilde{\boldsymbol{\epsilon}}^{\mathrm{P}}$ corresponds to the change of variable described in (25) in the case of non-associated plasticity. Appendix B details the change of variable one must handle in order to retrieve a second order cone complementarity problem using Voigt's notation. The complementarity problem (54) can be rewritten as

$$ \left\{ \begin{array}{l} Mv + hB^T \boldsymbol{\sigma} = hF + Mv^n \\ S\boldsymbol{\sigma} - hBv = -h\tilde{\boldsymbol{\epsilon}}^{\mathrm{P}} + S\boldsymbol{\sigma}^n \\ \mathcal{K}^* \ni \tilde{\boldsymbol{\epsilon}}^{\mathrm{P}} \perp \boldsymbol{\sigma} + \frac{k}{\eta} e \in \mathcal{K} \end{array} \right. \qquad (56) $$
By defining $q_v = hF + Mv^n + hB^T(k/\eta)e$, $q_\sigma = S\sigma^n + S(k/\eta)e$, and $z := -h\dot{\varepsilon}^P$ the incremental problem rewrites

$$ 
\left\{ 
\begin{array}{l}
Mv + hB^T\sigma = q_v \\
S\sigma - hBv = z + q_\sigma \\
\tilde{z} = \Psi(z) \\
\mathcal{K}^* \ni \tilde{z} \perp \sigma \in \mathcal{K}
\end{array}
\right. 
\qquad (57) 
$$

The transformation $\Psi$ corresponds to the change of variable due to the non-associated character of the plasticity. Note that the stress solving (57) corresponds to the translated stress taking into account the cohesion, and since $\mathcal{K}$ is a cone, $z/h \in \mathcal{K}$ is equivalent to $z \in \mathcal{K}$

**Remark 3.** In the incremental system (57), we use a fully implicit (backward Euler) scheme. It is also possible to use other integrators (θ-scheme, Newmark ...) for accuracy and energy considerations. The development for a θ-scheme would be similar to [1].

## 4.1 Solving of the complementarity conditions

The complementarity constraints can be expressed equivalently as finding the roots of a specific function. Following [41], we use the so called natural map to solve the plasticity conditions. It states that

$$ \mathcal{K}^* \ni \tilde{z} \perp \sigma \in \mathcal{K} \iff \forall \rho > 0, \sigma = P_{\mathcal{K}}(\sigma - \rho \tilde{z}), \qquad (58) $$

where the function $P_{\mathcal{K}}$ represents the projection onto the set $\mathcal{K}$ (which is well defined if $\mathcal{K}$ is convex). The parameter $\rho$ is a numerical parameter, which value is set to have relevant comparison between $\sigma$ and $\rho z$. The complementarity condition can then be solved by finding the roots of the natural map $F_{NM}$ defined as

$$ F_{NM}(\sigma, \tilde{z}) = \sigma - P_{\mathcal{K}}(\sigma - \rho \tilde{z}), \qquad (59) $$

The incremental problem presented in (57) is equivalent to finding the roots of the following function

$$ \mathcal{F} : \begin{pmatrix} v \\ \sigma \\ z \end{pmatrix} \mapsto \begin{pmatrix} Mv - hB^T - q_v \\ S\sigma - hBv - z - q_\sigma \\ F_{NM}(\sigma, \Psi(z)) \end{pmatrix} \qquad (60) $$

## Reduced system

The incremental problem presented in (57) can be reduced in the variable $v$ by multiplying the first row of the system (57) by $M^{-1}$ and replacing the expression of $v$ in the second row. This leads to the definition of the Schur complement (or sometimes called the Delassus operator) $W = S + h^2BM^{-1}B^T$, so that (57) becomes equivalent to

$$ 
\left\{ 
\begin{array}{l}
z = W\sigma + \hat{q} \\
\mathcal{K}^* \ni \Psi(z) \perp \sigma \in \mathcal{K},
\end{array}
\right. 
\qquad (61) 
$$

where $\hat{q} = -q_\sigma - hBM^{-1}q_v$. The system described in (61) is the root of the reduced function defined as

$$ \mathcal{F}^R : \begin{pmatrix} \sigma \\ z \end{pmatrix} \mapsto \begin{pmatrix} W\sigma - z + \hat{q} \\ F_{NM}(\sigma, \Psi(z)) \end{pmatrix}. \qquad (62) $$
For associated plasticity, the system described in eq. (61) becomes equivalent to the following minimisation problem:

$$ \min_{\sigma \in \mathcal{K}} \frac{1}{2} \sigma^T W \sigma + \hat{q}. \qquad (63) $$

The fact that $W$ is symmetric definite positive ensures good properties of the incremental problem such as existence and uniqueness as the root of the function described in (62) is the optimality condition of (63).

## 4.2 Dirichlet Boundary Conditions

### Dirichlet BC applied on grid nodes

The first way to impose Dirichlet BC is to impose them on grid nodes. We then choose a subset $G^{BC} \subset [1, n_g]$ of size $n^{BC}$ such that

$$ \forall i \in G^{BC}, \quad v_i = \bar{v}_i, \qquad (64) $$

where $\bar{v}_i$ is the desired velocity at node $i$. The condition in (64) can be written as a condition on the vector $v$

$$ Hv = \bar{v}, \qquad (65) $$

where $H$ is a selector operator, filled with 0 and 1 selecting only the nodes in $G^{BC}$, where the boundary conditions are applied.

### Dirichlet BC on material points

It is also possible to impose the Dirichlet BC directly on the material points. We select a set of $n^{BC}$ chosen material points $MP^{BC} \subset [1, n_p]$ such that we want to ensure that, at the end of the time step, the velocities at these nodes correspond to the prescribed values. It means that

$$ \forall p \in MP^{BC} \quad v_p^{n+1} = \bar{v}_p, \qquad (66) $$

where $\bar{v}_p$ is the desired velocity at the material point $p$. Once again, this can be written as a condition on $v$. If we note $N$ the matrix of interpolation such that, at step 4 of the MPM algorithm, the FLIP scheme writes as

$$ v^{n+1} = v^n + N(v - v^L). \qquad (67) $$

We can construct the selection matrix $A$ in the same manner as (65) such that (66) becomes

$$ A v^{n+1} = \bar{v}. \qquad (68) $$

Injecting (67) into (68) gives

$$ ANv = \bar{v} + A(v^n - Nv^L). \qquad (69) $$

By posing $H := AN$, we retrieve a similar form as (65). In all cases, Dirichlet boundary conditions are treated as equality constraints ensured by an equation of the form

$$ Hv = q_{BC}. \qquad (70) $$

After integrating the Dirichlet boundary conditions in (70) as equality constraints, the system to be solved is
$$ egin{cases} Mv + hB^T\sigma = q_v + H^T\lambda \\ S\sigma - hBv = z + q_\sigma \\ Hv = q_{BC} \\ \tilde{z} = \Psi(z) \\ \mathcal{K}^* \ni \tilde{z} \perp \sigma \in \mathcal{K} \end{cases} \tag{71} $$

where the vector $\lambda$ is the Lagrange multiplier associated to the constraint $Hv = q_{BC}$.

The method developed here is very close to the method presented in the recent work of [42]. The differences lie in the fact that they consider a variational formulation of the boundary conditions by integrating them over the boundary domain. In our case, we impose a strong formulation of the boundary conditions, which avoids having to introduce boundary particles. However, our formulation comes directly from the discretisation and does not derive from a variational principle. In order to retrieve a system similar to system (56), we can reduce the problem by expressing the variable $\lambda$ depending on $v$ and $\sigma$, we then have

$$ egin{cases} Mv + h\tilde{B}^T\sigma = \tilde{q}_v \\ S\sigma - hBv = z + q_\sigma \\ \tilde{z} = \Psi(z) \\ \mathcal{K}^* \ni \tilde{z} \perp \sigma \in \mathcal{K} \end{cases} \tag{72} $$

This operation is possible by introducing the matrix $L = HM^{-1}H^T$ and replacing

$$ B \rightarrow \tilde{B} = (I - M^{-1}H^T L^{-T} H)B \\ q_v \rightarrow \tilde{q}_v = q_v - H^T L^{-1} H M^{-1} q_v + H^T L^{-1} q_{BC}. $$

As for (56), solving (72) is equivalent to finding the zeros of the following function

$$ \tilde{F}: \begin{pmatrix} v \\ \sigma \\ z \end{pmatrix} \mapsto \begin{pmatrix} Mv - \tilde{B}^T\sigma - \tilde{q}_v \\ S\sigma - Bv - z - q_\sigma \\ F_{NM}(\sigma, \Psi(z)) \end{pmatrix}. \tag{73} $$

Note that, at first sight, we see a loss of symmetry in the linear part of system (72) induced by the differences between $\tilde{B}$ and $B$. However if we reduce again the system by replacing the velocity, we have a new system, similar to the reduced system (61)

$$ egin{cases} z = \tilde{W}\sigma + \hat{q} \\ \mathcal{K}^* \ni \Psi(z) \perp \sigma \in \mathcal{K}, \end{cases} \tag{74} $$

where $\tilde{W} = S + h^2 BM^{-1} \tilde{B}^T$. With the definition of $\tilde{B}$, we have

$$ \tilde{W} = S + h^2 BM^{-1}B^T - h^2 BM^{-1}H^T L^{-1} H M^{-1} B^T, \tag{75} $$

which is symmetric since $L$ and $M$ are symmetric.
## 4.3 Non smooth Newton's method

In order to compute the zeros of the function described in (62), (60) and (73), we use the Newton's method. As the function $F_{NM}$ is not differentiable everywhere (as it includes the projection on the stress admissible set), we have to use the so-called non-smooth or semi-smooth Newton's method, see [41] for more details. The main difference with the classic Newton's method lies in the fact that we have to consider a generalised jacobian to compute the new iterate:

$$x_{k+1} = x_k - V^{-1}F_{NM}(x_k), \quad V \in \partial F_{NM}(x_k). \qquad (76)$$

The generalised jacobian $\partial F_{NM}(x_k)$ is described in appendix C. For associated plasticity it can be shown [43] that the method converges locally for any choice of the parameter $\rho$. No theoretical results exist yet for non-associated plasticity but in practice, the method converges with the same parameters than for associated plasticity. As done classically, the algorithm stops when $\|F_{NM}(x_k)\|$ is smaller than a given tolerance.
# 5 Results and discussion

All the simulations presented have been done with our python implementation of the material point method. When specified, the incremental problem are solved with the non-smooth Gauss-Seidel solver included in the Siconos package. We refer to [32, 30] for a presentation of this solver.

## 5.1 Footing

We first illustrate our method on simulations of footing problems. This configuration has been extensively studied in soil mechanics and there exists in the literature plenty of analytical solutions based on limit analysis, in particular [44, 45, 46]. During a footing experiment or simulation, a soil mass is subjected to localized loading at its surface to evaluate its bearing capacity (Figure 3). The loading conditions consist of either applying a load and measuring the resulting displacement of the loading plate, or prescribing a displacement $u_d$ to the plate and measuring the corresponding reaction forces. In the simulations, the second option has been chosen. Consequently, Dirichlet boundary conditions on either grid nodes or material points (section 4.2) have been applied with a prescribed velocity set at 0.02 m s⁻¹.

The ultimate bearing capacity $Q$, defined as the limit load reached in footing tests, is used for comparison with literature values. In the literature, $Q$ is classically calculated by summing three distinct contributions, $Q_c$, $Q_\gamma$ and $Q_q$. The first one considers the effect of cohesion, the second one the effect of soil weight and the third one the effect of the overload. In this example, we consider no overload and a weightless soil, the ultimate bearing capacity can thus be calculated accounting only for the effect of cohesion ($Q = Q_c$). It can be computed by limit load analysis

$$Q_c = \frac{c}{\tan \phi} \left( e^{\pi \tan(\phi)} \tan \left( \frac{\pi}{4} + \frac{\phi}{2} \right)^2 - 1 \right). \qquad (77)$$

The calculation of $Q_c$ from (77) does not account for soil inertia and elasticity and only considers associated plasticity. It gives a lower bound for simulations with associated materials. The mechanical and numerical parameters used in the footing simulations are listed in Table 1. In order to minimise the error appearing when solving the incremental problem, the tolerance of the non-smooth Newton's method has been set to $10^{-10}$.
<figure><img src="image_4.png" /><figcaption>Figure 3: Schematic of the footing simulation</figcaption></figure>

We first performed simulations for both associated and non-associated material with boundary conditions applied on chosen grid nodes, as described in Figure 7. The range of material parameters explored in the simulations (Table 1) has been set according to the literature [38, 29]).

In Figure 5, we observe, in accordance with the literature, that the plastic deformations are localized around the loading zone and delimited first by a straight line located under the footing and starting from the edge of the footing making an angle of $\frac{\pi}{4} - \frac{\phi}{2}$ with the vertical axis. From the lower extremity of this line, the plastic deformations zone is located along an arc first directed towards the bottom of the soil sample and then towards the surface until this arc crosses a line starting from the footing edge, pointing towards the exterior of the footing and making an angle of $\frac{\pi}{4} + \frac{\phi}{2}$ with the vertical axis. Finally, the plastic deformation zone is located along a straight line tangent to the previously described, thus directed toward soil surface, see for example [38, 44] for more details.

A typical load response of the plate is represented in Figure 4. As expected, the load progressively increases until it reaches a peak value and then slightly decreases before stabilizing at a near-plateau. The limit load, corresponding to the near-plateau, is around 8.3 kN, which is in agreement with the analytical value obtained from (77).
<table><caption>Table 1: Ranges of parameters used for footing simulations</caption><tbody><tr><td rowspan="2">Geometry</td><td>Length L<sub>1-2</sub> (m)</td><td>2</td></tr><tr><td>Width b (m)</td><td>0.2</td></tr><tr><td rowspan="7">Material</td><td>Density ρ (kg·m<sup>-3</sup>)</td><td>2·10<sup>3</sup></td></tr><tr><td>Young's modulus E (MPa)</td><td>1·10<sup>8</sup></td></tr><tr><td>Poisson's ratio ν</td><td>0.3</td></tr><tr><td>Friction angle φ (°)</td><td>[20, 35]</td></tr><tr><td>Dilatancy angle θ (°)</td><td>[0, 35]</td></tr><tr><td>Cohesion c (kPa)</td><td>[2, 5]</td></tr><tr><td>boundary velocity u<sub>d</sub> (m.s<sup>-1</sup>)</td><td>2·10<sup>-2</sup></td></tr><tr><td rowspan="6">Numerical</td><td>Number of particles N<sub>p</sub></td><td>6400</td></tr><tr><td>Cell size dx (m)</td><td>[0.05, 1]</td></tr><tr><td>Time step h (s)</td><td>1·2<sup>-2</sup></td></tr><tr><td>Simulation time T (s)</td><td>[1, 2]</td></tr><tr><td>Interpolation function</td><td>Cubic</td></tr><tr><td>Interpolation scheme</td><td>FLIP</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: Force against time in a footing simulation for a cohesive soil with comparison of output of simulation with an analytical solution for c = 2kPa and φ = 25°. In blue : the associated case (θ = φ) and in orange the strongly non-associated case (θ = 0°) Boundary conditions applied on grid nodes.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: a) Plastic strain rate at the material points at the end of a footing simulation with $c = 2\,\text{kPa}$ and $\theta = \phi = 25^\circ$. b) Plastic strain rate at the material points at the end of a footing simulation with $c = 2\,\text{kPa}$, $\phi = 25^\circ$ and $\theta = 0^\circ$.</figcaption></figure>

For a fixed value of cohesion ($c = 2\,\text{kPa}$), the limit load increases for increasing friction angles (from $\phi = 20^\circ$ to $\phi = 35^\circ$). This result is in accordance with the analytical solutions although we observe a small relative difference ($\le 5\%$) between the maximum load in the simulation and the solution from limit analysis. As expected, the limit loads from the simulations are always greater than the solutions from limit analysis.
<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison of simulation results and analytical solutions for the lineic bearing capacity in kN/m with a cohesion $c = 2\text{kPa}$ and for different friction angles.</figcaption></figure>

The effect of dilatancy is highlighted by performing simulations with different dilatancy angles, ranging from the associated case ($\phi = \theta$) and a strongly non-associated case ($\theta = 0^\circ$). For a strongly non-associated material ($\theta = 0^\circ$ - Figure 5.b), we observe that the plastic deformation zone is smaller than for an associated material ($\phi = \theta$ - Figure 5.a), with a shape similar to the one obtained in [38]. In addition, the increase in volume is significantly smaller at the top of the plasticity zone, in particular near the load plate, because of significantly smaller soil dilatancy. The load response for a strongly non-associated material (Figure 6) is characterized by a more progressive load increase compared to an associated material, followed by a more clearly defined plateau, slightly smaller than the analytical value obtained from (77).
Challenges related to the boundary conditions modelling

<figure><img src="image_8.png" /><figcaption>Figure 7: Schematic of the boundary conditions applied in footing simulation with MPM. On top left are represented the BC on material points and on bottom left the BC on the grid nodes.</figcaption></figure>
As explained in section 4.2, the boundary conditions can be either imposed on grid nodes or directly imposed on the material points. For footing simulations, the differences between the two types of Dirichlet boundary conditions are illustrated in Figure 7. For large values of cohesion $c$ or frictional angle $\phi$, the footing displacement needed to reach the ultimate bearing capacity increases. When boundary conditions are imposed on grid nodes, this results in re-meshing issues that produce sudden decreases in the apparent bearing force, as illustrated in Figure 8 at time $t = 0.5$ s.

<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison between the simulations and the analytical solution for bearing capacit, with $c = 4$ kPa, $E = 10^2$ MPa, $h = 0.05$ s, $dx = 0.1$ m. boundary condition applied on grid nodes</figcaption></figure>

In the case of boundary conditions imposed on material points (Figure 9), the load increases more rapidly than when boundary conditions are imposed on grid nodes (Figure 8). In addition, the peak values are significantly more pronounced and followed by either decrease in the load or oscillations before reaching a plateau. These differences can be attributed to dynamical effects related to soil inertia that are more pronounced because significant velocity is applied directly on the material points at the beginning of the simulation.

The boundary conditions applied to material points can handle larger deformations than boundary conditions applied to grid nodes. This is well illustrated by the constant plateau on Figure 9 compared to the sudden decrease displayed in Figure 8. However, due to the large support of the shape functions, material points close to the loading zone are affected by the imposed velocity, as illustrated in Figure 10. The material points where boundary conditions are applied are framed in red and the material points in the grid cell adjacent to the loading are framed in green. The latter have a vertical velocity similar to the imposed vertical velocity $\bar{v}$. This has for consequence that more material points are "driven" that expected. In Figure 9, we compare the Prandtl's solution with the length $B$ and $B + dx$ and we observe a decent agreement by taking into account $B + dx$ in the analytical solution. Note that the difference between the black and the blue curves in Figure 9 is heavily marked due to the fact that the simulation is performed
with a coarse discretisation ($\frac{B}{dx} = 2$). The decrease in the size of $dx$ significantly decreases this difference (Figure 9), which illustrates that this limitation can be corrected by increasing mesh resolution.

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

Figure 9: Comparison between the simulations and the analytical solutions considering as loading length $B$ (in black) or $B + dx$ (in blue and red for different value of $dx$), with $c = 4\text{kPa}$, $E = 10^2\text{MPa}$, $h = 0.2\text{s}$, $dx = 10\text{ cm}$ and $dx = 5\text{ cm}$. Boundary conditions applied on material points.
<figure><img src="image_11.png" /><figcaption>Figure 10: Vertical velocities of the material points at the end of footing simulation with boundary conditions applied on material points. The loading zone is delimited by *B*. In black is displayed the background grid.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Force against displacement in footing simulation with boundary conditions on the material points. Comparisons of the outputs of simulation in associated and non-associated cases for a $c = 4$ kPa, $\phi = 25^\circ$, a time step $h = 0.02$ s and a cell size of $dx = 0.1$ m. In blue $\theta = \phi$, in orange $\theta = \frac{\phi}{2}$ and in green $\theta = 0^\circ$.</figcaption></figure>

In order to illustrate the capacities of the method proposed combined with Dirichlet boundary conditions applied on material points, we performed simulations for different values of the dilatancy angle (Figure 11). As expected from the literature, we observe that, as the value of the dilatancy angle decreases, the bearing capacity (represented by the value of the plateau of the curve in Figure 11) decreases. It is also worth noting that, even for a nul dilatancy angle and large time steps as regards to the discretisation ($h = \bar{v}/dx$), the solver manages to find the solution of the incremental problem in few iterations.

## Comparison with Finite Element Method simulations

During the footing simulations, only small strains of the soil occur. Consequently, the simulations obtained with our MPM approach can be compared with more classical FEM simulations. An incremental problem having a similar structure as the incremental problem presented in (57) can be written using the FEM [1]. This problem can be solved using the same solver as for our MPM simulations. One can note that, for the FEM simulations, we used a third order triangular element with an exact integration in order to prevent any locking phenomena, as suggested in [38].

Figure 12 shows the distribution of the plastic strain rate, averaged on each element at the end of the footing simulation. We observe a good agreement with the expected shear band from limit analysis and with the outputs from MPM simulations displayed in Figure 5. The load-displacement curve obtained is displayed in Figure 13. We observe a perfect agreement between the converged value of the force and the Prandlt's solution from limit analysis. As for MPM simulations with boundary conditions applied on material points, the force oscillations before reaching the plateau can be attributed to dynamical effects. Footing simulations were performed with different values of dilatancy angle, from associated case to strongly non-associated case.
Figure 13 represents the force against the displacement for three values of $\theta$. As expected from [38], we observe that the limit value of the force decreases as the dilatancy angle decreases.

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

Figure 12: Plastic strain rate averaged on the elements at the end of footing simulation with order 3 FEM for $c = 4\text{kPa}$, $\phi = 25^\circ$. The lines in black represents the expected borders of the plastic zone.
<figure><img src="image_14.png" /><figcaption>Figure 13: Force against displacement in footing simulations with third order FEM. Comparison of outputs of simulation for several values of the dilatancy angle, for $c = 4$ kPa with the analytical Prandtl's solution. In blue $\theta = \phi$, in orange $\theta = \phi/2$ and in green $\theta = 0^{\circ}$.</figcaption></figure>

## Solver performance

The solver's performance can be assessed on the basis of its accuracy, robustness and computation time. The tolerance of the solver is set at $10^{-10}$, guaranteeing a high precision of the computed solution. The first iterate, computed either by fixed point or with non-smooth Gauss-Seidel algorithm is used to find an approximation of the solution and use it as first iterate of the semi-smooth Newton's method, as the convergence of such method with a second order rate is only local. In the scope of the paper, the solver is evaluated based on the precision of the solution and the number of iterations of the semi-smooth Newton's method.

The robustness of the solver is illustrated by the fact that by keeping $h \sim \bar{v}/dx$, the solver manage to find the solution in few iterations independently of the value of $dx$. This condition comes from the fact that the material points must move no more than one cell size during one time step and that the maximal velocity of the material points is around $\bar{v}$ in footing simulations. The logarithm of the error versus the number of iteration is displayed in Figure 14 for associated and non-associated simulations. We observe as expected from [41], a local convergence with a second order rate. The yield criteria and the plastic flow rule are satisfied at each material points within the solver precision. Moreover, at the end of the Newton's method, the residual scalar product between the two complementary variables in (71) is also of the order of the solver precision.
<figure><img src="image_15.png" /><figcaption>Figure 14: Illustration of the solver performance. Top left: $\sqrt{J_2}(\sigma)$ as a function of $I_1(\sigma)$, each line corresponds to a material point through time. Top right: $\sqrt{J_2}(\dot{\epsilon}^P)$ as a function of $I_1(\dot{\epsilon}^P)$, each line corresponds to a material point through time. Bottom left: evolution at different time steps of the scalar product $\tilde{\epsilon}^P \cdot (\sigma + \frac{k}{\eta}e)$, each line corresponds to a material point. Bottom right: logarithm of the error of the semi-smooth Newton's method at each iteration, each line corresponds to a time step.</figcaption></figure>

## 5.2 Slope stability

To illustrate the method and highlight the benefits of the hybrid nature of MPM, we apply the method on a well known problem from geomechanics: the slope stability analysis. We considered for these simulations the geometry described in Figure 15. The simulations are performed by applying the gravity at the first
time step. The mechanical parameters considered are an internal friction angle $\phi = 25^\circ$, a Young's modulus $E = 100$ MPa, a density $\rho = 2 \times 10^3$ kg m$^{-3}$, a cohesion $c = 5$ kPa, a mesh cell size of $dx = 50$ cm and a time step of $h = 0.25$ s. The total simulation time is $T = 5$ s.

<figure><img src="image_16.png" /><figcaption>Figure 15: Schematic and dimensions for slope stability simulations.</figcaption></figure>

We first performed the simulations considering that the slope material is associated. Figure 16 displays the cumulative plastic strains and the positions of the material points for 4 different steps of the simulation. We observe, as expected from the literature, the apparition of a shear band which corresponds to a slip surface (or slip line in plane strain), clearly identified in Figure 16 c) and d). The cumulative plastic strains reach 200% at the end of the simulation, illustrating the fact that the MPM can process large deformations. The material in the upper part of the shear band does not undergo plastic transformations and behaves similarly to a rigid body. The rest of the material beneath the shear bands remains in the elastic domain. Due to material, we observe that the width of the shear band increases with time, leading to an unrealistic behaviour.
<figure><img src="image_17.png" /><figcaption>Figure 16: Cumulative plastic strains for an associated material ($\theta = \phi$) for different simulation times: a) $t = 0$ s, b) $t = 1.7$ s, c) $t = 3.3$ s, d) $t = 5$ s.</figcaption></figure>

The same simulations but with for a strongly non-associated material ($\theta = 0^\circ$) are performed. Figure 17 displays the cumulative plastic strains and the position of the material points for 4 different steps of the simulation. We observe again a shear band corresponding to a slip line delimiting a static part and a moving part with almost uniform velocity. The cumulative plastic strains reach more than 300% during the simulation. In the non-associated case, the shear band goes more deeply into the bulk of the material, leading to a higher volume of material being moved than in the associated case. This difference is highlighted by comparing the length of the static part at the top of the material in Figures 16 d) and 17 d). As expected the non-associated character of the material leads to no increase of volume within the shear band.
<figure><img src="image_18.png" /><figcaption>Figure 17: Cumulative plastic strains for a strongly non-associated material ($\theta = 0^{\circ}$) and for different simulation times: a) $t = 0$ s, b) $t = 1.7$ s, c) $t = 3.3$ s, d) $t = 5$ s.</figcaption></figure>

During the slope stability simulations, there is a transfer from elastic energy to kinetic energy. In a MPM simulation at time $t^k$, the kinetic energy $E_{kin}^k$, the elastic energy $E_{el}^k$ and the potential energy $E_{pot}^k$ can be computed as

$$E_{kin}^k = \frac{1}{2} \sum_{p=1}^{n_p} m_p |v_p^k|^2 \quad E_{el} = \frac{1}{2} \sum_{p=1}^{n_p} \sigma_p^k : S : \sigma_p^k \quad E_{pot} = E_0 + \sum_{p=1}^{n_p} m_p g h_p^k, \qquad (78)$$

where $h_p^k$ is the height of the material point $p$ at time $t^k$ and $E_0$ is a reference energy set such as $E_{pot} = 0$ at the end of the simulation. The mechanical energy $E_{mec}^k$ is defined as the sum of the potential energy, the elastic energy and the potential energy. The plastic energy dissipation $E_{dis}^k$ is computed as

$$E_{dis}^k = \sum_{l=0}^{k} h \sum_{p=1}^{n_p} \sigma_p : \dot{\varepsilon}_p^P, \qquad (79)$$
where *h* is the time step. In this framework, regardless the value of the dilatancy angle, the dissipation must remain non negative [29]. We can define the total energy as the sum of the mechanical energy and the dissipation $E_{tot} = E_{mec} + E_{dis}$

<figure><img src="image_19.png" /><figcaption>Figure 18: Evolution of the energy balance through time for a simulation of slope stability for an associated material. In blue the kinetic energy $E_{kin}$, in orange the elastic energy $E_{el}$, in green the potential energy $E_{pot}$, in red the mechanical energy $E_{mec}$, in purple the plastic energy dissipation $E_{dis}$ and in brown the total energy $E_{tot}$</figcaption></figure>

Figure 18 shows the evolution of the different energies during a simulation involving an associated material. At the first time step, we observe sudden increase of the elastic energy and a decrease of the potential energy, this is due to the fact that the material point are initiated to nil stress and gravity is applied at the first time step. During slope failure, the slope height decreases leading to a decrease of the potential energy. During shear band formation, the material is progressively moving along the slip line and thus the dissipation due to plastic strain increases. The evolution of the kinetic energy indicates that the slipping part of the slope has a maximum velocity at time $t = 1.5$ s and that, at the end of the simulation, a stable state is reached as kinetic energy is equal to zero.
<figure><img src="image_20.png" /><figcaption>Figure 19: Evolution of the energy balance through time for a simulation of slope stability for a non-associated material. In blue the kinetic energy $E_{kin}$, in orange the elastic energy $E_{el}$, in green the potential energy $E_{pot}$, in red the mechanical energy $E_{mec}$, in purple the plastic dissipation $E_{dis}$ and in brown the total energy $E_{tot}$</figcaption></figure>

The energy balance is also displayed for a non-associated material in Figure 19 where the evolution of the different energies evolve globally in a similar manner as in the simulation for an associated material. Two main differences are to be noticed. First, the variation of the potential energy is five time higher in the non-associated case. This is due to the fact that in the associated case, as there is a lot of dilatancy, the mean height used to compute the reference energy $E_0$ at the end of the simulation is significantly smaller. Second, the dissipation in the non-associated case is close to ten times greater than the dissipation in the associated case. This is due to the fact that in compression and shear, the normalized scalar product between $\dot{\epsilon}^P$ and $\sigma$ is larger in the non-associated case, as illustrated in Figure 2. What is the most important to notice is that in both non-associated and associated case the total energy $E_{tot}$ is non increasing during the simulation illustrating the fact that the presented method is dissipative regardless of the value of the dilatancy angle $\theta$, even with large time step ($h = 0.25$ s). The decrease of the total energy is less than 5% of its initial value. The main source of dissipation are the fact that we use fully implicit scheme instead of a $\theta$-scheme [1] and that we consider the lumped mass matrix instead of the consistent mass matrix [47, 14]. Having theoretical results in order to quantify the energy loss in the simulation will be the topic of further studies. Overall, the presented method is particularly well-suited for this type of geomechanical application. The hybrid nature of the MPM effectively handles large deformations. The solver converges even with large time steps with a relatively fine discretisation, while maintaining a good energy balance and producing a physically realistic solution.
# 6 Conclusion

This paper presents a new solver inspired from non-smooth dynamics and second order cone programming to handle non-associated plasticity. The plastic criteria and the flow rule are written as a cone complementarity problem. This solver is used to solve the incremental problem that comes from either MPM or FEM. Dirichlet boundary conditions, applied either directly on the material point or the grid nodes are considered as a equality constraints on the velocity field and are solved simultaneously with the plastic laws. The solver is illustrated on the well known footing problem, for which it has been compared with analytical solution coming from limit analysis and shows good agreement with the analytical solutions. The solver coupled with MPM has been shown to be precise and robust even with large time steps. The solver can also used with FEM to show similar results. The method is further validated on a standard geomechanics problem : the slope stability analysis. In addition, the presented method is proving to be dissipative with little loss of energy. The future of this work includes an extension of the solver to solve simultaneously the Signorini contact, the coulomb friction, the non-associated plasticity and the boundary conditions in one monolithic solver. Future works include also a theoretical and numerical study of the energy conservation of the method.
References

[1] V. Acary, F. Bourrier, and B. Viano, "Variational approach for nonsmooth elasto-plastic dynamics with contact and impacts," *Comput. Methods Appl. Mech. Eng.*, vol. 414, p. 116156, Sep. 2023.

[2] D. Sulsky, Z. Chen, and H. L. Schreyer, "A particle method for history-dependent materials," *Comput. Methods Appl. Mech. Eng.*, vol. 118, no. 1, pp. 179-196, Sep. 1994.

[3] M. W. Evans and F. Harlow, "The Particle-In-Cell Method for Hydrodynamic Calculations," *Scientific Laboratory Reports, Los Alamos*, Jan. 1957.

[4] J. U. Brackbill, D. B. Kothe, and H. M. Ruppel, "Flip: A low-dissipation, particle-in-cell method for fluid flow," *Comput. Phys. Commun.*, vol. 48, no. 1, pp. 25-38, Jan. 1988.

[5] K. Abe and K. Konagai, "Numerical simulation for runout process of debris flow using depth-averaged material point method," *Soils and Foundations*, vol. 56, no. 5, pp. 869-888, Oct. 2016.

[6] J. Gaume, T. Gast, J. Teran, A. van Herwijnen, and C. Jiang, "Dynamic anticrack propagation in snow," *Nat. Commun.*, vol. 9, no. 3047, pp. 1-10, Aug. 2018.

[7] B. Trottet, R. Simenhois, G. Bobillier, B. Bergfeld, A. van Herwijnen, C. Jiang, and J. Gaume, "Transition from sub-Rayleigh anticrack to supershear crack propagation in snow avalanches," *Nat. Phys.*, vol. 18, pp. 1094-1098, Sep. 2022.

[8] L. Guillet, L. Blatny, B. Trottet, D. Steffen, and J. Gaume, "A Depth-Averaged Material Point Method for Shallow Landslides: Applications to Snow Slab Avalanche Release," *J. Geophys. Res. Earth Surf.*, vol. 128, no. 8, p. e2023JF007092, Aug. 2023.

[9] H. Rousseau, J. Gaume, L. Blatny, and M. P. Lüthi, "Transition Between Mechanical and Geometric Controls in Glacier Crevassing Processes," *Geophys. Res. Lett.*, vol. 51, no. 9, p. e2024GL108206, May 2024.

[10] A. Huth, R. Duddu, and B. Smith, "A Generalized Interpolation Material Point Method for Shallow Ice Shelves. 2: Anisotropic Nonlocal Damage Mechanics and Rift Propagation," *J. Adv. Model. Earth Syst.*, vol. 13, Aug. 2021.

[11] G. Daviet and F. Bertails-Descoubes, "A semi-implicit material point method for the continuum simulation of granular materials," *ACM Trans. Graphics*, vol. 35, no. 4, pp. 1-13, Jul. 2016.

[12] C. Jiang, C. Schroeder, and J. Teran, "An angular momentum conserving Affine-Particle-In-Cell method," *arXiv*, Mar. 2016.

[13] E. Love and D. L. Sulsky, "An unconditionally stable, energy-momentum consistent implementation of the material-point method," *Comput. Methods Appl. Mech. Eng.*, vol. 195, pp. 3903-3925, Jul. 2006.

[14] G. Pretti, W. M. Coombs, C. E. Augarde, B. Sims, M. Marchena Puigvert, and J. A. R. Gutiérrez, "A conservation law consistent updated Lagrangian material point method for dynamic analysis," *J. Comput. Phys.*, vol. 485, p. 112075, Jul. 2023.
[15] A. Baumgarten and K. Kamrin, "Analysis and mitigation of spatial integration errors for the material point method," *Int. J. Numer. Methods Eng.*, Jan. 2023.

[16] M. Berzins, "Computational error estimation for the Material Point Method," *Comp. Part. Mech.*, vol. 10, no. 4, pp. 865-886, Aug. 2023.

[17] S. Duverger, "A multi-scale, MPMDEM, numerical modelling approach for geotechnical structures under severe loading," Ph.D. dissertation, Aix-Marseille Université (AMU), Feb. 2023. [Online]. Available: https://theses.hal.science/tel-04101270v1

[18] D. Sulsky and A. Kaul, "Implicit dynamics in the material-point method," *Computer Methods in Applied Mechanics and Engineering - COMPUT METHOD APPL MECH ENG*, vol. 193, pp. 1137-1170, Mar. 2004.

[19] R. E. Bird, G. Pretti, W. M. Coombs, C. E. Augarde, Y. U. Sharif, M. J. Brown, G. Carter, C. Macdonald, and K. Johnson, "An implicit material point-to-rigid body contact approach for large deformation soil-structure interaction," *Comput. Geotech.*, vol. 174, p. 106646, Oct. 2024.

[20] Z. Zong, C. Jiang, and X. Han, "A Convex Formulation of Frictional Contact for the Material Point Method and Rigid Bodies," in *International Conference on Intelligent Robots and Systems (IROS)*. IEEE, 2024, pp. 14-18.

[21] J. C. Simo, J. G. Kennedy, and R. L. Taylor, "Complementary mixed finite element formulations for elastoplasticity," *Computer Methods in Applied Mechanics and Engineering*, vol. 74, no. 2, pp. 177-206, 1989.

[22] P. W. Christensen, "A nonsmooth Newton method for elastoplastic problems," *Comput. Methods Appl. Mech. Eng.*, vol. 191, no. 11, pp. 1189-1219, Jan. 2002.

[23] K. Krabbenhoft, M. R. Karim, A. V. Lyamin, and S. W. Sloan, "Associated computational plasticity schemes for nonassociated frictional materials," *Int. J. Numer. Methods Eng.*, vol. 90, no. 9, pp. 1089-1117, Apr. 2012.

[24] K. Krabbenhøft, A. V. Lyamin, and S. W. Sloan, "Formulation and solution of some plasticity problems as conic programs," *Int. J. Solids Struct.*, vol. 44, no. 5, pp. 1533-1549, Mar. 2007.

[25] X.-W. Zhou, Y.-F. Jin, K.-Y. He, Z.-Y. Yin, and F.-T. Liu, "A novel implicit FEM-MPM coupling framework using convex cone programming for elastoplastic problems," *Comput. Methods Appl. Mech. Eng.*, vol. 429, p. 117153, Sep. 2024.

[26] R. Hill, "A VARIATIONAL PRINCIPLE OF MAXIMUM PLASTIC WORK IN CLASSICAL PLASTICITY," *Q. J. Mech. Appl. Math.*, vol. 1, no. 1, pp. 18-28, Jan. 1948.

[27] B. Halphen and Q. S. Nguyen, "Sur les matériaux standard généralisés," *Journal de Mécanique*, vol. 14, no. 1, pp. 39-63, 1975. [Online]. Available: https://hal.science/hal-03600755

[28] J. J. Moreau, "On Unilateral Constraints, Friction and Plasticity," in *New Variational Techniques in Mathematical Physics*. Berlin, Germany: Springer, 2011, pp. 171-322.
[29] A. Berga and G. De Saxcé, “Elastoplastic Finite Element Analysis of Soil Problems with Implicit Standard Material Constitutive Laws,” *Rev. Européenne Élém. Finis*, vol. 3, no. 3, pp. 411–456, Jan. 1994.

[30] V. Acary and B. Brogliato, *Numerical Methods for Nonsmooth Dynamical Systems*. Berlin, Germany: Springer, 2008. [Online]. Available: https://link.springer.com/book/10.1007/978-3-540-75392-6

[31] V. Acary, F. Cadoux, C. Lemaréchal, and J. Malick, “A formulation of the linear discrete Coulomb friction problem via convex optimization,” *Z. angew. Math. Mech.*, vol. 91, no. 2, pp. 155–175, Feb. 2011.

[32] F. Jourdan, P. Alart, and M. Jean, “A Gauss-Seidel like algorithm to solve frictional contact problems,” *Comput. Methods Appl. Mech. Eng.*, vol. 155, no. 1, pp. 31–47, Mar. 1998.

[33] Facchinei and Pang, *Finite-Dimensional Variational Inequalities and Complementarity Problems*. Springer, 2003.

[34] E. A. de Souza Neto, D. Perić, and D. R. J. Owen, *Computational Methods for Plasticity*. Wiley, Oct. 2008.

[35] J. C. Simo and T. J. R. Hughes, *Computational Inelasticity*. New York, NY, USA: Springer, 1998. [Online]. Available: https://link.springer.com/book/10.1007/b98904

[36] D. Drucker and W. Prager, “Soil mechanics and plastic analysis on limit design,” *Q. Appl. Math.*, vol. 10, no. 2, pp. 157–165, Jul. 1952, [Online; accessed 4. Mar. 2025]. [Online]. Available: https://www.jstor.org/stable/43633942

[37] M. Wojciechowski, “A note on the differences between Drucker-Prager and Mohr-Coulomb shear strength criteria,” *Studia Geotechnica et Mechanica*, vol. 40, no. 3, pp. 163–169, Nov. 2018.

[38] D. Loukidis and R. Salgado, “Bearing capacity of strip and circular footings in sand using finite elements,” *Comput. Geotech.*, vol. 36, no. 5, pp. 871–879, Jun. 2009.

[39] V. P. Nguyen, A. d. Vaucorbeil, and S. Bordas, *The Material Point Method*. Springer International Publishing, 2022.

[40] M. Cortis, W. Coombs, C. Augarde, M. Brown, A. Brennan, and S. Robinson, “Imposition of essential boundary conditions in the material point method,” *Int. J. Numer. Methods Eng.*, vol. 113, no. 1, pp. 130–152, Jan. 2018.

[41] F. Facchinei and J.-S. Pang, *Finite-Dimensional Variational Inequalities and Complementarity Problems*. New York, NY, USA: Springer, 2003. [Online]. Available: https://link.springer.com/book/10.1007/b97543

[42] V. Singer, T. Teschemacher, A. Larese, R. Wüchner, and K.-U. Bletzinger, “Lagrange multiplier imposition of non-conforming essential boundary conditions in implicit material point method,” *Comput. Mech.*, vol. 73, no. 6, pp. 1311–1333, Jun. 2024.
[43] C. Kanzow, I. Ferenczi, and M. Fukushima, “On the Local Convergence of Semismooth Newton Methods for Linear and Nonlinear Second-Order Cone Programs Without Strict Complementarity,” *SIAM J. Optim.*, Apr. 2009.

[44] W.-K. Chen, *Limit Analysis and Soil Plasticity*. Elsevier Science Ltd, Jul. 1975.

[45] J. Salencon, *Plasticite pour la Mecanique des Sols*. Springer, Jan. 1978.

[46] M. Hjiaj, A. V. Lyamin, and S. W. Sloan, “Numerical limit analysis solutions for the bearing capacity factor Nγ,” *Int. J. Solids Struct.*, vol. 42, no. 5, pp. 1681–1704, Mar. 2005.

[47] E. Love and D. L. Sulsky, “An energy-consistent material-point method for dynamic finite deformation plasticity,” *Int. J. Numer. Methods Eng.*, vol. 65, no. 10, pp. 1608–1638, Mar. 2006.
A Yield criterion in plane strain plasticity

Without any loss of generality, let us consider a frame coherent with the material principal directions and consider that the z direction corresponds to the larger dimension. The strain tensor becomes

$$ \epsilon = \begin{pmatrix} \epsilon_{11} & \epsilon_{12} & 0 \\ \epsilon_{12} & \epsilon_{22} & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (80) $$

The plane strain assumption also assumes that $\sigma_{23} = \sigma_{13} = 0$ but $\sigma_{33}$ is not null. It also entails that $0 = \epsilon_{33} = \epsilon_{33}^E = \epsilon_{33}^P$. We make the hypothesis that there is no plastic deformation in the third direction, meaning that $0 = \epsilon_{33} = \epsilon_{33}^E = \epsilon_{33}^P$. If we consider linear isotropic elasticity, we have the relation

$$ \sigma_{33} = \nu(\sigma_{11} + \sigma_{22}), \qquad (81) $$

where $\nu$ is the Poisson ratio. The stress invariants become

$$ I_1(\boldsymbol{\sigma}) = (1 + \nu)(\sigma_{11} + \sigma_{22}), \qquad (82) $$

$$ J_2(\boldsymbol{\sigma}) = \frac{1}{6}(\sigma_{11} - \sigma_{22})^2 + (\sigma_{11} - \nu(\sigma_{11} + \sigma_{22}))^2 + (\sigma_{22} - \nu(\sigma_{11} + \sigma_{11}))^2 + \sigma_{12}^2 \qquad (83) $$

$$ = \frac{1}{4}(\sigma_{11} - \sigma_{22})^2 + \frac{1}{3}(\nu - \frac{1}{2})^2(\sigma_{11} + \sigma_{22})^2 + \sigma_{12}^2. \qquad (84) $$

It is common [37] to set the Poisson's ratio equals to $\nu = \frac{1}{2}$ to define the stress invariant under the plane strain hypothesis.

$$ I_1(\boldsymbol{\sigma}) = \frac{3}{2}(\sigma_{11} + \sigma_{22}) \qquad (85) $$

$$ J_2(\boldsymbol{\sigma}) = \frac{1}{4}(\sigma_{11} - \sigma_{22})^2 + \sigma_{12}^2 = \left(\frac{\sigma_{11} - \sigma_{22}}{2}\right)^2 + \sigma_{12}^2. \qquad (86) $$

The parameters $\eta$ and $k$ of the Drucker-Prager yield criterion are then computed to match the plane strain Mohr-Coulomb criterion. If we denote $\sigma_1$ and $\sigma_2$, the principal stress components, the Mohr-Coulomb criterion can be written under the following form [44]

$$ \frac{\sigma_1 - \sigma_2}{2} \le \sin(\phi) \frac{\sigma_1 + \sigma_2}{2} + c \cos \phi. \qquad (87) $$

In the principal stress directions, the stress invariants are written

$$ I_1(\boldsymbol{\sigma}) = \frac{3}{2}(\sigma_1 + \sigma_2) \qquad (88) $$

$$ J_2(\boldsymbol{\sigma}) = \left(\frac{\sigma_1 - \sigma_2}{2}\right)^2. \qquad (89) $$

Therefore in order to match the Mohr-Coulomb criterion, the Drucker-Prager parameters $\eta$ and $k$ can be expressed depending on $\phi$ and $c$

$$ \eta = \frac{\sin \phi}{3}, \quad \text{and} \quad k = c \cos(\phi). \qquad (90) $$
Note that it is also possible to derive a plane strain model by considering only plastic deformations [36]. For non-associated plasticity, the coefficient $\eta_{\theta}$ will be defined as

$$ \eta_{\theta} = \frac{\sin \theta}{3} \qquad (91) $$

With $\theta$ the dialatency angle.

## B Retrieving a SOC using Voigt's notation

In order to use a solver derived from the convex cone programming, we express the yield criterion using Voigt's notation, associating $\sigma$ and $\varepsilon$ as vector of $\mathbb{R}^3$. We then denotes again $\sigma$ and $\dot{\varepsilon}^P$ the following vector of $\mathbb{R}^3$

$$ \sigma = \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{pmatrix} \qquad \dot{\varepsilon}^P = \begin{pmatrix} \dot{\varepsilon}_{11}^P \\ \dot{\varepsilon}_{22}^P \\ 2\dot{\varepsilon}_{12}^P \end{pmatrix} \qquad (92) $$

such that

$$ \sigma : \dot{\varepsilon}^P = \sigma \cdot \dot{\varepsilon}^P \qquad (93) $$

We then express the plasticity criteria and the plastic flow rule as second order cone (SOC) inclusion. Let us denote $n$ the normalised direction

$$ n = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix} \qquad (94) $$

the main idea is to express the two invariants $J_2(\sigma)$ and $I_1(\sigma)$ as a normal part and a tangential decomposition of $\sigma$. We first notice that

$$ \sigma_n := \sigma \cdot n = \frac{1}{\sqrt{2}} I_1 \sigma \qquad (95) $$

For the tangential part, a subtlety occurs, we have in deed $\sigma_t := \sigma - \sigma_n n$ thus

$$ \sigma_t = \begin{pmatrix} \sigma_{11} - \frac{1}{2}(\sigma_{11} + \sigma_{22}) \\ \sigma_{22} - \frac{1}{2}(\sigma_{11} + \sigma_{22}) \\ \sigma_{12} \end{pmatrix} = \frac{1}{2} \begin{pmatrix} \sigma_{11} - \sigma_{22} \\ \sigma_{22} - \sigma_{11} \\ 2\sigma_{12} \end{pmatrix} \qquad (96) $$

Then

$$ |\sigma_t|^2 = \frac{1}{4}(2(\sigma_{11} - \sigma_{22})^2 + 4\sigma_{12}^2) \qquad (97) $$

$$ = \frac{1}{2}(\sigma_{11} - \sigma_{22})^2 + \sigma_{12}^2 \qquad (98) $$

We notice that we don't have $|\sigma_t| \neq \sqrt{J_2(\sigma)} = \frac{1}{4}(\sigma_{11} - \sigma_{22})^2 + \sigma_{12}^2$. In order to deal with SOC, it is necessary to introduced a new variable

$$ \hat{\sigma} = \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ 2\sigma_{12} \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 2 \end{pmatrix} \sigma \qquad (99) $$
such that $\hat{\sigma}_n = \sigma_n$ and $|\hat{\sigma}_t|^2 = 2J_2(\sigma)$. we can now establish the equivalency

$$ \sigma \in \mathcal{C} \iff \hat{\sigma} \in \mathcal{K} \qquad (100) $$

$$ \text{where } \mathcal{K} = \{x \in \mathbb{R}^3, |x_t| \le 2\sqrt{2\eta}x_n\} $$

Regarding the plastic trian rate, no change of variable is needed, due to Voigt's notation, we directly have $\dot{\varepsilon}_n^P = \frac{1}{\sqrt{2}}I_1(\dot{\varepsilon}_n^P)$ and $|\dot{\varepsilon}_n^P|^2 = 2\sqrt{J_2(\dot{\varepsilon}_n^P)}$, which gives the following inclusion for the plastic flow rule

$$ \dot{\varepsilon}^P \in \mathcal{K}_\theta := \{x \in \mathbb{R}^3, |x_t| \le \frac{1}{2\sqrt{2\eta_\theta}} x_n\} \qquad (101) $$

## C Sub-gradient of natural map

in order to apply the semi-smooth Newton's method on the natural map, we have to compute its sub gradient with respect to $\sigma$ and $z$. We recall

$$ F_{NM}(\sigma, \tilde{z}) = \sigma - P_{\mathcal{K}}(\sigma - \rho \tilde{z}) \qquad (102) $$

and

$$ \Psi(z) = z + \mu|z_t|, \quad z_t = z - z_n \cdot n \qquad (103) $$

where $z_t, z_n, n$ are defined in appendix B. As $\Psi$ is also an non-smooth function and in order not to write the generalise jacobian of the composition of two non-smooth functions, we write the generalised jacobian of the function $G$ defined as

$$ G: \begin{pmatrix} \sigma \\ \tilde{z} \end{pmatrix} \mapsto \begin{pmatrix} \Psi(W\sigma + q) - \tilde{z} \\ F_{NM}(\sigma, \tilde{z}) \end{pmatrix} \qquad (104) $$

We then have that each element of $\partial G(\sigma, z)$ writes

$$ \begin{pmatrix} VW & -I \\ I - J & \rho J \end{pmatrix}, \qquad (105) $$

where $V$ is an element of $\partial\Psi(W\sigma + q)$ and $J$ is an element of $\partial P_{\mathcal{K}}(\sigma - \rho\tilde{z})$ and

$$ \partial\Psi(y) = \begin{cases} I + \frac{\mu}{|z_t|}nz_t^T, & \text{if } |z_t| > 0 \\ I + B(0, 1), & \text{else.} \end{cases} \qquad (106) $$

In equation (106), $B(0, 1)$ corresponds to the ball centred in 0 of radius 1. The first case corresponds to a singloton but for the sake of simplicity we identify it to its only element. The projection on a second order cone of opening $\eta$ and normal direction $n$ is given by

$$ P_{\mathcal{K}}(x) = \begin{cases} x, & \text{if } x \in \mathcal{K} \\ \left( \frac{\eta}{1+\eta^2} \left( \eta + \frac{x_n}{|x_t|} \right) \right) \left( \frac{1}{\eta} |x_t| n + x_t \right), & \text{if } x \notin \mathcal{K} \cup \mathcal{K}^* - \{0\} \\ 0, & \text{if } x \in \mathcal{K}^* \end{cases} \qquad (107) $$
on its regular part, the jacobian is given by

$$ \partial P_{\mathcal{K}}(x) = \begin{cases} I, & \text{if } x \in \mathring{\mathcal{K}} \\ \left( \frac{\eta}{1+\eta^2} \left( \frac{1}{|x_t|} n x_t^T + \frac{1}{\eta} n n^T - \frac{x_n}{|x_t|^3} x_t x_t^T + \frac{1}{|x_t|} x_t n^T + \left( \eta + \frac{x_n}{|x_t|} \right) (I - n n^T) \right), & \text{if } x \notin \mathring{\mathcal{K}} \cup \mathring{\mathcal{K}}^* \\ 0, & \text{if } x \in \mathring{\mathcal{K}}^* \end{cases} \tag{108} $$

where $\mathring{\mathcal{K}}$ denotes the interior of the set $\mathcal{K}$. When $x$ is at the border of the cone $\partial \mathcal{K}$ or its dual, the generalised jacobian $\partial P_{\mathcal{K}}(x)$ is the convex hull of the value given in case 1 and case 2 of equation (108)