发烟罐烟雾质量浓度时空分布特性

2020-11-24杨尚贤陈慧敏高丽娟马超齐斌邓甲昊

兵工学报 2020年9期

杨尚贤，陈慧敏，高丽娟，马超，齐斌，邓甲昊

(北京理工大学机电动态控制重点实验室,北京 100081)

0 引言

烟雾能够在一定时间和空间范围内对可见光、激光、毫米波等形成有效干扰和衰减。同时，人工烟雾释放操作简单、成本低廉，在现代战争领域得到了迅猛发展与广泛应用。研究烟雾的释放、扩散规律和浓度分布状态，为获取激光在烟雾环境下的传输特性，提高对地激光类武器系统的抗烟雾干扰性能，具有重要意义[1-3]。

发烟罐的作用原理是罐中发烟剂(以磷基发烟剂为代表)与氧气作用而燃烧，生成具有高分散度的氧化物粒子并扩散形成烟雾。烟雾较长时间悬浮在大气中，实现遮蔽、迷盲或干扰制导的目的[4]。为了得到烟雾的扩散规律、探索烟雾质量浓度的时空分布特性，近年来学者们做了大量研究。文献[5]采用数值模拟和野外测试相结合的方法，将烟幕初始云团的爆炸分散过程分为内部压力衰减为0之前和之后两个阶段，得到两个阶段烟幕初始云团的最大半径。文献[6]应用计算流体力学仿真软件Fluent研究了风速和风向对爆炸式催泪弹烟雾扩散的影响。文献[7]采用瞬时点源模拟爆炸烟雾，基于高斯烟团模式进行建模仿真，得到了烟雾浓度分布规律。文献[8]基于解析几何分析构建了垂直烟幕空间尺度计算模型，并基于高斯烟羽模式，应用牛顿下山法构建了单点释放垂直烟幕计算模型。文献[9]分析了发烟弹扩爆和烟雾在大气中扩散两个过程，提出了基于发烟弹扩爆后相关参数及气象条件来估算烟雾浓度的方法。文献[10]以湍流大涡算法对设置挡烟垂壁的防护工程烟雾扩散情况进行了仿真研究。文献[11]对地下娱乐建筑发生火灾时的烟雾流动进行模拟，分析了烟雾流动状态并与机械排风进行了比较。文献[12]基于气体- 固体两相流喷射的机理，数值模拟了微重力、高真空环境下的空间烟雾扩散。

以上文献应用不同方法针对烟雾扩散过程进行了研究，但仍存在以下问题：1)没有对烟雾的扩散进行分阶段、分时刻的详细研究，对烟雾扩散时烟团形状的变化过程没有展开描述；2)基于经验公式构建的烟雾扩散模型，以及基于烟雾释放参数和气象条件估计的烟雾空间质量浓度分布，数据误差较大；3)建筑通道对烟雾的扩散有阻碍作用，与发烟罐烟雾在自由空间中的扩散有明显差异。另外，大部分文献没有分析不同环境因素对烟雾扩散的影响规律。

为得到不同因素下烟雾扩散的时空分布特性，本文基于计算流体力学，应用Fluent软件对一种小型发烟罐烟雾的扩散过程进行仿真模拟。该方法能够直观地呈现计算结果，从而可以对烟雾的流动结构进行细致的分析。通过分析扩散时间、风速、质量流率等因素对烟雾形状、烟雾质量浓度分布状态的影响规律，为人工烟雾释放提供理论指导，并为进一步研究激光在烟雾中的传输特性奠定基础。

1 烟雾扩散数值模型

烟雾扩散过程属于典型的离散相扩散模型，建模时要同时考虑气流流动和烟雾颗粒的扩散运动。

1.1 气流流动数值模型

烟雾在上升过程中受到迎面空气的作用(包括空气介质的阻力、能量的交换)与上升流本身的扩散作用。设定流场速度u在3个坐标轴方向的分量为ux、uy、uz，温度场参数为T，流场压力为p.烟雾扩散过程遵循质量守恒定律、能量守恒定律和动量守恒定律。另外，烟雾的流动状态属于湍流流动，工程中常用雷诺时均方程来描述湍流状态，基本思想是通过k-ε(k为湍流动能；ε为耗散率)双方程模型表示时均方程中流体的瞬时脉动量[13]。相关的数学模型有如下5个。

质量守恒方程：

(1)

式中：ρ为空气密度(kg/m3)；xi为湍流模型的张量表示形式，i=x，y，z；ui(i=x，y，z)为流场速度在x轴、y轴、z轴方向上的分量；源项Sm为持续流入的气体质量(kg)。

