鄭丹丹 羅建軍 ,2) 張仁勇 劉 磊

?(西北工業(yè)大學(xué)航天學(xué)院，西安710072)

?(航天飛行動(dòng)力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安710072)

??(中國科學(xué)院空間應(yīng)用工程與技術(shù)中心，北京100094)

??(北京飛行控制中心，北京100094)

基于混合Lie算子辛算法的不變流形計(jì)算1)

鄭丹丹?,?羅建軍?,?,2)張仁勇??劉 磊?,??

?(西北工業(yè)大學(xué)航天學(xué)院，西安710072)

?(航天飛行動(dòng)力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安710072)

??(中國科學(xué)院空間應(yīng)用工程與技術(shù)中心，北京100094)

??(北京飛行控制中心，北京100094)

平動(dòng)點(diǎn)附近周期軌道的不變流形因其在低能軌道轉(zhuǎn)移中起著重要作用而受到廣泛關(guān)注.在設(shè)計(jì)低能軌道過程中不變流形要實(shí)時(shí)進(jìn)行能量匹配，但利用傳統(tǒng)數(shù)值積分方法進(jìn)行積分時(shí)能量會(huì)耗散.顯式辛算法具有比隱式辛算法計(jì)算效率高的優(yōu)勢，但其要求Hamilton系統(tǒng)必須分成兩個(gè)可積的部分，而旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題是不可分的，因而顯式辛算法難以用于求解旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題.本文通過引入混合Lie算子，成功實(shí)現(xiàn)了帶三階導(dǎo)數(shù)項(xiàng)的力梯度辛算法對圓型限制性三體問題的求解，并將基于混合Lie算子的帶三階導(dǎo)數(shù)項(xiàng)的辛算法與Runge-Kutta78算法和Runge-Kutta45算法進(jìn)行仿真對比，仿真結(jié)果表明基于混合Lie算子的含有三階導(dǎo)數(shù)項(xiàng)的辛算法位置精度高、能量誤差小且計(jì)算效率高.利用基于混合Lie算子的帶三階導(dǎo)數(shù)項(xiàng)的辛算法計(jì)算不變流形，可以實(shí)現(xiàn)低能軌道轉(zhuǎn)移過程中軌道拼接點(diǎn)的能量精準(zhǔn)匹配.

圓型限制性三體問題，不變流形，辛算法，混合Lie算子

引言

平動(dòng)點(diǎn)[1]是限制性三體問題中的動(dòng)平衡點(diǎn)，在隨兩個(gè)天體一起旋轉(zhuǎn)的參考系中處于靜止?fàn)顟B(tài).它們是幾何意義上的點(diǎn)，其本身的應(yīng)用價(jià)值十分有限.真正引人關(guān)注的是其周圍大量存在的周期軌道[2]，以及與之相聯(lián)的管狀不變流形[3](包括穩(wěn)定流形和不穩(wěn)定流形).周期軌道是完成某些特殊任務(wù)的理想場所，而不變流形則提供了往返于主天體和周期軌道之間的低能耗轉(zhuǎn)移途徑[49].太陽系中不同三體系統(tǒng)平動(dòng)點(diǎn)周期軌道的不變流形還存在相交的情形，從而拓展了系統(tǒng)間的深空轉(zhuǎn)移軌道設(shè)計(jì)范圍，因此尋找航天器從地球停泊軌道到平動(dòng)點(diǎn)附近周期軌道的轉(zhuǎn)移軌道顯得尤為重要.

