C RICJ2 SOURCE PV 21/11/04 21:15:15 11174 C sub ricjoi2d C RICJOI2D SOURCE CHAT 05/01/12 22:24:09 5004 SUBROUTINE RICJ2(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI, & DEPST,IFOUR,XMAT,NMATT,IVAL,DD,SIGF,DEFP,VARF,KERRE) C----------------------------------------------------------------------- C Modele RICJOI en version element joint 2D 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 DN,NCOEF,KCOEF DN = 0.D0 FVARLIM = 2.0D-5 C C-----COMPTEUR POUR LE NOMBRE DE PASSAGE DAS LA BOUCLE C COUN = VAR0(13) C C-----ETAT PRECEDENT---------------------------------------------------- C D = VAR0(1) EPSF = VAR0(2) XEC = VAR0(3) EPST = VAR0(4) EPSN = VAR0(5) CINE = VAR0(6) EPSR = VAR0(7) IFLAG = nint(VAR0(8)) FVAR = VAR0(9) SIGM = VAR0(10) XKS = VAR0(11) XKN = VAR0(12) 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 APP2(TCCOR,YR,15) C C-----ACTUALISATION DE LA DEFORMATION----------------------------------- C EPST = EPST + DEPST(1) EPSN = EPSN + DEPST(2) IF (NGONF.EQ.1) THEN GOTO 123 ENDIF C C-----CALCUL DE L ENERGIE----------------------------------------------- C EPSNP = 0.5D0*(ABS(EPSN)+EPSN) YN = 0.5D0*XKN*EPSNP*EPSNP YD = 0.5D0*XKS*EPST*EPST C C-----CALCUL DU SEUIL D ENDOMMAGEMENT----------------------------------- C XFD = XAL*YN+YD-(XY0+XEC+YR) C C-----TEST SUR LE SEUIL D ENDOMMAGEMENT--------------------------------- C IF (XFD.GT.0.0D0) THEN Dn = 1.0D0-1.0D0/(1+XAD*(XAL*YN+YD-XY0-YR)) XEC = XAL*YN+YD-XY0-YR ENDIF IF (Dn.GT.D) THEN D = Dn ENDIF IF (D.GE.1.0D0) THEN D = 0.999D0 ENDIF C C-----SUBSTEPPING------------------------------------------------------- C NENTIER = 1 DO IQ = 1,NENTIER EPST = VAR0(4)+DEPST(1)*(IQ/NENTIER) EPSN = VAR0(5)+DEPST(2)*(IQ/NENTIER) XSTE = XKN*EPSN C C-----FROTTEMENT--------------------------------------------------------- C SPI = D*XKS*(EPST-EPSF) CRITF = ABS(SPI-CINE)+XDILAT*0.5D0*(-ABS(XSTE)+XSTE) IF (CRITF.GT.(0.0D0)) THEN XF0 = ABS(SPI-CINE)+XDILAT*0.5D0*(-ABS(XSTE)+XSTE) XDLAM = 1.0D0 DO I=1,500 SEUIL = (ABS(SPI-CINE)+ & XDILAT*0.5D0*(-ABS(XSTE)+XSTE))/XF0 IF ((SEUIL.LE.1.0D-5).OR.(XDLAM.LE.1.0D-10)) THEN GOTO 900 ELSE XSG = ABS(SPI-CINE)/(SPI-CINE) XDLAM = (ABS(SPI-CINE)+ & XDILAT*0.5D0*(-ABS(XSTE)+XSTE))/ & (XKS*D-XSG*XGA*(-XSG+XAA*CINE)) CINE = CINE-XGA*XDLAM*(-XSG+XAA*CINE) SPI = SPI-D*XKS*XDLAM*XSG ENDIF EPSF = -SPI/(D*XKS)+EPST ENDDO ENDIF 900 CONTINUE ENDDO 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 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 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(2)).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(2) 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)*EPST + SPI IF (EPSN.GT.0.0D0) THEN SIGF(2) = XKN*(1-D)*EPSN ELSE SIGF(2) = XKN*(EPSN-EPSR) ENDIF C C-----STOCKAGE EN SORTIE DES VARIABLES INTERNES------------------------- C VARF(1 ) = D VARF(2 ) = EPSF VARF(3 ) = XEC VARF(6 ) = CINE VARF(7 ) = EPSR VARF(8 ) = IFLAG VARF(9 ) = FVAR VARF(10) = SIGM VARF(11) = XKSACT VARF(12) = XKNACT VARF(13) = COUN + 1 123 CONTINUE IF (NGONF.EQ.1) THEN SIGF(1) = XMAT(1)*EPST SIGF(2) = XMAT(2)*EPSN ENDIF VARF(4 ) = EPST VARF(5 ) = EPSN RETURN END