# Semi-Implicit Algorithm for Elastoplastic Damage Models Involving Energy Integration

Authors: Ji Zhang; Jie Li
Journal: Advances in Materials Science and Engineering
 (2016)

---

## Abstract

This study aims to develop a semi-implicit constitutive integration algorithm for a class of elastoplastic damage models where calculation of damage energy release rates involves integration of free energy. The constitutive equations with energy integration are split into the elastic predictor, plastic corrector, and damage corrector. The plastic corrector is solved with an improved format of the semi-implicit spectral return mapping, which is characterized by constant flow direction and plastic moduli calculated at initial yield, enforcement of consistency at the end, and coordinate-independent formulation with an orthogonally similar stress tensor. The tangent stiffness consistent with the updating algorithm is derived. The algorithm is implemented with a recently proposed elastoplastic damage model for concrete, and several typical mechanical tests of reinforced concrete components are simulated. The present semi-implicit algorithm proves to achieve a balance between accuracy, stability, and efficiency compared with the implicit and explicit algorithms and calculate free energy accurately with small time steps.


---

## Body

## 1. Introduction

Due to its accuracy, fast convergence, and unconditional stability, the (fully) implicit backward Euler return mapping method [1] has become the most favorable one for numerically integrating inelastic constitutive models. In this method, quite large time steps can be used, hence far fewer global (finite element) iterations.However, sufficiently small time steps are warranted to accurately calculate free energy for a class of elastoplastic damage models. In these constitutive theories, a damage energy release rate is derived from a Helmholtz free energy potential. Unlike the elastic free energy (i.e., strain energy), which is determined from the current state, the plastic free energy (a part of the plastic work [2]) has to be integrated along the plastic loading history [3, 4].With very small time steps and frequent iterations, the implicit return mapping method is too expensive computationally, especially for quasi-brittle material models, which usually employ complex hardening laws and flow rules. The extremely complex partial derivatives of the plastic moduli and flow direction with respect to the stress and hardening variables are required for the process of return mapping, and the inversion of certain matrices containing these derivatives needs to be carried out to construct the consistent tangent stiffness.Return mapping can be made much simpler with the so-called semi-implicit method [5], where some quantities, including the plastic moduli and flow direction, are calculated at the beginning of a time step and remained unchanged throughout the process of return mapping, and some, including the plasticity consistency parameter, are determined at the end of the time step. For a single small time step, this approach is much easier to implement and more efficient than the implicit method; at the same time, the stress is return-mapped effectively onto the subsequent yield surface by enforcing the consistency condition, and the problem of drift in the explicit forward Euler method is overcome. Thus, the semi-implicit return mapping method is an ideal choice for the abovementioned plastic free energy integration.It is the purpose of the present study to develop a semi-implicit constitutive integration algorithm for a class of elastoplastic damage models where calculation of damage energy release rates involves integration of free energy. It is intended that, with small time steps, free energy can be integrated accurately, and the numerical algorithm remains efficient and stable.


## 2. Constitutive Framework

The class of elastoplastic damage models under discussion are mainly for quasi-brittle materials, which exhibit distinct behavior in tension and compression. The general constitutive framework is presented as follows.


### 2.1. Overall State Description

Two scalarsd
+ and d
-, the tensile and compressive damage variables, are used to describe a damage state as (1)
σ
=
1
-
d
+
σ
-
+
+
1
-
d
-
σ
-
-
,where σ is the nominal (Cauchy) stress and σ
-
± is the tensile/compressive part of the effective stress σ
-: (2)
σ
-
±
=
±
∑
i
±
σ
-
i
n
→
i
⊗
n
→
i
,in which σ
-
i are the principal effective stresses, n
→
i is the orthonormal basis along the principal directions, and x
=
(
x
+
x
)
/
2.Elastoplasticity is formulated in the effective stress space:(3)
σ
-
=
C
-
:
ε
-
ε
p
,where C
- is the effective elastic stiffness, ε is the strain, and ε
p is the plastic strain.


### 2.2. Damage Evolution: Drive Chain

The tensile/compressive Helmholtz free energyψ
± is split into its elastic part ψ
e
± and plastic part ψ
p
±, mapped into their effective values ψ
-
e
± and ψ
-
p
±, and calculated as (4)
ψ
±
=
ψ
e
±
+
ψ
p
±
=
1
-
d
±
ψ
-
e
±
+
ψ
-
p
±
=
1
-
d
±
1
2
σ
-
±
:
ε
e
+
∫
0
κ
±
Z
-
±
d
κ
±
,where ε
e is the elastic strain, Z
-
± is the tensile/compressive effective plastic energy storage rate [4], and κ
± is the tensile/compressive hardening variable.The tensile/compressive damage energy release rate is defined and calculated as(5)
Y
±
=
-
∂
ψ
±
∂
d
±
=
1
2
σ
-
±
:
ε
e
+
∫
0
κ
±
Z
-
±
d
κ
±
.The tensile/compressive damage threshold is defined as(6)
R
±
=
max
⁡
R
0
±
,
m
a
x
τ
∈
0
,
t
Y
τ
±
,where R
0
± is the initial tensile/compressive damage threshold, τ is any instant throughout the loading history, and t is the current instant.Damage evolution is formulated as a function ofR
0
± and R
±: (7)
d
±
=
G
±
R
0
±
,
R
±
.


### 2.3. Plastic Evolution: Intertwined Components

Two scalarsκ
+ and κ
-, the tensile and compressive hardening variables, are used to characterize plasticity, which includes the yield criterion, (8)
f
p
σ
-
,
κ
+
,
κ
-
=
0
,and the flow rule(9)
ε
˙
p
=
λ
˙
∂
g
p
σ
-
,
κ
+
,
κ
-
∂
σ
-
,where λ is the plasticity consistency parameter and g
p is the plastic potential, and the evolution law for κ
± is as follows:(10)
κ
˙
±
=
λ
˙
h
±
σ
-
,
κ
+
,
κ
-
,where h
± is the tensile/compressive plastic modulus.Finally, with the Karush–Kuhn-Tucker loading/unloading conditions,(11)
f
p
≤
0
,
λ
˙
≥
0
,
λ
˙
f
p
=
0
,and the consistency condition, (12)
λ
˙
f
˙
p
=
0
,the plastic relation is complete.


## 2.1. Overall State Description

