不同风场下TNT 炸药爆炸烟云的扩散模型及特性*

段中山，龚朋彬，袁 伟，过惠平，罗永锋，罗昆升

（1. 火箭军工程大学核工程学院，陕西 西安 710025；2. 中国人民解放军陆军勤务学院环境工程教研室，重庆 401331；3. 火箭军研究院，北京 100094）

炸药爆炸烟云抬升高度变化与爆炸当量密切相关，建立烟云扩散高度变化模型可用于反演未知爆炸事故类型、当量或推测战场武器型号、威力及毁伤效应，实现未知爆炸源实时快速侦测。而对于脏弹、核武器炸药化学爆炸或其他“生化”爆炸，此类爆炸往往具有高度隐蔽性、烟雨危害性和应急时效性，事故初期难以通过调查现场炸坑痕迹、玻璃受冲击痕迹和爆炸震动记录等传统反演方法获取其爆炸当量及烟云高度[1]。由于移动摄像、视频监控、地面遥感、卫星观测等影像获取手段的使用与普及，使得利用爆炸视频中烟云高度变化信息快速反演炸药当量、预测污染烟云最终高度成为可能。总之，研究烟云扩散高度变化规律对战场爆炸远程评测、爆炸事故非现场侦测、有害物质爆轰污染快速预警与评估等都具有重要价值。

为获取复杂风场下爆炸烟云扩散高度变化过程及最终高度，本文拟采用仿真与实验相结合的方法开展研究。首先采用CFD 方法计算烟云时空分布模型；然后利用烟云时空分布实验结果验证建模方法和参数设置的正确性，确定一套可表征烟云扩散时空分布的建模与参数设置方法；最后建立不同风速下烟云扩散高度变化模型和最终高度计算方法并讨论风对烟云高度影响过程与机理。

1 爆炸烟云扩散理论

1.1 爆炸烟云扩散物理过程

炸药爆炸烟云扩散过程可分为4 个阶段[4]。阶段1：爆轰阶段，炸药引爆后爆轰产物不断膨胀直至地面烟团直径达到最大值，炸药参数、装置结构和地面环境决定了初始烟团形态、尺寸、密度、温度及垂直动量等参数，该过程几乎不受大气条件影响。阶段2：热抬升阶段，地面烟团在大气浮力作用下将脱离地面开始热抬升，此阶段烟云平均温度高于环境温度，大多存在化学燃烧和地面反冲冲量，持续时间几秒到十几秒。阶段3：持续上升阶段，浮力和惯性作用将烟云以扩散与湍流的形式抬升至稳定状态，持续时间达数十秒至上百秒。阶段2 和阶段3 可统称爆炸烟云主动抬升过程，过程中水平风和大气稳定度影响烟云与大气热交换速度、空气混入速度等，进而改变烟云高度分布[11]。阶段4：被动扩散阶段，烟云稳定后几乎不再上升，随着大气风场被动扩散并辐射下游区域。

1.2 爆炸烟云扩散理论方程

爆炸烟云主动抬升过程中主要受重力、上升阻力和大气密度差带来的浮力等影响，其扩散过程可以用N-S（奈维-斯托克斯）方程描述[11]：

式中：ρ 为流体密度，烟云简化为不可压缩时可认为ρ 为常数；u 为流体速度；p 为压强；µ为流体黏性系数；f 为外力项。

式(1)为流体质量守恒方程；式(2)为流体动量守恒方程，右边四项分别为外力项、压强项、黏性项、对流项。

式中：a、b 为系数项，T 为烟云温度，Tatm为环境温度，y 垂直向上向量。

烟云密度场和温度场受到烟云速度场的传送作用，可用如下关系表示：

2）其他土地利用类型都有不同程度的减少，其中水域和耕地减少的最为明显，1985－2000年和2000－2016年水域的减少速度分别为0.08%和0.19%，耕地的减少速度分别为0.12%和0.14%；

烟云扩散过程中存在热量交换，遵循能量守恒传递方程：

式中：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 实验用于验证无风与典型风场下的烟云扩散时空分布参数及模型。

选择阴天傍晚时刻进行实验以保证大气处于稳定状态，对TNT 裸装药进行地面引爆。实验布局如图1 所示，主要过程可分为4 步：(1) 设计1、16、62 kg 的TNT 炸药化学爆炸装置，为尽可能长时间地记录烟云，在实验爆炸装置外围均匀布置少量烟云示踪剂；(2) 对于1 kg TNT 爆炸，通过26 m 吊塔钢丝线直接标定烟云高度，对于16 和62 kg TNT 爆炸，则通过释放绳索牵引的高空气球后，采用激光测距仪测量气球高度，结合角度测量仪标定高空视野；(3) 根据测风仪测风结果选择风速到达实验要求时刻起爆装置，同时记录爆炸烟团初始形态及时空分布；(4) 根据视场标定结果，利用专业像素分析软件对实验烟云图片信息开展计算，获取爆炸烟云扩散时空分布数据。

