VUMAT debuging

C
C Coding for Isotropic Hardening Plasticity VUMAT
C

       SUBROUTINE VUMAT(
C Read only -
     1 NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, LANNEAL,
     2 STEPTIME, TOTALTIME, DT, CMNAME, COORDMP, CHARLENGTH,
     3 PROPS, DENSITY, STRAININC, RELSPININC,
     4 TEMPOLD, STRETCHOLD, DEFGRADOLD, FIELDOLD,
     5 STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,
     6 TEMPNEW, STRETCHNEW, DEFGRADNEW, FIELDNEW,
C Write only -
     7 STRESSNEW, STATENEW, ENERINTERNNEW, ENERINELASNEW)
C
       INCLUDE 'VABA_PARAM.INC'
C
       DIMENSION PROPS(NPROPS), DENSITY(NBLOCK), COORDMP(NBLOCK),
     1 CHARLENGTH(NBLOCK), STRAININC(NBLOCK, NDIR+NSHR),
     2 RELSPININC(NBLOCK, NSHR), TEMPOLD(NBLOCK),
     3 STRETCHOLD(NBLOCK, NDIR+NSHR),DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR),
     4 FIELDOLD(NBLOCK, NFIELDV), STRESSOLD(NBLOCK, NDIR+NSHR),
     5 STATEOLD(NBLOCK, NSTATEV), ENERINTERNOLD(NBLOCK),
     6 ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),
     7 STRETCHNEW(NBLOCK, NDIR+NSHR),DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),
     8 FIELDNEW(NBLOCK, NFIELDV), STRESSNEW(NBLOCK,NDIR+NSHR),
     9 STATENEW(NBLOCK, NSTATEV), ENERINTERNNEW(NBLOCK),
     1 ENERINELASNEW(NBLOCK)

       CHARACTER*80 CMNAME

C LOCAL ARRAYS
C ----------------------------------------------------------------
C EELAS - ELASTIC STRAINS
C EPLAS - PLASTIC STRAINS
C FLOW - DIRECTION OF PLASTIC FLOW
C ----------------------------------------------------------------
C
       DIMENSION EELAS(6),EPLAS(6),FLOW(6), HARD(3)
C
       PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, SIX=6.D0,
     1 ENUMAX=.4999D0, NEWTON=10, TOLER=1.0D-6,HALF=0.5D0)

C ----------------------------------------------------------------
C UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC MISES PLASTICITY
C CANNOT BE USED FOR PLANE STRESS
C ----------------------------------------------------------------
C PROPS(1) - E
C PROPS(2) - NU
C PROPS(3..) - SYIELD AN HARDENING DATA
C CALLS UHARD FOR CURVE OF YIELD STRESS VS. PLASTIC STRAIN
C ----------------------------------------------------------------
C
C ELASTIC PROPERTIES
C
       EMOD=PROPS(1)
       ENU=MIN(PROPS(2), ENUMAX)
       EBULK3=EMOD/(ONE-TWO*ENU)
       EG2=EMOD/(ONE+ENU)
       EG=EG2/TWO
       EG3=THREE*EG
       ELAM=(EBULK3-EG2)/THREE
       NVALUE=NPROPS/2-1
C
C CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN
C
       IF ( STEPTIME .EQ. ZERO ) THEN
       DO K = 1, NBLOCK
       TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
       STRESSNEW(K,1) = STRESSOLD(K,1)
     *                  + EG2 * STRAININC(K,1) + ELAM * TRACE
       STRESSNEW(K,2) = STRESSOLD(K,2)
     *                  + EG2 * STRAININC(K,2) + ELAM * TRACE
       STRESSNEW(K,3) = STRESSOLD(K,3)
     *                  + EG2 * STRAININC(K,3) + ELAM * TRACE
       STRESSNEW(K,4)=STRESSOLD(K,4) + EG2 * STRAININC(K,4)
       IF ( NSHR .GT. 1 ) THEN
       STRESSNEW(K,5)=STRESSOLD(K,5) + EG2 * STRAININC(K,5)
       STRESSNEW(K,6)=STRESSOLD(K,6) + EG2 * STRAININC(K,6)
       END IF
       END DO
       ELSE
C
C111111111111
C
       DO  K = 1, NBLOCK
C
       TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
       STRESSNEW(K,1) = STRESSOLD(K,1)
     *                  + EG2 * STRAININC(K,1) + ELAM * TRACE
       STRESSNEW(K,2) = STRESSOLD(K,2)
     *                  + EG2 * STRAININC(K,2) + ELAM * TRACE
       STRESSNEW(K,3) = STRESSOLD(K,3)
     *                  + EG2 * STRAININC(K,3) + ELAM * TRACE
       STRESSNEW(K,4)=STRESSOLD(K,4) + EG2 * STRAININC(K,4)
       IF ( NSHR .GT. 1 ) THEN
       STRESSNEW(K,5)=STRESSOLD(K,5) + EG2 * STRAININC(K,5)
       STRESSNEW(K,6)=STRESSOLD(K,6) + EG2 * STRAININC(K,6)
       END IF
C
       SMISES=(STRESSNEW(K,1)-STRESSNEW(K,2))**2 +
     1        (STRESSNEW(K,2)-STRESSNEW(K,3))**2 +
     2        (STRESSNEW(K,3)-STRESSNEW(K,1))**2
       DO 90 K1=NDIR+1,NDIR+NSHR
         SMISES=SMISES+SIX*STRESSNEW(K,K1)**2
90     CONTINUE
       SMISES=SQRT(SMISES/TWO)
C
C RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE FORWARD
C ALSO RECOVER EQUIVALENT PLASTIC STRAIN
C
C      CALL ROTSIG(STATEOLD(K,1), DROT, EELAS, 2, NDIR, NSHR)  
C      CALL ROTSIG(STATEOLD(K,1+NDIR+NSHR), DROT, EPLAS, 2, NDIR, NSHR)

       DO K1=1, NDIR+NSHR
       EELAS(K1)=STATEOLD(K,K1)
       EPLAS(K1)=STATEOLD(K,K1+NDIR+NSHR)
       END DO
       EQPLAS=STATEOLD(K,1+2*(NDIR+NSHR))!!!
C
C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
C

       CALL VUHARD(SYIEL0, HARD, EQPLAS, EQPLASRT,TOTALTIME,DT,TEMP,
     1             DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,
     2             NUMFIELDV,NVALUE,PROPS(3))
C
C DETERMINE IF ACTIVELY YIELDING
C
       IF (SMISES.GT.(ONE+TOLER)*SYIEL0) THEN
C
C ACTIVELY YIELDING
C SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS
C CALCULATE THE FLOW DIRECTION
C
       SHYDRO=(STRESSNEW(K,1)+STRESSNEW(K,2)+STRESSNEW(K,3))/THREE
       DO K1=1,NDIR
       FLOW(K1)=(STRESSNEW(K,K1)-SHYDRO)/SMISES
       END DO
       DO K1=NDIR+1,NDIR+NSHR
       FLOW(K1)=STRESSNEW(K,K1)/SMISES
       END DO
C
C SOLVE FOR EQUIVALENT VON MISES STRESS
C AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATION
C
       SYIELD=SYIEL0
       DEQPL=ZERO
       DO KEWTON=1, NEWTON
       RHS=SMISES-EG3*DEQPL-SYIELD
       DEQPL=DEQPL+RHS/(EG3+HARD(1))
       CALL VUHARD(SYIELD,HARD,EQPLAS+DEQPL,EQPLASRT,TOTALTIME,DT,TEMP,
     1            DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NUMFIELDV,
     2            NVALUE,PROPS(3))
       IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10
       END DO
C
C WRITE WARNING MESSAGE TO THE .MSG FILE
C
       WRITE(7,2) NEWTON
    2  FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT',
     1 'CONVERGE AFTER ',I3,' ITERATIONS')
   10  CONTINUE
C
C UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND
C EQUIVALENT PLASTIC STRAIN
C
       DO K1=1,NDIR
       STRESSNEW(K,K1)=FLOW(K1)*SYIELD+SHYDRO
       EPLAS(K1)=EPLAS(K1)+THREE/TWO*FLOW(K1)*DEQPL
       EELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPL
       END DO
       DO K1=NDIR+1,NDIR+NSHR
       STRESSNEW(K,K1)=FLOW(K1)*SYIELD
       EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL
       EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL
       END DO
       EQPLAS=EQPLAS+DEQPL
C
C CALCULATE PLASTIC DISSIPATION
C
       SPD=DEQPL*(SYIEL0+SYIELD)/TWO
  ENDIF
C
C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS
C IN STATE VARIABLE ARRAY
C
       DO K1=1, NDIR+NSHR
       STATENEW(K,K1)=EELAS(K1)
       STATENEW(K,K1+NDIR+NSHR)=EPLAS(K1)
       END DO
       STATENEW(K,1+2*(NDIR+NSHR)) = STATEOLD(K,1+2*(NDIR+NSHR)) + DEQPL
C
C UPDATE THE SPECIFIC INTERNAL ENERGY -
C
       IF ( NSHR .EQ. 1 ) THEN
       STRESSPOWER = HALF * (
     * ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
     * ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
     * ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3) ) +
     * ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4)
       ELSE
       STRESSPOWER = HALF * (
     * ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
     * ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
     * ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3) ) +
     * ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4) +
     * ( STRESSOLD(K,5) + STRESSNEW(K,5) ) * STRAININC(K,5) +
     * ( STRESSOLD(K,6) + STRESSNEW(K,6) ) * STRAININC(K,6)
       END IF
       ENERINTERNNEW(K) = ENERINTERNOLD(K) + STRESSPOWER / DENSITY(K)
C
C UPDATE THE DISSIPATED INELASTIC SPECIFIC ENERGY -
C
       PLASTICWORKINC = HALF * ( SYIEL0 + SYIELD ) * DEQPL
       ENERINELASNEW(K) = ENERINELASOLD(K) + PLASTICWORKINC / DENSITY(K)
       ENDDO  
       END IF
C
       RETURN
       END

       SUBROUTINE VUHARD(SYIELD,HARD,EQPLAS,EQPLASRT,TIME,DTIME,TEMP,
     1                   DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,
     2                   NUMFIELDV, NVALUE,TABLE)
       INCLUDE 'VABA_PARAM.INC'
       CHARACTER*80 CMNAME
       DIMENSION HARD(3)
C
       DIMENSION TABLE(2, NVALUE)
C
       PARAMETER(ZERO=0.D0)

C
C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
C
       SYIELD=TABLE(1, NVALUE)
       HARD(1)=ZERO
C
C IF MORE THAN ONE ENTRY, SEARCH TABLE
C
       IF(NVALUE.GT.1) THEN
       DO K1=1, NVALUE-1
       EQPL1=TABLE(2,K1+1)
       IF(EQPLAS.LT.EQPL1) THEN
       EQPL0=TABLE(2, K1)
       IF(EQPL1.LE.EQPL0) THEN
       WRITE(7, 1)
    1  FORMAT(//, 30X, '***ERROR - PLASTIC STRAIN MUST BE ',
     1 'ENTERED IN ASCENDING ORDER')
       CALL EXIT
       ENDIF
C
C CURRENT YIELD STRESS AND HARDENING
C
       DEQPL=EQPL1-EQPL0
       SYIEL0=TABLE(1, K1)
       SYIEL1=TABLE(1, K1+1)
       DSYIEL=SYIEL1-SYIEL0
       HARD(1)=DSYIEL/DEQPL
       SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD(1)
       GOTO 10
       ENDIF
       END DO
   10  CONTINUE
       ENDIF
       RETURN
       END

Dear Dr.wei

    Thanks for your suggestion. I have find the bug in my VUMAT subroutine. Your suggestion will also be very important to me to debugging the more complicated VUMAT.

   Thanks again.