# Elasto-plastic And Damage Modeling Of Reinforced Concrete

- 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

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.

- 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

ÏƒÏƒ Ïƒ Ïƒ