绝热指数γ 对平面爆轰过程中不同复合波区的参数特性影响分析*

2021-11-15李晓彬

爆炸与冲击 2021年10期

张娅，李晓彬，彭帅，施锐

（1. 武汉理工大学交通学院，湖北武汉 430063；2. 中国船舶工业系统工程研究院，北京，100036）

爆炸冲击波的传播与衰减一直是爆炸领域研究的重点问题[1]。当炸药爆轰后，爆炸气体产物内部形成的稀疏波会迅速衰减爆炸高压，其中不同方向的稀疏波会在中心汇聚，形成相互作用的稀疏波复合波区，稀疏效应的叠加加快了爆炸高压的衰减速率，且穿过复合波区的稀疏波会快速追赶上空气冲击波阵面，对该阵面的压缩特性进行衰减[2]。因此，在很多冲击波问题的研究中，为在实验段得到稳定的冲击超压，通常会在激波管实验中通过加长爆轰驱动段[3]，来减缓反射稀疏波到达激波阵面的时间。可以看出，该稀疏复合波区的特性直接影响了爆炸过程的衰减特征。除稀疏复合波区外，随着流场运动，爆炸区域还会形成其他类型的相交复合波区，且每个波区的传播和衰减特征均不同。

在平面爆轰流场中，该类复合波区均可以简化为两相向扰动波相交作用的特性耦合。在应用特征线法的爆炸理论中，已提出了两相向扰动波交汇后流场的一般解形式[4]，但由于相交特性复杂，并没有通用的解析解，仅在假定爆炸气体绝热指数γ=3 时有确定的特性，即相交后两簇波在x-t平面上的特征线仍均为直线。张守中[5]基于该γ=3 时的特殊解形式，对多个典型爆轰产物流场中稀疏复合波区的数值解进行了计算，并将计算解代入到后续的流场分析。可以看出，对于复合波区而言，该特殊解形式简单，方便用于流场参数的快速估算。因此，当对计算精度要求不高时，对高密度凝聚炸药常近似取γ=3.0。

但是在实际过程中，爆炸气体内部这种从极稠密的高压状态变化到较稀疏的低压状态，过程是极其复杂的，尤其是过程中爆炸气体的绝热指数γ 是持续变化的。常用高密度凝聚炸药的初始γ0测定值一般在2.2～3.50 之间[6-8]，在根据γ 律等熵状态方程进行膨胀衰减计算时，常采用在高压γ0到低压（1.1～1.4）分段取值计算的方法[9-10]，所以γ≠3 是这些复合波区的常态。

为研究γ=3 和γ≠3 不同条件下爆炸不同复合波区的特性差异，本文中基于特征线法和理想气体的γ律等熵状态方程，对一平面爆轰过程中不同复合波区的波系相交特性进行理论规律分析。由于在一般解形式中，γ≠3 的复合波区内相交的特征线为迹线未定的曲线，难以推导理论的位移-时间公式（即x-t解），本文中利用MATLAB 对其平面爆轰过程进行流场模拟，验证并分析不同复合波区流场内的参数变化特性。

1 平面爆轰过程中复合波区的形成

在如图1 所示的截面积为1 m2的圆形无限长管内，T0时刻在x(-1,1)区间内设有TNT 固体装药。基于瞬时爆轰模型，在炸药原体积空间内，爆炸气体产物的密度取为固体装药的密度，初始压力取pCJ爆压乘以等效系数0.434 8[11]，初始质点速度为零。

图1 平面爆轰过程中复合波区的形成及相对位置（+右行，-左行）Fig. 1 Formation and relative positions of complex wave zones in plane detonation(+ showing a going-right wave, - showing a going-left wave)

