形状不规则双小行星系统相互作用高精度高效率建模方法与流程

1.本发明涉及一种双小行星系统相互作用建模方法，尤其适用于涉及双小行星系统演化预测的双小行星系统相互作用高精度高效率建模方法，属于航空航天技术领域。


背景技术：

2.小行星是指太阳系内类似行星环绕太阳运动，但是体积和质量比行星小得多，无法清除轨道附近区域的天体。据估计，太阳系内有超过127万颗小行星，其中约15％是双小行星系统。如何计算双小行星系统内部的相互作用力和力矩，是研究双小行星系统演化的关键技术。双小行星系统内部的相互作用是目前科研人员关注的热点问题之一。
3.在已发展的研究双小行星系统内部相互作用力和力矩的方法中，在先技术[1](shi,yu,yue wang,and shijie xu."mutual gravitational potential,force,and torque of a homogeneous polyhedron and an extended body:an application to binary asteroids."celestial mechanics and dynamical astronomy 129.3(2017):307-320.)，将相互势考虑为一个小行星在另一个小行星引力场中的质量积分(即次小行星在主小行星引力场中的质量积分)。将相互势的表达式在次小行星的质心位置进行泰勒展开，利用主小行星引力势能函数的有限阶偏导数、次小行星的惯性积分和两个小行星之间的相对位置、姿态表示双小行星系统的相互势。该方法的缺点是当两个小行星的距离较近时，泰勒级数的截断误差将显著增大，且无法完全描述小行星不均匀的质量分布带来的影响。
[0004]
在先技术[2](hou,xiyun,daniel j.scheeres,and xiaosheng xin."mutual potential between two rigid bodies with arbitrary shapes and mass distributions."celestial mechanics and dynamical astronomy 127.3(2017):369-395.)，根据相互势的基本定义计算双小行星系统的相互势。将两个小行星质量微元之间的距离项按照勒让德级数展开，使用广义惯性积分来描述小行星不规则几何构型带来的质量分布不均。该方法的缺点是受限于勒让德级数的项数，存在截断误差，且无法完全描述小行星内部质量不均匀分布所带来的的影响。
[0005]
在先技术[3](yu,yang,et al."a finite element method for computational full two-body problem:i.the mutual potential and derivatives over bilinear tetrahedron elements."celestial mechanics and dynamical astronomy 131.11(2019):1-21.)，根据相互势、力、力矩的定义计算两个小行星之间的相互作用力、力矩。基于有限元划分方法将两个小行星建立为有限元模型，并根据双线性形状函数插值的方法，给出了相互势、力、力矩的数值计算形式。该方法的缺点是有限元方法需要大量的域元才能够保证精度，使得该方法在实际操作中具有较大的计算负担。


技术实现要素：

[0006]
本发明公开的形状不规则双小行星系统相互作用高精度高效率建模方法要解决
的技术问题是：提供一种双小行星系统内部相互作用力和力矩的高精度高效率建模方法，基于复合引力场模型并结合有限元方法实现相互作用的快速建模，描述两个小行星不均匀的质量分布，且能够控制结果的截断误差，保证预测精度。所述双小行星系统相互作用指双小行星系统内部相互作用力和力矩。
[0007]
本发明的目的是通过下述技术方案实现的：
[0008]
本发明公开的形状不规则双小行星系统相互作用高精度高效率建模方法，建立主小行星的复合引力场模型，通过数值插值方法，提高双小行星系统相互作用预测效率；建立次小行星的有限元模型，计算域元节点对相应域元的插值系数，域元划分的精度能够有效控制截断误差；通过复合引力场模型，确定域元节点在主小行星引力场中的引力场参数并保存；基于插值系数和前述引力场参数，确定域元受到主小行星的作用，有限元模型在预测相互作用的过程中能够充分利用主小行星引力场，提高计算效率，实现两个小行星之间相互作用力和力矩的快速建模；根据当前双小行星系统的状态进行相互作用快速建模，将建模结果代入运动学方程，获得下一时刻系统的状态，从而预测双小行星系统演化的过程。
[0009]
本发明公开的形状不规则双小行星系统相互作用高精度高效率建模方法，包括如下步骤：
[0010]
步骤一，建立主小行星的复合引力场模型，复合引力场模型通过标称引力场和残差引力场的叠加描述真实引力场。因为残差引力场使用数值插值方法，所以复合引力场模型具有较高的数值计算效率，且能够体现主小行星质量分布不均匀对引力场的影响，提高预测精度。通过建立复合引力场描述主小行星的引力场，可降低计算负担，提高双小行星系统相互作用建模效率。
[0011]
复合引力场模型能够高效地描述单颗小行星的引力场，当小行星的质量分布不均匀时，复合引力场仍是有效的。复合引力场模型通过标称引力场vn和残差引力场t的叠加描述真实引力场v
[0012]
v＝vn t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0013]
标称引力场vn是一个三轴椭球引力场，用于描述单颗小行星引力场的主要特征。
[0014]
复合引力场通过楔形网格的分层堆叠有效地覆盖小行星附近的空间区域。每一个楔形网格包含六个顶点，第k层网格的第i个顶点的位置矢量与描述小天体几何形状的多面体三角形表面顶点位置矢量ai的关系为
[0015][0016]
其中，θ是整个多面体的表面的平均固体角。
[0017]
在建立复合引力场的楔形网格结构后，需要计算在每一个网格节点上真实引力场v与标称引力场vn的残差值ω，即引力势能函数差u∈r和引力势能梯度差g＝[g
x
,gy,gz]∈r1×3。定义描述主小行星形状的多面体模型共有n
facet
个表面，则当楔形网格共有n
layers
层时，共包含n
laypers
×nfacet
个网格，若每一个楔形网格为n
node
点楔形单元，则共包含n
laypers
×nfacet
×nnode
个节点。将所有的残差值保存为n
laypers
×nfacet
×nnode
行4列的矩阵，所述矩阵的每一行即为节点上的残差值ω的行数据，即ω＝[u,g
x
,gy,gz]。
[0018]
残差引力场t通过数值插值方法表示，能够有效地确定楔形网格单元局部区域内的残差值
[0019][0020]
其中，系数ωi和pi分别是该楔形网格单元内第i个插值节点上的真实残差值和基函数，l是指定场点在当前楔形网格内的自然坐标，n是当前楔形单元。复合引力场模型包括残差引力场，因此，复合引力场模型能够快速计算出真实引力场由于小行星质量分布不均匀带来的引力场变化。
[0021]
为了进一步提高预测效率，作为优选，采用二阶二次球谐函数描述所述三轴椭球的引力场
[0022][0023]
其中，μ是椭球体的引力常数，是场点与椭球中心的距离，c
20
和c
22
是球谐系数，r＝[x,y,z]
t
是场点的位置矢量。
[0024]
步骤二，建立次小行星的有限元模型。将次小行星划分为有限个域元，并计算域元节点对相应域元的插值系数。根据域元划分的精度，能够有效控制双小行星系统相互作用随小天体之间距离减小而增大的截断误差。
[0025]
预设次小行星的有限元划分精度，利用有限元划分算法将目标小行星离散为有限个域元的集合，定义有限元模型包含n
fem
个域元，n
femvert
个域元节点。在每一个域元内，计算对于当前域元，域元节点的插值系数
[0026]
κ＝f
κ
(x,y,z)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0027]
其中，f
κ
是域元类型相应的插值系数计算方法，x,y,z是当前域元中的域元顶点位置参数。
[0028]
作为优选，设定次小行星的有限元划分精度，利用有限元划分算法将目标多面体离散为有限个四面体域元的集合，在每一个四面体域元内，计算对于当前域元，域元顶点的插值系数
[0029][0030]
其中，xi,yi,zi,i＝1,2,3,4是四面体域元的第i个节点的位置坐标。
[0031]
步骤三，根据步骤一中所述复合引力场模型，确定每一个域元节点在主小行星引力场中的引力场参数，所述引力场参数包括引力势能、引力势能梯度值。并将前述引力势能、引力势能梯度值保存，为双小行星系统相互作用的建模作准备，从而避免域元节点引力场参数的重复计算，提高双小行星系统相互作用建模效率。
[0032]
定义任一域元节点在主小行星固连坐标系中的位置矢量为r，则该域元节点在复合引力场模型中的引力势能为
[0033][0034]
引力势能梯度为
[0035][0036]
其中，un(r)是标称引力场的引力势能，ui是该域元节点所在楔形单元的第i个插值节点上引力势能的残差值，pi(r)是该楔形单元的基函数，是标称引力场的引力势能梯度，gi＝[g
xi
,g
yi
,g
zi
]
t
是该有限元节点所在楔形单元的第i个插值节点上引力势能梯度的残差值，n是该域元节点所在的楔形单元。
[0037]
将每一个域元节点上的引力势能值u(r)∈r，引力势能梯度值保存为n
femvert
行4列的矩阵，所述矩阵的每一行即为域元节点上的引力场参数的行数据，即f＝[f
p
,f
x
,fy,fz]。将前述引力势能、引力势能梯度值保存，为双小行星系统相互作用的建模作准备，从而避免域元节点引力场参数的重复计算。
[0038]
步骤四，基于步骤二得到的插值系数和步骤三保存的引力场参数，确定每一个次小行星域元受到主小行星的作用，每一个域元受到的作用都能够反映次小行星在该域元处不均匀的密度分布，从而提高双小行星系统相互作用预测精度。通过建立次小行星的有限元模型，在计算相互作用力和力矩的过程中能够充分地利用主小行星引力场，提高计算效率。综合考虑全部次小行星域元受到主小行星的作用，实现两个小行星之间相互作用力和力矩的快速建模。
[0039]
双小行星系统的相互作用包括双小行星系统中相互势p、主小行星对次小行星的力f、主小行星对次小行星的力矩μ1和次小行星对主小行星的力矩μ2。根据步骤二建立的有限元模型，两个小行星之间的相互作用为每一个次小行星域元受到主小行星的作用之和
[0040][0041][0042][0043]
μ2＝-r
×
f-μ1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0044]
其中，是有限元域元，是次小行星有限元模型有限元域元集合，d是次小行星域元节点的位置矢量在次小行星固连坐标系中的投影，κ是域元节点在当前域元内的插值系数，σ是次小行星在有限元域元节点处的密度值，第一个求和符号是指在有限元域元集合中对所有有限元域元求和，第二个求和符号是指在单个域元中对所有域元节点求和，第三个求和符号是指在单个楔形网格单元内对所有插值节点求和。
[0045]
由于在步骤三中已经计算过节点上的引力场参数值，在步骤四中仅需要根据不同的次小行星域元节点密度值计算节点参数。复合引力场所保存的引力场参数能够描述主小行星的质量不均匀分布，步骤三中域元节点参数值能够反映次小行星的质量不均匀分布，在此基础上，能够描述双小行星系统中主小行星和次小行星的质量不均匀分布，进而提高双小行星系统相互作用力和力矩的建模精度。
[0046]
遍历次小天体所有的域元并求和，即得到双小行星系统中相互势p、主小行星对次小行星的力f、主小行星对次小行星的力矩μ1和次小行星对主小行星的力矩μ2，从而实现双小行星系统相互作用的快速预测。
[0047]
步骤五，根据当前时刻双小行星系统的状态代入步骤四进行双小行星系统相互作用力和力矩的快速建模及预测，将步骤四得到的双小行星系统相互作用力和力矩建模结果代入双小行星系统相对运动微分方程，获得下一时刻双小行星系统的状态。
[0048]
双小行星系统相对运动微分方程为
[0049][0050]
其中，p是系统的总线动量，r是由小天体2质心指向小天体1质心的位置矢量，γi是小天体的角动量，ωi是小天体的角速度，a是由小天体2固连坐标系转至小天体1固连坐标系的姿态矩阵，a2是由小天体2固连坐标系转至惯性坐标系的姿态矩阵。f是小天体2对小天体1的力，μ1(μ2)是小天体2(小天体1)对小天体1(小天体2)的力矩。
[0051]
步骤六，重复步骤五，预测双小行星系统演化的过程，根据双小行星系统演化结果解决双小行星系统相关工程技术问题。
[0052]
有益效果：
[0053]
1、本发明公开的形状不规则双小行星系统相互作用高精度高效率建模方法，通过建立主小行星的复合引力场模型，利用残差引力场的数值插值方法，提高双小行星系统相互作用建模效率。
[0054]
2、本发明公开的形状不规则双小行星系统相互作用高精度高效率建模方法，通过建立次小行星的有限元模型，在计算相互作用力和力矩的过程中充分地利用主小行星的引力场，提高计算效率。
[0055]
3、本发明公开的形状不规则双小行星系统相互作用高精度高效率建模方法，基于复合引力场模型和有限元模型，能够同时描述双小行星系统中主小行星和次小行星的质量不均匀分布，提高双小行星系统相互作用预测精度。
[0056]
4、本发明公开的形状不规则双小行星系统相互作用高精度高效率建模方法，结合有限元方法能够有效控制双小行星系统相互作用随小天体之间距离减小而增大的截断误差，进而提高双小行星系统相互作用建模及预测精度。
附图说明
[0057]
图1是本发明的形状不规则双小行星系统相互作用高精度高效率建模方法流程图；
[0058]
图2是双小行星系统相互作用示意图；
3.8524
×
10-9
km2/s2,引力加速度为[-1.4043
×
10-8
,1.6442
×
10-9
,1.4062
×
10-9
]km/s2。复合引力场计算所得的引力势能为-1.6157
×
10-7
km2/s2，引力加速度为[-1.7325
×
10-7
,8.66598767621932
×
10-10
,2.2087
×
10-9
]km/s2。
[0070]
步骤四，基于步骤二得到的插值系数和步骤三保存的引力场参数，确定每一个次小行星域元受到主小行星的作用，每一个域元受到的作用都能够反映次小行星在该域元处不均匀的密度分布，从而提高双小行星系统相互作用建模精度。通过建立次小行星的有限元模型，在计算相互作用力和力矩的过程中能够充分地利用主小行星引力场，提高计算效率。综合考虑全部次小行星域元受到主小行星的作用，实现两个小行星之间相互作用力和力矩的快速建模。
[0071]
若主小行星相对于惯性空间的姿态角为(27.04
°
，10
°
，-83.93
°
)，从星相对于惯性空间的姿态角为(0
°
，0
°
,180
°
)，从星相对于主星的距离为2.5659km，则相互势为-8176kgkm2s-2
，相互力为3208kn，主星对从星的力矩1.8121knkm，从星对主星的力矩为22.3937knkm。
[0072]
步骤五，根据当前时刻双小行星系统的状态代入步骤四进行双小行星系统相互作用力和力矩的快速预测，将步骤四得到的双小行星系统相互作用力和力矩预测结果代入双小行星系统相对运动微分方程，获得下一时刻双小行星系统的状态。
[0073]
以步骤四双小行星系统的状态为例，系统的总线动量p＝[-2.6933,-1.7503,-0.0590]
t
×
107kgkm/s，由小天体2质心指向小天体1质心的位置矢量r＝[-1.3840,2.1147,0.4431]
t
km，小天体的角动量γ1＝[0.2969,-0.5818,3.7043]
×
105kgkm2/s，γ2＝[0,0，2.9123]
×
108kgkm2/s，小天体的角速度ω1＝[0,0,0.0001]rad/s，ω2＝[0,0,0.0006]rad/s，由小天体2固连坐标系转至小天体1固连坐标系的姿态矩阵由小天体2固连坐标系转至惯性坐标系的姿态矩阵
[0074]
步骤六，重复步骤五，预测双小行星系统演化的过程。
[0075]
经过200h后，系统的总线动量p＝[3.1103,0.7523,0.2859]
t
×
107kgkm/s，由小天体2质心指向小天体1质心的位置矢量r＝[0.5457,-2.4787,0.3788]
t
km，小天体的角动量γ1＝[-0.4878,0.4835,3.7297]
×
105kgkm2/s，γ2＝[-0.0031,-0.0061，2.9122]
×
108kgkm2/s，小天体的角速度ω1＝[0.0005,-0.0018,0.1010]rad/s，ω2＝[-0.0008,-0.0015,0.6313]rad/s，由小天体2固连坐标系转至小天体1固连坐标系的姿态矩阵由小天体2固连坐标系转至惯性坐标系的姿态矩阵
[0076]
bia-1999kw4双小行星系统中，次小行星相对于主小行星轨道运动的半长轴随时间的变化如附图3所示。
[0077]
bia-1999kw4双小行星系统中，次小行星相对于主小行星轨道运动的偏心率随时间的变化如附图4所示。
[0078]
bia-1999kw4双小行星系统中能量转化随时间的变化图如附图5所示。
[0079]
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。