umat二次開發(fā)超彈性本構(gòu)

1、APPENDIXNeo-Hookean Hyperelatic Material User SubroutineThis 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 s

2、train 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,

3、stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,C Write only (modifiable) variables -7 stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude vaba_param.incCdimension props(nprops), density(nblock), coordMp(nblock,*), 1 charLength(nblock), strainInc(nblock,ndi

4、r+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 de

5、fgradNew(nblock,ndir+nshr+nshr),9 fieldNew(nblock,nfieldv),1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),2 enerInternNew(nblock), enerInelasNew(nblock)Ccharacter*80 cmnameCif (cmname(1:6) .eq. VUMAT0) thencall VUMAT0(nblock, ndir , nshr, nstatev, nfieldv, nprops, lanneal, 2 stepTime, total

6、Time, 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.

7、VUMAT1) thencall 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, defgra

8、dNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)end ifendCsubroutine vumat0 (C Read only -* nblock, ndir, nshr , nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fiel

9、dOld,* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude vaba_param.incCdimension coordMp(nblock,*), charLength(nblock), props(nprops), 1 density(nblock), strainInc(nblock,ndir+nsh

10、r),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 defgrad

11、New(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = 0.D00, one = 1.D00, t

12、wo = 2.D00, three = 3.D00, * four = 4.D00, half = 0.5D0).real C10,D1,ak,twomu,amu,alamda,hydro,vonMises, maxShear,1 midStrain, maxPrincipalStrain118Cintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,Subroutine VUMAT is implemented /* only for plane strain and axis

13、ymmetric 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_exitend ifCC10 = props(1)D1 = props(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10twomu=four*C10alamda=(three*ak-twomu)/threeCC if stepTime equals

14、zero, assume pure elastic material and use initial elastic modulusCif(stepTime .EQ. zero) thendo k=1,nblocktrace1 = 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) +

15、alamda * trace1 stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(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 doelsedo k=1,nblock.CC JACOBIAN OF STRETCH

16、 TENSOR (U is symmetric and in local axis)Cdet=stretchNew(k, 3)*1191 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)*two)scale=det*(-ONE/THREE)stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretch

17、New(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B issymmetric)CBBar(k,1)=stretchNewBar(k, 1)*two+stretchNewBar(k, 4)*twoBBar(k,2)=stretchNewBar(k, 2)*two+stretchNewBar(k, 4)*twoBBar(k,3)=stretchNewBar(k, 3)*twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewB

18、ar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PRstressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PRstressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PRstressNew(k,4)=EG* BBar(k,4)CC Update the specific inte

19、rnal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4)enerInternNew(k) = enerInternOld(k)1 + stressPower

20、 / density(k)CC Strains under corotational coordinatesC.stateNew(k,1) = stateOld(k,1) + strainInc(k,1)stateNew(k,2) = stateOld(k,2) + strainInc(k,2)stateNew(k,3) = stateOld(k,3) + strainInc(k,3)stateNew(k,4) = stateOld(k,4) + strainInc(k,4)CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)

21、+1201 stressNew(k,3)/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydroend dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)*2end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)*2end dovonMises = sqrt(3./2*vonMises)C us

22、e 3/2 will get 2 (int) !CC write(6,*) totalTime,defgradNew(k, 4),stretchNew(k,4)C 1 ,defgradNew(k,3),defgradNew(k,4),defgradNew(k,5)C ,det,TRBBarC 1 ,stressNew(k,1),stressNew(k,2),stressNew(k,3),stressNew(k,4)CC Failure CriteriaCmidStrain = stateNew(k,1) + stateNew(k,2)maxShear = sqrt(stateNew(k,1)

23、- midStrain)*2. +1 stateNew(k,4)*2.)if (midStrain .GE. 0.) thenmaxPrincipalStrain = midStrain + maxShearelsemaxPrincipalStrain = maxShear - midStrainend ifif (vonMises .GE. 10.8565e6) thenstateNew(k,5) = 0end if.end doend ifreturnendCsubroutine vumat1 (C Read only -* nblock, ndir, nshr , nstatev, nf

24、ieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,121* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, en

25、erInternNew, enerInelasNew )Cinclude vaba_param.incdimension 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),

26、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), enerIn

27、elasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (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,vonMisesCintv(1) = ndirintv(2) = nshrCif (ndir

28、.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_exitend ifCC10 = props(1)D1 = props

29、(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10122twomu=four*C10alamda=(three*ak-twomu)/threeCC if stepTime equals zero, assume pure elastic material and use initial elastic modulusCif(stepTime .EQ. zero) thendo k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)stressNew(k,1) = s

30、tressOld(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 * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNe

31、w(k, 1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in local axis)Cdet=stretchNew(k, 3)*1 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)*two)scale=det*(-ONE/THREE).stretchNewBar(k, 1)=stretchNew(k, 1)*s

32、calestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B issymmetric)CBBar(k,1)=stretchNewBar(k, 1)*two+stretchNewBar(k, 4)*twoBBar(k,2)=stretchNewBar(k, 2)*two+stretchNewBar(k, 4)*twoBB

33、ar(k,3)=stretchNewBar(k, 3)*twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)123EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PRstressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PRstressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PRstressNew(k,4)=EG* BBar(k,4)CC Update the spec