..A.%...,‘

.1.

c ‘11-’25;

.. .4:

5.5:,

.2”,

1.}.
a

fi
35¢sz

‘.-l

{'54le
\ : ’
1..

.3 I I
3.353.;

. U (
._ Qur. 1 1:1. .zri i. w .r: :3. 3w. V . 4%. Us

_
adj. . ‘1

 

TH'Lz‘bIS

.9321
[9037 8455

This is to certify that the
thesis entitled

THEORETICAL DERIVATION AND NUMERICAL
IMPLEMENTATION OF CONTINUUM DAMAGE BASED
CONSTITUTIVE EQUATIONS

presented by
Nima Salajegheh

has been accepted towards fulfillment
of the requirements for the

MS. degree in Mechanical Engineering

 

 

aUW

Major Professors Signature
8 / 18’ [o‘z‘

Date

MSU Is an Affirmative Action/Equal Opportunity Institution

._.-.-.--

.0-o-a-l-o-l-0-I-I-I-I-I-o-a.l-l-A-l-0-l-e-n-I-n-a-I-o-a-¢-o-l-I-I-I-c-a“a.

v V ‘VV V‘

 

 

LIBRARY
Michigan State
University

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c:/ClFiC/DateDue.p85p15

THEORETICAL DERIVATION AND NUMERICAL IMPLEMENTATION OF
CONTINUUM DAMAGE BASED CONSTITUTIVE EQUATIONS

BY

Nima Salajegheh

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Mechanical Engineering

2004

ABSTRACT

THEORETICAL DERIVATION AND NUMERICAL IMPLEMENTATION OF
CONTINUUM DAMAGE BASED CONSTITUTIVE EQUATIONS

BY

Nima Salajegheh

In this work, the isotropic damage is assumed and the theoretical
formulations of damage-coupled material behavior are presented
based on the internal variable approach. Additionally, several methods
for damage measurement and specification of material damage
parameters are introduced and an experiment is performed to
measure damage parameters qualitatively. Finally, the constitutive and
evolution equations are implemented into a commercial finite element
code, Ls-Dyna, and the results are compared with the non-damage
based analysis. Throughout this work, other common approaches to
Continuum Damage Mechanics and their differences from the current

approach are mentioned briefly.

To my Family

ACKNOWLEDGMENTS

I would like to thank my advisor, Dr. Farhang Pourboghrat, who
supported me throughout my study. I would also like to thank other
faculty members of the department of Mechanical Engineering at
Michigan State University who were helpful and willing to help me with

my research.

TABLE OF CONTENTS

LIST OF TABLES ....................................................................... viii
LIST OF FIGURES ....................................................................... ix
LIST OF ABBREVIATIONS ............................................................ xi
CHAPTER 1: Continuum Damage Mechanics .................................... 1
1.1 Introduction ..................................................................... 1
1.2 Phenomenological Aspects of Damage ................................. 3
1.3 Atoms, Elasticity and Damage ............................................ 3
1.4 Plasticity and Irreversible Strains ........................................ 5
1.5 Scale of the Phenomena of Strain and Damage ..................... 6
1. 6 Different Manifestations of Damage ..................................... 7
1. 6. 1 Brittle damage ............................................................ 7
1. 6. 2 Ductile damage ........................................................... 7
1. 6. 3 Creep damage ............................................................ 8
1. 6. 4 Low cycle fatigue damage ............................................ 8
1.6.5 High cycle fatigue ........................................................ 9
1.7 Mechanical representation of Damage ............................... 10
1.7.1 One-Dimensional Surface Damage Variable (L.M. Kachanov
[1]) 10

1. 7. 2 Three-Dimensional tensorial representation of the Damage
parameter (Murakami and Ohno, 1978 [8]; Murakami, 1988) 12

1.8 Effective Stress Concept (Y.N. Rabotnov, 1968) .................. 14
1.9 Equivalence Principles ..................................................... 16
1.10 Strain Equivalence Principle (1. Lemaitre, 1971) ................. 18
1.11 Coupling between Strains and Damage .............................. 19
1.11.1 Elasticity law .......................................................... 19
1.11.2 Plasticity Law ......................................................... 20
1.12 Rupture criterion: ........................................................... 21
1.13 Damage initial threshold .................................................. 22
CHAPTER 2: Measurement of Damage ([15] & [18]) ....................... 24
2.1 Direct Measurements ....................................................... 24
2.2 Variation of the Elasticity Modulus ..................................... 24
2.3 Ultrasonic wave propagation ............................................ 26
2.4 Variation of the Micro hardness ......................................... 27
2.5 Other Methods ................................................................ 29
2.5.1 Variation of density ................................................... 29
2.5.2 Variation of electrical resistance (potential drop) ........... 3O

2.5.3 Variation of the cyclic plasticity response (stress amplitude

drop) 31
2.5.4 Tertiary creep response .............................................. 33
CHAPTER 3: Thermodynamics Approach ....................................... 37
3.1 How Thermodynamics Approaches Damage ........................ 37
3.2 Thermodynamical Variables, State Potential ....................... 37
3.3 Potential of Dissipation, Damage Evolution Equation ............ 43
3.4 Plasticity ........................................................................ 46
3.5 Damage Coupled Constitutive Equations ............................ 46
3.6 Three Dimensional Damage Threshold ............................... 52
3.7 Three-Dimensional Rupture Criterion ................................. 55
CHAPTER 4: Identification of Material Parameters .......................... 59
4.1 Identification of the Material Parameters by standard tensile
test for Al3003 ........................................................................ 60
4.1.1 Tensile-Test Specimen Specifications ........................... 61
4.2 A Simpler representation of the Lemaitre’s Damage Evolution
Equation ................................................................................ 65
CHAPTER 5: Numerical implementations ....................................... 67
5.1 Uncoupled Analysis of Crack Initiation ............................... 67
5.1.1 Integration of the Kinetic Damage Law ......................... 68
5.1.2 Uncoupled analysis gives more conservative results ....... 70
5.2 Coupled analysis: ............................................................ 73
5.2.1 Strain coupled algorithm with the assumption of perfect
plasticity ............................................................................. 73
5.2.2 Strain coupled algorithm with the assumption of using the
damage value at the beginning of the increment ...................... 84
5.3 Numerical Results ........................................................... 89

vi

LIST OF TABLES

Table 1 : State variable and their associated variables .................... 39
Table 2 : Flux and their dual variables .......................................... 44
Table 3 : Tensile specimen dimensions in the English units ............. 62
Table 4 : Young’s modulus data ................................................... 63
Table 5 : Comparison of the time to failure between uncoupled and
coupled analysis for monotonic loading ................................... 72
Table 6 : Material properties for numerical calculations ................... 89

vii

LIST OF FIGURES

Figure 1 : Plastic strain due to dislocation movement ....................... 4
Figure 2 : Elementary damage by nucleation of microcrack due to an
accumulation of dislocations (after D. Krajcinovic) ...................... 5
Figure 3 : Cyclic tension-compression curves for low cycle fatigue of
stainless steel (after J. Duffaily) ............................................... 9
Figure 4 : Cyclic tension-compression curves for low cycle fatigue of
stainless steel (after J. Duffaily) ............................................... 9
Figure 5 : Isotropic definition of damage parameter (Apple
representation is after J. Lemaitre, 1975) ................................ 11
Figure 6 : Evolution of damage during low cycle fatigue .................. 33
Figure 7 : Evolution of creep damage in IN 100 super alloy (after H.
Poheefla) ............................................................................. 35
Figure 8 : More number of stars specifies a better method for damage
measurement ....................................................................... 36
Figure 9 : Damage evolution and material degredation ................... 47
Figure 10 : Measurement of damage by elasticity modulus .............. 59
Figure 11 : Stress-strain curve for the loading-unloading test .......... 62
Figure 12 : Plot of Damage versus strain ...................................... 63
Figure 13 : The front view of the uniaxial test model ...................... 90
Figure 14 : The front view of the biaxial test model ........................ 91
Figure 15 : The isoperimetric view of the biaxial test model ............. 91
Figure 16: Mesh configuration and location of the critical element for
single element analysis in the uniaxial test model ..................... 93
Figure 17 : Mesh configuration and location of the critical element for

single element analysis in the biaxial test model ....................... 94

viii

Figure 18 : Comparison of von-Mises stress for without damage and
coupled damage analyses for element number 700 ................... 95

Figure 19 : Damage evolution and comparison of effective (total)
plastic strain for analysis cases of without damage and coupled

damage ............................................................................... 96
Figure 20 : Conservative prediction of failure from uncoupled analysis

.......................................................................................... 97
Figure 21 : Initiation of Damage for the uniaxial model. (The right

figure shows the effective plastic strain) .................................. 99

Figure 22 : The Initiation of Damage for the biaxial model. (The right
figure shows the effective plastic strain) ................................ 100

Figure 23 : Damage and effective plastic strain start to grow from the
edges ................................................................................ 102

Figure 24 : Localization of damage and its comparison with effective
plastic strain for the uniaxial model ...................................... 104

Figure 25 : Localization of damage and its comparison with effective
plastic strain for the biaxial model ........................................ 105

Figure 26 : Evolution of damage (black) and effective plastic strain
(blue) of the element shown in Figure 17 for the uniaxial model106

Figure 27 : Evolution of damage (black) and effective plastic strain
(blue) of the element shown in Figure 18 for the biaxial model. 106

Figure 28 : Distribution of the von-Mises stress for the uniaxial model.
(The right figure shows the damage model result and the left
figure shows the non-damage model result) .......................... 108

Figure 29 : Distribution of the von-Mises stress for the biaxial model.
(The right figure shows the damage model result and the left
figure shows the non-damage model result) .......................... 109

Figure 30 : Damage model result for von-Mises stress in element
shown in Figure 17 for the uniaxial model .............................. 110

Figure 31 : Non-Damage model result for von-Mises stress in element
shown in Figure 17 for the uniaxial model .............................. 110

Figure 32 : Damage model result for von-Mises stress in element
shown in Figure 18 for the biaxial model ............................... 111

Figure 33 : Non-Damage model result for von-Mises stress in element
shown in Figure 18 for the biaxial model ............................... 111

LIST OF ABBREVIATIONS

0'- Cauchy stress
6- Effective stress acting on the resisting area
a” - Deviatoric stress

0,, - Hydrostatic stress

ay - Yield stress

a, - Ultimate stress

a,- Fatigue limit stress

0. - Yield stress for elastic perfectly plastic materials
aeq- von Mises equivalent stress

a— Cauchy strain
8‘ - Elastic part of Cauchy strain
eP- Plastic part of Cauchy strain

p - Equivalent plastic strain

0- Damage parameter

DC- Critical value of damage parameter in three dimensions
Dk- Critical value of damage parameter in pure tension

8- Damage strength material parameter

epo- Damage plastic strain threshold in pure tension

£R- Strain to rupture in pure tension

xi

Y- Strain energy density release rate

I7 - The associated variable of damage parameter

Y; Critical strain energy density release rate at rupture

,0- Mass density

,6- Mass density for the damaged material

m - Mass

E- Young’s modulus

E- Young’s modulus for the damaged material

v- Poison’s ratio

23- Poison’s ratio for the damaged material

A, ,u- Lame coefficient

2°.- Plastic multiplier

a, R, R”, 7, - Nonlinear kinematic hardening material parameter
K, n- Nonlinear isotropic hardening material parameters

r, R, R”, b- Nonlinear isotropic hardening material parameters
we - Elastic strain energy density

w*,- Complementary elastic strain energy density

t- Time

T- Temperature

s- Entropy density

v- Wave speed

xii

i3- Wave speed for the damaged material
V - Electrical voltage and volume

r- Resistivity

R; Triaxiality function

W' Helmholtz specific free energy

F - Function potential of dissipation

FD- Damage potential function

NR - Number of cycles to rupture

C- Elastic stiffness tensor and

5 - Effective elastic stiffness tensor for the damaged material
M - Damage effect tensor

H - Hardness

H - Step function

k’- Coefficient of proportionality for hardness test

R, r-Radius of cavities and radius of the yield surface

1- Length

i- Intensity of the electrical current

KP, M - Material parameters for fatigue
K,, N - Material parameters for creep
6,].- Kronecker delta

21’ - Heat flux vector

xiii

I, - Time to failure

xiv

CHAPTER 1: Continuum Damage Mechanics

1.1 Introduction

A new branch of science usually develops thus. Somebody publishes
the basic ideas. Hesitatingly at first, then little by little, other original
contributions appear, until a certain threshold is reached. Then,
overview articles are printed, conferences are held, and a first mention
is made in textbooks, until specialized monographs are written.
Continuum damage mechanics has reached that status now.

Professor Horst Lippmann
Lehrstuhl A fUr Mechanik
Technische Universitiit Miinchen,

Germany

Ductile fracture has been the subject of active research in the last
decades because of its importance to practical engineering applications
where the forming of a part introduces large deformations beyond the
elastic region. Because of its success in the analysis of brittle fracture,
a major effort has been devoted to the application of fracture
mechanics to characterize ductile fracture. This effort has led to the

introduction of different ductile fracture criteria which notably include

J-integral and COO (Crack opening displacement). These criteria offer
a tool for predicting the crack progression at the macroscale but they
fail to take into account the continuous deterioration of material
properties as the effect of microcrack nucleation and accumulation.
Moreover, fracture mechanics is based on the analysis of existing
cracks which might be too late to prevent a disaster. Therefore, the
nucleation stage, which consists of the evolution of internal damage
before macrocracks become visible, caught attention. As a result it
was highly appreciated by the scientific community when L.M.
Kachanov published in 1958 a simple model of material damage which
subsequently could be extended to brittle elastic, plastic or viscous
materials under all conditions of uniaxial or multiaxial, simple or cyclic
loadings. Recently, the theory of continuum damage mechanics has
been developed to a state ready for engineering applications. The
theory has been applied to solve a number of important engineering
problems including low-cycle fatigue in metals, high cycle fatigue,
creep-fatigue interaction, damage in composites, creep damage, and
ductile rupture. Also, the theory has been extended by the
development of both isotropic and anisotropic models of continuum

damage mechanics.

1.2 Phenomenological Aspects of Damage

The damage of materials is the progressive process which degrades
the material parameters continuously and causes them to fail. At the
microscale level this is the breakage of bonds between atoms due to
concentration of stress in the neighborhood of defects or interfaces. At
the mesoscale level of the representative volume element (RVE), this
is the growth and the coalescence of microcracks or microvoids which
initiate a crack and at the macroscale level this is the growth of that
crack. The two first stages may be studied by means of damage
variables of the mechanics of continuous media defined at the
mesoscale level which is the goal of Continuum Damage Mechanics
approach. The third stage is usually studied using fracture mechanics
with variables defined at the macroscale level. Continuum mechanics
and the thermodynamics of irreversible processes model the material
behavior without detailed reference to the complexity of their physical

microstructures.
1.3 Atoms, Elasticity and Damage

All materials are composed of atoms, which are held together by
bonds resulting from the interaction of electromagnetic fields. Elasticity

is directly related to the relative movement of atoms. When debonding

occurs, this is the beginning of the damage process. For example,
metals are organized in crystals or grains, a regular array of atoms
except on lines of dislocations where atoms are missing; if an
adequate shear stress is applied, the dislocations will move, thus
creating a plastic strain by slip without any debonding as shown in the

Figure 1.

 

 

 

 

 

 

 

 

 

 

 

II LII
ll
\
l

 

 

 

 

 

 

 

 

Figure 1 : Plastic strain due to dislocation movement

If the dislocation is stopped by a defect or a stress concentration site,
it creates a zone in which another dislocation may be stopped. The
accumulation of dislocations creates a debonding damage as shown in

Figure 2 and eventually a microcrack will form.

 

 

 

 

\\\ ///

Figure 2 : Elementary damage by nucleation of microcrack due to an accumulation of
dislocations (after D. Krajcinovic)

There are some other damage mechanisms in metals such as
intergranular debonding and decohesion between inclusions and the
matrix all of which create plastic strains.

In all cases, elasticity is directly influenced by the damage, since the
number of atomic bonds responsible for elasticity decreases with
damage. This coupling, which occurs at the level of the state of the

material, is called a state coupling.
1.4 Plasticity and Irreversible Strains

Plasticity is directly related to slips due to the movement, climbing,
or twinning of dislocations. Damage influences plastic (irreversible)
strain only because the elementary area of resistance decreases as the
number of bonds decreases. However, damage does not directly

influence the mechanism of slip itself; that is, there is no state

coupling. The indirect coupling owing to an increase in the effective
stress arises only in the rate constitutive equations and it is called

kinetic coupling.
1.5 Scale of the Phenomena of Strain and Damage

- Elasticity takes place at the level of atoms.

- Plasticity is governed by slips at the level of crystals or grains.

- Damage is debonding from the level of atoms to the mesoscale level
for crack initiation and it can be viewed as a bridge between
microscale and macroscale analysis.

Continuum mechanics deals with quantities defined at a geometrical
point. From the physical point of view, these quantities represent
averages on a certain volume, the Representative Volume Element
(RVE). The RVE must be small enough to capture high gradients but
large enough to represent an average of the processes at the
microscale[15]. To summerize:

- The microscale is the scale of the mechanisms of strains and
damage;

- The mesoscale is the scale at which the constitutive equations for
numerical analysis are written;

- The macroscale is the scale of engineering structures;

1.6 Different Manifestations of Damage

Even though the damage at the microscale is governed by one
general mechanism of debonding, at the mesoscale it can be seen in
various ways depending upon the nature of the material, the type of

loading, and the temperature [15]. Some of these manifestations are:
1.6.1 Brittle damage

The damage is called brittle when a crack is initiated at the
mesoscale without a large amount of plastic strain. This means that
the cleavage forces are below the forces that could produce slips but
are higher than the debonding forces. The degree of localization is

high.
1.6.2 Ductile damage

On the other hand, the damage is called ductile when it occurs
simultaneously with plastic deformations, larger than a certain

threshold p0. Damage is due to the nucleation of cavities followed by

their growth and their coalescence. As a consequence, the degree of

localization of ductile damage is comparable to that of plastic strain.

1.6.3 Creep damage

When a metal is loaded at elevated temperatures, for instance a
temperature above 1/3 of the melting temperature, the plastic strain
occurs at a stress lower than the yield stress. When the strain is large
enough, there are intergranular decohesions which produce damage
and an increase of the strain rate through the period of tertiary creep
[1]. Like ductile damage, the gradients of creep damage are similar to

the viscoplastic strain gradients.

1.6.4 Low cycle fatigue damage

When a material is subjected to a cyclic loading at high values for
stress or strain, damage develops together with cyclic plastic strain
after a period of saturation preceding the phases of nucleation and
propagation of microcracks. The degree of damage localization is
higher than those of ductile or creep damage. Because of the high
values for the stress, the low cycle fatigue is characterized by low

values of the number of cycles to rupture, NR < 10000. If the material is
strain loaded, the damage induces a drop of the stress amplitude as
shown in Figure 3 for two stress-strain loops corresponding to the

initial cycles and a cycle close to the rupture.

O'IMPoI h 0’
3“) -l

     

 

I

200 «

 

 

 

 

560 1066 . 1500 20m

 

 

. . . -- -.- . —-_-
m -n~——-.—c w.- ‘Iv-‘*-*-—- ~A
l-Q‘HQ‘ u. n...-

 

 

 

 

 

 

 

 

 

 

 

_........ —“""—"———--_—‘—T-——'——_""]

___——"‘J

...-_..1.. .-~_ .—.-._
--<~

 

Figure 3 : Cyclic tension-compression curves for low cycle fatigue of stainless steel
(after J. Duffaily)

1.6.5 High cycle fatigue

When a material is loaded with lower values for stress or strain, the
plastic strain at the mesolevel remains small and negligible. The

number of cycles to failure may be very large, NR >100000 and the

localization of damage is high and similar to that of brittle fracture.

Also, a drop of stress occurs similar to low cycle fatigue.

 

 

 

I .

 

 

Figure 4 : Cyclic tension-compression curves for low cycle fatigue of stainless steel
(after J. Duffaily)

Note that for brittle damage and high cycle fatigue damage where
damage localization is high, a stress-strain curve obtained from a
classical tension-compression test at the macroscale does not
represent the true behavior for strain and damage because the
localization induces plastic and damage zones much smaller than those
of the specimen. Nevertheless, it is used because mechanical
experiments at the microscale are difficult to perform; the results are
averages of nonuniform quantities over a mesovolume. The
microhardness test may help to characterize a microvolume as it
involves a size of the order of microns but its state of stress is

complex.
1.7 Mechanical representation of Damage

1.7.1 One-Dimensional Surface Damage Variable (L.M. Kachanov [1])

At the microscale, damage may be interpreted as the creation of
discontinuous surfaces, breaking of atomic bonds, and enlargement of
microcavities. At the mesoscale, the number of broken bonds may be
approximated in any plane by the area of the intersections of all the
flaws with that plane [1]. In order to work with a dimensionless
quantity, this area is scaled by the size of the representative volume
element, which is of great importance in the definition of a continuous

variable in the sense of continuum mechanics. At one point, it must be

10

the representative effect of microdefects over the mesoscale volume
element. It is similar to plasticity where the plastic strain

5,. represents, at one point, the average of many slips. If we consider a

Representative Volume Element (RVE) oriented along the direction ii
and a plane passing through the RVE with the same orientation (Figure
5), we can define:

- as as the area of the intersection of the plane with the RVE;

- 85,, as the area of the intersections of all microcracks or

microcavities which lie in as;

 

Figure 5 : Isotropic definition of damage parameter (Apple representation is after J.
Lemaitre, 1975)

- The value of the damage D(M,fi)attached to the point M in the

direction iiis defined as:

as
D M,” = D
( n) as

 

11

The damage parameter at point M is the maximum value of

D(M,ii) for all possible orientations,ii .

It follows from this definition that the value of the scalar variable D
(which depends upon the point and the direction considered) is
bounded by 0 and 1:

D = 0 —-> Undamaged RVE material;

D =1—->Fully broken RVE material in two parts.

In fact, the failure occurs for D <1 through a process of instability.

1.7.2 Three-Dimensional tensorial representation of the Damage

parameter (Murakami and Ohno, 1978 [8]; Murakami, 1988)

The theory of continuum damage mechanics was first developed for
material deterioration in the process of creep. Kachanov in 1958 first
proposed a phenomenological theory of creep damage by introducing a
scalar damage variable D to characterize the material degradation
which is responsible for tertiary creep and rupture [1]. The underlying
physical suggestion behind the scalar characterization is that creep
damage in a continuum is assumed to be everywhere the same and
independent of the specific orientation chosen. Thereafter, Rabotnov
interpreted that the damage variable represents the fraction of

damaged material, thus reducing the effective resisting area of the

12

local material element from S to So =(l—D)S. In the one-dimensional

case, the effective stress was thus defined as:

F = a
80—0) 0—0)

 

