Elastoplastic Integration Method of Mohr-Coulomb Criterion

: A new method for implicit integration of the Mohr-Coulomb non-smooth multisurface plasticity models is presented, and Koiter’s requirements are incorporated exactly within the proposed algorithm. Algorithmic and numerical complexities are identiﬁed and introduced by the nonsmooth intersections of the Mohr-Coulomb surfaces; then, a projection contraction algorithm is applied to solve the classical Kuhn–Tucker complementary equations which provide the only characterization of possible active yield surfaces as a special class of variational inequalities, and the actual active yield surface is further determined by iteration. The basic idea is to calculate derivatives of the yield and potential functions with the expressions in the principal stresses and perform the return manipulations in the general stress space. Based on the principal stress characteristic equation, partial derivatives of principal stresses are calculated. The proposed algorithm eliminates the error caused by smoothing the corner of Mohr-Coulomb surfaces, avoids the numerical singularity at the intersections in the general stress space, and does not require the stress transformation needed in the principal stress space method. Lastly, several numerical examples are given to verify the validity of the proposed method.


Introduction
In geotechnical engineering, plastic models possessing non-smooth multiple yield surfaces are widely used in engineering applications, among which the Mohr-Coulomb criterion is the most typical. As is well known, the Mohr-Coulomb yielding surface is an irregular hexagonal cone in the principal stress space, which is composed of six facets with six edges and one apex. When the orthogonal law is applied at a point on an edge of the yielding surface, the normal is not unique and the final active yield surface is unknown, which introduces severe algorithmic and numerical complexities.
Huge efforts have been paid to the constitutive integration of plasticity with Mohr-Coulomb criterion. By applying the loading criterion for each yield surface of the singular point separately, Koiter proposed the evolution of plastic strain for multiple yield surfaces [1], which sets the foundation for the calculation of non-smooth multisurface plasticity. In order to overcome the serious numerical singularity at the corner of the Mohr-Coulomb criterion, the methods mentioned in the existing literature can be divided into three categories: (1) substitution method, Zienkiewicz et al. suggested that the suitable smooth yield criterion should be used to replace the yield criterion with singularity [2], that is, the Drucker-Prager criterion should be used to replace the Mohr-Coulomb criterion; near the corner of Mohr-Coulomb criterion, Owen and Hinton used the Drucker-Prager criterion to calculate the increment of plastic strain [3], so as to cut the corner of the Mohr-Coulomb yield surface. Marque [4] adopted a similar processing method, but such substitution would result in inaccurate numerical calculation results and affect the convergence speed of the overall equilibrium equation. (2) corner rounding method: Zienkiewicz and Pande [5] modified the Mohr-Coulomb yield surfaces and obtained a new yield criterion in the form of a curve, which eliminates the corners of Mohr-Coulomb on the meridional plane and π plane; Sloan and Booker [6] proposed a similar modified Mohr-Coulomb criterion for rounding corners, which is continuous and differentiable for all stress components; Abbo and Sloan [7] also smoothed the corners of Mohr-Coulomb, and the approximate yield function obtained is not only continuous but also second-order differentiable. Although these methods smoothed the Mohr-Coulomb yield surfaces approximately, the rapid change of derivative gradient near the corner point would lead to the instability of the numerical calculation. (3) principal stress space method: Larsson and Runesson [8] may be the first to perform the stress return in the principal stress space to avoid the singularity of Mohr-Coulomb functions expressed by invariants; then, Perić and de Souza Neto [9] extended this method to incorporate hardening; in the works of de Souza Neto et al. [10], the above method is further elaborated and extended to various common yield criteria and advanced yield criteria, such as modified Cam-clay and capped Drucker-Prager models.
Meanwhile, in order to handle the problem of unknown final active yield surfaces, the methods adopted in the literature can also be divided into three categories: (1) stress indicator method, De Borst proposed the concept of stress indicator factor to determine which group of yield surfaces are finally activated according to its sign [11]. Pankaj and Bićanić further developed this technology [12]. The generality of this method is limited because stress indication factors corresponding to different yield functions are different.
(2) geometric method, according to the geometric properties of Mohr-Coulomb yield surfaces in the principal stress space, Clausen et al. [13] derived the corresponding stress return geometry method, and this idea was later extended to other plastic models [14].
(3) Kuhn-Tucker method, Kuhn-Tucker complementarity theory has been widely used in an elastic-plastic constitutive integral [15,16]. Simo et al. adopted the mathematical idea of convex programming and, combined with the nearest point projection algorithm, the final active surface can be determined by iteration [17]. This method was continuously improved and applied to more complex elastic-plastic models [18][19][20]. The last method is the most widely used method at present because it provides a general and robust framework.
Based on Koiter's rule and the Kuhn-Tucker complementarity condition, a constitutive integral method suitable for the Mohr-Coulomb criterion is proposed in this paper. Koiter's rule determines the evolution of plastic strain and Kuhn-Tucker conditions characterize the possible active yield surfaces. The Kuhn-Tucker complementary equation is treated as a special class of variational inequalities and solved by a projection contraction algorithm.
Partial derivatives of principal stresses are calculated on the basis of the characteristic equation; then, the stress is returned in the general stress space, which not only avoids the singularity at the corner, but there is also no need for spectral decomposition and coordinate transformation. The algorithm can accommodate various widely used non-smooth multisurface plasticity models; single surface plasticity is also included. To avoid unnecessary algorithmic complexity, we restrict ourselves to the case of perfect plasticity. The generalization to isotropic hardening (or softening) can be made without any conceptual difficulties.

