一种碎片对航天器结构破坏效应的粒子评估方法

摘要

本发明涉及一种碎片对航天器破坏作用的计算评估方法，包括：建立工程计算模型，确定模型的几何尺寸参数；收集步骤1的计算模型中碎片及航天器材料的物理参数；以航天器的大小作为计算域，并采用正交网格离散计算域，网格步长由碎片及航天器中的最小部件尺寸确定，并在网格内加入PIC粒子；获取步骤3所确定的各网格和PIC粒子影响碰撞毁伤效果的物理量；连续输出步碰撞过程中变化的物理量，采用可视化技术，直观地得到碎片对航天器破坏效果图。该方法保证侵彻计算过程中弹靶接触区域的计算精度，减少计算量，不仅适用于碎片对航天器结构破坏的评估，还适用于冲击碰撞问题的分析评估。

说明书

技术领域

本发明涉及一种碎片对航天器结构破坏效应的粒子评估 方法，属于航天技术领域。

背景技术

随着人类在太空活动的日益增多，散布在太空中的大量 空间碎片对航天器的结构安全构成很大的威胁。因此了解碎 片对航天器的破坏作用并进行评估，将有助于航天器的优化 设计。

碎片对航天器的超高速碰撞问题是复杂的非线性动态响 应过程，既有结构发生大位移时所产生的几何非线性，又有 材料发生大变形时所表现出的物理非线性，还存在着复杂的 运动非线性以及复杂的接触和摩擦问题等。

现有的各种理论和实验研究方法都有各自不同的应用范 围与条件。随着计算机技术的发展及数值方法的发展，可以 采用数值手段将碎片对航天器的破坏的进行计算评估，从而 较全面地反应碰撞过程中航天器及碎片参数和物理量的变 化，直观演示航天器的破坏过程。

二维轴对称柱坐标系下，动能弹侵彻问题可用如下非守 恒形式的方程组形式表示：

$\frac{\partial \rho }{\partial t}+u_z\frac{\partial \rho }{\partial z}+u_r\frac{\partial \rho }{\partial r}+\rho (\frac{{\partial u}_z}{\partial z}+\frac{{\partial \mathrm{ru}}_r}{r\partial r})=0---\left(1\right)$

$\rho (\frac{{\partial u}_z}{\partial t}+u_z\frac{{\partial u}_z}{\partial z}+u_r\frac{{\partial u}_z}{\partial r})=-\frac{\partial p}{\partial z}+\frac{{\partial S}_{\mathrm{zz}}}{\partial z}+\frac{\partial \left({\mathrm{rS}}_{\mathrm{rz}}\right)}{r\partial r}---\left(2\right)$

$\rho (\frac{{\partial u}_r}{\partial t}+u_z\frac{{\partial u}_r}{\partial z}+u_r\frac{{\partial u}_r}{\partial r})=-\frac{\partial P}{\partial r}+\frac{{\partial S}_{\mathrm{rz}}}{\partial z}+\frac{\partial \left({\mathrm{rS}}_{\mathrm{rr}}\right)}{r\partial r}+\frac{S_{\mathrm{rr}}+S_{\mathrm{zz}}}{r}---\left(3\right)$

$\rho (\frac{\partial e}{\partial t}+u_z\frac{\partial e}{\partial z}+u_r\frac{\partial e}{\partial r})=-p(\frac{\partial \left({\mathrm{ru}}_r\right)}{r\partial r}+\frac{{\partial u}_z}{\partial z})$

$+S_{\mathrm{zz}}\frac{{\partial u}_z}{\partial z}+S_{\mathrm{rr}}\frac{{\partial u}_r}{\partial r}+S_{\mathrm{rz}}(\frac{{\partial u}_z}{\partial r}+\frac{{\partial u}_r}{\partial z})-\frac{u_r(S_{\mathrm{rr}}+S_{\mathrm{zz}})}{r}---\left(4\right)$

其中，r和z为空间坐标，ρ为密度，u_z为轴向速度，u_r为径 向速度，e为比内能，P为静水压力，S_rr和S_zz分别为r向和z向的 偏应力。