~ F
0':—:
50

Although the idea of scalar isotropic damage parameter is intuitively
appealing and has been universally accepted, it needs careful re-
examination. The main drawback of this isotropic model is that the
model predicts a constant value of the Poisson's ratio not affected by
material damage while the Young's modulus E is reduced to

E(1—D)after damage has occurred. However, experimental

observations have shown a change in Poisson's ratio. ([9], [11], and

[12])

On the other hand, despite the complexity of their mathematical
structure, tensors of second ranks can describe the spatial distribution
of microcavities more accurately, and hence, considerable efforts have
already been made to develop damage models based on second order

tensors [9], [10]. The second order damage field _D_(x) has the following

form in the principal directions:

D(x) =

cop
oNDo
poo

DUDZ, and D3 are simply the ratios of the damaged surface area over

the total surface area for any point in three orthogonal directions. It is

13

assumed that the principal directions of 2(x) coincide with the principal

directions of stress field at the specified material element because the
dislocations tend to line up and form the crack in the plane, which is
perpendicular to the principal stress direction[11].

In addition to scalars and second order tensors, vectors and forth
order tensors have also been used to reperesent the damage
parameter. ([3], [4], [5], and [6])

However since the differentiation with respect to damage parameter
will be necessary in the derivation of constitutive laws and tensors will
result in additional terms and calculations; it is very beneficial to avoid
tensorial representation whenever the material is not highly
anisotropic. Besides the increase in the cost of analysis, anisotropic
approaches introduce additional material properties which should be
determined from classical experiments such as tensile loading-
unloading test and measurements of Young's modulus and Poison’s
ratio along other orientations [11]. Throughout this work, the isotropic

