See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/253182180

Finite element analysis of reinforced concrete structures under monotonic loads

Article · January 1990

CITATIONS READS
132 2,496

2 authors, including:

Filip C. Filippou
University of California, Berkeley
87 PUBLICATIONS   3,768 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Decision-Making Guidelines for Overcoming Non-Convergence in Nonlinear Dynamic Analysis of Reinforced Concrete Bridges View project

Mixed formulations for nonlinear static and dynamic analysis of cables View project

All content following this page was uploaded by Filip C. Filippou on 20 July 2015.

The user has requested enhancement of the downloaded file.

REPORT NO. STRUCTURAL ENGINEERING
UCB/SEMM-90/14 MECHANICS AND MATERIALS

FINITE ELEMENT ANALYSIS OF

REINFORCED CONCRETE STRUCTURES
UNDER MONOTONIC LOADS

BY

HYO-GYOUNG KWAK

FILIP C. FILIPPOU

NOVEMBER 1990 DEPARTMENT OF CIVIL ENGINEERING

UNIVERSITY OF CALIFORNIA
BERKELEY, CALIFORNIA
 FINITE ELEMENT ANALYSIS OF
REINFORCED CONCRETE STRUCTURES
UNDER MONOTONIC LOADS

by

H. G. Kwak
Visiting Scholar
Korea Advanced Institute of Science and Technology

and

Filip C. Filippou
Associate Professor
Department of Civil Engineering

A Report on Research Conducted under

Grant RTA-59M848 from the
California Department of Transportation

Report No. UCB/SEMM-90/14

Structural Engineering, Mechanics and Materials
Department of Civil Engineering
University of California, Berkeley

November 1990
 ABSTRACT

This study deals with the finite element analysis of the monotonic behavior of
reinforced concrete beams, slabs and beam-column joint subassemblages. It is assumed that
the behavior of these members can be described by a plane stress field. Concrete and
reinforcing steel are represented by separate material models which are combined together
with a model of the interaction between reinforcing steel and concrete through bond-slip to
describe the behavior of the composite reinforced concrete material. The material behavior of
concrete is described by two failure surfaces in the biaxial stress space and one failure surface
in the biaxial strain space. Concrete is assumed as a linear elastic material for stress states
which lie inside the initial yield surface. For stresses outside this surface the behavior of
concrete is described by a nonlinear orthotropic model, whose axes of orthotropy are parallel
to the principal strain directions. The concrete stress-strain relation is derived from equivalent
uniaxial relations in the axes of orthotropy. The behavior of cracked concrete is described by
a system of orthogonal cracks, which follow the principal strain directions and are thus
rotating during the load history. Crushing or cracking of concrete takes place when the strains
lie outside the ultimate surface in the biaxial strain space.
A new smeared finite element model is proposed based on an improved cracking
criterion, which is derived from fracture mechanics principles. This model retains objectivity
of the results for very large finite elements, since it considers cracking to be concentrated
over a small region around the integration point and not over the entire finite element, as do
previous models.
A new reinforcing steel model which is embedded inside a concrete element, but
accounts for the effect of bond-slip is developed. This model results in significant savings in
the number of nodes needed to account for the effect of bond-slip, particularly, in three
dimensional finite element models. A new nonlinear solution scheme is developed in
connection with this model.
Finally, correlation studies between analytical and experimental results and several
parameter studies are conducted with the objective to establish the validity of the proposed
models and identify the significance of various effects on the local and global response of
reinforced concrete members. These studies show that the effects of tension-stiffening and
bond-slip are very important and should always be included in finite element models of the

i
response of reinforced concrete members. On the other hand, parameters, such as the tensile
strength of concrete and the value of the cracked shear constant, do not seem to affect the
response of slender beams in bending.

ii
 ACKNOWLEDGEMENTS

This report is part of a larger study on the seismic behavior of reinforced concrete
highway structures supported by Grant No. RTA-59M848 from the California Department of
Transportation (CALTRANS). This support is gratefully acknowledged. Any opinions
expressed in this report are those of the authors and do not reflect the views of the sponsoring
agency.

iii
 TABLE OF CONTENTS

ABSTRACT ....................................................................................................................................... i
ACKNOWLEDGEMENTS ...............................................................................................................ii
TABLE OF CONTENTS.................................................................................................................. iii
LIST OF TABLES ..............................................................................................................................iv
LIST OF FIGURES .......................................................................................................................... vii
CHAPTER 1 INTRODUCTION.................................................................................................... 1
1.1 General .............................................................................................................. 1

1.2 Literature Survey............................................................................................... 3

1.3 Objectives and Scope ........................................................................................ 8

CHAPTER 2 MATERIAL MODELING ...................................................................................... 11

2.1 General ............................................................................................................ 11

2.2 Concrete .......................................................................................................... 13

2.2.1 Behavior of Concrete........................................................................... 13
2.2.2 Concrete Models ................................................................................. 13
2.2.3 Behavior of Cracked Concrete ............................................................ 20
2.2.3.1 Description of a Cracked Section............................................ 20
2.2.3.2 Crack Models .......................................................................... 23
2.2.3.3 Proposed Model....................................................................... 28
2.2.3.4 Concrete Material Matrix ........................................................ 31

2.3 Reinforcing Steel............................................................................................. 34

2.3.1 Behavior of Reinforcing Steel............................................................. 34
2.3.2 Steel Material Matrix .......................................................................... 37

2.4 Bond between Reinforcing Steel and Concrete............................................... 43

2.4.1 Bond Behavior..................................................................................... 43
2.4.2 Bond Stiffness Matrix ......................................................................... 47

iv
 2.5 Discrete Reinforcing Steel Model with Bond-Slip.......................................... 49

CHAPTER 3 FINITE ELEMENT MODELING OF RC BEAMS AND SLABS ......................55

3.1 General ............................................................................................................ 55

3.2 Theoretical Considerations.............................................................................. 56

3.3 Finite Element Representation ........................................................................ 59

3.3.1 Reinforced Concrete Beams................................................................ 61
3.3.2 Reinforced Concrete Slabs .................................................................. 62

3.4 Numerical Implementation.............................................................................. 67

3.4.1 Iteration Method .................................................................................. 67
3.4.2 Solution Algorithm.............................................................................. 69
3.4.3 Solution Algorithm for Reinforcing Bar with Bond-Slip.................... 74
3.4.4 Convergence Criterion ........................................................................ 78

CHAPTER 4 APPLICATIONS.....................................................................................................81
4.1 Introduction ..................................................................................................... 81

4.2 Anchored Reinforcing Bar Under Monotonic and Cyclic Loads .................... 81

4.3 Reinforced Concrete Beams............................................................................ 88

4.4 Beam to Column Joint Subassemblage ........................................................... 99

4.5 Reinforced Concrete Slabs ............................................................................ 102

CHAPTER 5 CONCLUSIONS ...................................................................................................107

REFERENCES..................................................................................................................................109

v
 CHAPTER 1
INTRODUCTION

1.1 General

Reinforced concrete (RC) has become one of the most important building materials
and is widely used in many types of engineering structures. The economy, the efficiency, the
strength and the stiffness of reinforced concrete make it an attractive material for a wide
range of structural applications. For its use as structural material, concrete must satisfy the
following conditions:

(1) The structure must be strong and safe. The proper application of the fundamental
principles of analysis, the laws of equilibrium and the consideration of the mechanical
properties of the component materials should result in a sufficient margin of safety
against collapse under accidental overloads.

(2) The structure must be stiff and appear unblemished. Care must be taken to control
deflections under service loads and to limit the crack width to an acceptable level.

(3) The structure must be economical. Materials must be used efficiently, since the
difference in unit cost between concrete and steel is relatively large.

The ultimate objective of design is the creation of a safe and economical structure.
Advanced analytical tools can be an indispensable aid in the assessment of the safety and the
serviceability of a proposed design. This is, especially, true for many complex modern
structures such as nuclear power plants, bridges, off-shore platforms for oil and gas
exploration and underground or underwater tunnels, which are subjected to very complex
load histories. The safety and serviceability assessment of these structures necessitates the
development of accurate and reliable methods and models for their analysis. In addition, the
rise in cost of structures encourages engineers to seek more economical alternative designs
often resorting to innovative construction methods without lowering the safety of the
structure. Intimately related to the increase in scale of modern structures is the extent and
impact of disaster in terms of human and economical loss in the event of structural failure. As
a result, careful and detailed structural safety analysis becomes more and more necessary. The

1
2 CHAPTER 1

objective of such an analysis is the investigation of the behavior of the structure under all
possible loading conditions, both, monotonic and cyclic, its time-dependent behavior, and,
especially, its behavior under overloading.

Reinforced concrete structures are commonly designed to satisfy criteria of

serviceability and safety. In order to ensure the serviceability requirement it is necessary to
predict the cracking and the deflections of RC structures under service loads. In order to
assess the margin of safety of RC structures against failure an accurate estimation of the
ultimate load is essential and the prediction of the load-deformation behavior of the structure
throughout the range of elastic and inelastic response is desirable.

Within the framework of developing advanced design and analysis methods for
modern structures the need for experimental research continues. Experiments provide a firm
basis for design equations, which are invaluable in the preliminary design stages.
Experimental research also supplies the basic information for finite element models, such as
material properties. In addition, the results of finite element models have to be evaluated by
comparing them with experiments of full-scale models of structural subassemblages or, even,
entire structures. The development of reliable analytical models can, however, reduce the
number of required test specimens for the solution of a given problem, recognizing that tests
are time-consuming and costly and often do not simulate exactly the loading and support
conditions of the actual structure.

The development of analytical models of the response of RC structures is complicated

by the following factors:

• Reinforced concrete is a composite material made up of concrete and steel, two

materials with very different physical and mechanical behavior;
• Concrete exhibits nonlinear behavior even under low level loading due to nonlinear
material behavior, environmental effects, cracking, biaxial stiffening and strain
softening;

• Reinforcing steel and concrete interact in a complex way through bond-slip and
aggregate interlock.

These complex phenomena have led engineers in the past to rely heavily on empirical
formulas for the design of concrete structures, which were derived from numerous
experiments. With the advent of digital computers and powerful methods of analysis, such as
the finite element method, many efforts to develop analytical solutions which would obviate
CHAPTER 1 3

the need for experiments have been undertaken by investigators. The finite element method
has thus become a powerful computational tool, which allows complex analyses of the
nonlinear response of RC structures to be carried out in a routine fashion. With this method
the importance and interaction of different nonlinear effects on the response of RC structures
can be studied analytically. The present study is part of this continuing effort and concerns
the analysis of reinforced concrete beams, slabs, and beam-to-column subassemblages under
monotonic loads. A follow-up study will address the response of these structures under cyclic
load reversals.

1.2 Literature Survey

A brief review of previous studies on the application of the finite element method to
the analysis of reinforced concrete structures is presented is this section. A more detailed
description of the underlying theory and the application of the finite element method to the
analysis of linear and nonlinear reinforced concrete structures is presented in excellent state-
of-the-art reports by the American Society of Civil Engineers in 1982 (ASCE 1982) and 1985
(Meyer and Okamura, eds. 1985).

The earliest publication on the application of the finite element method to the analysis
of RC structures was presented by Ngo and Scordelis (1967). In their study, simple beams
were analyzed with a model in which concrete and reinforcing steel were represented by
constant strain triangular elements, and a special bond link element was used to connect the
steel to the concrete and describe the bond-slip effect. A linear elastic analysis was performed
on beams with predefined crack patterns to determine principal stresses in concrete, stresses
in steel reinforcement and bond stresses. Since the publication of this pioneering work, the
analysis of reinforced concrete structures has enjoyed a growing interest and many
publications have appeared. Scordelis et al. (1974) used the same approach to study the effect
of shear in beams with diagonal tension cracks and accounted for the effect of stirrups, dowel
shear, aggregate interlock and horizontal splitting along the reinforcing bars near the support.

Nilson (1972) introduced nonlinear material properties for concrete and steel and a
nonlinear bond-slip relationship into the analysis and used an incremental load method of
nonlinear analysis. Four constant strain triangular elements were combined to form a
quadrilateral element by condensing out the central node. Cracking was accounted for by
stopping the solution when an element reached the tensile strength, and reloading
4 CHAPTER 1

incrementally after redefining a new cracked structure. The method was applied to concentric
and eccentric reinforced concrete tensile members which were subjected to loads applied at
the end of the reinforcing bars and the results were compared with experimental data.

Franklin (1970) advanced the capabilities of the analytical method by developing a

nonlinear analysis which automatically accounted for cracking within finite elements and the
redistribution of stresses in the structure. This made it possible to trace the response of two-
dimensional systems from initial loading to failure in one continuous analysis. Incremental
loading with iterations within each increment was used to account for cracking in the finite
elements and for the nonlinear material behavior. Franklin used special frame-type elements,
quadrilateral plane stress elements, axial bar members, two-dimensional bond links and tie
links to study reinforced concrete frames and RC frames coupled with shear walls.

Plane stress elements were used by numerous investigators to study the behavior of
reinforced concrete frame and wall systems. Nayak and Zienkiewicz (1972) conducted two-
dimensional stress studies which include the tensile cracking and the elasto-plastic behavior
of concrete in compression using an initial stress approach. Cervenka (1970) analyzed shear
walls and spandrel beams using an initial stress approach in which the elastic stiffness matrix
at the beginning of the entire analysis is used in all iterations. Cervenka proposed a
constitutive relationship for the composite concrete-steel material through the uncracked,
cracked and plastic stages of behavior.

For the analysis of RC beams with material and geometric nonlinearities Rajagopal
(1976) developed a layered rectangular plate element with axial and bending stiffness in
which concrete was treated as an orthotropic material. RC beam and slab problems have also
been treated by many other investigators (Lin and Scordelis 1975; Bashur and Darwin 1978;
Rots et al. 1985; Barzegar and Schnobrich 1986; Adeghe and Collins 1986; Bergmann and
Pantazopoulou 1988; Cervenka et al. 1990; Kwak 1990) using similar methods.

Selna (1969) analyzed beams and frames made up of one-dimensional elements with
layered cross sections which accounted for progressive cracking and changing material
properties through the depth of the cross section as a function of load and time. Significant
advances and extensions of the finite element analysis of reinforced concrete beams and
frames to include the effects of heat transfer due to fire, as well as the time-dependent effects
of creep and shrinkage, were made by Becker and Bresler (1974).
CHAPTER 1 5

Two basically different approaches have been used so far for the analysis of RC slabs
by the finite element method: the modified stiffness approach and the layer approach. The
former is based on an average moment-curvature relationship which reflects the various
stages of material behavior, while the latter subdivides the finite element into imaginary
concrete and steel layers with idealized stress-strain relations for concrete and reinforcing
steel.

Experimental and analytical studies of RC slabs were conducted by Joffriet and

McNeice (1971). The analyses were based on a bilinear moment-curvature relation which
was derived from an empirically determined effective moment of inertia of the cracked slab
section including the effect of tension stiffening. The change in bending stiffness of the
elements due to cracking normal to the principal moment direction is accounted for by
reducing the flexural stiffness of the corresponding element.

Dotroppe et al. (1973) used a layered finite element procedure in which slab elements
were divided into layers to account for the progressive cracking through the slab thickness.
Scanlon and Murray (1974) have developed a method of incorporating both cracking and
time-dependent effects of creep and shrinkage in slabs. They used layered rectangular slab
elements which could be cracked progressively layer by layer, and assumed that cracks
propagate only parallel and perpendicular to orthogonal reinforcement. Lin and Scordelis
(1975) utilized layered triangular finite elements in RC shell analysis and included the
coupling between membrane and bending effects, as well as the tension stiffening effect of
concrete between cracks in the model.

The finite element analysis of an axisymmetric solid under axisymmetric loading can
be readily reduced to a two-dimensional analysis. Bresler and Bertero (1968) used an
axisymmetric model to study the stress distribution in a cylindrical concrete specimen
reinforced with a single plain reinforcing bar. The specimen was loaded by applying tensile
loads at the ends of the bar.

In one of the pioneering early studies Rashid (1968) introduced the concept of a
"smeared" crack in the study of the axisymmetric response of prestressed concrete reactor
structures. Rashid took into account cracking and the effects of temperature, creep and load
history in his analyses. Today the smeared crack approach of modeling the cracking behavior
of concrete is almost exclusively used by investigators in the nonlinear analysis of RC
structures, since its implementation in a finite element analysis program is more
straightforward than that of the discrete crack model. Computer time considerations also
6 CHAPTER 1

favor the smeared crack model in analyses which are concerned with the global response of
structures. At the same time the concerted effort of many investigators in the last 20 years has
removed many of the limitations of the smeared crack model (ASCE 1982; Meyer and
Okamura, eds. 1985).

Gilbert and Warner (1978) used the smeared crack model and investigated the effect
of the slope of the descending branch of the concrete stress-strain relation on the behavior of
RC slabs. They were among the first to point out that analytical results of the response of
reinforced concrete structures are greatly influenced by the size of the finite element mesh
and by the amount of tension stiffening of concrete. Several studies followed which
corroborated these findings and showed the effect of mesh size (Bazant and Cedolin 1980;
Bazant and Oh 1983; Kwak 1990) and tension stiffening (Barzegar and Schnobrich 1986;
Leibengood et al. 1986) on the accuracy of finite element analyses of RC structures with the
smeared crack model. In order to better account for the tension stiffening effect of concrete
between cracks some investigators have artificially increased the stiffness of reinforcing steel
by modifying its stress-strain relationship (Gilbert and Warner 1977). Others have chosen to
modify the tensile stress-strain curve of concrete by including a descending post-peak branch
(Lin and Scordelis 1975; Vebo and Ghali 1977; Barzegar and Schnobrich 1986; Abdel
Rahman and Hinton 1986).

In the context of the smeared crack model two different representations have
emerged: the fixed crack and the rotating crack model. In the fixed crack model a crack forms
perpendicular to the principal tensile stress direction when the principal stress exceeds the
concrete tensile strength and the crack orientation does not change during subsequent
loading. The ease of formulating and implementing this model has led to its wide-spread used
in early studies (Hand et al. 1973; Lin and Scordelis 1975). Subsequent studies, however,
showed that the model is associated with numerical problems caused by the singularity of the
material stiffness matrix. Moreover, the crack pattern predicted by the finite element analysis
often shows considerable deviations from that observed in experiments (Jain and Kennedy
1974).

The problems of the fixed crack model can be overcome by introducing a cracked
shear modulus, which eliminates most numerical difficulties of the model and considerably
improves the accuracy of the crack pattern predictions. The results do not seem to be very
sensitive to the value of the cracked shear modulus (Vebo and Ghali 1977; Barzegar and
Schnobrich 1986), as long as a value which is greater than zero is used, so as to eliminate the
CHAPTER 1 7

singularity of the material stiffness matrix and the associated numerical instability. Some
recent models use a variable cracked shear modulus to represent the change in shear stiffness,
as the principal stresses in the concrete vary from tension to compression (Balakrishnan and
Murray 1988; Cervenka et al. 1990).

de Borst and Nauta (1985) have proposed a model in which the total strain rate is
additively decomposed into a concrete strain rate and a crack strain rate. The latter is, in turn,
made up of several crack strain components. After formulating the two-dimensional concrete
stress-strain relation and transforming from the crack direction to the global coordinate
system of the structure, a material matrix with no coupling between normal and shear stress is
constructed. In spite of its relative simplicity and ease of application, this approach still
requires the selection of a cracked shear modulus of concrete.

In the rotating crack model proposed by Cope et al. (1980) the crack direction is not
fixed during the subsequent load history. Several tests by Vecchio and Collins (1982) have
shown that the crack orientation changes with loading history and that the response of the
specimen depends on the current rather than the original crack direction. In the rotating crack
model the crack direction is kept perpendicular to the direction of principal tensile strain and,
consequently, no shear strain occurs in the crack plane. This eliminates the need for a cracked
shear modulus. A disadvantage of this approach is the difficulty of correlating the analytical
results with experimental fracture mechanics research, which is at odds with the rotating
crack concept. This model has, nonetheless, been successfully used in analytical studies of
RC structures whose purpose is to study the global structural behavior, rather than the local
effects in the vicinity of a crack (Gupta and Akbar 1983; Adeghe and Collins 1986).

While the response of lightly reinforced beams in bending is very sensitive to the
effect of tension stiffening of concrete, the response of RC structures in which shear plays an
important role, such as over-reinforced beams and shear walls, is much more affected by the
bond-slip of reinforcing steel than the tension stiffening of concrete. To account for the bond-
slip of reinforcing steel two different approaches are common in the finite element analysis of
RC structures. The first approach makes use of the bond link element proposed by Ngo and
Scordelis (1967). This element connects a node of a concrete finite element with a node of an
adjacent steel element. The link element has no physical dimensions, i.e. the two connected
nodes have the same coordinates.

The second approach makes use of the bond-zone element developed by de Groot et
al. (1981). In this element the behavior of the contact surface between steel and concrete and
8 CHAPTER 1

of the concrete in the immediate vicinity of the reinforcing bar is described by a material law
which considers the special properties of the bond zone. The contact element provides a
continuous connection between reinforcing steel and concrete, if a linear or higher order
displacement field is used in the discretization scheme. A simpler but similar element was
proposed by Keuser and Mehlhorn (1987), who showed that the bond link element cannot
represent adequately the stiffness of the steel-concrete interface.

Even though many studies of the bond stress-slip relationship between reinforcing
steel and concrete have been conducted, considerable uncertainty about this complex
phenomenon still exists, because of the many parameters which are involved. As a result,
most finite element studies of RC structures do not account for bond-slip of reinforcing steel
and many researchers express the opinion that this effect is included in the tension-stiffening
model.

Very little work has been done, so far, on the three-dimensional behavior of
reinforced concrete systems using solid finite elements, because of the computational effort
involved and the lack of knowledge of the material behavior of concrete under three-
dimensional stress states. Suidan and Schnobrich (1973) were the first to study the behavior
of beams with 20-node three-dimensional isoparametric finite elements. The behavior of
concrete in compression was assumed elasto-plastic based on the von-Mises yield criterion. A
coarse finite element mesh was used in these analyses for cost reasons.

In spite of the large number of previous studies on the nonlinear finite element
analysis of reinforced concrete structures, only few conclusions of general applicability have
been arrived at. The inclusion of the effects of tension stiffening and bond-slip is a case in
point. Since few rational models of this difficult problem have been proposed so far, it is
rather impossible to assess exactly what aspects of the behavior are included in each study
and what the relative contribution of each is. Similar conclusions can be reached with regard
to other aspects of the finite element analysis. Even though the varying level of sophistication
of proposed models is often motivated by computational cost considerations, the multitude of
proposed approaches can lead to the conclusion that the skill and experience of the analyst is
the most important aspect of the study and that the selection of the appropriate model
depends on the problem to be solved.

Recognizing that many of the previously proposed models and methods have not been
fully verified so far, it is the intent of this study to address some of the model selection issues,
in particular, with regard to the effects of tension-stiffening and bond-slip.
CHAPTER 1 9

1.3 Objectives and Scope

The present investigation of the nonlinear response to failure of RC structures under

short term monotonic loads was initiated with the intent to investigate the relative importance
of several factors in the nonlinear finite element analysis of RC structures: these include the
effect of tension-stiffening and bond-slip and their relative importance on the response of
under- and over-reinforced beams and slabs, the effect of size of the finite element mesh on
the analytical results and the effect of the nonlinear behavior of concrete and steel on the
response of under- and over-reinforced beams and slabs. In the progress of this study
improved material models were developed and included in the analysis. These are presently
extended to study the behavior of RC structures under cyclic loads, which will be the subject
of a future report.

The main objectives of this study are:

• To develop improved analytical models for the study of the nonlinear behavior of RC
beams, slabs and beam-column joints under short term monotonic loads. These
models should be simple enough to allow extension to cyclic loading conditions
without undue computational cost.

• To investigate the relative importance of the nonlinear behavior of concrete,

reinforcing steel and their interaction through bond-slip on the response of under- and
over-reinforced concrete beams and slabs.

• To investigate the effect of size of the finite element mesh on the results of the
nonlinear analysis of RC structures and to develop an improved criterion for reducing
the numerical error associated with large finite element mesh size.
• To investigate the relative contribution of bond-slip and tension stiffening to the post-
cracking stiffness of under- and over-reinforced beams and beam-column joints.

• To investigate the effect of material and numerical parameters, such as concrete

tensile strength and shear retention factor, on the response of reinforced concrete
beams.

Following the introduction and a brief review of previous studies in Chapter 1

Chapter 2 deals with the description of the material behavior of concrete and reinforcing steel
and their interaction through bond-slip. The behavior of concrete and reinforcing steel is
modeled independently. A new discrete, embedded steel model is developed which accounts
10 CHAPTER 1

for the relative slip between reinforcing steel and concrete. The behavior of concrete under
biaxial loading conditions is described by a nonlinear orthotropic model, in which the axes of
orthotropy coincide with the principal strain directions (rotating crack model). The effect of
size of the finite element mesh is discussed in connection with a new smeared crack model
and an improved criterion is derived from fracture mechanics considerations in order to
reduce the numerical error associated with large size finite elements.

Chapter 3 discusses the inclusion of the material models in a finite element for plane
stress analysis of RC beams. A new plate bending element is also developed based on the
Mindlin plate theory and the layer concept. The element accounts for shear deformations and
is capable of modeling the gradual propagation of cracks through the depth of the slab.
Chapter 3 concludes with a discussion of computational aspects of the solution. The iteration
scheme and the resulting nonlinear solution procedure are described along with the associated
convergence criteria. A special purpose nonlinear algorithm is developed for the solution of
the embedded discrete bar model with bond-slip.

The validity of the proposed models is established by comparing the analytical

predictions with results from experimental and previous analytical studies in Chapter 4.
Finally, Chapter 5 presents the conclusions of this study.
 CHAPTER 2
