1. 简介

在Abaqus/Standard模块中，用户可以利用子程序UEL来开发自定义单元，用以实现一些通过Abaqus内置单元无法实现的功能。如果编写恰当，用户自定义单元可以正常使用Abaqus/Standard的大部分功能，并且在用户子程序UEL中可以编写多个自定义单元，这些单元可以同时使用。

相比于开发完整的有限元程序，使用Abaqus来开发单元的优势是非常明显的：首先，Abaqus提供了大量的结构单元、分析程序以及前处理工具，这些功能可以与Abaqus UEL配合使用；其次，Abaqus中可以完成自定义单元的部分后处理工作；并且，开发子程序的效率也远高于开发完整的有限元程序。

尽管用户子程序UEL为在Abaqus中进行有限元分析提供了非常强大和灵活的工具，但由于开发自定义单元需要较高的数学和力学理论基础，加之前、后处理过程不便，因此相比于诸如UMAT之类的子程序，可供参考的资料非常少。恰好笔者在完成毕设的过程中需要使用到用户子程序UEL进行结构单元的开发，因此借由此机会和大家分享一些子程序UEL的开发经验以及存在的问题，希望可以帮助到有这方面需求的同学。

为了理解用户子程序UEL的工作流程，必须熟悉有限元分析求解的整个过程。因此，本文以最为简单的杆单元为例，介绍杆结构的有限元分析流程，随后利用子程序UEL开发自定义杆单元，并验证计算结果的准确性。

(为了便于排版，部分内容采用图片格式)

附录1：自定杆单元UEL用户子程序

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC     程序用于定义2节点杆单元                                              CC     单元节点数:2   节点自由度数:1     单元自由度数:2×1=2                  CC     材料参数: (1)弹性模量E (2)截面积AREA                                 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      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(*)CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      C                       子程序主要变量说明                                 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC     需要更新的变量CC     AMATRX:单元的刚度矩阵 所有非零元素都必须定义C     RHS:单元刚度方程中的右端向量 变量中包含有残余载荷(不平衡力)向量C     SVARS:求解状态变量 状态变量的个数由NSVARS确定(可以不更新)C     ENERGY:用户自定义单元的能量 共有8个分量(可以不更新)C     传入模型的信息变量(不可修改)C     C     NNODE:自定义单元的节点个数C     JTYPE:单元类型 UnC     NDOFEL:自定义单元的自由度个数C     JELEM:用户指定的自定义单元号C     NSVARS:用户自定义状态变量的个数C     PROPS:单元材料参数实数数组 包含有NPROPS个实数参数C     JPROPS:单元材料参数整数数组 包含有NJPROPS个整数参数C     COORDS:坐标数组 COORDS(K1,K2)为第K2个节点的第K1个坐标C     U:单元计算中的自由度(本单元中为位移)C     DU:位移的增量值CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC                         变量定义及声明                                  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC          REAL::L              !杆单元长度      REAL::E              !材料弹性模量      REAL::AREA           !杆单元截面积      REAL::EPS            !单元应变      REAL::SIGMA          !单元应力      INTEGER::I,J         !循环索引      CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC                     定义刚度矩阵及右端矢量                               CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC     获取材料参数      E=PROPS(1)      AREA=PROPS(2)C     通过单元节点坐标计算杆单元长度C     节点坐标COORDS(K1,K2)表示当前单元第K2个节点第K1个坐标      L=ABS(COORDS(1,2)-COORDS(1,1))  C     构造刚度矩阵AMATRX      AMATRX(1,1)=E*AREA/L      AMATRX(1,2)=-1*E*AREA/L      AMATRX(2,1)=-1*E*AREA/L      AMATRX(2,2)=E*AREA/LC     构造右端矢量RHSC     RHS=-K・U      DO I=1,NDOFEL            !NDOFEL为单元的自由度数  本单元中为2          RHS(I,1)=0           !初始化右端矢量                DO J=1,NDOFEL              RHS(I,1)=RHS(I,1)-AMATRX(I,J)*U(J)    !按照矩阵相乘的公式计算RHS          END DO      END DO            CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC                        计算单元应力和应变                                CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC     计算单元应变      EPS=(U(2)-U(1))/LC     计算单元应力      SIGMA=E*EPSC     将计算结果储存置状态变量      SVARS(1)=EPS      SVARS(2)=SIGMA      RETURN      END

附录2：UEL用户子程序简介