Sources Esope | Cast3M

cmicri
C CMICRI    SOURCE    OF166741  25/11/04    21:15:28     12349          
 
      SUBROUTINE CMICRI (WRK52,WRK53,WRK54,NVARI,iecou)
*
*   modele d'endommagement microplan isotrope couple a la plasticite
*   C. La Borderie + S. Fichant Oct. 95
*   routines utilisees:
*   micro1: plasticite nadai
*            IDECAL=3 DANS LE CAS ISO IDECAL=8 DANS LE CAS ANISO
*   jacob3: diagonalisation:
*            attention jacob3 modifie la matrice a diagonaliser!!
*   prodt et prodt2
*   attention prodt2 ne fonctionne qu'avec la matrice des V. P. !!
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
      INTEGER NVARI
      INTEGER ISTRS,I,J
      REAL*8 YOUNG,XNU,EPSD0,BT,LAMB,DEUXMU,ALFA,MP,BT1
      REAL*8 DOM,SIGAN(6),TRSIG,TRSIG33,DEF33(3,3),EPSIPP(3),VECP(3,3)
      REAL*8 D1,DEFPT(6),DEFTOT(6),EPSITOP(3),EPSE1,EPSE
      REAL*8 SIGPP(3),SIGPM(3),SIG33(3,3),SIG33P(3,3),SIG33M(3,3)
      REAL*8 LAMBDAP,LAMBDAM
      INTEGER IDECAL
 
-INC PPARAM
-INC CCOPTIO
-INC DECHE
-INC TECOU
 
      SEGMENT WRKK1
        REAL*8 DEFELA(NSTRS1)
      ENDSEGMENT
*
*   on recupere les variable materielles
*
      YOUNG=XMAT(1)
      XNU=XMAT(2)
      EPSD0=XMAT(5)
      BT=XMAT(6)
      MP=XMAT(8)
      BT1=(1.D0-(YOUNG/(YOUNG+MP)))
      ALFA=XMAT(9)
      DEUXMU=YOUNG/(1.D0+XNU)
      LAMB=XNU*DEUXMU/(1.D0-2.D0*XNU)
*
*     recuperation des variables internes d'endommagement
*
      DOM=VAR0(3)
*
*     on ecoule plastiquement sur la contrainte effective
*
      IDECAL=3
      NSTRS1 =iecou.nstrss
 
*               PRINT*,'ON ECOULE'
      CALL CMICR1(WRK52,WRK53,WRK54,NSTRS1,NVARI,IDECAL,
     &             .FALSE.,DEFPT,EPSE1,EPSE)
      IF (KERRE .NE. 0) THEN
           print*,'on n''a pas converge dans micro1'
       CALL CMICR1(WRK52,WRK53,WRK54,NSTRS1,NVARI,IDECAL,
     &             .FALSE.,DEFPT,EPSE1,EPSE)
         RETURN
      ENDIF
*
*     on ecoule en endommagement sur les deformations elastiques
 
      SEGINI WRKK1
*
*     calcul des deformations elastiques
*
*      print*,'-----sigf------'
*      print*,sigf(1),sigf(2),sigf(3)
      TRSIG = SIGF(1)+SIGF(2)+SIGF(3)
      DO ISTRS=1,3
        DEFELA(ISTRS)=( (1.D0+XNU)*SIGF(ISTRS)-XNU*TRSIG)/YOUNG
        DEFTOT(ISTRS)=DEFELA(ISTRS)+DEFPT(ISTRS)
      END DO
*      print*,'-----defela------'
*      print*,defela(1),defela(2),defela(3)
*      print*,'-----defpt------'
*      print*,defpt(1),defpt(2),defpt(3)
         DO ISTRS=4,NSTRS1
           DEFELA(ISTRS)= (1.D0+XNU)*SIGF(ISTRS)/YOUNG
           DEFTOT(ISTRS)=DEFELA(ISTRS)+DEFPT(ISTRS)
         END DO
*
*     on met les deformations sous forme de matrice 3x3
*     pour calculer les valeurs propres
*
         CALL ENDOCA (DEFELA,DEF33,1)
*         print*,'deformations elastiques dans rpg(3x3)'
*         print*,def33
*      print*,'prodt defrdpe'
         CALL JACOB3 (DEF33,3,EPSIPP,VECP)
*         print*,'deformations principales'
*         print*,(epsipp(i),i=1,3)
*
*      PAREIL POUR DEFTOTAL
*
         CALL ENDOCA (DEFTOT,DEF33,1)
         CALL JACOB3 (DEF33,3,EPSITOP,VECP)