能量守恒方程：

(2)

式中：cp为空气的定压比热容(J/(kg·℃))；λ为气体的导热系数(W/(m·℃))；源项ST为持续流入的气体热量(J)。

动量守恒方程：

(3)

式中：uj为流场速度在x轴、y轴、z轴3个方向的分量，j=x，y，z；xj为湍流模型的张量表示形式；μ为层流黏度系数(Pa·s)；μt为湍流黏度系数(Pa·s)，μt=ρCμk2/ε，Cμ=0.09.

k方程：

(4)

ε方程：

(5)

式中：C1ε=1.44；C2ε=1.9；σε=1.3.

1.2 烟雾颗粒数值模型

应用Fluent软件中嵌套的离散相模型模拟烟雾颗粒的扩散运动。离散相模型是遵循欧拉-拉格朗日法的两相流数值模型，该模型将流体视为连续相，在拉格朗日坐标系中对颗粒作用力微分方程进行积分来求解颗粒的运动轨道。在笛卡尔坐标系中，烟雾颗粒的作用力平衡方程和轨迹方程(x轴方向)描述如下：

(6)

(7)

式中：u1、up分别为x轴方向流体相速度和颗粒速度(m/s)；FD(ul-up)为颗粒的单位质量曳力函数；ρp为颗粒堆积密度(kg/m3)；gx(ρp-ρ)/ρp为颗粒的单位质量重力与浮力之差；Fx为单位质量的其他作用力。对(6)式在离散的时间步长上逐步积分，即得到颗粒轨道上每一个位置上的颗粒速度；对(7)式沿每个坐标轴方向积分求解，即得到颗粒相的轨迹。

2 物理模型和基本参数设置

数值求解前建立物理模型并对数值模拟的参数进行设置，物理模型和参数设置的准确性决定了数值模拟结果的可靠性。

2.1 物理模型与网格划分

烟雾生成过程是发烟材料燃烧的过程，而燃烧是复杂的化学反应。本文重点研究发烟罐烟雾在自由空间的扩散与质量浓度分布状态，对发烟材料的燃烧过程不做分析。因此，仿真模型简化为带有一定初始速度的高温烟雾从发烟罐底部涌出，并在自由空间扩散。构建的发烟罐物理模型和计算域如图1所示。

图1 发烟罐物理模型和计算域Fig.1 Physical model and computational domain of smoke pot

由图1可见，发烟罐模型为直径70 mm、高度400 mm的圆柱，上端开口，下端封闭并置于地面上。计算域为10 m×10 m×10 m的立方体。参考坐标系Oxyz如图1所示，发烟罐底部中心为坐标原点O.

利用ICEM软件对上述几何模型进行四面体网格划分，兼顾仿真精度和收敛速度。由于发烟罐表面附近烟雾的速度、温度梯度大，对其网格进行了加密；远离发烟罐的区域对计算精度要求不高，网格以一定增长率逐渐变大，以减小计算量。另外，在Fluent软件中将四面体网格转化为多面体网格，并对网格进行光滑操作[14]。计算域网格划分结果如图2所示，最终计算域网格数量为28 627个，其中网格质量优于0.71和0.66的网格分别占95%和99%.

图2 计算域网格划分结果Fig.2 Partitioning result of computational domain mesh

2.2 基本参数设置

发烟材料燃烧后生成烟雾颗粒，其粒径分布是烟雾的基本特性参数之一。将烟雾颗粒尺寸划分成一系列尺寸带，每个尺寸带中烟雾颗粒占颗粒总量的百分数就是烟雾的粒径分布。采集小型发烟罐燃烧释放后的烟雾颗粒样品，使用日本HORIBA公司生产的LA-950型激光散射粒度分布分析仪对粒径进行测试分析，得到烟雾颗粒的最小粒径为3 μm，最大粒径为26 μm，平均粒径为9 μm.烟雾的粒径分布如图3所示。

图3 烟雾颗粒的粒径分布Fig.3 Particle size distribution of smoke particles

对烟雾粒径进行最小二乘回归分析和显著性检验，得到烟雾颗粒的粒径服从罗辛- 拉姆勒(Rosin-Rammler)分布。罗辛- 拉姆勒分布又称R-R分布，是最常用的粉尘颗粒分布函数[15-16]，其表达式为

F(d)=1-exp[-βdn],

(8)

式中：β为特征参数；d为粒径尺寸；n为分布指数；F(d)为粒径小于d的累积百分比。根据烟雾粒径分布数据，可求解出n=3.45，β=3.09×10-4.

发烟材料的成分主要为磷与硝酸钾、硝酸钡、二氧化锰、三氧化二铁、高锰酸钾等构成的混合物，其燃烧释放的热量使得烟雾的初始温度可达到400～500 ℃[17].设定仿真条件如表1所示。

表1 仿真条件设置Tab.1 Simulation conditions

3 仿真结果和分析

3.1 无风情况下烟雾的扩散分析

首先分析无风情况下烟雾释放后的扩散状态，选择瞬态计算模型进行仿真。将发烟罐烟雾的扩散过程分成两个阶段：0～10 s，烟雾持续从发烟罐底部喷涌；10～100 s，停止喷涌，烟雾在自由空间中弥散。

沿z=0 m截面观察不同时刻的烟雾状态，烟雾质量浓度在0～10 s的分布如图4所示。由图4可以看出，随着烟雾的持续喷涌，烟羽体积逐渐增大，外形呈伞状。这是因为除迎面空气阻力外，烟雾没有受到障碍物和风力等干扰，烟羽中轴周围压力保持平衡，烟羽中轴烟雾质量浓度明显比周围质量浓度高。

图4 烟雾质量浓度在0～10 s的分布Fig.4 Distributions of smoke mass concentration at 0-10 s

为了研究烟雾温度对烟雾扩散速度的影响规律，分析3 s、5 s、8 s、10 s 4个典型时刻的空间温度场分布和垂直方向(y轴方向)空间速度场分布，结果分别如图5、图6所示。图5、图6中，垂直方向空间速度场表征了烟雾扩散速度，烟雾上升速度越快，扩散速度越快。

由图5、图6可以看出，高温烟雾从发烟罐出口喷出后与冷空气进行热量的传递，温度降低。对比同时刻的烟雾质量浓度分布和温度场分布可以发现，高温烟雾主要集中于烟雾中轴，而从中轴附近扩散开后温度基本下降到与环境温度相同。发烟罐出口正上方，持续的烟雾供给为烟羽中轴部位提供热量，因此中轴部位的烟羽温度高于环境温度并能够以较快的速度向上喷涌。烟雾从中轴附近扩散开后不能得到持续的高温烟雾供给，因此降温较快，上升速度也明显减慢。

图5 典型时刻的空间温度场分布Fig.5 Spatial temperature field distribution at the characteristic moment

图6 典型时刻的空间速度场分布(垂直方向)Fig.6 Distribution spatial velocity field in vertical direction at characteristic typical moment

取样分析3 s、5 s、8 s、10 s 4个典型时刻烟羽中轴温度、上升速度分别与高度的变化关系，绘制曲线分别如图7、图8所示。由图7可以看出：4个时刻的烟羽中轴温度随高度增加先迅速下降，然后下降速度趋缓，最后达到稳定状态；4个时刻的烟羽长度依次递增，造成4个时刻烟羽中轴温度随高度的下降趋势依次放缓。对比图7、图8可以看出，烟羽中轴温度、上升速度随高度的变化趋势基本同步。

图7 典型时刻烟羽中轴温度与高度的变化关系Fig.7 Relation between temperature and height in the middle axis of plume at a typical moment

图8 典型时刻烟羽中轴上升速度与高度的变化关系Fig.8 Relation between rising speed and height in the middle axis of plume at typical moment

对每个取样点的温度和对应的上升速度进行曲线拟合，得到相关性曲线如图9所示。由图9可以看出，烟雾温度与上升速度有明显的线性关系。

图9 温度与上升速度的相关性曲线Fig.9 Correlation curve between temperature and rising speed

烟雾质量浓度在10～100 s的分布如图10所示。由图10可以看出，烟雾停止喷涌后，烟羽持续扩散并呈团状。烟团中轴部分上升速度仍然较快，但是没有了烟雾补充，导致中间烟雾质量浓度低而周围烟雾质量浓度高，周围烟雾围绕中轴形成一个明显的烟环。随着时间的延续，烟团体积持续增大并呈云状。整体来看，烟云上升速度呈逐渐放缓的趋势。

图10 烟雾质量浓度在10～100 s的分布Fig.10 Distributions of smoke mass concentration at 10-100 s

3.2 风速对烟雾扩散的影响

