# Single Spring Joint Element Based on the Mixed Coordinate System

Authors: Lanhao Zhao; Wei Zhang; Xin Bai; Tianyou Yan; Tongchun Li
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

As a FEM for reinforced concrete bond-slip problems, one important feature of the typical double spring joint element method is the selection of the normal stiffness, which may cause the mutual embedding problem and bring challenges to the calculation. In this paper, a novel single spring joint element method based on the mixed coordinate system is proposed to simulate the interaction of two materials. Instead of choosing the normal stiffness arbitrarily, the proposed method makes DOFs of two materials in the normal direction equal to ensure deformation compatibility. And its solid elements for the concrete are solved in global coordinate system, while the beam elements for the steel bar are solved in local coordinate system. In addition, the proposed method can also be applied to RC structures with irregular arrangements of steel bars. Numerical examples demonstrate the validity and accuracy of the proposed approach. Furthermore, the bond model is applied to RC beams with the description of the damage process.


---

## Body

## 1. Introduction

Reinforced concrete (RC) structures are widely used for engineering and construction purposes. In certain circumstances, knowledge about the damage behavior of RC structures can be essential for real applications. For example, in the case of RC tower for wind turbines [1] or nuclear reactor containment [2], predicting the damage area of RC structures is one of the key points to safety and security.Finite element method (FEM) has become the most commonly adopted method to predict the strength and the deformation of composites, namely, concrete and reinforcing bar. Typically, based on the methods used, nonlinear FEM of RC can be divided into three categories:discrete model,embedded model, andsmeared model [3–8]. In the discrete model (Figure 1(a)), a joint element is inserted to simulate the bond-slip performance of the concrete and the reinforcing bar. However, in the embedded model and smeared model, two materials are considered to bond together completely [9], so these two models do not focus on the description of the influence of the bond stress. To account for the joint element in the discrete model, several different approaches can be used indouble spring element (Figure 4(a)),quadrilateral bond element without thickness,bond zone element, and so forth [10–12]. Double spring joint element method with zero thickness is proposed firstly by Ngo and Scordelis [13], while bond zone element method (Figure 1(b)) is firstly proposed by de Groot et al. [14], assuming that there is a cylinder slip layer composed of a series of parallel corpus vertebrae with thickness outside the steel bar, to simulate softened layer after internal cracks in concrete, reflecting wedge action of deformed bars.Figure 1
(a) Discrete model of nonlinear FEM for bond-slip problems (e.g., a joint element with zero thickness is inserted to the interface between two materials, which has elastic stiffness but no physical dimensions; i.e., it is not a real component but a mechanical model). (b) Bond zone element (hypothetical bond-slip layer with thicknessβ
γ, composed of a series of parallel corpus vertebrae outside the steel bar to simulate softened layer after internal cracks in concrete, reflecting wedge action of deformed bars).Recently, improved methods have attracted increasing research interest, which can serve as alternatives to the above-mentioned methods for the analysis of RC members. Among these methods, one approach is to improve the embedded model with consideration of bond-slip effect. Ibrahimbegovic et al. [15], among others, developed embedded elements and described the bond-slip effect through an enrichment of the degrees of freedom. Other hybrid approaches also predict the mechanical behaviors of RC. For example, Wang et al. [16] proposed the FEM combined with the hybrid boundary node method (HdBNM) as a truly boundary-type meshless method with only discrete nodes distributed on the faces of the crack. Ismail et al. [17] studied the mechanical behavior of unidirectional fibre-reinforced polymer composites subjected to transverse tension using a two-dimensional discrete element method (DEM). This method, unlike the conventional FEM, is convenient to deal with local behavior of a material by defining local models or parameters for the specified particles and contacts. And both of these two new methods are using double spring element to simulate the bond stress. Furthermore, in order to take the steel strain into account in the bond-slip relationship, Wu and Gilbert [18] developed and implemented a modified CEB-FIP bond model (Figure 2) considering the effects of the concrete damage due to primary cracking and bond degradation under increasing steel stress and incorporated it into an existing FEM of reinforced concrete. Meanwhile, Santos and Henriques [19] presented a new bond element with modified material constitutive laws for finite element packages, which consists of an orthotropic four-node plane stress element.Figure 2
CEB-FIP bond model: at service load levels, when the steel stress-strain relationship is linear-elastic, slip normally does not exceedS
1. The basic bond-slip relationship (ignoring cracking and steel stress state) is therefore based on the ascending part of CEB-FIP bond model, which can be expressed as in the figure.Although these models can give appropriate solutions to mechanical analysis of certain bond-slip problems, there are also some drawbacks. For example, quadrilateral bond element without thickness is difficult to implement in 3D structures. Bond zone element has a heavy computational cost due to its actual size and also increases the level of difficulties in the meshing process. In studies where the overall structural behavior is of primary interest, bond zone elements may not provide a reasonable balance between accuracy and computational efficiency [11].Notably, when double spring element method is used to describe the interaction between the steel bar and the concrete, two parameters need to be predetermined, namely, normal stiffnessK
V and tangential stiffness K
H. K
H is determined by the relationship between bond stress τ and slip s, getting from the experimental data, such as pull-out test [20]. Until now, there is still a lack of research on the value of K
V. To ensure the deformation coordination and reflect the interaction between the two materials, a large enough value in arbitrariness choice is usually adopted to keep a small enough penetration between the two materials. However, overshoot values will bring challenges to the numerical calculation, while oversmall values will lead to mutual embedding problem between the steel elements and the concrete elements.In this paper, a novel single spring joint element method based on mixed coordinate system is proposed. In this method [21], the difficulty of choosing the normal stiffness coefficient has been avoided. Firstly, a pair of nodes are set to the interface between the two materials. Then, a single spring joint element (Figures 3 and 4(b)) is established between the two nodes. It is worth noting that the single spring is set along the tangent direction of the reinforcing bar to simulate the tangential interaction between two materials with its tangential stiffness determined by the relationship of the bond force and slip, while, in the normal direction, DOFs between two nodes are forced to be equal to ensure the deformation compatibility of two materials. The solid elements for the concrete are solved in a global coordinate system, while the beam elements for the steel bar are solved in the local coordinate system to establish the relation of DOFs in the normal direction. This method can also be used in beam elements of the steel bar with irregular arrangements.Figure 3
Single spring joint method model: interpolation nodej interpolated from solid element for concrete has the same coordinates with steel node i and single spring is set along the tangent direction of steel bar between node i and node j to simulate the tangential interaction while DOFs in normal direction are forced to be equal.Figure 4
Joint element model for plane problem: (a) typical double spring joint element, (b) single spring joint element.This paper is organized as follows: in Section2, typical double spring joint element method is introduced and the theoretical formulation of single spring joint element method is discussed in Section 3, followed by the damage model used for numerical simulation in Section 4. Validation cases including an unbonded prestressed reinforced concrete and a Shima’s #10 specimen are presented in Section 5, which proves the validity and accuracy of the proposed numerical scheme. Finally, the proposed single spring joint method is applied to analyze isotropic damage of reinforced concrete with slender beams.


## 2. Review of Typical Double Spring Joint Element Method

