Notifications
Clear all

VUMAT for Hyperelastic Materials - Update to Bergstrom's Example  


twhohen
Posts: 9
(@twhohen)
Active Member
Joined: 2 years ago

Hello,

Several posts discuss VUMATs for hyperelastic materials. Bergstrom (2015) provides an example for the Neo-Hookean strain-energy function (also copied to several forums) that was written for older versions of Abaqus. An updated, complete subroutine is published here:

Hohenberger et al (2019). A constitutive model for both low and high strain nonlinearities in highly filled elastomers and implementation with user-defined material subroutines in Abaqus. Rubber Chemistry and Technology, 92 (4), 653-686.

It uses a modified (generalized) Yeoh strain-energy function and is confirmed to work in Abaqus v6.14 and R2018. I am sharing this here since many posts inquiring about hyperelastic VUMATs are now closed.

Regards,
Travis

5 Replies
Jorgen
Posts: 3874
Moderator
(@jorgen)
Member
Joined: 10 months ago

Hi Travis, that I great. Are you able to also share the source code here?

-Jorgen

Reply
1 Reply
twhohen
(@twhohen)
Joined: 2 years ago

Active Member
Posts: 9

@jorgen

I will check with the publisher. In the meantime, anybody with immediate interest in the code can contact the corresponding author, James Busfield. He is easy to find with a web search.

-th

Reply
twhohen
Posts: 9
(@twhohen)
Active Member
Joined: 2 years ago

Hello,

I have confirmed that I can post the code. For a detailed explanation of the statements, refer to the article above which summarizes the finite elasticity theory that is necessary to build a hyperelastic VUMAT. A discussion of implementation of the theory in Abaqus is given in Appendix C. The article was built largely from the more detailed accounts in the books by Holzapfel (2000) and Bergström (2015).

Some pertinent updates to the Neo-Hookean subroutine that is commonly shared online are:

- updated names for Abaqus variables.
- a gen-Yeoh SEF is used. Set:
----- (K2,K3) = (0,0) for Neo-Hookean SEF.
----- (m,p,q) = (1,2,3) for Yeoh SEF.
----- (K2,K3) = (0,2) for Davies-De-Thomas SEF (1994, Rubber Chem. & Tech. 67 (4), 716).
- includes linear elastic initialization step required per Abaqus VUMAT subroutine template.
- includes logic to handle indeterminate stress condition that can arise with the gen-Yeoh SEF.

The code only handles 3D elements. A 2D plane-strain or axisymmetric (but not plane-stress) version can be written by removing references to the 5th and 6th elements in the stress, strain, and strain increment vectors. Alternatively, the logic in Bergström's Neo-Hookean example can be included to handle 2D and 3D cases.

Regards,
Travis