The Causes of Corner Problem of Mohr-Coulomb Criterion
The classical Mohr-Coulomb criterion can be expressed in the following form: where τ represents the shear stress; σ n represents the normal stress; c is cohesion; φ is the internal friction angle. The geometric interpretation of Equation (1) is shown in Figure 1, and the Mohr-Coulomb yield function can also be expressed as follows: which, rearranged, gives  Letting σ 1 ≥ σ 2 ≥ σ 3 , Equation (3) can be written as Complete yield surface can be obtained after considering all possible stress combinations. Generally, the Mohr-Coulomb yield function expressed by the invariant instead of the principal stress is used in calculation, namely in which According to Koiter's law, when the Mohr-Coulomb criterion is used to calculate plastic strain, we have Here, γ α represents the plastic multiplier corresponding to the yield surface α. A partial derivative of the Mohr-Coulomb yield function to stress component can be calculated based on the chain derivation rule: in which As we can see, when θ σ = ± π 6 (the corner of Mohr-Coulomb yield surface), tan 3θ σ and 1 cos 3θ σ in Equations (13) and (14) are undefined, that is, singularity arises from undefined partial derivatives. On the other hand, when the Mohr-Coulomb criterion is used for calculation, the final active yield surface at its corner is unknown. In the incremental form of the elastoplastic constitutive integral, after the strain increment is applied, whether the stress point is in elastic state or plastic loading state is judged in the form of test stress based on the following conditions: f trial α,n+1 ≤ 0 for all α ∈ (1, 2, · · · , 6) ⇒ elastic state f trial α,n+1 > 0 for some α ∈ (1, 2, · · · , 6) ⇒ plastic state (15) where f trial α,n+1 represents the test state of yield surface α corresponding to the strain increment, and the specific calculation method is introduced later.
It should be noted that, if only one yield surface is activated, namely, only one f trial α,n+1 > 0, then the yield surface is the final active yield surface, that is, f α,n+1 > 0. However, when multiple yield surfaces are activated, f trial α,n+1 > 0 cannot guarantee that the yield surface is the final active yield surface, namely, there may be f trial α,n+1 > 0 but f α,n+1 < 0. As we can see in Figure 2, N 1 and N 2 are the normal direction of yield surface f 1 and f 2 , C 1, C 2 , and C 12 form the corner cone of f 1 and f 2 together. When σ trial n+1 ∈ C 12 , we have f trial 1,n+1 > 0 and f trial 2,n+1 > 0; at the same time, the plastic multiplier ∆γ α > 0,α = 1, 2; therefore, the activated yield surfaces f 1 and f 2 both are the final active yield surfaces; when σ trial n+1 ∈ C 1 or σ trial n+1 ∈ C 2 , although f trial 1,n+1 > 0 and f trial 2,n+1 > 0 can be obtained, it can be seen that ∆γ 2 < 0 in C 1 and ∆γ 1 < 0 in C 2 , that is, although the activated yield surfaces are two, there is only one final active surface in the end.

