Numerical Algorithms for Elastoplacity: Finite Elements Code Development and Implementation of the Mohr–Coulomb Law

: Two numerical algorithms for solving elastoplastic problems with the ﬁnite element method are presented. The ﬁrst deals with the implementation of the return mapping algorithm and is based on a ﬁxed-point algorithm. This method rewrites the system of elastoplasticity non-linear equations in a form adapted to the ﬁxed-point method. The second algorithm relates to the computation of the elastoplastic consistent tangent matrix using a simple ﬁnite difference scheme. A ﬁrst validation is performed on a nonlinear bar problem. The results obtained show that both numerical algorithms are very efﬁcient and yield the exact solution. The proposed algorithms are applied to a two-dimensional rockﬁll dam loaded in plane strain. The elastoplastic tangent matrix is calculated by using the ﬁnite difference scheme for Mohr–Coulomb’s constitutive law. The results obtained with the developed algorithms are very close to those obtained via the commercial software PLAXIS. It should be noted that the algorithm’s code, developed under the Matlab environment, offers the possibility of modeling the construction phases (i.e., building layer by layer) by activating the different layers according to the imposed loading. This algorithmic and implementation framework allows to easily integrate other laws of nonlinear behaviors, including the Hardening Soil Model. and the results compared to those obtained using the commercial software PLAXIS. The results of the different approaches remain in agreement with those given by PLAXIS, with deviations of around 2% and 0.1% in displacements for ux and uy, and of a maximum of 11% in stresses. These results are very encouraging, especially since we do not need to analytically calculate the often-tedious tangent matrix, thereby conﬁrming the effectiveness of the numerical models developed in this work.


Introduction
The use of numerical codes that account for plasticity is essential when designing geotechnical structures, as they are functional decision-making tools. However, the implementation of a good plasticity model requires the development of powerful algorithms to resolve the numerical difficulties arising from complex plasticity laws such as those used for soil modeling. For example, de Souza et al. [1] present several algorithms for solving plasticity problems (that are incorporated into the HYPLAS calculation code written in FORTRAN 77). The classical Newton-Raphson scheme is very often used to solve the nonlinearities of the plasticity models since it converges at the second order, but under the condition that the initial solution is close to the true solution. Its implementation requires a rewriting of the load function as a function of the plastic deformation increment. Furthermore, it should be ensured that during the iterative procedure, this plastic increment is not less than zero, which could happen for certain plasticity models. However, the Newton-Raphson scheme can be improved by using the Homotopy-Newton-CPPM or CG-Newton-CPPM algorithms proposed by Dajiang et al. [2]. Scherzinger et al. [3] proposed a line-search algorithm with the Newton-Raphson method making it possible to force the convergence when the first guessed solution is far from the sought numerical solution. Similarly, Simo and Hughes [4] described the theoretical foundations of several formulations and algorithms for plasticity based on the Newton-Raphson procedure. In particular, they show that the construction of the elastoplastic tangent matrix, which is closely related to the rate of convergence in the Finite Element Method (FEM) and guarantees the consistency of the algorithmic integration process [5], is a crucial step for the robustness of these algorithms, hence the difficulty of implementing certain constitutive laws (such as the Hardening Soil Model law [6]). However, for the perfectly plastic behavior laws, particularly Mohr-Coulomb's law, some authors have been able to establish this matrix analytically, as with de Souza et al. [1], who determined this matrix in the principal basis before returning it to the global basis (Cartesian) using spectral projection. Indeed, the arising functions and their derivatives expressed in the global basis are complicated, and so the eigenvectors must therefore be determined numerically [7]. This difficulty can be avoided by keeping the stresses in the main base when the plasticity or failure criterion has a very simple shape compared to its formulation in the general stress space [8]. Another approach was proposed in [9], where the authors preferred to express the yield function in a non-invariant form. This enables the derivatives and the plastic matrix to be expressed in a fairly simple algebraic form. In general, great care should be taken to avoid errors in deriving the consistent elastoplastic tangent analytically, as the calculations can become very complicated. Recently, another approach for solving local constitutive equations using conic optimization is discussed [10]. This method avoids blind guessing in the search for the elastic stress tensor (trial stresses) for the constitutive scheme. Indeed, the solution of the problem can be searched as the optimal pair (isotropic and linear behavior and linear hardening) for the convex optimization problem. Moreover, this procedure also allows to lighten the construction of the tangent operator. This operator can also be decomposed into elastoplastic and elastic tangent operators [11].
Commercial finite element software systems are available to solve geotechnical problems, including PLAXIS, ZSOIL and FLAC. However, coupling them with calibration or uncertainty analysis algorithms can be difficult due to the impossibility of accessing source codes. In addition, these software packages are expensive, and academic versions often do not offer a complete library that can meet specific needs; hence the need to develop an in-house code for research purposes. Most software packages were initially designed to meet a specific purpose before being extended to other areas. For example, the PLAXIS software we use to validate our results was initially designed for the analysis of river embankments on earthen soils before being extended to most other geotechnical fields [12].
We develop a finite element code for elastoplasticity calculations with the application of the Mohr-Coulomb soil law. This law is widely used in geotechnical engineering to model the behavior of soils, particularly those of earth dams. We advocate the use of a 'purely' numerical strategy to compute the consistent tangent matrix. For this purpose, we evaluate a fixed point and a finite difference scheme. These schemes are simple alternatives to existing methods for solving the return mapping. The resulting finite element code offers the flexibility to integrate other much more complex constitutive laws, particularly those whose plasticity criteria evolve during loading (such as the Hardening Soil Model [6]). In addition, it can be easily coupled with parameter calibration algorithms.
The organization of the rest of this paper is as follows: Section 2 presents the formulation of the plasticity problem, as well as the return mapping algorithm for 1D and 2D models. Section 3 presents and discusses numerical examples, and we draw our conclusions in Section 4.