scalar representation of damage is used.
1.8 Effective Stress Concept (Y.N. Rabotnov, 1968)

Under the loading ofF, the usual uniaxial stress is:

14

If all the defects are open and no force is carried by the broken area

~

represented by SD, an effective stress a is introduced which is related

to the surface that effectively resists the load, namely(S-SD):

 

 

 

.. F
0' :
S-SD
. . . . _. BSD
Introducmg the IsotropIc damage variable D(M,n) = as ,
.. F F
0' : S S = S
" D S l__?_
( S )
Or 6 = -——a
1—D

This definition is the effective stress on the material in tension. In
compression, some defects close and the surface that effectively

resists the load is larger than(S—SD), In particular, if all the defects

close, the effective stress in compression is equal to the usual stress.
Moreover, the damage evolution is almost zero and the damage
parameter remains constant. This effect is called the crack closure
effect and is taken into account by introducing a factor in the
numerical implementation of damage evolution equation. This factor
equals to zero when the hydrostatic stress is negative corresponding
to the material in compression and equals to one for positive values of
hydrostatic stress. However since the damage evolution in

compression may not be negligible for some materials and

15

temperatures, one can introduce a material parameter ranging from

zero to one which of course should be determined experimentally.

1.9 Equivalence Principles

A way to avoid micromechanical analysis for each type of defect
and each type of mechanism of damage is to postulate a principle at
the mesoscale. Equivalence principles are used to relate the damaged
state to the undamaged state. Several equivalence principles have
been postulated and their validity has been investigated [20]. The
choice of an equivalence principle is the first step in the derivation of
damage-coupled equations and somehow depends on the problem we
are looking at. For instance, if the material is highly anisotropic, the
tensorial representation of damage has to be considered and elastic
strain energy equivalence criterion has shown better agreement with
experimental observations. Whereas, strain equivalence criterion is
theoretically more reasonable in the case of scalar representation of
damage corresponding to isotropic materials [19].

Two widely used equivalence principles are elastic strain energy
equivalence and strain equivalence. The elastic strain energy
equivalence states that (Sidoroff, 1982; Cordebois, 1983; Chow and

Wang, 1987; Ju, 1989):

16

There exists a pseudo-undamaged material made of the virgin
material in the sense that its elastic strain energy is equal to that of
the damaged material except that the stress is replaced by the
effective stress in the elastic strain energy formulation.

The complementary elastic energy of an undamaged material (D = 0)is:
1 T -1
we(0',0)=-2-0' :C :0'

Then the complimentary energy of a damaged material (D¢0)is

expressed as:

l

w*e(0',D)=w*e(6,0)= 6'7:C":&=30'T:(MT:C":M):0

l
2
Where C is the elastic stiffness tensor and 5 can be defined as
effective elastic stiffness tensor such that:
C" =(MT :C'1 :M)
Mis a fourth order tensor called the damage effect tensor and is
defined such that:
5' =M :0
For the general case of anisotropic damage, damage effect tensor is

equal to [10]:

17

 

 

 

 

 

 

 

 

 

—l— 0 0 0 O 0
1—Dl
O l 0 O O O
I—D2
O O l O O O
1—D3
M = 0 0 0 1 O O
Ja—D2XITDs)
l
O O 0 0 O
\Kl—Djl—Ds)
l
O O O 0 O
L J(l—DZX1—Dl)d

Throughout this work, the strain equivalence postulate is used which is
explained hereafter. Reader can find the details about elastic strain
energy equivalence criterion and other equivalence criteria in [11] and

[20].
1.10 Strain Equivalence Principle (J. Lemaitre, 1971)

At the microscale, this postulate assumes that the constitutive
equations for the strain of an element are not modified by a
neighboring element containing a microcrack [15]. At the mesoscale,
this means that the strain constitutive equations written for the

surface (S -S,,)are not modified by the damage or that the true stress

loading on the material is the effective stress 6 and no longer 0.

Consequently, any strain constitutive equation for a damaged material

18

may be derived in the same way as for a virgin material except that

the usual stress is replaced by the effective stress [15].

8(0', D) = £(6,0)
a

6':—
l-D
This principle is applied to both elasticity and plasticity.

1.11 Coupling between Strains and Damage

Applying the strain equivalence principle, we can write the uniaxial

laws of elasticity and plasticity of a damaged material as following:

1.1 1.1 Elasticity law

Using the concept of effective stress:

as
II
DulQI

And strain in other directions for isotropic damage:

C

822 =€33e=—U€IIC
E and v are the Young's modulus and the Poisson's ratio of the

undamaged material. As a result, the elasticity modulus of the

damaged material can be defined by the ratio% which is E“ = E(l—D)
8

This inspires a method for evaluating the damage parameter based on

the change in Young's modulus:

19

1.1 1.2 Plasticity Law

In order to model plasticity two kinds of strain hardening are
usually considered:
-The isotropic hardening which indicates the size of the yield function.
-The kinematic hardening or the back stress which indicates the center
of the yield function.

If try is the yield stress, R the stress due to isotropic hardening and

X the back stress, both functions of the plastic strain, the one-
dimensional plasticity criterion defining the current threshold of yield
limit is:
0‘ = 0y + R+ X
Or:
f =|0—X|—0’y —R=O
f is the yield function from which the kinetic constitutive equation for

plastic strain is derived.

f=0
35?" #01f and

f=0

And:

20

f <0
é‘”=0if or
f <0
When damage occurs, according to the principle of equivalence, the

yield function f must be written as:

0'
=———X—a —R=0
f l—D I ’

1.12 Rupture criterion:

The rupture at the mesoscale is a crack formation which occupies
the whole surface of the RVE; that is,D =1. In many cases this is
caused by a process of instability which suddenly causes the
separation of remaining area. This point of instability corresponds to a

critical value of damage D, <1 which depends upon the material and

the conditions of loading.
The final decohesion of atoms is characterized by a critical value of the

effective stress acting on the resisting area. If the ultimate stress a; is

taken as the critical value of the effective stress, we have:

 

.. a
0': =0,
l—Dc
Then:
D z1_£_

21

This gives the critical value of the damage at a mesocrack initiation

occurring for the one-dimensional stress, a. The ultimate stress 0'“
being identified as a material characteristic, D, may vary between
D, zo for pure brittle fracture to DC =1 for pure ductile fracture but
usually D, remains in the order of 0.2 to 0.5.

This relation, applied to the pure monotonic tension test, which is

taken as a reference, defines the corresponding critical damage DI.

considered as a material characteristic:

Where a, is the stress at rupture.
1.13 Damage initial threshold

Before the microcracks are initiated, creating the damage modeled
byD, they must nucleate by the accumulation of dislocations. In the
pure tension case, this corresponds to a certain value 8p» of the plastic
strain below which no damage occurs:

8” < app —> D = 0
Finally, the four main relations which comprise the basis of Continuum

Damage Mechanics are:

-s‘ = a for elasticity
E(l—D)

 

22

‘f =|—1—0—D— XI —c7y -R = 0 as the plastic yield criterion
-6” < app —> D = 0 as the damage threshold

- D = D, —> crack initiation

23

CHAPTER 2: Measurement of Damage ([15] & [18])

2.1 Direct Measurements

This method consists in the evaluation of the total crack areas (18,,

lying on a surface as at mesoscale. This can be done by observing
microscopic pictures. The main drawback of this method is that it is a
destructive method and one sample is needed for each data point.
Also, it is tedious to practice since the sample needs to be taken out at

various strain values and observed under a microscope.
2.2 Variation of the Elasticity Modulus

This is a non-direct measurement based on the elasticity law:

0 = a
E(l—D)

 

8

It assumes uniform homogeneous damage in the specimen gauge
section.

If E=E(l—D)is considered as the effective elasticity modulus of the
damaged material, the values of the damage may be derived from
measurements ofE, provided that Young's modulus Eis known:

~

0:1—5
E

24

This very useful method requires accurate strain measurements. Strain

gauges are commonly used and E is most accurately measured during
unloading [18]. This technique may be used for any kind of damage as
long as the damage is uniformly distributed in the volume on which
the strain is measured which is the main limitation of the method. If
the damage is greatly localized, as for high cycle fatigue of metals, for
example, another method must be used.
At the beginning and at the end of the unloading paths in the
plane(a,e), there are small nonlinearities, owing to viscous or
hardening effects and also due to the experimental devices as well. It
is best to ignore them and to identify E from the slope of the middle
points. Also, it is most important to always use the same procedure to
evaluate Eand the evolution of E.
For damage in metals, an early decrease of if at low strain levels may
happen which is due to movements of dislocations, and to texture
development, but not to damage [15]. As this phenomenon rapidly
saturates, it is common to consider a damage threshold such that:

8” < epD —> D = 0
Instead of strain gage, one may use grid patterns and observe the
deformation of the grid pattern. After the deformation, a circular
pattern may look like an oval whose major and minor diameters and

their orientations can be used to obtain principal strains and their

25

orientations respective to the loading direction. To avoid taking out the
sample at each unloading increment in order to measure the strains,
one may take pictures from the grid pattern while the sample is still in

the MTS machine and analyze the pictures when the test is complete.
2.3 Ultrasonic wave propagation

This method is based on the effect of change in elasticity modulus
on the speed of ultrasonic waves. For frequencies higher than 200 kHz

the longitudinal wave speed v, and the transversal wave speed v, in a

linear isotropic elastic cylindrical medium are [18]:

2 _£ 1"”
’” p(l+v)(l—22))

 

, E 1
v 1' =—
p2(l+v)

 

Where E is the Young's modulus, p the density, and v Poisson's ratio.

A measurement of the longitudinal wave speed of a damaged material

gives:

l—v
(1+ v)(l- 21))

~2

E
VL=-:
p

 

Where Eand ,5 are the actual damaged elastic modulus and density.

Poisson's ratio does not vary with damage if elasticity and damage are
isotropic. If anisotropic damage is considered, the Poisson's ratio for

the damaged material is different. The damage is calculated by:

26

If the damage consists mainly of microcracks or if small cavitation is

considered, 2 land:

‘0 I‘m

 

To measure the speed v,, ultrasonic transducers need to be placed on

both sides of the sample and damage parameter is assumed to be
uniform throughout the sample. Thus, the sample's size must be of an
order of magnitude coherent with that of the RVE and for metals, this
size is too small in comparison with ultrasonic transducers size and
accuracy of time measurements. Also, this method is destructive
because we need to cut the material into very tiny parts. Nevertheless,
this method gives good results for materials with larger RVE size such

as concrete [15].

2.4 Variation of the Micro hardness

This method is a non-direct measurement based on the influence of
damage on the plasticity criterion. In the uniaxial state of stress, the

yield function is written as:

 

f: —X—0'),—R=O

 

27

As microhardness is a process of very small indentation, this method
may be considered as practically nondestructive. The test consists in
inserting a diamond indenter in the material and the hardness His
defined by the mean stress:

Hzazi
S

The load F on the indenter is adjusted so as to obtain a projected
indented area S of the same order of magnitude as that of the RVE.
Theoretical analyses and many experimental results prove a linear
relationship between H and the plasticity threshold 0,, k’ being the
coefficient of proportionality:
H = k’as

This threshold corresponds to the actual yield stress,

a, =(c7y +R+XX1—D)
Then:

H =k’(a'y +R+X)(1-D)
In fact, the microhardness test itself increases the strain hardening by

an amount which corresponds to a plastic strain a” of the order of 5 to

8% [18]. His then always related to (away), a” being the current

plastic strain.

28

If H‘ =k’(ory +R+ X) is the microhardness of the material which would
exist without any damage f0r(£”+£,,”), and H =k'(0‘y +R+XX1-D) the

actual micro-hardness for (5" +e,,"), then:

D=l- H,
H

 

2.5 Other Methods

Several other methods have been introduced based on the influence
of damage on some physical or mechanically measurable properties

such as:
2.5.1 Variation of density

In the case of pure ductile damage, the defects are cavities which
can be assumed to be roughly spherical. This means that the volume
increases with damage and the corresponding decrease of density is

measurable with devices based on the Archimedean principle.

If (igihs the relative variation of density between the damaged state

,b' and the initial non-damaged state ,0, it is easy, considering a

spherical cavity of radius rin a spherical RVE of initial radius R and
mass m, to derive the following relationships between the surface

damage D and the variation of density or porosity:

29

 

 

 

 

p:
in?
3
.. m
p=4
-—7r(R3+r3)
3
fi_p R3 __r3
= 3 3_1_ 3 3
p R+r R+r

 

2.5.2 Variation of electrical resistance (potential drop)

Ohm's law for a non-damaged element of length 1, area S and

resistivity r is written as:
V = r—i

Where i is the intensity of the electrical current. Whereas, for a

damaged element of the same size:

 

 

Where 2":1—l— is the effective intensity of the electrical current

defined in the same way as the effective stress was defined, using the

surface definition of damage.

7 is the resistivity affected by the damage. Bridgman's law [18] gives

~

r as:

30

7=r[I+K—p'p]=r(1+KD%j
p

K is a coefficient which is approximately 2 for metals.

If the same intensity 1‘ is considered for the non damaged and the
damaged elements, the damage D may be derived from the two

expressions:

 

 

. F
For small values ofD, the correction term — due to the volume change
r

~

is close to 1 (for instance, D = 0.1 —> I— = 1.064); then:
r

ozI—L
V

This method is known as the "potential drop method".
2.5.3 Variation of the cyclic plasticity response (stress amplitude drop)

The influence of damage on plasticity may be used to measure the

low cycle fatigue damage in which the material undergoes plastic cyclic

loading.

31

The one-dimensional law of cyclic plasticity at stabilization may be

written as a power relationship between the amplitude of stress A0

M
and the amplitude of plastic strain Ag” as A5” =[%] for a non-

P

M
damaged material, and as A8” =[ A0 J for a damaged material

Kp(l—D
according to the principle of strain equivalence, where KP and M are

material parameters.
Considering a test at constant plastic strain amplitude, if A0* is the
stress amplitude before the beginning of the damage process and A0

is the stress amplitude after the beginning of the damage process, we

p _ Aa* M _ A0 M
A8 [76'] ‘I—K.<1—DII

From which it follows that:

have:

A0
110*

D=l—

 

This method successfully identifies the evolution of damage during low

cycle fatigue of metals. An example is given in Figure 6.

32

12 .-

 

I023 —— 10 k

Ad alOZIMPo

 

 

 

10 2 15 20
N =10

 

I
ll
’ l
'T— one mesocrack

/ develops
I —>

 

 

IO 15 20

Figure 6 : Evolution of damage during low cycle fatigue

2.5.4 Tertiary creep response

Creep damage occurs in metals loaded at temperatures above

approximately one third of the melting temperature. If we apply the

principle of strain equivalence to Norton's law of tertiary creep.

 

N
éPflt: 0‘
K.

33

K, and N being temperature dependent material parameters.

Assuming that the damage process begins at the end of the secondary

creep, during tertiary creep one may write:

Him—:71)"

w... %v
Mr)
8

Where .6” * is the minimum creep rate.

From which one derives:

 

This method yields good results which are in accordance with those
obtained by measuring the variation of the elasticity modulus. An

example is given in Figure 7.

34

a. We)

 

 

tIhI

 

 

ll ;
1.0r- ...;..........l
' I
0.8 ~— I:
O : I
0.1. "" I I
o 2 I. 4— one mesocrack
' develops
l i E I >
0 10 20 30

NM ‘R
Figure 7 : Evolution of creep damage in IN 100 super alloy (after H. Poheella)

As a conclusion, Figure 8 offers some advices for choosing the proper
method of damage measurement, depending on the kind of damage

involved. [18]

35

 

 

 

Low High
cycle cycle
Damage Brittle Ductile Creep fatigue fatigue
, , USO
MIcrography I) = ‘iS‘ * M H * *
r .
i) 2/3
Density D : (I — -) *3“ t a
P
EIBSIICII)’ D Z I _ E ** $4,... *4”; #1::
modulus E
. ‘wz
UIII‘BMImC D : l _ 1:17;- *** *1“ *g s: *
waves If
Cyclic slress D = l _ A0 a, a: 32* =l=
amplItudc 430"
_ _. l/N
Ten'dry D _____ I __ :13 :I: *3“: as
creep 5p
_ H
MIcro-hardness D = l —- I—F ** *** H *M *
Electrical D = 1 __ L ,.. a... n * *
reSIstance V

Figure 8 : More number of stars specifies a better method for damage measurement

36

CHAPTER 3: Thermodynamics Approach

3.1 How Thermodynamics Approaches Damage

The thermodynamics approach to obtain damage equations in three
dimensions is to postulate the existence of energy potentials from
which one can derive the state laws and the kinetic constitutive
equations. In the thermodynamics of irreversible processes, two
potentials are introduced [7]:
1-The state potential, a function of the state variables, which defines
the state laws and the variables associated with the state variables to
define the power involved in each physical process. For the damage
variable, an energy damage criterion is derived from the damaged
elasticity potential.
2-The potential of dissipation, which is a function of the associated
variables accounts for the kinetic laws of evolution of the state
variables including damage. A constitutive equation for the damage

gives the damage rate as a function of its associated variable.
3.2 Thermodynamical Variables, State Potential

The state variables can be categorized as follows:

- Observable variables such as:

37

8- The strain tensor of components suassociated with the Cauchy

stress tensoro

T- The temperature associated with the entropy density s

a- The back strain tensor, whose associated variable is the back
stress X D and X” represents the kinematic hardening, translation of
the center of the yield surface in the deviatoric space

r- Which is equal to effective plastic strain when damage is absent
and its associated variable is R representing the isotropic hardening.
D- The damage variable. If the damage is considered to be isotropic,

65,,
as

 

it has the same value in all directions, and the scalar D=

characterizes completely the three-dimensional state of damage in the
RVE at the point considered. 17 is the associated variable of D which
will be derived from the state potential. As D is dimensionless and
since the product - 17D is the power involved in the process of damage,
it can be seen that 17 is a volume energy density.

The following chart summarizes the variables involved in the

continuous mechanics of elasticity, plasticity and damage:

38

 

 

 

 

 

 

 

Associated
State Variables
Variables
Observable Internal

8 0’
T s

r R

a X D

D I7

 

 

 

 

 

Table 1 : State variable and their associated variables

It is postulated that the state laws are derived from a state potential;
a continuous scalar function, concave with the temperature, convex
with the other state variables and containing the origin [7]. Taking the
Helmholtz free energy:
‘P = ‘Y(£,T,£‘,£",r,a,D)
The strains act only through their difference 5‘ = 23—5”:
‘I’=‘I’(£‘,T,r,a,D)

The associated variables can be derived as:

39

B‘I’
a:p3€
_8‘P

5?
_ 8‘1’
7?

X0:

 

R

ii
paa

_ 6‘?
Y: —
pan

The analytical expression for ‘Pis chosen considering the experimental
observations and:

- The linear isotropic elasticity;

- The state coupling of damage with elastic strain and the principle of
strain equivalence together with the concept of effective stress written

for the three dimensional case as:

.. _ 0.,(t)
Mai-0(1)

For anisotropic damage, the damage variable is no longer a scalar and
the effective stress needs the application of damage effect tensor. The
difference of the behavior in tension and compression is not taken into
account here either. For isothermal loadings:

l l . . 1 — r X»?
=;[3Cmd5 I18 “(l—DH'R... [’“I'geI b )]+_3_aijau]

The law of elasticity coupled with damage is:

8‘? ,
0' =pa—8e-z CUHE “(l—D)

40

Where C is the fourth order elasticity stiffness tensor or, by inversion

for the isotropic case,

 

E being Young's modulus, v Poisson's ratio and 5,, the Kronecker

delta. The isotropic strain hardening scalar stress is expressed as:

R = p%\§-= R,[1—e“’”’]

Rwand b are two parameters which characterize the isotrOpic strain

hardening for each material. Xmand y are two parameters of nonlinear
kinematic hardening written with the factor if to ensure a simple

expression in one dimension. The associated variable for D is defined
as:

—_ 3‘P_ , ,
p _-_ I'jldgijekl

— BD 2
In order to work with a positive quantity,Y is introduced as:
Y=—Y
A relationship can be found between Y and the elastic strain energy

density we defined as:
dwe = aqde‘.)

Integrating with the law of elasticity, and assuming no variation of

damage, that is, D = const. , yields:

41

e C l C C
we = ICWE “(I—pm ,1. =§CW£ Us “(l—D)

This shows that:

Y: we
l-D

 

Yis also equal to one half the variation of the strain energy density

with damage at constant stress, do, =0, [7]. Starting with the law of

 

elasticity:
do, = CW [(1— D)d£‘,e, — e‘udD] = 0
Or:
dEeu = Eeu d0
1— D
Together with:
. . dD . , dD
dwe(a=ronsr) : Gilda i} = 0178 Id :3 = ("111:1E I] (l — D>8 kl i_—
Or:
dwc(a=romt) e e
T: gut" .78 It
And from the definition of Y =—17
Y = _1_ dwe(a’=ronsl)
2 dD

Because of this equation, Y is called the strain energy density release
rate. This is the energy released by loss of stiffness, damage

softening, in the RVE when damage occurs.

42

3.3 Potential of Dissipation, Damage Evolution Equation

Having defined all the state and associated variables, a second
potential is needed to describe the evolution of the phenomena. First
we investigate the second principle of thermodynamics to ensure a
non-negative dissipation. Starting from the Clausius-Duhem inequality

[21].

| NJ

O'ejéej—p(‘I’+sT)—qe. ‘I 20

"I

E; is the heat flux vector associated with the temperature gradient.

If we take the derivative of free energy with respect to the state

variables, we have:

a? é‘--+3\£T+-a—q:r'+ 3‘? ct +£D
age. " aT 8r 3a,, ‘1' an

 

 

Together with 5' = é‘ + 3” it becomes:

3‘? 3‘11 - 3‘? 3‘? 3‘1’ - T.
0'.. - é‘r - +— T+ ”6"”;- - — '- ——d.. - —D- .420

 

Or, with the definition of the associated variables,

., . D . —- T,-
O'e-j-E .j-Rr-X gag—YD—qe?20
To always satisfy this inequality of a positive dissipation and
particularly for an isothermal process where the plastic dissipation, the

first term, is negligible, we must have:

—YD20

43

As —17 is a positive quadratic function,—172.0, and the damage rate D
is non-negative, the Clausius-Duhem inequality is satisfied. Notice that
the dissipation inequality is the sum of the products of rate or flux
variables (with the negative sign) multiplied by their dual variables as

shown by Table 2:

 

 

 

 

 

 

 

flux variables dual variables
6" a
—Ri' R
_a X0
-1) — 17
- T,
q 7'

 

 

Table 2 : Flux and their dual variables

It is postulated that the kinetic laws are derived from a potential of
dissipation, a scalar continuous and convex function F, of the dual
variables and the state variables may act as parameters [15]. For the
isothermal case:

F(0,R,XD,Y:(£‘,r,a,D))
The laws of plasticity coupled with damage are derived from this
potential by means of a scalar multiplier which is always positive. This

ensures the normality condition of yielding for plasticity.

44

 

 

 

. e, z __ II.
" 80
. 3F - . . . . .
r = “SR—xi i plasticzty constitutive equations
a _ _ 8F
ax”

And for damage:

15:15.21
BY

In this work, we use the von-Mises criterion in the yield function which

identifies the equivalent stress as:

3 x
0' =(EO'DUO'DIJ‘)

“1

Together with kinematic hardening, the von-Mises criterion is applied
to define the size of the yield surface regardless of the