The typical double spring joint element method firstly proposed by Ngo and Scordelis [13] is composed of two orthogonal springs for plane problem (Figure 4(a)) or three orthogonal springs for space problems. They can pass the extrusion force and shear force between one node (node i in Figure 3) in the solid element for concrete and an adjacent node (node j in Figure 3) in the bar element for steel. The joint element is not a physical component but a mechanical model, which has elastic stiffness but no physical dimensions; that is, the two connected nodes have the same coordinates. The form of this model is simple and it is widely used in finite element simulation of the bond-slip relationship between the steel bar and the concrete.As shown in Figure4, the local coordinate system H
-
V along the reinforcing bar is established. At the same time, spring stiffness is defined in two directions in the local coordinate system H
-
V: normal stiffness K
V and tangential stiffness K
H. Using the principle of virtual work, the finite element balance equation is established between the nodal force of the double spring joint element and node displacement [12, 22]:(1)
F
e
=
T
T
K
H
0
0
K
V
T
δ
e
,where T is the coordinate transformation matrix between the global coordinate system X
-
Y and the local coordinate system H
-
V. In the reinforced concrete structures, the relationship between the steel and the concrete is governed by a bond stress distributed along the interface. The evolution of this bond stress is generally studied in the form of a bond stress-relative displacement (slip between the steel bar and the concrete) curve. Therefore, the bond behavior can be represented by a bond stress-slip curve with K
H standing for the tangential stiffness obtained from the slope of bond-slip curve, that is, d
τ
/
d
s
A. Equation (1) includes force instead of stress, so tangential stiffness is multiplied by area A. Generally, a large enough value for K
V is arbitrarily chosen to ensure the penetration of the steel and the concrete is small enough.


## 3. Proposed Analytical Model

### 3.1. General Description

As mentioned previously, the form of typical double spring element is simple and its mechanics concept is clear. In the cases where the overall structural behavior is of primary interest, the double spring element provides a reasonable balance between accuracy and computational efficiency and is still widely used in finite element simulation of the bond-slip relationship between steel and concrete. However, the use of the typical double spring element in the analysis of RC structures imposes the following restrictions:(1)
There is certain arbitrariness in the choice of normal stiffness, which varies from case to case.
(2)
Overshoot value will bring challenges to calculation, while too small value will lead to mutual embedding problem in elements between steel and concrete.To address some of these limitations, an analytical model-single spring joint element method is proposed in this study. Single spring joint element based on mixed coordinate system has only one spring set along the tangent direction of steel bar between nodei and node j, while degree of freedoms in normal direction are forced to be equal. The displacement of hypothetical node j is obtained by interpolation and then passes to the corresponding node i of the steel bar. With the mesh of RC structure, the solid elements for the concrete are solved in the global coordinate system, while the beam elements of the steel bar are solved in its local coordinate system. The displacement and the force can be expressed in different coordinate systems by coordinate transformation. Then, the establishment of finite element equation for single spring joint element method, related to normal degree of freedom constraint equations and interaction of the concrete and the steel, is specifically addressed.


### 3.2. Coordinate Transformation

Generally, stiffness matrix of the beam or bar element is calculated in its local coordinate system. To establish the whole balance equation, we customarily transform degrees of freedom in beam or bar element in the local coordinate system to global coordinate system. The solution of finite element equation does not depend on the choice of coordinate system. Different parts of the computing system can choose different reference coordinate system.As shown in Figure5, three coordinate systems are defined: global coordinate system x
y
z, local coordinate system x
′
y
′
z
′ for beam element, nodal local coordinate system x
∗
y
∗
z
∗ for steel node in beam element. In the coordinate system x
′
y
′
z
′, x
′ is arranged in axial direction of beam element. Considering a steel bar may not necessarily form a straight line in practical problems, coordinate system x
∗
y
∗
z
∗ is set up in each steel node in beam element. Node A (Figure 5) is the intersection of beam element 1 and beam element 2. However, these two beams have different axial directions. So axial direction of node A is taken as the mean axial direction of beam element 1 and beam element 2 using vector addition method, avoiding the difficulty of analyzing the axial direction of these nodes.Figure 5
Mixed coordinate system for single spring joint element (R represents coordinate transformation matrix between x
y
z and x
′
y
′
z
′, while r represents coordinate transformation matrix between x
y
z and x
∗
y
∗
z
∗. Node A is the intersection of the beam element 1 and beam element 2).Transformation relationship between the displacementu, u
′, u
∗ and the force F, F
′, F
∗ in three coordinate systems is established using the following:(2)
u
′
=
R
u
,
F
′
=
R
F
,
(3)
u
∗
=
r
u
,
F
∗
=
r
F
,where R is the coordinate transformation matrix between x
y
z and x
′
y
′
z
′ and r is the coordinate transformation matrix between x
y
z and x
∗
y
∗
z
∗.Usually, for space problems rotating around thez-axis, R can be written in matrix form as the following, where θ is the rotation angle:(4)
R
=
cos
⁡
θ
sin
⁡
θ
0
-
sin
⁡
θ
cos
⁡
θ
0
0
0
1
.Based on the FEM theory, finite element balance equation of beam element in local coordinate system can be established as(5)
K
′
u
′
=
F
′
.Accordingly, finite element balance equation in coordinate systemx
∗
y
∗
z
∗ is given by transformation equations (2) and (3) and substituted into FEM equation (5):(6)
K
∗
u
∗
=
F
∗
,where K
∗
=
r
K
r
T
=
r
R
T
K
′
R
r
T and K
∗ is defined as the stiffness matrix in the coordinate system x
∗
y
∗
z
∗.At this point, we establish contacts between three coordinate systems. Satisfactorily, stiffness matrixes in coordinate systemx
∗
y
∗
z
∗ can be expressed by variables in the other coordinate system.


### 3.3. Interpolation and Normal DOF Constraint Equations

In order to establish the relationship between the steel bar and the concrete in the normal direction, we also need to set a pair of nodes where the two materials intersect (Figures3 and 6). Unlike the double spring joint element method, there is no need to define a normal spring. Using the principle of compatible displacement method, the pair of nodal degrees of freedom is set to be equal to each other in the normal direction, and the relationship between the steel bar and the concrete in the normal direction is established.Figure 6
Steel bar in concrete with irregular arrangements.TakingX direction as the axial direction of steel bar and Y direction and Z direction as the normal direction, the normal displacement of the interpolation node in the concrete can be obtained as (3), following (7)
u
c
y
∗
=
∑
j
=
1
n
r
y
j
u
c
j
,
u
c
z
∗
=
∑
j
=
1
n
r
z
j
u
c
j
,where u
c is the displacement of surrounding concrete node in solid element in the global coordinate, u
c
∗ is the displacement of the interpolation node in Figure 6, r
j represents the interpolated coefficient, and r
j is the element of coordinate transform matrix. n represents two-dimensional or three-dimensional problems.For the problems with larger number and densely distributed bars, keeping all pairs of nodes’ positions coincident with the concrete elements and steel elements may introduce additional difficulties in meshing process. The method proposed in this paper can also be used in beam elements of the steel bar with irregular arrangements. Therefore, the nodes belonging to the concrete can be expressed by linear combinations of the surrounding concrete nodes (i.e., there are eight nodes in Figure6):(8)
u
c
j
=
∑
k
=
1
8
N
k
u
c
j
k
,where N
k is the shape functions of finite element method.K
H can relate the degrees of freedom of the concrete and the steel in the tangential direction of the reinforcement. To close the system, a perfect relation is imposed for the displacement in the other directions. Therefore, interpolation equation of normal displacement for the steel bar can be established as follows:(9)
u
s
y
∗
=
u
c
y
∗
=
∑
j
=
1
n
∑
k
=
1
8
r
y
j
N
k
u
c
j
k
,
u
s
z
∗
=
u
c
z
∗
=
∑
j
=
1
n
∑
k
=
1
8
r
z
j
N
k
u
c
j
k
,where u
s
∗ is the displacement of the steel node which has the same coordinate with the interpolation node.


### 3.4. Interaction of Concrete and Steel

