1. 简介
在Abaqus/Standard模块中,用户可以利用子程序UEL来开发自定义单元,用以实现一些通过Abaqus内置单元无法实现的功能。如果编写恰当,用户自定义单元可以正常使用Abaqus/Standard的大部分功能,并且在用户子程序UEL中可以编写多个自定义单元,这些单元可以同时使用。
相比于开发完整的有限元程序,使用Abaqus来开发单元的优势是非常明显的:首先,Abaqus提供了大量的结构单元、分析程序以及前处理工具,这些功能可以与Abaqus UEL配合使用;其次,Abaqus中可以完成自定义单元的部分后处理工作;并且,开发子程序的效率也远高于开发完整的有限元程序。
尽管用户子程序UEL为在Abaqus中进行有限元分析提供了非常强大和灵活的工具,但由于开发自定义单元需要较高的数学和力学理论基础,加之前、后处理过程不便,因此相比于诸如UMAT之类的子程序,可供参考的资料非常少。恰好笔者在完成毕设的过程中需要使用到用户子程序UEL进行结构单元的开发,因此借由此机会和大家分享一些子程序UEL的开发经验以及存在的问题,希望可以帮助到有这方面需求的同学。
为了理解用户子程序UEL的工作流程,必须熟悉有限元分析求解的整个过程。因此,本文以最为简单的杆单元为例,介绍杆结构的有限元分析流程,随后利用子程序UEL开发自定义杆单元,并验证计算结果的准确性。
(为了便于排版,部分内容采用图片格式)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 程序用于定义2节点杆单元 C
C 单元节点数:2 节点自由度数:1 单元自由度数:2×1=2 C
C 材料参数: (1)弹性模量E (2)截面积AREA C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
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 子程序主要变量说明 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 需要更新的变量
C
C AMATRX:单元的刚度矩阵 所有非零元素都必须定义
C RHS:单元刚度方程中的右端向量 变量中包含有残余载荷(不平衡力)向量
C SVARS:求解状态变量 状态变量的个数由NSVARS确定(可以不更新)
C ENERGY:用户自定义单元的能量 共有8个分量(可以不更新)
C 传入模型的信息变量(不可修改)
C
C NNODE:自定义单元的节点个数
C JTYPE:单元类型 Un
C NDOFEL:自定义单元的自由度个数
C JELEM:用户指定的自定义单元号
C NSVARS:用户自定义状态变量的个数
C PROPS:单元材料参数实数数组 包含有NPROPS个实数参数
C JPROPS:单元材料参数整数数组 包含有NJPROPS个整数参数
C COORDS:坐标数组 COORDS(K1,K2)为第K2个节点的第K1个坐标
C U:单元计算中的自由度(本单元中为位移)
C DU:位移的增量值
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 变量定义及声明 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
REAL::L !杆单元长度
REAL::E !材料弹性模量
REAL::AREA !杆单元截面积
REAL::EPS !单元应变
REAL::SIGMA !单元应力
INTEGER::I,J !循环索引
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 定义刚度矩阵及右端矢量 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C 获取材料参数
E=PROPS(1)
AREA=PROPS(2)
C 通过单元节点坐标计算杆单元长度
C 节点坐标COORDS(K1,K2)表示当前单元第K2个节点第K1个坐标
L=ABS(COORDS(1,2)-COORDS(1,1))
C 构造刚度矩阵AMATRX
E*AREA/L =
-1*E*AREA/L =
-1*E*AREA/L =
E*AREA/L =
C 构造右端矢量RHS
C RHS=-K・U
DO I=1,NDOFEL !NDOFEL为单元的自由度数 本单元中为2
0 !初始化右端矢量 =
DO J=1,NDOFEL
RHS(I,1)-AMATRX(I,J)*U(J) !按照矩阵相乘的公式计算RHS =
END DO
END DO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 计算单元应力和应变 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C 计算单元应变
EPS=(U(2)-U(1))/L
C 计算单元应力
SIGMA=E*EPS
C 将计算结果储存置状态变量
EPS =
SIGMA =
RETURN
END