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

不等待


即关注

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

【简述ABAQUS中UEL子程序】

ABAQUS作为成熟的商用有限元软件,可为高级用户提供特定的分析需求。ABAQUS常见的二次开发子程序包括:UMAT、VUMAT、UGENS、UEL和VUEL等。其中UEL/VUEL分别适用于ABAQUS的Standard/Explicit求解器。只有清楚有限元分析的基本原理,才能够较好地了解其分析的力学原理,才能对特定的分析需求编写合适的分析单元。
本文对ABAQUS/Standard中UEL子程序的有限元原理及编写原则进行初步梳理,并编写了平面三角形单元静力分析的UEL子程序,以加深对UEL的理解与认识。
【JY】 ABAQUS子程序UEL的有限元原理与应用的图3
【JY】 ABAQUS子程序UEL的有限元原理与应用的图4

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

当用户需要用ABAQUS/Standard单元库中没有的单元进行分析时,可通过ABAQUS提供的UEL子程序接口进行二次开发,编写适用于特定分析的单元。下面以平面三角形单元为例,给出编写UEL单元的基本流程。
【JY】 ABAQUS子程序UEL的有限元原理与应用的图5
系统的势能Π由下式给出
【JY】 ABAQUS子程序UEL的有限元原理与应用的图6

其中【JY】 ABAQUS子程序UEL的有限元原理与应用的图7为单元的应变能。

将应变-位移关系ε = BqUe中,由此对于势能方法形成的刚度矩阵如下:

【JY】 ABAQUS子程序UEL的有限元原理与应用的图8
【JY】 ABAQUS子程序UEL的有限元原理与应用的图9
其中 t e 为单元厚度、 A e 为单元面积。
有限元分析的第一步是将实体离散化为多个单元,此后构造出各个单元的单元刚度矩阵,在将单元刚度矩阵集成整体的刚度矩阵,最后通过整体刚度矩阵建立平衡方程从而求解各个节点的位移、应力、应变等响应。因此,根据有限元的分析原理,编写UEL的最终目的即是形成目标单元的单元刚度矩阵。以下对平面常应变三角形单元(CST)进行有限元分析,依据有限元分析思路编写UEL子程序,以初步了解有限单元的分析原理及UEL编写流程。
【JY】 ABAQUS子程序UEL的有限元原理与应用的图10

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

上图所示为平面三角形单元,由图示信息可得其节点坐标列阵 u 与位移列阵 q 分别为
【JY】 ABAQUS子程序UEL的有限元原理与应用的图11

1.形函数

设三角形单元节点1、2、3对应的形函数分别为N1、N2、N3,形函数满足条件:N1+N2+N3=1(在点i处,形函数Ni值为1,其余形函数值为0),形函数采用自然坐标 ξ η 描述,有
【JY】 ABAQUS子程序UEL的有限元原理与应用的图12

2.高斯积分点及权重

此平面三角形单元所选取的高斯积分点及权重如下:
表1高斯积分点及权重
【JY】 ABAQUS子程序UEL的有限元原理与应用的图13

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

得到形函数与单元节点位移后,根据等参元表示方法,单元内任意一点的位移都可用形状函数和未知结点位移场进行表示,有
【JY】 ABAQUS子程序UEL的有限元原理与应用的图14
可表示为u= N q ,其中N
【JY】 ABAQUS子程序UEL的有限元原理与应用的图15
对于三角形单元,可将x、y坐标表示为节点坐标的形式,可由此得到单元坐标的插值关系为:
【JY】 ABAQUS子程序UEL的有限元原理与应用的图16
其中
【JY】 ABAQUS子程序UEL的有限元原理与应用的图17

4.雅各比矩阵J

由节点坐标x、y与ξ、η关系,在根据链式求偏导法则,可得到平面三角形单元的雅可比矩阵
【JY】 ABAQUS子程序UEL的有限元原理与应用的图18

5.应变矩阵B 

根据应变ε与节点坐标u、节点位移q的关系可得应变计算公式为
【JY】 ABAQUS子程序UEL的有限元原理与应用的图19
 其中xij与yij计算公式同上。写成矩阵形式为
【JY】 ABAQUS子程序UEL的有限元原理与应用的图20
由此可导出单元应变-位移矩阵B,其表达式为
【JY】 ABAQUS子程序UEL的有限元原理与应用的图21
由此可看出B矩阵中所有元素均为常数,且每个常数均是通过节点坐标表示的。

6.弹性矩阵D

弹性矩阵D根据所研究问题为平面应力或平面应变问题有不同形式,其取值需根据所用材料确定,具体如下:
【JY】 ABAQUS子程序UEL的有限元原理与应用的图22
【JY】 ABAQUS子程序UEL的有限元原理与应用的图23

7.应变-位移矩阵

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

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

根据势能方法导出得到三角形单元的单元刚度矩阵 k e ,有
【JY】 ABAQUS子程序UEL的有限元原理与应用的图24
其中 t e 为单元厚度
最后,通过单元的边界条件可得单元的残余力向量。至此,UEL单元的编写完成,可进行初步调试直至结果无误。
【JY】 ABAQUS子程序UEL的有限元原理与应用的图25
【JY】 ABAQUS子程序UEL的有限元原理与应用的图26

【 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

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子程序应用于较多单元的分析案例中,可保证分析结果的准确性。
【JY】 ABAQUS子程序UEL的有限元原理与应用的图27
【JY】 ABAQUS子程序UEL的有限元原理与应用的图28

【在ABAQUS中的UEL单元研究】

根据以上平面三角形单元的有限元分析思路,编写对应的平面二维三角形单元UEL子程序,并通过两个不同的有限元分析算例验证该子程序的有效性。算例中采用的三角形单元为ABAQUS单元库中的二维三角形单元CPS3。
【JY】 ABAQUS子程序UEL的有限元原理与应用的图29
【JY】 ABAQUS子程序UEL的有限元原理与应用的图30

算例一:平面三角形单元

一、三角形单元有限元计算参数

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

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

三角形单元的相关计算参数如表2所示:
表2 三角形单元计算参数
【JY】 ABAQUS子程序UEL的有限元原理与应用的图32
该三角形单元共有三节点(如图1所示),每个节点有2个自由度
二、边界条件及载荷设置
将节点1与节点3的位移自由度施加约束,在节点2的水平x方向施加100N的节点力。
三、计算结果对比
在相同边界条件、载荷条件下,二维三角形单元的UEL子程序计算结果与有限元ABAQUS计算结果如下表3及图4-5所示,从图中可看出,子程序的位移计算结果与abaqus中CPS3单元位移计算结果一致。
表3 UEL子程序与abaqus有限元最大位移对比
【JY】 ABAQUS子程序UEL的有限元原理与应用的图33
【JY】 ABAQUS子程序UEL的有限元原理与应用的图34

(a)CPS3加载时程

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

  (b)UEL加载时程

图4两种不同单元的节点位移时程曲线

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

图5 ABAQUS杆单元计算结果与UEL单元计算结果对比

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

算例二:平面三角形单元组合

一、三角形单元有限元计算参数

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

图6平面三角形结构

由上图6可知该结构为四个相同的二维三角形单元组成。该杆件结构的相关参数如下表4所示:

表4 三角形单元参数
【JY】 ABAQUS子程序UEL的有限元原理与应用的图41
该平面三角形共有六个节点(如图4所示),共由三个平面二维三角形单元组成,每个节点有2个自由度。
二、边界条件设置
在三角形节点3施加100N水平力、节点4施加20N的水平力,并约束节点1、节点2的所有自由度,其余节点均为自由边界。
三、计算结果对比
在相同边界条件、相同载荷条件下,所编写的子程序位移计算结果与abaqus中的CPS3单元位移计算结果如下表4及图5-6所示,可从下图表得知,二者计算结果一致。
表5 UEL子程序与abaqus有限元最大位移对比
【JY】 ABAQUS子程序UEL的有限元原理与应用的图42
【JY】 ABAQUS子程序UEL的有限元原理与应用的图43

(a)CPS3加载时程-U1

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

(b)UEL加载时程-U1

图7-1两类单元节点的U1位移时程曲线

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

(a)CPS3加载时程-U2

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

(b)UEL加载时程-U2

图7-2两类单元节点的U2位移时程曲线

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

图8 CPS3单元计算结果与UEL单元计算结果对比

通过以上两个有限元算例可知,采用所编写的UEL二维三角形单元子程序所计算得到的位移及时程结果与有限元软件ABAQUS中采用CPS3单元对应计算结果一致,说明该三角形单元子程序可靠且合理。

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

【参考文献】

[1]:Abaqus6.14 User Subroutines Reference Guide

[2]: 工程中的有限元方法(第四版)

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

概念为先,机理为本,期待下篇!

 往期推荐 ·


#性能分析

【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


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

 ~关注未来更精彩~


(3条)
默认 最新
哈哈哈哈哈
评论 点赞
感谢分享
评论 点赞

查看更多评论 >

点赞 13 评论 3 收藏 29
关注