CFD理論過渡到編程的傻瓜入門教程

CFD理論過渡到編程的傻瓜入門教程（注：這是一篇不懂得誰寫的介紹一維無粘可壓縮Euler方程，以及如何編程實現(xiàn)求解該方程的論文。作者從最基本的概念出發(fā)，進(jìn)一步淺出的解說了控制方程，有限體積格式，MSUCL辦法，限制器，Roe格式等有關(guān)知識。這篇論文我覺得有助于大家學(xué)習(xí)CFD編程的有關(guān)知識，因此推薦給大家。文章的背面附有我寫的程序（C語言），用于求解一維激波管問題，大家有愛好能夠看看（程序中加了注釋闡明）胡偶）借寶地寫幾個小短文，介紹CFD的某些實際的入門知識。重要是由于這里支持Latex，寫起來比較方便。CFD，計算流體力學(xué)，是一種挺難的學(xué)科，涉及流體力學(xué)、數(shù)值分析和計算機(jī)算法，尚有計算機(jī)圖形學(xué)的某些知識。特別是有關(guān)偏微分方程數(shù)值分析的東西，不是那么容易入門。大多數(shù)圖書，片中數(shù)學(xué)原理而不重實際動手，由于作者都把讀者當(dāng)做已經(jīng)掌握基礎(chǔ)知識的科班學(xué)生了。因此數(shù)學(xué)基礎(chǔ)不那么好的讀者往往看得很吃力，看了還不懂得怎么實現(xiàn)。本人當(dāng)年雖說是學(xué)航天工程的，但是那時本科教育已經(jīng)退步，基礎(chǔ)的流體力學(xué)課被砍得只剩余一維氣體動力學(xué)了，因此自學(xué)CFD的時候也是頭暈眼花。不懂得怎么實現(xiàn)，也很難找到教學(xué)代碼——那時候網(wǎng)絡(luò)還不發(fā)達(dá)，只在教研室的故紙堆里搜羅到某些完全沒有注釋，編程風(fēng)格也不好的冗長代碼，硬著頭皮分析。后來網(wǎng)上淘到某些代碼研讀，結(jié)合書籍論文才慢慢入門。能夠說中間沒有老師教，后來賭博士為了混學(xué)分上過CFD專門課程，但是那時候我已經(jīng)都掌握課堂上那些了?；叵胱约喝腴T艱辛，不免有一種想法——寫點通俗易懂的CFD入門短文給師弟師妹們。本人不打算搞得很系統(tǒng)，而是但愿能結(jié)合實際，闡明某些最基本的概念和手段，其中某些復(fù)雜的道理只是點到為止?，F(xiàn)在也沒有具體的計劃，想到哪里寫到哪里，因此可能會很零碎。但是我爭取讓初學(xué)CFD的人能夠理解某些基本的東西，看過之后，會懂得一種CFD代碼怎么煉成的（這“煉”字仿佛很流行?。?。歡迎大家提出意見，這樣我盡量的能夠追加某些修改和解釋。言歸正傳，第一部分，我打算介紹一種最基本的算例，一維激波管問題。說白了就是一根兩端封閉的管子，中間有個隔板，隔板左邊和右邊的氣體狀態(tài)（密度、速度、壓力）不同，忽然把隔板抽去，管子內(nèi)面的氣體怎么運(yùn)動。這是個一維問題，被稱作黎曼間斷問題，仿佛是黎曼最初研究雙曲微分方程的時候提出的一種問題，用一維無粘可壓縮Euler方程就能夠描述了。這里這個方程就是描述的氣體密度、動量和能量隨時間的變化（）與它們各自的流量（密度流量，動量流量，能量流量）隨空間變化（）的關(guān)系。在CFD中普通把這個方程寫成矢量形式這里進(jìn)一步能夠?qū)懗缮⒍刃问揭欢ㄒ煜み@種矢量形式以上是控制方程，下面說說求解思路?？蓧嚎s流動計算中，有限體積（FVM）是最廣泛使用的算法，其它算法多多少少都和FVM有些聯(lián)系或者共通的思路。理解的FVM，學(xué)習(xí)其它高級點的算法（例如現(xiàn)在比較熱門的間斷有限元、譜FVM、譜FDM）就好說點了。針對一種微元控制體，把Euler方程在空間積分用微積分知識能夠得到也就是說控制體內(nèi)氣體狀態(tài)平均值的變化是控制體界面上流通量的成果。因此我們要計算的演化，核心問題是計算控制體界面上的。FVM就是以這個積分關(guān)系式出發(fā)，把整個流場劃分為許多小控制體，每個控制體和周邊相鄰的某個控制體共享一種界面，通過計算每個界面上的通量來得到相鄰控制體之間的影響，一旦每個控制體的變化得到，整個流場的變化也就懂得了。因此，再強(qiáng)調(diào)一次，核心問題是計算控制體界面上的。初學(xué)者會說，這個不難，把界面上的插值得到，然后就能夠計算。有道理！咱們畫個圖，有三個小控制體i-1到i+1，中間的“|”表達(dá)界面，控制體i右邊的界面用表達(dá)，左邊的就是。|i-1|i|i+1|好下個問題：每個小控制體長度都是如何插值計算界面上的？最自然的想法就是：取兩邊的平均值唄，但是很不幸，這是不行的。那么換個辦法？直接平均得到？還是很不行，這樣也不行。我靠，這是為什么？這明明是符合微積分里面的知識??？這個道理有點復(fù)雜，說開了去能夠講一本書，能夠說從50年代到70年代，CFD科學(xué)家就在琢磨這個問題。這里，初學(xué)者只需要記住這個結(jié)論：對于流動問題，不能夠這樣簡樸取平均值來插值或者差分。如果你非要想懂得這個終究，我現(xiàn)在也不想給你講清晰，由于我眼下的目的是讓你快速上手，并且該不刨根問底的時候就不要刨根問底，這也是初學(xué)階段一種重要的學(xué)習(xí)辦法。好了，既然目的只是為了求，我在這里，只告訴你一種計算辦法，也是非常重要、非常流行的一種辦法。簡樸的說，就是假設(shè)流動狀態(tài)在界面是不持續(xù)的，先計算出界面兩邊的值，和，再由它們用某種辦法計算出。上述辦法是非常重要的，是由一種蘇聯(lián)人Godunov在50年代首創(chuàng)的，后來被發(fā)展成為通用Godunov辦法，知名的ENO/WENO就是其中的一種。好了，現(xiàn)在的問題是：1怎么擬定和2怎么計算對于第一種問題，Godunov在他的論文中，是假設(shè)每個控制體中是均勻分布的，因此第二個問題，Godunov采用了精確黎曼解來計算。什么是“精確黎曼解”，就是計算這個激波管問題的精確解。既然有精確解，那還費功夫搞這些FVM算法干什么？由于只有這種簡樸一維問題有精確解，稍微復(fù)雜一點就不行了。精確解也比較麻煩，要分析5種狀況，用牛頓法迭代求解（牛頓法是什么？看數(shù)值計算的書去，哦，算了，現(xiàn)在臨時能夠不必看）。這是最初Godunov的辦法，后來在這個思想的基礎(chǔ)上，多種變體都出來了。也但是是在這兩個問題上做文章，怎么擬定，怎么計算。Godunov假設(shè)的是每個小控制體內(nèi)是均勻分布，也就是所謂分段常數(shù)(piecewiseconstant)，因此后來有分段線性(picewiselinear)或者分段二次分布(picewiseparabolic)，到后來ENO/WENO出來，那這個假設(shè)的多項式次數(shù)就繼續(xù)往上走了。都是用多項式近似的，這是數(shù)值計算中的一種強(qiáng)大工具，你能夠在諸多地方看到這種近似。Godunov用的是精確黎曼解，太復(fù)雜太慢，也不必要，因此后來就有多種近似黎曼解，最有名的是Roe求解器、HLL求解器和Osher求解器，都是對精確黎曼解做的簡化。這個多項式的階數(shù)是和計算精度親密有關(guān)的，階數(shù)越高，誤差就越小。但是普通來說，分段線性就能得到不錯的成果了，因此工程中都是用這個，F(xiàn)luent、Fastran以及NASA的CFL3D、OverFlow都是用這個。而黎曼求解器對精度的影響不是那么大，但是對整個算法的物理合用性有影響，也就是說某種近似黎曼求解器可能對某些流動問題不適宜，例如單純的Roe對于鈍頭體的脫體激波會算得亂七八糟，后來加了熵修正才算搞定。上次（）說到了求解可壓縮流動的一種重要算法，通用Godunov辦法。其兩個重要環(huán)節(jié)就是1怎么擬定和2怎么計算這里我們給出第一點一種具體的實現(xiàn)辦法，就是基于原始變量的MUSCL格式（下列簡稱MUSCL）。它是一種很簡樸的格式，并且含有足夠的精度，NASA知名的CFL3D軟件就是使用了這個格式，大家能夠去它的主頁（）上看手冊，里面空間離散那一章清晰的寫著。MUSCL假設(shè)控制體內(nèi)原始變量（就是）的分布是一次或者二次多項式，如果得到了這個多項式，就能夠求出控制體左右兩個界面的一側(cè)的值和。我們以壓力p為例來闡明怎么構(gòu)造這個多項式。這里我只針對二次多項式來解說，你看完之后必定能自己推導(dǎo)出一次多項式的成果（如果你搞不定，那我對你的智商表達(dá)懷疑）。OK，開始假設(shè)，這個假設(shè)不影響最后結(jié)論，由于你總能夠把一種區(qū)間線性的變換到長度為1的區(qū)間。假設(shè)壓力p在控制體i內(nèi)部的分布是一種二次多項式，控制體i的中心處在處，左右兩個界面就是和。這里先強(qiáng)調(diào)一種問題，在FVM中，每個控制體內(nèi)的求解出來的變量事實上是這個控制體內(nèi)的平均值。因此，。我們懂得，和，等距網(wǎng)格狀況下和處的導(dǎo)數(shù)能夠近似表達(dá)為那么（這里錯了，應(yīng)當(dāng)是2ax+b）由上述三個有關(guān)a，b和c的方程，我們能夠得到這樣就能夠得到f(x)的體現(xiàn)式了，由此能夠算出和普通MUSCL格式寫成以下形式對應(yīng)我們的推導(dǎo)成果（二次多項式假設(shè)）。但是這不是最后形式。如果直接用這個公式，就會造成流場在激波（間斷）附近的振蕩。由于直接用二次多項式去逼近一種間斷，會造成這樣的效果。因此科學(xué)家們提出要對間斷附近的斜率有所限制，因此引入了一種非常重要的修改——斜率限制器。加入斜率限制器后，上述公式就有點變化。這里是VanAlbada限制器是一種小數(shù)（），以避免分母為0。密度和速度通過同樣的辦法來搞定。密度、速度和壓力被稱作原始變量，因此上述辦法是基于原始變量的MUSCL。另外尚有基于特性變量的MUSCL，要復(fù)雜一點，但是被認(rèn)為適合更高精度的格式。然而普通計算中，基于原始變量的MUSCL由于含有足夠的精度、簡樸的形式和較低的代價而被廣泛使用。OK，搞定了。下面進(jìn)入第二點，怎么求。有關(guān)這一點，我不打算做具體介紹了，直接使用現(xiàn)有的近似黎曼解就能夠了，都是通過和計算得到。例如Roe由于形式簡樸，而非常流行。在CFL3D軟件主頁（）上看手冊，附錄C的C.1.3。想了一下，還是把Roe求解器稍微說說吧，力求比較完整。但是不要指望我把Roe求解器解釋清晰，由于這個不是很容易三言兩語說清的。Roe求解器的數(shù)學(xué)形式是這樣的顯然這個公式的第一項是一種中心差分形式，先前說過簡樸的中心差分不可行，因素是耗散局限性造成振蕩，說得通俗點就像一種彈簧，如果缺少耗散（阻尼）它就會始終振蕩。“耗散”這個術(shù)語在激波捕獲格式中是最常見的。第二項的作用就是提供足夠的耗散了。這里和已經(jīng)用MUSCL求得了，的定義在第一講中已經(jīng)介紹了。只有是還沒說過的。這個矩陣能夠?qū)懗商匦跃仃嚭吞匦韵蛄烤仃嚨男问蕉?，和的具體體現(xiàn)式在許多書上都有，并且這里的矩陣體現(xiàn)有問題，因此就不寫了。是由、和代入計算得到。而、和采用所謂Roe平均值這才是Roe求解器核心的地方！總結(jié)一下，就是用Roe平均計算界面上的氣體狀態(tài)，然后計算得到，這樣就能夠得到了。如果有時間，我背面會找一種代碼逐句分析一下。總之，計算還是很不直接的。構(gòu)造近似黎曼解是挺有學(xué)問的，需要對氣體動力學(xué)的物理和數(shù)學(xué)方面有較深的理解。普通，如果不是做基礎(chǔ)研究，你只需要懂得它們的特點，會用它們就能夠了，而不必深究它們怎么推導(dǎo)出來的。附錄程序：（新建一種.c類型文獻(xiàn)，將下面的程序復(fù)制粘貼到里面，就能夠運(yùn)行了）/****************************************************本程序用于求解一維無粘可壓縮歐拉方程（激波管問題）運(yùn)用DummyCell解決邊界條件;通量計算方式:AUSMScheme;重構(gòu)辦法：MUSCL辦法限制器：VanAlbada限制器時間離散：四步Runge-Kutta辦法****************************************************/#include"stdio.h"#include"conio.h"#include"malloc.h"#include"stdlib.h"#include"math.h"#include"string.h"#defineh(1/400.0)//網(wǎng)格步長#defineNc404//網(wǎng)格總數(shù)：與h之間的關(guān)系Nc=1/h+4#definePI3.1415927#defineIt1000#definegama1.4//氣體比熱比doubleKAKA=0.0;//限制器控制參數(shù)doubleXS1,XS2;//計算域的兩個端點doubledt=2.5e-5;//時間步長doubletimesum;//總的計算時間//函數(shù)聲明voidoutput();voidSolveWtoU(doubleW[3],doubleU[3]);voidSolveUtoW(doubleW[3],doubleU[3]);//基本變量和守恒變量之間的轉(zhuǎn)換函數(shù)//前后各留兩個網(wǎng)格單元作為虛擬網(wǎng)格單元計算網(wǎng)格單元從2到NOc-3structcell{ intflag;//網(wǎng)格點的類型doublexc;doubleW[3],Wp[3];//conservationvaraibledoubleU[3];//jibenbianliangdoubleR[3];doubleS;//entropy};structcellcell[Nc];//網(wǎng)格生成及流場初始化voidinitialsolve(){inti;doublex,xi,xe;XS1=-0.0;XS2=1.0;xi=XS1-2*h;xe=XS2+2*h; for(i=0;i<Nc;i++){cell[i].xc=(xi+h/2)+i*h;//網(wǎng)格中心坐標(biāo)cell[i].flag=0; x=cell[i].xc;if(x<0.5){cell[i].U[0]=0.445;cell[i].U[1]=0.698;cell[i].U[2]=3.528; }else{cell[i].U[0]=0.5;cell[i].U[1]=0.0;cell[i].U[2]=0.571;}SolveUtoW(cell[i].W,cell[i].U);// printf("x=%f,d=%f,u=%f,p=%f\n",cell[i].xc,cell[i].U[0],cell[i].U[1],cell[i].U[2]);getchar();}}//邊界條件：DummyCell辦法voidboundarycondition(){inti;doubleU[3];//直接賦遠(yuǎn)場值for(i=0;i<Nc;i++){if(cell[i].flag==1){cell[i].U[0]=0.445;cell[i].U[1]=0.698;cell[i].U[2]=3.528; SolveUtoW(cell[i].W,U);}elseif(cell[i].flag==2){cell[i].U[0]=0.5;cell[i].U[1]=0.0;cell[i].U[2]=0.571;SolveUtoW(cell[i].W,U);}}}//重構(gòu)：MUSCL辦法+VanAlbada限制器voidSolveReconstruction(doubleUm1[],doubleUI[],doubleUp1[],doubleUp2[],doubleUL[],doubleUR[]){inti;doubleepsilon,aR[3],aL[3],bL[3],bR[3],sL[3],sR[3];epsilon=1.0e-5*h; for(i=0;i<3;i++){ aR[i]=Up2[i]-Up1[i];//Dealta+U(I+1) bR[i]=Up1[i]-UI[i];//Dealta_U(I+1) aL[i]=Up1[i]-UI[i];//Dealta+U(I) bL[i]=UI[i]-Um1[i];//Dealta_U(I)}for(i=0;i<3;i++){sR[i]=(2.0*aR[i]*bR[i]+1.0e-6)/(aR[i]*aR[i]+bR[i]*bR[i]+1.0e-6);sL[i]=(2.0*aL[i]*bL[i]+1.0e-6)/(aL[i]*aL[i]+bL[i]*bL[i]+1.0e-6);}for(i=0;i<3;i++){UR[i]=Up1[i]-0.25*sR[i]*((1+KAKA*sR[i])*bR[i]+(1-KAKA*sR[i])*aR[i]);UL[i]=UI[i]+0.25*sL[i]*((1+KAKA*sL[i])*aL[i]+(1-KAKA*sL[i])*bL[i]);}}//守恒變量轉(zhuǎn)化為基本變量voidSolveWtoU(doubleW[3],doubleU[3]){ U[0]=W[0];//d U[1]=W[1]/W[0];//u U[2]=(gama-1.0)*(W[2]-0.5*U[0]*U[1]*U[1]);}//基本變量轉(zhuǎn)化為守恒變量voidSolveUtoW(doubleW[3],doubleU[3]){ W[0]=U[0];//d W[1]=U[1]*U[0];//u W[2]=U[2]/(gama-1.0)+0.5*U[0]*U[1]*U[1];}//格式：AUSM格式voidSolveAUSMFlux(doubleUL[],doubleUR[],doubleFc[])//UL左狀態(tài)，UR右狀態(tài)，F(xiàn)C通量{doubleWL[3]={0.0},WR[3]={0.0}; doubledeta,fMn;doubleML,cl,MR,cr,Ml,Mr,Mn;doublepl,pr,pn,PL,PR; deta=0.25;SolveUtoW(WL,UL);SolveUtoW(WR,UR);PL=UL[2];//計算左單元的壓力cl=sqrt(gama*PL/UL[0]);PR=UR[2];//計算右單元的壓力cr=sqrt(gama*PR/UR[0]);ML=UL[1]/cl;MR=UR[1]/cr;if(ML>=1.0) {Ml=ML;pl=PL;} elseif(fabs(ML)<1.0) {Ml=0.25*(ML+1)*(ML+1);pl=0.25*PL*(ML+1)*(ML+1)*(2.0-ML);}else {Ml=0.0;pl=0.0;}if(MR>=1.0) {Mr=0.0;pr=0.0;} elseif(fabs(MR)<1.0) {Mr=-0.25*(MR-1)*(MR-1);pr=0.25*PR*(MR-1)*(MR-1)*(2.0+MR);} else {Mr=MR;pr=PR;}Mn=Ml+Mr; pn=pl+pr;if(fabs(Mn)>deta)fMn=fabs(Mn); elsefMn=0.5*(fabs(Mn)*fabs(Mn)+deta*deta)/deta;Fc[0]=0.5*Mn*(WL[0]*cl+WR[0]*cr)-0.5*fMn*(WR[0]*cr-WL[0]*cl);Fc[1]=0.5*Mn*(WL[1]*cl+WR[1]*cr)-0.5*fMn*(WR[1]*cr-WL[1]*cl)+pn; Fc[2]=0.5*Mn*((WL[2]+PL)*cl+(WR[2]+PR)*cr)-0.5*fMn*((WR[2]+PR)*cr-(WL[2]+PL)*cl);}//殘差計算RvoidSolveResidual(){inti,j;doubleUL[3],UR[3],Fp[3],Fm[3];for(i=0;i<Nc;i++){if(cell[i].flag==0){//i+1/2SolveReconstruction(cell[i-1].U,cell[i].U,cell[i+1].U,cell[i+2].U,UL,UR);//MSUCL重構(gòu)SolveAUSMFlux(UL,UR,Fp);//二階格式// SolveAUSMFlux(cell[i].U,cell[i+1].U,Fp);//一階格式//i-1/2SolveReconstruction(cell[i-2].U,cell[i-1].U,cell[i].U,cell[i+1].U,UL,UR);//MSUCL重構(gòu) SolveAUSMFlux(UL,UR,Fm);//二階格式// SolveAUSMFlux(cell[i-1].U,cell[i].U,Fm);//一階格式for(j=0;j<3;j++){cell[i].R[j]=-(Fp[j]-Fm[j])/h;}}}}//流場值更新voidSolveNextstep(doublear[],intir){inti,j;for(i=0;i<=Nc;i++){if(cell[i].flag==0){ if(ir==0) { for(j=0;j<3;j++) cell[i].Wp[j]=cell[i].W[j]; } for(j=0;j<3;j++) cell[i].W[j]=cell[i].Wp[j]+ar[ir]*dt*cell[i].R[j];SolveWtoU(cell[i].W,cell[i].U);}}}//Runge-Kutta辦法voidSolveRungeKutta(){doublear[4]={1/4.0,1/3.0,0.5,1.0};intit,ir;timesum=0.0;for(it=0;;it++)//迭代步數(shù){ for(ir=0;ir<4;ir++)//四步Rouger-Kutta法{ boundarycondition();SolveResidual(); SolveNextstep(ar,ir); }timesum+=dt;printf("it=%d,timesum=%f\n",it,timesum);if(timesum>=0.16){output();getchar();printf("Pleaseenteranykeytocontinue...");break;}}}//成果輸出voidoutput(){inti;FILE*fpd,*fpu,*fpp;if((fpd=fopen("Density.plt","w"))==NULL) {printf("connotopeninfile\n");return; }fprintf(fpd,"TITLE=\"TestCase\"\n");fprintf(fpd,"VARIABLES=\"x\",\"density\"\n");fprintf(fpd,"ZONET=\"OnlyZone\",I=%d,F=POINT\n",Nc);if((fpp=fopen("Pressure.plt","w"))==NULL) {printf("connotopeninfile\n");return; }fprintf(fpp,"TITLE=\"TestCase\"\n");fp