FORTRAN 中的 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
变量说明:
1、abaqus 传入的模型信息的变量:
(1)数组:
PROPS(材料参数)
浮点数(实数)数组,包含用户单元 NPROPS 个待用实数型常量。NPROPS 为用户指
定的单元材料参数中的实数常量个数。
JPROPS(材料参数)
整型数组,包含用户单元 NJPROP 个待用整数型常量。NJPROP 是用户指定单元材料参
数中的整数常量个数。
COORDS
包含单元结点原始坐标的数组。COORDS(K1,K2)表示单元结点 K2 的 K1 坐标值。
U,DU,V,A
包含当前增量步结束时单元结点上基本求解变量当前估计值的数组(依赖于位移,旋转,
温度等自变量)。值如下:
U(K1) 变量全部的值(如位移)。如果是线性摄动分析步,基本状态时的值。
DU(K1,KRHS)
KRHS 时当前增量步的变量增加量(如增量位移)。如果是特征值提
取分析步,表示特征向量 KRHS 的特征向量大小。对于稳态动力,KRHS=1 表示摄动位移
的实部,KRHS=2 表示摄动位移的虚部。
V(K1) 变 量 的 时 间 一 次 导 数 ( 速 度 , 转 速 )。 只 能 在 隐 式 动 力 分 析 中 定 义
(LFLAGS(1)=11or12)。
A(K1)
变 量 的 时 间 二 次 导 数 ( 加 速 度 )。 只 能 在 隐 式 动 力 分 析 中 定 义
(LFLAGS(1)=11or12)。
JDLTYP
用于为单元定义分布式荷载类型的整数数组。荷载类型 Un 通过 JDLTYP 中的整数 n 来
识别;荷载类型 UnNU 通过 JDLTYP 中的负整数-n 来识别。JDLTYP(K1,K2)表示在第 K2 加
载情况下第 K1 种分布式荷载。对于一般非线性分析步,K2 通常为 1.
ADLMAG
对于一般非线性分析步,ADLMAG(K1,1)表示对分布荷载类型 Un 在当前增量结束时第
K1 种分布荷载的合力大小。对于分布荷载类型 UnNU,荷载大小在 UEL 中定义;因此,
ADLMAG 中相应项为零。对于线性摄动分析步,ADLMAG(K1,1)表示类型 Un 应用于基本
状 态 的 第 K1 种 分 布 荷 载 的 合 力 。 基 本 状 态 荷 载 类 型 UnNU 必 须 在 UEL 中 处 理 。
ADLMAG(K1,2),ADLMAG(K1,3)等目前没有使用。
DDLMAG
对于一般非线性分析步,DDLMAG 表示单元上当前活跃的分布荷载类型为 Un 的分布
荷载的大小值的增量。DDLMAG(K1,1)表示当前时间增量的荷载大小增量。在计算外部作业
贡献时需要荷载大小增量。对于分布荷载类型为 UnNU,荷载大小在 UEL 定义;因此,
DDLMAG 中相应项为零。对于线性摄动分析步,DDLMAG(K1,K2)表示对分布荷载类型 Un,
当前单元上活跃的分布荷载大小的摄动。K1 单元上活跃的第 K1 种摄动荷载。K2 通常为 1.
对于稳态动力分析例外。K2=1 表示真实荷载,K2=2 表示假想荷载。摄动荷载类型 UnNU
必须在 UEL 中处理。
PREDEF
包含预定义场变量值的数组,如在非耦合应力位移分析中单元结点上的温度。
数组第一个索引 K1 是 1 或 2,(1 表示在增量步结束时场变量值,2 表示场变量增量)。
第二个索引 K2 表示变量:温度相当于索引 1,预定义场变量相当于索引 1 和 2.在温度没有
定义的情形下,预定义场变量以索引 1 开始。第三个索引 K3 表示单元的局部结点号。
PREDEF(K1,1,K3) 温度
PREDEF(K1,2,K3) 第一个预定义场变量
PREDEF(K1,3,K3) 第二个预定义场变量
…
PREDEF(K1,K2,K3) 第 K2 个预定义场变量在单元第 K3 个结点上的全部或增量值
PREDEF(1,K2,K3) 在当前增量步结束时变量值
PREDEF(2,K2,K3) 对应于当前时间增量的增量值
PARAMS
存储 abaqus 内部与分析求解程序相关的参数的数组。数组中各项依赖于 UEL 调用时当
前求解程序,通过数组 LFLAGS 中各项来识别(如下)。
LFLAGS
非常重要的标示变量数组。反应当前分析步骤的类型,以及单元计算需要进行的内容,
这些内容必须在 UEL 中完成。Abaqus/Standard 程序的详细需求在本章早期定义。
定义程序类型
对于隐式动力(LFLAGS(1)=11or12))PARAMS 表示如下:
LFLAGS(1)
LFLAGS(1)=1
LFLAGS(1)=2
LFLAGS(1)=61
LFLAGS(1) =
分析过程为静态分析
分析过程为流体渗透/应力耦合分析
分析过程为静力分析,自动时间增量步长
分析过程为静力分析,固定时间增量步长
小位移分析
大位移分析
通用隐式求解时间增量过程,必须在 RHS 中定义残余向量,在
AMATRX 中定义雅克比矩阵
只定义当前刚度矩阵
只定义当前阻尼矩阵
只定义当前质量矩阵 ,Abaqus/Standard 总是在分析开始时要
求初始质量矩阵
只定义当前残余向量或荷载向量
为初始加速度计算(或撞击后的加速度计算)定义当前质量矩
阵和残余向量
定义输出扰动量
分析步为一般分析步
分析步为线性摄动分析步
当前近似值,取决于牛顿修正值
当前近似值通过之前增量外推得到
62~65
LFLAGS(2)=0
LFLAGS(2)=1
LFLAGS(3)=1
LFLAGS(3)=2
LFLAGS(3)=3
LFLAGS(3)=4
LFLAGS(3)=5
LFLAGS(3)=6
LFLAGS(3)=100
LFLAGS(4)=0
LFLAGS(4)=1
LFLAGS(5)=0
LFLAGS(5)=1
TIME(1)
当前子步时间
TIME(2)
当前全局时间
(2)标量参数
DTIME 时间增量
PERIOD 当前分析步历时
NDOFEL 单元自由度号
MLVARX 标注参数(在多位移或右手边向量被使用的情况下使用)
NRHS 荷载向量数。在大多数非线性问题中为 1;在修正 Risk 静力程序中为 2,在一
些线性分析程序和子结构中大于 1
NSVARS 用户定义的与单元相关的求解依赖状态变量个数
NPROPS 用户定义的与单元相关的实常数个数
NJPROP 用户定义的与单元相关的整数型常数个数
MCRD 用户定义的结点所需的最大坐标号和用户单元最大激活自由度,通常小于等
于 3. 例如,如果指定最大坐标号为 1 和用户单元激活自由度为 2,3 和 6,那么 MCRD 将
为 3. 如果指定最大坐标号为 2 和用户单元激活自由度为 11 和 12,那么 MCRD 将为 2.
NNODE 单元上用户指定的结点个数。
JTYPE 定义单元类型的整数。在单元类型 Un 中指定整数值 n。
KSTEP 当前子步号。
KINC 当前增量号。
JELEM 用户分配的单元号。
NDLOAD 当前单元上激活的分布荷载或流动的标识数
MDLOAD 定义在单元上的分布荷载和(或)流动总数
NPREDF 预定义场变量号,包括温度。Abaqus/Standard 为用户单元每个场变量每个
结点使用一个值。
2、(必须更新的量)依赖于数组 LFLAGS
值的矩阵
RHS
该数组包含用户自定义单元对总体系统方程右手边(right-hand-side)向量的贡献。对
于多数非线性分析,NRHS=1,RHS 包含残余载荷(不平衡力)向量。在修正的 Risks 静态
分析时例外,此时 NRHS=2,RHS 第一列包含残余向量,第二列包含作用于单元上的外荷
载的增量。
AMATRX
该数组包含了用户单元对雅克比刚度矩阵或其他的总体系统方程矩阵的贡献。何时为特
殊的矩阵依赖于在数组 LFLAGS 中的项数。无论 AMATRX 是否是对称的,都需要定义数组
AMATRX 中的非零项。如果在定义用户单元时没有指定矩阵非对称,Abaqus/Standard 将默
认为对称阵,即式中[A]子程序中定义的矩阵 AMATRX,并存储 1/2([A]+ [A]T)。如果指定
为非对称的,Abaqus/Standard 将直接用 AMATRX。
SVARS
该数组包含用户自定义单元求解所依赖的状态变量值。这些变量个数用 NSVARS 表示。
由用户定义这些变量的意义。
对一般非线性分析步,当前增量开始时该数组被传递给包含这些变量值的 UEL 子程序。
在增量结束时,这些值应该被更新,除非程序在调用 UEL 子程序时没有要求更新。这依赖
于 LFLAGS 数组中的项。对线性摄动分析步,在基本状态时,该数组被传递给包含这些变
量值的 UEL 子程序。如果你希望输出这些工作量,包含摄动值的数组就应该返回。
当 KINC 等于零时,调用 UEL 用于零增量输出。在这种情况下,返回值将只用于输出
用途,并且不再更新。
ENERGY
对一般非线性分析步,该数组包含用户单元的能量值。当调用 UEL 时数组中的值就是
当前增量起始时的单元能量。当前增量结束时值应该被更新。对线性摄动分析步,在基本状
态时将包含能量的数组传递给 UEL。如果希望输出这些工作量,含有摄动值的量应该被返
回。不定义不影响计算结果。数组中的项目如下:
ENERGY(1)
ENERGY(2)
ENERGY(3)
ENERGY(4)
ENERGY(5)
ENERGY(6)
Kinetic energy.动能
Elastic strain energy.弹性应变能
Creep dissipation. 蠕变损耗
Plastic dissipation. 塑性损耗
Viscous dissipation. 粘性损耗
“Artificial strain energy” associated with such effects as artificial stiffness
introduced to control hourglassing or other singular modes in the element.
用于沙漏控制或其他单一模型
ENERGY(7)
ENERGY(8)
Electrostatic energy. 静电能
Incremental work done by loads applied within the user element.
施加到自定义单元内的载荷的功的增量
3、需要更新的变量:
PNEWDT
建 议 新 的 时 间 增 量 与 当 前 被 用 时 间 增 量 (DTIME ) 的 比 值 。 该 变 量 允 许 用 户 在
Abaqus/Standard 中为自动时间增量算法提供输入(选择自动时间增量时)。该值只对平衡迭
代(正常时间增量,指明 LFLAGS(3)=1)起作用。在严重不连续迭代中,该值将被忽略,
该分析步中指定 CONVERT SDI=YES。用法如下:
在每次调用 UEL 前,PNEWDT 值被设为很大的值。
如果 PNEWDT 重定义的值小于 1.0,那么 Abaqus/Standard 必须放弃时间增量,并且尝
试 用 更 小 的 时 间 增 量 再 试 一 次 。 提 供 给 自 动 时 间 整 合 算 法 的 新 的 建 议 时 间 增 量 为
PNEWDT*DTIME。此处 PNEWDT 是所有调用用户子程序(允许在这次迭代中重新定义
PNEWDT)的最小值。
如果当前迭代用户子程序 PNEWDT 给定的值大于 1.0,并且迭代中增量收敛,那么
Abaqus/Standard 可能会增加时间增量。提供给自动时间整合算法的新的建议时间增量为
PNEWDT*DTIME。此处 PNEWDT 是所有调用用户子程序(允许在这次迭代中重新定义
PNEWDT)的最小值。
如果在分析程序中没有选择自动时间增量,PNEWDT 大于 1.0 将被忽略,PNEWDT 小
于 1.0 则将终止计算。
INP 中的单元说明:
需要在 abaqus 中修改 inp 文件,以调用上述 uel 程序。在此 inp 文件中需定义单元材料参数
个数、坐标、状态变量个数(应力、应变、应力增量、应变增量…)、单元编号、节点个数、
编号、自由度等参数。如:
*USER ELEMENT,NODES=2,TYPE=U1,PROPERTIES=3,COORDINATES=3 单元节点数、单元类
型、材料参数个数、坐标
1,2,3,8
*ELEMENT,TYPE=U1,elset=pvd 指定节点、单元类型、单元集
1001,253,248
1002,248,243
1003,243,238
1004,238,233
1005,233,228
1006,228,223
1007,223,218
1008,218,213
1009,213,208
1010,208,203
1011,203,198
1012,198,193
1013,193,188
1014,188,183
1015,183,178
1016,178,173
*UEL PROPERTY,ELSET=pvd 定义材料参数
0.0314,10,1e-2 面积、弹模、渗透系数