1 : C PPLAST SOURCE PV 11/03/07 21:17:49 6885 2 : 3 : SUBROUTINE PPLAST(XMAT,NMATT,VAR0,VARF,NVARI,SIG0, 4 : & SIGF,DEPST,NSTRS,KERRE) 5 : 6 : *************************************** 7 : * Routine de traitement en plasticité * 8 : *************************************** 9 : 10 : * Entrées / Sorties 11 : 12 : * Entrées 13 : 14 : INTEGER NMATT, NVARI, NSTRS 15 : * NMATT : Nombre de parametres materiaux 16 : * NVARI : Nombre de variables internes 17 : * NSTRS : Nombre de composantes 18 : 19 : REAL*8 XMAT(NMATT), VAR0(NVARI), SIG0(NSTRS) 20 : REAL*8 DEPST(NSTRS) 21 : * XMAT : Paramètres matériaux 22 : * VAR0 : Variables internes 23 : * SIG0 : Contraintes reelles 24 : * DEPST : Increment de deformations totales 25 : 26 : * Sorties 27 : 28 : REAL*8 SIGF(NSTRS), VARF(NVARI) 29 : * SIGF : Contraintes effective 30 : * VARF : Variables internes en sortie de plasticite 31 : 32 : * Variables materiau 33 : REAL*8 E,NU 34 : REAL*8 AC, AT, BC, BT, EPSD0, FC, FT 35 : REAL*8 P, AH, BH, CH, GAMMA, ALFA 36 : REAL*8 A, K0 37 : REAL*8 DEUXMU, LAMB, TREPS 38 : 39 : * Données intermédiaires et autres 40 : INTEGER ITER, I , J, K, KGLOB, KTEST 41 : LOGICAL CONV 42 : REAL*8 TEST,F0 43 : REAL*8 LIM, DEPST2(6) 44 : REAL*8 SIGE0(6), SIGE0P(3), TOL, TOLDYN 45 : REAL*8 VECPC(3,3), NORM, DSIGT(6) 46 : REAL*8 DFDS(6), H, SIG02(6), VAR02(9) 47 : REAL*8 FTB, FH, FM(6), C 48 : REAL*8 TRS, KH, KH0, KH2, TR(3,3) 49 : REAL*8 DMP, RES(8), RESUM(6), RESU(6), KHM, ELAS(6,6) 50 : REAL*8 JAC(7,7) 51 : REAL*8 DFDK, N(6), DHDK, KH1, DHDS(6), DMDK(6) 52 : REAL*8 FHP(6), FMM1(6) 53 : REAL*8 FHM(6) 54 : REAL*8 FMP1(6) 55 : REAL*8 DMDS(6,6), FHPH, FHMH, ETA, ETAP 56 : REAL*8 FHM1 57 : REAL*8 FHP1, SIGP(6), SIGMM(6), DELTA(8) 58 : REAL*8 SIGMP1(6),DSDE(7,7) 59 : REAL*8 M 60 : LOGICAL PLASTI,CONV1,CONV2,TEST5 61 : INTEGER FLAG 62 : 63 : REAL*8 ELASM(6,6) 64 : REAL*8 ES(7), NUMF, NUM(7), DENOF 65 : REAL*8 DENO(7), MAXI, MAXI2, SIGE00(6) 66 : 67 : 68 : 69 : *************************************** 70 : * Récupération des variables matériau * 71 : *************************************** 72 : 73 : E = XMAT(1) 74 : NU = XMAT(2) 75 : AC = XMAT(5) 76 : AT = XMAT(6) 77 : BC = XMAT(7) 78 : BT = XMAT(8) 79 : EPSD0 = XMAT(9) 80 : FC = XMAT(10) 81 : FT = XMAT(11) 82 : P = XMAT(12) 83 : AH = XMAT(13) 84 : BH = XMAT(14) 85 : CH = XMAT(15) 86 : GAMMA = XMAT(16) 87 : ALFA = XMAT(17) 88 : A = XMAT(18) 89 : K0 = XMAT(19) 90 : 91 : DEUXMU=E/(1.D0+NU) 92 : LAMB=NU*DEUXMU/(1.D0-2.D0*NU) 93 : 94 : 95 : 96 : ******************************* 97 : * Initialisationdes variables * 98 : ******************************* 99 : 100 : * Calcul de la matrice élastique (une seule fois !) 101 : 102 : DO I = 1,6 103 : DO J = 1,6 104 : ELAS(I,J) = 0.d0; 105 : ELASM(I,J) = 0.d0; 106 : END DO 107 : END DO 108 : 109 : ELAS(1,1)=LAMB+DEUXMU 110 : ELAS(2,2)=LAMB+DEUXMU 111 : ELAS(3,3)=LAMB+DEUXMU 112 : ELAS(1,2)=LAMB 113 : ELAS(2,1)=LAMB 114 : ELAS(1,3)=LAMB 115 : ELAS(3,1)=LAMB 116 : ELAS(2,3)=LAMB 117 : ELAS(3,2)=LAMB 118 : ELAS(4,4)=DEUXMU 119 : ELAS(5,5)=DEUXMU 120 : ELAS(6,6)=DEUXMU 121 : 122 : ELASM(1,1)=1.D0/E 123 : ELASM(2,2)=1.D0/E 124 : ELASM(3,3)=1.D0/E 125 : ELASM(1,2)=(-NU)/E 126 : ELASM(2,1)=(-NU)/E 127 : ELASM(1,3)=(-NU)/E 128 : ELASM(3,1)=-NU/E 129 : ELASM(2,3)=-NU/E 130 : ELASM(3,2)=-NU/E 131 : ELASM(4,4)=(1.D0+NU)/E 132 : ELASM(5,5)=(1.D0+NU)/E 133 : ELASM(6,6)=(1.D0+NU)/E 134 : 135 : 136 : 137 : ****************** 138 : * Initialisation * 139 : ****************** 140 : 141 : * Variables tampons 142 : 143 : DO I=1,NSTRS 144 : DEPST2(I) = DEPST(I) 145 : SIG02(I) = SIG0(I) 146 : END DO 147 : 148 : DO I=1,NVARI 149 : VAR02(I) = VAR0(I) 150 : END DO 151 : 152 : KGLOB = 1 153 : KTEST = 1 154 : 155 : * Test d'iterations globales 156 : TEST5 = .TRUE. 157 : 158 : *********************** 159 : * Iterations globales * 160 : *********************** 161 : 162 : DO WHILE (TEST5) 163 : 164 : * On saisit une petite valeur de KH et on saisit le plus grand à chaque itération... 165 : KH0 = 1.d-3 166 : KH = MAX(VAR0(9), KH0) 167 : 168 : * Initialisation de la variable KH à l'instant t 169 : KHM = KH 170 : 171 : * Initialisation de l'incrément du multiplicateur plastique 172 : DMP = 0.d0 173 : 174 : 175 : * Variables de test 176 : TEST = 1.d-6 177 : FLAG = 0 178 : 179 : 180 : 181 : **************************************************** 182 : * Calcul de l'incrément de contraintes (élastique) * 183 : **************************************************** 184 : 185 : * Prediction elastique 186 : TREPS = DEPST(1) + DEPST(2) + DEPST(3) 187 : DO K=1,3 188 : DSIGT(K) = LAMB*TREPS + DEUXMU*DEPST(K) 189 : END DO 190 : DO K=4,NSTRS 191 : DSIGT(K) = DEUXMU*DEPST(K) 192 : END DO 193 : 194 : * Initialisation des variables plastiques 195 : CONV=.FALSE. 196 : 197 : DO I=1, NSTRS 198 : DFDS(I) = 0.d0 199 : END DO 200 : 201 : DO K = 1, NSTRS 202 : SIGE0(K) = VAR0(K+2) 203 : END DO 204 : 205 : ITER=0 206 : 207 : * Prediction élastique (NSTRS : nombre de composantes de SIGMA) 208 : DO ISTRS = 1,NSTRS 209 : 210 : * La prédiction élastique est donnée par la contrainte réelle SIG0 211 : * à laquelle on somme la différence entre contrainte plastique et 212 : * contrainte réelle stockée en variable interne et la prédiction 213 : * élastique 214 : 215 : * Initialisation pour une itération lorsque l'incrément est divisé : 216 : * - Si c'est la première itération, on reprend la valeur initiale de 217 : * la contrainte plastique pour créer la prédiction 218 : * - Si c'est l'itération suivante, on travaille directement sur 219 : * la contrainte plastique obtenue à l'itération précédente et 220 : * stockée dans SIG0. On ne recalcule pas la contrainte plastique. 221 : 222 : IF (KTEST.EQ.1) THEN 223 : 224 : * SigE0 représente la différence entre la valeur de la contrainte dans le cas réel 225 : * et la valeur de la contrainte dans le cas plastique à l'instant t 226 : SIGMP1(ISTRS) = SIG0(ISTRS) + SIGE0(ISTRS) 227 : SIGE00(ISTRS) = SIG0(ISTRS) + SIGE0(ISTRS) 228 : SIGE0(ISTRS) = SIG0(ISTRS) + SIGE0(ISTRS) + DSIGT(ISTRS) 229 : 230 : ELSE 231 : SIGMP1(ISTRS) = SIG0(ISTRS) 232 : SIGE00(ISTRS) = SIG0(ISTRS) 233 : SIGE0(ISTRS) = SIG0(ISTRS) + DSIGT(ISTRS) 234 : 235 : END IF 236 : 237 : END DO 238 : 239 : * Apres, SIGE0 represente la solution (contrainte effective) 240 : 241 : ************************************* 242 : * Itération de calcul de plasticité * 243 : ************************************* 244 : 245 : DO WHILE (ITER .LE. 20 .AND. .NOT. CONV) 246 : ITER=ITER+1 247 : 248 : * Critere sur la prediction permettant de ramener la solution 249 : * a un domaine dans lequel les equations du modele sont definies 250 : * Note : Il a été modifié de telle sorte que l'on applique le critère jusqu'à 251 : * ce que la prédiction soit acceptable (10 itérations max modifiable) 252 : TRS = 0.D0 253 : DO I = 1,3 254 : TRS = TRS + SIGE0(I) 255 : END DO 256 : 257 : FTB = FT / FC 258 : M = 3.D0 * (1.D0-(FTB)**(2.D0/GAMMA))/(FTB + 259 : & 2.D0*(FTB)**(1.D0/GAMMA)) 260 : LIM = 3.D0*FC/M 261 : 262 : DO WHILE (TRS.GT.LIM .AND. FLAG.LT.10) 263 : 264 : FLAG = FLAG + 1 265 : 266 : IF (FLAG .LT. 10) THEN 267 : DO I =1,NSTRS 268 : SIGE0(I) = SIGE00(I) 269 : DSIGT(I) = DSIGT(I)/2.d0 270 : SIGE0(I) = SIGE00(I) + DSIGT(I) 271 : END DO 272 : DMP = 0.d0 273 : KH = max(var0(9),KH0) 274 : ITER = 2 275 : 276 : TRS = 0.D0 277 : 278 : DO I = 1,3 279 : TRS = TRS + SIGE0(I) 280 : END DO 281 : 282 : IF(TRS.GT.LIM) THEN 283 : write(6,*) 'ERREUR' 284 : END IF 285 : 286 : ELSE 287 : 288 : DO I =1,6 289 : SIGE0(I) = SIGE00(I) 290 : END DO 291 : 292 : END IF 293 : 294 : END DO 295 : 296 : * Calcul de FM et FH 296 : CALL FLOWM(FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH, C, 297 : & SIGE0, KH, F0, FM, ETA) 298 : 298 : CALL FLOWH(SIGE0, KH, AH, BH, CH, FC, FM, FH) 299 : 300 : * Calcul des résidus + vérification du critère de CV 301 : DO I=1,6 302 : RESU(I) = 0.D0 303 : RESUM(I) = 0.D0 304 : DO J=1,6 305 : RESU(I) = RESU(I) + ELASM(I,J)*SIGMP1(J) 306 : RESUM(I) = RESUM(I) + ELASM(I,J)*SIGE0(J) 307 : END DO 308 : END DO 309 : 310 : DO I = 1,6 311 : RES(I) = RESUM(I) + DMP*FM(I) - DEPST(I) - RESU(I) 312 : END DO 313 : 314 : RES(7) = KH - DMP*FH - KHM 315 : RES(8) = F0 316 : 317 : IF (KH .GE. 1.D0) THEN 318 : RES(7) = 0.D0 319 : KH = 1.d0 320 : END IF 321 : 322 : * Convergence 323 : * Deux criteres : un absolu, un relatif... Critere relatif conservé 324 : * Les variations de contraintes sont normées avec le max des contraintes 325 : * en fin d'itération et celles de KH et lambda sont normées par les valeurs 326 : * associées en sortie d'itération 327 : 328 : IF (((ABS(RES(8)).LT.TEST) .OR. (RES(8).LE.0.D0)) 329 : & .AND. ITER.EQ.1 .AND. (FLAG.EQ.0)) THEN 330 : CONV = .TRUE. 331 : ELSE 332 : CONV1 = .TRUE. 333 : MAXI = ABS(RES(1)) 334 : MAXI2 = ABS(SIGE0(1)) 335 : 336 : DO I=1,8 337 : MAXI = MAX(ABS(RES(I)),MAXI) 338 : END DO 339 : DO I=1,6 340 : MAXI2 = MAX(ABS(SIGE0(I)),MAXI2) 341 : END DO 342 : 343 : IF (MAXI .LT. TEST) THEN 344 : CONV1 = .TRUE. 345 : ELSE 346 : CONV1 = .FALSE. 347 : END IF 348 : 349 : CONV2 = .TRUE. 350 : DO I=1,NSTRS 351 : IF ((ABS(DELTA(I)) .LT. MAXI2*TEST) 352 : & .AND. (CONV2)) THEN 353 : CONV2 = .TRUE. 354 : ELSE 355 : CONV2 = .FALSE. 356 : END IF 357 : END DO 358 : 359 : IF ((ABS(DELTA(7)) .LT. ABS(KH)*TEST) .AND. (CONV2)) THEN 360 : CONV2 = .TRUE. 361 : ELSE 362 : CONV2 = .FALSE. 363 : END IF 364 : 365 : IF ((ABS(DELTA(8)) .LT. ABS(DMP)*TEST) .AND. (CONV2)) THEN 366 : CONV2 = .TRUE. 367 : ELSE 368 : CONV2 = .FALSE. 369 : END IF 370 : 371 : IF (CONV2) THEN 372 : CONV = .TRUE. 373 : END IF 374 : 375 : IF (ITER .EQ. 1) THEN 376 : CONV = .FALSE. 377 : END IF 378 : 379 : END IF 380 : 381 : IF (.NOT. CONV) THEN 382 : 383 : PLASTI = .TRUE. 384 : 385 : 386 : 387 : * Calcul de la dérivée de m : différentiation numérique 388 : * Note : le choix de h joue un role non négligeable 389 : * dans l'efficacité de l'algorithme 389 : CALL ENDOCB(SIGE0,TR,1,2) 389 : CALL JACOB3(TR,3,SIGE0P,VECPC) 390 : 391 : NORM= 0.D0 392 : DO I = 1,3 393 : NORM=NORM +(SIGE0P(I)**2) 394 : ENDDO 395 : NORM = DSQRT(NORM) 396 : 397 : H= NORM*2.d0**(-12) 398 : 399 : 400 : 401 : * Différenciation numérique 402 : DO J = 1,6 403 : 404 : DO I = 1,6 405 : SIGP(I) = SIGE0(I) 406 : SIGMM(I)= SIGE0(I) 407 : ENDDO 408 : 409 : SIGP(J) = SIGP(J) + H 410 : SIGMM(J) = SIGMM(J) - H 411 : 411 : CALL FLOWM (FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH, 412 : & C, SIGP, KH, F0, FMP1, ETAP) 412 : CALL FLOWM (FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH, 413 : & C, SIGMM, KH, F0, FMM1, ETAP) 414 : 414 : CALL FLOWH(SIGP, KH, AH, BH, CH, FC, FMP1, FHP1) 414 : CALL FLOWH(SIGMM, KH, AH, BH, CH, FC, FMM1, FHM1) 415 : 416 : 417 : 418 : * Et calcul de DMDS(i,j) = DM(i)/DS(J) 419 : 420 : DO I = 1,6 421 : DMDS(I,J) = (FMP1(I)-FMM1(I))/(2.D0*H) 422 : DHDS(J) = (FHP1-FHM1)/(2.D0*H) 423 : END DO 424 : 425 : END DO 426 : 427 : H=KH*2.d0**(-12) 428 : 429 : KH1= KH+H 430 : 431 : IF (KH1.GT.1.D0) THEN 432 : KH1 = KH 433 : END IF 434 : 434 : CALL FLOWM(FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH, C, 435 : & SIGE0, KH1, F0, FHP, ETAP) 435 : CALL FLOWH(SIGE0, KH1, AH, BH, CH, FC, FHP, FHPH) 436 : 437 : 438 : KH1 = KH -H 439 : 440 : * Traitement du cas kh = 1 : d2FdK2 = 0 441 : 442 : IF (KH.GE.1.D0) THEN 443 : KH1=KH 444 : END IF 445 : 445 : CALL FLOWM(FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH, C, 446 : & SIGE0, KH1, F0, FHM, ETAP) 447 : 447 : CALL FLOWH(SIGE0, KH1, AH, BH, CH, FC, FHM, FHMH) 448 : 449 : 450 : DO I=1,6 451 : DMDK(I) = (FHP(I) - FHM(I))/(2.D0*H) 452 : END DO 453 : 454 : DHDK = (FHPH-FHMH)/(2.D0*H) 455 : 456 : 457 : 458 : * Calcul de N = DFDS 459 : 460 : DO I = 1,6 461 : N(I) = FM(I) 462 : END DO 463 : 464 : DFDK = ETA 465 : 466 : DO I=1,6 467 : DO J=1,6 468 : JAC(I,J)= ELASM(I,J) + DMP*DMDS(I,J) 469 : ENDDO 470 : JAC(I,7) = DMP*DMDK(I) 471 : END DO 472 : 473 : DO J=1,6 474 : JAC(7,J) = -DMP*DHDS(J) 475 : END DO 476 : 477 : JAC(7,7) = 1.D0-DMP*DHDK 478 : 478 : CALL INVALM(JAC,7,7,KERRE,1.d-20) 479 : 480 : IF (KERRE.EQ.1) THEN 481 : WRITE(6,*) 'Erreur d inversion de matrice' 482 : END IF 483 : 484 : 485 : 486 : **************************************************************** 487 : * Calcul des variations en sortie d'itération (Newton Raphson) * 488 : **************************************************************** 489 : 490 : DO I=1,7 491 : DENO(I) = 0.d0 492 : NUM(I) = 0.d0 493 : END DO 494 : 495 : DENOF=0 496 : NUMF=0 497 : 498 : DO I = 1,8 499 : DELTA(I) = 0.d0 500 : END DO 501 : 502 : DO I=1,7 503 : DO J=1,6 504 : DENO(I) = DENO(I)+JAC(I,J)*FM(J) 505 : END DO 506 : DENO(I) = DENO(I)-JAC(I,7)*FH 507 : END DO 508 : 509 : DO I=1,6 510 : DENOF = DENOF+N(I)*DENO(I) 511 : END DO 512 : 513 : DENOF = DENOF+DFDK*DENO(7) 514 : 515 : DO I =1,7 516 : DO J = 1,7 517 : NUM(I)= NUM(I)+JAC(I,J)*RES(J) 518 : END DO 519 : END DO 520 : 521 : DO I =1,6 522 : NUMF = NUMF + N(I)*NUM(I) 523 : END DO 524 : 525 : NUMF = NUMF + DFDK*NUM(7) 526 : NUMF = RES(8)-NUMF 527 : 528 : DELTA(8) = NUMF/DENOF 529 : 530 : DO I =1,6 531 : ES(I) = -RES(I)-DELTA(8)*FM(I) 532 : END DO 533 : 534 : ES(7) = -RES(7)+DELTA(8)*FH 535 : 536 : DO I = 1,7 537 : DO J = 1,7 538 : DELTA(I)=DELTA(I)+JAC(I,J)*ES(J) 539 : END DO 540 : END DO 541 : 542 : 543 : * B est un vecteur à 8 composantes composées des DELTA (S, kh et LAMBDA) 544 : 545 : DO I = 1, 6 546 : SIGE0(I) = SIGE0(I) + DELTA(I) 547 : END DO 548 : 549 : KH2 = KH 550 : KH = KH + DELTA(7) 551 : 552 : 553 : 554 : * Test de KH < 1 555 : 556 : IF ((KH2.GE. 1.D0).OR.(KH.GE.1.D0)) THEN 557 : KH = 1.D0 558 : END IF 559 : 560 : 561 : 562 : * KH ne diminue pas... 563 : 564 : IF (KH .LT. KHM) THEN 565 : KH = KHM 566 : END IF 567 : 568 : DMP = DMP + DELTA(8) 569 : 570 : END IF 571 : 572 : END DO 573 : 574 : 575 : 576 : ******************* 577 : * Non convergence * 578 : ******************* 579 : 580 : IF ((ITER .EQ. 21) .AND. (.NOT. CONV) 581 : & .AND. (KGLOB.EQ.32)) THEN 582 : WRITE(6,*) '**************************************' 583 : WRITE(6,*) '* Divergence de l algorithme interne *' 584 : WRITE(6,*) '**************************************' 585 : WRITE(6,*) 'DEPST :',DEPST2 586 : WRITE(6,*) 'DELTA :',DELTA 587 : WRITE(6,*) 'MAXSIG,KH,DMP :',MAXI2,KH,DMP 588 : WRITE(6,*) 'VAR0 :',VAR02 589 : WRITE(6,*) 'SIG0 :',SIG02 590 : WRITE(6,*) 'SIGF :',SIGF 591 : * KERRE permet de stopper CAST3M et de reprendre la main 592 : * lorsqu'il est non nul 593 : KERRE = 1 594 : END IF 595 : 596 : 597 : 598 : ***************************************************** 599 : * Traitement des données en fin d'itération globale * 600 : ***************************************************** 601 : 602 : DO I=1,NSTRS 603 : SIG0(I) = SIGE0(I) 604 : END DO 605 : 606 : VAR0(9) = KH 607 : 608 : * Test de convergence 609 : IF (KGLOB.EQ.KTEST .AND. CONV) THEN 610 : TEST5 = .FALSE. 611 : ELSE IF (CONV) THEN 612 : KTEST = KTEST + 1 613 : END IF 614 : 615 : IF (.NOT.CONV) THEN 616 : KGLOB = KGLOB*2 617 : KTEST = 1 618 : DO I = 1,NSTRS 619 : SIG0(I) = SIG02(I) 620 : DEPST(I) = DEPST(I)/2.d0 621 : END DO 622 : VAR0(9) = VAR02(9) 623 : END IF 624 : 625 : * Limite de la division de l'intervalle 626 : * On redivise par 2 pour l'afficher en sortie 627 : IF (KGLOB.GT.32) THEN 628 : TEST5=.FALSE. 629 : KGLOB = KGLOB/2 630 : END IF 631 : 632 : END DO 633 : 634 : 635 : 636 : ******************************************* 637 : * FIN DE LA BOUCLE GLOBALE * 638 : * --> Sorties de la routine de plasticite * 639 : ******************************************* 640 : 641 : * Les résultats finaux sont stockés dans SIG0 et VAR0 642 : * Les valeurs en entrée sont stockés dans SIG02 et VAR02 643 : DO I = 1,NSTRS 644 : SIGF(I) = SIG0(I) 645 : SIG0(I) = SIG02(I) 646 : DEPST(I) = DEPST2(I) 647 : END DO 648 : 649 : VARF(9) = VAR0(9) 650 : VAR0(9) = VAR02(9) 651 : 652 : * On affiche s'il y a eu division de l'incrément de déformations 653 : IF (KGLOB.GT.1) THEN 654 : * WRITE(*,*) 'Division de l increment de deformations par',KGLOB 655 : * WRITE(*,*) 'D :',VAR02(2),' et KH :',KH 656 : END IF 657 : 658 : RETURN 659 : 660 : END 661 : 662 : 663 :
© Cast3M 2003 - All rights reserved.
Disclaimer