基于参激共振的受摄小行星悬停轨道设计方法1)

司震 钱霙婧 杨晓东 张伟

(北京工业大学材料与制造学部,机械结构非线性振动与强度北京市重点实验室,北京 100124)

引言

小行星形成于太阳系早期,在几十亿年的演化进程中,完整保留了太阳系早期物质等原始信息,不仅有助于人类探究太阳系起源及生命起源等谜题,也便于人类采取有效的应对措施来防御、减缓或消除碰撞威胁[1].如今,各国正逐步开展小行星的探测采样任务,JASA 的隼鸟号(Hayabusa)探测器于2005 年登陆小行星25143 Itokawa,成为人类第一个小行星采样返回任务[2].自此以后,JASA 的隼鸟II 号(Hayabusa II) 探测器、NASA 的OSIRIS-Rex 探测器相继对小行星162173 Ryugu、101955 Bennu 开展探测采样任务[3-4].

目前,小行星的探测采样任务主要采用悬停探测轨道,但由于探测器所处的动力学环境的差异性,针对小行星探测任务类型或任务的不同阶段,需要对各阶段的动力学进行分析进而设计相应的悬停轨道[5-6].隼鸟号探测器在小行星Itokawa 附近执行任务的大部分时间里都采用惯性系悬停方式,即探测器跟随小行星进行公转运动,在小行星轨道系中进行位置保持.而在小行星表面附近进行采样时,则采取本体系悬停操作,即相对于小行星表面某位置点进行位置保持,以便于在采样返回任务中进行下降和上升机动[7-9].

太阳系中已知小行星自身的引力场普遍呈现十分复杂的不规则特性,其附近探测器的运行轨道受不规则引力场影响[10-11].Scheeres 最早提出了考虑不规则引力场影响下本体系悬停轨道的受控运动[12].Sawai 考虑小行星引力场不规则特性,应用球谐、椭球谐模型、多面体模型等研究了均匀自旋和非均匀自旋小行星附近探测器的本体系悬停轨道[13-14].Broschart 在不考虑第三体引力及外界干扰情况下,对探测器施加控制推力来抵消小行星引力及惯性力,研究了本体系悬停和惯性系悬停的稳定区域[9,15].Zhang 等[16]参照限制性三体问题,考虑了小行星不规则引力对探测器轨道运动的影响,研究了翻滚小行星附近的惯性系悬停控制问题.Zeng 等[17]从非理想太阳帆模型的角度,分析了小行星不规则引力分布对悬停轨道的影响.

由于小行星普遍呈现质量小、引力场弱的特点,研究表明对于在小行星附近悬停运动的航天器而言,太阳辐射光压以及太阳引力摄动不可忽略[18-21].在光压摄动方面,Morrow 等[22]利用太阳光压辐射设计并分析了小行星附近的太阳帆悬停轨道.Williams等[23]在Morrow 等的基础上,根据小行星的球形简化模型,研究了太阳光压辐射影响下惯性系和本体系悬停探测轨道设计问题.Giancotti 和Funase[24]对太阳帆在特洛伊小行星附近的人工平衡点以及附近轨迹进行了可行性研究.Feng 等[25]研究了非球面项C20和太阳辐射压力对小行星附近航天器轨道的影响.Xin 等[26-27]以椭球体小行星模型,考虑太阳光压影响,研究了小行星悬停轨道及其稳定性.在太阳引力摄动方面,刘莹莹等[28]研究了非球形摄动、太阳光压摄动以及太阳引力摄动对绕飞小行星航天器轨道的影响.倪彦硕等[29]研究了考虑太阳引力摄动的小行星附近轨道动力学,结果表明太阳引力对平衡点附近航天器运动的影响较大[29].

由于太阳的周期运动,光压摄动与引力摄动在系统动力学方程中均呈现为周期性激励[30-33].根据文献分析可知,光压摄动呈现外激励的形式,即在系统动力学中表现为非齐次项.然而,引力摄动则呈现参数激励的形式,即在系统动力学中表现为周期性刚度项系数.通常,发生外激励共振时,外激频率需要接近系统的固有频率;而发生参激共振时,参数激励频率需要接近系统的两倍固有频率或与两倍固有频率相关的加减组合关系.由非线性振动理论可知,两种形式的共振并不能同时发生[34-35].因此,本文着重研究小行星不规则引力场与太阳引力摄动下,系统参数激励共振轨道设计问题,为小行星平衡点附近的悬停轨道设计提出新思路.

