Can anyone help me check my VUMAT code for Neo-Hookean material?

Hi,

I am new to writing VUMAT. I tested my Neo-Hookean hyperelastic material and the build-in Abaqus material on single element model but got quite different results. Can anyone take a look at my code and let me know if there is anything wrong?

The code was based on a previous example of Jorgen's VUMAT_NH, but I used a different constitutive model to match the build-in abaqus material.

Any comments would be greatly appreciated!

 

-------------------------------------------------------------

       subroutine vumat(nblock, ndi, nshr, nstatev, nfieldv, nprops,

     .     lanneal, stepTime, totTime, dt, cmname, coordMp,

     .     charLen, props, density, Dstrain, rSpinInc, temp0,

     .     U0, F0, field0, stressVec0, state0,

     .     intEne0, inelaEn0, temp1, U1,

     .     F1, field1, stressVec1, state1, intEne1, inelaEn1)

      implicit none

      integer nblock, ndi, nshr, nstatev, nfieldv, nprops, lanneal

      real stepTime, totTime, dt

      character*80 cmname

      real coordMp(nblock,*)

      integer charLen

      real props(nprops), density(nblock),

     .     Dstrain(nblock,ndi+nshr), rSpinInc(nblock,nshr),

     .     temp0(nblock), U0(nblock,ndi+nshr),

     .     F0(nblock,ndi+nshr+nshr), field0(nblock,nfieldv),

     .     stressVec0(nblock,ndi+nshr), state0(nblock,nstatev),

     .     intEne0(nblock), inelaEn0(nblock), temp1(nblock),

     .     U1(nblock,ndi+nshr), F1(nblock,ndi+nshr+nshr),

     .     field1(nblock,nfieldv), stressVec1(nblock,ndi+nshr),

     .     state1(nblock,nstatev), intEne1(nblock), inelaEn1(nblock)

C

C     local variables

      real F(3,3), B(3,3), J, t1, t2, t3, C10, D1

      integer i

C

      C10 = props(1)

      D1  = props(2)

C

C     loop through all blocks

      do i = 1, nblock

C        setup F (upper diagonal part)

         F(1,1) = U1(i,1)

         F(2,2) = U1(i,2)

         F(3,3) = U1(i,3)

         F(1,2) = U1(i,4)

         if (nshr .eq. 1) then

            F(2,3) = 0.0

            F(1,3) = 0.0

         else

            F(2,3) = U1(i,5)

            F(1,3) = U1(i,6)

         end if

C

C        calculate J

         t1 = F(1,1) * (F(2,2)*F(3,3) - F(2,3)**2)

         t2 = F(1,2) * (F(2,3)*F(1,3) - F(1,2)*F(3,3))

         t3 = F(1,3) * (F(1,2)*F(2,3) - F(2,2)*F(1,3))

         J = t1 + t2 + t3

         t1 = J**(-2.0/3.0)

C

C        Bstar = J^(-2/3) F Ft   (upper diagonal part)

         B(1,1) = t1*(F(1,1)**2 + F(1,2)**2 + F(1,3)**2)

         B(2,2) = t1*(F(1,2)**2 + F(2,2)**2 + F(2,3)**2)

         B(3,3) = t1*(F(1,3)**2 + F(2,3)**2 + F(3,3)**2)

         B(1,2) = t1*(F(1,1)*F(1,2)+F(1,2)*F(2,2)+F(1,3)*F(2,3))

         B(2,3) = t1*(F(1,2)*F(1,3)+F(2,2)*F(2,3)+F(2,3)*F(3,3))

         B(1,3) = t1*(F(1,1)*F(1,3)+F(1,2)*F(2,3)+F(1,3)*F(3,3))

C

C        Stress = 2/J*C10* Dev(Bbar) + 2/D1*(J-1)*Eye

         t1 = (B(1,1) + B(2,2) + B(3,3)) / 3.0

         t2 = 2.0 * C10 / J

         t3 = 2.0 * (J-1) / D1

         StressVec1(i,1) = t2 * (B(1,1) - t1) + t3

         StressVec1(i,2) = t2 * (B(2,2) - t1) + t3

         StressVec1(i,3) = t2 * (B(3,3) - t1) + t3

         StressVec1(i,4) = t2 * B(1,2)

         if (nshr .eq. 3) then

            StressVec1(i,5) = t2 * B(2,3)

            StressVec1(i,6) = t2 * B(1,3)

         end if

C

      end do

      return

      end