ABAQUS非线性有限元分析与实例umat

1、第19章ABAQUS用户材料子程序（UMAT）虽然ABAQUS为用户提供了大量的单元库和求解模型，使用户能够利用这些处理绝大多数的问题；但是现实世界毕竟十分复杂*ABAQUS不可能把所有可能的问题都包含进去。所以ABAQUS提供了大量的用户子程序（user觇子程序允许用户在找不到合适模型的情况下自行定义符合自己问题的模型*这些用程序涵盖了建模从载荷到单元的几乎各个部分。ABAQUS为用户提供的这个接口，允许用户通过自定义的子程序定制ABAQ以实现特定的功能。用户子程序具有以下的功能和特点；1）如果ABAQUS的一些固有选项模型功能有限，用户子程序可以提高ABA，中这些选项的功能。2）通常用户子

2、程序是用Fortran语言的代码写成。3）它可以以几种不同的方式包含在模型中。4）由子它们没有存储在restart文件中，如果需要的话，可以在重新开始运行I改它。5）在某些情况下它可以利用ABAQUS允许的已有程序要在模型中包含用户子程序，可以利用ABAQUS执行程序，在ABAQUS执;序中应用user选项指明包含这些子程序的Fortran源程序或者目标程序的名字提示：ABAQUS的输入文件除了可以通过ABAQUS/CAE的作业摸块提交运古还可以在ABAQUSCommand口中输入ABAQUS执疔程序直接运疔：ABAQUSjob二输入文件名user-用户予程序的Fortran文件名ABAQUS

3、/Standard和ABAQUS/Explicit都支持用户子程序功能，但是它们所的用户子程序种类不尽相同，读者在使用时请注意查询手册。在接下来的两章中，我们将讨论两种常用的用户子程序用户材料子程序和单元子程序。本章将通过在ABAQUS/Standard中创建Johnson-Cook的材料模型，对编写Sdard的用户材料子程序UMAT进行简单介绍。ABAQUS/Explic让中的用户材料子1VUMAT的思想与之相似，但是由于隐式和显式两种方法本身的差异”它们之间一些不同，请读者在自己具体使用前首先仔细查阅ABAQUS手册中的相关内容。19-1引言用户材料子程序是ABAOUS握供给阳户白宗义材料

