Numerical simulation of peak pressure in near-field underwater explosion
-
摘要:
为了研究水下爆炸近场内的压力状态,使用SPH (smoothed particle hydrodynamics)方法,采用C-J爆轰模型,对水下爆炸过程进行了数值模拟。通过与峰值压力规律和中场内经验公式的比较,验证了程序的可靠性。对水下爆炸过程进行了波系分析,与不同维数水下爆炸的数值模拟结果进行了对比研究。结果表明:比距离R/a=6是波形变化的一个分界点,在R/a<6的近场范围内,峰值压力的拟合分为两段更合适。另外,对数值结果的lnPm-ln(R/a)曲线进行了分段幂次拟合,得到了与数值模拟结果非常吻合的拟合曲线。
Abstract:In this work, to obtain the pressure state of underwater explosion in near-field, we numerically simulated the whole process of underwater explosion using the smoothed particle hydrodynamics and adopting the C-J detonation model, and following the empirical formulas, confirmed the laws of the peak pressure, thereby verifying the effectiveness of the numerical program. We also analyzed the waves' propagation in underwater explosion and compared it with the numerical results of underwater explosion in various dimensions. The results show that the distance ratio R/a=6 is a demarcation point in the waves structure, and in the R/a<6 near-field range the fitted peak pressure curve should be divided into two sections. Further, we performed segmented fitting of the numerical results with power function, and found the fitted curve in good agreement with the numerical results.
-
Key words:
- underwater explosion /
- peak pressure /
- near-field /
- smoothed particle hydrodynamics /
- dimension effect /
- power fitting
-
压力状态研究是水下爆炸研究中的重要内容。以往水下爆炸压力状态研究的距离范围主要集中于中远场,现有的中远场水下爆炸研究结果与实验符合较好。水下爆炸近场压力则具有典型的高压高频特点,压力峰值至少在GPa量级,持续时间则为ms量级,给实验测量造成很大困难[1]。近期,学者们对近场峰值压力进行了试验和数值研究[2-8],认为在R/a≥1的一定范围内峰值压力的变化更符合指数规律,而在1≤R/a≤6范围进行拟合时则都使用了幂次形式,最后的拟合公式也不尽相同。这些试验和计算工作都证实,水下爆炸冲击波中存在相似原理,因此用数值模拟对水下爆炸近场中不同波系的影响进行一些规律性的分析是可行的。
本文中利用SPH方法进行水下爆炸近场压力的数值模拟,对模拟结果进行基于波系传播规律的定性分析,总结近场内峰值压力的变化规律,以期提高对水下爆炸近场压力变化规律的认识。
1. SPH方法
SPH方法[9-11]对物理量进行局部空间的积分近似和粒子近似,是一种拉氏无网格粒子类方法,在水下爆炸领域得到了越来越广泛的应用。本文计算中忽略黏性、热传导,使用的SPH控制方程为欧拉方程组:
dρidt=ρiN∑j=1mjρjvβij⋅∂Wij∂xβi (1) dvαidt=−N∑j=1mjpi+pjρiρj∂Wij∂xβi (2) deidt=12N∑j=1mjpi+pjρiρjvβij⋅∂Wij∂xβi (3) 式中:
ρ 、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取为最常用的三次样条函数形式:
W(R,h)=bd×{2/3−R2+R3/20≤R<1(2−R)3/61≤R<20R≥2 (4) 式中:
bd 在不同维数中取不同的值,一维b1=1/h ,二维b2=15/(7πh2) ,三维b3=3/(2πh3) 。在轴对称的柱坐标系统中,欧拉方程组的相应形式[11]为:
∂ηi∂t=N∑j=1mj(˙si−˙sj)⋅∇iWij (5) ¨ri=2πpiηi−2πN∑j=1mj(piriη2i+pjrjη2j)∂Wij∂ri (6) ¨zi=−2πN∑j=1mj(piriη2i+pjrjη2j)∂Wij∂zi (7) deidt=−2πpiηi˙ri+πN∑j=1mj(piriη2i+pjrjη2j)(˙si−˙sj)⋅∇iWij (8) 式中:
s(r,z) 为粒子在柱坐标中的位置,˙s=v ,η=2πrρ 。计算中没有考虑重力,为了模拟水深,在水的初始状态中增加了初始压力,比如若水深H=0,则取初始压力
p0=1.01×105 Pa。2. CJ爆轰模型
在计算炸药爆轰过程时,使用了经典的CJ爆轰理论模型[12]。假定炸药微元在通过爆轰波阵面时瞬间转变为爆炸产物并完成放能,在这种假定下,反应区的宽度变为零,成为以爆速传播的强间断面。假设波阵面的传播速度为D,产物的粒子速度为u1,在自持定常爆轰状态下,产物离开波阵面的速度D−u1等于当地声速c1:
D=u1+c1 (9) 通常认为爆速D保持不变,需要由实验测定。
3. 状态方程
众多学者对水下爆炸中适用的水的状态方程进行了数值模拟与试验的对比研究,认为在PH<32.0 GPa,
ρ <2.02 g/cm3范围内[13],多项式形式的状态方程能够适用。本文中计算的TNT药包水下爆炸算例中,水中的压力峰值也在此范围内,因此水的状态方程选用了多项式形式:p={A1μ+A2μ2+A3μ3+(B0+B1μ)ρ0EMμ≥0T1μ+T2μ2+B0ρ0EMμ<0 (10) 式中:
ρ0 为初始密度,μ 为压缩度,μ=ρ/ρ0−1 ,EM是单位质量内能增量,A1、A2、A3、T1、T2为压强量纲拟合系数,B0、B1为无量纲系数。其具体参数如下:A1=2.2 GPa,A2=9.54 GPa,A3=14.57 GPa,T1=2.2 GPa,T2=0,B0=0.28,B1=0.28。
爆炸产物的状态方程选取JWL形式:
p=A(1−ωR1V)e−R1V+B(1−ωR2V)e−R2V+ωVE (11) 式中: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]在大量理论研究与试验数据的基础上,提出了自由场中水下爆炸气泡脉动第一个周期和最大半径的经验公式:
TCole=k1W1/3(H+10.33)5/6 (12) RColem=k2W1/3(H+10.33)1/3 (13) 式中: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炸药,气泡第一次振动的周期TCole=0.152 s,最大半径
RColem =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的结果进行了深入和扩充,给出的峰值压力公式为:
Pzm={44.1×106×(W1/3R)1.56<R/a≤1252.4×106×(W1/3R)1.1312<R/a<240 (14) 式中:a为球形炸药的初始半径,R为测点到炸药中心的距离。
本文中采用SPH二维轴对称格式,取H=0,对球形装药进行爆轰模拟。比较6≤R/a≤12内的计算结果与经验公式,令ρ=1 630 kg/m3,将经验公式中球形药包的质量与半径进行了换算,如图2所示,两者的结果最大差别不超过2%。
从对自由场水下爆炸的计算结果看,本文使用的SPH方法得到的计算结果较好地与经验公式相符合,使用SPH方法进行水下爆炸近场内的压力研究是可行的。
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。
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的平面一维、平面二维峰值压力计算结果进行统计拟合,可以发现,的确有这样每点误差都很小的区间:
平面二维:
Pm=7.2×107×(W1/2/R)1.06≤R/a≤12 (15) 平面一维:
Pm=5.4×107×(W/R)0.8510≤R/a≤30 (16) 将平面一维的拟合范围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表示曲线越近似为直线。
测定系数
R2=1−∑[yi−(r+sxi)]2/∑(yi−¯y)2 (17) 式中:样本平均值
¯y=∑yi/n ,n为样本数。通过计算可得到lnPm-ln(R/a)的拟合直线R2=0.983,lnPm-R/a的拟合直线R2=0.901。lnPm-ln(R/a)的R2更接近于1,说明在R/a<6的近场范围内,采用幂次形式的拟合公式能够比指数形式的拟合公式产生更小的偏差。其拟合直线为:lnPm=22.92−2.08ln(R/a),或Pm=9.00(a/R)2.08 (18) 式中:峰值压力单位为GPa,拟合结果与文献结果[4]是相近的。
在图5中lnPm-ln(R/a)计算曲线与拟合直线相比具有一定的下凸,这也意味着总体的衰减速度要快于幂次规律。对lnPm-ln(R/a)曲线进行幂次拟合能得到更加接近于计算曲线的拟合结果,图6中新的曲线形式为:
lnPm=23.37−2.61(lnRa)0.761<R/a<6 (19) 前文已经指出,R/a<6范围内的拟合曲线需要分为两段,而图6中的新结果在ln(R/a)<0.5范围与数值结果曲线有较大差异。图7中分别为lnPm-ln(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内的拟合形式保持一致,因而没有采用指数形式。
lnPm=23.28−2.76(lnRa)0.921<R/a<1.5lnPm=23.74−3.01(lnRa)0.661.5≤R/a<6 (20) 7. 结 论
用SPH方法进行了水下爆炸的数值模拟,通过自由场气泡、相似率以及中场峰值压力与试验结果的符合,验证了程序的有效性。进行了水下爆炸近场内波的传播分析,通过波系与数值结果的对比,分析了水下爆炸近场内的压力状态并采用新的形式进行了公式拟合。得到以下结论:
(1) R/a=6是水下爆炸近场波形变化的一个分界点。
(2) R/a<6区域的峰值压力衰减速度分为明显的两段,R/a<1.5的范围内,水下爆炸峰值压力降低的速度非常快。
(3) 使用幂次函数对1<R/a<6的lnPm-ln(R/a)曲线进行分段拟合,能够在水下爆炸近场峰值压力公式拟合中得到与数值结果更加符合的曲线。
-
-
[1] 库尔. 水下爆炸[M]. 罗耀杰, 韩润泽, 官信, 等译. 北京: 国防工业出版社, 1960. [2] DORSETT H, CLIFF M D. Detonation front curvature measurements and aquarium tests of tritonal variants: ARML, DSTO-TR-1411 [R]. 2003. [3] BENTEROU J, BENNETT C V, COLE G, et al. Internal detonation velocity measurements inside high Explosives: LLNL-PROC-409969 [R]. 2009. [4] 池家春, 马冰. TNT/RDX(40/60)炸药球水中爆炸波研究 [J]. 高压物理学报, 1999, 13(3): 199–204. DOI: 10.11858/gywlxb.1999.03.008.CHI Jiachun, MA Bing. Underwater explosion wave by a spherical charge of composition B-3 [J]. Chinese Journal of High Pressure Physics, 1999, 13(3): 199–204. DOI: 10.11858/gywlxb.1999.03.008. [5] 赵继波, 谭多望, 李金河, 等. TNT药柱水中爆炸近场压力轴向衰减规律 [J]. 爆炸与冲击, 2008, 28(6): 539–543. DOI: 10.3321/j.issn:1001-1455.2008.06.010.ZHAO Jibo, TAN Duowang, LI Jinhe, et al. Axial pressure damping of cylindrical tnt charges in the near underwater-explosion field [J]. Explosion and Shock Waves, 2008, 28(6): 539–543. DOI: 10.3321/j.issn:1001-1455.2008.06.010. [6] 师华强, 宗智, 贾敬蓓. 水下爆炸冲击波的近场特性 [J]. 爆炸与冲击, 2009, 29(2): 125–130. DOI: 10.3321/j.issn:1001-1455.2008.06.010.SHI Huaqiang, ZONG Zhi, JIA Jingbei. Short-range characters of underwater blast waves [J]. Explosion and Shock Waves, 2009, 29(2): 125–130. DOI: 10.3321/j.issn:1001-1455.2008.06.010. [7] 张远平, 李金河, 龚晏青, 等. 水下爆炸近场冲击波压力测试研究 [J]. 仪器仪表学报, 2009, 30(6): 58–61.ZHANG Yuanping, LI Jinhe, GONG Yanqing, et al. Measuring study on shock wave pressure at near-field during underwater explosion [J]. Chinese Journal of Scientific Instrument, 2009, 30(6): 58–61. [8] 李晓杰, 李现远, 张程娇, 等. 水下爆炸近场冲击波速度连续测试 [C] // 中国力学大会, 2013.LI Xiaojie, LI Xianyuan, ZHANG Chengjiao, et al. Continuous velocity measurement of underwater explosion shock wave in the near field [C] // The Chinese Congress of Theoretical and Applied Mechanics(CCTAM2013), 2013. [9] LIU G R, LIU M B. Smoothed particle hydrodynamics: a meshfree particle method [M]. translated by HAN Xu, YANG Gang, QIANG Hongfu, et al. Changsha: Hunan University Press, 2005. [10] ZHANG Aman, YANG Wenshan, HUANG Chao, et al. Numerical simulation of column charge underwater explosion based on SPH and BEM combination [J]. Computers & Fluids, 2013(7): 169–178. DOI: 10.1016/j.compfluid.2012.10.012. [11] BROOKSHAW L. Smooth particle hydrodynamics in cylindrical coordinates [J]. Australian and New Zealand Industrial and Applied Mathematics Journal, 2003, 44(E): C114–C139. DOI: 10.21914/anziamj.v44i0.675. [12] 李维新. 一维不定常流与冲击波[M]. 北京: 国防工业出版社, 2003. [13] 李晓杰, 张程娇, 王小红, 等. 水的状态方程对水下爆炸影响的研究 [J]. 工程力学, 2014, 31(8): 46–52. DOI: 10.6052/j.issn.1000-4750.2013.03.0180.LI Xiaojie, ZHANG Chengjiao, WANG Xiaohong, et al. Numerical study on the effect of equations of state of water on underwater explosions [J]. Engineering Mechanics, 2014, 31(8): 46–52. DOI: 10.6052/j.issn.1000-4750.2013.03.0180. [14] ZAMYSHLYAYEV B V, YAKOVLEV Y S. Dynamic loads in underwater explosion[M]. Washington, D. C. AD-757183,1973. 期刊类型引用(12)
1. 刘先辉,王高辉,卢文波. 水下双弹爆炸冲击波传播特性研究. 武汉大学学报(工学版). 2024(01): 28-36 . 百度学术
2. 严泽臣,岳松林,邱艳宇,王建平,赵跃堂,施杰,李旭. 水下爆炸冲击波反射压力计算方法的改进. 兵工学报. 2024(04): 1196-1207 . 百度学术
3. 宗周红,甘露,院素静,李明鸿,单玉麟,林津,夏梦涛,陈振健. 桥梁结构抗爆安全防护研究综述. 中国公路学报. 2024(05): 1-37 . 百度学术
4. 王丕光,卢冉冉,闫秋实,李述涛,杜修力. 水下爆炸作用下基于声学的爆源子结构输入方法. 力学学报. 2023(04): 915-924 . 百度学术
5. 任凯,周洪景,杨晨. 船体水下近距非接触爆炸损伤计算之两步迭代法. 爆炸与冲击. 2023(04): 90-101 . 本站查看
6. 刘哲函,王晓明,王燕,李健,南德,刘泽玉,曾志. WIGWAM水下核爆炸气泡脉动理论计算及数值模拟对比分析. 现代应用物理. 2023(03): 244-252 . 百度学术
7. 尹国福,陈建华,任西,任炜,李蛟,纪向飞. 反蛙人榴弹引战系统集成设计与毁伤规律研究. 火工品. 2022(01): 1-5 . 百度学术
8. 杨文,岳彩新,宋家良,吴志超,孙晨. 工业电子雷管抗冲击性能试验研究. 火工品. 2022(02): 16-19 . 百度学术
9. 杨文,吴竞,宋家良,夏光,黄孟文,陈以钻. 两种电子雷管抗冲击性能对比试验研究. 爆破器材. 2022(05): 38-42 . 百度学术
10. 岳彩新. 两种提高工业电子雷管抗冲击波性能缓冲材料的实验研究. 煤矿爆破. 2021(03): 1-5 . 百度学术
11. 林尚剑,王金相,马腾,秦健,李恒,黄瑞源,刘亮涛. 水下多点爆炸冲击波叠加效应研究. 兵工学报. 2020(S1): 39-45 . 百度学术
12. 邱清水,陈莹玉,古滨,李炳南,姚熊亮,王志凯. 水下近场爆炸载荷数值预报研究. 四川轻化工大学学报(自然科学版). 2020(05): 44-50 . 百度学术
其他类型引用(15)
-