Definef
∗ as the force applied by the concrete on the steel bar in nodal local coordinate x
∗
y
∗
z
∗; according to Newton’s third law, we know that -
r
T
f
∗ is the force applied by the steel bar on the concrete in the global coordinate x
y
z; reinforced concrete incremental finite element balance equations can be written, for concrete as(10)
K
Δ
u
c
=
Δ
F
-
r
T
Δ
f
∗
,and, correspondingly, for steel bar as(11)
K
∗
Δ
u
s
∗
=
Δ
F
∗
+
Δ
f
∗
.Due to the fact that the normal freedom of nodes is forced to be equal, normal force between the concrete and the steel bar has become internal force and will not be present in the balance equation; then,f
∗ is only the tangential bond force between the concrete and the steel bar:(12)
Δ
f
∗
=
K
s
d
s
=
K
s
r
Δ
u
c
-
Δ
u
s
∗
,
(13)
k
s
=
∂
τ
∂
s
π
D
l
d
s
,
K
s
=
k
s
0
0
0
0
0
0
0
0
,where D represents the diameter of steel bar, l is the length of steel bar and d
s is slip of steel bar, and τ
-
s represents the bond-slip curve obtained from experiments, and, for example, a nonlinear bond-slip constitutive relationship is represented below, proposed by Shima:(14)
τ
ε
s
,
s
=
0.73
f
c
′
ln
⁡
1
+
5000
s
D
3
1
1
+
1
0
5
ε
s
,where ε
s is the steel strain, f
c
′ is the compressive strength of concrete, D is the steel diameter and s is the slip, and τ is the bond stress.Finally, by substituting (12) into (10) and (11) and in matrix form,(15)
K
+
r
T
K
s
r
-
r
T
K
s
-
K
s
r
K
∗
+
K
s
Δ
u
c
Δ
u
s
∗
=
Δ
F
Δ
F
∗
.Equation (15) is the final finite element equation of single spring joint element method based on the mixed coordinate system.


## 3.1. General Description

As mentioned previously, the form of typical double spring element is simple and its mechanics concept is clear. In the cases where the overall structural behavior is of primary interest, the double spring element provides a reasonable balance between accuracy and computational efficiency and is still widely used in finite element simulation of the bond-slip relationship between steel and concrete. However, the use of the typical double spring element in the analysis of RC structures imposes the following restrictions:(1)
There is certain arbitrariness in the choice of normal stiffness, which varies from case to case.
(2)
Overshoot value will bring challenges to calculation, while too small value will lead to mutual embedding problem in elements between steel and concrete.To address some of these limitations, an analytical model-single spring joint element method is proposed in this study. Single spring joint element based on mixed coordinate system has only one spring set along the tangent direction of steel bar between nodei and node j, while degree of freedoms in normal direction are forced to be equal. The displacement of hypothetical node j is obtained by interpolation and then passes to the corresponding node i of the steel bar. With the mesh of RC structure, the solid elements for the concrete are solved in the global coordinate system, while the beam elements of the steel bar are solved in its local coordinate system. The displacement and the force can be expressed in different coordinate systems by coordinate transformation. Then, the establishment of finite element equation for single spring joint element method, related to normal degree of freedom constraint equations and interaction of the concrete and the steel, is specifically addressed.


## 3.2. Coordinate Transformation

Generally, stiffness matrix of the beam or bar element is calculated in its local coordinate system. To establish the whole balance equation, we customarily transform degrees of freedom in beam or bar element in the local coordinate system to global coordinate system. The solution of finite element equation does not depend on the choice of coordinate system. Different parts of the computing system can choose different reference coordinate system.As shown in Figure5, three coordinate systems are defined: global coordinate system x
y
z, local coordinate system x
′
y
′
z
′ for beam element, nodal local coordinate system x
∗
y
∗
z
∗ for steel node in beam element. In the coordinate system x
′
y
′
z
′, x
′ is arranged in axial direction of beam element. Considering a steel bar may not necessarily form a straight line in practical problems, coordinate system x
∗
y
∗
z
∗ is set up in each steel node in beam element. Node A (Figure 5) is the intersection of beam element 1 and beam element 2. However, these two beams have different axial directions. So axial direction of node A is taken as the mean axial direction of beam element 1 and beam element 2 using vector addition method, avoiding the difficulty of analyzing the axial direction of these nodes.Figure 5
Mixed coordinate system for single spring joint element (R represents coordinate transformation matrix between x
y
z and x
′
y
′
z
′, while r represents coordinate transformation matrix between x
y
z and x
∗
y
∗
z
∗. Node A is the intersection of the beam element 1 and beam element 2).Transformation relationship between the displacementu, u
′, u
∗ and the force F, F
′, F
∗ in three coordinate systems is established using the following:(2)
u
′
=
R
u
,
F
′
=
R
F
,
(3)
u
∗
=
r
u
,
F
∗
=
r
F
,where R is the coordinate transformation matrix between x
y
z and x
′
y
′
z
′ and r is the coordinate transformation matrix between x
y
z and x
∗
y
∗
z
∗.Usually, for space problems rotating around thez-axis, R can be written in matrix form as the following, where θ is the rotation angle:(4)
R
=
cos
⁡
θ
sin
⁡
θ
0
-
sin
⁡
θ
cos
⁡
θ
0
0
0
1
.Based on the FEM theory, finite element balance equation of beam element in local coordinate system can be established as(5)
K
′
u
′
=
F
′
.Accordingly, finite element balance equation in coordinate systemx
∗
y
∗
z
∗ is given by transformation equations (2) and (3) and substituted into FEM equation (5):(6)
K
∗
u
∗
=
F
∗
,where K
∗
=
r
K
r
T
=
r
R
T
K
′
R
r
T and K
∗ is defined as the stiffness matrix in the coordinate system x
∗
y
∗
z
∗.At this point, we establish contacts between three coordinate systems. Satisfactorily, stiffness matrixes in coordinate systemx
∗
y
∗
z
∗ can be expressed by variables in the other coordinate system.


## 3.3. Interpolation and Normal DOF Constraint Equations

In order to establish the relationship between the steel bar and the concrete in the normal direction, we also need to set a pair of nodes where the two materials intersect (Figures3 and 6). Unlike the double spring joint element method, there is no need to define a normal spring. Using the principle of compatible displacement method, the pair of nodal degrees of freedom is set to be equal to each other in the normal direction, and the relationship between the steel bar and the concrete in the normal direction is established.Figure 6
Steel bar in concrete with irregular arrangements.TakingX direction as the axial direction of steel bar and Y direction and Z direction as the normal direction, the normal displacement of the interpolation node in the concrete can be obtained as (3), following (7)
u
c
y
∗
=
∑
j
=
1
n
r
y
j
u
c
j
,
u
c
z
∗
=
∑
j
=
1
n
r
z
j
u
c
j
,where u
c is the displacement of surrounding concrete node in solid element in the global coordinate, u
c
∗ is the displacement of the interpolation node in Figure 6, r
j represents the interpolated coefficient, and r
j is the element of coordinate transform matrix. n represents two-dimensional or three-dimensional problems.For the problems with larger number and densely distributed bars, keeping all pairs of nodes’ positions coincident with the concrete elements and steel elements may introduce additional difficulties in meshing process. The method proposed in this paper can also be used in beam elements of the steel bar with irregular arrangements. Therefore, the nodes belonging to the concrete can be expressed by linear combinations of the surrounding concrete nodes (i.e., there are eight nodes in Figure6):(8)
u
c
j
=
∑
k
=
1
8
N
k
u
c
j
k
,where N
k is the shape functions of finite element method.K
H can relate the degrees of freedom of the concrete and the steel in the tangential direction of the reinforcement. To close the system, a perfect relation is imposed for the displacement in the other directions. Therefore, interpolation equation of normal displacement for the steel bar can be established as follows:(9)
u
s
y
∗
=
u
c
y
∗
=
∑
j
=
1
n
∑
k
=
1
8
r
y
j
N
k
u
c
j
k
,
u
s
z
∗
=
u
c
z
∗
=
∑
j
=
1
n
∑
k
=
1
8
r
z
j
N
k
u
c
j
k
,where u
s
∗ is the displacement of the steel node which has the same coordinate with the interpolation node.


