初识ABAQUS UMAT二次开发（二）
大家好，这一期小编接着上一期（初识ABAQUS UMAT二次开发）继续向大家介绍ABAQUS软件的UMAT二次开发功能。上一期我们说过ABAQUS是一款强大的非线性有限元软件。使用UMAT一般是为了实现ABAQUS材料库中不包含的本构模型。在这一期本来小编想介绍如何使用UMAT实现弹塑性本构关系，可是后来再百度中发现这些东西已经烂大街了。小编索性就自己臆想了一个本构关系，即应力-应变关系曲线符合二次抛物线函数。这样既有了上升段，也有下降段。通俗的讲就是可以实现随着荷载递增结构内力先增大，其后降低直至失去承载能力。
下面小编就带大家一起看看这种本构关系是如何在ABAQUS中用UMAT实现的。
模型简介
为了让大家更好的理解UMAT的使用方法。这一期的模型跟上一期一样采用最简单的桁架结构，如图1所示。模型由三根杆件组成，一根竖杆、两根斜杆。竖杆长为30，与斜杆的夹角均为45度，截面面积为1。下面三个节点均施加x、y向约束，上部节点施加一个向上的节点位移荷载0.0017。单元仍旧采用Truss单元，每根杆件划分为1个单元。
图1 有限元模型及边界条件示意图
材料属性设置：在User Material中定义两个参数，一个值是-7244290000，一个值是412000，这两个数分别是应变-应变二次抛物线函数的两个系数a和b。通过变量props传递给UAMT子程序。应力-应变关系曲线如图2所示：
图2 本构关系曲线（应力-应变关系曲线）
分析步和求解器设置：小编在试算的时候发现，如果将荷载定义为节点力，那么使用abaqus/standard Static General程序求解收敛性不稳定，使用Dynamic/implicit的拟静态程序求解收敛性会好很多。本例施加的是节点位移荷载，我们选用Static General即可。在增量设置里选用Fixed类型，最大增量步设置为100步，增量步长设置为0.01。求解器类型选用迭代法，其它保持默认即可。
有关abaqus/standard求解非线性问题的相关介绍或知识，大家可以查看abaqus的帮助文档或者庄茁老师的《基于abaqus的有限元分析和应用》的第八章，介绍的非常详细。
UMAT子程序设置
本例用到的UMAT相关知识与上一期基本相同。唯一不同的是多使用了一个变量STRAN。STRAN（NTENS）：应变矩阵（或向量），包括NDI个直接分量，NSHR个剪切分量，与应力变量STRESS相对应。整个UAMT子程序如下：
1 SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
2 1 RPL,DDSDDT,DRPLDE,DRPLDT,
3 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
4 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
5 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
6C
7 INCLUDE 'ABA_PARAM.INC'
8C
9 CHARACTER*80 CMNAME
10 DIMENSION STRESS(NTENS),STATEV(NSTATV),
11 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
12 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
13 3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
14 4 JSTEP(4)
15C
16 REAL*8 emod
17c ----------------------------------------------------------------
18c umat for nonlinear problems
19c be used for Truss element
20c ----------------------------------------------------------------
21c props(1) - a
22c props(2) - b
23c ----------------------------------------------------------------
24c
25 if (ndi.ne.1) then
26 write (7,*) 'This umat may only be used for elements with one direct stress components'
27 call xit
28 endif
29c
30c elastic properties
31 emod = 2*props(1)*STRAN(1)+props(2)
32c
33c elastic stiffness
34c
35 ddsdde(1, 1)=emod
36c
37c calculate stress
38c
39 stress(1) = stress(1)+ ddsdde(1,1)* dstran(1)
40c
41 RETURN
42 END
代码1~15行是UMAT的固定代码。
第16行定义了一个实数变量emod用于存储剪切刚度。
第25~28行，用于单元类型判断。当单元节点自由度数不等于1时，报错。
第31行计算剪切刚度，即应力-应变关系函数的一阶导函数。不清楚的可以看一下abaqus帮助文档中关于非线性求解器的介绍，或者百度一下牛顿迭代法的相关内容。
第35~42行，与上一期相同，都是用于根据应变增量更新应力。
计算结果
荷载作用下结构各单元的内力变化示意图：
图3 结构的内力变化示意图
由图3可知，竖杆的内力变化明显经历了两个过程，从随荷载的增大而增大到随荷载的增大而降低。竖杆的应力随荷载增量的变化曲线，如下图所示。
图4 竖杆的应力随荷载增量的变化曲线
可以从abaqus计算结果中导出竖杆的应力-应变关系曲线，如下图所示。
图5 从abaqus计算结果中导出的竖杆应力-应变曲线
对比图5和图2可知，模型中竖杆的应变-应变曲线与输入的本构关系完全一致，可见UMAT程序正确。
（本示例UMAT源码、CAE模型均可在公众号QQ群文件中免费获得）
结语
Abaqus UMAT用户子程序接口功能非常强大，编写一个漂亮的UMAT子程序需要一定的力学基础（尤其是有限元、连续介质力学）。本文旨在让大家对UMAT有一个初步的了解，以及使用UMAT需要用到哪些基本知识，好在今后的学习中有侧重点。后续的推文中，小编还会进一步结合我们有限元理论模块的进度，进一步讲解Abaqus UMAT的使用。
最后欢迎批评指正，小编水平有限，希望大家跟我一起完善这部教材，造福更多的人，谢谢！
陪伴是最长情的告白
学而不思则罔
思而不学则殆
欢迎大家跟我一起学习CAE知识！
QQ群号：517718332
注：相关资料可在QQ群文件中下载
微信
QQ
微博