Constitutive Law of Plasticity
In linear elasticity, the constitutive law is given by: where in plane strain case, σ = σ xx σ yy σ xy T is the stress vector, ε e = ε e xx ε e yy ε e xy T the elastic strain vector and D the matrix of elasticity material properties. As soon as the yield strength is exceeded during loading, and under the assumption of small strains, permanent strains (or plastic strains) ε p occur. The strain tensor ε is decomposed as: Appl. Sci. 2021, 11, 4637 3 of 27 and the stress-strain relationship is written as: This relationship can also be written in terms of tensor rates (or in incremental form linking the increments strains to those of the stresses): where . γ represents the plastic multiplier and the dots symbol means the pseudo-time derivative of the quantity.
The in-plane strain, the elasticity tensor linking stresses and strains in the principal basis are expressed as [1,13]: where K, G > 0 represents the elastic and shear modulus, respectively, I = 1 2 δ ik δ jl + δ il δ jk e i ⊗ e j ⊗ e k ⊗ e l is the fourth symmetric order identity tensor and 1 = δ ij e i ⊗ e j is the second order of identity tensor.

Return Mapping Algorithm
In plasticity problems, the constitutive laws linking stresses to strains, as well as the internal hardening law or that of the plastic multiplier (Equations (4) and (5)) are given as rates and therefore depend on the loading history. Thus, all these laws will have to be integrated over time with a numerical scheme. However, there are several numerical integration methods, including implicit and explicit Euler methods and Runge-Kutta methods. The choice of one of these methods should be simple; any one of them should provide accurate and robust results. The implicit Euler method also called the return-mapping method for elastoplastic problems and is often selected because of its stability [1,4,14,15]. This integration technique was first developed by Wilkins [16] for plane strain and elastoplasticity problems for low loads. Its extension to problems of large strains has been proposed by various authors, including Nagtegaal [17] and Simo and Taylor [18]. The iterative procedure often adopted for solving the constitutive problem is the classic Newton-Raphson scheme, which proves to be an optimal choice because of its quadratic convergence rate, generally leading to very efficient procedures in terms of the number of iterations required for convergence [1]. However, before applying this procedure, it is necessary, from the consistency condition, to express an equation dependent on the incremental plastic multiplier ∆γ, as presented in [1,4]. The implementation of the return mapping resolution algorithm is laborious and often analytically impossible when the yield function f is nonlinear. Therefore, we are developing a method for solving a nonlinear system based on the fixed-point algorithm. This system (Equation (7)) is generally defined by [1]: where at time n + 1, ε e n+1 represents the elastic strain tensor, ε e trial n+1 is the elastic strain test tensor; α n+1 is the internal hardening, N = ∂ f ∂σ the plasticity flux vector, A n+1 is all the combined thermodynamic forces associated with α, H = ∂ f ∂A is the generalized hardening module f represents the yield function and ∆γ the incremental plastic multiplier.

