umat二次开发超弹性本构.docx

umat二次开发超弹性本构.docx

《umat二次开发超弹性本构.docx》由会员分享，可在线阅读，更多相关《umat二次开发超弹性本构.docx（15页珍藏版）》请在冰点文库上搜索。

umat二次开发超弹性本构

APPENDIX

Neo-HookeanHyperelaticMaterialUserSubroutine

Thisprogramisbasedonthederivationofhyperelasticmaterialconstitutivemodelin

Section4.4.Astressandstrainrelationshipwasderivedfromtheneo-Hookean

hyperelasticmaterialconstitutivemodelthatisnormallyrepresentedasthestrainenergy

withstraininvariants.

subroutinevumat（

CReadonly（unmodifiable）variables-

1nblock,ndir,nshr,nstatev,nfieldv,nprops,lanneal,

2stepTime,totalTime,dt,cmname,coordMp,charLength,

3props,density,strainInc,relSpinInc,

4tempOld,stretchOld,defgradOld,fieldOld,

5stressOld,stateOld,enerInternOld,enerInelasOld,

6tempNew,stretchNew,defgradNew,fieldNew,

CWriteonly（modifiable）variables-

7stressNew,stateNew,enerInternNew,enerInelasNew）

C

include'vaba_param.inc'

C

dimensionprops（nprops）,density（nblock）,coordMp（nblock,*）,

1charLength（nblock）,strainInc（nblock,ndir+nshr）,

2relSpinInc（nblock,nshr）,tempOld（nblock）,

3stretchOld（nblock,ndir+nshr）,

4defgradOld（nblock,ndir+nshr+nshr）,

5fieldOld（nblock,nfieldv）,stressOld（nblock,ndir+nshr）,

6stateOld（nblock,nstatev）,enerInternOld（nblock）,

7enerInelasOld（nblock）,tempNew（nblock）,

8stretchNew（nblock,ndir+nshr）,

8defgradNew（nblock,ndir+nshr+nshr）,

9fieldNew（nblock,nfieldv）,

1stressNew（nblock,ndir+nshr）,stateNew（nblock,nstatev）,

2enerInternNew（nblock）,enerInelasNew（nblock）

C

character*80cmname

C

if（cmname（1:

6）.eq.'VUMAT0'）then

callVUMAT0（nblock,ndir,nshr,nstatev,nfieldv,nprops,lanneal,

2stepTime,totalTime,dt,cmname,coordMp,charLength,

3props,density,strainInc,relSpinInc,

4tempOld,stretchOld,defgradOld,fieldOld,

5stressOld,stateOld,enerInternOld,enerInelasOld,

6tempNew,stretchNew,defgradNew,fieldNew,

7stressNew,stateNew,enerInternNew,enerInelasNew）

117

elseif（cmname（1:

6）.eq.'VUMAT1'）then

callVUMAT1（nblock,ndir,nshr,nstatev,nfieldv,nprops,lanneal,

2stepTime,totalTime,dt,cmname,coordMp,charLength,

3props,density,strainInc,relSpinInc,

4tempOld,stretchOld,defgradOld,fieldOld,

5stressOld,stateOld,enerInternOld,enerInelasOld,

6tempNew,stretchNew,defgradNew,fieldNew,

7stressNew,stateNew,enerInternNew,enerInelasNew）

endif

end

C

subroutinevumat0（

CReadonly-

*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,

CWriteonly-

*stressNew,stateNew,enerInternNew,enerInelasNew）

C

include'vaba_param.inc'

C

dimensioncoordMp（nblock,*）,charLength（nblock）,props（nprops）,

1density（nblock）,strainInc（nblock,ndir+nshr）,

2relSpinInc（nblock,nshr）,tempOld（nblock）,

3stretchOld（nblock,ndir+nshr）,

4defgradOld（nblock,ndir+nshr+nshr）,

5fieldOld（nblock,nfieldv）,stressOld（nblock,ndir+nshr）,

6stateOld（nblock,nstatev）,enerInternOld（nblock）,

