# Influence of Opening Locations on the Structural Response of Shear Wall

Authors: G. Muthukumar; Manoj Kumar
Journal: Journal of Nonlinear Dynamics
 (2014)

---

## Abstract

Shear walls have been conferred as a major lateral load resisting element in a building in any seismic prone zone. It is essential to determine behavior of shear wall in the preelastic and postelastic stage. Shear walls may be provided with openings due to functional requirement of the building. The size and location of opening may play a significant role in the response of shear walls. Though it is a well known fact that size of openings affects the structural response of shear walls significantly, there is no clear consensus on the behavior of shear walls under different opening locations. The present study aims to study the dynamic behavior of shear walls under various opening locations using nonlinear finite element analysis using degenerated shell element with assumed strain approach. Only material nonlinearity has been considered using plasticity approach. A five-parameter Willam-Warnke failure criterion is considered to define the yielding/crushing of the concrete with tensile cutoff. The time history responses have been plotted for all opening cases with and without ductile detailing. The analysis has been done for different damping ratios. It has been observed that the large number of small openings resulted in better displacement response.


---

## Body

## 1. Introduction

Tall reinforced concrete buildings are subjected to lateral loads due to wind and earthquake. In order to resist these lateral loads, shear walls are provided in the framed structure as a lateral load resisting element [1–3]. Shear walls possess sufficient strength and stiffness under any loading conditions. The importance of shear wall in mitigating the damage to reinforced concrete structures is well documented in the literature [3, 4]. Shear walls are generally classified on the basis of aspect ratio (height/width ratio). Shear walls with aspect ratio between 1 and 3 are generally considered to be of squat type and shear walls with aspect ratio greater than 3 are considered to be of slender type. In general, the structural response of shear walls depends strongly on the type of loading, aspect ratio of shear wall, and size and location of the openings in the shear walls. Squat shear walls generally fail in shear mode whereas slender shear walls fail in a flexural mode. The presence of openings in shear walls makes the behavior of shear wall slightly vulnerable under dynamic loading conditions. The structural analysis of the shear walls with openings becomes complex due to the stress concentration near the openings [5]. Various experimental investigations have been performed on shear walls with and without openings subjected to severe dynamic earthquake loading conditions [6]. Neuenhofer [5], in his study on shear wall with opening, has observed that, for the same opening area, the reduction in stiffness for squat and slender shear walls is 50% and 20%, respectively [5]. Thus, the aspect ratio becomes critical for squat shear walls. Few analytical studies have been made on the response of shear wall with openings [5, 7–10]. Rosman [8] developed an approximate linear elastic approach using laminar analysis, based on different assumptions to analyze the shear wall with one row and two rows of openings. Schwaighofer used this approach to analyze the shear wall with three rows of openings and observed that Rosman’s theory predicts the behavior of shear wall with three rows of openings also with sufficient accuracy [9]. Though reasonable studies have been made on the response analysis of shear walls with openings, very few literatures exist on the influence of opening locations on the structural response of shear walls [11, 12]. Hence, it is essential to study the influence of opening locations on the structural response of RC shear walls. The present study analyses the response of RC shear walls for seven different opening locations for both squat and slender shear walls with and without ductile detailing [13] for different damping ratios. Moreover, it was found that the use of conventional methods for the analysis of shear walls with openings resulted in remarkably poor results especially for the squat shear walls where the mode of failure is predominantly shear. It was also shown in literature that the hand calculation either underestimates or overestimates the response in predicting the response of shear walls with openings [5]. The finite element analysis has been the most versatile and successfully employed method of analysis in the past to accurately predict the structural behavior of reinforced concrete shear wall in linear as well as in nonlinear range under any severe loading conditions [14]. With the advent in computing facilities, finite element method has gained an enormous popularity among the structural engineering community, especially in the nonlinear dynamic analysis. The nonlinearity of RC shear wall may be due to geometry or due to material. Since shear wall is a huge structure, the deformation of shear wall has been assumed to be in control and hence the geometric nonlinearity has not been considered. The next section describes the formulation of degenerated shell element formulation.


## 2. Degenerated Shell Element Formulation

The displacement based finite element method has been considered to be the most popular choice because of its simplicity and ease with which the computations can be performed. The use of shell element to model moderately thick structures like shear wall is well documented in the literature [15, 16]. In this section, the underlying basic ideas in the formulation of the degenerated curved shell element are described. Two assumptions are made in the formulation of the curved shell element which is degenerated from three-dimensional solid. First, it is assumed that, even for thick shells, the normal to the middle surface of the element remains straight after deformation. Second, the strain energy corresponding to stresses perpendicular to the middle surface is disregarded; that is, the stress component normal to the shell mid-surface is constrained to be zero. Five degrees of freedom are specified at each nodal point, corresponding to its three translations and two rotations of the normal at each node. The independent definition of the translational and rotational degrees of freedom permits the transverse shear deformation to be taken into account during the formulation of the element stiffness, since rotations are not necessarily normal to the slope of the mid-surface.Coordinate Systems. The geometry of the shell can be represented by the coordinates and normal vectors of its middle surface as Figure 1. The geometry of the degenerated shell element and kinematics of deformation are described by using four different coordinate systems, that is, global, natural, local, and nodal coordinate systems. The global coordinate system (
x
,
y
,
z
) is used to define the shell geometry. The shape functions are expressed in natural curvilinear coordinates ξ
,
η
,
ζ. In order to easily deal with the thin shell assumption of zero normal stress in the z direction, the strain components are defined in terms of local coordinate set of axes x
′
,
y
′
,
z
′. At each node of shell element, the nodal coordinate set (V
¯
1
k
,
V
¯
2
k
,
V
¯
3
k) with unit vectors is defined. The four coordinate sets employed in the present formulation are now described.Figure 1
Geometry of 9-noded degenerated shell element.Global Coordinate Set
(
x
,
y
,
z
)
. This is a Cartesian coordinate system, freely chosen, in relation to which the geometry of the structure is defined in space. Nodal coordinates and displacements, global stiffness matrix, and applied force vector are referred to this system. The displacements corresponding to x, y, and z directions are u, v, and w, respectively.Nodal Coordinate Set
(
V
¯
1
k
,
V
¯
2
k
,
V
¯
3
k
)
. A nodal coordinate system is defined at each nodal point with origin at the reference surface (mid-surface). The vector V
3
k is constructed from nodal coordinates at top and bottom surfaces at node k and is expressed as
(1)
V
3
k
x
V
3
k
y
V
3
k
z
=
x
3
k
top
-
x
3
k
bot
y
3
k
top
-
y
3
k
bot
z
3
k
top
-
z
3
k
bot
.
V
3
k defines the direction of the normal at any node “k” which is not necessarily perpendicular to the mid-surface. The major advantage of the definition of V
3
k with normal not necessary to be perpendicular to mid-surface is that there are no gaps or overlaps along element boundaries.The vectorV
1
k is constructed perpendicular to V
3
k and parallel to the global x
z plane. Hence
(2)
V
1
k
x
=
V
3
k
z
,
V
1
k
y
=
0.0
,
V
1
k
z
=
-
V
3
k
x
.
Alternatively, if the vector V
3
k is in the y direction (V
3
k
x
=
V
3
k
y
=
0.0), the following expressions are assumed: V
1
k
x
=
-
V
3
k
y, V
1
k
y
=
V
1
k
z
=
0.0 (x direction). The superscripts refer to the vector components in the global coordinate system.The vectorV
2
k is constructed perpendicular to the plane defined by V
1
k and V
3
k.HenceV
2
k
=
V
1
k
×
V
3
k. The unit vectors in the directions of V
1
k, V
2
k, and V
3
k are represented by V
¯
1
k, V
¯
2
k, and V
¯
3
k, respectively. The vectors V
¯
1
k, V
¯
2
k define the rotations (β
2
k  and  β
1
k, resp.) of the corresponding normal.Curvilinear Coordinate Set
(
ε
,
η
,
ζ
)
. In this system, ξ, η are the two curvilinear coordinates in the middle plane of the shell element and ζ is a linear coordinate in the thickness direction. It is assumed that ξ
,
η
,
ζ vary between −1 and +
1 on the respective faces of the elements. The relations between curvilinear coordinates and global coordinates are defined later in (4). It should also be noted that ζ direction is only approximately perpendicular to the shell-surface, since ζ is defined as a function of V
¯
3
k.Local Coordinate Set
(
x
′
,
y
′
,
z
′
)
. This is the Cartesian coordinate system defined at the sampling points wherein stresses and strains are to be calculated. The direction z
′ is taken perpendicular to the surface ξ
= constant, being obtained by the cross product of the ξ and η directions. The direction x
′ can be taken tangent to the ξ direction at the sampling point. The direction y
′ is defined by the cross product of the z
′and x
′ directions
(3)
z
′
=
∂
x
∂
ξ
∂
y
∂
ξ
∂
z
∂
ξ
×
∂
x
∂
η
∂
y
∂
η
∂
z
∂
η
;
x
′
=
∂
x
∂
ξ
∂
y
∂
ξ
∂
z
∂
ξ
;
y
′
=
z
′
×
x
′
.Element Geometry. The global coordinates of pairs of points on the top and bottom surface at each node are usually input to define the element geometry. Alternatively, the mid-surface nodal coordinates and the corresponding directional thickness can be furnished. In isoparametric formulation, coordinates of a point within an element are obtained by interpolating the nodal coordinates through element shape functions and are expressed as
(4)
x
y
z
=
∑
k
=
1
9
N
k
ξ
,
η
x
k
mid
y
k
mid
z
k
mid
︸
mid-surface
only
+
∑
k
=
1
9
N
k
ξ
,
η
ζ
h
k
2
V
3
k
x
V
3
k
y
V
3
k
z
︸
effect
of
shell
thickness
,
where x
k
mid, y
k
mid, and z
k
mid are the coordinates of the shell mid-surface and h
k is the shell thickness at node k. In the above expression N
k
ξ
,
η are the element shape functions at the point considered within the element (
ξ
,
η
) and ζ tells the position of the point in the thickness direction. The unit vector in the directions of V
3
k is represented by V
¯
3
k. The element shape functions are calculated in the natural coordinate system as
(5)
N
1
=
1
4
ξ
1
+
ξ
η
1
+
η
;
N
2
=
1
2
1
+
ξ
1
-
ξ
η
1
+
η
;
N
3
=
-
1
4
1
-
ξ
η
1
+
η
;
N
4
=
-
1
2
ξ
1
-
ξ
1
+
η
1
-
η
;
N
5
=
1
4
1
-
ξ
η
1
-
η
;
N
6
=
-
1
2
1
+
ξ
1
-
ξ
η
1
-
η
;
N
7
=
-
1
4
ξ
1
+
ξ
η
1
-
η
;
N
8
=
1
2
ξ
1
+
ξ
1
+
η
1
-
η
;
N
9
=
1
+
ξ
1
-
ξ
1
+
η
1
-
η
.
Based on two assumptions of the degeneration process previously described, the element displacement field can then be expressed by the five degrees of freedom at each node. The global displacements are determined from mid-surface nodal displacements u
k
mid, v
k
mid, and w
k
mid and the relative displacements are caused by the two rotations of the normal as
(6)
u
v
w
=
∑
k
=
1
n
N
k
u
k
mid
v
k
mid
w
k
mid
+
∑
k
=
1
n
N
k
ζ
h
k
2
v
1
k
x
-
v
2
k
x
v
1
k
y
-
v
2
k
y
v
1
k
z
-
v
2
k
z
β
1
k
β
2
k
,
where β
1
k and β
2
k are the rotations of the normals which results in the relative displacements and V
1
k and V
2
k are the unit vectors defined at each node and n is the number of nodes. At any point on the mid-surface of the nodes, an orthogonal set of local coordinates V
¯
1
k, V
¯
2
k, and V
¯
3
k is constructed. V
¯
1
k and V
¯
2
k are constructed in the following manner: V
¯
1
k
=
i
¯
×
V
¯
3
k; V
1
k
=
V
1
i
×
V
3
k. The vectors V
¯
1
k, V
¯
2
k, and V
¯
3
k are mutually perpendicular.


