不等待
即关注
【简述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代入Ue中，由此对于势能方法形成的刚度矩阵如下：
其中te为单元厚度、Ae为单元面积。
有限元分析的第一步是将实体离散化为多个单元，此后构造出各个单元的单元刚度矩阵，在将单元刚度矩阵集成整体的刚度矩阵，最后通过整体刚度矩阵建立平衡方程从而求解各个节点的位移、应力、应变等响应。因此，根据有限元的分析原理，编写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.单元刚度矩阵与残余力向量
根据势能方法导出得到三角形单元的单元刚度矩阵ke，有
其中te为单元厚度。
最后，通过单元的边界条件可得单元的残余力向量。至此，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
ENDUEL的编写可分为以下几个部分内容：
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
E = PROPS(2) ! Young's modulus
MIU = PROPS(3) ! Poisson's ratio
RHO = PROPS(4) ! Density
3.变量计算
①根据前述有限元基本原理，计算各个变量，并形成刚度矩阵
②在编写UEL过程中，可通过输出每一步的计算变量，以方便后期debug判断错误所在的位置，每个输出变量可在.log文件中查看。采用FORTRAN77 输出变量的方式：
（假设变量为雅各比矩阵Jacobi）PRINT *, 'Jacobi=' ,Jacobi
4.编写完成，开始debug
编写完UEL子程序之后，先将UEL子程序应用于单个单元进行测试（从静力、线性摄动、动力时程等各个方面逐一测试），在保证单个单元的计算结果正确无误的前提下，再将UEL子程序应用于较多单元的分析案例中，可保证分析结果的准确性。
【在ABAQUS中的UEL单元研究】
根据以上平面三角形单元的有限元分析思路，编写对应的平面二维三角形单元UEL子程序，并通过两个不同的有限元分析算例验证该子程序的有效性。算例中采用的三角形单元为ABAQUS单元库中的二维三角形单元CPS3。
算例一：平面三角形单元
一、三角形单元有限元计算参数
图3 平面三角形单元示意图
三角形单元的相关计算参数如表2所示：
表2 三角形单元计算参数
该三角形单元共有三节点（如图1所示），每个节点有2个自由度
二、边界条件及载荷设置
将节点1与节点3的位移自由度施加约束，在节点2的水平x方向施加100N的节点力。
三、计算结果对比
在相同边界条件、载荷条件下，二维三角形单元的UEL子程序计算结果与有限元ABAQUS计算结果如下表3及图4-5所示，从图中可看出，子程序的位移计算结果与abaqus中CPS3单元位移计算结果一致。
表3 UEL子程序与abaqus有限元最大位移对比
(a)CPS3加载时程
(b)UEL加载时程
图4两种不同单元的节点位移时程曲线
图5 ABAQUS杆单元计算结果与UEL单元计算结果对比
算例二：平面三角形单元组合
一、三角形单元有限元计算参数
图6平面三角形结构
由上图6可知该结构为四个相同的二维三角形单元组成。该杆件结构的相关参数如下表4所示：
表4 三角形单元参数
该平面三角形共有六个节点（如图4所示），共由三个平面二维三角形单元组成，每个节点有2个自由度。
二、边界条件设置
在三角形节点3施加100N水平力、节点4施加20N的水平力，并约束节点1、节点2的所有自由度，其余节点均为自由边界。
三、计算结果对比
在相同边界条件、相同载荷条件下，所编写的子程序位移计算结果与abaqus中的CPS3单元位移计算结果如下表4及图5-6所示，可从下图表得知，二者计算结果一致。
表5 UEL子程序与abaqus有限元最大位移对比
(a)CPS3加载时程-U1
(b)UEL加载时程-U1
图7-1两类单元节点的U1位移时程曲线
(a)CPS3加载时程-U2
(b)UEL加载时程-U2
图7-2两类单元节点的U2位移时程曲线
图8 CPS3单元计算结果与UEL单元计算结果对比
通过以上两个有限元算例可知，采用所编写的UEL二维三角形单元子程序所计算得到的位移及时程结果与有限元软件ABAQUS中采用CPS3单元对应计算结果一致，说明该三角形单元子程序可靠且合理。
【参考文献】
[1]：Abaqus6.14 User Subroutines Reference Guide
[2]: 工程中的有限元方法（第四版）
概念为先，机理为本，期待下篇！
往期推荐 ·
#性能分析
【JY】基于性能的抗震设计浅析（一）
【JY】基于性能的抗震设计浅析（二）
【JY】浅析消能附加阻尼比
【JY】近断层结构设计策略分析与讨论
【JY】浅析各动力求解算法及其算法数值阻尼(人工阻尼)
理念
【JY|体系】结构概念设计之(结构体系概念)
【JY|理念】结构概念设计之(设计理念进展)
【JY】有限单元分析的常见问题及单元选择
【JY】结构动力学之显隐式
【JY】浅谈结构设计
【JY】浅谈混凝土损伤模型及Abaqus中CDP的应用
【JY】浅谈混凝土结构/构件性能试验指标概念（一）
【JY】浅谈混凝土结构/构件性能试验指标概念（二）
【JY】橡胶系支座/摩擦系支座全面解析
#概念机理
【JY】基于Ramberg-Osgood本构模型的双线性计算分析
【JY】结构动力学初步-单质点结构的瞬态动力学分析
【JY】从一根悬臂梁说起
【JY】反应谱的详解与介绍
【JY】结构瑞利阻尼与经济订货模型
【JY】主成分分析与振型分解
【JY】浅谈结构多点激励之概念机理（上）
【JY】浅谈结构多点激励之分析方法（下）
【JY】板壳单元的分析详解
【JY】橡胶支座的简述和其力学性能计算
【JY】振型求解之子空间迭代
【JY】橡胶支座精细化模拟与有限元分析注意要点
【JY】推开土木工程振型求解之兰索斯法(Lanczos法)的大门
【JY】基于OpenSees和Sap2000静力动力计算案例分析
【JY】建筑结构施加地震波的方法与理论机理
【JY】力荐佳作《结构地震分析编程与应用》
#软件讨论
【JY】复合材料分析利器—内聚力单元
【JY】SDOF计算教学软件开发应用分享
【JY】Abaqus案例—天然橡胶隔震支座竖(轴)向力学性能
【JY】Abaqus6.14-4如何关联fortran？
【JY】如何利用python来编写GUI？
【JY】如何解决MATLAB GUI编程软件移植运行问题？
【JY】浅谈结构分析与设计软件
【JY|STR】求解器之三维结构振型分析
【JY】SignalData软件开发应用分享
【JY】基于Matlab的双线性滞回代码编写教程
【JY】动力学利器 —— JYdyn函数包分享与体验
【JY】混凝土分析工具箱：CDP模型插件与滞回曲线数据
【JY】结构工程分析软件讨论（上）
【JY】结构工程分析软件讨论（下）
#YJK前处理参数详解
【JY】YJK前处理参数详解及常见问题分析（一）
【JY】YJK前处理参数详解及常见问题分析：控制信息（二）
【JY】YJK前处理参数详解及常见问题分析：刚度系数（三）
【JY】YJK前处理参数详解及常见问题分析：二阶效应和分析求解（四）
【JY】YJK前处理参数详解及常见问题分析（五）：风荷载信息
【JY】YJK前处理参数详解及常见问题分析（六）：地震信息
#其他
【JY】位移角还是有害位移角？
【JY】如何利用python来编写GUI？
【JY】今日科普之BIM
~关注未来更精彩~