表 1 实验分组表Table 1 Experiment group table

图 1 实验场地布局Fig. 1 Experimental site layout

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 左右烟云基本稳定后几乎不再上升。计算结果则显示由于烟云密度梯度差产生的浮力将烟团由高斯形态浓度分布转变成由两个反漩涡环流组成的涡流分布，横截面上展现出由两个反漩涡环流所组成的“肾形状”模式，也就是常见的“蘑菇云”现象。实验和仿真结果可见，实验烟云扩散过程中具体形态不及仿真烟云对称和规整，但两者的形态演变过程和时空分布参数基本一致，验证了爆炸烟云低密度浮力烟团扩散机理和双涡环反漩涡环流扩散方式的正确性。

图 2 实验与仿真烟云时空形态对比Fig. 2 Comparison of time and space patterns between experimental and simulated clouds

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 实验与仿真烟云高度对比Fig. 3 Comparison of heights between experimental and simulated clouds

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]。

图 4 16 kg TNT 爆炸的烟云高度变化Fig. 4 Change of cloud heights for the explosion of 16 kg TNT

图 5 62 TNT 爆炸的烟云高度变化Fig. 5 Change of cloud heights for the explosion of 62 kg TNT

图 6 高度变化参数(b)拟合直线Fig. 6 Fitting line of heights variation parameter b

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 公式无法准确计算爆高，而该模型则可用于不同风速下爆高计算。复杂扩散条件下烟云爆高计算模型的建立可使爆高计算误差更小，进而提高烟云污染评估结果准确度。

图 7 16 kg TNT 爆炸的烟云最终高度Fig. 7 Final cloud height for 16 kg TNT explosion

图 8 62 kg TNT 烟云最终高度Fig. 8 Final cloud height for 62 kg TNT explosion

3.2.3 最终高度计算模型使用误差分析

烟云最终高度模型误差主要有三个来源：实验估计烟云最终高度误差，线性拟合R 值过大和烟云扩散计算模型适应条件。单纯实验估计最终测量高度定量误差在10%～15%[4]，由于仿真计算能清晰显示烟云顶高边界，极大降低了这一误差。模型线性拟合精度通过拟合参数R 评估，高度变化模型拟合度R＞0.99，最终高度模型拟合度R＞0.92，最终高度模型误差在10%内，该精度用于烟云污染评估可接受。最终高度模型忽略了较小地面因素及爆轰差异影响，计算模型适应条件为稳定大气，相比不稳定大气下模型烟云高度整体偏低[4]，该模型计算结果可认为是误差较小的爆炸烟云扩散高度下限值，烟云越低其计算所得危害结果却代表近场最大危害结果[13]，在烟云污染应急评估中采用高度下限值划分的应急区域是面积较大的保守结果，符合事故初期源项不确定情况下安全评价与应急救援基本原则。

图 9 最终高度参数d 的拟合直线Fig. 9 Fitting line of final height parameter (d)

3.3 不同风速下烟云扩散参数变化规律

针对前期不同风速下出现的高度变化差异，图10 和11 提取了烟云扩散参数云图中的最高温度、最小密度和最大速度进行数值分析。烟云中最高温度出现在左右涡环中心处，数据显示温度随时间呈反比例衰减，热释放过程大多在5～10 s，当量越大下降至常温的速度越慢，风速越高下降至常温越快。烟云密度分布呈现从涡环中心到外环逐步增大的规律，实验中也观测到烟云浓度外围大于中心处现象，烟云最小密度出现在涡环中心处，变化呈现前期快后期慢的指数函数增长规律，风速越大密度越快趋近空气密度，无风条件下16 和62 kg 烟云最小密度增至1 kg/m3需2 和30 s，而风场条件下仅10 和12 s，这是导致烟云高度差距的主要原因。烟云上升过程中最大速度并非上升速度而是涡环中心处的漩涡速度，3 秒左右烟云涡环翻滚至最大速度之后逐步衰减，风速越大衰减速度越快，因其先加速后衰减的特性导致烟云前期尽管与空气快速混合但高度差距较小，中后期速度和密度已接近空气导致其提前进入稳定状态，烟云最终高度也因此出现较大差距。

图 10 16 kg TNT 爆炸烟云扩散参数变化Fig. 10 Change of cloud diffusion parameters for 16 kg TNT explosions

图 11 62 kg TNT 爆炸烟云扩散参数变化Fig. 11 Change of cloud diffusion parameters for 62 kg TNT explosions

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) 爆炸烟云仿真结果显示烟云流场由左右两个反漩涡环流组成，温度、速度最大值和密度最小值均出现在左右涡环中心处。实验和仿真结果均显示风将加快烟云与空气混合速度，风速越大烟云温度和速度衰减越快，密度加速趋于空气。风场下烟云高度在前期热抬升阶段差距较小，但风对烟云抬升中后期瓦解作用明显导致高度差距明显，风速越大烟云抬升至稳定的时间越短、烟云最终高度越低。