1 : C BETOCY SOURCE PV 11/03/07 21:15:06 6885 2 : SUBROUTINE BETOCY(WRK0,WRK1,WRK2,NCURVT,NCURVC,KERRE) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : C 6 : C MODELE DE MACONNERIE ENDOMMAGEABLE EN CYCLIQUE 7 : C 8 : C======================================================================= 9 : C CETTE ROUTINE EST APPELE DANS ECOUL2 10 : C 11 : C ENTREES: 12 : C ------- 13 : C SIG0 = CONTRAINTES AU DEBUT DU PAS D'INTEGRATION 14 : C DSIGT = INCREMENT DE CONT. CALCULE ELASTIQUEMENT A PARTIR DE 15 : C L'INCREMENT DES DEFORMATIONS TOTALES 16 : C VAR0 = VARIABLES INTERNES AU DEBUT DU PAS D'INTEGRATION 17 : C 18 : C XMAT = CARACTERISTIQUES MECANIQUES DU MATERIAU 19 : C 20 : C SORTIES: 21 : C ------- 22 : C SIGF = CONTR. A LA FIN DU PAS D'INTEGRATION 23 : C VARF = VARIABLES INTERNES A LA FIN DU PAS D'INTEGRATION 24 : C DEFP = INCREMENT DES DEFORM. PLASTIQUES A LA FIN DU PAS 25 : C D'INTEGRATION 26 : C KERRE = INDICE QUI REGIT LES ERREURS 27 : C 28 : C========================================================================= 29 : C VARIABLES NECESSAIRES 30 : C========================================================================= 31 : -INC CCOPTIO 32 : C 33 : C Segment pour le Broyden 34 : C 35 : SEGMENT QUASIN 36 : REAL*8 XH0(MN,MN) 37 : REAL*8 WUP(MN,ITM),YUP(MN,ITM) 38 : INTEGER IT 39 : REAL*8 R(MN) 40 : REAL*8 D(MN),DES(MN) 41 : ENDSEGMENT 42 : C Segment des materiaux 43 : SEGMENT WRK0 44 : REAL*8 XMAT(NCXMAT) 45 : ENDSEGMENT 46 : C Segment des courbes d'ecrouissage 47 : SEGMENT WRK2 48 : REAL*8 TRAC(LTRAC) 49 : ENDSEGMENT 50 : C Segment des entrees/sortie de ecoul2 51 : C Les contraintes sont sous la forme [SXX,SYY,0,SXY] 52 : C Vecteurs de contrainte totale 53 : C 4 composantes pour les contraintes SIGT 54 : SEGMENT WRK1 55 : REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS) 56 : REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI) 57 : REAL*8 DEFP(NSTRS),XCAR(ICARA) 58 : ENDSEGMENT 59 : DIMENSION DSIGT(4) 60 : DIMENSION SIGT(4), SIGI(4), DSIGI(4), DSIGIDI(4) 61 : C Matrice de compliance (C-1) 62 : DIMENSION CINV(3) 63 : C Variables d'ecrouissage isotrope 64 : DIMENSION XKISO(2) 65 : C Travaux plastiques cycliques et Deformations plastiques equivalentes 66 : C associes au mecanisme de compression de la grande surf. 67 : REAL*8 WORK 68 : DIMENSION EPSI(2),EPSI0(2),EPSII(2) 69 : C Contraintes principales et variables d'ecrouissage cinematique 70 : DIMENSION SIGTP(2), SIGHP(2),SIGHP0(2) 71 : C Cos et Sin de l'angle definissant la premiere direction principale 72 : REAL*8 CPHI, CPHI0, CPHI00, SPHI, SPHI0, SPHI00 73 : C Indicateur du contact 74 : DIMENSION KONTAC(4),KONTA0(4) 75 : C Vecteurs contrainte transitoire 76 : DIMENSION SITPRO(4) 77 : C Jacobien initial 78 : DIMENSION XJ00(3,3) 79 : C Variable phi apparaissant au denominateur (calcul de la normale et de A) 80 : REAL*8 PHI 81 : C Matrices contenant A-1 82 : C AINV --> Inverse de A pour la grande surface 83 : C [ A B 0] A=AINV(1) 84 : C AINVGRANDESURFACE = [-B A 0] B=AINV(2) 85 : C [ 0 0 C] C=AINV(3) 86 : DIMENSION AINV(3) 87 : C Vecteurs normaux aux surface 88 : C XNOR1 --> 1eme mecanisme 89 : C XNOR2 --> 2eme mecanisme 90 : DIMENSION XNOR1(4), XNOR2(4) 91 : C 92 : C Multiplicateurs plastiques 2 mecanismes + 2 surfaces dans GAMBDA 93 : C Valeurs des criteres dans FFF 94 : C Valeurs des limites dans RESIST 95 : C Ordre pour les 3 tableaux: 96 : C 1 --> Traction petite surface direction 1 97 : C 2 --> Traction petite surface direction 2 98 : C 3 --> Compression petite surface direction 1 99 : C 4 --> Compression petite surface direction 2 100 : C 5 --> Traction grande surface 101 : C 6 --> Compression grande surface 102 : DIMENSION GAMBDA(6), FFF(6), RESIST(6) 103 : C Valeur necessaire pour le calcul des residus de la grande surface 104 : DIMENSION XKISO0(2),XKISOI(2) 105 : C 106 : C Separation de l'increment en deux ou trois parties dans le cas du contact 107 : C 108 : REAL*8 SPLITT 109 : REAL*8 XFFFS 110 : C 111 : C Pour differencier les differents cas 112 : C ICAS --> NUMERO DU MULTIPLICATEUR PLASTIQUE A METTRE A JOUR 113 : C IHARD --> NUMERO DE LA VARIABLE D'ECROUISSAGE CINEMATIQUE 114 : C ISGNTC --> Signe indiquant si on est en traction ou compression 115 : C ISGN2 --> Signe indiquant si on prend la plus grande ou la plus 116 : C petite contrainte principale 117 : C NCORNER --> Revient-on sur la pointe ? 118 : INTEGER ICAS, IHARD, ISGNTC, ISGN2 119 : INTEGER ICASB, IHARDB, ISGNTCB, ISGN2B 120 : INTEGER NCORNER 121 : C 122 : C Variable intermediaire valant la resistance de la grande surface 123 : C 124 : DIMENSION FG(2) 125 : C 126 : C Critere de convergence 127 : C 128 : DIMENSION EPSIL(6) 129 : CHARACTER*8 CMATE 130 : C======================================================================== 131 : C Parametres utiles 132 : C======================================================================== 133 : PARAMETER (XZER=0.D0,UN=1.D0,DEUX=2.D0,UNDEMI=.5D0) 134 : PARAMETER (UNRADE=0.707106781D0) 135 : PARAMETER (EPSILO=1.D-8,EPSILO2=-1.D-10) 136 : C======================================================================== 137 : C 138 : C CALCUL DE DSIGT 139 : C 140 : C calcul de la matrice elastique 141 : C 142 : CMATE = 'ISOTROPE' 143 : KCAS=1 143 : CALL DOHMAS(XMAT,CMATE,IFOUR,4,KCAS,DDHOOK,IRTD) 144 : DO IE1=1,4 145 : DSIGT(IE1)=XZER 146 : DO IE2=1,4 147 : DSIGT(IE1)=DSIGT(IE1)+DDHOOK(IE1,IE2)*DEPST(IE2) 148 : ENDDO 149 : ENDDO 150 : DSIGT(3)=XZER 151 : C========================================================================= 152 : C LECTURE DES PARAMETRES DU MATERIAU 153 : C 154 : C========================================================================= 155 : C Module d'Young YOUN, Poisson XNU, Shear GG 156 : C Ecrouissage cinematique XHH1 , 157 : YOUN = XMAT(1) 158 : XNU = XMAT(2) 159 : GGG = YOUN/(2*(UN+XNU)) 160 : XHH1 = XMAT(5) 161 : CPI = YOUN/(UN-XNU) 162 : C 163 : RESIST(1) = XMAT(6) 164 : RESIST(2) = XMAT(6) 165 : RESIST(3) = XMAT(7) 166 : RESIST(4) = XMAT(7) 167 : RESIST(5) = XMAT(8) 168 : RESIST(6) = XMAT(9) 169 : C Deformation plastique equivalente et travail plastique de ref. 170 : C pour la dissipation cyclique 171 : EPS0 = XMAT(10) 172 : WOR0 = XMAT(11) 173 : DEGRAD = EPS0/WOR0 174 : XLCAT = XMAT(14) 175 : XLCAC = XMAT(15) 176 : C 177 : C 178 : C Evolution des variables d'ecrouissage isotrope 179 : C 180 : C Voila les positions des points d'entree des deux courbes dans TRAC 181 : IPT = 1 182 : IPC = IPT+2*NCURVT 183 : C NCURVT : nombre de point de la courbe de traction 184 : C NCURVC : nombre de point de la courbe de compression 185 : C=========================================================== 186 : C Fin lecture des caracteristiques 187 : C On passe aux contraintes et variables d'etat 188 : C=========================================================== 189 : C 190 : C CONTRAINTES ELASTIQUES 191 : C 192 : SIGI(1) = SIG0(1) 193 : SIGI(2) = SIG0(2) 194 : SIGI(4) = SIG0(4) 195 : C DELTA SIGMA TRIAL 196 : DSIGI(1) = DSIGT(1) 197 : DSIGI(2) = DSIGT(2) 198 : DSIGI(4) = DSIGT(4) 199 : C 200 : C ECROUISSAGE CINEMATIQUE 201 : C 202 : SIGHP0(1) = VAR0(1) 203 : SIGHP0(2) = VAR0(2) 204 : C 205 : C ECROUISSAGE ISOTROPE 206 : C 207 : EPSI0(1) = VAR0(3) 208 : EPSI0(2) = VAR0(4) 209 : C 210 : CPHI00 = VAR0(5) + 1.D0 211 : CPHI0 = CPHI00 212 : SPHI00 = VAR0(6) 213 : SPHI0 = SPHI00 214 : C 215 : KONTA0(1) = nint(VAR0( 7)) 216 : KONTA0(2) = nint(VAR0( 8)) 217 : KONTA0(3) = nint(VAR0( 9)) 218 : KONTA0(4) = nint(VAR0(10)) 219 : C 220 : C Variables de contact 221 : C 222 : SPLITT = 0.D0 223 : NCON = 0 224 : IDICHO = 0 225 : NDICHO = 0 226 : NPASS = 0 227 : DSIGIDI(1) = XZER 228 : DSIGIDI(2) = XZER 229 : DSIGIDI(4) = XZER 230 : C=============================================== 231 : C Fin de la lecture des contraintes et des variables d'etat 232 : C================================================ 233 : C 234 : C Premier calcul des residus puis DRIVER 235 : C================================================ 236 : C 237 : C Entree quand on revient des routines, dernier test avant de finir 238 : C 239 : 10 CONTINUE 240 : SIGT(1) = SIGI(1) + DSIGI(1) 241 : SIGT(2) = SIGI(2) + DSIGI(2) 242 : SIGT(4) = SIGI(4) + DSIGI(4) 242 : CALL CONPRI (SIGT,CPHI0,SPHI0,SIGTP,CPHI,SPHI) 243 : EPSI(1) = EPSI0(1) 244 : EPSI(2) = EPSI0(2) 244 : CALL VALEUR(EPSI(1),TRAC(IPT),NCURVT,XLCAT, 245 : &XKISO(1),DUMM,KERRE) 245 : CALL VALEUR(EPSI(2),TRAC(IPC),NCURVC,XLCAC, 246 : &XKISO(2),DUMM,KERRE) 247 : SIGHP(1) = SIGHP0(1) 248 : SIGHP(2) = SIGHP0(2) 249 : NPASS2 = 0 250 : KONTAC(1) = KONTA0(1) 251 : KONTAC(2) = KONTA0(2) 252 : KONTAC(3) = KONTA0(3) 253 : KONTAC(4) = KONTA0(4) 254 : C 255 : C Calcul des residus initiaux 256 : C 257 : NMEC = 0 258 : FFF(1) = XZER 259 : FFF(2) = XZER 260 : FFF(3) = XZER 261 : FFF(4) = XZER 262 : FFF(5) = XZER 263 : FFF(6) = XZER 264 : FINT1 = XZER 265 : FINT2 = XZER 266 : GAMBDA(1) = XZER 267 : GAMBDA(2) = XZER 268 : GAMBDA(3) = XZER 269 : GAMBDA(4) = XZER 270 : GAMBDA(5) = XZER 271 : GAMBDA(6) = XZER 272 : ICAS = 0 273 : ICASB = 0 274 : C 275 : C Critere de convergence 276 : C 277 : EPSIL(1) = ABS(RESIST(1)*EPSILO) 278 : EPSIL(2) = ABS(RESIST(2)*EPSILO) 279 : EPSIL(3) = ABS(RESIST(3)*EPSILO) 280 : EPSIL(4) = ABS(RESIST(4)*EPSILO) 281 : EPSIL(5) = ABS(RESIST(5)*EPSILO) 282 : EPSIL(6) = ABS(RESIST(6)*EPSILO) 283 : C 284 : C 285 : IF (KONTAC(1).EQ.0) THEN 285 : CALL CRITPE (SIGTP(1), SIGHP0(1), RESIST(1), 1, FFF(1)) 286 : FINT1 = 0. 287 : IF (FFF(1).GE.EPSIL(1)) THEN 288 : NMEC = NMEC + 1 289 : IF (NMEC.EQ.1) THEN 290 : ICAS = 1 291 : ISGNTC = 1 292 : IHARD = 1 293 : ELSE 294 : ICASB = 1 295 : ISGNTCB = 1 296 : IHARDB = 1 297 : ENDIF 298 : ENDIF 299 : ELSE 300 : FG(1) = XKISO(1)*RESIST(5) 300 : CALL CRITPE (SIGTP(1), XZER, FG(1), 1, FINT1) 301 : ENDIF 302 : C 303 : IF (KONTAC(2).EQ.0) THEN 303 : CALL CRITPE (SIGTP(2), SIGHP0(2), RESIST(2), 1, FFF(2)) 304 : FINT2 = 0. 305 : IF (FFF(2).GE.EPSIL(2)) THEN 306 : NMEC = NMEC + 1 307 : IF (NMEC.EQ.1) THEN 308 : ICAS = 2 309 : ISGNTC = 1 310 : IHARD = 2 311 : ELSE 312 : ICASB = 2 313 : ISGNTCB = 1 314 : IHARDB = 2 315 : ENDIF 316 : ENDIF 317 : ELSE 318 : FG(1) = XKISO(1)*RESIST(5) 318 : CALL CRITPE (SIGTP(2), XZER, FG(1), 1, FINT2) 319 : ENDIF 320 : C 321 : IF ((FINT1.GE.EPSIL(5)).OR.(FINT2.GE.EPSIL(5))) THEN 322 : NMEC = NMEC + 1 323 : IF (FINT1.GE.FINT2) THEN 324 : FFF(5) = FINT1 325 : IF (NMEC.EQ.1) THEN 326 : IHARD = 1 327 : ICAS = 5 328 : ISGNTC = 1 329 : ELSE 330 : IHARDB = 1 331 : ICASB = 5 332 : ISGNTCB = 1 333 : ENDIF 334 : ELSE 335 : FFF(5) = FINT2 336 : IF (NMEC.EQ.1) THEN 337 : IHARD = 2 338 : ICAS = 5 339 : ISGNTC = 1 340 : ELSE 341 : IHARDB = 2 342 : ICASB = 5 343 : ISGNTCB = 1 344 : ENDIF 345 : ENDIF 346 : ENDIF 347 : C 348 : IF (KONTAC(3).EQ.0) THEN 348 : CALL CRITPE (SIGTP(1), SIGHP0(1), RESIST(3), -1, FFF(3)) 349 : FINT1 = 0. 350 : IF (FFF(3).GE.EPSIL(3)) THEN 351 : NMEC = NMEC + 1 352 : IF (NMEC.EQ.1) THEN 353 : ICAS = 3 354 : ISGNTC = -1 355 : IHARD = 1 356 : ELSE 357 : ICASB = 3 358 : ISGNTCB = -1 359 : IHARDB = 1 360 : ENDIF 361 : ENDIF 362 : ELSE 363 : FG(1) = XKISO(2)*RESIST(6) 363 : CALL CRITPE (SIGTP(1), XZER, FG(1), -1, FINT1) 364 : ENDIF 365 : C 366 : IF (KONTAC(4).EQ.0) THEN 366 : CALL CRITPE (SIGTP(2), SIGHP0(2), RESIST(4), -1, FFF(4)) 367 : FINT2 = 0. 368 : IF (FFF(4).GE.EPSIL(4)) THEN 369 : NMEC = NMEC + 1 370 : IF (NMEC.EQ.1) THEN 371 : ICAS = 4 372 : ISGNTC = -1 373 : IHARD = 2 374 : ELSE 375 : ICASB = 4 376 : ISGNTCB = -1 377 : IHARDB = 2 378 : ENDIF 379 : ENDIF 380 : ELSE 381 : FG(1) = XKISO(2)*RESIST(6) 381 : CALL CRITPE (SIGTP(2), XZER, FG(1), -1, FINT2) 382 : ENDIF 383 : C 384 : IF ((FINT1.GE.EPSIL(6)).OR.(FINT2.GE.EPSIL(6))) THEN 385 : NMEC = NMEC + 1 386 : IF (FINT1.GE.FINT2) THEN 387 : FFF(6) = FINT1 388 : IF (NMEC.EQ.1) THEN 389 : IHARD = 1 390 : ICAS = 6 391 : ISGNTC = -1 392 : ELSE 393 : IHARDB = 1 394 : ICASB = 6 395 : ISGNTCB = -1 396 : ENDIF 397 : ELSE 398 : FFF(6) = FINT2 399 : IF (NMEC.EQ.1) THEN 400 : IHARD = 2 401 : ICAS = 6 402 : ISGNTC = -1 403 : ELSE 404 : IHARDB = 2 405 : ICASB = 6 406 : ISGNTCB = -1 407 : ENDIF 408 : ENDIF 409 : ENDIF 410 : C 411 : C 412 : C Le tableau FFF contient les residus 413 : C 414 : C================================================ 415 : C DRIVER 416 : C================================================ 417 : IF (NMEC.EQ.0) GOTO 9999 418 : C 419 : IF (NMEC.EQ.1) GOTO 1000 420 : C 421 : IF (NMEC.EQ.2) GOTO 4000 422 : C 423 : GOTO 9999 424 : C================================================== 425 : C CAS 1: 426 : C Un seul mecanisme est active 427 : C================================================== 428 : 1000 CONTINUE 429 : C 429 : CALL CONPRI (SIGT,CPHI0,SPHI0,SIGTP,CPHI,SPHI) 430 : C 431 : C Mis a jour de isgn2 432 : C 433 : IF (IHARD.EQ.1) THEN 434 : IF (SIGTP(1).GE.SIGTP(2)) THEN 435 : ISGN2 = 1 436 : ELSE 437 : ISGN2 = -1 438 : ENDIF 439 : HARDC = 2 440 : ELSE 441 : IF (IHARD.EQ.2) THEN 442 : IF (SIGTP(2).GE.SIGTP(1)) THEN 443 : ISGN2 = 1 444 : ELSE 445 : ISGN2 = -1 446 : ENDIF 447 : ENDIF 448 : HARDC = 1 449 : ENDIF 450 : C 451 : IF (IHARD.EQ.1) THEN 452 : IHARDC = 2 453 : ELSE 454 : IHARDC = 1 455 : ENDIF 456 : C 457 : NCORNER = 0 458 : C 459 : 1010 CONTINUE 460 : C 461 : C Definition du segment 462 : C 463 : ITM = 30 464 : MN = 1 465 : SEGINI,QUASIN 466 : IT = -1 467 : C 468 : C Calcul de xhh2 et xkiso 469 : C 470 : IF (ICAS.GE.5) THEN 471 : IF (ICAS.EQ.5) THEN 471 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT, 472 : &XKISO(1),XHH2,KERRE) 472 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC, 473 : &XKISO(2),DUMM,KERRE) 474 : ELSE 474 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT, 475 : &XKISO(1),DUMM,KERRE) 475 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC, 476 : &XKISO(2),XHH2,KERRE) 477 : ENDIF 478 : EPSI(1) = EPSI0(1) 479 : XKISO0(1) = XKISO(1) 480 : EPSI(2) = EPSI0(2) 481 : XKISO0(2) = XKISO(2) 482 : XHH2I = XHH2 483 : EPSII(1) = EPSI0(1) 484 : XKISOI(1) = XKISO(1) 485 : EPSII(2) = EPSI0(2) 486 : XKISOI(2) = XKISO(2) 487 : ENDIF 488 : C 489 : DO 1011,IE1=1,6 490 : FFF(IE1) = XZER 491 : 1011 CONTINUE 492 : C 493 : IF (ICAS.LE.4) THEN 493 : CALL CRITPE (SIGTP(IHARD),SIGHP0(IHARD),RESIST(ICAS), 494 : & ISGNTC,FFF(ICAS)) 495 : ELSE 496 : FG(1) = XKISOI(ICAS-4)*RESIST(ICAS) 496 : CALL CRITPE (SIGTP(IHARD),XZER,FG(1),ISGNTC,FINT1) 496 : CALL CRITPE (SIGTP(IHARDC),XZER,FG(1),ISGNTC,FINT2) 497 : FFF(ICAS) = MAX(FINT1,FINT2) 498 : ENDIF 499 : C 500 : C Calcul de la normale et du jacobien 501 : C 502 : IF (ICAS.LE.4) THEN 502 : CALL NORMBL(SIGT,RESIST(ICAS),ISGNTC,ISGN2,1,XNOR1) 503 : XJ00(1,1) = -XXCCYY(XNOR1, XNOR1, YOUN, XNU)-XHH1 504 : ELSE 504 : CALL NORMBL(SIGT,RESIST(ICAS),ISGNTC,ISGN2,2,XNOR1) 505 : XJ00(1,1) = -XXCCYY(XNOR1,XNOR1,YOUN,XNU) 506 : & -ISGNTC*RESIST(ICAS)*XHH2 507 : ENDIF 508 : C 509 : R(1) = + FFF(ICAS) 510 : GAMBDA(1)=XZER 511 : GAMBDA(2)=XZER 512 : GAMBDA(3)=XZER 513 : GAMBDA(4)=XZER 514 : GAMBDA(5)=XZER 515 : GAMBDA(6)=XZER 516 : C 517 : XJ0INV = UN/XJ00(1,1) 518 : C 518 : CALL ZERO(D,1,MN) 519 : XH0(1,1) = XJ0INV 520 : C 521 : C On rentre dans les iterations internes 522 : C 523 : DO I=0,ITM 524 : C 525 : C Appel de Broyden 526 : C 526 : CALL BROYDE (QUASIN) 527 : C 528 : GAMBDA(ICAS) = D(1) 529 : C 530 : C Calcul du phi et des termes de la matrice A et A-1 531 : C 532 : SITPRO(1) = SIGI(1) + DSIGI(1) 533 : & - UNDEMI*(GAMBDA(1) + GAMBDA(2) + GAMBDA(5))*CPI 534 : & + UNDEMI*(GAMBDA(3) + GAMBDA(4) + GAMBDA(6))*CPI 535 : SITPRO(2) = SIGI(2) + DSIGI(2) 536 : & - UNDEMI*(GAMBDA(1) + GAMBDA(2) + GAMBDA(5))*CPI 537 : & + UNDEMI*(GAMBDA(3) + GAMBDA(4) + GAMBDA(6))*CPI 538 : SITPRO(4) = SIGI(4) + DSIGI(4) 539 : C 540 : C Mis a jour des variables d'ecrouissage 541 : C 542 : SIGHP(1) = SIGHP0(1) + XHH1*(GAMBDA(1)-GAMBDA(3)) 543 : SIGHP(2) = SIGHP0(2) + XHH1*(GAMBDA(2)-GAMBDA(4)) 544 : EPSI(1) = EPSI0(1) + GAMBDA(5) 545 : IF (ICAS.LE.4) THEN 546 : WORK = GAMBDA(ICAS)*RESIST(ICAS)*ISGNTC 547 : ELSE 548 : WORK = XZER 549 : ENDIF 550 : EPSI(2) = EPSI0(2) + GAMBDA(6) + DEGRAD*WORK 551 : C 552 : IF (ICAS.LE.4) THEN 552 : CALL VALEUR(EPSI(1),TRAC(IPT),NCURVT,XLCAT, 553 : &XKISO(1),DUMM,KERRE) 553 : CALL VALEUR(EPSI(2),TRAC(IPC),NCURVC,XLCAC, 554 : &XKISO(2),DUMM,KERRE) 555 : ELSE 556 : IF (ICAS.EQ.5) THEN 556 : CALL VALEUR(EPSI(1),TRAC(IPT),NCURVT,XLCAT, 557 : &XKISO(1),XHH2,KERRE) 557 : CALL VALEUR(EPSI(2),TRAC(IPC),NCURVC,XLCAC, 558 : &XKISO(2),DUMM,KERRE) 559 : ELSE 559 : CALL VALEUR(EPSI(1),TRAC(IPT),NCURVT,XLCAT, 560 : &XKISO(1),DUMM,KERRE) 560 : CALL VALEUR(EPSI(2),TRAC(IPC),NCURVC,XLCAC, 561 : &XKISO(2),XHH2,KERRE) 562 : ENDIF 563 : ENDIF 564 : C 565 : C Calcul des phi 566 : C 567 : IF (ICAS.LE.4) THEN 568 : PHI = ISGNTC*(-UNDEMI*(SITPRO(1) + SITPRO(2)) + SIGHP(IHARD) 569 : & + RESIST(ICAS)) 570 : ELSE 571 : PHI = ISGNTC*(-UNDEMI*(SITPRO(1)+SITPRO(2)) 572 : & +(XKISOI(ICAS - 4)+(EPSI(ICAS-4) - EPSII(ICAS-4))*XHH2I) 573 : & *RESIST(ICAS)) 574 : ENDIF 575 : C 576 : C Calcul de AINV 577 : C 578 : SOMLAM = GAMBDA(1) + GAMBDA(2) + GAMBDA(3) + GAMBDA(4) 579 : & + GAMBDA(5) + GAMBDA(6) 580 : AINV(1) = (DEUX*PHI+GGG*SOMLAM)/(DEUX*PHI+DEUX*GGG*SOMLAM) 581 : AINV(2) = GGG*SOMLAM/(DEUX*PHI + DEUX*GGG*SOMLAM) 582 : AINV(3) = PHI/(PHI + (GGG*SOMLAM)) 583 : C 584 : C Calcul des contraintes a l'aide de AINV 585 : C 586 : SIGT(1) = AINV(1)*SITPRO(1) + AINV(2)*SITPRO(2) 587 : SIGT(2) = AINV(2)*SITPRO(1) + AINV(1)*SITPRO(2) 588 : SIGT(4) = AINV(3)*SITPRO(4) 589 : C 590 : C Nouvelles contraintes principales 591 : C 591 : CALL CONPRI (SIGT,CPHI0,SPHI0,SIGTP,CPHI,SPHI) 592 : C 593 : C 594 : IF (ICAS.LE.4) THEN 594 : CALL CRITPE (SIGTP(IHARD),SIGHP(IHARD),RESIST(ICAS), 595 : & ISGNTC,FFF(ICAS)) 596 : R(1)= + FFF(ICAS) 597 : ELSE 598 : FG(1) = (XKISOI(ICAS-4)+(EPSI(ICAS-4)-EPSII(ICAS-4))*XHH2I) 599 : & *RESIST(ICAS) 600 : C 600 : CALL CRITPE (SIGTP(IHARD),XZER,FG(1),ISGNTC,FINT1) 600 : CALL CRITPE (SIGTP(IHARDC),XZER,FG(1),ISGNTC,FINT2) 601 : FFF(ICAS) = MAX (FINT1,FINT2) 602 : C 603 : R(1) = + FFF(ICAS) 604 : & - ISGNTC*RESIST(ICAS) 605 : & * (XKISO(ICAS-4) - XKISOI(ICAS-4) 606 : & -(EPSI(ICAS-4) - EPSII(ICAS-4))*XHH2I) 607 : IF (ABS(EPSI(ICAS-4) - EPSII(ICAS-4)).LT.EPSILO2) THEN 608 : XHH2I=XHH2 609 : ELSE 610 : XHH2I = (XKISO(ICAS-4) - XKISOI(ICAS-4))/ 611 : & (EPSI(ICAS-4) - EPSII(ICAS-4)) 612 : ENDIF 613 : XKISOI(ICAS-4) = XKISO(ICAS-4) 614 : EPSII(ICAS-4) = EPSI(ICAS-4) 615 : ENDIF 616 : C 617 : C Test de convergence 618 : C 619 : C IF (ABS(FFF(ICAS)).LT.EPSIL(ICAS)) GOTO 1020 620 : IF ((ABS(FFF(ICAS)).LT.EPSIL(ICAS)) 621 : & .AND.(ABS(R(1)).LT.EPSIL(ICAS))) GOTO 1020 622 : ENDDO 623 : C 624 : C Pas de convergence 625 : C 626 : SEGSUP,QUASIN 627 : GOTO 6000 628 : C 629 : 1020 CONTINUE 630 : SEGSUP,QUASIN 631 : C 632 : IF (GAMBDA(ICAS).LT.EPSILO2) THEN 633 : WRITE (IOIMP,*)'BETOCY: Multiplicateurs plastiques < 0' 634 : KERRE = 2 635 : RETURN 636 : ELSE 637 : CONTINUE 638 : ENDIF 639 : C 640 : C On verifie qu'on ne viole pas les autres criteres 641 : C 642 : IF (IHARD.EQ.2) THEN 643 : IHARDC=1 644 : ELSE 645 : IHARDC=2 646 : ENDIF 647 : C 648 : IF (ICAS.LE.2) THEN 649 : FG(1) = XKISO(1)*RESIST(5) 650 : FG(2) = XKISO(2)*RESIST(6) 651 : ELSE 652 : IF (ICAS.LE.4) THEN 653 : FG(1) = XKISO(2)*RESIST(6) 654 : FG(2) = XKISO(1)*RESIST(5) 655 : ELSE 656 : IF (ICAS.EQ.5) THEN 657 : FG(1) = XKISO(1)*RESIST(5) 658 : FG(2) = XKISO(2)*RESIST(6) 659 : ELSE 660 : FG(1) = XKISO(2)*RESIST(6) 661 : FG(2) = XKISO(1)*RESIST(5) 662 : ENDIF 663 : ENDIF 664 : ENDIF 665 : C 666 : C Test GS 667 : C 668 : IF (ICAS.GE.5) THEN 669 : C Regles de suivi si on est sur la grande surface 670 : C 1- Direction de l'ecoulement 671 : SIGHP(IHARD) = FG(1) -RESIST(2*ICAS-10+IHARD) 672 : C 2- Direction perpendiculaire a l'ecoulement 673 : C (si on se trouve dans le coin) 674 : XLIMIT = FG(1) -RESIST(2*ICAS-10+IHARDC) 675 : IF (ISGNTC*(SIGHP(IHARDC)-XLIMIT).GT.XZER) THEN 676 : SIGHP(IHARDC) = XLIMIT 677 : KONTAC(2*ICAS-10+IHARDC)=1 678 : ENDIF 679 : ENDIF 680 : C 681 : C 682 : C On teste les autres mecanismes 683 : C 684 : C Test des contacts lateraux 685 : C Meme signe signtc 686 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.3)) THEN 687 : ICASB = ICAS + 1 688 : ELSE 689 : IF ((ICAS.EQ.2).OR.(ICAS.EQ.4)) THEN 690 : ICASB = ICAS -1 691 : ELSE 692 : ICASB = 2*ICAS-10+IHARDC 693 : ENDIF 694 : ENDIF 695 : IF (KONTAC(ICASB).EQ.0) THEN 695 : CALL CRITPE(SIGTP(IHARDC), SIGHP(IHARDC), RESIST(ICASB), 696 : &ISGNTC,XFFFS) 697 : IF (((XFFFS.GE.EPSIL(5)).AND.(ICASB.LE.2)).OR. 698 : &((XFFFS.GE.EPSIL(6)).AND.(ICASB.GE.3))) THEN 699 : IHARDB = IHARDC 700 : ISGNTCB = ISGNTC 701 : GOTO 1100 702 : ELSE 703 : IF (((ABS(XFFFS).LT.EPSIL(5)).AND.(ICASB.LE.2)).OR. 704 : &((ABS(XFFFS).LT.EPSIL(6)).AND.(ICASB.GE.3))) THEN 705 : KONTAC(ICASB)=1 706 : ENDIF 707 : ENDIF 708 : ELSE 708 : CALL CRITPE(SIGTP(IHARDC), XZER, FG(1), ISGNTC,XFFFS) 709 : IF (ICASB.LE.2) THEN 710 : ICASB = 5 711 : ELSE 712 : ICASB = 6 713 : ENDIF 714 : IF (XFFFS.GE.EPSIL(ICASB)) THEN 715 : IHARDB = IHARDC 716 : ISGNTCB = ISGNTC 717 : GOTO 1100 718 : ENDIF 719 : ENDIF 720 : C Autre signe signtc 721 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.3)) THEN 722 : ICASB = ICAS + 1 723 : ELSE 724 : IF ((ICAS.EQ.2).OR.(ICAS.EQ.4)) THEN 725 : ICASB = ICAS - 1 726 : ELSE 727 : ICASB = 2*ICAS-10+IHARDC 728 : ENDIF 729 : ENDIF 730 : IF (ICASB.LE.2) THEN 731 : ICASB = ICASB +2 732 : ELSE 733 : ICASB = ICASB -2 734 : ENDIF 735 : IF (KONTAC(ICASB).EQ.0) THEN 735 : CALL CRITPE(SIGTP(IHARDC), SIGHP(IHARDC), RESIST(ICASB), 736 : &-ISGNTC,XFFFS) 737 : IF (((XFFFS.GE.EPSIL(5)).AND.(ICASB.LE.2)).OR. 738 : &((XFFFS.GE.EPSIL(6)).AND.(ICASB.GE.3))) THEN 739 : IHARDB = IHARDC 740 : ISGNTCB = -ISGNTC 741 : GOTO 1100 742 : ELSE 743 : IF (((ABS(XFFFS).LT.EPSIL(5)).AND.(ICASB.LE.2)).OR. 744 : &((ABS(XFFFS).LT.EPSIL(6)).AND.(ICASB.GE.3))) THEN 745 : KONTAC(ICASB)=1 746 : ENDIF 747 : ENDIF 748 : ELSE 748 : CALL CRITPE(SIGTP(IHARDC), XZER, FG(2), -ISGNTC,XFFFS) 749 : IF (ICASB.LE.2) THEN 750 : ICASB = 5 751 : ELSE 752 : ICASB = 6 753 : ENDIF 754 : IF (XFFFS.GE.EPSIL(ICASB)) THEN 755 : IHARDB = IHARDC 756 : ISGNTCB = -ISGNTC 757 : GOTO 1100 758 : ENDIF 759 : ENDIF 760 : C Test du contact frontal 761 : IF (ICAS.LE.4) THEN 761 : CALL CRITPE(SIGTP(IHARD), XZER, FG(1), ISGNTC,XFFFS) 762 : IF (((XFFFS.GE.EPSIL(5)).AND.(ICAS.LE.2)).OR. 763 : &((XFFFS.GE.EPSIL(6)).AND.(ICAS.GE.3))) THEN 764 : GOTO 1040 765 : ELSE 766 : IF (((ABS(XFFFS).LT.EPSIL(5)).AND.(ICAS.LE.2)).OR. 767 : &((ABS(XFFFS).LT.EPSIL(6)).AND.(ICAS.GE.3))) THEN 768 : KONTAC(ICAS) = 1 769 : ENDIF 770 : ENDIF 771 : ENDIF 772 : C 773 : C Il faut faire suivre la petite surface si la grande se retrecit trop vite 774 : C 775 : IF (ICAS.LE.4) THEN 776 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.3)) THEN 777 : ICASB = ICAS +1 778 : ELSE 779 : ICASB = ICAS -1 780 : ENDIF 781 : XLIMIT = FG(1) -RESIST(ICASB) 782 : IF (ISGNTC*(SIGHP(IHARDC)-XLIMIT).GT.XZER) THEN 783 : SIGHP(IHARDC) = XLIMIT 784 : KONTAC(ICASB) = 1 785 : ENDIF 786 : IF (ICASB.LE.2) THEN 787 : ICASB = ICASB +2 788 : ELSE 789 : ICASB = ICASB -2 790 : ENDIF 791 : XLIMIT = FG(2) -RESIST(ICASB) 792 : IF (-ISGNTC*(SIGHP(IHARDC)-XLIMIT).GT.XZER) THEN 793 : SIGHP(IHARDC) = XLIMIT 794 : KONTAC(ICASB) = 1 795 : ENDIF 796 : ENDIF 797 : C 798 : C Mis a jour des indicateurs de contact 799 : C 800 : IF ((XKISO(1)*RESIST(5)-XKISO(2)*RESIST(6)).GT. 801 : &(RESIST(1)-RESIST(3))) THEN 802 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) KONTAC(ICAS+2)=0 803 : IF ((ICAS.EQ.3).OR.(ICAS.EQ.4)) KONTAC(ICAS-2)=0 804 : ELSE 805 : KONTAC(1) = 1 806 : KONTAC(2) = 1 807 : KONTAC(3) = 1 808 : KONTAC(4) = 1 809 : ENDIF 810 : C 811 : ICASB = ICAS 812 : GOTO 9999 813 : C 814 : C Cas ou on trouve le contact frontal 815 : C 816 : 1040 CONTINUE 817 : SIGHP(1) = SIGHP0(1) 818 : SIGHP(2) = SIGHP0(2) 819 : SIGT(1) = SIGI(1) + DSIGI(1) 820 : SIGT(2) = SIGI(2) + DSIGI(2) 821 : SIGT(4) = SIGI(4) + DSIGI(4) 822 : GOTO 3000 823 : C 824 : C Cas ou on viole un critere lateral 825 : C 826 : 1100 CONTINUE 827 : SIGHP(1) = SIGHP0(1) 828 : SIGHP(2) = SIGHP0(2) 829 : SIGT(1) = SIGI(1) + DSIGI(1) 830 : SIGT(2) = SIGI(2) + DSIGI(2) 831 : SIGT(4) = SIGI(4) + DSIGI(4) 832 : IF (ICAS.EQ.ICASB) THEN 833 : IF (ICAS.GE.5) THEN 834 : GOTO 1000 835 : ELSE 836 : KERRE = 2 837 : RETURN 838 : ENDIF 839 : ELSE 840 : GOTO 4000 841 : ENDIF 842 : C================================================================ 843 : C CAS 2 844 : C Un seul mecanisme de chaque surface 845 : C mais entree en contact des 2 mecanismes 846 : C================================================================ 847 : 3000 CONTINUE 848 : C 848 : CALL CONPRI (SIGT,CPHI0,SPHI0,SIGTP,CPHI,SPHI) 849 : C 850 : C 851 : C Calcul de la normale a l'estimation elastique et du jacobien initial 852 : C 853 : IF (IHARD.EQ.1) THEN 854 : IF (SIGTP(1).GE.SIGTP(2)) THEN 855 : ISGN2 = 1 856 : ELSE 857 : ISGN2 = -1 858 : ENDIF 859 : ELSE 860 : IF (IHARD.EQ.2) THEN 861 : IF (SIGTP(2).GE.SIGTP(1)) THEN 862 : ISGN2 = 1 863 : ELSE 864 : ISGN2 = -1 865 : ENDIF 866 : ENDIF 867 : ENDIF 868 : C 869 : 3010 CONTINUE 870 : C 871 : C Definition du segment 872 : C 873 : ITM = 30 874 : MN = 2 875 : SEGINI,QUASIN 876 : IT = -1 876 : CALL ZERO(D,1,MN) 877 : C 878 : C Calcul de xhh2 et xkiso 879 : C 880 : IF (ICAS.LE.2) THEN 880 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT, 881 : &XKISO(1),XHH2,KERRE) 881 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC, 882 : &XKISO(2),DUMM,KERRE) 883 : ELSE 883 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT, 884 : &XKISO(1),DUMM,KERRE) 884 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC, 885 : &XKISO(2),XHH2,KERRE) 886 : ENDIF 887 : EPSI(1) = EPSI0(1) 888 : XKISO0(1) = XKISO(1) 889 : EPSI(2) = EPSI0(2) 890 : XKISO0(2) = XKISO(2) 891 : XHH2I = XHH2 892 : EPSII(1) = EPSI0(1) 893 : XKISOI(1) = XKISO(1) 894 : EPSII(2) = EPSI0(2) 895 : XKISOI(2) = XKISO(2) 896 : C 897 : DO 3011,IE1=1,6 898 : FFF(IE1) = XZER 899 : 3011 CONTINUE 900 : C 901 : IF (ISGNTC.GT.0) THEN 902 : FG(1) = XKISO(1)*RESIST(5) 903 : ELSE 904 : FG(1) = XKISO(2)*RESIST(6) 905 : ENDIF 905 : CALL CRITPE (SIGTP(IHARD),SIGHP0(IHARD) 906 : & ,RESIST(ICAS),ISGNTC,FFF(ICAS)) 906 : CALL CRITPE (SIGTP(IHARD), XZER, FG(1), ISGNTC, XFFFS) 907 : C 907 : CALL NORMBL(SIGT,RESIST(ICAS),ISGNTC,ISGN2,1,XNOR1) 908 : XJ00(1,1) = -XXCCYY(XNOR1, XNOR1, YOUN, XNU)-XHH1 909 : XJ00(2,2) = -XXYY(XNOR1,DSIGI) 910 : XJ00(1,2) = -XXYY(XNOR1,DSIGI) 911 : XJ00(2,1) = -XXCCYY(XNOR1, XNOR1, YOUN, XNU) 912 : IF (ISGNTC.LT.0) THEN 913 : XJ00(2,1)=XJ00(2,1)+RESIST(6)*XHH2*DEGRAD*ISGNTC*RESIST(ICAS) 914 : ENDIF 915 : R(1) = FFF(ICAS) 916 : R(2) = XFFFS 917 : GAMBDA(1)=XZER 918 : GAMBDA(2)=XZER 919 : GAMBDA(3)=XZER 920 : GAMBDA(4)=XZER 921 : GAMBDA(5)=XZER 922 : GAMBDA(6)=XZER 923 : C 924 : C Inverse du jacobien 925 : C 926 : XH0(1,1) = XJ00(1,1) 927 : XH0(1,2) = XJ00(1,2) 928 : XH0(2,1) = XJ00(2,1) 929 : XH0(2,2) = XJ00(2,2) 929 : CALL INVALM (XH0, MN , MN , IRD, 1.D-12) 930 : IF (IRD.NE.0) THEN 931 : SEGSUP,QUASIN 932 : GOTO 6000 933 : ENDIF 934 : C 935 : C On rentre dans les iterations internes 936 : C 937 : DO I=0,ITM 938 : C 939 : C Appel de Broyden 940 : C 940 : CALL BROYDE (QUASIN) 941 : C 942 : GAMBDA(ICAS) = D(1) 943 : SPLITT = D(2) 944 : C 945 : C Calcul du phi et des termes de la matrice A et A-1 946 : C 947 : SITPRO(1) = SIGI(1) + ((UN-SPLITT)*DSIGI(1)) 948 : & - UNDEMI*(GAMBDA(1) + GAMBDA(2))*CPI 949 : & + UNDEMI*(GAMBDA(3) + GAMBDA(4))*CPI 950 : SITPRO(2) = SIGI(2) + ((UN-SPLITT)*DSIGI(2)) 951 : & - UNDEMI*(GAMBDA(1) + GAMBDA(2))*CPI 952 : & + UNDEMI*(GAMBDA(3) + GAMBDA(4))*CPI 953 : SITPRO(4) = SIGI(4) + (UN-SPLITT)*DSIGI(4) 954 : C 955 : SIGHP(1) = SIGHP0(1) + XHH1*(GAMBDA(1)-GAMBDA(3)) 956 : SIGHP(2) = SIGHP0(2) + XHH1*(GAMBDA(2)-GAMBDA(4)) 957 : EPSI(1) = EPSI0(1) 958 : WORK = GAMBDA(ICAS)*ISGNTC*RESIST(ICAS) 959 : EPSI(2) = EPSI0(2) + DEGRAD*WORK 959 : CALL VALEUR(EPSI(2),TRAC(IPC),NCURVC,XLCAC, 960 : &XKISO(2),XHH2,KERRE) 961 : C 962 : PHI = ISGNTC*(-UNDEMI*(SITPRO(1)+SITPRO(2))+SIGHP(IHARD) 963 : & + RESIST(ICAS)) 964 : C 965 : C 966 : SOMLAM = GAMBDA(1) + GAMBDA(2) + GAMBDA(3) + GAMBDA(4) 967 : AINV(1) = (DEUX*PHI+GGG*SOMLAM)/(DEUX*PHI+DEUX*GGG*SOMLAM) 968 : AINV(2) = GGG*SOMLAM/(DEUX*PHI + DEUX*GGG*SOMLAM) 969 : AINV(3) = PHI/(PHI + GGG*SOMLAM) 970 : C 971 : C Calcul des contraintes a l'aide de AINV 972 : C 973 : SIGT(1) = AINV(1)*SITPRO(1)+AINV(2)*SITPRO(2) 974 : SIGT(2) = AINV(2)*SITPRO(1)+AINV(1)*SITPRO(2) 975 : SIGT(4) = AINV(3)*SITPRO(4) 976 : C 977 : C Nouvelles contraintes principales 978 : C 978 : CALL CONPRI (SIGT,CPHI0,SPHI0,SIGTP,CPHI,SPHI) 979 : C 980 : NMEC = 0 981 : C 982 : IF (ISGNTC.GT.0) THEN 983 : FG(1) = XKISO(1)*RESIST(5) 984 : ELSE 985 : FG(1) = (XKISOI(2) + XHH2I*(EPSI(2)-EPSII(2)))*RESIST(6) 986 : ENDIF 986 : CALL CRITPE (SIGTP(IHARD),SIGHP(IHARD) 987 : & ,RESIST(ICAS),ISGNTC,FFF(ICAS)) 988 : R(1) = + FFF(ICAS) 989 : IF (ABS(FFF(ICAS)).GE.EPSIL(ICAS)) THEN 990 : NMEC = NMEC + 1 991 : ENDIF 991 : CALL CRITPE (SIGTP(IHARD), XZER, FG(1), ISGNTC, XFFFS) 992 : R(2) = XFFFS 993 : IF (ICAS.GE.3) THEN 994 : R(2) = R(2) + RESIST(6)*(XKISO(2) - XKISOI(2) 995 : & - (EPSI(2) - EPSII(2))*XHH2I) 996 : IF (ABS(EPSI(2) - EPSII(2)).LT.EPSILO2) THEN 997 : XHH2I=XHH2 998 : ELSE 999 : XHH2I = (XKISO(2) - XKISOI(2))/ 1000 : & (EPSI(2) - EPSII(2)) 1001 : ENDIF 1002 : EPSII(2) = EPSI(2) 1003 : XKISOI(2) = XKISO(2) 1004 : ENDIF 1005 : IF (((ABS(XFFFS).GE.EPSIL(5)).AND.(ICAS.LE.2)).OR. 1006 : &((ABS(XFFFS).GE.EPSIL(6)).AND.(ICAS.GE.3))) THEN 1007 : NMEC = NMEC + 1 1008 : ENDIF 1009 : C 1010 : C Test de convergence 1011 : C 1012 : IF (NMEC.EQ.0) GOTO 3020 1013 : ENDDO 1014 : C 1015 : C Pas de convergence 1016 : C 1017 : SEGSUP,QUASIN 1018 : GOTO 6000 1019 : C 1020 : 3020 CONTINUE 1021 : C 1022 : C On verifie que la convergence ne s'est pas faite sur un point debile 1023 : C 1024 : SEGSUP,QUASIN 1025 : IF (((SPLITT-1.).LE.EPSILO).AND.(SPLITT.GE.(-1.*EPSILO))) THEN 1026 : C 1027 : IF (IHARD.EQ.2) THEN 1028 : IHARDC=1 1029 : ELSE 1030 : IHARDC=2 1031 : ENDIF 1032 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.3)) ICASC = ICAS + 1 1033 : IF ((ICAS.EQ.2).OR.(ICAS.EQ.4)) ICASC = ICAS - 1 1034 : C 1035 : IF (ICAS.LE.2) THEN 1036 : FG(1) = XKISO(1)*RESIST(5) 1037 : FG(2) = XKISO(2)*RESIST(6) 1038 : ELSE 1039 : FG(1) = XKISO(2)*RESIST(6) 1040 : FG(2) = XKISO(1)*RESIST(5) 1041 : ENDIF 1042 : C 1043 : C Test des contacts lateraux 1044 : C Meme signe signtc 1045 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.3)) THEN 1046 : ICASB = ICAS + 1 1047 : ELSE 1048 : ICASB = ICAS -1 1049 : ENDIF 1050 : IF (KONTAC(ICASB).EQ.0) THEN 1050 : CALL CRITPE(SIGTP(IHARDC), SIGHP(IHARDC), 1051 : & RESIST(ICASB), ISGNTC,XFFFS) 1052 : IF (((XFFFS.GE.EPSIL(5)).AND.(ICASB.LE.2)).OR. 1053 : &((XFFFS.GE.EPSIL(6)).AND.(ICASB.GE.3))) THEN 1054 : IHARDB = IHARDC 1055 : ISGNTCB = ISGNTC 1056 : GOTO 3100 1057 : ELSE 1058 : IF (((ABS(XFFFS).LT.EPSIL(5)).AND.(ICASB.LE.2)).OR. 1059 : &((ABS(XFFFS).LT.EPSIL(6)).AND.(ICASB.GE.3))) THEN 1060 : KONTAC(ICASB)=1 1061 : ENDIF 1062 : ENDIF 1063 : ELSE 1063 : CALL CRITPE(SIGTP(IHARDC), XZER, FG(1), ISGNTC,XFFFS) 1064 : IF (ICASB.LE.2) THEN 1065 : ICASB = 5 1066 : ELSE 1067 : ICASB = 6 1068 : ENDIF 1069 : IF (XFFFS.GE.EPSIL(ICASB)) THEN 1070 : IHARDB = IHARDC 1071 : ISGNTCB = -ISGNTC 1072 : GOTO 3100 1073 : ENDIF 1074 : ENDIF 1075 : C Autre signe signtc 1076 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.3)) THEN 1077 : ICASB = ICAS + 1 1078 : ELSE 1079 : ICASB = ICAS -1 1080 : ENDIF 1081 : IF (ICASB.LE.2) THEN 1082 : ICASB = ICASB +2 1083 : ELSE 1084 : ICASB = ICASB -2 1085 : ENDIF 1086 : IF (KONTAC(ICASB).EQ.0) THEN 1086 : CALL CRITPE(SIGTP(IHARDC), SIGHP(IHARDC),RESIST(ICASB), 1087 : & -ISGNTC,XFFFS) 1088 : IF (((XFFFS.GE.EPSIL(5)).AND.(ICASB.LE.2)).OR. 1089 : &((XFFFS.GE.EPSIL(6)).AND.(ICASB.GE.3))) THEN 1090 : IHARDB = IHARDC 1091 : ISGNTCB = ISGNTC 1092 : GOTO 3100 1093 : ELSE 1094 : IF (((ABS(XFFFS).LT.EPSIL(5)).AND.(ICASB.LE.2)).OR. 1095 : &((ABS(XFFFS).LT.EPSIL(6)).AND.(ICASB.GE.3))) THEN 1096 : KONTAC(ICASB)=1 1097 : ENDIF 1098 : ENDIF 1099 : ELSE 1099 : CALL CRITPE(SIGTP(IHARDC), XZER, FG(2), -ISGNTC,XFFFS) 1100 : IF (ICASB.LE.2) THEN 1101 : ICASB = 5 1102 : ELSE 1103 : ICASB = 6 1104 : ENDIF 1105 : IF (XFFFS.GE.EPSIL(ICASB)) THEN 1106 : IHARDB = IHARDC 1107 : ISGNTCB = -ISGNTC 1108 : GOTO 3100 1109 : ENDIF 1110 : ENDIF 1111 : C 1112 : C Il faut faire suivre la petite surface si la grande se retrecit trop vite 1113 : C 1114 : IF (ICAS.LE.4) THEN 1115 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.3)) THEN 1116 : ICASB = ICAS +1 1117 : ELSE 1118 : ICASB = ICAS -1 1119 : ENDIF 1120 : XLIMIT = FG(1) -RESIST(ICASB) 1121 : IF (ISGNTC*(SIGHP(IHARDC)-XLIMIT).GT.XZER) THEN 1122 : SIGHP(IHARDC) = XLIMIT 1123 : KONTAC(ICASB) = 1 1124 : ENDIF 1125 : IF (ICASB.LE.2) THEN 1126 : ICASB = ICASB +2 1127 : ELSE 1128 : ICASB = ICASB -2 1129 : ENDIF 1130 : XLIMIT = FG(2) -RESIST(ICASB) 1131 : IF (-ISGNTC*(SIGHP(IHARDC)-XLIMIT).GT.XZER) THEN 1132 : SIGHP(IHARDC) = XLIMIT 1133 : KONTAC(ICASB) = 1 1134 : ENDIF 1135 : ENDIF 1136 : C 1137 : C 1138 : C Mis a jour des indicateurs de contact 1139 : C 1140 : KONTAC(ICAS) = 1 1141 : IF ((XKISO(1)*RESIST(5)-XKISO(2)*RESIST(6)).GT. 1142 : &(RESIST(1)-RESIST(3))) THEN 1143 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) KONTAC(ICAS+2)=0 1144 : IF ((ICAS.EQ.3).OR.(ICAS.EQ.4)) KONTAC(ICAS-2)=0 1145 : ELSE 1146 : KONTAC(1) = 1 1147 : KONTAC(2) = 1 1148 : KONTAC(3) = 1 1149 : KONTAC(4) = 1 1150 : ENDIF 1151 : SIGI(1) = SIGT(1) 1152 : SIGI(2) = SIGT(2) 1153 : SIGI(4) = SIGT(4) 1154 : DSIGI(1) = DSIGI(1)*SPLITT 1155 : DSIGI(2) = DSIGI(2)*SPLITT 1156 : DSIGI(4) = DSIGI(4)*SPLITT 1157 : SIGT(1) = SIGI(1) + DSIGI(1) 1158 : SIGT(2) = SIGI(2) + DSIGI(2) 1159 : SIGT(4) = SIGI(4) + DSIGI(4) 1160 : SIGHP0(1) = SIGHP(1) 1161 : SIGHP0(2) = SIGHP(2) 1162 : EPSI0(1) = EPSI(1) 1163 : EPSI0(2) = EPSI(2) 1164 : CPHI0 = CPHI 1165 : SPHI0 = SPHI 1166 : NPASS2 = 0 1167 : IF (ICAS.LE.2) THEN 1168 : ICAS = 5 1169 : ELSE 1170 : ICAS = 6 1171 : ENDIF 1172 : GOTO 1000 1173 : C 1174 : C 1175 : ELSE 1176 : GOTO 6000 1177 : ENDIF 1178 : C 1179 : C Cas ou on viole un critere lateral 1180 : C 1181 : 3100 CONTINUE 1182 : SIGHP(1) = SIGHP0(1) 1183 : SIGHP(2) = SIGHP0(2) 1184 : SIGT(1) = SIGI(1) + DSIGI(1) 1185 : SIGT(2) = SIGI(2) + DSIGI(2) 1186 : SIGT(4) = SIGI(4) + DSIGI(4) 1186 : CALL CONPRI ( SIGT ,CPHI0,SPHI0,SIGTP,CPHI,SPHI) 1186 : CALL VALEUR(EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1), 1187 : &DUMM,KERRE) 1187 : CALL VALEUR(EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2), 1188 : &DUMM,KERRE) 1188 : CALL CRITPE(SIGTP(IHARD), SIGHP(IHARD), RESIST(ICAS), 1189 : & ISGNTC, FFF(ICAS)) 1190 : IF (ICASB.LE.4) THEN 1190 : CALL CRITPE(SIGTP(IHARDB), SIGHP(IHARDB), RESIST(ICASB), 1191 : & ISGNTCB, FFF(ICASB)) 1192 : ELSE 1193 : FG(1) = XKISO(ICASB-4)*RESIST(ICASB) 1193 : CALL CRITPE (SIGTP(IHARDB), XZER, FG(1), 1194 : & ISGNTCB,FFF(ICASB)) 1195 : ENDIF 1196 : IF (ICAS.LE.2) THEN 1197 : FG(1) = XKISO(1)*RESIST(5) 1198 : ELSE 1199 : FG(1) = XKISO(2)*RESIST(6) 1200 : ENDIF 1200 : CALL CRITPE (SIGTP(IHARD), XZER, FG(1), ISGNTC, XFFFS) 1201 : GOTO 5000 1202 : C================================================================ 1203 : C CAS 4 1204 : C Deux mecanismes 1205 : C Petite et Grande surfaces (18/07) 1206 : C================================================================ 1207 : 4000 CONTINUE 1208 : C 1208 : CALL CONPRI ( SIGT,CPHI0,SPHI0,SIGTP,CPHI,SPHI) 1209 : C 1210 : C 1211 : C Calcul de la normale a l'estimation elastique et du jacobien initial 1212 : C 1213 : IF (IHARD.EQ.1) THEN 1214 : IF (SIGTP(1).GE.SIGTP(2)) THEN 1215 : ISGN2 = 1 1216 : ELSE 1217 : ISGN2 = -1 1218 : ENDIF 1219 : ELSE 1220 : IF (IHARD.EQ.2) THEN 1221 : IF (SIGTP(2).GE.SIGTP(1)) THEN 1222 : ISGN2 = 1 1223 : ELSE 1224 : ISGN2 = -1 1225 : ENDIF 1226 : ENDIF 1227 : ENDIF 1228 : C 1229 : IF (IHARDB.EQ.1) THEN 1230 : IF (SIGTP(1).GE.SIGTP(2)) THEN 1231 : ISGN2B = 1 1232 : ELSE 1233 : ISGN2B = -1 1234 : ENDIF 1235 : ELSE 1236 : IF (IHARDB.EQ.2) THEN 1237 : IF (SIGTP(1).GE.SIGTP(2)) THEN 1238 : ISGN2B = -1 1239 : ELSE 1240 : ISGN2B = 1 1241 : ENDIF 1242 : ENDIF 1243 : ENDIF 1244 : C 1245 : C 1246 : C 1247 : 4010 CONTINUE 1248 : C 1249 : C Definition du segment 1250 : C 1251 : ITM = 30 1252 : MN = 2 1253 : SEGINI,QUASIN 1254 : IT = -1 1254 : CALL ZERO(D,1,MN) 1255 : C 1256 : C XKISO et XHH2 1257 : C 1258 : IF (ICAS.LE.4) THEN 1259 : IF (ICASB.LE.4) THEN 1259 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT, 1260 : &XKISO(1),DUMM,KERRE) 1260 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC, 1261 : &XKISO(2),DUMM,KERRE) 1262 : ELSE 1263 : IF (ICASB.EQ.5) THEN 1263 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1264 : & ,XHH2B,KERRE) 1264 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1265 : & ,DUMM,KERRE) 1266 : ELSE 1266 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1267 : & ,DUMM,KERRE) 1267 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1268 : & ,XHH2B,KERRE) 1269 : ENDIF 1270 : ENDIF 1271 : ELSE 1272 : IF (ICASB.LE.4) THEN 1273 : IF (ICAS.EQ.5) THEN 1273 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1274 : & ,XHH2,KERRE) 1274 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1275 : & ,DUMM,KERRE) 1276 : ELSE 1276 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1277 : & ,DUMM,KERRE) 1277 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1278 : & ,XHH2,KERRE) 1279 : ENDIF 1280 : ELSE 1281 : IF (ICASB.EQ.5) THEN 1281 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1282 : & ,XHH2B,KERRE) 1282 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1283 : & ,XHH2,KERRE) 1284 : ELSE 1284 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1285 : & ,XHH2,KERRE) 1285 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1286 : & ,XHH2B,KERRE) 1287 : ENDIF 1288 : ENDIF 1289 : ENDIF 1290 : EPSI(1) = EPSI0(1) 1291 : EPSI(2) = EPSI0(2) 1292 : XHH2I = XHH2 1293 : XHH2IB = XHH2B 1294 : XKISO0(1) = XKISO(1) 1295 : XKISO0(2) = XKISO(2) 1296 : EPSII(1) = EPSI0(1) 1297 : EPSII(2) = EPSI0(2) 1298 : XKISOI(1) = XKISO(1) 1299 : XKISOI(2) = XKISO(2) 1300 : C 1301 : DO 4011,IE1=1,6 1302 : FFF(IE1) = XZER 1303 : 4011 CONTINUE 1304 : C 1305 : IF (ICAS.LE.4) THEN 1305 : CALL CRITPE (SIGTP(IHARD),SIGHP0(IHARD),RESIST(ICAS), 1306 : & ISGNTC,FFF(ICAS)) 1307 : ELSE 1308 : FG(1) = XKISOI(ICAS-4)*RESIST(ICAS) 1308 : CALL CRITPE (SIGTP(IHARD),XZER,FG(1),ISGNTC,FFF(ICAS)) 1309 : ENDIF 1310 : C 1311 : IF (ICASB.LE.4) THEN 1311 : CALL CRITPE (SIGTP(IHARDB),SIGHP0(IHARDB),RESIST(ICASB), 1312 : & ISGNTCB,FFF(ICASB)) 1313 : ELSE 1314 : FG(1) = XKISOI(ICASB-4)*RESIST(ICASB) 1314 : CALL CRITPE (SIGTP(IHARDB),XZER,FG(1),ISGNTCB,FFF(ICASB)) 1315 : ENDIF 1316 : C 1317 : C Jcobien initial 1318 : C 1318 : CALL NORMBL(SIGT,RESIST(ICAS),ISGNTC,ISGN2,1,XNOR1) 1318 : CALL NORMBL(SIGT,RESIST(ICASB),ISGNTCB,ISGN2B,1,XNOR2) 1319 : C 1320 : XJ00(1,1) = -XXCCYY(XNOR1, XNOR1, YOUN, XNU) 1321 : XJ00(1,2) = -XXCCYY(XNOR1, XNOR2, YOUN, XNU) 1322 : XJ00(2,1) = -XXCCYY(XNOR2, XNOR1, YOUN, XNU) 1323 : XJ00(2,2) = -XXCCYY(XNOR2, XNOR2, YOUN, XNU) 1324 : C 1325 : C Premier Mecanisme 1326 : C 1327 : IF (ICAS.LE.4) THEN 1328 : XJ00(1,1) = XJ00(1,1) - XHH1 1329 : ELSE 1330 : XJ00(1,1) = XJ00(1,1)-ISGNTC*RESIST(ICAS)*XHH2 1331 : IF ((ICASB.LE.4).AND.(ICAS.EQ.6)) THEN 1332 : XJ00(1,2) = XJ00(1,2) 1333 : & +RESIST(6)*XHH2*DEGRAD*ISGNTCB*RESIST(ICASB) 1334 : ENDIF 1335 : ENDIF 1336 : C 1337 : C Second Mecanisme 1338 : C 1339 : IF (ICASB.LE.4) THEN 1340 : XJ00(2,2) = XJ00(2,2) - XHH1 1341 : ELSE 1342 : XJ00(2,2) = XJ00(2,2)-ISGNTCB*RESIST(ICASB)*XHH2B 1343 : IF ((ICAS.LE.4).AND.(ICASB.EQ.6)) THEN 1344 : XJ00(2,1) = XJ00(2,1) 1345 : & +RESIST(6)*XHH2B*DEGRAD*ISGNTC*RESIST(ICAS) 1346 : ENDIF 1347 : ENDIF 1348 : C 1349 : C 1350 : R(1) = + FFF(ICAS) 1351 : R(2) = + FFF(ICASB) 1352 : GAMBDA(1) = XZER 1353 : GAMBDA(2) = XZER 1354 : GAMBDA(3) = XZER 1355 : GAMBDA(4) = XZER 1356 : GAMBDA(5) = XZER 1357 : GAMBDA(6) = XZER 1358 : C 1359 : C Inverse du jacobien 1360 : C 1361 : XH0(1,1) = XJ00(1,1) 1362 : XH0(1,2) = XJ00(1,2) 1363 : XH0(2,1) = XJ00(2,1) 1364 : XH0(2,2) = XJ00(2,2) 1364 : CALL INVALM (XH0, MN , MN , IRD, 1.D-12) 1365 : IF (IRD.NE.0) THEN 1366 : SEGSUP,QUASIN 1367 : GOTO 6000 1368 : ENDIF 1369 : C 1370 : C On rentre dans les iterations internes 1371 : C 1372 : DO I=0,ITM 1373 : C 1374 : C Appel de Broyden 1375 : C 1375 : CALL BROYDE (QUASIN) 1376 : C 1377 : GAMBDA(ICAS) = D(1) 1378 : GAMBDA(ICASB) = D(2) 1379 : C 1380 : C Calcul du phi et des termes de la matrice A et A-1 1381 : C 1382 : SITPRO(1) = SIGI(1) + DSIGI(1) 1383 : & - UNDEMI*(GAMBDA(1) + GAMBDA(2) + GAMBDA(5))*CPI 1384 : & + UNDEMI*(GAMBDA(3) + GAMBDA(4) + GAMBDA(6))*CPI 1385 : SITPRO(2) = SIGI(2) + DSIGI(2) 1386 : & - UNDEMI*(GAMBDA(1) + GAMBDA(2) + GAMBDA(5))*CPI 1387 : & + UNDEMI*(GAMBDA(3) + GAMBDA(4) + GAMBDA(6))*CPI 1388 : SITPRO(4) = SIGI(4) + DSIGI(4) 1389 : C 1390 : C PHI et AINV 1391 : C 1392 : SIGHP(1) = SIGHP0(1) + XHH1*(GAMBDA(1)-GAMBDA(3)) 1393 : SIGHP(2) = SIGHP0(2) + XHH1*(GAMBDA(2)-GAMBDA(4)) 1394 : EPSI(1) = EPSI0(1) + GAMBDA(5) 1395 : IF (ICAS.LE.4) THEN 1396 : WORK = GAMBDA(ICAS)*RESIST(ICAS)*ISGNTC 1397 : ELSE 1398 : WORK = XZER 1399 : ENDIF 1400 : IF (ICASB.LE.4) THEN 1401 : WORK = WORK + GAMBDA(ICASB)*RESIST(ICASB)*ISGNTCB 1402 : ENDIF 1403 : EPSI(2) = EPSI0(2) + GAMBDA(6) + DEGRAD*WORK 1404 : C 1405 : IF (ICAS.LE.4) THEN 1406 : IF (ICASB.LE.4) THEN 1406 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT, 1407 : &XKISO(1),DUMM,KERRE) 1407 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC, 1408 : &XKISO(2),DUMM,KERRE) 1409 : ELSE 1410 : IF (ICASB.EQ.5) THEN 1410 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1411 : & ,XHH2B,KERRE) 1411 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1412 : & ,DUMM,KERRE) 1413 : ELSE 1413 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1414 : & ,DUMM,KERRE) 1414 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1415 : & ,XHH2B,KERRE) 1416 : ENDIF 1417 : ENDIF 1418 : ELSE 1419 : IF (ICASB.LE.4) THEN 1420 : IF (ICAS.EQ.5) THEN 1420 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1421 : & ,XHH2,KERRE) 1421 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1422 : & ,DUMM,KERRE) 1423 : ELSE 1423 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1424 : & ,DUMM,KERRE) 1424 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1425 : & ,XHH2,KERRE) 1426 : ENDIF 1427 : ELSE 1428 : IF (ICASB.EQ.5) THEN 1428 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1429 : & ,XHH2B,KERRE) 1429 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1430 : & ,XHH2,KERRE) 1431 : ELSE 1431 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1432 : & ,XHH2,KERRE) 1432 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1433 : & ,XHH2B,KERRE) 1434 : ENDIF 1435 : ENDIF 1436 : ENDIF 1437 : C 1438 : C Calcul des phi 1439 : C 1440 : IF (ICAS.LE.4) THEN 1441 : PHI1=ISGNTC*(-UNDEMI*(SITPRO(1) + SITPRO(2)) 1442 : & + SIGHP(IHARD) + RESIST(ICAS)) 1443 : ELSE 1444 : PHI1= ISGNTC*(-UNDEMI*(SITPRO(1)+SITPRO(2)) 1445 : & +(XKISOI(ICAS-4)+XHH2I*(EPSI(ICAS-4)-EPSII(ICAS-4))) 1446 : & *RESIST(ICAS)) 1447 : ENDIF 1448 : C 1449 : IF (ICASB.LE.4) THEN 1450 : PHI2=ISGNTCB*(-UNDEMI*(SITPRO(1) + SITPRO(2)) 1451 : & + SIGHP(IHARDB)+ RESIST(ICASB)) 1452 : ELSE 1453 : PHI2=ISGNTCB*(-UNDEMI*(SITPRO(1)+SITPRO(2)) 1454 : & +(XKISOI(ICASB-4)+XHH2IB*(EPSI(ICASB-4)-EPSII(ICASB-4))) 1455 : & *RESIST(ICASB)) 1456 : ENDIF 1457 : C 1458 : C Test sur phi1 et phi2 1459 : C 1460 : IF ((ABS(PHI2)).GT.(ABS(EPSILO*RESIST(ICASB)))) THEN 1461 : SOMLAM = GAMBDA(ICAS) + GAMBDA(ICASB)*PHI1/PHI2 1462 : PHI = PHI1 1463 : ELSE 1464 : IF ((ABS(PHI1)).GT.(ABS(EPSILO*RESIST(ICAS)))) THEN 1465 : SOMLAM = GAMBDA(ICASB) + GAMBDA(ICAS)*PHI2/PHI1 1466 : PHI = PHI2 1467 : ELSE 1468 : SOMLAM = GAMBDA(ICASB) + GAMBDA(ICAS) 1469 : PHI = PHI1 1470 : ENDIF 1471 : ENDIF 1472 : C 1473 : AINV(1) = (DEUX*PHI+GGG*SOMLAM)/(DEUX*PHI+DEUX*GGG*SOMLAM) 1474 : AINV(2) = (GGG*SOMLAM)/(DEUX*PHI + DEUX*GGG*SOMLAM) 1475 : AINV(3) = PHI/(PHI + GGG*SOMLAM) 1476 : C 1477 : C Calcul des contraintes a l'aide de AINV 1478 : C 1479 : SIGT(1) = AINV(1)*SITPRO(1)+AINV(2)*SITPRO(2) 1480 : SIGT(2) = AINV(2)*SITPRO(1)+AINV(1)*SITPRO(2) 1481 : SIGT(4) = AINV(3)*SITPRO(4) 1482 : C 1483 : C Nouvelles contraintes principales 1484 : C 1484 : CALL CONPRI ( SIGT,CPHI0,SPHI0,SIGTP,CPHI,SPHI) 1485 : C 1486 : C 1487 : IF (ICAS.LE.4) THEN 1487 : CALL CRITPE (SIGTP(IHARD),SIGHP(IHARD),RESIST(ICAS), 1488 : & ISGNTC,FFF(ICAS)) 1489 : R(1)= + FFF(ICAS) 1490 : ELSE 1491 : FG(1) = (XKISOI(ICAS-4)+XHH2I*(EPSI(ICAS-4)-EPSII(ICAS-4))) 1492 : & *RESIST(ICAS) 1492 : CALL CRITPE (SIGTP(IHARD),XZER,FG(1),ISGNTC,FFF(ICAS)) 1493 : R(1) = + FFF(ICAS) 1494 : & - ISGNTC*RESIST(ICAS) 1495 : & * (XKISO(ICAS-4) - XKISOI(ICAS-4) 1496 : & -(EPSI(ICAS-4) - EPSII(ICAS-4))*XHH2I) 1497 : IF (ABS(EPSI(ICAS-4) - EPSII(ICAS-4)).LT.EPSILO2) THEN 1498 : XHH2I=XHH2 1499 : ELSE 1500 : XHH2I = (XKISO(ICAS-4) - XKISOI(ICAS-4))/ 1501 : & (EPSI(ICAS-4) - EPSII(ICAS-4)) 1502 : ENDIF 1503 : EPSII(ICAS-4) = EPSI(ICAS-4) 1504 : XKISOI(ICAS-4) = XKISO(ICAS-4) 1505 : ENDIF 1506 : C 1507 : IF (ICASB.LE.4) THEN 1507 : CALL CRITPE (SIGTP(IHARDB),SIGHP(IHARDB),RESIST(ICASB), 1508 : & ISGNTCB,FFF(ICASB)) 1509 : R(2)= + FFF(ICASB) 1510 : ELSE 1511 : FG(1) = (XKISOI(ICASB-4)+XHH2IB*(EPSI(ICASB-4) - 1512 : & EPSII(ICASB-4)))*RESIST(ICASB) 1512 : CALL CRITPE (SIGTP(IHARDB),XZER,FG(1),ISGNTCB,FFF(ICASB)) 1513 : R(2) = + FFF(ICASB) 1514 : & - ISGNTCB*RESIST(ICASB) 1515 : & * (XKISO(ICASB-4) - XKISOI(ICASB-4) 1516 : & -(EPSI(ICASB-4) - EPSII(ICASB-4))*XHH2IB) 1517 : IF (ABS(EPSI(ICASB-4) - EPSII(ICASB-4)).LT.EPSILO2) THEN 1518 : XHH2IB=XHH2B 1519 : ELSE 1520 : XHH2IB = (XKISO(ICASB-4) - XKISOI(ICASB-4))/ 1521 : & (EPSI(ICASB-4) - EPSII(ICASB-4)) 1522 : ENDIF 1523 : XKISOI(ICASB-4) = XKISO(ICASB-4) 1524 : EPSII(ICASB-4) = EPSI(ICASB-4) 1525 : ENDIF 1526 : C 1527 : C Test de convergence 1528 : C 1529 : IF ((ABS(FFF(ICAS)).LE.EPSIL(ICAS)) 1530 : &.AND.(ABS(FFF(ICASB)).LE.EPSIL(ICASB)).AND. 1531 : &(ABS(R(1)).LE.EPSIL(ICAS)).AND.(ABS(R(2)).LE.EPSIL(ICASB))) THEN 1532 : GOTO 4020 1533 : ENDIF 1534 : ENDDO 1535 : C 1536 : C Pas de convergence 1537 : C 1538 : SEGSUP,QUASIN 1539 : GOTO 6000 1540 : C 1541 : 4020 CONTINUE 1542 : SEGSUP,QUASIN 1543 : C 1544 : IF ((GAMBDA(ICAS).LT.EPSILO2).OR. 1545 : &(GAMBDA(ICASB).LT.EPSILO2)) THEN 1546 : 4030 CONTINUE 1547 : IF (NPASS.LE.5) THEN 1548 : NPASS = NPASS + 1 1549 : ELSE 1550 : GOTO 6000 1551 : ENDIF 1552 : IF ((GAMBDA(ICAS).LT.EPSILO2).AND. 1553 : &(GAMBDA(ICASB).LT.EPSILO2)) THEN 1554 : C Ce cas n'arrive que si les normales ont ete mal estimes: 1555 : C (increment elastique trop grand) --> on bissecte 1556 : GOTO 6000 1557 : ELSE 1558 : IF (GAMBDA(ICAS).LT.EPSILO2) THEN 1559 : ICAS = ICASB 1560 : IHARD = IHARDB 1561 : ISGNTC = ISGNTCB 1562 : ELSE 1563 : CONTINUE 1564 : ENDIF 1565 : ENDIF 1566 : SIGHP(1) = SIGHP0(1) 1567 : SIGHP(2) = SIGHP0(2) 1568 : SIGT(1) = SIGI(1) + DSIGI(1) 1569 : SIGT(2) = SIGI(2) + DSIGI(2) 1570 : SIGT(4) = SIGI(4) + DSIGI(4) 1571 : GOTO 1000 1572 : ENDIF 1573 : C On teste les mecanismes de la grande surface 1574 : C Si un des criteres est viole (et donc s'il n'y a pas contact) 1575 : C il faut trouver le point de contact 1576 : C 1576 : CALL CONPRI ( SIGT ,CPHI0,SPHI0,SIGTP,CPHI,SPHI) 1577 : C 1578 : NCON = 0 1579 : FINT1 = 0. 1580 : FINT2 = 0. 1581 : C 1582 : XFFFS = XZER 1583 : XFFFSB = XZER 1584 : IF (ICAS.LE.4) THEN 1585 : IF (ICAS.LE.2) THEN 1586 : FG(1) = XKISO(1)*RESIST(5) 1587 : ELSE 1588 : FG(1) = XKISO(2)*RESIST(6) 1589 : ENDIF 1589 : CALL CRITPE (SIGTP(IHARD),XZER,FG(1),ISGNTC,XFFFS) 1590 : IF (((XFFFS.GE.EPSIL(5)).AND.(ICAS.LE.2)).OR. 1591 : & ((XFFFS.GE.EPSIL(6)).AND.(ICAS.GE.3))) THEN 1592 : NCON = NCON +1 1593 : ELSE 1594 : IF (((ICAS.LE.2).AND.(ABS(XFFFS).LT.EPSIL(5))).OR. 1595 : &((ICAS.GE.3).AND.(ABS(XFFFS).LT.EPSIL(6)))) THEN 1596 : KONTAC(ICAS) = 1 1597 : ELSE 1598 : CONTINUE 1599 : ENDIF 1600 : ENDIF 1601 : ENDIF 1602 : C 1603 : IF (ICASB.LE.4) THEN 1604 : IF (ICASB.LE.2) THEN 1605 : FG(1) = XKISO(1)*RESIST(5) 1606 : ELSE 1607 : FG(1) = XKISO(2)*RESIST(6) 1608 : ENDIF 1608 : CALL CRITPE (SIGTP(IHARDB),XZER,FG(1),ISGNTCB,XFFFSB) 1609 : IF (((XFFFSB.GE.EPSIL(5)).AND.(ICASB.LE.2)).OR. 1610 : & ((XFFFSB.GE.EPSIL(6)).AND.(ICASB.GE.3))) THEN 1611 : NCON = NCON +1 1612 : ELSE 1613 : IF (((ICASB.LE.2).AND.(ABS(XFFFSB).LT.EPSIL(5))).OR. 1614 : &((ICASB.GT.3).AND.(ABS(XFFFSB).LT.EPSIL(6)))) THEN 1615 : KONTAC(ICASB) = 1 1616 : ELSE 1617 : CONTINUE 1618 : ENDIF 1619 : ENDIF 1620 : ENDIF 1621 : C 1622 : IF (NCON.EQ.1) THEN 1623 : IF (((XFFFS.GE.EPSIL(5)).AND.(ICAS.LE.2)).OR. 1624 : & ((XFFFS.GE.EPSIL(6)).AND.(ICAS.GE.3))) THEN 1625 : GOTO 4040 1626 : ELSE 1627 : ICASI = ICAS 1628 : ISGNTCI = ISGNTC 1629 : IHARDI = IHARD 1630 : ICAS = ICASB 1631 : ISGNTC = ISGNTCB 1632 : IHARD = IHARDB 1633 : ICASB = ICASI 1634 : ISGNTCB = ISGNTCI 1635 : IHARDB = IHARDI 1636 : GOTO 4040 1637 : ENDIF 1638 : ELSE 1639 : IF (NCON.EQ.2) THEN 1640 : GOTO 4040 1641 : ELSE 1642 : CONTINUE 1643 : ENDIF 1644 : CONTINUE 1645 : ENDIF 1646 : C 1647 : C Regles de suivi avant de sortir normalement 1648 : C Direction de l'ecoulement 1649 : IF (ICAS.GE.5) THEN 1650 : SIGHP(IHARD) = XKISO(ICAS-4)*RESIST(ICAS) 1651 : & -RESIST(2*ICAS-10+IHARD) 1652 : ENDIF 1653 : IF (ICASB.GE.5) THEN 1654 : SIGHP(IHARDB) = XKISO(ICASB-4)*RESIST(ICASB) 1655 : & -RESIST(2*ICASB-10+IHARDB) 1656 : ENDIF 1657 : C 1658 : C Mis a jour des indicateurs de contact 1659 : C 1660 : IF ((XKISO(1)*RESIST(5)-XKISO(2)*RESIST(6)).GT. 1661 : &(RESIST(1)-RESIST(3))) THEN 1662 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) KONTAC(ICAS+2)=0 1663 : IF ((ICAS.EQ.3).OR.(ICAS.EQ.4)) KONTAC(ICAS-2)=0 1664 : IF ((ICASB.EQ.1).OR.(ICASB.EQ.2)) KONTAC(ICASB+2)=0 1665 : IF ((ICASB.EQ.3).OR.(ICASB.EQ.4)) KONTAC(ICASB-2)=0 1666 : ELSE 1667 : KONTAC(1) = 1 1668 : KONTAC(2) = 1 1669 : KONTAC(3) = 1 1670 : KONTAC(4) = 1 1671 : ENDIF 1672 : C 1673 : GOTO 9999 1674 : C 1675 : C On va au contact 1676 : 4040 CONTINUE 1677 : SIGHP(1) = SIGHP0(1) 1678 : SIGHP(2) = SIGHP0(2) 1679 : SIGT(1) = SIGI(1) + DSIGI(1) 1680 : SIGT(2) = SIGI(2) + DSIGI(2) 1681 : SIGT(4) = SIGI(4) + DSIGI(4) 1682 : EPSI(1) = EPSI0(1) 1683 : EPSI(2) = EPSI0(2) 1684 : GOTO 5000 1685 : C================================================================ 1686 : C CAS 5 1687 : C Deux mecanismes avec 1688 : C Entree en contact de 2 mecanismes 1689 : C================================================================ 1690 : 5000 CONTINUE 1691 : C 1691 : CALL CONPRI ( SIGT , CPHI0,SPHI0,SIGTP,CPHI,SPHI) 1692 : C 1693 : C ISGN2 et ISGN2B 1694 : IF (IHARD.EQ.1) THEN 1695 : IF (SIGTP(1).GE.SIGTP(2)) THEN 1696 : ISGN2 = 1 1697 : ELSE 1698 : ISGN2 = -1 1699 : ENDIF 1700 : ELSE 1701 : IF (IHARD.EQ.2) THEN 1702 : IF (SIGTP(1).GE.SIGTP(2)) THEN 1703 : ISGN2 = -1 1704 : ELSE 1705 : ISGN2 = 1 1706 : ENDIF 1707 : ENDIF 1708 : ENDIF 1709 : C 1710 : IF (IHARDB.EQ.1) THEN 1711 : IF (SIGTP(1).GE.SIGTP(2)) THEN 1712 : ISGN2B = 1 1713 : ELSE 1714 : ISGN2B = -1 1715 : ENDIF 1716 : ELSE 1717 : IF (IHARDB.EQ.2) THEN 1718 : IF (SIGTP(1).GE.SIGTP(2)) THEN 1719 : ISGN2B = -1 1720 : ELSE 1721 : ISGN2B = 1 1722 : ENDIF 1723 : ENDIF 1724 : ENDIF 1725 : C 1726 : 5010 CONTINUE 1727 : C 1728 : C Definition du segment 1729 : C 1730 : ITM = 30 1731 : MN = 3 1732 : SEGINI,QUASIN 1733 : IT = -1 1733 : CALL ZERO(D,1,MN) 1734 : C XKISO et XHH2 1735 : IF (ICASB.LE.4) THEN 1736 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) THEN 1736 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT, 1737 : &XKISO(1),XHH2,KERRE) 1737 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC, 1738 : &XKISO(2),DUMM,KERRE) 1739 : ELSE 1739 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT, 1740 : &XKISO(1),DUMM,KERRE) 1740 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC, 1741 : &XKISO(2),XHH2,KERRE) 1742 : ENDIF 1743 : ELSE 1744 : IF (ICASB.EQ.5) THEN 1745 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) THEN 1745 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1746 : & ,XHH2,KERRE) 1747 : XHH2B=XHH2 1747 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1748 : & ,DUMM,KERRE) 1749 : ELSE 1749 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1750 : & ,XHH2B,KERRE) 1750 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1751 : & ,XHH2,KERRE) 1752 : ENDIF 1753 : ELSE 1754 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) THEN 1754 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1755 : & ,XHH2,KERRE) 1755 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1756 : & ,XHH2B,KERRE) 1757 : ELSE 1757 : CALL VALEUR (EPSI0(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1758 : & ,DUMM,KERRE) 1758 : CALL VALEUR (EPSI0(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1759 : & ,XHH2,KERRE) 1760 : XHH2B = XHH2 1761 : ENDIF 1762 : ENDIF 1763 : ENDIF 1764 : C 1765 : XKISO0(1)=XKISO(1) 1766 : EPSI(1)=EPSI0(1) 1767 : XKISO0(2)=XKISO(2) 1768 : EPSI(2)=EPSI0(2) 1769 : XHH2I = XHH2 1770 : XHH2IB = XHH2B 1771 : XKISOI(1)=XKISO(1) 1772 : XKISOI(2)=XKISO(2) 1773 : EPSII(1)=EPSI0(1) 1774 : EPSII(2)=EPSI0(2) 1775 : C 1776 : DO 5011,IE1=1,6 1777 : FFF(IE1) = XZER 1778 : 5011 CONTINUE 1779 : C 1779 : CALL CRITPE (SIGTP(IHARD),SIGHP0(IHARD),RESIST(ICAS), 1780 : & ISGNTC,FFF(ICAS)) 1781 : C 1782 : IF (ICASB.LE.4) THEN 1782 : CALL CRITPE (SIGTP(IHARDB),SIGHP0(IHARDB),RESIST(ICASB), 1783 : & ISGNTCB,FFF(ICASB)) 1784 : ELSE 1785 : FG(1) = XKISOI(ICASB -4)*RESIST(ICASB) 1785 : CALL CRITPE (SIGTP(IHARDB), XZER, FG(1),ISGNTCB,FFF(ICASB)) 1786 : ENDIF 1787 : C 1788 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) THEN 1789 : FG(1) = XKISO(1)*RESIST(5) 1789 : CALL CRITPE (SIGTP(IHARD), XZER, FG(1), ISGNTC,XFFFS) 1790 : ELSE 1791 : FG(1) = XKISO(2)*RESIST(6) 1791 : CALL CRITPE (SIGTP(IHARD),XZER,FG(1),ISGNTC,XFFFS) 1792 : ENDIF 1793 : C 1793 : CALL NORMBL(SIGT,RESIST(ICAS),ISGNTC,ISGN2,1,XNOR1) 1793 : CALL NORMBL(SIGT,RESIST(ICASB),ISGNTCB,ISGN2B,1,XNOR2) 1794 : C 1795 : C Jacobien 1796 : C 1797 : XJ00(1,1) = -XXCCYY(XNOR1, XNOR1, YOUN, XNU) - XHH1 1798 : XJ00(1,2) = -XXCCYY(XNOR1, XNOR2, YOUN, XNU) 1799 : XJ00(1,3) = -XXYY(XNOR1 , DSIGI) 1800 : XJ00(2,1) = -XXCCYY(XNOR2, XNOR1, YOUN, XNU) 1801 : XJ00(2,2) = -XXCCYY(XNOR2, XNOR2, YOUN, XNU) 1802 : XJ00(2,3) = -XXYY(XNOR2,DSIGI) 1803 : XJ00(3,1) = -XXCCYY(XNOR1,XNOR1,YOUN,XNU) 1804 : XJ00(3,2) = -XXCCYY(XNOR1,XNOR2,YOUN,XNU) 1805 : XJ00(3,3) = -XXYY(XNOR1,DSIGI) 1806 : C 1807 : IF ((ICAS.EQ.3).OR.(ICAS.EQ.4)) THEN 1808 : XJ00(3,1)=XJ00(3,1)+RESIST(6)*XHH2*DEGRAD*ISGNTC*RESIST(ICAS) 1809 : ENDIF 1810 : C 1811 : IF (ICASB.LE.4) THEN 1812 : XJ00(2,2) = XJ00(2,2) - XHH1 1813 : IF ((ICAS.EQ.3).OR.(ICAS.EQ.4)) THEN 1814 : XJ00(3,2)=XJ00(3,2)+RESIST(6)*XHH2 1815 : &*DEGRAD*ISGNTCB*RESIST(ICASB) 1816 : ENDIF 1817 : ELSE 1818 : XJ00(2,2) = XJ00(2,2) - RESIST(ICASB)*XHH2B*ISGNTCB 1819 : ENDIF 1820 : C 1821 : C 1822 : R(1) = FFF(ICAS) 1823 : R(2) = FFF(ICASB) 1824 : R(3) = XFFFS 1825 : C 1826 : SPLITT = XZER 1827 : GAMBDA(1) = XZER 1828 : GAMBDA(2) = XZER 1829 : GAMBDA(3) = XZER 1830 : GAMBDA(4) = XZER 1831 : GAMBDA(5) = XZER 1832 : GAMBDA(6) = XZER 1833 : C 1834 : C Inverse du jacobien 1835 : C 1836 : XH0(1,1) = XJ00(1,1) 1837 : XH0(1,2) = XJ00(1,2) 1838 : XH0(1,3) = XJ00(1,3) 1839 : XH0(2,1) = XJ00(2,1) 1840 : XH0(2,2) = XJ00(2,2) 1841 : XH0(2,3) = XJ00(2,3) 1842 : XH0(3,1) = XJ00(3,1) 1843 : XH0(3,2) = XJ00(3,2) 1844 : XH0(3,3) = XJ00(3,3) 1844 : CALL INVALM (XH0, MN , MN , IRD, 1.D-12) 1845 : IF (IRD.NE.0) THEN 1846 : SEGSUP,QUASIN 1847 : GOTO 6000 1848 : ENDIF 1849 : C 1850 : C On rentre dans les iterations internes 1851 : C 1852 : DO I=0,ITM 1853 : C 1854 : C Appel de Broyden 1855 : C 1855 : CALL BROYDE (QUASIN) 1856 : C 1857 : GAMBDA(ICAS) = D(1) 1858 : GAMBDA(ICASB) = D(2) 1859 : SPLITT = D(3) 1860 : C 1861 : C Calcul du phi et des termes de la matrice A et A-1 1862 : C 1863 : SITPRO(1) = SIGI(1) + (UN-SPLITT)*DSIGI(1) 1864 : & - UNDEMI*(GAMBDA(1) + GAMBDA(2) + GAMBDA(5))*CPI 1865 : & + UNDEMI*(GAMBDA(3) + GAMBDA(4) + GAMBDA(6))*CPI 1866 : SITPRO(2) = SIGI(2) + (UN-SPLITT)*DSIGI(2) 1867 : & - UNDEMI*(GAMBDA(1) + GAMBDA(2) + GAMBDA(5))*CPI 1868 : & + UNDEMI*(GAMBDA(3) + GAMBDA(4) + GAMBDA(6))*CPI 1869 : SITPRO(4) = SIGI(4) + (UN-SPLITT)*DSIGI(4) 1870 : C 1871 : C PHI et AINV 1872 : C 1873 : SIGHP(1) = SIGHP0(1) + XHH1*(GAMBDA(1)-GAMBDA(3)) 1874 : SIGHP(2) = SIGHP0(2) + XHH1*(GAMBDA(2)-GAMBDA(4)) 1875 : EPSI(1) = EPSI0(1) + GAMBDA(5) 1876 : WORK = ISGNTC*RESIST(ICAS)*GAMBDA(ICAS) 1877 : IF (ICASB.LE.4) THEN 1878 : WORK = WORK + ISGNTCB*RESIST(ICASB)*GAMBDA(ICASB) 1879 : ENDIF 1880 : EPSI(2)=EPSI0(2)+ GAMBDA(6) + DEGRAD*WORK 1881 : C 1882 : IF (ICASB.LE.4) THEN 1883 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) THEN 1883 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT, 1884 : &XKISO(1),XHH2,KERRE) 1884 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC, 1885 : &XKISO(2),DUMM,KERRE) 1886 : ELSE 1886 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT, 1887 : &XKISO(1),DUMM,KERRE) 1887 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC, 1888 : &XKISO(2),XHH2,KERRE) 1889 : ENDIF 1890 : ELSE 1891 : IF (ICASB.EQ.5) THEN 1892 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) THEN 1892 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1893 : & ,XHH2,KERRE) 1894 : XHH2B=XHH2 1894 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1895 : & ,DUMM,KERRE) 1896 : ELSE 1896 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1897 : & ,XHH2B,KERRE) 1897 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1898 : & ,XHH2,KERRE) 1899 : ENDIF 1900 : ELSE 1901 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) THEN 1901 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1902 : & ,XHH2,KERRE) 1902 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1903 : & ,XHH2B,KERRE) 1904 : ELSE 1904 : CALL VALEUR (EPSI(1),TRAC(IPT),NCURVT,XLCAT,XKISO(1) 1905 : & ,DUMM,KERRE) 1905 : CALL VALEUR (EPSI(2),TRAC(IPC),NCURVC,XLCAC,XKISO(2) 1906 : & ,XHH2,KERRE) 1907 : XHH2B = XHH2 1908 : ENDIF 1909 : ENDIF 1910 : ENDIF 1911 : C 1912 : C Calcul des phi 1913 : C 1914 : PHI1=ISGNTC*(-UNDEMI*(SITPRO(1) + SITPRO(2)) 1915 : & + SIGHP(IHARD) + RESIST(ICAS)) 1916 : C 1917 : IF (ICASB.LE.4) THEN 1918 : PHI2=ISGNTCB*(-UNDEMI*(SITPRO(1) + SITPRO(2)) 1919 : & + SIGHP(IHARDB)+ RESIST(ICASB)) 1920 : ELSE 1921 : PHI2=ISGNTCB*(-UNDEMI*(SITPRO(1)+SITPRO(2)) 1922 : & + (XKISOI(ICASB-4)+XHH2IB*(EPSI(ICASB-4)-EPSII(ICASB-4))) 1923 : & *RESIST(ICASB)) 1924 : ENDIF 1925 : C 1926 : C Test sur phi1 et phi2 1927 : C 1928 : IF ((ABS(PHI2)).GT.(ABS(EPSILO*RESIST(ICASB)))) THEN 1929 : SOMLAM = GAMBDA(ICAS) + GAMBDA(ICASB)*PHI1/PHI2 1930 : PHI = PHI1 1931 : ELSE 1932 : IF ((ABS(PHI1)).GT.(ABS(EPSILO*RESIST(ICAS)))) THEN 1933 : SOMLAM = GAMBDA(ICASB) + GAMBDA(ICAS)*PHI2/PHI1 1934 : PHI = PHI2 1935 : ELSE 1936 : SOMLAM = GAMBDA(ICASB) + GAMBDA(ICAS) 1937 : PHI = PHI1 1938 : ENDIF 1939 : ENDIF 1940 : C 1941 : AINV(1) = (DEUX*PHI+GGG*SOMLAM)/(DEUX*PHI+DEUX*GGG*SOMLAM) 1942 : AINV(2) = (GGG*SOMLAM)/(DEUX*PHI + DEUX*GGG*SOMLAM) 1943 : AINV(3) = PHI/(PHI + GGG*SOMLAM) 1944 : C 1945 : C Calcul des contraintes a l'aide de AINV 1946 : C 1947 : SIGT(1) = AINV(1)*SITPRO(1)+AINV(2)*SITPRO(2) 1948 : SIGT(2) = AINV(2)*SITPRO(1)+AINV(1)*SITPRO(2) 1949 : SIGT(4) = AINV(3)*SITPRO(4) 1950 : C 1951 : C Nouvelles contraintes principales 1952 : C 1952 : CALL CONPRI ( SIGT , CPHI0,SPHI0,SIGTP,CPHI,SPHI) 1953 : C 1953 : CALL CRITPE (SIGTP(IHARD),SIGHP(IHARD),RESIST(ICAS), 1954 : & ISGNTC,FFF(ICAS)) 1955 : R(1)= + FFF(ICAS) 1956 : C 1957 : IF (ICASB.LE.4) THEN 1957 : CALL CRITPE (SIGTP(IHARDB),SIGHP(IHARDB),RESIST(ICASB), 1958 : & ISGNTCB,FFF(ICASB)) 1959 : R(2)= + FFF(ICASB) 1960 : ELSE 1961 : FG(1) = (XKISOI(ICASB -4)+XHH2IB*(EPSI(ICASB-4) 1962 : & - EPSII(ICASB-4)))*RESIST(ICASB) 1962 : CALL CRITPE (SIGTP(IHARDB), XZER, FG(1),ISGNTCB,FFF(ICASB)) 1963 : R(2) = + FFF(ICASB) 1964 : & - ISGNTCB*RESIST(ICASB) 1965 : & * (XKISO(ICASB-4) - XKISOI(ICASB-4) 1966 : & -(EPSI(ICASB-4) - EPSII(ICASB-4))*XHH2IB) 1967 : IF (ABS(EPSI(ICASB-4) - EPSII(ICASB-4)).LT.EPSILO2) THEN 1968 : XHH2IB=XHH2B 1969 : ELSE 1970 : XHH2IB = (XKISO(ICASB-4) - XKISOI(ICASB-4))/ 1971 : & (EPSI(ICASB-4) - EPSII(ICASB-4)) 1972 : ENDIF 1973 : ENDIF 1974 : C 1975 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) THEN 1976 : FG(1) = XKISO(1)*RESIST(5) 1976 : CALL CRITPE (SIGTP(IHARD), XZER, FG(1), ISGNTC,XFFFS) 1977 : R(3) = XFFFS 1978 : ELSE 1979 : FG(1) = (XKISOI(2)+XHH2I*(EPSI(2)-EPSII(2)))*RESIST(6) 1979 : CALL CRITPE (SIGTP(IHARD),XZER,FG(1),ISGNTC,XFFFS) 1980 : R(3) = XFFFS + RESIST(6)* 1981 : & (XKISO(2) - XKISOI(2)-(EPSI(2) - EPSII(2))*XHH2I) 1982 : IF (ABS(EPSI(2) - EPSII(2)).LT.EPSILO2) THEN 1983 : XHH2I=XHH2 1984 : ELSE 1985 : XHH2I = (XKISO(2) - XKISOI(2))/ 1986 : & (EPSI(2) - EPSII(2)) 1987 : ENDIF 1988 : ENDIF 1989 : XKISOI(1) = XKISO(1) 1990 : EPSII(1) = EPSI(1) 1991 : XKISOI(2) = XKISO(2) 1992 : EPSII(2) = EPSI(2) 1993 : C 1994 : C Test de convergence 1995 : C 1996 : IF (ICAS.LE.2) THEN 1997 : IF ((ABS(FFF(ICAS)).LT.EPSIL(ICAS)).AND. 1998 : & (ABS(R(1)).LT.EPSIL(ICAS)).AND. 1999 : & (ABS(FFF(ICASB)).LT.EPSIL(ICASB)).AND. 2000 : & (ABS(R(2)).LT.EPSIL(ICASB)).AND. 2001 : & (ABS(XFFFS).LT.EPSIL(5)).AND. 2002 : & (ABS(R(3)).LT.EPSIL(5))) THEN 2003 : GOTO 5020 2004 : ENDIF 2005 : ELSE 2006 : IF ((ABS(FFF(ICAS)).LT.EPSIL(ICAS)).AND. 2007 : & (ABS(R(1)).LT.EPSIL(ICAS)).AND. 2008 : & (ABS(FFF(ICASB)).LT.EPSIL(ICASB)).AND. 2009 : & (ABS(R(2)).LT.EPSIL(ICASB)).AND. 2010 : & (ABS(XFFFS).LT.EPSIL(6)).AND. 2011 : & (ABS(R(3)).LT.EPSIL(6))) THEN 2012 : GOTO 5020 2013 : ENDIF 2014 : ENDIF 2015 : ENDDO 2016 : C 2017 : C Pas de convergence 2018 : C 2019 : SEGSUP,QUASIN 2020 : GOTO 6000 2021 : C 2022 : 5020 CONTINUE 2023 : SEGSUP,QUASIN 2024 : C 2025 : C 2026 : IF (((SPLITT-1.).LE.EPSILO).AND.(SPLITT.GE.(-1.*EPSILO)).AND. 2027 : &(GAMBDA(ICAS).GE.EPSILO2) 2028 : &.AND.(GAMBDA(ICASB).GE.EPSILO2)) THEN 2029 : C 2030 : C On a les contraintes, les variables d'ecrouissage cinemat. et tgphi 2031 : C Mis a jour des variables d'ecrouissage isotrope ou regles de suivi 2032 : C 2033 : 2034 : IF (NCON.EQ.2) THEN 2035 : IF (ICASB.LE.4) THEN 2036 : C 2037 : IF (ICASB.LE.2) THEN 2038 : FG(1) = XKISO(1)*RESIST(5) 2039 : ELSE 2040 : FG(1) = XKISO(2)*RESIST(6) 2041 : ENDIF 2041 : CALL CRITPE (SIGTP(IHARDB),XZER,FG(1),ISGNTCB,XFFFSB) 2042 : C 2043 : C On regarde si on ne s'est pas trompe d'ordre 2044 : C 2045 : IF (((XFFFSB.GE.EPSIL(5)).AND.(ICASB.LE.2)).OR. 2046 : & ((XFFFSB.GE.EPSIL(6)).AND.(ICASB.GE.3))) THEN 2047 : ICASI = ICAS 2048 : ISGNTCI = ISGNTC 2049 : IHARDI = IHARD 2050 : ICAS = ICASB 2051 : ISGNTC = ISGNTCB 2052 : IHARD = IHARDB 2053 : ICASB = ICASI 2054 : ISGNTCB = ISGNTCI 2055 : IHARDB = IHARDI 2056 : GOTO 4040 2057 : ELSE 2058 : IF ((ICASB.LE.2).AND.(ABS(XFFFSB).LT.EPSIL(5))) THEN 2059 : KONTAC(ICASB) = 1 2060 : ICASB=5 2061 : ENDIF 2062 : IF (((ICASB.EQ.3).OR.(ICASB.EQ.4)) 2063 : & .AND.(ABS(XFFFSB).LT.EPSIL(6))) THEN 2064 : KONTAC(ICASB) = 1 2065 : ICASB=6 2066 : ENDIF 2067 : ENDIF 2068 : C 2069 : ELSE 2070 : CONTINUE 2071 : ENDIF 2072 : ELSE 2073 : CONTINUE 2074 : ENDIF 2075 : C 2076 : C Le contact est realise 2077 : C 2078 : KONTAC(ICAS) = 1 2079 : IF ((XKISO(1)*RESIST(5)-XKISO(2)*RESIST(6)).GT. 2080 : &(RESIST(1)-RESIST(3))) THEN 2081 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) KONTAC(ICAS+2)=0 2082 : IF ((ICAS.EQ.3).OR.(ICAS.EQ.4)) KONTAC(ICAS-2)=0 2083 : IF ((ICASB.EQ.1).OR.(ICASB.EQ.2)) KONTAC(ICASB+2)=0 2084 : IF ((ICASB.EQ.3).OR.(ICASB.EQ.4)) KONTAC(ICASB-2)=0 2085 : ELSE 2086 : KONTAC(1) = 1 2087 : KONTAC(2) = 1 2088 : KONTAC(3) = 1 2089 : KONTAC(4) = 1 2090 : ENDIF 2091 : C 2092 : SIGI(1) = SIGT(1) 2093 : SIGI(2) = SIGT(2) 2094 : SIGI(4) = SIGT(4) 2095 : DSIGI(1) = DSIGI(1)*SPLITT 2096 : DSIGI(2) = DSIGI(2)*SPLITT 2097 : DSIGI(4) = DSIGI(4)*SPLITT 2098 : SIGT(1) = SIGI(1) + DSIGI(1) 2099 : SIGT(2) = SIGI(2) + DSIGI(2) 2100 : SIGT(4) = SIGI(4) + DSIGI(4) 2101 : SIGHP0(1) = SIGHP(1) 2102 : SIGHP0(2) = SIGHP(2) 2103 : EPSI0(1) = EPSI(1) 2104 : EPSI0(2) = EPSI(2) 2105 : NPASS2 = 0 2106 : CPHI0 = CPHI 2107 : SPHI0 = SPHI 2108 : IF ((ICAS.EQ.1).OR.(ICAS.EQ.2)) THEN 2109 : ICAS = 5 2110 : ELSE 2111 : ICAS = 6 2112 : ENDIF 2113 : IF (ICASB.LE.4) THEN 2114 : GOTO 4000 2115 : ELSE 2116 : IF (ICASB.EQ.ICAS) THEN 2117 : C On passe de 2 mecanismes a 1 mecanisme 2118 : GOTO 1000 2119 : ELSE 2120 : GOTO 4000 2121 : ENDIF 2122 : ENDIF 2123 : C 2124 : C 2125 : ELSE 2126 : 5040 CONTINUE 2127 : IF ((GAMBDA(ICASB).LT.EPSILO2).AND.(NPASS.LE.3)) THEN 2128 : SIGHP(1) = SIGHP0(1) 2129 : SIGHP(2) = SIGHP0(2) 2130 : SIGT(1) = SIGI(1) + DSIGI(1) 2131 : SIGT(2) = SIGI(2) + DSIGI(2) 2132 : SIGT(4) = SIGI(4) + DSIGI(4) 2133 : IF (NPASS.LT.5) THEN 2134 : NPASS = NPASS + 1 2135 : GOTO 3000 2136 : ELSE 2137 : GOTO 6000 2138 : ENDIF 2139 : ELSE 2140 : IF ((GAMBDA(ICAS).LT.EPSILO2).AND.(NPASS2.EQ.0)) THEN 2141 : NPASS2 = 1 2142 : SIGHP(1) = SIGHP0(1) 2143 : SIGHP(2) = SIGHP0(2) 2144 : SIGT(1) = SIGI(1) + DSIGI(1) 2145 : SIGT(2) = SIGI(2) + DSIGI(2) 2146 : SIGT(4) = SIGI(4) + DSIGI(4) 2147 : ICAS = ICASB 2148 : IHARD = IHARDB 2149 : ISGNTC = ISGNTCB 2150 : IF (NPASS.LT.5) THEN 2151 : NPASS = NPASS + 1 2152 : GOTO 1000 2153 : ELSE 2154 : GOTO 6000 2155 : ENDIF 2156 : ELSE 2157 : GOTO 6000 2158 : ENDIF 2159 : ENDIF 2160 : ENDIF 2161 : C 2162 : C On bissecte 2163 : C 2164 : 6000 CONTINUE 2165 : NDICHO = NDICHO + 1 2166 : IF (NDICHO.LE.90) THEN 2167 : IDICHO = 1 2168 : DSIGIDI(1) = DSIGIDI(1) + UNDEMI*DSIGI(1) 2169 : DSIGIDI(2) = DSIGIDI(2) + UNDEMI*DSIGI(2) 2170 : DSIGIDI(4) = DSIGIDI(4) + UNDEMI*DSIGI(4) 2171 : DSIGI(1) = DSIGI(1)*UNDEMI 2172 : DSIGI(2) = DSIGI(2)*UNDEMI 2173 : DSIGI(4) = DSIGI(4)*UNDEMI 2174 : CPHI0 = CPHI00 2175 : SPHI0 = SPHI00 2176 : GOTO 10 2177 : ELSE 2178 : KERRE = 2 2179 : RETURN 2180 : ENDIF 2181 : C 2182 : C================================================================= 2183 : C MIS A JOUR DE TOUTES LES VARIABLES ET CONTRAINTES AVANT 2184 : C DE SORTIR 2185 : C================================================================= 2186 : 9999 CONTINUE 2187 : C Cas de la bissection 2188 : IF (IDICHO.EQ.1) THEN 2189 : IDICHO = 0 2190 : NDICHO = 0 2191 : NPASS = 0 2192 : DSIGI(1) = DSIGIDI(1) 2193 : DSIGI(2) = DSIGIDI(2) 2194 : DSIGI(4) = DSIGIDI(4) 2195 : SIGI(1) = SIGT(1) 2196 : SIGI(2) = SIGT(2) 2197 : SIGI(4) = SIGT(4) 2198 : DSIGIDI(1) = XZER 2199 : DSIGIDI(2) = XZER 2200 : DSIGIDI(4) = XZER 2201 : SIGHP0(1) = SIGHP(1) 2202 : SIGHP0(2) = SIGHP(2) 2203 : EPSI0(1) = EPSI(1) 2204 : EPSI0(2) = EPSI(2) 2205 : CPHI0 = CPHI 2206 : CPHI00 = CPHI0 2207 : SPHI0 = SPHI 2208 : SPHI00 = SPHI0 2209 : KONTA0(1) = KONTAC(1) 2210 : KONTA0(2) = KONTAC(2) 2211 : KONTA0(3) = KONTAC(3) 2212 : KONTA0(4) = KONTAC(4) 2213 : GOTO 10 2214 : ENDIF 2215 : C STRESS 2216 : SIGF(1) = SIGT(1) 2217 : SIGF(2) = SIGT(2) 2218 : SIGF(4) = SIGT(4) 2219 : SIGF(3) = XZER 2220 : C PLASTIC STRAIN 2221 : CINV(1) = UN/YOUN 2222 : CINV(2) = -XNU/YOUN 2223 : CINV(3) = UN/GGG 2224 : DEFP(1) =CINV(1)*(SIG0(1) + DSIGT(1) - SIGT(1)) 2225 : & +CINV(2)*(SIG0(2) + DSIGT(2) - SIGT(2)) 2226 : DEFP(2) =CINV(2)*(SIG0(1) + DSIGT(1) - SIGT(1)) 2227 : & +CINV(1)*(SIG0(2) + DSIGT(2) - SIGT(2)) 2228 : DEFP(4) =CINV(3)*(SIG0(4) + DSIGT(4) - SIGT(4)) 2229 : DEFP(3) = XZER 2230 : C 2231 : C Cas ou la grande surface est a l'interieur de la petite 2232 : C 2233 : IF ((XKISO(1)*RESIST(5)-XKISO(2)*RESIST(6)).LE. 2234 : &(RESIST(1)-RESIST(3))) THEN 2235 : KONTAC(1) = 1 2236 : KONTAC(2) = 1 2237 : KONTAC(3) = 1 2238 : KONTAC(4) = 1 2239 : ENDIF 2240 : VARF(1) = SIGHP(1) 2241 : VARF(2) = SIGHP(2) 2242 : VARF(3) = EPSI(1) 2243 : VARF(4) = EPSI(2) 2244 : VARF(5) = CPHI - 1.D0 2245 : VARF(6) = SPHI 2246 : VARF(7) = KONTAC(1) 2247 : VARF(8) = KONTAC(2) 2248 : VARF(9) = KONTAC(3) 2249 : VARF(10) = KONTAC(4) 2250 : VARF(11) = ICAS 2251 : VARF(12) = ICASB 2252 : VARF(13) = GAMBDA(ICAS) 2253 : VARF(14) = GAMBDA(ICASB) 2254 : C===================================================== 2255 : C FIN DE LA ROUTINE DU MODELE 2256 : C====================================================== 2257 : RETURN 2258 : END 2259 : 2260 : 2261 :
© Cast3M 2003 - All rights reserved.
Disclaimer