For your prompt reply and advice. I realized that material properties contribute to the excessive distortion. There is a pair of properties that work well for both the inbuilt and the VUMAT producing similar results whereas others only work for the inbuilt model but not VUMAT.

I am working on understanding more about this and improving the model.

Thank you so much.

]]>Abaqus/Explicit expects the stresses to be returned in a co-rotational coordinate system. You can achieve this by either calculating the Cauchy stress based on the deformation gradient and then rotating the stress, or you can achieve the same effect (but more efficiently) if you simply calculate the stress using the left Cauchy-Green tensor U.

-Jorgen

]]>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

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

]]>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 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

]]>