何粲 邢建文,?,1) 歐陽浩 鄧維鑫 肖保國,?
* (中國空氣動力研究與發展中心空天技術研究所,四川綿陽 621000)
? (中國空氣動力研究與發展中心高超聲速沖壓發動機技術重點實驗室,四川綿陽 621000)
飛得更高更快一直是人類不斷追尋的目標,作為飛行器的心臟,發動機可以說是制約其速度的核心部件.多年來,吸氣式超燃沖壓發動機憑借其簡單構型、不需攜帶氧化劑等諸多優點一直是高超聲速飛行器的理想推進系統之一.目前國內針對超燃沖壓發動機的研究按照飛行馬赫數范圍可以分為兩部分,首先是飛行馬赫數4~ 7 的雙模態超燃沖壓發動機,其次是飛行馬赫數Ma≥7 的高馬赫數超燃沖壓發動機[1].
經過60 多年的研究,雙模態超燃沖壓發動機關鍵技術逐漸被突破,在飛行演示試驗中已能克服燃燒室的內阻和飛行器的外阻,產生正推力,即將進行工程應用[2-4].高馬赫數超燃沖壓發動機可用作大氣層內超高速、高機動飛行器和可重復使用、低成本空天運輸系統動力系統的一部分,有助于進一步提高高超聲速武器的突防能力,具有極高的軍事和民用價值.但與低馬赫數超燃沖壓發動機不同,高馬赫數超燃沖壓發動機仍然面臨諸多關鍵科學和工程技術難題.因此,超燃沖壓發動機下一步的重點研究方向之一是在現有超燃沖壓發動機的研究基礎上將其工作速域向上拓展,開展高馬赫數超燃沖壓發動機研究.
高馬赫數超燃沖壓發動機與馬赫數4~ 7 發動機相比有以下幾個特點:(1)燃燒室入口氣流速度更高,通常在1500 m/s 以上,更高的流動速度給燃料混合、點火和燃燒組織帶來了更大的難度;(2)阻力更大,更難實現推阻平衡,導致對混合和燃燒效率有更高的要求;(3)來流總溫更高,燃燒室的熱負荷急劇上升,發動機熱防護技術更加困難.這些特點使得高馬赫數超燃沖壓發動機的研究極具挑戰性,需要數值模擬、地面試驗以及飛行試驗的強有力結合才可能取得突破性進展.
以日本M12 系列及澳大利亞RESTM12 模型發動機[5]為代表,目前國外對高馬赫數的發動機已經開展了一系列試驗[6]與計算研究[7].其中M12 系列的研制目的是獲得設計點為Ma12,工作范圍為Ma10~ 15 的超燃沖壓發動機.M12 系列發動機均為2D 構型,其研制經歷了一個不斷優化改進的過程,相較于M12-01[8],M12-02 做了兩方面的改進[9].一是增加進氣道壓縮比,以提升燃燒室入口氣流的靜溫(提升至1400 K)和密度;二是應用hypermixer 進行燃料噴注,以增強燃氣摻混,改善發動機的穩焰特性和點火能力[10].
M12-02 開展了總焓為4,5,7,9 MJ/kg,分別對應飛行馬赫數9,10,12,14,4 種工況條件下的試驗.試驗結果表明M12-02可以獲得穩定且強烈的燃燒,燃燒帶來的壓升幅度增加且強燃燒時間延長,在高焓條件下的燃燒室性能相較于M12-01 得到顯著改善.燃燒室性能隨著氣流總焓增加而提高,在設計點7 MJ/kg 時獲得最佳性能.研究結果為改善高焓條件下超燃沖壓發動機燃燒室性能提供了參考與支撐.但需要注意的是,結果也顯示隨著來流速度的進一步增加,M12-02 燃燒室性能將急劇下降[10].同時CFD 預測在M12-02 燃燒室出口燃氣溫度超過2600 K 時,由燃氣熱離解導致的凈釋熱損失很顯著[11].
近年來,國內針對飛行馬赫數7 以上的高馬赫數發動機也開展了一些地面試驗及計算的研究,姚軒宇等[12]基于力學所JF12 長試驗時間激波風洞開展了Ma7.0和Ma9.5 的氫燃料點火和燃燒試驗.盧洪波等[13]在航天十一院FD-21 高能脈沖風洞完成了模擬Ma8,高度31 km 飛行條件的超燃試驗.周建興[14]等構造了一種圓截面高馬赫數發動機并評估了其在Ma7~ 10 內的性能.Zhang和Li[15]針對M12-02 發動機開展了計算與分析[16].
中國空氣動力研究與發展中心針對高馬赫數發動機開展了建模與計算研究[17],以及基于不同燃料的一系列點火和燃燒地面試驗[18-22].從以往研究可知,對于高馬赫數發動機內流動及燃燒特性徹底理解尚有距離,計算與試驗能力均有待進一步提升.本文將基于AHL3D 軟件平臺開展針對高馬赫數發動機的計算方法改進與驗證,并針對M12-02 發動機[23]不同狀態進行三維數值模擬,驗證改進后的方法對于該類高馬赫數發動機的模擬能力,分析高馬赫發動機內波系、參數分布及燃燒性能的特征.
本文基于AHL3D 軟件對發動機開展三維定常計算.AHL3D 軟件平臺可以模擬二維或三維、定常或非定常、完全氣體或化學非平衡流動[24],可以對激波邊界層干擾[25]、穩焰、燃燒[26]等復雜現象進行模擬,其對超燃沖壓發動機的計算能力得到了大量算例的驗證.這里簡單介紹使用的控制方程與求解方法.
采用求解直角坐標系下的三維N-S 方程,形式如下

式中,Q=(ρ,ρu,ρv,ρw,ρEt,ρYi)T,E,F,G表示無黏通量,Fv,Gv,Ev表示黏性通量,S為源項,u,v,w為x,y,z方向速度,ρYi表示氣體的密度和組份的質量分數,氣體的總內能Et=e+(u2+v2+w2)/2,其中e表示熱力學內能.
求解三維化學反應控制方程采用隱式有限體積法離散,無黏通量為AUSM PW+格式,黏性通量的計算方法采用Gauss 定理構造方法.湍流模型采用考慮可壓縮性修正的wilcox2006 兩方程湍流模型[27].
流體的可壓縮性對湍流的演變有很大的影響,尤其是脈動壓力和脈動速度散度之間的相關(或稱壓力-散度項)、可壓縮性引起的湍流動能的耗散率(或稱可壓縮性耗散率)的影響最為重要,但在常規兩方程模型中未考慮這兩項的影響,本文對計算軟件進行了基于wilcox2006 模型[27]的可壓縮性修正.
考慮可壓縮修正的wilcox2006 兩方程湍流模型為

考慮膨脹耗散和壓力膨脹的Sarkar可壓縮修正時,對應的參數

算例是Waltrup 經典圓截面隔離段試驗[28]中的一個模型,其直徑69.85 mm,長度578 mm.計算采用的網格如圖1 中所示,網格總量288 萬,壁面網格距離1.0 μm.

圖1 隔離段網格Fig.1 The grid of the isolator
計算狀態為:來流馬赫數2.6,總溫289 K,進口總壓為3.06 atm (1 atm=101.3 kPa),進口靜壓為0.15 atm,出口靜壓為1 atm.分別采用原始wilcox2006 湍流模型及經過可壓縮修正的湍流模型對隔離段進行計算,將計算所得壁壓與試驗壓力對比如圖2 所示.

圖2 采用經過可壓縮修正的湍流模型與原始模型計算所得壓力與試驗對比Fig.2 Comparison of simulations to experimental pressures with compressible modified turbulence model and initial model
可見湍流模型經過Sarkar可壓縮修正后,壓力分布與試驗結果更加吻合,經過可壓縮修正的湍流模型更加準確地預測了的激波串位置.
基于AHL3D 軟件平臺,選用進行可壓縮修正后的湍流模型對日本M12-02 高馬赫數超燃沖壓發動機開展冷態與燃燒狀態的三維數值模擬,并對不同狀態下發動機內流場結構、燃燒釋熱、推阻特性等進行分析.
本文針對圖3 所示日本的M12-02 發動機開展研究.發動機總長2.9 m,由二維側壓式進氣道、平行段、等直燃燒室及擴張噴管組成.進氣道入口截面尺寸為0.25 m × 0.2 m.進氣道與擴張噴管的頂部與底部均為平板,側壁與來流夾角分別為7.4°及8.9°.等直燃燒室長1.7 m,寬0.070 7 m.

圖3 M12-02 發動機構型[23]Fig.3 Configuration of M12-02 scramjet model[23]
燃料為氫氣,為增強混合采用流向渦支板進行燃料噴注,噴注系統詳細結構如圖4 所示[23],噴注裝置位于進氣道之后平行段的側壁上.模擬氫氣燃燒時采用13 組分33 方程簡化模型[29].壁面采用無滑移壁面條件,考慮為等溫壁,壁溫為300 K.

圖4 M12-02 噴射模塊結構[23]Fig.4 Injector structure of M12-02 scramjet[23]
由于所研究的高馬赫數發動機構型沿展向(z方向)對稱,為減少計算量,采用半模計算,總網格量約3000 萬,如圖5 所示.為能準確模擬邊界層流動,壁面法向第一層網格高度為1 μm,以確保y+≤ 1.來流邊界固定為來流參數,出口采用外推.計算過程收斂判斷準則為:流場結構不再明顯改變,密度最大殘差下降3 個量級如圖6 所示,入口流量+噴油量與出口流量差別小于2%,繼續計算2 萬步后流量變化量小于0.2%.

圖5 M12-02 模型計算網格Fig.5 The mesh topology of M12-02 model

圖6 密度殘差收斂曲線Fig.6 Convergence curve of density residual
針對上文介紹的發動機模型開展三維定常數值計算,計算狀態與M12-02 發動機地面試驗狀態保持一致,進氣道入口馬赫數為6.72,速度為3480 m/s,靜壓4.0 kPa,靜溫677 K,氮氣與氧氣質量分數分別為0.203 1 與0.796 9.分別采用未經過可壓縮修正的wilcox2006 模型與經過可壓縮修正的wilcox2006 模型對發動機的冷態與燃燒狀態(當量比(equivalence ratio,ER) 0.5和1.0)進行三維數值模擬,計算所得壁面壓力與文獻中的試驗數據進行對比,如圖7 所示.
從圖7(a)中冷態下的壓力分布可以看出,采用經過可壓縮修正的湍流模型后,軟件對該類問題的模擬能力明顯增強,計算所得壁面壓力波動位置及幅度均與試驗值吻合良好.圖7(b)中則給出了當量比1 條件下氫氣燃燒狀態下的壁壓對比,可見增加了可壓縮修正后,燃燒室后部的壓力與試驗值吻合得相對更好.

圖7 采用經過可壓縮修正的湍流模型與原始模型計算所得壓力與試驗對比Fig.7 Comparison of simulations to experimental pressures with compressible modified turbulence model and initial model
圖8 中對稱面波系分布更直觀地展示了原始模型與修正后模型計算所得激波及其反射波系的差別.可見經過可壓縮修正后計算所得激波角更大,且差異在經過多次激波反射后被放大,燃燒室后部的反射激波位置與原始模型結果相比已存在不可忽略的明顯區別,修正后的模型計算所得激波位置與試驗吻合更好.同時對比兩種模型所得激波后高壓區范圍,可見頭兩道激波后高壓范圍相差不大,但對于后續反射激波,修正后的湍流模型計算得到的高壓區域更大,反映了反射激波強度更大,激波帶來的壓升與試驗值吻合也更好.

圖8 采用經過可壓縮修正的湍流模型與原始模型計算所得壓力云圖對比Fig.8 Comparison of numerical pressure with compressible modified turbulence model to initial model
總的來說,計算結果驗證了經過可壓縮修正的程序對該類高馬赫數問題具有較好的模擬能力.
對于固定來流條件的高馬赫數超燃沖壓發動機,未注油的冷態情形下發動機內激波、膨脹波波系主要由流道結構決定.如圖9 中對稱面馬赫數及波系云圖可見,冷態時進氣道入口處形成上下對稱的激波,進氣道與平行段轉折處形成膨脹波,入口激波與注油位處收縮轉折帶來的激波一同在上下壁面之間反射,等直燃燒室內形成貫穿流道的反射波系,直至噴管擴張處形成膨脹波.

圖9 不同狀態下對稱面馬赫數及波系云圖Fig.9 Centerline planes of Mach number and shock systems for different cases
圖10 中一維平均馬赫數、靜溫、靜壓沿程分布反映了冷態流場參數的整體趨勢,可見平均馬赫數從入口處逐漸下降,直到降至燃燒室出口處達到最小值(3.3),在擴張噴管內不斷增加至出口處的4.18,而平均靜溫趨勢相反,從入口開始逐漸升高至1789 K,從噴管開始下降至出口處的1284 K,反射波系的存在使得平均靜壓有一定的波動,壓力峰值在注油位下游的激波交匯處,靜壓沿燃燒室流向小幅增加.

圖10 不同狀態下一維質量平均馬赫數/靜溫/靜壓Fig.101 D mass averaged Mach number and static pressure and temperature for different cases
區別于雙模態超燃沖壓發動機燃燒的強弱可以改變原有波系,形成激波串等新的流動結構[18],對于工作在Ma12 的高馬赫數發動機而言,燃料注入并點火后,盡管化學反應帶來的熱釋放會引起波系及流動參數變化,但并未改變激波與膨脹波系貫穿流道的基本結構.當量比0.5 時,激波角相較冷態時明顯增大,反射激波數量增多,燃燒室及尾噴管段馬赫數明顯下降.等直燃燒室出口處,平均馬赫數達到最小值2.5,平均靜壓與靜溫增至最大值,分別為68 kPa及2630 K.
當量比1.0 時,激波角進一步增大,激波在上下壁面間反射次數增多,燃燒室內馬赫數進一步下降,溫度整體提升,平均馬赫數最小為2.3,平均靜溫最高為2730 K,均位于燃燒室出口(x=2.435 m)處,但平均靜壓峰值前移至x=2 m 處,為82.15 kPa.對比平均靜溫曲線,可知在燃燒室前半段當量比0.5 的平均靜溫略高于當量比1.0,燃燒室后半段(x> 1.76)當量比1.0 才展現出更高的溫升.
為便于對不同當量比下高馬赫數發動機燃燒特性有更綜合直觀的認識,圖11 中給出了當量比0.5 與1.0 條件下4 種組分(OH,H2O,O 及H)的一維質量加權平均質量分數的分布曲線.

圖11 不同狀態下OH/H2O/O/H 一維質量平均質量分數Fig.111 D mass averaged mass fraction of OH/H2O/O/H for different cases
OH 質量分數的分布規律與上文中一維平均靜溫一致,當量比0.5 時燃燒室前半段OH 更多,燃燒室后半段OH 變化較小,直到噴管處OH 不斷下降,而當量比1.0 時OH 從注油位后不斷攀升直至燃燒室出口/尾噴管入口處.如果說OH 的分布反映了點火與燃燒,那完全產物H2O 則反映了H2與O2反應是否充分完成,燃料熱值能否完全釋放出來.H2O 的分布并非單調遞增,當量比0.5 時在燃燒室中后段(x=1.88 m) H2O 的質量分數達到峰值后在燃燒室尾部略微降低,直到尾噴管中又小幅度增大.當量比1.0 時這種變化更明顯,H2O 質量分數一直增長到燃燒室后部(x=2 m)處,開始明顯下降,直至尾噴管中回升.分析H2O 質量分數的下降主要與燃燒室后部的靜溫水平有關,平均靜溫高達2730 K,局部靜溫更高,高溫使得H2O 更易離解.
當量比的增加更明顯地體現在H 原子質量分數的差異之上,高當量比時H 原子整體更多.而分析O 原子沿程分布規律,可見與當量比1.0 時O 原子在整個燃燒室內逐漸增多相比,當量比0.5 時O 原子生成主要發生在燃燒室前半段,生成速率更快.兩種狀態下,O 原子復合均主要發生在尾噴管內,O 質量分數明顯減小且幅度相當.
圖12 中進一步給出了不同當量比下三維燃燒流場中靜溫、OH 以及H2O 的分布,當量比0.5 時,燃燒室前段溫升更明顯且均勻,OH 也更多,H2O 的生成主要發生在前段,可見貧油狀態下,燃燒距離相對較短,大部分反應在燃燒室前段完成.而當量比1.0 時,燃燒室前段溫升與貧油時相比相對略低,溫升與OH 生成反映了燃燒反應在燃燒室后段一直持續,燃燒距離較長.

圖12 不同狀態下流場中靜溫、OH 及H2O 分布Fig.12 Distribution of static temperature,OH and H2O in flow field for different cases
為進一步增加對燃燒及火焰結構的理解,引入火焰指數GFO[30]定義如下

式中,YO為氧氣質量分數,YF為燃料質量分數,GFO可以區分預混火焰及擴散火焰(非預混),GFO為正時代表此處為預混火焰,GFO為負時則為擴散火焰(非預混).
從圖13 中當量比0.5 及1.0 狀態下流場中火焰指數分布可知,兩種情形下流場中絕大部分的燃燒均屬于非預混燃燒,火焰為擴散火焰.僅在注油裝置hypermixer 沿流向不遠處有少部分預混燃燒區域,且預混火焰的位置恰好在注油孔位置的下游.與當量比0.5 時相比,當量比1.0 時更充分的燃料使得預混燃燒區域相對更多.對燃燒形式的認識可以為下一步高馬赫數發動機數值模擬中模型的選取與進一步優化提供參考.

圖13 不同狀態下流場中火焰指數分布Fig.13 Flame index distribution in flow field for different cases
燃燒反應永遠伴隨著能量的轉化,對于絕熱系統而言,化學反應前后系統總焓值恒定不變,總焓由生成焓與總顯焓(靜顯焓與動能之和)組成.燃燒是將燃料的生成焓轉化為總顯焓的過程.燃燒系統的很大部分能量在化學反應之前都以燃料生成焓的形式存在,化學反應使分子重組.燃燒改變了系統的組分,使反應物向生成物轉化.產物的生成焓小于反應物的生成焓,所減少的生成焓提升了系統的總顯焓.因此對于絕熱系統而言,增加的總顯焓等于減小的生成焓(燃燒釋熱),總顯焓增量此時就是燃燒釋熱量.能量的改變和轉移表現為燃燒釋熱,引起了溫度和壓力的增加,改變了發動機的流動速度.
針對本文研究的等溫壁燃燒系統,化學反應前后系統總焓不再恒定不變,總焓減小量等于壁面傳走熱量,此時總顯焓變化量等于減小的生成焓(即燃燒釋熱)減去壁面傳走熱量.冷態時燃燒釋熱為零,總顯焓減小量等于壁面傳走的熱量,本文稱為“壁面熱損失”.熱態時總顯焓變化量等于燃燒釋熱減去壁面傳走的熱量,本文稱為流道內可利用的“有效釋熱”.其中總顯焓通過以下公式計算

上式中積分沿流向的每個截面進行,ci指的是第i個組分的濃度,Hi(T)與Hi(298.0) 分別指的是溫度為T與298.0 K 時第i個組分的靜焓,|V|2/2 是氣流的動能.V為計算單元表面的流動速度,ds為表面微元的面積向量.
將釋熱量Hrx對x進行求導,可獲得沿流向的釋熱率HRx,如下式所示

計算得到不同狀態下發動機有效釋熱量與釋熱率沿流向分布如圖14 所示.可見冷態時發動機壁面熱損失沿流向基本呈線性分布,壁面傳熱帶來均勻的熱損失.

圖14 不同狀態下燃燒有效釋熱量及釋熱變化率分布Fig.14 Distribution of effective heat release and heat rate for different cases
當量比0.5 與1.0 時,注油位后,燃燒釋熱開始增加,經過一小段距離的累計后可以抵消壁面熱損失,隨后有效釋熱量沿流向逐漸增大至燃燒室后部又開始減小,可見燃燒室后部氫燃燒效果較差.該區域靜溫明顯超過Kutschenreuter[31]提出的2500 K 界限,即靜溫超過2500 K 時,H2與O2反應可獲得的凈燃燒釋熱開始迅速減少,性能顯著下降.與尾噴管中O 原子復合,完全產物H2O 質量分數增加相符,尾噴管有效釋熱略有提升.當量比0.5 與1.0 狀態下,沿流向釋熱率均出現波動的情形,釋熱率較高的波峰處與流場中激波交匯位置較為吻合,可見激波交匯帶來的溫升與壓升有利于燃燒釋熱.
將燃燒釋熱(有效釋熱與壁面傳走熱量之和)用氫燃料熱值無量綱化得到發動機燃燒效率(熱值轉化率)曲線如圖15 所示.可見當量比0.5 時燃料少,更利于摻混,燃燒效率高,發動機出口處約為90%,當量比1.0 時發動機出口燃燒效率約為60%.且根據燃燒釋熱計算的燃燒效率在燃燒室后部也出現減小的情況,這與H2O 離解吸熱相符.

圖15 不同狀態下燃燒效率分布Fig.15 Distribution of combustion efficiency for different cases
在分析發動機性能時,壁面傳熱導致的能量損失是不可忽略的一部分,影響著燃燒熱的利用率,圖16以當量比1.0 狀態為例給出了無量綱壁面熱流分布(以進氣道入口段壁面熱流無量綱),反映了壁面熱流的相對變化量.圖中同時給出了對稱面波系分布,可見激波與壁面交匯區域壁面熱流明顯升高,且隨著反射激波沿流向減弱,所帶來的壁面熱流升高逐漸減少,激波在流道中央交匯的位置壁面熱流則相對略小.高馬赫數發動機內激波與反射波系常貫穿于流道中,發動機總體設計中必須考慮激波與壁面交匯帶來的高熱流區域.

圖16 當量比1.0 狀態下無量綱壁面熱流分布Fig.16 Nondimensional wall heat flux distribution for equivalence ratio 1.0 case
最后,分析計算所得M12-02 發動機各部件的推阻力情形,采用動壓與面積的乘積對推力進行無量綱,得到推力系數

其中 ρ∞為來流密度,u∞為來流速度,A為進氣道入口的橫截面積.以同樣的方式對摩阻進行無量綱,給出各分部件摩阻系數CFD.計算所得整機及各分部件推力系數與摩阻系數如表1 所示.

表1 M12-02 發動機不同部件推力系數及摩阻系數Table 1 Thrust and friction coefficients of different components for M12-02 scramjet
可見不同狀態下燃燒室均為阻力部件,燃燒與冷態相比燃燒室阻力變小,但當量比從0.5 提升至1.0 燃燒室阻力變化較小,總推力系數的差異主要由尾噴管貢獻.不同狀態下進氣道摩阻基本不變,尾噴管摩阻差異較小,燃燒會導致燃燒室摩阻及整機總摩阻減小.
針對飛行Ma12 條件下M12-02 超燃沖壓發動機開展計算方法的改進以及多狀態下精細三維數值模擬,分析發動機內激波與膨脹波波系、參數分布以及燃燒性能特征.結果表明:
(1) 基于AHL3D 軟件,對計算方法完成了基于wilcox2006 模型的可壓縮性修正.修正后的方法計算所得激波位置及強度與M12-02 發動機試驗值吻合良好,在激波串模擬、高馬赫數發動機模擬上均展現了更優的能力.
(2) 區別于雙模態發動機燃燒的強弱可以改變原有波系并形成激波串等新的流動結構,該高馬赫數發動機內形成貫穿流道的激波與反射波系,燃燒熱釋放并未改變波系貫穿流道的基本結構,且隨著當量比增加,激波角增大,反射激波數量增多.
(3) 當量比0.5 時,燃燒室前半段溫升、OH 及O 原子質量分數比當量比1.0 時更高,燃燒反應主要發生在燃燒室前部;當量比1.0 時,OH 與靜溫在整個燃燒室內沿流向不斷升高,反應距離更長.O 原子復合主要發生在擴張噴管中.
(4) 流場中絕大部分的燃燒均屬于非預混燃燒,火焰為擴散火焰,可為下一步高馬赫數發動機數值模擬中模型的選取與進一步優化提供參考.
(5) 燃燒室后段平均靜溫超過2500 K,完全產物H2O 減少,H2與O2燃燒效果變差.發動機可利用的有效釋熱在燃燒室前段增加,在燃燒室后段減小.當量比0.5 與1.0 下燃燒室阻力差值較小,總推力系數差異主要由尾噴管貢獻.不同狀態下進氣道摩阻基本不變,尾噴管摩阻差異較小,燃燒會導致燃燒室摩阻及整機總摩阻減小.
(6) 激波與壁面交匯區域壁面熱流明顯升高,且隨著反射激波沿流向減弱,所帶來的壁面熱流升高逐漸減少,發動機總體設計中必須考慮激波帶來高熱流區域的影響.