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