C LABORD    SOURCE    CB215821  17/11/30    21:16:43     9639           
C     ==================================================================
C     =       SUBROUTINE INTEGRANT LE MODELE D'ENDOMMAGEMENT           =
C     =     UNILATERAL EN UNIAXIAL AVEC PILOTAGE EN DEFORMATION        =
C     =              FG + FR (LMT-Cachan, fevrier 2001)                =
C     ==================================================================
      SUBROUTINE LABORD(XMAT,DEPS,SIG0,VAR0,SIGF,VARF)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (a-h,o-z)

      REAL*8 XMAT(13),DEPS(3),SIG0(3),VAR0(8),SIGF(3),VARF(8)

C---------------------------------------------------------------
C     MODELE DE CLB (CLB_UNI) SPECIALISE POUR LE MODELE A FIBRE
C---------------------------------------------------------------
C
C     variables en entree
C
C     XMAT : Caracteristique du materiaux
C
C     EPS : Deformation totale
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
      logical change

C     ==================================================================
C     =     Initialisation des Constantes                              =
C     ==================================================================
      E=XMAT(1)
      XNU=XMAT(2)
      Y01=XMAT(5)
      Y02=XMAT(6)
      A1=XMAT(7)
      A2=XMAT(8)
      B1=XMAT(9)
      B2=XMAT(10)
      BETA1=XMAT(11)
      BETA2=XMAT(12)
      SIGF1=XMAT(13)
C     ==================================================================
C     =     Initialisation des VARIABLES HISTORIQUES                   =
C     ==================================================================
      ESEC = VAR0(1)
      D1 = VAR0(2)
      D2 = VAR0(3)
      YLIM1 = VAR0(4)
      YLIM2 = VAR0(5)
      EPS10 = VAR0(6)
C
      EPS = EPS10 + DEPS(1)
      y1=YLIM1-1.D0
      y2=YLIM2-1.D0
C     ==================================================================
***
      change=.true.
      do while (change)
***
      X1=(BETA1*D1/(1.D0-D1)+BETA2*D2/(1.D0-D2))/E
      X2=(BETA2*D2-SIGF1)/(1.D0-D2)/E
      if ( X1 .LE. Eps ) then
        EPSF=Eps-BETA2*D2/(1.0D0-D2)/E
        CALL CAS13(EPSF,SIG,ESEC,E,D1,BETA1,A1,B1,Y1,YLIM1,Y01)
        icas1=1
      else if ( Eps .LE. X2 ) then
        CALL CAS13(Eps,SIG,ESEC,E,D2,BETA2,A2,B2,Y2,YLIM2,Y02)
        icas1=3
      else
        call CAS2(Eps,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ESEC)
        icas1=2
      endif
      X1=(BETA1*D1/(1.D0-D1)+BETA2*D2/(1.D0-D2))/E
      X2=(BETA2*D2-SIGF1)/(1.D0-D2)/E
      if ( X1 .LE. Eps ) then
        icas2=1
      else if ( Eps .LE. X2 ) then
        icas2=3
      else
        icas2=2
      endif
      if (icas1 .eq. icas2) then
         change=.false.
      endif
***
      enddo
***
c     Si on garde le module du beton constant
      ESEC=E

C     ==================================================================
C     =     reactualisation des seuils                                 =
C     ==================================================================
      YLIM1= MAX(Y1,YLIM1)
      YLIM2= MAX(Y2,YLIM2)

C     ==================================================================
C     =     reactualisation des VARIABLES HISTORIQUES                  =
C     ==================================================================
      VARF(1) = ESEC
      VARF(2) = D1
      VARF(3) = D2
      VARF(4) = YLIM1
      VARF(5) = YLIM2
      VARF(6) = EPS
      VARF(7) = VAR0(7) + DEPS(2)
      VARF(8) = VAR0(8) + DEPS(3)
C     ==================================================================
C     =     Calcul des contraintes finales                             =
C     ==================================================================
      SIGF(1) = SIG
      r_z = 0.5D0 * E / (1.D0+XNU)
      SIGF(2) = r_z * VARF(7)
      SIGF(3) = r_z * VARF(8)

      RETURN
      END

C     ==================================================================
C     =               CALCUL DE SIGMA DANS LE CAS 1 & 3                =
C     =             ( le calcul se fait par dichotomie )               =
C     ==================================================================
      SUBROUTINE CAS13(EPS,SIG,ESEC,E,D,BETA,A,B,Y,YLIM,Y00)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

      REAL*8 m0,m

      dt = 1.0d0
      m0 = 0.0d0
      alpha = 0.0d0

