卫星解体碎片生成数值模拟

2014-03-30张晓天贾光辉

实验流体力学 2014年4期

张晓天, 贾光辉

(北京航空航天大学宇航学院, 北京 100191)

0 引言

航天器解体碎片生成研究是空间碎片环境建模的关键。目前的主要研究手段是对解体碎片数据进行统计分析，建立经验公式。解体碎片数据来源于两个方面：轨道碎片雷达观测和地面解体实验。

美国太空观测网(SSN，Space Surveillance Network)对地球轨道的物体进行跟踪和编目。以目前的观测能力，能够对低地球轨道特征长度10cm以上的物体和高地球轨道特征长度1m以上的物体进行可靠跟踪，难以获得小尺寸碎片的数据。

在空间碎片群体中，小碎片的数目要远大于大碎片。为了研究小碎片的特性，NASA 开展了SOCIT(Satellite Orbital Debris Characterization Impact Test)地面解体实验[1]。解体实验中使用了OSCAR 小卫星，通过加速铝弹丸对其碰撞造成解体。日本T. Hanada 研究小组与NASA 约翰逊研究中心合作，开展了复合材料小卫星地面解体实验研究[2]。柳森、黄洁、兰胜威和李毅等开展了小型卫星模型的解体实验研究[3-5]，进行了三种卫星模型的撞击解体实验，并提出了新的解体阈值。王海福、余庆波和金学科等开展了圆柱壳结构爆炸解体实验和数值仿真研究[6-7]。地面解体实验受到加速设备能力限制，难以构造常规尺寸卫星解体的碰撞条件，地面解体实验难度很大。

本文通过超高速碰撞数值模拟技术对航天器解体过程进行研究，从而获得解体碎片的特性。数值模拟技术是一种相对成本低的技术，一旦建立了有效的模拟算法，就可以大量开展。由于目前广泛使用的SPH方法在单个碎片的识别和置信度方面存在不足，本文通过一种SPH和有限元结合的方法来实现对碎片的识别和统计，并与卫星模型地面解体实验结果进行对比。

1 卫星解体仿真方法

目前的超高速撞击数值模拟研究主要针对航天器防护结构分析，广泛使用的数值模拟方法是SPH(Smoothed Particle Hydrodynamics)方法[8]。但在航天器防护结构仿真中，由于撞击体较小，航天器一般不会发生解体；另外这种分析的关注点为防护结构是否被碎片穿透，并不关注撞击产生的二次碎片云中单个碎片的信息(如碎片数目、大小和速度矢量等)。而解体模型就是要根据解体条件给出解体碎片特征信息，解体碎片特性是解体动力学数值模拟最关心的结果。因此，将现有数值模拟算法成果应用于航天器解体动力学数值模拟，首先需要对单个碎片的识别性和置信度等方面进行分析。

对于超高速碰撞问题，SPH方法的仿真结果常表现为离散的粒子云，难以从中识别出哪些粒子构成一个碎片，也难以分辨破碎材料和未破碎材料的界面。虽然可以通过搜索算法直接将结果中邻近的粒子划分为一个碎片，但是这种方法存在几个方面的不足。第一，置信度问题。在初始模型中，粒子被以一定间距δ填充入物体所占据的空间。此时粒子作为1个群体代表了1个物体，其可信性较高。但是单个粒子代表了1块δ立方的物质微元，其作为孤立个体可信性较低。因为该物质微元有可能进一步分裂破碎，而一个粒子不可能再被分裂。所以仿真结果中的孤立粒子，甚至是少量几个粒子都不能作为置信的单个碎片存在。第二，流动性问题。超高速碰撞过程中，物质流动范围很大，初始不相邻的物质微元完全有可能在撞击过程中被“挤”到一起。因此仅仅根据某个时间点上粒子的空间坐标是否相邻来判断这两个物质微元是否属于同一个碎片是不够可靠的。第三，邻近阈值问题。判断粒子是否邻近，一般是判断两者间距是否小于给定阈值。这个阈值的给定缺乏力学依据，难以制定统一的规则，这会影响不同问题之间的一致性。

