馮偉業(yè),陳林泉,吳 秋,陳林君
(西安航天動力技術(shù)研究所,陜西 西安 710025)
固體火箭的姿態(tài)和軌道控制很大程度上依賴于內(nèi)彈道性能的預(yù)示精度,提高內(nèi)彈道的性能預(yù)示精度對固體火箭發(fā)動機(jī)設(shè)計(jì)和導(dǎo)彈制導(dǎo)律設(shè)計(jì)有一定意義。固體火箭發(fā)動機(jī)內(nèi)彈道學(xué)的基本任務(wù)是在發(fā)動機(jī)各種工作條件下計(jì)算燃燒室工作壓強(qiáng)隨時間和空間的變化規(guī)律[1-3]。為了提高內(nèi)彈道性能預(yù)示的精度,許多學(xué)者在內(nèi)彈道方向上做了深入研究。
王磊等[4]利用Pro/E軟件三維建模,以數(shù)字化設(shè)計(jì)模擬發(fā)動機(jī)燃燒,得到精準(zhǔn)燃面變化。田忠亮等[5]通過分析橫向過載帶來的燃燒偏心對發(fā)動機(jī)燃面影響,預(yù)示了橫向過載下固體火箭發(fā)動機(jī)的內(nèi)彈道性能。劉赟等[6]通過仿真與試驗(yàn)相結(jié)合的方法對點(diǎn)火藥量對內(nèi)彈道性能影響進(jìn)行了研究。張凌等[7]通過對試驗(yàn)數(shù)據(jù)進(jìn)行參數(shù)辨識,反算得到肉厚燃面曲線,又以遺傳算法對內(nèi)彈道參數(shù)進(jìn)行了二次辨識修正得到更精準(zhǔn)的預(yù)示曲線。
由于喉襯的燒蝕過程難以在試驗(yàn)中獲得準(zhǔn)確數(shù)據(jù),常以試驗(yàn)前后的喉徑計(jì)算其燒蝕率,在內(nèi)彈道性能預(yù)示中采用線性燒蝕對其進(jìn)行簡化。然而在發(fā)動機(jī)實(shí)際工作過程中,燒蝕并非是線性的[8-10]。
文中在零維內(nèi)彈道分析的基礎(chǔ)上提出了一種基于遺傳算法和LSTM神經(jīng)網(wǎng)絡(luò)的內(nèi)彈道性能預(yù)示方法,利用試驗(yàn)實(shí)測的壓強(qiáng)數(shù)據(jù),和Pro/E模型計(jì)算出的理論燃面-肉厚曲線[11],使用遺傳算法辨識出其喉徑在燒蝕過程中的變化曲線;再以此數(shù)據(jù)建立LSTM網(wǎng)絡(luò),以建立好的模型通過理論燃面和燃速參數(shù)等對后續(xù)發(fā)動機(jī)喉徑變化曲線進(jìn)行預(yù)測,并對內(nèi)彈道性能進(jìn)行預(yù)示。經(jīng)過對案例計(jì)算分析驗(yàn)證,證明該方法可以有效的對內(nèi)彈道性能進(jìn)行預(yù)示,為內(nèi)彈道性能預(yù)示提供了一種新方法。
在零維內(nèi)彈道計(jì)算中,燃面和喉徑變化曲線是兩項(xiàng)重要的輸入,直接影響著燃燒室壓力高低,在過往的內(nèi)彈道性能預(yù)示中,常以傳熱學(xué)計(jì)算喉襯的平均燒蝕率,再以此計(jì)算喉襯的線性燒蝕變化,代入到零維內(nèi)彈道計(jì)算中。但由于喉襯在發(fā)動機(jī)實(shí)際工作中是非線性燒蝕過程,就會產(chǎn)生較大誤差。為提高零維內(nèi)彈道預(yù)示的精度,文中基于神經(jīng)網(wǎng)絡(luò)建立了喉徑燒蝕預(yù)示模型,為內(nèi)彈道性能預(yù)示提供了非線性的喉徑燒蝕曲線。該神經(jīng)網(wǎng)絡(luò)以燃面、自由容積、肉厚為輸入,喉襯燒蝕曲線為輸出。為構(gòu)建神經(jīng)網(wǎng)絡(luò)模型,需要先獲得足夠的訓(xùn)練樣本,文中以發(fā)動機(jī)地試數(shù)據(jù)為基礎(chǔ),采用遺傳算法對喉徑進(jìn)行辨識,獲取喉徑變化曲線,以此作為訓(xùn)練神經(jīng)網(wǎng)絡(luò)的數(shù)據(jù)輸入。
完整的工作流程如圖1所示。
圖1 預(yù)示方法流程圖Fig.1 Predictive method flow chart
遺傳算法(GA)起源于對生物系統(tǒng)的計(jì)算機(jī)模擬研究,是一種隨即全局搜索優(yōu)化的方法,通過對待優(yōu)化參數(shù)編碼,使編碼模擬自然遺傳過程中的復(fù)制、交叉、變異等現(xiàn)象,由適應(yīng)度函數(shù)從變異中篩選出更優(yōu)的個體,使種群不斷向好的方向進(jìn)化,最終得到最優(yōu)個體以及問題優(yōu)質(zhì)解。
在固體火箭發(fā)動機(jī)參數(shù)辨識中,采用最小二乘法等傳統(tǒng)參數(shù)辨識方法,常會因?yàn)槌跏键c(diǎn)敏感、病態(tài)梯度等問題導(dǎo)致效果不佳或陷入局部收斂,遺傳算法因其固有的全局搜索性,初值穩(wěn)健性等特點(diǎn)更適合用于固體火箭發(fā)動機(jī)的參數(shù)辨識[12]。
由于喉襯的燒蝕會使喉部面積增大導(dǎo)致燃燒室壓強(qiáng)的下降和噴管膨脹面積比的減小,因而預(yù)示出噴管喉徑的非線性變化曲線可以提高預(yù)示發(fā)動機(jī)內(nèi)彈道性能準(zhǔn)確性。為建立喉徑預(yù)示神經(jīng)網(wǎng)絡(luò)模型,需先以發(fā)動機(jī)以往的地試數(shù)據(jù)為基礎(chǔ),辨識出喉徑變化曲線作為神經(jīng)網(wǎng)絡(luò)的輸入。
在喉徑辨識中,以多項(xiàng)式曲線擬合喉徑的非線性變化,可將求喉徑變化曲線轉(zhuǎn)化為求解多項(xiàng)式的各次項(xiàng)系數(shù)。通過對多項(xiàng)式各次項(xiàng)系數(shù)的二進(jìn)制編碼生成個體。
將個體代表的喉徑變化曲線與Pro/E模型計(jì)算出的理論燃面—肉厚曲線代入零維內(nèi)彈道公式:
(1)
在零維內(nèi)彈道公式迭代計(jì)算出的壓強(qiáng)曲線結(jié)果中取m個采樣點(diǎn),其適應(yīng)度函數(shù)為:
(2)
式中C為合適的較大正數(shù)。
采取適應(yīng)度比例法對每輪產(chǎn)生的變異個體進(jìn)行篩選以識別出最終的最優(yōu)個體。對案例中的每臺發(fā)動機(jī)試驗(yàn)數(shù)據(jù)進(jìn)行識別,得到如圖2所示的喉徑變化曲線。
圖2 喉徑變化曲線Fig.2 Throat diameter variation curve
圖2中擬合曲線的前端,喉徑出現(xiàn)不同程度的下降,這是由于推導(dǎo)的理論燃面與真實(shí)燃面有一定的差距,零維內(nèi)彈道計(jì)算中是忽略了燃?xì)饬鲃拥?在工作初期,由于燃?xì)馔ǖ澜孛娣e較小,燃?xì)庠谒幹ǖ乐辛魉佥^高,從而發(fā)生侵蝕燃燒[13]。案例中的發(fā)動機(jī)是增面燃燒的藥柱設(shè)計(jì),使得實(shí)際燃面比理論燃面稍高。因此在將理論燃面當(dāng)作真實(shí)燃面并以零維內(nèi)彈道預(yù)示殘差平方和作為目標(biāo)函數(shù)進(jìn)行辨識后,這部分誤差就轉(zhuǎn)移到了喉徑變化曲線上,使得喉徑變化與真實(shí)值產(chǎn)生了偏離。但在神經(jīng)網(wǎng)絡(luò)的預(yù)示中,預(yù)示值同樣會產(chǎn)生偏離。通過代入零維內(nèi)彈道公式對內(nèi)彈道性能進(jìn)行預(yù)示,即可將誤差補(bǔ)償,得到較為準(zhǔn)確的壓力曲線。
以上節(jié)辨識出的喉徑和發(fā)動機(jī)地試數(shù)據(jù)建立喉徑燒蝕預(yù)示模型。喉徑燒蝕曲線需要以一種帶有時序記憶功能的建模方法,文中采用長短效記憶神經(jīng)網(wǎng)絡(luò),其獨(dú)特的遺忘門結(jié)構(gòu),使得它能夠控制累計(jì)記憶的信息多少,關(guān)注短期內(nèi)的有效信息,避免長時間記憶積累帶來的長時依賴問題,避免梯度爆炸或梯度消失對預(yù)測精度的影響。
試驗(yàn)得到的壓強(qiáng)時間曲線如圖3所示,摘取工作段部分對其進(jìn)行處理,以獲得建立模型的輸入數(shù)據(jù)。
圖3 壓強(qiáng)-時間曲線Fig.3 Pressure-time curve
利用燃速壓強(qiáng)公式和試驗(yàn)測得壓力-時間曲線,將肉厚-燃面曲線轉(zhuǎn)化為肉厚時間曲線、燃面時間曲線,自由容積時間曲線,如圖4~圖6所示。
圖4 肉厚-時間曲線Fig.4 Propellant thickness-time curve
圖5 燃面-時間曲線Fig.5 Burning surface-time curve
圖6 自由容積-時間曲線Fig.6 Free volume-time curve
由于不同的輸入具有不同的量綱,數(shù)值間差距可能很大,直接輸入網(wǎng)絡(luò)會導(dǎo)致網(wǎng)絡(luò)無法全面考量每個輸入對結(jié)果的影響,需要對數(shù)據(jù)進(jìn)行規(guī)范化,采用標(biāo)準(zhǔn)差標(biāo)準(zhǔn)化方法:
(3)
選用序列到序列的LSTM網(wǎng)絡(luò),以肉厚、燃面面積、自由容積、時間序列作為輸入,喉徑變化作為輸出。建立6層神經(jīng)網(wǎng)絡(luò),除輸入輸出層外,加入LSTM層和Dropout層,以及2個常規(guī)的全連接層。LSTM層可實(shí)現(xiàn)長短期記憶功能,只對臨近當(dāng)前時刻的前序狀態(tài)進(jìn)行記憶,Dropout層則可以在訓(xùn)練時隨即屏蔽部分神經(jīng)元節(jié)點(diǎn),使得神經(jīng)網(wǎng)絡(luò)適應(yīng)性更強(qiáng),同時也可以避免出現(xiàn)過擬合現(xiàn)象。神經(jīng)網(wǎng)絡(luò)的超參數(shù)設(shè)置如表1所示。
表1 神經(jīng)網(wǎng)絡(luò)超參數(shù)設(shè)置Table 1 Neural network hyperparameter settings
對神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練,圖7為訓(xùn)練過程的均方根誤差RMSE隨訓(xùn)練輪次收斂情況。
圖7 均方根誤差隨訓(xùn)練收斂圖Fig.7 RMSE varies with training
由圖7可見,模型訓(xùn)練情況較好,預(yù)測誤差收斂較快,模型建立完成。
為了對發(fā)動機(jī)的內(nèi)彈道性能進(jìn)行預(yù)示,首先要確定其工作時間。由于待預(yù)示的發(fā)動機(jī)不能提供壓強(qiáng)時間曲線,在準(zhǔn)備輸入數(shù)據(jù)時,需要以工作段燃速恒定對肉厚燃面曲線進(jìn)行處理,預(yù)測出壓強(qiáng)時間曲線后再對其進(jìn)行修正。在預(yù)示內(nèi)彈道性能前,需先對待預(yù)示發(fā)動機(jī)的燃速進(jìn)行預(yù)測。影響燃速的因素分為兩類:1)燃燒的環(huán)境,包括壓強(qiáng)、初溫等[14-15];2)推進(jìn)劑本身性質(zhì),如藥條燃速、燃速壓強(qiáng)指數(shù)等。通過建立BP神經(jīng)網(wǎng)絡(luò),對其中的函數(shù)關(guān)系進(jìn)行擬合,可以對發(fā)動機(jī)的燃速進(jìn)行較為準(zhǔn)確的預(yù)示[16-17]。
以全尺寸發(fā)動機(jī)燃燒室壓強(qiáng)、初溫,直徑為127 mm燃速測試發(fā)動機(jī)的燃速、燃燒室壓強(qiáng)、初溫,藥條燃速和燃速壓強(qiáng)指數(shù)、溫度敏感系數(shù)作為預(yù)測全尺寸發(fā)動機(jī)的基本數(shù)據(jù)集。建立3層BP神經(jīng)網(wǎng)絡(luò)。
在預(yù)示燃燒室喉徑燒蝕時,由于尚未獲得壓強(qiáng)時間曲線,需先按恒定燃速計(jì)算出肉厚時間曲線、燃面時間曲線、自由容積時間曲線,并利用神經(jīng)網(wǎng)絡(luò)對喉徑變化進(jìn)行求解。在得到喉徑變化曲線后,根據(jù)零維內(nèi)彈道公式進(jìn)行求解,獲得壓強(qiáng)時間曲線。然后對燃速進(jìn)行修正,以此壓強(qiáng)曲線和燃速壓強(qiáng)公式再對肉厚燃面進(jìn)行處理,并重新對喉徑及壓強(qiáng)進(jìn)行預(yù)示。修正過程可進(jìn)行多次。由圖3可知,案例中的發(fā)動機(jī),工作段壓強(qiáng)變化幅度較小,燃速變化不大,接近于直線,修正1次即可。預(yù)測得到喉襯變化曲線如圖8所示。
圖8 喉徑變化預(yù)測曲線Fig.8 Throat diameter change prediction curve
由于LSTM序列預(yù)測中,每個點(diǎn)的預(yù)測需要以前序狀態(tài)的輸出作為輸入,在前幾個預(yù)測點(diǎn)因?yàn)槿鄙偾靶驙顟B(tài),容易產(chǎn)生較大的突變,需要對明顯的離群點(diǎn)進(jìn)行處理,使曲線光滑。
預(yù)測曲線中的上升段使用工程計(jì)算法,對點(diǎn)火啟動段壓強(qiáng)建立過程進(jìn)行簡化計(jì)算。假設(shè)主裝藥點(diǎn)火瞬時完成,且忽略火焰?zhèn)鞑ズ腿細(xì)饬鲃?由此可得:
(4)
式中Pig為點(diǎn)火壓強(qiáng),即主裝藥剛開始點(diǎn)火時,燃燒室內(nèi)部的由點(diǎn)火裝置燃燒產(chǎn)生的壓強(qiáng)。
預(yù)測曲線中的下降段,由于有余藥剩余,余藥的燃燒彌補(bǔ)了燃燒室內(nèi)燃?xì)庠趪姽芘蛎涍^程中溫度的下降,近似看做等溫膨脹。由此可得:
(5)
分離變量求積分得:
(6)
以零維內(nèi)彈道計(jì)算公式對壓強(qiáng)進(jìn)行預(yù)測,可得到壓強(qiáng)時間曲線。
為了驗(yàn)證該方法的內(nèi)彈道性能預(yù)示的準(zhǔn)確度,添加使用傳統(tǒng)的燃面反算方法進(jìn)行預(yù)示的對照組,得到的預(yù)示結(jié)果如圖9和表2所示。
表2 部分采樣點(diǎn)預(yù)測結(jié)果表Table 2 Prediction results of some sampling points MPa
圖9 壓強(qiáng)-時間預(yù)示結(jié)果曲線Fig.9 Prediction curve of pressure-time
由圖9可見,工作段的誤差主要來源于前期,預(yù)測壓強(qiáng)比實(shí)測壓強(qiáng)要高,這部分誤差很明顯的來自于實(shí)測壓強(qiáng)的點(diǎn)火初始壓強(qiáng)峰和回落,這是由于點(diǎn)火藥粒度太細(xì)、燃燒時間過短,未能與主裝藥的點(diǎn)火和火焰?zhèn)鞑チ己玫你暯訉?dǎo)致[6]。文中預(yù)示內(nèi)彈道性能時,主要關(guān)注工作段的壓強(qiáng)預(yù)示,簡化了點(diǎn)火過程,認(rèn)為主裝藥點(diǎn)火和火焰?zhèn)鞑ニ查g完成,導(dǎo)致難以對初始壓強(qiáng)峰進(jìn)行預(yù)示。
為評估模型預(yù)測精度,引入相對誤差作為評價函數(shù),相對誤差按M=(b/a)×100%計(jì)算,其中:a為實(shí)際值的平均值;b為預(yù)測值與實(shí)際值的均方根誤差。在工作段均勻設(shè)置采樣點(diǎn),計(jì)算預(yù)測值與實(shí)際值的均方根誤差。部分采樣點(diǎn)的預(yù)測結(jié)果如表2所示。
對包含上述采樣點(diǎn)的所有采樣點(diǎn)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析,得到使用文中方法的壓強(qiáng)預(yù)示曲線均方根誤差為0.118 MPa,相對誤差為1.35%,效果總體較好。傳統(tǒng)零維內(nèi)彈道預(yù)示得到的均方根誤差為0.133 MPa,相對誤差為1.50%。為對比預(yù)測精度,對曲線進(jìn)行進(jìn)一步處理,并得到表3結(jié)果。
表3 內(nèi)彈道計(jì)算結(jié)果對比表Table 3 Comparison of internal ballistic calculation results
基于文中方法預(yù)測的最大壓強(qiáng)的相對誤差為0.62%,對應(yīng)時間相差0.04 s。基于傳統(tǒng)反算燃面的內(nèi)彈道分析方法預(yù)測的最大壓強(qiáng)的相對誤差為0.42%,對應(yīng)時間相差0.93 s。對最大壓強(qiáng)的預(yù)示中,雖然壓強(qiáng)峰值的預(yù)示效果相差0.2%,但由于峰值壓強(qiáng)對應(yīng)時間擬合偏差較小,整體曲線更貼近于實(shí)測曲線。燃燒時間的預(yù)示上,文中方法預(yù)測相對誤差為1.09%,反算燃面方法預(yù)測的相對誤差為1.44%。綜合表2~表3可見,文中所用方法對內(nèi)彈道性能的預(yù)示精度較高。
文中提出了一種基于遺傳算法和LSTM神經(jīng)網(wǎng)絡(luò)的固體火箭發(fā)動機(jī)內(nèi)彈道性能預(yù)示方法,經(jīng)案例驗(yàn)證,該方法對內(nèi)彈道壓強(qiáng)時間曲線的預(yù)示均方根誤差為0.118 MPa,相對誤差為1.35%。對最大壓強(qiáng)的預(yù)示偏差為0.054 MPa,相對誤差為0.42%,最大壓強(qiáng)響應(yīng)時間誤差為0.2 s,燃燒時間相對誤差為1.09%,能夠較為準(zhǔn)確的預(yù)測其工作段壓強(qiáng)曲線。與傳統(tǒng)的利用試驗(yàn)數(shù)據(jù)反算燃面的內(nèi)彈道預(yù)示方法相比,精度有所提升。
該方法以LSTM神經(jīng)網(wǎng)絡(luò)進(jìn)行建模并實(shí)現(xiàn)對喉徑變化的預(yù)測,需要較多的試驗(yàn)數(shù)據(jù)對模型進(jìn)行訓(xùn)練,這是該方法的局限性。若訓(xùn)練樣本較少,則準(zhǔn)確度不佳,但若是能夠?qū)⒆銐蚨喔S富的樣本試驗(yàn)數(shù)據(jù)用于訓(xùn)練深度神經(jīng)網(wǎng)絡(luò)模型,使該模型對更多工況進(jìn)行適應(yīng),會進(jìn)一步提升其準(zhǔn)確性。