基于線性引力場(chǎng)的運(yùn)載火箭彈道設(shè)計(jì)研究

胡冬生,童科偉,張 烽,劉丙利,李 爍

(中國(guó)運(yùn)載火箭技術(shù)研究院 研究發(fā)展部,北京 100076)

運(yùn)載火箭彈道設(shè)計(jì)是火箭總體設(shè)計(jì)的基礎(chǔ)工作之一,對(duì)火箭總體方案、運(yùn)載性能和任務(wù)設(shè)計(jì)等都有著重要作用。在總體設(shè)計(jì)中,由于各種工程因素的制約,火箭發(fā)動(dòng)機(jī)推力并不總是按彈道計(jì)算的理論最優(yōu)推力來(lái)進(jìn)行配置,某些運(yùn)載火箭末級(jí)推重比偏小,飛行時(shí)間長(zhǎng),如阿里安5火箭和德爾它4火箭的末級(jí)起飛推重比約為0.23～0.25,連續(xù)動(dòng)力飛行時(shí)間在850 s以上。而末級(jí)點(diǎn)火時(shí)火箭剛脫離稠密大氣層,還需要繼續(xù)爬升和加速,這時(shí)規(guī)劃和設(shè)計(jì)一條最優(yōu)彈道就顯得尤為重要,這直接決定了火箭的性能或者能否將衛(wèi)星等有效載荷成功送入目標(biāo)軌道。

國(guó)內(nèi)工程單位在進(jìn)行彈道設(shè)計(jì)時(shí),在大氣層內(nèi)飛行段普遍是將飛行攻角表示成隨飛行時(shí)間指數(shù)變化的參數(shù)化經(jīng)驗(yàn)公式,在真空飛行段通常是基于平行常值引力場(chǎng)假設(shè),利用變分法或最優(yōu)控制理論推導(dǎo)得到推力方向的雙線性正切解析表達(dá)式,并進(jìn)一步簡(jiǎn)化為線性變化規(guī)律,進(jìn)而轉(zhuǎn)化為對(duì)線性斜率值和發(fā)動(dòng)機(jī)關(guān)機(jī)時(shí)間等少量參數(shù)的求解問(wèn)題,再應(yīng)用牛頓迭代法或優(yōu)化算法進(jìn)行求解,是一種工程化處理的次優(yōu)解。這種方法極大地滿足了工程需要,支撐了火箭型號(hào)方案論證,但在火箭推重比較小、飛行時(shí)間和飛行航程較長(zhǎng)的情況下,問(wèn)題求解對(duì)初值更加敏感,增加了收斂難度,且更加偏離平面地球常值引力場(chǎng)的假設(shè),因而也更加偏離最優(yōu)彈道。國(guó)內(nèi)高校在運(yùn)載火箭彈道設(shè)計(jì)方面也做了大量的研究,主要集中于智能優(yōu)化算法以及直接法中的偽譜法等。

國(guó)外在運(yùn)載火箭彈道設(shè)計(jì)方面基于打靶法、配點(diǎn)法、偽譜法等直接法開發(fā)了若干款成熟的軟件工具,在應(yīng)用間接法進(jìn)行彈道設(shè)計(jì)時(shí)則一般是與制導(dǎo)算法一體開展的,采用更加復(fù)雜的引力場(chǎng)模型和最優(yōu)化理論來(lái)進(jìn)行飛行程序角設(shè)計(jì)。DUKEMAN、BEREND等在進(jìn)行大氣層外彈道設(shè)計(jì)時(shí)采用了基于極大值原理的間接方法,將大氣層外制導(dǎo)代碼嵌入到彈道設(shè)計(jì)中,大氣層內(nèi)則采用包含2個(gè)設(shè)計(jì)變量的參數(shù)化設(shè)計(jì)方法,并應(yīng)用到戰(zhàn)神I和阿里安5火箭的彈道設(shè)計(jì)中;文獻(xiàn)[18-19]基于線性引力場(chǎng)開展真空段彈道和制導(dǎo)設(shè)計(jì),通過(guò)迭代求解含7個(gè)變量的兩點(diǎn)邊值問(wèn)題,可以滿足4～5個(gè)終端軌道參數(shù)要求。這些研究和算法對(duì)于不同推重比情況下的飛行任務(wù)具有較好的適應(yīng)能力,但由于需要同時(shí)對(duì)俯仰角和偏航角進(jìn)行優(yōu)化,當(dāng)射面與目標(biāo)軌道面有一定偏差時(shí)會(huì)產(chǎn)生較大的偏航角,從而影響火箭運(yùn)載能力,另外也未對(duì)迭代初值的設(shè)置規(guī)律進(jìn)行分析和研究。

此外,迭代制導(dǎo)也是一種通過(guò)最優(yōu)化原理推導(dǎo)出來(lái)的程序角生成方法,長(zhǎng)期應(yīng)用于國(guó)外運(yùn)載火箭飛行中,在我國(guó)火箭上也已成功應(yīng)用。但由于需要給出具體入軌位置等終端條件,且采用了基于當(dāng)前位置和入軌位置的平均常值引力場(chǎng)假設(shè),因此主要用于火箭在線制導(dǎo)和故障情況下彈道重構(gòu),無(wú)法直接應(yīng)用于彈道設(shè)計(jì),在小推重比情況下最優(yōu)性也欠佳。

綜上可知,引力場(chǎng)模型對(duì)彈道設(shè)計(jì)方法和結(jié)果有著重要的影響,甚至關(guān)系到小推重比情況下彈道優(yōu)化問(wèn)題的求解和性能。本文從線性引力場(chǎng)假設(shè)入手,通過(guò)簡(jiǎn)化偏航程序角來(lái)建立真空段彈道設(shè)計(jì)模型,將其轉(zhuǎn)換為可以迭代求解的兩點(diǎn)邊值問(wèn)題,從而為運(yùn)載火箭彈道設(shè)計(jì)提供一種新的解決思路,并優(yōu)化小推重比情況下的彈道性能。

1 彈道設(shè)計(jì)問(wèn)題的描述

1.1 運(yùn)動(dòng)方程

火箭在真空段飛行時(shí),所受的力主要是地球引力和推力,其在發(fā)射慣性坐標(biāo)系的運(yùn)動(dòng)方程為

(1)

式中:為火箭相對(duì)地心的位置矢量,為速度矢量;()為引力加速度矢量;為常值推力;為發(fā)動(dòng)機(jī)推力方向的單位矢量,假定推力方向與火箭縱軸方向重合;()為火箭當(dāng)前質(zhì)量。

1.2 線性引力場(chǎng)模型

引力加速度由位置矢量決定,引力加速度的大小僅僅取決于地心距的大小。一般情況下,地心距在火箭真空段飛行過(guò)程中變化不大,因此在設(shè)計(jì)中,引力方向的重要性要遠(yuǎn)大于引力大小的重要性。為了克服傳統(tǒng)彈道設(shè)計(jì)中平面地球常值引力場(chǎng)假設(shè)帶來(lái)的不足,本文基于二體模型(不考慮地球引力攝動(dòng)),采用線性引力場(chǎng)假設(shè)來(lái)推導(dǎo)和求解,認(rèn)為引力加速度是位置矢量的線性函數(shù):

(2)

式中:為地球引力常數(shù),為當(dāng)前時(shí)刻的地心距。在每個(gè)程序角計(jì)算周期中都進(jìn)行更新,這樣隨著火箭不斷地飛向入軌點(diǎn),該假設(shè)的誤差會(huì)越來(lái)越小。也可以借鑒迭代制導(dǎo)的計(jì)算模型,采用當(dāng)前地心距與入軌點(diǎn)地心距的平均值來(lái)替代。本文采用當(dāng)前時(shí)刻的地心距。

1.3 歸一化處理

(3)

1.4 邊界條件及性能指標(biāo)

(4)

在常值推力的情況下,真空段最優(yōu)飛行彈道為推進(jìn)劑消耗最少的飛行軌跡,即要求飛行時(shí)間最短,其性能指標(biāo)設(shè)置如下:

min=

(5)

式中:為歸一化的火箭動(dòng)力飛行時(shí)間。

2 彈道設(shè)計(jì)問(wèn)題的求解

2.1 彈道設(shè)計(jì)問(wèn)題的一般解

應(yīng)用最優(yōu)控制理論,哈密頓函數(shù)為

(6)

式中:,為協(xié)態(tài)變量;為拉格朗日乘子。

根據(jù)龐特里亞金極小值原理,最優(yōu)推力方向滿足:

(7)

協(xié)態(tài)變量滿足:

(8)

協(xié)態(tài)變量和狀態(tài)變量的解析解分別為

(9)

(10)

由協(xié)態(tài)變量和哈密頓函數(shù)的橫截條件,最終可得:

(11)

(12)

式中:,f,,f為協(xié)態(tài)變量的終端值。

2.2 程序角和協(xié)態(tài)變量的處理

(13)

對(duì)于彈道設(shè)計(jì)來(lái)說(shuō),一般設(shè)計(jì)偏航程序角=0,軌道傾角約束通過(guò)對(duì)發(fā)射方位角的迭代來(lái)滿足,這是與制導(dǎo)設(shè)計(jì)的重要區(qū)別。由軌道力學(xué)可知,飛行器速度越小,改變速度方向和射面所付出的代價(jià)也越小,因此對(duì)偏航程序角的這種處理可以減小由于真空段偏航機(jī)動(dòng)帶來(lái)的運(yùn)載能力損失。

故式(13)變?yōu)?/p>

(14)

()≡0

(15)

結(jié)合式(15)由協(xié)態(tài)方程(8)可知:

()≡0

(16)

即矢量()和()的方向分量均為0,相當(dāng)于減少了2個(gè)待定的值。由此,5個(gè)未知量,0=(,0,00),，0=(,0,00),與式(4)、式(11)、式(12)中的5個(gè)終端約束形成了一個(gè)非線性方程組,可通過(guò)迭代求解。

同時(shí),由式(14)可得,俯仰程序角公式變?yōu)?/p>

(17)

因此,每一次迭代均能給出當(dāng)前時(shí)刻的最優(yōu)俯仰程序角和火箭飛行時(shí)間,在火箭飛行過(guò)程中不斷迭代,就可計(jì)算出滿足入軌約束要求的最優(yōu)彈道。

2.3 迭代算法

本文采用Levenberg-Marquardt算法對(duì)該非線性方程組進(jìn)行求解。Levenberg-Marquardt算法是使用最廣泛的非線性最小二乘算法,其通過(guò)調(diào)節(jié)阻尼因子在高斯-牛頓法和最速下降法之間實(shí)現(xiàn)某種插值,兼具2種方法的優(yōu)點(diǎn),對(duì)迭代初值不敏感且求解速度快,因而在最優(yōu)設(shè)計(jì)、管理優(yōu)化等領(lǐng)域得到了廣泛應(yīng)用。

2.4 迭代初值的設(shè)置

在求解非線性方程組時(shí),需要設(shè)置迭代初值,進(jìn)而通過(guò)迭代算法使方程組的解逐漸收斂。迭代初值的設(shè)置對(duì)于方程組的求解具有重要影響。

運(yùn)載火箭一般為垂直發(fā)射,真空飛行段開始時(shí)火箭距離發(fā)射點(diǎn)不遠(yuǎn),俯仰程序角通常為正值,由式(7)和式(14)可知,初始的和應(yīng)為負(fù)值;但隨著火箭持續(xù)飛行,飛行速度方向逐漸接近水平,同時(shí)受到地球曲率的影響,到后半程一般會(huì)出現(xiàn)俯仰程序角為負(fù)值的情況,此時(shí)<0,>0。因此可以認(rèn)為是遞增的,由式(8)可以得到應(yīng)小于0;而由于一直小于0,由式(8)可得呈遞減變化,但的符號(hào)無(wú)法確定。據(jù)此,可設(shè)置(,0,0,0,0)的迭代初值為(-05 -05 -05 -05)或(05 -05 -05 -05)。

火箭飛行時(shí)間的迭代初值則可采用齊奧爾科夫斯基公式進(jìn)行估算。若初始飛行速度為,則由13節(jié)邊界條件可知火箭飛行入軌所需的理想速度增量為

(18)

由齊奧爾科夫斯基公式可推導(dǎo)得到對(duì)應(yīng)的推進(jìn)劑消耗量為

(19)

式中:為初始火箭質(zhì)量,為發(fā)動(dòng)機(jī)比沖。

(20)

(21)

由12節(jié)中時(shí)間的歸一化系數(shù)可將火箭飛行時(shí)間的初值設(shè)置為

(22)

需要說(shuō)明的是,由于火箭真空段飛行過(guò)程中重力和攻角會(huì)引起速度損失,根據(jù)式(22)計(jì)算出的迭代初值一般偏小,但在量級(jí)上不會(huì)存在差異,足以適應(yīng)迭代算法的收斂。

3 基于線性引力場(chǎng)的彈道設(shè)計(jì)思路

3.1 真空段彈道設(shè)計(jì)思路

根據(jù)前文推導(dǎo),對(duì)基于線性引力場(chǎng)的真空段彈道設(shè)計(jì)思路進(jìn)行整理如下:每一個(gè)程序角計(jì)算周期均根據(jù)當(dāng)前位置更新,并設(shè)置(,0,0,0,0)的迭代初值,對(duì)非線性方程組進(jìn)行求解,將求解出的程序角代入運(yùn)動(dòng)方程(1)進(jìn)行積分運(yùn)算(此時(shí)引力采用隨位置矢量變化和考慮項(xiàng)攝動(dòng)的精確模型,為地球引力二階帶諧系數(shù)),進(jìn)而開展下一個(gè)周期的計(jì)算,不斷循環(huán)?；鸺l(fā)動(dòng)機(jī)按軌道半長(zhǎng)軸關(guān)機(jī),以保證軌道半長(zhǎng)軸的設(shè)計(jì)精度。此外,第一個(gè)計(jì)算周期設(shè)置的初值至關(guān)重要,若能夠迭代成功,則將該求解出的值設(shè)為下一個(gè)計(jì)算周期迭代的初值,即能保證后續(xù)飛行時(shí)間內(nèi)均能迭代成功。

具體流程見圖1。圖中,為軌道半長(zhǎng)軸約束,為當(dāng)前時(shí)刻對(duì)應(yīng)的半長(zhǎng)軸,為一小量,根據(jù)設(shè)計(jì)精度要求進(jìn)行設(shè)定。因此,通過(guò)該流程計(jì)算出的真空段飛行彈道是自然滿足半長(zhǎng)軸、偏心率等軌道面內(nèi)的入軌要求的。

圖1 真空段彈道設(shè)計(jì)流程

3.2 全程飛行彈道設(shè)計(jì)思路

運(yùn)載火箭上升段全程飛行彈道設(shè)計(jì)一般包括一級(jí)大氣層內(nèi)飛行段和二級(jí)真空飛行段。大氣層內(nèi)飛行段仍然按照攻角的參數(shù)化經(jīng)驗(yàn)公式進(jìn)行設(shè)計(jì),設(shè)計(jì)變量為最大攻角值。由于真空段彈道僅通過(guò)迭代和積分給出,因而對(duì)于全程飛行彈道來(lái)說(shuō),設(shè)計(jì)變量只有和發(fā)射方位角,其中,主要影響射面內(nèi)的飛行彈道,通過(guò)調(diào)節(jié)來(lái)滿足一子級(jí)落點(diǎn)約束或使得火箭入軌消耗推進(jìn)劑最少;主要影響射面的方向,用于滿足軌道傾角要求?？赏ㄟ^(guò)內(nèi)外2層迭代或優(yōu)化進(jìn)行求解,其中內(nèi)層用于求解,外層用于求解。

若對(duì)一子級(jí)落點(diǎn)無(wú)約束(如海上發(fā)射),設(shè)計(jì)目標(biāo)是火箭消耗推進(jìn)劑最少,則全程飛行彈道設(shè)計(jì)具體流程見圖2。圖中,為軌道傾角約束;為當(dāng)前時(shí)刻對(duì)應(yīng)的軌道傾角;為一小量,根據(jù)設(shè)計(jì)精度要求進(jìn)行設(shè)定;的初值可由軌道傾角約束、發(fā)射點(diǎn)緯度通過(guò)球面幾何公式得出。若對(duì)一子級(jí)落點(diǎn)航程有約束,則僅需將圖中的內(nèi)層優(yōu)化改為通過(guò)大氣層內(nèi)彈道積分來(lái)迭代求解,真空段彈道設(shè)計(jì)不參與內(nèi)層迭代。

圖2 全程飛行彈道設(shè)計(jì)流程

在真空段采用線性引力場(chǎng)進(jìn)行彈道計(jì)算的起始點(diǎn),迭代得到的俯仰程序角相對(duì)于大氣層內(nèi)終端程序角可能會(huì)有一定差異甚至跳躍。為避免這種跳躍對(duì)火箭姿態(tài)控制帶來(lái)的不利影響,應(yīng)對(duì)迭代計(jì)算獲得的程序角速率加以限制,本文在計(jì)算中限制程序角速率不超過(guò)2(°)/s。

4 仿真分析

本節(jié)以美國(guó)SpaceX公司獵鷹9火箭為仿真對(duì)象,驗(yàn)證基于線性引力場(chǎng)的彈道設(shè)計(jì)方法的正確性,以及其在小推重比情況下的優(yōu)勢(shì),并給出該設(shè)計(jì)思路在運(yùn)載火箭上升段全程彈道設(shè)計(jì)中的應(yīng)用。

4.1 真空段彈道設(shè)計(jì)

以獵鷹9火箭二級(jí)飛行進(jìn)入200 km近地軌道為例進(jìn)行仿真。相關(guān)設(shè)計(jì)輸入見表1。為驗(yàn)證本文所提算法的魯棒性,在4.1節(jié)和4.2節(jié)的所有算例中,均按2.4節(jié)的方法來(lái)設(shè)置第一個(gè)計(jì)算周期中(,0,0,0,0)的迭代初值,結(jié)果表明均能在8次左右收斂,后續(xù)的迭代則一般在5次以內(nèi)收斂。

表1 真空段彈道設(shè)計(jì)輸入初始參數(shù)

在標(biāo)準(zhǔn)推力下(初始推重比0.716),分別采用基于線性變化率的傳統(tǒng)設(shè)計(jì)方法、基于平均常值引力模型的迭代制導(dǎo)方法和基于線性引力場(chǎng)的彈道設(shè)計(jì)方法進(jìn)行仿真計(jì)算,結(jié)果見表2及圖3、圖4。由表2可知,3種方法得到的飛行時(shí)間和對(duì)應(yīng)的推進(jìn)劑消耗量均差異很小,速度損失接近。由圖3和圖4可見,3種方法下的飛行地心距和俯仰程序角均基本重合,且俯仰程序角呈典型的近線性變化,驗(yàn)證了本文所推導(dǎo)方法的正確性。

表2 真空段彈道設(shè)計(jì)結(jié)果(標(biāo)準(zhǔn)推力)

圖3 火箭二級(jí)飛行地心距

圖4 火箭二級(jí)飛行俯仰程序角

當(dāng)推力降低至53%時(shí)(初始推重比0.379),采用3種方法得到的設(shè)計(jì)結(jié)果見表3及圖5、圖6,飛行彈道和俯仰程序角的整體變化趨勢(shì)較為一致。由表3可知,基于線性引力場(chǎng)的設(shè)計(jì)結(jié)果要比其他2種方法的結(jié)果分別少消耗推進(jìn)劑3 525.7 kg、2 509.5 kg,相當(dāng)于節(jié)省推進(jìn)劑2.9%、2.1%,速度損失分別少1 054.345 m/s、715.317 m/s,意味著飛行彈道更優(yōu),體現(xiàn)出了本文所用方法的優(yōu)勢(shì),而迭代制導(dǎo)仿真相比于傳統(tǒng)設(shè)計(jì)結(jié)果的優(yōu)勢(shì)并不明顯。

表3 真空段彈道設(shè)計(jì)結(jié)果(推力降低)

圖5 火箭二級(jí)飛行地心距(推力降低)

圖6 火箭二級(jí)飛行俯仰程序角(推力降低)

此外,在小推重比情況下采用單個(gè)線性段的傳統(tǒng)設(shè)計(jì)方法已不能迭代出飛行彈道,需要分2個(gè)乃至更多個(gè)線性段進(jìn)行處理,見圖5、圖6。

由圖5、圖6中的線性引力場(chǎng)設(shè)計(jì)結(jié)果可見,在二級(jí)飛行的前500 s時(shí)間里,火箭以50°左右的較大程序角飛行且近似呈線性緩慢減小,直到關(guān)機(jī)前300多秒才逐漸快速地減小至-60°左右,因而,在小推重比情況下飛行程序角已不適宜用簡(jiǎn)單的線性變化來(lái)近似。這也使得整個(gè)彈道的飛行高度呈現(xiàn)跳躍、起伏的情況,中間高度最小約為133 km。當(dāng)推力進(jìn)一步降低時(shí),有可能出現(xiàn)中間的飛行高度回落到大氣層內(nèi)的情況,這在設(shè)計(jì)中是應(yīng)該避免的,這時(shí)需要犧牲一定的優(yōu)化性能來(lái)優(yōu)先保證火箭不再入大氣層,是未來(lái)應(yīng)予以研究的一個(gè)方向。而當(dāng)推力降低到一定程度時(shí),火箭將沒(méi)有足夠大的推力加速度來(lái)擺脫地球引力,將再入大氣燒毀,無(wú)法完成入軌任務(wù)。

4.2 全程飛行彈道設(shè)計(jì)

同樣以獵鷹9火箭在卡納維拉爾角執(zhí)行28.5°傾角、200 km近地軌道的發(fā)射任務(wù)為例進(jìn)行全程彈道的仿真,二級(jí)推力采用標(biāo)準(zhǔn)值。設(shè)定一子級(jí)落點(diǎn)約束,經(jīng)過(guò)迭代的一級(jí)飛行段最大攻角為1.369 9°,發(fā)射方位角為92.66°,總飛行時(shí)間為525.873 9 s,入軌時(shí)刻火箭質(zhì)量為28 687.9 kg,扣除入軌時(shí)刻火箭二子級(jí)質(zhì)量后,衛(wèi)星質(zhì)量為22 745 kg,與獵鷹9火箭公布的運(yùn)載能力(22.8 t)極為接近。地心距和俯仰程序角曲線見圖7和圖8。由圖8可知,二級(jí)飛行初始俯仰程序角在短時(shí)間內(nèi)有一定幅度的增加,文中對(duì)角速率進(jìn)行了限幅,確保程序角速率不超過(guò)2(°)/s。

圖7 火箭全程飛行地心距

圖8 火箭全程飛行俯仰程序角

5 結(jié)束語(yǔ)

本文基于線性引力場(chǎng)假設(shè),通過(guò)簡(jiǎn)化偏航程序角和協(xié)態(tài)變量建立了真空段彈道優(yōu)化設(shè)計(jì)模型,將最優(yōu)彈道的計(jì)算轉(zhuǎn)換為兩點(diǎn)邊值問(wèn)題的求解,并給出了設(shè)計(jì)思路和設(shè)計(jì)流程,以及協(xié)態(tài)變量和火箭飛行時(shí)間迭代初值的設(shè)置方法。以獵鷹9火箭為算例進(jìn)行了仿真,結(jié)果表明,該方法應(yīng)用于某給定小推重比情況下彈道設(shè)計(jì)時(shí)具有更好的優(yōu)化效果,相比傳統(tǒng)設(shè)計(jì)方法速度損失減小1 054.345 m/s,節(jié)省推進(jìn)劑2.9%,相比迭代制導(dǎo)仿真速度損失減小715.317 m/s,節(jié)省推進(jìn)劑2.1%;與一級(jí)飛行彈道參數(shù)化設(shè)計(jì)相結(jié)合可實(shí)現(xiàn)上升段全程彈道設(shè)計(jì)和優(yōu)化。

此外,該方法和思路在運(yùn)載火箭在線軌跡規(guī)劃、發(fā)動(dòng)機(jī)推力下降故障情況下彈道重構(gòu)以及上面級(jí)軌跡規(guī)劃等方面均具有應(yīng)用潛力。