熱完全空氣模型及其應用

熱完全空氣模型及其應用

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