燃?xì)獍l(fā)生器循環(huán)變推力發(fā)動機(jī)頻率響應(yīng)分析

黃夕軒,黃敏超,胡小平

(國防科技大學(xué) 空天科學(xué)學(xué)院,湖南 長沙 421000)

0 引言

液體火箭發(fā)動機(jī)推力可變技術(shù)是航天液體動力的重要發(fā)展范疇。采用變推力技術(shù)的液體火箭發(fā)動機(jī)可以實(shí)現(xiàn)航天運(yùn)輸系統(tǒng)的最佳推力控制,登月軟著陸的實(shí)現(xiàn)離不開變推力發(fā)動機(jī)[1-4],空間姿軌控任務(wù)使用變推力發(fā)動機(jī)可以提高控制系統(tǒng)的靈活性。對于導(dǎo)彈系統(tǒng),變推力發(fā)動機(jī)可以提高其突防能力。目前擁有高關(guān)注度的可重用發(fā)動機(jī)技術(shù)也要求具有推力可調(diào)性能的發(fā)動機(jī)[5]。

隨著可變推力發(fā)動機(jī)的發(fā)展,追求更高推力比,實(shí)現(xiàn)穩(wěn)定高效的變推過程成為未來的發(fā)展趨勢[6-7]。對于變推發(fā)動機(jī)的穩(wěn)定性能研究多以數(shù)值仿真進(jìn)行。陳啟智[8]等人對帶有可調(diào)噴注器的雙組元變推力液體火箭發(fā)動機(jī)出現(xiàn)的低頻振蕩現(xiàn)象進(jìn)行了研究,分析了設(shè)計(jì)參數(shù)對穩(wěn)定性的影響。文獻(xiàn)[9]通過動力學(xué)模型的方法分別描述了以流量信號反饋和位移信號反饋兩種方式的擠壓式推進(jìn)劑供應(yīng)變推力液體火箭發(fā)動機(jī)的傳遞函數(shù)以及狀態(tài)空間。在狀態(tài)空間的基礎(chǔ)上對兩種反饋方式的穩(wěn)定條件進(jìn)行了分析,并且通過試車數(shù)據(jù)驗(yàn)證了計(jì)算得到的系統(tǒng)的穩(wěn)定區(qū)域與不穩(wěn)定區(qū)域,證明了采用動力學(xué)模型的方法尋找系統(tǒng)穩(wěn)定性區(qū)間的可行性。Rodden[10]等人對可調(diào)噴注器進(jìn)行了數(shù)學(xué)建模并繪制了噴注器調(diào)節(jié)系統(tǒng)流程圖,結(jié)合頻域分析,給出了噴注器的設(shè)計(jì)建議。文獻(xiàn)[11]對流量調(diào)節(jié)器及泵壓式供應(yīng)系統(tǒng)建立了集中參數(shù)線性化動力學(xué)模型,研究了系統(tǒng)頻率響應(yīng)特性以及系統(tǒng)的固有穩(wěn)定性,通過時域與頻域的對比驗(yàn)證了頻域分析的正確性,因此本文使用的集中參數(shù)線性化動力學(xué)模型是可以應(yīng)用于頻域分析的。劉景華[12]采用混合仿真方法,建立了過氧化氫發(fā)動機(jī)穩(wěn)定性分析模型,得到脈沖擾動下發(fā)動機(jī)響應(yīng)特性以及燃燒時滯、噴注壓降等參數(shù)對穩(wěn)定性影響?？傮w看來,對于動態(tài)特性的研究主要是通過數(shù)學(xué)方法建立狀態(tài)空間的模型,進(jìn)而分析頻率響應(yīng)特性,且大部分工作是對于擠壓式系統(tǒng)進(jìn)行的。目前對于泵壓式變推力發(fā)動機(jī)的變推力過程的穩(wěn)定性研究較少。

