Sources Esope | Cast3M

increp
C INCREP    SOURCE    CB215821  17/11/30    21:16:27     9639           
C INCREP    SOURCE    BV      03/11/04    21:29:05     1995
       SUBROUTINE INCREP(SIG,VAR,EPSVPT,VARPT,XMAT,NSTRS0,
     &                                     MFR,NVARI,NCOMAT)
       IMPLICIT REAL*8(A-H,O-Z)
       DIMENSION SIG(*),VAR(*),EPSVPT(*),VARPT(*),XMAT(*)
       DIMENSION XX(6),EPS0(6),XP(3,3)
       DIMENSION ST(3,3),D(3),ZMU(3,3)
       DIMENSION TEMP(3),EPSD(3),VPT(6)
       PARAMETER (AMAX = 1.0D20 , AMIN = 1.D-10)
-INC CCREEL
       DETIER = 2.0D0/3.0D0
       ROOT = SQRT(DETIER)
C-------------------------------------------------------------------|
C*******  EVALUATION OF J2 SIGMA                                    |
C-------------------------------------------------------------------I
       TRACE =(SIG(1)+SIG(2)+SIG(3))/3.0D0
       DO 70 I=1,NSTRS0
         A = 0.0D0
         IF (I.LE.3) A=1.0D0
         XX(I) = SIG(I)-A*TRACE
  70   CONTINUE
       AJ2 = PROCON (XX,XX,NSTRS0)
       AJ2 = SQRT(1.5D0*AJ2)
C------------------------------------------------------
C*******  EVALUATION OF J0 CONTRAINTE PRINCIPALE MAX
C------------------------------------------------------
       CALL ZERO(ST,3,3)
       DO 50 I=1,3
         ST(I,I) = SIG(I)
50     CONTINUE
       ST(1,2) = SIG(4)
       IF(NSTRS0.NE.6) GOTO 54
       ST(1,3) = SIG(5)
       ST(2,3) = SIG(6)
54     CONTINUE
       ST(2,1) = ST(1,2)
       ST(3,1) = ST(1,3)
       ST(3,2) = ST(2,3)
       CALL ZERO(XP,3,3)
*
       CALL JACOB3(ST,3,D,XP)
C
       AJ0 = D(1)
*       AJ0 =  MAX(D(1),D(2),D(3))
C--------------------------------------------------------
C*******  CALCUL DE LA CONTRAINTE EQUIVALENTE AU SENS
C*******  DE L'ENERGIE DE CISAILLEMENT (VON MISES)
C--------------------------------------------------------
       DD = VAR(3)
       DC = VAR(4)
       IF (DD.GT.0.98D0) THEN
          DD =  0.98D0
       ENDIF
       IF (DC.GT.0.98D0) THEN
          DC =  0.98D0
       ENDIF
       ZMU(1,1) = 2.*(1.+(XMAT(2)*(1.-DC)))
       ZMU(1,2) = (DC / 2.0D0)-(1.0D0+XMAT(2)*(1.0D0-DC))
       ZMU(1,3) = ZMU(1,2)
       ZMU(2,2) = 2.0D0*(1.0D0+XMAT(2))*(1.0D0-DC)
       ZMU(2,3) = -1.0D0*(1.0D0+XMAT(2))*(1.0D0-DC)
       ZMU(3,3) = ZMU(2,2)
       ZMU(2,1) = ZMU(1,2)
       ZMU(3,1) = ZMU(1,3)
       ZMU(3,2) = ZMU(2,3)
       DO 18 J=1,3
         DO 17 I=1,3
            ZMU(I,J) =  ZMU(I,J)/(2.0D0*(1.0D0+XMAT(2)))
  17     CONTINUE
  18   CONTINUE
       TEMP(1) = 0.0D0
       TEMP(2) = 0.0D0
       TEMP(3) = 0.0D0
       TEMP(1) = ZMU(1,1)*D(1)+ZMU(1,2)*D(2)+ZMU(1,3)*D(3)
       TEMP(2) = ZMU(2,1)*D(1)+ZMU(2,2)*D(2)+ZMU(2,3)*D(3)
       TEMP(3) = ZMU(3,1)*D(1)+ZMU(3,2)*D(2)+ZMU(3,3)*D(3)
       SEQE = TEMP(1)*D(1)+TEMP(2)*D(2)+TEMP(3)*D(3)
       SEQE = SEQE/((1.0D0-DD)*(1.0D0-DC))
       SEQE = SQRT(SEQE)
       SEQE=max(xpetit,seqe)
C--------------------------------------------------------------------|
C*******  CALCUL OF SIGV/K                                           |
C -------------------------------------------------------------------I
       RR = VAR(1)
       SK = (SEQE - RR-XMAT(5))/( XMAT(8) )
C--------------------------------------------------------------------|
C*******  CALCULATION OF EFFECTIVE INELASTIC STRAIN INCREMENT (P)    |
C -------------------------------------------------------------------I
       IF (SK.GT.0.0D0) THEN
         RPT= SK**XMAT(9)
       ELSE
         RPT= 0.0D0
       ENDIF
