# Implicit Damping Iterative Algorithm to Solve Elastoplastic Static and Dynamic Equations

Authors: Huaifa Ma; Jikai Zhou; Guoping Liang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

This paper presents an implicit damping iterative algorithm to simultaneously solve equilibrium equations, yield function, and plastic flow equations, without requiring an explicit expression of elastoplastic stiffness matrices and local iteration for “return mapping” stresses to the yield surface. In addition, a damping factor is introduced to improve the stiffness matrix conformation in the nonlinear iterative process. The incremental iterative scheme and whole amount iterative scheme are derived to solve the dynamical and static and dynamical elastoplastic problems. To validate the proposed algorithms, computation procedures are designed and the numerical tests are implemented. The computational results verify the correctness and reliability of the proposed implicit iteration algorithms.


---

## Body

## 1. Introduction

The solution of the nonlinear problem using FEM eventually boils down to solving discretized nonlinear equations. FEM uses a series of modified linear approximate solutions to approach the solution of a nonlinear problem based on an iterative process. There are many methods available for solving nonlinear equations [1, 2], which may be divided into two categories: the direct iteration method and the Newton-Raphson method, or simply Newton method. The convergence rate of the direct iteration method is highly dependent on the choice of initial values. For problems with many degrees of freedom, instability may occur. In addition, for problems related to deformation history, the application of this method is quite limited. For these reasons, the direct method is rarely used. The Newton-Raphson method is probably the most popular method for solving nonlinear equations, and extensive researches and investigations have been performed. There are many derivatives of these methods, including the modified Newton-Raphson method, Quasi-Newton method, and incremental method (which can be regarded as the incremental form of the Newton-Raphson method). Researches on the Newton-Raphson method are focused on two aspects, namely, the computational efficiency and the stability of the solution. The Aitken acceleration method [3] and linear-search method [4, 5] are used in conjunction with Newton’s method to reduce the number of computational iterations. Wempner [6] and Riks [7] proposed the arc-length method, and Forde and Stiemer [8], Müller [9], and others improved it. These improvements made it possible to analyze and solve the ultimate load of structure as well as the material weakening problems.The elastoplastic problem is a typical nonlinear problem. In the loading and unloading processes, elastoplastic materials show different deformation characteristics: plastic hardening or weakening when loading, but elastic deformation when unloading. Two basic issues must be addressed for solving elastoplastic problems using FEM, namely, the linear scheme of nonlinear equations and its solution algorithm and the material constitutive relationship and its integration. The linear scheme of nonlinear equations and its solution process often depend on the material properties, load magnitude, loading history, and loading method. Therefore, the method of incremental stress-strain constitutive equations combined with iterative schemes is widely used, but the iterative operations will inevitably cause the effective stress to deviate from the yield surface, that is, the “drift” phenomenon.At present, the Newton-Raphson method is still an effective way to solve plasticity problems. In particular for serious nonlinear processes involving material softening, it is more effective to combine the Newton-Raphson method with the arc-length method or displacement controlled method. In order to describe the complex deformation processes of elastoplastic problems, a variety of constitutive equation integration methods [10] have been proposed, and the most widely used one is the return-mapping algorithm [11–16]. The return-mapping algorithm is actually an elastic forecast and plastic amendment process which requires the local iteration to correct the plastic parameters in the iterative process for solving nonlinear equations and return the forecast stress to the yield surface. This paper proposes a new nonlinear iteration method, that is, the implicit damping iterative method. Different from existing methods, this new method is designed to simultaneously solve equilibrium equations, yield functions, and plastic flow equations of elastoplastic static and dynamical problems. In the iterative process of certain nonlinear problems, the global stiffness matrix of a finite element formulation tends to be ill-conditioned. In order to avoid the singularity of the coefficient matrix, the authors introduce a damping factor in the numerical iterative process, which is stable, fast, and easy to use and involves no local iteration for “return mapping” stress to the yield surfaces.


## 2. Basic Concept of the Implicit Damping Iterative Algorithm

