關鍵詞：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，確定出哪些節點在左側，再識別出它們屬于哪些單元的哪個邊。這是因為壓力邊界條件三角形單元三個邊是有區別的。

在生成這個案例網格的時候，也需要編一個小程序，把單元的面積、每個邊長度、每個邊法向等等信息都事先定義好。方便求解器做前處理的時候識別。

很多人在開發求解器做案例測試的時候，前處理都是手動敲單元，麻煩不說，這也說明了程序通用性差，不利于后期集成。

效果

從速度分布結果可以看出，自研線性三角形單元的結果和文獻一致。

自研求解器結果：頂部拖曳+壓力效應的速度分布

文獻結果：頂部拖曳+壓力效應的速度分布

不足之處

根據得到的速度，帶回單元中，按照下式求解壓力。

結果如下，很明顯這個壓力分布完全不對。還是前面的原因造成的：壓力是速度的降階，應該更換單元類型進行壓力處理。

 使用高階單元

改用三角形二次單元作為速度單元，三角形線性單元作為壓力的單元，并且不使用罰函數，速度和壓力耦合求解，結果就對了。

改進后速度、壓力結果

文獻結果