CFD理論過渡到編程的傻瓜入門教程.doc
《CFD理論過渡到編程的傻瓜入門教程.doc》由會員分享,可在線閱讀,更多相關《CFD理論過渡到編程的傻瓜入門教程.doc(19頁珍藏版)》請在裝配圖網上搜索。
CFD理論過渡到編程的傻瓜入門教程 (注:這是一篇不知道誰寫的介紹一維無粘可壓縮Euler方程,以及如何編程實現求解該方程的論文。作者從最基本的概念出發(fā),深入淺出的講解了控制方程,有限體積格式,MSUCL方法,限制器,Roe格式等相關知識。這篇論文我覺得有利于大家學習CFD編程的相關知識,所以推薦給大家。文章的后面附有我寫的程序(C語言),用于求解一維激波管問題,大家有興趣可以看看(程序中加了注釋說明)胡偶2011) 借寶地寫幾個小短文,介紹CFD的一些實際的入門知識。主要是因為這里支持Latex,寫起來比較方便。 CFD,計算流體力學,是一個挺難的學科,涉及流體力學、數值分析和計算機算法,還有計算機圖形學的一些知識。尤其是有關偏微分方程數值分析的東西,不是那么容易入門。大多數圖書,片中數學原理而不重實際動手,因為作者都把讀者當做已經掌握基礎知識的科班學生了。所以數學基礎不那么好的讀者往往看得很吃力,看了還不知道怎么實現。本人當年雖說是學航天工程的,但是那時本科教育已經退步,基礎的流體力學課被砍得只剩下一維氣體動力學了,因此自學CFD的時候也是頭暈眼花。不知道怎么實現,也很難找到教學代碼——那時候網絡還不發(fā)達,只在教研室的故紙堆里搜羅到一些完全沒有注釋,編程風格也不好的冗長代碼,硬著頭皮分析。后來網上淘到一些代碼研讀,結合書籍論文才慢慢入門。可以說中間沒有老師教,后來賭博士為了混學分上過CFD專門課程,不過那時候我已經都掌握課堂上那些了。 回想自己入門艱辛,不免有一個想法——寫點通俗易懂的CFD入門短文給師弟師妹們。本人不打算搞得很系統(tǒng),而是希望能結合實際,闡明一些最基本的概念和手段,其中一些復雜的道理只是點到為止。目前也沒有具體的計劃,想到哪里寫到哪里,因此可能會很零散。但是我爭取讓初學CFD的人能夠了解一些基本的東西,看過之后,會知道一個CFD代碼怎么煉成的(這“煉”字好像很流行?。g迎大家提出意見,這樣我盡可能的可以追加一些修改和解釋。 言歸正傳,第一部分,我打算介紹一個最基本的算例,一維激波管問題。說白了就是一根兩端封閉的管子,中間有個隔板,隔板左邊和右邊的氣體狀態(tài)(密度、速度、壓力)不一樣,突然把隔板抽去,管子內面的氣體怎么運動。這是個一維問題,被稱作黎曼間斷問題,好像是黎曼最初研究雙曲微分方程的時候提出的一個問題,用一維無粘可壓縮Euler方程就可以描述了。 這里 這個方程就是描述的氣體密度、動量和能量隨時間的變化()與它們各自的流量(密度流量,動量流量,能量流量)隨空間變化()的關系。 在CFD中通常把這個方程寫成矢量形式 這里 進一步可以寫成散度形式 一定要熟悉這種矢量形式 以上是控制方程,下面說說求解思路??蓧嚎s流動計算中,有限體積(FVM)是最廣泛使用的算法,其他算法多多少少都和FVM有些聯系或者共通的思路。了解的FVM,學習其他高級點的算法(比如目前比較熱門的間斷有限元、譜FVM、譜FDM)就好說點了。 針對一個微元控制體,把Euler方程在空間積分 用微積分知識可以得到 也就是說控制體內氣體狀態(tài)平均值的變化是控制體界面上流通量的結果。因此我們要計算的演化,關鍵問題是計算控制體界面上的。 FVM就是以這個積分關系式出發(fā),把整個流場劃分為許多小控制體,每個控制體和周圍相鄰的某個控制體共享一個界面,通過計算每個界面上的通量來得到相鄰控制體之間的影響,一旦每個控制體的變化得到,整個流場的變化也就知道了。 所以,再強調一次,關鍵問題是計算控制體界面上的。 初學者會說,這個不難,把界面上的插值得到,然后就可以計算。有道理! 咱們畫個圖,有三個小控制體 i-1到i+1,中間的“|”表示界面,控制體i右邊的界面用表示,左邊的就是。 | i-1 | i | i+1 | 好下個問題:每個小控制體長度都是如何插值計算界面上的? 最自然的想法就是:取兩邊的平均值唄, 但是很不幸,這是不行的。 那么換個方法?直接平均得到? 還是很不行,這樣也不行。 我靠,這是為什么?這明明是符合微積分里面的知識?。? 這個道理有點復雜,說開了去可以講一本書,可以說從50年代到70年代,CFD科學家就在琢磨這個問題。這里,初學者只需要記住這個結論:對于流動問題,不可以這樣簡單取平均值來插值或者差分。如果你非要想知道這個究竟,我現在也不想給你講清楚,因為我眼下的目的是讓你快速上手,而且該不刨根問底的時候就不要刨根問底,這也是初學階段一種重要的學習方法。 好了,既然目的只是為了求,我在這里,只告訴你一種計算方法,也是非常重要、非常流行的一種方法。簡單的說,就是假設流動狀態(tài)在界面是不連續(xù)的,先計算出界面兩邊的值,和,再由它們用某種方法計算出。上述方法是非常重要的,是由一個蘇聯人Godunov在50年代首創(chuàng)的,后來被發(fā)展成為通用Godunov方法,著名的ENO/WENO就是其中的一種。 好了,現在的問題是: 1 怎么確定和 2 怎么計算 對于第一個問題,Godunov在他的論文中,是假設每個控制體中是均勻分布的,因此 第二個問題,Godunov采用了精確黎曼解來計算。什么是“精確黎曼解”,就是計算這個激波管問題的精確解。既然有精確解,那還費功夫搞這些FVM算法干什么?因為只有這種簡單一維問題有精確解,稍微復雜一點就不行了。精確解也比較麻煩,要分析5種情況,用牛頓法迭代求解(牛頓法是什么?看數值計算的書去,哦,算了,現在暫時可以不必看)。 這是最初Godunov的方法,后來在這個思想的基礎上,各種變體都出來了。也不過是在這兩個問題上做文章,怎么確定,怎么計算。 Godunov假設的是每個小控制體內是均勻分布,也就是所謂分段常數(piecewise constant),所以后來有分段線性(picewise linear)或者分段二次分布(picewise parabolic),到后來ENO/WENO出來,那這個假設的多項式次數就繼續(xù)往上走了。都是用多項式近似的,這是數值計算中的一個強大工具,你可以在很多地方看到這種近似。 Godunov用的是精確黎曼解,太復雜太慢,也不必要,所以后來就有各種近似黎曼解,最有名的是Roe求解器、HLL求解器和Osher求解器,都是對精確黎曼解做的簡化。 這個多項式的階數是和計算精度密切相關的,階數越高,誤差就越小。不過一般來說,分段線性就能得到不錯的結果了,所以工程中都是用這個,Fluent、Fastran以及NASA的CFL3D、OverFlow都是用這個。而黎曼求解器對精度的影響不是那么大,但是對整個算法的物理適用性有影響,也就是說某種近似黎曼求解器可能對某些流動問題不合適,比如單純的Roe對于鈍頭體的脫體激波會算得亂七八糟,后來加了熵修正才算搞定。 上次(http://gezhi.org/node/399)說到了求解可壓縮流動的一個重要算法,通用Godunov方法。其兩個主要步驟就是 1 怎么確定和 2 怎么計算 這里我們給出第一點一個具體的實現方法,就是基于原始變量的MUSCL格式(以下簡稱MUSCL)。它是一種很簡單的格式,而且具有足夠的精度,NASA著名的CFL3D軟件就是使用了這個格式,大家可以去它的主頁(http://cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html)上看手冊,里面空間離散那一章清楚的寫著。 MUSCL假設控制體內原始變量(就是)的分布是一次或者二次多項式,如果得到了這個多項式,就可以求出控制體左右兩個界面的一側的值和。 我們以壓力p為例來說明怎么構造這個多項式。這里我只針對二次多項式來講解,你看完之后肯定能自己推導出一次多項式的結果(如果你搞不定,那我對你的智商表示懷疑)。 OK,開始 假設,這個假設不影響最終結論,因為你總可以把一個區(qū)間線性的變換到長度為1的區(qū)間。 假設壓力p在控制體i內部的分布是一個二次多項式,控制體i的中心處于處,左右兩個界面就是和。 這里先強調一個問題,在FVM中,每個控制體內的求解出來的變量實際上是這個控制體內的平均值。 所以, 。 我們知道,和,等距網格情況下和處的導數可以近似表示為 那么 (這里錯了,應該是2ax+b) 由上述三個有關a,b和c的方程,我們可以得到 這樣就可以得到f(x)的表達式了,由此可以算出和 通常MUSCL格式寫成如下形式 對應我們的推導結果(二次多項式假設)。 但是這不是最終形式。如果直接用這個公式,就會導致流場在激波(間斷)附近的振蕩。因為直接用二次多項式去逼近一個間斷,會導致這樣的效果。所以科學家們提出要對間斷附近的斜率有所限制,因此引入了一個非常重要的修改——斜率限制器。加入斜率限制器后,上述公式就有點變化。 這里是Van Albada限制器 是一個小數(),以防止分母為0。 密度和速度通過同樣的方法來搞定。 密度、速度和壓力被稱作原始變量,所以上述方法是基于原始變量的MUSCL。此外還有基于特征變量的MUSCL,要復雜一點,但是被認為適合更高精度的格式。然而一般計算中,基于原始變量的MUSCL由于具有足夠的精度、簡單的形式和較低的代價而被廣泛使用。 OK,搞定了。下面進入第二點,怎么求。關于這一點,我不打算做詳細介紹了,直接使用現有的近似黎曼解就可以了,都是通過和計算得到。比如Roe因為形式簡單,而非常流行。在CFL3D軟件主頁(http://cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html)上看手冊,附錄C的C.1.3。 想了一下,還是把Roe求解器稍微說說吧,力求比較完整。但是不要指望我把Roe求解器解釋清楚,因為這個不是很容易三言兩語說清的。 Roe求解器的數學形式是這樣的 顯然這個公式的第一項是一個中心差分形式,先前說過簡單的中心差分不可行,原因是耗散不足導致振蕩,說得通俗點就像一個彈簧,如果缺乏耗散(阻尼)它就會一直振蕩?!昂纳ⅰ边@個術語在激波捕捉格式中是最常見的。第二項的作用就是提供足夠的耗散了。 這里和已經用MUSCL求得了,的定義在第一講中已經介紹了。只有是還沒說過的。 這個矩陣可以寫成特征矩陣和特征向量矩陣的形式 而 ,,和的具體表達式在許多書上都有,而且這里的矩陣表達有問題,所以就不寫了。 是由、和代入計算得到。而、和采用所謂Roe平均值 這才是Roe求解器關鍵的地方! 總結一下,就是用Roe平均計算界面上的氣體狀態(tài),然后計算得到,這樣就可以得到了。如果有時間,我后面會找一個代碼逐句分析一下。 總之,計算還是很不直接的。構造近似黎曼解是挺有學問的,需要對氣體動力學的物理和數學方面有較深的理解。通常,如果不是做基礎研究,你只需要知道它們的特點,會用它們就可以了,而不必深究它們怎么推導出來的。 附錄程序:(新建一個.c類型文件,將下面的程序復制粘貼到里面,就可以運行了) /**************************************************** 本程序用于求解一維無粘可壓縮歐拉方程(激波管問題) 運用Dummy Cell處理邊界條件; 通量計算方式: AUSM Scheme; 重構方法:MUSCL方法 限制器:Van Albada限制器 時間離散:四步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) //網格步長 #define Nc 404 //網格總數:與h之間的關系 Nc=1/h+4 #define PI 3.1415927 #define It 1000 #define gama 1.4 //氣體比熱比 double KAKA=0.0;//限制器控制參數 double XS1,XS2;//計算域的兩個端點 double dt=2.5e-5;//時間步長 double timesum;//總的計算時間 //函數聲明 void output(); void SolveWtoU(double W[3],double U[3]); void SolveUtoW(double W[3],double U[3]);//基本變量和守恒變量之間的轉換函數 //前后各留兩個網格單元作為虛擬網格單元 計算網格單元從2到NOc-3 struct cell { int flag;//網格點的類型 double xc; double W[3],Wp[3];//conservation varaible double U[3];//jibenbianliang double R[3]; double S;//entropy }; struct cell cell[Nc]; //網格生成及流場初始化 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;i- 配套講稿:
如PPT文件的首頁顯示word圖標,表示該PPT已包含配套word講稿。雙擊word圖標可打開word文檔。
- 特殊限制:
部分文檔作品中含有的國旗、國徽等圖片,僅作為作品整體效果示例展示,禁止商用。設計者僅對作品中獨創(chuàng)性部分享有著作權。
- 關 鍵 詞:
- CFD 理論 過渡 編程 傻瓜 入門教程
裝配圖網所有資源均是用戶自行上傳分享,僅供網友學習交流,未經上傳用戶書面授權,請勿作他用。
鏈接地址:http://www.3dchina-expo.com/p-1544534.html