随着精确制导武器命中精度的提高，大威力武器可能在坑道口部发生爆炸，高强冲击波进入坑道内传播^[1]。在坑道壁面的约束下，爆炸冲击波压力衰减缓慢，正压持续时间增长^[2]，对坑道内人员和设备构成了极大的威胁。近年来，关于如何衰减坑道内的冲击波，国内外学者开展了一系列研究工作。一种方法是增加坑道的长度，但过长的坑道防护性能弱，维护投入大；另一种方法是改变坑道的构型，在主坑道的基础上增加分叉^[3]和穿廊^[4]等结构，但容易受到地理环境限制，建造成本高。研究表明^[5]，在原有坑道结构基础上增设扰流板可以加速冲击波的衰减。扰流板是指安装在坑道壁面上的一系列突出平板，平时收起贴于壁面，不影响坑道的正常使用，发生爆炸时迅速打开。

坑道爆炸实验具有成本高、可重复性差、破坏性大等缺点，因而目前主要采用数值仿真和管道冲击波绕流障碍物的模型实验进行研究。Sha等^[6]针对障碍物形状对冲击波衰减的影响开展了仿真研究，发现迎风面与冲击波来流夹角为钝角的障碍物对冲击波的衰减效果最明显。Takayama^[7]通过全息干涉法得出，阻挡面积相同时，冲击波流经球形障碍物比圆柱形障碍物衰减更快。Berger等^[8]通过冲击波绕流一个或多个障碍物的模型实验得出：只有一个障碍物时，阻挡面积是影响冲击波衰减的主要因素；存在多个障碍物时，倾角是影响冲击波衰减的主要因素。朱建等^[9]采用数值方法研究了扰流板数量增加对冲击波衰减的影响。杨科之等^[5]研究了扰流板倾角对冲击波衰减的影响。

坑道和管道内冲击波流经扰流板的作用机理基本相同，但管道内障碍物形状可以任意设置，坑道内扰流板的设置则要考虑工程实际，因此，文献[6-8]中结论的适用性受到限制。杨科之等^[5]和朱建等^[9]分别研究了扰流板倾角和数量对冲击波衰减的影响，但缺乏对扰流板几何尺寸、间距等参数影响的分析。本工作以某发射井的矩形坑道为研究对象，考虑到坑道设备的运输功能，设扰流板宽度不超过坑道宽度的1/2；同时，为保证扰流板的最大阻挡面积，将扰流板形状设置为矩形；然后基于有限元模拟软件ANSYS/LS-DYNA探究扰流板厚度、倾角、间距对冲击波衰减规律的影响。

1.   仿真模型验证

1957年，Schardin利用阴影摄像技术清晰地记录了平面激波在激波管内绕流三角形楔体的物理现象，该现象称为Schardin问题^[10]。按照Schardin的实验设计，基于ANSYS/LS-DYNA有限元软件，建立等比模型，如图 1所示。

图  1  Schardin问题的几何模型

Figure  1.  Geometric model of Schardin problem

下载: 全尺寸图片 幻灯片

几何模型由空气域和三角楔组成，空气采用MAT_NULL材料模型及线性多项式状态方程进行描述，部分参数设置见表 1，其中ρ为密度，C₁~C₆为状态方程参数，E₀为初始内能密度，V₀为相对体积。三角楔采用Johnson-Cook材料模型和Grüneisen状态方程进行描述^[11]，参数设置如表 2所示^[11]，其中G为剪切模量，E为弹性模量，ν为泊松比，A_J-C、B_J-C、C_J-C、N为本构方程参数，T_m为熔化温度。

表  1  空气的材料参数

Table  1.  Material parameters of air

ρ/(kg·m-3)	C1	C2	C3	C4	C5	C6	E0/MPa	V0
1.29	0	0	0	0.4	0.4	0	0.25	1.0

下载: 导出CSV 

| 显示表格

表  2  三角楔的材料参数

Table  2.  Material parameters of triangular wedge

ρ/(g·cm-3)	G/MPa	E/MPa	ν	AJ-C/MPa	BJ-C/MPa	CJ-C	N	Tm/K
7.83	0.30	0	0.25	496	434	0.014	0.26	1 788

下载: 导出CSV 

| 显示表格

在空气域入口施加Ma=1.34的冲击波载荷；在空气域出口施加无反射边界条件，防止冲击波发生反射重新进入空气域；空气域的其他边界面施加无位移约束，以模拟刚性管道壁面。为保证计算效率，三角楔网格尺寸取为0.05 mm，空气域网格尺寸取为2 mm，计算时间步长为0.01 ms。

图 2为不同时刻实验和仿真得到的压力流线。由图 2可知：t=0.07 ms时，入射激波先在三角楔前方发生正面碰撞，形成向上游传播的圆弧形反射波；接着，入射波继续沿三角楔表面向下游传播，绕过三角楔的上、下顶点，部分波阵面传播方向发生改变，形成衍射波；t=0.09 ms时，在尾部上、下角边缘分别形成了顺时针与逆时针运动的两个涡流，涡流不断发展变大；随着入射角的增大，下游两侧壁面发生马赫反射，两个涡流逐渐融合，马赫杆随之发生弯曲变形；t=1.25 ms时，三波点在激波管的中间轴汇合。对比发现，仿真得到的波系结构演化与实验结果吻合较好，验证了仿真模型的有效性。

图  2  实验结果与仿真结果对比

Figure  2.  Comparison of experimental results and simulation results

下载: 全尺寸图片 幻灯片

2.   扰流板结构参数对冲击波衰减影响

数值模拟过程中，设某发射井坑道为长度s=82 m、横截面宽度l=3.5 m、高度h=8.5 m的矩形。几何模型由空气域、坑道及扰流板3种材料组成，如图 3所示。空气采用MAT_NULL材料模型及线性多项式状态方程进行描述。坑道壁面为钢筋混凝土材料，可以采用MAT_BRTTLE_DAMAGE材料模型进行描述，参数见表 3，其中R_m为抗拉强度，σ_c为剪切强度，η为动力黏度，σ_s为屈服极限。扰流板采用Johnson-Cook材料模型和Grüneisen状态方程进行描述^[11]。爆炸发生在距坑道入口一定距离处，将传播到坑道入口的爆炸载荷简化为峰值超压为0.4 MPa、正压持续时间为13.2 ms的压力载荷，并沿冲击波的传播方向施加到坑道入口面的法线上。为保证计算效率，取空气域网格尺寸为40 cm，坑道网格尺寸为10 cm，扰流板网格尺寸为5 cm，计算时间步长为0.2 ms。

图  3  设置扰流板的坑道模型

Figure  3.  Tunnel model with spoilers

下载: 全尺寸图片 幻灯片

表  3  壁面的材料参数

Table  3.  Material parameters of the wall surface

ρ/(g·cm-3)	E/GPa	ν	Rm/MPa	σc/MPa	η/(Pa·s)	σs/MPa
2.70	25	0.30	3.10	14.48	0	28.96

下载: 导出CSV 

| 显示表格

建模时，空气采用欧拉单元，坑道和扰流板采用拉格朗日单元，两者通过共用节点耦合起来发生作用；忽略钢筋混凝土壁面所吸收的爆炸能量，通过约束壁面质点法线方向的运动形成刚性边界；由于坑道壁面的粗糙程度难以准确量化，将其近似为光滑壁面^[12]。

2.1   扰流板厚度影响

图 4为坑道Section Ⅰ段(见图 3)的俯视图。设扰流板宽度l_s=1.75 m，与坑道入口距离z=5 m，y方向高度与主坑道相同。为探究扰流板厚度对冲击波衰减的影响，保持扰流板宽度l_s、高度h_s不变，在10~80 cm之间每隔10 cm改变一次扰流板厚度d_s，共进行8组数值仿真。

图  4  坑道Section Ⅰ段的俯视图

Figure  4.  Top view of tunnel Section Ⅰ

下载: 全尺寸图片 幻灯片

选取d_s为10、40、80 cm 3组算例进行重点分析，图 5给出了6.5 ms时扰流板附近的压力云图，其中I_i为入射波，R_i为反射波，D_i为衍射波，E_i为膨胀波，M_i为马赫波(i取值为1、2、3，分别对应于算例a、算例b、算例c)。入射波在扰流板正面发生碰撞后，一部分反射波向坑道上游运动，另一部分则继续向下游追赶入射波。当入射波绕过扰流板的两个直角时，其传播方向发生改变，发生衍射现象。算例a中，反射波尚未追上入射波，入射波I₁与衍射波D₁相连。而算例b和算例c中，反射波追赶上入射波，此时，入射角大于临界角，发生马赫反射。入射波与反射波合成马赫波M₂、M₃，并分别与衍射波D₂、D₃相连；同时，扰流板直角拐弯处产生顺时针运动的涡流，由内部低压区和扇形膨胀波构成。扰流板厚度较小时，只在扰流板后方产生一个涡流；随着扰流板厚度的增加，扰流板附近产生一前一后两个涡流。

图  5  6.5 ms时不同厚度扰流板附近的压力云图

Figure  5.  Pressure nephogram near spoilers with different thickness at 6.5 ms

下载: 全尺寸图片 幻灯片

不同算例中坑道中轴线上波阵面的峰值超压随与坑道入口距离z的变化曲线如图 6所示。由图 6可知：当2.8 m≤z≤5.0 m时，波阵面超压均呈现增大趋势。这是由于入射波在扰流板正前方发生反射，反射波向上游运动提高了波阵面压力。当z=5.0 m时，波阵面压力开始迅速衰减，此时冲击波在扰流板直角处发生绕流，形成一系列膨胀波，膨胀波与入射波间的相互作用导致波阵面压力减小。随着扰流板厚度的增加，冲击波与扰流板的作用时间变长，算例b和算例c中衍射现象较为明显，膨胀波对入射波压力的衰减作用更加突出，但在一前一后两个涡流之间，波阵面压力会有小幅度回升。在扰流板后方，波阵面压力呈现随距离增大逐渐衰减的趋势。观察局部放大图可以得出，波阵面流经d_s=80 cm的扰流板时，与膨胀波发生两次相互作用，峰值超压衰减最为明显。

图  6  波阵面峰值超压随距离变化的曲线

Figure  6.  Curves of peak overpressure versus distance

下载: 全尺寸图片 幻灯片

在距坑道入口10 m处取A₁、B₁、C₁ 3个单元，如图 4所示。单元A₁靠近无扰流板一侧壁面，单元B₁位于坑道中轴线上，单元C₁靠近扰流板正后方的壁面。冲击波流经不同厚度的扰流板后，波阵面在3个单元处的峰值超压见表 4。由表 4可知，到达同一位置处的波阵面峰值超压随着扰流板厚度的增加呈减小趋势，即扰流板厚度越大，其对冲击波的衰减作用越明显。另外，单元A₁和单元C₁的峰值超压均高于单元B₁，这是因为冲击波在两侧壁面发生反射，波阵面压力局部升高。

表  4  不同扰流板厚度下单元的峰值超压

Table  4.  Peak overpressure of units with different spoiler thickness

ds/cm	pA1/MPa	pB1/MPa	pC1/MPa
10	0.457	0.323	0.597
20	0.451	0.320	0.588
30	0.442	0.317	0.581
40	0.434	0.315	0.577
50	0.425	0.313	0.570
60	0.421	0.310	0.564
70	0.416	0.309	0.562
80	0.413	0.307	0.558

下载: 导出CSV 

| 显示表格

2.2   扰流板倾角影响

将扰流板与冲击波来流方向的夹角α定义为扰流板的倾角。为了进一步分析扰流板倾角对冲击波衰减的影响，固定扰流板厚度d_s=40 cm，并保持扰流板与坑道壁面的垂直距离b=1.75 m不变，扰流板高度与主坑道相同，在0°~165°之间每隔15°选取1组扰流板倾角，共进行12组数值仿真。几何模型如图 7所示。

图  7  设置倾斜扰流板的坑道模型

Figure  7.  Tunnel model with inclined spoiler

下载: 全尺寸图片 幻灯片

冲击波流经不同倾角的扰流板后，计算得到波阵面在单元A₁、B₁、C₁处的峰值超压和冲量，并将其与坑道内无扰流板的工况进行对比。定义超压比ε为不同扰流板倾角的峰值超压与无扰流板工况下峰值超压的比值，冲量比η为不同扰流板倾角的冲量与无扰流板工况下冲量的比值。单元A₁、单元B₁的超压比和冲量比随倾角的变化曲线如图 8所示，单元C₁的超压比和冲量比随倾角的变化曲线如图 9所示。

图  8  单元A₁和单元B₁处超压比(a)和冲量比(b)与倾角的关系

Figure  8.  Overpressure ratio (a) and impulse ratio (b) vesus the inclination angle at position A₁ and B₁

下载: 全尺寸图片 幻灯片

图  9  单元C₁处超压比(a)和冲量比(b)与倾角的关系

Figure  9.  Overpressure ratio (a) and impulse ratio (b) vesus the inclination angle at position C₁

下载: 全尺寸图片 幻灯片

由图 8可知，单元A₁和单元B₁处超压比和冲量比均小于1，其原因是冲击波与扰流板相互作用的过程中，波阵面破碎重组，消耗部分能量，加上膨胀波对冲击波的衰减作用，使单元A₁和单元B₁处的超压和冲量小于无扰流板的工况。由图 9可知，单元C₁处超压比和冲量比均大于1，其原因是在扰流板后方发生马赫反射，形成局部高压区，使单元C₁处的峰值超压和冲量大大增强。另外，从图 8和图 9还可以看出，单元A₁、单元B₁、单元C₁处超压比和冲量比均随着倾角的增大上下波动，并且均满足扰流板倾角α=105°时，超压比和冲量比最小。

为进一步从波系结构角度分析扰流板倾角对冲击波衰减的影响，图 10给出了t=7.0、7.7、8.5 ms时α=60°和α=120°扰流板附近的压力云图。由图 10可知：入射波传播到倾角α=60°的扰流板时，首先在扰流板前表面发生反射，反射波阵面几乎与扰流板倾斜方向平行，并以垂直于扰流板的速度v向上游运动；随后，入射波在拐角处发生衍射，波阵面弯曲变形，形成顺时针的涡流；当t=7.7 ms时，入射角大于临界角，在坑道上壁面发生马赫反射，入射波I₁与反射波R₁合成马赫波M₁，马赫波另一端与衍射波D₁相连。当入射波流经α=120°扰流板时，波的反射和衍射几乎同时发生，反射波R₂、衍射波D₂均与入射波I₂相连，此时，大部分反射波在板前表面和下壁面之间，而衍射产生的涡流则沿着板的后表面逐渐扩展，扰流板的前、后表面分别形成反射高压区和膨胀低压区；在此期间，反射波与衍射波能量不断消耗；最后，当t=8.5 ms时，反射波与衍射波脱离入射波，而膨胀波E₂紧跟在入射波I₂后，使入射波强度逐渐衰减。

图  10  α=60°和α=120°扰流板附近的压力云图

Figure  10.  Pressure nephogram of spoilers with α=60° and α=120°

下载: 全尺寸图片 幻灯片

不同扰流板倾角情况下，冲击波与扰流板、壁面之间的相互作用极其复杂。无法量化冲击波的反射方位、衍射角度及叠加过程，也难以估计涡流的运动方向与膨胀程度，并且波阵面之间还会相互影响、彼此制约，因此冲击波超压、冲量与扰流板倾角之间没有呈现出特定的、相对简单的规律。

2.3   扰流板间距影响

为分析扰流板间距对坑道内冲击波衰减情况的影响，本节以装有6个扰流板的坑道为例进行仿真计算。设扰流板在坑道一侧均匀排列，相邻扰流板的间距为i，宽度和高度保持不变，厚度d_s=40 cm，倾角α=105°，通过计算得到了坑道内不同时刻的压力流场，如图 11所示。可以看出，冲击波遇到扰流板的阻挡后，在扰流板正面发生第一次反射，反射波折转向坑道入口运动。同时，另一部分入射波试图绕流，波阵面弯曲变形，产生圆弧状的衍射波。由于扰流板是非流线型，在边缘处会形成很薄的边界层，容易发生分离脱落，形成大小各异的旋涡。前驱波阵面独自向前传播，与下一个扰流板碰撞后，又会发生波的反射、绕流和衍射等一系列现象，继续重复上述过程。值得注意的是，在运动过程中反射波仍然会在壁面与地面之间反射叠加，不断提高波阵面的强度，相邻扰流板之间生成了局部高压区。为了降低局部高压对坑道壁面和扰流板产生的破坏作用，工程上可以采用选取高强度材料、增加坑道壁面厚度或者在壁面和扰流板表面贴上吸能材料等防护措施。

图  11  不同时刻坑道内压力流场

Figure  11.  Pressure flow field in tunnel at different times

下载: 全尺寸图片 幻灯片

为获得最有利于冲击波衰减的扰流板间距，保持扰流板数量和形状不变，在1~10 m之间每1 m选取1组扰流板间距i，共进行10组数值仿真。同时，在扰流板后方，从距坑道出口20 m开始每隔2 m取1个单元，共10个单元，依次记为单元A~单元J，如图 11(a)所示。通过仿真计算可得，各单元峰值超压随扰流板间距的变化曲线如图 12所示。

图  12  单元超压随扰流板间距和距离变化的三维曲线

Figure  12.  Three dimensional curve of over pressure changing with interval and distance

下载: 全尺寸图片 幻灯片

由图 12可知，单元A~单元J的超压均随着扰流板间距的增大呈现先下降后上升的趋势。当扰流板间距i=6 m时，超压曲线达到最小值，再继续增大扰流板间距，超压反而迅速升高，表明扰流板间距i=6 m时最有利于冲击波的衰减。分析出现上述现象的原因：一方面，扰流板前后间距过小，冲击波在绕流时，涡流发展不充分，膨胀波与入射波的相互作用时间较短，对波阵面强度的削弱效果不明显；另一方面，扰流板前后间距过大，膨胀波与入射波脱离后，入射波继续向前运动，在坑道壁面和地面之间来回反射，使得波阵面超压在减小之后又有所增大，冲击波的衰减过程相对缓慢。将扰流板间距i=6 m时各个单元的峰值超压连接成一条曲线，即图 12中的超压衰减曲线。冲击波经过扰流板后峰值超压不断减小，随着传播距离的增加，衰减趋势逐渐变缓。

3.   结论

针对爆炸冲击波在装有扰流板坑道内的传播进行了仿真研究，在扰流板宽度一定的条件下，探究了扰流板厚度、倾角和间距对冲击波衰减规律的影响，得到以下结论：

(1) 扰流板厚度d_s越大，绕流产生的膨胀波对冲击波的衰减作用越明显；

(2) 冲击波与扰流板、壁面之间的相互作用极其复杂，扰流板倾角与超压、冲量之间没有特定的规律；当扰流板厚度d_s=40 cm，扰流板倾角105°时，冲击波的超压比和冲量比最小；

(3) 当扰流板厚度为40 cm，倾角为105°时，冲击波超压随着扰流板间距增大呈现先降后升的趋势，扰流板间距6 m超压达到最小值。