# Numerical Implementation of Spatial Elastoplastic Damage Model of Concrete in the Framework of Isogeometric Analysis Approach

Authors: Cheng Ma; Wei-zhen Chen; Jian-yuan Sun
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

This paper is a study of the numerical implementation of the spatial elastoplastic damage model of concrete by isogeometric analysis (IGA) method from three perspectives: the geometric modeling and the numerical formulation via IGA method, the constitutive model of concrete, and the solution algorithms for the local and global problems. The plasticity of concrete is considered on the basis of a nonassociated flow rule, where a three-parameter Barcelona yield surface and a modified Drucker-Prager plastic potential are used. The damage evolution of concrete driven by the internal variables is expressed by a piecewise function. In the study, the return-mapping algorithm and the substepping strategy are used for stress updating, and a new dissipation-based arc-length method with constraint path that considers the combined contribution of plasticity and damage to the energy dissipation is employed to trace the equilibrium path. After comparisons between simulation results and experimental data, the use of the elastoplastic damage model in the framework of IGA approach is proven to be practical in reflecting material properties of concrete.


---

## Body

## 1. Introduction

Armed with the progress of computer science and the in-depth study of computational theory, the numerical simulation technology of concrete structures has been developing towards highly advanced analysis that requires comprehensive capture of mechanical responses, especially the postcrack performances in three-dimensional (3D) space.In the traditional finite element modeling, computer-aided design (CAD) and computer-aided engineering (CAE) are constructed on two different platforms based on one-way information transmission which is not only complex but also time-consuming. Therefore, Hughes et al. [1, 2] integrated CAD and CAE and proposed the philosophy of IGA, which directly uses the geometric languages in CAD instead of the commonly used Lagrange functions to construct the analysis domain.The modeling of concrete mechanism is a crux in the advanced analysis due to high material nonlinearity. In continuum mechanics, concrete behavior can be simulated by elastic, elastic-plastic, or elastoplastic damage models. But the elastic damage model [3, 4] cannot consider the unrecoverable deformation of concrete, and the elastic-plastic model [5, 6] ignores the influence of damage. In comparison, the elastoplastic damage model, a combination of the first two models, is more reasonable because it can truly reflect concrete properties and therefore prevails in ductile failure [6–8], strong coupling [9, 10], or effective stress approaches [3, 11–16], among which the effective stress approach is the most widely used.The solution algorithm is also a highly nonlinear issue that can be separated into the solutions of local and global problems to be solved at the level of Gaussian point as well as the level of structure, respectively. The difficulty in solving local problem exists in updating concrete stress. To facilitate convergence, Pérez-Foguet et al. [17] presented substepping algorithm combined with return-mapping algorithm. For global problem, the arc-length method is better than the traditional force control or displacement control methods. The classical arc-length method was proposed by Riks [18] and developed by Ramm [19] and Crisfield [20]. But the method was poor in solving material nonlinear problems. Gutiérrez [21], Verhoose et al. [22], and van der Meer et al. [23] established the constraint path of arc-length on the amount of the released energy and proved convergent capability in brittle and ductile analysis.At present, the use of IGA has been widely explored in geometric modeling and mechanical engineering, but the extension use of it to structural analysis is rare. In this study, the knowledge for the numerical implementation of spatial elastoplastic damage model of concrete via IGA is introduced in detail from three aspects, that is, the geometric modeling and the numerical formulation of IGA, the constitutive model of concrete, and the solution algorithms for the local and global problems.


## 2. Basic Theory of IGA Solid Modeling

### 2.1.B-Spline Basis Function and NURBS Basis Function

The basic analyzed object of IGA [24] in 1D, 2D, and 3D space is the NURBS (Nonuniform Rational B-Spline), NURBS surface, and NURBS solid, respectively. A pth-order B-spline C
B
ξ can be constructed by n
B-spline basis functions N
i
,
p
ξ defined on the knot sequence Θ
=
ξ
1
,
ξ
2
,
…
,
ξ
n
+
p
+
1 and their corresponding control points A
i:(1)CBξ=∑i=1nNi,pξAi.The coordinates of knots should be monotonically nondecreasing (ξ
1
≤
⋯
≤
ξ
i
≤
ξ
i
+
1
≤
⋯
≤
ξ
n
+
p
+
1), with two nonoverlapped neighboring knots constructing an element interval.Forp
=
0, N
i
,
0
ξ is defined as(2)Ni,0ξ=1ifξi≤ξ≤ξi+10otherwise.Forp
≥
1, N
i
,
p
ξ can be calculated according to the Cox-de Boor recursion formula:(3)Ni,pξ=ξ-ξiξi+p-ξiNi,p-1ξ+ξi+p+1-ξξi+p+1-ξi+1Ni+1,p-1ξ.Apth-order NURBS basis function of R
i
,
p
ξ is defined by giving each N
i
,
p
ξ a proper weight, and a pth-order NURBS is constructed by n
pth-order NURBS basis functions and their corresponding control points A
i:(4)CNξ=∑i=1nRi,pξAi,Ri,pξ=Ni,pξωiWξ=Ni,pξωi∑i^=1nNi^,pξωi^,where ω
i
i
=
1,2
,
…
,
n is the weight for N
i
,
p
ξ. NURBS basis function degenerates to B-spline basis function when all the weights take the same value.The NURBS solid is a tensor product ofN
i
,
p
ξ, M
j
,
q
η, and P
k
,
r
ζ:(5)Vξ,η,ζ=∑i=1n∑j=1m∑k=1lRi,j,kp,q,rξ,η,ζAi,j,k,Ri,j,kp,q,rξ,η,ζ=Ni,pξMj,qηPk,rζωi,j,k∑i^=1n∑j^=1m∑k^=1lNi^,pξMj^,qηPk^,rζωi^,j^,k^,where n, m, and l refer to the number of B-spline basis functions of N
i
,
p
ξ, M
j
,
q
η, and P
k
,
r
ζ in the directions of u, v, and w. A
i
,
j
,
k is the control point of R
i
,
j
,
k
p
,
q
,
r and ω
i
,
j
,
k is the weight.


### 2.2. Numerical Formulation for NURBS Solid Model

The numerical integral of a given functionB (i.e., stiffness matrix K) on a domain Ω consists of n
e NURBS solids that can be achieved via the philosophy of isoparametric transformation. According to the coordinate systems where the NURBS solids locate in the process of transformation, three spaces, that is, the parent space Ω
~
e
l
e, the parametric space Ω
^
e
l
e, and the physical space Ω
e
l
e, are defined by the local, knot, and Cartesian coordinates, respectively. Hence,(6)Ω~ele=ξ~i,ξ~i+1⊗η~j,η~j+1⊗ζ~k,ζ~k+1=-1,1⊗-1,1⊗-1,1,(7)Ω^ele=ξi,ξi+1⊗ηj,ηj+1⊗ζk,ζk+1,ξi≠ξi+1,ηj≠ηj+1,ζk≠ζk+1,(8)Ωele=xi,xi+1⊗yj,yj+1⊗zk,zk+1,xi≠xi+1,yj≠yj+1,zk≠zk+1.As shown in Figure1, we define the mappings of parametric space → parent space and physical space → parametric space as mapping 1 and mapping 2, respectively. For a NURBS solid, if the coordinates of analyzed point in the parent, parametric, and physical spaces are ξ
~
,
η
~
,
ζ
~, ξ
,
η
,
ζ, and x
,
y
,
z, respectively, then the Jacobian matrices for mapping 1 and mapping 2 are(9)
J
1
=
∂
ξ
,
η
,
ζ
∂
ξ
~
,
η
~
,
ζ
~
=
ξ
i
+
1
-
ξ
i
2
0
0
0
η
j
+
1
-
η
j
2
0
0
0
ζ
k
+
1
-
ζ
k
2
,
J
2
=
∂
x
,
y
,
z
∂
ξ
,
η
,
ζ
=
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ξ
x
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ξ
y
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ξ
z
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
η
x
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
η
y
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
η
z
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ζ
x
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ζ
y
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ζ
z
i
p
t
,where n
p
t is the total number of control points for a single NURBS solid and x
i
p
t
,
y
i
p
t
,
z
i
p
t is the Cartesian coordinates of control point i
p
t.Figure 1
Schematic and mappings ofΩ
~
e
l
e, Ω
^
e
l
e, and Ω
e
l
e in two-dimensional space.The integral of functionB on Ω can be derived as(10)∫ΩBx,y,zdΩ=∑ele=1ne∫ΩeleBx,y,zdΩele=∑ele=1ne∫Ω~eleBξ~,η~,ζ~JdΩ~elewith determinant J
=
J
2
·
J
1.