Two scalarsd
+ and d
-, the tensile and compressive damage variables, are used to describe a damage state as (1)
σ
=
1
-
d
+
σ
-
+
+
1
-
d
-
σ
-
-
,where σ is the nominal (Cauchy) stress and σ
-
± is the tensile/compressive part of the effective stress σ
-: (2)
σ
-
±
=
±
∑
i
±
σ
-
i
n
→
i
⊗
n
→
i
,in which σ
-
i are the principal effective stresses, n
→
i is the orthonormal basis along the principal directions, and x
=
(
x
+
x
)
/
2.Elastoplasticity is formulated in the effective stress space:(3)
σ
-
=
C
-
:
ε
-
ε
p
,where C
- is the effective elastic stiffness, ε is the strain, and ε
p is the plastic strain.


## 2.2. Damage Evolution: Drive Chain

The tensile/compressive Helmholtz free energyψ
± is split into its elastic part ψ
e
± and plastic part ψ
p
±, mapped into their effective values ψ
-
e
± and ψ
-
p
±, and calculated as (4)
ψ
±
=
ψ
e
±
+
ψ
p
±
=
1
-
d
±
ψ
-
e
±
+
ψ
-
p
±
=
1
-
d
±
1
2
σ
-
±
:
ε
e
+
∫
0
κ
±
Z
-
±
d
κ
±
,where ε
e is the elastic strain, Z
-
± is the tensile/compressive effective plastic energy storage rate [4], and κ
± is the tensile/compressive hardening variable.The tensile/compressive damage energy release rate is defined and calculated as(5)
Y
±
=
-
∂
ψ
±
∂
d
±
=
1
2
σ
-
±
:
ε
e
+
∫
0
κ
±
Z
-
±
d
κ
±
.The tensile/compressive damage threshold is defined as(6)
R
±
=
max
⁡
R
0
±
,
m
a
x
τ
∈
0
,
t
Y
τ
±
,where R
0
± is the initial tensile/compressive damage threshold, τ is any instant throughout the loading history, and t is the current instant.Damage evolution is formulated as a function ofR
0
± and R
±: (7)
d
±
=
G
±
R
0
±
,
R
±
.


## 2.3. Plastic Evolution: Intertwined Components

Two scalarsκ
+ and κ
-, the tensile and compressive hardening variables, are used to characterize plasticity, which includes the yield criterion, (8)
f
p
σ
-
,
κ
+
,
κ
-
=
0
,and the flow rule(9)
ε
˙
p
=
λ
˙
∂
g
p
σ
-
,
κ
+
,
κ
-
∂
σ
-
,where λ is the plasticity consistency parameter and g
p is the plastic potential, and the evolution law for κ
± is as follows:(10)
κ
˙
±
=
λ
˙
h
±
σ
-
,
κ
+
,
κ
-
,where h
± is the tensile/compressive plastic modulus.Finally, with the Karush–Kuhn-Tucker loading/unloading conditions,(11)
f
p
≤
0
,
λ
˙
≥
0
,
λ
˙
f
p
=
0
,and the consistency condition, (12)
λ
˙
f
˙
p
=
0
,the plastic relation is complete.


## 3. Operator Split

In accordance with the concept and methods of product formulas and operator split [3, 6], the constitutive equations presented in the previous section are replaced here with the elastic, plastic, and damage parts in Table 1, which can serve as a mathematical foundation for developing first-order accurate algorithms.Table 1
Split of rate constitutive equations.
Total
=Elastic part
+Plastic part
+Damage part
ε
˙
=
∇
s
u
˙
→
=
u
˙
→
⊗
∇
→
+
∇
→
⊗
u
˙
→
2
ε
˙
=
∇
s
u
˙
→
ε
˙
p
=
λ
˙
∂
g
p
∂
σ
-
ε
˙
p
=
λ
˙
∂
g
p
∂
σ
-
κ
˙
±
=
λ
˙
h
±
κ
˙
±
=
λ
˙
h
±
σ
-
˙
=
C
-
:
ε
˙
-
ε
˙
p
σ
-
˙
=
C
-
:
ε
˙
σ
-
˙
=
-
C
-
:
ε
˙
p
Z
-
˙
±
=
σ
-
˙
±
∂
Z
-
±
∂
σ
-
Z
-
˙
±
=
σ
-
˙
±
e
l
a
s
t
i
c
p
a
r
t
+
σ
-
˙
±
∂
Z
-
±
∂
σ
-
ψ
-
˙
e
±
=
σ
-
±
:
ε
˙
-
ε
˙
p
ψ
-
˙
e
±
=
σ
-
±
:
ε
˙
ψ
-
˙
e
±
=
-
σ
-
±
:
ε
˙
p
ψ
-
˙
p
±
=
Z
-
±
κ
˙
±
ψ
-
˙
p
±
=
Z
-
±
κ
˙
±
Y
˙
±
=
ψ
-
˙
e
±
+
ψ
-
˙
p
±
Y
˙
±
=
ψ
-
˙
e
±
e
l
a
s
t
i
c
p
a
r
t
+
ψ
-
˙
e
±
p
l
a
s
t
i
c
p
a
r
t
+
ψ
-
˙
p
±
p
l
a
s
t
i
c
p
a
r
t
R
˙
±
=
Y
˙
±
,
Y
±
=
R
±
,
Y
˙
±
>
0
0
,
Y
±
<
R
±
o
r
Y
˙
±
≤
0
R
˙
±
=
Y
˙
±
,
Y
±
=
R
±
,
Y
˙
±
>
0
0
,
Y
±
<
R
±
o
r
Y
˙
±
≤
0
d
˙
±
=
R
˙
±
∂
G
±
∂
R
±
d
˙
±
=
R
˙
±
∂
G
±
∂
R
±
σ
˙
=
1
-
d
+
σ
-
˙
+
+
1
-
d
-
σ
-
˙
-
-
d
˙
+
σ
-
+
-
d
˙
-
σ
-
-
σ
˙
=
1
-
d
+
σ
-
˙
+
+
1
-
d
-
σ
-
˙
-
σ
˙
=
1
-
d
+
σ
-
˙
+
+
1
-
d
-
σ
-
˙
-
σ
˙
=
-
d
˙
+
σ
-
+
-
d
˙
-
σ
-
-These three parts are taken as split problems in a series: the first problem provides the initial condition for the second one, the second one for the third one, and the third one for the first one in the next time step; and all the split problems add up to the original elastoplastic damage problem. In this manner, the elastoplastic damage constitutive equations can be integrated incrementally over a sequence of time stepst
n
,
t
n
+
1
,
n
=
0,1
,
2
….


### 3.1. Elastic Predictor

