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

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

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

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

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

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

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

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

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

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

10、流特征線，需規(guī)定全部四個(gè)條件 （51） 其中的W表示邊上的守恒變量，表示來(lái)流值。七、計(jì)算結(jié)果 本文計(jì)算了三個(gè)算例，一個(gè)是攻角，馬赫數(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ù)分布圖 等壓線圖 等馬赫?qǐng)D 壓力云圖 馬赫數(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)格歐拉方程定常解的問(wèn)題。采用的是有限體積空間離散方法，單元形狀可以是任意多邊形。使用顯式多步推進(jìn)過(guò)程，求解非定常方程得到定常解。一些標(biāo)準(zhǔn)的加速技術(shù)可以加快收斂速度。 Jameson格式是學(xué)習(xí)CFD編程的入門(mén)格式，其精度和收斂速度都不高，但通過(guò)實(shí)際編寫(xiě)程序可以學(xué)到很多方法和技巧。因此，學(xué)習(xí)Jameson格式為以后編寫(xiě)其他格式如Roe，AUSM等等以及求解三維方程和求解N-S方程都打下了良好的基礎(chǔ)。本次編程采用的是C語(yǔ)言，并且存在著一個(gè)比較大的問(wèn)題，CFL數(shù)不能取到大于1.2的數(shù)值，在超音速計(jì)算時(shí)候CFL數(shù)能取的大值變的更小，估計(jì)是時(shí)間步長(zhǎng)上的問(wèn)題但是目前還沒(méi)有找到問(wèn)題所在。忘老師指點(diǎn)。參考文獻(xiàn)1、L. STOLCIS* and L. J. JOHNSTON*, Solution of the Euler equations on unstructured grids for two-dimensional compressible flow, Aeron