C OUGLFI    SOURCE    CB215821  16/04/21    21:17:55     8920
C sub ouglofibc
C====&===1=========2=========3=========4=========5=========6=========7==
C Commentaires : Subroutine permettant de mettre en oeuvre le
C                    modele d'OUGLOVA d'acier corrodé
C Traits       : - Endommagement isotrope couplé à la plasticité
C                - Monotone et cyclique
C Auteurs      : B. Richard
C Date         : Février 2008
C====&===1=========2=========3=========4=========5=========6=========7==
      SUBROUTINE OUGLFI(XMAT,DEPST,SIG0,VAR0,SIGA,VARF)

      REAL*8 XMAT(20),VAR0(25),VARF(25),DEPST(3)
        REAL*8 SIG0(3),SIGA(3)
        REAL*8 ES,NUS,DCS,SIGMAES,TCS,MS,KS
        REAL*8 ENDS,PVAR,DEPP,RVAR,SIGS,EPSP
        REAL*8 SIGMA,XLAM,EPTO,EPSR,PR,PD
        REAL*8 CRITPLAS,J2,CRITPLAS0,CRITP,DELTAP,CRITENDO
        REAL*8 TAU,DFDSIG,DGDSIG,DFDR

        INTEGER INDITER,ITER,FLAG,I1,ITERMAX,I2

C
C OPTIONS INTERNES
C
        ITERMAX = 100000
C  ---------------------------
C |     ACIER CORRODE         |
C  ---------------------------
C      XMAT( 1) ES
C      XMAT( 2) NUS
C      XMAT( 3)
C      XMAT( 4)
C      XMAT( 5)
C      XMAT( 6) DCS
C      XMAT( 7) SIGMAES
C      XMAT( 8) TCS
C      XMAT( 9) MS
C      XMAT(10) KS
C
C  ---------------------------
C |     ACIER CORRODE         |
C  ---------------------------
C      VAR0( 1) VARF( 1) ENDS
C      VAR0( 2) VARF( 2) PVAR
C      VAR0( 3) VARF( 3) DEPP
C      VAR0( 4) VARF( 4) RVAR
C      VAR0( 5) VARF( 5) SIGS
C      VAR0( 6) VARF( 6) EPSP
C      VAR0( 7) VARF( 7)
C      VAR0( 8) VARF( 8)
C      VAR0( 9) VARF( 9)
C      VAR0(10) VARF(10)
C
C  ---------------------------
C |          AUTRES           |
C  ---------------------------
C      VAR0(21) VARF(21) XLAM
C      VAR0(22) VARF(22) EPTO
C      VAR0(23) VARF(23) INDITER

C
C LES PARAMETRES MATERIAUX
C
        ES      = XMAT( 1)
      NUS     = XMAT( 2)
      DCS     = XMAT( 6)
      SIGMAES = XMAT( 7)
      TCS     = XMAT( 8)
      MS      = XMAT( 9)
      KS      = XMAT(10)
C
C LES VARIABLES INTERNES
C
      ENDS = VAR0( 1)
      PVAR = VAR0( 2)
      DEPP = VAR0( 3)
      RVAR = VAR0( 4)
      SIGS = VAR0( 5)
      EPSP = VAR0( 6)

        XLAM    = 1.0D0

        EPTO    = VAR0(22)
        INDITER = VAR0(23)
C
C AVANCE EN DEFORMATION
C
        EPTO = (EPTO + DEPST(1))*XLAM
C
C CALCUL DE QUELQUES GRANDEURS
C
        IF (TCS.LE.15.0D0) THEN
                EPSR = -0.0111D0*TCS+0.2345D0
        ELSE
                EPSR = -0.0006D0*TCS+0.0510D0
        END IF

        PR = EPSR

        PD = 0.8D0*PR
C
C PREDICTION ELASTIQUE
C
        SIGMA = (1.0D0-ENDS)*ES*(EPTO-EPSP)
