APPENDIX
Neo-Hookean Hyperelatic Material User Subroutine
This program is based on the derivation of hyperelastic material constitutive model in Section 4.4. A stress and strain relationship was derived from the neo-Hookean
hyperelastic material constitutive model that is normally represented as the strain energy with strain invariants. 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
if (cmname(1:6) .eq. 'VUMAT0') then
call VUMAT0(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) 117
else if (cmname(1:6) .eq. 'VUMAT1') then
call VUMAT1(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) end if end C
subroutine vumat0 ( C Read only -
* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
* stepTime, totalTime, dt, cmname, coordMp, charLength, * props, density, strainInc, relSpinInc,
* tempOld, stretchOld, defgradOld, fieldOld,
* stressOld, stateOld, enerInternOld, enerInelasOld, * tempNew, stretchNew, defgradNew, fieldNew, C Write only -
* stressNew, stateNew, enerInternNew, enerInelasNew ) C
include 'vaba_param.inc' C
dimension coordMp(nblock,*), charLength(nblock), props(nprops), 1 density(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),
9 defgradNew(nblock,ndir+nshr+nshr), 1 fieldNew(nblock,nfieldv),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), 3 enerInternNew(nblock), enerInelasNew(nblock) C
dimension devia(nblock,ndir+nshr),
1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2) C
character*80 cmname
parameter (zero = 0.D00, one = 1.D00, two = 2.D00, three = 3.D00, * four = 4.D00, half = 0.5D0)
real C10,D1,ak,twomu,amu,alamda,hydro,vonMises, maxShear, 1 midStrain, maxPrincipalStrain 118 C
intv(1) = ndir intv(2) = nshr C
if (ndir .ne. 3 .or. nshr .ne. 1) then
call xplb_abqerr(1,'Subroutine VUMAT is implemented '// * 'only for plane strain and axisymmetric cases '// * '(ndir=3 and nshr=1)',0,zero,' ')
call xplb_abqerr(-2,'Subroutine VUMAT has been called '// * 'with ndir=%I and nshr=%I',intv,zero,' ') call xplb_exit end if C
C10 = props(1) D1 = props(2)
C C10=1.11619E6 D1=4.48E-8 C
ak=two/D1 amu=two*C10 twomu=four*C10
alamda=(three*ak-twomu)/three C
C if stepTime equals zero, assume pure elastic material and use initial elastic modulus C
if(stepTime .EQ. zero) then do k=1,nblock
trace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3) stressNew(k,1) = stressOld(k,1)
* + twomu * strainInc(k,1) + alamda * trace1 stressNew(k,2) = stressOld(k,2)
* + twomu * strainInc(k,2) + alamda * trace1 stressNew(k,3) = stressOld(k,3)
* + twomu * strainInc(k,3) + alamda * trace1 stressNew(k,4) = stressOld(k,4) * + twomu * strainInc(k,4)
C write(6,*) totalTime,k,defgradNew(k, 1),stretchNew(k,1), C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4) end do else
do k=1,nblock