The first problem, usually referred to as the elastic predictor, is simply the update of the strain(13)
ε
n
+
1
=
ε
n
+
∇
s
Δ
u
→
n
+
1and the effective stress (14)
σ
-
n
+
1
t
r
i
a
l
=
C
-
:
ε
n
+
1
-
ε
n
passuming purely elastic behavior.


### 3.2. Plastic Corrector

The second problem, the plastic corrector, is to relax the plastic strain incrementΔ
ε
n
+
1
p and direct the excursive effective stress σ
-
n
+
1
t
r
i
a
l back to the updated yield surface: (15)
σ
-
n
+
1
=
σ
-
n
+
1
t
r
i
a
l
-
C
-
:
Δ
ε
n
+
1
pif σ
-
n
+
1
t
r
i
a
l violates the yield criterion (16)
f
p
σ
-
n
+
1
t
r
i
a
l
,
κ
n
+
,
κ
n
-
>
0
.Various types of algorithms have been developed for this problem, including the explicit forward Euler method, implicit backward Euler method, and semi-implicit method, which differ from one another in their specific (iterative) updating schemes for (15). In a broader sense, all these methods can be regarded as “return mapping”; many researchers mean only those that accurately enforce consistency by this popular term.


### 3.3. Damage Corrector

The third problem, the damage corrector, is to evaluate damage evolution:(17)
σ
n
+
1
=
1
-
d
n
+
1
+
σ
-
n
+
1
+
+
1
-
d
n
+
1
-
σ
-
n
+
1
-
.Like the elastic predictor, it is also straightforward and involves no iteration. However, its formulation plays an important role in the construction of the algorithmic consistent tangent stiffness as demonstrated in the following Section 5.


## 3.1. Elastic Predictor

The first problem, usually referred to as the elastic predictor, is simply the update of the strain(13)
ε
n
+
1
=
ε
n
+
∇
s
Δ
u
→
n
+
1and the effective stress (14)
σ
-
n
+
1
t
r
i
a
l
=
C
-
:
ε
n
+
1
-
ε
n
passuming purely elastic behavior.


## 3.2. Plastic Corrector

The second problem, the plastic corrector, is to relax the plastic strain incrementΔ
ε
n
+
1
p and direct the excursive effective stress σ
-
n
+
1
t
r
i
a
l back to the updated yield surface: (15)
σ
-
n
+
1
=
σ
-
n
+
1
t
r
i
a
l
-
C
-
:
Δ
ε
n
+
1
pif σ
-
n
+
1
t
r
i
a
l violates the yield criterion (16)
f
p
σ
-
n
+
1
t
r
i
a
l
,
κ
n
+
,
κ
n
-
>
0
.Various types of algorithms have been developed for this problem, including the explicit forward Euler method, implicit backward Euler method, and semi-implicit method, which differ from one another in their specific (iterative) updating schemes for (15). In a broader sense, all these methods can be regarded as “return mapping”; many researchers mean only those that accurately enforce consistency by this popular term.


## 3.3. Damage Corrector

The third problem, the damage corrector, is to evaluate damage evolution:(17)
σ
n
+
1
=
1
-
d
n
+
1
+
σ
-
n
+
1
+
+
1
-
d
n
+
1
-
σ
-
n
+
1
-
.Like the elastic predictor, it is also straightforward and involves no iteration. However, its formulation plays an important role in the construction of the algorithmic consistent tangent stiffness as demonstrated in the following Section 5.


## 4. Semi-Implicit Return Mapping

For the plastic corrector problem, the semi-implicit method [5] is adopted and improved in the present study in order to integrate free energy accurately and efficiently with small time steps t
n
,
t
n
+
1.


### 4.1. General Format

Different from the implicit backward Euler method, the flow direction∂
g
p
/
∂
σ
- and plastic moduli h
+
,
h
- at the beginning of the elastoplastic loading range t
n
+
e are used in the whole elastoplastic range t
n
+
e
,
t
n
+
1, where e
∈
0,1 is the relative length of the purely elastic range t
n
,
t
n
+
e, so that the iterations become much easier and simpler [7]. Different from the explicit forward Euler method, the consistency condition is enforced at the end of the step t
n
+
1 so that drift no longer occurs.The set of nonlinear equations (8)–(10) are cast in a semi-implicit format of their residuals: (18)
f
n
+
1
p
=
f
p
σ
-
n
+
1
,
κ
n
+
1
+
,
κ
n
+
1
-
,
R
n
+
1
ε
p
=
-
ε
n
+
1
p
+
ε
n
p
+
Δ
λ
n
+
1
∂
g
n
+
e
p
∂
σ
-
n
+
e
R
n
+
1
κ
±
=
-
κ
n
+
1
±
+
κ
n
±
+
Δ
λ
n
+
1
h
n
+
e
±
,where (19)
g
n
+
e
p
=
g
p
σ
-
n
+
e
,
κ
n
+
,
κ
n
-
,
h
n
+
e
±
=
h
±
σ
-
n
+
e
,
κ
n
+
,
κ
n
-
.To solve (18), the Newton-Raphson iterations (20)
f
n
+
1
p
k
+
∂
f
n
+
1
p
k
∂
σ
-
n
+
1
k
:
δ
σ
-
n
+
1
k
+
∂
f
n
+
1
p
k
∂
κ
n
+
1
+
k
δ
κ
n
+
1
+
k
+
∂
f
n
+
1
p
k
∂
κ
n
+
1
-
k
δ
κ
n
+
1
-
k
=
0
,
R
n
+
1
k
ε
p
+
C
-
-
1
:
δ
σ
-
n
+
1
k
+
∂
g
n
+
e
p
∂
σ
-
n
+
e
δ
Δ
λ
n
+
1
k
=
0
,
R
n
+
1
k
κ
±
-
δ
κ
n
+
1
±
k
+
h
n
+
e
±
δ
Δ
λ
n
+
1
k
=
0are used, where the superscript (k) denotes the kth iteration within a time step, leading to (21)
Δ
λ
n
+
1
k
+
1
=
Δ
λ
n
+
1
k
+
δ
Δ
λ
n
+
1
k
,
κ
n
+
1
±
k
+
1
=
κ
n
+
1
±
k
+
δ
κ
n
+
1
±
k
,
σ
-
n
+
1
k
+
1
=
σ
-
n
+
1
k
+
δ
σ
-
n
+
1
k
,
ε
n
+
1
p
k
+
1
=
ε
n
+
1
p
k
-
C
-
-
1
:
δ
σ
-
n
+
1
k
,where δ
Δ
λ
n
+
1
k, δ
κ
n
+
1
±
k, and δ
σ
-
n
+
1
k are obtained by solving (20).Remark 1.
It is noted that whene
=
0, that is, the stress state is on the current yield surface at the beginning of the time step t
n, the above procedure reduces to the existing semi-implicit format [5].


### 4.2. Spectral Version

The plastic constitutive equations for quasi-brittle materials often involve terms of principal stresses. As a consequence, eigenvalue calculation needs to be performed for each iteration (k) in the return mapping process developed in the previous subsection, which substantially increases the computational cost. For these cases, the spectral version [8, 9] of the return mapping methods is much more efficient. It is adopted and formulated rigorously in the present subsection.The matrixσ
-
n
+
1
t
r
i
a
l corresponding to the tensor σ
-
n
+
1
t
r
i
a
l in the structural coordinate system is spectrally decomposed as (22)
σ
-
n
+
1
t
r
i
a
l
=
Q
σ
-
^
n
+
1
t
r
i
a
l
Q
T
,where Q
=
n
→
1
n
→
2
n
→
3 and σ
-
^
n
+
1
t
r
i
a
l
=
d
i
a
g
σ
-
1
t
r
i
a
l
,
σ
-
2
t
r
i
a
l
,
σ
-
3
t
r
i
a
l, which corresponds to σ
-
n
+
1
t
r
i
a
l in the principal coordinate system.In this study, a new tensorσ
-
^
n
+
1
t
r
i
a
l corresponding to the diagonal matrix σ
-
^
n
+
1
t
r
i
a
l in the structural coordinate system is defined, and (22) is expressed in a tensor form: (23)
σ
-
n
+
1
t
r
i
a
l
=
Q
·
σ
-
^
n
+
1
t
r
i
a
l
·
Q
T
,where Q is the tensor corresponding to Q in the structural coordinate system. Because Q is an orthogonal matrix, Q is an orthogonal tensor. By definition, the newly introduced tensor σ
-
^
n
+
1
t
r
i
a
l and the traditional tensor σ
-
n
+
1
t
r
i
a
l are orthogonally similar.It has been proved [8, 9] that, for an isotropic material with coaxial plastic flow, that is, a plastic potential g
p exists such that (9) holds, σ
-
n
+
1 after return mapping in (15) has the identical eigenspace with σ
-
n
+
1
t
r
i
a
l. Thus, σ
-
n
+
1 corresponding to σ
-
n
+
1 in the structural coordinate system can be spectrally decomposed with the same eigenvector matrix Q as (24)
σ
-
n
+
1
=
Q
σ
-
^
n
+
1
Q
T
,where σ
-
^
n
+
1
=
d
i
a
g
σ
-
1
,
σ
-
2
,
σ
-
3, which corresponds to σ
-
n
+
1 in the principal coordinate system.In a similar way,σ
-
^
n
+
1 corresponding to σ
-
^
n
+
1 in the structural coordinate system is defined, and (24) is expressed in a tensor form: (25)
σ
-
n
+
1
=
Q
·
σ
-
^
n
+
1
·
Q
T
,which indicates that σ
-
^
n
+
1 and σ
-
n
+
1 are orthogonally similar.The return mapping procedure expressed by (18)–(21) is carried out here on the orthogonally similar tensor σ
-
^ as (26)
f
^
n
+
1
p
=
f
^
p
σ
-
^
n
+
1
,
κ
n
+
1
+
,
κ
n
+
1
-
=
f
p
σ
-
n
+
1
,
κ
n
+
1
+
,
κ
n
+
1
-
,
R
n
+
1
ε
p
=
-
ε
n
+
1
p
+
ε
n
p
+
Δ
λ
n
+
1
∂
g
^
n
+
e
p
∂
σ
-
^
n
+
e
,
R
n
+
1
κ
±
=
-
κ
n
+
1
±
+
κ
n
±
+
Δ
λ
n
+
1
h
^
n
+
e
±
,
(27)
g
^
n
+
e
p
=
g
^
p
σ
-
^
n
+
e
,
κ
n
+
,
κ
n
-
=
g
p
σ
-
n
+
e
,
κ
n
+
,
κ
n
-
,
h
^
n
+
e
±
=
h
^
±
σ
-
^
n
+
e
,
κ
n
+
,
κ
n
-
=
h
±
σ
-
n
+
e
,
κ
n
+
,
κ
n
-
,
(28)
f
^
n
+
1
p
k
+
∂
f
^
n
+
1
p
k
∂
σ
-
^
n
+
1
k
:
δ
σ
-
^
n
+
1
k
+
∂
f
^
n
+
1
p
k
∂
κ
n
+
1
+
k
δ
κ
n
+
1
+
k
+
∂
f
^
n
+
1
p
k
∂
κ
n
+
1
-
k
δ
κ
n
+
1
-
k
=
0
,
R
n
+
1
k
ε
p
+
C
-
-
1
:
δ
σ
-
^
n
+
1
k
+
∂
g
^
n
+
e
p
∂
σ
-
^
n
+
e
δ
Δ
λ
n
+
1
k
=
0
,
R
n
+
1
k
κ
±
-
δ
κ
n
+
1
±
k
+
h
^
n
+
e
±
δ
Δ
λ
n
+
1
k
=
0
,
(29)
Δ
λ
n
+
1
k
+
1
=
Δ
λ
n
+
1
k
+
δ
Δ
λ
n
+
1
k
,
κ
n
+
1
±
k
+
1
=
κ
n
+
1
±
k
+
δ
κ
n
+
1
±
k
,
σ
-
^
n
+
1
k
+
1
=
σ
-
^
n
+
1
k
+
δ
σ
-
^
n
+
1
k
,
ε
n
+
1
p
k
+
1
=
ε
n
+
1
p
k
-
C
-
-
1
:
δ
σ
-
^
n
+
1
k
,where eigenvalue calculation needs to be performed only for the first iteration (k
=
1), and then σ
-
i are iterated upon with n
→
i fixed. The existence of the functions f
^
p
,
g
^
p, and h
^
± in (26) and (27) is based on the assumption that the material is isotropic [10].Remark 2.
It can be seen from (28) and (29) that the introduction of the orthogonally similar tensor σ
-
^ is necessary for the coordinate-independent implementation of the double dot operation between higher-order tensors.


## 4.1. General Format