## 3.4. Interaction of Concrete and Steel

Definef
∗ as the force applied by the concrete on the steel bar in nodal local coordinate x
∗
y
∗
z
∗; according to Newton’s third law, we know that -
r
T
f
∗ is the force applied by the steel bar on the concrete in the global coordinate x
y
z; reinforced concrete incremental finite element balance equations can be written, for concrete as(10)
K
Δ
u
c
=
Δ
F
-
r
T
Δ
f
∗
,and, correspondingly, for steel bar as(11)
K
∗
Δ
u
s
∗
=
Δ
F
∗
+
Δ
f
∗
.Due to the fact that the normal freedom of nodes is forced to be equal, normal force between the concrete and the steel bar has become internal force and will not be present in the balance equation; then,f
∗ is only the tangential bond force between the concrete and the steel bar:(12)
Δ
f
∗
=
K
s
d
s
=
K
s
r
Δ
u
c
-
Δ
u
s
∗
,
(13)
k
s
=
∂
τ
∂
s
π
D
l
d
s
,
K
s
=
k
s
0
0
0
0
0
0
0
0
,where D represents the diameter of steel bar, l is the length of steel bar and d
s is slip of steel bar, and τ
-
s represents the bond-slip curve obtained from experiments, and, for example, a nonlinear bond-slip constitutive relationship is represented below, proposed by Shima:(14)
τ
ε
s
,
s
=
0.73
f
c
′
ln
⁡
1
+
5000
s
D
3
1
1
+
1
0
5
ε
s
,where ε
s is the steel strain, f
c
′ is the compressive strength of concrete, D is the steel diameter and s is the slip, and τ is the bond stress.Finally, by substituting (12) into (10) and (11) and in matrix form,(15)
K
+
r
T
K
s
r
-
r
T
K
s
-
K
s
r
K
∗
+
K
s
Δ
u
c
Δ
u
s
∗
=
Δ
F
Δ
F
∗
.Equation (15) is the final finite element equation of single spring joint element method based on the mixed coordinate system.


## 4. Concrete Damage Model

A large number of damage models may be found in the literature [23–28] that the concrete exhibits significant strain softening beyond the peak or initial failure stress, and the failure criterion of concrete cannot be easily determined under complex stress states [29]. Four-parameter equivalent strain for the isotropic damage model which is modified Mazars [30] model incorporates this behavior and also reflects the entire failure process. This equivalent strain, based on the idea of Hsieh et al. [31] strength failure criterion in the stress space, can depict the mechanical behavior of the concrete under uniaxial stress states or multiaxial stress states. In this work, single spring joint method is applied to analyze the isotropic damage of reinforced concrete with this model.As we all know, the damage parameter can be calculated by the following function:(16)
d
=
1
-
σ
ε
,
ε
p
,
d
E
0
ε
,where d (ranging from 0 to 1) is a scalar value representing the local damage parameter, σ and ε are, respectively, the stress and the strain tensors, and E
0 is the modulus of elasticity.In four-parameter equivalent strain for the isotropic damage model, a four-parameter criterion in strain space is developed:(17)
F
I
1
′
,
J
2
′
,
ε
1
=
A
J
2
′
ε
0
+
B
J
2
′
+
C
ε
1
+
D
I
1
′
-
ε
0
,where ε
1, ε
2, and ε
3 express the three principal stresses, respectively. ε
0 is the peak strain corresponding to the strength limitation f
c under the uniaxial compression, ε
0
=
C
t
f
c
/
E
=
f
t
/
E, f
t is the tensile strength, f
c is the compressing strength, J
2
′ represents the second deviatoric strain invariant, I
1
′ represents the first spherical strain invariant, and A, B, C, and D are the material parameters determined by available experimental data. The criterion in strain space is also formulated by other researchers with a similar formula as used in the current work but the definition of equivalent strain is different.For the multiaxial stress, the strainε should be changed as the equivalent strain ε
∗. A four-parameter equivalent strain based on the basic idea of Hsieh-Ting-Chen strength failure criterion is proposed as follows:(18)
ε
∗
=
A
J
2
′
ε
0
+
B
J
2
′
+
C
ε
1
+
D
I
1
′
,where I
1
′
=
ε
i
i, J
2
′
=
e
i
j
e
i
j
/
2, ε
1
=
2
/
3
J
2
′
sin
⁡
(
θ
+
2
/
3
π
)
+
1
/
3
I
1
′, Lode’s angle θ
=
1
/
3
s
i
n
-
1
(
-
3
3
J
3
′
/
2
J
2
′
3
), and third deviatoric strain invariant J
3
′
=
e
i
j
e
j
k
e
k
l
/
3.The equivalent strain in all kinds of stress states can be solved from (18); then, substitute ε
∗ into (16), using Ghrib and Tinawi [32] damage model; the damage parameter can be solved:(19)
d
=
0
ε
∗
≤
ε
0
d
=
1
-
f
t
E
ε
∗
2
e
-
b
ε
∗
-
ε
0
-
e
-
2
b
ε
∗
-
ε
0
ε
∗
>
ε
0
,where b
=
3
/
ε
0
(
2
G
f
E
/
l
c
h
f
t
2
-
1
), G
f is the fracture energy and the commonly used formula for G
f is G
f
=
(
2.72
+
0.0214
f
t
)
f
t
2
d
m
a
x
/
E
0, l
c
h is the characteristic length and its expression is l
c
h
=
E
G
f
/
f
t
2, and d
m
a
x is the maximum particle size for concrete.Validation of this model will be demonstrated in Sections5.3 and 5.4.


## 5. Numerical Simulations

In this section, the method is applied to the numerical simulations of an unbonded prestressed reinforced concrete and a Shima’s #10 specimen. The first simulation is studied to verify the accuracy and capability of the proposed single spring joint method, and the latter one is to evaluate the influence of the bond effects. And the results are compared with existing published experimental and numerical data in literature. Then, the single spring joint method is applied to analyze isotropic damage of reinforced concrete with slender beams, including three- and four-point bending beams, showing the midspan deflection and the damage distribution of the slender concrete beams in comparison with experiments by other researchers [33–36].


### 5.1. Unbonded Prestressed RC Structure

