详解有限元法中K、M和C矩阵的集成

零概述分享一篇关于ABAQUS二次开发的文章，这是比例边界有限法与商业软件结合的第一篇文章，来自于CMAME期刊，很有意义，文章配套的有开源程序，方便学习。文章和开源程序的链接为https://www.sciencedirect.com/science/article/pii/S004578252100102X?via%3Dihubhttps://github.com/ShukaiYa/SBFEM-UEL壹sbfem嵌入到abaqus中比例边界有限元法基本理论在这里不做过多介绍，可参考相关文献。这里重点讲解比例边界有限元法如何通过UEL等子程序嵌入到ABAQUS中。有限元求解过程中，单元刚度矩阵的求解和组装是重中之重，这些过程经过商业软件的gui封装，我们已经看不到，只需要点击鼠标就能完成，看起来非常的简单，但这种简单是封装的结果，其内部程序是非常复杂的，abaqus软件主要将fem的算法封装到软件中，但随着有限元理论的发展，各种新型理论层出不穷。作为商业程序，abaqus需要迎合市场需求，因此，abaqus向从业者提供了一系列的二次开发接口，用户可以值只关注自己的科学问题，在编程的过程中能够大幅度的避免重复造轮子，将更多的精力用于优化自己专业部分的理论和程序。其中uel可用于实现刚度矩阵等计算，当abaqus主程序计算单元刚度矩阵的时候，abaqus不再调用自己内部的程序，而是调用用户开发的子程序，子程序执行结束会向主程序返回刚度矩阵、质量矩阵和阻尼矩阵等等，abaqus用于迭代求解。下面是文章给出的带有子程序的abaqus执行流程图sbfem的优势之一是可以直接处理任意多面体的单元，如下面的单元除了编写uel子程序，还需要手动修改inp文件，对于上面的多面体单元，inp文件中的代码为其中，userelement用于声明一个用户自定义单元，nodes指明了单元的节点数目，type相当于给这种单元指明了单元类型，u8中的数字8会以参数的形式传进子程序里面，可用于识别不同的单元类型，properties指的是单元的浮点型材料属性个数，coordinates代表该单元激活的自由度数目，variables是单元的状态变量个数，与单元的节点数目有关，详细的计算方法可查看文章正文。element用于声明单元的编号和单元节点信息，需要特别注意的是，后面的节点排序是非常有讲究的，我在另一篇帖子详细写过，具体可查看https://mp.weixin.qq.com/s?__biz=MzkwNTc1MDM1MA==&mid=2247483683&idx=1&sn=0750ac09ccc3e840b53be8aaf6213deb&chksm=c0f3b302f7843a14a3253499b92793560c8c467bdd762a7b78c58384fc110f061a53594ccfbf#rd以六面体单元为例，在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模型为例，CT模型的分辨率为200×200×200，共计539170个多面体单元，208880个cohesive单元，909402个节点，给模型顶部施加荷载，计算结果为可以看到，在用户自定义单元中嵌入的cohesive单元成功参与计算，界面上清晰地产生了损伤。陆人体膝盖接触分析算例文章将子程序应用到STL模型中，人体膝盖的STL模型为例采用八叉树单元离散上述几何模型，模型的单元数目为253079，节点数目为361296，绘制的网格为在力的作用下，上部骨头向下移动，计算的位移云图为直到与下部接触，中间的骨头在上下部的挤压下发生变形，中部骨头的应力为柒结论将基于SBFEM的多面体单元作为用户自定义单元在商业有限元软件ABAQUS中实现，实现过程涉及用户子程序UEXTERNALDB和UEL，分别用于存储多面体单元几何拓扑结构和执行多面体单元的计算。详细介绍了存储多面体网格和多面体单元几何拓扑所需的数据结构。结合比例边界有限元理论和abaqus计算框架，对UEL子程序计算流程进行了全面的解释。此外，提出了一种将abaqus内置单元叠加在SBFEM用户自定义单元上的改进方法来创建基于单元的面，有助于后续采用自定义单元研究界面接触问题。文章给出了一系列的算例，包括基准测试、动荷载计算、非匹配界面接触计算、基于CT和STL模型的大规模计算，最终的计算结果均与理论解或abaqus保持一致，表明了开发的子程序有良好的计算精度和计算稳定性。点击卡片关注我们来源：有限元先生