水下爆炸冲击波数值仿真精度研究

2024-03-14敖启源姜智雅康珀阁

水下无人系统学报 2024年1期

敖启源 ,卢熹 *,姜智雅 ,康珀阁

(1.沈阳理工大学装备工程学院,辽宁沈阳,110159;2.山西江阳化工有限公司,山西太原,030041)

0 引言

水下武器作为舰船生命力的主要威胁之一,其爆炸冲击波及气泡载荷会对舰船造成严重的局部和总体破坏[1]。各国相继开展了大量水下爆炸试验,但实弹实测试验安全风险较高、成本高昂、重复性低且观测范围有限[2]。随着计算机性能和仿真技术的发展,数值仿真以其较高计算精度、低成本和高可重复性等优点成为继实验、理论研究后第3 种水下爆炸研究方法[3]。因此开展水下爆炸数值仿真精度研究具有十分重要的意义。

数值仿真研究水下爆炸冲击波载荷问题时,为了处理冲击波的强间断面,抑制波阵面前后的数值振荡,引入了人工粘性。由于人工粘性的引入会在几个网格宽度上光滑冲击波波阵面,因此水下爆炸冲击波载荷计算结果直接受人工粘性的影响。同时，人工粘性的引入要求网格尺寸不能过大,否则计算过程中网格会忽略部分波阵面信息,致使峰值过低影响计算精度。Huang 等[4]通过典型三硝基甲苯(Trinitrotolue,TNT)炸药水下爆炸数值分析,探讨了一次、二次粘性系数对数值计算结果的影响,并给出了一定比例爆距范围内的建议值。Wang 等[5]研究了不同装药质量下网格尺寸对冲击波峰值的影响,并引入与装药半径和单元边长相关的无因次量表征网格尺寸。胡亮亮等[6]以常规TNT 水下爆炸为例,对水的状态方程、人工粘性系数和网格尺寸对于数值仿真结果的影响进行了研究。张社荣等[7]基于有限元软件AUTODYN建立了不同炸药量的水下爆炸数值模型,对比分析了网格尺寸对不同爆距处冲击波峰值压力及比冲量的影响。此外还有其他学者[8-13]讨论了网格尺寸及粘性系数对水下爆炸计算结果的影响,但现有研究无法在预定精度下快速确定网格尺寸和人工粘性系数。因此文中以TNT 水下自由场爆炸数值计算为例,探究网格尺寸和一次项人工粘性系数对水下爆炸冲击波峰值压力的影响,综合分析不同工况下网格密度因子和一次项系数与峰值压力平均误差间的关系,构建出普适性较高的水下爆炸数值误差预估模型,为预定精度的仿真模型设计提供依据。

1 水下爆炸数值仿真

1.1 数值模型

基于文献[14]中开展的爆炸水井试验,建立二维轴对称计算模型。试验水域尺寸4.5 m×9 m,空气域尺寸4.5 m×0.1 m,炸药为直径2 cm 的等高药柱,炸药质量78 g,放置在水深4 m 位置处,采用中心起爆方式。网格的排列走向和过渡方式会影响计算结果,因此采用均匀网格划分方式。在距装药中心0.4～2.8 m/kg1/3比例爆距内选定16 个观测点,采用关键字*INITIAL_HYDROSTATIC_ALE初始化静水压力,设置水域压力梯度,以模拟真实条件下的水下压强环境。炸药、空气和水介质均选用ALE（Arbitrary Lagrangian-Eulerian）算法,数值计算模型如图1 所示。

图1 数值计算模型Fig.1 Numerical calculation model

1.2 数值模型

装药选用典型的单质TNT 炸药,采用JWL(Jones Wilkins Lee)状态方程描述其爆炸过程,具体形式为

式中:V为相对体积;A、B、R1、R2为常数,取值如表1[15]所示;P为压力;ω为药量;E为单位体积内能。

表1 TNT 状态方程参数Table 1 State equation parameters for TNT

当水介质处于压缩状态时,其状态方程为

当水介质处于膨胀状态时,其状态方程为

式中:µ为水的压缩比;C为水中声速;S1、S2和S3为常数;γ0为GRUNEISEN 系数;α为体积修正系数;V0为初始相对体积。以上参数取值如表2所示。

表2 水状态方程参数Table 2 State equation parameters for water

空气使用POLYNOMIAL 状态方程进行描述,其形式为

式中:C0、C1、C2、C3、C4、C5和C6为常数。取值如表3 所示。

表3 空气状态方程参数Table 3 State equation parameters for air

1.3 冲击波峰值压力计算结果

根据文献[10],长径比为1∶1 的柱形装药可以近似为球形装药,对于球形装药水下自由场爆炸冲击波的传播,Cole[16]通过大量试验标定了水下爆炸冲击波相似律的公式系数,获得了不考虑水深影响的TNT 炸药水下爆炸冲击波峰值压力计算公式,Zamyshlyaev[17]在其基础上将经验公式修正为