20世紀(jì)90年代，G′omez等[10]引入非線性動(dòng)力系統(tǒng)理論，利用穩(wěn)定流形將航天器送到了平動(dòng)點(diǎn)附近，是轉(zhuǎn)移軌道設(shè)計(jì)[11]的重大突破.這一發(fā)現(xiàn)不但節(jié)省了能量也為星際超級公路理論的誕生奠定了基礎(chǔ).但是不論是日--地系統(tǒng)還是地--月系統(tǒng)，不變流形距離地球都比較遠(yuǎn)，需要在推力或引力輔助[12]作用下將航天器從地球停泊軌道轉(zhuǎn)移到穩(wěn)定流形上.由于在優(yōu)化拼接點(diǎn)時(shí)需要多次計(jì)算不變流形，如何選取不變流形的拼接點(diǎn)是進(jìn)行低能軌道轉(zhuǎn)移設(shè)計(jì)的關(guān)鍵.已經(jīng)有學(xué)者對不變流形的計(jì)算進(jìn)行研究，Zhang等[13]利用三次回旋插值法計(jì)算圓型限制性三體問題不穩(wěn)定平動(dòng)點(diǎn)周期軌道的不變流形，Lei等[1417]利用Lindstedt-Poincar′e方法求解了限制性三體和四體平動(dòng)點(diǎn)周期軌道不變流形的高階近似解析解，Dellnitz等[18]提出了方向集數(shù)值計(jì)算方法，Howell等[19]利用胞映射法[20]描述和存儲(chǔ)不變流形的數(shù)據(jù).這些計(jì)算方法提高了不變流形的計(jì)算速度，但是卻沒有關(guān)注不變流形的能量變化.由于圓型限制性三體問題是一個(gè)非線性自治哈密頓系統(tǒng)，沒有嚴(yán)格的解析解，而且由于強(qiáng)非線性其對初值誤差十分敏感，較小的計(jì)算誤差將導(dǎo)致微分方程的解較大地偏離實(shí)際情況.因此需要尋找一種保能量的不變流形計(jì)算方法.

普通數(shù)值解法具有人為耗散性，長時(shí)間計(jì)算會(huì)導(dǎo)致系統(tǒng)總能量的耗散隨時(shí)間的平方增長.在研究三體問題的過程中，對辛結(jié)構(gòu)的破壞無疑會(huì)對動(dòng)力學(xué)特性造成很大影響，使得系統(tǒng)的長期演化不能真實(shí)反映客觀事實(shí).辛積分方法[2122]由于具有保辛結(jié)構(gòu)和能量守恒兩個(gè)重要特性，近年來得到了快速發(fā)展[2326].顯式辛算法具有比隱式辛算法計(jì)算效率高的優(yōu)勢，但其要求Hamilton系統(tǒng)必須可以分成兩個(gè)可積的部分.旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題由于受Coriolis力的影響，不能像慣性系下的系統(tǒng)一樣可分，因此，從形式上看顯式差分格式對其不適用.廖新浩等[27]將圓型限制性三體問題的Hamilton函數(shù)分成動(dòng)能與勢能的和以及由坐標(biāo)旋轉(zhuǎn)產(chǎn)生的非慣性力兩部分，然后利用算法合成構(gòu)造差分結(jié)構(gòu)，計(jì)算復(fù)雜.Sun等[28]分析了帶有三階導(dǎo)數(shù)項(xiàng)的力梯度辛算法在保能量上的優(yōu)勢，但是沒有將該算法應(yīng)用在旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題中，也沒有將辛算法用在不變流形的求解上.

鑒于顯式辛算法不適合求解旋轉(zhuǎn)坐標(biāo)系下圓型限制性三體問題，本文通過重新定義類動(dòng)能的Lie算子，嚴(yán)格證明基于混合Lie算子的含有三階導(dǎo)數(shù)項(xiàng)的顯式辛算法可以求解旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題.首先，給出文中所用的動(dòng)力學(xué)模型和不變流形的計(jì)算方法；然后，建立混合Lie算子并研究顯式辛算法在求解Hamilton系統(tǒng)類動(dòng)能不具有標(biāo)準(zhǔn)二次型的圓型限制性三體問題時(shí)的可行性；最后，將本文改進(jìn)算法與 Runge-Kutta78(RK78)和 Runge-Kutta45(RK45)[29]算法進(jìn)行仿真對比.

1 空間圓型限制性三體問題

1.1 動(dòng)力學(xué)模型

圓型限制性三體問題描述的是質(zhì)量為m3的航天器P3在兩個(gè)繞著共同質(zhì)心做勻速圓周運(yùn)動(dòng)的主天體P1和P2的引力場下運(yùn)動(dòng)，其質(zhì)量分別為m1和m2，且m3?m1?m2，以P1和P2的質(zhì)心為原點(diǎn)，P1指向P2的方向?yàn)閤軸，y軸在兩個(gè)主引力體旋轉(zhuǎn)平面上，z軸垂直于xy平面，滿足右手法則，如圖1所示.

圖1 圓型限制性三體問題Fig.1 Circular restricted three-body problem

假設(shè)P3在oxy平面運(yùn)動(dòng)，為了簡化模型，對單位進(jìn)行無量綱化處理.令引力常量G、主天體P1和P2之間的距離、旋轉(zhuǎn)角速度ω、質(zhì)量和均為1，因此 P1的質(zhì)量為1?μ=m1/(m1+m2)，位置坐標(biāo)為(?μ,0,0)，P2的質(zhì)量為 μ =m2/(m1+m2)，位置坐標(biāo)為(1?μ,0,0).文中所有量如無特殊說明都是無量綱的.在此旋轉(zhuǎn)坐標(biāo)系下航天器的運(yùn)動(dòng)方程和Largrange函數(shù)[30]分別為

和

?(x,y,z)是系統(tǒng)的勢函數(shù)，通常表示為

其中r1,r2分別表示航天器與P1,P2之間的距離

通過Legendre變換

其中qi,pi(i=1,2,3)分別為航天器的坐標(biāo)q=(x,y,z)和動(dòng)量p=(px,py,pz)，可以得到圓型限制性三體問題的Hamilton函數(shù)

該系統(tǒng)具有雅可比能量積分

由Hamilton函數(shù)(2)，可推導(dǎo)出圓型限制性三體問題的正則方程為

1.2 不變流形

不變流形是與平動(dòng)點(diǎn)周期軌道光滑連接的一簇空間軌道，航天器在不變流形上可以無動(dòng)力飛行.因此，不變流形可以作為低能軌道轉(zhuǎn)移的通道，在深空探測軌道轉(zhuǎn)移設(shè)計(jì)中有著重大應(yīng)用價(jià)值.

其中A(t)為系統(tǒng)(1)的雅可比矩陣，?ˉx為相對于不動(dòng)點(diǎn)狀態(tài)的偏移量.由ˉx0點(diǎn)出發(fā)，積分一個(gè)周期后的狀態(tài)所對應(yīng)的狀態(tài)轉(zhuǎn)移矩陣Φ(0,T)稱為單值矩陣[31].周期軌道上不同的點(diǎn)對應(yīng)不同的單值矩陣，每個(gè)單值矩陣有一個(gè)不穩(wěn)定特征根λs(λs＞1)及一個(gè)穩(wěn)定特征根 λu=1/λs(λu＞ 1)，它們分別對應(yīng)著特征向量和，它們包含了穩(wěn)定流形和不穩(wěn)定流形的方向信息.在點(diǎn)處沿特征向量方向增加一個(gè)微小狀態(tài)擾動(dòng)，則可以得到不變流形的積分初值.對于不變流形的計(jì)算有

一般都采用數(shù)值積分來計(jì)算不變流形.不變流形計(jì)算的關(guān)鍵是獲得積分狀態(tài)初值，而此初值和Halo軌道上的點(diǎn)有密切關(guān)系.由于圓型限制性三體問題的強(qiáng)非線性使得微分方程組(1)的解對初始值十分敏感，可能較小的計(jì)算誤差都將導(dǎo)致不變流形的能量誤差巨大.辛算法因具有保能量的特點(diǎn)而受到廣泛關(guān)注，而顯式辛算法具有比隱式辛算法效率高的優(yōu)勢，因此本文考慮利用含有三階導(dǎo)數(shù)項(xiàng)的力梯度辛算法來計(jì)算不變流形.力梯度辛算法是將Hamilton系統(tǒng)分解成動(dòng)能T(p)(僅關(guān)于動(dòng)量p=(px,py,pz)的正定二次函數(shù))和勢能V(q)(關(guān)于坐標(biāo)q=(x,y,z)的函數(shù))兩個(gè)可積部分，然后分別求解.圓型限制性三體問題的Hamilton函數(shù)有多種分解形式，最簡單的是分解成類動(dòng)能T(p)和勢能V(q)兩部分

其中，ypx?xpy是由于Coriolis力的影響而產(chǎn)生的偏移部分.從式(3)可以看出，類動(dòng)能T(p)不具有標(biāo)準(zhǔn)的二次型形式，因此從形式上看力梯度辛算法不適合求解圓型限制性三體問題.后文試圖通過改進(jìn)算法來解決這一難題.

2 基于混合Lie算子的力梯度辛算法

考慮可以分解為動(dòng)能T(p)和勢能V(q)的Hamilton系統(tǒng)

其中，動(dòng)能T(p)僅僅是關(guān)于動(dòng)量p的正定二次函數(shù),即T(p)=p2/2，勢能V(q)是關(guān)于坐標(biāo)q的函數(shù).分別定義T(p)和V(q)的Lie算子[33-34]

由于圓型限制性三體問題Hamilton函數(shù)分解中類動(dòng)能T(p)不是動(dòng)量p的標(biāo)準(zhǔn)二次型，因此具有形式(4)的Lie算子不再適用，考慮將類動(dòng)能T(p)的Lie算子變?yōu)槲恢门c動(dòng)量的混合型算子

它作用于方程組(3)第一個(gè)方程的位置坐標(biāo)和動(dòng)量得到微分方程組

該方程組的分析解為

其中,(x0,y0,z0,px0,py0,pz0)為初始狀態(tài)，t是積分時(shí)間，于是新定義的混合算子ˉA可積，說明重新定義的混合Lie算子是合理的.容易驗(yàn)證算子B也可積.

利用混合Lie算子ˉA和B構(gòu)造復(fù)雜算子

算子D是含有一階、二階和三階導(dǎo)數(shù)項(xiàng)的算子，可以與混合Lie算子ˉA和Lie算子B一起構(gòu)造高階算法

其中,h是積分步長，各積分系數(shù)分別為

其差分格式如下

這里，力 f= ?V/?q.

由

具體證明過程見附錄A,得到

式中，E和F的具體表達(dá)形式見附錄B.這表明算子D是兩主天體的引力梯度，而不是旋轉(zhuǎn)坐標(biāo)系所產(chǎn)生的非慣性力與主天體引力的混合梯度.

因此，通過引入混合Lie算子，分解成形式(3)的圓型限制性三體問題可以用改進(jìn)的顯式辛算法進(jìn)行求解.

3 仿真算例

仿真使用Matlab R2008b軟件，計(jì)算機(jī)為Windows 7系統(tǒng)，配置Intel(R)Core(TM)2 Duo CPU E7500處理器，主頻2.93 GHz，內(nèi)存2.00GB，32位操作系統(tǒng).利用第2節(jié)所得到的改進(jìn)顯式辛算法MTGS研究地--月圓型限制性三體問題并與RK78和RK45法進(jìn)行仿真對比，初始值參數(shù)如表1所示.

表1 初始值參數(shù)Table 1 The initial value of the parameter

圖2 三種算法積分兩個(gè)周期得到的軌道Fig.2 The orbits of three algorithms integral two orbit periods

積分步長取0.001，積分兩個(gè)周期得到的軌道如圖2所示，可以看出本文改進(jìn)的算法和RK78算法得到的軌道仍然能夠閉合，而通過RK45算法得到的軌道開始偏離周期軌道.當(dāng)積分時(shí)長為3個(gè)周期時(shí)，軌道如圖3所示，改進(jìn)的顯式辛算法、RK78算法與RK45算法都不能保證軌道閉合，RK45算法發(fā)散速度最快，改進(jìn)的算法發(fā)散速度最慢.雖然Halo軌道產(chǎn)生偏移的快慢程度與其自身穩(wěn)定性也有一定關(guān)系，但是，在同樣時(shí)間內(nèi)，以同樣步長計(jì)算，不同算法發(fā)生偏移的快慢反映了算法的精確度.

圖3 三種算法積分3個(gè)周期得到的軌道Fig.3 The orbits of three algorithms integral three orbit periods

積分100000步時(shí)絕對能量誤差如圖4所示，可以看出當(dāng)步長取固定值0.001，積分100000步時(shí)改進(jìn)的顯式辛算法的能量誤差最大為10?9量級.雖然RK78算法的初始能量誤差比較小，但隨著時(shí)間推移，能量誤差積累變得越來越大，出現(xiàn)驟增的情況，而RK45算法的能量誤差始終很大.與之相比，顯式辛算法的能量誤差始終保持在一定范圍內(nèi)，能量誤差突然增大是因?yàn)楹教炱髋c主天體P2的距離突然變得非常小.

圖4 三個(gè)算法積分100000步的絕對能量誤差Fig.4 The absolute energy errors of 100000 steps with three algorithms

利用改進(jìn)的顯式辛算法、RK78和RK45算法分別計(jì)算由表1所示初始參數(shù)得到的Halo軌道對應(yīng)的一條不變流形，綜合考慮積分時(shí)間和積分精度，積分步長取0.0001，積分100000步時(shí)三種算法得到的穩(wěn)定流形如圖5所示，不穩(wěn)定流形與其關(guān)于y軸對稱，3種算法得到的不變流形從形狀來看差別無幾，但是不重合.從圖6所示能量相對誤差分析，可知利用改進(jìn)的算法得到的不變流形，可以保證在長時(shí)間積分過程中能量誤差在很小范圍內(nèi)，RK78算法初始能量誤差比較小，但呈增大趨勢，而RK45的能量誤差一直都很大.這3種算法所耗時(shí)長和相對能量誤差的最大值如表2所示，改進(jìn)的顯式辛算法耗時(shí)最短，計(jì)算效率約是RK78的12.3084倍，且相對能量誤差的最大值最小，因此本文改進(jìn)的顯式辛算法在低能軌道轉(zhuǎn)移過程中可以實(shí)現(xiàn)軌道拼接點(diǎn)精準(zhǔn)能量匹配.

圖5 三種算法得到的一條穩(wěn)定流形軌道Fig.5 Stable manifold of three algorithms

圖6 三種算法積分100000步的相對能量誤差Fig.6 The relative energy error of 100000 steps with three algorithms

表2 三種算法的積分效率Table 2 Eefficiency of three algorithms

4 結(jié)論

在基于流形拼接法設(shè)計(jì)航天器低能轉(zhuǎn)移軌道的過程中，需要實(shí)時(shí)進(jìn)行不變流形的能量計(jì)算.本文針對傳統(tǒng)數(shù)值積分能量耗散問題，研究了顯式辛算法在旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題不變流形的保能量問題.顯式辛算法長時(shí)間積分時(shí)沒有長期的變化，能夠保持Hamilton系統(tǒng)的辛結(jié)構(gòu)，本文通過改變圓型限制性三體問題分解形式中類動(dòng)能(動(dòng)能部分不是嚴(yán)格的正定二次型)Lie算子的形式，證明含有一階導(dǎo)數(shù)項(xiàng)、二階導(dǎo)數(shù)項(xiàng)和三階導(dǎo)數(shù)項(xiàng)的改進(jìn)顯式辛算法可以求解旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題，這也說明改進(jìn)的Lie算子形式合理.最后，應(yīng)用改進(jìn)的顯式辛算法、RK78和RK45算法分別計(jì)算了地--月圓型限制性三體問題L1平動(dòng)點(diǎn)附近的周期軌道和不變流形，并對其能量誤差和計(jì)算效率進(jìn)行了分析.仿真結(jié)果顯示，本文改進(jìn)的顯式辛算法在精度和能量誤差兩方面都具有較好的優(yōu)勢.

1 Sezebehely V.Theory of Orbits:The Restricted Problem of Three Bodies.New York:Academic Press，1967

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

3雷漢倫.平動(dòng)點(diǎn)、不變流形及低能軌道.[博士論文].南京：南京大學(xué)，2015(Lei Hanlun.Equilibrium point，invariant manifold and low energy trajectory.[PhD Thesis].Nanjing:Nanjing University，2015(in Chinese))

4 Farquhar RW.The role of Sun-Earth collinear libration points in future space exploration//AAS Annual Meeting,1999

5 Dunham DW，F(xiàn)arquhar RW.Libration point missions，1978-2002//Libration Point Orbits and Applications.Singapore:World Scienti fi c,2003

6 Condon GL,Pearson DP.The role of humans in libration point missionswithspeci fi capplicationtoanEarth-Moonlibrationpointgateway station//AAS 01-307，AAS/AIAA Astrodynamics Specialist Coference,2001

7 Lo MW.The inter-planetary superhighway and the origins program//IEEE Aerospace Conference Proceedings，2002

8 Baoyin HX，McInnes CR.Solar sail Halo orbits at the Sun–Earth arti fi cial L1 point.Celestial Mechanics and Dynamical Astronomy,2006，94:155-17

9 Xu M.Application of Hamiltonian structure-preserving control to formation fl ying on a J2-perturbed mean circular orbit.Celestial Mechanics and Dynamical Astronomy,2012,113(4):403-433

10 G′omez G，Jorba A，Masdenmont J，et al.Study of the transfer from the Earth to a halo orbit around equilibrium point L1.Celestial Mechanics，1993，56:541-562

11袁建平,孫沖,方群.基于虛擬中心引力場方法的航天器轉(zhuǎn)移軌道設(shè)計(jì).力學(xué)學(xué)報(bào),2015,47(1):180-184(Yuan Jianping,Sun Chong,Fang Qun.Spacecraft’s transfer orbit design based on the virtual central gravity?eld method.Chinese Journal of Theoretical and Applied Mechanics,2015,47(1):180-184(in Chinese))

12賈建華,呂敬,王琪.帶脈沖的三維引力輔助變軌研究.力學(xué)學(xué)報(bào),2016,48(2):437-446(Jia Jianhua,L¨u Jing,Wang Qi.Powered gravity assist in three dimensions.Chinese Journal of Theoretical and Applied Mechanics,2016,48(2):437-446(in Chinese))

13 Zhang RY，Luo JJ.Attainable sets approach for low-energy，lowthrustInterplanetarytransfers//InternationalAstronauticalCongress.Beijing，China,2013

14 Lei HL，Xu B，Hou XY，et al.High-order solutions around triangular libration points in the elliptic circular restricted three-body problem and applications to low energy transfers.Communications in nonlinear science and Numerical Simulation，2014，19:3374-3398

15 Lei HL，Xu B.High-order analytical solutions around triangular libration points in the circular restricted three-body problem.Monthly Notices of the Royal Astronomical Society，2013，434(2):1376-1386

16 Lei HL，Xu B，Hou XY，et al.High-order solutions of invariant manifolds associated with libration point in the elliptic circular restricted three-body system.Celestial Mechanics and Dynamical Astronomy，2013，117(4):349-384

17 Lei HL，Xu B.Analytiacl study on the motions around equilibrium points of restricted Four-body problem.Compute Nonlinear Science Number Simulat,2015，29:441-458

18 Dellnitz M，Junge O，Post M.On Target for venus:Set oriented computation of energy efficient low thrust trajectories.Celestial Mechanics and Dynamical Astronomy,2006，95:357-370

19 Howell K，Beckman M，Patterson C，et al.Representations of invariant manifolds for applications in three-body system//14th AAS/AIAA Space Flight Mechanics Conference，Maui，Hawaii，2004

20 Hsu CS.A theory of cell-to-cell mapping dynamical systems.Applied Mechanics，1980,147:931-939

21 Ruth RD．A Canonical integration technique.IEEE Transactions on Nuclear Science，1983，30：2669-2671

22李淵,鄧子辰,葉學(xué)華等.基于辛理論的載流碳納米管能帶分析.力學(xué)學(xué)報(bào),2016,48(1):135-139(Li Yuan,Deng Zichen,Ye Xuehua,et al.Analysing the wave scattering in single-walled carbon nanotube conveying fl uid based on the symplectic theory.Chinese Journal of Theoretical and Applied Mechanics,2016,48(1):135-139(in Chinese))

23 Pauli Pihajoki.Explicit methods in extended phase space for inseparable Hamiltonian problems.Celestial Mechanics and Dynamical Astronomy,2015，121:211-231

24 Liu L，Wu X,Huang GQ,et al.Higher order explicit symmetric integrators for inseparable forms of coordinates and momenta MNRAS459.MonthlyNoticesoftheRoyalAstronomicalSociety,2016 25 Hu WP，Deng ZC，Han SM，et al.Generalized multi-symplectic Integrators for a class of Hamiltonian nonlinear wave PDEs.Journal of Computational Physics，2013，235:394-406

26 Mei LJ，Wu X，Liu FY.On preference of Yoshida construction over Forest–Ruth fourth-order symplectic algorithm.The Europe Physical Journal C,2013,73:2413

27廖新浩，劉林.辛算法在限制性三體問題數(shù)值研究中的應(yīng)用.計(jì)算物理，1995，12(1):102-108(Liao Xinhao，Liu Lin.Applications of symplectic algorithms to the numerical researches of restricted three-body problem.Physic Computation,1995，12(1):102-108(in Chinese))

28 Sun W，Wu X，Hang GQ.Symplectic integrators with potential derivatives to third order. Research in Astronomy and Astrophysics，2011，11:353-368

29 Hairer E，Wanner G.Solving Ordinary Di ff erential Equation.[s.l.]:Springer Verlag，1991

30李言俊，張科，呂梅柏等.利用拉格朗日點(diǎn)的深空探測技術(shù).西安：西北工業(yè)大學(xué)出版社，2014(Li Yanjun，Zhang Ke，L¨u Meibo.et al.Deep Space Exploration Using Lagrange Equilibrium Point.Xi’an:Northwestern Polytechnical Chivesity Press,2014(in Chinese))

31 Bernelli-Zazzera F，Topputo F，Massari M.Assessment of Mission Design Including Utilization of Libration Points and Weak Stability Boundaries.[s.l.]:Politecnico di Milano，2004

32 G′omez G，Jorba A，Masdemont J，et al.Study re fi nement of semianalytical halo orbit theory.Final Report，ESOC Contract No:8625/89/D/MD(SC)，1991

33 Omelyan IP，Mryglod IM，F(xiàn)olk R.Optimized Forest-Ruth-and Suzuki-like algorithms for integration of motion in many-body systems.Computer Physics Communications,2002，146:188-202

34 Omelyan IP，Mryglod IM，F(xiàn)olk R.Symplectic analytically integrable decomposition algorithms:Classi fi cation，derivation,and application to molecular dynamics,quantum and celestial mechanics simulations.Computer Physics Communications，2003，151:272-314

COMPUTATION OF INVARIANT MANIFOLD BASED ON SYMPLECTIC ALGORITHM OF MIXED LIE OPERATOR1)

Zheng Dandan?,?Luo Jianjun?,?,2)Zhang Renyong??Liu Lei?,???(School of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China)
?(National Key Laboratory of Aerospace of Flight Dynamics，Xi’an 710072，China)
??(Technology and Engineering Center for Space Utilization，Chinese Academy of Sciences，Beijing 100094，China)
??(Beijing Aerospace Control Center，Beijing 100094，China)

Invariant manifolds of periodic orbit near the libration points attract a lot of attentions due to their importance in the low-energy orbits transfer problem.In the process of low-energy orbit design,the energy of the invariant manifolds must be matched,but the energy is dissipated when integrating with traditional numerical integration method.The explicit symplectic algorithm with energy conservation is more efficient than the implicit symplectic algorithm,but it requires the Hamiltonian system to be divided into two integral parts,while the circular restricted three-body problem in the rotating coordinate system being inseparable.It is difficult to solve the circular restricted three-body problem in the rotating coordinate system by explicit symplectic algorithm.In this paper,the mixed Lie derivative operator of kinetic energy is used to solve the circular restricted three-body problem in the rotating coordinate system,and the e ff ectiveness of this explicit symplectic algorithm with the third derivation in dealing with this problem has been showed.Compared with the Runge-Kutta45 method and Runge-Kutta78 method,the symplectic algorithm with the third-order derivative term not only has high precision but also the smallest energy error and the highest efficiency.Finally,the invariant manifolds are calculated by the symplectic algorithm with the third derivative term,the patched point can match accurately with this method.

circular restricted three-body problem，invariant manifold，symplectic algorithm，mixed Lie operator

V412.4

A

10.6052/0459-1879-17-079

2017–01–11收稿，2017–07–18 錄用,2017–07–27 網(wǎng)絡(luò)版發(fā)表.

1)國家自然科學(xué)基金重大項(xiàng)目(61690210，61690211)和航天飛行動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目(2015afdl014,2015afd1017)資助.

2)羅建軍，教授，主要研究方向：航天飛行動(dòng)力學(xué)與控制.E-mail:jjluo@nwpu.edu.cn

鄭丹丹,羅建軍,張仁勇,劉磊.基于混合Lie算子辛算法的不變流形計(jì)算.力學(xué)學(xué)報(bào),2017,49(5):1126-1134

Zheng Dandan,Luo Jianjun,Zhang Renyong,Liu Lei.Computation of invariant manifold based on symplectic algorithm of mixed Lie operator.Chinese Journal of Theoretical and Applied Mechanics,2017,49(5):1126-1134

附錄A

附錄B