Sources Esope | Cast3M

ouglof
C OUGLOF    SOURCE    BR232186  16/09/12    12:43:41     9078           
      SUBROUTINE OUGLOF(XMAT,DEPS,SIG0,VAR0,SIGF,VARF)
 
C
 
C====&===1=========2=========3=========4=========5=========6=========7==
 
C Commentaires : Subroutine permettant de mettre en oeuvre le
 
C                modele OUGLOVA pour representer
 
C                le comportement de l'acier corrodé
 
C
 
C
 
C Auteurs      : R. Paredes
 
C              : CEA-DEN/DANS/DM2S/SEMT/EMSI
 
C              : Romili.Paredes@cea.fr
 
C
 
C Date         : Aout 2016
 
C====&===1=========2=========3=========4=========5=========6=========7==
 
C
 
C-----DECLARATION GENERALE----------------------------------------------
 
C
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
C
 
C-----DECLARATION DES VARIABLES-----------------------------------------
 
C
 
      DIMENSION XMAT(*),DEPS(*),SIG0(*),VAR0(*),SIGF(*),VARF(*)
 
 
 
C     Parametre pour le nombre d iteration locales internes
 
      IMAX = 1000
 
C
 
C-----DONNEES MATERIAUX-------------------------------------------------
 
C
 
      XYG   = XMAT(1)
 
      XNU   = XMAT(2)
 
      XGC   = XYG/(2.0D0*(XNU + 1.0D0))
 
      XSIGY = XMAT(5)
 
      XK    = XMAT(6)
 
      XM    = XMAT(7)
 
      XTC   = XMAT(8)
 
      XDC   = XMAT(9)
 
C
 
C-----SEUIL DE CORROSION------------------------------------------------
 
C
 
C-----Deformation plastique a rupture
 
      IF (XTC.LE.15.0D0) THEN
 
         EPSR = -0.0111D0*XTC + 0.2345D0
 
      ELSE 
 
         EPSR = -0.0006D0*XTC + 0.051D0
 
      ENDIF
 
 
 
C-----Seuil de rupture et d'endommagement
 
      XPR = EPSR
 
      XPD = 0.8D0*XPR
 
C
 
C-----ETAT INITIAL------------------------------------------------------
 
C
 
      XD    = VAR0(1)
 
      XR    = VAR0(2)
 
      XP    = VAR0(3)
 
      XZT   = VAR0(4)
 
      XNRUP = VAR0(5)
 
      XEPSP = VAR0(6)
 
      XEPSF = DEPS(1) + VAR0(10)
 
C
 
C-----CALCUL------------------------------------------------------------
 
C      
 
 
 
C-----Prediction elastique
 
      SIG    = (1.0D0 - XD)*XYG*(XEPSF - XEPSP)
 
      SIGEQ  = SIG
 
 
 
C-----Boucle d'endommagement      
 
      DO I = 1,IMAX
 
C--------Evaluation de la fonction seuil
 
         FP  = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY)
 
         FP0 = FP
 
 
 
C--------Boucle de plasticite          
 
         IF (FP.GT.1.0D0) THEN
 
 
 
            DO J = 1,IMAX
 
C--------------Calcul des derives
 
               DFDSG  = SIG/(SIGEQ*(1.0D0 - XD))
 
               DPHDSG = DFDSG
 
 
 
C--------------Calcul du multiplicateur plastique
 
               TAU     = XYG*DPHDSG
 
               AUX1    = DFDSG*TAU
 
               AUX2    = SIGEQ/(XK*(1.0D0 - XD))
 
               AUX3    = XSIGY/XK
 
               AUX4    = 1.0D0 - XM
 
               AUX5    = XK/XM
 
               XLAMBDA = FP/(AUX1 + AUX5*((AUX2 - AUX3)**AUX4))
 
 
 
C--------------Mise à jour des variables
 
               SIG    = SIG - XLAMBDA*TAU
 
               SIGEQ  = SIG
 
               XEPSP  = XEPSP + XLAMBDA*DPHDSG
 
               XP     = XEPSP
 
               XR     = XK*(XP**(1.0D0/XM))
 
 
 
C--------------Reevaluation de la fonction seuil
 
               FP  = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY)
 
 
 
               CRIT1 = DABS(FP/FP0)
 
               CRIT2 = DABS(XLAMBDA)
 
 
 
               IF ((CRIT1.LE.1.0D-8).OR.(CRIT2.LE.1.0D-10)) THEN 
 
                   GOTO 910
 
               ENDIF
 
 
 
            ENDDO
 
 
 
C-----Sortie de la boucle de plasticité
 
 910  CONTINUE
 
 
 
C-----------Evaluation du seuil d'endommagement
 
            FD = XP - (XPD + XZT)
 
            IF (I.EQ.1) THEN
 
                FD0 = FD
 
            ENDIF
 
 
 
            CRIT = DABS(FD/FD0)
 
 
 
            IF ((FD.GT.0.0D0).AND.(CRIT.GT.1.0D-5)) THEN
 
 
 
C--------------Calcul de l'endommagement
 
               XD  = (XDC/(XPR - XPD))*(XP - XPD)
 
               XZT = XP - XPD
 
 
 
C--------------Endommagement critique                
 
               IF (XD.GT.XDC) THEN
                  XD = XDC
 
               ENDIF
 
 
 
            ELSE
 
                GOTO 900
 
            ENDIF
 
 
 
            IF (CRIT.LE.1.0D-5) THEN
 
                GOTO 900
 
            ENDIF
 
 
 
         ELSE
 
            GOTO 900
 
         ENDIF
 
 
 
      ENDDO
 
 
 
C-----Sortie de la boucle d'endommagement
 
 900  CONTINUE
 
 
 
C-----Rupture
 
      IF (XP.GE.XPR) THEN
 
         XNRUP = 1
 
         XP  = XPR
 
         SIG = 1.0D0
 
      ENDIF
 
C
 
C-----SOCKAGE EN SORTIE-------------------------------------------------
 
C
 
C-----Les variables internes
 
      VARF(1)  = XD
 
      VARF(2)  = XR
 
      VARF(3)  = XP
 
      VARF(4)  = XZT
 
      VARF(5)  = XNRUP
 
      VARF(6)  = XEPSP
 
      VARF(10) = XEPSF
 
 
 
C-----Les contraintes
 
      SIGF(1) = SIG
 
      SIGF(2) = SIG0(2) + XGC*DEPS(2)
 
      SIGF(3) = SIG0(3) + XGC*DEPS(3)
 
 
 
C-----Fin de l integration
 
      RETURN
 
      END