Sources Esope | Cast3M

ricbetf1
C RICBETF1  SOURCE    CB215821  16/04/21    21:18:16     8920
      SUBROUTINE RICBETF1(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 RICBET pour representer
C                le comportement simplifie d'un beton ordinaire sous
C                sollicitations cycliques
C
C Traits       : - Endommagement scalaire
C                - Glissement frottant
C                - multifibre
C                - Monotone et cyclique
C                - Effet unilatteral
C                - Comportement non lineaire en compression
C
C Auteurs      : B. Richard
C              : CEA-DEN/DANS/DM2S/SEMT/EMSI
C              : Benjamin.Richard@cea.fr
C
C Date         : Mars 2011
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(15),DEPS(3),SIG0(3),VAR0(10),SIGF(3),VARF(10)
      REAL*8 VAR1(10),DEPS1(3)
 
C     Parametre pour le nombre d iteration locales internes
      IMAX = 500
 
C     Parametre pour le nombre de sous pas lie au sub-steeping
      NENTIER = 50
 
C
C-----DONNEES MATERIAUX-------------------------------------------------
C
      XYG  = XMAT(1)
      XGC  = XMAT(1)/2.0D0/(XMAT(2)+1.0D0)
      XNU  = XMAT(2)
      XY0  = 0.5D0*XYG*(XMAT(5)/XMAT(1))**2.0D0
 
      XAT  = XMAT(6)
      XGAM = XMAT(7)
      XAA  = XMAT(8)
      XSIGF= XMAT(9)
 
      XFC = XMAT(10)
      XAF = XMAT(11)
      XAG = XMAT(12)
      XAC = XMAT(13)
      XBC = XMAT(14)
      XSIGU = XMAT(15)
 
      SXX = SIG0(1)
C
C-----PREPARATION SUBSTEPPING-------------------------------------------
C
      DO I=1,10
 
         VAR1(I) = VAR0(I)
 
      ENDDO
 
      DO I=1,3
 
         DEPS1(I) = DEPS(I)/NENTIER
 
      ENDDO
C
C-----LANCEMENT DU SUBSTEPPING------------------------------------------
C
      DO IQ=1,NENTIER
C
C-----VARIABLES INTERNES------------------------------------------------
C
      EPSXX = VAR1(1)+DEPS1(1)
      EPSXY = (VAR1(2)+DEPS1(2))/2.0D0
      EPSXZ = (VAR1(3)+DEPS1(3))/2.0D0
 
C
      XD = VAR1(4)
      XZ = VAR1(5)
      EPSPI = VAR1(6)
      XEC = VAR1(7)
C
      RECO = VAR1(8)
      DEFN = VAR1(9)
      DEFP = VAR1(10)
 
      IF (SXX.GE.0.0D0) THEN
C
C-----ENDOMMAGEMENT-----------------------------------------------------
C
 
C     Partie positive de la deformation
      EPSXXP = 0.5D0*(ABS(EPSXX-DEFN)+(EPSXX-DEFN))
 
C     Seuil en energie
      XYE=0.5D0*XYG*EPSXXP*EPSXXP
 
C     Calcul du seuil
      XENDO = XYE - (XY0 + XZ)
 
      IF ((XENDO.GT.0.0D0).AND.(SXX.GE.0.0D0)) THEN
         XD=1.0D0-1.0D0/(1.0D0+XAT*(XYE-XY0))
         XZ=XYE-XY0
      ENDIF
C
C-----GLISSEMENT INTERNE------------------------------------------------
C
 
C     La contrainte de frottement
      SPI = XD*XYG*((EPSXX-DEFN)-EPSPI)
 
C     Verification du seuil de glissement interne
C      IF ((ABS(SPI-XEC).GT.0.0D0).AND.(SXX.GE.0.0D0)) THEN
 
      IF (ABS(SPI-XEC).GT.0.0D0) THEN
C       Valeur initiale du seuil pour critere relatif
        XF0=ABS(SPI-XEC)
 
C       Initialisation du multiplicateur de Lagrange
        XDLAM=1.0D0
 
C       Debut des iterations internes
        DO I=1,IMAX
 
C          Calcul du critere relatif
           SEUIL=ABS(SPI-XEC)/XF0
 
C          Calcul de la derive du critere
           XSG=ABS(SPI-XEC)/(SPI-XEC)
 
C          Calcul du multiplicateur de Lagrange
           XDLAM=ABS(SPI-XEC)/(XYG*XD-XSG*XGAM*(-XSG+XAA*XEC))
 
C          Actualisation des variables internes
           XEC=XEC-XGAM*XDLAM*(-XSG+XAA*XEC)
 
           SPI=SPI-XD*XYG*XDLAM*XSG
 
C          Calcul de la deformation de glissement
           EPSPI=-SPI/(XD*XYG)+(EPSXX-DEFN)
 
C          Recalcul du seuil pour eventuel arret
           SEUIL=ABS(SPI-XEC)/XF0
 
C          Critere d arret
           IF ((SEUIL.LE.1.0D-8).OR.(XDLAM.LE.1.0D-10)) THEN
 
              GOTO 900
 
           ENDIF
 
        ENDDO
 
      ENDIF
 
C     Balise de sortie
 900  CONTINUE
 
      SXX=XYG*(EPSXX-DEFN-XD*EPSPI)
 
      GOTO 761
 
      ELSEIF (SXX.LE.XSIGF) THEN
 
C        Calcul de la contrainte effective
         SXX=XYG*(EPSXX-DEFN)
 
C        Calcul du seuil plastique
         XFP = XAF*ABS(SXX) - (RECO + XFC)
 
C        On test le seuil plastique
         IF (XFP.GT.0.0D0) THEN
 
C           On retient le seuil pour le critere relatif
            XFP0 = XFP
 
C           Debut des iterations internes
            DO NT=1,IMAX
 
               XFP = XAF*ABS(SXX) - (RECO + XFC)
 
               DXFP = XAF*ABS(SXX)/SXX
 
               DXGP = XAG*ABS(SXX)/SXX
 
               XTEMP1 = XYG*DXGP
 
               XTEMP2 = EXP(-XBC*DEFP)*(XAC-(XAC*DEFP+XFC+XSIGU)*XBC)
 
               XLAMP = XFP/(XTEMP1*DXFP+XTEMP2)
 
               RECO = RECO + XTEMP2*XLAMP
 
               DEFP = DEFP + XLAMP
 
               DEFN = DEFN + XLAMP*DXGP
 
               SXX = SXX - XLAMP*XTEMP1
 
               XFP = XAF*ABS(SXX) - (RECO + XFC)
 
             IF ((ABS(XFP/XFP0).LE.1.0D-8).OR.(XLAMP.LE.1.0D-10)) THEN
 
                  GOTO 666
 
               ENDIF
 
            ENDDO
 
         ENDIF
 
 666     CONTINUE
         GOTO 761
 
      ELSE
 
         SXX = SXX + DEPS1(1)*XYG/
     &         (1.0D0-(XYG*XD*EPSPI)/(XSIGF))
         GOTO 761
 
      ENDIF
 
 761  CONTINUE
 
C
C-----MISE A JOUR DE VAR1-----------------------------------------------
C
      VAR1(1 ) = EPSXX
      VAR1(2 ) = 2.0D0*EPSXY
      VAR1(3 ) = 2.0D0*EPSXZ
      VAR1(4 ) = XD
      VAR1(5 ) = XZ
      VAR1(6 ) = EPSPI
      VAR1(7 ) = XEC
      VAR1(8 ) = RECO
      VAR1(9 ) = DEFN
      VAR1(10) = DEFP
C
C-----FIN DU SUBSTEPPING------------------------------------------------
C
      ENDDO
C
C-----SOCKAGE EN SORTIE-------------------------------------------------
C
C     Les variables internes
      VARF(1) = EPSXX
      VARF(2) = 2.0D0*EPSXY
      VARF(3) = 2.0D0*EPSXZ
      VARF(4) = XD
      VARF(5) = XZ
      VARF(6) = EPSPI
      VARF(7) = XEC
      VARF(8) = RECO
      VARF(9) = DEFN
      VARF(10) = DEFP
 
C     Les contraintes
      SIGF(1) = SXX
      SIGF(2) = SIG0(2)+XGC*DEPS(2)
      SIGF(3) = SIG0(3)+XGC*DEPS(3)
 
C     Fin de l integration
      RETURN
      END