Different from the implicit backward Euler method, the flow direction∂
g
p
/
∂
σ
- and plastic moduli h
+
,
h
- at the beginning of the elastoplastic loading range t
n
+
e are used in the whole elastoplastic range t
n
+
e
,
t
n
+
1, where e
∈
0,1 is the relative length of the purely elastic range t
n
,
t
n
+
e, so that the iterations become much easier and simpler [7]. Different from the explicit forward Euler method, the consistency condition is enforced at the end of the step t
n
+
1 so that drift no longer occurs.The set of nonlinear equations (8)–(10) are cast in a semi-implicit format of their residuals: (18)
f
n
+
1
p
=
f
p
σ
-
n
+
1
,
κ
n
+
1
+
,
κ
n
+
1
-
,
R
n
+
1
ε
p
=
-
ε
n
+
1
p
+
ε
n
p
+
Δ
λ
n
+
1
∂
g
n
+
e
p
∂
σ
-
n
+
e
R
n
+
1
κ
±
=
-
κ
n
+
1
±
+
κ
n
±
+
Δ
λ
n
+
1
h
n
+
e
±
,where (19)
g
n
+
e
p
=
g
p
σ
-
n
+
e
,
κ
n
+
,
κ
n
-
,
h
n
+
e
±
=
h
±
σ
-
n
+
e
,
κ
n
+
,
κ
n
-
.To solve (18), the Newton-Raphson iterations (20)
f
n
+
1
p
k
+
∂
f
n
+
1
p
k
∂
σ
-
n
+
1
k
:
δ
σ
-
n
+
1
k
+
∂
f
n
+
1
p
k
∂
κ
n
+
1
+
k
δ
κ
n
+
1
+
k
+
∂
f
n
+
1
p
k
∂
κ
n
+
1
-
k
δ
κ
n
+
1
-
k
=
0
,
R
n
+
1
k
ε
p
+
C
-
-
1
:
δ
σ
-
n
+
1
k
+
∂
g
n
+
e
p
∂
σ
-
n
+
e
δ
Δ
λ
n
+
1
k
=
0
,
R
n
+
1
k
κ
±
-
δ
κ
n
+
1
±
k
+
h
n
+
e
±
δ
Δ
λ
n
+
1
k
=
0are used, where the superscript (k) denotes the kth iteration within a time step, leading to (21)
Δ
λ
n
+
1
k
+
1
=
Δ
λ
n
+
1
k
+
δ
Δ
λ
n
+
1
k
,
κ
n
+
1
±
k
+
1
=
κ
n
+
1
±
k
+
δ
κ
n
+
1
±
k
,
σ
-
n
+
1
k
+
1
=
σ
-
n
+
1
k
+
δ
σ
-
n
+
1
k
,
ε
n
+
1
p
k
+
1
=
ε
n
+
1
p
k
-
C
-
-
1
:
δ
σ
-
n
+
1
k
,where δ
Δ
λ
n
+
1
k, δ
κ
n
+
1
±
k, and δ
σ
-
n
+
1
k are obtained by solving (20).Remark 1.
It is noted that whene
=
0, that is, the stress state is on the current yield surface at the beginning of the time step t
n, the above procedure reduces to the existing semi-implicit format [5].


## 4.2. Spectral Version

The plastic constitutive equations for quasi-brittle materials often involve terms of principal stresses. As a consequence, eigenvalue calculation needs to be performed for each iteration (k) in the return mapping process developed in the previous subsection, which substantially increases the computational cost. For these cases, the spectral version [8, 9] of the return mapping methods is much more efficient. It is adopted and formulated rigorously in the present subsection.The matrixσ
-
n
+
1
t
r
i
a
l corresponding to the tensor σ
-
n
+
1
t
r
i
a
l in the structural coordinate system is spectrally decomposed as (22)
σ
-
n
+
1
t
r
i
a
l
=
Q
σ
-
^
n
+
1
t
r
i
a
l
Q
T
,where Q
=
n
→
1
n
→
2
n
→
3 and σ
-
^
n
+
1
t
r
i
a
l
=
d
i
a
g
σ
-
1
t
r
i
a
l
,
σ
-
2
t
r
i
a
l
,
σ
-
3
t
r
i
a
l, which corresponds to σ
-
n
+
1
t
r
i
a
l in the principal coordinate system.In this study, a new tensorσ
-
^
n
+
1
t
r
i
a
l corresponding to the diagonal matrix σ
-
^
n
+
1
t
r
i
a
l in the structural coordinate system is defined, and (22) is expressed in a tensor form: (23)
σ
-
n
+
1
t
r
i
a
l
=
Q
·
σ
-
^
n
+
1
t
r
i
a
l
·
Q
T
,where Q is the tensor corresponding to Q in the structural coordinate system. Because Q is an orthogonal matrix, Q is an orthogonal tensor. By definition, the newly introduced tensor σ
-
^
n
+
1
t
r
i
a
l and the traditional tensor σ
-
n
+
1
t
r
i
a
l are orthogonally similar.It has been proved [8, 9] that, for an isotropic material with coaxial plastic flow, that is, a plastic potential g
p exists such that (9) holds, σ
-
n
+
1 after return mapping in (15) has the identical eigenspace with σ
-
n
+
1
t
r
i
a
l. Thus, σ
-
n
+
1 corresponding to σ
-
n
+
1 in the structural coordinate system can be spectrally decomposed with the same eigenvector matrix Q as (24)
σ
-
n
+
1
=
Q
σ
-
^
n
+
1
Q
T
,where σ
-
^
n
+
1
=
d
i
a
g
σ
-
1
,
σ
-
2
,
σ
-
3, which corresponds to σ
-
n
+
1 in the principal coordinate system.In a similar way,σ
-
^
n
+
1 corresponding to σ
-
^
n
+
1 in the structural coordinate system is defined, and (24) is expressed in a tensor form: (25)
σ
-
n
+
1
=
Q
·
σ
-
^
n
+
1
·
Q
T
,which indicates that σ
-
^
n
+
1 and σ
-
n
+
1 are orthogonally similar.The return mapping procedure expressed by (18)–(21) is carried out here on the orthogonally similar tensor σ
-
^ as (26)
f
^
n
+
1
p
=
f
^
p
σ
-
^
n
+
1
,
κ
n
+
1
+
,
κ
n
+
1
-
=
f
p
σ
-
n
+
1
,
κ
n
+
1
+
,
κ
n
+
1
-
,
R
n
+
1
ε
p
=
-
ε
n
+
1
p
+
ε
n
p
+
Δ
λ
n
+
1
∂
g
^
n
+
e
p
∂
σ
-
^
n
+
e
,
R
n
+
1
κ
±
=
-
κ
n
+
1
±
+
κ
n
±
+
Δ
λ
n
+
1
h
^
n
+
e
±
,
(27)
g
^
n
+
e
p
=
g
^
p
σ
-
^
n
+
e
,
κ
n
+
,
κ
n
-
=
g
p
σ
-
n
+
e
,
κ
n
+
,
κ
n
-
,
h
^
n
+
e
±
=
h
^
±
σ
-
^
n
+
e
,
κ
n
+
,
κ
n
-
=
h
±
σ
-
n
+
e
,
κ
n
+
,
κ
n
-
,
(28)
f
^
n
+
1
p
k
+
∂
f
^
n
+
1
p
k
∂
σ
-
^
n
+
1
k
:
δ
σ
-
^
n
+
1
k
+
∂
f
^
n
+
1
p
k
∂
κ
n
+
1
+
k
δ
κ
n
+
1
+
k
+
∂
f
^
n
+
1
p
k
∂
κ
n
+
1
-
k
δ
κ
n
+
1
-
k
=
0
,
R
n
+
1
k
ε
p
+
C
-
-
1
:
δ
σ
-
^
n
+
1
k
+
∂
g
^
n
+
e
p
∂
σ
-
^
n
+
e
δ
Δ
λ
n
+
1
k
=
0
,
R
n
+
1
k
κ
±
-
δ
κ
n
+
1
±
k
+
h
^
n
+
e
±
δ
Δ
λ
n
+
1
k
=
0
,
(29)
Δ
λ
n
+
1
k
+
1
=
Δ
λ
n
+
1
k
+
δ
Δ
λ
n
+
1
k
,
κ
n
+
1
±
k
+
1
=
κ
n
+
1
±
k
+
δ
κ
n
+
1
±
k
,
σ
-
^
n
+
1
k
+
1
=
σ
-
^
n
+
1
k
+
δ
σ
-
^
n
+
1
k
,
ε
n
+
1
p
k
+
1
=
ε
n
+
1
p
k
-
C
-
-
1
:
δ
σ
-
^
n
+
1
k
,where eigenvalue calculation needs to be performed only for the first iteration (k
=
1), and then σ
-
i are iterated upon with n
→
i fixed. The existence of the functions f
^
p
,
g
^
p, and h
^
± in (26) and (27) is based on the assumption that the material is isotropic [10].Remark 2.
It can be seen from (28) and (29) that the introduction of the orthogonally similar tensor σ
-
^ is necessary for the coordinate-independent implementation of the double dot operation between higher-order tensors.


