水下爆炸近场峰值压力的数值模拟*

宫翔飞，刘文韬，张树道，杨基明

（1.中国科学技术大学近代力学系，安徽 合肥 230027；2.北京应用物理与计算数学研究所，北京 100094）

压力状态研究是水下爆炸研究中的重要内容。以往水下爆炸压力状态研究的距离范围主要集中于中远场，现有的中远场水下爆炸研究结果与实验符合较好。水下爆炸近场压力则具有典型的高压高频特点，压力峰值至少在GPa量级，持续时间则为ms量级，给实验测量造成很大困难[1]。近期，学者们对近场峰值压力进行了试验和数值研究[2-8]，认为在R/a≥1的一定范围内峰值压力的变化更符合指数规律，而在1≤R/a≤6范围进行拟合时则都使用了幂次形式，最后的拟合公式也不尽相同。这些试验和计算工作都证实，水下爆炸冲击波中存在相似原理，因此用数值模拟对水下爆炸近场中不同波系的影响进行一些规律性的分析是可行的。

本文中利用SPH方法进行水下爆炸近场压力的数值模拟，对模拟结果进行基于波系传播规律的定性分析，总结近场内峰值压力的变化规律，以期提高对水下爆炸近场压力变化规律的认识。

1 SPH方法

SPH方法[9-11]对物理量进行局部空间的积分近似和粒子近似，是一种拉氏无网格粒子类方法，在水下爆炸领域得到了越来越广泛的应用。本文计算中忽略黏性、热传导，使用的SPH控制方程为欧拉方程组：

式中：ρ、p、e分别为粒子密度、静水压和单位质量的内能，m为 粒子质量，v为粒子速度，下标i、j为粒子编号，vij=vi-vj，N为粒子总数目，x表示坐标轴，α ,β=1,2,3，xα或xβ分别表示不同的坐标轴方向,光滑函数Wij(R,h)有两个参数，h为光滑长度，R为i、j粒子间的距离与h的比值。

光滑函数W取为最常用的三次样条函数形式：

式中：bd在不同维数中取不同的值，一维b1=1/h，二维b2=15/(7πh2)，三维b3=3/(2πh3)。

在轴对称的柱坐标系统中，欧拉方程组的相应形式[11]为：

式中：s(r,z)为粒子在柱坐标中的位置，。

计算中没有考虑重力，为了模拟水深，在水的初始状态中增加了初始压力，比如若水深H=0，则取初始压力p0=1.01×105Pa。

2 CJ爆轰模型

在计算炸药爆轰过程时，使用了经典的CJ爆轰理论模型[12]。假定炸药微元在通过爆轰波阵面时瞬间转变为爆炸产物并完成放能，在这种假定下，反应区的宽度变为零，成为以爆速传播的强间断面。假设波阵面的传播速度为D，产物的粒子速度为u1，在自持定常爆轰状态下，产物离开波阵面的速度D-u1等于当地声速c1：

通常认为爆速D保持不变，需要由实验测定。

3 状态方程

众多学者对水下爆炸中适用的水的状态方程进行了数值模拟与试验的对比研究，认为在PH＜32.0 GPa，ρ＜2.02 g/cm3范围内[13]，多项式形式的状态方程能够适用。本文中计算的TNT药包水下爆炸算例中，水中的压力峰值也在此范围内，因此水的状态方程选用了多项式形式：

式中：ρ0为初始密度，µ为压缩度， µ= ρ/ρ0-1，EM是单位质量内能增量，A1、A2、A3、T1、T2为压强量纲拟合系数，B0、B1为无量纲系数。其具体参数如下：

爆炸产物的状态方程选取JWL形式：

式中：E为比内能，ρ0为炸药初始密度，V=ρ0/ρ，A、B、R1、R2、ω为常系数。对TNT炸药，常用值为A=373.8 GPa、B=3.747 GPa、R1=4.15、R2=0.9、ω=0.3。炸药起爆方式使用时间起爆，爆速取D=6 930 m/s。

4 自由场爆炸验证

先将自由场水下爆炸中气泡脉动周期、最大半径和中场峰值压力的计算结果与经验公式相比较，验证SPH计算中相似率的符合，结果表明了本文数值计算的可靠性。

Cole[1]在大量理论研究与试验数据的基础上，提出了自由场中水下爆炸气泡脉动第一个周期和最大半径的经验公式：

式中：T为气泡第一次脉动周期，s；Rm为气泡最大半径，m；W为炸药质量，kg；H为爆炸深度，m。对TNT炸药，一般取系数值k1=2.11、k2=3.8。为了加快计算速度，取水深H=40 m，根据经验公式，初始半径a=0.1 m(W=6.83 kg)的球形TNT炸药，气泡第一次振动的周期最大半径1.94 m。

