基于ABAQUS UEL二次开发的比例边界有限法

分享一篇关于ABAQUS二次开发的文章，这是比例边界有限法与商业软件结合的第一篇文章，来自于CMAME期刊，很有意义，文章配套的有开源程序，方便学习。文章和开源程序的链接为

比例边界有限元法基本理论在这里不做过多介绍，可参考相关文献。这里重点讲解比例边界有限元法如何通过UEL等子程序嵌入到ABAQUS中。

有限元求解过程中，单元刚度矩阵的求解和组装是重中之重，这些过程经过商业软件的gui封装，我们已经看不到，只需要点击鼠标就能完成，看起来非常的简单，但这种简单是封装的结果，其内部程序是非常复杂的，abaqus软件主要将fem的算法封装到软件中，但随着有限元理论的发展，各种新型理论层出不穷。作为商业程序，abaqus需要迎合市场需求，因此，abaqus向从业者提供了一系列的二次开发接口，用户可以值只关注自己的科学问题，在编程的过程中能够大幅度的避免重复造轮子，将更多的精力用于优化自己专业部分的理论和程序。

其中uel可用于实现刚度矩阵等计算，当abaqus主程序计算单元刚度矩阵的时候，abaqus不再调用自己内部的程序，而是调用用户开发的子程序 ，子程序执行结束会向主程序返回刚度矩阵、质量矩阵和阻尼矩阵等等，abaqus用于迭代求解。下面是文章给出的带有子程序的abaqus执行流程图

sbfem的优势之一是可以直接处理任意多面体的单元，如下面的单元

除了编写uel子程序，还需要手动修改inp文件，对于上面的多面体单元，inp文件中的代码为

其中，user element用于声明一个用户自定义单元，nodes指明了单元的节点数目，type相当于给这种单元指明了单元类型，u8中的数字8会以参数的形式传进子程序里面，可用于识别不同的单元类型，properties指的是单元的浮点型材料属性个数，coordinates代表该单元激活的自由度数目，variables是单元的状态变量个数，与单元的节点数目有关，详细的计算方法可查看文章正文。

element用于声明单元的编号和单元节点信息，需要特别注意的是，后面的节点排序是非常有讲究的，我在另一篇帖子详细写过，具体可查看

以六面体单元为例，在ABAQUS中，前四个节点和后四个节点分别组成了六面体的两个对面，其中第15、26、37和48节点是在同一条棱上的，而且要求第一个面的法向是指向单元内部的，那么后一个面的法向就是指向单元外部的，这样依靠节点的排序就可以重现单元的几何拓扑信息。

但是在任意多面体中，上述排布方式就不能完全重现单元的几何拓扑信息，那就需要补充外部文件来导入数据，在编程的时候读取这个文件，为此，文章设计了一种txt文件用于补充说明信息，

这张图片就是上面两个多面体单元的txt文件，txt文件包括四种信息。

首先是模型的节点坐标信息，这部分数据的第一行是模型的总体节点数目n，后面是n行节点的坐标。

第二部分是模型的面信息，这部分数目的第一行是模型的总面数目，假如是一个六面体，那这个数目就是6，后面6行就是六面体的面节点信息，但是面节点信息是有顺序的，采用面的法向表示，即第四部分数据。

第四部分数据的第一行是模型的单元总数n，后面是n行数据，每一行的第一个数字是该单元的面数目，后面的数字，是面在第三部分数据中的索引，前面的正负号则代表了面的法向，负号代表面法向是指向单元外部，反之则是指向单元内部。

文章通过子程序UEXTERNALDB将上述txt数据导入到程序中，UEXTERNALDB子程序可以在分析的不同位置被调用，如分析开始前、分析进行中和分析结束等止时刻被调用。

因为软件的设计机制问题，目前的用户自定义单元并不能无缝的与abaqus自带功能相结合，如用户自定义单元无法施加面荷载、无法定义面，也无法完成单元的可视化和接触属性的定义。

根本的原因是abaqus主程序无法得知用户自定义单元的形函数，因此面荷载无法插值到节点上面，也就无法定义面，则abaqus所有涉及到面的功能均无法在用户自定义单元上面使用，接触模型、界面模型也无法使用，为了解决这个问题，文章采用了蒙皮单元技术。

