CFD理論過渡到編程的最簡單教程(附源程序)

CFD理論過渡到編程的最簡單教程(附源程序)CFD理論過渡到編程的最簡單教程(附源程序)CFD理論過渡到編程的最簡單教程(附源程序)CFD理論過渡到編程的最簡單教程(附源程序)編制僅供參考審核批準(zhǔn)生效日期地址：電話：傳真:郵編:CFD理論過渡到編程的傻瓜入門教程（注：這是一篇不知道誰寫的介紹一維無粘可壓縮Euler方程，以及如何編程實(shí)現(xiàn)求解該方程的論文。作者從最基本的概念出發(fā)，深入淺出的講解了控制方程，有限體積格式，MSUCL方法，限制器，Roe格式等相關(guān)知識。這篇論文我覺得有利于大家學(xué)習(xí)CFD編程的相關(guān)知識，所以推薦給大家。文章的后面附有我寫的程序（C語言），用于求解一維激波管問題，大家有興趣可以看看（程序中加了注釋說明）胡偶2011）借寶地寫幾個小短文，介紹CFD的一些實(shí)際的入門知識。主要是因?yàn)檫@里支持Latex，寫起來比較方便。CFD，計(jì)算流體力學(xué)，是一個挺難的學(xué)科，涉及流體力學(xué)、數(shù)值分析和計(jì)算機(jī)算法，還有計(jì)算機(jī)圖形學(xué)的一些知識。尤其是有關(guān)偏微分方程數(shù)值分析的東西，不是那么容易入門。大多數(shù)圖書，偏重?cái)?shù)學(xué)原理而不重實(shí)際動手，因?yàn)樽髡叨及炎x者當(dāng)做已經(jīng)掌握基礎(chǔ)知識的科班學(xué)生了。所以數(shù)學(xué)基礎(chǔ)不那么好的讀者往往看得很吃力，看了還不知道怎么實(shí)現(xiàn)。本人當(dāng)年雖說是學(xué)航天工程的，但是那時本科教育已經(jīng)退步，基礎(chǔ)的流體力學(xué)課被砍得只剩下一維氣體動力學(xué)了，因此自學(xué)CFD的時候也是頭暈眼花。不知道怎么實(shí)現(xiàn)，也很難找到教學(xué)代碼——那時候網(wǎng)絡(luò)還不發(fā)達(dá)，只在教研室的故紙堆里搜羅到一些完全沒有注釋，編程風(fēng)格也不好的冗長代碼，硬著頭皮分析。后來網(wǎng)上淘到一些代碼研讀，結(jié)合書籍論文才慢慢入門?？梢哉f中間沒有老師教，后來讀博士為了混學(xué)分上過CFD專門課程，不過那時候我已經(jīng)都掌握課堂上那些了?；叵胱约喝腴T艱辛，不免有一個想法——寫點(diǎn)通俗易懂的CFD入門短文給師弟師妹們。本人不打算搞得很系統(tǒng)，而是希望能結(jié)合實(shí)際，闡明一些最基本的概念和手段，其中一些復(fù)雜的道理只是點(diǎn)到為止。目前也沒有具體的計(jì)劃，想到哪里寫到哪里，因此可能會很零散。但是我爭取讓初學(xué)CFD的人能夠了解一些基本的東西，看過之后，會知道一個CFD代碼怎么煉成的（這“煉”字好像很流行?。?。歡迎大家提出意見，這樣我盡可能的可以追加一些修改和解釋。言歸正傳，第一部分，我打算介紹一個最基本的算例，一維激波管問題。說白了就是一根兩端封閉的管子，中間有個隔板，隔板左邊和右邊的氣體狀態(tài)（密度、速度、壓力）不一樣，突然把隔板抽去，管子內(nèi)的氣體怎么運(yùn)動。這是個一維問題，被稱作黎曼間斷問題，好像是黎曼最初研究雙曲微分方程的時候提出的一個問題，用一維無粘可壓縮Euler方程就可以描述了。這里這個方程就是描述的氣體密度、動量和能量隨時間的變化（）與它們各自的流量（密度流量，動量流量，能量流量）隨空間變化（）的關(guān)系。在CFD中通常把這個方程寫成矢量形式

