C RICJ3 SOURCE PV 21/11/04 21:15:16 11174 C sub ricjoi3d C RICJOI3D SOURCE CHAT 05/01/12 22:24:09 5004 SUBROUTINE RICJ3(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI, & DEPST,IFOUR,TETA1,TETA2,XMAT, & NMATT,IVAL,DD,SIGF,DEFP,VARF,KERRE) C----------------------------------------------------------------------- C Modele RICJOI en version element joint 3D C----------------------------------------------------------------------- C Développé par : Benjamin Richard C C Supervisé par : Frederic Ragueneau C Lucas Hector Adelaide C Christian Cremona C Jean Louis Tailhan C----------------------------------------------------------------------- C C-----DECLARATION GENERALE---------------------------------------------- C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) C C-----APPEL AUX INCLUDES------------------------------------------------ C -INC CCREEL C C-----PARAMETRES C PARAMETER(XDILAT = 80.0D0) C C-----DIMENSION--------------------------------------------------------- C DIMENSION SIG0(*),EPIN0(*),VAR0(*),DEPST(*),XMAT(*) DIMENSION DD(NSTRS,NSTRS) DIMENSION SIGF(*),DEFP(*),VARF(*) REAL*8 NCOEF,KCOEF FVARLIM = 2.0D-5 C C-----COMPTEUR POUR LE NOMBRE DE PASSAGE DAS LA BOUCLE C COUN = VAR0(17) C C-----ETAT PRECEDENT---------------------------------------------------- C D = VAR0(1 ) EPSF1 = VAR0(2 ) EPSF2 = VAR0(3 ) XEC1 = VAR0(4 ) XEC2 = VAR0(5 ) EPST1 = VAR0(6 ) EPST2 = VAR0(7 ) EPSN = VAR0(8 ) CIN1 = VAR0(9 ) CIN2 = VAR0(10) EPSR = VAR0(11) IFLAG = nint(VAR0(12)) FVAR = VAR0(13) SIGM = VAR0(14) XKS = VAR0(15) XKN = VAR0(16) crit2tol=1.001d0 C C-----INITIALISATION SPECIFIQUE POUR LES VARIABLES INTERNES DIFF. DE ZERO C IF (COUN.EQ.0.0D0) THEN FVAR = 0.2D0 SIGM = XMAT(13) XKS = XMAT(1) XKN = XMAT(2) ENDIF C C-----CONSTANTES MATERIAUX---------------------------------------------- C XKS0 = XMAT(1) XKN0 = XMAT(2) XAD = XMAT(5) XY0 = XMAT(6) XAL = XMAT(7) XGA = XMAT(8) XAA = XMAT(9) Q1COE = XMAT(10) Q2COE = XMAT(11) Q3COE = XMAT(12) SYCOE = XMAT(13) NCOEF = XMAT(14) KCOEF = XMAT(15) TCCOR = XMAT(16) NGONF = nint(XMAT(17)) CALL APP3(TCCOR,YR,15) C C-----ACTUALISATION DE LA DEFORMATION----------------------------------- C EPST1 = EPST1 + DEPST(1) EPST2 = EPST2 + DEPST(2) EPSN = EPSN + DEPST(3) IF (NGONF.EQ.1) THEN GOTO 123 ENDIF XSTE = XKN*EPSN XSAD = XDILAT*0.5D0*(-ABS(XSTE)+XSTE) C C-----CALCUL DE L ENERGIE----------------------------------------------- C EPSNP = 0.5D0*(ABS(EPSN)+EPSN) YN = 0.5D0*XKN*EPSNP*EPSNP YD1 = 0.5D0*XKS*EPST1*EPST1 YD2 = 0.5D0*XKS*EPST2*EPST2 C C-----CALCUL DU SEUIL D ENDOMMAGEMENT 1--------------------------------- C XFD1 = XAL*YN+YD1-(XY0+XEC1+YR) C C-----TEST SUR LE SEUIL D ENDOMMAGEMENT--------------------------------- C IF (XFD1.GT.0.0D0) THEN Dn = 1.0D0-1.0D0/(1+XAD*(XAL*YN+YD1-XY0-YR)) XEC1 = XAL*YN+YD1-XY0-YR IF (Dn.GT.D) THEN D = Dn ENDIF ENDIF IF (D.GE.1.0D0) THEN D = 0.999D0 ENDIF C C-----CALCUL DU SEUIL D ENDOMMAGEMENT 2--------------------------------- C XFD2 = XAL*YN+YD2-(XY0+XEC2+YR) C C-----TEST SUR LE SEUIL D ENDOMMAGEMENT--------------------------------- C IF (XFD2.GT.0.0D0) THEN Dn = 1.0D0-1.0D0/(1+XAD*(XAL*YN+YD2-XY0-YR)) XEC2 = XAL*YN+YD2-XY0-YR IF (Dn.GT.D) THEN D = Dn ENDIF ENDIF IF (D.GE.1.0D0) THEN D = 0.999D0 ENDIF C C-----FROTTEMENT 1------------------------------------------------------- C SPI1 = D*XKS*(EPST1-EPSF1) CRITF1= ABS(SPI1-CIN1)+XSAD IF (CRITF1.GT.(0.0D0)) THEN XF0 = CRITF1 XDLAM = 1.0D0 DO I=1,500 SEUIL = (ABS(SPI1-CIN1)+XSAD)/XF0 IF ((SEUIL.LE.1.0D-5).OR.(XDLAM.LE.1.0D-10)) THEN GOTO 900 ELSE XSG = ABS(SPI1-CIN1)/(SPI1-CIN1) XDLAM = (ABS(SPI1-CIN1)+XSAD)/ & (XKS*D-XSG*XGA*(-XSG+XAA*CIN1)) CIN1 = CIN1-XGA*XDLAM*(-XSG+XAA*CIN1) SPI1 = SPI1-D*XKS*XDLAM*XSG ENDIF EPSF1 = -SPI1/(D*XKS)+EPST1 ENDDO ENDIF 900 CONTINUE C C-----FROTTEMENT 2------------------------------------------------------- C SPI2 = D*XKS*(EPST2-EPSF2) CRITF2= ABS(SPI2-CIN2)+XSAD IF (CRITF2.GT.(0.0D0)) THEN XF0 = CRITF2 XDLAM = 1.0D0 DO I=1,500 SEUIL = (ABS(SPI2-CIN2)+XSAD)/XF0 IF ((SEUIL.LE.1.0D-5).OR.(XDLAM.LE.1.0D-10)) THEN GOTO 901 ELSE XSG = ABS(SPI2-CIN2)/(SPI2-CIN2) XDLAM = (ABS(SPI2-CIN2)+XSAD)/ & (XKS*D-XSG*XGA*(-XSG+XAA*CIN2)) CIN2 = CIN2-XGA*XDLAM*(-XSG+XAA*CIN2) SPI2 = SPI2-D*XKS*XDLAM*XSG ENDIF EPSF2 = -SPI2/(D*XKS)+EPST2 ENDDO ENDIF 901 CONTINUE C C-----ANELASTICITE-------------------------------------------------------- C IF (EPSN.GT.0.0D0) THEN SIGN = (1-D)*XKN*EPSN ELSE SIGN = XKN*(EPSN-EPSR) ENDIF C C-----CONTRAINTE DUE AUX EFFETS INELASTIQUE DE LA ROUILLE----------------- C sigr=sign IF (EPSN.LE.0.0D0) THEN SIGR = -1.0D0*SIGN ENDIF C C-----CRITERE D ANELASTICITE DE LA ROUILLE-------------------------------- C IF (IFLAG.EQ.1) THEN GOTO 20 ELSE CRIT2 = 2.0D0*Q1COE*FVAR*DCOSH(Q2COE*SIGR/2.0D0/SYCOE)- & (1.0D0+(Q3COE*FVAR)**2.0D0) ENDIF * print*,'crit2=',crit2 C C SI LE CRITERE EST DEPASSE ON DEMARRE LE RETURN MAPPING C IF ((CRIT2.GT.0.0D0).AND. 1 (EPSN .LE.0.0D0).AND. 2 (TCCOR.GT.0.0D0).AND. 3 (FVAR.GT.FVARLIM).AND. 4 (IFLAG.EQ.0).AND. 5 ((EPSN*DEPST(3)).GT.0.0D0)) THEN CRIT20 = CRIT2 DLAM2 = 1.0D0 GOTO 10 ELSE IF (FVAR.EQ.FVARLIM) THEN IFLAG = 1 GOTO 20 ELSE GOTO 30 END IF C C-----DEBUT DES ITERATIONS INTERNES C 10 DO I=1,10000 C C EVALUATION DU CRITERE C CRIT2 = 2.0D0*Q1COE*FVAR*DCOSH(Q2COE*SIGR/2.0D0/SYCOE)- & (1.0D0+(Q3COE*FVAR)**2.0D0) C C TEST DE CONVERGENCE SUR LA VALEUR DU CRITERE RELATIVE C IF ((ABS(CRIT2/CRIT20).LE.CRIT2TOL).OR. & (DLAM2.LT.0.0D0).OR. & (CRIT2.LT.0.0D0)) THEN GOTO 20 END IF C C CALCUL DES DERIVEES DU CRITERES C C dFr/dSIGR : DERIV21 C dFr/dSIGM : DERIV22 C dFr/dFVAR : DERIV23 C DERIV21 = SINH(Q2COE*SIGR/2.0D0/SIGM)* & Q1COE*Q2COE*FVAR/SIGM DERIV22 = SINH(Q2COE*SIGR/2.0D0/SIGM)* & (-1.0D0)*Q1COE*Q2COE*FVAR*SIGR/SIGM**2.0D0 DERIV23 = 2.0D0*Q1COE*DCOSH(Q2COE*SIGR/2.0D0/SIGM)- & 2.0D0*(Q3COE**2.0D0)*FVAR C C CALUL DU MODULE TANGENT COHERENTS ENTRE DSIGM ET DEPSM C ETAN = XKN/NCOEF*(SYCOE/SIGM)**(NCOEF-1.0D0) ESTA = 1.0D0/ETAN - 1.0D0/XKN ESTA = 1.0D0/ESTA C C CALCUL DU MULTIPLICATEUR DE LAGRANGE C DLAM2 = -CRIT2/(XKN*(DERIV21)**2.0D0- & DERIV23*KCOEF*FVAR*(1.0D0-FVAR)*DERIV21- & ESTA/((1.0D0-FVAR)*SIGM)*DERIV21*SIGR*DERIV22) C C ACTUALISATION DES VARIABLES FORCES C SIGR = SIGR + DLAM2 * XKN * DERIV21 FVAR = FVAR - DLAM2 * KCOEF * FVAR * (1.0D0-FVAR) * DERIV21 SIGM = SIGM - DLAM2 * SIGR * DERIV21 * & ESTA/(1.0D0-FVAR)/SIGM C C SECURITE SUR LA VALEUR DE FVAR (CAS OU ELLE DEVIENT TROP FAIBLE) C IF (FVAR.LE.FVARLIM) THEN FVAR = FVARLIM MODR = nint(DLAM2 * XKN * DERIV21/DEPST(2) ) END IF C C -------------------------------------------- C RETURN MAPPING : FIN DES ITERATIONS INTERNES C -------------------------------------------- C END DO C C BALISE DE SORTIE DU RETURN MAPPING AVANT ACUALISATION C 20 CONTINUE C C CAS OU LA FRACTION VOLUMIQUE DES PORES EST EGALE A 0.02 C IF ((FVAR.EQ.FVARLIM).AND. & (IFLAG.EQ.1)) THEN SIGR = SIGR - MODR * DEPST(3) END IF C C CALCUL DES QUANTITES RESULTATS (EPSR, VKACT ET XKN) C SIGN = -SIGR EPSR = -(EPSN - SIGR/XKN) C C BALISE DE SORTIE DU RETURN MAPPING APRES ACUALISATION C 30 CONTINUE C C ACTUALISATION DES PARAMETRES ELASTIQUES C XKNACT = 4.0D0*XKN0*XKS0*(1.0D0-FVAR)/ & (4.0D0*XKS0+3.0D0*XKN0*FVAR) XKSACT = XKS0*(1.0D0-FVAR)/ & (1.0D0+((6.0D0*XKN0+12.0D0*XKS0)/ & (9.0D0*XKN0+8.0D0*XKS0))*FVAR) XKN = XKNACT XKS = XKSACT C C-----CALCUL DES CONTRAINTES FINALES------------------------------------ C SIGF(1) = XKS*(1-D)*EPST1 + SPI1 SIGF(2) = XKS*(1-D)*EPST2 + SPI2 IF (EPSN.GT.0.0D0) THEN SIGF(3) = XKN*(1-D)*EPSN ELSE SIGF(3) = XKN*(EPSN-EPSR) ENDIF C C-----STOCKAGE EN SORTIE DES VARIABLES INTERNES------------------------- C VARF(1 ) = D VARF(2 ) = EPSF1 VARF(3 ) = EPSF2 VARF(4 ) = XEC1 VARF(5 ) = XEC2 VARF(9 ) = CIN1 VARF(10) = CIN2 VARF(11) = EPSR VARF(12) = IFLAG VARF(13) = FVAR VARF(14) = SIGM VARF(15) = XKS VARF(16) = XKN VARF(17) = COUN + 1 123 CONTINUE IF (NGONF.EQ.1) THEN SIGF(1) = XMAT(1)*EPST1 SIGF(2) = XMAT(1)*EPST2 SIGF(3) = XMAT(2)*EPSN ENDIF VARF(6 ) = EPST1 VARF(7 ) = EPST2 VARF(8 ) = EPSN RETURN END