基于多项式展开的三角平动点垂直周期轨道解析构建方法研究*

翟冠峤 张伟 钱霙婧

(北京工业大学机电学院,北京 100124)

引言

对于三体问题的研究已经进行几百年,最早可以追述到Newton于1687年发表的Principia中.平动点是圆型限制性三体问题中的五个平衡解,包括3个共线平动点L1、L2、L3以及2个三角平动点L4、L5.其中,三角平动点具有“中心×中心”的动力学特性,其附近存在着大量的周期轨道,可以被用于构建空间中转站,编队导航等[1-4].研究这些轨道在深空探测中具有理论价值及工程意义.

为了更加有效地利用三角平动点的动力学特性,在其邻近定点小型探测器,需要在真实力学模型下构造一个稳定的目标轨道,减少探测器在工作过程中的轨道保持需求.显然,以圆型限制性三体问题模型下的线性解提供构造目标轨道的初值选择是不理想的,应用此初值会引起解在积分过程中发散,从而无法保持轨道的周期性,这就需要我们给出更高阶的解.在高阶解的研究中,已有学者对此进行了大量的研究.Richardson[5]应用Lindstedt-Poincaré(L-P)法给出了圆型限制性三体问题共线平动点附近Halo周期轨道的三阶解析解.Erdi[6]和Zagouras[7]基于小参数展开法分别推导了三角平动点附近周期轨道的三阶和四阶解析解.这为动平衡点附近轨道的分析和研究奠定了基础.近期,Lei和Xu[8,9]则通过L-P法构建了三角平动点附近周期轨道的任意阶解析解.但是,这些传统摄动方法主要着重于修正线性条件下的振幅与频率,使其更加接近非线性条件下的真实运动,但是却很少关注运动中各维度之间的联系以及它们分别对系统非线性动力学特性的贡献.

在振动理论最新的进展中,由Shaw等[10-12]基于模态分析的思想提出了一种基于多项式展开理论的求解方法,为周期轨道求解问题提供了新的思想.这种多项式展开方法定义了一种不变的相空间关系[13],从而得到两自由度之间的多项式关系.能为数值求解真实力学模型下的周期轨道提供满足物理规律的约束条件[14].我们之前有过关于三角平动点附近周期轨道的工作,那主要是考虑了x-y平面内长周期或短周期的展开式,并未考虑z轴方向的运动[15].本文则是采用多项式展开的方法求解运动方程,得到圆型限制性三体问题三角平动点附近周期轨道三个维度之间的运动关系.

文中所提出的采用多项式展开方法得到的运动关系可以清晰地反映圆型限制性三体问题模型中三角平动点附近周期轨道三个自由度之间的关系,为分析其轨道动力学特性提供理论依据,并且可以为数值迭代求解周期轨道提供约束关系,为设计真实力学模型下的飞行器轨道提供借鉴.

1 基本动力学模型

在圆型限制性三体问题中,一个质量相对无限小的第三体在两个围绕其公共质心做圆周运动的主天体的引力作用下做运动.

假设质量较大的主天体P1质量为m1,质量较小的主天体P2质量为m2.两个主天体绕其共同的质心C做匀速圆周运动.选取质心会合坐标系进行问题的研究,记为C-XYZ.其原点为C,X-Y平面为两个主天体相对运动平面,X轴由主天体P1指向主天体P2,Z轴垂直于X-Y平面.如图1所示.

图1 坐标系示意图Fig.1 Schematic for coordinate systems

为了计算方便,通常将运动方程无量纲化,取两个主天体质量之和为质量量纲,两个主天体间的距离为长度量纲,即定义μ=m2/(m1+m2)为质量参数.则主天体P1质量表示为1-μ,坐标为(-μ,0,0).主天体P2质量为表示μ,坐标为(1-μ,0,0).小天体在此会合坐标系中的运动方程为:

(1)

式中Ω为系统中的拟势能函数,表示为[16]:

(2)

其中,R1与R2分别代表小天体到主天体P1与P2的距离.

(3)

(4)

2 三角平动点附近的运动方程展开

为了更方便地描述三角平动点附近的运动,将坐标系原点移动到三角平动点.本文以L4点为例,将坐标系原点移动到L4点,新的坐标系的坐标轴与原坐标系的坐标轴平行,如图1所示.在此坐标系下将原运动方程(1)按Legendre展开,可以表示为[17]:

(5)

虽然可以在L4-xyz坐标系中直接求解运动方程的解析解,了解三角平动点附近的运动特征,然而相应的几何特征在该坐标系中并不十分明显,这是由于式(5)的等号左边存在x-y平面的线性耦合项,不利于计算,所以为了更清晰地体现三角平动点附近运动的几何特征,我们选择将原L4-xyz坐标系在x-y平面内绕z轴旋转θ角,得到一个新坐标系L4-ξηζ[16],如图1所示.L4-ξηζ为这一系统的主坐标(principle coordinate system).

引入新变量(ξ,η,ζ)代替(x,y,z),其关系可以表示为:

x=ξcosθ-ηsinθ

y=ξsinθ+ηcosθ

z=ζ

(6)

其中θ满足:

(7)

在新的坐标下,式(5)可以表示为:

(8)

3 周期运动的多项式展开分析

本节使用多项式展开的方法研究三个方向上运动之间的关系,得到它们之间的解析关系,为分析其三个维度分别对系统动力学特性的影响以及三角平动点附近周期运动的运动形式和动力学特性提供参照.这种方法的核心思想是首先选取一个方向为基方向,基上方向的运动状态(位置和速度)为一组空间的二维基状态,将其他方向的运动描述为与基方向状态相关的多项式形式[10-12].通过求解多项式系数的方式寻找多个方向运动之间的关系.

首先,将式化简,可得到运动方程的形式为:

(9)

其中,i,j,k从0开始,n代表截断的阶数,ε代表小参量,

g=2

(10)

令S=sin(θ),C=cos(θ),则:

(11)

其中,Oijk、Pijk、Qijk的具体表达式见附录.

选取ζ方向的位置和速度为周期运动时的基状态,即令:

ζ=u

(12)

(13)

根据文献[10][11]和[12],ξ和η方向的运动可以被描述为与ζ方向相关的如下多项式形式:

ξ=a1u+a2v+a3uv+a4u2+a5v2+

a6u2v+a7uv2+a8u3+a9v3+…

(14)

b6u2v+b7uv2+b8u3+b9v3+…

(15)

η=c1u+c2v+c3uv+c4u2+c5v2+

c6u2v+c7uv2+c8u3+c9v3+…

(16)

d6u2v+d7uv2+d8u3+d9v3+…

(17)

其中,ai、bi、ci、di是待定系数.通过对这些系数的求解,可以得到ξ和η两个方向上位移与速度的与ζ方向上位移与速度的关系.

将式(14)、(15)、(16)与(17)分别带入到式(9)的三个方程中,同时只保留到3次项,可以得到:

(gd3-k1a3-2α200a1a2-2α020c1c2-

α110a1c2-α110a2c1)uv+(gd4-k1a4-α002-

(gd6-k1a6)u2v+(gd7-k1a7)uv2+

(gd8-k1a8)u3+(gd9-k1a9)v3

(18)

(-gb3-k2c3-2β200a1a2-2β020c1c2-

β110a1c2-β110a2c1)uv+(-gb4-k2c4-β002-

(-gb6-k2c6)u2v+(-gb7-k2c7)uv2+

(-gb8-k2c8)u3+(-gb9-k2c9)v3

(19)

(γ101a2+γ011c2)uv-γ003u3

(20)

分别将式(14)、(15)、(16)、(17)对时间求导,可得:

(3a8-2a7)u2v+(2a6-3a9)uv2-a6u3+a7v3

(21)

(-γ101a2b2-γ011b2c2+2b4-2b5)uv+b3v2+

(3b8-2b7)u2v+(2b6-3b9)uv2-b6u3+b7v3

(22)

(3c8-2c7)u2v+(2c6-3c9)uv2-c6u3+c7v3

(23)

(-γ101a2d2-γ011c2d2+2d4-2d5)uv+d3v2+

(3d8-2d7)u2v+(2d6-3d9)uv2-d6u3+d7v3

(24)

分别对比式(15)与式(21)、式(17)与式(23)、式(18)与式(22)、式(19)与式(24)的一次项系数可得到:

u:gd1-k1a1=-b2-a2=b1

-gb1-k2c1=-d2-c2=d1

v:gd2-k1a2=b1a1=b2

-gb2-k2c2=d1c1=d2

(25)

对比二次项系数可得到:

gd3-k1a3-2α200a1a2-2α020c1c2-α110a1c2-

α110a2c1+γ101a2b2+γ011b2c2-2b4+2b5=0

-gb3-k2c3-2β200a1a2-2β020c1c2-β110a1c2-

β110a2c1+γ101a2d2+γ011c2d2-2d4+2d5=0

u2: -γ101a1a2-γ011a2c1-a3-b4=0

-γ101a1c2-γ011c1c2-c3-d4=0

γ101a1b2+γ011b2c1+b3=0

γ101a1d2+γ011c1d2+d3=0

v2:a3-b5=0

c3-d5=0

(26)

对比三次项系数可得到:

u2v:b6=3a8-2a7gd6-k1a6-3b8+2b7=0

d6=3c8-2c7-gb6-k2c6-3d8+2d7=0

uv2:b7=2a6-3a9gd7-k1a7-2b6+3b9=0

d7=2c6-3c9-gb7-k2c7-2d6+3d9=0

u3:b8+a6=0gd8-k1a8+b6=0

d8+c6=0gb8+k2c8-d6=0

v3:b9-a7=0gd9-k1a9-b7=0

d9-c7=0gb9+k2c9+d7=0

(27)

求解式(25)、式(26)、式(27)可以得到线性项系数:

a1=0,a2=0,b1=0,b2=0

c1=0,c2=0,d1=0,d2=0

(28)

二次项系数的表达式为:

(29)

(30)

三次项系数的表达式为:

a7=0,a8=0,a9=0

b7=0,b8=0,b9=0

c7=0,c8=0,c9=0

d7=0,d8=0,d8=0

(31)

其中A=-4k1-4k2+k1k2.

由此得到圆型限制性三体问题中三角平动点附近运动三个自由度之间关系表达式中所有系数,同样,我们也可以将系数求解至更高次项,使结果更加精确.三角点附近周期运动的三个方向ξ、η和ζ之间关系满足如下方程:

(32)

4 周期运动方程解析求解

通过前文对三自由度运动之间关系的求解,可以将运动方程转化为一个自由度的运动方程,将式(32)带入到式(9)的第三个方程中,并忽略3次以上项,可以得到:

(33)

这样,一个三自由度的运动方程就转化为单一自由度的振动方程,我们可以通过传统的摄动方法求解其解析解.本节使用多尺度法进行解析求解.

引入小参数ε,并设方程的解为:

u=u0(T0,T1)+εu1(T0,T1)+…

(34)

其中,T0=t,T1=εt.

微分算子为:

(35)

(36)

将式(34)、式(35)和式(36)带入到式(33)中,令方程两边ε同阶项系数相等,则得微分方程如下:

ε0阶:

D02u0+u0=0

(37)

ε1阶:

(38)

设方程(37)有如下形式的解:

u0=A(T1)eiω0T0+cc

(39)

其中,cc代表共轭项,且ω0=1.

将式(39)带入到式(38)中得:

(40)

消除式(40)中得共振项,可得:

(41)

取:

(42)

则:

(43)

将式(42)和式(43)代入到式(41)中并整理得到:

(44)

分离实虚部得到:

(45)

解得:

(46)

则:

(47)

消除共振项之后的式(40)变为:

(48)

设式(48)的解有如下形式:

(49)

带入到式(48)中可得:

(50)

令式(50)等号左右两边相同项系数相等,可得:

(51)

将式(51)带入到式(49)中即可得到:

(52)

最后可以得到式(33)的解为:

ξ=u

(53)

(54)

(55)

α02a3ω2cos(2ωt+2β0)-

(56)

(57)

α02c3ω2cos(2ωt+2β0)-

(58)

(59)

(60)

5 数值仿真与对比

本节采用地-月-飞行器圆形限制性三体问题中三角平动点附近的周期轨道为例进行仿真计算,其质量参数为μ=0.012150568.同时与传统的摄动方法L-P法得到的解析解[18]进行对比.选定一组参数为幅值α0=0.05,相角β=0.图2(a)～(c)所示为两种方法得到的垂直周期轨道在各个平面上的投影对比,(d)所示为三维坐标下周期轨道,其中实线代表本文方法得到的周期轨道,星线代表L-P法得到的周期轨道.

图2 垂直周期轨道图Fig.2 Vertical Periodic Orbits

通过图2的对比可以发现,在当前无量纲化的方程形式下,两种方法得到的解析解精度接近.且同时得到了垂直周期轨道三个自由度运动的解析关系,通过对比验证了此方法的正确性.

6 结论

本文提出了基于多项式展开法构造圆型限制性三体问题三角平动点附近周期轨道三阶解析解的方法.通过多项式展开方法得到了其运动3个自由度之间的解析关系,为分析其轨道动力学特性提供了理论依据,同时揭示了周期运动的物理规律.

本文所得到的周期运动三阶解析解精度适合,但同时也可继续求解更高阶数、更高精度的解析解,并可以被拓展至椭圆型限制性三体问题模型中,同时也可求解共线平动点的解析解.

1Brouwer D. Clemennce G M. Methods of celestial mechanics(2nd). USA:Academic Press, 1985

2Beuler G. Methods of celestial mechanics, Springer-Verlag Berlin:Heideberg, 2005

3刘林. 人造地球卫星轨道力学. 北京:高等教育出版社, 1992 (Liu L. Orbit Dynamics of the artificial earth satellite. Beijing:The Higher Education Press, 1992 (in Chinese))

4刘林. 航天器轨道理论. 北京:国防工业出版社, 2000 (Liu L. The spacecraft orbit theory. Beijing:National Defence Industry Press, 2000 (in Chinese))

5Richardson D L. Analytic construction of periodic-orbits about the collinear points.CelestialMechanics, 1980,22(3):241～253

6Erdi B. 3-Dimensional motion of trojan asteroids.CelestialMechanics, 1978,18(2):141～161

7Zagouras C G. 3-Dimensional periodic-orbits about the triangular equilibrium points of the restricted problem of 3 bodies.CelestialMechanics, 1985,37(1):27～46

8Lei H L, Xu B. High-order solutions around triangular libration points in the elliptic restricted three-body problem and applications to low energy transfers.CommunicationsinNonlinearScienceandNumericalSimulation, 2014,19(9):3374～3398

9Zhao L. Quasi-periodic solutions of the spatial lunar three-body problem.CelestialMechanics&DynamicalAstronomy, 2014,119(1):91～118

10 Shaw S W, Pierre C. Normal-modes for nonlinear vibratory-systems.JournalofSoundandVibration, 1993,164(1):85～124

11 Shaw S W, Pierre C. Normal-modes of vibration for nonlinear continuous systems.JournalofSoundandVibration, 1994,169(3):319～347

12 Shaw S W. An invariant manifold approach to nonlinear normal-modes of oscillation.JournalofNonlinearScience, 1994,4(5):419～448

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

14 Vakakis A F. Normal Modes and Localization in Nonlinear Systems.SpringerNetherlands, 2001

15 钱霙婧,翟冠峤,张伟.基于多项式约束的三角平动点平面周期轨道设计方法研究. 力学学报, 2017,49(1):154～164 (Qian Y J, Zhai G Q, Zhang W. Planar periodic orbit construction around the triangular libration points based on polynomial constraints.ChineseJournalofTheoreticalandAppliedMechanics, 2017,49(1):154～164(in Chinese))

16 Szebehely V. Theory of orbits. New York and London:Academic Press, 1967

17 Jorba A, Masdemont J. Dynamics in the centre manifold of the collinear points of the restricted three body problem.PhysicaD, 1999,132:189～213

18 刘林,侯锡云. 深空探测器轨道力学. 北京:电子工业出版社, 2012(Liu L, Hou X Y. Orbital mechanics of thedeep space probe. Beijing:Publishing House of Electronics Industry, 2012 (in Chinese))