Computer methods in applied mechanics and englneerlng

ELSEVIER Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232

Damage 90: a post processor for crack initiation

Issam Doghri
of

~niversi~ Received 10 December

of

1992

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.

During thepast decades, the problem of macrocrack growth up to failure of structures has received much attention through the tremendous development of fracture mechanics. Crack initiation, which requires a knowledge of what happens before a mesocrack breaks the representative volume element (RVE), cannot be treated by classical fracture mechanics dealing with pre-existing meso or macro cracks. The usual engineering method of predicting the crack initiation consists of using some phenomenological or empirical criteria or relations to define the conditions of local rupture of which the following are examples: - a critical value of the maximum principal stress; - the von Mises equivalent stress or better the damage equivalent stress equal to some critical value for brittle fracture [l]; -the 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]; - the Manson-Coffin law for low cycle fatigue [4,5] - the Woehler-Palmgreen Miner rule for high cycle fatigue [6]; - for creep, a relation between the stress and the time to rupture corresponding to an infinite strain flow in a creep constitutive equation 173 or a stress criterion based on a linear combination of the stress tensor invariants [8]. The poor accuracy of such predictions in complex cases comes from the inability of those models to

* Corresponding

author.

1Visiting Professor at University of California, Santa Barbara.

take into account the history of loadings and for most of them, the difficulty of being written in three-dimensions. The tool for modeling the evolution of the progressive deterioration of materials up to a mesocrack initiation is the damage mechanics which was developed after fracture mechanics but which is now in its stage of maturity. A continuous damage variable is introduced as an internal variable. It is of tensorial nature, but it reduces to a scalar if the damage is isotropic [9] which is the hypothesis considered here. It is the maximum equivalent area density of microcracks or microcavities which lies in any cross plane that can be defined in a RVE [lo]:

where &S is the area of the most damaged cross section of the RVE and &S, is the equivalent area of the microcracks and microcavities in 6s. The concept of effective stress Gii = aij/( 1 - D) [ll] associated with the principle of strain equivalence [12] allows one to write the Helmholtz free energy t,!tin the framework of the thermodynamics of irreversible processes. The couplings or uncouplings between elasticity, plasticity and damage are identified and modeled in accordance with the state kinetic coupling theory [13] which defines a logical procedure to choose the analytical expressions of the state potential r(l and the potential of dissipation F from which are derived the kinetic constitutive equations. If small deformations and linear isotropic elasticity are assumed, I+G =5 [A &Q&z + (1 + v)Ll _ 2v) ciil (1 - D) + terms for plasticity ,

where &piis the elastic strain tensor, E and Y are the Youngs modulus and Poissons ratio and p is the mass density considered as a constant. The law of elasticity coupled to damage is derived from q,=p$=E(l-D) or l+V &?. =11 E 6, - z iQij E operator , 6, = 1 if i = j, 6, = 0 if i #j. The associated variable to D is defined by

dk4,
(1+ v)(l - 2v)

8, is the Kroneckers

Y is the strain energy density release rate [14] defining the power dissipated in the damaging process as +,=YD>O, D=$, tbeingthetime.

It can be written in a close form as

4qR,
Y=2E(1-D)2 where Use is the von Mises stress C_ = ($~ycry)~~, CJ~ is the stress deviator a: = aii - a,?jij, a, is the hydrostatic stress an = 3 wkk, and R, is the triaxiality function which depends on the triaxiality ratio u,/u_i, i.e. R, =;(I + v) + 3(12~)(?).
eq

This allows the definition of the damage equivalent stress u* as the one-dimensional stress which gives, for the same value of the damage, the same elastic strain energy density W, as the three-dimensional state [l]. Since Y = w,/(l - D), then U* = creqRt". 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, 161 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 (FEM) 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 uy, ultimate stress a, and fatigue limit of (o, < my,). - The behavior of 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 uy and equal to the fatigue limit of if not known by another consideration

Fig. 1. Locally coupled analysis of crack initiation.

200

c-

Fig. 2. Two scale representative

volume element.

Below the fatigue limit ur no damage should occur. The fatigue limit of the inclusion is supposed to be reduced in the same proportion:

The second assumption which makes the calculation simple is the Lin-Taylors 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, u&,,*, = Max(u(*,))+ 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 u* for a long time may be more damaging than a large u* for a short time! There is no rigorous way to select the critical point(s) but an intelligent look at the evolution of u* 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 6A = d. 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,

But -6 W/2 is also the energy dissipated in the inclusion by the damaging process at constant stress, DC
d31 0

YdD
6A

G (6A) =

Assuming a constant strain energy density release rate, Y = Y, yields

Y,D,d .

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,

where G, is the toughness of the material.

It yields

Comparing the two relations for shows that GcsA)= 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.

The following hypotheses are made: - small deformation; - uncoupled partitioning of the total strain in elastic strain and plastic strain eij = eQ + E:; - 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 a, is bounded by the fatigue limit and the ultimate stress

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:

The plastic strain rate derives from

f through the normality rule of standard materials [19],

where h: is the plastic multiplier.

This shows that

where 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 in the framework of the state kinetic coupling theory D31.
Y2

d&
= 2E(1-Dj2 =e

+2S(l-D) where S is a damage strength,

a constant characteristic

of each material,

/i

or

ifpsp, in Section 2).

(mesocrack

initiation if

as demonstrated

202

J. Lemaitre, 1. Doghri I Comput. Methods Appl. Mech. Engrg. 11.5 (1994) 197-232

Fig. 3. Modelling a stress-strain

curve

pb 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 oS and a fatigue limit of = of/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 v;p. Then the stored energy density is
w, =

uii

ds; - U;P plastic constitutive equation,

From the perfectly

w, = (a, - cr!)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 up,,, a plastic threshold equal to the ultimate stress mUand a fatigue limit Up, a (Us cT;)PD = tu gff)EPD and PD = EPD 3 s Uf

with of = a:/~~.

If piecewise perfect plasticity of several thresholds - urP

oSsiis considered

wS = C usi(Pi -Pi-l) and pD is defined by PgOusi(Pi

-Pi-l)

urPD

tau

a;IEPD .

The critical value of damage at mesocrack initiation D, corresponds to an instability [20]. Nevertheless, it can be related to a certain amount of energy dissipated in the damaging process D, I
0

Y dD = constant at failure

Assuming for that calculation CT2R edD=eD,.

a proportional CT2R

loading for which R, = const,

