Experimental and simulation of diffusion model and characteristics of explosive cloud at different wind fields
-
摘要: 为获取不同风场下TNT爆炸烟云扩散时空分布规律与高度变化模型,本文理论描述了爆炸烟云扩散过程与机理,开展了不同水平风速下烟云扩散的计算流体力学(computational fluid dynamics, CFD)仿真和外场时空分布实验,建立了不同水平风速下烟云高度随时间变化模型及烟云最终高度计算模型,分析了烟云扩散过程中形态、温度、密度、速度变化规律。研究结果显示:CFD方法仿真烟云分布结果与实验结果基本一致,大气稳定且无风条件下烟云高度随时间呈指数0.5的幂函数关系,最终高度与爆炸当量可拟合为指数0.47的幂函数模型;水平风会加快烟云与空气混合的速度,导致幂函数模型中指数参数随风速变大而呈线性减小规律,风速越大烟云上升速度衰减越快、上升时间越短、最终高度越低。Abstract: In order to obtain the spatial and temporal distribution law and height variation model of TNT explosion cloud diffusion under different wind Fields, this paper theoretically describes the diffusion process and mechanism of the explosion cloud, and carries out the computational fluid dynamics (CFD) simulation and the field-time distribution experiment of the cloud diffusion under different horizontal wind speeds, and analyzes the diffusion process of the cloud. Morphology, temperature, density and speed change law were established, and the variation model of cloud height with time at different u-wind speeds and the final height calculation model of cloud were established. The results show that the CFD method simulation cloud diffusion results are consistent with the experimental results. The cloud height has a power function with an exponent of 0.5 under windless conditions. The final height and explosive equivalent can be fitted to a power function model with an index of 0.47. The horizontal wind will speed up the mixing of the cloud and the air, causing the exponential parameter of the power function model to decrease linearly with the wind speed becoming larger. The higher the wind speed, the faster the decay rate of the cloud, the shorter the rise time, and the lower the final height.
-
炸药爆炸烟云抬升高度变化与爆炸当量密切相关,建立烟云扩散高度变化模型可用于反演未知爆炸事故类型、当量或推测战场武器型号、威力及毁伤效应,实现未知爆炸源实时快速侦测。而对于脏弹、核武器炸药化学爆炸或其他“生化”爆炸,此类爆炸往往具有高度隐蔽性、烟雨危害性和应急时效性,事故初期难以通过调查现场炸坑痕迹、玻璃受冲击痕迹和爆炸震动记录等传统反演方法获取其爆炸当量及烟云高度[1]。由于移动摄像、视频监控、地面遥感、卫星观测等影像获取手段的使用与普及,使得利用爆炸视频中烟云高度变化信息快速反演炸药当量、预测污染烟云最终高度成为可能。总之,研究烟云扩散高度变化规律对战场爆炸远程评测、爆炸事故非现场侦测、有害物质爆轰污染快速预警与评估等都具有重要价值。
爆炸烟云抬升主要受爆轰条件及扩散条件影响。目前研究主要集中于获取爆炸烟云最终高度用于污染评估,利用烟云信息反演当量或推测烟云高度的研究还鲜为报导。Church[2]根据53 kg以上TNT爆炸及核爆炸烟云高度拟合得到爆高公式(H=76M0.25),该公式无法描述烟云高度变化规律,且计算小当量炸药时,爆高误差较大。Thielen等[3]和Sharon等[4]根据不同当量及大气条件下炸药爆炸实验结果拟合了小于100 kgTNT烟云爆高模型,由于实验条件控制较差,实验拟合精度低且重复性较差。单纯依靠实验方法开展烟云研究成本高且风速、大气稳定度等随机因素影响较大,烟云后期由于高空视野消失更是难以准确测量,要获取复杂条件下烟云扩散定量模型需结合可表征烟云时空分布的数值仿真方法。Makhviladze等[5]、Kansa[6]和Kanarska等[7]采用多相流理论计算了爆炸烟云涡环运动,郑毅[8]和李晓丽等[9]利用变密度法开展了中小尺度烟云时空分布数值仿真,以上均未讨论风对烟云扩散高度的影响,Zhang等[10]通过理论计算证明了风对爆炸烟云高度有重要影响,未对其作用机理和定量模型开展相关研究。
为获取复杂风场下爆炸烟云扩散高度变化过程及最终高度,本文拟采用仿真与实验相结合的方法开展研究。首先采用CFD方法计算烟云时空分布模型;然后利用烟云时空分布实验结果验证建模方法和参数设置的正确性,确定一套可表征烟云扩散时空分布的建模与参数设置方法;最后建立不同风速下烟云扩散高度变化模型和最终高度计算方法并讨论风对烟云高度影响过程与机理。
1. 爆炸烟云扩散理论
1.1 爆炸烟云扩散物理过程
炸药爆炸烟云扩散过程可分为4个阶段[4]。阶段1:爆轰阶段,炸药引爆后爆轰产物不断膨胀直至地面烟团直径达到最大值,炸药参数、装置结构和地面环境决定了初始烟团形态、尺寸、密度、温度及垂直动量等参数,该过程几乎不受大气条件影响。阶段2:热抬升阶段,地面烟团在大气浮力作用下将脱离地面开始热抬升,此阶段烟云平均温度高于环境温度,大多存在化学燃烧和地面反冲冲量,持续时间几秒到十几秒。阶段3:持续上升阶段,浮力和惯性作用将烟云以扩散与湍流的形式抬升至稳定状态,持续时间达数十秒至上百秒。阶段2和阶段3可统称爆炸烟云主动抬升过程,过程中水平风和大气稳定度影响烟云与大气热交换速度、空气混入速度等,进而改变烟云高度分布[11]。阶段4:被动扩散阶段,烟云稳定后几乎不再上升,随着大气风场被动扩散并辐射下游区域。
1.2 爆炸烟云扩散理论方程
爆炸烟云主动抬升过程中主要受重力、上升阻力和大气密度差带来的浮力等影响,其扩散过程可以用N-S(奈维-斯托克斯)方程描述[11]:
∂ρ∂t+∇⋅(ρu)=0 (1) ∂u∂t=f−(u⋅∇)u−1ρ∇p+μ∇2u (2) 式中:ρ为流体密度,烟云简化为不可压缩时可认为ρ为常数;u为流体速度;p为压强;μ为流体黏性系数;f为外力项。
式(1)为流体质量守恒方程;式(2)为流体动量守恒方程,右边四项分别为外力项、压强项、黏性项、对流项。
f=−aρy+b(T−Tatm)y (3) 式中:a、b为系数项,T为烟云温度,Tatm为环境温度,y垂直向上向量。
烟云密度场和温度场受到烟云速度场的传送作用,可用如下关系表示:
{∂ρ∂t=−u⋅∇ρ∂T∂t=−u⋅∇T (4) 烟云扩散过程中存在热量交换,遵循能量守恒传递方程:
∂∂t(ρE)+∂∂x[u(ρE+p)]=∂∂x(keff∂T∂x)−∑j′hj′Jj′+uj(τj)eff+Sh (5) 式中:E为流体微团总能量;
hj′ 、Jj′ 分别表示组分j′ 的焓和扩散通量;keff为有效热传导系数;Sh为热源项,包括了化学反应热及其他体积热源项。式(5)右边前3项分别为热传导、组分扩散和黏性耗散能量输运。烟云扩散N-S方程待求解变量为速度u、密度ρ、温度T和压强p。CFD计算方法一般通过划分待求解区域网格来求解N-S微分方程组,计算后可得各网格点上各时刻物理变量的分布场。
2. 爆炸烟云扩散实验与仿真
2.1 烟云扩散实验设计与实施
实验选取1、16、62 kg TNT先后开展5组外场烟云扩散实验研究。实验分组如表1所示,其中1 kg TNT无风实验用于近距离观察烟云扩散形态演化过程及其特征,16和62 kg实验用于验证无风与典型风场下的烟云扩散时空分布参数及模型。
表 1 实验分组表Table 1. Experiment group table分组 M/kg v/(m∙s−1) 分组 M/kg v/(m∙s−1) 1 1 <0.5 (无风) 4 16 3±0.5(典型风) 2 16 <0.5 (无风) 5 62 3±0.5(典型风) 3 62 <0.5 (无风) 注:M为TNT当量,v为风速。 选择阴天傍晚时刻进行实验以保证大气处于稳定状态,对TNT裸装药进行地面引爆。实验布局如图1所示,主要过程可分为4步:(1) 设计1、16、62 kg的TNT炸药化学爆炸装置,为尽可能长时间地记录烟云,在实验爆炸装置外围均匀布置少量烟云示踪剂;(2) 对于1 kg TNT爆炸,通过26 m吊塔钢丝线直接标定烟云高度,对于16和62 kg TNT爆炸,则通过释放绳索牵引的高空气球后,采用激光测距仪测量气球高度,结合角度测量仪标定高空视野;(3) 根据测风仪测风结果选择风速到达实验要求时刻起爆装置,同时记录爆炸烟团初始形态及时空分布;(4) 根据视场标定结果,利用专业像素分析软件对实验烟云图片信息开展计算,获取爆炸烟云扩散时空分布数据。
2.2 烟云扩散仿真思路与方法
烟云抬升演化过程主要指阶段2和3,阶段1可采用实验经验公式结合有限元软件AUTODYN计算共同判定地面膨胀最大时刻非均匀的初始烟云形态、尺寸、密度、温度及垂直动量,将初始烟云参数优化后耦合到可编程软件FLUENT中构建烟云扩散仿真模型[12-14]。
初始烟云参数:地面烟团最大膨胀半径符合实验经验公式r=1.93M0.32/ (T1/3 600)1/3,T1为爆温(2 861 K)[9];AUTODYN爆轰计算模型显示几十毫秒内初始烟云达到最大膨胀体积,由于地面影响初始烟团可简化为半椭球形态耦合到CFD模型中,1、16、62 kg TNT的地面烟团直径约为4.17、10.12、15.61 m,高约3.13,、7.59、11.71 m,火球膨胀后压力与空气压力相当,整体符合ρ1=ρ0(T0/ T1)[9],其中ρ1为烟团密度,ρ0与T0为空气密度与温度,ρ1均值取0.12 kg/m3,由于地面反冲效应导致烟云存在约5 m/s初始垂直动量,由于温度迅速下降取整体烟团均值温度约1 000 K,烟团组分按TNT爆炸产物体积分数设定,假设初始烟团参数分布均匀且忽略颗粒携带、炸药抛洒和后续化学燃烧[15-16]。
边界与网格:在二维空间中采用从中间到两边、从下到上递增的步进网格划分方式,这样中间和下部网格密集利于关注烟云扩散;计算中发现边界条件设置对烟云分布有影响,综合分析与试算后将文献[12]中边界设置变为左边界为Inflow、上边界与右边界均为Pressure Outlet、底部边界Wall计算更科学;为使风场更真实,采用UDF编程风速随高度变化规律来设置风场,风速u=U10(z/10)a,U10为地面10 m处风速,大小取0~6 m/s(0~4级),z为高度,a为地面粗糙系数(取0.3),在空气入口处加载UDF velocity[17],5级及以上风速属较极端气象条件暂不做讨论。
算法与参数:采用混合物耦合计算、湍流标准k-ε模型、随时间变化的非定常流计算[17-18];计算中对x向和y向速度、能量、各物相体积分数及k值等进行收敛设定和残差监视。
3. 结果分析
3.1 烟云扩散仿真与实验结果分析
3.1.1 烟云形态时空分布结果分析
图2为1 kg TNT爆炸后不同时刻烟云分布实验与仿真结果,其中实验吊塔高度26 m,仿真空间高度为40 m。实验较完整记录了烟云时空分布及形态变化过程,前期烟云由于地面反冲动量及浮力作用加速上升、宽高比例较小,后期在浮力涡环主导下烟云扩散整体呈现“球形态”。烟团上端边界分明,下端烟柱明显,由于大气不确定性、示踪剂非均匀抛洒等影响导致后期烟云左侧出现少许不规则形态,30 s左右烟云基本稳定后几乎不再上升。计算结果则显示由于烟云密度梯度差产生的浮力将烟团由高斯形态浓度分布转变成由两个反漩涡环流组成的涡流分布,横截面上展现出由两个反漩涡环流所组成的“肾形状”模式,也就是常见的“蘑菇云”现象。实验和仿真结果可见,实验烟云扩散过程中具体形态不及仿真烟云对称和规整,但两者的形态演变过程和时空分布参数基本一致,验证了爆炸烟云低密度浮力烟团扩散机理和双涡环反漩涡环流扩散方式的正确性。
3.1.2 烟云扩散高度分布实验与仿真结果
实验烟云早期能清楚显示上边界,由于空气稀释中后期烟云上边界模糊,特别是接近顶高时更是无法表征,除1 kg TNT能测量实验烟云最终高度外,16和62 kg烟云后期均难以准确测量,只能通过仿真表征其高度。图3结果显示仿真与实验烟云高度分布基本一致,实验结果由于风场不确定性和视野误差出现少量拐点,仿真曲线平滑性和规律性更好。实验烟云前期(爆后10 s)高度略高于其仿真值,考虑为实验烟云扩散前期存在地面反冲和燃烧反应,中后期实验与仿真规律几乎一致,平均相对误差约5%说明仿真结果可信度高。实验和仿真均显示,典型风场条件下烟云高度会降低,风对烟云第2阶段影响较小,16 kg TNT爆炸前15 s和62 kg TNT爆炸前25 s无风和典型风下烟云高度差距15%内,考虑为前期快速加速和逐步减速过程中,尽管风加速烟云与空气混合,但加速和减速过程总路程(高度)差距较小。风场加速空气与烟云混合导致第3阶段烟云高度差距扩大,烟云快速瓦解、后期上升动力不足,典型风场烟云上升时间和最终高度明显小于无风烟云,风场烟云将较快趋于稳定。
3.2 不同风速下烟云扩散高度模型
3.2.1 烟云扩散过程中高度变化模型
不同风速下的烟云高度分布数据如图4~5所示,可以看出,不同风速下烟云高度变化呈前快后慢的幂函数增长规律,遵循浮力加速与重力、大气阻力减速物理模型。采用Matlab幂函数数学模型H(t)=atb进行高度变化规律拟合(H单位为m,t单位为s),同公斤级下初始烟云H(1)高度相同[12],a=H(1)=(6.3±1)M(0.29±0.03),计算所得b值如图6所示,发现b值与风速呈线性关系,拟合后得到b=(0.5±0.003)−(0.024±0.001)v,拟合度R为99%,故烟云高度随时间变化模型H(t)=(6.3±1)M(0.29±0.03)t(0.5±0.003)−(0.024±0.001)v。该扩散规律可用于非现场方式反演爆炸当量并推测污染烟云高度,直接实现爆炸危害快速侦测与预警,也可结合其他现场痕迹计算方法精确反演爆炸当量[1]。
3.2.2 烟云扩散最终高度计算模型
不同风速下的爆炸烟云最终高度数据如图7~8所示,可以看出最终高度与风速大小呈线性下降趋势,高度与当量整体符合Hmax=cMd模型(Hmax单位为m,M单位为kg)[12],由实验可得当M=1时,c值区间为33~36 m,高度模型可转变为Hmax(M,v)=(34.5±1.5)Md(v),取33、34.5、36分别代入16和62 kg最大高度下求解d值,所得d值如图9所示,可以看出d值随风速呈线性下降,拟合得到d(v)=(0.47±0.01)−(0.038±0.002)v,拟合度R为92.6%,Hmax(M,v)=(34.5±1.5)M(0.47±0.01)−(0.038±0.002)v。该爆高计算模型与Church公式相比,当v=0时该模型在计算小于50 kg当量爆高时准确度更高,在50 kg以上与Church计算结果差距不大;当v≠0时,Church公式无法准确计算爆高,而该模型则可用于不同风速下爆高计算。复杂扩散条件下烟云爆高计算模型的建立可使爆高计算误差更小,进而提高烟云污染评估结果准确度。
3.2.3 最终高度计算模型使用误差分析
烟云最终高度模型误差主要有三个来源:实验估计烟云最终高度误差,线性拟合R值过大和烟云扩散计算模型适应条件。单纯实验估计最终测量高度定量误差在10%~15%[4],由于仿真计算能清晰显示烟云顶高边界,极大降低了这一误差。模型线性拟合精度通过拟合参数R评估,高度变化模型拟合度R>0.99,最终高度模型拟合度R>0.92,最终高度模型误差在10%内,该精度用于烟云污染评估可接受。最终高度模型忽略了较小地面因素及爆轰差异影响,计算模型适应条件为稳定大气,相比不稳定大气下模型烟云高度整体偏低[4],该模型计算结果可认为是误差较小的爆炸烟云扩散高度下限值,烟云越低其计算所得危害结果却代表近场最大危害结果[13],在烟云污染应急评估中采用高度下限值划分的应急区域是面积较大的保守结果,符合事故初期源项不确定情况下安全评价与应急救援基本原则。
3.3 不同风速下烟云扩散参数变化规律
针对前期不同风速下出现的高度变化差异,图10和11提取了烟云扩散参数云图中的最高温度、最小密度和最大速度进行数值分析。烟云中最高温度出现在左右涡环中心处,数据显示温度随时间呈反比例衰减,热释放过程大多在5~10 s,当量越大下降至常温的速度越慢,风速越高下降至常温越快。烟云密度分布呈现从涡环中心到外环逐步增大的规律,实验中也观测到烟云浓度外围大于中心处现象,烟云最小密度出现在涡环中心处,变化呈现前期快后期慢的指数函数增长规律,风速越大密度越快趋近空气密度,无风条件下16和62 kg烟云最小密度增至1 kg/m3需2和30 s,而风场条件下仅10和12 s,这是导致烟云高度差距的主要原因。烟云上升过程中最大速度并非上升速度而是涡环中心处的漩涡速度,3秒左右烟云涡环翻滚至最大速度之后逐步衰减,风速越大衰减速度越快,因其先加速后衰减的特性导致烟云前期尽管与空气快速混合但高度差距较小,中后期速度和密度已接近空气导致其提前进入稳定状态,烟云最终高度也因此出现较大差距。
4. 结 论
(1) 爆炸烟云仿真与实验结果显示烟云时空分布过程与参数几乎一致,说明基于CFD建模的烟云扩散仿真方法可较好表征烟云扩散过程、形态及物理参数。
(2) 大气稳定不同风速条件下烟云可见阶段高度分布拟合函数为H(t)=(6.3±1.0)M(0.29±0.03)t(0.5±0.003)−(0.024±0.001)v,该规律可初步用于实现非现场方式快速反演爆炸当量并推测烟云污染最终高度。风场下烟云扩散最终高度理论模型为Hmax(M,v)=(34.5±1.5)M(0.47±0.01)−(0.038±0.002)v,该模型可较准确计算复杂风场条件下事故污染烟云稳定高度值。
(3) 爆炸烟云仿真结果显示烟云流场由左右两个反漩涡环流组成,温度、速度最大值和密度最小值均出现在左右涡环中心处。实验和仿真结果均显示风将加快烟云与空气混合速度,风速越大烟云温度和速度衰减越快,密度加速趋于空气。风场下烟云高度在前期热抬升阶段差距较小,但风对烟云抬升中后期瓦解作用明显导致高度差距明显,风速越大烟云抬升至稳定的时间越短、烟云最终高度越低。
-
表 1 实验分组表
Table 1. Experiment group table
分组 M/kg v/(m∙s−1) 分组 M/kg v/(m∙s−1) 1 1 <0.5 (无风) 4 16 3±0.5(典型风) 2 16 <0.5 (无风) 5 62 3±0.5(典型风) 3 62 <0.5 (无风) 注:M为TNT当量,v为风速。 -
[1] 周宣赤. 基于爆炸现场痕迹反演爆源参数方法及应用[D]. 北京: 北京理工大学, 2014: 3−14. [2] CHURCH H W. Cloud rise from high explosive detonations: TID-4500 [R]. USA: Sandia National Laboratories (SNL), 1969. [3] THIELEN H, SCHRODL E. Blast experiments for the derivation of initial cloud dimensions after a “Dirty Bomb” event: Koln- 50667 [R]. Germany, 2004. [4] SHARON A, HALEVY I, SATTINGER D, et al. Cloud rise model for radiological dispersal devices events [J]. Atmospheric Environment, 2012, 54: 603–610. DOI: 10.1016/j.atmosenv.2012.02.050. [5] MAKHVILADZE G M, ROBERTS J P, YAKUSH S E. Modelling of atmospheric pollution by explosions [J]. Environmental Software, 1995, 10(2): 117–127. DOI: 10.1016/0266-9838(94)00005-R. [6] KANSA E J. A time dependant buoyant puff model for explosion sources: UCRL-ID-128733 [R]. USA: Lawrence Livermore National Laboratory, 1997. [7] KANARSKA Y, LOMOV I, GLENN L, et al. Numerical simulation of cloud rise phenomena associated with nuclear bursts [J]. Annals of Nuclear Energy, 2009, 36(10): 1475–1483. DOI: 10.1016/j.anucene.2009.08.009. [8] 郑毅. 瞬时热源(爆炸烟云)浮力涡环研究[D]. 北京: 清华大学, 2008: 12−40. [9] 李晓丽, 郑毅, 刘伟, 等. 爆炸烟云运动的试验与数值模拟研究初探 [J]. 核电子学与探测技术, 2011, 31(2): 131–135. DOI: 10.3969/j.issn.0258-0934.2011.02.002.LI X L, ZHENG Y, LIU W, et al. Numerical modeling and experimental research on the movement of the explosion clouds [J]. Nuclear Electronics and Detection Technology, 2011, 31(2): 131–135. DOI: 10.3969/j.issn.0258-0934.2011.02.002. [10] ZHANG X, MOUSSA N A, GROSZMANN D E, et al. A simple analytical buoyant puff rise model [C] // Proceedings of JANNAF Safety and Environmental Protection Subcommittee Meeting. Tampa, Florida, 1995. [11] BROWN R C, KOLB C E, CONANT J A, et al. Source characterization model (SCM) [R]. Washington: United States Department of Energy, 2004. [12] 段中山, 过惠平, 冯孝杰, 等. 爆炸烟云扩散的时空分布模型及特性 [J]. 爆炸与冲击, 2019, 39(5): 054202. DOI: 10.11883/bzycj-2017-0380.DUAN Z S, GUO H P, FENG X J, et al. Temporal and spatial distribution models of explosive cloud diffusion and their characteristics [J]. Explosion and Shock Waves, 2019, 39(5): 054202. DOI: 10.11883/bzycj-2017-0380. [13] LEBEL L, BOURGOUIN P, CHOUHAN S, et al. The sensitivity of atmospheric dispersion calculations in near-field applications: modeling of the FullScale RDD experiments with operational models in Canada, part I [J]. Health Physics, 2016, 110(5): 499–517. DOI: 10.1097/HP.0000000000000365. [14] KWAK H Y, KANG K M, KO I, et al. Fire-ball expansion and subsequent shock wave propagation from explosives detonation [J]. International Journal of Thermal Sciences, 2012, 59: 9–16. DOI: 10.1016/j.ijthermalsci.2012.04.022. [15] ABDUL-KARIM N, BLACKMAN C S, GILL P P, et al. The spatial distribution patterns of condensed phase post-blast explosive residues formed during detonation [J]. Journal of Hazardous Materials, 2016, 316: 204–213. DOI: 10.1016/j.jhazmat.2016.04.081. [16] FEDINA E, FUREBY C. Investigating ground effects on mixing and afterburning during a TNT explosion [J]. Shock Waves, 2013, 23(3): 251–261. DOI: 10.1007/s00193-012-0420-9. [17] 曹毅. 基于Fluent的大气雾霾流场研究[D]. 山东: 青岛科技大学, 2018: 31−43. [18] MISHRA K B, WEHRSTEDT K D, KREBS H. Boiling liquid expanding vapour explosion (BLEVE) of peroxy-fuels: experiments and computational fluid dynamics (CFD) simulation [J]. Energy Procedia, 2015, 66: 149–152. DOI: 10.1016/j.egypro.2015.02.082. -