写出一个有限元仿真求解器的方法和具体工作

详细信息

作者优秀: 优秀教师/意见领袖/博士学历/特邀专家/独家讲师
平台推荐: 内容稀缺

2年前浏览2888

导读：本文来自803教研室主编、仿真秀专栏作者刘畅老师的原创文章《如何写出第一个有限元求解器》，刘老师毕业于南京航空航天大学，主要从事工业软件开发、飞行器结冰设计、先进复合材料结构设计、图像识别等领域研究。真心希望能够引发正在学习进阶二次开发、国产软件和CAE软件程序开发等专业技能的理工科学子和研发工程师的共鸣，起到抛砖引玉目的，如有不当欢迎批评指正和共同进步。也强烈推荐大家关注仿真秀平台的CAE软件程序开发公开课，详情见后文。

一、写在前面

仿真求解器距离可用的仿真软件之间，差异很大。仿真软件的开发是一场完整的战役，开发者作为“统帅”，必须统筹考虑软件逻辑、通用的前后处理、可能的bug、用户体验等方方面面。相比之下，单纯的求解器开发更像一场战斗，百万军中取敌酋首级，干就完了。

笔者以前学有限元的时候，写过简单的二维杆、二维梁、平面应力单元的有限元程序。后来研究需要，还写过一个梁单元的屈曲特征值程序。再后来在做二维结冰算法的时候，刚开始采用有限元法计算流场，所以也做了一段时间的计算流体力学有限元程序。

所以对我来说，这并不能算是我严格意义上的第一个有限元求解器。但是，从实用性的角度、程序通用性的层面，这次做的有限元程序要远高于以往。

二、想清楚你要什么样的求解器

我这次想做的是一个针对四面体单元、线弹性、静力学的有限元求解器。为了便于后面描述，我们给它取个名字：TetSolver。为了增强趣味性，还可以给它设计个LOGO。几乎所有复杂的结构都可以用四面体单元进行离散，因此从应用覆盖面来说，TetSolver是可以的。

除了上述考量，基于通用性的考虑，我们还要搞清楚TetSolver的输入和输出。参考ABAQUS软件，比较方便的输入，就是定义统一的inp文件，把材料参数、节点坐标、单元信息、边界条件、载荷都放在一个文本中，供程序调用，然后直接进行计算。

为了尽快开展实际工作，减少后期程序文本处理的时间，我们分别定义五个txt文本用于存储上述信息。

另一个方面，为了保证通用性，还要定义好每个输入文本的格式，并且在程序的开头就把文本都读入进来，统一命名，后面的所有程序都是基于这些命名开展，这样即便更换输入模型，也能开展计算。

凡是设计到对输入文本进行参数、大小的输入，均自动化。说了这么多，一言以蔽之，即便更换输入模型，我们也仅仅需要输入五个txt文本的名称就能自动计算。

现在，可以定义TetSolver的初步框架了：

接下来我们对TetSolver计算功能做一个规划，这是第一个四面体求解器，为了减少开发难度，我们规定：

（1）位移边界条件均为固支；
（2）载荷为力载荷；
（3）暂时不开展网格离散工具和后处理工具的开发；
（4）处理线弹性静力学问题。

最后一个很重要的方面是，模块化编程。我们把大大小小各类计算、输入读取、处理等等，尽可能写成一个个自定义函数，而在主程序中近保留必要的流程顺序代码。这样做，处于三个考虑：

增强程序的可移植性；
便于后续改进；
便于程序调试。

到这里，基本上我们把TetSolver规划的差不多了，总结下，TetSolver 1.0版本的概念设计中，我们考虑了开发难度、通用性、可移植性、应用范围等因素。

三、选择一个合适的程序语言

从个人使用习惯来说，我更倾向于用MATLAB来做。并且，有限元程序中主要的工作是矩阵运算，MATLAB处理起来更有优势，编程工作量相对小。

最后出于练习Python的考虑，我这次决定用Python来做。

四、准备工作

一些准备工作是必须的。首先就是找到合适的参考资料。目前关于有限元入门编程的书籍很多，不少的书籍上还会附上代码，这些代码可以帮助我们快速开展编程工作。我这次主要参考的资料有：

《MATLAB有限元分析与应用》，这是翻译过来的书籍，作者Peter I. Kattan，这本书是一本很棒的有限元入门教材，从简单的弹簧元，到梁、平面、四面体单元均有深入的讲解，并且提供了源代码。
《Python与有限元》，作者是湖北工业大学讲师裴尧尧博士，这本书不局限于讲授有限元程序，而是从有限元框架讲起来，很有利于做更深入的有限元编程研究。

