Numerical simulation of intense blast wave reflected on rigid ground
-
摘要: 为准确预测空中强爆炸产生的冲击波载荷分布,基于欧拉坐标系建立了能够模拟具有高密度比、高压力比的强激波问题的二维多介质流体数值方法。结合网格自适应技术,对1 kt TNT当量的空中强爆炸在不同爆炸高度下的冲击波地面反射过程进行了数值模拟,并考虑了真实气体状态方程和空气随高度不均匀分布的影响。计算得到了地面上距爆心投影点大尺度范围内的反射超压和冲量等冲击波载荷分布,并给出了冲击波载荷随爆高的变化规律。Abstract: In order to predict the intense blast wave in air, a multi-material numerical scheme is proposed in two-dimensional cylindrical coordinates on Eulerian grids, which can handle the blast wave problems with high ratio in initial density and pressure. Combined with the adaptive mesh technique, the propagation of blast wave produced by a 1 kt TNT intense explosion is simulated, and the effects of real gas equation of state and the nonuniform atmosphere are taken into account. The calculated blast wave parameters on the ground, such as peak overpressures and impulses, agree well with the experimental data in a wide space range, and the influences of heights of burst are analyzed.
-
Key words:
- air blast /
- multi-material flow /
- mesh adaption /
- shock reflection /
- height of burst
-
炸药在距地面有限高度内爆炸时,冲击波将会在地面、建筑等障碍物的表面发生反射和绕射现象,传播规律变得更加复杂。相比于化学爆炸,强爆炸产生的初始压力可达上万吉帕,初始温度高达数十万摄氏度[1-2],空气在冲击波作用下将发生剧烈的压缩和温升,热力学性质发生显著变化,需要建立高温高压下真实气体的状态方程。此外,由于强爆炸冲击波的传播距离远,其在传播过程中会受到空气不均匀性的影响,空气的初始压强和密度越低,冲击波的衰减越明显[2-3]。空中强爆炸问题具有的上述特点,对数值计算的精度和效率提出了很高的要求。
空中强爆炸问题属于可压缩多介质大变形问题的范畴,一般会存在多种介质的相互作用。由于不同介质的状态方程和初始条件通常存在较大的差异,会给数值计算特别是介质界面附近的数值模拟带来很大的困难。目前针对空中爆炸冲击波问题的主流数值方法是基于Euler坐标系下的多介质流体数值方法,其包含两个关键步骤:界面追踪和介质间相互作用的处理。其中,常见的界面追踪方法包括VOF (volume-of-fluid)方法[4]、Level Set方法[5]和Front Tracking方法[6]等。当确定介质界面的位置后,需进一步准确描述不同介质间的相互作用,常见的处理方法包括虚拟流体类方法(ghost fluid method)[7-8]和切割单元法(cut cell method)[9]。
在爆炸冲击波地面反射规律研究方面,Crowl[10]给出了TNT等化学炸药在不同爆高下的地面反射冲击波载荷分布规律,乔登江[1]和Glasstone等[2]依据实验结果分别给出了标准大气状态下1 kt TNT当量的空中强爆炸冲击波在大尺度范围内的等超压曲线。目前公开发表的文献大多针对TNT等化学爆炸产生的冲击波在局部范围内的传播[11-15],较少见到能够完整给出空中强爆炸冲击波在地面反射过程的大尺度范围数值模拟。
针对空中强爆炸问题的强对流性,本文在前期工作[16]的基础上进一步改进基于欧拉坐标系的多介质流体数值方法,并考虑高温、高压下的真实气体状态方程和空气随高度不均匀分布的影响[17],结合网格自适应技术,计算1 kt TNT当量空中强爆炸产生的冲击波在不同爆高下的地面反射完整过程,得到地面上距爆心投影点5 km范围内的冲击波载荷分布。通过和实测结果进行比对,验证计算结果的正确性,并在此基础上研究地面反射冲击波载荷随爆高的变化规律。
1. 物理模型和数值方法
1.1 物理模型
根据空中强爆炸问题的对称性,采用二维柱对称计算模型,控制方程组为:
∂U∂t+∂F(U)∂r+∂G(U)∂z+H(U)=0 (1) U=(ρρuρvE),F=(ρuρu2+pρuv(E+p)u),G=(ρvρuvρv2+p(E+p)v),H=(ρu/rρu2/rρuv/r(E+p)u/r) (2) 式中:
ρ 为密度,u 、v 分别为r,z 方向的速度,E 为体积总能量,p 为压力。由于强爆炸伴随着强烈的放热、电离等效应,爆炸产物迅速膨胀至气体状态,且周围空气在受到冲击波、热辐射等因素的综合作用下发生剧烈的加热和压缩,物理状态极其复杂。本文中采用基于Saha电离平衡理论给出的、以数值表形式给出的状态方程[1]来描述强爆炸产物的热力学行为,并采用真实气体状态方程来描述空气:
p=(γ−1)ρe (3) 式中:
γ 为绝热指数,e 为比内能。Symbalisty等[18]给出了8组参考密度下(10−6~10 kg/m3)的温度与内能、温度与绝热指数之间关系的实测结果,如图1所示。通过对参考曲线进行双对数插值来获得各热力学变量之间的关系:首先,根据已知的比内能e,对图1(a)进行对数插值,得到相应的温度T;然后,对图1(b)进行插值,获得相应的绝热指数γ−1 ;最后,将γ−1 代入式(3)进行求解,最终完成整个状态方程的完整计算。实际空气的初始压力和密度随着高度的增加逐渐降低,且满足下列关系[1, 17]:
p0=pref(1−bTrefz)μgRb (4) ρ0=ρref(1−bTrefz)μgRb−1 (5) 式中:
Tref=288.16 K,pref=1.01325×105 Pa,ρref=1.225 kg/m3分别为标准大气温度、压强和密度;μ=28.966 g/mol,g=9.8 kg/m2,R=8.31 m2·s2/K;z为垂直高度,b为拟合系数。1.2 数值方法
采用基于欧拉坐标系、具有互不相溶特征的二维柱对称多介质流体计算模型来进行数值离散。在任意时刻,整个计算区域
Ω 可以分成2个子区域Ω+(t) 和Ω−(t) ,且满足:Ω+(t)∪Ω−(t)=Ω,Ω+(t)∩Ω−(t)=Γ(t) (6) 式中:
Γ(t) 为2种互不相溶流体之间的物质界面。在前期工作[16]的基础上,采用Level Set方法捕捉爆炸产物与空气的物质界面,根据Level Set函数分片线性的特征清晰重构出具有分片线性的物质界面,并通过在界面的法向上求解多介质Riemann问题的精确解来构造爆炸产物与空气之间的数值通量。整个计算流程如下。
1.2.1 物质界面的追踪和重构
Level Set方法[5, 19-20]把随时间运动的物质界面
Γ(t) 定义为距离函数φ(x,t) 的零等值面。利用特征线方法求解Level Set函数的演化方程:∂φ∂t+˜u⋅∇φ=0 (7) 式中:
˜u 为由物质界面的运动速度延拓得到的速度场[21]。根据特征线方法的思想,只要知道
tn+1 时刻网格顶点Xi 上的流体在tn 时刻的位置x(tn+1) ,就可以得到新的Level Set函数在Xi 上的值φ(x,tn+1) 。对每个顶点Xi ,利用特征线方法可得:x(tn+1)=x(tn)−(tn+1−tn)˜un (8) 再利用代数平均,可得新的Level Set函数:
φ(x,tn+1)=∑Xi∈τ(n)iφ(x,tn)/∑Xi∈τ(n)i1 (9) 式中:
τ(n)i 表示顶点Xi 所处于的单元。利用显式正系数格式[21]求解重新初始化方程:
{φt+sgn(φ0)(|∇φ|−1)=0φ(x,0)=φ0 (10) sgn(φ0)=φ0√φ20+ε2 (11) 以保持Level Set函数的符号距离性质,其中,
φ0 为重新初始化前φ(x,t) 的值,ε 为一很小的正数(例如10−6)。显示正系数格式的求解过程[21]比较繁锁,本文不再说细展开。当得到
tn 时刻的Level Set函数φn 后,根据其分片线性的特征重构出界面单元内分片线性流形的物质界面,如图2所示,图中nKi,n 为物质界面的单位法向。假设单元Ki,n 与单元Kj,n 具有共边Sij,n ,且包含物质界面ΓKi,n 。ΓKi,n 将Ki,n 和Sij,n 依次分为2个部分K±i,n 和S±ij,n ,此时该单元包含2种流体,且每种流体的守恒量分别定义为:U±Ki,n=(ρ±Ki,nu±Ki,nE±Ki,n) 式中:
ρ±Ki,n 、u±Ki,n 和E±Ki,n 分别为每种流体的密度、速度和体积比总能。1.2.2 单元边界数值通量
单元边界的数值通量是指单元边界上同种流体之间的数值通量,且满足:
ˆF±ij,n=Δtn|S±ij,n|ˆF(U±Ki,n,U±Kj,n;n±ij,n) (12) 式中:
Δtn 为时间步长,n±ij,n 为单元边界上的单位法向,ˆF(U±Ki,n,U±Kj,n;n±ij,n) 为相容、单调的数值通量函数,|S±ij,n| 表示单元的边长。本文中采用local Lax-Friedrich数值通量函数:
ˆF(U±Ki,n,U±Kj,n;n±ij,n)=12(F(U±Ki,n)+F(U±Kj,n))⋅n±ij,n−12λ(U±Kj,n−U±Kj,n) (13) 式中:
λ=max{λKi,n,λKj,n} ,λKi,n 、λKj,n 分别为U±Ki,n 和U±Kj,n 的最大信号速度。1.2.3 物质界面数值通量
物质界面数值通量是物质界面两侧不同流体之间的数值通量,且满足:
ˆF±Ki,n=Δtn|ΓKi,n|[0p∗Ki,nnKi,np∗Ki,nu∗Ki,n] (14) 式中:
|ΓKi,n| 为物质界面的长度;p∗Ki,n 和u∗Ki,n 分别为物质界面上的压力和法向速度,且通过在物质界面法向上求解局部一维多介质Riemann问题的精确解来获得,具体求解过程见文献[22]。1.2.4 守恒量的更新
当计算得到
Ki,n 的单元边界数值通量和物质界面数值通量后,可对该单元中每种流体的守恒量进行更新:U±Ki,n+1={0K±i,n+1=01|K±i,n+1|(|K±i,n|U±Ki,n+∑S±ij,n⊆K±i,nˆF±ij,n+ˆF±Ki,n)K±i,n+1≠0 (15) 式中:
U±Ki,n 、U±Ki,n+1 和K±i,n 、K±i,n+1 分别表示tn 和tn+1 时刻单元内每种流体的守恒量和体积。2. 空中强爆炸冲击波地面反射
2.1 空中强爆炸冲击波反射数值模拟
利用数值模拟程序对当量为1 kt TNT、爆高为100 m的空中强爆炸冲击波传播过程进行数值模拟。其中,强爆炸产物采用等温等压球模型,取爆炸总当量的85%作为爆炸产物的力学初始能量[1],产物质量设为50 kg,初始半径设为30 cm。空气采用真实气体状态方程,取水平面处的初始压力为101.3 kPa,初始密度为1.29 kg/m3,并按照式(4)和(5)对不同高度处空气的初始压力和密度进行设置。
计算区域在
r=0 和z=0 设为固壁反射边界条件,其余边界设为无反射边界条件。数值计算中采用h-网格自适应技术[23],根据相邻网格的压力梯度来进行网格的局部加密或稀疏,得到典型时刻的冲击波压力等值线以及相应的网格自适应图,如图3所示。由图3可知,当空中强爆炸产生的冲击波传播到地面时,最早在爆心投影点发生正反射,然后以逐渐增大的入射角发生斜反射,产生双激波结构的规则反射。随着冲击波向外继续传播,入射角不断增大。当入射角增大到临界角附近时将发生马赫反射,此时反射波阵面和入射波阵面的交点离开地面,形成垂直于地面的马赫杆。随着冲击波的进一步传播,马赫杆迅速增长。图4给出了数值计算得到的地面冲击波载荷分布(峰值超压和冲量)与实测结果[2]的对比情况,其中在距爆心投影点1 km以内,数值计算得到的峰值超压与实测结果符合较好,而在1 km以外则比实测结果略低,这是由于随着计算区域的尺度增大,网格自适应后的单元尺寸不够精细,使得数值计算对峰值超压的捕捉能力减弱导致。此外,数值计算得到的冲量比实测结果略大,但最大误差不超过50%。
2.2 不同爆炸高度的冲击波反射规律
利用数值模拟程序进一步对当量为1 kt TNT,爆高
H 依次为20、50、100、200、400和800 m的空中强爆炸冲击波传播过程进行了数值模拟,得到了地面上不同距离处的冲击波峰值超压和冲量,如图5所示。图5的计算结果表明,对1 kt TNT当量的空中强爆炸,在距爆心投影点300 m范围以内,地面上的峰值超压和冲量随着爆高的增大而逐渐减小。在距爆心投影点300 m范围以外,随着距爆心投影点的距离的增大,不同爆高条件下地面反射冲击波载荷的差距逐渐变小。其中,爆高约为200 m的空中强爆炸产生的地面峰值超压和冲量处于最大值,其他爆高条件下的冲击波载荷随着爆高的进一步增大或减小均呈现逐渐减弱的趋势,其原因主要是由于空中强爆炸地面反射而引发的复杂冲击波传播引起的。当空中强爆炸形成的冲击波达到地面时,随着入射角的增大,会依次在地面形成规则反射和马赫反射等复杂的波系结构。其中,在规则反射区,反射冲击波载荷由入射波强度和入射角决定。对于距爆心投影点相同距离的观察点,爆高越大,冲击波载荷越弱。
当入射角增大到临界角附近时,冲击波将会发生马赫反射。马赫反射的特点是,反射波阵面和入射波阵面的交点离开地面,形成垂直于地面的马赫杆。马赫反射对地面冲击波载荷分布起到了双重效果,一方面马赫反射的出现使得冲击波载荷的强度局部增大,另一方面由于马赫杆高度的不断增大,加快了冲击波波阵面能量的扩散速度,使得马赫反射区地面的冲击波载荷衰减速度加快。因此,空中强爆炸冲击波在地面反射后产生的大尺度空间载荷分布与爆炸当量以及爆炸高度密切相关,在实际应用中需要根据预定的目的,正确地选择爆炸方式,可以有效地扩大冲击波破坏区域。
3. 结 论
本文中建立了基于欧拉坐标系的互不相溶、具有清晰锐利界面的可压缩多介质流体数值方法,并研制了适用于模拟空中强爆炸冲击波传播的数值计算程序。利用Level Set方法捕捉爆炸产物与空气的运动界面,并清晰重构出分片线性流形的物质界面。采用有限体积方法求解每种物质的控制方程组,并通过在物质界面法向求解局部一维多介质Riemann问题的精确解来计算物质界面的数值通量。数值计算中采用了网格自适应技术,在捕捉激波峰值的前提下,有效提高了计算效率。
利用计算程序首先对空中强爆炸冲击波遇到刚性地面时的反射过程进行了数值模拟,得到了与实测数据一致的数值结果,验证了数值模拟程序的可靠性。在此基础上进一步对不同爆高条件下强爆炸冲击波的反射过程进行了数值模拟,给出了不同爆高条件下地面反射冲击波峰值超压和冲量在不同反射区域的分布情况,得到了冲击波载荷随爆高的变化规律,并基于数值模拟结果给出了强爆炸冲击波对典型目标结构毁伤破坏的最佳爆炸高度,为空中强爆炸冲击波的杀伤破坏效果预测以及爆炸方式的选取提供了研究基础。
-
-
[1] 乔登江. 核爆炸物理概论[M]. 北京: 国防工业出版社, 2003: 51−55. [2] GLASSTONE S, DOLAN P J. The effects of nuclear weapons [R]. USA: Defense Technical Information Center, 1977: 453−501. DOI: 10.21236/ada087568. [3] 段晓瑜, 崔庆忠, 郭学永, 等. 炸药在空气中爆炸冲击波的地面反射超压实验研究 [J]. 兵工学报, 2016, 37(12): 2277–2283. DOI: 10.3969/j.issn.1000-1093.2016.12.013.DUAN Xiaoyu, CUI Qingzhong, GUO Xueyong, et al. Experimental investigation of ground reflected overpressure of shock wave in air blast [J]. Acta Armamentarii, 2016, 37(12): 2277–2283. DOI: 10.3969/j.issn.1000-1093.2016.12.013. [4] HIRT C W, AMSDEN A A, COOK J L. An arbitrary Lagrangian-Eulerian computing method for all flow speeds [J]. Journal of Computational Physics, 1974, 14(3): 227–253. DOI: 10.1016/0021-9991(74)90051-5. [5] OSHER S, SETHIAN J A. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations [J]. Journal of Computational Physics, 1988, 79(1): 12–49. DOI: 10.1016/0021-9991(88)90002-2. [6] TRYGGVASON G, BUNNER B, ESMAEELI A, et al. A front-tracking method for the computations of multiphase flow [J]. Journal of Computational Physics, 2001, 169(2): 708–759. DOI: 10.1006/jcph.2001.6726. [7] FEDKIW R P, ASLAM T, MERRIMAN B, et al. A non-oscillatory Eulerian approach to interfaces in multimaterial flows: the ghost fluid method [J]. Journal of Computational Physics, 1999, 152(2): 457–492. DOI: 10.1006/jcph.1999.6236. [8] LIU T G, KHOO B C, WANG C W. The ghost fluid method for compressible gas-water simulation [J]. Journal of Computational Physics, 2005, 204(1): 193–221. DOI: 10.1016/j.jcp.2004.10.012. [9] SCHOCH S, NORDIN-BATES K, NIKIFORAKIS N. An Eulerian algorithm for coupled simulations of elastoplastic-solids and condensed-phase explosives [J]. Journal of Computational Physics, 2013, 252: 163–194. DOI: 10.1016/j.jcp.2013.06.020. [10] CROWL W K. Structures to resist the effects of accidental explosions [M]. USA: US Army, Navy and Air Force, US Government Printing Office, 1969: 205−315. [11] 徐维铮, 吴卫国. 爆炸波高精度数值计算程序开发及应用 [J]. 中国舰船研究, 2017, 12(3): 64–74. DOI: 10.3969/j.issn.1673-3185.2017.03.010.XU Weizheng, WU Weiguo. Development of in-house high-resolution hydrocode for assessment of blast waves and its application [J]. Chinese Journal of Ship Research, 2017, 12(3): 64–74. DOI: 10.3969/j.issn.1673-3185.2017.03.010. [12] TÜRKER L. Thermobaric and enhanced blast explosives (TBX and EBX) [J]. Defence Technology, 2016, 12(6): 423–445. DOI: 10.1016/j.dt.2016.09.002. [13] 张洪武, 何扬, 张昌权. 空中爆炸冲击波地面荷载的数值模拟 [J]. 爆炸与冲击, 1992, 12(2): 156–165.ZHANG Hongwu, HE Yang, ZHANG Changquan. Numerical simulation on ground surface loading of shock wave from air explosions [J]. Explosion and Shock Waves, 1992, 12(2): 156–165. [14] 赵海涛, 王成. 空中爆炸问题的高精度数值模拟研究 [J]. 兵工学报, 2013, 34(12): 1536–1546. DOI: 10.3969/j.issn.1000-1093.2013.12.008.ZHAO Haitao, WANG Cheng. High resolution numerical simulation of air explosion [J]. Acta Armamentarii, 2013, 34(12): 1536–1546. DOI: 10.3969/j.issn.1000-1093.2013.12.008. [15] BAKER W E. Explosions in air [M]. USA: University of Texas Press, 1973: 55−95. [16] 姚成宝, 李若, 田宙, 等. 空气自由场中强爆炸冲击波传播二维数值模拟 [J]. 爆炸与冲击, 2015, 35(4): 585–590. DOI: 10.11883/1001-1455(2015)04-0585-06.YAO Chengbao, LI Ruo, TIAN Zhou, et al. Two dimensional simulation for shock wave produced by strong explosion in free air [J]. Explosion and Shock Waves, 2015, 35(4): 585–590. DOI: 10.11883/1001-1455(2015)04-0585-06. [17] 姚成宝, 浦锡锋, 寿列枫, 等. 强爆炸冲击波在不均匀空气中传播数值模拟 [J]. 计算力学学报, 2015, 32(S1): 6–9.YAO Chengbao, PU Xifeng, SHOU Liefeng, et al. Numeircal simulation of blast wave propagation in nonuniform air [J]. Chinese Journal of Computational Mechanics, 2015, 32(S1): 6–9. [18] SYMBALISTY E M D, ZINN J, WHITAKER R W. RADFLO physics and algorithms: LA-12988-MS [R]. USA: Los Alamos National Lab, 1995. DOI: 10.2172/110714. [19] SETHIAN J A. Evolution, implementation, and application of level set and fast marching methods for advancing fronts [J]. Journal of Computational Physics, 2001, 169(2): 503–555. DOI: 10.1006/jcph.2000.6657. [20] SUSSMAN M, SMEREKA P, OSHER S. A level set approach for computing solutions to incompressible two-phase flow [J]. Journal of Computational Physics, 1994, 114(1): 146–159. DOI: 10.1006/jcph.1994.1155. [21] DI Yana, LI Ruo, TANG Tao, et al. Level set calculations for incompressible two-phase flows on a dynamically adaptive grid [J]. Journal of Scientific Computing, 2007, 31(1/2): 75–98. DOI: 1007/s10915-006-9119-3. [22] TORO E F. Riemann solvers and numerical methods for fluid dynamics [M]. Berlin, Heidelberg: Springer Berlin Heidelberg, 2009: 102-200. DOI: 10.1007/b79761. [23] LI R, WU S N. h-adaptive mesh method with double tolerance adaptive strategy for hyperbolic conservation laws [J]. Journal of Scientific Computing, 2013, 56(3): 616–636. DOI: 10.1007/s10915-013-9692-1. 期刊类型引用(10)
1. 陆曼君,解利军,祁佳晨,王攀,姚成宝. 面向爆炸压力场的体绘制自动传输函数设计. 计算机应用研究. 2024(02): 602-608 . 百度学术
2. 王振宁,尹建平,伊建亚,李旭东. 柱形装药近地动爆冲击波周向传播规律研究. 爆炸与冲击. 2023(06): 90-107 . 本站查看
3. 高永红,王子聪,张伟,辛凯,王巍,段亚鹏. 强冲击荷载作用下钢结构建筑毁伤后果研究. 结构工程师. 2023(03): 141-149 . 百度学术
4. 崔元博,孔德仁,张学辉,张逸飞. 典型炸药爆炸过程中电磁辐射特性分析. 国防科技大学学报. 2022(06): 70-80 . 百度学术
5. 姚成宝,付梅艳,韩峰,闫凯,雷雨. 基于多介质Riemann问题的流体-固体耦合数值方法及其在爆炸与冲击问题中的应用. 兵工学报. 2021(02): 340-355 . 百度学术
6. 张琳,张泽,姚李刚,杨剑波,秦将,潘登,王军龙. 防爆舱设计及抗爆炸冲击波超压试验研究. 科学技术与工程. 2021(29): 12383-12389 . 百度学术
7. 任会兰,储著鑫,栗建桥,马天宝. B炸药爆炸过程中电磁辐射研究. 力学学报. 2020(04): 1199-1210 . 百度学术
8. 杨科之,刘盛. 空气冲击波传播和衰减研究进展. 防护工程. 2020(03): 1-10 . 百度学术
9. 曹涛,孙浩,周游,罗赓,孙吉伟. 近地爆炸冲击波传播特性数值模拟与应用. 兵器装备工程学报. 2020(12): 187-191 . 百度学术
10. 王良全,商飞,孔德仁. 静动爆冲击波数值仿真分析. 兵器装备工程学报. 2020(12): 208-213 . 百度学术
其他类型引用(11)
-