1 : C COUL2 SOURCE AM 08/08/29 21:15:05 6145 2 : SUBROUTINE COUL2(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI, 3 : . DEPST,IFOUR,XMAT,NMATT,IVAL,DD,SIGF,DEFP,VARF,KERRE) 4 : C-------------------------------------------------------------- 5 : c 6 : C 7 : C PLASTICITE MODELE 2D MOHR COULOMB ELEMENTS JOINTS 8 : C 9 : C ENTREES 10 : C IB = NUMERO DE L'ELEMENT 11 : C IGAUSS = NUMERO DU POINT DE GAUSS 12 : C NSTRS = NOMBRE DE COMPOSANTES DE CONTRAINTES 13 : C SIG0 = CONTRAINTES INITIALES (AU PAS PRECEDENT) 14 : C EPIN0 = DEFORMATIONS INITIALES (AU PAS PRECEDENT) 15 : C VAR0(NVARI) = VARIABLES INTERNES DEBUT 16 : C VAR0(1) = DEFORMATION PLASTIQUE EQUIVALENTE 17 : C VAR0(2) = EPOUN (DEFO. PLASTIQUE DUE A L'OUVERTURE SEULE) 18 : C VAR0(3) = STAT (ETAT DU JOINT) 19 : C VAR0(4) = LAM1 ( MULTIPLICATEUR PLASTIQUE EN TRACTION ) 20 : C DSIGT(NSTRS) = INCREMENT DE CONTRAINTES TOTALES 21 : C XMAT(NCOMAT) = COMPOSANTES DE MATERIAU 22 : C NMATT = NOMBRE DE COMPOSANTES DE MATERIAU 23 : C IVAL(NCOMAT) = INDICE DES COMPOSANTES DE MATERIAU 24 : C DD(NSTRS,NSTRS) = MATRICE DE HOOK 25 : C 26 : C SORTIES 27 : C SIGF(NSTRS) = CONTRAINTES FINALES 28 : C DEFP(NSTRS) = INCREMENTS DE DEFORMATIONS INELASTIQUES FINALES 29 : C VARF(NVARI) = VARIABLES INTERNES FINALES 30 : C KERRE = INDICE D'ERREUR 31 : C = 60 PAS DE CONVERGENCE APRES 10 ITERATIONS 32 : C INTERNES (VARIABLE ISSU DANS COULM) 33 : C = 61 34 : C 35 : C----------------------------------------------------------------------- 36 : C 37 : IMPLICIT INTEGER(I-N) 38 : IMPLICIT REAL*8(A-H,O-Z) 39 : C Include contenant quelques constantes dont XPI : 40 : -INC CCREEL 41 : C 42 : DIMENSION SIG0(*),EPIN0(*),VAR0(*),DEPST(*),XMAT(*) 43 : DIMENSION DD(NSTRS,NSTRS) 44 : DIMENSION SIGF(*),DEFP(*),VARF(*) 45 : KERRE = 0 46 : C 47 : C 48 : C 49 : C ----------PARAMETRES POUR LE TEST DE NULLITE DES CONTRAINTES ET 50 : C ----------DEFORMATIONS 51 : C 52 : PRESP = +1.D-5 53 : PRESM = -1.D-5 54 : ccc PREEP = +1.D-10 55 : ccc PREEM = -1.D-10 56 : PREEP = +1.D-4 57 : PREEM = -1.D-7 58 : c 59 : C 60 : C ----------PARAMETRE D'IMPPRESSION 61 : C 62 : IMPRI = 0 63 : C 64 : C ----------NOMBRE MAXIMUM D'ITERATIONS INTERNES AUTORISEES 65 : C 66 : MAXITE = 1000 67 : C 68 : C 69 : C------CHARGEMENT DU EPSILON DE PRECISION DES ITERATIONS INTERNES 70 : C 71 : PRECIS = 1.D-6 72 : PRECIN = -1.D-6 73 : C 74 : C 75 : C------CHARGEMENT DE LA MATRICE DE HOOK INVERSE 76 : C------POUR LE CALCUL DE L'INCREMENT DE DEFORMATION 77 : C------A PARTIR DE L'INCREMENT DE CONTRAINTE 78 : C 79 : DD(1,1) = 1.0D0/XMAT(1) 80 : DD(1,2) = 0.0D0 81 : DD(2,1) = 0.0D0 82 : DD(2,2) = 1.0D0/XMAT(2) 83 : c 84 : c contrainte en fin de pente ec 85 : c 86 : sigsei = -xmat(6)*xmat(2) 87 : c 88 : C 89 : C-----------INCREMENT DE DEFORMATION 90 : C 91 : DEPS1 = DEPST(1) 92 : DEPSN = DEPST(2) 93 : C 94 : C TRAITEMENT DU PASSAGE DIRECT ENTRE UN ETAT OUVERT 95 : C (I.E. EXISTENCE D'UN JEU AU TEMPS N) A UN ETAT 96 : C COMPRIME AU TEMPS N+1 97 : C 98 : DEFPN=0.0D0 99 : INDIC=0 100 : IFIN =0 101 : STOTAN = 0.0D0 102 : DEFP11 = 0.0D0 103 : DEFP12 = 0.0D0 104 : ALP = 1.0D0 105 : S1INT = 0.0D0 106 : ZJEU = EPIN0(2) 107 : C CAS DU RAPPROCHEMENT DES LEVRES 108 : C 109 : C AM CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION 110 : C 111 : C IF (DEPSN.LT.-1.D-4) THEN 112 : IF (DEPSN.LT.-1.D-4.AND.SIG0(2).GT.PRESP) THEN 113 : C IL EXISTE UN JEU 114 : IF (ZJEU.GT.1.D-3) THEN 115 : C PASSAGE DIRECT D'UN ETAT OUVERT 116 : C A UN ETAT COMPRIME 117 : valdif=(ABS(DEPSN) - ZJEU)/ABS(DEPSN) 118 : IF (valdif.gt.1.D-3) THEN 119 : C 120 : ALP = ZJEU/(ABS(DEPSN)) 121 : STOTAN= DEPS1 122 : DEPS1 = ALP * DEPS1 123 : STOCK = DEPSN 124 : DEPSN = ALP * DEPSN 125 : INDIC = 1 126 : GO TO 99 127 : ENDIF 128 : C 129 : ENDIF 130 : C 131 : ENDIF 132 : C 133 : 98 CONTINUE 134 : C 135 : C IL RESTE A PARCOURIR LE TRAJET 136 : C (1-ALP) * DEPSN 137 : C 138 : IF (INDIC.EQ.1) THEN 139 : DEPS1=(1.0D0-ALP)*STOTAN 140 : DEPSN=(1.0D0-ALP)*STOCK 141 : INDIC=0 142 : IFIN =1 143 : ENDIF 144 : C 145 : 99 CONTINUE 146 : C 147 : C 148 : C ----CAS GENERAL SANS TRANSITION DIRECTE D'UN ETAT OUVERT A UN ETAT 149 : C ----A UN ETAT COMPRIME 150 : C 151 : C 152 : C 153 : C 154 : IF ((INDIC.EQ.1).OR.(IFIN.EQ.0)) THEN 155 : C 156 : C--------------------CALCUL DE L ETAT DE CONTRAINTES, DE DEFORMATIONS 157 : C--------------------DES VARIABLES INTERNES AU PAS PRECEDENT 158 : C 159 : C------CHARGEMENT DES VALEURS AU PAS PRECEDENT 160 : C 161 : C 162 : C----------DEFORMATIONS 163 : C 164 : C DEFORMATIONS INELASTIQUES 165 : C 166 : EPSS1 = EPIN0(1) 167 : EPSN1 = EPIN0(2) 168 : C 169 : C DEFORMATIONS ELASTIQUES 170 : C 171 : DEFS10 = DD(1,1)*SIG0(1) 172 : c 173 : if (sig0(2).lt.sigsei) then 174 : DEFN0 = ( -xmat(6) ) + ( (sig0(2) - sigsei) / xmat(5) ) 175 : else if ( (sig0(2).ge.sigsei) .and. (sig0(2).le.presp) ) then 176 : DEFN0 = DD(2,2)*SIG0(2) 177 : else if (sig0(2).gt.presp) then 178 : C 179 : C AM CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION 180 : C DEFN0 = 0.0D0 181 : DEFN0 = DD(2,2)*SIG0(2) 182 : endif 183 : C 184 : C DEFORMATIONS TOTALES AU PAS PRECEDENT N 185 : C 186 : EPS1N = EPSS1 + DEFS10 187 : EPSNN = EPSN1 + DEFN0 188 : C 189 : ENDIF 190 : C 191 : C 192 : C 193 : C 194 : C ----CAS DU 2ND PASSAGE : D'UN ETAT FERME ON PASSE A UN ETAT COMPRIME 195 : C 196 : C 197 : C 198 : C 199 : IF (IFIN.EQ.1) THEN 200 : C 201 : C--------------------CALCUL DE L ETAT DE CONTRAINTES, DE DEFORMATIONS 202 : C--------------------DES VARIABLES INTERNES AU PAS PRECEDENT 203 : C 204 : C------CHARGEMENT DES VALEURS AU PAS PRECEDENT 205 : C 206 : C 207 : C----------DEFORMATIONS 208 : C 209 : C DEFORMATIONS INELASTIQUES 210 : C 211 : EPSS1 = DEFP11 212 : EPSN1 = 0.0D0 213 : C 214 : C DEFORMATIONS ELASTIQUES 215 : C 216 : DEFS10 = DD(1,1)*S1INT 217 : DEFN0 = 0.0D0 218 : C 219 : C DEFORMATIONS TOTALES AU PAS PRECEDENT N 220 : C 221 : EPS1N = EPSS1 + DEFS10 222 : EPSNN = EPSN1 + DEFN0 223 : C 224 : ENDIF 225 : C 226 : C ESTIMATION DE LA DEFO. TOTALE AU PAS N+1 227 : C 228 : EP1NP1 = EPS1N + DEPS1 229 : EPNNP1 = EPSNN + DEPSN 230 : C 231 : C----------EPSILON ELASTIQUE (1 2 N) AU PAS N+1 A L'ITERATION 0 232 : C 233 : EE1P10 = DEFS10 + DEPS1 234 : EENP10 = DEFN0 + DEPSN 235 : C 236 : C----------SIGMA (1 2 N) AU PAS N+1 A L'ITERATION 0 237 : C 238 : S1NP10 = XMAT(1) * EE1P10 239 : c 240 : x6 = -1.0d0 * xmat(6) 241 : c 242 : c AM CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION 243 : c 244 : c if (depsn.lt.preem) then 245 : if (depsn.lt.preem.and.sig0(2).le.presp) then 246 : c 247 : c -----------fermeture du joint 248 : c 249 : if (epnnp1.le.preep) then 250 : c 251 : c --------------fermeture du joint avec compression des levres 252 : c 253 : if ( (eenp10.ge.x6) .and. (eenp10.le.preep) ) then 254 : c 255 : c ----------------compression : pente EC 256 : c 257 : SNNP10 = XMAT(2) * EENP10 258 : endif 259 : c 260 : if (eenp10.lt.x6) then 261 : c 262 : c ----------------compression : pente Ef 263 : c 264 : snnp10= ( xmat(2) * x6 ) + 265 : . ( (eenp10 + xmat(6)) * xmat(5) ) 266 : endif 267 : c 268 : endif 269 : c 270 : if ((epnnp1.gt.preem).or.(indic.eq.1)) then 271 : c 272 : c --------------rapprochement des levres sans compression 273 : c --------------ou fermeture des levres sans compression 274 : c 275 : SNNP10 = 0.0D0 276 : c 277 : endif 278 : c 279 : 280 : c 281 : else 282 : c 283 : c -----------ouverture du joint ou fermeture si joint pas casse 284 : c 285 : SNNP10 = XMAT(2) * EENP10 286 : c 287 : endif 288 : C 289 : C-----------VARIABLES INTERNES 290 : C 291 : C DEFORMATION PLASTIQUE EQUIVALENTE 292 : C 293 : EPSEQ=VAR0(1) 294 : C 295 : C DEFORMATION PLASTIQUE DUE A L'OUVERTURE SEULE 296 : C 297 : EPOUN=VAR0(2) 298 : C 299 : C ETAT DU JOINT 300 : C 301 : C STAT = 0 LEVRES EN CONTACT, AVEC OU SANS CISAILLEMENT 302 : C STAT = 1 JOINT EN OUVERTURE PENTE EC 303 : C STAT = 2 JOINT EN FERMETURE PENTE EC 304 : C STAT = 3 JOINT EN FERMETURE PENTE EF 305 : C 306 : C 307 : C 308 : STAT =VAR0(3) 309 : C 310 : C MULTIPLICATEUR PLASTIQUE EN TRACTION 311 : C 312 : XLAM1=VAR0(4) 313 : C 314 : C 315 : C---- MOHR-COULOMB NORMAL 316 : C 317 : C = XMAT(7) 318 : PHI = XMAT(8) 319 : FTRC = XMAT(4) 320 : C 321 : C AM 28/7/95 322 : C FTRC EST EN 4 ET NON EN 9 A CAUSE DE RHO ET ALPH 323 : C VOIR COMMENTAIRE DANS ECOUL2 324 : C 325 : C LA LIMITE EN TRACTION EST MISE BRUTALEMENT A 0 326 : C DES QU'ELLE EST DEPASSEE. 327 : C 328 : C 329 : PHI = (PHI*XPI)/180.D0 330 : ccc ZMU = XMAT(8) 331 : ZMU = 0.0D0 332 : ZMU = (ZMU*XPI)/180.D0 333 : C 334 : C MISE A 0 DE LA LIMITE EN TRACTION SI DEJA CASSE OU SI JEU 335 : C 336 : RTRAC=FTRC 337 : IF(XLAM1.GT.0.D0) RTRAC=0.D0 338 : C 339 : C-----------ETUDE DE L'ECOULEMENT SELON LES DIFFERENTS CRITERES 340 : C 341 : S2 = 0.0D0 342 : CRICIS = CALG(S1NP10,S2,SNNP10,PHI,C) 343 : CRINOR = SNNP10 - RTRAC 344 : CRICIP = CALN(S1NP10,S2,SNNP10,PHI,C) 345 : ABSTAU = SQRT( (S1NP10 ** 2) + (S2 ** 2) ) 346 : C 347 : C----------AIGUILLAGE SUIVANT LA ZONE 0, 1, 2, 3 348 : C 349 : IF ((CRICIS.LE.PRECIS).AND.(CRINOR.LE.PRECIS)) THEN 350 : GO TO 10 351 : ENDIF 352 : C 353 : IF ((CRICIS.GT.PRECIS).AND.(CRICIP.LE.PRECIS)) THEN 354 : GO TO 40 355 : ENDIF 356 : C 357 : IF ((CRICIP.GT.PRECIS).AND.(ABSTAU.GT.C)) THEN 358 : GO TO 20 359 : ENDIF 360 : C 361 : IF ((ABSTAU.LE.C).AND.(CRINOR.GT.PRECIN)) THEN 362 : GO TO 30 363 : ENDIF 364 : C 365 : C 366 : C-----------CAS DE LA CONVERGENCE A L ITERATION 0 : 367 : C 368 : 10 CONTINUE 369 : C 370 : C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS, 371 : C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE 372 : C 373 : C 374 : IF (IMPRI.EQ.1) THEN 375 : WRITE(6,*) 'PASSAGE N. 1' 376 : WRITE(6,*) 'LES DEUX CRITERES SONT VERIFIES-PAS D ECOULEMENT' 377 : WRITE(6,*) 'CONVERGENCE A L ITERATION 0' 378 : ENDIF 379 : C 380 : C 381 : C 382 : IF (INDIC.EQ.1) THEN 383 : C FERMETURE DU JOINT SANS COMPRESSION 384 : S1INT = S1NP10 385 : DEFP11= 0.0D0 386 : cc DEFPN = 0.0D0 387 : DEFPN = -epin0(2) 388 : GO TO 98 389 : C 390 : ELSE 391 : C 392 : IF (IFIN.EQ.0) THEN 393 : C 394 : DEFP(1) = 0.0D0 395 : C 396 : if ( abs(snnp10).le.presp) then 397 : c 398 : c deformation plastique du pas (n+1) moins 399 : c deformation plastique du pas (n) 400 : c 401 : defp(2) = epnnp1 - epsn1 402 : else 403 : c 404 : c le joint se deforme toujours elastiquement en fermeture 405 : c 406 : DEFP(2) = 0.0D0 407 : endif 408 : c 409 : ELSE 410 : c 411 : c cas du passage direct d'un etat ouvert a un etat 412 : c comprime 413 : c 414 : DEFP12 = 0.0D0 415 : DEFP(1) = 0.0d0 416 : cc DEFP(2) = 0.0d0 417 : DEFP(2) = defpn 418 : IFIN = 0 419 : ENDIF 420 : C 421 : C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN 422 : C 423 : c 424 : c calcul du delta epsilon plastique entre 425 : c le pas actuel et le pas precedent 426 : c 427 : delep1 = defp(1) 428 : delepn = defp(2) 429 : c 430 : c produit contracte pour le calcul de la defo. plasti. cumulee 431 : c 432 : pro1 = ( delep1 * delep1 ) + ( delepn * delepn ) 433 : depeq = sqrt( (2.0D0/3.0D0) * pro1 ) 434 : VARF(1) = EPSEQ + DEPEQ 435 : VARF(4) = XLAM1 436 : c 437 : c calcul de la deformation normale due a l'ouverture seule 438 : c 439 : if ( abs(snnp10).le.presp) then 440 : epoun = defp(2) 441 : VARF(2) = EPOUN 442 : else 443 : epoun = 0.0d0 444 : VARF(2) = EPOUN 445 : endif 446 : c 447 : c etat du joint stat = 1 cisaillement pur, levres en contact 448 : c ou ouverture du joint pente EC 449 : c stat = 2 compression du joint pente EC 450 : c stat = 3 compression du joint pente EF 451 : c 452 : if ( abs(snnp10).le.presp) then 453 : stat = 1.D0 454 : if ( abs(epnnp1).le.1.D-8 ) then 455 : stat = 0.D0 456 : endif 457 : VARF(3) = STAT 458 : endif 459 : if ( (snnp10.gt.sigsei).and.(snnp10.lt.presm) ) then 460 : stat = 2.D0 461 : VARF(3) = STAT 462 : endif 463 : if (snnp10.le.sigsei) then 464 : stat = 3.D0 465 : VARF(3) = STAT 466 : endif 467 : C 468 : SIGF(1) = S1NP10 469 : SIGF(2) = SNNP10 470 : C 471 : RETURN 472 : c 473 : ENDIF 474 : C 475 : C 476 : C----------CRITERES DE L'EFFORT NORMAL ET DE CISAILLEMENT 477 : C----------TOUS LES DEUX VIOLES EN MEME TEMPS 478 : C 479 : 20 CONTINUE 480 : C 481 : C ECOULEMENT SELON LES DEUX CRITERES A L'ITERATION ITER 482 : C 483 : ITER = 1 484 : DELAM1=0.D0 485 : C 486 : 2 CONTINUE 487 : C 488 : IF (ITER.GT.MAXITE) THEN 489 : WRITE(6,*) ' LES DEUX CRITERES SONT VIOLES ' 490 : WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES' 491 : KERRE = 2 492 : RETURN 493 : ENDIF 494 : C 495 : IF (IMPRI.EQ.1) THEN 496 : write(6,*) 'PASSAGE N. 2' 497 : write(6,*) 'LES DEUX CRITERES SONT VIOLES' 498 : write(6,*) 'ITERATION',ITER 499 : ENDIF 500 : C 501 : C 502 : C CALCUL DES CORRECTEURS PLASTIQUES 503 : C 504 : C 505 : C DERIVEE DU CRITERE F SUIVANT SIGMA 506 : C 507 : DFIDS1 = 0.0D0 508 : DFIDSN = 1.0D0 509 : C 510 : C DERIVEE DU CRITERE G SUIVANT SIGMA 511 : C 512 : if ( s1np10.lt.precis ) then 513 : dgids1 = -1.d0 514 : else 515 : dgids1 = 1.d0 516 : endif 517 : c 518 : DZGIDS = TAN(ZMU) 519 : DGIDSN = TAN(PHI) 520 : C 521 : C MATRICE DE HOOK * DFI/DSIGMA 522 : C 523 : HODF1 = XMAT(1) * DFIDS1 524 : HODFN = XMAT(2) * DFIDSN 525 : C 526 : C MATRICE DE HOOK * DGI/DSIGMA 527 : C 528 : HODG1 = XMAT(1) * DGIDS1 529 : HODGN = XMAT(2) * DZGIDS 530 : C 531 : C FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA 532 : C 533 : FIFI = (DFIDS1 * HODF1) + (DFIDSN * HODFN) 534 : C 535 : C FIGI = DFI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA 536 : C 537 : FIGI = (DFIDS1 * HODG1) + (DFIDSN * HODGN) 538 : C 539 : C GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA 540 : C 541 : GIFI = (DGIDS1 * HODF1) + (DGIDSN * HODFN) 542 : C 543 : C GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA 544 : C 545 : GIGI = (DGIDS1 * HODG1) + (DGIDSN * HODGN) 546 : C 547 : C RESOLUTION DU SYSTEME POUR TROUVER LES CORRECTEURS 548 : C PLASTIQUES DLAM1 ET DLAM2 549 : C 550 : DET = (FIFI * GIGI) - (FIGI * GIFI) 551 : C 552 : IF (ABS(DET).LE.PRECIS) THEN 553 : WRITE(6,*) ' SYSTEME INDETERMINE (F G): ITERATION ',ITER 554 : KERRE = 2 555 : RETURN 556 : ENDIF 557 : C 558 : DLAM1 = (CRINOR * GIGI) - (CRICIS * FIGI) 559 : DLAM1 = DLAM1 / DET 560 : DELAM1 = DELAM1 + DLAM1 561 : RTRAC=0.D0 562 : C 563 : DLAM2 = (CRICIS * FIFI) - (CRINOR * GIFI) 564 : DLAM2 = DLAM2 / DET 565 : C 566 : C CALCUL DU NOUVEL ETAT DE CONTRAINTES 567 : C 568 : S1NP11=S1NP10 - (XMAT(1) * ( (DLAM1*DFIDS1) + (DLAM2*DGIDS1))) 569 : SNNP11=SNNP10 - (XMAT(2) * ( (DLAM1*DFIDSN) + (DLAM2*DZGIDS))) 570 : C 571 : C TEST DE LA CONVERGENCE 572 : C 573 : S2 = 0.0D0 574 : CRICIS = CALG(S1NP11,S2,SNNP11,PHI,C) 575 : CRINOR = SNNP11 - RTRAC 576 : C 577 : IF ((ABS(CRICIS).LE.PRECIS).AND.(ABS(CRINOR).LE.PRECIS)) THEN 578 : C 579 : C----------LES DEUX CRITERES SONT VERIFIES 580 : C 581 : C 582 : IF (IMPRI.EQ.1) THEN 583 : write(6,*) 'LES DEUX CRITERES SONT VERIFIES' 584 : write(6,*) 'CONVERGENCE A L ITERATION',ITER 585 : ENDIF 586 : C 587 : C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS, 588 : C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE 589 : C 590 : C 591 : SIGF(1) = S1NP11 592 : SIGF(2) = SNNP11 593 : C 594 : DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1 595 : c 596 : if ( abs(sigf(2)) .le. presp ) then 597 : c 598 : c deformation plastique du pas (n+1) moins 599 : c deformation plastique du pas (n) 600 : c 601 : DEFP(2) = EPNNP1 - EPSN1 602 : else 603 : c 604 : c le joint se deforme toujours elastiquement en fermeture 605 : c 606 : DEFP(2) = 0.0D0 607 : c 608 : endif 609 : C 610 : C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN LAM1 611 : C 612 : c 613 : c calcul du delta epsilon plastique entre 614 : c le pas actuel et le pas precedent 615 : c 616 : delep1 = defp(1) 617 : delepn = defp(2) 618 : c 619 : c produit contracte pour le calcul de la defo. plasti. cumulee 620 : c 621 : pro1 = ( delep1 * delep1 ) + ( delepn * delepn ) 622 : depeq = sqrt( (2.0D0/3.0D0) * pro1 ) 623 : VARF(1) = EPSEQ + DEPEQ 624 : VARF(4) = XLAM1 + DELAM1 625 : c 626 : c calcul de la deformation normale due a l'ouverture seule 627 : c 628 : if ( abs(snnp11).le.presp ) then 629 : epoun = delepn 630 : VARF(2) = EPOUN 631 : else 632 : epoun = 0.0d0 633 : VARF(2) = EPOUN 634 : endif 635 : c 636 : c etat du joint 637 : c stat = 0 levres en contact, avec ou sans cisaillement 638 : c stat = 1 ouverture du joint pente EC 639 : c stat = 2 compression du joint pente EC 640 : c stat = 3 compression du joint pente EF 641 : c 642 : if ( abs(snnp11).le.presp ) then 643 : stat = 1.D0 644 : if ( abs(epnnp1).le.1.D-8 ) then 645 : stat = 0.D0 646 : endif 647 : VARF(3) = STAT 648 : endif 649 : if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then 650 : stat = 2.D0 651 : VARF(3) = STAT 652 : endif 653 : if (snnp11.le.sigsei) then 654 : stat = 3.D0 655 : VARF(3) = STAT 656 : endif 657 : c 658 : RETURN 659 : ELSE 660 : C 661 : S1NP10=S1NP11 662 : SNNP10=SNNP11 663 : S2 = 0.0D0 664 : CRICIS = CALG(S1NP10,S2,SNNP10,PHI,C) 665 : CRINOR = SNNP10 - RTRAC 666 : ITER = ITER + 1 667 : C 668 : GO TO 2 669 : ENDIF 670 : C 671 : C 672 : C____________________________________________________________________ 673 : C 674 : C 675 : C----------CRITERE DE L'EFFORT NORMAL VIOLE MAIS 676 : C----------CRITERE DE CISAILLEMENT VERIFIE 677 : C 678 : 30 CONTINUE 679 : C 680 : C ECOULEMENT SELON LE CRITERE F A L'ITERATION ITER 681 : C 682 : ITER = 1 683 : DELAM1=0.D0 684 : C 685 : 3 CONTINUE 686 : C 687 : IF (ITER.GT.MAXITE) THEN 688 : WRITE(6,*) ' LE CRITERE NORMAL EST VIOLE ' 689 : WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES' 690 : KERRE = 2 691 : RETURN 692 : ENDIF 693 : C 694 : IF (IMPRI.EQ.1) THEN 695 : write(6,*) 'PASSAGE N. 3' 696 : write(6,*) 'CRITERE DE L EFFORT NORMAL VIOLE' 697 : write(6,*) 'ITERATION',ITER 698 : ENDIF 699 : C 700 : C 701 : C CALCUL DU CORRECTEUR PLASTIQUE 702 : C 703 : C 704 : C DERIVEE DU CRITERE F SUIVANT SIGMA 705 : C 706 : DFIDS1 = 0.0D0 707 : DFIDSN = 1.0D0 708 : C 709 : C MATRICE DE HOOK * DGI/DSIGMA 710 : C 711 : HODF1 = XMAT(1) * DFIDS1 712 : HODFN = XMAT(2) * DFIDSN 713 : C 714 : C FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA 715 : C 716 : FIFI = (DFIDS1 * HODF1) + (DFIDSN * HODFN) 717 : C 718 : C RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR 719 : C PLASTIQUE DLAM1 720 : C 721 : IF (ABS(FIFI).LE.PRECIS) THEN 722 : WRITE(6,*) ' SYSTEME INDETERMINE (F) : ITERATION ',ITER 723 : KERRE = 2 724 : RETURN 725 : ENDIF 726 : C 727 : DLAM1 = CRINOR / FIFI 728 : DELAM1 = DELAM1 + DLAM1 729 : RTRAC =0.D0 730 : C 731 : C CALCUL DU NOUVEL ETAT DE CONTRAINTES 732 : C 733 : S1NP11=S1NP10 - (XMAT(1) * (DLAM1*DFIDS1) ) 734 : SNNP11=SNNP10 - (XMAT(2) * (DLAM1*DFIDSN) ) 735 : C 736 : C LE TEST DE LA CONVERGENCE 737 : C 738 : S2 = 0.0D0 739 : CRICIS = CALG(S1NP11,S2,SNNP11,PHI,C) 740 : CRINOR = SNNP11 - RTRAC 741 : C 742 : IF (ABS(CRINOR).LE.PRECIS) THEN 743 : C 744 : C----------LES DEUX CRITERES SONT VERIFIES 745 : C 746 : C 747 : IF (IMPRI.EQ.1) THEN 748 : write(6,*) 'LES DEUX CRITERES SONT VERIFIES' 749 : write(6,*) 'CONVERGENCE A L ITERATION',ITER 750 : ENDIF 751 : C 752 : C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS, 753 : C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE 754 : C 755 : C 756 : SIGF(1) = S1NP11 757 : SIGF(2) = SNNP11 758 : C 759 : DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1 760 : c 761 : if ( abs(sigf(2)) .le. presp ) then 762 : c 763 : c deformation plastique du pas (n+1) moins 764 : c deformation plastique du pas (n) 765 : c 766 : DEFP(2) = EPNNP1 - EPSN1 767 : else 768 : c 769 : c le joint se deforme toujours elastiquement en fermeture 770 : c 771 : DEFP(2) = 0.0D0 772 : c 773 : endif 774 : C 775 : C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN XLAM1 776 : C 777 : c 778 : c calcul du delta epsilon plastique entre 779 : c le pas actuel et le pas precedent 780 : c 781 : delep1 = defp(1) 782 : delepn = defp(2) 783 : c 784 : c produit contracte pour le calcul de la defo. plasti. cumulee 785 : c 786 : pro1 = ( delep1 * delep1 ) + ( delepn * delepn ) 787 : depeq = sqrt( (2.0D0/3.0D0) * pro1 ) 788 : VARF(1) = EPSEQ + DEPEQ 789 : VARF(4) = XLAM1 + DELAM1 790 : c 791 : c calcul de la deformation normale due a l'ouverture seule 792 : c 793 : if ( abs(snnp11).le.presp ) then 794 : epoun = delepn 795 : VARF(2) = EPOUN 796 : else 797 : epoun = 0.0d0 798 : VARF(2) = EPOUN 799 : endif 800 : c 801 : c etat du joint 802 : c stat = 0 levres en contact, avec ou sans cisaillement 803 : c stat = 1 ouverture du joint pente EC 804 : c stat = 2 compression du joint pente EC 805 : c stat = 3 compression du joint pente EF 806 : c 807 : if ( abs(snnp11).le.presp ) then 808 : stat = 1.D0 809 : if ( abs(epnnp1).le.1.D-8 ) then 810 : stat = 0.D0 811 : endif 812 : VARF(3) = STAT 813 : endif 814 : if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then 815 : stat = 2.D0 816 : VARF(3) = STAT 817 : endif 818 : if (snnp11.le.sigsei) then 819 : stat = 3.D0 820 : VARF(3) = STAT 821 : endif 822 : c 823 : RETURN 824 : ELSE 825 : C 826 : S1NP10=S1NP11 827 : SNNP10=SNNP11 828 : S2 = 0.0D0 829 : CRICIS = CALG(S1NP10,S2,SNNP10,PHI,C) 830 : CRINOR = SNNP10 - RTRAC 831 : ITER = ITER + 1 832 : C 833 : GO TO 3 834 : ENDIF 835 : C 836 : C 837 : C____________________________________________________________________ 838 : C 839 : C 840 : C----------CRITERES DE CISAILLEMENT VIOLE MAIS 841 : C----------CRITERE DE L'EFFORT NORMAL VERIFIE 842 : C 843 : 40 CONTINUE 844 : C 845 : C ECOULEMENT SELON LE CRITERE G A L'ITERATION ITER 846 : C 847 : ITER = 1 848 : C 849 : 4 CONTINUE 850 : C 851 : IF (ITER.GT.MAXITE) THEN 852 : WRITE(6,*) ' LE CRITERE DE CISAILLEMENT EST VIOLE ' 853 : WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES' 854 : KERRE = 2 855 : RETURN 856 : ENDIF 857 : C 858 : IF (IMPRI.EQ.1) THEN 859 : write(6,*) 'PASSAGE N. 4' 860 : write(6,*) 'CRITERE DE CISAILLEMENT VIOLE' 861 : write(6,*) 'ITERATION',ITER 862 : ENDIF 863 : C 864 : C 865 : C CALCUL DU CORRECTEUR PLASTIQUE 866 : C 867 : C 868 : C DERIVEE DU CRITERE G SUIVANT SIGMA 869 : C 870 : if ( s1np10.lt.precis ) then 871 : dgids1 = -1.d0 872 : else 873 : dgids1 = 1.d0 874 : endif 875 : c 876 : DZGIDS = TAN(ZMU) 877 : DGIDSN = TAN(PHI) 878 : C 879 : C MATRICE DE HOOK * DGI/DSIGMA 880 : C 881 : HODG1 = XMAT(1) * DGIDS1 882 : HODGN = XMAT(2) * DZGIDS 883 : C 884 : C GIGI = DGI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA 885 : C 886 : GIGI = (DGIDS1 * HODG1) + (DGIDSN * HODGN) 887 : C 888 : C RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR 889 : C PLASTIQUE DLAM2 890 : C 891 : IF (ABS(GIGI).LE.PRECIS) THEN 892 : WRITE(6,*) ' SYSTEME INDETERMINE (G) : ITERATION ',ITER 893 : KERRE = 2 894 : RETURN 895 : ENDIF 896 : C 897 : DLAM2 = CRICIS / GIGI 898 : C 899 : C CALCUL DU NOUVEL ETAT DE CONTRAINTES 900 : C 901 : S1NP11=S1NP10 - (XMAT(1) * (DLAM2*DGIDS1) ) 902 : SNNP11=SNNP10 - (XMAT(2) * (DLAM2*DZGIDS) ) 903 : c 904 : if (SNNP11.GE.0.0d0) THEN 905 : SNNP11 = 0.0d0 906 : ENDIF 907 : cccssszzzeeerrr if (abs(zmu).le.presp) then 908 : cccssszzzeeerrr differ = abs( depsn ) - epsn1 909 : cccssszzzeeerrr differ = abs(differ)/abs( depsn ) 910 : cccssszzzeeerrr if (depsn.lt.preem) then 911 : cccssszzzeeerrr if (epnnp1.gt.preem) then 912 : cccssszzzeeerrr SNNP11 = 0.0d0 913 : cccssszzzeeerrr endif 914 : cccssszzzeeerrr else 915 : cccssszzzeeerrrc 916 : cccssszzzeeerrr SNNP11 = 0.0d0 917 : cccssszzzeeerrr endif 918 : cccssszzzeeerrr endif 919 : C 920 : C TEST DE LA CONVERGENCE 921 : C 922 : S2 = 0.0D0 923 : CRICIS = CALG(S1NP11,S2,SNNP11,PHI,C) 924 : CRINOR = SNNP11 - RTRAC 925 : C 926 : IF (ABS(CRICIS).LE.PRECIS) THEN 927 : C 928 : IF (INDIC.EQ.1) THEN 929 : C 930 : C----------1ER TRAJET : ON PARCOURT ALP * DEPSN 931 : C 932 : DEFPN = -epin0(2) 933 : S1INT = S1NP11 934 : DEFP11 = EP1NP1 - ( DD(1,1) * S1NP11 ) 935 : GO TO 98 936 : ELSE 937 : IF (IFIN.EQ.0) THEN 938 : C 939 : C----------LES DEUX CRITERES SONT VERIFIES 940 : C 941 : C 942 : IF (IMPRI.EQ.1) THEN 943 : write(6,*) 'LES DEUX CRITERES SONT VERIFIES' 944 : write(6,*) 'CONVERGENCE A L ITERATION',ITER 945 : ENDIF 946 : C 947 : C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS, 948 : C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE 949 : C 950 : C 951 : SIGF(1) = S1NP11 952 : SIGF(2) = SNNP11 953 : C 954 : DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1 955 : c 956 : if ( abs(sigf(2)) .le. presp ) then 957 : c 958 : c deformation plastique du pas (n+1) moins 959 : c deformation plastique du pas (n) 960 : c 961 : DEFP(2) = EPNNP1 - EPSN1 962 : else 963 : c 964 : c le joint se deforme toujours elastiquement en fermeture 965 : c 966 : DEFP(2) = 0.0D0 967 : c 968 : endif 969 : ELSE 970 : C 971 : C----------CAS DU PASSAGE DIRECT D'UN ETAT OUVERT A UN ETAT COMPRIME 972 : C 973 : DEFP12 = EP1NP1 - (DD(1,1)*S1NP11) 974 : DEFP(1) = DEFP12 - epin0(1) 975 : DEFP(2) = -epin0(2) 976 : SIGF(1) = S1NP11 977 : SIGF(2) = SNNP11 978 : IFIN = 0 979 : ENDIF 980 : C 981 : ENDIF 982 : C 983 : C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN XLAM1 984 : C 985 : c 986 : c calcul du delta epsilon plastique entre 987 : c le pas actuel et le pas precedent 988 : c 989 : delep1 = defp(1) 990 : delepn = defp(2) 991 : c 992 : c produit contracte pour le calcul de la defo. plasti. cumulee 993 : c 994 : pro1 = ( delep1 * delep1 ) + ( delepn * delepn ) 995 : depeq = sqrt( (2.0D0/3.0D0) * pro1 ) 996 : VARF(1) = EPSEQ + DEPEQ 997 : VARF(4) = XLAM1 998 : c 999 : c calcul de la deformation normale due a l'ouverture seule 1000 : c 1001 : if ( abs(snnp11).le.presp ) then 1002 : epoun = delepn 1003 : VARF(2) = EPOUN 1004 : else 1005 : epoun = 0.0d0 1006 : VARF(2) = EPOUN 1007 : endif 1008 : c 1009 : c etat du joint 1010 : c stat = 0 levres en contact, avec ou sans cisaillement 1011 : c stat = 1 ouverture du joint pente EC 1012 : c stat = 2 compression du joint pente EC 1013 : c stat = 3 compression du joint pente EF 1014 : c 1015 : if ( abs(snnp11).le.presp ) then 1016 : stat = 1.D0 1017 : if ( abs(epnnp1).le.1.D-8 ) then 1018 : stat = 0.D0 1019 : endif 1020 : VARF(3) = STAT 1021 : endif 1022 : if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then 1023 : stat = 2.D0 1024 : VARF(3) = STAT 1025 : endif 1026 : if (snnp11.le.sigsei) then 1027 : stat = 3.D0 1028 : VARF(3) = STAT 1029 : endif 1030 : c 1031 : RETURN 1032 : ELSE 1033 : C 1034 : S1NP10=S1NP11 1035 : SNNP10=SNNP11 1036 : S2 = 0.0D0 1037 : CRICIS = CALG(S1NP10,S2,SNNP10,PHI,C) 1038 : CRINOR = SNNP10 - RTRAC 1039 : ITER = ITER + 1 1040 : C 1041 : GO TO 4 1042 : ENDIF 1043 : C 1044 : C 1045 : C____________________________________________________________________ 1046 : C 1047 : RETURN 1048 : END 1049 : 1050 : 1051 : 1052 : 1053 : 1054 : 1055 :
© Cast3M 2003 - All rights reserved.
Disclaimer