C
C CRITERE DE PLASTICITE ET INITIALISATIONS
C
        J2        = ABS(SIGMA)
        CRITPLAS  = J2/(1.0D0-ENDS)-(RVAR+SIGMAES)
        CRITPLAS0 = CRITPLAS
        CRITP     = 1.0D0
        DELTAP    = 20.0D0
        ITER      = 0

        IF ((CRITPLAS.LE.0.0D0).OR.((DEPST(1)*EPTO).LE.0.0D0)) THEN
                FLAG = 1
                DEPP = 0.0D0
                GOTO 30
        ELSE IF (CRITPLAS.GT.0.0D0) THEN
                FLAG = 0
                GOTO 10
        END IF
C
C DEBUT DES ITERATIONS INTERNES
C
 10     DO I1=1,ITERMAX
C
C ON VERIFIE LE CRITERE
C

                DO I2=1,ITERMAX
                        IF ((CRITP.GT.1.0D-5).AND.
     &            ((DELTAP.GT.1.0D-15).OR.(CRITPLAS.GT.0.0D0)).AND.
     &                  (ITER.LT.50)) THEN

                               ITER = ITER +1

                               DFDSIG = (ABS(SIGMA)/SIGMA)/(1.0D0-ENDS)
                               DGDSIG = DFDSIG
                               DFDR   = -1.0D0

                               TAU    = DGDSIG*ES
                               DELTAP = CRITPLAS/(DFDSIG*TAU-
     &                     DFDR*KS/MS*(J2/(KS*(1.0D0-ENDS))-
     &                                   SIGMAES/KS)**(1.0D0-MS))

                                IF (DELTAP.LT.0.0D0) THEN
                                        DELTAP = 0.0D0
                                END IF
C
C ACTUALISATION
C
                        SIGMA = SIGMA - DELTAP*TAU
                        EPSP  = EPSP  + DELTAP*DGDSIG
                        DEPP  = DELTAP*DGDSIG
                        PVAR  = PVAR  + DELTAP/(1.0D0-ENDS)
                        RVAR  = KS*(PVAR)**(1.0D0/MS)
C
C ACTUALISATION DU SEUIL
C
                        J2        = ABS(SIGMA)
                        CRITPLAS  = J2/(1.0D0-ENDS)-(RVAR+SIGMAES)
                        CRITP     = CRITPLAS/CRITPLAS0

                        ELSE

                                GOTO 15

                        END IF
                END DO
C
C BALISE DE SORTIE DE LA SECONDE BOUCLE D ITERATIONS INTERNES
C
 15   CONTINUE
                RVAR = J2/(1.0D0-ENDS)-SIGMAES

                GOTO 20

 20   CONTINUE
                CRITENDO = PVAR - PD

                IF (CRITENDO.LT.0.0D0) THEN
                        FLAG = 1

                        GOTO 30

                ELSE
                        ENDS = DCS/(PR-PD)*(PVAR-PD)

                  FLAG = 1

                        IF (ENDS.GT.0.99D0) THEN
                                ENDS = 0.99D0
                        END IF

                        GOTO 30

           END IF

      END DO
C
C BALISE DE SORTIE FINALE
C
 30     CONTINUE

        IF (PVAR.GT.PR) THEN
                PVAR  = PR
            SIGMA = 0.0D0
        ELSE
                SIGMA = (1.0D0-ENDS)*ES*(EPTO-EPSP)
        END IF
C
C LES SORTIES
C
        VARF( 1) = ENDS
      VARF( 2) = PVAR
      VARF( 3) = DEPP
      VARF( 4) = RVAR
      VARF( 5) = SIGMA
      VARF( 6) = EPSP
      VARF( 9) = EPTO

      VARF(22) = EPTO
      VARF(23) = INDITER

C
C LES CONTRAINTES
C
        SIGA(1)  = SIGMA
        SIGA(2)  = SIG0(2)+(ES/(2.0D0*(1.0D0+NUS)))*DEPST(2)*XLAM
        SIGA(3)  = SIG0(3)+(ES/(2.0D0*(1.0D0+NUS)))*DEPST(3)*XLAM

        RETURN
        END