C CONCYC2 SOURCE CB215821 16/04/21 21:15:55 8920 SUBROUTINE CONCYC2(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 CONCYC (M. Vassaux & B. Richard) pour representer C le comportement 2D/3D d'un beton ordinaire sous C sollicitations cycliques C C Traits : - Endommagement scalaire C - Boucles hysteretiques C - Deformations permanentes C - Refermeture lineaire des fissures C - Effet unilatteral complet C - Non localite portant sur le seuil d'endommagement (Fd) C C Auteur : M. Vassaux C C Co-auteur : B. Richard C C Date : 2010 - 2011 C====&===1=========2=========3=========4=========5=========6=========7== C C----DECLARATION GENERALES---------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) C C-----DECLARATION GENERALE---------------------------------------------- C REAL*8 UNIT(3,3) REAL*8 SIGTOT(3,3),SIGMAT(3,3),SIGFIS(3,3),QF(3,3),SIGFISEFF(3,3) REAL*8 EPSTOT(3,3),EPSTOTPREC(3,3) REAL*8 EPSFIS(3,3),DEPSFIS(3,3),EPSFISMAX(3,3),EPSFISP(3,3) REAL*8 J11(3,3),J22(3,3),J12(3,3),J21(3,3),DSFIEFIJ12(3,3) REAL*8 SIGFISEFFDEV(3,3),J12DEV(3,3) REAL*8 EPSTOTDEV(3,3),EPSFISDEV(3,3),DEPSFISDEV(3,3) REAL*8 EPSTOTPOS(3,3),EPSTOTNEG(3,3) REAL*8 DI2SIGFISEFFDEV(3,3),EPSIPP(3),VALP33(3,3) REAL*8 SIGELA(3,3), SIGELADEV(3,3),SIGIPP(3),VALS33(3,3) REAL*8 XMAT(*),SIG0(*),SIGF(*),VAR0(*),VARF(*),DEPST(*) C C-----DEFINITION DE LA MATRICE UNITE D ORDRE 3------------------------------ 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-----DEFINITION DE PARAMETRES------------------------------------------ C IMAX = 100 C C-----CHARGEMENT DES PARAMETRES MATERIAUX---------------------------------------------------- C XE = XMAT(1) XNU = XMAT(2) XK = XE/((1.0D0-2.D0*XNU)*3.0D0) XG = XE/(2.0D0*(1.D0+XNU)) XEND = XMAT(5) XSIGT = XMAT(6) ATRA = XMAT(7) BTRA = XMAT(8) XQP = XMAT(9) XCF = XMAT(10) C Ecrouissage cinematique XHF = 0.0 C Non-linearite de l'ecrouissage cinematique XAF = 0.0 C C-----CHARGEMENT DES VARIABLES INTERNES---------------------------------------------------- C C Deformations totales EPSTOT(1,1) = VAR0(1) EPSTOT(2,2) = VAR0(2) EPSTOT(3,3) = VAR0(3) EPSTOT(1,2) = VAR0(4) EPSTOT(1,3) = VAR0(5) EPSTOT(2,3) = VAR0(6) EPSTOT(2,1) = VAR0(4) EPSTOT(3,1) = VAR0(5) EPSTOT(3,2) = VAR0(6) C Contrainte initiales SIGTOT(1,1) = SIG0(1) SIGTOT(2,2) = SIG0(2) SIGTOT(3,3) = SIG0(3) SIGTOT(1,2) = SIG0(4) SIGTOT(1,3) = SIG0(5) SIGTOT(2,3) = SIG0(6) SIGTOT(2,1) = SIG0(4) SIGTOT(3,1) = SIG0(5) SIGTOT(3,2) = SIG0(6) C Contrainte fissures initiales SIGFIS(1,1) = VAR0(7) SIGFIS(2,2) = VAR0(8) SIGFIS(3,3) = VAR0(9) SIGFIS(1,2) = VAR0(10) SIGFIS(1,3) = VAR0(11) SIGFIS(2,3) = VAR0(12) SIGFIS(2,1) = VAR0(10) SIGFIS(3,1) = VAR0(11) SIGFIS(3,2) = VAR0(12) C Tenseur d'ecrouissage cinematique QF(1,1) = VAR0(13) QF(2,2) = VAR0(14) QF(3,3) = VAR0(15) QF(1,2) = VAR0(16) QF(1,3) = VAR0(17) QF(2,3) = VAR0(18) QF(2,1) = VAR0(16) QF(3,1) = VAR0(17) QF(3,2) = VAR0(18) C Endommagement XDT = VAR0(19) XDT0= VAR0(19) XZT = VAR0(20) C Fonctions de refermeture associees a la trace et au deviateur XPHI = VAR0(21) C Deformations associées a la fissuration maximale enregistrees EPSFISMAX(1,1) = VAR0(22) EPSFISMAX(2,2) = VAR0(23) EPSFISMAX(3,3) = VAR0(24) EPSFISMAX(1,2) = VAR0(25) EPSFISMAX(1,3) = VAR0(26) EPSFISMAX(2,3) = VAR0(27) EPSFISMAX(2,1) = VAR0(25) EPSFISMAX(3,1) = VAR0(26) EPSFISMAX(3,2) = VAR0(27) C Deformations plastiques associées au glissement EPSFISP(1,1) = VAR0(28) EPSFISP(2,2) = VAR0(29) EPSFISP(3,3) = VAR0(30) EPSFISP(1,2) = VAR0(31) EPSFISP(1,3) = VAR0(32) EPSFISP(2,3) = VAR0(33) EPSFISP(2,1) = VAR0(31) EPSFISP(3,1) = VAR0(32) EPSFISP(3,2) = VAR0(33) C Stockage des variables au pas precedent XDTPREC = XDT DO I=1,3 DO J=1,3 EPSTOTPREC(I,J) = EPSTOT(I,J) ENDDO ENDDO C Fonctions de refermeture associees a la trace et au deviateur XEDISS = VAR0(34) C C-----ACTUALISATION DES DEFORMATIONS------------------------------------ C EPSTOT(1,1) = EPSTOT(1,1) + DEPST(1) EPSTOT(2,2) = EPSTOT(2,2) + DEPST(2) EPSTOT(3,3) = EPSTOT(3,3) + DEPST(3) EPSTOT(1,2) = EPSTOT(1,2) + (DEPST(4)/2.0D0) EPSTOT(1,3) = EPSTOT(1,3) + (DEPST(5)/2.0D0) EPSTOT(2,3) = EPSTOT(2,3) + (DEPST(6)/2.0D0) EPSTOT(2,1) = EPSTOT(2,1) + (DEPST(4)/2.0D0) EPSTOT(3,1) = EPSTOT(3,1) + (DEPST(5)/2.0D0) EPSTOT(3,2) = EPSTOT(3,2) + (DEPST(6)/2.0D0) CALL JACOB3(EPSTOT,3,EPSIPP,VALP33) EPSTIL=MAX(0.0D0,EPSIPP(1))**2.0D0 + 1 MAX(0.0D0,EPSIPP(2))**2.0D0 + 2 MAX(0.0D0,EPSIPP(3))**2.0D0 EPSTIL=SQRT(EPSTIL) XEPSTOTTRAC = EPSTOT(1,1)+EPSTOT(2,2)+EPSTOT(3,3) CALL DEVIAT(EPSTOT,EPSTOTDEV) DO I=1,3 DO J=1,3 SIGELA(I,J) = (XK*XEPSTOTTRAC*UNIT(I,J) & + 2.0D0*XG*EPSTOTDEV(I,J)) ENDDO ENDDO CALL JACOB3(SIGELA,3,SIGIPP,VALP33) C-----Modele endommagement RICRAG classique ou avec consolidation modifiee IF (XEND.EQ.1.OR.XEND.EQ.2) THEN YP = 0.5D0*XSIGT*EPSTIL C-----Modele endommagement RICRAG avec critere modifie ELSEIF (XEND.EQ.3) THEN XSI1 = MIN(0.0D0,SIGIPP(1)) + MIN(0.0D0,SIGIPP(2)) & + MIN(0.0D0,SIGIPP(3)) EPSTIL=EPSTIL + BTRA*MIN(0.0D0,XSI1)/(3.0*XK) YP = 0.5D0*XSIGT*EPSTIL ENDIF C C-----NON LOCAL---------------------------------------------------------- C IF (ISTEP.EQ.0) THEN SEUIL1=YP VARF(35)=YP ELSE IF (ISTEP.EQ.1) THEN VARF(35)=YP DO I=1,34 VARF(I)=VAR0(I) ENDDO GOTO 2000 ELSE IF (ISTEP.EQ.2) THEN SEUIL1=VAR0(35) VARF(35)=SEUIL1 ENDIF C C-----EVOLUTION DE L ENDOMMAGEMENT--------------------------------------- C C Modele endommagement RICRAG classique ou avec critere modifie IF (XEND.EQ.1.OR.XEND.EQ.3) THEN XFD = SEUIL1 - ((XSIGT**2)/(2.0D0*XE) + XZT) XYY0 = (XSIGT**2)/(2.0D0*XE*SEUIL1) IF (XFD.GT.0.0D0) THEN XDT = 1.0D0 1 - XYY0*EXP(-ATRA*(SEUIL1-((XSIGT**2)/(2.0D0*XE)))) XZT = SEUIL1 - ((XSIGT**2)/(2.0D0*XE)) ENDIF C Modele endommagement RICRAG avec consolidation modifiee ELSEIF (XEND.EQ.2) THEN XSIGEELA = SIGIPP(1)**2.0D0 1 + SIGIPP(2)**2.0D0 2 + SIGIPP(3)**2.0D0 XSIGEELANEG = MIN(0.0D0,SIGIPP(1))**2.0D0 1 + MIN(0.0D0,SIGIPP(2))**2.0D0 2 + MIN(0.0D0,SIGIPP(3))**2.0D0 XCONF = (1 + BTRA*(XSIGEELANEG/XSIGEELA)**0.5) XYY0 = (XSIGT**2)/(2.0D0*XE*SEUIL1) XFD = SEUIL1 - ((XSIGT**2)/(2.0D0*XE) + XZT) IF (XFD.GT.0.0D0) THEN XDT = 1.0D0 & - XYY0*EXP(-(ATRA/XCONF)*(SEUIL1-((XSIGT**2)/(2.0D0*XE)))) XZT = SEUIL1 - ((XSIGT**2)/(2.0D0*XE)) ENDIF ENDIF C C-----ACTUALISATION DE LA CONTRAINTE DANS LA MATRICE FISSUREE------------ C C Calcul de la trace de EPSTOT XEPSTOTTRAC = EPSTOT(1,1)+EPSTOT(2,2)+EPSTOT(3,3) C Calcul du deviateur de EPSTOT CALL DEVIAT(EPSTOT,EPSTOTDEV) DO I=1,3 DO J=1,3 SIGMAT(I,J) = (1.0D0-XDT)*(XK*XEPSTOTTRAC*UNIT(I,J) & + 2.0D0*XG*EPSTOTDEV(I,J)) ENDDO ENDDO XI1SIGMAT = (SIGMAT(1,1)+SIGMAT(2,2)+SIGMAT(3,3))/3.0D0 C C-----ACTUALISATION DE LA DEFORMATION ASSOCIEE AUX FISSURES------------- C DO I=1,3 DO J=1,3 EPSFIS(I,J) = EPSTOT(I,J)*XDT ENDDO ENDDO DO I=1,3 DO J=1,3 DEPSFIS(I,J) = EPSTOT(I,J)*XDT & - EPSTOTPREC(I,J)*XDTPREC ENDDO ENDDO C Calcul de la trace de EPSFIS XEPSFISTRAC = EPSFIS(1,1) + EPSFIS(2,2) + EPSFIS(3,3) IF (ABS(XEPSFISTRAC).LE.3.0D-12) THEN XEPSFISTRAC = 3.0D-12 ENDIF C XEPSFISTRAC = EPSFIS(3,3) C Calcul du deviateur de DEPSFIS CALL DEVIAT(EPSFIS,EPSFISDEV) C Calcul de la trace de EPSFIS XEPSFISMAXTRAC = EPSFISMAX(1,1) + EPSFISMAX(2,2) + EPSFISMAX(3,3) C Securite... IF (XEPSFISMAXTRAC.LE.0.0D0) THEN XEPSFISMAXTRAC = 3.0D-12 ENDIF C Mise a jour de EPSFISMAX IF (XEPSFISTRAC.GT.XEPSFISMAXTRAC) THEN DO I=1,3 DO J=1,3 EPSFISMAX(I,J) = EPSFIS(I,J) ENDDO ENDDO ENDIF C C-----CALCUL DES FONCTIONS DE REFERMETURE DE FISSURE-------------------- C C Calcul de la trace de EPSFISMAX XEPSFISMAXTRAC = EPSFISMAX(1,1) + EPSFISMAX(2,2) + EPSFISMAX(3,3) C Securite... IF (XEPSFISMAXTRAC.LE.0.0D0) THEN XEPSFISMAXTRAC = 3.0D-12 ENDIF C Calcul de la trace de EPSFISP XEPSFISPTRAC = EPSFISP(1,1) + EPSFISP(2,2) + EPSFISP(3,3) C Calcul de la fonction de refermeture XPHI = 1.0D0 - 1.0D0/(1.0D0 + EXP(-1.0D0*XQP* & XEPSFISTRAC/XEPSFISMAXTRAC)) C C-----CALCUL DE LA CONTRAINTE DANS LES FISSURES: TEST-------------------- C C Calcul de la trace de EPSFIS XEPSFISTRAC = EPSFIS(1,1) + EPSFIS(2,2) + EPSFIS(3,3) C Calcul de la trace de DEPSFIS XDEPSFISTRAC = DEPSFIS(1,1)+DEPSFIS(2,2)+DEPSFIS(3,3) C Calcul du deviateur de DEPSFIS CALL DEVIAT(DEPSFIS,DEPSFISDEV) DO I=1,3 DO J=1,3 SIGFIS(I,J) = SIGFIS(I,J) & + XPHI*(XK*XDEPSFISTRAC)*UNIT(I,J) & + XPHI*(2.0D0*XG*DEPSFISDEV(I,J)) ENDDO ENDDO C C-----SEUIL GLISSEMENT INTERNE FISSURES-------------------------------- C C Calcul du I1 de SIGFIS XI1SIGFIS = (SIGFIS(1,1)+SIGFIS(2,2)+SIGFIS(3,3))/3.0D0 XI1SIGTOT = XI1SIGMAT + XI1SIGFIS C Calcul de la quantite (SIGFIS - QF) DO I=1,3 DO J=1,3 SIGFISEFF(I,J) = SIGFIS(I,J) - QF(I,J) ENDDO ENDDO C Calcul de I2(DEV(SIGFIS - QF)) CALL DEVIAT(SIGFISEFF,SIGFISEFFDEV) CALL I2(SIGFISEFFDEV,XI2SIGFISEFFDEV) C Calcul de la non-linearite de la dependance au confinement XFPI = XI2SIGFISEFFDEV + XCF*(XI1SIGFIS) XTEDISS = 0.0D0 C Test sur la positivite du seuil de glissement interne IF (XFPI.GT.0.0D0) THEN XFPI0 = XFPI XLAM = 1.0D0 C Debut du "Return-Mapping" 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 C Si le test n est pas concluant... ELSE C Calcul de dI2(DEV(SIGFIS - QF)) IF (XI2SIGFISEFFDEV.EQ.0.0D0) THEN DO I=1,3 DO J=1,3 DI2SIGFISEFFDEV(I,J) = 0.0D0 ENDDO ENDDO ELSE DO I=1,3 DO J=1,3 DI2SIGFISEFFDEV(I,J) = & 1.5D0*SIGFISEFFDEV(I,J)/XI2SIGFISEFFDEV ENDDO ENDDO ENDIF C Calcul des jacobiens DO I=1,3 DO J=1,3 J11(I,J) = DI2SIGFISEFFDEV(I,J) & + XCF*(UNIT(I,J)/3.0D0) J12(I,J) = DI2SIGFISEFFDEV(I,J) J21(I,J) = -1.0D0*DI2SIGFISEFFDEV(I,J) J22(I,J) = -1.0D0*DI2SIGFISEFFDEV(I,J) & + XAF*QF(I,J) ENDDO ENDDO C Calcul de la trace de DEPSFIS XJ12TRAC = J12(1,1)+J12(2,2)+J12(3,3) C Calcul du deviateur de DEPSFIS CALL DEVIAT(J12,J12DEV) C Calcul de (dsigfis/depsfis):J12 DO I=1,3 DO J=1,3 DSFIEFIJ12(I,J)= & +1.0D0*((XPHI)*XJ12TRAC*UNIT(I,J) & + 2.0D0*XG*(XPHI)*J12DEV(I,J)) ENDDO ENDDO C Calcul de produits contractes de matrices XJ11JXY = 0.0D0 XJ21J22 = 0.0D0 DO I=1,3 DO J=1,3 XJ11JXY = XJ11JXY + J11(I,J)*DSFIEFIJ12(I,J) XJ21J22 = XJ21J22 + J21(I,J)*J22(I,J) ENDDO ENDDO XDFPI = (XJ11JXY + XHF*XJ21J22) C Calcul du multiplicateur plastique IF(XDFPI.NE.0.0D0) THEN XLAM = XFPI/XDFPI ELSE GOTO 400 ENDIF C Actualisation de SIGFIS et QF DO I=1,3 DO J=1,3 SIGFIS(I,J)= SIGFIS(I,J) - XLAM*DSFIEFIJ12(I,J) QF(I,J) = QF(I,J) - XLAM*(XHF*J22(I,J)) EPSFISP(I,J) = EPSFISP(I,J) + XLAM*J12(I,J) ENDDO ENDDO C Calcul du I1 de SIGFIS XI1SIGFIS = & (SIGFIS(1,1)+SIGFIS(2,2)+SIGFIS(3,3))/3.0D0 XI1SIGTOT = XI1SIGMAT + XI1SIGFIS C Calcul de la quantite (SIGFIS - QF) DO I=1,3 DO J=1,3 SIGFISEFF(I,J) = SIGFIS(I,J) - QF(I,J) ENDDO ENDDO C Calcul de I2(DEV(SIGFIS - QF)) CALL DEVIAT(SIGFISEFF,SIGFISEFFDEV) CALL I2(SIGFISEFFDEV,XI2SIGFISEFFDEV) C Calcul du critere de glissement XFPI = XI2SIGFISEFFDEV + XCF*(XI1SIGFIS) C Calcul de l'increment d'energie dissipee XTEDISS = 0.0D0 DO I=1,3 DO J=1,3 XTEDISS = XTEDISS & + ABS(SIGFISEFFDEV(I,J)*XLAM*J12(I,J)) ENDDO ENDDO ENDIF ENDDO 400 CONTINUE ENDIF C C-----CALCUL DE LA CONTRAINTE TOTALE C DO I=1,3 DO J=1,3 SIGTOT(I,J) = SIGMAT(I,J) + SIGFIS(I,J) ENDDO ENDDO C PRINT *,'----------------FIN ITERATION------------------' C C-----CALCUL DE L ENERGIE DISSIPEE C XEDISS = XEDISS + XTEDISS C C-----STOCKAGE DES DEFORMATIONS TOTALES C VARF(1)=EPSTOT(1,1) VARF(2)=EPSTOT(2,2) VARF(3)=EPSTOT(3,3) VARF(4)=EPSTOT(1,2) VARF(5)=EPSTOT(1,3) VARF(6)=EPSTOT(2,3) C C-----STOCKAGE DES CONTRAINTES AVEC SYMMETRISATION POUR ROBUSTESSE C SIGF(1)= SIGTOT(1,1) SIGF(2)= SIGTOT(2,2) SIGF(3)= SIGTOT(3,3) SIGF(4)= (SIGTOT(1,2)+SIGTOT(2,1))/2.0D0 SIGF(5)= (SIGTOT(1,3)+SIGTOT(3,1))/2.0D0 SIGF(6)= (SIGTOT(2,3)+SIGTOT(3,2))/2.0D0 C C-----STOCKAGE DES CONTRAINTES DANS LES FISSURES C VARF(7)=SIGFIS(1,1) VARF(8)=SIGFIS(2,2) VARF(9)=SIGFIS(3,3) VARF(10)=(SIGFIS(1,2)+SIGFIS(2,1))/2.0D0 VARF(11)=(SIGFIS(1,3)+SIGFIS(3,1))/2.0D0 VARF(12)=(SIGFIS(2,3)+SIGFIS(3,2))/2.0D0 C C-----STOCKAGE DE L ECROUISSAGE CINEMATIQUE DE GLISSEMENT C VARF(13) =QF(1,1) VARF(14) =QF(2,2) VARF(15) =QF(3,3) VARF(16)=(QF(1,2)+QF(2,1))/2.0D0 VARF(17)=(QF(1,3)+QF(3,1))/2.0D0 VARF(18)=(QF(2,3)+QF(3,2))/2.0D0 C C-----STOCKAGE DE L ENDOMMAGEMENT C VARF(19) = XDT VARF(20) = XZT C C-----STOCKAGE DES FONCTIONS DE REFERMETURE C VARF(21) = XPHI C C-----STOCKAGE DEFORMATIONS MAXIMALES ASSOCIEES AU FISSURES C VARF(22) = EPSFISMAX(1,1) VARF(23) = EPSFISMAX(2,2) VARF(24) = EPSFISMAX(3,3) VARF(25) = (EPSFISMAX(1,2)+EPSFISMAX(2,1))/2.0D0 VARF(26) = (EPSFISMAX(1,3)+EPSFISMAX(3,1))/2.0D0 VARF(27) = (EPSFISMAX(2,3)+EPSFISMAX(3,2))/2.0D0 C-----STOCKAGE DEFORMATIONS PLASTIQUES ASSOCIEES AU GLISSEMENT C VARF(28) = EPSFISP(1,1) VARF(29) = EPSFISP(2,2) VARF(30) = EPSFISP(3,3) VARF(31) = (EPSFISP(1,2)+EPSFISP(2,1))/2.0D0 VARF(32) = (EPSFISP(1,3)+EPSFISP(3,1))/2.0D0 VARF(33) = (EPSFISP(2,3)+EPSFISP(3,2))/2.0D0 C C-----STOCKAGE DE L ENERGIE DISSIPEE PAR FRICTION C VARF(34) = XEDISS C C-----BALISE DE SORTIE POUR LE NON LOCAL C 2000 CONTINUE RETURN END