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