translation X ”and it acts upon the difference between avg-and X D.

Furthermore, In the presence of damage, the coupling between the
damage and the plastic strain is written in accordance with the
principle of strain equivalence. The yield criterion is written, in the
same way as for a non-damaged material except that the stress is
replaced by the effective stress, which, for isotropic damage is:

0' ~
_=0'
l-D

Then, the yield function f is written as:

45

 

De
Th

 

 

 

 

 

 

f=|00_xveq_x_a,=o
With:
V
lO'D-XD 232:”; gpgl_xi.0 2
«I 21—D l-D ’
3.4 Plasticity

Plastic strain occurs when the yield criterion, f =0, is satisfied. The
plastic strain continues to grow if the yield criterion is continuously
satisfied, that is, if f = 0. Then plasticity deals with these two

conditions, which define loading, or unloading with f =0:

i=0”
and Waco
f°=0.
f<OI
or >é”=0
f.<OJ

 

xi is deduced from these two conditions:

f=0andf=0.

3.5 Damage Coupled Constitutive Equations

The choice of an analytical expression for the two potentials and
particularly for the potential of dissipation is of great importance.

Thermodynamics provides the general framework and some

46

restrictions on the functions that can be used, but only experiments
can give the details and verify the model. As the constitutive equations
must be as general as possible and cover a wide range of materials,
the general trends of basic experiments are considered and only the
value of a few parameters are taken to be as material parameters.

A schematic test in tension with some unloadings and compressive

loadings is shown in Figure 9:

 

 

 

 

‘0’
Im
Il'T-T'I"“ \
0’ I I I II\
I I I “y i E<E I .
El I._e__I_____~ [I \
o . I I Ix I .1 >8
I l i A0 /
Ii 1 I I
l l I I x
I; ,. -4
0 as?
EPD

 

Figure 9 : Damage evolution and material degredation

Below a certain value of the plastic strain, a threshold app, no damage

occurs. This is implemented in the code by comparing the equivalent

plastic to 5pc at the end of each increment. The kinematic back

stressX, defined as the locus of the center of the elastic domain,

47

increases with the plastic strain, is nonlinear with the plastic strain,

and tends to saturate to some value X“. The evolution of X is chosen

as:
X = Xw[1—eI'”PI]

Xe, Ancl y are parameters to be identified for each material and each
temperature. The isotropic hardening streSSR, identified as:

R = a- 0‘, - X
- Increases with the plastic strain;

- is nonlinear with the plastic strain;

- tends to also saturate to some value R...

Like XThe evolution of R is chosen as:
R = Rm|:1-eI’b‘PI]

If we use the exponential representation for hardening components,
the potential of dissipation is proposed to be [15]:

3

D D D D
F=|U -X —R—0'y+4—X—Xij Xi} +FD(Y:(r.D))

 

‘4

However, since the last term,F,,, is used to derive the damage

evolution law, other hardening equations such as the power law can be

used in the numerical procedure as well. The first three terms

48

correspond to the yield function f. The fourth term Zi—XUDX ”is

I}

responsible for the nonlinear kinematic hardening.

The choice of the function FD is of course the key to representing the

damage evolution and is based on physical observations and
experimental data.

Damage is always related to some irreversible strain. This property is
taken into account in the damage law by the multiplier/i, which is
proportional to the accumulated plastic strain.

- 3F - 6F
D=__£,t=—D ' l—D
air 81' p( I

The variable pshows the irreversible nature of the damage, as p is
always positive or zero.

As the accumulated plastic strain increases from zero the damage
remains equal to zero during the nucleation of microcracks and
accumulation of dislocations. A one-dimensional damage threshold

related to the plastic strain seep has already been introduced. As the

equation of damage is governed by the accumulated plastic strain and

as pzs" in one-dimensional loading, it is needed to introduce a

threshold on the variable p such that:

If prD:
- 3F - 3F
D=——_2,t=—D'1—D
air 61/ p( I

49

If p<pD:
D=o

This allows us to introduce a step function in the potentialFD:

PZPD_)1
H(p-p):
D P<PD—)0

In the next section, an approach is introduced to identify prrom the
uniaxial damage threshold app.
On the basis of a thermodynamical analysis, the dual variable for the
damage is the strain energy density release rateY. Hence, FD must be
a function on:

FD = FD(Y,...)
Another important feature of fracture is the influence of the triaxiality

0H

 

ratio , (0,, is the hydrostatic stress, 0... and is the Von Mises

at
equivalent stress). The modeling of this effect is contained in the
expression of Y by the triaxialitv factorRe :

2
Re — 2(l+v)+3(l—Zv)[a”]
U

3 eq

 

~2
0' .4

Y:
2E R”
_ 0'2“!
25(1—0)2 '

 

50

In order to choose the proper and simplest expression forFD, a general
trend in the damage evolution obtained by micromechanics for
particular mechanisms is considered as a guide [15] which show:
[)sz
Thus:
FD = Y2
Here, as in every constitutive equation, a scale factor must be
introduced. 80—0) is taken as the scale factor where Sis a material
constant. The term (1-D)is considered here to be cancelled, with
(l—D) coming from2= p(1—D) because experiments show a non-

decreasing damage rate when Yand pare constant:

2
FD = Y
S(l-D)

 

Finally, according to the quantitative properties listed above, the
damage potential is written as:

Y2
meo» = mHip-p.)

The factor V2 is used here to avoid 2 in the derivation:

BFD 1 er), . Y

D:- _
aY {W 50-0)

 

 

13(1 _ D)H(P'PD)
Thus:

D : —pH(P‘PD)

51

With the rupture condition for crack initiation,
D = D,
Three material parameters are introduced to characterize the damage
evolution:
-S the damage strength,

- p0 the initial damage threshold,
-D, the critical damage,

The effects of temperature can be taken into account by the variation

of these coefficients and yield function coefficients with temperature.
3.6 Three Dimensional Damage Threshold

The damage threshold pp in three dimensions which is equivalent

to spDin one dimension corresponds to the initial threshold of damage

evolution. There has not been a universally acceptable way for

derivation of po and in many cases 5pc is used with very satisfactory

results. Nonetheless, one of the most dependable criteria is developed
by Lemaitre [17] and is cited in this section which is related to the
amount of energy stored in the material as a result of plastic
deformation.

The total plastic strain energy dissipated may reach tremendous

values before failure but it is postulated that the stored energy

52

remains constant at the beginning of damage evolution. This stored
energy is the result of stress concentrations which develop in the

neighborhood of dislocations in metals. For a unit volume, it is equal to

the difference between the total plastic strain energy javépudtand the
0

energy dissipated in heat given by the Clausius-Duhem inequality of
the second principle of thermodynamics. The dissipated power for an
isothermal deformation is:

above”, -Rp-XD.-,-d,j 20
By replacing the rate variables with their definition from the potential

of dissipation:

 

8F.-
-=___g=g
p 312

_ 9}: 31: DnaF .>
¢—[0”80+R8R+X .,——axD)p_O

53

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

( “0'00 -X D ( O'Dij __ “D W
(p: 003 1;D " -R+X”.,— _3 1;D ” f 3 Xfio p
2 0' ., -X..D 2 0' ij -X..D 2X,

k l-D " (q k l—D " (q U
00,-,- D 3

 

 

 

D D .
¢=[1_D U —R+2—X-—XU Xij ]p
«I

3 D o .
¢=[0'y+§—X—XU Xij JP

Then, the stored energy was a function of time is:

' . ' 3 o D .
W"'):£0"8P'jdt—£[a’+§X—XU Xij )pdt
If we neglect the effect of the kinematic hardening,

0..
x”,- =0—>é”.-,- =%IE-,l=3‘I " p

20'“,

 

w :IMpdt_JO-ypdt
0 0

Or with:

It is assumed that this energy is equal to that of a perfectly plastic

material of plastic threshold 0', Zaywhose stored energy is a function

of the difference between a,and the fatigue limita, [17]. Then, as

0'“, = const.

54

W5 = (deg _0f )p
If we put this energy equal to that of pure tension reference case

having a damage threshold spDand a plastic threshold on, the ultimate

stress:

(Uta ‘01)170 :(au ‘0'] )EpD

 

0' - 0'
pD : EpD u _ f
0:4 0!

We have to remember that in this formula the material is considered

to be perfectly plastic. Then for a varying loadingama), a more
accurate calculation consists in performing the integration forw,. After

all these assumptions, it is fair to say that the derivation of three
dimensional damage threshold calls for more investigation.

Nevertheless, since acqis close to a, but slightly smaller, taking p0 = app

is a good approximation and will be used in the numerical analysis.

3.7 Three-Dimensional Rupture Criterion

Similar to the three dimensional damage threshold, the expansion
of one dimensional critical value of damage to three dimensions is
subjective and there is no universally acceptable method to obtain it. A
one-dimensional rupture criterion is obtained from the tensile loading-
unloading test. The damage parameter is calculated from the change

in the elasticity modulus and the critical value of damage is the final

55

value of damage parameter. For the general case of loading, a critical
three dimensional damage value is needed for numerical prediction of
rupture. In so many numerical analyses, the critical value of damage is
taken to be one in both one dimensional and three dimensional loading
situations. Nevertheless, experimental data show that the material
fails at the critical value much less than one and this value decreases
as the ductility decreases. One of the approaches for obtaining a
reasonable value for three-dimensional critical damage has been
proposed by Lemaitre [15] and is cited hereafter.

This approach assumes that the rupture moment corresponding to the
critical value of damage is governed by the amount of energy

dissipated in damage growth:

DC
I YdD = const.
0

For the case of perfect plasticity in proportional loading considered

here for simplicity:

 

 

Dc Dc 0-2 Dc 0.2 0.2
ijD=I—-i—2-RvdD=I ‘RvdD= ‘RVDC
0 0 2150-1)) 0 2E 25

This quantity should be equal to that of the uniaxial reference case

which is:

Y _ 02R _ 0214
2150—1)“)2 2E

 

56

 

 

 

 

 

 

 

And since a, z a“:

 

This formula shows that the critical value of the damage in three

0'”

 

dimensions decreases as the triaxiality ratio, contained in the

‘4

RV expression, increases. Since, R, is usually very close to one, taking
0,:ch is a good approximation and will be used throughout the

numerical part of this work.
The set of damage equations used in this work can be summarized as

follows:

D=D

C

D" —> Crack

The material parameters are:

- S, apDand chwhich must be determined from damage

measurements in tensile loading-unloading test;

57

- 0'“ and a, which are classical characteristics and easy to find in

handbooks or to identify by tensile and fatigue tests.

58

CHAPTER 4: Identification of Material Parameters

The material damage parameterss , epDand D,,which characterize

the damage must be measured for each material and temperature by
the means of tensile loading-unloading test. Let us assume that a good
tensile test has been performed with measurement of the damage
during unloading by elasticity change and the following curves have

been obtained:

 

 

 

 

Figure 10 : Measurement of damage by elasticity modulus

For the one-dimensional case:

2

Y=—3‘——2 as R, =1
25(1—0)

(w

Thus:

59

 

= msg— D)2 IEPIHV'TP")
The parameter S is determined from the slope of the curve, D versus
the plastic strain 5p:
= “2 2M
2ES(l—-D)
Or:

dD _ 0'2
dEP 2ES(1- D)2

 

 

At each point of the curve, D is known, ais known from the stress

strain curve, Edgis estimated and E is known from a previous
6‘

idenfiflcaflon:

0.2

2E(l-D)

S:

 

2 d0
d8”

 

The main difficulty is obtaining a good stress strain curve in the
softening range where necking occurs and strains must be measured

locally by a small strain gauge.

4.1 Identification of the Material Parameters by standard

tensile test for A|3003

In the following experiment, an extensometer was used to measure

the strain values. Extensometers assume uniform distribution of strain

60

as well as damage over their gage length and can not capture the
localization and necking of the specimen during which almost the
entire process of damage evolution occurs. To obtain accurate data,
small strain gages should be used to obtain strains at the much
localized points. Nonetheless, the method for derivation of the material
damage parameters is explained hereafter using the data acquired
from the extensometer keeping in mind that the calculated values are
just for the purpose of describing the procedure and are not
quantitatively accurate. For numerical analyses, more accurate