为了考虑不规则小行星引力场和太阳引力,本文采用受摄粒杆模型,即将小行星视为由无质量杆连接的质量粒子组成[36-39],并加入引力摄动天体部分.本文把太阳引力视为受摄小行星系统的组成部分,而非单纯的摄动源,采用非线性振动理论中的参激共振概念,推导了小行星平衡点附近稳定的参数共振轨道.研究了受摄小行星系统的幅频响应曲线、力频响应曲线以及长短周期运动之间的能量转移现象等非线性特性.此外,由于粒杆模型将小行星视为由无质量杆连接的质量粒子组成,本文的研究成果可以拓展现有小行星附近周期轨道族.

1 受摄不规则小行星系统建模

本文采用粒杆模型模拟小行星的不规则引力场,即根据小行星形状特点及质量分布特征,将小行星系统等效看作由3 个质量粒子(m1,m2和m3) 组成,并假定m3m1=m2.3 个质量粒子由两个恒定距离的无质量杆连接,在空间上呈等腰三角形构型.小行星系统的质心围绕摄动天体M4做圆形轨道运动,且M4到小行星的距离远大于小行星内部每个质量粒子之间的距离.假定所有天体都在同一个平面内运动,且探测器质量忽略不计.

以小行星系统的质心为原点O,建立小行星本体坐标系(O-xyz),受摄粒杆模型示意图如图1 所示.Oz轴沿着小行星的自转角速度Ω 方向,Ox轴平行于从粒子m1指向粒子m2的连线,Oy轴由右手坐标系确定.将系统无量纲化,即设单元长度为粒子m1与m2之间的距离(表示为L),粒子m3与m1的距离在Oy轴上的投影为h,单位质量为M=m1+m2+m3.引入µ=m1/M作为系统的质量参数,则各粒子m1,m2,m3和M4的质量可以表示为

图1 受摄粒杆模型示意图Fig.1 Perturbed particle-linkage model

设定长度比σ=h/L为距离参数.当粒子m3在小行星的质心上方时,σ 数值为正,否则为负.由于小行星的质心以恒定距离Lso绕M4做圆形轨道运动,因此摄动天体M4的位置可以写为(Lsocos θ,Lsosin θ).其中,θ 是Ox轴与摄动天体M4和小行星质心的连线之间的角度(逆时针测量).探测器(x,y)在小行星本体坐标系O-xy中的动力学方程可以表示为

其中,k是无量纲参数,代表小行星表面引力和离心力之比

各质量粒子到探测器的距离ri可以展开为

定义摄动天体的特征参数

其中,以无量纲参数k=1 为边界,k>1 时,引力大于离心力,小行星处于聚集状态;k< 1 时,引力小于离心力,小行星处于离散状态.

观察式(6)可知,摄动天体M4对探测器的引力由两部分组成:由于平均引力而产生的直接引力部分(以线性刚度项表示为kβx/2 和kβy/2),以及因摄动天体位置的周期性变化激励而产生的参数激励部分.相应地,类似著名的Mathieu 方程[30],式(6) 可被分为两部分:第一部分,非扰动系统即小行星系统,表示为

其中,有效势V表示为

以及第二部分,参数激励部分,即式(6) 中包含具有非自治正弦函数和余弦函数的项.值得注意的是,无扰动系统(7)与经典粒杆模型系统并不相同.无扰动系统的有效势函数包含摄动天体的特征参数β,而β最终会影响系统的固有频率.

在本研究中,摄动小行星系统的平衡点坐标(xo,yo)与小行星的质量参数µ,距离参数σ,自旋参数k和摄动参数β 均相关.在给定参数[µ,σ,k,β]的情况下,在xy平面上生成Vx=∂V/∂x和Vy=∂V/∂y的零值曲线,从中可以得出交点的初始值.再通过Newton-Raphson 迭代法得到平衡点的精确位置[40].