A reinforced concrete block structure with U-shaped unbonded prestressed steel bar is considered. The radius of bending steel is equal to 2 m, and the bottom of the concrete block is fixed. The diameter of the reinforcing steel is 36 mm with material parameters listed in Table1 and its end is subjected to a concentrated load equal to 1000 kN. For this simulation, bend steel bar is meshed by beam elements. FEM model and its dimensions are illustrated in Figure 7 To highlight the effectiveness of the single spring joint method in describing the normal interaction of the concrete and steel bar; the bond stress between steel bar and concrete is not considered; that is, τ
=
0.Table 1
Material parameters for example 1.
Schedule
Properties
Level
Concrete
Elastic modulus
25.5 GPa
Poisson’s ratio
0.167
Reinforced steel
Elastic modulus
200 GPa
Poisson’s ratio
0.23Figure 7
Finite element method (FEM) model for example 1.Figure8(a) represents the distribution of force along the steel bar. The value is from 499.905 kN/m to 500.015 kN/m, which is close to the analytical solution 500 kN/m obtained from (20) and illustrated in Figure 8(b). Although distribution of force is close between the two methods for this simple test case, there is no normal interpenetration amount between concrete and steel with the single spring joint element method as compared with the typical double spring joint element method. The relationship between normal stiffness K
V and normal embedment amount is given in Table 2, and the order of normal embedment amount is inversely proportional to the order of K
V. When K
V is 1 kPa, the normal embedment amount is 0.1013 × 10−1 m. While K
V is 1 × 107 kPa, the embedment amount is reduced to 0.13 × 10−8 m, and when K
V continues to increase, the convergence deteriorates or even fails:(20)
∫
σ
d
r
=
∫
0
π
σ
r
d
θ
=
2
F
.Table 2
The relationship between normal stiffness and normal embedment amount.
Normal stiffnessK
V (KPa)
Normal embedment amount (m)
1
0.1013
×
10
-
1
1
×
10
2
0.1013
×
10
-
3
1
×
10
4
0.1013
×
10
-
5
1
×
10
7
0.1013
×
10
-
8
1
×
10
12
0.1013
×
10
-
13Figure 8
(a) Distribution of force along the steel bar. (b) The analytical solution of example 1.
(a)
(b)


### 5.2. Shima Tensile Specimen

In order to obtain bond-slip relationship suitable for different thicknesses of the concrete cover, anchor length, and different parameters of materials, Shima et al. [37], among others, did a series of pull tests and tensile tests of cylindrical specimens. With cylindrical specimens broad applicability, their tests became the few of the best known cases for testing the bond-slip relationship.The proposed joint element method in this paper is used to simulate the behavior of steel bar and concrete by using Shima’s #10 specimen. The three-dimensional FEM model is illustrated in Figure9. Steel bar with 10 D (steel diameter) length is left at free load-side end in order to prevent the concrete splitting failure and make the stress distribution more uniform on the load side. The steel strains in the end of steel bar in the five load stages are 0.786 × 10−3, 1.471 × 10−3, 2.286 × 10−3, 2.886 × 10−3, and 3.586 × 10−3. The problem is numerically simulated with 4700 points and 4003 mesh elements made of 8 node isoparametric elements. Other material parameters are given in Table 3 in detail and the model is represented by a nonlinear bond-slip constitutive relationship, proposed by Shima et al., with its expression shown in (14).Table 3
Material parameters for example 2.
Schedule
Properties
Level
Concrete
Elastic modulus
33 GPa
Poisson’s ratio
0.167
Uniaxial compressive strength
21.6 MPa
Reinforced steel
Elastic modulus
190 GPa
Poisson’s ratio
0.23
Yielding stress
428 MPaFigure 9
Finite element method (FEM) model of Shima #10 specimen.The distribution of steel strain along the steel bar in the first four load steps, obtained by this method, is shown in Figure10. The experiment values obtained by Shima are also plotted for comparison. The calculated results by single spring joint element method match well with the experiment values, especially at low stress levels. With increasing tensile load, the steel strain obtained by the proposed method exhibits a similar distribution discipline but is numerically smaller than the experiment value. This may be attributed to small differences in models between the computational one and the experimental one. On the other hand, at high stress level, the applicability of the bond-slip constitutive relationship can be affected by other factors.Figure 10
Distribution of steel strain along the bar in different load step: (a) first load step, (b) second load step, (c) third load step, and (d) fourth load step.
(a)
First load step
(b)
Second load step
(c)
Third load step
(d)
Fourth load stepFigure10 also shows the calculated results simulated by the double spring joint element model under the same circumstances. Because the specimen in this case is the axial tensile member, there is no problem in the choice of normal stiffness and mutual embedding. These two methods are almost completely consistent, which in return verifies the validity and effectiveness of the proposed method.


### 5.3. Three-Point Bending RC Beam

To demonstrate the advantage of the model proposed in this paper, the numerical simulation of the three-point bending RC beam test conducted in the Stevin Laboratory at Delft University of Technology by Slobbe et al. [35] is presented. The shear critical reinforced concrete beam without shear reinforcement showed an extremely brittle failure behavior. Figure 11 sketches the geometrical dimensions and the type of loading. The 3D model beam has a thickness of 200 mm. The span ratio (a
/
d) is 2.93. The beam is reinforced with a single layer of three longitudinal steel bars, of diameter d
=
20 mm each. The ratio between the area of the steel cross section and the concrete cross section is equal to 1.15%. The fine discretization is required by both the characteristic length of the concrete and the small steel cross section leading to a mesh of 9120 concrete elements and 300 steel elements in FEM calculation. With respect to the adopted solution strategy of the NLFE analyses, the analyses have been run in force control with a small step size. Maximal 50 iterations per load step are used. A force tolerance of 10−5 is used as convergence criterion.Figure 11
Slender concrete beams with steel bars (mm).The material of concrete is described by the damage model of Section4 and the constitutive relation curves for concrete are shown in Figure 12. The mechanical properties of concrete and steel bars are given in Table 4. The adopted constitutive relation curves for concrete and material properties have an excellent agreement with the simulation by Slobbe et al. [35].Table 4
Mechanical properties of concrete and reinforced bar for the three-point bending test.
Schedule
Properties
Level
Concrete
Elastic modulus
33.551 GPa
Poisson’s ratio
0.2
Uniaxial compressive strength
38 MPa
Uniaxial tensile strength
3.5 MPa
Fracture energy
76.5 N/m
Reinforced steel
Elastic modulus
220 GPa
Poisson’s ratio
0.3
Yielding stress
587 MPaFigure 12
The adopted constitutive relation curves for concrete.Figure13 gives a comparison of load-displacement curves obtained in the test with their numerical counterparts by Slobbe et al. and the current simulation using the single spring joint element method. Embedded reinforcement model was used in Slobbe et al.’s simulation. As referred to in Introduction, perfect bond was assumed between the concrete and reinforcement in the embedded reinforcement model. Load-deflection curve showed two clear stages: before crack appeared (I) and after crack appeared (II). When the second stage began, cracks appearance corresponds to the degradation process of the concrete and the elastic behavior of the steel bars. In consequence, the slope of the structural curve in this branch depends strongly on the steel participation. The same trend of load-deflection curve can be observed between the numerical simulation and the experimental test. While the numerical results showed a stronger and slightly stiffer behavior, further, it is observed that the initial stiffness of the beam is overestimated in the numerical analyses. This could be due to an overestimation of the input concrete material parameter E
0. Apart from that, in general, the numerical results revealed the same type of behavior with the experimental curve. More details of the damage behavior of the RC beam models at different load levels are provided in Figure 14, confirming that the crack is at the same position as shown in the experimental observation.Figure 13
Three-point bending test: load-deflection curves between the numerical analysis and experiment.Figure 14
(a) Crack propagation from FE analysis compared with experiment. (b) Damage evolution at different characteristic load levels: 8.4 kN, 10.5 kN, 14.7 kN, and 17.22 kN.
(a)
(b)Figure15 shows load-deflection curves of five different kinds of normal stiffness in double spring joint element method (DSEM) compared with single spring-joint element method, while Figure 16 shows mutual embedment amount by different stiffness parameters with the increasing load levels (take a point on the intermediate joint element of reinforced concrete beams). At the first stage, the load-deflection curves are almost the same for the two methods. This is because the change in the beam’s deflection depends largely on the parameters of concrete in this stage, while the bond effect shows small effect. However, at the first stage, when the normal stiffness is taken to be 102 and 104, the mutual embedment amount is relatively large and shows an unrealistic behavior. At the second stage, as normal stiffness increases, the deflection curve from DSEM gets closer to the curve simulated with the single spring joint element method. When the normal stiffness is taken to be 102 and 104, the mutual embedment amount achieves 1.5 mm. However, the maximum displacement of the beam is only 2.5 mm, which is also not realistic. When normal stiffness is larger than 106, the penetration is inversely proportional to the normal stiffness. When the normal stiffness reaches 108 and 1010, the embedment amount has relatively small effect on the results, which can be ignored. On the other hand, at the same load levels, with the increase of the normal stiffness, the iterative convergence of the finite element calculation becomes worse.Figure 15
Load-deflection curves of five different kinds of normal stiffness in double spring joint element method (DSEM) compared with single spring joint element.Figure 16
Load-deflection curves of five different kinds of normal stiffness in double spring joint element method (DSEM) compared with single spring joint element.Overall, the results demonstrate the ability of the isotropic damage model to describe the bending of reinforcement concrete in multiaxial stress state and the advantages of the proposed model over the typical DSEM. Next, this model will be applied to another four-point bending RC beam to see the development of crack in reinforcement concrete; the results using different bond-slip relationships are also showed.


### 5.4. Four-Point Bending RC Beam

The four-point bending test considers the reinforced concrete beam of Figure17, which was experimentally loaded up until failure by Walraven [38] and also simulated by Scotta et al. [39] who proposed a scalar damage model with a shear retention factor for the analysis of reinforced concrete structures. The beam has no shear reinforcement. The experimental observation showed that failure was imputable to the shear force, after the developing of diagonal cracks which reduces the available area of concrete in the vicinity of the load points. In this test, the single spring joint method and the isotropic damage model with four-parameter equivalent strain are applied to analyze the isotropic damage of reinforced concrete with this model.Figure 17
Slender concrete beams with steel bars (cm).The problem is numerically simulated with an encryption grid with 346 × 20 × 8 concrete elements and 1632 reinforcing bar elements in 3D calculation (Figure17). Three steel bars were set along the concrete. The properties of the concrete and of the reinforcing steel used in the test are summarized in Table 5.Table 5
Mechanical properties for the four-point bending test.
Schedule
Properties
Level
Concrete
Elastic modulus
28 GPa
Poisson’s ratio
0.2
Uniaxial compressive strength
30 MPa
Uniaxial tensile strength
2.5 MPa
Fracture energy
10 N/m
Reinforced steel
Elastic modulus
206 GPa
Poisson’s ratio
0.3
Yielding stress
440 MPaThe resulting numerical load-deflection curve is plotted in Figure18, compared with the experimental one and the numerical one simulated by Scotta et al. [39]. The load-deflection curve showed two stages obviously: elastic stage (I) and stage after cracks appear (II). In the first stage before concrete cracked, the stress of the reinforcing bar and concrete on the cross section is not large, the deformation is essentially elastic, and thus load-deflection remains approximately in linear relationship with increasing load. Tension is shared between the two materials, both of about the same strain. In the second phase, load continues to increase and tensile strain of concrete over the tensile edges up to the ultimate deformation. Cracks appear in the concrete tension zone, and the concrete is out of work. Then, stiffness of the cross section decreased significantly. So growth of the deflection is faster than stage I. After reaching the cracking load (about 23 kN), tension is almost entirely borne by the steel. Hereafter, there is an increase in deflection as in Figure 18. The experimental observation was that vertical cracks formed firstly in the central zone. On subsequent loading, the stresses rotated and new inclined cracks formed in the region between the point load and the supports. At the moment of failure, both the steel and the compressed concrete are far below their limit strength, which indicates that the failure of the beam is without doubt imputable to the shear force. The damage distributions at different load levels plotted in Figure 20 confirm the experimental observation and testify the presence at incipient failure of a diagonal band at about 45°, in the same position as the experimental shear cracks. The damage penetrates under the compression zone, reducing the resistant area, finally leading to the failure.Figure 18
Four-point bending test: comparison between the numerical and experimental load-deflection curves.In the experiment by Walraven, the reported failure load is about 70 kN, while the numerical model predicts a failure load of about 40 kN by Scotta et al. and 60 kN in this paper. This may be attributed to the inability of the numerical model without considering the various mechanisms of shear transfer. As this is not the issue of interest at the current stage of research, detailed review and discussion of this problem can be found in Roberto’s paper. In the rear portion of the second stage, the curve from the current simulation is closer to the experimental values. This may be attributed to the damage model employed. When the normal stiffness is chosen to be 108, the results are almost the same by the single spring element method or double spring element method. However, there is normal embedment amount between concrete and steel bar in the double spring element as example 1 in Section 5.1. The embedment amount of about 6 × 10−8 m in stage (I) and 8 × 10−7 m in stage (II) in the middle of the beam is relatively large than the ends of beam. When larger normal stiffness is chosen, the calculation does not converge. Figure 19 shows three different curves, simulated by bond-slip curves proposed by three different researchers: Nilson [10], Mirza and Houde [40], and Di [41]. It can be seen that the results obtained using the different curves are almost the same, which implies that the impact of reinforced concrete structures subjected to bending deflection curve is relatively small. The three curves are almost the same in stage (I) implying that concrete plays a major role in resisting moment. Additionally, small differences in stage (II) imply reinforcement is also involved in the role.Figure 19
Calculated vertical forceF versus deflection u with various bond definitions.Figure 20
Damage evolution and crack propagation from FE analysis at different characteristic load levels: (a) Scotta et al.’s simulation, at 21, 30, 34, and 38 kN; (b) this paper, at 23, 34, 43, and 59 kN; and (c) three-dimensional damage pattern, at 59 kN.
(a)
(b)
(c)Figure21 illustrates the longitudinal force distribution of the steel bars along axial direction; it corresponds to the different loading stages. Note that before concrete cracked, the force in steel bar increased as load step increased, and the distribution of steel bar element force is smooth; the force is almost the same between two loading points, while the force gradually increased from displacement constraint points to loading points. After cracks appeared, uneven distributions are shown between loading points. This is apparently attributed to the softening property of concrete. And we can see a larger force at loading points.Figure 21
Longitudinal bar element force distribution of steel rebar.


## 5.1. Unbonded Prestressed RC Structure

A reinforced concrete block structure with U-shaped unbonded prestressed steel bar is considered. The radius of bending steel is equal to 2 m, and the bottom of the concrete block is fixed. The diameter of the reinforcing steel is 36 mm with material parameters listed in Table1 and its end is subjected to a concentrated load equal to 1000 kN. For this simulation, bend steel bar is meshed by beam elements. FEM model and its dimensions are illustrated in Figure 7 To highlight the effectiveness of the single spring joint method in describing the normal interaction of the concrete and steel bar; the bond stress between steel bar and concrete is not considered; that is, τ
=
0.Table 1
Material parameters for example 1.
Schedule
Properties
Level
Concrete
Elastic modulus
25.5 GPa
Poisson’s ratio
0.167
Reinforced steel
Elastic modulus
200 GPa
Poisson’s ratio
0.23Figure 7
Finite element method (FEM) model for example 1.Figure8(a) represents the distribution of force along the steel bar. The value is from 499.905 kN/m to 500.015 kN/m, which is close to the analytical solution 500 kN/m obtained from (20) and illustrated in Figure 8(b). Although distribution of force is close between the two methods for this simple test case, there is no normal interpenetration amount between concrete and steel with the single spring joint element method as compared with the typical double spring joint element method. The relationship between normal stiffness K
V and normal embedment amount is given in Table 2, and the order of normal embedment amount is inversely proportional to the order of K
V. When K
V is 1 kPa, the normal embedment amount is 0.1013 × 10−1 m. While K
V is 1 × 107 kPa, the embedment amount is reduced to 0.13 × 10−8 m, and when K
V continues to increase, the convergence deteriorates or even fails:(20)
∫
σ
d
r
=
∫
0
π
σ
r
d
θ
=
2
F
.Table 2
The relationship between normal stiffness and normal embedment amount.
Normal stiffnessK
V (KPa)
Normal embedment amount (m)
1
0.1013
×
10
-
1
1
×
10
2
0.1013
×
10
-
3
1
×
10
4
0.1013
×
10
-
5
1
×
10
7
0.1013
×
10
-
8
1
×
10
12
0.1013
×
10
-
13Figure 8
(a) Distribution of force along the steel bar. (b) The analytical solution of example 1.
(a)
(b)