### 2.1. Strain Displacement Relationship

The Mindlin and Reissner type assumptions are used to derive the strain components defined in terms of the local coordinate system of axesx
′
-
y
′
-
z
′ where z
′ is perpendicular to the material surface layer. For the small deformations and neglecting the strain energy associated with stresses perpendicular to the local x
′
-
y
′ surface, the strain components may be written as
(7)
ε
=
ε
f
′
ε
s
′
=
ε
x
′
ε
y
′
γ
x
′
y
′
γ
x
′
z
′
γ
y
′
z
′
=
∂
u
′
∂
x
′
∂
v
′
∂
y
′
∂
u
′
∂
y
′
+
∂
v
′
∂
x
′
∂
u
′
∂
z
′
+
∂
w
′
∂
x
′
∂
v
′
∂
z
′
+
∂
w
′
∂
y
′
.
ε
f
′ and ε
s
′ are the in-plane and transverse shear strains, respectively, and u
′, v
′, and w
′are the displacement components in the local system x
′, y
′, z
′. Assuming the shell mid-surface tangential to x
′
-
y
′, the in-plane shear strains ε
f
′ can further be divided into membrane strains ε
m
′ and bending strains ε
b
′ as
(8)
ε
f
′
=
ε
m
′
+
ε
b
′
,
where
(9)
ε
m
′
=
∂
u
′
∂
x
′
∂
v
′
∂
y
′
∂
u
′
∂
y
′
+
∂
v
′
∂
x
′
,
ε
b
′
=
z
′
∂
θ
x
′
∂
x
′
z
′
∂
θ
y
′
∂
y
′
z
′
∂
θ
x
′
∂
y
′
+
∂
θ
y
′
∂
x
′
.
The transformation matrix [T] is used to convert the strains in local coordinate system into strains in global coordinate system as
(10)
∂
u
′
∂
x
′
∂
v
′
∂
x
′
∂
w
′
∂
x
′
∂
u
′
∂
y
′
∂
v
′
∂
y
′
∂
w
′
∂
y
′
∂
u
′
∂
z
′
∂
v
′
∂
z
′
∂
w
′
∂
z
′
=
T
T
∂
u
∂
x
∂
v
∂
x
∂
w
∂
x
∂
u
∂
y
∂
v
∂
y
∂
w
∂
y
∂
u
∂
z
∂
v
∂
z
∂
w
∂
z
T
,
where transformation matrix [T] is given by
(11)
T
=
∂
x
∂
x
′
∂
x
∂
y
′
∂
x
∂
z
′
∂
y
∂
y
′
∂
y
∂
y
′
∂
y
∂
z
′
∂
z
∂
z
′
∂
z
∂
y
′
∂
z
∂
z
′
.
The derivatives of displacements with respect to Cartesian coordinate system into derivatives of displacements with respect to natural coordinate system are transformed as
(12)
∂
u
∂
x
∂
v
∂
x
∂
w
∂
x
∂
u
∂
y
∂
v
∂
y
∂
w
∂
y
∂
u
∂
z
∂
v
∂
z
∂
w
∂
z
=
J
-
1
∂
u
∂
ξ
∂
v
∂
ξ
∂
w
∂
ξ
∂
u
∂
η
∂
v
∂
η
∂
w
∂
η
∂
u
∂
ζ
∂
v
∂
ζ
∂
w
∂
ζ
,
where [J] is the Jacobian matrix defined as
(13)
J
=
∂
x
∂
ξ
∂
y
∂
ξ
∂
z
∂
ξ
∂
x
∂
η
∂
y
∂
η
∂
z
∂
η
∂
x
∂
∂
y
∂
ς
∂
z
∂
ς
.
The strain displacement matrix [B] relates the strain components and the nodal variables as
(14)
ε
=
B
δ
,
where
(15)
δ
=
u
v
w
β
1
β
2
T
.
The layered element formulation [17] allows the integration through the element thicknesses, which are divided into several concrete and steel layers. Each layer is assumed to have one integration point at its mid-surface. The steel layers are used to model the in-plane reinforcement only. The strain displacement matrixB and the material stiffness matrixD are evaluated at the midpoint of each layer and for all integration points in the plane of the layer. The element stiffness matrix K
e is defined using numerical integration as follows:
(16)
K
e
=
∭
B
T
D
B
d
V
,
where the integration is made over the volume of the element.In the Euclidean space, the volume element is given by the product of the differentials of the Cartesian coordinates and is expressed asd
V
=
d
x
d
y
d
z. Using numerical integration, the volume integration is converted into area integration using Jacobian and is expressed as
(17)
K
e
=
∬
B
T
D
B
J
d
ζ
d
A
.
Similarly, the internal force vector is expressed as {
f
e
} as
(18)
f
e
=
∬
B
T
σ
J
d
ζ
d
A
.
The element stiffness matrix relates the force vector with the displacement vector as
(19)
f
=
K
δ
,
where
(20)
∫
d
A

=

J
∬
-
1
+
1

d
ξ
d
η
Integration
on
⁡
layer
mid-surface
.
Once the displacements are determined, the strains and stresses are calculated using strain displacement matrix and material constitutive matrix, respectively. The formulation of degenerated shell element is completely described in Huang [18].


### 2.2. Assumed Strain Approach

Nevertheless, the general shell theory based on the classical approach has been found to be complex in the finite element formulation. On the other hand, the degenerated shell element [19, 20] derived from the three-dimensional element has been quite successful in modeling moderately thick structures because of their simplicity and circumvents the use of classical shell theory. The degenerated shell element is based on assumption that the normal to the mid-surface remain straight but not necessarily normal to the mid-surface after deformation. Also, the stresses normal to the mid-surface are considered to be negligible. However, when the thickness of element reduces, degenerated shell element has suffered from shear locking and membrane locking when subjected to full numerical integration. The shear locking and membrane locking are the parasitic shear stresses and membrane stresses present in the finite element solution. In order to alleviate locking problems, the reduced integration technique has been suggested and adopted by many authors [21, 22]. However, the use of reduced integration resulted in spurious mechanisms or zero energy modes in some cases. The reduced integration ignores the high ranked terms in interpolated shear strain by numerical integration, thus introducing the chance of development of spurious or zero energy modes in the element. The selective integration, wherein different integration orders are used to integrate the bending, shear, and membrane terms of stiffness matrix, avoids the locking in most of the cases.The assumed strain approach has been successfully adopted by many researchers [23, 24] as an alternative to avoid locking. In the assumed strain based degenerated shell elements, the transverse shear strain and membrane strains are interpolated from the assumed sampling points obtained from the compatibility requirement between flexural and shear strain fields, respectively. The assumed transverse shear strain fields, interpolated at the six appropriately located sampling points, as shown in Figure 2 are
(21)
γ
-
ξ
ζ
=
∑
i
=
1
3
∑
j
=
1
2
P
i
η
·
Q
j
ξ
γ
ξ
ζ
i
j
,
γ
-
η
ζ
=
∑
i
=
1
3
∑
j
=
1
2
P
i
ξ
·
Q
j
η
γ
η
ς
i
j
.
γ
ξ
ς
i
j and γ
η
ζ
i
j are the shear strains obtained from Lagrangian shape functions. The interpolating functions P
i
z and Q
j
z are
(22)
P
1
z
=
z
2
z
+
1
,
P
2
z
=
1
-
z
2
,
P
3
z
=
z
2
z
-
1
,
Q
1
z
=
1
2
1
+
3
z
,
Q
2
z
=
1
2
1
-
3
z
.
Hence, it can be observed that γ
-
ξ
ζ is linear in ξ direction and quadratic in η direction, while γ
-
η
ζ is linear in η direction and quadratic in ξ direction. The polynomial terms for curvature of nine-node Lagrangian elements, κ
ξ and κ
η, are the same as the assumed shear strain, as given by
(23)
κ
ξ
=
∂
θ
ξ
∂
ξ
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
,
η
2
,
ξ
η
2
,
ξ
2
η
2
,
κ
ξ
=
κ
ξ
1
,
ξ
,
η
,
ξ
η
,
η
2
,
ξ
η
2
,
κ
η
=
∂
θ
η
∂
η
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
,
η
2
,
ξ
η
2
,
ξ
2
η
2
,
κ
η
=
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
2
,
γ
-
ξ
ζ
=
γ
-
ξ
ζ
κ
ξ
=
γ
-
ξ
ζ
1
,
ξ
,
η
,
ξ
η
,
η
2
,
ξ
η
2
,
γ
-
η
ζ
=
γ
-
η
ζ
κ
η
=
γ
-
η
ζ
1
,
ξ
,
η
,
ξ
η
,
η
2
,
ξ
η
2
.
The original shear strains obtained from the Lagrange shape functions γ
ξ
ζ and γ
η
ζ are
(24)
γ
ξ
ζ
=
θ
ξ
+
∂
w
∂
ξ
=
γ
ξ
ζ
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
,
η
2
,
ξ
η
2
,
ξ
2
η
2
,
γ
η
ζ
=
θ
η
+
∂
w
∂
η
=
γ
η
ζ
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
,
η
2
,
ξ
η
2
,
ξ
2
η
2
.
The total potential energy expression has the form
(25)
π
-
=
π
+
∫
λ
13
γ
-
ξ
ζ
-
γ
ξ
ζ
d
V
+
∫
λ
23
γ
-
η
ζ
-
γ
η
ζ
d
V
.
λ
13 and λ
23 are Lagrangian multipliers and are independent functions. The terms γ
ξ
ζ and γ
η
ζ are the transverse shear strains evaluated from the displacement field.Sampling point locations for assumed shear/membrane strains.
(a)
(b)The assumed shear strain fields are chosen as(26)
γ
-
ξ
ζ
=
∑
i
=
1
n
R
i
ξ
,
η
γ
ξ
ζ
i
;
γ
-
η
ζ
=
∑
i
=
1
n
S
i
ξ
,
η
γ
η
ζ
i
.
The Lagrangian multipliers are taken as
(27)
λ
13
=
∑
i
=
1
n
λ
i
13
δ
ξ
~
i
-
ξ
δ
η
~
i
-
η
,
λ
23
=
∑
i
=
1
n
λ
i
23
δ
ξ
-
i
-
ξ
δ
η
-
i
-
η
.
By substituting, the following equation can be obtained:
(28)
γ
-
ξ
ζ
ξ
~
i
,
η
~
i
=
γ
ξ
ζ
ξ
~
i
,
η
~
i
,
γ
-
η
ζ
ξ
-
i
,
η
-
i
=
γ
η
ζ
ξ
-
i
,
η
-
i
.
θ
ξ and θ
η are the rotating normal and w is the transverse displacement of the element. It can be clearly seen that the original shear strain and assumed shear strain are not compatible and hence the shear locking exists for very thin shell cases. The appropriately chosen polynomial terms and sampling points ensure the elimination of risk of spurious zero energy modes. The assumed strain can be considered a special case of integration scheme, wherein for function γ
-
ξ
ζ full integration is employed in η direction and reduced integration is employed in ξ direction. On the other hand, for function γ
-
η
ζ reduced integration is employed in η direction and full integration is employed in ξ direction. The membrane and shear strains are interpolated from identical sampling points even though the membrane strains are expressed in orthogonal curvilinear coordinate system and transverse shear strains are expressed in natural coordinate system. The flow chart explaining the formulation of strain displacement matrix using assumed strain approach is shown in Figure 3.Figure 3
Determination of strain displacement matrix using the assumed strain approach.