当β=0 时,系统退化为经典的粒杆模型.根据Yang 等[36]的研究,可得小行星433 Eros 和243 Ida 的参数,按照上述迭代方法求解平衡点,如图2(a)和图2(b)所示.小行星433 Eros 和243 Ida 都有4 个外部平衡点.当考虑摄动天体的引力即β0 时,小行星433 Eros 和243 Ida 的平衡点的位置将略有改变.

为了分析式(6)中的探测器动力学行为,有必要将坐标系的原点从原有的小行星本体坐标系移动到无扰动系统的平衡点处.相应的,单位长度调整为平衡点与其最接近的质量粒子(m1,m2或m3)之间的距离,表示为γo.这里,选择任意一平衡点Lo作为原点,在新的坐标系中,ξ-和η-的方向与小行星本体坐标系的方向平行,参见图1.需要注意的是,选择作为新坐标系原点的平衡点应该是稳定的.对于以Lo−ξη(o=1,2,3,···)表示的以稳定平衡点Lo为中心的坐标系和小行星本体坐标系O–xy之间的转换关系可以表示为

图2 平衡点随参数变化的分布示意图Fig.2 Distributions of equilibrium points with variations of parameters

根据式(9),将各粒子到探测器的距离ri(i=1,2,3)进行Legendre 展开到四阶,引入非线性项,得到非线性系统.进而可将式(6)表示为

其中,有效势函数V的第二阶导数由下标ξ 和η 表示,上下标中的o 表示求导过程取在平衡点处.Ei,Fi,Di(i=1,2,3)定义为新坐标系下平衡点与各质量粒子的位置关系,其表达式列于附录A 中.由于图1 中系统模型关于y轴对称,对于平衡点L2点和L4点,有D1=D2,E1=−E2,E3=0,F1=F2.M1,M2,···,M7,N1,N2,···,N7是有效势函数V经Legendre 展开获得的多项式系数,其表达式见附录A.ε 为小量符号,用于标记非线性项的量级,在计算中取值为1.

2 参数激励共振条件分析

根据Nayfeh 和Mook[34]的观点,参激共振发生有两个重要的前提条件.首先,参激共振周期轨道设计问题应在稳定平衡点附近进行分析.平衡点的稳定性与其周围轨道的稳定性密切相关,平衡点的稳定性决定了其周围不变流形的稳定性.其次,参数激励频率与系统固有频率应存在相应二倍关系.结合这两个前提条件,可以找到发生参数共振的参数区域.下面依次分析两条件.

条件一:平衡点稳定性条件.对于ξ–η 分量,等式(10)中系统线性部分的特征方程为

其中

如果特征方程(12)具有纯虚部或具有负实部的复数根,则平衡点将是稳定的.因此,需满足以下条件

条件二:参数激励频率与固有频率关系条件.根据式(12),得到平衡点Lo附近的非扰动系统的固有频率为

此外,M4绕小行星系统质心的圆周运动是二体问题,且满足开普勒第三定律.在惯性坐标系中,M4的圆周运动频率Ωs可表示为

由于θ 是图1 中x轴与从质心到M4的线之间的夹角,因此,在平衡点坐标系中,θ 相对于x轴旋转的频率为

这样,可以将式(10)中的θ 改写为θ=ω0t,即

其中,ω=2ω0.

通常,二维系统存在两个固有频率ω1和ω2(假设ω2> ω1),分别对应于长周期和短周期的轨道运动.因此,在理论上可能存在第一类和第二类参激的两种主共振类型,即满足频率关系

以及和型参激共振和差型参激共振的两种组合共振,即满足频率关系

由于平衡点的位置、稳定性以及固有频率均与参数[µ,σ,k,β]有关,而参数激励频率与k,β 参数相关,因此,通过条件一和条件二,可获得满足参激共振条件的参数空间.

为了方便观察摄动参数β 对系统的影响,在仿真中固定参数σ 和k,研究系统参数µ与β 的取值范围.图3 为σ=0.8,k=0.9 时,满足条件一和条件二的µ与β 取值范围.在图3 中,灰色区域表示满足平衡点稳定条件的取值区域,红线表示满足第二类参激主共振条件的取值区域,蓝线则表示满足和型参激共振条件的取值区域.

图3 满足参激共振条件的参数取值区域σ=0.8,k=0.9Fig.3 Parametric resonance region when σ=0.8,k=0.9

图4 不同系统参数下满足参激共振条件的参数取值Fig.4 Parametric resonance regions with different k and σ

图4 为更多组合取值变化条件下的分析结果,图中红线和蓝线分别为满足第二类参激主共振条件以及和型参激共振条件的共振线.为了观察小行星自旋特性对系统共振区域的影响,固定参数σ,改变参数k,得到满足不同参激共振情况的参数取值范围,如图4(a)所示.由图4(a)可知,当k=0.6 时,系统只能发生第二类参激主共振,随着参数k增大,系统同时出现第二类参激主共振以及和型参激共振.

同样,为了观察小行星形状对系统共振区域的影响,固定参数k,改变参数σ,得到满足不同参激共振情况的参数取值范围,如图4(b).在图4(b)中,随着参数σ 减小,满足共振条件所对应的系统µ值增大.因此,对于不规则形状的小行星,形状越接近细长型,且处于聚集状态时,即形状参数σ 越小,自旋参数k越大,系统更容易找到满足参激共振的周期轨道族.需要说明的是,在数值模拟的过程中,没有发现满足第一类参激主共振和差型参激共振关系的参数取值区域.

结合图3 和图4 的分析可知,对于摄动小行星系统,普遍存在第二类参激主共振,在一定条件下,存在和型参激共振,但是普遍不存在第一类参激主共振和差型参激共振.

此外,本文分析了系统存在内共振的可能性,即系统固有频率之间存在关系

同样,以σ=0.8 和k=0.9 时为例,结果如图3 所示.图中,绿线和黑线分别表示满足1:2 和1:3 内共振条件的参数范围.

根据分析可知,系统存在同时发生参激共振以及内共振的可能性.在此区域内,参数激励共振使得系统在微小引力摄动激励条件下产生较大幅值的响应,即参激共振周期轨道,而内共振则可能使得系统整体的能量在长短周期的轨道运动之间转换,由此产生更为复杂的非线性动力学行为.

3 参数激励共振轨道求解

基于上述分析可知,受摄小行星系统普遍存在内共振以及第二类参激主共振,本节以第二类参激主共振和1:3 内共振为例,对参数激励共振轨道进行求解.

为了确定非线性和参数激励对系统振幅和相位的组合影响,采用多尺度方法,假定式(18)具有如下形式的解

其中,上述解均与快时间尺度T0=t以及慢时间尺度T2=ε2t有关.约定算子

把式(22)代入式(18),比较ε 的同次幂,得到

不难发现,式(24)为齐次方程.因此,式(24)的解可以表示为

其中,A1和A2分别为长短周期轨道振幅,cc表示其前面各项的共轭

将式(27)代入式(24),可以获得式(25)中ξ2和η2的解.考虑参激频率和固有频率间的关系以及1:3 的内共振情况,引入调谐参数τ 和κ,将频率之间的关系表示为

由此可得

为了确定式(26)的可解条件,寻求ξ3和η3的特解,以此来消除长期项,可设

然后,将ξ1,η1和ξ2,η2代入式(26)右端,式(31)代入到式(26)左端,并使两侧exp(iω1T0)和exp(iω2T0)的系数相等,得到

且W11,W22,S11,S22的表达式列于附录B 中,上标“—”表示相应参数的共轭.

由此,确定式(27) 的可解性条件的问题被简化为确定线性代数方程(32)的系数行列式为零的条件,考虑式(28)中的关系,可得

把式(33)和式(34)代入式(35),得到

式中,G20与Gr1,Gr3(r=1,2)都是关于系统非线性系数M1,M2,...,M7,N1,N2,...,N7的复数,分离其实虚部为R20与Rr1,Rr2,Rr3,Ir1,Ir2,Ir3,(r=1,2)表达式列于附录B 中,上标“—”表示相应参数的共轭.

将式(36) 中的A1和A2表示为实数ar,βr的形式

结合欧拉公式,分离实部和虚部,并引入新的相位角

把关系式(39)代入式(36),得到一次近似解式(27)中振幅和相位应满足的自治微分方程

其中,方程(40)和(41)明确了系统非线性的影响.非线性项不仅可以直接影响振幅,还可以通过改变相位间接影响振幅.由此,可将系统的一次近似解表示为形式

并根据由式(28)、式(30) 和式(39) 所给出的相位关系

得到参激共振运动的一次近似解的通解形式

其中,ω 是参数激励的频率.需要注意的是,与式(27)不同,式(44)并不是系统的线性解,其包含瞬态和稳态两种运动.振幅ar和相位φr(r=1,2)的取值必须严格遵循方程(40) 和(41),而式(27) 中系统线性解的振幅可以选择为任意值.

此外,由式(44) 可知,长周期轨道运动频率为ω/6,而短周期轨道运动频率为ω/2.由此可见,由于系统非线性的影响,使得短周期轨道的运动频率与参数激励频率之间呈现精确的二倍关系,与此同时,长短周期运动之间呈现精确的1:3 内共振.虽然共振响应轨道同时包含长短周期分量,但是由于频率之间的谐振关系依然呈现周期性.

4 参数激励共振轨道特性分析

4.1 系统稳态解及其稳定性判断

由于通解式(44) 中包含有周期性的稳态运动,显然,考虑到非线性项的影响后,根据式(40) 和式(41) 可知,系统可能存在定常振幅==0 和定常相位==0,即在摄动天体的影响下,参数激励共振响应的振幅不再是无边界,而是存在实现稳态周期轨道运动的可能性.该轨道可为小行星探测任务的轨道设计提供参考.

根据式(40) 和式(41),稳态运动的定常解振幅ar0和相位φr0应满足的代数方程

其中,a10,a20,φ10,φ20表示a1,a2,φ1,φ2对应的稳态解,由于模型的对称性,I20,I21,I22,I23取值为0.

按照a10是否为零,上述方程组的解可分为平凡解和非平凡解两种情况.当a10取平凡解时,系统退化为只与ω2相关的短周期运动.本文考虑一般情况,即a10取非平凡解情况,a100,a200.由式(45)和式(47)可知

即系统存在的4 种相位组合情形.

给定系统参数,采用Newton-Raphson 法针对不同相位组合情况可求得相应振幅稳态解a10,a20,从而得到振幅稳态解与调谐参数τ 之间的关系图,称之为幅频响应曲线.由式(29)可知,调谐参数τ 可以用于表征参数激励频率的大小,由此,通过幅频响应曲线可以分析得知不同参数激励频率所对应的系统振幅稳态解.

此外,为了确保参激共振轨道的稳定性,需要对其长短周期运动进行分析.根据式(45) 至式(48) 可知,长周期运动稳态解的稳定性由式(45)和式(46)的奇点性质判定,短周期运动稳态解的稳定性由式(47)和式(48)的奇点性质判定.考虑到式(49),参激共振轨道的稳定性可以通过式(46)和式(48)的奇点性质来判断.对式(46)和式(48)的奇点性质进行分析,得到用于判断长短周期振幅a10,a20稳定性的Jacobi 矩阵为

矩阵A的特征值可以表示为

其中,p是A的迹,q是A的行列式

当p2−4q> 0 时,λ1,2为互异实特征值.若q< 0,则λ1λ2< 0,奇点为鞍点.若q> 0,p< 0,则λ1λ2>0,奇点为渐近稳定结点.当p2−4q< 0 时,λ1,2为共轭复特征值.若p< 0,奇点为渐近稳定焦点.若p> 0,奇点为不稳定焦点.综上所述,得到奇点渐近稳定的判别条件,即同时满足参激共振轨道长短周期振幅稳定性的条件如下

下面通过数值模拟,对系统稳态解的性质进一步分析,µ=0.001 13,σ=0.8,k=0.9,β=0.014 为例,选取稳定的平衡点L4,对4 种相位组合情况分别讨论,绘制幅频响应曲线,如图5 和图6 所示.

