小白最近對流體計算的收斂產生了困惑。以前在學習高等數學的時候,小白接觸過了級數的收斂,由於當時貪玩,並未將其放在心上,因此大學結束了小白也只是記住有這么一個名詞罷了。現如今在利用CFD的過程中,小白又一次碰到了“收斂”這一名詞。小白找了很多的資料,然而資料中關於收斂的介紹,無一不是一大堆的數學公式,小白覺得頭很疼。
“出來混,總是要還的。”小白的心情很復雜。
“流體計算為什么要收斂?收斂是什么意思?不收斂又有什么后果?如何判斷是會否收斂?如果不收斂該采取何種措施使其收斂?”小白帶着滿腔的疑惑,找到了小牛師兄。
“師兄,最近在看資料的時候,老是碰到收斂這個名詞,我想問一下,什么叫收斂呢?我最近看了很多的資料,都是用一堆數學公式來表達,看了很久概念也很模糊。”小白問道。
“也好,既然你注意到了收斂,今天我們就來聊一聊這個話題。”小牛師兄說。
計算機解決工程問題
“在談收斂之前,我們還是先來談一談數值方法。關於數值計算,我想你也了解了不少。你說計算機解決這些工程問題到底是如何實現的呢?”小牛師兄問道。
“我看了一些資料。利用計算機來解決工程問題,其核心是利用計算機求解物理模型,或者說是利用計算機求解數學方程。”小白說。
“基本上是這樣。我們知道描述現實世界物理現象的數學模型大部分都是微分方程,利用計算機並不能直接求解微分方程,計算機通常需要利用數學方法將這些微分方程轉化為代數方程,通過求解代數方程獲取原微分方程的解。那么問題來了,計算機是如何將微分方程轉化為代數方程的呢?”小牛師兄說。
“轉化的方法其實有很多種,比如我們常說的有限差分法、有限元法、有限體積法等都是用來干這種工作的。為了實現這種軟件,引入了網格。采用網格的目的是將連續的空間和時間轉化為離散的空間與時間,從而在每一個網格單元或網格節點上獲得代數方程。所有網格上的代數方程集合在一起,就構成了整個計算域上的代數方程組,求解方程組就可以得到每一個網格單元或網格節點上的物理量。”
“當然,求解方程組得到的物理量還是離散的,為了得到節點之間或單元之間的物理量分布,采用了一種稱之為插值的數學方法。利用插值可以得到近似連續的物理量分布。當計算網格或計算時間間隔足夠小的時候,通過插值可以得到相當精確的計算結果。”小牛師兄說。
“因此我們可以說,利用計算機求解工程問題,實際上歸根結底是利用計算機求解工程問題所涉及的物理模型轉化而成的代數方程。那么問題來了,利用計算機求解代數方程是怎么實現的呢?”小牛師兄說。
“我記得線
性代數里面有講到,可以采用矩陣求逆、高斯消去法、三角分解法、克拉默法則等。”小白說。
“嗯,這些方法的確可以用於計算機求解線性方程組,這些方法可以統稱為直接求解法。然而在工程應用上是存在很大問題的。我們知道,為了得到比較精確的計算結果,網格數量常常非常多,直接后果就是最終的代數方程組包含的方程數量非常非常多。打個比方,如果網格數量是100萬個,那么最終求解的代數方程組就有100萬個方程,在利用直接法求解的過程中,需要采用100萬維的矩陣進行存儲,如此龐大的矩陣,消耗的內存是非常巨大的,而且其求逆計算或三角分解計算都是非常困難的,計算效率是極其低下的。因此,工程上求解這類數量龐大的方程組通常采用迭代法。”小牛師兄說。
“迭代法其實也有很多種,比較常見的如雅克比迭代、高斯-賽德爾迭代、超松弛迭代等等。關於迭代法求解線性方程組的具體方法,可以找一本數值計算的數看看。收斂就是迭代法求解線性方程組中的一個概念。”小牛師兄最后說。
關於迭代殘差
為了更好的說明收斂的問題,這里先舉個例子來說明什么叫殘差。
比如說要求解下面的線性方程組:
可將(1)式改寫為迭代計算的形式:

給定一個初值,如x(0)=(0,0,0)T,代入到式(2)中,可以得到x(1)=(2.5,3,3)T
反復迭代,

