1 : C GURSO2 SOURCE CHAT 05/01/13 00:22:19 5004 2 : SUBROUTINE GURSO2(DEPST,NSTRS,MFR1,IB,IGAU, 3 : 1 DSIGT,NCOMAT,SIG0,VAR0,XMAT,XCAR,NVARI,ICARA, 4 : 2 SIGF,VARF,DEFP,TRAC,KERRE,necou) 5 : * 6 : *_________________________________________________________________ 7 : * 8 : * 9 : * ENTREES : 10 : * --------- 11 : * 12 : * DEPST = INCREMENT DE DEFORMATIONS TOTALES 13 : * NSTRS = NBRE DE COMPOSANTES DES DEFORMATIONS 14 : * NCOMAT= NBRE DE CARACTERISTIQUES MECANIQUES DU MATERIAU 15 : * MFR1 = NUMERO DE LA FORMULATION 16 : * IB = NUMERO DE L ELEMENT COURANT 17 : * IGAU = NUMERO DU POINT COURANT 18 : * NVARI = NBRE DE VARIABLES INTERNES 19 : * SIG0(NSTRS) = CONTRAINTES AU DEBUT DU PAS D'INTEGRATION 20 : * VAR0(NVARI) = VARIABLES INTERNES AU DEBUT DU PAS DE TEMPS 21 : * XMAT(NCOMAT) = CARACTERISTIQUES MECANIQUES DU MATERIAU 22 : * ICARA = NBRE DE CARACTERISTIQUES GEOMETRIQUES DES ELEMENTS FINIS 23 : * XCARA(ICARA) = CARACTERISTIQUES GEOMETRIQUES DES ELEMENTS FINIS 24 : * TRAC = COURBE DE TRACTION NON ENDOMMAGEE 25 : * 26 : * SORTIE : 27 : * -------- 28 : * 29 : * SIGF(NSTRS) = CONTRAINTES FINALES 30 : * VARF(NVARI) = VARIALES INTERNES A LA FIN DU PAS D'INTEGRATION 31 : * DEFP(NSTRS) = INCREMENT DE DEFORMATION PLASTIQUE A LA FIN 32 : * DU PAS D'INTEGRATION 33 : * ============================================================ 34 : * ICI IL FAUT PROGRAMMER EN FORTRAN PUR 35 : *============================================================= 36 : * 37 : IMPLICIT INTEGER(I-N) 38 : IMPLICIT REAL*8(A-H,O-Z) 39 : * 40 : -INC CCREEL 41 : -INC CCOPTIO 42 : * 43 : SEGMENT NECOU 44 : * COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO, 45 : INTEGER NCOURB,IPLAST,IT,IMAPLA,ISOTRO, 46 : . ITYP,IFOURB,IFLUAG, 47 : . ICINE,ITHER,IFLUPL,ICYCL,IBI, 48 : . JFLUAG,KFLUAG,LFLUAG, 49 : . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF 50 : ENDSEGMENT 51 : C COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO, 52 : C . ITYP,IFOURB,IFLUAG, 53 : C . ICINE,ITHER,IFLUPL,ICYCL,IBI, 54 : C . JFLUAG,KFLUAG,LFLUAG, 55 : C . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF 56 : * 57 : DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*),SIGF(*), 58 : & VARF(*),DEFP(*),DSIGT(*),TRAC(*) 59 : DIMENSION RDEPS0(6) 60 : DIMENSION RSIG0(6),RDEPS(6),RSIGF(6),RDEFP(6),RDIGT(6) 61 : DIMENSION DEVT(6),SIGT(6),DEFT(6),SIGT1(6) 62 : * 63 : * Adaptation de l'option de calcul vers le 3D massif de SIG0 a RSIG0 64 : *==================================================================== 65 : * 66 : IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN 67 : * 68 : *---> 1 formulation massive 69 : *---> 2 formulation quasi incompressible 70 : *---> MASSIF 3D 71 : * 72 : IF (NSTRS .EQ. 6) THEN 73 : DO 110 I=1,NSTRS 74 : RSIG0(I)=SIG0(I) 75 : RDEPS0(I)=DEPST(I) 76 : 110 CONTINUE 77 : ELSE IF ( NSTRS .EQ. 4 .AND. ((IFOUR .EQ. 0) 78 : & .OR.(IFOUR .EQ. -1).OR.(IFOUR.EQ.-2))) THEN 79 : * 80 : *---> Calcul en mode deformations planes ou axisymetrique 81 : * ou contraintes planes 82 : * 83 : DO 115 I=1,NSTRS 84 : RSIG0(I)=SIG0(I) 85 : RDEPS0(I)=DEPST(I) 86 : 115 CONTINUE 87 : RSIG0(5)=0.D0 88 : RSIG0(6)=0.D0 89 : RDEPS0(5)=0.D0 90 : RDEPS0(6)=0.D0 91 : ENDIF 92 : ELSE 93 : KERRE = 99 94 : RETURN 95 : ENDIF 96 : * 97 : * 98 : * Passage des deformations de cisaillement exprimées 99 : * en GAMA aux déformations de cisaillement exprimées 100 : * en déformations 101 : * 102 : DO 116 I=1,6 103 : A=1.D0 104 : IF (I.GT.3) A=2.D0 105 : RDEPS0(I)=RDEPS0(I)/A 106 : RDEPS(I)=0.D0 107 : 116 CONTINUE 108 : * 109 : * Données du materiau 110 : *=========================================================== 111 : * 112 : YOUNG=XMAT(1) 113 : XNU=XMAT(2) 114 : Q0=XMAT(6) 115 : FU0=XMAT(7) 116 : FF0=XMAT(8) 117 : FC0=XMAT(9) 118 : FNS0=XMAT(10) 119 : FNE0=XMAT(11) 120 : SNS0=XMAT(12) 121 : SNE0=XMAT(13) 122 : SIGN0=XMAT(14) 123 : EPSN0=XMAT(15) 124 : F0=XMAT(16) 125 : XSRMA=XMAT(17) 126 : Fmin0=1.D-8 127 : IF ((FF0.LE.Fmin0).OR.(FF0.GT.1.D0)) FF0=1.D0 128 : IF ((F0.LE.Fmin0).OR.(F0.GE.1.D0)) F0=0.D0 129 : C RMIN0=(1.D0-FF0)/(1.D0-F0) 130 : * 131 : * Calculs préliminaires 132 : *=================================== 133 : * 134 : *---> Module de cisaillement G0 et de pression hydro XK0 135 : * 136 : G0=YOUNG/(2.D0*(1.D0+XNU)) 137 : XK0=YOUNG/(3.D0*(1.D0-(2.D0*XNU))) 138 : * 139 : * Début des sous incrémentations 140 : *=================================== 141 : * 142 : recal0=-100.D0 143 : iter2=0 144 : nmax0=1 145 : * 146 : *---> Variables internes test 147 : * 148 : * - déformation plastique équivalente 149 : * - fraction volumique de cavités 150 : * 151 : 93 EPSPT=VAR0(1) 152 : FT=VAR0(2) 153 : SIGMT=VAR0(3) 154 : EPSMT=VAR0(4) 155 : RHOT=VAR0(5) 156 : FNST=VAR0(6) 157 : FNET=VAR0(7) 158 : VARF(3)=VAR0(3) 159 : VARF(4)=VAR0(4) 160 : VARF(6)=VAR0(6) 161 : VARF(7)=VAR0(7) 162 : IF (RHOT.LE.1.D-30) THEN 163 : RHOT=1.D0 164 : VAR0(5)=1.D0 165 : ENDIF 166 : IF (FT.LE.Fmin0) FT=0.D0 167 : * 168 : * Contraintes planes 169 : * - on stocke les variables internes au début de 170 : * chaque sous incrémentation 171 : * - on stocke les déformations plastiques au début de 172 : * chaque sous incrémentation 173 : * 174 : IF (IFOUR.EQ.-2) THEN 175 : EPSPT1=EPSPT 176 : FT1=FT 177 : EPSMT1=EPSMT 178 : SIGMT1=SIGMT 179 : RHOT1=RHOT 180 : FNST1=FNST 181 : FNET1=FNET 182 : def1=0.D0 183 : def2=0.D0 184 : def3=0.D0 185 : rdeps3=0.D0 186 : ENDIF 187 : * 188 : * Contraintes planes 189 : * - on stocke les contraintes au début de 190 : * chaque sous incrémentation 191 : * 192 : IF (IFOUR.EQ.-2) THEN 193 : DO 171 I=1,6 194 : SIGT1(I)=RSIG0(I) 195 : 171 CONTINUE 196 : ENDIF 197 : DO 172 I=1,6 198 : SIGT(I)=RSIG0(I) 199 : 172 CONTINUE 200 : SMT0=(RSIG0(1)+RSIG0(2)+RSIG0(3))/3.D0 201 : SMT=SMT0 202 : IF (IFOUR.EQ.-2) SMT1=SMT0 203 : * 204 : * On divise le chargement pour les sous incrementations 205 : * 206 : 95 DO 118 I=1,6 207 : RDEPS(I)=RDEPS0(I)/nmax0 208 : 118 CONTINUE 209 : iter2=iter2+1 210 : * 211 : * Initialisations 212 : *=================================== 213 : * 214 : iter1=0 215 : * 216 : *---> paramètre de recalcul éventuel de la nucléation 217 : * il est ré-initialisé à chaque sous itération 218 : * 219 : recal1=-100.D0 220 : * 221 : *---> Calcul de la trace de l'incrément de déformation 222 : * 223 : IF (IFOUR.EQ.-2) THEN 224 : XLAM0=(1.D0-FT)*XK0-(2.D0*G0/3.D0) 225 : RDEPS(3)=-1.D0*XLAM0/(2.D0*G0+XLAM0)*(RDEPS(1)+RDEPS(2)) 226 : ENDIF 227 : 98 treps0=RDEPS(1)+RDEPS(2)+RDEPS(3) 228 : * 229 : * 230 : * Calculs des grandeurs élastiques test 231 : *==================================================== 232 : * 233 : *---> Variables internes test: cas contraintes planes 234 : * 235 : * - déformation plastique équivalente 236 : * - fraction volumique de cavités 237 : * 238 : IF (IFOUR.EQ.-2) THEN 239 : EPSPT=EPSPT1 240 : FT=FT1 241 : EPSMT=EPSMT1 242 : SIGMT=SIGMT1 243 : RHOT=RHOT1 244 : FNST=FNST1 245 : FNET=FNET1 246 : SMT=SMT1 247 : ENDIF 248 : * 249 : *---> Calcul des contraintes test élastiques 250 : * et de l'incrément des déformations plastiques test au 251 : * cours du pas 252 : * 253 : DO 101 I=1,6 254 : A=0.D0 255 : IF (I.LE.3) A=1.D0 256 : IF (IFOUR.EQ.-2) SIGT(I)=SIGT1(I) 257 : DEVT(I)=SIGT(I)-(A*SMT) 258 : DEVT(I)=DEVT(I)+(2.D0*G0*(RDEPS(I)-(A*treps0/3.D0))) 259 : DEFT(I)=0.D0 260 : 101 CONTINUE 261 : * 262 : *---> Contrainte équivalente STEST 263 : * 264 : STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0 265 : IF (STEST2.LE.1.D-10) STEST2=0.D0 266 : Stest=STEST2**0.5D0 267 : * 268 : *---> Limite d'élasticité lue sur la courbe d'écrouissage 269 : * 270 : * - Recherche des valeurs de déformations plastiques équivalentes 271 : * encadrant la valeur courante 272 : * 273 : EPSM0=TRAC(2*NCOURB) 274 : IF (EPSMT.GE.EPSM0) THEN 275 : Y1=TRAC(2*NCOURB-1) 276 : Y2=TRAC(2*NCOURB-1) 277 : EPS1=EPSM0 278 : EPS2=2.D0*EPSM0 279 : ELSE 280 : DO 300 I=1,(NCOURB-1) 281 : EPS1=TRAC(2*I) 282 : EPS2=TRAC(2*(I+1)) 283 : IF ((EPSMT.GE.EPS1).AND.(EPSMT.LT.EPS2)) THEN 284 : Y1=TRAC(2*I-1) 285 : Y2=TRAC(2*(I+1)-1) 286 : GOTO 96 287 : ENDIF 288 : 300 CONTINUE 289 : ENDIF 290 : * 291 : * - Limite d'élasticité 292 : * 293 : 96 H1=(Y2-Y1)/(EPS2-EPS1) 294 : Ytest=(H1*(EPSMT-EPS2))+Y2 295 : * 296 : *---> Fraction de densité test 297 : * 298 : RHOT=1.D0/((1.D0/RHOT)+treps0) 299 : IF (RHOT.LT.1.D-10) RHOT=1.D-10 300 : ** IF ((FT.GT.Fmin0).AND.(RHOT.LT.RMIN0)) RHOT=RMIN0 301 : * 302 : *---> Coefficients de nucléation des cavités 303 : * 304 : dif0=VARF(3)-Ytest-SMT 305 : IF ((dif0.LE.1.D-5).AND.(treps0.GE.0.D0)) THEN 306 : BB0=EXP(-.5D0*((VARF(3)-SIGN0)/SNS0)**2.D0) 307 : BB0=BB0*FNS0/(SNS0*((2.D0*XPI)**.5D0)) 308 : ELSE 309 : BB0=0.D0 310 : ENDIF 311 : * 312 : *---> Cas de la nucléation pilotée en contrainte 313 : * 314 : IF ((FNS0.GT.Fmin0).AND.(recal1.LT.0.D0)) THEN 315 : * 316 : * Si on se trouve dans la "fenetre" où la 317 : * nucléation apparait, on sous incrémente le pas 318 : * de calcul afin de suivre les variations de BB0 319 : * 320 : dif0=Ytest+SMT 321 : IF ((dif0.GE.(SIGN0-SNS0*2.D0)).AND. 322 : . (VARF(3).LE.(SIGN0+SNS0*2.D0))) THEN 323 : xmax1=200.D0*(dif0-VARF(3))/(4.D0*SNS0) 324 : nmax1=NINT(xmax1)*nmax0 325 : IF (nmax1.GT.2000) nmax1=2000 326 : IF (nmax1.GT.nmax0) THEN 327 : nmax0=nmax1 328 : iter2=0 329 : recal1=100.D0 330 : GOTO 93 331 : ENDIF 332 : ENDIF 333 : * 334 : ENDIF 335 : * 336 : *---> Fraction de cavité test 337 : * On résoud une équation du second degré 338 : * 339 : AA1=BB0*XK0/((1.D0-F0)*RHOT) 340 : BB1=BB0*XK0-1.D0 341 : CC1=BB0*SMT+1.D0-FT 342 : delta0=BB1*BB1+4.D0*AA1*CC1 343 : * 344 : IF (AA1.GT.1.D-5) THEN 345 : * 346 : * Le terme en contrainte de la nucléation n'est pas négligeable 347 : * 348 : XX1=(BB1+(delta0**.5D0))/(2.D0*AA1) 349 : FT=1.D0-XX1 350 : ELSE 351 : * 352 : * Le terme en contrainte de la nucléation est négligé 353 : * 354 : BB0=0.D0 355 : ENDIF 356 : FT=MAX(FT,0.D0) 357 : * 358 : *---> Contrainte totale et contrainte moyenne test 359 : * 360 : FNST=FNST-BB0*SMT 361 : SMT=(1.D0-FT)*XK0*((1.D0-FT+XSRMA*(FNET+FNST))/(1.D0-F0)/RHOT- 362 : . 1.D0) 363 : FNST=FNST+BB0*SMT 364 : DO 200 I=1,6 365 : A=0.D0 366 : IF (I.LE.3) A=1.D0 367 : SIGT(I)=DEVT(I)+(A*SMT) 368 : 200 CONTINUE 369 : * 370 : *---> Fonction d'endommagement 371 : * 372 : IF (FT.LE.FC0) THEN 373 : FST0=FT 374 : ddel0=1.D0 375 : ELSE 376 : FST0=FC0+(FU0-FC0)/(FF0-FC0)*(FT-FC0) 377 : ddel0=(FU0-FC0)/(FF0-FC0) 378 : ENDIF 379 : IF (FST0.GT.FU0) FST0=FU0 380 : * 381 : *---> Terme d'endommagement 382 : * 383 : SMTM=MAX(SMT,0.D0) 384 : SMTM=SMT 385 : Gtest=1.D0+(Q0*FST0)*(Q0*FST0)- 386 : . (2.D0*Q0*FST0*COSH(3.D0*SMTM/(2.D0*Ytest))) 387 : * 388 : *---> Fonction de charge test PHIT 389 : * 390 : PHIT=(Stest*Stest)-(Ytest*Ytest)*Gtest 391 : * 392 : *---> Cas du matériau complètement endommagé 393 : * 394 : 91 IF (FT.GE.FF0) THEN 395 : * 396 : * Les cavités 397 : * 398 : FT=FF0+Fmin0 399 : * 400 : * Les contraintes sont mises a zero 401 : * Les deformations sont plastiques 402 : * 403 : SMT=0.D0 404 : Stest=0.D0 405 : treps1=(SMT/(1.D0-FF0)-SMT0/(1.D0-VAR0(2)))/XK0 406 : DO 12 I=1,6 407 : SIGT(I)=0.D0 408 : A=0.D0 409 : IF (I.LE.3) A=1.D0 410 : depel0=SIGT(I)-RSIG0(I) 411 : DEVT(I)=(depel0-(SMT-SMT0)*A)/(2.D0*G0) 412 : depel0=DEVT(I)+(treps1*A/3.D0) 413 : DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0 414 : IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN 415 : DEFT(3)=rdeps3+RDEPS(3)-depel0 416 : ENDIF 417 : 12 CONTINUE 418 : * 419 : * Deformation plastique equivalente 420 : * 421 : treps1=DEFT(1)+DEFT(2)+DEFT(3) 422 : DO 14 I=1,6 423 : A=0.D0 424 : IF (I.LE.3) A=1.D0 425 : DEVT(I)=DEFT(I)-treps1*A/3.D0 426 : 14 CONTINUE 427 : IF (treps1.LE.0.D0) THEN 428 : DO 16 I=1,3 429 : DEFT(I)=DEVT(I) 430 : 16 CONTINUE 431 : ENDIF 432 : dlam0=(2.D0*PROCON(DEVT,DEVT,6)/3.D0)**.5D0 433 : EPSPT=EPSPT+dlam0 434 : ** RHOT=RMIN0 435 : * 436 : Gtest=0.D0 437 : PHIT=0.D0 438 : * 439 : ENDIF 440 : * 441 : * 442 : * Vérification du critère de plasticité 443 : *========================================================= 444 : * 445 : *---> Erreur admise 446 : * 447 : ERR0=1.D-7*(ABS(PHIT)) 448 : ERR0=MAX(ERR0,1.D-8*Ytest*Ytest) 449 : * 450 : *---> Critère de plasticité 451 : * 452 : PHI0=PHIT 453 : iter0=0 454 : * 455 : 99 IF (PHI0.LE.ERR0) THEN 456 : * 457 : * On est élastique 458 : *========================================================= 459 : * 460 : *-------------------------------------------------------------- 461 : * Cas particulier des contraintes planes 462 : *-------------------------------------------------------------- 463 : * 464 : IF (IFOUR.EQ.-2) THEN 465 : * 466 : *---> On vérifie qu'on est en contraintes planes 467 : * 468 : iter1=iter1+1 469 : SIG3=ABS(SIGT(3)) 470 : SIG30=MAX(Stest*1.D-5,1.D-3) 471 : IF (SIG3.GT.SIG30) THEN 472 : XLAM0=(1.D0-FT)*XK0-(2.D0*G0/3.D0) 473 : RDEPS(3)=-1.D0*(XLAM0/(2.D0*G0+XLAM0))* 474 : . (RDEPS(1)-DEFT(1)+RDEPS(2)-DEFT(2)+def1+def2) 475 : RDEPS(3)=RDEPS(3)+DEFT(3)-def3 476 : IF (iter1.LE.200)THEN 477 : GOTO 98 478 : ELSE 479 : KERRE=460 480 : ENDIF 481 : ELSE 482 : EPSPT1=EPSPT 483 : FT1=FT 484 : SIGMT1=SIGMT 485 : EPSMT1=EPSMT 486 : RHOT1=RHOT 487 : FNST1=FNST 488 : FNET1=FNET 489 : SMT1=SMT 490 : SIGT1(1)=SIGT(1) 491 : SIGT1(2)=SIGT(2) 492 : SIGT1(3)=SIGT(3) 493 : SIGT1(4)=SIGT(4) 494 : SIGT1(5)=SIGT(5) 495 : SIGT1(6)=SIGT(6) 496 : rdeps3=rdeps3+RDEPS(3) 497 : def1=DEFT(1) 498 : def2=DEFT(2) 499 : def3=DEFT(3) 500 : ENDIF 501 : ENDIF 502 : * 503 : *----------------------------------------------------- 504 : * Fin du cas particulier des contraintes planes 505 : *----------------------------------------------------- 506 : * 507 : DO 400 I=1,6 508 : RSIGF(I)=SIGT(I) 509 : RDEFP(I)=DEFT(I) 510 : RDIGT(I)=RSIGF(I)-RSIG0(I) 511 : 400 CONTINUE 512 : VARF(1)=EPSPT 513 : VARF(2)=FT 514 : VARF(3)=MAX(VARF(3),Ytest+SMT) 515 : VARF(4)=MAX(VARF(4),EPSMT) 516 : VARF(5)=RHOT 517 : VARF(6)=FNST 518 : VARF(7)=FNET 519 : * 520 : ELSE 521 : * 522 : * On plastifie 523 : *========================================================= 524 : * 525 : *---> Valeur minimale de Stest au dessous de laquelle 526 : * les contraintes sont considérées comme étant 527 : * purement hydrostatiques 528 : * 529 : Gmin0=2.D0*Q0*FST0*COSH(3.D0*SMTM/(2.D0*Ytest)) 530 : Gmin0=Gmin0-((Q0*FST0)*(Q0*FST0)) 531 : Gmin0=(ABS(Gmin0))**.5D0 532 : Gmin0=Gmin0*Ytest*1.D-6 533 : * 534 : *---> Coefficient de nucléation 535 : * 536 : dif0=VARF(4)-EPSMT 537 : IF (dif0.LE.1.D-5) THEN 538 : DD0=EXP(-.5D0*((VARF(4)-EPSN0)/SNE0)**2.D0) 539 : DD0=DD0*FNE0/(SNE0*((2.D0*XPI)**.5D0)) 540 : ELSE 541 : DD0=0.D0 542 : ENDIF 543 : * 544 : *---> Condition de consistance 545 : * 546 : IF (Stest.LT.Gmin0) THEN 547 : tr0=3.D0*Q0*FST0*Ytest*SINH(3.D0*SMTM/(2.D0*Ytest)) 548 : DE0=SMT*tr0/((1.D0-FT+XSRMA*(FNET+FNST))*Ytest) 549 : ELSE 550 : tr0=3.D0*Q0*FST0*Ytest*SINH(3.D0*SMTM/(2.D0*Ytest)) 551 : tr0=tr0/(2.D0*Stest) 552 : DE0=(SMT*tr0+Stest)/((1.D0-FT+XSRMA*(FNET+FNST))*Ytest) 553 : ENDIF 554 : * 555 : H0=H1*DE0 556 : DS1=XK0*(2.D0*(1.D0-FT+XSRMA*(FNET+FNST))/(1.D0-F0)/RHOT-1.D0) 557 : DF0=((1.D0-FT+XSRMA*(FNET+FNST))*tr0+DD0*DE0+BB0*H0)/ 558 : . (1.D0+BB0*DS1) 559 : * 560 : DCOS0=(DS1*DF0/(2.D0*Ytest))+(SMT*H0)/(2.D0*Ytest*Ytest) 561 : DCOS0=DCOS0*3.D0*SINH(3.D0*SMTM/(2.D0*Ytest)) 562 : * 563 : DG0=2.D0*Q0*Q0*ddel0*DF0*FST0 564 : DG0=DG0+2.D0*Q0*FST0*DCOS0 565 : DG0=DG0-2.D0*Q0*ddel0*DF0*COSH(3.D0*SMTM/(2.D0*Ytest)) 566 : * 567 : IF (Stest.LT.Gmin0) THEN 568 : denom0=Ytest*Ytest*DG0 569 : ELSE 570 : denom0=6.D0*G0*Stest+Ytest*Ytest*DG0 571 : denom0=denom0+2.D0*H0*Gtest*Ytest 572 : ENDIF 573 : * 574 : IF (denom0.GT.1.D-5) THEN 575 : dlam0=PHIT/denom0 576 : ELSE 577 : IF (recal0.LE.0.D0) THEN 578 : recal0=100.D0 579 : nmax0=2000 580 : iter2=0 581 : GOTO 93 582 : ENDIF 583 : ENDIF 584 : * 585 : * 586 : * 587 : * Vérification des hypothèses de calcul 588 : *=============================================================== 589 : * 590 : IF (recal0.LE.0.D0) THEN 591 : * 592 : xnum0=10.D0 593 : * 594 : xlim0=1.D-1 595 : xmax1=0.D0 596 : xmax2=0.D0 597 : xmax3=0.D0 598 : xmax4=0.D0 599 : xmax5=0.D0 600 : xmax6=0.D0 601 : xmax7=0.D0 602 : * 603 : *---> Condition 3.D0*G*dlam0/Stest << 1 604 : * 605 : IF (Stest.GE.Gmin0) THEN 606 : cri01=ABS(3.D0*G0*dlam0/Stest) 607 : IF (cri01.GT.xlim0) THEN 608 : xmax1=xnum0*cri01+5.D0 609 : ENDIF 610 : ENDIF 611 : * 612 : *---> Condition H0*dlam0/Ytest << 1 613 : * 614 : cri01=ABS(H0*dlam0/Ytest) 615 : IF (cri01.GT.xlim0) THEN 616 : xmax2=xnum0*cri01+5.D0 617 : ENDIF 618 : * 619 : *---> Condition DCOS0*dlam0/COSH << 1 620 : * 621 : cri01=ABS(DCOS0*dlam0/COSH(3.D0*SMTM/(2.D0*Ytest))) 622 : IF (cri01.GT.xlim0) THEN 623 : xmax3=xnum0*cri01+5.D0 624 : ENDIF 625 : * 626 : *---> Condition DF0*dlam0/F << 1 627 : * 628 : Fmax0=MAX(FNS0,FNE0) 629 : Fmax0=MAX(Fmax0,F0) 630 : Fmax0=MAX(Fmax0,FT) 631 : Fmax0=MAX(Fmax0,Fmin0) 632 : cri01=ABS(DF0*dlam0/Fmax0) 633 : IF (cri01.GT.xlim0) THEN 634 : xmax4=xnum0*cri01+5.D0 635 : ENDIF 636 : * 637 : *---> Condition DG0*dlam0/G << 1 638 : * 639 : Gmax0=MAX(ABS(Gtest),1.D0) 640 : cri01=ABS(DG0*dlam0/Gmax0) 641 : IF (cri01.GT.xlim0) THEN 642 : xmax5=xnum0*cri01+5.D0 643 : ENDIF 644 : * 645 : *---> Condition DS1*dlam0/SMT << 1 646 : * 647 : SMTmax=MAX(ABS(SMT),(1.D-3*TRAC(1))) 648 : cri01=ABS(DS1*DF0*dlam0/SMTmax) 649 : IF (cri01.GT.xlim0) THEN 650 : xmax6=xnum0*cri01+5.D0 651 : ENDIF 652 : * 653 : *---> Condition tr0*dlam0/tr(EPSP) << 1 654 : * 655 : tremp0=(FT-F0)/(1.D0-F0)/RHOT 656 : IF (ABS(tremp0).GT.1.D-5) THEN 657 : cri01=ABS(tr0*dlam0/tremp0) 658 : IF (cri01.GT.xlim0) THEN 659 : xmax7=xnum0*cri01+5.D0 660 : ENDIF 661 : ENDIF 662 : * 663 : * 664 : IF (xmax1.GT.xmax2) THEN 665 : xmax00=xmax1 666 : ELSE 667 : xmax00=xmax2 668 : ENDIF 669 : IF (xmax3.GT.xmax00) xmax00=xmax3 670 : IF (xmax4.GT.xmax00) xmax00=xmax4 671 : IF (xmax5.GT.xmax00) xmax00=xmax5 672 : IF (xmax6.GT.xmax00) xmax00=xmax6 673 : IF (xmax7.GT.xmax00) xmax00=xmax7 674 : * 675 : nmax00=NINT(xmax00) 676 : recal0=100.D0 677 : IF (nmax00.GT.nmax0) nmax0=nmax00 678 : IF (nmax0.GT.2000) nmax0=2000 679 : * 680 : *---> On recommence le calcul en sous incrémentant le 681 : * pas de temps 682 : * 683 : IF (nmax0.GE.2) THEN 684 : iter2=0 685 : GOTO 93 686 : ENDIF 687 : * 688 : ENDIF 689 : * 690 : * Cas particulier de la nucléation 691 : *============================================================ 692 : * 693 : IF (recal1.LT.0.D0) THEN 694 : * 695 : *---> Condition sur la nucléation pilotée en contrainte 696 : * 697 : IF (FNS0.GT.Fmin0) THEN 698 : * 699 : * Si on se trouve dans la fenetre, on sous incremente 700 : * 701 : dif0=Ytest+H0*dlam0+SMT-DS1*DF0*dlam0 702 : IF ((dif0.GE.(SIGN0-SNS0*2.D0)).AND. 703 : . (VARF(3).LE.(SIGN0+SNS0*2.D0))) THEN 704 : xmax1=200.D0*(dif0-VARF(3))/(4.D0*SNS0) 705 : nmax1=NINT(xmax1)*nmax0 706 : IF (nmax1.GT.2000) nmax1=2000 707 : IF (nmax1.GT.nmax0) THEN 708 : nmax0=nmax1 709 : iter2=0 710 : recal1=100.D0 711 : GOTO 93 712 : ENDIF 713 : ENDIF 714 : * 715 : ENDIF 716 : * 717 : *---> Condition sur la nucléation pilotée en deformation 718 : * 719 : IF (FNE0.GT.Fmin0) THEN 720 : * 721 : * Si on se trouve dans la fenetre, on sous incremente 722 : * 723 : dif0=EPSMT+DE0*dlam0 724 : IF ((dif0.GE.(EPSN0-SNE0*2.D0)).AND. 725 : . (VARF(4).LE.(EPSN0+SNE0*2.D0))) THEN 726 : xmax1=200.D0*(dif0-VARF(4))/(4.D0*SNE0) 727 : nmax1=NINT(xmax1)*nmax0 728 : IF (nmax1.GT.2000) nmax1=2000 729 : IF (nmax1.GT.nmax0) THEN 730 : nmax0=nmax1 731 : iter2=0 732 : recal1=100.D0 733 : GOTO 93 734 : ENDIF 735 : ENDIF 736 : * 737 : ENDIF 738 : * 739 : ENDIF 740 : * 741 : * 742 : * Mise à jour des grandeurs après écoulement plastique 743 : *============================================================ 744 : * 745 : *---> Variables internes mises à jour 746 : * 747 : IF (Stest.GT.Gmin0) EPSPT=EPSPT+dlam0 748 : EPSMT=EPSMT+DE0*dlam0 749 : FNST=FNST+BB0*(H0-DS1*DF0)*dlam0 750 : FNET=FNET+DD0*DE0*dlam0 751 : FT=FT+DF0*dlam0 752 : IF (FT.LT.Fmin0) FT=0.D0 753 : * 754 : *---> Contraintes mises à jour 755 : * 756 : SMT=(1.D0-FT)*XK0*((1.D0-FT+XSRMA*(FNET+FNST))/RHOT/(1.D0-F0)- 757 : . 1.D0) 758 : DO 601 I=1,6 759 : A=0.D0 760 : IF (I.LE.3) A=1.D0 761 : IF (Stest.LE.Gmin0) THEN 762 : DEVT(I)=0.D0 763 : ELSE 764 : DEVT(I)=DEVT(I)-3.D0*G0*DEVT(I)*dlam0/Stest 765 : ENDIF 766 : SIGT(I)=DEVT(I)+(SMT*A) 767 : 601 CONTINUE 768 : * 769 : *---> Incrément de déformations plastiques mise à jour 770 : * 771 : treps1=(SMT/(1.D0-FT)-SMT0/(1.D0-VAR0(2)))/XK0 772 : DO 610 I=1,6 773 : A=0.D0 774 : IF (I.LE.3) A=1.D0 775 : depel0=(DEVT(I)-RSIG0(I)+(SMT0*A))/(2.D0*G0) 776 : depel0=depel0+(treps1*A/3.D0) 777 : DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0 778 : IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN 779 : DEFT(3)=rdeps3+RDEPS(3)-depel0 780 : ENDIF 781 : 610 CONTINUE 782 : * 783 : *---> Contrainte équivalente Stest mise à jour 784 : * 785 : STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0 786 : IF (STEST2.LE.1.D-10) STEST2=0.D0 787 : Stest= STEST2 ** 0.5D0 788 : * 789 : *---> Limite d'élasticité lue sur la courbe d'écrouissage 790 : * mise à jour 791 : * 792 : * - Recherche des valeurs de déformations plastiques équivalentes 793 : * encadrant la valeur courante 794 : * 795 : EPSM0=TRAC(2*NCOURB) 796 : IF (EPSMT.GE.EPSM0) THEN 797 : Y1=TRAC(2*NCOURB-1) 798 : Y2=TRAC(2*NCOURB-1) 799 : EPS1=EPSM0 800 : EPS2=2.D0*EPSM0 801 : ELSE 802 : DO 700 I=1,(NCOURB-1) 803 : EPS1=TRAC(2*I) 804 : EPS2=TRAC(2*(I+1)) 805 : IF ((EPSMT.GE.EPS1).AND.(EPSMT.LT.EPS2)) THEN 806 : Y1=TRAC(2*I-1) 807 : Y2=TRAC(2*(I+1)-1) 808 : GOTO 97 809 : ENDIF 810 : 700 CONTINUE 811 : ENDIF 812 : * 813 : * - Limite d'élasticité mise à jour 814 : * 815 : 97 H1=(Y2-Y1)/(EPS2-EPS1) 816 : Ytest=(H1*(EPSMT-EPS2))+Y2 817 : * 818 : *---> Mise à jour des coefficients de nucléation 819 : * 820 : IF ((Ytest+SMT).GE.VARF(3)) THEN 821 : BB0=EXP(-.5D0*((Ytest+SMT-SIGN0)/SNS0)**2.D0) 822 : BB0=BB0*FNS0/(SNS0*((2.D0*XPI)**.5D0)) 823 : ELSE 824 : BB0=0.D0 825 : ENDIF 826 : IF (EPSMT.GE.VARF(4)) THEN 827 : DD0=EXP(-.5D0*((EPSMT-EPSN0)/SNE0)**2.D0) 828 : DD0=DD0*FNE0/(SNE0*((2.D0*XPI)**.5D0)) 829 : ELSE 830 : DD0=0.D0 831 : ENDIF 832 : * 833 : *---> Fonction d'endommagement 834 : * 835 : IF (FT.LE.FC0) THEN 836 : FST0=FT 837 : ddel0=1.D0 838 : ELSE 839 : FST0=FC0+(FU0-FC0)/(FF0-FC0)*(FT-FC0) 840 : ddel0=(FU0-FC0)/(FF0-FC0) 841 : ENDIF 842 : IF (FST0.GT.FU0) FST0=FU0 843 : * 844 : *---> Terme d'endommagement 845 : * 846 : SMTM=MAX(SMT,0.D0) 847 : SMTM=SMT 848 : Gtest=1.D0+(Q0*FST0)*(Q0*FST0)- 849 : . (2.D0*Q0*FST0*COSH(3.D0*SMTM/(2.D0*Ytest))) 850 : * 851 : *---> Fonction de charge test PHIT mise à jour 852 : * 853 : PHIT=(Stest*Stest)-(Ytest*Ytest)*Gtest 854 : PHI0=ABS(PHIT) 855 : * 856 : *---> Cas particulier des matériaux totalement endommagés 857 : * 858 : 92 IF (FT.GE.FF0) THEN 859 : * 860 : FT=FF0+Fmin0 861 : * 862 : SMT=0.D0 863 : treps1=(SMT/(1.D0-FF0)-SMT0/(1.D0-VAR0(2)))/XK0 864 : DO 13 I=1,6 865 : SIGT(I)=0.D0 866 : A=0.D0 867 : IF (I.LE.3) A=1.D0 868 : depel0=SIGT(I)-RSIG0(I) 869 : DEVT(I)=(depel0-((SMT-SMT0)*A))/(2.D0*G0) 870 : depel0=DEVT(I)+(treps1*A/3.D0) 871 : DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0 872 : IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN 873 : DEFT(3)=rdeps3+RDEPS(3)-depel0 874 : ENDIF 875 : 13 CONTINUE 876 : treps1=DEFT(1)+DEFT(2)+DEFT(3) 877 : DO 15 I=1,6 878 : A=0.D0 879 : IF (I.LE.1) A=1.D0 880 : DEVT(I)=DEFT(I)-(DEFT(1)+DEFT(2)+DEFT(3))*A/3.D0 881 : 15 CONTINUE 882 : IF (treps1.LE.0.D0) THEN 883 : DO 17 I=1,3 884 : DEFT(I)=DEVT(I) 885 : 17 CONTINUE 886 : ENDIF 887 : dlam0=(2.D0*PROCON(DEVT,DEVT,6)/3.D0)**.5D0 888 : EPSPT=EPSPT+dlam0 889 : ** RHOT=RMIN0 890 : * 891 : Gtest=0.D0 892 : PHIT=0.D0 893 : PHI0=0.D0 894 : * 895 : ENDIF 896 : * 897 : *---> Test de convergence ou itération suivante 898 : * 899 : iter0=iter0+1 900 : IF (iter0.LT.200) THEN 901 : GOTO 99 902 : ELSE 903 : PHI0=0.D0 904 : KERRE=460 905 : GOTO 99 906 : ENDIF 907 : * 908 : ENDIF 909 : * 910 : * 911 : * Vérification des sous incrémentations 912 : * 913 : IF ((iter2.LT.nmax0).AND.(KERRE.EQ.0)) THEN 914 : GOTO 95 915 : ENDIF 916 : * 917 : * Passage des deformations de cisaillement exprimées 918 : * en déformations aux déformations de cisaillement exprimées 919 : * en GAMA 920 : * 921 : DO 117 I=1,6 922 : A=1.D0 923 : IF (I.GT.3) A=2.D0 924 : RDEFP(I)=RDEFP(I)*A 925 : 117 CONTINUE 926 : * 927 : * 928 : * Passage a l'option de calcul pour les contraintes 929 : *========================================================= 930 : * 931 : IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN 932 : IF (NSTRS .EQ. 6) THEN 933 : * 934 : *---> MASSIF 3D 935 : * 936 : DO 170 I=1,NSTRS 937 : SIGF(I)=RSIGF(I) 938 : DEFP(I)=RDEFP(I) 939 : DSIGT(I)=RDIGT(I) 940 : 170 CONTINUE 941 : ELSE IF ( NSTRS .EQ. 4 ) THEN 942 : * 943 : *---> Calcul axisymétrique ou contraintes planes 944 : * 945 : DO 180 I=1,NSTRS 946 : SIGF(I)=RSIGF(I) 947 : DEFP(I)=RDEFP(I) 948 : DSIGT(I)=RDIGT(I) 949 : 180 CONTINUE 950 : ENDIF 951 : ENDIF 952 : RETURN 953 : * 954 : END 955 : 956 : 957 : 958 : 959 : 960 : 961 : 962 : 963 : 964 : 965 : 966 : 967 : 968 : 969 : 970 :
© Cast3M 2003 - All rights reserved.
Disclaimer