UMAT for Kinematic hardening

I just followed the code that ABAQUS presented (Copied below), but I'm getting very strange results comparing with ABAQUS GUI, what's wrong with that? 

Please help if you have experience 

 

 

      SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD, RPL,

     1 DDSDDT, DRPLDE, DRPLDT, STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP,

     2 PREDEF, DPRED, CMNAME, NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS,

     3 COORDS, DROT, PNEWDT, CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER,

     4 KSPT, KSTEP, KINC)

C

      INCLUDE 'ABA_PARAM.INC'

C

      CHARACTER*8 CMNAME

C

      DIMENSION STRESS(NTENS), STATEV(NSTATV), DDSDDE(NTENS, NTENS),

     1 DDSDDT(NTENS), DRPLDE(NTENS), STRAN(NTENS), DSTRAN(NTENS),

     2 PREDEF(1), DPRED(1), PROPS(NPROPS), COORDS(3), DROT(3, 3),

     3 DFGRD0(3, 3), DFGRD1(3, 3)

C

C     LOCAL ARRAYS

C -----------------------------------------------------------------

C     EELAS - ELASTIC STRAINS

C     EPLAS - PLASTIC STRAINS 

C     ALPHA - SHIFT TENSOR

C     FLOW  - DIRECTION OF PLASTIC FLOW

C     OLDS  - STRESS AT START OF INCREMENT

C     OLDPL - PLASTIC STRAIN AT START OF INCREAMENT

C -----------------------------------------------------------------

C      

      DIMENSION EELAS(6),EPLAS(6), ALPHA(6), FLOW(6), OLDS(6), OLDPL(6)

C    

      PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, SIX=6.D0,

     1           ENUMAX=0.4999D0, TOLER=1.0D-6)

C

C ----------------------------------------------------------------

C    UMAT FOR ISOTRPIC ELASTICITY AND ISOTROPIC MISES PLASTICITY 

C    WITH KINEMATIC HARDENING - CANNOT BE USED FOR PLANE STRESS

C ----------------------------------------------------------------

C    PROPS(1) - E

C    PROPS(2) - NU

C    PROPS(3)- SYIELD 

C    PROPS(3)- HARD

C ----------------------------------------------------------------

C

C    ELASTICS PROPORTIES

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

C

C     ELASTIC STIFFNESS

C

      DO K1=1, NDT

          DO K2=1, NDI

              DDSDDE(K2,K1)=ELAM

          END DO

          DDSDDE(K1,K1)=EG2+ELAM

      END DO

      DO K1=NDI+1, NTENS

          DDSDDE(K1,K1)=EG

      END DO

C     

C     RECOVER ELASTIC STRAIN, PLASTIC STRAIN AND SHIFT TENSOR AND ROTATE

C     NOTE: USE CODE 1 FOR (TENSOR) STRESS , CODE 2 FOR (ENGINEERING) STRAIN

C

      CALL ROTSIG(STATEV(        1), DROT, EELAS, 2, NDI, NSHR)

      CALL ROTSIG(STATEV(  NTENS+1), DROT, EPLAS, 2, NDI, NSHR)

      CALL ROTSIG(STATEV(2*NTENS+1), DROT, ALPHA, 1, NDI, NSHR)

C

C     SAVE STRESS AND PLASTIC STRAINS AND 

C     CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN

C

      DO K1=1, NTENS

          OLDS(K1)=STRESS(K1)

          OLDPL(K1)=EPLAS(K1)

          EELAS(K1)=EELAS(K1)+DSTRAN(K1)

          DO K2=1, NTENS

              STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)

          END DO

      END DO

C

C     CALCULATE EQUIVALENT VON MISES STRESS

C

       SMISES= (STRESS(1)-ALPHA(1)-STRESS(2)+ALPHA(2))**2

     1       + (STRESS(2)-ALPHA(2)-STRESS(3)+ALPHA(3))**2

     2       + (STRESS(3)-ALPHA(3)-STRESS(1)+ALPHA(1))**2

      DO K1=NDI+1, NTENS

          SMISES=SMISES+SIX*(STRESS(K1)-ALPHA(K1))**2

      END DO

      SMISES=SQRT(SMISES/TWO)