泵壓式液體火箭發(fā)動機(jī)是使用渦輪泵在較低的壓力下從推進(jìn)劑貯箱中接收液體推進(jìn)劑,并以所需的流量和噴射壓力將其供給發(fā)動機(jī)燃燒室。渦輪的動力來自高壓氣體的膨脹,而高壓氣體通常是由渦輪泵輸送的推進(jìn)劑的燃燒產(chǎn)物。泵壓式火箭發(fā)動機(jī)應(yīng)用范圍廣泛:如20世紀(jì)60年代為土星五號運(yùn)載火箭開發(fā)的F-1和J-2發(fā)動機(jī); 20世紀(jì)70年代航天飛機(jī)使用的主引擎(SSME),使航天器首次實(shí)現(xiàn)可重復(fù)使用[13]。因此,對于泵壓式液體火箭發(fā)動機(jī)頻率響應(yīng)方面的研究進(jìn)行完善是必要的。

通過進(jìn)行線性分析,可以了解系統(tǒng)狀態(tài)空間的特征,可以反映時域仿真中不能明確體現(xiàn)的頻率響應(yīng)特性,進(jìn)一步了解系統(tǒng)高、低頻率下的幅相頻特性。因此,本文在模擬了發(fā)動機(jī)啟動過程及變推力過程的主要部件時域響應(yīng)之外,用快速傅里葉方法提取了發(fā)動機(jī)的變推力過程推力室壓力的主要頻率,對于加入該頻率擾動信號的推力室壓力進(jìn)行了線性分析,得到其頻率響應(yīng)特性。

1 發(fā)動機(jī)主要組件建模

文獻(xiàn)[14]中對液氧煤油發(fā)動機(jī)起動仿真結(jié)果與試車數(shù)據(jù)進(jìn)行了對比,發(fā)現(xiàn)2組數(shù)據(jù)的參數(shù)變化趨勢與過渡時間吻合性較好,驗(yàn)證了模型的正確性與可行性。本文使用的物理模型與文獻(xiàn)[14]相同,故不再進(jìn)行驗(yàn)證。

1.1 仿真系統(tǒng)示意圖

本文所研究燃?xì)獍l(fā)生器循環(huán)發(fā)動機(jī)的系統(tǒng)如圖1所示,該系統(tǒng)由燃料貯箱、氧化劑貯箱、渦輪、離心泵、燃?xì)獍l(fā)生器、混合比調(diào)節(jié)器、推力室閥門、管道及節(jié)流孔板所構(gòu)成。

圖1 泵壓式燃?xì)獍l(fā)生器循環(huán)發(fā)動機(jī)系統(tǒng)示意圖Fig.1 Schematic diagram of gas generator cycle engine system

1.2 管路模型

發(fā)動機(jī)管路采用集中參數(shù)法構(gòu)建模型。集中參數(shù)法是將管路按照流阻、流感和流容3類基本元件進(jìn)行有限元分割的方法。在低頻動力學(xué)范圍內(nèi)廣泛使用,因該方法的前提假設(shè)為小擾動、線性化,故不適合于中高頻的震蕩特性的研究[15]。

在該方法中,將主要考慮流體慣性的單元稱為慣性單元(I),主要考慮局部阻力的單元稱為阻性單元(R),主要考慮流體壓縮性的單元設(shè)置為容性單元(C),同時考慮慣性和阻性的單元稱為運(yùn)動單元(IR)。

1.3 離心泵模型

建立發(fā)動機(jī)離心泵模型過程中使用了泵揚(yáng)程及阻力矩的性能數(shù)據(jù)。離心泵性能數(shù)據(jù)由文獻(xiàn)[16-17]中實(shí)驗(yàn)數(shù)據(jù)擬合得到。泵的特性曲線是將一維數(shù)據(jù)表擬合成關(guān)于自變量θ的函數(shù)而繪成,其中,θ為反映泵工作特性的特征角,其定義為

θ=π+arctan(ν/n)

(1)

式中ν和n分別為流量和轉(zhuǎn)速的無量綱特征數(shù),其定義為

(2)

揚(yáng)程和轉(zhuǎn)矩的無量綱參數(shù)h和β定義如下

(3)

式中：H為實(shí)際揚(yáng)程;HR為參考揚(yáng)程, m;ηR為泵參考效率;ρ為泵入口和出口的工質(zhì)平均密度,kg/m3。

1.4 燃燒室模型

由于燃燒過程對發(fā)動機(jī)低頻特性影響較小,燃燒零維模型被廣泛用于發(fā)動機(jī)動態(tài)特性分析[20-21]。燃燒室模型采用基于燃燒時滯的燃燒動力學(xué)模型,且假設(shè)燃燒時滯是一個常數(shù),燃?xì)庠谌紵沂蔷鶆蚍植嫉?各參數(shù)不隨空間坐標(biāo)變化。

質(zhì)量微分方程

(4)

能量微分方程

(5)

混合比微分方程

(6)

式中：ρ為燃燒室平均密度;p為燃燒室室壓;kc為燃燒室混合比(O/F);γ為定壓比熱與定容比熱之比,是關(guān)于時間的函數(shù);qmf為燃料質(zhì)量流量;qmo為氧化劑質(zhì)量流量;V為燃燒室體積;hf和ho分別為燃料和氧化劑的比焓;R為燃?xì)鈿怏w常數(shù);T為燃?xì)鉁囟取?/p>

2 結(jié)果與討論

2.1 工況說明

本文設(shè)置了2個工況：其中每個工況對應(yīng)的燃料及氧化劑入口流量如表1所示,入口流量調(diào)節(jié)是通過調(diào)節(jié)氧化劑主閥與燃料主閥開度實(shí)現(xiàn)的。

表1 設(shè)置工況說明

本文設(shè)置了2個工況轉(zhuǎn)換點(diǎn),分別位于仿真時間的4 s由高工況轉(zhuǎn)變?yōu)榈凸r和7 s由低工況轉(zhuǎn)變?yōu)楦吖r。切換工況時推進(jìn)劑主閥開度改變,同時主系統(tǒng)工況切換閥和副系統(tǒng)工況切換閥切換工作狀態(tài)。

2.2 發(fā)動機(jī)啟動及變工況仿真結(jié)果

圖2～圖4分別給出啟動后經(jīng)過工況轉(zhuǎn)換的發(fā)動機(jī)系統(tǒng)主要部件主要參數(shù)隨時間變化曲線。推力室室壓在啟動時出現(xiàn)超調(diào),在高工況時穩(wěn)定在4.2 MPa,低工況時穩(wěn)定在1.7 MPa,在2個工況切換點(diǎn)均未出現(xiàn)超調(diào)。離心泵轉(zhuǎn)速啟動后在升至36 000 r/min中無超調(diào),在轉(zhuǎn)換為低工況后有極小的超調(diào),之后穩(wěn)定在37 500 r/min。渦輪功率在兩個工況切換點(diǎn)均出現(xiàn)了超調(diào),且超調(diào)量較大,反映了副系統(tǒng)受變工況的影響較大。

圖2 推力室室壓隨時間變化曲線Fig.2 Thrust chamber pressure vs time

圖3 離心泵轉(zhuǎn)速隨時間變化曲線Fig.3 Centrifugal pump speed vs time

圖4 渦輪功率隨時間變化曲線Fig.4 Turbine power vs time

2.3 系統(tǒng)特征頻率與阻尼系數(shù)分析

在燃?xì)獍l(fā)生器循環(huán)液體火箭發(fā)動機(jī)系統(tǒng)中,燃?xì)獍l(fā)生器及推力室參數(shù)被作為主要被控對象。將燃?xì)獍l(fā)生器與推力室壓力等參數(shù)設(shè)置為觀測變量,分別在2種工況達(dá)到穩(wěn)定狀態(tài)時考察系統(tǒng)的特征頻率及阻尼系數(shù)。根據(jù)各特征頻率下觀測變量的模態(tài)振型,系統(tǒng)的特征頻率可分為2類：對應(yīng)于燃?xì)獍l(fā)生器的特征頻率與對應(yīng)于4臺并聯(lián)的推力室的特征頻率。由表2可以看出,4臺并聯(lián)的推力室特征頻率非常接近。

表2 高工況系統(tǒng)阻尼系數(shù)及頻率

Tab.2 Eigenvalues and frequency of system at high operating condition

特征頻率/Hz阻尼比阻尼系數(shù)實(shí)部阻尼系數(shù)虛部6.476 01.00-40.691 40.0056.860 21.00-357.264 00.0056.861 01.00-357.268 00.0056.861 51.00-357.271 00.0056.861 61.00-357.272 00.00110.777 01.00-696.036 00.00110.778 01.00-696.042 00.00110.779 01.00-696.045 00.00110.779 01.00-696.046 00.00

表2所示的高工況9個特征頻率中,對應(yīng)于燃?xì)獍l(fā)生器壓力的特征頻率是6.476 Hz。4個在56.86 Hz附近的和4個在110.78 Hz附近的特征頻率對應(yīng)4臺推力室的特征頻率。高工況下系統(tǒng)所有阻尼系數(shù)均為實(shí)數(shù),且實(shí)部為負(fù)數(shù)。這表明系統(tǒng)是穩(wěn)定的。

表3中所示的低工況9個特征頻率中,燃?xì)獍l(fā)生器特征頻率為17.997 Hz。其中4個為46.563 1 Hz和4個在94.896 4 Hz附近的特征頻率為4臺推力室的特征頻率。相對于高工況,低工況下系統(tǒng)所有特征頻率及阻尼系數(shù)均發(fā)生了改變,且阻尼系數(shù)全部為負(fù)實(shí)數(shù)。這是由于變工況時,主閥及主、副系統(tǒng)工況切換閥工作狀態(tài)改變,導(dǎo)致通流面積改變,因而引起管道局部阻力的變化,而流阻的變化必將改變系統(tǒng)的動力學(xué)方程參數(shù),進(jìn)而影響系統(tǒng)狀態(tài)空間,最終表現(xiàn)為特征頻率與阻尼系數(shù)的改變。

表3 低工況系統(tǒng)阻尼系數(shù)及頻率

Tab.3 Eigenvalues and frequency of system atlow operating condition

特征頻率/Hz阻尼比阻尼系數(shù)實(shí)部阻尼系數(shù)虛部17.997 01.00-113.078 30.0046.563 11.00-292.565 00.0046.563 11.00-292.565 00.0046.563 11.00-292.565 00.0046.563 11.00-292.565 00.0094.896 41.00-596.252 00.0094.896 41.00-596.252 00.0094.896 41.00-596.252 00.0094.896 41.00-596.252 00.00

總體看來,系統(tǒng)一階特征頻率對應(yīng)的阻尼系數(shù)最小,表明系統(tǒng)對該頻率較敏感。該特征頻率是燃?xì)獍l(fā)生器壓力的特征頻率,說明由于副系統(tǒng)對低頻信號更敏感。

由阻尼系數(shù)變化可知系統(tǒng)由高工況變?yōu)榈凸r后系統(tǒng)對低頻敏感性減小,對高頻敏感性增加,但影響系統(tǒng)穩(wěn)定性的主要信號仍然是低頻信號。

2.4 頻率響應(yīng)分析

本文對第一個變工況點(diǎn)后的變推力段推力室室壓進(jìn)行快速傅里葉變換,如圖5所示,得到其最主要的頻率成分在10 Hz以內(nèi),在5.34 Hz達(dá)到了峰值,驗(yàn)證了上文對特征頻率的分析。為了研究該頻率下的入口壓力震蕩對系統(tǒng)穩(wěn)定性的影響,將擾動設(shè)置為一正弦波,因此推進(jìn)劑入口壓力波方程為

(7)

式中：A=10%PIn0；f=5.34 Hz；PIn0為推進(jìn)劑入口壓力；φ為相位差；PIn0為0.49 MPa；φ=0°。本文選取4.01 s,4.02 s,4.03 s,4.04 s,4.05 s 5個時間點(diǎn),將推力室壓力設(shè)為輸出信號,繪制以入口壓力為輸入信號的Nyquist圖,結(jié)果如圖6～圖10所示。

圖5 變推力過程室壓FFT分析 Fig.5 FFT analysis of chamber pressure of thrustvarying process

圖6 4.01 s推力室壓力的Nyquist曲線Fig.6 Nyquist curve of thrust chamber pressure at 4.01 s

圖7 4.02 s推力室壓力的Nyquist圖Fig.7 Nyquist curve of thrust chamber pressure at 4.02 s

圖8 4.03 s推力室壓力的Nyquist圖Fig.8 Nyquist curve of thrust chamber pressure at 4.03 s

圖9 4.04 s推力室壓力的Nyquist圖Fig.9 Nyquist curve of thrust chamber pressure at 4.04 s

圖10 4.05 s推力室壓力的Nyquist圖Fig.10 Nyquist curve of thrust chamber pressure at 4.05 s

Nyquist圖反映了相對于控制信號(入口壓力),輸出信號(推力室壓力)的相位差以及幅值比,可以觀察到在高頻區(qū)曲線的行為趨向于二階震蕩環(huán)節(jié),但由于幅值較小沒有影響系統(tǒng)穩(wěn)定性。

根據(jù)Nyquist判據(jù),在開環(huán)系統(tǒng)穩(wěn)定時,Nyquist曲線不包圍(-1,j0)點(diǎn),閉環(huán)系統(tǒng)也穩(wěn)定。圖6～圖10中Nyquist曲線均未包圍(-1,j0)點(diǎn)。因此將該系統(tǒng)在所設(shè)置擾動下變工況過程穩(wěn)定。

系統(tǒng)的幅相頻特性曲線的最大幅值小于1時可知系統(tǒng)的相位不會改變系統(tǒng)的穩(wěn)定性,系統(tǒng)的幅值裕度可以由曲線與負(fù)實(shí)軸交界處幅值求得。各測點(diǎn)的穩(wěn)定裕度如表4所示。在系統(tǒng)變工況過程中, 4.03 s低頻段出現(xiàn)了幅值過大的現(xiàn)象,且在相角滯后大于30.856°后會出現(xiàn)系統(tǒng)不穩(wěn)定現(xiàn)象。4.05 s低頻段時出現(xiàn)了較低的穩(wěn)定裕度,目前系統(tǒng)的開環(huán)增益增大3.404 4倍則系統(tǒng)失穩(wěn)。由此可知變工況的過程中低頻擾動信號對于穩(wěn)定性影響不容忽略。

表4 各觀測點(diǎn)穩(wěn)定裕度

3 結(jié)論

本文對泵壓式雙組元燃?xì)獍l(fā)生器循環(huán)液體火箭發(fā)動機(jī)系統(tǒng)變工況頻率響應(yīng)特性進(jìn)行了研究,可以得到以下一些結(jié)論：

1)本文所研究的發(fā)動機(jī)系統(tǒng)在高工況與低工況下的特征頻率存在明顯差別：系統(tǒng)的一階特征頻率(6.476 Hz和17.996 Hz)均為燃?xì)獍l(fā)生器壓力對應(yīng)的特征頻率,屬于低頻范圍,是影響系統(tǒng)穩(wěn)定性的主要頻率。并聯(lián)的4臺推力室特征頻率十分接近,特征頻率相對較高,且阻尼系數(shù)較大,對系統(tǒng)影響相對較小。

2)由高工況下轉(zhuǎn)變?yōu)榈凸r過程中,該發(fā)動機(jī)系統(tǒng)阻尼系數(shù)改變。由阻尼系數(shù)變化可以得出：系統(tǒng)由高工況變?yōu)榈凸r后系統(tǒng)對低頻敏感性減小,對高頻敏感性增加,但影響系統(tǒng)穩(wěn)定性的主要信號仍然是低頻信號。因此對于該發(fā)動機(jī)系統(tǒng)推力較大的工況,應(yīng)當(dāng)特別注意頻率為10 Hz以內(nèi)的信號的影響。

3)變工況過程中, Nyquist曲線反映了高頻段系統(tǒng)偏向于二階震蕩環(huán)節(jié)但幅值小不影響系統(tǒng)穩(wěn)定性。低頻段曲線反映出的加入低頻壓力擾動信號可以造成推力室壓力的穩(wěn)定裕度變化,出現(xiàn)了較小的幅值裕度(3.404 4)。由于出現(xiàn)相角裕度小于40°的點(diǎn),說明系統(tǒng)若存在相角延遲容易造成系統(tǒng)失穩(wěn)的現(xiàn)象,進(jìn)一步說明變工況過程容易受到低頻壓力入口擾動影響而失穩(wěn)。系統(tǒng)在工作時,為了維持更高的系統(tǒng)穩(wěn)定性裕度應(yīng)該盡可能避免低頻入口壓力擾動信號的出現(xiàn)。