Sources Esope | Cast3M

cmicr1
C CMICR1    SOURCE    PV        17/12/08    21:16:05     9660           
C MICRO1    SOURCE    AM        00/12/13    21:40:25     4045
      SUBROUTINE CMICR1(WRK52,WRK53,WRK54,NSTRS1,NVARI,IDECAL
     &      ,LIMPR,DEFPT,EPSE1,EPSE)
      IMPLICIT INTEGER(I-N)
 
-INC PPARAM
-INC CCOPTIO
-INC DECHE
      INTEGER NSTRS1,NVARI,ISTRS,ITER,NSIG,IDECAL
      LOGICAL CONV,LCP,LIMPR
      REAL*8 SEQ,TRDFDS,YOUNG,XNU,DEUXMU,LAMB,R0,DR,DR0,RINF,RI,FTR
      REAL*8 FPC,MP,EPSE,R,TEST,TRSIG,F0,DELTL,DELTLT,EPSE1,EPSD0
      REAL*8 DEPS3,EPS3,DFDELTL,DFDEPS3,DS3DELTL,DS3DEPS3
      REAL*8 A,B,C,A1,B1,C1,A2,B2,C2,DEFPT(6)
*
*
      SEGMENT WRKK1
        REAL*8 SIGE0(NSTRS1),DFDS(NSTRS1)
      END SEGMENT
*
      SEGINI WRKK1
      YOUNG=XMAT(1)
      XNU=XMAT(2)
      FPC=XMAT(7)
      MP=XMAT(8)
      EPSD0=XMAT(5)
      EPSE=VAR0(1)
      RI=VAR0(2)
      R0=EPSD0*YOUNG
*      R0=0.4D0*FPC
      RINF=0.6d0*FPC
      LCP=IFOUR.EQ.-2
      DEUXMU=YOUNG/(1.D0+XNU)
      LAMB=XNU*DEUXMU/(1.D0-2.D0*XNU)
      TEST=1.D-6*R0
      DR=0.D0
      KERRE=0
 
 
C
C     calcul de la matrice elastique
C
      CMATE = 'ISOTROPE'
      KCAS=1
      CALL DOHMAS(XMAT,CMATE,IFOUR,NSTRS1,KCAS,DDHOOK,IRTD)
      IF ( IRTD .EQ. 0) THEN
        KERRE=123
        return
      end if
C
C      calcul de l'increment de contrainte
C
      CALL MATVE1 (DDHOOK,DEPST,NSTRS1,NSTRS1,DSIGT,1)
*      print*,'IFOUR=',IFOUR
*      print*,'DDHOOK'
*      do i=1,3
*      print*,(ddhook(i,j),j=1,3)
*      end do
*      do i=1,3
*      print*,'DEPST(',I,')=',DEPST(i)
*      print*,'DSIGT(',I,')=',DSIGT(I)
*      end do
*
*     calcul des parametres du critere de NADAI
*
C-
C- A,B,C : PARAMETRES DU CRITERE DE PLASTICITE (NADAI)
C- ILS SONT DIFFERENTS SUIVANT LE CHARGEMENT
C-
 
       CONV=.FALSE.
*       CONV=.TRUE.
      CALL ZDANUL (DEFP,NSTRS1)
      ITER=0
      EPS3=0.d0
      DELTLT=0.D0
      R=RI
      A1=0.4d0
      B1=(1.d0 - A1) / (1.D0 + A1)
      C1=(2.D0 * A1) / (1.D0 + A1)
      A2=1.16d0
*     B2=(A2 - 1.D0) / (2.D0 * A2 - 1.D0)
      B2=0.2
*     C2=A2 / (2.D0 * A2 - 1.D0)
      C2=0.8
*
*     NSIG est la dimension du tenseur des contraintes:
*     NSIG=3 sauf 2 en contraintes planes
*
      IF (LCP)THEN
        NSIG=2
      ELSE
        NSIG=3
      END IF
*
*     on recupere les contraintes anelastiques due a l'endommagement
*
      DO ISTRS=1,NSTRS1
         SIGE0(ISTRS)=VAR0(ISTRS+IDECAL)
      END DO
*
*     prediction elastique de la contrainte effective
*     contrainte + contrainte anelastique + increment de contrainte elastique
*
*     l'increment de contrainte est decoupe:
*     1) point de contact avec la surface seuil ???
*     2) decoupage du reste de facon a ce que le pas en contrainte
*        equivalente soit inferieur a disons pour le moment E*1.e-04
*
*     en deux coups les gros
*
*
*      DS=0.d0
*      do istrs=1,nstrs
*        ds=ds + dsigt(istrs)**2
*      end do
*      nstep=int(1.e+4*sqrt(ds)/youn)
      DO ISTRS=1,NSTRS1
         SIGE0(ISTRS)=SIG0(ISTRS)+SIGE0(ISTRS)+DSIGT(ISTRS)
      END DO
      if(limpr)then
       DO ISTRS=1,NSTRS1
         PRINT*,'SIG0(',ISTRS,')=',SIGE0(ISTRS)
         PRINT*,'DSIGT(',ISTRS,')=',DSIGT(ISTRS)
       END DO
      endif
 
      DO WHILE (ITER .LE. 100 .AND. .NOT. CONV)
        ITER=ITER+1
*
*       on calcule la contrainte equivalente de Mises (SEQ)
*               et la trace des contraintes (TRSIG)
*
*       on utilise la fonction modulo:
*       MOD( 1,3)= 1
*       MOD( 2,3)= 2
*       MOD( 3,3)= 0
*       MOD( 4,3)= 1
*       MOD( 5,3)= 2
*       MOD( 6,3)= 0
*
        SEQ=0.D0
        TRSIG=0.d0
* ici NSIG change en 3 le 6/11/95 (clb)
        DO ISTRS=1,3
          TRSIG=TRSIG + SIGE0(ISTRS)
          SEQ=SEQ+SIGE0(ISTRS)*(SIGE0(ISTRS)-SIGE0(MOD(ISTRS,3)+1))
        END DO
        DO ISTRS=4,NSTRS1
          SEQ=SEQ + 3.D0*SIGE0(ISTRS)**2
        END DO
        if (limpr) then
          print*,'VMIS**2=',SEQ,'TRSIG=',TRSIG
        end if
        SEQ=SQRT(SEQ)
*
*       en fonction de la valeur de la trace de sigma, on determine les
*       coef b et c a applique a la surface seuil
*
        IF (MOD(ITER-1,3).EQ.0 )then
          IF ( TRSIG .GT. (4.D0*(R+R0)*(C2-C1)/(B2-B1)))THEN
*       PRINT*,'OUI TRSIG .GT. (R+R0)'
*       PRINT*,'(R+R0)*(C2-C1)/(B2-B1)',4.D0*(R+R0)*(C2-C1)/(B2-B1)
*           B=B1
*           C=C1
            B=1.D0
                   C=2.D0
*             C=1.D0
*             B=((((0.4D0*FPC)/FTR)*C)-1)
*       PRINT*,'B',B
*       PRINT*,'C',C
          ELSE
*       PRINT*,'NON TRSIG .GT. (R+R0)'
*       PRINT*,'(R+R0)*(C2-C1)/(B2-B1)',4*(R+R0)*(C2-C1)/(B2-B1)
*           B=B2
*           C=C2
*           B=-1.D0
*           C=8.0D0
            B=0.25D0
                   C=3.0D0
*       PRINT*,'B',B
*       PRINT*,'C',C
          ENDIF
        ENDIF
*       B=1.D0
*       C=2.D0
*
*       fonction seuil
*
          F0=(SEQ + B*TRSIG)/C - (R0 + R)
        if (limpr) then
          PRINT*,'B=',B,'C=',C,'F0=',F0
          PRINT*,'SIGMAEQ=',(SEQ + B*TRSIG)/C
          PRINT*,'R0',R0
          PRINT*,'(1.D0+B)/C',(1.D0+B)/C
          PRINT*,'B', B
          PRINT*,'C', C
          PRINT*,'TRSIG',TRSIG,'SEQ',SEQ
          PRINT*,'F0',F0
        end if
*
*       le critere de convergence est verifie si la valeur absolue
*       de la fonction seuil est plus petite que le critere ou si
*       la fonction seuil est negative pour la premiere iteration
*
        CONV=(ABS(F0).LE.TEST .AND.
     1       (.NOT. LCP .OR. ABS(SIGE0(3)) .LE. TEST ) )
     2        .OR. (F0 .LE. 0.D0 .AND. ITER .EQ. 1)
*
*       dans le cas des contraintes planes on rajoute aussi
*       que la valeur absolue de sigma33 est inferieure au critere
*
*       PRINT*,'CONV',CONV
*       PRINT*,'F(0)',F0
        IF (.NOT. CONV) THEN
*        PRINT*,'On ecoule'
*
*       calcul de dF/dSigma
*
          DO ISTRS=1,3
            DFDS(ISTRS)=(2.D0*SIGE0(ISTRS)-SIGE0(MOD(ISTRS,3)+1)
     &              -SIGE0(MOD(ISTRS+1,3)+1) )/(2.D0*C*SEQ) +B/C
          END DO
*       PRINT*,'DFDS(11)',DFDS(1)
*       PRINT*,'DFDS(22)',DFDS(2)
*       PRINT*,'DFDS(33)',DFDS(3)
          DO ISTRS=4,NSTRS1
            DFDS(ISTRS)=3.D0*SIGE0(ISTRS)/(2.d0*C*SEQ)
          END DO
          TRDFDS=0.D0
          DO ISTRS=1,NSIG
            TRDFDS=TRDFDS+DFDS(ISTRS)
          END DO
          if (limpr) then
*           on verifie que les valeur ppales sont bonnes
            do istrs=1,nstrs
              print*,'dfds(',istrs,')=',dfds(istrs)
            end do
            print*,'trdfds=',trdfds
          end if
