Numérotation des lignes :

C DESMOR    SOURCE    PV        17/12/08    21:17:20     9660           C DESMORAT (V2)      SUBROUTINE DESMOR(WRK52,WRK53,WRK54,NVARI,Iecou)CC     variables en entreeCC     WRK52, WRK53, WRK54 pointeur actif a l'entreeC     NVARI      nombre de variables internes (doit etre egal a ?)CC     ISTEP      flag utilise pour separer les etapes dans un calcul non localC                ISTEP=0 -----> calcul localC                ISTEP=1 -----> calcul non local etape 1 on calcule les seuilsC                ISTEP=2 -----> calcul non local etape 2 on continue le calculC                               a partir des seuils moyennesCC     variables en sortieCC     VARINF      variables internes finalesCC     SIGMAF      contraintes finalesCC      IMPLICIT REAL*8(A-H,O-Z)C-INC CCOPTIO-INC DECHECC      DIMENSION XD0(3,3),DF(3,3),POSEPSI2(3,3),EPS33(3,3)      DIMENSION SIGTILDE33(3,3),SIGDEV33(3,3),SIGD33(3,3),SIGTILDE(6)      DIMENSION UNMOINSDVP(3),HM1(3,3),DFD(3,3),UNMOINSD(3,3),DFVP(3)      DIMENSION ROT(3,3),ROTT(3,3),VECP33T(3,3),TRAV(3,3),VECP33(3,3)      DIMENSION EPSIPOS(3),HM1DIAG(3),EPS(2,2),DELTAF22(2,2),EPSIPP(3)      DIMENSION VPDF(2),VECPDF(2,2),EPSTF(6)C      PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)C      SEGMENT IECOU*      COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,       INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,C      INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM,  NYKK,     1  icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,C    .  NYALF1,NYBET1,NYR,  NYA,   NYRHO,NSIGY,  NNKX,  NYKX,  IND,     2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,C    .  NSOM,  NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT,  NCXMAT,     3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA,C    .  LTRAC, MFR,  IELE,  NHRM,   NBNN, NBELEM,ICARA,     4  icow32,icow33,NSTRS1,MFR1,icow36,icow37,icow38,C    .  LW2,   NDEF,  NSTRSS,MFR1,  NBGMAT,NELMAT,MSOUPA,     5  icow39,icow40,icow41,icow42,icow43,icow44C    .  NUMAT1,LENDO, NBBB,  NNVARI,KERR1, MELEME        INTEGER icow45,icow46,icow47,icow48,icow49,icow50,     .          icow51,icow52,icow53,icow54,icow55,icow56     .          icow57,icow58      ENDSEGMENTCCC Recuperation des variables initiales dans les tableauxC      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)CC Recuperation de l'endommagement initialC      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)CC Recuperation de la base d'endommagementC      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)CC Recuperation de l'etat d'endommagement critiqueC      IFAILURE=nint(VAR0(17))CC  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)=GEGECCC Deformation totaleC      DO ISTRS=1,NSTRS         EPSTF(ISTRS)=VAR0(17+ISTRS)+DEPST(ISTRS)      ENDDOC - Trace de la deformation      TRAEPS=EPSTF(1)+EPSTF(2)+EPSTF(3)CC  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) THENc Over nonlocal, coef 2C         EPSTIL=DEUX*VAR0(1)-EPSTILNLOC         EPSTIL=VAR0(1)      ENDIFCC**********************************************************C Debut de l'intégration numerique du comportementC**********************************************************CC     calcul des seuilsCC  partie positive de la matrice des deformations      CALL PRODT2(POSEPSI2,EPSIPOS,VECP33,3)CC calcul de la trace de D au pas n      TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)CC  calcul du seuil      XKAPPA=A*TAN(TRAD0/(A*AA)+ATAN(XK0/A))      FCRIT=EPSTIL+Xk*TRAEPS-XKAPPACC  TEST sur le seuil      IF(FCRIT.LT.XZERO) THENC       on est elastiqueC       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      ELSECC       debut de la correction d'endom.CC       Calcul de Kappa^(-1)=Tr Dn+1C        TRADF=A*AA*(ATAN((EPSTIL+Xk*TRAEPS)/A)-ATAN(XK0/A))CC       test sur l'eventuelle presence d'une direction d'endommagement fixeeC        IF(IFAILURE.EQ.0)THENC         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          ENDDOC         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=1C           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        ENDIFC        IF(IFAILURE.EQ.1)THENC       une direction d'endom est fixee a DcCC       On calcule le multiplicateur d'endommagement dans la bonne baseC       On met les deformations positives dans le repère principalC                 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          ENDIFCC         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)CC         On verifie la deuxieme direction d'endommagement          CALL JACOD3(DF,3,DFVP)          IF(DFVP(2).GE.Dc)THEN            IFAILURE=2            PRINT*, 'IFAILURE=', IFAILUREC           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        ENDIFC        IF(IFAILURE.EQ.2) THENC         deux directions d'endom sont fixee a Dc          CALL TRSPOD(ROT,3,3,ROTT)C          DF33=TRADF-DEUX*DcC         print*, 'DF33=',DF33          IF(DF33.LT.XZERO) THEN            DF33=XPETIT          ENDIFC         On corrige l'endommagement          CALL ZERO(DF,3,3)          DF(1,1)=Dc          DF(2,2)=Dc          DF(3,3)=DF33C          CALL MATMAT(DF,ROTT,3,3,3,TRAV)          CALL MATMAT(ROT,TRAV,3,3,3,DF)CC         On verifie la troisieme direction d'endommagement          IF(DF33.GE.Dc)THEN            IFAILURE=3C           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        ENDIFC        IF(IFAILURE.EQ.3)THENC         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)=DcC          CALL MATMAT(DF,ROTT,3,3,3,TRAV)          CALL MATMAT(ROT,TRAV,3,3,3,DF)        ENDIFCC   Fin de la correction d'endommagementCC        DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0C      ENDIFCC  Fin des differents cas (elasticite ou endommagement)CC     Calcul de la contrainte finaleC      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)CC    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      ENDIFC      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)CCC stockage des variables du problemeC      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)CC flag pour l'endommagement maximum      VARF(17)=IFAILURECC  Deformations totales finales      DO ISTRS=1,NSTRS       VARF(17+ISTRS)=EPSTF(ISTRS)      ENDDOC2000  CONTINUECC      RETURN      ENDC

© Cast3M 2003 - Tous droits réservés.
Mentions légales