CFD理論過渡到編程的傻瓜入門教程.doc
《CFD理論過渡到編程的傻瓜入門教程.doc》由會(huì)員分享,可在線閱讀,更多相關(guān)《CFD理論過渡到編程的傻瓜入門教程.doc(19頁珍藏版)》請?jiān)谘b配圖網(wǎng)上搜索。
CFD理論過渡到編程的傻瓜入門教程(注:這是一篇不知道誰寫的介紹一維無粘可壓縮Euler方程,以及如何編程實(shí)現(xiàn)求解該方程的論文。作者從最基本的概念出發(fā),深入淺出的講解了控制方程,有限體積格式,MSUCL方法,限制器,Roe格式等相關(guān)知識(shí)。這篇論文我覺得有利于大家學(xué)習(xí)CFD編程的相關(guān)知識(shí),所以推薦給大家。文章的后面附有我寫的程序(C語言),用于求解一維激波管問題,大家有興趣可以看看(程序中加了注釋說明)胡偶2011)借寶地寫幾個(gè)小短文,介紹CFD的一些實(shí)際的入門知識(shí)。主要是因?yàn)檫@里支持Latex,寫起來比較方便。CFD,計(jì)算流體力學(xué),是一個(gè)挺難的學(xué)科,涉及流體力學(xué)、數(shù)值分析和計(jì)算機(jī)算法,還有計(jì)算機(jī)圖形學(xué)的一些知識(shí)。尤其是有關(guān)偏微分方程數(shù)值分析的東西,不是那么容易入門。大多數(shù)圖書,片中數(shù)學(xué)原理而不重實(shí)際動(dòng)手,因?yàn)樽髡叨及炎x者當(dāng)做已經(jīng)掌握基礎(chǔ)知識(shí)的科班學(xué)生了。所以數(shù)學(xué)基礎(chǔ)不那么好的讀者往往看得很吃力,看了還不知道怎么實(shí)現(xiàn)。本人當(dāng)年雖說是學(xué)航天工程的,但是那時(shí)本科教育已經(jīng)退步,基礎(chǔ)的流體力學(xué)課被砍得只剩下一維氣體動(dòng)力學(xué)了,因此自學(xué)CFD的時(shí)候也是頭暈眼花。不知道怎么實(shí)現(xiàn),也很難找到教學(xué)代碼那時(shí)候網(wǎng)絡(luò)還不發(fā)達(dá),只在教研室的故紙堆里搜羅到一些完全沒有注釋,編程風(fēng)格也不好的冗長代碼,硬著頭皮分析。后來網(wǎng)上淘到一些代碼研讀,結(jié)合書籍論文才慢慢入門??梢哉f中間沒有老師教,后來賭博士為了混學(xué)分上過CFD專門課程,不過那時(shí)候我已經(jīng)都掌握課堂上那些了?;叵胱约喝腴T艱辛,不免有一個(gè)想法寫點(diǎn)通俗易懂的CFD入門短文給師弟師妹們。本人不打算搞得很系統(tǒng),而是希望能結(jié)合實(shí)際,闡明一些最基本的概念和手段,其中一些復(fù)雜的道理只是點(diǎn)到為止。目前也沒有具體的計(jì)劃,想到哪里寫到哪里,因此可能會(huì)很零散。但是我爭取讓初學(xué)CFD的人能夠了解一些基本的東西,看過之后,會(huì)知道一個(gè)CFD代碼怎么煉成的(這“煉”字好像很流行啊)。歡迎大家提出意見,這樣我盡可能的可以追加一些修改和解釋。言歸正傳,第一部分,我打算介紹一個(gè)最基本的算例,一維激波管問題。說白了就是一根兩端封閉的管子,中間有個(gè)隔板,隔板左邊和右邊的氣體狀態(tài)(密度、速度、壓力)不一樣,突然把隔板抽去,管子內(nèi)面的氣體怎么運(yùn)動(dòng)。這是個(gè)一維問題,被稱作黎曼間斷問題,好像是黎曼最初研究雙曲微分方程的時(shí)候提出的一個(gè)問題,用一維無粘可壓縮Euler方程就可以描述了。這里這個(gè)方程就是描述的氣體密度、動(dòng)量和能量隨時(shí)間的變化()與它們各自的流量(密度流量,動(dòng)量流量,能量流量)隨空間變化()的關(guān)系。在CFD中通常把這個(gè)方程寫成矢量形式這里進(jìn)一步可以寫成散度形式一定要熟悉這種矢量形式以上是控制方程,下面說說求解思路??蓧嚎s流動(dòng)計(jì)算中,有限體積(FVM)是最廣泛使用的算法,其他算法多多少少都和FVM有些聯(lián)系或者共通的思路。了解的FVM,學(xué)習(xí)其他高級(jí)點(diǎn)的算法(比如目前比較熱門的間斷有限元、譜FVM、譜FDM)就好說點(diǎn)了。針對(duì)一個(gè)微元控制體,把Euler方程在空間積分用微積分知識(shí)可以得到也就是說控制體內(nèi)氣體狀態(tài)平均值的變化是控制體界面上流通量的結(jié)果。因此我們要計(jì)算的演化,關(guān)鍵問題是計(jì)算控制體界面上的。FVM就是以這個(gè)積分關(guān)系式出發(fā),把整個(gè)流場劃分為許多小控制體,每個(gè)控制體和周圍相鄰的某個(gè)控制體共享一個(gè)界面,通過計(jì)算每個(gè)界面上的通量來得到相鄰控制體之間的影響,一旦每個(gè)控制體的變化得到,整個(gè)流場的變化也就知道了。所以,再強(qiáng)調(diào)一次,關(guān)鍵問題是計(jì)算控制體界面上的。初學(xué)者會(huì)說,這個(gè)不難,把界面上的插值得到,然后就可以計(jì)算。有道理!咱們畫個(gè)圖,有三個(gè)小控制體 i-1到i+1,中間的“|”表示界面,控制體i右邊的界面用表示,左邊的就是。| i-1 | i | i+1 |好下個(gè)問題:每個(gè)小控制體長度都是如何插值計(jì)算界面上的?最自然的想法就是:取兩邊的平均值唄,但是很不幸,這是不行的。那么換個(gè)方法?直接平均得到?還是很不行,這樣也不行。我靠,這是為什么?這明明是符合微積分里面的知識(shí)啊?這個(gè)道理有點(diǎn)復(fù)雜,說開了去可以講一本書,可以說從50年代到70年代,CFD科學(xué)家就在琢磨這個(gè)問題。這里,初學(xué)者只需要記住這個(gè)結(jié)論:對(duì)于流動(dòng)問題,不可以這樣簡單取平均值來插值或者差分。如果你非要想知道這個(gè)究竟,我現(xiàn)在也不想給你講清楚,因?yàn)槲已巯碌哪康氖亲屇憧焖偕鲜?,而且該不刨根問底的時(shí)候就不要刨根問底,這也是初學(xué)階段一種重要的學(xué)習(xí)方法。好了,既然目的只是為了求,我在這里,只告訴你一種計(jì)算方法,也是非常重要、非常流行的一種方法。簡單的說,就是假設(shè)流動(dòng)狀態(tài)在界面是不連續(xù)的,先計(jì)算出界面兩邊的值,和,再由它們用某種方法計(jì)算出。上述方法是非常重要的,是由一個(gè)蘇聯(lián)人Godunov在50年代首創(chuàng)的,后來被發(fā)展成為通用Godunov方法,著名的ENO/WENO就是其中的一種。好了,現(xiàn)在的問題是:1 怎么確定和2 怎么計(jì)算對(duì)于第一個(gè)問題,Godunov在他的論文中,是假設(shè)每個(gè)控制體中是均勻分布的,因此第二個(gè)問題,Godunov采用了精確黎曼解來計(jì)算。什么是“精確黎曼解”,就是計(jì)算這個(gè)激波管問題的精確解。既然有精確解,那還費(fèi)功夫搞這些FVM算法干什么?因?yàn)橹挥羞@種簡單一維問題有精確解,稍微復(fù)雜一點(diǎn)就不行了。精確解也比較麻煩,要分析5種情況,用牛頓法迭代求解(牛頓法是什么?看數(shù)值計(jì)算的書去,哦,算了,現(xiàn)在暫時(shí)可以不必看)。這是最初Godunov的方法,后來在這個(gè)思想的基礎(chǔ)上,各種變體都出來了。也不過是在這兩個(gè)問題上做文章,怎么確定,怎么計(jì)算。Godunov假設(shè)的是每個(gè)小控制體內(nèi)是均勻分布,也就是所謂分段常數(shù)(piecewise constant),所以后來有分段線性(picewise linear)或者分段二次分布(picewise parabolic),到后來ENO/WENO出來,那這個(gè)假設(shè)的多項(xiàng)式次數(shù)就繼續(xù)往上走了。都是用多項(xiàng)式近似的,這是數(shù)值計(jì)算中的一個(gè)強(qiáng)大工具,你可以在很多地方看到這種近似。Godunov用的是精確黎曼解,太復(fù)雜太慢,也不必要,所以后來就有各種近似黎曼解,最有名的是Roe求解器、HLL求解器和Osher求解器,都是對(duì)精確黎曼解做的簡化。這個(gè)多項(xiàng)式的階數(shù)是和計(jì)算精度密切相關(guān)的,階數(shù)越高,誤差就越小。不過一般來說,分段線性就能得到不錯(cuò)的結(jié)果了,所以工程中都是用這個(gè),F(xiàn)luent、Fastran以及NASA的CFL3D、OverFlow都是用這個(gè)。而黎曼求解器對(duì)精度的影響不是那么大,但是對(duì)整個(gè)算法的物理適用性有影響,也就是說某種近似黎曼求解器可能對(duì)某些流動(dòng)問題不合適,比如單純的Roe對(duì)于鈍頭體的脫體激波會(huì)算得亂七八糟,后來加了熵修正才算搞定。上次(http:/gezhi.org/node/399)說到了求解可壓縮流動(dòng)的一個(gè)重要算法,通用Godunov方法。其兩個(gè)主要步驟就是1 怎么確定和2 怎么計(jì)算這里我們給出第一點(diǎn)一個(gè)具體的實(shí)現(xiàn)方法,就是基于原始變量的MUSCL格式(以下簡稱MUSCL)。它是一種很簡單的格式,而且具有足夠的精度,NASA著名的CFL3D軟件就是使用了這個(gè)格式,大家可以去它的主頁(http:/cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html)上看手冊,里面空間離散那一章清楚的寫著。MUSCL假設(shè)控制體內(nèi)原始變量(就是)的分布是一次或者二次多項(xiàng)式,如果得到了這個(gè)多項(xiàng)式,就可以求出控制體左右兩個(gè)界面的一側(cè)的值和。我們以壓力p為例來說明怎么構(gòu)造這個(gè)多項(xiàng)式。這里我只針對(duì)二次多項(xiàng)式來講解,你看完之后肯定能自己推導(dǎo)出一次多項(xiàng)式的結(jié)果(如果你搞不定,那我對(duì)你的智商表示懷疑)。OK,開始假設(shè),這個(gè)假設(shè)不影響最終結(jié)論,因?yàn)槟憧偪梢园岩粋€(gè)區(qū)間線性的變換到長度為1的區(qū)間。假設(shè)壓力p在控制體i內(nèi)部的分布是一個(gè)二次多項(xiàng)式,控制體i的中心處于處,左右兩個(gè)界面就是和。這里先強(qiáng)調(diào)一個(gè)問題,在FVM中,每個(gè)控制體內(nèi)的求解出來的變量實(shí)際上是這個(gè)控制體內(nèi)的平均值。所以,。我們知道,和,等距網(wǎng)格情況下和處的導(dǎo)數(shù)可以近似表示為那么 (這里錯(cuò)了,應(yīng)該是2ax+b)由上述三個(gè)有關(guān)a,b和c的方程,我們可以得到這樣就可以得到f(x)的表達(dá)式了,由此可以算出和通常MUSCL格式寫成如下形式對(duì)應(yīng)我們的推導(dǎo)結(jié)果(二次多項(xiàng)式假設(shè))。但是這不是最終形式。如果直接用這個(gè)公式,就會(huì)導(dǎo)致流場在激波(間斷)附近的振蕩。因?yàn)橹苯佑枚味囗?xiàng)式去逼近一個(gè)間斷,會(huì)導(dǎo)致這樣的效果。所以科學(xué)家們提出要對(duì)間斷附近的斜率有所限制,因此引入了一個(gè)非常重要的修改斜率限制器。加入斜率限制器后,上述公式就有點(diǎn)變化。這里是Van Albada限制器是一個(gè)小數(shù)(),以防止分母為0。密度和速度通過同樣的方法來搞定。密度、速度和壓力被稱作原始變量,所以上述方法是基于原始變量的MUSCL。此外還有基于特征變量的MUSCL,要復(fù)雜一點(diǎn),但是被認(rèn)為適合更高精度的格式。然而一般計(jì)算中,基于原始變量的MUSCL由于具有足夠的精度、簡單的形式和較低的代價(jià)而被廣泛使用。OK,搞定了。下面進(jìn)入第二點(diǎn),怎么求。關(guān)于這一點(diǎn),我不打算做詳細(xì)介紹了,直接使用現(xiàn)有的近似黎曼解就可以了,都是通過和計(jì)算得到。比如Roe因?yàn)樾问胶唵?,而非常流行。在CFL3D軟件主頁(http:/cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html)上看手冊,附錄C的C.1.3。想了一下,還是把Roe求解器稍微說說吧,力求比較完整。但是不要指望我把Roe求解器解釋清楚,因?yàn)檫@個(gè)不是很容易三言兩語說清的。Roe求解器的數(shù)學(xué)形式是這樣的顯然這個(gè)公式的第一項(xiàng)是一個(gè)中心差分形式,先前說過簡單的中心差分不可行,原因是耗散不足導(dǎo)致振蕩,說得通俗點(diǎn)就像一個(gè)彈簧,如果缺乏耗散(阻尼)它就會(huì)一直振蕩?!昂纳ⅰ边@個(gè)術(shù)語在激波捕捉格式中是最常見的。第二項(xiàng)的作用就是提供足夠的耗散了。這里和已經(jīng)用MUSCL求得了,的定義在第一講中已經(jīng)介紹了。只有是還沒說過的。這個(gè)矩陣可以寫成特征矩陣和特征向量矩陣的形式而,和的具體表達(dá)式在許多書上都有,而且這里的矩陣表達(dá)有問題,所以就不寫了。是由、和代入計(jì)算得到。而、和采用所謂Roe平均值這才是Roe求解器關(guān)鍵的地方!總結(jié)一下,就是用Roe平均計(jì)算界面上的氣體狀態(tài),然后計(jì)算得到,這樣就可以得到了。如果有時(shí)間,我后面會(huì)找一個(gè)代碼逐句分析一下。總之,計(jì)算還是很不直接的。構(gòu)造近似黎曼解是挺有學(xué)問的,需要對(duì)氣體動(dòng)力學(xué)的物理和數(shù)學(xué)方面有較深的理解。通常,如果不是做基礎(chǔ)研究,你只需要知道它們的特點(diǎn),會(huì)用它們就可以了,而不必深究它們怎么推導(dǎo)出來的。附錄程序:(新建一個(gè).c類型文件,將下面的程序復(fù)制粘貼到里面,就可以運(yùn)行了)/*本程序用于求解一維無粘可壓縮歐拉方程(激波管問題)運(yùn)用Dummy Cell處理邊界條件;通量計(jì)算方式: AUSM Scheme;重構(gòu)方法:MUSCL方法限制器:Van Albada限制器 時(shí)間離散:四步Runge-Kutta方法 */#include stdio.h#include conio.h#include malloc.h#include stdlib.h#include math.h#include string.h#define h (1/400.0) /網(wǎng)格步長#define Nc 404 /網(wǎng)格總數(shù):與h之間的關(guān)系 Nc=1/h+4#define PI 3.1415927 #define It 1000#define gama 1.4 /氣體比熱比double KAKA=0.0;/限制器控制參數(shù)double XS1,XS2;/計(jì)算域的兩個(gè)端點(diǎn)double dt=2.5e-5;/時(shí)間步長double timesum;/總的計(jì)算時(shí)間/函數(shù)聲明void output();void SolveWtoU(double W3,double U3);void SolveUtoW(double W3,double U3);/基本變量和守恒變量之間的轉(zhuǎn)換函數(shù)/前后各留兩個(gè)網(wǎng)格單元作為虛擬網(wǎng)格單元 計(jì)算網(wǎng)格單元從2到NOc-3struct cellint flag;/網(wǎng)格點(diǎn)的類型double xc;double W3,Wp3;/conservation varaibledouble U3;/jibenbianliangdouble R3;double S;/entropy;struct cell cellNc;/網(wǎng)格生成及流場初始化void initialsolve()int i;double x,xi,xe;XS1=-0.0;XS2=1.0;xi=XS1-2*h;xe=XS2+2*h;for(i=0;iNc;i+)celli.xc=(xi+h/2)+i*h;/網(wǎng)格中心坐標(biāo)celli.flag=0;x=celli.xc; if(x0.5)celli.U0=0.445; celli.U1=0.698;celli.U2=3.528;elsecelli.U0=0.5;celli.U1=0.0;celli.U2=0.571;SolveUtoW(celli.W,celli.U);/printf(x=%f,d=%f,u=%f,p=%fn,celli.xc,celli.U0,celli.U1,celli.U2);getchar();/邊界條件:Dummy Cell方法void boundarycondition()int i;double U3;/直接賦遠(yuǎn)場值for(i=0;iNc;i+)if(celli.flag=1)celli.U0=0.445; celli.U1=0.698;celli.U2=3.528;SolveUtoW(celli.W,U);else if(celli.flag=2)celli.U0=0.5;celli.U1=0.0;celli.U2=0.571;SolveUtoW(celli.W,U);/重構(gòu):MUSCL方法+Van Albada限制器void SolveReconstruction(double Um1,double UI,double Up1,double Up2,double UL,double UR)int i;double epsilon,aR3,aL3,bL3,bR3,sL3,sR3;epsilon=1.0e-5*h;for(i=0;i3;i+)aRi=Up2i-Up1i;/Dealta+U(I+1)bRi=Up1i-UIi;/Dealta_U(I+1)aLi=Up1i-UIi;/Dealta+U(I)bLi=UIi-Um1i;/Dealta_U(I) for(i=0;i3;i+)sRi=(2.0*aRi*bRi+1.0e-6)/(aRi*aRi+bRi*bRi+1.0e-6);sLi=(2.0*aLi*bLi+1.0e-6)/(aLi*aLi+bLi*bLi+1.0e-6);for(i=0;i=1.0) Ml=ML;pl=PL;else if(fabs(ML)=1.0) Mr=0.0;pr=0.0;else if(fabs(MR)deta) fMn=fabs(Mn);else fMn=0.5*(fabs(Mn)*fabs(Mn)+deta*deta)/deta;Fc0=0.5*Mn*(WL0*cl+WR0*cr)-0.5*fMn*(WR0*cr-WL0*cl);Fc1=0.5*Mn*(WL1*cl+WR1*cr)-0.5*fMn*(WR1*cr-WL1*cl)+pn;Fc2=0.5*Mn*(WL2+PL)*cl+(WR2+PR)*cr)-0.5*fMn*(WR2+PR)*cr-(WL2+PL)*cl);/殘差計(jì)算Rvoid SolveResidual()int i,j;double UL3,UR3,Fp3,Fm3;for(i=0;iNc;i+)if(celli.flag=0)/i+1/2SolveReconstruction(celli-1.U, celli.U, celli+1.U, celli+2.U, UL, UR); /MSUCL重構(gòu)SolveAUSMFlux(UL,UR,Fp);/二階格式/SolveAUSMFlux(celli.U,celli+1.U,Fp);/一階格式/i-1/2SolveReconstruction(celli-2.U, celli-1.U, celli.U, celli+1.U, UL, UR); /MSUCL重構(gòu)SolveAUSMFlux(UL,UR,Fm); /二階格式 /SolveAUSMFlux(celli-1.U,celli.U,Fm);/一階格式for(j=0;j3;j+)celli.Rj=-(Fpj-Fmj)/h;/流場值更新void SolveNextstep(double ar,int ir)int i,j;for(i=0;i=Nc;i+)if(celli.flag=0)if(ir=0)for(j=0;j3;j+)celli.Wpj=celli.Wj;for(j=0;j3;j+)celli.Wj=celli.Wpj+arir*dt*celli.Rj;SolveWtoU(celli.W,celli.U);/Runge-Kutta方法void SolveRungeKutta()double ar4=1/4.0,1/3.0,0.5,1.0;int it,ir; timesum=0.0; for(it=0;it+)/迭代步數(shù)for(ir=0;ir=0.16)output();getchar();printf(Please enter any key to continue.);break;/結(jié)果輸出void output()int i;FILE *fpd,*fpu,*fpp;if(fpd=fopen(Density.plt,w)=NULL) printf(connot open infilen);return;fprintf(fpd,TITLE = TestCasen); fprintf(fpd,VARIABLES = x, densityn);fprintf(fpd,ZONE T=Only Zone, I=%d, F=POINTn,Nc);if(fpp=fopen(Pressure.plt,w)=NULL) printf(connot open infilen);return;fprintf(fpp,TITLE = TestCasen); fprintf(fpp,VARIABLES = x, pressuren);fprintf(fpp,ZONE T=Only Zone, I=%d, F=POINTn,Nc);if(fpu=fopen(Velocity.plt,w)=NULL) printf(connot open infilen);return;fprintf(fpu,TITLE = TestCasen); fprintf(fpu,VARIABLES = x, velocityn);fprintf(fpu,ZONE T=Only Zone, I=%d, F=POINTn,Nc);for(i=0;iNc;i+)fprintf(fpd,%ft%fn,celli.xc,celli.U0);fprintf(fpu,%ft%fn,celli.xc,celli.U1);fprintf(fpp,%ft%fn,celli.xc,celli.U2);/printf(x=%f,d=%f,u=%f,p=%fn,celli.xc,celli.U0,celli.U1,celli.U2);getchar();fclose(fpd);fclose(fpp);fclose(fpu);/網(wǎng)格類型判斷void SetMeshFlag()int i;for(i=0;i=XS1&celli.xcXS2)celli.flag=2;elsecelli.flag=1;/printf(cell%d.flag=%dn,i,celli.flag);getchar();/熵的計(jì)算void SovleEntropy()int i;for(i=0;iNc;i+)celli.S=celli.U2/pow(celli.U0,gama);/主函數(shù)void main(void)initialsolve();SetMeshFlag(); SolveRungeKutta();output();- 1.請仔細(xì)閱讀文檔,確保文檔完整性,對(duì)于不預(yù)覽、不比對(duì)內(nèi)容而直接下載帶來的問題本站不予受理。
- 2.下載的文檔,不會(huì)出現(xiàn)我們的網(wǎng)址水印。
- 3、該文檔所得收入(下載+內(nèi)容+預(yù)覽)歸上傳者、原創(chuàng)作者;如果您是本文檔原作者,請點(diǎn)此認(rèn)領(lǐng)!既往收益都?xì)w您。
下載文檔到電腦,查找使用更方便
32 積分
下載 |
- 配套講稿:
如PPT文件的首頁顯示word圖標(biāo),表示該P(yáng)PT已包含配套word講稿。雙擊word圖標(biāo)可打開word文檔。
- 特殊限制:
部分文檔作品中含有的國旗、國徽等圖片,僅作為作品整體效果示例展示,禁止商用。設(shè)計(jì)者僅對(duì)作品中獨(dú)創(chuàng)性部分享有著作權(quán)。
- 關(guān) 鍵 詞:
- CFD 理論 過渡 編程 傻瓜 入門教程
鏈接地址:http://appdesigncorp.com/p-1544534.html