航天器无控再入解体非规则碎片散布范围分析研究

2020-09-02蒋新宇党雷宁李志辉李四新唐小伟

载人航天 2020年4期

蒋新宇，党雷宁，李志辉*，李四新，唐小伟

（1.中国空气动力研究与发展中心超高速空气动力研究所，绵阳621000）2.国家计算流体力学实验室，北京100191）

1 引言

航天器被送入太空，在航天器服役期满或寿命末期，部分大型航天器（如和平号空间站等）会采取主动离轨再入策略；部分航天器（如同步轨道卫星等）会进入坟场轨道，腾出宝贵的轨道资源，以免对正常轨道上的卫星构成威胁；而绝大部分失效航天器则是沿着螺旋型椭圆轨道衰降无控飞行，最终再入大气层解体坠毁。失效航天器和因空间活动（如太空行走、航天器相互碰撞、燃料泄漏等）产生的小型零部件、残骸等，都统称为空间碎片。空间碎片会威胁在轨航天器的安全运行，需要对这类航天器做好监测防护或主动预报［1-3］，必要时进行主动轨道机动来躲避高风险碎片［4］。随着空间碎片的增多，碎片清理也成为一个越来越受到重视的研究方向。Bérend等［5］给出了一种空间碎片清理策略的双目标优化方法；Forshaw等［6］完成了一次在轨碎片清理试验，分别对编网捕获、鱼叉捕获、视觉导航交会、阻力帆技术进行了验证，展示了未来可能的轨道碎片清理技术。受大气阻力、空间环境和引力等因素影响，空间碎片的轨道高度都会逐渐衰降，特别是低轨道碎片受大气阻力影响明显，800 km以下的空间碎片一般会在十余年内再入大气层［7］。目前每年都会发生100次左右的航天器再入事件，除去有特殊热防护设计的返回式卫星/飞船等，大部空间物体会在再入环境的强气动力/热作用下熔融/烧蚀/解体，大部分质量熔融/烧蚀殆尽，但一般仍有10%～40%的残余质量会到达地面，并可能对地面造成危害［8］。

个别较特殊的再入事件会人为主动处理，如美国军方声称为避免有毒材料危害地球，采用标准-3导弹在约250 km高度将USA-193卫星击毁。公开资料显示大部分碎片会在24～48 h内坠落，余下的碎片则在40天内全部坠落。Pardini等［9］采用公开的观测数据对此事件进行了模拟，研究了卫星被摧毁后的碎片云散布范围和演化过程，并希望有空间监测能力的国家、机构对此类信息做到公开透明，以确保公共空间安全。

对于空间物体再入的研究主要分为2个方向，一个方向是再入环境下的材料毁坏，如Prevereaud等［10］对航天器再入过程的烧蚀破坏作了数值分析，并采用试验方法研究了合金材料的氧化/熔融过程；另一个方向是复杂外形物体的弹道计算与落点散布分析，如Mehta等［11］建立了一种考虑参数不确定性的高维模拟方法，其模拟的落点较常规高斯分布和椭圆分布相去甚远；Falsone等［12］提出了一种基于仿真的碎片落点散布分析方法，结果优于之前使用的协方差传播方法，并认为以解体点为模拟起始点是一个有趣的研究方向。

目前主要的空间物体再入预测软件包括美国的DAS［13］、ORSAT［14］和欧洲的SCARAB［15］等。胡锐锋等［16-17］也开发了空间碎片方面的工程预测软件DRAPS。空间物体再入预测软件按分析方法可以分为面向物体法和面向航天器法，除欧洲的SCARAB外，上述软件均采用了面向物体法［18］。面向物体法将航天器以及解体后产生的部件外形认为由简单的基本几何模型组成，且其几何信息可用若干特征参数描述；面向航天器法则尽可能模拟真实的航天器外形，采用基于表面网格的飞行器和部件模型，理论上具有更高的预测精度，但建模较复杂并需要较长的计算时间。以上软件主要是模拟了航天器部件的再入情况，分析可能到达地面的残骸，分析对象都是已知且提前预设好的航天器结构部件，暂时还没有考虑航天器在再入过程中因熔融、烧蚀、爆炸等因素突然产生的大量微小碎片。微小碎片因为数量大，散布范围广，可能造成更大的地面伤害概率，因此也需要研究。

本文旨在开展航天器再入解体后产生的大量非规则微小碎片地面散布范围可计算建模，并研究其分布规律。根据碎片几何形状与尺度，将碎片分为有限类，利用空间分组策略实现对碎片运动随机性的模拟，使用当地化快速方法计算不规则碎片的空气动力特性，采用六阶显式Adams积分方法求解弹道方程，最终获得碎片的地面散布范围。

2 非规则微小碎片散布范围可计算建模

2.1 碎片的几何分类与质量分布模型

采用面向物体法的思想，将碎片抽象为具有简单几何形状的有限类，对各类碎片进行弹道模拟，落点相近的类可以整合。根据大量模拟结果，本文将碎片按几何形状分为块状、片状和杆状等5类。同时注意到碎片面质比对弹道特性影响很大，所以在确定了几何分类方案后，还要考虑碎片尺度的影响。本文将碎片按尺度分为3类，最终确定的碎片分类方案见表1所示。某类碎片的单个质量越小，则该类碎片的数量就越大，且质量和数量在取对数后基本呈线性关系［19］。为此，可设置n为某类碎片的数量，如式（1）所示：

式中，n为该类碎片的单个质量，C基于质量约束条件确定。k根据解体类型不同而不同，在本文中取为0.553。

表1 碎片有限分类方案*Table 1 Finite classifying scheme for debris

结合表1与式（1），可确定各类碎片的权重。例如假设碎片总质量为200 kg，3个尺度取为0.1m，0.01 m和0.001 m，则各类碎片数量与质量的关系如图1所示。

2.2 碎片运动随机性的模拟

服役期满航天器再入多次解体生成大量碎片之前，航天器作为一个整体运动，具有确定的空间位置、速度、姿态角。在熔融、烧蚀、爆炸等因素下突然生成碎片的时间点，产生大量的不同形状、不同尺度的碎片，在按形状、尺度将碎片分为有限类以后，各类碎片的数目依然众多。采用李志辉等［20］研究建立的箔条云随机统计模拟方法来处理每一类碎片。同一类中的各碎片具有随机的飞行姿态角，以相同的分离速度（绝对值）沿不同的空间方位角散开［21］，不同种类碎片的分布相互独立。为实现对同类碎片的跟踪模拟，可采用如下规则：

图1 各类碎片的数量与质量关系Fig.1 Relationship between quantity and mass of various debris

1）将整个空间方位角均匀划分为ng=Mφ×Mθ个组，每组具有确定的方位角 ( φ，θ)；

2）每组根据该组的方位角 ( φ，θ)和相对分离速度得到本组碎片生成后的初速度，初速度保证了各碎片组之间具有相互散开的趋势；

3）对于具有N1个碎片的种类1，通过生成N1个随机数R1，再按规则映射到某个方位角中去即可。同理进行其它种类碎片的映射；

4）若某个种类的碎片数目Ni≫ng，由于在空间方位上是均匀分布，为减少随机数的生成和映射，可以先将大部分数目（比如80%）的碎片平均分配到各空间方位角，再将余下部分按规则（3）进行映射。

为实时跟踪有限的Mφ×Mθ组碎片运动飞行情况，需要将Mφ×Mθ组碎片一一映射分布到具有（αL，βM）飞行迎角与侧滑角的各小组Mα×Mβ中，使得每组碎片具有若干碎片个体，且在每个d t时间内飞行姿态相对不变［21］。

依次处理每个组［21-22］：先生成一个随机数，再把随机数按规则映射到一个角度值，如果该角度值未使用过，就进行随机筛选；如果已经使用过，则需要重新生成随机数。针对大型航天器再入多次解体形成碎片问题，需要重复多次生成合适随机数，采用如下规则进行碎片组飞行姿态角的随机筛选确定：

1）一次性生成Mφ个均匀随机数，存放在数组Rφ中；

2）对Rφ中的元素进行排序，记录每个元素所对应的序号，存放在数组Iφ中；

3）将飞行迎角均匀划分成Mφ个组，Iφ即可根据序号映射到不同的飞行迎角；

4）同理，可对飞行测滑角进行类似的映射；

5）分别从2个方向完成对Mα×Mβ碎片组的上述循环，即可确定各碎片组的飞行迎角和侧滑角。

使用上述随机动力学筛选方法，可将Mφ×Mθ组碎片一一映射到具有（αi，βj）飞行迎角与侧滑角的各小组Mα×Mβ中，使得同一类碎片能被模拟成有限的碎片组，各碎片组在每个d t时间内具有确定的碎片个体数Nij，以确定的飞行速度Vij、飞行迎角αi和侧滑角βj飘落飞行，于是可研究发展相关空气动力学理论、方法，计算确定此d t时间内碎片个体所受的空气动力学特性［23-25］。

2.3 再入跨流域气动特性的当地化快速算法

常用的数值方法如直接模拟蒙特卡罗（DSMC）方法、N-S方程解算器与求解玻尔兹曼（Boltzmann）模型方程、气体动理论统一算法（GKUA）等，由于计算量太大，不能满足弹道计算中的气动特性快速获取要求，本文采用了当地化快速算法［26］。

当地化快速算法以稀薄气体高超声速绕流当地化面元法计算理论为基础，将非规则解体物面划分成若干块小曲面，对于每个小曲面，选用一个小的平面三角形面元或四边形面元来代替，利用面元逼近复杂解体物形。这样，计算小曲面上的气动力就转换成计算面元上的气动力，将这些面元上的气动力累加起来就可以得到整个非规则解体物的气动特性。当面元分得极其细小时，由此引起的误差就会变得很小。

当地面元所受气动力系数依赖于来流和当地性质，如当地迎角、当地表面作用等。对于给定的入射角，面元上经归一化的压力和摩擦力系数表达为式（2）：

式中，CP和Cf分别是过渡流区域的压力和摩擦力系数，下标FM和Cont分别表示自由分子流和连续流；FP，b和Ff，b分别是压力和摩擦力桥函数，下标b表示桥函数，依赖于独立参数当地努森数Kn、壁温总温比TW/T0和入射角θ。桥函数的意义在于使气动系数在自由分子流区和连续流区分别趋近于理论值，而在中间过渡流区实现光滑过渡。

分别使用气体动理论统一算法（GKUA）［23-24］、DSMC方法［25，27］或N-S方程解算器对残骸碎片简化外形典型绕流状态进行精细计算，获取典型状态的精确气动数据。调试修正非规则碎片的桥函数关联参数，使得当地化快速算法的结果在典型状态与数值算法吻合后，则其它任意状态的气动数据可由式（2）直接计算得出。

2.4 碎片运动弹道计算

2.2 节的随机分组方法已经考虑了飞行姿态的影响，则碎片运动的弹道计算中只需计算三自由度弹道。矢量形式的弹道方程如式（3）所示：

式中，n为碎片的瞬时质量；V为碎片的瞬时速度矢量；R为碎片空气动力；G为地球引力，G=n▽U；U为地球的引力势；Fe为离心惯性力，Fe=n rΩ×( Ω×r)；Ω为地球的自转角速度矢量，r是地心惯性坐标系中质心的矢径；Fco为科氏惯性力，Fco=2 n( Ω×V)。

为求解弹道方程（3），采用六阶显式Adams-Bashforth格式，在J2000.0地心惯性坐标系下进行数值积分。对于微分方程式（4）

其六阶显式Adams-Bashforth格式为式（5）：

式中，h为时间步长，un+1为待求的下一时间步函数值，un为当前时刻函数值，fn，fn-1，…，fn-5为当前时刻及当前时刻之前1～5时间步的右函数值。在计算启动时刻，由于时间步结果不足，无法直接采用此数值格式，需要采用龙格-库塔等其它数值方法来计算初始的几个时间步。

3 碎片散布范围计算分析

根据第2节航天器再入多次解体后产生的大量碎片存活坠落地面散布范围可计算模型方法，对航天器解体后的碎片进行气动融合弹道的全程模拟。本文以数值预报某大型航天器再入为例，其解体点位置为北纬42°，东经0°，海拔高度100 km；解体点速度为7500 m/s，朝向正东；弹道倾角为0.1°；碎片总质量为200 kg。对该航天器无控再入解体碎片散布范围计算分析，图2、图3分别绘出碎片在地面散布的数量与质量分布。可以看出，在该给定初始状态下，解体碎片的纵向散布范围可达5 000多公里，而且碎片质量越大，其数量也就越少，航程也就越远。这是因为在再入环境主要受重力和空气动力影响，空气动力与尺度平方（横截面积）正相关，而重力与尺度立方（体积）正相关。在尺度减小时，重力下降得更快，相对而言，空气动力对碎片运动的影响就逐渐增大。图4绘出特征尺度为0.1 m的片状碎片在地面的分布范围。计算表明，由于考虑到了飞行姿态的随机性，碎片在横向也会散开，但横向航程较纵向航程小很多，宽度约为100 km。

图2 碎片在地面散布的数量分布（φ=0.1°）Fig.2 Quantity distribution of debris on the ground（φ=0.1°）

图5绘出特征尺度为0.01 m的立方体状碎片飞行轨迹与碎片云质心弹道的关系，其中蓝色点划线为碎片云质心运动弹道，彩色圆圈为该组碎片不同时刻的空间位置，其颜色代表飞行马赫数。由图中可看出碎片飞行轨迹与质心运动弹道的一致性很好。

图3 解体碎片在地面的质量分布（φ=0.1°）Fig.3 M ass distribution of debris on the ground（φ=0.1°）

图4 特征尺度为0.1m的片状碎片在地面的分布Fig.4 Distribution of sheet-like debris（characteristic scale is 0.1 m）on the ground

图5 特征尺度0.01 m立方体状碎片飞行轨迹与质心弹道Fig.5 Flight path of cube-like debris w ith the characteristic scale of 0.01 m

为了研究航天器再入解体弹道倾角对落区散布影响，采用前述拟定计算状态，仅将弹道倾角变更为1.1°，所得落区散布见图6、图7所示。对比分析图2、图3可以发现，碎片的纵向散布范围大大减小，约为2500 km。该结果表明弹道倾角对落区散布范围影响很大，弹道倾角越大，再入飞行时间越短，落区散布范围越小，在后续服役期满大型航天器受控再入解体数值预报与离轨控制策略的制定与分析中，应特别注意这方面的研究。