MATERIAL MODELING

2.1 General

Reinforced concrete structures are made up of two materials with different

characteristics, namely, concrete and steel. Steel can be considered a homogeneous material
and its material properties are generally well defined. Concrete is, on the other hand, a
heterogeneous material made up of cement, mortar and aggregates. Its mechanical properties
scatter more widely and cannot be defined easily. For the convenience of analysis and design,
however, concrete is often considered a homogeneous material in the macroscopic sense.

III

II
Range I: Elastic
LOAD

Range II: Cracking

Range III: Steel Yielding or
Concrete Crushing

DEFLECTION

FIGURE 2.1 TYPICAL LOAD-DISPLACEMENT RESPONSE OF RC ELEMENT

11
12 CHAPTER 2

The typical stages in the load-deformation behavior of a reinforced concrete simply

supported beam are illustrated in Figure 2.1. Similar relations are obtained for other types of
reinforced concrete structural elements. This highly nonlinear response can be roughly
divided into three ranges of behavior: the uncracked elastic stage, the crack propagation and
the plastic (yielding or crushing) stage.

The nonlinear response is caused by two major effects, namely, cracking of concrete
in tension, and yielding of the reinforcement or crushing of concrete in compression.
Nonlinearities also arise from the interaction of the constituents of reinforced concrete, such
as bond-slip between reinforcing steel and surrounding concrete, aggregate interlock at a
crack and dowel action of the reinforcing steel crossing a crack. The time-dependent effects
of creep, shrinkage and temperature variation also contribute to the nonlinear behavior.
Furthermore, the stress-strain relation of concrete is not only nonlinear, but is different in
tension than in compression and the mechanical properties are dependent on concrete age at
loading and on environmental conditions, such as ambient temperature and humidity. The
material properties of concrete and steel are also strain-rate dependent to a different extent.

Because of these differences in short- and long-term behavior of the constituent

materials, a general purpose model of the short- and long-term response of RC members and
structures should be based on separate material models for reinforcing steel and concrete,
which are then combined along with models of the interaction between the two constituents
to describe the behavior of the composite reinforced concrete material. This is the approach
adopted in this study. The assumptions made in the description of material behavior are
summarized below:

• The stiffness of concrete and reinforcing steel is formulated separately. The results are
then superimposed to obtain the element stiffness;
• The smeared crack model is adopted in the description of the behavior of cracked
concrete;

• Cracking in more than one direction is represented by a system of orthogonal cracks;

• The crack direction changes with load history (rotating crack model);

• The reinforcing steel is assumed to carry stress along its axis only and the effect of
dowel action of reinforcement is neglected;
CHAPTER 2 13

• The transfer of stresses between reinforcing steel and concrete and the resulting bond-
slip is explicitly accounted for in a new discrete reinforcing steel model, which is
embedded in the concrete element.

In the following the behavior of each constituent material and the derivation of the
corresponding material stiffness matrix is discussed separately. This is followed by the model
of the interaction between reinforcing steel and concrete through bond. The superposition of
the individual material stiffness matrices to form the stiffness of the composite reinforced
concrete material and the numerical implementation of this approach in the nonlinear analysis
of beams and slabs is discussed in the next chapter.

2.2 Concrete

2.2.1 Behavior of Concrete

Concrete exhibits a large number of microcracks, especially, at the interface between

coarser aggregates and mortar, even before subjected to any load. The presence of these
microcracks has a great effect on the mechanical behavior of concrete, since their propagation
during loading contributes to the nonlinear behavior at low stress levels and causes volume
expansion near failure. Many of these microcracks are caused by segregation, shrinkage or
thermal expansion of the mortar. Some microcracks may develop during loading because of
the difference in stiffness between aggregates and mortar. Since the aggregate-mortar
interface has a significantly lower tensile strength than mortar, it constitutes the weakest link
in the composite system. This is the primary reason for the low tensile strength of concrete.
The response of a structure under load depends to a large extent on the stress-strain
relation of the constituent materials and the magnitude of stress. Since concrete is used
mostly in compression, the stress-strain relation in compression is of primary interest. Such a
relation can be obtained from cylinder tests with a height to diameter ratio of 2 or from strain
measurements in beams.

The concrete stress-strain relation exhibits nearly linear elastic response up to about
30% of the compressive strength. This is followed by gradual softening up to the concrete
compressive strength, when the material stiffness drops to zero. Beyond the compressive
strength the concrete stress-strain relation exhibits strain softening until failure takes place by
crushing.
14 CHAPTER 2

2.2.2 Concrete Models

Many mathematical models of the mechanical behavior of concrete are currently in

use in the analysis of reinforced concrete structures. These can be divided into four main
groups: orthotropic models, nonlinear elasticity models, plastic models and endochronic
models (Chen 1976, ASCE 1982, Meyer and Okamura, eds. 1985).

The orthotropic model is the simplest. It can match rather well experimental data
under proportional biaxial loading and approximates the concrete behavior under general
biaxial loading adequately. The model was also found capable of representing the hysteretic
behavior of concrete under cyclic loading (Darwin and Pecknold 1977). It is particularly
suitable for the analysis of reinforced concrete beams, panels and shells, since the stress state
of these structures is predominantly biaxial, and the model can be calibrated against an
extensive experimental data base. The equivalent uniaxial strain model was recently extended
to monotonic triaxial behavior (Chen 1976; ASCE 1982).

The nonlinear elasticity model is based on the concept of variable moduli and matches
well several available test data. In the pre-failure regime unique approximate relationships
have been established between hydrostatic and volumetric strain and between deviatoric
stress and strain. From these relationships expressions for the tangent bulk and shear modulus
can be derived. Thus, the nonlinear response of concrete is simulated by a piecewise linear
elastic model with variable moduli. The model is, therefore, computationally simple and is
particularly well suited for finite element calculations. When unloading takes place, the
behavior can be approximated by moduli which are different from those under loading
conditions. The model exhibits, however, continuity problems for stress paths near neutral
loading. As a result, the variable moduli model is unable to describe accurately the behavior
of concrete under high stress, near the compressive strength and in the strain softening range.

The plastic model, especially, the strain hardening plastic model can be considered as
a generalization of the previous models. The formulation of the constitutive relations in the
strain hardening plastic model is based on three fundamental assumptions: (1) the shape of
the initial yield surface; (2) the evolution of the loading surface, i.e. the hardening rule; and
(3) the formulation of an appropriate flow rule. Even though this model represents
successfully the behavior in the strain hardening region, the strain softening behavior of
concrete beyond the peak stress cannot be described adequately by the classical theory of
work-hardening plasticity which is based on Drucker's postulate of material stability. It is,
CHAPTER 2 15

therefore, not appropriate to use this model in the analysis of reinforced concrete structures
which experience strain softening. The model is, nevertheless, used extensively in the study
of concrete behavior, since the introduction of additional assumptions renders it capable of
simulating the behavior of concrete with sufficient accuracy (Chen 1976; Arnesen et al.
1980).

The endochronic theory of plasticity is based on the concept of intrinsic or

endochronic time. The intrinsic time is used to measure the extent of damage in the internal
structure of the concrete material under general deformation histories. Many features of
concrete behavior may be represented by this theory without need for loading-unloading
conditions. The introduction of loading criteria becomes, however, necessary for an accurate
material representation. The most logical way of accomplishing this is by introducing loading
surfaces and plasticity hardening rules. Recent applications have clearly demonstrated the
power of the endochronic approach (ASCE 1982). Further research, is, however, needed to
refine the theory, and to reduce the number of material constants which describe the material
behavior.

σ2

ft σ1
C
TI
AS
EL

0.6 fc
R
EA

feq
N
LI

INITIAL YIELD SURFACE

ORTHOTROPIC
fc
ULTIMATE LOAD NONLINEAR ELASTIC
SURFACE

FIGURE 2.2 STRENGTH FAILURE ENVELOPE OF CONCRETE

A very promising model has been recently proposed by Bazant and Ozbolt (1989).
The so-called microplane model seems to be capable of representing quite well several
16 CHAPTER 2

features of the monotonic and triaxial behavior of concrete and is thus, particularly, well
suited for local analyses of reinforced concrete structures. It is, however, very expensive
computationally and its use in the analysis of large scale structures does not appear feasible at
present.

The orthotropic model is adopted in this study for its simplicity and computational
efficiency. The ultimate objective of this work is the response analysis of large structures
under cyclic loads for which the endochronic and the microplane model are prohibitively
expensive.

Under combinations of biaxial compressive stress concrete exhibits strength and

stress-strain behavior which is different from that under uniaxial loading conditions. Figure
2.2 shows the biaxial strength envelope of concrete under proportional loading (Kupfer et al.
1969; Tasuji et al. 1978). Under biaxial compression concrete exhibits an increase in
compressive strength of up to 25% of the uniaxial compressive strength, when the stress ratio
61 62 is 0.5. When the stress ratio 61 62 is equal to 1 the concrete compressive strength is
approximately 1.16 ⋅ f c , where f c is the uniaxial concrete compressive strength. Under
biaxial tension concrete exhibits constant or perhaps slightly increased tensile strength
compared with that under uniaxial loading. Under a combination of tension and compression
the compressive strength decreases almost linearly with increasing principal tensile stress
(Figure 2.2).

The principal stress ratio also affects the stiffness and strain ductility of concrete.
Under biaxial compression concrete exhibits an increase in initial stiffness which may be
attributed to Poisson's effect. It also exhibits an increase in strain ductility indicating that less
internal damage takes place under biaxial compression than under uniaxial loading. More
details are presented by Kupfer et al. (1969), Chen (1976), and Tasuji et al. (1978). Even
though the principal compressive strain and the principal tensile strain at failure decrease
with increasing tensile stress in the compression-tension region, failure basically takes place
by cracking and the principal compressive stress and strain remain in the ascending branch of
the concrete stress-strain relation.

Since most of a beam or slab which is subjected to bending moments, experiences

biaxial stress combinations in the tension-tension or compression-compression region and
only small portions near the supports lie in the compression-tension region, the proposed
concrete model uses a different degree of approximation in each region of the biaxial stress
and strain space.
CHAPTER 2 17

The behavior of the model depends on the location of the present stress state in the
principal stress space in Figure 2.2. In the biaxial compression region the model remains
linear elastic for stress combinations inside the initial yield surface in Figure 2.2. Both the
initial yield and the ultimate load surface are described by the expression proposed by Kupfer
et al. (1969) (Figure 2.2)

( 61 + 62 )
2

F= − A ⋅ fc = 0 (2.1)
6 2 + 3.6561
where 61 and 62 are the principal stresses, f c is the uniaxial compressive strength and A is a
parameter. A=0.6 defines the initial yield surface, while A=1.0 defines the ultimate load
surface under biaxial compression.

For stress combinations outside the initial yield surface but inside the ultimate failure
surface the behavior of concrete is described by a nonlinear orthotropic model. This model
derives the biaxial stress-strain response from equivalent uniaxial stress-strain relations in the
axes of orthotropy.

σ σ

ε ε
(a) (b)
σ σ

ε ε
(c) (d)

FIGURE 2.3 IDEALIZATION OF STRESS-STRAIN CURVE OF CONCRETE

18 CHAPTER 2

In describing the uniaxial stress-strain behavior of concrete many empirical formulas

have been proposed. These are summarized by Popovics (1969) and ASCE (1982). Figure
2.3a shows the simplest of the nonlinear models, the linearly elastic-perfectly plastic model,
which was used by Lin and Scordelis (1975) in a study of reinforced concrete slabs and walls.
Figure 2.3b shows the inelastic-perfectly plastic model proposed by the European Concrete
Committee (CEB 1978) made up of a parabola and a horizontal line. The model proposed by
Hognestad (1951) is shown in Figure 2.3c. This model is capable of representing quite well
the uniaxial stress-strain behavior of a wide range of concretes. Finally, Figure 2.3d shows a
piecewise linear model in which the nonlinear stress-strain relation is approximated by a
series of straight line segments. Although this is the most versatile model capable of
representing a wide range of stress-strain curves, its use is restricted to cases in which
experimental data for the uniaxial concrete stress-strain relation are available.

σip
0.85 σip

0.6 σip

crushing

εo
ε
ε ip ε iu

f eq

FIGURE 2.4 STRESS-STRAIN RELATION OF CONCRETE

In the present study the model of Hognestad (1951) is used after some modifications.
These modifications are introduced in order to increase the computational efficiency of the
model, and in view of the fact that the response of typical reinforced concrete structures is
much more affected by the tensile than by the compressive behavior of concrete. This stems
CHAPTER 2 19

from the fact that the concrete tensile strength is generally less than 20% of the compressive
strength. In typical reinforced concrete beams and slabs which are subjected to bending, the
maximum compressive stress at failure does not reach the compressive strength. This means
that the compressive stresses in most of the member reach a small fraction of the compressive
strength at failure. The behavior of these members is, therefore, dominated by crack
formation and propagation, and the yielding of reinforcing steel.
The equivalent concrete compressive strength in each axis of orthotropy 6ip is
determined from the biaxial failure surface of concrete, where i is equal to 1 or 2. In order to
simplify the concrete material model the stress-strain relation in compression is assumed
piecewise linear with three branches. The material remains linear elastic up to a stress of
0.6 ⋅ 6ip (Figure 2.4), since at a stress between 50% to 70% of 6ip cracks at nearby aggregate
surfaces start to bridge in the form of mortar cracks and other bond cracks continue to grow
slowly. Beyond this stress the behavior is assumed linear up to the compressive strength 6ip
followed by a linear descending branch which represents strain softening. The strain .ip at
the compressive strength 6ip in Figure 2.4 is determined so as to match the strain energy in
compression of the experimental uniaxial stress-strain curve.

When the biaxial stresses exceed Kupfer's failure envelope, concrete enters into the
strain softening range of behavior where an orthotropic model describes the biaxial behavior
(Darwin and Pecknold 1977). In this region failure occurs by crushing of concrete when the
principal compressive strain exceeds the limit value .iu in Figure 2.4. In order to define the
crushing of concrete under biaxial compressive strains a strain failure surface is used as
ultimate failure criterion. This surface, which is shown in Figure 2.5, is described by the
following equation in complete analogy to Kupfer's stress failure envelope in Eq. 2.1:

( .1 + . 2 )
2

C= − . cu = 0 (2.2)
. 2 + 3.65.1
where .1 and . 2 are the principal strains and . cu is the ultimate strain of concrete in uniaxial
loading. For strain combinations outside the strain failure surface the element is assumed to
lose its strength completely and is not able to carry any more stress.

The equations which define the parameters of the equivalent uniaxial stress-strain
relation in the main axes of orthotropy can now be summarized for the compression-
compression region of the principal stress space. Denoting the principal stress ratio by α
20 CHAPTER 2

61
*= where 61 ≥ 62 (2.3)
62
yields the following equations
1 + 3.65*
62 p = ⋅ fc (2.4)
(1 + * ) 2

æ 6 ö
. 2 p = . co ⋅ ç 3 2 p − 2 ÷ (2.5)
è fc ø
61 p = * ⋅ 6 2 p (2.6)

é æ 61 p ö
3
æ 61 p ö
2
æ 61 p ö ù
.1 p = .co ⋅ ê −1.6 ⋅ ç ÷ + 2.25 ⋅ ç ÷ + 0.35 ⋅ ç ÷ú (2.7)
êë è fc ø è fc ø è f c ø úû
where . co is the strain corresponding to the concrete compressive strength f c under uniaxial
stress conditions.

ε2

ε1
εo

CRACKED

ε cu

CRUSHED

FIGURE 2.5 STRAIN FAILURE ENVELOPE OF CONCRETE

In the biaxial compression-tension and tension-tension region the following

assumptions are adopted in this study: (1) failure takes place by cracking and, therefore, the
tensile behavior of concrete dominates the response; (2) the uniaxial tensile strength of
CHAPTER 2 21

concrete is reduced to the value f eq , as shown in Figs. 2.2 and 2.4 to account for the effect of
the compressive stress; in the tension-tension region the tensile strength remains equal to the
uniaxial tensile strength f t (Figure 2.2); (3) the concrete stress-strain relation in compression
is the same as under uniaxial loading and does not change with increasing principal tensile
stress. Thus 6ip in Figure 2.4 is equal to f c and .iu is equal to . cu . The last assumption holds
true in the range of compressive stresses which is of practical interest in typically reinforced
beams and slabs.
The proposed model assumes that concrete is linear elastic in the compression-tension
and the biaxial tension region for tensile stresses smaller than f eq in Figure 2.4. Beyond the
tensile strength the tensile stress decreases linearly with increasing principal tensile strain
(Figure 2.4). Ultimate failure in the compression-tension and the tension-tension region takes
place by cracking, when the principal tensile strain exceeds the value .o in Figure 2.4. This
results in the strain failure surface of Figure 2.5. The value of .o is derived from fracture
mechanics concepts, as will be described in the following section. When the principal tensile
strain exceeds .o , the material only loses its tensile strength normal to the crack, while it is
assumed to retain its strength parallel to the crack direction.

2.2.3 Behavior of Cracked Concrete

2.2.3.1 Description of a Cracked Section

The nonlinear response of concrete is often dominated by progressive cracking which

results in localized failure. Figure 2.6a depicts part of a reinforced concrete member in
flexure. The member has cracked at discrete locations where the concrete tensile strength is
exceeded.
At the cracked section all tension is carried by the steel reinforcement. Tensile
stresses are, however, present in the concrete between the cracks, since some tension is
transferred from steel to concrete through bond. The magnitude and distribution of bond
stresses between the cracks determines the distribution of tensile stresses in the concrete and
the reinforcing steel between the cracks.
22 CHAPTER 2

Additional cracks can form between the initial cracks, if the tensile stress exceeds the
concrete tensile strength between previously formed cracks. The final cracking state is
reached when a tensile force of sufficient magnitude to form an additional crack between two

cracks (a)
M1 M2

M1 M
M2 (b)

u
(c)

ft (d)

fs
(e)

EI (f)

FIGURE 2.6 EFFECT OF CRACKING IN REINFORCED CONCRETE BEAM

(A) PORTION OF BEAM
(B) BENDING MOMENT DISTRIBUTION
(C) BOND STRESS DISTRIBUTION
(D) CONCRETE TENSILE STRESS DISTRIBUTION
(E) STEEL TENSILE STRESS DISTRIBUTION
(F) FLEXURAL STIFFNESS DISTRIBUTION IN ELASTIC RANGE
CHAPTER 2 23

existing cracks can no longer be transferred by bond from steel to concrete. Figs. 2.6c, 2.6d
and 2.6e show the idealized distribution between cracks of bond stress, concrete tensile stress
and steel tensile stress, respectively. Because concrete is carrying some tension between the
cracks, the flexural rigidity is clearly greater between the cracks than at the cracks, as shown
in Figure 2.6f (Park and Paulay 1975).
In order to improve the accuracy of finite element models in representing cracks and,
in some cases, in order to improve the numerical stability of the solution the tension
stiffening effect was introduced in several models. The physical behavior in the vicinity of a
crack can be inferred from Figure 2.6d and Figure 2.6e. As the concrete reaches its tensile
strength, primary cracks form. The number and the extent of cracks are controlled by the size
and placement of the reinforcing steel. At the primary cracks the concrete stress drops to zero
and the steel carries the entire tensile force. The concrete between the cracks, however, still
carries some tensile stress, which decreases with increasing load magnitude. This drop in
concrete tensile stress with increasing load is associated with the breakdown of bond between
reinforcing steel and concrete. At this stage a secondary system of internal cracks, called
bond cracks, develops around the reinforcing steel, which begins to slip relative to the
surrounding concrete.
Since cracking is the major source of material nonlinearity in the serviceability range
of reinforced concrete structures, realistic cracking models need to be developed in order to
accurately predict the load-deformation behavior of reinforced concrete members. The
selection of a cracking model depends on the purpose of the finite element analysis. If overall
load-deflection behavior is of primary interest, without much concern for crack patterns and
estimation of local stresses, the "smeared" crack model is probably the best choice. If detailed
local behavior is of interest, the adoption of a "discrete" crack model might be necessary.
Unless special connecting elements and double nodes are introduced in the finite
element discretization of the structure, the well established smeared crack model results in
perfect bond between steel and concrete, because of the inherent continuity of the
displacement field. In this case the steel stress at the cracks will be underestimated.
One way of including the tension stiffening effect in the smeared crack model is to
increase the average stiffness of the finite element which contains the crack. Considering that
the finite element has relatively large dimensions compared to the size of the cracked section
two methods have been proposed. In the first method the tension portion of the concrete
stress-strain curve is assigned a descending branch. In this case the tension stiffening effect is
24 CHAPTER 2

represented as either a step-wise reduction of tensile stress or as a gradually unloading model

(Lin and Scordelis 1975; Gilbert and Warner 1977; Barzegar and Schnobrich 1986). Even in
terms of plain concrete this is not completely unrealistic, since the stiffness of the test loading
apparatus affects the strain softening stiffness in tension and disagreement exists on the
ability of concrete to carry stress beyond the ultimate tensile strength. In the second tension
stiffening model the steel stiffness is increased (Gilbert and Warner 1977, Cervenka et al.
1990). The additional stress in the reinforcing steel represents the tensile force carried by the
concrete between the cracks. For reasons of computational convenience it is assumed that the
orientation of this additional stress coincides with the orientation of the reinforcing steel.

2.2.3.2 Crack Models

The first reinforced concrete finite element model which includes the effect of
cracking was developed by Ngo and Scordelis (1967), who carried out a linear elastic
analysis of beams with predefined crack patterns. The cracks were modeled by separating the
nodal points of the finite element mesh and thus creating a discrete crack model (Figure
2.7a). With the change of topology and the redefinition of nodal points the narrow band width
of the stiffness matrix is destroyed and a greatly increased computational effort results in this
model. Moreover, the lack of generality in crack orientation has made the discrete crack
model unpopular. In spite of these shortcomings, the use of discrete crack models in finite
element analysis offers certain advantages over other methods. For those problems that
involve a few dominant cracks, the discrete crack approach offers a more realistic description
of the cracks, which represent strain discontinuities in the structure. Such discontinuities are
correctly characterized by the discrete crack model.
The need for a crack model that offers automatic generation of cracks and complete
generality in crack orientation, without the need of redefining the finite element topology, has
led the majority of investigators to adopt the smeared crack model. Rather than representing a
single crack, as shown in Figure 2.7a, the smeared crack model represents many finely spaced
cracks perpendicular to the principal stress direction, as illustrated in Figure 2.7b. This
approximation of cracking behavior of concrete is quite realistic, since the fracture behavior
of concrete is very different from that of metals. In concrete fracture is preceded by
microcracking of material in the fracture process zone, which manifests itself as strain
softening. This zone is often very large relative to the cross section of the member due to the
large size of aggregate (Figure 2.8a). In a steel member fracture is preceded by yielding of
material in the process zone which is concentrated near the crack tip and has a relatively
CHAPTER 2 25

small size (Figure 2.8b). In this case a discrete crack model is a more realistic representation
of actual behavior.

node

element

crack

(a) (b)

FIGURE 2.7 CRACKING MODELS:

(A) DISCRETE CRACK,
(B) SMEARED CRACK

The smeared crack model first used by Rashid (1968) represents cracked concrete as
an elastic orthotropic material with reduced elastic modulus in the direction normal to the
crack plane. With this continuum approach the local displacement discontinuities at cracks
are distributed over some tributary area within the finite element and the behavior of cracked
concrete can be represented by average stress-strain relations. In contrast to the discrete crack
concept, the smeared crack concept fits the nature of the finite element displacement method,
since the continuity of the displacement field remains intact.
Although this approach is simple to implement and is, therefore, widely used, it has
nevertheless a major drawback, which is the dependency of the results on the size of the finite
element mesh used in the analysis (Vebo and Ghali 1977; Bazant and Cedolin 1980). When
large finite elements are used, each element has a large effect on the structural stiffness.
When a single element cracks, the stiffness of the entire structure is greatly reduced. Higher
order elements in which the material behavior is established at a number of integration points
do not markedly change this situation, because, in most cases, when a crack takes place at
one integration point, the element stiffness is reduced enough so that a crack will occur at all
other integration points of the element in the next iteration. Thus, a crack at an integration
point does not relieve the rest of the material in the element, since the imposed strain
continuity increases the strains at all other integration points. The overall effect is that the
formation of a crack in a large element results in the softening of a large portion of the
26 CHAPTER 2

structure. The difficulty stems from the fact that a crack represents a strain discontinuity
which cannot be modeled correctly within a single finite element in which the strain varies
continuously. Many research efforts have been devoted to the solution of this problem based,
in particular, on fracture mechanics concepts (Hillerborg et al. 1976; Bazant and Cedolin
1980, 1983).

concrete steel

F Y

(a) (b)

FIGURE 2.8 RELATIVE SIZE OF CRACK PROCESS ZONE

(A) FRACTURE ZONE (F) OF CONCRETE
(B) YIELDING ZONE (Y) OF STEEL

The success fracture mechanics theory (Broek 1974) had in solving different types of
cracking problems in metals, ceramics and rocks has lead to its use in the finite element
analysis of reinforced concrete structures. If it is accepted that concrete is a notch-sensitive
material, it can be assumed that a cracking criterion which is based on tensile strength may be
dangerously unconservative and only fracture mechanics theory provides a more rational
approach to the solution of the problem. In its current state of development, however, the
practical applicability of fracture mechanics to reinforced concrete is still in question and
much remains to be done. Intensive research in this area is presently undertaken by several
researchers (Hillerborg et al. 1976; Bazant and Cedolin 1980, 1983; Jenq and Shah 1986).
In order to define the strain softening branch of the tensile stress-strain relation of
concrete by fracture mechanics concepts three important parameters need to be defined: (1)
the tensile strength of concrete at which a fracture zone initiates; (2) the area under the stress-
strain curve; and (3) the shape of the descending branch (Reinhardt 1986). Among these
parameters, the first two can be considered as material constants, while the shape of the
descending branch varies in the models that have been proposed (Bazant and Oh 1983).
CHAPTER 2 27