C  

C     GET YIELD STRESS AND HARDENING MODULUS

C

      SYIELD=PROPS(3)

      HARD=PROPS(4)

C     

C     DETERMINE IF ACTIVELY YIELDING 

C

      IF(SMISES.GT.(ONE+TOLAR)*SYIELD) THEN

C

C     ACTIVELY YIELDING

C     SEPARATED THE HYDORSTATIC FROM THE DEVATORIC STRESS

C     CALCULATE THE FLOW DIRECTION 

C

      SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE

      DO K1=1, NDI

          FLOW(K1)=(STRESS(K1)-ALPHA(K1)-SHYDRO)/SMISES

      END DO

      DO K1=NDI+1, NTENS

          FLOW(K1)=(STRESS(K1)-ALPHA(K1))/SMISES

      END DO

C

C     SLOVE FOR EQUIVALENT PLASTIC STRAIN INCREAMENT 

C

      DEQPL=(SMISES-SYIELD)/(EG3+HARD)

C

C     UPDATE SHIFT TENSOR (BACK STRESS), ELASTIC AND PLASTIC STRAIN AND STRESS

C

      DO K1=1, NDI

          ALPHA(K1)=ALPHA(K1)+HARD*FLOW(K1)*DEQPL

          EPLAS(K1)=EPLAS(K1)+THREE/TWO*FLOW(K1)*DEQPL

          EELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPL

          STRESS(K1)=ALPHA(K1)+FLOW(K1)*SYIELD+SHYDRO

      END DO

      DO K1=NDI+1, NTENS

           ALPHA(K1)=ALPHA(K1)+HARD*FLOW(K1)*DEQPL

          EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL

          EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL

          STRESS(K1)=ALPHA(K1)+FLOW(K1)*SYIELD

      END DO

C

C     CALCULATE PLASTIC DISSIPATION 

C

      SPD=ZERO

      DO K1=1,NTENS 

          SPD=SPD+(STRESS(K1)+OLDS(K1))*(EPLAS(K1)-OLDPL(K1))/TWO

      END DO

C

C     FORMULATE THE JACOBIAN (MATERIAL TANGENT)

C     FIRST CALCUALTE EFFECTIVE MODULI

C

      EFFG=EG*(SYIELD+HARD*DEQPL)/SMISES

      EFFG2=TWO*EFFG

      EFFG3=THREE*EFFG

      EFFLAM=(EBULK3-EFFG2)/THREE

      EFFHRD=EG3*HARD/(EG3+HARD)-EFFG3 

      DO K1=1,NDI

          DO K2=1,NDI

              DDSDDE(K2,K1)=EFFLAM

          END DO

          DDSDDE(K1,K1)=EFFG2+EFFLAM

      END DO

      DO K1=NDI+1, NTENS

          DDSDDE(K1,K1)=EFFG

      END DO

      DO K1=1, NTENS

          DO K2=1,NTENS

              DDSDDE(K2,K1)=DDSDDE(K2,K1)+EFFHRD*FLOW(K2)*FLOW(K1)

          END DO

      END DO 

      ENDIF

C

C     STORE ELASTIC STRAINS, PLASTIC STRAINS AND SHIFT TENSOR

C     IN STATE VARIABLE ARRAY

C

      DO K1=1, NTENS

          STATEV(K1)=EELAS(K1)

          STATEV(K1+NTENS)=EPLAS(K1)

          STATEV(K1+2*NTENS)=ALPHA(K1)

      END DO

C     

C     OUTPUT REQUEST

      IF (NOEL.EQ.1 .AND. NPT.EQ.1) THEN

          WRITE(*,*)

          WRITE(*,*) '################################################'

          WRITE(*,*) 'ALPHA',ALPHA

          WRITE(*,*) 'EPLAS',EPLAS

          WRITE(*,*) 'EELAS',EELAS

          WRITE(*,*) 'DDSDDE',DDSDDE

          WRITE(*,*) 'FLOW',FLOW

          WRITE(*,*) '################################################'

          ENDIF

      RETURN

      END