李 柯,聶萬勝,馮必鳴

(裝備學(xué)院,北京 101416)

助推-滑翔飛行器具有遠(yuǎn)程快速精確打擊、大范圍區(qū)域覆蓋、機(jī)動性好等優(yōu)點(diǎn),已成為國內(nèi)外研究的熱點(diǎn)。其整個飛行過程可描述為:固體火箭從地面垂直發(fā)射,經(jīng)過垂直上升和轉(zhuǎn)彎過程后發(fā)動機(jī)關(guān)機(jī),飛行器與火箭和整流罩分離,這段過程稱為主動段;飛行器與箭體分離后,依靠氣動控制實(shí)現(xiàn)遠(yuǎn)距離無動力滑翔,完成對目標(biāo)的打擊,這段過程稱為滑翔段。

整個飛行過程中,飛行器可以達(dá)到的極限距離為最大射程,是評價(jià)助推-滑翔飛行器性能的一項(xiàng)重要指標(biāo)。文獻(xiàn)[1]利用解析法近似分析了最大射程和升阻比及關(guān)機(jī)點(diǎn)速度的關(guān)系。文獻(xiàn)[2]利用Gauss偽譜法對 CAV再入滑翔段的軌跡進(jìn)行了優(yōu)化,突出了偽譜法在計(jì)算效率和精度上的優(yōu)勢。文獻(xiàn)[3]基于Gauss偽譜法對固體火箭上升段進(jìn)行了優(yōu)化研究。文獻(xiàn)[4-5]設(shè)定少量的約束條件,通過序列二次規(guī)劃法對整體過程進(jìn)行了優(yōu)化?？偟膩碇v,目前對全程彈道優(yōu)化研究的仍然很少,尤其對于主動段關(guān)機(jī)點(diǎn)指標(biāo)的選取仍不明確。本文在上述文獻(xiàn)的基礎(chǔ)上,利用Radau偽譜法分別對主動段和滑翔段進(jìn)行優(yōu)化,然后將不同優(yōu)化方案求得的最大射程進(jìn)行了對比分析,明確了主動段的性能指標(biāo),給出了最佳方案。

1 數(shù)學(xué)模型

1.1 動力學(xué)模型

在方案論證階段可忽略次要影響因素,因此不考慮地球自轉(zhuǎn)、扁率和彈體自身滾轉(zhuǎn),其簡化的主動段運(yùn)動方程為

其中,V為飛行速度, T為發(fā)動機(jī)推力,D為阻力, L為升力,α為攻角,θ為彈道傾角,s為航程,h為高度, m為火箭總質(zhì)量,R0為地球平均半徑,g為重力加速度,I為發(fā)動機(jī)比沖。飛行器處于滑翔段時(shí),自身不攜帶發(fā)動機(jī),質(zhì)量為定值,所以將T＝0,d m / d t ＝ 0 代入式(1)即得到滑翔段的運(yùn)動方程。

式中,密度ρ、重力加速度g,阻力D和升力L可由式(2)計(jì)算得出:

其中,0ρ為海平面大氣密度,H為常量,μ為引力常數(shù),CD為阻力系數(shù),CL為升力系數(shù),S為氣動參考面積。

1.2 氣動模型

固體火箭采用軸對稱結(jié)構(gòu),其氣動參數(shù)可按照文獻(xiàn)[6]選取:

飛行器采用升力體構(gòu)型,氣動外形如圖1所示。以飛行高度30km為例,壓強(qiáng)1170Pa,溫度216.5K,大氣密度 1.9×10-5kg/m2,計(jì)算得到不同馬赫數(shù)狀態(tài)下,升力系數(shù)、阻力系數(shù)和升阻比隨攻角的變化,如圖 2至圖4所示。

圖1 飛行器的氣動外形

圖2 升力系數(shù)隨攻角變化情況

圖3 阻力系數(shù)隨攻角變化情況

圖4 升阻比隨攻角變化情況

由圖2、圖3可知,升力系數(shù)隨攻角的變化曲線近似為一次函數(shù),阻力系數(shù)隨攻角的變化曲線近似為二次函數(shù)。在馬赫數(shù)高于9時(shí),氣動特性曲線非常接近,這與高馬赫數(shù)無關(guān)原理相吻合。因此,氣動系數(shù)可擬合成攻角的函數(shù),即

1.3 優(yōu)化模型

不失一般性,考慮 Bolza形式的最優(yōu)控制問題,設(shè)狀態(tài)變量 x （ τ）∈ Rn,控制變量 u （τ）∈ Rm,初始時(shí)間t0,終端時(shí)間tf,指標(biāo)函數(shù)為[7]

受材料和自身結(jié)構(gòu)限制,必須使駐點(diǎn)熱流、動壓和過載始終處于飛行器允許的范圍之內(nèi),即路徑約束:

狀態(tài)變量的初始值和終端值要根據(jù)初始條件和末端指標(biāo)來確定,即邊界條件:

其中包括:初始參數(shù)h0、0θ、s0、V0、m0,末端參數(shù)hf、fθ、sf、Vf、mf。

攻角作為控制量,應(yīng)處于飛行器所允許的最大可控范圍之內(nèi),即控制量約束:

2 Radau優(yōu)化算法

Radau偽譜法將狀態(tài)變量和控制變量在一系列LGR(Legendre-Gauss-Radau)點(diǎn)上離散,并以離散點(diǎn)為節(jié)點(diǎn)構(gòu)造Largrange插值多項(xiàng)式來擬合狀態(tài)變量和控制變量。通過對全局插值多項(xiàng)式求導(dǎo)來近似狀態(tài)變量對時(shí)間的導(dǎo)數(shù),將動力學(xué)微分方程約束轉(zhuǎn)化成一組代數(shù)約束,將軌跡優(yōu)化問題最終轉(zhuǎn)化為非線性規(guī)劃問題[8]。

由于LGR點(diǎn)所定義的時(shí)間區(qū)間 τ ＝ [-1 ,1),若要把動力學(xué)模型在LGR點(diǎn)上進(jìn)行離散,需要把優(yōu)化問題的時(shí)間區(qū)間 t ＝[t0, tf)通過式(11)轉(zhuǎn)換:

轉(zhuǎn)換后,τ取代t成為獨(dú)立變量,則指標(biāo)函數(shù),微分方程,約束條件和邊界條件分別為

狀態(tài)變量和控制變量用Lagrange插值多項(xiàng)式近似表示:

經(jīng)過上述離散化,指標(biāo)函數(shù)、微分方程、路徑約束和邊界條件分別為

其中, wk,（ k ＝ 1 ,… … ,N -1）為LGR權(quán)重。

最終將連續(xù)問題轉(zhuǎn)化為非線性規(guī)劃問題進(jìn)行求解。

3 仿真及分析

在主動段,初始條件設(shè)為h0=0km,V0=0m/s,θ0= 90°,由于整流罩的存在,可不考慮熱流密度約束,僅考慮動壓和法向過載約束,qmax＝ 1 00kPa,nymax＝ 1 .5。洲際固體導(dǎo)彈的攻角一般不超過7°,所以控制量約束設(shè)為 - 7°≤ α ≤ 7°。假設(shè)發(fā)動機(jī)關(guān)機(jī)時(shí)速度為5000m/s,以此設(shè)計(jì)出兩級固體火箭的性能參數(shù)如表1所示,具體設(shè)計(jì)方法參考文獻(xiàn)[9]。

由于發(fā)動機(jī)的關(guān)機(jī)點(diǎn)參數(shù)(即滑翔段的初始參數(shù))直接決定了滑翔彈道的形狀和射程,因此,根據(jù)不同的性能指標(biāo),主動段采取不同的方案進(jìn)行優(yōu)化,然后從中挑選出最佳方案。

方案1:主動段飛行距離最遠(yuǎn) J ＝ m axsf;

方案3:關(guān)機(jī)速度最大 J ＝ m axVf。

不同方案優(yōu)化出的關(guān)機(jī)點(diǎn)參數(shù)如表2所示。

表1 固體火箭性能參數(shù)

表2 關(guān)機(jī)點(diǎn)參數(shù)

在滑翔段,設(shè)飛行器氣動參考面積 S ＝ 0 .48387m2,質(zhì)量m ＝ 9 07.18kg,將主動段的關(guān)機(jī)點(diǎn)參數(shù)作為初始值,仿真終止條件為 hf=0km 。需要考慮路徑約束包括動壓約束、過載約束和駐點(diǎn)熱流密度約束,即控制量約 束 為0≤ α ≤ 3 0°。 性 能 指 標(biāo) 為 滑 翔 距 離 最 遠(yuǎn)

如圖5和圖6所示,3種方案的彈道形狀有很大差異。方案 1中,滑翔段的初始傾角較小,彈道平滑,跳躍幅度小,飛行高度都保持在 30～50km,該高度為大部分防空導(dǎo)彈的攔截盲區(qū),利于突防,攻角變化較為平滑,易于控制,但射程最小;方案2中,滑翔段的初始傾角較大,彈道的跳躍幅度較大,自由飛行彈道(近似于彈道導(dǎo)彈的中段飛行)較長,不利于突防,而且在200～543s這段時(shí)間內(nèi),攻角出現(xiàn)劇烈波動,不利于控制,射程比方案1略大;方案3中,彈道跳躍幅度適中,攻角最大值為 16.8o,沒有出現(xiàn)劇烈波動的情況,且射程最遠(yuǎn)。由此可見,主動段的終端傾角決定滑翔段的跳躍幅度,主動段的終端速度決定最大射程。

圖5 不同方案的彈道曲線

圖6 不同方案的控制量變化曲線

圖7 主動段路徑約束變化情況

圖8 滑翔段路徑約束變化情況

表3 路徑約束最大值

由表3可知,3種方案的路徑約束最大值都滿足本文所設(shè)的約束條件,驗(yàn)證了方案的可行性。對比方案2和方案3,在熱流密度最大值相同的情況下,方案3的射程最遠(yuǎn),此外,方案2在滑翔段的路徑約束最大值更接近極限值,這對飛行器的自身承受能力要求較高,因此方案3優(yōu)于方案2。對比方案1和方案3,方案1在滑翔段的路徑約束最大值皆小于方案3,但在主動段的法向最大過載接近極限值,這對固體火箭的結(jié)構(gòu)強(qiáng)度要求較高,同時(shí)射程也小于方案3,可見方案3優(yōu)于方案1。綜合對比發(fā)現(xiàn),方案3射程最遠(yuǎn),實(shí)現(xiàn)難度適中,可選為最佳方案。

文獻(xiàn)[1]利用解析法,得到的最大射程計(jì)算公式為

式中,L/ D為飛行器的升阻比,Vf為主動段的關(guān)機(jī)點(diǎn)速度,Vs為第一宇宙速度,可見最大射程與飛行器的升阻比成正比。令 L / D ＝ 3 .681,即最大升阻比,求得方案3的最大射程解析值如表4所示。

由表 4可知,通過本文的設(shè)計(jì)方案,優(yōu)化后最大射程比解析值增加了9.8%,優(yōu)化方案可行。

4 結(jié)束語

本文利用Radau偽譜法,將軌跡優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃問題,分別優(yōu)化了主動段和滑翔段的彈道,仿真計(jì)算了不同方案時(shí)的飛行器最大射程,得出了以下結(jié)論:

1)性能指標(biāo)選取為主動段關(guān)機(jī)點(diǎn)速度最大和滑翔段射程最遠(yuǎn)的組合方式,可以使總射程達(dá)到最大,實(shí)現(xiàn)難度適中;

2)主動段的終端傾角決定滑翔段的跳躍幅度,主動段的終端速度決定最大射程;

3)采用分段優(yōu)化的方法處理助推滑翔問題,便于分析各段的彈道特征,為技戰(zhàn)術(shù)指標(biāo)的制定提供參考。

[1]Eggers A J, Allen J H. A comparative analysis of the performance of long-range hypervelocity vehicles[R].NACA Technical Report, March l954.

[2]雍恩米.高超聲速滑翔式再入飛行器軌跡優(yōu)化與制導(dǎo)方法研究[D].長沙:國防科技大學(xué), 2008.

[3]楊希祥.基于Gauss偽譜法的固體運(yùn)載火箭上升段軌跡快速優(yōu)化研究[J].宇航學(xué)報(bào),2011,32(1):15-21.

[4]李瑜.助推-滑翔導(dǎo)彈彈道優(yōu)化研究[J].宇航學(xué)報(bào),2008,29(1):66-71.

[5]李瑜.助推-滑翔導(dǎo)彈最大射程優(yōu)化[J].彈道學(xué)報(bào),2008,20(4):53-56.

[6]何麟書.彈道導(dǎo)彈和運(yùn)載火箭設(shè)計(jì)[M].北京:北京航空航天大學(xué)出版社,2002.

[7]胡壽松.自動控制原理[M].北京:科學(xué)出版社,2004.

[8]Divya Garg, Michael A Patterson. Direct trajectory optimization and costate estimation of general optimal problems using a radau pseudospectral method[C].AIAA Guidance, Navigation, and Control Conference.AIAA, August 2009.

[9][蘇]A.M.西紐科夫.固體彈道式導(dǎo)彈[M].北京:國防工業(yè)出版社,1984.