406 天文学进展 36 卷摄动运动方程描述, 还是采用运动方程描述, 两者对坐标系选取与对方程的展开和解算方法等都存在差异这些模型在一般的教科书中都有介绍, 如由 Brouwer 和 Clemence 所著的经典专著 Methods of Celestial Mechanics [1],

Size: px

Start display at page:

Download "406 天文学进展 36 卷摄动运动方程描述, 还是采用运动方程描述, 两者对坐标系选取与对方程的展开和解算方法等都存在差异这些模型在一般的教科书中都有介绍, 如由 Brouwer 和 Clemence 所著的经典专著 Methods of Celestial Mechanics [1],"

治屠
5 years ago
Views:

1 第 36 卷第 4 期天文学进展 Vol. 36, No 年 11 月 PROGRESS IN ASTRONOMY Nov., 2018 doi: /j.issn 限制性三体问题中摄动运动方程的坐标系选择沈欣和 1,2,3, 王文磊 1,2,3, 许雪晴 2,3, 周永宏 1,2,3, 廖新浩 2,3 (1. 中国科学院大学, 北京 ; 2. 中国科学院行星科学重点实验室, 南京 ; 3. 中国科学院上海天文台, 上海 ) 摘要 : 针对限制性三体问题, 分别选取以中心天体和摄动体质心为坐标原点的惯性系, 及以中心天体为坐标原点的非惯性系, 讨论了不同坐标系下天体运动轨道描述的异同利用运动天体轨道能量 E 的大小, 可以确定受摄运动方程采用椭圆轨道根数还是采用双曲线轨道根数进行描述为此, 推导出一个关于轨道半长径和偏心率满足的临界关系判别式结果表明, 在摄动天体质量较大的情况下, 非惯性系中存在大量轨道, 这些轨道在原惯性坐标系中是稳定的椭圆轨道, 转换到非惯性系中后却无法用椭圆轨道根数进行描述只能引入双曲线轨道根数来描述轨道, 由此将产生非惯性系下摄动运动方程轨道根数类型选择问题最后, 指出选择雅可比坐标系可以避免上述问题, 并推导出适用于任意运动区域的具有统一形式的摄动函数展开式关键词 : 限制性三体问题 ; 椭圆运动 ; 雅可比坐标系 ; 摄动运动方程中图分类号 : P133 文献标识码 : A 1 引言对自然界中一切运动的描述, 都必须在某一特定的参考空间中进行, 如果选用的参考系不同, 对同一运动的描述也会不同我们知道, 牛顿力学的动力学方程仅适用于惯性参考系对非惯性参考系, 需考虑非惯性力效应, 才能给出正确的动力学描述换言之, 无论在惯性或非惯性系, 只要我们全面考虑惯性力和非惯性力, 无论采用何种坐标系, 我们都可以得到正确的动力学方程基于太阳系或其他天体系统的分布特性, 在限制性三体问题的理论分析研究中, 涉及最多的是两类动力学模型 : 摄动运动模型和运动模型这两类模型的主要差别体现在采用收稿日期 : ; 修回日期 : 资助项目 : 国家自然科学基金 ( , , ) 通讯作者 : 沈欣和,shenxinhe@shao.ac.cn

2 406 天文学进展 36 卷摄动运动方程描述, 还是采用运动方程描述, 两者对坐标系选取与对方程的展开和解算方法等都存在差异这些模型在一般的教科书中都有介绍, 如由 Brouwer 和 Clemence 所著的经典专著 Methods of Celestial Mechanics [1], 由 Murray 和 Dermott 所著的新出版的教材 Solar System Dynamics [2] ; 而专门介绍运动模型的书籍是 Szebehely 所著的 Theory of Orbits [3], 以及由 Marchal 等人所著的 The Three-Body Problem [4] 为了便于理论分析, 更加直观地研究轨道的动力学特性, 人们惯于采用摄动运动方程, [5] 例如经典的拉格朗日行星摄动运动方程和高斯型受摄运动方程等人们通常选取椭圆的轨道根数 a, e, i, Ω, ω 和 M 作为方程的基本变量, 它们分别表示椭圆轨道半长径偏心率倾角升交点角近点角距和平近点角因为抛物线和双曲线轨道根数对应不稳定的运动, 所 [6] 以在对各类行星轨道的实际研究中很少被采用, 这类轨道根数可用于彗星轨道的研究显然, 对于行星动力学的研究, 这里采用的椭圆轨道根数, 是以大质量中心天体作为椭圆的焦点, 同时该中心天体也是所选取坐标系的原点由此可见, 摄动运动方程对应的坐标系一般都是非惯性的, 其采用的椭圆轨道根数也是非惯性系中几何变量当受摄变量的运动方程能够恰当地描述天体的运动时, 惯性系和非惯性系中建立的两组模型给出的结论是一致的, 那么无论采用哪种模型, 对天体进行运动分析都是合理的问题是惯性系中椭圆轨道在非惯性系中是否仍是椭圆的? 如果答案是否定的, 那么就会出现这样一种情况 : 惯性空间的初始瞬时椭圆运动, 因为采用了非惯性系, 在新的坐标系中变为瞬时双曲线或抛物线运动, 所以椭圆受摄运动方程对其不再适用如果发生这样的情况, 那么在使用非惯性系时, 可能会在轨道的定性分析中得出错误的结论因此, 我们将分析不能采用非惯性系的摄动运动模型本文在第 2 章中基于限制性三体问题模型理论分析了在惯性和非惯性坐标系下描述椭圆运动的异同, 推导出一个临界关系判别式为了解决限制性三体问题的坐标系选择问题, 第 3 章在雅可比坐标系下推导出了适用于任意运动区域的受摄二体分析的新模型, 并对相应摄动函数进行了分析展开最后本文对所研究的内容进行了总结 2 惯性与非惯性坐标系下椭圆轨道描述的异同考虑一个三体系统, 由主天体 S, 摄动天体 P 和小天体 A 组成,3 个天体的质量分别记为 m 1, m 2 和 m, 且满足 m 1 > m 2 > 0 和 m = 0 其中主天体 S 和摄动天体 P 在惯性坐标系中始终保持圆周运动取两者的质量中心 O 为坐标原点, 建立一个惯性坐标系, 记为 O-XY Z 同时, 也可以建立一个非惯性系, 取主天体 S 为坐标原点 o, 将其记为 o-xyz, 如图 1 所示令 m 1 + m 2 为质量单位, 摄动天体 P 在坐标系 o-xyz 中的轨道半长径为长度单位设 m 2 = µ, 则 m 1 = 1 µ, 万有引力常数 G = 1, 则时间单位可以直接导出为了定性研究坐标系转换造成轨道描述的问题, 我们首先考虑最简单的小天体初始轨道为平面圆轨道情况, 然后得到主天体和摄动天体在惯性系中的直角坐标的表示形式 ( 假设时间 t = 0 时,S 和 P 位于 x 轴上 ), 对应的位置矢量分别为 R 1 和 R 2, 分别表示为 :

3 4 期沈欣和, 等 : 限制性三体问题中摄动运动方程的坐标系选择 407 注 : 位置矢量 r 表示小天体 A 与非惯性系坐标原点 o ( 主天体 S) 的相对位置, 位置矢量 R 1, R 2 和 R 分别表示天体 S, P 和 A 与惯性坐标系原点 O ( 主天体 S 与摄动天体 P 的质量中心 ) 的相对位置,ψ 为 OP 与 OA 之间的夹角图 1 三体系统坐标示意图 R 1 = R 2 = µ cos t µ sin t 0 (1 µ) cos t (1 µ) sin t 0 = X 1 Y 1 Z 1 = X 2 Y 2 Z 2, (1). (2) 为了研究惯性系与非惯性系下天体轨道运动的异同, 在惯性系中假设小天体 m 的初始状态是在惯性坐标系中作半径为 R 的平面圆周运动, 对应的能量为 1 2 R 令其平近点角 M = nt + M 0,n 为平均运动角速度, 则相应的直角坐标可以表示为 : { X = R cos M Y = R sin M. (3) 将上述惯性系 O-XY Z 中的小天体的运动转换到非惯性系 o-xyz, 相应的坐标可以通过如下坐标变换获得, 即 { x = X X1 y = Y Y 1, (4) 其中,(X 1, Y 1 ) 为中心天体的坐标, 速度 v 只要对时间求导数即可获得惯性坐标系中的椭圆运动被转换到非惯性系中后, 三体系统的能量为 : E = 1 2 υ2 1 µ r. (5) 图 2 图 3 的小天体轨道均对应惯性坐标系中的椭圆轨道, 从图中可以看出, 当 µ 取值不同时, 惯性系中的椭圆轨道转换到非惯性系后, 其轨线形状和能量也不同当 µ = 10 3

4 408 天文学进展 36 卷时, 能量 E 随时间变化均小于 0, 因此可以用变化的椭圆描述其运动 ; 但 µ = 0.4 时, 出现了能量大于 0 的情况, 因此不能用变化的椭圆去描述其运动了出现此差别的原因可以用图 2 和 3 的轨线性质给出清楚的解释图 2 的轨线并没有像图 3 中轨线出现的回转, 这里的回转代表轨道的能量由正变到负的情况, 详见图 4 从图中可以看出, 惯性系中稳定的椭圆轨道转换到非惯性系中后, 如果要使用轨道根数对轨道进行描述, 轨道将在椭圆运动和双曲线运动中不停地转换, 因此即便在某一时刻轨道变为双曲线运动, 也无法判定轨道是否已经失去稳定性上述惯性和非惯性系中对应轨道描述所出现的差异说明, 当参数 µ 较大时, 非惯性坐标系中轨道描述会出现非椭圆情况, 因此在用椭圆摄动运动方程表示小天体运动时可能会出现一些定性分析上的错误注 :a) 轨道示意图 ;b) 能量 E 随时间变化, 单位时间为 G(m1 + m 2 )/(R 1 + R 2 ) 3,µ = 10 3, 初始轨道半长径 R 0 = 10 图 2 非惯性系下试验小天体轨道注 :µ = 0.4, R 0 = 10 图 3 非惯性系下试验小天体轨道

5 4 期沈欣和, 等 : 限制性三体问题中摄动运动方程的坐标系选择 409 注 : 红色点对应能量 E < 0 的部分, 即运动可以用椭圆轨道根数描述蓝色点对应能量 E > 0 部分, 即运动不可以用椭圆轨道根数描述图 4 非惯性系下试验小天体轨道 ( 图 3 的局部放大 ) 为了对该问题有一个全面和清晰的理解, 我们对其进行更详细的理论分析首先设定在惯性系中某一时刻, 做椭圆运动的小天体坐标和速度可以表示为 : X = R[cos(f + ω) cos i sin Ω sin(ω + f)] 其中, Y = R[sin(f + ω) + cos i cos Ω sin(ω + f)] Z = R sin i sin(ω + f) 1 Ẋ = [cos Ω sin(f + ω) + cos i sin Ω cos(ω + f)+ p e(cos Ω sin ω + cos i sin Ω cos ω)] 1 Ẏ = [sin Ω sin(f + ω) + cos i cos Ω cos(ω + f)+ p e(sin Ω sin ω + cos i cos Ω cos ω)] 1 Ż = [sin i sin ω sin f sin i cos ω(cos f + e)] p R =, (6) p (1 + e cos f), p = a(1 e2 ). (7) 变换到非惯性系中 ( 计算中仍使用惯性坐标系下的轨道根数 ), 可得 : { x = X X1, y = Y Y 1, z = Z ẋ = Ẋ X 1, ẏ = Ẏ Y 1, ż = Ż. (8) 由此可得能量为 : E = 1 2 υ2 (1 µ) r, (9)

6 410 天文学进展 36 卷其中, υ 2 = ẋ 2 + ẏ 2 + ż 2, (10) r = x 2 + y 2 + z 2. (11) 设 µ < 1 2, a > 1 µ, 舍去高阶小量项, 并令 f = 0, α = Ω t, 可以得到 : E = 1 2a + µ p µ(1 + e)[cos(ω + α) cos ω cos α(1 cos i)] + µ R, (12) 对于式中 i 取极大值, 可得到 : 令式 (13) 中 E 0, 并对其中 a 1 µ 2 [ E = 1 2a + µ R 1 p µ(1 + e) + µ R. (13) 取极大值, 则 R = a(1 e), 最后可以解得 : 3 + e 2µ 2 2 (1 µ)(1 + e) 1 e ]. (14) 上述关系式给出了在非惯性系中, 可以用椭圆表示的受摄二体问题对应的原惯性系中 a, e 和摄动体质量 µ 应满足的限制条件, 该条件与轨道倾角 i 无关当满足关系式时, 在以中心天体为原点的惯性坐标系下, 小天体的运动在特定时刻会变为双曲轨道, 不再适用椭圆型受摄方程此关系式是通过删掉 µ 的三次方及以上高阶小量后推导出来的, 图 5 表示当 µ 取值不同时, 判别式的误差情况图中红色区域表示运动可以用椭圆轨道根数描述, 蓝色区域表示运动不可以用椭圆轨道根数描述, 它们均为精确的轨道积分数值计算结果数值计算过程对 ω 和 f 在区间 [0, 2π],Ω 在区间 [0, π] 扫描取值, 在轨道积分中只要出现小天体的能量大于 0 的时刻, 便认为此运动不可用椭圆轨道根数描述而黑色曲线对应临界判别式结果, 曲线直接由判别式得到图 5 给出临界判别式给出的限制条件, 其可以很好地揭示在惯性系与非惯性系中对同一小天体的运动轨道描述上的差异, 并且通过对比不同轨道倾角, 进一步验证了轨道倾角对于临界判别式没有影响从图 5 中可以看出, 随着 µ 的增大, 原惯性系下的椭圆轨道转移到非惯性系中后, 能用椭圆轨道根数描述的运动区域在不断减小, 所以在非惯性系中, 当摄动天体质量增大时, 距离主天体较远处轨道均不能用稳定的椭圆来描述而在实际问题的研究中, [7] 此时若转而使用抛物线或双曲线描述轨道, 问题会变得十分复杂考虑由太阳 - 木星 - 小天体组成的限制性三体问题, 它对应的 µ 10 3, 由判别式可得, 在距离太阳 500 AU 范围内, 对稳定的运动, 这类小天体的轨道在两种坐标系下均可以使用椭圆轨道根数描述判别式同样可以表明, 对于大多数长周期彗星, 坐标系选择不重要, 但对于由双星 - 实验粒子组成的限制性三体问题, 如果 µ = 0.15, 在非惯性坐标系下可以用椭圆轨道根数的范围则缩小为小于 10 倍双星间距离所覆盖的区域, 此时不宜在以中心天体为坐标原点的非惯性系中建立摄动运动方程

7 4期沈欣和等限制性三体问题中摄动运动方程的坐标系选择 411 注 a) µ = 0.01 b) µ = 0.03 c) µ = 0.05 d) µ = 0.07 图为数值轨道积分与推导出的临界关系判别式比较结果红色区域表示运动可以用椭圆轨道根数描述蓝色区域表示运动不可以用椭圆轨道根数描述均为精确的惯性系中椭圆轨道根数转换到非惯性系中的数值计算结果图5 三体椭圆轨道半长径随偏心率的演化 3 雅可比坐标系下受摄二体问题摄动运动方程对于 N 体问题雅可比坐标系是近似于惯性的坐标系特别当 N = 3 时采用雅可比坐标系可以使得问题对应的哈密顿正则运动方程呈现很好的对称性根据三 [8] 体问题在惯性系下运动方程可以导出相应的雅可比坐标系下运动方程图 6 中位置矢量 r1 表示天体 C 与主天体 B 的相对位置位置矢量 r2 表示天体 D 和 B, C 的质心 O 的相对位置 mb 为主天体 B 的质量 mc 为天体 C 的质量 md 为天体 D 的质量约等于 0 设图6 雅可比坐标系下三体问题示意图

8 412 天文学进展 36 卷 r 2 = r 2 > r 1 = r 1, 则 C 和 D 的运动方程为 : r 1 = G(m B + m C )r 1 r Gm D(m B + m C ) m C R + r 1 + Gm D(m B + m C ) m B R r 1, (15) 其中, 或写为 : r 2 = Gm B(m B + m C + m D ) R + m B + m C r 2 + Gm C(m B + m C + m D ) m B + m C R 1 = m B = 1 [( ) n m B r 1 P n (cos ψ)] r 1 r 2 r 2 m n=0 B + m C r 2 m B + m C R + 1 = m C = 1 [( m ) n C r 1 P n (cos ψ)] r 1 + r 2 r 2 m n=0 B + m C r 2 m B + m C R r 2, (16), (17). (18) r 1 = G(m B + m C )r 1 r R 0 r 1, (19) r 2 = G(m B + m C + m D )r 2 r R 1 r 2, (20) R 0 = 1 r 2 R 1 = 1 r 2 n=2 n=2 [ ( ) n r1 Gm D α n P n (cos ψ)] r 2 [ ( ) n r1 Gm C β n P n (cos ψ)] r 2, (21), (22) 其中, α = m B n 1 + ( m C ) n 1 (m B + m C ) n 1, β = m B(m B + m C + m D ) (m B + m C ) 2 α n, n 2. (23) 由式 (19) 和式 (20) 可知, 右端第一项均对应二体问题运动对 C, 二体是由 B 和 C 构成, 摄动函数为 R 0 ; 而对 D, 二体却是由 B 和 C 的质心与 D 构成的, 摄动函数为 R 1 当 m C = 0, m D 0 时,m D 为摄动体 D 的质量,R 1 = 0, R 0 0, 运动体在摄动体轨道内绕中心天体运动, 称为外摄情形 ; 当 m D = 0, m C 0 时,m C 为摄动体的质量,R 0 = 0, R 1 0, 运动体在摄动体轨道外绕中心天体和摄动体的质心运动, 称为内摄情形选择这样的坐标系, 对限制性三体问题将始终给出完整且正确的运动描述对内摄和外摄情形,n 的定义不同当 m C = 0, m D 0 时,n 2 a 3 = Gm B ; 当 m D = 0, m C 0 时,n 2 a 3 = G(m B + m C ) 特别需要指出的是, 上述摄动函数的展开式, 无论是内摄还是外摄情况, 均不含 n = 1 项, 两种情况在形式上得以统一, 这在采用其他坐标系情况下难以达到而这种形式上的统一, 使得对于轨道运动的后续分析更加便捷

9 4 期沈欣和, 等 : 限制性三体问题中摄动运动方程的坐标系选择总结本文讨论的内容是在限制三体问题模型下对小天体轨道进行摄动理论分析过程中派生出来的进行轨道摄动分析首先要解决的问题就是如何选取轨道根数 ( 或积分常数 ), 即选择椭圆抛物线或是双曲线模型描述小天体轨道, 而模型选择是否合适简便, 与采用的坐标系直接相关尽管选择不同的坐标系不会造成小天体轨道本身动力学性质的改变, 但它将决定理论分析过程的复杂或简便程度, 不恰当的坐标系有时也可能导致错误的结论如本文描述的在惯性坐标系中的椭圆轨道, 转换到非惯性系中就变成了双曲线轨道这不仅带来了摄动运动方程的转换, 而且还面临着摄动函数对双曲根数的展开问题更有甚者, 当轨道从椭圆变化到双曲线时, 人们往往容易产生错觉, 认为运动本身变成不稳定了本文讨论的坐标系选择问题, 与小轨道偏心率和小轨道倾角导致的摄动运动方程奇点问题非常类似, 都不是运动本身产生的问题, 而是变量选取问题而本文详细讨论了惯性系和非惯性系中小天体轨道运动描述问题, 建立了两种坐标系下椭圆运动的判别关系式, 并经过数值计算, 验证了它的有效性为了解决坐标系选择造成的问题, 对于限制性三体问题内摄和外摄的情况, 本文建议采用雅可比坐标系从本文推导出的摄动运动方程可以看出, 选择雅可比坐标系对于第三体摄动可以建立形式统一的拉格朗日天体摄动运动方程 ( 质量因子有所不同 ), 在此基础上可以进一步建立无奇点摄动运动方程和哈密顿正则运动方程参考文献 : [1] Brouwer D, Clemence G M. Methods of Celestial Mechanics. New York: Academic Press, 1961: 465 [2] Murray C D. Dermott S F. Solar System Dynamics. Cambridge: Cambridge University Press, 1999:225 [3] Szebehely V. Theory of Orbits. London: Academic Press, 1967: 14 [4] Marchal C, Yoshida J, Sun Yi-Sui. Celestial mechanics, 1984, 34: 65 [5] 刘林. 人造地球卫星轨道力学. 北京 : 高等教育出版社, 1992: 91 [6] Kuznetsov V B. Solar System Research, 2016, 50: 224 [7] 刘林, 王歌. 深空探测器的双曲线轨道及其变化规律. 北京 : 中国科学出版社, 2002: 1127 [8] 孙义燧, 周济林. 现代天体力学导论. 北京 : 高等教育出版社, 2008: 29

10 414 天文学进展 36 卷 Selection of Coordinate System for Perturbed Equations of Motion in Restricted Three-body Problems SHEN Xin-he 1,2,3, WANG Wen-lei 1,2,3, XU Xue-qing 2,3, ZHOU Yong-hong 1,2,3, LIAO Xin-hao 2,3 (1.University of Chinese Academy of Sciences, Beijing , China; 2. Key Laboratory for Planetary Sciences of CAS, Nanjing , China; 3. Shanghai Astronomical Observatory,CAS, Shanghai , China) Abstract: This investigation discusses different expressions of motion orbit against the restricted three-body problems (RTBP), with two-selected coordinate system: (1) Coordinate system with central celestial body and mass center of perturbed body as the origin; (2) Noninertial coordinate system with central celestial body as the origin. As the orbital energy E of moving celestial body is positive or negative, we can determine the perturbed functions should be described by elliptical or hyperbolic orbit element. For convenience, this work give a critical relation discriminant with the semi-major axis and eccentricity. Results show that the stable elliptical orbit of inertial coordinate system, could not be described by the elliptical orbit element of non-inertial coordinate system, in case of big mass of the perturbed celestial body. So the elliptical orbit of inertial coordinate system should be represented by hyperbolic orbit element. Finally, this study find the up problem could be solved in the Jacobi coordinates, and a perturbation function expansion formula for any moving area is derived. Key words: restricted three-body problem; elliptical motion; Jacobi coordinate system; perturbation equations of motion

5 551 [3-].. [5]. [6]. [7].. API API. 1 [8-9]. [1]. W = W 1) y). x [11-12] D 2 2πR = 2z E + 2R arcsin D δ R z E = πr 1 + πr ) 2 arcsin

5 551 [3-].. [5]. [6]. [7].. API API. 1 [8-9]. [1]. W = W 1) y). x [11-12] D 2 2πR = 2z E + 2R arcsin D δ R z E = πr 1 + πr ) 2 arcsin 38 5 216 1 1),2) 163318) 163318). API. TE256 A doi 1.652/1-879-15-298 MODE OF CASING EXTERNA EXTRUSION BASED ON THE PRINCIPE OF VIRTUA WORK 1) ZHAO Wanchun,2) ZENG Jia WANG Tingting FENG Xiaohan School