c     Erreur relative par rapport a D
      Erreur_max=1.0D-6
      Nbpas=50

c     Dans le cas No 1 ou 3, on recherche D13 dans [0,1]
      X1= d-0.01D0
      X2= 1.0D0

      pas=(X2-X1)/Nbpas
      m=m0
      if ( EPS .ge. 0.D0) then
         m=m0/650.D0
      endif

C     ------------------------------------------------------------------
C     -   Avant de faire les calculs, on verifie que l'on depasse le   -
C     -   seuil YLIM                                                   -
C     ------------------------------------------------------------------
      D13=D
      y0=y00
      YY = (E*EPS+BETA)*(E*EPS+BETA)-BETA*BETA/(1.0D0-D13)/(1.0D0-D13)
      YY = YY/(2.0D0*E)

***   Prise en compte de la viscosite

****  coefficient de viscosite
c      m=3.d4
c      alpha=0.18d0
****
**      Si yy est inferieur a y0 on sort sans evolution de d
      If ( YY .LE. Y0 ) Then
        y=yy
c       Module Secant
        ESEC=E*(1.0D0-D)
c       Valeur de la contrainte
        SIG=EPS*ESEC-BETA*D
        return
      endif

***********
**      Calcul du seuil
c      y=(E*EPS+BETA)*(E*EPS+BETA)/(2.D0*E)

**      Calcul de l endommagement
      Erreur=1.D0
      XG=d
      Xd=X1
      Do while ( (Erreur .GT. Erreur_max) .and. (pas .ge. 1.d-8) )

c       Valeur de la fonction a la borne de droite D13=XD
        D13=XD
        YY = (E*EPS+BETA)*(E*EPS+BETA)/(2.0D0*E)
     &      -(BETA*BETA/(1.D0-D13)/(1.D0-D13))/(2.D0*E)
c        if (yy .lt. ylim) then
c          yy=ylim
c        endif
        y0=y00+m*( (abs(d13-d)+(d13-d))/(2.D0*dt) )**alpha
        If ( YY .GT. Y0 ) Then
**      si Dconverge est > a D13 alors fctd est < 0
          Fctd =  1.0D0 - (1.0D0+(A*(YY-Y0))**B)*(1.0D0-D13)
        Else
          Fctd =  1.d0
        Endif

        If ( Fctd .lt. 0.D0) then
          Xg=Xd
          Xd=Xd+pas
        Else
          Xd=Xg
          Pas=Pas * 0.5D0
          Xd=Xd+pas
        Endif

        If ( Xg .gt. X2 ) Then
          Stop ' . DANS CAS13 . Pas de Zero.'
        Endif

        Erreur=abs(Fctd)
      Enddo

c     On verifie que d13 > d . Si ce n'est pas le cas rien n'evolue.
c      print*, (abs(d13-d)+(d13-d))/(2.D0*dt)
***
      y=y00+1.D0/a*(d13/(1.D0-d13))**(1.D0/b)
***
      If ( d13 .Ge. d ) then
        D = D13
      Endif
      if (d .gt. 0.9999999D0) then
         d=0.9999999D0
      endif

c     Module Secant
      ESEC=E*(1.0D0-D)
c     Valeur de la contrainte
      SIG=EPS*ESEC-BETA*D

      Return
      End

C     ==================================================================
C     =                 CALCUL DE SIGMA DANS LE CAS 2                  =
C     =             ( le calcul se fait par dichotomie )               =
C     ==================================================================
      SUBROUTINE CAS2(EPS,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ESEC)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

      fprim=  (e*eps*(1.D0-d2)-beta2*d2+sigf1)
     &      / (sigf1+beta1*d1*(1.D0-d2)/(1.D0-d1))

c     La solution est la borne de droite.
      SIG=e*eps*(1.D0-d2)-beta2*d2-beta1*d1*(1.D0-d2)*fprim/(1.D0-d1)

c     ESEC est le module du beton de compression
      ESEC= E*(1.0D0-D2)

      Return
      End

C=======================================================================
C REMARQUE : LES SOUS-PROGRAMMES RAIDTANi (avec i=1,2,3) NE SONT
C            PAS UTILISES !
C=======================================================================
      SUBROUTINE RAIDTAN1 (EPS,SIG,ETAN,E,D1,BETA1,D2,BETA2,A,B,Y,Y0)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