Before discussing two of the most prominent models, a relation between the area under the
tensile stress-crack strain diagram in Figure 2.9a and the fracture energy G f is needed. This
relation can be readily derived by the following procedure.

σ nn

smeared
ft

gf

ε cr
nn
(a)
σnn

discrete
ft

Gf

w
(b)

FIGURE 2.9 STRAIN SOFTENING BEHAVIOR OF CONCRETE

(A) TENSILE STRESS VS. CRACK STRAIN RELATION
(B) TENSILE STRESS VS. CRACK OPENING DISPLACEMENT RELATION

The area g f under the curve in Figure 2.9a can be expressed as:

g f = ò 6nn ⋅ d .nn
cr
(2.8)
The fracture energy G f is defined as the amount of energy required to crack one unit
of area of a continuous crack and is considered a material property. This definition results in
the following expression for the fracture energy G f

G f = ò 6 nn ⋅ dw (2.9)
28 CHAPTER 2

where w represents the sum of the opening displacements of all microcracks within the
fracture zone. Eq. 2.9 is schematically shown in Figure 2.9b.
In the smeared crack model w is represented by a crack strain which acts over a
certain width within the finite element called the crack band width b. Since w is the
accumulated crack strain, this is represented by the following relation

w = ò .crnn ⋅ dn (2.10)

Assuming that the microcracks are uniformly distributed across the crack band width,
Eq. 2.10 reduces to:

w = b ⋅ .crnn (2.11)
The combination of Eq. 2.11 with Eqs. 2.8 and 2.9 yields the relation between G f and g f :

Gf = b ⋅ g f (2.12)

The simplicity of Eq. 2.12 is misleading, because the actual size of the crack band
width b depends on the selected element size, the element type, the element shape, the
integration scheme and the problem type to be solved.

σ
σ

ft ft

gf gf

1
- f
3 t

ε ε
0 εo 0 2/9 εo εo
(a) (b)

FIGURE 2.10 TYPICAL SHAPES OF SOFTENING BRANCH OF CONCRETE

(A) BAZANT AND OH MODEL
(B) HILLERBORG MODEL

Two widely used models of the strain softening behavior of concrete in tension are
those of Bazant and Hillerborg. Bazant and Oh (1983) introduced the "crack band theory" in
CHAPTER 2 29

the analysis of plain concrete panels. This model is one of the simplest fictitious crack
models. The two basic assumptions of the model are that the width of the fracture zone wc is
equal to three times the maximum aggregate size (approximately 1 inch) and that the concrete
strains are uniform within the band. In this case the final equation for determining the tensile
fracture strain .o takes the form (Figure 2.10a)

2G f
.o = (2.13)
ft ⋅ b
where b is the element width and G f is the fracture energy required to form a crack.

After an extensive experimental study Hillerborg et. al. (1976) proposed a bilinear
descending branch for the tensile strain softening behavior of concrete (Figure 2.10b). Using
the assumption that the microcracks are uniformly distributed over the crack band width and
combining the area g f with the fracture energy G f according to Eq. 2.12 the following
equation is derived for the tensile fracture strain

18 G f
.o = (2.14)
5 ft ⋅ b
Both models have been extensively used in the analysis of RC members and yield
very satisfactory results when the size of the finite element mesh is relatively small. The
analytical results, however, differ significantly from the experimental data when the finite
element mesh size becomes very large. This happens because both models assume a uniform
distribution of microcracks over a significant portion of a relatively large finite element while
the actual microcracks are concentrated in a much smaller cracked region of the element.
Thus Eqs. 2.13 and 2.14 cannot be directly applied to the numerical analysis of RC structures
with relatively large finite elements.

2.2.3.3 Proposed Model

Fracture and crack propagation in concrete depends to a large extent on the material
properties in tension and the post-cracking behavior. Experimental studies (Welch and
Haisman 1969; Bedard and Kotsovos 1986) indicate that the behavior of concrete after
cracking is not completely brittle and that the cracked region exhibits some ductility. As the
applied loads are increased the tensile stress in the critical cross section of the member
reaches the tensile strength f t . At this stage microcracks develop and form a fracture zone.
This process is characterized by the strain softening behavior of the section which ends when
30 CHAPTER 2

the microcracks coalesce to form one continuous macrocrack and stresses in the section
reduce to zero.
In order to account for the fact that microcracks are concentrated in a fracture process
zone which may be small compared to the size of the finite element mesh a distribution
function for the microcracks across the element width is introduced in this study. The
distribution function is exponential, so that it can represent the concentration of microcracks
near the crack tip when the finite element mesh size becomes fairly large (Figure 2.11).

f ( x ) = * ⋅ e +x (2.15)
in which * and + are constants to be determined.

1 b=3in
f(x)

b=6in

h=1 b=9in

b=12in
3
b b=92in
b
x 0 -
2
b

FIGURE 2.11 ASSUMED DISTRIBUTION OF MICROCRACKS IN AN ELEMENT

Using the boundary conditions that f (0) = 1.0 and f (b 2) = 3 b into Eq. 2.10 yields
the following equation for the distribution function

f ( x ) = e −2 b⋅ln b 3⋅ x (2.16)
where b is the element width. The condition that f (b 2) = 3 b ensures that, when the finite
element mesh size is equal to 3 inches, i.e. three times the maximum aggregate size (Bazant
and Oh 1983), the proposed distribution function reduces to Eq. 2.13 of the crack band
theory. As shown in Figure 2.12 the microcrack distribution f ( x ) is uniform for a finite
element mesh size smaller than 3 inches.
The fracture energy G f is defined as the product of the area under the equivalent
uniaxial stress-strain curve g f and the fracture zone. It can, therefore, be expressed as:
CHAPTER 2 31

b2
1
G f = .o ⋅ f t ⋅ 2 ⋅
2 ò
0
f ( x ) ⋅ dx (2.17)

where f t is the tensile strength of concrete, .o is the fracture tensile strain which
characterizes the end of the strain softening process when the microcracks coalesce into a
continuous crack and G f is the fracture energy which is dissipated in the formation of a
crack of unit length per unit thickness and is considered a material property.

0
-2b 0
-2b 0
-b2
(a) (b) (c)

FIGURE 2.12 DISTRIBUTION FUNCTION OF MICROCRACKS

RELATIVE TO FINITE ELEMENT MESH SIZE
(A) SMALL FINITE ELEMENT MESH SIZE
(B) AVERAGE FINITE ELEMENT MESH SIZE
(C) LARGE FINITE ELEMENT MESH SIZE

The experimental study by Welch and Haisman (1969) indicates that for normal
strength concrete the value of G f f t is in the range of 0.005-0.01 mm. If G f and f t are
known from measurements, then .o can be determined from

Gf
.o = b2 (2.18)
ft ⋅ ò
0
f ( x ) ⋅ dx

After substituting the function f ( x ) from Eq. 2.16 and integrating the following relation
results

2G f ⋅ ln ( 3 b )
Fo = (2.19)
f t ⋅ (3 − b)
which clearly shows that .o depends on the finite element mesh size. This approach of
defining .o renders the analytical solution insensitive to the mesh size and guarantees the
32 CHAPTER 2

objectivity of the results. At the same time this approach allows for the realistic
representation of the microcrack concentration near the tip of the crack in the case of large
finite elements, in which case Eq. 2.13 yields unsatisfactory results. With this approach large
finite elements can be used in the modeling of RC structures without loss of accuracy.

2.2.3.4 Concrete Material Matrix

In the analysis of reinforced concrete structures plane stress problems make up a large
majority of practical cases. In the following the constitutive relation for plane stress problems
is derived for the concrete model of this study.
For stress combinations inside the initial yield surface in Figure 2.2 concrete is
assumed to be a homogeneous, linear isotropic material. Thus, the stress-strain relation for
plane stress problems has the simple form

é ù
ì Tx ü ê1 O 0 ú ì Fx ü
ï ï E ê ú ï ï
í Ty ý = ⋅ êO 1 0 ú ⋅ í Fy ý (2.20)
ïU ï 1 − O
2

î xy þ
ê 1 − O ú ïî H xy ïþ
ê 0 0 ú
ë 2 û
where E is the initial elastic modulus of concrete and ν is Poisson's ratio.
Once the biaxial stress combination exceeds the initial yield surface in the
compression-compression region of Figure 2.2 concrete is assumed to behave as an
orthotropic material. The same assumption holds for stress combinations outside the ultimate
loading surface indicating that the material has entered into the strain softening range. With
reference to the principal axes of orthotropy the incremental constitutive relationship can be
expressed (Desai and Siriwardance 1972; Chen 1976)

ì d T11 ü é E1 O E1E2 0 ù ì dF ü
ï ï 1 ê ú ï 11 ï
íd T22 ý = ⋅ ê O E1E2 E2 0 ú ⋅ í d F 22 ý (2.21)
ï d U ï (1 − O ) ê 0
2

î 12 þ 0 (1 − O )G ú ïî d H12 ïþ
2
êë úû
where E1 and E2 are the secant moduli of elasticity in the direction of the axes of orthotropy,
which are oriented perpendicular and parallel to the crack direction, ν is Poisson's ratio and
(
(1 − O2 )G = 0.25 ⋅ E1 + E2 − 2O E1E2 .)
CHAPTER 2 33

The main advantage of this model is its simplicity and ease of calibration with
uniaxial concrete test data. Since the model is derived from the observation that the biaxial
stress-strain relation is not very sensitive to the compression stress ratio, it is mainly
applicable to plane stress problems such as beams, panels and thin shells where the stress
field is predominantly biaxial. By contrast, it is well known from experimental evidence that
hydrostatic pressure markedly influences the behavior of concrete under three-dimensional
stress states.
The use of the incremental orthotropic model under general stress histories, in which
the principal stress directions rotate during the loading process, has met with strong criticism,
both, on physical and theoretical grounds. This aspect of the model does not, however, appear
to affect its practical usefulness. The model was found capable of modeling satisfactorily the
concrete behavior under cyclic as well as monotonic loads. This model and its adaptations
have been applied to a wide variety of practical finite element problems with quite good
agreement between theoretical and experimental results in most cases.
When the principal tensile strain exceeds .o in Figure 2.4 a crack forms in a direction
perpendicular to the principal strain. The stress normal to the crack is zero and the shear
modulus needs to be reduced to account for the effect of cracking on shear transfer. In this
case the incremental constitutive relation takes the form

ì d T11 ü é E1 0 0 ù ì d F11 ü
ï ï 1 ê ú ⋅ ïdF ï
íd T22 ý = ⋅ 0 0
ê
0
ú í 22 ý
(2.22)
ï d U ï (1 − O ) ê 0 0 M (1 − O2 )G ú ïd H ï
2

î 12 þ ë û î 12 þ
where 1 and 2 are the directions parallel and perpendicular to the crack, respectively, and λ is
a cracked shear constant.
The use of the cracked shear constant λ not only solves most of the numerical
difficulties associated with a singular stiffness matrix, but also improves the representation of
the concrete cracking phenomenon in finite element analysis. The shear factor may also be
used as a way of suppressing the resulting singularity when all elements surrounding a
particular node crack in the same direction. The value of λ has a lower bound, which is
greater than zero and depends on the type of structure, the type of load and the accuracy of
the numerical representation. The effect of dowel action of the reinforcement and the
aggregate interlock tend to make the determination of an effective shear modulus rather
complex. According to previous studies (ASCE 1982) a crack shear constant value of 0.5 was
34 CHAPTER 2

used in the analysis of shear panels and deep coupling beams, a value of 0.25 was used in the
analysis of deep beams and a value of 0.125 in the analysis of shear walls and shear-wall
frame systems. Hand, Pecknold and Schnobrich (1973) used a value of 0.4 in the analysis of
RC plates and shells, Lin and Scordelis (1975) used a fixed value of 1.0 and Gilbert and
Warner (1977) used a fixed λ value of 0.6 in their analysis of RC slabs. From these and other
studies (Vebo and Ghali; Bashur and Darwin 1978; Barzegar and Schnobrich 1986) it appears
that the value of λ is not critical to the accuracy of the final results. This fact is also
corroborated by sensitivity analyses which are presented in Chapter 4, where it is shown that
the value of the cracked shear constant does not affect the response of beams in bending. The
cracked shear value λ in Eq. 2.22 is, consequently, fixed in the present study at a value of 0.4.
The use of the orthotropic constitutive relation in Eq. 2.21 to represent cracked
concrete may not be totally realistic. In the case of a real crack the crack surface is rough and
any sliding parallel to the crack will generate some local stresses or movement normal to the
crack. To properly represent this type of behavior the off-diagonal terms of the material
matrix which relate shear strain with normal stress should not be zero. The relative
magnitude of these off-diagonal terms decreases as the crack widens. However, this effect
may not be significant in a study which focuses attention on overall member behavior and
most researchers have neglected it (Lin and Scordelis 1975; Bashur and Darwin 1978).
The proposed concrete model accounts for progressive cracking and changes in the
crack direction by assuming that the crack is always normal to the principal strain direction.
In contrast to the model used by Hand et. al. (1973) and Lin and Scordelis (1975), the
material axes are not fixed after formation of the initial crack, but their orientation is
determined from the direction of principal strains at the beginning of each iteration.
In developing a numerical algorithm for the rotating crack concept Gupta and Akbar
(1983) obtained the rotating crack material matrix as the sum of the conventional fixed crack
material matrix in Eq. 2.22 and a contribution which reflects the change in crack direction.
This is expressed by the following equation

[ DLO ]*c = [ DLO ]c + [G ] (2.23)

where [ DLO ]c is given by Eq. 2.22 and [G ] reflects the change in crack direction and is given
by
CHAPTER 2 35

ésin 2 2+ − sin 2 2+ − sin 2+ cos 2+ ù

T1 ê ú
[G ] = − ⋅ê . sin 2 2+ sin 2+ cos 2+ ú (2.24)
2 F 2x − F 2x + H 2xy ê sym
ë . cos2 2+ úû

While Eq. 2.24 is theoretically correct within the assumptions of the rotating crack model,
any suitable incremental material stiffness matrix can be used in the context of an iterative
nonlinear solution algorithm. It is, therefore, possible to neglect the rotating crack
contribution [G ] provided that the change in crack direction is accounted for in the
orientation of the material axes and in the transformation from material to element coordinate
axes. Milford and Schnobrich (1984) found that neglecting the rotating crack contribution
[G ] in Eq. 2.23 only rarely increased the number of iterations and did not introduce any
numerical instabilities.
Since the material matrix is defined with reference to crack (principal strain)
directions, it must be transformed to the global coordinate system before all element stiffness
matrices can be assembled. This is accomplished by the following transformation.

[ DGL ]c = [T ]T ⋅ [ DLO ]c ⋅ [T ] (2.25)

where

é cos 2 + sin 2 + sin + cos + ù

ê ú
[T ] = ê sin + 2
cos +2
− sin + cos + ú (2.26)
ê −2sin + cos + 2sin + cos + cos 2 + − sin 2 + ú
ë û
and
H xy
tan 2+ = (2.27)
Fx − F y

+ is the angle between the direction normal to the crack and the global x-axis.

During the state determination process the location of the present state in the biaxial principal
stress and strain space needs to be determined. In this phase of the solution process the
principal strains and stresses are obtained from the strains and stresses in the global
coordinate system by the following relationships (Chia 1978):
36 CHAPTER 2

ì Fx ü
ì F1 ü é cos2 + sin 2 + sin + cos + ù ï ï
í ý=ê 2 ú ⋅ í Fy ý (2.28)
î F 2 þ ë sin + cos + − sin + cos + û ï H ï
2

î xy þ
ì Tx ü
ì T1 ü é cos2 + sin 2 + 2 sin + cos + ù ï ï
í ý=ê 2 ú ⋅ íTy ý (2.29)
îT 2 þ ë sin + cos + −2 sin + cos + û ïU ï
2

î xy þ

2.3 Reinforcing Steel

2.3.1 Behavior of Reinforcing Steel

The properties of reinforcing steel, unlike concrete, are generally not dependent on
environmental conditions or time. Thus, the specification of a single stress-strain relation is
sufficient to define the material properties needed in the analysis of reinforced concrete
structures.

σ
Es 2

σy 1

Es1

1
-ε u
ε
0 ε
u

-σ
y

FIGURE 2.13 STEEL STRESS-STRAIN RELATION

Typical stress-strain curves for reinforcing steel bars used in concrete construction are
obtained from coupon tests of bars loaded monotonically in tension. For all practical
purposes steel exhibits the same stress-strain curve in compression as in tension. The steel
CHAPTER 2 37

stress-strain relation exhibits an initial linear elastic portion, a yield plateau, a strain-
hardening range in which stress again increases with strain and, finally, a range in which the
stress drops off until fracture occurs. The extent of the yield plateau is a function of the
tensile strength of steel. High-strength, high-carbon steels, generally, have a much shorter
yield plateau than relatively low-strength, low-carbon steels.
Since steel reinforcement is used in concrete construction in the form of reinforcing
bars or wire, it is not necessary to introduce the complexities of three-dimensional
constitutive relations for steel. For computational convenience it even often suffices to
idealize the one dimensional stress-strain relation for steel. Three different idealizations,
shown in Fig. 2.13, are commonly used depending on the desired level of accuracy (ASCE
1982; Cervenka et al. 1990).
The first idealization neglects the strength increase due to strain hardening and the
reinforcing steel is modeled as a linear, perfectly plastic material, as shown in 2.13a. This
assumption underlies the design equations of the ACI code. If the strain at the onset of strain
hardening is much larger than the yield strain, this approximation yields very satisfactory
results. This is the case for low-carbon steels with low yield strength.
If the steel hardens soon after the onset of yielding, this approximation underestimates
the steel stress at high strains. In several instances it is necessary to evaluate the steel stress at
strains higher than yield to more accurately assess the strength of members at large
deformations. This is, particularly, true in seismic design, where assessing the available
ductility of a member requires that the behavior be investigated under strains many times the
yield strain. In this case more accurate idealizations which account for the strain hardening
effect are required, as shown in Figs. 2.13b and 2.13c. The parameters of these models are the
stress and strain at the onset of yielding, the strain at the onset of strain hardening and the
stress and strain at ultimate. These parameters can be derived from experimentally obtained
stress-strain relations.
In this study the reinforcing steel is modeled as a linear elastic, linear strain hardening
material with yield stress 6 y , as shown in Figure 2.13. The reasons for this approximation
are: (1) the computational convenience of the model; (2) the behavior of RC members is
greatly affected by the yielding of reinforcing steel when the structure is subjected to
monotonic bending moments. Yielding is accompanied by a sudden increase in the
deformation of the member. In this case the use of the elastic-perfectly plastic model in 2.13a
leads to numerical convergence problems near the ultimate member strength. It is, therefore,
38 CHAPTER 2

advisable to take advantage of the strain-hardening behavior of steel in improving the

numerical stability of the solution. The assumption of a linear strain hardening behavior
immediately after yielding of the reinforcement does not adversely affect the accuracy of the
results, as long as the slope of the strain hardening branch is determined so that the strain
energy of the model is equal to the strain energy of the experimental steel stress-strain
relation (Figure 2.13). Such a model has been successfully used in many analyses of
reinforced concrete structures (Ngo and Scordelis 1967; Vebo and Ghali 1977; Bashur and
Darwin 1978).

2.3.2 Steel Material Matrix

Three different approaches can be used in the modeling of reinforcing steel in a

reinforced concrete structure: (1) the distributed or smeared model; (2) the embedded model;
and (3) the discrete model (ASCE 1982).
In the distributed model the reinforcing steel is assumed to be distributed over the
concrete element at a certain orientation angle + . A composite concrete-steel constitutive
relation is used in this case, which, however, requires that perfect bond be assumed between
concrete and steel.
An embedded steel model is useful in connection with higher order isoparametric
concrete elements. The reinforcing steel is considered as an axial member built into the
isoparametric element such that its displacements are consistent with those of the element.
Such a model again implies perfect bond between concrete and steel.
The most widely used model represents the reinforcement with discrete one-
dimensional truss elements, which are assumed to be pin connected and possess two degrees
of freedom at each node. Alternatively, beam elements can also be used, in which case, three
degrees of freedom are allowed at each end of the bar element. In either case the one-
dimensional reinforcing bar elements can be easily superimposed on the two-dimensional
concrete element mesh. A significant advantage of the discrete representation, in addition to
its simplicity, is that it can include the slip of reinforcing steel with respect to the surrounding
concrete (Cervenka et al. 1990).
In this study a new discrete reinforcing steel model is used in the analysis of plane
stress problems. The model is embedded in the concrete element, but includes the effect of
bond-slip. In the analysis of slabs, where the effect of bond-slip is presently excluded from
CHAPTER 2 39

the model, it is computationally convenient to use the distributed steel model, in which the
reinforcing mesh can be represented by an equivalent steel layer, as will be described later.
In the proposed discrete model the reinforcing steel is represented by a one-
dimensional truss element which is embedded in the concrete element, as shown in Figure
2.14. The nodes of the steel element do not need to coincide with the nodes of the concrete
element. In this section the end displacements of the steel element are assumed to be
compatible with the boundary displacements of the concrete element, so that perfect bond is
implied. The discrete embedded steel element is then extended to include the effect of bond-
slip between reinforcing steel and concrete in Section 2.5.
Even though the end displacements of the steel element are compatible with the
concrete displacements, displacement compatibility within the concrete element is, generally,
not satisfied when the two node steel element is embedded in the 8-node serendipity element.
This fact, however, does not seem to affect the accuracy of results of the global behavior of
RC structures.
For the one-dimensional truss element with constant strain the stiffness matrix is
given by

ì P1 ü AE é 1 −1ù ì d1 ü
í ý= ⋅í ý (2.30)
î P2 þ L êë −1 1 úû î d 2 þ

or

{P} = [ K LO ]c ⋅ {d } (2.31)

where E is the modulus of elasticity, A the cross-sectional area and L the length of the bar
element. P1 and P2 are axial end forces and d1 and d 2 are axial end displacements of the
reinforcing bar. E is equal to Es1 before yielding of the reinforcement and equal to Es 2 after
yielding (Figure 2.13). denotes the local stiffness matrix of the axially loaded reinforcing bar.
Eq. 2.31 can also be expressed relative to the global coordinate system by applying a
rotation with angle + , which is the angle between the axis of the reinforcing bar and the
global x-axis of the structure. This can be formally expressed by the following relation:
40 CHAPTER 2

ì P1x ü ì d1 x ü
ïP ï ïd ï
ï 1y ï ï 1y ï
í ý = [T1 ] ⋅ [ K LO ]s ⋅ [T1 ] ⋅ í ý
T
(2.32)
ï P2 x ï ïd2 x ï
ïî P2 y ïþ ïî d 2 y ïþ

where

é cos + sin + 0 0 ù
[T1 ] = ê (2.33)
ë 0 0 cos + sin + úû
1 and 2 refer to the end points of the steel element and the subscripts x and y indicate that
steel forces and displacements are defined with reference to the global coordinate system of
the structure.
Since the end points of the reinforcing bar element do not generally coincide with the
nodes of the concrete element in Figure 2.14, Eq. 2.32 has to undergo another transformation
before it can be assembled together with the concrete element stiffness matrix. This can be
formally expressed by the following relation:

[ KGL ]s = [T2 ]T ⋅ [T1 ]T ⋅ [ K LO ]s ⋅ [T1 ] ⋅ [T2 ] (2.34)

SIDE 3
steel element
SIDE 2 nodes
SIDE 4

SIDE 1

5 5
7 6 7 6
j j
8 4 l2 8 4 l2
l1 c2 c2
i
c1 1 2 3 1 i 2 3
c1
l1

FIGURE 2.14 STEEL ELEMENT EMBEDDED IN CONCRETE ELEMENT

CHAPTER 2 41

Transformation matrix [T2 ] can be derived with the procedure used to establish the consistent
nodal forces of the finite element method. The consistent nodal forces for node i are:

æ Ix Iy ö
1
Pxi = òN
−1
i ⋅ ç Pt ⋅
è I;
− Pn ⋅ ÷ d ;
I; ø
(2.35)

æ Ix Iy ö
1
Pyi = òN
−1
i ⋅ ç Pn ⋅
è I;
− Pt ⋅ ÷ d ;
I; ø
(2.36)

where N i is the shape function for node i and Pn , Pt are the normal and tangential boundary
loads, respectively. If a concentrated load acts at a point ; = -i , Eqs. 2.35 and 2.36 reduce to
the following:
Pxi = N i (- i ) ⋅ Px (2.37)
Pyi = N i (- i ) ⋅ Py (2.38)

When the 8-node isoparametric element is used in the two-dimensional mesh

representation of the member, the shape functions for nodes i, j and k in Figure 2.15 are:
1
Ni = ; ⋅ ( ; − 1) (2.39)
2
N j = (1 − ;) ⋅ (1 + ;) (2.40)

1
Nk =
; ⋅ (; + 1) (2.41)
2
To obtain N i (- i ) in Eqs. 2.37 and 2.38 -i = (2 ci li − 1) is substituted for ξ in Eqs.
2.39-2.41 where ci characterizes the position of the steel node relative to the concrete
element boundary, as shown in Figs. 2.15 and 2.16. By setting ci li = ri Eqs. 2.39-2.41
simplify to the following relations which are used to construct the transformation matrix
[T2 ] .
N i = ( ri − 1) ⋅ (2ri − 1) (2.42)
N j = 4 ri ⋅ (1 − ri ) (2.43)

N k = ri ⋅ (2ri − 1) (2.44)

If the reinforcing bar element crosses the concrete element boundary on sides 2 and 4
in Figure 2.14, nodes i, j and k correspond to node numbers 1, 8 and 7, on side 4, and node
42 CHAPTER 2

numbers 3, 4 and 5 on side 2, respectively. With the notation of Figure 2.14 the
transformation matrix [T2 ] has the following form:

é A1 0 0 0 0 0 A3 A2 ù
[T2 ] = ê (2.45)
ë 0 0 B1 B2 B3 0 0 0 úû

where

é2 p2 − 3 p + 1 0 ù é −4 p 2 + 4 p 0 ù é2 p2 − 2 p 0 ù
A1 = ê ú A = ê ú A = ê ú
2 p 2 − 3 p + 1û −4 p 2 + 4 p û 2 p2 − 2 pû
2 3
ë 0 ë 0 ë 0
(2.46)
é 2q 2 − 3q + 1 0 ù é −4q 2 + 4q 0 ù é 2q 2 − 2q 0 ù
B1 = ê ú B = ê ú B = ê ú
2q − 3q + 1û −4q + 4q û 2q − 2q û
2 2 2 3 2
ë 0 ë 0 ë 0
(2.47)
p = c1 l1 , q = c2 l2 and 0 is the 2x2 null matrix.

in Eq. 2.45 the position of submatrices A and B within the transformation matrix [T2 ]
is related to the side of the concrete element which the reinforcing bar element crosses. If the
reinforcing bar element crosses the concrete element boundary on sides 1 and 2 in Figure
2.14, the transformation matrix [T2 ] takes the form

é A1 A2 A3 0 0 0 0 0ù
[T2 ] = ê (2.48)
ë0 0 B1 B2 B3 0 0 0úû
noting that the submatrices A and B are the same as before, but c1 is now defined, as shown
on the right hand side of Figure 2.14.
When the 4-node isoparametric element is used in the two-dimensional mesh
representation of the member, the shape functions for nodes i, and k change and the
transformation matrix [T2 ] simplifies considerably. For example, if the reinforcing bar
element crosses the concrete element boundary on sides 2 and 4 in Figure 2.14, the
transformation matrix [T2 ] takes the form

é1 − p 0 0 0 0 0 p 0ù
ê 0 1− p 0 0 0 0 0 pú
[T2 ] = ê ú (2.49)
ê 0 0 1− q 0 q 0 0 0ú
ê ú
ë 0 0 0 1− q 0 q 0 0û
where p and q are defined as in Eqs. 2.46 and 2.47.
CHAPTER 2 43

i 1
Ni
ci

j steel Nj

li nodes Nk
1

k
-1 0 1

(a) (b)

FIGURE 2.15 PARAMETERS AND SHAPE FUNCTIONS OF ISOPARAMETRIC

CONCRETE AND EMBEDDED STEEL ELEMENT

The steel stiffness matrix [ KGL ]s in Eq. 2.34 can now be assembled together with the
concrete element stiffness matrix to form the total stiffness of the structure.
In the analysis of slabs with the layer model the effect of bond-slip is not taken into
account. In this case it is computationally convenient to use the distributed steel model. In
this model the reinforcing bars inside a concrete element are replaced by an equivalent steel
element with distributed uniaxial material properties in each reinforcing direction. The
equivalent steel element has the dimensions of the concrete element and thickness ts which is
determined from the following relation:
As
ts =
= 5s ⋅ d e (2.50)
b
where As is the cross sectional area of one reinforcing bar in the particular direction, b is the
spacing of reinforcing bars, 5 s is the reinforcing ratio and d e is the effective depth of the
member.
Since the equivalent steel element has uniaxial properties in the direction parallel to
the axis of the reinforcing bars, the constitutive material model takes the simple form
44 CHAPTER 2

y x

x
0

FIGURE 2.16 COORDINATE TRANSFORMATION FOR STEEL

ì 61 ü é Es1 0 0ù ì .1 ü
ï ï ê ú ï ï
í 6 2 ý = ê 0 0 0ú ⋅ í . 2 ý (2.51)
ï7 ï ê 0 0 0ú ï 0 ï
î 12 þ ë û î 12 þ
where Es1 is Young's modulus of reinforcing steel. Es1 in Eq. 2.51 is replaced by the strain
hardening modulus Es 2 , when the reinforcing steel yields.

The local stiffness matrix in Eq. 2.51 needs to be transformed to the global x-y
coordinate system before it can be assembled into the structure stiffness matrix. This is
accomplished exactly as in Eqs. 2.25 and 2.26 for concrete, where + now represents the
angle between the axis of the reinforcing bars and the global x-axis (Figure 2.16).

2.4 Bond between Reinforcing Steel and Concrete

2.4.1 Bond Behavior

Bond is the interaction between reinforcing steel and surrounding concrete. The force
transfer from steel to concrete can be attributed to three different phenomena: (1) chemical
adhesion between mortar paste and bar surface; (2) friction and wedging action of small
dislodged sand particles between the bar and the surrounding concrete; and (3) mechanical
CHAPTER 2 45

interaction between concrete and steel. Bond of plain bars derives primarily from the first two
mechanisms, even though there is some mechanical interlocking caused by the roughness of
the bar surface. Deformed bars have better bond than plain bars, because most of the steel
force is transferred through the lugs to concrete. Friction and chemical adhesion forces are
not negligible, but secondary and tend to decrease as the reinforcing bars start to slip.
Since bond stresses in reinforced concrete members arise from the change in the steel
force along the length, the effect of bond becomes more pronounced at end anchorages of
reinforcing bars and in the vicinity of cracks.

τ2
Eb2
τ1 1

Eb1

- dsu 0 1
ds
ds1 d s2 d su

−τ
1

−τ
2

FIGURE 2.17 BOND STRESS-SLIP RELATION FOR PLANE STRESS PROBLEMS

In the simplified analysis of reinforced concrete structures complete compatibility of

strains between concrete and steel is usually assumed, which implies perfect bond. This
assumption is realistic only in regions where negligible stress transfer between the two
components takes place. In regions of high transfer stresses along the interface between
reinforcing steel and surrounding concrete, such as near cracks, the bond stress is related to
the relative displacement between reinforcing steel and concrete. The assumption of perfect
bond near crack zones requires infinitely high strains to explain the existence of a finite crack
width. In reality there is no strain compatibility between reinforcing steel and surrounding
46 CHAPTER 2

concrete near cracks. This incompatibility and the crack propagation give rise to relative
displacements between steel and concrete, which are known as bond-slip.
Two types of experiments are used in the determination of the bond stress-slip
relation. In an anchorage test, such as the standard ASTM pull-out test, the force is applied at
the projecting end of a bar which is embedded in a concrete cylinder. This force is gradually
transferred to the concrete. In the transfer test a self-equilibrating pair of forces are applied at
the two ends of a bar projecting from a concentrically cast concrete block or cylinder. Tests
of this type are intended to simulate conditions in the tension zone of a concrete beam
between primary flexural cracks.
In either type of test the bond stress is determined from the change in steel stress over
a certain measurement length, which is usually taken equal to five bar diameters, and the
relative slip is determined externally or internally. It is, therefore, practically impossible to
establish a "local" bond stress-slip relation, since the measured bond stress-slip relation
generally represents the average relation over the measurement length. Moreover, since the
bond stress is derived from the change in steel stress, the result is very sensitive to
experimental error. The bond-slip relation also depends on the position of the bars, the bar
surface conditions, the loading state, the boundary conditions and the anchorage length of
bars. In spite of these difficulties, several experimental bond-stress slip relations have been
proposed (ASCE 1982; Eligehausen et al. 1983; Hayashi and Kokusho 1985).
In this study the simple trilinear bond stress-slip model in Figure 2.17 is adopted in
the finite element analysis of plane stress problems. The parameters of the model are derived
from the material properties of each specimen in the correlation studies of Chapter 4. This
model is a good approximation of the actual behavior in cases, which do not exhibit
significant bond-slip and associated bond damage. Under monotonic loading this holds true
in all RC members which do not experience anchorage failure.
Since it is the objective of this study to investigate the bond-slip behavior of
reinforcing steel in more detail, a more sophisticated bond-slip model is used in the
correlation studies of anchored reinforcing bars under monotonic and cyclic loads in Chapter
4. A brief description of this model follows, while more details are given elsewhere
(Eligehausen et al. 1983; Zulfiqar and Filippou 1990).
The bond stress-slip relation between reinforcing bars and surrounding concrete
consists of the following parts:
CHAPTER 2 47

a,b MONOTONIC ENVELOPES

τ
c UNLOADING
d RELOADING τ1
e REDUCED ENVELOPES
f FIRST RELOADING a

e
τ3
d
τf
u u u u
1 2 3
c
f

e
b

FIGURE 2.18 BOND STRESS-SLIP RELATION FOR ANCHORAGE PROBLEMS

Two monotonic envelopes, one in tension and one in compression, which are updated
at each slip reversal to reflect incurred damage (curves (a) and (b) in Figure 2.18).
A typical unloading-reloading path described by the unloading curve (c), the current
frictional bond resistance 7 f and a reloading curve (d), along with a set of rules for
unloading and reloading for incomplete cycles (Figure 2.18).
A set of relations for updating the monotonic envelope values and the frictional bond
resistance as a function of incurred damage.
The monotonic envelope consists of an initial nonlinear relation 7 = 71 ⋅ ( u u1 ) , valid
*

for u ≤ u1 , followed by a plateau 7 = 71 for u1 ≤ u ≤ u2 . For u ≥ u2 , 7 decreases linearly to

the value of ultimate frictional bond resistance 7 3 at a slip value of u3 which is assumed to
be equal to the clear distance between the lugs of deformed bars.
When the load is reversed at some slip value, unloading takes place along a steep
straight line until the frictional bond resistance 7 f is reached (curve c Figure 2.18).

The reloading curve is described by a fourth degree polynomial (curve d in Figure

2.18) starting from the frictional bond resistance value and leading to the intersection
48 CHAPTER 2

between reduced envelope and previous unloading curve. In case that no slip has been
previously imposed in the same loading direction, reloading takes place along a horizontal
line until reaching the reduced envelope (curve f in Figure 2.18).
During each half cycle following the first unloading, the monotonic envelopes are
updated (curve e Figure 2.18) by reducing the characteristic bond stress values 71 and 7 3 by
a factor, which is a function of the "damage parameter" d. The damage parameter has the
form (Eligehausen et al. 1983)
71 ( N ) = 71 ⋅ (1 − d ) (2.52)
where 71 is the characteristic value of the virgin envelope curve and 71 ( N ) is the
corresponding value after N cycles. The damage parameter depends on the total energy
dissipated by the bond-slip process and is given by
−1.2( E Eo )
1.1

d = 1− e (2.53)
where E is the total dissipated energy. Eo is the energy absorbed under monotonically
increasing slip up to the value u3 and is used as a normalization parameter.
The frictional bond resistance 7 f depends on the maximum value of previous slip
umax in either direction of loading and the ultimate bond resistance 7 3 ( N ) of the
corresponding reduced envelope curve. If cyclic reversals are performed under fixed values
of slip, 7 f is further reduced by multiplying its initial value with a factor which depends on
the energy dissipated by friction alone. Explicit expressions for these relations are presented
by Eligehausen et al. (1983).
In the nonlinear bond stress-slip model all expressions are cast in dimensionless form
except for the characteristic values of the monotonic envelope curve. The model can thus be
readily extended to describe the local bond stress-slip relation of reinforcing bars under
different bond and loading conditions by only changing the characteristic values of the
monotonic envelope curve. These characteristic values depend on several parameters such as
bar diameter, concrete strength, bar deformation pattern, clear spacing between bars and
transverse pressure due to loading or confinement. Ideally, the characteristic values should be
derived from experiments which simulate the actual geometric and loading conditions of the
structure under investigation. In the absence of experimental data the characteristic values of
the envelope curve can be derived from the material properties of the structure based on
recommendations of Eligehausen et al. (1983).
CHAPTER 2 49

2.4.2 Bond Stiffness Matrix

Two basically different elements have been proposed to date for including the bond-
slip effect in the finite element analysis of RC structures. The first is the bond link element by
Ngo and Scordelis (1967) shown in Figure 2.19. It consists of two orthogonal springs which
connect and transmit shear and normal forces between a reinforcing steel node and an
adjacent concrete node. Since the link element has no physical dimensions, the two connected
nodes originally occupy the same location in the finite element mesh of the undeformed

concrete node

r
θ

steel node

FIGURE 2.19 BOND LINK ELEMENT

structure.
The second bond-slip element, called the bond zone element, is significantly different
from the bond link element, the most important difference being its finite dimension. In the
bond zone element by de Groot et al. (1981) the contact surface between reinforcing steel and
concrete and the concrete in the immediate vicinity of the reinforcing bar are modeled by a
material law, which represents the special properties of the bond zone. In this model, which is
shown in Figure 2.20, the bond stress is the sum of the slip resistance and the stress due to
mechanical interlocking and can be determined from
7 = 7 o + 7 k = S ⋅ u − k ⋅ 6 r (2.54)
50 CHAPTER 2

where S is the slip modulus, u is the relative slip, k is the rib factor which represents the
effect of the bar deformations and 6r is the radial stress.

v v
node

u u

concrete bond-slip
element element
(a)

1 2 = translation
2 3 = rotation

3
concrete 1
2
bond-slip
layer
1

reinforcement

(b)

FIGURE 2.20 BOND ZONE ELEMENT

The contact element describes the behavior of the layer between reinforcing steel and
concrete and provides a continuous connection between the reinforcing steel and the
surrounding concrete element, if a linear or higher order displacement field is used in the
discretization scheme. A simpler but similar element was proposed by Keuser and Mehlhorn
(1987). They introduced higher order displacement interpolation functions and showed
through energy considerations that the bond link element cannot adequately represent the
stiffness of the steel-concrete interface. Finally, Yankelevsky (1985) derived closed form
solutions of the stiffness matrix of a reinforcing bar element with bond-slip on the
assumption that the steel remains linear elastic and the local bond-stress slip relation is
approximated with a piecewise linear curve.
In studies where detailed local behavior is of interest, the continuous bond models are
most appropriate. In cases, however, where the overall structural behavior is of primary
interest, the bond link element provides a reasonable compromise between accuracy and
CHAPTER 2 51

computational efficiency. It is, therefore, selected for representing the bond-slip effect in this
study.
The element stiffness matrix relates shear and normal force to the corresponding
nodal displacements by the following relation (ASCE 1982):

ì Fr ü é K r 0 ù ìd r ü
í ý=ê ⋅í ý (2.55)
î Fs þ ë 0 K s úû î d s þ
where K r represents the stiffness due to dowel action and may be assigned a very low value,
if dowel action is not accounted for. K s is the shear stiffness of the interface and can be
derived from the measured bond stress-slip relation according to the following equation
K s = Eb1 ⋅ A (2.56)
where Eb1 is the initial slip modulus and is replaced by Eb 2 when the bond-slip d s exceeds
the value d s1 (Figure 2.17). A is the bar circumferential area tributary to one bond link
element and is determined from the following relation:
m3d bl
A= (2.57)
2b
where m is the number of bars of diameter d b in the cross section of the member, l is the
spacing of the bond links along the reinforcing bar and b is the width of the member cross
section. The factor 2 appears in the denominator to account for the fact that it is usually
convenient to place bond link elements at, both, the top and the bottom of the reinforcing bar
element.
The local stiffness matrix of the bond link element in Eq. 2.55 is transformed to
global coordinates by a rotation matrix. This can be expressed formally as

[ KGL ]B = [T ]T ⋅ [ K LO ]B ⋅ [T ] (2.58)

where

é Kr 0 −Kr 0 ù
ê 0 Ks 0 − Ks ú
[ K LO ]B =ê ú (2.59)
ê− Kr 0 Kr 0 ú
ê ú
ë 0 − Ks 0 Ks û

and
52 CHAPTER 2

é c s 0 0ù
ê − s c 0 0ú
[T ] = ê ú (2.60)
ê 0 0 c sú
ê ú
ë 0 0 −s cû
with c = cos + , s = sin + and + as defined in Figure 2.19.

The global bond link stiffness matrix in Eq. 2.59 takes the following explicit form

é Kr ⋅ c2 + Ks ⋅ s2 K r ⋅ sc − K s ⋅ sc − K r ⋅ c 2 − K s ⋅ s 2 − K r ⋅ sc + K s ⋅ sc ù
ê ú
K r ⋅ s 2 + K s ⋅ c 2 − K r ⋅ sc + K s ⋅ sc − K r ⋅ c 2 − K s ⋅ s 2 ú
[ KGL ]B = ê
ê
(2.61)
Kr ⋅ c2 + K s ⋅ s2 K r ⋅ sc − K s ⋅ sc ú
ê ú
ëê sym K r ⋅ s 2 + K s ⋅ c 2 ûú

2.5 Discrete Reinforcing Steel Model with Bond-Slip

The bond link element connects one steel node with a corresponding concrete node
which occupies the same physical location in the undeformed configuration of the structure.
Consequently, the use of this element in the finite element analysis of RC structures imposes
the following restrictions: (a) the finite element mesh must be arranged so, that a reinforcing
bar is located along the edge of a concrete element and (b) a double node is required to
represent the relative slip between reinforcing steel and concrete. These restrictions arise
from the fact that the stiffness of the bond link element is associated with the relative
displacement between steel and concrete. Consequently, the total displacement of, both,
reinforcing steel and concrete is required at each node of the finite element mesh, so that the
relative displacement and, consequently, the bond stress between steel and concrete can be
determined. In a complex structure, particularly in three-dimensional models, these
requirements lead to a considerable increase in the number of degrees of freedom, not only
because of the doubling of the number of nodes along the reinforcing steel bars, but also
because the mesh has to be refined, so that the bars pass along the edges of concrete
elements. The complexity of mesh definition and the large number of degrees of freedom has
discouraged researchers from including the bond-slip effect in many studies to date.
To address some of these limitations of the bond link element a new discrete
reinforcing steel model which includes the bond-slip deformation is proposed in this study. In
this model the reinforcing bar is assumed to be embedded inside the concrete element, as
shown in Figure 2.21a, so that the analyst can choose the finite element mesh configuration
CHAPTER 2 53

independently of the location of the reinforcing bars. At the same time the relative slip
between reinforcing steel and concrete is explicitly taken into account in the model. Since the
finite element model only includes the concrete displacement degrees of freedom, the degrees
of freedom which are associated with the reinforcing steel need to be condensed out from the
element stiffness matrix, before it is assembled into the structure stiffness matrix. A second
pass, however, is now required in order to satisfy the equilibrium of each reinforcing bar and
determine the steel forces. The formal derivation of this model is presented in the following,
while the numerical implementation is deferred to Chapter 3.
A reinforcing steel element which is embedded in a concrete element is shown in
Figure 2.21a. A convenient free body diagram is selected which isolates the steel element
with the bond link elements attached at its end points. Figure 2.21b shows this element before
and Figure 2.21c after deformation. i and j denote the end points of the element, points 1 and
3 are associated with the concrete and points 2 and 4 are associated with the reinforcing steel
at ends i and j, respectively. The corresponding degrees of freedom of the reinforcing steel
and concrete at each end are connected by the bond link element whose stiffness depends on
the relative displacement between steel and concrete. With this assumption the stiffness
matrix which relates the end displacements along the axis of the reinforcing bar with the
corresponding forces is given by

ì P1 ü é kbi − kbi 0 0 ù ì d1 ü
ï P ï ê −k k s + kbi − kbj ú ïïd 2 ïï
ï 2 ï ê bi 0
ú⋅í ý
í ý= (2.62)
ï P3 ï ê 0 0 kbj −k s ú ï d 3 ï
ïî P4 ïþ êêë 0 − kbj −ks
ú
k s + kbj úû ïîd 4 ïþ
where k s = EA L is the steel bar stiffness, kb is the stiffness of the bond link parallel to the
bar axis at the corresponding end of the steel element and the dowel action is neglected in Eq.
2.62 so that term k r in Eq. 2.59 is equal to zero. In order to include the effect of dowel action
Eq. 2.62 has to be expanded to a relation between eight displacement degrees of freedom and
the corresponding forces. Since the following derivation procedure is not affected by this
change, the dowel action is neglected here for the sake of simplicity.
By rearranging rows and columns in Eq. 2.62 the following relation results
54 CHAPTER 2

ì P1 ü é kbi 0 − kbi 0 ù ì d1 ü
ïP ï ê 0 − kbj ú ïï d 3 ïï
ï 3ï ê kbj 0
ú⋅í ý
í ý= (2.63)
ï P2 ï ê − kbi 0 k s + kbi −ks ú ïd 2 ï
ê ú
îï P4 þï êë 0 − kbj −ks k s + kbj úû îïd 4 þï

or

ì Pc ü é K cc K cs ù ì d c ü
í ý=ê ⋅í ý (2.64)
î Ps þ ë K cs K ss úû î d s þ

where the subscript c stands for concrete and s stands for reinforcing steel.

steel element
j
concrete element

i
(a) embedded steel element

2 4
j (b) before deformation
i
1 3
steel element

2 4
(c) after deformation

1 3
1 , 3 : concrete d.o.f
bond element

FIGURE 2.21 DISCRETE REINFORCING STEEL ELEMENT WITH BOND-SLIP

Since the finite element model only includes the concrete displacement degrees of
freedom, the degrees of freedom which are associated with the reinforcing steel need to be
condensed out from the element stiffness matrix, before it is assembled into the structure
stiffness matrix. By applying static condensation of the steel degrees of freedom in Eq. 2.64
the following relation between concrete displacements and corresponding forces results

{P } = éë K ùû ⋅ {d }
c
* *
cc c (2.65)
CHAPTER 2 55

where

{P } = {P } − [ K ] ⋅ [ K ] ⋅ {P }
c
*
c cs ss
−1
s (2.66)

éë K cc* ùû = [ K cc ] − [ K cs ] ⋅ [ K ss ] ⋅ [ K cs ]
−1
(2.67)

After some calculations for evaluating the inverse and carrying out the multiplications Eq.
2.67 reduces to:

k s ⋅ kbi ⋅ kbj é 1 −1ù

éë K cc* ùû = ⋅ê ú = éë K eq ùû s (2.68)
k s ⋅ ( kbi + kbj ) + kbi ⋅ kbj ë −1 1 û

which is the local stiffness matrix of the reinforcing steel element including the effect of
bond slip and is equivalent to the local stiffness matrix in Eq. 2.30. The coefficient of the
equivalent stiffness matrix can be rewritten as
k s ⋅ kbi ⋅ kbj ks
⋅= (2.69)
k s ⋅ ( kbi + kbj ) + kbi ⋅ kbj æ 1 1 ö
1 + ks ⋅ ç +
ç kbi kbj ÷÷
è ø
It is now readily apparent from Eq. 2.69 that bond slip reduces the stiffness of the reinforcing
steel element. In case of perfect bond the bond stiffness terms kbi and kbj become infinitely
large and the expression in Eq. 2.69 reduces to k s . In this case the stiffness matrix in Eq. 2.68
reduces to the local stiffness matrix of the embedded steel model with perfect bond in
Eq. 2.30.
The process of transforming the local stiffness matrix of the reinforcing steel element
to the global coordinate system is identical to that presented in Section 2.3.2 for the
embedded steel element with perfect bond. Since the end points of the reinforcing steel
element do not coincide with the nodes of the concrete element mesh, two transformation
matrices are needed. By replacing [ K LO ]s in Eq. 2.34 with éë K eq ùû s the transformation to

global coordinates of the local stiffness matrix of the reinforcing steel element with bond slip
takes the form

[ KGL ]s = [T2 ]T ⋅ [T1 ]T ⋅ éë K eq ùû s ⋅ [T1 ] ⋅ [T2 ] (2.70)

The steel stiffness matrix [ KGL ]s in Eq. 2.70 can now be assembled together with the
concrete element stiffness matrix to form the total stiffness of the structure.
56 CHAPTER 2

Since the degrees of freedom associated with the reinforcing steel are condensed out
of the stiffness matrix in Eq. 2.70, an additional step now becomes necessary for determining
the deformations and forces in the reinforcing steel. Once the displacement increments at the
nodes of the concrete finite elements are determined for the current load increment, they can
be transformed by matrices [T1 ] and [T2 ] to yield the concrete displacements {d c } at the ends
of the steel element. The vector {d c } has only two components which are the end
displacements parallel to the axis of the bar. Using the second row of the matrix relation in
Eq. 2.64

{Ps } = [ K cs ] ⋅ {d c } + [ K ss ] ⋅ {d s } (2.71)
and solving for {d s } yields

[ K ss ] ⋅ {d s } = ({Ps } − [ K cs ] ⋅ {d c }) (2.72)

Eq. 2.72 expresses the equilibrium of the reinforcing bar element. After assembling
the steel element stiffness matrices on the left hand side and the forces on the right hand side
of the equation and imposing appropriate boundary conditions at the ends of the reinforcing
bar the simultaneous solution of the equations of equilibrium yields the steel deformations
and forces. Alternatively, it is possible to rewrite Eq. 2.72 so as to yield the steel force and
deformation at one end of the steel element as a function of the force and deformation at the
other end. The successive application of this relation along with appropriate conditions for
the transition from one element to the next results in a transfer matrix solution method, which
offers certain computational advantages over the direct solution method. The numerical
implementation of these methods and their advantages and disadvantages will be discussed in
the following chapter.
 CHAPTER 3
FINITE ELEMENT MODELING OF RC BEAMS AND SLABS

3.1 General

The finite element method is a general method of structural analysis in which the
solution of a problem in continuum mechanics is approximated by the analysis of an
assemblage of finite elements which are interconnected at a finite number of nodal points and
represent the solution domain of the problem.

The finite element method is now well accepted as the most powerful general
technique for the numerical solution of a variety of engineering problems. Applications range
from the stress analysis of solids to the solution of acoustical phenomena, neutron physics
and fluid dynamic problems. Indeed the finite element method is now established as a general
numerical method for the solution of partial differential equations subject to known boundary
and initial conditions.