4、属件的Fnrtrn稈序榕口ao=A+B(19ao=A+B(19Eee=E(Eep)(过ABAQUS/Sundard的用户材料子程序UMAT编程实现。在UMAT编程中使用相关塑性理论以及完全隐式的应力更新算法。19.2模型的数学描述19*2.1Johnson-Cook强化模型简介Johnson-Cook模型用来模拟在冲击载荷作用下的变形。Johnson-Cook强化模型(JC)表示为三项的乘积，分别反映了应变硬化、应变化和温度软化。这里使用JC模型的修正形式：(1a=(A+Ben)1+Cln1+(1一模型中包含A.n.C、椚五个参数，需要通过实验来确定。使参考应变率&这样公式中的A即为材料的静态

5、屈服应力。公式中的:T为无量纲化的温度：T.r-rrTm-TT其中丁为室温,Tm为材料的熔点。Johnson-Cook模型在温度从室温到材料熔点的范围内都是有效的。高应变率的变形经常伴有温升现象，这是因为材料变形过程中塑性功转化为?对于大多数金属，90%-100%的塑性变形将耗散为热量。所以JC模型中温度的喪以用如下的公式计算：丁=2血任(1其中，丁为温度的增量；盘为塑性耗散比，表示塑性功转化为热量的比例；为的比热，p为材料的密度口公式(19-2)考虑的是一个绝热过程，即认为温度的？|全起因于塑性耗散。19.2.2率相关塑性的基本公式JohnsomCook本构模型考虑率相关塑性，塑性变形是关联

6、的，即塑性流动沿看面的法线方向，并采用Mises屈服面。将应变的增量分解为弹性部分和塑性部分：de=dee+dep(1将上式两端同时对时间的增量df微分得到率形式：右二于+尹(1在率相关塑性中，材料的塑性反应取决于加载率，以率的形式给岀材料的弾应；上式即为流动法则，其中誓为塑性流动势能，入为塑性率参数*“为运动硬化时的力。对于各项同性硬化，不存在背应力，因此有0=0,此时有(r誓hI7I二crsign(d-)(I1与率无关塑性不同的是,而是直接通过经验定律给出，sign(7)率相关塑性中等效塑性应变率不能通过一致性条件老成为过应力模型：A(F式中0是过应力，可为粘性。在过应力模型中，等效塑性应

7、变率取决于超过了多少应力。上而一维的率相关塑性公式可以很方便地推广到三维情况。对于小应变的情形力度量之间无须区分，这里采用Cauchy应力6塑性率参数由应力和内变量的经:数给出。对照一维情况，三维情况下分解应变率为弹性和塑，性部分：(7=了).WF)(19s=ee+(19应力率和弹性应变率之间的关系为F=0,则=/;如果/0,则(f)二0击为Mises等:力,云为等效应变，神为率敏感系数。对于Johnson-Cook模型，可以得到等效塑性应变率的表达式为(19ao=A+B(19ao=A+B(19Eee=E(Eep)(ao=A+B(19ao=A+B(19Eee=E(Eep)(其中为静态屈眼应力:

8、Eee=E(Eep)(TOC o 1-5 h zsn+l二兔+壮(1屠+1=叮十為+1G+1(1盛+Q=骑+九+1.九十1(1%+1二C：(畸+1礁十1)(1fni-/(rn+in+i)=0(1在时刻肚给出一组（馮，?,qn）和应变增量爲公式（19-17）是一组关（為+“瑞务I）的非线性代数方程。将公式（19-17b）代入（19-17d）得到為十1=C：（程n+lr畔+1）=C：务+醇一码j二C：（en一咄）+C：Ac口公十1=(r+C：Ac)-C：典+itrailj-6+1一二典+t=AA(0)=O,(f(o)=Cj略+护)应力在第k次迭代时为ff(A=rt-co)一AA(i)C护(励且在整

9、个算法的塑性修正状态过程中始终保持不变，因此塑性应变的更新是久的參函数。在k次迭代时将检查屈服条件：(19-：(19；严）二存的一aY（&）二（笄）-3#AA）却3）若收敛，则迭代完毕，增量步结束。否则将计算塑性参数的增量：加_育一卄；（）令（巨）并对塑性应变和内变量进一步更新:(19-2：(19-2；(19-2：应+3二Fa)-2阙叫号nFu+n=W+加AAU+1)=3+虻=ptk)+厶严)（19-2；（19-2：然后将更新的变量返回屈服条件进行检査，整个过程将重复直至收敛为止。这就是埠步中应力更新的过程。19.3ABAQUS用户材料子程序用户材料子程序（user-definedmateri

10、almechanicalbehavior简称UMAT）通迥ABAQUS主求解程序的接口实现与ABAQUS的数据交流。在输入文件中，使用关铤*USERMATERIAL表示定义用户材料属性口19.3.1子程序概况与接口UMAT子程序具有强大的功能使用UMAT子程序有以下几点需注意:1）可以定义材料的本构关系，使用ABAQUS材料库中没有包含的材料进行计扩充程序功能。2）几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋ABAQUS中的任何单元。3）必须在UMAT中提供材料本构模型的雅可比（Jacobian）矩阵，即应力增量应变增量的变化率。4）可以和用户子程序“USDFIJT联合使用

11、，通过USDFUT重新定义单元一物质点上传递到UMAT中场变量的数值。由于主程序与UMAT之间存在数据传递*甚至共用一些变量，因此必须遵守有TTMAT園爭宜構才TTMAT中誉田的恥畐肚方枠#乩丰刘审暫滴直施倉吗.3NDI,NSHR,WTENS,NSTATV,PROPS.NPROPS,COORDS,IWT.PNEDT,LENT,DFGRDO,HFGRDl,NOEL,NET,LAYER,KSPT,KSTEP,KING)INCLUDE1ABA_PARAM.INCCHARACTER*80CMNAMEDIMENSIONSTRESS(NTENS),STATEV(NSTATV),IXJSDDE(WTENS,

12、tTTEHS),MSDDT(1WNS),DRPLDE(NTENS)fSTRAN(MENS),DSTRAN(WIENS),TIME(2),PREDEF(1)TDPRED，PROPS(NPROPS),COORDS(3),DROT(3,3),LFGRDO(3.3),HFGRDl(3f3)usercoding七0defineDDSDDE,STRESS,STATEV,SSE,SPD、SCDand,ifnecessary,RPL,DDSDDT,DRHBE,imOT.PNEBJCTRTURHUMAT中的应力矩阵、应变矩阵以及矩阵DDSDDE、DDSDDT.DRPLDE是直接分量存储在前，剪切分量存储在后。直

13、接分量有个DNI个，剪切分量有个。各分量之间的顺序根据单元自由度的不同有一些差异，所以编写UMAT时到所使用单元的类别。下面对UMAT中用到的一些变量进行说明。DDSDDE(NTENS,NTENS):是一个NTEN5维的方阵，称作雅可比矩P尬是应为的增量，虚是应变的增量，DDSDDE(/,J)表示增量2时第J个应变分量的改变引起的第个应力分量的变化通常雅可比是一个对称除非在“诜USERMATERIAL”语句中加入了UNSYMM参数乜STRESS(NTENS)：应力张量矩阵，对应NDI个直接分量和NSHR个剪切在增量步的开始，应力张量矩阵中的数值通过UMAT和主程序之间的接口申UMAT中；在增量

14、步的结束，UMAT将对应力张量矩阵更新。对于包含刚体转Z限应变问题，一个增量步调用UMAT之前就已经对应力张量的进行了刚体转动,在UMAT中只需处理应力张量的共旋部分。UMAT中应力张量的度量为柯西(应力。STATEV(NSTATEV):用于存储状态变量的矩阵，在增量步开始时将数勺到UMAT中。也可在于程序USDFLD或UEXPAN中先更新数据，然后在增量戈时将更新后的数据传递到UMAT中。在增量步结束时必须更新状态变量矩阵匚据。和应力张量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的数寸以外，状态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运2NSTATEV：状态变量矩

15、阵的维数，尊于关键字rDEPVAR中定义的数STANTS”常数设定的值PROPS(NPROPS)：材料常数矩阵，矩阵中元素的数值对应于关键字u*US：MATERIAL下面的数据行。SSE.SPD.SCD：分别是义每一增量步的弹性应变能、塑性耗散和蠕变耗散。们对计算结果没有影响，仅仅作为能量输出。其他变量：STRAN(NTENS)：应变矩阵。DSTRAN(NTENS)：应变增量矩阵。DTIME：增量步的时间增量。NDI：直接应力分量的个数。NSHR；剪切应力分量的个数。NTENS：总应力分量的个数，NTENS=NDI+NSHRo使用UMAT时需要注意单元的沙漏控制刚度和横向剪切別度。通常减缩积分

16、単的沙漏控制刚度和板、壳、梁单元的横向剪切刚度是通过材料属性中的弹性性质定的。这些刚度基于材料初始剪切模量的值，通常在材料定义中通过ELASTIC，项定义。但是使用UMAT时，ABAQUS对程序输入文件进行预处理的时候得不到變模量的数值，所以这时候用户必须使用JHOURGLASSSTIFFNESS选项来定5!有沙漏模式的单元的沙漏控制刚度，使用*TRANSVERSESHEARSTIFFNESS项来定义板、壳、梁单元的模向剪切刚度。19.3.2编程基于上面所述的率相关材料公式和应力更新算法参照ABAQUS用户材料子稍的接口规范*进行UMAT的编程。有限元模拟结果将在下一节给岀，在最后一节中给出了

17、相应的程序源代码。由于UMAT在单元的积分点上调用，增量步开始时，主程序路径将通过UM的接口进入UMAX单元当前积分点必要变量的初始值将随之传递给UMAT的相応量。在UMAT结束时，变量的更新值将通过接口返回主程序。整个UMAT的流砖图19-2所示。一共有8个材料常数需要给定，并申请了一个13维的状态变量矩阵，它们表示物理含义如表19-1所示。19-1UMAT材料常数PROPS12345678物理性质杨氏模量泊松比塑性耗散比ABnCmSTATEV1671213变量意义弹性应变塑性应变等效鳖性应变512512$11ra19-2UMATJt稈囹19.4SHPB实验的有限元模拟F面将建立SHPB实验

18、的有限元模型.井把前面所建立的UMAT接入ABAQUS/Standard进行有PH元棋拟，进行有眼元模拟的目的不绘单地为了禅现SHPB实验的过程.同时也足为了对选择的ttffim和建立的IJMAT进行评价，19.4.1分离式Hopkhison压杆(SHPB)奕監分离式Hopkinson压杆(splitHopkirvonpressurebar.面称SHPB)实验是从经典blopkinson实脸系砒之上发展而来的一种实验技术.用来测供材料的动念应力应变行为。该实验技术的理论慕础是堆应力波理论.通过测傲两根压杆上的应变来推寻试件上的应力-应变关系。分离式Hopkizn压杆实验的示jg图见图193。I

19、U19-35/：HcvkinwimiF*W19.4.2有限元連模有限元模熨主要杲倉照前血介绍的SHPB实验装程.if过載荷、边界条件爭的定义在有限元中棋拟SHPB实脸的歼境.尽萤在较少的机时耗费F达到更髙的榊度.L援型的尙亿与有限无岡站为f不使棋熨过干庞大对箱型进行了一些简化。首先.改变入力杆和出力杆的尺寸.长度由K(来的MMUrwn减小为1(MWmm,直位堆加到25nun,试件的长度和点径分别变化为22mmfU18mmv这样不仅优化了网格的.还成倍地咸小了復中的規模其带来的负面影响就是试件能达到的应变将降低另外.由于撞击杆仅仅起到产生应力脉神的作用庄数值模型中没必姿有虑曲击杆.取代的方法足点

20、接在入力杆的输入端魏加均布的应力關冲著虑到实够裝世的对称性.也作了一觀简化。菸个实验装匱以及戟荷弄柿是关于杆的中心线辅对称的所以可以使用轴对称血元进行二维分析。刃外也建辽了四分之一横嚴面的维模型作为补充一蠢轴对休模型和三维棋熨分别如图194、图195所示。在棋乜中对试件以及入力杆、出力杼和试件接粧的部分进行了局85阿格加密.这样的网格IM分可以取得比校经济的结果。團194元禅51fflN5堆有阪元锲型$11单元类型上，选择一阶常规单元，由于没有使用减缩积分单元，所以使用时无须指定单元的沙漏控制刚度。最后的模型中，二维网格单元的总数为122C模型网格的单元总数为17160.关于单元的详细信息如表

21、19-2所示。表1头2模型信息-型模L”X-单-一数个元单总二籐型杆力入-一Z-I一5-n41-1-31杆力出三養型一杆力入r勺r-5007-一711一-杆力出-750在二维模型局部网格存在疏密连接的部位，一个单元边要同时和两个单元逍在通常的有限元网格中是不易实现的。本例在这里使用了ABAQU5的多点约(multi-pointsconstrain,简称MFC)来解决，线性多点约束方式如图19-6所示口图中点p使用线性多点约束后，其节点自由度均由旁边的节点抚和方线性到。所以使用多点约束方式可以很好地连接模型中网格疏密不同的部位，划分出细的网格。2材料定义入力杆和岀力杆使用线弹性材料，弹性模量和泊

22、松比分别为200GPa和0二为7,85xio3试件采用用户在UMAT中的自定义材料*材料参数如表表193试件的材料定义性质密度/(kg/m5)杨氏模量/MPa泊松比Johnson-Cook模型参数A/MPaB/MPanCm数值27*068.0X1030.3366x562108.8530.2380.0290J3边界条件为了保证SHPB实验的要求，在二维模型和三维模型中均施加了必要的边界条在对称轴或对称面上施加了对称性边界条件，同时保证压杆和试件可以沿轴线方向E无约束地运动。压杆和试件之间的接触为硬接触，光滑无摩擦*为了确定输入应力脉冲的时间，进行了简单的计算.弹性材料中纵波波速的计爭式为$11其

23、中E为材料弹性模量，卩为材料密度。由此可以计算输入应力波在压杆中的传*度为Cd=5048m/so要求在入力杆应力波的输入端不能出现入射波和反射波的重叠，也就是说在输丿力脉冲的时间内，应力波的传播距离不应超过两倍的杆长，即TsAtJ4=TEtef*/Cr.上茲丄urznixL-tAstiz上八宀ZnI+rfIAtJ4=TEtef*/Cr.上茲丄ur$11图1916四种应变率下的应变-时间曲线70伽顶ss应应应应二二J-J:-80604020-04mm0.00000.0001O.OOO20,00030.000400005图19J7四种应变率下的应力时间曲线生，引起单元温度的变化，但是相邻单元之间没

24、有热量的传递。在应变率250s下件最终的温度场如图19-21所示。取塑性耗散比为0.95,即95%的塑性耗散将转化为热量。这一比例对于大多数属是合适的.从图19.21中看出，早期应力波没有到达试件时.试件上的温度没有:$11T应变率100应变寧250图19-18四种应变率下的应力-应变曲线Edw、R圏图19-19四种设计应变率下的真实应变率-时间曲线将考虑温度效应时单元的应力应变与前面没有考虑温度效应的情况作一下比件表面中间单元的温度曲线如图1922所示。曲线中发生塑性变形后，单元的温上升。右祖詹隔毘响卞单三萌庙*审亦兩由站纬伽園IQg霸条肚闫酚苗木窘$11z1a*锲sJu魚篆：SWM一z-6

25、-s#但qsY2xofss*0仝君?-*辺&Hs-E:*tiqe(设I8-$11sVIS*整連驴裂ST尝1VM355SE写猱J国$&班七SH徉另Bgm亠89U?_EOWW鬼壬誉整房艺糊廻注艺ww9sfe潑生E4姿題如犖aJYlte5-$1128$11$11270耳000000.0001o6oo200003o.6oo4o.doos时间虫图19-22单元温度的变化$11$110.04-003-制圏0020.01aoo0,0000o.oooi$11（町应变时闾曲钱$11$11时同沧800800524urr1-r-t1*r-OQOO0850.010001$n.(i200.02500300035(XH0

26、OXHS800800524800800524图19-24专球催度效应与不(Sfi度效应的比较19.4.4三维动态分析本例IS进行的二维动态分析也相当成功试件往康拟过程中各处应力都比较均匀議大应力和廉小应力的建别不超过IMPao图19-25为试件的Mises应力分布。800800524800800524800800524变形过穆中单元的应力、应变以及应变來历史桶与二维的悄况非常吻合如图1926所示。图1927所示为应力应变曲线.这是对编写的UMAT用于三维实体单元的-个笑证。UMATM然見从-tftJohn-son-Cook施熨中建立起来的但是在UMAT中率相关也性公式将它扩展列了三维情况.所以

27、能用于三维实体单元也在意料之中。虽怡取得的结果是一致的.但是在计算时间上却冇很大差别在个人计算机上使用轴对称-.维网格完成一次计篦只需要10min.然后完成-次潍网恪的计韓需翌3h。图1926三维有限元模拟的应力、应变率历史图19-27三维有限元模拟中的应力-应变曲线相比而言使用轴对称二维模型要经济得多。195UMAT的Fortran程序19.5.1UMATSUBROUTINEUMAT(STRESS,STATEV,DDSME,SSE,SEO,9CD,RPL,EDSDOT,DRPLEffi.EREWT,STRAN,DSTFiAN,524INCLUDE7ABAPARAM-INCCHARACTERB

28、OMATERLDIMENSIONSTBESS(WIE阳)，STAIEV旳/TV.DDSXe(NTENS,!nENS),DDSDET(inENS,DRElDE(iniMS人STRAMdnENShDSTRMinENShTIMEdhPFEWdhEREDd),PK5(NPROPS.COORDS.lWr(32),DFDO(3,3.DFGRD1(33)DIMEWSIONLAS，EELAS(6),F1QW524524PARAMETER(ONE=l.OIX),TTO=2,0D0tTHREE=3.ODO,SIX=+ODO,HALF=0.5d0)DATANEKTON,TOLER/40覧,96/CCUMATFORJ

29、OHNSON-COCKMODELCCPROPSYANGSMCDULUSCPFCPSPOISSONRATIOCFROPS(3)INELASTICHEATFRACTIONPARAMETOFJOHNSOJUDCKMCQEL：524524CPROPSATOC o 1-5 h zCPP0PS(5)七BCPROPS-nCFR0PS7)-CCPROPS(0)mIF(NDLNE.3)THENWRTrE(6,l)FQEMAT(/r3QXSl粉ERROR-THISUMATMAYCMLYHEUSEDFOR1ELEMEKTSWTIHTHREEDIRECTSTRESSCCMPONEJnS，)ENDIFCCELASTIC

30、PROPERTIESCEtW=FRQPS(l)ENUPE?OPS(2)EG3=THREE*EGELAM二(EBUTO-EG2)/THREECCELASTICSTIFFNESSCDO20Kl=l,NTEKSDO10K2=1,NTENSDDSM(K23KI)=0R10COHTINUE20OGMTINUECDO40Kl=l.NDIDO30K2=1?NDIDDSDDE(K2,K1)二ELAM30COKTINUEEDSCCEtKl,KI)=+ELAM40CONTINUEDO50KI=NDI+l.NTENS50CONTINUECCCALCULMESTRESSFROMELASTICSTRAINSCDO70K

31、l=lsNTENSDO60K2=1,NTEHSSTRESS(K2)=STRESS(K2)+DDSEC(K2,Kl)并DSIRAN(Kl)60CONTINUE70CONTINUECCREOOVERELASTICANDPLASTICSTRAINSCDO80Kll.NTENSLRS(Kl)-STATEV(Kl)+DSTRfiN(Kl)EPLAS(Kl)=STATEV(K1+WIENS)80(XMINUEEQPLRS=STATEVC1+2供NTENS)CCALCULATEMISESSTRESSDO90KI=NDI+l.NTENSSHISES=SMISES+SIX*STRESS(K1)*STRESS(K

32、1)90OCSTINUESHISES=SQRTfSMlSES/TWO)CCCALLUSEPHAK3SUBROUTINE,GETHARDENINGRATEANDYIELDSTRESSCCCALLUSERHMD(SYIELO,HMD,EQPLAS,H?0PS)CDEIERHINEIFACTIVELYYIELDINGCIF(SMISES.GT_(1.0+TOLER)*SYIELO)THENcCMATERIALRESPONSEISPLASTIC,DETERMINEFLOWDIRECTIONCSHDRO=(STRESS(l)+STRESS(2)+STRESS(3)/THREEONESY=ONE/SMIS

33、ESDO110KI=1,NDIFTjCW(KI)=ONESY*(STRESS(K1)-SHYDRO)110CONTOJUEDO120KI-NDI+1,NTENSFLOW(Kl)=STRESS(K1)*ONESY120CONTINUElaREADPARAMETERSCFJOHNS(潞COOKFCA=TOPS(4)B=PROPS(5)EN=PR0PS(6)C=PR0PS(7)EM=PR0PS(6)CCNEWTONITERATIONCSyiELD=SYIELODECPL=(SMISES-SYIEED)/ES3DSTRES=IOLER*SYIELO/EG3DEQMIN=HALF*nriME*EXP(1

34、.OD-4/C)DO130KEWTON=1FNEWTONTTfTtviriThlTtvirnrnrETtr524HARD1=HM3TVRL+SYIELX/EEQPLSYIELD=SYIELD*TVP1TTfTtviriThlTtvirnrnrETtr524TTfTtviriThlTtvirnrnrETtr524RHS-SMISES一EG3#FL-SYIELDTTfTtviriThlTtvirnrnrETtr524EEQPL二DEQRL+RHS/(EG3+HARD1)IF(ABS(RHS/EG3).LE.DSTRES)GOTO140130CONTINUEWRITE(6,2)NEWTON2FORM

35、AT/,30X,*WARNINGPLASTICITYALGORITHMDIDNOT1CONVERTAFTER13,ITERATICSS)140CONTINUEEFFHRD=EG3*HARD1/(EG3+HARD1)CCCALCULATESTRESSANDUPDATESTRAINSDO150KI=1#NDISTRESS(Kl)=FLOW(K1)*SYIELD+SHYDROEPLfiS(Kl)=EPLRS(K1)+IHREE餐FLO锹K1”DEQPL/TO150CONTINUEDO160K1=NDI+1,NTENSSTRESS(Kl)=FLOW(SYIELO+SYIELD)/IWORPL=PROPS星SPD/DTIMECCJACOBIANEFFG=EOSYIELD/31ISESEETG2=TH(E