# [Solved] VUMAT EXAMPLES for Hyperelastic materials

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

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

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

- VUMAT for Hyperelastic Materials - Update to Bergstrom's Example2 years ago
- Solution Dependent Variables (*DEPVAR) in VUMAT3 years ago
- Deformation Gradient in VUMAT7 years ago
- help with vumat8 years ago
- Material damage model applied in ductile /brittle packaging material11 years ago

- 21 Forums
- 3,870 Topics
- 13.2 K Posts
- 4 Online
- 29.4 K Members