In the realm of linear analysis the finite element method is now widely used as a
design tool. A similar acceptance in nonlinear analysis problems depends on two major
factors. First, the increase in computational effort which is required for nonlinear problems
necessitates that considerable computing power be available at low cost to the designer.
Developments in the last two decades have ensured that high-speed digital computers have
gradually become available to the average designer and indications are that reductions in unit
computing costs will continue at an accelerating pace.

The second major factor is related to the level of complexity of nonlinear analysis.
Before application of nonlinear methods can become commonplace in design situations, the
accuracy and reliability of the proposed models has to be established beyond doubt. The
development of improved element characteristics and more efficient nonlinear solution
algorithms as well as the experience gained in their application to engineering problems have
ensured that nonlinear finite element analysis can now be performed with some confidence.
Thus, barriers to the wide use of nonlinear finite element techniques are gradually removed.

57
58 CHAPTER 3

Nevertheless, difficulties still abound whose solution will require much effort on the part of
researchers and designers.

The behavior of RC structures is distinctly nonlinear, because of several factors:

(1) nonlinear material behavior of concrete and steel and their interaction through bond and
dowel action; (2) cracking of concrete; and (3) time dependent effects such as creep,
shrinkage, temperature and load history. In dealing with these problems in an analytical
setting several nonlinear solution algorithms are now available.

This chapter first presents the underlying theoretical assumptions of the RC beam and
slab model used in this study. This is followed by the solution of the problem in the finite
element context. Following a brief description of the nonlinear solution algorithm a detailed
step-by-step organization of the nonlinear solution strategy is presented.

3.2 Theoretical Considerations

For the analysis of RC beams the Timoshenko beam theory is used in the present
study. Since this theory is well established and widely used in the analysis of beams, attention
is focused below on some theoretical aspects of the plate bending problem followed by the
finite element implementation of the material models of Chapter 2 in the analysis of RC
beams and slabs.

θx
dw
dx φ
x X

ACTUAL
MID-SURFACE
DEFORMATION

NORMAL TO MID-SURFACE
ASSUMED DEFORMATION

FIGURE 3.1 ASSUMED DEFORMATION OF MINDLIN PLATE THEORY

CHAPTER 3 59

In the analysis of thin plates it is generally assumed that the plate thickness is small
relative to the dimensions in the x-y plane, which is the plane of the plate. When plates are
subjected to in-plane loads the resulting stresses are assumed to be constant through the
thickness of the plate and stresses 6 z , 7 xz and 7 yz are ignored. For out-of-plane loading the
Mindlin plate bending theory, which includes the effect of transverse shear deformations, is
adopted in this study (Owen and Hinton 1978). The main assumptions of the theory are:

(1) transverse displacements are small relative to the plate thickness,

(2) the stress normal to the mid-surface of the plate is negligible, and

(3) normals to the mid-surface of the slab before deformation remain straight, but not
necessary normal to the mid-surface after deformation, as shown in Fig. 3.1. The force
resultants and corresponding deformations of a typical Mindlin plate are shown in Fig.
3.2.

The deformation field of the problem consists of three components:

{d } = {w + y}
T
+x (3.1)

Qx
z,w
M xy
Qy
Mx
y
θx
M yx
My
My
M yx
x
Mx

Qy
θy M xy
Qx

FIGURE 3.2 TYPICAL MINDLIN PLATE

where w is the transverse plate displacement normal to the xy plane and + x and + y are the
rotations about the y-and the x-axis, respectively (Fig. 3.2). Under the assumptions of the
Mindlin plate theory it holds that
60 CHAPTER 3

sw
+x = − Kx (3.2)
sx
sw
+y = − Ky (3.3)
sy
where K x and K y are transverse shear rotations (Fig. 3.1).

Making use of the basic assumptions of the theory the displacements of any point in
the plate can be expressed by:

u( x, y , z ) = uo ( x, y ) + [d nx ( x, y ) − z ] ⋅ + x (3.4)

v ( x, y , z ) = vo ( x, y ) + éë d ny ( x, y ) − z ùû ⋅ + y (3.5)

w( x, y , z ) = wo ( x, y ) (3.6)
where d nx and d ny is the distance from the top surface of the plate to a point in the plate with
in-plane deformations of uo and vo , respectively. It should be noted that d nx and d ny are not
necessarily equal.
The bending strains . x , . y and 0 xy vary linearly through the plate thickness and can
be derived by taking derivatives of Eqs. 3.4-3.6:
su suo s+ s+
Fx = = + ( d nx − z ) ⋅ x = Fox + ( d nx − z ) ⋅ x (3.7)
sx sx sx sx
sv svo s+ s+
Fy = = + ( d ny − z ) ⋅ y = Foy + ( d ny − z ) ⋅ y (3.8)
sy sy sy sy
su sv suo svo s+ s+
H xy = + = + + ( d nx − z ) ⋅ x + ( d ny − z ) ⋅ y
sy sx sy sx sy sx
(3.9)
s+ s+
= H oxy + ( d nx − z ) ⋅ x + ( d ny − z ) ⋅ y
sy sx
or

{. 0 xy } = {. p }
T
x .y (3.10)
Since 0 xy varies linearly with distance z, Eq. 3.9 can be simplified to the following relation

æ s+ s+ y ö
H xy = H oxy + ( d nxy − z ) ⋅ ç x + (3.11)
è sy sx ÷ø
CHAPTER 3 61

The transverse shear strains, on the other hand, are constant through the thickness and can be
expressed as follows:
sw
H zx = − +x (3.12)
sx
sw
H zy = − +y (3.13)
sy
or

{0 xz 0 zy } = {. t }
T
(3.14)

CONCRETE LAYER
j-th layer
dn

Zj

NEUTRAL AXIS

EFFECTIVE STEEL LAYER

FIGURE 3.3 LAYER SECTION

Finally, d nx , d ny and d nxy , which represent the arbitrary depths of .ox , .oy , and 0 oxy ,
respectively, (Figure 3.3) can be determined from the condition that
ò 6 x dz = ò 6 y dz = ò 7 xy dz = 0 , which decouples the in-plane from the out-of-plane forces.
Assuming that the shear modulus G is constant through the depth of the slab leads to the
conclusion that d nxy is equal to d/2, where d is the slab thickness.

The bending strain components of Eqs. 3.7-3.9 can be expressed as follows:

62 CHAPTER 3

ì s+ ü
ï ( d nx − z ) ⋅ x ï
ì Fx ü ï sx ï
ï s+ ïï
{F p } = ïí F y ïý = ïí ( d ny − z ) ⋅ syy ý (3.15)
ïH ï ï ï
î xy þ ï
æ s+ s+ y öï
ï( d nxy − z ) ⋅ ç x + ÷ï
îï è sy sx ø þï
The effect of 7 zx and 7 zy was ignored in this derivation. This is a reasonable
approximation, since, in the absence of in-plane forces, the in-plane stresses reach a
maximum at the extreme fibers, where the transverse shear stresses are at their lowest value,
and reach a minimum at mid-depth of the slab, where the transverse shear stresses are
highest. Thus, in the absence of in-plane forces the interaction between transverse shear
stresses and in-plane stresses is relatively small and can be neglected.

3.3 Finite Element Representation

The finite element analysis of a continuum starts with the subdivision of the physical
system into an assemblage of discrete elements. The displacement vector {d } at any point
within a particular structure is approximated by interpolating functions associated with
generalized coordinates d i , which are the displacements of the nodes of the finite element
discretization of the structure. For a finite element this approximation can be formally written
as

ì d1 ü
ïMï
ï ï
{d } = [ N ] ⋅ {d } = [ N1 K N i K] ⋅ í ý
ele
(3.16)
ïdi ï
ïî M ïþ

where the components of [ N ] are, in general, functions of position and {d }

ele
is the vector of
the node displacements of a particular element. Once the displacements are known at any
point within the element, the strains can be determined from the strain-displacement relation.
In the discrete problem this can be formally written as

{.} = [ B ] ⋅ {d }
ele
(3.17)
CHAPTER 3 63

The matrix [ B ] will be derived later for beams and slabs separately. The stresses can now be
determined from the material constitutive law

{6} = [ D ] ⋅ ({.} − {.o }) + {6o } (3.18)

where [ D ] is the element material matrix, {.o } is the initial strain vector and {6o } is the
initial stress vector.
By applying the virtual work principle or the theorem of minimum potential energy to
the assemblage of discrete elements the following equilibrium equations result

[ K ] ⋅ {d } + {F }p + {F }g + {F }. + {F }6 − {R} = 0
o o
(3.19)

The terms in Eq. 3.19 are derived as follows:

the stiffness matrix [ K ] ,

[ K ] = å ò [ B ]T ⋅ [ D ] ⋅ [ B ] dV (3.20)
ele

the nodal forces due to surface traction,

{F }p = −å ò [ N ]T ⋅ { p}dV (3.21)
ele

the nodal forces due to body forces,

{F }g = − å ò [ N ]T ⋅ {g }dV (3.22)
ele

the nodal forces due to initial strains,

{F }. = − å ò [ B ] ⋅ [ D ] ⋅ {.o } dV
T
(3.23)
o
ele

and the nodal forces due to initial stresses

{F }6 = − å ò [ B ] ⋅ [ D ] ⋅ {6o } dV
T
(3.24)
o
ele

In Eqs. 3.19-3.24 {d } is the vector of node displacements, {R} is the vector of applied nodal
forces, { p} is the vector of surface forces and {g} is the vector of body forces.

The node displacements {d } can be determined from the solution of the system of
simultaneous algebraic equations in Eq. 3.19, whereupon the strains and stresses at any point
of the structure can be obtained from Eqs. 3.17 and 3.18. In a nonlinear problem the stiffness
64 CHAPTER 3

matrix [ K ] depends on the displacement vector {d } and the nonlinear system of algebraic
equations in Eq. 3.19 has to be solved by one of the nonlinear solution methods in Section
3.4.

3.3.1 Reinforced Concrete Beams

The finite element idealization of reinforced concrete beams is accomplished by two-

dimensional elements which lie in the x-y plane of the beam elevation. Following the
assumptions of the Timoshenko beam theory the resulting stresses are assumed to be constant
through the thickness of each element and stresses 6 z , 7 xz and 7 yz are ignored. These
assumptions lead to a plane stress field. The displacement field is described by u and v, which
represent the node displacements in the x and y-direction, respectively, and are in general
functions of x and y. These functions must satisfy C0 continuity. The process of the finite
element discretization of this displacement field can be formally written as

ìu ü n é N j 0 ù ìu ü
{d } = í ý=å ê ⋅í ý
N j úû î v þ j
(3.25)
î v þ j =1 ë 0
where n is the number of nodes of the finite element discretization of the structure and N j is
the shape function for j-th node.
The application of the strain-displacement relation leads to

é IN ù
ê 0 ú
ì . x ü ê Ix ú ele
ï ï ê IN ú ì u ü
í .y ý = ê 0 ⋅í ý (3.26)
ï0 ï ê Iy ú î v þ
ú
î xy þ ê IN IN ú
ê Iy Ix úû
ë
or, in compact form to

{. } = [ B ] ⋅ {d }
p
ele
(3.27)

The stiffness matrix of the composite reinforced concrete finite element is arrived at
by superposition of the concrete and reinforcing steel stiffness matrices derived in Chapter 2.
The steel stiffness matrix of the discrete element also includes the effect of bond-slip, as
described in Section 2.5. This process can be formally expressed as
CHAPTER 3 65

[ K ]ele = [ K ]c + å ([ KGL ]s )
n
(3.28)
i =1

n is the number of steel elements embedded in the concrete element and [ KGL ]s is the global
stiffness matrix of the discrete steel element without bond-slip (Eq. 2.34) or with bond-slip
(Eq. 2.70). The concrete element stiffness [ K ]c is determined from

[ K ]c = ò [ B ]T ⋅ [ DGL ]c ⋅ [ B ] dV (3.29)
V

where [ DGL ]s is the concrete material matrix in global coordinates in Eq. 2.25. The assembly
of the element stiffness matrices [ K ]ele results in the structural stiffness matrix [ K ] .

3.3.2 Reinforced Concrete Slabs

The Mindlin plate element requires only C0 continuity of the transverse displacement
function w and the two independent rotation functions + x and + y . By contrast, the plate
element which is based on the classical Kirchhoff thin plate theory requires C1 continuity,
since the derivatives sw and Iw of the displacement function w as well as w itself
sx Iy
must be continuous across the element interfaces. The proposed Mindlin plate element,
therefore, not only accounts for the transverse shear deformations, but is also simpler to
formulate.
Based on the assumptions of the Mindlin plate theory and the continuity requirements
of the finite element method the displacement field in Eq. 3.1 can be expressed in the
following matrix form

ìwü éN j 0 0 ù ìwü
{d } = ïí + x ïý = å ê ú ï ï
n

ê0 Nj 0 ú ⋅ í+ x ý (3.30)
ï+ ï j =1 êë 0 0 N j úû ïî+ y ïþ j
î yþ
where n is the number of nodes of the finite element discretization of the structure and N j is
the shape function for j-th node.
Using Eq. 3.15 the relation between in-plane strains and displacements becomes
66 CHAPTER 3

é sN ù
ê0 − 0 ú
ì F x ü é( z − d x ) 0 0 ù ê sx ú ìwü
ele

ï ï ê ú ê sN ú ï ï
í Fy ý = ê 0 (z − d y ) 0 ú ⋅ ê0 0 − ⋅ í+ x ý (3.31)
ï ï sy ú ï ï
î H xy þ ëê 0 0 ( z − d xy ) ûú ê ú +y
ê sN sN ú î þ
ê 0 − sy −
sx úû
ë
or, in compact form

{. } = [d ] ⋅ éë B ùû ⋅ {d }
p n p
ele
(3.32)

From Eqs. 3.12 and 3.13 the relation between transverse shear strains and displacements
becomes

é sN ù ele
−N 0 ú ìwü
ì zx ü ê sx
H ï ï
í ý=ê ú ⋅ í+ x ý (3.33)
î H yz þ ê sN 0 − N ú ïî+ y ïþ
ëê sx ûú
or, in compact form

{.t } = [ Bt ] ⋅ {d }
ele
(3.34)

The strain vector results from the combination of Eq. 3.33 with Eq. 3.34
ele
ìwü
ïì{F p }ïü é[d n ] ⋅ ëé B p ûù ù ï ï
{F p } = í ý = ê ú ⋅ í+ x ý (3.35)
{F }
îï t þï ë ê [ B t ] úû ï+ ï
î yþ
With the decomposition of the strain field in Eq. 3.35 the Mindlin plate element
stiffness matrix consists of two parts

[ K ] = ò éë Bp ùû ⋅ [d n ]T ⋅ éë D p ùû ⋅ [d n ] ⋅ éë B p ùû dV + ò [ Bt ]T ⋅ [ Dt ] ⋅ [ Bt ] dV
T

V V (3.36)
= éë K p ùû + [ K t ]

where
CHAPTER 3 67

éë K p ùû = ò éë B p ùû ⋅ [d n ] ⋅ éë D p ùû ⋅ [d n ] ⋅ éë B p ùû dV
T T

= ò éë B p ùû ⋅
A
T
( ò [d ]
n
T
)
⋅ éë D p ùû ⋅ [d n ] dz ⋅ éë B p ùû dA (3.37)

T
 p ù ⋅ é B p ù dA
= ò éë B p ùû ⋅ éë D û ë û
A

and

[ Kt ] = ò [ Bt ]T ⋅ [ Dt ] ⋅ [ Bt ] dV
V
(3.38)
= ò [ Bt ] ⋅
A
T
( ò [ D ] dz ) ⋅ [ B ] dA = ò [ B ]
t t
A
t
T
t ù ⋅ [ Bt ] dA
⋅ éë D û

In Eqs. 3.36-3.38 éë D p ùû is the bending and [ Dt ] is the shear component of the material
matrix, respectively.
The integration through the depth of the slab, which is placed in parentheses in
Eqs. 3.37 and 3.38, is performed by subdividing the slab into imaginary steel and concrete
layers, as shown in Fig. 3.3. The transverse displacement field is assumed to be continuous in
each slab element so that there are no gaps between layers. Each layer can have different
material properties, which are, however, assumed to be constant through the thickness of the
layer. The integrals through the depth of the slab can be evaluated by summation over all
layers
 ù = [d ] ⋅ é D ù ⋅ [d ] dz
éë D ò n ë pû n
T
pû

nc z j +1 ns
=å ò z ⋅ éë D pc ùû dz + å zi2 ⋅ [ Ds ]i ti
2
(3.39)
j
j =1 z j i =1

1 nc 3 ns
= å ( z j +1 − z j ) ⋅ éë D pc ùû + å zi2 ⋅ [ Ds ]i ti
3

3 j =1 j
i =1

and
nc
t ù = [ Dt ] ⋅ dz = å ( z j +1 − z j ) ⋅ [ Dtc ]
éë D û ò j
(3.40)
j =1

where z = z − ( d nx + d ny ) 2 for simplicity and nc and ns is the total number of concrete and
steel layers, respectively. éë D pc ùû j is the plane stress material matrix of the j-th concrete layer,
[ Dtc ] j is the transverse shear material matrix of the j-th concrete layer and [ Ds ]i is the
68 CHAPTER 3

material matrix of the i-th steel layer whose thickness is established from Eq. 2.50 based on
the distributed steel model. The reinforcing steel only contributes to the in-plane stiffness of
the slab, since the dowel action of the reinforcement is neglected in this study.
Once the node displacements of the finite element model are determined, Eqs. 3.31
and 3.33 yield the in-plane and the transverse shear strains of each layer, respectively. These
relations are summarized below

ì é ù s+ ü
F xj = ê d nx − ( z j +1 + z j ) ú ⋅ x
1
ï ï
ï ë 2 û sx ï
ï é ù s+ ï
F yj = ê d ny − ( z j +1 + z j ) ú ⋅ y
1
ï ï
ï ë 2 û sy ï
ï s+ y öï
{F c } j = í H xyj = é d − 1 ( z j +1 + z j ) ù ⋅ çæ s+ x + ÷ý (3.41)
êë 2 2 úû sy sx
ï è øï
ï sw ï
ï H xzj = − +x ï
ï sx ï
ï sw ï
ï H yzj = − +y ï
îï sy þï
ì s+ x ü
ïï ( d − z ) ⋅
sx ïï
nx i

{F s }i =í
s+ ý
(3.42)
ï( d ny − zi ) ⋅ y ï
ïî sy ïþ
The stresses of each layer are determined from the corresponding material stress-
strain relation

{6c } j = [ Dc ] j ⋅ {. c } j (3.43)

{6s }i = [ Ds ]i ⋅ {. s }i (3.44)
where {6s }i and {. s }i denote the stress and strain at mid-depth zi of the i-th steel layer, and
{6c } j and {. c } j denote the stress and strain at mid-depth of the j-th concrete layer,

respectively.
In the nonlinear analysis algorithm used in this study, the nonlinear constitutive
relations are satisfied by iterative successive corrections. In each iteration the stresses are
determined from the corresponding strain increments at the end of the previous iteration step,
while the material properties are only updated at the beginning of the next load step.
CHAPTER 3 69

Since in a nonlinear problem the stresses which are determined during the iterative
phase of the algorithm are not, generally, in equilibrium with the applied loads, unbalanced
nodal forces result. These are corrected during the subsequent iteration until a specified
tolerance is satisfied. The unbalanced nodal forces are the difference between applied and
resisting or equivalent forces

{R}unbalanced = {R}applied − {R}equivalent (3.45)

The resisting nodal forces are statically equivalent to the stress field which results
from the current deformation state of the finite element model. The resisting nodal forces are
determined from

{R}equivalent = ò [ B ] ⋅{6}dV = ò éë Bp ùû ⋅[d n ] ⋅ {6 p } dV + ò [ Bt ] ⋅{6t }dV

T T T T

V V V

A
T
= ò éë B p ùû ⋅ ( ò [d ] ⋅ {6 }dz ) dA + ò [ B ] ⋅ ( ò {6 }dz ) dA
n
T
p
A
t
T
p (3.46)

= ò éë B p ùû ⋅ {6
 p } dA + ò [ Bt ] ⋅ {6
 t }dA
T T

A A

and the integration through the depth of the slab is performed by summation over all
imaginary concrete and steel layers

ì nc ns ü
 x = å é ( z j +1 − d nx ) − ( z j − d nx ) ù ⋅ 6 x + å ( zi − d nx ) ⋅ 6 x ti ï
ï6 1 2 2

ï j =1 2
ëê ûú i =1 ï
ï ï
ï ï
nc ns
 y = å ( z j +1 − d ny ) − ( z j − d ny ) ⋅ 6 y + å ( zi − d ny ) ⋅ 6 y ti ý
1é 2 2
ù
 p ùû = í6
éë6 (3.47)
j =1 2
ê
ë ú
û
ï i =1 ï
ï ï
1 éæ dö ù
2 2
nc
dö æ
ï7 =å ê z − − z − ú⋅7 ï
ï xy j =1 2 ëê çè j +1 2 ÷ø çè j 2 ÷ø ûú xy ï
î þ
ì nc
ü
ï zx å ( z j +1 − z j ) ⋅ 7 zx ï

7 =
{6 t } = ïí ï
j =1
nc ý (3.48)
ï7 = ( z j +1 − z j ) ⋅ 7 zy ï
ïî zy åj =1
ïþ

Details of the numerical implementation along with a summary of the algorithm and a
discussion of the convergence criterion are presented in the following section.
70 CHAPTER 3

3.4 Numerical Implementation

3.4.1 Iteration Method

The numerical implementation of the finite element model requires the solution of Eq.
3.19. This is a system of simultaneous nonlinear equations, since the stiffness matrix [ K ] , in
general, depends on the displacement vector {d } . The solution of this system of nonlinear
equations is typically accomplished with an iterative method. The load vector {R} is
subdivided into a number of sufficiently small load increments, which are successively
applied (Fig. 3.4).

F4

F3

F2

Actual

F1 Numerical

0 r

FIGURE 3.4 INCREMENTAL LOAD METHOD WITHOUT CORRECTION

At each load step a linear approximation of the stiffness matrix [ K ] is established and
the resulting system of linear equilibrium equations is solved for the displacement increments
which correspond to the applied load increments. Since the stiffness matrix [ K ] changes
under these displacement increments, the resisting forces of the structure do not equilibrate
the applied loads and unbalanced loads result (Eq. 3.45). In the subsequent correction phase
the displacement increments are iteratively improved, until a specified convergence criterion
is satisfied. If no correction phase is included in the nonlinear analysis algorithm, the
CHAPTER 3 71

numerical error grows from one load step to the next and the numerical solution drifts away
from the actual response, as shown in Fig. 3.4.
Depending on how the stiffness
matrix [ K ] is updated during the
correction phase, the iterative method
F
can be classified into three broad
categories even though variations of
these schemes are also possible: the
initial or constant stiffness method, the
tangent stiffness method and the secant
stiffness method (Taylor and Hinton
0 r
(a) 1980).
The tangent stiffness method
F
requires the smallest number of iterations
to arrive at the solution, but has the
disadvantage that the stiffness matrix
[ K ] needs to be reformed and
triangularized at each iteration. The
0 r initial stiffness method, on the other
(b)
extreme, requires the largest number of
F
iterations, but the stiffness matrix [ K ] is
only formed and triangularized once at
the beginning of the load step. Generally,
the tangent stiffness method is more
efficient than the initial stiffness method
for models with a small number of
0 r degrees of freedom, even though the
(c)
selection of one method over another
FIGURE 3.5 ITERATIVE METHOD
often depends on numerical stability
(A) INITIAL STIFFNESS METHOD considerations. The secant method, for
(B) TANGENT STIFFNESS METHOD example, is well suited for structures
(C) SECANT STIFFNESS METHOD which exhibit softening response.
Even though it is possible to use
72 CHAPTER 3

any stiffness matrix in the first iteration of a new load increment, which constitutes the
advancing phase of the solution algorithm, it is quite common to use the tangent stiffness
matrix for this purpose. The nonlinear solution scheme selected in this study uses the tangent
stiffness matrix at the beginning of the load step in combination with a constant stiffness
matrix during the subsequent correction phase, as shown in Fig. 3.6.

3.4.2 Solution Algorithm

Every nonlinear analysis algorithm consists of four basic steps: the formation of the
current stiffness matrix, the solution of the equilibrium equations for the displacement
increments, the state determination of all elements in the model and the convergence check.
These steps are presented in some detail in the flow diagram of Fig. 3.7 for the plane stress
and the plate bending problem.
A summary of the steps of the nonlinear solution algorithm along with the
corresponding operations is presented below and is also identified in Fig. 3.7.
1. Form the tangent stiffness matrix [ K ] at the current displacement state of the finite
element model.
(a) In the case of RC beams use Eqs. 3.28 and 3.29 and 2.34
(b) In the case of RC slabs use Eqs. 3.36-3.38
2. Given the current load increment {R} calculate the node displacement increments
from Eq. 3.19, which is recast in incremental form

[ K ] ⋅ {d } + {F }p + {F }g + {F }. + {F }6 − {R} = 0
o o

and determine the current displacement increments {d } and displacements {d } .

3. Apply the strain-displacement relation to determine the concrete strains.

(a) In the case of RC beams use Eq. 3.27.

{. } = [ B ] ⋅ {d }
p
ele

(b) In the case of RC slabs use Eq. 3.41

CHAPTER 3 73

