Elasto-plastic And Damage Modeling Of Reinforced Concrete

    Reinforced concrete sample abaqus. 3.3.3 examples of using nonlinear equilibrium algorithms to solvemodeling the mechanical behavior of reinforced concrete (rc) is still one of the software abaqus via user defined material subroutine umat to analyze and better. etd.lsu.edu.

  • Id: # ca835ac
  • File Type: pdf
  • Author: etd.lsu.edu
  • File Size 1.49 MB
  • Read by User: 90 Times
  • Published: Thursday, April 3, 2014
  • index: Reinforced Concrete Sample Abaqus

Rating

  • Read Online
ELASTO-PLASTIC AND DAMAGE MODELING OF REINFORCED
CONCRETE
A Dissertation
Submitted to the Graduate Faculty of the
Louisiana State University and
Agricultural and Mechanical College
in partial fulfillment of the
requirements for the degree of
Doctor of Philosophy
in
The Department of Civil & Environmental Engineering
by
Ziad N. Taqieddin
B.S., Applied Science University, Amman, Jordan, 2001
M.S., Louisiana State University, Baton Rouge, Louisiana, 2005
August 2008
ii
TO THOSE WHO TAUGHT ME EVERYTHING THAT
MATTERS
TO MY BELOVED MOTHER & FATHER
iii
ACKNOWLEDGEMENTS
I would like to gratefully and sincerely thank Dr. George Z. Voyiadjis for his
impeccable guidance, knowledge, understanding, patience, and most importantly, his
friendship during my graduate studies at LSU. His mentorship was paramount in
providing me with a well rounded experience throughout the years. Prof. Voyiadjis
encouraged me not only to grow as a scientist but also as an instructor and an
independent thinker. His enthusiasm and unlimited zeal have been major driving forces
throughout my graduate career. For everything you’ve done to help me, Dr. Voyiadjis, I
can’t thank you enough.
I would also like to extend my thanks to my exceptional doctoral committee members,
Dr. Steve Cai, Dr. Su-Seng Pang, Dr. Eurico D’Sa and especially Dr. Suresh Moorthy for
teaching me what Finite Element Analysis is all about.
I would also like to thank all of the members of the advanced Computational Solid
Mechanics (CSM) Laboratory. Through five years of Masters and PhD, I developed a
very strong bond of friendship with a dozen of the CSM Tigers (as we like to call
ourselves). These friends were from all over the world yet they were my family in Baton
Rouge. I attended the graduation of many of them, just as many of them will attend the
graduation of mine. My experience with them is one of the memories that I will hold dear
to my heart forever. To all of you CSM Tigers, thank you. For you, my dear friends, have
provided much needed humor and entertainment in what could have otherwise been a
stressful and boring concentrated research environment.
Finally, and most importantly, I would like to express my deepest gratitude to my
family members:
• My parents: Dr. Noureddin and Mona. Thank you for being the infinite source of
love and support that enlightened my path from six thousand miles away. Thank
you for always believing in me and for encouraging me to be strong and focused
in the times I needed encouragement most. Thank you for welcoming my dreams
and making sure they become true throughout the years.
• My brother and sisters. Thank you for all the joy that you bring into my life.
Thank you for sending me pictures of your new born babies and for allowing me
to be a part of their daily life as an uncle even though I’m distantly far away.
Thank you for being the wonderful people that you are.
• My wife and best friend, Rawan. Your support, encouragement, quiet patience
and unwavering love were undeniably the bedrock upon which the past two years
of my life have been built. Your sprit lifting and cheering smile conquered all my
fears. Your tolerance of my occasional stress-outbursts is a testament in itself of
your unyielding love and devotion. You are the reason why my glass is always
half full, and your happiness and welfare will be the primary goal in my life.
iv
TABLE OF CONTENTS
DEDICATION . . . . . . . . . ii
ACKNOWLEDGMENTS . . . . . . . . iii
ABSTRACT . . . . . . . . . . vii
CHAPTER 1: INTRODUCTION . . . . . . . 1
1.1 Scope and Objectives . . . . . . . . 3
1.2 Outline of the Dissertation . . . . . . . 4
CHAPTER 2: LITERATURE REVIEW . . . . . . 6
2.1 Concrete Materials . . . . . . . . 8
2.1.1 Concrete Plasticity . . . . . . . . 13
2.1.2 Concrete Damage Mechanics . . . . . . 14
2.2 Steel Reinforcement . . . . . . . . 19
2.3 Bond and Interaction at Steel – Concrete Interface . . . . 21
2.4 Methodology . . . . . . . . . 24
2.4.1 Concrete Plasticity . . . . . . . . 24
2.4.2 Concrete Damage Mechanics . . . . . . . 27
2.4.3 Steel Reinforcement . . . . . . . . 30
2.4.4 Steel-Concrete Bond and Interaction . . . . . . 31
2.4.5 Numerical Implementation of the Model . . . . . 31
CHAPTER 3: REINFORCING STEEL MATERIAL MODEL . . . 33
3.1 Introduction . . . . . . . . . 33
3.2 The Application of the FE Method to Nonlinear Continuum Mechanics . 37
3.2.1 FE Analysis Procedure for Elastic Problems . . . . . 37
3.2.2 FE Analysis Procedure for Elastic-Plastic Problems . . . . 40
3.3 Algorithms for Nonlinear Global/Equilibrium Iterations . . . . 42
3.3.1 The Newton-Raphson Method . . . . . . 42
3.3.2 The Modified Newton-Raphson Method . . . . . 44
3.3.3 Examples of Using Nonlinear Equilibrium Algorithms to Solve
FE Elastic-Plastic Problems . . . . . . . 45
3.3.3.1 Solution Using Newton-Raphson Scheme . . . . 46
3.3.3.2 Solution Using the Modified Newton-Raphson Scheme . . . 48
3.4 Algorithms for Nonlinear Local Iterations . . . . . 52
3.4.1 The Incremental Constitutive Theory of Metal Plasticity . . . 52
3.4.2 Integrating the Incremental Constitutive Equations . . . . 59
3.4.2.1 Explicit (Forward-Euler) Integration Algorithm . . . . 60
3.4.2.2 Implicit Integration Algorithm: Radial Return Method . . . 63
3.4.2.3 Other Forms of Strain Hardening: Nonlinear Strain Hardening
and Linear Strain Hardening Following Perfectly-Plastic Behavior . 67
3.5 Implementation and Verification of the Integration Schemes . . . 69
v
CHAPTER 4: CONCRETE MATERIAL MODEL . . . . 76
4.1 Introduction . . . . . . . . . 76
4.2 Framework for the Elastic-Plastic-Damage Model . . . . 78
4.3 Consistent Thermodynamic Formulation . . . . . . 80
4.4 The Helmholtz Free Energy Function . . . . . . 84
4.4.1 The Elastic/Damage Part of the Helmholtz Free Energy Function . . 84
4.4.2 The Plastic Part of the Helmholtz Free Energy Function . . . 89
4.5 Plasticity Formulation . . . . . . . . 91
4.5.1 Plasticity Yield Surface and Hardening Functions . . . . 91
4.5.2 Plasticity Non-associative Flow Rule and Consistency Condition . . 93
4.6 Damage Formulation . . . . . . . . 94
4.6.1 Tensile and Compressive Damage Surfaces and Hardening Functions . 94
4.6.2 Damage Consistency Conditions . . . . . . 96
4.7 Numerical Integration of the Constitutive Model . . . . . 96
4.7.1 The Effective (Undamaged) Elastic-Plastic Steps . . . . 97
4.7.1.1 The Elastic – Predictor . . . . . . . 97
4.7.1.2 The Plastic – Corrector . . . . . . . 98
4.7.1.3 The Effective Configuration Integration Algorithm . . . 103
4.7.2 The Degradation (Damage) Step . . . . . . 107
4.7.2.1 The Damage - Corrector Step . . . . . . 107
4.8 The Consistent Elastic – Plastic – Damage Tangent Operator . . . 109
4.9 Implementation and Verification of the Integration Scheme . . . 114
4.9.1 Identification of the Proposed Model’s Parameters . . . . 114
4.9.2 Monotonic Uniaxial Tensile Test . . . . . . 115
4.9.3 Monotonic Uniaxial Compressive Test . . . . . 117
4.9.4 Monotonic Biaxial Tests . . . . . . . 119
4.9.5 Three Point Bending Test of a Notched Beam . . . . 126
CHAPTER 5: FINITE ELEMENT ANALYSIS OF REINFORCED
CONCRETE BEAMS . . . . . . 129
5.1 Introduction . . . . . . . . . 129
5.2 Bond-Slip Effect on the Stress – Strain Relation of the Steel
Reinforcement . . . . . . . . . 130
5.3 Applications of the Constitutive Model to the FE Analysis of RC Beams . 132
CHAPTER 6: SUMMARY, CONCLUSIONS AND FUTURE WORK . . 143
6.1 Summary and Conclusions . . . . . . . 143
6.2 Future Work . . . . . . . . . 146
BIBLIOGRAPHY . . . . . . . . . 148
APPENDIX A: NUMERICAL PROCEDURE FOR THE SPECTRAL
DECOMPOSITION CONCEPT USING THE
MATHEMATICAL SOFTWARE MAPLE
(CHAPTER 4) . . . . . . . 163
vi
APPENDIX B: UPDATING THE EFFECTIVE STRESS USING THE
RETURN MAPPING EQUATION (CHAPTER 4) . . 165
VITA . . . . . . . . . . . 166
vii
ABSTRACT
Modeling the mechanical behavior of Reinforced Concrete (RC) is still one of the
most difficult challenges in the field of structural engineering. The Nonlinear Finite
Element Analysis (NFEA) and modeling of the behavior of RC members are the primary
goals of this study. The macroscopic components of RC, Concrete material and
reinforcing steel, are represented herein by separate material models. These material
models are combined together using a model that describes the global effect of
interaction between reinforcing steel and concrete in order to simulate the behavior of the
composite RC material.
A thermodynamically consistent constitutive model for concrete that incorporates
concrete-plasticity and fracture-energy-based continuum damage mechanics is presented.
An effective stress space plasticity yield criterion, with multiple hardening functions and
a non-associative plasticity flow rule, is used simultaneously with two (tensile and
compressive) isotropic damage criteria. The spectral decomposition of the stress tensor
into tensile and compressive components is utilized in all criteria in order to simulate
different responses of the material under various loading patterns. The damage criteria
are based on the hydrostatic-deviatoric sensitive damage energy release rates in tension
and compression derived from the Helmholtz free energy function. Three dissipation
mechanisms are defined, one for plasticity and two for damage, to control the dissipation
process of the material model.
Elastic-plastic models that account for isotropic perfectly-plastic and plastic-strain-
hardening (linear, bilinear and nonlinear) of the steel reinforcement are provided as well.
The global effect of bond-slip is incorporated into the stress-strain diagram of the
reinforcing bars in an attempt to describe this interaction phenomenon in a stress-strain
driven environment.
The Numerical implementation and application are important parts of this study. A
suitable elastoplasticity-implicit/damage-explicit scheme is adapted here for the
integration of the incremental constitutive equations. The elastic-predictor, plastic-
corrector and damage-corrector steps are used to facilitate the integration procedure. The
constitutive approach is implemented, through numerical algorithms; in the advanced FE
software ABAQUS via user defined material subroutine UMAT to analyze and better
describe the overall behavior of such a composite material. Concrete and RC beams
subjected to static-short-term-monotonic loading are analyzed in an assumed isothermal
environment. The simulated results are compared to experimental studies conducted by
other researchers.
1
CHAPTER 1
INTRODUCTION
Reinforced Concrete (RC) is one of the most commonly used building materials
nowadays. It is a composite material made of plain concrete, which possesses relatively
high compressive strength but low tensile strength, and steel bars embedded in the
concrete, which can provide the needed strength in tension. The economy, efficiency,
strength and stiffness of RC make it an attractive material for a wide range of structural
engineering applications, such as nuclear power-plants, bridges, cooling towers and
offshore platforms. For RC to be used as a structural material, it should satisfy special
criteria including:
• Strength and Stiffness
• Safety and Appearance
• Economy
By applying the principles of structural analysis, the laws of equilibrium and the
consideration of the mechanical properties of the components studied; RC design
procedure should yield a sufficient margin of safety against collapse under ultimate loads.
Serviceability analysis is conducted to control the deflections under service loads and to
limit the crack width to an acceptable level for the structural component to perform and
appear safe and inhabitable for the human eye. Economical considerations are satisfied
by optimizing the usage of steel/concrete quantities to account for the difference in unit
costs of steel and concrete.
The ultimate objective of design is the safety and economy of the RC structural
member. The design process is usually based on a linear elastic analysis to calculate the
internal forces in the member which are then used to design the reinforcement and the
details of the member using some code provision. Codes are usually based on empirical
approaches that utilize experimental data and provide design rules to satisfy safety and
serviceability requirements. Although the design of RC structures based on linear-elastic
stress analysis is adequate and reliable in many cases, the extent and impact of a disaster
in terms of human and economical losses in the event of structural failure of large scale
modern structures necessitate more careful and detailed structural safety analysis. Thus,
Nonlinear Finite Element Analysis (NFEA) is often required to obtain detailed
information regarding the ultimate loading capacity and the post-failure behavior of RC
structures. The importance and interaction of different nonlinear effects on the response
of RC structures can be studied analytically using NFEA.
The complex behavior of concrete, which arises from the composite nature of the
material, is characterized by a reduction of the load carrying capacity with increasing
deformations after reaching a certain limit load. This global behavior is usually caused by
a material behavior which is described as strain softening and occurs in tension and in
compression. This necessitates the development of appropriate constitutive models to
describe such behavior.
2
In RC, the response of the structure is even more complicated. In general a number of
cracks will develop in the structure due to the bond action between concrete and
reinforcement. This results in a redistribution of the tensile loads from concrete to the
reinforcement. This phenomenon is called tension-stiffening, because the response is
stiffer than the response with a brittle fracture approach.
The behavior of RC is highly nonlinear which is caused by mechanisms such as
cracking, crushing, creep and shrinkage of concrete, but also caused by the interaction
between reinforcement and concrete, where the load transferring mechanism of the
interface between concrete and reinforcement plays an important role. Because all these
mechanisms are interacting, it is not realistic to try to formulate a constitutive model
which incorporates all these mechanisms, but a model has to be formulated to adequately
describe the behavior of a structure within the range of application which has been
restricted in advance. Although the constitutive models which are developed within this
phenomenological approach are usually simplified representations of the real behavior of
the material, it is believed that more insight can be gained by tracing the entire response
of a structure in this manner, than modeling a structure with highly sophisticated material
models which do not result in a converged solution after failure load and are
computationally expensive and complicated.
A large variety of models have been proposed to characterize the stress-strain relation
and failure behavior of RC materials. All these models have certain inherent advantages
and disadvantages which depend to a large degree on their particular application and
complexity. Macroscopic constitutive studies have been conducted with different levels
of complexity and applicability in order to address the different aspects of the concrete
material behavior. On the other hand, microscopic modeling and multi-scale modeling
offer useful ways to model the material behavior, but their applicability to full-scale
structural problems is still problematic, due to their requirement for huge amounts of
computer resources. Therefore, further development in the macroscopic constitutive
modeling of concrete is justified and needed, with the motivation of incorporating
contemporary experimentally observed features of the material behavior in the modeling.
In this study, constitutive models representing the macro constituents of RC (Concrete
and Steel Reinforcement) are developed with emphasis on the rigor and consistency in
formulation and implementation into the Finite Element Analysis (FEA) commercial
software ABAQUS.
Concrete and reinforcing steel are represented herein by separate material models
which are combined together using a model that describes the interaction between
reinforcing steel and concrete to simulate the overall behavior of the composite RC
material. An elasto-plastic-damage constitutive model is used to describe the behavior of
concrete, while steel reinforcement is modeled as an elastoplastic material with strain
hardening using the classical von Mises plasticity. Bond considerations are accounted for
within the steel reinforcement model. Coupling between damage and plasticity in the
constitutive model is employed to capture the observed phenomenological behavior of
concrete. In this combined approach, damage theory is used to model the material
deterioration, while the permanent deformation and some other behavioral features of
3
concrete can be captured using plasticity theory. All features of the two theories can be
incorporated in this combined approach, making it very promising for use in constitutive
modeling or RC structures.
1.1 Scope and Objectives
It is the purpose of this study to introduce a thermodynamically consistent model for
the elasto-plastic-damage NFEA of RC beams. Developing a better understanding and
representation of the behavior of RC beam structures subjected to short term static
loading intensity will be the primary goal. The loading regime will be such that the
rotation of the direction of the principal strain vector remains moderate, allowing the use
of small strain theory in the Finite Elements (FE) simulations. The mechanical behavior
of RC will be studied while isothermal conditions are assumed throughout this work.
The model will be translated into algorithms that will simulate the nonlinear material
behavior of concrete, steel reinforcement, and their interaction. They will also facilitate
the reproduction of experimentally observed load carrying capacity curves of RC beams.
These algorithms will be incorporated into the FEA software ABAQUS via a user-
defined material subroutine (UMAT). While ABAQUS performs the standard FE
procedure using standard types of finite elements, the UMAT will govern the behavior of
these materials during different loading stages, i.e., elastic, inelastic, failure, post-failure
loads. The model is intended to be robust, efficient and reliable but by no means vague
and complicated.
The complexity of the behavior of RC, and the scatter of the experimental data
associated with machine precision, variations in testing techniques, and statistical
distributions of material properties from one sample to another is one of the main factors
enforcing the notion that the primary goal of any constitutive model should be set in the
prediction of essential features of experimentally observed behavior, rather than in
exactly replicating the entire history of stress-strain curves. Along this line, it should be
emphasized that numerical implementation of a proposed constitutive model into a
computer code is almost as important an issue to consider as the model itself. A literature
survey can easily reveal models that are mathematically very elegant, but pose
overwhelming computational difficulties. It is thus important that a constitutive model,
although rigorous in theory, should also be suitable for use in computation and should
lend itself well to an efficient implementation in computer codes.
The objective of this study can be cast into the following points:
• Capturing the elasto-plastic-damage behavior of concrete under monotonic loads
in tension, compression, tension-tension, compression-compression, and
compression-tension stress states. This incorporates defining the continuous
damage mechanism in concrete that will represent the strain softening in the post-
peak regions and the degradation of elastic stiffness.
• Modeling the behavior of structural steel - embedded in the concrete - as an
elastoplastic material with isotropic hardening.
4
• Accounting for the effect of bond between steel and concrete and their interaction
on the overall short term material (stress-strain) behavior of RC.
• Verifying the applicability of the proposed model by comparing the predicted
behaviors with those observed in experimental results obtained by other
researchers.
1.2 Outline of the Dissertation
This study starts with a literature review on the physical behaviors of RC and its
macroscopic constituents, i.e., steel and concrete, as well as the constitutive models
applied to describe such behaviors, all of which are presented in Chapter 2. Emphasis
here is placed on the experimental work done by various researchers to study concrete
and RC in structural members and on the constitutive models applied to capture the
important features of the experimentally observed material behaviors. This emphasis
leads to the application of combined approaches employing both damage mechanics and
plasticity theory to describe the concrete behavior and classical hardening plasticity
theory to describe that of steel reinforcement. The bond effect and material interactions as
well as modeling of these phenomena are also addressed.
Chapter 3 addresses the von Mises plasticity-based material constitutive models used
to describe the physical behavior of steel reinforcement. This includes perfectly plastic
and (linear, bilinear and nonlinear) hardening models. The chapter starts with an
introduction discussing the application of the FE method to nonlinear continuum
mechanics. The fundamentals of FE analysis procedures for elastic-plastic problems are
considered, with algorithms required for convergence at the local (Gauss point) and
global (equilibrium) levels. This is followed by a discussion of the incremental
constitutive theory of metal plasticity and the numerical techniques employed to integrate
the constitutive equations. Examples are provided to verify the implemented algorithms.
This chapter is intended to serve as an introduction to the chapters that follow in terms of
theoretical formulation and computational implementation methods.
In Chapter 4, the material model for plain concrete is discussed. Rigorous consistent
thermodynamic formulation is employed to derive the framework of the elastic-plastic-
damage constitutive model. The Helmholtz free energy function is discussed and the
dissipation potentials for plastic and damage processes are postulated. To capture the
different responses in tension and compression, the approach makes use of the separation
of tensile and compressive behaviors, obtained through the decomposition of stress tensor
and the introduction of the damage variables, all of which are integrated into the
thermodynamic framework. The dissipation process therefore consists of three separate
dissipation mechanisms: tensile and compressive damage coupled with plasticity.
Detailed description of the theoretical formulation of the constitutive model in terms of
the theories of plasticity with multiple hardening rules and continuum damage mechanics
is then provided, followed by a step-by-step detailing of the numerical scheme applied to
integrate the incremental constitutive equations. Non associative plasticity - with multiple
hardening rules - is combined with continuum damage mechanics, where the accumulated
damage in the concrete material is represented by two internal damage parameters, one in
5
tension and one in compression. The constitutive model is formulated as a relation
between the undamaged stress and an internal accumulated damage parameter defined
using the tensile and compressive damage parameters. The integration scheme of elastic
predictor followed by plastic and damage correctors is adopted in order to overcome the
difficulties that arise from the combined plastic-damage approach. Several verification
examples are provided in order to test the model’s predictions under uniaxial and biaxial
stress states, as well as under three point bending test of a single-edge-notched concrete
beam. The numerical results are compared to the experimental ones in order to evaluate
the performance of the proposed concrete model and to asses the model’s ability to
capture the experimentally observed behaviors of concrete materials.
Chapter 5 addresses modeling of RC beams. The bond effect on RC is discussed and
a simple experimentally based methodology is applied to account for the effect of steel-
concrete interaction on the reinforcement stress level. This discussion is followed by two
RC comprehensive examples provided to demonstrate the applicability of the applied
constitutive models combined to the NFEA of RC beams. The model is tested using two
RC beams with different geometries and reinforcement ratios. The predicted (numerical)
results are compared to their experimental counterparts obtained by other researchers.
Concrete damage distributions as well as RC stress-strain curves are presented and
compared to other works.
Conclusions and further studies are proposed in Chapter 6. A summary of the
constitutive model and results is introduced followed by a brief discussion of the merits
and weaknesses of the currently proposed model. Results of the current investigation
suggest that additional fundamental research is required if computer simulation is to be a
viable tool for future research and design of RC structures. A discussion of additional
research needs in the area of characterization of material response through experimental
investigation, material modeling and non-linear analysis is presented.
Throughout this work, (+) and (
−
) superscripts are used to indicate tension or
compression, respectively. All symbols with an over bar, e.g.
σ
, are considered to be in
the effective fictitious undamaged configuration. On the other hand, the absence of the
over bar indicates the actual damaged configuration.
6
CHAPTER 2
LITERATURE REVIEW
RC structures are made up of two materials with different characteristics, namely,
concrete and steel. Steel can be considered as a homogeneous material with generally
well defined material properties. Concrete, on the other hand, is a heterogeneous material
made up of cement, mortar and aggregates. Its mechanical properties are widely scattered
and cannot be defined easily. For the convenience of analysis and design, however,
concrete is often considered a homogeneous material at the macroscopic scale.
The typical stages in the load-deformation behavior of a RC simply supported beam
are illustrated in Fig. 2.1. Similar relations are obtained for other types of RC structural
elements. This 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 (Chen, 1982).
Figure 2.1 Typical load-displacement curve of RC beams
The nonlinear response is caused by three major effects, namely, cracking of concrete
in tension, yielding of the reinforcement or crushing of concrete in compression, and the
interaction of the constituents of RC. Interaction includes 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
7
on environmental conditions, such as ambient temperature and humidity. The material
properties of concrete and steel are also strain-rate dependent to some extent.
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 to determine principal stresses in concrete, stresses in
steel reinforcement and bond stresses. Ngo and Scordelis (1967) reported that one of the
main difficulties in constructing an analytical model for RC member is due to the
composite action of steel and concrete. Prefect bonding between steel and concrete can
only exist at an early stage under low load intensity. As the load is increased, cracking as
well as breaking of bond inevitably occurs, and a certain amount of bond slip will take
place in the beam, all of which will in turn affect the stress distributions in concrete and
steel.
Since Ngo and Scordelis published their landmark paper in 1967, the analysis of RC
structures has enjoyed a growing interest and many publications have appeared.
Important progress has also been made in the finite-element-based numerical analysis of
plain and RC structures (recent examples: Fantilli et. al., 2002; Sumarac et. al., 2003;
Marfia et. al., 2004; Phuvoravan and Sotelino, 2005; Junior and Venturini, 2007; just to
mention a few). However, despite this progress, modeling the mechanical behavior of
concrete is still one of the most difficult challenges in the field of structural engineering.
This is due to the inherent complexity and uncertainty concerning the properties of
concrete, which make it excessively difficult to develop accurate constitutive models or
algorithms that are sufficiently robust to obtain reliable and converged solutions in
numerical analyses.
An adequate numerical analysis of the nonlinear behavior of RC structures is based on
the coupled modeling of different inelastic processes in concrete and in reinforcement.
Macroscopic representation of crystal dislocation in steel reinforcement within the
framework of elastoplasticity yields a reliable prediction of deformation history of
reinforcement. The realistic constitutive behavior of concrete is, however, more complex.
It has resulted in the appearance of many different concepts to its theoretical description.
Concrete failure is usually characterized by many macroscopic cracks. If the discrete
cracks are considered, the necessary adaptation of the FE mesh to the trajectory of each
crack under new load step makes the FE analysis cumbersome. Even in the two-
dimensional case, it limits the applicability of conventional fracture mechanics to simple
specimens with one or two cracks. Therefore, the use of models based on continuum
damage mechanics and elastoplasticity have found large application in the numerical
modeling of concrete fracture.
The development of analytical models for the response of RC structures is
complicated due to the following factors (Kwak and Fillipou, 1997):
8
• RC is a composite material made up of concrete and steel, two materials with very
different physical and mechanical behavior.
• Concrete exhibits nonlinearities even under low level of 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.
Because of these factors and the differences in short- and long-term behaviors of the
constituent materials, it is common practice among researches to model the short- and
long-term response of RC members and structures based on separate material models for
reinforcing steel and concrete, which are then combined along with models of interaction
between the two constituents to describe the behavior of the composite RC material. This
will be the approach adopted in this study of short-term behavior of RC beams. In the
following, literature review of the research done to describe the behavior of concrete,
reinforcing steel, and their bond-interaction analysis is presented.
2.1 Concrete Materials
Concrete by itself is a composite material. It is made of cement, mortar, and
aggregates. The thermo-chemical interaction between these constituents results in a
unique building material. One of the most important characteristics of concrete is low
tensile strength, which results in tensile cracking at a very low stress compared with
compressive stresses. The tensile cracking reduces the stiffness of the concrete
component.
Several issues in the current practice of constitutive modeling of concrete material
need to be addressed. First, concrete is a non-homogeneous and anisotropic material, the
mechanical behavior of which is nonlinear (Kupfer et. al., 1969; Kotsovos and Newman,
1977). Its compressive strength increases as it is loaded in a biaxial compressive state, but
decreases as the tensile stress is increased under biaxial compression–tension. Moreover,
the ductility of concrete under biaxial stresses is also dependent on the stress state.
Concrete exhibits a large number of micro-cracks, especially at the interface between
coarser aggregates and mortar, even before the application of any external loads. The
presence of these micro-cracks has a great effect on the mechanical behavior of concrete,
since their propagation (concrete damage) during loading contributes to the nonlinear
behavior at low stress levels and causes volume expansion near failure. Many of these
micro-cracks are initially caused by segregation, shrinkage or thermal expansion of the
mortar. Some micro-cracks 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 the mortar; it constitutes the weakest link in the
composite system. This is the primary reason for the low tensile strength of concrete.
Furthermore, the relation between the microstructure and mechanical behavior of
concrete is quite complex because of the considerable heterogeneity of the distinct phases
9
of the material. Concrete may be treated as a composite material and it may contain
porosity in its matrix. The porosity in the matrix is not homogenous and a strong porosity
gradient is observed around the inclusions formed by the aggregates (Panoskaltsis and
Lubliner, 1994; Ollivier et. al., 1995). This area in the matrix affected by the surface of
the aggregate is known as “transition zone”. This transition zone has a damaging effect
on the mechanical behavior of concrete. It is clear from the above discussion that the size
and texture properties of the aggregates will also have a significant effect on the
mechanical behavior of concrete under various types of loading (Chen, 1982).
The nonlinear material behavior of concrete can be attributed to two distinct material
mechanical processes; plasticity and damage mechanisms. The cracking process in
concrete is distinguished from cracking of other materials in that it is not a sudden onset
of new free surfaces but a continuous forming and connecting of micro-cracks. The
formation of micro-cracks is presented macroscopically as softening behavior of the
material, which causes localization and redistribution of strains in the structure. This
phenomenological behavior at the macroscopic level can be modeled by classical
plasticity (Pramono and Willam, 1989). On the other hand, the micro-processes such as,
micro-cracking, micro-cavities, and their nucleation and coalescence, also cause stiffness
degradation, which is difficult to represent with classical plasticity (Lee and Fenves,
1998). This introduces the need for continuum damage mechanics, where the stiffness
degradation can be modeled by making use of the effective stress concept (Kachanov,
1958) as will be shown later. Damage mechanics can also be used to represent the post-
peak softening behavior of concrete materials; a behavior that cannot be addressed by
classical plasticity theory.
The failure behavior of concrete is governed by complex degradation processes
starting within the aggregate-matrix interface. These processes are shown in Fig. 2.2. The
aggregate-matrix interface contains micro-cracks that exist even before any load is
applied to concrete (Fig. 2.2a). The formation of these micro-cracks is primarily due to
stress and strain concentrations resulting from the incompatibility of the elastic moduli of
the aggregate and cement paste components. Strain concentrations at the aggregate-
mortar interface may also occur as a result of volume changes in the concrete due to
shrinkage and/or thermal effects resulting from the difference in thermal coefficients of
various constituents. Additional micro-cracks can be initiated when concrete is subjected
to external loads with magnitudes beyond the micro-crack initiation threshold (Fig. 2.2b).
These micro-cracks spread and grow under the effect of continuous loading until they
merge into the matrix after a certain threshold is reached (Fig. 2.2c). Cracks in the matrix
grow in size and coalesce with each other to form major cracks that eventually lead to
failure (Fig. 2.2d).
The nonlinear stress-strain behavior under uniaxial compression of concrete is shown
in Fig. 2.3. Investigators (e.g. Kotsovos and Newman, 1977) have shown that concrete
compression behavior and fracture characteristics may be explained by the creation and
propagation of micro-cracks inside the concrete. It is observed that under different
magnitudes of the applied load, the concrete behavior can be summarized in four stages
shown in Fig. 2.3. The first stage is observed during 30-60% of the ultimate strength
10
(shown as 45% in Fig. 2.3). In this initial stage, one can observe the highest tensile strain
concentration at particular points where new micro-cracks are initiated as shown in Fig.
2.2b. At this load state, localized cracks are initiated, but they do not propagate
(stationary cracks). Hence, the stress-strain behavior is linearly elastic. Therefore, 0.3
'
c
f
is usually proposed as the limit of elasticity. Beyond this limit, the stress-strain curve
begins to deviate from a straight line. Stresses up to 70-90% of the ultimate strength
(shown as 85% in Fig. 2.3) characterize the second stage. In this stage, as the applied load
is progressively increased, the crack system multiplies and propagates as shown in Fig.
2.2c. The increase of the internal damage in concrete, revealed by deviation from the
linear elastic behavior, reduces the material stiffness and causes irrecoverable
deformation in unloading. Although the relief of strain concentration continues during
this stage, void formation causes the rate of increase of the tensile strain in the direction
normal to that of branching to increase with respect to the rate of increase of the strain in
the direction of branching (Kotsovos and Newman, 1977). The start of such deformation
behavior is called “onset of stable fracture propagation” (OSFP). In this load stage, the
mortar cracks tend to bridge the aggregate-matrix bond cracks.
Figure 2.2 Aggregate-matrix interface: a) prior to loading, b) 65% of ultimate load, c)
85% of ultimate load, d) failure load, (Kotsovos and Newman, 1977)
A third stage shown in Fig. 2.3 extends up to the ultimate strength. Interface micro-
cracks are linked to each other by mortar cracks as shown in Fig. 2.2c, and void
formation (dilation) begins to have its effect on deformation at this stage. The start of this
stage is called “onset of unstable fracture propagation” (OUFP). This level is easily
11
defined since it coincides with the level at which the overall volume of the material
becomes a minimum. In this stage, the progressive failure of concrete is primarily caused
by cracks through the mortar. These cracks merge with bond cracks at the surface of
nearby aggregates and form crack zones of internal damage. Following that, a smoothly
varying deformation pattern may change and further deformation may be localized.
A fourth stage defines the region beyond the ultimate strength. In this region, the
energy released by the propagation of a crack is greater than the energy needed for
propagation. Therefore, the cracks become unstable and self-propagating until complete
disruption and failure occurs. In this stage, the major cracks form parallel to the direction
of the applied load, causing failure of the concrete. The volume of voids increases
dramatically causing a rapid dilation of the overall volume of concrete as shown in Fig.
2.2d.
Figure 2.3 Uniaxial compression stress-strain relation of concrete, (Chen, 1982).
All the above mentioned stages are for the uniaxial compression case. Stages I, (II and
III), and IV could be categorized into the linear elastic, inelastic, and the localized stages
respectively. Understanding these stages is crucial for the development of any concrete
model.
Figure 2.4 shows a typical uniaxial tension stress-elongation curve. In general the
limit of elasticity is observed to be about 60-80% of the ultimate tensile strength. Above
this level, the aggregate-matrix interface micro-cracks start to grow. As the uniaxial
tension state of stress tends to arrest the cracks much less frequently than the compressive
stage of stress, one can expect the interval of stable crack propagation to be quite short,
and the unstable crack propagation to occur much earlier. That is why the deformation
behavior of concrete in tension is quite brittle in nature. In addition, the aggregate-matrix
12
interface has a significantly lower tensile strength than the matrix, which is the primary
reason for the low tensile strength of concrete.
Nonlinearities in concrete behavior are well documented and arise from two distinct
micro-structural changes that take place in the material: one is the plastic flow; the other
is the development of micro-cracks and micro-voids. From a plasticity point of view, the
number of bonds between atoms during the plastic flow process is hardly altered;
therefore, the elastic compliances remain insensitive to this mode of micro-structural
change. On the other hand, micro-cracking destroys the bond between material grains,
affects the elastic properties, and may also result in permanent deformations, which can
be modeled by damage mechanics.
Figure 2.4 Uniaxial tensile stress-strain curve (Peterson, 1981)
The idea of combined plasticity and damage mechanics theories through the
description of plasticity and damage surfaces has been explored and used in the past
(Oritz, 1985; Lubliner et. al. 1989). Many researches attempted to expand the application
of plasticity and damage theories to concrete (Chen and Chen, 1975; Lubliner et. al.,
1989; Yazdani and Schreyer, 1990; Abu-Lebdeh and Voyiadjis, 1993; Voyiadjis and
Abu-Lebdeh, 1994; Luccioni et. al., 1996; Lee and Fenves, 1998, 2001; Faria et. al.,
1998; Hansen et. al., 2001; Nechnech et. al., 2002; Gatuingt and Pijaudier-Cabot, 2002;
Kratzig and Polling, 2004; Salari et. al., 2004; Shen et. al., 2004; Jankowiak and
Lodygowski, 2004; Luccioni and Rougier, 2005, Rabczuk et. al., 2005; Wu et. al., 2006;
Grassl and Jirasek, 2006; Nguyen and Korsunsky, 2006; Shao et. al., 2006; Jason et. al.,
2006; Contrafatto and Cuomo, 2006; Nguyen and Houlsby, 2007; Mohamad-Hussein and
Shao, 2007; Ananiev and Ožbolt, 2007; Voyiadjis et. al., 2008a,b; Yu et. al., 2008; and
others).
Plasticity by itself fails to address the softening behavior of concrete under tension and
compression caused by damage propagation due to micro-cracking in the strained
13
material. On the other hand, damage mechanics is only concerned with the description of
this progressive weakening of solids due to the development of micro-cracks and micro-
voids (Loland, 1980; Oritz and Popov, 1982; Krajcinovic, 1985; Simo and Ju, 1987a,
1987b; Ju et. al., 1989; Voyiadjis and Kattan, 1989, 1999, 2006). Therefore, a
constitutive model should equally address these two distinct physical modes of
irreversible changes and should satisfy the basic postulates of mechanics and
thermodynamics governing these phenomena. Concrete plasticity and damage will be
further discussed in what follows.
2.1.1 Concrete Plasticity
Plasticity theory successfully treated concrete problems in which the material is
subjected to primary compressive loads. In situations where tension-compression plays a
significant role, plasticity theory is applied to model the compression zones while
damage or fracture mechanics is used to model the tensile zones (Lubliner et. al., 1989).
The uniaxial behaviors of plain concrete under tension and compression up to tensile
and compressive failure are shown in Fig. 2.5. For tensile failure, the behavior is
essentially linear elastic up to the failure load, followed by a strain softening response
which was generally neglected or idealized in the past. For compression failure, the
material initially exhibits almost linear behavior up to the proportional limit at point A,
after which the material is progressively weakened by internal micro-cracking up to the
end of the perfectly plastic flow region CD at point D. The nonlinear deformations are
basically plastic, since upon unloading only the elastic portion
e
ε
can be recovered from
the total deformation
. It is clear that the phenomenon in regions AC and CD
corresponds exactly to the behavior of a work-hardening elastoplastic and elastic
perfectly plastic solid, respectively. As can be seen from Fig. 2.5, the total strain
in a
plastic material can be considered as the sum of the reversible elastic strain
e
ε
and the
permanent plastic strain
p
ε
. A material is called perfectly plastic or work-hardening
according to whether it does or does not admit changes of permanent stain under constant
stress (Chen, 1982).
The plastic response of concrete exhibits some characteristics that the classical theory
of plasticity can not describe. It was shown experimentally that there is a lack in
simulating the normality rule (Adenaes et. al., 1977). Furthermore, the descending branch
of the uniaxial stress-strain diagram of concrete resembles a violation of the Drucker’s
stability postulate. On the other hand, and from a macroscopic point of view; classical
plasticity can simulate the behavior of concrete particularly in the pre-peak regime such
as the nonlinearity of the stress-strain curve and the significant irreversible strain upon
loading. Therefore, the plasticity theory can be used in the modeling of strain hardening
behavior of concrete. Many works were presented by researchers to modify the classical
theory of plasticity in order to make it more suitable for concrete materials (Feenstra and
de Borst, 1996; Bicanic and Pearce, 1996; Grassl et. al., 2002; Park and Kim, 2005; and
others).
14
Figure 2.5 Uniaxial stress-strain curve, pre- and post failure regimes, (Chen, 1982)
The main characteristics of the plasticity models (Chen and Schnobrich, 1981; Han
and Chen, 1985; Oritz, 1985; Lubliner et. al., 1989; Voyiadjis and Abu-Lebdeh, 1994;
Lee and Fenves, 1998, 2001; Grassl and Jirásek, 2006; and others) used to describe the
behavior of concrete include: pressure and path sensitivity, non-associative flow rules,
work or strain hardening, and limited tensile strength. Many of those models have been
developed for the finite element analysis of structural elements. Some of theses models
are associated with high mathematical complexity which renders them undesirable for
many engineering applications, especially the analysis and design of simple structural
elements, such as beams and columns.
A drawback of the plasticity-based approach is that the stiffness degradation due to
progressive damage is not modeled. However, some researcher advocated that
experimental evidence (Willam et. al., 1986; Hordijk, 1991) shows that the stiffness
degradation due to tensile cracking is substantial only when the tensile cracking has
developed fully, and the stiffness degradation due to compressive loading is even less
pronounced than the stiffness degradation due to tensile loading. Therefore, they argued
that under monotonic loading where only local unloading occurs, neglecting the
degradation of the elastic stiffness does not seem to entail major errors (Feenstra and de
Borst, 1996).
2.1.2 Concrete Damage Mechanics
Concrete contains numerous micro-cracks even before the application of any external
loads. These inherent micro-cracks are mainly present at the aggregate-cement interface
as a result of shrinkage, and thermal expansion in the cement paste or aggregates.
Damage in concrete is primarily caused by the propagation and coalescence of these
micro-cracks. The growth of these micro-cracks during loading causes reduction in
strength and deterioration in the mechanical properties of the concrete material.
15
Therefore, modeling of crack initiation and propagation is very important in the failure
analysis of concrete structures.
To model damage in concrete, various types of constitutive laws have been presented
including different approaches such as the endochronic theory (Bazant, 1978) the plastic
fracturing theory (Dragon and Mroz, 1979) the total strain models (Kotsovos 1980),
plasticity with decreasing yield limit (Wastiels 1980), microplane models (Bazant and
Ozbolt, 1990), and the bounding surface concept (Voyiadjis and Abu-Lebdeh, 1993).
Continuum damage mechanics was first applied to metals and later modified to model
different materials. The term damage mechanics has been conventionally used to refer to
models that are characterized by a loss of stiffness or a reduction of the secant
constitutive modulus. It was first introduced by Kachanov (1958) for creep-related
problems and then later applied to the description of progressive failure of metals and
composites and to represent the material behavior under fatigue (Kachanov, 1986). The
use of continuum damage mechanics in concrete began in 1980’s. Damage models were
used to describe the strain-softening behavior of concrete. Since then, researchers
developed different damage models to represent concrete (Lemaitre and Mazars, 1982;
Krajcinovic, 1986; Mazars and Pijaudier-Cabot, 1989; Voyiadjis and Abu Al-Lebdeh,
1992; and others).
Damage theory in concrete materials can represent the post-peak region (Krajcinovic,
1983b; Mazars and Pijaudier-Cabot, 1989; and others.). The use of fracture mechanics in
concrete materials, on the other hand, was debated extensively (Mindess, 1983;
Krajcinovic and Fanella, 1986). The debate of the use of fracture mechanics for concrete
materials was supported by some experimental evidence (Pak and Trapeznikov, 1981),
which showed that energy dissipation and planar crack idealization in concrete make the
application of fracture mechanics in concrete materials very arguable. Therefore,
continuum damage mechanics is applied for the determination of macroscopic damaging
variables and material properties (Krajcinovic, 1979; Dragon and Mroz, 1979; among
others).
There are several approaches of how the stiffness degradation is included in a model.
Some researchers coupled damage with elastic analysis only (Mazars, 1984; Willam et.
al., 2001; Tao and Phillips, 2005; Labadi and Hannachi, 2005; Junior and Venturini,
2007; Khan et. al., 2007), while others coupled damage with plasticity. In the plastic-
damage model (Simo and Ju, 1987a,b; Ju, 1989; Lubliner et. al., 1989; Voyiadjis and
Kattan, 1989; Luccioni et. al., 1996; Armero and Oller, 2000; Salari et. al., 2004;
Voyiadjis et. al., 2008a) stiffness degradation is embedded in a plasticity model. These
models introduce the use of elastic and plastic damage variables to represent the stiffness
degradation. The damage variables are coupled with the plastic deformation in the
constitutive formulations which provide help for calibrating the parameters with the
experimental results. Yet, the coupled relations are complex and result in an unstable
numerical algorithm. This kind of algorithm may cause unrealistic representation of the
plastic behavior of the concrete during numerical implementation and iteration
procedures (Lee and Fenves, 1994).
16
In the effective-plasticity and damage model, plasticity is formulated in the effective
stress space to model the plastic irreversible phenomena while continuum damage
mechanics is adopted to represent stiffness degradation (Simo and Ju, 1987a,b; Ju, 1989;
Hansen and Schreyer, 1992; Yazdani and Schreyer, 1990; Fichant et. al., 1999; Voyiadjis
and Kattan, 1999; Mazars and Pijaudier-Cabot, 1989; Salari et. al., 2004; Shen et. al.,
2004; Lee and Fenves, 1998; Faria et. al., 1998; Jefferson, 2003; Voyiadjis and Kattan,
2006; Jason et. al., 2006; Voyiadjis et. al., 2008b). It has the advantage to decouple the
stiffness degradation from the plastic deformation by linearizing the evolution equations
(e.g., Jason et. al. 2006).
Others considered the interaction of thermal and mechanical damage processes in
concrete and concrete-like heterogeneous materials (Nechnech et. al., 2002; Willam et.
al., 2003; and others). They studied the interaction of thermal and mechanical damage at
the mesomechanical level of observations, when volumetric and deviatoric degradation
take place simultaneously. They also studied the effect of thermal expansion and
shrinkage in the two-phase concrete material when thermal softening of the elastic
properties leads to massive degradation of the load resistance.
Quasi-brittle material under different kinds of loading undergoes several damage
states, such as tensile cracking, compressive failure, and stiffness degradation. Some
researchers recently accounted for all of these effects in a concrete model using a single
(scalar) damage variable
Ï•
(Luccioni et. al., 1996; Lee and Fenves, 1998, 2001; Willam
et. al., 2001; Ferrara and di Prisco, 2001; Nechnech et. al., 2002; Soh et. al., 2003; Shen
et. al., 2004; Jankowiak and Lodygowski, 2004; Salari et. al., 2004; Luccioni and
Rougier, 2005; Labadi and Hannachi, 2005; Kukla et. al., 2005; Shao et. al., 2006; Jason
et. al., 2006; Grassl and Jirásek, 2006; Nguyen and Korsunsky, 2006; Junior and
Venturini, 2007; Sain and Kishen, 2007; Mohamad-Hussein and Shao, 2007; and others).
In order to account for different responses of concrete under various loadings, other
researchers used multiple hardening scalar damage variables
Ï•
±
(Mazars, 1986; Mazars
and Pijaudier-Cabot, 1989; Faria et. al., 1998; Comi and Perego, 2001; Gatuingt and
Pijaudier-Cabot, 2002; Willam et. al., 2003; Jirásek, 2004; Marfia et. al., 2004; Tao and
Phillips, 2005; Wu et. al., 2006; Contrafatto and Cuomo, 2006; Ananiev and Ožbolt,
2007; Nguyen and Houlsby 2008a,b; and others). It was argued that isotropic continuum
damage mechanics models with multiple damage variables cannot represent the entire
spectrum of damage effects on biaxial tensile and compressive strength of the material
because the damage variables eventually contribute to the same isotropic evolution in
both strengths. This introduced the need for concrete anisotropic damage models. Yet
isotropic damage is very widely used in different types of combinations with plasticity
models due to the simplicity of its implementation in a material model for concrete
material.
There are models that present thermodynamic theories of anisotropic damage mainly
by the use of tensor damage variables (Krajcinovic and Lemaitre, 1986; Krajcinovic and
Fonseka, 1981; Krajcinovic, 1983b; Voyiadjis and Kattan, 1989, 1999, 2006). Chow and
his co-workers (Chow and Wang, 1988; Chow and Lu, 1989) proposed an energy-based
elastic-plastic damage model in order to describe the difference in the observed failure
17
modes of geological materials under compression and tension, by using a damage tensor
identified by the fourth rank order and fourth rank projection tensors. The coupling
between elastic-plastic deformation and anisotropic damage using symmetric second
order damage tensors is widely presented by several authors (Simo and Ju, 1987a,b; Ju,
1989; Voyiadjis and Kattan, 1989, 1990, 1992a,b; Hansen et. al., 2001; Tikhomirov and
Stein, 2001; Voyiadjis et. al., 2008a,b; and others).
In the case of anisotropic damage, different levels of damage are related to different
principal directions. This can be done using a symmetric second-order damage tensor,
ij
Ï•
(Murakami and Ohno, 1981; Murakami, 1983; Ortiz, 1985; Murakami, 1988; Voyiadjis
and Kattan, 1992a,b, 1999, 2006; Voyiadjis and Abu-Lebdeh, 1993; Voyiadjis and
Venson, 1995; Voyiadjis and Park, 1997, 1999; Voyiadjis and Deliktas, 2000, Voyiadjis
et. al., 2001, Hansen et. al., 2001; Tikhomirov and Stein, 2001;Voyiadjis et. al., 2008a,b;
and others), or fourth order damage tensor
ijkl
Ï•
(Chaboche, 1993, Chaboche et. al., 1995,
Taqieddin et. al., 2005, 2006; Voyiadjis et. al., 2007a,b; and others).
The variety in all of these anisotropic models is somewhat puzzling because a) the
relation between these models is difficult to establish (except maybe in the case of the
isotropic damage) and b) the comparison of damage-induced anisotropy with
experimental data is difficult and therefore the characterization of the damage-induced
anisotropy of the material requires three-dimensional experimental facilities. Because of
this difficulty in the experimental simulation of concrete specimens, together with the
inconvenience related to computational aspects, many researchers avoided these
complexities by resorting to simplified (single or multiple) isotropic damage models.
Concrete damage can be characterized as a reduction in the material stiffness. As
shown in Fig. 2.6 and 2.7, stiffness reduction in tensile loading is less than that under
Normalized strain, ε
nom
= ε / ε
c0
Normalized stress f
nom
= f / f
c0
Figure 2.6 Concrete response to monotonic and cyclic compression load
(data from Bahn and Hsu, 1998)
18
compressive loading and therefore, more damage occurs in tension than in compression.
Thus, the characterization of material damage in tension and in compression should be
considered when modeling the response of plain concrete, justifying the use of multiple
damage variables.
Development of a damage-based model requires the definition of a damage rule that
characterizes the rate at which material damage is accumulated. The identification of this
damage rule may also include the definition of a damage surface that defines an initial
elastic domain. Various proposed damage models differ in the definition of the damage
surface and the corresponding damage rules. Some of the first constitutive relationships
for damage characterization were proposed for the isotropic damage case based on the
assumption that one defines an effective stress that is larger than the Cauchy stress and
accounts for the reduction in the material’s area that results from micro-cracking
(Kachanov, 1958).
In order to solve the problem of mesh non-objectivity (sensitivity) in the FE analysis,
which is often encountered when using models that are based on the strength theory, the
fracture energy method is often incorporated in damage mechanics to give an improved
mesh-size independent description of the post-peak behavior of concrete (Faria et. al.,
1998; Lee and Fenves, 1998, 2001; Salari et. al., 2004; Kratzig and Polling, 2004;
Rabczuk et. al., 2005; Luccioni and Rougier, 2005; He et. al., 2006; Contrafatto and
Cuomo, 2006; Wu et. al., 2006; Nguyen and Houlsby 2008a,b; and others). The concept
of energy dissipation, i.e. fracture energy, has been used extensively in finite element
calculations and can be considered as accepted in the engineering community (Feenstra
and de Borst, 1996). With the assumption that the fracture energy is uniformly dissipated
in a representative area of the structure, the equivalent length, the finite element
calculations will lead to results that are insensitive with regard to the global structural
response upon mesh refinement, at least below a certain level of refinement.
Axial deformation, µm
Axial stress (MPa)
Figure 2.7 Stress-deformation history for concrete subjected to cyclic tensile
loading (data from Reinhardt, 1984)
19
2.2 Steel Reinforcement
Reinforcement comes in different types and shapes. Those most commonly used are
the deformed circular cross-sectional bars. The spiral deformation pattern on the bars
strengthens the mechanical bond between the bars and concrete. 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 RC structures.
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 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 (Wang and Salmon,
1998).
Two different idealizations, shown in Fig. 2.8, are commonly used depending on the
desired level of accuracy (ASCE 1982). The first idealization neglects the strength
increase due to strain hardening and the reinforcing steel is modeled as a linear elastic,
perfectly plastic material, as shown in 3.8a (Mazars and Pijaudier-Cabot, 1989; Sumarac
et. al., 2003; Kratzig and Polling, 2004; Luccioni and Rougier, 2005; and Wu et. al.,
2006). 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.
Figure 2.8 Idealizations of the steel stress-strain curves.
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
20
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 Fig. 2.8b for the case of bilinear stress-
strain models (Salari and Spacone, 2001; Limkatanyu and Spacone, 2003). 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 (Fig. 2.9). These
parameters can be derived from experimentally obtained stress-strain relations.
Figure 2.9 Linear elastic, linear strain hardening steel stress-strain relation
In this study the reinforcing steel is modeled as a linear elastic, linear strain hardening
material with yield stress
y
σ
, as shown in Fig. 2.9. 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 Fig. 2.8a leads to numerical
convergence problems near the ultimate member strength. It is, therefore,
advisable to take advantage of the strain-hardening behavior of steel to improve
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 (Fig. 2.9). Such a model has been
successfully used for the analyses of RC structures (Ngo and Scordelis, 1967;
Feenstra, 1993; Kwak and Fillippou, 1997; Lowes, 1999; Tikhomirov and Stein,
21
2001; Rabczuk et. al., 2005; Phuvoravan and Sotelino, 2005; He et. al., 2006;
Junior and Venturini, 2007, and others).
Different approaches can be used for modeling reinforcing steel in RC structure using
the finite element method. Some researchers use two dimensional elements to represent
steel reinforcement in their finite element meshes (Ngo and Scordelis, 1967; Mazars and
Pijaudier-Cabot, 1989; Pijaudier-Cabot et. al., 1991; Luccioni and Rougier, 2005; and
Wu et. al., 2006; and others). On the other hand, one dimensional representation of steel
reinforcement is more widely used since it is unnecessary to introduce the complexities
of multiaxial constitutive relationships for structural steel (Chen, 1982; Feenstra, 1993;
Kwak and Fillippou, 1997; Faria et. al., 1998; Lowes, 1999; Soh et. al., 2003; Sumarac et.
al., 2003; Kratzig and Polling, 2004; Rabczuk et. al., 2005; Phuvoravan and Sotelino,
2005; He et. al., 2006; and others). The discrete model, the distributed or smeared model,
and the embedded model are examples of one dimensional representation of steel
reinforcement in a RC finite element mesh. The most widely used approach is the
discrete one dimensional truss element, which is assumed to be pin connected to the
concrete elements and posses two degrees of freedom at each node (Feenstra, 1993; Faria
et. al., 1998; Kratzig and Polling, 2004; Phuvoravan and Sotelino, 2005; He et. al., 2006;
and others). Alternatively, beam elements can also be used where three degrees of
freedom are allowed at each end of the bar element (Rabczuk et. al., 2005). In the
distributed model, the reinforcing steel is assumed to be distributed over the concrete
elements at a certain orientation angle, and in the embedded steel model the reinforcing
steel is considered as an axial member built into the isoparametric elements representing
the concrete. A significant advantage of the discrete representation, in addition to the
simplicity, is that it can include the slip of reinforcing steel with respect to the
surrounding concrete in contrast with the other 1D models where perfect bond between
concrete and steel is usually assumed (e.g. de Borst et. al., 2004).
2.3 Bond and Interaction at Steel – Concrete Interface
Utilization of RC as a structural material is derived from the combination of concrete
and reinforcing steel into structural elements. Concrete is strong and relatively durable in
compression while reinforcing steel is strong and ductile in tension and in compression.
Maintaining this composite action requires transfer of load between concrete and steel.
This load transfer is referred to as bond and is idealized as a continuous stress field that
develops in the vicinity of the steel-concrete interface. For RC structures subjected to
moderate loading, the bond stress capacity of the system exceeds the demand allowing
steel and the surrounding concrete to behave as a unit (full bond). However, further
application of external loads results in an increase in the stresses in the interface between
concrete and steel. It also results in localized bond demand that exceeds the bond
capacity, resulting in the deterioration of that capacity, a localized damage that gradually
spreads to the surrounding material, and a significant movement between the reinforcing
steel and the surrounding concrete. Therefore, the properties of the interface between
concrete and steel is very important in the analysis of RC structures.
22
The main stress transfer mechanisms between concrete and steel in RC elements are
represented by adhesion, mechanical interaction, and friction (Luccioni et. al., 2005).
Adhesion is constituted by chemical bonds and stresses developed during the curing
process of concrete. This transfer mechanism, schematically represented in Fig. 2.10a, is
prevailing in the case of bars of smooth surface and its failure is characterized by the
initiation and propagation of cracks in the concrete/steel interface. In the case of
corrugated bars, the described mechanism is secondary and the stress transfer is mainly
due to the interaction between ribs and the surrounding concrete. Adhesion is relatively
soon exhausted in the global response and consequently, the transfer force is transmitted
by friction and mechanical interaction between the ribs and the adjacent concrete.
As the force in the reinforcing bar is increased, the transfer forces are dominated by
the mechanical interaction concentrating at the faces of the ribs. In this state, the term
“adhesion stress” refers to the mean force per unit of surface. At increased loading, the
concrete begins to fail near the ribs with two different modes of failure; by failure of the
concrete adjacent to the contact area, as illustrated in Fig. 2.10b and by transverse
cracking.
Figure 2.10 a) Schematic presentation of adhesion between concrete and steel and its
failure, b) Stress concentration: Transverse cracking
Transverse cracking initiates at the ribs presenting a characteristic cone shape, also
called secondary cracking or adhesion cracking. The bond zone, constituted by a
damaged zone with finite depth surrounding the steel bar, is defined by the extension of
the transverse cracking (Fig. 2.11a). With the stable propagation of transverse cracking,
the concrete next to the bar seems to form inclined struts that are known as compression
cones (Fig. 2.11b). The adhesion stiffness is usually characterized by the stiffness of
these struts (Luccioni et. al., 2005).
When the load is further increased, radial splitting forces can be developed. This
phenomenon is due to the rotation of the inclined struts (Fig. 2.11b) that produces a
considerable radial component of the contact force, the increase of the radial force that is
produced by the increase of the effective contact angle between ribs and concrete due to
the deposition of the crushed concrete at the faces of the ribs, the wedge action of the
ribs, and the radial component of the contact forces. Without an accurate confinement, a
splitting failure can occur, spreading the effect of bond outside the bond zone.
23
Figure 2.11 a) Failure of concrete in contact with face of rib, b) Change in direction of
struts: Longitudinal cracking
The bond zone dilation, due to the longitudinal cracking, takes place when the
adhesion stress reaches values next to the limit one. Following this moment, a stress
softening is observed in the response. This stress softening, characteristic of the bond–
slip behavior, is frequently interpreted as a progressive shear failure of concrete between
the ribs. When the confinement stresses are low, the geometric variation in the contact
occurring between the ribs and concrete can also contribute to the stress softening. Both
the progressive shear failure of concrete and the contact zone geometric variation are
stimulated by the reduction of the confinement stresses produced by the propagation of
longitudinal cracking. In both cases, the stress softening reveals a strong discontinuity
between the reinforcing bar and the surrounding concrete.
Different approaches have been adopted to incorporate the effect of steel–concrete
interaction into the analysis of RC structures. One of the earliest attempts to describe the
bond in RC beams is that of Ngo and Scordelis, (1967). They used a link element which
can be conceptually thought of as consisting of two linear springs parallel to a set of
orthogonal axes. These axes can be in sync with the steel bar longitudinal and tangential
directions. The link element has no physical dimensions, and only its mechanical
properties are of importance, i.e., longitudinal and tangential stiffness. Therefore, link
elements can be placed anywhere in the beam without disturbing the beam geometry.
This link element represents the bond between concrete and steel and can permit a certain
amount of slip to take place during stress transfer through that bond.
Based on the bond zone mentioned above, some researchers combined 2D interface
elements with regular 2D elements representing steel and concrete (Rots, 1988; Clement,
1987). The nonlinear response of concrete near the steel bar is lumped into a fictitious
interface that has special constitutive equations. In addition some mechanisms such as the
friction that depends on the shape of the steel bars (deformed or not) have been included
into the behavior of interface elements. The implementation of such relations raises the
problem of the identification of additional material properties. Others lumped the damage
zone characteristics into an interface of zero thickness. In the work of (Dragosavic and
Groneveld, 1984), it is assumed that the thickness of the damaged layer around the
reinforcement is equal to the radius of the bar. This technique is justified using
homogenization in (Clement et. al., 1985). Thus omitting interface elements and
considering that concrete is progressively damaged around the steel bar is strictly
24
equivalent to using interface elements if the bar presents surface deformations, i.e., if
displacements are continuous at the interface. In finite element applications, omitting the
interface elements should certainly reduce the computational cost. However, it requires
the implementation of accurate constitutive equations for concrete which can be a very
sensitive subject.
Pijaudier-Cabot et. al. (1991) introduced a non-local damage mechanics approach to
analyze the bond between steel and concrete. A scalar damage variable was used to relate
the stresses in the concrete to the strains. Critical state of concrete stress at the interface
will initiate that damage process which will eventually expand to propagate damage to
the surrounding concrete. They compared their results to experiments of pull-out test.
Many researchers recently addressed the problem of steel-concrete interaction (Cox
and Herrmann, 1998, 1999; Ghandehari et. al., 1999; Soh et. al., 1999; Monti and
Spacone, 2000; Ayoub and Filippou, 2000; Salari and Spacone, 2001; Ben Romdhane
and Ulm, 2002; Limkatanyu and Spacone, 2003; Spacone and El-Tawil, 2004; Luccioni
et. al., 2005; Liang et. al., 2005; Rabczuk et. al., 2005; Kwak and Kim, 2006; Ragueneau
et. al., 2006; Wu et. al., 2006; and others). On the other hand, many researchers assumed
full bonding between the reinforcing steel and concrete in order to allow for detailed
modeling of other phenomena observed in RC (Bazant and Pijaudier-Cabot, 1987;
Feenstra, 1993; He et. al., 2006; Junior and Venturini, 2007; and others).
2.4 Methodology
It was mentioned earlier that because of the many factors affecting the nonlinear
behavior of RC and the differences in short- and long-term behavior of the constituent
materials, it is a common practice among researches to model the short- and long-term
response of RC members and structures 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 RC material. This
will be the adopted approach in this study of short-term behavior of RC beams under
monotonic loads. Time dependent effects such as creep, shrinkage, and temperature
change will not be considered here. Thermodynamically consistent derivation of the
formulations will be shown throughout this work. The RC behavior predicted by the
proposed model will be eventually compared to results obtained through experiments
conducted by other researchers. In the following, the methodologies used to describe the
behavior of concrete, reinforcing steel, and their bond-interaction analysis is presented.
2.4.1 Concrete Plasticity
The plasticity concrete model selected for this study is the model presented originally
by Lubliner et. al., (1989) (also known as the Barcelona model) and later modified by Lee
and Fenves, (1998) and Wu et. al., (2006). The model is based on an internal variable-
formulation of plasticity theory for the nonlinear analysis of concrete. The model makes
use of the fact that concrete eventually exhibits strain-softening in tension and
compression, leading to complete loss of strength. In this regard concrete resembles such
25
materials as cohesive soils, and may be classified with them as frictional materials with
cohesion, where the eventual loss of strength may be thought of as the vanishing of the
cohesion (Lubliner et. al., 1989). For such a model to be capable of representing the
behavior of concrete materials, the yield criterion must be of the form in which the
concept of cohesion is clear, and the hardening rule should be such that it will eventually
lead to vanishing of the cohesion (total damage).
Any plasticity model must include three basic components: an initial yield surface that
defines the stress level at which plastic deformation begins, a hardening rule that defines
the change of loading surface as well as the change of the hardening properties of the
material during the course of plastic flow, and a flow rule which gives an incremental
plastic stress-strain relation to allow for the plastic strain evolution during the course of
loading.
The adopted yield surface accounts for plasticity in tension and compression, and is
observed as a successful yield criterion in simulating the behavior of concrete under
uniaxial, biaxial, multiaxial, and cyclic loadings (Lee and Fenves, 1998, 2001; We et. al.,
2006). Multiple hardening variables are introduced through this yield surface to account
for different hardening. The uniaxial strength functions are factored into two parts,
corresponding to the effective stress and the degradation of elastic stiffness. The
constitutive relations for elastoplastic responses are decoupled from the degradation
damage response (discussed in the next section), which provides advantages in the
numerical implementation. This criterion is given in the effective (undamaged)
configuration as follows:
(
)
21 maxmax
ˆˆ
3()H()1()0fJI c
αβκ σσ ακ
±−
=++ −− =
−
(2.1)
where
2
/2
ij ij
Jss= is the second-invariant of the effective deviatoric stress
/3
ij ij kk ij
s
σ
σδ
=− ,
1 kk
I
σ
= is the first-invariant of the effective stress tensor
ij
σ
,
κ
±
are
the equivalent plastic strains,
max
ˆ
H( )
σ
is the Heaviside step function (H = 1 for
max
ˆ
0
σ
>
and H = 0 for
max
ˆ
0
σ
< ), and
max
ˆ
σ
is the maximum principal stress. The parameters
α
and
β
were originally defined by Lubliner et. al. (1989) as dimensionless constants. Lee
and Fenves, (1998) modified
β
to be a function of the tensile and compressive
cohesions:
00
00
(/)1
2( / ) 1
b
b
ff
ff
α
−−
−−
−
=
−
(2.2)
()
()(1 ) (1 )
()
c
c
κ
β
κα α
κ
−
±
+
=
−−+
−
+
(2.3)
where
0b
f
−
and
0
f
−
are the initial biaxial and uniaxial compressive yield stresses,
respectively,
−
c and
+
c are the tensile and compressive cohesions, respectively.
26
Figure 2.12 Concrete behavior under uniaxial loading in, a) tension, b) compression
(ABAQUS, 2004)
The tensile and compressive cohesions
±
c in Eq. (2.3), which are functions of the
tensile and compressive equivalent plastic strains
κ
±
, should be scaled so that their initial
27
values are
0
±
f
, the initial yield strength in uniaxial tension and compression, respectively,
and their final values vanishes when the damage parameters reach their maximum.
However, unlike the usual plasticity models with isotropic hardening, the cohesions
±
c
themselves are assumed to be internal variables governed by rate equations that simulate
their evolution.
It should be noted here that for tensile loading, damage and plasticity are initiated
when the equivalent applied stress reaches the uniaxial tensile strength
o
f
+
as shown in
Fig. 2.12a. However, under compressive loading, damage is initiated at a different stage
than plasticity. Once the equivalent applied stress reaches
o
f
−
(i.e. when nonlinear
behavior starts) damage is initiated, whereas plasticity occurs once
u
f
−
is reached (Fig.
2.12b). Therefore, generally
ou
f
f
+
+
= for tensile loading, but this is not true for
compressive loading (i.e.
ou
f
f
−−
≠ ).
The flow rule gives the relation between the plastic flow direction and the plastic
strain rate. In contrast with metals, the nonassociative flow rule is necessary to control the
dilatancy in modeling frictional materials (Chen and Han, 1988). Because the yield
surface in Eq. (2.1) is a combined geometric shape from two different Drucker-Prager
type functions, a Drucker-Prager type function is used as the plastic potential function
(Lee and Fenves, 1998, 2001; Wu et. al., 2006).
21
3
pp
FJI
α
=+
(2.4)
2
3
2
3
p
ij
p
ij
ij
s
F
J
α
δ
σ
∂
=+
∂
(2.5)
The flow rule also connects the loading surface/function and the stress-strain relation.
When the current yield surface
f
is reached, the material is considered to be in plastic
flow state upon increase of the loading. The flow rule is presented as follows:
p
p
ij
ij
F
ελ
σ
∂
=
∂


