用有限體積方法求解歐拉方程

1、精選優(yōu)質(zhì)文檔-傾情為你奉上有限體積法求解二維可壓縮Euler方程計算流體力學(xué)課程大作業(yè)老師： 夏健、劉學(xué)強(qiáng) 學(xué)生： 徐錫虎 學(xué)號： SQ 日期： 2010年2月5日 目 錄一、內(nèi)容摘要（2）二、流動控制方程（2）三、有限體積法的空間離散（2）四、人工耗散（3）五、時間離散（4）六、邊界條件（5）七、計算結(jié)果（8）八、結(jié)論與展望（11）參考文獻(xiàn)（11）一、內(nèi)容摘要本文通過運(yùn)用JAMESON有限體積法求解了二維定常和非定?？蓧嚎sEuler方程。程序?qū)崿F(xiàn)語言為C+。其中，使用的網(wǎng)格是三角形非結(jié)構(gòu)網(wǎng)格。在時間推進(jìn)上使用的是四步龍庫塔推進(jìn)格式。推進(jìn)的時間步長取的是當(dāng)?shù)氐臅r間步長。為了消除迭代誤差、rou

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

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

4、信息中，包含了這條邊的邊號，左右單元號和邊的端點。在計算通量時采用按邊循環(huán)的方式： 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，具體表現(xiàn)在上交的程序中。在計算時間步長、人工耗散項等也可用象這樣按邊循環(huán)，從此處我們可以看出求解時與單元的形狀無關(guān)。四、人工耗散人工粘性模型對方法的成功應(yīng)用起

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

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

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

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

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

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

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