1 : C STEINB SOURCE CHAT 05/01/13 03:25:08 5004 2 : SUBROUTINE STEINB(DEPST,NSTRS,MFR1,IB,IGAU, 3 : 1 DSIGT,NCOMAT,SIG0,VAR0,XMAT,XCAR,NVARI,ICARA, 4 : 2 SIGF,VARF,DEFP,TETA1,TETA2,KERRE) 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 : * 27 : * SORTIE : 28 : * -------- 29 : * 30 : * SIGF(NSTRS) = CONTRAINTES FINALES 31 : * VARF(NVARI) = VARIALES INTERNES A LA FIN DU PAS D'INTEGRATION 32 : * DEFP(NSTRS) = INCREMENT DE DEFORMATION PLASTIQUE A LA FIN 33 : * DU PAS D'INTEGRATION 34 : * ============================================================ 35 : * ICI IL FAUT PROGRAMMER EN FORTRAN PUR 36 : *============================================================= 37 : * 38 : IMPLICIT INTEGER(I-N) 39 : IMPLICIT REAL*8(A-H,O-Z) 40 : * 41 : -INC CCOPTIO 42 : DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*),SIGF(*), 43 : & VARF(*),DEFP(*),DSIGT(*) 44 : DIMENSION RSIG0(6),RDEPS(6),RSIGF(6),RDEFP(6),RDIGT(6) 45 : DIMENSION DEVT(6),ELT(6),SIGT(6),ST(6),DEFT(6) 46 : DIMENSION DIJ(6),CIJ(6),CIJ2(6) 47 : * 48 : * Adaptation de l'option de calcul vers le 3D massif de SIG0 a RSIG0 49 : *==================================================================== 50 : * 51 : IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN 52 : * 53 : *---> 1 formulation massive 54 : *---> 2 formulation quasi incompressible 55 : *---> MASSIF 3D 56 : * 57 : IF (NSTRS .EQ. 6) THEN 58 : DO 110 I=1,NSTRS 59 : RSIG0(I)=SIG0(I) 60 : VARF(I+1)=VAR0(I+1)+DEPST(I) 61 : RDEPS(I)=DEPST(I) 62 : 110 CONTINUE 63 : ELSE IF ( NSTRS .EQ. 4 .AND. ((IFOUR .EQ. 0) 64 : & .OR.(IFOUR.EQ.-1).OR.(IFOUR.EQ.-2))) THEN 65 : * 66 : *---> Calcul en mode deformations planes ou axisymetrique 67 : * 68 : DO 115 I=1,NSTRS 69 : RSIG0(I)=SIG0(I) 70 : VARF(I+1)=VAR0(I+1)+DEPST(I) 71 : RDEPS(I)=DEPST(I) 72 : 115 CONTINUE 73 : RSIG0(5)=0.D0 74 : RSIG0(6)=0.D0 75 : VARF(6)=VAR0(6) 76 : VARF(7)=VAR0(7) 77 : RDEPS(5)=0.D0 78 : RDEPS(6)=0.D0 79 : ENDIF 80 : ELSE 81 : KERRE = 99 82 : RETURN 83 : ENDIF 84 : * 85 : * Passage des deformations de cisaillement exprimées 86 : * en GAMA aux déformations de cisaillement exprimées 87 : * en déformations 88 : * 89 : DO 116 I=1,6 90 : A=1.D0 91 : IF (I.GT.3) A=2.D0 92 : RDEPS(I)=RDEPS(I)/A 93 : VAR0(I+1)=VAR0(I+1)/A 94 : VARF(I+1)=VARF(I+1)/A 95 : VAR0(I+8)=VAR0(I+8)/A 96 : VARF(I+8)=VARF(I+8)/A 97 : 116 CONTINUE 98 : * 99 : * Données du materiau 100 : *=========================================================== 101 : * 102 : YOUNG=XMAT(1) 103 : XNU=XMAT(2) 104 : Y0=XMAT(5) 105 : BETA=XMAT(6) 106 : Xm=XMAT(7) 107 : EPSI=XMAT(8) 108 : G0=YOUNG/(2.D0*(1.D0+XNU)) 109 : Gp1=XMAT(9) 110 : Gt1=XMAT(10) 111 : YMAX=XMAT(11) 112 : TM0=XMAT(12) 113 : XMU0=XMAT(13) 114 : * 115 : * Calculs préliminaires 116 : *=================================== 117 : * 118 : *---> Nombre maximal d'itérations 119 : * 120 : NMAX=200 121 : iter00=0 122 : * 123 : *---> Trace de la déformation totale 124 : * 125 : IF (IFOUR.EQ.-2) THEN 126 : VARF(4)=VAR0(11)-(XNU/(1.D0-XNU) 127 : . *(VARF(2)-VAR0(9)+VARF(3)-VAR0(10))) 128 : ENDIF 129 : 98 treps0=VARF(2)+VARF(3)+VARF(4) 130 : * 131 : *---> Calcul de la compression finale ETAnew 132 : * 133 : ETAnew=1.D0/(1.D0+treps0) 134 : * 135 : *---> Temperature de fusion a la fin (TM2) 136 : * du pas de temps 137 : * 138 : TM2=TM0*EXP(2.D0*XMU0*(1.D0-(1.D0/ETAnew))) 139 : TM2=TM2/(ETAnew**(2.D0/3.D0)) 140 : * 141 : * 142 : * Calculs des grandeurs élastiques test 143 : *==================================================== 144 : * 145 : *---> Déformations 146 : * -élastiques totales test ELT 147 : * -trace des déformations élastiques totales TRELT 148 : * -déviatoires élastiques test DEVT 149 : * -incrément de déformation plastique DEFT 150 : * -déformation plastique équivalente test EPSPT 151 : * 152 : DO 100 I=1,6 153 : ELT(I)=VARF(I+1)-VAR0(8+I) 154 : 100 CONTINUE 155 : * IF (IFOUR.EQ.-2) 156 : * . ELT(3)=-1.D0*XNU*(ELT(1)+ELT(2))/(1.D0-XNU) 157 : TRELT=ELT(1)+ELT(2)+ELT(3) 158 : DO 200 I=1,6 159 : A=0.D0 160 : IF (I.LE.3) A=1.D0 161 : DEVT(I)=ELT(I)-(1.D0*TRELT*A/3.D0) 162 : DEFT(I)=0.D0 163 : 200 CONTINUE 164 : EPSPT=VAR0(1) 165 : * 166 : *---> Coefficient proportionnel entre le module de 167 : * cisaillement et le module de pression hydrostatique 168 : * 169 : E1=2.D0*(1.D0+XNU)/(3.D0*(1.D0-(2.D0*XNU))) 170 : * 171 : *---> Données matériaux Gtest et XKtest 172 : * 173 : pente1=E1*Gp1*TRELT/(ETAnew**(1.D0/3.D0)) 174 : Gtest=(G0+Gt1)/(1.D0+pente1) 175 : XKtest=E1*Gtest 176 : * 177 : *---> Si on dépasse la température de fusion, 178 : * le solide devient liquide 179 : * -les caractéritiques matériaux sont constantes 180 : * -G=0, donc le déviateur élastique est nul 181 : * -la limite d'élasticité est constante et nulle donc: 182 : * -la trace de la déformation plastique est nulle 183 : * 184 : IF (TETA2.GT.TM2) THEN 185 : Gtest=0.D0 186 : DEVT(1)=0.D0 187 : DEVT(2)=0.D0 188 : DEVT(3)=0.D0 189 : DEVT(4)=0.D0 190 : DEVT(5)=0.D0 191 : DEVT(6)=0.D0 192 : TREPT=RDEPS(1)+RDEPS(2)+RDEPS(3) 193 : TRELT=TREPT 194 : DEFT(1)=RDEPS(1)-(1.D0*TREPT/3.D0) 195 : DEFT(2)=RDEPS(2)-(1.D0*TREPT/3.D0) 196 : DEFT(3)=RDEPS(3)-(1.D0*TREPT/3.D0) 197 : DEFT(4)=RDEPS(4) 198 : DEFT(5)=RDEPS(5) 199 : DEFT(6)=RDEPS(6) 200 : XKtest=VAR0(8) 201 : DEPS=PROCON(DEFT,DEFT,6) 202 : DEPS=(2.D0*DEPS/3.D0)**.5D0 203 : EPSPT=EPSPT+DEPS 204 : IF (IFOUR.EQ.-2) THEN 205 : KERRE=99 206 : RETURN 207 : ENDIF 208 : ENDIF 209 : * 210 : *---> Contraintes test 211 : * 212 : Ptest=-1.D0*XKtest*TRELT 213 : DO 300 I=1,6 214 : ST(I)=2.D0*Gtest*DEVT(I) 215 : A=0.D0 216 : IF (I.LE.3) A=1.D0 217 : SIGT(I)=ST(I)-(Ptest*A) 218 : 300 CONTINUE 219 : * 220 : *---> Contrainte équivalente test STEST 221 : * 222 : STEST2=(SIGT(1)*SIGT(1))+(SIGT(2)*SIGT(2)) 223 : STEST2=STEST2+(SIGT(3)*SIGT(3)) 224 : STEST2=STEST2-(SIGT(1)*SIGT(2))-(SIGT(2)*SIGT(3)) 225 : STEST2=STEST2-(SIGT(3)*SIGT(1)) 226 : STEST3=(SIGT(4)*SIGT(4))+(SIGT(5)*SIGT(5)) 227 : STEST3=STEST3+(SIGT(6)*SIGT(6)) 228 : STEST2=STEST2+(3.D0*STEST3) 229 : IF (STEST2.LT.0.D0) STEST2=0.D0 230 : STEST=(STEST2)**(0.5D0) 231 : * 232 : *---> Limite d'élasticité test 233 : * 234 : Ytest=Y0*((1.D0+(BETA*(EPSPT+EPSI)))**Xm) 235 : IF (Ytest.GT.YMAX) Ytest=YMAX 236 : Ytest=Ytest*Gtest/G0 237 : IF (Ytest.LE.0.D0) Ytest=0.D0 238 : * 239 : *---> Fonction de charge test PHIT 240 : * 241 : PHIT=STEST-Ytest 242 : IF (TETA2.GT.TM2) PHIT=0.D0 243 : * 244 : * 245 : * Vérification du critère de plasticité 246 : *========================================================= 247 : * 248 : *---> Erreur admise 249 : * 250 : ERR0=1.D-7*ABS(PHIT) 251 : ERR0=MAX(ERR0,1.D-8*Y0) 252 : * 253 : *---> Critère de plasticité 254 : * 255 : PHI0=PHIT 256 : iter0=0 257 : * 258 : 99 IF (PHI0.LE.ERR0) THEN 259 : * 260 : * On est élastique 261 : *========================================================= 262 : * 263 : *-------------------------------------------------------------- 264 : * Cas particulier des contraintes planes 265 : *-------------------------------------------------------------- 266 : * 267 : IF (IFOUR.EQ.-2) THEN 268 : * 269 : *---> On vérifie qu'on est en contraintes planes 270 : * 271 : iter00=iter00+1 272 : SIG3=ABS(SIGT(3)) 273 : IF (SIG3.GT.1.D-2) THEN 274 : ELT(3)=-1.D0*XNU*(ELT(1)+ELT(2))/(1.D0-XNU) 275 : VARF(4)=VAR0(11)+ELT(3)+DEFT(3) 276 : IF (iter00.LE.NMAX)THEN 277 : GOTO 98 278 : ELSE 279 : KERRE=460 280 : ENDIF 281 : ENDIF 282 : ENDIF 283 : * 284 : *----------------------------------------------------- 285 : * Fin du cas particulier des contraintes planes 286 : *----------------------------------------------------- 287 : * 288 : DO 400 I=1,6 289 : RSIGF(I)=SIGT(I) 290 : RDEFP(I)=DEFT(I) 291 : RDIGT(I)=RSIGF(I)-RSIG0(I) 292 : VARF(I+8)=DEFT(I)+VAR0(I+8) 293 : 400 CONTINUE 294 : VARF(1)=EPSPT 295 : VARF(8)=XKtest 296 : DEPST(3)=VARF(4)-VAR0(4) 297 : * 298 : * 299 : ELSE 300 : * 301 : * On plastifie 302 : *========================================================= 303 : * 304 : *---> Coefficient donnant la trace de l'incrément de 305 : * déformation plastique 306 : * 307 : tr0=(1.D0+(BETA*(EPSPT+EPSI)))**Xm 308 : tr0=tr0*Y0*Gp1/(G0*(ETAnew**(1.D0/3.D0))) 309 : * 310 : *---> Coefficient donnant la variation de pression due 311 : * à l'incrément de déformation plastique 312 : * 313 : coeff0=E1*Gp1*TRELT/(ETAnew**(1.D0/3.D0)) 314 : coeff0=XKtest*tr0/(1.D0+coeff0) 315 : * 316 : *---> Module d'écrouissage H = - d PHI / d EPSeq 317 : * 318 : H0=Y0*Xm*BETA*((1.D0+(BETA*(EPSPT+EPSI)))**(Xm-1.D0)) 319 : H0=H0*Gtest/G0 320 : * 321 : *---> Matrice de la direction d'écoulement DIJ = d PHI / d SIG 322 : * 323 : DO 500 I=1,6 324 : A=0.D0 325 : IF (I.LE.3) A=1.D0 326 : DIJ(I)=3.D0*ST(I)/(STEST*2.D0) 327 : DIJ(I)=DIJ(I)+(A*tr0/3.D0) 328 : 500 CONTINUE 329 : * 330 : *---> Matrice de la variation de contrainte due à l'incrément de 331 : * déformation plastique CIJ = - ( Hooke . d PHI / d SIG) 332 : * 333 : DO 600 I=1,6 334 : A=0.D0 335 : IF (I.LE.3) A=1.D0 336 : CIJ(I)=3.D0*Gtest*ST(I)/STEST 337 : CIJ(I)=CIJ(I)-(2.D0*Gp1*coeff0*DEVT(I)/(ETAnew**(1.D0/3.D0))) 338 : CIJ(I)=CIJ(I)+(coeff0*A) 339 : CIJ2(I)=3.D0*Gp1*coeff0*ST(I)/((ETAnew**(1.D0/3.D0))*STEST) 340 : 600 CONTINUE 341 : * 342 : *---> Calcul du produit DIJ : CIJ 343 : * 344 : PROD0=PROCON(DIJ,CIJ,6) 345 : * 346 : *---> Calcul de DIJ : CIJ2 347 : * 348 : PROD1=PROCON(DIJ,CIJ2,6) 349 : PROD11=ABS(PROD1) 350 : * 351 : *---> Calcul de l'incrément du paramètre plastique 352 : * On résoud une équation du second degré 353 : * Si le terme (PROD0+H0) est positif, dlam0 est 354 : * clairement défini par la solution positive de l'équation 355 : * Sinon, on prend la solution qui est de meme signe que 356 : * PHIT (une solution est nécessairement 357 : * de signe opposé a PHIT). Si aucune solution ne convient 358 : * on prend celle qui a le plus petit module afin d'éviter de 359 : * diverger. 360 : * 361 : delta0=(PROD0+H0)*(PROD0+H0) 362 : delta0=delta0+(4.D0*PHIT*PROD1) 363 : IF ((delta0.LT.0.D0).OR.(PROD11.LT.1.D-15)) THEN 364 : dlam0=PHIT/(PROD0+H0) 365 : ELSE 366 : IF ((PROD0+H0).GE.0.D0) THEN 367 : dlam0=(delta0**(.5D0))/(2.D0*PROD1) 368 : dlam0=dlam0-((PROD0+H0)/(2.D0*PROD1)) 369 : ELSE 370 : dlam0=(delta0**(.5D0))/(2.D0*PROD1) 371 : dlam0=dlam0-((PROD0+H0)/(2.D0*PROD1)) 372 : IF ((dlam0*PHIT).LE.0.D0) THEN 373 : dlam1=-1.D0*(delta0**(.5D0))/(2.D0*PROD1) 374 : dlam1=dlam1-((PROD0+H0)/(2.D0*PROD1)) 375 : IF (ABS(dlam1).LT.ABS(dlam0)) dlam0=dlam1 376 : ENDIF 377 : ENDIF 378 : ENDIF 379 : * 380 : * 381 : * Mise à jour des grandeurs après écoulement plastique 382 : *============================================================ 383 : * 384 : *---> Déformation plastique équivalente 385 : * 386 : EPSPT=EPSPT+dlam0 387 : * 388 : *---> Pression mise à jour 389 : * 390 : Ptest=Ptest+(coeff0*dlam0) 391 : * 392 : *---> Contraintes 393 : * 394 : DO 800 I=1,6 395 : A=0.D0 396 : IF (I.LE.3) A=1.D0 397 : SIGT(I)=SIGT(I)-(CIJ(I)*dlam0)-(CIJ2(I)*dlam0*dlam0) 398 : ST(I)=SIGT(I)+(Ptest*A) 399 : 800 CONTINUE 400 : * 401 : *---> Contrainte équivalente mise à jour STEST 402 : * 403 : STEST2=(SIGT(1)*SIGT(1))+(SIGT(2)*SIGT(2)) 404 : STEST2=STEST2+(SIGT(3)*SIGT(3)) 405 : STEST2=STEST2-(SIGT(1)*SIGT(2))-(SIGT(2)*SIGT(3)) 406 : STEST2=STEST2-(SIGT(3)*SIGT(1)) 407 : STEST3=(SIGT(4)*SIGT(4))+(SIGT(5)*SIGT(5)) 408 : STEST3=STEST3+(SIGT(6)*SIGT(6)) 409 : STEST2=STEST2+(3.D0*STEST3) 410 : IF (STEST2.LE.0.D0) STEST2=0.D0 411 : STEST=(STEST2)**(0.5D0) 412 : * 413 : *---> Calcul des données matériaux mises à jour 414 : * 415 : Gtest=G0+(Gp1*Ptest/(ETAnew**(1.D0/3.D0)))+Gt1 416 : XKtest=E1*Gtest 417 : * 418 : * 419 : * Mise à jour des déformations après écoulement plastique 420 : *============================================================ 421 : * 422 : *---> Trace des déformations élastiques mises à jour 423 : * 424 : TRELT=-1.D0*Ptest/XKtest 425 : * 426 : *---> Déviateur des déformations élastiques 427 : * 428 : DO 900 I=1,6 429 : A=0.D0 430 : IF (A.LE.1) A=1.D0 431 : DEVT(I)=ST(I)/(2.D0*Gtest) 432 : IF (IFOUR.EQ.-2) ELT(I)=DEVT(I)+(TRELT*A/3.D0) 433 : 900 CONTINUE 434 : * 435 : *---> Incrément de déformation plastique au cours du pas 436 : * 437 : DO 700 I=1,6 438 : A=0.D0 439 : IF (I.LE.3) A=1.D0 440 : DEFT(I)=VARF(I+1)-VAR0(I+8) 441 : DEFT(I)=DEFT(I)-DEVT(I)-(1.D0*TRELT*A/3.D0) 442 : 700 CONTINUE 443 : * 444 : *---> limite d'élasticité mise à jour 445 : * 446 : Ytest=Y0*((1.D0+(BETA*(EPSPT+EPSI)))**Xm) 447 : IF (Ytest.GT.YMAX) Ytest=YMAX 448 : Ytest=Ytest*Gtest/G0 449 : * 450 : *---> Nouvelle fonction de charge 451 : * 452 : PHIT=STEST-Ytest 453 : PHI0=ABS(PHIT) 454 : * 455 : *---> Test de convergence ou itération suivante 456 : * 457 : iter0=iter0+1 458 : IF (iter0.LT.NMAX) THEN 459 : GOTO 99 460 : ELSE 461 : KERRE=460 462 : ENDIF 463 : * 464 : ENDIF 465 : * 466 : * 467 : * Passage des deformations de cisaillement exprimées 468 : * en déformations aux déformations de cisaillement exprimées 469 : * en GAMA 470 : * 471 : DO 117 I=1,6 472 : A=1.D0 473 : IF (I.GT.3) A=2.D0 474 : RDEFP(I)=RDEFP(I)*A 475 : VAR0(I+1)=VAR0(I+1)*A 476 : VARF(I+1)=VARF(I+1)*A 477 : VAR0(I+8)=VAR0(I+8)*A 478 : VARF(I+8)=VARF(I+8)*A 479 : 117 CONTINUE 480 : * 481 : * Passage a l'option de calcul pour les contraintes 482 : *========================================================= 483 : * 484 : IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN 485 : IF (NSTRS .EQ. 6) THEN 486 : * 487 : *---> MASSIF 3D 488 : * 489 : DO 170 I=1,NSTRS 490 : SIGF(I)=RSIGF(I) 491 : DEFP(I)=RDEFP(I) 492 : DSIGT(I)=RDIGT(I) 493 : 170 CONTINUE 494 : ELSE IF ( NSTRS .EQ. 4 ) THEN 495 : * 496 : *---> Calcul axisymétrique ou contraintes planes 497 : * 498 : DO 180 I=1,NSTRS 499 : SIGF(I)=RSIGF(I) 500 : DEFP(I)=RDEFP(I) 501 : DSIGT(I)=RDIGT(I) 502 : 180 CONTINUE 503 : ENDIF 504 : ENDIF 505 : * 506 : * 507 : RETURN 508 : * 509 : END 510 : 511 : 512 : 513 : 514 : 515 : 516 : 517 :
© Cast3M 2003 - All rights reserved.
Disclaimer