SPH方法在航天器解体碎片仿真问题中存在上述不足；另一方面有限元方法材料界面清晰，碎片效果直观，单个碎片置信度高，但是在大变形问题中会发生网格畸变，使得计算无法进行。将SPH与有限元方法相结合的方法，称为有限元重构方法，提高了对单个碎片的识别性和置信度。在有限元重构方法中首先通过立方体填充，将任意形状的几何体离散为粒子-形心重合模型；进而通过显式积分迭代步中的SPH积分、有限元重构与失效分析三个步骤，获得包含FE单元与SPH粒子的混合仿真结果。

1.1粒子-形心重合模型

图1是一个离散化的实例。A代表一个任意形状的物体，B表示使用立方体FE单元对物体区域进行填充，C表示在每个单元的形心处创建一个SPH粒子，相应的粒子间距等于FE的单元边长。此时的模型称为“粒子-形心重合模型”，包括FE网格与SPH粒子两部分。

A B C D E

1.2显式积分迭代步

对于粒子-形心重合模型中的粒子部分，按照SPH控制方程进行显示时间积分迭代，在每个时间步后能够获得所有粒子的位置、速度信息。

AEi={EH|Ni∈EH}

(1)

在立方体填充模型中，根据节点被多少个单元所共有——也就是AEi集合中单元的个数——共分为8种情况。如果1个节点被8个单元共有，则称其为内点，否则称其为边界点。有限元重构分为3个步骤依次执行：内点重构、第一类边界点重构，第二类边界点重构。

(1)内点重构

用IN表示立方体填充网格中的内点集合，它包含了所有被8个单元共享的节点。内点的重构坐标定义为这8个单元对应的SPH粒子当前坐标的算术平均：

(2)

其中，pH表示EH形心对应的SPH粒子的坐标向量，ri表示Ni的坐标向量。

(2)第一类边界点重构

定义一个与Ni关联的单元标识集合：

(3)

第一类边界点集合定义为：

(4)

第一类边界点的重构坐标按下式进行计算：

(5)

(3)第二类边界点重构

定义一个与Ni关联的单元标识集合：

(6)

第二类边界点集合定义为：

(7)

第二类边界点的坐标按下式进行计算：

(8)

如果一个边界点所有的背点都是第二类边界点，将陷入一个死循环，这些点无法被更新。为避免这种极端情况，在几何体离散的过程中需要对最薄的部位采取两层以上的单元来填充。

在大变形问题中，SPH粒子流动范围很大，重构后的单元不一定能维持常规的形状，很多单元会发生畸变，此时失效单元将被删除，单元形心对应的SPH粒子被保留；未失效单元被保留，对应的SPH粒子被删除。本文采用Johnson-Cook失效准则。FE重构分析和失效分析按照SPH计算的时间步逐步进行，第n+1个时间步的分析继承第n个时间步的结果。在结果中FE单元代表未失效材料，SPH粒子代表失效材料。

2 两体碰撞解体算例

使用上述方法对实心球体撞击实心立方体解体问题进行仿真，两物体相对位置关系如图2所示，从右上角局部放大图中可以看出粒子和单元形心的重合关系。球体的直径为1m，立方体边长1m。两物体的材料密度均为0.2g/cm3，球体质量为104.72kg，立方体质量为200kg。图中的撞击速度v=7.6km/s，相对撞击速度为15.2km/s。建模中有限元单元边长定为2cm，因此在仿真分析中可以得到边长大于2cm的置信碎片结果。

图2 球体与立方体碰撞解体仿真模型

图3给出了350μs时刻碎片云的对比。左图给出了SPH仿真结果，从图中可以看出，SPH的结果表现为粒子云。由局部放大图中难以分辨出哪些粒子构成一个碎片，破碎材料和未破碎材料的界面难以清晰定义。中图为有限元重构的结果，结果由SPH粒子和有限元单元混合表达。在局部放大图中，由单元表达的单个碎片清晰可见。右图为有限元重构结果中的单元部分，也就是置信碎片部分。SPH粒子代表了充分破碎的材料，作为整体表达了碎片云的形状、演化与材料分布，但是作为单个碎片是不够可信的。

图3 球体与立方体碰撞解体数值模拟结果

固定单元是否失效的状态不变，将节点和粒子坐标反推回0时刻，就得到了反推图。图4所示为350μs的反推图。其中由单元表示的部分在撞击过程结束后会转变为置信碎片，或者说是较大的碎片；由粒子表示的部分在撞击结束后会转变为非置信碎片，这些部分会充分破碎为小碎片云。由图中可见，撞击区域内的材料大部分都充分破碎了，大碎片主要由远离撞击区的材料形成。

图4 球体与立方体碰撞结果反推图

有限元重构方法识别了SPH仿真结果中的碎片，确定了清晰的破碎与未破碎材料的界面，对SPH仿真结果中单个碎片的置信部分和非置信部分进行了清晰的划分。首先，在解体过程中如果单元失效，它将被转换为粒子；而粒子继续受变形破裂，就无法描述，因为其不可再分。因此在解体数值模拟结果中，以单元形式存在的物质就构成了置信碎片；而以粒子形式存在的失效单元则属于非置信碎片。对于单个碎片的分析应该基于置信碎片或者说基于单元部分进行。第二，由于引入了单元，刻画了单元之间的物质连通性，初始相邻的单元如果解体结束后仍然相邻则属于同一个碎片；如果初始不相邻，而在解体结束后相邻，则说明它们的邻近是由于材料的运动引起的，因此不构成同一个碎片。这样，哪些单元构成一个碎片将由单元部分的拓扑结构决定。第三，单元转换粒子判据是材料失效或相变，而不是人为指定几何阈值，因此置信与否的分界线从根本上是由力学机理决定的。

3 碎片统计方法

在有限元重构的仿真结果中，置信碎片由单元部分表达，碎片的总数以及碎片的单元构成是由单元部分的拓扑结构决定的。碎片统计就是要设计一种自动的算法，能够统计出碎片的数目、碎片的单元构成以及每个碎片的尺寸、质量和速度矢量等信息。

把置信部分的单元看作顶点，如果两个单元之间存在一个公共节点，则在两单元之间定义一条边，这样顶点和边的集合可以表示为：

(9)

〈Vtx,Edg〉构成一个无向图。如果该图中存在一个首尾相接的边序列，能够把图中两个顶点连接起来，则这两个顶点是连通的。如果〈Vtx,Edg〉的一个子集中任意两个顶点是连通的，则称这个子集为该无向图的一个连通分支。每一个连通分支就对应了碎片云中的一个碎片。因此，碎片的统计与特征量计算就是对连通分支的分析。

本文使用了基于广度优先的连通分支遍历算法，时间复杂度为o(n2)。对图4中的碎片进行分析，图5给出了连通分支搜寻的结果，其中每一种颜色代表了一个连通分支，也就是一个碎片。该算例中共产生了大于2cm的碎片562个，其中178来自球体，384个来自立方体。

在获得了碎片的单元构成之后，可以对每个碎片的特征长度进行计算：

(10)

图5 碎片统计结果

4 卫星模型解体碎片数值模拟

利用上述方法对文献[3]开展的地面卫星模型撞击解体实验进行数值模拟。

图6给出了卫星解体仿真模型，图7为Case B的仿真结果。从图7中可以看出，虽然弹丸尺寸较小，但是造成了很大的破坏。由于卫星模型在撞击法平面内的两个轴向并不对称，所以结构中的开孔也不对称，而是呈椭圆形。仿真结果中，卫星结构并没有像实验中一样板结构之间互相脱离。这是因为在实验结构中，相邻板结构采用了L型钣金配合螺钉连接，连接强度较弱；在仿真中，连接形式进行了近似处理，直接采用了相邻边焊接的方式。

图6 卫星模型撞击解体仿真建模