C-----------------------------------------------------------------|
C*******  CALCULATION OF DR                                       |
C-----------------------------------------------------------------I
       DR= XMAT(6)*(XMAT(7)-VAR(1))*RPT
       VARPT(1)=DR
C---------------------------------------------------------------------|
C*******  CALCULATION OF INELASTIC STRAIN INCREMENTS  (EPSVPT)        |
C---------------------------------------------------------------------I
C  Calcul des deformations principales
       EPSD(1) = ZMU(1,1)*D(1)+ZMU(1,2)*D(2)+ZMU(1,3)*D(3)
       EPSD(1) = EPSD(1)*RPT/(SEQE*(1.0D0-DD)*(1.0D0-DC))
       EPSD(2) = ZMU(2,1)*D(1)+ZMU(2,2)*D(2)+ZMU(2,3)*D(3)
       EPSD(2) = EPSD(2)*RPT/(SEQE*(1.0D0-DD)*(1.0D0-DC))
       EPSD(3) = ZMU(3,1)*D(1)+ZMU(3,2)*D(2)+ZMU(3,3)*D(3)
       EPSD(3) = EPSD(3)*RPT/(SEQE*(1.0D0-DD)*(1.0D0-DC))
C
       IF (SEQE.EQ.0.0D0) THEN
         PPT = 0.0D0
       ELSE
         PPT = DETIER*(EPSD(1)*EPSD(1)+EPSD(2)*EPSD(2)+EPSD(3)*EPSD(3))
         PPT = SQRT(PPT)
       ENDIF
C
       VARPT(2) = PPT
       DO 11 I=1,NSTRS0
         EPS0(I)=0.0D0
  11   CONTINUE
C
       EPS0(1)=XP(1,1)*XP(1,1)*EPSD(1)
       EPS0(1)=EPS0(1)+XP(1,2)*XP(1,2)*EPSD(2)
       EPS0(1)=EPS0(1)+XP(1,3)*XP(1,3)*EPSD(3)
       EPS0(2)=XP(2,1)*XP(2,1)*EPSD(1)
       EPS0(2)=EPS0(2)+XP(2,2)*XP(2,2)*EPSD(2)
       EPS0(2)=EPS0(2)+XP(2,3)*XP(2,3)*EPSD(3)
       EPS0(3)=XP(3,1)*XP(3,1)*EPSD(1)
       EPS0(3)=EPS0(3)+XP(3,2)*XP(3,2)*EPSD(2)
       EPS0(3)=EPS0(3)+XP(3,3)*XP(3,3)*EPSD(3)
       EPS0(4)=XP(1,1)*XP(2,1)*EPSD(1)
       EPS0(4)=EPS0(4)+XP(1,2)*XP(2,2)*EPSD(2)
       EPS0(4)=EPS0(4)+XP(1,3)*XP(2,3)*EPSD(3)
       IF(NSTRS0.NE.6) GOTO 56
         EPS0(5)=XP(1,1)*XP(3,1)*EPSD(1)
         EPS0(5)=EPS0(5)+XP(1,2)*XP(3,2)*EPSD(2)
         EPS0(5)=EPS0(5)+XP(1,3)*XP(3,3)*EPSD(3)
         EPS0(6)=XP(2,1)*XP(3,1)*EPSD(1)
         EPS0(6)=EPS0(6)+XP(2,2)*XP(3,2)*EPSD(2)
         EPS0(6)=EPS0(6)+XP(2,3)*XP(3,3)*EPSD(3)
  56   CONTINUE
C
       DO 47 I=1,NSTRS0,1
         IF (RPT.EQ.0.0) THEN
            EPSVPT(I)=0.0D0
         ELSE
            EPSVPT(I)=EPS0(I)
         ENDIF
  47   CONTINUE
C--------------------------------------------------------
C*******  CALCUL de D ductile AVEC SEUIL D ENDOMMAGEMENT
C--------------------------------------------------------
       YD = (D(1)*D(1)/(1.0D0-DC))+(D(2)*D(2))+(D(3)*D(3))
       YD = YD-2.0D0*XMAT(2)*(D(1)*D(2))+(D(2)*D(3)+D(1)*D(3))
       YD = YD/(2.0D0*XMAT(1)*(1.0D0-DD)*(1.0D0-DD))
C
       P = VAR(2)
       SE1 = XMAT(12)
       IF ((P.GT.SE1).AND.(DC.LT.0.98)) THEN
          DDD =  ((YD/XMAT(10))**XMAT(11))*RPT
       ELSE
          DDD = 0.0D0
       ENDIF
       VARPT(3)= DDD
C       write(6,*)'DDD',DDD
C-----------------------------------------------------
C*******  CALCUL de D creep AVEC SEUIL D ENDOMMAGEMENT
C-----------------------------------------------------
       YC = (D(1)*D(1)/((1.0D0-DC)**2))
       YC = YC/(2.0D0*XMAT(1)*(1.0D0-DD))
C
       SE2 = XMAT(15)
       IF ((P.GT.SE2).AND.(DC.LT.0.98)) THEN
          DDC =  (YC/XMAT(13))**XMAT(14)
       ELSE
          DDC = 0.0D0
       ENDIF
       VARPT(4)= DDC
C       write(6,*)'DDC',DDC
C
       RETURN
       END