唐桂琴，高發(fā)寶，王 坤，邱林飛

(揚(yáng)州大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，江蘇 揚(yáng)州 225002)

木星-火星-航天器系統(tǒng)的零速度曲面和轉(zhuǎn)移軌道的數(shù)值研究

唐桂琴，高發(fā)寶，王 坤，邱林飛

(揚(yáng)州大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，江蘇 揚(yáng)州 225002)

針對(duì)木星-火星-航天器系統(tǒng)，根據(jù)歐拉-拉格朗日方程建立了航天器的動(dòng)力學(xué)方程，并通過數(shù)值模擬，研究了系統(tǒng)的零速度曲面和Jacobi常數(shù)的幾何結(jié)構(gòu)的關(guān)系，設(shè)計(jì)了一條可用于多任務(wù)軌道的從火星至木星的轉(zhuǎn)移軌道。

木星；火星；零速度曲面；轉(zhuǎn)移軌道；數(shù)值模擬

0 引言

在深空探測(cè)任務(wù)中，假設(shè)其中一個(gè)質(zhì)量體的質(zhì)量與其他兩個(gè)天體的質(zhì)量相比時(shí)，小到可以忽略，這樣的三體問題就稱為限制性三體問題。這個(gè)小質(zhì)量的天體一般稱為無限小質(zhì)量體，把兩個(gè)大質(zhì)量的天體稱為有限質(zhì)量體。由于無限小質(zhì)量體的質(zhì)量可以看成無限小，所以就可以不考慮它對(duì)另外兩個(gè)有限質(zhì)量體的吸引，也就是說，它不影響這兩個(gè)有限質(zhì)量體的運(yùn)動(dòng)。如果三體中的這兩個(gè)有限質(zhì)量體以一定的角速度繞其公共質(zhì)心作勻速圓周運(yùn)動(dòng)，則此三體問題就稱為圓型限制性三體問題[1-5]。隨著第二輪深空探測(cè)熱潮的到來，火星、木星等天體正逐漸成為各航天大國的探測(cè)重點(diǎn)[6-7]。自1976年火星著陸器第一次成功登陸火星以來，人類已經(jīng)發(fā)射多個(gè)火星探測(cè)器，對(duì)火星的探測(cè)積累了豐富的經(jīng)驗(yàn)。由于地球距離木星要比火星遠(yuǎn)得多，因此人類對(duì)木星的探測(cè)也少得多。但如果借助探索火星的成功經(jīng)驗(yàn)，設(shè)計(jì)一條火星至木星的轉(zhuǎn)移軌道，就可以在探測(cè)火星的同時(shí)探測(cè)木星。當(dāng)然，由于限制性三體問題存在相應(yīng)的零速度曲面，使得無限小質(zhì)量體不能無限制地在宇宙空間中任意飛行。在實(shí)際探測(cè)任務(wù)中，零速度曲面具有極其重要的應(yīng)用價(jià)值，利用零速度曲面不但可以確定火箭飛向月球的最小速度，在討論運(yùn)動(dòng)區(qū)域時(shí)也具有重要的意義，近年來還被用來研究雙星的演化。

本文將深入分析由火星-木星-航天器組成的限制性三體系統(tǒng)的零速度曲面的幾何結(jié)構(gòu)，討論多種情形下Jacobi常數(shù)對(duì)航天器飛行區(qū)域的影響，同時(shí)針對(duì)其中最重要的一種情形，數(shù)值模擬出一條航天器可在火星和木星之間飛行的轉(zhuǎn)移軌道。該軌道可用于多任務(wù)軌道的設(shè)計(jì)，相比獨(dú)立探測(cè)火星和木星來說，將會(huì)減少大量的能源消耗并節(jié)省可觀的費(fèi)用。

1 圓型限制性三體模型

首先介紹兩個(gè)常用坐標(biāo)系：質(zhì)心慣性坐標(biāo)系OXYZ和質(zhì)心轉(zhuǎn)動(dòng)坐標(biāo)系oxyz。質(zhì)心慣性坐標(biāo)系OXYZ中由航天器M，木星M1和火星M2組成的系統(tǒng)，坐標(biāo)系的原點(diǎn)O在M1和M2的質(zhì)心(如圖1所示)，其中X-Y平面是木星-火星系統(tǒng)的軌道平面，M1和M2圍繞它們的公共質(zhì)心作圓周運(yùn)動(dòng)，航天器M只受到M1和M2的引力作用。X軸由M1指向M2，Y軸與之垂直并滿足右手坐標(biāo)系(默認(rèn)Z軸垂直于軌道平面X-Y)，M，M1和M2的質(zhì)量分別為m，m1和m2，且m?m1，m2。

圖1 慣性和轉(zhuǎn)動(dòng)坐標(biāo)系

方便起見，假設(shè)[m]，[L]和[T]分別為質(zhì)量單位、長(zhǎng)度單位和時(shí)間單位，LM1M2為火星和木星的平均距離，ω為火星和木星相對(duì)運(yùn)動(dòng)的角速度，通常按照[m]=m1+m2，[L]=LM1M2，[T]=[LM1M2/(G(m1+m2))]1/2=1/ω對(duì)系統(tǒng)進(jìn)行無量綱化。圖1中θ為M1和M2的連線繞其質(zhì)心轉(zhuǎn)過的角度，系統(tǒng)的質(zhì)量比率μ=m2/(m1+m2)，火星和木星到其公共質(zhì)心的無量綱距離分別為1-μ和μ。質(zhì)心轉(zhuǎn)動(dòng)坐標(biāo)系oxyz的原點(diǎn)o和質(zhì)心慣性坐標(biāo)系的原點(diǎn)O重合，其坐標(biāo)軸x和y以單位角速度圍繞其質(zhì)心作逆時(shí)針轉(zhuǎn)動(dòng)(相對(duì)于慣性坐標(biāo)系中的X軸和Y軸)，并假設(shè)兩個(gè)坐標(biāo)系在時(shí)間t=0時(shí)重合，則有θ=t。如果(x，y，z)表示小天體M在質(zhì)心轉(zhuǎn)動(dòng)坐標(biāo)系中的位置，則有

(1)

將方程(1)兩邊對(duì)時(shí)間t求導(dǎo)，得到

(2)

于是在質(zhì)心轉(zhuǎn)動(dòng)坐標(biāo)系下，系統(tǒng)的拉格朗日函數(shù)為：

(3)

將方程(3)代入歐拉-拉格朗日方程[3]，得到質(zhì)心轉(zhuǎn)動(dòng)坐標(biāo)系oxyz下航天器運(yùn)動(dòng)方程：

(4a)

(4b)

(4c)

定義廣義勢(shì)能Ω(x，y，z)為：

(5)

則質(zhì)心轉(zhuǎn)動(dòng)坐標(biāo)系oxyz中航天器運(yùn)動(dòng)的無量綱方程可改寫為：

(6a)

(6b)

(6c)

2 木星-火星-航天器系統(tǒng)的零速度曲面分析

(7a)

(7b)

(7c)

將方程(7a)-(7c)3式的兩邊分別相加，得

(8)

對(duì)方程(8)兩邊積分得到航天器的運(yùn)動(dòng)狀態(tài)流形：

2Ω(x，y，z)-v2=C，

(9)

由于航天器的飛行區(qū)域由v2≥0決定，因此航天器的運(yùn)動(dòng)區(qū)域?qū)⒂纱肆餍螞Q定。方程(7)在三維空間中所張成的曲面稱為零速度曲面或Hill曲面，在二維平面上則表現(xiàn)為一系列的曲線，稱為零速度曲線。二者的幾何結(jié)構(gòu)將隨著Jacobi常數(shù)C的變化而變化。根據(jù)美國國家航空航天局的官方數(shù)據(jù)，m1=1898.3e+24 kg，m2=0.64174e+24 kg，那么系統(tǒng)的質(zhì)量比率μ=3.379 461 236e-04。