式中:Pm为冲击波峰值压力;R为爆距;W为装药质量;Re为装药半径。

实时监测柴油发电机的运行状态、报警信息及运行参数，并对柴油发电机的输出电源的质量进行监测。在动力机房室，通过发电机厂家提供的专用智能通讯接口，及相应的通讯协议，软件工程师将根据此编集驱动，导入系统平台进行实时监测；实时地监视柴油机的起停状态，各相电压，电流，频率等参数进行监测。

根据式(5)可计算出水中一定范围内的冲击波峰值压力,即表4 中的经验值。为获得水下爆炸冲击波传播演化规律,文献[14]在爆炸水井中开展了水下爆炸冲击波试验,获得了78 g TNT 装药不同爆距Z处的冲击波峰值压力,即表4 中的试验值。基于爆炸水井试验,利用LS-DYNA 有限元软件对78 g TNT 炸药水下爆炸过程进行仿真,网格尺寸0.25 cm,一次项系数取值0.06,获得不同爆距处的冲击波峰压力,即数值解。将试验值、数值解及经验值进行对比。

表4 不同爆距处峰值压力对比Table 4 Comparison of peak pressure at different scaled blast distances MPa

对比表4 数据,可以发现数值解与试验值最小误差仅2.6%,平均误差6.9%;与经验值最小误差4%,平均误差9.2%。因此网格尺寸和一次项粘性系数取值合理时,数值模型可获得较高的计算精度。经验值与试验值峰值压力平均误差7.2%,表明经验公式可以较为准确地预估水下爆炸冲击波峰值压力。因此,在探究网格尺寸及人工粘性对数值计算精度的影响时,仅比较数值解与经验值。

2 网格尺寸及人工粘性对Pm 的影响

在有限元计算中为了反映冲击波波阵面的强间断,引入了人工粘性来光滑冲击波,这使得数值计算峰值压力低于真实值。一次项人工粘性系数对冲击波峰值压力影响较大,二次项人工粘性系数主要用于抑制冲击波衰减过程中的虚假振荡,但对水中爆炸数值仿真计算中冲击波峰值影响较小。参阅文献[8],二次项人工粘性系数取定值1.0。LS-DYNA 中人工粘性形式为

式中:ρ为材料密度;Q1为二次项人工粘性系数;Q2为一次项人工粘性系数;l为特征长度;C为当地声速;ε为体积变化率。

在水下爆炸数值计算中,网格密度和排列方式对计算结果影响很大,过大的网格尺寸在计算过程中会忽略冲击波波阵面信息,冲击波爬升至峰值所需时间变长;而过小的网格尺寸对计算资源带来的压力也不可忽视。同时,水域网格长宽比尽量接近于1,尤其是炸药附近的网格,否则爆炸冲击波易出现失真[3]。在探究网格尺寸对数值计算结果的影响时,引入与装药半径R0和网格尺寸L0相关的无因次量,因此文中为研究网格尺寸对数值计算结果的影响[5,7-8],其具体形式为

为探究水下爆炸数值仿真中网格尺寸和Q2对冲击波峰值压力的影响,对78 g TNT 水下自由场爆炸过程进行数值仿真计算。当Q2＞0.1 时,计算峰值压力误差较大[8],因此调整分别取Q2=0.02、0.03、0.04、0.05、0.06、0.07、0.08 和0.10。对于二维模型而言,网格密度因子λ=8 时便满足大部分二维数值模型,因此取λ分别取λ=1、2、3、4、5、6、7 和8,共计算64 个工况。为了从全局意义上分析网格尺寸和Q2对Pm的影响,引入平均误差,定义为比例爆距0.375～2.8 m/kg1/3范围内选定观测点峰值压力相对于经验公式误差的平均值。

网格尺寸及Q2对Pm结果如图2～图6 所示。由图2 和图3 可知,当λ较小时,Pm较低,与经验公式值偏差较大,随着λ的增加,Pm显著增大且近场峰值压力与经验值吻合较好。但随着λ增加,Pm趋于稳定,不会随着λ的变化出现显著变化,这表明随着网格密度的增加,Pm对网格的敏感性越来越低。与此同时,网格尺寸对冲击波超压爬升至峰值所需时间影响很大,λ越小,冲击波超压爬升至峰值时间越长且峰值压力越低,如图4 所示。

图2 Q2=0.08 时,不同λ 时Pm 随比例爆距变化曲线Fig.2 Variation of Pm with scaled blast distance for different λ at Q2=0.08

图3 Q2=0.08 时,不同比例爆距处Pm 随λ 变化曲线Fig.3 Variation of Pm with λ for different scaled blast distances at Q2=0.08

图4 λ 对Pm 上升速度影响Fig.4 Effect of λ on the rate of increase of Pm

