随着航天技术的高速发展，大量空间碎片遗留在地球轨道并处于高速运动状态，时刻威胁着在轨航天器的安全^[1]。在航天器舱壁外加装防护屏的被动防护技术作为一种有效的防护方式越来越受到重视，这种结构最先由F.L.Whipple^[2]提出，防护机理为空间碎片超高速碰撞防护屏后，生成大量更为细小的碎片和颗粒，少部分碎片向相反方向运动，大部分碎片向运动前方抛出。目前，研究碎片运动的方法主要包括实验，理论计算以及数值模拟方法。其中，A.J.Piekutowski^[3]和管公顺等^[4]进行了大量的实验，得出了弹丸直径、薄板厚度以及碎片云各特征点速度的变化规律，但研究存在着成本高、周期长、得到的实验数据有限等问题；H.F.Swift等^[5]在实验的基础上，作出若干适当假设，结合经验公式对碎片云特性进行研究，提出了球形弹丸高速碰撞薄板的碎片云模型，思想具有开创性意义，模型虽然不断进行细节改进，但仍存在着对碎片云描述过于简化，与实际情况差距较大的局限性；数值模拟具有快速、直观的特点，不仅能弥补实验研究的不足并为实验提供指导，还可以充分了解超高速碰撞全过程和内在机理。目前，基于网格的计算方法对弹塑性材料中低速变形比较可靠，但超高速碰撞问题涉及了移动的边界和运动的物质交界面，当采用Lagrange有限元方法计算超高速碰撞问题时，需要引入网格重分和界面重构等算法进行处理，将引入难以忽略的数值误差，且算法复杂，计算量大。虽然有限差分(FD)、有限体积(FVM)等Euler方法不存在计算大畸变问题的困难，但是难以准确描述各类界面，并且无法追踪单碎片的运动轨迹。

光滑粒子流体动力学方法(smoothed particle hydrodynamics, SPH)^[6-7]作为一种完全Lagrange无网格粒子法，相对于传统网格方法具有精确追踪自由表面、变形边界和材料交界面等优势，非常适合模拟穿透、层裂、剥落及碎片等现象，在处理多介质、大变形破坏问题上具有网格方法无法比拟的优势，被广泛应用于高速冲击，爆炸和断裂等问题。其中，L.D.Libersky等^[8]率先将材料强度效应引入SPH方法，成功地开展了高速碰撞数值模拟计算。G.R.Johnson等^[9]将SPH方法和有限元方法结合开展侵彻方面的数值模拟，取得一些进步。S.Hiermaier等^[10]在对铝球超高速碰撞不同材料的薄靶实验的基础上，使用传统SPH方法对其作了数值模拟，得到的弹坑直径、碎片云形状等结果与实验结果相似。徐金宏等^[11]将改进的SPH方法应用于超高速碰撞问题，研究了靶板厚度、弹丸速度和形状对碎片云的影响。C.E.Zhou等^[12]使用传统SPH方法对铝球超高速撞击铝薄板的实验进行了三维数值模拟。G.R.Liu等^[13]使用ASPH方法求解弹丸超高速碰撞薄板问题，相比传统SPH方法增加了空间求解能力。综上所述，SPH算法经过不断完善与发展，成为近年来超高速碰撞数值计算中一种非常有吸引力和有前途的新的数值方法。但目前国内外学者在求解超高速碰撞问题时，对碰撞后的物质大多采用Grüneisen和Tillotson状态方程进行描述，通过将“碎片云”中颗粒的运动类比于“流体”的行为进行相态划分，仿真结果仅能达到形态上相似、定量相近求解精度较低，无法得到碰撞后碎片运动理想效果。同时，由于碎片云在膨胀过程中碎片与碎片之间以及碎片与空气间的相互作用起着至关重要的作用，碎片表现出典型的离散颗粒的性质。而传统SPH方法大多仅用于离散求解连续介质，对于与连续介质完全不同的离散颗粒来说传统求解连续性物质的数值方法不再有效。因此对于弹丸超高速碰撞形成碎片云的数值模拟需要采用适用于离散颗粒相求解的理论进行补充研究。