## 2.1. Strain Displacement Relationship

The Mindlin and Reissner type assumptions are used to derive the strain components defined in terms of the local coordinate system of axesx
′
-
y
′
-
z
′ where z
′ is perpendicular to the material surface layer. For the small deformations and neglecting the strain energy associated with stresses perpendicular to the local x
′
-
y
′ surface, the strain components may be written as
(7)
ε
=
ε
f
′
ε
s
′
=
ε
x
′
ε
y
′
γ
x
′
y
′
γ
x
′
z
′
γ
y
′
z
′
=
∂
u
′
∂
x
′
∂
v
′
∂
y
′
∂
u
′
∂
y
′
+
∂
v
′
∂
x
′
∂
u
′
∂
z
′
+
∂
w
′
∂
x
′
∂
v
′
∂
z
′
+
∂
w
′
∂
y
′
.
ε
f
′ and ε
s
′ are the in-plane and transverse shear strains, respectively, and u
′, v
′, and w
′are the displacement components in the local system x
′, y
′, z
′. Assuming the shell mid-surface tangential to x
′
-
y
′, the in-plane shear strains ε
f
′ can further be divided into membrane strains ε
m
′ and bending strains ε
b
′ as
(8)
ε
f
′
=
ε
m
′
+
ε
b
′
,
where
(9)
ε
m
′
=
∂
u
′
∂
x
′
∂
v
′
∂
y
′
∂
u
′
∂
y
′
+
∂
v
′
∂
x
′
,
ε
b
′
=
z
′
∂
θ
x
′
∂
x
′
z
′
∂
θ
y
′
∂
y
′
z
′
∂
θ
x
′
∂
y
′
+
∂
θ
y
′
∂
x
′
.
The transformation matrix [T] is used to convert the strains in local coordinate system into strains in global coordinate system as
(10)
∂
u
′
∂
x
′
∂
v
′
∂
x
′
∂
w
′
∂
x
′
∂
u
′
∂
y
′
∂
v
′
∂
y
′
∂
w
′
∂
y
′
∂
u
′
∂
z
′
∂
v
′
∂
z
′
∂
w
′
∂
z
′
=
T
T
∂
u
∂
x
∂
v
∂
x
∂
w
∂
x
∂
u
∂
y
∂
v
∂
y
∂
w
∂
y
∂
u
∂
z
∂
v
∂
z
∂
w
∂
z
T
,
where transformation matrix [T] is given by
(11)
T
=
∂
x
∂
x
′
∂
x
∂
y
′
∂
x
∂
z
′
∂
y
∂
y
′
∂
y
∂
y
′
∂
y
∂
z
′
∂
z
∂
z
′
∂
z
∂
y
′
∂
z
∂
z
′
.
The derivatives of displacements with respect to Cartesian coordinate system into derivatives of displacements with respect to natural coordinate system are transformed as
(12)
∂
u
∂
x
∂
v
∂
x
∂
w
∂
x
∂
u
∂
y
∂
v
∂
y
∂
w
∂
y
∂
u
∂
z
∂
v
∂
z
∂
w
∂
z
=
J
-
1
∂
u
∂
ξ
∂
v
∂
ξ
∂
w
∂
ξ
∂
u
∂
η
∂
v
∂
η
∂
w
∂
η
∂
u
∂
ζ
∂
v
∂
ζ
∂
w
∂
ζ
,
where [J] is the Jacobian matrix defined as
(13)
J
=
∂
x
∂
ξ
∂
y
∂
ξ
∂
z
∂
ξ
∂
x
∂
η
∂
y
∂
η
∂
z
∂
η
∂
x
∂
∂
y
∂
ς
∂
z
∂
ς
.
The strain displacement matrix [B] relates the strain components and the nodal variables as
(14)
ε
=
B
δ
,
where
(15)
δ
=
u
v
w
β
1
β
2
T
.
The layered element formulation [17] allows the integration through the element thicknesses, which are divided into several concrete and steel layers. Each layer is assumed to have one integration point at its mid-surface. The steel layers are used to model the in-plane reinforcement only. The strain displacement matrixB and the material stiffness matrixD are evaluated at the midpoint of each layer and for all integration points in the plane of the layer. The element stiffness matrix K
e is defined using numerical integration as follows:
(16)
K
e
=
∭
B
T
D
B
d
V
,
where the integration is made over the volume of the element.In the Euclidean space, the volume element is given by the product of the differentials of the Cartesian coordinates and is expressed asd
V
=
d
x
d
y
d
z. Using numerical integration, the volume integration is converted into area integration using Jacobian and is expressed as
(17)
K
e
=
∬
B
T
D
B
J
d
ζ
d
A
.
Similarly, the internal force vector is expressed as {
f
e
} as
(18)
f
e
=
∬
B
T
σ
J
d
ζ
d
A
.
The element stiffness matrix relates the force vector with the displacement vector as
(19)
f
=
K
δ
,
where
(20)
∫
d
A

=

J
∬
-
1
+
1

d
ξ
d
η
Integration
on
⁡
layer
mid-surface
.
Once the displacements are determined, the strains and stresses are calculated using strain displacement matrix and material constitutive matrix, respectively. The formulation of degenerated shell element is completely described in Huang [18].


## 2.2. Assumed Strain Approach

Nevertheless, the general shell theory based on the classical approach has been found to be complex in the finite element formulation. On the other hand, the degenerated shell element [19, 20] derived from the three-dimensional element has been quite successful in modeling moderately thick structures because of their simplicity and circumvents the use of classical shell theory. The degenerated shell element is based on assumption that the normal to the mid-surface remain straight but not necessarily normal to the mid-surface after deformation. Also, the stresses normal to the mid-surface are considered to be negligible. However, when the thickness of element reduces, degenerated shell element has suffered from shear locking and membrane locking when subjected to full numerical integration. The shear locking and membrane locking are the parasitic shear stresses and membrane stresses present in the finite element solution. In order to alleviate locking problems, the reduced integration technique has been suggested and adopted by many authors [21, 22]. However, the use of reduced integration resulted in spurious mechanisms or zero energy modes in some cases. The reduced integration ignores the high ranked terms in interpolated shear strain by numerical integration, thus introducing the chance of development of spurious or zero energy modes in the element. The selective integration, wherein different integration orders are used to integrate the bending, shear, and membrane terms of stiffness matrix, avoids the locking in most of the cases.The assumed strain approach has been successfully adopted by many researchers [23, 24] as an alternative to avoid locking. In the assumed strain based degenerated shell elements, the transverse shear strain and membrane strains are interpolated from the assumed sampling points obtained from the compatibility requirement between flexural and shear strain fields, respectively. The assumed transverse shear strain fields, interpolated at the six appropriately located sampling points, as shown in Figure 2 are
(21)
γ
-
ξ
ζ
=
∑
i
=
1
3
∑
j
=
1
2
P
i
η
·
Q
j
ξ
γ
ξ
ζ
i
j
,
γ
-
η
ζ
=
∑
i
=
1
3
∑
j
=
1
2
P
i
ξ
·
Q
j
η
γ
η
ς
i
j
.
γ
ξ
ς
i
j and γ
η
ζ
i
j are the shear strains obtained from Lagrangian shape functions. The interpolating functions P
i
z and Q
j
z are
(22)
P
1
z
=
z
2
z
+
1
,
P
2
z
=
1
-
z
2
,
P
3
z
=
z
2
z
-
1
,
Q
1
z
=
1
2
1
+
3
z
,
Q
2
z
=
1
2
1
-
3
z
.
Hence, it can be observed that γ
-
ξ
ζ is linear in ξ direction and quadratic in η direction, while γ
-
η
ζ is linear in η direction and quadratic in ξ direction. The polynomial terms for curvature of nine-node Lagrangian elements, κ
ξ and κ
η, are the same as the assumed shear strain, as given by
(23)
κ
ξ
=
∂
θ
ξ
∂
ξ
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
,
η
2
,
ξ
η
2
,
ξ
2
η
2
,
κ
ξ
=
κ
ξ
1
,
ξ
,
η
,
ξ
η
,
η
2
,
ξ
η
2
,
κ
η
=
∂
θ
η
∂
η
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
,
η
2
,
ξ
η
2
,
ξ
2
η
2
,
κ
η
=
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
2
,
γ
-
ξ
ζ
=
γ
-
ξ
ζ
κ
ξ
=
γ
-
ξ
ζ
1
,
ξ
,
η
,
ξ
η
,
η
2
,
ξ
η
2
,
γ
-
η
ζ
=
γ
-
η
ζ
κ
η
=
γ
-
η
ζ
1
,
ξ
,
η
,
ξ
η
,
η
2
,
ξ
η
2
.
The original shear strains obtained from the Lagrange shape functions γ
ξ
ζ and γ
η
ζ are
(24)
γ
ξ
ζ
=
θ
ξ
+
∂
w
∂
ξ
=
γ
ξ
ζ
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
,
η
2
,
ξ
η
2
,
ξ
2
η
2
,
γ
η
ζ
=
θ
η
+
∂
w
∂
η
=
γ
η
ζ
1
,
ξ
,
η
,
ξ
η
,
ξ
2
,
ξ
2
η
,
η
2
,
ξ
η
2
,
ξ
2
η
2
.
The total potential energy expression has the form
(25)
π
-
=
π
+
∫
λ
13
γ
-
ξ
ζ
-
γ
ξ
ζ
d
V
+
∫
λ
23
γ
-
η
ζ
-
γ
η
ζ
d
V
.
λ
13 and λ
23 are Lagrangian multipliers and are independent functions. The terms γ
ξ
ζ and γ
η
ζ are the transverse shear strains evaluated from the displacement field.Sampling point locations for assumed shear/membrane strains.
(a)
(b)The assumed shear strain fields are chosen as(26)
γ
-
ξ
ζ
=
∑
i
=
1
n
R
i
ξ
,
η
γ
ξ
ζ
i
;
γ
-
η
ζ
=
∑
i
=
1
n
S
i
ξ
,
η
γ
η
ζ
i
.
The Lagrangian multipliers are taken as
(27)
λ
13
=
∑
i
=
1
n
λ
i
13
δ
ξ
~
i
-
ξ
δ
η
~
i
-
η
,
λ
23
=
∑
i
=
1
n
λ
i
23
δ
ξ
-
i
-
ξ
δ
η
-
i
-
η
.
By substituting, the following equation can be obtained:
(28)
γ
-
ξ
ζ
ξ
~
i
,
η
~
i
=
γ
ξ
ζ
ξ
~
i
,
η
~
i
,
γ
-
η
ζ
ξ
-
i
,
η
-
i
=
γ
η
ζ
ξ
-
i
,
η
-
i
.
θ
ξ and θ
η are the rotating normal and w is the transverse displacement of the element. It can be clearly seen that the original shear strain and assumed shear strain are not compatible and hence the shear locking exists for very thin shell cases. The appropriately chosen polynomial terms and sampling points ensure the elimination of risk of spurious zero energy modes. The assumed strain can be considered a special case of integration scheme, wherein for function γ
-
ξ
ζ full integration is employed in η direction and reduced integration is employed in ξ direction. On the other hand, for function γ
-
η
ζ reduced integration is employed in η direction and full integration is employed in ξ direction. The membrane and shear strains are interpolated from identical sampling points even though the membrane strains are expressed in orthogonal curvilinear coordinate system and transverse shear strains are expressed in natural coordinate system. The flow chart explaining the formulation of strain displacement matrix using assumed strain approach is shown in Figure 3.Figure 3
Determination of strain displacement matrix using the assumed strain approach.


