ABAQUS UEL 静/动力隐式二次开发

概述：结合HHT时程积分法，推导了ABAQUS 静/动力隐式算法中的UEL关键矩阵表达式，并将公式应用到自编CPS4/CPE4、C3D8 BBAR和C3D20用户自定义单元中，计算结果均与ABAQUS自带单元保持一致。其中，静力计算中，关键矩阵AMARTX和RHS等可直接按照刚度矩阵和方程右端不平衡力输出。动力隐式计算中，这两者的输出较为复杂，需要结合HHT时程积分法进行推导，将刚度矩阵、质量矩阵和阻尼矩阵依据LFLAGS数组的数值进行组合，RHS同样需要进行推导计算，并以合适的方法将解相关的状态量储存在SVARS中，供后面的增量步调用。

关于UEL的程序设计，只支持静力通用计算分析步的资料有很多，本帖子内容适不但适用于静力通用，而且适用于动力隐式、频率分析等分析步计算，采用模块化程序设计，所有的矩阵求解均被封装，调用方便，读者可以针对自己的需求对相应函数进行改编即可，尤其是涉及到动力隐式计算部分，适用于任何运动方程的动力隐式求解，可直接移植使用。

----------------------------

（）UEL接口基本参数介绍

----------------------------

用户自定义单元(USER DEFINED ELEEMNT, UEL)适合进阶的工程师/学者使用，u、UEL的目的是实现一个单元的力学行为，即力-位移关系，这部分的功能需要的编程工作量和理论功底都比较高， 涉及到的知识大致包括：

有限元理论：形函数插值、应力-应变关系，刚度矩阵组装，质量矩阵求解，阻尼矩阵求解，

数值算法：高斯积分（全积分、缩减积分）、非线性方程组迭代求解流程（增量迭代，常刚度迭代等）、雅各比矩阵求解

弹/塑性力学：应力-应变关系，位移-应变关系

熟练的FORTRAN编程操作是基础技能。

UEL接口参数众多，包括自己本身的参数和与其他子程序联合的参数，下面介绍我比较熟悉的参数。

<>需要编程定义的参数

RHS（right hand side）：这个命名是从方程组的角度来的，顾名思义，他就是方程的右端量，其本质是：外力-内力，外力部分程序的编写涉及到与其他子程序的联合使用，包括DLOAD和UTRACLOAD等等，这部分内容是给用户自定义单元施加复杂的广义力，UEL接口为其提供了相应的参数，如JDLTYP、NDLOAD等等，目前没有做过尝试。至于内力的求解，在静力线性计算中，数值上等于-KU（后面有这个公式的推导），即刚度与位移乘积负数。在动力隐式分析中，这个矩阵需要结合HHT时程积分法推导具体的表达式。

AMATRX：这个参数往往被认为是刚度，但并不这样，他只有在静力计算的时候才是刚度。他的具体取值依据分析类型的不同而不同，具体的表达形式不唯一。在动力隐式分析中，这个矩阵需要结合HHT时程积分法推导具体的表达式。

SVARS：这个参数官方文档说是取决于结果的状态量，他的具体意义和数据由我们自己确定，而我把它理解为一个小仓库，可以存放我们的数据。需要注意的是里面的数据可以在不同的增量步之间传输，就是说，里面的数据会传到下一个增量步，只要不更新他，他可以一直被使用。这个数据我们自己决定更新与否。

ENERGY：与用户自定义单元相关的一些动量等能量，也是依据不同的分析类型而不同。

<>可以被编程更新的参数

PNEWDT：这个参数可以用来修改当前增量步长，

NDOFEL：单元的自由度数目，如平面四节点单元，自由度为8。这个参数是主程序传给我们用于计算的，我们直接使用，不更新。

NRHS：荷载向量的数目，对一般的分析，数值位移，对特殊的分析类型，数值具体而定。

CROODS，MCRD和NNODE：描述单元的坐标等信息，croos储存单元的初始坐标，mcrd可简单的理解为单元节点的自由度数目，nnode为单元的节点总数。

U，DU，V和A：这一系列参数是ABAQUS对该数据在当前增量步结束时候的估计。

TIME（1）：为当前分析步时间，TIME(2)为计算总时间。DTIME为时间增量步长。

KSTEP和KINC：分别为当前分析步和增量步个数。

JELEM：为用户自定义单元编号。

<>传入UEL作为信息的参数

PROPS：包含INP文件中定义的浮点数组，这些数据需要我们提前在inp文件中定义好，然后在UEL是使用，他的大小为NPROPS。

JPROPS：包含整型单元属性的数组，他的大小为NJPROP。

COORDS：单元节点坐标，列代表节点，行代表维度

U、DU、V、A：位移等数据量

U（K1）：第K1个自由度的位移

DU（K1，KRHS）：没用过，不翻译

V（K1）：第K1个自由度的速度

A（K1）：第K1个自由度的加速度

JDTYPE：没用过，不翻译。只知道这个参数是用于给用户自定义单元施加复杂荷载的

ADMAG：没用过，不翻译。

DDLMAG：没用过，不翻译。

PREDEF：没用过，不翻译。

PREDEF（K1，1，K3）：没用过，不翻译。

PREDEF（K1，2，K3）：没用过，不翻译。

PREDEF（K1，3，K3）：没用过，不翻译。

PREDEF（K1，K2，K3）：没用过，不翻译。

PREDEF（1，K2，K3）：没用过，不翻译。

PREDEF（1，K2，K3）：没用过，不翻译。

PARAMS：求解运动方程，积分法里面用到的参数

PARAMS（1）：α

PARAMS（2）：β

PARAMS（3）：γ

LFLAGS：一个至关重要的数组，在动力UEL编程中，向ABAQUS主程序输出的关键矩阵，就是依据这个数组。该数组内容较多，详细的内容放在后面讲解。

TIME（1）：当前分析步时间。

TIME（2）当前分析总时间

<>标量参数

DTIME：增量步长

PERIOD：当前分析步时间

NDOFEL：单元的自由度数目

MLVARX：没用过，不翻译

NRHS：右端量数目

NSVARS：用户自定义的解状态量数目，与USER ELEMENT关键字后面的参数对应

NPROPS：单元属性参数个数，针对浮点数

NJPROP：单元属性参数个数，针对整数

MCRD：用户自定义单元节点维度个数

NNODE：用户自定义单元节点数目

JTYPE：用户自定义单元合在标识符，用于给用户自定义单元施加复杂荷载（这部分应用较复杂）

KSTEP：当前分析步数目

KINC：当前增量步数目

JELEM：用户自定义单元编号

NDLOAD：自定义合在标识符

MDLOAD：自定义荷载总数目

NPREDF：自定义场变量编号

 () CPS4/CPE4公式推导

等参单元中，母单元向笛卡尔坐标的转换为：（注意：这里的箭头标反了）

要再现位移场，可假设单元内部位移为坐标的函数，即：

将四个节点的位移和坐标代入上式，有：

如此，求解系数，并用节点位移表示单元内部位移：

上式即为单元节点位移与坐标的插值关系。考虑到等参单元，即位移的插值函数与母单元-笛卡尔单元变换采用相同插值函数。则等参变化的插值函数为：

上式的矩阵表达为：

其中：

由弹性力学几何方程，应变矩阵为：

合并为：

JOCABIN矩阵为：

其中：

则JOCABIN的具体表达式为：

刚度矩阵的求解：

质量矩阵分为两种，集中质量矩阵和协调质量矩阵，下面给出的是一致质量矩阵（又叫协调质量矩阵）矩阵，该种质量矩阵因其采用了与位移插值相同的插值函数而被称为一致质量矩阵。但是协调质量矩阵为满阵，不利于数值计算，所以有限元中一般采用集中质量矩阵，即只有对角元素不为零。集中质量矩阵通常由协调质量矩阵处理而来，比如将协调质量矩阵每一行的质量集中到对角元素，然后将非对角元素数值置为零。

阻尼矩阵，采用比例阻尼：

UEL主程序采用常刚度法设计，刚度矩阵和质量矩阵只在第一个增量步计算一次，然后储存在变量SVARS中，后面一直到分析结束之间的所有增量步不再计算，均直接读取SVARS数组中的数据并输出会ABAQUS主程序。该主程序将各个关键量的求解分为不同的子程序，这些子程序基本是所有用户自定义单元都会涉及的，因此该主程序的可移植性很强，在线弹性问题中，即可用常刚度迭代法求解的运动方程中，读者完全可以按照这个直接编写自己的程序。

不想写文字了，麻烦，直接把付费内容放在附件里面了：

（1）CPE4/CPS4单元UEL，适用于静/动力隐式计算

（2）C3D8 BBAR单元UEL，，适用于静/动力隐式计算

（3）C3D20 单元UEL，适用于静/动力隐式计算

（4）WORD文档，推导、讲解动力隐式计算中的AMARTX、RHS

赠送内容：

（5）收集的一些UEL-UMAT程序