由图5 和图6 所见,随着Q2逐渐变小,Pm逐渐增大,且近场冲击波峰值压力受Q2影响较大。由于Q2取值范围较小,不同比例爆距处Pm随Q2的变化曲线与文献[8]中有所不同。在实际计算中发现,当网格密度较大时,过小的Q2反而导致Pm数值解与经验值的偏差增大,如表5 所示。通过计算,Q2=0.02 时不同爆距处的平均误差为11.2%,而Q2=0.06 时平均误差仅4.2%,因此有必要探究λ和Q2对峰值压力平均误差的影响。

表5 λ=6 时,不同比例爆距处冲击Pm 对比Table 5 Different scaled blast distances peak pressure of shock waves at λ=6 MPa

图5 λ=3 时,不同Q2 时Pm 随比例爆距变化曲线Fig.5 Variation of Pm with scaled blast distance for different Q2 at λ=3

图6 λ=3 时,不同比例爆距处Pm 随Q2 变化曲线Fig.6 Variation of Pm with Q2 for different blast scaled distances at λ=3

对比数值解和经验值在不同爆距处的冲击波峰值压力,获得不同Q2下峰值压力平均误差EP随λ的变化关系,如图7 所示。整体来看,随着λ的增大和Q2的减小,EP逐渐降低。但以Q2=0.10为例,随着λ的增大,EP先减小后增大,这是因为网格密度较大时,过小的Q2会加大伪振荡,致使Pm数值解与经验的相对误差增大,因此在数值计算中不能一味地减小网格尺寸和一次项系数。同时,当λ=1 时,EP很大,即使调整Q2也未能使EP满足工程精度。因此在水下爆炸数值计算中应首先确定λ,同时调整Q2方能得到精度较高的计算结果。

图7 峰值压力平均误差Fig.7 Average errors of peak pressure

3 计算误差预估模型

为了便于应用,实现在预定精度下快速确定网格尺寸和一次项系数,需构建出关于λ和Q2的误差EP预估模型。考虑到实际工程需要,误差应控制在20%以内。同时,由前面分析可知,较大的网格密度以及过小的一次项系数均可能导致计算误差增大。因此,将误差预估模型的变量区间范围限定在0.03≤Q2≤0.10 和3≤λ≤8,则图7 中的数据在该区间范围显示为图8 中曲线。

图8 限定区间后峰值压力平均误差曲线Fig.8 Average errors of peak pressure curves after limited interval

将图8 中的EP和λ取对数得到lgEP和lgλ的关系曲线如图9 所示,k为斜率,b为截距。可以看到,各曲线呈近似平行的线性关系,因此可构造线性表达式为

图9 lgEP 和lgλ 关系曲线Fig.9 The relationship curves between lgEP and lgλ

利用式(8)对图9 中的数据进行拟合得到拟合系数如表6 所示。其中,复合相关系数R1均值0.979,决定系数R2均值0.959,可见数据线性拟合精度较高。表6 中各曲线斜率k值比较接近,为了得到归一化的误差预估模型,可以取k值平均值,并且将截距b看作关于Q2的函数。如图10 所示为截距b与lgQ2的关系,可以看到二者近似呈线性关系,因此可构建线性关系式

表6 lgEP 关于lgλ 的拟合参数Table 6 Fitting parameters for λ in EP

图10 截距-对数粘性系数线性拟合Fig.10 Linear fitting of intercept-logarithm viscosity coefficient

利用式(9)对图10 中数据进行线性拟合,结果如表7 所示,可以看出数据具有较高拟合精度。将式(9)代入式(8),整理得

表7 b 关于Q2 的拟合参数Table 7 Fitting parameters for b in Q2

从式(10)可以看到,误差EP是关于的函数。利用图8 中的数据绘制EP与的关系曲线如图11 所示。可以看到,EP与近似呈线性关系。因此,为了进一步提高误差预估模型的拟合精度,对图11 中数据重新进行线性拟合,其中,拟合模型的复合相关系数R=0.998,决定系数R2=0.996,说明拟合精度较高,则式(11)为最终得到的误差预估模型

图11 峰值压力平均误差拟合结果Fig.11 Fitting results of peak pressure average error

为验证误差预估模型具有较高的普适性,需要对不同药量和装药形状的炸药进行仿真计算,同时结合研究背景,对0.2、5、500、1 500 和5 000 kg TNT 柱形装药(长径比为1)和球形装药水下爆炸进行计算。在距装药中心0.327～2.8 m/kg1/3比例爆距间设置16 个观测点,取Q2=0.06,λ=6。不同工况下的峰值压力平均误差及预估误差如表8 所示。可以看到,预估误差与实际计算误差相近,因此误差预估模型对于不同装药量和装药形状也具有很高的适用性,有助于建立水下爆炸模型对整体计算精度和网格数量进行综合分析,为预定精度的仿真模型设计提供依据。