張 凌,王德友,李世鵬,王寧飛

(北京理工大學(xué) 宇航學(xué)院,北京 100081)

0 引言

固體火箭發(fā)動機內(nèi)彈道學(xué)的基本任務(wù)是在發(fā)動機各種工作條件下計算燃燒室壓強隨時間和空間的變化規(guī)律。零維內(nèi)彈道計算方法假定燃燒室內(nèi)參數(shù)均勻一致,與一維及多維內(nèi)彈道計算方法相比更加簡便、直觀,是工程上應(yīng)用較多的求解方法[1-2]。發(fā)動機的內(nèi)彈道特性直接決定著整機性能,進而影響著火箭導(dǎo)彈的飛行特性。因此,進一步提高內(nèi)彈道預(yù)示的精度,對保障飛行器工作可靠性有著重要意義。

工程上通常根據(jù)單臺或多臺發(fā)動機的實測數(shù)據(jù)和燃面幾何退移規(guī)律,對燃速特性等參數(shù)進行辨識,以獲得更精確的內(nèi)彈道性能預(yù)示結(jié)果。目前,業(yè)內(nèi)學(xué)者針對發(fā)動機內(nèi)彈道性能精確預(yù)示問題,從理論算法修正、參數(shù)辨識方法改進等角度開展了研究。如劉楊等[3]基于對發(fā)動機實測數(shù)據(jù)處理誤差的理論分析,提出了修正的燃速系數(shù)處理方法,該方法可提高具有較高壓強比的發(fā)動機內(nèi)彈道預(yù)示精度。李春艷等[4]針對短時間工作固體發(fā)動機特點,從理論角度對現(xiàn)有內(nèi)彈道計算方法進行了適應(yīng)性修改,所得預(yù)示結(jié)果與實測數(shù)據(jù)相比具有較高精度。

在辨識參數(shù)選擇方面,CAVALLINI等[5-6]通過研究Vega固體運載火箭發(fā)動機實際飛行性能影響因素,發(fā)現(xiàn)噴管喉部燒蝕率對推力和燃燒室壓強的影響較大。BIANCHI等[7]采取數(shù)值仿真方法研究噴管喉部燒蝕率和燃速變化對燃燒室壓強的影響,發(fā)現(xiàn)這兩個參數(shù)的改變對壓強-時間曲線影響十分明顯。李曉斌等[8]基于發(fā)動機工作過程的內(nèi)彈道數(shù)學(xué)模型及實測燃燒室壓強時變數(shù)據(jù),通過參數(shù)辨識方法可對燃速、喉襯燒蝕率等參數(shù)進行計算。

在參數(shù)辨識方法改進方面,由于最小二乘法、最大似然法等傳統(tǒng)辨識方法常存在初始點敏感等問題,因而如遺傳算法等具有初值穩(wěn)健性和全局搜索性的方法在內(nèi)彈道參數(shù)辨識領(lǐng)域逐漸得到深入應(yīng)用。樊超等[9]將遺傳算法用于固體火箭發(fā)動機參數(shù)辨識,獲得了推進劑燃速模型和喉徑變化模型的全局最優(yōu)辨識值,所得內(nèi)彈道性能預(yù)示結(jié)果與試驗測量值吻合良好。莊建華等[10]對發(fā)動機貯存老化內(nèi)彈道性能的影響因素進行參數(shù)建模,基于拉丁超立方設(shè)計方法進行了瞬態(tài)參數(shù)辨識,結(jié)合非定常內(nèi)彈道模型對發(fā)動機貯存老化內(nèi)彈道進行了預(yù)示。蒲曉航等[11]基于部分實測壓強數(shù)據(jù),通過遺傳算法逆向反算獲得預(yù)示參數(shù),再結(jié)合正向計算實現(xiàn)了內(nèi)彈道的預(yù)示自修正,對復(fù)雜環(huán)境剖面下的發(fā)動機內(nèi)彈道性能進行了準確預(yù)示。孫瑞陽等[12]建立了基于經(jīng)典AlexNet卷積神經(jīng)網(wǎng)絡(luò)的內(nèi)彈道參數(shù)辨識方法,對長期貯存發(fā)動機推進劑裝藥的燃速特性參數(shù)進行了辨識。張小平等[13]建立了高能固體推進劑的配方表征方法,利用遺傳神經(jīng)網(wǎng)絡(luò)(GA-BP)對推進劑高壓燃燒性能及其變化規(guī)律進行了預(yù)示。

傳統(tǒng)發(fā)動機的實測數(shù)據(jù)處理常通過燃燒時間內(nèi)的平均壓強和藥柱平均燃速計算獲得燃速系數(shù)等參數(shù)并反算燃面,而實際上平均壓強與平均燃速間并不符合燃速模型對應(yīng)關(guān)系[3],因而計算存在一定誤差。而目前的研究工作大多采用遺傳算法和神經(jīng)網(wǎng)絡(luò)等方法對關(guān)鍵內(nèi)彈道參數(shù)進行辨識,并沒有根據(jù)實測數(shù)據(jù)對藥柱燃面以及內(nèi)彈道參數(shù)進行修正與二次辨識,內(nèi)彈道參數(shù)辨識以及性能預(yù)示的精度還可以進一步提高。因此,有必要對關(guān)鍵內(nèi)彈道參數(shù)開展二次辨識,得到更加精確的辨識結(jié)果,進而提高發(fā)動機內(nèi)彈道預(yù)示精度。

本文提出了一種發(fā)動機內(nèi)彈道多參數(shù)辨識及性能精準預(yù)示方法,根據(jù)發(fā)動機實測壓強數(shù)據(jù),通過遺傳算法對燃速系數(shù)、壓強指數(shù)、推進劑密度和噴管喉襯燒蝕速率系數(shù)開展初步辨識,基于辨識結(jié)果修正燃面-肉厚關(guān)系,在此基礎(chǔ)上開展二次參數(shù)辨識,最終確定關(guān)鍵內(nèi)彈道參數(shù)。根據(jù)第二次參數(shù)辨識結(jié)果和修正后的燃面-肉厚關(guān)系,通過正向預(yù)示可獲得發(fā)動機的內(nèi)彈道性能。以美國C1xb型發(fā)動機[14-15]和國內(nèi)某型長通道發(fā)動機[16]為案例進行計算分析,驗證了該方法內(nèi)彈道性能預(yù)示的精度和有效性,相關(guān)結(jié)果可為固體火箭發(fā)動機的優(yōu)化設(shè)計提供技術(shù)參考。

1 內(nèi)彈道預(yù)示模型

本文以零維內(nèi)彈道模型[17-18]為基礎(chǔ)進行發(fā)動機性能正向預(yù)示,該模型有如下假設(shè):(1)燃氣為理想氣體,氣體流動為零維定常流動;(2)裝藥點燃過程瞬時完成;(3)噴管喉襯燒蝕速率與燃燒室壓強正相關(guān)。在考慮點火藥和噴管喉襯燒蝕作用對內(nèi)彈道計算影響的情況下,結(jié)合燃燒室內(nèi)質(zhì)量守恒定律、能量守恒定律、燃燒室內(nèi)體積守恒方程、燃速方程以及氣體狀態(tài)方程,推導(dǎo)出的內(nèi)彈道微分方程組如下:

(1)

其中,

(2)

(3)

(4)

(5)

2 多參數(shù)辨識及精準預(yù)示方法

考慮到發(fā)動機試驗成本問題,基于一次或幾次試車數(shù)據(jù)開展發(fā)動機內(nèi)彈道關(guān)鍵參數(shù)的精確辨識研究,對降低固體發(fā)動機設(shè)計成本和提高設(shè)計效率有著重要意義?；诠こ讨械膬?nèi)彈道參數(shù)辨識方法,本文提出一種發(fā)動機內(nèi)彈道多參數(shù)辨識及預(yù)示方法,可提高內(nèi)彈道預(yù)示精度。

2.1 工程中內(nèi)彈道參數(shù)辨識方法

(6)

(7)

(8)

式中W為裝藥肉厚;t為燃燒時間;Ab為燃面面積;At為噴管喉部面積;C*為特征速度。

2.2 多參數(shù)辨識及精準預(yù)示方法流程

上述平均燃速和平均壓強并不具有實際對應(yīng)關(guān)系,計算存在一定誤差,在內(nèi)彈道預(yù)示精度上仍有較大優(yōu)化空間。本文提出了一種內(nèi)彈道多參數(shù)辨識及精準預(yù)示方法,該方法是對工程辦法的改進。以單臺發(fā)動機的實測壓強數(shù)據(jù)為基礎(chǔ),在發(fā)動機穩(wěn)定工作段選取部分時間段進行計算,基于所建立辨識修正方法獲得燃速模型參數(shù)和燃面-肉厚數(shù)據(jù)等內(nèi)彈道預(yù)示所需全部關(guān)鍵參數(shù),進行內(nèi)彈道性能預(yù)示,該方法能提高內(nèi)彈道預(yù)示精度,具體流程如圖1所示。

圖1 多參數(shù)辨識及精準預(yù)示方法流程Fig.1 Multi-parameter identification and prediction method process

從發(fā)動機實測壓強數(shù)據(jù)出發(fā),基于GA遺傳算法對燃速系數(shù)a、壓強指數(shù)n、推進劑密度ρp及噴管喉襯燒蝕速率系數(shù)as進行第一次參數(shù)辨識。通過參數(shù)辨識結(jié)果和燃面-壓強公式對燃面-肉厚曲線進行反算,進而更新多參數(shù)辨識的計算模型,再通過所得燃面數(shù)據(jù)和遺傳算法開展第二次參數(shù)辨識。根據(jù)二次辨識結(jié)果對初次辨識參數(shù)進行修正,在第一章介紹的內(nèi)彈道預(yù)示模型的基礎(chǔ)上,根據(jù)所得參數(shù)及燃面數(shù)據(jù)對發(fā)動機內(nèi)彈道性能開展正向預(yù)示。

2.3 基于遺傳算法的內(nèi)彈道多參數(shù)辨識

內(nèi)彈道多參數(shù)辨識主要根據(jù)發(fā)動機工作時間內(nèi)的實測壓強數(shù)據(jù)進行逆向反算,該工作基于GA遺傳算法完成,對燃速系數(shù)a、壓強指數(shù)n、推進劑密度ρp及噴管喉襯燒蝕速率系數(shù)as四個參數(shù)進行辨識,具體步驟如下:

(1)根據(jù)實測數(shù)據(jù)信息及設(shè)計經(jīng)驗,選定時間間隔Δt,作為計算兩條曲線擬合程度的最小時間微元;

(2)給定參數(shù)初值范圍,其中燃速系數(shù)a和壓強指數(shù)n還需滿足燃速模型關(guān)系,故需施加額外限制條件;

(3)在給定初值范圍內(nèi),隨機給出N組數(shù)據(jù)(a,n,ρp,as)作為模型初始群體,每組數(shù)據(jù)為初始群體的個體;

(4)根據(jù)零維內(nèi)彈道理論建立正向預(yù)示數(shù)學(xué)模型,基于任意一組數(shù)據(jù)可進行一次內(nèi)彈道性能預(yù)示,即N組數(shù)據(jù)可獲得N條壓強曲線;

(5)由于發(fā)動機穩(wěn)定工作段數(shù)據(jù)更具代表性,故在穩(wěn)定工作段內(nèi)以Δt為間隔總計選取M個節(jié)點,建立內(nèi)彈道預(yù)示模型的適應(yīng)度函數(shù),其為M個節(jié)點處預(yù)示壓強與實測壓強的差值組成的行向量,該行向量的范數(shù)作為衡量該組數(shù)據(jù)適應(yīng)度高低的判據(jù);

(6)淘汰適應(yīng)度低的初始個體,保留適應(yīng)度高的初始個體,在保留個體間進行數(shù)據(jù)的交換和變異,形成初始群體的子代,再次通過適應(yīng)度函數(shù)對新群體進行選擇;

(7)重復(fù)上步操作直至結(jié)果收斂,將適應(yīng)度最高的數(shù)據(jù)代入內(nèi)彈道預(yù)示模型中,即可獲得與實測數(shù)據(jù)吻合度最高的壓強曲線。

上述為初次參數(shù)辨識的步驟?；诔醮螀?shù)辨識結(jié)果進行燃面反算,根據(jù)反算燃面數(shù)據(jù)更新內(nèi)彈道預(yù)示數(shù)學(xué)模型,再進行二次參數(shù)辨識修正。第二次辨識的參數(shù)初值范圍可適當縮小,以加快收斂速度并提高效率。

3 案例計算分析

為驗證內(nèi)彈道性能預(yù)示方法的準確性和有效性,本文選取C1xb型發(fā)動機和一種長通道發(fā)動機作為研究案例。

3.1 C1xb型發(fā)動機案例

C1xb型發(fā)動機的裝藥組成為79% AP和21% HTPB的復(fù)合推進劑,該發(fā)動機結(jié)構(gòu)參數(shù)和試車數(shù)據(jù)較為詳盡[15],許多關(guān)于發(fā)動機燃燒穩(wěn)定性和內(nèi)流場特性的研究均選取該發(fā)動機作為驗證案例。該型發(fā)動機的具體參數(shù)如表1所示,結(jié)構(gòu)如圖2所示。

表1 C1xb發(fā)動機參數(shù)Table 1 Parameters of C1xb SRM

圖2 C1xb發(fā)動機結(jié)構(gòu)圖Fig.2 Schematic of C1xb SRM

3.1.1 第一次參數(shù)辨識

(1)取最小時間微元Δt=0.01 s進行曲線擬合度對比計算;

(2)在穩(wěn)定工作段內(nèi),以Δt為時間間隔選取300個點作為曲線擬合度計算時間節(jié)點;

(3)計算所有節(jié)點上預(yù)測壓強值與實測壓強值的差向量,將該行向量的Euclid范數(shù)作為適應(yīng)度函數(shù)。

C1xb發(fā)動機實測壓強曲線存在初始尖峰,根據(jù)工程經(jīng)驗可知其為點火藥與裝藥共同燃燒導(dǎo)致產(chǎn)氣量劇增所致[17-18]。點火過程影響因素較為復(fù)雜,但其對內(nèi)彈道曲線整體趨勢影響并不顯著。辨識計算中取發(fā)動機穩(wěn)定工作部分時間段作為計算區(qū)間,避開初始點火過程和壓強下降拖尾段,以提高多參數(shù)辨識的精度。

在第一次多參數(shù)辨識的過程中,設(shè)定4項待辨識參數(shù)的初始取值范圍分別為3.5 mm/(s·MPan)≤a≤5 mm/(s·MPan),0.4≤n≤0.5,1550 kg/m3≤ρp≤1650 kg/m3,0.1 mm/(s·MPa)≤as≤0.2 mm/(s·MPa)?；谶z傳算法進行多參數(shù)辨識,在經(jīng)過多次迭代計算后結(jié)果收斂,求解過程如圖3所示。初次參數(shù)辨識結(jié)果為:燃速系數(shù)a=4.48 mm/(s·MPan),壓強指數(shù)n=0.487 6,推進劑密度ρp=1556 kg/m3,噴管喉襯燒蝕速率系數(shù)as=0.138 7 mm/(s·MPa)。

圖3 第一次遺傳算法優(yōu)化求解過程Fig.3 Calculation progress of first multi- parameter identification

由于原始文獻未給出發(fā)動機點火裝置參數(shù),且點火過程對宏觀內(nèi)彈道特性影響較小,故將點火藥假定為球形顆粒黑火藥,再基于遺傳算法和考慮點火過程的內(nèi)彈道預(yù)示模型對點火藥質(zhì)量和顆粒粒徑進行辨識。該過程所需其余參數(shù)采用初次參數(shù)辨識的結(jié)果,最終獲得點火藥質(zhì)量為2.5 g,粒徑為1.5 mm。根據(jù)參數(shù)辨識結(jié)果進行內(nèi)彈道預(yù)示,可獲得初步預(yù)示的壓強-時間曲線,如圖4所示。

圖4 第一次辨識計算結(jié)果與實測數(shù)據(jù)對比Fig.4 Comparison between first prediction and measured pressure data

3.1.2 燃面反算及第二次參數(shù)辨識

基于初次參數(shù)辨識結(jié)果,根據(jù)燃面-壓強公式進行燃面反算。為避免點火過程對反算燃面數(shù)據(jù)產(chǎn)生影響,在計算時不考慮點火藥燃燒的影響,可獲得如圖5所示的反算燃面-肉厚曲線。

圖5 反算燃面-肉厚曲線Fig.5 Inverse calculation of burning area-web curve

在初次辨識基礎(chǔ)上開展二次辨識,基于反算燃面數(shù)據(jù)對已有數(shù)學(xué)模型進行更新。為加速收斂,第二次辨識的參數(shù)初值范圍相較第一次辨識適當縮小。設(shè)定4項待辨識參數(shù)的初始取值范圍分別為4.2 mm/(s·MPan)≤a≤4.8 mm/(s·MPan),0.45≤n≤0.5,1550 kg/m3≤ρp≤1600 kg/m3,0.12 mm/(s·MPa)≤as≤0.2 mm/(s·MPa)。如圖6所示,經(jīng)過多次迭代計算后結(jié)果收斂,第二次參數(shù)辨識結(jié)果為:燃速系數(shù)a=4.648 mm/(s·MPan),壓強指數(shù)n=0.455 1,推進劑密度ρp=1589 kg/m3,噴管喉襯燒蝕率系數(shù)as=0.142 3 mm/(s·MPa)。根據(jù)反算燃面數(shù)據(jù)和第二次參數(shù)辨識結(jié)果再次進行內(nèi)彈道預(yù)示,可獲得圖7所示的修正壓強-時間曲線。由圖7可知,第二次性能預(yù)示精度較第一次更高,且預(yù)示結(jié)果與實測數(shù)據(jù)高度吻合。

圖6 第二次遺傳算法優(yōu)化求解過程Fig.6 Calculation progress of second multi- parameter identification

圖7 第二次辨識計算結(jié)果與實測數(shù)據(jù)對比Fig.7 Comparison between second prediction and measured pressure data

3.2 某型長通道發(fā)動機案例

測試案例二為某型長通道固體發(fā)動機,該發(fā)動機裝藥為端面包覆的內(nèi)外燃管形裝藥,采用雙基推進劑并以黑火藥為點火藥,其相關(guān)參數(shù)較為詳盡[16],具體結(jié)構(gòu)參數(shù)如表2所示,結(jié)構(gòu)圖如圖8所示。

表2 長通道發(fā)動機參數(shù)Table 2 Parameters of the long passage SRM

圖8 長通道發(fā)動機結(jié)構(gòu)圖Fig.8 Schematic of the long passage SRM

應(yīng)用上述介紹的內(nèi)彈道多參數(shù)辨識方法對該長通道發(fā)動機進行第一次參數(shù)辨識,第一次參數(shù)辨識結(jié)果為:燃速系數(shù)a=3.01 mm/(s·MPan),壓強指數(shù)n=0.444 1,推進劑密度ρp=1563 kg/m3,噴管喉襯燒蝕率系數(shù)as=0.112 2 mm/(s·MPa)。用第一次多參數(shù)辨識結(jié)果進行燃面反算,根據(jù)內(nèi)彈道預(yù)示模型進行性能預(yù)示,計算得到如圖9所示的燃燒室壓強預(yù)示曲線和如圖10所示反算燃面曲線。由于燃面數(shù)據(jù)計算采用式(8)所示的工程方法,故上升段和下降段與實際燃面變化的物理過程存在一定差異。

圖9 長通道發(fā)動機第一次辨識預(yù)示曲線Fig.9 First prediction of the long passage SRM

圖10 長通道發(fā)動機反算燃面-肉厚曲線Fig.10 Inverse calculation of burning area-web curve of the long passage SRM

根據(jù)反算燃面數(shù)據(jù)更新模型后,進行第二次多參數(shù)辨識,得到最終參數(shù)辨識結(jié)果為:燃速系數(shù)a=3.09 mm/(s·MPan),壓強指數(shù)n=0.440 5,推進劑密度ρp=1564 kg/m3,噴管喉襯燒蝕率系數(shù)as=0.113 1 mm/(s·MPa)。由于辨識計算過程的主要依據(jù)為發(fā)動機穩(wěn)定工作段壓強數(shù)據(jù),故所得辨識參數(shù)結(jié)果具有較好的合理性。把反算的燃面數(shù)據(jù)以及修正后的參數(shù)代入內(nèi)彈道正向預(yù)示模型進行計算,可獲得如圖11所示的燃燒室壓強曲線。由圖11可知,第二次性能預(yù)示精度較第一次更高,且預(yù)示結(jié)果與實測數(shù)據(jù)高度吻合。

圖11 長通道發(fā)動機第二次辨識預(yù)示曲線Fig.11 Second prediction of the long passage SRM

4 結(jié)論

(1)提出了固體火箭發(fā)動機內(nèi)彈道多參數(shù)辨識及精準預(yù)示方法,并以C1xb型發(fā)動機和某型長通道發(fā)動機為研究案例進行計算分析,基于二次參數(shù)辨識修正和燃面反算,獲得了與實測壓強數(shù)據(jù)高度吻合的預(yù)示結(jié)果,驗證了該方法對提高內(nèi)彈道性能預(yù)示精度的有效性。

(2)除本文案例參數(shù)外,還可針對發(fā)動機使用環(huán)境剖面下具有較強散布特性的多項預(yù)示參數(shù)開展辨識修正,相關(guān)方法也適用于具有復(fù)雜裝藥構(gòu)型和燃面退移規(guī)律的固體發(fā)動機,其應(yīng)用具有較強的普遍性。

(3)基于該內(nèi)彈道參數(shù)辨識與預(yù)示方法,可充分利用有限次發(fā)動機試車實測數(shù)據(jù)開展內(nèi)彈道性能精準預(yù)示,掌握復(fù)雜環(huán)境剖面下的內(nèi)彈道性能散布情況,提高發(fā)動機設(shè)計效率并降低成本。