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