## 5.2. Shima Tensile Specimen

In order to obtain bond-slip relationship suitable for different thicknesses of the concrete cover, anchor length, and different parameters of materials, Shima et al. [37], among others, did a series of pull tests and tensile tests of cylindrical specimens. With cylindrical specimens broad applicability, their tests became the few of the best known cases for testing the bond-slip relationship.The proposed joint element method in this paper is used to simulate the behavior of steel bar and concrete by using Shima’s #10 specimen. The three-dimensional FEM model is illustrated in Figure9. Steel bar with 10 D (steel diameter) length is left at free load-side end in order to prevent the concrete splitting failure and make the stress distribution more uniform on the load side. The steel strains in the end of steel bar in the five load stages are 0.786 × 10−3, 1.471 × 10−3, 2.286 × 10−3, 2.886 × 10−3, and 3.586 × 10−3. The problem is numerically simulated with 4700 points and 4003 mesh elements made of 8 node isoparametric elements. Other material parameters are given in Table 3 in detail and the model is represented by a nonlinear bond-slip constitutive relationship, proposed by Shima et al., with its expression shown in (14).Table 3
Material parameters for example 2.
Schedule
Properties
Level
Concrete
Elastic modulus
33 GPa
Poisson’s ratio
0.167
Uniaxial compressive strength
21.6 MPa
Reinforced steel
Elastic modulus
190 GPa
Poisson’s ratio
0.23
Yielding stress
428 MPaFigure 9
Finite element method (FEM) model of Shima #10 specimen.The distribution of steel strain along the steel bar in the first four load steps, obtained by this method, is shown in Figure10. The experiment values obtained by Shima are also plotted for comparison. The calculated results by single spring joint element method match well with the experiment values, especially at low stress levels. With increasing tensile load, the steel strain obtained by the proposed method exhibits a similar distribution discipline but is numerically smaller than the experiment value. This may be attributed to small differences in models between the computational one and the experimental one. On the other hand, at high stress level, the applicability of the bond-slip constitutive relationship can be affected by other factors.Figure 10
Distribution of steel strain along the bar in different load step: (a) first load step, (b) second load step, (c) third load step, and (d) fourth load step.
(a)
First load step
(b)
Second load step
(c)
Third load step
(d)
Fourth load stepFigure10 also shows the calculated results simulated by the double spring joint element model under the same circumstances. Because the specimen in this case is the axial tensile member, there is no problem in the choice of normal stiffness and mutual embedding. These two methods are almost completely consistent, which in return verifies the validity and effectiveness of the proposed method.


## 5.3. Three-Point Bending RC Beam

To demonstrate the advantage of the model proposed in this paper, the numerical simulation of the three-point bending RC beam test conducted in the Stevin Laboratory at Delft University of Technology by Slobbe et al. [35] is presented. The shear critical reinforced concrete beam without shear reinforcement showed an extremely brittle failure behavior. Figure 11 sketches the geometrical dimensions and the type of loading. The 3D model beam has a thickness of 200 mm. The span ratio (a
/
d) is 2.93. The beam is reinforced with a single layer of three longitudinal steel bars, of diameter d
=
20 mm each. The ratio between the area of the steel cross section and the concrete cross section is equal to 1.15%. The fine discretization is required by both the characteristic length of the concrete and the small steel cross section leading to a mesh of 9120 concrete elements and 300 steel elements in FEM calculation. With respect to the adopted solution strategy of the NLFE analyses, the analyses have been run in force control with a small step size. Maximal 50 iterations per load step are used. A force tolerance of 10−5 is used as convergence criterion.Figure 11
Slender concrete beams with steel bars (mm).The material of concrete is described by the damage model of Section4 and the constitutive relation curves for concrete are shown in Figure 12. The mechanical properties of concrete and steel bars are given in Table 4. The adopted constitutive relation curves for concrete and material properties have an excellent agreement with the simulation by Slobbe et al. [35].Table 4
Mechanical properties of concrete and reinforced bar for the three-point bending test.
Schedule
Properties
Level
Concrete
Elastic modulus
33.551 GPa
Poisson’s ratio
0.2
Uniaxial compressive strength
38 MPa
Uniaxial tensile strength
3.5 MPa
Fracture energy
76.5 N/m
Reinforced steel
Elastic modulus
220 GPa
Poisson’s ratio
0.3
Yielding stress
587 MPaFigure 12
The adopted constitutive relation curves for concrete.Figure13 gives a comparison of load-displacement curves obtained in the test with their numerical counterparts by Slobbe et al. and the current simulation using the single spring joint element method. Embedded reinforcement model was used in Slobbe et al.’s simulation. As referred to in Introduction, perfect bond was assumed between the concrete and reinforcement in the embedded reinforcement model. Load-deflection curve showed two clear stages: before crack appeared (I) and after crack appeared (II). When the second stage began, cracks appearance corresponds to the degradation process of the concrete and the elastic behavior of the steel bars. In consequence, the slope of the structural curve in this branch depends strongly on the steel participation. The same trend of load-deflection curve can be observed between the numerical simulation and the experimental test. While the numerical results showed a stronger and slightly stiffer behavior, further, it is observed that the initial stiffness of the beam is overestimated in the numerical analyses. This could be due to an overestimation of the input concrete material parameter E
0. Apart from that, in general, the numerical results revealed the same type of behavior with the experimental curve. More details of the damage behavior of the RC beam models at different load levels are provided in Figure 14, confirming that the crack is at the same position as shown in the experimental observation.Figure 13
Three-point bending test: load-deflection curves between the numerical analysis and experiment.Figure 14
(a) Crack propagation from FE analysis compared with experiment. (b) Damage evolution at different characteristic load levels: 8.4 kN, 10.5 kN, 14.7 kN, and 17.22 kN.
(a)
(b)Figure15 shows load-deflection curves of five different kinds of normal stiffness in double spring joint element method (DSEM) compared with single spring-joint element method, while Figure 16 shows mutual embedment amount by different stiffness parameters with the increasing load levels (take a point on the intermediate joint element of reinforced concrete beams). At the first stage, the load-deflection curves are almost the same for the two methods. This is because the change in the beam’s deflection depends largely on the parameters of concrete in this stage, while the bond effect shows small effect. However, at the first stage, when the normal stiffness is taken to be 102 and 104, the mutual embedment amount is relatively large and shows an unrealistic behavior. At the second stage, as normal stiffness increases, the deflection curve from DSEM gets closer to the curve simulated with the single spring joint element method. When the normal stiffness is taken to be 102 and 104, the mutual embedment amount achieves 1.5 mm. However, the maximum displacement of the beam is only 2.5 mm, which is also not realistic. When normal stiffness is larger than 106, the penetration is inversely proportional to the normal stiffness. When the normal stiffness reaches 108 and 1010, the embedment amount has relatively small effect on the results, which can be ignored. On the other hand, at the same load levels, with the increase of the normal stiffness, the iterative convergence of the finite element calculation becomes worse.Figure 15
Load-deflection curves of five different kinds of normal stiffness in double spring joint element method (DSEM) compared with single spring joint element.Figure 16
Load-deflection curves of five different kinds of normal stiffness in double spring joint element method (DSEM) compared with single spring joint element.Overall, the results demonstrate the ability of the isotropic damage model to describe the bending of reinforcement concrete in multiaxial stress state and the advantages of the proposed model over the typical DSEM. Next, this model will be applied to another four-point bending RC beam to see the development of crack in reinforcement concrete; the results using different bond-slip relationships are also showed.


