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,
represents the elastic strain tensor,
is the elastic strain test tensor;
is the internal hardening,
the plasticity flux vector,
is all the combined thermodynamic forces associated with
,
is the generalized hardening module
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:
where
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:
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. |
|
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
at the instant
from a given strain increment
and a state
at instant
. The initial values for
are
. The iteration is stopped for a time step
when
is tess then a fixed threshold. Then, once
is calculated it is tested on the plasticity criterion as presented in [
1]. If
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. |
else EndIf 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
and its increment
. It also involves the shear stress τ and normal stress
. For this criterion, the plastic strain begins when the shear stress τ, as well as the normal stress
become critical. These two stresses are linked by the relationship:
The plastic function
f is expressed in terms of principal stresses by:
with the convention
where
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 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
(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
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:
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
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
given by:
- ▪
For two activated mechanisms, the plastic flow rule is written:
with
the normal at the intersection of the two planes
and
or
and
, 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]):
where
represents the
i-th volumetric component of the updated plastic flow vector, tr[ ] the trace of the second-order tensor and
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
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:
with
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
. In this algorithm, the plastic multiplier ∆γ is obtained by substituting the principal trial stresses
in the consistency equation for this activated mechanism Equation (28). Moreover, from Equation (28) we deduce
.
Algorithm 3: Fixed point algorithm applied to the main plane mechanism. |
Update EndWhile
|
- (b)
Return mapping to the right edge
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:
with
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
and
, defined by:
or
with
Remark 1. Regardless of the «physical» values taken by parameters K, G, ϕ and ψ,.
(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 (
).
Algorithm 4: Fixed point method applied to the right edge mechanism. |
Update EndWhile
|
- (c)
Return mapping to the left edge
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:
with
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
and
are:
or
with
Remark 2.
Similarly, for «physical» values of parameters K, G, ϕ and ψ,(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 (
).
Algorithm 5: Fixed point algorithm applied to the left edge mechanism. |
EndWhile
|
- (d)
Return mapping to the apex
Finally, in the case of the Apex, the stress is located on the hydrostatic axis. The hydrostatic pressure is then written:
where
represents the volumetric plastic strain increment and is given by:
The updated stress and the accumulated plastic strain are given by [
1]:
with
.
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
and
The volumetric plastic strain increment
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. |
Update EndWhile
|