C RICBET3D SOURCE CB215821 16/04/21 21:18:15 8920 SUBROUTINE RICBET3D(XMAT,SIG0,SIGF,VAR0,VARF,DEPST,ISTEP) C C====&===1=========2=========3=========4=========5=========6=========7== C Commentaires : Subroutine permettant de mettre en oeuvre le C modele RICBET (B. Richard) pour representer C le comportement 3D d'un béton ordinaire sous C sollicitations cycliques C C Traits : - Endommagement scalaire C - Glissement frottant C - 3D C - Monotone et cyclique C - Effet unilatteral C - Comportement lineaire en compression C - Non localite portant sur le seuil d'endommagement (Fd) C C Auteurs : B. Richard (Dr - Ing.) & F. Ragueneau (Pr) - LMT/ENS C Date : Mars 2011 C====&===1=========2=========3=========4=========5=========6=========7== C C-----DECLARATION GENERALE---------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) C C-----DECLARATION GENERALE---------------------------------------------- C REAL*8 X(3,3),UNIT(3,3),EPSPI(3,3),EPSI(3,3),SIGMA(3,3) REAL*8 EPSDEV(3,3),SIGMA1(3,3),EPSPIDEV(3,3) REAL*8 EPSIPOS(3,3),EPSINEG(3,3),SIGPI(3,3),DEPSN(3,3) REAL*8 TERME1(3,3),TERME2(3,3),TERME3(3,3),TERME4(3,3) REAL*8 SIGMOX(3,3),TRAV(3,3),SIG1DEV(3,3),TERME5(3,3) REAL*8 XDT,XTRAV2(3,3),XTRAV3(3,3),DEFN(3,3),EPSIT(3,3) REAL*8 XMAT(*),SIG0(*),SIGF(*),VAR0(*),VARF(*),DEPST(*) REAL*8 XDFP(3,3),XDGP(3,3),XDJ2P(3,3),DEVGP(3,3) C C-----DEFINITION DE PARAMETRES------------------------------------------ C PARAMETER (IMAX = 1000) C C-----MISE EN DONNEE---------------------------------------------------- C XE = XMAT(1) XNU = XMAT(2) XFT = XMAT(5) ALDIR = XMAT(6) GAMMA1 = XMAT(7) A1 = XMAT(8) SIGRF = XMAT(9) XAF = XMAT(10) XBF = XMAT(11) XAG = XMAT(12) XBG = XMAT(13) XAC = XMAT(14) XBC = XMAT(15) XSIGU = XMAT(16) XFC = XMAT(17) XK = XE/((1.0D0-2.D0*XNU)*3.0D0) XG = XE/(2.0D0*(1.D0+XNU)) C C-----DEFINITION DE LA MATRICE UNITE D ORDRE 2------------------------------ C DO I=1,3 DO J=1,3 IF (I.EQ.J) THEN UNIT(I,J) = 1.0D0 ELSE UNIT(I,J) = 0.0D0 ENDIF ENDDO ENDDO C C-----VARIABLES INTERNES---------------------------------------------------- C C Tenseur de glissement interne EPSPI(1,1) = VAR0(2) EPSPI(2,2) = VAR0(3) EPSPI(3,3) = VAR0(4) EPSPI(1,2) = VAR0(5) EPSPI(1,3) = VAR0(6) EPSPI(2,3) = VAR0(7) EPSPI(2,1) = VAR0(5) EPSPI(3,1) = VAR0(6) EPSPI(3,2) = VAR0(7) C Tenseur de la contrainte de rappel X(1,1) = VAR0(8) X(2,2) = VAR0(9) X(3,3) = VAR0(10) X(1,2) = VAR0(11) X(1,3) = VAR0(12) X(2,3) = VAR0(13) X(2,1) = VAR0(11) X(3,1) = VAR0(12) X(3,2) = VAR0(13) C Endommagement XDT = VAR0(14) C Force associee a l ecrouissage isotrope XFFT = VAR0(15) C Deformations totales EPSI(1,1) = VAR0(16) EPSI(2,2) = VAR0(17) EPSI(3,3) = VAR0(18) EPSI(1,2) = VAR0(19) EPSI(1,3) = VAR0(20) EPSI(2,3) = VAR0(21) EPSI(2,1) = VAR0(19) EPSI(3,1) = VAR0(20) EPSI(3,2) = VAR0(21) C Tenseur des deformations plastiques DEFN(1,1) = VAR0(35) DEFN(2,2) = VAR0(36) DEFN(3,3) = VAR0(37) DEFN(1,2) = VAR0(38) DEFN(1,3) = VAR0(39) DEFN(2,3) = VAR0(40) DEFN(2,1) = VAR0(38) DEFN(3,1) = VAR0(39) DEFN(3,2) = VAR0(40) C Ecrouissage plastique XECO = VAR0(41) C Deformation plastique cumulee XDPP = VAR0(42) C Contraintes tests SIGMA1(1,1) = VAR0(22) SIGMA1(2,2) = VAR0(23) SIGMA1(3,3) = VAR0(24) SIGMA1(1,2) = VAR0(25) SIGMA1(1,3) = VAR0(26) SIGMA1(2,3) = VAR0(27) SIGMA1(2,1) = VAR0(25) SIGMA1(3,1) = VAR0(26) SIGMA1(3,2) = VAR0(27) C Indicateur XBEN = VAR0(22)+VAR0(23)+VAR0(24) C Contrainte initiales SIGMA(1,1) = SIG0(1) SIGMA(2,2) = SIG0(2) SIGMA(3,3) = SIG0(3) SIGMA(1,2) = SIG0(4) SIGMA(1,3) = SIG0(5) SIGMA(2,3) = SIG0(6) SIGMA(2,1) = SIG0(4) SIGMA(3,1) = SIG0(5) SIGMA(3,2) = SIG0(6) C C-----ACTUALISATION DES DEFORMATIONS------------------------------------ C EPSI(1,1) = EPSI(1,1) + DEPST(1) EPSI(2,2) = EPSI(2,2) + DEPST(2) EPSI(3,3) = EPSI(3,3) + DEPST(3) EPSI(1,2) = EPSI(1,2) + (DEPST(4)/2.0D0) EPSI(1,3) = EPSI(1,3) + (DEPST(5)/2.0D0) EPSI(2,3) = EPSI(2,3) + (DEPST(6)/2.0D0) EPSI(2,1) = EPSI(2,1) + (DEPST(4)/2.0D0) EPSI(3,1) = EPSI(3,1) + (DEPST(5)/2.0D0) EPSI(3,2) = EPSI(3,2) + (DEPST(6)/2.0D0) C C-----CALCUL DE QUELQUES QUANTITES-------------------------------------- C C Calcul de la trace de EPSI TRA = EPSI(1,1)+EPSI(2,2)+EPSI(3,3) C Calcul du deviateur de EPSI CALL DEVIAT(EPSI,EPSDEV) C Calcul de la trace de EPSPI TRAPI = EPSPI(1,1)+EPSPI(2,2)+EPSPI(3,3) C Calcul du deviateur de EPSPI CALL DEVIAT(EPSPI,EPSPIDEV) C Calcul du deviateur de DEFN CALL DEVIAT(DEFN,DEPSN) C Calcul de la trace de DEFN TRADN = DEFN(1,1)+DEFN(2,2)+DEFN(3,3) C C-----PARTIE POSITIVES ET NEGATIVES DE (EPS)----------------------------- C DO I=1,3 DO J=1,3 EPSIT(I,J) = EPSI(I,J)-DEFN(I,J) ENDDO ENDDO CALL POSNEG(EPSIT,EPSIPOS,EPSINEG) C C-----SEUIL ENDOMMAGEMENT------------------------------------------------ C C Calcul de l enerigie Y+ YP = 0.0D0 DO I=1,3 DO J=1,3 YP = YP + EPSIPOS(I,J)*EPSIPOS(I,J) ENDDO ENDDO YP = 0.5D0*YP*XE C C-----NON LOCAL---------------------------------------------------------- C IF (ISTEP.EQ.0) THEN SEUIL1=YP VARF(1)=YP ELSE IF (ISTEP.EQ.1) THEN VARF(1)=YP DO I=2,42 VARF(I)=VAR0(I) ENDDO GOTO 2000 ELSE IF (ISTEP.EQ.2) THEN SEUIL1=VAR0(1) VARF(1)=SEUIL1 ENDIF C C-----EVOLUTION DE L ENDOMMAGEMENT--------------------------------------- C C Calcul du seuil XFD = SEUIL1 - ((XFT**2)/(2.0D0*XE) + XFFT) C Evolution explicite de l'endommagement IF ((XFD.GT.0.0D0).AND.(XBEN.GE.0.0D0)) THEN XDT = 1.0D0-1.0D0/(1.0D0+ALDIR*(SEUIL1-((XFT**2)/(2.0D0*XE)))) XFFT = SEUIL1 - ((XFT**2)/(2.0D0*XE)) ENDIF C C-----SEUIL GLISSEMENT INTERNE------------------------------------------- C C Calcul de sigma_pi DO I=1,3 DO J=1,3 SIGPI(I,J) = XDT*XK*(TRA-TRAPI-TRADN)*UNIT(I,J)+ & 2.0D0*XG*XDT*(EPSDEV(I,J)-EPSPIDEV(I,J)-DEPSN(I,J)) ENDDO ENDDO C Calcul de la quantite (sigma_pi - X) DO I=1,3 DO J=1,3 SIGMOX(I,J) = SIGPI(I,J)-X(I,J) ENDDO ENDDO C Calcul de I2(sigma_pi-X) CALL I2(SIGMOX,XFPI) C Sigma_test = Sigma_finale DO I=1,3 DO J=1,3 SIGMA(I,J)=SIGMA1(I,J) ENDDO ENDDO C Mise a 1 du multiplicateur plastique XLAM = 1.0D0 C CAlcul de la trace de Sigma_finale TRASIG1 = SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3) C Test sur la positivite du seuil de glissement interne IF ((XFPI.GT.0.0D0).AND.(XBEN.GE.0.D0)) THEN C On retient la valeur initiale du seuil XFPI0 = XFPI C On itere en bloquant le nombre d iteration max a IMAX DO NI=1,IMAX C On realise un test de sortie IF ((ABS(XFPI/XFPI0 ).LE.1.0D-8).OR. & (ABS(XLAM).LE.(1.0D-10))) THEN GOTO 400 ELSE C Calcul des derivees dfpi/dsigpi et dphipi/dsigpi CALL DI2(SIGMOX,TERME1) CALL DI2(SIGMOX,TERME2) C Calcul des derivees dfpi/dX et dphipi/dX DO I=1,3 DO J=1,3 TERME3(I,J) = 0.0D0 TERME4(I,J) = 0.0D0 TERME3(I,J) = -1.0D0*TERME2(I,J) TERME4(I,J) = -1.0D0*TERME2(I,J)+A1*X(I,J) ENDDO ENDDO C On calcul le terme C:dfpi/dsigpi TRA_SE = TERME2(1,1)+TERME2(2,2)+TERME2(3,3) CALL DEVIAT(TERME2,TRAV) DO I=1,3 DO J=1,3 TERME5(I,J) = 0.0D0 TERME5(I,J) = XK*TRA_SE*UNIT(I,J)+ & 2.0D0*XG*TRAV(I,J) ENDDO ENDDO C On calcul le denominateur du multiplicateur plastique XT1 = 0.0D0 DO I=1,3 DO J=1,3 XT1 = XT1 + TERME1(I,J)*TERME5(I,J) ENDDO ENDDO XT2 = 0.0D0 DO I=1,3 DO J=1,3 XT2 = XT2 + TERME3(I,J)*TERME4(I,J) ENDDO ENDDO C Calcul du multiplicateur plastique XLAM = XFPI/(XDT*XT1 + GAMMA1*XT2) C On actualisatise X et sigpi DO I=1,3 DO J=1,3 SIGPI(I,J)=SIGPI(I,J)-XDT*XLAM*TERME5(I,J) X(I,J)=X(I,J)-GAMMA1*XLAM*TERME4(I,J) ENDDO ENDDO C On recalcule le nouveau seuil DO I=1,3 DO J=1,3 SIGMOX(I,J) = 0.0D0 SIGMOX(I,J) = SIGPI(I,J)-X(I,J) ENDDO ENDDO CALL I2(SIGMOX,XFPI) ENDIF ENDDO 400 CONTINUE C On calcul epspi TRA_SE=SIGPI(1,1)+SIGPI(2,2)+SIGPI(3,3) CALL DEVIAT(SIGPI,TRAV) TRAPI = (TRA-TRADN) - (TRA_SE/(XDT*3.0D0*XK)) DO I=1,3 DO J=1,3 EPSPIDEV(I,J) = (EPSDEV(I,J)-DEPSN(I,J))- & (1.0D0/(2.0D0*XG*XDT))*TRAV(I,J) ENDDO ENDDO DO I=1,3 DO J=1,3 EPSPI(I,J) = (TRAPI/3.0D0)*UNIT(I,J)+EPSPIDEV(I,J) ENDDO ENDDO ENDIF C C-----CONTRAINTES ACTUALISEE C TRA1 = SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3) IF (XBEN.GE.0.0D0) THEN DO I=1,3 DO J=1,3 SIGMA1(I,J)=XK*(TRA-XDT*TRAPI-TRADN)*UNIT(I,J)+ & 2.0D0*XG* & (EPSDEV(I,J)-XDT*EPSPIDEV(I,J)-DEPSN(I,J)) ENDDO ENDDO ELSEIF (XBEN.LE.SIGRF) THEN C La contrainte effective DO I=1,3 DO J=1,3 SIGMA1(I,J)=XK*(TRA-TRADN)*UNIT(I,J)+ & 2.0D0*XG*(EPSDEV(I,J)-DEPSN(I,J)) ENDDO ENDDO C On verifie le seuile de plasticite CALL J2(SIGMA1,TEMP1) TEMP2=SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3) TEMP2 = TEMP2/3.00 XFP = XAF*TEMP1+XBF*TEMP2 - (XFC + XECO) C On verifie la positivite du seuil IF (XFP.GT.0.0D0) THEN C On retient la valeur initiale pour le critere relatif XFP0 = XFP C On commence les choses serieuses DO NP=1,IMAX C La derivee du critere CALL DI2(SIGMA1,XDJ2P) DO I=1,3 DO J=1,3 XDFP(I,J) = XAF*XDJ2P(I,J)+XBF*UNIT(I,J)/3.0D0 XDGP(I,J) = XAG*XDJ2P(I,J)+XBG*UNIT(I,J)/3.0D0 ENDDO ENDDO C Le hessien elasique CALL DEVIAT(XDGP,DEVGP) XTGP = XDGP(1,1)+XDGP(2,2)+XDGP(3,3) DO I=1,3 DO J=1,3 TRAV(I,J)=0.0D0 TRAV(I,J)=XK*XTGP*UNIT(I,J)+2.0D0*XG*DEVGP(I,J) ENDDO ENDDO C La derivee de l ecrouissage XTEMP= EXP(-1.0D0*XBC*XDPP)*(XAC-(XAC*XDPP+XFC+XSIGU)*XBC) C Le critere actuel CALL J2(SIGMA1,TEMP1) TEMP2=SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3) TEMP2 = TEMP2/3.00 XFP = XAF*TEMP1+XBF*TEMP2 - (XFC + XECO) C Le multiplicateur plastique CALL DBLECONT(XDFP,TRAV,TEMP9) XLAMP = XFP/(TEMP9+XTEMP) C Actualisation des variables internes XECO = XECO + XLAMP*XTEMP XDPP = XDPP + XLAMP DO I=1,3 DO J=1,3 DEFN(I,J)=DEFN(I,J) + XLAMP*XDGP(I,J) ENDDO ENDDO DO I=1,3 DO J=1,3 SIGMA1(I,J)=SIGMA1(I,J) - XLAMP*TRAV(I,J) ENDDO ENDDO C On calcul le seuil pour un eventuel arret CALL J2(SIGMA1,TEMP1) TEMP2=SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3) TEMP2 = TEMP2/3.00 XFP = XAF*TEMP1+XBF*TEMP2 - (XFC + XECO) C Le critere d arret des iterations internes IF ((ABS(XFP/XFP0).LE.(1.0D-8)).OR.(NP.GE.IMAX)) THEN GOTO 294 ENDIF ENDDO 294 CONTINUE ENDIF ELSE C Calcul du deviateur de DEFN CALL DEVIAT(DEFN,DEPSN) C Calcul de la trace de DEFN TRADN = DEFN(1,1)+DEFN(2,2)+DEFN(3,3) TRASIG1 = 3.0D0*XK*(TRA-XDT*TRAPI-TRADN)/ & (1.0D0-(3.0D0*XK/SIGRF)*XDT*TRAPI) DO I=1,3 DO J=1,3 SIG1DEV(I,J)=2.00D0*XG*(EPSDEV(I,J)-XDT*EPSPIDEV(I,J)* & (1.0D0-(TRASIG1/SIGRF))-DEPSN(I,J)) ENDDO ENDDO DO I=1,3 DO J=1,3 SIGMA1(I,J)=TRASIG1*UNIT(I,J)/3.0D0+SIG1DEV(I,J) ENDDO ENDDO ENDIF C C-----CALCUL DE L ENDOMMAGEMENT EFFECTIF (GRANDEUR INDICATIVE) C-----PAS DE SENS THERMO C XFSIG = 1.0D0-(SIGMA1(1,1)+SIGMA1(2,2)+SIGMA1(3,3))/SIGRF IF ((XFSIG.LE.1.0D0).AND.(XFSIG.GE.0.0D0)) THEN XDEFF = XDT*XFSIG ELSE XDEFF = XDT ENDIF C C-----STOCKAGE DES CONTRAINTES AVEC SYMMETRISATION POUR ROBUSTESSE C SIGF(1)= SIGMA1(1,1) SIGF(2)= SIGMA1(2,2) SIGF(3)= SIGMA1(3,3) SIGF(4)= (SIGMA1(1,2)+SIGMA1(2,1))/2.0D0 SIGF(5)= (SIGMA1(1,3)+SIGMA1(3,1))/2.0D0 SIGF(6)= (SIGMA1(2,3)+SIGMA1(3,2))/2.0D0 C C-----STOCKAGE DES DEFORMATIONS DE GLISSEMENT C VARF(2)=EPSPI(1,1) VARF(3)=EPSPI(2,2) VARF(4)=EPSPI(3,3) VARF(5)=EPSPI(1,2) VARF(6)=EPSPI(1,3) VARF(7)=EPSPI(2,3) C C-----STOCKAGE DE L ECROUISSAGE CINEMATIQUE C VARF(8)=X(1,1) VARF(9)=X(2,2) VARF(10)=X(3,3) VARF(11)=X(1,2) VARF(12)=X(1,3) VARF(13)=X(2,3) C C-----STOCKAGE DES CONTRAINTES TESTS C VARF(22)=SIGMA1(1,1) VARF(23)=SIGMA1(2,2) VARF(24)=SIGMA1(3,3) VARF(25)=SIGMA1(1,2) VARF(26)=SIGMA1(1,3) VARF(27)=SIGMA1(2,3) C C-----TENSEUR DES DEFORMATIONS PLASTIQUES C VARF(35) = DEFN(1,1) VARF(36) = DEFN(2,2) VARF(37) = DEFN(3,3) VARF(38) = DEFN(1,2) VARF(39) = DEFN(1,3) VARF(40) = DEFN(2,3) C C-----Ecrouissage plastique C VARF(41) = XECO C C-----Deformation plastique cumulee C VARF(42) = XDPP C C-----STOCKAGE DES CONTRAINTES TESTS C VARF(29)=SIGPI(1,1) VARF(30)=SIGPI(2,2) VARF(31)=SIGPI(3,3) VARF(32)=SIGPI(1,2) VARF(33)=SIGPI(1,3) VARF(34)=SIGPI(2,3) C C-----ENDOMMAGEMENT EFFECTIF C VARF(28) = XDEFF C C-----STOCKAGE DE L ENDOMMAGEMENT C VARF(14) = XDT C C-----STOCKAGE DES FORCES THERMODYNAMIQUES ASSOCIEES A L ECROUISSAGE ISOTROPE C VARF(15) = XFFT C C-----STOCKAGE DES DEFORMATIONS TOTALES C VARF(16)=EPSI(1,1) VARF(17)=EPSI(2,2) VARF(18)=EPSI(3,3) VARF(19)=EPSI(1,2) VARF(20)=EPSI(1,3) VARF(21)=EPSI(2,3) C C-----BALISE DE SORTIE POUR LE NON LOCAL C 2000 CONTINUE RETURN END