三維引力輔助機理分析

賈建華 王 琪

(北京航空航天大學 航空科學與工程學院，北京100191)

引力輔助技術(shù) (gravity-assist)，又稱借力飛行或引力甩擺，是一種航天任務(wù)中經(jīng)常用來減少燃料消耗的方法.目前已經(jīng)有多個任務(wù)利用引力輔助技術(shù)設(shè)計了行星際探測軌道，完成了太陽系八大行星的探測任務(wù)［1］.

引力輔助技術(shù)起源于18世紀，Laplace、leverrier、Tisserand用此來解釋行星對彗星軌道的攝動，Lawden首先提出可將引力輔助技術(shù)應(yīng)用于飛行器［2］.目前對引力輔助技術(shù)的研究主要集中在應(yīng)用方面，如引力輔助序列的選?。?－4］、小推力引力輔助的研究［5－6］、氣動引力輔助的研究［7－8］，而對引力輔助機理本身的探討和研究較少.1988年文獻［9］對二維引力輔助機理進行了分析，給出了前側(cè)飛越使飛行器能量減少、后側(cè)飛越使飛行器能量增加等重要結(jié)論，揭開了人類對引力輔助機理認識的新篇章，但討論僅限于平面圓型限制性三體問題.1997年文獻 ［10］采用平面橢圓型限制性三體模型對引力輔助機理進行了深入研究，結(jié)果表明，偏心率和真近點角對引力輔助軌道具有重要影響.文獻 ［11－12］也對基于平面限制性三體問題模型的引力輔助機理進行了一些討論.以上的研究都只限于二維引力輔助，飛行器運動被限制在借力天體的運動平面內(nèi)，這顯然不能滿足實際任務(wù)的需要.目前，對三維引力輔助機理的研究尚不多見，文獻［13］初步討論了面外角對引力輔助軌道的影響，但是，對面外角的研究較為簡略，而且對飛行器近心點速度有一定限制，因此需要對三維引力輔助機理進行進一步的研究.

為此，本文在前人研究工作的基礎(chǔ)上，通過引入兩個參數(shù)，將對引力輔助機理的探討從平面橢圓型限制性三體模型拓展到三維橢圓型限制性三體模型，重點討論了這兩個參數(shù)對引力輔助前后的飛行器軌道類型、軌道能量和軌道傾角的影響.最后，以地月系統(tǒng)為例，針對不同任務(wù)要求，給出了引力輔助參數(shù)的選擇區(qū)域，并對引力輔助參數(shù)進行了優(yōu)化.

1 問題描述

在三維橢圓型限制性三體模型下，三維引力輔助可通過6個自由參數(shù)來描述 (如圖1所示):①Vp:飛越借力天體時近心點處的速度;②rp:飛越借力天體時近心點距離;③ν:借力天體的真近點角;④ψ:rp在x-y平面的投影與M1，M2連線的夾角;⑤α:rp與x-y平面的夾角;⑥β:VP與x-y平面的夾角.

圖1 三維橢圓型限制性三體模型引力輔助示意圖

2 三維橢圓型限制性三體模型

假設(shè)兩個主天體分別為M1和M2，繞其公共質(zhì)心做橢圓運動，第三體M3(本文中為飛行器)的質(zhì)量遠遠小于兩個主天體的質(zhì)量，在兩個主天體的作用下運動，而不影響兩個主天體的運動.M3的運動軌跡不限于兩個主天體的運動平面內(nèi)，可以是三維的.

為了討論和計算上的方便，采用無量綱形式，定義:以M1和M2相對軌道的半長軸為單位長度1;平均角速度為1;單位質(zhì)量為m1+m2，其中，m1，m2分別是 M1和 M2的真實質(zhì)量，這樣，M2的質(zhì)量可以用μ=m2/(m1+m2)表示，M1的質(zhì)量為 (1－μ);定義單位時間，使得主天體的運動周期為2π;M1-M2系統(tǒng)的引力常數(shù) (萬有引力常數(shù)G=6.67×10－11m3·kg－1·s－2與系統(tǒng)總質(zhì)量的乘積)定義為1.

用于分析此模型的參照系有很多，本文采用慣性坐標系.在慣性系中，坐標原點為M1和M2的公共質(zhì)心;橫軸x是初始時刻連接M1和M2的連線;縱軸y過原點垂直于x軸并且在兩個主天體的運動平面內(nèi);z軸符合右手定則.在此參考系中，M1和M2的位置分別為

其中，r=(1－e2)/(1+ecosν)，為M1和M2間的距離;ν為M2的真近點角;e為橢圓軌道的偏心率.

M3的動力學方程可表示為

其中，r1，r2分別為M3到M1和M2的距離.

此外，時間和真近點角的關(guān)系為

其中，p=a(1－e2)，為橢圓的半正交弦.

可以求解出此慣性系下，M3的能量EI和角動量CI分別為

其中，RI是M3的矢徑;VI是M3的速度向量.則M3的軌道傾角i=arccos(CIz/|CI|)，其中，CIz為CI的z軸分量.

3 數(shù)值算法

由于橢圓型限制性三體問題目前無法得到軌道變化的解析表達式，所以只能通過數(shù)值積分進行分析.數(shù)值算法流程如下:

1)給參數(shù)Vp，rp，ν，ψ，α，β，e，μ賦初值.

2)計算M3的初始位置和初始速度:

3)在初始條件下，對動力學方程進行向前積分，直到M3與M2的距離大于指定距離d時停止.計算出引力輔助之后飛行器的能量E+和角動量C+.

4)在初始條件下，對動力學方程進行向后積分，直到M3與M2的距離大于指定距離d時停止.計算出引力輔助之前飛行器的能量E－和角動量C－.

本文所有數(shù)值模擬，指定距離d都為0.5，通過變化參數(shù)的初值，可以研究參數(shù)對引力輔助軌道的影響.

4 數(shù)值仿真

根據(jù)軌道動力學的基本定義:軌道能量為“正”表示雙曲線軌道， “負”表示橢圓軌道;角動量z軸分量為“正”表示順行軌道，“負”表示逆行軌道.根據(jù)得到的E－，C－，E+，C+，可以將引力輔助軌道劃分為16種類型，分別用前16個英文字母表示，具體見表1.

表1 英文字母與軌道類型的對應(yīng)關(guān)系

關(guān)于引力輔助參數(shù)Vp，rp，ν，ψ對飛行器軌道的影響，文獻 ［9－13］已經(jīng)進行了討論，在此不再贅述.本文將重點討論參數(shù)α和β對引力輔助前后的飛行器軌道類型、軌道能量和軌道傾角的影響.

4.1 參數(shù)α的分析

4.1.1 α對引力輔助前后軌道類型的影響

本文對大量不同的α取值進行了仿真，特選取了其中較有代表性的情況，見圖2.圖2所示是 α 分別為0°，60°，90°，120°時的 ψ-Vp圖，圖中字母代表對應(yīng)的軌道類型.其他參數(shù)取值分別為:μ=0.012 1，rp=0.004 95，e=0，ν=0，β=0.此組參數(shù)下，0＜ψ＜π圖形與π＜ψ＜2π圖形關(guān)于ψ=π鏡像 (鏡像規(guī)律:B和E，C和I，D和M，G和J，H和N，L和O互為鏡像，A，F(xiàn)，K，P不變)，因此只給出ψ為π～2π時情況;α與360°－α的圖形完全相同，因此只給出α為0°～180°情況.

圖2 α對引力輔助軌道類型的影響

當α=β=0°時，三維橢圓型限制性三體問題就退化為平面橢圓型限制性三體問題，三維引力輔助退化為二維引力輔助，本文結(jié)果與文獻［10］的結(jié)果完全相同，說明本文對三維引力輔助機理的研究可以很好地退化到二維情況.

當α≠0°時，引力輔助為三維情況.從圖2可以看到:隨著α的變化，各種軌道類型區(qū)域均發(fā)生變化.現(xiàn)將其變化規(guī)律總結(jié)如下:

1)橢圓軌道區(qū)域 (引力輔助前后都為橢圓軌道，包含類型A，B，E，F(xiàn))總是出現(xiàn)在左下方.α從0°～90°，整體區(qū)域增大，A，F(xiàn)增大，B減小;α=90°時B消失;α從90°～180°，整體區(qū)域減小，A，F(xiàn)減小，E(B的鏡像類型)增大.

2)雙曲軌道區(qū)域 (引力輔助前后都為雙曲軌道，包含類型K，L，O，P)總是出現(xiàn)在右上方.α從0°～90°，整體區(qū)域增大，K，P增大，L減小;α=90°時L消失;α從90°～180°，整體區(qū)域減小，K、P減小，O(L的鏡像類型)增大.

3)逃逸軌道區(qū)域 (引力輔助前后分別為橢圓軌道和雙曲軌道，包含類型I，J，M，N)總是夾在橢圓軌道區(qū)域和雙曲軌道區(qū)域之間.α從0°～90°，整體區(qū)域縮小，I，J，N均減小;α=90°時均消失;α為90°～180°時，不出現(xiàn).

4)俘獲軌道區(qū)域 (引力輔助前后分別為雙曲軌道和橢圓軌道，包含類型C，D，G，H)總是夾在橢圓軌道區(qū)域和雙曲軌道區(qū)域之間.α為0°～90°時，不出現(xiàn);α從90°～180°，整體區(qū)域增大，C，G，H均增大.

5)逆行軌道區(qū)域 (引力輔助前后都為逆行軌道，包含類型F，H，N，P)總是出現(xiàn)在左上方.α從0°～90°，整體區(qū)域略有增大，F(xiàn)，P增大，N減小;α=90°時N消失;α從90°～180°，整體區(qū)域略有減小，F(xiàn)，P減小，H(N的鏡像類型)增大.

6)順行軌道區(qū)域 (引力輔助前后都為順行軌道，包含類型A，C，I，K)總是出現(xiàn)在右半部和底部 ，一直占據(jù)大部分.α從0°～90°，整體區(qū)域增大，A，K增大，I減小;α=90°時，I消失;α從90°～180°，整體區(qū)域減小，A，K減小，C(I的鏡像類型)增大.

7)逆行變順行軌道區(qū)域 (引力輔助前后分別為逆行軌道和順行軌道，包含類型B，D，J，L)總是夾在逆行軌道區(qū)域和順行軌道區(qū)域之間.α從0°～90°，整體區(qū)域縮小，B，J，L均減小;α=90°時均消失;α為90～180°時，不出現(xiàn).

8)順行變逆行軌道區(qū)域 (引力輔助前后分別為順行軌道和逆行軌道，包含類型:E，G，M，O)總是夾在逆行軌道區(qū)域和順行軌道區(qū)域之間.α為 0°～90°時，不出現(xiàn);α從 90°～180°，整體區(qū)域增大，E，G，O均增大.

4.1.2 α對飛行器能量的影響

圖3所示為α對飛行器能量變化的影響.其中參數(shù)設(shè)置如下:μ=0.012 1，rp=0.004 95，e=0，ν=0，β=0，ψ =20°，Vp=2.6.

圖3 α對軌道能量變化的影響

可以看到，曲線呈鐘形，關(guān)于 α=180°對稱.α=0°或180°時，引力輔助前后飛行器軌道能量改變最大.當α=90°或270°時，Δe=0.本文對其他rp值進行了模擬，所得曲線仍呈鐘形.

圖4 α對軌道傾角變化的影響

4.1.3 α對飛行器軌道傾角的影響

圖4所示為α對飛行器軌道傾角變化的影響.參數(shù)設(shè)置與圖3相同.

可以看到，圖形關(guān)于α=180°對稱.α=0°～180°曲線形狀類似于正弦曲線.α =0°，90°，180°，270°時，Δi=0.本文對其他 rp值進行了模擬，所得曲線形狀類似.

4.2 參數(shù)β的分析

4.2.1 β對引力輔助前后軌道類型的影響

本文對大量的不同β取值進行了仿真，特選取了其中較有代表性的情況，見圖5.圖5所示是β分別為 0°，67.5°，75°，90°，105°，112.5°，180°時的ψ-Vp圖，其他參數(shù)取值分別為:μ=0.0121，rp=0.00495，e=0，ν=0，α=0.此組參數(shù)下，0＜ψ＜π圖形與π＜ψ＜2π圖形關(guān)于ψ=π鏡像 (鏡像規(guī)則與4.1.1節(jié)中相同)，因此只給出ψ為π～2π時情況;β與360°－β的圖形完全相同，因此只給出β為0°～180°情況.

從圖5可以看到:隨著β的變化，各種軌道類型區(qū)域均發(fā)生變化.現(xiàn)將其變化規(guī)律總結(jié)如下:

1)橢圓軌道區(qū)域總是出現(xiàn)在左下方、右下方.β從0°～180°，左下方區(qū)域減小，右下方區(qū)域增大;B，F(xiàn)初始在左下角，不斷減小至消失，后又從右下角出現(xiàn)，并不斷增大.

2)雙曲軌道區(qū)域總是出現(xiàn)在上部.β從0°～90°，整體區(qū)域增大，P，L先增大后減小;β=75°時 P，L消失，只剩 K;β從 90°～180°，整體區(qū)域減小;P，L又從右上角出現(xiàn).

3)逃逸軌道區(qū)域總是出現(xiàn)在橢圓軌道區(qū)域和雙曲軌道區(qū)域之間.β從0°～67.5°，J，N均減小;β=67.5°時J，N消失，只剩下I;β從112.5°～180°，J，N再次出現(xiàn)并不斷增大.

4)俘獲軌道區(qū)域沒有出現(xiàn)在圖中，在ψ為0°～180°中會出現(xiàn).

5)逆行軌道區(qū)域，β為0°～67.5時，出現(xiàn)在左上方;β為0°～67.5°時，出現(xiàn)在右上方.β從0°～67.5°，整體區(qū)域減小;β=67.5°時消失;β從112.5～180°，F(xiàn)，N，P都再次出現(xiàn)并不斷增大.

圖5 β對引力輔助軌道類型的影響

6)順行軌道區(qū)域一直占據(jù)大部分.β從0°～75°，整體區(qū)域增大，從右側(cè)不斷向左側(cè)擴張(主要是K在擴張);β為75°～105°時獨占整個平面;β從105°～180°，整體區(qū)域減小，向左側(cè)縮小 (主要是K在縮小).

7)逆行變順行軌道區(qū)域總是出現(xiàn)在逆行軌道區(qū)域和順行軌道區(qū)域之間.β從0°～75°，向上部移動，β=75°時均消失;β為75°～105°時，不出現(xiàn);β為105°～180°時，B，J，L 再次出現(xiàn).

8)順行變逆行軌道區(qū)域沒有出現(xiàn)在圖中，在ψ為0°～180°中會出現(xiàn).

4.2.2 β對飛行器能量的影響

圖6所示為β對飛行器能量變化的影響.參數(shù)設(shè)置如下:μ=0.0121，rp=0.004 95，e=0，ν=0，α =0，ψ =20°，Vp=2.6.

可以看到，曲線呈鐘形，關(guān)于β=180°對稱;曲線變化幅度很小，說明不同的β引起的軌道能量差別很小.本文對其他rp值進行了模擬，所得曲線仍呈鐘形.

4.2.3 β對飛行器軌道傾角的影響

圖6 β對軌道能量變化的影響

圖7和圖8所示為不同rp下β對飛行器軌道傾角變化的影響.其余參數(shù)設(shè)置與圖6相同.

圖7 rp=0.008時β對軌道傾角變化的影響

圖8 rp=0.00495時β對軌道傾角變化的影響

從圖7可以看出:曲線關(guān)于β=180°對稱.軌道傾角變化的兩個峰值約在β=150°和β=210°處取得，而當β=0°或180°時，軌道傾角變化為0.從圖8可以看出:此參數(shù)下，軌道傾角變化較大，β=180°時，軌道傾角變化為180°，說明軌道發(fā)生順逆轉(zhuǎn)換.本文對其他rp值進行了模擬，發(fā)現(xiàn)β=180°時，使軌道發(fā)生順逆轉(zhuǎn)換的臨界rp值約為0.00536.

綜合圖7和圖8，可以看到:①圖形總是關(guān)于β=180°對稱;②β=180°時，可能的Δi只有兩個值:0°或 180°.

5 以地月系統(tǒng)為例進行參數(shù)優(yōu)化

上述研究對航天器軌道設(shè)計具有重要指導意義.在特定任務(wù)中，往往需要獲得特定軌道類型，如俘獲與逃逸、順行與逆行等，而且往往有一定參數(shù)限制，特別對Vp，rp會有一定限制.如果任務(wù)目標是對M2進行數(shù)據(jù)采集，則往往需要設(shè)計使Vp和rp最小;如果需要躲避M2的大氣層、射線等，則需要使Vp和rp最大.

本文以地月系統(tǒng)為例進行說明.地月系統(tǒng)中，e=0.005 49，μ=0.012 15.限制 rp=0.00504(距月球表200 km)，ν=90°，α =60°，β=30°，限制 ψ 范圍為 180°～360°，分別求逃逸、俘獲的最大Vp.

從圖9可以看出:想要獲得逃逸軌道 (I，J，M，N)，最大 Vp為2.75，軌道類型為 I.想要獲得俘獲軌道 (C，D，G，H)，最大Vp為3.17，軌道類型為H.

圖9 ψ-Vp圖

6 結(jié)論

為了研究三維引力輔助技術(shù)，本文通過引入?yún)?shù)α和β，將平面橢圓型限制性三體模型拓展到三維橢圓型限制性三體模型.深入討論了α和β對飛行器軌道類型、軌道能量和軌道傾角的影響.研究結(jié)果表明:

