錢霙婧 翟冠嶠 張 偉

(北京工業(yè)大學(xué)機電學(xué)院，北京100124)

基于多項式約束的三角平動點平面周期軌道設(shè)計方法1)

錢霙婧 翟冠嶠 張 偉2)

(北京工業(yè)大學(xué)機電學(xué)院，北京100124)

平動點是圓型限制性三體問題中的五個平衡解.其中，三角平動點在平面問題中具有“中心×中心”的動力學(xué)特性，其附近存在著大量的周期軌道，研究這些周期軌道的構(gòu)建方法在深空探測中具有理論及工程意義.本文從振動角度分析周期軌道，通過多項式展開法構(gòu)建出主坐標(biāo)下周期軌道兩個運動方向之間的漸近關(guān)系.從新的角度分析了系統(tǒng)的動力學(xué)特性和平面周期運動兩個方向內(nèi)在關(guān)聯(lián)以及物理規(guī)律.這種多項式形式的關(guān)系式，可以作為約束條件用于數(shù)值微分修正算法中，通過迭代的方式尋找周期軌道.數(shù)值仿真算例驗證了方法的正確性及精確性.文章從振動的角度對周期軌道進行分析，改進了微分修正算法.提出的方法可以被拓展至圓型/橢圓型限制性三體問題的三維周期軌道構(gòu)建中.

平面圓型限制性三體問題，平動點，多項式展開法

引言

平動點是圓型限制性三體問題中的 5個平衡解，包括3個共線平動點L1,L2和L3以及2個三角平動點L4和L5.其中，三角平動點具有“中心×中心”的動力學(xué)特性，其附近存在著大量的周期軌道，可以被用于構(gòu)建空間中轉(zhuǎn)站，編隊導(dǎo)航等[1-4].研究這些軌道對深空探測具有理論價值及工程意義.

目前，已有學(xué)者對于平動點附近的周期軌道構(gòu)建問題進行了解析以及數(shù)值研究.在解析求解方面，Richardson[5]應(yīng)用Lindstedt-Poincar′e(L-P)法給出了圓型限制性三體問題共線平動點附近Halo周期軌道的三階解析解.Erdi[6]和Zagouras[7]基于小參數(shù)展開法分別推導(dǎo)了三角平動點附近周期軌道的三階和四階解析解.這為動平衡點附近軌道的分析和研究奠定了基礎(chǔ).近期，Lei等[8]和Zhao[9]通過L-P法構(gòu)建了三角平動點附近周期軌道的任意階解析解.在發(fā)展解析解的同時，應(yīng)用數(shù)值方法設(shè)計平動點附近周期軌道的思想也在飛速發(fā)展.Goodrich[10]，Bray等[11]，Zagouras等[12]在解析解的基礎(chǔ)上，使用微分修正的方法分別求得了太陽系中不同系統(tǒng)下平動點附近的軌道數(shù)值解.Howell[13]則詳細推導(dǎo)論述了圓型限制性三體問題中擬周期軌道的數(shù)值求解方法，即多步打靶法(multiple shooting).這種方法適用性強，可以根據(jù)不同的終端約束需求做出改動[14-15].隨后，Andreu等[16-17]應(yīng)用多步打靶法，通過初值積分的方式來得到多步打靶的拼接點，最終完成軌道設(shè)計.許多其他學(xué)者也研究了通過利用數(shù)值方法來設(shè)計軌道[18-20].

綜上所述，在解析分析的研究中，傳統(tǒng)攝動方法著重于修正線性條件下的振幅與頻率，使其更加接近非線性條件下的真實運動，但是卻很少關(guān)注運動中各維度之間的聯(lián)系.數(shù)值分析的研究則主要以解析解為軌道的理想猜測初值，通過微分修正以及軌道本身的對稱性來得到周期軌道，或者是以解析解軌道上的點為拼接點(patch points)為約束條件采用多步打靶法來求解連續(xù)軌道.然而，拼接點的求解本身就是迫切需要解決的難點問題[21-22].

在振動理論最新的進展中，由Shaw等[23-25]基于模態(tài)分析的思想提出了的一種基于多項式展開理論的求解方法，為周期軌道求解問題提供了新的思想.這種多項式展開方法定義了一種不變的相空間關(guān)系[26]，從而得到兩自由度之間的多項式關(guān)系.能為數(shù)值求解真實力學(xué)模型下的周期軌道提供滿足物理規(guī)律的約束條件[27].受啟發(fā)于這種思想，本文首先采用多項式展開求解運動方程，得到平面圓型限制性三體問題三角平動點附近周期軌道兩個維度之間的運動關(guān)系.利用微分修正的思想，采用兩個維度之間的運動關(guān)系作為約束條件，進行多次迭代得到滿足誤差精度的周期軌道.

文中所提出的采用多項式展開方法得到的運動關(guān)系可以清晰地反映平面圓型限制性三體問題模型中三角平動點附近周期軌道兩個自由度之間的關(guān)系，為分析其軌道動力學(xué)特性提供理論依據(jù)，并且可以為數(shù)值迭代求解周期軌道提供約束關(guān)系，為設(shè)計真實力學(xué)模型下的飛行器軌道提供借鑒.

1 基本動力學(xué)模型

1.1 平面圓型限制性三體問題

在平面圓型限制性三體問題中，一個質(zhì)量相對無限小的第三體在兩個圍繞其公共質(zhì)心做圓周運動的主天體的引力作用下做運動.

假設(shè)質(zhì)量較大的主天體P1質(zhì)量為m1，質(zhì)量較小的主天體P2質(zhì)量為m2.兩個主天體繞其共同的質(zhì)心C做勻速圓周運動.選取質(zhì)心會合坐標(biāo)系進行問題的研究，記為C-XY.其原點為C，XY平面為兩個主天體相對運動平面，X軸由主天體P1指向主天體P2.如圖1所示.

圖1 坐標(biāo)系示意圖Fig.1 Schematic for coordinate systems

為了計算方便，通常取兩個主天體質(zhì)量之和、兩個主天體間的距離分別為質(zhì)量與長度的單位，即定義μ=m2/(m1+m2)為質(zhì)量參數(shù).則主天體P1質(zhì)量為1-μ，坐標(biāo)為(-μ,0).主天體P2質(zhì)量為μ，坐標(biāo)為(1-μ,0).小天體在此會合坐標(biāo)系中的運動方程為

Ω為系統(tǒng)中的擬勢能函數(shù)，通常表示為[28]

其中，R1與R2分別代表小天體到主天體P1與P2的距離

1.2 三角平動點附近的運動方程展開

為了描述三角平動點附近的運動，將坐標(biāo)系移動到三角平動點.本文以L4點為例，將坐標(biāo)系原點移動到L4，新的坐標(biāo)軸x與y和原坐標(biāo)系X與Y軸平行，如圖1所示.在此坐標(biāo)系下將原運動方程(1)按Legendre展開，可以表示為[29]

式中，Pn為n階的Legendre多項式，

由于式(4)的左邊存在線性耦合項，不利于后續(xù)計算，同時為了更清晰地體現(xiàn)三角平動點附近運動的幾何特征，選擇將原L4-xy坐標(biāo)系旋轉(zhuǎn)θ角，得到一個新坐標(biāo)系L4-ξη[28]，如圖1所示.L4-ξη為這一系統(tǒng)的主坐標(biāo)(principle coordinate system).

引入新變量(ξ,η)代替(x,y)，其關(guān)系可以表示為在新坐標(biāo)下，式(4)可以表示為

由此，消除了方程(4)左端的線性耦合項，線性橢圓運動軌道的長軸和短軸就分別位于這一新坐標(biāo)系的ξ軸和η軸上.

2 周期運動的多項式展開分析

本節(jié)使用多項式展開法研究兩個方向運動之間的關(guān)系，得到它們對系統(tǒng)動力學(xué)特性的影響，為分析三角平動點附近周期運動的運動形式以及動力學(xué)特性提供參照.這種方法的核心思想在于選取一個方向為基方向，基方向上的運動狀態(tài)(位置和速度)為一組空間的二維基狀態(tài)，將另一方向的運動描述為與基方向狀態(tài)相關(guān)的多項式形式[23-25].通過求解多項式系數(shù)的方式尋找兩個方向運動之間的關(guān)系.

將式(7)Legendre展開后表示成如下形式

這里的ξ與η分別代表兩個運動方向的位移.為了得到比較精確的結(jié)果，將式(8)保留前3次項，可將原運動方程改寫為

其中，A2i,B2i,A3i和B3i的具體表達式見附錄.

選取ξ方向的位置和速度為周期運動時的基狀態(tài)，即令

根據(jù)文獻[23-25]，η方向的運動可以被描述為與ξ方向相關(guān)的如下多項式形式

其中,ai與bi是待定系數(shù).通過對這些系數(shù)的求解，可以得到兩個方向上位移與速度的關(guān)系.

將式(15)與式(16)分別代入式(9)的兩個方程中得到

將式(15)對時間求導(dǎo)，可得

將式(16)對時間求導(dǎo)，可得

對比式(16)與式(20)、式(18)與式(22)的一次項系數(shù)可得到

對比二次項系數(shù)可得到

對比三次項系數(shù)可得到

2.1 線性項關(guān)系

求解式(23)，可以得到兩組線性化運動方程的關(guān)系，得到

上述方程表明，對于平面問題，飛行器在三角平動點附近的線性運動是橢圓形的.其一個方向上的位移與另一個方向上的速度成比例，即存在π/2的相位差，當(dāng)其中一個方向上的速度達到最大值時，另一個方向則剛好經(jīng)過其零點.即如果取初始條件滿足式(28)時，會使得原運動方程的其中一個運動分量為零，而只有另外一個運動分量，此時xy平面內(nèi)的運動對應(yīng)為周期運動.而式(26)與式(27)所表示的兩個模態(tài)就對應(yīng)兩類周期軌道：長周期軌道和短周期軌道.

將式(28)代入到式(9)中，可以得到運動方程的線性化頻率

這里的負號對應(yīng)式(26)，正號對應(yīng)式(27).根據(jù)平動點附件周期運動線性穩(wěn)定性分析[30]，在質(zhì)量參數(shù)μ小于μc=0.038520896(Routh極限值)時，上述兩個頻率分別代表了三角平動點附近運動的短周期和長周期運動.

2.2 二次非線性項關(guān)系

繼續(xù)求解式(24)，可以得到運動方程二次非線性系數(shù)的解析表達式.令

2.3 三次非線性項關(guān)系

解得運動方程三次非線性系數(shù)方程(25)的解析形式，可得

其中，S1～S7的表達式見附錄.

由此，用多項式展開的方法得到了平面限制性三體問題中三角平動點附近運動的兩個自由度之間關(guān)系表達式中所有系數(shù).三角點附近周期運動的兩個方向ξ和η滿足如下方程

式(33)和式(34)是描述平面周期運動兩個方向內(nèi)在關(guān)聯(lián)物理規(guī)律的方程，在后文中將用于數(shù)值求解中微分修正的約束條件.

3 微分修正方法

本節(jié)主要介紹軌道微分修正方法，通過前文求得的運動方程兩自由度之間的關(guān)系修正初值，為后文高精度軌道的求解提供理論基礎(chǔ).

微分修正法在本質(zhì)上是一種迭代的打靶法[31]，是對廣義約束和自由變量理論的應(yīng)用.對于這樣的問題，存在著控制變量和約束變量，而這兩種變量則是通過狀態(tài)轉(zhuǎn)移矩陣聯(lián)系起來的.微分修正方法通過狀態(tài)轉(zhuǎn)移矩陣描述約束變量相對控制變量微小改變的敏感性，并通過調(diào)整控制變量使約束變量達到期望值.

求解式(7)的線性部分，容易得到其在L4-ξη坐標(biāo)系下線性解的形式

其中，ω如式(29)所示，α0和β分別代表幅值與相角.ξ和η代表兩個方向上的位移，vξ和vη代表兩個方向上的速度.

確定幅值和相角后，由式(35)得到L4-ξη坐標(biāo)系下微分修正的迭代初值，記為固定時間進行數(shù)值積分得到終端狀態(tài)值運用前文求得的兩個自由度之間的關(guān)系式(33)和式(34)以確定終端目標(biāo)的狀態(tài)值為即

為了保證軌道設(shè)計過程中軌跡滿足設(shè)計目的，Xf應(yīng)該等于的值.即需要滿足如下4個約束條件

對于初始猜測的特解X，可以將F(X)在X0附近進行泰勒展開，并保留其一階項之后表示為

其中，?F(Xj)/?X為狀態(tài)轉(zhuǎn)移矩陣，δX0則為初始猜測X0需要進行修正的量.對于絕大多數(shù)非線性問題而言，需要通過多次迭代的方式得到最終的修正量，因而，聯(lián)立式(37)和式(38)，并用j表示迭代的次數(shù)，可以得出

直到Xj+1=X?或者||F(Xj+1)||≤ε而終止，其中ε為需要的收斂精度.

顯然約束變量等于控制變量，迭代方程有唯一解，可以得到修正量δXj的表達式

這樣，通過多次迭代修正初值量X使其滿足精度要求.

4 數(shù)值仿真與對比

本節(jié)采用地--月--飛行器圓型限制性三體問題中三角平動點附近的周期軌道為例進行仿真計算.其質(zhì)量參數(shù)為μ=0.012150568.取線性解式(35)中，幅值α0=0.05，相角β=π.以此為初值，在平面限制性三體問題全模型，即式(1)條件下，進行微分修正.

4.1 短周期軌道

首先對短周期軌道的情況進行算例仿真，即取ω=0.9545009306377使用線性化模型得到的狀態(tài)初值進行數(shù)值積分，固定積分時間為5個單位化時間，得到積分終值.迭代8次后積分初值滿足精度，迭代完畢.圖2所示為C-XY坐標(biāo)系下，8次修正過程積分初值點的變化情況(僅選取第1次、第3次、第6次、第8次修正過程示意).曲線代表線性解析解軌道，而點劃線表示數(shù)值修正初值所得軌道，星號代表積分初值點.

分別對修正前以及修正后的初值進行50個單位時間積分.圖3(a)所示為線性解積分后得到的軌道圖，圖3(b)所示為微分修正后初值所積分得到的軌道圖.通過對比可以發(fā)現(xiàn)，由線性解積分得到的軌道不能保持周期特性，但由于三角平動點的“中心”特性，其軌跡不會遠離三角平動點，而是呈現(xiàn)擬周期形式的運動.而由修正后的初值進行積分所得到的軌道則可以保持周期運動.數(shù)值仿真驗證了文中方法的有效性.

圖2 短周期軌道微分修正結(jié)果圖Fig.2 Dif f erential correction for the short-period periodic motion

圖3 微分修正前后短周期軌道圖Fig.3 Comparison of the trajectories numerically integrates with the linearized initial condition and corrected initial condition for the short-period motion

4.2 長周期軌道

與短周期軌道的驗證方式相同，通過式(29)得到長周期軌道頻率ω=0.2982079365337，使用線性解(35)得到的狀態(tài)初值.對于長周期軌道，仿真中取固定積分時間為15個單位化時間，得到積分終值.迭代8次后積分初值滿足精度.圖4所示為C-XY坐標(biāo)系下，8次修正過程積分初值點的變化情況(僅選取第1次、第3次、第6次、第8次修正過程示意).曲線代表線性解析解軌道，而點劃線表示數(shù)值修正初值所得軌道，星號代表積分初值點.

圖4 長周期軌道微分修正結(jié)果圖Fig.4 Dif f erential correction for the long-period periodic motion

分別對修正前以及修正后的初值進行50個單位時間積分.圖5(a)所示為線性解積分后得到的軌道圖，圖5(b)所示為微分修正后初值所積分得到的軌道圖.與短周期軌道類似，通過對比可以發(fā)現(xiàn)，由線性解積分得到的軌道不能保持周期特性，但由于三角平動點的“中心”特性，其軌跡不會遠離三角平動點，而是呈現(xiàn)擬周期形式的運動.而由修正后的初值進行積分所得到的軌道則可以保持周期運動.數(shù)值仿真驗證了文中方法的有效性.

圖5 微分修正前后長周期軌道圖Fig.5 Comparison of the trajectories numerically integrates with the linearized initial condition and corrected initial condition for the long-period motion

圖5 微分修正前后長周期軌道圖(續(xù))Fig.5 Comparison of the trajectories numerically integrates with the linearized initial condition and corrected initial condition for the long-period motion(continued)

值得注意的是，文中方法的主要目標(biāo)是尋找平面運動兩個方向上的非線性關(guān)系，是一類幾何約束.這一類約束可以維持周期軌道的幾何特征，但卻無法約束軌道的振幅.長/短周期軌道的仿真中，修正后的軌道振幅都相比較修正前放大.這一特性可以通過添加振幅約束得以改善.

5 結(jié)論

本文提出了基于多項式約束構(gòu)建平面圓形三角平動點附件周期軌道的方法.通過多項式展開方法得到的運動關(guān)系可以清晰地反映平面圓型限制性三體問題模型中三角平動點附近周期軌道兩個自由度之間的運動關(guān)系.為分析其軌道動力學(xué)特性提供理論依據(jù)，并揭示了周期運動的物理規(guī)律.文章創(chuàng)新性的將運動規(guī)律作為微分修正的約束條件，通過多次迭代構(gòu)建平動點附近周期軌道.

本文所使用的研究方法同樣適用于共線平動點附近的軌道構(gòu)建問題，并可以被拓展至圓型/橢圓型限制性三體問題的三維周期軌道構(gòu)建中.

1 Brouwer D,Clemennce GM.Methods of Celestial Mechanics.2nd edn.Academic Press,1985

2 Beuler G.Methods of Celestial Mechanics.Berlin,Heideberg：Springer-Verlag,2005

3 劉林.人造地球衛(wèi)星軌道力學(xué).北京：高等教育出版社,1992(Liu Lin.Orbit Dynamics of the Artificia Earth Satellite.Beijing：The Higher Education Press,1992(in Chinese))

4 劉林.航天器軌道理論.北京：國防工業(yè)出版社,2000(Liu Lin. The spacecraft orbit theory.Beijing：National Defence Industry Press,2000(in Chinese))

5 Richardson DL.Analytic construction of periodic-orbits about the collinear points.Celestial Mechanics,1980,22(3)：241-253

6 Erdi B.3-dimensional motion of trojan asteroids.Celestial Mechanics,1978,18(2)：141-161

7 Zagouras CG.3-dimensional periodic-orbits about the triangular equilibrium points of the restricted problem of 3 bodies.Celestial Mechanics,1985,37(1)：27-46

8 Lei HL,Xu B.High-order solutions around triangular libration points in the elliptic restricted three-body problem and applications to low energy transfers.Communications in Nonlinear Science and Numerical Simulation,2014,19(9)：3374-3398

9 Zhao L.Quasi-periodic solutions of the spatial lunar three-body problem.Celestial Mechanics&Dynamical Astronomy,2014, 119(1)：91-118

10 GoodrichE.Numericaldeterminationofshort-periodtrojanorbitsin the restricted three-body problem.The Astronomical Journal,1966, 71(2)：88-93

11 Bray T,Goudas L.Doubly symmetric orbits about the collinear Lagrangian points.The Astronomical Journal,1967,72(2)：202-213

12 Zagouras C,Kazantzis P.Three-dimensional periodic oscillations generating from plane periodic ones around the collinear Lagrangian points.Astrophysics and Space Science,1979,61(2)：389-409

13 Howell KC.Three-dimensional periodic“halo”orbits.Celestial Mechanics,1984,32(1)：53-71

14 Hughes S,Cooley D,Guzman JA.direct method for fuel optimal maneuversofdistributedspacecraftinmultiplefligh regimes//Space Flight Mechanics Meeting,Copper Mountain,Colorado,2005

15 Marchand B,Howell KC.A spherical formations near the libration points in the sun-earth/moon ephemeris system//14th AAS/AIAA Space Flight Mechanics Conference,Maui,Hawaii,2004

16 Andreu MA,Simo C.Translunar halo orbits in the quasi-bicircular problem//NATO ASI,1997：309-314

17 Andreu MA.Dynamic in the center manifold around L2in quasibicircular problem.Celestial Mechanics and Dynamical Astronomy, 2002,84(2)：105-133

18 Parker J,Born GH.Direct lunar halo orbit transfers//In 17th AAS/AIAA Space Flight Mechanics Meeting,2007,January 28-February 1(AAS)：07-229

19 Perozzi E,Salvo A.Novel spaceways for reaching the moon：an assessment for exploration.Celestial Mechanics&Dynamical Astronomy,2008,102(1-3)：207-218

20 Parker J,Martin L.Shoot the Moon 3D//AAS/AIAA Space Flight Mechanics Meeting,2005,August 7-11(Paper AAS-05-383)

21 Qian YJ,Liu Y,Zhang W,et al.Station keeping strategy for quasiperiodic orbit around Earth-Moon L2 point//Proceedings of the Institution of Mechanical Engineers,Part G：Journal of Aerospace Engineering,2016

22 侯錫云.平動點的動力學(xué)特征及其應(yīng)用.[博士論文].南京：南京大學(xué),2008(Hou Xiyun.Dynamic characteristics and applications of the libration point.[PhD Thesis].Nanjing：Nanjing University, 2008(in Chinese))

23 Shaw SW,Pierre C.Normal-modes for nonlinear vibratory-systems.Journal of Sound and Vibration,1993,164(1)：85-124

24 ShawSW,PierreC.Normal-modesofvibrationfornonlinearcontinuous systems.Journal of Sound and Vibration,1994,169(3)：319-347

25 Shaw SW.An invariant manifold approach to nonlinear normalmodesofoscillation.JournalofNonlinearScience,1994,4(5)：419-448

26 Arquier R.Two methods for the computation of nonlinear modes of vibrating systems at large amplitudes.Computers&Structures, 2006,84(24-25)：1565-1576

27 Vakakis AF,Manevitch L,Mikhlin Y,et al.Normal Modes and Localization in Nonlinear Systems.New York：John Wiley,1996

28 Szebehely V.Theory of Orbits.New York and London：Academic Press,1967

29 Jorba A,Masdemont J.Dynamics in the centre manifold of the collinear points of the restricted three body problem.Physica D, 1999,132：189-213

30 劉林,侯錫云.深空探測器軌道力學(xué).北京：電子工業(yè)出版社,2012 (Liu Lin,Hou Xiyun.Orbital Mechanics of the Deep Space Probe. Beijing：Publishing House of Electronics Industry,2012(in Chinese))

31 Keller HB.Numerical solution of two point boundary value problems.Society for Industrial and Applied Mathematics,1976

附錄

PLANAR PERIODIC ORBIT CONSTRUCTION AROUND THE TRIANGULAR LIBRATION POINTS BASED ON POLYNOMIAL CONSTRAINTS1)

Qian Yingjing Zhai Guanqiao Zhang Wei2)
(College of Mechanical Engineering，Beijing University of Technology，Beijing100124，China)

Libration points are the fve equilibrium solutions in the circular restricted three-body problem(CRTBP). The linearized motions around triangular libration points are typical center×center type.Studies about probes moving around orbits in the vicinity of the libration points have theoretical significance From the vibrational point of view,the polynomial series are used to derive approximately the relations in dif f erent directions during periodic motions,which provides a new point of view to exploring the dynamics and analyzing the overall characteristics of the whole system with general rules.The nonlinear relations in polynomial form between the directions of the planar motions can be treated as constraints to obtain the solutions by numerical integration.Numerical simulations verify the efficiency of the proposed method.The methodology of deriving topological relations has the potential to be extended to circular/elliptical R3BP in three dimensional cases.

planar circular restricted three-body problem，libration point，polynomial expansion method

V412.4

A doi：10.6052/0459-1879-16-215

2016-08-01收稿，2016-10-11錄用,2016-10-18網(wǎng)絡(luò)版發(fā)表.

1)國家自然科學(xué)基金資助項目(11402007).

2)張偉，教授，主要研究方向：非線性動力學(xué).E-mail：sandyzhang0@163.com

錢霙婧,翟冠嶠,張偉.基于多項式約束的三角平動點平面周期軌道設(shè)計方法.力學(xué)學(xué)報,2017,49(1)：154-164

Qian Yingjing,Zhai Guanqiao,Zhang Wei.Planar periodic orbit construction around the triangular libration points based on polynomial constraints.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1)：154-164