## 5. Consistent Tangent Stiffness

The differentiation with respect to the strain incrementΔ
ε
n
+
1 of the resulting σ
-
n
+
1 and σ
n
+
1 gives, respectively, the tangent stiffness consistent with the updating algorithms in elastic-plastic parts and elastic-plastic-damage parts: (30)
C
-
n
+
1
t
a
n
=
∂
σ
-
n
+
1
∂
ε
n
+
1
,
C
n
+
1
t
a
n
=
∂
σ
n
+
1
∂
ε
n
+
1which lead to (31)
C
n
+
1
t
a
n
=
∂
σ
n
+
1
∂
σ
-
n
+
1
:
C
-
n
+
1
t
a
n
.On the basis of the continuum elastoplastic tangent stiffness expression(32)
C
-
t
a
n
=
C
-
-
∂
g
p
/
∂
σ
-
:
C
-
⊗
C
-
:
∂
f
p
/
∂
σ
-
∂
g
p
/
∂
σ
-
:
C
-
:
∂
f
p
/
∂
σ
-
-
∂
f
p
/
∂
κ
+
h
+
-
∂
f
p
/
∂
κ
-
h
-replacing the flow direction and plastic moduli with their values at t
n
+
e and the other quantities at t
n
+
1, one obtains(33)
C
-
n
+
1
t
a
n
=
C
-
-
∂
g
n
+
e
p
/
∂
σ
-
n
+
e
:
C
-
⊗
C
-
:
∂
f
n
+
1
p
/
∂
σ
-
n
+
1
∂
g
n
+
e
p
/
∂
σ
-
n
+
e
:
C
-
:
∂
f
n
+
1
p
/
∂
σ
-
n
+
1
-
∂
f
n
+
1
p
/
∂
κ
n
+
1
+
h
n
+
e
+
-
∂
f
n
+
1
p
/
∂
κ
n
+
1
-
h
n
+
e
-
.∂
σ
n
+
1
/
∂
σ
-
n
+
1 in (31) needs to be derived in accordance with specific damage formulation. From (1), one can obtain (34)
∂
σ
∂
σ
-
=
1
-
d
+
∂
σ
-
+
∂
σ
-
+
1
-
d
-
∂
σ
-
-
∂
σ
-
-
σ
-
+
⊗
∂
d
+
∂
σ
-
-
σ
-
-
⊗
∂
d
-
∂
σ
-
,where (35)
∂
σ
-
+
∂
σ
-
=
∑
i
H
σ
-
i
n
→
i
⊗
n
→
i
⊗
n
→
i
⊗
n
→
i
+
2
∑
i
<
j
σ
-
i
-
σ
-
j
σ
-
i
-
σ
-
j
sym
⁡
n
→
i
⊗
n
→
j
⊗
s
y
m
n
→
i
⊗
n
→
j
,
∂
σ
-
-
∂
σ
-
=
∑
i
H
-
σ
-
i
n
→
i
⊗
n
→
i
⊗
n
→
i
⊗
n
→
i
-
2
∑
i
<
j
-
σ
-
i
-
-
σ
-
j
σ
-
i
-
σ
-
j
sym
⁡
n
→
i
⊗
n
→
j
⊗
sym
⁡
n
→
i
⊗
n
→
jas derived in [11], and (36)
∂
d
±
∂
σ
-
=
∂
d
±
∂
R
±
∂
R
±
∂
σ
-
.In (36), ∂
d
±
/
∂
R
± depend on damage evolution functions and ∂
R
±
/
∂
σ
- on definitions of Y
±. If the formulation of Z
± in [4] is adopted, the specific expressions read (37)
∂
R
+
∂
σ
-
=
1
2
σ
-
:
C
-
-
1
:
∂
σ
-
+
∂
σ
-
+
σ
-
+
:
C
-
-
1
+
1
w
Δ
κ
+
H
σ
-
m
a
x
n
→
m
a
x
⊗
n
→
m
a
x
,
∂
R
-
∂
σ
-
=
1
2
σ
-
:
C
-
-
1
:
∂
σ
-
-
∂
σ
-
+
σ
-
-
:
C
-
-
1
-
1
1
-
w
Δ
κ
-
H
-
σ
-
m
i
n
n
→
m
i
n
⊗
n
→
m
i
n
,where (38)
w
=
∑
i
σ
-
i
∑
i
σ
-
i
.


