# Semi-smooth Newton method for nonassociative plasticity using the bi-potential approach

Louis Guillet, Vincent Acary, Franck Bourrier, Olivier Goury

► To cite this version:

Louis Guillet, Vincent Acary, Franck Bourrier, Olivier Goury. Semi-smooth Newton method for nonassociative plasticity using the bi-potential approach. 16ème Colloque National en Calcul de Structures (CSMA 2024), CNRS; CSMA; ENS Paris-Saclay; CentraleSupélec, May 2024, Hyères, France. hal-04611037

HAL Id: hal-04611037

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

Submitted on 3 Dec 2024

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
# Semi-smooth Newton method for nonassociative plasticity using the bi-potential approach

L. Guillet<fnref n="1" />, V. Acary<fnref n="1" />, F. Bourrier <fnref n="1" />,<fnref n="2" />, O. Goury <fnref n="1" />

**Résumé** — This paper presents an implicit solver for non-associative plasticity problems based on the semi-smooth Newton method. The method is derived from the Implicit Standard Material and is easily compatible with various space discretization techniques, particularly the Material Point Method and the Finite Element Method. The solver converges quadratically, even for large time steps, although theoretical results have only been demonstrated for restricted cases. The method is exemplified through a footing simulation.

**Mots clés** — Non-Smooth Dynamics, Non-Associative Plasticity, Complementarity Based Solvers

## 1 Introduction and motivations

From footing simulation to slope stability analysis, geomaterials are often described by an elasto-plastic model. In order to simulate their correct behaviour, a non-associative flow rule must be added to the model. There are already several approaches and numerical methods to solve this model. The most common ones are the use of a Finite Element Method (FEM) [17] for the spatial discretisation and the Return Map Algorithm to solve the conservation equation and the constitutive laws. The long term goal of this work, in continuity with [2], is to develop a monolithic solver capable of solving the constitutive equation, the potential hardening and the Coulomb friction, in order to simulate complex phenomena such as debris flows or landslide initiation. This framework is compatible with various space discretisations, including the Material Point Method (MPM), which has been widely used in geomaterial simulations over the last two decades. Inspired by the Implicit Generalised Materials (IGM) [4], the constitutive equations are written as differential inclusion and a semi-smooth Newton method [3, 6, 10] is presented to solve the constitutive law in an implicit time scheme. For now, this work is limited to the simulation of non-associative elasto-plastic materials. In the first part, the IGM framework is developed, then the spatial discretisation (FEM or MPM) is derived before the non-smooth solver is presented.

## 2 Governing equations Standard Generalized Materials (SGM) and Implicit Generalized Materials (IGM)

If we restrict ourselves to the framework of stability analysis, we assume the small deformation hypothesis, where the strain tensor is :

$$ \varepsilon = \frac{1}{2} D(u), \qquad (1) $$

where the operator $D$ is the symmetric part of the gradient and $u$ is the displacement with respect to the position and time. The small deformation hypothesis is often combined with an additive elasto-plastic model with linear elasticity. The strain is then decomposed as

$$ \varepsilon = \varepsilon^e + \varepsilon^p \qquad (2) $$

$$ \sigma = \mathcal{D} : \varepsilon^e, \qquad (3) $$

where $\sigma$ is the Cauchy stress tensor and $\mathcal{D}$ is the fourth order elastic tensor. To complete the constitutive equations, we need to consider a yield criterion based on critical state theory and the plastic flow rule, which satisfies the second principle of thermodynamics, which states that the dissipation $d = \dot{\varepsilon}^p : \sigma$ must

<footnote n="1">¹ Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, Institute of Engineering, LJK, 38000, Grenoble, {louis.gillet,vincent.acary,olivier.goury}@inria.fr</footnote>

<footnote n="2">² Univ. Grenoble Alpes, INRAE, CNRS, IRD, Grenoble INP, IGE, 38000 Grenoble, France, franck.bourrier@inrae.fr</footnote>
always be non negative. The constitutive laws are written in the form of a differential inclusion or a variational inequality.

A standard framework for associative plasticity is the Generalized Standard Material (GSM) approach developed in [11] which assumes the existence of a convex pseudo-potential of dissipation V and its Fenchel transform W such that the dissipation inequality is

