## 2 Plasticity Formulation
### Return Mapping Algorithm
#### 2.2.1 Return-Mapping Scheme for the 1D Model
The volumetric plastic strain increment \(\Delta\varepsilon_{v}^{p}\) can still be written as a function of the hydrostatic pressure p using equations Equation (15) and Equation (48). ### The Elastoplastic Consistent Tangent Operator

The elastoplastic consistent tangent operator is closely related to the algorithmic form for updating the stress tensor. Indeed, it is obtained by performing a linearization of the stress updating algorithm [1; 4; 20]. Thus, in this section, we present a numerical method using a first-order finite difference scheme for calculating this tangent operator. However, the analytical method is presented in Section 2.3.1 for 1D case and in the Appendix A for 2D case. Analytical calculation is often laborious, and sometimes even impossible to perform exactly. In the numerical approach, the return mapping algorithm is used twice. The additional numerical calculations come with more time consumption, but the complexities associated with the analytical calculations are overcome. Furthermore, the computational overload due the finite difference algorithm can be substantially reduced using parallel programing. The general expression of the consistent elastoplastic tensor is given by:

\[\mathbf{D}^{ep}=\frac{d\mathbf{\sigma}_{n+1}}{d\mathbf{\varepsilon}_{n+1}^{\text{ \emph{trial}}}} \tag{51}\]

where \(\sigma_{n+1}\) represents the stress tensor at time \(t_{n+1}\). #### 2.3.1 1D Illustration

Using the consistency condition \(\left(\dot{f}=0\right)\), the expression of \(\dot{\sigma}\) is given by:

\[\dot{\sigma}=E\bigg{(}\dot{\varepsilon}-\frac{E}{(E+H+K)}\dot{\varepsilon} \bigg{)}, \tag{52}\]

and so, the relationship of the constitutive law linking \(\dot{\sigma}\) to \(\dot{\varepsilon}\) becomes:

\[\dot{\sigma}=\left\{ \begin{array}{cc}\frac{E(H+K)}{(E+H+K)}\dot{\varepsilon} &if\ \ \dot{\gamma}>0\\ E\dot{\varepsilon}&if\ \ \dot{\gamma}=0\end{array} \right. \tag{53}\]We can therefore deduce the tangent modulus \(D^{\mathit{elp}}\) obtained analytically for the 1D model:

\[D^{\mathit{elp}}=\frac{E(H+K)}{(E+H+K)} \tag{54}\]

The tangent module is in general not easy to determine. Indeed, in this case, it was only made possible because of the linear shape of the yield loading function. However, if this function is non-linear, and for higher dimension problems, the calculations become tedious or even impossible and a numerical method may be necessary. Thus, at a time \(t_{n+1}\) corresponding to the load increment and iteration \(k\) of the Newton-Raphson algorithm solving the equations of equilibrium, the elastoplastic tangent modulus \(D_{n+1}^{(k)}\) is:

\[D_{n+1}^{\mathit{elp}\;(k)}=\frac{\partial\sigma_{n+1}^{(k)}}{\partial\epsilon _{n+1}^{(k)}} \tag{55}\]

In the 1D case, the analytical derivation leads to Equation (54). However, this result is not general [4]. The tangent elastoplastic tensor can be obtained in the triaxial case from the yield function expressed in terms of principal stresses. It is therefore necessary to bring it back to the Cartesian basis. The following approximation for a small \(\mathit{de}_{n+1}^{(k)}\) can be used:

\[\frac{\partial\sigma_{n+1}^{(k)}}{\partial\epsilon_{n+1}^{(k)}}\approx\ \frac{ \sigma_{n+1}^{(k)}\Big{(}\epsilon_{n+1}^{(k)\;\mathit{per}}\Big{)}-\sigma_{



#### 2.3.2 Case of Mohr-Coulomb's Law

The numerical method of the elastoplastic tangent matrix presented in the 1D case is generalized for the in-plane strain and for Mohr-Coulomb's law. The elastoplastic tangent matrix is approximated using a first order finite difference scheme for each component of:

\[\left(\mathbf{D}^{elp}\right)_{(n+1)}^{(k)}=\left(