Sources Esope | Cast3M

lcgldm
C LCGLDM    SOURCE    CB215821  16/04/21    21:17:40     8920
      SUBROUTINE LCGLDM (XMAT,VAR0,VARF,SIG0,SIGF,DEPST,XCAR)
C
C====&===1=========2=========3=========4=========5=========6=========7==
C Commentaires :   Subroutine permettant de mettre en oeuvre le
C                  modele LGCR_DM (EDF R&D/AMA) pour representer
C                  le comportement cyclique du béton armé sous
C                  sollicitations cycliques
C
C Traits       : - Endommagement scalaire
C                - 2D (DKT/COQUES/PLAQUES MINCES) pré-intégrée
C                - Monotone et cyclique
C                - Effet unilateral (retour à l'origine)
C
C Auteur       : B. Richard (Dr - Ing.) - CEA/DEN/DANS/DM2S/SEMT/EMSI
C====&===1=========2=========3=========4=========5=========6=========7==
C
C----DECLARATION GENERALES----------------------------------------------
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
C----DECLARATION SEGMENTS-----------------------------------------------
C
      REAL*8 EPSIM(3,3),EPSIF(3,3),UNIT(2,2),SIGMAM(3,3)
      REAL*8 SIGMAF(3,3),VAPM(3),VEPM(3,3),VAPF(3),VEPF(3,3)
      REAL*8 XMAT(*),SIG0(*),SIGF(*),VAR0(*),VARF(*),DEPST(*),XCAR(*)
      REAL*8 XLAMM,XLAMF,XMUM,XMUF,XD1,XD2,XTOL
      REAL*8 XTRAV1(3,3),XTRAV2(3,3),VEPMT(3,3),VEPFT(3,3)
      REAL*8 XTRAV3(3,3),XTRAV4(3,3),XIMI(2),XIFI(2),BETAF(3,3)
      REAL*8 XTRAV5(2,2),XTRAV6(2,2),BETAM(3,3),SGEFF(3)
C
C----DECLARATION PARAMETRES---------------------------------------------
C
      PARAMETER (UN=1.D0)
      PARAMETER (XTOL=1.D-8)
 
      CALL ZERO(EPSIM,3,3)
      CALL ZERO(EPSIF,3,3)
 
      DO I=1,2
          DO J=1,2
              IF (I.EQ.J) THEN
                  UNIT(I,J) = 1.0D0
              ELSE
                  UNIT(I,J) = 0.0D0
              ENDIF
          ENDDO
      ENDDO
C
C----MISE EN DONNEES----------------------------------------------------
C
      YOUN = XMAT(1)
      XNU  = XMAT(2)
 
      XYOUNM = XMAT(1)
      XNUM  = XMAT(2)
      XYOUNF = XMAT(5)
      XNUF  = XMAT(6)
 
      EPAI = XCAR(1)
 
      XLAMM = XNUM*XYOUNM*EPAI/((1.0D0+XNUM)*(1.0D0-2.0D0*XNUM))
      XMUM  = XYOUNM*EPAI/(2.0D0*(1.0D0+XNUM))
 
      XLAMF = (XNUF*XYOUNF*(EPAI**(3.0D0)))/
     &        (12.0D0*(1.0D0-(XNUF**(2.0D0))))
      XMUF  = (XYOUNF*(EPAI**(3.0D0)))/
     &        (24.0D0*(1.0D0+XNUF))
 
      XGAMMT = XMAT(7)
      XGAMMC = XMAT(8)
      XGAMF  = XMAT(9)
      XSEUI0 = XMAT(10)
      XALFA  = XMAT(11)
C
C----VARIABLES INTERNES-------------------------------------------------
C
      EPSIM(1,1) = VAR0(1)
      EPSIM(2,2) = VAR0(2)
      EPSIM(1,2) = VAR0(3)/2.0D0
      EPSIM(2,1) = VAR0(3)/2.0D0
 
      EPSIF(1,1) = VAR0(4)
      EPSIF(2,2) = VAR0(5)
      EPSIF(1,2) = VAR0(6)/2.0D0
      EPSIF(2,1) = VAR0(6)/2.0D0
 
      XD1 = VAR0(9)
      XD2 = VAR0(10)
 
      XD10 = VAR0(9)
      XD20 = VAR0(10)
 
C         print*,'DEPST',(DEPST(I),I=1,3)
 
C
C----ACTUALISATION DEFORMATION------------------------------------------
C
      EPSIM(1,1) = EPSIM(1,1) + DEPST(1)
      EPSIM(2,2) = EPSIM(2,2) + DEPST(2)
      EPSIM(1,2) = EPSIM(1,2) + DEPST(3)/2.0D0
      EPSIM(2,1) = EPSIM(2,1) + DEPST(3)/2.0D0
 
      EPSIF(1,1) = EPSIF(1,1) + DEPST(4)
      EPSIF(2,2) = EPSIF(2,2) + DEPST(5)
      EPSIF(1,2) = EPSIF(1,2) + DEPST(6)/2.0D0
      EPSIF(2,1) = EPSIF(2,1) + DEPST(6)/2.0D0
C
C----PREDICTION ELASTIQUE-----------------------------------------------
C
 
C     Deformation hors plan
      CALL CALCEP33(EPSIM,XGAMMT,XGAMMC,XD1,XD2,XMUM,XLAMM,XEPS33)
C
C----CALCUL DES SEUILS YM1 ET YM2---------------------------------------
C
      CALL CALCYMI(EPSIM,XGAMMT,XGAMMC,XLAMM,XMUM,XD1,XEPS33,YM1)
      CALL CALCYMI(EPSIM,XGAMMT,XGAMMC,XLAMM,XMUM,XD2,XEPS33,YM2)
C
C----CALCUL DES SEUILS YF1 ET YF2---------------------------------------
C
 
C     Calcul de YF01 et YF02
      CALL CALCYF01(EPSIF,XGAMF,XLAMF,XMUF,YF01)
      CALL CALCYF02(EPSIF,XGAMF,XLAMF,XMUF,YF02)
 
C     Calcul de YF1 et YF2
      YF1 = (XALFA/(XD1+XALFA)**2.0D0)*YF01
      YF2 = (XALFA/(XD2+XALFA)**2.0D0)*YF02
C
C----CALCUL DES SEUILS Y1 ET Y2-----------------------------------------
C
      XSEUIL1 = YM1 + YF1 - XSEUI0
      XSEUIL2 = YM2 + YF2 - XSEUI0
C
C----TEST SEUIL 1-------------------------------------------------------
C
      IF (XSEUIL1.GT.0.0D0) THEN
 
C     Debut des iterations internes
      DO NINTER=1,100
 
C     Calcul du residu 1
      CALL CALCYMI(EPSIM,XGAMMT,XGAMMC,XLAMM,XMUM,XD1,XEPS33,YM1)
 
      CALL CALCYF01(EPSIF,XGAMF,XLAMF,XMUF,YF01)
 
      YF1 = (XALFA/(XD1+XALFA)**2.0D0)*YF01
 
      XRES1 = YM1 + YF1 - XSEUI0
 
C     Calcul de l'operateur tangent 1
      CALL CALCRTG1(EPSIM,EPSIF,XGAMMT,XGAMMC,XGAMF,XALFA,XLAMM,
     &              XMUM,XLAMF,XMUF,XD1,XD2,XEPS33,DR1DD1)
 
C     Corrections de l'endommagement
      XDD1 = -1.0D0*(1.0D0/DR1DD1)*XRES1
 
C     Critere d arret des iterations internes
      XCRI1 = ABS(XDD1*DR1DD1/XSEUI0)
 
      IF (XCRI1.LE.XTOL) THEN
 
         GOTO 111
 
      ELSE
 
         XD1 = XD1 + XDD1
 
         CALL CALCEP33(EPSIM,XGAMMT,XGAMMC,XD1,XD2,XMUM,XLAMM,
     &                 XEPS33)
 
      ENDIF
 
      ENDDO
      ENDIF
 
 111  CONTINUE
 
C
C----TEST SEUIL 2-------------------------------------------------------
C
      IF (XSEUIL2.GT.0.0D0) THEN
 
C     Debut des iterations internes
      DO NINTER=1,100
 
C     Calcul du residu 1
      CALL CALCYMI(EPSIM,XGAMMT,XGAMMC,XLAMM,XMUM,XD2,XEPS33,YM2)
 
      CALL CALCYF02(EPSIF,XGAMF,XLAMF,XMUF,YF02)
 
      YF2 = (XALFA/(XD2+XALFA)**2.0D0)*YF02
 
      XRES2 = YM2 + YF2 - XSEUI0
 
C     Calcul de l'operateur tangent 1
      CALL CALCRTG2(EPSIM,EPSIF,XGAMMT,XGAMMC,XGAMF,XALFA,XLAMM,
     &              XMUM,XLAMF,XMUF,XD1,XD2,XEPS33,DR2DD2)
 
C     Corrections de l'endommagement
      XDD2 = -1.0D0*(1.0D0/DR2DD2)*XRES2
 
C     Critere d arret des iterations internes
      XCRI2 = ABS(XDD2*DR2DD2/XSEUI0)
 
      IF (XCRI2.LE.XTOL) THEN
 
         GOTO 112
 
      ELSE
 
         XD2 = XD2 + XDD2
 
         CALL CALCEP33(EPSIM,XGAMMT,XGAMMC,XD1,XD2,XMUM,XLAMM,
     &                 XEPS33)
 
      ENDIF
 
      ENDDO
      ENDIF
 
 112  CONTINUE
 
      XTRA = EPSIM(1,1)+EPSIM(2,2)
      CALL CALCXIM (XTRA,XGAMMT,XGAMMC,XD1,XD2,XIMF)
C
C-----CALCUL DE CONTRAINTES---------------------------------------------
C
 
C     Test sur les endommagements pour obtenir une dissipation positive
      IF (XD1.LE.XD10) THEN
         XD1 = XD10
      ENDIF
 
      IF (XD2.LE.XD20) THEN
         XD2 = XD20
      ENDIF
 
C     Cacul du repere propre
      CALL JACOB4(EPSIM,2,VAPM,VEPM)
C         print*,'EPSIM=',((EPSIM(I,J),I=1,2),J=1,2)
C         print*,'VAPM=',(VAPM(I),I=1,3)
C         print*,'VEPM=',((VEPM(I,J),I=1,2),J=1,2)
 
          CALL JACOB4(EPSIF,2,VAPF,VEPF)
 
C         print*,'EPSIF=',((EPSIF(I,J),I=1,3),J=1,3)
C         print*,'VAPF=',(VAPF(I),I=1,3)
C         print*,'VEPF=',((VEPF(I,J),I=1,2),J=1,2)
 
C     Calcul de l'endommagement xim, xim1 et xim2
      XTRA = EPSIM(1,1)+EPSIM(2,2)
      CALL CALCXIM(XTRA,XGAMMT,XGAMMC,XD1,XD2,XIM)
      CALL CALCXIM(VAPM(1),XGAMMT,XGAMMC,XD1,XD2,XIM1)
      CALL CALCXIM(VAPM(2),XGAMMT,XGAMMC,XD1,XD2,XIM2)
 
      XIMI(1) = XIM1
      XIMI(2) = XIM2
 
C     Calcul de l'endommagement xif, xif1 et xif2
      XTRAK = EPSIF(1,1)+EPSIF(2,2)
      CALL CALCXIF(XTRAK,XALFA,XGAMF,XD1,XD2,XIF)
      CALL CALCXIF(VAPF(1),XALFA,XGAMF,XD1,XD2,XIF1)
      CALL CALCXIF(VAPF(2),XALFA,XGAMF,XD1,XD2,XIF2)
 
      XIFI(1) = XIF1
      XIFI(2) = XIF2
 
C     Calcul du repere propre des deformations generalisees
      CALL TRSPOD(VEPM,3,3,VEPMT)
      CALL TRSPOD(VEPF,2,2,VEPFT)
 
C     Calcul des contraintes generalisees dans le repere propre
      XTRA = EPSIM(1,1)+EPSIM(2,2)
      XTRAK = EPSIF(1,1)+EPSIF(2,2)
 
      CALL ZERO(SIGMAM,3,3)
      CALL ZERO(SIGMAF,3,3)
 
      XTRAK = EPSIF(1,1) + EPSIF(2,2)
 
      DO I = 1,2
         SIGMAM(I,I) = XIM*XLAMM*(XTRA+XEPS33)+
     &                 2.0D0*XMUM*VAPM(I)*XIMI(I)
 
C Abandon de la non linearite en partie flexion
C         SIGMAF(I,I) = XIF*XLAMF*XTRAK+
C     &                 2.0D0*XMUF*VAPF(I)*XIFI(I)
      ENDDO
 
C     Calcul des contraintes generalisees dans le repere global
      CALL ZERO(XTRAV5,3,3)
      CALL ZERO(XTRAV6,3,3)
 
      CALL MATMAT(VEPM,SIGMAM,3,3,3,XTRAV5)
      CALL ZERO(SIGMAM,3,3)
      CALL MATMAT(XTRAV5,VEPMT,3,3,3,SIGMAM)
 
C Abandon de la non linearite en partie flexion
C      CALL MATMAT(VEPF,SIGMAF,3,3,3,XTRAV6)
C      CALL ZERO(SIGMAF,3,3)
C      CALL MATMAT(XTRAV6,VEPFT,3,3,3,SIGMAF)
C
C----CONTRAINTES-------------------------------------------------------
C
      SIGF(1) = SIGMAM(1,1)
      SIGF(2) = SIGMAM(2,2)
      SIGF(3) = SIGMAM(1,2)
 
 
      FAC  = (EPAI**3)/12.D0
      AUX  = FAC*XYOUNF/(1.D0-XNUF*XNUF)
      AUX1 = FAC*XYOUNF*.5D0/(1.D0+XNUF)
      AUX2 = EPAI*XYOUNF*.5D0/(1.d0+XNUF)/1.2d0
 
      SIGF(4) = SIG0(4)+AUX*(DEPST(4)+XNUF*DEPST(5))
      SIGF(5) = SIG0(5)+AUX*(DEPST(5)+XNUF*DEPST(4))
      SIGF(6) = SIG0(6)+AUX1*DEPST(6)
 
C         print*,'SIGF(I)',(SIGF(I),I=4,6)
 
C Cisaillement élastique
      SIGF(7) = SIG0(7)+AUX2*DEPST(7)
      SIGF(8) = SIG0(8)+AUX2*DEPST(8)
 
C
C----VARIABLES INTERNES-------------------------------------------------
C
      VARF(1) = EPSIM(1,1)
      VARF(2) = EPSIM(2,2)
      VARF(3) = EPSIM(1,2)*2.0D0
 
      VARF(4) = EPSIF(1,1)
      VARF(5) = EPSIF(2,2)
      VARF(6) = EPSIF(1,2)*2.0D0
 
      VARF(9)  = XD1
      VARF(10) = XD2
 
      VARF(11) = 1.0D0-XIM
      VARF(12) = 1.0D0-XIF
 
      RETURN
      END