当我们手里有了很好的参考文献的时候，开始会满怀期待，以为很容易就能把TetSolver开发出来。实际上，教材处于方便讲解的考虑造成代码不连贯、印刷错误、书本往实际应用转换理解错误等等，只要你想按照自己的框架开发有限元程序，都会遇到各种各样的困难，所有不可过分依赖参考资料，要独立思考自己遇到的问题。

五、编程与验证

终于到了编程这一步。按照既定的策略，我开始如下流程的编写：

1、单元刚度矩阵代码的编写。这个直接参考《MATLAB有限元分析与应用》；

2、总体刚度矩阵组装代码的编写。《MATLAB有限元分析与应用》在这一部分，采用的是笨方法，一个元素一个元素的定义，这样不利于我们后面引入其他单元。所以这里参考《Python与有限元》给出的思路，进行自动组装程序的编写。

3、根据位移边界条件约束总体刚度矩阵。这里采用置大数法，参考《飞行器结构力学》；

4、载荷列阵生成。

5、文本输入方式涉及的材料参数读入代码、节点/单元/边界条件/载荷读入调用代码；

6、输出代码编写；

7、梳理代码，搭建框架：

INPUT，输入；
LOAD，加载；
SOLVE，求解
OUTPUT；输出。

验证是个至关重要的一步。对于写大型程序来说，一定避免一次性验证，这样很难定位问题。TetSolver采用的策略是：

（1）每一步都验证，每个自定义函数或者关键模块都配一个对应的test程序。完善一个才进行下一个。

（2）找到一个信息完整的算例。《MATLAB有限元分析与应用》给出一个悬臂梁案例，每个单元都给出来刚度矩阵结果、总体刚度矩阵结果。比较有利于我们对比、查找错误。

（3）多尝试几个验证案例。这个便于了解我们代码中bug，以及熟悉我们自己编写的求解器的特点。

这里贴出我们最后成稿的TetSolver主程序框架：

为了便于查看计算结果，增强信心，我把之前用MATLAB写的云图处理程序拿过来做了可视化：

《MATLAB有限元分析与应用》中的案例复现

自验证的三点弯曲案例

六、写在最后

好啦，至此我们完成了TetSolver 1.0的开发。这个过程中遇到的比较大的问题主要有三个：

1、Python的语法。尽管之前学过一段时间，这几天用它写程序还是非常不习惯。特别是矩阵处理很容易出错，list、float、string。

2、总体刚度矩阵的集成。刚度矩阵集成本来就很绕，再用Python写，简直灾难，特别是Python数组初始位为0，极容易出错。

3、个人不细心造成的问题。这个就太多了，可以说时刻都在发生。

如果要进一步研究开发，实际上我们有两条路可以走：

横向。继续做其他方向的初级求解器，比如屈曲、模态分析等等。
纵向。在静力学这个方向上，进一步开发非线性求解器，并谋求在某个特定领域的应用，作为解决某个特殊问题的专用工具。积累一定经验后，再开始往支持各向异性、支持多种本构、支持多种单元这个方向走。然后再往动力学、传热方面走。当然条路很难也很远。

具体怎么走？要看个人兴趣，以及能力基础。

如果你的是对标ABAQUS、NASTRAN，那就是要打正规战和阵地战，最好背靠官方组织。大牛牵头，团队作战，做科学计算的、静力学的、动力学的、热学的、电磁学的、声学的、隐式的、显式的、做网格的、做后处理的、做本构的、做接触的、做交互的、做验证的、做集成的，拧成一股绳，十年一剑。

时代变了，尽管人们津津乐道于两三个年轻人创办ABAQUS这样的故事，但是实际上在ABAQUS核心能力形成的关键阶段，已经是有大量的经费、成熟的团队了。NASTRAN则一开始就是官方组织的。现阶段要想三两人做出目前这个状态的ABAQUS，几乎是不可能的。

如果你的目标是先做一个求解器，先用起来解决某一问题，形成一个特色。那么最好农村包围城市，需求牵引着，先打下一块根据地，活下来，再慢慢壮大。积少成多，由易到难。

（完）

作者简介：刘畅，803教研室主编，仿真秀专栏作者，毕业于南京航空航天大学，主要从事工业软件开发、飞行器结冰设计、先进复合材料结构设计、图像识别等领域研究。

声明：原创作品，部分图片和内容源自网络，如有不当请联系我们，欢迎分享，禁止私自转载，转载请联系我们。

来源：仿真秀App

Fluent 静力学复合材料非线性二次开发通用航空航天汽车 MATLAB 芯片声学材料

著作权归作者所有，欢迎分享，未经许可，不得转载

首次发布时间：2022-08-24