VUMAT EXAMPLES for Hyperelastic materials

[willsen]

Hello,

I am new with VUMAT. Recently I have finished reading VUMAT example for Neo-Hookean material from the "Mechanics of Solid Polymers" textbook. Could anyone share more VUMAT examples of Hyperelastic materials for learning purposes? Thank you very much.

Regards,

Willsen

[twhohen]

Hello,

I have written a hyperelastic VUMAT that can handle strain-energy functions of the form:

W = K1(I1-3)^m + K2(I1-3)^p + K3(I1-3)^q

which covers the Neo-Hookean, Yeoh, and (almost) Davies-De-Thomas (1994) models. I think it will be published in an appendix of a related work soon. If you want to take a look before then, let me know and I may be allowed to share it with you before it hits the press.

Regards,
Travis

[twhohen]

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

[willsen]

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

[sjafari63]

Hi,

I am new to VUMAT too and appreciate if anyone share me a VUMAT subroutine and its CAE file so I can go through it carefully.

Thanks so much,
Somaye

[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. Jorgen's 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

[sjafari63]

Thank you so much Willsen!
Can the Jorgen's VUMAT for Neo Hookean material or yours be run as .f or .for file in ABAQUS? It seems those formats are kind of different from the umat file I see here :
https://abaqus-docs.mit.edu/2017/Eng...ces/umathrt2.f

[sjafari63]

Hi There,
I am having a problem with inputting my Vumat into Abaqus explicit.
I have already implemented many subroutines into Abaqus and got results. However, my recent vumat gives this error :"problem during compilation"
my vumat is just simple neo-hookean model.
I appreciate if someone helps me with that. the code is written here:

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
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
parameter( zero = 0.D0, one = 1.D0, two = 2.D0, three = 3.D0,
1 third = one/three, half = 0.5D0, twoThirds = two/three,
2 threeHalfs = 1.5D0 )
C
C elastic constants
C
c10 = props(1)
d1 = props(2)
shrMod = two*c10
blkMod = two/d1
C
C
do k = 1,nblock
C
C calculate left cauchy-green tensor in terms of stretch tensor
C and jacobian
bxx =stretchNew(k,1) * stretchNew(k,1)
1 + stretchNew(k,4) * stretchNew(k,4)
byy =stretchNew(k,4) * stretchNew(k,4)
1 + stretchNew(k,2) * stretchNew(k,2)
bzz =stretchNew(k,3) * stretchNew(k,3)
bxy =stretchNew(k,1) * stretchNew(k,4)
1 + stretchNew(k,4) * stretchNew(k,2)
bxz = zero
byz = zero
detu=stretchNew(k,3)*(stretchNew(k,1)*stretchNew(k ,2)
1 - stretchNew(k,4)*stretchNew(k,4))
C
if ( nshr .gt. 1 ) then
bxx =bxx + stretchNew(k,6) * stretchNew(k,6)
byy =byy + stretchNew(k,5) * stretchNew(k,5)
bzz =bzz + stretchNew(k,6) * stretchNew(k,6)
1 + stretchNew(k,5) * stretchNew(k,5)
bxy =bxy + stretchNew(k,6) * stretchNew(k,6)
1 + stretchNew(k,4) * stretchNew(k,2)
bxz = stretchNew(k,1) * stretchNew(k,6)
1 + stretchNew(k,4) * stretchNew(k,5)
+ stretchNew(k,6) * stretchNew(k,3)
byz = stretchNew(k,4) * stretchNew(k,6)
1 + stretchNew(k,2) * stretchNew(k,5)
2 + stretchNew(k,5) * stretchNew(k,3)
detu = detu + stretchNew(k,6) *
1 ( stretchNew(k,4) * stretchNew(k,5)
2 - stretchNew(k 6) * stretchNew(k,2) )
3 - stretchNew(k,5) * ( stretchNew(k,1) * stretchNew(k,5)
4 - stretchNew(k,6) * stretchNew(k,4) )
end if
xpow = exp ( - log(detu) * two_thirds )
bxx = bxx * xpow
byy = byy * xpow
bzz = bzz * xpow
bxy = bxy * xpow
bxz = bxz * xpow
byz = byz * xpow
C BI1 ( first invariant of BIJ )
C
bi1 = bxx + byy + bzz
C
C BDIJ ( deviatoric part of BIJ )
C
bdxx = bxx - third * bi1
bdyy = byy - third * bi1
bdyy = byy - third * bi1
bdzz = bzz - third * bi1
bdxy = bxy
bdxz = bxz
bdyz = byz
C
duDi1 = half * shrMod
duDi3 = blkMod * ( detu - one )
C
C Calculate Cauchy (true) stresses
C
detuInv = one / detu
factor = two * duDi1 * detuInv
C
stressNew(k,1) = factor * bdxx + duDi3
stressNew(k,2) = factor * bdyy + duDi3
stressNew(k 3) = factor * bdzz + duDi3
stressNew(k,3) = factor * bdzz + duDi3
stressNew(k,4) = factor * bdxy
if ( nshr .gt. 1 ) then
stressNew(k,5) = factor * bdyz
stressNew(k,6) = factor * bdxz
end if
C
end do
C
return
end




Thanks!