C RIC2NL SOURCE PV 17/12/08 21:17:38 9660 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-----OUVERTURE DE SEGMENTS--------------------------------------------- C---------------------------------------------------------------------- C SEGMENT WRK0 C REAL*8 XMAT(NMATT) C ENDSEGMENT C C SEGMENT WRK1 C REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS) C REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI) C REAL*8 DEFP(NSTRS),XCAR(ICARA) C ENDSEGMENT C C SEGMENT WRK5 C REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS) C ENDSEGMENT C C SEGMENT WRK6 C REAL*8 SIG0S(NSTRS),DEPSTS(NSTRS) C END SEGMENT C CHARACTER*8 CMATE C INTEGER NSTRS,NVARI,NMATT,REMP,KIN,NENTIER,COCO 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 I=1,3 DO J=1,3 IF (I.EQ.J) THEN UNIT(I,J)=UN ELSE UNIT(I,J)=XZERO ENDIF ENDDO ENDDO 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**2.D0 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) IF (TRA.GE.XZERO) THEN c A1 = 7.0D-6 XPDC=1.0D0 ELSE c A1 = 7.0D-7 XPDC=XZERO ENDIF DO I=1,3 DO J=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 I=1,3 DO J=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 I=1,3 DO J=1,3 DONMDNL=DONMDNL+(K/6.0D0)*EPSDIR(I,J)*EPSDIR(I,J) DONMINL=DONMINL+(K/6.0D0)*EPSIND(I,J)*EPSIND(I,J) ENDDO ENDDO SEUIL1NL=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 IF (Dn.GE.D) THEN D=Dn ENDIF C C-----ON LIMITE LA VALEUR DE D C IF (D.GE.0.9999D0) THEN D=0.9999D0 ENDIF 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 I=1,3 DO J=1,3 X(I,J)=GAMMA1*ALPHA(I,J) EPSDEV(I,J)=EPSI(I,J)-(1.0D0/3.0D0)*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)/2.0D0 EPSI(1,3)=0.0D0 EPSI(2,3)=0.0D0 EPSI(2,1)=EPSILO(4)/2.0D0 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) TRAPI = EPSIPI(1,1)+EPSIPI(2,2)+EPSIPI(3,3) DO I=1,3 DO J=1,3 EPSDEV(I,J)=EPSI(I,J)-(1.0D0/3.0D0)*TRA*UNIT(I,J) DEVPI(I,J)=EPSIPI(I,J)-(1.0D0/3.0D0)*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 I=1,3 DO J=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 I=1,3 DO 201 J=1,3 SIGMOX(I,J)=XZERO SIGMOX(I,J)=SIGMPI(I,J)-X(I,J) 201 CONTINUE 101 CONTINUE TRASX=XZERO TRASX=SIGMOX(1,1)+SIGMOX(2,2)+SIGMOX(3,3) DO 102 I=1,3 DO 202 J=1,3 SIGMOD(I,J)=XZERO SIGMOD(I,J)=SIGMOX(I,J)-(1.0D0/3.0D0)*TRASX*UNIT(I,J) 202 CONTINUE 102 CONTINUE TERME4=XZERO DO 103 I=1,3 DO 203 J=1,3 TERME4=TERME4+SIGMOD(I,J)*SIGMOD(I,J) 203 CONTINUE 103 CONTINUE C C CALCUL DU J2(SIGMPI-X) C TERME4=SQRT(3.D0/2.D0*TERME4) SEUIL2=TERME4+XCOEF*(1.0D0/3.0D0)*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 I=1,3 DO J=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 I=1,3 DO J=1,3 TERME6(I,J)=(3.D0/2.D0)*SIGMOX(I,J)/TERME4+XCOEF*1.0D0/3.D0* & UNIT(I,J) ENDDO ENDDO C C CALCUL DE DF/DX.................................................... C DO I=1,3 DO J=1,3 TERME7(I,J)=-TERME6(I,J) ENDDO ENDDO C C CALCUL DE DG/DS.................................................... C DO I=1,3 DO J=1,3 TERME8(I,J)=(3.D0/2.D0)*SIGMOX(I,J)/TERME4+XCOEF*1.0D0/3.D0* & UNIT(I,J) ENDDO ENDDO C C CALCUL DE DG/DX.................................................... C DO I=1,3 DO J=1,3 TERME9(I,J)=-(3.D0/2.D0)*SIGMOX(I,J)/TERME4+A1*X(I,J)- & XCOEF*1.0D0/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) DO I=1,3 DO J=1,3 DEVDGS(I,J)=TERME8(I,J)-(1.0D0/3.0D0)*TRDGZ*UNIT(I,J) ENDDO ENDDO TERM10=XZERO DO I=1,3 DO J=1,3 TERM10=TERM10+TERME6(I,J)*((K/3.0D0)*TRDGZ*UNIT(I,J)+ & 2.0D0*G*DEVDGS(I,J)) ENDDO ENDDO TERM11=XZERO DO I=1,3 DO J=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 I=1,3 DO J=1,3 SIGMPI(I,J)=SIGMPI(I,J)-DLAM2*D*((K/3.0D0)*TRDGZ*UNIT(I,J)+ & 2.0D0*G*DEVDGS(I,J)) ENDDO ENDDO DO I=1,3 DO J=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 I=1,3 DO 403 J=1,3 ALPHA(I,J)=XZERO ALPHA(I,J)=X(I,J)/GAMMA1 403 CONTINUE 303 CONTINUE C C 2) EPSILON_PI C ---------- TROGO = TRA-(SIGMPI(1,1)+SIGMPI(2,2)+SIGMPI(3,3))/(D*K) DO I=1,3 DO J=1,3 DEVTRO(I,J)=EPSDEV(I,J)-(SIGMPI(I,J)-(1.0D0/3.0D0)* & (SIGMPI(1,1)+SIGMPI(2,2)+SIGMPI(3,3))*UNIT(I,J))/(2.0D0*G*D) ENDDO ENDDO DO I=1,3 DO J=1,3 EPSIPI(I,J)=(1.0D0/3.0D0)*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)/2.0D0 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