### 2.3. Comparison with Traditional FEM

Although the numerical formulations of IGA and FE models are both constructed on the philosophy of isoparametric transformation, the basic ideas of these two are intrinsically different. The FE model, whose basis functions are chosen as Lagrange functions, is only an approximation to the original geometric model. The IGA model, on the other hand, with its element bases directly adopting the NURBS basis functions that exactly describe the geometric model, has no loss of information between the analyzed model and the original one, providing an extremely apparent advantage for curved or complicated shapes, as shown in Figure2.Figure 2
Comparison on geometric modeling of a circle by FEM and IGA.
(a)
FEM modeling
(b)
IGA modelingConsidering a cubic as shown in Table1 and discretizing the cubic to n
·
n
·
n
pth-order elements, respectively, via FEM and IGA methods, we can calculate the total number of nodes in FE model (NFEM) and the total number of control points in IGA model (NIGA) as functions of n and p. If linear basis functions are used, then NFEM = NIGA. Once the order of basis functions is elevated, as seen in Figure 3, NIGA is apparently smaller than NFEM for the same value of n, indicating that the scale of the equilibrium equation of IGA model is significantly reduced. It is concluded, as for Figure 4, that NIGA is not sensitive to the order changes of basis functions. Even though quartic NURBS basis functions are used, the scale of the equilibrium equation of IGA model is smaller than that of quadratic FE model. It should also be mentioned that the conduction of mesh refinement in IGA model is more convenient than that in FE model, since order elevation and knot insertion can be directly carried out on its linear control model.Table 1
Comparison between NFEM and NIGA.
Analyzing model
Order of basis function
NFEM
NIGA
p
=
1
n
+
1
3
n
+
1
3
p
=
2
n
+
1
2
4
n
+
1
n
+
2
3
p
=
3
n
+
1
2
7
n
+
1
n
+
3
3
⋮
⋮
⋮
p
n
+
1
2
3
p
-
2
n
+
1
n
+
p
3Figure 3
Evolution curves of NFEM and NIGA with respect ton.
(a)
p
=
2
(b)
p
=
3Figure 4
Evolution curves of NIGA with different basis orders.


## 2.1.B-Spline Basis Function and NURBS Basis Function

The basic analyzed object of IGA [24] in 1D, 2D, and 3D space is the NURBS (Nonuniform Rational B-Spline), NURBS surface, and NURBS solid, respectively. A pth-order B-spline C
B
ξ can be constructed by n
B-spline basis functions N
i
,
p
ξ defined on the knot sequence Θ
=
ξ
1
,
ξ
2
,
…
,
ξ
n
+
p
+
1 and their corresponding control points A
i:(1)CBξ=∑i=1nNi,pξAi.The coordinates of knots should be monotonically nondecreasing (ξ
1
≤
⋯
≤
ξ
i
≤
ξ
i
+
1
≤
⋯
≤
ξ
n
+
p
+
1), with two nonoverlapped neighboring knots constructing an element interval.Forp
=
0, N
i
,
0
ξ is defined as(2)Ni,0ξ=1ifξi≤ξ≤ξi+10otherwise.Forp
≥
1, N
i
,
p
ξ can be calculated according to the Cox-de Boor recursion formula:(3)Ni,pξ=ξ-ξiξi+p-ξiNi,p-1ξ+ξi+p+1-ξξi+p+1-ξi+1Ni+1,p-1ξ.Apth-order NURBS basis function of R
i
,
p
ξ is defined by giving each N
i
,
p
ξ a proper weight, and a pth-order NURBS is constructed by n
pth-order NURBS basis functions and their corresponding control points A
i:(4)CNξ=∑i=1nRi,pξAi,Ri,pξ=Ni,pξωiWξ=Ni,pξωi∑i^=1nNi^,pξωi^,where ω
i
i
=
1,2
,
…
,
n is the weight for N
i
,
p
ξ. NURBS basis function degenerates to B-spline basis function when all the weights take the same value.The NURBS solid is a tensor product ofN
i
,
p
ξ, M
j
,
q
η, and P
k
,
r
ζ:(5)Vξ,η,ζ=∑i=1n∑j=1m∑k=1lRi,j,kp,q,rξ,η,ζAi,j,k,Ri,j,kp,q,rξ,η,ζ=Ni,pξMj,qηPk,rζωi,j,k∑i^=1n∑j^=1m∑k^=1lNi^,pξMj^,qηPk^,rζωi^,j^,k^,where n, m, and l refer to the number of B-spline basis functions of N
i
,
p
ξ, M
j
,
q
η, and P
k
,
r
ζ in the directions of u, v, and w. A
i
,
j
,
k is the control point of R
i
,
j
,
k
p
,
q
,
r and ω
i
,
j
,
k is the weight.


## 2.2. Numerical Formulation for NURBS Solid Model

The numerical integral of a given functionB (i.e., stiffness matrix K) on a domain Ω consists of n
e NURBS solids that can be achieved via the philosophy of isoparametric transformation. According to the coordinate systems where the NURBS solids locate in the process of transformation, three spaces, that is, the parent space Ω
~
e
l
e, the parametric space Ω
^
e
l
e, and the physical space Ω
e
l
e, are defined by the local, knot, and Cartesian coordinates, respectively. Hence,(6)Ω~ele=ξ~i,ξ~i+1⊗η~j,η~j+1⊗ζ~k,ζ~k+1=-1,1⊗-1,1⊗-1,1,(7)Ω^ele=ξi,ξi+1⊗ηj,ηj+1⊗ζk,ζk+1,ξi≠ξi+1,ηj≠ηj+1,ζk≠ζk+1,(8)Ωele=xi,xi+1⊗yj,yj+1⊗zk,zk+1,xi≠xi+1,yj≠yj+1,zk≠zk+1.As shown in Figure1, we define the mappings of parametric space → parent space and physical space → parametric space as mapping 1 and mapping 2, respectively. For a NURBS solid, if the coordinates of analyzed point in the parent, parametric, and physical spaces are ξ
~
,
η
~
,
ζ
~, ξ
,
η
,
ζ, and x
,
y
,
z, respectively, then the Jacobian matrices for mapping 1 and mapping 2 are(9)
J
1
=
∂
ξ
,
η
,
ζ
∂
ξ
~
,
η
~
,
ζ
~
=
ξ
i
+
1
-
ξ
i
2
0
0
0
η
j
+
1
-
η
j
2
0
0
0
ζ
k
+
1
-
ζ
k
2
,
J
2
=
∂
x
,
y
,
z
∂
ξ
,
η
,
ζ
=
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ξ
x
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ξ
y
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ξ
z
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
η
x
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
η
y
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
η
z
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ζ
x
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ζ
y
i
p
t
∑
i
p
t
=
1
n
p
t
∂
R
i
p
t
ξ
,
η
,
ζ
∂
ζ
z
i
p
t
,where n
p
t is the total number of control points for a single NURBS solid and x
i
p
t
,
y
i
p
t
,
z
i
p
t is the Cartesian coordinates of control point i
p
t.Figure 1
Schematic and mappings ofΩ
~
e
l
e, Ω
^
e
l
e, and Ω
e
l
e in two-dimensional space.The integral of functionB on Ω can be derived as(10)∫ΩBx,y,zdΩ=∑ele=1ne∫ΩeleBx,y,zdΩele=∑ele=1ne∫Ω~eleBξ~,η~,ζ~JdΩ~elewith determinant J
=
J
2
·
J
1.


## 2.3. Comparison with Traditional FEM

Although the numerical formulations of IGA and FE models are both constructed on the philosophy of isoparametric transformation, the basic ideas of these two are intrinsically different. The FE model, whose basis functions are chosen as Lagrange functions, is only an approximation to the original geometric model. The IGA model, on the other hand, with its element bases directly adopting the NURBS basis functions that exactly describe the geometric model, has no loss of information between the analyzed model and the original one, providing an extremely apparent advantage for curved or complicated shapes, as shown in Figure2.Figure 2
Comparison on geometric modeling of a circle by FEM and IGA.
(a)
FEM modeling
(b)
IGA modelingConsidering a cubic as shown in Table1 and discretizing the cubic to n
·
n
·
n
pth-order elements, respectively, via FEM and IGA methods, we can calculate the total number of nodes in FE model (NFEM) and the total number of control points in IGA model (NIGA) as functions of n and p. If linear basis functions are used, then NFEM = NIGA. Once the order of basis functions is elevated, as seen in Figure 3, NIGA is apparently smaller than NFEM for the same value of n, indicating that the scale of the equilibrium equation of IGA model is significantly reduced. It is concluded, as for Figure 4, that NIGA is not sensitive to the order changes of basis functions. Even though quartic NURBS basis functions are used, the scale of the equilibrium equation of IGA model is smaller than that of quadratic FE model. It should also be mentioned that the conduction of mesh refinement in IGA model is more convenient than that in FE model, since order elevation and knot insertion can be directly carried out on its linear control model.Table 1
Comparison between NFEM and NIGA.
Analyzing model
Order of basis function
NFEM
NIGA
p
=
1
n
+
1
3
n
+
1
3
p
=
2
n
+
1
2
4
n
+
1
n
+
2
3
p
=
3
n
+
1
2
7
n
+
1
n
+
3
3
⋮
⋮
⋮
p
n
+
1
2
3
p
-
2
n
+
1
n
+
p
3Figure 3
Evolution curves of NFEM and NIGA with respect ton.
(a)
p
=
2
(b)
p
=
3Figure 4
Evolution curves of NIGA with different basis orders.


## 3. Constitutive Model of Concrete

### 3.1. Basic Functions

The basic functions for elastoplastic damage model are(11)ε=εe+εp,σ=1-dσ-,σ-=Deεe,ε˙p=λ˙mσ-,κ,κ˙=λ˙hσ-,κ,where ε, ε
e, and ε
p are the total, the elastic, and the plastic strains, respectively; σ and σ
- are the nominal and the effective stresses, respectively; d is the damage; D
e is the initial elastic stiffness matrix; m is the flow vector that can be derived from the plastic potential Φ; h is the plastic modulus; κ is the set of internal variables; λ
˙ is the plastic multiplier that satisfies the Kuhn-Tucker condition(12)λ˙≥0,Fσ-,κ≤0withλ˙F=0,where F
σ
-
,
κ is the yield function.Two internal variables ofκ
t and κ
c are introduced to record the hardening/softening behaviors of concrete under tension and compression, respectively; that is, κ
=
κ
t
,
κ
c
=
ε
~
t
p
,
ε
~
c
p, where ε
~
t
p and ε
~
c
p are the tensile and compressive equivalent plastic strains, respectively.


### 3.2. Yield Function

The yield surface of concrete is described by the modified Barcelona function(13)
F
σ
-
=
1
1
-
α
F
α
F
I
-
1
+
3
J
-
2
+
β
F
κ
σ
-
^
m
a
x
-
γ
F
-
σ
-
^
m
a
x
-
c
-
c
κ
cwith(14)αF=θ-12θ-1,0≤θ≤5,βFκ=c-cκcc-tκt1-αF-1+αF,γF=31-Kc2Kc-1,where I
-
1 and J
-
2 are the first stress invariant and the second deviatoric stress invariant computed by the effective stress tensor σ
-, respectively; σ
-
^
m
a
x is algebraically the maximum principal effective stress with x
=
x
+
x
/
2, the Macaulay bracket function; θ
=
f
b
0
/
f
c
0 with f
b
0 and f
c
0 being the uniaxial and biaxial yield strengths of concrete in compression, respectively; c
-
t
κ
t and c
-
c
κ
c are the effective tensile and compressive cohesion stresses of concrete, respectively; γ
F is introduced only when concrete is subjected to triaxial compression with 0.5
<
K
c
<
1.0.


### 3.3. Plastic Potential

A hyperbolic form of the Drucker-Prager function is used as the plastic potential of concrete:(15)Φ=ρft0tan⁡φ2+3J-2+tan⁡φ3I-1,where φ is the dilation angle; f
t
0 is the uniaxial tensile yield stress; ρ is the eccentricity parameter.


### 3.4. Plastic Modulus

The plastic modulus is written as(16)hσ-^,κ=rσ-^0000-1-rσ-^·∂Φ∂σ-^with r
σ
-
^ being the weight function which can be computed as(17)rσ-^=0ifσ-^=0∑i=13σ-^i∑i=13σ-^iotherwise.


### 3.5. Damage Law

The total damaged can be seen as the combination of tensile and compressive damage, which has been discussed in detail by Lubliner et al. [25] and Lee and Fenves [13] and then furtherly developed in ABAQUS [26]:(18)d=1-1-stdcκc1-scdtκtwith(19)st=1-wtrσ-^,0≤wt≤1,sc=1-wc1-rσ-^,0≤wc≤1,where d
t
=
d
t
κ
t and d
c
=
d
c
κ
c are the tensile and compressive damage variables, respectively; w
t and w
c are the stiffness recovery factors reflecting the unilateral effects of concrete in tension and compression, respectively, which are usually set as w
t
=
0 and w
c
=
1.Based on the experimental data, Yu and Wu [27] assumed that the plastic strain of concrete under uniaxial loading satisfies(20)εℵp=αεℵ-ε0ℵβ,where ℵ
∈
t
,
c is the state variable with ℵ
=
t and ℵ
=
c referring to uniaxial tension and compression, respectively; α and β are calibrated parameters and are commonly set as α
=
1.79 and β
=
1.18; ε
ℵ is the total strain and ε
0
ℵ is the strain corresponding to the elasticity limit.For uniaxial loading,κ
ℵ
=
ε
ℵ
p, the effective stress of concrete can be calculated as(21)σ-ℵκℵ=Eℵ0εℵκℵ=0Eℵ0ε0ℵ+κℵα1/β-κℵκℵ>0.Yu and Wu [27] also recommended a piecewise function to describe the damage evolution:(22)dℵκℵ=d0ℵκℵ=0d0ℵ+A1ℵκℵαB1ℵ/βκℵ≤κfℵ1-1-dfℵΠℵB2ℵA2ℵΠℵ-1B3ℵ+Πℵκℵ>κfℵwithΠℵ=ε0ℵ+κℵ/α1/βεfℵ,where d
0
ℵ
=
0 is the initial damage; d
f
ℵ and κ
f
ℵ are the damage and internal variable of concrete at the peak strain ε
f
ℵ; A
2
ℵ and B
3
ℵ are parameters controlling the softening curve of concrete; A
1
ℵ, B
1
ℵ, and B
2
ℵ can be computed as below according to the continuity requirements:(23)A1ℵ=1-d0ℵ-λℵεfℵ-ε0ℵ-B1ℵ,B1ℵ=λℵσ-ℵ′εfℵεfℵ-ε0ℵσℵεfℵ1-d0ℵ-λℵ,B2ℵ=σℵεfℵ-εfℵσ-ℵ′εfℵσℵεfℵ,λℵ=σℵεfℵσ-ℵεfℵ=1-dfℵ,where λ
ℵ is the ratio of nominal stress to effective stress at ε
f
ℵ; σ
-
ℵ
′
ε
f
ℵ is the derivative of effective stress with respect to strain at ε
f
ℵ. The damage law is selected not only because of its applicability in recording both the tensile and the compressive degeneration mechanisms of concrete, but also because of its feasibility in adjusting the softening behaviors of concrete. The influences of A
2
ℵ and B
3
ℵ on the evolution discipline of d
ℵ are shown in Figure 5.Figure 5
Influences ofA
2
ℵ and B
3
ℵ on the evolution discipline of d
ℵ.
(a)
Influence ofA
2
t
(b)
Influence ofA
2
c
(c)
Influence ofB
3
t
(d)
Influence ofB
3
c


## 3.1. Basic Functions

The basic functions for elastoplastic damage model are(11)ε=εe+εp,σ=1-dσ-,σ-=Deεe,ε˙p=λ˙mσ-,κ,κ˙=λ˙hσ-,κ,where ε, ε
e, and ε
p are the total, the elastic, and the plastic strains, respectively; σ and σ
- are the nominal and the effective stresses, respectively; d is the damage; D
e is the initial elastic stiffness matrix; m is the flow vector that can be derived from the plastic potential Φ; h is the plastic modulus; κ is the set of internal variables; λ
˙ is the plastic multiplier that satisfies the Kuhn-Tucker condition(12)λ˙≥0,Fσ-,κ≤0withλ˙F=0,where F
σ
-
,
κ is the yield function.Two internal variables ofκ
t and κ
c are introduced to record the hardening/softening behaviors of concrete under tension and compression, respectively; that is, κ
=
κ
t
,
κ
c
=
ε
~
t
p
,
ε
~
c
p, where ε
~
t
p and ε
~
c
p are the tensile and compressive equivalent plastic strains, respectively.


## 3.2. Yield Function

The yield surface of concrete is described by the modified Barcelona function(13)
F
σ
-
=
1
1
-
α
F
α
F
I
-
1
+
3
J
-
2
+
β
F
κ
σ
-
^
m
a
x
-
γ
F
-
σ
-
^
m
a
x
-
c
-
c
κ
cwith(14)αF=θ-12θ-1,0≤θ≤5,βFκ=c-cκcc-tκt1-αF-1+αF,γF=31-Kc2Kc-1,where I
-
1 and J
-
2 are the first stress invariant and the second deviatoric stress invariant computed by the effective stress tensor σ
-, respectively; σ
-
^
m
a
x is algebraically the maximum principal effective stress with x
=
x
+
x
/
2, the Macaulay bracket function; θ
=
f
b
0
/
f
c
0 with f
b
0 and f
c
0 being the uniaxial and biaxial yield strengths of concrete in compression, respectively; c
-
t
κ
t and c
-
c
κ
c are the effective tensile and compressive cohesion stresses of concrete, respectively; γ
F is introduced only when concrete is subjected to triaxial compression with 0.5
<
K
c
<
1.0.


## 3.3. Plastic Potential

A hyperbolic form of the Drucker-Prager function is used as the plastic potential of concrete:(15)Φ=ρft0tan⁡φ2+3J-2+tan⁡φ3I-1,where φ is the dilation angle; f
t
0 is the uniaxial tensile yield stress; ρ is the eccentricity parameter.


## 3.4. Plastic Modulus

The plastic modulus is written as(16)hσ-^,κ=rσ-^0000-1-rσ-^·∂Φ∂σ-^with r
σ
-
^ being the weight function which can be computed as(17)rσ-^=0ifσ-^=0∑i=13σ-^i∑i=13σ-^iotherwise.


## 3.5. Damage Law

The total damaged can be seen as the combination of tensile and compressive damage, which has been discussed in detail by Lubliner et al. [25] and Lee and Fenves [13] and then furtherly developed in ABAQUS [26]:(18)d=1-1-stdcκc1-scdtκtwith(19)st=1-wtrσ-^,0≤wt≤1,sc=1-wc1-rσ-^,0≤wc≤1,where d
t
=
d
t
κ
t and d
c
=
d
c
κ
c are the tensile and compressive damage variables, respectively; w
t and w
c are the stiffness recovery factors reflecting the unilateral effects of concrete in tension and compression, respectively, which are usually set as w
t
=
0 and w
c
=
1.Based on the experimental data, Yu and Wu [27] assumed that the plastic strain of concrete under uniaxial loading satisfies(20)εℵp=αεℵ-ε0ℵβ,where ℵ
∈
t
,
c is the state variable with ℵ
=
t and ℵ
=
c referring to uniaxial tension and compression, respectively; α and β are calibrated parameters and are commonly set as α
=
1.79 and β
=
1.18; ε
ℵ is the total strain and ε
0
ℵ is the strain corresponding to the elasticity limit.For uniaxial loading,κ
ℵ
=
ε
ℵ
p, the effective stress of concrete can be calculated as(21)σ-ℵκℵ=Eℵ0εℵκℵ=0Eℵ0ε0ℵ+κℵα1/β-κℵκℵ>0.Yu and Wu [27] also recommended a piecewise function to describe the damage evolution:(22)dℵκℵ=d0ℵκℵ=0d0ℵ+A1ℵκℵαB1ℵ/βκℵ≤κfℵ1-1-dfℵΠℵB2ℵA2ℵΠℵ-1B3ℵ+Πℵκℵ>κfℵwithΠℵ=ε0ℵ+κℵ/α1/βεfℵ,where d
0
ℵ
=
0 is the initial damage; d
f
ℵ and κ
f
ℵ are the damage and internal variable of concrete at the peak strain ε
f
ℵ; A
2
ℵ and B
3
ℵ are parameters controlling the softening curve of concrete; A
1
ℵ, B
1
ℵ, and B
2
ℵ can be computed as below according to the continuity requirements:(23)A1ℵ=1-d0ℵ-λℵεfℵ-ε0ℵ-B1ℵ,B1ℵ=λℵσ-ℵ′εfℵεfℵ-ε0ℵσℵεfℵ1-d0ℵ-λℵ,B2ℵ=σℵεfℵ-εfℵσ-ℵ′εfℵσℵεfℵ,λℵ=σℵεfℵσ-ℵεfℵ=1-dfℵ,where λ
ℵ is the ratio of nominal stress to effective stress at ε
f
ℵ; σ
-
ℵ
′
ε
f
ℵ is the derivative of effective stress with respect to strain at ε
f
ℵ. The damage law is selected not only because of its applicability in recording both the tensile and the compressive degeneration mechanisms of concrete, but also because of its feasibility in adjusting the softening behaviors of concrete. The influences of A
2
ℵ and B
3
ℵ on the evolution discipline of d
ℵ are shown in Figure 5.Figure 5
Influences ofA
2
ℵ and B
3
ℵ on the evolution discipline of d
ℵ.
(a)
Influence ofA
2
t
(b)
Influence ofA
2
c
(c)
Influence ofB
3
t
(d)
Influence ofB
3
c


## 4. Solution of Local Problem