Taking again the pure tension test as a reference which gives the critical value of the damage at mesocrack initiation D,, related to the stress to rupture a, and the tensile decohesive ultimate stress oU together with R, = 1 gives the rupture criterion, (r2R --fgDc=$~,c,

that is

Collecting all the equations gives the set of constitutive

equations to be solved by a numerical analysis:

3 6;
&;= { 0,

. 2 as P,

iff=k

CT f=-&-us,

iff<O,

if pap,,
ifp<p,,

U - Uf *

Mesocrack initiation occurs if D = 1. As the plastic strain rate .$ is governed by p, and Ij 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 only as a function of the total strain, since U = E(l-D) &K7 l-2V eq = a,(1 - D) .
&H= E

+ .s&) and

0,

Then
UH -=-Ueq

E
l-2V

EH

US.
of

3.2.

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 handbook or the test results used correspond exactly to the material considered in the application. Those coefficients are: and v for elasticity; uf, uy , uu and the choice of a, for plasticity; S, cpD, D,, for damage. They can be easily obtained from a tensile stress-strain curve with unloadings or from a very low fatigue test at constant strain amplitude (NR s 10) as those of Fig. 4. A classical experimental procedure with strains measured by strain gauges permits the determination of Youngs modulus Poissons ratio v and also the damaged elasticity modulus i? as defined on the graphs of Fig. 4 and as a function of Fig. 4(a); 2N As, in Fig. 4(b) since 2u,lE( 1 2u, const, where N is the number of cycles. The yield stress cry and the ultimate stress a, are also of a classical procedure. The fatigue limit uf is the stress for which the number of cycles to failure NR in a high cycle fatigue test is very large: NR z lo6 to 10 cycles. Usually

The damaged parameters write [21]

are deduced from the measurement

of ,6 as the law of elasticity allows

204

A&

(4
Fig. 4. Identification (%,,). of elasto-plastic and damage parameters. (a) Tensile stress-strain

(b)
curve; (b) very low cycle fatigue test

D=l-g.
Then, the damage D is plotted versus the accumulated plastic strain as in to deduce: - the damage threshold defined by EPDIj = 0 ; since in one dimension -the damage strength S as S = - the critical value of the damage defined by Max(D). The plastic yield threshold a, is a matter of choice regarding the type calculation: - take crS= of for the highest bound on the time-like loading parameter to - take crS= gU for the lowest bound on the time-like loading parameter to Fig. 5 from which it is easy

1 and D = (afI2ES)@

of accuracy needed in the crack initiation; crack initiation.

Numerical procedure

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,+l). We show here that the general numerical procedure [22] becomes particularly simple and efficient because it is analytically derived in a closed form.

*D l-

Fig. 5. Identification of damaged parameters.

In the following, the subscript (n + 1) is omitted and all the variables which do not contain 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 ci:= A tr &I+ 2p(s - E:) , where A. and I_L are the Lame coefficients

the

and 1 is the second order identity tensor. All the other plastic variables arc equal to their values at t,_ If this elastic predictor satisfies the yield condition f Q 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 has to satisfy the following relations:
f = GeqCT, =

0 )

r7;=Atrel+2~(e-~~-AeP),A~*=NAp,

AD=:Ap,

and N = $(GD/GG,). If we replace Aep by its expression in the second where equation, we see that the problem isreduced to the first two equations for the two unknots, & and p. Let us find @ and p which satisfy f=K, -cls=O,

This non-linear
f

system is solved iteratively

by Newtons method. For each iteration

(s), we have

+~:c;=u, h+$:c,+$c,=o,
are taken at in+1 and at the iteration (s). The corrections C, and Cp = (P)~+~ - (P), s

where f, h and their partial derivatives and Cp are defined by C, = (G)s+l - (4,

The starting iteration (s = 0) corresponds to the elastic predictor. The set of equations for f and ii may be rewritten as
f +N:C,=u,

where 17 is the fourth order identity tensor and

Since N C J-N:h
P

and N: &Vi&% = 0, the system gives explicitly the correction

for p:

3P

We shall prove that the correction for & can also be found explicitly. First let us prove that C, is a deviator% tensor. From the second equation of the system and using the expression of aNlaG, one obtains tr(C&) = -tr h .

206

Using the definitions of h and C,, tr[( G),+,] = (3A + 2~) tr E = tr[( 6)J = constant . These relations being true for any iteration (s), then tr(C,) = 0 and the proof is achieved. Using this last result together with the expressions of aNla6 and C, yields for the second equation of the system,
A:C, = -h

-$f-N:~)N,
L+ApN@N.
eq

where
A=

One can also check the following result: the tensor A is invertable if and only if Gek, # 0, and in this case
A-=

1++p
eq

i7+$ApN@N
eq

I.
for G
1 h
-l

Using this expression

gives explicitly the correction

-$f-N:~)N

l+$Ap
eq

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 and 6 are found, E and D are calculated from their discretized constitutive equations and the stress components are given by cij = (1 - D)Gij. 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. 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 (i) the calculation is performed until a stabilized cycleN, is reached. Let be the increment of over this cycle. We assume that during a number AN of cycles, remains linear versus N with the slope This number AN is given by

where Afiis a given value which determines the accuracy on accumulated plastic strain. (ii) a jump AN of cycles is performed and is updated as +AN. to (i) (b) Following damage growth (i.e. when (i) the calculation is performed with a constant value of the damage (fi = 0) until a stabilized cycle N, is reached. Then a fully coupled elasto-plastic and damage computation is performed for the next cycle. Let and SD be the increments of and D for this cycle, respectively. We assume that during a number L\N of cycles, pnd D remain linear versus N with the slopes and respectively. This number AN is given by

where AD is a given value which determines the accuracy on the damage. (ii) A jump E of cycles is performed and p and D are updated as p(N,+l+AN)=p(N,+l)+AN+p D(N, + 1 +AN) = D(N, + 1) +[Link] Go to (i) The choice of & and AD is of first importance.
-

For this work, we chose

D,, S AD = 50 and Ap = -Y(G) AD, where As is the strain amplitude .

This method gives good results but as any heuristic method, sometimes it fails and the values of AD or G must be changed. 4.3. The post processor DAMAGE 90

A computer program called DAMAGE 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 eij. DAMAGE 90 gives as output the evolution of the damage D, the accumulated plastic strain p, and the stresses aii up to crack initiation. DAMAGE 90 distinguishes between two loading cases: (a) General loading history: in this case, the history is defined by the values of .sij at given time values. DAMAGE 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 cij varies linearly between two given peak values. DAMAGE 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. DAMAGE 90 is written in FORTRAN 77 as available on a CONVEX computer, and contains about 600 FORTRAN 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 CONVEX machine was between 1 and 15 s for each run. The questions asked by DAMAGE 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 DAMAGE 90.

5. Examples The following examples were performed with the DAMAGE 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, 300 MPa, a, = 500 MPa, S = 0.06 MPa, epD = 10%) D,, = 0.99. v = 0.3, a, = 200 MPa, ff,, =

Pure tension at microscale The pure tension case at microscale is obtained by giving pD = epD = lo%, D, = D,, = 0.99 and &11 as the main input. DAMAGE 90 computes the other strains

208

J. Lemaitre, I. Doghri ! Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232

which represents the elasto-plastic contraction resulting in the incompressible plastic flow hypothesis (& = &!& = -+&yl). For a better representation of the large strain hardening of that material, the piecewise perfect plasticity procedure is used with the following data: 0 us (MPa) 200 0.25 300 1.5 400 5 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 Those cases of two-dimensional loadings were performed in pure perfect plasticity with a plasticity threshold a, = a, = 500 MPa. The influence of the triaxiality (ratio an/a,,) upon the accumulated plastic strain to rupture pR was obtained by calculations for which in each of them mH/ueeqis fixed to a certain value. As from the this fixes the value of sI1 + sz2 = 3~~. Then, for constitutive equations, on /a,, = [E/(1 - 2v)](s,la,); each point sI1 was chosen and cz2 taken as cz2 = 3~~ - .Q. The points of the limit curve in strains corresponding also to rupture were obtained with proportional loading in strain histories: .sz2= (YEAS, Q = 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 ( pR - pD) + 0 [2] and the classical S shape of the limit curve of metal forming [25] here in plane strains. 5.2. Case brittle damage

The material data are those of a ceramic: E = 400 000 MPa, Y = 0.2, or = 200 MPa, uY = 250 MPa, a, = 300 MPa, s = 0.00012 MPa, &pD = 0, D,, = 0.05.

loo_

Fig. 6. Elastoplasticity

and ductile damage in pure tension at microscale. = l/3, epn = 19.6%.

Fig. 7. Effect of the triaxiality upon the strain to rupture. epR is the strain to rupture in pure tension ~,/a,,

J. Lemaitre, I. Doghri I 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, elr is the main input and the two other principal strain components are imposed as .ez2= es3 = - ~ei,. 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 as = 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. Low 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 damage rate equation, it is valid even when a cycle cannot be defined. The material data are those of an

0 0

I 1

I 2

I 3

I 4

I 5

I 6

I ?

I 8

E,,a O-9

Fig. 9. Stress-strain

and brittle damage at mesoscale in pure tension.

210

Fig. 10. Stress-strain

and damage

at microscale

when

pure

tension

is applied

at mesoscale.

aluminum alloy: E = 72 000 MPa, v = 0.32, a, = 303 MPa, a,, = 306 MPa, aU = 500 MPa, s = 6 MPa, &pn = lo%, D,, = 0.99. As for the example of Section 5.1, pure tension cases at microscale are considered by giving as input data eli, the other strains being computed by DAMAGE 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: 20.425 us (MPa) 303 kO.43 305 +0.47 308 +1 370 +3.5 440 +4.5 460

The number 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, (p =p,) is also reported. As shown by experiments, the ratio N,lN, increases as NR increases. The details of the results corresponding to the case AE~~= 7% are given in Fig. 12, showing the evolution of the strain &II and the stress-strain loops (ull, E11) as a function of the number of cycles, the equivalent stress ueeqand the accumulated plastic strain p together with the damage D.

.,
I II

10

10

lo

Id

N
at microscale.

Fig. 11. Fatigue

rupture

curve

in pure tension

600%
400.00

200.00 0.M) I -200.00 -400.00

1:i.L

40.00

50.00

Fig. 12. Results of a very low cycle fatigue in pure tension at microscale AE~~ = 7%, NR = 40 cycles.

It is interesting to compare this with the case of a pure elastic tension case at mesoscale inducing a hree-dimensional state of stress at microscale. The same input are used as for the case of Fig. 12 except hat the strains applied are &ii = +3.5%, &22= &33= - mll. The four similar graphs are shown in Fig. -3.

An important feature of fatigue damage is the non-linearity of the accumulation of damages due to oadings of different amplitudes. If a two level loading history is considered with N1 cycles at the strain amplitude A&i and N2 cycles at the strain amplitude AF* with N1 + N2 = NR the number of cycles to upture, in contrast with the Palmgreen-Miner rule, generally

Nl

N2

NR1(Aq) + NRZ(A~2) < if &I> A&z N1 N2 NR,(Ael) + NRZ(A&J if A&l< A&z

15Ot20

1000.00

500.00

0.00

-500.00

-1000.00

-1500.00

IO

Em3

-2000.00 -0.06

-0.04

-0.02

0.04

Fig. 13. Results of a very low cycle fatigue in pure tension at mesoscale inducing three-dimensional he,, = 7%, NR = 8.

fatigue at microscale

Nal and Na, being the number of cycles to failure corresponding to the constant amplitude cases of fatigue respectively at As, and AQ. This effect was studied with the same material data as those in Example 5.3 for a one-dimensional case at mesoscale .sz2= .s33= -VEX1. +0.47% for the highest level , NR = 7720 cycles NR = 109 570 cycles . diagram of Fig. 14; they are in good qualitative agreement

-0.425% for the lowest level , The results are given in the accumulation with experimental results.

Two cases of multi-axial fatigue were studied: biaxial 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 DAMAGE 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.

213

N2 z,

.a

.6

Fig. 14. Accumulation

diagram for two level tests in tension at mesoscale.

A plane state of strain is considered

with the following in phase strains: 0. to NR = 7800 cycles.

=+-Y 9

E33 =

Fig. 15 shows contours of cycles to failure corresponding

This is the same type of calculation with

= +x )

El2

2Y9

all other components = 0.

Fig. 15. Biaxial fatigue in plane strain.

214

i!v
2 .4_ , .I_

.2_

.l_

0 0

I .l

I .2

I .3

I .I

I .s

I .a 5% 2

Fig. 16. Fatigue in tension and shear strains imposed.

Fig. 16 summarizes

the results for the same number of cycles to failure.

ill J. Lemaitre and D. Baptiste, On damage criteria, Proc. NSF Workshop on Mechanics of Damage and Fracture, Atlanta, GA (1982). PI F. McClintock, A criterion for ductile fracture by the growth of holes, ASME J. Appl. Mech. (1968). [31 J. Rice and D. Tracey. On ductile enlargement of voids in triaxial stress fields, J. Mech. Phys. Solids 17 (1969). [41 S.S. Manson, Behavior of materials under conditions of thermal stress, NACA Tech. Rep. 1170, 1954. [51 F.F. Coffin, Study on the effects of cyclic thermal stresses in a ductile metal, Trans. ASME (1954) 76-931. bl M.A. Miner, Cumulative damage in fatigue, J. Appl. Mech. 12 (1945) no. 3. I71 N.J. Hoff, J. Appl. Mech. 20 (1953) 105-108. PI D. Hayhurst, J. Phys. Solids 20 (1972) 381. I91 D. Krajcinovic and G.U. Fonseka, The continuous damage mechanics of brittle materials - Parts I and II, J. Appl. Mech. 48 (1981) 809. [lOI L.M. Kachanov, On the creep rupture, Izv. SSSR, Otd. Tekhn. Nauk 8 (1958) 26. [Ill Y.N. Rabotonov, On the equations of state for creep, in: Progress in Applied Mechanics, Prager Anniversary Volume (Macmillan, New York, 1963) 307. t1-21 J. Lemaitre, Evaluation of dissipation and damage in metals submitted to dynamic loading, Proc. ICMl, Kyoto, Japan (1971). D. Marquis and J. Lemaitre, Modelling complex behaviors of metals by the State Kinetic Coupling theory, Proc. ASME Winter Meeting, San Francisco (1989). J.L. Chaboche, Description thermodynamique et phenomenologique de la viscoplasticite avec endommagement, These de Doctorat es Sciences, Universite Paris 6, 1978. C. Lienard, Plasticite couplee a Iendommagement en conditions quasi-unilaterales pour la prevision de Iamoqage des fissures, These, Universite Paris 6, 1989. [I61 J. Lemaitre, Micro-mechanics of crack initiation, Int. J. Fracture 42 (1990) 87-99. [I71 A. Benallal, Thermoviscoplasticite et endommagement des structures, These de Doctorat es Sciences, Universite Paris 6, 1989. [I81 G.I. Taylor, Plastic strains in metals, J. Inst. Met. (1938) 62-307. 1191 J. Lemaitre and J.L. Chaboche, Mechanics of Solid Materials (Cambridge University Press, 1990). WI 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. WI J. Lemaitre and J. Dufailly, Damage measurements, Eng. Fract. Mech., 28 (1987) 643-661. WI 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. v31 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. (241 R. Billardon, Etude de la rupture par la mtcanique de Iendommagement, These de Doctorat es Sciences, Universite Paris 6, 1989. WI J.-P. Cordebois, Criteres dinstabilities plastiques et endommagement ductile en grandes deformations. Applications a Iemboutissage, These de Doctorat es Sciences, Universite Paris 6, 1983.

J. Lemaitre, I. Doghri I Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232

215

Give material constants and the strains history * * * DAMAGE90 * * * will give you the damage growth up to crack initiation Elasticity. Give: Youngs modulus, Poissons 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: dD = (Y/S)dp. Give: S 6. Damage threshold dD = 0 if p < pD. Do you know pD? Gi:e EpD in tension: pD = EpD * (SIGu - SIGf)/ [SIGs - (SIGf * * 2/SIGy)] 10.E - 2 Crack initiation: D = DC. Do you know DC? n Give Dlc in tension: DC = Dlc* [(SIGu/SIGs) * *2]/Rnu 0.99 Initial conditions. Give: p0, DO o., 0. Is the stress state uniaxial? y or n n Is the strain history cyclic? y or n Y * ** YOUR LOADING 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 11-strain 0.425E - 2, -0.425E - 2 Give: 1st peak, 2nd peak of 22-strain -0.136E - 2, 0.136E - 2 Give: 1st peak, 2nd peak of 33-strain -0.136E - 2, 0.136E - 2 1Zstrain Give: 1st peak, 2nd peak of o., 0. Give: 1st peak, 2nd peak of 13-strain o., 0. 23-strain Give: 1st peak, 2nd peak of o., 0. Give the value of the plastic threshold stress SIGs for this block 303. Suggest a number of increments per cycle (minimum: 4) 4 * ** CRACK INITIATION. * ** THE JOB IS ENDED. YOUR RESULTS FILES ARE:

216

J.

[Link] STOP:

[Link]

[Link]

CYCLE [Link] +OO, 0.2.500E +OO, 0.5OOOE+OO, 0.7500E+OO, [Link] +Ol , o.l858E+04, O.l858E+04, 0.2785E+04, 0.2786E+04, 0.2786E+04, 0.2786E304, 0.3713E+04, 0.3714E+04, 0.3714E+04, 0.3714E+04, 0.4641E +04, 0.4642E+04, 0.4642E +04, 0.4642E+04, O.l086E+06, O.l086E+06, O.l086E+06, 0.1086E +06, 0.1086E +06,, O.l086E+06, O.l086E+06, O.l086E+06, O.l096E+06,

D [Link]+OO, [Link] + 00, [Link] + 00, [Link]+OO, [Link]+OO, [Link]+OO, [Link]+OO, [Link]+OO, [Link]+OO, [Link] + 00, [Link]+OO, [Link]+OO, [Link]+OO, [Link]+OO, [Link]+OO, [Link]+OO, [Link]+OO, [Link]+OO, 0.0000E+00, 0.9858E300, 0.9858E +OO, 0.9858E +OO, 0.9858E +OO, 0.9858E+OO, 0.9858E +OO, 0.9858E +OO, 0.9858E+OO,
[Link]+Ol ,

~.OOOOE+OO 0.3667E-04; 0.3667E-04, [Link]-03, [Link] -03, 0.2725E+OO, 0.2725E+OO, 0.4085E+OO, 0.4085E+OO, 0.4086E+OO, 0.4086E+OO, 0.5446E+OO, 0.5446E+OO, 0.5447E+OO, 0.5447E+OO, 0.6807E+OO, 0.6807E+OO, 0.6808E+00, 0.6808E+00, O.l593E+02, O.l593E+02, O.l593E+02, O.l593E+02, 0,1593E+02, 0.1593E302, O.l593E+02, O.l593E+02, 0.1607E +02,

MISESSTRESS [Link]+OO, 0.3030E+03, 0.3OOOE+Ol, 0.3030E+03, 0.3000E+Ol, 0.3030E +03, 0.3000E+Ol, 0.3030E +03, 0.3OOOE+Ol, 0.3030E-tO3, 0.3OOOE+Ol, 0.3030E-tO3, 0.3000E+Ol, 0.3030E-tO3, 0.30OOE +Ol, 0.3030E+03, 0.3000E+Ol, 0.3030E+03, 0.3000E+Ol, 0.4298E+Ol, 0.4255E -01, 0.4298E+Ol, 0.4255E -01, 0.4296E+Ol, 0.4253E -01, 0.4293E+Ol, 0.4251E-01, -0.9480E-01,

[Link] [Link]+OO, 0.3034E+03, 0.2814E+Ol, 0.3034E+03, 0.2814E+Ol, 0.3034E+03, 0.2814E+Ol, 0.3034E+03, 0.2814E +Ol, 0.3034E303, 0.2814E+Ol, 0.3034E303, 0.2814E +Ol, 0.3034E+03, 0.2814E+Ol, 0.3034E+03, 0.2814E+Ol, 0.3034E+03, 0.2814E+Ol, 0.4303E+Ol, 0.3992E -01, 0.4303E +Ol, 0.3992E-01, 0.4301E +Ol, 0.3990E -01, 0.4298E+Ol, 0.3987E-01, -01, -,0.9492E

Appendix B
__________________-_----~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~----------C* ____________________--~~-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~----------X _____-

C*** DAMAGE SO Elastic-perfectly plastic law coupled to a ductile damage model. C WI* Fully implicit integrationscheme. C*2* Jump in cycles procedure. C*** Written by Issam Doghri C*** Version: May 1990 c*====================================================================x====== implicit none integer ntens, nstatv, ms, mb, naxi, ncycle, nblock, ibl integer ist, ngiv, ig, ilook0, ipass, nul, nu2, igiv, icycle integer nrcycl, i, ncysum, ideltn, pslope, ncycl0, kk, j real*8 e0, xnu, sigs, smpd, SO, one, two, three, dcrit real*8 sigf, sigy, sigu, epd, dlc, smp0, do, energO, energb real*8 tolene, rac2, dtime, energl, time, tfin, strsb, star real*8 mu, energ2, dtmin, tper, dhtime, smpi, dhi, rtime real*8 xxnu, ehl, epsmin, smpf, dhf, ddmax0, yyy, trdeps, yyl real*8 yy2, dpmax, dpmax0, deltp, ddmax, dslope, deltd CHARACTER*1UNIAXI,ANSPD,ANSDC,CYCLIC,JUMP,CONVER,STAB,COUPL CHARACTER*20FILEI,FILE2 CHARACTER*6COMMENT PARAMETER(NTENS=G,NSTATV=ll) REAL*8 STRS(NTENS),STATEV(NSTATV),STATEVI(NSTATV), . STRAN(NTENS),DSTRAN(NTENS),SIGHI(6),SIGHF(6), . INCUSE(50>,TIM(200>,HIST(6,2OO),EPSB(6,5O),EPSS(6,5O), PERDIV(G),SYIEL(50) 'INTEGER ISLOPE(6),IPER(6),INDEX(G),NBCYCL(50) COMMON/ETIql/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=l.D+lO DCRIT=[Link] WRITE(MS,*)' WRITE(MS,*)'
WRITE(MS,*) **********************************************************

WRITE(MS,*)' Give material constants and the strain history' WRITE(MS,*)' WRITE(MS,*)' *** DAMAGE SO ***' WRITE(MS.*)' WRITE(MS,*)'will give you the damage growth up to crack initiation' NRITE(MS,*)'**********************************************************' WRITE(MS,*) WRITE(MS,*)'**ELASTICITY . Give :' WRITE(MS,*)' YOUNG"s modulus :'

READ(MB,*)EO WRITE(MS,*)' POISSON"s ratio : READ(#B,*)XNU WRITE(MS,*)'**PERFECT PLASTICITY : plastic threshold SIGs . given with loading . Give :' WRITE(IIS,*)' Fatigue limit SIGf : READ(MB,*)SIGF WRITE(MS,*)' Yield stress SIGy :' READ(HB,*)SIGY WRITE(MS,*)' Ultimate stress SIGu :' READ(MB,*)SIGU WRITE(MS,*)'**DAHAGE EVOLUTION : dD = (Y/S) dp. Give S' READ(MB,*)SO FORMAT(A) h!RITE(MS,*)'** DAMAGE THRESHOLD : dD=O if p<pD. Do you know
. pD? JaYJ) or NpJJ

READ(MB,I)ANSPD IF(([Link].'y').OR.([Link].'Y'))TIiEN WRITE(HS,*)' Give the value of pD :' READ(MB,*)SMPD IF([Link].O.)SMPD=l.D-6 ELSE IF(([Link].'n').OR.([Link].'N'))TAEN WRITE(MS,I) 2' Give EpD In Tension : 3 pD=EpD*(SIGu-SIGf)/[SIGs-(SIGf**2/SIGy)]' READ(MB,*)EPD IF([Link].O.)EPD=i.D-6 ELSE STOP'*** WRONG DATA. GOOD BYE.' END IF WRITE(MS,*)'**CRACK INITIATION : D=Dc . Do you know DC ? "Y" READ(MB,l)ANSDC IF(([Link].'y').OR.([Link].'Y'))THEN WRITE(MS,*)' Give the value of DC (remember:O<Dc<i) :' READ(WB,*)DCRIT ELSEIF(([Link].'nJ).0R.(ANSDC.E9.'1'))THEN WRITE(MS,*)' Give Dlc in tension : . DC = Dlc * [(SIGu/SIGs)**2] / Rnu READ(MB,*)Dlc ELSE STOP'*** WRONG DATA. GOOD BYE.' END IF WRITE(MS,*)'**INITIAL CONDITIONS . Give :' WRITE(MS,*)' The value of PO :' READ(KB,*)SMPO WRITE(MS,*)' The value of Do : READ(MB,*)DO WRITE(HS,*)'**LOADING' WRITE(MS,*)' Is the stress state uniaxial? "Y" or READ(MB,I)UNIAXI IF([Link].'y'.[Link].'n'.AND.

or "N"'

[Link].'Y'.[Link].'N')STOP'*** WRONG DATA. GOOD BYE.***' 'IF(("[Link]&'y').OR.([Link].'Y')) NAXI=l URITE(HS,*)' Is the strain history cyclic? "Y" or "I"' READ(MB,I)CYCLIC IF([Link].'y'.[Link].'n'.AND. [Link].'Y'.[Link].'N')STOP'***WRONG DATA. . . GOOD BYE.***' IF(([Link].'y').OR.([Link].'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'J' READ(MB,I)JUMP IF([Link].'y'.[Link].'n'.AND. [Link].'Y'.[Link].'N')STOP'*** WRONG DATA. . GOOD BYE.***' WRITE(MS,*)' Give the number of blocks of constant amplitude . (50 max)' READ(MB,*)NBLOCK DO IBL=l,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) 3 FORMAT('Give:1st 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 b SIGs for this block : READ(MB,*)SYIEL(IBL) WRITE(MS,*)' Suggest a number of incrementsper cycle . (minimum : 4) : ) READ(MB,*)INCUSE(IBL) END DO ELSE IF(([Link].'n').OR.([Link].'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=I,NAXI WRITE(MS,2)INDEX(IST) 2 FORMAT0 Give the values of ',12,' -strain at these' . times :'> READ(MB,*)(HIST(IST,IG),IG=i,NGIV)

END DO WRITE(MS,*)' Give the values of the plastic threshold stress S &IGs at these times : READ(MB,*)(SYIEL(IG),IG=l,NGIV) WRITE(MS,*)' Suggest an initial time increment < to Lintervallebetueen first 2 points, READ(MB,*)DTIME END IF C***[Link]= r,p,D,dD,(pstrn>,dp DO IST=l,NSTATV STATEVI(IST)=O. END DO STATEVI(2)=SHPO STATEVI(3)=DO STATEVI(4)=0. DO IST=l,6 STRAN(IST)=O. DSTRAN(IST)=O. END DO IF(([Link].'y').OR.([Link].'Y'))EPD=SMPD ENERGl=(SIGU-SIGF)*EPD ENERGO=O. ENERGS=O. TOLENE=[Link]-2 ILOOKO=O COUPL='y' IPASS=O ONE=l.D+OO TWO=2.D+00 THREE=3.D+00 RACZ=DSQRT(TWO) NUl=ll NU2=12 FILEl=,[Link]' FILEZ=,[Link], OPEN(UNIT=NUl,file=FILEl) OPEN(UNIT=NU2,file=FILEZ) COMMENT=' TIME' IF(([Link].'y').OR.([Link].'Y'))COMMENT=' CYCLE' WRITE(NUl,312)COMMENT URITE(NU2,312)COMMENT WRITE(90,3lS)COMMENT 312 FORMAT(2X,A6,2lX,'STRAINS',27X,'STRESSES',l6X,'DAMAGE', . 9X,'p',8X, 'MISESJ,2X,'[Link]') WRITE(NUl,313) WRITE(NU2,314) 313 FORMAT(l9X,'li',9X,'22',9X,'33',lOX,'ll',9X,'22',9X,'33', llX,,D,,23X,,SIG eq,,SX,,SIG*,) 314 FORnAT(l9X,,l2',9X,'l3',9X,'23',lOX,'l2',9X,'i3',9X,'23', llX,,D,,23X,'SIG eq',SX,'SIG*,) 315 FORMAT(SX,A6,lOX,'DAMAGE',8X,'p',SX,'MISES [Link]'

221 [Link]') C*** IF TEE LOADING IS NOT CYCLIC IF(([Link].'n').or.([Link].'N'))THEN JUMP='n' TIME=TIM(l) DTMIN=DTIME/l0000. TFIN=TIM(NGIV) IGIV=I CALL OUTPUT(NUI,NU2,TIME,STRAN,STRS,STATEV,NTENS,NSTATV, . STRSB,STAR) DOUHILE(TIME+[Link]) DO IST=I,NAXI DSTRAN(IST)=DTIME*(HIST(IST,IGIV+1)-HIST(IST,IGIV~~ & /(TI~(IGIV+I)-TIM(IGIV)) END DO CONVER='y' DO IST=I,NSTATV STATEV(IST)=STATEVI(IST) END DO SIGS=SYIEL(IGIV) IF([Link])STOP'*** PLEASE REVIEW THE & 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([Link].'y'.[Link](ll).GT.O.)THEN ENERG2=ENERG3+(SIGS-SIGF*(SIGF/SIGY))*STATEV(ll) IF([Link].'n')DCRIT=DlC*((SIGU/SICS)**2~/RNU IF(ILOOKO.E~.[Link].(ANSPD.E~.'y'.[Link].E~.'Y') .[Link](2).[Link]*[Link]. ILOOKO.E~.[Link].(ANSPD.E~.'n'.[Link].E~.'N'~.A~. [Link].I.05*[Link]. STATEV(3).GT.l.O5*DCRIT)THEN CONVER='n' END IF END IF IF([Link].'n')THEN IPASS=IPASS+l DTIME=DTIME/2. C*** IF CONVERGENCE ELSE IF(([Link].'n'.[Link].'N') .[Link].0.99DO)DCRIT=O.99DO IPASS=O DO IST=l,NSTATV STATEVI(IST)=STATEV(IST) END DO DO IST=I,6 STRAN(IST)=STRAN(IST)+DSTRAN(IST) END DO TIME=TIME+DTIME

CALL OUTPUT(NUl,NU2,TIHE,STRAN,STRS,STATEV,NTENS,NSTATV, . STRSB,STAR) IF(STATEV(3).[Link])THEN WRITE(HS,*)' ' WRITE(HS,*)'***CRACK INITIATION.' GO TO 308 END IF IF(STATEV(3).GE.l.l)STOP'*** DAMAGE CANNOT EXCEED I.' C*** FIND pD IF(ILOOKO.E~.[Link](ii).GT.O.)THEI ENERG3=ENERG2 IF(([Link].'n'.[Link].'N') .[Link])TIiEN SMPD=STATEV(2) ILOOKO=I ELSE IF(([Link].'y').OR.([Link].'Y').AND. STATEV(2).[Link]>THEN ILOOKO=I END IF END IF IF(DABS(TIME-TIM(IGIV+I)).[Link])THEN IGIV=IGIV+l END IF END IF DTIME=DTIME*l.I IF(TIME+[Link](IGIVtl)[Link]~[Link])THEN DTIME=TIM(IGIV+I)-TIME END IF IF([Link])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(([Link].'y').OR.([Link].'Y'))THEN TPER=l. DTMIN=TPER/40000. TIME=O. ICYCLE=I NRCYCL=ICYCLE DHTIME=O. DO 1=1,6 SIGHI(I)=O. END DO SMPI=SMPO DHI=DO IF(([Link].'y').OR.([Link].'Y'))THEN COUPL='n' OPEN (UNIT=80,file=,[Link],)

D P' TIME WRITE(80,*)'REAL TIME WRITE(80,*)'******************************' WRITE(80,*)'BEGINNING OF THE CYCLE' =',URCYCL WRITE(80,*)'CYCLE(REAL) WRITE(80,*)'CYCLE(MACHINE) =',ICYCLE WRITE(80,400)RTIME,TIME,DHI,SMPI END IF CALL OUTPUT(NUi,NU2,TIME,STRAN,STRS,STATEV,NTENS,NSTATV, . STRSB'STAR) IBL=l DTIME=TPER/INCUSE(IBL) DO IST=l,6 IF(EPSB(IST,IBL)*EPSS(IST,IBL).LT.O.)THEN PERDIV(IST)=TPER/4. ELSE PERDIV(IST)=TPER/2. END IF IPER(IST)=l ISLOPE(IST)=l END DO NCYSUM=O DOWHILE([Link]) DO IST=I,NAXI IF(TIME+[Link](IST))THEN DSTRAN(IST)=EPSB(IST,IBL)*DTIME/PERDIV(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=l,NSTATV STATEV(IST)=STATEVI(IST) END DO SIGS=SYIEL(IBL) EHl=SIGS/lOOO. EPSMIN=SIGS/E0/10000. IF([Link])STOP'*** PLEASE REVIEW THE & VALUES OF SIGf SIGu and SIGs . CALL INTEGR(STRAN,DSTRAN,NTENS,NSTATV,[Link], b CONVER,STRS,STATEV,STRSB,STAR,RNU) C*** IF NO CONVERGENCE,DIVIDE TIME INC. BY 2 IF(([Link].'n'.[Link].'U').AND. STATEV(Il).GT.O.) DCRIT=DlC*((SIGU/SIGS)**2)/RNU IF([Link].'n')THEN IPASS=IPASS+l DTIME=DTIME/2. C*** IF CONVERGENCE ELSE IF(([Link].'n'.[Link].'W')

224

.[Link].0.99DO)DCRIT=O.99DO IPASS=O DO IST=l,NSTATV STATEVI(IST)=STATEV(IST) END DO DO IST=1,6 STRAN(IST)=STRAN(IST)+DSTRAN(IST) END DO TIME=TIME+DTIME RTIME=TIHE+DHTIME CALL OUTPUT(NUI,NU2,RTIn,STRAI,STRS,STATEV,NTENS,NSTATV, L STRSB,STAR) IF(STATEV(3).[Link])THEN WRITE(HS,*)' ' WRITE(MS.*)'***CRACK INITIATION.' GO TO 308 END IF IF(STATEV(3).GE.l.l)STOPJ*** DAMAGE CANNOT EXCEED 1.' DO IST=I,6 IF(DABS(TIME-IPER(IST)*PERDIV(IST)).[Link])THEN IPER(IST)=IPER(IST)+l DTIME=I./INCUSE(IBL) END IF IF(ISLOPE(IST).[Link](STRAN(IST)-EPSB(IST,IBL)).LE. . [Link](IST).[Link](STRN(IST)-EPSS(IST,IBL)) . .[Link])THEN ISLOPE(IST)=-ISLOPE(IST) END IF END DO C*** END OF A CYCLE IF(DABS(TIME-TPER*ICYCLE).[Link])THEN IDELTN=O DO 1=1,6 SIGHF(I)=STRS(I) END DO SMPF=STATEV(2) DBF=STATEV(3) DELTP=SMPF-SMPI DELTD=DBF-DHI C*** FIND pD IF([Link].O)THEN ENERG3=ENERGO+(SIGS-SIGF*(SIGF/SIGY))* . DELTP*(NRCYCL-NCYSUM) IF(([Link].'n'.[Link].'N') .[Link])THEN SMPD=SMPF NCYCLO=NRCYCL-NCYSUM ILOOKO=l ELSE IF(([Link].'y'.[Link].'Y').AND. STATEV(2).[Link])THEN NCYCLO=NRCYCL-NCYSUM

&

ILOOKO=l END IF END IF C*** JUMP IN CYCLES PROCEDURE IF(([Link].'y').OR.([Link].'Y'))THEB WRITE(80,*)'ENDOF THE CYCLE' WRITE(80,400)RTIME,TIHE,DHF,SMPF DDMAXO=DCRIT/SO. IF(([Link].'n'>.OR.([Link].'N')) DDMAXO=DlC/50. IF(([Link].'y').OR.([Link].'Y'))THEN YYY=SIGS*SIGS/2./EO ELSE TRDEPS=O. DO IST=1,3 TRDEPS=TRDEPS+EPSB(IST,IBL)-EPSS(IST,IBL) END DO YYl=(I.+XNU)*SIGS*SIGS/3./EO YY2=EO*TRDEPS*TRDEPS/(l.-2.*XNU)/6. YYY=YYl+YYZ END IF DPMAXO=SO*DDMAXO/YYY DPMAX=DPMAXO PSLOPE=NBCYCL(IBL)*DELTP IF([Link])DPMAX=PSLOPE WRITE(80,*)'ILOOKO=',ILOOKO,'NCYCLO=',NCYCLO,'pD=',SMPD IF([Link].'[Link](2).[Link])THEN STAB='y' DO KK=I,6 IF(DABS(SIGHF(KK)-SIGHI(KK)).[Link])STAB='n' END DO WRITE(80,*)'ENERGl=',ENERGl,'ENERG3=',ENERG3 IF([Link].'y') THEN WRITE(80,*)'STABILIZED CYCLE' IF(STATEV(2>.[Link])THEN C*** JUMP OF CYCLES BEFORE DAMAGE GROWTH IDELTN=IDINT(DPMAX/DELTP) IF(NRCYCL+[Link]+NBCYCL(IBL)) IDELTN=NCYSUM+NBCYCL(IBL)-NRCYCL ENERG2=ENERGO+(SIGS-SIGY)*DELTP*(NRCYCL+IDELTN-NCYS~) IF([Link].ENERG1*1.05) IDELTN=(ENERGl-ENERG3)/(SIGS-SIGY)/DELTP WRITE(80,*)'***JUMP OF ',IDELTN,'CYCLES' WRITE(80,*)'DPMAX=',IDELTN*DELTP, 'YYY=',YYY ELSE WRITE(80,*)'COUPLED COMPUTATIONFOR NEXT CYCLE COUPL='y' END IF END IF ELSE C*** JUMP OF CYCLES AFTER DAMAGE GROWTH

DDMAX=DDHAXO DSLOPE=(NBCYCL(IBL)-NCYCLO)*DELTD IF([Link])DDMAX=DSLOPE NCYCLO=O IDELTN=IDINT(DDMAX/DELTD) IF(NRCYCL+[Link](IBL)) . IDELTN=NCYSUMtNBCYCL(IBL)-NRCYCL IF(DBFtDELTD*[Link]*1.O5) . IDELTN=IDINT((DCRIT-DHF)/DELTD) IF([Link](DPMAX/DELTP))IDELTN=IDINT(DPMAX/DELTP) WRITE(80,*)'***JUMP OF ',IDELTN,'CYCLES' WRITE(80,*)'DDHAX=',IDELTN*DELTD, 'DPMAX=',IDELTN*DELTP COUPL=,n, END IF IF([Link])WRITE(80,*)'*** PROBLEM : DPMAXO IS TOO SMALL.' IF([Link])YRITE(80,*)'*** PROBLEM : . DDMAXO IS TOO SMALL.' IF([Link])IDELTN=O IF([Link].O.)STOP'*** PROBLEM : NEGATIVE JUMP OF CYCLES.' END IF DO J=1,6 SIGHI(J)=SIGHF(J) END DO SMPI=SMPFtDELTP*IDELTN DHI=DHFtDELTD*IDELTN STATEVI(2)=SMPI STATEVI(3)=DHI DHTIME=DHTIMEtTPER*IDELTN NRCYCL=NRCYCLtIDELTN RTIME=TIMEtDHTIHE IF(([Link].'y').OR.([Link].'Y'))THEN WRITE(80,*)'-------------------------' WRITE(80,*)'BEGINNING OF THE CYCLE' WRITE(80,*)'CYCLE(REAL) =',NRCYCLtI WRITE(80,*)'CYCLE(MACHINE) =',ICYCLEtl WRITE(80,400)RTIME,TIME,DHI,SMPI END IF IF([Link](IBL))THEN ENERGO=ENERGO+(SIGS-SIGF*(SIGF/SIGY))*DELTP*(NRCYCL-NCYSUM) NCYSUM=NCYSUMtNBCYCL(IBL) IBL=IBL+l END IF ICYCLE=ICYCLE+l NRCYCL=NRCYCLtl END IF FORMAT(IX,4(El2.6,',')) END IF DO IST=1,6 IF([Link](IST)*PERDIV(IST)tDTMIN)THEN DTIME=(IPER(IST)*PERDIV(IST))-TIME

400

.I. Lemaitre, I. Doghri I Comput. Methods Appl. Mech. Engrg. 115 (1994) 197-232

221

END IF END DO IF([Link])THEN WRITE(HS,*)'***NO CONVERGENCE' WRITE(MS,*)'DTIME=',DTIME,'IPASS=',IPASS GO TO 308 END IF END DO END IF WRITE(MS,*)'***THE JOB IS ENDED. YOUR RESULT FILES ARE :' 308 WRITE(MS,*)' IF(([Link].'n'>.OR.([Link].'N'))WRITE(MS,*~FILEl,FILE2 IF(([Link]).'y').OR.([Link].'Y')) . WRITE(MS,*)FILE1,FILE2,'[Link] STOP END C*____________________________________________________________________X______ _______------_______-------_~~~~~~~~~-~--~~~~~~~~~~~~~~~~~~~~~ ______ SUBROUTINE INTEGR(STRAN,DSTRAN,NTENS,NSTATV,UNIAXI,COUPL, % CONVER,[Link],STRSB,STAR,RNU) C*____________________________________________________________________ ____________________~~~~~~~~~~~~~~~~~~~~~~~_~~~~~~~~~~~~~~~~~~~~~~~~ ______ X_____implicit none integer ntens, nstatv, niss, i, j, nfois, ntest real*8 e0, xnu, sigs, smpd, SO, one, two, three, strsb, star real*8 mu, 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), % STRN(S),ID(S,S),AS(S>, % STRS(G),ESTRS(G), CESTRS(6),KESTRS(6), % N(S),PSTRN(S),DPSTRN(6) REAL*8 LAMDA,NN LOGICAL CRITERE CHARACTER*1CONVER,COUPL,UNIAXI COMMON/ETIql/EO,XNU,SIGS,SMPD,SO COMMON/ETIq2/0NE,TWO,TIiREE C NN=THREE/TWO RAC2=DSqRT(TWO) TOLRES=SIGS*l.D-6 NISS=NTENS C*** IDENTITY MATRIX DO I=l,NISS-1 DO J=I+l,NISS ID(I,J)=O. ID(J,l)=O. END DO ID(I,I)=ONE END DO ID(NISS,NISS)=ONE C*** IDENTITY VECTOR DO I=l,NISS

IF([Link].3)AS(I)=ONE IF([Link].4)AS(I)=O. END DO C*** LAME COEFF. XMU=EO/(ONE+XNU)/TUO LAMDA=EO*XNU/(ONE-TUO*XNU)/(ONE+XNU)

SMR=STATEV(l) SMP=STATEV(2) D=STATEV(3) DO I=l,NISS PSTRN(I)=STATEV(I+4) END DO C*** TOTAL STRAIN DO I=I,NISS STRN(I)=STRAN(I)+DSTRAN(I) IF([Link].4)STRN(I>=STRN(I)*RAC2 END DO C*** ELASTIC PREDICTOR DSMP=O. IF(([Link].'y').OR.([Link].'Y'))THEN STRN(2)=-XNU*STRN(l>-(0.5-XNU)*PSTRN(l) STRN(3)=STRN(2) END IF C . . . . Trace of strain tensor TREPS=STRN(l)+STRN(2>+STRN(3) C . . . . Effective stress DO I=l,NISS ESTRS(I)=LAMDA*TREPS*AS(I)+2.*XMU*(ST~(I)-PSTRN(I)) END DO CALL INVAR(ESTRS,ESTRSH,ESTRSB,NISS) SIGN=l. IF(ESTRS(l>.LT.O.)SIGN=-1. C*** TEST THE YIELD CONDITION F=ESTRSB-SIGS IF([Link].O.)THEN GO TO 500 ELSE C*** IF THE ELASTIC PREDICTOR DOES NOT VERIFY THE YIELD CONDITION DO I=l,NISS N(I>=NN*(ESTRS(I)-ESTRSH*AS(I))/ESTRSB END DO C*** Correctionsover the elastic predictor c....CorrectionOVER P DENOM=3.*XMU CSMP=F/DENOM C . . . . Correction over effective stress DO I=l,NISS CESTRS(I)=-2.*F*N(I)/3. END DO C*** LOOP ON THE PLASTIC CORRECTIONS

NFOIS=O CRITERE=.TRUE. DO WHILE(CRITERE) NFOIS=NFOIS+l DIVl=ESTRSB NTEST=2 IF([Link].l)NTEST=l C*** UPDATE THE STATE C . . . . Update inc. of p DSMP=DSMP+CSMP C . . . . Update effective stress DO I=l,NISS ESTRS(I)=ESTRS(I)+CESTRS(I) END DO CALL INVAR(ESTRS,ESTRSH,ESTRSB,NISS) DO I=l,NISS N(I>=NN*(ESTRS(I)-ESTRSIi*AS(I))/ESTRSB END DO IF([Link].2)THEN DIVZ=ESTRSB REL=(DIVZ-DIVI)/DIVI END IF main C . . . . If we converge too slowly,or we [Link] C . . . . program will propose a smaller "time" increment. IF([Link].50.0R. & ([Link].I5.D-02))THEN CONVER='n' GO TO 500 END IF C*** Compute residual functions F=ESTRSB-SIGS IF(([Link].'y').OR.([Link].'Y'))THEN STRN(2)=-O.5*STRN(1)+(O.5-XNU)*SIGN*SIGS/EO STRN(3)=STRN(2) TREPS=STRN(l)+STRN(2)+STRN(3) END IF DO I=i,NISS KESTRS(I)=ESTRS(I)-LAHDA*TREPS*AS(I)& 2.*XMU*(STRN(I)-PSTRN(I)-DSMP*N(I)) END DO C . . . . Max. residuals RESTRS=O. DO I=I,NISS IF(DABS(KESTRS(I)).[Link])RESTRS=DABS(KESTRS(I)) END DO C*** ITERATIVE TEST ON THE YIELD CONDITION IF(DABS(F).[Link])THEN CRITERE=.FALSE. ELSE C*** PLASTIC CORRECTIONS

CALL VTRANVl(N,KESTRS,XNKES,NISS) C . . . . Correction over p DENOM=3.*XMU CSHP=(F-XNKES)/DENOM C . . . . Correction over effective stress DENOM=DENOM*DSMP/ESTRSB+(I.) DO I=I,NISS CESTRS(I)=(2./3.)*(XNKES-DENOM*F)*N(I)-KESTRS(1) 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([Link].'y')THEN IF(([Link].'y'>.OR.([Link].'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)*DSHP END DO C .... Damage increment DD=O. CALL DAMAGE(TREPS,ESTRSB,Y,RNU) IF([Link].O..[Link]+[Link].E~.'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(l)=SMR+DSMR STATEV(2)=SMP+DSMP STATEV(3)=D+DD STATEV(4)=DD DO I=l,NISS STATEV(I+4)=PSTRN(I)+DPSTRN(I) END DO STATEV(ll)=DSMP END .IF RETURN

231 END C*__-_______________________________________________________--________ ____________________________________________________________________X====== SUBROUTINE INVAR(V,VH,VB,NISS) C*===_____________________________________________-_________------____ ____-______________________________________________-__________-__X====== C ... 1st and 2nd stress invariant8 implicit none integer niss, i real*8 one, two, three, vh, vb, const REAL*8 V(6) COMMON/ETIq2/0NE,TWO,THREE C VH=(V(I)+V(2)+V(3))/THREE CONST=O. DO I=I,NISS IF([Link].3)CONST=CONST+(V(I)-VH)*(V(I)-VH) IF([Link].4>CONST=CONST+V(I)*V(I) END DO VB=DSqRT(THREE*CONST/TWO) RETURN END C*________-___________________________--_______---______----_____----_____--____----_____~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~--~~~~ X___--SUBROUTINEDAMAGE(TREPS,ESTRSB,Y,RNU) C*~~________________________--______----____------______----____-----_____--___----______~-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-~~~ ____-X___--implicit none real*8 e0, xnu, sigs, smpd, SO, treps, estrsb, y, mu, term1 real*8 term2 COMMON/ETI~l/EO,XNU,SIGS,SMPD,SO C TERMl=ESTRSB*ESTRSB*(I.+XNU)/EO TERM2=EO*TREPS*TREPS/(l.-2.*XNU)/2. Y=(TERMl+TERM2)/3. RNU=2.*EO*Y/ESTRSB/ESTRSB RETURN END C*_________________________---_____----______----______---_-_____---__ _____--____--_______~~~~~~~~~~~~~~~~~~~~~~-~~~~~~~~~~-~~~~~~~~--~~~~ X____-SUBROUTINE VTRANVl(V,Vl,VTVI,NISS) C*___--___----______-_____________________________________________________--____---__________________________~--~~~~~~~~~~-~~~~~~~~--~~~~X~~~~~~ C . . . Inner product of 2 symmetric 2nd order tensors implicit none integer niss, i real*8 vtvl REAL*8 V(S),Vi(S> C VTVl=O. DO I=[Link] VTVl=VTVI+V(I)*Vl(I) END DO RETURN END C*~~_________________________-______----______----________--__________ _____--___---___________________________~~~~~~~~~~~~~~~~~~~~~~~~~~X~~~~~~ SUBROUTINE OUTPUT(NUl,NU2,TIME,STRAN,STRS,STATEV,

232

NTENS,NSTATV,STRSB,STAR) implicit none integer nul, nu2, ntens, nstatv, i real*8 time, strsb, star REAL*8 STRAN(NTENS),STRS(NTENS),STATEV(NSTATV)
C

WRITE(NU1,130)TIME,(STRAN(I~,I=I,3~,~STRS~I~,I=1,3~, b STATEV(3),STATEV(2),STRSB,STAR WRITE(NU2,130)TIME,(STRAN(I),I=4,6),(STRS(I),I=4,6), & STATEV(3),STATEV(2),STRSB,STAR WRITE(90,444)TIME,STATEV(3),STATEV(2),STRSB,STAR 130 FORMAT(3X,ll(ElO.4,',')) FORMAT(3X,5(ElO.4,',')) 444 RETURN END