1 : C FATI SOURCE CHAT 07/12/04 21:15:55 5985 2 : C FATI 3 : SUBROUTINE FATI (WRK0,WRK1,WRKK2,WRK5,NSTRS,NVARI, 4 : 1 NMATT,ISTEP,ICARA,JDIM,IFOUR2,DT) 5 : C MAZARS SOURCE AM 00/12/13 21:39:37 4045 6 : C SUBROUTINE MAZARS (WRK0,WRK1,WRKK2,WRK5,NSTRS,NVARI, 7 : C 1 NMATT,ISTEP,ICARA,JDIM,IFOUR2) 8 : *====================================================================== 9 : * DB 9 : * FOR THE NONLOCAL CASE CALLED by ECOU40 --> FATTT --> FATI 10 : * New source: one nonlocal damage material model added: 11 : * Fatigue damage 12 : * This routine is inspired based on mazars.eso 13 : * 14 : C Goal: allowing cyclic cumulative damage calculation under simple fatigue loading 15 : C 16 : C Description: 17 : C 1. The threshold is maintiened allowing an endurance limit but not tested 18 : C 2. The equivalent strain is changed replaced by the geometric average of the positive stresses 19 : C induced strain ( no damage growth due to compression induced tension - Poisson effect) 20 : C 3. The damage rate given in [1] is integrated over a cycle with 21 : C dD/dN = f(D) eps_a^(beta+1)/(beta+1) 22 : C 4. Two expression of the f(D) funcion are implemented 23 : C a) LOI = 2. power law as chosen in [2] 24 : C b) LOI = 3. power law including a damage rate acceleration as describe in [1] 25 : C 26 : C References: 27 : C [1] D. Bodin, G. Pijaudier-Cabot, C. de La Roche, J.-M Piau and A. Chabot, (2004) 28 : C A Continuum Damage Approach to Asphalt Concrete Fatigue Modelling, Journal of 29 : C Engineering Mechanics, ASCE, vol. 130 (6), pp. 700-708. 30 : C 31 : C [2] Paas, R. H. J. W., Scheurs, P. J. G., and Brekelmans, W. A. M. (1993). ‘A 32 : C continuum approach to brittle and fatigue damage: Theory and numerical 33 : C procedures.’ Int. J. Solids Struct., 30~4!, 579–599. 34 : C 35 : * DB 36 : * 37 : *====================================================================== 38 : C 39 : C 40 : C 41 : C 42 : C variables en entree 43 : C 44 : C 45 : C WRK0 pointeur sur un segment deformation au pas precedent 46 : C 47 : C WRK1 pointeur sur un segment increment de deformation 48 : C 49 : C WRKK2 pointeur sur un segment variables internes au pas precedent 50 : C 51 : C WRK5 pointeur sur un segment de deformations inelastiques 52 : C 53 : C XMATER constantes du materiau 54 : C 55 : C NSTRS nombre de composantes dans les vecteurs des contraintes 56 : C et les vecteurs des deformations 57 : C 58 : C NVARI nombre de variables internes (doit etre egal a 2) 59 : C 60 : C NMATT nombre de constantes du materiau 61 : C 62 : C ISTEP flag utilise pour separer les etapes dans un calcul non local 63 : C ISTEP=0 -----> calcul local 64 : C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils 65 : C ISTEP=2 -----> calcul non local etape 2 on continue le calcul 66 : C a partir des seuils moyennes 67 : C 68 : C Modif L.Bode - 14/10/92 69 : C Nouveaux parametres en entree 70 : C JDIM Dimension de travail 71 : C ( Coques JDIM =2 , Massifs JDIM = IDIM ) 72 : C IFOUR2 Type de formulation 73 : C ( Coques IFOUR2 = -2 => contraintes planes , 74 : C Massifs IFOUR2 = IFOUR) 75 : C 76 : C variables en sortie 77 : C 78 : C VARINF variables internes finales 79 : C 80 : C SIGMAF contraintes finales 81 : C 82 : C C. LA BORDERIE MARS 1992 83 : C declaration des variables 84 : C 85 : C 86 : IMPLICIT REAL*8(A-H,O-Z) 87 : -INC CCOPTIO 88 : SEGMENT WRK0 89 : REAL*8 XMAT(NMATT) 90 : ENDSEGMENT 91 : * 92 : SEGMENT WRK1 93 : REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS) 94 : REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI) 95 : REAL*8 DEFP(NSTRS),XCAR(ICARA) 96 : ENDSEGMENT 97 : SEGMENT WRKK2 98 : REAL*8 EPSILI(NSTRS),DSIGT(NSTRS) 99 : ENDSEGMENT 100 : C Modif L.Bode - 14/10/92 101 : SEGMENT WRK3 102 : REAL*8 EPSILO(NSTRS) 103 : ENDSEGMENT 104 : C Definition dynamique de EPSILO 105 : C Fin modif L.Bode 106 : SEGMENT WRK5 107 : REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS) 108 : ENDSEGMENT 109 : C 110 : CHARACTER*8 CMATE 111 : INTEGER NSTRS,NVARI,NMATT,ISTEP 112 : REAL*8 EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3) 113 : REAL*8 SIGP(3),SIGPT(3),SIGPC(3),TRSIGT,TRSIGC 114 : C REAL*8 YOUN,XNU,EPSD0,ACOM,BCOM,ATRA,BTRA,BETA 115 : REAL*8 YOUN,XNU,EPSD0 116 : C ajout bodin 14 02 2006 117 : REAL*8 BETA,LOI 118 : C pour la loi L2R 119 : REAL*8 C,ALFA 120 : C pour la loi L3R 121 : REAL*8 ALFA1,ALFA2,ALFA3 122 : C fin ajout bodin 14 02 2006 123 : INTEGER ISTRS,JSTRS,KCAS,IRTD 124 : REAL*8 XZERO,UN,DEUX,TROIS 125 : REAL*8 EPSTL0,EPSTM0 126 : REAL*8 DINI,D,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA 127 : REAL*8 dt1 128 : C Fin modif Bodin - 26-09-2005 129 : PARAMETER (XZERO=0.D0,UN=1.D0,DEUX=2.D0,TROIS=3.D0) 130 : 131 : C 132 : C 133 : C recuperation des variables initiales dans les tableaux 134 : C 135 : C 136 : N=NSTRS 137 : CMATE = 'ISOTROPE' 138 : YOUN = XMAT(1) 139 : XNU = XMAT(2) 140 : EPSD0= XMAT(5) 141 : BETA = XMAT(6) 142 : LOI = XMAT(7) 143 : C = XMAT(8) 144 : ALFA = XMAT(9) 145 : ALFA1 = XMAT(10) 146 : ALFA2 = XMAT(11) 147 : ALFA3 = XMAT(12) 148 : C ajout bodin 17-11-2000 149 : C introduction de la variable EPSTIMI déformation équivalent initiale 150 : EPSTL0 = VAR0(1) 151 : C fin ajout bodin 152 : DINI = VAR0(2) 153 : C 154 : C 155 : C print*,'Je suis dans mazars' 156 : C 157 : C ajout bodin 26 09 2005 158 : C print*,'mazars Variable DT ',DT 159 : dt1 =DT 160 : C print*,'mazars variable dt1',dt1 161 : C fin ajout bodin 26 09 2005 162 : 163 : C 164 : C calcul des seuils 165 : C 166 : C 167 : C calcul de la deformation totale 168 : C 169 : SEGINI WRK3 170 : DO 100 ISTRS=1,NSTRS 171 : EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS) 172 : 100 CONTINUE 173 : C 174 : C calcul des deformations principales 175 : C 176 : C 177 : C on reecrit les deformations sous forme matricielle 178 : C 179 : C Modif L.Bode - 14/10/92 180 : C Rajout de IFOUR2 en argument de ENDOCA 181 : * print*,'on appelle ENDOCB' 181 : CALL ENDOCB (EPSILO,EPS33,2,IFOUR2) 182 : * print*,'apres endocb' 183 : C Fin modif L.Bode 184 : C 185 : C et on diagonalise 186 : C 187 : C Modif L.Bode - 14/10/92 188 : C Pour les elts Coques, on travaille en contraintes planes => JDIM =2 189 : C Pour les elts Massifs JDIM =IDIM 190 : * print*,'avant JACOB3' 190 : CALL JACOB3 (EPS33,JDIM,EPSIPP,VALP33) 191 : * print*,'apres JACOB3' 192 : C Fin modif L.Bode 193 : IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN 194 : C 195 : C on calcule la matrice de hooke et les contraintes ppales 196 : C 197 : C CMATE = 'ISOTROPE' 198 : C KCAS=1 199 : C Modif L.Bode - 14/10/92 200 : C IFOUR --> IFOUR2 dans appel DOHMAS 201 : C* print*,'avant dohmas' 201 : C CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS,KCAS,DDHOOK,IRTD) 202 : C* print*,'apres dohmas' 203 : C Fin modif L.Bode 204 : C DO 200 ISTRS=1,3 205 : C SIGP(ISTRS)= XZERO 206 : C DO 210 JSTRS=1,3 207 : C SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS) 208 : C210 CONTINUE 209 : C200 CONTINUE 210 : END IF 211 : C 212 : C on complete la deformation dans le cas des contraintes planes 213 : C 214 : C Modif L.Bode - 14/10/92 215 : C IFOUR remplace par IFOUR2 216 : IF (IFOUR2.EQ. -2) THEN 217 : EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU) 218 : END IF 219 : C Fin modif L.Bode 220 : C 221 : C on calcule le epsilontild 222 : Cmodif bodin 6 juillet 2001 pour le calcul de epsilon equivalente 223 : C on calcule la matrice de hooke et les contraintes ppales 224 : C 225 : CMATE = 'ISOTROPE' 226 : KCAS=1 227 : C Modif L.Bode - 14/10/92 228 : C IFOUR --> IFOUR2 dans appel DOHMAS 229 : * print*,'avant dohmas' 229 : CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS,KCAS,DDHOOK,IRTD) 230 : * print*,'apres dohmas' 231 : C Fin modif L.Bode 232 : 233 : DO 200 ISTRS=1,3 234 : SIGP(ISTRS)= XZERO 235 : DO 210 JSTRS=1,3 236 : SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS) 237 : 210 CONTINUE 238 : 200 CONTINUE 239 : C 240 : C on calcule le signe des contraintes elastiques 241 : C 242 : DO 300 ISTRS=1,3 243 : IF (SIGP(ISTRS).LT. XZERO) THEN 244 : SIGPC(ISTRS)=SIGP(ISTRS) 245 : SIGPT(ISTRS)=XZERO 246 : ELSE 247 : SIGPT(ISTRS)=SIGP(ISTRS) 248 : SIGPC(ISTRS)=XZERO 249 : END IF 250 : 300 CONTINUE 251 : TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3) 252 : TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3) 253 : C 254 : C on calcule les deformations dues aux contraintes positives 255 : C 256 : DO 400 ISTRS=1,3 257 : EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN 258 : 400 CONTINUE 259 : C 260 : EPSTIL=MAX( XZERO , EPSILT(1) )**2 + 261 : 1 MAX( XZERO , EPSILT(2) )**2 + 262 : 2 MAX( XZERO , EPSILT(3) )**2 263 : EPSTIL=SQRT (EPSTIL) 264 : C 265 : C EPSTIL=MAX( XZERO , EPSIPP(1) )**2 + 266 : C 1 MAX( XZERO , EPSIPP(2) )**2 + 267 : C 2 MAX( XZERO , EPSIPP(3) )**2 268 : C EPSTIL=SQRT (EPSTIL) 269 : C fin modif bodin 6 juillet 2001 pour le calcul de epsilon equivalente 270 : IF (ISTEP .EQ. 0) THEN 271 : EPSTIM=EPSTIL 272 : EPSTM0=EPSTL0 273 : VARF(1)=EPSTIL 274 : ELSE IF (ISTEP .EQ. 1) THEN 275 : VARF(2)=DINI 276 : VARF(1)=EPSTIL 277 : C ajout bodin 13-03-2001 pour une nouvelle variable interne 278 : C VARx(3) introduite pour garder la mémoire de la déf eq initiale 279 : C marche en non local si le champ de EPSTL0 est non local par exemple identiquement nul 280 : VARF(3)=EPSTL0 281 : C fin ajout bodin 13-03-2001 282 : ELSE IF (ISTEP .EQ. 2) THEN 283 : C ajout bodin 13-03-2001 pour une nouvelle variable interne 284 : C VARx(3) introduite pour garder la mémoire de la déf eq initiale 285 : EPSTM0=VAR0(3) 286 : C fin ajout bodin 13-03-2001 287 : EPSTIM=VAR0(1) 288 : VARF(1)=EPSTIM 289 : ELSE 290 : PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP 291 : END IF 292 : IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN 293 : C 294 : C on calcule l'endommagement et les contraintes 295 : C 296 : C 297 : C on calcule ALFAT ALFAC DT et DC puis D 298 : C dans le cas ou le seuil initial est depasse 299 : C 300 : C 301 : IF ( EPSTIM .GT. EPSD0) THEN 302 : C 303 : C calcul de l'endommagement 304 : C 305 : C 306 : C on calcule le signe des contraintes elastiques 307 : C 308 : C DO 300 ISTRS=1,3 309 : C IF (SIGP(ISTRS).LT. XZERO) THEN 310 : C SIGPC(ISTRS)=SIGP(ISTRS) 311 : C SIGPT(ISTRS)=XZERO 312 : C ELSE 313 : C SIGPT(ISTRS)=SIGP(ISTRS) 314 : C SIGPC(ISTRS)=XZERO 315 : C END IF 316 : C300 CONTINUE 317 : C TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3) 318 : C TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3) 319 : C 320 : C on calcule les deformations dues aux contraintes positives 321 : C 322 : C DO 400 ISTRS=1,3 323 : C EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN 324 : C400 CONTINUE 325 : C 326 : C on en deduit ALFAT et ALFAC 327 : C 328 : C ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) + 329 : C 1 MAX(XZERO,EPSIPP(2))*EPSILT(2) + 330 : C 2 MAX(XZERO,EPSIPP(3))*EPSILT(3) 331 : C ALFAT = ALFAT/(EPSTIL*EPSTIL) 332 : C ALFAC = UN - ALFAT 333 : C PRINT*,'alfat = ',ALFAT,' et alfac = ',ALFAC 334 : C PRINT*,'alfa1=',A99,'beta1 =',B99,'C1 = ',C99 335 : C 336 : C modification pour la bi ou tricompression 337 : C 338 : C IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN 339 : C GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+ 340 : C 1 SIGPC(3)*SIGPC(3) 341 : C GAMMA=-SQRT(GAMMA)/TRSIGC 342 : C EPSTIM=EPSTIM*GAMMA 343 : C END IF 344 : C 345 : C amelioration de la reponse en cisaillement pour beta > 1. 346 : C 347 : C IF (BETA .GT. UN) THEN 348 : C IF ( ALFAT .GT. XPETIT ) THEN 349 : C ALFAT=ALFAT**BETA 350 : C END IF 351 : C IF ( ALFAC .GT. XPETIT ) THEN 352 : C ALFAC=ALFAC**BETA 353 : C END IF 354 : C END IF 355 : C 356 : C 357 : C on calcule DT et DC puis D 358 : C dans le cas ou le seuil initial est depasse 359 : C 360 : C on est oblige de verifier car on a pu multiplier par gamma 361 : C 362 : C modif bodin 26-02-2001 363 : C modif bodin 26-02-2001 364 : C calcul incrémental pour le pas EPSTIM-EPSTIMI de déformation 365 : C équivalente 366 : C 367 : C 368 : C 369 : C 370 : C 371 : C 372 : C PRINT*,'Je suis dans mazars_bod et je suis content !!' 373 : C PRINT*,'La def initiale est',EPSTM0 374 : C PRINT*,'La def finale est',EPSTIM 375 : 376 : IF (LOI.EQ.DEUX) THEN 377 : D=((DINI**(1-ALFA))+ 378 : 1 ((C*(1-ALFA)/(BETA+1))* 379 : 2 ((EPSTIM**(BETA+1))-(EPSD0**(BETA+1)))*(dt1*2) 380 : 3 ))**(1/(1-ALFA)) 381 : C 382 : ELSE IF (LOI.EQ.TROIS) THEN 383 : C on traite la loi L3R 384 : C calcul d'un term intermediaire 385 : TERM_A= (((ALFA1*(UN -EXP(-1*((DINI/ALFA2)**ALFA3))))+ 386 : 1 ((((EPSTIM**(BETA+1))-(EPSD0**(BETA+1))) 387 : 2 /(BETA+UN))*(dt1*2)) 388 : 3 )/ALFA1) 389 : C 390 : IF (TERM_A .LT. UN) THEN 391 : D= ALFA2 *((-1*(LOG(UN - TERM_A)))**(UN/ALFA3)) 392 : ELSE 393 : D=UN 394 : END IF 395 : C sinon roblème de définition de LOI 396 : ELSE 397 : WRITE(*,*) 'LOI doit etre egal à 2. ou 3.' 398 : STOP 399 : END IF 400 : C 401 : C on borne la valeur de D a 0.99999 402 : C 403 : D=MIN ( D , UN-1.D-05 ) 404 : ELSE 405 : C on n'a pas dépassé le seuil 406 : D=XZERO 407 : END IF 408 : C 409 : C on teste la croissance de D 410 : C 411 : D=MAX ( D , DINI ) 412 : C 413 : C on le stocke dans les variables internes finales 414 : C 415 : VARF(2)= D 416 : C 417 : C on calcule les contraintes finales 418 : C 418 : CALL MATVE1 (DDHOOK,EPSILO,NSTRS,NSTRS,SIGF,2) 419 : DO 500 ISTRS=1,NSTRS 420 : SIGF(ISTRS)=SIGF(ISTRS)*(UN-D) 421 : 500 CONTINUE 422 : C 423 : C et les deformations inelastiques finales 424 : C 425 : DO 600 ISTRS=1,NSTRS 426 : EPINF(ISTRS)=EPSILO(ISTRS)*D 427 : 600 CONTINUE 428 : 429 : END IF 430 : SEGSUP WRK3 431 : RETURN 432 : END 433 :
© Cast3M 2003 - All rights reserved.
Disclaimer