這里

進(jìn)一步可以寫成散度形式

一定要熟悉這種矢量形式以上是控制方程，下面說說求解思路。可壓縮流動計(jì)算中，有限體積（FVM）是最廣泛使用的算法，其他算法多多少少都和FVM有些聯(lián)系或者共同的思路。了解了FVM，學(xué)習(xí)其他高級點(diǎn)的算法（比如目前比較熱門的間斷有限元、譜FVM、譜FDM）就好說點(diǎn)了。針對一個微元控制體，把Euler方程在空間積分

用微積分知識可以得到

也就是說控制體內(nèi)氣體狀態(tài)平均值的變化是控制體界面上流通量的結(jié)果。因此我們要計(jì)算的演化，關(guān)鍵問題是計(jì)算控制體界面上的。FVM就是以這個積分關(guān)系式出發(fā)，把整個流場劃分為許多小控制體，每個控制體和周圍相鄰的某個控制體共享一個界面，通過計(jì)算每個界面上的通量來得到相鄰控制體之間的影響，一旦每個控制體的變化得到，整個流場的變化也就知道了。所以，再強(qiáng)調(diào)一次，關(guān)鍵問題是計(jì)算控制體界面上的。初學(xué)者會說，這個不難，把界面上的插值得到，然后就可以計(jì)算。有道理！咱們畫個圖，有三個小控制體i-1到i+1，中間的“|”表示界面，控制體i右邊的界面用表示，左邊的就是。

|i-1|i|i+1|好，下個問題：每個小控制體長度都是。如何插值計(jì)算界面上的最自然的想法就是：取兩邊的平均值唄，

但是很不幸，這是不行的。那么換個方法直接平均得到？

還是很不行，這樣也不行。我靠，這是為什么這明明是符合微積分里面的知識啊這個道理有點(diǎn)復(fù)雜，說開了去可以講一本書，可以說從50年代到70年代，CFD科學(xué)家就在琢磨這個問題。這里，初學(xué)者只需要記住這個結(jié)論：對于流動問題，不可以這樣簡單取平均值來插值或者差分。如果你非要想知道這個究竟，我現(xiàn)在也不想給你講清楚，因?yàn)槲已巯碌哪康氖亲屇憧焖偕鲜郑以摬慌俑鶈柕椎臅r候就不要刨根問底，這也是初學(xué)階段一種重要的學(xué)習(xí)方法。好了，既然目的只是為了求，我在這里，只告訴你一種計(jì)算方法，也是非常重要、非常流行的一種方法。簡單的說，就是假設(shè)流動狀態(tài)在界面是不連續(xù)的，先計(jì)算出界面兩邊的值，和，再由它們用某種方法計(jì)算出。上述方法是非常重要的，是由一個蘇聯(lián)人Godunov在50年代首創(chuàng)的，后來被發(fā)展成為通用Godunov方法，著名的ENO/WENO就是其中的一種。好了，現(xiàn)在的問題是：

1怎么確定和

2怎么計(jì)算對于第一個問題，Godunov在他的論文中，是假設(shè)每個控制體中是均勻分布的，因此