The solution of local problem, also known as stress updating, can be decomposed into three steps according to the concept of operator split [15, 16, 28], that is, the elastic predictor, the plastic corrector, and the damage corrector. The elastic predictor and the plastic corrector are executed to obtain the effective stresses, the internal variables, and the plastic multiplier that satisfy the updated yield surface, and the damage corrector is used to compute damage and nominal stresses.


### 4.1. Elastic Predictor and Plastic Corrector

Assuming that the state with variables{
ε
n
,
ε
p
n
,
κ
n
} is already known at time t
n, then the effective stress σ
-
n and damage d
n can be computed as(24)σ-n=De·εen=Deεn-εpnwithdn=dκn.For time incrementΔ
t, we give t
n
+
1
=
t
n
+
Δ
t. Let the strain increment of concrete within Δ
t be Δ
ε; the total strain and effective stress of concrete at time t
n
+
1 then satisfy(25)εn+1=εn+Δεσ-n+1=σ-trialn+1-λn+1Demσ-n+1,κn+1withσ-trialn+1=Deεn+1-εpn,where σ
-
t
r
i
a
l
n
+
1 is the trial stress at time t
n
+
1. If F
(
σ
-
t
r
i
a
l
n
+
1
,
κ
n
)
≤
0, then concrete still remains in elasticity and the trial values are the final values of t
n
+
1; that is, σ
-
n
+
1
=
σ
-
t
r
i
a
l
n
+
1, κ
n
+
1
=
κ
n, and λ
n
+
1
=
0; if F
(
σ
-
t
r
i
a
l
n
+
1
,
κ
n
)
>
0, then plasticity develops and all trial values should be recalculated. Let all the unknowns to be solved at t
n
+
1 be x
n
+
1
=
σ
-
n
+
1
κ
n
+
1
λ
n
+
1
T; we then have(26)fxn+1=σ-n+1+λn+1Demσ-n+1,κn+1κn+1-λn+1hσ-n+1,κn+1Fσ-n+1,κn+1=Inσ-0nσ-,nκ0nσ-,10nκ,nσ-Inκ0nκ,101,nσ-01,nκ0σ-nκnλn+De0nκ,nσ-01,nσ-Δε,where n
σ
-
(
=
6
) and n
κ
(
=
2
) are the dimensions of the stress tensor and the number of the internal variables, respectively; I
∗ is the identity matrix of order ∗. A backward Euler scheme is used in the iterative solution of (26). Given a specific iteration step k within Δ
t, that is, x
k
n
+
1
=
σ
-
k
n
+
1
κ
k
n
+
1
λ
k
n
+
1
T, for the next iteration step k
+
1, x
k
+
1
n
+
1 can be derived as(27)xk+1n+1=xkn+1+Δxkn+1withΔxkn+1=-JLocalkn+1-1Exkn+1,where E
x
k
n
+
1 and J
L
o
c
a
l
k
n
+
1 are the residual and Jacobian matrices of (26) with the following expressions:(28)
E
x
k
n
+
1
=
E
σ
-
k
n
+
1
E
κ
k
n
+
1
E
λ
k
n
+
1
=
σ
-
k
n
+
1
+
λ
k
n
+
1
D
e
m
σ
-
k
n
+
1
,
κ
k
n
+
1
-
σ
-
t
r
i
a
l
n
+
1
κ
k
n
+
1
-
λ
k
n
+
1
h
σ
-
k
n
+
1
,
κ
k
n
+
1
-
κ
n
F
σ
-
k
n
+
1
,
κ
k
n
+
1
,
J
k
n
+
1
=
d
f
d
σ
-
k
n
+
1
,
κ
k
n
+
1
,
λ
k
n
+
1
=
I
n
σ
-
+
λ
D
e
∂
m
∂
σ
-
λ
D
e
∂
m
∂
κ
D
e
m
-
λ
∂
h
∂
σ
-
I
n
κ
-
λ
∂
h
∂
κ
-
h
∂
F
∂
σ
-
T
∂
F
∂
κ
T
0
k
n
+
1
,where T refers to the transposition.For the first iteration step, the initial solution of all the unknowns of (26) is usually chosen from the elastic trial state; that is, x
0
n
+
1
=
σ
-
t
r
i
a
l
n
+
1
κ
n
0
T. Convergence is satisfied when the Euclidean norm of the residuals is less than given tolerance; that is, E
x
k
n
+
1
⩽
T
O
L
L
o
c
a
l.


### 4.2. Substepping Strategy

The convergence of (26) may not be easily made in some cases, especially when the strain increment Δ
ε within Δ
t is large. Therefore, an adaptive substepping strategy [17, 29] is employed if the convergence cannot be made for given tolerance or a maximum number of iterations in the solution. The total strain increment Δ
ε from t
n to t
n
+
1 can then be divided into m parts according to (not necessary equal)(29)Δε=Δε1+Δε2+⋯+Δεm=∑k=1mΔεk=∑k=1mαkΔε,where α
k
=
Δ
ε
k
/
Δ
ε.


### 4.3. Damage Corrector

Once the solutions of elastic predictor and plastic corrector have been finished, the nominal stress can be updated as(30)σn+1=1-dn+1σ-n+1,dn+1=dκn+1.


### 4.4. Consistent Tangent Stiffness

The consistent tangent stiffness at timet
n
+
1 can be calculated by differentiating both sides of (11) with respect to the strain increment:(31)dσn+1dΔε=-σ-n+1ddn+1dκn+1dκn+1dΔε+1-dn+1dσ-n+1dΔε.If substepping strategy is used, the computations ofd
κ
n
+
1
/
d
Δ
ε and d
σ
-
n
+
1
/
d
Δ
ε should be derived step by step recursively. For substep k
(
k
≠
1
), (26) can be written as(32)fσ-n+k/m,κn+k/m,λn+k/m=Inσ-0nσ-,nκ0nσ-,10nκ,nσ-Inκ0nκ,101,nσ-01,nκ0σ-n+k-1/mκn+k-1/mλn+k-1/m+De0nκ,nσ-01,nσ-Δεk.Let(33)dfdσ-n+k/m,κn+k/m,λn+k/m=JLocaln+k/mdΔεkdΔε=αkInσ-;and then differentiating both sides of (32) with respect to Δ
ε we obtain (34)
d
σ
-
n
+
k
/
m
,
κ
n
+
k
/
m
,
λ
n
+
k
/
m
d
Δ
ε
=
A
n
+
k
/
m
d
σ
-
n
+
k
-
1
/
m
,
κ
n
+
k
-
1
/
m
,
λ
n
+
k
-
1
/
m
d
Δ
ε
+
α
k
P
D
e
,where A
n
+
(
k
/
m
)
=
J
L
o
c
a
l
n
+
(
k
/
m
)
-
1
I
n
σ
-
+
n
κ
0
n
σ
-
+
n
κ
,
1
0
1
,
n
σ
-
+
n
κ
0 and P
T
=
(
I
n
σ
-
,
0
n
σ
-
,
n
κ
+
1
).For substepk
=
1,(35)dσ-n+1/m,κn+1/m,λn+1/mdΔε=α1APDen+1/m.Combining (34) and (35), we get the following equation according to the recursive relation: (36)
d
σ
-
n
+
1
,
κ
n
+
1
,
λ
n
+
1
d
Δ
ε
=
A
n
+
1
A
n
+
m
-
1
/
m
⋯
α
2
P
D
e
+
α
1
A
P
D
e
n
+
1
/
m
⋯
+
α
m
-
1
P
D
e
+
α
m
P
D
e
=
Ξ
σ
-
Ξ
κ
•
,where uninterested submatrix in (36) is represented by bold dots •. Obviously, d
σ
-
n
+
1
/
d
Δ
ε
=
Ξ
σ
- and d
κ
n
+
1
/
d
Δ
ε
=
Ξ
κ.


