Abaqus 自定义单元 UEL 子程序实例研究(附源码)

2年前浏览4237

Abaqus 的UEL子程序的出现，使得大家可以自定义单元来满足解决实际科研、工程问题的需要，比如自定义粘结单元(cohesive elements)等。但另一方面，作为 Abaqus 最难的子程序，自定义单元UEL子程序需要较高的有限元力学基础、灵活的编程能力和持久的调试耐心。并且，网上涉及UEL子程序的教程较少，还主要局限在 Abaqus 在线帮助手册上对于UEL的解释内容，未能有相关结合具体问题的代码实例来对UEL子程序做一个研究学习。

下面结合一个经典的4节点平面应力单元的自定义UEL子程序来进行深入explore：

并且UEL的建模一般没有GUI（点点点式）这种方式，都是通过inp文件直接导入相关单元几何节点坐标信息，下列inp代码第5-9行。前3行为一般性的inp开头行。

代码段1：

*Heading** Job name: 4 nodes quad Plane Model name: UEL*Preprint, echo=YES, model=yes, history=yes, contact=NO
*Node      1,-1.,-1.      2,1.,-1.       3,1.,1.       4,-1.,1.
*User element, nodes=4, type=U1001, properties=3, coordinates=2, variables=71,2*Element, type=U1001 1,1,2,3,4*Elset, elset=El 1*Uel property, elset=El  1.0, 2.0E11,  0.3*Nset, nset=Node, generate  1, 4, 1*Nset, nset=N1  1, 2*Nset, nset=N3  3,4
*Step, name=Force*Static0.01, 1., 1e-05, 0.01*BoundaryN1, 2, 2N1, 1, 1*CloadN3, 2, -1e 7*Restart, write, frequency=0*Output, field*Node OutputCF, RF, U*Element Output, directions=YESLE, PE, PEEQ, PEMAG, PS, S, VE, VEEQ, VS,SDV*Output, history, variable=PRESELECT
*EL PRINTSDV1,SDV2,SDV3,SDV4,SDV5,SDV6,SDV7,SDV8,SDV9
*End Step

代码段1的第5-9行表示用inp输入的节点信息，即1号节点的坐标(-1,-1)，2号节点的坐标(1,-1)，3号节点的坐标(1,1)，4号节点坐标为(-1,1)。

UEL子程序接口

代码段1第11行到第18行，为UEL子程序特殊的声明段，其接口定义如图1所示。

图1 UEL Interface定义形式

我们先学习 Abaqus 关于UEL的官方解释内容（见图2）：*USER ELEMENT这个选项必须出现来定义UEL子程序。*USER ELEMENT后面的TYPE为Un，Un为单元类型，必须以U开头，可以为U1,U2等；NODES为每个单元具有的节点数量；COORDINATES为每个节点具有的坐标数；PROPERTIES为单元浮点数材料属性的参数个数，比如弹性模量、泊松比等；I PROPERTIES为单元整数材料属性的参数个数；VARIABLES为用户自定义结果变量的数目（Number of SDVs）；UNSYMM决定雅克比矩阵是否为非对称矩阵。

因此，对于代码段1中的*User element, nodes=4, type=U1001, properties=3, coordinates=2, variables=7，则表示一个单元有4个节点，单元类型定义为U1001，单元的浮点数材料属性参数有3个，每个节点有两个坐标，有7个结果状态变量（SDV）。

图2 UEL各变量解释

如图3，*USER ELEMENT下的Data line(s)这里输入激活的节点自由度，比如dof_1表示节点的第一个激活自由度，dof_2表示节点的第二个激活自由度。需注意，如果所有自定义单元拥有相同的激活自由度，则不需要进行进一步的数据输入。对于输入的自由度数据，1-3表示位移，4-6表示旋转。

而对于代码段1中*User element下的第12行1,2则表示每个节点有x,y位移两个激活的自由度。

图3 UEL中的Data Lines

代码段1的第13-14行表示建立单元类型为U1001的1号单元，该单元由1,2,3,4四个节点组成；代码第15行表示将1号单元建立单元集合El（element set）；而*UEL PROPERTY表示赋予自定义单元材料属性，如图4，通过*UEL PROPERTY该方式，即便是同一UEL子程序也可以赋予自定义单元不同属性，并且在*UEL PROPERTY中定义好了相关属性，可不必要在UEL子程序中定义单元材料属性。因此代码段1的第17到18行表示赋予自定义单元集合的材料属性为1.0，2.0E11和0.3。

图4 *UEL PROPERTY的注解

代码段1的第19-20行表示通过generate(除了Set名，inp文件一般不区分大小写)来生成节点集合，具体解释如图5。因此第19-20行表示从1号节点开始，4号节点结束，每次递增1；这些生成的所有节点放在Node这个Set里。而代码段1的第21-24行跟前面建立单元集合的方式一致，即节点编号1,2放在N1节点集(NodeSet)合里，节点编号3,4放在N3节点集合(NodeSet)里。

图5 INP文件中关于generate的注解

这里还有个不好理解的VARIABLES？其实如图6，VARIABLES表示每个单元的SDVs（结果状态变量数目）。其中，如果单元采用数值积分，VARIABLES应该等于单元积分点数目乘以每点的SDVs数目。可采用代码段1的第42-43行来输出SDVs。UEL子程序的SDVs（SDV1,SDV2）只能输出到data(.dat)，results(.fil)，output database(.odb)文件，也可以用于X-Y plots绘图。

图6 VARIABLES的注解

Monitor界面即可查看data(.dat)文件，如图7，可看到其中输出的SDVs的数据。

图7 查看data(.dat)文件中的结果

建立分析步

代码段1的第26行表明建议名为Force的一个Step，Step的类型为静力分析（Static），结合图8，则代码段1的第33行表示初始时间增量为0.01，时间总步长为1.0，最小时间步长为1e-05，最大时间步长为0.01。

图7 *STATIC下时间步长注释

添加荷载边界条件

结合图8，*BOUNDARY后分别是节点编号或者节点集合，第一个自由度，最后一个自由度，对应自由度的大小。代码段1的第30、31行表示约束住N1节点集合的1、2自由度，即x和y方向的平动自由度。

图8 *BOUNDARY关键字注解

结合图9，*Cload后分别是节点编号或节点集合名称，自由度，荷载大小。代码段1的第33行表示对N3节点集合沿着2自由度方向（Y方向）施加-1e 7N的集中力。

图9 *Cload关键字注解

代码段1的第34行表示重启动分析，如图10。*restart重启动分析表示对之前的分析进行再次读写，但基本模型数据是不能改动的，比如单元、材料和节点，能够改动的只能是单元集合（单元set）和节点集合（节点set）。这里第34行表示重新write，frequency=0表示打断重分析数据写入的连续性，如图11。