1)α從0°～90°，橢圓軌道、雙曲軌道、逆行軌道、順行軌道整體區(qū)域都增大，逃逸軌道、俘獲軌道、逆行變順行軌道、順行變逆行軌道整體區(qū)域都減小;α從90～180°，變化趨勢相反.

2)α=90°時，只存在A，F(xiàn)，K，P四種軌道類型，意味著引力輔助前后軌道類型保持不變.α=90°或270°時，引力輔助前后飛行器軌道能量、軌道傾角均保持不變.α=0°或180°時，引力輔助前后飛行器軌道能量改變最大.

3)限定π ＜ψ ＜2π，則只有當0°＜α＜90°時，才存在逃逸軌道和逆行變順行軌道;只有當90°＜α＜180°時，才存在俘獲軌道和順行變逆行軌道.0＜ψ＜π時的規(guī)律，可根據(jù)鏡像規(guī)則得到.

4)β從0°～75°，逆行軌道、逆行變順行軌道、順行變逆行軌道整體區(qū)域都減小，順行軌道整體區(qū)域增大;β從150°～180°，變化趨勢相反.

5)β為75°～105°時，只存在順行軌道，說明引力輔助前后的軌道只可能是順行軌道.

6)不論β取何值，只有當π＜ψ＜2π時，才存在逃逸軌道和逆行變順行軌道;只有當0＜ψ＜π時，才存在俘獲軌道和順行變逆行軌道.

以上的研究和分析揭示了參數(shù)α和β對飛行器軌道特性的影響.通過選取適當?shù)膮?shù)α和β，可以獲得特定任務(wù)所需要的軌道類型、軌道能量和軌道傾角.這些研究可以為引力輔助軌道的設(shè)計提供依據(jù).

References)

［1］Curtis H D.軌道力學 ［M］.周建華，等譯.北京:科學出版社，2009:308－309

Curtis H D.Orbital mechanics for engineering students［M］.Translated by Zhou Jianhua，et al.Beijing:Science Press，2009:308－309(in Chinese)

［2］Flandro G A.From instrumented comets to grand tours:on the history of gravity assist［R］.AIAA 2011-0176，2001

［3］Strange N J，Longuski J M.A graphical method for gravity-assist trajectory design ［J］.Journal of Spacecraft and Rockets，2002，39(1):9－16

［4］Chen K J，Kloster K W，Longuski J M.A graphical method for preliminary design of low-thrust gravity-assist Trajectories［R］.AIAA 2008-6952，2008

［5］Carnelli I，Dachwald B，Vasile M.Evolutionary neurocontrol:a novel method for low-thrust gravity-assist trajectory optimization［J］.Journal of Guidance，Control and Dynamics，2009，32(2):615－624

［6］Pascale P，De Vasile M.Preliminary design of low-thrust multiple gravity-assist trajectories［J］.Journal of Spacecraft and Rockets，2006，43(5):1065－1076

［7］Sims J A，Longuski J M，Patel M R.Aerogravity-assist trajectories to the outer planets and the effect of drag［J］.Journal of Spacecraft and Rockets，2000，37(1):49 －55

［8］Armellin R，le Lavagna M，Ercoli-Finzi A.Aero-gravity assist maneuvers:coupled trajectory and vehicle shape optimization［J］.Celestial Mechanics and Dynamical Astronomy，2006，95:391－405

［9］Broucke R A.The celestial mechanics of gravity assist［C］//AIAA/AAS Astrodynamics Conference.Washington D C:American Institute of Aeronautics and Astronautics，1988:69－78

［10］Prado A F B A.Close-approach trajectories in the elliptic restricted problem ［J］.Journal of Guidance，Control and Dynamics，1997，20(4):797 － 802

［11］喬棟，崔平遠，崔祜濤.基于圓型限制性三體模型的借力飛行機理研究［J］.宇航學報，2009，30(1):82－87

Qiao Dong，Cui Pingyuan，Cui Hutao.Research on gravity-assist mechanism in circular restricted three-body model［J］.Journal of Astronautics，2009，30(1):82－87(in Chinese)

［12］喬棟，崔平遠，尚海濱.基于橢圓型限制性三體模型的借力飛行機理研究［J］.宇航學報，2010，31(1):36－43

Qiao Dong，Cui Pingyuan，Shang Haibin.Research on gravity-assist mechanism in elliptic restricted three-body model［J］.Journal of Astronautics，2010，31(1):36－43(in Chinese)

［13］Felipe G，Prado A F B A.Classification of out-of-plane swing-by trajectories［J］.Journal of Guidance，Control and Dynamics，1999，22(5):643 －649