## 4.1. Elastic Predictor and Plastic Corrector

Assuming that the state with variables{
ε
n
,
ε
p
n
,
κ
n
} is already known at time t
n, then the effective stress σ
-
n and damage d
n can be computed as(24)σ-n=De·εen=Deεn-εpnwithdn=dκn.For time incrementΔ
t, we give t
n
+
1
=
t
n
+
Δ
t. Let the strain increment of concrete within Δ
t be Δ
ε; the total strain and effective stress of concrete at time t
n
+
1 then satisfy(25)εn+1=εn+Δεσ-n+1=σ-trialn+1-λn+1Demσ-n+1,κn+1withσ-trialn+1=Deεn+1-εpn,where σ
-
t
r
i
a
l
n
+
1 is the trial stress at time t
n
+
1. If F
(
σ
-
t
r
i
a
l
n
+
1
,
κ
n
)
≤
0, then concrete still remains in elasticity and the trial values are the final values of t
n
+
1; that is, σ
-
n
+
1
=
σ
-
t
r
i
a
l
n
+
1, κ
n
+
1
=
κ
n, and λ
n
+
1
=
0; if F
(
σ
-
t
r
i
a
l
n
+
1
,
κ
n
)
>
0, then plasticity develops and all trial values should be recalculated. Let all the unknowns to be solved at t
n
+
1 be x
n
+
1
=
σ
-
n
+
1
κ
n
+
1
λ
n
+
1
T; we then have(26)fxn+1=σ-n+1+λn+1Demσ-n+1,κn+1κn+1-λn+1hσ-n+1,κn+1Fσ-n+1,κn+1=Inσ-0nσ-,nκ0nσ-,10nκ,nσ-Inκ0nκ,101,nσ-01,nκ0σ-nκnλn+De0nκ,nσ-01,nσ-Δε,where n
σ
-
(
=
6
) and n
κ
(
=
2
) are the dimensions of the stress tensor and the number of the internal variables, respectively; I
∗ is the identity matrix of order ∗. A backward Euler scheme is used in the iterative solution of (26). Given a specific iteration step k within Δ
t, that is, x
k
n
+
1
=
σ
-
k
n
+
1
κ
k
n
+
1
λ
k
n
+
1
T, for the next iteration step k
+
1, x
k
+
1
n
+
1 can be derived as(27)xk+1n+1=xkn+1+Δxkn+1withΔxkn+1=-JLocalkn+1-1Exkn+1,where E
x
k
n
+
1 and J
L
o
c
a
l
k
n
+
1 are the residual and Jacobian matrices of (26) with the following expressions:(28)
E
x
k
n
+
1
=
E
σ
-
k
n
+
1
E
κ
k
n
+
1
E
λ
k
n
+
1
=
σ
-
k
n
+
1
+
λ
k
n
+
1
D
e
m
σ
-
k
n
+
1
,
κ
k
n
+
1
-
σ
-
t
r
i
a
l
n
+
1
κ
k
n
+
1
-
λ
k
n
+
1
h
σ
-
k
n
+
1
,
κ
k
n
+
1
-
κ
n
F
σ
-
k
n
+
1
,
κ
k
n
+
1
,
J
k
n
+
1
=
d
f
d
σ
-
k
n
+
1
,
κ
k
n
+
1
,
λ
k
n
+
1
=
I
n
σ
-
+
λ
D
e
∂
m
∂
σ
-
λ
D
e
∂
m
∂
κ
D
e
m
-
λ
∂
h
∂
σ
-
I
n
κ
-
λ
∂
h
∂
κ
-
h
∂
F
∂
σ
-
T
∂
F
∂
κ
T
0
k
n
+
1
,where T refers to the transposition.For the first iteration step, the initial solution of all the unknowns of (26) is usually chosen from the elastic trial state; that is, x
0
n
+
1
=
σ
-
t
r
i
a
l
n
+
1
κ
n
0
T. Convergence is satisfied when the Euclidean norm of the residuals is less than given tolerance; that is, E
x
k
n
+
1
⩽
T
O
L
L
o
c
a
l.


## 4.2. Substepping Strategy

The convergence of (26) may not be easily made in some cases, especially when the strain increment Δ
ε within Δ
t is large. Therefore, an adaptive substepping strategy [17, 29] is employed if the convergence cannot be made for given tolerance or a maximum number of iterations in the solution. The total strain increment Δ
ε from t
n to t
n
+
1 can then be divided into m parts according to (not necessary equal)(29)Δε=Δε1+Δε2+⋯+Δεm=∑k=1mΔεk=∑k=1mαkΔε,where α
k
=
Δ
ε
k
/
Δ
ε.


## 4.3. Damage Corrector

Once the solutions of elastic predictor and plastic corrector have been finished, the nominal stress can be updated as(30)σn+1=1-dn+1σ-n+1,dn+1=dκn+1.


## 4.4. Consistent Tangent Stiffness

The consistent tangent stiffness at timet
n
+
1 can be calculated by differentiating both sides of (11) with respect to the strain increment:(31)dσn+1dΔε=-σ-n+1ddn+1dκn+1dκn+1dΔε+1-dn+1dσ-n+1dΔε.If substepping strategy is used, the computations ofd
κ
n
+
1
/
d
Δ
ε and d
σ
-
n
+
1
/
d
Δ
ε should be derived step by step recursively. For substep k
(
k
≠
1
), (26) can be written as(32)fσ-n+k/m,κn+k/m,λn+k/m=Inσ-0nσ-,nκ0nσ-,10nκ,nσ-Inκ0nκ,101,nσ-01,nκ0σ-n+k-1/mκn+k-1/mλn+k-1/m+De0nκ,nσ-01,nσ-Δεk.Let(33)dfdσ-n+k/m,κn+k/m,λn+k/m=JLocaln+k/mdΔεkdΔε=αkInσ-;and then differentiating both sides of (32) with respect to Δ
ε we obtain (34)
d
σ
-
n
+
k
/
m
,
κ
n
+
k
/
m
,
λ
n
+
k
/
m
d
Δ
ε
=
A
n
+
k
/
m
d
σ
-
n
+
k
-
1
/
m
,
κ
n
+
k
-
1
/
m
,
λ
n
+
k
-
1
/
m
d
Δ
ε
+
α
k
P
D
e
,where A
n
+
(
k
/
m
)
=
J
L
o
c
a
l
n
+
(
k
/
m
)
-
1
I
n
σ
-
+
n
κ
0
n
σ
-
+
n
κ
,
1
0
1
,
n
σ
-
+
n
κ
0 and P
T
=
(
I
n
σ
-
,
0
n
σ
-
,
n
κ
+
1
).For substepk
=
1,(35)dσ-n+1/m,κn+1/m,λn+1/mdΔε=α1APDen+1/m.Combining (34) and (35), we get the following equation according to the recursive relation: (36)
d
σ
-
n
+
1
,
κ
n
+
1
,
λ
n
+
1
d
Δ
ε
=
A
n
+
1
A
n
+
m
-
1
/
m
⋯
α
2
P
D
e
+
α
1
A
P
D
e
n
+
1
/
m
⋯
+
α
m
-
1
P
D
e
+
α
m
P
D
e
=
Ξ
σ
-
Ξ
κ
•
,where uninterested submatrix in (36) is represented by bold dots •. Obviously, d
σ
-
n
+
1
/
d
Δ
ε
=
Ξ
σ
- and d
κ
n
+
1
/
d
Δ
ε
=
Ξ
κ.


## 5. Solution of Global Problem