Constitutive Integral for Mohr-Coulomb Criterion
The elastoplastic constitutive equation can be written as the following form: Here,  ,  , and p  represent stress, strain, and plastic strain, respectively, and D is the elastic matrix. The plastic strain p  is determined according to Koiter's rule in Equation (11), and the plastic multiplier satisfies the Kuhn-Tucker complementarity equation: For the Mohr-Coulomb model with multiple yield surfaces, it is necessary to define the elastic space, which can be represented by E and defined as:

Constitutive Integral for Mohr-Coulomb Criterion
The elastoplastic constitutive equation can be written as the following form: Here, σ, ε, and ε p represent stress, strain, and plastic strain, respectively, and D is the elastic matrix. The plastic strain ε p is determined according to Koiter's rule in Equation (11), and the plastic multiplier satisfies the Kuhn-Tucker complementarity equation: For the Mohr-Coulomb model with multiple yield surfaces, it is necessary to define the elastic space, which can be represented by E σ and defined as: Furthermore, Equation (18) can be analytically written in terms of the principal stresses σ 1 , σ 2 , σ 3 as: Considering a time discretization of interest, and letting σ n , ε n , ε p n be the initial stress and strain, then, given an incremental strain ∆ε and applying an implicit backward Euler difference to Equations (11) and (16), we have Accordingly, the discrete form of Equation (17) is expressed as: Finally, by setting ∆γ α = 0 in Equations (20) and (21), the trial state is obtained formally: A yield surface f α,n+1 is termed active if the condition f α,n+1 > 0 holds; once the trial state is obtained, an initial set of trial constraints of Mohr-Coulomb model is defined as: If J trial act = ∅, it indicates that it is currently in an elastic state, and the test state is the real stress state. Then, assuming plastic loading, it holds that J trial act = ∅; the trial stress state σ trial n+1 lays beyond the elastic region E σ and is considered inadmissible. The solution σ n+1 to the discrete problem can be expressed as: Substituting (24) into Kuhn-Tucker conditions (21), we have Equation (25)  Compute elastic predictor Check for plastic process Evaluate stress at iteration (k) act , Exit Else continue 5 Computate plastic multiplier In elastoplastic computations using the finite element method, the load is applied incrementally. By performing the above iteration, the stress state after applying the strain increment is obtained.

Partial Derivatives of Principal Stresses with Regard to Stress Components
For Mohr-Coulomb surfaces, which have sharp corners and can be explicitly expressed in terms of principal stresses,σ i , i = 1, 2, and 3. While calculating plastic strains, knowing partial derivatives of principal stresses with regard to stress components will give us great convenience.
Principal stress σ i satisfies the following characteristic equation: where the coefficients I 1 , I 2 , and I 3 are functions of the six stress components σ x , σ y , . . . , τ xy , defined as Differentiating Equation (26) with regard to one of σ x , σ y , . . . , τ xy , we have with the single quote associated with a quantity denoting the derivative of the quantity with regard to one of σ x , σ y , . . . , τ xy , or, equivalently, Thus, we have Similarly, we can obtain the second derivative of the principal stress with regard to σ x , σ y , . . . , τ xy .

Algorithm PCA
Given a subset K of Euclidean n-dimensional space R n and a mapping F : K → R n , the variational inequality finds a vector x ∈ K such that Given a mapping F : R n + → R n , the nonlinear complementarity problem is to find a vector x ∈ R n satisfying 0 ≤ x⊥F(x) ≥ 0 In fact, when the subset K in the variational inequality is non-negative, the nonlinear complementarity problem is equivalent to the variational inequality [21]. As a method for solving convex optimization problems under the framework of variational inequalities, just as the name implies, the projection contraction algorithm is a contraction algorithm based on projection. Here, we only refer to the algorithm part; for more details, see [22]. Rewrite complementary equation system (21) to vector form: in which ∆γ = (∆γ 1 , . . . , ∆γ m ), F(∆γ) = (− f 1 (∆γ 1 ), . . . , − f m (∆γ m )). As we can see from Equation (35), the complementary system in (36) can be summarized as a typical nonlinear complementarity problem, which can be solved by the projection contraction method. The projection contraction algorithm (PCA) is invoked in this way: The input arguments of ∆γ 0 is the initial guess of plastic multiplier ∆γ. In general, ∆γ 0 = 0.
The pseudocode of PCA is listed as follows, and He proposed [22] three parameters in PCA as ξ = 0.9, η = 0.4, and λ = 1.9, which is also adopted in this paper. It should be noted that the parameters have little effect on the result.