利用特征波系的运动分析该爆炸流场，当t＞0 时，阶段1：空气受压缩形成冲击波阵面，在内部中心稀疏波的作用下爆炸气体压力迅速衰减；阶段2：原爆炸气体产物（③区）全部被稀疏扰动，两相向中心稀疏波相遇并相交作用，形成稀疏波-稀疏波复合波区（⑤区）；阶段3：稀疏波波头完全扰动④区，继续向前扰动②区（扰动后形成稀疏⑥区），在到达与空气的交界面后，由于爆炸气体的波阻抗大于空气的波阻抗，会在空气中形成透射稀疏波，并向爆炸气体中反射一个压缩波，该压缩波又与原稀疏波相交作用，形成稀疏波-压缩波复合波区（⑦区）。

为分析阶段2、阶段3 中生成的两复合波区的参数特性，首先利用特征线法确定流场的初始波动特征（阶段1）。为便于分析流场中的能量流动，计算中假定高压气体和空气均为无黏理想气体，满足理想气体状态方程和γ 律等熵方程。绝热指数γ 不同，整个流场的参数分布不同。

对于TNT 炸药的pCJ爆压测定值一般在19～21 GPa，以pCJ=21 GPa 时为例，爆热e为4 187 kJ/kg，则γ=2.34。根据类激波管问题的特征线理论，可以得到图1 中阶段1 时段无限长管内各区域的状态参量，见表1，且①区前端的冲击波阵面速度vs=5 266.14 m/s。当左右两稀疏波在原点相遇时，阶段1 结束，记此时刻为T1，此时无限长管内各区域的状态参量分布，如图2 所示，图2(a)、(b)、(d)为对称半流场，图2(d)中Ed为线能量密度，Ed-int、Ed-kin分别代表内能密度、动能密度，沿管长度方向积分可得该时刻系统的内能及动能。

表1 阶段1 各区域的状态参量Table 1 Characteristic parameters of each region in stage 1

图2 两稀疏波在原点相遇时流场内的特征参数分布Fig. 2 Distribution of characteristic parameters when two rarefaction waves meet at the origin

2 当γ=3 与γ≠3 时复合波区相交特性的差异

由特征线理论可知，两个相向行波相交满足黎曼不变特性，即右行波的每一条C+特征线，沿该线的黎曼不变量J+为定值；左行波的每一条C-特征线上，沿该线的黎曼不变量J-为定值。对于理想气体，其J±=u±2c/(γ-1)。由此可得任一条C+和C-特征线相交后，交点位置M的u、c参数满足：

由于中心稀疏波特有的中心发散特性，其运动（④区）有通解方程。此时若一条C+或C-特征线穿过常态区与一相向中心稀疏波相遇，其迹线的解析方程为：

式中：xc为中心稀疏波发散中心，x0为相向特征线出发点，ccons为常态区气体声速。

图3 为图1 所示流场的x-t平面示意图，图中⑤区就是两稀疏波特征线的相交区域，在交汇前各单波系的特征线运动均是直线。根据公式（2）可以得出⑤区边界线即图3 中左右波头特征线C+0和C-0的x-t迹线。当γ=3 时，式（2）变化为一直线形式。当γ≠3 时，该x-t解为一明显的多项式曲线形式，可以看出当γ＞3 时C+0或C-0的波速是不断减小的，曲线向内凹；γ＜3 时C+0或C-0的波速是不断增大的，曲线向外凹。

图3 两中心稀疏波相交的x-t 平面Fig. 3 x-t plane of intersection of two central rarefaction waves

而对于复合波区内部任意两条相交特征线的x-t方程，其一般解形式[2]为：

式中：F+、F-分别为与特征线波速u+c、u-c相关的系数。由式（1）可以看出，特征线相交会影响其运动波速，在复合波区内一条特征线会受N条相向特征线的影响，其x-t迹线变为斜率不断变化的曲线，所以式(3)很难得到直接的解析解。但由于在γ=3 条件下，dx/dt=u±c=J+(或J-)=常数，所有特征线均沿直线向前推进，所以可以通过左右行波的初始条件确定F+、F-。

如图3 中两中心稀疏波相遇时，由右行稀疏波t=0 时，x=-1；左行稀疏波t=0 时，x=1；可得F+=-1，F-=1。可解得复合波⑤区任一交点的u、c参数：

