張 皓 師 鵬 趙育善 李保軍

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

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

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

1 基本方程

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

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

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

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

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

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

上述的兩個積分可以表示為

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

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

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

2 運動有界性分析

如果要利用徑向推力逃逸中心體的中心引力場，需對運動的有界性進行研究.

2.1 運動的邊界以及臨界條件

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 邊界上的運動分析

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

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

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

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

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

3 運動周期性分析

3.1 運動的擬周期性

其中，dr＞0時取正號，dr＜0時取負(fù)號.從式(18)容易分析運動時間和極角的變化情況.

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

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

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

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

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

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

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

此時有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性，即運動軌跡關(guān)于近地點方向是對稱的，同樣可知其關(guān)于遠(yuǎn)地點方向也是對稱的.

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

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

圖3 運動軌跡

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

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

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

3.2 周期軌道

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

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

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

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

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

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

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

圖5 迭代出界示意圖

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

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

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

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

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

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

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

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

圖6 多周期軌道仿真算例

圖7 單周期軌道仿真算例

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

4 結(jié)束語

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

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］.宇航學(xué)報，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］張皓，李保軍，趙育善.常值徑向推力下飛行器運動的可行域研究［C］//袁建平.空間操作自主控制專題研討會論文集.長沙:國防科技大學(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)霞.空間操作與非開普勒運動［J］.宇航學(xué)報，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