Return-Mapping Scheme for the 1D Model
To illustrate the proposed approach, we consider a 1D model, where we define the yield function f as obeying the mixed hardening law: where σ Y represents the yield strength, K the plastic modulus and α the internal hardening variable, which has the evolution equation: γ (9) and q is the back stress defined by Ziegler's rule as: with H representing the kinematic hardening modulus. For a 1D problem, the return mapping algorithm is explicitly presented in the work of Simo and Hughes [4]. It is reproduced here (Algorithm 1) with a reformulation as an algebraic system (Equation (11)) in order to highlight our numerical method. Algorithm 1: Return-mapping algorithm for 1D.
1. ε p n , α n and q n are known in step n 2.
Compute elastic trial stress σ trial n+1 and test for plastic loading 4.
Elseplastic step, solve the system : The analytical resolution of the system Equation (11) gives [4]: While there are several numerical methods for solving non-linear equation systems, we chose a simple fixed-point method. This choice is motivated by the fact that it does not require the calculation of the derivatives. Integrating Equation (11) amounts to calculating the mechanical state Z n+1 = σ n+1 , α n+1 , ∆γ, q n+1 , ε p n+1 at the instant t n+1 from a given strain increment ∆ε and a state Z n at instant t n . The initial values for Z 0 are σ trial n+1 , α n , 0, q n , ε p n . The iteration is stopped for a time step n + 1 ≤ N max when Z n+1 − Z 0 / Z 0 is tess then a fixed threshold. Then, once Z n+1 is calculated it is tested on the plasticity criterion as presented in [1]. If f trial n+1 ≤ 0 then the predicted state is the solution of the problem. Otherwise, the predicted state is not admissible, and a plastic correction step is performed. Algorithm 2 illustrates this iterative process for the 1D model: Algorithm 2: Fixed point method applied to the return mapping algorithm.