目前基于拟流体模型的用于描述气体-颗粒两相流动问题的双流体模型理论逐渐成熟，陈福振等^[14]已成功将颗粒动理学理论引入到传统SPH方法之中，建立了适用于气粒两相流问题求解的光滑离散颗粒流体动力学方法(smoothed discrete particle hydrodynamics, 简称SDPH)。该方法不仅突破了传统SPH方法材料属性固定、粒子仅作为几何质点的局限，将SPH粒子拓展至具有物理质点特性的范畴，还拓展了SPH方法的应用范围，现已成功应用于风沙运动过程^[14]、气-粒两相间传热蒸发^[15]和燃料抛撒成雾及其燃烧爆炸^[16]等问题。本文中在此基础上，将新方法引入到超高速碰撞领域，建立基于拟流体模型的SPH新方法，将超高速碰撞后产生的物质离散为SPH粒子，将处于损伤状态的颗粒相等效于拟流体，同时在描述其运动过程中引入气体相的作用，通过施加空气外力和计算碎片阻力的方式用来表征空气对碎片云的影响；对球形弹丸超高速碰撞薄靶板和Whipple防护结构进行三维数值模拟，分析靶板开坑直径、碎片云膨胀过程、形态及不同撞击速度下Whipple结构的损伤特性。

1.   数学模型

在球形弹丸超高速碰撞铝薄板数值模拟过程中，未处于失效状态的铝球和铝薄板属于SPH方法中连续相的离散求解问题。超高速碰撞后处于失效状态的碎片属于离散颗粒相求解问题，因此本文在Johnson-Cook损伤模型的基础上引入SDPH模型，采用适用于离散颗粒相求解的方法求解碎片的运动过程，将碎片等效为拟流体，引入碎片间相互作用和气体相对碎片的作用。

1.1   SPH基本方程

在SPH方法中，采用粒子近似方法对核函数进行插值离散并对粒子i处场函数进行粒子估计有：

式中：i、j为里编号，N为紧支域内的粒子总数，f为任意的连续光滑函数，〈〉为近似符号，m_i、ρ_i和r_i分别表示粒子i的质量、密度和位置矢量，W_ij=W(r_i-r_j, h)为核函数，h为定义核函数影响区域的光滑长度。对式(1)求导可得：

未处于失效状态的铝球和铝薄板属于SPH方法中连续相的离散求解问题，不考虑热传导过程，Lagrange框架下的控制方程组表述为：

式中：ρ、v、P、x分别为流体的密度、速度、压力、位置，τ为偏应力张量，g为质量体积力。

最后，考虑人工粘性的影响，对式(3)~(5)进行SPH离散，得到适用于超高速碰撞问题的SPH方程组：

式中：v_ij=v_i－v_j，$\mathit{\boldsymbol{\sigma }} = - \nabla P + \nabla \cdot \mathit{\boldsymbol{\tau }}$，Π为人工黏性，α和β表示坐标方向。

1.2   弹丸和靶板的本构模型及状态方程

为描述靶板材料的屈服应力及损伤演化，这里引入Johnson-Cook损伤模型^[17]，该模型中将材料的屈服强度表示为损伤变量、等效应变、等效应变率和温度的函数：

式中：A、B、C、n、m是材料常数，D为损伤变量，D=0表示材料没有损伤，D=1表示材料完全失效。r是累积损伤塑性应变，r=(1－D)ε，${\dot r^ * } = \left( {1 - D} \right){\dot \varepsilon ^ * }, {\dot \varepsilon ^ * } = \dot \varepsilon /{\dot \varepsilon _0}$，ε是累积塑性应变，${\dot \varepsilon _0}$是自定义参考应变率。温度T^*=(T-T₀)/(T_m-T₀)，T₀是室温，T_m是材料熔点。

损伤变D是累积塑性应变ε的函数，当D=1时发生损伤破坏：

式中：ε_f是断裂塑性应变，与材料的应力三轴度、应变率和温度相关。本构模型中的剪切损伤演化模型将ε_f描述如下：

式中：D₁~D₅为材料常数，σ^*=σ_m/σ_eq为应力三轴度，σ_m=(σ_x+σ_y+σ_z)/3为平均正应力。

超高速碰撞问题中最常用的两个状态方程为Tillotson状态方程和Grüneisen状态方程，Tillotson状态方程形式如下：

式中：P为压力；a、b、A_T为拟合常数；B_T为调节参数；α_T、β_T为材料常数；E为比内能；E₀为初始比内能；V为比容；μ=ρ/ρ₀－1，η=ρ/ρ₀，ρ为材料的实时密度，ρ₀为材料的初始密度。

Grüneisen状态方程形式如下：

式中：系数Γ=1.99，C_S=3940，S_S=1.489；ρ为材料的实时密度；ρ₀为材料的初始密度；e为材料的内能。

1.3   SDPH模型方程

超高速碰撞后处于失效状态的碎片属于离散颗粒相求解问题。本节中列出了SDPH模型中颗粒相质量、动量、拟温度守恒方程以及颗粒相的SPH离散方程等，具体方法及各参数表达式参见文献[14]。

同时，碎片区别于具有不同粒径、任意分布的矿石、风沙等颗粒，相对于风沙运动过程等二维模型中单位面积的颗粒数，球形弹丸超高速碰撞铝薄板产生的碎片数目相对较少、计算量较小，本文SDPH单粒子不必为减少计算量而表征具有一系列粒径分布的颗粒群。由于SDPH单粒子与单个碎片属于几何近似关系，因此SDPH单粒子代表一个真实碎片，携带质量、密度、速度、拟温度以及压力等参量。但此时的SDPH粒子不等同于传统SPH方法中用于求解连续相流体问题的连续性介质的一个质点粒子，而是作为离散颗粒的形式用来求解拟流体模型和描述颗粒性质的变化历程。

固体颗粒相的控制方程中，连续性方程、动量方程、拟温度守恒方程分别为：

式中：ρ_s和v_s分别为颗粒相密度和速度，P为连续相压力，P_s为离散相压力，τ_s为颗粒相粘性应力张量，ρ_sg为外部体积力，R_sg为相间相互作用力；T_s为颗粒相拟温度，I为单位矩阵，$\left( { - {p_{\rm{s}}}\mathit{\boldsymbol{I + }}{\mathit{\boldsymbol{\tau }}_{\rm{s}}}} \right):\nabla {\mathit{\boldsymbol{v}}_{\rm{s}}}$为由颗粒相应力产生的能量，${k_{{T_{\rm{s}}}}}\nabla {T_{\rm{s}}}$为能量耗散项，k_Ts为能量耗散系数，γT_s为颗粒间碰撞产生的能量耗散项，${\phi _{{\rm{gs}}}}$为连续相与颗粒相间的能量交换，具体参数表达式详见文献[14]。

固体颗粒相SPH离散方程组中，连续性方程、动量方程、拟温度守恒方程分别为：

式中：应力σ=-p_sI+τ_s，ρ_i为SDPH粒子i的密度(即颗粒相有效密度)；ρ_s为颗粒的实际密度；速度矢量v_ij=v_i-v_j，R′_gs^sph为作用于SPH粒子上的单位质量曳力，f_i^bp为作用于粒子i上的壁面力；拟温度梯度$\nabla {{\rm{T}}_s}$的SPH离散公式为：

式中：ρ_ij=(ρ_i+ρ_j)/2。

2.   算例验证分析

本章算例中室温及周围环境的初始温度均设为273K，时间积分采用蛙跳格式，时间步长为0.1μs。假设碰撞过程中，靶板不受侧面边界约束，光滑长度采用1.5倍粒子间距。

2.1   铝球超高速碰撞铝薄板

因为航天器的防护结构多为铝合金，所以本文参照Hiermaier等^[10]的实验中铝球超高速碰撞铝薄板的过程进行数值计算，并结合Zhou等^[12]和Liu等^[13]采用传统SPH方法及ASPH方法的计算结果进行对比验证。

如图 1所示，铝球直径为10mm，铝靶板尺寸为40mm×40mm×4mm。铝球的撞击速度为6.18km/s。SPH粒子直径为0.67mm，整个模型共有23391个粒子，其中铝球离散为1791个粒子，铝薄板离散为21600个粒子。人工粘性参数同Hiermaier^[10]和Zhou保持一致。由于Hiermaier并没有列出所用铝合金型号，本文根据密度、屈服应力等强度模型参数确定其类别，其他参数参考同类材料^[18-19]，具体材料参数参见表 1，其中：C_V为定容比热。状态方程选用同Hiermaier等^[10]和Zhou等^[12]相同的Tillotson状态方程，具体参数见表 2。

图  1  算例模型结构

Figure  1.  Experimental model

下载: 全尺寸图片 幻灯片

表  1  本构模型参数

Table  1.  Parameters of constitutive model

A/MPa	B/MPa	n	C	m	Tm/K	${\dot \varepsilon _0}$/s-1
300	426	0.34	0.015	1.0	775	1.0
CV/(J·kg-1·C-1)		D1	D2	D3	D4	D5
875		0.12	0.13	-1.5	0.0175	0

下载: 导出CSV 

| 显示表格

表  2  状态方程参数

Table  2.  Parameters of state equation

ρ/(kg·m-3)	a	b	AT/GPa	BT/GPa	αT/GPa	βT/GPa	E0/(MJ·kg-1)
2790	0.5	1.63	75	65	5	5	5

下载: 导出CSV 

| 显示表格

图 2为采用SPH新方法计算得到的碎片分布与Hiermaier实验^[10]、Hiermaier二维数值模拟结果^[10]、Zhou等采用传统SPH方法的三维数值模拟结果^[12]和Liu等采用ASPH方法的二维数值模拟结果^[13]对比图。从图 2中可以看出，采用SPH新方法和传统SPH方法中的内核碎片云均位于外泡碎片云的前方且占了整个碎片云质量的大部分。但通过对比，可以发现存在两处不同：

图  2  不同SPH方法数值模拟结果与实验结果对比

Figure  2.  Comparison of different SPH algorithm's simulation results and the experimental result

下载: 全尺寸图片 幻灯片

第一处不同为采用传统SPH方法的内核碎片云成圆锥型，随着外泡碎片云的不断膨胀锥面不断弯曲，分布比较集中。而采用SPH新方法的将处于失效状态的碎片等效于拟流体，充分考虑碎片间相互作用及数值模拟中经常被忽略的气体相对碎片的影响，故内核碎片云呈“圆球”型，分布比较分散，与实验吻合较好；

第二处是采用SPH新方法的计算结果中反溅碎片云比采用传统SPH方法的反溅碎片云膨胀的距离和宽度都小，这是因为SPH方法引入Johnson-Cook损伤模型后，薄板屈服应力要小于未引入损伤模型的屈服应力，如式(9)所示；根据反溅碎片云形成过程研究^[20]，当铝球与薄板接触点后移速度高于薄板材料反溅速度向后的分量时，铝球阻挡薄板该部分材料继续反溅，反之无法阻止薄板材料反溅；而薄板的屈服应力直接决定了其反溅速度，所以采用SPH新方法的计算结果中反溅碎片云比采用传统SPH方法中的反溅碎片云膨胀的距离和宽度都小；但从反溅碎片云形态来看，新方法计算得到的反溅碎片云运动形态与实验更为接近，而传统SPH方法中反溅碎片云的计算结果明显过大。

表 3给出了采用SPH新方法得到的弹坑直径、外泡碎片云和内核碎片云的形状、分布及其与使用传统SPH方法^[10]、ASPH方法的模拟结果^[13]的对比，其中：d₁和d₂分别为不包括和包括弹坑边缘的弹坑直径，l为碎片云的膨胀距离，w为碎片云宽度，Δ为l/w的相对误差。从表 3可以看出，采用SPH新方法得到的结果与实验吻合较好。与使用传统SPH方法、ASPH方法的模拟结果进行对比，新方法在内核碎片云形状和分布上计算结果更加准确，与实验更加符合。表明SPH新方法适合模拟铝球冲击铝薄板等超高速碰撞问题。

表  3  超高速碰撞数值模拟结果对比

Table  3.  Comparison of high-velocity impact simulation results

方法	D1	D2	l/mm	w/mm	l/w	Δ/%
SPH新方法	29.4	33.7	105.7	81.4	1.30	6.5
Hiermaier实验[10]	27.5	34.5	-	-	1.39	-
Hiermaier模拟[10]	35.0	-	-	-	1.11	20.1
Liu模拟[13]	28.9	-	105.1	86.1	1.22	12.2
Zhou模拟[12]	31.6	35.3	102.8	75.5	1.36	2.2

下载: 导出CSV 

| 显示表格

为形象展示了碎片云形成、膨胀的过程，图 3~4分别以俯视、左视两种视角展示了不同时刻铝球超高速撞击铝薄板三维数值模拟过程，其中图 4以左视角展示了“碎片云”形成、发展及最终形态，并与Zhou采用传统SPH方法的三维数值模拟结果^[12]进行对比。

图  3  铝弹超高速碰撞铝薄板俯视图

Figure  3.  Top view of hypervelocity impact of projectile on thin plates

下载: 全尺寸图片 幻灯片

图  4  SPH新方法与传统SPH方法计算结果对比

Figure  4.  Comparison of the new SPH algorithm and traditional SPH algorithm's simulation results

下载: 全尺寸图片 幻灯片

2.2   Whipple防护屏超高速碰撞数值模拟研究

柳森等^[21]进行了一系列有关Whipple防护屏的超高速碰撞实验，本文中选取正碰撞04-0080，04-0090两组典型过程^[21]对其进行数值模拟研究。图 5为Whipple屏结构示意图。为了减小计算量，在不影响防护屏弹孔尺寸及舱壁弹坑分布的前提下，分别将防护屏和舱壁长宽减少至50和80mm。SPH粒子直径为0.3mm，弹丸直径为5mm，共离散为2440个粒子，撞击速度为5.29和6.15km/s。防护屏和舱壁厚度均为2mm，分别离散为169344和427734个粒子。弹丸、防护屏及舱壁材料为LY12铝合金。Johnson-Cook损伤模型的具体参数^[22]见表 4。由于本算例的撞击速度处于超高速撞击范畴中相对较低的撞击速度范围，碰撞过程中所有材料处于固体状态，且Grüneisen状态方程膨胀态实际上是固相材料的线性膨胀，且冲击绝热线关系与高压固体状态方程之间联系密切，冲击绝热关系的实验数据较多^[23]。因此本算例中状态方程采用Grüneisen状态方程，LY12铝合金式中各系数为ρ=2790kg/m³，Γ=2.0，C_S=5328，S_S=1.338。

图  5  Whipple防护屏结构

Figure  5.  Construction of Whipple shield

下载: 全尺寸图片 幻灯片

表  4  LY12铝合金本构模型参数

Table  4.  Parameters of LY12 aluminium alloy of constitutive model

A/MPa	B/MPa	n	C	m	Tm/K	${\dot \varepsilon _0}$/s-1
369	684	0.73	0.083	1.7	775	0.1
CV/(J·kg-1·C-1)		D1	D2	D3	D4	D5
875		0.13	0.13	-1.5	0.011	0

下载: 导出CSV 

| 显示表格