## 5.4. Four-Point Bending RC Beam

The four-point bending test considers the reinforced concrete beam of Figure17, which was experimentally loaded up until failure by Walraven [38] and also simulated by Scotta et al. [39] who proposed a scalar damage model with a shear retention factor for the analysis of reinforced concrete structures. The beam has no shear reinforcement. The experimental observation showed that failure was imputable to the shear force, after the developing of diagonal cracks which reduces the available area of concrete in the vicinity of the load points. In this test, the single spring joint method and the isotropic damage model with four-parameter equivalent strain are applied to analyze the isotropic damage of reinforced concrete with this model.Figure 17
Slender concrete beams with steel bars (cm).The problem is numerically simulated with an encryption grid with 346 × 20 × 8 concrete elements and 1632 reinforcing bar elements in 3D calculation (Figure17). Three steel bars were set along the concrete. The properties of the concrete and of the reinforcing steel used in the test are summarized in Table 5.Table 5
Mechanical properties for the four-point bending test.
Schedule
Properties
Level
Concrete
Elastic modulus
28 GPa
Poisson’s ratio
0.2
Uniaxial compressive strength
30 MPa
Uniaxial tensile strength
2.5 MPa
Fracture energy
10 N/m
Reinforced steel
Elastic modulus
206 GPa
Poisson’s ratio
0.3
Yielding stress
440 MPaThe resulting numerical load-deflection curve is plotted in Figure18, compared with the experimental one and the numerical one simulated by Scotta et al. [39]. The load-deflection curve showed two stages obviously: elastic stage (I) and stage after cracks appear (II). In the first stage before concrete cracked, the stress of the reinforcing bar and concrete on the cross section is not large, the deformation is essentially elastic, and thus load-deflection remains approximately in linear relationship with increasing load. Tension is shared between the two materials, both of about the same strain. In the second phase, load continues to increase and tensile strain of concrete over the tensile edges up to the ultimate deformation. Cracks appear in the concrete tension zone, and the concrete is out of work. Then, stiffness of the cross section decreased significantly. So growth of the deflection is faster than stage I. After reaching the cracking load (about 23 kN), tension is almost entirely borne by the steel. Hereafter, there is an increase in deflection as in Figure 18. The experimental observation was that vertical cracks formed firstly in the central zone. On subsequent loading, the stresses rotated and new inclined cracks formed in the region between the point load and the supports. At the moment of failure, both the steel and the compressed concrete are far below their limit strength, which indicates that the failure of the beam is without doubt imputable to the shear force. The damage distributions at different load levels plotted in Figure 20 confirm the experimental observation and testify the presence at incipient failure of a diagonal band at about 45°, in the same position as the experimental shear cracks. The damage penetrates under the compression zone, reducing the resistant area, finally leading to the failure.Figure 18
Four-point bending test: comparison between the numerical and experimental load-deflection curves.In the experiment by Walraven, the reported failure load is about 70 kN, while the numerical model predicts a failure load of about 40 kN by Scotta et al. and 60 kN in this paper. This may be attributed to the inability of the numerical model without considering the various mechanisms of shear transfer. As this is not the issue of interest at the current stage of research, detailed review and discussion of this problem can be found in Roberto’s paper. In the rear portion of the second stage, the curve from the current simulation is closer to the experimental values. This may be attributed to the damage model employed. When the normal stiffness is chosen to be 108, the results are almost the same by the single spring element method or double spring element method. However, there is normal embedment amount between concrete and steel bar in the double spring element as example 1 in Section 5.1. The embedment amount of about 6 × 10−8 m in stage (I) and 8 × 10−7 m in stage (II) in the middle of the beam is relatively large than the ends of beam. When larger normal stiffness is chosen, the calculation does not converge. Figure 19 shows three different curves, simulated by bond-slip curves proposed by three different researchers: Nilson [10], Mirza and Houde [40], and Di [41]. It can be seen that the results obtained using the different curves are almost the same, which implies that the impact of reinforced concrete structures subjected to bending deflection curve is relatively small. The three curves are almost the same in stage (I) implying that concrete plays a major role in resisting moment. Additionally, small differences in stage (II) imply reinforcement is also involved in the role.Figure 19
Calculated vertical forceF versus deflection u with various bond definitions.Figure 20
Damage evolution and crack propagation from FE analysis at different characteristic load levels: (a) Scotta et al.’s simulation, at 21, 30, 34, and 38 kN; (b) this paper, at 23, 34, 43, and 59 kN; and (c) three-dimensional damage pattern, at 59 kN.
(a)
(b)
(c)Figure21 illustrates the longitudinal force distribution of the steel bars along axial direction; it corresponds to the different loading stages. Note that before concrete cracked, the force in steel bar increased as load step increased, and the distribution of steel bar element force is smooth; the force is almost the same between two loading points, while the force gradually increased from displacement constraint points to loading points. After cracks appeared, uneven distributions are shown between loading points. This is apparently attributed to the softening property of concrete. And we can see a larger force at loading points.Figure 21
Longitudinal bar element force distribution of steel rebar.


## 6. Conclusions

In this paper, a novel single spring joint element method based on mixed coordinate system for bond-slip problems in RC structures has been proposed. Since typical double spring joint element method needs to choose normal stiffness, which may cause mutual embedding problems and bring additional challenges to calculation, the single spring joint element has the advantages in terms of accuracy and computational efficiency. In the proposed method, a pair of nodes are set along the tangent direction of beam elements for steel bar, making DOFs of two materials in the normal direction equal to ensure deformation compatibility instead of choosing normal stiffness arbitrarily. The tangential stiffness of joint element is determined by the relationship of the bond force and slip. Its solid elements for concrete are solved in the global coordinate system, while the beam elements for the steel bar are solved in the local coordinate system to establish the relation between DOFs in the normal direction. In addition, the proposed method can also apply to RC structures with irregular arrangements of steel bars.The proposed method is then applied to an unbonded prestressed reinforcement concrete. In the typical double spring joint element method, there is normal embedment amount between concrete and steel bar. And the normal embedment amount is inversely proportional to the chosen normal stiffness, while, in the proposed method, there is no normal embedment amount and the distribution of force between steel and concrete corresponds qualitatively to analytical solution.In this case of a typical Shima’s #10 specimen, the distribution of steel strain along the steel bar in the first four load steps is shown. The calculated results by two methods coincide well with the experiment values under the same circumstances, especially at low stress levels, verifying the validity and effectiveness of the proposed method.Then, the single spring joint method is applied to analyze the isotropic damage of reinforced concrete with a three-point bending beam, showing consistent midspan deflection and damage distribution of the slender concrete beams in comparison with experiment observations. The results are also compared with double spring joint method choosing different normal stiffness.Finally, the single spring joint method is applied to another four-point bending beam. Four different bond-slip curves proposed by different researchers are employed, and almost identical load-deflection curves are obtained. Damage evolution and crack propagation are observed which are similar to experimental observation and simulations results.



---

*Source: 979678-2015-12-02.xml*