三維空間中，航天器的禁飛區(qū)域或零速度面(2.999 67≤C≤3.023)如圖2所示。當(dāng)Jacobi常數(shù)C的取值逐漸減小時(shí)，表明系統(tǒng)消耗的能量逐漸增大，那么零速度曲面所張成的范圍越來越小，此時(shí)航天器的飛行區(qū)域逐漸開闊；當(dāng)Jacobi常數(shù)C的取值逐漸增大時(shí)，此時(shí)系統(tǒng)消耗的能量逐漸減小，航天器沒有足夠的能量，就像陷入了一個(gè)勢(shì)阱區(qū)域，其飛行的區(qū)域也就越來越小，最終甚至飛不出圖中間木星所在的“大圓”區(qū)域或者火星所在的“小圓”區(qū)域。

圖2 航天器的零速度曲面C∈[2.999 67，3.023]

二維平面中，當(dāng)C=3.023時(shí)，圖2中最上面的區(qū)域由三個(gè)“圓”組成，其中最外圍的“大圓”圍繞著中間的木星所在區(qū)域的“大圓”以及火星所在區(qū)域的“小圓”，且三者呈相離狀態(tài)，說明此時(shí)航天器僅能圍繞木星或者火星飛行，卻不能在兩星之間作穿越飛行，如圖3所示。

圖3 航天器的零速度曲線C=3.023

當(dāng)C由3.023逐漸減小時(shí)，圖3中木星和火星所在的區(qū)域逐漸相互靠近。當(dāng)C減小到3.019 9時(shí)，兩個(gè)區(qū)域連接成了“∞”形。此時(shí)，雖然航天器的飛行區(qū)域有所增加，但“∞”形腰部的那個(gè)交叉點(diǎn)(0.950 359 873 1，0)，卻成了航天器在木星和火星之間作穿越飛行的“要塞”，如圖4所示。圖4局部放大的部分如圖5所示。

圖4 航天器的零速度曲線C=3.019 9

當(dāng)C<3.019 9且逐漸減小時(shí)，連接木星和火星之間飛行區(qū)域的“要塞”被打開，如圖6所示，此時(shí)航天器不僅可以在木星和火星各自的區(qū)域飛行，還可以通過轉(zhuǎn)移軌道從其中一個(gè)星球穿越飛行到另一個(gè)星球。

圖5 航天器的局部零速度曲線C=3.019 9

圖6 航天器的零速度曲線C=3.019 8

當(dāng)C減小到3.019 42時(shí)，木星和火星之間“保齡球”狀飛行區(qū)域與最外圍的“大圓”在臨界點(diǎn)(1.049 041 965，0)處連接在一起，演化出新的“要塞”，如圖7所示，航天器僅能在“保齡球”狀飛行區(qū)域內(nèi)部自由活動(dòng)，但不能穿越飛行到外圍的宇宙空間。

圖7 航天器的零速度曲線C=3.019 42

當(dāng)C繼續(xù)逐漸減小時(shí)，“保齡球”狀飛行區(qū)域與外圍宇宙空間之間新的“要塞”(1.049 041 965，0)被打開。當(dāng)C=3.019時(shí)，探測(cè)器不僅可以在木星和火星之間來回飛行，還可通過打開的新“要塞”去探索外圍的宇宙空間，如圖8所示。

圖8 航天器的零速度曲線C=3.019

當(dāng)C再次逐漸減小至3.000 4時(shí)，航天器的禁飛區(qū)域已退縮至由兩塊彎弧狀在臨界點(diǎn)(-0.995 600 046 3，0)處連接成的區(qū)域，如圖9所示。此時(shí)，航天器雖然可以通過新打開的“要塞”缺口飛向外圍宇宙空間，但卻不能通過兩塊彎弧狀連接點(diǎn)處飛向外圍宇宙空間。隨著C越來越小，航天器的禁飛區(qū)域也逐漸變小，如圖10和11所示。

圖9 航天器的零速度曲線C=3.000 4

圖10 航天器的零速度曲線C=3.000 1

圖11 航天器的零速度曲線C=2.999 7

通過上述數(shù)值模擬結(jié)果，我們發(fā)現(xiàn)當(dāng)Jacobi常數(shù)C的值較大時(shí)，它描繪出一條遠(yuǎn)離原點(diǎn)且形狀近于圓形的閉曲線以及兩條分別圍繞木星和火星的閉曲線；當(dāng)C的值逐漸減小時(shí)，木星和火星附近的兩條閉曲線卻逐漸擴(kuò)大，隨后在交點(diǎn)處相遇，由于在交點(diǎn)處曲線的法線方向不確定，也就是奇點(diǎn)的情況，那么這個(gè)交點(diǎn)實(shí)際上就是平動(dòng)點(diǎn)L1。類似地，當(dāng)C繼續(xù)減小時(shí)，兩條閉曲線在平動(dòng)點(diǎn)L1處融合成一條不斷膨脹的閉曲線并與不斷縮小的外面曲線相遇于平動(dòng)點(diǎn)L2；當(dāng)C再繼續(xù)減小，里外兩曲線變成一條在平動(dòng)點(diǎn)L3處自交的閉曲線；最后，當(dāng)C再度減小時(shí)，曲線分裂成上下兩塊，最終收縮為兩個(gè)平動(dòng)點(diǎn)L4和L5。L1，L2和L3三個(gè)平動(dòng)點(diǎn)，它們?cè)诨鹦呛湍拘堑倪B線上，也稱為共線平動(dòng)點(diǎn)，而L4和L5這兩個(gè)平動(dòng)點(diǎn)關(guān)于坐標(biāo)軸對(duì)稱且分別與火星和木星組成等邊三角形，也稱為三角平動(dòng)點(diǎn)[1，4，8]。

3 木星-火星之間轉(zhuǎn)移軌道的數(shù)值模擬

在實(shí)際的深空探測(cè)任務(wù)中，圖8是我們最為關(guān)注的情形，最具應(yīng)用價(jià)值。因?yàn)槿绻鸍acobi常數(shù)C的取值比圖7中的C還大的話，航天器沒有足夠的能量，脫離不了火星或者木星引力的束縛；如果Jacobi常數(shù)C的取值比圖9中的值還小的話，這不但意味著航天器系統(tǒng)需要消耗大量的能量，而且對(duì)運(yùn)載火箭等也提出了更高的要求，甚至能否成功地將航天器從地球發(fā)射至火星都是個(gè)未知數(shù)。只有當(dāng)Jacobi常數(shù)C的取值介于上述兩種情形之間時(shí)，航天器既能飛離木星-火星系統(tǒng)，又不需攜帶額外的能量。

我們選取適當(dāng)?shù)某跏紬l件[0.999 6，0，0，-39，-0.9，-0.12]，通過數(shù)學(xué)軟件Matlab 2011a對(duì)系統(tǒng)(6)進(jìn)行數(shù)值模擬，得到一條從火星(0.999 6，0，0)出發(fā)，最終到達(dá)木星周圍(-0.000 334 5，0.002 237，-0.002 645)的轉(zhuǎn)移軌道，如圖12所示。

圖12 火星到木星的轉(zhuǎn)移軌道

4 結(jié)論

本文深入討論了火星-木星-航天器系統(tǒng)的零速度曲線、零速度曲面的幾何結(jié)構(gòu)和Jacobi常數(shù)C的關(guān)系。結(jié)果表明：隨著C的取值逐漸增大，雖然此時(shí)系統(tǒng)需要消耗的能量越來越小，但探測(cè)器能夠飛行的區(qū)域也隨之逐漸縮?。环粗?，隨著C的取值逐漸減小，探測(cè)器能夠飛行的區(qū)域也越來越大，連接火星和木星之間的禁飛區(qū)域也將首次在平動(dòng)點(diǎn)L1處打開，此時(shí)環(huán)繞火星(木星)飛行的航天器就可以被木星(火星)捕獲；或者說，火星(木星)航天器將會(huì)逃逸到環(huán)繞木星(火星)飛行的軌道上。由于C的取值越小，系統(tǒng)需要消耗的能量越大，對(duì)運(yùn)載火箭系統(tǒng)的要求也越來越高，使得很多的探測(cè)任務(wù)不能實(shí)施，所以我們最希望關(guān)注的還是二者的最優(yōu)情形，即根據(jù)探測(cè)器實(shí)際探測(cè)任務(wù)的飛行范圍和火箭的實(shí)際運(yùn)載能力，選取適當(dāng)?shù)腏acobi常數(shù)C，使得航天器的飛行范圍盡可能大的同時(shí)系統(tǒng)攜帶的能量盡可能地小。因?yàn)橄牡娜剂显缴伲w行器的總重量也就越輕，飛行器的重量每減少1kg，不但可以使運(yùn)載火箭減輕500kg，還能拓寬發(fā)射窗口，延長(zhǎng)航天器運(yùn)行的壽命。

此外，對(duì)于火星-木星-航天器系統(tǒng)，本文還數(shù)值模擬出一條火星和木星之間的轉(zhuǎn)移軌道，該轉(zhuǎn)移軌道可用于多任務(wù)軌道的設(shè)計(jì)，由于人類探測(cè)火星所積累的經(jīng)驗(yàn)越來越多，因此將大大提高探測(cè)木星的成功率，在探測(cè)火星的同時(shí)探測(cè)木星，相對(duì)于獨(dú)立探測(cè)火星和木星來說，可節(jié)約大量的財(cái)力和物力。

[1] Marchal C. The three-body problem： the restricted problem of three bodies[M]. Amsterdam： Elsevier Publishing Company，1990.

[2] Gao F B， Zhang W. A study on periodic solutions for the circular restricted three-body problem[J]. Astronomical Journal，2014，148(6)：116.

[3] Koon W S， Lo M W， Marsden J E， et al. Dynamical systems， the three-body problem and space mission design[M]. Springer，2011：32-34.

[4] 劉林，侯錫云.深空探測(cè)器軌道力學(xué)[M].北京：電子工業(yè)出版社，2012.

[5] 劉暾，趙鈞.空間飛行器動(dòng)力學(xué)[M].哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2003.

[6] 陳楊，寶音賀西，李俊峰.木星探測(cè)軌道分析與設(shè)計(jì)[J].天文學(xué)報(bào)，2012(2)：106-118.

[7] 李俊峰，寶音賀西.深空探測(cè)中的動(dòng)力學(xué)與控制[J].力學(xué)與實(shí)踐，2007(4)：1-8.

[8] 徐明.平動(dòng)點(diǎn)軌道的動(dòng)力學(xué)與控制研究綜述[J].宇航學(xué)報(bào)，2009(4)：1299-1313.

(責(zé)任編校：夏玉玲)

A Numerical Study of the Zero-Velocity Surface and Transfer Orbit of the Jupiter-Mars-Spacecraft System

TANG Gui-qin， GAO Fa-bao， WANG Kun， QIU Lin-fei

(College of Mathematical Science， Yangzhou University， Yangzhou 225002， China)

The authors of this paper， based on the Jupiter-Mars-Spacecraft system and the Euler Lagrange equation， established a dynamical equation of the spacecraft， carried out numerical simulation， studied the relationship between the zero-velocity surface of the system and the geometric structure of the Jacobi constant and designed a multi-task transfer orbit from Mars to Jupiter.

Jupiter； Mars； zero-velocity surface； transfer orbit； numerical simulation

O242；P135

A

1672-349X(2015)03-0008-04

10.16160/j.cnki.tsxyxb.2015.03.003