基于線性引力場的運載火箭彈道設計研究

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

(中國運載火箭技術研究院 研究發展部,北京 100076)

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

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

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

此外,迭代制導也是一種通過最優化原理推導出來的程序角生成方法,長期應用于國外運載火箭飛行中,在我國火箭上也已成功應用。但由于需要給出具體入軌位置等終端條件,且采用了基于當前位置和入軌位置的平均常值引力場假設,因此主要用于火箭在線制導和故障情況下彈道重構,無法直接應用于彈道設計,在小推重比情況下最優性也欠佳。

綜上可知,引力場模型對彈道設計方法和結果有著重要的影響,甚至關系到小推重比情況下彈道優化問題的求解和性能。本文從線性引力場假設入手,通過簡化偏航程序角來建立真空段彈道設計模型,將其轉換為可以迭代求解的兩點邊值問題,從而為運載火箭彈道設計提供一種新的解決思路,并優化小推重比情況下的彈道性能。

1 彈道設計問題的描述

1.1 運動方程

火箭在真空段飛行時,所受的力主要是地球引力和推力,其在發射慣性坐標系的運動方程為

(1)

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

1.2 線性引力場模型

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

(2)

式中:為地球引力常數,為當前時刻的地心距。在每個程序角計算周期中都進行更新,這樣隨著火箭不斷地飛向入軌點,該假設的誤差會越來越小。也可以借鑒迭代制導的計算模型,采用當前地心距與入軌點地心距的平均值來替代。本文采用當前時刻的地心距。

1.3 歸一化處理

(3)

1.4 邊界條件及性能指標

(4)

在常值推力的情況下,真空段最優飛行彈道為推進劑消耗最少的飛行軌跡,即要求飛行時間最短,其性能指標設置如下:

min=

(5)

式中:為歸一化的火箭動力飛行時間。

2 彈道設計問題的求解

2.1 彈道設計問題的一般解

應用最優控制理論,哈密頓函數為

(6)

式中:,為協態變量;為拉格朗日乘子。

根據龐特里亞金極小值原理,最優推力方向滿足:

(7)

協態變量滿足:

(8)

協態變量和狀態變量的解析解分別為

(9)

(10)

由協態變量和哈密頓函數的橫截條件,最終可得:

(11)

(12)

式中:,f,,f為協態變量的終端值。

2.2 程序角和協態變量的處理

(13)

對于彈道設計來說,一般設計偏航程序角=0,軌道傾角約束通過對發射方位角的迭代來滿足,這是與制導設計的重要區別。由軌道力學可知,飛行器速度越小,改變速度方向和射面所付出的代價也越小,因此對偏航程序角的這種處理可以減小由于真空段偏航機動帶來的運載能力損失。

故式(13)變為

(14)

()≡0

(15)

結合式(15)由協態方程(8)可知:

()≡0

(16)

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

同時,由式(14)可得,俯仰程序角公式變為

(17)

因此,每一次迭代均能給出當前時刻的最優俯仰程序角和火箭飛行時間,在火箭飛行過程中不斷迭代,就可計算出滿足入軌約束要求的最優彈道。

2.3 迭代算法

本文采用Levenberg-Marquardt算法對該非線性方程組進行求解。Levenberg-Marquardt算法是使用最廣泛的非線性最小二乘算法,其通過調節阻尼因子在高斯-牛頓法和最速下降法之間實現某種插值,兼具2種方法的優點,對迭代初值不敏感且求解速度快,因而在最優設計、管理優化等領域得到了廣泛應用。

2.4 迭代初值的設置

在求解非線性方程組時,需要設置迭代初值,進而通過迭代算法使方程組的解逐漸收斂。迭代初值的設置對于方程組的求解具有重要影響。

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

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

(18)

由齊奧爾科夫斯基公式可推導得到對應的推進劑消耗量為

(19)

式中:為初始火箭質量,為發動機比沖。

(20)

(21)

由12節中時間的歸一化系數可將火箭飛行時間的初值設置為

(22)

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

3 基于線性引力場的彈道設計思路

3.1 真空段彈道設計思路

根據前文推導,對基于線性引力場的真空段彈道設計思路進行整理如下:每一個程序角計算周期均根據當前位置更新,并設置(,0,0,0,0)的迭代初值,對非線性方程組進行求解,將求解出的程序角代入運動方程(1)進行積分運算(此時引力采用隨位置矢量變化和考慮項攝動的精確模型,為地球引力二階帶諧系數),進而開展下一個周期的計算,不斷循環。火箭發動機按軌道半長軸關機,以保證軌道半長軸的設計精度。此外,第一個計算周期設置的初值至關重要,若能夠迭代成功,則將該求解出的值設為下一個計算周期迭代的初值,即能保證后續飛行時間內均能迭代成功。

具體流程見圖1。圖中,為軌道半長軸約束,為當前時刻對應的半長軸,為一小量,根據設計精度要求進行設定。因此,通過該流程計算出的真空段飛行彈道是自然滿足半長軸、偏心率等軌道面內的入軌要求的。

圖1 真空段彈道設計流程

3.2 全程飛行彈道設計思路

運載火箭上升段全程飛行彈道設計一般包括一級大氣層內飛行段和二級真空飛行段。大氣層內飛行段仍然按照攻角的參數化經驗公式進行設計,設計變量為最大攻角值。由于真空段彈道僅通過迭代和積分給出,因而對于全程飛行彈道來說,設計變量只有和發射方位角,其中,主要影響射面內的飛行彈道,通過調節來滿足一子級落點約束或使得火箭入軌消耗推進劑最少;主要影響射面的方向,用于滿足軌道傾角要求。可通過內外2層迭代或優化進行求解,其中內層用于求解,外層用于求解。

若對一子級落點無約束(如海上發射),設計目標是火箭消耗推進劑最少,則全程飛行彈道設計具體流程見圖2。圖中,為軌道傾角約束;為當前時刻對應的軌道傾角;為一小量,根據設計精度要求進行設定;的初值可由軌道傾角約束、發射點緯度通過球面幾何公式得出。若對一子級落點航程有約束,則僅需將圖中的內層優化改為通過大氣層內彈道積分來迭代求解,真空段彈道設計不參與內層迭代。

圖2 全程飛行彈道設計流程

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

4 仿真分析

本節以美國SpaceX公司獵鷹9火箭為仿真對象,驗證基于線性引力場的彈道設計方法的正確性,以及其在小推重比情況下的優勢,并給出該設計思路在運載火箭上升段全程彈道設計中的應用。

4.1 真空段彈道設計

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

表1 真空段彈道設計輸入初始參數

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

表2 真空段彈道設計結果(標準推力)

圖3 火箭二級飛行地心距

圖4 火箭二級飛行俯仰程序角

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

表3 真空段彈道設計結果(推力降低)

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

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

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

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

4.2 全程飛行彈道設計

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

圖7 火箭全程飛行地心距

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

5 結束語

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

此外,該方法和思路在運載火箭在線軌跡規劃、發動機推力下降故障情況下彈道重構以及上面級軌跡規劃等方面均具有應用潛力。