对于理想气体而言，该点的压力、密度均可通过该气体声速得到（其中p3、ρ3、c3为③区参数）：

所以在γ=3 时复合波⑤区内各处气体的声速c、压力p、密度ρ 都是一个均布值，与气体分布位置无关，仅随运动时间t变化，并且随着t增长，c、p、ρ 逐渐减小，而气体的质点速度u同时受位置和时间的影响。可以看出该特殊解形式简单，且方程可解析。但当γ≠3 时，复合波区内的各特征线均为变化曲线，任一交点的u、c参数无通解。图4 显示了近似利用γ=3 的c值均布特性求解后原图2 中γ=2.34 时流场的特性分布变化，其中E为能量，Eint为内能，Ekin为动能。

图4 应用c 值均布特性后系统流场特性变化Fig. 4 Changes of flow field characteristics after the application of c-value uniform distribution

从图4（a）可以看出，复合波⑤区内除气体质点速度外，内部的其他参数都是均布状态。图4（b）计算了两稀疏波从初始T0、相遇T1到④区完全被扰动为复合波区的T2时刻过程中各区域的能量变化，可以发现，系统内的总能量在复合波⑤区形成后逐渐减少，到T2时刻总能量已减少了37.53%。这说明，当γ≠3 时其复合波区的特性应与γ=3 时有较大差异。

从波系相交特性的角度分析，如图5 所示，当γ≠3 时选取3 条右行C+特征线与左行中心稀疏波相交，C+0线按已有解析方程式（2）画出。由于图1 中③区为常态区，根据特征线理论，该左行波C-族线的右向黎曼不变量J+与C+0线的J+0相同，所以在C+0线通过后，各C-线的u、c不变，仍按原斜率u-c运动，此时左行波仍保持中心发散特性。逐一相交迭代计算。同理可推，当γ≠3 时在图3 中相向运动的两个中心稀疏波相交后，C+族线和C-族线都不再具有中心特性，变为非中心稀疏波，即其波系发散特性发生了改变，因此难以确定其各特征线相交的具体位置。

图5 中心稀疏波的发散特性Fig. 5 Divergence characteristics of central rarefaction waves

由上可知在x-t平面上无法直接对比其特性差异，改从能量角度分析，理想气体的内能与压力p相关，动能与密度ρ、质点速度u相关，结合式（5）可知，该系统的能量特性主要受u、c变化的影响，所以为分析其特性差异，可以从复合波区的u-c特征进行切入。

3 当γ=3 时复合波区的u-c 平面特性

由特征线理论可知，对于每一个不变量J+/J-，在u-c平面上必有一条曲线 Γ+/Γ-与之对应，且曲线Γ+和 Γ-的交点处的(u,c)值与C+和C-特征线相应交点的（x,t）值具有一一对应关系。由于理想气体的等熵特性，u-c平面上的 Γ+/Γ-曲线均为直线。

（1）复合波⑤区的u-c平面特性

图6 为两同等强度中心稀疏波相交的x-t平面和u-c平面。假定tk时刻，x-t平面上左行波波头C-0与C+k线相交(A点)，由对称可知，此时右行波波头C+0与C-k线相交(B点)，线AB为复合波区的长度，对应到u-c平面，可以找到相应的复合波区边界点A、B，当γ=3 时，x-t平面上特征线为直线，复合波区中点位置为P点，u=0。根据式（4）已知c与分布位置无关，为均值，因此在u-c平面，AB之间的路径为直线A-P-B。

图6 两同等强度稀疏波相交的平面特性Fig. 6 Plane characteristics of the intersection of two rarefaction waves with equal intensity

（2）复合波⑦区的u-c平面特性

图1 中复合波⑦区的类型是稀疏波和压缩波的相互作用。为分析⑦区的特性，首先了解两同等强度的压缩波相遇的复合波特性。在图7（a）中分别在长管内-x0、x0处设置活塞，管内气体静止，初始声速c0。利用活塞相向匀加速运动，加速度为a，压缩管内气体，生成两相向等强度压缩波。两压缩波运动到中点（x=0）处相遇，形成压缩波-压缩波的复合波区。

