ABAQUS-VUMAT初学者用户子程序小例子

1上一次发过《ABAQUS初学者用户子程序小例子》，給学习UMAT的初学者带来了一定的帮助。现在用到VUMAT，发现网上这种小例子很少，关于VUMAT的资料也不多。摸索了一天，做个VUMAT的小例子供大家分享。实例：简单的平面平板拉伸，材料本构模型采用随动强化模型，E＝210000MPa,ν=0.3，fy=200MPa,E’=10000MPa。左端约束，右端施加位移载荷V＝0.2mm。Author:[e-mail]1材料本构模型：在Property中定义材料时，在General下选中UserMaterial，输入210000、0.3、200、10000，以上两个数值代表E、ν，、fy、E’。在用户子程序中代表着PROPS(1)、PROPS(2)、PROPS(3)、PROPS(4)。然后在General下选中Depvar，由于该例子中没有使用状态变量，因此在第一项中输入大于0的数值即可。另外还需要输入密度。2建模大家都会，故省略2ABAQUS中STEP的设置由于VUMAT需要用到Explicit求解，因此需要在step步骤中设置explicit选项，如下图，其设置可以用默认设置。33ABAQUS调用VUMAT用户子程序同UMAT用户子程序的调用方法。在JobManager中点击Edit选项，在General选项的最后一项中选择自己建立好的用户子程序文件。（注：用户子程序文件可以使用文本编辑器进行编辑，当然也可以用Fortran编译器，如果对自己的用户子程序文件的语法不放心，可以先用Fortan编译器进行编译，不过编译前要先建立Project，关于Fortran编译，这里不再介绍）。完成后，点击submit即可进行分析。4结果以下两张图分布是用户子程序和ABAQUS自带的材料模型（Standard求解器）得到的应力云图，可以看出两种图形基本相同。45VUMAT子程序对于初学者来说，需要注意的是，FORTRAN对于程序语言格式上的要求。例如，对于FORTRAN语言，前六个字符必须空出来，等等。检查语法最好的方法就是在FORTRAN编译器上进行编译。SUBROUTINEVUMAT(CReadonly-1nblock,ndir,nshr,nstatev,nfieldv,nprops,lanneal,2stepTime,totalTime,dt,cmname,coordMp,charLength,3props,density,strainInc,relSpinInc,4tempOld,stretchOld,defgradOld,fieldOld,3stressOld,stateOld,enerInternOld,enerInelasOld,6tempNew,stretchNew,defgradNew,fieldNew,CWriteonly-5stressNew,stateNew,enerInternNew,enerInelasNew)Cinclude'vaba_param.inc'CCJ2MisesPlasticitywithkinematichardeningforplaneCstraincase.CElasticpredictor,radialcorrectoralgorithm.CCThestatevariablesarestoredas:CSTATE(*,1)=backstresscomponent11CSTATE(*,2)=backstresscomponent22CSTATE(*,3)=backstresscomponent33CSTATE(*,4)=backstresscomponent125CSTATE(*,5)=equivalentplasticstrainCCCAllarraysdimensionedby(*)arenotusedinthisalgorithmdimensionprops(nprops),density(nblock),1coordMp(nblock,*),2charLength(*),strainInc(nblock,ndir+nshr),3relSpinInc(*),tempOld(*),4stretchOld(*),defgradOld(*),5fieldOld(*),stressOld(nblock,ndir+nshr),6stateOld(nblock,nstatev),enerInternOld(nblock),7enerInelasOld(nblock),tempNew(*),8stretchNew(*),defgradNew(*),fieldNew(*),9stressNew(nblock,ndir+nshr),stateNew(nblock,nstatev),1enerInternNew(nblock),enerInelasNew(nblock)Ccharacter*80cmnameCparameter(zero=0.,one=1.,two=2.,three=3.,1third=one/three,half=.5,twoThirds=two/three,2threeHalfs=1.5)Ce=props(1)xnu=props(2)yield=props(3)hard=props(4)Ctwomu=e/(one+xnu)thremu=threeHalfs*twomusixmu=three*twomualamda=twomu*(e-twomu)/(sixmu-two*e)term=one/(twomu*(one+hard/thremu))con1=sqrt(twoThirds)Cdo100i=1,nblockCCTrialstresstrace=strainInc(i,1)+strainInc(i,2)+strainInc(i,3)sig1=stressOld(i,1)+alamda*trace+twomu*strainInc(i,1)sig2=stressOld(i,2)+alamda*trace+twomu*strainInc(i,2)sig3=stressOld(i,3)+alamda*trace+twomu*strainInc(i,3)sig4=stressOld(i,4)+twomu*strainInc(i,4)CCTrialstressmeasuredfromthebackstress6s1=sig1-stateOld(i,1)s2=sig2-stateOld(i,2)s3=sig3-stateOld(i,3)s4=sig4-stateOld(i,4)CCDeviatoricpartoftrialstressmeasuredfromthebackstresssmean=third*(s1+s2+s3)ds1=s1-smeands2=s2-smeands3=s3-smeanCCMagnitudeofthedeviatorictrialstressdifferencedsmag=sqrt(ds1**2+ds2**2+ds3**2+2.*s4**2)CCCheckforyieldbydeterminingthefactorforplasticity,Czeroforelastic,oneforyieldradius=con1*yieldfacyld=zeroif(dsmag-radius.ge.zero)facyld=oneCCAddaprotectiveadditionfactortopreventadividebyzeroCwhendsmagiszero.Ifdsmagiszero,wewillnothaveexceededCtheyieldstressandfacyldwillbezero.dsmag=dsmag+(one-facyld)CCCalculatedincrementingamma(thisexplicitlyincludestheCtimestep)diff=dsmag-radiusdgamma=facyld*term*diffCCUpdateequivalentplasticstraindeqps=con1*dgammastateNew(i,5)=stateOld(i,5)+deqpsCCDividedgammabydsmagsothatthedeviatoricstressesareCexplicitlyconvertedtotensorsofunitmagnitudeintheCfollowingcalculationsdgamma=dgamma/dsmagCCUpdatebackstressfactor=hard*dgamma*twoThirdsstateNew(i,1)=stateOld(i,1)+factor*ds1stateNew(i,2)=stateOld(i,2)+factor*ds2stateNew(i,3)=stateOld(i,3)+factor*ds37stateNew(i,4)=stateOld(i,4)+factor*s4CCUpdatethestressfactor=twomu*dgammastressNew(i,1)=sig1-factor*ds1stressNew(i,2)=sig2-factor*ds2stressNew(i,3)=sig3-factor*ds3stressNew(i,4)=sig4-factor*s4CCUpdatethespecificinternalenergy-stressPower=half*(1(stressOld(i,1)+stressNew(i,1))*strainInc(i,1)1+(stressOld(i,2)+stressNew(i,2))*strainInc(i,2)1+(stressOld(i,3)+stressNew(i,3))*strainInc(i,3)1+two*(stressOld(i,4)+stressNew(i,4))*strainInc(i,4))CCenerInternNew(i)=enerInternOld(i)C1+stressPower/density(i)CCUpdatethedissipatedinelasticspecificenergy-plasticWorkInc=dgamma*half*(1(stressOld(i,1)+stressNew(i,1))*ds11+(stressOld(i,2)+stressNew(i,2))*ds21+(stressOld(i,3)+stressNew(i,3))*ds31+two*(stressOld(i,4)+stressNew(i,4))*s4)CenerInelasNew(i)=enerInelasOld(i)C1+plasticWorkInc/density(i)100continueCreturnend6INP输入文件*Heading**Jobname:Job-1Modelname:Model-1**Generatedby:Abaqus/CAE6.10-1*Preprint,echo=NO,model=NO,history=NO,contact=NO****PARTS***Part,name=Part-1*Node1,10.,5.2,8.48214245,5.3,6.96428585,5.4,5.44642878,5.5,3.92857146,5.6,2.41071439,5.87,0.892857134,5.8,-0.625,5.9,-2.14285707,5.10,-3.66071439,5.11,-5.17857122,5.12,-6.69642878,5.13,-8.21428585,5.14,-9.73214245,5.15,-11.25,5.16,-12.7678576,5.17,-14.2857141,5.18,-15.8035717,5.19,-17.321