## 6. Implementation

The complete algorithm of the present study is summarized as follows.Summary of integration algorithm is as follows:(39)
ε
n
+
1
=
ε
n
+
∇
s
Δ
u
→
n
+
1
I
f
f
p
σ
-
n
+
1
t
r
i
a
l
,
κ
n
+
,
κ
n
-
>
T
O
L
Δ
λ
n
+
1
k
+
1
=
Δ
λ
n
+
1
k
+
δ
Δ
λ
n
+
1
k
ε
n
+
1
p
k
+
1
=
ε
n
p
+
Δ
λ
n
+
1
k
+
1
∂
g
^
n
+
e
p
∂
σ
-
^
n
+
e
κ
n
+
1
±
k
+
1
=
κ
n
±
+
Δ
λ
n
+
1
k
+
1
h
^
n
+
e
±
σ
-
n
+
1
t
r
i
a
l
=
C
-
:
ε
n
+
1
-
ε
n
p
σ
-
^
n
+
1
k
+
1
=
σ
-
^
n
+
1
t
r
i
a
l
-
Δ
λ
n
+
1
k
+
1
C
-
:
∂
g
^
n
+
e
p
∂
σ
-
^
n
+
e
I
f
f
^
p
σ
-
^
n
+
1
k
+
1
,
κ
n
+
1
+
k
+
1
,
κ
n
+
1
-
k
+
1
≤
T
O
L
ψ
-
n
+
1
e
±
=
1
2
σ
-
n
+
1
±
:
ε
n
+
1
-
ε
n
+
1
p
ψ
-
n
+
1
p
±
=
ψ
-
n
p
±
+
1
2
Z
-
n
±
+
Z
-
n
+
1
±
Δ
κ
n
+
1
±
Y
n
+
1
±
=
ψ
-
n
+
1
e
±
+
ψ
-
n
+
1
p
±
R
n
+
1
±
=
max
⁡
R
n
±
,
Y
n
+
1
±
d
n
+
1
±
=
G
±
R
0
±
,
R
n
+
1
±
σ
n
+
1
=
1
-
d
n
+
1
+
σ
-
n
+
1
+
+
1
-
d
n
+
1
-
σ
-
n
+
1
-
C
n
+
1
t
a
n
=
1
-
d
n
+
1
+
∂
σ
-
n
+
1
+
∂
σ
-
n
+
1
+
1
-
d
n
+
1
-
∂
σ
-
n
+
1
-
∂
σ
-
n
+
1
-
∂
d
n
+
1
+
∂
R
n
+
1
+
σ
-
n
+
1
+
⊗
∂
R
n
+
1
+
∂
σ
-
n
+
1
-
∂
d
n
+
1
-
∂
R
n
+
1
-
σ
-
n
+
1
-
⊗
∂
R
n
+
1
-
∂
σ
-
n
+
1
:
C
-
-
∂
g
n
+
e
p
/
∂
σ
-
n
+
e
:
C
-
⊗
C
-
:
∂
f
n
+
1
p
/
∂
σ
-
n
+
1
∂
g
n
+
e
p
/
∂
σ
-
n
+
e
:
C
-
:
∂
f
n
+
1
p
/
∂
σ
-
n
+
1
-
∂
f
n
+
1
p
/
∂
κ
n
+
1
+
h
n
+
e
+
-
∂
f
n
+
1
p
/
∂
κ
n
+
1
-
h
n
+
e
-
.As numerical examples, this algorithm is implemented with the specific elastoplastic damage model for concrete presented in [4]. Several typical mechanical tests of reinforced concrete components are simulated with the commercial finite element program Abaqus and its user material subroutine UMAT. Solid elements are used for concrete and truss elements for steel rebars, between which perfect bonding is assumed.The uniaxial mechanical properties for concrete are first calibrated to fit the stress-strain curves of uniaxial tension and compression recommended in the Chinese code [12]. The tension stiffening effects due to the presence of reinforcement are accounted for by using the descending stress-strain curve by Stevens et al. [13] and the confining effects resulting from closely spaced stirrups by fitting the stress-strain curve for core concrete by Kent and Park [14] and Scott et al. [15]. Besides, the common values for the equibiaxial strength ratio f
b
0
-
/
f
0
-
=
1.2, Poisson’s ratio ν
=
0.2, and dilation parameter α
p
=
0.2 are set for multiaxial behavior. For reinforcing steel, uniaxial elastoplastic models are employed.


### 6.1. Slab

For application in reinforced concrete slabs, the test by McNeice [16] (Figure 1(a)) is simulated. The two-way slab has a bottom reinforcement ratio of 0.85% in each direction. In the simulation, Young’s modulus E
s
=
200
×
1
0
3
N
/
m
m
2 and yield point f
s
=
345
N
/
m
m
2 are set for steel, and effective Young’s modulus E
-
=
35.1
×
1
0
3
N
/
m
m
2, initial damage stresses in uniaxial tension and compression f
0
+
=
3.17
N
/
m
m
2, f
0
-
=
21.0
N
/
m
m
2, and initial effective yield stresses in uniaxial tension and compression f
-
y
+
=
3.17
N
/
m
m
2, f
-
y
-
=
37.92
N
/
m
m
2 for concrete. The calculated load-deflection curve is shown along with the test data in Figure 1(b).Figure 1
Reinforced concrete slab under out-of-plane loading [16]. (a) Setup (dimensions in mm). (b) Results.
(a)
(b)For comparison, the implicit backward Euler return mapping algorithm is also implemented, and the calculated load-deflection curve is precisely the same as in Figure1(b). The times of calculation consumed with the semi-implicit and implicit methods are 179 s and 804 s, respectively, at 2000 displacement steps.


### 6.2. Wall