$$V(\sigma') + W(\dot{\varepsilon}^p) \geq \sigma' : \dot{\varepsilon}^p, \quad \text{for all } (\sigma', \dot{\varepsilon}^p) \text{ admissible.} \quad (4)$$

Using convex analysis to link the generalized forces (stress in this case) and the generalized velocity (plastic strain rate), we obtain $\sigma \in \partial W(\dot{\varepsilon}^p)$ and $\dot{\varepsilon}^p \in \partial V(\sigma)$ which satisfies the equality in (4).

While the GSM framework is very well suited for associative plasticity where it can be used in simple cases to solve non associative plasticity problems, a more general framework, the Implicit Generalized Materials (ISM), introduced in [4], is an efficient and elegant approach. It has already been used to handle this type of problems with a return map algorithm in [19]. The pseudo-potential V and its Fenchel dual W are replaced by a bipotential b, convex with respect to each variable separately, such that for all $(\dot{\varepsilon}^p, \sigma')$ admissible

$$b(\dot{\varepsilon}^p, \sigma') \leq \dot{\varepsilon}^p : \sigma'. \quad (5)$$

The equality in (5) is reached when

$$\sigma \in \partial b_{\dot{\varepsilon}^p}(\dot{\varepsilon}^p, \sigma) \text{ and } \dot{\varepsilon}^p \in \partial b_{\sigma'}(\dot{\varepsilon}^p, \sigma). \quad (6)$$

This formulation allows some coupling terms between $\sigma$ and $\dot{\varepsilon}^p$ in the bipotential formulation and thus allows the modeling of non associative plasticity flow rules.

In the remaining, we consider the Drucker-Prager yield criterion. The latter can be written as $\sigma \in \mathcal{K}$ where $\mathcal{K}$ is the translated closed convex cone defined by

$$\mathcal{K} = \{ \sigma = (\sigma^h, \operatorname{dev}(\sigma)) \text{ s.t. } \frac{1}{k_d} \|\operatorname{dev}(\sigma)\| + \sigma^h \tan\phi \le c \} = \{ \sigma \text{ s.t. } f(\sigma) \le 0 \}, \quad (7)$$

where we decompose any tensor field $\tau$ into its hydrostatic part $\tau^h$ and its deviatoric dev($\tau$) part. Since the real $k_d$ is only a normalization factor depending on the dimension, the yield criterion is determined by two parameters, the friction angle $\phi$ and the cohesion $c$. To simulate realistic behaviours of geomaterials, we need to be able to adjust the dilatancy angle $\theta$. The plastic strain rate $\dot{\varepsilon}^p = (\dot{\varepsilon}^{p,h}, \operatorname{dev}(\dot{\varepsilon}^p))$ and the stress must satisfy the differential inclusion

$$(\dot{\varepsilon}^{p,h} + k_d(\tan\phi - \tan\theta)\|\operatorname{dev}(\dot{\varepsilon}^p)\|, \operatorname{dev}(\dot{\varepsilon}^p)) \in \mathcal{N}_{\mathcal{K}}(\sigma^h, \operatorname{dev}(\sigma)) \quad (8)$$

where $\mathcal{N}_{\mathcal{K}}(x)$ correspond to the normal cone of C taken at the point x. To define the bipotential associated with the problem, we need to specify the admissible rate of plastic strain. Eq. (8) implies that the admissible set is

$$\mathcal{K}_{\varepsilon} = \{ \dot{\varepsilon}^{p,h} \le k_d \tan\theta \|\operatorname{dev}(\dot{\varepsilon}^p)\| \}. \quad (9)$$

The function

$$b(\dot{\varepsilon}^p, \sigma) := \psi_{\mathcal{K}_{\varepsilon}}(\dot{\varepsilon}^p) + \psi_{\mathcal{K}}(\sigma) + \frac{\dot{\varepsilon}^{p,h}c}{\tan\phi} + k_d \left( \sigma^h - \frac{c}{\tan\phi} \right) (\tan\theta - \tan\phi) \|\operatorname{dev}(\dot{\varepsilon}^p)\| \quad (10)$$

defines the bipotential corresponding to the plasticity model. The first two terms represent the fact that the couple $(\dot{\varepsilon}^p, \sigma)$ must be admissible (feasible) and the coupled term encodes the non-associative plastic flow rule.

## 3 Space Discretization

**Material Point Method** The FEM is appropriate for stability simulations under the small perturbation hypothesis, it may not accurately simulate large deformations caused by mesh distortion. The spatial discretisation can be performed using the Material Point Method (MPM). The MPM is a hybrid approach first introduced in [18], which couples Lagrangian particles with a background Eulerian mesh.
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Non associative flow rule with $s = \text{dev}(\sigma)$, taken from [4]</figcaption></figure>

The conservation equations are solved at the mesh nodes and the material quantities, including mass, velocity, and stress, are stored in the particles. This allows, due to the Lagrangian nature of the material particles, not to consider the convective term in the momentum conservation equation. An interpolation scheme, with the same type of shape function used in FEM, is used to transfer the quantities from the material points to the grid nodes and from the grid nodes to the particles. This method has been broadly used for applications in geosciences [9, 7, 16].

As described in [15], the MPM algorithm can be divided into three different steps. We denote by a subscript *i* the quantities over the grid and by a subscript *p* the quantities over the material points (or particles). The quantities are interpolated with the "shape function" of the fine support, centred at each grid node and denoted $N_i$. The subscript *n* corresponds to the known approximation at time $t_n$ and *h* is the fixed time step. We refer the reader to [15] or [13] for further details.

1. **Particles to Grid interpolation.** During this first phase, the mass and momentum are interpolated from the material particles to the grid nodes with

$$ (Mv)_i = \sum_p N_i(x_p) m_p v_p \quad (11) $$

The mass on grid nodes $m_i$ are also computed using a lumped mass matrix on the grid

2. **Mesh Dynamics.** During this phase the momentum conservation equation are solved, we then have

$$ \frac{(Mv)_{n+1} - (Mv)_n}{h} = B^T \sigma_p + f_{n+1} \quad (12) $$

where B is a sort of discrete divergence operator in the same sense as in [5]. The equation (12) is coupled with the constitutive law and plastic flow rule to solve the whole system with the stress at the particles $\sigma_p$. The method and the full system are described in details below.

3. **Particles update.** The velocity and the position of the particles are then updated with the new grid quantities obtained by solving (12). Several schemes can be used to update the particles (PIC, FLIP, APIC ...), the most common being a combination between the PIC and FLIP schemes. For all particles $p$, we obtain

$$ v_{p,n+1} = \alpha \sum_i N_i(x_p) v_{i,n+1} + (1-\alpha) \left( v_{p,n} + h \sum_i N_i(x_p) (v_{i,n+1} - v_{i,n}) \right) \quad (13) $$

$$ x_{p,n+1} = x_{p,n} + h v_{p,n+1} \quad (14) $$

with $\alpha$ being the so-called PIC-FLIP ratio often chosen between 0.95 and 0.99.
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Schematics representing differences between FEM and MPM discretization</figcaption></figure>

In the realm of finite strain theory, research by Love et al. [15] has shown that if the plastic law satisfies certain discrete dissipation inequalities, the numerical method becomes dissipative. The mass conservation equation is solved at the end of each time step.

**Resolution of the mesh dynamics** By introducing the slack variable $z = -\dot{\epsilon}^p$ the mesh dynamics and the constitutive equation become equivalent to : find $z, \sigma, v$ such that

$$ egin{cases} M \frac{v - v_n}{h} + B^T \sigma = f_{n+1} \\ S \frac{\sigma - \sigma_n}{h} = Bv + z \\ G(-z, \sigma) = 0. \end{cases} \qquad (15) $$

In the previous system $S$ is the inverse of the linear elastic operator and the function $G$ is defined by $G(\dot{\epsilon}^p, \sigma) = 0$ if and only if $(\dot{\epsilon}^p, \sigma)$ satisfies the plastic flow rule described in the previous section. By setting $W = S + h^2 BM^{-1}B^T$ and $q = -(S\sigma_n - hBv_n + h^2 BM^{-1}f_{n+1})$, the system reduces to

$$ \left\{ \begin{array}{l} hz = W\sigma + q \\ G(-z, \sigma) = 0. \end{array} \right. \qquad (16) $$

**Associative plasticity.** Assuming $\theta = \phi$, the system is equivalent to :

$$ -(W\sigma + q) \in \mathcal{N}_c(\sigma) \qquad (17) $$

which is the first order optimality condition of the following quadratic convex minimization problem :

$$ \min_{\sigma} \quad \sigma^T W \sigma + q^T \sigma \\ \sigma \in \mathcal{C} \qquad (18) $$

since $\mathcal{C}$ is convex and $W$ is positive definite (with is always the case with a sufficiently small time step), that ensures that the discrete minimization problem (18) is well posed.
As the vector field $\sigma$ is the concatenation of the stress at each interpolation point (Gauss point or material point) denoted $\sigma^p$, the set $\mathcal{C}$ can be seen as the intersection of all admissible sets, we have thus $\sigma \in \mathcal{C}$ if and only if for all interpolation points, $\sigma^p \in \mathcal{K}$. Another equivalent formulation of (17), based on the natural map [8] is

$$ \sigma = P_c(\sigma - \rho(W\sigma + q)), \quad \forall \rho > 0, \qquad (19) $$

where $P_c$ is the projection operator onto the set $\mathcal{C}$. This will be used for solving the system with semi-smooth Newton method.

**Non associative case** Even though the problem can not be written on the form of an optimality condition, the solution of problem in (8) will verify a problem similar to (19) where we need to find $z, \sigma$ such that

$$ \left\{ \begin{array}{l} z = W\sigma + q \\ \sigma = P_c\left(\sigma - \rho\left(z^h + \alpha\|\mathrm{dev}(z)\|, \mathrm{dev}(z)\right)\right) \end{array} \right. \qquad (20) $$

The parameter $\alpha$ is equal to $\tan(\phi) - \tan(\theta)$ and quantifies the dilatancy.

**Resolution of the system** Both systems (19) and (20) can be put into the form : find $\sigma$ such that $F(\sigma) = 0$ and can therefore be solved using the semi-smooth Method as described in [8] or [12]. Supposing that we found the k-th iterate $\sigma_k$, the descent direction $d_k$ is computed as

$$ d_k = -V_k^{-1}\sigma_k, \quad V \in \partial F(\sigma_k) \qquad (21) $$

$$ \sigma_{k+1} = \sigma_k + d_k, \qquad (22) $$

where $V$ is a non singular element of the Clarke derivative of $F$ at the point $\sigma_k$, denoted $\partial F(\sigma_k)$. This set includes the convex hull of all directional derivatives. The method converges if $\|F(\sigma_k)\|$ reaches a given user tolerance. This method has the advantage of having a quadratic convergence but the Clarke derivative of the semi-smooth function needs to be known. Fortunately, when $\mathcal{K}$ is a translated second order cone as in Drucker-Prager plasticity, the explicit form of the projection and its gradient are known. If we denote $\mathcal{K}_\eta$ the cone of angle arctan($\eta$) and oriented with the normal $n$, the projection onto the cone is given by

$$ P_{\mathcal{K}_\eta}(x) = \begin{cases} x, & \text{if } x \in \mathcal{K}_\eta \\ \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}_\eta \cup \mathcal{K}_\eta^* - \{0\} \\ 0, & \text{if } x \in \mathcal{K}_\eta^* \end{cases} \qquad (23) $$

where $x_n = x \cdot n$, $x_t = x - x_n n$ and $\mathcal{K}_\eta^*$ is the dual cone of $\mathcal{K}_\eta$. The projection is then given by

$$ JP_{\mathcal{K}_\eta}(x) = \begin{cases} I, & \text{if } x \in \mathcal{K}_\eta \\ \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 \mathcal{K}_\eta \cup \mathcal{K}_\eta^* - \{0\} \\ 0, & \text{if } 0 \in \mathcal{K}_\eta^* \end{cases} \qquad (24) $$

The Jacobian of the function for which we are searching the zero values is then given in the associative case by

$$ JF(\sigma) = I - JP_c(\sigma - \rho(W\sigma + q))(I - \rho W) \qquad (25) $$

Furthermore, existing local convergence results can be used in the associative scenario [14], whereby it is demonstrated that for a chosen $\rho$, all elements of $\partial F(\sigma^*)$ are non-singular at the solution $\sigma^*$. However, note that this convergence is solely local, meaning that the chosen initial value $\sigma_0$ must be sufficiently close to the solution $\sigma^*$ for the method to converge. For the non-associative scenario, it is possible to demonstrate that there exists a small value $\delta > 0$, where the condition $\alpha \le \delta$ is adequate to achieve local convergence. Nevertheless, in the non-convex case, our simulation indicates that the method converges with sufficiently small time steps for any $\alpha$ value in our simulation.
## 4 Applications : Footing simulation

The presented method can be illustrated on a small example of footing simulation using a FEM discretization. The schematics of the simulation and the parameter are presented in Fig. 3. The set of parameters is inspired from [4]. The dilatancy angle $\theta$ is chosen such that $\tan \theta = \frac{\tan \phi}{2}$.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Left : schematics of the simulation. Right : parameters chosen for the simulation.</figcaption></figure>

The simulation results are shown in Fig.4a, displaying the cumulative plastic strains norm. The observed outcome is an expected plastic shear band under the loading zone.

Additionally, Fig.4b highlights the force plotted against displacement, where the force increases linearly before reaching a plateau caused by plastic deformation. It is important to note that the simulation exaggerates the displacement to showcase the shear phenomena. The Newton method converges towards a tolerance of $10^{-12}$ in a maximum of five iterations.

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Footing simulation</figcaption></figure>
# 5 Conclusion and outlooks

In this paper, we introduce an implicit time solver for non-associative plasticity that is compatible with several spatial discretisations, including MPM and FEM. To illustrate the solver, we apply it to a FEM footing simulation and observe a quadratic convergence, even for large time steps. We expect that this solver will serve as a useful tool for characterising MPM in small deformations and for performing larger-scale simulations of landslide initiation using this discretisation.

Using the same method as [2], we can incorporate isotropic and/or kinematic hardening into our model and maintain its well-posedness by creating an internal variable, $\alpha$, for complementary conditions. The plastic flow rule outlined here is highly similar to the Coulomb friction problem, as explained in [1], so the method can be easily modified for coupled elasto-plasticity with frictional contact. As in [2], the aim is to incorporate one-sided contact conditions, which in the dynamics of dimensional systems means including a plastic impact law as well as Coulomb friction. The idea is to provide a monolithic solver and a time integration scheme that guarantees the energetic aspects in discrete time. Finally, the aim of this work is to integrate this numerical method into the MPM with a strong coupling to the non-smooth DEM (NCSD), in order to simulate elasto-plastic flows with rigid particles within a single solver.

# Références

[1] V. ACARY et al. "A formulation of the linear discrete Coulomb friction problem via convex optimization". In: *Z. angew. Math. Mech.* 91.2 (fév. 2011), p. 155-175. ISSN: 0044-2267. DOI: 10.1002/zamm.201000073.

[2] Vincent ACARY, Franck BOURRIER et Benoit VIANO. *Variational approach for nonsmooth elasto-plastic dynamics with contact and impacts*. [Online; accessed 4. Jul. 2023]. Fév. 2023. URL: https://cnrs.hal.science/LJK/hal-03978387v1.

[3] P. ALART et A. CURNIER. "A mixed formulation for frictional contact problems prone to Newton like solution methods". In: *Comput. Methods Appl. Mech. Eng.* 92.3 (nov. 1991), p. 353-375. ISSN: 0045-7825. DOI: 10.1016/0045-7825(91)90022-X.

[4] Abdelmajid BERGA et Géry DE SAXCÉ. "Elastoplastic Finite Element Analysis of Soil Problems with Implicit Standard Material Constitutive Laws". In: *Rev. Européenne Élém. Finis* 3.3 (jan. 1994), p. 411-456. ISSN: 1250-6559. DOI: 10.1080/12506559.1994.10511137.

[5] Jeremy BLEYER. "Automating the formulation and resolution of convex variational problems : applications from image processing to computational mechanics". In: *arXiv* (nov. 2019). DOI: 10.48550/arXiv.1911.13185. eprint : 1911.13185.

[6] Peter W. CHRISTENSEN. "A nonsmooth Newton method for elastoplastic problems". In: *Comput. Methods Appl. Mech. Eng.* 191.11 (jan. 2002), p. 1189-1219. ISSN: 0045-7825. DOI: 10.1016/S0045-7825(01)00321-8.

[7] Sachith DUNATUNGA et Ken KAMRIN. "Continuum modelling and simulation of granular flows through their many phases". In: *arXiv* (nov. 2014). DOI: 10.1017/jfm.2015.383. eprint : 1411.5447.

[8] Francisco FACCHINEI et Jong-Shi PANG. *Finite-Dimensional Variational Inequalities and Complementarity Problems*. New York, NY, USA: Springer, 2003. ISBN: 978-0-387-21814-4. URL: https://link.springer.com/book/10.1007/b97543.

[9] J. GAUME et al. "Dynamic anticrack propagation in snow". In: *Nat. Commun.* 9.3047 (août 2018), p. 1-10. ISSN: 2041-1723. DOI: 10.1038/s41467-018-05181-w.

[10] Corinna HAGER et B. I. WOHLMUTH. "Nonlinear complementarity functions for plasticity problems with frictional contact". In: *Comput. Methods Appl. Mech. Eng.* 198.41 (sept. 2009), p. 3411-3427. ISSN: 0045-7825. DOI: 10.1016/j.cma.2009.06.021.

[11] Bernard HALPHEN et Quoc Son NGUYEN. "Sur les matériaux standard généralisés". In: *Journal de Mécanique* 14.1 (1975), p. 39-63. URL: https://hal.science/hal-03600755.
[12] M. HINTERMÜLLER, K. ITO et K. KUNISCH. "The Primal-Dual Active Set Strategy as a Semismooth Newton Method". In : *SIAM J. Optim*. (2003). URL : https://epubs.siam.org/doi/10.1137/S1052623401383558.

[13] Chenfanfu JIANG, Craig SCHROEDER et Joseph TERAN. "An angular momentum conserving Affine-Particle-In-Cell method". In : *arXiv* (mars 2016). DOI : 10.1016/j.jcp.2017.02.050. eprint : 1603.06188.

[14] Christian KANZOW, Izabella FERENCZI et Masao FUKUSHIMA. "On the Local Convergence of Semismooth Newton Methods for Linear and Nonlinear Second-Order Cone Programs Without Strict Complementarity". In : *SIAM J. Optim*. (avr. 2009). URL : https://epubs.siam.org/doi/10.1137/060657662.

[15] E. LOVE et D. L. SULSKY. "An energy-consistent material-point method for dynamic finite deformation plasticity". In : *Int. J. Numer. Methods Eng.* 65.10 (mars 2006), p. 1608-1638. ISSN : 0029-5981. DOI : 10.1002/nme.1512.

[16] Vinh Phu NGUYEN, Alban de VAUCORBEIL et Stephane BORDAS. *The Material Point Method*. Cham, Switzerland : Springer International Publishing, 2023. ISBN : 978-3-031-24070-6. URL : https://link.springer.com/book/10.1007/978-3-031-24070-6.

[17] J. C. SIMO, J. G. KENNEDY et R. L. TAYLOR. "Complementary Mixed Finite Element Formulations for Elastoplasticity". In : *Computer Methods in Applied Mechanics and Engineering* 74.2 (1er sept. 1989), p. 177-206. ISSN : 0045-7825. DOI : 10.1016/0045-7825(89)90102-3. URL : https://www.sciencedirect.com/science/article/pii/0045782589901023 (visité le 23/09/2021).

[18] D. SULSKY, Z. CHEN et H. L. SCHREYER. "A particle method for history-dependent materials". In : *Comput. Methods Appl. Mech. Eng.* 118.1 (sept. 1994), p. 179-196. ISSN : 0045-7825. DOI : 10.1016/0045-7825(94)90112-0.

[19] Yang-Jing ZHOU et al. "Application of the Bipotential Theory to a Nonassociated Drucker-Prager Model". In : *Adv. Civ. Eng.* 2018 (mai 2018). ISSN : 1687-8086. DOI : 10.1155/2018/4534750.