有限體積法下的控制方程

02-04

看著題目，肯定自忖又是老生常談的東西，事實上並非如此。

如果大家文獻看的多了，肯定能夠發現一個有趣的現象——不同文獻的有限體積控制方程有的相似有的差別較大。

這裡我選出兩種有代表性的例子，以後的講解也基於這兩種形式的控制方程。

第一種：

$nfrac{partial }{partial t} int_{Omega }^{} bar{W} dOmega + oint_{partial Omega }^{} [(bar{bar{F_{c} } } - bar{bar{F_{v} } })cdot bar{n} ]dS=int_{Omega }^{} bar{Q } dOmega$

書籍：《COMPUTATIONAL FLUID DYNAMICS: Principles and Applications》

這種形式方程文獻中見的最多，早期的大牛paper中方程都長這樣，這本書也很有名，開計算流體力學課的高校應該多用這本書。

第二種：

$nfrac{partial rho phi }{partial t} + nablacdot (rho bar{v} phi ) =nablacdot (Gamma ^{phi } nablaphi )+Q^{phi }$

$nfrac{partial }{partial t} int_{Vc }^{} {rho phi } dV + oint_{partial S}^{} J^{phi,C } cdot bar{n} dS= oint_{partial S}^{} J^{phi,D} cdot bar{n} dS+int_{Vc}^{} Q^{phi } dV$

書籍：《The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM? and Matlab?》

這種形式的方程不能小覷，因為OpenFOAM用的就是這個形式。

我個人更喜歡第二種形式中的微分方程，非常優美，雖然其割裂了某些東西（下面會講到）。

首先明確一點，上述兩種形式是數學上是殊途同歸的，但在FVM中有些許不同。目前暫且講不到具體的離散過程，但是可以初步比較一下個中差別。

第一種方程初期是用於處理可壓流的，所以其變數寫成守恆形式，且方程用統一格式表達，因為可壓流流場變數全部耦合，每個變數的求解是耦合式求解（每個點流場變數同時求出）。第二種方程初期是用於不可壓流的（OpenFOAM早期只針對不可壓流），所以喜歡寫成微分方程形式，雖然也用統一格式表達，但是方程並非全部耦合，所以常常用分離式方法（和耦合式相反）求解流場變數。
兩種方式中定義的通量不同。兩個區別，第一個很淺顯，第一種形式為通量沿界面法向，第二種形式通量沿速度方向。第二個區別不太直觀——兩種方式通量的含義不同。

首先是對流通量，第一種形式中，壓力 $p$ 放入到了 $bar{bar{F_{c} } }$ 中，而第二種形式則純為對流通量。原因也在於第一種形式起初是針對可壓流，所有特徵速度都是有限值，壓力波也不例外，也有迎風性。求解對流通量，第一種形式中常常出現 JST和Roe這些格式，而第二種沒有這些東西，原因也在於可壓縮和不可壓縮。
其次是第二個通量，為什麼我叫它第二個通量，因為他們名字確實不一樣。第一種形式中 $bar{bar{F_{v} } }$ 叫粘性通量，第二種形式中 $J^{phi ,D}$ 叫做擴散通量。其最根本的區別在於粘性項的處理，粘性力 $tau _{ij} =mu [nablabar{v} +(nablabar{v} )^{T} ]+lambda (nablacdot bar{v} )bar{bar{I} }$ 。第一種形式中這些項全部在粘性通量中，而第二種形式對比可以發現，粘性力只有第一項在擴散通量中（剩下的在源項中）。

兩種形式的最終離散目標都是化為： $Ax=b$ ，問題在於怎麼獲得 $A ,b$ 。

首先必須明確，現今採用的有限體積法都是時空分別離散（不耦合），英文名：method of lines。

第一種形式——整體考慮：

定義殘差（residual）： $nR = oint_{partial Omega }^{} [(bar{bar{F_{c} } } -bar{bar{F_{v} } })cdot bar{n} ]dS -int_{Omega }^{} bar{Q } dOmega$ 。

原方程化為： $frac{dOmega W}{dt} =-R$ ，可以看到這裡用了常微分方程的符號，method of lines的原理就是把空間先離散把偏微分方程化為常微分方程，因為針對常微分方程的數值方法很多。可以看到 $W,R$ 我沒有標上標，因為不同的上標代表求解方式是隱式還是顯式的。

顯式（explicit）： $frac{Omega }{ Delta t}Delta W^{n} =-R^{n}$ ， $Delta W=W^{n+1} -W^{n}$ 。用上一時間步計算下一時間步，每個點和其他點不需聯立求解，也即 $A$ 為對角陣。

隱式（implicit）： $frac{Omega }{ Delta t}Delta W^{n} =-R^{n+1}$

注意，關鍵來了，這裡要用到一個必須採用的技巧，第二種方式依舊使用的——線化。

泰勒展開保留一階導數項， $R^{n+1} = R^{n}+frac{partial R}{partial W} Delta W^{n}$ ，上式化為： $[frac{Omega }{ Delta t} +frac{partial R}{partial W} ]Delta W^{n} =-R^{n}$ ，方括弧內的為 $A$ ，並非對角陣，流場中離散點需要耦合求解。其中 $R$ 對 $W$ 的jacobi矩陣的計算構成了有限體積法的一個難點，其對隱式格式時間步長提出了限制（隱式格式理論上無時間步長限制）。

第二種形式——每項分別考慮：

每一項上（比如對流項、擴散項）分別能夠產生部分的 $A$ 和 $b$ ，把它們加起來獲得最終的 $A$ 和 $b$ 。了解過OpenFOAM朋友應該知道，OpenFOAM中有兩個名稱空間fvm和fvc，這兩個東西作用還不太相同，採用fvm作用某一項叫隱式處理，為什麼？因為其產生 $A$ 和 $b$ ，處理思路和第一種形式方程隱式處理基本一致，也是採用線化，一部分跑到 $A$ 中，一部分到 $b$ 中。採用fvc則叫顯式處理，它不產生 $A$ 。例子我就不舉了，有興趣可以看我提供的那本書。

-------------------------------------------------------------------------------------------------

這篇寫得有點長，之後的文章不定期更新。

我有兩篇回答，對初學者可能有些幫助，推薦看一下，自學是最重要的。

有哪些適合cfd初學者練習的題目？n以及用什麼工具求解這些題目？ - HangZS 的回答

空氣動力學/流體力學的big picture和自學指南？ - HangZS 的回答

本來我準備在我的個人網站（CFD之旅）上寫文章，奈何主機不行，半天反應不過來急死人，知乎這個平台確實不錯，以後就在這上發布了。