For application in reinforced concrete walls, the tests by Lefas et al. [17] (Figure 2(a)) are simulated. In one of the two selected tests (SW22 and SW24), a constant vertical load is held while the horizontal load increases. These loads are transferred through the top spreader beams, which are cast integrally with the walls. The vertical and horizontal reinforcement ratios in the walls are 2.49% and 0.82%, respectively, and the ratio of closed stirrups to confined concrete of 140 mm-wide concealed columns at both edges is 0.52%. In the simulations, E
s
=
210
×
1
0
3
N
/
m
m
2, f
s
,
v
=
470
N
/
m
m
2 (vertical rebars), f
s
,
h
=
520
N
/
m
m
2 (horizontal rebars), and f
s
,
s
=
420
N
/
m
m
2 (stirrups) are set for steel, E
-
=
35.4
×
1
0
3
N
/
m
m
2, f
0
+
=
3.33
N
/
m
m
2, f
0
-
=
21.3
N
/
m
m
2, f
-
y
+
=
3.33
N
/
m
m
2, and f
-
y
-
=
38.5
N
/
m
m
2 for concrete of SW24, and E
-
=
36.7
×
1
0
3
N
/
m
m
2, f
0
+
=
3.42
N
/
m
m
2, f
0
-
=
22.0
N
/
m
m
2, f
-
y
+
=
3.42
N
/
m
m
2, and f
-
y
-
=
40.8
N
/
m
m
2 for concrete of SW22. The calculated load-deflection curves are displayed along with the test data in Figures 2(b) and 2(c), respectively, where f
c
u is the cubic strength of concrete.Figure 2
Reinforced concrete walls under in-plane loading [17]. (a) Setup (dimensions in mm). (b) Results for SW24 (N
=
0, f
c
u
=
48.3N/mm2). (c) Results for SW22 (N
=
182 kN, f
c
u
=
50.6 N/mm2).
(a)
(b)
(c)For comparison, the implicit return mapping algorithm is also implemented, and the calculated load-deflection curves are precisely the same as in Figures2(b) and 2(c). The times of calculation consumed with the semi-implicit and implicit methods are 495 s/508 s and 1870 s/1922 s, respectively, at 2000 displacement steps.


## 6.1. Slab

For application in reinforced concrete slabs, the test by McNeice [16] (Figure 1(a)) is simulated. The two-way slab has a bottom reinforcement ratio of 0.85% in each direction. In the simulation, Young’s modulus E
s
=
200
×
1
0
3
N
/
m
m
2 and yield point f
s
=
345
N
/
m
m
2 are set for steel, and effective Young’s modulus E
-
=
35.1
×
1
0
3
N
/
m
m
2, initial damage stresses in uniaxial tension and compression f
0
+
=
3.17
N
/
m
m
2, f
0
-
=
21.0
N
/
m
m
2, and initial effective yield stresses in uniaxial tension and compression f
-
y
+
=
3.17
N
/
m
m
2, f
-
y
-
=
37.92
N
/
m
m
2 for concrete. The calculated load-deflection curve is shown along with the test data in Figure 1(b).Figure 1
Reinforced concrete slab under out-of-plane loading [16]. (a) Setup (dimensions in mm). (b) Results.
(a)
(b)For comparison, the implicit backward Euler return mapping algorithm is also implemented, and the calculated load-deflection curve is precisely the same as in Figure1(b). The times of calculation consumed with the semi-implicit and implicit methods are 179 s and 804 s, respectively, at 2000 displacement steps.


## 6.2. Wall

For application in reinforced concrete walls, the tests by Lefas et al. [17] (Figure 2(a)) are simulated. In one of the two selected tests (SW22 and SW24), a constant vertical load is held while the horizontal load increases. These loads are transferred through the top spreader beams, which are cast integrally with the walls. The vertical and horizontal reinforcement ratios in the walls are 2.49% and 0.82%, respectively, and the ratio of closed stirrups to confined concrete of 140 mm-wide concealed columns at both edges is 0.52%. In the simulations, E
s
=
210
×
1
0
3
N
/
m
m
2, f
s
,
v
=
470
N
/
m
m
2 (vertical rebars), f
s
,
h
=
520
N
/
m
m
2 (horizontal rebars), and f
s
,
s
=
420
N
/
m
m
2 (stirrups) are set for steel, E
-
=
35.4
×
1
0
3
N
/
m
m
2, f
0
+
=
3.33
N
/
m
m
2, f
0
-
=
21.3
N
/
m
m
2, f
-
y
+
=
3.33
N
/
m
m
2, and f
-
y
-
=
38.5
N
/
m
m
2 for concrete of SW24, and E
-
=
36.7
×
1
0
3
N
/
m
m
2, f
0
+
=
3.42
N
/
m
m
2, f
0
-
=
22.0
N
/
m
m
2, f
-
y
+
=
3.42
N
/
m
m
2, and f
-
y
-
=
40.8
N
/
m
m
2 for concrete of SW22. The calculated load-deflection curves are displayed along with the test data in Figures 2(b) and 2(c), respectively, where f
c
u is the cubic strength of concrete.Figure 2
Reinforced concrete walls under in-plane loading [17]. (a) Setup (dimensions in mm). (b) Results for SW24 (N
=
0, f
c
u
=
48.3N/mm2). (c) Results for SW22 (N
=
182 kN, f
c
u
=
50.6 N/mm2).
(a)
(b)
(c)For comparison, the implicit return mapping algorithm is also implemented, and the calculated load-deflection curves are precisely the same as in Figures2(b) and 2(c). The times of calculation consumed with the semi-implicit and implicit methods are 495 s/508 s and 1870 s/1922 s, respectively, at 2000 displacement steps.


## 7. Conclusions

History-dependent energy calculation makes the quality of stability in large time steps of the popular implicit backward Euler return mapping algorithm no longer attractive for some inelastic constitutive models. More flexible and efficient integrating schemes are favored.The semi-implicit constitutive integration algorithm developed for a class of elastoplastic damage models that involve integration of free energy in the present study achieves a balance between accuracy, stability, and efficiency. In comparison with the implicit backward Euler algorithm, it is by far less challenging mathematically and less time-consuming computationally. While it is less stable, the problem can be fixed with small time steps, which are necessary for free energy integration. Another drawback is the need to find the stress where the yield surface is first reached within a time step. In comparison with the explicit forward Euler algorithm, it is much more accurate and stable at the cost of iterations for consistency.The splitting of operators and the elastic-plastic-damage updating algorithm on its basis presented in this paper are particularly formulated for damage energy release rates that are derived from elastoplastic free energy potentials. Nevertheless, the basic idea and method are also applicable to other types of constitutive models that include history-dependent calculation of certain quantities.



---

*Source: 5289642-2016-03-23.xml*