1 uel简要介绍

1.1 uel变量

在uel中，主要需要更新amatrix（雅克比矩阵）和rhs（右端残值项），如果有状态变量的话，也需要更新状态变量svars.

1.2 inp文件修改

并且要成功使用uel的话，需要对inp文件做以下修改，一一说明。

（1） 需要首先定义uel，加入以下keyword以及data line(s)

各变量的含义如下：

该keyword下的data line(s)表示uel激活的自由度。

（2） 定义完uel之后，可以使用该uel，keyword以及data line(s)如下：

（3） 接着需要给uel赋予性质

以C3D8单元计算一个各向同性材料为例，uel需要的性质为材料的杨氏模量和泊松比。

1.3 uel结果的可视化

由于Abaqus不支持uel结果的可视化，因此我们可以利用umat来辅助进行可视化。使用Abaqus的标准单元覆盖uel，即标准单元与uel共享节点，但是材料使用umat（设置一个很小的刚度），关键步骤在于使用common block将uel中的svars数据传递到umat中的statev中。

2 有限元离散方程

2.1 形函数

本部分记录八节点实体单元（C3D8）的刚度矩阵和内力列阵等 ，熟悉此部分可跳过。

该单元的节点编号顺序示意图如下：

单元的形函数矩阵为：

参数坐标系下的形函数为：

因此形函数对参数坐标的导数为：

2.2 雅克比矩阵

利用形函数对参数坐标的导数矩阵和节点的物理坐标信息可以求得雅克比矩阵：

具体公式为：

2.3 B矩阵

利用雅克比矩阵可以求得形函数对物理坐标的导数：

从而可以组装B矩阵为：

2.4 单元刚度矩阵和内力列阵

单元的刚度矩阵和内力公式为：

uel的主要目的就是更新这两项，其中和的计算需要用到具体的本构方程。在此使用各向同性线弹性本构。

2.5 高斯积分

上述积分的计算需要使用数值积分，一维高斯积分点的坐标以及权系数如下：

采用2x2x2的高斯积分点。

3 B bar 算法

B-bar算法将应变修正为:

式中:

体应变插值B矩阵为:

记为:

也有:

因此将B矩阵修正为:

4 测试

4.1 一个单元测试

对一个单元进行单轴拉伸测试，借助umat进行可视化时，在代码中设置参数为：

需要根据job的单元总数手动修改该参数。

位移结果为：

应变结果为：

应力结果为：

4.2 带孔板单轴拉伸测试

位移结果对比为（左列为Abaqus计算结果，右列为uel计算结果）：

应变的对比结果为（左列为Abaqus计算结果，右列为uel计算结果）：

应力的对比结果为（左列为Abaqus计算结果，右列为uel计算结果）：

uel源代码

2024/03/27更新:已加入B-bar算法进行刚度矩阵的修正,可看到uel的计算结果已于Abaqus计算结果一致