发烟罐在户外开放空间使用，环境中的风极易对烟雾的扩散状态造成影响，进而影响发烟罐的遮蔽效果。对不同风速下的烟雾分布状态进行仿真分析，分别设定风速为0.5 m/s、1.0 m/s、1.5 m/s、2.0 m/s、2.5 m/s、3.0 m/s、4.0 m/s、5.0 m/s、10.0 m/s、15.0 m/s，风流方向沿x轴正方向。烟雾喷涌时间为10 s，仿真总时间为10 s.取不同风速下第5 s和第10 s典型时刻的烟雾质量浓度分布状态进行对比，沿z=0 m截面观察，结果如图11、图12所示。

图11 不同风速下5 s时刻的烟雾质量浓度分布Fig.11 Distributions of smoke mass concentration at 5 s at different wind speeds

图12 不同风速下10 s时刻的烟雾质量浓度分布Fig.12 Distributions of smoke mass concentration at 10 s at different wind speeds

由图11和图12可以看到，环境中的风使烟羽发生了倾斜，且随着风速的增加烟羽倾斜幅度增加。有利于烟雾较大范围弥漫于近地面附近，有效提高对地面保护目标的遮蔽效果。当风速小于2.5 m/s时，10 s时刻的烟羽明显长于5 s时刻的烟羽，且扩散范围更大，而当风速超过2.5 m/s后则差异不大，即5 s时刻烟雾已经达到了稳定扩散状态。这是因为风加快了烟雾的扩散，使烟雾在较短时间内达到稳定扩散状态，快速增大遮蔽范围。但是，当风速增加到一定程度后又会造成烟雾被急剧吹散，烟雾质量浓度太低，遮蔽效果下降。对比风速为4 m/s、5 m/s、10 m/s、15 m/s时的烟雾质量浓度分布，可以发现烟雾整体上的空间质量浓度依次降低。风速在2.0 m/s、2.5 m/s、3.0 m/s时，发烟罐释放的烟雾能在较短时间内达到稳定扩散状态，且在稳定扩散时具有较大的遮蔽范围和较高的空间浓度。因此，初步认为发烟罐的最佳使用风速在2～3 m/s附近。

在计算域空间中选择典型位置进行特征点采样，分析烟雾质量浓度与风速之间的关系。在z=0 m截面上沿烟雾扩散方向选取P(1.0 m，1.0 m，0 m)、P(2.0 m，1.0 m，0 m)、P(3.0 m，1.0 m，0 m)、P(4.0 m，1.0 m，0 m)4个采样点，得到5 s时刻和10 s时刻质量浓度与风速的关系曲线，如图13所示。

图13 不同时刻4个采样点的质量浓度与风速关系Fig.13 Relationship between mass concentration and wind speed at 4 sampling points at different times

总体来看，随着风速增大各采样点的浓度值先增大后减小。5 s时每个位置的质量浓度峰值，主要集中在风速为1～2 m/s时。10 s时每个位置的质量浓度峰值，主要集中在风速为2～3 m/s时。这也从数值上验证了发烟罐的最佳使用风速在2～3 m/s附近的结论，此风速下发烟罐烟雾能在较短的时间内达到较大的遮蔽范围和较高的空间质量浓度。

另外，对比风速为3 m/s、4 m/s、5 m/s、10 m/s、15 m/s 5种情况下烟雾已经达到稳定扩散状态的采样点质量浓度与风速的关系，为负相关，近似呈反比关系。这一结果符合文献[8]引用的高斯烟羽模式中下风向任一点的烟雾浓度与风速大小呈反比的表述。

3.3 质量流率对烟雾扩散的影响

z=0 m截面不同质量流率下10 s时刻的烟雾质量浓度分布状态如图14所示。由图14可以看到，随着质量流率的增加，烟雾的空间质量浓度明显增加。选取z=0 m和z=1.0 m两个截面，在每个截面上选取6个采样点，分析发烟罐质量流率变化对烟雾空间质量浓度造成的影响。

图14 不同质量流率下烟雾质量浓度分布状态Fig.14 Smoke mass concentration distributions at different mass flow rates

为与发烟罐的使用场景相结合，采样点选在计算域内烟雾扩散范围较大的区域。如图15所示分别是z=0 m和z=1.0 m截面上(4.0 m，1.0 m，0 m)、(4.0 m，2.0 m，0 m)、(4.5 m，2.0 m，0 m)、(4.5 m，1.0 m，0 m)、(4.0 m，1.5 m，1.0 m)、(4.5 m，1.5 m，0 m)6个采样点质量浓度值随质量流率的变化曲线。由图15可以明显看出，采样点的烟雾质量浓度值随发烟罐质量流率的增加而线性增加。由此可以认为在1～20 mg/s的质量流率范围内，发烟罐烟雾的空间质量浓度与质量流率呈线性关系。