**      derive de d par rapport a y
        dddy= b*(a**b)*(y-y0)**(b-1.D0)/(1.D0+(a*(y-y0))**b)**2

**      derive de d par rapport a epsilon
        dddeps=dddy*( e*eps+beta1-beta2*d2/(1.D0-d2) )
        dddeps=dddeps/( 1.D0+dddy*beta1**2/(e*(1.D0-d1)**3) )

**      derive de sigma par rapport a epsilon -- calcul de E tangent
        etan=e*(1.D0-d1)-e*eps*dddeps-beta1*dddeps+beta2*d2*dddeps/
     $  (1.D0-d2)

        return
        end
C---------
      SUBROUTINE RAIDTAN2 (EPS,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ETAN)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

        etan=E*(1.D0-d2)/( sigf1+beta1*d1*(1.D0-d2)/(1.D0-d1) )
        etan=etan*beta1*d1*(1.D0-d2)/(1.D0-d1)
        etan=e*(1.D0-d2)-etan

        return
        end
C---------
      SUBROUTINE RAIDTAN3 (EPS,SIG,ETAN,E,D1,BETA1,D2,BETA2,A,B,Y,Y0)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

**      derive de d par rapport a y
        dddy= b*(a**b)*(y-y0)**(b-1.D0)/(1.D0+(a*(y-y0))**b)**2

**      derive de d par rapport a epsilon
        dddeps=dddy*(e*eps+beta2)
        dddeps=dddeps/(1.D0+dddy*beta2**2/(e*(1.D0-d2)**3))

**      derive de sigma par rapport a epsilon -- calcul de E tangent
        etan=e*(1.D0-d2)-e*eps*dddeps-beta2*dddeps

        return
        end

C     ==================================================================
C     =                 CALCUL DE SIGMA DANS LE CAS 2                  =
C     =             ( le calcul se fait par dichotomie )               =
C     ==================================================================
      SUBROUTINE CAS21(EPS,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ESEC)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

c     Erreur relative par rapport a Sigma
      Erreur_max=1.0D-6
      Nbpas=50
      itermax=1000
      iter=0
c     Dans le cas No 2, on recherche Y dans [-1,0]   Y=sig/sigf1
      X1= -1.0D0
      X2=  0.0D0

      pas=(X2-X1)/Nbpas

c     Valeur initiale des bornes de gauche et droite
      XG=X1
      XD=XG+PAS

c     Valeur de la fonction a la borne de gauche y=XG
      Y=XG
      Fctg    =  EPS*E - BETA1*D1*(1.0D0-3.0D0*Y*Y-2.0D0*Y*Y*Y)/
     $  (1.0D0-D1)
     &        - BETA2*D2/(1.0D0-D2) - Y*SIGF1/(1.0D0-D1)
     &        - SIGF1*(1.0D0/(1.0D0-D1)-1.0D0/(1.0D0-D2))*(2.0D0+Y)*Y*Y
      If (Fctg .GE. 0.0D0) then
        isigg=1
      else
        isigg=0
      endif
      Erreur=abs(Fctg)
C
      Do while ( (Erreur .GT. Erreur_max) .and. (iter .le. itermax) )
         iter=iter+1
c       Valeur de la fonction a la borne de droite y=XD
        Y=XD
        Fctd    =  EPS*E - BETA1*D1*(1.0D0-3.0D0*Y*Y-2.0D0*Y*Y*Y)/
     &   (1.0D0-D1)    - BETA2*D2/(1.0D0-D2) - Y*SIGF1/(1.0D0-D1)
     &       - SIGF1*(1.0D0/(1.0D0-D1)-1.0D0/(1.0D0-D2))*(2.0D0+Y)*Y*Y
        If (Fctd .GE. 0.0D0) then
          isigd=1
        else
          isigd=0
        endif
        Erreur=abs(Fctd)
        If ( isigg .EQ. isigd) then
          Fctg=Fctd
          isigg=isigd
          Xg=Xd
          Xd=Xd+pas
        Else
          Pas=Pas * 0.5D0
          Xd=Xg+pas
        Endif
        If ( Xg .gt. X2 ) Then
          Stop ' . DANS CAS2 . Pas de Zero.'
        Endif
      Enddo
c     La solution est la borne de droite.
      SIG=Y*Sigf1
c     Si on fait varier lineairement ESEC en fonction de Sig
      ESEC= E*(1.0D0 - D1 + Y*(D2-D1))
      Return
      End