碎片对航天器的破坏的评估计算就是结合材料的应力应 变关系及状态方程来求解偏微分守恒方程组，按照计算时所 采用的坐标系可分为Euler方法Lagrange方法。超高速碰撞问 题是一个高速瞬态大变形问题，对于该类问题的计算，目前 Lagrange法的研究已经比较成熟。但是，尽管Lagrange法引入 了一些新的算法，如网格重分，在一定程度上解决了网格畸 变等问题，与Euler法相比，仍不可能像Euler法那样自然的反 映大变形。采用Euler方法，当所研究的系统中含有多种物质 的时候，便会出现混合网格，此时混合网格处理方法的优劣 就显得至关重要。如何确定混合网格中的物质界面位置，如 何计算混合网格的力学量以及混合网格向周围网格的物质输 运量，一直是Euler方法中的难题。

混合网格处理方法研究从20世纪五六十年代有了实质性 的发展，其中美国Los Alamos的科学家们提出和发展的格子类 方法(Cell-type method)，应用效果非常好。著名的格子类方法 有PIC(Particle in Cell)方法、MAC(Marker and Cell)方法以及 FLIC(Fluid In Cell)方法等。为了计算多种物质流动的问题， Kershner和Mader提出了2DE方法，以及我国的徐国荣等人提 出了多流体网格方法。1981年，Hirt和Nichols提出著名的VOF (volume of fluid)论文，对于运动界面追踪问题的数值模拟研 究作出了开创性的贡献。目前，比较有代表性的VOF方法有： FLAIR方法及Youngs界面重构方法等。

Youngs于1982年发表了其著名的关于Youngs界面重构技 术的论文。由于其重构精度高且算法简单有效，目前已在国 外一些大型数值仿真软件中得到了广泛的应用，如 AUTODYN、CTH等。其基本思想是将混合网格中的物质界面 近似成一直线，四周8个网格中物质的体积份额用来确定该直 线的法线方向，而混合网格本身的体积份额用来确定直线的 位置。设f为混合格中介质A的体积份额函数，f_E、f_W、f_N和 f_S为网格四边上介质A的体积份额，按下式进行计算：

$f_E=\frac{f_{i-1,j+1}+\alpha f_{i,j+1}+f_{i+1,j+1}}{2+\alpha },$$f_W=\frac{f_{i-1,j-1}+\alpha f_{i,j-1}+f_{i+1,j-1}}{2+\alpha }$

$f_N=\frac{f_{i+1,j-1}+\alpha f_{i+1,j}+f_{i+1,j+1}}{2+\alpha },$$f_S=\frac{f_{i-1,j-1}+\alpha f_{i-1,j}+f_{i-1,j+1}}{2+\alpha }---\left(5\right)$

式中，i、j是网格编号，α为调节系数，通常取α＝2。

对界面斜率k有如下定义：

$\frac{\partial f}{\partial r}=\frac{f_E-f_W}{2\mathrm{\Delta r}},$$\frac{\partial f}{\partial z}=\frac{f_N-f_S}{2\mathrm{\Delta z}}---\left(6\right)$

$k=-\frac{\partial f}{\partial r}/\frac{\partial f}{\partial z}---\left(7\right)$

保持k不变，移动界面线，直到该直线段将网格分割的两 部分体积份额与网格中两物质实际份额相一致，此时直线段 的位置就是所求的物质分界面。在计算出物质界面后，再根据 混合网格的输运体积，便可计算出每种介质相应的输运体积， 进而可以完成质量、动量、能量的输运计算。

质点网格法(particle in cell，PIC方法)，首先由Evan和 F.H.Harlow在1955年提出，Amsden(1966)做了详细总结。该方 法把连续介质(可看作无穷质点系)凝聚成有限的质点系， 质点有质量，参与守恒运算，通过对质点的运算和追踪来实 现流场数值模拟和界面显示。主要的特点是：(1)采用欧拉 网格；(2)拉氏质点是带质量的，可以采用不同的物质符号， 以计算多种物质同时存在的系统；(3)流体状态完全由质点 及分布决定。

PIC方法网格内的物质被看成是若干个带质量的质点，如 图1所示，假设一个网格内分布K_max个质点，那么第k个质点p_k的 坐标(r_k，z_k)可用下面的公式求得：

$r_k=i+(\frac{1}{2}+a)·\frac{\mathrm{\Delta r}}{\sqrt{K_{\max }}}$$a=k\%\sqrt{K_{\max }}$