parameters are used which have been obtained by Lemaitre [14].
4.1.1 Tensile-Test Specimen Specifications

The tensile loading-unloading test was carried out for Aluminum

alloy 3003. Table 3 specifies the dimensions of the Dog-:Bone shape

 

 

 

 

 

 

specimen.

G - Gage Length 2.000 :0.005 in

W —- Width 0.500 i0.010 in

T - Thickness Thickness of Material
R - Radius of fillet 1/2 in

L - Overall length 8 in

A - Length of reduced section 2 1/4 in

 

 

 

 

61

$0.

 

 

 

 

B - Length of grip section 2 in

 

 

C - Width of grip section 3/4 in

 

 

Table 3 : Tensile specimen dimensions in the English units

The unloading of the specimen was done at the strain

5% and the following stress-strain curve was obtained:

increments of

 

Tensile Data

 

 

 

—_————

 

 

Tomlin Stress (Wu)

 

 

 

 

 

 

 

—— —————

-—-—__q

———————

———————

 

 

Tensile Strain ($6)

 

 

Figure 11 : Stress-strain curve for the loading-unloading test

The Young's modulus values are calculated using the linear portion of
the loading part after each unloading. The reason for taking the

loading portion is that we want to follow the same procedure as we

used for measuring the initial Young’s modulus.

parameter is calculated based on the strain equivalence principle by

using the following formula:

62

The damage

 

 

Where E is for the damaged material and measured at each unloading

increment. Eis for the very first loading step and its variation from the

data indicated in the handbooks shows an initial damage in the

specimen. The following data is obtained:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Unloading step 1 2 3 4 5 6 Efor the input material
E(Gpmor unloading data 47.771 42.85 39.72 35.29 34.29 33.1078 5143552201
E(Gpajfor loading data 50.225 47 43.79 40.8 37.15
Strain quLemaitrejforloading 0.0235 0.086 0.149 0.207 0.278 0.337 {—ch
Strain °/o 5.07 9.97 14.84 19.6 24.48 29.3
Table 4 : Young’s modulus data

Since 5‘ << 5", a values can be used to plot D versus a”.

,g 0.4 —

g 0.35 44— —+ , r 4 ,_MW 4 4 . «~—

,3; 0.3 41-——————— ‘/

: j‘.1

- 0.25 4 — rrrrr - w

8..§ l ./

E g 0.2 fi 4,3

a h

e '3 0.15 «i— /

O

s 0.1 +___ ,/

"3 / ‘1

E 0.05 7/

m o I I 7 r V l

5.07 9.97 14.84 19.6 24.48 29.3
Straln%

 

 

 

Figure 12 : Plot of Damage versus strain

This plot is used to calculate the material parameters. The variation of
damage shows a linear behavior as obtained by Lemaitre[14]. S is
determined by measuring the slope of the line using the following

formula:

63

7

0...

2E(l—D)

S:

 

2 dB
de”

We also use the Kuhn-Tucker condition for damage-coupled plasticity

which for perfectly plastic material is written as:

0'
:——0' :0
f l—D ’

Where a, is the yield value for the perfectly plastic material and is

usually taken to be equal to ultimate stress, on.

So:

0' 2
{—1
l-D

Thus, the equation for S can be written as:

 

 

2 2
s = a“ = 100 ..... 7.4891MPa
215111 2x51435x0.01298
d6”

law is obtained by intersecting the line with the horizontal axis:
em = 3.357
And finally, D1, is obtained by finding the value of the fitted line at the

failure strain:

0,, z 0.337

The mechanical properties of the specimen can also be found from the

tensile test. They are measures as:

64

0', == 4OMPa
0'“ == lOOMPa
E =- 51436MPa

4.2 A Simpler representation of the Lemaitre’s Damage

Evolution Equation

Let’s take the power law for hardening written as:

0' n
“=_= =K

If we use this relation in the definition of damage evolution equation,
we have:

_ 02R. ._K2(P)2"Rv.
-2ES(l—D)2 1" p

 

2E5

In the case of radial loading where the directions of principal stresses

remain unchanged, Rvis constant throughout the loading and we can

integrate the above equation to obtain a relation between damage

parameter and equivalent plastic strain:

 

_ ” K2(p)2"R. ._ K212. P 2,. ._K2Rv((P)2""-(Po)z"”)
0‘! 2ES p“ 2195 £07) ‘0' 2ES(2n+l)

PD

p = p, corresponds to D = Dcwhich is evaluated as:

_ K215 ((pR)2n+1 _(pD)2n+1)
” — 2ES(2n+1)

 

65

Since nis very small, (2n+1) is of order one and the above equation

can be written as:

 

= K2R.(pk -po)
C 2155
Or:
D, _ KZK

 

(PR—I’D) 2E5

For the one-dimensional tensile case, R, =1:

D,, _ K2

8R -£pD _ 2ES
2

If we replace 21:35 in the damage evolution law by its relation from the

 

above equation, we have:

D=[——D“ )MP)“ 1b

8R _€pD

This equation will be used in the next chapter for the coupled analysis

of damage evolution. D.” 8,, and 5,0 values for some materials can be

found in [14].

66

CHAPTER 5: Numerical implementations

In previous chapters, the isotropic formulations of damage coupled
constitutive equations were developed and the corresponding damage
material parameters were introduced. These results can be
implemented as a material model or in conjunction with the available
material models as a failure criterion. In this chapter, the constitutive
and rate equations are utilized in two different ways, uncoupled and
coupled, to analyze two simple problems under uniaxial and biaxial

loadings.

5.1 Uncoupled Analysis of Crack Initiation

Damage analysis of a model subjected to a given history of loading
consists in the calculation of the evolution of the damage as a function
of the time at the critical point(s) and the critical time at which the
damage reaches its critical value corresponding to a crack initiation.
An assumption which simplifies the analysis consists in neglecting the
coupling between the damage and the strains. In a first step the stress
and strain field histories are- calculated by a commercial finite element

software. Let us suppose that the results for point M are:

£(M,t)
0'(M ,t)

67

The second step consists in determining the point(s) where the
damage has its maximum value. There is no exact way to select the
critical point(s) except to calculate the evolution of the damage field in
each point. However, an "intelligent" look at the evolution of o as a
function of space and time will restrict the number of dangerous areas

where the damage must be calculated.
5.1.1 Integration of the Kinetic Damage Law

At the critical point, the structural calculation gives:

£,j(M*,t)
£‘r(M*,t)
£".;(M*,t)
0,}. (M*,t)

It is easy to deduce:
- The accumulated plastic strain rate, (It’s automatically calculated in

some finite element commercial software)
. 2 . . K
P =[§€pij€p0)

- The strain energy density release rare,

2
0' eq
Y =____
"’ 213(1—D)2 R”

From:

3 )2
0' :[EUDUUDU]

“I

68

 

R, =.:_(1+v)+3(1—2v)["”]

«I

And to write the general kinetic damage law as:

= Y0.D) . H _ 0.402190) .

1') _
S p (lip-€120) 2ES(l—D)2p (92-920)

 

The initial condition for the damage evolution is the end of the period
of accumulation of dislocations and the initiation of the first mesocrack

[15] corresponding to the value p0 of the accumulated plastic strain.

The initial conditions are:

The damage evolution is given by the integral:

D t
_ 2 LL 2 -
g (1 D )dD 2ES i 0,,(1) Rv(t)p(t)dt

.1._ _ 3 =;' 2 -
3( (1 D) +1) 255 {0.0) R.(t)p(t)dt

. 1%
D =1—[l——3—ja (t)2Rv(t)p(t)dt
2E5 ‘4 J
L to
The critical time t, at which a crack is initiated is reached when the

damage itself reaches its critical value given by the rupture criterion:

D,=D,,s1

If the material is highly ductile:

69

D=l

(‘

And t, can be calculated by the expression:

'R %

3 .
1:1— 1-3-Em a,,(z)2 Rv(t)p(t)dt

Or:
IR 2
j a,,(z)2R,(r)p(r)dt =-—ES
to 3
Where 1, is given by the elastoplastic analysis: to =t(p = pD)

5.1.2 Uncoupled analysis gives more conservative results

Uncoupled analysis is not an accurate analysis for determination of
the failure time. However, it gives more conservative results meaning
that the design is on the safe side but not optimized as opposed to the
strain-damage coupled calculation. In order to demonstrate this
analytically, we assume the material to be elastic-perfectly plastic and

the loading to be proportional:

= 0', = const

0'“,
R, =1
The last expression becomes:

’ . _2 ES
“(Mt—3&0,

to S

 

70

The strain history is given and this allows us to take [)(I) as a given

function, which is particularly simple to calculate if the elastic strain is

neglected:

. 2.. }5

Then t, may be determined.

For the case of damage-coupled analysis of the same material loaded
under the same conditions:

- For p< pD no damage occurs; the same calculation gives the same
result fort0 ,

t(p=PD)=to
- The material is perfectly plastic with the same threshold 0, which

allows us to write the coupled plasticity criterion as:

0'“,

—— = 0', = const.
(1 - D)

- The loading is proportional:

R, = const

- The same strain history is imposed and the elastic strain is again

neglected, thus the function [7(t) is the same.
The critical time for crack initiation in the coupled case r, is calculated

from the same kinetic damage law:

71

. Y, , 0 (0212.0). 0212.0).
D: ('D’ H = e" H = S H
s p (Sp-920) 2ES(l—D)2 p (92-900) 2ES p (Ep-Epo)

 

 

Where (1—D)2disappears due to the coupling in the plasticity criterion.

The integration is obvious:
402. ‘R .
l)=-———- tcfl
2155 JP”
‘0
Taking again D, =1as the critical condition:
D=D, =l—->t=t'R

"R . 2155
I p(t)dt = RN,

to S

 

Comparison with the uncoupled case shows that:

t' 'R

f p(:)d:=3j‘ pmdi

to to
Which implies that:

tR <t'R

In the particular case of loading where p =c0nst.

 

 

 

 

Uncoupled case Coupled case
r -t +———2ES t' -t + 2E5
R ° 3Rvazsp R ° K030

 

 

Table 5 : Comparison of the time to failure between uncoupled and coupled analysis
for monotonic loading

To summarize, the uncoupled calculation is always conservative. It is a

lower bound, but usually far from the coupled solution. In other terms,

72

 

component design may benefit from coupled calculations, which can
prove the enhanced safety of components or indicate how light-weight

economical components could be built.

5.2 Coupled analysis:

In this section, two algorithms are developed for the strain coupled
damage analysis and implemented in parallel with the finite element
software. In each algorithm, a simplifying assumption or condition is
taken into account to avoid excessive complexity in the formulation of
the algorithm. In this work, first the strain components were obtained
from Ls-Dyna and the first algorithm is utilized in Matlab and then the
second algorithm was written in Fortran and coupled as a material

model with Ls-Dyna.
5.2.1 Strain coupled algorithm with the assumption of perfect plasticity

By assuming the material to be perfectly plastic and by neglecting
the microcrack closure effect. The following equations should be

solved :

 

73

 

 

5'”,- =0 Otherwise.

2

pr>pD,then D=ZESR”p and if p<pD,then D=0.

 

Crack initiation if D, = D, is reached.

These equations may be used for piecewise perfect plasticity by

considering several values of the plastic threshold 0, as the loading or

the time like parameter vary. Then, the material parameters must be
considered as follows:

-Eand v for elasticity;

- f, a, and a, for plasticity

-a,=a,for perfect plasticity. Since, in reality, no material has the
perfect plasticity behavior, the choice of a, is somehow subjective.

However, choosing or, is a conservative choice which is taken in this

work.

- S,£,,D,D,C for damage

The input of the calculation is the time history of the strain

components 8,.(t) which may come from the result of a finite element

structural calculation. The outputs are:

74

- The time like parameter t,

- The damage evolution D(r), the last point corresponding to D=D,,

that is, crack initiation,

- The accumulated plastic strain evolution p(t) ;

- The evolution of the von-Mises equivalent stress 0,,(1‘),

The numerical procedure is a strain-driven incremental time like
algorithm using an elastic predictor and a plastic corrector. The
hypothesis of perfect plasticity allows us to explicitly formulate this
plastic correction and helps us see the softening effect of damage. The
constitutive equations are written in an incremental form
corresponding to a fully implicit integration scheme having
unconditional stability; but the damage parameter is updated
explicitly.

After an elastic prediction obtained from the law of elasticity, the

plasticity criterion fso is checked. If f>0, the plasticity corrections

are obtained by Newton's iterative procedure applied to a system of
two equations deduced from the constitutive equations.
At each time increment, the values of the stress and the other state

variables at the beginning of the increment (t,)are known and they
have to be updated at the end of the incrementum). In the following
formulation, tensorial notation is used for convenience.

We first assume that the increment is entirely elastic. So, we have:

75

62,“, = (l — Dom )0,” + Atr(A8)I + 2,u(A8)
Where xland ,uare the Lame coefficients and I is the second order

identity tensor.

 

 

)1: DE
(l+v)(l—20)
E
”—0— 2(l+v)