In the failure analysis of concrete material, convergence is a challenging problem in both stress update as mentioned above and capture of global responses. The traditional force control method finds it difficult to determine the exact load limit and is unable to trace the softening curve. The displacement control method, performing better than force control, has a strong dependency on selecting the proper loading freedoms. The classical arc-length method, which is able to trace “snap-through” and “snap-back” phenomena, is poor in localization problems of concrete. To help convergence, local arc-length method was proposed, but it requires the failure pattern of concrete to be a priori predictable. Gutiérrez [21], Verhoose et al. [22], and van der Meer et al. [23] developed a new arc-length method whose constraint path is constructed on the energy dissipation of a system. The method was proven to have better applicability than other methods in analyzing brittle/ductile failure. However, all the existing dissipation-based arc-length methods depend on simple constitutive models such as elastic damage model and elastic-plastic model. In this section, we have derived a constraint path, considering the combined contribution of plasticity and damage to the amount of released energy.The basic function by use of the arc-length method in global solution of a domainΩ is(37)fintugu,λF=λFf^0,where f
i
n
t
u is the equivalent internal force; u is the displacement vector; λ
F is the external load factor; f
^ is a prescribed loading pattern; g
u
,
λ
F is the constraint path.


### 5.1. Energy Release Rate

The energy release rate ofΩ can be computed as(38)G=P-V˙with(39)P=u˙Tfext=λFu˙Tf^(40)V=12∫ΩεeTσdΩ=12∫ΩεTσdΩ︷V1-12∫ΩεpTσdΩ︷V2,where P and V
˙ are the exerted power and rate of elastic energy, respectively; f
e
x
t
=
λ
F
f
^ is the vector of applied external forces; u
˙ is the nodal velocity.Differentiating (40) to time, for the item V
1, we can obtain(41)V˙1=12λ˙FuTf^+12λFu˙Tf^.The rate of the second itemV
2 can be derived as follows:(42)V˙2=12∫Ωε˙pTσ+εpTσ˙dΩ=12∫Ωε˙TΛTσ+σ˙TεpdΩ=12u˙T∫ΩBTΛTσ+CTεpdΩ=12u˙Tf∗,where Λ
=
d
ε
p
/
d
ε; C
=
d
σ
/
d
ε which can be obtained from (31).Equations (41) and (42) are combined to obtain(43)V˙=12λ˙FuTf^+12λFu˙Tf^-12u˙Tf∗.Substituting (39) and (43) into (38) yields(44)G=12u˙TλFf^+f∗-λ˙FuTf^.A forward Euler discretization is used to obtain the incremental path-following constraint as a function of the path parameterΔ
τ:(45)g=12ΔuTλ0Ff^+f0∗-ΔλFu0Tf^-Δτ,where Δ
u and Δ
λ
F are the displacement increment vector and load increment factor with regard to current load step, respectively; u
0, λ
0
F, and f
0
∗ are the converged displacement vector, converged load factor, and f
∗. The dissipation increment is illustrated as the shaded area in Figure 6.Figure 6
Energy dissipation increment for elastoplastic damage model under uniaxial loading.


### 5.2. Solution Scheme

The solution algorithm can be applied with the help of Sherman-Morrison formula, where the unknowns to be solved at iteration stepk are(46)dudλF=Δuk+1ΔλFk+1-ΔukΔλFk=dI-gk-1vTdII-wvTdI+gkdII-vTdI-gk1+vTdII-wwith(47)KdI=rk,KdII=-f^,K=∂fint∂u,r=λFf^-fint,v=∂g∂u,w=∂g∂λF.According to (45), we have(48)∂g∂u=12λF0f^+f∗0,∂g∂λF=-12u0Tf^.In our study, the global solution scheme is carried out in a hybrid scheme, as shown as follows.The flowchart of the global solution scheme is presented as follows:(1)
Initialization:λ
F
=
0; set n
d, Δ
λ
F, Δ
τ
m
i
n, and Δ
τ
m
a
x.
(2)
Solve a load step fromλ
0
F to λ
F (using force control):
(a)
Solvef
i
n
t
=
λ
F
f
^ by using Newton-Raphson method and record the iteration numbers n.
(b)
Checkn
>
n
d. If not, do λ
0
F
=
λ
F and λ
F
=
λ
F
+
Δ
λ
F; then go to the next load step (2(a)); else, do the following.
(c)
Compute the released energyG according to (44).
(d)
ComputeΔ
τ
=
0.5
r
·
G, with r
=
0.25
·
(
n
-
n
d
).
(e)
CheckΔ
τ
m
i
n
≤
Δ
τ
≤
Δ
τ
m
a
x. If not, do Δ
τ
=
Δ
τ
m
i
n or Δ
τ
=
Δ
τ
m
a
x; then go to step (3).
(3)
Switch to dissipation-based arc-length method:
(a)
UpdateK, f
i
n
t from local computations.
(b)
Updatef
∗, v, and w.
(c)
Start iteration:(49)r=λFf^-fintdI=K-1·r,dII=-K-1·f^.
Computed
u and d
λ
F according to (46):(50)u=u+du,λF=λF+dλF.
(d)
Check convergence. If not, go back to (3(a)); else, record iteration numbersn.
(e)
ComputeΔ
τ
=
0.5
r
·
G, with r
=
0.25
·
(
n
-
n
d
).
(f)
Updateλ
0
F
=
λ
F; then go to next load step (3(a)).In the above, Δ
τ
m
i
n and Δ
τ
m
a
x are thresholds that control the step length. The global solution is regarded as converged once the residual of the load satisfies(51)rkλFkf^≤TOLGlobal.


## 5.1. Energy Release Rate

The energy release rate ofΩ can be computed as(38)G=P-V˙with(39)P=u˙Tfext=λFu˙Tf^(40)V=12∫ΩεeTσdΩ=12∫ΩεTσdΩ︷V1-12∫ΩεpTσdΩ︷V2,where P and V
˙ are the exerted power and rate of elastic energy, respectively; f
e
x
t
=
λ
F
f
^ is the vector of applied external forces; u
˙ is the nodal velocity.Differentiating (40) to time, for the item V
1, we can obtain(41)V˙1=12λ˙FuTf^+12λFu˙Tf^.The rate of the second itemV
2 can be derived as follows:(42)V˙2=12∫Ωε˙pTσ+εpTσ˙dΩ=12∫Ωε˙TΛTσ+σ˙TεpdΩ=12u˙T∫ΩBTΛTσ+CTεpdΩ=12u˙Tf∗,where Λ
=
d
ε
p
/
d
ε; C
=
d
σ
/
d
ε which can be obtained from (31).Equations (41) and (42) are combined to obtain(43)V˙=12λ˙FuTf^+12λFu˙Tf^-12u˙Tf∗.Substituting (39) and (43) into (38) yields(44)G=12u˙TλFf^+f∗-λ˙FuTf^.A forward Euler discretization is used to obtain the incremental path-following constraint as a function of the path parameterΔ
τ:(45)g=12ΔuTλ0Ff^+f0∗-ΔλFu0Tf^-Δτ,where Δ
u and Δ
λ
F are the displacement increment vector and load increment factor with regard to current load step, respectively; u
0, λ
0
F, and f
0
∗ are the converged displacement vector, converged load factor, and f
∗. The dissipation increment is illustrated as the shaded area in Figure 6.Figure 6
Energy dissipation increment for elastoplastic damage model under uniaxial loading.


## 5.2. Solution Scheme

