C FIBMAZ    SOURCE    MB234859  24/07/09    21:15:03     11959          
      SUBROUTINE FIBMAZ (XMAT,DEPS,SIG0,VAR0,SIGF,VARF)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XMAT(*),DEPS(3),SIG0(3),VAR0(5),SIGF(3),VARF(5)
C---------------------------------------------------------------
C     MODELE DE MAZARS (MAZARS_FIB) SPECIALISE POUR LE MODELE A FIBRE
C     Didier Combescure et Pierre Pegon Nov. 93
C---------------------------------------------------------------
C
C
C     variables en entree
C
C     XMAT : Caracteristique du materiaux
C
C     EPS : Deformations initiales
C
C     DEPS : Increment de deformations
C
C     SIG0 : Contraintes initiales
C
C     VAR0 : Variables internes initiales
C
C
C     variables en sortie
C
C     SIGF : Contraintes finales
C
C     VARF : Variables internes finales
C
C     declaration des variables
C
C
      REAL*8 EPS33(3,3),EPSIPP(3),EPSILC(3),VALP33(3,3)
      REAL*8 SIGP(3),SIGPC(3)
      PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)

C
C
C     recuperation des variables initiales dans les tableaux
C
C
      YOUN = XMAT(1)
      XNU  = XMAT(2)
      EPSD0= XMAT(5)
      ACOM = XMAT(6)
      BCOM = XMAT(7)
      ATRA = XMAT(8)
      BTRA = XMAT(9)
      BETA = XMAT(10)
      DINI = VAR0(2)
      EPSX = VAR0(3)
      EPSY = VAR0(4)
      EPSZ = VAR0(5)

C
C     calcul des seuils
C
C
C        calcul de la deformation totale
C
      DO 100 ISTRS=1,3
         VARF(ISTRS+2)=VAR0(ISTRS+2)+DEPS(ISTRS)
100   CONTINUE
      EPS33(1,1) = VARF(3)
      EPS33(1,2) = (VARF(4)/2.D0)
      EPS33(2,1) = EPS33(1,2)
      EPS33(1,3) = (VARF(5)/2.D0)
      EPS33(3,1) = EPS33(1,3)
      EPS33(2,2) = ((-1.D0*XNU)*VARF(3))
      EPS33(3,3) = EPS33(2,2)
      EPS33(3,2) = XZERO
      EPS33(2,3) = XZERO

C
C        calcul des deformations principales
C
C         on diagonalise
C
      CALL JACOB3 (EPS33,3,EPSIPP,VALP33)

C
C       on calcule les contraintes ppales
C
C
      SIGP(1)=((UN-XNU)*EPSIPP(1)+
     1         (XNU*EPSIPP(2))+
     2         (XNU*EPSIPP(3)))
      SIGP(2)=((UN-XNU)*EPSIPP(2)+
     1         (XNU*EPSIPP(1))+
     2         (XNU*EPSIPP(3)))
      SIGP(3)=((UN-XNU)*EPSIPP(3)+
     1         (XNU*EPSIPP(2))+
     2         (XNU*EPSIPP(1)))
      DO 200 ISTRS=1,3
         SIGP(ISTRS)=(SIGP(ISTRS)*YOUN/((UN+XNU)*(UN-2.D0*XNU)))
200   CONTINUE

C
C        on calcule le epsilontild
C
      EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
     1       MAX( XZERO , EPSIPP(2) )**2 +
     2       MAX( XZERO , EPSIPP(3) )**2
      EPSTIL=SQRT (EPSTIL)
C
C      on est toujours dans le cas istep=0
C
         EPSTIM=EPSTIL
         VARF(1)=EPSTIL
C
C      on calcule l'endommagement et les contraintes
C
C
C        on calcule ALFAT ALFAC DT et DC puis D
C        dans le cas ou le seuil initial est depasse
C
C        on calcule le signe des contraintes elastiques
C
            DO 300 ISTRS=1,3
                  SIGPC(ISTRS)=MIN(XZERO,SIGP(ISTRS))
300         CONTINUE
            TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
C
C       ainsi que GAMMA et EPSTIM si BETA est different de zero
C
CCC         IF (BETA.GT.1.D-5) THEN
CCC            IF (ABS(TRSIGC).GT.XPETIT) THEN
CCC              GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
CCC  1                                   SIGPC(3)*SIGPC(3)
CCC              GAMMA= SQRT(GAMMA)/ABS(TRSIGC)
CCC              EPSTIM=EPSTIM*GAMMA
CCC            END IF
CCC         END IF
C
         IF ( EPSTIM .GT. EPSD0) THEN
C
C      calcul de l'endommagement
C
C        on calcule les deformations dues aux contraintes negatives
C
            DO 400 ISTRS=1,3
               EPSILC(ISTRS)=(SIGPC(ISTRS)*(UN+XNU)-TRSIGC*XNU)/YOUN
400         CONTINUE
C
C        on en deduit ALFAT et ALFAC
C
            ALFAC = MAX(XZERO,EPSIPP(1))*EPSILC(1) +
     1              MAX(XZERO,EPSIPP(2))*EPSILC(2) +
     2              MAX(XZERO,EPSIPP(3))*EPSILC(3)
            ALFAC = ALFAC/(EPSTIL*EPSTIL)
            ALFAT = UN - ALFAC
C
C        amelioration de la reponse en cisaillement pour beta > 1.
C
            IF (BETA .GT. UN) THEN
               IF ( ALFAT .GT. XPETIT ) THEN
                  ALFAT=ALFAT**BETA
               END IF
               IF ( ALFAC .GT. XPETIT ) THEN
                  ALFAC=ALFAC**BETA
               END IF
            END IF

C
C        on calcule DT et DC puis D
C        dans le cas ou le seuil initial est depasse
C
               DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
     1            ATRA*EXP(BTRA*(EPSD0-EPSTIM))
               DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
     1            ACOM*EXP(BCOM*(EPSD0-EPSTIM))
               D = ALFAT*DT + ALFAC*DC
C
C           on borne la valeur de D a 0.99999
C
            D=MIN ( D , UN-1.D-05 )
         ELSE
            D=XZERO
         END IF
C
C           on teste la croissance de D
C
         D=MAX ( D , DINI )
C
C           on le stocke dans les variables internes finales
C
         VARF(2)= D
C
C           on calcule les contraintes finales
C
         SIGF(1) = YOUN*VARF(3)
         SIGF(2) = YOUN*VARF(4)/2.D0/(UN+XNU)
         SIGF(3) = YOUN*VARF(5)/2.D0/(UN+XNU)
         DO 500 ISTRS=1,3
            SIGF(ISTRS)=SIGF(ISTRS)*(UN-D)
500      CONTINUE

      RETURN
      END