VUMAT implementation error in ABAQUS:Excessive Distortion of Elements

 Hi, 

I am trying to implement probabilistic damage model of the dynamic fragmentation of brittle materials,Hild,Advanaces in applied mechanics,2010 
I have got excessive distortion error.Please anyone suggest me how to rectify the error.Please see my subroutine if any error correct me.I also have a doubt how to update stress in VUMAT my model includes anisotropic damage(D1,D2,D3 in principal directions which are found in each time step,alpha is known).Please tell with one example how to update stress. I got the following error: excessive deformation of elements.

      subroutine vumat(   

  subroutine vumat(

     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,

     3  stressOld, stateOld, enerInternOld, enerInelasOld,

     6  tempNew, stretchNew, defgradNew, fieldNew,

     5  stressNew, stateNew, enerInternNew, enerInelasNew )

C

      INCLUDE 'VABA_PARAM.INC'

C

C

C All arrays dimensioned by (*) are not used in this algorithm

      dimension props(nprops), density(nblock),

     1  coordMp(nblock,*),

     2  charLength(*), strainInc(nblock,ndir+nshr),

     3  relSpinInc(*), tempOld(*),

     4  stretchOld(*), defgradOld(*),

     5  fieldOld(*), stressOld(nblock,ndir+nshr),

     6  stateOld(nblock,nstatev), enerInternOld(nblock),

     7  enerInelasOld(nblock), tempNew(*),

     8  stretchNew(*), defgradNew(*), fieldNew(*),

     9  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),

     1  enerInternNew(nblock), enerInelasNew(nblock)

C

      character*80 cmname

C

C

***********************************************************************************

**  VUMAT,PROBABILISTIC DAMAGE MODEL FOR ABAQUS/STANDARD INCORPORATING            **

**  ELASTIC PROPERTIES AND  WEIBULL PARAMETERS                                   **

**                                                                               **

***********************************************************************************

***********************************************************************************

**

**

**

*USER SUBROUTINE

      PARAMETER (M=3,N=3,ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,

     +          SIX=6.D0, NINE=9.D0,XALPHA=0.31,S=4.18)

C

      DIMENSION XIDEN(M,N),XC(6,6),XS(6,6),

     +          STR(M,N),XDAM(nblock,M),AJAC(6,6),STRN(M,N),

     +          XPRIN(nblock,M),XDEF(nblock,M),XPRINN(nblock,M),

     +          XFAIL(nblock),XPRINR(nblock,M)

      dimension  eigVal(nblock,M), eigVec(nblock,M)

      real V,B,RKO,RKT,RKTH,RKF,RK

      integer seed

       seed=7654321

C

C     XDAM=DAMAGE VALUE CORRESPONDING TO EACH PRICIPAL DIRECTION,XPRIN=PRINCIPAL STRESS

C     XDEF=MODIFIED GROWTH DEFECT DENSITY CORRESPONDING TO EACH PRINCIPAL STRESS

C     XPRINN-PRINCIPALSTRESS NEW,XPRIN-PRICIPAL STRESS OLD

C     XPRINR-PRINCIPAL STRESS RATE

C

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

C

C     SPECIFY MATERIAL PROPERTIES

C     ALL DIMENSIONS ARE IN meter

C

       E=props(1)

       XNUE=props(2)

       XDENSITY=props(3)

       XPOROSITY=props(4)

       XM=props(5)

       XMEANFS=props(6)

       XEFFVOL=props(7)

       ALAMBDA=E*XNUE/(ONE+XNUE)/(ONE-TWO*XNUE)

       AMU=E/(ONE+XNUE)/2

       V=SQRT(E/XDENSITY)

       DO I=1,M

       DO J=1,N

        STR(I,J)=ZERO

       END DO

      END DO

      NTENS=ndir+nshr

       DO I=1,NTENS

        DO J=1,NTENS

            AJAC(I,J)=ZERO

        ENDDO

       ENDDO

C     STATE VARIABLES:1 TO 3 ARE DAMAGE VALUE VARYING FROM 0 TO 1

C     STATE VARIABLES:4-FAILURE STRESS

C     STATE VARIABLES:5,6,7-PRINCIPAL STRESSES

C     STATE VARIABLES:8,9,10-DEFECT DENSITY CORRESPONDING TO EACH PRINCIPAL STRESS

       if (stepTime.eq.zero) then

        do 101 i = 1,nblock   

        B=((XEFFVOL)*((XMEANFS)**XM))

        PK=ran(seed)    

      Q=ONE/(ONE-PK)

      XFAIL(i)=(XMEANFS)*((LOG(Q))**(ONE/XM))     

       AJAC(1,1)= (ALAMBDA+TWO*AMU)

       AJAC(2,2)= (ALAMBDA+TWO*AMU)

       AJAC(3,3)= (ALAMBDA+TWO*AMU)

       AJAC(4,4)= AMU

       AJAC(5,5)= AMU

       AJAC(6,6)= AMU

       AJAC(1,2)= ALAMBDA

       AJAC(1,3)= ALAMBDA

       AJAC(2,1)= ALAMBDA

       AJAC(2,3)= ALAMBDA

       AJAC(3,1)= ALAMBDA

       AJAC(3,2)= ALAMBDA

       

        DO j=1,NTENS

        stressNew(i,j)=stressOld(i,j)+

     2   (DOT_PRODUCT(AJAC(j,:),strainInc(i,:)))

       ENDDO      

      stateNew(i,4)=XFAIL(i)     

      seed=seed+1

101     continue

      

      endif

      do 100 i = 1,nblock

C

C     WRITE STRESSES IN MATRIX FORM

C

      DO j = 1,3

      STR(j,j) = stressOld(i,j)

      END DO

      STR(1,2) = stressOld(i,4)

      STR(2,1) = stressOld(i,4)

      STR(1,3) = stressOld(i,5)

      STR(3,1) = stressOld(i,5)

      STR(2,3) = stressOld(i,6)

      STR(3,2) = stressOld(i,6)

      

c     FIND THE EIGEN VALUES OF THE STRESSOLD MATRIX

 

      CALL VSPRINC( NBLOCK, STR, eigVal, ndir, nshr )

C     EVR IS THE PRICIPAL STRESSES

      DO j=1,3

        XPRIN(i,j)=eigVal(i,j)

      ENDDO

       print*, 'old principal stress'

      print*,  XPRIN(i,:) 

C

C     WRITE STRESSES IN MATRIX FORM

C

      DO j = 1,3

      STRN(j,j) = stressNew(i,j)

      END DO

      STRN(1,2) = stressNew(i,4)

      STRN(2,1) = stressNew(i,4)

      STRN(1,3) = stressNew(i,5)

      STRN(3,1) = stressNew(i,5)

      STRN(2,3) = stressNew(i,6)

      STRN(3,2) = stressNew(i,6)     

c     FIND THE EIGEN VALUES OF THE STRESSNEW MATRIX

 

      CALL VSPRINC( NBLOCK, STRN, eigVal, ndir, nshr )

C     EVR IS THE PRICIPAL STRESSES

      DO j=1,3

        XPRINN(i,j)=eigVal(i,j)

      ENDDO

      print*, 'new principal stress'

      print*,  XPRINN(i,:) 

C     FINDING THE PRINCIPAL STRESS RATE

      DO j=1,3

        XPRINR(i,j)=(XPRINN(i,j)-XPRIN(i,j))/dt

      ENDDO

      print*, 'principal stress rate'

      print*,  XPRINR(i,:) 

C     FINDING THE DEFECT DENSITY CORRESPONDING TO EACH PRICIPAL STRESS

C     FINDING THE DAMAGE VALUE CORRESPONDING TO EACH PRINCIPAL STRESS

C     REFER EQUATION 83

 

      DO j=1,3

       IF (XPRINR(i,j).LE.ZERO) THEN

          stateNew(i,j)=stateOld(i,j)

       ELSE

        IF ((XPRIN(i,j).GT.ZERO).AND.(XPRINR(i,j).GT.ZERO))THEN

         IF (XPRIN(i,j).LE.XFAIL(i)) THEN

                stateNew(i,j+7)=ZERO

                stateNew(i,j)=ZERO

            ELSE

             stateNew(i,j+7)=K*((XPRIN(i,j))**XM)

C     XBETA IS THE VALUE OF THE RIGHT HAND SIDE IN THE EQUATION 81

             XBETA=6*S*((0.33*V)**3)*XDEF(i,j)

             RKO=((1-XDAM(i,j))*(XBETA*(totalTime**2)))/2

             RKT=((1-(XDAM(i,j)+((dt*RKO)/2)))

     1         *(XBETA*((totalTime+dt/2)**2))/2)

             RKTH=(1-(XDAM(i,j)+(dt*RKT)/2))

     1         *(XBETA*((totalTime+dt/2)**2))/2

             RKF=(1-(XDAM(i,j)+(dt*RKTH)))

     1         *(XBETA*((totalTime+dt)**2))/2

             RK=RKO+2*RKT+2*RKTH+RKF

             stateNew(i,j)=stateOld(i,j)+((dt/6)*RK)

             ENDIF

            ELSE

                stateNew(i,j+7)=ZERO

                stateNew(i,j)=ZERO

           ENDIF

        ENDIF

      ENDDO