图 6给出了撞击速度为5.29km/s时Whipple屏损伤的实验和数值模拟结果，图 7给出了撞击速度为6.15km/s时Whipple屏损伤的数值模拟结果。可以看出，舱壁的损伤分布集中在撞击中心，损伤有弹坑、斑点和侵蚀三种类型。撞击中心处的中心损伤区域弹坑密集重叠，损伤程度较重，主要由内核碎片云撞击形成，中心损伤区外围出现了由外泡碎片云撞击形成若干个环形分布的圆形撞击坑群，损伤程度较轻。将撞击速度5.29和6.15km/s的舱壁损伤情况进行对比，撞击速度值越大，弹丸破碎程度更大，弹坑及斑点数量越多，撞击中心所在的中心损伤区的损伤值越小。这表明在该速度范围内，撞击速度值越大，对舱壁的损伤破坏越小，Whipple防护屏防护效果越好，这符合Whipple防护结构的典型撞击极限曲线中粉碎段(中速段)的变化规律。由计算结果与实验结果对比可见，防护屏弹孔尺寸、舱壁中心撞击坑损伤区的尺寸和舱壁弹坑形状、分布与实验吻合较好，验证了新方法对超高速碰撞问题数值模拟的有效性，具体数值参见表 5。

图  6  撞击速度为5.29km/s时Whipple屏损伤情况

Figure  6.  Damage characteristics of the Whipple shield at impact velocity of 5.29 km/s

下载: 全尺寸图片 幻灯片

图  7  撞击速度为6.15km/s时Whipple屏损伤情况

Figure  7.  Damage characteristics of the Whipple shield at impact velocity of 6.15 km/s

下载: 全尺寸图片 幻灯片

表  5  采用SPH新方法的数值模拟结果与实验结果比较

Table  5.  Comparison of the new SPH algorithm and experimental results

编号	弹丸直径/mm	撞击角度/(°)	撞击速度/(km·s-1)	防护屏弹孔直径/mm	中心损伤区域宽度/mm	舱壁损伤情况
实验04-0090	5	0	5.29	11.5	约50	3处剥落，无穿孔
仿真04-0090	5	0	5.29	12.1	54	无剥落、穿孔
实验04-0080	5	0	6.15	12.6	约50	无剥落、穿孔
仿真04-0080	5	0	6.15	12.6	62	无剥落、穿孔

下载: 导出CSV 

| 显示表格

由于实验中存在靶板安装和撞击点散布误差引起的碰撞角度误差，防护屏弹孔直径的计算结果与实验结果的误差处在合理区间；此外，实验中并没有直接给出中心损伤区域宽度，本文中对实验结果中进行短尺测量估计，计算与测量结果在容许误差内；因为数值模拟层裂破坏存在一定困难，在计算过程中只考虑了不穿孔和穿孔两种损伤状态，实验中的剥落等层裂损伤状态的数值模拟有待进一步研究。

3.   结论

引入颗粒动力学理论(拟流体模型)建立适用于超高速碰撞的SPH新方法，将SDPH方法的应用领域拓展至超高速碰撞领域，对球形弹丸超高速碰撞薄板形成碎片云的过程和Whipple防护屏超高速碰撞问题进行了数值模拟，得到以下结论：

(1) 对球形弹丸超高速碰撞薄板形成碎片云的基础算例进行了数值模拟，得到的弹坑直径、碎片云形态及内核碎片云的形状和分布均与实验吻合较好。与使用传统SPH方法、自适应光滑粒子流体动力学(ASPH)方法的仿真结果进行对比，新方法的计算结果在内核碎片云形状和分布上更加准确；验证了新方法在求解超高速碰撞问题的可行性和准确性；

(2) 通过Whipple屏超高速碰撞的算例，模拟了不同撞击速度下防护屏弹坑尺寸及舱壁损伤情况等内容；计算结果与实验结果对比吻合较好，符合Whipple防护结构的典型撞击极限曲线，验证了新方法对球形弹丸超高速碰撞Whipple防护结构问题的有效性；为下一步研究弹丸斜撞击Whipple屏全过程以及撞击损伤的内在机理等打下了工作基础。