*
*        calcul de dsigma/deltl=C * depsp
*        attention !! dans le cas des contraintes planes:epsp33 n'est
*        pas une inconnue directe mais c'est depse33
*
*
*         R=MP*p  -----> DR=MP * Dp  =MP*DELTL
*         R=Rinf(1-exp(-bR))  ---> dR=b(Rinf-R)deltl
*
*          DR = (RINF - R)*MP
*          DR=MP
*          R=MP*(p)**0.5 ---> DR=0.5 * MP *P**(-05)dp
*         Si On prend un écrouissage linéaire
*         R=MP*(p)**1.D0 ---> DR= MP dp
*
 
*          IF ((EPSE+DELTLT).GE. 1.D-20 )THEN
*             DR=0.5*MP*(EPSE+DELTLT)**(-0.5)
              DR=MP
*       PRINT*,'OUI (EPSE+DELTLT).GE. 1.D-20 )',DR
*          ELSE
*             DR=0.D0
*          END IF
          DFDELTL=-DR
          DO ISTRS=1,NSIG
            DFDELTL=DFDELTL-DFDS(ISTRS)*(DFDS(ISTRS)*DEUXMU+TRDFDS*LAMB)
          END DO
          DO ISTRS=4,NSTRS1
            DFDELTL=DFDELTL - 2.d0*DEUXMU * DFDS(ISTRS)**2
          END DO
          IF (LCP) THEN
            DFDELTL=DFDELTL-DFDS(3)*TRDFDS*LAMB
            DFDEPS3=LAMB*TRDFDS + (LAMB+DEUXMU)*DFDS(3)
            DS3DELTL=-LAMB*TRDFDS
            DS3DEPS3=LAMB + DEUXMU
            DELTL=(SIGE0(3)*DFDEPS3 - F0*DS3DEPS3)/
     &            (DS3DEPS3*DFDELTL - DS3DELTL*DFDEPS3)
            DEPS3=-(SIGE0(3)-DS3DELTL*DELTL)/DS3DEPS3
            EPS3=EPS3+DEPS3
            DELTLT=DELTLT+DELTL
            DO ISTRS=1,NSIG
               SIGE0(ISTRS)=SIGE0(ISTRS)+LAMB*DEPS3
            END DO
            SIGE0(3)=SIGE0(3)+(LAMB+DEUXMU)*DEPS3-LAMB*DELTL*TRDFDS
          ELSE
            DELTL=-F0/DFDELTL
            DELTLT=DELTLT+DELTL
          END IF
          DO ISTRS=1,NSIG
            SIGE0(ISTRS)=SIGE0(ISTRS)
     &               -DELTL*(DEUXMU*DFDS(ISTRS)+LAMB*TRDFDS)
          END DO
          DO ISTRS=4,NSTRS1
            SIGE0(ISTRS)=SIGE0(ISTRS)-DELTL*DEUXMU*DFDS(ISTRS)
          END DO
        ELSE
                DELTL=0.D0
        END IF
        R=R + DR*DELTL
        DO ISTRS=1,3
          DEFP(ISTRS)=DEFP(ISTRS)+DELTL*DFDS(ISTRS)
        END DO
        DO ISTRS=4,NSTRS1
          DEFP(ISTRS)=DEFP(ISTRS)+2.D0*DELTL*DFDS(ISTRS)
        END DO
         if (limpr) then
           print*,'ITER=',ITER,'F0=',F0,'DELTL=', DELTL
           print*,'DELTLT=',DELTLT,'DR=',DR,'R',R
           do istrs=1,nstrs
             print*,'SIGMA(',ISTRS,')=', SIGE0(ISTRS)
             print*,'DEFP(',ISTRS,')=', DEFP(ISTRS)
           end do
         end if
*       PRINT*,'EPSE+DELTLT',EPSE+DELTLT
      END DO
*      EN CAS DE NON CONVERGENCE: ON MET KERRE a 1
      IF (.NOT. CONV) THEN
              KERRE=1
      END IF
*      PRINT*,'CONVERGENCE OBTENUE APRES ',ITER, 'ITERATIONS'
*
*  on remplit les variables finales
 
**
* Je rajoute ça pour les déformations plastiques
*
        DO ISTRS=1,NSTRS1
        DEFPT(ISTRS)=DEFP(ISTRS)
        ENDDO
      VARF(1)=EPSE+DELTLT
      VARF(2)=R
        EPSE1=EPSE*((1.D0+B)/C)
*       EPSE=EPSE1*(C/(1.D0+B))
*
*  ON REND LES CONTRAINTES
*
      DO ISTRS=1,NSTRS1
        SIGF(ISTRS)=SIGE0(ISTRS)
*        PRINT*,'SIGF(',ISTRS,')=',SIGF(ISTRS)
      END DO
 
* ON REND L'INCREMENT DE DEFORMATIONS PLASTIQUES
*
 
      SEGSUP WRKK1
      RETURN
      END
*