熱噴干擾氣體模型對飛行器氣動特性影響分析1)

傅楊奧驍 劉慶宗 丁明松 江 濤 李 鵬 董維中 許 勇 高鐵鎖

(中國空氣動力研究與發展中心,計算空氣動力研究所,四川綿陽 621000)

引言

噴流反作用控制系統(reaction control system,RCS)是目前很多高超聲速飛行器進行跨流域飛行時所采用的氣動控制方法,相比于常規的氣動舵面控制技術,具有響應快、不受空域限制的優點[1-2],可明顯提升高超聲速飛行器跨流域飛行的控制和機動能力.

當噴流反作用控制系統工作時,發動機產生的噴流將與來流發生強烈的相互作用,產生一系列復雜的流場結構,這種現象被稱為噴流干擾效應[3-4],會對飛行器氣動力熱特性產生很大影響.在實際工程應用中,噴流反作用控制系統一般采用小型火箭發動機,噴管出口為高溫燃氣,一般稱之為“熱噴”,與之相對的是在風洞試驗中,出于成本和可實現性的考慮,常常采用空氣作為噴流氣體進行試驗,由于其溫度一般較低,因此被稱為“冷噴”.噴流燃氣的氣體性質與常溫空氣有較大差別,同時,噴流燃氣中未完全燃燒的組分進入主流后還會發生二次燃燒現象[5-6],使噴流干擾問題更加復雜,這些現象可統稱為熱噴干擾效應.

目前,已有一些學者針對熱噴干擾效應開展了相關研究.在這些研究中,很大一部分采用了二元等效比熱比異質流模型進行計算:如文獻[7]采用二元等效比熱比異質流模型,研究了尖錐外形的熱噴干擾問題;Shingo 等[8]針對火星著陸器外形,采用此類模型研究了著陸器RCS 噴流干擾對氣動力特性的影響;趙弘睿等[9]針對攔截彈外形,采用類似方法研究了高空環境下的側向噴流干擾問題.還有部分學者采用化學反應凍結流對熱噴干擾問題進行了研究:如Votta 等[10]采用化學反應凍結流數值模擬了火星著陸器的噴流干擾問題;DeSpirito[11]采用此模型研究了熱噴干擾問題中湍流模型的影響.由于在風洞中采用空氣進行冷噴試驗成本低、重復性好、易于實現,部分學者嘗試采用空氣噴流來模擬燃氣熱噴流,如林薄希等[12]采用空氣完全氣體模型,研究了噴流干擾氣動熱數值模擬的相似準則問題.以上研究中采用的多種氣體模型,均對實際的燃氣熱噴干擾問題進行了不同程度的簡化,可以明顯提高計算效率,但在一定條件下,同樣可能引入誤差,不能準確反映真實的物理情況.因此,一些學者采用了較為精細的氣體模型進行研究,如Li 等[13-14]采用化學反應流模型對火箭級間熱分離流場進行了數值模擬,結果顯示噴流燃氣的化學反應對噴流干擾區范圍有很大影響,但文中同樣指出,采用該種模型計算需要巨大的計算資源;Dong 等[15]采用類似模型對錐柱裙外形的熱噴干擾流場進行了數值模擬,結果顯示燃氣化學反應對分離區長度具有重要影響;趙法明等[16]基于類似模型對鈍錐外形的噴流干擾流場進行數值模擬,結果顯示熱噴燃燒效應會引起噴流干擾流場結構變化,進而影響飛行器的氣動力特性.

可見,目前對于熱噴干擾問題的研究,從提高計算效率的角度出發而采用各種簡化模型,可以解決大規模計算的問題,但在某些條件下可能引入誤差;從更加全面地考慮問題的角度出發而采用精細模型計算時,其計算效率往往很低.在各種條件下,綜合考慮計算效率和計算精度而對氣體模型做出選擇,是工程應用中所面臨的問題,因此有必要系統地考察不同氣體模型對熱噴干擾流場數值模擬結果的影響,為實際工程應用提供參考.

本文針對典型外形的RCS 熱噴干擾問題,基于自主研發的氣動物理流場計算軟件AEROPHFlow 以及噴流燃氣物理化學模型,分別采用化學反應流、反應凍結流、二元異質流以及空氣噴流四種氣體模型進行了數值模擬,分析了不同氣體模型對熱噴干擾流場結構及飛行器氣動力熱特性的影響,同時考察了不同馬赫數、飛行高度下的變化規律.

1 數值計算方法

流場控制方程采用包含化學反應源項的三維Navier-Stokes 方程,其無量綱化形式如下[17]

式中Q是守恒變量,完全氣體模型Q=(ρ,ρu,ρv,ρw,ρE)T,化學反應流、反應凍結流及二元異質氣體模型.其中ρ是混合氣體總密度,ρi是組分i的密度,u,v,w為直角坐標下三個方向的速度,Et為氣體內能,Re是雷諾數,F,G,H和FV,GV,HV分別對應了三個方向的對流項和黏性項,W為化學反應源項.

采用結構網格的有限差分方法離散控制方程(1),對流項采用AUSMPW+格式離散,黏性項采用中心差分格式離散,湍流模型采用k-ωSST 模型[18],時間離散采用LU-SGS 隱式方法.具體處理方法詳見文獻[17,19].

2 氣體模型

2.1 化學反應流模型

RCS 發動機產生的高溫燃氣進入外流空氣流場后,首先需要考慮高溫下空氣組分發生的離解、復合和置換等化學反應,同時,由于燃氣中含有一些易燃組分,其進入空氣后還會發生一系列復雜的二次燃燒反應,而本文采用的化學反應流模型則考慮了上述復雜化學反應.

由于噴流燃氣中的液體相雜質占比很小,對干擾流場的影響很小,因此這里忽略燃氣噴流中少量液體相雜質的影響.采用的化學反應模型如表1 所示,其中包含了空氣以及燃氣組分在高溫下發生的離解、復合、置換等化學反應,表中M 代表碰撞體.各反應的反應速率通過Arrhenius 公式計算,系數常數取自文獻[20-25].

表1 化學反應模型Table 1 Chemical reaction model

流動過程中,化學反應的模擬通過流動方程與化學反應源項的強耦合自動實現,源項中i組分生成源項wi可寫為

式中,Nr為模型中化學反應個數,Mi和Qj分別為第i組分的分子量和第j化學反應的生成源項,和γij分別為第j反應中第i組分的生成物系數和反應物系數.

2.2 化學反應凍結流模型

本文采用的化學反應凍結流凍結了燃氣組分相關的化學反應,即將噴流燃氣視為不反應的多組分化學惰性氣體,只考慮來流空氣組分在高溫下發生的各種化學反應,以簡化問題并提高計算速率,這種處理方法在熱噴干擾研究中比較常見[10-11].

在計算中,該模型的輸運系數計算方法與化學反應流一致:混合氣體的輸運系數采用Wilke 半經驗公式計算,單個組分的輸運系數采用氣體動理論計算[25].該模型的狀態方程及能量關系式同樣與化學反應流一致,考慮了氣體在高溫下的熱力學激發,包括各組分的束縛電子激發效應和分子組分振動能量激發效應[26].具體處理方法詳見文獻[17,19,27].

2.3 二元異質流模型

二元異質流模型中,將來流空氣與RCS 噴流工質當作兩種具有不同分子量和比熱比的完全氣體,流場介質可看作兩組分量熱完全氣體混合物,它們之間無化學反應,但存在擴散現象.相比于化學反應凍結流,該模型的處理更加簡化,因此計算效率更高,在熱噴干擾研究中被廣泛應用[7-9].

在該模型中,來流空氣及噴流工質均采用了量熱完全氣體假設,輸運系數采用Sutherland 公式計算,氣體內能只是溫度的線性關系式,即

式中,e為氣體內能,Cv為氣體定容比熱,T為氣體溫度,γ為氣體比熱比,p為壓力,ρ為密度.本文中噴流工質比熱比為1.27,來流空氣為1.4.

2.4 空氣噴流模型

由于在風洞試驗中采用空氣噴流重復性好、容易實現,因此有許多學者試圖通過一些相似準則,以實現空氣噴流對燃氣噴流的替代模擬[12],這種模型同樣采用量熱完全氣體假設,并且不需要考慮噴流和來流氣體性質差異,因此可以極大提升計算效率.

本文采用與多數文獻中相同的做法[12,28-29],即在空氣、燃氣噴流的噴管出口形狀與面積保持一致的條件下,保證二者的壓力比、動量比和質量流量比一致,如下

式中,p為靜壓,ρ 為密度,V為速度,下標j表示噴管出口參數,air 代表空氣噴流,gas 代表燃氣噴流.

3 計算方法的驗證

這里采用文獻[30]中的側噴干擾風洞試驗,來對本文的計算方法進行驗證,計算外形如圖1 所示,該外形由頭部、身部和尾裙三部分組成,身部直徑D=90 mm,長度為3.2D,頭部長2.8D,尾部長3D,底部直徑為1.66D,噴管出口直徑4.6 mm,位于X/D=4.3 處.

圖1 試驗模型外形Fig.1 Configuration of the test model

該風洞試驗包括冷噴和熱噴兩組狀態,兩組狀態下的風洞來流均為空氣,來流靜壓0.923 bar (1 bar=1 × 105Pa),來流靜溫105 K,來流馬赫數3.0,攻角0°,噴流總壓120 bar,噴流壓比130.冷噴時噴流氣體為空氣,噴流總溫293 K,噴流出口溫度244 K;熱噴時噴流氣體為燃氣,噴流總溫2300 K,噴流出口溫度2058 K,噴流出口比熱比1.235,噴口處組分及質量分數為:CO2(38.21%),H2(1.73%),H2O(10.47%),CO(35.47%),N2(14.12%).

冷噴狀態采用完全氣體模型計算(比熱比1.4),熱噴狀態采用前述的化學反應流模型計算.壁面處取絕熱壁、無滑移、零壓力梯度、非催化條件.

圖2 給出了流場參數分布云圖(其中流場對稱面為Ma云圖,模型表面為Cp云圖),圖3 給出了計算得到的模型表面壓差系數分布與文獻中冷噴試驗數據的對比(ΔCP=CP,jet-CP,nojet),可以看出,本文的計算結果與試驗數據吻合較好,計算結果準確地捕捉到了模型表面流向和周向的流動分離位置,對表面整個噴流干擾區域的壓力系數分布預測也較為準確.

圖2 流場參數分布云圖Fig.2 Flow field parameters distribution contour

圖3 模型表面表面壓力系數分布對比Fig.3 Comparison of surface pressure coefficient distribution

4 計算結果與分析

4.1 研究對象及網格無關性驗證

本文的研究對象為典型的錐柱裙外形,計算網格如圖4 所示,該外形總長為2250 mm,身部直徑386 mm,底部直徑502 mm,RCS 系統噴管位于質心,距頭部1177 mm.

圖4 計算網格Fig.4 Sketch of computation grid

RCS 系統發動機為典型的小型火箭發動機,推進劑為N2O4/MMH,燃燒室總溫3200 K,總壓4 MPa,總推力為6800 N,噴管出口直徑78 mm,噴管出口馬赫數Ma=2.796,出口靜壓P=133.3 kPa,出口靜溫T=1596 K,出口燃氣比熱比y=1.27,氣體常數R=400 J/(kg·K),組分質量分數見表2.

表2 噴管出口燃氣組分質量分數Table 2 Hot gas species mass fraction at nozzle exit

為了驗證網格收斂性,分別采用稀疏(網格高度5 μm,結果用coarse 表示),中等加密(網格高度1 μm,結果用medium 表示),加密網格(網格高度0.5 μm,結果用dense 表示)三套網格進行對比計算,計算狀態為H=20 km,Ma=5,攻角a=0°的熱噴干擾流場.表3 給出了氣動力系數對比,其中CA為軸向力系數,CN為法向力系數,CMZ為俯仰力矩系數,圖5 給出了噴口前方高熱流區域的熱流分布對比.可以看出,三套網格得到的氣動力系數相差不大;稀疏網格得到的熱流峰值偏低,加密網格與中等加密網格的熱流分布差異較小,滿足網格無關性要求,為了兼顧計算效率和精度,本文采用中等加密網格進行計算.

圖5 不同網格得到的熱流分布對比Fig.5 Comparison of heat flux

表3 不同網格計算得到的氣動力系數對比Table 3 Comparison of aerodynamic characteristics computed by different grid

本文采用工程上常用的噴流干擾附加力、力矩以及放大因子來表征噴流控制效果.定義無噴時的氣動力Fjetoff,力矩Mjetoff,噴流開啟時的氣動力Fjeton,力矩Mjeton,則噴流干擾附加力Fji、附加力矩Mji為

在本文條件下,噴流沿法向噴出,因此定義推力放大因子

由于本文選取的計算外形,噴管中心軸線位于飛行器質心位置,因此采用文獻[31]中的方法定義力矩放大因子

式中,Fjet是噴管推力,D是彈體身部直徑.

在本文計算中,氣動力參考面積取飛行器最大截面積Sref=0.198 m2,參考長度lref=1.0 m,力矩參考點為質心位置Xref=1.177 m,Yref=0.0 m,Zref=0.0 m,俯仰力矩以抬頭為正.

4.2 氣體模型對熱噴干擾流場的影響

這里首先采用前述的化學反應流(reacting)、反應凍結流(froz)、二元異質流(binary)以及空氣噴流(air)模型對本文第3 節中的熱噴試驗狀態進行了數值模擬,并與試驗數據進行了對比.

圖6 給出了計算得到的模型表面壓差系數分布與文獻中熱噴試驗數據的對比,可見,化學反應流得到的表面壓力分布與試驗數據吻合較好,其余三種模型預測的分離區范圍均明顯小于試驗數據.這說明在該試驗條件下,為了反映真實物理環境,有必要更加全面地考慮干擾流場中的詳細化學反應機制,采用化學反應流模型進行計算.

圖6 模型表面壓力系數分布對比Fig.6 Comparison of surface pressure coefficient distribution

為了進一步分析氣體模型對熱噴干擾流場數值模擬的影響,這里針對前述的錐柱裙外形,選取H=20 km,Ma=5,a=0°狀態,開展了不同氣體模型的對比分析.在計算中,采用k-ωSST 湍流模型,對壁面處取等溫壁300 K、無滑移、非催化、零壓力梯度條件,噴口處取為噴管出口條件(表3).

表4 給出了不同模型得到的飛行器表面氣動力結果對比.從表中可以看出,在該計算狀態下,不同計算模型得到的軸向力系數差別較小,最大差別在1.3%以內,但不同模型得到的噴流附加推力和附加力矩差異較大,從而導致法向力系數和俯仰力矩系數差異較大.具體來看,當前計算狀態下的噴流干擾效應會產生向下的噴流附加推力以及附加低頭力矩,與化學反應流的結果相比,其余三種模型的附加推力明顯更小,附加力矩明顯更大,下面將分析產生這種現象的原因.

表4 不同計算模型得到的氣動力特性對比Table 4 Comparison of aerodynamic characteristics by different gas model

圖7 給出了飛行器表面不同位置的壓力分布對比,從圖中可以看出,化學反應流相比于其他三種模型,得到的表面噴流干擾區范圍增大,同時干擾區內壓力明顯升高,這是由于噴流燃氣中含有未反應完全的CO 和H2等組分,其進入主流后會與來流發生二次燃燒反應,如圖8 所示,在噴口附近區域,流場中的O2組分質量分數下降,CO2和H2O 組分質量分數升高,而CO2和H2O 正是二次燃燒反應的主要產物,由于這些反應均為放熱反應,其釋放的化學能使得噴流氣體進一步膨脹,進而導致噴流干擾區范圍變大,干擾區壓力升高.

表5 給出了飛行器的局部區域氣動力系數對比,表中下標up 表示彈體上表面的氣動力積分結果,下標down 表示下表面,upstream 表示噴口上游部分的氣動力積分結果,downstream 表示噴口下游部分.結合圖7 可知,化學反應流得到的噴流干擾區范圍變大,干擾區壓力升高,由于噴流干擾區主要位于彈體上表面,從而導致彈體上表面的氣動力積分結果CN,up絕對值增大,同時,由于噴流干擾效應對下表面壓力分布影響較小,各氣體模型得到的下表面積分結果CN,down相差不大,因此,整體的積分效果表現為化學反應流的法向力系數CN絕對值增大,噴流附加推力增大.從表5 還可以看出,噴口上下游部分的積分效果均表現為低頭力矩,但由于噴流干擾區大部分位于噴口下游,因此力矩主要由下游部分貢獻,從圖8(a)和8(b)可知,噴口下游上表面對稱線附近存在較強的二次燃燒反應,導致這一區域壓力明顯升高(如圖7(b)所示),這一部分的壓力積分表現為一個抬頭力矩,因此化學反應流得到的CMZ,downstream絕對值明顯偏小,從而導致其附加低頭力矩明顯偏小.

