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

PAGE1-有限體積法求解二維可壓縮Euler方程——計(jì)算流體力學(xué)課程大作業(yè)老師：夏健、劉學(xué)強(qiáng)學(xué)生：徐錫虎學(xué)號(hào)：SQ09018013018日期：2010年目錄一、內(nèi)容摘要……………………（2）二、流動(dòng)控制方程………………（2）三、有限體積法的空間離散……（2）四、人工耗散……………………（3）五、時(shí)間離散……………………（4）六、邊界條件……………………（5）七、計(jì)算結(jié)果……………………（8）八、結(jié)論與展望…………………（11）參考文獻(xiàn)…………（11）在有限體積法中，耗散項(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表示邊上的值，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ì)與無(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提出的遠(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)的特征值為：，因此，上式的兩族特征線(xiàn)為：(34)(35)(34)為特征線(xiàn)，(35)為特征線(xiàn)。沿，給出了眾所周知的Riemann不變量，：(36)(37)這里不變量，沿入流特征線(xiàn)是常數(shù)，可以用來(lái)流條件計(jì)算得到；沿出流特征線(xiàn)是常數(shù)，可以用流場(chǎng)內(nèi)部向外插值計(jì)算：(38)(39)上式中下標(biāo)“”表示來(lái)流值，下標(biāo)“”表示以計(jì)算域內(nèi)部參數(shù)外插獲得的值。通過(guò)Riemann不變量，的加減，可獲得遠(yuǎn)場(chǎng)的法向速度和音速：(40)(41)根據(jù)Riemann不變量，按邊界附近信息傳播的性質(zhì)把遠(yuǎn)場(chǎng)邊界條件分成以下四種情況：亞音速入流條件，它有三條入流特征線(xiàn)，需規(guī)定三個(gè)條件：(42)(43)(44)(45)其中下標(biāo)t表示切向，n代表方向，代表自由流，代表從流場(chǎng)內(nèi)到邊界的外插值，上式的右端皆為已知，可解出邊界上的值，再由和求出和。亞音速出流條件，它有一條入流特征線(xiàn)，需規(guī)定一個(gè)條件(46)(47)(48) (49)C、超音速出流條件，無(wú)入流特征線(xiàn)，不需在邊界上規(guī)定邊界條件（50）其中的W表示邊上的守恒變量，表示與此邊相鄰元素的守恒變量值。D、超音速入流條件，它有四條入流特征線(xiàn)，需規(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ù)分布等壓線(xiàn)等馬赫線(xiàn)壓力云圖馬赫數(shù)云圖升力系數(shù)CL=-3.27408E-005阻力系數(shù)CD=1.14998670E-0022.上下表面壓強(qiáng)系數(shù)分布圖等壓線(xiàn)圖等馬赫?qǐng)D壓力云圖馬赫數(shù)云圖升力系數(shù)為CL＝2.741130755E-001阻力系數(shù)為CD=2.15905376E-0023.上下表面壓強(qiáng)系數(shù)分布?jí)毫υ茍D馬赫數(shù)云圖升力系數(shù)為CL=0.1352582564阻力系數(shù)為CD=0.1037842242八、結(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*andL.J.JOHNSTON*,SolutionoftheEulerequationsonunstructuredgridsfortwo-d