## 3. Material Modeling

The modeling of material may play a crucial role in achieving the correct response. The presence of nonlinearity may add another dimension of complexity to it. The nonlinearities in the structure may accurately be estimated and incorporated in the solution algorithm. The accuracy of the solution algorithm depends strongly on the prediction of second-order effects that cause nonlinearities, such as tension stiffening, compression softening, and stress transfer nonlinearities around cracks.These nonlinearities are usually incorporated in the constitutive modeling of the reinforced concrete. In order to incorporate geometric nonlinearity, the second-order terms of strains are to be included. In this study, only material nonlinearity has been considered. The subsequent sections describe the modeling of concrete in compression and tension, modeling of steel.


### 3.1. Concrete Modeling in Tension

The presence of crack in concrete has much influence on the response of nonlinear behavior of reinforced concrete structures. The crack in the concrete is assumed to occur when the tensile stress exceeds the tensile strength. The cracking of concrete results in the loss of continuity in the load transfer and hence the stresses in both concrete and steel reinforcement differ significantly. Hence the analysis of concrete fracture has been very important in order to predict the response of structure precisely. The numerical simulation of concrete fracture can be represented either by discrete crack, proposed by Ngo and Scordelis [25], or by smeared crack, proposed by Rashid [26]. The objective of discrete crack is to simulate the initiation and propagation of dominant cracks present in the structure. In the case of discrete crack approach, nodes are disassociated due to the presence of cracks and therefore the structure requires frequent renumbering of nodes, which may render the huge computational cost. Nevertheless, when the structure’s behavior has been dominated by only few dominant cracks, the discrete modeling of cracking seems the only choice. On the other hand, the smeared crack approach smears out the cracks over the continuum and captures the deterioration process through the constitutive relationship and reduces the computational cost and time drastically.Crack modeling has gone through several stages due to the advancement in technology and computing facilities. Earlier research work indicates that the formation of crack results in the complete reduction in stresses in the perpendicular direction, thus neglecting the phenomenon called tension stiffening. With the rapid increase in extensive experimental investigations as well as in computing facilities, many finite element codes have been developed for the nonlinear finite element analysis, which incorporates the tension stiffening effect. The first tension stiffening model using degraded concrete modulus was proposed by Scanlon and Murray [27] and subsequently many analytical models have been developed such as Lin and Scordelis [28] model, Vebo and Ghali model, Gilbert and Warner model [29], and Nayal and Rasheed model [30]. The cracks are always assumed to be formed in the direction perpendicular to the direction of the maximum principal stress. These directions may not necessarily remain the same throughout the analysis and loading and hence the modeling of orientation of crack plays a significant role in the response of structure. Still, due to simplicity, many investigations have been performed using fixed crack approach, wherein the direction of principal strain axes may remain fixed throughout the analysis. In this study also, the direction of crack has been considered to be fixed throughout the duration of the analysis. However, the modeling of aggregate interlock has not been taken very seriously. The constant shear retention factor or the simple function has been employed to model the shear transfer across the cracks. Apart from the initiation of crack, the propagation of crack also plays a crucial role in the response of structure. The prediction of crack propagation is a very difficult phenomenon due to scarcity and confliction of test results. Nevertheless, the propagation of cracks plays a crucial role in the response of nonlinear analysis of RC structures. The plain concrete exhibits softening behavior and reinforced concrete exhibits stiffening behavior due to the presence of active reinforcing steel. A gradual release of the concrete stress is adopted in this present study as shown in Figure 4 [31]. The reduction in the stress is given by the following expression:
(29)
E
i
=
α
f
t
′
1
-
ε
i
ε
m
1
ε
i
;
ε
t
≤
ε
i
≤
ε
m
.
α and ε
m are the tension stiffening parameters. ε
m is the maximum value reached by the tensile strain at the point considered. ε
i is the current tensile strain in material direction i. The coefficient depends on the percentage of steel in the section. In the present study, values of α and ε
m are taken as 0.5 and 0.0020, respectively. It has also been reported that the influence of the tension stiffening constants on the response of the structures is generally small and hence the constant value is justified in the analysis [31]. Generally, the cracked concrete can transfer shear forces through dowel action and aggregate interlock. The magnitude of shear moduli has been considerably affected because of extensive cracking in different directions.Figure 4
Tension stiffening effect of cracked concrete.Thus, the reduced shear moduli can be put to incorporate the aggregate interlock and dowel action. In the plain concrete, aggregate interlock is the major shear transfer mechanism and for reinforced concrete, dowel action is the major shear transfer mechanism, with reinforcement ratio being the critical variable. In order to incorporate the aggregate interlock and dowel action, the appropriate value of cracked shear modulus [32] has been considered in the material modeling of concrete.Cracked in One Direction. The stress-strain relationship for cracked concrete where cracking is assumed to take place in only one direction is given as
(30)
σ
1
σ
2
τ
12
τ
13
τ
23
=
0
0
0
0
0
0
E
0
0
0
0
0
G
12
c
0
0
0
0
0
G
13
c
0
0
0
0
0
G
23
ε
1
ε
2
γ
12
γ
13
γ
23
,
G
12
c
=
0.25
×
G
1
-
1
0.004
if
ε
1
<
0.004
,
G
12
c
=
0
if
ε
1
≥
0.004
,
G
13
c
=
G
12
c
;
G
23
=
5
G
6
.Cracked in Two Directions. The stress-strain relationship for cracked concrete where cracking is assumed to take place in both directions is given as
(31)
σ
1
σ
2
τ
12
τ
13
τ
23
=
0
0
0
0
0
0
0
0
0
0
0
0
G
12
c
2
0
0
0
0
0
G
13
c
0
0
0
0
0
G
23
c
ε
1
ε
2
γ
12
γ
13
γ
23
,
G
13
c
=
0.25
×
G
1
-
ε
1
0.004
if
ε
1
<
0.004
,
G
13
c
=
0
if
ε
1
≥
0.004
,
G
23
c
=
0.25
×
G
1
-
ε
2
0.004
,
G
12
c
=
0.5
×
G
13
c
if
G
23
c
<
G
13
c
.
It has also been mentioned by Hinton and Owen [33] that the tensile strength of concrete is a relatively small and unreliable quantity which is not highly influential to the response of structures. In the above stress-strain relationship, the cracked shear modulus (G
c) is assumed to be a function of the current tensile strain. G is the uncracked concrete shear modulus. If the crack closes, the uncracked shear modulus G is assumed in the corresponding direction. Even after the formation of initial cracks, the structure can often deform further without further collapse. In addition to the formation of new cracks, there may be a possibility of crack closing and opening of the existing cracks. If the normal strain across the existing crack becomes greater than that just prior to crack formation, the crack is said to have opened again; otherwise it is assumed to be closed. Nevertheless, if all cracks are closed, then the material is assumed to have gained the status equivalent to that of noncracked concrete with linear elastic behavior.


### 3.2. Concrete Modeling in Compression

The theory of plasticity has been used in the compression modeling of the concrete. The failure surface or bounding surface has been defined to demarcate plastic behavior from the elastic behavior. Failure surface is the important component in the concrete plasticity. Sometimes, the failure surface can be referred to as yield surface or loading surface. The material behaves in the elastic fashion as long as the stress lies below the failure surface. Several failure models have been developed and reported in the literature [34]. Nevertheless, the five-parameter failure model proposed by Willam and Warnke [35] seems to possess all inherent properties of the failure surface. The failure surface is constructed using two meridians, namely, compression meridian and tension meridian. The two meridians are pictorially depicted in a meridian plane and cross section of the failure surface is represented in the deviatoric plane.The variations of the average shear stressesτ
m
t and τ
m
c along tensile (θ
=
0
°) and compressive (θ
=
6
0
°) meridians, as shown in Figure 5, are approximated by second-order parabolic expressions in terms of the average normal stresses σ
m, as follows:
(32)
τ
m
t
f
c
′
=
ρ
t
5
f
c
′
=
a
0
+
a
1
σ
m
f
c
′
+
a
2
σ
m
f
c
′
2
θ
=
0
°
,
τ
m
c
f
c
′
=
ρ
c
5
f
c
′
=
b
0
+
b
1
σ
m
f
c
′
+
b
2
σ
m
f
c
′
2
θ
=
6
0
°
.Willam-Warnke failure model.
(a)
Deviatoric plane in octahedral stress
(b)
Meridian plane in principal stressThese two meridians must intersect the hydrostatic axis at the same pointσ
m
/
f
c
′
=
ξ
-
0 (corresponding to hydrostatic tension); the number of parameters that need to be determined is reduced to five. The five parameters a
0
or
b
0
,
a
1
,
a
2
,
b
1
,
b
2 are to be determined from a set of experimental data, with which the failure surface can be constructed using second-order parabolic expressions. The failure surface is expressed as
(33)
f
σ
m
,
τ
m
,
θ
=
5
τ
m
ρ
σ
m
,
θ
-
1
=
0
,
ρ
θ
=
4
ρ
c
2
-
ρ
t
2
cos
⁡
2
θ
+
5
ρ
t
2
-
4
ρ
t
ρ
c
1
/
2
2
ρ
c
ρ
c
2
-
ρ
t
2
cos
⁡
θ
+
ρ
c
2
ρ
t
-
ρ
c
ρ
c
2
1
/
2
g
g
g
×
4
ρ
c
2
-
ρ
t
2
cos
⁡
2
θ
+
5
ρ
t
2
-
4
ρ
t
ρ
c
1
/
2
×
4
ρ
c
2
-
ρ
t
2
cos
⁡
2
θ
+
ρ
c
-
2
ρ
t
2
-
1
.
The formulation of Willam-Warnke five-parameter material model is described in Chen [34]. Once the yield surface is reached, any further increase in the loading results in the plastic flow. The magnitude and direction of the plastic strain increment are defined using flow rule, which is described in the next section.