可以得到一系列x1,x2及x3的值。
迭代次數 | x1 | x2 | x3 |
---|---|---|---|
0 | 0 | 0 | 0 |
1 | 2.5 | 3 | 3 |
2 | 2.875 | 2.363636364 | 1 |
3 | 3.136363636 | 2.045454545 | 0.971590909 |
4 | 3.024147727 | 1.947830579 | 0.920454545 |
5 | 3.000322831 | 1.983987603 | 1.000968492 |
6 | 2.993753228 | 1.999970652 | 1.003841684 |
7 | 2.999028573 | 2.002620797 | 1.003130723 |
8 | 3.000200118 | 2.000637857 | 0.999830514 |
9 | 3.000281568 | 1.999911822 | 0.999740477 |
10 | 3.000031814 | 1.999874019 | 0.999881261 |
11 | 2.999982442 | 1.999977637 | 1.000015588 |
12 | 2.999987717 | 2.000007802 | 1.00001437 |
13 | 2.999999333 | 2.000005773 | 1.000004191 |
14 | 3.000001117 | 2.000000623 | 0.99999889 |
15 | 3.000000511 | 1.999999493 | 0.999999286 |
16 | 2.999999988 | 1.999999749 | 0.999999871 |
17 | 2.999999938 | 1.999999992 | 1.000000068 |
18 | 2.99999998 | 2.000000029 | 1.000000033 |
19 | 3.000000003 | 2.00000001 | 1.000000003 |
20 | 3.000000003 | 1.999999999 | 0.999999996 |
可以看出,隨着迭代計算的進行,計算結果越來越接近解析解[3,2,1]。
我們可以將某一物理量兩次迭代計算值的差值稱之為該物理量的迭代殘差。比如上面的迭代計算x1的第一次迭代計算殘差為2.5,第二次迭代計算后殘差變為0.375。
當然上面提到的殘差為絕對殘差,在實際計算過程中,有時也常常采用相對殘差,即物理量隨迭代進行的變化量。可以看到,當殘差非常小時,可以認為不用再往下計算了,或者說再往下算已經沒有太大意義了。
CFD計算中的殘差
在CFD計算中,每一個網格上都會存儲眾多物理量,因此每一個網格上的任一個物理量在計算迭代過程中都會存在一個殘差,這意味着在一次迭代過程中,同一物理量在不同的計算網格上有不同的計算殘差,而實際上我們在進行CFD計算時,每一個迭代步只對應着一個殘差值。
CFD中殘差分為幾種:
最大殘差:在一次迭代中,取所有網格中的殘差值的最大值作為本次計算的殘差。
平均殘差:在一次迭代中,計算所有網格中的算術平均值作為本次迭代的殘差
均方根殘差:在一次計算,計算所有網格中殘差值的均方根作為本次迭代計算的殘差
在CFD計算中,常常采用均方根殘差(RMS)作為殘差值。
CFD中的收斂
所謂迭代收斂,簡單來說,就是在迭代計算過程中,物理量趨於某一值的情況。CFD計算中判斷收斂通常有三種方法。
收斂判斷規則之一:殘差達到某一設定標准時可以認為迭代計算達到收斂。
這條規則最簡單,在實際工程應用中也最常用。通常設定某一標准,當迭代計算過程中殘差值低於此標准時則認為計算收斂。這也是幾乎所有CFD軟件用於判斷收斂的基本方法。
然而此規則在實際工程應用中常常無法使用,有些復雜的問題,無論你怎么計算,其殘差也不會下降,甚至有時候殘差會出現周期性震盪。
殘差穩定在某一位置不下降的原因有很多,常見的原因包括:
- 計算區域中存在低質量的網格。低質量的網格會造成計算殘差增大及殘差震盪
- 邊界條件設置有誤。錯誤的邊界條件或邊界類型搭配都會導致計算殘差震盪。
- 利用穩態求解器計算瞬態問題也會造成殘差的震盪。
在殘差無法達到設定標准的時候該如何判斷收斂呢?需要注意的是,此時的收斂並非數學意義上的收斂了,而只是意味着我們可以停止計算。
收斂判斷規則之二:進出口物理量通量達到平衡
最常用的是判斷進出口質量是否相等,這實際上是判斷連續性條件是否達到滿足的。實際上還有很多,如計算域中包含化學反應時,判斷進出口組分是否守恆;如計算域中包含多相流時,判斷進出口各相質量是否守恆等。此規則是一種非常弱的規則,實際上只是收斂的一個必要條件而已。但是在第一條判斷規則無法達到時,也常常采用此規則來判斷。
收斂判斷規則三:計算域中的物理量隨迭代進行不再發生變化
這條規則在實際應用中也很常用,甚至比第二條規則更常用。在實際工程中,經常監測某些敏感位置的物理量,當隨着迭代進行,監測的量不再發生變化時,基本可以認為計算達到收斂。
注意后兩種方法只是在殘差無法達到標准時才采用的判斷方法,並不意味着計算就收斂了。
穩態及瞬態計算的收斂
在CFD迭代計算過程中,穩態殘差曲線與瞬態殘差曲線的形態有很大的不同。
如下圖所示為穩態計算殘差,每一個物理量都有一條殘差曲線,當殘差曲線低於設定的標准時認為計算收斂,當所有的殘差曲線低於設定的標准時,計算結束。
下圖所示為典型的瞬態計算殘差。可以看到其形狀特征與穩態殘差曲線有很大的不同。瞬態計算要求在每一個計算時間步內達到收斂。若不能達到收斂,則迭代次數達到所設定的內迭代次數后進入下一個時間步,重新開啟迭代計算。因此書台計算殘差呈現出下圖所示的鋸齒狀。
注意:瞬態計算要求每一個時間步內均達到收斂。