1 : C PRESTO SOURCE CHAT 05/01/13 02:25:58 5004 2 : SUBROUTINE PRESTO(DEPST,NSTRS,MFR1,IB,IGAU, 3 : 1 DSIGT,NCOMAT,SIG0,VAR0,XMAT,XCAR,NVARI,ICARA, 4 : 2 SIGF,VARF,DEFP,TETA1,TETA2,KERRE,DT0) 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 : * TETA1 = TEMPERATURE AU DEBUT DU PAS DE TEMPS 25 : * TETA2 = TEMPERATURE A LA FIN DU PAS DE TEMPS 26 : * DT0 = PAS DE TEMPS DU CALCUL 27 : * 28 : * SORTIE : 29 : * -------- 30 : * 31 : * SIGF(NSTRS) = CONTRAINTES FINALES 32 : * VARF(NVARI) = VARIALES INTERNES A LA FIN DU PAS D'INTEGRATION 33 : * DEFP(NSTRS) = INCREMENT DE DEFORMATION PLASTIQUE A LA FIN 34 : * DU PAS D'INTEGRATION 35 : * ============================================================ 36 : * ICI IL FAUT PROGRAMMER EN FORTRAN PUR 37 : *============================================================= 38 : * 39 : IMPLICIT INTEGER(I-N) 40 : IMPLICIT REAL*8(A-H,O-Z) 41 : * 42 : -INC CCREEL 43 : -INC CCOPTIO 44 : DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*),SIGF(*), 45 : & VARF(*),DEFP(*),DSIGT(*) 46 : DIMENSION RSIG0(6),RDEPS(6),RSIGF(6),RDEFP(6),RDIGT(6) 47 : DIMENSION DEVT(6),ELT(6),SIGT(6),ST(6),DEFT(6) 48 : DIMENSION DIJ(6),CIJ(6),CIJ2(6) 49 : * 50 : * Adaptation de l'option de calcul vers le 3D massif de SIG0 a RSIG0 51 : *==================================================================== 52 : * 53 : IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN 54 : * 55 : *---> 1 formulation massive 56 : *---> 2 formulation quasi incompressible 57 : *---> MASSIF 3D 58 : * 59 : IF (NSTRS .EQ. 6) THEN 60 : DO 110 I=1,NSTRS 61 : RSIG0(I)=SIG0(I) 62 : VARF(I+1)=VAR0(I+1)+DEPST(I) 63 : RDEPS(I)=DEPST(I) 64 : 110 CONTINUE 65 : ELSE IF (NSTRS.EQ.4.AND.((IFOUR.EQ.0) 66 : & .OR.(IFOUR.EQ.-1).OR.(IFOUR.EQ.-2))) THEN 67 : * 68 : *---> Calcul en mode deformations planes ou axisymetrique 69 : * ou en contraintes planes 70 : * 71 : DO 115 I=1,NSTRS 72 : RSIG0(I)=SIG0(I) 73 : VARF(I+1)=VAR0(I+1)+DEPST(I) 74 : RDEPS(I)=DEPST(I) 75 : 115 CONTINUE 76 : RSIG0(5)=0.D0 77 : RSIG0(6)=0.D0 78 : VARF(6)=VAR0(6) 79 : VARF(7)=VAR0(7) 80 : RDEPS(5)=0.D0 81 : RDEPS(6)=0.D0 82 : ENDIF 83 : ELSE 84 : KERRE = 99 85 : GOTO 98 86 : ENDIF 87 : * 88 : * Passage des deformations de cisaillement exprimées 89 : * en GAMA aux déformations de cisaillement exprimées 90 : * en déformations 91 : * 92 : DO 116 I=1,6 93 : A=1.D0 94 : IF (I.GT.3) A=2.D0 95 : RDEPS(I)=RDEPS(I)/A 96 : VAR0(I+1)=VAR0(I+1)/A 97 : VARF(I+1)=VARF(I+1)/A 98 : VAR0(I+8)=VAR0(I+8)/A 99 : VARF(I+8)=VARF(I+8)/A 100 : 116 CONTINUE 101 : * 102 : * Données du materiau 103 : *=========================================================== 104 : * 105 : * 106 : YOUNG=XMAT(1) 107 : XNU=XMAT(2) 108 : G00=YOUNG/(2.D0*(1.D0+XNU)) 109 : IF (IFOUR.EQ.-2) THEN 110 : RHO0=XMAT(21) 111 : ELSE 112 : RHO0=XMAT(20) 113 : ENDIF 114 : TAU0=XMAT(5) 115 : P0=XMAT(6) 116 : S0=XMAT(7) 117 : SINF0=XMAT(8) 118 : XK00=XMAT(9) 119 : G0=XMAT(10) 120 : Y0=XMAT(11) 121 : YINF0=XMAT(12) 122 : Y1=XMAT(13) 123 : Y2=XMAT(14) 124 : BETA0=XMAT(15) 125 : Gp1=XMAT(16) 126 : Gt1=XMAT(17) 127 : XMU0=XMAT(18) 128 : TM0=XMAT(19) 129 : * 130 : * Calculs préliminaires 131 : *=================================== 132 : * 133 : *---> Nombre d'iterations internes maximales 134 : * 135 : NMAX=200 136 : iter00=0 137 : * 138 : *---> Trace de l'incrément de déformation totale 139 : * 140 : IF (IFOUR.EQ.-2) THEN 141 : VARF(4)=VAR0(11)-(XNU/(1.D0-XNU) 142 : . *(VARF(2)-VAR0(9)+VARF(3)-VAR0(10))) 143 : RDEPS(3)=VARF(4)-VAR0(4) 144 : ENDIF 145 : 96 treps1=RDEPS(1)+RDEPS(2)+RDEPS(3) 146 : * 147 : *---> Vitesse de deformation equivalente 148 : * 149 : * DT0=5.D-7 150 : DEPS2=RDEPS(1)*RDEPS(1)+RDEPS(2)*RDEPS(2)+RDEPS(3)*RDEPS(3) 151 : DEPS2=DEPS2-RDEPS(1)*RDEPS(2)-RDEPS(2)*RDEPS(3) 152 : DEPS2=DEPS2-RDEPS(3)*RDEPS(1) 153 : DEPS3=RDEPS(4)*RDEPS(4)+RDEPS(5)*RDEPS(5)+RDEPS(6)*RDEPS(6) 154 : DEPS2=DEPS2+3.D0*DEPS3 155 : DEPS0=(DEPS2**(0.5D0))*2.D0/3.D0 156 : DEPS0=DEPS0+(treps1/3.D0) 157 : IF (ABS(DT0).LT.1.D-15) THEN 158 : DEPS=0.D0 159 : ELSE 160 : DEPS=DEPS0/DT0 161 : END IF 162 : * 163 : *---> Trace de la déformation totale 164 : * 165 : treps0=VARF(2)+VARF(3)+VARF(4) 166 : * 167 : *---> Calcul de la compression a la fin du pas de temps 168 : * 169 : RHO1=RHO0/(1.D0+treps0) 170 : ETAnew=RHO1/RHO0 171 : * 172 : *---> Temperature de fusion a la fin (TM2) 173 : * du pas de temps 174 : * 175 : TM2=TM0*EXP(2.D0*XMU0*(1.D0-(1.D0/ETAnew))) 176 : TM2=TM2/(ETAnew**(2.D0/3.D0)) 177 : * 178 : *---> On divise le coefficient variable avec la température 179 : * par la température de fusion 180 : * 181 : XK0=XK00/TM2 182 : * 183 : * Calculs des grandeurs élastiques test 184 : *==================================================== 185 : * 186 : *---> Déformations 187 : * -élastiques totales test ELT 188 : * -trace des déformations élastiques totales TRELT 189 : * -déviatoires élastiques test DEVT 190 : * -incrément de déformation plastique DEFT 191 : * -déformation plastique équivalente test EPSPT 192 : * 193 : DO 100 I=1,6 194 : ELT(I)=VARF(I+1)-VAR0(8+I) 195 : 100 CONTINUE 196 : TRELT=ELT(1)+ELT(2)+ELT(3) 197 : DO 200 I=1,6 198 : A=0.D0 199 : IF (I.LE.3) A=1.D0 200 : DEVT(I)=ELT(I)-(1.D0*TRELT*A/3.D0) 201 : DEFT(I)=0.D0 202 : 200 CONTINUE 203 : EPSPT=VAR0(1) 204 : * 205 : *---> Coefficient proportionnel entre le module de 206 : * cisaillement et le module de pression hydrostatique 207 : * 208 : E1=2.D0*(1.D0+XNU)/(3.D0*(1.D0-(2.D0*XNU))) 209 : * 210 : *---> Données matériaux Gtest et XKtest 211 : * 212 : pente1=E1*Gp1*TRELT/(ETAnew**(1.D0/3.D0)) 213 : Gtest=(G00+Gt1)/(1.D0+pente1) 214 : XKtest=E1*Gtest 215 : * 216 : *---> Si on dépasse la température de fusion, 217 : * le solide devient liquide 218 : * -les caractéritiques matériaux sont constantes 219 : * -G=0, donc le déviateur élastique est nul 220 : * -la limite d'élasticité est constante et nulle donc: 221 : * -la trace de la déformation plastique est nulle 222 : * 223 : IF (TETA2.GT.TM2) THEN 224 : Gtest=0.D0 225 : DEVT(1)=0.D0 226 : DEVT(2)=0.D0 227 : DEVT(3)=0.D0 228 : DEVT(4)=0.D0 229 : DEVT(5)=0.D0 230 : DEVT(6)=0.D0 231 : TREPT=RDEPS(1)+RDEPS(2)+RDEPS(3) 232 : TRELT=TREPT 233 : DEFT(1)=RDEPS(1)-(1.D0*TREPT/3.D0) 234 : DEFT(2)=RDEPS(2)-(1.D0*TREPT/3.D0) 235 : DEFT(3)=RDEPS(3)-(1.D0*TREPT/3.D0) 236 : DEFT(4)=RDEPS(4) 237 : DEFT(5)=RDEPS(5) 238 : DEFT(6)=RDEPS(6) 239 : XKtest=VAR0(8) 240 : DEPS0=PROCON(DEFT,DEFT,6) 241 : DEPS0=(2.D0*DEPS0/3.D0)**.5D0 242 : EPSPT=EPSPT+DEPS0 243 : ENDIF 244 : * 245 : *---> Contraintes test 246 : * 247 : Ptest=-1.D0*XKtest*TRELT 248 : DO 300 I=1,6 249 : ST(I)=2.D0*Gtest*DEVT(I) 250 : A=0.D0 251 : IF (I.LE.3) A=1.D0 252 : SIGT(I)=ST(I)-(Ptest*A) 253 : 300 CONTINUE 254 : * 255 : *---> Contrainte équivalente test STEST 256 : * 257 : STEST2=(SIGT(1)*SIGT(1))+(SIGT(2)*SIGT(2)) 258 : STEST2=STEST2+(SIGT(3)*SIGT(3)) 259 : STEST2=STEST2-(SIGT(1)*SIGT(2))-(SIGT(2)*SIGT(3)) 260 : STEST2=STEST2-(SIGT(3)*SIGT(1)) 261 : STEST3=(SIGT(4)*SIGT(4))+(SIGT(5)*SIGT(5)) 262 : STEST3=STEST3+(SIGT(6)*SIGT(6)) 263 : STEST2=STEST2+(3.D0*STEST3) 264 : IF (STEST2.LT.0.D0) STEST2=0.D0 265 : STEST=STEST2**(0.5D0) 266 : * 267 : *---> Calcul de la pulsation de Debye: W0 et de DTETA0 268 : * 269 : W0=(Gtest/RHO1)**(.5D0) 270 : * 271 : DTETA0=(1.D0/6.D0)*((4.D0/XPI)**(.5D0))*W0 272 : * 273 : *---> Calcul de la contrainte de saturation YSAT0 274 : * 275 : IF (ABS(DT0).LT.1.D-15) THEN 276 : YS1=SINF0 277 : ELSE 278 : FONC0=XK0*LOG(G0*DTETA0/DEPS) 279 : ERF0=ERGAUS(FONC0) 280 : YS1=S0-((S0-SINF0)*ERF0) 281 : ENDIF 282 : YS2=S0*((DEPS/(G0*DTETA0))**BETA0) 283 : YSAT0=MAX(YS1,YS2) 284 : * 285 : *---> Calcul de la limite élastique YLIM0 286 : * 287 : IF (ABS(DT0).LT.1.D-15) THEN 288 : YLIM1=YINF0 289 : ELSE 290 : FONC0=XK0*LOG(G0*DTETA0/DEPS) 291 : ERF0=ERGAUS(FONC0) 292 : YLIM1=Y0-((Y0-YINF0)*ERF0) 293 : ENDIF 294 : YLIM2=Y1*((DEPS/(G0*DTETA0))**Y2) 295 : YLIM0=MAX(YLIM1,MIN(YLIM2,YS2)) 296 : * 297 : *---> Limite d'élasticité test : CAS CUBIQUE CENTRE C.C. 298 : * 299 : DYSL0=ABS(YSAT0-YLIM0) 300 : DYSL1=ABS(S0-YLIM0) 301 : IF (P0.EQ.0.D0) THEN 302 : * 303 : IF (DYSL0.LT.1.D-15) THEN 304 : Yel0=YSAT0 305 : ELSE 306 : Yel0=-1.D0*TAU0*EPSPT 307 : Yel0=EXP(Yel0/(YSAT0-YLIM0)) 308 : Yel0=YSAT0-((YSAT0-YLIM0)*Yel0) 309 : ENDIF 310 : Ytest=Yel0*Gtest 311 : IF (Ytest.LE.0.D0) Ytest=0.D0 312 : * 313 : ELSE 314 : * 315 : *---> Limite d'élasticité test : AUTRES CAS 316 : * 317 : IF (DYSL0.LT.1.D-15) THEN 318 : Yel0=YLIM0 319 : ELSE IF (DYSL1.LT.1.D-15) THEN 320 : Yel0=YLIM0 321 : ELSE 322 : Yel1=-1.D0*TAU0*EPSPT*P0 323 : denom0=P0*(YSAT0-YLIM0)/(S0-YLIM0) 324 : denom1=EXP(denom0)-1.D0 325 : denom1=denom1*(S0-YLIM0) 326 : Yel1=EXP(Yel1/denom1) 327 : denom2=1.D0-(EXP(-1.D0*denom0)) 328 : Yel0=LOG(1.D0-(denom2*Yel1)) 329 : coeff0=(S0-YLIM0)/P0 330 : Yel0=YSAT0+(coeff0*Yel0) 331 : ENDIF 332 : Ytest=Yel0*Gtest 333 : IF (Ytest.LE.0.D0) Ytest=0.D0 334 : * 335 : ENDIF 336 : * 337 : *---> Fonction de charge test PHIT 338 : * 339 : PHIT=STEST-Ytest 340 : IF (TETA2.GT.TM2) PHIT=0.D0 341 : * 342 : * 343 : * Vérification du critère de plasticité 344 : *========================================================= 345 : * 346 : *---> Erreur admise 347 : * 348 : ERR0=1.D-7*ABS(PHIT) 349 : * 350 : *---> Critère de plasticité 351 : * 352 : PHI0=PHIT 353 : iter0=0 354 : * 355 : * 356 : 97 IF (PHI0.LE.ERR0) THEN 357 : * 358 : * On est élastique 359 : *========================================================= 360 : * 361 : *-------------------------------------------------------------- 362 : * Cas particulier des contraintes planes 363 : *-------------------------------------------------------------- 364 : * 365 : IF (IFOUR.EQ.-2) THEN 366 : * 367 : *---> On vérifie qu'on est en contraintes planes 368 : * 369 : iter00=iter00+1 370 : SIG3=ABS(SIGT(3)) 371 : IF (SIG3.GT.1.D-2) THEN 372 : ELT(3)=-1.D0*XNU*(ELT(1)+ELT(2))/(1.D0-XNU) 373 : VARF(4)=VAR0(11)+ELT(3)+DEFT(3) 374 : RDEPS(3)=VARF(4)-VAR0(4) 375 : IF (iter00.LE.NMAX) THEN 376 : GOTO 96 377 : ELSE 378 : KERRE=460 379 : GOTO 98 380 : ENDIF 381 : ENDIF 382 : ENDIF 383 : * 384 : *----------------------------------------------------- 385 : * Fin du cas particulier des contraintes planes 386 : *----------------------------------------------------- 387 : * 388 : DO 400 I=1,6 389 : RSIGF(I)=SIGT(I) 390 : RDEFP(I)=DEFT(I) 391 : RDIGT(I)=RSIGF(I)-RSIG0(I) 392 : VARF(I+8)=DEFT(I)+VAR0(I+8) 393 : 400 CONTINUE 394 : VARF(1)=EPSPT 395 : VARF(8)=XKtest 396 : DEPST(3)=VARF(4)-VAR0(4) 397 : * 398 : * 399 : ELSE 400 : * 401 : * On plastifie 402 : *========================================================= 403 : * 404 : *---> Coefficient donnant la trace de l'incrément de 405 : * déformation plastique 406 : * 407 : tr0=Yel0*Gp1/(ETAnew**(1.D0/3.D0)) 408 : * 409 : *---> Coefficient donnant la variation de pression due 410 : * à l'incrément de déformation plastique 411 : * 412 : coeff1=E1*Gp1*TRELT/(ETAnew**(1.D0/3.D0)) 413 : coeff1=XKtest*tr0/(1.D0+coeff1) 414 : * 415 : *---> Module d'écrouissage H = - d PHI / d EPSeq 416 : * 417 : IF (P0.EQ.0.D0) THEN 418 : IF (DYSL0.LT.1.D-15) THEN 419 : H0=0.D0 420 : ELSE 421 : H0=TAU0*EXP(-1.D0*TAU0*EPSPT/(YSAT0-YLIM0)) 422 : H0=H0*Gtest 423 : ENDIF 424 : ELSE 425 : IF ((DYSL0.LT.1.D-15).OR.(DYSL1.LT.1.D-15)) THEN 426 : H0=0.D0 427 : ELSE 428 : H0=coeff0*denom2*P0*TAU0 429 : fac0=denom2*EXP(-1.D0*P0*TAU0*EPSPT/denom1) 430 : H0=H0/(denom1*(1.D0-fac0)) 431 : H0=H0*Gtest 432 : ENDIF 433 : ENDIF 434 : * 435 : *---> Matrice de la direction d'écoulement DIJ = d PHI / d SIG 436 : * 437 : DO 500 I=1,6 438 : A=0.D0 439 : IF (I.LE.3) A=1.D0 440 : DIJ(I)=3.D0*ST(I)/(STEST*2.D0) 441 : DIJ(I)=DIJ(I)+(A*tr0/3.D0) 442 : 500 CONTINUE 443 : * 444 : *---> Matrice de la variation de contrainte due à l'incrément de 445 : * déformation plastique CIJ = - ( Hooke . d PHI / d SIG) 446 : * 447 : DO 600 I=1,6 448 : A=0.D0 449 : IF (I.LE.3) A=1.D0 450 : CIJ(I)=3.D0*Gtest*ST(I)/STEST 451 : CIJ(I)=CIJ(I)-(2.D0*Gp1*coeff1*DEVT(I)/(ETAnew**(1.D0/3.D0)) 452 : $ ) 453 : CIJ(I)=CIJ(I)+(coeff1*A) 454 : CIJ2(I)=3.D0*Gp1*coeff1*ST(I)/((ETAnew**(1.D0/3.D0))*STEST) 455 : 600 CONTINUE 456 : * 457 : *---> Calcul du produit DIJ : CIJ 458 : * 459 : PROD0=PROCON(DIJ,CIJ,6) 460 : * 461 : *---> Calcul de DIJ : CIJ2 462 : * 463 : PROD1=PROCON(DIJ,CIJ2,6) 464 : PROD11=ABS(PROD1) 465 : * 466 : *---> Calcul de l'incrément du paramètre plastique 467 : * On résoud une équation du second degré 468 : * Si le terme (PROD0+H0) est positif, dlam0 est 469 : * clairement défini par la solution positive de l'équation 470 : * Sinon, on prend la solution qui est de meme signe que 471 : * PHIT (une solution est nécessairement 472 : * de signe opposé a PHIT). Si aucune solution ne convient 473 : * on prend celle qui a le plus petit module afin d'éviter de 474 : * diverger. 475 : * 476 : delta0=(PROD0+H0)*(PROD0+H0) 477 : delta0=delta0+(4.D0*PHIT*PROD1) 478 : IF ((delta0.LT.0.D0).OR.(PROD11.LT.1.D-15)) THEN 479 : dlam0=PHIT/(PROD0+H0) 480 : ELSE 481 : IF ((PROD0+H0).GE.0.D0) THEN 482 : dlam0=(delta0**(.5D0))/(2.D0*PROD1) 483 : dlam0=dlam0-((PROD0+H0)/(2.D0*PROD1)) 484 : ELSE 485 : dlam0=(delta0**(.5D0))/(2.D0*PROD1) 486 : dlam0=dlam0-((PROD0+H0)/(2.D0*PROD1)) 487 : IF ((dlam0*PHIT).LE.0.D0) THEN 488 : dlam1=-1.D0*(delta0**(.5D0))/(2.D0*PROD1) 489 : dlam1=dlam1-((PROD0+H0)/(2.D0*PROD1)) 490 : IF (ABS(dlam1).LT.ABS(dlam0)) dlam0=dlam1 491 : ENDIF 492 : ENDIF 493 : ENDIF 494 : * 495 : * 496 : * Mise à jour des grandeurs après écoulement plastique 497 : *============================================================ 498 : * 499 : *---> Déformation plastique équivalente 500 : * 501 : EPSPT=EPSPT+dlam0 502 : * 503 : *---> Pression mise à jour 504 : * 505 : Ptest=Ptest+(coeff1*dlam0) 506 : * 507 : *---> Contraintes 508 : * 509 : DO 800 I=1,6 510 : A=0.D0 511 : IF (I.LE.3) A=1.D0 512 : SIGT(I)=SIGT(I)-(CIJ(I)*dlam0)-(CIJ2(I)*dlam0*dlam0) 513 : * SIGT(I)=SIGT(I)-(CIJ(I)*dlam0) 514 : ST(I)=SIGT(I)+(Ptest*A) 515 : 800 CONTINUE 516 : * 517 : *---> Contrainte équivalente mise à jour STEST 518 : * 519 : STEST2=(SIGT(1)*SIGT(1))+(SIGT(2)*SIGT(2)) 520 : STEST2=STEST2+(SIGT(3)*SIGT(3)) 521 : STEST2=STEST2-(SIGT(1)*SIGT(2))-(SIGT(2)*SIGT(3)) 522 : STEST2=STEST2-(SIGT(3)*SIGT(1)) 523 : STEST3=(SIGT(4)*SIGT(4))+(SIGT(5)*SIGT(5)) 524 : STEST3=STEST3+(SIGT(6)*SIGT(6)) 525 : STEST2=STEST2+(3.D0*STEST3) 526 : IF (STEST2.LE.0.D0) STEST2=0.D0 527 : STEST=STEST2**(0.5D0) 528 : * 529 : *---> Calcul des données matériaux mises à jour 530 : * 531 : Gtest=G00+(Gp1*Ptest/(ETAnew**(1.D0/3.D0)))+Gt1 532 : XKtest=E1*Gtest 533 : * 534 : * 535 : * Mise à jour des déformations après écoulement plastique 536 : *============================================================ 537 : * 538 : *---> Trace des déformations élastiques mises à jour 539 : * 540 : TRELT=-1.D0*Ptest/XKtest 541 : * 542 : *---> Déviateur des déformations élastiques 543 : * 544 : DO 900 I=1,6 545 : A=0.D0 546 : IF (I.LE.3) A=1.D0 547 : DEVT(I)=ST(I)/(2.D0*Gtest) 548 : IF (IFOUR.EQ.-2) ELT(I)=DEVT(I)+(A*TRELT/3.D0) 549 : 900 CONTINUE 550 : * 551 : *---> Incrément de déformation plastique au cours du pas 552 : * 553 : DO 700 I=1,6 554 : A=0.D0 555 : IF (I.LE.3) A=1.D0 556 : DEFT(I)=VARF(I+1)-VAR0(I+8) 557 : DEFT(I)=DEFT(I)-DEVT(I)-(1.D0*TRELT*A/3.D0) 558 : 700 CONTINUE 559 : * 560 : *---> Mise à jour de la pulsation de Debye: W0 et de DTETA0 561 : * 562 : W0=(Gtest/RHO1)**(.5D0) 563 : * 564 : DTETA0=(1.D0/6.D0)*((4.D0/XPI)**(.5D0))*W0 565 : * 566 : *---> Mise à jour de la contrainte de saturation YSAT0 567 : * 568 : IF (ABS(DT0).LT.1.D-15) THEN 569 : YS1=SINF0 570 : ELSE 571 : ERF0=ERGAUS(FONC0) 572 : YS1=S0-((S0-SINF0)*ERF0) 573 : ENDIF 574 : YS2=S0*((DEPS/(G0*DTETA0))**BETA0) 575 : YSAT0=MAX(YS1,YS2) 576 : * 577 : *---> Mise à jour de la limite élastique YLIM0 578 : * 579 : IF (ABS(DT0).LT.1.D-15) THEN 580 : YLIM1=YINF0 581 : ELSE 582 : FONC0=XK0*LOG(G0*DTETA0/DEPS) 583 : ERF0=ERGAUS(FONC0) 584 : YLIM1=Y0-((Y0-YINF0)*ERF0) 585 : ENDIF 586 : YLIM2=Y1*((DEPS/(G0*DTETA0))**Y2) 587 : YLIM0=MAX(YLIM1,MIN(YLIM2,YS2)) 588 : * 589 : *---> limite d'élasticité mise à jour 590 : * 591 : DYSL0=ABS(YSAT0-YLIM0) 592 : DYSL1=ABS(S0-YLIM0) 593 : IF (P0.EQ.0.D0) THEN 594 : IF (DYSL0.LT.1.D-15) THEN 595 : Yel0=YSAT0 596 : ELSE 597 : Yel0=-1.D0*TAU0*EPSPT 598 : Yel0=EXP(Yel0/(YSAT0-YLIM0)) 599 : Yel0=YSAT0-((YSAT0-YLIM0)*Yel0) 600 : ENDIF 601 : Ytest=Yel0*Gtest 602 : IF (Ytest.LE.0.D0) Ytest=0.D0 603 : ELSE 604 : IF (DYSL0.LT.1.D-15) THEN 605 : Yel0=YLIM0 606 : ELSE IF (DYSL1.LT.1.D-15) THEN 607 : Yel0=YLIM0 608 : ELSE 609 : Yel1=-1.D0*TAU0*EPSPT*P0 610 : denom0=P0*(YSAT0-YLIM0)/(S0-YLIM0) 611 : denom1=EXP(denom0)-1.D0 612 : denom1=denom1*(S0-YLIM0) 613 : Yel1=EXP(Yel1/denom1) 614 : denom2=1.D0-(EXP(-1.D0*denom0)) 615 : Yel0=LOG(1.D0-(denom2*Yel1)) 616 : coeff0=(S0-YLIM0)/P0 617 : Yel0=YSAT0+(coeff0*Yel0) 618 : ENDIF 619 : Ytest=Yel0*Gtest 620 : ENDIF 621 : * 622 : *---> Nouvelle fonction de charge 623 : * 624 : PHIT=STEST-Ytest 625 : PHI0=ABS(PHIT) 626 : * 627 : *---> Test de convergence ou itération suivante 628 : * 629 : iter0=iter0+1 630 : * 631 : IF (iter0.LT.NMAX) THEN 632 : GOTO 97 633 : ELSE 634 : KERRE=460 635 : GOTO 98 636 : ENDIF 637 : * 638 : ENDIF 639 : * 640 : * Passage des deformations de cisaillement exprimées 641 : * en déformations aux déformations de cisaillement exprimées 642 : * en GAMA 643 : * 644 : DO 117 I=1,6 645 : A=1.D0 646 : IF (I.GT.3) A=2.D0 647 : RDEFP(I)=RDEFP(I)*A 648 : VAR0(I+1)=VAR0(I+1)*A 649 : VARF(I+1)=VARF(I+1)*A 650 : VAR0(I+8)=VAR0(I+8)*A 651 : VARF(I+8)=VARF(I+8)*A 652 : 117 CONTINUE 653 : * 654 : * Passage a l'option de calcul pour les contraintes 655 : *========================================================= 656 : * 657 : IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN 658 : IF (NSTRS .EQ. 6) THEN 659 : * 660 : *---> MASSIF 3D 661 : * 662 : DO 170 I=1,NSTRS 663 : SIGF(I)=RSIGF(I) 664 : DEFP(I)=RDEFP(I) 665 : DSIGT(I)=RDIGT(I) 666 : 170 CONTINUE 667 : ELSE IF ( NSTRS .EQ. 4 ) THEN 668 : * 669 : *---> Calcul axisymétrique ou contraintes planes 670 : * 671 : DO 180 I=1,NSTRS 672 : SIGF(I)=RSIGF(I) 673 : DEFP(I)=RDEFP(I) 674 : DSIGT(I)=RDIGT(I) 675 : 180 CONTINUE 676 : ENDIF 677 : ENDIF 678 : * 679 : * 680 : 98 CONTINUE 681 : RETURN 682 : * 683 : END 684 : 685 : 686 : 687 : 688 : 689 : 690 : 691 : 692 : 693 : 694 : 695 : 696 :
© Cast3M 2003 - All rights reserved.
Disclaimer