#### 3.2.1. Flow Rule

In this method, associated flow rule is employed because of the lack of experimental evidence in nonassociated flow rule. The plastic strain increment expressed in terms of current stress increment is given as(34)
d
ε
i
j
p
=
d
λ
∂
f
σ
∂
σ
i
j
.
d
λ determines the magnitude of the plastic strain increment. The gradient ∂
f
σ
/
∂
σ
i
j defines the direction of plastic strain increment to be perpendicular to the yield surface; f
σ is the loading condition or the loading surfaces.


#### 3.2.2. Hardening Rule

The relationship between loading surfaces (or effective stress) and the plastic work (accumulated plastic strain) is represented by a hardening rule. The “Madrid parabola” is used to define the hardening rule. The isotropic hardening is adopted in the present study, as shown in Figure6
(35)
σ
=
E
0
ε
-
1
2
E
0
ε
0
ε
2
.
E
0 is the initial elasticity modulus, ε is the total strain, and ε
0 is the total strain at peak stress f
c
′. The total strain can be divided into elastic and plastic components as
(36)
ε
=
ε
e
+
ε
p
,
σ
˙
=
D
e
ε
˙
-
ε
˙
p
,
ε
˙
p
=
λ
˙
a
;
a
=
∂
F
∂
σ
.
F is the yield function and λ
˙ is the consistency parameter, which defines the magnitude of the plastic flow.Figure 6
Isotropic hardening with expanding yield surfaces and the corresponding uniaxial stress-strain curve.The loading, unloading conditions (Kuhn-Tucker conditions) can be stated as(37)
λ
˙
≥
0
;
F
≤
0
;
λ
˙
F
=
0
.
The first of these Khun-Tucker conditions indicates that the consistency parameter is nonnegative; the second condition implies that the stress states must lie on or within the yield surface. The third condition ensures that the stresses lie on the yield surface during the plastic loading
(38)
a
T
σ
˙
=
0
;
a
T
D
e
ε
˙
-
ε
˙
p
=
0
;
a
T
D
e
ε
˙
-
λ
˙
a
=
0
,
λ
˙
=
a
T
D
e
ε
˙
a
T
D
e
a
.
The elastoplastic constitutive matrix is given by the following expression:
(39)
D
e
p
=
D
-
D
a
a
T
D
H
+
a
T
D
a
a
=
flow vector, defined by the stress gradient of the yield function; D = constitutive matrix in elastic range. The second term in (39) represents the effect of degradation of material during the plastic loading.


### 3.3. Modeling of Reinforcement in Tension and Compression

Reinforcing bars in structural concrete are generally assumed one-dimensional elements without transverse shear stiffness or flexural rigidity. The reinforcing bar can generally be treated as either discrete or smeared. The major advantage of discrete representation of reinforcing bar is existence of one-to-one correspondence between the real structure and model. In the smeared reinforcement, the average stress-strain relationship is calculated for an element area and incorporated directly as part of the overall concrete element stiffness matrix. In the present investigation, the smeared layered approach is adopted. The bilinear stress-strain curve with linear elastic and strain hardening region is adopted in this study as shown in Figure7. Sometimes, the trilinear idealization has also been adopted as the stress-strain curve in tension. Typically, the hardening strain modulus is assumed to be 1% of initial elasticity modulus. The position and thickness of steel layers are to be defined as input parameters along with the elasticity modulus and hardening modulus. The direction of steel (horizontal or vertical) can be set up by defining the angle with respect to local x-axis. There can only be two states of stress for the reinforcing bar, namely, elastic and linear strain hardening.Figure 7
Stress-strain curve of steel.


## 3.1. Concrete Modeling in Tension

The presence of crack in concrete has much influence on the response of nonlinear behavior of reinforced concrete structures. The crack in the concrete is assumed to occur when the tensile stress exceeds the tensile strength. The cracking of concrete results in the loss of continuity in the load transfer and hence the stresses in both concrete and steel reinforcement differ significantly. Hence the analysis of concrete fracture has been very important in order to predict the response of structure precisely. The numerical simulation of concrete fracture can be represented either by discrete crack, proposed by Ngo and Scordelis [25], or by smeared crack, proposed by Rashid [26]. The objective of discrete crack is to simulate the initiation and propagation of dominant cracks present in the structure. In the case of discrete crack approach, nodes are disassociated due to the presence of cracks and therefore the structure requires frequent renumbering of nodes, which may render the huge computational cost. Nevertheless, when the structure’s behavior has been dominated by only few dominant cracks, the discrete modeling of cracking seems the only choice. On the other hand, the smeared crack approach smears out the cracks over the continuum and captures the deterioration process through the constitutive relationship and reduces the computational cost and time drastically.Crack modeling has gone through several stages due to the advancement in technology and computing facilities. Earlier research work indicates that the formation of crack results in the complete reduction in stresses in the perpendicular direction, thus neglecting the phenomenon called tension stiffening. With the rapid increase in extensive experimental investigations as well as in computing facilities, many finite element codes have been developed for the nonlinear finite element analysis, which incorporates the tension stiffening effect. The first tension stiffening model using degraded concrete modulus was proposed by Scanlon and Murray [27] and subsequently many analytical models have been developed such as Lin and Scordelis [28] model, Vebo and Ghali model, Gilbert and Warner model [29], and Nayal and Rasheed model [30]. The cracks are always assumed to be formed in the direction perpendicular to the direction of the maximum principal stress. These directions may not necessarily remain the same throughout the analysis and loading and hence the modeling of orientation of crack plays a significant role in the response of structure. Still, due to simplicity, many investigations have been performed using fixed crack approach, wherein the direction of principal strain axes may remain fixed throughout the analysis. In this study also, the direction of crack has been considered to be fixed throughout the duration of the analysis. However, the modeling of aggregate interlock has not been taken very seriously. The constant shear retention factor or the simple function has been employed to model the shear transfer across the cracks. Apart from the initiation of crack, the propagation of crack also plays a crucial role in the response of structure. The prediction of crack propagation is a very difficult phenomenon due to scarcity and confliction of test results. Nevertheless, the propagation of cracks plays a crucial role in the response of nonlinear analysis of RC structures. The plain concrete exhibits softening behavior and reinforced concrete exhibits stiffening behavior due to the presence of active reinforcing steel. A gradual release of the concrete stress is adopted in this present study as shown in Figure 4 [31]. The reduction in the stress is given by the following expression:
(29)
E
i
=
α
f
t
′
1
-
ε
i
ε
m
1
ε
i
;
ε
t
≤
ε
i
≤
ε
m
.
α and ε
m are the tension stiffening parameters. ε
m is the maximum value reached by the tensile strain at the point considered. ε
i is the current tensile strain in material direction i. The coefficient depends on the percentage of steel in the section. In the present study, values of α and ε
m are taken as 0.5 and 0.0020, respectively. It has also been reported that the influence of the tension stiffening constants on the response of the structures is generally small and hence the constant value is justified in the analysis [31]. Generally, the cracked concrete can transfer shear forces through dowel action and aggregate interlock. The magnitude of shear moduli has been considerably affected because of extensive cracking in different directions.Figure 4
Tension stiffening effect of cracked concrete.Thus, the reduced shear moduli can be put to incorporate the aggregate interlock and dowel action. In the plain concrete, aggregate interlock is the major shear transfer mechanism and for reinforced concrete, dowel action is the major shear transfer mechanism, with reinforcement ratio being the critical variable. In order to incorporate the aggregate interlock and dowel action, the appropriate value of cracked shear modulus [32] has been considered in the material modeling of concrete.Cracked in One Direction. The stress-strain relationship for cracked concrete where cracking is assumed to take place in only one direction is given as
(30)
σ
1
σ
2
τ
12
τ
13
τ
23
=
0
0
0
0
0
0
E
0
0
0
0
0
G
12
c
0
0
0
0
0
G
13
c
0
0
0
0
0
G
23
ε
1
ε
2
γ
12
γ
13
γ
23
,
G
12
c
=
0.25
×
G
1
-
1
0.004
if
ε
1
<
0.004
,
G
12
c
=
0
if
ε
1
≥
0.004
,
G
13
c
=
G
12
c
;
G
23
=
5
G
6
.Cracked in Two Directions. The stress-strain relationship for cracked concrete where cracking is assumed to take place in both directions is given as
(31)
σ
1
σ
2
τ
12
τ
13
τ
23
=
0
0
0
0
0
0
0
0
0
0
0
0
G
12
c
2
0
0
0
0
0
G
13
c
0
0
0
0
0
G
23
c
ε
1
ε
2
γ
12
γ
13
γ
23
,
G
13
c
=
0.25
×
G
1
-
ε
1
0.004
if
ε
1
<
0.004
,
G
13
c
=
0
if
ε
1
≥
0.004
,
G
23
c
=
0.25
×
G
1
-
ε
2
0.004
,
G
12
c
=
0.5
×
G
13
c
if
G
23
c
<
G
13
c
.
It has also been mentioned by Hinton and Owen [33] that the tensile strength of concrete is a relatively small and unreliable quantity which is not highly influential to the response of structures. In the above stress-strain relationship, the cracked shear modulus (G
c) is assumed to be a function of the current tensile strain. G is the uncracked concrete shear modulus. If the crack closes, the uncracked shear modulus G is assumed in the corresponding direction. Even after the formation of initial cracks, the structure can often deform further without further collapse. In addition to the formation of new cracks, there may be a possibility of crack closing and opening of the existing cracks. If the normal strain across the existing crack becomes greater than that just prior to crack formation, the crack is said to have opened again; otherwise it is assumed to be closed. Nevertheless, if all cracks are closed, then the material is assumed to have gained the status equivalent to that of noncracked concrete with linear elastic behavior.


## 3.2. Concrete Modeling in Compression

The theory of plasticity has been used in the compression modeling of the concrete. The failure surface or bounding surface has been defined to demarcate plastic behavior from the elastic behavior. Failure surface is the important component in the concrete plasticity. Sometimes, the failure surface can be referred to as yield surface or loading surface. The material behaves in the elastic fashion as long as the stress lies below the failure surface. Several failure models have been developed and reported in the literature [34]. Nevertheless, the five-parameter failure model proposed by Willam and Warnke [35] seems to possess all inherent properties of the failure surface. The failure surface is constructed using two meridians, namely, compression meridian and tension meridian. The two meridians are pictorially depicted in a meridian plane and cross section of the failure surface is represented in the deviatoric plane.The variations of the average shear stressesτ
m
t and τ
m
c along tensile (θ
=
0
°) and compressive (θ
=
6
0
°) meridians, as shown in Figure 5, are approximated by second-order parabolic expressions in terms of the average normal stresses σ
m, as follows:
(32)
τ
m
t
f
c
′
=
ρ
t
5
f
c
′
=
a
0
+
a
1
σ
m
f
c
′
+
a
2
σ
m
f
c
′
2
θ
=
0
°
,
τ
m
c
f
c
′
=
ρ
c
5
f
c
′
=
b
0
+
b
1
σ
m
f
c
′
+
b
2
σ
m
f
c
′
2
θ
=
6
0
°
.Willam-Warnke failure model.
(a)
Deviatoric plane in octahedral stress
(b)
Meridian plane in principal stressThese two meridians must intersect the hydrostatic axis at the same pointσ
m
/
f
c
′
=
ξ
-
0 (corresponding to hydrostatic tension); the number of parameters that need to be determined is reduced to five. The five parameters a
0
or
b
0
,
a
1
,
a
2
,
b
1
,
b
2 are to be determined from a set of experimental data, with which the failure surface can be constructed using second-order parabolic expressions. The failure surface is expressed as
(33)
f
σ
m
,
τ
m
,
θ
=
5
τ
m
ρ
σ
m
,
θ
-
1
=
0
,
ρ
θ
=
4
ρ
c
2
-
ρ
t
2
cos
⁡
2
θ
+
5
ρ
t
2
-
4
ρ
t
ρ
c
1
/
2
2
ρ
c
ρ
c
2
-
ρ
t
2
cos
⁡
θ
+
ρ
c
2
ρ
t
-
ρ
c
ρ
c
2
1
/
2
g
g
g
×
4
ρ
c
2
-
ρ
t
2
cos
⁡
2
θ
+
5
ρ
t
2
-
4
ρ
t
ρ
c
1
/
2
×
4
ρ
c
2
-
ρ
t
2
cos
⁡
2
θ
+
ρ
c
-
2
ρ
t
2
-
1
.
The formulation of Willam-Warnke five-parameter material model is described in Chen [34]. Once the yield surface is reached, any further increase in the loading results in the plastic flow. The magnitude and direction of the plastic strain increment are defined using flow rule, which is described in the next section.


### 3.2.1. Flow Rule

In this method, associated flow rule is employed because of the lack of experimental evidence in nonassociated flow rule. The plastic strain increment expressed in terms of current stress increment is given as(34)
d
ε
i
j
p
=
d
λ
∂
f
σ
∂
σ
i
j
.
d
λ determines the magnitude of the plastic strain increment. The gradient ∂
f
σ
/
∂
σ
i
j defines the direction of plastic strain increment to be perpendicular to the yield surface; f
σ is the loading condition or the loading surfaces.


### 3.2.2. Hardening Rule

The relationship between loading surfaces (or effective stress) and the plastic work (accumulated plastic strain) is represented by a hardening rule. The “Madrid parabola” is used to define the hardening rule. The isotropic hardening is adopted in the present study, as shown in Figure6
(35)
σ
=
E
0
ε
-
1
2
E
0
ε
0
ε
2
.
E
0 is the initial elasticity modulus, ε is the total strain, and ε
0 is the total strain at peak stress f
c
′. The total strain can be divided into elastic and plastic components as
(36)
ε
=
ε
e
+
ε
p
,
σ
˙
=
D
e
ε
˙
-
ε
˙
p
,
ε
˙
p
=
λ
˙
a
;
a
=
∂
F
∂
σ
.
F is the yield function and λ
˙ is the consistency parameter, which defines the magnitude of the plastic flow.Figure 6
Isotropic hardening with expanding yield surfaces and the corresponding uniaxial stress-strain curve.The loading, unloading conditions (Kuhn-Tucker conditions) can be stated as(37)
λ
˙
≥
0
;
F
≤
0
;
λ
˙
F
=
0
.
The first of these Khun-Tucker conditions indicates that the consistency parameter is nonnegative; the second condition implies that the stress states must lie on or within the yield surface. The third condition ensures that the stresses lie on the yield surface during the plastic loading
(38)
a
T
σ
˙
=
0
;
a
T
D
e
ε
˙
-
ε
˙
p
=
0
;
a
T
D
e
ε
˙
-
λ
˙
a
=
0
,
λ
˙
=
a
T
D
e
ε
˙
a
T
D
e
a
.
The elastoplastic constitutive matrix is given by the following expression:
(39)
D
e
p
=
D
-
D
a
a
T
D
H
+
a
T
D
a
a
=
flow vector, defined by the stress gradient of the yield function; D = constitutive matrix in elastic range. The second term in (39) represents the effect of degradation of material during the plastic loading.


## 3.2.1. Flow Rule

In this method, associated flow rule is employed because of the lack of experimental evidence in nonassociated flow rule. The plastic strain increment expressed in terms of current stress increment is given as(34)
d
ε
i
j
p
=
d
λ
∂
f
σ
∂
σ
i
j
.
d
λ determines the magnitude of the plastic strain increment. The gradient ∂
f
σ
/
∂
σ
i
j defines the direction of plastic strain increment to be perpendicular to the yield surface; f
σ is the loading condition or the loading surfaces.


## 3.2.2. Hardening Rule

The relationship between loading surfaces (or effective stress) and the plastic work (accumulated plastic strain) is represented by a hardening rule. The “Madrid parabola” is used to define the hardening rule. The isotropic hardening is adopted in the present study, as shown in Figure6
(35)
σ
=
E
0
ε
-
1
2
E
0
ε
0
ε
2
.
E
0 is the initial elasticity modulus, ε is the total strain, and ε
0 is the total strain at peak stress f
c
′. The total strain can be divided into elastic and plastic components as
(36)
ε
=
ε
e
+
ε
p
,
σ
˙
=
D
e
ε
˙
-
ε
˙
p
,
ε
˙
p
=
λ
˙
a
;
a
=
∂
F
∂
σ
.
F is the yield function and λ
˙ is the consistency parameter, which defines the magnitude of the plastic flow.Figure 6
Isotropic hardening with expanding yield surfaces and the corresponding uniaxial stress-strain curve.The loading, unloading conditions (Kuhn-Tucker conditions) can be stated as(37)
λ
˙
≥
0
;
F
≤
0
;
λ
˙
F
=
0
.
The first of these Khun-Tucker conditions indicates that the consistency parameter is nonnegative; the second condition implies that the stress states must lie on or within the yield surface. The third condition ensures that the stresses lie on the yield surface during the plastic loading
(38)
a
T
σ
˙
=
0
;
a
T
D
e
ε
˙
-
ε
˙
p
=
0
;
a
T
D
e
ε
˙
-
λ
˙
a
=
0
,
λ
˙
=
a
T
D
e
ε
˙
a
T
D
e
a
.
The elastoplastic constitutive matrix is given by the following expression:
(39)
D
e
p
=
D
-
D
a
a
T
D
H
+
a
T
D
a
a
=
flow vector, defined by the stress gradient of the yield function; D = constitutive matrix in elastic range. The second term in (39) represents the effect of degradation of material during the plastic loading.


## 3.3. Modeling of Reinforcement in Tension and Compression

Reinforcing bars in structural concrete are generally assumed one-dimensional elements without transverse shear stiffness or flexural rigidity. The reinforcing bar can generally be treated as either discrete or smeared. The major advantage of discrete representation of reinforcing bar is existence of one-to-one correspondence between the real structure and model. In the smeared reinforcement, the average stress-strain relationship is calculated for an element area and incorporated directly as part of the overall concrete element stiffness matrix. In the present investigation, the smeared layered approach is adopted. The bilinear stress-strain curve with linear elastic and strain hardening region is adopted in this study as shown in Figure7. Sometimes, the trilinear idealization has also been adopted as the stress-strain curve in tension. Typically, the hardening strain modulus is assumed to be 1% of initial elasticity modulus. The position and thickness of steel layers are to be defined as input parameters along with the elasticity modulus and hardening modulus. The direction of steel (horizontal or vertical) can be set up by defining the angle with respect to local x-axis. There can only be two states of stress for the reinforcing bar, namely, elastic and linear strain hardening.Figure 7
Stress-strain curve of steel.


## 4. Dynamic Analysis of RC Shear Wall

The dynamic analysis of structure can be performed by three ways, namely, (i) equivalent lateral force method, (ii) response spectrum method, and (iii) time history method. The equivalent lateral force method determines the equivalent dynamic effect in the static manner. The response spectrum method aims at determining the maximum response quantity of the structure. For tall and irregular buildings, dynamic analysis by response spectrum method seems to be a popular choice among designers. The time history analysis of the structure has been successfully used to analyze the structure especially of huge importance. Even though time history analysis consumes time, it is the only method capable of giving results closer to the actual one especially in the nonlinear regime. In the dynamic analysis, the loads are applied over a period of time and the response is obtained at different time intervals. The equation of dynamic equilibrium at any time “t” is given by (1):
(40)
M
U
¨
t
+
C
U
˙
t
+
K
U
t
=
R
t
.
M, C, and K are the mass, damping, and stiffness matrices, respectively. The mass matrix can be formulated either by using consistent mass approach or by using lumped mass approach. Since damping cannot be precisely determined analytically, the damping can be considered proportional to mass or stiffness or both depending on the type of the problem. The direct time integration [24] of the equation of motion can be performed using explicit (central difference scheme) and implicit (Houbolt method, Newmark Beta method, and Wilson Theta method) time integration. In the explicit time integration, the formation of complete stiffness matrix of the structure is not required and hence saves a lot of computer time and money in storing and saving of those data. Moreover, in the case of all explicit time integration schemes, the iterations are not required as the equilibrium at time t
+
Δ
t depends on the equilibrium at time t. Nevertheless, the major drawback of explicit time integration is that the time step (Δ
t) used for calculation of response has to be smaller than the critical time step (Δ
t
cr) to ensure the stable solution
(41)
Δ
t
≤
Δ
t
cr
=
T
n
π
=
2
ω
.
On the other hand, implicit time integration requires the iterations to be carried out within the time step as the solution at time t
+
Δ
t involves the equilibrium equation at t
+
Δ
t. The Newmark β method converges to various implicit and explicit schemes for different values of Beta, called the stability parameter. In this study, for β
=
0.25, the Newmark β method converges to the constant acceleration implicit method, known as trapezoidal rule. The trapezoidal rule is unconditionally stable and hence allows larger time step to be used in the calculation of response. Nevertheless, the time step can be made smaller from the accuracy point of view. The formulation of implicit Newmark Beta method (trapezoidal rule) is mentioned in the literature [24].


### 4.1. Formulation of Mass Matrix

In a dynamic analysis, a correct estimate of mass matrix is very important in predicting the dynamic response of RC structures. There are two different ways, namely, (i) consistent approach and (ii) lumped approach, by which the element mass matrix can be developed. In the case of consistent approach, the masses are assumed to be distributed over the entire finite element mesh. In this approach, the shape functions (N
i) used for the computation of mass matrix are the same shape functions used for the development of stiffness matrix and hence the name “consistent” approach. The mass matrix developed using the consistent approach is known as consistent mass matrix. The consistent element mass matrix (M
e) is given by
(42)
M
e
=
∭
V
ρ
m
N
i
N
i
d
V
,
where ρ
m is the mass density, “i” is the node number, and the integration is performed over the entire volume of the element. The consistent mass matrix contains off-diagonal terms and hence is computationally expensive.On the other hand, the lumped mass matrix is purely diagonal and hence computationally cheaper than the consistent mass matrix. Nevertheless, the diagonalization of the mass matrix from the full mass matrix results in the loss of information and accuracy [18]. Nodal quadrature, row sum, and special lumping are the three lumping procedures available to generate the lumped mass matrices. All the three methods of lumping lead to the same mass matrix for nine-node rectangular elements. Nevertheless, one of the most efficient means of lumping is to distribute the element mass in proportion to the diagonal terms of consistent mass matrix [36] and also discarding the off-diagonal elements. This way of lumping has been successfully used in many finite element codes in practice.The advantage of this special lumping scheme is the assurance of positive definiteness of mass matrix. The use of lumped mass matrix is mostly employed in lower order elements. For higher order elements, the use of lumped mass matrix may not be an appropriate option and hence the present study uses only consistent mass matrix. Moreover, the lumped mass matrix may be an ideal option in the case of RC framed structures in which the masses can be lumped at floor level. As RC shear wall is the concrete structure, it may be appropriate to use consistent mass matrix. The element mass matrices for consistent and lumped mass matrices are as mentioned below:(43)
M
e
c
=
I
1
0
0
0
I
2
0
I
1
0
-
I
2
0
0
0
I
1
0
0
0
-
I
2
0
I
3
0
I
2
0
0
0
I
3
,
(44)
M
e
L
=
I
L
1
0
0
0
0
0
I
L
1
0
0
0
0
0
I
L
1
0
0
0
0
0
I
L
2
0
0
0
0
0
I
L
2
.
Equation (43) is very general and can be used to develop the consistent mass matrix for any displacement based finite element. In (43), “I
1” represents the contribution of mass to resisting the linear or translational motion and is defined as I
1
=
∫
ρ
m
d
z. Hence, masses corresponding to translational degrees of freedom are represented by “I
1.” On the other hand, “I
3” represents the contribution of mass to resisting the change in the rotatory motion and is expressed mathematically as I
3
=
∫
ρ
m
z
×
z
d
z. I
2
=
∫
ρ
m
z
d
z. The total mass matrix M is the sum of the element mass matrices M
e. z is the position of layer middle surface from shell middle surface. Equation (44) represents the lumped mass matrix which is not used in the present analytical study.


### 4.2. Formulation of Damping Matrix

Mass and stiffness matrices can be represented systematically by overall geometry and material characteristics. However, damping can only be represented in a phenomenological manner thus making the dynamic analysis of structures in a state of uncertainty. The quantification and representation of damping is certainly complicated by the relationship between its mathematical representation and the physical sources. The damping may be assumed to be contributed through friction, hysteretic, and viscous characteristics. There is no single universally accepted methodology for representing damping because of the nature of the state variables, which control damping. Nevertheless, several investigations have been done in making the representation of damping in a simplistic yet logical manner [37]. Only for the mathematical convenience, the damping has been modeled as equivalent viscous damping, represented as the percentage of critical damping. The governing equation of motion second-order differential equation with constant coefficients is rewritten as
(45)
M
u
¨
+
C
u
˙
+
K
u
=
R
t
.
The trial solution is given by
(46)
u
=
c
e
s
t
.On substituting the trial solution and simplifying, the roots of quadratic equation are as(47a)
s
=
-
c
2
m
±
c
2
m
2
-
ω
2
,
c
2
m
2
-
ω
2
>
0
⟹
2
real
roots
over
damped
⟶
s
1
=
-
ξ
ω
+
ω
ξ
2
-
1
;
s
2
=
-
ξ
ω
-
ω
ξ
2
-
1
,
c
2
m
2
-
ω
2
=
0
⟹
1
real
root
critically
damped
⟶
s
1
=
-
ξ
ω
,
c
2
m
2
-
ω
2
<
0
⟹
2
complex
roots
under
damped
⟶
s
1
=
-
ξ
ω
+
i
ω
1
-
ξ
2
;
s
2
=
-
ξ
ω
-
i
ω
1
-
ξ
2
.
In the above equation, the damping ratio ξ is given by
(47b)
ξ
=
c
2
m
ω
⟹
c
2
m
=
ξ
ω
;
ξ
=
c
c
cr
.Over damped system does not vibrate it all. The classical example is automatic door closing. Critical damping has the linear exploding function and hence the amplitude is higher than the over damped system, followed by exponentially exploding function resulting in the fast movement over the time. The bottom line is that both over damped and critically damped systems do not vibrate at all. Nevertheless, in a building structure, critically damped and over damped situation may not arise. Damping matrix can be formulated analogous to mass and stiffness matrices [38, 39]. It is also important to note that the damping matrix should be formulated from damping ratio and not from the member sizes. Rayleigh dissipation function assumes that the dissipation of energy takes place and can be idealized as the function of velocity. When Rayleigh damping is used, the resultant damping matrix is of the same size as stiffness matrix. Rayleigh damping is being used conveniently because of its versatility in segregating each mode independently. The damping can be defined as the linear combination of mass and stiffness matrices as
(48)
C
=
α
M
+
β
K
,
(49)
ς
i
=
α
2
ω
i
+
β
ω
i
2
.
It is to be noted that the damping is controlled by only two parameters (Figure 8). From (49), it is observed that if β is zero, the higher modes of the structure will be assigned very little damping. When α alpha is zero, the higher modes will be heavily damped, as the damping ratio is directly proportional to circular frequency (ω) [40]. Thus, the choice of damping is problem dependent. Hence, it is inevitable to perform modal analysis to determine the different frequencies for different modes. However, in the present study, the fundamental natural period (T) and fundamental natural frequency (f) have been calculated using the formulas as mentioned in Goel and Chopra [41]. The coupling of the modes usually can be avoided easily in the case of undamped free vibration. The same is not true for damped vibration. Hence in order to represent the equation of motion in uncoupled form, it is suggested to have a damping matrix proportional to uncoupled mass and stiffness matrices. Thus, Rayleigh’s proportional damping has the specific advantage that the equation of motion can be uncoupled when it is proportional to mass and stiffness matrices. Thus, it is proposed to use Rayleigh damping in this study.Figure 8
Variation of damping with circular frequency.


### 4.3. Nonlinear Solution

The numerical procedure for nonlinear analysis employs the iterative procedure to satisfy the equilibrium at the end of the load step. Once the convergence of the solution is achieved, the algorithm proceeds to the next step. It is always desirable to keep the load step very small especially after the onset of nonlinear behavior. The stiffness matrix is updated at the beginning of each load step. The convergence is said to be achieved if the out of balance forces, calculated as follows, are less than the specified tolerance:(50)
ψ
i
n
=
f
n
-
p
i
n
=
f
n
-
∫
V
B
T
σ
i
n
d
V
<
Tolerance
0.0025
.


## 4.1. Formulation of Mass Matrix

In a dynamic analysis, a correct estimate of mass matrix is very important in predicting the dynamic response of RC structures. There are two different ways, namely, (i) consistent approach and (ii) lumped approach, by which the element mass matrix can be developed. In the case of consistent approach, the masses are assumed to be distributed over the entire finite element mesh. In this approach, the shape functions (N
i) used for the computation of mass matrix are the same shape functions used for the development of stiffness matrix and hence the name “consistent” approach. The mass matrix developed using the consistent approach is known as consistent mass matrix. The consistent element mass matrix (M
e) is given by
(42)
M
e
=
∭
V
ρ
m
N
i
N
i
d
V
,
where ρ
m is the mass density, “i” is the node number, and the integration is performed over the entire volume of the element. The consistent mass matrix contains off-diagonal terms and hence is computationally expensive.On the other hand, the lumped mass matrix is purely diagonal and hence computationally cheaper than the consistent mass matrix. Nevertheless, the diagonalization of the mass matrix from the full mass matrix results in the loss of information and accuracy [18]. Nodal quadrature, row sum, and special lumping are the three lumping procedures available to generate the lumped mass matrices. All the three methods of lumping lead to the same mass matrix for nine-node rectangular elements. Nevertheless, one of the most efficient means of lumping is to distribute the element mass in proportion to the diagonal terms of consistent mass matrix [36] and also discarding the off-diagonal elements. This way of lumping has been successfully used in many finite element codes in practice.The advantage of this special lumping scheme is the assurance of positive definiteness of mass matrix. The use of lumped mass matrix is mostly employed in lower order elements. For higher order elements, the use of lumped mass matrix may not be an appropriate option and hence the present study uses only consistent mass matrix. Moreover, the lumped mass matrix may be an ideal option in the case of RC framed structures in which the masses can be lumped at floor level. As RC shear wall is the concrete structure, it may be appropriate to use consistent mass matrix. The element mass matrices for consistent and lumped mass matrices are as mentioned below:(43)
M
e
c
=
I
1
0
0
0
I
2
0
I
1
0
-
I
2
0
0
0
I
1
0
0
0
-
I
2
0
I
3
0
I
2
0
0
0
I
3
,
(44)
M
e
L
=
I
L
1
0
0
0
0
0
I
L
1
0
0
0
0
0
I
L
1
0
0
0
0
0
I
L
2
0
0
0
0
0
I
L
2
.
Equation (43) is very general and can be used to develop the consistent mass matrix for any displacement based finite element. In (43), “I
1” represents the contribution of mass to resisting the linear or translational motion and is defined as I
1
=
∫
ρ
m
d
z. Hence, masses corresponding to translational degrees of freedom are represented by “I
1.” On the other hand, “I
3” represents the contribution of mass to resisting the change in the rotatory motion and is expressed mathematically as I
3
=
∫
ρ
m
z
×
z
d
z. I
2
=
∫
ρ
m
z
d
z. The total mass matrix M is the sum of the element mass matrices M
e. z is the position of layer middle surface from shell middle surface. Equation (44) represents the lumped mass matrix which is not used in the present analytical study.


## 4.2. Formulation of Damping Matrix

Mass and stiffness matrices can be represented systematically by overall geometry and material characteristics. However, damping can only be represented in a phenomenological manner thus making the dynamic analysis of structures in a state of uncertainty. The quantification and representation of damping is certainly complicated by the relationship between its mathematical representation and the physical sources. The damping may be assumed to be contributed through friction, hysteretic, and viscous characteristics. There is no single universally accepted methodology for representing damping because of the nature of the state variables, which control damping. Nevertheless, several investigations have been done in making the representation of damping in a simplistic yet logical manner [37]. Only for the mathematical convenience, the damping has been modeled as equivalent viscous damping, represented as the percentage of critical damping. The governing equation of motion second-order differential equation with constant coefficients is rewritten as
(45)
M
u
¨
+
C
u
˙
+
K
u
=
R
t
.
The trial solution is given by
(46)
u
=
c
e
s
t
.On substituting the trial solution and simplifying, the roots of quadratic equation are as(47a)
s
=
-
c
2
m
±
c
2
m
2
-
ω
2
,
c
2
m
2
-
ω
2
>
0
⟹
2
real
roots
over
damped
⟶
s
1
=
-
ξ
ω
+
ω
ξ
2
-
1
;
s
2
=
-
ξ
ω
-
ω
ξ
2
-
1
,
c
2
m
2
-
ω
2
=
0
⟹
1
real
root
critically
damped
⟶
s
1
=
-
ξ
ω
,
c
2
m
2
-
ω
2
<
0
⟹
2
complex
roots
under
damped
⟶
s
1
=
-
ξ
ω
+
i
ω
1
-
ξ
2
;
s
2
=
-
ξ
ω
-
i
ω
1
-
ξ
2
.
In the above equation, the damping ratio ξ is given by
(47b)
ξ
=
c
2
m
ω
⟹
c
2
m
=
ξ
ω
;
ξ
=
c
c
cr
.Over damped system does not vibrate it all. The classical example is automatic door closing. Critical damping has the linear exploding function and hence the amplitude is higher than the over damped system, followed by exponentially exploding function resulting in the fast movement over the time. The bottom line is that both over damped and critically damped systems do not vibrate at all. Nevertheless, in a building structure, critically damped and over damped situation may not arise. Damping matrix can be formulated analogous to mass and stiffness matrices [38, 39]. It is also important to note that the damping matrix should be formulated from damping ratio and not from the member sizes. Rayleigh dissipation function assumes that the dissipation of energy takes place and can be idealized as the function of velocity. When Rayleigh damping is used, the resultant damping matrix is of the same size as stiffness matrix. Rayleigh damping is being used conveniently because of its versatility in segregating each mode independently. The damping can be defined as the linear combination of mass and stiffness matrices as
(48)
C
=
α
M
+
β
K
,
(49)
ς
i
=
α
2
ω
i
+
β
ω
i
2
.
It is to be noted that the damping is controlled by only two parameters (Figure 8). From (49), it is observed that if β is zero, the higher modes of the structure will be assigned very little damping. When α alpha is zero, the higher modes will be heavily damped, as the damping ratio is directly proportional to circular frequency (ω) [40]. Thus, the choice of damping is problem dependent. Hence, it is inevitable to perform modal analysis to determine the different frequencies for different modes. However, in the present study, the fundamental natural period (T) and fundamental natural frequency (f) have been calculated using the formulas as mentioned in Goel and Chopra [41]. The coupling of the modes usually can be avoided easily in the case of undamped free vibration. The same is not true for damped vibration. Hence in order to represent the equation of motion in uncoupled form, it is suggested to have a damping matrix proportional to uncoupled mass and stiffness matrices. Thus, Rayleigh’s proportional damping has the specific advantage that the equation of motion can be uncoupled when it is proportional to mass and stiffness matrices. Thus, it is proposed to use Rayleigh damping in this study.Figure 8
Variation of damping with circular frequency.


## 4.3. Nonlinear Solution

The numerical procedure for nonlinear analysis employs the iterative procedure to satisfy the equilibrium at the end of the load step. Once the convergence of the solution is achieved, the algorithm proceeds to the next step. It is always desirable to keep the load step very small especially after the onset of nonlinear behavior. The stiffness matrix is updated at the beginning of each load step. The convergence is said to be achieved if the out of balance forces, calculated as follows, are less than the specified tolerance:(50)
ψ
i
n
=
f
n
-
p
i
n
=
f
n
-
∫
V
B
T
σ
i
n
d
V
<
Tolerance
0.0025
.


## 5. Development of Computer Program

In the present study an analysis module, NLDAS, was developed using Fortran 77 and used to perform the nonlinear dynamic finite element analysis of RC shear walls. The computer codes developed by Huang [23] and Owen and Hinton [31] for the static and dynamic elastoplastic analysis of RC structures based on simple Owen-Figurious yield/failure criterion have been taken as the base programs in this study. These programs have been merged and modified to include state-of-the-art five-parameter yield/failure model and concrete cracking. A modular approach has been adopted for the program development. To this end several new subroutines have been incorporated and few subroutines have been enhanced.The program consists of 19 subroutines which are developed to perform various operations. In the program NLDAS, though the input files for static and dynamic analysis are different, all the data sets have to be read in at the beginning of the program. The program mainly consisted of the following subroutines: input, loading, incremental loading, stiffness, mass and damping matrices development and assembly, solution of equations, residual force calculations, and convergence check and output results, in addition to modules for storage of global arrays such as nodal coordinates, element connectivity, material properties, and boundary conditions. Figure9 shows the program layout explaining the process of the program.Figure 9
Program layout, NLDAS.


## 6. Displacement Time History Response of Shear Wall with Different Opening Locations

In order to identify the safe regions where the openings can be provided in a shear wall, two representative problems of ten storeys (35 m high, 8 m wide, and 0.3 thick) and five storeys (17.5 m high, 8 m wide, and 0.3 thick) as shown in Figures10(a) and 10(b), behaviorally slender and squat type, are chosen and analyzed for dynamic loading condition subjected to Figure 10(c), using the finite element analysis. The degenerated shell element of size 1 m × 0.5 m has been used to discretize the shear wall geometry. The opening location is varied keeping all practical positions. In total, there are seven cases considered as the possible opening locations prevailed in practice, namely, (i) door cum window, (ii) two doors, (iii) two windows, (iv) three windows, (v) four windows, (vi) staggered two windows, and (vii) staggered four windows. The size of the opening in all cases is kept at 14%. The openings are provided uniformly in all the storeys. For simplicity, only a typical storey is represented in Figure 10. The material properties used for the material modeling are as mentioned in Table 1.Table 1
Material property of RC shear wall.
Material property
Magnitude
Young’s modulus of concrete
2.7 × 1010 kN/m2
Young’s modulus of steel
2.0 × 1011 kN/m2
Poisson’s ratio of concrete
0.100
Uniaxial compressive strength of concrete
30 × 106 kN/m2
Tensile strength of concrete
3 × 106 kN/m2
Ultimate crushing strain of concrete
0.0035
Yield stress of steel
50 × 107 kN/m2
Tension stiffening coefficient (α
ts
)
0.6
Tension stiffening constant (ε
ts)
0.0020(a) Slender shear wall with openings. (b) Squat shear wall with openings. (c) Dynamic ground acceleration applied at the base of the structure.
(a)
(b)
(c)The analysis has been done for all practical damping cases, namely, (i) no damping, (iii) 2% damping, (iii) 5% damping, and (iv) 10% damping. The reinforcement ratio of 0.25% is adopted for both vertical and horizontal reinforcement. The strengthening of shear walls around openings was considered as per IS: 13920 requirements. The simulated EL-Centro earthquake, shown in Figure8, is applied as the input ground acceleration with maximum amplitude of 1.05 g. The response of the structure is traced for 1.5 seconds of duration. Several investigators have also adopted this way of response calculation by predicting the response only for the most intense earthquake period in order to simplify the computation. The total duration of the ground motion has been taken as 1.5 seconds. The analysis has been carried out with the time step of 0.015 seconds. The undamped top displacement time histories of shear walls with different opening locations have been presented in Table 2. Though the response analyses of shear walls have been conducted for different damping ratios, only undamped responses have been plotted for brevity.Table 2
Displacement time history responses of shear walls with different door window opening locations.
A typical storey
Slender shear wall (10 storeys)
Squat shear wall (5 storeys)Tables3 and 4 show the influence of opening locations on the maximum displacement response of slender and squat shear walls, respectively.Table 3
Influence of opening locations on the maximum displacement of slender shear walls.
Case
Maximum displacement of shear wall at
0% damping
2% damping
5% damping
10% damping
Door cum window
135
80.04
3.77
3.76
2.55
2.55
1.60
1.60
Two doors
114
80
3.04
3.03
2.12
2.12
1.37
1.37
Two windows
167.8
83.23
3.90
3.90
2.58
2.59
1.63
1.63
Three windows
185.1
87.44
3.78
3.78
2.73
2.73
1.77
1.77
Four windows
99.46
61.68
3.02
3.02
2.13
2.13
1.39
1.39
Staggered two windows
161.3
16.7
3.86
3.87
2.48
2.48
1.55
1.55
Staggered four windows
126.5
56.83
3.47
3.46
2.37
2.37
1.50
1.50Table 4
Influence of opening locations on the maximum displacement of squat shear walls.
Case
Maximum displacement of shear wall at
0% damping
2% damping
5% damping
10% damping
Door cum window
16.20
11.84
1.38
1.39
0.64
0.64
0.39
0.39
Two doors
16.20
8.14
1.31
1.31
0.53
0.53
0.34
0.34
Two windows
35.44
12.91
1.59
1.59
0.65
0.65
0.40
0.40
Three windows
24.68
12.27
1.66
1.66
0.69
0.69
0.43
0.43
Four windows
8.28
8.06
1.30
1.30
0.53
0.53
0.34
0.34
Staggered two windows
32.92
31.05
1.52
1.52
0.61
0.61
0.38
0.38
Staggered four windows
9.02
9.00
1.41
1.41
0.59
0.59
0.37
0.37
Not strengthened  strengthened.As evidenced through Tables3 and 4, the displacement has been found to be higher in the case of slender shear walls than squat shear walls. It was inferred from Tables 2 and 3 that the maximum displacement responses of slender shear walls have been found to be in the case of three window openings followed by two window openings. Nevertheless, the strengthening around openings reduces the displacement response to 53%. The influence of strengthening has been considered massive in the case of staggered openings (two windows and four windows), as seen in Table 3. In general, the strengthening has resulted in better behavior of the shear wall.On the other hand, in the case of squat shear walls (Tables2 and 4), the displacement response has been found to be least in the case of shear wall with four windows, regular and staggered. Also from Table 4, it is further observed that the squat shear walls undergo high displacement in the presence of two window openings. Hence, it is better to provide large number of small openings to have better structural performance.Table5 shows the maximum displacement response on the RC shear wall with different opening locations for different damping ratios. It is observed that the influence of strengthening has been found to be significant in case of shear wall with no damping. On the other hand, the influence of strengthening has been offset by the presence of damping as seen in Table 2 for 2%, 5%, and 10%.Table 5
Maximum displacement response on RC shear walls with different opening locations for different damping ratios.
Case
Slender shear wall (10 storeys)
Squat shear wall (5 storeys)
No damping
2% damping
5% damping
10% damping


## 7. Conclusions

On the basis of displacement time history responses of shear walls with different opening locations and with various damping ratios, it has been concluded that shear walls are penetrated by large number of small openings than small number of large openings. Moreover, the influence of strengthening has been considered essential especially for undamped shear wall. The shear wall with four windows has been considered best for both slender and squat shear walls. The strengthening (ductile detailing) has been considered significant in the case of slender shear wall with staggered openings. However, for higher damping, the strengthening has not been considered essential.



---

*Source: 858510-2014-11-25.xml*