ì é ù s+ ü
F xj = ê d nx − ( z j +1 + z j ) ú ⋅ x
1
ï ï
ï ë 2 û sx ï
ï é ù s+ ï
F yj = ê d ny − ( z j +1 + z j ) ú ⋅ y
1
ï ï
ï ë 2 û sy ï
ï s+ y öï
{F c } j = í H xyj = é d − 1 ( z j +1 + z j ) ù ⋅ çæ s+ x + ÷ý
êë 2 2 úû sy sx
ï è øï
ï sw ï
ï H xzj = − +x ï
ï sx ï
ï sw ï
ï H yzj = − +y ï
îï sy þï
4. Determine the new crack direction (Eq. 2.27)
H xy
tan 2+ =
Fx − F y

5. Calculate the principal concrete strains (Eq. 2.28)

ì F1 ü é cos2 + sin 2 + sin + cos + ù ì T1 ü

í ý=ê 2 ú⋅í ý
î F 2 þ ë sin + cos + − sin + cos + û îT 2 þ
2

6. Use the concrete stress-strain relation to calculate the stresses in the global coordinate

F4

F3

F2
Actual
Numerical
F1

0 r

FIGURE 3.6 MODIFIED NEWTON-RAPHSON ALGORITHM

74 CHAPTER 3

system
(a) In the case of RC beams use the following equation.

ì 6x ü ì .x ü
ï ï ï ï
í 6 y ý = [ DGL ]c ⋅ í . y ý
ï7 ï ï0 ï
î xy þ î xy þ
(b) In the case of RC slabs use Eq. 3.43.

{6c } j = [ Dc ] j ⋅ {. c } j
7. Determine the principal stresses from the stresses in the global coordinate system using
Eq. 2.29.

ì Tx ü
ì T1 ü é cos2 + sin 2 + 2 sin + cos + ù ï ï
í ý=ê 2 ú ⋅ íTy ý
îT 2 þ ë sin + cos + −2 sin + cos + û ïU ï
2

î xy þ
8. Use the principal stresses to locate the present state of concrete in the biaxial principal
stress space of Fig. 2.2 and, if required, determine new principal stress values from the
principal strain values of Step 5. The process of state determination is different in each
region of the principal stress space.
(a) In the compression-compression region:
(1) Unless concrete remains inside the initial yield surface, in which case the stresses
are those of Steps 6 and 7, determine the principal stress ratio α from Eq. 2.3.
61
*= where 61 ≥ 62
62
(2) The ultimate strengths 61 p , 6 2 p and the corresponding strains .1 p , . 2 p of the
stress-strain relations in the principal strain directions are recalculated with Eqs.
2.4-2.7.
1 + 3.65*
62 p = ⋅ fc
(1 + * ) 2

æ 6 ö
. 2 p = . co ⋅ ç 3 2 p − 2 ÷
è fc ø
CHAPTER 3 75

61 p = * ⋅ 6 2 p

é æ 61 p ö
3
æ 61 p ö
2
æ 61 p ö ù
.1 p = .co ⋅ ê −1.6 ⋅ ç ÷ + 2.25 ⋅ ç ÷ + 0.35 ⋅ ç ÷ú
êë è fc ø è fc ø è f c ø úû

(3) The principal stresses are determined from the corresponding principal strains
using the modified uniaxial stress-strain relation in Fig. 2.4. If the principal strain
exceeds the crushing strain, the corresponding stress is set equal to zero.
(b) In the compression-tension and the tension-tension region:
(1) The principal stresses can be directly determined from the uniaxial stress-strain
relation in Fig. 2.4. It is assumed that the presence of a tensile stress does not
modify the uniaxial stress-strain relation in compression. The maximum tensile
strength f eq is reduced in accordance with the stress failure envelope in Fig. 2.2.

9. Determine the new stresses in the global coordinate system from the principal stresses.

ì T x ü é cos 2 + sin 2 + ù
ï ï ê ú ì T1 ü
í T y ý = ê sin + cos2 + ú ⋅ í ý
2

ïU ï êsin + cos + − sin + cos + ú îT 2 þ

î xy þ ë û
10. Determine the equivalent nodal forces.
11. For RC slabs repeat Steps 3-10 for each concrete layer.
12. Repeat steps 3-11 for each integration point.
13. The calculations for the reinforcing steel depend on the type of reinforcing steel model
used in the analysis.
(a) In the analysis of RC slabs with the distributed steel model:
(1) Determine the steel strains from the strain-displacement relation with Eq. 3.42

ì s+ x ü
ïï ( d nx − zi ) ⋅ sx ïï
{F s }i =í
s+ ý
ï( d ny − zi ) ⋅ y ï
îï sy þï
(2) Determine the new stress for the distributed steel model (Eq. 2.51).
76 CHAPTER 3

START

Read and transform input

Update and assemble the element Tangent

stiffness matrices in the total Step 1 stiffness
stiffness matrix of the structure
Iteration
method
Assemble the incremental load vector Step 2

Solve equilibrium equations.

Add displacement increments to the Step 2
current displacements
For all slab layers
For all elements

Determine the element strains

Step 3

Initial stiffness
for the current displacements
For all load increments

Use the material constitutive model

to determine the internal stresses Steps 4-13
and the internal element nodal forces

Determine the neutral axis location

Step 14
(this step is for RC slabs only)

Determine the difference between

external and internal nodal forces Steps 16 and 17

no
Convergence

yes
Print displacements, forces and stresses

no All load
steps
yes

STOP

FIGURE 3.7 OUTLINE OF SOLUTION ALGORITHM

CHAPTER 3 77

(3) Transform the local steel stress to the global coordinate system.
(4) Determine the equivalent nodal forces.
(5) Repeat Steps (1)-(4) for each steel layer.
(6) Repeat steps (1)-(5) for each integration point.
(b) In the analysis of RC beams with a reinforcing steel model with perfect bond:
(1) Determine the displacement increments at the ends of the reinforcing bar element
(Eq. 3.49).

{d s } = [T1 ] ⋅ [T2 ] ⋅ {d }

ele

(2) Calculate the steel forces according to the discrete reinforcing steel model with
perfect bond with Eq. 2.30.
(3) Transform the local steel forces to the global coordinate system to obtain the
contribution of the steel model to the equivalent nodal forces.
(4) Repeat steps (1)-(3) for each embedded reinforcing steel element.
(c) In the analysis of RC beams with a reinforcing steel model with bond-slip a
special procedure is implemented which is described in detail in the following
section. This procedure determines the contribution of steel and bond to the
equivalent nodal forces at the concrete degrees of freedom from Eq. 2.66.

{P } = {P } − [ K ] ⋅ [ K ] ⋅ {P }
c
*
c cs ss
−1
s

This process is repeated for each embedded reinforcing steel element.

14. For RC slabs only calculate the new neutral axis location using the conditions:

ò 6 dz = ò 6 dz = ò 7
x y xy dz = 0

Repeat for each integration point. The new neutral axis location is only used in the next
load increment.
15. Assemble the equivalent nodal forces.
16. Repeat Steps 3-15 for all elements.
17. Calculate the unbalanced nodal forces from Eq. 3.45
78 CHAPTER 3

{R}unbalanced = {R}applied − {R}equivalent

18. Repeat Steps 2-17 until convergence occurs.
19. Increment the load and go back to Step 1.

3.4.3 Solution Algorithm for Reinforcing Bar with Bond-Slip

In Step 2 of the solution algorithm the concrete displacement increments at the nodes
of the finite element mesh are determined. If a reinforcing steel model with perfect bond is
used, the displacement increments at the ends of the bar element are equal to the
corresponding displacement increments of the quadrilateral concrete element, in which the
bar element is embedded. In this case the displacement increments at the end of the
reinforcing bar element can be determined from

{d s } = [T1 ] ⋅ [T2 ] ⋅ {d }

ele
(3.49)
where the transformation matrices [T1 ] and [T2 ] are given by Eqs. 2.33 and 2.45 or 2.49,
respectively. {d } is the vector of concrete displacement increments and {d s } is the
ele

vector of displacement increments at the ends of the reinforcing element parallel to the axis
of the bar. If the reinforcing steel model with bond-slip is used in the finite element analysis,
then the relation between displacement increments at the ends of the reinforcing bar element
and the corresponding concrete displacement increments is provided by Eq. 2.64. Since the
terms of the stiffness matrix in Eq. 2.64 which correspond to the steel degrees of freedom are
removed by static condensation before assembly into the structure stiffness matrix, the
determination of the steel displacements and the corresponding forces now requires an
additional step.
Once the displacement increments at the nodes of the concrete finite elements {d }
ele

are determined for the current load increment, they can be transformed by matrices [T1 ] and
[T2 ] to yield the concrete displacement increments {d c } at the ends of the steel element and
parallel to the axis of the bar

{d c } = [T1 ] ⋅ [T2 ] ⋅ {d }

ele
(3.50)

Using the second row of the matrix relation in Eq. 2.64, which is rewritten in terms of
displacement increments
CHAPTER 3 79

{Ps } = [ Kcs ] ⋅ {d c } + [ K ss ] ⋅ {d s } (3.51)

and solving for {d s } yields

[ K ss ] ⋅ {d s } = ({Ps } − [ K cs ] ⋅ {d c }) (3.52)

Eq. 3.52 expresses the condition of equilibrium of the reinforcing bar element. The
steel force and deformation increments can now be determined by assembling the steel
element matrices on the left hand side and the forces on the right hand side of the equation
and imposing appropriate boundary conditions at the ends of the entire reinforcing bar. It
should be noted that the only forces acting on the reinforcing bar are the bond (transfer)
forces and any applied end forces, so that upon completion of the assembly process only the
term − [ K cs ] ⋅ {d c } and one concentrated force at each end remain on the right hand side of
Eq. 3.52. Since the concrete displacement increments {d c } are known from Eq. 3.50, this
system of linear equations can now be solved for the unknown steel displacement increments
{d s } . Once the steel displacement increments are determined, the state determination of the
steel and bond elements is performed and the corresponding forces are determined. These are
substituted into Eq. 2.66 to yield the contribution of the reinforcing steel model to the
equivalent nodal forces at the global degrees of freedom. These forces are assembled for all
steel and concrete elements (Step 16) and subtracted from the applied load increments
(Step 16) to yield the unbalanced forces. The process continues until the convergence
criterion is satisfied, as described previously. It should be noted that the bond and steel
stiffness is not updated at every iteration, but only at the beginning of a new load increment,
in accordance with the initial stiffness method which is used in this study. Consequently, the
equilibrium of the reinforcing bar and the calculation of the equivalent nodal forces requires
that the stiffness matrix of the entire reinforcing bar be triangularized only once at the
beginning of a new load step.
An alternative solution to the direct stiffness method is the transfer matrix method.
This method was already used by Ciampi et al. (1982), Filippou et al. (1983) and Filippou
(1986) in the solution of the hysteretic behavior of anchored reinforcing bars. In the transfer
matrix method the solution process starts at one end of the reinforcing bar where one
boundary condition is known. The second condition is assumed, so that the solution can
advance from one element to the next until reaching the far end of the reinforcing bar. The
boundary condition which is known at the far end of the reinforcing bar provides the equation
for determining the assumed boundary condition at the starting end. The implementation of
80 CHAPTER 3

the transfer matrix method in the solution of the equilibrium of the embedded reinforcing bar
with bond-slip offers certain computational advantages over the direct stiffness method, as
will be discussed in more detail later. Details of the implementation are presented in the
following.
The reinforcing bar is subdivided into n-elements, as shown in Fig. 3.8. Eq. 2.64
represents the equilibrium equations of a steel element. Focusing attention on the k-th
element Eq. 3.51 yields

{Ps } = [ K cs ]k ⋅ {d c } + [ K ss ]k ⋅ {d s }

k k k
(3.53)

with reference to Eq. 2.63 in incremental form Eq. 3.53 can be explicitly written (Fig. 3.8)

ìP2 ü é k s + kbi − k s ù ìd 2 ü ì kbi ⋅ d1 ü

k k k k

í ý =ê ⋅í ý −í ý
k s + kbj ûú îd 4 þ îkbj ⋅ d 3 þ
(3.54)
îP4 þ ë − k s
where it should be noted that the concrete displacement increments d1 and d 3 are known
from Eq. 3.50 and Step 2 of the global solution algorithm. Solving Eq. 3.54 for the force and
displacement increment at node 4 yields

1st ...... (k-1)-th k-th (k+1)-th ...... n-th

(a)

2 k-th 4 2 (k+1)-th 4

1 3 1 3
(b) (c)

Imaginary node

Element

FIGURE 3.8 MODELING OF REINFORCING BAR FOR TRANSFER MATRIX METHOD

CHAPTER 3 81

k s + kbj
k
é ù
1 é − ( k s + kbi ) k s ⋅ ( kbi + kbj ) + kbi ⋅ kbj ù
k
ì P2 ü ê kbi ⋅
k k k
ì P4 ü kbj ú ì d1 ü
í ý = ê ú ⋅í ý − ks ⋅í ý (3.55)
îd 4 þ k s êë −1 k s + kbi úû îd 2 þ êê ú îd 3 þ
ë kbi 0 úû

which is written in compact form

k k k
ì P4 ü k ì P2 ü k ì  d1 ü
í ý = [Q ] ⋅ í ý − [ R] ⋅ í ý (3.56)
îd 4 þ îd 2 þ îd 3 þ
For the transition from the k-th to the (k+1)-th element the following relations hold (Fig. 3.8)
k +1 k k k
ì P2 ü ì −P4 ü é −1 0ù ì P4 ü ì P4 ü
í ý =í ý =ê ú ⋅ í ý = [ S ] ⋅ í ý (3.57)
îd 2 þ î d 4 þ ë 0 0 û îd 4 þ îd 4 þ
The substitution of Eq. 3.56 into Eq. 3.57 yields
k +1 k k
ì P2 ü ì P ü k ì d1 ü
= [ S ] ⋅ [Q ] ⋅ í 2 ý − [ S ] ⋅ [ R ] ⋅ í
k
í ý ý (3.58)
îd 2 þ îd 2 þ îd 3 þ
After realizing that the multiplication by matrix [ S ] amounts to changing the signs of the
first column of matrices [Q ] and [ R ] Eq. 3.58 is simplified to
k +1 k k
ì P2 ü k ì P ü k ì d1 ü
í ý = éëQ ùû ⋅ í 2 ý − éë R ùû ⋅ í ý (3.59)
îd 2 þ îd 2 þ îd 3 þ
where the change of sign of the first column of [Q ] and [ R ] results in éëQ ùû and éë R ùû ,
respectively.
Eq. 3.59 relates the force and displacement increment at the beginning of steel element k+1
with the force and displacement increment at the beginning of steel element k. By applying
Eq. 3.58 successively to elements k-1, k-2, ..., 2, 1 and combining the results the following
transfer matrix relation results
82 CHAPTER 3

k +1 k k
ì P2 ü k ì P ü k ì d ü
í ý = éëQ ùû ⋅ í 2 ý − éë R ùû ⋅ í 1 ý
îd 2 þ îd 2 þ îd 3 þ
k −1 k −1 k
k k −1 ì P2 ü k k −1 ì d1 ü k ì d1 ü
= éëQ ùû ⋅ éëQ ùû ⋅í ý − éëQ ùû ⋅ éë R ùû ⋅í ý − éë R ùû ⋅ í ý
îd 2 þ îd 3 þ îd 3 þ
1 1 k
k k −1 1 ì P ü k k −1 2 1 ì  d1 ü k ì  d1 ü
= éëQ ùû ⋅ éëQ ùû K éëQ ùû ⋅ í 2 ý − éëQ ùû ⋅ éëQ ùû K éëQ ùû ⋅ éë R ùû ⋅ í ý − K − éë R ùû ⋅ í ý
îd 2 þ î d 3 þ îd 3 þ
(3.60)
After replacing k+1 with n in Eq. 3.60 and applying Eq. 3.56 for element n the following
relation between the forces and displacements at the two ends of the reinforcing bar results
n 1 1 k
ì P4 ü 1 ì P2 ü 1 ì  d1 ü k ì d1 ü
ý = [Q ] ⋅ éëQ ùû K éëQ ùû ⋅ í ý − [Q ] ⋅ éëQ ùû K éëQ ùû ⋅ éë R ùû ⋅ í
n −1 n −1 2
ý − K − éë R ùû ⋅ í
n n
í ý
îd 4 þ îd 2 þ îd 3 þ îd 3 þ
(3.61)
In Eq. 3.61 the concrete displacement increments are known from Eq. 3.50.
One boundary condition is known at each end of the reinforcing bar. After assuming
the second boundary condition at the starting end, Eq. 3.61 yields the force and displacement
increment at the other end of the reinforcing bar. Since one condition is known at that end,
the initial assumption needs to be corrected until the known boundary condition at the far end
is satisfied.
It should be noted that the transfer matrix method is applied during the correction
phase of the global solution algorithm. Consequently, the transfer matrix relations in Eq. 3.61
are linear, since no updating of the bond and steel stiffness takes place during the correction
phase which is based on the initial stiffness method. Thus, Eq. 3.61 can be solved very
rapidly through a series of multiplications and no iterations are required. An advantage of the
transfer matrix method is its numerical stability in connection with the cyclic bond-slip
model of Section 2.4.1. This will be discussed further in the correlation studies of Chapter 4.
The satisfaction of equilibrium of the reinforcing bar by the transfer matrix method of
Eq. 3.61 yields the steel displacement increments. Since the concrete displacement
increments at the ends of the bar are also known from Eq. 3.50, the relative displacements
can be readily determined. The state determination of the steel and bond elements can now be
undertaken yielding the new steel and bond forces and the updated stiffness matrices. The
latter are only needed at the beginning of a new load step, when the stiffness matrix of the
structure is updated. The substitution of the new steel and bond forces in Eq. 2.66 yields the
CHAPTER 3 83

equivalent nodal forces at the global degrees of freedom. These forces are subtracted from the
applied load increments to yield the unbalanced forces. The process continues until the
convergence criterion is satisfied.

3.4.4 Convergence Criterion

The criterion for measuring the convergence of the iterative solution is based on the
accuracy of satisfying the global equilibrium equations or on the accuracy of determining the
total displacements. The accuracy of satisfying the global equilibrium equations is controlled
by the magnitude of the unbalanced nodal forces. The accuracy of the node displacements
depends on the magnitude of the additional displacement increment after each iteration. The
latter convergence criterion is used in this study. This can be expressed as
1
é i 2ù
( )
2

êå  d j ú
Ed = ë j û ≤ TOLER (3.62)
1
é i 2ù
êå ( d j ) ú
2

ë j û
where the summation extends over all degrees of freedom j, d j is the displacement of degree
of freedom j, d ij is the corresponding increment after iteration i and TOLER is the
specified tolerance.
In the nonlinear analysis of RC structures the load step size must be small enough so
that unrealistic "numerical cracking" does not take place. These spurious cracks can
artificially alter the load transfer path within the structure and result in incorrect modes of
failure. Crisfield (1982) has shown that such numerical disturbance of the load transfer path
after initiation of cracking can give rise to alternative equilibrium states and, hence, lead to
false ultimate strength predictions. In order to avoid such problems after crack initiation the
load is increased in steps of 2.5-5.0% of the ultimate load of the member.
The failure load is assumed to occur at a load level for which a large number of
iterations are required for convergence. This means that very large strain increments take
place during this step and that equilibrium cannot be satisfied under the applied loads.
Obviously, the maximum number of iterations depends on the problem and the specified
tolerance, but a maximum of 30 iterations seems adequate for a tolerance of 1%. This is the
limit in the number of iterations selected in this study.
 CHAPTER 4
APPLICATIONS

4.1 Introduction

A number of correlation studies are conducted in this chapter with the objective of
establishing the ability of the proposed model to simulate the response of reinforced concrete
beams, slabs and beam-column joint subassemblages. In order to independently test the
reinforcing steel model with bond-slip, the response of anchored reinforcing bars under
monotonic pull-out and cyclic push-pull loads is also studied. Even though cyclic loading is
outside the scope of this study, the cyclic behavior of reinforcing bars is presented in order to
demonstrate the capabilities of the proposed numerical solution scheme and in order to
highlight the pronounced effect of cyclic bond-deterioration on the hysteretic response of RC
members, which is the subject of a follow-up study.

4.2 Anchored Reinforcing Bar Under Monotonic and Cyclic Loads

Viwathanatepa, Popov and Bertero (1979b) tested several anchored reinforcing bars
simulating anchorage and loading conditions in interior beam-column joints of moment
resisting frames which are subjected to a combination of gravity and high lateral loads. #6, #8
and #10 reinforcing bars were anchored in well confined concrete blocks and were subjected
to monotonic pull-out at one end, monotonic pull-out at one end with simultaneous push-in at
the other (called push-pull) and cyclic push-pull.

Two specimens are selected for comparison with the proposed reinforcing steel model
with bond-slip. The first specimen is an anchored #8 bar in a well confined block of 25 in.
width, which corresponds to an anchorage length of 25 bar diameters. This specimen was
subjected to a monotonic pull-out under displacement control at one end only. The second
specimen also involves a #8 reinforcing bar with identical dimensions which was subjected to
a cyclic push-pull loading with gradually increasing end slip value. Both specimens have
been the subject of earlier analytical correlation studies by Viwathanatepa et al. (1979b),

85
86 CHAPTER 4

Ciampi et al. (1982), Yankelevky (1985) and Filippou (1986). The material properties of
concrete and reinforcing steel, as reported in Table 2.1 of the report by Viwathanatepa et al.
(1979b), are as follows: the concrete cylinder strength is 4,700 psi for the specimen under
monotonic pull-out and 4,740 psi for the specimen under cyclic push-pull. The yield strength
of the reinforcing steel is 68 ksi and the yield strain is 0.23%.
100

80
STEEL STRESS [ksi]

60

40
25 in

20

0
0.00 0.02 0.04 0.06 0.08 0.10

SLIP [in]
FIGURE 4.1 STRESS-SLIP RESPONSE OF ANCHORED REINFORCING BAR
UNDER MONOTONIC PULL-OUT LOADING

In the monotonic pull-out test the bilinear steel model in Fig. 2.14 of Section 2.3 is
used. The strain hardening modulus is equal to 411 ksi, which corresponds to 1.4% of
Young's modulus. In the cyclic push-pull test the model of Menegotto-Pinto, as modified by
Filippou et al. (1983), is used. In both specimens the refined bond-slip model in Fig. 2.19 of
Section 2.4 is used. The parameters of the model are equal to those used by Filippou et al.
(1983) and Filippou (1986) in earlier investigations, namely: u1 = 0. 02756 in ,
u2 = 0. 07874 in, u3 = 0. 2756 in , τ1 = 2350 psi and τ 3 = 870 psi (Fig. 2.19). In the study by
Ciampi et al. (1982) the bond-slip relation was modified in the outer unconfined portions of
the anchorage length. By contrast, in the present study the same bond stress-slip relation is
used along the entire anchorage length for the sake of simplicity, since the hysteretic behavior
is not the focus of the present study. Under cyclic loading conditions this assumption leads to
underestimation of the bond resistance at the push-in end of the reinforcing bar. 25 steel
elements of 1 in length each were used in modeling the anchored reinforcing bar.
CHAPTER 4 87

20 40

Experiment Experiment
Present study Present study
16 Yankelevsky
32 Yankelevsky
STEEL STRESS (ksi)

Viwathanatepa et al. Viwathanatepa et al.

STEEL STRESS (ksi)

12 24

8 16
25 in 25 in

4 8

0 0
0 5 10 15 20 25 0 5 10 15 20 25

DISTANCE (in) DISTANCE (in)

(a) End Stress=20 ksi (b) End Stress=40 ksi
60 70
Experiment
Experiment
Present study
Present study
48 Yankelevsky 56
Viwathanatepa et al.
STEEL STRESS (ksi)
STEEL STRESS (ksi)

36 42

25 in 25 in
24 28

12 14

0 0
0 5 10 15 20 25 0 5 10 15 20 25
DISTANCE (in) DISTANCE (in)
(c) End Stress=60 ksi (d) End Stress=70 ksi

FIGURE 4.2 STRESS DISTRIBUTION ALONG ANCHORED REINFORCING BAR

Fig. 4.1 shows the relation between pull-out slip and steel stress at one end of the
reinforcing bar under monotonic pull-out loading. Figs. 4.2a-d show the distribution of steel
stresses along the anchorage length of the reinforcing bar at different loading stages. The
experimental results are compared with the analytical results of Viwathanatepa et al. (1979b),
Yankelevsky (1985) and those of the present study. The results of the present study show the
best agreement, particularly, with increasing load. It should be noted that the model of
Yankelevsky (1985) does not allow for yielding of the reinforcing steel. The results by
Viwathanatepa et al. (1979b) are from a linear finite element analysis, since no stress
distributions are presented for the nonlinear model proposed in that study.
88 CHAPTER 4

100

80

60
STEEL STRESS (ksi)

40 25
in
20

-20

-40

-60

-80

-100
-0.04 -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 0.04

SLIP (in)

(a) Load Cycles Before Yielding of Reinforcing Steel

100

80

60
STEEL STRESS (ksi)

40 25
in
20

-20

-40

-60

-80

-100
-0.10 -0.08 -0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0.10

SLIP (in)

(b) Load Cycles After Yielding of Reinforcing Steel

FIGURE 4.3 STRESS-SLIP RESPONSE OF ANCHORED REINFORCING BAR

UNDER CYCLIC PUSH-PULL. ANALYTICAL RESULTS.

CHAPTER 4 89

(a) Load Cycles Before Yielding of Reinforcing Steel

(b) Load Cycles After Yielding of Reinforcing Steel

FIGURE 4.4 STRESS-SLIP RESPONSE OF ANCHORED REINFORCING BAR

UNDER CYCLIC PUSH-PULL. EXPERIMENTAL RESULTS.

90 CHAPTER 4

12 24
Experiment
Experiment Present study
Yankelevsky
8 Present study 16 Viwathanatepa et al.
Yankelevsky

STEEL STRESS (ksi)

STEEL STRESS (ksi)

Viwathanatepa et al.
4 8

0 0

