Sources Esope | Cast3M

endothm
C ENDOTHM   SOURCE    FD218221  26/02/13    21:15:14     12462          
      SUBROUTINE ENDOTHM(EPSILO,NSTRS,NVARI,CMAT,DEP,SIGF,VAR0,
     *                   VARF,NMAT,IFOUR)
C
C_________________________________________________________________
C
C  CALCUL DE L ENDOMMAGEMENT ET DES CONTRAINTES APPARENTES FINALES
C  SOUS-PROGRAMME APPELE PAR CAMECTHM
C
C  Attention la regularisation par l'energie de fracture est prise
C  en compte (lignes 43-44)
C_________________________________________________________________
C
C     declaration des variables
C
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3)
      DIMENSION SIGP(3),SIGPT(3),SIGPC(3)
      DIMENSION EPSILO(NSTRS),SIGF(NSTRS),CMAT(NMAT)
      DIMENSION DEP(NSTRS,NSTRS),VARF(NVARI),VAR0(NVARI)
C
      XZERO=0.D0 
      UN=1.D0 
      XPETIT=1.D-12
C
C RECUPERATION DES PARAMETRES MATERIAUX
C         ET DES VARIABLES INTERNES
C
      YOUN = CMAT(1)
      XNU  = CMAT(2)
      EPSD0= CMAT(5)
      ATRA = CMAT(6)
      FRAE = CMAT(7)
      ACOM = CMAT(8)
      BCOM = CMAT(9)
      EQUL = CMAT(16)
      DINI = VAR0(2)
C
C REGULARIZATION CALCUL BTRA  (MODIFICA 01-11-14)
C        
      FTRA =  YOUN * EPSD0
      BTRA = (FTRA * EQUL) / (FRAE - ((FTRA**2.D0)*(EQUL/(2.D0*YOUN))))
C
C CHANGEMENT EPSILO
C  CONTRAINTES PLANES
C
      IF (IFOUR.EQ.-2) THEN
        temp = EPSILO(3)
        EPSILO(3) = EPSILO(4)
        EPSILO(4) = temp
      ENDIF
C
C REECRITURE DU VECTEUR DEFORMATION SOUS FORME MATRICIELLE
C   ENTREE : EPSILO -> VECTEUR DEFORMATION ELASTIQUE
C   SORTIE : EPS33 -> MATRICE DEFORMATION ELASTIQUE
C      Le 2 correspond au cisaillement
C
       CALL ENDOCB (EPSILO,EPS33,2,IFOUR)
C
C DIAGONALISATION DE LA MATRICE DEFORMATION ELASTIQUE EPS33
C    SORTIE : EPSIPP(3)   - VALEURS PROPRES
C             VALP33(3,3) - VECTEURS PROPRES
C
       IF (IFOUR.EQ.2) THEN
         CALL JACOB3 (EPS33,3,EPSIPP,VALP33)
       ELSE
         CALL JACOB3 (EPS33,2,EPSIPP,VALP33)
       ENDIF
C
C CALCUL DES CONTRAINTES PRINCIPALES
C
       DO 10 I=1,3
            SIGP(I)= XZERO
            DO 11 J=1,3
               SIGP(I)=SIGP(I)+DEP(I,J)*EPSIPP(J)
  11        CONTINUE
  10   CONTINUE
C
C ON COMPLETE LA DEFORMATION DANS LE CAS DES
C  CONTRAINTES PLANES 
C
      IF (IFOUR.EQ. -2) THEN
         EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU)
      ENDIF
C
C CALCUL DE EPSILON TILDA (DEFORMATION POSITIVE)
C
      EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
     *       MAX( XZERO , EPSIPP(2) )**2 +
     *       MAX( XZERO , EPSIPP(3) )**2
      EPSTIL=SQRT (EPSTIL)
C
C TEST EPSTIL>EPSD0
C
      IF ( EPSTIL .GT. EPSD0) THEN
C
C       CALCUL DE L ENDOMMAGEMENT
C
C        CALCUL DE LA TRACE DES CONTRAINTES
C          NEGATIVE (COMPRESSION) : TRSIGC
C          POSITIVE (TRACTION)    : TRSIGT
C
        DO 20 I=1,3
          IF (SIGP(I).LT. XZERO) THEN
            SIGPC(I) = SIGP(I)
            SIGPT(I) = XZERO
          ELSE
            SIGPT(I) = SIGP(I)
            SIGPC(I) = XZERO
          ENDIF
  20    CONTINUE
C
        TRSIGT = SIGPT(1)+SIGPT(2)+SIGPT(3)
        TRSIGC = SIGPC(1)+SIGPC(2)+SIGPC(3)
C
C     CALCUL DES DEFORMATIONS DUES AU CONTRAINTES POSITIVES
C
        DO 30 I=1,3
          EPSILT(I) = (SIGPT(I)*(UN+XNU)-TRSIGT*XNU)/YOUN
  30    CONTINUE
C
C     CALCUL DE ALFAT ET ALFAC
C
        ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
     *          MAX(XZERO,EPSIPP(2))*EPSILT(2) +
     *          MAX(XZERO,EPSIPP(3))*EPSILT(3)
        ALFAT = ALFAT/(EPSTIL*EPSTIL)
        ALFAC = UN - ALFAT
C
C     AMELIORATION DE LA REPONSE EN CISAILLEMENT
C        POUR BETA > 1
C
        BETA = UN
        IF (BETA .GT. UN) THEN
           IF ( ALFAT .GT. XPETIT ) THEN
             ALFAT=ALFAT**BETA
           ENDIF
           IF ( ALFAC .GT. XPETIT ) THEN
             ALFAC=ALFAC**BETA
           ENDIF
         ENDIF
C
c      CALCUL DE DT, DC ET D
C
         DT = UN - (EPSD0/EPSTIL)*(
     *             (UN + ATRA)*EXP(-BTRA*(EPSTIL-EPSD0)) -
     *             ATRA*EXP(-2.D0*BTRA*(EPSTIL-EPSD0)))
         DC = UN - (EPSD0/EPSTIL)*(
     *             (UN + ACOM)*EXP(-BCOM*(EPSTIL-EPSD0)) -
     *             ACOM*EXP(-2.D0*BCOM*(EPSTIL-EPSD0)))
         D = ALFAT*DT + ALFAC*DC
C
C      ON BORNE LA VALEUR DE D a 0.99999
C
         D = MIN ( D , UN-1.D-05 )
C         
      ELSE
         D = XZERO
      ENDIF 
C
C    ON TESTE LA CROISSANCE DE D 
C
      D=MAX ( D , DINI )
C
C STOCKAGE DES VARIABLES INTERNES FINALES
C 
      VARF(2) = D
C
C CONTRAINTES APPARENTES
C
      DO 40 I=1,NSTRS
        SIGF(I) = SIGF(I)*(1.D0-D)
   40 CONTINUE
C
C CHANGEMENT EPSILO
C  CONTRAINTES PLANES
C
      IF (IFOUR.EQ.-2) THEN
        temp = EPSILO(3)
        EPSILO(3) = EPSILO(4)
        EPSILO(4) = temp
      ENDIF
C
      END