$z_k=j+(\frac{1}{2}+b)·\frac{\mathrm{\Delta z}}{\sqrt{K_{\max }}}$$b=k/\sqrt{K_{\max }}---\left(8\right)$

其中，Δr和Δz为网格步长。质点p_k的速度u和v由相邻四个 网格的速度按面积加权得到，其计算公式如下：

$\left(\begin{array}{c}u\\ v\end{array}\right)=\frac{1}{\mathrm{\Delta r\Delta z}}\left(\begin{array}{c}{A}_3{\left(\begin{array}{c}u\\ v\end{array}\right)}_{i+1,j+1}\end{array},{+A_2\left(\begin{array}{c}u\\ v\end{array}\right)}_{i,j+1}+A_1{\left(\begin{array}{c}u\\ v\end{array}\right)}_{i,j}+A_4{\left(\begin{array}{c}u\\ v\end{array}\right)}_{i+1,j}\right)---\left(9\right)$

加权面积采用下式计算：

$A_3=[z_k-(j+\frac{1}{2})\mathrm{\Delta z}][r_k-(i+\frac{1}{2})\mathrm{\Delta r}]$

$A_2=[z_k-(j+\frac{1}{2})\mathrm{\Delta z}][(i+\frac{3}{2})\mathrm{\Delta r}-r_k]$

$A_1=[(j+\frac{3}{2})\mathrm{\Delta z}-z_k][(i+\frac{3}{2})\mathrm{\Delta r}-r_k]---\left(10\right)$

$A_4=[(j+\frac{3}{2})\mathrm{\Delta z}-z_k][r_k-(i+\frac{1}{2})\mathrm{\Delta r}]$

若p_k不位于右上角区域，取对应的相邻网格进行加权， 其计算类似。

质点p_k的n+1时刻的坐标为

$r_k^{n+1}=r_k^n+u·\mathrm{\Delta t}$

$z_k^{n+1}=z_k^n+v·\mathrm{\Delta t}---\left(11\right)$

若越过了所在(i，j)网格的格边坐标，那么该质 点就输运了出去，同时带走相应的质量、动量、能量及比能 等。

PIC方法能够通过质点位置的变化来反映计算过程中不 同介质的运动规律，特别适合于模拟材料或结构的断裂、破 碎问题如侵彻问题、高速碰撞问题等，但PIC方法数值波动较 大，计算精度不高。

发明内容

为了解决上述问题，本发明的目的在于提供一种碎片对 航天器破坏作用的粒子评估方法，保证侵彻计算过程中弹靶 接触区域的计算精度。

本发明提供一种碎片对航天器破坏作用的粒子评估方 法，该方法包括以下步骤：

步骤1：建立工程计算模型，确定模型的几何尺寸参数；

步骤2：收集步骤1的计算模型中碎片及航天器材料的物 理参数；

步骤3：根据步骤1的计算模型，以航天器的大小作为计 算域，并采用正交网格离散计算域，网格步长由碎片及航天 器中的最小部件尺寸确定，并在网格内加入PIC粒子；

步骤4：获取步骤3所确定的各网格和PIC粒子影响碰撞毁 伤效果的物理量；

步骤5：连续输出碰撞过程中变化的物理量，采用可视化 技术，直观地得到碎片对航天器破坏效果图。

根据本发明的碎片对航天器破坏作用的粒子评估方法， 不仅适用于碎片对航天器结构破坏的评估，还适用于冲击碰 撞问题的分析评估。本发明在PIC方法的基础上，采用质点映 射输运方法来计算质点在网格间的输运，提高了计算精度， 减小了数值振荡。

附图说明

图1示出PIC方法网格示意图。

图2示出质点p移动后，与相邻网格产生新的重叠影响 区域B₁～B₉。

图3示出实施例中计算模型示意图。

图4示出实施例中试验结果与数值模拟结果的对比。

具体实施方式

下面结合附图，通过具体实施例对本发明进行详细说明。 本发明的特点和优点将随着这些说明变得更为清楚、明确。

在根据本发明的碎片对航天器破坏作用的计算评估方法 中，在步骤1中，建立工程计算模型，确定模型的几何尺寸参 数，即根据工程方案，建立二维轴对称计算模型。

在根据本发明的碎片对航天器破坏作用的计算评估方法 中，在步骤2中，收集上述步骤1的计算模型中碎片及航天器 材料的物理参数，例如包括碎片及航天器材料的初始密度、 状态方程参数、本构方程参数等。

在根据本发明的碎片对航天器破坏作用的计算评估方法 中，在步骤3中，根据步骤1的计算模型，计算域的尺寸取为 碰撞毁伤区宽度的4倍以上，进一步优选地，即计算域的尺寸 取为碰撞毁伤区宽度的4-6倍。采用正交网格离散计算域，网 格步长由碎片及航天器中的最小部件尺寸确定，优选地，保 证最小部件上至少划分5个网格，即网格步长最大为碎片及航 天器中的最小部件尺寸的1/5，进一步优选地，即网格步长为 碎片及航天器中的最小部件尺寸的1/10-1/5，并在网格内加入 PIC粒子。

在根据本发明的碎片对航天器破坏作用的计算评估方法 中，在步骤4中，获取步骤3所确定的各网格和PIC粒子影响碰 撞毁伤效果的物理量，所述物理量例如包括质量、动量、能 量、密度、速度等物理量。

具体而言，获取步骤3所确定的各网格和PIC粒子影响碰 撞毁伤效果的物理量通过采用算子分裂法求解偏微分守恒方 程组来实现。

将质量、动量、能量守恒方程统一写成如下形式：

$\frac{\partial \phi }{\partial t}+u·▿\phi =H---\left(12\right)$

式中，φ代表物理量，例如密度(ρ)、能量(e)或速度 (u)等，表示物理量φ在网格上的变化率；为对流项， 即物理量φ在单位空间网格边界上的通量；等式右边H为源项。

物理效应意义上的算子分裂法将方程(12)分裂成如下两 个方程，同时在数值计算时也分成相应的两步完成：

$\frac{\partial \phi }{\partial t}=H---\left(13\right)$

$\frac{\partial \phi }{\partial t}+u·▿\phi =0---\left(14\right)$

在优选的实施方式中，步骤4包括以下子步骤：

子步骤4.1：求解方程(13)，即Lagrange步，不考虑对流项 的影响，只考虑源项H(压力和偏应力的梯度效应)的作 用，得到网格中各物理量的中间值；

子步骤4.2：在子步骤4.1更新网格和粒子的物理量后，考 虑方程(14)中对流项的影响，即Euler输运步。Euler输运步 是按照质量、动量及能量守恒的原则对质量、动量及能量进 行输运，其中，采用质点映射输运方法，对步骤3中的粒子进 行精确计算。

其中，在子步骤4.1中，具体求解方法为：对方程(13)离 散，时间采用向前差分，空间采用一阶中心差分计算，得到 碰撞过程中网格和粒子的物理量如质量、动量、能量、密度、 速度的更新值。

在子步骤4.2中，本发明采用质点映射输运方法，对步骤3 中的粒子进行精确计算。

子步骤4.2的具体实现过程如下：

(1)采用上述公式(9)、(10)将网格物理量映射到质点；

(2)采用上述公式(11)计算质点新时刻运动到的位置；

(3)经过以上两步，质点运动到了新的位置，改变了质点的 影响域与相邻网格之间的关系，本发明采用质点映射输运方 法将运动后质点的信息再映射回网格。

在质点映射输运方法中，需要记录质点移动之前网格边界 lineⁿ在影响域中的位置，lineⁿ与lineⁿ⁺¹共同划分了质点P的影 响域，构成B₁～B₉。对于图1和图2所示情况，影响域面积关 系有：

A₁＝B₁+B₂+B₄+B₅

A₂＝B₇+B₈                                      (15)

A₃＝B₉

A₄＝B₃+B₆

式中，A₁-A₄为图1所示的质点移动前在各网格内的影响域，B₁-B₉为 图2所示的质点移动后的影响域。

对处于不同位置的质点，其影响域与网格的关系也是不同 的，可以按照类似于上述的对应关系来分析其他情况。分析所 有情况可知，质点影响域与网格总共存在九种位置关系，网格 边界lineⁿ与lineⁿ⁺¹的位置关系又可以分为四种，对于上述各种 位置关系，均可采用如下算法进行计算。

针对每一个质点p，分析影响域B_m(m＝1，2，…，9)与质点移 动前的影响域A_n(n＝1，2，3，4)的对应关系，再将影响域A_n所 分配的物理量U_An按比例分配给影响域B_m，最后将影响域B_m对应的物理量映射回网格。

对于图2所示的情况，质点p影响域与相邻网格的位置关 系为式，质点p影响域A_n，所分配的物理量为U_An，可由式(16) 计算得到：

$U_{A_n}=\frac{A_n}{\Sigma A_{i,j}}{U}_{i,j}$

式中，U_i，j表示影响域A_n所处的(i，j)网格的物理量，∑A_i，j表 示(i，j)网格中所有质点的影响域面积之和。

求得每个影响域A_n的物理量后，按照影响域面积大小把物 理量分配给B_m，并加入到所属网格中，对于图2所示情况，每 个质点对四个相邻网格的输运贡献量为：

${\mathrm{\Delta U}}_{i,j}=\frac{B_1}{A_1}{U}_{A_1}$

${\mathrm{\Delta U}}_{i+1,j}=\frac{B_2}{A_1}{U}_{A_1}+\frac{B_3}{A_4}{U}_{A_4}$

${\mathrm{\Delta U}}_{i+1,j+1}=\frac{B_5}{A_1}{U}_{A_1}+\frac{B_6}{A_4}{U}_{A_4}+\frac{B_8}{A_2}{U}_{A_2}+\frac{B_9}{A_3}{U}_{A_3}$

${\mathrm{\Delta U}}_{i,j+1}=\frac{B_4}{A_1}{U}_{A_1}+\frac{B_7}{A_2}{U}_{A_2}$

对于所有质点，都采用上述过程进行计算。

经过以上三个步骤，通过质点的移动完成了子步骤4.2网 格间的物理量输运。

在根据本发明的碎片对航天器破坏作用的计算评估方法 中，在步骤5中，子步骤4.1和子步骤4.2不断迭代循环，碰撞 过程中变化的物理量得以连续输出。采用可视化技术，可以 直观地得到碎片对航天的器破坏效果图。

实施例

以下通过范例性实施例进一步描述本发明。

该实例是球形碎片高速撞击薄靶。

(1)将球形碎片对薄靶的碰撞视为正碰撞，建立二维轴 对称计算模型，如图3所示。

其中，铝制小球直径R为9.5mm，铝制薄靶的厚度为 2.2mm。

(2)确定球形碎片及靶板的材料参数。具体如下表1中 所示。

表1铝制球形碎片及靶板材料参数列表

  ρ0(g/cm³)   E(GPa)   c₀(m/s)   γ₀  a   s₁  s₂  s₃  2.7   26.4   5100   2   0.43   1.339   0   0

其中ρ₀为材料密度，c₀为材料中的声速，E为弹性模量，γ₀、 a、s₁、s₂、s₃为材料状态方程中的无量纲系数。

(3)确定计算域。

计算域尺寸取为80mm×60mm，根据球形碎片及靶板的尺 寸，两个方向的网格步长均为0.2mm，采用正交网格离散计算 域，并在每个网格内加入9个粒子。

(4)按上述步骤4及其子步骤4.1和4.2，编程实现各网格 和粒子的Lagrange步和Euler输运步的求解。输入球形碎片及 靶板的上述几何参数和材料参数，启动计算。

(5)本实施例共进行了20000步计算，每100步输出一个 可视化数值模拟结果。

图4示出了试验结果与数值模拟结果的对比，其中a表示实 验结果，b表示数值模拟结果。

从图4中可以看出，数值模拟的结果与试验结果是一致 的。具体而言，球形碎片在撞击过程中，受到薄靶的挤压， 快速加载的应力波使得小球头部产生变形、流动，后续卸载 波使得小球向四周扩张、分裂。靶板在受到小球的挤压，靶 板背部凸起，拉伸形变，随后分靶板受撞击部分材料以碎片 云的方式从靶的后表面喷发出来。数值模拟的仿真结果表明， 该发明较好地模拟了超高速碰撞中碎片云的形成，能够自然 地处理计算过程中产生的大量物质界面。表明本发明可以用 于碎片对航天器的破坏评估。

以上通过具体实施方式和范例性实例已对本发明进行详 细说明，不过这些实施方式和实例仅是说明性的，并不对本 发明的保护范围构成任何限制，在不偏离本发明精神和范围 的情况下，本领域技术人员能对本发明及其实施方式进行多 种改进、等价替换或修改，这些均落入本发明的保护范围内。 本发明的保护范围以所附权利要求书为准。