張 皓，黃 萌，鄧 恒，田小濤，顏 密，賈勝錫

(1 西安現(xiàn)代控制技術(shù)研究所，陜西 西安 710065；2 現(xiàn)代控制技術(shù)重點實驗室，陜西 西安 710065)

0 引言

作為高超聲速巡航飛行器潛在動力系統(tǒng)之一,固體燃料超燃沖壓發(fā)動機(以下簡稱：固體超燃)具有比沖高、結(jié)構(gòu)簡單、成本低、維護方便等優(yōu)點。與固體火箭發(fā)動機不同,由于固體超燃具有進氣道部件,發(fā)動機內(nèi)部工作效率與外部飛行環(huán)境強相關(guān),導(dǎo)致采用固體超燃的飛行器在總體設(shè)計時無法將發(fā)動機作為獨立部件模塊化處理[1]。采用固體超燃的高超聲速巡航飛行器總體設(shè)計高度依賴發(fā)動機-飛行器機體一體化設(shè)計技術(shù)。內(nèi)彈道計算與設(shè)計作為固體超燃高超聲速飛行器一體化設(shè)計的關(guān)鍵技術(shù),其準確計算與預(yù)估至關(guān)重要。與固體火箭發(fā)動機不同,固體超燃內(nèi)彈道與外界環(huán)境相關(guān)性強,且內(nèi)彈道的計算涉及到高超聲速流場中激波的產(chǎn)生、相交、反射,火焰穩(wěn)定區(qū)低速回流燃燒放熱,固體顆粒兩相流燃燒等物理化學現(xiàn)象,一般需要開展二維/三維數(shù)值計算進行預(yù)測。二維/三維數(shù)值計算雖然可以較為準確預(yù)測內(nèi)彈道性能,但計算時間較長。在進行內(nèi)彈道計算與設(shè)計時,設(shè)計變量較多,往往需要開展多工況計算,計算時間成本較大,不利于飛行器一體化總體設(shè)計的快速迭代。一維或準一維內(nèi)彈道計算模型在允許的精度范圍內(nèi),簡化各個物理化學過程,考慮主要因素,能極大減少計算時間,快速對發(fā)動機內(nèi)彈道性能進行評估與預(yù)測。

目前,超燃沖壓發(fā)動機的準一維計算模型主要有兩類：一類是以空間為基礎(chǔ)的常微分方程組穩(wěn)態(tài)求解法[2];另一類是以時間為基礎(chǔ)的偏微分方程組的非穩(wěn)態(tài)求解法[3]。二者都考慮了超聲速流動的主要物理過程。在此基礎(chǔ)上,張鵬等[4]對幾種公開的典型燃燒室流場建立了一維模型,完善了相關(guān)實驗分析。范學軍等[5]利用一維模型開展了普通煤油、加熱煤油和超臨界煤油在超聲速流場中的燃燒效率研究,認為超臨界煤油燃燒效率高于冷態(tài)煤油。陳強[6-7]通過實驗數(shù)據(jù)對一維模型的釋熱模塊進行參數(shù)優(yōu)化,并用該一維模型對燃燒室結(jié)構(gòu)進行了優(yōu)化。王厚慶等[8]通過一維計算模型開展了超燃內(nèi)彈道計算,并以該計算結(jié)果為邊界條件輸入,對碳化硅陶瓷復(fù)合材料的燃燒室冷卻結(jié)構(gòu)進行了優(yōu)化設(shè)計。Srinivasan等[9]對一維計算模型進行了補充,主要考慮了液滴蒸發(fā)、有限速率化學反應(yīng)、傳熱與摩擦等因素對計算結(jié)果的影響。Micka等[10]基于化學熒光測量法測試了燃燒室釋熱分布,并用壁面壓力修正了一維計算模型的相關(guān)參數(shù),提高了一維模型對該類燃燒室構(gòu)型的預(yù)估精度。Tomioka等[11-12]在一維模型中加入了燃油噴注簡化模型,并利用該模型,以燃燒效率、推力、摩擦阻力作為目標參數(shù),研究了燃油噴注位置對超燃性能的影響。Kliche等[13]建立了基于Euler方程的一維非定常計算模型,該模型考慮了有限速率化學反應(yīng)模型、平衡點火模型的影響。

目前為止,對超聲速流場的一維計算模型已開展了較為豐富的研究,模型功能也較為豐富。但由于固體超燃燃燒室內(nèi)流場中普遍存在亞聲速-跨聲速-超聲速流動現(xiàn)象[14],一維計算模型在進行此類問題計算時,控制方程會在聲速臨界點存在數(shù)學奇異性。此外,在進行部分工況計算時,由于方程組參數(shù)數(shù)量級相差過大,會引起嚴重的數(shù)值剛性問題。在進行飛行器總體設(shè)計時,外彈道計算希望將發(fā)動機看做一個“黑盒”部件,確定的輸入對應(yīng)一個確定的輸出,對發(fā)動機內(nèi)部工作細節(jié)并不關(guān)心。因此,文中對一維模型進一步簡化,基于能量守恒方程,建立一個以推力為目標的零維模型,研究燃料質(zhì)量流量與空氣流量變化對貼壁式固體超燃沖壓發(fā)動機推力的影響。

1 推力零維模型

1.1 基本假設(shè)

基于典型燃燒室構(gòu)型(圖1),建立燃料供應(yīng)量對貼壁式固體超燃沖壓發(fā)動機推力的影響模型。為建立燃料供應(yīng)量與推力之間的零維模型,現(xiàn)假設(shè)如下：

圖1 固體超燃沖壓發(fā)動機典型構(gòu)型Fig.1 Typical configuration of sold fuel scramjet

1)不考摩擦的影響。由于燃燒室長度較短,而固體燃料貼壁式超燃沖壓發(fā)動機內(nèi)氣流流動速度為1 000 m/s以上的,如此短的行程導(dǎo)致摩擦的影響較小。

2)認為燃燒室內(nèi)的氣體為理想氣體。

3)不考慮發(fā)動機的點火過程,認為燃燒室內(nèi)已經(jīng)建立起穩(wěn)定工作狀態(tài)。

4)不考慮燃料的熱解過程,認為固體燃料瞬間熱解并向燃燒室內(nèi)注入燃料氣體。

1.2 模型建立

基于以上假設(shè),根據(jù)能量守恒原則,可以列出方程：

(1)

根據(jù)質(zhì)量守恒,可以列出方程：

(2)

式中：ρe為出口氣流密度;ue為出口氣流速度;Ae為燃燒室出口截面積。

根據(jù)理想氣體定律：

Pe=ρeRTe

(3)

式中：Pe為出口靜壓;R為氣體常數(shù);Te為出口靜溫。

根據(jù)總溫與靜溫的等熵關(guān)系式,有如下方程：

(4)

式中：k為比熱比;M為出口馬赫數(shù)。

馬赫數(shù)可以通過式(5)計算：

(5)

聯(lián)立求解式(1)～式(5),可得：

(6)

按照文獻[15]中的方程計算地面直連試驗的發(fā)動機推力為：

(7)

式中Pa為出口反壓,即大氣壓。

聯(lián)立式(5)和式(7),可得出口馬赫數(shù)的計算方程：

(8)

1.3 模型驗證

對比Ben-Yakar的試驗數(shù)據(jù)[15]與上節(jié)推力零維模型的計算結(jié)果以驗證該模型的正確性。根據(jù)Yakar的實驗條件,取計算條件如表1所示。

表1 根據(jù)Yakar實驗確定的計算條件Table 1 Computational condition determined based on Yakar′s experiment

推力零維模型計算所得的理論推力曲線與Yakar的實驗結(jié)果對比如圖2所示。從圖中可以發(fā)現(xiàn),理論推力與Yakar的實驗推力結(jié)果總體吻合的較好。理論推力與實驗結(jié)果相比,理論推力整體偏低,且隨著發(fā)動機工作過程的進行,理論推力與實驗結(jié)果偏差逐漸增大。這可能是由于零維推力模型計算中只考慮了平均燃面退移速率,且使用了平均燃燒效率造成的。

圖2 推力對比Fig.2 Thrust contrast

零維推力模型的馬赫數(shù)和出口速度與Yakar的實驗結(jié)果對比如圖3和圖4所示。從圖中可以看出,出口馬赫數(shù)的理論結(jié)果較實驗值整體偏低,而出口速度的計算結(jié)果與實驗中吻合的較好。文中的零維推力模型主要用于研究空氣流量與燃料供應(yīng)量對推力的影響,因此在建立模型時未考慮燃料的熱解行為,而是將燃料供應(yīng)量作為輸入量用于發(fā)動機推力的計算。零維推力模型未考慮摩擦的影響,導(dǎo)致計算的理論結(jié)果偏低。但從圖2和圖4可以看出,理論結(jié)果與實驗結(jié)果偏差較小,且趨勢一致,因此,建立的零維推力模型是考慮了固體超燃工作過程中的主要因素,具有一定的精度,可以用來初步研究空氣流量與燃料供應(yīng)量對發(fā)動機推力的影響。

圖3 馬赫數(shù)對比Fig.3 Mach contrast

圖4 出口速度對比Fig.4 Outlet velocity contrast

2 結(jié)果與討論

將式(6)代入式(7)中可得：

(9)

式中：De,0為發(fā)動機出口初始直徑;而Ae可表示為：

(10)

當t=0時,

(11)

將式(11)代入式(9)中可得：

(12)

對式(12)求關(guān)于De,0的導(dǎo)數(shù)可得：

(13)

整理式(13)得：

(14)

在式(14)中,顯然有：

(15)

那么對于任意的De,0,式(16)成立：

(16)

由式(16)可知,初始出口直徑越大,推力越小。該結(jié)論在李彪[16]的數(shù)值仿真計算結(jié)果中也有所體現(xiàn)。李彪的計算模型與文中討論的固體燃料超燃沖壓發(fā)動機典型構(gòu)型與圖1相同,因此,擴張半角越大,出口的初始直徑越大。推力與比沖隨著出口的初始直徑下降,這也符合式(16)的結(jié)論。

為研究在貼壁式固體燃料超燃沖壓發(fā)動機燃燒室工作過程中空氣流量與燃料質(zhì)量流量的變化對推力的影響。對式(12)求時間t的導(dǎo)數(shù)：

(17)

其中,

(18)

只要求出式(17)的零點,即式(19)的解,便可以了解推力的變化過程。

(19)

(20)

對式(20)求時間t的導(dǎo)數(shù)：

(21)

假設(shè)入口空氣質(zhì)量流量為：

(22)

則有：

(23)

系數(shù)a,b,c的物理含義分別為：系數(shù)a與燃燒過程中燃面退移速率與燃面面積相關(guān);系數(shù)b為燃燒室初始時刻燃面面積與燃面退移速率的乘積;系數(shù)c為空氣質(zhì)量流量。通過設(shè)計不同的燃燒室燃面構(gòu)型及選用不同燃面退移速率的固體燃料,可以得到不同的a,b值。需要注意的是,系數(shù)a有負值,而系數(shù)b和系數(shù)c必須是正值才有意義。

對式(19)使用4-5階龍格庫塔法迭代求解,所得零點為t*,稱之為零點時間。圖5所示為零點時間t*隨系數(shù)a的變化曲線。從圖中可以看出,當設(shè)定b=0.04和c=0.5時,不同系數(shù)a下的零點時間t*不同,且t*與系數(shù)a呈現(xiàn)正相關(guān)特性。當t*<0時,則在t>0的范圍內(nèi)不存在零點時間,即推力F在t>0的時間段內(nèi)是單調(diào)的;而當t*>0時,在t>0的時間段內(nèi)存在零點;當0t*時,推力開始下降,當t=t*時,推力到達理論峰值。

圖5 零點時間隨系數(shù)a的變化曲線Fig.5 Zero-point time varies with a

設(shè)b=0.04,c=0.5。當系數(shù)a=0時的推力曲線如圖6所示。a=0表示燃料質(zhì)量流量不隨時間變化,此時零點時間t*<0,從圖中可以看出,這種情況下推力會隨時間逐漸降低,且降低幅度明顯,在20 s時,理論推力就已經(jīng)降低了30%左右。到50 s時,理論推力下降程度已經(jīng)超過了初始推力的一半。

圖6 系數(shù)a=0時的推力曲線Fig.6 Thrust curve when a=0

當系數(shù)a=0.002時的推力曲線如圖7所示。

圖7 系數(shù)a=0.002時的推力曲線Fig.7 Thrust curve when a=0.002

由圖5可知,此時零點時間t*>0。則在t>0的時間段內(nèi)存在推力變化率為0的點。從圖7可以明顯看出,推力隨時間先增大,在20 s左右到達推力的峰值,然后推力開始下降。雖然這種情況下的推力有增加的趨勢,但初始推力與峰值推力相差不大。在50 s 的工作時間內(nèi),推力整體變化的幅值不大。

當系數(shù)a=0.003時的推力曲線如圖8所示。在這種情況下,推力在前50 s內(nèi)會一直增加,且推力增加幅度較a=0.002情況有明顯增大。雖然推力在前50 s內(nèi)一直增大,但從圖5中可以看出,這種情況下存在零點時間t*,即存在一個推力的峰值,在到達峰值前,推力會一直上升。

圖8 系數(shù)a=0.002時的推力曲線Fig.8 Thrust curve when a=0.003

從上述關(guān)于系數(shù)a對推力影響的討論中,可以了解到在質(zhì)量流量恒定的情況下,推力會隨著時間一直降低;選取一個合適的系數(shù)a可以讓推力在工作期間的變化幅度小很多,而選取較大的系數(shù)a會使得推力在發(fā)動機工作期間迅速增大,且系數(shù)a越大,推力的峰值到來越晚。

下面討論系數(shù)b對推力的影響。在系數(shù)a=0.002,c=0.5的計算條件下,不同系數(shù)b下的零點時間t*如圖9所示。

圖9 不同系數(shù)b下的零點時間t*Fig.9 Zero-point time varies with b

零點時間t*隨系數(shù)b的變化趨勢與零點時間t*隨系數(shù)a的變化趨勢相反。系數(shù)b的取值較小時t*>0,即推力存在一個峰值,使得推力先增大后減小。當b的取值較大時,t*<0,此時推力隨時間變化沒有峰值,推力隨時間變化呈現(xiàn)單調(diào)性。

當b=0.03時,推力隨時間的變化曲線如圖10所示,此時t*>0,即推力存在一個峰值。從圖中可以得知,推力在初始時刻的大小為760 N,在25 s前推力隨時間緩慢增加,在25 s左右時達到推力峰值,隨后開始緩慢下降。從0～25 s和25～50 s的推力曲線對比可以看出,推力上升的幅度要高于推力下降的幅度。

圖10 系數(shù)b=0.03時的推力曲線。Fig.10 Thrust curve when b=0.03

當b=0.1時的推力隨時間變化曲線如圖11所示,該種情況下t*<0,即在發(fā)動機工作過程中不存在推力峰值,推力隨時間單調(diào)變化。從圖中可知,推力在初始時刻的大小約為1 260 N,隨著時間的退移,推力一直下降,到50 s時,推力已經(jīng)下降到1 120 N左右。通過對比圖10和圖11中推力的變化幅度,可以發(fā)現(xiàn)推力在單調(diào)變化時變化幅度較大。

圖11 系數(shù)b=0.1的推力曲線圖Fig.11 Thrust curve when b=0.1

為考察系數(shù)a,b,c對推力變化幅度的影響,定義推力變化率如下：

(24)

式中：Fmax為工作時間內(nèi)的最大推力;Fmin為工作時間內(nèi)的最小推力;F|t=0為初始時刻的推力。

當系數(shù)b=0.04,c=0.5時推力變化率隨系數(shù)a的變化曲線如圖12所示。推力變化率ΔF為1的情況對應(yīng)于推力在工作期間下降為0的理論極限情況。從圖中可以看出,當系數(shù)a=0.002時,推力變化率ΔF最小;當a<0.002時,推力變化率隨著系數(shù)a的增大而降低;當系數(shù)a>0.002時,推力變化率ΔF會隨著系數(shù)a的增大而增大。

圖12 推力變化率隨系數(shù)a的變化曲線Fig.12 Curve of thrust change rate with a

當系數(shù)a=0.002,c=0.5時推力變化率ΔF隨系數(shù)b的變化曲線如圖13所示。當b=0.04時,推力變化率ΔF最小;當b<0.04時,推力變化率ΔF隨系數(shù)b的增加而降低;而當b>0.04后,推力變化率ΔF隨著系數(shù)b的增加先增加,然后趨于平穩(wěn)。對比圖12和圖13可以發(fā)現(xiàn),系數(shù)a對推力變化率ΔF的影響要遠遠強于系數(shù)b。

圖13 推力變化率隨系數(shù)b的變化曲線Fig.13 Curve of thrust change rate with b

當系數(shù)c發(fā)生變化時,零點時間隨系數(shù)a和系數(shù)b的改變而變化的曲線如圖14和圖15所示。圖中a*和b*表示零點時間t*=0時的系數(shù)a和系數(shù)b,稱之為系數(shù)a,b的臨界值。從圖中可以發(fā)現(xiàn),不同的系數(shù)c下,臨界值a*和b*幾乎不發(fā)生變化。且在ab*的區(qū)域內(nèi),4種情況的曲線幾乎重合。而在ab*區(qū)域內(nèi),推力是單調(diào)遞減的。說明增大空氣質(zhì)量流量并不會讓固體燃料超燃沖壓發(fā)動機推力變化的趨勢發(fā)生變化。系數(shù)c的變化在區(qū)域a>a*與b

圖14 不同系數(shù)c下的零點時間t*隨系數(shù)a增加而變化的曲線(b=0.04)Fig.14 Zero-point time varies with a under different c(b=0.04)

圖15 不同系數(shù)c下的零點時間t*隨系數(shù)b增加而變化的曲線(a=0.002)Fig.15 Zero-point time varies with b under different c(a=0.002)

當系數(shù)c發(fā)生變化時,推力變化率ΔF隨系數(shù)a和系數(shù)b的改變而變化的曲線如圖16和圖17所示。

圖16 不同系數(shù)c下的推力變化率隨系數(shù)a增加而變化的曲線(b=0.04)Fig.16 Thrust change rate varies with a under different c(b=0.04)

圖17 不同系數(shù)c下的推力變化率隨系數(shù)b增加而變化的曲線(a=0.002)Fig.17 Thrust change rate varies with b under different c(a=0.002)

從圖中可以看出,當系數(shù)c增大時,在相同的系數(shù)a或系數(shù)b處,推力變化率ΔF都降低了,這說明增大空氣質(zhì)量流量可以使固體燃料超燃沖壓發(fā)動機推力更加穩(wěn)定。且從圖中還可以看出,在推力變化率較大的系數(shù)a和系數(shù)b處,增大系數(shù)c降低推力變化率的效果更加明顯。

3 結(jié)論

建立了固體超燃沖壓發(fā)動機的零維推力模型,系統(tǒng)研究了燃料質(zhì)量流量與空氣流量對發(fā)動機推力的影響。

固體超燃沖壓發(fā)動機在工作過程中燃料質(zhì)量流量特性對推力的穩(wěn)定性有顯著影響,通過燃面退移速率與燃燒室構(gòu)型的匹配,理論上可以構(gòu)造出使推力變化幅度較小的發(fā)動機構(gòu)型。增大空氣流量更有利于固體超燃沖壓發(fā)動機推力的穩(wěn)定。然而,影響推力的因素不僅僅為燃料質(zhì)量流量和空氣流量,燃燒室內(nèi)流場結(jié)構(gòu)同樣對其有影響。因此,文中僅從理論上對燃料質(zhì)量流量和空氣流量對固體超燃沖壓發(fā)動機推力的影響有一個定性認識,即合適的燃面退移速率特性搭配合理的燃面構(gòu)型能夠讓固體燃料超燃沖壓發(fā)動機在規(guī)定工作時間內(nèi)推力變化幅度較小,推力更加穩(wěn)定。