Computer methods
in applied
mechanics and
engineering
EI_SEVIER
Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
Damage 90: a post processor for crack initiation
Jean
Lemaitre*'t
Prof. Universit~ Paris 6, LMT Cachan, 61 Avenue du President Wilson, 94230 Cachan, France
Issam Doghri
Department of Mechanical and Environmental Engineering, University of California, Santa Barbara, CA 93106, USA
Received 10 December 1992
Abstract
A post processor is fully described which allows the calculation of the crack initiation conditions from the history of strain
components taken as the output of a finite element calculation. It is based upon damage mechanics using coupled strain damage
constitutive equations for linear isotropic elasticity, perfect plasticity and a unified kinetic law of damage evolution. The
localization of damage allows this coupling to be considered only for the damaging point for which the input strain history is taken
from a classical structure calculation in elasticity or elastoplasticity. The listing of the code, a 'friendly' code, with less than 600
FORTRAN instructions is given and some examples show its ability to model ductile failure in one or multi dimensions, brittle
failure, low and high cycle fatigue with the non-linear accumulation, and multi-axial fatigue.
1. Introduction
D u r i n g t h e ' p a s t decades, the p r o b l e m of m a c r o c r a c k growth up to failure of structures has received
m u c h attention t h r o u g h the t r e m e n d o u s d e v e l o p m e n t of fracture mechanics. Crack initiation, which
requires a k n o w l e d g e of what happens before a mesocrack breaks the representative v o l u m e e l e m e n t
( R V E ) , c a n n o t be treated by classical fracture mechanics dealing with pre-existing m e s o or m a c r o
cracks. T h e usual engineering m e t h o d of predicting the crack initiation consists of using s o m e
p h e n o m e n o l o g i c a l or empirical criteria or relations to define the conditions of local rupture of which the
following are examples:
- a critical value of the m a x i m u m principal stress;
- the von Mises equivalent stress or better the d a m a g e equivalent stress equal to some critical value for
brittle fracture [1];
- t h e von Mises equivalent plastic strain equal to some critical value for ductile failure or better the
micro-void growth models up to instability [2, 3];
- t h e M a n s o n - C o f f i n law for low cycle fatigue [4, 5]
- t h e W o e h l e r - P a l m g r e e n Miner rule for high cycle fatigue [6];
- for creep, a relation b e t w e e n the stress and the time to rupture c o r r e s p o n d i n g to an infinite strain flow
in a creep constitutive e q u a t i o n [7] or a stress criterion based on a linear c o m b i n a t i o n of the stress
tensor invariants [8].
T h e p o o r accuracy of such predictions in complex cases comes from the inability of those models to
* Corresponding author.
Visiting Professor at University of California, Santa Barbara.
0045-7825/94/$07.00 t~) 1994 Elsevier Science B.V. All rights reserved
J. Lemaitre. 1. Doghri I Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
198
take into account the history of loadings and for most of t h e m , the difficulty of being written in
three-dimensions.
T h e tool for modeling the evolution of the progressive deterioration of materials up to a m e s o c r a c k
initiation is the d a m a g e mechanics which was d e v e l o p e d after fracture mechanics but which is n o w in its
stage of maturity. A continuous d a m a g e variable is introduced as an internal variable. It is of tensorial
nature, but it reduces to a scalar if the d a m a g e is isotropic [9] which is the hypothesis considered here.
It is the m a x i m u m equivalent area density of microcracks or microcavities which lies in any cross plane
that can be defined in a R V E [10]:
D-
6Sp
6S ' 0 ~ < D ~ < I '
w h e r e 8S is the area of the most d a m a g e d cross section of the R V E and 8S o is the equivalent area of
the microcracks and microcavities in 6S.
T h e concept of effective stress ~q = ~rq/(1 - D ) [11] associated with the principle of strain equivalence [12] allows one to write the H e l m h o l t z free energy qJ in the f r a m e w o r k of the t h e r m o d y n a m i c s of
irreversible processes. T h e couplings or uncouplings b e t w e e n elasticity, plasticity and d a m a g e are
identified and m o d e l e d in accordance with the 'state kinetic coupling t h e o r y ' [13] which defines a logical
p r o c e d u r e to choose the analytical expressions of the state potential q, and the potential of dissipation F
from which are derived the kinetic constitutive equations. If small d e f o r m a t i o n s and linear isotropic
elasticity are a s s u m e d ,
E[
1
~ ~
v
q / = ~ p [-]--+~v e'ie'J + (1 + v)~l - 2v)
~.2]
ekk_ (1 - D ) + terms for plasticity,
e
w h e r e % is the elastic strain tensor, E and v are the Y o u n g ' s m o d u l u s and Poisson's ratio and p is the
mass density considered as a constant.
T h e law of elasticity coupled to d a m a g e is derived from
o-, =p 3eyj -E(1 - D ) L-TT-~v + (1 + v-~]C 2v)
or
l+v
v
6 0 is the K r o n e c k e r ' s o p e r a t o r
6q -- 1 if i=j, 6q = 0 if i # j . T h e associated variable to D is defined by
t2-
-Y=P
aD-
2 LI+v
+ (l+e)(1-2v)
"
Y is the strain energy density release rate [14] defining the p o w e r dissipated in the d a m a g i n g process as
d~D=YL)~>0,
dD
L ) = dt '
t being the time .
It can be written in a close form as
y -
°r~q R,,
2E( 1 - D ) 2 ,
D
w h e r e ~%q is the yon Mises stress ~rq = ,~, 3 o D # cr#D x)l / 2 , crqD is the stress deviator ¢z~j
= o-a - 6rH6q, O"H is the
hydrostatic stress crH = ~crkk, and R,, is the triaxiality function which d e p e n d s on the triaxiality ratio
cru/crq, i.e.
2
R,, = ~ ( 1 + v ) + 3 ( 1 -
2v)( °'n 12 .
\ O-~qI
J. Lemaitre, I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
199
This allows the definition of the damage equivalent stress ~r* as the one-dimensional stress which gives,
for the same value of the damage, the same elastic strain energy density w e as the three-dimensional
state [1]. Since Y = W e / ( 1 - D), then t r * = O-eqRl,/2.
These basic concepts of damage mechanics are used in Section 3 to derive the set of constitutive
equations.
2. Locally coupled analysis of damage in structures
In most applications, the damage is very localized in such a way that the damaged material occupies a
volume small in comparison to the macroscale of the structural component and even to the mesoscale of
the RVE. This is due to the high sensitivity of damage to stress concentrations at the macroscale and to
defects at the microscale. This allows us to consider that the effect of the damage on the sate of stress
and strain occurs only in very small damaged regions. In other words, the coupling between damage
and strains may be neglected everywhere in the structure except in the RVE(s) where the damage
develops. This is the principle of the locally coupled analysis [15, 16] where the procedure may be split
into the following two steps as shown in Fig. 1:
- a classical structure calculation in elasticity or elastoplasticity by the finite element method ( F E M ) or
any other method to obtain the fields of strain and stress;
- a local analysis at the critical point(s) only dealing with the elasto-plastic constitutive equations
coupled with the kinetic law of damage evolution, that is a set of differential equations.
This method is much simpler and saves a lot of computer time in comparison to the fully coupled
analysis which takes into account the coupling between damage and strain in the whole structure. The
fully coupled method must be used when the damage is not localized but diffused in a large region as in
some areas of ductile and creep failures [17].
The locally coupled analysis is most suitable for those cases of brittle or fatigue failures when the
structure remains elastic everywhere except in the critical point(s).
If the damage localizes, this is of course because stress concentrations may occur but also, because
some weakness in strength always exists at the microscale. Let us consider a micro-mechanics model of
a RVE at mesoscale made of a matrix containing a micro-element weaker by its yield stress, all other
material characteristics being the same in the inclusion and in the matrix (Fig. 2).
- The behavior of the matrix is elastic or elasto-plastic with the following characteristics of the material:
yield stress o-v, ultimate stress o-u and fatigue limit o-f (of < O-y).
- The b e h a v i o r o f the inclusion is elastic-perfectly plastic and can suffer damage. Its weakness is due to
a plastic threshold which is taken to be lower than o-y and equal to the fatigue limit crt if not known by
another consideration
O.p.
s ~ O ' y / . L Off.
POST PROCESSOR
El~to(pl~tic)~
Constitutive
Eq..~tions
Ehk~oplL~icityg~lplcd
I
todmgecomtitutiveJ
equaaom
I
CRACK
INITIATION
CONDITIONS
Dlmlng¢mec.hani~
Cootinuum mechanics
FEM Calculations
Fig. 1. Locally coupled analysis of crack initiation.
200
J. Lemaitre, I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
).
J
3t
~
oTyr=0-~
Fig. 2. Two scale representative volume element.
Below the fatigue limit ort no damage should occur. The fatigue limit of the inclusion is supposed to be
reduced in the same proportion:
orf
O'f
~
O'f-.
O" v
The second assumption which makes the calculation simple is the Lin-Taylor's strain compatibility
hypothesis [18] which states that the state of strain at microscale is equal to the state of strain imposed
at macroscale.
Then, there is no boundary value problem to be considered. Only the set of coupled constitutive
equations must be solved for the given history of the mesostrains at the critical point(s).
The determination of this critical point is the bridge between the FEM calculation and the post
processor. If the loading is proportional, that is constant principal directions of the stresses, this is the
point M* where at any time the damage equivalent stress is maximum,
O-[M. ) = Max(o-[M))---> M* .
If the loading is non-proportional, the damage equivalent stress may vary differently at different points
as functions of the time-like parameter defining the history. A small cr* for a long time may be more
damaging than a large o-* for a short time! There is no rigorous way to select the critical point(s) but an
'intelligent' look at the evolution of or* as a function of the space and the time may restrict the number
of the dangerous areas to a few points at which the calculation may be performed.
Another point of interest is the mesocrack initiation criterion. The result of the calculation is a crack
initiation at microscale that is a crack of the size of the inclusion; its area is 8A = d 2. This corresponds
to a value of the strain energy release rate at mesoscale G that may be calculated from the variation of
the strain energy W at constant stress,
G
~
lOW
2 OA '
_ _ - - - -
But - 8 W / 2 is also the energy dissipated in the inclusion by the damaging process at constant stress,
Dc
d3/
G(~. 4 ) _
o
Y dD
8A
Assuming a constant strain energy density release rate, Y - Yc yields
G~A) =
YcDcd.
Another interesting relation comes from the matching between damage mechanics and fracture
mechanics. Writing the equality between the energy dissipated by damage at microcrack initiation and
the energy dissipated to create the same crack through a brittle fracture mechanics process,
201
J. Lernaitre. I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
ycDcd 3 = God 2 ,
where Gc is the toughness of the material. It yields
YcD¢d = G~ .
Comparing the two relations for YcD~d shows that G(~A) = G~ is the instability of the meso RVE. Then,
the criterion of microcrack initiation is also the criterion for a brittle crack instability at mesoscale.
3. Constitutive equations
3.1. Formulation
The following hypotheses are made:
- small deformation;
- uncoupled partitioning of the total strain in elastic strain and plastic strain % = e,~ + e~P;
-linear isotropic elasticity;
- perfect plasticity.
The latter hypothesis is justified by the fact that in most cases damage develops when the
accumulated plastic strain is large enough to consider the strain hardening is saturated. Nevertheless,
the code described in Section 4.3 allows one to consider stepwise perfect plasticity in order to take into
consideration high values of strain hardening and the cyclic stress-strain curve for multilevel fatigue
processes.
The threshold of plasticity o-s is bounded by the fatigue limit and the ultimate stress
crf ~< Crs ~ ~ u .
The plastic constitutive equations are derived in the most classical way from the yield function in
which the kinetic coupling with damage obeys the strain equivalence principle applied on the plastic
yield function:
~
~
f=O'q-O-s=0,
Oreq
witho-~.- l - D "
The plastic strain rate derives from f through the normality rule of standard materials
[19],
D
t~P = Ocr;j
20req 1 - D '
if
0,
where A is the plastic multiplier. This shows that
~D
1 -
-
(2
\t/2
/~iP~ ip)
= [~'
where p is the accumulated plastic strain.
The kinetic law of damage evolution, valid for any kind of ductile, brittle, low or high cycle fatigue
damage derives from the potential of dissipation F in the framework of the state kinetic coupling theory
[13].
y2
F = f + 2S(1 - D) '
2 Rv
0req
2
o'sRv
Y - 2e(1 - D) 2 - 2-----E- '
where S is a damage strength, a constant characteristic of each material,
OF.
Y A
D=-~-yA-s
1-D
Y
or
D=-~p,
ifp>~pD
(mesocrack initiation if D = D c as demonstrated in Section 2).
J. Lemaitre. 1. Doghri / Comput. Methods Appl, Mech. Engrg. 115 (1994) 197-232
202
cc--;,
"
/
/
.--.------:-
"Si-o re,d' ~_r~cr~,,.,
-7--)--:
Fig. 3. Modelling a stress-strain curve.
PD iS a damage threshold corresponding to microcrack nucleation; it is related to the energy stored in
the material. Consider the material of the inclusion having a plastic threshold stress tr~ and a fatigue
limit o "r = o~/o:, considered as the 'physical yield stress' (Fig. 3).
The second principle of thermodynamics states that in the absence of any damage, the energy
dissipated in heat is o-~p. Then the stored energy density is
w =f%de
-o
p.
From the perfectly plastic constitutive equation,
w, = (os - o4')p.
Considering that the threshold PD is governed by a critical value of this energy identified from the
unidimensional reference tension test having a damage threshold epD, a plastic threshold equal to the
ultimate stress o-u and a fatigue limit of.
O u -- Of
(o~ - o~)pD = (Ou -- o f ) e p o
and
P o = epD o~ - cr--~f
with o- "f = o-~/o-v. If piecewise perfect plasticity of several thresholds oi~i is considered
w =2
-
and PD is defined by
PD
E
- Or" PD = ( o u -- °r)ePD
p-:0
The critical value of damage at mesocrack initiation D c corresponds to an instability [20].
Nevertheless, it can be related to a certain amount of energy dissipated in the damaging process
Dc
Y d D = constant at failure.
0
Assuming for that calculation a proportional loading for which R,, = const,
Dc
Dc
YdD=
0
~ d D =
2E
De.
0
Taking again the pure tension test as a reference which gives the critical value of the damage at
mesocrack initiation Die related to the stress to rupture o-R and the tensile decohesive ultimate stress o-u
together with R,, = 1 gives the rupture criterion,
o'eR,,
o'2
2E D~ = ~ - D i c
,
J. Lemaitre, 1. Doghri I Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
203
that is
2
o-u
D¢ = Dl~ o-~R, "
Collecting all the equations gives the set of constitutive equations to be solved by a numerical analysis:
gij
c
eij -{- e p ,
eu-
l +v %
E 1 -D
v o-k~DSU
E 1'
~D
"3 o-o .
4=
p ,
O,
£)=
if f = O,
o-~q
f=l-D
o'~,
iff<O,
~ E s R , ,,6'
if p
0,
ifp <PD ,
o-u - - o-f
~>PD ,
PD = epD
-
"
Mesocrack initiation occurs if D = Dtc(o-z,/o-~R, ) <- 1.
As the plastic strain rate ~P is governed by/5, and/~ is related to the damage rate D, in this model it
is the damage which governs the plastic strain.
Note also that the perfect plasticity hypothesis allows one to write R~ only as a function of the total
strain, since
_
o-H
E(1 - D ) ¢
1 - 2v eH '
el. | =
e H
as
g'H =
-'3(ekk + EPk) and eP, = O,
o-¢q = o-~(1 - D ) .
Then
Or H
E
e H
~q-
1-2v
o-,
3.2. Identification of the material parameters
In order to be able to perform numerical calculations, it is an important task to determine the
material coefficients in each case of application. This is always difficult because one is never sure that
the h a n d b o o k or the test results used correspond exactly to the material considered in the application.
Those coefficients are: E and v for elasticity; o-f, o-y, o-, and the choice of ~ for plasticity; S, epD, D~¢
for damage.
T h e y can be easily obtained from a tensile stress-strain curve with unloadings or from a very low
fatigue test at constant strain amplitude (N R ~< 10) as those of Fig. 4.
A classical experimental procedure with strains measured by strain gauges permits the determination
of Young's modulus E, Poisson's ratio u and also the damaged elasticity modulus E as defined on the
graphs of Fig. 4 and as a function of p:
- - p = ep in Fig. 4(a);
- p = 2N Aep in Fig. 4(b) since Ae~ = 2o-M/E(1 -- D ) = 2o-,/E = const, where N is the n u m b e r of cycles.
T h e yield stress % and the ultimate stress or, are also of a classical procedure. The fatigue limit o-f is the
stress for which the n u m b e r of cycles to failure N R in a high cycle fatigue test is very large:
N R ~ 106 to 107 cycles. Usually
The damaged p a r a m e t e r s are deduced from the measurement- of E as the law of elasticity allows us to
write [21]
204
J. Lemaitre. I. Doghri I Comput. Methods Appl. Mech. Engrg. I15 (1994) 197-232
'O--
Uo-~
,Z
.
,-&
~p~
E
I
•
~I'.D
A'Er
~.
A~
(a)
(b)
Fig. 4. Identification of elasto-plastic and damage parameters. (a) Tensile stress-strain curve; (b) very low cycle fatigue test
(Ae,,°,).
D=I---
E"
Then, the damage D is plotted versus the accumulated plastic strain p as in Fig. 5 from which it is easy
to deduce:
- the damage threshold epD, defined by ep < e p D ~ 10 = 0 ;
- t h e damage strength S as S = ( t r ~ / 2 E ) ( 6 p / 6 D ) since in one dimension R,, = 1 and /3 = ( t r ~ / 2 E S ) p ;
- the critical value of the damage D~c defined by Die = Max(D).
The plastic yield threshold o-s is a matter of choice regarding the type of accuracy needed in the
calculation:
- t a k e tr~ = o f for the highest bound on the time-like loading p a r a m e t e r to crack initiation;
- t a k e tr~ = t r u for the lowest bound on the time-like loading p a r a m e t e r to crack initiation.
4. N u m e r i c a l p r o c e d u r e
4.1. M e t h o d o f integration
The method used is a strain driven algorithm: at each 'time' (or time-like parameter) increment, and
for a given total strain increment, one knows the values of the stresses and the other material model
variables at the beginning of the increment (t,) and they have to be updated at the end of the increment
(t,+t). We show here that the general numerical procedure [22] becomes particularly simple and
efficient because it is analytically derived in a closed form.
D
o
~I~D
Fig. 5. Identification of damaged parameters.
J. Lemaitre, I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
205
In the following, the subscript (n + 1) is omitted and all the variables which do not contain the
subscript (n + 1) are computed at t,,+~. Tensorial notation is used for convenience.
We first assume that the increment is entirely elastic. So, we have
& = A tr el + 2/~(e - e,P,),
where A and /~ are the Lamd coefficients
Ev
A-(l_2v)(l+v),
E
~=2(1+v)
and 1 is the second order identity tensor.
All the other 'plastic' variables are equal to their values at t,,. If this 'elastic predictor' satisfies the
yield condition f ~ 0, then the assumption is valid 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 rate relations of Section 3 are discretized in an incremental form corresponding to a fully implicit
integration scheme which presents the advantage of unconditional stability [23]. Then the solution at
t,,+t has to satisfy the following relations:
~
f = °'eq
- ~s = 0 ,
& = 3, tr el + 2/~(e - %P
-- AEP),
AEp = N A p '
Y
AD = - -s Ap '
~
where A ( e ) = ( o ) , + l - ( e ) ,
and N = ~ 3( ~ ~ D /O-eq).
If we replace Ae p by its expression in the second
equation, we see that the problem is reduced to the first two equations for the two unknowns, & and p.
Let us find ~ and p which satisfy
f=
~eq -- O's = 0 ,
h = ~-
A tr el - 2 ~ ( e -
e~) + 2 ~ N A p = 0 .
This non-linear system is solved iteratively by Newton's method. For each iteration (s), we have
f+
:Ca=O,
Oh
Oh
h +-~-~: Ca +-Tpp Cp = 0 ,
where f, h and their partial derivatives are taken at t,,+t and at the iteration (s). The 'corrections' Ca
and Cp are defined by
C a = (¢7)~+, - (&)s
and
Cp = (P)~+t - (P)s.
The starting iteration (s = O) corresponds to the elastic predictor.
The set of equations for f and h may be rewritten as
f +N:Ca=O,
E
h+ H+2txAp-~--~ :Ca+21~NCp=O,
where H is the fourth order identity tensor and
ON
_1 [ 3 ( H _ I _ ~ I ® I ) _ N ~ N
0{T
O'eq
1.
Since N : N -- - 73 and N:ON/Oo" = 0, the system gives explicitly the correction for p:
Cp -
f -N:h
3~
We shall prove that the correction for & can also be found explicitly. First let us prove that Ca is a
deviatoric tensor.
From the second equation of the system and using the expression of ON~06", one obtains
tr(Ca) = - tr h .
J. Lernaitre, I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
206
Using the definitions of h and C a,
tr[(tT)~, l] = (3A + 2~) tr e = tr[(~)0 ] = constant.
These relations being true for any iteration (s), then tr(Ca) = 0 and the proof is achieved.
Using this last result together with the expressions of ON~aS" and Cp yields for the second equation of
the system,
2
A :C a = - h - ~ ( f - N : h ) N ,
where
3-~ Ap) 1 1 - 2 P ' A p N ® N .
a=(l+
U~q
/
tr~q
One can also check the following result: the tensor A is invertable if and only if t~q ~ 0, and in this case
A- ~
1
3~z
-
[
11+ 2~
~ ApN®N ] .
1+~
Ap
O-~q
O-~q
Using this expression gives explicitly the correction for t~
C a = - - 2~ ( f - N'h)N
l
1 + 3/z Ap
Crq
[h + 2tZAp(N:h)N ] .
O-~q
To summarize, this method is a fully implicit scheme with the advantage of an explicit scheme; no linear
system is to be solved and the unknown are updated explicitly.
Once p and ~7 are found, e p and D are calculated from their discretized constitutive equations and
the stress components are given by % = (1 - D ) ~ i r
Let us remark that the method described above can be applied to the case of non-linear isotropic and
kinematic hardenings, with or without damage.
4.2. Jump in cycles procedure
For periodic loadings of fatigue with large number of cycles, the computation step by step in time
becomes prohibitive and the complete computation is out of the reasonable occupation of any
computer. For that reason, a simplified method has been proposed [24] and a slightly different version
is used in this work. It allows one to 'jump' large numbers of cycles for a given approximation to the
final solution.
(a) Before any damage growth (i.e. when p < P o ) :
(i) the calculation is performed until a stabilized cycle N~ is reached. Let 8p be the increment of p
over this cycle. We assume that during a number AN of cycles, p remains linear versus N with
the slope 8p/8N. This number AN is given by
AN=
ap
8p"
where Ap is a given value which determines the accuracy on accumulated plastic strain.
(ii) a jump AN of cycles is performed and p is updated as p(N~ +AN)=p(Ns)+AN.Sp. Go to (i)
(b) Following damage growth (i.e. when p ~>PD):
(i) the calculation is performed with a constant value of the damage (/) = 0) until a stabilized cycle
N~ is reached. Then a fully coupled elasto-plastic and damage computation is performed for the
next cycle. Let 8p and 8D be the increments of p and D for this cycle, respectively. We assume
that during a number AN of cycles, p and D remain linear versus N with the slopes 8p/SN and
8D/fN, respectively. This number AN is given by
J. Lemaitre. I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
207
AN = min{ Ap
z~D} ,
~p, ~-~
where AD is a given value which determines the accuracy on the damage.
(ii) A jump AN of cycles is performed and p and D are updated as
p ( N S + 1 +AN) = p ( N s + 1) + A N . 6 p
D(Ns + 1 +AN) = D(N., + 1) + A N . 6D
Go to (i)
The choice of Ap and AD is of first importance. For this work, we chose
Dl~
-S
AD = - ~ - and ~p = y(A----e)~D, where Ae is the strain amplitude.
This method gives good results but as any heuristic method, sometimes it fails and the values of AD
or Ap must be changed.
4.3. The post processor D A M A G E 90
A computer program called D A M A G E 90 has been written as a 'friendly'. code on the basis of the
method developed. The input data are the material parameters and the history of the total strain
components eq. D A M A G E 90 gives as output the evolution of the damage D, the accumulated plastic
strain p, and the stresses o-q up to crack initiation. D A M A G E 90 distinguishes between two loading
cases:
(a) General loading history: in this case, the history is defined by the values of eq at given 'time' values.
D A M A G E 90 assumes a linear history between two consecutive given values.
(b) Piecewise periodic history: in this case, the loading is defined by blocks of cycles. In each block,
each total strain % varies linearly between two given 'peak' values. D A M A G E 90 asks the user if
he wishes to perform a complete calculation, or a simplified one using the jump in cycles procedure
described in Section 4.2 if the number of cycles is large.
D A M A G E 90 is written in F O R T R A N 77 as available on a C O N V E X computer, and contains about
600 F O R T R A N instructions. It is designed to be used as a post processor after a FEM code. However,
it may be used in an interactive way. In all the examples of Section 5, the average computing time on
the C O N V E X machine was between 1 and 15 s for each run.
The questions asked by D A M A G E 90 for a case corresponding to an example of Section 5.4 are
listed in Appendix A together with the results. Appendix B is the complete Fortran listing of the
program D A M A G E 90.
5. Examples
The following examples were performed with the D A M A G E 90 code in order to point out the main
properties of the constitutive equations used together with the locally coupled method.
5.1. Case of ductile damage
The material data are those of a stainless steel: E --- 200 000 MPa, ~, = 0.3, o-t = 200MPa, O-y=
300 MPa, o"u = 500 MPa, S = 0.06 MPa, epo = 10%, Die = 0.99.
Pure tension at microscale
The pure tension case at microscale is obtained by giving PD
the main input. D A M A G E 90 computes the other strains
=
ePD
=
10%, D c = D l c = 0.99 and ell as
J. Lemaitre, L Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
208
1
E22 = 633 = -- ~ e t t
+
-- v
){6,
E
D)
'
which represents the elasto-plastic contraction resulting in the incompressible plastic flow hypothesis
--
33
~_ E l l ) "
For a better representation of the large strain hardening of that material, the piecewise perfect
plasticity procedure is used with the following data:
ell ( % )
0
tr~ (MPa)
0.25
200
1.5
300
5
400
500
The results are shown in Fig. 6 where the strain softening and the damage are linear functions of the
strain as it is introduced in the constitutive equations.
Plane strain at meso (or micro)scale
Those cases of two-dimensional loadings were performed in pure perfect plasticity with a plasticity
threshold trs -- ~. = 500 MPa.
The influence of the triaxiality (ratio trH/treq) upon the accumulated plastic strain to rupture PR was
obtained by calculations for which in each of them trH/treq is fixed to a certain value. As from the
constitutive equations, trn/C%q = [E/(1 - 2u)](e./o's); this fixes the value of e~l + e22 = 3e H. Then, for
each point el, was chosen and e22 taken as e22 = 3 e H - - e l l .
The points of the limit curve in strains corresponding also to rupture were obtained with proportional
loading in strain histories: t~22= aett, 633 = 0. The results are plotted in Figs. 7 and 8, showing the
classical strong effect of the triaxiality which makes the rupture more and more brittle ( P R - - P D ) - - - ~ 0 [2]
and the classical 'S' shape of the limit curve of metal forming [25] here in plane strains.
5.2. Case of brittle damage
The material data are those of a ceramic: E = 400 000 MPa, u = 0.2, trf = 200 MPa, try = 250 MPa,
tru = 300 MPa, s = 0.00012 MPa, epD = 0, D l c = 0.05.
O"
M~
400
300
200
~.-~o
/
500
I-
S
.8
.5.
.6
100
.4
.2
0
10
15
20
0
I
1
I
2
T
3
4
~
a.
o',q
Fig. 6. Elastoplasticity and ductile damage in pure tension at microscale.
Fig. 7. Effect of the triaxiality upon the strain to rupture, epR is the strain to rupture in pure tension tru/trcq = 1/3, ePR = 19.6%.
J. Lemaitre. I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
209
,,
Fig. 8. Crack initiation limit curve in plane strains.
An elastic tension case is considered at mesoscale, e,~ is the main input and the two other principal
strain c o m p o n e n t s are imposed as e22 = ~33 = - I " E º I .
At microscale, this induces a pure tension in the elastic range but a three-dimensional state of stress
occurs as soon as the plastic threshold is reached due to the difference of the elasto-plastic contraction
of the inclusion and the elastic contraction imposed by the matrix.
The plastic threshold is taken as trs = 200 MPa.
The results are shown in Fig. 9 for the behavior at mesoscale where no appreciable plastic strain
appears and in Fig. 10 at microscale where the effect of damage is sensitive.
5.3. L o w and high cycle fatigue
In the damage model, there is no difference in the equation between low and high cycle fatigue,
which obey the same energetic mechanisms, according to the hypothesis. Furthermore, written as a
d a m a g e rate equation, it is valid even when a cycle cannot be defined. The material data are those of an
0"11
300 _
200
tO0
I
1
I
2
i
3
f
4
f
5
I
6
I
7
i
8
~1 i (1 if.a)
Fig. 9. Stress-strain and brittle damage at mesoscale in pure tension.
210
J. Lemaitre, 1. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
O',q M/:b
200
-
IOO
,0B
0
0
1
2
3
,4
5
6
7
~ t( I 0-")
Fig. 10. Stress-strain and damage at microscale when pure tension is applied at mesoscale.
aluminum alloy: E = 72 000 MPa, v = 0.32, crf = 303 MPa, G = 306 MPa, cru = 500 MPa, s = 6 MPa,
epD = 10%, Die = 0.99.
As for the example of Section 5.1, pure tension cases at microscale are considered by giving as input
data etl, the other strains being computed by D A M A G E 90.
In order to take into account the cyclic strain hardening the following plastic thresholds taken from a
cyclic stress-strain curve are considered for the different constant strain amplitudes imposed:
e~, ( % )
+-0.425
+-0.43
+-0.47
+-i
+3.5
+-4.5
cr (MPa)
303
305
308
370
440
460
The n u m b e r of cycles to rupture obtained are given as a function of the strain amplitude imposed in
Fig. 11 in which the number of cycles to damage initiation N o ( p = P D ) is also reported. As shown by
experiments, the ratio N o / N R increases as N R increases.
The details of the results corresponding to the case Ae~ = 7% are given in Fig. 12, showing the
evolution of the strain ~:tt and the stress-strain loops (~t~, eta) as a function of the n u m b e r of cycles,
the equivalent stress ~ q and the accumulated plastic strain p together with the damage D.
As,,~
\
\
\\ No
\
\
\
\
\,,,,
,
IO
i
,
IO~
1o ~
llo,
i
Io ~
Fig. 11. Fatigue rupture curve in pure tension at microscale.
J. Lemaitre, I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
0.0 e
2,
211
600.00
400.00
0.0:
200.00
0.0(
0.00
-200.00
-0.0
-400.00
-0.O
-600,00
-800.00
-C.vu
-0.0(
o.q
500.00
-UmU~ -u.u~
u.uo
0.02
0.04
/~
5,00
400.00
4,00 I
300.00
3.00 ~
200.00
2.00 '
tO0.O0
1.00
0.00
O,uu
! O,0O
20.00
30.00
40,00
,50.00
/v
0.00
0.00
P
tO.O0
20.00
30.00
......... )
4 0 . 0 0 50.00"
N
Fig. 12. Results of a very low cycle fatigue in pure tension at microscale Ae~ = 7%, N R = 40 cycles.
It is i n t e r e s t i n g to c o m p a r e this with the case of a p u r e elastic t e n s i o n case at m e s o s c a l e i n d u c i n g a
h r e e - d i m e n s i o n a i state o f stress at m i c r o s c a l e . T h e s a m e input are u s e d as for the case o f Fig. 12 e x c e p t
h a t t h e strains a p p l i e d are et~ = ---3.5%, e22 = e33 = - v e l t . T h e f o u r similar g r a p h s a r e s h o w n in Fig.
3.
~.4. Non-linear accumulation in fatigue
A n i m p o r t a n t f e a t u r e of fatigue d a m a g e is the n o n - l i n e a r i t y o f the a c c u m u l a t i o n o f d a m a g e s d u e to
o a d i n g s o f d i f f e r e n t a m p l i t u d e s . If a two level l o a d i n g history is c o n s i d e r e d with N l cycles at the strain
l m p l i t u d e Aet and N 2 cycles at the strain a m p l i t u d e Ae 2 with N t + N 2 = N R the n u m b e r o f cycles to
u p t u r e , in c o n t r a s t with the P a l m g r e e n - M i n e r rule, g e n e r a l l y
NI
NRt(Ael)
N1
NR, (AEI)
+
+
N,
NR2(Ae2)
N,
NR,(Ae_,)
< 1,
i f A e t > Ag 2 ,
>1,
i f A e t < A e 2,
J. Lemaitre. 1. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
212
~q
nRI
H~
500.00
1500.00
1000.00
400.00
[500.00
500.00
pJ
0.00
-500.00
200.00
,d
f
-1000.00
100.00
-1500.00
o.o~.~6 .
. . . . . . . .
2.00 . . . . . '4:66 . . . . . . '6.66
'"
-2000.00 ,,,
-0.06
a:60
N
004
0.80
0.02 l
o60
-0,04
-0.02
0.00
0.02
0.04
P
0.00
-0.02
0.20
-0.04
-0.06
0.00
. . . . . . . .
~ . . . . . . . . .
2.00
, . . . . . . . . .
4.00
, . . . . . . . . .
6.00
~
0.00
B.00
.......
0.00
,,
.........
2.00
i .........
4.00
i .........
6.00
i
8.00
iV'
Results of a very low cycle fatigue in pure tension at mesoscale inducing three-dimensional fatigue at microscale
Fig.
13.
Aet~
= 7%,
N R = 8.
NR~ and NR: being the number of cycles to failure corresponding to the constant amplitude cases of
fatigue respectively at Ae~ and Ae 2.
This effect was studied with the same material data as those in Example 5.3 for a one-dimensional
case at mesoscale • 2 2 = ~ ' 3 3 = - - V E I l "
e~ = -+0.47%
for the highest level,
e~t = ---0.425% for the lowest level,
N n = 7720 cycles,
NR
=
109 570 cycles,
The results are given in the accumulation diagram of Fig. 14; they are in good qualitative agreement
with experimental results.
5.5. Multi-axial fatigue
Two cases of multi-axial fatigue were studied: bi-axial fatigue and fatigue in tension and shear. Both
were performed with the material data of Examples 5.3 and 5.4. The game played with D A M A G E 90
consisted of finding the state of strains which produces a given number of cycles to failure. This was
done by repeated trial from many calculations.
J. Lernaitre. I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
1%
vv!N
1.
¢ r=
.2
i
.6
I °
.4
i ~
.8
1
Fig. 14. Accumulation diagram for two level tests in tension at mesoscale.
Bi-axial fatigue in plane strain at meso (or micro)scale
A plane state of strain is considered with the following in phase strains:
6"11 "= ± X ,
622 = - - Y ,
fi'33 = 0 .
Fig. 15 shows contours of cycles to failure corresponding to N R = 7800 cycles•
Fatigue in tension and shear strains imposed at meso (or micro)scale
This is the same type of calculation with
e~ = - x ,
e~2 = -+y,
all other components = 0.
"~, s
.7
/
/
oS
/
.4
NR =
/
.3
.2
/
.I
/.
/
0
i
i
!
.!
.2
.3
.4
.S
.li
.7
Fig. 15. Biaxial fatigue in plane strain.
7(300~.;,,~-1~
213
214
J. Lemaitre, 1. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
,41_
.3.
Na= 7800
¢.7¢.1~s
,2_
,1_
0
,
.t
,
.2
.3
.4
.s
.e z~Si ~$
2
Fig. 16. Fatigue in tension and shear strains imposed.
Fig. 16 summarizes the results for the same number of cycles to failure.
References
[1[ J. Lemaitre and D. Baptiste, On damage criteria, Proc. NSF Workshop on Mechanics of Damage and Fracture, Atlanta, GA
(1982).
[2] F. McClintock, A criterion for ductile fracture by the growth of holes, ASME J. Appl. Mech. (1968).
[3] J. Rice and D. Tracey. On ductile enlargement of voids in triaxial stress fields, J. Mech. Phys. Solids 17 (1969).
[4] S.S. Manson, Behavior of materials under conditions of thermal stress, NACA Tech. Rep. 1170, 1954.
[5] F.F. Coffin, Study on the effects of cyclic thermal stresses in a ductile metal, Trans. ASME (1954) 76-931.
[6] M.A. Miner, Cumulative damage in fatigue, J. Appl. Mech. 12 (1945) no. 3.
[7l N.J. Hoff, J. Appl. Mech. 20 (1953) 105-108.
[8] D. Hayhurst, J. Phys. Solids 20 (1972) 381.
[9] D. Krajcinovic and G.U. Fonseka, The continuous damage mechanics of brittle materials - Parts I and II, J. Appl. Mech. 48
( 1981 ) 809.
[10] L.M. Kachanov, On the creep rupture, Izv. SSSR, Otd. Tekhn. Nauk 8 (1958) 26.
[11] Y.N. Rabotonov, On the equations of state for creep, in: Progress in Applied Mechanics, Prager Anniversary Volume
(Macmillan, New York, 1963) 307.
[12] J. Lemaitre, Evaluation of dissipation and damage in metals submitted to dynamic loading, Proc. ICMI, Kyoto, Japan
(1971).
[13] D. Marquis and J. Lemaitre, Modelling complex behaviors of metals by the "State Kinetic Coupling" theory, Proc. ASME
Winter Meeting, San Francisco (1989).
[14] J.L. Chaboche, Description thermodynamique et ph6nom6nologique de la viscoplasticit6 avec endommagement, These de
Doctorates Sciences, Universit~ Paris 6, 1978.
[15] C. Lienard, Plasticit6 coupl6e a I'endommagement en conditions quasi-unilat6rales pour la pr6vision de I'amorqage des
fissures, These, Universit6 Paris 6, 1989.
[16] J. Lemaitre, Micro-mechanics of crack initiation, Int. J. Fracture 42 (1990) 87-99.
[17] A. Benallal, Thermoviscoplasticit6 et endommagement des structures, These de Doctorat es Sciences, Universit6 Paris 6,
1989.
[18] G.I. Taylor, Plastic strains in metals, J. Inst. Met. (1938) 62-307.
[19] J. Lemaitre and J.L. Chaboche, Mechanics of Solid Materials (Cambridge University Press, 1990).
[20] R. Billardon and I. Doghri, Localization bifurcation analysis for damage softening elasto-plastic materials, in: J. Mazars and
Z.P. Bazant, eds., Strain Localization and Size Effect due to Cracking and Damage (Elsevier, Amsterdam, 1989) 295-307.
[21] J. Lemaitre and J. Dufailly, Damage measurements, Eng. Fract. Mech., 28 (1987) 643-661.
[22] A. Benallal, R. Billardon and I. Doghri, An integration algorithm and the corresponding consistent tangent operator for
fully coupled elastoplastic and damage equations, Comm. Appl. Numer. Methods 4 (1988) 731-740.
[23] M. Ortiz and E.P. Popov, Accuracy and stability of integration algorithms for elastoplastic constitutive equations, Internat.
J. Numer. Methods Engrg. 21 (1985) 1561-1576.
[24] R. Billardon, Etude de la rupture par la m6canique de l'endommagement, These de Doctorates Sciences, Universit6 Paris 6,
1989.
[25] J.-P. Cordebois, Criteres d'instabilities plastiques et endommagement ductile en grandes d6formations. Applications a
l'emboutissage, These de Doctorates Sciences, Universit6 Paris 6, 1983.
J. Lemaitre, 1. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
Appendix A
Input data for Example 5.4
Give material constants and the strains history
* * * DAMAGE90 * * *
will give you the damage growth up to crack initiation
Elasticity. Give: Young's modulus, Poisson's ratio
72.E + 3, 0.32
Perfect plasticity: plastic threshold stress SIGs given with loading.
Give: fatigue limit SIGf, yield stress SIGy, ultimate stress SIGu
303., 306., 500.
Damage evolution: d D = (Y/S)dp. Give: S
6.
Damage threshold dD = 0 if p < pD. Do you know pD?
n
Give EpD in tension: pD = EpD * (SIGu - SIGf)/[SIGs - ( S I G f * * 2/SIGy)]
10.E-2
Crack initiation: D = Dc. Do you know Dc?
n
Give D l c in tension: Dc = D l c * [(SIGu/SIGs) * *2]/Rnu
0.99
Initial conditions. Give: p0, DO
0.,0.
Is the stress state uniaxial? 'y' or 'n'
n
Is the strain history cyclic? 'y' or 'n'
Y
• ** Y O U R L O A D I N G IS CYCLIC.
Do you wish the jump in cycles procedure for large N?
Y
Give the number of blocks of constant amplitude (50 max)
1
Give the number of cycles for block: 1
120000
Give: 1st peak, 2nd peak of
ll-strain
0.425E - 2, - 0 . 4 2 5 E - 2
Give: 1st peak, 2nd peak of
22-strain
- 0 . 1 3 6 E - 2, 0 . 1 3 6 E - 2
Give: 1st peak, 2nd peak of
33-strain
- 0 . 1 3 6 E - 2, 0.136E - 2
Give: 1st peak, 2nd peak of
12-strain
0.,0.
Give: 1st peak, 2nd peak of
13-strain
0.,0.
Give: 1st peak, 2nd peak of
23-strain
0.,0.
Give the value of the plastic threshold stress SIGs for this block
303.
Suggest a number of increments per cycle (minimum: 4)
4
*** C R A C K I N I T I A T I O N .
*** T H E JOB IS E N D E D . Y O U R R E S U L T S FILES ARE:
215
J. Lemaitre, 1. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
216
direct, out
STOP:
she ar. out
fatigue, out
Results of Example 5.4
CYCLE
0.0000E
0.2500E
0.5000E
0.7500E
0.1000E
+
+
+
+
+
00,
00,
00,
00,
01,
D
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
+
+
+
+
+
00,
00,
00,
00,
00,
0.1858E
0.1858E
0.2785E
0.2786E
0.2786E
0.2786E
0.3713E
0.3714E
0.3714E
0.3714E
0.4641E
0.4642E
0.4642E
0.4642E
+
+
+
+
+
+
+
+
+
+
+
+
+
+
04
04
04
04
04
04
04
04
04
04
04
04
04
O4
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
0.0000E
+
+
+
+
+
+
+
+
+
+
+
+
+
+
0.1086E
0.1086E
0.1086E
0.1086E
0.1086E
0.1086E
0.1086E
0.1086E
0.1096E
+
+
+
+
+
+
+
+
+
06,
06,
06,
06,
06,,
06,
06,
06,
06,
0.9858E
0.9858E
0.9858E
0.9858E
0.9858E
0.9858E
0.9858E
0.9858E
0.1000E
+
+
+
+
+
+
+
+
+
0.1100E
-
03,
MISES S T R E S S
0.0000E + 0 0 ,
0.3030E + 0 3 ,
0.3000E + 0 1 ,
0.3030E + 0 3 ,
0.3000E + 0 1 ,
O0
O0
O0
O0
O0
O0
O0
O0
O0
O0
O0
O0
O0
O0
0.2725E
0.2725E
0.4085E
0.4085E
0.4086E
0.4086E
0.5446E
0.5446E
0.5447E
0.5447E
0.6807E
0.6807E
0.6808E
0.6808E
+
+
+
+
+
+
+
+
+
+
+
+
+
+
00
00
00
00
00
00
00
00
00
00
00
00
00
00
0.3030E
0.3000E
0.3030E
0.3000E
0.3030E
0.3000E
0.3030E
0.3000E
0.3030E
0.3000E
0.3030E
0.3000E
0.3030E
0.3000E
+
+
+
+
+
+
+
+
+
+
+
+
+
+
03
01
03
01
03
01
03
01
03
01
03
01
03
01
0.3034E
0.2814E
0.3034E
0.2814E
0.3034E
0.2814E
0.3034E
0.2814E
0.3034E
0.2814E
0.3034E
0.2814E
0.3034E
0.2814E
+
+
+
+
+
+
+
+
+
+
+
+
+
+
03,
01,
03,
01,
03
01
03
01
03
01
03
01
03
01
O0
O0
O0
O0
O0
O0
O0
O0
01
0.1593E
0.1593E
0.1593E
0.1593E
0.1593E
0.1593E
0.1593E
0.1593E
0.1607E
+
+
+
+
+
+
+
+
+
02
02
02
02
02
02
02
02
02
0.4298E
0.4255E
0.4298E
0.4255E
0.4296E
0.4253E
0.4293E
0.4251E
-0.9480E
+
+
+
+
-
01
01
01
01
01
01
01
01
01
0.4303E
0.3992E
0.4303E
0.3992E
0.4301E
0.3990E
0.4298E
0.3987E
-0.9492E
+
+
+
+
-
01
01
01
01
01
01
01
01,
01 ,
P
0.0000E + 0 0 ,
0.3667E - 0 4 ,
0.3667E - 0 4 ,
0.1100E - 0 3 ,
DAM. EQU. STRESS
0.0000E + 0 0 ,
0.3034E + 0 3 ,
0.2814E + 0 1 ,
0.3034E + 0 3 ,
0.2814E + 0 1 ,
J. Lemaitre, I. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
Appendix B
C.========
-X......
C*** DAMAGE 90
C
C.1.
C.2.
C***
C***
Elastic-perfectly plastic law coupled to a ductile damage model.
Fully implicit integration scheme.
Jump in cycles procedure.
Written by Issam Doghri
Version: May 1990
............. X ......
implicit none
integer ntens, nstatv, ms, mb, naxi, ncycle, nblock, ibl
integer ist, ngiv, ig, ilookO, ipass, nul, nu2, igiv, icycle
integer nrcycl, i, ncysum, ideltn, pslope, ncyclO, kk, j
real*8 eO, xnu, sigs, smpd, sO, one, two, three, dcrit
real*8 sigf, sigy, sigu, epd, dlc, smpO, dO, energO, energ3
real*8 tolene, rac2, dtime, energl, time, tfin, strsb, star
real*8 rnu, energ2, dtmin, tper, dhtime, smpi, dhi, rtime
real*8 xxnu, ehl, epsmin, smpf, dhf, ddmaxO, yyy, trdeps, yyl
real*8 yy2, dpmax, dpmaxO, deltp, ddmax, dslope, deltd
CHARACTER*I UNIAXI,ANSPD,ANSDC,CYCLIC,JUMP,CONVER,STAB,COUPL
CHARACTER*20 FILEI,FILE2
CHARACTER*6 COMMENT
PARAMETER(NTENS=6,NSTATV=11)
REAL*8 STRS(NTENS),STATEV(NSTATV),STATEVI(NSTATV),
STRAN(NTENS),DSTRAN(NTENS),SIGHI(6),SIGHF(6),
INCUSE(50),TIM(2OO),HIST(6,2OO),EPSB(6,50),EPSS(6,50),
PER/)IV(6),SYIEL(50)
INTEGER ISLOPE(6),IPER(6),INDEX(6),NBCYCL(50)
COMMON/ETIQ1/EO,XNU,SIGS,SMPD,SO
COMMON/ETIQ2/ONE,TWO,THREE
DATA INDEX /11,22,33,12,13,23/
C*** READ DATA
MS=6
MB=5
NAXI=6
NCYCLE=O
SMPD=I.D+IO
DCRIT=O.99DO
WRITE(MS*)'
WRITE(MS*)
**********************************************************'
WRITE(MS*)
WRITE(MS*)
Give material constants and the strain history'
WRITE(MS*)
WRITE(MS*)
*** DAMAGE 90 ***'
WRITE(MS*)
WRITE(MS *) will give you the damage growth up to crack initiation'
WRITE(MS*)
WRITE(MS *)
WRITE(MS*) '** ELASTICITY . Give :'
WRITE(MS*) '
YOUNG''s modulus :'
217
218
J. Lemaitre. I. Doghri/ Comput. Me~o~ Appl. Mech. Eng~.115(19~)197-232
READ(MB,*)EO
WRITE(MS,*)'
POISSON''s ratio :'
READ(MB,*)XNU
WRITE(MS,*)'** PERFECT PLASTICITY : plastic threshold SIGs
given with loading . Give :'
WRITE(MS,*)'
Fatigue limit SIGf
:'
READ(MB,*)SIGF
WRITE(MS,*)'
Yield stress SIGy
:'
READ(MB,*)SIGY
WRITE(MS,*)'
Ultimate stress SIGu :'
READ(MB,*)SIGU
WRITE(MS,*)'** DAMAGE EVOLUTION :
dD = (Y/S) dp. Give S'
READ(MB,*)SO
FORMAT(A)
WRITE(MS,*)'** DAMAGE THRESHOLD : dD=O if p<pD. Do you know
pD?
''Y'' or ''N'''
READ(MB,I)ANSPD
IF((ANSPD.EQ.'y').OR.(ANSPD.EQ.'Y'))THEN
WRITE(MS,*)'
Give the value of pD :'
READ(MB,*)SMPD
IF(SMPD.EO.O.)SMPD=I.D-6
ELSE IF((ANSPD.EQ.'n').OR.(ANSPD.EQ.'N'))THEN
WRITE(MS,I)
2 '
Give EpD In Tension :
3 pD=EpD*(SIGu-SIGf)/[SIGs-(SIGf**2/SIGy)]'
READ(MB,*)EPD
IF(EPD.EO.O.)EPD=I.D-6
ELSE
STOP'*** WRONG DATA. GOOD BYE.'
END IF
WRITE(MS,*)'** CRACK INITIATION : D=Dc . Do you know Dc ? ''Y'' or ''N'''
READ(MB,I)ANSDC
IF((ANSDC.EQ.'y').OR.(ANSDC.EQ.'Y'))THEN
WRITE(MS,*)'
Give the value of Dc (remember: O<Dc<l) :'
READ(MB,*)DCRIT
ELSEIF((ANSDC.EQ.'n').OR.(ANSDC.EQ.'N'))THEN
WRITE(MS,*)'
Give Dlc in tension :
Dc = D l c * [(SIGu/SIGs)**2] / Rnu '
READ(MB,*)DIc
ELSE
STOP'*** WRONG DATA. GOOD BYE.'
END IF
WRITE(MS,*)'** INITIAL CONDITIONS . Give :'
WRITE(MS,*)'
The value of Po
:'
READ(MB,*)SMPO
WRITE(MS,*)'
The value of Do
:'
READ(MB,*)DO
WRITE(MS,*)'** LOADING'
WRITE(MS,*)'
Is the stress state uniaxial?
''Y'' or ''N'''
READ(MB,I)UNIAXI
IF(UNIAXI.NE.'y'.AND.UNIAXI.NE.'n'.AND.
J. Lemaitre, 1. Doghri / Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232
UNIAXI.NE.'Y'.AND.UNIAXI.NE.'N')STOP'*** WRONG DATA.
GOOD BYE.***'
IF((UNIAXI.EQ.'y').OR.(UNIAXI.Eq.'Y')) NAXI=I
WRITE(MS,*)'
Is t h e strain history cyclic? ''Y'' or'''N'''
READ(MB,I)CYCLIC
IF(CYCLIC.NE.'y'.AND.CYCLIC.NE.'n'.AND.
CYCLIC.NE.'Y'.AND.CYCLIC.NE.'N')STOP'***WRONG DATA.
• GOOD BYE.***'
IF((CYCLIC.EQ.'y').OR.(CYCLIC.EO.'Y'))THEN
WRITE(MS,*)'
'
WRITE(MS,*)'***** YOUR LOADING IS CYCLIC.'
WRITE(MS,*)'
'
WRITE(MS,*)'
Do you wish the jump in cycles procedure for
large N ? ''Y'' or ''N'''
READ(MB,I)JUMP
IF(JUMP.NE.'y'.AND.JUMP.NE.'n'.AND.
JUMP.NE.'Y'.AND.JUMP.NE.'N')STOP'*** WRONG DATA.
GOOD BYE.***'
WRITE(MS,*)'
Give the number of blocks of constant amplitude
(50 max)'
READ(MB,*)NBLOCK
DO IBL=I,NBLOCK
WRITE(MS,*)'
Give the number of cycles for block :',IBL
READ(MB,*)NBCYCL(IBL)
NCYCLE=NCYCLE+NBCYCL(IBL)
DO IST=I,NAXI
WRITE(MS,3)INDEX(IST)
FORMAT('Give: ist peak , 2nd peak
of ',12,'-strain')
READ(MB,*)EPSB(IST,IBL),EPSS(IST,IBL)
END DO
WRITE(MS,*)'
Give the value of the plastic threshold stress
& SIGs for this block :'
READ(MB,*)SYIEL(IBL)
WRITE(MS,*)'
Suggest a number of increments per cycle
(minimum : 4) : '
READ(MB,*)INCUSE(IBL)
END DO
ELSE IF((CYCLIC.EQ.'n').OR.(CYCLIC.EQ.'N'))then
WRITE(MS,*)'
'
WRITE(Ms,*)'***** YOUR LOADING IS NOT CYCLIC.'
WRITE(MS,*)'
'
WRITE(MS,*)'
Give the number of points which define the his
.tory (200 max) :'
READ(MB,*)NGIV
WRITE(MS,*)'
Give the values of time at these points :'
READ(MB,*)(TIM(IG),IG=I,NGIV)
DO IST=%,NAXI
WRITE(MS,2)INDEX(IST)
FORMAT('
Give the values of ',I2,' -strain at these'
' times :')
READ(MB,*)(HIST(IST,IG),IG=I,NGIV)
219
J. Lemaitre, l. Doghri/ Comput. Metho~ Appl. Me~. Eng~.115(1994) 197-232
220
END DO
WRITE(MS,*)'
Give the
~IGs a t t h e s e t i m e s : '
values
of the plastic
threshold
READ(MB,*)(SYIEL(IG),IG=I,NGIV)
WRITE(MS,*)'
Suggest an initial time increment
~intervalle between first 2 points'
READ(MB,*)DTIME
END IF
C***INITIALIZE. STATEV= r,p,D,dD,{pstrn},dp
DO IST=I,NSTATV
STATEVI(IST)=O.
END DO
STATEVI(2)=SMPO
STATEVI(3)=DO
STATEVI(4)=O.
DO IST=I,6
STRAN(IST)=O.
DSTRAN(IST)=O.
END DO
IF((ANSPD.EQ.'y').OR.(ANSPD.EQ.'Y'))EPD=SMPD
ENERGI=(SIGU-SIGF)*EPD
ENERGO=O.
ENERG3=O.
TOLENE=2.5D-2
ILOOKO=O
COUPL='y'
IPASS=O
ONE=I.D+O0
TWO=2.D+O0
THREE=3.D+O0
RAC2=DSQRT(TWO)
NUI=II
<
stress
to
NU2=12
312
313
314
315
FILEl='direct.out'
FILE2='shear.out'
OPEN(UNIT=NUt,file=FILEt)
OPEN(UNIT=NU2,file=FILE2)
COMMENT=' TIME'
IF((CYCLIC.EQ.'y').OR.(CYCLIC.EQ.'Y'))COMMENT=' CYCLE'
WRITE(NU1,312)COMMENT
WRITE(NU2,312)COMMENT
WRITE(90,315)COMMENT
FORMAT(2X,A6,21X,'STRAINS',27X,'STRESSES',16X,'DAMAGE',
• 8X,'p',8X, 'MISES',2X,'DAM.EQ.STRESS')
WRITE(NU1,313)
WRITE(NU2,314)
FORMAT(19X,'11',gx,'22',gx,'33',IOX,'II',9X,'22',9X,'33',
11X,'D',23X,'SIG eq',SX,'SIG*')
FORMAT(19X,'I2',gx,'13',gx,'23',10X,'I2',gx,'13',9X,'23',
11X,'D',23X,'SIG eq',5X,'SIG*')
FORMAT(SX,A6,1OX,'DAMAGE',8X,'p',5X,'MISES EQ.STRESS'
S
J. Lemaitre, l. Doghn / Comput. Me~o~ Appl. Mech. E n g ~ . l 1 5 ( 1 9 ~ ) 1 9 7 - 2 3 2
,2X,'DAMAGE EO.STRESS')
C*** IF THE LOADING IS NOT CYCLIC
IF((CYCLIC.EQ.'n').or.(CYCLIC.EQ.'N'))THEN
JUMP='n'
TIME=TIN(1)
DTMIN=DTIME/IO000.
TFIN=TIM(NGIV)
IGIV=I
CALL OUTPUT(NUI,NU2,TIME,STRAN,STBS,STATEV,NTENS,NSTATV,
• STRSB,STAR)
DOWRILE(TIME+DTIME.LE.TFIN)
DO IST=I,NAXI
DSTRAN(IST)=DTIME.(HIST(IST,IGIV+I)-HIST(IST,IGIV))
/(TIM(IGIV+I)-TIM(IGIV))
END DO
CONVER='y'
DO IST=I,NSTATV
STATEV(IST)=STATEVI(IST)
END DO
SIGS=SYIEL(IGIV)
IF(SIGS.LT.SIGF.OR.SIGS.GT.SIGU)STOP'*** PLEASE REVIEW THE
k
VALUES OF SIGf SIGu
and SIGs .'
CALL INTEGR(STRAN,DSTRAN,NTENS,NSTATV,UNIAXI,COUPL,
&
CONVER,STRS,STATEV,STRSB,STAR,RNU)
C*** IF NO CONVERGENCE, DIVIDE TIME INC. BY 2
IF(CONVER.EQ.'y'.AND.STATEV(II).GT.O.)THEN
ENERG2=ENERG3+(SIGS-SIGF*(SIGF/SIGY))*STATEV(11)
IF(ANSDC.EQ.'n')DCRIT=DIC*((SIGU/SIGS)**2)/RI~U
C***
IF(ILOOKO.EQ.O.AND.(ANSPD.EQ.'y'.OR.ANSPD.EQ.'Y')
.AND.STATEV(2).GT.I.O5*SMPD.OR.
ILOOKO.EQ.O.AND.(ANSPD.EQ.'n'.OR.ANSPD.EO.'N').AND.
ENERG2.GT.I.O5*ENERGI.OR.
STATEV(3).GT.I.05*DCRIT)THEN
CONVER='n'
END IF
END IF
IF(CONVER.EQ.'n')THEN
IPASS=IPASS+I
DTIME=DTIME/2.
IF CONVERGENCE
ELSE
IF((ANSDC.EO.'n'.OR.ANSDC.EO.'N')
.AND.DCRIT.GT.O.99DO)DCRIT=O.99DO
IPASS=O
DO IST=I,NSTATV
STATEVI(IST)=STATEV(IST)
END DO
DO IST=I,6
STRAN(IST)=STRAN(IST)+DSTRAN(IST)
END DO
TIME=TIME+DTIME
221
J. Lemaitre. l. Doghn / Comput. M e ~ o ~ Appl. Me~. E n g ~ . l 1 5 ( 1 9 ~ ) 1 9 7 - 2 3 2
222
CALL OUTPUT(NU1,NU2,TIME,STRAN,STRS,STATEV,NTENS,NSTATV,
STRSB,STAR)
IF(STATEV(3).GE.DCRIT)THEN
WRITE(MS,*)'
'
WRITE(MS,*)'*** CRACK INITIATION.'
GO TO 308
END IF
IF(STATEV(3).GE.I.1)STOP'*** DAMAGE CANNOT EXCEED i.'
C*** FIND pD
IF(ILOOKO.EO.O.AND.STATEV(ll).GT.O.)THEN
ENERG3=ENERG2
IF((ANSPD.EO.'n'.OR.ANSPD.EQ.'N')
.AND.ENERG3.GE.ENERG1)THEN
SMPD=STATEV(2)
ILOOKO=I
ELSE IF((ANSPD.EQ.'y').OR.(ANSPD.EQ.'Y').AND.
STATEV(2).GE.SMPD)THEN
ILOOKO=I
END IF
END IF
IF(DABS(TIME-TIM(IGIV+I)).LE.DTMIN)THEN
IGIV=IGIV+I
END IF
END IF
DTIME=DTIME*I.I
IF(TIME+DTIME.GT.TIM(IGIV+I)+DTMIN.AND.IGIV+I.LE.NGIV)THEN
DTIME=TIM(IGIV+I)-TIME
END IF
IF(IPASS.EQ.13.0R.DTIME.LT.DTMIN)THEN
WRITE(MS,*)'*** NO CONVERGENCE ~
WRITE(MS,*)'DTIME=',DTIME,'IPASS=',IPASS
GO TO 308
END IF
END DO
END IF
C*** IF THE LOADING IS CYCLIC
IF((CYCLIC.EQ.'y').OR.(CYCLIC.EQ.'Y'))THEN
TPER=I.
DTMIN=TPER/40000.
TIME=O.
ICYCLE=I
NRCYCL=ICYCLE
DHTIME=O.
DO I=1,6
SIGHI(I)=O.
END DO
SMPI=SMPO
DHI=DO
IF((JUMP.EQ.'y').OR.(JUMP.EQ.'Y'))THEN
COUPL='n'
OPEN (UNIT=80,file='fatigue.out')
J. Lemaitre. l. Doghn / Comput. M e ~ o ~ Appl. Mech. E n g ~ . l 1 5 ( 1 9 9 4 ) 197-232
WRITE(80,*)'KEAL TIME
TIME
D
P'
WRITE(B0,*)'******************************'
WRITE(80,*)'BEGINNING OF THE CYCLE'
WRITE(80,*)'CYCLE (REAL)
=',NRCYCL
WRITE(BO,*)'CYCLE (MACHINE) =',ICYCLE
WRITE(80,4OO)RTIME,TIME,DHI,SMPI
END IF
CALL OUTPUT(NUI,NU2,TIME,STRAN,STRS,STATEV,NTENS,NSTATV,
STRSB,STAR)
IBL=I
DTIME=TPER/INCUSE(IBL)
DO IST=I,6
IF(EPSB(IST,IBL)*EPSS(IST,IBL).LT.O.)THEN
PERDIV(IST)=TPER/4.
ELSE
PERDIV(IST)=TPER/2.
END IF
IPER(IST)=I
ISLOPE(IST)=I
END DO
NCYSUM=O
DOWHILE(NRCYCL.LE.NCYCLE)
DO IST=I,NAXI
IFCTIME+DTMIN.LE.PER~IV(IST))THEN
DSTRAN(IST)=EPSB(IST,IBL)*DTIME/PEP~IV(IST)
ELSE
xxnu= epsb(ist,ibl)-epss(ist,ibl)
DSTRAN(IST)=ISLOPE(IST)*xxnu*2.*DTIME/TPER
END IF
END DO
CONVER='y'
DO IST=I,NSTATV
STATEV(IST)=STATEVI(IST)
END DO
SIGS=SYIEL(IBL)
EHI=SIGS/IO00.
EPSMIN=SIGS/EO/IO000.
IF(SIGS.LT.SIGF.OR.SIGS.GT.SIGU)STOP'*** PLEASE REVIEW THE
k
VALUES OF SIGf SIGu
and SIGs .'
CALL INTEGR(STRAN,DSTRAN,NTENS,NSTATV,UNIAXI,COUPL,
CONVER,STRS,STATEV,STRSB,STAR,RNU)
C*** IF NO CONVERGENCE, DIVIDE TIME INC. BY 2
IF((ANSDC.EQ.'n'.OR.ANSDC.EQ.'N').AND.
STATEV(II).GT.O.)
DCRIT=DIC*((SIGU/SIGS)**2)/RNU
IF(CONVER.EQ.'n')THEN
IPASS=IPASS+I
DTIME=DTIME/2.
C*** IF CONVERGENCE
ELSE
IF((ANSDC.EQ.'n'.OR.ANSDC.EQ.'N')
223
224
J. Lemaitre, l. Doghn / Comput. M e ~ o ~ Appl. Mech. Eng~.115(1994) 197-232
.AND.DCRIT.GT.O.99DO)DCRIT=O.99DO
IPASS=O
DO IST=I,NSTATV
STATEVI(IST)=STATEV(IST)
END DO
DO IST=I,6
STRAN(IST)=STRAN(IST)+DSTRAN(IST)
END DO
TIME=TIME+DTIME
RTIME=TIME+DHTIME
CALL OUTPUT(NUI,NU2,RTIME,STKAN,STRS,STATEV,NTENS,NSTATV,
STRSB,STAR)
IF(STATEV(3).GE.DCRIT)THEN
WRITE(MS,*)'
'
WRITE(MS,*)'*** CRACK INITIATION.'
GO TO 308
END IF
IF(STATEV(3).GE.I.I)STOP'*** DAMAGE CANNOT EXCEED I.'
DO IST=I,6
IF(DABS(TIME-IPER(IST)*PERDIV(IST)).LE.DTMIN)THEN
IPER(IST)=IPER(IST)+I
DTIME=I./INCUSE(IBL)
END IF
IF(ISLOPE(IST).GT.O.AND.DABS(STRAN(IST)-EPSB(IST,IBL)).LE.
• EPSMIN.OR.ISLOPE(IST).LT.O.AND.DABS(STRAN(IST)-EPSS(IST,IBL))
.LE.EPSMIN)THEN
ISLOPE(IST)=-ISLOPE(IST)
END IF
END DO
C*** END OF A CYCLE
IF(DABS(TIME-TPER*ICYCLE).LE.DTMIN)THEN
IDELTN=O
DO I = 1 , 6
SIGHF(I)=STRS(I)
END DO
SMPF=STATEV(2)
DHF=STATEV(3)
DELTP=SMPF-SMPI
DELTD=DHF-DHI
C * * * FIND pD
IF(ILOOKO.EQ.O)THEN
ENERG3=ENERGO+(SIGS-SIGF*(SIGF/SIGY))*
DELTP*(NRCYCL-NCYSUM)
IF((ANSPD.EQ.'n'.OR.ANSPD.EQ.'N')
.AND.ENERG3.GE.ENERGI)THEN
SMPD=SMPF
NCYCLO=NRCYCL-NCYSUM
ILOOKO=I
ELSE IF((ANSPD.EQ.'y'.OR.ANSPD.EQ.'Y').AND.
STATEV(2).GE.SNPD)THEN
NCYCLO=NRCYCL-NCYSUM
J. Lemaitreol. Doghd / Comput. Me&odsAppl. Mech. Eng~.115 (1994) 197-232
ILOOKO=I
END IF
END IF
C*** JUMP IN CYCLES PROCEDURE
IF((JUNP.EQ.'y').OR.(JUMP.EQ.~Y'))THEN
WRITE(80~*)'END OF THE CYCLE'
WRITE(80,4OO)RTIME,TIME,DHF,SMPF
DDMAXO=DCRIT/50.
IF((ANSDC.Eq.'n').OR.(ANSDC.EQ.'N'))
DDMAXO=DIC/50.
IF((UNIAXI.EQ.'y').OR.(UNIAXI.EQ.'Y'))THEN
YYY=SIGS*SIGS/2./EO
ELSE
TRDEPS=O.
DO IST=I,3
TRDEPS=TRDEPS+EPSB(IST,IBL)-EPSS(IST,IBL)
END DO
YYI=(I.+XNU)*SIGS*SIGS/3./EO
YY2=EO*TRDEPS*TRDEPS/(I.-2.*XNU)/6.
YYY=YY1+YY2
END I F
DPMAXO=SO*DDMAXO/YYY
DPMAX=DPMAXO
PSLOPE=NBCYCL(IBL)*DELTP
IF(DPMAX.GT.PSLOPE)DPNAX=PSLOPE
NRITE(80,*)'ILOOKO=J,ILOOKO,'NCYCLO=',NCYCLO,'pD=',SMPD
IF(COUPL.EQ.'n'.OR.STATEV(2).LT.SNPD)THEN
STAB=,y *
DO KK=I,6
IF(DABS(SIGHF(KK)-SIGHI(KK)).GT.EHI)STAB='n'
END DO
WRITE(80,*)'ENERGI=',ENERGI,'ENERG3=',ENERG3
IF(STAB.EQ.'y') THEN
WRITE(80,*)'STABILIZED CYCLE'
IF(STATEV(2).LT.SMPD)THEN
C*** JUMP OF CYCLES BEFORE DAMAGE GROWTH
IDELTN=IDINT(DPMAX/DELTP)
IF(NRCYCL+IDELTN.GT.NCYSUM+NBCYCL(IBL))
IDELTN=NCYSUM+NBCYCL(IBL)-NRCYCL
ENERG2=ENERGO+(SIGS-SIGY)*DELTP*(NRCYCL+IDELTN-NCYSUM)
IF(ENERG2.GT.ENERGI*I.05)
IDELTN=(ENERGI-ENERG3)/(SIGS-SIGY)/DELTP
WRITE(80,*)'*** JUMP OF ',IDELTN,' CYCLES'
WRITE(80,*)'DPMAX=',IDELTN*DELTP, 'YYY=',YYY
ELSE
NRITE(80,*)'COUPLED COMPUTATION FOR NEXT CYCLE '
COUPL='y'
END IF
END IF
ELSE
C*** JUMP OF CYCLES AFTER DAMAGE GROWTH
225
J. Lemaitre. I. Doghn / Comput. Metho~ Appl. Mech. E n g ~ . l l 5 (19~)197-232
226
DDMAX=DDMAXO
DSLOPE=(NBCYCL(IBL)-NCYCL0)*DELTD
IF(DDMAX.GT.DSLOPE)DDMAX=DSLOPE
NCYCLO=O
IDELTN=IDINT(DDMAX/DELTD)
IF(NRCYCL+IDELTN.GT.NCYSUM+NBCYCL(IBL))
IDELTN=NCYSUM+NBCYCL(IBL)-NRCYCL
IF(DHF+DELTD*IDELTN.GT.DCRIT*I.05)
• IDELTN=IDINT((DCRIT-DHF)/DELTD)
IF(IDELTN.GT.IDINT(DPMAX/DELTP))IDELTN=IDINT(DPMAX/DELTP)
WRITE(80,*)'*** JUMP OF ',IDELTN,' CYCLES'
WRITE(80,*)'DDMAX=',IDELTN*DELTD, 'DPNAX=',IDELTN*DELTP
COUPL='n'
END IF
IF(DELTP.GT.DPMAX)WRITE(80,*)'*** PROBLEM :
DPMAXO IS TOO SMALL.'
IF(DELTD.GT.DDMAX)WRITE(80,*)'*** PROBLEM :
DDMAXO IS TOO SMALL.'
IF(DELTD.GT.DDMAX.OR.DELTP.GT.DPMAX)IDELTN=O
IF(IDELTN.LT.O.)STOP'*** PROBLEM : NEGATIVE JUMP OF CYCLES.'
END IF
DO J=l,6
SIGHI(J)=SIGHF(J)
END DO
SMPI=SMPF+DELTP*IDELTN
DHI=DHF+DELTD*IDELTN
STATEVI(2)=SMPI
STATEVI(3)=DHI
DHTIME=DHTIME+TPER*IDELTN
NRCYCL=NRCYCL+IDELTN
RTIME=TIME+DHTIME
IF((JUMP.EQ.'y').OR.(JUMP.EQ.'Y'))THEN
WRITE(80,*)'
WRITE(80,*)'BEGINNING OF THE CYCLE'
WRITE(80,*)'CYCLE (REAL)
=',NRCYCL+I
WRITE(80,*)'CYCLE (MACHINE) =',ICYCLE+I
WRITE(80,4OO)RTIME,TIME,DHI,SMPI
END IF
IF(NRCYCL.EQ.NCYSUM+NBCYCL(IBL))THEN
ENERGO=ENERGO+(SIGS-SIGF*(SIGF/SIGY))*DELTP*(NRCYCL-NCYSUM)
NCYSUM=NCYSUM+NBCYCL(IBL)
IBL=IBL+I
END IF
ICYCLE=ICYCLE+I
NRCYCL=NRCYCL+l
END IF
•
400
FORMAT(1X,4(E12.6,','))
END IF
DO IST=I,6
IF(TIME+DTIME.GT.IPER(IST)*PERDIV(IST)+DTNIN)THEN
DTIME=(IPER(IST)*PERDIV(IST))-TIME
J. Lemaitreol. Doghn / Comput. M e ~ o ~ Appl. Mech. E n g . . I f 5 ( 1 9 9 4 ) 197-232
END IF
END DO
IF(IPASS.EQ.13.0R.DTIME.LT.DTMIN)THEN
WRITE(MS,*)'*** NO CONVERGENCE'
WRITE(MS,*)'DTIME=',DTIME,'IPASS=',IPASS
GO
TO
308
END IF
END DO
END IF
308
WRITE(MS,*)'*** THE JOB IS ENDED• YOUR RESULT FILES ARE :'
WRITE(MS,*)'
'
IF((JUMP.EQ.'n').OR.(JUMP.EQ.'N'))WRITE(MS,*)FILEI,FILE2
IF((JUMP.EQ.'y').OR.(JUMP.EQ.'Y'))
• WRITE(MS,*)FILEI,FILE2,'fatigue.out'
STOP
END
C.==============
---X ......
SUBROUTINE INTEGK(STRAN,DSTRAN,NTENS,NSTATV,UNIAXI,COUPL,
&
CONVER,STRS,STATEV,STRSB,STAR,RNU)
C.=======:
-'~X......
implicit none
integer ntens, nstatv, niss, i, j, nfois, ntest
real*8 eO, xnu, sigs, smpd, sO, one, two, three, strsb, star
real*8 rnu, rac2, smr, smp, d, dsmp, treps, estrsh, estrsb
real*8 sign, f, denom, xmu, divl, div2, xnkes, dd, y
real*8 tolres, csmp, rel, restrs, dsmr
REAL*8 STATEV(NSTATV),STRAN(NTENS),DSTRAN(NTENS),
k
STRN(6),ID(6,6),AS(6),
STKS(6),ESTRS(6), CESTRS(6),KESTRS(6),
N(6),PSTRN(6),DPSTRN(6)
REAL*8 LAMDA,NN
LOGICAL CRITERE
CHARACTER*I CONVER,COUPL,UNIAXI
COMMON/ETIQI/EO,XNU,SIGS,SMPD,SO
COMMON/ETIQ2/ONE,TWO,THREE
NN=THKEE/TWO
RAC2=DSORT(TWO)
TOLRES=SIGS*I.D-6
NISS=NTENS
C*** IDENTITY MATRIX
DO I=I,NISS-1
DO J=I+I,NISS
ID(I,J)=O.
ID(J,I)=O.
END DO
ID(I,I)=ONE
END DO
ID(NISS,NISS)=ONE
C*** IDENTITY VECTOR
DO I=I,NISS
227
228
J. Lemaitre, l. Doghri/ Comput. Metho~ Appl. Mech. Eng~.l15(1994) 197-232
IF(I.LE.3)AS(I)=ONE
IF(I.GE.4)AS(I)=O.
END DO
C*** LAME COEFF.
XMU=EO/(ONE+XNU)/TWO
LAMDA=EO*XNU/(ONE-TWO*XNU)/(ONE+XNU)
SMR=STATEV(1)
SMP=STATEV(2)
D=STATEV(3)
DO I = I , N I S S
PSTRN(I)=STATEV(I+4)
END DO
C*** TOTAL STRAIN
DO I = I , N I S S
STRN(I)=STRAN(I)+DSTRAN(I)
IF(I.GE.4)STRN(I)=STRN(I)*RAC2
END DO
C*** ELASTIC PREDICTOR
DSMP=O.
IF((UNIAXI.EQ.'y').OR.(UNIAXI.EQ.'Y'))THEN
STRN(2)=-XNU*STRN(1)-(O.S-XNU)*PSTRN(1)
STRN(3)=STRN(2)
END IF
C . . . . T r a c e of strain tensor
TREPS=STRN(1)+STRN(2)+STRN(3)
C .... Effective stress
DO I=I,NISS
ESTRS(I)=LAMDA*TREPS*AS(I)+2.*XMU*(STRN(I)-PSTRN(I))
END DO
CALL INVAR(ESTRS,ESTRSH,ESTRSB,NISS)
SIGN=I.
IF(ESTRS(1).LT.O.)SIGN=-I.
C*** TEST THE YIELD CONDITION
F=ESTRSB-SIGS
IF(F.LE.O.)THEN
GO TO 500
ELSE
C*** IF THE ELASTIC PREDICTOR DOES NOT VERIFY THE YIELD CONDITION
DO I=I,NISS
N(I)=NN*(ESTRS(I)-ESTRSH*AS(I))/ESTRSB
END DO
C*** Corrections over the elastic predictor
c .... Correction OVER P
DENOM=3.*XMU
CSMP=F/DENOM
C .... Correction over effective stress
DO I=I,NISS
CESTRS(I)=-2.*F*N(I)/3.
END DO
C*** LOOP ON THE PLASTIC CORRECTIONS
]. Lemaitre, 1. Doghri/ Comput. Me~o~ Appl. Mech. Eng~.115(19~)197-232
NFOIS=O
CRITERE=.TRUE.
DO WHILE(CRITERE)
NFOIS=NFOIS+I
DIVI=ESTRSB
NTEST=2
IF(NFOIS.EQ.I)NTEST=I
C*** UPDATE THE STATE
C .... Update inc. of p
DSMP=DSMP+CSMP
C .... Update effective stress
DO I=I,NISS
ESTRS(I)=ESTRS(I)+CESTRS(I)
END DO
CALL INVAR(ESTRS,ESTRSH,ESTRSB,NISS)
DO I=I,NISS
N(I)=NN*(ESTRS(I)-ESTRSH*AS(I))/ESTRSB
END DO
IF(NTEST.EQ.2)THEN
DIV2=ESTRSB
REL=(DIV2-DIV%)/DIVI
END IF
C .... If we converge too slowly,or we diverge,exit.The main
C .... program will propose a smaller "time" increment.
IF(NFOIS.EQ.50.OR.
(NTEST.EQ.2.AND.REL.GE.15.D-O2))THEN
CONVER='n'
GO TO 500
END IF
C*** Compute residual functions
F=ESTRSB-SIGS
IF((UNIAXI.EQ.'y').OR.(UNIAXI.EQ.'Y'))TBEN
STRN(2)=-O.5*STRN(1)+(O.5-XNU)*SIGN*SIGS/EO
STRN(3)=STRN(2)
TREPS=STRN(1)+STRN(2)+STRN(3)
END IF
DO I=I,NISS
KESTRS(I)=ESTRS(I)-LAMDA*TREPS*AS(I)2.*XMU*(STRN(1)-PSTRN(1)-DSMP*N(I))
END DO
C .... Max. residuals
RESTRS=O.
DO I=I,NISS
IF(DABS(KESTRS(I)).GT.RESTRS)RESTKS=DABS(KESTRS(I))
END DO
C*** ITERATIVE TEST ON THE YIELD CONDITION
IF(DABS(F).LT.TOLRES.AND.RESTRS.LT.TOLRES)THEN
CRITERE=.FALSE.
ELSE
C*** PLASTIC CORRECTIONS
229
230
J. Lemaitre, l. Doghn / Comput. Metho~ Appl. Mech. E n g ~ . 1 1 5 ( 1 9 ~ ) 1 9 7 - 2 3 2
CALL VTRANVI(N,KESTRS,XNKES,NISS)
C .... Correction over p
DENOM=3.*XMU
CSMP=(F-XNKES)/DENOM
C .... Correction over effective stress
DENOM=DENOM*DSMP/ESTRSB+(1.)
DO I=I,NISS
CESTRS(I)=(2./3.)*(XNKES-DENOM*F)*N(I)-KESTRS(I)
CESTRS(I)=CESTRS(I)/DENOM
END DO
C*** END OF THE ITERATIVE TEST ON THE YIELD CONDITION
END IF
C*** END OF THE LOOP ON THE PLASTIC CORRECTIONS
END DO
C*** END OF TEST ON THE YIELD CONDITION
END IF
500
CONTINUE
IF(CONVER.EQ.'y')THEN
IF((UNIAXI.EQ.'y').OR.(UNIAXI.EQ.'Y'))THEN
DSTRAN(2)=STRN(2)-STRAN(2)
DSTRAN(3)=DSTRAN(2)
END IF
C .... plastic strain increment
DO I=I,NISS
DPSTRN(I)=N(I)*DSMP
END DO
C .... Damage increment
DD=O.
CALL DAMAGE(TREPS,ESTRSB,Y,RNU)
IF(DSMP.GT.O..AND.SMP+DSMP.GE.SMPD.AND.COUPL.EQ.'y')THEN
DD=Y*DSMP/SO
END IF
STAR=(ONE-D-DD)*DSQRT(2.*EO*Y)
C .... plastic multiplier increment
DSMR=(ONE-D-DD)*DSMP
C .... Compute Stresses
DO I=I,NISS
STRS(I)=ESTRS(I)*(ONE-D-DD)
END DO
STRSB=ESTRSB*(ONE-D-DD)
C*** STORE STATEV AT THE END OF THE INCREMENT
STATEV(1)=SMR+DSMR
STATEV(2)=SMP+DSMP
STATEV(3)=D+DD
STATEV(4)=DD
DO I=I,NISS
STATEV(I+4)=PSTRN(I)+DPSTRN(I)
END DO
STATEV(II)=DSMP
END IF
RETURN
J. Lemaitre, l. Doghri/ Comput. Metho~ Appl. Mech. Eng~.115(19~)197-232
END
C*=================
SUBROUTINE INVAR(V,VH,VB,NISS)
C,=======:
C... ist and 2nd stress invariants
implicit none
integer
niss, i
real*8
one, two, three, vh, vb, const
REAL*8 V(6)
CONMON/ETIQ2/ONE,TWO,THREE
-'X ......
"-~=========X
VH=(V(1)+V(2)+V(3))/THREE
CONST=O.
DO I=I,NISS
IF(I.LE.3)CONST=CONST+(V(1)-VH)*(V(1)-VH)
IF(I.GE.4)CONST=CONST+V(I)*V(I)
END DO
VB=DSQRT(THREE*CONST/TWO)
RETURN
END
C *===============~ ...................................................
SUBROUTINE DAMAGE(TREPS,ESTRSB,Y,RNU)
implicit none
real*8
eO, xnu, sigs, smpd, sO, treps,
real*8
term2
COMMON/ETIQI/EO,XNU,SIGS,SMPD,SO
estrsb,
......
~X:
y, rnu, terml
TERMI=ESTRSB*ESTRSB*(1.+XNU)/EO
TERN2=EO*TREPS*TREPS/(1.-2.*XNU)/2.
Y=(TERMI+TERM2)/3.
RNU=2.*EO*Y/ESTRSB/ESTRSB
RETURN
END
"-X......
C~=================:
SUBROUTINE
VTRANVI(V,VI,VTV1,NISS)
::::::::::::::::::::::
C... Inner product of 2 symmetric
implicit none
integer
niss, i
real*8
vtvl
REAL*8 V(6),Vl(6)
...... -X ......
2nd order tensors
VTVI=O.
DO I=I,NISS
VTVI=VTVI+V(1)*VI(1)
END DO
RETURN
END
SUBROUTINE
OUTPUT(NUI,NU2,TIME,STRAN,STRS,STATEV,
231
232
J. Lemaitre. I. Doghri / Comput. Methot~ Appl. Mech. Eng,. 115 (19~) 197-232
NTENS,NSTATV,STRSB,STAR)
. . . . . . . . . .
C,=
implicit none
integer nul, nu2, ntens, nstatv, i
real*8
time, strsb, star
REAL*8 STRAN(NTENS),STRS(NTENS),STATEV(NSTATV)
WRITE(NUI,130)TIME,(STRAN(I),I=I,3),(STKS(I),I=I,3),
STATEV(3),STATEV(2),STRSB,STAR
WRITE(NU2,130)TIME,(STRAN(I),I=4,6),(STRS(I),I=4,6),
130
444
&
STATEV(3),STATEV(2),STRSB,STAR
WRITE(90,444)TIME,STATEV(3),STATEV(2),STRSB,STAR
FORMAT(3X,11(EIO.4,','))
FORMAT(3X,S(EIO.4,','))
RETURN
END
X ~ . . . . .
© Copyright 2024 ExpyDoc