The solution algorithm can be applied with the help of Sherman-Morrison formula, where the unknowns to be solved at iteration stepk are(46)dudλF=Δuk+1ΔλFk+1-ΔukΔλFk=dI-gk-1vTdII-wvTdI+gkdII-vTdI-gk1+vTdII-wwith(47)KdI=rk,KdII=-f^,K=∂fint∂u,r=λFf^-fint,v=∂g∂u,w=∂g∂λF.According to (45), we have(48)∂g∂u=12λF0f^+f∗0,∂g∂λF=-12u0Tf^.In our study, the global solution scheme is carried out in a hybrid scheme, as shown as follows.The flowchart of the global solution scheme is presented as follows:(1)
Initialization:λ
F
=
0; set n
d, Δ
λ
F, Δ
τ
m
i
n, and Δ
τ
m
a
x.
(2)
Solve a load step fromλ
0
F to λ
F (using force control):
(a)
Solvef
i
n
t
=
λ
F
f
^ by using Newton-Raphson method and record the iteration numbers n.
(b)
Checkn
>
n
d. If not, do λ
0
F
=
λ
F and λ
F
=
λ
F
+
Δ
λ
F; then go to the next load step (2(a)); else, do the following.
(c)
Compute the released energyG according to (44).
(d)
ComputeΔ
τ
=
0.5
r
·
G, with r
=
0.25
·
(
n
-
n
d
).
(e)
CheckΔ
τ
m
i
n
≤
Δ
τ
≤
Δ
τ
m
a
x. If not, do Δ
τ
=
Δ
τ
m
i
n or Δ
τ
=
Δ
τ
m
a
x; then go to step (3).
(3)
Switch to dissipation-based arc-length method:
(a)
UpdateK, f
i
n
t from local computations.
(b)
Updatef
∗, v, and w.
(c)
Start iteration:(49)r=λFf^-fintdI=K-1·r,dII=-K-1·f^.
Computed
u and d
λ
F according to (46):(50)u=u+du,λF=λF+dλF.
(d)
Check convergence. If not, go back to (3(a)); else, record iteration numbersn.
(e)
ComputeΔ
τ
=
0.5
r
·
G, with r
=
0.25
·
(
n
-
n
d
).
(f)
Updateλ
0
F
=
λ
F; then go to next load step (3(a)).In the above, Δ
τ
m
i
n and Δ
τ
m
a
x are thresholds that control the step length. The global solution is regarded as converged once the residual of the load satisfies(51)rkλFkf^≤TOLGlobal.


## 6. Model Validation

To validate the applicability and effectiveness of the model, a cubic plain concrete specimen with size of 100 × 100 × 100 (unit: mm) is simulated and computed under different loading conditions. The numerical results are compared with experiment data. Material parameters for the yield function and the plastic potential are listed in Table2. The convergence tolerance values for the local and global solutions are set as T
O
L
L
o
c
a
l
=
T
O
L
G
l
o
b
a
l
=
0.1%.Table 2
Material parameters for the yield function and plastic potential.
θ
γ
ρ
φ
1.15
3.0
0.1
π/6


### 6.1. Uniaxial Loading

Two experiments [30, 31] are simulated to investigate the applicability of the model in describing the uniaxial tensile and compressive loading behaviors of concrete. Corresponding parameters are listed in Tables 3 and 4, respectively. As shown in Figure 7, it is seen that the computed stress-strain curve agrees well with experimental data in both ascending and descending branches, indicating that the proposed model is able to capture the uniaxial hardening/softening properties of concrete in both tension and compression.Table 3
Material parameters for uniaxial tension.
E
0
t (MPa)
v
ε
0
t
ε
f
t
f
f
t (MPa)
d
0
t
A
2
t
B
3
t
Δ
λ
F
n
d
Δ
τ
m
i
n (N⋅m)
Δ
τ
m
a
x (N⋅m)
3.1
×
10
4
0.18
0.000105
0.00012
3.48
0.0
0.7f
f
t
2
1.5
1.0
3
0.0001
1.0Table 4
Material parameters for uniaxial compression.
E
0
c (MPa)
v
ε
0
c
ε
f
c
f
f
c (MPa)
d
0
c
A
2
c
B
3
c
Δ
λ
F
n
d
Δ
τ
m
i
n (N⋅m)
Δ
τ
m
a
x (N⋅m)
3
.
1
×
10
4
0.18
0.00067
0.002
-
27.6
0.0
0.75
2.0
1.0
3
0.0001
5.0Figure 7
Numerical results compared with experimental data: (a) uniaxial tension and (b) uniaxial compression.
(a)
(b)


### 6.2. Biaxial Loading

The same parameters in Tables3 and 4 are used to observe the performance of the proposed model under biaxial loadings [32], with exception of f
f
c. In Figure 8, the comparisons of numerical results obtained from different biaxial stress ratios with experimental data are shown. All the values of σ
1 are regularized by β
p
=
f
f
c, with f
f
c taking the values of −29.5 MPa and −32.8 MPa for biaxial tensile and compressive loadings, respectively. Good agreements are seen in these two cases.Figure 8
Numerical results compared with experimental data: (a) biaxial tension and (b) biaxial compression.
(a)
(b)


## 6.1. Uniaxial Loading

Two experiments [30, 31] are simulated to investigate the applicability of the model in describing the uniaxial tensile and compressive loading behaviors of concrete. Corresponding parameters are listed in Tables 3 and 4, respectively. As shown in Figure 7, it is seen that the computed stress-strain curve agrees well with experimental data in both ascending and descending branches, indicating that the proposed model is able to capture the uniaxial hardening/softening properties of concrete in both tension and compression.Table 3
Material parameters for uniaxial tension.
E
0
t (MPa)
v
ε
0
t
ε
f
t
f
f
t (MPa)
d
0
t
A
2
t
B
3
t
Δ
λ
F
n
d
Δ
τ
m
i
n (N⋅m)
Δ
τ
m
a
x (N⋅m)
3.1
×
10
4
0.18
0.000105
0.00012
3.48
0.0
0.7f
f
t
2
1.5
1.0
3
0.0001
1.0Table 4
Material parameters for uniaxial compression.
E
0
c (MPa)
v
ε
0
c
ε
f
c
f
f
c (MPa)
d
0
c
A
2
c
B
3
c
Δ
λ
F
n
d
Δ
τ
m
i
n (N⋅m)
Δ
τ
m
a
x (N⋅m)
3
.
1
×
10
4
0.18
0.00067
0.002
-
27.6
0.0
0.75
2.0
1.0
3
0.0001
5.0Figure 7
Numerical results compared with experimental data: (a) uniaxial tension and (b) uniaxial compression.
(a)
(b)


## 6.2. Biaxial Loading

The same parameters in Tables3 and 4 are used to observe the performance of the proposed model under biaxial loadings [32], with exception of f
f
c. In Figure 8, the comparisons of numerical results obtained from different biaxial stress ratios with experimental data are shown. All the values of σ
1 are regularized by β
p
=
f
f
c, with f
f
c taking the values of −29.5 MPa and −32.8 MPa for biaxial tensile and compressive loadings, respectively. Good agreements are seen in these two cases.Figure 8
Numerical results compared with experimental data: (a) biaxial tension and (b) biaxial compression.
(a)
(b)


## 7. Conclusion

This paper introduces a numerical framework of implementing the spatial elastoplastic damage model of concrete via the recently developed IGA method. For the material model, a nonassociated flow rule is considered to describe the plasticity of concrete, where a three-parameter Barcelona yield surface and a modified Drucker-Prager plastic potential are used. The damage of concrete is considered as a combination of tensile and compressive damage variables, with the evolution of each driven by the internal variables and expressed by a piecewise function. The main procedure is divided into two kinds of solutions: local and global. For the local kind, the decoupling of plasticity and damage is achieved via the philosophy of operator split and both return-mapping algorithm and substepping strategy are used to help convergence. For the global kind, a dissipation-based arc-length method with constraint path is developed, considering contribution of plasticity and damage to the released amount of energy. Through comparisons between numerical simulations and experimental data, the approach thus employed has been proven to be reliable in reflecting the concrete properties and can therefore be further applied to the advanced analysis of concrete structures.



---

*Source: 4273024-2016-06-22.xml*