Can someone 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?

I knew there was an early example wrote by Jorgen Bergstrom. But the consititutive model was somehow different from the Abaqus build-in one. I would really like to develop my own code-writing skills and confirm it by comparing with Abaqus build-in materials, so I am more confident in writing more challenging code in my future application.

Any comments would be greatly appreciated!

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

     subroutine vumat(

C Read only (unmodifiable)variables -

     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 (modifiable) variables -

     7  stressNew, stateNew, enerInternNew, enerInelasNew )

C

      include 'vaba_param.inc'

C

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),

     4  defgradOld(nblock,ndir+nshr+nshr),

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

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

     7  enerInelasOld(nblock), tempNew(nblock),

     8  stretchNew(nblock,ndir+nshr),

     8  defgradNew(nblock,ndir+nshr+nshr),

     9  fieldNew(nblock,nfieldv),

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

     2  enerInternNew(nblock), enerInelasNew(nblock)

C

      character*80 cmname

C 

C neo-Hookean hyperelasticity for 3D cases 

C 

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

C props(1) - C10 

C props(2) - D1 

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

C     local variables

      real F(3,3), B(3,3), J, scale, trace, C10, D1

C      

      C10    = props(1)

      D1     = props(2)

C

      do i = 1,nblock

C        setup the deformation gradient tensor F

         F(1,1) = defgradNew(i,1)

         F(2,2) = defgradNew(i,2)

         F(3,3) = defgradNew(i,3)

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

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

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

         F(2,1) = defgradNew(i,7)

         F(3,2) = defgradNew(i,8)

         F(1,3) = defgradNew(i,9)

C

C        calculate J

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

     *          - F(1,2)*F(2,1)*F(3,3) + F(1,2)*F(2,3)*F(3,1) 

     *          + F(1,3)*F(2,1)*F(3,2) - F(1,3)*F(2,2)*F(3,1)

         scale  = J**(-2.0/3.0)

C         

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

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

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

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

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

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

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

C      

C        Calculate Tr(Bbar) and Cauchy Stress

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

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

         stressNew(i,1) =  (B(1,1)-trace)*2.0*C10/J + (J-1.0)*2.0/D1

         stressNew(i,2) =  (B(2,2)-trace)*2.0*C10/J + (J-1.0)*2.0/D1

         stressNew(i,3) =  (B(2,2)-trace)*2.0*C10/J + (J-1.0)*2.0/D1

         stressNew(i,4) =          B(1,2)*2.0*C10/J

         stressNew(i,5) =          B(2,3)*2.0*C10/J

         stressNew(i,6) =          B(1,3)*2.0*C10/J 

C          

      end do

C

      return

      end