C **********************************************************************
C VUMAT (3D only) for gen-Yeoh Strain-Energy Function
C Compatibility confirmed for Abaqus v6.14 & R2018
C ---------------------------------------------------
C Authors: Travis Hohenberger & Richard Windslow
C Date: 2019-07-11
C E-mail: twhohen@gmail.com
C Source: Hohenberger et al (2019). Rubber Chem. & Tech., 92 (4), 653
C **********************************************************************
C
C Strain-energy function:
C
C W = K1*(I1-3)^m + K2*(I1-3)^p + K3*(I1-3)^q + (1/D1)*(J-1)^2
C
C where K1, K2, K3, m, p, q are distortional fitting parameters and D1
C is a volumetric fitting parameter. I1 is the first invariant of the
C modified stretch tensor. J is the volumetric ratio.
C
C **********************************************************************
C
SUBROUTINE VUMAT(
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,
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),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
8 stretchNew(nblock,ndir+nshr),
9 defgradNew(nblock,ndir+nshr),
1 fieldNew(nblock,nfieldv),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
3 enerInternNew(nblock), enerInelasNew(nblock)
C
CHARACTER*80 cmname
C
C PARAMETERS
C ----------
REAL*8 oneThrd, half, twoThrd, one, two, three, thresh
PARAMETER(oneThrd=1.d0/3.d0, half=0.5d0, twoThrd=2.d0/3.d0,
$ one=1.d0, two=2.d0, three=3.d0, thresh=10.d0**-12.d0)
C
C LOCAL VARIABLES
C ---------------
REAL*8 g0 , k0 , twoG , lmda , trace , d1 ,
$ k1 , k2 , k3 , em , pe , qu ,
$ Bxx , Byy , Bzz , Bxy , Bxz , Byz ,
$ BbarXX , BbarYY , BbarZZ , BbarXY , BbarXZ , BbarYZ ,
$ dBbarXX , dBbarYY, dBbarZZ, dBbarXY, dBbarXZ, dBbarYZ ,
$ J , J23 , duDi1 , duDi3 , I1 , g1 ,
$ p0 , j1 , u1
C
k1 = props(1)
k2 = props(2)
k3 = props(3)
em = props(4)
pe = props(5)
qu = props(6)
d1 = props(7)
C
g0 = two * k1
k0 = two / d1
C
C ******************************************************************
C ------------ INITIALIZE MATERIAL AS LINEARLY ELASTIC -------------
C ******************************************************************
C
twoG = two * g0
lmda = k0 - twoG * oneThrd
C
IF (totalTime.EQ.0.0) THEN
C
DO k = 1,nblock
trace = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
stressNew(k,1) = stressOld(k,1) + twoG*strainInc(k,1) +
$ lmda*trace
stressNew(k,2) = stressOld(k,2) + twoG*strainInc(k,2) +
$ lmda*trace
stressNew(k,3) = stressOld(k,3) + twoG*strainInc(k,3) +
$ lmda*trace
stressNew(k,4) = stressOld(k,4) + twoG*strainInc(k,4)
stressNew(k,5) = stressOld(k,5) + twoG*strainInc(k,5)
stressNew(k,6) = stressOld(k,6) + twoG*strainInc(k,6)
END DO
C
RETURN
C
END IF
C
C ******************************************************************
C ----------- START LOOP FOR MATERIAL POINT CALCULATIONS -----------
C ******************************************************************
C
DO k = 1,nblock
C
C CALCULATE LEFT CAUCHY-GREEN STRAIN TENSOR, B = U*U
C --------------------------------------------------
Bxx = stretchNew(k,1) * stretchNew(k,1) +
$ stretchNew(k,4) * stretchNew(k,4) +
$ stretchNew(k,6) * stretchNew(k,6)
Byy = stretchNew(k,2) * stretchNew(k,2) +
$ stretchNew(k,4) * stretchNew(k,4) +
$ stretchNew(k,5) * stretchNew(k,5)
Bzz = stretchNew(k,3) * stretchNew(k,3) +
$ stretchNew(k,5) * stretchNew(k,5) +
$ stretchNew(k,6) * stretchNew(k,6)
Bxy = stretchNew(k,1) * stretchNew(k,4) +
$ stretchNew(k,4) * stretchNew(k,2) +
$ stretchNew(k,6) * stretchNew(k,5)
Bxz = stretchNew(k,1) * stretchNew(k,6) +
$ stretchNew(k,4) * stretchNew(k,5) +
$ stretchNew(k,6) * stretchNew(k,3)
Byz = stretchNew(k,4) * stretchNew(k,6) +
$ stretchNew(k,2) * stretchNew(k,5) +
$ stretchNew(k,5) * stretchNew(k,3)
C
C CALCULATE J = |F| = |U|
C -----------------------
C
J = stretchNew(k,1) *
$ ( stretchNew(k,2) * stretchNew(k,3) -
$ stretchNew(k,5) * stretchNew(k,5) ) +
$ stretchNew(k,4) *
$ ( stretchNew(k,5) * stretchNew(k,6) -
$ stretchNew(k,3) * stretchNew(k,4) ) +
$ stretchNew(k,6) *
$ ( stretchNew(k,4) * stretchNew(k,5) -
$ stretchNew(k,2) * stretchNew(k,6) )
C
C CALCULATE MODIFIED STRAIN TENSOR, Bbar = J^(-2/3)*B
C ---------------------------------------------------
J23 = J**(-twoThrd)
C
BbarXX = J23 * Bxx
BbarYY = J23 * Byy
BbarZZ = J23 * Bzz
BbarXY = J23 * Bxy
BbarXZ = J23 * Bxz
BbarYZ = J23 * Byz
C
C FIRST INVARIANT OF Bbar = tr(Bbar)
C ----------------------------------
I1 = BbarXX + BbarYY + BbarZZ
C
C DEVIATORIC PART OF Bbar
C -----------------------
p0 = oneThrd * I1
C
dBbarXX = BbarXX - p0
dBbarYY = BbarYY - p0
dBbarZZ = BbarZZ - p0
dBbarXY = BbarXY
dBbarXZ = BbarXZ
dBbarYZ = BbarYZ
C
C DERIVATIVES OF STRAIN-ENERGY FUNCTION
C -------------------------------------
j1 = I1 - three
C
IF (j1.LT.thresh) THEN
duDi1 = zero
ELSE
duDi1 = em * k1 * j1**(em-one) +
$ pe * k2 * j1**(pe-one) +
$ qu * k3 * j1**(qu-one)
END IF
C
duDi3 = two/d1 * (J - one)
C
C COROTATIONAL CAUCHY (TRUE) STRESSES
C -----------------------------------
g1 = two/J * duDi1
C
stressNew(k,1) = g1 * dBbarXX + duDi3
stressNew(k,2) = g1 * dBbarYY + duDi3
stressNew(k,3) = g1 * dBbarZZ + duDi3
stressNew(k,4) = g1 * dBbarXY
stressNew(k,5) = g1 * dBbarYZ
stressNew(k,6) = g1 * dBbarXZ
C
C UPDATE SPECIFIC INTERNAL ENERGY
C -------------------------------
u1 = half * ( (stressOld(k,1)+stressNew(k,1))*strainInc(k,1) +
$ (stressOld(k,2)+stressNew(k,2))*strainInc(k,2) +
$ (stressOld(k,3)+stressNew(k,3))*strainInc(k,3) +
$ two * ( (stressOld(k,4) + stressNew(k,4))*
$ strainInc(k,4) +
$ (stressOld(k,5) + stressNew(k,5))*
$ strainInc(k,5) +
$ (stressOld(k,6) + stressNew(k,6))*
$ strainInc(k,6) ) )
C
enerInternNew(k) = enerInternOld(k) + u1 / density(k)
C
END DO
C
RETURN
C
END SUBROUTINE VUMAT
Reply
1 Reply
fat32ntfs
(@fat32ntfs)
Joined: 3 years ago

New Member
Posts: 2

@twhohen hi, i think you may have forgot to define zero after the line:

 IF (j1.LT.thresh) THEN
Reply
fat32ntfs
Posts: 2
(@fat32ntfs)
New Member
Joined: 3 years ago

So it seems I cannot edit my post, so I apologize for the double post. I wanted to add the following general question:

In my VUMATs, I use matrix multiplications rather than computing the components individually, for example:

C	F^eE = U
U(1,1) = stretchNew(k,1)
U(2,2) = stretchNew(k,2)
U(3,3) = stretchNew(k,3)
U(1,2) = stretchNew(k,4)
U(2,3) = stretchNew(k,5)
U(1,3) = stretchNew(k,6)
U(2,1) = U(1,2)
U(3,1) = U(1,3)
U(3,2) = U(2,3)
C
C B^ee = F^eE * F^eE^{T} = U^2
B = zero
B = matmul(U,U)

Is there any reason to choose one over the other?

Reply
Share: