[Solved] VUMAT EXAMPLES for Hyperelastic materials

I forgot to add the compressibility term which is required for VUMAT. The form is:

W = K1(I1-3)^m + K2(I1-3)^p + K3(I1-3)^q + 1/D1(J-1)^2

Hi,

Thank you very much for your kind offer. I am quite confident with VUMAT now so all good!
Congratulation with your publication!

Regards,

Willsen

Hi Somaye,

Recalling my experience on how to learn VUMAT, I did a lot of in depth study on the basic of Continuum Mechanics. Through a proper understanding of deformation gradient, different strain and stress measures in large deformation, it becomes really straight forward to understand the different physical quantities available in VUMAT. Jorgens VUMAT for Neo Hookean material was super valuable to me.
My suggestion will be to understand his VUMAT in depth and apply it in a single element test. Do an identical test using ABAQUS in-built Neo-Hookean material model and compare the results.
Build upon that. Try to write a more complicated VUMAT to describe Mooney-Rivlin material and compare it with the equivalent mateiral model from ABAQUS, and so on.

Attached is an example of VUMAT for Mooney Rivlin model:

! This is a VUMAT subroutine for Mooney-Rivlin

subroutine vumat(nblock, ndi, nshr, nstatev, nfieldv, nprops, &
lanneal, stepTime, totTime, dt, cmname, coordMp, &
charLen, props, density, strainInc, relSpinInc, temp0, &
U0, F0, field0, stress0, state0, &
intEne0, inelaEn0, temp1, U1, &
F1, field1, stress1, state1, intEne1, inelaEn1)
implicit none

integer :: nblock, ndi, nshr, nstatev, nfieldv, nprops, lanneal
real :: stepTime, totTime, dt
character*80 :: cmname
real :: coordMp(nblock,*)
real :: charLen, props(nprops), density(nblock), &
strainInc(nblock,ndi+nshr), relSpinInc(nblock,nshr), &
temp0(nblock), U0(nblock,ndi+nshr), &
F0(nblock,ndi+nshr+nshr), field0(nblock,nfieldv), &
stress0(nblock,ndi+nshr), state0(nblock,nstatev), &
intEne0(nblock), inelaEn0(nblock), temp1(nblock), &
U1(nblock,ndi+nshr), F1(nblock,ndi+nshr+nshr), &
field1(nblock,nfieldv), stress1(nblock,ndi+nshr), &
state1(nblock,nstatev), intEne1(nblock), inelaEn1(nblock)

! local variables
real :: J, U(3,3), U2dev(3,3), U4dev(3,3), C10, C01, D1 ,traceU2dev, traceU4dev , constant1, constant2
integer :: i
C10 = props(1)
C01 = props(2)
D1 = props(3)

! loop through all blocks
do i = 1, nblock
! setup F (upper diagonal part)
U(1,1) = U1(i,1)
U(2,2) = U1(i,2)
U(3,3) = U1(i,3)
U(1,2) = U1(i,4)
U(2,1) = U(1,2)
if (nshr .eq. 1) then
U(2,3) = 0.0
U(1,3) = 0.0
else
U(2,3) = U1(i,5)
U(1,3) = U1(i,6)
end if
U(3,2)=U(2,3)
U(3,1)=U(1,3)
! calculate J
J = (U(1,1)*(U(2,2)*U(3,3)-U(2,3)**2)) + (U(1,2)*(U(2,3)*U(1,3)-U(1,2)*U(3,3))) + (U(1,3)*(U(1,2)*U(2,3)-U(2,2)*U(1,3)))

! Define the square of the deviatoric stretch tensor
U2dev=(J**(-2.0/3.0))*matmul(U,U)

! Define the forth power of the deviatoric stretch tensor
U4dev=matmul(U2dev,U2dev)

! Define some constants to simplify the stress expression
traceU2dev=(U2dev(1,1)+U2dev(2,2)+U2dev(3,3))
traceU4dev=(U4dev(1,1)+U4dev(2,2)+U4dev(3,3))
constant1=2.0*C10/J
constant2=2.0*C01/J

!The corotational stress
stress1(i,1) = (constant1+constant2*traceU2dev)*(U2dev(1,1)-1.0/3.0*traceU2dev) - constant2*(U4dev(1,1)-1.0/3.0*traceU4dev) + 2.0/D1*(J-1.0)
stress1(i,2) = (constant1+constant2*traceU2dev)*(U2dev(2,2)-1.0/3.0*traceU2dev) - constant2*(U4dev(2,2)-1.0/3.0*traceU4dev) + 2.0/D1*(J-1.0)
stress1(i,3) = (constant1+constant2*traceU2dev)*(U2dev(3,3)-1.0/3.0*traceU2dev) - constant2*(U4dev(3,3)-1.0/3.0*traceU4dev) + 2.0/D1*(J-1.0)
stress1(i,4) = (constant1+constant2*traceU2dev)*U2dev(1,2) - constant2*U4dev(1,2)
if (nshr .eq. 3) then
stress1(i,5) = (constant1+constant2*traceU2dev)*U2dev(2,3) - constant2*U4dev(2,3)
stress1(i,6) = (constant1+constant2*traceU2dev)*U2dev(3,1) - constant2*U4dev(3,1)
end if
end do
end subroutine vumat

I appreciate this forum.

I am new to mechanics of Polymers but I have to say that I'm catching up quite fast thanks to Jorgen's book, Mechanics of Solid Polymers.

I am trying to implement the Ne0-Hookean VUMAT model in the book and compare with the inbuilt one. I am facing some problems.

1. There is excessive distortion of elements. At the moment, I have no idea about what might be causing the error. Can anyone advise?

2. Jorgen mentions that the stress has to be returned in a co-rotational frame as given in equation 10.3 but the stress formula in the code uses the deviatoric part of the Left Cauchy-Green deformation tensor. Please let me understand this.

If anyone can share a working Neo-Hookean VUMAT, I'd really appreciate.

Thank you

Melly

Just as a reference, here is the VUMAT that I wrote for the Neo-Hookean material model.

c FILE:
c vumat_nh.f
c AUTHOR:
c Jorgen Bergstrom (jorgen@polymerFEM.com)
c CONTENTS:
c A single-precision Abaqus VUMAT subroutine for the
c Neo-Hookean (NH) model. The subroutine is an example
c of how to write a VUMAT, and not is not designed to
c be a commercial quality implementation.
c CREATED:
c 2004-02-17
c USAGE:
c 2 material constants:
c 1: mu (shear modulus)
c 2: kappa (bulk modulus)
c
c |<- column 1 begins here
c |
c *User material, constants=2
c **..:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
c ** mu, kappa
c 5.0, 100.0
c *Density
c 1000e-12
c
c NOTES:
c The subroutine only works for plane-strain, axisymmetric,
c and 3D-elements.
c
c LICENSE:
c
c This program is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 2 of the License, or
c (at your option) any later version.
c
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c GNU General Public License for more details.
c
c You should have received a copy of the GNU General Public License
c along with this program; if not, write to the Free Software
c Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
c
c Copyright 2004 Jorgen Bergstrom
c
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, mu, kappa
integer i

mu = props(1)
kappa = props(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
t1 = (B(1,1) + B(2,2) + B(3,3)) / 3.0
t2 = mu/J
t3 = kappa * log(J) / J
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