从式（11）可以看出cM/cB≠1 ，因此在该压缩波复合波区的u-c平面上AB路径并不是一条直线。将a、c0、tk赋值后发现， (t-M+t+M) 在中心P点最大，在两端A、B处最小，因此cP＞cB。由于在u-c平面上压缩波和稀疏波波头波尾走向相反，因此该波区的u-c平面如图7（b）所示，其AB路径为曲线A-P-B。

图7 两同等强度压缩波相交的平面特性Fig. 7 Plane characteristics of the intersection of two compression waves with equal intensity

将该分析方法延伸至图1 复合波⑦区的类型分析中，在长管内利用图8（a）中活塞向左匀加速运动生成不同强度（ -a1～-a5）的压缩波与迎面稀疏波相交，两相向波在x-t平面上的初始位置分别为( -x0,0)、（x0,0）。当γ=3 时结合式（4）和式（6）的分析，可以得到图8 波区交点的F+、F-解为：

将F-i展开，由此推出该波区任意时刻的声速解c：

由式（13）可以看出，在指定a和时刻t后，该波区范围内c仅与各压缩线的初始形成时间t-i有关，因此t时刻由压缩波波头到压缩波波尾c逐渐增大。但对于单条压缩线来讲，在复合波区内其c值随t的增长逐渐减小，这就是稀疏波对压缩波稀疏效应的体现。

对于该路径上u的特性，根据稀疏波方程x=(u+c)t+F+可得：

由于在t时刻，(x+x0)/t项随x增大而增大，-c值会随x增大而减小，而c值变化速率主要受a的影响，因此该路径上u的特性也主要受活塞加速度a的影响，即受压缩波强度的影响。假定tk时刻，该复合波区边界交点为A、B，改变活塞初始加速度（a1～a5），以式（13）～（14）对应计算图8（b）u-c平面上的AB路径。计算得出，当a较小（a1～a3）时，-c项的影响较小，AB路径上u基本呈单调递增的状态，但当a继续增大（a4～a5）时，-c项的影响开始增大，AB路径上u开始有明显的先减小后增大的趋势，若继续增大a，可预见递减会成为主要趋势。

图8 不同强度压缩波和稀疏波相交的平面特性Fig. 8 Plane characteristics of the intersection of a compression wave and a rarefaction wave with different intensities

综上分析可以看出，当γ=3 时各类复合波区的c值均只跟时间有关，但在不同类型复合波区中分布形式不同，而u值与x和t均相关，在同类型波相交的复合波区内是单调递增的，若相交波为等强度波，对称面u=0。在不同类型波相交中u值特性与参与波的强度直接相关。u、c特性的组合就是复合波区的主要特性。

4 当γ≠3 时复合波区的u-c 平面特性

虽然γ≠3 时依据式（3）流场难以求解，但对于具有对称特性的复合波⑤区，仍可近似借鉴γ=3 时的直线特性，假定其特征线也是沿直线运动，同样得到式（4）的c值均布特性（当然此假定过程并不等熵），然后放在原γ≠3 的x-t和u-c平面上进行对比和特征规律分析，如图9 所示。

图9 γ≠3 时两同等强度稀疏波相交的平面特性示意图Fig. 9 Plane characteristics of intersection of two rarefaction waves with equal intensity at γ≠3

在图9 中假定tk时刻，A0-P0-B0为假定特征线沿直线运动时的c值均布直线。而其实际情况由式（2）及图3 可知，当γ≠3 时，x-t平面上复合波区内的特征线均为曲线，此时需分两种情况：

（1）当γ＞3 时，波头C+0/C-0迹线的波速是不断减小的，所以波头C+0/C-0会比直线运动时更晚与C-k/C+k线相交（tk1＞tk），x-t平面上线A1B1为复合波区的长度，中点位置为P1点，u=0。但在u-c平面上复合波区边界点仍为A0、B0两点，而P1点在x-t平面上高于P0点，也就是在u-c平面上其相交特征线的编号高于P0点，因此当γ＞3 时，在u-c平面上A1B1之间的路径为曲线A0-P1-B0，即此时复合波⑤区的特性参数变化规律为：c（或p、ρ）在复合波区中心处最小，向着复合波区边界逐渐增大。u为矢量，在对称中心处均为0，数值上向着复合波区边界逐渐增大。

（2）当γ＜3 时，波头C+0/C-0迹线的波速是不断增大的，所以波头C+0/C-0会比直线运动时更早与C-k/C+k线相交（tk2＜tk），x-t平面上线A2B2为复合波区的长度，中点位置为P2点，u=0。但在u-c平面上复合波区边界点仍为A0、B0两点，而P2点在x-t平面上低于P0点，也就是在u-c平面上其相交特征线的编号低于P0点，因此当γ＜3 时，在u-c平面上A2B2之间的路径为曲线A0-P2-B0，即此时复合波⑤区的特性参数变化规律为：c（或p、ρ）在复合波区中心处最大，向着复合波区边界逐渐减小。u为矢量，在对称中心处均为0，数值上向着复合波区边界逐渐增大。

从上述分析可以看出，当γ＞3 和γ＜3 时复合波⑤区的u-c路径特性是相反的。

为验证上述结论，也为得到非对称复合波⑦区在γ≠3 时的u-c平面特性，本文利用MATLAB 基于特征线法对图1 所示的爆轰流场进行模拟，首先假定左右中心稀疏波④区各有N条特征线从中心发出，并在每个区域预先设置特征空矩阵，以便于下一区域直接调用上一区域的边界作为初始条件。为保证计算精度，计算中采用尽可能多的特征线充满流场，取N=10 000。

图1 中阶段3 有2 种不同的复合波形式，对于稀疏波相交的复合波⑤区，以公式（1）为迭代准则，假定在每个迭代胞格内特征线为直线，计算胞格各交点的（x,t）及（u,c），再推进到下一个胞格，如图10(a)所示；对于在介质分界面处的复合波⑦区，其入射波、反射波和透射波的关系[12]按下式计算：

图10 对复合波区的迭代思路Fig. 10 The iterative method for the complex wave zone

为充分讨论γ≠3 的各种情况，本文中通过调节初始爆热，得到7 组不同γ 值：3.5、3.2、2.7、2.34、2.0、1.7、1.4。由于该流场程序具有通用性，因此增加γ=3 工况作为对比。为验证该流场程序的合理性，以γ=2.34 时为例，MATLAB 计算生成的特征线流场如图11 所示。每个特征线上节点位置处都包含了6 个特征参数：压力、密度、质点速度、气体声速、内能密度、动能密度。选定时刻，即可确定该时刻内流场各处的参数分布。计算过程在Tn=5.93 ms 时刻报错，显示⑦区N7=173 与N7=172 反射压缩特征线在26.569 5 m处相交。

图11 MATLAB 计算生成的特征线流场Fig. 11 Characteristic line flow field generated by MATLAB calculation

首先选定3 个特征时刻进行流场对比：T2时刻（右行稀疏波波头刚穿出复合波⑤区）、T3时刻（右行稀疏波波头刚穿出②区，即到达介质分界面）、Tn时刻（计算结束时刻）。流场参数分布如图12 所示，由于特征线流场左右对称，只显示右半流场。

由图12（d）各区域随时间的能量变化曲线可以看出，Tn时刻之前流场的总能量保持守恒状态，说明该MATLAB 迭代流场是合理可行的，若忽略Tn时刻报错继续计算，流场总能量开始增大。

图12 T2、T3、Tn 时刻的流场参数分布Fig. 12 Distribution of flow field parameters at times T2, T3 and Tn

（1）复合波⑤区的u-c平面特性

从图11 的x-t平面可以看出，中心稀疏波相交作用的复合波⑤区在T1时刻开始形成。图13（a）为该区与图11（x,t）交点对应的u-c平面，从图中不同时刻的u-c路径及图12（a）～（c）中不同时刻⑤区内p、ρ、c的分布规律，可以看出在该波区内c（或p、ρ）均在波区中点处最大，向着两侧波区边界逐渐减小。这与上述γ＜3 时复合波⑤区的u-c特性分析是一致的。图13（b）为γ 不同时各流场中复合波⑤区在各自T2时刻的u-c路径，从对比中可以看出：不同流场的u值特性是一致的，均是沿着波区长度方向逐渐增大；而其c值特性则随着γ 增大，其沿波区的路径逐渐由上凸曲线（γ＜3，中心c值最大）-直线（γ=3，c值均布）-下凹曲线（γ＞3，中心c值最小）变化，且离γ=3 值越远时，曲线特征更明显。当γ 值越接近3 时，其曲线路径越接近直线特性。

图13 爆轰流场中复合波⑤区的u-c 特性Fig. 13 u-c characteristics of complex wave zone ⑤ in the post-detonation flow field

由于在爆轰流场中，复合波⑤区范围最广，且前期衰减速率及衰减跨度大，如对比γ=2.34 时流场中图2 和图12，波区中点p值从T1时刻9 130.8 MPa 快速衰减到了T2时刻的44.07 MPa，该时段内中心ρ 值从1 630 kg/m3衰减到了166.8 kg/m3，因此该c值特性不同影响的p、ρ 参数对系统的内能、动能均产生了较大影响。

（2）复合波⑦区的u-c平面特性

对于在T3时刻开始形成的复合波⑦区，通过图11 中Tn时刻压缩波特征线交点的位置可以看出，相交位置位于靠近介质分界面，也就是靠近反射压缩波波尾处。若将介质分界面DL看成一个可透射的活塞，若活塞向左做均匀加速、且波前静止的情况，间断出现的位置一般在波头上，而该⑦区流场中间断率先发生在了波尾处，说明后生成的反射压缩线追赶速率不断增大，也就是可透射的活塞DL向左加速度是不断增大的，间断将由波后向波前推进。在其他7 组不同γ 的流场中均出现了类似的间断位置，而且随着γ 增大间断位置出现的越早。因此可以得出，该平面爆轰过程中的复合波⑦区是由右行稀疏波（其中γ≠3 时为非中心稀疏波，γ=3 时为中心稀疏波）和左行变加速反射压缩波相交形成的。

图14（a）为γ=2.34 时该波区与图11（x,t）交点对应的u-c平面，从c值特性看，在不同时刻该波区范围内c值均是沿压缩波波头向波尾（介质分界面）逐渐增大的，这与γ=3 时的规律特性是一致的。从u值特性看，在T7时刻前u值沿着波区主要呈递增的趋势，但随着时间增长，u的增幅越来越小，且结合图12（c）中⑦区的u-x曲线，可以看出在Tn时刻该波区的u值出现了先减小后增大的变化趋势，由此也可以看出该反射压缩波的强度随生成时间是不断增大的。不过从整体上该波区内各位置的u、c还是随时间不断衰减的。

图14（b）为γ 不同时各流场中复合波⑦区在各自Tn时刻的u-c路径，从对比中可以看出：不同流场的c值特性是一致的，均是沿着波区长度方向逐渐增大，而其u值特性随γ 增大逐渐由单调递增变化为先减小后增大，该先减小的特性也间接表现了右行稀疏波的强度，说明随γ 增大，流场生成的稀疏波强度越大。

图14 爆轰流场中复合波⑦区的u-c 特性Fig. 14 u-c characteristics of complex wave zone ⑦ in the post-detonation flow field

可以看出，在该平面爆轰过程中不同区域都有其特定的衰减特性，根据其特性可以快速了解整个流场的参数分布情况，如复合波⑤区的边界点是流场压力最低的位置。同时若将图11 模型继续计算，在Tn时刻之后约1 ms 的时间，⑧区的右行透射稀疏波将与冲击波波阵面交汇，还会形成稀疏波和冲击波的复合阵面，此时空气超压和冲击波速开始被衰减。