Return-Mapping Scheme for the Mohr-Coulomb (MC) Law
The MC model generally depends on two parameters: the internal friction angle φ of the material and the cohesion c, which can also depend on the accumulated plastic strain ε p n and its increment ∆ε p (c ε p n , ∆ε p . It also involves the shear stress τ and normal stress σ n . For this criterion, the plastic strain begins when the shear stress τ, as well as the normal stress σ n become critical. These two stresses are linked by the relationship: The plastic function f is expressed in terms of principal stresses by: Appl. Sci. 2021, 11, 4637 6 of 27 with the convention σ 1 ≥ σ 2 ≥ σ 3 where σ j (j = 1, 2, 3) is the j-th principal stress. The representation of its yield surface is a hexagonal pyramid whose axis is colinear to the hydrostatic pressure ( Figure 1) defined by: Appl. Sci. 2021, 11, x FOR PEER REVIEW 6 of 29 ( , ) = 1 − 3 + ( 1 + 3 ) sin − 2 cos (14) with the convention 1 ≥ 2 ≥ 3 where ( = 1,2,3) is the j-th principal stress. The representation of its yield surface is a hexagonal pyramid whose axis is colinear to the hydrostatic pressure ( Figure 1) defined by: The two-by-two intersections of these load surfaces and the vertex form singularities that pose some problems for numerical integration. Each of the six load surfaces has a plastic flow rule associated with it. In addition, the behavior of the soil during loading involves friction, and so non-associated plastic flow rules must be used to avoid excessive expansion of the material [9]. For this purpose, an expansion angle, noted as ψ, that allows the introduction of six plastic flow potentials (I = 1, …, 6) associated with each yield function . These plastic potentials are written in the same way as their respective yield surfaces (the relationships in (Equation (16)).
For example, for the yield surface 1 ( , ), we have: The plastic flow rule is then written as: The six edges of this criterion in three-dimensional space, allow it to be expressed by six load functions given by [1]: The two-by-two intersections of these load surfaces and the vertex form singularities that pose some problems for numerical integration. Each of the six load surfaces has a plastic flow rule associated with it. In addition, the behavior of the soil during loading involves friction, and so non-associated plastic flow rules must be used to avoid excessive expansion of the material [9]. For this purpose, an expansion angle, noted as ψ, that allows the introduction of six plastic flow potentials g i (I = 1, . . . , 6) associated with each yield function f i . These plastic potentials are written in the same way as their respective yield surfaces (the relationships in (Equation (16)).
For example, for the yield surface f 1 (σ, c), we have: The plastic flow rule is then written as: where we recall that m is the number of active mechanisms. Thus, with the principal stresses ordered as follows: σ 1 ≥ σ 2 ≥ σ 3 , the flow rule can be formulated in the edge space located in the same plane of the principal stresses space (the sextant plane) ( Figure 2). There are then four distinct possibilities for defining the plastic flow rule: If the final stress is within the load surface, then the point is regular and a mechanism is activated, m = 1.
If the final stress is located on the right edge of the cone relative to the sextant plan, then the point is singular and two mechanisms are activated, m = 2.
If the final stress is located on the left edge of the cone relative to the sextant plan, then the point is singular and also two mechanisms are activated, m = 2.
If the final stress is not located inside the load surface or on an edge, then it is projected to the top of the cone, the point is singular and six mechanisms are activated, m = 6.
(the sextant plane) ( Figure 2). There are then four distinct possibilities for defining the plastic flow rule:  If the final stress is within the load surface, then the point is regular and a mechanism is activated, m = 1.
 If the final stress is located on the right edge of the cone relative to the sextant plan, then the point is singular and two mechanisms are activated, m = 2.
 If the final stress is located on the left edge of the cone relative to the sextant plan, then the point is singular and also two mechanisms are activated, m = 2.
 If the final stress is not located inside the load surface or on an edge, then it is projected to the top of the cone, the point is singular and six mechanisms are activated, m = 6. Figure 2 illustrates the plastic flow rule of Mohr-Coulomb's law. Figure 2a shows how, according to the condition 1 ≥ 2 ≥ 3 linking the principal stresses, the main plane and the two possible configurations for two mechanisms are activated simultaneously: the right-edge mechanism or the left-edge mechanism and the apex (Figure 2b).


For a single activated mechanism, the plastic flow rule is written: with the flow vector normal to the corresponding plane defined by 1 = 0, given by:  For two activated mechanisms, the plastic flow rule is written:   Figure 2a shows how, according to the condition σ 1 ≥ σ 2 ≥ σ 3 linking the principal stresses, the main plane and the two possible configurations for two mechanisms are activated simultaneously: the right-edge mechanism or the left-edge mechanism and the apex (Figure 2b).
For a single activated mechanism, the plastic flow rule is written: γ N a (19) with N a the flow vector normal to the corresponding plane defined by f 1 = 0, given by: For two activated mechanisms, the plastic flow rule is written: with N b the normal at the intersection of the two planes f 1 = 0 and f 6 = 0 or f 1 = 0 and f 2 = 0, given by: For six activated mechanisms, the plastic flow rule is written The updated stress tensor in the return mapping algorithm is generally given by equation Equation (24). For an isotropic model, and taking into account the elastic rigidity matrix linking the principal stresses and strains given to Equation (5), Mohr-Coulomb's constitutive law written in principal stresses is ( [1,19]): For a single activated mechanism and taking into account its flow rule, the principal updated stresses give: with the accumulated plastic strain: By combining Equation (25a,c) we have: The consistency condition in this case is given by: By substituting Equation (27) in Equation (28) and after a few transformations, the analytical form of ∆γ is obtained by: In the case where the material properties are non-linear, we adopt a fixed-point algorithm to complete the resolution. The fixed point algorithm applied to the return mapping of the main plane mechanism (Algorithm 3) is based on the constitutive law of the activated mechanism (Equation (25)), with the initial conditions taken in the elastic domain , 0, ε p n = 0 . In this algorithm, the plastic multiplier ∆γ is obtained by substituting the principal trial stresses σ trial 1 , σ trial 2 , σ trial 3 in the consistency equation for this activated mechanism Equation (28). Moreover, from Equation (28) we deduce σ 1 .
Similarly, for two mechanisms activated on the right edge, the principal updated stresses are: with the accumulated plastic strain given by: The combination of equations Equation (30a,c) gives: By substitution of Equation (32) into Equation (28) and after some transformations, we obtain: For the second mechanism, the consistency condition is given by: The combination of equations Equation (30a,b) gives: By substitution of Equation (36) into Equation (35) and after some transformations, we obtain: From Equations (33) and (37), we can deduce the scalar form of ∆γ a1 and ∆γ b1 , defined by: with Regardless of the «physical» values taken by parameters K, G, φ and ψ, Det 1 = 0. (see Figure 3).  The fixed-point algorithm applied to the return mapping of the right edge mechanism (Algorithm 4) is based on the constitutive law given in Equation (30). Its initial conditions are also taken in the elastic domain ( 1 , 2 , 3 , 0,0, ̅ = 0).
For two mechanisms activated on the left edge, principal updated stresses are given by: Similarly, the accumulated plastic strain is given by: The combinations of equations Equation (39a,c) gives: By substituting Equation (41) into Equation (28) and making some transformations, we obtain: For the second activated mechanism, the consistency condition is given by: The combination of equations Equation (39a,b) gives: By substituting Equation (45) into Equation (44) and making some transformations, we obtain: From Equations (42) and (46), the scalar forms of ∆γ a2 and ∆γ b2 are: with
Finally, in the case of the Apex, the stress is located on the hydrostatic axis. The hydrostatic pressure is then written: where ∆ε p v represents the volumetric plastic strain increment and is given by: The updated stress and the accumulated plastic strain are given by [1]: with β ≡ cos φ sin ψ . The fixed-point algorithm applied to the return mapping to the Apex mechanism (Algorithm 6) is also based on the constitutive law of the activated mechanism (Equation (48)). Its initial conditions, taken in the elastic domain, are p trial n+1 , 0 and ε p n = 0. The volumetric plastic strain increment ∆ε p v 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: where σ n+1 represents the stress tensor at time t n+1 .

1D Illustration
Using the consistency condition . f = 0 , the expression of . σ is given by: .
and so, the relationship of the constitutive law linking . σ to . ε becomes: .
We can therefore deduce the tangent modulus D el p obtained analytically for the 1D model: 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 (k) n+1 is: 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 dε can be used: In Equation (57), ω is an appropriate small number (e.g., 10 −6 ) that represents the amount of perturbation. The perturbed strain ε per is: (58) and the strain increment is ∆ε For a given ε n+1 is obtained by the return mapping method presented in Algorithm 2. The perturbed stress σ (k) n+1 per is calculated from the perturbed strain ε (k) n+1 per using the same procedure. Finally, the tangent modulus is determined by the relationship: The illustration of this numerical method is presented in Algorithm 7.

Algorithm 7:
Numerical computation of the elastoplastic tangent modulus.
1. ε p n , α n and q n are known in step n and iteration k of Newton-Raphson 2.
Compute elastic trial stress σ trial n+1 and test for plastic loading 5.

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: Each component of the strain vector is thus perturbed according to the same principle stated in Section 2.3.1, and so allows an approximation of the elastoplastic tangent matrix to be calculated for each activated mechanism.

Numerical Study
In this section, we show the efficiency of the proposed numerical methods. For this purpose, two examples are considered. The first example concerns the solution of a bar (under a traction load) problem proposed by Kim [14]. This is a classical benchmark test for which there exists an analytical solution. Therefore, any numerical algorithm in plasticity must pass it.
The second example deals with a symmetrical rockfill dam modeled in plane strains according to Mohr-Coulomb's law.

Example 1: Nonlinear Bar Problem
The bar considered is uniform with a cross-sectional area A = 1 × 10 −4 m and a unit length L = 1 m. A traction force of magnitude F = 50 kN is applied at the end of the bar ( Figure 5).
to be calculated for each activated mechanism.

Numerical Study
In this section, we show the efficiency of the proposed numerical methods. For this purpose, two examples are considered. The first example concerns the solution of a bar (under a traction load) problem proposed by Kim [14]. This is a classical benchmark test for which there exists an analytical solution. Therefore, any numerical algorithm in plasticity must pass it.
The second example deals with a symmetrical rockfill dam modeled in plane strains according to Mohr-Coulomb's law.

Example 1: Nonlinear Bar Problem
The bar considered is uniform with a cross-sectional area = 1 × 10 −4 and a unit length = 1 . A traction force of magnitude = 50 is applied at the end of the bar ( Figure 5). The yield strength is = 400 . In the elastic region, the bar's Young's modulus is = 200 and in the plastic region, its plastic modulus is zero. Its yield function obeys Equation (7). The yield strength is σ Y = 400 MPa. In the elastic region, the bar's Young's modulus is E = 200 GPa and in the plastic region, its plastic modulus K is zero. Its yield function obeys Equation (7). Figure 6 shows the stress-strain for a single element of the 1D bar. Ten equal load levels are applied to the bar. In the legend, Ana-NR represents the results from the analytical calculation of D el p , Num-NR represents the numerical calculation of D el p , and (Kim) represents the exact solution [14].
Appl. Sci. 2021, 11, x FOR PEER REVIEW 18 of 29 Figure 6 shows the stress-strain for a single element of the 1D bar. Ten equal load levels are applied to the bar. In the legend, Ana-NR represents the results from the analytical calculation of , Num-NR represents the numerical calculation of , and (Kim) represents the exact solution [14]. The results obtained for both models and for the exact solution (Ana-NR, Num-NR and Kim) are virtually identical ( Figure 6). For loads below 40 kN, all models converge (Table 1)  the Euclidean norm) after at most two iterations. All the results show that the proposed numerical model for the calculation of the tangent module (Num-NR) is sufficiently accurate and efficient. In summary, the results prove that the exact solution is obtained no matter how many elements are used. The case of 10 elements is to show that the code behaves correctly for any mesh.  The results obtained for both models and for the exact solution (Ana-NR, Num-NR and Kim) are virtually identical ( Figure 6). For loads below 40 kN, all models converge (Table 1) in a single iteration (elastic domain). Above 40 kN (the plastic domain), the Ana-NR and Num-NR models converge (i.e., meet the convergence criterion which is F ext −F int F ext ≤ 10 −6 where F ext represents the external force, F int the internal force and · the Euclidean norm) after at most two iterations. All the results show that the proposed numerical model for the calculation of the tangent module (Num-NR) is sufficiently accurate and efficient. In summary, the results prove that the exact solution is obtained no matter how many elements are used. The case of 10 elements is to show that the code behaves correctly for any mesh.  Figure 6 also shows that the yield strength is reached at 400 MPa, and beyond this value, the bar becomes in plastic mode. To better observe these areas, we divided the bar into ten identical elements of two nodes each. Ten successive loads in 5 kN steps were applied. The evolution of plasticity in the bar during these loads for the Num-NR model is presented in Figure 7.

Example 2: Rockfill Dam
To test the efficiency of the algorithm in the case of a more complex behavior, we consider an example of a rockfill dam. This rockfill dam is assumed to be dry, homogeneous and symmetrical, and subject to its own weight. This dam has a width L of 320 m, a height h equal to 100 m and a crest length of 9 m. The depth of the dam axis is assumed to be infinite, which makes it possible to carry out a study of in-plane strain. Vertical and horizontal movements are zero at the base. Hooke's law is used in cases where the material has an elastic behavior. In the plastic case, the Mohr-Coulomb model is used in its nonassociated version to avoid excessive expansion of the material [9]. The material properties correspond to those used in [21] of a rockfill dam: Young's modulus E = 17 × 10 7 Pa , Poisson coefficient ν = 0.33, density ρ = 2370 Kg/m 3 , material cohesion c = 0, internal friction angle φ = 45 • and material dilatancy angle ψ = 15 • . We propose to compare the results obtained with our code developed in the Matlab environment and those obtained using the commercial software PLAXIS for the mesh sizes presented in Figure 8. The construction of the dam is simulated by adding successive layers, with 20 layers in total. The vertical and horizontal displacements, as well as the principal stresses in the first and second direction for the entire domain, are presented in Figure 9. Although the comparison of isovalues for the Matlab and PLAXIS models must be done with care, the displacements and calculated main stresses show a satisfactory correspondence. Indeed, for displacements, we note a relative difference of only 2% and 0.1% respectively for the maximum values of u x and u y . For the maximum stresses, there is a relative difference of 4% and 11% for the principal stress 1 and principal stress 2, respectively.

Example 2: Rockfill Dam
To test the efficiency of the algorithm in the case of a more complex behavior, we consider an example of a rockfill dam. This rockfill dam is assumed to be dry, homogeneous and symmetrical, and subject to its own weight. This dam has a width L of 320 m, a height h equal to 100 m and a crest length of 9 m. The depth of the dam axis is assumed to be infinite, which makes it possible to carry out a study of in-plane strain. Vertical and horizontal movements are zero at the base. Hooke's law is used in cases where the material has an elastic behavior. In the plastic case, the Mohr-Coulomb model is used in its non-associated version to avoid excessive expansion of the material [9]. The material properties correspond to those used in [21] of a rockfill dam: Young's modulus = 17 × The vertical and horizontal displacements, as well as the principal stresses in the first and second direction for the entire domain, are presented in Figure 9. Although the comparison of isovalues for the Matlab and PLAXIS models must be done with care, the displacements and calculated main stresses show a satisfactory correspondence. Indeed, for displacements, we note a relative difference of only 2% and 0.1% respectively for the maximum values of ux and uy. For the maximum stresses, there is a relative difference of 4% and 11% for the principal stress 1 and principal stress 2, respectively. Comparisons of the values obtained for this approach, particularly in the column and row under consideration (Figure 10), show that the vertical (u y ) and horizontal (u x ) displacements, as well as the normal stresses in the X and Y directions, are also in agreement with those obtained with PLAXIS ( Figures 11 and 12, respectively).
Appl. Sci. 2021, 11, x FOR PEER REVIEW 22 of 29 Comparisons of the values obtained for this approach, particularly in the column and row under consideration (Figure 10), show that the vertical (uy) and horizontal (ux) displacements, as well as the normal stresses in the X and Y directions, are also in agreement with those obtained with PLAXIS ( Figures 11 and 12, respectively).    Comparisons of the values obtained for this approach, particularly in the column and row under consideration (Figure 10), show that the vertical (uy) and horizontal (ux) displacements, as well as the normal stresses in the X and Y directions, are also in agreement with those obtained with PLAXIS ( Figures 11 and 12, respectively).   However, some small differences are observed. On the one hand, these could be attributed to the mesh sizes and to the shape of the elements, as those developed for our Matlab code differ from those used for the commercial software PLAXIS. On the other hand, these differences may also be due to how the loading was applied. It should be noted, however, that the results are very encouraging, especially since we do not need to calculate the tangent matrix analytically with our numerical model for more complex constitutive laws.
The code developed for Matlab also offers the possibility of studying the problem layer by layer according to the imposed load, and by refining the mesh size. These loads are represented by the different layers of the domain. Table 2 shows the number of iterations and the convergence calculated according to the loading for the numerical model (i.e., the elastoplastic tangent matrix is calculated by the finite difference method). For the three meshes (20 × 10, 60 × 20 and 30 × 30), the convergence tolerance in the equilibrium iterations has been set to 10 −10 . The maximum number of iterations in any load step (layer) is 8. However, some small differences are observed. On the one hand, these could be attributed to the mesh sizes and to the shape of the elements, as those developed for our Matlab code differ from those used for the commercial software PLAXIS. On the other hand, these differences may also be due to how the loading was applied. It should be noted, however, that the results are very encouraging, especially since we do not need to calculate the tangent matrix analytically with our numerical model for more complex constitutive laws.
The code developed for Matlab also offers the possibility of studying the problem layer by layer according to the imposed load, and by refining the mesh size. These loads are represented by the different layers of the domain. Table 2 shows the number of iterations and the convergence calculated according to the loading for the numerical model (i.e., the elastoplastic tangent matrix is calculated by the finite difference method). For the three meshes (20 × 10, 60 × 20 and 30 × 30), the convergence tolerance in the equilibrium iterations has been set to 10 −10 . The maximum number of iterations in any load step (layer) is 8.

Conclusions
Two methods for solving elastoplastic problems have been developed in this work. The first is related to the return-mapping scheme and is based on the fixed-point algorithm. The second is the calculation of the consistent tangent matrix and is based on finite differences. These two methods were implemented on an in-house finite element Matlab code. The algorithms were first validated by applying them to a 1D problem. The numerical results obtained during the simulations proved to be in excellent agreement with the exact solution [14]. They were then extended to the problem of a 2D dam subjected to plane strains and using Mohr-Coulomb's law. Two numerical methods for solving the system of non-linear return-mapping equations and the calculation of the tangent elastoplastic matrix were presented. Using these approaches, a set of simulations was performed, and the results compared to those obtained using the commercial software PLAXIS. The results of the different approaches remain in agreement with those given by PLAXIS, with deviations of around 2% and 0.1% in displacements for ux and uy, and of a maximum of 11% in stresses. These results are very encouraging, especially since we do not need to analytically calculate the often-tedious tangent matrix, thereby confirming the effectiveness of the numerical models developed in this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
For the analytical method, this matrix will first be expressed in the principal basis, as the yield functions have already been expressed. Then, by means of the chain rule derivation, it will be returned to the global basis using the DGISO2 procedure proposed by de Souza Neto et al. [1] for calculating the tensor derivative of a 2D isotropic function.
Knowing that dσ tr j = Ddε j (A2) and because the trial stress is assumed to be elastic, Equation (A2) becomes: representing the normal vector applied to the flow potential of the main plane, and the differentiation of the consistency condition given in Equation (28) is defined by d f tr 1 σ tr , c = 2V cm dε j , j = 1, 2, 3 (A5) representing the normal vector applied to the yield function of the main plane, so that dσ j = D − 4 a V pm ·V T cm dε j , j = 1, 2, 3 Finally, the consistent tangent operator in the principal basis and for the main plane D ep p,m obtained analytically can be written as: where V pr represents the normal vector applied to the flow potential of the right edge, and In this last relationship, we still have to express the incremental plastic multiplier as a function of the incremental plastic strain. Therefore, after differentiating the consistency conditions given in Equations (28) the increment of the j principal stress is then written as: The consistent tangent operator in the main base for the right edge D ep p,r obtained analytically is therefore: Remark A1. It should be noted that the analytical calculation of the elastoplastic tangent tensor in this case is only valid for constant parameters (without any material non-linearity as in the Duncan-Chang hyperbolic constitutive law [22] or the HS law [6]).