All the other plastic variables are equal to their values at(t,.). if this
elastic predictor satisfies the yield condition fso, then the elastic

assumption is correct and the computation for this time increment

ends. If f>0, this elastic state is corrected as explained in the

following in order to find the plastic solution.

The solution at (rm) has to satisfy the following equations:

(0...). __

 

= 0' =0
f l_ Dn+l 5
am: 0"“ +Azr(Ae)1+2p(Ae)-2p(AeP)
(1—D,,,,)
A8” = NAp
Y
AD=-—-
5 AP
3 0'” ,1 .
Where A(x)=x,,+, —x,andN =——-—"—. If we replace A8” by its

2 (UM-l )eq

expression in the second equation, we see that the problem is reduced

to the first two equations for the two unknowns, 6",, andAp.

76

We need to find 6,,,and prhich satisfy:
f :(6-n+l)eq -Us :0

h = a — 0" -4tr(A8)I — 20015) + 2,uNAp = 0

(1-0.)

This nonlinear system is solved iteratively by Newton's method. For

 

each iteration (s), we have:
f(5'n+ls TCa) = (6M1: +C6‘)eq ‘Us : 0
mam-i +C,-,,Ap+Cp) =0
So:
f +%3Ca =0

ah ah
h+—:C +—-C =0
35' 3 3p p

Where f, hand their partial derivatives are taken at (t,,,)and at the
iteration (s). (The starting point is the elastic predictor)

In other words, the integration is done explicitly over the time and

implicitly at each time over fiandp. The corrections Caand C,are

defined by:

C6 = 6-3-14 -6-3

Cp = ps+l __ ps
The starting iteration (s=0)corresponds to the elastic predictor. The

set of equations for fand h may be rewritten as:

77

Because:

And:

ah em
__= H 2 __
aa [ + ”A‘aa]

 

Where II is the fourth order identity tensor and: N =%

 

~ (3 ~D ~D ]% E(l—D)(3 e D e D]%
0- 2 ..

= —81'j £1]
1+0 2

 

M 3 1
—= — H——I®I —N®N
aa "' [2( 3 ) ]

0'“,

N :N=% and N :a—N=0 because (II—-:-I®I)is a deviatoric projector

35
because if you take its inner product with a 2nd-rank deviatoric tensor

N , the results is still N :

(fl—%I®I):N=N

So:

78

N: a——]§:=N: 0'3[21(n—11®1)_ N®N]=:l—[§N—(N®N):N]
0‘ 3 a 2

 

a 0,,
#l N (N NW] 3 [5-9.le
0,, 2 2
If we multiply the second equation of the

system,h+[n+2pApg—I;]:Ca+2#NCP=0, by N, the system gives

explicitly the correction for p :

f—N:h
3!!

C:

P

 

~

Now, the correction for a should also be found explicitly. First we

prove that C, is a deviatoric tensor. From the second equation of the

system and using the expression ofg—IY, one obtains:
0'

all=~;[-:—(H—%I®I)—N®N]

80 0,,

(H—%I®1)is a deviatoric projector because if you take its inner

product with a 2nd-rank tensor (a) you find the deviator of (a):

(H—%I®I):a=ab

(N ®N):a = (N:a)N this is collinear with N and since N is deviatoric, the

result is also deviatoric:

8N
— :=a deviatoric

85'

Therefore taking the trace on the second equation of the system gives:

79

3N
h+|:H+2;1Ap5-g]:C5 +ZpNCp =0
8N
tr(h) +tr(C&) + ”OWE; : C5) +tr(2,uNCp) = O

tr(2,uApiIZ:C5)=O because g—N;:C& is deviatoric.
0' 0'

tr(2,uNC,)=0 because N is deviatoric. So:
tr(h) + tr(C6) = 0
Using the definition of 6",, , and noting that:
6",, = 6', +ltr(A£)I + 20(A£)— 2,u(A£" )

tr(£pn) = O
tr(A£p) =0
tr(I)=3

We have:
”(6,,1) = (321 + 2,u)tr(£,,,,)
Therefore using the definition ofh:
h = 5i... - 6,, —/ltr(A8)I — 2,u(A£)+ szAp
And noting that tr(N) =0, we find that:
tr(C,) = —tr(h) = 0
These relations being true for any iteration (s), then tr(C,,)=Oand the

proof is achieved. Using this last result together with the expressions

0f %and C, yields for the second equation:

80

f+NzC3=0

h+[n+2mp%1§]zc,+2mc,=0
0'

 

~

31;: 1 [3(n—11®1)—N®N]
30' a 2 3

ea

=f—N:h

Where:

A=[H+2,w_\p%]
1

0'“,

 

+1.2... [;(N-.;.I®N-N®N]]

In calculatingAzCa, the term ’ [%I®I]:Ca =(%I:Ca)1=[%tr(Ca)]I =0

because C5. is deviatoric.

Therefore :

l

