C DESMNL SOURCE CB215821 16/04/21 21:16:24 8920 C DESMORAT SUBROUTINE DESMNL(WRK0,WRK1,WRK5,NSTRS,NVARI,LHOOK,ISTEP, & NMATT,ICARA,KERRE,MFR) C C variables en entree C C WRK0, WRK1 pointeur actif a l'entree C NVARI nombre de variables internes (doit etre egal a ?) C NMATT nombre de constantes du materiau C C ISTEP flag utilise pour separer les etapes dans un calcul non local C ISTEP=0 -----> calcul local C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils C ISTEP=2 -----> calcul non local etape 2 on continue le calcul C a partir des seuils moyennes C C variables en sortie C C VARINF variables internes finales C C SIGMAF contraintes finales C C IMPLICIT REAL*8(A-H,O-Z) C -INC PPARAM -INC CCOPTIO -INC SMINTE -INC SMCOORD -INC SMELEME C C CHARACTER*8 CMATE REAL*8 XD0(3,3),DF(3,3),POSEPSI2(3,3),EPSIN(6),EPS33(3,3) REAL*8 SIGTILDE33(3,3),SIGDEV33(3,3),SIGD33(3,3),SIGTILDE(6) REAL*8 UNMOINSDVP(3),HM1(3,3),DFD(3,3),UNMOINSD(3,3),DFVP(3) REAL*8 ROT(3,3),ROTT(3,3),VECP33T(3,3),TRAV(3,3),VECP33(3,3) REAL*8 EPSIPOS(3),HM1DIAG(3),EPS(2,2),DELTAF22(2,2),EPSIPP(3) REAL*8 VPDF(2),VECPDF(2,2),EPSTF(6),DHOOK(6,6) C PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12) C SEGMENT WRK0 REAL*8 XMAT(NMATT) ENDSEGMENT C SEGMENT WRK1 REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS) REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI) REAL*8 DEFP(NSTRS),XCAR(ICARA) ENDSEGMENT C SEGMENT WRK5 REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS) ENDSEGMENT C C C Recuperation des variables initiales dans les tableaux C CMATE= 'ISOTROPE' YOUN = XMAT(1) XNU = XMAT(2) XK0 = XMAT(5) AA = XMAT(6) A = XMAT(7) etaC = XMAT(8) etaT = XMAT(9) Xk = XMAT(10) DC=XMAT(11) C C Recuperation de l'endommagement initial C XD0(1,1)=VAR0(2) XD0(2,2)=VAR0(3) XD0(3,3)=VAR0(4) XD0(1,2)=VAR0(5) XD0(1,3)=VAR0(6) XD0(2,3)=VAR0(7) XD0(2,1)=XD0(1,2) XD0(3,1)=XD0(1,3) XD0(3,2)=XD0(2,3) C C Recuperation de la base d'endommagement C ROT(1,1)=VAR0(8) ROT(1,2)=VAR0(9) ROT(1,3)=VAR0(10) ROT(2,1)=VAR0(11) ROT(2,2)=VAR0(12) ROT(2,3)=VAR0(13) ROT(3,1)=VAR0(14) ROT(3,2)=VAR0(15) ROT(3,3)=VAR0(16) C C Recuperation de l'etat d'endommagement critique C IFAILURE=nint(VAR0(17)) C C Matrice de raideur elastique CALL ZERO(DDHOOK,6,6) AUX0=YOUN/((UN+XNU)*(UN-XNU-XNU)) AUX=AUX0*(UN-XNU) AUX1=AUX0*XNU GEGE=0.5D0*YOUN/(UN+XNU) DDHOOK(1,1)=AUX DDHOOK(2,1)=AUX1 DDHOOK(3,1)=AUX1 DDHOOK(1,2)=AUX1 DDHOOK(2,2)=AUX DDHOOK(3,2)=AUX1 DDHOOK(1,3)=AUX1 DDHOOK(2,3)=AUX1 DDHOOK(3,3)=AUX DDHOOK(4,4)=GEGE DDHOOK(5,5)=GEGE DDHOOK(6,6)=GEGE C C C Deformation totale C DO ISTRS=1,NSTRS EPSTF(ISTRS)=VAR0(17+ISTRS)+DEPST(ISTRS) ENDDO C - Trace de la deformation TRAEPS=EPSTF(1)+EPSTF(2)+EPSTF(3) C C calcul des deformations principales CALL ENDOCB(EPSTF,EPS33,2,IFOUR) CALL JACOB3(EPS33,3,EPSIPP,VECP33) EPSIPOS(1)= MAX( XZERO , EPSIPP(1) )**2 EPSIPOS(2)= MAX( XZERO , EPSIPP(2) )**2 EPSIPOS(3)= MAX( XZERO , EPSIPP(3) )**2 EPSTILNLOC=SQRT(EPSIPOS(1)+EPSIPOS(2)+EPSIPOS(3)) C IF (ISTEP.EQ.0) THEN EPSTIL=EPSTILNLOC VARF(1)=EPSTIL ELSE IF (ISTEP.EQ.1) THEN VARF(1)=EPSTILNLOC DO I=2,23 VARF(I)=VAR0(I) ENDDO GOTO 2000 ELSE IF (ISTEP.EQ.2) THEN c Over nonlocal, coef 2 C EPSTIL=DEUX*VAR0(1)-EPSTILNLOC EPSTIL=VAR0(1) ENDIF C C********************************************************** C Début de l'intégration numérique du comportement C********************************************************** C C calcul des seuils C C partie positive de la matrice des deformations CALL PRODT2(POSEPSI2,EPSIPOS,VECP33,3) C C calcul de la trace de D au pas n TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) C C calcul du seuil XKAPPA=A*TAN(TRAD0/(A*AA)+ATAN(XK0/A)) FCRIT=EPSTIL+Xk*TRAEPS-XKAPPA C C TEST sur le seuil IF(FCRIT.LT.XZERO) THEN C on est elastique C endommagement et deformations permanentes inchanges DO I=1,3 DO J=1,3 DF(I,J)=XD0(I,J) ENDDO ENDDO DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0 ELSE C C debut de la correction d'endom. C C Calcul de Kappa^(-1)=Tr Dn+1 C TRADF=A*AA*(ATAN((EPSTIL+Xk*TRAEPS)/A)-ATAN(XK0/A)) C C test sur l'eventuelle presence d'une direction d'endommagement fixee C IF(IFAILURE.EQ.0)THEN C actualisation de l'endommagement, Dc pas encore atteint (cas "normal") IF (EPSTIL.EQ.XZERO) THEN DLAMBDA=XZERO ELSE DLAMBDA=(TRADF-TRAD0)/(EPSTIL**2) ENDIF IF(DLAMBDA.LT.XZERO) THEN DLAMBDA=XZERO ENDIF CALL ZERO(DF,3,3) DO I=1,3 DO J=1,3 DF(I,J)=XD0(I,J)+DLAMBDA*POSEPSI2(I,J) ENDDO ENDDO C On borne l'endommagement D(1) a Dc si D>Dc et Fixed-crack CALL JACOD3(DF,3,DFVP) IF(DFVP(1).GE.DC) THEN IFAILURE=1 C mise a jour de XD0 pour traiter IFAILURE=1 dans le meme pas DO I=1,3 DO J=1,3 XD0(I,J)=DF(I,J) ENDDO ENDDO TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3) CALL JACOB3(XD0,3,DFVP,ROT) ENDIF ENDIF C IF(IFAILURE.EQ.1)THEN C une direction d'endom est fixee a Dc C C On calcule le multiplicateur d'endommagement dans la bonne base C On met les deformations positives dans le repère principal C de l'endommagement CALL JACOD3(XD0,3,DFVP) CALL TRSPOD(ROT,3,3,ROTT) CALL MATMAT(POSEPSI2,ROT,3,3,3,TRAV) CALL MATMAT(ROTT,TRAV,3,3,3,POSEPSI2) C IF(ABS(POSEPSI2(2,2)+POSEPSI2(3,3)).GT.XPETIT) THEN DLAMBDA=(TRADF-TRAD0)/(POSEPSI2(2,2)+POSEPSI2(3,3)) ELSE DLAMBDA=XZERO ENDIF IF(DLAMBDA.LT.XZERO) THEN DLAMBDA=XZERO ENDIF C C On corrige l'endommagement CALL ZERO(DF,3,3) DF(1,1)=Dc DF(2,2)=DFVP(2)+DLAMBDA*POSEPSI2(2,2) DF(3,3)=DFVP(3)+DLAMBDA*POSEPSI2(3,3) DF(2,3)=DLAMBDA*POSEPSI2(2,3) DF(3,2)=DF(2,3) C CALL MATMAT(DF,ROTT,3,3,3,TRAV) CALL MATMAT(ROT,TRAV,3,3,3,DF) C C On verifie la deuxieme direction d'endommagement CALL JACOD3(DF,3,DFVP) IF(DFVP(2).GE.Dc)THEN IFAILURE=2 C mise a jour de XD0 pour traiter IFAILURE=2 dans le meme pas DO I=1,3 DO J=1,3 XD0(I,J)=DF(I,J) ENDDO ENDDO TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3) CALL JACOB3(XD0,3,DFVP,ROT) ENDIF ENDIF C IF(IFAILURE.EQ.2) THEN C deux directions d'endom sont fixee a Dc CALL TRSPOD(ROT,3,3,ROTT) C DF33=TRADF-DEUX*Dc C print*, 'DF33=',DF33 IF(DF33.LT.XZERO) THEN DF33=XPETIT ENDIF C On corrige l'endommagement CALL ZERO(DF,3,3) DF(1,1)=Dc DF(2,2)=Dc DF(3,3)=DF33 C CALL MATMAT(DF,ROTT,3,3,3,TRAV) CALL MATMAT(ROT,TRAV,3,3,3,DF) C C On verifie la troisieme direction d'endommagement IF(DF33.GE.Dc)THEN IFAILURE=3 C mise a jour de XD0 pour traiter IFAILURE=3 dans le même pas DO I=1,3 DO J=1,3 XD0(I,J)=DF(I,J) ENDDO ENDDO TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) CALL JACOB3(XD0,3,DFVP,ROT) ENDIF ENDIF C IF(IFAILURE.EQ.3)THEN C On est totalement rompu, endommagement isotrope egal a DC CALL TRSPOD(ROT,3,3,ROTT) C On corrige l'endommagement CALL ZERO(DF,3,3) DF(1,1)=Dc DF(2,2)=Dc DF(3,3)=Dc C CALL MATMAT(DF,ROTT,3,3,3,TRAV) CALL MATMAT(ROT,TRAV,3,3,3,DF) ENDIF C C Fin de la correction d'endommagement C C DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0 C ENDIF C C Fin des differents cas (elasticite ou endommagement) C C Calcul de la contrainte finale C DO I=1,3 DO J=1,3 UNMOINSD(I,J)=-UN*DF(I,J) ENDDO ENDDO UNMOINSD(1,1)=UN-DF(1,1) UNMOINSD(2,2)=UN-DF(2,2) UNMOINSD(3,3)=UN-DF(3,3) CALL JACOB3(UNMOINSD,3,UNMOINSDVP,VECP33) HM1DIAG(1)=SQRT(UNMOINSDVP(1)) HM1DIAG(2)=SQRT(UNMOINSDVP(2)) HM1DIAG(3)=SQRT(UNMOINSDVP(3)) CALL PRODT2(HM1,HM1DIAG,VECP33,3) C C contrainte effective CALL MATVE1(DDHOOK,EPSTF,NSTRS,NSTRS,SIGTILDE,2) CALL ENDOCB(SIGTILDE,SIGTILDE33,1,IFOUR) CALL MATMAT(HM1,SIGTILDE33,3,3,3,TRAV) CALL MATMAT(TRAV,HM1,3,3,3,SIGD33) XTRAV=XZERO DO I=1,3 DO J=1,3 XTRAV=XTRAV+UNMOINSD(I,J)*SIGTILDE33(I,J) ENDDO ENDDO DO I=1,3 DO J=1,3 SIGD33(I,J)=SIGD33(I,J)-(XTRAV/(3.D0*(UN-DH)))*UNMOINSD(I,J) ENDDO ENDDO TRASIGTILDE=SIGTILDE(1)+SIGTILDE(2)+SIGTILDE(3) C IF(TRASIGTILDE.GE.XZERO) THEN COEF=UN-etaT*DH IF(COEF.LE.XZERO) THEN COEF=UN-Dc ENDIF TRASIG=COEF*TRASIGTILDE ELSE TRASIG=(UN-etaC*DH)*TRASIGTILDE ENDIF C SIGF(1)=SIGD33(1,1)+TRASIG/3.d0 SIGF(2)=SIGD33(2,2)+TRASIG/3.d0 SIGF(3)=SIGD33(3,3)+TRASIG/3.d0 SIGF(4)=SIGD33(1,2) SIGF(5)=SIGD33(1,3) SIGF(6)=SIGD33(2,3) C C C stockage des variables du probleme C VARF(2)=DF(1,1) VARF(3)=DF(2,2) VARF(4)=DF(3,3) VARF(5)=DF(1,2) VARF(6)=DF(1,3) VARF(7)=DF(2,3) VARF(8)=ROT(1,1) VARF(9)=ROT(1,2) VARF(10)=ROT(1,3) VARF(11)=ROT(2,1) VARF(12)=ROT(2,2) VARF(13)=ROT(2,3) VARF(14)=ROT(3,1) VARF(15)=ROT(3,2) VARF(16)=ROT(3,3) C C flag pour l'endommagement maximum VARF(17)=IFAILURE C C Deformations totales finales DO ISTRS=1,NSTRS VARF(17+ISTRS)=EPSTF(ISTRS) ENDDO C C 2000 CONTINUE C RETURN END C