第二個問題，Godunov采用了精確黎曼解來計(jì)算。什么是“精確黎曼解”，就是計(jì)算這個激波管問題的精確解。既然有精確解，那還費(fèi)功夫搞這些FVM算法干什么因?yàn)橹挥羞@種簡單一維問題有精確解，稍微復(fù)雜一點(diǎn)就不行了。精確解也比較麻煩，要分析5種情況，用牛頓法迭代求解（牛頓法是什么看數(shù)值計(jì)算的書去，哦，算了，現(xiàn)在暫時可以不必看）。這是最初Godunov的方法，后來在這個思想的基礎(chǔ)上，各種變體都出來了。也不過是在這兩個問題上做文章，怎么確定，怎么計(jì)算。Godunov假設(shè)的是每個小控制體內(nèi)是均勻分布，也就是所謂分段常數(shù)(piecewiseconstant)，所以后來有分段線性(picewiselinear)或者分段二次分布(picewiseparabolic)，到后來ENO/WENO出來，那這個假設(shè)的多項(xiàng)式次數(shù)就繼續(xù)往上走了。都是用多項(xiàng)式近似的，這是數(shù)值計(jì)算中的一個強(qiáng)大工具，你可以在很多地方看到這種近似。Godunov用的是精確黎曼解，太復(fù)雜太慢，也不必要，所以后來就有各種近似黎曼解，最有名的是Roe求解器、HLL求解器和Osher求解器，都是對精確黎曼解做的簡化。這個多項(xiàng)式的階數(shù)是和計(jì)算精度密切相關(guān)的，階數(shù)越高，誤差就越小。不過一般來說，分段線性就能得到不錯的結(jié)果了，所以工程中都是用這個，F(xiàn)luent、Fastran以及NASA的CFL3D、OverFlow都是用這個。而黎曼求解器對精度的影響不是那么大，但是對整個算法的物理適用性有影響，也就是說某種近似黎曼求解器可能對某些流動問題不合適，比如單純的Roe對于鈍頭體的脫體激波會算得亂七八糟，后來加了熵修正才算搞定。上次（node/399）說到了求解可壓縮流動的一個重要算法，通用Godunov方法。其兩個主要步驟就是

1怎么確定和

2怎么計(jì)算這里我們給出第一點(diǎn)一個具體的實(shí)現(xiàn)方法，就是基于原始變量的MUSCL格式（以下簡稱MUSCL）。它是一種很簡單的格式，而且具有足夠的精度，NASA著名的CFL3D軟件就是使用了這個格式，大家可以去它的主頁（）上看手冊，里面空間離散那一章清楚的寫著。MUSCL假設(shè)控制體內(nèi)原始變量（就是）的分布是一次或者二次多項(xiàng)式，如果得到了這個多項(xiàng)式，就可以求出控制體左右兩個界面的一側(cè)的值和。我們以壓力p為例來說明怎么構(gòu)造這個多項(xiàng)式。這里我只針對二次多項(xiàng)式來講解，你看完之后肯定能自己推導(dǎo)出一次多項(xiàng)式的結(jié)果（如果你搞不定，那我對你的智商表示懷疑）。OK，開始

假設(shè)，這個假設(shè)不影響最終結(jié)論，因?yàn)槟憧偪梢园岩粋€區(qū)間線性的變換到長度為1的區(qū)間。

假設(shè)壓力p在控制體i內(nèi)部的分布是一個二次多項(xiàng)式，控制體i的中心處于處，左右兩個界面就是和。這里先強(qiáng)調(diào)一個問題，在FVM中，每個控制體內(nèi)的求解出來的變量實(shí)際上是這個控制體內(nèi)的平均值。

所以，

。我們知道，和，等距網(wǎng)格情況下和處的導(dǎo)數(shù)可以近似表示為

那么

（這里錯了，應(yīng)該是2ax+b）

由上述三個有關(guān)a，b和c的方程，我們可以得到

這樣就可以得到f(x)的表達(dá)式了，由此可以算出和通常MUSCL格式寫成如下形式

對應(yīng)我們的推導(dǎo)結(jié)果（二次多項(xiàng)式假設(shè)）。但是這不是最終形式。如果直接用這個公式，就會導(dǎo)致流場在激波（間斷）附近的振蕩。因?yàn)橹苯佑枚味囗?xiàng)式去逼近一個間斷，會導(dǎo)致這樣的效果。所以科學(xué)家們提出要對間斷附近的斜率有所限制，因此引入了一個非常重要的修改——斜率限制器。加入斜率限制器后，上述公式就有點(diǎn)變化。

這里是VanAlbada限制器

