umat二次开发超弹性本构.docx
《umat二次开发超弹性本构.docx》由会员分享，可在线阅读，更多相关《umat二次开发超弹性本构.docx（13页珍藏版）》请在冰点文库上搜索。
umat二次开发超弹性本构
APPENDIX
Neo-HookeanHyperelaticMaterialUserSubroutine
Thisprogramisbasedonthederivationofhyperelasticmaterialconstitutivemodelin
Section.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''
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''
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=,one=,two=,three=,
*four=,half=
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'Q.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.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''
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=,one=,two=,three=,
*four=,half=
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'Q.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）*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）
123
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
升级会员 