图5 幅频响应曲线(情况1、情况3)Fig.5 Frequency-amplitude response curves(Case 1,Case 3)

图6 幅频响应曲线(情况2、情况4)Fig.6 Frequency-amplitude response curves(Case 2,Case 4)

图6 幅频响应曲线(情况2、情况4)(续)Fig.6 Frequency-amplitude response curves(Case 2,Case 4)(continued)

图5 为情况1(φ10=0,φ20=0)和情况3(φ10=π,φ20=0)的幅频响应曲线.在图5 中,两种情况的曲线重合,黑线为长周期振幅稳态解a10,蓝线为短周期振幅a20.不难发现,稳态解a10,a20均呈现随调谐参数τ 的增大而单调递增的趋势,即一个参数激励频率只对应于一组稳态解.当调谐参数τ=0.000 7 时,系统存在一组稳态解a10=0.014 2,a20=0.000 3,如图5中红点所示.由式(52)和式(54)可得,a10=0.014 2,a20=0.000 3 时,p=−0.010 6 < 0,q=−0.004 8 < 0,λ1=0.004 3,λ2=−0.074 7,该组稳态解是不稳定的.通过式(54) 验证稳态解的稳定性,发现在幅频响应曲线中,任何一组稳态解所对应的q值均小于0,稳态解具有不稳定性.因此,在设计参激共振周期轨道时,情况1 和情况3 将不予以考虑.

图6 为相位组合关系为情况2(φ10=0,φ20=π)和情况4(φ10=π,φ20=π)的幅频响应曲线.当调谐参数τ=0.000 7,系统存在一组稳态解a10=0.013 9,a20=0.000 4,如图6 中红点所示.基于式(52) 和式(54)可得,p=−0.114 1<0,q=0.003 1>0,p2−4q>0,λ1=−0.056 7,λ2=−0.070 0,该组稳态解是稳定的.为了进一步描述调谐参数τ 对稳态解的稳定性的影响,将幅频响应曲线分为3 个区域,标记为I,II,III,对长短周期振幅变化轨迹分别进行标记.

当τ 从0 增大到0.02 时,长周期振幅a10经过A1至A4点,轨迹方向由黑色箭头表示;短周期振幅a20经过B1至B4点,轨迹方向由蓝色箭头表示,如图6(a)所示.在A1、B1点(τ=0)到A2,B2点(τ=0.010 4)的区域I 内,系统只存在一组非平凡解,且通过式(46)和式(48)可分别求出.基于式(54)可得,p<0,q>0,即区域I 中的各组稳态解均是稳定的,如图6 区域I中的实线所示.在A2,B2点到A3,B3点(τ=0.012 3)的区域II 中,长周期振幅a10和短周期振幅a20分别存在两组非平凡解.一组非平凡解所对应的q值小于0,为不稳定稳态解,如图6 区域II 中的虚线所示.一组非平凡解的p< 0,q> 0,则该组稳态解是稳定的,如图6 区域II 中的实线所示.但奇点种类有所不同,当τ ∈(0.010 4,0.011 8)时,p2−4q> 0,λ1,λ2为负实数,呈现“稳定结点”特性,当τ ∈(0.011 8,0.012 3)时,p2−4q< 0,λ1,λ2为实部为负的复数,呈现“稳定焦点”特性.特别注意的是,在A3,B3点,两个非平凡解不再以原轨迹变化,出现“跳跃现象”,a10从A3点跳跃到A4点,a20从B3点跳跃到B4点.在区域III 内,即τ ∈(0.012 3,+∞)时,p< 0,q< 0,呈现“鞍点”特征,该组稳态解是不稳定的,如图6 区域III 中的虚线所示.

当τ 从0.02 减小到0 时,长周期振幅a10经过C1至C4点,轨迹方向由黑色箭头表示; 短周期振幅a20经过D1至D4点,轨迹方向由蓝色箭头表示,如图6(b)所示.在τ 经过C2,D2点到达C3,D3点之前,系统存在两组非平凡解,之后可以观察到跳跃到C4,D4点的现象.与上述τ 增大过程类似,从C4,D4点之后,唯一可实现的稳定的非平凡解由式(46) 和式(48) 给出.除此之外,系统不存在稳定非平凡解.因此,调谐参数的变化方向不同,使得系统稳定非平凡解所对应的调谐参数区域也有所不同.

结合上述长短周期振幅变化的特点,发现了一个有趣的现象.在稳定区域I 和区域II 中,随着调谐参数τ 的增大,长周期振幅a10先增大后减小,而短周期振幅a20不断增大,直至超过长周期振幅a10.根据式(29) 得知,参数激励ω 与系统固有频率ω2之间存在共振关系,通过参数激励,能量不断输入系统,激发系统短周期运动(振幅a20相关运动).在区域I的初始阶段,内共振调谐参数κ < 0.005,使得在区域I 内共振的影响显著,作用于短周期运动的能量通过内共振不断被转移到长周期运动,进而出现长周期振幅a10为主的稳态运动.随着调谐参数τ 不断增大,内共振调谐参数κ 也随之增大,内共振能量转移效应减弱,从而在区域II 内短周期运动的振幅a20不断增大直至超过长周期振幅a10.

为了进一步描述参数激励能量(引力摄动β)对系统的影响程度,以µ=0.001 13,σ=0.8,k=0.9为例,固定调谐参数τ=0.000 7,绘制参数激励振幅1.5kβ 和参激响应振幅ar0的力频关系曲线,如图7所示.

图7 力频响应曲线Fig.7 Amplitude of the excitation–amplitude of response curves

在情况1 和情况3 中,β ∈[0,0.016 0],稳态解随参激振幅的增大而减小.p< 0,q< 0,稳态解呈现“鞍点”特性,该组稳态解是不稳定的,如图7(a)中虚线所示.力频关系分析结果和幅频响应曲线结果一致,因此,在设计参激共振周期轨道时,情况1 和情况3 将不予以考虑.

在情况2 和情况4 中,类似于幅频响应曲线,根据稳态解的个数和稳定性,可以将力频响应曲线分为3 个区域,如图7(b)中虚线所示.其中,实线表示稳定的稳态解,虚线表示不稳定的稳态解.区域边界(红线)为β=0.003 8 (1.5kβ=0.005 1)和β=0.006 9(1.5kβ=0.009 3).当β ∈[0,0.003 8) 时,长周期振幅a10随参激振幅的增大而减小,短周期振幅a20则呈现不断增大,p< 0,q> 0,该稳态解是稳定的.当β ∈[0.003 8,0.006 9] 时,出现两组稳态解,一组稳态解p<0,q<0,是不稳定的;一组稳态解p<0,q>0,是稳定的.当β ∈(0.006 9,0.016 0]时,短周期振幅a20随参激振幅的增大而减小,长周期振幅a10则呈现先增大后减小的特点,p< 0,q> 0,对应的稳态解是不稳定的.力频关系分析结果和幅频响应曲线结论一致.

4.2 参数激励共振响应轨道

通过上述分析判断,根据参激共振响应(44) 以及所求稳定的稳态解,最终得到了满足第二主参数共振和1:3 内共振条件下稳定的参激周期轨道.由于在稳定区域I 和区域II 中,探测器的动力学行为有所差别.因此,将用两组数据分别对参激共振轨道进行模拟求解.

以参数µ=0.001 13,σ=0.8,k=0.9,β=0.014为例,得到调谐参数τ=0.000 7,κ=0.004 1,参激周期轨道长短周期的振幅a10=0.013 9,a20=0.000 4.

图8(a)中红线和蓝线分别表示长短周期运动振幅随时间的变化情况,黑线表示长短周期叠加运动振幅随时间变化曲线.不难发现,由于非线性效应,长短周期轨道频率之间呈现精确的1:3 内共振比例关系.

图8(b)中,蓝线表示摄动天体随时间t在x方向的运动分量,黑线表示参数激励响应即探测器周期轨道随时间t在x方向的运动分量.由于在该组参数条件下,内共振引起的能量转移效应显著,长周期轨道振幅远远大于短周期轨道振幅.虽然共振响应中同时包含有长短周期的分量,但是参激周期轨道几乎完全呈现为长周期运动,如图8(c)所示.此外,由于非线性因素影响,摄动天体频率之间精确保持为1:6的关系.

以参数µ=0.001,σ=0.8,k=0.9,β=0.012 为例,得到调谐参数τ=0.000 6,κ=0.076 7,以及两组稳态解,(1)a10=0.023 5,a20=0.036 9,p< 0,q> 0,p2−4q< 0,稳态解呈稳定特性; (2)a10=0.013 2,a20=0.049 9,p< 0,q< 0,p2−4q> 0,稳态解呈现不稳定特性.取稳定稳态解a10=0.023 5,a20=0.036 9为参激周期轨道长短周期的振幅,绘制共振轨道如图9 所示.图9(a)中红线和蓝线分别表示长短周期运动振幅随时间的变化情况,黑线表示长短周期叠加运动振幅随时间变化曲线.由图可知,该组参数条件下长短周期轨道振幅量级相当,且由于非线性效应,长短周期轨道频率之间呈现精确的1:3 内共振比例关系.图9(b) 中蓝线和黑线分别表示摄动天体和探测器周期轨道随时间t在x方向的运动分量.显然共振响应呈现为长短周期叠加运动,且由于1:3 内共振的影响,使得探测器周期轨道振幅呈现亚谐共振状态.此外由于长短周期轨道频率之间的谐振关系,最终呈现的响应依然呈现周期性,如图9(c)所示.

图8 当σ=0.8,k=0.9,µ=0.001 13,β=0.014 时参激共振轨道示意图Fig.8 Stable parametric resonance orbit when σ=0.8,k=0.9,µ=0.001 13,β=0.014

图9 当σ=0.8,k=0.9,µ=0.001,β=0.012 时参激共振轨道示意图Fig.9 Stable parametric resonance orbit when σ=0.8,k=0.9,µ=0.001,β=0.012

通过上述模拟分析,不难发现,参激周期轨道与系统参数选取密切相关.最终参激周期轨道的形状以及周期关系将严格遵循式(44)的表达形式.

值得一提的是:µ,σ,k是与小行星自身特征相关的参数,β 与摄动天体的引力扰动有关.当参数满足参激共振条件一和条件二时,系统可得到一个第二类参激主共振轨道或一个和型参激共振轨道.此外,参激共振的条件颇为苛刻,并非所有受摄小行星系统都能满足该条件.本文所提出的方法为轨道设计中处理引力摄动提供了新的思路.

5 结论

本文将太阳引力视为受摄小行星系统的组成部分,采用非线性振动理论中的参激共振概念,创新性地设计了小行星平衡点附近稳定的参数共振轨道.

(1)通过引入质量参数µ,形状参数σ,自旋参数k以及摄动参数β,对受摄小行星系统进行建模.研究发现,系统平衡点位置以及稳定性、平衡点附近运动的固有频率均与4 个参数相关.

(2)通过固定参数σ 和k,定性分析了摄动参数β对系统的影响.对于摄动小行星系统而言,普遍存在第二类参激主共振,在一定条件下,存在和型参激共振,但是普遍不存在第一类参激主共振和差型参激共振,且系统普遍存在1:2 和1:3 内共振.

(3)以第二类参激主共振和1:3 内共振为例,采用多尺度方法得到参激共振响应的一次近似解,求解了振幅和相位的稳态解并分析其稳定性,得到由参激共振引起的稳定周期轨道.通过受摄小行星系统的幅频响应曲线和力频响应曲线分析了系统非线性特性.通过分析可知,由于系统非线性的影响,短周期轨道的运动频率与参数激励频率之间呈现精确的二倍关系,与此同时,长短周期运动之间呈现精确的1:3 内共振.虽然共振响应轨道同时包含长短周期分量,但是由于频率之间的谐振关系,使得轨道依然呈现周期性.

此外,由于粒杆模型将小行星视为由无质量杆连接的质量粒子组成,因此本文提出的参激共振轨道设计方法同样适用于双小行星系统或三小行星系统中参数共振轨道动力学问题.本文的研究成果可以拓展现有小行星附近周期轨道族,为轨道设计中处理引力摄动提供了新的思路.

附录A

附录B