C LCGLDM SOURCE BR232186 12/10/04 21:15:18 7520 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 = DABS(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 = DABS(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