7enerInelasOld（nblock）,tempNew（nblock）,

8stretchNew（nblock,ndir+nshr）,

9defgradNew（nblock,ndir+nshr+nshr）,

1fieldNew（nblock,nfieldv）,

2stressNew（nblock,ndir+nshr）,stateNew（nblock,nstatev）,

3enerInternNew（nblock）,enerInelasNew（nblock）

C

dimensiondevia（nblock,ndir+nshr）,

1BBar（nblock,4）,stretchNewBar（nblock,4）,intv

（2）

C

character*80cmname

parameter（zero=0.D00,one=1.D00,two=2.D00,three=3.D00,

*four=4.D00,half=0.5D0）

realC10,D1,ak,twomu,amu,alamda,hydro,vonMises,maxShear,

1midStrain,maxPrincipalStrain

118

C

intv

（1）=ndir

intv

（2）=nshr

C

if（ndir.ne.3.or.nshr.ne.1）then

callxplb_abqerr（1,'SubroutineVUMATisimplemented'//

*'onlyforplanestrainandaxisymmetriccases'//

*'（ndir=3andnshr=1）',0,zero,''）

callxplb_abqerr（-2,'SubroutineVUMAThasbeencalled'//

*'withndir=%Iandnshr=%I',intv,zero,''）

callxplb_exit

endif

C

C10=props

（1）

D1=props

（2）

CC10=1.11619E6D1=4.48E-8

C

ak=two/D1

amu=two*C10

twomu=four*C10

alamda=（three*ak-twomu）/three

C

CifstepTimeequalszero,assumepureelasticmaterialanduseinitialelastic

modulus

C

if（stepTime.EQ.zero）then

dok=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）

Cwrite（6,*）totalTime,k,defgradNew（k,1）,stretchNew（k,1）,

C1stressNew（k,2）,stressNew（k,3）,stressNew（k,4）

enddo

else

dok=1,nblock

C

CJACOBIANOFSTRETCHTENSOR（Uissymmetricandinlocal

axis）

C

det=stretchNew（k,3）*

119

1（stretchNew（k,1）*stretchNew（k,2）-stretchNew（k,4）**two）

scale=det**（-ONE/THREE）

stretchNewBar（k,1）=stretchNew（k,1）*scale

stretchNewBar（k,2）=stretchNew（k,2）*scale

stretchNewBar（k,3）=stretchNew（k,3）*scale

stretchNewBar（k,4）=stretchNew（k,4）*scale

C

CCALCULATELEFTCAUCHY-GREENTENSOR（Bis

symmetric）

C

BBar（k,1）=stretchNewBar（k,1）**two+stretchNewBar（k,4）**two

BBar（k,2）=stretchNewBar（k,2）**two+stretchNewBar（k,4）**two

BBar（k,3）=stretchNewBar（k,3）**two

BBar（k,4）=stretchNewBar（k,1）*stretchNewBar（k,4）+

1stretchNewBar（k,2）*stretchNewBar（k,4）

C

CCALCULATESTRESStensor

C

TRBBar=BBar（k,1）+BBar（k,2）+BBar（k,3）

EG=two*C10/det

PR=two/D1*（det-one）

stressNew（k,1）=EG*（BBar（k,1）-TRBBar/Three）+PR

stressNew（k,2）=EG*（BBar（k,2）-TRBBar/Three）+PR

stressNew（k,3）=EG*（BBar（k,3）-TRBBar/Three）+PR

stressNew（k,4）=EG*BBar（k,4）

C

CUpdatethespecificinternalenergy

C

stressPower=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/density（k）

C

CStrainsundercorotationalcoordinates

C

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）

C

CCalculatevonMises

C

hydro=（stressNew（k,1）+stressNew（k,2）+

120

1stressNew（k,3））/3.

dok1=1,ndir

devia（k,k1）=stressNew（k,k1）-hydro

enddo

dok1=ndir+1,ndir+nshr

devia（k,k1）=stressNew（k,k1）

enddo

vonMises=0.

dok1=1,ndir

vonMises=vonMises+devia（k,k1）**2

