User login

Navigation

You are here

VUMAT implementation error in ABAQUS:Excessive Distortion of Elements

KOSURI BHANU's picture
AttachmentSize
Image icon COMPLIANCE TENSOR USED27.1 KB

 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

 

 

 

 

 

 

 

Assuming that you are updating the stress correctly, the cause for an "Excessive Distortion of Elements" error is likely that your stiffness is too small. Because you are dealing with damage, it seems likely that as your element damages the stiffness reduces and you get excessive distortion. If you expect that much distortion in your simulation then maybe you should look into element erosion algorithms to remove the elements with too much damage. Another possible fix is to refine your mesh.

For debuging your VUMAT, I would suggest setting the parameters such that no damage occurs and you are left with a linear elastic, isotropic material. Run the simulation and see if you get excessive distortion. Run the simulation with the built-in isotropic, linear elastic material and see if you get the same results. If not, keep debuging.

You asked for an example VUMAT. The best document I've ever found that is an introduction to Abaqus VUMATs/UMATs can be found here:

http://imechanica.org/files/Writing%20User%20Subroutines%20with%20ABAQUS...

Good luck!

Subscribe to Comments for "VUMAT implementation error in ABAQUS:Excessive Distortion of Elements"

More comments

Syndicate

Subscribe to Syndicate