CN116306185A - 基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统 - Google Patents
基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统 Download PDFInfo
- Publication number
- CN116306185A CN116306185A CN202310159735.8A CN202310159735A CN116306185A CN 116306185 A CN116306185 A CN 116306185A CN 202310159735 A CN202310159735 A CN 202310159735A CN 116306185 A CN116306185 A CN 116306185A
- Authority
- CN
- China
- Prior art keywords
- model
- discrete element
- irregular discrete
- meteorite
- representing
- 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.)
- Granted
Links
- 230000001788 irregular Effects 0.000 title claims abstract description 225
- 238000000034 method Methods 0.000 title claims abstract description 56
- 230000033001 locomotion Effects 0.000 title claims abstract description 42
- 239000002245 particle Substances 0.000 claims abstract description 211
- 230000008859 change Effects 0.000 claims abstract description 47
- 238000004880 explosion Methods 0.000 claims abstract description 39
- 230000008569 process Effects 0.000 claims abstract description 24
- 238000004364 calculation method Methods 0.000 claims abstract description 10
- 206010034719 Personality change Diseases 0.000 claims abstract description 9
- 238000005474 detonation Methods 0.000 claims description 37
- 230000009471 action Effects 0.000 claims description 34
- 238000004088 simulation Methods 0.000 claims description 29
- 238000002679 ablation Methods 0.000 claims description 26
- 239000013598 vector Substances 0.000 claims description 18
- 238000013016 damping Methods 0.000 claims description 12
- 230000007704 transition Effects 0.000 claims description 9
- 230000001133 acceleration Effects 0.000 claims description 7
- 238000011084 recovery Methods 0.000 claims description 7
- 238000010276 construction Methods 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 108050006002 RNA polymerase sigma factor FliA Proteins 0.000 claims description 3
- 230000002146 bilateral effect Effects 0.000 claims description 3
- 230000004069 differentiation Effects 0.000 claims description 3
- 230000009545 invasion Effects 0.000 claims description 3
- 230000035939 shock Effects 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000012544 monitoring process Methods 0.000 abstract description 4
- 239000012634 fragment Substances 0.000 description 8
- 230000000694 effects Effects 0.000 description 6
- 230000003993 interaction Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 235000017899 Spathodea campanulata Nutrition 0.000 description 3
- 230000008878 coupling Effects 0.000 description 3
- 238000010168 coupling process Methods 0.000 description 3
- 238000005859 coupling reaction Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 239000000126 substance Substances 0.000 description 3
- 230000005457 Black-body radiation Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000011439 discrete element method Methods 0.000 description 2
- 239000002360 explosive Substances 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000002093 peripheral effect Effects 0.000 description 2
- 230000005855 radiation Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 230000008901 benefit Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000001808 coupling effect Effects 0.000 description 1
- 230000007123 defense Effects 0.000 description 1
- 230000008021 deposition Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000000375 direct analysis in real time Methods 0.000 description 1
- 238000012063 dual-affinity re-targeting Methods 0.000 description 1
- 239000000428 dust Substances 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000010438 heat treatment Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000006855 networking Effects 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 239000012798 spherical particle Substances 0.000 description 1
- 239000004575 stone Substances 0.000 description 1
- 238000013517 stratification Methods 0.000 description 1
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
- G06F30/25—Design optimisation, verification or simulation using particle-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种基于不规则离散元的陨星进入地球大气层后运动仿真方法，其包括：借助不规则离散元颗粒构建陨星模型，计算两两相邻不规则离散元颗粒间的接触力，计算陨星模型的形变；计算不规则离散元颗粒的轨迹，计算陨星模型的角速度和姿态变化，求解陨星模型的质量损失、温度变化和形状变化，计算陨星模型的结构变化；判断陨星模型在高度为0前的质量损失是否与质量相等，若是，则陨星模型的运动参数解算停止；否则，获得陨星模型的经纬度，当高度为0时确定陨星模型的落点范围。本发明利用可变形不规则离散元模型，解算陨星进入地球大气层后的运动演化过程，更符合陨星内部真实结构特征和爆炸扩张原理，对陨星监测和防御有现实意义。
Description
技术领域
本发明属于星际空间的颗粒分析技术领域，特别是一种基于不规则离散元的陨星进入地球大气层后运动仿真方法及其仿真系统。
背景技术
陨星(又称流星)是指来自星际空间的颗粒(如小行星或星体碎片)以超高速进入地球大气层后在剧烈气动作用下的产物。陨星进入大气层后会以12km/s～20km/s的极高速度下落，并发生烧蚀、解体和空爆等耦合事件。未被完全烧蚀的陨星着陆后形成陨石。
历史上对陨星的研究有两波高潮，首先是20世纪末至21世纪初期，美国、加拿大和欧洲等国建立了陨星的观测组网，希望在地球表面约100万平方公里范围的天空内实现陨星和流星观测全覆盖。陨星观测组网通过CCD相机拍摄陨星进入时的火球和明亮轨迹，进而确定撞击区域，并对陨石进行回收。这些组网包括美国中部的“草原网络”(PrairieNetwork,PN)、加拿大的“陨石观测与回收计划”(Meteorite Observation and RecoveryProgram,MORP)和在德国、捷克及西班牙等地设置的“欧洲火球网络”(European FireballNetwork,EN)。基于所设立的观测组网，研究者们收集了部分陨星进入地球大气层时的轨迹和速度变化数据，例如捷克的Pribram陨星、美国的Lost City陨星和加拿大的Innisfree陨星。2013年，在俄罗斯车里雅宾斯克(Chelyabinsk)上空，一枚直径20m的陨星以约20km/s的速度进入地球大气层，并在距离地面约30km高度处发生了剧烈空爆。爆炸当量可达500万吨TNT，造成1500多人受伤，多处房屋受损。这是继1908年通古斯陨星大爆炸(爆炸当量为1000～4000万吨TNT)以来最严重的一次陨星爆炸事故。Chelyabinsk陨星爆炸重新引发了国际上对大体积陨星进入事件的关注，开启了第二轮针对陨星探测与防御的研究热潮。美国NASA陆续建立和开展了“近地天体观测计划”(Near-Earth Object ObservationsProgram)、“近地天体研究中心”(Center for Near-Earth Object Studies,CNEOS)和“双小行星重定向测试任务”(Double Asteroid Redirection Test,DART)。2013年国际上建立了“国际小行星预警网”(International Asteroid Warning Network,IAWN)，各成员国之间共享可能有威胁的小行星信息。2018年我国正式加入该组织，南京紫金山天文台的近地天体望远镜是我国目前主要的观测设备。虽然我国对近地小行星和陨星观测的起步晚，但发展迅速。
与各向同性、结构规则的人造航天器不同，陨星大多来源于小行星或小行星母体碰撞解体后产生的碎片，而小行星多是由不规则物质堆积而成的碎石堆结构，因此陨星在初始进入地球大气层时通常有较大的孔隙率，而非独立的完整结构。经统计，陨星的质量分布范围广(从克到百万公斤)、进入速度极快、外形各异、结构分布极不均匀，更易发生烧蚀。受陨星的复杂结构和超高速进入条件限制，目前对陨星进入过程的运动规律认识较浅，特别是陨星在超高速空气动力学作用下的解体和爆炸过程，其中代表性研究结果有美国NASA艾姆斯研究中心(Ames Research Center)Lorien F.Wheeler等人提出的模拟陨星解体的碎片云(Fragment-cloud model,FCM)模型。该模型认为陨星在解体时会形成大量由小碎片和灰尘液滴组成的碎片云结构，并通过调整初始参数可与描述陨星动能损失的能量沉积曲线吻合。但目前的研究多将陨星简化为单独的球形物质，无法准确模拟陨星的不规则外形、内部结构的相互作用、爆炸成因以及横向碎片的弹出。通过借助不规则离散元方法可以仿真陨星的内部结构和不规则外形特征，不规则离散元方法是利用不规则外形的颗粒来仿真颗粒系统的动力学特性，从而解决球形离散元无法准确描述颗粒外形以及非球形颗粒间可能出现的结构互锁等现象的问题。然而，在陨星进入地球大气的过程中，由于超高的进入速度，气动力和气动热作用十分剧烈，且伴随质量烧蚀、解体和爆炸等一系列复杂过程。因此，为解决上述两个关键问题，寻求一种基于不规则离散元的陨星进入地球大气层后运动仿真方法，以有效利用不规则离散元模型模拟陨星烧蚀边界的后缩和空爆现象，进而解算陨星的轨迹、姿态和内部结构变化，并模拟陨星解体后碎片的弹出是十分迫切且必要的。
发明内容
本发明针对上述现有技术中的缺陷，提出一种基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统。该方法包括借助不规则离散元颗粒构建陨星模型，计算两两相邻不规则离散元颗粒间的接触力，计算陨星模型的形变；计算不规则离散元颗粒的轨迹，计算陨星模型的角速度和姿态变化，求解陨星模型的质量损失、温度变化和形状变化，计算陨星模型的结构变化；判断陨星模型在高度为0前的质量损失是否与质量相等，若是，则陨星模型的运动参数解算停止；否则，获得陨星模型的经纬度，当高度为0时确定陨星模型的落点范围。本发明利用可变形不规则离散元模型，解算陨星进入地球大气层后的运动演化过程，更符合陨星内部真实结构特征和爆炸扩张原理，对陨星监测和防御有现实意义。
本发明提供一种基于不规则离散元的陨星进入地球大气层后运动仿真方法，其包括以下步骤：
S1、借助不规则离散元颗粒构建陨星模型并初始化：在不规则离散元程序中，采用Delaunay三角网格构建不规则离散元颗粒，若干所述不规则离散元颗粒组合而成陨星模型，对所有所述不规则离散元颗粒赋予不同的密度、体积、外形和力学性质参数；
S2、计算两两相邻所述不规则离散元颗粒之间的接触力，并计算所述陨星模型的形变；
S3、计算所述不规则离散元颗粒的轨迹，并计算所述陨星模型的角速度和姿态变化；
S31、建立地心固连系O-XYZ和东北天坐标系o-xyz；
S32、计算所述不规则离散元颗粒的轨迹方程；
S33、计算所述陨星模型的角速度和姿态变化；
S4、求解所述陨星模型的质量损失、温度变化和形状变化；
S5、计算所述陨星模型的结构变化，即模拟所述陨星模型的解体和空爆过程；
S51、基于陨星强度，获得陨星在气动作用下的解体判据；
ρairV2＞S (18)
其中，S表示所述陨星模型的解体强度；V表示速度幅值；ρair表示空气密度；
S52、模拟所述陨星模型的空爆：
S521、设定爆轰产物压力的作用时间为To，并将所述作用时间To内的相对比容νp设为1，在所述作用时间To内，若所有所述不规则离散元颗粒全部被爆轰作用弹开，即爆炸后所有所述不规则离散元颗粒间失去接触，则爆轰压力作用结束，执行步骤S6；否则，执行步骤S522；
S522、爆轰驱动阶段采用JWL方程计算单位面积爆轰产物压力PJWL：
其中，E0表示单位体积爆轰产物内能；AJ、BJ、R1、R2和ω0分别表示第一、第二、第三、第四和第五无量纲参数；e表示以自然常数e＝2.71828为底的指数函数；
S523、求爆炸的所述不规则离散元颗粒与接触的所述不规则离散元颗粒的爆炸作用面积；
S524、基于所述单位面积爆轰产物压力PJWL和爆炸作用面积AP，计算接触的所述不规则离散元颗粒所受的爆轰接触力；
S6、判断所述陨星模型在高度为0前的质量损失是否与质量相等，即判断所述陨星模型的质量是否变为0，若是，则所述陨星模型完全烧蚀，此时所述陨星模型的运动参数解算停止；否则，即所述陨星模型在高度为0时仍有质量，所述陨星模型以陨石形态降落地球表面，获得所述陨星模型的经纬度，确定陨星模型的落点范围。
进一步，所述步骤S2具体包括以下步骤：
S21、采用多边形接触模型，依次判断两两相邻所述不规则离散元颗粒之间是否发生接触，若是，则执行步骤S22，否则执行步骤S3；
S22、通过层次包围盒和双边链表找到所述不规则离散元颗粒间发生接触时的活动区域；
S23、通过法向Hertz模型和切向无滑动Mindlin-Deresiewicz模型，计算相邻两个所述不规则离散元颗粒之间的接触力Fc：
其中，vi和vj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的质心速度；vr表示相对速度；vn和vt分别表示法向相对速度和切向相对速度；n表示法向量，由第i个所述不规则离散元颗粒的质心指向第j个所述不规则离散元颗粒的质心；ωi和ωj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的角速度；li和lj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的边界球的半径；δn表示法向相对侵入深度；δt表示切向侵入矢量；Fcn和Fct分别表示第i个所述不规则离散元颗粒所受的法向接触力和切向接触力；kn和kt分别表示法向刚度系数和切向刚度系数；Cn和Ct分别表示法向阻尼系数和切向阻尼系数；μ表示摩擦系数；
S24、计算所述不规则离散元颗粒的形变结果，得到所述陨星模型的形变结果。
可优选的，所述步骤S32具体包括以下步骤：
S321、设定所述不规则离散元颗粒的高度h、速度幅值V、经度θ、纬度φ、速度进入角γ和速度方位角ψ；
S322、在气动阻力作用下，构建所述不规则离散元颗粒的轨迹方程为：
其中，gr和gt分别表示径向重力加速度和切向重力加速度；Re表示地球半径；ωE表示地球自转角速度；Cd表示气动阻力系数；m表示陨星质量；A表示阻力参考面积；d表示微分；t表示时间；
S323、考虑所述不规则离散元颗粒之间的接触力时，式(6)中所述速度幅值V、速度进入角γ和速度方位角ψ的解算公式为：
其中，∑Fc表示所述不规则离散元颗粒所受的接触合力；fc1,fc2,fc3分别表示∑Fc在所述东北天坐标系o-xyz中x轴、y轴和z轴上的分量；()new表示迭代新值；
S324、解算所述不规则离散元颗粒的位置和速度，得到所述不规则离散元颗粒的轨迹；
所述步骤S33具体包括以下步骤：
S331、计算所述不规则离散元颗粒的气动阻力矩MD：
S3311、基于克努森数，计算所述陨星模型进入过程中经历自由分子流区时的第一表面压力pfm和第一剪切应力τfm：
其中，ρ∞和V∞分别表示来流密度和来流速度；α表示来流撞击角；
S3312、基于Lees修正的牛顿流模型，计算所述陨星模型进入过程中经历连续流区时的第二表面压力pcont和第二剪切应力τcont：
其中，M∞表示来流马赫数；Г表示气体比热比；CPmax表示正激波后驻点处压力系数；
S3313、基于sine-squared桥函数，计算所述陨星模型进入过程中经历过渡流区时的第三表面压力ptrans和第三剪切应力τtrans：
其中，Kn表示克努森数；lg表示以10为底的log函数；
S3314、根据不规则离散元颗粒中三角面距质心的距离计算力臂，结合自由分子流区、过渡流区和连续流区的第一表面压力pfm、第一剪切应力τfm、第二表面压力pcont、第二剪切应力τcont、第三表面压力ptrans和第三剪切应力τtrans，获得所述不规则离散元颗粒每个三角面相应的气动力矩，进而加和得到总气动力矩MD：
其中，numf表示所述不规则离散元颗粒的面数；si表示第i个三角面的面积；ai表示第i个三角面的力臂；v和ω分别表示所述不规则离散元颗粒的质心速度和相对所述地心固连系O-XYZ的角速度；ni表示第i个面的法向量；pi表示对应的第一表面压力pfm或第二表面压力pcont或第三表面压力ptrans；τi表示对应的第一剪切应力τfm或第二剪切应力τcont或第三剪切应力τtrans；
S332、计算所述陨星模型的角速度ω*和姿态q：
其中，I表示所述不规则离散元颗粒的转动惯量；q0,q1,q2,q3分别表示描述姿态的四元数中第一元素、第二元素、第三元素和第四元素；ωE表示所述地心固连系O-XYZ相对惯性系的角速度且取值为7.292×10-5rad/s；∑Mc表示所述不规则离散元颗粒间的接触力矩之和；-1表示逆运算；T表示转置运算。
可优选的，所述步骤S4具体包括以下步骤：
S41、计算所述陨星模型在进入大气层过程中烧蚀导致的质量损失；
其中，σ表示烧蚀系数且取值范围为3.5×10-10～7.0×10-8kg/J；
S42、计算所述陨星模型的温度变化；
其中，Tm表示所述陨星模型的温度；Tair表示空气温度；Λ表示热流系数；σB表示Stefan-Boltzmann常数；ε表示所述陨星模型的辐射系数；L表示烧蚀热；c表示所述陨星模型的比热容；
S43、基于式(14)，计算所述不规则离散元颗粒单位面积上的线烧蚀率
其中，AS表示所述不规则离散元颗粒的表面积，即所述不规则离散元颗粒中各个三角面的面积之和；Vai,Vbi,Vci分别表示第i个三角面逆时针排列的第一顶点、第二顶点和第三顶点；
S44、基于所述不规则离散元颗粒单位面积上的线烧蚀率
对所述不规则离散元颗粒的顶点进行边界后缩，获得边界后缩后所述不规则离散元颗粒的顶点位置Vj′；
其中，nj表示第j个顶点处相应的外法向量；numa表示包含第j个顶点的三角面的个数；Vak,Vbk,Vck分别表示逆时针排列的第k个包含顶点j的三角面的第一顶点、第二顶点和第三顶点；Vj分别表示边界后缩前所述不规则离散元颗粒的顶点坐标；
S45、依据所述不规则离散元颗粒顶点位置的变化得到所述不规则离散元颗粒的形状变化。
可优选的，所述步骤S23中所述法向刚度系数kn和切向刚度系数kt分别表示为：
其中，Ei和Ej分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的杨氏模量；ιi和ιj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的泊松比；
所述法向阻尼系数Cn和切向阻尼系数Ct分别表示为：
其中，mi和mj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的质量；εn表示法向恢复系数；εt表示切向恢复系数。
可优选的，所述步骤S323中fc1,fc2,fc3由不规则离散元程序在所述地心固连系中计算得出的接触力通过状态转移矩阵变换得到，即：
其中，Fc1、Fc2和Fc3分别表示不规则离散元程序中计算的接触力Fc在所述地心固连系O-XYZ中X轴、Y轴和Z轴上的分量。
可优选的，所述步骤S523中所述爆炸作用面积为接触的所述不规则离散元颗粒在爆炸的所述不规则离散元颗粒的质心指向接触的所述不规则离散元颗粒质心方向上的投影面积；基于Jarvis’March凸包算法将所述爆炸作用面积转化为多边形的面积，所述多边形的面积借助l个三角形面积之和求取，即所述爆炸作用面积AP为：
其中，D1表示所述多边形第一顶点；D(i+1)表示多边形第i+1个顶点且i取1～l；D(i+2)表示所述多边形第i+2个顶点；×表示叉乘运算。
可优选的，所述步骤S31中所述地心固连系O-XYZ与地球固连且原点O位于地球质心，Z轴垂直于地球赤道面指北，X轴为沿赤道平面与本初子午面的交线，Y轴按右手定则确定；所述东北天坐标系o-xyz的原点o为地心与陨星颗粒质心连线与地球表面的交点，z轴沿地球质心指向陨星质心，x轴指向当地铅锤面的东向，y轴根据右手定则指向当地铅锤面的北向；所述步骤S321中所述速度进入角γ为进入速度矢量v在oxy平面内的投影与v的夹角，所述速度方位角ψ为进入速度矢量v在oxy平面的投影与x轴的夹角。
可优选的，所述步骤S323中的∑Fc在爆炸时还包含爆轰产物压力；所述步骤S42中考虑到随着质量损失，陨星表面温度的变化受到陨星在与空气分子相互作用中吸收的热量、黑体辐射损失的热量和物质烧蚀损失的热量的影响，暂不考虑所述不规则离散元颗粒间的热量传递和所述不规则离散元颗粒内部的温度分层。
本发明的另一方面，提供一种利用前述的基于不规则离散元的陨星进入地球大气层后运动仿真方法的运动仿真系统，其包括陨星模型构建模块、颗粒碰撞仿真模块、轨迹和姿态仿真模块、形态仿真模块和落点分析预测模块，所述陨星模型构建模块借助不规则离散元颗粒构建陨星模型，所述颗粒碰撞仿真模块用于计算两两相邻不规则离散元颗粒间的接触力与陨星模型的形变，所述轨迹和姿态仿真模块用于计算不规则离散元颗粒的轨迹以及陨星模型的角速度和姿态变化，所述形态仿真模块用于仿真陨星模型的质量损失、温度变化、形状变化与结构变化，所述落点分析预测模块用于判断陨星模型在高度为0前的质量损失是否与质量相等，若是，则陨星模型的运动参数解算停止；否则，获得陨星模型的经纬度，当高度为0时确定陨星模型的落点范围。
与现有技术相比，本发明的技术效果为：
1、本发明提出的一种基于不规则离散元的陨星进入地球大气层后运动仿真方法，针对单独球体模拟陨星结构的不足，提供一种利用可变形不规则离散元模型，解算碎石堆陨星在超高速进入地球过程中运动演化的数值仿真方法，所提方法能够在陨星进入地球大气的运动解算中考虑不规则陨星内部结构之间的相互作用力，并基于JWL爆轰方程和不规则颗粒间的接触力传播，模拟陨星空爆时的膨胀过程，更符合碎石堆陨星真实的结构特征和爆炸扩张原理。
2、本发明提出的一种基于不规则离散元的陨星进入地球大气层后运动仿真方法，采用不规则离散元模型构建陨星的松散结构，通过拉伸不规则离散元颗粒的多面体网格模拟陨星在烧蚀作用下的边界后缩、体积和外形变化，计算陨星进入地球大气层后的运动和内部组成，模拟陨星解体后碎片的弹出，对于陨星的监测和防御都有着重要的现实与战略意义。
附图说明
通过阅读参照以下附图所作的对非限制性实施例所作的详细描述，本申请的其它特征、目的和优点将会变得更明显。
图1是本发明的基于不规则离散元的陨星进入地球大气层后运动仿真方法流程图；
图2是本发明的陨星进入地球过程中多学科的耦合分析示意图；
图3是本发明的地心固连坐标系与东北天坐标系；
图4是本发明的空爆流程图；
图5是本发明的不规则颗粒投影面积计算示意图。
具体实施方式
下面结合附图和实施例对本申请作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。
图1示出了本发明的基于不规则离散元的陨星进入地球大气层后运动仿真方法，该方法包括以下步骤：
陨星进入地球大气层的不规则离散元仿真中涉及离散元颗粒系统姿轨与结构动力学、超高速跨流域空气动力学、气动热力学、材料热力学、计算爆炸力学等多方面的耦合分析，如图2所示。陨星进入地球的运动参数包括陨星模型的轨迹、角速度变化、姿态变化、质量损失、温度变化、形状变化和结构变化等。
S1、借助不规则离散元颗粒构建陨星模型并初始化：在不规则离散元程序中，采用Delaunay三角网格构建不规则离散元颗粒，若干不规则离散元颗粒组合而成陨星模型，对所有不规则离散元颗粒赋予不同的密度、体积、外形和力学性质参数。
S2、计算两两相邻不规则离散元颗粒之间的接触力，并计算陨星模型的形变。
S21、采用多边形接触模型，依次判断两两相邻不规则离散元颗粒之间是否发生接触，若是，则执行步骤S22，否则执行步骤S3。
S22、通过层次包围盒和双边链表找到不规则离散元颗粒间发生接触时的活动区域。
S23、通过法向Hertz模型和切向无滑动Mindlin-Deresiewicz模型，计算相邻两个不规则离散元颗粒之间的接触力Fc：
其中，vi和vj分别表示第i个不规则离散元颗粒和第j个不规则离散元颗粒的质心速度；vr表示相对速度；vn和vt分别表示法向相对速度和切向相对速度；n表示法向量，由第i个不规则离散元颗粒的质心指向第j个不规则离散元颗粒的质心；ωi和ωj分别表示第i个不规则离散元颗粒和第j个不规则离散元颗粒的角速度；li和lj分别表示第i个不规则离散元颗粒和第j个不规则离散元颗粒的边界球的半径；δn表示法向相对侵入深度；δt表示切向侵入矢量；Fcn和Fct分别表示第i个不规则离散元颗粒所受的法向接触力和切向接触力；kn和kt分别表示法向刚度系数和切向刚度系数；Cn和Ct分别表示法向阻尼系数和切向阻尼系数；μ表示摩擦系数。
法向刚度系数kn和切向刚度系数kt分别表示为：
其中，Ei和Ej分别表示第i个不规则离散元颗粒和第j个不规则离散元颗粒的杨氏模量；ιi和ιj分别表示第i个不规则离散元颗粒和第j个不规则离散元颗粒的泊松比。
法向阻尼系数Cn和切向阻尼系数Ct分别表示为：
其中，mi和mj分别表示第i个不规则离散元颗粒和第j个不规则离散元颗粒的质量；εn表示法向恢复系数；εt表示切向恢复系数。
考虑陨星组成颗粒之间的接触力，从而在气动力基础上更准确地描述陨星的进入解体和爆炸演化过程。
S24、计算不规则离散元颗粒的形变结果，得到陨星模型的形变结果。
S3、计算不规则离散元颗粒的轨迹，并计算陨星模型的角速度和姿态变化。陨星进入过程中，忽略气动升力作用。
S31、建立地心固连系O-XYZ和东北天坐标系o-xyz：如图3所示，地心固连系O-XYZ与地球固连且原点O位于地球质心，Z轴垂直于地球赤道面指北，X轴为沿赤道平面与本初子午面的交线，Y轴按右手定则确定；东北天坐标系o-xyz的原点o为地心与陨星颗粒质心连线与地球表面的交点，z轴沿地球质心指向陨星质心(天向)，x轴指向当地铅锤面的东向，y轴根据右手定则指向当地铅锤面的北向。
S32、计算不规则离散元颗粒的轨迹方程。
S321、设定不规则离散元颗粒的高度h、速度幅值V、经度θ、纬度φ、速度进入角γ和速度方位角ψ。速度进入角γ为进入速度矢量v在oxy平面内的投影与v的夹角，速度方位角ψ为进入速度矢量v在oxy平面的投影与x轴的夹角。
S322、在气动阻力作用下，构建不规则离散元颗粒的轨迹方程为：
其中，gr和gt分别表示径向重力加速度和切向重力加速度；Re表示地球半径；ωE表示地球自转角速度；Cd表示气动阻力系数；m表示陨星质量；A表示阻力参考面积；d表示微分；t表示时间。
S323、考虑不规则离散元颗粒之间的接触力时，式(6)中速度幅值V、速度进入角γ和速度方位角ψ的解算公式为：
其中，∑Fc表示不规则离散元颗粒所受的接触合力，∑Fc在爆炸时还包含爆轰产物压力；( )new表示迭代新值；fc1,fc2,fc3分别表示∑Fc在东北天坐标系o-xyz中x轴、y轴和z轴上的分量，由不规则离散元程序在地心固连系中计算得出的接触力通过状态转移矩阵变换得到，即：
其中，Fc1、Fc2和Fc3分别表示不规则离散元程序中计算的接触力Fc在地心固连系O-XYZ中X轴、Y轴和Z轴上的分量。
S324、解算不规则离散元颗粒的位置和速度，得到不规则离散元颗粒的轨迹。
S33、计算陨星模型的角速度和姿态变化。
S331、计算不规则离散元颗粒的气动阻力矩MD。
S3311、为计算气动阻力矩，需要根据不同的分子流域求解陨星所受气动力。基于克努森数，陨星进入过程会依次经历自由分子流区、过渡流区和连续流区。计算陨星模型进入过程中经历自由分子流区时的第一表面压力pfm和第一剪切应力τfm：
其中，ρ∞和V∞分别表示来流密度和来流速度；α表示来流撞击角。
S3312、基于Lees修正的牛顿流模型，计算陨星模型进入过程中经历连续流区时的第二表面压力pcont和第二剪切应力τcont：
其中，M∞表示来流马赫数；Г表示气体比热比；CPmax表示正激波后驻点处压力系数。
S3313、基于sine-squared桥函数，计算陨星模型进入过程中经历过渡流区时的第三表面压力ptrans和第三剪切应力τtrans：
其中，Kn表示克努森数；lg表示以10为底的log函数。
S3314、已知各个流区的压力和剪切应力后，根据不规则离散元颗粒中三角面距质心的距离计算力臂，结合自由分子流区、过渡流区和连续流区的第一表面压力pfm、第一剪切应力τfm、第二表面压力pcont、第二剪切应力τcont、第三表面压力ptrans和第三剪切应力τtrans，获得不规则离散元颗粒每个三角面相应的气动力矩，进而加和得到总气动力矩MD：
其中，numf表示不规则离散元颗粒的面数；si表示第i个三角面的面积；ai表示第i个三角面的力臂；v和ω分别表示不规则离散元颗粒的质心速度和相对地心固连系O-XYZ的角速度；ni表示第i个面的法向量；pi表示对应的第一表面压力pfm或第二表面压力pcont或第三表面压力ptrans；τi表示对应的第一剪切应力τfm或第二剪切应力τcont或第三剪切应力τtrans。
S332、计算陨星模型的角速度ω*和姿态q：
其中，I表示不规则离散元颗粒的转动惯量；q0,q1,q2,q3分别表示描述姿态的四元数中第一元素、第二元素、第三元素和第四元素；ωE表示地心固连系O-XYZ相对惯性系的角速度且取值为7.292×10-5rad/s；∑Mc表示不规则离散元颗粒间的接触力矩之和；-1表示逆运算；T表示转置运算。
S4、求解陨星模型的质量损失、温度变化和形状变化。
S41、计算陨星模型在进入大气层过程中烧蚀导致的质量损失；
其中，σ表示烧蚀系数且取值范围为3.5×10-10～7.0×10-8kg/J；ρair表示空气密度。
S42、考虑到随着质量损失，陨星表面温度的变化受到陨星在与空气分子相互作用中吸收的热量、黑体辐射损失的热量和物质烧蚀损失的热量的影响，暂不考虑不规则离散元颗粒间的热量传递和不规则离散元颗粒内部的温度分层，计算陨星模型的温度变化；
其中，Tm表示陨星模型的温度；Tair表示空气温度；Λ表示热流系数；σB表示Stefan-Boltzmann常数；ε表示陨星模型的辐射系数；L表示烧蚀热；c表示陨星模型的比热容。
在已有不规则离散元模型中，不涉及不规则离散元颗粒的变形问题。但在陨星仿真中，需要模拟陨星质量烧蚀导致的边界后缩，也即对陨星组成颗粒的不规则离散元进行变形。
S43、基于式(14)，计算不规则离散元颗粒单位面积上的线烧蚀率
其中，AS表示不规则离散元颗粒的表面积，即不规则离散元颗粒中各个三角面的面积之和；Vai,Vbi,Vci分别表示第i个三角面逆时针排列的第一顶点、第二顶点和第三顶点。
S44、基于不规则离散元颗粒单位面积上的线烧蚀率
对不规则离散元颗粒的顶点进行边界后缩，获得边界后缩后不规则离散元颗粒的顶点位置Vj′；
其中，nj表示第j个顶点处相应的外法向量；numa表示包含第j个顶点的三角面的个数；Vak,Vbk,Vck分别表示逆时针排列的第k个包含顶点j的三角面的第一顶点、第二顶点和第三顶点；Vj分别表示边界后缩前不规则离散元颗粒的顶点坐标。
S45、依据不规则离散元颗粒顶点位置的变化得到不规则离散元颗粒的形状变化。
S5、计算陨星模型的结构变化，即模拟陨星模型的解体和空爆过程。
S51、松散结构的陨星在进入过程中很可能伴随一系列的解体和空爆反应，基于陨星强度，获得陨星在气动作用下的解体判据；
ρairV2＞S (18)
其中，S表示陨星模型的解体强度。
S52、除解体事件外，陨星在巨大压力和剧烈加热作用下还会发生空爆，释放出大量的光辐射和热辐射，模拟陨星模型的空爆。
将爆轰产物JWL状态方程与不规则离散元模型结合，模拟陨星的空爆现象。JWL状态方程描述了在爆轰C-J状态之后爆轰产物压力、体积和温度等物理量之间的关系式，反应了爆轰产物的做功能力，已被广泛运用在有限元仿真软件中，例如LS-DYNA等。在不规则离散元模型中应用JWL方程的空爆流程图如图4所示。首先根据温度判定是否有颗粒发生爆炸；之后运用JWL方程对与爆炸颗粒接触的颗粒施加爆轰产物压力。
S521、因无法实际测量陨星爆轰时相对比容的变化，设定爆轰产物压力的作用时间为To，并将作用时间To内的相对比容νp设为1，在作用时间To内，若所有不规则离散元颗粒全部被爆轰作用弹开，即爆炸后所有不规则离散元颗粒间失去接触，则爆轰压力作用结束，作用结束后，首次接受爆轰产物压力的颗粒会获得较大的加速度，颗粒之间通过接触力传递爆炸响应，执行步骤S6；否则，执行步骤S522。
S522、爆轰驱动阶段采用JWL方程计算单位面积爆轰产物压力PJWL：
其中，E0表示单位体积爆轰产物内能；e表示以自然常数e＝2.71828为底的指数函数；AJ、BJ、R1、R2和ω0分别表示第一、第二、第三、第四和第五无量纲参数，一般通过圆筒实验测量。在陨星爆轰中，因无量纲参数的具体数值无法进行实际测量，可根据仿真中所需的爆炸当量设定各参数值。
S523、求爆炸的不规则离散元颗粒与接触的不规则离散元颗粒的爆炸作用面积。爆炸作用面积为接触的不规则离散元颗粒在爆炸的不规则离散元颗粒的质心指向接触的不规则离散元颗粒质心方向上的投影面积；基于Jarvis’March凸包算法将爆炸作用面积转化为多边形的面积，多边形的面积借助l个三角形面积之和求取，即爆炸作用面积AP为：
其中，D1表示多边形第一顶点；D(i+1)表示多边形第i+1个顶点且i取1～l；D(i+2)表示多边形第i+2个顶点；×表示叉乘运算。
在一个具体实施例中，求解投影面积的示意图如图5所示，其中o1和o2分别为爆炸颗粒与其周围接触颗粒的质心。投影平面内圆点表示不规则接触颗粒顶点的投影点，其中蓝色圆点为根据Jarvis’March凸包算法得到的外围点，中间的绿色圆点为被包围点。将蓝色外围点围成的多边形划分成若干三角形，三角形面积之和即为所求的投影面积AP。
S524、基于单位面积爆轰产物压力PJWL和爆炸作用面积AP，计算接触的不规则离散元颗粒所受的爆轰接触力。
综上所述，利用可变形不规则离散元模拟陨星进入时在气动力、气动热、接触力和爆炸等多方耦合作用下，陨星轨迹、姿态、温度、质量、结构和外形的变化。进入过程中，陨星会经历剧烈的烧蚀和减速。
S6、判断陨星模型在高度为0前的质量损失是否与质量相等，即判断陨星模型的质量是否变为0，若是，则陨星模型完全烧蚀，此时陨星模型的运动参数解算停止；否则，即陨星模型在高度为0时仍有质量，则认为该陨星颗粒未被完全烧蚀，陨星模型以陨石形态降落地球表面，获得陨星模型的经纬度，确定陨星模型的落点范围，进而开展陨石回收等工作。
本发明的另一方面，提供一种利用前述的基于不规则离散元的陨星进入地球大气层后运动仿真方法的运动仿真系统，其包括陨星模型构建模块、颗粒碰撞仿真模块、轨迹和姿态仿真模块、形态仿真模块和落点分析预测模块，陨星模型构建模块借助不规则离散元颗粒构建陨星模型，颗粒碰撞仿真模块用于计算两两相邻不规则离散元颗粒间的接触力与陨星模型的形变，轨迹和姿态仿真模块用于计算不规则离散元颗粒的轨迹以及陨星模型的角速度和姿态变化，形态仿真模块用于仿真陨星模型的质量损失、温度变化、形状变化与结构变化，落点分析预测模块用于判断陨星模型在高度为0前的质量损失是否与质量相等，若是，则陨星模型的运动参数解算停止；否则，获得陨星模型的经纬度，当高度为0时确定陨星模型的落点范围。
本发明提出的一种基于不规则离散元的陨星进入地球大气层后运动仿真方法，针对单独球体模拟陨星结构的不足，提供一种利用可变形不规则离散元模型，解算碎石堆陨星在超高速进入地球过程中运动演化的数值仿真方法，所提方法能够在陨星进入地球大气的运动解算中考虑不规则陨星内部结构之间的相互作用力，并基于JWL爆轰方程和不规则颗粒间的接触力传播，模拟陨星空爆时的膨胀过程，更符合碎石堆陨星真实的结构特征和爆炸扩张原理；采用不规则离散元模型构建陨星的松散结构，通过拉伸不规则离散元颗粒的多面体网格模拟陨星在烧蚀作用下的边界后缩、体积和外形变化，计算陨星进入地球大气层后的运动和内部组成，模拟陨星解体后碎片的弹出，对于陨星的监测和防御都有着重要的现实与战略意义。
最后所应说明的是：以上实施例仅以说明而非限制本发明的技术方案，尽管参照上述实施例对本发明进行了详细说明，本领域的普通技术人员应当理解：依然可以对本发明进行修改或者等同替换，而不脱离本发明的精神和范围的任何修改或局部替换，其均应涵盖在本发明的权利要求范围当中。
Claims (9)
1.一种基于不规则离散元的陨星进入地球大气层后运动仿真方法，其特征在于，其包括以下步骤：
S1、借助不规则离散元颗粒构建陨星模型并初始化：在不规则离散元程序中，采用Delaunay三角网格构建不规则离散元颗粒，若干不规则离散元颗粒组合而成陨星模型，对所有所述不规则离散元颗粒赋予不同的密度、体积、外形和力学性质参数；
S2、计算两两相邻所述不规则离散元颗粒之间的接触力，并计算所述陨星模型的形变；
S3、计算所述不规则离散元颗粒的轨迹，并计算所述陨星模型的角速度和姿态变化；
S31、建立地心固连系O-XYZ和东北天坐标系o-xyz；
S32、计算所述不规则离散元颗粒的轨迹方程；
S33、计算所述陨星模型的角速度和姿态变化；
S4、求解所述陨星模型的质量损失、温度变化和形状变化；
S5、计算所述陨星模型的结构变化，即模拟所述陨星模型的解体和空爆过程；
S51、基于陨星强度，获得陨星在气动作用下的解体判据；
ρairV2＞S (18)
其中，S表示所述陨星模型的解体强度；V表示速度幅值；ρair表示空气密度；
S52、模拟所述陨星模型的空爆：
S521、设定爆轰产物压力的作用时间为To，并将所述作用时间To内的相对比容νp设为1，在所述作用时间To内，若所有所述不规则离散元颗粒全部被爆轰作用弹开，即爆炸后所有所述不规则离散元颗粒间失去接触，则爆轰压力作用结束，执行步骤S6；否则，执行步骤S522；
S522、爆轰驱动阶段采用JWL方程计算单位面积爆轰产物压力PJWL：
其中，E0表示单位体积爆轰产物内能；AJ、BJ、R1、R2和ω0分别表示第一、第二、第三、第四和第五无量纲参数；e表示以自然常数e＝2.71828为底的指数函数；
S523、求爆炸的所述不规则离散元颗粒与接触的所述不规则离散元颗粒的爆炸作用面积；
S524、基于所述单位面积爆轰产物压力PJWL和爆炸作用面积AP，计算接触的所述不规则离散元颗粒所受的爆轰接触力；
S6、判断所述陨星模型在高度为0前的质量损失是否与质量相等，即判断所述陨星模型的质量是否变为0，若是，则所述陨星模型完全烧蚀，此时所述陨星模型的运动参数解算停止；否则，即所述陨星模型在高度为0时仍有质量，所述陨星模型以陨石形态降落地球表面，获得所述陨星模型的经纬度，确定陨星模型的落点范围。
2.根据权利要求1所述的基于不规则离散元的陨星进入地球大气层后运动仿真方法，其特征在于，所述步骤S2具体包括以下步骤：
S21、采用多边形接触模型，依次判断两两相邻所述不规则离散元颗粒之间是否发生接触，若是，则执行步骤S22，否则执行步骤S3；
S22、通过层次包围盒和双边链表找到所述不规则离散元颗粒间发生接触时的活动区域；
S23、通过法向Hertz模型和切向无滑动Mindlin-Deresiewicz模型，计算相邻两个所述不规则离散元颗粒之间的接触力Fc：
其中，vi和vj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的质心速度；vr表示相对速度；vn和vt分别表示法向相对速度和切向相对速度；n表示法向量，由第i个所述不规则离散元颗粒的质心指向第j个所述不规则离散元颗粒的质心；ωi和ωj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的角速度；li和lj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的边界球的半径；δn表示法向相对侵入深度；δt表示切向侵入矢量；Fcn和Fct分别表示第i个所述不规则离散元颗粒所受的法向接触力和切向接触力；kn和kt分别表示法向刚度系数和切向刚度系数；Cn和Ct分别表示法向阻尼系数和切向阻尼系数；μ表示摩擦系数；
S24、计算所述不规则离散元颗粒的形变结果，得到所述陨星模型的形变结果。
3.根据权利要求1所述的基于不规则离散元的陨星进入地球大气层后运动仿真方法，其特征在于，所述步骤S32具体包括以下步骤：
S321、设定所述不规则离散元颗粒的高度h、速度幅值V、经度θ、纬度φ、速度进入角γ和速度方位角ψ；
S322、在气动阻力作用下，构建所述不规则离散元颗粒的轨迹方程为：
其中，gr和gt分别表示径向重力加速度和切向重力加速度；Re表示地球半径；ωE表示地球自转角速度；Cd表示气动阻力系数；m表示陨星质量；A表示阻力参考面积；d表示微分；t表示时间；
S323、考虑所述不规则离散元颗粒之间的接触力时，式(6)中所述速度幅值V、速度进入角γ和速度方位角ψ的解算公式为：
其中，∑Fc表示所述不规则离散元颗粒所受的接触合力；fc1,fc2,fc3分别表示∑Fc在所述东北天坐标系o-xyz中x轴、y轴和z轴上的分量；()new表示迭代新值；
S324、解算所述不规则离散元颗粒的位置和速度，得到所述不规则离散元颗粒的轨迹；
所述步骤S33具体包括以下步骤：
S331、计算所述不规则离散元颗粒的气动阻力矩MD：
S3311、基于克努森数，计算所述陨星模型进入过程中经历自由分子流区时的第一表面压力pfm和第一剪切应力τfm：
其中，ρ∞和V∞分别表示来流密度和来流速度；α表示来流撞击角；
S3312、基于Lees修正的牛顿流模型，计算所述陨星模型进入过程中经历连续流区时的第二表面压力pcont和第二剪切应力τcont：
其中，M∞表示来流马赫数；Г表示气体比热比；CPmax表示正激波后驻点处压力系数；
S3313、基于sine-squared桥函数，计算所述陨星模型进入过程中经历过渡流区时的第三表面压力ptrans和第三剪切应力τtrans：
其中，Kn表示克努森数；lg表示以10为底的log函数；
S3314、根据不规则离散元颗粒中三角面距质心的距离计算力臂，结合自由分子流区、过渡流区和连续流区的第一表面压力pfm、第一剪切应力τfm、第二表面压力pcont、第二剪切应力τcont、第三表面压力ptrans和第三剪切应力τtrans，获得所述不规则离散元颗粒每个三角面相应的气动力矩，进而加和得到总气动力矩MD：
其中，numf表示所述不规则离散元颗粒的面数；si表示第i个三角面的面积；ai表示第i个三角面的力臂；v和ω分别表示所述不规则离散元颗粒的质心速度和相对所述地心固连系O-XYZ的角速度；ni表示第i个面的法向量；pi表示对应的第一表面压力pfm或第二表面压力pcont或第三表面压力ptrans；τi表示对应的第一剪切应力τfm或第二剪切应力τcont或第三剪切应力τtrans；
S332、计算所述陨星模型的角速度ω*和姿态q：
其中，I表示所述不规则离散元颗粒的转动惯量；q0,q1,q2,q3分别表示描述姿态的四元数中第一元素、第二元素、第三元素和第四元素；ωE表示所述地心固连系O-XYZ相对惯性系的角速度且取值为7.292×10-5rad/s；∑Mc表示所述不规则离散元颗粒间的接触力矩之和；-1表示逆运算；T表示转置运算。
4.根据权利要求1所述的基于不规则离散元的陨星进入地球大气层后运动仿真方法，其特征在于，所述步骤S4具体包括以下步骤：
S41、计算所述陨星模型在进入大气层过程中烧蚀导致的质量损失；
其中，σ表示烧蚀系数且取值范围为3.5×10-10～7.0×10-8kg/J；
S42、计算所述陨星模型的温度变化；
其中，Tm表示所述陨星模型的温度；Tair表示空气温度；Λ表示热流系数；σB表示Stefan-Boltzmann常数；ε表示所述陨星模型的辐射系数；L表示烧蚀热；c表示所述陨星模型的比热容；
S43、基于式(14)，计算所述不规则离散元颗粒单位面积上的线烧蚀率
其中，AS表示所述不规则离散元颗粒的表面积，即所述不规则离散元颗粒中各个三角面的面积之和；Vai,Vbi,Vci分别表示第i个三角面逆时针排列的第一顶点、第二顶点和第三顶点；
S44、基于所述不规则离散元颗粒单位面积上的线烧蚀率
对所述不规则离散元颗粒的顶点进行边界后缩，获得边界后缩后所述不规则离散元颗粒的顶点位置Vj′；
其中，nj表示第j个顶点处相应的外法向量；numa表示包含第j个顶点的三角面的个数；Vak,Vbk,Vck分别表示逆时针排列的第k个包含顶点j的三角面的第一顶点、第二顶点和第三顶点；Vj分别表示边界后缩前所述不规则离散元颗粒的顶点坐标；
S45、依据所述不规则离散元颗粒顶点位置的变化得到所述不规则离散元颗粒的形状变化。
5.根据权利要求1所述的基于不规则离散元的陨星进入地球大气层后运动仿真方法，其特征在于，所述步骤S23中所述法向刚度系数kn和切向刚度系数kt分别表示为：
其中，Ei和Ej分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的杨氏模量；ιi和ιj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的泊松比；
所述法向阻尼系数Cn和切向阻尼系数Ct分别表示为：
其中，mi和mj分别表示第i个所述不规则离散元颗粒和第j个所述不规则离散元颗粒的质量；εn表示法向恢复系数；εt表示切向恢复系数。
6.根据权利要求1所述的基于不规则离散元的陨星进入地球大气层后运动仿真方法，其特征在于，所述步骤S323中fc1,fc2,fc3由不规则离散元程序在所述地心固连系中计算得出的接触力通过状态转移矩阵变换得到，即：
其中，Fc1、Fc2和Fc3分别表示不规则离散元程序中计算的接触力Fc在所述地心固连系O-XYZ中X轴、Y轴和Z轴上的分量。
7.根据权利要求1所述的基于不规则离散元的陨星进入地球大气层后运动仿真方法，其特征在于，所述步骤S523中所述爆炸作用面积为接触的所述不规则离散元颗粒在爆炸的所述不规则离散元颗粒的质心指向接触的所述不规则离散元颗粒质心方向上的投影面积；基于Jarvis’March凸包算法将所述爆炸作用面积转化为多边形的面积，所述多边形的面积借助l个三角形面积之和求取，即所述爆炸作用面积AP为：
其中，D1表示所述多边形第一顶点；D(i+1)表示多边形第i+1个顶点且i取1～l；D(i+2)表示所述多边形第i+2个顶点；×表示叉乘运算。
8.根据权利要求1所述的基于不规则离散元的陨星进入地球大气层后运动仿真方法，其特征在于，所述步骤S31中所述地心固连系O-XYZ与地球固连且原点O位于地球质心，Z轴垂直于地球赤道面指北，X轴为沿赤道平面与本初子午面的交线，Y轴按右手定则确定；所述东北天坐标系o-xyz的原点o为地心与陨星颗粒质心连线与地球表面的交点，z轴沿地球质心指向陨星质心，x轴指向当地铅锤面的东向，y轴根据右手定则指向当地铅锤面的北向；所述步骤S321中所述速度进入角γ为进入速度矢量v在oxy平面内的投影与v的夹角，所述速度方位角ψ为进入速度矢量v在oxy平面的投影与x轴的夹角。
9.一种根据权利要求1-8所述的基于不规则离散元的陨星进入地球大气层后运动仿真方法的运动仿真系统，其特征在于，其包括陨星模型构建模块、颗粒碰撞仿真模块、轨迹和姿态仿真模块、形态仿真模块和落点分析预测模块，所述陨星模型构建模块借助不规则离散元颗粒构建陨星模型，所述颗粒碰撞仿真模块用于计算两两相邻不规则离散元颗粒间的接触力与陨星模型的形变，所述轨迹和姿态仿真模块用于计算不规则离散元颗粒的轨迹以及陨星模型的角速度和姿态变化，所述形态仿真模块用于仿真陨星模型的质量损失、温度变化、形状变化与结构变化，所述落点分析预测模块用于判断陨星模型在高度为0前的质量损失是否与质量相等，若是，则陨星模型的运动参数解算停止；否则，获得陨星模型的经纬度，当高度为0时确定陨星模型的落点范围。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310159735.8A CN116306185B (zh) | 2023-02-14 | 2023-02-14 | 基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310159735.8A CN116306185B (zh) | 2023-02-14 | 2023-02-14 | 基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116306185A true CN116306185A (zh) | 2023-06-23 |
| CN116306185B CN116306185B (zh) | 2023-11-03 |
Family
ID=86821543
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310159735.8A Active CN116306185B (zh) | 2023-02-14 | 2023-02-14 | 基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116306185B (zh) |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| RU2388061C1 (ru) * | 2009-03-16 | 2010-04-27 | Феликс Герасимович Герасимато | Способ моделирования катастроф, вызванных падением метеоритов |
| CN105739537A (zh) * | 2016-03-29 | 2016-07-06 | 北京理工大学 | 一种小天体表面附着运动主动控制方法 |
| CN108871349A (zh) * | 2018-07-13 | 2018-11-23 | 北京理工大学 | 一种深空探测器光学导航位姿加权确定方法 |
| CN111241634A (zh) * | 2019-11-19 | 2020-06-05 | 中国空气动力研究与发展中心超高速空气动力研究所 | 一种航天器再入陨落的分析预报方法 |
| CN113656896A (zh) * | 2021-09-10 | 2021-11-16 | 中国空气动力研究与发展中心高速空气动力研究所 | 翻转平板在陨落过程中的气动力模型建立方法 |
| CN113722958A (zh) * | 2021-08-30 | 2021-11-30 | 北京理工大学 | 一种不规则形状小天体引力场高效建模方法 |
| CN114565742A (zh) * | 2022-02-16 | 2022-05-31 | 青岛科技大学 | 小天体表面动态模拟与着陆视景仿真系统及方法 |
| CN115292657A (zh) * | 2022-07-04 | 2022-11-04 | 中国科学院国家空间科学中心 | 一种求解撞击地球小天体的物理特性的方法 |
-
2023
- 2023-02-14 CN CN202310159735.8A patent/CN116306185B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| RU2388061C1 (ru) * | 2009-03-16 | 2010-04-27 | Феликс Герасимович Герасимато | Способ моделирования катастроф, вызванных падением метеоритов |
| CN105739537A (zh) * | 2016-03-29 | 2016-07-06 | 北京理工大学 | 一种小天体表面附着运动主动控制方法 |
| CN108871349A (zh) * | 2018-07-13 | 2018-11-23 | 北京理工大学 | 一种深空探测器光学导航位姿加权确定方法 |
| CN111241634A (zh) * | 2019-11-19 | 2020-06-05 | 中国空气动力研究与发展中心超高速空气动力研究所 | 一种航天器再入陨落的分析预报方法 |
| CN113722958A (zh) * | 2021-08-30 | 2021-11-30 | 北京理工大学 | 一种不规则形状小天体引力场高效建模方法 |
| CN113656896A (zh) * | 2021-09-10 | 2021-11-16 | 中国空气动力研究与发展中心高速空气动力研究所 | 翻转平板在陨落过程中的气动力模型建立方法 |
| CN114565742A (zh) * | 2022-02-16 | 2022-05-31 | 青岛科技大学 | 小天体表面动态模拟与着陆视景仿真系统及方法 |
| CN115292657A (zh) * | 2022-07-04 | 2022-11-04 | 中国科学院国家空间科学中心 | 一种求解撞击地球小天体的物理特性的方法 |
Non-Patent Citations (3)
| Title |
|---|
| ZIWEN LI, ET AL: "Numerical Comparison of Contact Force Models in the Discrete Element Method", AEROSPACE, pages 1 - 20 * |
| 崔平远;朱圣英;崔祜涛;: "小天体软着陆自主光学导航与制导方法研究", 宇航学报, no. 06, pages 2159 - 2164 * |
| 耿淑娟;周炳红;韩鹏;郑伟;李明涛;: "小天体撞击地球大气层的空爆问题研究", 空间碎片研究, no. 03, pages 35 - 42 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116306185B (zh) | 2023-11-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111241634B (zh) | 一种航天器再入陨落的分析预报方法 | |
| CN114580224B (zh) | 一种分布式气动融合轨道耦合姿态摄动分析方法 | |
| CN105205248A (zh) | 一种基于ode物理引擎的车辆地形通过性仿真分析组件的设计方法 | |
| CN114936471B (zh) | 一种基于并行计算的航天器碰撞预警分层快速筛选方法 | |
| CN104751012A (zh) | 沿飞行弹道的扰动引力快速逼近方法 | |
| CN116933487A (zh) | 航天器体系仿真碰撞损伤裁决系统及方法 | |
| CN116306185B (zh) | 基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统 | |
| Moss et al. | Survey of blunt body flows including wakes at hypersonic low-density conditions | |
| Schonberg et al. | Empirical hole size and crack length models for dual-wall systems under hypervelocity projectile impact | |
| Tewari | Entry trajectory model with thermomechanical breakup | |
| CN115292657B (zh) | 一种求解撞击地球小天体的物理特性的方法 | |
| Prevereaud et al. | Debris aerodynamic interactions during uncontrolled atmospheric reentry | |
| Schonberg et al. | RCS-based ballistic limit curves for non-spherical projectiles impacting dual-wall spacecraft systems | |
| Shuvalov et al. | Numerical simulation of the LCROSS impact experiment | |
| Szalai et al. | Mars exploration rover heatshield observation campaign | |
| Caldwell et al. | Exploring density and strength variations in asteroid 16 Psyche’s composition with 3D hydrocode modeling of its deepest impact structure | |
| Rakib et al. | A Review of Shielding Systems for Protecting Off-Earth Structures from Micrometeoroid and Orbital Debris Impact | |
| Huxley-Reynard | An airbag landing system for the Beagle 2 Mars probe | |
| Kaplan et al. | Drag analysis of an orbiting tumbling body at the onset of reentry | |
| Orlando | Modeling the Evolution from Massive Stars to Supernovae and Supernova Remnants | |
| Petrescu | Some Aspects of Modern Drones | |
| Brack et al. | Effects of Momentum Transfer Deflection Efforts on Small-Body Rotational State | |
| CN117828766A (zh) | 小行星防御撞击器轨道仿真系统 | |
| Prevereaud et al. | Debris Aerodynamic Interaction and its Effect on Reentry Risk Assessment | |
| Megliola | Simplified structural fragmentation analysis of space debris subjected to destructive re-entry in the atmosphere |
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 |