熱完全空氣模型及其應(yīng)用

熱完全空氣模型及其應(yīng)用

真實(shí)氣體效應(yīng)的影響空天飛機(jī)的動(dòng)態(tài)結(jié)構(gòu)是高超速飛機(jī)研究的第一個(gè)和核心技術(shù)。所有國(guó)家都非常重視這方面的研究。如美國(guó)的Trailblazer計(jì)劃、Hyper-X計(jì)劃、Loflyte計(jì)劃,俄羅斯CIAM與美國(guó)NASA的合作計(jì)劃,澳大利亞的Hyshot計(jì)劃等均十分重視高超聲速條件下飛行器的氣動(dòng)構(gòu)型研究。高超聲速飛行器氣動(dòng)問(wèn)題的研究最基本的問(wèn)題是高超聲速時(shí)空氣的氣動(dòng)熱力學(xué)特性發(fā)生了很大的變化。當(dāng)高超聲速氣流通過(guò)激波壓縮或粘性阻滯而減速時(shí),部分氣體將運(yùn)動(dòng)的動(dòng)能轉(zhuǎn)化為分子隨機(jī)運(yùn)動(dòng)的能量,從而使氣體達(dá)到很高的溫度。在這種極高的溫度下,非惰性氣體分子的振動(dòng)能和電子能將被激發(fā),甚至氣體之間會(huì)發(fā)生碰撞離解反應(yīng)、置換反應(yīng)、碰撞電離反應(yīng)等化學(xué)反應(yīng)。高超聲速飛行器的流場(chǎng)和氣動(dòng)特性(如升力、阻力、力矩等參數(shù))均受到真實(shí)氣體效應(yīng)的顯著影響。例如:阿波羅(Apollo)的飛行實(shí)驗(yàn)結(jié)果表明,其指揮艙的配平攻角比風(fēng)洞實(shí)驗(yàn)的預(yù)測(cè)值在高馬赫數(shù)時(shí)要大2～4°。而據(jù)最新報(bào)道,美國(guó)用無(wú)化學(xué)反應(yīng)氣體和有化學(xué)反應(yīng)氣體兩種模型分別計(jì)算了Apollo指揮艙的配平攻角,其差別和風(fēng)洞預(yù)測(cè)結(jié)果與飛行實(shí)驗(yàn)結(jié)果的差別一致,從而說(shuō)明了研究真實(shí)氣體效應(yīng)流動(dòng)的重要性。為了用NAPA軟件分析高超聲速飛行器氣動(dòng)問(wèn)題,本文首先用五次多項(xiàng)式擬合給出了溫度在50～3000K范圍內(nèi)的熱完全空氣的焓值與溫度之間的函數(shù)關(guān)系式,推出了其它熱力參數(shù)的內(nèi)能e,定壓比熱cp,定容比熱cv及比熱比γ的表達(dá)式。接著提出了熱完全空氣總溫、總壓的計(jì)算方法,并將其計(jì)算結(jié)果與量熱完全空氣的結(jié)果作了比較。將該熱完全空氣模型用于N-S方程求解,對(duì)三維粘性流場(chǎng)計(jì)算軟件(NAPA)進(jìn)行了改進(jìn),利用改進(jìn)后的軟件對(duì)高馬赫數(shù)鈍體繞流流場(chǎng)和乘波體流場(chǎng)進(jìn)行了計(jì)算。1熱完全空氣的概念常溫下氣體的熱力參數(shù)定壓比熱cp、定容比熱cv和比熱比γ等均為常數(shù)。但是,隨著馬赫數(shù)和溫度的提高,當(dāng)T>800K時(shí)空氣分子的振動(dòng)自由度的激發(fā),此時(shí)空氣的比熱僅為溫度的函數(shù),這種空氣稱為熱完全空氣;當(dāng)T>2500K時(shí)會(huì)發(fā)生化學(xué)反應(yīng),此時(shí)空氣的比熱不僅是溫度還是壓強(qiáng)的函數(shù)。1.1空氣力學(xué)參數(shù)的計(jì)算由于不考慮化學(xué)反應(yīng),空氣可看作是一種雙原子分子氣體。考慮到熱完全空氣的振動(dòng)自由度被激發(fā),于是根據(jù)統(tǒng)計(jì)熱力學(xué)和近代量子力學(xué)關(guān)于氣體熱力性質(zhì)與溫度的理論,單位質(zhì)量空氣的內(nèi)能由平動(dòng)能、轉(zhuǎn)動(dòng)能、振動(dòng)能組成,其表達(dá)式如下e=32RΤ+RΤ+hν(kΤ)ehν/(kΤ)-1RΤ(1)再根據(jù)焓與內(nèi)能的關(guān)系h=e+p/ρ=e+RT(2)可以得到單位質(zhì)量空氣的焓的表達(dá)式h=32RΤ+RΤ+hν/(kΤ)ehν/(kΤ)-1RΤ+RΤ(3)利用式(1～3)可以得到空氣的熱力學(xué)參數(shù)h,e,cp,cv和γ與溫度的理論函數(shù)關(guān)系式。但是由于式(1,3)中含有指數(shù)運(yùn)算和分式運(yùn)算,這使得求導(dǎo)計(jì)算非常復(fù)雜,會(huì)大大降低程序的運(yùn)行效率,增加計(jì)算機(jī)時(shí)?？紤]到指數(shù)和分?jǐn)?shù)函數(shù)可以展開(kāi)成冪級(jí)數(shù)多項(xiàng)式的形式,而多項(xiàng)式的積分和微分運(yùn)算簡(jiǎn)單,易用程序?qū)崿F(xiàn)。因此,可以采用溫度的冪級(jí)數(shù)多項(xiàng)式來(lái)描述空氣的熱力學(xué)參數(shù)。本文采用溫度的5次冪級(jí)數(shù)多項(xiàng)式來(lái)近似地表示空氣的焓與溫度的函數(shù)關(guān)系h=B0+B1T+B2T2+B3T3+B4T4+B5T5(4)根據(jù)文給出的50～3000K溫度范圍內(nèi)的空氣熱力學(xué)參數(shù)表,采用最小二乘法曲線擬合的方法,得到了式(4)的系數(shù)B0=1.817160×103B1=9.890504×102B2=-9.595592×10-3B3=1.041469×10-4B4=-4.433065×10-8B5=5.879263×10-12圖1給出了文的焓值與式(4)計(jì)算值的誤差。由圖可以看出在200～3000K溫度范圍內(nèi)擬合值與表中數(shù)據(jù)的誤差均在0%～0.5%之間,在低溫范圍內(nèi)(50～200K)誤差較大,但仍在5%之內(nèi)。1.2熱傳導(dǎo)系數(shù)的修正根據(jù)熱力學(xué)關(guān)系式,由式(4)可以求出單位質(zhì)量熱完全空氣的其它熱力學(xué)參數(shù)的內(nèi)能e、定壓比熱cp、定容比熱cv和比熱比γ與溫度的函數(shù)關(guān)系e=h-RΤ=B0+(B1-R)Τ+B2Τ2+B3Τ3+B4Τ4+B5Τ5(5)cp=dhdΤ=B1+2B2Τ+3B3Τ2+4B4Τ3+5B5Τ4(6)cv=dedΤ=cp-R=(B1-R)+2B2Τ+3B3Τ2+4B4Τ3+5B5Τ4(7)γ=cpcv=B1+2B2Τ+3B3Τ2+4B4Τ3+5B5Τ4(B1-R)+2B2Τ+3B3Τ2+4B4Τ3+5B5Τ4(8)同樣利用式(6,8)計(jì)算出各個(gè)溫度點(diǎn)下的定壓比熱cp和比熱比γ并與表中的數(shù)據(jù)進(jìn)行比較,結(jié)果表明在200～2800K溫度范圍內(nèi),計(jì)算誤差均在0%～1%之間,在低溫范圍內(nèi)(50～200K)誤差在4%之內(nèi)。另外,在數(shù)值求解N-S方程時(shí),還需知道另外兩個(gè)重要的參數(shù)即普朗特?cái)?shù)Pr和粘性系數(shù)μ。在T<3500K時(shí),空氣的粘性系數(shù)μ可以采用Sutherland公式來(lái)計(jì)算,即μ=1.4587×10-6Τ3/2Τ+110.4(9)而普朗特?cái)?shù)Pr的定義式為Ρr=cpμk(10)式中k為熱傳導(dǎo)系數(shù)。本文采用文給出的公式來(lái)計(jì)算k值k=kk0k0(11)k0=1.994×10-3Τ3/2Τ+112(12)式中k/k0為修正系數(shù)。文僅給出了一些離散溫度點(diǎn)上的k/k0值。為了便于在計(jì)算程序中使用,本文仍采用冪級(jí)數(shù)多項(xiàng)式最小二乘法擬合的方法對(duì)文中給出的500～3000K溫度范圍內(nèi)的k/k0值進(jìn)行擬合(圖2),得出了如下四次多項(xiàng)式擬合式kk0=A0Τ+A1Τ1+A2Τ2+A3Τ3+A4Τ4A0=1.2173336A1=-1.0664188×10-3A2=1.8580564×10-6A3=-1.1454078×10-9A4=2.3933339×10-13(13)2發(fā)動(dòng)機(jī)的材料選擇氣流的總溫、總壓計(jì)算對(duì)流場(chǎng)的氣動(dòng)特性、進(jìn)氣道的性能(如總壓恢復(fù)系數(shù)等)以及發(fā)動(dòng)機(jī)的強(qiáng)度和材料的選擇具有重要的意義。但是,對(duì)于熱完全空氣,由于焓h、定壓比熱cp、定容比熱cv和比熱比γ等均是溫度的函數(shù),不能像量熱完全空氣那樣將氣流總溫、總壓與馬赫數(shù)之間的關(guān)系用一個(gè)簡(jiǎn)單的解析函數(shù)表達(dá)式來(lái)描述。2.1和h=t#tcd對(duì)于絕熱流動(dòng),熱力學(xué)第一定律為h*=h+V22(14)式中h*,h分別為熱完全空氣的總焓和靜焓。由于熱完全空氣的cp和γ是溫度的函數(shù),因此h*和h之差必須通過(guò)積分來(lái)求出h*-h=∫T*TcpdT(15)文是利用上式并通過(guò)引入平均比熱和平均比熱比等平均的概念來(lái)求解熱完全空氣的總溫。這種算法要反復(fù)進(jìn)行積分,因此計(jì)算起來(lái)相當(dāng)復(fù)雜,不利于在CFD軟件中應(yīng)用。于是本文采用下面的方法來(lái)求總溫。由于靜溫T、靜壓P和馬赫數(shù)Ma是已知的,利用式(4)可以求出靜焓h。而V=Μa√γRΤ(16)將V和h代入式(14)則可求出總焓h*,最后利用牛頓迭代法由式(4)就可以求出T*。2.2滯止壓力和靜壓關(guān)系根據(jù)熱力學(xué)第一定律和第二定律Τds=dh-1ρdp(17)對(duì)于等熵過(guò)程有dpp=dhRΤ=cpRΤdΤ(18)積分上式前,引入相對(duì)壓強(qiáng)函數(shù)π0。定義如下lnπ0=1R∫ΤΤ01ΤcpdΤ(19)將式(6)代入式(19),得到了相對(duì)壓力π0的計(jì)算式為lnπ0=1R(B1lnΤ+2B2Τ+32B3Τ2+43B4Τ3+54B5Τ4+B6)(20)式中B6為待定常數(shù)。但由于本文所關(guān)心的是總壓與靜壓的比值,因此B6可以取任意值。熱完全氣體的等熵壓比關(guān)系為ln(p2p1)=ln(π20π10)=lnπ20-lnπ10(21)于是,利用上一節(jié)中計(jì)算出來(lái)的總溫度,通過(guò)式(20,21)可以算出熱完全空氣的滯止壓力與靜壓的比值。圖3,4分別給出了給定氣流靜溫和不同馬赫數(shù)的熱完全空氣的總溫、總壓及與對(duì)應(yīng)的量熱完全空氣的總溫、總壓的比較。圖中方框?yàn)榱繜嵬耆諝獾挠?jì)算結(jié)果,空心圓為熱完全空氣的計(jì)算結(jié)果,三角為靜溫。由圖可以看出,當(dāng)Ma數(shù)較低時(shí),熱完全空氣與量熱完全空氣的滯止參數(shù)差別很小;隨著Ma數(shù)的提高,空氣的振動(dòng)能、電子能被激發(fā),真實(shí)氣體效應(yīng)越來(lái)越顯著,于是兩種氣體的滯止參數(shù)產(chǎn)生了很大的差別。3時(shí)間/空間方向算法本文利用上述熱完全氣體的氣動(dòng)參數(shù)計(jì)算式對(duì)已有的NAPA軟件進(jìn)行了改進(jìn),用N-S方程進(jìn)行數(shù)值求解,在時(shí)間方向上采用了時(shí)間相關(guān)法和五步Runge-Kutta法推進(jìn)求解,在空間方向采用有限體積法離散,空間通量計(jì)算采用2階精度的AUSM+格式。為了減少計(jì)算機(jī)時(shí),還采用了多項(xiàng)加速收斂技術(shù)如多層網(wǎng)格、當(dāng)?shù)貢r(shí)間步長(zhǎng)和殘值光順等。3.1熱完全空氣實(shí)驗(yàn)高超聲速鈍體繞流在宇航飛行器的設(shè)計(jì)研究中有重要的意義。因此,本文首先采用熱完全空氣(變比熱)模型對(duì)一個(gè)圓柱繞流流場(chǎng)進(jìn)行了數(shù)值計(jì)算,然后與量熱完全空氣(定比熱)的計(jì)算結(jié)果進(jìn)行了比較。圓柱體半徑為r=0.1m,圓柱體前方自由來(lái)流的條件為Μa∞=8.03ρ∞=0.4135kg/m3p∞=26500ΡaΤ∞=223.25ΚRe=68.187×106/m圖5給出了圓柱繞流流場(chǎng)壓強(qiáng)分布的等值線圖。由圖可以看出,采用熱完全空氣模型計(jì)算出來(lái)的脫體激波(圖5A)比量熱完全空氣情況下的激波更貼近圓柱表面(圖5B)。這主要是由于熱完全空氣經(jīng)激波壓縮后的密度比量熱完全空氣的值高,流動(dòng)所需的面積變小,因此熱完全空氣的激波脫體距離比量熱完全空氣的小。這一點(diǎn)與文中給出的高超聲速圓柱繞流實(shí)驗(yàn)與理想氣體計(jì)算結(jié)果比較圖一致。圖6為沿圓柱固壁表面溫度、密度和壓比的分布曲線,圖中p0為來(lái)流靜壓強(qiáng)。從圖6(a)可以看出壁面溫度很高,氣動(dòng)加熱現(xiàn)象很明顯,但熱完全空氣模型的溫度值明顯低于量熱完全空氣的結(jié)果。圓柱前緣駐點(diǎn)溫度在熱完全空氣情況下為:Tb=2887.0K,而量熱完全空氣情況下Td=3368.5K,兩者相差14.3%。這主要是由于熱完全空氣經(jīng)過(guò)激波壓縮時(shí)受真實(shí)氣體效應(yīng)的影響,波前的動(dòng)能除了像量熱完全空氣那樣轉(zhuǎn)換成波后空氣分子的平動(dòng)能和轉(zhuǎn)動(dòng)能以外,還有一部分能量將被空氣分子振動(dòng)能和電子能的激發(fā)等過(guò)程吸收,因此它的溫度要明顯低于量熱完全空氣。從圖6(b)則可以看出,熱完全空氣壁面上的密度值要比量熱完全空氣情況下的結(jié)果高。如前所述,這是造成熱完全空氣的激波更貼近圓柱表面的原因之一。圖6(c)表明,熱完全空氣壁面上的壓強(qiáng)同來(lái)流壓強(qiáng)的比值與量熱完全空氣的差別不大,這是由于壓強(qiáng)主要取決于流體動(dòng)力學(xué)過(guò)程,受熱力學(xué)過(guò)程的影響較小的緣故。3.2雙波場(chǎng)激發(fā)對(duì)比本文選用了文中的錐型流乘波體。其表面形狀可用式(22)描述。θ/δ=1-0.1cos2?+(0.39+0.1cos2φ)×(tan(π-?)/2.75)7.96(22)生成其下表面的形狀。乘波體前方自由來(lái)流條件為Μa=5ρ∞=0.1755kg/m3Τ∞=288.9Κp∞=14552.64ΡaRe=16.67×106/m圖8給出了乘波體下四個(gè)站位上的激波,可以看出在乘波體的下表面覆蓋著一道很強(qiáng)的錐形激波。計(jì)算結(jié)果表明,量熱完全空氣和熱完全空氣兩種情況下激波的位置并沒(méi)有明顯的差別。這主要是由于乘波體是根據(jù)所需的激波形狀和位置通過(guò)流場(chǎng)反設(shè)計(jì)出來(lái)的,因此其激波的形狀和位置主要由其下表面形狀所決定,流體的熱力性質(zhì)對(duì)