The constitutive relation of elastoplastic material is usually given in the incremental form by the yield function and flow rule. Constitutive relations are generally expressed as an explicit form, namely,Δ
σ
=
D
e
p
Δ
ε, where Δ
σ and Δ
ε are the stress increment and the strain increment, respectively, and D
e
p is the elastoplastic matrix. The elastoplastic matrix will be adjusted according to the stress-strain state. Here the displacement (strain) will be obtained by simultaneously solving the equilibrium equations, yield function, and plastic flow equations. The explicit form of elastoplastic matrix is not required in solving the plastic strain based on the loading-unloading process and the current deformation.Based on the above concepts, the implicit damping iterative method is proposed to solve the general material nonlinear problem. Supposing thatε
n, ε
n
p, and σ
n are, respectively, the strain, plastic strain, and stress at the nth load step, and Δ
ε
n, Δ
ε
n
p, and Δ
σ
n are, respectively, the strain increment, plastic strain increment, and stress increment at the n
+
1th load step, then the total strain, plastic strain, and stress are, respectively, expressed as:(1)
ε
n
+
1
=
ε
n
+
Δ
ε
n
;
ε
n
+
1
p
=
ε
n
p
+
Δ
ε
n
p
;
σ
n
+
1
=
σ
n
+
Δ
σ
n
.In then
+
1th load step, the equations are(2a)
(
σ
n
+
1
,
δ
ε
)
=
(
F
n
+
1
,
δ
u
)
Equalibrium
equations
,
f
(
σ
n
+
1
,
κ
n
+
1
)
=
0
Yield
function
,
Δ
ε
n
+
1
p
=
Δ
λ
n
∂
f
∂
σ
Flow
rules
.In the above formulas, the plastic potential function is the same as the yield function. Assuming that the yield function is a function of the stress and internal variableκ, then F is the external distributed load, and f
n
+
1, λ
n
+
1, and κ
n
+
1, respectively, represent the yield function value, plastic factor, and internal plastic variable value at n
+
1th load step. The notation (
*
,
*
) in the first equation of (2a) and the formulas in the following section indicate the inner product of two functions, and the vectors and matrix are expressed in bold characters. In the iterative process, the value of yield function (f
n
≠
0) in the former load step is allowed to not be on the yield surface, that is, f
n
≠
0, but it can be automatically corrected by the second equation, so that the value of the yield function at n
+
1th load step meets the yield condition f
n
+
1
=
0 and returns to the yield surface.Equation (2a) is rewritten as follows:
(2b)
(
σ
n
+
Δ
σ
n
,
δ
ε
)
=
(
F
n
+
1
,
δ
u
)
,
f
(
σ
n
+
Δ
σ
n
,
κ
n
+
Δ
κ
n
)
=
0
,
Δ
ε
n
+
1
p
=
Δ
λ
n
∂
f
∂
σ
.LetΔ
κ
n
=
Δ
λ
n
H
(
σ
,
κ
), where H
(
σ
,
κ
) is defined as the hardening function; thus the stress increment is expressed as follows:
(3)
Δ
σ
n
=
D
e
(
Δ
ε
n
-
Δ
ε
n
p
)
=
D
e
Δ
ε
n
-
Δ
λ
n
D
e
∂
f
∂
σ
.HereD
e is the elastic matrix. By using these notations, we then extend f
(
σ
n
+
1
,
κ
n
+
1
) at (
σ
n
,
κ
n
) to obtain(4)
f
(
σ
n
+
Δ
σ
n
,
κ
n
+
Δ
κ
n
)
=
f
n
+
(
∂
f
∂
σ
)
T
Δ
σ
n
+
∂
f
∂
κ
Δ
κ
n
=
(
∂
f
∂
σ
)
T
(
D
e
Δ
ε
n
-
Δ
λ
n
D
e
∂
f
∂
σ
)
+
Δ
λ
n
∂
f
∂
κ
H
(
σ
,
κ
)
=
f
n
+
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
-
Δ
λ
n
{
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
-
H
(
σ
,
κ
)
∂
f
∂
κ
}
=
0
.Here and in the following sections the superscriptT in (
*
)
T indicates the transposition of a matrix:
(5)
Set
A
0
=
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
-
H
(
σ
,
κ
)
∂
f
∂
κ
.The plastic factor can be calculated by the following formula:(6)
Δ
λ
n
=
1
A
0
[
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
+
f
n
]
.Substituting (6) into (3), we obtain
(7)
Δ
σ
n
=
D
e
(
Δ
ε
n
-
Δ
ε
n
p
)
=
D
e
Δ
ε
n
-
Δ
λ
n
D
e
∂
f
∂
σ
=
D
e
Δ
ε
n
-
1
A
0
D
e
∂
f
∂
σ
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
+
f
n
)
=
[
D
e
-
1
A
0
D
e
∂
f
∂
σ
(
∂
f
∂
σ
)
T
D
e
]
Δ
ε
n
-
f
n
A
0
D
e
∂
f
∂
σ
.That is,(8)
Δ
σ
n
=
D
e
p
Δ
ε
n
-
f
n
A
0
D
e
∂
f
∂
σ
,
where
(9)
D
e
p
=
[
D
e
-
1
A
0
D
e
∂
f
∂
σ
(
∂
f
∂
σ
)
T
D
e
]
.
As usual, D
e
p is referred to as the elastoplastic matrix. If f
n
=
0, the stress increments are given by Δ
σ
n
=
D
e
p
Δ
ε
n.The first term in (2b) can be expressed as follows:
(10)
(
Δ
σ
n
,
δ
ε
)
=
(
F
n
+
1
,
δ
u
)
-
(
σ
n
,
δ
ε
)
.Substituting (7) into the above equation, we obtain
(11)
(
D
e
Δ
ε
n
,
δ
ε
)
-
(
Δ
λ
n
D
e
∂
f
∂
σ
,
δ
ε
)
=
(
F
n
+
1
,
δ
u
)
-
(
σ
n
,
δ
ε
)
.EliminatingΔ
λ
n by (6), the basic equation to solve the static elastoplastic problem can be expressed in the integral weak forms as follows:(12a)
(
D
e
Δ
ε
n
,
δ
ε
)
-
(
1
A
0
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
1
A
0
f
n
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
,
δ
ε
)
.Equation (12a) can be rewritten explicitly as follows:
(12b)
(
D
e
p
Δ
ε
n
,
δ
ε
)
=
(
f
n
A
0
(
∂
f
∂
σ
)
T
D
e
,
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
,
δ
ε
)
.In accordance with Newton iteration method, at thek
+
1th iterative step the nonlinear equation f
(
x
(
k
+
1
)
)
=
0, after linearization at x
(
k
), f
(
x
(
k
)
)
+
d
f
(
x
(
k
)
)
/
d
x
Δ
x
(
k
)
=
0.Herex
(
k
+
1
)
=
x
(
k
)
+
Δ
x
(
k
), when Δ
x
(
k
)
→
0, x
(
k
+
1
) is the solution of the problem. Substituting d
f
(
x
(
k
)
)
/
d
x by coefficients A, the equation f
(
x
(
k
)
)
+
d
f
(
x
(
k
)
)
/
d
x
Δ
x
(
k
)
=
0 will be changed into f
(
x
(
k
)
)
+
A
Δ
x
(
k
)
=
0, as long as Δ
x
(
k
)
→
0, and x
(
k
+
1
)
=
x
(
k
)
+
Δ
x
(
k
); x
(
k
+
1
) also tends to the solution of the problem. To avoid the singularity of the coefficient matrix of the equation, the damping factor μ is introduced, so that A
=
A
0
+
μ; that is,
(13)
A
=
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
-
∂
f
∂
κ
H
+
μ
.The value ofμ can be modified according to actual situations. We attempt to keep the absolute value of |
A
| far from zero but not so high that it will affect the iterative convergence rate.(1)
If(14)
A
0
=
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
-
∂
f
∂
κ
H
>
0
,
that is,(15)
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
>
∂
f
∂
κ
H
,
set
μ
=
C
1
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
;
(2)
if(16)
A
0
=
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
-
∂
f
∂
κ
H
<
0
;
that is,(17)
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
<
∂
f
∂
κ
H
,
set
μ
=
-
C
2
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
.HereC
1 and C
2 are all real numbers greater than zero. When A
0 is equal to or approach zero, then μ can be given a greater number. In this paper, C
1 and C
2 are taken as 2, and the damping factor μ is determined by the following formula:
(18)
μ
=
{
2
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
∂
f
∂
κ
H
<
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
-
2
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
∂
f
∂
κ
H
>
(
∂
f
∂
σ
)
T
D
e
∂
f
∂
σ
.The constants in (18) are dependent on actual cases. The forms of the damping factors may be different from the yield function and flow rules. It is not difficult to extend this implicit damping iterative method to solve general material nonlinear problems with stiff damage weakening (degradation), but the representation of the yield function and flow rule may be different.


## 3. Iterative Schemes

### 3.1. Construction of Iterative Formulas

The unknowns in (12a) are displacement increments. In the elastic loading or plastic unloading, the plastic factor Δ
λ
n is less than zero; then (12a) is transformed into the general linear elastic equation: (
D
e
Δ
ε
n
,
δ
ε
)
=
(
F
n
+
1
,
δ
u
)
-
(
σ
n
,
δ
ε
); when the stress-strain falls into the plastic zone, the plastic factor Δ
λ
n is greater than zero. Therefore, one may solve elastic equations if in the iterative process the plasticity actor value is negative or the plastic equation is positive.The following section introduces two iterative methods: one is an incremental iterative method, and the other is the full amount displacements of iteration.There are two algorithms for displacement increment iteration:(
1
) the iteration algorithm based on independent unknowns of strain increments and plastic factor and (
2
) the numerical solution based on the only unknowns of strain increment.


#### 3.1.1. Increment Iterative Scheme (I)

At each load step, the external load remains constant, by (12a) after inducing the damping factor at the kth iteration:
(19)
(
D
e
Δ
Δ
ε
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
Δ
Δ
ε
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
1
A
f
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.We then obtain the displacement subincrementsΔ
Δ
u
(
k
) by the iterative equation (19) at the kth iteration. In fact, (19) can be converted to the following equivalent form:
(20)
(
D
e
Δ
ε
n
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
D
e
Δ
ε
n
(
k
-
1
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
1
A
f
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.From iterative equation (20), we can directly obtain the total displacement increments Δ
u
n
(
k
) and Δ
ε
n
(
k
) at the current load step. We then compute stresses σ
n
+
1
(
k
) by means of the following formula:
(21)
σ
n
+
1
(
k
)
=
σ
n
+
D
e
Δ
ε
n
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
f
n
(
σ
n
,
κ
n
)
.Note that when calculating the stresses by (21), A
0 is used, rather than A, to compute the real stresses. We then compute ∂
f
/
∂
σ by difference, and let ∂
f
/
∂
σ
i
j
≈
[
f
(
σ
i
j
+
δ
)
-
f
(
σ
i
j
)
]
/
δ. In this paper δ takes the 10
-
3 of the material ultimate stress. When the iteration converges, Δ
ε
n
(
k
)
→
Δ
ε
n, ε
n
+
1
(
k
)
=
ε
n
+
1
(
k
-
1
)
+
Δ
ε
n
(
k
)
→
ε
n
+
1, and σ
n
+
1
(
k
)
→
σ
n
+
1, the external loads are in equilibrium with internal forces, and the equilibrium equations tend to the following form:
(22)
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
,
δ
ε
)
=
0
.The procedure of increment iterative scheme (I) is listed in Box1.Box 1:Increment iterative scheme (I).
(i) The given initial values includeu
0, Δ
u
0, σ
0, κ
0, material and iterative controlling parameters;
(ii) According to the acting loadsF
n
+
1, execute iteration within n
+
1th loading step, to find displacement increments
Δ
u
n
(
k
) by (20) until achieving convergence:
(1) For thekth iteration:
We assumeu
n
+
1
(
0
)
=
u
n; u
n
+
1
(
0
)
=
σ
n; and Δ
u
n
+
1
(
0
)
=
0, then execute the iterative loops by (20):
① Compute the strain incrementsΔ
ε
n
(
k
), stresses σ
n
+
1
(
k
), yield function value f
n
+
1
(
k
), ∂
f
/
∂
σ and the other coefficients.
② Compute the plastic parameterd
λ
=
(
∂
f
/
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
+
f
n
(
σ
n
,
κ
n
), and determine the stress state:
Ifd
λ
>
0, the stress is located in the plastic area; if d
λ
≤
0, then the stress is in the elastic area.
③ Update stiffness matrix by the following regulations:
(a) In the elastic area, eliminate the terms containingA in  (20).
(b) In the plastic area, compute A by  (13) and (18).
④ SolveΔ
u
n
(
k
) by  (20), calculate e
r
r
=
∥
Δ
u
n
(
k
)
-
Δ
u
n
(
k
-
1
)
∥.
⑤ Determine whether or noterr is less than α
∥
Δ
u
n
(
k
)
∥, (α
=
10
-
9 in this paper).
If not, repeat Steps ①–④, and conduct the k
+
1th iteration, until achieving convergence. Then return to Step (2).
(2) Re-compute the plastic parameterd
λ by the displacements at the final iteration, determine the stress state,
compute the total stressesσ
n
+
1 by  (21), as well as inner variable κ
n
+
1, plastic strains, and so on.
(iii) Repeat Step (ii), continue the computation of the next load step, determine whether or not failure takes place,
and obtain the  ultimate load.


#### 3.1.2. Increment Iterative Scheme (II)

Equation (6) is written in the form of weak form and associated with (11):
(23)
(
D
e
Δ
ε
n
,
δ
ε
)
-
(
Δ
λ
n
D
e
∂
f
∂
σ
,
δ
ε
)
=
(
F
n
+
1
,
δ
u
)
-
(
σ
n
,
δ
ε
)
,
(
A
0
Δ
λ
n
,
δ
λ
)
-
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
,
δ
λ
)
=
(
f
n
,
δ
λ
)
.The strain incrementsΔ
u and the plastic factor Δ
λ
n are taken as an independent unknown, and these are simultaneously solved. The iterative scheme’s integral weak form is as follows:
(24)
(
D
e
Δ
ε
n
(
k
)
,
δ
ε
)
-
(
Δ
λ
n
(
k
)
D
e
∂
f
∂
σ
,
δ
ε
)
+
(
A
Δ
λ
n
(
k
)
,
δ
λ
)
-
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
,
δ
λ
)
=
(
D
e
Δ
ε
n
(
k
-
1
)
,
δ
ε
)
-
(
Δ
λ
n
(
k
-
1
)
D
e
∂
f
∂
σ
,
δ
ε
)
+
(
A
Δ
λ
n
(
k
-
1
)
,
δ
λ
)
-
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
-
1
)
,
δ
λ
)
+
(
f
n
+
1
(
k
-
1
)
,
δ
λ
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.In the iterative process of (24), if the plastic factor Δ
λ
n
(
k
) is less than zero, let it be equal to zero, and then calculate the stress increments according to the elastic matrix. In this paper, if the calculated A is very small (e.g., 10
-
10), a large number is set (e.g., 1020). The main difference between the iteration scheme (II) and the program flow of the previous iteration (I) lies in the processing of the plastic factor. Corresponding to Box 1, in Box 2 we only list a portion of the nonlinear iterative process.Box 2:Increment iterative scheme (II).
For thekth iteration:
We assume thatu
n
+
1
(
0
)
=
u
n, σ
n
+
1
(
0
)
=
σ
n and Δ
u
n
+
1
(
0
)
=
0, then execute the iterative loops by (24):
① Compute the total strain incrementsΔ
ε
n
(
k
), stresses σ
n
+
1
(
k
), yield function value f
n
+
1
(
k
), ∂
f
/
∂
σ and the other coefficients.
② ComputeA by (13) and (18). If A is very small (e.g. <10
-
10), let A equal a high number (e.g. 10
20), then obtain
Δ
u
n
(
k
) and Δ
λ
n
(
k
) by (24).
③ Determine the stress state: IfΔ
λ
n
(
k
)
<
0, let Δ
λ
n
(
k
)
=
0.
④ Calculate the errors ofΔ
u
n
(
k
) and Δ
λ
n
(
k
), and see if they satisfy the tolerance error requirement.
⑤ If not, repeat Steps  ①–④, then conduct the k
+
1th iteration, until achieving convergence.Equation (20) is rewritten as the total strain (displacement) iterative formula for a structural system to withstand the current loads:
(25)
(
D
e
ε
n
+
1
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
n
+
1
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
D
e
ε
n
+
1
(
k
-
1
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
1
A
f
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.Clearly, (25) is essentially another form of the iteration equation (20). ε
n
+
1
(
k
-
1
) and ε
n
+
1
(
k
), respectively, represent the total strain components of the k
-
1th and kth with the current load F
n
+
1. The total displacements u
n
+
1
(
k
) can be solved at the kth iterative step by (25). We then compute the displacement increments Δ
u
n
(
k
)
=
u
n
+
1
(
k
)
-
u
n and strain increments Δ
ε
n
(
k
)
=
ε
n
+
1
(
k
)
-
ε
n at the current load step and then substitute Δ
ε
n
(
k
) into (21), after which the total stress components of the current iteration step can be obtained. When the iteration converges, that is, ε
n
+
1
(
k
)
→
ε
n
+
1
(
k
-
1
), f
n
+
1
(
k
)
→
0, and σ
n
+
1
(
k
)
→
σ
n
+
1, then the same equilibrium as that in (22) is reached. The nonlinear iteration procedure of whole amount iterative scheme is almost the same as that of the increment iterative scheme (I).In (20) the displacement increments are taken as independent unknowns, and in (24) the plastic factors also act as independent unknowns. However, both of these incremental iterative formulas are completely equivalent, except for some slight differences in their expressions. Equation (25), the total displacement iteration scheme of displacements, is also equivalent to the incremental iteration schemes, by which the total strains (displacements) can be directly obtained under the current loads.


### 3.2. FEPG Script Codes

In order to verify the above iteration algorithm, the FEM FORTRAN program for the above iterative algorithm is developed based on the Finite Element Program Generator (FEPG) [17, 18]. The essence of FEPG lies in adopting both component programming techniques and human intelligence technology to automatically generate finite element source code from the given partial differential equation (PDE) codes and algorithm expressions (nonlinear finite element (NFE)). FEPG is based on the unified mathematical theory for FEM and its intrinsic rule, and the program generation is similar to the deduction process of mathematical formulas. The finite element computation program is composed of five components, namely, START, BFT, SOLV, E, and U, among which E and U are automatically generated according to the expressions given by users, while the other three components, which do not change with the expressions, are given by the system.According to the language rules of FEPG, the integral weak equations (20), (24), and (25) are written by the PDE; the iterative algorithm and loops are written by NFE codes; and the coupling relation of the displacement with stress fields is described by the Generate Command stream file and NFE codes (GCN) and either Generate IO (GIO) or Multidisciplinary Input (MDI) codes. Running the GIO or MDI command can generate the source codes for the calculation procedures [17]. In the following section, these source codes are executed to assess the iterative procedures previously established.


### 3.3. Validation of Iteration Method

We take a circular tunnel excavation problem as a testing example. Tunnel excavation will cause rock stress redistribution. The tunnel surrounding may enter the plastic state. Similar to the thick-walled cylinder, the out boundary will be subject to uniform external pressure, as shown in Figure1. Adopting the Mohr-Coulomb rule, here we let the cohesive strength be equal to 1.0 kPa and let the circular orifice radius be 1.0 m, and then the expansion processes of the plastic zone are simulated. Taking the friction coefficients f as 0.2 and 0.4, the radius of the elastoplastic interface is computed. The curves of plastic radii via pressure are given in Figure 2. As a whole, the plastic zone radii increase with confining pressure and agree well with the theoretical solutions. When friction coefficient is 0.4, especially, the computational results coincide almost exactly with the theoretical solution.Figure 1
Plastic area.Figure 2
Plastic radii via pressure.It is worth noting that far away from the outer boundary of the rectangle area, the plastic area is in good agreement with the theoretical solutions, but when the plastic area reaches the outer edge, it no longer meets the infinite boundary conditions of the theoretical solution, which may result in a larger calculated error.


## 3.1. Construction of Iterative Formulas

The unknowns in (12a) are displacement increments. In the elastic loading or plastic unloading, the plastic factor Δ
λ
n is less than zero; then (12a) is transformed into the general linear elastic equation: (
D
e
Δ
ε
n
,
δ
ε
)
=
(
F
n
+
1
,
δ
u
)
-
(
σ
n
,
δ
ε
); when the stress-strain falls into the plastic zone, the plastic factor Δ
λ
n is greater than zero. Therefore, one may solve elastic equations if in the iterative process the plasticity actor value is negative or the plastic equation is positive.The following section introduces two iterative methods: one is an incremental iterative method, and the other is the full amount displacements of iteration.There are two algorithms for displacement increment iteration:(
1
) the iteration algorithm based on independent unknowns of strain increments and plastic factor and (
2
) the numerical solution based on the only unknowns of strain increment.


### 3.1.1. Increment Iterative Scheme (I)

At each load step, the external load remains constant, by (12a) after inducing the damping factor at the kth iteration:
(19)
(
D
e
Δ
Δ
ε
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
Δ
Δ
ε
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
1
A
f
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.We then obtain the displacement subincrementsΔ
Δ
u
(
k
) by the iterative equation (19) at the kth iteration. In fact, (19) can be converted to the following equivalent form:
(20)
(
D
e
Δ
ε
n
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
D
e
Δ
ε
n
(
k
-
1
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
1
A
f
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.From iterative equation (20), we can directly obtain the total displacement increments Δ
u
n
(
k
) and Δ
ε
n
(
k
) at the current load step. We then compute stresses σ
n
+
1
(
k
) by means of the following formula:
(21)
σ
n
+
1
(
k
)
=
σ
n
+
D
e
Δ
ε
n
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
f
n
(
σ
n
,
κ
n
)
.Note that when calculating the stresses by (21), A
0 is used, rather than A, to compute the real stresses. We then compute ∂
f
/
∂
σ by difference, and let ∂
f
/
∂
σ
i
j
≈
[
f
(
σ
i
j
+
δ
)
-
f
(
σ
i
j
)
]
/
δ. In this paper δ takes the 10
-
3 of the material ultimate stress. When the iteration converges, Δ
ε
n
(
k
)
→
Δ
ε
n, ε
n
+
1
(
k
)
=
ε
n
+
1
(
k
-
1
)
+
Δ
ε
n
(
k
)
→
ε
n
+
1, and σ
n
+
1
(
k
)
→
σ
n
+
1, the external loads are in equilibrium with internal forces, and the equilibrium equations tend to the following form:
(22)
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
,
δ
ε
)
=
0
.The procedure of increment iterative scheme (I) is listed in Box1.Box 1:Increment iterative scheme (I).
(i) The given initial values includeu
0, Δ
u
0, σ
0, κ
0, material and iterative controlling parameters;
(ii) According to the acting loadsF
n
+
1, execute iteration within n
+
1th loading step, to find displacement increments
Δ
u
n
(
k
) by (20) until achieving convergence:
(1) For thekth iteration:
We assumeu
n
+
1
(
0
)
=
u
n; u
n
+
1
(
0
)
=
σ
n; and Δ
u
n
+
1
(
0
)
=
0, then execute the iterative loops by (20):
① Compute the strain incrementsΔ
ε
n
(
k
), stresses σ
n
+
1
(
k
), yield function value f
n
+
1
(
k
), ∂
f
/
∂
σ and the other coefficients.
② Compute the plastic parameterd
λ
=
(
∂
f
/
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
+
f
n
(
σ
n
,
κ
n
), and determine the stress state:
Ifd
λ
>
0, the stress is located in the plastic area; if d
λ
≤
0, then the stress is in the elastic area.
③ Update stiffness matrix by the following regulations:
(a) In the elastic area, eliminate the terms containingA in  (20).
(b) In the plastic area, compute A by  (13) and (18).
④ SolveΔ
u
n
(
k
) by  (20), calculate e
r
r
=
∥
Δ
u
n
(
k
)
-
Δ
u
n
(
k
-
1
)
∥.
⑤ Determine whether or noterr is less than α
∥
Δ
u
n
(
k
)
∥, (α
=
10
-
9 in this paper).
If not, repeat Steps ①–④, and conduct the k
+
1th iteration, until achieving convergence. Then return to Step (2).
(2) Re-compute the plastic parameterd
λ by the displacements at the final iteration, determine the stress state,
compute the total stressesσ
n
+
1 by  (21), as well as inner variable κ
n
+
1, plastic strains, and so on.
(iii) Repeat Step (ii), continue the computation of the next load step, determine whether or not failure takes place,
and obtain the  ultimate load.


### 3.1.2. Increment Iterative Scheme (II)

Equation (6) is written in the form of weak form and associated with (11):
(23)
(
D
e
Δ
ε
n
,
δ
ε
)
-
(
Δ
λ
n
D
e
∂
f
∂
σ
,
δ
ε
)
=
(
F
n
+
1
,
δ
u
)
-
(
σ
n
,
δ
ε
)
,
(
A
0
Δ
λ
n
,
δ
λ
)
-
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
,
δ
λ
)
=
(
f
n
,
δ
λ
)
.The strain incrementsΔ
u and the plastic factor Δ
λ
n are taken as an independent unknown, and these are simultaneously solved. The iterative scheme’s integral weak form is as follows:
(24)
(
D
e
Δ
ε
n
(
k
)
,
δ
ε
)
-
(
Δ
λ
n
(
k
)
D
e
∂
f
∂
σ
,
δ
ε
)
+
(
A
Δ
λ
n
(
k
)
,
δ
λ
)
-
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
,
δ
λ
)
=
(
D
e
Δ
ε
n
(
k
-
1
)
,
δ
ε
)
-
(
Δ
λ
n
(
k
-
1
)
D
e
∂
f
∂
σ
,
δ
ε
)
+
(
A
Δ
λ
n
(
k
-
1
)
,
δ
λ
)
-
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
-
1
)
,
δ
λ
)
+
(
f
n
+
1
(
k
-
1
)
,
δ
λ
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.In the iterative process of (24), if the plastic factor Δ
λ
n
(
k
) is less than zero, let it be equal to zero, and then calculate the stress increments according to the elastic matrix. In this paper, if the calculated A is very small (e.g., 10
-
10), a large number is set (e.g., 1020). The main difference between the iteration scheme (II) and the program flow of the previous iteration (I) lies in the processing of the plastic factor. Corresponding to Box 1, in Box 2 we only list a portion of the nonlinear iterative process.Box 2:Increment iterative scheme (II).
For thekth iteration:
We assume thatu
n
+
1
(
0
)
=
u
n, σ
n
+
1
(
0
)
=
σ
n and Δ
u
n
+
1
(
0
)
=
0, then execute the iterative loops by (24):
① Compute the total strain incrementsΔ
ε
n
(
k
), stresses σ
n
+
1
(
k
), yield function value f
n
+
1
(
k
), ∂
f
/
∂
σ and the other coefficients.
② ComputeA by (13) and (18). If A is very small (e.g. <10
-
10), let A equal a high number (e.g. 10
20), then obtain
Δ
u
n
(
k
) and Δ
λ
n
(
k
) by (24).
③ Determine the stress state: IfΔ
λ
n
(
k
)
<
0, let Δ
λ
n
(
k
)
=
0.
④ Calculate the errors ofΔ
u
n
(
k
) and Δ
λ
n
(
k
), and see if they satisfy the tolerance error requirement.
⑤ If not, repeat Steps  ①–④, then conduct the k
+
1th iteration, until achieving convergence.Equation (20) is rewritten as the total strain (displacement) iterative formula for a structural system to withstand the current loads:
(25)
(
D
e
ε
n
+
1
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
n
+
1
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
D
e
ε
n
+
1
(
k
-
1
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
1
A
f
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.Clearly, (25) is essentially another form of the iteration equation (20). ε
n
+
1
(
k
-
1
) and ε
n
+
1
(
k
), respectively, represent the total strain components of the k
-
1th and kth with the current load F
n
+
1. The total displacements u
n
+
1
(
k
) can be solved at the kth iterative step by (25). We then compute the displacement increments Δ
u
n
(
k
)
=
u
n
+
1
(
k
)
-
u
n and strain increments Δ
ε
n
(
k
)
=
ε
n
+
1
(
k
)
-
ε
n at the current load step and then substitute Δ
ε
n
(
k
) into (21), after which the total stress components of the current iteration step can be obtained. When the iteration converges, that is, ε
n
+
1
(
k
)
→
ε
n
+
1
(
k
-
1
), f
n
+
1
(
k
)
→
0, and σ
n
+
1
(
k
)
→
σ
n
+
1, then the same equilibrium as that in (22) is reached. The nonlinear iteration procedure of whole amount iterative scheme is almost the same as that of the increment iterative scheme (I).In (20) the displacement increments are taken as independent unknowns, and in (24) the plastic factors also act as independent unknowns. However, both of these incremental iterative formulas are completely equivalent, except for some slight differences in their expressions. Equation (25), the total displacement iteration scheme of displacements, is also equivalent to the incremental iteration schemes, by which the total strains (displacements) can be directly obtained under the current loads.


## 3.1.1. Increment Iterative Scheme (I)

At each load step, the external load remains constant, by (12a) after inducing the damping factor at the kth iteration:
(19)
(
D
e
Δ
Δ
ε
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
Δ
Δ
ε
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
1
A
f
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.We then obtain the displacement subincrementsΔ
Δ
u
(
k
) by the iterative equation (19) at the kth iteration. In fact, (19) can be converted to the following equivalent form:
(20)
(
D
e
Δ
ε
n
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
D
e
Δ
ε
n
(
k
-
1
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
1
A
f
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.From iterative equation (20), we can directly obtain the total displacement increments Δ
u
n
(
k
) and Δ
ε
n
(
k
) at the current load step. We then compute stresses σ
n
+
1
(
k
) by means of the following formula:
(21)
σ
n
+
1
(
k
)
=
σ
n
+
D
e
Δ
ε
n
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
f
n
(
σ
n
,
κ
n
)
.Note that when calculating the stresses by (21), A
0 is used, rather than A, to compute the real stresses. We then compute ∂
f
/
∂
σ by difference, and let ∂
f
/
∂
σ
i
j
≈
[
f
(
σ
i
j
+
δ
)
-
f
(
σ
i
j
)
]
/
δ. In this paper δ takes the 10
-
3 of the material ultimate stress. When the iteration converges, Δ
ε
n
(
k
)
→
Δ
ε
n, ε
n
+
1
(
k
)
=
ε
n
+
1
(
k
-
1
)
+
Δ
ε
n
(
k
)
→
ε
n
+
1, and σ
n
+
1
(
k
)
→
σ
n
+
1, the external loads are in equilibrium with internal forces, and the equilibrium equations tend to the following form:
(22)
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
,
δ
ε
)
=
0
.The procedure of increment iterative scheme (I) is listed in Box1.Box 1:Increment iterative scheme (I).
(i) The given initial values includeu
0, Δ
u
0, σ
0, κ
0, material and iterative controlling parameters;
(ii) According to the acting loadsF
n
+
1, execute iteration within n
+
1th loading step, to find displacement increments
Δ
u
n
(
k
) by (20) until achieving convergence:
(1) For thekth iteration:
We assumeu
n
+
1
(
0
)
=
u
n; u
n
+
1
(
0
)
=
σ
n; and Δ
u
n
+
1
(
0
)
=
0, then execute the iterative loops by (20):
① Compute the strain incrementsΔ
ε
n
(
k
), stresses σ
n
+
1
(
k
), yield function value f
n
+
1
(
k
), ∂
f
/
∂
σ and the other coefficients.
② Compute the plastic parameterd
λ
=
(
∂
f
/
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
+
f
n
(
σ
n
,
κ
n
), and determine the stress state:
Ifd
λ
>
0, the stress is located in the plastic area; if d
λ
≤
0, then the stress is in the elastic area.
③ Update stiffness matrix by the following regulations:
(a) In the elastic area, eliminate the terms containingA in  (20).
(b) In the plastic area, compute A by  (13) and (18).
④ SolveΔ
u
n
(
k
) by  (20), calculate e
r
r
=
∥
Δ
u
n
(
k
)
-
Δ
u
n
(
k
-
1
)
∥.
⑤ Determine whether or noterr is less than α
∥
Δ
u
n
(
k
)
∥, (α
=
10
-
9 in this paper).
If not, repeat Steps ①–④, and conduct the k
+
1th iteration, until achieving convergence. Then return to Step (2).
(2) Re-compute the plastic parameterd
λ by the displacements at the final iteration, determine the stress state,
compute the total stressesσ
n
+
1 by  (21), as well as inner variable κ
n
+
1, plastic strains, and so on.
(iii) Repeat Step (ii), continue the computation of the next load step, determine whether or not failure takes place,
and obtain the  ultimate load.


## 3.1.2. Increment Iterative Scheme (II)

Equation (6) is written in the form of weak form and associated with (11):
(23)
(
D
e
Δ
ε
n
,
δ
ε
)
-
(
Δ
λ
n
D
e
∂
f
∂
σ
,
δ
ε
)
=
(
F
n
+
1
,
δ
u
)
-
(
σ
n
,
δ
ε
)
,
(
A
0
Δ
λ
n
,
δ
λ
)
-
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
,
δ
λ
)
=
(
f
n
,
δ
λ
)
.The strain incrementsΔ
u and the plastic factor Δ
λ
n are taken as an independent unknown, and these are simultaneously solved. The iterative scheme’s integral weak form is as follows:
(24)
(
D
e
Δ
ε
n
(
k
)
,
δ
ε
)
-
(
Δ
λ
n
(
k
)
D
e
∂
f
∂
σ
,
δ
ε
)
+
(
A
Δ
λ
n
(
k
)
,
δ
λ
)
-
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
)
,
δ
λ
)
=
(
D
e
Δ
ε
n
(
k
-
1
)
,
δ
ε
)
-
(
Δ
λ
n
(
k
-
1
)
D
e
∂
f
∂
σ
,
δ
ε
)
+
(
A
Δ
λ
n
(
k
-
1
)
,
δ
λ
)
-
(
(
∂
f
∂
σ
)
T
D
e
Δ
ε
n
(
k
-
1
)
,
δ
λ
)
+
(
f
n
+
1
(
k
-
1
)
,
δ
λ
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.In the iterative process of (24), if the plastic factor Δ
λ
n
(
k
) is less than zero, let it be equal to zero, and then calculate the stress increments according to the elastic matrix. In this paper, if the calculated A is very small (e.g., 10
-
10), a large number is set (e.g., 1020). The main difference between the iteration scheme (II) and the program flow of the previous iteration (I) lies in the processing of the plastic factor. Corresponding to Box 1, in Box 2 we only list a portion of the nonlinear iterative process.Box 2:Increment iterative scheme (II).
For thekth iteration:
We assume thatu
n
+
1
(
0
)
=
u
n, σ
n
+
1
(
0
)
=
σ
n and Δ
u
n
+
1
(
0
)
=
0, then execute the iterative loops by (24):
① Compute the total strain incrementsΔ
ε
n
(
k
), stresses σ
n
+
1
(
k
), yield function value f
n
+
1
(
k
), ∂
f
/
∂
σ and the other coefficients.
② ComputeA by (13) and (18). If A is very small (e.g. <10
-
10), let A equal a high number (e.g. 10
20), then obtain
Δ
u
n
(
k
) and Δ
λ
n
(
k
) by (24).
③ Determine the stress state: IfΔ
λ
n
(
k
)
<
0, let Δ
λ
n
(
k
)
=
0.
④ Calculate the errors ofΔ
u
n
(
k
) and Δ
λ
n
(
k
), and see if they satisfy the tolerance error requirement.
⑤ If not, repeat Steps  ①–④, then conduct the k
+
1th iteration, until achieving convergence.Equation (20) is rewritten as the total strain (displacement) iterative formula for a structural system to withstand the current loads:
(25)
(
D
e
ε
n
+
1
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
n
+
1
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
=
(
D
e
ε
n
+
1
(
k
-
1
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
1
A
f
n
+
1
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
n
+
1
,
δ
u
)
-
(
σ
n
+
1
(
k
-
1
)
,
δ
ε
)
.Clearly, (25) is essentially another form of the iteration equation (20). ε
n
+
1
(
k
-
1
) and ε
n
+
1
(
k
), respectively, represent the total strain components of the k
-
1th and kth with the current load F
n
+
1. The total displacements u
n
+
1
(
k
) can be solved at the kth iterative step by (25). We then compute the displacement increments Δ
u
n
(
k
)
=
u
n
+
1
(
k
)
-
u
n and strain increments Δ
ε
n
(
k
)
=
ε
n
+
1
(
k
)
-
ε
n at the current load step and then substitute Δ
ε
n
(
k
) into (21), after which the total stress components of the current iteration step can be obtained. When the iteration converges, that is, ε
n
+
1
(
k
)
→
ε
n
+
1
(
k
-
1
), f
n
+
1
(
k
)
→
0, and σ
n
+
1
(
k
)
→
σ
n
+
1, then the same equilibrium as that in (22) is reached. The nonlinear iteration procedure of whole amount iterative scheme is almost the same as that of the increment iterative scheme (I).In (20) the displacement increments are taken as independent unknowns, and in (24) the plastic factors also act as independent unknowns. However, both of these incremental iterative formulas are completely equivalent, except for some slight differences in their expressions. Equation (25), the total displacement iteration scheme of displacements, is also equivalent to the incremental iteration schemes, by which the total strains (displacements) can be directly obtained under the current loads.


## 3.2. FEPG Script Codes

In order to verify the above iteration algorithm, the FEM FORTRAN program for the above iterative algorithm is developed based on the Finite Element Program Generator (FEPG) [17, 18]. The essence of FEPG lies in adopting both component programming techniques and human intelligence technology to automatically generate finite element source code from the given partial differential equation (PDE) codes and algorithm expressions (nonlinear finite element (NFE)). FEPG is based on the unified mathematical theory for FEM and its intrinsic rule, and the program generation is similar to the deduction process of mathematical formulas. The finite element computation program is composed of five components, namely, START, BFT, SOLV, E, and U, among which E and U are automatically generated according to the expressions given by users, while the other three components, which do not change with the expressions, are given by the system.According to the language rules of FEPG, the integral weak equations (20), (24), and (25) are written by the PDE; the iterative algorithm and loops are written by NFE codes; and the coupling relation of the displacement with stress fields is described by the Generate Command stream file and NFE codes (GCN) and either Generate IO (GIO) or Multidisciplinary Input (MDI) codes. Running the GIO or MDI command can generate the source codes for the calculation procedures [17]. In the following section, these source codes are executed to assess the iterative procedures previously established.


## 3.3. Validation of Iteration Method

We take a circular tunnel excavation problem as a testing example. Tunnel excavation will cause rock stress redistribution. The tunnel surrounding may enter the plastic state. Similar to the thick-walled cylinder, the out boundary will be subject to uniform external pressure, as shown in Figure1. Adopting the Mohr-Coulomb rule, here we let the cohesive strength be equal to 1.0 kPa and let the circular orifice radius be 1.0 m, and then the expansion processes of the plastic zone are simulated. Taking the friction coefficients f as 0.2 and 0.4, the radius of the elastoplastic interface is computed. The curves of plastic radii via pressure are given in Figure 2. As a whole, the plastic zone radii increase with confining pressure and agree well with the theoretical solutions. When friction coefficient is 0.4, especially, the computational results coincide almost exactly with the theoretical solution.Figure 1
Plastic area.Figure 2
Plastic radii via pressure.It is worth noting that far away from the outer boundary of the rectangle area, the plastic area is in good agreement with the theoretical solutions, but when the plastic area reaches the outer edge, it no longer meets the infinite boundary conditions of the theoretical solution, which may result in a larger calculated error.


## 4. Solving the Nonlinear Dynamic Equations

### 4.1. Basic Equations of Dynamic Problems

For dynamic problems, the various physical quantities of the object are a function of time. At timet the strain, plastic strain, and stress are, respectively, expressed as ε
t, ε
t
p, σ
t, the strain increment Δ
ε
t, the plastic strain increment Δ
ε
t
p, and stress increment Δ
σ
t. At time t
+
Δ
t, a total strain, plastic strain, and stress are expressed as
(26)
ε
t
+
Δ
t
=
ε
t
+
Δ
ε
t
;
ε
t
+
Δ
t
p
=
ε
t
p
+
Δ
ε
t
p
;
σ
t
+
Δ
t
=
σ
t
+
Δ
σ
t
.By analogy with (2a) and (2b), the dynamic equations at time t
+
Δ
t are(27)
(
Δ
σ
t
,
δ
ε
)
+
(
ρ
u
¨
t
+
Δ
t
,
δ
u
)
+
(
μ
u
˙
t
+
Δ
t
,
δ
u
)
=
(
F
t
+
Δ
t
,
δ
u
)
-
(
σ
(
t
)
,
δ
ε
)
,
f
t
+
Δ
t
=
0
,
Δ
ε
p
=
Δ
λ
t
∂
f
∂
σ
.Hereρ is the mass density; μ is the damping coefficient, whose meaning is different from that of (13); u
¨ and u
˙, respectively, represent the components of the acceleration and velocities.


### 4.2. Iterative Methods of Dynamic Problems

The implicit damping iterative method proposed above can be easily expanded and applied to elastoplastic dynamic problems. By analogy with the total amount iterative algorithm (see (25)) of static problems, the total amount iterative formula of dynamic problems can be written as follows:
(28)
(
D
e
ε
t
+
Δ
t
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
t
+
Δ
t
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
ρ
u
¨
t
+
Δ
t
(
k
)
,
δ
u
)
+
(
μ
u
˙
t
+
Δ
t
(
k
)
,
δ
u
)
=
(
D
e
ε
t
+
Δ
t
(
k
-
1
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
t
+
Δ
t
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
1
A
f
t
+
Δ
t
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
t
+
Δ
t
,
δ
u
)
-
(
σ
t
+
Δ
t
(
k
-
1
)
,
δ
ε
)
.The strain incrementsΔ
ε
t
(
k
)
=
ε
t
+
Δ
t
(
k
)
-
ε
t at time t
+
Δ
t up to the kth iteration. The stress is calculated by the following:
(29)
σ
t
+
Δ
t
(
k
)
=
σ
t
+
D
e
Δ
ε
t
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
(
∂
f
∂
σ
)
T
D
e
Δ
ε
t
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
f
t
(
σ
t
,
κ
t
)
.The iterative process of (25) has been rewritten as a dynamic time-dependent loading process. When ε
t
+
Δ
t
(
k
)
→
ε
t
+
Δ
t
(
k
-
1
)
=
ε
t
+
Δ
t, f
t
+
Δ
t
(
k
)
→
0, (28) converges to the following kinetic equation:
(30)
(
ρ
u
¨
t
+
Δ
t
,
δ
u
)
+
(
μ
u
˙
t
+
Δ
t
,
δ
u
)
=
(
F
t
+
Δ
t
,
δ
u
)
-
(
σ
t
+
Δ
t
,
δ
ε
)
.The iterative equation (28) is discretized in the time domain using the basic assumptions of the Newmark integration method. The dynamic finite element equations are obtained as follows:
(31)
(
K
t
(
k
-
1
)
+
a
0
M
+
a
1
C
)
u
t
+
Δ
t
(
k
)
=
K
t
(
k
-
1
)
u
t
+
Δ
t
(
k
-
1
)
+
Δ
F
t
+
Δ
t
(
k
-
1
)
+
M
[
a
0
u
t
t
+
a
2
u
˙
t
t
+
a
3
u
¨
t
t
]
+
C
[
a
1
u
t
t
+
a
4
u
˙
t
t
+
a
5
u
¨
t
t
]
.HereK
t
(
k
), M
, and C
, respectively, represent the stiffness matrix, mass matrix, and damping matrix, and Δ
F
t
+
Δ
t
(
k
) is the load increment generated by the last three terms on the right side of (28). Note that the nodal forces produced by K
t
u
t
+
Δ
t are not necessarily equal to the nodal force caused by σ
t
+
Δ
t, due to possible plastic deformation. If not, K
t
u
t
+
Δ
t
+
Δ
F
t
+
Δ
t
=
F
t
+
Δ
t and (28) tend to be the full amount of the Newmark integral equation of a linear elastic problem.Carrying numerical iterations untilu
t
+
Δ
t
=
u
t
+
Δ
t
(
k
)
≈
u
t
+
Δ
t
(
k
-
1
) and Δ
u
t
=
u
t
+
Δ
t
-
u
t, by which one can obtain the acceleration u
¨
t
+
Δ
t
=
a
0
Δ
u
t
-
a
2
u
˙
t
-
a
3
u
¨
t at the time of t
+
Δ
t and velocities u
˙
t
+
Δ
t
=
u
˙
t
+
a
6
u
¨
t
+
a
7
u
¨
t
+
Δ
t, where a
0

~

a
7 indicate the Newmark integration constants, given as follows: a
0
=
1
/
γ
Δ
t
2, a
1
=
δ
/
γ
Δ
t, a
2
=
1
/
γ
Δ
t, a
3
=
1
/
2
γ
-
1, a
4
=
δ
/
γ
-
1, a
5
=
(
Δ
t
/
2
)
(
δ
/
γ
-
2
), a
6
=
Δ
t
(
1
-
γ
), a
7
=
γ
Δ
t, and δ
≥
0.5, γ
=
0.25
(
0.5
+
δ
)
2. Listed in Box 3 is the total amount iterative scheme of the elastoplastic dynamical equations.Box 3:Whole amount iterative scheme of elastoplastic dynamical equations.
(i) Given time intervalΔ
t, the integration constants; material and iterative controlling parameters; initial values,
includingu
0, u
˙
0, u
¨
0, σ
0, κ
0, and their increments.
(ii) According to acting loadsF
t
+
Δ
t or time-displacement history, execute iteration within and t
+
Δ
t, to find total
displacementsu
t
+
Δ
t
(
k
) by (31) until achieving convergence:
(1) For thekth nonlinear iteration (not update u
˙
t
+
Δ
t and u
¨
t
+
Δ
t for time being in nonlinear iterative loops):
Assumeu
t
+
Δ
t
(
0
)
=
u
t, σ
t
+
Δ
t
(
0
)
=
σ
t and Δ
u
t
+
Δ
t
(
0
)
=
0, then execute the iterative loops and determine u
t
+
Δ
t
(
k
) by (31):
① Compute the displacement incrementsΔ
u
t
(
k
)
=
u
t
+
Δ
t
(
k
)
-
u
t, strain increments Δ
ε
t
(
k
)
=
ε
t
(
k
)
-
ε
t, σ
t
+
Δ
t
(
k
), f
t
+
Δ
t
(
k
), ∂
f
/
∂
σ
and the other coefficients.
② Compute plastic parameterd
λ
=
(
∂
f
/
∂
σ
)
T
D
e
Δ
ε
t
(
k
)
+
f
t
(
σ
t
,
κ
t
), and determine the stress state: If d
λ
>
0, the stress is
located in the plastic area; ifd
λ
≤
0, the stress is in the elastic area.
③ Update the stiffness matrix by the following regulations:
(a) In the elastic area, eliminate the terms containing A in  (28).
(b) In the plastic area, compute A by (13) and (18).
④ Findu
t
+
Δ
t
(
k
) by (31), calculate Δ
u
t
(
k
)
=
u
t
+
Δ
t
(
k
)
-
u
t, e
r
r
=
∥
Δ
u
t
(
k
)
-
Δ
u
t
(
k
-
1
)
∥.
⑤ Determine whether or note
r
r is less than α
∥
Δ
u
t
(
k
)
∥; if not, repeat Steps ①–④, conduct the k
+
1th iteration until
achieving convergence, then return to Step (2).
(2) Updateu
t
+
Δ
t
=
u
t
+
Δ
t
(
k
), u
¨
t
+
Δ
t
=
a
0
Δ
u
t
(
k
)
-
a
2
u
˙
t
-
a
3
u
¨
t, and u
˙
t
+
Δ
t
=
u
˙
t
+
a
6
u
¨
t
+
a
7
u
¨
t
+
Δ
t.
(3) Re-compute the plastic parameterd
λ by the displacements at the final iteration, determine the stress state, and compute
the total stresses and inner variableκ
t
+
Δ
t, plastic strains and so on, for the next time step.
(iii) Repeat Step (ii), continue the computation of the next load step, determine whether or not failure take places, and obtain
the ultimate load.Shown in (28) is the implicit damping iteration scheme to solve the total displacements of elastoplastic dynamic problems. This creates the boundary conditions given in the form of whole amount displacements or loads, thus providing a great convenience to conduct the nonlinear response analyses.


### 4.3. Numerical Examples

The example is the elastoplastic bending of a simply supported beam subjected to a uniform load of intensityq with the rectangular cross-section (b
×
2
h), as shown in Figure 3. Using the perfectly elastoplastic model and Mises’ yield criterion, when h
e
=
h, the stresses on the upper and lower edges of the middle part of the beam just reach the plastic limit. If the cross-section is in the fully elastic state, then the elastic limit moment at the middle beam is M
e
=
(
2
/
3
)
b
h
2
σ
s; when h
e
=
0, the cross-section of the beam fully accesses the full plastic state, formatting a plastic hinge, and then the plastic limit moments are M
p
=
b
h
2
σ
s and M
p
=
1.5
M
e.Theoretical plastic area of flexural beam under static load.
(a)
Theoretical plastic area (M
p
=
1.5
M
e)
(b)
Stress in the cross-sectionThe above theoretical solutions are obtained without gravity weight. For the convenience of comparison with the theoretical solutions, this study does not consider the weight of the beam and instead considers the mass of the beam material in the dynamic calculation and takes the damping proportional to mass; that is,C
=
α
M, with the dynamic loading rate of 100.0 kN/s.We conduct static numerical simulations by exerting static distributed load with the increments of0.01
q
s, where q
s (100 kN/m2) denotes the plastic limit load (M
p
=
1.5
M
e). In order to verify the accuracy and efficiency of the proposed method, the numerical tests were executed, respectively, by using the Implicit Damping Iterative Algorithm (IDIA) and the Return-Mapping Algorithm provided by the ABAQUS software. The curves of deflection via the loading history are obtained as shown in Figure 4. It is shown that the two curves are almost identical. Then, by exerting static distributed load with the increments of the increments of 0.10
q
s, 0.05
q
s, and 0.01
q
s, respectively, we obtain the three curves using IDIA, as shown in Figure 5, which almost completely overlap with each other. The plastic zones as shown in Figure 6, which are calculated by using both IDIA and ABAQUS, under the plastic limit moments, are all very close to the theoretical solution. Furthermore, we noted that the plastic zone calculated under the plastic limit load is not affected by loading increments. Those numerical tests demonstrate that the proposed implicit iteration algorithm has both strong stability and reliability.Figure 4
Computational results by IDIA compared with ABAQUS.Figure 5
Stability analysis of IDIA.Plastic area under static load (M
p
=
1.5
M
e).
(a)
The result computed by IDIA
(b)
The result computed by ABAQUSThe dynamical tests were conducted with 0.001 s of time step and 100.0 kN/s of loading rate. The curves of deflection under dynamic load (α
=
30.0), obtained, respectively, by IDIA and ABAQUS, as shown in Figure 7, are completely overlapped. In addition, we have similar results when the damping coefficient α is, respectively, taken as 0 and 10. The distributions of the plastic zones in Figure 8 show that the dynamic flexural strength is higher than the static flexural strength, and as damping coefficient increases, the corresponding plastic zone becomes smaller and smaller. Figures 8(a)(A) and 8(b)(A) show the case with no damping (α
=
0) and the dynamic effect induced only by the inertia force; the plastic zone is a bit smaller than that of the static load. In general, the above computational results obtained from IDIA and ABAQUS completely agree well and in line with the general rule of the dynamic effect.Figure 7
Dynamical results computed, respectively, by IDIA and ABAQUS (α
=
30.0).Plastic area under dynamical load (M
p
=
1.5
M
e), respectively, by IDIA and ABAQUS.
(a)
The result computed by IDIA
(b)
The result computed by ABAQUSThese numerical experiments are carried out on a PC configured as Intel Core i5-2400 CPU@3.10 GHz, MEM 2.99 GB. In the examples given, it iterates only 3 or 4 times to converge to solution: the first example with 16,320 DOF, 80-step calculation, taking less than 20 minutes; the second example, with 4662 DOF, running 15 loading steps to form a plastic hinge, taking less than one minute. It is shown from these calculations that the iterative algorithm proposed has good computational efficiency.In addition, a series of numerical tests performed by using the above examples showed the magnitude of the constantsC
1 and C
2 in (18) affects little the iteration convergence rate and solution accuracy.


## 4.1. Basic Equations of Dynamic Problems

For dynamic problems, the various physical quantities of the object are a function of time. At timet the strain, plastic strain, and stress are, respectively, expressed as ε
t, ε
t
p, σ
t, the strain increment Δ
ε
t, the plastic strain increment Δ
ε
t
p, and stress increment Δ
σ
t. At time t
+
Δ
t, a total strain, plastic strain, and stress are expressed as
(26)
ε
t
+
Δ
t
=
ε
t
+
Δ
ε
t
;
ε
t
+
Δ
t
p
=
ε
t
p
+
Δ
ε
t
p
;
σ
t
+
Δ
t
=
σ
t
+
Δ
σ
t
.By analogy with (2a) and (2b), the dynamic equations at time t
+
Δ
t are(27)
(
Δ
σ
t
,
δ
ε
)
+
(
ρ
u
¨
t
+
Δ
t
,
δ
u
)
+
(
μ
u
˙
t
+
Δ
t
,
δ
u
)
=
(
F
t
+
Δ
t
,
δ
u
)
-
(
σ
(
t
)
,
δ
ε
)
,
f
t
+
Δ
t
=
0
,
Δ
ε
p
=
Δ
λ
t
∂
f
∂
σ
.Hereρ is the mass density; μ is the damping coefficient, whose meaning is different from that of (13); u
¨ and u
˙, respectively, represent the components of the acceleration and velocities.


## 4.2. Iterative Methods of Dynamic Problems

The implicit damping iterative method proposed above can be easily expanded and applied to elastoplastic dynamic problems. By analogy with the total amount iterative algorithm (see (25)) of static problems, the total amount iterative formula of dynamic problems can be written as follows:
(28)
(
D
e
ε
t
+
Δ
t
(
k
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
t
+
Δ
t
(
k
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
ρ
u
¨
t
+
Δ
t
(
k
)
,
δ
u
)
+
(
μ
u
˙
t
+
Δ
t
(
k
)
,
δ
u
)
=
(
D
e
ε
t
+
Δ
t
(
k
-
1
)
,
δ
ε
)
-
(
1
A
(
∂
f
∂
σ
)
T
D
e
ε
t
+
Δ
t
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
1
A
f
t
+
Δ
t
(
k
-
1
)
,
(
∂
f
∂
σ
)
T
D
e
δ
ε
)
+
(
F
t
+
Δ
t
,
δ
u
)
-
(
σ
t
+
Δ
t
(
k
-
1
)
,
δ
ε
)
.The strain incrementsΔ
ε
t
(
k
)
=
ε
t
+
Δ
t
(
k
)
-
ε
t at time t
+
Δ
t up to the kth iteration. The stress is calculated by the following:
(29)
σ
t
+
Δ
t
(
k
)
=
σ
t
+
D
e
Δ
ε
t
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
(
∂
f
∂
σ
)
T
D
e
Δ
ε
t
(
k
)
-
1
A
0
D
e
∂
f
∂
σ
f
t
(
σ
t
,
κ
t
)
.The iterative process of (25) has been rewritten as a dynamic time-dependent loading process. When ε
t
+
Δ
t
(
k
)
→
ε
t
+
Δ
t
(
k
-
1
)
=
ε
t
+
Δ
t, f
t
+
Δ
t
(
k
)
→
0, (28) converges to the following kinetic equation:
(30)
(
ρ
u
¨
t
+
Δ
t
,
δ
u
)
+
(
μ
u
˙
t
+
Δ
t
,
δ
u
)
=
(
F
t
+
Δ
t
,
δ
u
)
-
(
σ
t
+
Δ
t
,
δ
ε
)
.The iterative equation (28) is discretized in the time domain using the basic assumptions of the Newmark integration method. The dynamic finite element equations are obtained as follows:
(31)
(
K
t
(
k
-
1
)
+
a
0
M
+
a
1
C
)
u
t
+
Δ
t
(
k
)
=
K
t
(
k
-
1
)
u
t
+
Δ
t
(
k
-
1
)
+
Δ
F
t
+
Δ
t
(
k
-
1
)
+
M
[
a
0
u
t
t
+
a
2
u
˙
t
t
+
a
3
u
¨
t
t
]
+
C
[
a
1
u
t
t
+
a
4
u
˙
t
t
+
a
5
u
¨
t
t
]
.HereK
t
(
k
), M
, and C
, respectively, represent the stiffness matrix, mass matrix, and damping matrix, and Δ
F
t
+
Δ
t
(
k
) is the load increment generated by the last three terms on the right side of (28). Note that the nodal forces produced by K
t
u
t
+
Δ
t are not necessarily equal to the nodal force caused by σ
t
+
Δ
t, due to possible plastic deformation. If not, K
t
u
t
+
Δ
t
+
Δ
F
t
+
Δ
t
=
F
t
+
Δ
t and (28) tend to be the full amount of the Newmark integral equation of a linear elastic problem.Carrying numerical iterations untilu
t
+
Δ
t
=
u
t
+
Δ
t
(
k
)
≈
u
t
+
Δ
t
(
k
-
1
) and Δ
u
t
=
u
t
+
Δ
t
-
u
t, by which one can obtain the acceleration u
¨
t
+
Δ
t
=
a
0
Δ
u
t
-
a
2
u
˙
t
-
a
3
u
¨
t at the time of t
+
Δ
t and velocities u
˙
t
+
Δ
t
=
u
˙
t
+
a
6
u
¨
t
+
a
7
u
¨
t
+
Δ
t, where a
0

~

a
7 indicate the Newmark integration constants, given as follows: a
0
=
1
/
γ
Δ
t
2, a
1
=
δ
/
γ
Δ
t, a
2
=
1
/
γ
Δ
t, a
3
=
1
/
2
γ
-
1, a
4
=
δ
/
γ
-
1, a
5
=
(
Δ
t
/
2
)
(
δ
/
γ
-
2
), a
6
=
Δ
t
(
1
-
γ
), a
7
=
γ
Δ
t, and δ
≥
0.5, γ
=
0.25
(
0.5
+
δ
)
2. Listed in Box 3 is the total amount iterative scheme of the elastoplastic dynamical equations.Box 3:Whole amount iterative scheme of elastoplastic dynamical equations.
(i) Given time intervalΔ
t, the integration constants; material and iterative controlling parameters; initial values,
includingu
0, u
˙
0, u
¨
0, σ
0, κ
0, and their increments.
(ii) According to acting loadsF
t
+
Δ
t or time-displacement history, execute iteration within and t
+
Δ
t, to find total
displacementsu
t
+
Δ
t
(
k
) by (31) until achieving convergence:
(1) For thekth nonlinear iteration (not update u
˙
t
+
Δ
t and u
¨
t
+
Δ
t for time being in nonlinear iterative loops):
Assumeu
t
+
Δ
t
(
0
)
=
u
t, σ
t
+
Δ
t
(
0
)
=
σ
t and Δ
u
t
+
Δ
t
(
0
)
=
0, then execute the iterative loops and determine u
t
+
Δ
t
(
k
) by (31):
① Compute the displacement incrementsΔ
u
t
(
k
)
=
u
t
+
Δ
t
(
k
)
-
u
t, strain increments Δ
ε
t
(
k
)
=
ε
t
(
k
)
-
ε
t, σ
t
+
Δ
t
(
k
), f
t
+
Δ
t
(
k
), ∂
f
/
∂
σ
and the other coefficients.
② Compute plastic parameterd
λ
=
(
∂
f
/
∂
σ
)
T
D
e
Δ
ε
t
(
k
)
+
f
t
(
σ
t
,
κ
t
), and determine the stress state: If d
λ
>
0, the stress is
located in the plastic area; ifd
λ
≤
0, the stress is in the elastic area.
③ Update the stiffness matrix by the following regulations:
(a) In the elastic area, eliminate the terms containing A in  (28).
(b) In the plastic area, compute A by (13) and (18).
④ Findu
t
+
Δ
t
(
k
) by (31), calculate Δ
u
t
(
k
)
=
u
t
+
Δ
t
(
k
)
-
u
t, e
r
r
=
∥
Δ
u
t
(
k
)
-
Δ
u
t
(
k
-
1
)
∥.
⑤ Determine whether or note
r
r is less than α
∥
Δ
u
t
(
k
)
∥; if not, repeat Steps ①–④, conduct the k
+
1th iteration until
achieving convergence, then return to Step (2).
(2) Updateu
t
+
Δ
t
=
u
t
+
Δ
t
(
k
), u
¨
t
+
Δ
t
=
a
0
Δ
u
t
(
k
)
-
a
2
u
˙
t
-
a
3
u
¨
t, and u
˙
t
+
Δ
t
=
u
˙
t
+
a
6
u
¨
t
+
a
7
u
¨
t
+
Δ
t.
(3) Re-compute the plastic parameterd
λ by the displacements at the final iteration, determine the stress state, and compute
the total stresses and inner variableκ
t
+
Δ
t, plastic strains and so on, for the next time step.
(iii) Repeat Step (ii), continue the computation of the next load step, determine whether or not failure take places, and obtain
the ultimate load.Shown in (28) is the implicit damping iteration scheme to solve the total displacements of elastoplastic dynamic problems. This creates the boundary conditions given in the form of whole amount displacements or loads, thus providing a great convenience to conduct the nonlinear response analyses.


## 4.3. Numerical Examples

The example is the elastoplastic bending of a simply supported beam subjected to a uniform load of intensityq with the rectangular cross-section (b
×
2
h), as shown in Figure 3. Using the perfectly elastoplastic model and Mises’ yield criterion, when h
e
=
h, the stresses on the upper and lower edges of the middle part of the beam just reach the plastic limit. If the cross-section is in the fully elastic state, then the elastic limit moment at the middle beam is M
e
=
(
2
/
3
)
b
h
2
σ
s; when h
e
=
0, the cross-section of the beam fully accesses the full plastic state, formatting a plastic hinge, and then the plastic limit moments are M
p
=
b
h
2
σ
s and M
p
=
1.5
M
e.Theoretical plastic area of flexural beam under static load.
(a)
Theoretical plastic area (M
p
=
1.5
M
e)
(b)
Stress in the cross-sectionThe above theoretical solutions are obtained without gravity weight. For the convenience of comparison with the theoretical solutions, this study does not consider the weight of the beam and instead considers the mass of the beam material in the dynamic calculation and takes the damping proportional to mass; that is,C
=
α
M, with the dynamic loading rate of 100.0 kN/s.We conduct static numerical simulations by exerting static distributed load with the increments of0.01
q
s, where q
s (100 kN/m2) denotes the plastic limit load (M
p
=
1.5
M
e). In order to verify the accuracy and efficiency of the proposed method, the numerical tests were executed, respectively, by using the Implicit Damping Iterative Algorithm (IDIA) and the Return-Mapping Algorithm provided by the ABAQUS software. The curves of deflection via the loading history are obtained as shown in Figure 4. It is shown that the two curves are almost identical. Then, by exerting static distributed load with the increments of the increments of 0.10
q
s, 0.05
q
s, and 0.01
q
s, respectively, we obtain the three curves using IDIA, as shown in Figure 5, which almost completely overlap with each other. The plastic zones as shown in Figure 6, which are calculated by using both IDIA and ABAQUS, under the plastic limit moments, are all very close to the theoretical solution. Furthermore, we noted that the plastic zone calculated under the plastic limit load is not affected by loading increments. Those numerical tests demonstrate that the proposed implicit iteration algorithm has both strong stability and reliability.Figure 4
Computational results by IDIA compared with ABAQUS.Figure 5
Stability analysis of IDIA.Plastic area under static load (M
p
=
1.5
M
e).
(a)
The result computed by IDIA
(b)
The result computed by ABAQUSThe dynamical tests were conducted with 0.001 s of time step and 100.0 kN/s of loading rate. The curves of deflection under dynamic load (α
=
30.0), obtained, respectively, by IDIA and ABAQUS, as shown in Figure 7, are completely overlapped. In addition, we have similar results when the damping coefficient α is, respectively, taken as 0 and 10. The distributions of the plastic zones in Figure 8 show that the dynamic flexural strength is higher than the static flexural strength, and as damping coefficient increases, the corresponding plastic zone becomes smaller and smaller. Figures 8(a)(A) and 8(b)(A) show the case with no damping (α
=
0) and the dynamic effect induced only by the inertia force; the plastic zone is a bit smaller than that of the static load. In general, the above computational results obtained from IDIA and ABAQUS completely agree well and in line with the general rule of the dynamic effect.Figure 7
Dynamical results computed, respectively, by IDIA and ABAQUS (α
=
30.0).Plastic area under dynamical load (M
p
=
1.5
M
e), respectively, by IDIA and ABAQUS.
(a)
The result computed by IDIA
(b)
The result computed by ABAQUSThese numerical experiments are carried out on a PC configured as Intel Core i5-2400 CPU@3.10 GHz, MEM 2.99 GB. In the examples given, it iterates only 3 or 4 times to converge to solution: the first example with 16,320 DOF, 80-step calculation, taking less than 20 minutes; the second example, with 4662 DOF, running 15 loading steps to form a plastic hinge, taking less than one minute. It is shown from these calculations that the iterative algorithm proposed has good computational efficiency.In addition, a series of numerical tests performed by using the above examples showed the magnitude of the constantsC
1 and C
2 in (18) affects little the iteration convergence rate and solution accuracy.


## 5. Conclusion

This paper first describes the basic concept of the implicit damping iterative method and then gives the displacement incremental iterative scheme as well as the whole iterative scheme to solve the static and dynamical elastoplastic problems. At the same time, the authors present the corresponding computing procedures and script files of the iterative schemes in accordance with the FEPG language rules and generate the FORTRAN programs. The circular tunnel excavation and elastoplastic bending of simply supported beam problems are numerically calculated, and the results verify the correctness and reliability of the proposed implicit iteration algorithms.The elastoplastic damping implicit iterative algorithm will simultaneously solve equilibrium equations, yield functions, and plastic flow equations. The method does not require an explicit expression of the elastoplastic stiffness matrix and local iteration for “return mapping” stress to the yield surface. In addition, the damping factor introduced in the paper improves the stiffness matrix conformation. Although the numerical algorithm proposed is based on the elastoplastic problem, it can easily be expanded and applied to solve the general material nonlinearity problem. In particular, the whole amount implicit damping iterative scheme allows the provision of displacement or stress boundary conditions in the form of the whole amount. The method will bring great convenience to solving nonlinear dynamic response of high concrete dams with complex seismic ground motion inputs and other similar problems.The drawback of the implicit damping iterative method proposed in this paper is that it requires recalculation of the stiffness matrix in each iteration step. However, as computer processing power increases, especially with the development of high performance computers, computing speed and scale will no longer be the bottleneck of scientific computing. Therefore, the shortcomings of the proposed algorithm will not be the outstanding problems of the application. In addition, the authors only carried out the numerical experiments of the incremental loading process, but if under displacement load, the implicit damping iterative algorithm can also be employed to simulate material plastic softening processes without much additional effort.



---

*Source: 486171-2014-08-31.xml*