-4 -8

-8 25
-16 in
25 in

-12 -24
0 5 10 15 20 25 0 5 10 15 20 25
DISTANCE (in) DISTANCE (in)

(a) First Cycle (b) Fourth Cycle

40 70
Experiment
Present study
Yankelevsky
20 Viwathanatepa et al. 35 25 in

STEEL STRESS (ksi)

STEEL STRESS (ksi)

0 0

-20 -35

25 in Experiment
Present study

-40 -70
0 5 10 15 20 25 0 5 10 15 20 25
DISTANCE (in) DISTANCE (in)

(c) Seventh Cycle (d) First Post-Yield Cycle

10 24

1A 4C
STEEL STRESS (ksi)

5 12
STEEL STRESS (ksi)

1B 4F
0 0

1C 4I
-5 -12

25 in -24 25 in
-10
1A 24 4C
10
1B 4F
1 Experiment 4 Experiment
-10 1C Present study -24 4I Present study

0 5 10 15 20 25 0 5 10 15 20 25
DISTANCE (in) DISTANCE (in)

(e) Complete First Cycle (f) Complete Fourth Cycle

FIGURE 4.5 STRESS DISTRIBUTION ALONG ANCHORED REINFORCING BAR

UNDER CYCLIC PUSH-PULL LOADING

CHAPTER 4 91

Figs. 4.3a and b show the end stress-slip relation of the reinforcing bar under cyclic push-pull
loading. Fig. 4.3a shows the cycles before yielding of the reinforcement and Fig. 4.3b shows
the entire response. The corresponding experimental data are shown in Figs. 4.4a and b. The
steel stress distributions along the anchorage length of the reinforcing bar are shown in
Figs. 4.5a-f. The experimental results are compared with the analytical results of
Viwathanatepa et al. (1979b), Yankelevsky (1985) and those of the present study. The results
of the proposed model show excellent agreement in the early cycles in Figs. 4.5a-f. However,
the comparison of the overall response in Figs. 4.3 and 4.4 shows discrepancies between the
present model and the experimental data as the loading increases. Two factors contribute to
this: (a) the model is tested under load controlled conditions, while the specimen was
subjected to displacement controlled testing, and (b) the assumption that the bond stress-slip
relation is the same along the anchorage length of the bar is not reasonable for cycles which
induce significant bond damage.

EXAMPLES MATERIAL PROPERTIES

Ec Es fc fy ρ
[ksi] [ksi] [ksi] [ksi] [%]

J-4 3800 29500 4.82 44.9 0.99

BEAMS A-1 3367 31600 3.49 80.5 1.53

T1MA 3860 28200 4.60 46.0 0.62

BEAM-COLUMN JOINT 3900 28770 4.97 71.4 1.22

SLABS McNeice 4150 29000 5.50 40.0 0.85

B7 4000 30000 5.15 50.0 1.00

Note: The material properties of the reinforcing steel for beam A-1 and for the beam-
column joint are for #9 and #6 bars, respectively. Other material properties are as follows:
G f = 0. 5 lb / in, ft = 5 f c , ν c = 0.167 and ν s = 0. 333 .

TABLE 4.1 MATERIAL PROPERTIES USED IN APPLICATIONS

It is concluded that the proposed reinforcing steel model with bond-slip can describe
quite well the response of anchored reinforcing bars under monotonic and cyclic loading
conditions. It should be noted that the analytical results of the proposed model in Figs. 4.1-
4.5 are obtained with both solution methods of Chapter 2, namely, the direct solution method
and the transfer matrix method. The direct solution method is, in general, more economical,
92 CHAPTER 4

particularly, if advantage is taken of the tridiagonal stiffness matrix of the reinforcing bar
with bond-slip. The load step size required for convergence is, however, larger in the transfer
matrix method, so that both methods seem to require about the same amount of time for
determining the entire load-displacement response of Fig. 4.3b. Further comparison studies of
the computer time of these two methods are needed, before a definitive statement can be
made.

4.3 Reinforced Concrete Beams

Three simply supported reinforced concrete beams have been investigated. These
beams are specimen J-4 tested by Burns and Siess (1962), specimen A-1 tested by Bresler
and Scordelis (1963) and specimen T1MA tested by Gaston, Siess and Newmark (1972). In

0.5P
8in

18in
20in

6ft

FIGURE 4.6 CONFIGURATION OF BEAM J-4

these case studies the concrete was modeled by 8-node serendipity plane stress elements with
3x3 Gauss integration and the reinforcement was modeled by the 2-node truss element
described in Section 2.3. In all studies the bond-slip effect is taken into account with the bond
link element described in Section 2.4. The material properties of the three test specimens are
summarized in Table 4.1.
Specimen J-4 consists of a simply supported beam with a span of 12 ft (3.7 m) which
was subjected to a concentrated load at midspan. The geometry and the cross section of beam
J-4 are shown in Fig. 4.6 and the material properties are summarized in Table 4.1.
The finite element model in Fig. 4.7 represents only half of the structure taking
advantage of the symmetry in geometry and loading. The finite element discretization, the
arrangement of the reinforcement, the loading and support conditions are shown in Fig. 4.7.
CHAPTER 4 93

The parameters of the bond stress-slip relation in Fig. 2.17 are assumed as follows:
d s1 = 0.0007854 in , d s2 = 0.00589 in , τ1 = 286 psi and τ 2 = 572 psi. These values are derived
from the material properties in Table 4.1 based on the recommendations of Eligehausen et al.
(1983). The tension stiffening effect is taken into account by the model proposed in Section
2.2.3.3. The correlation between the measured load-displacement curve of the beam and the
analytical results which include the effects of tension stiffening and bond-slip is shown in
Fig. 4.8. The results are also compared with those of an earlier study by Barzegar and
Schnobrich (1986), who included only the tension stiffening effect by determining the slope
of the strain softening region according to Hillerborg's model (Hillerborg et al. 1976).

0.5P
NODES
ELEMENT NO.

4 8 12

3 7 11

2 6 10

1 5 9

FIGURE 4.7 FINITE ELEMENT IDEALIZATION OF BEAM J-4

In order to establish the contribution of the tension stiffening and the bond-slip effect
to the structural response a series of studies on beam J-4 were conducted. These studies also
aim at identifying the effect, if any, of finite element mesh size, tensile strength of concrete
and the cracked shear constant value. The results of these studies are discussed below. In
order to study the effect of finite element mesh size on the analytical results three different
mesh configurations were investigated, as shown in Fig. 4.9. In all configurations the number
of elements through the depth of the member remained the same and only the element size in
the span direction was progressively reduced.
The effect of finite element mesh size is studied for the three tension-stiffening
models discussed in Section 2.2.3.2. The results in Figs. 4.10a-c indicate that all tension-
94 CHAPTER 4

stiffening models exhibit satisfactory behavior and lead to response predictions which are
essentially independent from the finite element mesh size. The size dependence is more
pronounced in Hillerborg's model, while the crack band theory and the model proposed in
this study show a very satisfactory objectivity of the results. When comparing, however, the
results of the crack band theory and Hillerborg's model with those of the tension-stiffening
model proposed in this study, it is concluded that the load-deflection behavior of the former
models is stiffer than that of the present model, which as shown in Fig. 4.8, matches the
experimental behavior very well. This can be attributed to the assumption of those models
that the microcracks are uniformly distributed over the entire finite element.
36

30
TOTAL LOAD P (kips)

24

18

12

Experiment
6 Barzegar and Schnobrich
Present Study

0
0 0.1 0.2 0.3 0.4 0.5 0.6

CENTER DEFLECTION (in)

FIGURE 4.8 LOAD-DEFLECTION RELATIONSHIP OF BEAM J-4

If the tension stiffening effect is not included in the model, there is a marked
dependence of the analytical results on the finite element mesh size (Fig. 4.10d). The load-
deflection curves exhibit more flexible response with increasing grid refinement. Thus finite
element models which are based on the tensile strength of concrete and not on a fracture
energy criterion produce results which strongly depend on the analyst's choice of mesh
configuration and can lead to erroneous estimates of structural stiffness. This fact
corroborates earlier findings by several investigators (Bazant and Cedolin 1983; Reinhardt
1986).
CHAPTER 4 95

0.5P
16 ELEMENTS

0.5P
24 ELEMENTS

0.5P
48 ELEMENTS

FIGURE 4.9 FINITE ELEMENT MESH CONFIGURATIONS USED

The models of tension-stiffening and bond slip described in Chapter 2 reveal that
tension-stiffening has the opposite effect than bond-slip on the load-displacement response of
the finite element model of the structure. One way of including the effect of tension-
stiffening between cracks is to increase the stiffness of the reinforcing steel model. On the
other hand bond slip reduces the stiffness of the reinforcing steel model, as is clearly reflected
in the coefficient of the steel stiffness matrix in Eqs. 2.68 and 2.69. To identify the relative
contribution of each effect four different analyses are performed. In Fig. 4.11a the effect of
tension stiffening is included in both analytical results. The response which is depicted by the
long-dash line excludes the effect of bond-slip, while the response depicted with the dotted
line includes this effect. It is clear from the comparison of these results with the experimental
data shown by the solid line that the inclusion of both effects yields a very satisfactory
96 CHAPTER 4

agreement of the model with reality. Fig. 4.11b shows two analytical results which exclude
the effect of tension-stiffening. In this case the inclusion of bond-slip (dotted line) produces a
slightly more flexible response than the experiment, while the exclusion of bond-slip
produces a slighter stiffer response. The comparison of Fig. 4.11a with Fig. 4.11b illustrates
why some investigators have concluded that neither tension-stiffening nor bond-slip is
important in the analysis, provided that the right size of finite element mesh is used.
Fig. 4.12 shows the influence of the tensile strength of concrete on the response of the
beam specimen. In Fig. 4.12 the tensile strength of concrete, which is determined from the

36 36

30 30
TOTAL LOAD P (kips)

24 TOTAL LOAD P (kips) 24

18 18

12 12

16 ELEMENTS
16 ELEMENTS
6 24 ELEMENTS
24 ELEMENTS 6
48 ELEMENTS
48 ELEMENTS

0 0
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
CENTER DEFLECTION (in) CENTER DEFLECTION (in)
(a) Hillerborg's Model (b) Crack Band Model
36 36

30 30
TOTAL LOAD P (kips)

TOTAL LOAD P (kips)

24 24

18 18

12 12

16 ELEMENTS 16 ELEMENTS
6 24 ELEMENTS 6
24 ELEMENTS
48 ELEMENTS 48 ELEMENTS

0 0
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
CENTER DEFLECTION (in) CENTER DEFLECTION (in)

(c) Proposed Model (d) No Tension Stiffening

FIGURE 4.10 MESH SIZE EFFECT

 CHAPTER 4 97

36 36

30 30

TOTAL LOAD P (kips)

TOTAL LOAD P (kips)

24 24

18 18

12 12

Experiment Experiment
6 Perfect Bond 6 Perfect Bond
Bond-Slip Bond-Slip

0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.1 0.2 0.3 0.4 0.5 0.6

CENTER DEFLECTION (in) CENTER DEFLECTION (in)

(a) Tension Stiffening Effect (b) Bond-Slip Effect of Beam

FIG. 4.11 TENSION STIFFENING AND BOND-SLIP EFFECT FOR BEAM J-4

relation ft = α ⋅ f c , is varied between ft = 5 f c and ft = 7 f c . It is apparent that the

tensile strength only affects the load magnitude at initial cracking and the load magnitude at
the formation of a stable crack configuration. Since, the load magnitude at initial cracking is
relatively low, the effect of tensile strength on the overall response is negligible. Moreover, if
the tensile fracture strain ε o is adjusted to yield the same fracture energy according to
Eq. 2.14 of the proposed model, then the effect of tensile strength on the response is reduced
even further.

36 36

30 30
TOTAL LOAD P (kips)
TOTAL LOAD P (kips)

24 24

18 18

12 12
α = 5.0
α = 6.0
6
α = 7.0 6
λ =0.1=0.4=1.0

0 0
0 0.4 0.8 1.2 0 0.4 0.8 1.2
CENTER DEFLECTION (in) CENTER DEFLECTION (in)
Fig. 4.12 Effect of Concrete Tensile Strength Fig. 4.13 Effect of Cracked Shear Constant
98 CHAPTER 4

Fig. 4.13 shows that the effect of the cracked shear constant λ is negligible in the case
of slender members subjected to bending, since most of the member is in the stress state of
compression-compression or tension-tension and only a small portion of the structure near
the supports is in a stress state of tension-compression. It is only in the tension-compression
region of the biaxial stress space in Fig. 2.2 that the coupling between normal and shear
strains, which is represented by the cracked shear constant λ, plays an important role. While
this coupling effect is negligible in the case of slender beams in bending, it becomes
significant in the analysis of shear walls and deep beams.
The crack prediction and crack propagation in beam J-4 at two load stages is shown in
Fig. 4.14. It is evident that the response of this specimen is dominated by diagonal tension
cracking in good agreement with the experimental evidence.

7.2kips

14.4kips

FIGURE 4.14 CRACK PROPAGATION IN BEAM J-4

The next specimen which is investigated is beam A-1 tested by Bresler and Scordelis
(1963). The geometry and cross section dimensions of beam A-1 are presented in Fig. 4.15,
while Fig. 4.16 shows the finite element model of this structure. The material properties of
concrete and steel are summarized in Table 4.1. The bond stress-slip relation and any other
CHAPTER 4 99