蒙皮单元技术具体实行的操作非常的简单，只需要在用户自定义单元的节点上再嵌套一层abaqus自带的单元即可，称之为蒙皮单元，蒙皮单元属性需要设置为无限小，然后将面荷载定义在蒙皮单元上，或者在蒙皮单元上定义接触属性，后面的计算就可以顺利进行。

蒙皮单元的原理涉及到运动方程。运动方程等号左边的刚度矩阵、质量矩阵和阻尼矩阵通过uel计算，因为蒙皮单元的材料属性无限小，因为对方程的左端量贡献可以忽略，右端力向量通过abaqus蒙皮单元提供，abaqus主程序将蒙皮单元的力转换到节点上面，而节点不区分用户自定义单元和abaqus内置单元。

如下图所示

外部的六面体单元为用户自定义单元，内部的五面体单元为蒙皮单元，即abaqus自带的C3D5单元，这种单元从abaqus2017开始有，之前的版本没有，右面的代码为对应的inp文件，可以看到面定义是基于蒙皮单元的。

蒙皮单元的形式不仅仅只有五面体，理论上来讲，abaqus所有内置的单元均可以作为蒙皮单元，只要保证蒙皮单元的节点与用户自定义单元重合、并且材料属性无限小即可。

文章设计悬臂梁算例用与探究子程序在动力问题中的计算性能，悬臂梁一端为固定边界，另一端受动荷载，边界条件、荷载和网格示意图为

文章首先对上述模型进行动力特性分析，验证刚度矩阵和质量矩阵的计算精度，为后续的计算打基础，与abaqus对比，统计了自振频率为

可以发现，计算精度与abaqus相当，从右图发现sbfem的收敛速度是优于abaqus的。

再对悬臂梁的悬臂端施加动荷载，荷载示意图为

在统计了悬臂梁悬臂端的下部中点的计算计算结果

观测点的加载向位移、速度和加速度数据与abaqus计算保持一致，初步验证了开发的子程序在动力计算中的计算效果。

文章设计了非匹配网格的接触分析算例。采用蒙皮单元技术定义界面接触属性，界面的主从面网格非匹配，与abaqus对比计算结果，模型加载和网格示意图为

统计了界面的应力数据，与abaqus对比的结果为

二者计算一致，初步说明开发的子程序在非匹配界面接触问题中的计算精度与abaqus吻合。

设置以上算例的目标是探究子程序在小自由度问题中的计算精度，为了探究子程序在大规模计算中的稳定性，有必要进行进一步研究。

文章以沥青混凝土的CT模型为例，CT模型的分辨率为200×200×200，共计539170个多面体单元，208880个cohesive单元，909402个节点，给模型顶部施加荷载，计算结果为

可以看到，在用户自定义单元中嵌入的cohesive单元成功参与计算，界面上清晰地产生了损伤。

文章将子程序应用到STL模型中，人体膝盖的STL模型为例

采用八叉树单元离散上述几何模型，模型的单元数目为253079 ，节点数目为361296， 绘制的网格为

在力的作用下，上部骨头向下移动，计算的位移云图为

直到与下部接触，中间的骨头在上下部的挤压下发生变形，中部骨头的应力为

将基于SBFEM的多面体单元作为用户自定义单元在商业有限元软件ABAQUS中实现，实现过程涉及用户子程序UEXTERNALDB和UEL，分别用于存储多面体单元几何拓扑结构和执行多面体单元的计算。详细介绍了存储多面体网格和多面体单元几何拓扑所需的数据结构。结合比例边界有限元理论和abaqus计算框架，对UEL子程序计算流程进行了全面的解释。

此外，提出了一种将abaqus内置单元叠加在SBFEM用户自定义单元上的改进方法来创建基于单元的面，有助于后续采用自定义单元研究界面接触问题。

文章给出了一系列的算例，包括基准测试、动荷载计算、非匹配界面接触计算、基于CT和STL模型的大规模计算，最终的计算结果均与理论解或abaqus保持一致，表明了开发的子程序有良好的计算精度和计算稳定性。