流体有限元求解器开发-二维斯托克斯方程 原创

关键词：CFD，有限元，三角形单元，罚函数，粘性流动

最近工作室有流体有限元求解器的开发需求，我在前面讲飞机结冰的文章提到过，差不多10年前瞎捣鼓过这个东西。

好多东西都记不清了，先从一些简单的流动模型入手，做一些恢复性训练。考虑到我是结构力学出身，在进行流体有限元开发的时候，我会代入结构有限元的视角进行分析。

流体也好，固体也好，CFD也好，FEM也好，有很多开源工具、源代码可以用。但是如果想要开发自己的软件，并期待它未来解决工程的各种复杂艰深问题，那么从理论公式入手的推导，每一个细节开始自己编程，是非常重要的。

吸星大法得到的功力只能应付一时之急，时间长了后患无穷。尤其在AI能力越发强大的今天，在科学与工程问题上，如果没有足够的理论基础和工程经验，完全靠AI堆积的代码看似庞大，实则是无根之萍。

初中的时候，我们数学老师逼着我们自学。到高中的时候，我的自学能力已经颇为可观，基本暑假一个月左右的时间，我就把一本数学必修和一本物理必修学完了，还各自刷完一本练习册，做到通关的程度。当然自学经常遇到的问题是，很多曲折复杂的思路一时难以完全搞懂，需要花很多时间去想去练。当时有同学跟我说：你这自己学太累了，你花了几天时间搞懂的东西，老师课上讲一遍我们就听懂了，还是听课的效率高。

别人直接喂的东西，当然效率高，但是这玩意吃了只会长肥肉，它不长力气。事实也证明，老师一讲他就懂，一到考试就全不会。

控制方程

二维斯托克斯方程的控制方程装X写法如下：

朴实点的写法是这样：

和NS方程比，它没有对流项，适用于高粘度流体、低速流动等场景。

有限元思路

搞结构力学有限元和其他方向有限元最大的区别是：结构力学有限元发展的太成熟了，杆梁板壳，各种模型的刚度矩阵前辈都给你推导好了。我在开发结构力学有限元求解器的时候，都是先去查资料，直接就把单元刚度矩阵拿过来用。

但是到流体这就不能完全这么干，原因是：

（1） 未必能找到直接可用的单元矩阵；

（2） 流体的边界条件中，有很多第二类边界条件（压力、热流、对流），这些边界条件有的是放到右侧载荷项，需要推导。有的载荷甚至会影响左侧单元矩阵，部分放到左边去修正。结构力学里面大部分情况下，载荷就是在右侧，位移就是放到左边修正刚度矩阵，清晰明确。

所以，开发流体求解器的时候，还是要从有限元的基本方法入手。这里采用加权余量法进行处理。有限元的教材里面讲的很多了，这里简单说一下流程：

（1） 根据单元类型，确定插值函数。此时速度、压力等变量，都可以用权函数表达。

（2） 采用伽辽金方法，权函数=插值函数，控制方程与权函数相乘，积分取0。

（3） 在每个单元域内，方程转换为权函数的积分形式，最终形成单元矩阵。

单元方程

最终得到单元的方程形式如下：

类比到结构有限元，左侧第一大项，就是刚度矩阵。u、v、p相当于3个自由度，右侧就是载荷列阵。

我的本意是开发一个三角形单元的斯托克斯求解器。使用三角形线性单元对应的插值函数：

无论是CFD还是结构有限元，只要单元类型一致，插值函数都是一样的，区别只在单元方程。

但是这样直接求解结果是震荡的，原因是结构有限元中，节点的几个自由度本质是同一种物理量，它们是“平级的”。但是速度u、v与压力p不是同一种物理量，它们不平级。压力和速度的降阶是平级的。

一般的处理思路是，对速度用高阶单元，对压力用低阶单元。还有一种思路是使用罚函数方法，将压力用如下形式表达，只要λ足够大即可。这样就在原控制方程中消除了压力。

得到下式，这种方程的求解和结构有限元就完全一样了，基本等效成了每个节点两自由度的结构有限元。

约束与总体方程装配

对于指定的速度条件，它等同于结构有限元的位移约束，采用乘大数法，对左侧矩阵对角线元素进行处理，同时对右侧对应的载荷项进行处理。

对于指定压力边界条件，这个类似于结构力学里面的将边载荷移植到节点上的处理。

总体方程的装配，和有限元中刚度矩阵的装配方法完全一致。

编程方法

采用Python进行编程实现。只要完成理论推导，单元方程、总体方程矩阵组集都是很简单的。

比较麻烦的是边界条件的处理。

CFD中，域的概念很强，常见的inlet、outlet、wall等等条件都不是夹在单点上的，对于二维问题，一般都是加载在边上的。

以如下矩形域算例为例（粘度取1000Pa·s)：

用户在定义边界条件的时候，以左侧边为例，最方便的表达是：

{“x_min” :{“P” :1000 } }

到我们做前处理的时候，就需要根据x_min，确定出哪些节点在左侧，再识别出它们属于哪些单元的哪个边。这是因为压力边界条件三角形单元三个边是有区别的。

在生成这个案例网格的时候，也需要编一个小程序，把单元的面积、每个边长度、每个边法向等等信息都事先定义好。方便求解器做前处理的时候识别。

很多人在开发求解器做案例测试的时候，前处理都是手动敲单元，麻烦不说，这也说明了程序通用性差，不利于后期集成。

效果

从速度分布结果可以看出，自研线性三角形单元的结果和文献一致。

自研求解器结果：顶部拖曳+压力效应的速度分布

文献结果：顶部拖曳+压力效应的速度分布

不足之处

根据得到的速度，带回单元中，按照下式求解压力。

结果如下，很明显这个压力分布完全不对。还是前面的原因造成的：压力是速度的降阶，应该更换单元类型进行压力处理。

 使用高阶单元

改用三角形二次单元作为速度单元，三角形线性单元作为压力的单元，并且不使用罚函数，速度和压力耦合求解，结果就对了。

改进后速度、压力结果

文献结果