不等待

即关注

【简述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)
CINCLUDE 'ABA_PARAM.INC'
CDIMENSION 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
（此部分是我们所要编写的内容）RETURNEND

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

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

 ~关注未来更精彩~