【JY】 ABAQUS子程序UEL的有限元原理与应用

1年前浏览9209

【简述ABAQUS中UEL子程序】

ABAQUS作为成熟的商用有限元软件，可为高级用户提供特定的分析需求。ABAQUS常见的二次开发子程序包括：UMAT、VUMAT、UGENS、UEL和VUEL等。其中UEL/VUEL分别适用于ABAQUS的Standard/Explicit求解器。只有清楚有限元分析的基本原理，才能够较好地了解其分析的力学原理，才能对特定的分析需求编写合适的分析单元。

本文对ABAQUS/Standard中UEL子程序的有限元原理及编写原则进行初步梳理，并编写了平面三角形单元静力分析的UEL子程序，以加深对UEL的理解与认识。

【ABAQUS的UEL单元之有限元基本原理简述】

当用户需要用ABAQUS/Standard单元库中没有的单元进行分析时，可通过ABAQUS提供的UEL子程序接口进行二次开发，编写适用于特定分析的单元。下面以平面三角形单元为例，给出编写UEL单元的基本流程。

系统的势能Π由下式给出

其中为单元的应变能。

将应变-位移关系ε = Bq代入U_e中，由此对于势能方法形成的刚度矩阵如下：

其中t_e为单元厚度、A_e为单元面积。

有限元分析的第一步是将实体离散化为多个单元，此后构造出各个单元的单元刚度矩阵，在将单元刚度矩阵集成整体的刚度矩阵，最后通过整体刚度矩阵建立平衡方程从而求解各个节点的位移、应力、应变等响应。因此，根据有限元的分析原理，编写UEL的最终目的即是形成目标单元的单元刚度矩阵。以下对平面常应变三角形单元（CST）进行有限元分析，依据有限元分析思路编写UEL子程序，以初步了解有限单元的分析原理及UEL编写流程。

图2 平面三角形单元示意图

上图所示为平面三角形单元，由图示信息可得其节点坐标列阵u与位移列阵q分别为

1.形函数

设三角形单元节点1、2、3对应的形函数分别为N1、N2、N3，形函数满足条件：N1+N2+N3=1（在点i处，形函数Ni值为1，其余形函数值为0），形函数采用自然坐标ξ和η描述，有

2.高斯积分点及权重

此平面三角形单元所选取的高斯积分点及权重如下：

表1高斯积分点及权重

3.形函数与单元位移关系

得到形函数与单元节点位移后，根据等参元表示方法，单元内任意一点的位移都可用形状函数和未知结点位移场进行表示，有

可表示为u=Nq，其中N为

对于三角形单元，可将x、y坐标表示为节点坐标的形式，可由此得到单元坐标的插值关系为：

其中

4.雅各比矩阵J

由节点坐标x、y与ξ、η关系，在根据链式求偏导法则，可得到平面三角形单元的雅可比矩阵

5.应变矩阵B

根据应变ε与节点坐标u、节点位移q的关系可得应变计算公式为

其中xij与yij计算公式同上。写成矩阵形式为

由此可导出单元应变-位移矩阵B，其表达式为

由此可看出B矩阵中所有元素均为常数，且每个常数均是通过节点坐标表示的。

6.弹性矩阵D

弹性矩阵D根据所研究问题为平面应力或平面应变问题有不同形式，其取值需根据所用材料确定，具体如下：

7.应变-位移矩阵

由已有的应变矩阵B与应力应变关系，得应变位移矩阵σ = Dε=D B q。

8.单元刚度矩阵与残余力向量

根据势能方法导出得到三角形单元的单元刚度矩阵k^e，有

其中t^e为单元厚度。

最后，通过单元的边界条件可得单元的残余力向量。至此，UEL单元的编写完成，可进行初步调试直至结果无误。

【 UEL的编写基本原则】

ABAQUS的子程序可通过FORTRAN 77/95 进行编写，本文采用FORTRAN77进行编写。

在FORTRAN中编写UEL，需在以下框架中进行，才可被ABAQUS识别。以下框架也即是ABAQUS的UEL接口（见ABAQUS6.14 User Subroutines Reference Guide 1.1.28 UEL）。

SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
     1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,
     2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,
     3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),
     1 SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL),
     2 DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*),
     3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
     4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)

      user coding to define RHS, AMATRX, SVARS, ENERGY, and PNEWDT
（此部分是我们所要编写的内容）

      RETURN
      END

UEL的编写可分为以下几个部分内容：

1.定义变量

编写UEL子程序最重要的工作是形成单元的刚度矩阵AMATRX与残余力向量RHS。单元刚度矩阵AMATRX与残余力向量RHS的形成可按照上一小节给出的流程进行。UEL子程序框架中出现的如SVARS、ENERGY、COORDS等变量无需定义（UEL子程序框架中各变量的具体含义见帮助文档），可直接赋值。

2.变量赋值与信息传入

在开始采用所编写的UEL单元时，需要从inp文件中传入单元所用材料的基本属性，如密度、弹性模量及泊松比等。同时，在inp文件中需要指定UEL单元的单元结点数NODES、单元类型命名TYPE、基本属性个数PROPERTIES、空间维数COORDINATES、产生的变量个数VARIABLES、单元节点编号*ELEMENT等信息。

参考Abaqus6.14帮助文档中调用UEL子程序时inp文件所用的命令流：

*USER ELEMENT, NODES=2, TYPE=U1, PROPERTIES=4, COORDINATES=3,

VARIABLES=12

1, 2, 3

*ELEMENT, TYPE=U1

101, 101, 102

*ELGEN, ELSET=UTRUSS

101, 5

*UEL PROPERTY, ELSET=UTRUSS

0.002, 2.1E11, 0.3, 7200.

上述命令表示所采用的UEL单元为两点的线性杆件单元，其中ELGEN表示逐步增加单元数量（详见帮助文档关键词解释*ELGEN）。要将inp文件中的单元特性传入UEL，可在UEL中写入以下命令：

AREA= PROPS(1) ! Cross-sectional area