enddo

dok1=ndir+1,ndir+nshr

vonMises=vonMises+2*devia（k,k1）**2

enddo

vonMises=sqrt（3./2*vonMises）

Cuse3/2willget2（int）!

C

Cwrite（6,*）totalTime,defgradNew（k,4）,stretchNew（k,4）

C1,defgradNew（k,3）,defgradNew（k,4）,defgradNew（k,5）

C,det,TRBBar

C1,stressNew（k,1）,stressNew（k,2）,stressNew（k,3）,stressNew（k,4）

C

CFailureCriteria

C

midStrain=stateNew（k,1）+stateNew（k,2）

maxShear=sqrt（（stateNew（k,1）-midStrain）**2.+

1stateNew（k,4）**2.）

if（midStrain.GE.0.）then

maxPrincipalStrain=midStrain+maxShear

else

maxPrincipalStrain=maxShear-midStrain

endif

if（vonMises.GE.10.8565e6）then

stateNew（k,5）=0

endif

enddo

endif

return

end

C

subroutinevumat1（

CReadonly-

*nblock,ndir,nshr,nstatev,nfieldv,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,

CWriteonly-

*stressNew,stateNew,enerInternNew,enerInelasNew）

C

include'vaba_param.inc'

dimensioncoordMp（nblock,*）,charLength（nblock）,props（nprops）,

1density（nblock）,strainInc（nblock,ndir+nshr）,

2relSpinInc（nblock,nshr）,tempOld（nblock）,

3stretchOld（nblock,ndir+nshr）,

4defgradOld（nblock,ndir+nshr+nshr）,

5fieldOld（nblock,nfieldv）,stressOld（nblock,ndir+nshr）,

6stateOld（nblock,nstatev）,enerInternOld（nblock）,

7enerInelasOld（nblock）,tempNew（nblock）,

8stretchNew（nblock,ndir+nshr）,

9defgradNew（nblock,ndir+nshr+nshr）,

1fieldNew（nblock,nfieldv）,

2stressNew（nblock,ndir+nshr）,stateNew（nblock,nstatev）,

3enerInternNew（nblock）,enerInelasNew（nblock）

C

dimensiondevia（nblock,ndir+nshr）,

1BBar（nblock,4）,stretchNewBar（nblock,4）,intv

（2）

C

character*80cmname

parameter（zero=0.D00,one=1.D00,two=2.D00,three=3.D00,

*four=4.D00,half=0.5D0）

realC10,D1,ak,twomu,amu,alamda,hydro,vonMises

C

intv

（1）=ndir

intv

（2）=nshr

C

if（ndir.ne.3.or.nshr.ne.1）then

callxplb_abqerr（1,'SubroutineVUMATisimplemented'//

*'onlyforplanestrainandaxisymmetriccases'//

*'（ndir=3andnshr=1）',0,zero,''）

callxplb_abqerr（-2,'SubroutineVUMAThasbeencalled'//

*'withndir=%Iandnshr=%I',intv,zero,''）

callxplb_exit

endif

C

C10=props

（1）

D1=props

（2）

CC10=1.11619E6D1=4.48E-8

C

ak=two/D1

amu=two*C10

122

twomu=four*C10

alamda=（three*ak-twomu）/three

C

CifstepTimeequalszero,assumepureelasticmaterialanduseinitialelastic

modulus

C

if（stepTime.EQ.zero）then

dok=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）

Cwrite（6,*）totalTime,k,defgradNew（k,1）,stretchNew（k,1）,

C1stressNew（k,2）,stressNew（k,3）,stressNew（k,4）

enddo

else

dok=1,nblock

C

CJACOBIANOFSTRETCHTENSOR（Uissymmetricandinlocal

axis）

C

det=stretchNew（k,3）*

1（stretchNew（k,1）*stretchNew（k,2）-stretchNew（k,4）**two）

scale=det**（-ONE/THREE）

stretchNewBar（k,1）=stretchNew（k,1）*scale

stretchNewBar（k,2）=stretchNew（k,2）