高压气体载荷下预制破片与空气冲击波的运动关系

夏晓旭，宁建国，李 健

（北京理工大学爆炸科学与技术国家重点实验室，北京 100081）

在高效毁伤领域，杀爆战斗部一直备受关注。炸药爆炸后产生的破片和冲击波是两种最主要的毁伤元。在早期的毁伤理论研究中，一般把高速破片的穿甲作用和冲击波冲量效应进行单独考虑，再组合得到综合毁伤结果。近年来，人们意识到在近距空爆情形下，破片和冲击波对结构的破坏效果存在耦合效应，且大于二者单独作用时的效果之和。对此，学者们开展了大量研究工作。Marchand 等[1]通过试验研究了平板目标在破片、冲击波单独作用及联合作用下的响应。Nyström 等[2]对冲击波和破片的耦合毁伤进行了数值模拟研究，指出耦合效应对目标的破坏效果大于两种载荷单独作用时产生的破坏之和。Kong 等[3]、Li 等[4]、曹兵等[5]、张志倩等[6]对不同目标的耦合毁伤进行了研究。陈长海等[7-8]探讨了破片与冲击波联合毁伤对目标作用的时序问题，通过理论推导提出了耦合作用区间。龚超安等[9]、王庆[10]通过理论研究了临界爆距问题。陈兴等[11]通过实验和数值模拟研究了临界爆距问题。

破片与冲击波对目标不同时序的耦合毁伤与二者的动力学过程有关。目前对这方面的研究较少，更多的是通过半经验公式去研究整个过程，而忽略掉一些实际情况。李茂等[12]、郑红伟等[13-14]通过数值模拟分析了冲击波在传播过程中与破片发生的绕流和反射现象。由于爆轰波结构复杂，给绕流过程、波系结构以及破片运动到临界爆距处与冲击波之间的追赶问题的分析带来困难。因此采用简化模型，通过能量守恒，将TNT 炸药等效成高温高压气体，将其作为载荷加载到预制破片上，通过数值模拟的方法，对这一过程进行研究，旨在更加清晰地分析冲击波与破片之间的绕流作用。通过改变预制破片的数量、与高压气体间距，研究了不同工况下临界爆距问题。

要研究冲击波和破片的联合毁伤作用，必须考虑冲击波与破片的相互运动关系，本质上就是要考虑冲击波和爆轰产物对破片的驱动作用。真实的凝聚相炸药驱动预制破片问题在物理上是一个非常复杂的问题，直接数值模拟和实验虽然可以从整体上相对准确地预测预制破片速度，但是很难精细地描述冲击波和破片在运动初期的相互作用细节。基于此，本研究将复杂的问题简化为高温高压气体与圆形刚体破片的相互作用问题，以清楚地揭示决定破片-冲击波运动关系的物理机制，为两者的联合毁伤提供理论基础。

1 物理模型

1.1 控制方程

考虑无限长的圆柱形高温高压气体和刚体预制破片，忽略轴向差异。在此条件下，问题可以简化为二维模型，如图1 所示。

图1 模型示意图Fig. 1 Schematic diagram of the model

相应的控制方程为[15]

守恒分量U和通量F和G的定义为

1.2 状态方程

本研究中，将凝聚相炸药等效为高温高压气体，采用理想气体状态方程

式中：R为理想气体常数；T为温度；e为单位质量内能； γ为多方指数，本研究取1.4。

1.3 数值模拟方法

通量求解采用Roe-HLL 黎曼求解器[16]，可精确捕捉以接触间断为代表的线性波，混合求解鲁棒性更强。当通量的雅克比矩阵特征值过小时，会违反熵条件并产生非物理解，因此需要引入熵修正方法。本研究采用Haeten-Yee 型熵修正[17]，以解决该黎曼求解器计算出现的Carbuncle 不稳定现象[18]等非物理解。采用有限体积法进行离散，空间离散采用MUSCL-Hancock 方法[19]，该迎风格式重构具有二阶精度，并采用Minmod 耗散型限制器[20-21]，限制重构时的斜率。时间离散采用二阶Runge-Kutta 格式。计算时通过设置CFL 数小于0.9 来控制自适应时间步长。

1.4 圆形破片界面处理

本研究将圆形破片视为刚体，其边界由Level-set 函数确定，即边界为函数值为零的等值曲线。在数值模拟中，气体流动状态通过求解欧拉方程获得，同时也可以得到圆形破片边界上的压力。圆形破片边界曲线被分割成若干直线段，对曲线上的压力进行线积分可以得到破片整体的受力和加速度，在一个时间步长上进行积分，可以得到破片的速度，二次积分可以得到破片的位移。在下个时间步里，重新初始化Level-set 函数，完成对破片运动的求解。

2 理论分析

假定炸药爆炸后，炸药能量E0转化为破片动能Es、产物动能Eg和内能Ee，即

2.1 破片初速

求解式(5)可得

2.2 等效TNT 当量

2.3 空气冲击波传播规律

空气冲击波波阵面传播时间ts与传播距离Rs可由下式确定[7-8]

2.4 破片运动规律

破片飞行时间tk与飞行距离Rk可由下式确定[24]

式中：mk为破片质量，u0为破片初速。

本研究中，QTNT取4 200 J/g，炸药密度取1.6 g/cm3，给定炸药半径，可得炸药质量，由E0=mQTNT可得炸药所含总能量E0。本研究的物理模型忽略黏性效应，不考虑能量的耗散，则炸药的总能量可由式(4)具体表示，初始时刻的破片动能Es、爆轰产物动能Eg均为零，采用瞬时爆轰假定，炸药瞬间反应完毕，总能量E0全部转化为初始时刻爆轰产物的内能Ee，通过式(3)可以求得初始时刻爆轰产物压强。采用无反应Euler 物理模型，在数值模拟中不考虑化学反应和爆轰波的产生与传播，从能量守恒的角度，用高温高压气体等效TNT 炸药，在数值模拟设置中，初始时刻高温高压气体密度和炸药密度相等，其能量与炸药总能量相等，即与初始时刻爆轰产物的内能Ee相等，初始压强和初始时刻爆轰产物压强相等。

3 数值模拟结果与讨论

3.1 数值模拟设置

高温高压气团半径为R，密度为 ρ0，初始时刻速度为零。圆形刚体预制破片半径为r，密度为 ρs0，初始时刻速度为零。整体计算域为正方形，长度为l，边界均采用流出边界条件。具体参数由表1 给定。如图2 所示，当预制破片到高温高压气团的距离d为0.02 m 时，破片数量分别设置为24、16、8、4；距离d取0.04 m 时，破片数量分别设置为30、24、16、8、4。所有算例中破片均为中心对称分布。以破片数量24、相距0.02 m 为例，记为工况24-d0.02。

图2 计算域设置示意图Fig. 2 Schematic of the computational domain

表1 数值模拟初始参数Table 1 Initial parameters of numerical simulation

3.2 网格收敛性测试

以工况4-d0.02 为例，对网格收敛性进行了测试。网格分辨率分别设为1.92、3.84、7.68、15.36 pts/mm。从图3 可以看出，网格的最大分辨率为1.92 pts/mm 时，破片速度-位移（v-x）曲线最低；最大分辨率为3.84 pts/mm 时，曲线最高；最大分辨率为7.68、15.36 pts/mm 时，曲线重合，居于1.92 和3.84 pts/mm 时之间，已达到收敛，考虑到计算的经济性，本研究中网格的最大分辨率选择7.68 pts/mm。

图3 网格收敛性测试结果Fig. 3 Grid resolution test results

3.3 结果分析

图4 为工况24-d0.02 不同时刻的密度纹影图像，显示了破片与冲击波的整个运动过程。初始时刻，高温高压气体与周围空气存在一个强间断，这本质上是一个黎曼问题。高压气体迅速膨胀，压缩空气产生向外传播的冲击波，同时产生向内的稀疏波。在运动初期，冲击波强度高，波后的流场质点速度也很大，而破片的启动和加速需要时间，因此初期速度较小。冲击波穿过圆形破片之间的间隙发生绕射作用，冲击波后的超音速气流在圆柱后形成花瓣状的激波波系结构，同时也可以看到，冲击波绕射后相互碰撞导致局部扭曲的波阵面及其波后的流动不稳定性，如图4(b)所示。圆形破片在前后压力差作用下向前加速，冲击波及流场质点速度随着传播距离的增大而迅速减小，破片与空气冲击波波阵面的相对距离减小，如图4(c)所示。随着破片向外运动，与花瓣状的波系结构间的相对位置发生变化。花瓣状的波系结构是由激波之间的相互碰撞产生的，花瓣所在区域的介质质点速度较大，而在花瓣所在区域的前方，受到花瓣波系的作用较小，介质速度相对较低，当破片从介质质点速度较高的区域运动到较低的区域时，破片相对介质质点的运动速度大于当地声速，在其头部产生脱体激波，如图4(d)所示。随着破片与空气冲击波继续向外传播，二者的相对位置发生变化，破片追上并超过空气冲击波，如图4(e)所示。但是在很多算例中发现圆形破片并没有追上冲击波，因此两者之间的相对运动关系较复杂，取决于初始设置，特别是球的几何分布。后面将对这一问题进行详细的参数研究。

图4 工况24-d0.02 不同时刻密度纹影图Fig. 4 Schlieren diagram of density at different moments in case 24-d0.02

圆形破片与冲击波相互作用过程中的压力曲线如图5 所示。初始设置本质上是一维圆柱中心对称的黎曼问题，同时产生一道向外传播的冲击波以及一道向内传播的稀疏波，如图5(a)所示。向外传播的冲击波与圆形破片发生碰撞，发生规则反射以及其后的马赫反射。不同圆形破片上的反射波发生碰撞产生一道向内传播的冲击波，使得内部流场压力增加，破片在前后压差作用下启动并向前加速运动。外传冲击波透过圆形破片之间的通道后在凸面上发生衍射作用，形成多个前凸的冲击波，如图5(b)和图6(c)所示。冲击波绕射后相邻之间会发生碰撞，如图5(c)和图6(d)所示。初始黎曼问题形成的内传稀疏波，不断虚弱中心处的压力，随后在中心处发生反射，形成向外的稀疏波，继续削弱流场的压力，如图5(d)所示。圆形破片上形成的内传反射冲击波与向外的稀疏波发生碰撞，随后受到削弱的稀疏波继续向前传播，而同样受到削弱的冲击波继续向内传播到达中心处发生汇聚碰撞反射，并产生外传反射冲击波，如图5(e)所示。外传反射冲击波向外传播，被波阵面扫过的介质速度增大，向外迅速扩散，后方压力不断降低，见图5(f)；由于几何尺寸的稀疏作用，外传反射冲击波强度不断降低，见图5(g)。当中心区域附近压力不断降低时，空间中存在一个较大的压力梯度，此时产生一道内传冲击波并在中心处发生汇聚碰撞，之后反射回来一道外传冲击波，如图5(h)所示；最终压力场中有明显的内外两个压力波阵面，最外面为空气冲击波波阵面，里面为第2 次中心汇聚碰撞产生的向外传播的反射冲击波波阵面，如图5(i)所示。如图6 所示，运动初期，冲击波作用在破片上，发生反射和绕射，在破片之间的间隙处存在多种冲击波间相互作用（见图6(b)）；相邻的绕射冲击波之间会发生碰撞（见图6(c)～图6(d)）；最终形成向外传播的空气冲击波，其波阵面存在局部扭曲，但是整体上接近圆形（见图6(f)）。

图5 工况24-d0.02 不同时刻空间压力曲线Fig. 5 Pressure distribution at different moments in case 24-d0.02

图6 工况24-d0.02 不同时刻的密度纹影图Fig. 6 Local schlieren photography at different moments in case 24-d0.02

图7 为不同工况下破片和冲击波的v-t和x-t曲线。由于破片的阻碍作用，冲击波速度小于没有布置破片时的空爆冲击波速度，如图7(a)所示。破片和冲击波波阵面的相对位置关系可以分为3 种情况。(1)相遇两次：在近区第1 次相遇，破片追赶上冲击波，相遇距离小于20 m，相遇时间小于30 ms，随后在较远的位置第2 次相遇，冲击波追赶上破片，见图7(a)、图7(b)、图7(e)、图7(f)。(2)相遇一次：在中区破片追上冲击波，并在很短的时间内冲击波反超破片，可认为仅相遇一次，相遇距离大约在50 m 处，相遇时间约100 ms，相遇时二者速度大致相等，并接近声速，见图7(c)、图7(d)。(3)没有相遇：冲击波始终在破片之前，见图7(g)、图7(h)、图7(i)。存在一个临界速度，当破片初速大于临界速度时，冲击波和破片可以发生相遇，由于时序不同，对目标的耦合毁伤有两种情况：一是破片先对目标进行侵彻开孔，使目标结构出现弱点，提高易损性，随后冲击波再对目标进行毁伤；二是冲击波先作用于目标，使结构发生变形并产生应力集中，随后破片再对目标进行侵彻穿孔。当破片初速小于临界速度时，冲击波和破片不能相遇，对目标的耦合毁伤只能是冲击波先作用、破片后作用。

从图7(a)可以看出，冲击波和破片的传播在前期可以分为3 个阶段。(1) 冲击波快速膨胀，绕射通过破片向前传播，速度迅速降低，而破片以较低的初速向前运动，并以稳定速度衰减率减小，经过一段时间向前传播，二者速度相等，此时相距最远，冲击波在前，破片在后；(2) 冲击波速度衰减率大于破片，破片速度大于冲击波速度，破片向前追赶冲击波，二者相对距离减小，同时二者速度也在减小；(3) 破片运动在冲击波之前，但随着向前传播距离的增加，冲击波逐渐衰减为声波，可近似认为以稳定的声速向前传播[25]，而破片速度持续降低，之后小于波阵面传播速度，在较远的位置发生第2 次相遇，冲击波追赶上破片。当破片初速减小时，第1 次相遇距离增大，第2 次相遇距离减小，两个位置不断靠近，最终从相遇两次变为相遇一次。从时程曲线可以看到：初始时冲击波速度大约为2 860 m/s，由于波阵面尺寸膨胀效应，波阵面强度降低，在0～4 ms 内，冲击波速度随着时间迅速衰减；在4～9 ms 内，冲击波速度缓慢衰减；在9 ms 时，空气冲击波速度产生波动，这是由于破片追赶波阵面造成的。在图7(a)、图7(b)、图7(e)、图7(f) 4 种工况中，空气冲击波与破片相遇之后速度持续增加，之后随着传播距离的增加而减小，相遇时破片与波阵面速度相差越大，波阵面速度增加越显著。图8 为工况30-d0.04 破片追赶空气冲击波的局部压力，从图中可以清晰地看到高速破片产生了脱体冲击波，并在后方发生碰撞。破片向前追赶空气冲击波波阵面时，脱体冲击波随之与波阵面发生碰撞，波阵面强度增加，波速增大，原本光滑外凸的波阵面变得扭曲不规则，但随着传播距离的增大，空气冲击波波阵面逐渐变得光滑，速度仍然继续衰减。对于图7(b)、图7(f) 所示的2 种工况，空气冲击波与破片相遇之前速度有所增加，也是由于破片的前导脱体冲击波导致的。破片前方的脱体冲击波与空气冲击波波阵面发生碰撞，波阵面强度增加，速度增大，但很快衰减下去，波阵面强度越低，速度增加越明显。而破片速度则在很短的时间内加速到某一值，并在一段时间内缓慢地增加至最大初速，之后开始逐渐衰减。

图7 不同工况下冲击波和破片的速度及位移时程曲线Fig. 7 Time history curves of velocity and displacement of shock wave and fragment in different cases

图8 工况30-d0.04 中局部的压力云图Fig. 8 Local pressure contours in case 30-d0.04

图9 为不同工况下破片速度随距离的变化曲线。从图9 中可以看到，破片数量一定时，距离高压气体中心越远，其最大初速越小，但速度衰减快慢大致相同。这是因为距离越远，冲击波强度进一步衰减，当冲击波运动到破片所在位置时，超压相应降低，破片加速能力变弱，使得破片初速降低，高压气体的能量更多地用作空气冲击波的传播。当间距相同时，破片数量越多，初速越大，这是由破片与冲击波之间的相互作用更加剧烈所致。当冲击波绕射通过破片时，在相邻破片之间的间隙处，存在入射冲击波与壁面碰撞、反射冲击波与反射冲击波碰撞、反射冲击波与壁面碰撞等多种作用，过程较为复杂，且破片间隙越小，相互作用越剧烈，见图5(b)。破片后方会产生向后传播的反射冲击波，破片数量越多，反射冲击波之间的汇聚碰撞越剧烈，此时破片后的压力场在短时间内维持相对慢速衰减阶段，使得破片加速能力变强，初速增加，如图5(f)所示。从能量的角度讲，破片越多，通过破片间隙逃逸的气体能量越少，越多的能量用于破片的加速。

图9 不同工况下破片速度的空间分布曲线Fig. 9 Spatial distribution of fragment velocity under different working conditions

图10 为不同工况条件下冲击波与破片发生相遇的模拟结果。从图10 中可以看到，随着破片初速u0减小，相遇时间t延后，相遇距离越远，二者相遇时，冲击波速度us和破片速度uk均减小。对比数值模拟结果和理论值发现，速度比较吻合，这是因为理论速度的推导只涉及能量的计算，并考虑了气体的内能，与载荷是凝聚炸药还是高压气体无关。相遇时间与相遇距离差异较大，但趋势相同，即破片初速越低，相遇越晚，相遇距离越大。由此说明，本研究的计算结果是合理的。

图10 不同工况下相遇时n-t、n-x、n-v 曲线Fig. 10 n-t, n-x, n-v curves of encounter under different working conditions

图11 给出了统计出的所有算例中冲击波和破片是否存在相遇的情况。可以看出，初始时若破片与高温高压气体中心的间距为0.02 m，则必然存在破片追上冲击波的情况。若该间距增大，则破片不一定能够追上冲击波。另外，在本研究中，当n/(d/r)（破片数量除以间距与破片半径之比）大于4 时，破片能够追上冲击波并且超越。这一判据只考虑了初始间距和破片数量对后续运动规律的影响，并没有考虑破片的尺寸和质量以及高温高压区的尺寸。如果考虑所有影响因素，这一问题将更复杂。后续工作将对此问题进行更详细的研究。

图11 统计出的破片-冲击波的相遇情况Fig. 11 Encountering statistics of the fragment and shock wave

4 结 论

对于破片战斗部，在一些情况下，破片先行穿孔弱化目标结构，随后冲击波进一步作用，形成联合毁伤。能否形成有效耦合毁伤作用主要取决于冲击波与破片的相对运动关系，进一步则取决于初始时刻炸药的能量以及破片的质量、尺寸和几何分布等因素。采用基于网格自适应的二维流体力学Euler 程序，通过改变预制破片数量、破片与高温高压气团间距，分析了冲击波和破片在运动初期的相互作用过程和流场分布，以及破片与冲击波的追逐关系，得到以下结论。

(1) 数量一定时，破片距离高压气体中心越远，最大初速越小，但速度衰减快慢大致相同；当破片与高温高压气体中心的间距相同时，破片数量越多，反射和绕射作用越强，初速越大。

(2) 破片初速较大时，破片很快追赶上冲击波，存在二次相遇；当破片初速为562.1、552.5 m/s 时，破片与波阵面之间仅相遇一次；当破片初速进一步减小，则不能与波阵面相遇。

(3) 速度理论值和模拟结果较吻合，这是因为速度推导仅与能量有关。相遇时间和相遇距离差异较大，这与模拟初值条件与理论空气冲击波超压分布不同有关，但趋势相同，说明本研究的数值方法具有合理性。