图8给出了Case A和Case B的累计碎片数目对碎片尺寸的分布曲线对比。总体来说仿真结果与实验数据比较一致。仿真曲线横坐标值在0.001以下的缺失是由模型分辨率引起的，为了获得更小碎片的特性，需要进一步提高模型分辨率。另外，在小尺寸碎片区域，NASA模型的计算值明显高于实验值，这可能是由于NASA模型的建模数据源与文献[3]实验的差异造成的。NASA模型小碎片特性部分的建模数据主要来源于SOCITE4解体实验，该实验中用到的OSCAR卫星的内部结构与材料和文献[3]实验中用的卫星模型相差较大。这也体现了卫星的结构和材料对解体碎片特性的影响，表明了基于解体动力学过程分析、计算碎片特性的必要性。

图7 Case B数值模拟结果

(a) Case A

(b) Case B

5 结论

本文基于超高速碰撞数值模拟，建立了一种用于模拟卫星解体碎片生成的有限元重构方法，对卫星模型解体问题进行了仿真研究，得到如下结论：

(1) 有限元重构方法能够对卫星解体生成的碎片进行模拟，结合无向图连通域分析方法，能够获得碎片分布以及单个碎片的数据；

(2) 对地面卫星模型解体实验进行了数值模拟分析。对于Case A和Case B两种工况，模拟获得的累计碎片数对碎片尺寸的分布曲线与实验曲线基本一致。表明了本文方法的有效性。另外，通过和NASA模型计算结果对比，表明解体卫星的结构形式和材料对解体碎片特性有较大影响，通过对解体动力学过程的研究，可以更加有效的计算解体碎片特性。

参考文献:

[1]Krisko P H, Horstman M, Fudge M L. SOCIT4 collisional-breakup test data analysis: With shape and materials characterization[J]. Advances in Space Research, 2008, 41(7): 1138-1146.

[2]Hanada T, Liou J C, Nakajima T, et al. Outcome of recent satellite impact experiments[J]. Advances in Space Research, 2009, 44(5): 558-567.

[3]Lan S, Li Y, Huang J, et al. Fragment characteristic of simulated spacecraft under hypervelociy impact[C]//Proceeding of 62nd International Astronautical Congress. Cape Town, SA, 2011.

[4]兰胜威, 李毅, 黄洁, 等. 模拟航天器超高速撞击解体碎片特性研究[C]//第六届全国空间碎片学术交流会论文集, 四川成都, 2011.

Lan Shengwei, Li Yi, Huang Jie, et al. Investigation on the characteristics of fragments generated by the simulated spacecraft hypervelocity impact test[C]//Proceeding of 6thNational conference of space debris. Chengdu, Sichuan, 2011.

[5]李毅, 黄洁, 马兆侠, 等. 一种新的卫星超高速撞击解体阈值模型研究[C]//第六届全国空间碎片学术交流会论文集. 四川成都, 2011.

Li Yi, Huang Jie, Ma Zhaoxia, et al. Investigation on new

breakup threshold model of satellite breakup under hypervelocity impact[C]//Proceeding of 6thNational conference of space debris. Chengdu, Sichuan, 2011.

[6]余庆波, 金学科, 王勤智, 等. 结构强度对卫星爆炸解体影响研究[C]//第六届全国空间碎片学术交流会论文集. 四川成都, 2011.

Yu Qingbo, Jin Xueke, Wang Qinzhi, et al. The influence of structure strength on satellite breakup caused by explosion[C]//Proceeding of 6thNational conference of space debris. Chengdu, Sichuan, 2011.

[7]金学科, 余庆波, 王勤智, 等. 燃料量对卫星爆炸解体影响研究[C]//第六届全国空间碎片学术交流会论文集. 四川成都, 2011.

Jin Xueke, Yu Qingbo, Wang Qinzhi, et al. The influence of the amount of fuel on satellite breakup caused by explosion[C]//Proceeding of 6thNational conference of space debris. Chengdu, Sichuan, 2011.

[8]Johnson G R, Stryk R A, Beissel S R. SPH for high velocity impact computations[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 57(4): 573-580.

作者简介: