1、精選優質文檔-傾情為你奉上有限體積法求解二維可壓縮Euler方程計算流體力學課程大作業老師： 夏健、劉學強 學生： 徐錫虎 學號： SQ 日期： 2010年2月5日 目 錄一、內容摘要（2）二、流動控制方程（2）三、有限體積法的空間離散（2）四、人工耗散（3）五、時間離散（4）六、邊界條件（5）七、計算結果（8）八、結論與展望（11）參考文獻（11）一、內容摘要本文通過運用JAMESON有限體積法求解了二維定常和非定常可壓縮Euler方程。程序實現語言為C+。其中，使用的網格是三角形非結構網格。在時間推進上使用的是四步龍庫塔推進格式。推進的時間步長取的是當地的時間步長。為了消除迭代誤差、rou

2、nd-off等誤差，本文采用了添加人工耗散項的辦法。另外，本文計算了NACA0012翼型在跨音速下不同迎角的情況，并與fluent軟件的計算結果進行了比較，來驗證程序的準確性。二、流動控制方程守恒形式的Euler方程： （1）其中x和y代表笛卡兒坐標系。W是守恒變量。 (2)F,G表示通量, (3),P , H和E表示密度，壓強，單元總焓和單元總能量。U,V表示笛卡兒坐標系下的速度矢量。這些量由理想氣體的單位體積的總能量和總焓相互聯系。 (4) (5)三、有限體積法的空間離散 計算域被劃分為互不重疊的單元。在每個單元運用守恒形式的Euler方程。由于每個單元相對于時間都是不變的，所以等式（1）

3、可以寫成： （8） 其中和S是單元的體積和邊界。W是單元的平均值。在對上述方程進行時間離散前，先對空間進行離散，則方程（6）可以寫為： （9）其中表示第k個單元的體積，是第k個單元的守恒變量。表示第k個單元的通量。方程（7）的右邊項可以寫成： （10）其中 （11）（8）式中的求和是對第k個單元的所有邊進行的。守恒參數的量是單元中心值，在求通量時，第I條邊的守恒參數值是用左右單元的平均來表示的： （12）引入變量： （13）則第k單元的Euler方程可以寫為： （14）在本文中，采用的是JAMENSON有限體積法，為了減少存儲的相關信息的量，其存儲的方式選擇的是按邊存儲的方法。在存儲的每條邊的

4、信息中，包含了這條邊的邊號，左右單元號和邊的端點。在計算通量時采用按邊循環的方式： do I=1,nedge k=connmatrix(I,1) a=connmatrix(I,2) b=connmatrix(I,3) p=connmatrix(I,4) flux=function(k,a,b,p) sum(k)=sum(k)+flux sum(p)=sum(p)-flux end do這里給出的是FOTRAN語言的形式，我編寫采用的是C，具體表現在上交的程序中。在計算時間步長、人工耗散項等也可用象這樣按邊循環，從此處我們可以看出求解時與單元的形狀無關。四、人工耗散人工粘性模型對方法的成功應用起

5、著關鍵作用，人工粘性抑制解在激波附近的振蕩，又阻尼解在光滑區域內的高階誤差，對解的線性穩定和收斂于定態是很重要的。本文在方程（14）的右端加入了人工耗散項，如對于單元k，其表達式可以表示為： （15）在有限體積法中，耗散項的公式可以表示為： （16）其中： （17）其中的I表示單元k和p的公共邊，定義為： （18）上面的j表示與k相鄰的單元。 （19） （20）其中的量的范圍是：。 在計算時發現上面方法得到的人工耗散項并不太適合。其在光滑區域耗散項太大，而在大剃度區域又顯得太小，為了彌補上面的不足，作下面的修改： （21）自適應系數為： （22）尺度系數為： （23）其中的U,V表示邊上的值，

6、C表示當地聲速。五、時間離散方程最后的穩定解是通過時間上的迭代得到的，可以寫為： （24）右邊項的表達式為： （25）為了加速收斂，時間迭代使用的是4步龍庫塔推進格式。格式如下： （26） 其中的n表示的是當前的時間步，n+1表示的是新的時間步： （27） （28） 為了減小計算時間，人工耗散項的計算只在第一步進行，在下面幾步的迭代中保持不變。運用上面的方法計算，可以發現CFL數可以取到，本文中使用的是2.0。使用顯示格式迭代的主要缺點是由于穩定區域的限制，所以不能使用過大的時間步長。可以用近似的方法估算時間步長，對于任意形狀的網格，可以使用下面的方法： （29）六、邊界條件1 固面邊界條件對

7、與無粘流動，固面邊界條件無穿透條件，設其法向的速度通量為零，即。由于壓強項的影響，x向和y向的動量通量并不為零。固面的壓強近似的取為其相鄰的單元的單元中心壓強。廣泛的數值研究證明，如果貼近壁面的單元足夠小，并且人工耗散項運用正確，用這種方法取得的壓強對結果的精度不會產生太大的影響。2 遠場邊界條件本文提到的遠場，實際是人為的有界邊界，對于流場中的擾動會傳到很遠的地方。因而對于遠場邊界條件，情況比較復雜，它不能直接給定具體的流場值，需要與流場內的值來共同確定遠邊場的流場值。如果邊界取得過小，則通常采用環量修正。一般情況下，我們采用無反射邊界條件。為了保證擾動波不會反射回流場，應用A.Jamson

8、提出的遠場邊界法向一維特征分析方法，來建立無反射的遠場邊界條件。一維均熵流動的Euler方程可寫成： (30) (31)這里，動量方程中除去了壓強項，將上式寫成矩陣形式為： (32)這里：， (33)的特征值為：，因此，上式的兩族特征線為： (34) (35)(34)為特征線，(35)為特征線。沿，給出了眾所周知的Riemann不變量，： (36) (37)這里不變量，沿入流特征線是常數，可以用來流條件計算得到；沿出流特征線是常數，可以用流場內部向外插值計算： (38) (39)上式中下標“”表示來流值，下標“”表示以計算域內部參數外插獲得的值。通過Riemann不變量，的加減，可獲得遠場的法

9、向速度和音速： (40) (41)根據Riemann不變量，按邊界附近信息傳播的性質把遠場邊界條件分成以下四種情況：A、 亞音速入流條件，它有三條入流特征線，需規定三個條件： (42) (43) (44) (45)其中下標t表示切向，n代表方向，代表自由流，代表從流場內到邊界的外插值，上式的右端皆為已知，可解出邊界上的值，再由和求出和。B、 亞音速出流條件，它有一條入流特征線，需規定一個條件 (46) (47) (48) (49) C、 超音速出流條件，無入流特征線，不需在邊界上規定邊界條件 （50）其中的W表示邊上的守恒變量，表示與此邊相鄰元素的守恒變量值。 D、 超音速入流條件，它有四條入

10、流特征線，需規定全部四個條件 （51） 其中的W表示邊上的守恒變量，表示來流值。七、計算結果 本文計算了三個算例，一個是攻角，馬赫數0.80的情況，二是攻角，馬赫數0.80的情況，三是攻角2.5°，馬赫數1.5的情況 翼型網格示意圖 1 表面壓強系數分布 等壓線 等馬赫線 壓力云圖 馬赫數云圖升力系數 CL=-3.27408E-005阻力系數 CD= 1.E-0022. 上下表面壓強系數分布圖 等壓線圖 等馬赫圖 壓力云圖 馬赫數云圖升力系數為 CL2.E-001阻力系數為 CD=2.E-0023. 上下表面壓強系數分布 壓力云圖 馬赫數云圖升力系數為 CL= 0.阻力系數為 CD=

11、 0.八、結論與展望本文主要介紹了求解二維非結構網格歐拉方程定常解的問題。采用的是有限體積空間離散方法，單元形狀可以是任意多邊形。使用顯式多步推進過程，求解非定常方程得到定常解。一些標準的加速技術可以加快收斂速度。 Jameson格式是學習CFD編程的入門格式，其精度和收斂速度都不高，但通過實際編寫程序可以學到很多方法和技巧。因此，學習Jameson格式為以后編寫其他格式如Roe，AUSM等等以及求解三維方程和求解N-S方程都打下了良好的基礎。本次編程采用的是C語言，并且存在著一個比較大的問題，CFL數不能取到大于1.2的數值，在超音速計算時候CFL數能取的大值變的更小，估計是時間步長上的問題但是目前還沒有找到問題所在。忘老師指點。參考文獻1、L. STOLCIS* and L. J. JOHNSTON*, Solution of the Euler equations on unstructured grids for two-dimensional compressible flow, Aeron