C RIC2NL SOURCE OF166741 25/11/04 21:16:03 12349 C sub ricnl2D C====&===1=========2=========3=========4=========5=========6=========7== C Commentaires : Subroutine permettant de mettre en oeuvre le C modele RICRAG en 2D et en non local C Traits : - Endommagement anisotrope C - Monotone C Auteurs : B. Richard (doctorant) C Date : Février 2008 C====&===1=========2=========3=========4=========5=========6=========7== SUBROUTINE RIC2NL(wrk52,wrk53,wrk54,nvari,iecou) C---------------------------------------------------------------------- C-----DECLARATION GENERALE---------------------------------------------- C---------------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) C---------------------------------------------------------------------- C-----APPEL AUX LIBRAIRIES---------------------------------------------- C---------------------------------------------------------------------- -INC PPARAM -INC CCOPTIO -INC DECHE C---------------------------------------------------------------------- C-----DECLARATION PARTICULIERES----------------------------------------- C---------------------------------------------------------------------- REAL*8 UNIT(3,3),E,NU,FT,K,ALIND,ALDIR,GAMMA1,A1,D,Dn REAL*8 EPSI(3,3),EPSIPI(3,3),ALPHA(3,3),DEPSI(3,3) REAL*8 XLAM,XMU,EPSDIR(3,3),EPSIND(3,3),FFFZ,Y0,EPSILO(4) REAL*8 SIGMA(3,3),VPEPS(3),X(3,3),FORCED,FORCEI,EPSILI(4) REAL*8 EPSILB(4),SIGMP(4),EPSILC(4),EPSILP(4),SIGMPI(3,3) REAL*8 SIGMOX(3,3),SIGMOD(3,3),TERME4,TERME6(3,3),TERME7(3,3) REAL*8 TERME8(3,3),TERME9(3,3),DEVDGS(3,3),TERM11,DLAM2 REAL*8 DEVTRO(3,3),EPSDEV(3,3),DEVPI(3,3),TRA,TRAPI,TRASX REAL*8 SEUINI,SEUIL2,TRDGZ,TERM10,TROGO,XPDC C---------------------------------------------------------------------- C-----PARAMETRES DE CALCUL---------------------------------------------- C---------------------------------------------------------------------- PARAMETER (XZERO=0.0D0 , UN=1.0D0 , DEUX=2.0D0, XCOEF=0.09D0) DO J=1,3 DO I=1,3 UNIT(I,J)=XZERO ENDDO ENDDO UNIT(1,1)=UN UNIT(2,2)=UN UNIT(3,3)=UN C---------------------------------------------------------------------- C-----LES VARIABLES INTERNES D ENTREE----------------------------------- C---------------------------------------------------------------------- C C-----GLISSEMENT (EPSILON_PI) C EPSIPI(1,1)=VAR0(2) EPSIPI(2,2)=VAR0(3) EPSIPI(3,3)=VAR0(4) EPSIPI(1,2)=VAR0(5) EPSIPI(1,3)=VAR0(6) EPSIPI(2,3)=VAR0(7) EPSIPI(2,1)=VAR0(5) EPSIPI(3,1)=VAR0(6) EPSIPI(3,2)=VAR0(7) C C-----ECROUISSAGE CINEMATIQUE (ALPHA) C ALPHA(1,1)=VAR0(8) ALPHA(2,2)=VAR0(9) ALPHA(3,3)=VAR0(10) C ALPHA(1,2)=VAR0(11) ALPHA(1,3)=VAR0(12) ALPHA(2,3)=VAR0(13) C ALPHA(2,1)=VAR0(11) ALPHA(3,1)=VAR0(12) ALPHA(3,2)=VAR0(13) C C-----ENDOMMAGEMENT (D) C D=VAR0(14) C C-----FORCE ECROUISSAGE ISOTROPE C FFFZ=VAR0(17) C C-----DEFORMATION TOTALES C DO I=1,4 EPSILO(I)=VAR0(17+I)+DEPST(I) EPSILI(I)=VAR0(17+I) ENDDO C C-----ECCROUISSAGE C FORCED=VAR0(25) FORCEI=VAR0(26) C---------------------------------------------------------------------- C-----RENOMMAGE DES CONSTANTES ELASTIQUES------------------------------- C---------------------------------------------------------------------- E = XMAT(1) NU = XMAT(2) FT = XMAT(5) ALIND = XMAT(6) GAMMA1 = XMAT(7) A1 = XMAT(8) ALDIR = XMAT(9) EPS0 = FT/E K = E/(1.D0-2.D0*NU) G = E/(2.D0*(1.D0+NU)) Y0 = K/6.0D0*(EPS0*EPS0) C---------------------------------------------------------------------- C-----INDICATEURS POUR LES OPTIONS DE CALCUL---------------------------- C---------------------------------------------------------------------- NSTRSV = NSTRS IFOUR2 = IFOUR JDIM = IDIM C---------------------------------------------------------------------- C-----CALCUL DE LA CONTRAINTE ELASTIQUE--------------------------------- C---------------------------------------------------------------------- CMATE = 'ISOTROPE' KCAS=1 CALL DOHMAS (XMAT,CMATE,IFOUR2,NSTRSV,KCAS,DDHOOK,IRTD) CALL MATVE1 (DDHOOK,DEPST,NSTRSV,NSTRSV,DSIGT,1) PREC=1.D-08 CALL DOHMAS (XMAT,CMATE,IFOUR2,NSTRSV,2,DDHOOK,IRTD) CALL INVALM (DDHOOK,NSTRSV,NSTRSV,IRTD,PREC) CALL MATVE1 (DDHOOK,SIG0,NSTRSV,NSTRSV,EPSILI,1) CALL ENDOCB (EPSILO,EPSI,2,IFOUR2) CALL ENDOCB (DEPST,DEPSI,2,IFOUR2) KCAS=1 CALL DOHMAS (XMAT,CMATE,IFOUR2,NSTRSV,KCAS,DDHOOK,IRTD) CALL MATVE1 (DDHOOK,EPSILO,NSTRSV,NSTRSV,SIGF,2) SIGMA(1,1) = SIGF(1) SIGMA(1,2) = SIGF(4) SIGMA(1,3) = 0.0D0 SIGMA(2,1) = SIGF(4) SIGMA(2,2) = SIGF(2) SIGMA(2,3) = 0.0D0 SIGMA(3,1) = 0.0D0 SIGMA(3,2) = 0.0D0 SIGMA(3,3) = SIGF(3) C---------------------------------------------------------------------- C-----INTEGRATION DE L ENDOMMAGEMENT----------------------------------- C---------------------------------------------------------------------- C C-----CALCUL DES VALEURS PRINCIPALES C CALL JACOD3(EPSI,2,VPEPS) C C-----CALCUL DE EPSILON DIRECT C TRA = (EPSI(1,1) + EPSI(2,2) + EPSI(3,3))/3.D0 IF (TRA.GE.XZERO) THEN c A1 = 7.0D-6 XPDC=1.0D0 ELSE c A1 = 7.0D-7 XPDC=XZERO ENDIF DO J=1,3 DO I=1,3 IF (I.EQ.J) THEN EPSDIR(I,J)=0.5D0*(ABS(VPEPS(I))+VPEPS(I))*XPDC ELSE EPSDIR(I,J)=XZERO ENDIF ENDDO ENDDO C C-----CALCUL DE EPSILON INDUIT C DO J=1,3 DO I=1,3 IF (I.EQ.J) THEN EPSIND(I,J)=0.5D0*(ABS(VPEPS(I))+VPEPS(I))-EPSDIR(I,J) ELSE EPSIND(I,J)=XZERO ENDIF ENDDO ENDDO C C-----TAUX D ENERGIES CORRESPONDANTES C DONMDNL=XZERO DONMINL=XZERO DO J=1,3 DO I=1,3 DONMDNL=DONMDNL+EPSDIR(I,J)*EPSDIR(I,J) DONMINL=DONMINL+EPSIND(I,J)*EPSIND(I,J) ENDDO ENDDO SEUIL1NL= (K/6.0D0) * (DONMDNL*XPDC+DONMINL*(1.0D0-XPDC)) C C-----LE NON LOCAL C IF (ISTEP.EQ.0) THEN SEUIL1=SEUIL1NL VARF(1)=SEUIL1 ELSE IF (ISTEP.EQ.1) THEN VARF(1)=SEUIL1NL DO I=2,26 VARF(I)=VAR0(I) ENDDO GOTO 2000 ELSE IF (ISTEP.EQ.2) THEN SEUIL1=VAR0(1) VARF(1)=SEUIL1 ENDIF FFFZ=XZERO FFFZ=FORCED*XPDC+FORCEI*(1.0D0-XPDC) SEUILT=SEUIL1-(FFFZ+Y0) C C-----EVOLUTION ENDOMMAGEMENT C IF (SEUILT.GT.XZERO) THEN IF (TRA.GE.XZERO) THEN Dn = 1.0D0-1.0D0/(1.0D0+ALDIR*(SEUIL1-Y0)) FORCED = SEUIL1-Y0 ELSE Dn = 1.0D0-1.0D0/(1.0D0+ALIND*(SEUIL1-Y0)) FORCEI = SEUIL1-Y0 ENDIF C C-----ON VERIFIE LA CROISSANCE DE D C D = MAX(D,Dn) C C-----ON LIMITE LA VALEUR DE D C D = MIN(D,0.9999D0) C C-----ON BLOQUE LES ECROUISSAGES C FFFZ=FORCED*XPDC+FORCEI*(1.0D0-XPDC) ENDIF C C---------------------------------------------------------------------- C-----INTEGRATION DU FROTTEMENT----------------------------------------- C---------------------------------------------------------------------- C C-----CALCUL DE X ET DU DEVIATEUR DE EPSILON TOTAL C DO J=1,3 DO I=1,3 X(I,J)=GAMMA1*ALPHA(I,J) EPSDEV(I,J)=EPSI(I,J)-TRA*UNIT(I,J) ENDDO ENDDO C C-----PREPARATION AU SUBSTEPPING C NENTIER=1 DO KIN=1,NENTIER DO I=1,4 EPSILO(I)=VAR0(17+I)+KIN*DEPST(I)/NENTIER ENDDO C C-----ON PASSE EPSILO DANS EPSI C EPSI(1,1)=EPSILO(1) EPSI(2,2)=EPSILO(2) EPSI(3,3)=EPSILO(3) EPSI(1,2)=EPSILO(4)*0.5D0 EPSI(1,3)=0.0D0 EPSI(2,3)=0.0D0 EPSI(2,1)=EPSILO(4)*0.5D0 EPSI(3,1)=0.0D0 EPSI(3,2)=0.0D0 C C-----ON CALCUL LA TRACE DE EPSILON ET DE EPSILON_PI C TRA = (EPSI(1,1)+EPSI(2,2)+EPSI(3,3))/3.D0 TRAPI = (EPSIPI(1,1)+EPSIPI(2,2)+EPSIPI(3,3))/3.D0 DO J=1,3 DO I=1,3 EPSDEV(I,J)=EPSI(I,J)-TRA*UNIT(I,J) DEVPI(I,J)=EPSIPI(I,J)-TRAPI*UNIT(I,J) ENDDO ENDDO C C-----ON PASSE EPSILON_PI DANS EPSILP C EPSILP(1)=EPSIPI(1,1) EPSILP(2)=EPSIPI(2,2) EPSILP(3)=EPSIPI(3,3) EPSILP(4)=2.0D0*EPSIPI(1,2) C C-----ON FORME LES DEFORMATIONS ASSOCIEES A SIGMA_PI C DO I=1,4 EPSILC(I)=(EPSILO(I)-EPSILP(I))*D ENDDO C C-----CALCUL DE SIGMPI C KCAS=1 CALL DOHMAS (XMAT,CMATE,IFOUR2,NSTRSV,KCAS,DDHOOK,IRTD) CALL MATVE1 (DDHOOK,EPSILC,NSTRSV,NSTRSV,SIGMP,2) C C-----ON PASSE SIGMP DANS SIGMPI C DO J=1,3 DO I=1,3 SIGMPI(I,J)=XZERO ENDDO ENDDO SIGMPI(1,1)=SIGMP(1) SIGMPI(2,2)=SIGMP(2) SIGMPI(3,3)=SIGMP(3) SIGMPI(1,2)=SIGMP(4) SIGMPI(2,1)=SIGMP(4) C C-----ON FORME LES DEFORMATIONS ASSOCIEES A SIGMA C DO I=1,4 EPSILB(I)=(EPSILO(I)-EPSILP(I)*D) ENDDO C C-----ON CALCUL SIGMA EN VECTEUR STOCKEE DANS SIGF C KCAS=1 CALL DOHMAS (XMAT,CMATE,IFOUR2,NSTRSV,KCAS,DDHOOK,IRTD) CALL MATVE1 (DDHOOK,EPSILB,NSTRSV,NSTRSV,SIGF,2) C C-----ON PASSE SIGF DANS SIGMA C SIGMA(1,1)=SIGF(1) SIGMA(2,2)=SIGF(2) SIGMA(3,3)=SIGF(3) SIGMA(1,2)=SIGF(4) SIGMA(2,1)=SIGF(4) C C-----DEBUT DES ITERATIONS INTERNES C DO 500 IREMP=1,500 C C 1) EVALUATION DU SEUIL C ---------------------- DO 101 J=1,3 DO 201 I=1,3 SIGMOX(I,J)=SIGMPI(I,J)-X(I,J) 201 CONTINUE 101 CONTINUE TRASX=(SIGMOX(1,1)+SIGMOX(2,2)+SIGMOX(3,3))/3.D0 DO 102 J=1,3 DO 202 I=1,3 SIGMOD(I,J)=SIGMOX(I,J)-TRASX*UNIT(I,J) 202 CONTINUE 102 CONTINUE TERME4=XZERO DO 103 J=1,3 DO 203 I=1,3 TERME4=TERME4+SIGMOD(I,J)*SIGMOD(I,J) 203 CONTINUE 103 CONTINUE C C CALCUL DU J2(SIGMPI-X) C TERME4=SQRT(1.5D0*TERME4) SEUIL2=TERME4+XCOEF*TRASX IF (IREMP.EQ.1) THEN SEUINI=SEUIL2 DLAM2=XZERO ENDIF IF (SEUINI.GT.1.0D0) THEN CRIT=SEUIL2/SEUINI ELSE CRIT=XZERO ENDIF IF (IREMP.EQ.1) ICOCO=2 C C 2) ON VERIFIE LA VALEUR RELATIVE DU SEUIL C ----------------------------------------- IF ((CRIT.LE.1.0D-5).OR.(DLAM2.LE.1.0D-10.AND.IREMP.GT.1).OR. & (D.EQ.XZERO)) THEN GOTO 666 ELSE IF (IREMP.EQ.1) ICOCO=1 C C 3) CALCUL DES DERIVÉES POUR LA REDESCENTE AU GRADIENT C ----------------------------------------------------- DO J=1,3 DO I=1,3 TERME6(I,J)=XZERO TERME7(I,J)=XZERO TERME8(I,J)=XZERO TERME9(I,J)=XZERO ENDDO ENDDO C C CALCUL DE DF/DS.................................................... C DO J=1,3 DO I=1,3 TERME6(I,J)=1.5D0*SIGMOX(I,J)/TERME4+XCOEF/3.D0*UNIT(I,J) ENDDO ENDDO C C CALCUL DE DF/DX.................................................... C DO J=1,3 DO I=1,3 TERME7(I,J)=-TERME6(I,J) ENDDO ENDDO C C CALCUL DE DG/DS.................................................... C DO J=1,3 DO I=1,3 TERME8(I,J)=1.5D0*SIGMOX(I,J)/TERME4+XCOEF/3.D0*UNIT(I,J) ENDDO ENDDO C C CALCUL DE DG/DX.................................................... C DO J=1,3 DO I=1,3 TERME9(I,J)=-1.5D0*SIGMOX(I,J)/TERME4+A1*X(I,J) & -XCOEF/3.D0*UNIT(I,J) ENDDO ENDDO C C 4) CALCUL DU MULTIPLICATEUR DE GLISSEMENT C ----------------------------------------- TRDGZ=(TERME8(1,1)+TERME8(2,2)+TERME8(3,3))/3.D0 DO J=1,3 DO I=1,3 DEVDGS(I,J)=TERME8(I,J)-TRDGZ*UNIT(I,J) ENDDO ENDDO TERM10=XZERO DO J=1,3 DO I=1,3 TERM10=TERM10+TERME6(I,J)*(K*TRDGZ*UNIT(I,J)+ & 2.0D0*G*DEVDGS(I,J)) ENDDO ENDDO TERM11=XZERO DO J=1,3 DO I=1,3 TERM11=TERM11+TERME7(I,J)*TERME9(I,J) ENDDO ENDDO DLAM2=SEUIL2/(TERM10+GAMMA1*TERM11) C C 5) ACTUALISATION DES VARIABLES FORCES C ------------------------------------- DO J=1,3 DO I=1,3 SIGMPI(I,J)=SIGMPI(I,J)-DLAM2*D*(K*TRDGZ*UNIT(I,J)+ & 2.0D0*G*DEVDGS(I,J)) ENDDO ENDDO DO J=1,3 DO I=1,3 X(I,J)=X(I,J)-GAMMA1*DLAM2*TERME9(I,J) ENDDO ENDDO ENDIF 500 CONTINUE C C-----ACTUALISATION DES VARIABLES FLUX (S'IL Y A LIEU) C C 1) ALPHA C ------ 666 CONTINUE IF (ICOCO.EQ.1) THEN DO 303 J=1,3 DO 403 I=1,3 ALPHA(I,J)=X(I,J)/GAMMA1 403 CONTINUE 303 CONTINUE C C 2) EPSILON_PI C ---------- TRA_Z = (SIGMPI(1,1)+SIGMPI(2,2)+SIGMPI(3,3))/3.D0 TROGO = TRA-(TRA_Z/(D*K)) DO J=1,3 DO I=1,3 DEVTRO(I,J)= EPSDEV(I,J) & -(SIGMPI(I,J)-TRA_Z*UNIT(I,J))/(2.0D0*G*D) ENDDO ENDDO DO J=1,3 DO I=1,3 EPSIPI(I,J)=TROGO*UNIT(I,J)+DEVTRO(I,J) ENDDO ENDDO EPSILP(1)=EPSIPI(1,1) EPSILP(2)=EPSIPI(2,2) EPSILP(3)=EPSIPI(3,3) EPSILP(4)=2.0D0*EPSIPI(1,2) ENDIF ENDDO C C-----ACTUALISATION DES CONTRAINTES C 777 CONTINUE DO I=1,4 EPSILB(I)=(EPSILO(I)-D*EPSILP(I)) ENDDO CALL DOHMAS (XMAT,CMATE,IFOUR2,NSTRSV,KCAS,DDHOOK,IRTD) CALL MATVE1 (DDHOOK,EPSILB,NSTRSV,NSTRSV,SIGF,2) C C-----MISE DE SIGF DANS SIGMA C SIGMA(1,1) = SIGF(1) SIGMA(1,2) = SIGF(4) SIGMA(1,3) = 0.0D0 SIGMA(2,1) = SIGF(4) SIGMA(2,2) = SIGF(2) SIGMA(2,3) = 0.0D0 SIGMA(3,1) = 0.0D0 SIGMA(3,2) = 0.0D0 SIGMA(3,3) = SIGF(3) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C STOCKAGE EN SORTIE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C VECTEUR DES CONTRAINTES C SIGF(1) = SIGMA(1,1) SIGF(2) = SIGMA(2,2) SIGF(3) = SIGMA(3,3) SIGF(4) = SIGMA(1,2) C C GLISSEMENT (EPSILON_PI) C VARF(2)=EPSILP(1) VARF(3)=EPSILP(2) VARF(4)=EPSILP(3) VARF(5)=EPSILP(4)*0.5D0 VARF(6)=0.0D0 VARF(7)=0.0D0 C C ECROUISSAGE CINEMATIQUE (ALPHA) C VARF(8 )=ALPHA(1,1) VARF(9 )=ALPHA(2,2) VARF(10)=ALPHA(3,3) VARF(11)=ALPHA(1,2) VARF(12)=ALPHA(1,3) VARF(13)=ALPHA(2,3) C C ENDOMMAGEMENT (D) C VARF(14)=D C C FORCE LIEE A L ECROUISSAGE ISOTROPE C VARF(17)=FFFZ C C DEFORMATION TOTALE C DO I=1,4 VARF(I+17)=EPSILO(I) ENDDO C C ECCROUISSAGE C VARF(25)=FORCED VARF(26)=FORCEI 2000 CONTINUE RETURN END