Sources Esope | Cast3M

rotkta
C ROTKTA    SOURCE    PV        11/03/07    21:18:08     6885
      SUBROUTINE ROTKTA(XMAT,VARF,SIGF,XKTAT,KERRE)
C
C Rotating Crack model avec 2 directions de fissuration
C Calcul de la matrice tangente
C  D. Combescure 24/08/95 ELSA-Ispra
C
C ENTREES:
C -------
C     XMAT    = Caracteristiques materiaux
C
C     VARF    = Variables internes
C
C     SIGF    = Contraintes
C
C SORTIES:
C -------
C     XKTAT    = Matrice tangente
C
C================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
      DIMENSION XMAT(8),VARF(7),SIGF(4),XXX(2)
      DIMENSION DAM(2), XKTAT(4,4),XKTAT2(3,3),EPSF(2),XMINT(3,3)
      DIMENSION EPSN(2), EPS(2), XKTAP(3,3), XKINV(2,2)
      DIMENSION TP(3,3),TM(3,3),IMEC(2),EPSP(2), SIGM(2)
      PARAMETER (UN=1.D0,UNDEMI=0.5D0,EPSILO=1.D-20,XZER=0.D0)
      KERRE=0
C===============================================================
C   LECTURE DES CARACTERISTIQUES MAT. ET DES VAR. INT.
C===============================================================
      YOUN = XMAT(1)
      XNU  = XMAT(2)
      XFTR = XMAT(3)
      EPSR = XMAT(4)
      XFRE = XMAT(5)
      BETA = XMAT(6)
      XHH1 = (XFRE - XFTR)/EPSR
C--------------------------------------------------------------
C Variables internes->1 et 2- cos et sin de l'angle
C                     3 et 4- "Deformation" maximale de la fissure
C                     5 et 6- Deformations principales
C                     7     - Cas pour le calcul de Kta
C--------------------------------------------------------------
      KERRE   = 0
      CPHI    = VARF(1) + UN
      SPHI    = VARF(2)
      EPS(1)  = VARF(3)
      EPS(2)  = VARF(4)
      EPSN(1) = VARF(5)
      EPSN(2) = VARF(6)
      ICAS    = nint(VARF(7))
      CALL ZERO(XKTAT,4,4)
      CALL ZERO(XKTAT2,3,3)
C
C  Remplissage de la matrice de transformation
C (Ecrite pour le tenseur des deformations contenant 2*EPSXY)
C
      TP(1,1) = CPHI**2
      TP(1,2) = SPHI**2
      TP(1,3) = -SPHI*CPHI
      TP(2,1) = TP(1,2)
      TP(2,2) = TP(1,1)
      TP(2,3) = -TP(1,3)
      TP(3,1) = 2.*SPHI*CPHI
      TP(3,2) = -2.*SPHI*CPHI
      TP(3,3) = (CPHI**2) - (SPHI**2)
C
C  Ecrite pour les contraintes
C
      TM(1,1) = CPHI**2
      TM(1,2) = SPHI**2
      TM(1,3) = 2.*SPHI*CPHI
      TM(2,1) = TM(1,2)
      TM(2,2) = TM(1,1)
      TM(2,3) = -TM(1,3)
      TM(3,1) = -SPHI*CPHI
      TM(3,2) = SPHI*CPHI
      TM(3,3) = (CPHI**2) - (SPHI**2)
C
C
C================================================================
C   MATRICE TANGENTE DANS LE REPERE PRINCIPAL
C================================================================
C
C  Calcul des variables d'endommagement
      DO I=1,2
        IF (EPS(I).LE.EPSR) THEN
          SIGM(I) = XFTR - (XFTR - XFRE)*EPS(I)/EPSR
          EPSP(I) = BETA*EPS(I)
          XXX(I) = UN + YOUN*(UN - BETA)*EPS(I)/SIGM(I)
        ELSE
          SIGM(I) = XFRE
          EPSP(I) = BETA*EPSR
          XXX(I) = UN + YOUN*(UN - BETA)*EPSR/XFRE
     &+ YOUN*(EPS(I) - EPSR)/XFRE
        ENDIF
        XXX(I) = UN/XXX(I)
        DAM(I) = UN
      ENDDO
C
C
C
      IF (ICAS.EQ.1) THEN
        GOTO 1000
      ELSE
        IF (ICAS.GE.11) THEN
           IF (ICAS.EQ.11) DAM(1) = XXX(1)
           IF (ICAS.EQ.12) DAM(2) = XXX(2)
           IF (ICAS.EQ.13) THEN
             DAM(1) = XXX(1)
             DAM(2) = XXX(2)
           ENDIF
           GOTO 1000
        ELSE
          IF (ICAS.EQ.2) THEN
            IMEC(1) = 1
            IMEC(2) = 2
            GOTO 2000
          ELSE
            IF (ICAS.EQ.3) THEN
              IMEC(1) = 2
              IMEC(2) = 1
              GOTO 2000
            ELSE
              IF (ICAS.EQ.4) THEN
                GOTO 3000
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
C
C  Cas sans direction d'ecoulement
C
 1000 CONTINUE
       XKINV(1,1)=(UN/(YOUN*DAM(1)))
       XKINV(1,2)=((-1.)*XNU/YOUN)
       XKINV(2,1)=((-1.)*XNU/YOUN)
       XKINV(2,2)=(UN/(YOUN*DAM(2)))
      GOTO 5000
C
C  Cas avec 1 direction d'ecoulement
C
 2000 CONTINUE
      IF (EPS(IMEC(1)).LE.EPSR) THEN
        XKINV(IMEC(1),IMEC(1)) = UN/YOUN + UN/XHH1
        XKINV(1,2) = -XNU/YOUN
        XKINV(2,1) = -XNU/YOUN
        XKINV(IMEC(2),IMEC(2)) = UN/(DAM(IMEC(2))*YOUN)
        GOTO 5000
      ELSE
        XKTAP(IMEC(1),IMEC(1)) = XZER
        XKTAP(1,2) = XZER
        XKTAP(2,1) = XZER
        XKTAP(IMEC(2),IMEC(2)) = YOUN*DAM(IMEC(2))
        GOTO 9999
      ENDIF
C
C  Cas avec 2 directions d'ecoulement
C
C
 3000 CONTINUE
      IF ((EPS(1).LE.EPSR).AND.(EPS(2).LE.EPSR)) THEN
        XKINV(1,1) = UN/YOUN + UN/XHH1
        XKINV(1,2) = -XNU/YOUN
        XKINV(2,1) = -XNU/YOUN
        XKINV(2,2) = UN/YOUN + UN/XHH1
        GOTO 5000
      ELSE
        IF ((EPS(1).GT.EPSR).AND.(EPS(2).LE.EPSR)) THEN
           XKTAP(1,1) = XZER
           XKTAP(1,2) = XZER
           XKTAP(2,1) = XZER
           XKTAP(2,2) = UN/((UN/XHH1 + UN/YOUN))
           GOTO 9999
        ELSE
          IF ((EPS(1).LE.EPSR).AND.(EPS(2).GT.EPSR)) THEN
            XKTAP(1,1) = UN/((UN/XHH1 + UN/YOUN))
            XKTAP(1,2) = XZER
            XKTAP(2,1) = XZER
            XKTAP(2,2) = XZER
            GOTO 9999
          ELSE
            XKTAP(1,1) = XZER
            XKTAP(1,2) = XZER
            XKTAP(2,1) = XZER
            XKTAP(2,2) = XZER
            GOTO 9999
          ENDIF
        ENDIF
      ENDIF
C
C  On inverse la matrice si necessaire
C
 5000 CONTINUE
      XDET = XKINV(1,1)*XKINV(2,2) - XKINV(1,2)*XKINV(2,1)
      IF (ABS(XDET).GT.EPSILO) THEN
         XDET = UN/XDET
         XKTAP(1,1) = XDET*XKINV(2,2)
         XKTAP(2,2) = XDET*XKINV(1,1)
         XKTAP(1,2) = -XDET*XKINV(2,1)
         XKTAP(2,1) = -XDET*XKINV(1,2)
      ELSE
        KERRE = 2
        WRITE(6,*)'Impossible d inverser'
      ENDIF
C================================================================
C   PASSAGE DANS LE REPERE FIXE
C===============================================================
 9999  CONTINUE
      XKTAP(1,3) = XZER
      XKTAP(2,3) = XZER
      XKTAP(3,1) = XZER
      XKTAP(3,2) = XZER
      EPSF(1) = TP(1,1)*EPSN(1) + TP(1,2)*EPSN(2)
      EPSF(2) = TP(2,1)*EPSN(1) + TP(2,2)*EPSN(2)
      IF (ABS(EPSF(2)-EPSF(1)).GT.EPSILO) THEN
        XKTAP(3,3) = (SIGF(2)-SIGF(1))/(2.*(EPSF(2)-EPSF(1)))
      ELSE
        XKTAP(3,3) = XZER
      ENDIF
      CALL MATMAT(TM,XKTAP,3,3,3,XMINT)
      CALL MATMAT(XMINT,TP,3,3,3,XKTAT2)
      XKTAT(1,1)=XKTAT2(1,1)
      XKTAT(1,2)=XKTAT2(1,2)
      XKTAT(2,1)=XKTAT2(2,1)
      XKTAT(2,2)=XKTAT2(2,2)
      XKTAT(4,1)=XKTAT2(3,1)
      XKTAT(1,4)=XKTAT2(1,3)
      XKTAT(4,2)=XKTAT2(3,2)
      XKTAT(2,4)=XKTAT2(2,3)
      XKTAT(4,4)=XKTAT2(3,3)
      RETURN
      END