umat二次开发超弹性本构

APPENDIXNeo-HookeanHyperelaticMaterialUserSubroutineThisprogramisbasedonthederivationofhyperelasticmaterialconstitutivemodelinSection4.4.Astressandstrainrelationshipwasderivedfromtheneo-Hookeanhyperelasticmaterialconstitutivemodelthatisnormallyrepresentedasthestrainenergywithstraininvariants.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)Cinclude'vaba_param.inc'Cdimensionprops(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)Ccharacter*80cmnameCif(cmname(1:6).eq.'VUMAT0')thencallVUMAT0(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)117elseif(cmname(1:6).eq.'VUMAT1')thencallVUMAT1(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)endifendCsubroutinevumat0(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)Cinclude'vaba_param.inc'CdimensioncoordMp(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)Cdimensiondevia(nblock,ndir+nshr),1BBar(nblock,4),stretchNewBar(nblock,4),intv(2)Ccharacter*80cmnameparameter(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,maxPrincipalStrain118Cintv(1)=ndirintv(2)=nshrCif(ndir.ne.3.or.nshr.ne.1)thencallxplb_abqerr(1,'SubroutineVUMATisimplemented'//*'onlyforplanestrainandaxisymmetriccases'//*'(ndir=3andnshr=1)',0,zero,'')callxplb_abqerr(-2,'SubroutineVUMAThasbeencalled'//*'withndir=%Iandnshr=%I',intv,zero,'')callxplb_exitendifCC10=props(1)D1=props(2)CC10=1.11619E6D1=4.48E-8Cak=two/D1amu=two*C10twomu=four*C10alamda=(three*ak-twomu)/threeCCifstepTimeequalszero,assumepureelasticmaterialanduseinitialelasticmodulusCif(stepTime.EQ.zero)thendok=1,nblocktrace1=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)stressNew(k,1)=stressOld(k,1)*+twomu*strainInc(k,1)+alamda*trace1stressNew(k,2)=stressOld(k,2)*+twomu*strainInc(k,2)+alamda*trace1stressNew(k,3)=stressOld(k,3)*+twomu*strainInc(k,3)+alamda*trace1stressNew(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)enddoelsedok=1,nblockCCJACOBIANOFSTRETCHTENSOR(Uissymmetricandinlocalaxis)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)=stretchNew(k,4)*scaleCCCALCULATELEFTCAUCHY-GREENTENSOR(Bissymmetric)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)+1stretchNewBar(k,2)*stretchNewBar(k,4)CCCALCULATESTRESStensorCTRBBar=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)CCUpdatethespecificinternalenergyCstressPower=half*(1(