太空中微流星体和空间碎片虽然体形很小，但由于通常都具有超高的相对速度（2～15 km/s），会对航天器的安全造成巨大的威胁。例如，长时间暴露的航天器外部超过30 000个直径大于0.3 mm的陨石坑都是由微流星体或轨道碎片撞击形成^[1]。航天防护中常采用Whipple屏对特征尺寸1 cm以下的碎片进行破碎^[2-3]。Whipple屏在超高速撞击后的损伤状态对于Whipple屏的优化设计至关重要。

超高速撞击的物理过程非常复杂。撞击过程中，材料的惯性、可压缩性效应或相变效应比结构效应更显著，伴随着大变形、热流固耦合、相变（液化、气化、凝固等）、碎裂等现象，采用数值方法进行准确模拟是一项巨大的挑战^[4]。基于网格的超高速撞击模拟存在着网格畸变或者需要不断重新划分网格^[5]的缺点；而无网格法由于不需要进行网格离散及采用高阶插值函数，有利于解决大变形和流固耦合问题，在超高速撞击问题上应用更为广泛。目前常用的无网格方法有光滑粒子流体动力学法（smoothed particle hydrodynamics, SPH）和再生核粒子法（reproducing kernel particle method, RKPM）以及质量点法（material point method, MPM）等。有大量的研究应用以上的方法^[6-7]及其改进方法^[8]对Whipple屏的超高速撞击进行模拟，虽然取得了丰富研究成果，但由于算法本身的固有缺陷，仍然无法解决准确设置位移边界条件^[9]、计算效率不高^[10]、拉应力不稳定^[11]及有效计算带摩擦的动态接触等问题，同时也缺乏严密的收敛性与误差理论分析。

最优运输无网格方法（optimal transportation meshfree, OTM）是由Li等^[12]提出的一种基于最优运输理论和局部最大熵插值函数的无网格方法，该方法结合基于能量释放率的物质点失效方式，能够很好地解决现有Whipple屏超高速撞击模拟方法存在的问题。

本文首先采用OTM法模拟铝球超高速撞击单层铝板，通过与实验结果以及其他数值方法的计算结果的比较，验证OTM法在超高速撞击问题上的适用性；然后采用OTM法对Whipple屏超高速撞击问题进行模拟，研究不同速度不同撞击角度下的碎片云形状、缓冲墙弹孔尺寸以及后墙剥落穿透的损伤情况，并与文献[13]的实验结果进行对比。

1.   OTM理论

OTM 法的主要特点是采用局部最大熵无网格插值函数，克服了一般无网格法中插值函数不满足Kronecker delta 属性的本质缺陷，解决了传统无网格法难以准确施加位移边界条件的问题；采用物质点充当积分点，有效避免了计算结果在拉伸载荷下的不稳定性；采用最优运输理论对时间离散，保证了哈密顿作用量的时间离散形式满足动量守恒条件，且收敛性能得到严格的数学证明^[14]；采用基于能量释放率的物质点失效方式，能够很好地模拟材料的损伤情况，并且该方式已被证明可收敛到 Griffith 断裂准则^[15]。

1.1   时间离散

与时间相关的问题可以转化为一系列以最小化原理为特征的递增问题。OTM法采用最优运输理论离散流体或固体粒子流动的动能。最优运输理论通过最小化某一时间间隔$\left[{t}_{k},{t}_{k+1}\right]$内物体运动的哈密顿作用量$S\left[\rho ,v\right]$来得到$\left[\rho ,v\right]$(ρ和v分别为物体的密度和速度)，描述了质量运输问题的真实运动过程。时间间隔$\left[{t}_{k},{t}_{k+1}\right]$中初始质量密度${\rho }_{k}$和最终质量密度${\rho }_{k+1}$之间的Wasserstein距离给出了动能增量的确切最小值。Wasserstein距离d_w的定义为：

式中：T是对于${\rho }_{k}$到${\rho }_{k+1}$任何可能的映射。

设${t}_{0}=a{\text{＜}}{t}_{1}{\text{＜}}{t}_{2}{\text{＜}}\cdots {\text{＜}}{t}_{k}{\text{＜}}{t}_{k+1}{\text{＜}}\cdots {\text{＜}}{t}_{n-1}{\text{＜}}{t}_{n}=b$是时间间隔$\left[a,b\right]$的离散化，并假设在时间${t}_{k}$时，${\rho }_{k}$、${\varphi }_{k}$是已知的，在时间${t}_{k+1}$处的状态${\varphi }_{k+1}$可以通过求极限的方式近似得到。差分(${\varphi }_{k+1}$−${\varphi }_{k}$)/${\Delta }t$的极限即${\Delta }t={t}_{k+1}-{t}_{k}$趋于零时可以表示速率。因此采用Wasserstein距离表示的时间离散后的动能项为：

由此可得到时间离散后的哈密顿作用量形式

式中：$U\left(\varphi \right)$为自由能，代表与体积力、运动约束或边界压力相关的能量，其表达式为：

式中：B为积分区域，$f(\nabla \varphi )$是单位体积上的增量自由能密度，${\varphi }_{k\to k+1}$为增量传输映射，详见1.2节。

最优运输理论为OTM法中显式的时间积分算法提供了可靠的理论基础，如果将时间视为一个独立的广义坐标并且采用连续的时间步长，总能量也能精确守恒^[16]。

1.2   空间离散

OTM法将计算域分为两类点的集合。一类为节点，用下标a表示，包含了位置、速度、加速度的信息；一类为物质点，用下标p表示，包含了材料参数、质量、应力应变等信息，相当于一般有限元法中的积分点。定义在${t}_{k}$时刻节点的坐标为${x}_{a,k}$，物质点的坐标为${x}_{p,k}$，如图1所示。

图  1  空间离散示意图^[17]

Figure  1.  Spatial discrete diagram^[17]

下载: 全尺寸图片 幻灯片

从${t}_{k}$时刻到${t}_{k+1}$时刻的位置迁移由变形映射函数${\varphi }_{k\to k+1}\left(x\right)$表示

式中：$M$是${x}_{p,k+1}$邻域N_h(x_p,k)范围内节点的个数；插值函数${N}_{{a},k}\left(x\right)$为局部最大熵形函数，它是考虑了局部熵最大（也就是最随机的状态）、同时又考虑了邻域最小的优化问题的解^[14]，其表达式为

其中

式中：$\beta$是Pareto最优参数；与局部熵有关；${\lambda }^{*}\left(x\right)$是最优拉格朗日乘子；插值函数${N}_{{\rm{a}},k}\left(x\right)$的约束条件为

约束条件确保了所得到的形函数满足零阶和一阶连续性。

空间离散后可得到完全离散的作用量形式

1.3   失效准则

OTM法中通过删除失效物质点的方式来模拟裂纹的扩展以及破碎的形成。当某个物质点失效时，应力计算等过程中会将其忽略。对于失效的判定，OTM法采用能量释放率${G}_{p,k+1}$作为其判定准则。能量释放率${G}_{p,k+1}$通过局部能量平均进行预估，定义如下

式中：${B}_{\varepsilon }\left({x}_{p,k+1}\right)$是以${x}_{p,k+1}$为中心、$\varepsilon$为半径的球，如图2所示；${F_{{\rm{q}},k + 1}}$为变形梯度；${m}_{p,k+1}$是$\mathrm{\varepsilon }$邻域的质量：

图  2  评估能量释放率的局部邻域^[15]

Figure  2.  The local neighborhood used to estimate the energy-release rate^[15]

下载: 全尺寸图片 幻灯片

当某物质点的能量释放率(${G}_{p,k+1}$)大于临界能量释放率(${G}_{\mathrm{c}}$)，认为该物质点失效。临界能量释放率${G}_{c}$刻画了单位面积下材料表面发生破碎所需要的能量。

2.   铝球超高速撞击铝板算例分析

文献[8]采用基于拟流体模型的SPH新方法（拟流体SPH法）对铝球超高速撞击铝板进行了数值模拟，并与Hiermaier等的实验结果、传统SPH法和自适应光滑粒子流体动力学法（adaptive smoothed particle hydrodynamics, ASPH）的计算结果进行了对比。本节将采用OTM法对相同的算例进行模拟，并与文献[8]中的结果进行比较。

2.1   数值模拟模型

算例为直径10 mm的铝球和尺寸40 mm$\times$40 mm$\times$2 mm的铝靶板的超高速碰撞，撞击速度为6.18 km/s。OTM法的空间离散方式为：对几何模型划分网格，将网格节点作为OTM法中的节点，单元中心作为物质点。本算例中网格划分尺寸为0.67 mm，离散结果如图3所示，共有22 094个节点，100 484个物质点。

图  3  铝球撞击铝板离散模型

Figure  3.  Discrete model of aluminum ball impacting single aluminum plate

下载: 全尺寸图片 幻灯片

2.2   材料模型

超高速撞击需要考虑弹丸和靶材的应变率硬化和热软化等问题，因此选择能较好地模拟铝合金材料塑性响应的J2粘塑性模型（J2-viscoplasticity）。

J2粘塑性模型有效屈服应力为

式中：${\varepsilon }^{\rm{p}}$和${\dot{\varepsilon }}^{\rm{p}}$为有效塑性应变和应变率，$T$为绝对温度，${\sigma }_{0}$为参考温度下的初始准静态屈服应力，${\varepsilon }_{0}^{\rm{p}}$和${\dot{\varepsilon }}_{0}^{\rm{p}}$为参考有效塑性应变和应变率，${T}_{\rm{m}}$为熔化温度，${T}_{0}$为参考温度，$n$为强化指数，$m$为应变率指数，$q$为热软化指数，熔化温度${T}_{\rm{m}}$与变形相关，表达式为

式中：J=V/V₀，V为当前体积，V₀为初始体积；${T}_{{\rm{m}}0}$为参考熔化温度；${\gamma }_{0}$为参考体积下的Grüneisen参数，$a$为常数。铝合金型号为LY12，材料参数如表1所示。本算例中依据文献[8]提供的材料信息进行了部分参数修改。式(15)和式(16)中提到的各参考系数的取值如表2所示。高温高压下材料的变形与温度、压力的关系采用SESAME状态数据库描述^[18]。

表  1  LY12材料参数

Table  1.  Material parameters of LY12

密度/(kg∙m−3)	弹性模量/GPa	泊松比	比热容/(J∙kg−1∙K−1)
2 700	68.9	0.33	896

下载: 导出CSV 

| 显示表格

表  2  J2黏塑性模型参数

Table  2.  Parameters of J2 viscoplasticity model

σ0/MPa	${\varepsilon }_{0}^{\rm{p}}$	${\dot{\varepsilon } }_{0}^{\rm{p} }$	n	m	q	Tm0/K	$a$	${\gamma }_{0}$
276	$5\times {10}^{-4}$	1 000	0.075	0.08	0.5	925	1.5	1.97

下载: 导出CSV 

| 显示表格

2.3   结果分析

OTM法能够很好地模拟内核碎片云的位置、外泡碎片云的形态以及反溅碎片云的形态等特征信息。特别是反溅碎片云的膨胀距离和宽度，与实验很好吻合。一般的SPH法由于采用Johnson-Cook损伤模型而导致薄板屈服应力小于真实的屈服应力，外溅碎片云的反溅程度过大，与实验偏差较大（如图4所示）。

图  4  OTM法与各类SPH方法计算结果对比

Figure  4.  Comparison of OTM and various SPH methods’ simulation results

下载: 全尺寸图片 幻灯片

表3展示了相关实验和计算方法在弹孔直径$d$（不包括弹孔边缘）、弹孔直径误差$\varepsilon$、碎片云膨胀距离$l$，碎片云膨胀宽度$w$，以及$l/w$及其相对误差的对比结果。从该表中可以看出，OTM法预测的弹孔直径与真实实验中的弹孔直径最为相近，误差仅为4.7%；碎片云的膨胀距离与宽度之比的相对误差也是各类方法中最小的，说明了OTM法能很好地模拟出碎片云形态。

表  3  铝球超高速撞击铝板结果对比

Table  3.  Comparison of high-velocity impact results between aluminum projectile and plate

方法	d/mm	ε/%	l/mm	w/mm	l/w	Δ/%
Hiermaier实验	27.5	−	−	−	1.39	−
Hiermaier模拟	35.0	27.3	−	−	1.11	−
SPH法	31.6	14.9	102.8	75.5	1.36	2.2
ASPH法	28.9	5.1	105.1	86.1	1.22	12.2
拟流体SPH法	29.4	6.9	105.7	81.4	1.30	6.5
OTM法	26.2	4.7	104.2	76.7	1.36	2.2

下载: 导出CSV 

| 显示表格

3.   Whipple屏超高速撞击模拟

3.1   对比实验简介

文献[13]进行了一系列Whipple屏超高速撞击实验（如图5所示），球形弹丸直径0.4～0.5 cm，撞击速度4.47～6.15 km/s，撞击角度分为0°和45°两种。靶材间距为10 cm，厚度为0.192 cm。实验得到不同撞击速度和撞击角度下的弹孔尺寸、后墙损伤情况和碎片云激光阴影照片等结果。

图  5  实验模型示意图

Figure  5.  Schematic diagram of experimental model

下载: 全尺寸图片 幻灯片

3.2   数值模拟模型

本文按照对比实验建立模型，模拟了撞击角度为0°和45°两种情况。在划分网格时细化了缓冲墙和后墙的中心区域，最小网格尺寸为0.2 mm。在正撞模拟中，物质点共有340 382个，节点共有67 389个（如图6所示）；在斜撞模拟中，物质点共有340 473个，节点共有67 401个。

图  6  Whipple屏超高速撞击数值模拟模型

Figure  6.  The numerical simulation model of Whipple shield hypervelocity impact

下载: 全尺寸图片 幻灯片

弹丸和靶材的材料型号均为LY12，材料模型仍然使用J2粘塑性模型（材料参数见表1、表2）和SESAME状态方程。

3.3   结果分析

文献[13]进行了一系列不同撞击速度的Whipple屏正撞与斜撞实验。参数设置如表4所示。

表  4  实验参数设置

Table  4.  Parameters in experiments

实验	弹丸直径/mm	弹丸质量/g	缓冲墙厚度/mm	后墙厚度/mm	撞击速度/(km·s−1)	撞击角/(°)
04-0090	5.00	0.1797	1.92	1.94	5.29	0
04-0092	5.02	0.1826	1.94	1.90	5.52	0
04-0079	5.00	0.1810	1.94	1.92	6.08	0
04-0080	5.00	0.1811	1.92	1.90	6.15	0
04-0084	4.04	0.0972	1.92	1.90	5.95	45
04-0083	4.02	0.0960	1.92	1.94	6.02	45
04-0075	4.02	0.0958	1.92	1.90	4.47	45
04-0077	4.00	0.0940	1.92	1.94	4.74	45

下载: 导出CSV 

| 显示表格

OTM法采用相同的参数进行对应的数值模拟，实验与仿真中的缓冲墙弹孔尺寸对比结果如表5所示。其中后四组是斜撞实验，得到的弹孔呈现椭圆形，弹孔尺寸指的是椭圆的长轴和短轴。

表  5  缓冲墙弹孔尺寸对比

Table  5.  Bullethole size comparison of outer bumper

实验	撞击速度/(km∙s−1)	实验缓冲墙弹孔尺寸/mm	仿真缓冲墙弹孔尺寸/mm	相对误差
04-0090	5.29	11.5	10.5	8.69%
04-0092	5.52	11.7	10.9	6.84%
04-0079	6.08	12.4	11.2	9.68%
04-0080	6.15	12.6	11.8	6.35%
04-0075	4.47	10.6×8.5	10.9×8.99	2.83%×5.76%
04-0077	4.74	10.6×8.7	11.2×9.39	5.66%×7.93%
04-0084	5.95	11.6×10.2	12.3×10.1	6.03%×0.98%
04-0083	6.02	11.6×10.3	12.1×9.75	4.31%×5.34%

下载: 导出CSV 

| 显示表格

缓冲墙弹孔尺寸的模拟结果与实验吻合得较好，如在撞击速度5.29 km/s的实验中，弹孔直径为1.15 cm，而仿真的结果为1.05 cm，相对误差为8.69%（如图7所示）。

图  7  缓冲墙损伤对比图(撞击速度5.29 km/s)

Figure  7.  Damage characteristics comparison chart of outer bumper (Impact velocity 5.29 km/s)

下载: 全尺寸图片 幻灯片

后墙的损伤形式一般有成坑、产生鼓包、层裂、剥落和穿孔^[19]。文献[13]只关注剥落和穿透，如图8(a)所示。图8(b)为OTM仿真中的剥落和穿透。

图  8  实验与仿真中的剥落和穿透

Figure  8.  Definitions of spalling and penetration in experiments and simulations

下载: 全尺寸图片 幻灯片

对比仿真与实验中后墙损伤情况（如表6所示），可见正撞仿真中，撞击速度越大，后墙的损伤越小。这是由于速度大的弹丸被缓冲墙破碎得更充分，形成了更小的碎片，减轻了对后墙的破坏作用。这与文献[13]中的结论一致。

表  6  后墙损伤情况对比

Table  6.  Damagecomparison of spacecraft wall

实验	撞击速度/(km∙s−1)	实验后墙损伤情况	仿真后墙损伤情况
04-0090	5.29	3处剥落，无穿透	无剥落，2处穿透
04-0092	5.52	2处剥落，无穿透	5处剥落，2处穿透
04-0079	6.08	无剥落，无穿透	无剥落，无穿透
04-0080	6.15	无剥落，无穿透	无剥落，无穿透
04-0075	4.47	无剥落，2处穿透	无剥落，无穿透
04-0077	4.74	无剥落，无穿透	无剥落，无穿透
04-0084	5.95	1处剥落，1处穿透	无剥落，无穿透
04-0083	6.02	1处剥落，1处穿透	1处剥落，无穿透

下载: 导出CSV 

| 显示表格

仿真显示，弹丸碎片撞击到后墙上，残余应力约为150 MPa，而撞击较严重的区域残余应力达到300 MPa以上，如图9所示。

图  9  后墙损伤图(撞击速度5.29 km/s)

Figure  9.  Damage characteristics of spacecraft wall (impact velocity is 5.29 km/s)

下载: 全尺寸图片 幻灯片

目前超高速弹丸与Whipple防护屏撞击的数值模拟，很难精确模拟出后墙的损伤情况。大部分后墙的损伤数据都来自于实验。文献[8]用拟流体SPH法对Whipple屏撞击的研究中，也只给出了后墙的中心损伤区域，缺乏对于后墙剥落与穿透等损伤状态的探究。由此可看出，OTM法基于能量释放率的物质点失效的方式，在模拟材料的断裂损伤方面有着明显的优势。

OTM法对Whipple屏的超高速撞击模拟不仅能模拟出弹丸穿透缓冲墙形成碎片这一过程，还能很好地模拟出碎片云的形态，包括碎片云呈现出头部椭圆形、尾迹扇形并有扩散趋势、碎片大多集中在头部和缓冲墙穿孔处等细节特征，如图10～11所示。由于实验中的碎片云照片是采用激光阴影技术得到的，细微的碎片无法在照片中显示，因而实验中的碎片云轮廓较为清晰；而OTM计算得到的碎片云对比图中，所有碎片均有显示，轮廓有一定发散性。

图  10  正撞碎片云对比图

Figure  10.  Fragment cloud comparison chart of vertical impact

下载: 全尺寸图片 幻灯片

图  11  斜撞碎片云对比图

Figure  11.  Fragment cloud comparison chart of oblique impact

下载: 全尺寸图片 幻灯片

4.   结　论

本文采用OTM法对铝球超高速撞击铝板和Whipple防护屏超高速撞击进行了数值模拟。通过铝球超高速撞击铝板这一验证算例，可看出OTM法能够为超高速撞击问题提供有效的数值模拟手段。在Whipple屏的超高速撞击模拟中，OTM法能够很好地预测Whipple防护屏与弹丸撞击时缓冲墙和后墙的损伤情况，尤其在模拟后墙的剥落、穿透等损伤状态和破碎过程中碎片云的形态变化方面有着显著的优势，验证了其在超高速撞击数值模拟方面的适用性，为Whipple防护屏在航天防护方面的相关探究提供了有效的数值模拟手段。