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