material properties not specifically mentioned in Table 4.1 are the same as those used in the
analysis of beam J-4, since the concrete composition and strength of the two specimens is not
much different. Since Table 4.1 only contains the steel properties of the bottom longitudinal
steel of specimen A-1 (#9 bars), the following steel properties are assumed for the top steel
and the transverse reinforcement: Es = 29, 200 ksi and f y = 50.1 ksi for #4bars as top
longitudinal steel and Es = 27, 500 ksi and f y = 47. 2 ksi for #2 bars as transverse
reinforcement.

0.5P
305mm

2@#4
#2 TIES

2@63mm 553mm

1828.5mm 228mm

4@#9

FIGURE 4.15 CONFIGURATION OF BEAM A-1

0.5P CONCRETE NODES

STEEL ELEMENT NO.

5 10 15

4 9 14

3 8 13

2 7 12

1 6 11

STIRRUP

FIGURE 4.16 FINITE ELEMENT IDEALIZATION OF BEAM A-1

Fig. 4.17 compares the analytical results with the measured load-displacement
response of beam A-1. Very satisfactory agreement between analysis and experiment is
100 CHAPTER 4

observed. The analytical results by Adeghe and Collins (1986) are also shown in Fig. 4.17.
Since Adeghe and Collins did not account for the effect of bond slip in their study, the
comparison of analytical predictions demonstrates the significance of this effect in the load-
displacement response of a heavily reinforced specimen, such as beam A-1.

450

360
TOTAL LOAD P (KN)

270

180

Experiment
90 Adeghe and Collins
Present Study

0
0 3 6 9 12

CENTER DEFLECTION (mm)

FIGURE 4.17 LOAD-DEFLECTION RELATIONSHIP OF BEAM A-1

To identify the relative contribution of each effect four different analyses are also
performed for this specimen. In Fig. 4.18a the effect of tension stiffening is included in both
analytical results. The response which is depicted by the long-dash line excludes the effect of
bond-slip, while the response depicted with the dotted line includes this effect. It is clear
from the comparison of these results with the experimental data shown by the solid line that
the inclusion of both effects yields a very satisfactory agreement of the model with reality.
Fig. 4.18b shows two analytical results which exclude the effect of tension-stiffening. In this
case the inclusion of bond-slip (dotted line) produces a slightly more flexible response than
the experiment, while the exclusion of bond-slip produces a slighter stiffer response. Figs.
4.18a and b show that the contribution of bond-slip to the load-displacement response of the
specimen increases with the load. Near the ultimate strength of the beam the magnitude of the
CHAPTER 4 101

bond-slip contribution to the load-displacement response is almost twice that of the tension-
stiffening effect, which is of opposite sign. By contrast, in the lightly reinforced beam
specimen J-4 in Figs. 4.11a and b the magnitude of the bond-slip contribution to the load-
displacement response is comparable to that of the tension-stiffening effect, so that these two
effects practically cancel each other (long-dashed line in Fig. 4.11b).
450 450

360 360

TOTAL LOAD P (KN)

TOTAL LOAD P (KN)

270 270

180 180

Experiment Experiment
90 Perfect Bond 90 Perfect Bond
Bond-Slip Bond-Slip

0 0
0 3 6 9 12 0 3 6 9 12
CENTER DEFLECTION (mm) CENTER DEFLECTION (mm)

Fig. 4.18a Tension Stiffening and Fig. 4.18b Bond-Slip Effect without
Bond-Slip Effect in Beam A-1 Tension Stiffening in Beam A-1
The studies of beam specimens J-4 and A-1 show that the inclusion of, both, tension-
stiffening and bond-slip effect, yields excellent agreement of the analytical results with the
experiments. In lightly reinforced beams these two effects are of comparable importance and
can cancel each other at certain stages of the load history. This does not happen over the
entire response, since gradual bond damage leads to an increasing bond-slip contribution with
increasing load. The bond-slip effect is more marked in heavily reinforced beams, such as

0.5P 6in

3ft
10.72in
12in

4.5ft

FIGURE 4.19 CONFIGURATION OF BEAM T1MA

102 CHAPTER 4

specimen A-1, where it outweighs the effect of tension-stiffening.

0.5P NODES
ELEMENT NO.

4 8 12

3 7 11

2 6 10

1 5 9

FIGURE 4.20 FINITE ELEMENT IDEALIZATION OF BEAM T1MA

12

9
TOTAL LOAD P (kips)

3 Experiment
Bashur and Darwin
Present study

0
0 0.1 0.2 0.3 0.4 0.5

CENTER DEFLECTION (in)

FIGURE 4.21 LOAD-DEFLECTION RELATIONSHIP OF BEAM T1MA

CHAPTER 4 103

Finally, the predictions of the model are compared with the measured response of
specimen T1MA tested by Gaston et. al (1972). This specimen is also under-reinforced.
Fig. 4.19 shows the geometry, the dimensions of the cross section, the boundary conditions
and the loading arrangement of the specimen. Fig. 4.20 shows the finite element model of the
structure. As with the previous test specimens only half of the structure needs to be modeled,
because of symmetry in geometry and loading arrangement. The material properties are
summarized in Table 4.1. The bond stress-slip relation and any other material properties not
specifically mentioned in Table 4.1 are the same as those used in the analysis of beam J-4,
since the concrete composition and strength of the two specimens is not much different.
12 12

9 9
TOTAL LOAD P (kips)

TOTAL LOAD P (kips)

6 6

3 Experiment 3 Experiment
Perfect Bond Perfect Bond
Bond-Slip Bond-Slip

0 0
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
CENTER DEFLECTION (in) CENTER DEFLECTION (in)
Fig. 4.22a Bond-Slip and Tension Fig. 4.22b Bond-Slip Effect without
Stiffening Effect of Beam T1MA Tension Stiffening of Beam T1MA

Fig. 4.21 compares the analytical results with the measured load-deflection relation of
specimen T1MA. The results of a previous study by Bashur and Darwin (1978) are also
shown in Fig. 4.21. The results of the proposed model in Fig. 4.21 include the effects of
tension stiffening and bond slip. The initial discrepancy between analysis and experiment
stems from the fact that the specimen was probably extensively cracked before loading.
Otherwise the proposed model shows very satisfactory agreement with the measured
response.
Figs. 4.22a and b show the effect of tension stiffening and bond slip and the relative
contribution of each source of deformation to the load-displacement response of the
specimen. It is concluded that the contribution of tension stiffening and bond-slip is of
comparable magnitude but of opposite effect near the ultimate load of the specimen, in
104 CHAPTER 4

agreement with the results from specimen J-4. Thus the exclusion of both effects in Fig.
4.22b yields good overall agreement between analysis and experiment, even though the
analytical response is a little too flexible in the initial post-cracking stage and a little to stiff
near yielding of the specimen. This can be attributed to the load dependence of bond damage
and the bond-slip contribution.

4.4 Beam to Column Joint Subassemblage

Experimental studies have shown that bond slip has a pronounced effect on the global
deformation of beam-column subassemblages which are subjected to loading simulating the
effect of high lateral loads. This happens because the anchorage length of the reinforcing bars

4-8"

1"

3-11.5"

SHEAR
19.25" TRANSDUCER

2-7.5"

JOINT TIES(TOT. 7)
9" 17" 12@#6
REGULAR TIES(TOT. 10) 4@#6
AT EACH SIDE
16" 17"
COL TIES @ 1.6"(TOT. 13) 3@#5
AT EACH SIDE 0.75" #2
TIES
(a) (b) (c)

FIGURE 4.23 GEOMETRY AND REINFORCING DETAILS

OF BEAM TO COLUMN SUBASSEMBLAGE

in the joint does not suffice to transfer to the concrete the steel forces at the beam-column
interfaces of the joint. These forces cause the reinforcing bars to be pulled from one end of
the joint and pushed from the other, so that a force equal to almost twice the yield force of the
bar needs to be transferred to the concrete within the joint. When lateral load reversals occur,
the stress transfer problem is aggravated by bond deterioration under cyclic load reversals.
CHAPTER 4 105

V=470kips

FIGURE 4.24 FINITE ELEMENT IDEALIZATION OF BEAM TO COLUMN SUBASSEMBLAGE

In order to assess the ability of the proposed finite element model to simulate the
behavior of beam-column joint subassemblages, specimen BC4, which was tested by
Viwathanatepa et al. (1979a), is selected for further study. The subassemblage consisted of
two 9 in. by 16 in. girders framing into a 17 in. by 17 in. square column. The dimensions and
the reinforcing arrangement of the subassemblage are shown in Fig. 4.23. The column was
bolted at the top and bottom to steel clevises for mounting the specimen on the testing frame.
The main longitudinal reinforcement of the beams consisted of 4#6 bars at the top and 3#5
bars at the bottom of the section with #2 tied stirrups spaced at 3.5 in. as transverse
reinforcement (Fig. 4.23b). The longitudinal reinforcement of the column consisted of 12#6
bars with #2 ties spaced 1.6 in. center-to-center along column providing the transverse
reinforcement (Fig. 4.23c). Seven #2 ties were placed inside the beam-column joint region to
satisfy the confinement and shear resistance requirements of the building code. The material
properties of the specimen are summarized in Table 4.1. The bond stress-slip relation and any
other material properties not specifically mentioned in Table 4.1 are the same as those used in
106 CHAPTER 4

the analysis of beam J-4, since the concrete composition and strength of the two specimens is
not much different (Viwathanatepa et al. 1979a).
The beam-column subassemblage was subjected to a constant axial load P of 470 kips
which simulated the effect of gravity loads and a horizontal load H at the lower column end
which was cycled to simulate the effect of lateral loads on the subassemblage. Specimen BC4
was subjected to a very severe, pulse-type loading with a single load reversal in order to study
the monotonic behavior of the subassemblage and establish the load-displacement response
envelope. Other specimens were then cycled several times with gradually increasing
displacement ductility to simulate the effect of earthquake excitation on the subassemblage.
Since the load-displacement response of subassemblage BC4 during the first load cycle,
which represents a monotonic load test to near failure, was well established during the
experiment, it serves as an ideal case study for the proposed model.

36

30
TOTAL LOAD H (kips)

24

18

12

Experiment
6
Present Study

0
0 0.1 0.2 0.3 0.4 0.5 0.6

CENTER DEFLECTION (in)

FIGURE 4.25 LOAD-DISPLACEMENT RELATION OF

BEAM TO COLUMN SUBASSEMBLAGE

The finite element representation of the beam-column subassemblage is shown in Fig.

4.24. Concrete was modeled by 8-node isoparametric elements and the longitudinal and
transverse reinforcement was modeled by 2-node truss elements. The bond-slip effect is
included in the analysis with bond link elements, as described in Sections 2.4 and 2.5.
CHAPTER 4 107

36

30

TOTAL LOAD H (kips) 24

18

12

Experiment
6 Perfect Bond
Bond-Slip

0
0 0.1 0.2 0.3 0.4 0.5 0.6

CENTER DEFLECTION (in)

FIGURE 4.26 TENSION STIFFENING AND BOND-SLIP EFFECT

IN BEAM TO COLUMN SUBASSEMBLAGE

Fig. 4.25 compares the analytical results with the measured load-displacement
response of the subassemblage. With the effects of tension stiffening and bond-slip the
analysis shows excellent agreement with the experimental results. The lateral force of 36 kips
in Fig. 4.25 closely approximates the ultimate load of the specimen.
Fig. 4.26 shows that the effect of bond slip affects the behavior of the subassemblage
much more than tension stiffening. In this case the bond-slip of reinforcing bars in the joint
contributes approximately 33% of the total deformation of the subassemblage near the
ultimate load of 36 kips. The observation of the significance of bond-slip in the beam-column
subassemblage agrees with the results of the over-reinforced concrete beam A-1, where the
interaction of bond and shear plays a significant role in the response. This interaction also has
an important effect on the stress transfer in beam-column joints, particularly, when these are
subjected to loading simulating the effect of lateral loads. Even though the tension stiffening
effect plays a minor role in the response of this beam-column joint, it should not be
concluded that it can be excluded from the model, since it helps prevent numerical instability
problems in connection with crack formation and propagation.
108 CHAPTER 4

4.5 Reinforced Concrete Slabs

Tension stiffening has a significant effect in the analysis of RC slabs. In order to

investigate the validity of the proposed model a corner-supported two-way slab tested by
Jofriet and McNeice (1971) was used in the correlation studies. The slab is square in plan and
spans 36 in. with a thickness of 1.75 in. It is reinforced with wire mesh with a steel ratio of
0.85% in both directions and is subjected to a concentrated load at midspan, as shown in Fig.
4.27. The material properties of the specimen are summarized in Table 4.1. This slab was
previously studied by a number of investigators (Jofriet and McNeice 1971; Hand et al. 1973;
Lin and Scordelis 1975; Bashur and Darwin 1978) and it, therefore, provides a good case
study for the comparison of different models with measured data.

36in
2
3in

FIGURE 4.27 CONFIGURATION AND FINITE ELEMENT

DISCRETIZATION OF JOFRIET-MCNEICE SLAB

Fig. 4.27 shows the finite element representation of one quarter of the slab only,
because of symmetry considerations. Nine slab elements were used in the basic analysis with
eight concrete layers through the thickness of the slab for a total of 72 elements.
Fig. 4.28 compares the analytical load-deflection relation at node 2 in Fig. 4.27 with
measured data. The deflection at node 2 is obtained by interpolation of measured deflections
at nearby nodes. Fig. 4.28 shows that the predictions of the proposed tension-stiffening
model, which is based on a fracture energy criterion, are in very good agreement with
CHAPTER 4 109

experimental data. In the analysis of the same slab Joffriet and McNeice (1971) used the
effective stiffness approach and treated, both, steel and concrete as elastic materials; Lin and
Scordelis (1975) included the elasto-plastic behavior of reinforcing steel and concrete and
accounted for the coupling effect between membrane and bending action in the slab; Bashur
and Darwin (1978) used the analytical moment-curvature relation based on yield line theory
to analyze the slab; finally, Hand, Pecknold and Schnobrich (1973) considered both the shear
and the in-plane deformation of the slab, but neglected the tension-stiffening effect.

3.6

3
TOTAL LOAD P (kips)

2.4

1.8

1.2
Hand et al.
Jofriet and McNeice
Bashur and Darwin
0.6 Experiment
Present Study

0
0 4 8 12 16 20 24 28 32 36
-2
DEFLECTION AT NODE 2 ( x10 in )

FIGURE 4.28 LOAD-DEFLECTION RELATIONSHIP AT NODE 2 IN FIG. 4.27

In order to investigate the effect of finite element mesh size on the results, a nine and
a four element mesh is used in the finite element idealization of the slab and the results are
compared in Fig. 4.29. Solutions B and C which are obtained by including the tension
stiffening effect and the size of the element, as proposed in Section 2.2.3.3, are virtually
identical and agree very well with the experimental results. By contrast, solution E which
does not include the element size effect significantly overestimates the slab deflection. In this
case the crack band width extends over the entire element and the model is much more
flexible than the real structure. A much finer mesh is required to improve the accuracy of the
110 CHAPTER 4

results at the expense of computational efficiency. Fig. 4.29 also shows that neglecting the
strain softening of concrete in solution A results in a model which is stiffer than the actual
slab. Fig. 4.29 underlines the significance of tension stiffening in the analysis of slabs and
stresses the computational advantage of the proposed model, which is independent on finite
element size and thus computationally very efficient.
3.6

3
TOTAL LOAD P (kips)

2.4

1.8

1.2
SOLUTION A
SOLUTION B
SOLUTION C
0.6 SOLUTION D
SOLUTION E

0
0 4 8 12 16 20 24 28 32 36
-2
DEFLECTION AT NODE 2 [x10 in]

FIG. 4.29 EFFECT OF MESH SIZE ON LOAD-DEFLECTION RELATION AT NODE 2

A: NO STRAIN SOFTENING (9 ELEMENTS)
B: STRAIN SOFTENING AND ELEMENT SIZE EFFECT INCLUDED (9 ELEMENTS)
C: STRAIN SOFTENING AND ELEMENT SIZE EFFECT INCLUDED (4 ELEMENTS)
D: EXPERIMENT
E: CASE C WITHOUT ELEMENT SIZE EFFECT

Finally, the rectangular slab B7 which was tested by Cardenas and Sozen (1968) is
used as the last test case of this study. The slab was simply supported along two opposite
edges and free along the other two, as shown in Fig. 4.30. Uniformly distributed moments
were applied along the simply supported edges. The dimensions of the slab and the loading
arrangement are shown in Fig. 4.30. The slab was reinforced with wire mesh oriented at 45°
to the edges with equal amount of reinforcement in both directions.
CHAPTER 4 111

42in

θ =45

54in
SIMPLE SIMPLE
SUPPORT FREE EDGE
SUPPORT

FIGURE 4.30 CARDENAS RECTANGULAR SLAB B7

The slab is modeled with nine 8-node isoparametric elements using nine concrete
layers through the thickness of the slab for a total of 81 elements. The material properties
used in the analysis are summarized in Table 4.1. Fig. 4.31 compares the experimental
moment-curvature relation of the slab with the results of this study and an earlier study by
Vebo and Ghali (1977). The results of the present study which include the effects of tension
stiffening and finite element mesh size are in very satisfactory agreement with measured data.

6
UNIT MOMENT (in-kips/in)

Experiment
Vebo and Ghali
Present Study

0
0 4 8 12 16 20
-4 -1
CURVATURE ( x10 in )
FIGURE 4.31 MOMENT-CURVATURE RELATIONSHIPS OF B7
 CHAPTER 5
CONCLUSIONS

The objective of this study is to develop reliable and computationally efficient finite
element models for the analysis of reinforced concrete beams, slabs and beam-column joint
subassemblages under monotonic loading conditions. Since this is part of a larger study
which endeavors to study the hysteretic behavior of reinforced concrete members, the models
developed in this first part should also be general enough, so that they can be extended to
cyclic loading conditions.

It is assumed that the behavior of RC beams, slabs and beam-column joints can be
described by a plane stress field. Concrete and reinforcing steel are represented by separate
material models which are combined together with a model of the interaction between
reinforcing steel and concrete through bond-slip to describe the behavior of the composite
reinforced concrete material. The material behavior of concrete is described by two failure
surfaces in the biaxial stress space and one failure surface in the biaxial strain space.
Concrete is assumed as a linear elastic material for stress states which lie inside the initial
yield surface. For stresses outside this surface the behavior of concrete is described by a
nonlinear orthotropic model, whose axes of orthotropy are parallel to the principal strain
directions. The concrete stress-strain relation is derived from equivalent uniaxial relations in
the axes of orthotropy. The behavior of cracked concrete is described by a system of
orthogonal cracks, which follow the principal strain directions and are thus rotating during
the load history. Crushing or cracking of concrete takes place when the strains lie outside the
ultimate surface in the biaxial strain space.

A new smeared finite element model is proposed based on an improved cracking

criterion, which is derived from fracture mechanics principles. This model retains objectivity
of the results for very large finite elements, since it considers cracking to be concentrated
over a small region around the integration point and not over the entire finite element, as do
previous models.

A new reinforcing steel model which is embedded inside a concrete element, but
accounts for the effect of bond-slip is developed. This model results in significant savings in

113
114 CHAPTER 5

the number of nodes needed to account for the effect of bond-slip, particularly, in three
dimensional finite element models. A new nonlinear solution scheme is developed in
connection with this model.

The correlation studies between analytical and experimental results and the
parametric studies associated with them lead to the following conclusions:

• The effect of tension-stiffening is important in the analysis of RC beams, but more so

in the analysis of RC slabs.

• The inclusion of tension-stiffening is important for the independence of the analytical

results from the size of the finite element mesh, but also for avoiding numerical
problems in connection with crack formation and propagation.

• Present smeared models are too stiff in connection with large finite elements. A new
criterion which limits the effect of tension stiffening to the vicinity of the integration
point yields very satisfactory results.

• The effect of bond-slip is very important in the analysis of RC beams and beam-
column subassemblages even under monotonic loads. This effect is more pronounced
in heavily reinforced beams and in beam-column joints, which are subjected to loads
simulating the effect of lateral loads on the structure.

• Tension-stiffening and bond-slip cause opposite effects on the response of RC

members. While tension-stiffening, which accounts for the concrete tensile stresses
between cracks, increases the stiffness of the member, bond-slip leads to a stiffness
reduction. In lightly reinforced beams, these effects can cancel each other at certain
load stages, thus, leading to the erroneous impression that they can be neglected in the
analysis. Since bond-slip increases with loading, while tension-stiffening does not,
consistent results can only be obtained when both effects are included in the model.
Moreover, the effect of bond-slip clearly outweighs the contribution of tension-
stiffening in heavily reinforced beams and beam-column joint subassemblages. In
these cases the exclusion of the bond-slip effect can lead to significant overestimation
of the stiffness of the member.

• The tensile strength of concrete has insignificant effect on the load-displacement

response of RC beams. Crack formation and propagation is rather influenced by the
fracture energy of concrete.
CHAPTER 5 115

• The value of the cracked shear constant λ has negligible effect on the response of
slender RC beams in bending, since most of the member is in the stress state of
compression-compression or tension-tension and only a small portion of the structure
near the supports is in a stress state of tension-compression. It is only under a stress
state of tension-compression that the coupling between normal and shear strains
represented by the cracked shear constant λ plays an important role. While this
coupling effect is negligible in the case of slender beams in bending, it becomes
significant in the analysis of shear walls and deep beams.
 REFERENCES

Abdel Rahman, H.H. and Hinton, E. (1986). "Nonlinear Finite Element Analysis of Reinforced Concrete
Stiffened and Cellular Slabs". Computers & Structures, Vol. 23, No. 3, pp. 333-350.
Adeghe, L.N. and Collins, M.P. (1986). "A Finite Element Model for Studying Reinforced Concrete Detailing
Problems". Publication No. 86-12, Department of Civil Engineering, University of Toronto.
Arnesen, A., Sorensen, S. I. and Bergan, P.G. (1980). "Nonlinear Analysis of Reinforced Concrete". Computers
& Structures, Vol. 12, pp. 571-579.
ASCE Task Committee on Finite Element Analysis of Reinforced Concrete Structures. (1982). State-of-the-Art
Report on Finite Element Analysis of Reinforced Concrete, ASCE Special Publications.
Balakrishnan, S. and Murray, D.W. (1988). "Concrete Constitutive Model for NLFE Analysis of Structures".
Journal of Structural Engineering, ASCE, Vol. 114, No. 7, pp. 1449-1466.
Barzegar, F. and Schnobrich, W.C. (1986). "Nonlinear Finite Element Analysis of Reinforced Concrete under
Short Term Monotonic Loading". Civil Engineering Studies SRS No. 530, Univ. of Illinois at Urbana,
Illinois.
Bashur, F.K, and Darwin, D. (1978). "Nonlinear Model for Reinforced Concrete Slabs". Journal of Structural
Division, ASCE, Vol. 104, No. ST1, pp. 157-170.
Bazant, Z.P. and Cedolin, L. (1980). "Fracture Mechanics of Reinforced Concrete". Journal of the Engineering
Mechanics, ASCE, Vol. 106, No. EM6, pp. 1287-1306.
Bazant, Z.P. and Cedolin, L.. (1983). "Finite Element Modeling of Crack Band Propagation". Journal of
Structural Engineering, ASCE, Vol. 109, No. 1, pp. 69-93.
Bazant, Z.P. and Oh, B.H. (1983). "Crack Band Theory for Fracture of Concrete". Materials and Structures,
RILEM, Paris, Vol. 16, pp. 155-176.
Bazant, Z.P. and Ozbolt, J. (1989) "Nonlocal Microplane Model for Fracture, Damage, and Size Effect in
Structures", Journal of Engineering Mechanics, ASCE, Vol. 116, No. 11.
Becker, J.M. and Bresler, B.. (1974). "FIRES-RC-A Computer Program for the Fire Response of Structures-
Reinforced Concrete Frames". Report No. UCB/FRG 74-3, Department of Civil Engineering, University
of California, Berkeley.
Bedard, C. and Kotsovos, M.D.. (1986). "Fracture Processes of Concrete for NLFEA Methods". Journal of
Structural Engineering, ASCE, Vol. 112, No. 3, pp. 573-587.
Bergmann, R. and Pantazopoulou, V. A. (1988). "Finite Element for R/C Shear Walls Under Cyclic Loads".
Department of Civil Engineering, Report UCB/SEMM-88/09, University of California, Berkeley, 1988.
de Borst, R. and Nauta, P. (1985) "Non-Orthogonal Cracks in a Smeared Finite Element Model". Engineering
Computations, Vol. 2, pp. 35-46.
Bresler, B. and Scordelis, A.C. (1963). "Shear Strength of Reinforced Concrete Beams". Journal of ACI,
Vol. 60, No. 1, pp. 51-72.
Bresler, B. and Bertero, V.V. (1968). "Behavior of Reinforced Concrete Under Repeated Load". Journal of the
Structural Division, ASCE, Vol. 94, No. ST6, pp. 1567-1590.
Broek, D. (1974). Elementary Engineering Fracture Mechanics, Sijthoff & Noordoff, London.
Burns, N.H. and Siess, C.P. (1962). "Load-Deformation Characteristics of Beam-Column Connections in
Reinforced Concrete". Civil Engineering Studies, SRS No. 234, University of Illinois, Urbana.

117
118 REFERENCES

Cardenas A. and Sozen, M.A. (1968). "Strength and Behavior of Isotropally and Nonisotropally Reinforced
Concrete Slabs Subjected to Combinations of Flexural and Torsional Moments". Civil Engineering
Studies SRS No. 336, Univ. of Illinois at Urbana, Illinois.
Comité Euro-International du Beton (CEB) (1978). CEB-FIP Model Code of Concrete Structures. Bulletin
d'Information 124/125E, Paris, France.
Cervenka, V. (1970). "Inelastic Finite Element Analysis of Reinforced Concrete Panels". Ph.D. Dissertation,
University of Colorado, Boulder.
Cervenka, V., Eligehausen, R. and Pukl, R. (1990). "SBETA-Computer Program for Nonlinear Finite Element
Analysis of Reinforced Concrete Structures". Report 90/1, Institute of Building Materials, University of
Stuttgart.
Chang, T.Y., Taniguchi, H. and Chen, W. F. (1987). "Nonlinear Finite Element Analysis of Reinforced
Concrete Panels". Journal of Structural Engineering, ASCE, Vol. 113, No. 1, pp. 122-140.
Chen, W.F. (1976). Plasticity in Reinforced Concrete, McGraw-Hill, New York.
Chia, C.Y. (1978). Nonlinear Analysis of Plates, McGraw-Hill, New York.
Ciampi, V., Eligehausen, R., Bertero, V.V. and Popov, E.P. (1982). "Analytical Model for Concrete Anchorages
of Reinforcing Bars under Generalized Excitations". Report No EERC 82-23, Earthquake Engineering
Research Center, University of California, Berkeley.
Cope, R.J., Rao, P.V., Clark, L.A. and Norris, R. (1980). "Modeling of Reinforced Concrete Behavior for Finite
Element Analysis of Bridge Slabs". Numerical Methods for Nonlinear Problems , C. Taylor, E. Hinton
and D.R.J. Oden, eds., Pineridge Press, Swansea, pp. 457-470.
Crisfield, M.A. (1982). "Accelerated Solution Techniques and Concrete Cracking". Computer Methods in
Applied Mechanics and Engineering, Vol. 33, pp. 585-607.
Darwin, D. and Pecknold, D.A. (1977). "Analysis of Cyclic Loading of Plane R/C Structures". Computers &
Structures, Vol. 7, No. 1, pp. 137-147.
Desai, C.S. and Siriwardance, H.J. (1972). Constitutive Laws for Engineering Materials, Prentice-Hall,
Englewood Cliffs, New Jersey.
Dotroppe, J.C., Schnobrich, W.C. and Pecknold, D.A. (1973). "Layered Finite Element Procedure for Inelastic
Analysis of Reinforced Concrete Slabs". IABSE Publication, Vol. 33-11.
Dougill, J.W. (1976). "On Stable Progressively Fracturing Solids," Journal of Applied Mathematics and
Physics, Vol. 27, pp. 423-437.
Eligehausen, R., Popov, E.P. and Bertero, V.V. (1983). "Local Bond Stress-Slip Relationships of Deformed
Bars Under Generalized Excitations". Report No. UCB/EERC 83-23, Earthquake Engineering Research
Center, University of California, Berkeley.
Filippou, F.C., Popov, E.P., and Bertero, V.V. (1983). "Effects of Bond Deterioration on Hysteretic Behavior of
Reinforced Concrete Joints". Report No. UCB/EERC-83/19, Earthquake Engineering Research Center,
University of California, Berkeley.
Filippou, F.C. (1986). "A Simple Model for Reinforcing Bar Anchorages Under Cyclic Excitations". ASCE,
Journal of Structural Engineering, Vol. 112, No. 7, pp. 1639-1659.
Franklin, H.A. (1970). "Non-Linear Analysis of Reinforced Concrete Frames and Panels". Ph.D. Dissertation,
Division of Structural Engineering and Structural Mechanics, University of California, Berkeley,
SEMM 70-5.
Gaston, J.R., Siess, C.P. and Newmark, N.M. (1972). "A Layered Finite Element Non-Linear Analysis of
Reinforced Concrete Plates and Shells". Civil Engineering Studies, SRS No. 389, University of Illinois,
Urbana.
Gilbert, R.I. and Warner, R.F. (1978). "Tension Stiffening in Reinforced Concrete Slabs". Journal of Structural
Division, ASCE, Vol. 104, No. ST12, pp. 1885-1900.
REFERENCES 119

de Groot, A.K., Kusters, G.M.A. and Monnier, T. (1981). "Numerical Modeling of Bond-Slip Behavior".
Heron, Concrete Mechanics, Vol. 26, No. 1B.
Gupta, A.K. and Akbar, H. (1983). "A Finite Element for the Analysis of Reinforced Concrete Structures".
International Journal for Numerical Methods in Engineering, Vol. 19, pp. 1705-1712.
Hand, F.R., Pecknold, D.A. and Schnobrich, W.C. (1973). "Nonlinear Layered Analysis of RC Plates and
Shells". Journal of Structural Division, ASCE, Vol. 99, No. ST7, pp. 1491-1505.
Hayashi, S. and Kokusho, S. (1985). "Bond Behavior in the Neighborhood of the Crack". Proceedings of the
U.S.-Japan Joint Seminar on Finite Element Analysis of Reinforced Concrete, pp. 364-373, Tokyo,
Japan.
Hillerborg, A., Modeer, M. and Petersson, P.E. (1976). "Analysis of Crack Formation and Growth in Concrete
by Means of Fracture Mechanics and Finite Element". Cement and Concrete Research, Vol. 6, No. 6,
pp. 773-782.
Hognestad, E. (1951). "A Study of Combined Bending and Axial Load in Reinforced Concrete Members".
University of Illinois Engineering Experiment Station, Bulletin Series No. 399, Bulletin No. 1.
Jain, S.C. and Kennedy, J.B. (1974). "Yield Criterion for Reinforced Concrete Slabs," Journal of Structural
Division, ASCE, Vol. 100, No. ST3, pp. 631-644.
Jenq, Y.S. and Shah, S.P. (1986). "Crack Propagation in Fiber-Reinforced Concrete". Journal of Structural
Engineering, ASCE, Vol. 112, No. 1, pp. 19-34.
Jofriet, J.C. and McNeice, G.M. (1971). "Finite Element Analysis of RC Slabs". Journal of Structural Division,
ASCE, Vol. 97, No. ST3, pp. 785-806.
Keuser, M. and Mehlhorn, G. (1987). "Finite Element Models for Bond Problems". Journal of Structural
Engineering, ASCE, Vol. 113, No. 10, pp. 2160-2173.
Kupfer, H., Hilsdorf, H.K. and Rüsch, H. (1969). "Behavior of Concrete under Biaxial Stresses". ACI Journal,
Vol. 66, No. 66-62, pp. 656-666.
Kwak, H.G. (1990). "Material Nonlinear Finite Element Analysis and Optimal Design of Reinforced Concrete
Structures". Ph.D. Dissertation, Department of Civil Engineering, KAIST, Korea.
Leibengood, L.D., Darwin, D. and Dodds, R.H. (1986). "Parameters Affecting FE Analysis of Concrete
Structures," Journal of Structural Engineering, ASCE, Vol. 112, No. 2, pp. 326-341.
Lin, C.S. and Scordelis, A.C. (1975). "Nonlinear Analysis of RC Shells of General Form". Journal of Structural
Division, ASCE, Vol. 101, No. ST3, pp. 523-538.
Meyer, C. and Okamura, H., eds. (1985). "Finite Element Analysis of Reinforced Concrete Structures".
Proceedings of the US-Japan Joint Seminar on Finite Element Analysis of Reinforced Concrete, Tokyo,
Japan.
Milford, R.V. and Schnobrich, W.C. (1984). "Nonlinear Behavior of Reinforced Concrete Cooling Towers".
Civil Engineering Studies SRS No. 514, Univ. of Illinois at Urbana, Illinois.
Nayak, G.C. and Zienkiewicz, O.C. (1972). "Elasto-Plastic Stress Analysis". International Journal of Numerical
Methods in Engineering, Vol. 5, pp. 113-135.
Ngo, D. and Scordelis, A.C. (1967). "Finite Element Analysis of Reinforced Concrete Beams," Journal of ACI,
Vol. 64, No. 3, pp. 152-163.
Nilson, A.H. (1972). "Internal Measurement of Bond Slip". Journal of ACI, Vol. 69, Title No. 7, pp. 439-441.
Owen, J. and Hinton, E. (1978). Finite Elements in Plasticity, Pineridge Press, Swansea.
Park, P. and Paulay, T. (1975). Reinforced Concrete Structures John Wiley & Sons, New York.
Popovics, S. (1969). "A Review of Stress-Strain Relationships for Concrete". Journal of ACI, Vol. 66, No. 5,
pp. 756-764.
Rajagopal, K.R. (1976). "Nonlinear Analysis of Reinforced Concrete Beams, Beam-Columns and Slabs by
Finite Elements". Ph.D. Dissertation, Iowa State University.
 120 REFERENCES

Rashid, Y.R. (1968). "Analysis of Prestressed Concrete Pressure Vessels". Nuclear Engineering and Design,
Vol. 7, No. 4, pp. 334-344.
Reinhardt, H.W. (1986). "The Role of Fracture Mechanics in Rational Rules for Concrete Design". IABSE
Surveys, Vol. S-34, pp. 1-15.
Rots, J.G., Nauta, P., Kusters, G.M.A. and Blaauwendraad, J. (1985). "Smeared Crack Approach and Fracture
Localization in Concrete". HERON, Vol. 30, No. 1, pp. 3-48.
Scanlon, A. and Murray, D.W. (1974). "Time Dependent Reinforced Concrete Slab Deflections". Journal of the
Structural Division, ASCE, Vol. 100, No. ST9, pp. 1911-1924.
Scordelis, A.C., Ngo, D. and Franklin, H.A. (1974). "Finite Element Study of Reinforced Concrete Beams with
Diagonal Tension Cracks". Proceedings of Symposium on Shear in Reinforced Concrete, ACI
Publication SP-42.
Selna, L.G. (1969). "Creep, Cracking and Shrinkage in Concrete Frame Structures". Journal of the Structural
Division, ASCE, Vol. 95, No. ST12, pp. 2743-2761.
Suidan, M.T. and Schnobrich, W.C. (1973). "Finite Element Analysis of Reinforced Concrete". Journal of the
Structural Division, ASCE, Vol. 99, No. ST10, pp. 2109-2122.
Tasuji, M.E., Slate, F.O. and Nilson, A.H. (1978). "Stress-Strain Response and Fracture of Concrete in Biaxial
Loading". Journal of ACI, Vol. 75, No. 3, pp. 306-312.
Taylor, C. and Hinton, E. (1980). Numerical Methods for Nonlinear Problems, Pineridge Press, Swansea.
Vebo, A. and Ghali, A. (1977). "Moment-Curvature Relation of Reinforced Concrete Slabs". Journal of
Structural Division, ASCE, Vol. 103, No. ST3, pp. 515-531.
Vecchio, F. and Collins, M.P. (1982). "The Response of Reinforced Concrete to In-Plane Shear and Normal
Stress". Publication No. 82-03, Department of Civil Engineering, University of Toronto, Toronto,
Canada.
Viwathanatepa, S., Popov, E.P. and Bertero, V.V. (1979a). "Seismic Behavior of Reinforced Concrete Interior
Beam-Column Subassemblages". Report No. UCB/EERC-79/14, Earthquake Engineering Research
Center, University of California, Berkeley.
Viwathanatepa, S., Popov, E.P. and Bertero, V.V. (1979b). "Effects of Generalized Loadings on Bond of
Reinforcing Bars Embedded in Confined Concrete Blocks". Report No EERC 79-22, Earthquake
Engineering Research Center, University of California, Berkeley.
Wamba, T.E. (1980). "Nonlinear Analysis of Reinforced Concrete Beams by Finite Element Method". Ph.D.
Dissertation, Department of Civil Engineering, State University of New York at Buffalo.
Welch, G.B. and Haisman, B. (1969). "Fracture Toughness Measurements of Concrete". Report No. R42,
University of New South Wales, Sydney, Australia.
Wunderlich, W., Stein, E. and Bathe, K.J., eds. (1980). Nonlinear Finite Element Analysis in Structural
Mechanics, Proceedings of the Europe-US Workshop, Ruhr-University, Germany.
Yankelevsky, D.Z. (1985). "New Finite Element for Bond-Slip Analysis", Journal of Structural Engineering,
ASCE, Vol. 111, No. 7, pp. 1533-1542.
Yoshikawa, Y. and Tanabe, T. (1987). "An Analytical Model for Frictional Shear Slip of Cracked Concrete".
IABSE Colloquium Computational Mechanics of Concrete Structures, pp. 75-86, Delft University,
Netherlands.
Zulfiqar, N. and Filippou, F.C. (1990). "Models of Critical Regions in Reinforced Concrete Frames under
Seismic Excitations". Report No. UCB/EERC-90/06, Earthquake Engineering Research Center,
University of California, Berkeley.

View publication stats