C     STATE VARIABLES:1 TO 3 ARE DAMAGE VALUE VARYING FROM 0 TO 1

C     STATE VARIABLES:4-FAILURE STRESS

        stateNew(i,4)=XFAIL(i)

        stateNew(i,5)=XPRIN(i,1)

        stateNew(i,6)=XPRIN(i,2)

        stateNew(i,7)=XPRIN(i,3)         

C     UPDATING THE STRESS MATRIX:TO DO THAT UPDATE THE COMPLIANCE TENSOR

C     MATRIX AND INVERSE THE COMPLIANCE TENSOR TO OBTAIN THE CONSTITUTIVE TENSOR

C     XC=COMPLIANCE TENSSOR,XS=INVERSE OF COMPLIANCE TENSOR

      XC(1,1)=ONE/E/(1-stateNew(i,1))

      XC(1,2)=(-XNUE)

      XC(1,3)=(-XNUE)

      XC(1,4)=0.0

      XC(1,5)=0.0

      XC(1,6)=0.0

      XC(2,1)=(-XNUE)

      XC(2,2)=ONE/E/(1-stateNew(i,2))

      XC(2,3)=(-XNUE)

      XC(2,4)=0.0

      XC(2,5)=0.0

      XC(2,6)=0.0

      XC(3,1)=(-XNUE)

      XC(3,2)=(-XNUE)

      XC(3,3)=ONE/E/(1-stateNew(i,3))

      XC(3,4)=0.0

      XC(3,5)=0.0

      XC(3,6)=0.0

      XC(4,1)=0.0

      XC(4,2)=0.0

      XC(4,3)=0.0

      XC(4,4)=(ONE+XNUE)/((ONE-stateNew(i,2))**(XALPHA))

     3 /((ONE-stateNew(i,3))**(XALPHA))

      XC(4,5)=0.0

      XC(4,6)=0.0

      XC(5,1)=0.0

      XC(5,2)=0.0

      XC(5,3)=0.0

      XC(5,4)=0.0

      XC(5,5)=(ONE+XNUE)/((ONE-stateNew(i,3))**(XALPHA))

     4 /((ONE-stateNew(i,1))**(XALPHA))

      XC(5,6)=0.0

      XC(6,1)=0.0

      XC(6,2)=0.0

      XC(6,3)=0.0

      XC(6,4)=0.0

      XC(6,5)=0.0

      XC(6,6)=(ONE+XNUE)/((ONE-stateNew(i,1))**XALPHA)

     5 /((ONE-stateNew(i,2))**XALPHA)

C     FINDING XS=INVERSE OF COMPLIANCE TENSOR BY GUASS JORDAN METHOD

      CALL INVERSEA(XC,XS)         

        DO j=1,NTENS

      stressNew(i,j)=stressOld(i,j)+DOT_PRODUCT(XS(j,:),strainInc(i,:))

       ENDDO

100     continue

        return

        end

 

********************END OF VUMAT SUBROUTINE******************************

*************************************************************************

C **   INVERSE OF A 6X6 MATRIX BY GAUSS JORDAN METHOD

*************************************************************************

C **   INVERSE OF A 6X6 MATRIX BY GAUSS JORDAN METHOD

      SUBROUTINE INVERSEA(A,C)

C

      INCLUDE 'VABA_PARAM.INC'

C

      PARAMETER (M=6,N=6)

      DIMENSION A(M,N),B(6,12),C(M,N)

      INTEGER I,J,N

C **  BUILDING AUGUMENTED MATRIX B

      DO 40 J=1,N

        DO 30 I=1,N

            B(I,J)=A(I,J)

30    CONTINUE

40    CONTINUE

      DO I=1,N

        DO J=N+1,2*N

            IF ((J-I).EQ.N) THEN

                B(I,J)=1.0

                ELSE

                    B(I,J)=0.0

            ENDIF

        ENDDO

      ENDDO

      DO 160 I=1,N

        B(I,:)=B(I,:)/B(I,I)

        DO 150 J=1,N

            IF (I.NE.J) THEN

                B(J,:)=B(J,:)-B(J,I)*B(I,:)

            ENDIF

150        CONTINUE

160        CONTINUE

       DO I=1,N

        DO J=N+1,2*N

            C(I,J-6)=B(I,J)

        ENDDO

       ENDDO

       RETURN

       END