蒙皮单元在UEL二次开发中的作用

为了用户能够在abaqus中实现自己的有限元理论，abaqus提供了二次开发接口uel，uel子程序设计用来完成一个单元的核心计算，包括单元对运动方程左端贡献：刚度矩阵、质量矩阵和阻尼矩阵，以及对方程右端不平衡力的贡献，这里要注意区分“方程右端量”和“右端不平衡力”的区别。

uel是一款功能强大的子程序，但同时也是一把双刃剑。用户自定义单元无法直接与abaqus自带的功能完美结合，比如给用户自定义单元施加复杂荷载，面力、体力等，也无法给用户自定义单元定义surface，也就无法与abaqus自带的接触非线性模型相结合。

帖子分享了蒙皮单元在uel二次开发中几个常见的作用，包括可视化、施加复杂荷载和定义接触属性，重在讲解思路。

首先分析用户自定义单元为什么无法与abaqus某些功能直接结合。

从荷载讲起。在有限元计算中，任何荷载都是需要被转化为节点荷载的，无论是线荷载、面荷载还是体荷载，形函数在荷载转化过程中起到了至关重要的作用，这些节点力最终组装到运动方程的右端，参与迭代计算。

但是abaqus主程序是无法得知用户自定义单元的形函数的，因为uel只向主程序输出一些关键矩阵，这些关键矩阵直接被组装到总体模型矩阵中，这么一来，用户自定义单元的形函数只有写程序的人知道，那么复杂荷载也就无法转化到节点上了，因此，用户自定义单元只能定义基于节点的荷载，比如节点位移、节点力和节点速度等等。

与此同时，用户自定义单元的可视化也是无法进行的，因为位移、应变等数据的磨平也需要形函数，这就是用户自定义单元的odb文件中只有“X”形状的原因。如下图

通过以上分析，主要的问题是用户自定义单元无法提供方程的右端量，也就是节点数据。完成abaqus的所有功能，包括但不限于定义surface、体荷载和接触属性等。

那么，只需要在自定义单元原有的节点上面在嵌套一层abaqus内置的单元，被称为蒙皮单元，蒙皮单元的材料属性被设置的无限小，因此蒙皮单元对运动方程的左端量没有贡献，只对右端的荷载向量有贡献，与此同时，也能完成自定义单元的可视化。

设计悬臂梁顶端受节点力算例。悬臂梁尺寸为2×2×8，弹性模量1e10，密度2000，泊松比0.25，悬臂端受垂直轴向的节点力，幅值为1e6，另一端为固定端，边界条件和荷载示意图为

网格尺寸为0.5，共计单元数目256，节点数目为425，网格为

采用静力通用分析步，计算总时长为1，最小增量步长1e-5。（静力计算中的计算时长没有物理意义，只有数值意义）

采用自编的C3D8单元uel与abaqus做对比，设置三种对比工况，分别为：abaqus计算，不带蒙皮单元的uel计算，带蒙皮单元的uel计算，如下图

可以看到，三者计算结果一致。且不带蒙皮单元的云图显示效果很差，仅仅显示了“XX”形状，非常不利于结果展示。第三幅图为带有蒙皮单元的uel计算效果，云图同时显示了“XX”和abaqus自带的单元，把自定义单元隐藏之后的效果为

隐藏了自定义单元之后，显示的云图与abaqus保持一致。同样地，当隐藏了自定义单元以后，u1方向的位移云图为

u2方向的位移云图为

初步说明了蒙皮单元在可视化方面的应用是可行的。

设计悬臂梁一端受垂直于轴向面力算例。计算模型参数与分析步信息与上述相同。边界条件和荷载示意图为

这里为了显示面力，将模型视图调整为透明，可以看到施加的是垂直于轴向的面力。

对该计算模型设置两个对比工况，分别为abaqus计算和带有蒙皮单元的uel计算。因为不带有蒙皮单元无法施加面力，所以没有设置不带蒙皮单元的工况。

统计了加载向的位移云图，不隐藏自定义单元的时候，显示如下

隐藏了自定义单元之后的云图效果为

u1方向隐藏了自定义单元之后的云图

u2方向隐藏了自定义单元之后的云图

以上云图和计算效果说明了蒙皮单元施加面力荷载是可行的。

设计悬臂梁一端受轴向体力算例。计算模型参数与分析步信息与上述相同。边界条件和荷载示意图为

这里为了显示面力，将模型视图调整为透明，可以看到施加的是垂直于轴向的面力。

对该计算模型设置两个对比工况，分别为abaqus计算和带有蒙皮单元的uel计算。因为不带有蒙皮单元无法施加体力，所以没有设置不带蒙皮单元的工况。

统计了加载向的位移云图，不隐藏自自定义单元的时候，显示如下

隐藏了自定义单元之后的云图效果为

u1方向隐藏了自定义单元之后的云图

u2方向隐藏了自定义单元之后的云图

以上云图和计算效果说明了蒙皮单元施加体力荷载是可行的。

设计了滑块在斜面滑落的算例，滑块一共有两个运动过程：在斜面滑动和出斜面运动。该算例采用比例边界有限元法计算，涉及到给用户自定义单元定义重力荷载和接触属性，滑块与斜面的几何尺寸如下图：

滑块与斜面的网格图如下：

该算例中的三种方法均采用动力隐式算法，计算总时长为10s，设置固定增量步长为0.01s，增量步数共计1000。重力加速度为9.81m/s²。具体参数见附件。

滑块两个运动过程受力分析见下图：

上面左图为斜面内的受力分析，右图为斜面外的受力分析。

可知滑块沿斜面的分力与摩擦力大小比值有两种情况：

当滑块在斜面内运动的时，遵循如下控制方程：

通过上式子，可求得位移和速度的解析解。滑块运动水平向位移时程对比如图

误差分析：观察四个时程曲线，发现滑块在斜面上运动的过程中，四种计算结果均吻合良好。而在滑块离开斜面的第二运动过程中，三个数值解吻合较好，但均与解析解有一定的差别。这是因为解析求解是建立在质点的假设上，但是数值解是一个小方块，因此方块第一个角点离开斜面的时候，方块并未完全离开斜面，所以误差会从方块的第一个角点离开斜面开始。

在该算例中，完成荷载和接触定义的同时，附加的蒙皮单元兼具用户自定义单元的位移、速度和加速度可视化功能。

当蒙皮单元和用户自定义单元共存的时候，位移云图为

当只显示用户自定义单元的时候，ABAQUS显示如图

当只显示蒙皮单元的时候，ABAQUS显示如图

这个例子结合了蒙皮单元的所有作用：可视化、施加荷载和接触属性定义，验证了蒙皮单元在uel二次开发中的作用

当然，蒙皮单元也是有限制的。

1、蒙皮单元只适用于用户自定义单元和abaqus内置单元单元节点形状相同的时候。为了避免这种限制，可以将几何拓扑信息复杂的自定义单元拆分为最简单的四面体单元即可。

2、蒙皮单元只能给自定义单元的位移、速度和加速度进行可视化，不能对应变、应力等数据进行可视化。

应力应变等数据的可视化可以采用python处理到第三方软件进行可视化，如sufur、paraview和tecplot等软件。

也可以采用UMAT子程序进行可视化，但是采用UMAT同样要求自定义单元的形状与abauqs内置单元相同。

注：除用户子程序外，上述算例所有inp文件均来源，如有需要，可留言或私信。