計(jì)算流體力學(xué)過(guò)渡到編程的傻瓜入門教程

借寶地寫幾個(gè)小短文，介紹CFD的某些實(shí)際的入門知識(shí)。重要是由于這里支持Latex，寫起來(lái)比較方便。CFD，計(jì)算流體力學(xué)，是一種挺難的學(xué)科，涉及流體力學(xué)、數(shù)值分析和計(jì)算機(jī)算法，尚有計(jì)算機(jī)圖形學(xué)的某些知識(shí)。特別是有關(guān)偏微分方程數(shù)值分析的東西，不是那么容易入門。大多數(shù)圖書，片中數(shù)學(xué)原理而不重實(shí)際動(dòng)手，由于作者都把讀者當(dāng)做已經(jīng)掌握基礎(chǔ)知識(shí)的科班學(xué)生了。因此數(shù)學(xué)基礎(chǔ)不那么好的讀者往往看得很吃力，看了還不懂得怎么實(shí)現(xiàn)。本人當(dāng)年雖說(shuō)是學(xué)航天工程的，但是那時(shí)本科教育已經(jīng)退步，基礎(chǔ)的流體力學(xué)課被砍得只剩余一維氣體動(dòng)力學(xué)了，因此自學(xué)CFD的時(shí)候也是頭暈眼花。不懂得怎么實(shí)現(xiàn)，也很難找到教學(xué)代碼——那時(shí)候網(wǎng)絡(luò)還不發(fā)達(dá)，只在教研室的故紙堆里搜羅到某些完全沒(méi)有注釋，編程風(fēng)格也不好的冗長(zhǎng)代碼，硬著頭皮分析。后來(lái)網(wǎng)上淘到某些代碼研讀，結(jié)合書籍論文才慢慢入門。能夠說(shuō)中間沒(méi)有老師教，后來(lái)賭博士為了混學(xué)分上過(guò)CFD專門課程，但是那時(shí)候我已經(jīng)都掌握課堂上那些了。回想自己入門艱辛，不免有一種想法——寫點(diǎn)通俗易懂的CFD入門短文給師弟師妹們。本人不打算搞得很系統(tǒng)，而是但愿能結(jié)合實(shí)際，闡明某些最基本的概念和手段，其中某些復(fù)雜的道理只是點(diǎn)到為止?，F(xiàn)在也沒(méi)有具體的計(jì)劃，想到哪里寫到哪里，因此可能會(huì)很零碎。但是我爭(zhēng)取讓初學(xué)CFD的人能夠理解某些基本的東西，看過(guò)之后，會(huì)懂得一種CFD代碼怎么煉成的（這“煉”字仿佛很流行?。g迎大家提出意見(jiàn)，這樣我盡量的能夠追加某些修改和解釋。言歸正傳，第一部分，我打算介紹一種最基本的算例，一維激波管問(wèn)題。說(shuō)白了就是一根兩端封閉的管子，中間有個(gè)隔板，隔板左邊和右邊的氣體狀態(tài)（密度、速度、壓力）不同，忽然把隔板抽去，管子內(nèi)面的氣體怎么運(yùn)動(dòng)。這是個(gè)一維問(wèn)題，被稱作黎曼間斷問(wèn)題，仿佛是黎曼最初研究雙曲微分方程的時(shí)候提出的一種問(wèn)題，用一維無(wú)粘可壓縮Euler方程就能夠描述了。這里這個(gè)方程就是描述的氣體密度、動(dòng)量和能量隨時(shí)間的變化（）與它們各自的流量（密度流量，動(dòng)量流量，能量流量）隨空間變化（）的關(guān)系。在CFD中普通把這個(gè)方程寫成矢量形式這里進(jìn)一步能夠?qū)懗缮⒍刃问揭欢ㄒ煜み@種矢量形式以上是控制方程，下面說(shuō)說(shuō)求解思路?？蓧嚎s流動(dòng)計(jì)算中，有限體積（FVM）是最廣泛使用的算法，其它算法多多少少都和FVM有些聯(lián)系或者共通的思路。理解的FVM，學(xué)習(xí)其它高級(jí)點(diǎn)的算法（例如現(xiàn)在比較熱門的間斷有限元、譜FVM、譜FDM）就好說(shuō)點(diǎn)了。針對(duì)一種微元控制體，把Euler方程在空間積分用微積分知識(shí)能夠得到也就是說(shuō)控制體內(nèi)氣體狀態(tài)平均值的變化是控制體界面上流通量的成果。因此我們要計(jì)算的演化，核心問(wèn)題是計(jì)算控制體界面上的。FVM就是以這個(gè)積分關(guān)系式出發(fā)，把整個(gè)流場(chǎng)劃分為許多小控制體，每個(gè)控制體和周邊相鄰的某個(gè)控制體共享一種界面，通過(guò)計(jì)算每個(gè)界面上的通量來(lái)得到相鄰控制體之間的影響，一旦每個(gè)控制體的變化得到，整個(gè)流場(chǎng)的變化也就懂得了。因此，再?gòu)?qiáng)調(diào)一次，核心問(wèn)題是計(jì)算控制體界面上的。初學(xué)者會(huì)說(shuō)，這個(gè)不難，把界面上的插值得到，然后就能夠計(jì)算。有道理！咱們畫個(gè)圖，有三個(gè)小控制體i-1到i+1，中間的“|”表達(dá)界面，控制體i右邊的界面用表達(dá)，左邊的就是。|i-1|i|i+1|好下個(gè)問(wèn)題：每個(gè)小控制體長(zhǎng)度都是如何插值計(jì)算界面上的？最自然的想法就是：取兩邊的平均值唄，但是很不幸，這是不行的。那么換個(gè)辦法？直接平均得到？還是很不行，這樣也不行。我靠，這是為什么？這明明是符合微積分里面的知識(shí)??？這個(gè)道理有點(diǎn)復(fù)雜，說(shuō)開(kāi)了去能夠講一本書，能夠說(shuō)從50年代到70年代，CFD科學(xué)家就在琢磨這個(gè)問(wèn)題。這里，初學(xué)者只需要記住這個(gè)結(jié)論：對(duì)于流動(dòng)問(wèn)題，不能夠這樣簡(jiǎn)樸取平均值來(lái)插值或者差分。如果你非要想懂得這個(gè)終究，我現(xiàn)在也不想給你講清晰，由于我眼下的目的是讓你快速上手，并且該不刨根問(wèn)底的時(shí)候就不要刨根問(wèn)底，這也是初學(xué)階段一種重要的學(xué)習(xí)辦法。好了，既然目的只是為了求，我在這里，只告訴你一種計(jì)算辦法，也是非常重要、非常流行的一種辦法。簡(jiǎn)樸的說(shuō)，就是假設(shè)流動(dòng)狀態(tài)在界面是不持續(xù)的，先計(jì)算出界面兩邊的值，和，再由它們用某種辦法計(jì)算出。上述辦法是非常重要的，是由一種蘇聯(lián)人Godunov在50年代首創(chuàng)的，后來(lái)被發(fā)展成為通用Godunov辦法，出名的ENO/WENO就是其中的一種。好了，現(xiàn)在的問(wèn)題是：1怎么擬定和2怎么計(jì)算對(duì)于第一種問(wèn)題，Godunov在他的論文中，是假設(shè)每個(gè)控制體中是均勻分布的，因此第二個(gè)問(wèn)題，Godunov采用了精確黎曼解來(lái)計(jì)算。什么是“精確黎曼解”，就是計(jì)算這個(gè)激波管問(wèn)題的精確解。既然有精確解，那還費(fèi)功夫搞這些FVM算法干什么？由于只有這種簡(jiǎn)樸一維問(wèn)題有精確解，稍微復(fù)雜一點(diǎn)就不行了。精確解也比較麻煩，要分析5種狀況，用牛頓法迭代求解（牛頓法是什么？看數(shù)值計(jì)算的書去，哦，算了，現(xiàn)在臨時(shí)能夠不必看）。這是最初Godunov的辦法，后來(lái)在這個(gè)思想的基礎(chǔ)上，多個(gè)變體都出來(lái)了。也但是是在這兩個(gè)問(wèn)題上做文章，怎么擬定，怎么計(jì)算。Godunov假設(shè)的是每個(gè)小控制體內(nèi)是均勻分布，也就是所謂分段常數(shù)(piecewiseconstant)，因此后來(lái)有分段線性(picewiselinear)或者分段二次分布(picewiseparabolic)，到后來(lái)ENO/WENO出來(lái)，那這個(gè)假設(shè)的多項(xiàng)式次數(shù)就繼續(xù)往上走了。都是用多項(xiàng)式近似的，這是數(shù)值計(jì)算中的一種強(qiáng)大工具，你能夠在諸多地方看到這種近似。Godunov用的四精確黎曼解，太復(fù)雜太慢，也不必要，因此后來(lái)就有多個(gè)近似黎曼解，最有名的是Roe求解器、HLL求解器和Osher求解器，都是對(duì)精確黎曼解做的簡(jiǎn)化。這個(gè)多項(xiàng)式的階數(shù)是和計(jì)算精度親密有關(guān)的，階數(shù)越高，誤差就越小。但是普通來(lái)說(shuō)，分段線性就能得到不錯(cuò)的成果了，因此工程中都是用這個(gè)，F(xiàn)luent、Fastran以及NASA的CFL3D、OverFlow都是用這個(gè)。而黎曼求解器對(duì)精度的影響不是那么大，但是對(duì)整個(gè)算法的物理合用性有影響，也就是說(shuō)某種近似黎曼求解器可能對(duì)某些流動(dòng)問(wèn)題不適宜，例如單純的Roe對(duì)于鈍頭體的脫體激波會(huì)算得亂七八糟，后來(lái)加了熵修正才算搞定。上次（）說(shuō)到了求解可壓縮流動(dòng)的一種重要算法，通用Godunov辦法。其兩個(gè)重要環(huán)節(jié)就是1怎么擬定和2怎么計(jì)算這里我們給出第一點(diǎn)一種具體的實(shí)現(xiàn)辦法，就是基于原始變量的MUSCL格式（下列簡(jiǎn)稱MUSCL）。它是一種很簡(jiǎn)樸的格式，并且含有足夠的精度，NASA出名的CFL3D軟件就是使用了這個(gè)格式，大家能夠去它的主頁(yè)（）上看手冊(cè)，里面空間離散那一章清晰的寫著。MUSCL假設(shè)控制體內(nèi)原始變量（就是）的分布是一次或者二次多項(xiàng)式，如果得到了這個(gè)多項(xiàng)式，就能夠求出控制體左右兩個(gè)界面的一側(cè)的值和。我們以壓力p為例來(lái)闡明怎么構(gòu)造這個(gè)多項(xiàng)式。這里我只針對(duì)二次多項(xiàng)式來(lái)解說(shuō)，你看完之后必定能自己推導(dǎo)出一次多項(xiàng)式的成果（如果你搞不定，那我對(duì)你的智商表達(dá)懷疑）。OK，開(kāi)始假設(shè)，這個(gè)假設(shè)不影響最后結(jié)論，由于你總能夠把一種區(qū)間線性的變換到長(zhǎng)度為1的區(qū)間。假設(shè)壓力p在控制體i內(nèi)部的分布是一種二次多項(xiàng)式，控制體i的中心處在處，左右兩個(gè)界面就是和。這里先強(qiáng)調(diào)一種問(wèn)題，在FVM中，每個(gè)控制體內(nèi)的求解出來(lái)的變量事實(shí)上是這個(gè)控制體內(nèi)的平均值。因此，。我們懂得，和，等距網(wǎng)格狀況下和處的導(dǎo)數(shù)能夠近似表達(dá)為那么由上述三個(gè)有關(guān)a，b和c的方程，我們能夠得到這樣就能夠得到f(x)的體現(xiàn)式了，由此能夠算出和普通MUSCL格式寫成以下形式對(duì)應(yīng)我們的推導(dǎo)成果（二次多項(xiàng)式假設(shè)）。但是這不是最后形式。如果直接用這個(gè)公式，就會(huì)造成流場(chǎng)在激波（間斷）附近的振蕩。由于直接用二次多項(xiàng)式去逼近一種間斷，會(huì)造成這樣的效果。因此科學(xué)家們提出要對(duì)間斷附近的斜率有所限制，因此引入了一種非常重要的修改——斜率限制器。加入斜率限制器后，上述公式就有點(diǎn)變化。這里是VanAlbada限制器是一種小數(shù)（），以避免分母為0。密度和速度通過(guò)同樣的辦法來(lái)搞定。密度、速度和壓力被稱作原始變量，因此上述辦法是基于原始變量的MUSCL。另外尚有基于特性變量的MUSCL，要復(fù)雜一點(diǎn)，但是被認(rèn)為適合更高精度的格式。然而普通計(jì)算中，基于原始變量的MUSCL由于含有足夠的精度、簡(jiǎn)樸的形式和較低的代價(jià)而被廣泛使用。OK，搞定了。下面進(jìn)入第二點(diǎn)，怎么求。有關(guān)這一點(diǎn)，我不打算做具體介紹了，直接使用現(xiàn)有的近似黎曼解就能夠了，都是通過(guò)和計(jì)算得到。例如Roe由于形式簡(jiǎn)樸，而非常流行。在CFL3D軟件主頁(yè)（）上看手冊(cè)，附錄C的C.1.3。想了一下，還是把Roe求解器稍微說(shuō)說(shuō)吧，力求比較完整。但是不要指望我把Roe求解器解釋清晰，由于這個(gè)不是很容易三言兩語(yǔ)說(shuō)清的。Roe求解器的數(shù)學(xué)形式是這樣的顯然這個(gè)公式的第一項(xiàng)是一種中心差分形式，先前說(shuō)過(guò)簡(jiǎn)樸的中心差分不可行，因素是耗散局限性造成振蕩，說(shuō)得通俗點(diǎn)就像一種彈簧，如果缺少耗散（阻尼）它就會(huì)始終振蕩?！昂纳ⅰ边@個(gè)術(shù)語(yǔ)在激波捕獲格式中是最常見(jiàn)的。第二項(xiàng)的作用就是提供足夠的耗散了。這里和已經(jīng)用MUSCL求得了，的定義在第一講中已經(jīng)介紹了。只有是還沒(méi)說(shuō)過(guò)的。這個(gè)矩陣能夠?qū)懗商匦跃仃嚭吞匦韵蛄烤仃嚨男问蕉?，和的具體體現(xiàn)式在許多書上都有，并且這里的矩陣體現(xiàn)有問(wèn)題，因此就不寫了。