4:0, =[(1+;—”Ap)r1— 2” ApNeNJ:C,

 

~

“I “1

Although:

A¢[(l+;—#Ap)II—%£ApN®N]

eq eq

81

So if we cal|B=[(1+gfi-Ap)H--E—”—ApN®N], the tensor Bis invertible if
0'

eq eq

and only if 5,, =0, and in this case:

B" = ——31-—[n +g—‘E-ApN ®N:l
l+—~E—Ap ‘4

‘4

Using this expression gives explicitly the correction fora:

Ca =—-3-(f—N2h)N———3}-—[h+ 3'”
1.2—$513,) «1
eq

 

Ap(N : MN]

Once pand (fare found, spand D are calculated from their discretized

constitutive equations and the stress components are given by:
0.1] = (1— D)5:ij

The following flow chart summarizes the above-mentioned procedure.

82

 

( fl
[ Lu = Variable at t,

[ Lew = Variable at tn+1

A8 = 5",, — 80,,

All the variables are known at t- J

T
s = O & Ap = 0 (Elastic Predictor)
6;” = (1 — D,” )0,“ + Xtr(A£)I + 2,u(A£)

 

 

 

 

 

 

 

 

 

 

 

0’
h = 6:0, — ¢ — lumen - 2p(A£) + ZpNAp
(l—Dold)
_ ~3 D
CH N'h&1v=3 0m
34 2&3...)
k “1
C,=—3(f-N:h)N— 31 h+~32—'uAp(N:h)N
3 1+ I“ Ap (anewlq
(2:4)
«I
6:1:=a;,,+c,&Ap=Ap+C,&s=s+1

 

 

 

 

No

 

 

 

 

83

This algorithm was then written in Matlab and Fortran and paralleled

with Ls-Dyna as a user defined material to give stress distribution.

5.2.2 Strain coupled algorithm with the assumption of using the

damage value at the beginning of the increment

This time, we write the yield function in terms of the Cauchy stress
instead of the effective stress and we take into account the damage
effects into the yield stress. Thus, the von-Mises yield function is

written as:

f(aD)=%aD :a" -(1—D)20'y2 so
Radius of this yield surface is:

2
R = \[g(l—D)ay

And the normal unit vector of the yield surface is:

N=i=_2”_=9_"_ J24;

[0”] an :02 R 2 (l—D)a'y

The other equations used in this section are:
Total strain rate in its additive form:

é=é‘+£‘p

Flow rule:

Equivalent plastic strain rate:

84

2 2
'= —é‘”:é"= —°
p \(3 (3”

d- : (1—D)(/itrace(é‘)l+2#ée)

Stress rate:

Equivalent von-Mises stress:

0,, - 30'” 20'”
2
Isotropic hardening law:
0'y = K (80 + p)"

Lemaitre’s damage law:

D:[g 131; ][%(1+v)+3(1-2v)[g-”—]2](€o+P)2”P

The Computational procedure is as follows:

 

First, we assume that the increment is totally elastic and we calculate
the elastic predictor as the trial stress for the new increment:

0133’ = 00,, +(l— Dd,)(/1trace(A8)I + 2,uA£)

Then, if the elastic predictor, 0“” , is within the current yield surface,

MW

the elastic predictor is the actual value of stress for the current
increment:

_ rial
03W — 0”

new

85

If not, we need to subtract a plastic correction written as following:
anew : 0:21:21 "ZMEP (1_Dold) : 0:12:1_21‘1A7N(1—Dold)

Since the increment is no longer totally elastic, we need to update the

value of total equivalent plastic strain, p :

2
pnew : paid +J;A}’

And the new yield stress considering the plastic hardening is calculated

as:

n- 2
Gym = 0,0,, +hAp = 0,0,, +nK(.€0 +p) ngAy

prZé‘pD, the damage parameter has to be updated using the

Lemaitre's equation:

2
D 2 0 2n 2
D =D+—4——1+ +3l—2 ” 6+ —A
new old [£R_EPD][3( V) ( V)[0' J]( 0 p) J; 7

eq

 

Then, A7 is determined considering constant D during one step which

is justified in the explicit calculation because of very small increments.

At the end of new step:

2
UDnew : anN = £(1 - Dold )aynewN

m‘a 2
0.0M): —2/‘IA},(l—Dold)N=J—3-(l—Dold)[ayold +J:hA}/]N

00:11 = {E[0ydd +£hAy]+2pAy}(l—D,M)N = {€0,014 +(%h+2,u)Ay}(l—D,,d)N

86

DJN

 

JOBS: 30'0er : (l-DOM){‘/%0'yold +[%h+2fl)A}’}

W: 1.1 {1002140121 E(l—D.,.)a...}
2fl[1+§;](l_Dold)

 

 

The following flow chart summarizes the above-mentioned

computational procedure.

87

 

/

[ L, .1 = Variable at t,

[ l...

A.s=.€,,+1 —6

n

= Variable at tn+1

 

 

All the variables are known at t-

l

[ om"=0,“+(1-Ddd)(/1trace(A£)I+2/1A£) ]

\

 

 

 

   

f(o’,’,‘f;’)so =03

 

 

 

 

No

 

l

 

 

A?

1 .
= JUDmal

 

Dm'al
: 0'

"(W

 

new

(E (1 D...)a W}

241.2%)(1—0,»

2
pnew =pold +J:Ay
3
Bold =+R'%c[£

DE]( (1—+v)+3(1 2v)[

abm’al

N: new
Dm'al

I: ...

 

 

2
0'” 2n 2
-A
0:]:|(80+P) J; 7

—2:uA}/N(1_Dold)

-€p

_0.m'al_

new

Andom,

 

 

 

88

 

5.3 Numerical Results

In what follows, the results of numerical analysis are presented and
compared for three analysis cases of without damage, uncoupled
damage, and coupled damage material model. The case problems are
uniaxial and biaxial tensions and the material properties are those of

copper 99.9% obtained from [14]:

 

 

 

 

 

 

 

 

 

 

 

 

Material Copper 99.9%
Young’s modulus 98990 MPa
Poison’s ratio 0.34
Yield stress 90 MPa
Ultimate stress 500 MPa
K (Hardening coefficient) 491.3 MPa
n (Hardening exponent) 0.2459
Initial damage threshold 0.35
Final damage threshold 1.07
Critical damage 0.85
Damage strength 1.0696

 

 

Table 6 : Material properties for numerical calculations

The geometry of the models is shown in the following. The thicknesses
are 1mm for the uniaxial case and 10mm for the biaxial case. In each

case, some kind of defect or stress concentrator is needed to induce

89

 

localization. This is done in the uniaxial case by decreasing the width
at the center by 2mm and in the biaxial case by considering a hole in
the center. Due to symmetry, just a quarter of the model is considered
for the biaxial case and the symmetry boundary conditions are applied

to the edges adjacent to the hole.

‘ 1
50.0 mm

L
, :

200.0 mm 50,0 mm

“—20.0 mm

 

 

 

 

 

12.5

 

 

 

 

Figure 13 : The front view of the uniaxial test model

90

: 100 mm _

 

 

80 mm

 

 

 

 

 

Figure 14 : The front view of the biaxial test model

V
A

Figure 15 : The isoperimetric view of the biaxial test model

Since the code is developed for the general case, three dimensional

brick elements are used to mesh the model avoiding any plain strain or

91

 

plain stress assumption. The preprocessing is carried out in

Hypermesh and the model is then solved by Ls-Dyna.

92

Side View

 

 

    
 

    

.‘J-f HE
,/""‘ F 8‘ re Eie r \\
Velocity BC's f... 0r, mg men \,
Presa‘bed on both ,. Ana y31s ‘1
ends )/ ‘
.7 I
)i N
\ /
\ |_ Z .
\\ V ,./
\ \ ’//
\-‘7\._‘ M/
K, ___,_._-v-

Figure 16 : Mesh configuration and location of the critical element for single element
analysis in the uniaxial test model

93

fl! Velocity 805

 

 

  
 
 
 
       
    
        
 
 

Isoperimetnc View of /
Thickness

 

I My;
Symmetry 8C5

  

 

 

For Single Element

  
 
     
 

 

    

 
  

\ III

  

   
      

 
 

 

 

 

 

 

 

Analyses

 

 

 

\\.

Figure 17 : Mesh configuration and location of the critical element for single element

analysis in the biaxial test model

First, the algorithm mentioned in section 5.2.1 was written in Matlab

and the strain components for the element number 700 of the plate

specimen which is shown in Figure 17 were obtained from Ls-Dyna

with the assumption of perfect plasticity. The following plots show the

comparison of von—Mises stress, effective plastic strain, and damage

for three analysis cases of without damage, uncoupled damage, and

coupled damage.

94

Damage softening effects

1 I I I l

 

.0
O)
1
l

Th§CConventional FE without the damage effects_

.°
01
1

.°
.5
I
1

 

Equivalent van-Mises stress (GPa)

 

 

 

0'3 P The damage softening effects A
0.2 _ ~
0.1 f —
0 l l l l l
0 1 2 3 4 5 6
Time (ms)

Figure 18 : Comparison of von-Mises stress for without damage and coupled damage
analyses for element number 700

Damage softening effects are noticeable and indicate a lower load

bearing capability than what is calculated with a non-damage analysis.

95

Comparison between the evolution of Damage and Equivalent plastic strain

7 l I I T l

 

7 {Total plastic strain evolution withoUt damage

.5
N
1

Total plastic strain evolution for the Coupled analysisl

  

—l
l

Damage evolution for the Coupled analysis

.0
m
1

I
l.
,l
.l
-I
'l
I
I
I T .
'1
l
I
I
l
I
.l
I
l
l
I

.o
m
1

 

Equivalent plastic strain and Damage evolution
9
g
1

.0
N
1

 

. 1 1 1 1
0 1 2 3 4 5 6 7
Time (ms)

 

mL-t

 

Figure 19 : Damage evolution and comparison of effective (total) plastic strain for
analysis cases of without damage and coupled damage

Effective or total plastic strain is one of the most common failure
criteria and is commonly used in element deletion or node release
processes in the finite element codes. The above mentioned plot shows
that the damage parameter not only takes into account the continuous
degradation of the material properties but also the failure prediction
based on the critical value of damage is in good agreement with the
critical effective plastic strain criterion. That is, the damage parameter
reaches its critical value, 0.85, when the effective plastic strain is

approximately equal to its critical value, 1.07.

96

The following plot shows conservative prediction of failure from
uncoupled analysis, which was analytically derived in section 5.1.2.
The rate of damage evolution increases rapidly as damage increases

since there is a (1—D)2term in the denominator of damage rate

equafion.

The evolution of Damage and Equivalent plastic strain for Uncoupled analysis

1 1 1 1 1

 

.a
N
l

- --<

. Damageevolution for the Uncofipled, analyse

—L
l

otal plastic strain evolution . . ~ . —.

Q
G
T

Equivalent plastic strain and Damage evolution
9 o
h C)

C_
M
F

 

 

l— ~
.—

 

 

#
3

Hmemm)
Figure 20 : Conservative prediction of failure from uncoupled analysis

Figure 18, Figure 19, and Figure 20 are results for the biaxial model
and for a single element (The element is shown in Figure 17). In what
follows, the algorithm mentioned in section 5.2.2 is written in Fortran

and implemented as a user defined material model with Ls-Dyna. The

97

showing the continuous degradation of the material. Another point of
interest is the comparison between the prediction of failure by the
critical value of damage and by the critical value of effective plastic
strain. As shown in Figure 21 and Figure 22, damage occurs when the

effective plastic strain reaches the initial damage threshold:

98

 

1: ii
i

\ ' -,l:r-\;«Q i
' 4“,."

.h;

       

Figure 21 : Initiation of Damage for the uniaxial model. (The bottom figure shows
the effective plastic strain)

99

 

Figure 22 : The Initiation of Damage for the biaxial model. (The bottom figure shows
the effective plastic strain)

100

In the uniaxial model, both damage and effective plastic strain start to
grow from the edges where the three dimensionality thus the

triaxiality ratio are maximum. This is shown in Figure 23.

 

101

 

“I ...¢

I
I as.-
.i-In: B“

    

_
,
.!

inn-5, win-l!
In: Fifi—Eul-
--- ____=_------
-----—_=__-_------
---= -__—=_------
“Essen
assess...
In: :5: E..-

.
- = 1... =5-
rP— _._ E 2..-E!

Figure 23 : Damage and effective plastic strain start to grow from the edges (The top
picture shows the effective plastic strain)

102

Then, damage reaches its critical value when the effective plastic
strain is approximately equal to its critical value both at the same
location. This indicates that the critical value of damage can ,
satisfactorily predict failure. In Figure 24 and Figure 25, a comparison
is made between contours of damage and effective plastic strain when

damage reaches its critical value.

103

 

Figure 24 : Localization of damage and its comparison with effective plastic strain for
the uniaxial model (The top picture shows the effective plastic strain)

104

 

Figure 25 : Localization of damage and its comparison with effective plastic strain for
the biaxial model (The top picture shows the effective plastic strain)

105

 

 

 

 

 

 

2.6 2.8
Time (1115)

Figure 26 : Evolution of damage (black) and effective plastic strain (blue) of the
element shown in Figure 16 for the uniaxial model

1.2 i i i i

 

 

 

 

 

Time (ms)

Figure 27 : Evolution of damage (black) and effective plastic strain (blue) of the
element shown in Figure 17 for the biaxial model

106

And finally, the softening effects of damage can be seen in the
following figures and plots. The von-Mises equivalent stress is
compared between the damage model results and results from a

model with no damage consideration.

 

107

 

Figure 28 : Distribution of the von-Mises stress for the uniaxial model. (The bottom
figure shows the damage model result and the top figure shows the non-damage
model result)

108

 

.°
.15

 

.9
h.)

 

 

 

 

\

1-....-_--.._..-- 1

 

Effective Stress (v-m)

 

 

 

 

 

 

 

l L i l i L l l
T l
1.5 2 2.5 3
Time

9
.°
01
c—I

Figure 30 : Damage model result for von-Mises stress in element shown in Figure 16
for the uniaxial model

 

.9
c»

 

P
01

 

O
.5
fi

\

 

 

\

Effective Stress (v-m)
I

 

 

 

 

 

 

 

 

 

P
o —b
\
l—

 

0 0.5 1 1.5 2 2.5 3 3.5
Time

Figure 31 : Non-Damage model result for von-Mises stress in element shown in
Figure 16 for the uniaxial model

110

0.4 .

 

 

 

O
(.9

 

 

 

 

 

 

 

 

5.
3.. ’ i
g i
a) 02 .
a / i
o .
.3 i
g
a: 0.1 i
“J
l
l

c l L 1 i l l I.

0 1 2 3 4 5

Time

Figure 32 : Damage model result for von-Mises stress in element shown in Figure 17
for the biaxial model

0.5

 

O
A
1

 

" /

 

\

 

.0
nub

Effective Stress (v—m)
O
N

 

 

 

 

 

 

 

 

 

Time

Figure 33 : Non-Damage model result for von-Mises stress in element shown in
Figure 17 for the biaxial model

111

References

[1] L.M. Kachanov, Introduction to continuum damage
mechanics (Dordrecht ; Boston : M. Nijhoff, 1986).

[2] D. Krajcinovic, G.U. Fonseka, The continuous damage
theory of brittle materials-parts I & II, Journal of Applied
Mechanics, 48 (1981) 809-824.

[3] D. Krajcinovic, Continuous damage mechanics revisited:
basic concepts and definitions Journal of Applied Mechanics, 52
(1985) 829-834.

[4] F.A Leckie, E.T. Onat, Tensorial nature of damage
measuring internal variables, Proc. IUTAM Symp. on Physical
Nonlinearities in Structures, (1987) 547-558.

[5] ET. Onat, F.A. Leckie, Representation of mechanical
behavior in the presence of changing internal structure, Journal
of Applied Mechanics, 55 (1988) 1-10.

[6] J.L. Chaboche, Anisotropic creep damage in the framework
of continuum damage mechanics, Nuclear Engineering and
Design, 79 (1984) 309-319.

[7] J.L. Chaboche, Thermodynamically founded CDM models
for creep and other conditions, CISM Courses, "Creep and
Damage in Materials and Structures", H. Altenbach and 3.3.
Skrzypek eds., Springer, 209-283, 1999.

[8] S. Murakami, Mechanical modeling of material damage,
Journal of Applied Mechanics, 55 (1988) 280-286.

[9] CL. Chow, J. Wang, An anisotropic theory of elasticity for
continuum damage mechanics, International Journal of Fracture,
33 (1987) 3-16.

[10] CL. Chow, J. Wang, An anisotropic theory of continuum
damage mechanics for ductile fracture, Engineering Fracture
Mechanics, 27 (1987) 547-558.

[11] CL. Chow, T.J. Lu, On evolution laws of anisotropic

112

damage, Engineering Fracture Mechanics, 34 (1989) 679-701.

[12] CL. Chow, J. Wang, Ductile fracture characterization with
an anisotropic continuum damage theory, Engineering Fracture
Mechanics, 30 (1988) 547-563.

[13] CL. Chow, J. Wang, A finite element analysis of continuum
damage mechanics for ductile fracture, International Journal of
Fracture, 38 (1988) 83-102.

[14] J. Lemaitre, A continuous damage mechanics model for
ductile fracture. Journal of Engineering Materials and
Technology, 107 (1985) 83-89.

[15] J. Lemaitre, A Course on Damage Mechanics, Second
Edition (Springer: New York, 1996)

[16] J. Lemaitre, J.L. Chaboche, Mechanics of Solid Materials
(Cambridge University Press, 1990).

[17] J. Lemaitre, I. Doghri, Damage 90: a post processor for
crack initiation, Computer methods in applied mechanics and
engineering, 115 (1994) 197-232

[18] J.Lemaitre, J. Dufailly, Damage Measurements,
Engineering Fracture Mechanics, 28 (1987) 643-661

[19] J. Lemaitre, R. Desmorat, M. Sauzay, Anisotropic damage
law of evolution, European Journal of Mechanics - A/Solids , 19
(2000) 187-208

[20] 1.]. Skryzpek, Material damage models for creep failure
analysis and design of structures

[21] G.T. Mase, G.E. Mase, Continuum Mechanics for Engineers,
Second Edition (CRC Press, 1999)

113

1liliijiiiiiijji1111111111