Hi,

I have been recently learning VUMAT to implement hyperelastic materials. I found Dr.Bergstroms earlier post with the example VUMAT code vumat_nh.for for Neo-Hookean Material.

I modified this VUMAT by rewriting the Cauchy stress calculation in the following format, in order to keep it the same as the Abaqus built-in model for comparison:

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

However, when I tested this code on a C3D8R element (four nodes on the same surface fixed, instantaneous force applied to the rest four nodes, same setting for the Abaqus built-in material) , there was quite large a difference. My VUMAT element deformed very quickly and got blown up.

I have been stuck in this for quite a while. Could anyone please take a look at my code and kindly comment where I may get it wrong?

Any comments would be greatly appreciated!! Thank you!

Best,

Yun

Code attached (original copyright information not included)

___________________________________________________________

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 local variables

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

integer i

c mu = props(1)

c kappa = props(2)

C10 = props(1) ! Abaqus built-in NeoHookean model input 1

D1 = props(2) ! Abaqus built-in NeoHookean model input 2

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 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 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 Stress = mu/J * Dev(Bstar) + kappa*log(J)/J * Eye (Jorgen)

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

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

t2 = 2*C10/J

t3 = 2*(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

end do

return

end