CN108345754B - 一种彗尾数值仿真方法 - Google Patents

一种彗尾数值仿真方法 Download PDF

Info

Publication number
CN108345754B
CN108345754B CN201810182688.8A CN201810182688A CN108345754B CN 108345754 B CN108345754 B CN 108345754B CN 201810182688 A CN201810182688 A CN 201810182688A CN 108345754 B CN108345754 B CN 108345754B
Authority
CN
China
Prior art keywords
comet
dust
equation
comet dust
ordinary differential
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810182688.8A
Other languages
English (en)
Other versions
CN108345754A (zh
Inventor
王晓慧
毛李恒
岳宇贤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN201810182688.8A priority Critical patent/CN108345754B/zh
Publication of CN108345754A publication Critical patent/CN108345754A/zh
Application granted granted Critical
Publication of CN108345754B publication Critical patent/CN108345754B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)
  • Sampling And Sample Adjustment (AREA)

Abstract

本发明提供一种彗尾数值仿真方法，它包括以下几个步骤：步骤一：建立彗星尘埃尺寸分布模型；步骤二：建立彗星尘埃动力学模型；步骤三：彗星尘埃动力学方程转化；步骤四：动力学方程求解；通过以上步骤，能对彗尾中尘埃粒子的动力学进行精确数值仿真，确定彗尾的时间和空间分布，为人类了解彗尾的演化机制和彗星的动力学过程提供帮助，同时为彗星探测轨道设计提供重要参考。

Description

一种彗尾数值仿真方法
技术领域
本发明涉及一种彗尾数值仿真方法，可以对彗尾中尘埃粒子的动力学进行精确数值仿真，确定彗尾的时间和空间分布，属于航天技术领域。
背景技术
彗星回归时逐渐靠近太阳，在其处于太阳一定范围内时，受到太阳辐射加温的影响会在彗核周围因物质脱离而产生彗星尘埃，大量彗星尘埃逐步扩散到彗核附近区域时就形成了彗星独特的彗尾。由于每颗彗星的成分、轨道、彗星尘埃喷射速度都不同，因此彗星形成的彗尾也各不相同。目前对于彗尾的研究大都集中于观测数据分析，即通过地面或空间中的可见光、红外等设备对彗尾进行观测，然后对观测数据进行分析，确定彗尾的形状、成分等，目前还缺少对彗尾形状变化的预测模型的研究。因此，开展彗尾的数值仿真工作是十分迫切的，仿真结果不仅可以帮助解释彗尾的演化机制以及彗星的动力学过程，还可以对彗尾的分布范围进行精确预测，为彗星探测轨道设计提供重要参考。
发明内容
1、目的
本发明的目的是提供一种彗尾数值仿真方法，对彗尾中尘埃粒子的动力学进行精确数值仿真，确定彗尾的时间和空间分布，为人类了解彗尾的演化机制和彗星的动力学过程提供帮助，同时为彗星探测轨道设计提供重要参考。
2、技术方案
为实现上述发明目的，本发明采用以下技术方案。
本发明一种彗尾数值仿真方法，它包括以下四个步骤：
步骤一：建立彗星尘埃尺寸分布模型
所述彗星尘埃尺寸分布模型用于确定彗星尘埃尺寸满足的分布类型，然后根据待仿真的彗星尘埃总数，确定不同尺寸的彗星尘埃的个数；其建立的过程如下：
首先假设彗星尘埃尺寸满足的分布类型，一般假设为指数分布，即直径为a彗星尘埃的数目正比于a-p，其中p为常数，一般在3到4之间；然后根据计算机的计算性能以及数值仿真的精度，确定待仿真的彗星尘埃总数为N，利用前面假设的分布模型，即可求得各个区间范围内的彗星尘埃个数；
上述“建立彗星尘埃尺寸分布模型”，其建立过程的具体作法见“具体实施方式”内容；
步骤二：建立彗星尘埃动力学模型
所述建立彗星尘埃动力学模型是首先建立彗星尘埃逃离速度计算模型，该模型需通过地面观测数据进行参数拟合；然后建立彗星尘埃受力模型，分析彗星尘埃在空间中的受力状况；最后利用牛顿第二定律建立彗星尘埃动力学方程；
步骤三：彗星尘埃动力学方程转化
所述彗星尘埃动力学方程转化是指将步骤二中建立的二阶常微分动力学方程转化为一阶常微分方程组，为后续的数值求解奠定基础；其转化过程如下：
在步骤二中建立的动力学方程是关于彗星尘埃日心距D的二阶常微分方程，通过引入参数
Figure BDA0001589346320000021
可将二阶常微分方程转化为一阶常微分方程组，其中，
Figure BDA0001589346320000022
表示对彗星尘埃日心距D进行求导，v的物理含义是彗星尘埃在空间中的速度；
上述“彗星尘埃动力学方程转化”的具体作法见“具体实施方式”内容；
步骤四：动力学方程求解
所述动力学方程求解是指利用Runge-Kutta法对步骤三中建立的一阶常微分方程组进行数值求解，确定所有尘埃粒子在不同时刻的位置和速度；即首先利用Runge-Kutta法建立统一的迭代格式，然后利用计算机编程进行求解，确定所有尘埃粒子在不同时刻的位置和速度；
该Runge-Kutta法为通用的常微分方程数值解法。
上述“动力学方程求解”的具体作法见“具体实施方式”内容；
通过以上步骤，可以对彗尾中尘埃粒子的动力学进行精确数值仿真，确定彗尾的时间和空间分布，为人类了解彗尾的演化机制和彗星的动力学过程提供帮助，同时为彗星探测轨道设计提供重要参考。
3、优点及功效
本发明给出了一套完整的彗尾尘埃粒子动力学建模及数值求解方法，利用计算机编程可快速、准确求解彗尾的形成与演化过程，给出彗尾的时间和空间分布预报，对于人类了解彗尾的科学问题以及彗星探测轨道设计都有重要参考意义。
附图说明
图1为彗尾数值仿真方法流程图。
具体实施方式
以下结合附图对本发明作进一步描述。
由图1可以看出，本发明提供一种彗尾数值仿真方法，它包括以下四个步骤：
步骤一：建立彗星尘埃尺寸分布模型
所述建立彗星尘埃尺寸分布模型用于确定彗星尘埃尺寸满足的分布类型，然后根据待仿真的彗星尘埃总数，确定不同尺寸的彗星尘埃的个数；彗星尘埃尺寸一般假设为指数分布，即
n(a)＝ka-p (1)
式中，n(a)表示直径为a的彗星尘埃个数；k为比例系数，待定；p为常数，一般在3到4之间；
现假设待仿真的尘埃粒子总数为N，则可求得公式(1)中的比例系数为
Figure BDA0001589346320000031
式中，amin表示尘埃粒子的最小直径，amax表示尘埃粒子的最大直径。
根据公式(1)和公式(2)，可求得在任意区间[aL,aU](amin≤aL<aU≤amax)内的尘埃粒子数为
Figure BDA0001589346320000032
步骤二：建立彗星尘埃动力学模型
所述建立彗星尘埃动力学模型是指通过分析彗星尘埃在空间中的受力情况，利用牛顿第二定律建立彗星尘埃动力学方程。具体建模过程为：首先建立彗星尘埃逃离速度计算模型，该模型需通过地面观测数据进行参数拟合；然后建立彗星尘埃受力模型，分析彗星尘埃在空间中的受力状况；最后利用牛顿第二定律建立彗星尘埃动力学方程。相关建模方法已另申请发明专利“一种基于观测数据的彗星尘埃动力学建模方法”。
最终得到的彗星尘埃动力学方程为
Figure BDA0001589346320000041
式中，G为万有引力常量；MS为太阳的质量；m为彗星尘埃质量；D为彗星尘埃在日心黄道惯性坐标系下的位置向量；MC为彗核的质量；R为彗星质心在日心黄道惯性坐标系下的位置向量；S表示尘埃所处位置的能流密度；A表示彗星尘埃最大截面积；c表示光速；Qpr为辐射压力效率因子；vD为尘埃径向速度；
Figure BDA0001589346320000042
为单位向量，方向从太阳指向尘埃；
Figure BDA0001589346320000043
表示对彗星尘埃位置向量D求二阶导数，即为彗星尘埃加速度。
对公式(4)进行变形得
Figure BDA0001589346320000044
公式(5)即为待转化的二阶常微分动力学方程。
步骤三：彗星尘埃动力学方程转化
所述彗星尘埃动力学方程转化是指将步骤二中建立的二阶常微分动力学方程转化为一阶常微分方程组，为后续的数值求解奠定基础。
现引入参数
Figure BDA0001589346320000045
则公式(5)转化为
Figure BDA0001589346320000046
公式(6)为矢量形式的一阶常微分方程组，现将公式(6)中的向量都在日心惯性坐标系下进行投影，得到标量形式的六维一阶常微分方程组如公式(7)所示：
Figure BDA0001589346320000047
式中，带下标x，y，z的变量分别表示原向量在x，y，z轴的投影。公式(7)即为待求解的常微分方程组。
步骤四：动力学方程求解
所述动力学方程求解是指利用Runge-Kutta法对步骤三中建立的一阶常微分方程组进行数值求解，确定所有尘埃粒子在不同时刻的位置和速度。
所述Runge-Kutta法为通用的常微分方程数值解法。
公式(8)所示的是一般形式的常微分方程初值问题
Figure BDA0001589346320000051
式中，u表示函数自变量，t表示时间，t0表示初始时刻，T表示结束时刻。
对于该问题，经典的四级Runge-Kutta法的迭代格式为：
Figure BDA0001589346320000052
式中，tn表示当前时刻；un表示当前时刻自变量u的值；un+1表示下一时刻自变量u的值；h为积分步长；k1，k2，k3，k4均为迭代参数。
根据公式(9)所示的迭代格式，利用计算机编程就可对公式(7)所示彗星尘埃动力学微分方程组就行求解，得到所有尘埃粒子在任意时刻的位置和速度信息，从而确定彗尾的时间和空间分布。
运用以上方法，可以对彗尾中尘埃粒子的动力学进行精确数值仿真，确定彗尾的时间和空间分布，为人类了解彗尾的演化机制和彗星的动力学过程提供帮助，同时为彗星探测轨道设计提供重要参考。

Claims (5)

1.一种彗尾数值仿真方法，其特征在于：它包括以下几个步骤：
步骤一：建立彗星尘埃尺寸分布模型
所述彗星尘埃尺寸分布模型用于确定彗星尘埃尺寸满足的分布类型，然后根据待仿真的彗星尘埃总数，确定不同尺寸的彗星尘埃的个数；其建立的过程如下：
首先假设彗星尘埃尺寸满足的分布类型，假设为指数分布，即直径为a彗星尘埃的数目正比于a-p，其中p为常数，在3到4之间；然后根据计算机的计算性能以及数值仿真的精度，确定待仿真的彗星尘埃总数为N，利用前面假设的分布模型，即能求得各个区间范围内的彗星尘埃个数；
步骤二：建立彗星尘埃动力学模型
所述建立彗星尘埃动力学模型是首先建立彗星尘埃逃离速度计算模型，该模型需通过地面观测数据进行参数拟合；然后建立彗星尘埃受力模型，分析彗星尘埃在空间中的受力状况；最后利用牛顿第二定律建立彗星尘埃动力学方程；
步骤三：彗星尘埃动力学方程转化
所述彗星尘埃动力学方程转化是指将步骤二中建立的二阶常微分动力学方程转化为一阶常微分方程组，为后续的数值求解奠定基础；其转化过程如下：
在步骤二中建立的动力学方程是关于彗星尘埃日心距D的二阶常微分方程，通过引入参数
Figure FDA0002965778120000011
能将二阶常微分方程转化为一阶常微分方程组，其中，
Figure FDA0002965778120000012
表示对彗星尘埃日心距D进行求导，v的物理含义是彗星尘埃在空间中的速度；
步骤四：动力学方程求解
所述动力学方程求解是指利用Runge-Kutta法对步骤三中建立的一阶常微分方程组进行数值求解，确定所有尘埃粒子在不同时刻的位置和速度；即首先利用Runge-Kutta法建立统一的迭代格式，然后利用计算机编程进行求解，确定所有尘埃粒子在不同时刻的位置和速度；
该Runge-Kutta法为通用的常微分方程数值解法。
2.根据权利要求1所述的一种彗尾数值仿真方法，其特征在于：在步骤一中所述的“建立彗星尘埃尺寸分布模型”，其建立过程的具体作法如下：
所述建立彗星尘埃尺寸分布模型用于确定彗星尘埃尺寸满足的分布类型，然后根据待仿真的彗星尘埃总数，确定不同尺寸的彗星尘埃的个数；彗星尘埃尺寸假设为指数分布，即
n(a)＝ka-p·······················(1)
式中，n(a)表示直径为a的彗星尘埃个数；k为比例系数，待定；p为常数，在3到4之间；
现假设待仿真的尘埃粒子总数为N，则求得公式(1)中的比例系数为
Figure FDA0002965778120000021
式中，amin表示尘埃粒子的最小直径，amax表示尘埃粒子的最大直径；
根据公式(1)和公式(2)，求得在任意区间[aL,aU](amin≤aL<aU≤amax)内的尘埃粒子数为
Figure FDA0002965778120000022
3.根据权利要求1所述的一种彗尾数值仿真方法，其特征在于：在步骤二中所述的“建立彗星尘埃动力学模型”，其建立过程的具体作法如下：
所述建立彗星尘埃动力学模型是指通过分析彗星尘埃在空间中的受力情况，利用牛顿第二定律建立彗星尘埃动力学方程；具体建模过程为：首先建立彗星尘埃逃离速度计算模型，该模型需通过地面观测数据进行参数拟合；然后建立彗星尘埃受力模型，分析彗星尘埃在空间中的受力状况；最后利用牛顿第二定律建立彗星尘埃动力学方程；最终得到的彗星尘埃动力学方程为
Figure FDA0002965778120000023
式中，G为万有引力常量；MS为太阳的质量；m为彗星尘埃质量；D为彗星尘埃在日心黄道惯性坐标系下的位置向量；MC为彗核的质量；R为彗星质心在日心黄道惯性坐标系下的位置向量；S表示尘埃所处位置的能流密度；A表示彗星尘埃最大截面积；c表示光速；Qpr为辐射压力效率因子；vD为尘埃径向速度；
Figure FDA0002965778120000031
为单位向量，方向从太阳指向尘埃；
Figure FDA0002965778120000032
表示对彗星尘埃位置向量D求二阶导数，即为彗星尘埃加速度；
对公式(4)进行变形得
Figure FDA0002965778120000033
公式(5)即为待转化的二阶常微分动力学方程。
4.根据权利要求3所述的一种彗尾数值仿真方法，其特征在于：在步骤三中所述的“彗星尘埃动力学方程转化”，其转化的具体作法如下：
所述彗星尘埃动力学方程转化是指将步骤二中建立的二阶常微分动力学方程转化为一阶常微分方程组，为后续的数值求解奠定基础；
现引入参数
Figure FDA0002965778120000034
则上述公式(5)转化为
Figure FDA0002965778120000035
公式(6)为矢量形式的一阶常微分方程组，现将公式(6)中的向量都在日心惯性坐标系下进行投影，得到标量形式的六维一阶常微分方程组如公式(7)所示：
Figure FDA0002965778120000036
式中，带下标x，y，z的变量分别表示原向量在x，y，z轴的投影；公式(7)即为待求解的常微分方程组。
5.根据权利要求4所述的一种彗尾数值仿真方法，其特征在于：在步骤四中所述的“动力学方程求解”，其求解的具体作法如下：
所述动力学方程求解是指利用Runge-Kutta法对步骤三中建立的一阶常微分方程组进行数值求解，确定所有尘埃粒子在不同时刻的位置和速度；
所述Runge-Kutta法为通用的常微分方程数值解法；
下述公式(8)所示的是一般形式的常微分方程初值问题
Figure FDA0002965778120000041
式中，u表示函数自变量，t表示时间，t0表示初始时刻，T表示结束时刻；
对于该问题，经典的四级Runge-Kutta法的迭代格式为：
Figure FDA0002965778120000042
式中，tn表示当前时刻；un表示当前时刻自变量u的值；un+1表示下一时刻自变量u的值；h为积分步长；k1，k2，k3，k4均为迭代参数；
根据公式(9)所示的迭代格式，利用计算机编程就对公式(7)所示彗星尘埃动力学微分方程组就行求解，得到所有尘埃粒子在任意时刻的位置和速度信息，从而确定彗尾的时间和空间分布。
CN201810182688.8A 2018-03-06 2018-03-06 一种彗尾数值仿真方法 Active CN108345754B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810182688.8A CN108345754B (zh) 2018-03-06 2018-03-06 一种彗尾数值仿真方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810182688.8A CN108345754B (zh) 2018-03-06 2018-03-06 一种彗尾数值仿真方法

Publications (2)

Publication Number Publication Date
CN108345754A CN108345754A (zh) 2018-07-31
CN108345754B true CN108345754B (zh) 2021-05-25

Family

ID=62957838

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810182688.8A Active CN108345754B (zh) 2018-03-06 2018-03-06 一种彗尾数值仿真方法

Country Status (1)

Country Link
CN (1) CN108345754B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110955955B (zh) * 2019-11-05 2022-11-29 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种基于对数拟合的尘埃等离子体数值计算提速方法
CN114154309A (zh) * 2021-11-17 2022-03-08 北京航空航天大学 一种考虑冰含量的彗星准束缚态尘埃动力学建模方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
SU1165950A1 (ru) * 1983-11-23 1985-07-07 Osherov Ruslan S Нефелометр Ошерова
CN105760575A (zh) * 2016-01-17 2016-07-13 中国海洋大学 渤海海上溢油输移、扩展数值预报系统的建立方法
CN105865459A (zh) * 2016-03-31 2016-08-17 北京理工大学 一种考虑视线角约束的小天体接近段制导方法
CN105890577A (zh) * 2015-01-23 2016-08-24 北京空间飞行器总体设计部 一种适用于深空探测器在轨多个天体合影成像方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103419947A (zh) * 2013-08-21 2013-12-04 北京理工大学 微重力环境下自主着陆导航控制地面试验验证系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
SU1165950A1 (ru) * 1983-11-23 1985-07-07 Osherov Ruslan S Нефелометр Ошерова
CN105890577A (zh) * 2015-01-23 2016-08-24 北京空间飞行器总体设计部 一种适用于深空探测器在轨多个天体合影成像方法
CN105760575A (zh) * 2016-01-17 2016-07-13 中国海洋大学 渤海海上溢油输移、扩展数值预报系统的建立方法
CN105865459A (zh) * 2016-03-31 2016-08-17 北京理工大学 一种考虑视线角约束的小天体接近段制导方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Influence of the magnetic field on the density distribution of solar wind protons and cometary ions in the shock layer ahead of cometary ionospheres;Baranov, V. B 等;《ASTRONOMY LETTERS-A JOURNAL OF ASTRONOMY AND SPACE ASTROPHYSICS》;20170430;第43卷(第2期);126-133 *
Model of the dust tail of comet C/2011 L4 (PANSTARRS);S. V. Kharchuk 等;《Dynamics and Physics of Bodies of the Solar System》;20150101;第31卷;34-36 *
尘埃等离子体动力学及电磁特性研究;李辉;《中国博士学位论文全文数据库 基础科学辑》;20180115(第1期);A005-78 *

Also Published As

Publication number Publication date
CN108345754A (zh) 2018-07-31

Similar Documents

Publication Publication Date Title
Leung et al. Research needs and directions of regional climate modeling using WRF and CCSM
Žic et al. Heliospheric propagation of coronal mass ejections: Drag-based model fitting
CN108090621B (zh) 一种基于分阶段整体优化的短期风速预测方法与系统
CN106646645B (zh) 一种重力正演加速方法
CN105224715A (zh) 一种山区地貌下强风三维脉动风场综合模拟方法
CN108345754B (zh) 一种彗尾数值仿真方法
CN106681343B (zh) 一种航天器姿态跟踪低复杂度预设性能控制方法
Pardyjak et al. QUIC-URB v. 1.1: Theory and User’s Guide
CN113467241B (zh) 凸曲率着陆轨迹燃耗优化方法
CN108627667A (zh) 基于光度序列同时估计空间失稳目标进动和自旋速率方法
CN108181806A (zh) 基于采样输出的空间机器人位置与姿态自抗扰控制方法
CN114580310A (zh) 一种基于palm实现wrf模拟风场降尺度处理的方法
US20120150499A1 (en) Invertible contact model
McCann et al. Rigid body pose estimation on TSE (3) for spacecraft with unknown moments of inertia
Pfaff et al. The Spherical Grid Filter for Nonlinear Estimation on the Unit Sphere
CN109325288A (zh) 一种基于不确定性优化的固体运载器总体参数确定方法及系统
CN110231619B (zh) 基于恩克法的雷达交接时刻预报方法及装置
Waibel et al. Physics meets machine learning: Coupling FFD with regression models for wind pressure prediction on high-rise facades
CN111985050B (zh) 弱引力小天体表面着陆器安全投放窗口生成方法
CN108446466B (zh) 一种基于观测数据的彗星尘埃动力学建模方法
CN104657595B (zh) 一种单颗粒曳力模型系数标定方法
CN109063268B (zh) 人群疏散宏观模型迭代式仿真方法
CN110597088B (zh) 一种极坐标系下的车辆动力学仿真方法
CN107291986B (zh) 一种飞行器最优传感器选择方法
Mortezazadeh et al. Modelling urban airflows by a new parallel high-order semi-Lagrangian 3D fluid flow solver

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant