Next Article in Journal
Expressing Robot Personality through Talking Body Language
Next Article in Special Issue
Numerical Simulation of Flame Retardant Polymers Using a Combined Eulerian–Lagrangian Finite Element Formulation
Previous Article in Journal
Sargassum Influx on the Mexican Coast: A Source for Synthesizing Silver Nanoparticles with Catalytic and Antibacterial Properties
Previous Article in Special Issue
Consistent Implicit Time Integration for Viscoplastic Modelingof Subsidence above Hydrocarbon Reservoirs
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Gildas Yaovi Amouzou
and
Azzeddine Soulaïmani
*
Departement of Mechanical Engineering, École de Technologie Supérieure, 1100 Notre-Dame W., Montréal, QC H3C 1K3, Canada
*
Author to whom correspondence should be addressed.
Appl. Sci. 2021, 11(10), 4637; https://doi.org/10.3390/app11104637
Submission received: 19 April 2021 / Revised: 15 May 2021 / Accepted: 17 May 2021 / Published: 19 May 2021
(This article belongs to the Special Issue Element-Based Methods for the Solution of Engineering Problems)

Abstract

:
Two numerical algorithms for solving elastoplastic problems with the finite element method are presented. The first deals with the implementation of the return mapping algorithm and is based on a fixed-point algorithm. This method rewrites the system of elastoplasticity non-linear equations in a form adapted to the fixed-point method. The second algorithm relates to the computation of the elastoplastic consistent tangent matrix using a simple finite difference scheme. A first validation is performed on a nonlinear bar problem. The results obtained show that both numerical algorithms are very efficient and yield the exact solution. The proposed algorithms are applied to a two-dimensional rockfill dam loaded in plane strain. The elastoplastic tangent matrix is calculated by using the finite 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.

1. 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.

2. Plasticity Formulation

2.1. Constitutive Law of Plasticity

In linear elasticity, the constitutive law is given by:
σ = D : ε e
where in plane strain case, σ = [ σ x x   σ y y   σ x y ] T is the stress vector, ε e = [ ε x x e ε y y e ε x y e ] 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:
ε = ε e + ε p
and the stress-strain relationship is written as:
σ = D : ( ε ε p )
This relationship can also be written in terms of tensor rates (or in incremental form linking the increments strains to those of the stresses):
σ ˙ = D : ( ε ˙ ε ˙ p )
with
ε ˙ p = γ ˙ f σ
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]:
D = 2 G   I + ( K 2 3 G )   1 1  
where K ,   G > 0 represents the elastic and shear modulus, respectively, I = 1 2 [ δ i k δ j l + δ i l δ j k ] e i e j e k e l   is the fourth symmetric order identity tensor and 1 = δ i j e i e j is the second order of identity tensor.

2.2. 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]:
{ ε n + 1 e ε n + 1 e   t r i a l + Δ γ   N ( σ n + 1 , A n + 1 )       = 0 α n + 1 α n Δ γ   H ( σ n + 1 , A n + 1 )             = 0 f ( σ n + 1   , A n + 1     )                                         = 0
where at time n + 1, ε n + 1 e represents the elastic strain tensor, ε n + 1 e   t r i a l 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.

2.2.1. 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:
f ( σ , q ) = | σ q | ( σ Y + K α )
where σ Y   represents the yield strength, K the plastic modulus and α the internal hardening variable, which has the evolution equation:
α ˙ = γ ˙
and q is the back stress defined by Ziegler’s rule as:
q ˙ = γ ˙   H   f σ = γ ˙   H   sign ( σ q ) ,
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.
  • ε n p ,   α n   and   q n   are   known   in   step   n
  • Given   a   strain   ε n + 1 = ε n + Δ ε n
  • Compute   elastic   trial   stress   σ n + 1 t r i a l   and   test   for   plastic   loading  
  • σ n + 1 t r i a l = E ( ε n + 1 ε n p )
  • ξ n + 1 t r i a l = σ n + 1 t r i a l q n
  • f n + 1 t r i a l = | ξ n + 1 t r i a l | [ σ Y + K α n ]
  • if   f n + 1 trial 0   then
  • Elastic   step :   ( ) n + 1 = ( ) n + 1 trial
  • Else plastic   step ,   solve   the   system :  
    { σ n + 1 = σ n + 1 t r i a l Δ γ   E   sign ( ξ n + 1 )     ε n + 1 p = ε n p + Δ γ   sign ( ξ n + 1 ) α n + 1 = α n + Δ γ q n + 1 = q n + Δ γ   H   sign ( ξ n + 1 ) f n + 1 = | ξ n + 1 | [ σ Y + K α n + 1 ]
  • EndIf  
{ σ n + 1 = σ n + 1 t r i a l Δ γ   E   sign ( ξ n + 1 )     ε n + 1 p = ε n p + Δ γ   sign ( ξ n + 1 ) α n + 1 = α n + Δ γ q n + 1 = q n + Δ γ   H   sign ( ξ n + 1 ) f n + 1 = | ξ n + 1 | [ σ Y + K α n + 1 ]
The analytical resolution of the system Equation (11) gives [4]:
Δ γ = f n + 1 t r i a l E + [ K + H ] > 0 σ n + 1 = σ n + 1 t r i a l Δ γ   E   sign ( ξ n + 1 ) ε n + 1 p = ε n p + Δ γ   sign ( ξ n + 1 ) q n + 1 = q n + Δ γ   H   sign ( ξ n + 1 ) α n + 1 = α n + Δ γ
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 ,   ε n + 1 p ] 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 [ σ n + 1 t r i a l ,   α n , 0 ,   q n ,   ε n p ] . The iteration is stopped for a time step n + 1 N m a x 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 n + 1 trial 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.
  • Initialization :   Z 0 = [ σ n + 1 t r i a l ,   α n , 0 ,   q n ,   ε n p ]
  • t o l :   tolerance   and   N m a x   maximum   number   of   iterations   are   set
  • While   Z Z 0 / Z 0 > tol   and   n N max = 20
  • { Δ γ = ( σ n + 1 σ n + 1 t r i a l ) E sign ξ n + 1 ε n + 1 p = ε n p + Δ γ   sign ( ξ n + 1 ) α n + 1 = α n + Δ γ q n + 1 = q n + Δ γ   H   sign ( ξ n + 1 )
  • If   σ n + 1 q n + 1 0   then
  • σ n + 1 = q n + 1 + [ σ Y + K α n + 1 ]
  • else
  • σ n + 1 = q n + 1 [ σ Y + K α n + 1 ]
  • EndIf
  • Z = [ σ n + 1 ,   α n + 1 , Δ γ ,   q n + 1 ,   ε n + 1 p ]
  • Update   Z 0 = Z
  • EndWhile

2.2.2. 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 ε ¯ n p and its increment Δ ε ¯   p ( c ( ε ¯ n p , Δ ε ¯   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:
τ = c σ n   tan ( ϕ )
The plastic function f is expressed in terms of principal stresses by:
f ( σ , c ) = σ 1 σ 3 + ( σ 1 + σ 3 ) sin ϕ 2 c cos ϕ
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:
p = c cot ( ϕ )
The six edges of this criterion in three-dimensional space, allow it to be expressed by six load functions given by [1]:
f 1 ( σ , c ) = σ 1 σ 3 + ( σ 1 + σ 3 ) sin ϕ 2 c cos ϕ f 2 ( σ , c ) = σ 2 σ 3 + ( σ 2 + σ 3 ) sin ϕ 2 c cos ϕ f 3 ( σ , c ) = σ 2 σ 1 + ( σ 2 + σ 1 ) sin ϕ 2 c cos ϕ f 4 ( σ , c ) = σ 3 σ 1 + ( σ 3 + σ 1 ) sin ϕ 2 c cos ϕ f 5 ( σ , c ) = σ 3 σ 2 + ( σ 3 + σ 2 ) sin ϕ 2 c cos ϕ f 6 ( σ , c ) = σ 1 σ 2 + ( σ 1 + σ 2 ) sin ϕ 2 c cos ϕ
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:
g 1 ( σ , c ) = σ 1 σ 3 + ( σ 1 + σ 3 ) sin ψ 2 c cos ψ
The plastic flow rule is then written as:
ε ˙ p = i = 1 m γ i ˙ g i ( σ , c ) σ
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.
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:
ε ˙ p = γ   ˙ N a
with N a the flow vector normal to the corresponding plane defined by f 1 = 0 , given by:
N a = g 1 σ = ( 1 + sin ψ ) e 1 e 1 ( 1 sin ψ ) e 3 e 3
For two activated mechanisms, the plastic flow rule is written:
ε ˙ p = γ   ˙ a N a + γ   ˙ b N b
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:
N b = { N 6 = ( 1 + sin ψ ) e 1 e 1 ( 1 sin ψ ) e 2 e 2 N 2 = ( 1 + sin ψ ) e 2 e 2 ( 1 sin ψ ) e 3 e 3
For six activated mechanisms, the plastic flow rule is written
ε ˙ p = i = 1 6 γ i ˙ N i
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]):
σ j = σ j   t r i a l i = 1 6 Δ γ i ( 2 G [ N d i ] j + K   N v i ) ,     j = 1 , 2 , 3
where N v i tr [ N i ] represents the i-th volumetric component of the updated plastic flow vector, tr[ ] the trace of the second-order tensor and [ N d i ] j the j-th eigenvalue of its deviatoric projection.
From Equation (24), four return mappings of Mohr–Coulomb’s law according to the active mechanisms follow: to the smooth portion, to the right edge, to the left edge, and to the apex.
(a)
Return mapping to the smooth portion ( σ 1 > σ 2 > σ 3 )
For a single activated mechanism and taking into account its flow rule, the principal updated stresses give:
{ σ 1 = σ 1 t r i a l [ 2 G ( 1 + 1 3 sin ψ ) + 2 K sin ψ ] Δ γ                   ( a ) σ 2 = σ 2 t r i a l + ( 4 3 G 2 K ) Δ γ   sin ψ                       ( b )     σ 3 = σ 3 t r i a l + [ 2 G ( 1 1 3 sin ψ ) 2 K sin ψ ] Δ γ                       ( c )
with the accumulated plastic strain:
ε ¯ n + 1 p = ε ¯ n p + 2 Δ γ cos ϕ
By combining Equation (25a,c) we have:
{ σ 1 σ 3 = σ 1 t r i a l σ 3 t r i a l 4 G Δ γ σ 1 + σ 3 = σ 1 t r i a l + σ 3 t r i a l 4 ( G 3 + K ) Δ γ sin ψ
The consistency condition in this case is given by:
f 1   ( σ   , c ) = σ 1   σ 3   + ( σ 1   + σ 3   ) sin ϕ 2 c cos ϕ = 0
By substituting Equation (27) in Equation (28) and after a few transformations, the analytical form of Δ γ is obtained by:
Δ γ = f 1 t r i a l ( σ t r i a l , c ) a
with
a = 4 [ G + ( G 3 + K ) sin ψ sin ϕ ]
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 ( σ 1 t r i a l , σ 2 t r i a l , σ 3 t r i a l , 0 ,   ε ¯ n p = 0 ) . In this algorithm, the plastic multiplier ∆γ is obtained by substituting the principal trial stresses σ 1 trial , σ 2 trial , σ 3 trial in the consistency equation for this activated mechanism Equation (28). Moreover, from Equation (28) we deduce σ 1   .
Algorithm 3: Fixed point algorithm applied to the main plane mechanism.
  • Initialization :   Z 0 = [ σ 1 trial , σ 2 trial , σ 3 trial , 0 , ε ¯ n p ]
  •   tol :   tolerance   and   N max   maximum   number   of   iterations  
  • While   Z Z 0 Z 0 > tol   and   n N max = 100
  • a = 4 [ G + ( G 3 + K ) sin ψ sin ϕ ]  
  • Δ γ = ( σ 1 trial σ 3 trial ) + ( σ 1 trial + σ 3 trial ) sin ϕ 2   c ( ε ¯ n p , Δ ε ¯   p ) cos ϕ a
  • σ 1 = ( 1 sin ψ ) σ 3 + 2   c ( ε ¯ n p , Δ ε ¯   p )   cos ϕ   1 + sin ψ
  • σ 2 = σ 2 trial + Δ γ ( 4 G 3 2 K ) sin ψ
  • σ 3 = σ 3 trial + [ 2 G ( 1 1 3 sin ψ ) 2 K sin ψ ] Δ γ
  • ε ¯ n + 1 p = ε ¯ n p + 2 Δ γ cos ϕ
  • Update   c ,   K   and   G
  • Z = [ σ 1   , σ 2   , σ 3   , Δ γ , ε ¯ n + 1 p ]
  • Update Z 0 = Z
  • EndWhile
(b)
Return mapping to the right edge ( σ 1 > σ 2 = σ 3 )
Similarly, for two mechanisms activated on the right edge, the principal updated stresses are:
{ σ 1 = σ 1 t r i a l [ 2 G ( 1 + 1 3 sin ψ ) + 2 K sin ψ ] ( Δ γ a 1 + Δ γ b 1 )                                             ( a )   σ 2 = σ 2 t r i a l + ( 4 3 G 2 K ) Δ γ a 1   sin ψ +   [ 2 G ( 1 1 3 sin ψ ) 2 K sin ψ ] Δ γ b 1                     ( b ) σ 3 = σ 3 t r i a l + [ 2 G ( 1 1 3 sin ψ ) 2 K sin ψ ] Δ γ a 1 + ( 4 3 G 2 K ) Δ γ b 1   sin ψ               ( c )
with the accumulated plastic strain given by:
ε ¯ n + 1 p = ε ¯ n p + 2 ( Δ γ a 1 + Δ γ b 1 ) cos ϕ
The combination of equations Equation (30a,c) gives:
{ σ 1 σ 3 = σ 1 t r i a l σ 3 t r i a l 4 G Δ γ a 1 2 G ( 1 + sin ψ ) Δ γ b 1 σ 1 + σ 3 = σ 1 t r i a l + σ 3 t r i a l 4 ( G 3 + K ) Δ γ a 1 sin ψ   2 [ G ( 1 1 3 sin ψ ) + 2 K sin ψ ] Δ γ a 1  
By substitution of Equation (32) into Equation (28) and after some transformations, we obtain:
a Δ γ a 1 + b 1 Δ γ b 1 = f 1 t r i a l ( σ t r i a l , c )
with
b 1 = 2 G ( 1 + sin ψ + sin ϕ 1 3 sin ψ sin ϕ ) + 4 K sin ψ sin ϕ
For the second mechanism, the consistency condition is given by:
f 6   ( σ   , c ) = σ 1   σ 2   + ( σ 1   + σ 2   ) sin ϕ 2 c cos ϕ = 0
The combination of equations Equation (30a,b) gives:
{ σ 1 σ 2 = σ 1 t r i a l σ 2 t r i a l 4 G Δ γ b 1 2 G ( 1 + sin ψ ) Δ γ a 1 σ 1 + σ 2 = σ 1 t r i a l + σ 2 t r i a l 4 ( G 3 + K ) Δ γ b 1 sin ψ   2 [ G ( 1 1 3 sin ψ ) + 2 K sin ψ ] Δ γ a 1  
By substitution of Equation (36) into Equation (35) and after some transformations, we obtain:
a Δ γ b 1 + b 1 Δ γ a 1 = f 6 t r i a l ( σ t r i a l , c )
From Equations (33) and (37), we can deduce the scalar form of Δ γ a 1 and Δ γ b 1 , defined by:
{ Δ γ a 1 = a f 1 t r i a l b 1 f 6 t r i a l D e t 1 Δ γ b 1 = a f 6 t r i a l b 1 f 1 t r i a l D e t 1
or
{ Δ γ a 1 = f 1 t r i a l b 1 Δ γ b 1 a Δ γ b 1 = f 6 t r i a l b 1 Δ γ a 1 a
with
D e t 1 = a 2 b 1 2
Remark 1.
Regardless of the «physical» values taken by parameters K, G, ϕ and ψ, D e t 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 t r i a l , σ 2 t r i a l , σ 3 t r i a l , 0 , 0 ,   ε ¯ n p = 0 ).
Algorithm 4: Fixed point method applied to the right edge mechanism.
  • Initialization :   Z 0 = [ σ 1 trial , σ 2 trial , σ 3 trial , 0 ,   0 , ε ¯ n p ]
  • tol :   tolerance   and   N max   maximum   number   of   iterations  
  • While   Z Z 0 Z 0 > tol   and   n N max = 100
  • a = 4 [ G + ( G 3 + K ) sin ψ sin ϕ ]
  • b 1 = 2 G ( 1 + sin ψ + sin ϕ 1 3 sin ψ sin ϕ ) + 4 K sin ψ sin ϕ
  • Δ γ a 1 = ( σ 1 trial σ 3 trial ) + ( σ 1 trial + σ 3 trial ) sin ϕ 2   c ( ε ¯ n p , Δ ε ¯   p ) cos ϕ b 1 Δ γ b 1 a
  • Δ γ b 1 = ( σ 1 trial σ 2 trial ) + ( σ 1 trial + σ 2 trial ) sin ϕ 2   c ( ε ¯ n p , Δ ε ¯   p ) cos ϕ b 1 Δ γ a 1 a
  • σ 1 = σ 1 trial [ 2 G ( 1 + 1 3 sin ψ ) + 2 K sin ψ ] ( Δ γ a 1 + Δ γ b 1 )
  • σ 2 = σ 2 trial + ( 4 3 G 2 K ) Δ γ a 1   sin ψ + [ 2 G ( 1 1 3 sin ψ ) 2 K sin ψ ] Δ γ b 1
  • σ 3 = σ 3 trial + [ 2 G ( 1 1 3 sin ψ ) 2 K sin ψ ] Δ γ a 1 + ( 4 3 G 2 K ) Δ γ b 1   sin ψ
  • ε ¯ n + 1 p = ε ¯ n p + 2 ( Δ γ a 1 + Δ γ b 1 ) cos ϕ
  • Update   c ,   K   and   G
  • Z = [ σ 1   , σ 2   , σ 3   , Δ γ a 1 ,   Δ γ b 1 , ε ¯ n + 1 p ]
  • Update Z 0 = Z
  • EndWhile
(c)
Return mapping to the left edge ( σ 1 = σ 2 > σ 3 )
For two mechanisms activated on the left edge, principal updated stresses are given by:
{ σ 1 = σ 1 t r i a l [ 2 G ( 1 + 1 3 sin ψ ) + 2 K sin ψ ] Δ γ a 2 + ( 4 3 G 2 K ) Δ γ b 2   sin ψ             ( a ) σ 2 = σ 2 t r i a l + ( 4 3 G 2 K ) Δ γ a 2   sin ψ   [ 2 G ( 1 + 1 3 sin ψ ) + 2 K sin ψ ] Δ γ b 2           ( b ) σ 3 = σ 3 t r i a l + [ 2 G ( 1 1 3 sin ψ ) 2 K sin ψ ] ( Δ γ a 2 + Δ γ b 2 ) .                                         ( c )
Similarly, the accumulated plastic strain is given by:
ε ¯ n + 1 p = ε ¯ n p + 2 ( Δ γ a 2 + Δ γ b 2 ) cos ϕ .
The combinations of equations Equation (39a,c) gives:
{ σ 1 σ 3 = σ 1 t r i a l σ 3 t r i a l 4 G Δ γ a 2 2 G ( 1 sin ψ ) Δ γ b 2 σ 1 + σ 3 = σ 1 t r i a l + σ 3 t r i a l 4 ( G 3 + K ) Δ γ a 2 sin ψ + 2 [ G ( 1 + 1 3 sin ψ ) 2 K sin ψ ] Δ γ b 2 .
By substituting Equation (41) into Equation (28) and making some transformations, we obtain:
a Δ γ a 2 + b 2 Δ γ b 2 = f 1 t r i a l ( σ t r i a l , c )
with
b 2 = 2 G ( 1 sin ψ sin ϕ 1 3 sin ψ sin ϕ ) + 4 K sin ψ sin ϕ .
For the second activated mechanism, the consistency condition is given by:
f 2   ( σ   , c ) = σ 2   σ 3   + ( σ 2   + σ 3   ) sin ϕ 2 c cos ϕ = 0 .
The combination of equations Equation (39a,b) gives:
{ σ 2 σ 3 = σ 2 t r i a l σ 3 t r i a l 4 G Δ γ b 2 2 G ( 1 sin ψ ) Δ γ a 2 σ 2 + σ 3 = σ 2 t r i a l + σ 3 t r i a l 4 ( G 3 + K ) Δ γ b 2 sin ψ + 2 [ G ( 1 + 1 3 sin ψ ) 2 K sin ψ ] Δ γ a 2
By substituting Equation (45) into Equation (44) and making some transformations, we obtain:
a Δ γ b 2 + b 2 Δ γ a 2 = f 2 t r i a l ( σ t r i a l , c ) .
From Equations (42) and (46), the scalar forms of Δ γ a 2 and Δ γ b 2 are:
{ Δ γ a 2 = a f 1 t r i a l b 2 f 2 t r i a l D e t 2 Δ γ b 2 = a f 2 t r i a l b 2 f 1 t r i a l D e t 2
or
{ Δ γ a 2 = f 1 t r i a l b 2 Δ γ b 2 a Δ γ b 2 = f 2 t r i a l b 2 Δ γ a 2 a
with
D e t 2 = a 2 b 2 2 .
Remark 2.
Similarly, for «physical» values of parameters K, G, ϕ and ψ, D e t 2 0 (see Figure 4).
The fixed-point algorithm applied to the left edge mechanism (Algorithm 5) is based on the constitutive law of the activated mechanism (Equation (39)). Its initial conditions are taken in the elastic domain ( σ 1 t r i a l , σ 2 t r i a l , σ 3 t r i a l , 0 , 0 ,   ε ¯ n p = 0 ).
Algorithm 5: Fixed point algorithm applied to the left edge mechanism.
  • Initialization :   Z 0 = [ σ 1 trial , σ 2 trial , σ 3 trial , 0 ,   0 , ε ¯ n p ]
  • tol :   tolerance   and   N _ max   maximum   number   of   iterations  
  • W h i l e Z Z 0 Z 0 > tol   and   n N max = 100
  • a = 4 [ G + ( G 3 + K ) sin ψ sin ϕ ]
  • b 2 = 2 G ( 1 sin ψ sin ϕ 1 3 sin ψ sin ϕ ) + 4 K sin ψ sin ϕ
  • Δ γ a 2 = ( σ 1 trial σ 3 trial ) + ( σ 1 trial + σ 3 trial ) sin ϕ 2   c ( ε ¯ n p , Δ ε ¯   p ) cos ϕ b 2 Δ γ b 2 a
  • Δ γ b 2 = ( σ 2 trial σ 3 trial ) + ( σ 2 trial + σ 3 trial ) sin ϕ 2   c ( ε ¯ n p , Δ ε ¯   p ) cos ϕ b 2 Δ γ a 2 a
  • σ 1 = σ 1 trial [ 2 G ( 1 + 1 3 sin ψ ) + 2 K sin ψ ] Δ γ a 2 + ( 4 3 G 2 K ) Δ γ b 2   sin ψ
  • σ 2 = σ 2 trial + ( 4 3 G 2 K ) Δ γ a 2   sin ψ [ 2 G ( 1 + 1 3 sin ψ ) + 2 K sin ψ ] Δ γ b 2
  • σ 3 = σ 3 trial + [ 2 G ( 1 1 3 sin ψ ) 2 K sin ψ ] ( Δ γ a 2 + Δ γ b 2 )
  • ε ¯ n + 1 p = ε ¯ n p + 2 ( Δ γ a 2 + Δ γ b 2 ) cos ϕ
  • Update   c ,   K   and   G
  • Z = [ σ 1 ,   , σ 2   , σ 3   , Δ γ a 2 ,   Δ γ b 2 , ε ¯ n + 1 p ]
  • Update     Z 0 = Z
  • EndWhile
(d)
Return mapping to the apex ( σ 1 = σ 2 = σ 3 )
Finally, in the case of the Apex, the stress is located on the hydrostatic axis. The hydrostatic pressure is then written:
p n + 1 = p n + 1 t r i a l K   Δ ε v p ,
where Δ ε v p represents the volumetric plastic strain increment and is given by:
Δ ε v p = 2 sin ( ψ ) i = 1 6 Δ γ i .
The updated stress and the accumulated plastic strain are given by [1]:
  { σ n + 1 = p n + 1   1               ε ¯ n + 1 p = ε ¯ n p + β Δ ε v p .        
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 n + 1 t r i a l , 0 and ε ¯ n p = 0 . The volumetric plastic strain increment Δ ε v p can still be written as a function of the hydrostatic pressure p using equations Equation (15) and Equation (48).
Algorithm 6: Fixed point algorithm applied to the Apex mechanism.
  • Initialization :   Z 0 = [ p n + 1 trial , 0 , ε ¯ n p ]
  • tol   tolerance   and   N _ max   maximum   number   of   iterations  
  • W h i l e Z Z 0 Z 0 > tol   and   n N max = 100  
  • β = cot ϕ / sin ψ
  • p n + 1   = p n + 1 trial K Δ ε v p
  • Δ ε v p = ( p n + 1 trial c cot ϕ ) / K
  • ε ¯ n + 1 p = ε ¯ n p + β   Δ ε v p
  • Update   c ,   K   and   G
  • Z = [ p n + 1   , Δ ε v p , ε ¯ n + 1 p ]
  • Update Z 0 = Z
  • EndWhile

2.3. 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:
D e p = d σ n + 1 d ε n + 1 e   t r i a l
where σ n + 1 represents the stress tensor at time t n + 1 .

2.3.1. 1D Illustration

Using the consistency condition ( f ˙ = 0 ) , the expression of σ ˙ is given by:
σ ˙ = E ( ε ˙ E ( E + H + K ) ε ˙ ) ,
and so, the relationship of the constitutive law linking σ ˙ to ε ˙ becomes:
σ ˙ = { E ( H + K ) ( E + H + K ) ε ˙       i f     γ ˙ > 0           E ε ˙               i f     γ ˙ = 0
We can therefore deduce the tangent modulus D e l p obtained analytically for the 1D model:
D e l p = E ( H + K ) ( E + H + K )
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 e l p   ( k ) = σ n + 1 ( k ) ε n + 1 ( k )
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 ε n + 1 ( k ) can be used:
σ n + 1 ( k ) ε n + 1 ( k )   σ n + 1 ( k ) ( ε n + 1 ( k ) p e r ) σ n + 1 ( k ) ( ε n + 1 ( k ) ) d ε n + 1 ( k )
where:
d ε n + 1 ( k ) = ω | ε n + 1 ( k ) | .
In Equation (57), ω is an appropriate small number (e.g., 10 6 ) that represents the amount of perturbation. The perturbed strain ε n + 1 ( k ) p e r is:
ε n + 1 ( k ) p e r = ε n + 1 ( k ) + ω | ε n + 1 ( k ) |
and the strain increment is
Δ ε n ( k ) p e r = Δ ε n ( k ) + ω | ε n + 1 ( k ) | .
For a given ε n + 1 ( k ) , the value of σ n + 1 ( k ) is obtained by the return mapping method presented in Algorithm 2. The perturbed stress σ n + 1 ( k ) p e r is calculated from the perturbed strain ε n + 1 ( k ) p e r using the same procedure. Finally, the tangent modulus is determined by the relationship:
D n + 1 e l p   ( k ) σ n + 1 ( k ) p e r σ n + 1 ( k ) ε n + 1 ( k ) p e r ε n + 1 ( k ) .
The illustration of this numerical method is presented in Algorithm 7.
Algorithm 7: Numerical computation of the elastoplastic tangent modulus.
  • ε n p ,   α n and q n   are   known   in   step   n   and   iteration   k   of   Newton Raphson
  • Compute σ n + 1  
  • Given strain a ε n + 1 = ε n + Δ ε n
  • Compute elastic trial stress σ n + 1 trial   and   test   for   plastic   loading
  • σ n + 1 trial = E ( ε n + 1 ε n p )
  • ξ n + 1 trial = σ n + 1 trial q n
  • f n + 1 trial = | ξ n + 1 trial | [ σ Y + K α n ]
  • if f n + 1 trial 0   then
  • Elastic step: ( ) n + 1 = ( ) n + 1 trial
  • Else plastic step, go to Algorithm
  • EndIf
  • Compute σ n + 1   p e r  
  • ω : positive scalar representing the perturbation
  • Given the perturbed strain ε n + 1 per   = ε n + 1   + ω | ε n + 1   |
  • the strain increment is: Δ ε n per = Δ ε n ( k ) + ω | ε n + 1 ( k ) |
  • Compute elastic trial stress σ n + 1 trial ,   per   and   test   for   plastic   loading
  • σ n + 1 trial ,   per = σ n   + E   Δ ε n per  
  • ξ n + 1 trial ,   per = σ n + 1 trial ,   per q n
  • f n + 1 trial ,   per = | ξ n + 1 trial ,   per | [ σ Y + K α n ]
  • if f n + 1 trial ,   per 0   then
  • Elastic step: ( ) n + 1 = ( ) n + 1 trial ,   per
  • Else plastic step, go to Algorithm
  • EndIf
  • Compute D e l p
  • D e l p = σ n + 1   p e r   σ n + 1   ε n + 1 p e r ε n + 1  

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:
( D e l p ) ( n + 1 ) ( k ) = ( σ x x ε x x e   t r i a l σ x x ε x y e   t r i a l σ x x ε y y e   t r i a l σ x y ε x x e   t r i a l σ x y ε x y e   t r i a l σ x y ε y y e   t r i a l σ y y ε x x e   t r i a l σ y y ε x y e   t r i a l σ y y ε y y e   t r i a l ) ( n + 1 ) ( k )
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.

3. 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.

3.1. 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   k N is applied at the end of the bar (Figure 5).
The yield strength is σ Y = 400   M P a . In the elastic region, the bar’s Young’s modulus is E = 200   G P a 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 e l p , Num-NR represents the numerical calculation of D e l p , 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) 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 e x t F i n t F e x t 10 6 where F e x t represents the external force, F i n t 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.

3.2. 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 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 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 (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 (Figure 11 and Figure 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.

4. 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.

Author Contributions

Conceptualization: A.S.; methodology, A.S. and G.Y.A.; software, A.S. and G.Y.A.; validation, A.S. and G.Y.A.; investigation, G.Y.A.; resources, A.S.; writing—original draft preparation, G.Y.A.; writing—review and editing, G.Y.A.; visualization, G.Y.A.; supervision, A.S.; project administration, A.S.; funding acquisition, A.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Natural Sciences and Engineering Research Council of Canada and Hydro-Québec, grant number RDCPJ470297-14.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

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.
(a)
Main plane
By substituting Equation (29) into Equation (25), after differentiation we have:
{ d σ 1 d σ 2   d σ 3 } = { d σ 1 t r d σ 2 t r d σ 3 t r } 2 Δ γ { ( K + G 3 ) sin ψ + G ( K 2 G 3 ) sin ψ               ( K + G 3 ) sin ψ G }
Knowing that
d σ j t r = D d ε j
and because the trial stress is assumed to be elastic, Equation (A2) becomes:
d σ j   = D d ε j 2   d f 1 t r i a l a V p m
with
V p m = { ( K + G 3 ) sin ψ + G ( K 2 G 3 ) sin ψ               ( K + G 3 ) sin ψ G }
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 1 t r ( σ t r , c ) = 2 V c m d ε j ,     j = 1 , 2 , 3
with
V c m = { ( K + G 3 ) sin ϕ + G ( K 2 G 3 ) sin ϕ               ( K + G 3 ) sin ϕ G }
representing the normal vector applied to the yield function of the main plane, so that
d σ j   = ( D 4   a V p m   · V c m T ) d ε j ,     j = 1 , 2 , 3
Finally, the consistent tangent operator in the principal basis and for the main plane ( D p , m e p ) obtained analytically can be written as:
D p , m e p = D 4   a V p m   . V c m T
(b)
Right edge
We replace Equation (38) in Equation (30), and then differentiate to find:
{ d σ 1 d σ 2 d σ 3 }   =   { d σ 1 t r d σ 2 t r d σ 3 t r } 2   Δ γ a 1 { ( K + G 3 ) sin ψ + G ( K 2 3 G ) sin ψ ( K + G 3 ) sin ψ G } V p m 2   Δ γ b 1 { ( K + G 3 ) sin ψ + G ( K + G 3 ) sin ψ G ( K 2 3 G ) sin ψ } V p r
where V p r represents the normal vector applied to the flow potential of the right edge, and
d σ j   = D d ε j 2   d Δ γ a 1   V p m 2   d Δ γ b 1   V p r
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) and (35) defined by:
{ d f 1 t r ( σ t r , c ) = 2 V c m d ε j d f 6 t r i a l ( σ t r i a l , c ) = 2 V c r d ε j
with
V c r = { ( K + G 3 ) sin ϕ + G ( K + G 3 ) sin ϕ G ( K 2 G 3 ) sin ϕ }
representing the normal vector applied to the yield function of the right edge and the plastic multipliers defined by
{ d Δ γ a 1 = a   d f 1 t r i a l b 1   d f 6 t r i a l D e t 1 d Δ γ b 1 = a   d f 6 t r i a l b 1   d f 1 t r i a l D e t 1     ,
the increment of the j principal stress is then written as:
d σ j   = { D 4   D e t 1 [ a ( V p m   · V c m T + V p r   · V c r T ) b 1 ( V p m   · V c r T + V p r   · V c m T ) ] } d ε j
The consistent tangent operator in the main base for the right edge ( D p , r e p ) obtained analytically is therefore:
D p , r e p = D 4   D e t 1 [ a ( V p m   · V c m T + V p r   · V c r T ) b 1 ( V p m   · V c r T + V p r   · V c m T ) ]
(c)
Left edge
As before, we substitute Equation (47) into Equation (39) and then after differentiating, to obtain:
{ d σ 1 d σ 2 d σ 3 }   =   { d σ 1 t r d σ 2 t r d σ 3 t r } 2   Δ γ a 2 { ( K + G 3 ) sin ψ + G ( K 2 3 G ) sin ψ ( K + G 3 ) sin ψ G } V p m 2   Δ γ b 2 { ( K 2 3 G ) sin ψ ( K + G 3 ) sin ψ + G ( K + G 3 ) sin ψ G } V p l
and
d σ j   = D d ε j 2   d Δ γ a 2   V p m 2   d Δ γ b 2   V p l
In addition, taking into account the differentiation of the consistency conditions given in Equations (28) and (44) and defined by
{ d f 1 t r ( σ t r , c ) = 2 V c m d ε j d f 2 t r i a l ( σ t r i a l , c ) = 2 V c l d ε j
with
V c l = { ( K 2 G 3 ) sin ϕ ( K + G 3 ) sin ϕ + G ( K + G 3 ) sin ϕ G }
and plastic multipliers given as:
{ d Δ γ a 2 = a   d f 1 t r i a l b 2   d f 2 t r i a l D e t 2 d Δ γ b 2 = a   d f 2 t r i a l b 2   d f 1 t r i a l D e t 2     ,
the increment of the main stress j is then written as:
d σ j   = { D 4   D e t 2 [ a ( V p m   · V c m T + V p l   · V c l T ) b 2 ( V p m   · V c l T + V p l   · V c m T ) ] } d ε j
Finally, the consistent tangent operator in the principal basis and for the left edge ( D p , m e p ) obtained analytically is:
D p , l e p = D 4   D e t 2 [ a ( V p m   · V c m T + V p l   · V c l T ) b 2 ( V p m   · V c l T + V p l   · V c m T ) ]
(d)
Apex
By evaluating the pressure and strain equations Equations (15) and (50), respectively, we find:
d σ j   = 0
Thus, the consistent tangent operator in the principal basis for the Apex ( D p , a e p ) is
D p , a e p = 0
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]).

References

  1. Neto, E.A.d.; Peric, D.; Owen, D.R. Computational Methods for Plasticity: Theory and Applications; John Wiley & Sons: Hoboken, NJ, USA, 2008. [Google Scholar]
  2. Geng, D.; Dai, N.; Guo, P.; Zhou, S.; Di, H. Implicit numerical integration of highly nonlinear plasticity models. Comput. Geotech. 2021, 132, 103961. [Google Scholar] [CrossRef]
  3. Scherzinger, W.M. A return mapping algorithm for isotropic and anisotropic plasticity models using a line search method. Comput. Methods Appl. Mech. Eng. 2017, 317, 526–553. [Google Scholar] [CrossRef] [Green Version]
  4. Simo, J.; Hughes, T. Computational Inelasticity; Springer: Berlin/Heidelberg, Germany, 2006; Volume 7. [Google Scholar]
  5. Simo, J.C.; Taylor, R.L. Consistent tangent operators for rate-independent elastoplasticity. Comput. Methods Appl. Mech. Eng. 1985, 48, 101–118. [Google Scholar] [CrossRef]
  6. Schanz, T.; Vermeer, P.; Bonnier, P. The hardening soil model: Formulation and verification. In Beyond 2000 in Computational Geotechnics; Ten Years of PLAXIS International. Proceedings of the International Symposium; Routledge: London, UK, 1999; pp. 281–296. [Google Scholar]
  7. Meng, C.; Tang, Z.; Chen, M.; Peng, Q. Return mapping algorithm in principal space for general isotropic elastoplasticity involving multi-surface plasticity and combined isotropic-kinematic hardening within finite deformation framework. Finite Elem. Anal. Des. 2018, 150, 1–19. [Google Scholar] [CrossRef]
  8. Clausen, J.; Damkilde, L. An exact implementation of the Hoek–Brown criterion for elasto-plastic finite element calculations. Int. J. Rock Mech. Min. Sci. 2008, 45, 831–847. [Google Scholar] [CrossRef]
  9. Griffiths, D.; Willson, S. An explicit form of the plastic matrix for a Mohr–Coulomb material. Commun. Appl. Numer. Methods 1986, 2, 523–529. [Google Scholar] [CrossRef]
  10. Bruno, H.; Barros, G.; Menezes, I.F.; Martha, L.F. Return-mapping algorithms for associative isotropic hardening plasticity using conic optimization. Appl. Math. Model. 2020, 78, 724–748. [Google Scholar] [CrossRef]
  11. Čermák, M.; Sysala, S.; Valdman, J. Efficient and flexible MATLAB implementation of 2D and 3D elastoplastic problems. Appl. Math. Comput. 2019, 355, 595–614. [Google Scholar] [CrossRef] [Green Version]
  12. Brinkgreve, R.; Kumarswamy, S.; Swolfs, W.; Waterman, D.; Chesaru, A.; Bonnier, P. Plaxis 2016; PLAXIS bv: Delft, The Netherlands, 2016. [Google Scholar]
  13. Marc, K. Loi de Mohr-Coulomb- Code_Aster. 2016. Available online: https://www.code-aster.org/V2/doc/v13/fr/man_r/r7/r7.01.28.pdf (accessed on 1 February 2021).
  14. Kim, N.H. Introduction to Nonlinear Finite Element Analysis; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2014. [Google Scholar]
  15. Simo, J.C.; Govindjee, S. Non-linear B-stability and symmetry preserving return mapping algorithms for plasticity and viscoplasticity. Int. J. Numer. Methods Eng. 1991, 31, 151–176. [Google Scholar] [CrossRef]
  16. Wilkins, M.L. Calculation of elastic-plastic flow. In Methods in Computational Physics; Alder, B., Fernbach, S., Roteerg, M., Eds.; Academic Press: New York, NY, USA, 1964; Volume 3, pp. 211–263. [Google Scholar]
  17. Nagtegaal, J.C. On the implementation of inelastic constitutive equations with special reference to large deformation problems. Comput. Methods Appl. Mech. Eng. 1982, 33, 469–484. [Google Scholar] [CrossRef]
  18. Simo, J.C.; Taylor, R.L. A return mapping algorithm for plane stress elastoplasticity. Int. J. Numer. Methods Eng. 1986, 22, 649–670. [Google Scholar] [CrossRef]
  19. Sysala, S.; Čermák, M.; Ligurský, T. Subdifferential-based implicit return-mapping operators in Mohr-Coulomb plasticity. Zamm-J. Appl. Math. Mech. /Z. Für Angew. Math. Und Mech. 2017, 97, 1502–1523. [Google Scholar] [CrossRef] [Green Version]
  20. Semenov, A.S.; Liskowsky, A.; Balke, H. Return mapping algorithms and consistent tangent operators in ferroelectroelasticity. Int. J. Numer. Methods Eng. 2010, 81, 1298–1340. [Google Scholar] [CrossRef]
  21. Hamed, A.A. Predictive Numerical Modeling of the Behavior of Rockfill Dams. Ph.D. Thesis, École de Technologie Supérieure, Montréal, QC, Canada, 2017. [Google Scholar]
  22. Ni, P.; Mei, G.; Zhao, Y.; Chen, H. Plane strain evaluation of stress paths for supported excavations under lateral loading and unloading. Soils Found. 2018, 58, 146–159. [Google Scholar] [CrossRef]
Figure 1. Mohr–Coulomb yield surface in principal stress space [1].
Figure 1. Mohr–Coulomb yield surface in principal stress space [1].
Applsci 11 04637 g001
Figure 2. The Mohr–Coulomb flow rule; (a) faces and edges, and (b) apex [1].
Figure 2. The Mohr–Coulomb flow rule; (a) faces and edges, and (b) apex [1].
Applsci 11 04637 g002
Figure 3. Example of D e t 1 values for K = 1.67 × 108 and G = 6.39 × 105.
Figure 3. Example of D e t 1 values for K = 1.67 × 108 and G = 6.39 × 105.
Applsci 11 04637 g003
Figure 4. Example of D e t 2 values for K = 1.67 × 108 and G = 6.39 × 105.
Figure 4. Example of D e t 2 values for K = 1.67 × 108 and G = 6.39 × 105.
Applsci 11 04637 g004
Figure 5. Nonlinear bar problem [14].
Figure 5. Nonlinear bar problem [14].
Applsci 11 04637 g005
Figure 6. Stress-strain for a single element considered for ω = 10 6 and t o l = 10 6 .
Figure 6. Stress-strain for a single element considered for ω = 10 6 and t o l = 10 6 .
Applsci 11 04637 g006
Figure 7. Evolution of the plasticization in the bar with to the applied load.
Figure 7. Evolution of the plasticization in the bar with to the applied load.
Applsci 11 04637 g007
Figure 8. Mesh used (a) PLAXIS: 327 triangular elements with 15 nodes and (b) Matlab: 800 quadrilateral isoparametric bilinear quadrilateral elements (40 × 20) with 4 nodes.
Figure 8. Mesh used (a) PLAXIS: 327 triangular elements with 15 nodes and (b) Matlab: 800 quadrilateral isoparametric bilinear quadrilateral elements (40 × 20) with 4 nodes.
Applsci 11 04637 g008
Figure 9. Comparisons of horizontal (ux) and vertical (uy) displacements and normal stresses calculated for ω = 10 6 and t o l = 10 10 .
Figure 9. Comparisons of horizontal (ux) and vertical (uy) displacements and normal stresses calculated for ω = 10 6 and t o l = 10 10 .
Applsci 11 04637 g009
Figure 10. Choice of the studied areas for the comparison of variables: (a) column considered and (b) row considered.
Figure 10. Choice of the studied areas for the comparison of variables: (a) column considered and (b) row considered.
Applsci 11 04637 g010
Figure 11. Displacements and stresses in the column considered for Matlab and PLAXIS models according to mesh size.
Figure 11. Displacements and stresses in the column considered for Matlab and PLAXIS models according to mesh size.
Applsci 11 04637 g011
Figure 12. Displacements and stresses in the line considered for Matlab and PLAXIS models according to mesh size.
Figure 12. Displacements and stresses in the line considered for Matlab and PLAXIS models according to mesh size.
Applsci 11 04637 g012
Table 1. Convergence rate and number of iterations by loading: 1D case.
Table 1. Convergence rate and number of iterations by loading: 1D case.
Convergence: F e x t F i n t F e x t < 10 6 Number of Iterations
Load (kN)Ana-NRNum-NRAna-NRNum-NR
50011
1004.6 × 10−1111
151.2 × 10−163.1 × 10−1111
2002.7 × 10−1111
2504.9 × 10−1211
3009.9 × 10−1311
3502.1 × 10−1111
4003.8 × 10−1211
451.62 × 10−162.3 × 10−722
5004.8 × 10−1111
Table 2. Convergence rate and number of iterations by loading: 2D case.
Table 2. Convergence rate and number of iterations by loading: 2D case.
Number   of   Iterations   vs .   to   Number   of   Layer   ( with   the   Convergence   Criterion   Given   by   F e x t F i n t F e x t 10 10 )
Meshes
20 × 1060 × 2030 × 30
Number of LayerNumber of IterationsConvergenceNumber of IterationsConvergenceNumber of IterationsConvergence
111.692 × 10−1612.708 × 10−1611.653 × 10−16
218.636 × 10−1633.191 × 10−1115.335 × 10−16
346.068 × 10−1658.266 × 10−1611.512 × 10−15
453.192 × 10−1351.234 × 10−1546.654 × 10−16
551.590 × 10−1541.718 × 10−1141.099 × 10−15
652.110 × 10−1552.896 × 10−1549.045 × 10−12
751.712 × 10−1154.163 × 10−1553.325 × 10−14
854.503 × 10−1553.024 × 10−1452.198 × 10−15
975.917 × 10−1566.844 × 10−1553.132 × 10−15
1087.336 × 10−1578.043 × 10−1554.763 × 10−15
11 71.136 × 10−1451.078 × 10−14
12 79.497 × 10−1455.637 × 10−15
13 71.948 × 10−1456.700 × 10−15
14 74.625 × 10−1456.571 × 10−15
15 72.017 × 10−1157.162 × 10−15
16 75.878 × 10−1451.670 × 10−12
17 78.057 × 10−1159.223 × 10−11
18 84.806 × 10−1461.077 × 10−14
19 82.862 × 10−1351.517 × 10−11
20 81.182 × 10−1151.321 × 10−14
21 51.427 × 10−14
22 53.637 × 10−13
23 65.778 × 10−12
24 71.794 × 10−14
25 71.791 × 10−14
26 72.160 × 10−14
27 71.100 × 10−13
28 74.871 × 10−12
29 82.473 × 10−14
30 82.938 × 10−14
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Amouzou, G.Y.; Soulaïmani, A. Numerical Algorithms for Elastoplacity: Finite Elements Code Development and Implementation of the Mohr–Coulomb Law. Appl. Sci. 2021, 11, 4637. https://doi.org/10.3390/app11104637

AMA Style

Amouzou GY, Soulaïmani A. Numerical Algorithms for Elastoplacity: Finite Elements Code Development and Implementation of the Mohr–Coulomb Law. Applied Sciences. 2021; 11(10):4637. https://doi.org/10.3390/app11104637

Chicago/Turabian Style

Amouzou, Gildas Yaovi, and Azzeddine Soulaïmani. 2021. "Numerical Algorithms for Elastoplacity: Finite Elements Code Development and Implementation of the Mohr–Coulomb Law" Applied Sciences 11, no. 10: 4637. https://doi.org/10.3390/app11104637

APA Style

Amouzou, G. Y., & Soulaïmani, A. (2021). Numerical Algorithms for Elastoplacity: Finite Elements Code Development and Implementation of the Mohr–Coulomb Law. Applied Sciences, 11(10), 4637. https://doi.org/10.3390/app11104637

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop