張 皓 師 鵬 趙育善 李保軍

(北京航空航天大學(xué) 宇航學(xué)院，北京100191)

常值徑向連續(xù)推力機(jī)動(dòng)是常用的軌道機(jī)動(dòng)方式，文獻(xiàn)［1－7］對(duì)這類問題進(jìn)行了研究.早在1953年，文獻(xiàn)［1］就研究了二體模型下使用常值徑向推力逃離中心引力場(chǎng)的問題.后來文獻(xiàn)［2］用橢圓積分求解了初始軌道為圓軌道的情形.1998年，文獻(xiàn)［3］引入勢(shì)能阱的概念，研究了從圓軌道逃逸的條件和不發(fā)生逃逸時(shí)飛行器所能達(dá)到的最大半徑.文獻(xiàn)［4］研究了利用常值徑向推力軌道來構(gòu)建非開普勒［8］意義下的圓軌道的方法，并與開普勒?qǐng)A軌道進(jìn)行了對(duì)比.然而對(duì)初始軌道為橢圓的研究文獻(xiàn)還很少.文獻(xiàn)［5］利用動(dòng)力系統(tǒng)理論，通過數(shù)值積分得到了一些特殊的周期軌道，并分析了其穩(wěn)定性.文獻(xiàn)［6］通過幾何作圖法，得到了橢圓軌道下逃離中心引力場(chǎng)時(shí)所需的最小推力.針對(duì)橢圓初始軌道情況的解析性質(zhì)研究，目前還很少有文獻(xiàn)提及，此類研究對(duì)將來的空間應(yīng)用具有重要的意義，例如運(yùn)動(dòng)可達(dá)區(qū)［9］的分析和深空軌道的設(shè)計(jì)等.

本文考慮初始軌道為橢圓，施加常值徑向連續(xù)推力后飛行器的運(yùn)動(dòng)軌跡，并從解析角度對(duì)這一特殊軌跡的有界性和周期性進(jìn)行研究.首先建立了運(yùn)動(dòng)微分方程，并通過首次積分，將運(yùn)動(dòng)寫成求積的形式;然后通過研究一個(gè)三次多項(xiàng)式的根，分析了運(yùn)動(dòng)的有界性;最后用橢圓積分研究了運(yùn)動(dòng)的周期性并給出了周期軌道的求法.

1 基本方程

定義航天器施加推力前所處的開普勒軌道為停泊軌道，設(shè)該停泊軌道半長(zhǎng)軸為a，偏心率為e.初始時(shí)刻t0航天器的真近點(diǎn)角為f，從此刻開始施加常值徑向推力.不考慮航天器質(zhì)量變化，設(shè)推力加速度為ar.由于僅受徑向力的作用，航天器的軌道角動(dòng)量守恒，運(yùn)動(dòng)將局限在停泊軌道所在平面內(nèi)［3］.為方便研究，在此停泊軌道平面內(nèi)建立極坐標(biāo)系，以引力中心為坐標(biāo)原點(diǎn)，以停泊軌道近地點(diǎn)方向?yàn)闃O軸方向，航天器運(yùn)動(dòng)方向?yàn)闃O角增加的方向，如圖1所示.

圖1 坐標(biāo)系定義

選取極角θ和地心距(極距)r為運(yùn)動(dòng)的廣義坐標(biāo)，建立運(yùn)動(dòng)的極坐標(biāo)方程為

式中，μ為中心體的引力常數(shù);ar沿徑向向外為正.

由于常值徑向推力可以看作引力勢(shì)的一部分，即

故該系統(tǒng)為保守系統(tǒng)，存在能量積分E.另一方面，由于原點(diǎn)在推力方向上，不改變系統(tǒng)的角動(dòng)量，因此存在角動(dòng)量積分H.

上述的兩個(gè)積分可以表示為

由于該積分為常值，可以用初始條件表示:

式中，E0=－μ/2a和H0=［μa(1－e2)］1/2分別為原停泊軌道的能量和角動(dòng)量;r0=a(1－e2)/(1+ecosf)為初始時(shí)刻的地心距.利用積分可將運(yùn)動(dòng)方程(1)化成一階微分的形式:

對(duì)式(5)分離變量后直接積分可以得到r關(guān)于t的表達(dá)式，然而在積分過程中需要用橢圓函數(shù)［10］，且涉及變量代換，形式比較復(fù)雜，為更好地得到常值徑向推力下的軌跡特性，本文避免冗長(zhǎng)的解析表達(dá)，而對(duì)其進(jìn)行定性研究.

2 運(yùn)動(dòng)有界性分析

如果要利用徑向推力逃逸中心體的中心引力場(chǎng)，需對(duì)運(yùn)動(dòng)的有界性進(jìn)行研究.

2.1 運(yùn)動(dòng)的邊界以及臨界條件

在式(5)的第1式中存在根式，因而要使其有意義，需滿足如下關(guān)系［11］:

式(6)為一元三次不等式，通過對(duì)該式的分析，可以得到運(yùn)動(dòng)的范圍.令不等式(6)對(duì)應(yīng)的一元三次多項(xiàng)式為F(r)，稱之為基本多項(xiàng)式［7］:

如果F(r)=0僅有一個(gè)實(shí)根r1，那么不等式(6)的解為r＞r1，從而易知其運(yùn)動(dòng)是無界的.對(duì)于有界運(yùn)動(dòng)，F(xiàn)(r)=0必然有多于一個(gè)實(shí)根.不失一般性，假設(shè)方程F(r)=0有3個(gè)不等實(shí)根，從小到大依次為r1，r2和r3，那么不等式(6)的解可寫成:

若r0∈［r1，r2］，則其運(yùn)動(dòng)有界，且運(yùn)動(dòng)邊界分別為r1和r2，運(yùn)動(dòng)軌跡限制在半徑分別為r1和r2的同心圓環(huán)內(nèi);否則其運(yùn)動(dòng)無界.從式(4)可以看到，當(dāng)初始a，e和f確定后，方程F(r)=0的根只取決于ar，它決定了方程實(shí)根的個(gè)數(shù)以及在數(shù)軸上的分布，從而影響運(yùn)動(dòng)的有界性.本文將有界運(yùn)動(dòng)和無界運(yùn)動(dòng)的分界推力稱為臨界推力，用arc表示.

雖然式(6)的根可以利用求根公式［10］求解，但求根公式太過復(fù)雜，不易看出規(guī)律.當(dāng)初始位置比較特殊時(shí)，式(6)可以分解因式，從而可以得到運(yùn)動(dòng)邊界以及臨界推力的簡(jiǎn)單表達(dá)式［7］.例如，當(dāng)初始位置位于停泊軌道的近地點(diǎn)時(shí)，運(yùn)動(dòng)的臨界推力為

當(dāng)ar＞arc時(shí)運(yùn)動(dòng)無界，否則運(yùn)動(dòng)有界，且其邊界由下式確定:

當(dāng)初始的位置位于停泊軌道的遠(yuǎn)地點(diǎn)，也存在類似的結(jié)論，此時(shí)arc的表達(dá)式還和e有關(guān):

同樣，當(dāng)ar＞arc時(shí)運(yùn)動(dòng)無界，否則運(yùn)動(dòng)有界，且其邊界由下式確定:

對(duì)于一般的初始位置 0＜f＜π 或 π＜f＜2π，不能得到類似的簡(jiǎn)單表達(dá)式，但仍然可以用類似的方法進(jìn)行分析.通過研究式(6)的根的情況，可以得到一般情況下arc的表達(dá)式［7］:

式中，ε是方程(14)的最小實(shí)根.

有界運(yùn)動(dòng)的運(yùn)動(dòng)邊界需求解F(r)=0，此時(shí)方程有3個(gè)實(shí)根，邊界由其最小的兩個(gè)實(shí)根決定，即運(yùn)動(dòng)范圍r∈［r1，r2］.

另外，隨著推力ar的變化，方程F(r)=0的根也會(huì)發(fā)生變化.可以證明［7］，隨著ar的增大，r1和r2增大而r3減小，且當(dāng)ar=arc時(shí)，方程的根r2和r3重合，但都大于r1，如圖2所示.本文算例都使用無量綱的引力常數(shù)μ=1，后面不再說明.圖2取a=1，e=0.8，f=170°.

圖2 方程F(r)=0的根隨ar的變化趨勢(shì)

2.2 邊界上的運(yùn)動(dòng)分析

由前面的分析可知，運(yùn)動(dòng)的有界性主要由r的可行范圍決定，下面分析r達(dá)到邊界之后的運(yùn)動(dòng)趨勢(shì).將式(5)的第1式對(duì)時(shí)間求導(dǎo)并利用F(r)=0的實(shí)根，得到r相對(duì)于時(shí)間的二階導(dǎo)數(shù):

在內(nèi)邊界處，有r=r1，代入式(15)和式(5)的第1式，得到

在外邊界處，有r=r2，代入式(15)和式(5)的第1式，得到

從式(16)和式(17)可以看出:當(dāng)運(yùn)動(dòng)到達(dá)內(nèi)邊界時(shí)，r相對(duì)于時(shí)間的一階導(dǎo)數(shù)為0，二階導(dǎo)數(shù)為正，因此其向著r增大的趨勢(shì)運(yùn)動(dòng).當(dāng)運(yùn)動(dòng)到達(dá)外邊界時(shí)，r相對(duì)于時(shí)間的一階導(dǎo)數(shù)為0，二階導(dǎo)數(shù)為負(fù)，因此向著r減小的趨勢(shì)運(yùn)動(dòng).

特殊情況，當(dāng)ar=arc時(shí)，r1＜r2=r3.由式(17)可知，當(dāng)航天器運(yùn)動(dòng)到外邊界r2時(shí)，和都為0，r不再增大或者減小，運(yùn)動(dòng)的軌跡形成一個(gè)圓形.這一圓形是典型的非開普勒軌道［8］，文獻(xiàn)［4］對(duì)停泊軌道為圓時(shí)的這類軌道進(jìn)行了分析，并且和開普勒?qǐng)A軌道進(jìn)行了對(duì)比.

3 運(yùn)動(dòng)周期性分析

3.1 運(yùn)動(dòng)的擬周期性

其中，dr＞0時(shí)取正號(hào)，dr＜0時(shí)取負(fù)號(hào).從式(18)容易分析運(yùn)動(dòng)時(shí)間和極角的變化情況.

定義徑向從某次內(nèi)邊界運(yùn)動(dòng)到下一次內(nèi)邊界的運(yùn)動(dòng)時(shí)間為徑向周期，以Tr表示.這一概念類似于橢圓開普勒軌道的軌道周期.Tr可以表示為

式(19)可以用橢圓積分［10］表示:

注意，當(dāng)初始位置固定時(shí)，Tr僅決定于ar，因此可以看成ar的單變量函數(shù)，即Tr(ar).

描述軌道運(yùn)動(dòng)的另外一個(gè)重要參數(shù)是θ.由前可知，運(yùn)動(dòng)在內(nèi)外邊界循環(huán)一次，其時(shí)間為一個(gè)徑向周期，由式(18)的第2式可知這一過程中極角的變化Δθr(稱為極角轉(zhuǎn)動(dòng))為

式(21)同樣可以用橢圓積分表示:

這里的特征數(shù)n＜0，為利用文獻(xiàn)［10］給出的計(jì)算方法，需要將式(22)的橢圓積分轉(zhuǎn)化為標(biāo)準(zhǔn)形式［10］(n位于0和1 之間)，即

式中，系數(shù)c1，c2和新特征數(shù)N分別為

此時(shí)有k2＜N＜1，滿足第3類完全橢圓積分標(biāo)準(zhǔn)形的要求.

由式(22)可以看到，每完成一次徑向循環(huán)航天器轉(zhuǎn)過的極角Δθr也是相同的.當(dāng)初始位置確定，Δθr僅決定于ar，因此可以看成ar的單變量函數(shù)，即 Δθr(ar).

考慮某次內(nèi)邊界r1前后到任意極距r的極角變化，分別用Δθf和Δθb表示:

容易看到

因此軌跡具有徑向?qū)ΨQ性，即運(yùn)動(dòng)軌跡關(guān)于近地點(diǎn)方向是對(duì)稱的，同樣可知其關(guān)于遠(yuǎn)地點(diǎn)方向也是對(duì)稱的.

航天器的整體運(yùn)動(dòng)可以看成徑向運(yùn)動(dòng)和極角轉(zhuǎn)動(dòng)的合成，表現(xiàn)為一個(gè)進(jìn)動(dòng)的橢圓.橢圓的近地點(diǎn)和遠(yuǎn)地點(diǎn)各自在半徑為r1和r2的圓上漂移.如果一個(gè)徑向周期內(nèi)，航天器的極角轉(zhuǎn)動(dòng)滿足:

式中，p，q分別是互質(zhì)的整數(shù)，那么其運(yùn)動(dòng)是周期性的，否則其運(yùn)動(dòng)為擬周期的，其軌跡將會(huì)布滿整個(gè)運(yùn)動(dòng)可行區(qū)域.這類似于平面上的振子，當(dāng)各個(gè)方向振動(dòng)頻率不同時(shí)，其運(yùn)動(dòng)的軌道會(huì)布滿整個(gè)相平面［12］.圖3給出了一個(gè)運(yùn)動(dòng)的實(shí)例，參數(shù)分別為 a=1，e=0.5，f=120°，ar=0.08，從中可以看出運(yùn)動(dòng)的徑向?qū)ΨQ性和擬周期性.

圖3 運(yùn)動(dòng)軌跡

另一方面，Tr和Δθr都可以看成ar為自變量的單變量函數(shù).計(jì)算可知，隨著ar的增加，Tr和Δθr都增大.當(dāng)ar=0時(shí)，運(yùn)動(dòng)軌跡就是停泊軌道，此時(shí)Tr為停泊軌道周期T0，Δθr=0;當(dāng)ar=arc時(shí)，航天器一旦到達(dá)外邊界，之后一直沿著該邊界運(yùn)動(dòng)，因此可看作Tr和Δθr為無窮大.因此對(duì)于有界運(yùn)動(dòng)，Tr和Δθr的取值范圍是

取 a=1，e=0.5，f=120°，Tr和 Δθr隨著 ar的變化如圖4所示.

圖4 Tr和Δθr隨推力大小的變化

3.2 周期軌道

從3.1節(jié)的分析可知，當(dāng)Δθr滿足:

運(yùn)動(dòng)具有周期性，且運(yùn)動(dòng)的周期為qTr.由于整體運(yùn)動(dòng)的周期僅與q有關(guān)，稱q為周期特征數(shù).通過上式求反函數(shù)，可得所需的推力ar，即

求解式(28)會(huì)遇到以下的問題，下面逐一說明.

1)雖然利用式(28)解析求解ar是可行的，然而Δθr的表達(dá)式(22)用到了橢圓積分，使得ar具體的形式過于復(fù)雜，因此本文將用迭代方法求解.

2)迭代時(shí)需知道Δθr對(duì)ar的導(dǎo)數(shù)，而ar對(duì)Δθr的部分影響是通過方程F(r)=0的根 r1，r2和r3實(shí)現(xiàn)的，因此求導(dǎo)時(shí)不可避免要應(yīng)用鏈?zhǔn)椒▌t.文中r1，r2和r3和橢圓積分的參數(shù)有關(guān).橢圓積分本質(zhì)是一個(gè)無窮限的廣義積分，對(duì)其求導(dǎo)會(huì)出現(xiàn)無界函數(shù)，在數(shù)值計(jì)算中會(huì)出現(xiàn)困難，因此本文采用差分導(dǎo)數(shù).

3)由于周期軌道是有界運(yùn)動(dòng)的特例，因此ar的求解范圍應(yīng)限制在區(qū)間［0，arc］上.迭代過程中可能會(huì)有某次迭代值超出區(qū)間［0，arc］，如圖5所示.為避免這一現(xiàn)象對(duì)于求解的不利影響，考慮使用如下的單調(diào)映射:

將求解范圍從［0，arc］變換到(－∞，+∞)，將求解變量從ar轉(zhuǎn)化為ur.通過數(shù)值驗(yàn)證發(fā)現(xiàn)，采用這樣的代換還可以加快收斂速度.

圖5 迭代出界示意圖

4)迭代初值對(duì)于算法的收斂性影響很大.經(jīng)過式(29)變換后，收斂性可以得到較好的保證.為簡(jiǎn)單起見，文中選擇ar0=arc/2，對(duì)應(yīng)于變換后的變量ur0=0.

算例:為使仿真圖形簡(jiǎn)單，本文只考慮簡(jiǎn)單分?jǐn)?shù)的情形，即p和q至少有一個(gè)為1.

1)多周期軌道.指 p=1，q＞1的情況，即 Δθr完成一次循環(huán)的同時(shí)徑向運(yùn)動(dòng)要完成q次循環(huán).這里的多周期是相對(duì)于徑向運(yùn)動(dòng)的周期而言的.

給定如下仿真參數(shù):設(shè) a=1.41，e=0.418，f=60°.當(dāng) p=1，q=3 時(shí)，ar=0.045579211004;當(dāng) p=1，q=2 時(shí)，ar=0.046 327 987 567.迭代精度取10－12，仿真中積分器選用 MATLAB?ODE45，仿真精度為 10－12.

仿真結(jié)果如圖6，每個(gè)圖都含有4種曲線，即停泊軌道、運(yùn)動(dòng)軌跡、內(nèi)邊界和外邊界.圖6a對(duì)應(yīng)p=1，q=3的情況，可以看到徑向運(yùn)動(dòng)經(jīng)過3個(gè)周期，極角轉(zhuǎn)動(dòng)完成一次循環(huán);圖6b對(duì)應(yīng)p=1，q=2的情況，即徑向運(yùn)動(dòng)經(jīng)過兩個(gè)周期，極角轉(zhuǎn)動(dòng)完成一次循環(huán).

2)單周期軌道.指p≥1，q=1的情況，即徑向完成一次循環(huán)，極角已完成了多個(gè)循環(huán).

取a=1.41，e=0.418，f=180°.當(dāng)p=1，q=1 時(shí)，ar=0.100000000057;當(dāng) p=2，q=1 時(shí)，ar=0.103953507 988.迭代精度取10－12，仿真精度為10－12.

仿真結(jié)果如圖7，圖7a對(duì)應(yīng)p=1，q=1的情況，可以看到徑向運(yùn)動(dòng)經(jīng)過一個(gè)周期，極角轉(zhuǎn)動(dòng)也剛好完成一次循環(huán);圖7b對(duì)應(yīng)p=2，q=1的情況，徑向運(yùn)動(dòng)經(jīng)過一個(gè)周期，極角完成兩次循環(huán).

圖6 多周期軌道仿真算例

圖7 單周期軌道仿真算例

可以用類似的方法構(gòu)建不同周期的軌道，這里不再贅述.

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

本文分別從有界性和周期性研究了橢圓停泊軌道上常值徑向推力下飛行器的運(yùn)動(dòng)軌跡.研究表明運(yùn)動(dòng)存在有界和無界兩種情況，其分界由臨界推力確定.對(duì)某些特殊的情況，可以寫出運(yùn)動(dòng)邊界和臨界推力的簡(jiǎn)單表達(dá)式.應(yīng)用橢圓積分研究了有界運(yùn)動(dòng)的周期性，結(jié)果表明運(yùn)動(dòng)具有擬周期性，即徑向運(yùn)動(dòng)和極角的變化都有周期性，當(dāng)這兩個(gè)周期通約時(shí)，可以得到周期軌道.文中最后給出了一種獲取周期軌道的數(shù)值算法.本文的結(jié)果可以用來進(jìn)行深空小推力軌道的初步設(shè)計(jì)等.

References)

［1］Tsien H S.Take-off from satellite orbit［J］.Journal of the American Rocket Society，1953，24(4):233－236

［2］Battin R H.An introduction to the mathematics and methods of astrodynamics［M］.Revised Edition.Reston，VA:AIAA，1999:408－415

［3］Prussing J E，Coverstone-Carroll V.Constant radial thrust acceleration redux［J］.Journal of Guidance Control and Dynamics，1998，21(3):516－518

［4］周姜濱，袁建平，羅建軍.空間飛行器連續(xù)徑向推力機(jī)動(dòng)軌道研究［J］.宇航學(xué)報(bào)，2009，30(1):67－71

Zhou Jiangbin，Yuan Jianping，Luo Jianjun.Spacecraft orbits under continuous radial thrust［J］.Journal of Astronautics，2009，30(1):67－71(in Chinese)

［5］Akella M R，Broucke R A.Anatomy of the constant radial thrust problem［J］.Journal of Guidance Control and Dynamics，2002，26(3):563－570

［6］Mengali G，Quarta A.Escape from elliptic orbit using constant radial thrust［J］.Journal of Guidance Control and Dynamics，2009，32(3):1018－1022

［7］張皓，李保軍，趙育善.常值徑向推力下飛行器運(yùn)動(dòng)的可行域研究［C］//袁建平.空間操作自主控制專題研討會(huì)論文集.長(zhǎng)沙:國(guó)防科技大學(xué)，2009:509－516

Zhang Hao，Li Baojun，Zhao Yushan.Feasible region of spacecraft with a constant radial thrust［C］//Yuan Jianping.Seminar on Space Operation and Autonomous Control.Changsha:National University of Defense Technology，2009:509－516(in Chinese)

［8］袁建平，朱戰(zhàn)霞.空間操作與非開普勒運(yùn)動(dòng)［J］.宇航學(xué)報(bào)，2009，30(1):42－46

Yuan Jianping，Zhu Zhanxia.Space operations and non-keplerian orbit motion［J］.Journal of Astronautics，2009，30(1):42－46(in Chinese)

［9］Xue Dan，Li Junfeng，Baoyin Hexi，et al.Reachable domain for spacecraft with a single impulse［J］.Journal of Guidance Control and Dynamics，2010，33(3):934－942

［10］Abramowits M，Stegun I A.Handbook of mathematical functions with formulas，graphs，and mathematical tables［M］.Ninth Printing.New York:Dover，1972:17，567－684

［11］Vinti J P.Orbital and celestial mechanics［M］.Reston，VA:AIAA，1998:58

［12］Hirsch M W，Smale S，Devaney R L.Differential equations，dynamical systems，and an introduction to chaos［M］.Second Edition.San Diego，CA:Academic Press，2004:114－119