Numerical Examples
The practical application of the proposed algorithm is illustrated in this section by two different numerical examples. For all the examples, four node isoparametric elements with four Gaussian points are adopted, and the error tolerance for unbalanced force in the equilibrium iteration is 0.1%.

Strip Footing Collapse
The bearing capacity of a strip footing is calculated under undrained conditions by the proposed method. The soil is assumed weightless, isotropic, and is modeled as a Mohr-Coulomb perfectly plastic material. There is no friction at the footing/soil interface. The computational model and discretized mesh are shown in Figure 3. The material parameters are assumed as follows: Young's modulus: E = 107 MPa, Poisson's ratio: v = 0.48, cohesion c = 490 kPa, and internal friction angle φ = 20 • . A plane strain state is assumed.
The strip footing is subjected to a vertical pressure P and the aim of present analysis is to obtain the limit pressure P lim . A solution to this problem has been derived by Prandtl and Hill [23]: P lim = c e π tan ϕ tan 2 (45 The pressure P is increased incrementally to failure and corresponding load/settlement (average) curves obtained are plotted in Figure 4.
Mohr-Coulomb perfectly plastic material. There is no friction at the footing/soil interface. The computational model and discretized mesh are shown in Figure 3 The strip footing is subjected to a vertical pressure P and the aim of present analysis is to obtain the limit pressure Plim. A solution to this problem has been derived by Prandtl and Hill [23]: The pressure P is increased incrementally to failure and corresponding load/settlement (average) curves obtained are plotted in Figure 4.   Figure 5. A from Figure 5, at the moment that the foundation enters flow, the global sh the foundation develops downward to a certain depth and extends to the failure mechanism captured in the finite element analysis is in good agreem The equivalent plastic strain at the failure load is shown in Figure 5. As we can see from Figure 5, at the moment that the foundation enters flow, the global shear failure of the foundation develops downward to a certain depth and extends to the ground The failure mechanism captured in the finite element analysis is in good agreement with the slip-line field of strip footing as shown in Figure 6 [24]. the foundation develops downward to a certain dept failure mechanism captured in the finite element analy slip-line field of strip footing as shown in Figure 6 [24]

Slope Stability
In this example, the plane strain analysis of a slope subjected to self-weigh out. To assess the safety of the slope shown in Figure 7, the strength reduction adopted. Table 1 lists the mechanical parameters. The materials comply with th Mohr-Coulomb criterion. Horizontal displacements are fixed for nodes along th

Slope Stability
In this example, the plane strain analysis of a slope subjected to self-weight is carried out. To assess the safety of the slope shown in Figure 7, the strength reduction method is adopted. Table 1 lists the mechanical parameters. The materials comply with the perfectly Mohr-Coulomb criterion. Horizontal displacements are fixed for nodes along the left and right boundaries while both horizontal and vertical displacements are fixed along the bottom boundary.
adopted. Table 1 lists the mechanical parameters. The materials comply with the pe Mohr-Coulomb criterion. Horizontal displacements are fixed for nodes along the le right boundaries while both horizontal and vertical displacements are fixed along th tom boundary.   For slopes, the factor of safety Fs is traditionally defined as the ratio of the actual shear strength to the minimum shear strength required to prevent failure. An obvious way of computing Fs with a finite element or a finite difference program is simply to reduce the shear strength until collapse occurs. This technique was used as early as 1975 by Zienkiewicz et al. [25], and has since been applied by Naylor [26], Matsui and San [27], Ugai and Leschinsky [28], Dawson et al. [29], and many others. There are a variety of criteria of the limit equilibrium state, see [30]. The differences of the factors of safety based on different criteria are small.
The initial reduction factor is set as 1, then increases gradually until the slope collapse. It is noteworthy that Young's modulus E and Poission's ratio, v are rectified based on [31] while using the strength reduction technique.
The horizontal displacement at the top corner of the slope with the reduction factor is shown in Figure 8; when the reduction factor F = 1.35, there is an obvious inflection point on the displacement curve. At this time, the distribution of the plastic zone is shown in Figure 9, the plastic zone is just go-through.
The horizontal displacement at the top corner of the slope with the red is shown in Figure 8; when the reduction factor F = 1.35, there is an obvio point on the displacement curve. At this time, the distribution of the plastic z in Figure 9, the plastic zone is just go-through. If the reduction factor at the inflection point is taken as the safety fact 1.35; if the reduction factor when the calculation fails to converge is taken factor, then Fs = 1.38. The gap between the two results is very small. Figure 9 once an earth slope is led to the limit equilibrium state by means of the streng technique, a plastic zone will go through the slope from the toe to the top.
The slope safety factor obtained by Spencer is Fs = 1.33, and the correspo slip surface is shown in Figure 10, in which green line represents the slip surf colors top-down represents different soils in Table 1; due to the obeserva critical slip surface within the plastic zone will consist of the points at which th plastic strain takes the local maximum in the vertical direction [32], we can  Geotechnics 2022, 2, FOR PEER REVIEW obtained critical sliding surface by the Spencer method is in good agreement w failure mechanism shown in Figure 9.  If the reduction factor at the inflection point is taken as the safety factor, then Fs = 1.35; if the reduction factor when the calculation fails to converge is taken as the safety factor, then Fs = 1.38. The gap between the two results is very small. Figure 9 shows that, once an earth slope is led to the limit equilibrium state by means of the strength reduction technique, a plastic zone will go through the slope from the toe to the top.
The slope safety factor obtained by Spencer is Fs = 1.33, and the corresponding critical slip surface is shown in Figure 10, in which green line represents the slip surface and three colors top-down represents different soils in Table 1; due to the obeservation that the critical slip surface within the plastic zone will consist of the points at which the equivalent plastic strain takes the local maximum in the vertical direction [32], we can see that the obtained critical sliding surface by the Spencer method is in good agreement with the failure mechanism shown in Figure 9.  A comparison about computational efficiency is made between the PC sic mapping return method in [17]. Integral points that are undergoin mation are selected. A comparison is made based on the computation stress returns, and is implemented in Matlab on a computer with Intel(R 3230M 2.60 GHz processor. Table 2 lists the results. It can be seen that the present PCA is substantially more efficient mapping return method, especially when more than one yield surface is e

Conclusions
Based on the complementarity theory, the elasto-plastic complementa established. The corner problems of multiple surfaces are unified into a set tary equations by using Koiter's law, and then solved by using the projec algorithm. The theory and examples show that: the projection contractio the complementary equation with multiple yield surfaces avoids the acc vergence problems caused by corner smoothing; the stress return operatio in the general stress space, and the partial derivative of the correspondin A comparison about computational efficiency is made between the PCA and the classic mapping return method in [17]. Integral points that are undergoing plastic deformation are selected. A comparison is made based on the computation time for 10,000 stress returns, and is implemented in Matlab on a computer with Intel(R) Core(TM) i5-3230M 2.60 GHz processor. Table 2 lists the results. It can be seen that the present PCA is substantially more efficient than the classic mapping return method, especially when more than one yield surface is ever activated.

Conclusions
Based on the complementarity theory, the elasto-plastic complementarity equation is established. The corner problems of multiple surfaces are unified into a set of complementary equations by using Koiter's law, and then solved by using the projection contraction algorithm. The theory and examples show that: the projection contraction algorithm for the complementary equation with multiple yield surfaces avoids the accuracy and convergence problems caused by corner smoothing; the stress return operation is carried out in the general stress space, and the partial derivative of the corresponding stress component is calculated in the principal stress space, which not only avoids the numerical singularity at the corner, but also eliminates the stress transformation of the stress return in the principal stress space.

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