計算流體力學

流體力學和其他學科一樣，是通過理論分析和實驗研究兩種手段發展起來的。很早就已有理論流體力學和實驗液體力學兩大分支。理論分析是用數學方法求出問題的定量結果。但能用這種方法求出結果的問題畢竟是少數，計算流體力學正是為彌補分析方法的不足而發展起來的。

早在20世紀初，理查德就已提出用數值方法來解流體力學問題的思想。但是由於這種問題本身的複雜性和當時計算工具的落後，這一思想並未引起人們重視。自從40年代中期電子計算機問世以來，用電子計算機進行數值模擬和計算才成為現實。1963年美國的F.H.哈洛和J.E.弗羅姆用當時的IBM7090計算機，成功地解決了二維長方形柱體的繞流問題並給出尾流渦街的形成和演變過程，受到普遍重視。1965年，哈洛和弗羅姆發表“流體動力學的計算機實驗”一文，對計算機在流體力學中的巨大作用作了引人注目的介紹。從此，人們把60年代中期看成是計算流體力學興起的標誌。

計算流體力學的歷史雖然不長，但已廣泛深入到流體力學的各個領域，相應地也形成了各種不同的數值解法。就目前情況看，主要是有限差分方法和有限元法。有限差分方法在流體力學中已得到廣泛應用。而有限元法是從求解固體力學問題發展起來的。近年來在處理低速流體問題中，已有相當多的應用，而且還在迅速發展中。

為了説明計算流體力學主要方法，需先了解流體力學運動的基本方程的性質和分類。流體力學的基本方程是在19世紀上半葉由C.-L.-M.-H.納維和G.G.斯托克斯等人建立的，稱為納維-斯托克斯方程，簡稱N-S方程 ^[1]  ，二維非定常不可壓縮流體的N-S方程為：

式中u、v為沿着x、y方向上的速度分量；t為時間；p為壓力；ρ為密度；ν為運動粘性係數。在不同條件下，N-S方程的數學性質也不一樣。

①N-S方程描述粘性流體隨時間而變的非定常運動。時間項和方程右邊的高階導數項決定方程的性質。它同二維熱傳導方程類似，屬於拋物型方程。

②粘性流體的定常運動是將原方程中的時間項省去。此時N-S方程的性質，取決於它的高階導數項，和拉普拉斯方程一樣，為橢圓型方程。

③無粘流的歐拉方程是將N-S方程的右邊粘性項略去而得。它也適用於可壓縮流體。從形式上不容易判斷歐拉方程的性質。因多數無粘流動皆為無旋流動，故如將歐拉方程改用速度勢ψ表示，則二維定常可壓縮氣流的方程為：

式中c為聲速。此式是二階偏微分方程

的一般形式，其性質要看B²-AC

0而定。在超聲速區，B²-AC

0，即

，上式類似於波動方程，為雙曲型；在亞聲速區，B²-AC

0，即

，上式便與拉普拉斯方程相同，為橢圓型。總之，流體力學的運動方程是極其複雜的非線性偏微分方程，具有各種不同的類型，而且往往還是混合型的。要全面描述流體的運動，還必須同時考慮其他方程，如連續性方程、能量方程和狀態方程等。所以計算流體力學在很大程度上就是針對不同性質的偏微分方程採用和發展相應的數值解方法。

低速無粘流動數值解 在無旋條件下，低速流動的速度勢滿足拉普拉斯方程或泊松方程。很多平面問題利用複變函數和保角映射可以求得解析解，這是經典流體力學的重要內容。但對幾何形狀比較複雜的物體，必須用下述的數值解法。

這是用逐步近似求解聯立方程的方法，也是橢圓型微分方程的主要數值解法。此法程序簡單，存儲量與運算量均比較小，一般先假定一組初值，然後求每個網點上的新值。以五點格式為例，網點上的新值是鄰近四點初值的平均。新值求出後，舊值還要保留，以便計算其他各點的新值。這種簡單迭代收斂很慢，現已很少使用。但若稍加改進，用算出的新值沖掉舊值，並引進一個松馳因子，以加速收斂，將算出來的新值與原來的舊值加權平均，就成為50年代發展起來的逐次超松馳法。

這是用非定常方程求解定常問題的方法，常用於求解N-S方程和歐拉方程等。雖然用的是非定常方程，但所解的並不是非定常問題。根據給定的初始條件以及隨時間改變的約束條件，非定常問題是研究流動隨時間的演變過程。這種非定常行為和給出的初值很有關係。然而時間相關法的初值，原則上是隨意選取的，只是須滿足定常問題所規定的邊界條件。在求解過程中，流動隨時間的變化並不代表真實的物理過程。當時間足夠長後，未知函數值逐步與時間無關，便漸近趨於定常解。所以時間相關法實際上也是一種迭代法，時間變量只不過是用來記錄迭代的次數而已。

流體力學的應用問題，往往是二維和三維的空間問題。由於穩定性的要求，時間步長受維數的限制，維數愈高，要求時間步長愈小，計算工作量也愈大。50年代中期D.W.裴斯曼和J.道格拉斯等人提出所謂交替方向隱式法，以加快計算速度。如在二維非定常方程中，第一步先對x的導數用隱式差分，而y方向的導數則用前一個的數值。第二步對y的導數用隱式差分，x方向的導數則用第一步算出來的數值。這一方法的優點是穩定性好，有足夠的二階精度，所產生的差分方程是三對角矩陣方程，便於求解。

解位勢流動的一種數值方法。航空工業中的低速飛機設計採用位勢理論計算各種氣動力參數，就是求解二維或三維拉普拉斯方程。在經典流體力學中，用基本解的疊加來解拉普拉斯方程的做法是很成功的。這種方法的要點是，用源、匯、偶極子的分佈代替機翼和機身對流場的影響。它們的強度由邊界條件確定，結果需要求解積分方程。對一些簡單情況可以求解，對一般情況則比較困難。高速電子計算機的出現使這種積分方程的數值解法也有了突破。其主要思想是把積分方程離散化，積分方程代表源、匯等奇點在空間連續分佈的總和。例如，若把機翼和機身表面，分割成若干個小單元，每個單元上的奇點強度取平均值。把這些奇點的總和疊加起來，就得出流場總的效應。因此，它用有限項的求和來代替積分，而最後要解的是一組代數方程。由於基本解都是具有奇點的函數，所以這種方法又稱為有限奇點法或鱗片法。(見有限基本解法)

跨聲速流動的流場是既有亞聲速區又有超聲速區的一種混合流場。在不考慮粘性影響和小擾動的情況下，定常二維速度勢方程是混合型的，即Vφ_xx+φ_yy=0，式中V是來流馬赫數※與φ_x的複雜函數。V>0是亞聲速區(橢圓型)，而V<0為超聲速區(雙曲型)。美國的E.M.穆曼和J.D.科爾在1971年首先採用混合差分格式，並運用鬆弛法成功地解出定常小擾動速度勢方程。混合差分格式就是在亞聲速區用中心差分格式，所有鄰近網點上的條件都會影響計算點，而在超聲速區，則用迎風格式，因為上游迎風網點正好是雙曲型波動方程的依賴區。(見跨聲速流數值計算)

在超聲速流動中，主要問題是如何處理激波。用數值方法處理超聲速流場中的激波現有兩種方法。一是激被捕捉法，另一是激泣裝配法。激波捕捉法對激波本身並不需作任何特殊處理，只是在計算公式中，直接或間接地引進“粘性”項，自動算出激波的位置和強度，以“捕捉”激波。其中又有所謂人工粘性和格式粘性兩種方法。人工粘性方法是J.von諾伊曼和R.D.裏希特邁爾於1950年首先提出的，它是以真實粘性流體的物理理論為基礎的一種自動處理激波近似方法。該法是在激波層內，人為地加入粘性項，使激波間斷變成光滑的過渡區。近年來，在超聲速流動中得到廣泛的應用。格式粘性是通過某種差分格式間接地引入粘性項拉克斯格式。拉克斯-文德霍夫格式和麥克馬克格式都具有類似的效果。激波裝配法是把激波仍當作間斷面來處理，激波前後要滿足激波跳躍條件。但是在普通座標中，它的實現很困難。一般採用座標變換，使激波位置(此時是未知的)和一個座標軸重合，然後把激波看作內邊界。這種處理是比較精確的，但也是很麻煩租不方便的。最好的辦法是把激波捕捉法和激波裝配法結合起來。例如在流場外圍的離體激波用激波裝配法，在流場內的激波用激波捕捉法。(見超聲速無粘繞流數值解)

粘性流功數值解法可參見納維-斯托克斯方程數值解、邊界層方程數值解法和湍流數值計算等。 ^[1]