(2.6)
where
λ

is the plastic loading factor known as the Lagrangian multiplier.
2.4.2 Concrete Damage Mechanics
Damage mechanics can be illustrated using the effective stress concept proposed first
by Kachanov (1958) as explained below: Consider a uniform bar subjected to a uniaxial
uniform tensile stress,
σ
, as shown in Fig. 2.13a. The cross-sectional area of the bar in
the stressed configuration is
A and it is assumed that both voids and cracks appear as
damage in the bar. The uniaxial tensile force
T , acting on the bar is expressed using the
relation TA
σ
= . In order to use the principles of continuum damage mechanics, one
28
considers a fictitious undamaged configuration (effective configuration) of the bar as
shown in Fig. 2.13b. In this configuration all types of damage, including voids and
cracks, are removed from the bar. The effective stressed cross-sectional area of the bar in
this configuration is denoted by
A and the effective uniaxial stress is
ij
σ
. The bars in
both the damaged configuration and the effective undamaged configuration are subjected
to the same tensile force
T
. Therefore, considering the effective undamaged
configuration, one obtains the relation
TA
σ
= . Equating the two expressions of T
obtained from both configurations, the following expression is derived:
A
A
σ
σ
= (2.7)
Moreover, as it is seen from Fig. 2.13, the effective area
A is obtained from A by
removing the surface intersections of the micro-cracks and cavities (Voyiadjis and
Kattan, 1999, 2006) and correcting for the micro-stress concentrations in the vicinity of
discontinuities and for the interactions between these effects. Therefore, the damage
parameter
Ï•
in case of uniaxial loading can be defined as follows:
1
A
A
Ï•
=
− (2.8)
In the above equation, in the case of no damage (effective state) in the material the
damage parameter is equal to zero (i.e.
0
Ï•
=
for AA
=
). The critical damage
cr
Ï•
Ï•
=
corresponds to the rupture of the element. Lemaitre (1984) showed that the damage
parameter value ranges between 0.2 and 0.8 (
0.2 0.8
Ï•
≤
≤ ) for metals. The theoretical
value of damage parameter,
Ï•
, for the general case lies in the range 01
Ï•
≤≤.
T
A
T
φ
Remove Both
Voids and Cracks
T
A
T
σ
Effective Undamaged
Configuration
(
b
)
Damaged
Configuration
(a)
Figure 2.13 A cylindrical bar subjected to uniaxial tension: both voids and cracks
are removed (Voyiadjis and Kattan, 1999, 2006)
29
The effective area
A can be obtained through mathematical homogenization
techniques (Voyiadjis and Kattan, 1999, 2006). Homogenization techniques can be used
when the shape and the size of the defects are known which is somewhat difficult to
obtain even with electron microscopes.
Making use of Eqs. (2.7) and (2.8), one obtains the following expression for the
effective uniaxial stress
σ
(Kachanov, 1958):
1
σ
σ
Ï•
=
−
(2.9)
It should be noted that the undamaged (effective) stress
σ
can be considered as a
fictitious stress acting on an undamaged equivalent area
A .
For a three-dimensional state of stress, Eq. (2.9) can be generalized for isotropic
damage as follows:
1
ij
ij
σ
σ
Ï•
=
−
(2.10)
where
ij
σ
and
ij
σ
are the stress tensors in the damaged and effective configurations,
respectively.
The transformation equations from the nominal configuration to the effective one can
be obtained using either the strain energy equivalence hypothesis or the strain
equivalence hypothesis (Voyiadjis and Kattan, 1990, 1999, 2006). The strain equivalence
hypothesis states that the strain in the effective configuration is equal to the strain in the
nominal configuration such that in the constitutive equations one can simply replace the
nominal strain by the corresponding effective strain. On the other hand, the strain energy
equivalence hypothesis states that the elastic strain energy density in the damaged
configuration is equal to the elastic strain energy density in the effective (undamaged)
configuration.
The constitutive equation for the damaged material is written in terms of that of the
virgin material; that is, the damaged material is modeled using the constitutive laws of
the effective undamaged material in which the Cauchy stress tensor
ij
σ
is replaced by the
effective stress tensor,
ij
σ
(Murakami and Ohno, 1981):
ij ijkl kl
M
σ
σ
=
(2.11)
where
ijkl
M
is the fourth-order damage effect tensor which has many definitions in
literature (isotropic or anisotropic). One of these definitions is given by Voyiadjis and
30
Venson (1995) and Voyiadjis et. al. (2008a,b) for anisotropic damage as follows (in an
effort to symmetrize the stress tensor):
(
)
()
1
2
ijkl ij ij kl ij kl kl
M
δϕδ δδϕ
−
⎡
⎤
=−+−
⎣
⎦
(2.12)
In this study, two isotropic damage variables are used to model the tensile and
compressive damage phenomena in concrete material. The constitutive equation is
written in terms of an accumulated damage variable that is function of the tensile and
compressive damage variables and the stress state using the spectral decomposition of the
stress tensor. Detailed formulation is presented in Chapter 4.
2.4.3 Steel Reinforcement
In this study the reinforcing steel is modeled as a linear elastic, linear strain hardening
material with yield stress
y
σ
, as shown in Fig. 2.9. A
2
J elasto-plasticity model with
linear hardening will be adopted to describe the behavior or steel reinforcement. The von
Mises yield criterion, associative flow rule and isotropic hardening are suitable for
modeling structural steel. The von Mises yield criterion can be written as:
(,) () 0
eq y
FR Rp
σ
σσ
=
−− ≤ (2.13)
where
y
σ
is the yield stress,
()
R
p
is the isotropic hardening stress (linear function of the
accumulated plastic strain
R
kp= ), and
eq
σ
is the von Mises equivalent stress defined
by:
3
2
eq ij ij
SS
σ
= (2.14)
where
ij
S is the deviatoric part of the Cauchy stress tensor:
1
3
ij ij mm ij
S
σ
σδ
=− (2.15)
The plastic flow rule that governs the evolution of the plastic strain is given as follows:
p
ij
ij
f
ελ
σ
∂
=
∂