圖7 表面壓力分布對比Fig.7 Comparison of surface pressure distribution

圖8 噴口附近的組分質量分數分布對比Fig.8 Species mass fraction distribution near nozzle outlet

表5 飛行器局部區域氣動力系數對比Table 5 Comparison of part aerodynamic coefficient

從表4 和圖7 還可以看出,三種簡化模型之間的表面壓力分布和氣動力特性差異相對較小,說明在當前計算狀態下,氣體熱力學激發和噴流工質比熱比差別的影響相對較小.

圖9 給出了上表面對稱線上的熱流分布對比,可見三種簡化模型得到的熱流峰值差異不大,但與化學反應流的結果相比,其熱流峰值低約22%,結合圖8 可知,該區域是生成CO2和H2O 反應較劇烈的區域,二次燃燒反應釋放的化學能造成了該區域表面熱流升高,這也說明若采用簡化模型進行計算,將低估飛行器表面的熱流值,對防熱系統設計不利.

圖9 上表面對稱線上熱流分布對比Fig.9 Comparison of heat flux distribution along upper side symmetric line

圖10 給出了不同氣體模型的殘差收斂曲線和CPU 計算時間對比,綜合來看,采用空氣噴流模型的計算效率明顯最高,二元異質流次之,而由于采用了多組分氣體描述噴流燃氣,化學反應流和反應凍結流不僅收斂較慢、在同樣步數下的CPU 計算時間明顯更長,計算效率較低.

圖10 不同氣體模型的計算效率對比Fig.10 Comparison of efficiency for different model

總的來看,在當前計算條件下,采用簡化模型進行熱噴干擾流場的數值模擬時,飛行器的氣動力熱特性預測將出現一定偏差;若采用更加復雜的化學反應流模型進行計算,其計算效率較低.

更進一步來看,當飛行條件變化時,流場中化學反應和熱力學激發等物理化學現象的劇烈程度也會發生改變,在一些條件下,采用合適的氣體模型可以實現計算效率和精度的兼顧.因此,本文接下來將開展不同飛行條件下氣體模型對熱噴干擾流場數值模擬的影響分析.

圖11 給出了高度H=20 km,攻角a=0°,馬赫數Ma=4～ 10 狀態下的氣動力熱特性對比.可見,隨著馬赫數增大,化學反應流與簡化模型得到的結果差異逐漸增大,同時,在高馬赫數條件下,三種簡化模型相互之間的差異也明顯增大.以下從三個方面分析其原因:(1)隨著飛行馬赫數增大,來流的氧氣質量流量增大,使得二次燃燒反應更加充分,同時,高馬赫數下噴流與來流相互干擾作用更強,類似于激波誘導燃燒,在更強的激波作用影響下會誘導發生更強的二次燃燒反應,因此化學反應流與簡化模型得到的結果差異逐漸增大;(2)隨著馬赫數增大,噴流與來流相互干擾作用增強,造成噴流干擾區內溫度升高(如圖12 所示),由于凍結流考慮了多組分氣體內能各個模態的激發,在高溫狀態下,將更加偏離完全氣體狀態,此時雖然噴口處燃氣比熱比是一致的,但相比于低馬赫數狀態,高馬赫數下整個干擾區內氣體比熱比出現了明顯變化(如圖13 所示),而這將導致氣體熱力學特性發生變化,從而導致凍結流與二元異質流的差異增大;(3)在高溫條件下,由噴流工質的比熱比差別所導致的氣體膨脹特性差異將會更加明顯,因此二元異質流與空氣噴流的差異隨著飛行馬赫數增大而增大.

圖11 不同馬赫數下的氣動力熱特性對比(H=20 km,a=0°)Fig.11 Comparison of aerodynamic and aerothermal characteristics for different flight Mach number (H=20 km,a=0°)

圖12 反應凍結流溫度云圖對比(H=20 km,a=0°)Fig.12 Comparison of specific heat ratio contour for chemical frozen flow (H=20 km,a=0°)

圖13 反應凍結流比熱比云圖對比(H=20 km,a=0°)Fig.13 Comparison of specific heat ratio contour for chemical frozen flow (H=20 km,a=0°)

總的來看,隨著飛行馬赫數增加,簡化模型對熱噴干擾氣動力熱特性預估的誤差增大,同時不同簡化模型之間的差異也將增大,此時有必要更加全面地考慮問題的物理本質,采用化學反應流模型進行計算.

圖14 給出了高度H=20～ 50 km,馬赫數Ma=5,攻角0°狀態下的氣動力熱特性對比.可見,隨著飛行高度增加,化學反應流與簡化模型之間的差異迅速減小,在高度30 km 以上時,不同計算模型得到的飛行器氣動力熱特性差別已經較小.下面分析其原因:隨著飛行高度增加,來流密度迅速降低(從高度20 km 到50 km,來流密度下降了兩個量級),因此來流氧氣質量流量迅速降低,當來流氧氣質量流量過低時,由于缺乏氧化物,流場中二次燃燒反應將明顯減弱,因此化學反應流與簡化模型的差異明顯減小;同時,高空條件下,噴流與來流相互干擾作用將會減弱,噴流干擾區內溫度降低,此時氣體內能模態激發以及噴流工質比熱比差別的影響將減小,因此三種簡化模型相互之間的差別進一步減小.

圖14 不同高度下的氣動力熱特性對比(Ma=5,a=0°)Fig.14 Comparison of aerodynamic and aerothermal characteristics for different flight altitude (Ma=5,a=0°)

圖15 給出了高度H=20～ 50 km,馬赫數Ma=10,攻角0°狀態下的氣動力熱特性對比.可見,即使在較高馬赫數條件下,隨著高度增加,氣體模型的影響同樣迅速減小,說明此時飛行高度仍然是主要影響因素,類似地,在30 km 以上時,不同計算模型得到的飛行器氣動力熱特性差別已經較小.需要補充說明的是在Ma=10 狀態下,由于飛行速度很高、噴流干擾作用很強,在低高度條件下俯仰力矩出現了反向的現象(由低頭變為抬頭),因此圖15(b)中的變化趨勢與圖14(b)有所不同.

圖15 不同高度下的氣動力熱特性對比(Ma=10,a=0°)Fig.15 Comparison of aerodynamic and aerothermal characteristics for different flight altitude (Ma=10,a=0°)

總的來看,簡化模型對飛行器熱噴干擾氣動力熱特性預估的偏差,在低飛行高度狀態下較明顯,當飛行高度較高時,這種偏差明顯減小,此時可以采用簡化模型進行計算,在提高計算效率的同時也可以保持較高的預測精度.

5 結論

本文針對錐柱裙外形,分別采用化學反應流、反應凍結流、二元異質流以及空氣噴流四種氣體模型開展了RCS 熱噴干擾流場的數值模擬,開展了不同氣體模型對熱噴干擾流場結構、飛行器氣動力熱特性的影響特性和規律研究,得到以下結論.

(1)從與風洞試驗數據的對比來看,化學反應流模型與試驗數據的吻合程度優于三種簡化模型,要準確反映熱噴干擾流場特性,應盡量采用化學反應流模型進行計算.

(2)在本文條件下(高度H=20～ 50 km,Ma=4～10),飛行高度較低時(約30 km 以下),采用簡化模型進行熱噴干擾流場數值模擬,會低估分離區大小,使飛行器氣動力特性預測出現偏差,同時也會低估表面熱環境,對防熱系統設計不利;隨著馬赫數增加,這種偏差將會進一步增大,同時不同簡化模型之間的結果差異也進一步增大.

(3)飛行高度較高時(約30 km 以上),簡化模型對飛行器熱噴干擾氣動力熱特性預估的偏差明顯減小,此時可采用簡化模型進行計算,在保證精度的同時大幅提升計算效率.

本文選取的外形和計算狀態仍相對有限,在更高飛行高度及更高馬赫數狀態下,流場中會進一步出現離解/電離效應、熱化學非平衡效應、稀薄滑移效應等復雜的物理化學現象,因此在后續的研究中,有必要進一步探索多種外形、多種飛行狀態下熱噴干擾效應對飛行器氣動特性的影響特征及規律.