是一個小數(shù)（），以防止分母為0。密度和速度通過同樣的方法來搞定。密度、速度和壓力被稱作原始變量，所以上述方法是基于原始變量的MUSCL。此外還有基于特征變量的MUSCL，要復(fù)雜一點(diǎn)，但是被認(rèn)為適合更高精度的格式。然而一般計(jì)算中，基于原始變量的MUSCL由于具有足夠的精度、簡單的形式和較低的代價而被廣泛使用。OK，搞定了。下面進(jìn)入第二點(diǎn)，怎么求。關(guān)于這一點(diǎn)，我不打算做詳細(xì)介紹了，直接使用現(xiàn)有的近似黎曼解就可以了，都是通過和計(jì)算得到。比如Roe因?yàn)樾问胶唵?，而非常流行。在CFL3D軟件主頁（）上看手冊，附錄C的。想了一下，還是把Roe求解器稍微說說吧，力求比較完整。但是不要指望我把Roe求解器解釋清楚，因?yàn)檫@個不是很容易三言兩語說清的。Roe求解器的數(shù)學(xué)形式是這樣的

顯然這個公式的第一項(xiàng)是一個中心差分形式，先前說過簡單的中心差分不可行，原因是耗散不足導(dǎo)致振蕩，說得通俗點(diǎn)就像一個彈簧，如果缺乏耗散（阻尼）它就會一直振蕩?！昂纳ⅰ边@個術(shù)語在激波捕捉格式中是最常見的。第二項(xiàng)的作用就是提供足夠的耗散了。這里和已經(jīng)用MUSCL求得了，的定義在第一講中已經(jīng)介紹了。只有是還沒說過的。這個矩陣可以寫成特征矩陣和特征向量矩陣的形式

而

，，和的具體表達(dá)式在許多書上都有，而且這里的矩陣表達(dá)有問題，所以就不寫了。是由、和代入計(jì)算得到。而、和采用所謂Roe平均值

這才是Roe求解器關(guān)鍵的地方！總結(jié)一下，就是用Roe平均計(jì)算界面上的氣體狀態(tài)，然后計(jì)算得到，這樣就可以得到了。如果有時間，我后面會找一個代碼逐句分析一下?？傊?，計(jì)算還是很不直接的。構(gòu)造近似黎曼解是挺有學(xué)問的，需要對氣體動力學(xué)的物理和數(shù)學(xué)方面有較深的理解。通常，如果不是做基礎(chǔ)研究，你只需要知道它們的特點(diǎn)，會用它們就可以了，而不必深究它們怎么推導(dǎo)出來的。附錄程序：（新建一個.c類型文件，將下面的程序復(fù)制粘貼到里面，就可以運(yùn)行了）/****************************************************本程序用于求解一維無粘可壓縮歐拉方程（激波管問題）運(yùn)用DummyCell處理邊界條件;通量計(jì)算方式:AUSMScheme;重構(gòu)方法：MUSCL方法限制器：VanAlbada限制器時間離散：四步Runge-Kutta方法****************************************************/#include""#include""#include""#include""#include""#include""#defineh(1/c=(xi+h/2)+i*h;lag=0; x=cell[i].xc; if(x< { cell[i].U[0]=;cell[i].U[1]=;cell[i].U[2]=; } else { cell[i].U[0]=;cell[i].U[1]=;cell[i].U[2]=; } SolveUtoW(cell[i].W,cell[i].U);c,cell[i].U[0],cell[i].U[1],cell[i].U[2]);getchar(); }}lag==1) { cell[i].U[0]=;cell[i].U[1]=;cell[i].U[2]=; SolveUtoW(cell[i].W,U); } elseif(cell[i].flag==2) { cell[i].U[0]=;cell[i].U[1]=;cell[i].U[2]=; SolveUtoW(cell[i].W,U); } }}URUR { ,cell[i].U,cell[i+1].U,cell[i+2].U,UL,UR);,cell[i+1].U,Fp);,cell[i-1].U,cell[i].U,cell[i+1].U,UL,UR);,cell[i].U,Fm);[j]=-(Fp[j]-Fm[j])/h; } } }}lag==0) { if(ir==0) { for(j=0;j<3;j++) cell[i].Wp[j]=cell[i].W[j]; } for