(2.16)
The scalar
λ

is the plastic multiplier which is equal in this case to the rate of the
accumulated plastic strain:
31
2
3
pp
ij ij
p
λ
εε
==



(2.17)
where
p
ij
ε

is the plastic strain rate.
2.4.4 Steel-Concrete Bond and Interaction
The importance of bond consideration in analysis of RC structures is shown through
the reinforcing bar stress-strain diagram in Fig. 2.14. Significant difference in the initial
behavior of the bar can be observed depending on the degree of bond. The initial stiffness
in both the normal and weak bond cases is smaller than for the theoretical full bond case
(Monti and Spacone, 2000). As for the post-yield phase, the normal bond case is
practically identical to the full bond case. It is obvious from the figure that the nature of
bond failure drastically changes the stress-strain behavior of the bar.
Figure 2.14 Monotonic stress-strain response of reinforcing bar fiber with different
degrees of bond (Monti and Spacone, 2000)
The effect of bond between steel and concrete will be incorporated in this study
through the description of the steel constitutive model. The global effect of bond and
interaction on the stress-strain diagram will be accounted for by reducing the yield stress
and the elastic and hardening moduli of the steel bars according to the parametric studies
conducted by Belarbi and Hsu (1994), Kwak and Kim (2006), and the references therein.
Further discussion is provided in Chapter 5.
2.4.5 Numerical Implementation of the Model
Once the material constitutive models for concrete, steel reinforcement, and bond-
interaction are developed, they will be all implemented into the advanced finite element
analysis software ABAQUS via a user defined material subroutine (UMAT). While
ABAQUS performs the standard finite element procedure using standard types of finite
elements, the UMAT will govern the behavior of these materials during different loading
32
stages, i.e., elastic, inelastic, failure, post-failure loads. Combining a UMAT subroutine
with standard finite element procedure will manifest the simplicity and applicability of
the proposed model into any engineering problem where the average user will not be
exposed to the complexities associated with introducing non standard finite elements.
33
CHAPTER 3
REINFORCING STEEL MATERIAL MODEL
3.1 Introduction
Steel reinforcement is used in RC structures to provide the tensile strength that
concrete lacks. The properties of reinforcing steel, unlike concrete, are generally well
known and do not dependent on environmental conditions or time (creep). Reinforcement
is usually classified on the basis of geometrical properties, such as size and surface
characteristics, and on the mechanical properties, such as characteristic yield stress and
ductility.
Typical stress-strain curves for reinforcing steel bars used in RC construction are
obtained from coupon tests of bars loaded monotonically in tension. For all practical
purposes, reinforcing steel exhibits the same stress-strain curve in compression as in
tension. Thus, the specification of a single stress-strain relation is sufficient to define the
material properties needed in the analysis of RC structures. The steel 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
(softens) until fracture occurs (Fig. 3.1).
Figure 3.1 Experimental Stress-Strain curve for reinforcing steel
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 (Wang and Salmon, 1998). The experimentally obtained
34
stress-strain diagram for the reinforcing steel is usually replaced in research by an
idealized characteristic diagram (Fig. 3.2).
Figure 3.2 Idealizations of the steel stress-strain curves.
Two different idealizations, shown in Fig. 3.2, are commonly used in literature
depending on the desired level of accuracy (ASCE 1982). The first idealization neglects
the strength increase due to strain hardening and the reinforcing steel is modeled as a
linear elastic, perfectly plastic material, as shown in Fig. 3.2a. 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. At several instances it is necessary to evaluate the steel
stress at strains much 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. The second idealization - the case of the bilinear
stress strain models - would be more appropriate to model the steel behavior. As shown
in Fig. 3.2b, the reinforcing steel is modeled as a linear elastic, linear-plastic-hardening
material. The parameters of these models are the stress and strain at the onset of yielding
and the stress and strain at the ultimate load (Fig. 3.3). These parameters can be derived
from experimentally obtained stress-strain relations.
The first idealization (Fig. 3.2a) is frequently used by researches (Mazars and
Pijaudier-Cabot, 1989; Sumarac et. al., 2003; Kratzig and Polling, 2004; Luccioni and
Rougier, 2005; and Wu et. al., 2006) because of its simplicity, especially in the cases
where most of the interesting activities in the RC member occur before the strain-
hardening of reinforcing steel starts (Vecchio and Collins, 1982). The second idealization
(Fig. 3.2b) is also used in many studies (Ngo and Scordelis, 1967; Feenstra, 1993; Kwak
and Fillippou, 1997; Lowes, 1999; Tikhomirov and Stein, 2001; Rabczuk et. al., 2005;
35
Phuvoravan and Sotelino, 2005; He et. al., 2006; and Junior and Venturini, 2007). It is
claimed by some researchers (e.g. Kwak and Fillippou, 1997 and the references therein)
that adopted the second idealization that the use of the elastic-perfectly plastic model
shown in Fig. 3.2a leads to numerical convergence problems near the ultimate member
strength in situations where yielding is accompanied by a sudden increase in the
deformation of the member under monotonic bending moments. In these situations, the
RC member is greatly affected by the hardening of reinforcing steel, the fact that
encourages the use of bilinear stress-strain models to improve the numerical stability of
the solution. They also claimed that 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 (Fig. 3.3). Despite the discussion, both idealizations have been reported to
be successfully implemented into FE codes used to model RC structures.
Figure 3.3 Linear elastic, linear strain hardening steel stress-strain relation
Some researchers extended the bilinear stress-strain representations to models
involving linear elastic behavior followed by nonlinear strain hardening (e.g., Marfia et.
al., 2004). A recent and more accurate idealization of the reinforcing steel behavior is the
one that accounts for the different phenomena/stages governing the steel behavior
(elasticity, perfect-plasticity, and strain-hardening-plasticity) as shown in Fig. 3.4. The
material parameters of this idealization are the stress and strain at the onset of yielding,
the strain at the onset of strain hardening, and the stress and strain at the ultimate load
(Selby and Vecchio, 1997; Fantilli et. al., 2002; and Hoehler and Stanton, 2006). Again,
these parameters can be derived from experimentally obtained stress-strain relations. This
36
idealization requires the development of a new criterion to determine the strains at which
hardening starts, a procedure that increases the complexity of this idealization technique.
Different approaches can be used to model reinforcing steel in RC structure using the
FE method. Some researchers chose to use two dimensional elements to represent steel
reinforcement in their FE meshes (Ngo and Scordelis, 1967; Mazars and Pijaudier-Cabot,
1989; Pijaudier-Cabot et. al., 1991; Luccioni and Rougier, 2005; and Wu et. al., 2006).
On the other hand, one dimensional FE representation of steel reinforcement is more
widely used to avoid introducing the complexities associated with multiaxial constitutive
relationships for structural steel (Chen, 1982; Feenstra, 1993; Kwak and Fillippou, 1997;
Faria et. al., 1998; Lowes, 1999; Soh et. al. 2003; Sumarac et. al., 2003; Kratzig and
Polling, 2004; Rabczuk et. al., 2005; Phuvoravan and Sotelino, 2005; and He et. al.,
2006).
Figure 3.4 Linear elastic, yield plateau, plastic hardening steel stress-strain relation
The introduction of the steel reinforcing elements into RC FE analysis can be done in
many different ways. The discrete element model, the distributed or smeared model, and
the embedded model are examples of representing steel reinforcement in an RC FE mesh.
The most widely used approach is the discrete element (e.g., Wu et. al., 2006), where
steel and concrete are modeled as two distinguished elements with different material
properties and behaviors. Whereas, in the distributed model, the reinforcing steel is
assumed to be distributed over the concrete elements at a certain orientation angle, and in
the embedded steel model the reinforcing steel is considered as an axial member built
into the isoparametric elements representing the concrete. A significant advantage of the
discrete representation, in addition to its simplicity, is that it eliminates the need to
develop more sophisticated types of finite elements needed by the other models (e.g.,
Soh et. al., 2003), and therefore facilitates the use of classical finite elements available to
a wider spectrum of engineers and researchers.
In this study, the steel reinforcing bars are modeled by two-dimensional constitutive
relations to comply with the primary objective of this work which is the two-dimensional
FE analysis of RC structures. The elastic-plastic with linear strain hardening model (i.e.,
37
the second idealization, Fig. 3.2b) will be eventually used in the RC beam analysis in a
subsequent chapter. Nevertheless, numerical algorithms facilitating the use of different
idealizations will be discussed and presented in this chapter.
The following sections describe how metal plasticity may be applied to solve practical
boundary-value problems. Emphasis is given to the development of numerical procedures
used in solving elastic-plastic boundary-value problems by the nonlinear FE method and
the theory of plasticity. The application of the FE method to nonlinear mechanics will be
presented first, followed by the discussion of the aspects related to the constitutive
equations. Implementation of incremental elastic-plastic constitutive relations in a
numerical process where strain and stress increments are no longer infinitesimal but with
a finite size, will also be discussed and presented.
3.2 The Application of the FE Method to Nonlinear Continuum Mechanics
The basic formulations and procedures of the FE method for elastic analysis are
described briefly as an introduction, followed by a description of the formulations and
procedures for elastic-plastic analysis, where elastic-plastic FE analysis procedure results
in a set of simultaneous nonlinear equations, the solution of which will also be discussed.
3.2.1 FE Analysis Procedure for Elastic Problems
Recalling the classical statement of a boundary-value problem in quasi-static
materials, any arbitrary solid body can be considered. The body before deformation
occupies a volume V and is subjected to forces per unit volume
i
q , traction forces acting
per unit surface area
i
T , and displacement boundary conditions imposed on certain parts
of the body
i
u . The general governing equation of the FE method for a static elastic
analysis may be derived from the principle of virtual work (weak form) as follows:
ij ij i i i i
VAV
dV T u dA q u dV
σδε δ δ
=+
∫∫∫
(3.1)
where
ij
σ
and
ij
ε
are the stress and strain tensors, respectively,
i
u
δ
and
ij
δ
ε
are virtual
displacement and virtual strain increments, respectively, where these increments form a
compatible set of deformations, and
ij
σ
with
i
T and
i
q form an equilibrium set. In matrix
form, which is adopted more frequently for FE formulations than tensorial notations, Eq.
(3.1) can be represented as:
{}{} {}{} {}{}
TTT
VAV
dV u T dA u q dV
δε σ δ δ
=+
∫∫∫
(3.2)
where the vectors for displacement
{
}
u , strain
{
}
ε
(similarly
{
}
δ
ε
and
{
}
u
δ
) and stress
{
}
σ
are defined in general (3D) terms as follows:
38
{
}
{
}
{}
{}
,,
,,
T
xyz
T
xyz
uuuu
uuuu
δδδδ
=
=
(3.3)a,b
{
}
{
}
{}
{}
,,, , ,
,,, , ,
T
xyzyzzxxy
T
xyzyzzxxy
εεεεγγγ
δ
εδεδεδεδγδγδγ
=
=
(3.4)a,b
{}
{
}
,,,,,
T
xyzyzzxxy
σ
σσστ τ τ
= (3.5)
For a small-deformation analysis, the following well-established strain displacement
relationship is obtained:
{
}
[
]
{
}
{}
[]
{}
Lu
Lu
ε
δ
εδ
=
=
(3.6)a,b
where
[
]
L
is the differential operator matrix defined as:
[]
000
000
00 0
T
x
zy
L
yzx
zyx
⎡
⎤
∂
∂∂
⎢
⎥
∂
∂∂
⎢
⎥
⎢
⎥
∂
∂∂
=
⎢
⎥
∂
∂∂
⎢
⎥
⎢
⎥
∂∂∂
⎢
⎥
∂∂∂
⎣
⎦
(3.7)
In the displacement-based FE method, the body being analyzed is approximated by an
assemblage of discrete finite elements interconnected at nodal points on their boundaries.
The nodal points are numbered from 1 to N, and the elements are numbered from 1 to M.
The displacement of the body measured in a given coordinate system is approximated as
a piece-wise continuous function over the body such that the function is continuous
inside each element and the continuity at element boundaries is assumed to be a certain
degree. The displacement at the nodal points in an FE system defines the displacement
vector
{
}
U
as follows:
{}
{
}
11 1,2 2 2
, , , , ,.......,
T
n
Uuvwuvw w= (3.8)
Within an element, say element m, the displacement is approximated as:
{
}
[
]
{
}
m
m
uNU= (3.9)
where
[
]
m
N
is a matrix of 3 by 3N, and is the matrix of displacement interpolation
function or the shape function of the element (m). For example, if the element (m) has a
39
four nodal points, say points i, j, k, l, the shape function matrix will have non-zero sub
matrices of 3 by 3 only at columns i, j, k, l, and may be expressed as:
[
]
0,...,0, , 0,...,0, , 0,...,0, , 0,...,0, ,0,..., 0
ijkl
m
N NINININI
⎡⎤
=
⎣⎦
(3.10)
where
[
]
I
is a 3 by 3 unit (Identity) matrix, and
[
]
0 is a 3 by 3 zero matrix. Substituting
Eq. (3.9) into Eqs. (3.6)a,b, one obtains the strain vectors for the element (m) as:
{
}
[
]
{
}
{}
[]
{}
m
m
m
m
B
U
B
U
ε
δ
εδ
=
=
(3.11)a,b
where
[
]
m
B
is called the strain-displacement matrix, and defined as:
[
]
[
]
[
]
mm
B
LN= (3.12)
Substituting Eqs. (3.11)a,b and (3.9) into Eq. (3.2), one obtains the governing equation
for small-deformation FE analysis as:
[
]
{}
[
]
{}
[
]
{}
mmm
TTT
mmm
mmm
VAV
B
dV N T dA N q dV
σ
=+
∑∑∑
∫∫∫
(3.13)
where m = 1,2,….,M, or summation is over every element. This equation may also be
written in a simplified form as:
[
]
{}
[
]
{}
[
]
{}
TTT
VAV
B
dV N T dA N q dV
σ
=+
∫∫∫
(3.14)
where
[
]
B
and
[
]
N take the values of
[
]
m
B
and
[
]
m
N , respectively, when the
integrations are performed in the element (m). Next, let:
{}
[
]
{}
[
]
{}
TT
AV
R
N T dA N q dV=+
∫∫
(3.15)
denote the equivalent external force acting on the nodal points, then accordingly Eq.
(3.14) can be reduced to:
[
]
{} {}
T
V
B
dV R
σ
=
∫
(3.16)
For an elastic analysis, the stress-strain relationship may be generally written as:
{
}
[
]
{
}
C
σ
ε
= (3.17)
40
where
[
]
C is the elasticity matrix. The governing Eq. (3.16) may be further written as:
[
]
{
}
{
}
[] [][][]
T
V
KU R
KBCBdV
=
=
∫
(3.18)a,b
where
[
]
K is termed the stiffness matrix of the FE system. Equation (3.18)a represents a
set of linear simultaneous equations. The displacement vector
{
}
U may be determined by
solving this set of equations, and the strain and stress of each element can then be
determined by Eqs. (3.11)a and (3.17).
It should be noted here that in Eqs. (3.9), (3.11), and (3.12), the element shape
function matrix
[
]
m
N and the strain-displacement matrix
[
]
m
B
are expressed in terms of
the global displacement vector
{
}
U . More compact form of these two matrices may be
obtained if they are expressed in terms of the element displacement vector
{
}
m
U . For a
four-node element with nodal points i, j, k, and l, one can write:
{}
{
}
,,,
T
TTTT
ijkl
m
Uuuuu= (3.19)
where each sub-vector in the equation represents the displacement vector at a nodal point.
In terms of
{
}
m
U , Eqs. (3.9) and (3.10) may be rewritten as:
{
}
[
]
{
}
mm
m
uNU= (3.20)
[
]
,,,
ijkl
m
N NINININI
⎡
⎤
=
⎣
⎦
(3.21)
In an actual FE representation/program, such compact form is usually used. The
stiffness matrix and the external equivalent force vector are computed individually for
each element, and are then assembled to the global stiffness matrix and the global force
vector, respectively.
3.2.2 FE Analysis Procedure for Elastic-Plastic Problems
For most elastic-plastic problems, closed-form solutions are not possible and
numerical solutions are sought via the FE method. Due to the nonlinear relationship
between the stress
{
}
σ
and the strain
{
}
ε
, the governing equation, Eq. (3.16), is a
nonlinear equation of strains, and thus, a nonlinear equation of nodal displacements
{
}
U .
Iterative methods are therefore necessary to solve this equation for a given set of external
forces. Moreover, because of the deformation history dependence of an elastic-plastic
constitutive relation, an incremental analysis following the actual variation of the external
41
forces must be used to trace the history of the displacements, strains and stresses along
with the applied external forces.
In an incremental analysis, the external force history may be expressed as a
progressive accumulation of external force increments in certain load steps. At the
(n+1)th step, the external force may be expressed as:
{} {} { }
1nn
R
RR
+
=
+∆ (3.22)
where the right superscript (n) is used to refer to the n-th incremental step, and
{
}
R
∆ is
the force increment from step (n) to step (n+1). Assuming that the solution at the n-th
step,
{}
n
U ,
{}
n
σ
and
{}
n
ε
are known, and at the (n+1)th step, corresponding to the load
increment
{
}
R
∆ , one obtains:
{} {} { }
{} {} { }
1
1
nn
nn
UUU
σ
σσ
+
+
=
+∆
=+∆
(3.23)a-b
{} {}
11nn
FR
+
+
= (3.24)
Using these equations, Eq. (3.16) may be rewritten as:
{}
[
]
{}
11
T
nn
V
FB dV
σ
++
=
∫
(3.25)
where
{}
1n
F
+
is the stress equivalent force acting on the nodal points. Substituting Eq.
(3.23)b into Eq. (3.22), the following can be obtained:
[
]
{} {}
[
]
{}
1
TT
nn
VV
B
dV R B dV
σσ
+
∆=−
∫∫
(3.26)
which is another way of noting that the stiffness matrix, at any (n+1) step, is a function of
the displacement. Therefore, the structural equation, Eq. (3.18)a, which now can be
written as:
(
)
{
}
{
}
KU U R=⎡⎤
⎣⎦
(3.27)
cannot be immediately solved for the nodal displacement
{
}
U because information
needed to construct the stiffness matrix
(
)
KU
⎡
⎤
⎣
⎦
is not known in advance. An iterative
process is required to obtain
{
}
U and its associated
(
)
KU
⎡
⎤
⎣
⎦
such that the product of
(){
}
KU U⎡⎤
⎣⎦
is in equilibrium with
{
}
R
.
42
Equation (3.24) in fact represents the equilibrium of the external force
{}
1n
R
+
with the
internal force
{}
1n
F
+
. Equation (3.26) is the governing equation of the incremental FE
formulations. Two types of numerical algorithms are involved in solving this equation for
the displacement increment
{
}
U∆
and the stress increment
{
}
σ
∆
. One is the algorithm
used for solving the nonlinear simultaneous equations generated from Eqs. (3.24) or
(3.26) for the displacement increment
{
}
U∆ (global/equilibrium iterations). Another is
the algorithm used to determine the stress increment
{
}
σ
∆
corresponding to a strain
increment
{
}
ε
∆
, which is computed from
{
}
U∆
, at a given stress state and a given
deformation history (local iterations) . These two algorithms constitute the nonlinear FE
procedure for elastic-plastic analysis. Many algorithms exist for solving the nonlinear
equations shown above. What follows is a description of two equation-solving techniques
applicable to the time-independent nonlinear equations.
3.3 Algorithms for Nonlinear Global/Equilibrium Iterations
The governing equation of an elastic-plastic FE incremental analysis may be generally
written in terms of the displacement
{
}
U at the incremental step (
1n
+
) as:
{}
()
{}
(
)
{
}
{}
1
1
11
n
n
nn
UFU R
+
+
++
∏= − (3.28)
This equation represents an equilibrium of the external force,
{}
1n
R
+
, with the internal
force,
{}
1n
F
+
. The iterative methods for solving this equation are therefore referred to as
(equilibrium iterative methods). They are also known as global iterations because they
are defined at the level of the entire body or structure (Doghri, 2000). This concept is
very important when dealing with the user defined material subroutine (UMAT)
associated with ABAQUS nonlinear FE analysis; this will be discussed later on. If the
weak form of equilibrium is satisfied at a time step
1n
t
+
, then the solution at
1n
t
+
has been
found, and the process can move on to the next time interval [
1n
t
+
,
2n
t
+
]. Otherwise, a new
iteration within the same time interval [
n
t ,
1n
t
+
] starts, that is, a new approximation to the
nodal displacements at
1n
t
+
is proposed by solving the nonlinear system of equations
again.
3.3.1 The Newton-Raphson Method
Assume that the (
1i − )th approximation of the displacement in the ( 1n
+
)th
incremental step has already been obtained, and is denoted by
{}
1
1
n
i
U
+
−
. Expanding
{
}
()
1n
U
+
∏ using the Taylor series expansion at
{}
1
1
n
i
U
+
−
and neglecting all higher-order
terms than the linear term, one obtains:
43
{}
()
{}
{} {}
(
)
111
11
0
nnn
ii
UUU
U
+++
−−
∂
∏
∏+ −=
∂
(3.29)
where
{}
U
∂∏
∂
is evaluated at
{}
1
1
n
i
U
+
−
, which can be written in terms of the external and
internal forces as:
{}
{}
{}
()
{}
11
1
0
nn
ii
F
FUR
U
++
−
∂
+
∆− =
∂
(3.30)
where
{}
F
U
∂
∂
is evaluated at
{}
1
1
n
i
U
+
−
,
{ } {} {}
11
1
nn
ii
UU U
+
+
−
∆= − (3.31)
and
{} {}
(
)
{
}
[]
{}
1
11 1
11 1
n
T
nn n
ii i
V
FFU B dV
σ
+
++ +
−− −
==
∫
(3.32)
Recalling that
{}
F
U
∂
∂
evaluated at
{}
1
1
n
i
U
+
−
is equal to:
[
]
[
]
[
]
1
1
nT
ep
i
V
KBCBdV
+
−
⎡⎤
=
⎣⎦
∫
(3.33)
where
ep
C
⎡⎤
⎣⎦
is evaluated at
{}
1
1
n
i
U
+
−
and is the elastic-plastic matrix, and
[
]
1
1
n
i
K
+
−
is the
tangential stiffness matrix of the structure. Therefore, the iteration scheme of a Newton-
Raphson algorithm is obtained as follows:
[
]
{}{} {}
{} {} { }
1
11
1
1
11
1
n
nn
ii
i
nn
ii i
KUR F
UU U
+
+
+
−
−
++
−
∆= −
=+∆
(3.34)a,b
and the iteration scheme starts at:
{} {}
1
0
nn
UU
+
= ,
[
]
[
]
1
0
nn
KK
+
= , and
{} {}
1
0
nn
FF
+
= . This
iteration continues until a proper convergence criterion is satisfied. Convergence criteria
will be discussed and demonstrated through the iteration procedure of a one-degree-of
freedom nonlinear system later on in this chapter (section 3.3.3). The Newton-Raphson
algorithm has a high quadratic convergence rate. However, it should be noted from Eq.
(3.34)a that the tangential stiffness matrix,
[
]
1
1
n
i
K
+
−
, is evaluated and factorized in every
44
iteration step. This can be prohibitively expensive for a large scale system. Moreover for
a perfectly-plastic or a strain-softening material, the tangential matrix may become
singular or ill-conditioned during the iteration. This may cause difficulty in finding the
solution of the nonlinear equation system. Modification of the Newton-Raphson
algorithm might therefore become necessary.
3.3.2 The Modified Newton-Raphson Method
To reduce the expensive operations of stiffness matrix evaluation and factorization,
one of the modifications of the Newton-Raphson algorithm is to replace the tangential
stiffness matrix
[
]
1
1
n
i
K
+
−
in Eq. (3.34)a by a matrix
[
]
m
K , which is the stiffness matrix
evaluated at a load step
m
, where
1mn
<
+
. If the matrix is evaluated only at the
beginning of the first lead step, the initial elastic matrix
[
]
0
K
is used throughout all load
steps. Such a method is referred to as the initial stress method. Most commonly, when the
stiffness matrix is evaluated at the beginning of each load step, or for the (
1n + )th step,
the following stiffness matrix is used:
[
]
[
]
[
]
1
0
mn n
KK K
+
== (3.35)
The iteration scheme for the Modified Newton Raphson algorithm is expressed as:
[
]
{}{} {}
{} {} { }
11
1
11
1
n
nn
ii
nn
ii i
KU R F
UU U
+
+
−
++
−
∆= −
=+∆
(3.36)a,b
and the iteration scheme starts with:
{} {}
1
0
nn
UU
+
=
, and
{} {}
1
0
nn
FF
+
=
. Similar to the
Newton-Raphson scheme, this iteration procedure continues until a proper convergence
criterion is satisfied. This modified iteration procedure is illustrated for a one-degree-of-
freedom nonlinear system in the example discussed in section (3.3.3).
The Modified Newton-Raphson algorithm involves less stiffness matrix evaluation
and factorization operations than the Newton-Raphson algorithm. This leads to a
significantly reduced computational effort in one iteration cycle for a large scale system.
However, this modified algorithm converges linearly and, in general, more slowly than
the original algorithm. Thus, for a specific nonlinear system, more iterations are needed
to reach a convergence when the Modified Newton-Raphson algorithm is used. In some
situations, such as in the analysis of a stain softening material, it may become
prohibitively slow. The convergence rate of this algorithm depends to a large extent on
the number of times the stiffness matrix is updated. The more frequently the stiffness
matrix is updated, the lesser the number of iterations needed to reach a convergence.
Moreover, the problem that the stiffness matrix may become singular or ill-conditioned
still exists.
45
Another problem associated with the Modified Newton-Raphson algorithm is that if a
change in the external load causes an unloading in the structure analyzed, this algorithm
may not result in a convergent iteration, unless the stiffness matrix is re-evaluated once
an unloading is detected. This problem increases the programming complexity of the
implementation of the Modified Newton-Raphson.
3.3.3 Examples of Using Nonlinear Equilibrium Algorithms to Solve FE Elastic-
Plastic Problems
An FE example is presented here to demonstrate the application of the first algorithm
of the two algorithms constituting the nonlinear FE procedure for elastic-plastic analysis,
i.e., using global iterations to solve for the displacement increment
{
}
U∆ . A simple
structure consisting of two horizontal bars is shown in Fig. 3.5. Bar (1) and (2) are made
of elastic-linear strain hardening materials with stress-strain relationships shown in the
figure. The two materials have the same yield stress
y
σ
and modulus of elasticity E .
t
E
is the tangential (elastic-plastic) modulus of the material.
1
/4
t
EE
=
is the tangential
modulus of material of bar (1), while
2
/2
t
EE
=
is the tangential modulus of material of
bar (2). A horizontal force 3
y
R
A
σ
= is applied at the joint connecting the two bars,
where
A
is the cross-sectional area of the bars. Two one-dimensional bar elements with
linear shape function are used for the analysis to avoid becoming consumed with the FE
formulations of more sophisticated elements, and thus concentrate on the concept of
using the nonlinear algorithms. The matrix of the shape function of such elements is
given as:
[]
11
(1 ), (1 )
22
Nrr
⎡
⎤
=− +
⎢
⎥
⎣
⎦
(3.37)
Figure 3.5 Numerical Example 1: Elastic-plastic analysis of a two bar structure
where r is the natural coordinate with origin at the center of the bar element. The strain-
displacement matrix of the bar element is then obtained as:
[] []
1
1, 1B
L
=− (3.38)
46
The element stiffness matrices of the two-bar elements are:
[] []
12
12
11 11
() ()
11 11
EuA EuA
KK
LL
−
−
⎡⎤ ⎡⎤
==
⎢⎥ ⎢⎥
−−
⎣⎦ ⎣⎦
(3.39)a,b
and the global stiffness matrix, using Eq. (3.18)b, is given as:
[] [][][]
11
112 2
22
() () 0
() () () ()
0()()
T
L
Eu Eu
A
KBCBAdx EuEuEuEu
L
Eu Eu
−
⎡
⎤
⎢
⎥
==−+−
⎢
⎥
⎢
⎥
−
⎣
⎦
∫
(3.40)
Using the governing equation, Eq. (3.25), and the boundary condition
13
0uu==, the
governing equation at (n+1)th step becomes:
() () 0fu Fu R
=
−= (3.41)
where
2
uu= is the displacement of node 2 (see Fig. 3.5), and ( )Fu is the stress
equivalent force given by:
(
)
12
() () ()
F
uA u u
σσ
=−
(3.42)
Setting up
[
]
{
}
{
}
KU R= at any given step yields the following:
11 11
12
22
2
33
12
() () 0
() ()
0()
()
(
(
)
)
A
Eu Eu u R
L
Eu Eu u R
Eu Eu
Eu Eu u R
⎡⎤⎧⎫⎧⎫
⎪
⎪⎪⎪
⎢⎥
+=
⎨
⎬⎨⎬
⎢⎥
⎪
⎪⎪⎪
⎢⎥
⎣⎦⎩⎭⎩
−
−−
⎭
−
(3.43)
After applying the loads and boundary conditions and removing the subscripts of the load
and displacement (only one load and one displacement remains), Eq. (3.43) becomes:
()()
12
() ()
A
Eu Eu u R
L
+
= (3.44)
This is a nonlinear equation of the displacement
2
uu
=
. Once this displacement is
determined, the stress and strain of the two bars may be obtained.
3.3.3.1 Solution Using Newton-Raphson Scheme
The tangential stiffness matrix of the two-bar structure may be expresses as:
47
()
12
(evaluated at ) ( ) ( )
ii i i
FA
uK EuEu
uL
∂
== +
∂
(3.45)
where
12
12
12
,,
() and ()
,,
42
yy
yy
EE
Eu Eu
EE
ε
εεε
ε
εεε
<<
⎧⎫ ⎧⎫
⎪⎪ ⎪⎪
==
⎨⎬ ⎨⎬
>>
⎪⎪ ⎪⎪
⎩⎭ ⎩⎭
(3.46)
The scheme of the iteration is as follows:
10
() , , 0
iiii
KuRFu R uu u u
−
∆= − =∆ = +∆ = (3.47)
The iteration converges in three steps (
3i
=
). These steps are shown in details below:
Step 1
Start with
0
0u = ,
()
01020
() () () 0Fu A u u
σσ
=−=,
0
()3 03
yy
R
RFu A A
σ
σ
∆
=− = −=
when
0
0u = ,
1
ε
and
2
ε
are both equal to zero, i.e.,
y
ε
ε
<
which means that
10
()Eu =
20
()Eu = E , into
0
()
i
KuRFu R∆= − =∆ gives
()
3
y
A
EE u R A
L
σ
+∆=∆=
which yields
3
2
y
L
u
E
σ
∆= . The next step starts
Step 2
10
3
2
y
L
uu u
E
σ
=+∆= , by rearranging
1
12
3
3
22
y
y
u
LE
σ
ε
εε
=− = = = which means that both
materials have yielded, and the tangential matrix must be used for the two bar elements
after yielding.
11
()
4
t
E
Eu=
,
21
()
t
Eu=
2
E
, therefore,
1
1111
( ) ( ) ( ( ) 2.375
42 22 8 4
yyyy yyyy y
EE
Fu A A A
σ
εσ ε σσσσ σ
⎛⎞⎛⎞
=+ −−− =+++=
⎜⎟⎜⎟
⎝⎠⎝⎠
1
()
i
KuRFu R∆= − =∆ or 3 2.375 0.625
42
yy y
AE E
uA A A
L
σ
σσ
⎛⎞
+∆= − =
⎜⎟
⎝⎠
which gives
0.8333
y
L
u
E
σ
∆= , the third and final step starts
Step 3
21
2.333
y
L
uu u
E
σ
=+∆= , by rearranging
1
12
2.333 2.333
y
y
u
LE
σ
ε
εε
=− = = = which
means that both materials have yielded, and the tangential matrix must be used for the
48
two bar elements after yielding.
12
()
4
t
E
Eu
=
,
22
()
t
Eu=
2
E
into
2
()Fu gives
2
4412
() ( )( ( ) 3
43 23 3 3
yyyy yyyyy
EE
Fu A A A
σ
εσ ε σσσσ σ
⎛⎞⎛⎞
=+ −−− =+++=
⎜⎟⎜⎟
⎝⎠⎝⎠
() 3 3 0
iyy
RRFu A A
σ
σ
∆=−=−=, which means that
2
()
i
KuRFu R
∆
=− =∆=0
0u∆=
end of solution procedure.
The solution procedure using Newton-Raphson scheme was shown here in full details
because it consists of only three steps and it forms the basis for all the successive
discussions related to solving the nonlinear equilibrium equations. Note that the solution
procedure terminates when
0u∆=, which originates from 0R
∆
= when the material
stiffness in not equal to zero. The term
R
∆
is defined as the out-of-balance, or the
difference between the external force
R
and the stress equivalent force F . This term,
R
∆ , can be used as a convergence criterion to end the iterative procedure.
3.3.3.2 Solution Using the Modified Newton-Raphson Scheme
Using the Modified Newton-Raphson method, the tangential stiffness matrix is
evaluated only once in the first iteration step as follows:
()
001020
(evaluated at ) ( ) ( )
FA
KuEuEu
uL
∂
==+
∂
(3.48)
The scheme of the iteration then proceeds as follows:
010
() and where 0
iii
KuRFu R u u u u
−
∆= − =∆ = +∆ = (3.49)
Only the first three iterations will be shown here to describe the procedure. The results
will be shown afterwards:
Step 1
Start with
0
0u = ,
()
01020
() () () 0Fu A u u
σσ
=−=,
0
()3 03
yy
R
RFu A A
σ
σ
∆
=− = −=
when
0
0u = ,
1
ε
and
2
ε
are both equal to zero, i.e.,
y
ε
ε
<
which means that:
10
()Eu =
20
()Eu =
E
, into
0
()
i
KuRFu R
∆
=− =∆ gives
()
3
y
A
EE u R A
L
σ
+∆=∆=
which yields
3
2
y
L
u
E
σ
∆= , where
0
2EA
K
L
= . The next step starts
Step 2
10
3
2
y
L
uu u
E
σ
=+∆= , by rearranging
1
12
3
3
22
y
y
u
LE
σ
ε
εε
=− = = = which means that both
49
materials have yielded, and the tangential matrix must be used for the two bar elements
after yield.
11
()
4
t
E
Eu= ,
21
()
t
Eu=
2
E
, therefore,
1
1111
( ) ( ) ( ( ) 2.375
42 22 8 4
yyyy yyyy y
EE
Fu A A A
σ
εσ ε σσσσ σ
⎛⎞⎛⎞
=+ −−− =+++=
⎜⎟⎜⎟
⎝⎠⎝⎠
, the
tangent stiffness matrix
0
K is not updated and remains constant
0
2EA
K
L
=
0
()
i
KuRFu R∆= − =∆ or
2
3 2.375 0.625
yy y
EA
uA A A
L
σ
σσ
∆= − = , which gives
0.3125
y
L
u
E
σ
∆= , the third step starts
Step 3
21
3
0.3125 1.8125
2
yyy
LLL
uu u
EEE
σ
σσ
=+∆= + = , by rearranging:
1
12
1.8125 1.8125
y
y
u
LE
σ
ε
εε
=− = = = which means that both materials have yielded, and
the tangential matrix must be used for the two bar elements after yielding.
12
()
4
t
E
Eu =
,
22
()
t
Eu=
2
E
into
2
()Fu gives:
2
( ) (0.8125 ) ( (0.8125 ) 2.6094
42
yyyy y
EE
Fu A A
σ
εσ ε σ
⎛⎞
=+ −−− =
⎜⎟
⎝⎠
, again the tangent
stiffness matrix
0
K is not updated and remains constant
0
2EA
K
L
=
( ) 3 2.6094 0.3906
iy y y
R
RFu A A A
σ
σσ
∆= − = − =
, which means that:
0
()
i
KuRFu R∆= − =∆, or
2
0.3906
y
EA
uA
L
σ
∆=
, which gives
0.1953
y
L
u
E
σ
∆= , the
fourth step then starts.
By observing the value of
R
∆ through the first three steps, we realized that as long as
the solution is converging, the value for
R
∆
is decreasing. This iterative procedure is
repeated until the convergence criterion
R
∆
reaches a practically small value (in this
case, when
0.0001R∆≤ ) . Twenty one iterations are performed to complete the analysis.
The results of the two iterative procedures used to solve the one-degree-of-freedom
nonlinear problem are shown in Fig. 3.6, where the stress equivalent force is normalized
in terms of
y
A
σ
, i.e., the ordinate is /
y
FA
σ
.
The results for the Modified Newton-Raphson method showing the variation of the
element normalized stresses in bar (1) and bar (2) along with the iteration steps are
demonstrated in Fig. 3.7. In the case of the Newton-Raphson method; only three points
would be plotted for each stress.
50
Newton-Raphson .vs. Modified Newton-Raphson
0
0.5
1
1.5
2
2.5
3
0 5 10 15 20
Iteration number
Stress equivalent forc
e
Modified NR
NR
Figure 3.6 Example 1: Stress equivalent force convergence iterations
Modified Newton-Raphson
-2
-1.5
-1
-0.5
0
0.5
1
1.5
0 5 10 15 20
Iteration number
Stress
Sigma bar 1
Sigma bar 2
Figure 3.7 Example 1: Element stress variation along with iterations
Figure 3.8 Numerical Example 2: Elastic-plastic analysis of a two bar structure
If the same FE example of the one-degree-of-freedom nonlinear analysis is carried out
again but with the material of bar (1) being a strain softening material,
1
/4
t
EE=− , as
shown in Fig. 3.8, while the material of bar (2) remains the same as in the first example
51
(strain hardening material with
2
/2
t
EE
=
), the following results are obtained for the
Newton-Raphson and the Modified Newton-Raphson methods. Again, only three steps
were required for the Newton-Raphson method to converge. On the other hand, the
iterative procedure converges after 70 iterations (
0.0001R
∆
≤ ) for the Modified Newton-
Raphson method (see Fig. 3.9).
Newton-Raphson .vs. Modified Newton-Raphson
0
0.5
1
1.5
2
2.5
3
0204060
Iteration number
Stress equivalent forc
e
Modified NR
NR
Figure 3.9 Example 2: Stress equivalent force convergence iterations
Figure 3.10 shows the normalized stresses in bar (1) and bar (2) along with iterations
for the Modified Newton-Raphson method. Note that as bar (1) fails to carry anymore
load, the entire load will be carried out by bar (2) in compression.
As can be seen from the previous two examples, the numerical algorithms/procedures
discussed above are applied to satisfy the weak form of the equilibrium equation. In each
step, the stresses are found and used to calculate the residual
R
∆
which in turn is used to
calculate a new
u∆ until the equilibrium equation,
0Ku
∆
=
, is satisfied. If the
equilibrium equation is not satisfied, a new step is applied where a new approximation to
the nodal displacement is proposed by solving the linear system. The procedure is
repeated until the equilibrium equation is satisfied.
Up to this point, only the global/equilibrium iterations - required to solve the nonlinear
governing equation at a specific time step for the entire structure - were discussed. The
stresses were calculated using very simple elastic E and elastic-plastic (tangent)
t
E
moduli. This is straight forward for a one-degree-of-freedom nonlinear problem, but
when discussing more complicated problems, the elastic-plastic tangential operator
depends on the incremental properties of the material as does the stresses. The
displacement increment
u∆ will be used to calculate the strain increment
ε
∆ using an
incremental format of Eq. (3.11)a, and the latter,
ε
∆
, will then be used to update the
stress. This introduces the need for the incremental theory of plasticity and the integration
of the constitutive equations to yield the required solution parameters at the local level.
52
Modified Newton-Raphson
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
0204060
Iteration number
Stress
Sigma 1
Sigma 2
Figure 3.10 Example 2: Element stress variation along with iterations
3.4 Algorithms for Nonlinear Local Iterations
What follows is a discussion of the local or Gaussian iterations, i.e., how to integrate
the plasticity constitutive model to update the quantities at the point level. The
constitutive equations and computational algorithms for rate-independent elastic-plastic
analysis of metals are presented. Everything here will be demonstrated for isotropic metal
plasticity, which will be used to model the steel reinforcement in a RC beam. Different
forms of metal strain-hardening plasticity will be presented.
3.4.1 The Incremental Constitutive Theory of Metal Plasticity
When the stress tensor,
ij
σ
, and the strain tensor,
ij
ε
, are related by a strain-dependent
matrix rather than a matrix of constants, the material studied is said to be nonlinear, i.e.,
material nonlinearity is observed. Computational challenges arise from the fact that the
equilibrium equations (discussed earlier) must be written using material properties that
depend on strains, but strains are not known in advance. The loading history and
geometry, support conditions, and material properties are assumed to be known (as
discussed in the examples of section 3.3.3). The deformations and stresses in the body as
a function of the applied load are being sought.
In the theory of plasticity - the mathematical theory of time-independent irreversible
deformations - plastic flow is the cause of material nonlinearity, the physical basis of
which involves the movement of dislocations without the influence of viscous
phenomena or presence of decohesion which damages the material (Lemaitre and
Chaboche, 1990). When the material behavior is nonlinear, material properties in a finite
element are dictated by material properties at a finite number of sampling points in each
element. Typically these points are quadrature stations of a numerical integration rule, i.e.
element centroids or Gauss points of isoparametric elements. At each point one must
keep a record of the material’s history and update the record in each computational cycle.
This is accomplished during an FE analysis using ABAQUS through storing internal
variables in the material subroutine to represent the material history parameters. The
53
number of sampling points must be small to reduce computational expense. Accordingly,
simple elements should be used. A contrary argument is that many sampling points are
needed to accurately capture the spread of yielding in individual elements. In practice, the
choice should be made between many simple elements and a small number of more
sophisticated elements.
An incremental constitutive relation for an elastic-plastic material is presented next.
Tensorial notation will be adopted in this section to comply with the tensorial nature of
the UMAT files developed in this study. For the case of isotropic elastic-plastic analysis
of metals, one should understand the general features of metal inelastic response, and
then address the key concepts in modeling this plastic behavior, which are:
•
The decomposition of strain into elastic and plastic parts.
•
The yield criterion that predicts whether the metal responds elastically or
plastically.
•
The strain hardening rule that controls the shape of the stress-strain curve in the
plastic regime.
•
The plastic flow rule that determines the relationship between stress and plastic
strain rate.
•
The plastic (loading) and elastic (unloading) condition.
The total strain tensor
ij
ε
is well accepted to be composed of two parts: a small
recoverable (reversible) elastic strain
e
ij
ε
and a large irreversible plastic strain
p
ij
ε
such
that:
ep
ij ij ij
ε
εε
=
+ (3.50)
This decomposition which is a basic assumption in the classical theory of rate-
independent plasticity is justified by the physics of the theory. An elastic deformation
corresponds to a variation in the inter-atomic distances without changes of place while
plastic deformation implies slip movements with modification of inter-atomic bonds
(Lemaitre and Chaboche, 1990). The elastic (reversible) part is related to the stress
through the linear elastic constitutive model, where the Cauchy stress tensor
ij
σ
is given
in terms of the elastic strain tensor
e
kl
ε
by:
e
ij ijkl kl
C
σ
ε
=
(3.51)
where
ijkl
C
is the fourth order Hooke’s elastic operator given as:
2
dev
ijkl ijkl ij kl
CGIK
δ
δ
=+
(3.52)
The constants
G
and K are the shear and bulk moduli, respectively given by:
54
2
and
(1 ) 3(1 2 )
EE
GK
vv
==
+−
(3.53)a,b
where
E
and v are the Young’s modulus and Poisson’s ratio, respectively, and
dev
ijkl
I
is
the deviatoric part of the fourth-rank identity tensor
ijkl
I
given as:
11
()
23
dev
ijkl ik jl il jk ij kl
I
δ
δδδ δδ
=+− (3.54)
In order to develop a description of the relation between the plastic stresses and strains,
one needs to establish a yield criterion that will predict the onset of inelastic deformation.
This is done by adopting a yield function ( , )
f
R
σ
in terms of the stress (
σ
) and the
hardening equivalent stress R (linear strain-hardening is assumed).
The yield function defines a yield surface ( , ) = 0fR
σ
(a space of admissible stresses).
Satisfying the yield surface equation indicates that the material is yielding. Drifting away
from the yield surface either indicates elastic deformation ( , ) < 0fR
σ
or progression of
plastic deformation. Progression of plastic deformations alters the strain hardening
parameters (hardening state variable
p
and its corresponding stress equivalent
R
) and
therefore modifies (expands) the yield surface enforcing ( , ) = 0fR
σ
at the updated
values of ( , )
R
σ
. This is the numerical interpretation of a physical phenomenon.
Hardening is due to an increase in the dislocation density. Higher density of dislocations
leads to more intermingle and interlock between dislocations, which causes dislocations
to block each other. Once these dislocations interlock, and even if the material is
unloaded, the updated value of the yield stress
y
R
σ
+
is the new yield stress of the
material. Thus, the situation where
(,)>0fR
σ
is not physically possible, yet, it will be
encountered during the numerical iterative procedure used to solve for material
nonlinearity. Special techniques are then required to bring the stresses to an updated yield
surface, i.e.
(,)=0fR
σ
.
By studying the general nature of polycrystalline solids, these solids are assumed to be
isotropic and their yield criteria can be assumed to be independent of the hydrostatic
pressure. Therefore, the yield criterion here will only depend on the deviatoric
components of stresses and furthermore, because of isotropy, the yield criterion will only
depend on the magnitudes (not the directions) of the deviatoric stresses (
2
J
elastoplasticity or
2
J flow theory). This leads to the von Mises yield criterion (maximum
distortion-energy criterion) where the yield criterion is a function of the invariants of the
deviatoric stresses, given here as:
(,) () 0
eq y
fR Rp
σ
σσ
=
−− ≤ (3.55)
55
where
y
σ
is the yield stress, ()
R
p is the isotropic hardening stress (linear function of the
accumulated plastic strain
p
), and
eq
σ
is the von Mises equivalent stress defined by:
3
2
eq ij ij
SS
σ
=
(3.56)
where
ij
S
is the deviatoric part of the Cauchy stress tensor, given as:
1
3
ij ij mm ij
S
σ
σδ
=− (3.57)
and
1
3
mm
σ
is the mean stress and
ij
δ
is the Kronecker delta.
Isotropic strain hardening can be modeled by relating the size of the yield surface to
plastic strain in some appropriate way. One of the most commonly used forms of
hardening stress ( )
R
p for steel reinforcement is (Ngo and Scordelis, 1967; Feenstra,
1993; Lowes, 1999; Tikhomirov and Stein, 2001; Phuvoravan and Sotelino, 2005;
Rabczuk et. al., 2005; He et. al., 2006; and Junior and Venturini, 2007):
()
R
pkp
=
(3.58)
which is known as linear isotropic hardening law,
k is the isotropic hardening constant,
and ( )
y
R
p
σ
+ represents the radius of the yield surface (a cylinder) in the space of
principal stresses. As the accumulated plastic strain ( )
p
increases, the radius of the yield
surface increases.
In order to calculate the plastic strains induced by loading beyond yield, we note that
the magnitude of the plastic strain can be determined from the hardening behavior of the
material. Since the stress must be on the yield surface at all times (the consistency
condition), and the radius of the yield surface is related to the magnitude of the
accumulated plastic strain, the magnitude of the increment of the plastic strain must be
related to the stress increment. This leads to the definition of a plastic flow rule. A plastic
flow rule derived from the yield surface (associative plastic flow rule) governs the
evolution of the plastic strain, and is given as follows:
p
ij ij
N
ε
λ
=


(3.59)
where the scalar
λ

is called the plastic multiplier and
ij
ij
f
N
σ
∂
=
∂
is the gradient of the
yield function with respect to the stress tensor. The geometric interpretation of the flow
rule is very simple. The gradient
ij
N
is perpendicular to the yield surface, and therefore,
56
the flow rule is also referred to as the normality rule. It should be noted here that the
plastic flow rule is incremental in nature.
Using the equation for the yield criterion given in Eq. (3.55), the gradient
ij
N can be
evaluated for a von Mises material as:
3
2
ij
ij
ij eq
S
f
N
σ
σ
∂
==
∂
(3.60)
and the flow rule can now be written as:
3
2
ij
p
ij
eq
S
λ
ε
σ
=


(3.61)
This is known as the Levy-Mises flow rule.
If the scalar ( )
p
, which is the accumulated plastic strain, is defined as the integration
of the rate of the accumulated plastic strain during an iterative procedure, i.e.:
()
t
o
p
tpd
Ï„
=
∫

(3.62)
where
Ï„
is an integration operator, and the rate itself is defined as:
2
3
pp
ij ij
p
ε
ε
=


(3.63)
where
p
ij
ε

is the plastic strain rate, the plastic multiplier can then be easily proven to be
equal to the rate of the accumulated plastic strain, using Eqs. (3.56), (3.61) and (3.63):
2
23 3 3
()()
32 2 2
ij ij ij ij
eq eq eq eq
SS SS
p
p
λλ
λ
σ
σσσ
λ
==
=





(3.64)
It should be noted here that the plastic multiplier
λ

is a positive scalar ( 0
λ
≥

)
determined through enforcing the consistency condition. Also known as the Kuhn-Tucker
complementary or the loading/unloading condition, the consistency condition can be
written as:
0, ( , ) 0, ( , ) 0,which gives ( , ) 0
fR fR fR
λ σ λσ λσ
≥≤ = =

 
(3.65)
57
leading to the following possible situations in the evaluation of plastic multiplier
λ

:
( , ) 0 0 (elastic)
0 0 (elastic unloading)
( , ) 0 0 0 (neutral loading)
0 0 (plastic loading)
fR
f
fR f
f
σλ
λ
σλ
λ
⎡⎤
<⇒=
⎢⎥
⎡
⎤
<⇒ =
⎢⎥
⎢
⎥
⎢⎥
=⇒ =⇒ =
⎢
⎥
⎢⎥
⎢
⎥
⎢⎥
=⇒ >
⎣
⎦
⎣⎦







(3.66)
If the updated yield function
1n
f
+
(a scalar) at the end of a step is written in terms of
its initial value at the beginning of that step
n
f
plus the rate of change df (which can be
considered equal to
f
∆ in a numerical procedure), one can write:
1
0
nn
ffdf
+
=
+= (3.67)
but 0
n
f = , which means that 0df
=
is the only possible outcome. Taking the rate of the
yield function, Eq. (3.55), results in the following:
0
ij ij
ij ij
ffffR
df p p
pRp
σσ
σσ
∂∂∂∂∂
=+=+ =
∂∂∂∂∂
  
(3.68)
and making use of Eqs. (3.50) and (3.51) along with substitutions as
1
f
R
∂
=−
∂
and
R
k
p
∂
=
∂
, we obtain the following expression for the plastic multiplier (
p
λ
=


):
ij ijkl kl
ij ijkl kl
NC
NC N k
ε
λ
=
+


(3.69)
Beyond the elastic limit, plastic flow rules are incremental in nature and stresses are
related to strains by means of an incremental constitutive relation obtained by considering
the rate of Hook’s law (Eq. (3.51)):
e
ij ijkl kl
C
σ
ε
=


(3.70)
as well as the rate of the strain additive decomposition equation, Eq. (3.50):
ep
ij ij ij
ε
εε
=
+

(3.71)
Applying Eq. (3.71) into Eq. (3.70), one obtains:
()
p
ij ijkl kl kl
C
σ
εε
=−


(3.72)
58
and substituting for the increment of the plastic strain,
p
kl
ε

, using the expression obtained
from the normality rule, Eq. (3.59), along with the expression for the plastic multiplier,
Eq. (3.69), the following equation can be obtained:
()
mn mnpq pq
ij ijkl kl kl
ab abcd cd
NC
CN
NC N k
ε
σε
=−
+



(3.73)
In order to find an operator relating
ij
σ

and
kl
ε

, a number of tensorial manipulations and
rearrangements (shown below) leads to the following incremental constitutive relation:
()
()
pq pqrs rs
ij ijmn mn mn
ab abcd cd
pq pqrs kr ls kl
ij ijmn km ln kl mn
ab abcd cd
ep
ij ijkl kl
NC
CN
NC N k
NC
CN
NC N k
C
ε
σε
δδε
σδδε
σε
=−
+
=−
+
=








(3.74)
where
ijmn mn pq pqkl
ep
ijkl ijkl
ab abcd cd
CNNC
CC
NC N k
=−
+
(3.75)
ep
ijkl
C is the fourth-order tensor known as the elastoplastic tangent operator. It has the same
symmetries as
ijkl
C , but unlike
ijkl
C ,
ep
ijkl
C
is not constant; it depends on the deviatoric
stress and the hardening parameter (Doghri, 2000).
If the expression in Eq. (3.75) is simplified to compare with the previously discussed
examples of the one-degree-of-freedom FE elastoplastic problem,
ep
ijkl
C simplifies to
t
E (section 3.3.3) which is given by:
1
t
EEk
EE
Ek Ek
⎛⎞
=− =
⎜⎟
+
+
⎝⎠
(3.76)
For the first example in section 3.3.3:
1
/4
t
EE= , but
1
1
1
t
Ek
E
Ek
=
+
which gives
1
/3kE
=
(hardening occurs)
2
/2
t
EE= but
2
2
2
t
Ek
E
Ek
=
+
which gives
2
kE
=
(hardening occurs)
For the second example in section 3.3.3:
59
1
/4
t
EE=− , but
1
1
1
t
Ek
E
Ek
=
+
which gives
1
/5kE
=
− (softening occurs)
2
/2
t
EE= but
2
2
2
t
Ek
E
Ek
=
+
which gives
2
kE
=
(hardening occurs)
where
k for a given material is a constant property known as the strain-hardening
(softening) parameter or the plastic modulus.
The elastic-plastic incremental constitutive equation, Eq. (3.74), relates an
infinitesimal stress increment with an infinitesimal strain increment at a given stress state
and plastic deformation history. However, the load increment resulting from the FE
equilibrium iterative procedure is of a finite magnitude rather than an infinitesimal one.
The resulting increments of stress and strain have finite sizes too. Therefore, the
incremental constitutive equations need be integrated numerically to compute the stress
increment. The numerical algorithm required to integrate the constitutive equations, along
with the equilibrium iterative procedure discussed earlier, are the core of the elastic-
plastic FE analysis.
It should be noted here that in ABAQUS nonlinear FE analysis associated with the
UMAT subroutine, ABAQUS carries out the global/equilibrium iterations while UMAT
integrates the material incremental constitutive model (local iterations). Yet, satisfying
the equilibrium equations (i.e., convergence) depends on the tangent operator passed back
to ABAQUS from the UMAT subroutine. This was illustrated previously in the one-
degree-of-freedom examples (3.2.2.2). Recalling that the value of tangential stiffness
matrix
i
K was calculated based on the tangential modulus
t
E , for both bar elements,
obtained in each step. In a multi-degree-of-freedom problem,
t
E will become the elastic-
plastic tangent operator
ep
C
⎡⎤
⎣⎦
calculated at each integration point and passed back to
ABAQUS where it will be used to calculate the stiffness matrix
[
]
K shown in Eq. (3.18).
3.4.2 Integrating the Incremental Constitutive Equations
The purpose here is to devise accurate and efficient algorithms for the integration of
the constitutive relations governing the material behavior. In the context of the FE
analysis, the integration of constitutive equations is carried out at the Gauss points for a
given deformation increment. The unknowns to be found are the updated stresses and
plastic variables.
Efficient algorithms should satisfy three basic requirements: consistency with the
constitutive equations to be integrated, numerical stability, and incremental plastic
consistency. The soundness and efficiency of the algorithm ensure the effectiveness of
the numerical procedure. An improper algorithm may lead not only to an inaccurate stress
solution, but may also delay the convergence of the equilibrium iterations or even lead to
divergence of the iteration.
60
In the scope of time-independent plasticity, it is common to use fictitious time
increments (pseudo-time) merely as a way of counting successive increments of loads,
stresses and strains. Therefore, we denote all quantities at time
n
t with a subscript ( n )
while all quantities at time
1n
t
+
are denoted by a subscript ( 1n
+
). Increments within the
time step [
n
t ,
1n
t
+
] (
1nn
tt t
+
−=∆) are left with no subscripts. Tensorial notation will be
dropped whenever an algorithm is introduced to reduce the complication of indices and
subscripts. Each term will be written to correspond to its original format shown in the
previous sections. The reader can refer to the original format of each equation in the
algorithms to figure out what each symbol represents (scalar, vector, or tensor). The
procedure starts at time
n
t with knowledge of the stresses
n
σ
and the total strain
increment
ε
∆ passed to the UMAT file from the previous equilibrium iteration carried
out by ABAQUS. The goal now is to find the updated stresses
1n
σ
+
, plastic strains
1
p
n
ε
+
,
and the equivalent plastic strain
1n
p
+
.
Two schemes for the integration of the constitutive equations will be presented next:
The Explicit (Forward-Euler) scheme and a Radial Return form of the Implicit
(Backward-Euler) scheme. Other integration methods exist and are numerous with
different levels of complexity. But due to the advantage of the presence of high level
processors, complicating the integration procedure will not result in any better outputs
than those obtained using these two schemes with a reasonable number of time
increments.
In both integration schemes, the first step is to use an elastic procedure to update the
stresses. If these updated stresses are found to lie within the yield surface, the material at
the Gauss point is assumed to have remained elastic. In this case, there is no need to
integrate the rate equation; elastic analysis can resume. However, if the elastic stresses
are outside the yield surface, an integration scheme is adopted to bring the stresses back
to the yield surface.
3.4.2.1 Explicit (Forward-Euler) Integration Algorithm
This integration scheme is famous for its simplicity and being straightforward to
implement. By calculating the gradient of the yield function, Eq. (3.60), and thus the
plastic multiplier, Eq. (3.69) at the beginning of the increment, the Explicit integration
scheme uses the parameters evaluated at the previous time step
n
t to calculate the
updated parameters at the end of the current time step
1n
t
+
. The following is a block
representation of the Explicit (Forward-Euler) integration scheme:
Given
n
σ
,
p
n
ε
,
n
p
and
ε
∆
If ( , ) 0
nn
fR
σ
< then exit the integration scheme-End If
If ( , ) 0
nn
fR
σ
> then
61
Calculate the plastic multiplier
n
nn
NC
NCN k
ε
λ
∆
∆=
+
Calculate the increment of plastic strain
p
n
N
ε
λ
∆=∆
Calculate the increment of elastic strain using
ε
∆
ep
ε
εε
∆
=∆ −∆
Calculate the increment of stress
e
C
σ
ε
∆
=∆
Update all the quantities
Stresses
1nn
σ
σσ
+
=+∆
Strains
1
pp p
nn
ε
εε
+
=+∆
Equivalent plastic strain
1nn
pp
λ
+
=+∆
Calculate the elastic-plastic tangent operator
ep
nn
nn
CN N C
CC
NCN k
=−
+
End If
Associated with its simplicity, however, are three of its disadvantages that need to be
discussed: a) Because it is an Explicit procedure, it is conditionally stable. b) The
accuracy of the integration depends on the increment size chosen. c) The plastic
multiplier was obtained in such a way that the yield condition is satisfied at time
n
t .
Satisfying the yield condition at time
1n
t
+
was not checked. Therefore, the solution at
time
1n
t
+
may drift away from the yield surface over many time steps. This drift can be
reduced by decreasing the size of the increment. Furthermore, if the total strain increment
1
ε
∆ (see Fig. 3.11) at time
n
t causes the stress to cross the yield surface, then
(,)0
nn
fR
σ
< , and the increment will be treated elastically, which will over predict the
location of the yield stress from
y
σ
to
1n
σ
+
. Then the new total strain increment
2
ε
∆ and
all the consequent increments will result in a shifted stress strain diagram as shown in
Fig. 3.11. This leads to accumulation of errors and overestimation of the computed
ultimate (collapse) load.
Figure 3.11 Explicit integration without correction
62
In order to minimize the drift from the yield surface, the total strain increment
1
ε
∆
causing the stress to cross the yield surface for the first time needs to be resolved into two
parts (Fig. 3.11). One part of
1
ε
∆ is elastic (
1
e
ε
∆
) taking the stress to the first yield
y
σ
while the other part is plastic (
1
p
ε
∆
). Only the plastic part needs to be passed into the
integration algorithm shown above. The elastic part needs to be considered such that:
1
e
ny
C
σ
εσ
+∆ = and
11
eee
nn
ε
εε
+
=+∆. The remainder plastic part,
111
pe
ε
εε
∆=∆−∆,
becomes the new total strain increment passed to the integration algorithm with the
known value of the stress being equal to the yield stress, i.e.,
1
p
ε
ε
∆
=∆ and
ny
σ
σ
= .
Different methods are suggested in literature for finding a scalar scaling factor
β
such
that the elastic and plastic part of the strain increment
1
ε
∆
can be found:
11
e
ε
βε
∆=∆
and
11
(1 )
p
ε
βε
∆=−∆. In terms of stresses, the yield function should be satisfied at
1
e
nn y
C
σ
βσ σ ε σ
+∆ = +∆ = , i.e.,
()0
y
f
σ
=
. Note that ( ) 0
n
f
σ
<
as shown in Fig. 3.11.
When
β
reaches its maximum value ( 1
β
=
), the elastically predicted stress would be
1nn
C
σ
σσ ε
+∆ = + ∆ , which is equal to
1n
σ
+
. At a given value of
β
, say
i
β
the stress
would be equal to
ini
σ
σβσ
=+∆ and the yield function becomes
() ( )
ini
ff
σ
σβσ
=+∆. By using a truncated Taylor series, the value of the yield function
()
i
f
σ
can be given as: () () () ( ) 0
i
in n
ii i
ff
ff f
σ
σ σ δβ σ σ δβ
σβ σ
∂
∂
∂
=+ =+∆=
∂∂
. Such a
scheme might start with an initial estimate
1
()
()()
n
o
nn
f
ff
σ
β
σ
σ
+
−
=
−
to calculate
ono
σ
σβσ
=+∆, evaluate
o
f
σ
∂
∂
and then calculate
o
δ
β
using the truncated Taylor series,
() ( ) 0
no
o
f
f
σσδβ
σ
∂
+∆ =
∂
. The next iteration then starts with
1 oo
β
βδβ
=
+ to obtain
11o
σ
σβσ
=+∆. By calculating
i
f
σ
∂
∂
, and using the truncated series again,
1
1
() ( ) 0
o
f
f
σσδβ
σ
∂
+∆ =
∂
, the next estimate of
i
δ
β
can be found. The iterations continue
until the value of
i
δ
β
practically approaches zero.
Once
β
has been evaluated, the new strain increment
1
p
ε
ε
∆
=∆ (see Fig. 3.11) will
be passed into the Forward-Euler Explicit integration algorithm shown above to complete
the analysis of this yield point defining increment. A block representation of the
correction to the Explicit (Forward-Euler) integration algorithm is shown below:
Given
n
σ
, 0
p
n
ε
= , 0
n
p = and
ε
∆
If ( ) 0
n
f
σ
< then
63
Calculate a trial elastic stress
1
trial
nn
C
σ
σε
+
=+∆
If
1
()0
trial
n
f
σ
+
<
then exit integration scheme- End If
If
1
()0
trial
n
f
σ
+
> then solve for
β
as shown above
Store the new value of
n
σ
as
ny
C
σ
βεσ
+∆=
Store the new value of
ε
∆ as
ε
βε
∆−∆
Proceed with the Forward-Euler Explicit integration
Algorithm by calculating the plastic multiplier
λ
∆
using the new value of
ε
∆
End If
End If
3.4.2.2 Implicit Integration Algorithm: Radial Return Method
This method is popular because, for the von Mises yield criterion, it takes a
particularly simple form. At a certain step in time,
n
t , applying the strain increment
ε
∆
takes the elastically updated stress
1
trial
nn
C
σ
εσ
+
+∆= outside of the yield surface (see Fig.
3.12). This trial stress is known as the elastic predictor. The stress is then updated with a
plastic corrector
p
C
ε
∆
to bring it back onto the yield surface at the end of the current
time step [
n
t ,
1n
t
+
] as follows:
(
)
1
ep p
nn n n
CC CC
σ
σεσ εεσεε
+
=
+∆ = + ∆−∆ = +∆−∆
1
trial p
n
C
σ
ε
+
=−∆. This plastic correction involves calculating the plastic multiplier
λ
∆
using the trial stress
1
trial
n
σ
+
, followed by the calculation of the increment of the plastic
strain
p
ε
∆ .
Figure 3.12 Implicit integration scheme
64
In contrast to the Explicit scheme, all quantities are written here at the end of the time
increment, ensuring that the yield condition is satisfied at the end of the time increment,
and therefore, avoiding the drift from the yield surface which was observed in the
Explicit scheme. This scheme also leads to a faster (larger time increments) solution.
In the deviatoric stress space, the plane stress von Mises ellipse becomes a circle, and
the plastic correction term is always directed towards the center of the yield surface
(because of the normality of the flow rule). This fact gives the technique its name, i.e.,
the radial return method, which can be derived as follows:
Starting with the following expression: (the elastic predictor
1
trial
n
σ
+
will be referred to
as
trial
ij
σ
):
1n trial p
ij ij ijkl kl
C
σ
σε
+
=
−∆ (3.77)
Expanding the plastic strain increment, using Eq. (3.59), one obtains:
11ntrial n
ij ij ijkl kl
CN
σσ λ
+
+
=−∆
(3.78)
Only the deviatoric part
ij
S
of the stress tensor
ij
σ
will affect the plastic analysis merely
due to the fact that the hydrostatic pressure remains constant. By obtaining the deviatoric
parts of each term in Eq. (3.78), the following expression is obtained:
11n trial n
ij ij ijkl kl
SS CN
λ
+
+
=−∆
(3.79)
However, using Eqs. (3.52) and (3.60):, the following can be derived in details for the last
term in Eq. (3.79):
()()
1
11
1
11
1
11
1
11
1
3
22
2
33
2
22
3
22
2
n
ndev ndev
kl
ijkl kl ijkl ij kl kl ijkl ij kl
n
eq
nn
ndev
kl kl
ijkl kl ijkl ij kl
nn
eq eq
n
ij
nn
ijkl kl ij
n
eq
S
C N GI K N GI K
SS
CN GI K
S
CN G GN
δδ δδ
σ
δδ
σσ
σ
+
++
+
++
+
++
+
++
+
⎛⎞
=+ =+
⎜⎟
⎜⎟
⎝⎠
=+
⎛⎞
==
⎜⎟
⎜⎟
⎝⎠
(3.80)
Note that
1n
kl
S
+
is deviatoric, therefore,
11dev n n
ijkl kl ij
I
SS
+
+
=
and
1
0
n
kk
S
+
=
(incompressible
plasticity). Substituting Eq. (3.80) back into Eq. (3.79), one obtains the following:
1
11
1
23
n
ij
ntrial ntrial
ij ij ij ij
n
eq
S
SS GNS G
λλ
σ
+
++
+
=−∆ =−∆
(3.81)
65
Multiplying
1n
ij
S
+
(Eq. (3.81)) by itself gives:
()
()
11
11
11
111
2
11
2
1
1
1
11
33
69
6
nn
ij ij
n n trial trial
ij ij ij ij
nn
eq eq
ntrial nn
ij ij ij ij
nn trialtrial
ij ij ij ij
n
n
eq
eq
n
ij i
n n trial trial
ij ij ij ij
SS
SSSGSG
SS SS
SS S S G G
SS
SS S S G
λλ
σσ
λλ
σ
σ
λ
++
++
++
+
++
++
+
+
+
++
⎛⎞⎛⎞
=−∆ −∆
⎜⎟⎜⎟
⎜⎟⎜⎟
⎝⎠⎝⎠
=−∆ +∆
=−∆
()
2
1
6
trial
j
n
eq
G
λ
σ
+
+∆
(3.82)
but using Eq. (3.79), one obtains:
111
111 11
11
111 1
33
2
nnn
ij ij ij
ntrial n n nn
ij ij ij ij ij ij
nn
eq eq
ntrial nn n
ij ij ij ij eq
SSS
SS S S G SS G
SS SS G
λλ
σσ
λσ
+
++
+++ ++
++
+++ +
⎛⎞
=+∆=+∆
⎜⎟
⎜⎟
⎝⎠
=+∆
(3.83)
Substituting Eq. (3.83) back into Eq. (3.82), one obtains:
(
)
()
()()
()
11 1
2
11
1
22
11 1
2
11 1
2
66
4126
46
nn n
ij ij eq
n n trial trial
ij ij ij ij
n
eq
n n trial trial n
ij ij ij ij eq
nn trialtrial n
ij ij ij ij eq
SS G
SS S S G G
SS S S G G G
SS S S G G
λσ
λλ
σ
λσ λ λ
λσ λ
++ +
++
+
++ +
++ +
+∆
=−∆ +∆
=−∆−∆+∆
=−∆−∆
(3.84)
Rearranging the previous equation gives the following:
()
2
11 1
46
trial trial n n n
ij ij ij ij eq
SS SS G G
λσ λ
++ +
=+∆+∆ (3.85)
Multiplying all the terms in Eq. (3.85) by
3
2
and taking the square root of both sides
gives:
()
()
()
()
()
2
11 1
2
2
11
1
33
69
22
23 3
3
trial trial n n n
ij ij ij ij eq
trial n n
eq eq eq
trial n
eq eq
SS SS G G
GG
G
λσ λ
σσ λσ λ
σσ λ
++ +
++
+
=+∆+∆
=+∆ +∆
=+∆
(3.86)
and by multiplying all terms of Eq. (3.81) by
1n
eq
σ
+
, and rearranging on obtains:
66
()
11 1 1
111
3
3
nn ntrial n
eq ij eq ij ij
nnntrial
eq ij eq ij
SSGS
GS S
σσ λ
σλ σ
++ + +
+++
=−∆
+∆ =
(3.87)
Note that the term between brackets in Eq. (3.87) is equal to
trial
eq
σ
. This is shown in Eq.
(3.86). Therefore, Eq.(3.87) can now be written as:
11trial n n trial
eq ij eq ij
SS
σσ
++
= (3.88)
and recalling that
1
,
trial n
eq eq
σ
σ
+
are scalars, multiplying both sides by
3
2
and rearranging
terms, the above relation reduces to the following:
1
1
1
33
or
22
ntrial
ij ij
ntrial
ij ij
ntrial
eq eq
SS
NN
σσ
+
+
+
==
(3.89)
which is the proof of the radial return using the von Mises yield criterion.
When the yield criterion, Eq. (3.55), is satisfied, i.e.,
11
()
nn
eq y
Rp
σσ
+
+
=+ , Eq. (3.86)
can be rewritten as:
11
1
3()3
()3 0
trial n n
eq eq y
ntrial
yeq
GRp G
Rp G
σ
σλσ λ
σλσ
++
+
=+∆=+ +∆
++∆−=
(3.90)
This equation can be used to solve for the plastic multiplier
λ
∆
. For linear hardening, the
hardening equivalent stress
1
()
n
Rp
+
, which is a function of the equivalent plastic strain
1n
p
+
, can be linearly decomposed into the following:
(
)
11
()
nn n n
R
pkpkppkpkp
++
=
=+∆=+∆ (3.91)
which facilitates rewriting Eq. (3.90) in the following format:
30
ntrial
yeq
kp k p G
σλσ
++∆+∆− = (3.92)
Rearranging the previous equation, and substituting
λ
∆
for
p
∆
(Eqs. (3.62) to (3.64)),
the following expression for the plastic multiplier
λ
∆
is derived:
33
trial n
trial
eq y
kp
f
kG kG
σσ
λ
−+
∆= =
++
(3.93)
A block representation of the Radial Return integration scheme is shown below:
67
Given
n
σ
,
p
n
ε
,
n
p
and
ε
∆
Calculate a trial elastic stress
1
trial
nn
C
σ
σε
+
=+∆
If
1
(,)0
trial
nn
fR
σ
+
< then exit integration scheme-End If
If
1
(,)0
trial
nn
fR
σ
+
> then
Calculate the plastic multiplier
3
trial
f
Gk
λ
∆=
+
Calculate the increment of plastic strain
ptrial
N
ελ
∆=∆
Calculate the increment of elastic strain using
ε
∆
ep
ε
εε
∆
=∆ −∆
Calculate the increment of stress
e
C
σ
ε
∆
=∆
Update all the quantities
Stresses
1nn
σ
σσ
+
=+∆
Strains
1
pp p
nn
ε
εε
+
=+∆
Equivalent plastic strain
1nn
pp
λ
+
=+∆
Calculate the elastic-plastic tangent operator
11
11
nn
ep
nn
CN N C
CC
NCN k
++
++
=−
+
End If
In FE nonlinear material analysis, and when the return mapping algorithm has
converged but the weak form of equilibrium, Eq. (3.1), is not satisfied, a new global
iteration is needed in order to propose new approximation to the nodal displacements
used to calculate the strain field at time
1n
t
+
during the same time increment
[
]
1
,
nn
tt
+
.
When the Newton-Raphson method, Eq. (3.29), is used to iterate on the global level, the
so called consistent tangent operator
δ
σ
δ
ε
has been reported to preserve a quadratic rate of
convergence, i.e., significantly faster convergence rate than that of the classical tangent
operator
ep
C . This results in the substitution of
ep
C calculated above by a type of
material consistent tangent operator
σ
ε
∂
∂
or one of its approximations, in order to achieve
higher convergence rates. Simo and Taylor (1985) demonstrated the importance of the
consistent linearization in preserving a quadratic rate of convergence when using the
Newton-Raphson’s method.
3.4.2.3 Other Forms of Strain Hardening: Nonlinear Strain Hardening and Linear
Strain Hardening Following Perfectly-Plastic Behavior
For nonlinear hardening, the breakdown of
1
()
n
Rp
+
, Eq. (3.91), is no longer valid, and
the above algorithms (Forward-Euler or Radial Return) need to be modified. Equations
(3.68) for Explicit and (3.90) for Radial Return, are now nonlinear functions of the
plastic multiplier
λ
∆
. These equations are solved using a local Newton-Raphson
iteration scheme. The iteration scheme is termed local to distinguish it from the
gobal/equilibrium Newton-Raphson iterations.
68
The algorithms can be easily altered to account for a nonlinear isotropic hardening
rule, such as a power hardening law, ( )
n
R
pkp= , where (
n
) is a material constant that
can be calibrated to match experimental results. The expressions given by the algorithms
for
λ
∆ , i.e. Eq. (3.69) for Explicit and Eq. (3.93) for Radial Return, are replaced by the
following expressions given as:
2
and
33
trial
ij ij
GN
f
R
R
GG
p
p
ε
λλ
∆
∆= ∆=
∂
∂
++
∂
∂
(3.94)a,b
where
1n
R
nk p
p
−
∂
=
∂
is no longer a constant. This means that the above equations will
have to be solved iteratively for
λ
∆
using a nonlinear equation solver scheme. The
scheme selected here is the Newton-Raphson iterative procedure. Once the plastic
multipliers
λ
∆
(for both algorithms) are obtained through nonlinear iterations, the rest of
the algorithms proceed as shown above for the case of linear hardening except for the
elastoplastic tangent operators
ep
C , where the linear hardening parameter k in the
expressions for
ep
C needs to be replaced with
R
p
∂
∂
.
In order to introduce an algorithm to represent the elastic, perfectly plastic followed
by strain hardening behavior of steel, the perfectly plastic behavior is discussed first
followed by the introduction of a mechanism to define the strains at the onset of
hardening. By considering the simple case of one dimensional analysis for this
demonstration, the yield criterion, Eq. (3.55), can be written as 0
y
f
σ
σ
=
−≤, and the
consistency condition (Eq. (3.68) without the hardening term) gives
p
ε
ε
=

, which
means that once plasticity starts, the strain increment is entirely plastic, which eliminates
the need to update the stress
e
E
σ
ε
=


. This explains the meaning of perfect plasticity. In
a general 2 or 3 dimensional analysis, however, this argument has to be accounted for
numerically.
Beyond the yield point, perfect plasticity dominates the analysis until the onset of
strain hardening. At that point and beyond, the strain increment is divided into elastic and
plastic parts, where the elastic part is used to update the stress. The issue here is to
determine the point at which the strain hardening starts. In a one dimensional analysis,
the onset of strain hardening is the experimentally obtained strain value. However, in a
two dimensional analysis, a more general term has to be developed to determine the onset
of hardening. In this work, the maximum strain criterion is used to define a value for the
maximum principal strain in steel, where the latter is used as an indicator of whether the
onset of hardening has been reached or not through comparing it to the experimental
value of strain at the onset of hardening. This is accomplished through adding IF
statements to the UMAT FORTRAN file after the yield condition to determine when the
stress is updated. The additional IF statements acquire the value of the maximum
69
principal strain and allow hardening and stress update to take place if the maximum strain
criterion is satisfied.
The following is a block representation of the Explicit (Forward-Euler) integration
Algorithm used to model the three phase behavior (elastic, perfectly plastic followed by
strain hardening). The same logic can be applied to the Implicit (Radial Return)
integration scheme:
Given
n
σ
,
p
n
ε
,
n
p
and
ε
∆
If ( , ) 0
nn
fR
σ
< then exit the integration scheme-End If
If ( , ) 0
nn
fR
σ
> then
Update the strain tensor
1nn
ε
εε
+
=+∆
Calculate the maximum principal strain
max
p
ε
If
max
p
ε
≥ uniaxial strain hardening then allow hardening
Else use perfectly plastic analysis
End If
Calculate the plastic multiplier
n
nn
NC
NCN k
ε
λ
∆
∆=
+
Calculate the increment of plastic strain
p
n
N
ε
λ
∆=∆
Calculate the increment of elastic strain using
ε
∆
ep
ε
εε
∆
=∆ −∆
Calculate the increment of stress
e
C
σ
ε
∆
=∆
Update all the quantities
Stresses
1nn
σ
σσ
+
=+∆
Strains
1
pp p
nn
ε
εε
+
=+∆
Equivalent plastic strain
1nn
pp
λ
+
=+∆
Calculate the elastic-plastic tangent operator
ep
nn
nn
CN N C
CC
NCN k
=−
+
End If
3.5 Implementation and Verification of the Integration Schemes
The two integration schemes discussed above were implemented into the FE
commercial code ABAQUS where their applicability, with different post yield behaviors,
was verified. The Implicit code (ABAQUS Standard) will be used. In such an Implicit
code, provision of both the integration scheme of the plasticity constitutive equations
(whether Implicit or Explicit) along with the material tangent operator/matrix is
necessary. These will be written using FORTRAN into a material behavior subroutine
called UMAT. This UMAT is then linked to ABAQUS during the execution of the
command that runs the analysis using the data stored in the input file. While ABAQUS
performs the standard FE procedure along with the nonlinear global/equilibrium
iterations, the UMAT file will govern the behavior of the materials during different
loading stages; elastic and plastic.
70
Figure 3.13 Schematic diagram of a steel plate modeled using plane stress elements
A wide range of information and parameters are passed in and out the subroutine
UMAT during FE analysis at the beginning and end of each time increment. Data passed
into UMAT at the beginning of the time increment would include, but is not limited to,
the values of stresses, strains, state variables as well as the strain increments. Data passed
back to ABAQUS would include the updated stresses, state variables and material
tangent operator (the Jacobian). The Jacobian is required for the global iterative
procedure used to minimize the force residual. If convergence occurs after a given
number of iterations, the type of the Jacobian does not influence the accuracy of the
solution as much as it influences the rate of convergence (Dunne and Petrinic, 2005).
The verification of the UMAT file is carried out under the plane stress condition
shown in Fig. 3.13. The plate is subjected to uniaxial displacement control at one of the
vertical ends while pinned at the other. Four ABAQUS CPS4R (4 nodded Continuum
Plane Stress Reduced Integration) elements are used. The analytical solution to this
problem provides guidelines that help in assessing the applicability of the UMAT
subroutines. Using the parameters given in Fig. 3.13, the strain at yield will be equal to
y
ε
= 0.00207.
At the beginning of each time increment, ABAQUS passes the stress, strain and strain
increment into UMAT. These are all passed in vector form, which can then be transferred
into tensors using proper functions/subroutines in UMAT file. The material model is then
used to update the stresses and return them back to ABAQUS in vector format. The
tangent operator is returned in the form of a matrix.
The results obtained by running the problem shown above with different UMAT files
are shown below for the case of linear hardening first, followed by other hardening laws.
The number of increments used is indicated on each figure.
Decreasing the size of the total strain increment
1
ε
∆
(Fig. 3.11) - causing the stress to
cross the yield surface for the first time - will improve the results of the Explicit
integration scheme (Fig. 3.14), eliminating the need for a correction procedure. The
correction applied to the Explicit integration will have a pronounced effect at a lesser
Plane stress
Steel plate 0.2x0.2x0.02m
3
Left edge is pinned
Right edge subjected to uniaxial
displacement control of 0.01m applied
through different time increments.
y
σ
= 414 MPa,
u
σ
= 720 MPa
E = 200 GPa ,
v = 0.3
k = 5.3 GPa
71
number of time increments (Fig. 3.15) rather than at high number of increments where
the correction would be rendered unnecessary.
Explicit Integration Scheme
0
100
200
300
400
500
600
700
800
0 0.01 0.02 0.03 0.04 0.05 0.06
Total Strain
Stress , M Pa
Explicit 100 Inc.
Explicit 1000 Inc.
Figure 3.14 Explicit integration scheme using 100 and 1000 increments
Corrected Explicit .vs. Explicit Integration Schemes
0
100
200
300
400
500
600
700
800
0 0.01 0.02 0.03 0.04 0.05 0.06
Total Strain
Stress , M Pa
Corrected Explicit 100 Inc.
Standard Explicit 100 Inc.
Figure 3.15 Explicit .vs. Corrected Explicit at 100 increments
As was discussed earlier, the Implicit integration scheme requires lesser amount of
increments to give good results than the Explicit integration scheme. The differences
between the 100 and 1000 increments curves shown below, Fig. 3.16, for the Implicit
integration scheme are the smoothness of each line, and the smoothness of the region of
transition from elastic to plastic behavior.
Fig. 3.17 shows the difference between Explicit and Implicit integrations at different
time increments. It shows the difference between the Implicit and Explicit integrations at
100 increments. It also shows that Implicit and Explicit integrations will eventually give
the same result at increased numbers of time increments.
72
Implicit Integration Scheme
0
100
200
300
400
500
600
700
800
0 0.01 0.02 0.03 0.04 0.05 0.06
Total Strain
Stress , M Pa
Implicit 100 Inc.
Implicit 1000 Inc.
Figure 3.16 Explicit integration scheme using 100 and 1000 increments
Explicit .vs. Implicit Integration Schemes
0
100
200
300
400
500
600
700
800
0 0.010.020.030.040.050.06
Total Strain
Stress , M Pa
Explicit 100 Inc.
Implicit 100 Inc.
Implicit 1000 Inc.
Explicit 1000 Inc.
Figure 3.17 Explicit .vs. Implicit integration schemes
If the results of the algorithms of linear strain hardening are compared to a
mathematical representation of the experimental results of Grade 60 steel (
y
σ
= 414 MPa,
u
σ
= 720 MPa), Fig. 3.18, it can be seen that the linear hardening parameter
k
was
chosen such that the ultimate stress of the steel is reached.
If the linear hardening rule is replaced by a nonlinear hardening rule, such as a power
hardening law, ( )
n
R
pkp= , the following figures (Fig. 3.19 and Fig. 3.20) can be
obtained. One should note the difference in the values of the yield stress between the
Explicit and Implicit integrations at 100 increments. Reducing the size of these
increments results in better agreement between the two integration schemes. The Explicit
integration results can be further enhanced by incorporating the correction procedure
discussed above.
73
Numerical .vs. Theoretical Stress-Strain Curves
0
100
200
300
400
500
600
700
800
0 0.05 0.1 0.15
Strain
Stress, M Pa
Theoretical
Implicit 1000 inc.
Explicit 1000 inc.
Implicit 100 inc.
Explicit 100 inc.
Figure 3.18 Numerical .vs. theoretical stress strain curve for Reinforcing steel
Explicit Integration 100 inc.
0
500
1000
1500
0 0.01 0.02 0.03 0.04 0.05 0.06
Total Strain
Stress , M Pa
n = 0.9
n = 0.8
n = 0.6
n = 0.5
Figure 3.19 Explicit integration of nonlinear hardening plasticity
It is worthy to note that as the power (n) increases from 0.5 to 0.9, the solution
approaches that of a linear hardening law. Substituting unity for the value of (n) is not
admissible. In addition, increasing the number of increments to 1000 will give identical
results for both integration schemes, as shown in Fig. 3.21.
A UMAT file was also used to produce an elastic perfectly plastic behavior where the
plastic slip in the reinforcing steel equals the applied strain rate which renders the
increment of the stress to be equal to zero, resulting with a plateau beyond yielding, see
Fig. 3.22. As discussed earlier, this algorithm neglects the strain hardening effect of the
steel and therefore, might not be suitable for large deformation analyses such as seismic
analysis, etc.
74
Implicit Integration 100 inc.
0
200
400
600
800
1000
1200
1400
1600
0 0.01 0.02 0.03 0.04 0.05 0.06
Total Strain
Stress , M Pa
n = 0.9
n = 0.8
n = 0.6
n = 0.5
Figure 3.20 Implicit integration of nonlinear hardening plasticity
Compare Explicit and Implicit Integrations 1000 inc.
0
200
400
600
800
1000
1200
1400
1600
0 0.01 0.02 0.03 0.04 0.05 0.06
Total Strain
Stress , M Pa
Ex n = 0.9
Ex n = 0.8
Ex n = 0.6
Ex n = 0.5
Im n = 0.9
Im n = 0.8
Im n = 0.6
Im n = 0.5
Figure 3.21 Comparison of integration schemes at 1000 increments
Elastic Perfectly Plastic Stress Strain
0
100
200
300
400
500
0 0.01 0.02 0.03 0.04 0.05 0.06
Total Strain
Stess, MPa
Corrected Explicit
100 increments
Figure 3.22 Corrected Explicit elastic-perfectly plastic analysis
75
Elastic - Perfectly Plastic - Linear Hardening
0
100
200
300
400
500
600
700
0 0.01 0.02 0.03 0.04 0.05 0.06
Total Strain
Stress, MPa
Corrected Explicit Full Range 100 inc.
Implicit 100 Full Range 100 inc.
Figure 3.23 Elastic-perfectly plastic-strain hardening analysis
Another modification was applied to the UMAT file in order to obtain a refined
description of the stress-strain behavior of reinforcing steel (Fig. 3.23). Here the analysis
beyond the yield point presumes under the perfectly plastic condition, until the point
where the strain reaches the strain hardening value of the reinforcing steel, where
hardening behavior starts. The analysis changes from a perfectly plastic to a plastic
hardening analysis leading to the stress-strain curve shown in Fig. 3.23. It should be
mentioned here that this algorithm is sensitive to increment size and sometimes it leads to
divergence of the solution. Further work is needed in order to enhance the performance of
such an algorithm and to substitute the linear hardening stage with nonlinear hardening to
better describe the steel behavior.
76
CHAPTER 4
CONCRETE MATERIAL MODEL
4.1 Introduction
The analysis and design of a concrete structure requires prior knowledge of its
mechanical properties. When continuum mechanics is considered, elastic damage models
or elastic plastic constitutive laws are generally the standard approaches to describe the
behavior of concrete. In the first case, the mechanical effect of the progressive
microcracking and strain softening are represented by a set of internal state variables
which act on the elastic behavior (decrease of the stiffness) at the macroscopic level (e.g.
Mazars, 1984; Simu and Ju, 1987a,b; Mazars and Pijaudier-Cabot, 1989; Labadi and
Hannachi, 2005; Tao and Phillips, 2005; Junior and Venturini, 2007; Khan et. al., 2007).
In plasticity models, softening is directly included in the expression of a plastic yield
surface by means of a hardening–softening function (e.g. Feenstra and de Borst, 1996;
Bicanic and Pearce, 1996; Grassl et. al., 2002; Park and Kim, 2005).
In concrete material analysis, it is very important to capture the variations
(degradation) of the elastic stiffness of the material upon mechanical loading, which
cannot be captured by plasticity-based approaches (Feenstra and de Borst, 1996).
Continuum damage mechanics is the appropriate theoretical framework for that in order
to capture material degradation. However, continuum damage models cannot capture
alone the irreversible (plastic) deformations that the material undergoes during loading.
Therefore, the combined use of elastic-plastic constitutive equations along with
continuum damage mechanics is vital to better describe the mechanical behavior of
concrete.
There are several possibilities for coupling plasticity and damage effects in a single
constitutive relation. Historically, damage has first been coupled to plasticity (Lemaitre
and Chaboche, 1984) in the so-called ductile failure approaches for metal alloys. The
underlying assumption was that void nucleation is triggered by plastic strains.
Applications to concrete were proposed among others, (e.g. Oller et. al., 1990; Voyiadjis
and Abu-Lebdeh, 1993, 1994; Abu-Lebdeh and Voyiadjis, 1993; Kratzig and Polling,
2004). In these models, damage growth is a function of the plastic strains. There is a
difficulty, however. In uniaxial tension there is little plasticity and quite a lot of damage
while in uniaxial compression, the picture is reversed with little damage and important
plastic strains. Furthermore, it can be hardly explained how plastic strain may develop in
concrete prior to microcracking. A common assumption is that irreversible strains are due
to microcrack sliding and internal friction. Such a process requires the prior formation of
internal surfaces (microcracks).
The second approach, that is more suited to both tension and compression responses,
uses the effective stress. The plastic yield function is written in the effective
configuration pertaining to the stresses in the undamaged material. Many authors (Simo
77
and Ju, 1987a,b; Ju, 1989; Mazars and Pijaudier-Cabot, 1989; Yazdani and Schreyer,
1990; Hansen and Schreyer, 1992; Lee and Fenves, 1998; Faria et. al., 1998; Fichant et.
al., 1999; Voyiadjis and Kattan, 1999, 2006; Jefferson, 2003; Salari et. al., 2004; Shen et.
al., 2004; Jason et. al., 2006, Voyiadjis et. al., 2008b) applied this approach to isotropic
and anisotropic damage coupled to elasto-plasticity. It has been extended to other sources
of damage, for instance to thermal damage by Nechnech et. al. (2002) and Willam et. al.
(2002).
A last possibility is what could be called the strong coupled approach. As opposed to
the above where the plastic yield function is written in term of the effective stress, the
applied stress appears in the plastic process, which becomes coupled to damage. Luccioni
et. al. (1996), Armero and Oller (2000) and Voyiadjis et. al. (2008a) provided the
thermodynamic consistent backgrounds of such a model.
In this chapter, an elastic plastic damage formulation is proposed to model the
nonlinear behavior of concrete materials. The model is intended to circumvent the
disadvantages of pure plastic and pure damage approaches applied separately. It is based
on an isotropic damage model, with tensile and compressive damage criteria, combined
with a plasticity yield criterion with multiple hardening rules. The isotropic damage is
responsible for the softening response and the decrease in the elastic stiffness, while
hardening plasticity accounts for the development of irreversible strains and volumetric
compressive behavior within the effective configuration.
The effective stress approach has been chosen because it provides a simple way to
separate the damage and plastic processes. Plastic effects, driven by the effective stresses,
can be described independently from damage ones and vice versa. One of the main
interests is to ease the numerical implementation which is Implicit/Explicit. The plastic
part is Implicit and the damage part is Explicit, same as in classical continuum damage
computations. As a consequence, existing robust algorithms for integrating the
constitutive relations can be implemented. The calibration of the material parameters is
also easier to handle as a consequence of the separation of damage and plasticity
processes.
In this contribution, the damage process is (elastic) strain controlled. The isotropic
damage model proposed by Tao and Phillips (2005) will be adopted here to describe the
damage behavior of concrete. While the Tao and Phillips (2005) model incorporated
strain-softening in an elastic-damage framework, it is used in this work simultaneously
with the effective stress space plasticity in order to describe the damage and irreversible
phenomena in concrete materials under tension and compression. The plastic process
shall be described using a yield function inspired from Lubliner et. al. (1989) and later
modified by Lee and Fenves (1998) and Wu et. al. (2006). It is an isotropic hardening
process in the present model. Softening is controlled by damage, while plasticity controls
hardening, in tension and compression. Fracture energy related coefficients have been
defined and incorporated in order to achieve a reasonable degree of discretization
insensitivity in numerical calculations (Feenstra and de Borst, 1996; Lee and Fenves,
1998; Wu et. al., 2006).
78
In the following sections, the general framework of the elastic plastic damage model is
discussed first, followed by a consistent thermodynamic formulation. The Helmholtz free
energy function is then introduced with specific forms to be used in the FE code. The
effective stress space plasticity is then introduced followed by the discussion of the
damage yield criteria and evolution laws for concrete under tension and compression.
Elements of validation and comparisons between the model and existing damage and
damage plasticity approaches are also presented. The behavior of the model is then tested
using four elementary loading cases: simple tension, simple compression, biaxial tension,
and biaxial compression. The model is also used to reproduce the load capacity of a
notched beam subjected to three point bending test.
4.2 Framework for the Elastic-Plastic-Damage Model
The model presented in this work is thermodynamically consistent and comes from a
generalization of the effective space plasticity theory and isotropic damage theory applied
simultaneously under the assumptions of small strains and isothermal conditions.
The transformation from the effective (undamaged) configuration to the damaged one
can be done by utilizing the strain equivalence hypothesis, which basically states that the
strains in the undamaged (effective) configuration are equal to the strains in the damaged
configuration. This hypothesis is commonly applied to the coupling of plasticity and
continuum damage mechancis (Lemaitre and Chaboche, 1990; Steinmann et. al., 1994;
Lammer and Tsakmakis, 2000; Menzel et. al., 2005; Voyiadjis and Kattan, 2006). Using
the additive decomposition of the total strain tensor
ij
ε
(=
ij
ε
) into elastic (reversible)
e
ij
ε
(=
e
ij
ε
) and plastic (irreversible)
p
ij
ε
(=
p
ij
ε
) strain tensors, along with the strain equivalence
hypothesis, the following arrangements can be assumed:
ep
ij ij ij
ij ij
ep
ij ij ij
εεε
ε
ε
εεε
⎧⎫
=+
⎪⎪
⇒=
⎨⎬
=+
⎪⎪
⎩⎭
(4.1)
The equivalence of the elastic strains will be used to obtain an expression for the
elasticity tensor ( )
ijkl
E Φ in the damaged configuration as well as the damage
thermodynamic conjugate forces Y
±
used in the damage yield criteria g
±
. The
equivalence of the plastic strains, on the other hand, will be justified through the use of
the effective stress space plasticity (Ju, 1989) and the definition of the plastic Helmholtz
free energy function. Both equivalences will be discussed later in this chapter.
By taking the time derivative of the arrangements in Eq. (4.1), the following strain rate
equations necessary for the plastic-damage incremental procedure are obtained:
ep
ij ij ij
ij ij
ep
ij ij ij
εεε
ε
ε
εεε
⎧⎫
=+
⎪⎪
⇒=
⎨⎬
=+
⎪⎪
⎩⎭




(4.2)
79
The effective stress tensor (stresses in the undamaged configuration) can now be
written in terms of the strain equivalence hypothesis and using Hook’s law as:
()
ep
ij ijkl kl ijkl kl kl
EE
σ
εεε
== − (4.3)
where
ijkl
E
is the fourth-order undamaged isotropic elasticity tensor, also known as the
undamaged elastic operator, given as:
2
dev
ijkl ijkl ij kl
EGIK
δ
δ
=+ (4.4)
where
1
3
dev
ijkl ijkl ij kl
II
δ
δ
=− is the deviatoric part of the fourth-order identity tensor
1
2
()
ijkl ik jl il jk
I
δ
δδδ
=+
, and G and K are the effective shear and bulk moduli,
respectively. The tensor
ij
δ
is the Kronecker delta, and is equal to one,
1
ij
δ
=
when ij
=
or zero, 0
ij
δ
= when ij≠ .
Taking the time derivative of Eq. (4.3), the rate of the constitutive equation in the
effective configuration can be obtained as:
()
ep
ij ijkl kl ijkl kl kl
EE
σ
εεε
== −


(4.5)
The damage configuration counterpart of Eq. (4.3), i.e. the stress tensor for the
damaged material, can be written as follows:
() ()( )
ep
ij ijkl kl ijkl kl kl
EE
σ
εεε
=Φ=Φ−
(4.6)
where
ijkl
E
is the fourth-order elasticity tensor in the damaged configuration.
The stress-strain behavior is affected by the development of micro and macro cracks in
the material body. Concrete contains a large number of micro cracks, especially at
interfaces between coarse aggregates and mortar, even before the application of external
loads. These initial microcracks are caused by segregation, shrinkage, or thermal
expansion in the cement paste. Under applied loading, further micro-cracking may occur
at the aggregate-cement paste interface, which is the weakest link in the composite
system. These microcracks which are initially small (invisible), will eventually lead to
visible cracks that extend as the applied external loads are increased. These cracks
contribute to the generally obtained nonlinear stress-strain behavior. Since a
phenomenological continuum approach is followed in this work, these effects are
smeared out (i.e. averaged) throughout the body where the material is considered as a
mechanical continuum with degraded (damaged) properties.
In this work, the stress-strain relation involves a scalar (isotropic) damage variable
Φ
which is a weighted average of the tensile and compressive damage scalar variables,
Ï•
+
80
and
Ï•
−
. It was shown in Chapter 2 that the Cauchy stress tensor
ij
σ
is related to the
effective stress tensor
ij
σ
by:
(1 )
ij ij
σ
σ
=
−Φ (4.7)
where Φ is a dimensionless scalar (i.e. isotropic) damage variable interpreted here as
averaged the crack density. It is clear from Eq. (4.7) that when the material is in the
virgin state (undamaged),
0Φ= , the effective stress
ij
σ
is equivalent to the Cauchy
stress,
ij
σ
. In the case of the damaged material, the effective stress is more representative
than the Cauchy stress because it acts on the effective area that is resisting the external
loads. Furthermore, the scalar damage variable
Φ
is still used in order to represent the
macroscopic effect of the material microdamage mechanism discussed above.
By substituting Eqs. (4.3) and (4.6) into Eq. (4.7) one obtains the following relations:
(1 )
ijkl ijkl
EE=−Φ (4.8)
(1 ) ( )
p
ij ijkl kl kl
E
σ
εε
=−Φ −
(4.9)
The expression for the fourth order elasticity tensor in the damaged configuration
ijkl
E
in
terms of its effective counterpart
ijkl
E
will also be derived from the elastic dissipation
potential later on. The damage variable
Φ
has values from zero to one. The value
0
Φ
=
corresponds to the undamaged (effective) material and the value 1
Φ
= corresponds to the
fully damaged material. Damage associated with the failure mechanisms of the concrete
(cracking and crushing) results in a reduction in the elastic stiffness (Eq. (4.8)). Within
the context of the scalar-damage theory, the stiffness degradation is isotropic (i.e. the
same damage evolution is assumed in different directions) and represented by a single
degradation value Φ . The time derivative can now be applied to Eq. (4.9) to obtain the
following:
(1 )
ee
ij ijkl kl ijkl kl
EE
σ
εε
=−Φ −Φ



(4.10)
which represents the constitutive relation for the elastic-plastic-damage model used in
this work.
4.3 Consistent Thermodynamic Formulation
In this section, a general thermodynamic framework of the elastic-plastic-damage
formulation for concrete is developed. Isothermal conditions and rate independence are
assumed throughout this work. Irreversible thermodynamic following the internal
variable procedure of Coleman and Gurtin (1967) will be applied. The internal variables
and potentials used to describe the thermodynamic processes are introduced. The
Lagrange minimization approach (calculus of functions of several variables) is used later
81
on to derive the evolution equations for the proposed model. The constitutive equations
are derived from the second law of thermodynamics, the expression of Helmholtz free
energy, the additive decomposition of the total strain rate into elastic and plastic
components, the Clausius-Duhem inequality, and the maximum dissipation principle.
The Helmholtz free energy can be expressed in terms of a suitable set of internal state
variables that characterizes the elastic, plastic, and damage behaviors of concrete. In this
work, the following internal state variables are assumed to satisfactorily characterize the
behavior of concrete both in tension and compression: the elastic strain tensor
e
ij
ε
, a set
of plastic hardening variables (
κ
+
,
κ
−
) defined as the equivalent plastic strains in tension
and compression, respectively, and the scalar damage variables (
Ï•
+
and
Ï•
−
) representing
the damage densities in the material under tension or compression, respectively, such
that:
(, , , , )
e
ij
ψψ
εκκ
Ï•Ï•
+
−+−
=
(4.11)
The constitutive model proposed here is based on the hypothesis of uncoupled
elasticity (e.g. Lubliner, 1990; Luccioni et. al., 1996; Faria et. al., 1998; Nechnech et. al.,
2002; Salari et. al., 2004; Kratzig and Polling, 2004; Luccioni and Rougier, 2005; Shao
et. al., 2006; Wu et. al., 2006). According to this hypothesis, the total free energy density
per unit volume
ψ
can be assumed to be formed by two independent parts: an elastic part
e
ψ
and a plastic part
p
ψ
, corresponding to the elastic and plastic processes respectively
(both dissipative). Therefore, the Helmholtz free energy is given as:
(, , ) ( , )
ee p
ij
ψ
ψεϕϕ ψκκ
+
−+−
=+ (4.12)
It can be noted from the above decomposition that damage affects only the elastic
properties and not the plastic ones. This can be justified by the following: once micro-
cracks are initiated during loading of a concrete material, local stresses are redistributed
to undamaged material micro-bonds over the effective area (
A shown in Fig. 2.13,
Chapter 2). Thus, effective stresses of undamaged material points are higher than nominal
stresses. Accordingly, it appears reasonable to state that the plastic flow occurs only in
the undamaged material micro-bounds by means of effective quantities (Ju, 1989). The
plastic response is therefore characterized in the effective stress space and the yield
function is no longer written in term of the applied stress, rather, it is a function of the
effective stress, i.e., the stress in the undamaged material in between the microcracks.
This approach, which is more suited for brittle materials like concrete, has been
extensively used by researchers (Simo and Ju, 1987a,b; Ju, 1989; Mazars and Pijaudier-
Cabot, 1989; Yazdani and Schreyer, 1990; Hansen and Schreyer, 1992; Lee and Fenves,
1998; Faria et. al., 1998; Fichant et. al., 1999; Salari et. al., 2004; Jefferson, 2003; Jason
et. al., 2006; and others).
In the following, the thermodynamic conjugate forces associated with the internal
state variables in Eq. (4.12) are derived based on the second law of thermodynamics. For
82
isothermal behavior, the second-law of thermodynamics states that the rate of change in
the internal energy is less than or equal to the external expenditure of power such that:
dv
ext
v
P
ρψ
≤
∫

(4.13)
where
ext
P is the external power which according to the principle of virtual power should
be equal to the internal power such that:
int
dv
ext ij ij
v
PP
σε
==
∫

(4.14)
Substituting Eq. (4.14) into Eq. (4.13), one obtains the following:
dv dv 0 ( ) dv 0
ij ij ij ij
vv v
ρψ σ ε ρψ σ ε
−≤⇔ −≤
∫∫ ∫
 
(4.15)
In a stepwise sense, the Clausius-Duhem inequality can be inferred from Eq. (4.15) as
follows:
0
ij ij
σ
ερψ
−
≥

(4.16)
Taking the time derivative of Eq. (4.12), the following expression can be written:
+
κκ
ee e p p
ep e
ij
e
ij
ψψ ψ ψ ψ
ψ
ψψ ε ϕ ϕ κ κ
εϕϕ
+
−+−
+− −
∂∂ ∂ ∂ ∂
=+= + + + +
∂∂ ∂ ∂ ∂
      
(4.17)
By substituting the rate of the Helmholtz free energy density, Eq. (4.17), into the
Clausius-Duhem inequality, Eq. (4.16), one can write the following relation:
0
eeepp
pe
ij ij ij ij
e
ij
ψψψψψ
σε σ ρ ε ρ ϕ ρ ϕ ρ κ ρ κ
εϕϕκκ
+− + −
+− + −
⎛⎞
∂∂∂∂∂
+− − − − − ≥
⎜⎟
⎜⎟
∂∂∂∂∂
⎝⎠

(4.18)
The above equation is valid for any admissible internal state variable such that the
Cauchy stress tensor can be define as:
e
ij
e
ij
ψ
σρ
ε
∂
=
∂
(4.19)
as well as the non-negativeness of intrinsic dissipation:
0
p
ij ij
YYcc
σε ϕ ϕ κ κ
++ −− ++ −−
++−−≥

(4.20)
83
where the damage and plasticity conjugate forces that appear in the above expression are
defined as follows:
e
Y
ψ
ρ
Ï•
+
∂
=−
∂
(4.21)
e
Y
ψ
ρ
Ï•
−
∂
=−
∂
(4.22)
p
c
ψ
ρ
κ
+
+
∂
=
∂
(4.23)
p
c
ψ
ρ
κ
−
−
∂
=
∂
(4.24)
The mechanical dissipation must satisfy the first (Clausius-Duhem) inequality of
thermodynamics and can be decomposed in two parts: one part due to the plastic process
p
Π and the other due to the damage process
d
Π
. The mechanical dissipation energy
function
Π can therefore be written as follows:
0
pd
Π
=Π +Π ≥
(4.25)
The plasticity and damage dissipation potentials are given, respectively, as follows:
0
pp
ij ij
cc
σε κ κ
++ −−
Π= − − ≥

(4.26)
0
d
YY
Ï•Ï•
++ −−
Π
=+≥

(4.27)
The rate of the internal variables associated with plastic and damage deformations are
obtained by utilizing the calculus of functions of several variables with the plasticity and
damage Lagrange multipliers
p
λ

and
d
λ
±

, respectively. Thus the following general
objective function can be defined:
0
p
dd
Fg g
λλ λ
++ −−
Ω=Π− − − ≥
 
(4.28)
where
F and g
±
are the plastic potential function and the tensile and compressive
damage potential functions, respectively, to be defined later.
Use is now made of the well known maximum dissipation principle (Simo and
Honein, 1990; Simo and Hughes, 1998), which describes the actual state of the
thermodynamic forces (
ij
σ
, Y
±
, c
±
) as the state that maximizes the dissipation function
over all other possible admissible states. Hence, one can maximize the objective function
Ω by using the following necessary conditions:
84
0, 0, 0
ij
Yc
σ
±±
∂Ω ∂Ω ∂Ω
=
==
∂∂∂
(4.29)
Substituting Eq. (4.28) into Eqs. (4.29) along with Eqs. (4.26) and (4.27) yields the
following corresponding thermodynamic laws:
p
pp
ij
ij
F
ελ
σ
∂
=
∂


(4.30)
d
g
Y
ϕλ
+
++
+
∂
=
∂


(4.31)
d
g
Y
ϕλ
−
−−
−
∂
=
∂


(4.32)
p
F
c
κλ
+
∂
=
∂


(4.33)
p
F
c
κλ
−
−
∂
=
∂


(4.34)
Note that Eq. (4.30) is defined in terms of a plastic potential
p
F different from F to
indicate the use of a non-associative flow rule. It is also worthy to note that in this work,
the damage criteria are characterized with scalar quantities – scalar thermodynamic
conjugate forces and scalar damage parameter - therefore, the above general
thermodynamic evolution laws for damage, Eqs. (4.31) and (4.32), will be greatly
simplified in the implementation procedure.
4.4 The Helmholtz Free Energy Function
Based on the additive decomposition of the Helmholtz free energy function into
elastic-damage and plastic parts discussed earlier, Eq. (4.12), this section introduces
specific forms for the elastic-damage and plastic parts of the Helmholtz free energy
function adopted in this work. The elastic-damage part of the Helmholtz free energy
function will be defined first, followed by a definition for the plastic part.
4.4.1 The Elastic/Damage Part of the Helmholtz Free Energy Function
In order to define the elastic-damage part of the Helmholtz free energy, the
spectral decomposition procedure for obtaining the tensile and compressive parts of the
Cauchy stress tensor, as well as the combined scalar damage variable,
Φ , are to be
defined first.
To account for the different effects of damage mechanisms on the nonlinear
performance of concrete under tension and compression, spectral decomposition of the
effective stress tensor
ij
σ
into positive and negative components (
ij
σ
+
,
ij
σ
−
) is performed
85
(e.g. Ortiz, 1985; Ju, 1989; Lubliner et. al., 1989; Faria et. al., 1998; Lee and Fenves,
1998; Wu et. al., 2006) such that:
ij ij ij
σ
σσ
+
−
=
+ (4.35)
The total effective stress tensor
ij
σ
can be written in terms of its principal values
()
ˆ
k
σ
and their corresponding principal directions
()k
i
n ( k =1, 2, 3) as follows:
3
() () () (1) (1) (1) (2) (2) (2) (3) (3) (3)
1
ˆˆˆˆ
kkk
ij i j i j i j i j
k
nn nn nn nn
σσ σ σ