使用二维轴对称SPH算法，水中初始压力为50.33 m高水柱在重力作用下产生的压力值，进行自由场球形药包中心起爆水下爆炸气泡的计算。得到水下爆炸气泡第一次振动的周期和最大半径分别为T=0.159 s、Rm=1.98 m。计算结果与经验公式相比的误差分别为4.6%和2.1%，是较为接近的。

Cole[1]还进行了水下爆炸中相似率的理论论证，并用实验数据作了简单而充分的证明。根据这一理论，以(R/a,t/a)作为参数，不同装药量的同种炸药具有相同的压力分布。

本文中对半径为0.1 m和0.2 m(W=54.62 kg)的球形TNT炸药进行计算，取相同比距离R/a=10处的压力曲线，a=0.2 m的模型中，炸药半径为另一模型的2倍，根据相似率其时间尺度也是另一模型的2倍。从图1中可见，以t/a作为参数时，两种计算模型的压力衰减曲线基本吻合，不同计算模型中，相同R/a处的峰值压力也是基本一致的，说明相似率符合很好。

Cole[1]总结了经验峰值压力公式，Zamyshlyaev等[14]后来对Cole的结果进行了深入和扩充，给出的峰值压力公式为：

图1 不同模型压力的相似率Fig.1 Similarity of pressure in various models

式中：a为球形炸药的初始半径，R为测点到炸药中心的距离。

本文中采用SPH二维轴对称格式，取H=0，对球形装药进行爆轰模拟。比较6≤R/a≤12内的计算结果与经验公式，令ρ=1 630 kg/m3，将经验公式中球形药包的质量与半径进行了换算，如图2所示，两者的结果最大差别不超过2%。

从对自由场水下爆炸的计算结果看，本文使用的SPH方法得到的计算结果较好地与经验公式相符合，使用SPH方法进行水下爆炸近场内的压力研究是可行的。

图2 峰值压力计算结果与经验公式的比较Fig.2 Comparison of peak pressure between calculated results and empirical formula

5 水下爆炸近场中的波系

爆炸产物内的波系对水下爆炸近场中的压力有重要的影响，不同波的传播速度和衰减规律并不相同，以比距离(R/a)作为参数，同一波形在不同范围区间上的规律也是有差异的，多种复杂的衰减规律综合在一起就造成了水下爆炸近场压力状态的复杂性。

水下爆炸近场内的波系结构主要决定于炸药内传播的爆轰波，使用CJ爆轰模型计算半径0.1 m (W=6.83 kg)的TNT球形炸药在中心起爆的波系，爆轰波由前导冲击波和之后的稀疏波区组成，由计算结果结合波的传播理论，得到图3中的波系结构。需要注意的是，图中以实线表示的激波、爆炸产物内的稀疏波区在向左的中心稀疏波到达之前，以及反射稀疏波波头IU线为计算所得，其他虚线表示的稀疏波区则为示意图。也就是说，图3中I、M、U、S的位置为计算所得，爆炸产物和水界面上的K、L点的确切位置则是未知的，而J点的实际横坐标R/a＜1.1，为了从波系图中清晰描述反射压缩波的汇聚过程，才标示在图中位置。冲击波OI在I点(爆炸产物和水的界面)形成向水内的透射冲击波IC和向爆炸产物传播的反射中心稀疏波区IUS，能够汇聚到O点的左传稀疏波会形成反射稀疏波向右传播。爆炸产物内的稀疏波区OIK在界面IK段上都会形成向水内的透射稀疏波，而由于中心稀疏波区IUS的作用，界面左侧爆炸产物内的声阻抗迅速降低，由高于界面右侧水的声阻抗降低到低于水的声阻抗，因此界面IK段会分为两段，IJ段上反射波为压缩波，JK段上反射波为稀疏波。IJ段上的反射压缩波区在追赶稀疏波区IUS时，压缩波的波头会不断进入稀疏波的波尾形成一条包络线IS，同时压缩波的波尾特征线会不断汇聚到前方的压缩波特征线上也形成一条包络线JS，IS和JS的距离(压缩波宽度)越来越小，形成二次冲击波，二次冲击波汇聚到O点后反射形成反射冲击波SM。

图3 水下爆炸近场内的波系Fig.3 Waves in near-field underwater explosion

6 峰值压力公式

本文中取水深H=0，使用平面一维、平面二维和二维轴对称格式的SPH方法计算半径a=0.1 m的炸药水下爆炸，平面一维、平面二维和二维轴对称不同坐标系中水下爆炸的峰值压力曲线如图4所示。虽然维数不同，冲击波传入水中的初始压力较为相近，都在15 GPa左右，峰值压力随距离增大而降低，峰值压力在同等距离上降低的数值大体上随距离增大而减小。由于几何效应，峰值压力衰减在高维中比低维更快。

对于6≤R/a≤12范围的峰值压力，本文中已经将二维轴对称峰值压力的计算结果与经验公式进行了比较，符合幂次关系。由于维数效应，越高维中压力衰减越快，各波形影响区域越窄，三维峰值压力中以幂次变化的区域起点R/a=6，R/a跨度6～12，因此低维中相应波形影响的峰值压力幂次变化区域起点R/a≥6，跨度不小于6。对R/a≥6的平面一维、平面二维峰值压力计算结果进行统计拟合，可以发现，的确有这样每点误差都很小的区间：

图4 峰值压力曲线Fig.4 Peak pressure curves

平面二维：

平面一维：

将平面一维的拟合范围10≤R/a≤30与图4进行对比可知，R/a=10正对应于压力峰值曲线中平台的末端，也就是图3中的L点，据此认为在平面二维、二维轴对称中，幂次拟合范围的起点R/a=6对应于L点也是合理的。由此可进行以下推论：在平面二维、二维轴对称水下爆炸问题中，R/a＜6的近场区域内，存在1＜β＜6，在1＜R/a＜β的范围内，压力峰值受IJ段影响，降低速度要快于6＜R/a＜12区域，β＜R/a＜6范围，则符合维数效应的幂次关系。

赵继波等[6]的研究结果认为，在1＜R/a＜4.7(140 mm/30 mm)的近场内采用指数形式能够较好地描述峰值压力的规律，与本文中认为的在1＜R/a＜β的一定范围内，峰值压力的降低速度要高于幂次形式是相符合的。而在1＜R/a＜6的区域上采用统一形式进行分段拟合，还需要考虑到β＜R/a＜6(文献[6]中β=4.7)范围内幂次形式与指数形式的较大差别。

本文中二维轴对称计算结果使用线性回归方法拟合曲线，根据最小二乘法，得到lnPm-ln(R/a)曲线拟合直线y=r+sx的系数估计r=22.92，s=-2.08，同样方法得到lnPm-R/a拟合直线的系数估计r=23.48，s=-0.84。在线性回归理论中，用R2表示直线拟合曲线的符合程度，R2越接近于1表示曲线越近似为直线。

测定系数

式中：样本平均值为样本数。通过计算可得到lnPm-ln(R/a)的拟合直线R2=0.983，lnPm-R/a的拟合直线R2=0.901。lnPm-ln(R/a)的R2更接近于1，说明在R/a＜6的近场范围内，采用幂次形式的拟合公式能够比指数形式的拟合公式产生更小的偏差。其拟合直线为：

式中：峰值压力单位为GPa，拟合结果与文献结果[4]是相近的。

在图5中lnPm-ln(R/a)计算曲线与拟合直线相比具有一定的下凸，这也意味着总体的衰减速度要快于幂次规律。对lnPm-ln(R/a)曲线进行幂次拟合能得到更加接近于计算曲线的拟合结果，图6中新的曲线形式为：

图5 近场中的lnPm-ln(R/a)和lnPm-R/a曲线Fig.5 lnPm-ln(R/a) curves and lnPm-R/a curves in near-field

前文已经指出，R/a＜6范围内的拟合曲线需要分为两段，而图6中的新结果在ln(R/a)＜0.5范围与数值结果曲线有较大差异。图7中分别为lnPmln(R/a)曲线的一阶和二阶导数，从曲线中可以看到在ln(R/a)=0.45，也就是R/a=1.57附近为明显的拐点，这与文献[4]中选择的R/a=1.52作为分界点是较为相近的。

图8以R/a=1.5为界，进行了分段拟合，最后的拟合形式能够使拟合曲线与计算曲线更加符合。对1＜R/a＜1.5范围的拟合，一方面计算中这部分的误差较大，计算结果只能作为参考，另一方面在拟合形式上也是为了与1.5＜R/a＜6内的拟合形式保持一致，因而没有采用指数形式。

图6 近场中lnPm-ln(R/a)的拟合曲线Fig.6 Fitted curve of lnPm-ln(R/a) in near-field

图7 近场中lnPm-ln(R/a)曲线的一阶和二阶导数Fig.7 First- and second-order derivatives of lnPm-ln(R/a) curve in near-field

图8 近场中lnPm-ln(R/a)的分段幂次拟合曲线Fig.8 Fitted curves of piecewise power of lnPm-ln(R/a) in near-field

7 结 论

用SPH方法进行了水下爆炸的数值模拟，通过自由场气泡、相似率以及中场峰值压力与试验结果的符合，验证了程序的有效性。进行了水下爆炸近场内波的传播分析，通过波系与数值结果的对比，分析了水下爆炸近场内的压力状态并采用新的形式进行了公式拟合。得到以下结论：

(1)R/a=6是水下爆炸近场波形变化的一个分界点。

(2)R/a＜6区域的峰值压力衰减速度分为明显的两段，R/a＜1.5的范围内，水下爆炸峰值压力降低的速度非常快。

(3) 使用幂次函数对曲线进行分段拟合，能够在水下爆炸近场峰值压力公式拟合中得到与数值结果更加符合的曲线。