Download ricbet2d.eso

Back to the list

   1 : C RICBET2D  SOURCE    CHAT      12/08/06    21:15:28     7479           
   2 :           SUBROUTINE RICBET2D(XMAT,SIG0,SIGF,VAR0,VARF,DEPST,ISTEP)
   3 : C
   4 : C====&===1=========2=========3=========4=========5=========6=========7==
   5 : C Commentaires : Subroutine permettant de mettre en oeuvre le 
   6 : C                modele RICBET (B. Richard) pour representer
   7 : C                le comportement 2D d'un béton ordinaire sous                 
   8 : C                sollicitations cycliques
   9 : C
  10 : C Traits       : - Endommagement scalaire 
  11 : C                - Glissement frottant
  12 : C                - 2D
  13 : C                - Monotone et cyclique
  14 : C                - Effet unilatteral
  15 : C                - Comportement linŽaire en compression
  16 : C                - Non localite portant sur le seuil d'endommagement (Fd)
  17 : C
  18 : C Auteurs      : B. Richard (Dr - Ing.) & F. Ragueneau (Pr) - LMT/ENS
  19 : C Date         : Mars 2011
  20 : C====&===1=========2=========3=========4=========5=========6=========7==
  21 : C
  22 : C-----DECLARATION GENERALE----------------------------------------------
  23 : C
  24 :       IMPLICIT REAL*8(A-H,O-Z)
  25 : C
  26 : C-----LISTE DES INCLUDES A CHARGER--------------------------------------
  27 : C
  28 : C -INC CCOPTIO 
  29 : C -INC DECHE
  30 : C
  31 : C-----DECLARATION GENERALE----------------------------------------------
  32 : C
  33 :       REAL*8 X(3,3),UNIT(3,3),EPSPI(3,3),EPSI(3,3),SIGMA(3,3)
  34 :       REAL*8 EPSDEV(3,3),SIGMA1(3,3),EPSPIDEV(3,3)
  35 :       REAL*8 EPSIPOS(3,3),EPSINEG(3,3),SIGPI(3,3),DEPSN(3,3)
  36 :       REAL*8 TERME1(3,3),TERME2(3,3),TERME3(3,3),TERME4(3,3)
  37 :       REAL*8 SIGMOX(3,3),TRAV(3,3),SIG1DEV(3,3),TERME5(3,3)
  38 :       REAL*8 XDT,XTRAV2(3,3),XTRAV3(3,3),DEFN(3,3),EPSIT(3,3)     
  39 :       REAL*8 XMAT(*),SIG0(*),SIGF(*),VAR0(*),VARF(*),DEPST(*)
  40 :       REAL*8 XDFP(3,3),XDGP(3,3),XDJ2P(3,3),DEVGP(3,3)
  41 :       REAL*8 DEPSI(3,3)
  42 : C
  43 : C-----DEFINITION DE PARAMETRES------------------------------------------
  44 : C
  45 :       PARAMETER (IMAX = 1000)
  46 :       PARAMETER (IMAX2 = 10)
  47 : C
  48 : C-----MISE EN DONNEE----------------------------------------------------
  49 : C             
  50 :       XE = XMAT(1)
  51 :       XNU = XMAT(2)
  52 :       XFT = XMAT(5)
  53 :       ALDIR = XMAT(6)
  54 :       GAMMA1 = XMAT(7)
  55 :       A1 = XMAT(8)
  56 :       SIGRF = XMAT(9)
  57 :       
  58 :       XAF = XMAT(10)
  59 :       XBF = XMAT(11)
  60 :       XAG = XMAT(12)
  61 :       XBG = XMAT(13)
  62 :       XAC = XMAT(14)
  63 :       XBC = XMAT(15)
  64 :       XSIGU = XMAT(16)
  65 :       XFC = XMAT(17)
  66 : 
  67 :       XK = XE/((1.0D0-2.D0*XNU)*3.0D0)
  68 :       XG = XE/(2.0D0*(1.D0+XNU))
  69 : C
  70 : C-----DEFINITION DE LA MATRICE UNITE D ORDRE 2--------------------------
  71 : C
  72 :           DO I=1,3
  73 :          DO J=1,3
  74 :             IF (I.EQ.J) THEN
  75 :                UNIT(I,J) = 1.0D0
  76 :             ELSE
  77 :                UNIT(I,J) = 0.0D0
  78 :             ENDIF
  79 :          ENDDO
  80 :       ENDDO
  81 : C
  82 : C-----ACTUALISATION DES DEFORMATIONS------------------------------------
  83 : C       
  84 : 
  85 : C     Indicateur
  86 :       XBEN = VAR0(22)+VAR0(23)+VAR0(24)*0.0D0
  87 : 
  88 : C     Boucle pour les contraintes planes
  89 :       DO NIDX=1,IMAX2
  90 :       
  91 :       IF (NIDX.EQ.IMAX) THEN
  92 : C      print*,'non converge'
  93 : C      stop
  94 :       ENDIF
  95 :       
  96 :       IF (DABS(DEPST(3)).GE.1.0D10) THEN
  97 : C         print*,'pb depst'
  98 :       ENDIF
  99 :       
 100 : C     Contraintes tests      
 101 :       SIGMA1(1,1) = VAR0(22)
 102 :       SIGMA1(2,2) = VAR0(23)
 103 :       SIGMA1(3,3) = VAR0(24)
 104 :       SIGMA1(1,2) = VAR0(25)
 105 :       SIGMA1(1,3) = VAR0(26)
 106 :       SIGMA1(2,3) = VAR0(27)
 107 :       SIGMA1(2,1) = VAR0(25)
 108 :       SIGMA1(3,1) = VAR0(26)
 109 :       SIGMA1(3,2) = VAR0(27)
 110 :       
 111 : C     Tenseur de glissement interne
 112 :       EPSPI(1,1) = VAR0(2)
 113 :       EPSPI(2,2) = VAR0(3)
 114 :       EPSPI(3,3) = VAR0(4)
 115 :       EPSPI(1,2) = VAR0(5)
 116 :       EPSPI(1,3) = VAR0(6)
 117 :       EPSPI(2,3) = VAR0(7)
 118 :       EPSPI(2,1) = VAR0(5)
 119 :       EPSPI(3,1) = VAR0(6)
 120 :       EPSPI(3,2) = VAR0(7)
 121 : 
 122 : C     Tenseur des contraintes de rappel
 123 :       X(1,1) = VAR0(8)
 124 :       X(2,2) = VAR0(9)
 125 :       X(3,3) = VAR0(10)
 126 :       X(1,2) = VAR0(11)
 127 :       X(1,3) = VAR0(12)
 128 :       X(2,3) = VAR0(13)
 129 :       X(2,1) = VAR0(11)
 130 :       X(3,1) = VAR0(12)
 131 :       X(3,2) = VAR0(13)
 132 : 
 133 : C     Tenseur des deformations plastiques
 134 :       DEFN(1,1) = VAR0(35)
 135 :       DEFN(2,2) = VAR0(36)
 136 :       DEFN(3,3) = VAR0(37)
 137 :       DEFN(1,2) = VAR0(38)
 138 :       DEFN(1,3) = VAR0(39)
 139 :       DEFN(2,3) = VAR0(40)
 140 :       DEFN(2,1) = VAR0(38)
 141 :       DEFN(3,1) = VAR0(39)
 142 :       DEFN(3,2) = VAR0(40)
 143 : 
 144 : C     Ecrouissage plastique
 145 :       XECO = VAR0(41)
 146 : 
 147 : C     Deformation plastique cumulee
 148 :       XDPP = VAR0(42)
 149 : 
 150 : C     Endommagement
 151 :       XDT = VAR0(14)
 152 : 
 153 : C     Force associee a l ecrouissage isotrope
 154 :       XFFT = VAR0(15)
 155 : 
 156 : C     Deformations totales      
 157 :       EPSI(1,1) = VAR0(16)
 158 :       EPSI(2,2) = VAR0(17)
 159 :       EPSI(3,3) = VAR0(18)
 160 :       EPSI(1,2) = VAR0(19)
 161 :       EPSI(1,3) = VAR0(20)
 162 :       EPSI(2,3) = VAR0(21)
 163 :       EPSI(2,1) = VAR0(19)
 164 :       EPSI(3,1) = VAR0(20)
 165 :       EPSI(3,2) = VAR0(21)
 166 : 
 167 :       EPSI(1,1) = EPSI(1,1) + DEPST(1)
 168 :       EPSI(2,2) = EPSI(2,2) + DEPST(2)
 169 :       EPSI(3,3) = EPSI(3,3) + DEPST(3)
 170 :       EPSI(1,2) = EPSI(1,2) + (DEPST(4)/2.0D0)
 171 :       EPSI(1,3) = 0.0D0
 172 :       EPSI(2,3) = 0.0D0
 173 :       EPSI(2,1) = EPSI(2,1) + (DEPST(4)/2.0D0)
 174 :       EPSI(3,1) = 0.0D0
 175 :       EPSI(3,2) = 0.0D0
 176 : 
 177 :       DEPSI(1,1) = DEPST(1)
 178 :       DEPSI(2,2) = DEPST(2)
 179 :       DEPSI(3,3) = DEPST(3)
 180 :       DEPSI(1,2) = (DEPST(4)/2.0D0)
 181 :       DEPSI(1,3) = 0.0D0
 182 :       DEPSI(2,3) = 0.0D0
 183 :       DEPSI(2,1) = (DEPST(4)/2.0D0)
 184 :       DEPSI(3,1) = 0.0D0
 185 :       DEPSI(3,2) = 0.0D0
 186 : C
 187 : C-----CALCUL DE QUELQUES QUANTITES--------------------------------------
 188 : C
 189 : 
 190 : C     Calcul de la trace de EPSI
 191 :       TRA = EPSI(1,1)+EPSI(2,2)+EPSI(3,3)
 192 : 
 193 : C     Calcul du deviateur de EPSI
 193 :       CALL DEVIAT(EPSI,EPSDEV)
 194 :           
 195 : C     Calcul de la trace de EPSPI
 196 :       TRAPI = EPSPI(1,1)+EPSPI(2,2)+EPSPI(3,3)
 197 :       
 198 : C     Calcul du deviateur de EPSPI
 198 :       CALL DEVIAT(EPSPI,EPSPIDEV)
 199 : 
 200 : C     Calcul de la trace de SIGMA1
 201 :       TRAS1 = SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3)
 202 : 
 203 : C     Calcul du deviateur de DEFN
 203 :       CALL DEVIAT(DEFN,DEPSN)
 204 : 
 205 : C     Calcul de la trace de DEFN
 206 :       TRADN = DEFN(1,1)+DEFN(2,2)+DEFN(3,3)
 207 : 
 208 : C
 209 : C-----PARTIE POSITIVES ET NEGATIVES DE (EPS-EPSP)------------------------
 210 : C
 211 :       DO I=1,3
 212 :          DO J=1,3
 213 :             EPSIT(I,J) = EPSI(I,J)-DEFN(I,J)
 214 :          ENDDO
 215 :       ENDDO
 216 : 
 216 :       CALL POSNEG(EPSIT,EPSIPOS,EPSINEG)
 217 : C
 218 : C-----SEUIL ENDOMMAGEMENT------------------------------------------------
 219 : C
 220 : 
 221 : C     Calcul de l enerigie Y+
 222 :       YP = 0.0D0
 223 :       
 224 :       DO I=1,3
 225 :          DO J=1,3
 226 :             YP = YP + EPSIPOS(I,J)*EPSIPOS(I,J)
 227 :          ENDDO
 228 :       ENDDO
 229 : 
 230 :       YP = 0.5D0*YP*XE      
 231 : C
 232 : C-----NON LOCAL----------------------------------------------------------
 233 : C
 234 :       IF (ISTEP.EQ.0) THEN
 235 :          SEUIL1=YP
 236 :          VARF(1)=YP
 237 :       ELSE IF (ISTEP.EQ.1) THEN
 238 :          VARF(1)=YP
 239 :          DO I=2,42
 240 :             VARF(I)=VAR0(I)
 241 :          ENDDO
 242 :          GOTO 2000
 243 :       ELSE IF (ISTEP.EQ.2) THEN
 244 :          SEUIL1=VAR0(1)
 245 :          VARF(1)=SEUIL1
 246 :       ENDIF 
 247 : C
 248 : C-----EVOLUTION DE L ENDOMMAGEMENT---------------------------------------
 249 : C
 250 : 
 251 : C     Calcul du seuil
 252 :       XFD = SEUIL1 - ((XFT**2)/(2.0D0*XE) + XFFT)
 253 : 
 254 :       TRAS1 = 1.0D0
 255 :       
 256 : C     Evolution explicite de l'endommagement
 257 :       IF ((XFD.GT.0.0D0).AND.(XBEN.GE.0.0D0)) THEN
 258 :          XDT  = 1.00D0-1.00D0/
 259 :      &   (1.0D0+ALDIR*(SEUIL1-((XFT**2)/(2.0D0*XE))))
 260 :          XFFT = SEUIL1 - ((XFT**2)/(2.0D0*XE))
 261 :       ENDIF   
 262 :       
 263 : C
 264 : C-----SEUIL GLISSEMENT INTERNE-------------------------------------------
 265 : C
 266 : 
 267 : C     Calcul de sigma_pi
 268 :       DO I=1,3
 269 :          DO J=1,3
 270 :             SIGPI(I,J) = XDT*XK*(TRA-TRAPI-TRADN)*UNIT(I,J)+
 271 :      &      2.0D0*XG*XDT*(EPSDEV(I,J)-EPSPIDEV(I,J)-DEPSN(I,J))
 272 :          ENDDO
 273 :       ENDDO
 274 : 
 275 : C     Calcul de la quantite (sigma_pi - X)
 276 :       DO I=1,3
 277 :          DO J=1,3
 278 :             SIGMOX(I,J) = SIGPI(I,J)-X(I,J)
 279 :          ENDDO
 280 :       ENDDO
 281 : 
 282 : C     Calcul de I2(sigma_pi-X)
 282 :       CALL I2(SIGMOX,XFPI)
 283 : 
 284 : C     Mise a 1 du multiplicateur plastique
 285 :       XLAM = 1.0D0
 286 : 
 287 : C     CAlcul de la trace de Sigma_finale  
 288 :       TRASIG1 = SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3)
 289 : 
 290 : C     Test sur la positivite du seuil de glissement interne
 291 :       IF ((XFPI.GT.0.0D0).AND.(XBEN.GE.0.D0)) THEN   
 292 :           
 293 : C         On retient la valeur initiale du seuil
 294 :           XFPI0 = XFPI
 295 : 
 296 : C         On itere en bloquant le nombre d iteration max a IMAX
 297 :           DO NI=1,IMAX
 298 : 
 299 : C             On realise un test de sortie
 300 :               IF ((DABS(XFPI/XFPI0 ).LE.(1.0D-8)).OR.
 301 :      &            (DABS(XLAM).LE.(1.0D-20))) THEN
 302 : 
 303 :                  GOTO 400
 304 : 
 305 :               ELSE
 306 : 
 307 : C                Calcul des derivees dfpi/dsigpi et dphipi/dsigpi
 307 :                  CALL DI2(SIGMOX,TERME1)
 307 :                  CALL DI2(SIGMOX,TERME2)
 308 : 
 309 : C                Calcul des derivees dfpi/dX et dphipi/dX
 310 :                  DO I=1,3
 311 :                     DO J=1,3
 312 :                        TERME3(I,J) = 0.0D0
 313 :                        TERME4(I,J) = 0.0D0
 314 : 
 315 :                        TERME3(I,J) = -1.0D0*TERME2(I,J)
 316 :                        TERME4(I,J) = -1.0D0*TERME2(I,J)+A1*X(I,J)
 317 :                     ENDDO
 318 :                  ENDDO
 319 : 
 320 : C                On calcul le terme C:dfpi/dsigpi
 321 :                  TRA_SE = TERME2(1,1)+TERME2(2,2)+TERME2(3,3)
 322 :                  
 322 :                  CALL DEVIAT(TERME2,TRAV)
 323 : 
 324 :                  DO I=1,3
 325 :                     DO J=1,3
 326 :                        TERME5(I,J) = 0.0D0
 327 :                        TERME5(I,J) = XK*TRA_SE*UNIT(I,J)+
 328 :      &                 2.0D0*XG*TRAV(I,J)
 329 :                     ENDDO
 330 :                  ENDDO
 331 : 
 332 : C                On calcul le denominateur du multiplicateur plastique
 333 :                  XT1 = 0.0D0
 334 :                  DO I=1,3
 335 :                     DO J=1,3
 336 :                        XT1 = XT1 + TERME1(I,J)*TERME5(I,J)
 337 :                     ENDDO
 338 :                  ENDDO
 339 :                  
 340 :                  XT2 = 0.0D0                 
 341 :                  DO I=1,3
 342 :                     DO J=1,3
 343 :                        XT2 = XT2 + TERME3(I,J)*TERME4(I,J)
 344 :                     ENDDO
 345 :                  ENDDO
 346 : 
 347 : C                Calcul du multiplicateur plastique
 348 :                  XLAM = XFPI/(XDT*XT1 + GAMMA1*XT2)
 349 : 
 350 : C                On actualisatise X et sigpi
 351 :                  DO I=1,3
 352 :                     DO J=1,3
 353 :                        SIGPI(I,J)=SIGPI(I,J)-XDT*XLAM*TERME5(I,J)
 354 :                        X(I,J)=X(I,J)-GAMMA1*XLAM*TERME4(I,J)
 355 :                     ENDDO
 356 :                  ENDDO
 357 : 
 358 : C                On recalcule le nouveau seuil
 359 :                  DO I=1,3
 360 :                     DO J=1,3
 361 :                        SIGMOX(I,J) = 0.0D0
 362 :                        SIGMOX(I,J) = SIGPI(I,J)-X(I,J)
 363 :                     ENDDO
 364 :                  ENDDO
 365 : 
 365 :                  CALL I2(SIGMOX,XFPI)
 366 : 
 367 :               ENDIF
 368 : 
 369 :           ENDDO
 370 :           
 371 :  400      CONTINUE 
 372 : 
 373 : 
 374 : C         On calcul epspi
 375 :           TRA_SE = SIGPI(1,1)+SIGPI(2,2)+SIGPI(3,3)
 375 :           CALL DEVIAT(SIGPI,TRAV)
 376 : 
 377 :           TRAPI = (TRA-TRADN) - (TRA_SE/(XDT*3.0D0*XK))
 378 : 
 379 :           DO I=1,3
 380 :              DO J=1,3
 381 :                 EPSPIDEV(I,J) = (EPSDEV(I,J)-DEPSN(I,J))-
 382 :      &          (1.0D0/(2.0D0*XG*XDT))*TRAV(I,J)
 383 :              ENDDO
 384 :           ENDDO
 385 : 
 386 :           DO I=1,3
 387 :              DO J=1,3
 388 :                 EPSPI(I,J) = (TRAPI/3.0D0)*UNIT(I,J)+EPSPIDEV(I,J)
 389 :              ENDDO
 390 :           ENDDO
 391 : 
 392 :       ENDIF
 393 : C
 394 : C-----CONTRAINTES ACTUALISEE
 395 : C
 396 :       TRA1 = SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3)
 397 :             
 398 :       IF (XBEN.GE.0.0D0) THEN
 399 : 
 400 :  752  CONTINUE         
 401 :       
 402 :          DO I=1,3
 403 :             DO J=1,3
 404 :                SIGMA1(I,J)=XK*(TRA-XDT*TRAPI-TRADN)*UNIT(I,J)+
 405 :      &          2.0D0*XG*
 406 :      &          (EPSDEV(I,J)-XDT*EPSPIDEV(I,J)-DEPSN(I,J))
 407 :             ENDDO
 408 :          ENDDO
 409 :             
 410 :       ELSEIF (XBEN.LE.SIGRF) THEN
 411 : 
 412 : C        La contrainte effective
 413 :          DO I=1,3
 414 :             DO J=1,3
 415 :                SIGMA1(I,J)=XK*(TRA-TRADN)*UNIT(I,J)+
 416 :      &         2.0D0*XG*(EPSDEV(I,J)-DEPSN(I,J))
 417 :             ENDDO
 418 :          ENDDO
 419 : 
 420 : C        On verifie le seuile de plasticite
 420 :          CALL J2(SIGMA1,TEMP1)
 421 : 
 422 :          TEMP2 = SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3)
 423 : 
 424 :          TEMP2 = TEMP2/3.00
 425 :  
 426 :          XFP = XAF*TEMP1+XBF*TEMP2 - (XFC + XECO)
 427 : 
 428 : C        On verifie la positivite du seuil
 429 :          IF (XFP.GT.0.0D0) THEN
 430 :          
 431 : C           On retient la valeur initiale pour le critere relatif
 432 :             XFP0 = XFP
 433 : 
 434 : C           On commence les choses serieuses
 435 :             DO NP=1,IMAX
 436 :             
 437 : C           La derivee du critere
 437 :             CALL DI2(SIGMA1,XDJ2P)
 438 :             
 439 :             DO I=1,3
 440 :             DO J=1,3
 441 :             XDFP(I,J) = XAF*XDJ2P(I,J)+XBF*UNIT(I,J)/3.0D0
 442 :             XDGP(I,J) = XAG*XDJ2P(I,J)+XBG*UNIT(I,J)/3.0D0
 443 :             ENDDO
 444 :             ENDDO
 445 : 
 446 : C           Le hessien elasique             
 446 :             CALL DEVIAT(XDGP,DEVGP)
 447 : 
 448 :             XTGP = XDGP(1,1)+XDGP(2,2)+XDGP(3,3)
 449 :             
 450 :             DO I=1,3
 451 :             DO J=1,3
 452 :             TRAV(I,J)=0.0D0
 453 :             TRAV(I,J)=XK*XTGP*UNIT(I,J)+2.0D0*XG*DEVGP(I,J)
 454 :             ENDDO
 455 :             ENDDO
 456 : 
 457 : C           La derivee de l ecrouissage            
 458 :             XTEMP= DEXP(-1.0D0*XBC*XDPP)*(XAC-(XAC*XDPP+XFC+XSIGU)*XBC)
 459 : 
 460 : C           Le critere actuel 
 461 :            
 461 :             CALL J2(SIGMA1,TEMP1)
 462 : 
 463 :             TEMP2 = SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3)
 464 : 
 465 :             TEMP2 = TEMP2/3.00
 466 :  
 467 :             XFP = XAF*TEMP1+XBF*TEMP2 - (XFC + XECO)
 468 : 
 469 : C           Le multiplicateur plastique
 469 :             CALL DBLECONT(XDFP,TRAV,TEMP9)
 470 : 
 471 :             XLAMP = XFP/(TEMP9+XTEMP)
 472 : 
 473 : C           Actualisation des variables internes
 474 :             XECO = XECO + XLAMP*XTEMP
 475 :             
 476 :             XDPP = XDPP + XLAMP
 477 :             
 478 :             DO I=1,3
 479 :             DO J=1,3
 480 :             DEFN(I,J)=DEFN(I,J) + XLAMP*XDGP(I,J)
 481 :             ENDDO
 482 :             ENDDO
 483 : 
 484 :             DO I=1,3
 485 :             DO J=1,3
 486 :             SIGMA1(I,J)=SIGMA1(I,J) - XLAMP*TRAV(I,J)
 487 :             ENDDO
 488 :             ENDDO
 489 :             
 490 : C           On calcul le seuil pour un eventuel arret
 490 :             CALL J2(SIGMA1,TEMP1)
 491 :             TEMP2 = SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3)
 492 :             TEMP2 = TEMP2/3.00
 493 :  
 494 :             XFP = XAF*TEMP1+XBF*TEMP2 - (XFC + XECO)
 495 :             
 496 : C           Le critere d arret des iterations internes            
 497 :             IF ((DABS(XFP/XFP0).LE.(1.0D-8)).OR.(NP.GE.IMAX)) THEN
 498 :                
 499 :                GOTO 294
 500 :                
 501 :             ENDIF
 502 :             
 503 :             
 504 :             ENDDO
 505 :          
 506 :  294     CONTINUE        
 507 :          
 508 :          ENDIF
 509 : 
 510 :       ELSE
 511 : 
 512 : C        Calcul du deviateur de DEFN
 512 :          CALL DEVIAT(DEFN,DEPSN)
 513 : 
 514 : C        Calcul de la trace de DEFN
 515 :          TRADN = DEFN(1,1)+DEFN(2,2)+DEFN(3,3)
 516 : 
 517 :          TRASIG1 = 3.0D0*XK*(TRA-XDT*TRAPI-TRADN)/
 518 :      &             (1.0D0-(3.0D0*XK/SIGRF)*XDT*TRAPI)
 519 :               
 520 :          DO I=1,3
 521 :             DO J=1,3
 522 :                SIG1DEV(I,J)=2.00D0*XG*(EPSDEV(I,J)-XDT*EPSPIDEV(I,J)*
 523 :      &         (1.0D0-(TRASIG1/SIGRF))-DEPSN(I,J))
 524 :             ENDDO
 525 :          ENDDO
 526 :          
 527 :          DO I=1,3
 528 :             DO J=1,3
 529 :                SIGMA1(I,J)=TRASIG1*UNIT(I,J)/3.0D0+SIG1DEV(I,J)
 530 :             ENDDO
 531 :          ENDDO
 532 :          
 533 :       ENDIF 
 534 :       
 535 :  101  CONTINUE     
 536 :       
 537 :       
 538 :       IF (DABS(SIGMA1(1,1)).GE.1.0D50) THEN
 539 :          print*,'TRASIG1',TRASIG1
 540 :          print*,'QUANT=',(1.0D0-(3.0D0*XK/SIGRF)*XDT*TRAPI)
 541 :          print*,'HEHEF=',(TRA-XDT*TRAPI-TRADN)
 542 :          print*,'XEXEF=',TRAPI
 543 :          print*,'XEXEF=',TRADN
 544 :          print*,'FEFEF=',XDT
 545 :          print*,'TRA=',TRA
 546 :          print*,'DEPST',(DEPST(I),I=1,6)
 547 :       ENDIF
 548 : C
 549 : C-----CALCUL DE L ENDOMMAGEMENT EFFECTIF (GRANDEUR INDICATIVE)
 550 : C-----PAS DE SENS THERMO
 551 : C
 552 :       XFSIG = 1.0D0-(SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3))/SIGRF
 553 :       
 554 :       IF ((XFSIG.LE.1.0D0).AND.(XFSIG.GE.0.0D0)) THEN
 555 :           XDEFF = XDT*XFSIG
 556 :       ELSE
 557 :           XDEFF = XDT
 558 :       ENDIF  
 559 :       
 560 :  444  CONTINUE          
 561 :           
 562 :       IF (DABS(SIGMA1(3,3)).LE.1.0D0) THEN
 563 : 
 564 :          GOTO 678
 565 : 
 566 :       ELSE
 567 : 
 568 :          DEPST(3) = -SIGMA1(3,3)/XE
 569 :          IF (DABS(DEPST(3)).GE.1.0D10) THEN
 570 :             print*,'SIGMA1(3,3)=',SIGMA1(3,3)
 571 :             print*,'soucis d actualisation'
 572 :             print*,'XE=',XE
 573 :             print*,'DPEST3',DEPST(3)
 574 :             
 575 :             
 576 :          ENDIF
 577 :          
 578 :          VAR0(18) = EPSI(3,3)
 579 :          
 580 : C     Tenseur de glissement interne
 581 :       VAR0(2) = EPSPI(1,1)
 582 :       VAR0(3) = EPSPI(2,2)
 583 :       VAR0(4) = EPSPI(3,3)
 584 :       VAR0(5) = EPSPI(1,2)
 585 :       VAR0(6) = EPSPI(1,3)
 586 :       VAR0(7) = EPSPI(2,3)
 587 : 
 588 : C     Tenseur de la contrainte de rappel
 589 :       VAR0(8) = X(1,1)
 590 :       VAR0(9) = X(2,2)
 591 :       VAR0(10) = X(3,3)
 592 :       VAR0(11) = X(1,2)
 593 :       VAR0(12) = X(1,3)
 594 :       VAR0(13) = X(2,3)
 595 : 
 596 : C     Tenseur des deformations plastiques
 597 :       VAR0(35) = DEFN(1,1)
 598 :       VAR0(36) = DEFN(2,2)
 599 :       VAR0(37) = DEFN(3,3)
 600 :       VAR0(38) = DEFN(1,2)
 601 :       VAR0(39) = DEFN(1,3)
 602 :       VAR0(40) = DEFN(2,3)
 603 : 
 604 : C     Endommagement
 605 :       VAR0(14) = XDT
 606 : 
 607 : C     Force associee a l ecrouissage isotrope
 608 :       VAR0(15) =  XFFT
 609 : 
 610 : C     Contraintes tests      
 611 :       VAR0(22) =  SIGMA1(1,1)
 612 :       VAR0(23) =  SIGMA1(2,2)
 613 :       VAR0(24) =  SIGMA1(3,3)
 614 :       VAR0(25) =  SIGMA1(1,2)
 615 :       VAR0(26) =  SIGMA1(1,3)
 616 :       VAR0(27) =  SIGMA1(2,3)
 617 : 
 618 : C     Ecrouissage plastique
 619 :       VAR0(41) = XECO
 620 : 
 621 : C     Deformation plastique cumulee
 622 :       VAR0(42) = XDPP
 623 : 
 624 :       ENDIF
 625 :       
 626 :       ENDDO
 627 :       
 628 :  678  CONTINUE
 629 :          
 630 : C
 631 : C-----STOCKAGE DES CONTRAINTES AVEC SYMMETRISATION POUR ROBUSTESSE
 632 : C-----ON MET SIG_33 BIEN QU ELLE SOIT CALCULEE A UNE TOL. PRES
 633 : C
 634 :       SIGF(1)= SIGMA1(1,1)
 635 :       SIGF(2)= SIGMA1(2,2)
 636 :       SIGF(3)= SIGMA1(3,3)*0.0D0
 637 :       SIGF(4)= (SIGMA1(1,2)+SIGMA1(2,1))/2.0D0
 638 : 
 639 : C
 640 : C-----STOCKAGE DES DEFORMATIONS DE GLISSEMENT
 641 : C
 642 :       VARF(2)=EPSPI(1,1)
 643 :       VARF(3)=EPSPI(2,2)
 644 :       VARF(4)=EPSPI(3,3)
 645 :       VARF(5)=EPSPI(1,2)
 646 :       VARF(6)=EPSPI(1,3)
 647 :       VARF(7)=EPSPI(2,3)
 648 : C
 649 : C-----TENSEUR DES DEFORMATIONS PLASTIQUES
 650 : C     
 651 :       VARF(35) = DEFN(1,1)
 652 :       VARF(36) = DEFN(2,2)
 653 :       VARF(37) = DEFN(3,3)
 654 :       VARF(38) = DEFN(1,2)
 655 :       VARF(39) = DEFN(1,3)
 656 :       VARF(40) = DEFN(2,3)
 657 : C      
 658 : C-----Ecrouissage plastique
 659 : C
 660 :       VARF(41) = XECO
 661 : C
 662 : C-----Deformation plastique cumulee
 663 : C
 664 :       VARF(42) = XDPP
 665 : C
 666 : C-----STOCKAGE DE L ECROUISSAGE CINEMATIQUE
 667 : C
 668 :       VARF(8)=X(1,1)
 669 :       VARF(9)=X(2,2)
 670 :       VARF(10)=X(3,3)
 671 :       VARF(11)=X(1,2)
 672 :       VARF(12)=X(1,3)
 673 :       VARF(13)=X(2,3)
 674 : C
 675 : C-----STOCKAGE DES CONTRAINTES TESTS
 676 : C
 677 :       VARF(22)=SIGMA1(1,1)
 678 :       VARF(23)=SIGMA1(2,2)
 679 :       VARF(24)=SIGMA1(3,3)
 680 :       VARF(25)=SIGMA1(1,2)
 681 :       VARF(26)=SIGMA1(1,3)
 682 :       VARF(27)=SIGMA1(2,3)
 683 : C
 684 : C-----STOCKAGE DES CONTRAINTES TESTS
 685 : C
 686 :       VARF(29)=SIGPI(1,1)
 687 :       VARF(30)=SIGPI(2,2)
 688 :       VARF(31)=SIGPI(3,3)
 689 :       VARF(32)=SIGPI(1,2)
 690 :       VARF(33)=SIGPI(1,3)
 691 :       VARF(34)=SIGPI(2,3)
 692 : C
 693 : C-----ENDOMMAGEMENT EFFECTIF
 694 : C
 695 :       VARF(28) = XDEFF
 696 : C
 697 : C-----STOCKAGE DE L ENDOMMAGEMENT
 698 : C
 699 :       VARF(14) = XDT
 700 : C
 701 : C-----STOCKAGE DES FORCES THERMODYNAMIQUES ASSOCIEES A L ECROUISSAGE ISOTROPE
 702 : C
 703 :       VARF(15) = XFFT
 704 : C
 705 : C-----STOCKAGE DES DEFORMATIONS TOTALES
 706 : C
 707 :       VARF(16)=EPSI(1,1)
 708 :       VARF(17)=EPSI(2,2)
 709 :       VARF(18)=EPSI(3,3)
 710 :       VARF(19)=EPSI(1,2)
 711 :       VARF(20)=EPSI(1,3)
 712 :       VARF(21)=EPSI(2,3)
 713 : C
 714 : C-----BALISE DE SORTIE POUR LE NON LOCAL
 715 : C
 716 :  2000 CONTINUE
 717 : 
 718 :       RETURN
 719 :           END
 720 : 
 721 : 
 722 : 
 723 : 
 724 : 
 725 : 
 726 : 
 727 : 
 728 : 
 729 : 
 730 : 
 731 : 
 732 : 
 733 : 
 734 : 
 735 : 
 736 : 
 737 : 
 738 : 
 739 :  
 740 :  

© Cast3M 2003 - All rights reserved.
Disclaimer