空氣渦輪火箭發(fā)動機風(fēng)車狀態(tài)數(shù)值仿真研究

張留歡，逯婉若

（西安航天動力研究所，陜西西安710100）

0 引言

空氣渦輪火箭發(fā)動機是航空渦輪噴氣發(fā)動機與火箭發(fā)動機的有機組合。與航空發(fā)動機空氣路與主燃燒室相通不同，ATR發(fā)動機利用流路及供應(yīng)獨立的火箭發(fā)動機（或燃氣發(fā)生器）工作產(chǎn)生的高溫、高壓燃氣驅(qū)動渦輪，進而帶動壓氣機對來流空氣進行增壓，之后空氣進入燃燒室 (類似渦噴加力燃燒室)與渦輪出口燃氣摻混、燃燒，并通過尾噴管高速排出產(chǎn)生推力，見圖1。ATR發(fā)動機燃氣發(fā)生器氣體流路與空氣流路并聯(lián)，在進入燃燒室前互不影響。渦輪前燃氣參數(shù) (總溫等)不受發(fā)動機來流空氣參數(shù)影響。這在一定程度上拓寬了ATR發(fā)動機的穩(wěn)定工作包線 (馬赫數(shù)范圍可達 0～6)[1-2]。

圖1 ATR發(fā)動機（燃氣發(fā)生器循環(huán)）Fig.1 Gas generator cycle of ATR engine

以ATR發(fā)動機為動力的飛行器在飛行過程中，可能遇到無動力下壓滑翔狀態(tài)。此時，發(fā)動機關(guān)機，壓氣機轉(zhuǎn)子在來流作用下做自轉(zhuǎn)運動，發(fā)動機處于風(fēng)車狀態(tài)，發(fā)動機的內(nèi)流道阻力（以下簡稱內(nèi)阻）等參數(shù)在很大程度上影響著整個飛行器下壓飛行的時間、距離。因此，獲得發(fā)動機風(fēng)車狀態(tài)參數(shù)（包含風(fēng)車轉(zhuǎn)速、內(nèi)阻等參數(shù)）對飛行彈道規(guī)劃有重要意義。

在國內(nèi)外研究方面，尤·阿·李特維諾夫指出在雷諾數(shù)自行?；母黠w行狀態(tài)范圍內(nèi)，如果轉(zhuǎn)子沒有消耗功率，任何結(jié)構(gòu)形式的發(fā)動機的所有換算參數(shù)實際上都是單值的，取決于設(shè)計馬赫數(shù)，并由此給出了計算風(fēng)車轉(zhuǎn)速和內(nèi)阻的計算公式[3]；王占學(xué)等人在尤·阿·李特維諾夫方法的基礎(chǔ)上，發(fā)展了計算燃氣渦輪發(fā)動機風(fēng)車特性的理論方法[4]；劉志友等人利用高空臺開展了某型渦扇發(fā)動機風(fēng)車試驗，并依據(jù)試驗數(shù)據(jù)推導(dǎo)了針對特定發(fā)動機構(gòu)型的風(fēng)車內(nèi)阻計算經(jīng)驗公式[5]；陳建民等針對航空發(fā)動機風(fēng)車狀態(tài)起動特性開展了部分試驗和理論研究[6-9]。美國飛行推進研究試驗中 心 （Flight Propulsion Research Laboratory）在20世紀中葉開展了大量的燃氣渦輪發(fā)動機風(fēng)車特性試驗研究，獲得了豐富的試驗數(shù)據(jù)，不過試驗成本較高[10-14]。

本文采用CFD（Computational Fluid Dynamics）方法，對特定構(gòu)型的ATR發(fā)動機開展了冷態(tài)條件下發(fā)動機內(nèi)流道三維流場數(shù)值仿真，獲得了一定來流條件下ATR發(fā)動機的風(fēng)車轉(zhuǎn)速及對應(yīng)的流道內(nèi)阻，可為后續(xù)發(fā)動機試驗及飛行器總體設(shè)計工作的開展提供借鑒。

1 結(jié)構(gòu)及氣動參數(shù)

本文研究的某構(gòu)型ATR發(fā)動機設(shè)計轉(zhuǎn)速30 000 r/min，最大外徑300 mm，長度約1.5 m。發(fā)動機沿氣流流向依次包含斜流壓氣機（空氣路）、燃氣發(fā)生器和沖擊式渦輪（燃氣路）、圓形燃燒室及收擴尾噴管等組件。

ATR發(fā)動機處于風(fēng)車狀態(tài)時，進入發(fā)動機的空氣經(jīng)過壓氣機后，直接進入圓形等徑燃燒室，最后通過收擴尾噴管排出。此時，燃氣發(fā)生器未工作，渦輪流道內(nèi)不存在燃氣。

本文選取了飛行過程中的某一工況點（飛行高度20 km，速度3 Ma），開展了ATR發(fā)動機風(fēng)車狀態(tài)特性計算。飛行工況氣動參數(shù)具體見表1。

表1 飛行工況參數(shù)Tab.1 Parameters of flight condition

表1中的入口為壓氣機入口截面。其中，考慮進氣道總壓損失，給定壓氣機入口總壓0.15 MPa。

2 數(shù)值仿真

ATR發(fā)動機空氣與燃氣流路并聯(lián)，且在關(guān)機時，燃氣路不存在燃氣。在風(fēng)車狀態(tài)時，壓氣機轉(zhuǎn)子、渦輪轉(zhuǎn)子及其軸系在轉(zhuǎn)動周向上所受合力為零，且合力對轉(zhuǎn)軸的扭矩為0。由于軸系摩擦力、渦輪轉(zhuǎn)子流路反竄空氣對轉(zhuǎn)子系統(tǒng)的影響較小，當(dāng)不考慮以上因素時，可認為來流空氣對壓氣機轉(zhuǎn)子在周向上的作用合力為0，壓氣機轉(zhuǎn)子葉片所受的氣體作用力對轉(zhuǎn)軸的扭矩為0。

在開展冷態(tài)條件下ATR發(fā)動機內(nèi)流道流場數(shù)值仿真過程中，本文首先給定發(fā)動機一系列轉(zhuǎn)速（3 000～9 000 r/min），計算收斂后提取壓氣機轉(zhuǎn)子扭矩，將其作為判定參數(shù)，當(dāng)扭矩為0即為風(fēng)車狀態(tài)，對應(yīng)的轉(zhuǎn)速即為風(fēng)車轉(zhuǎn)速。在給定的風(fēng)車轉(zhuǎn)速下，獲得壁面沿程壓力和粘性摩擦力，并由公式（1）計算得到內(nèi)阻值：

2.1 計算模型

根據(jù)ATR發(fā)動機構(gòu)型，采用ICEM軟件建立了三維流場計算模型（如圖2所示）。

圖2 ATR發(fā)動機三維計算網(wǎng)格Fig.2 3D mesh on ATR engine

模型采用混合網(wǎng)格劃分，其中燃燒室前的空氣流域采用非結(jié)構(gòu)網(wǎng)格劃分，燃燒室及尾噴管流域采用結(jié)構(gòu)網(wǎng)格劃分。發(fā)動機入口設(shè)為壓力入口，尾噴管出口設(shè)置計算遠場，遠場邊界設(shè)為壓力遠場，遠場沿流向出口設(shè)為壓力出口。對壓氣機葉柵等流動復(fù)雜區(qū)域及壁面附近網(wǎng)格進行加密，整個模型網(wǎng)格量約250萬。

2.2 計算方法

數(shù)值仿真采用Fluent軟件進行。計算采用有限體積法，求解方法為基于密度的隱式求解法，采用Sutherland公式計算分子粘性系數(shù)，壁面取絕熱無滑移和固體邊界條件，壁面附近采用非平衡壁面函數(shù)處理，湍流模型采用RNG（Re-Normalization Group）k-ε湍流模型。由于壓氣機轉(zhuǎn)子流域為三維旋轉(zhuǎn)流動，因此，計算中使用多參考坐標(biāo)系模型，給定壓氣機轉(zhuǎn)子流域流體轉(zhuǎn)動角速度以及對應(yīng)的靜子壁面相對角速度。計算過程中，監(jiān)測質(zhì)量流量、動量、能量等殘差數(shù)量級降至10-3以下，認為迭代達到收斂。

2.3 結(jié)果及分析

數(shù)值仿真獲得了冷態(tài)條件下發(fā)動機內(nèi)流道流場及氣動參數(shù)結(jié)果。圖3給出了轉(zhuǎn)速n為6 900 r/min時，ATR發(fā)動機不同截面流場的馬赫數(shù)分布情況。

圖3 不同截面流場馬赫數(shù)分布云圖（n=6 900 r/min）Fig.3 Mach number distribution of flowfield at different sections（n=6 900 r/min）

在ATR發(fā)動機中，壓氣機擴壓器葉柵通道入口幾何面積小于噴管喉部面積，為整個空氣流道的最小截面。在飛行高度20 km、飛行速度Ma 3、入口總壓0.15 MPa的條件下，發(fā)動機入口馬赫數(shù)約0.16。經(jīng)過壓氣機時，流通面積相對發(fā)動機入口有所減小，氣流速度增加。在擴壓器一級靜子葉柵通道入口，流通面積進一步縮小，在一級靜子葉片葉背處局部氣流首先達到聲速。之后，葉柵通道面積增加，且沿流向超聲速區(qū)域逐漸增大，在二級靜子葉柵通道出口截面氣流馬赫數(shù)局部達到Ma 2。擴壓器后，流道突然擴張，超聲速氣流與葉片尾跡區(qū)低速流相互作用，氣流動能逐漸耗散，降至亞聲速。在燃燒室內(nèi)，流道面積再次增加，入口速度約Ma 0.8的空氣氣流與渦輪盤后的低速氣流（約Ma 0.1）沿流向逐漸摻混，在燃燒室出口處氣流均勻，馬赫數(shù)約為0.5。氣流到達尾噴管時，由于流道面積收縮，氣流速度增加，且尾噴管落壓比達到12，噴管喉部達到聲速。

同時，沿氣流流向，由于壁面粘性摩擦、氣流分離、氣流剪切摻混等損失的影響，氣流總壓逐漸降低。發(fā)動機進出口總壓損失約61%。其中，壓氣機流道進出口總壓損失達到了32.6%，發(fā)動機不同截面流場的總壓分布云圖見圖4。

圖4 不同截面流場總壓分布云圖（n=6 900 r/min）Fig.4 Total pressure distribution of flowfield at different sections（n=6 900 r/min）

圖5給出了不同給定轉(zhuǎn)速下ATR發(fā)動機入口氣動參數(shù)的變化曲線。隨著給定壓氣機轉(zhuǎn)子轉(zhuǎn)速的不斷升高，發(fā)動機通流能力逐漸增強，發(fā)動機入口氣流靜壓p1逐漸降低，速度V1逐漸增加。

圖6給出了不同給定轉(zhuǎn)速下ATR發(fā)動機轉(zhuǎn)子扭矩變化曲線。圖6中，扭矩正值表示沿流向看去，轉(zhuǎn)子周向所受氣動力為逆時針方向；反之，為順時針方向?？梢钥吹?，隨著發(fā)動機轉(zhuǎn)速的逐漸升高，來流空氣對壓氣機轉(zhuǎn)子作用扭矩由正值逐漸變?yōu)樨撝?，即轉(zhuǎn)子周向受力由逆時針變?yōu)轫槙r針方向。轉(zhuǎn)速約為6 900 r/min時，轉(zhuǎn)子扭矩接近0，表明此時發(fā)動機處于風(fēng)車狀態(tài)。根據(jù)數(shù)值仿真結(jié)果，計算發(fā)動機內(nèi)流道壁面受力在軸向上的分量，得到發(fā)動機軸向內(nèi)阻約為2 170 N。

圖5 不同轉(zhuǎn)速下入口參數(shù)變化曲線Fig.5 Parameters at entrance at different speeds

圖6 不同轉(zhuǎn)速下轉(zhuǎn)子扭矩變化曲線Fig.6 Curve of rotor torque at different speeds

3 結(jié)論

1）在飛行高度20 km，飛行速度Ma 3及入口總壓0.15 MPa的條件下，ATR發(fā)動機風(fēng)車轉(zhuǎn)速約為6 900 r/min，內(nèi)阻約2 170 N。

2）ATR發(fā)動機處于風(fēng)車狀態(tài)時，沿氣流流向，氣流總壓逐漸降低。在飛行高度20 km、飛行速度Ma 3、入口總壓0.15 MPa的條件下，處于風(fēng)車狀態(tài)的ATR發(fā)動機進出口總壓損失約61%。其中，壓氣機流道進出口總壓損失達到了32.6%。

3）冷態(tài)條件下，隨著給定壓氣機轉(zhuǎn)子轉(zhuǎn)速的不斷升高，發(fā)動機通流能力逐漸增強，發(fā)動機入口氣流靜壓逐漸降低，速度逐漸增加。

實際飛行過程中，由于軸系摩擦、氣流分離等因素的影響，風(fēng)車轉(zhuǎn)速可能與數(shù)值仿真結(jié)果有所出入。后續(xù)工作需優(yōu)化計算模型以盡可能接近實際工況，提高計算結(jié)果的可信度。

[1]TAKESHI K.Combined-cycle engines,encyclopedia of aerospace engineering[M].[S.l.]：John Wiley&Sons,2010.

[2]NOBUHIRO T.Development study on airtur-boramjet,developments in high-speed-vehicle propulsion systems[M].[S.l.]：Progress in Astronautics and Aeron-autics,1996.

[3]李特維諾夫尤·阿.航空渦輪噴氣發(fā)動機的特性和使用性能[M].國防工業(yè)出版社,1986.

[4]王占學(xué),劉增文.某型燃氣渦輪發(fā)動機風(fēng)車狀態(tài)內(nèi)阻力的計算[J].燃氣渦輪試驗與研究,2006,19(3)：8-10.

[5]劉志友,侯敏杰,文剛.航空發(fā)動機風(fēng)車阻力的試驗研究[J].航空動力學(xué)報,2006,21(2)：391-395.

[6]田琳,陶冶,申曉霞.民用渦扇發(fā)動機空中啟動風(fēng)車轉(zhuǎn)速研究[J].工程與試驗,2013,53(1)：6-9.

[7]楊欣毅,沈偉,劉海峰,董可海.一種彈用渦噴發(fā)動機風(fēng)車起動數(shù)值仿真方法[J].航空動力學(xué)報,2015,25(8)：1776-1782.

[8]張東方,李應(yīng)紅,尉詢楷.計算渦扇發(fā)動機風(fēng)車起動特性的辨識模型[J].航空動力學(xué)報,2007,22(8)：1320-1324.

[9]萬照云.微型渦輪發(fā)動機風(fēng)車起動特性研究[D].南京航空航天大學(xué),2011.

[10]CONRAD E W,DURHAM J D.Preliminary results of an altitude-wind-tunnel investigation of a TG-10CA gas turbine-propeller engineⅡ-Windmilling characteristics,NACA RM No.E7G25[R].USA：NACA,1947

[11]CONRAD E W.,DURHAM J D.Preliminary results of an altitude-wind-tunnel investigation of an axial-flow gas turbine-propeller engine II-Windmilling characteristics,NACA RM No.E8F10a[R].USA：NACA,1948.

[12]FLEMING W A,DIETZ R O,Jr.Altitude-wind-tunnel investigation ofWestinghouse 19B-2,19B-8,and 19XB-1 jet-propulsion engines III performance and windmilling drag characteristics,NACA RM No.E8J28b[R].USA：NACA,1948

[13]SOBOLEWAKI A E,FARLEY J M.Steady-state engine windmilling and engine speed decay characteristics of an axial-flow turbojetengine,NACA RM No.E51I06[R].USA：NACA,1951.

[14]VINCENT K R,HUNTLEY S C.,WILSTED H D.Comparison of locked-rotor and windmilling drag characteristics of an axial-flow-compressor type turbojet engine,NACARMNo.E51K15[R].USA：NACA,1952.

[15]南向誼,王拴虎,李平.空氣渦輪火箭發(fā)動機研究的進展及展望[J].火箭推進,2008,34(6)：31-35.

NAN Xiangyi,WANG Shuanhu,LI Ping.Investigation on status and prospect of air turbine rocket[J],Journal of Rocket Propulsion,2008,34(6)：31-35.

[16]李文龍,郭海波,南向誼.空氣渦輪火箭發(fā)動機熱力循環(huán)特性分析[J].火箭推進,2015,41(4)：48-54.

LI Wen-long,GUO Hai-bo,NAN Xiang-yi.Analysis on thermodynamic cycle characteristics of air-turbo-rocket engine[J].Journal of Rocket Propulsion,2015,41(4)：48-54.