*         print*,'deformations TOTAL principales'
*         print*,(epsitop(i),i=1,3)
*
*     on calcule l'endommagement resultant
*
*         print*,'BT=',BT,'EPSD0=',EPSD0,'EPSIPP(1)',EPSIPP(1)
         IF ( (EPSIPP(1)) .GT. (EPSD0) ) THEN
*         IF ( EPSITOP(1) .GT. (EPSD0) ) THEN
*       PRINT*,'OUI On calcul l endommagement'
*       PRINT*,'EPSIPP(1)/BT1',EPSIPP(1)/BT1
*       PRINT*,'EPSITOP(1)',EPSITOP(1)
*        PRINT*,'BT1',BT1
      D1=1.D0-(((EPSD0)/EPSIPP(1))*EXP(BT*(EPSD0- EPSIPP(1))))
*          D1=1.D0-EPSD0/EPSIPP(1)*EXP(BT*(EPSD0 - (EPSIPP(1)+EPSE1)))
*           D1=1.D0-EXP(-BT*(EPSE))
*               D1=0.D0
         ELSE
           D1=0.D0
         END IF
*         PRINT*,'D1=',D1
*         PRINT*,'EPSIPP(1)=',EPSIPP(1)
*         PRINT*,'EPSE1=',EPSE1
*         PRINT*,'EPSE',EPSE
*
*     et on en l'endommagement final
*
         IF(d1.gt.dom)then
            dom=d1
         endif
*         print*,'DOM=',DOM
*
*     on separe les contraintes effectives en + et - dans rpsigma
*
         CALL ENDOCA (SIGF,SIG33,1)
         CALL JACOB3 (SIG33,3,SIGPP,VECP)
*         print*,'contraintes ppales'
*         print*,sigpp
         DO I=1,3
           IF (SIGPP(I) .LT. 0.D0)THEN
              SIGPM(I)=SIGPP(I)
              SIGPP(I)=0.D0
           ELSE
              SIGPM(I)=0.D0
           END IF
         END DO
         CALL PRODT2(SIG33P,SIGPP,VECP,3)
         CALL PRODT2(SIG33M,SIGPM,VECP,3)
*         print*,'contraintes dans rpg'
*         print*,sig33p
*         print*,sig33m
*
       LAMBDAP=1.D0-DOM
       LAMBDAM=1.D0-DOM**ALFA
       SIG33(1,1)=LAMBDAP*SIG33P(1,1)+LAMBDAM*SIG33M(1,1)
       SIG33(1,2)=LAMBDAP*SIG33P(1,2)+LAMBDAM*SIG33M(1,2)
       SIG33(1,3)=LAMBDAP*SIG33P(1,3)+LAMBDAM*SIG33M(1,3)
       SIG33(2,1)=SIG33(1,2)
       SIG33(2,2)=LAMBDAP*SIG33P(2,2)+LAMBDAM*SIG33M(2,2)
       SIG33(2,3)=LAMBDAP*SIG33P(2,3)+LAMBDAM*SIG33M(2,3)
       SIG33(3,1)=SIG33(1,3)
       SIG33(3,2)=SIG33(2,3)
       SIG33(3,3)=LAMBDAP*SIG33P(3,3)+LAMBDAM*SIG33M(3,3)
*
*         Modif Mohammed calcul des def total
         TRSIG33=SIG33(1,1)+SIG33(2,2)+sig33(3,3)
         DEFTOT(1)=( (1.D0+XNU)*SIG33(1,1)-XNU*TRSIG33)/YOUNG
*
*     on rend les contraintes et les variables internes finales
*
      SIGAN(1)=SIGF(1)-SIG33(1,1)
      SIGF(1)=SIG33(1,1)
      VARF(3)=MAX(DOM,0.d0)
*      VARF(3)=DEFELA(1)
      SIGAN(2)=SIGF(2)-SIG33(2,2)
      SIGF(2)=SIG33(2,2)
      SIGAN(3)=SIGF(3)-SIG33(3,3)
      SIGF(3)=SIG33(3,3)
      SIGAN(4)=SIGF(4)-SIG33(1,2)
      SIGF(4)=SIG33(1,2)
      IF(IFOUR.GE.1.OR.IFOUR.LE.-3) THEN
        SIGAN(5)=SIGF(5)-SIG33(1,3)
        SIGF(5)=SIG33(1,3)
        SIGAN(6)=SIGF(6)-SIG33(2,3)
        SIGF(6)=SIG33(2,3)
      ELSE
        SIGAN(5)=0.D0
        SIGAN(6)=0.D0
      END IF
      DO ISTRS=1,6
         VARF(ISTRS+3)=SIGAN(ISTRS)
      END DO
      SEGSUP WRKK1
 
      RETURN
      END