1 引言

烟火药通常是由可燃剂、氧化剂、粘合剂等粉末材料经过机械混合、造粒、压制等工艺制成的复合型含能材料, 是一种典型的非均匀多孔介质, 这种含能材料的燃烧具有高温、多相和复杂的中间成分, 因此建立烟火药火焰结构的研究历来是个难题^[1-3]。为此, 不少学者对于该问题进行了深入的研究与探讨, 埃利·弗里曼研究认为烟火药燃烧是从初始温度开始到转变生成火焰温度时的气体、液体和固体燃烧产物的复杂过程^[4]; 希洛夫认为烟火药燃烧火焰结构分还原层、完全燃烧层、热辐射层、氧气层和烟层等^[5]; 还有柏木(Kashiwagi)建立的均匀物系固体燃料的燃烧模型以及正田强半经验式延期药燃烧模型^[6], 都是经验的物理模型, 缺乏具体实施数据描述, 忽略了烟火药这种粉剂性材料燃烧时火焰中的燃烧粒子, 有很大的局限性。

近些年, 随着科学仪器的发展进步, 人们对烟火药燃烧的认识发生较大改变。一方面, 使用了大量计算软件, 以化学反应和热力学计算出发描述烟火药性能, 如Farid等^[7]通过不同侧面利用热力学理论对镁粉和聚四氟乙烯体系的燃烧机理和红外辐射效应进行了研究; 另一方面, 将流体力学两相流动理论引入烟火药燃烧火焰的计算与分析, 并借助粒子图像速度场仪测试设备(Particle Image Velocimetry，PIV)与高速摄影仪等设备, 发现烟火药燃烧火焰由气体流场与离散燃烧粒子组成, 且离散燃烧粒子是不可忽略的重要因素^[8-9]。改变了原有观点, 即认为烟火药燃烧火焰是一种与碳氢燃料一样的连续火焰, 分成内外焰结构, 炙热烟尘的影响很难分清火焰的观点。火焰结构诊断设备的不断更新使火焰结构研究取得了突破性进展, 如电荷藕合器件图像传感器(Charge Coupled Device, CCD)、平面激光诱导荧光(Planar Laser Induced Fluorescence, PLIF)、可调谐半导体激光吸收光谱(Tunable Diode Laser Absorption Spectroscopy, TDLAS)等的使用。但烟火药燃烧产生的大量烟尘, 明显干扰这些设备对烟火火焰诊断。感光设备所记录的图像及其强度与曝光时间直接相关。曝光时间长, 对于不动点光源, 体现为光强累加, 对于动点光源, 体现为迹线效应, 即在图像上表现为一个光线段(在曝光时间内动点光源移动距离的体现)。烟火药燃烧火焰含有大量高速运动的燃烧粒子, 在普通光学设备所成图像表现为连续的。先进高速摄影仪, 可以设置超短曝光时间, 不仅可以滤掉不发光的燃烧产物颗粒, 而且可消除动光源的迹线效应^[9]。尽管采用高速摄影仪(High-Speed Camera, HSC)可以拍摄火焰的大量燃烧粒子图像, 但这种图像并不能体现出粒子的速度与空间分布等火焰流场的重要参数。

为此, 本研究对烟火药火焰中的正在燃烧粒子进行实验与统计, 通过图片的二值化处理得出正在燃烧粒子的坐标位置, 通过连续图片的粒子追踪法得到速度分布, 获得其行为轨迹, 进而确立火焰粒子速度分布的数学模型。

2 实验 2.1 药剂与设备

所用的烟火药剂由高氯酸钾(KClO₄)、镁粉(Mg)和硝化棉组成。其中, KClO₄为分析纯, 粒径为80~100目, 分解温度为400 ℃, 国药集团生产; 可燃剂为雾化镁粉, 粒径为200~325目, 熔点648.5 ℃, 沸点1090 ℃, 唐山威豪公司生产; 硝化棉, 四川泸州245厂。

采用丙酮将硝化棉溶解, (按照质量比为8:92混合); 将高氯酸钾和镁粉混合均匀后, 与硝化棉溶液湿混, 并经过造粒、干燥, 称量后用30吨油压机分别压制在内径为18 mm和10 mm、高度为20 mm的单向开口的柱状铁壳内用于实验。其中高氯酸钾50%~55%, 镁粉35% ~40%, 硝化棉约10%。

HG100K高速摄影系统, 柯达REDLAKE公司, 170万像素CMOS感光器, 后触发连续记录, 可提供25~100000幅/s的帧比率, 图像最大分辨率为1504×1128。

2.2 实验

实验在暗室中进行, 装置如图 1所示。实验通过调整HSC与燃烧火焰距离, 保证燃烧火焰的燃烧质点充满高速摄影仪视场。根据燃烧质点的速度、亮度等因素, 确定曝光时间为15 μs, 高速摄影仪的帧频设置为500帧/s。

实验采用普通摄像机(最低照度0.1Lux)与高速摄像仪对直径为18 mm的烟火药柱的燃烧火焰进行研究, 对比火焰图像中粒子的解析效果, 结果如图 2所示。由图 2a可以看出, 普通摄像机的曝光时间处在十几毫秒到几百毫秒之间, 能记录火焰气体的光辐射与发光粒子辐射, 以及那些通过反射和散射的不发光粒子, 如图像中烟缕与背景, 记录的是整个火焰的外形图像, 但难以分辨火焰中的发光粒子。由图 2b~图 2e可以看出, 随着曝光时间的缩短, 不发光物质与发光粒子的迹线效应被消除, 烟火药燃烧火焰逐渐表现出高密度粒子特征。

图 2中烟火药柱直径为18 mm, 燃烧粒子密集度偏高, 为方便图像的处理, 选择直径为10 mm药柱进行处理, 取一组曝光时间为15 μs的图像作为示例进行研究分析, 结果如图 3所示。

从图 2b~图 2e与图 3可以看出, 曝光时间为15 μs时, 火焰大部分为粒子簇, 仅在喷口高温区域仍有显示一点火焰形状, 对其二值化处理, 适当选取阈值后, 这部分小火焰也变为粒子图像。

3 数学模型

对图 3进行分析可以发现, 图 3中的亮点部分(连续喷射的过程图片)为烟火药燃烧粒子。通过高速摄影仪慢放观察可知, 烟火药火焰中心轴线且在喷口位置的粒子较快, 周边速度随着直径的增大而逐渐降低。随着燃烧的进行, 这些正在燃烧的粒子将向上运动, 呈喷射状, 当粒子上升速度为零时, 由于重力的作用开始逐渐下降。由于高速摄影仪每张图片只能够给出每个燃烧粒子二维平面坐标和形状信息, 要想获得燃烧粒子的速度、轨迹及其流场, 还需要根据给出图面尺寸和连续拍摄时间间隔, 通过建立粒子追踪模型和轨迹模型, 可以求出粒子运动速度的变化过程, 进而建立该燃烧过程的数学模型。

3.2 建模思路

建模思路如图 4所示。

3.3 问题分析

(1) 二值化处理:依据文献[10]的运动轨迹提取算法对图片进行二值化处理, 高通滤波之后, 选择0.9为二值化的阈值比较合适, 影像的灰度高于0.9的令为1, 其余为0。直接使用MATLAB中DIP工具箱函数im2bw, 使用阈值变换法把灰度图像转换成二值图像, 为下一步用灰度值加权质心计算法^[10]求燃烧粒子的坐标奠定基础。

(2) 灰度值加权质心计算法:本研究的粒子定位采用灰度值加权质心计算法, 该方法不仅适用于对称图像的中心计算, 而且其优点还在于能够充分利用图像中每一点的灰度信息。加权质心计算公式如下:

$ \overline{x}=\frac{\sum\limits_{i=1~}^{m}{{{x}_{i}}\cdot f({{x}_{i}}, {{y}_{i}})}}{\sum\limits_{i=1}^{m}{f\left( {{x}_{i}}, {{y}_{i}} \right)}}, \overline{y}=\frac{\sum\limits_{i=1~}^{m}{{{y}_{i}}\cdot f({{x}_{i}}, {{y}_{i}})}}{\sum\limits_{i=1}^{m}{f({{x}_{i}}, {{y}_{i}})}} $

(1)

式中, x和y分别为质心的横坐标和纵坐标; m为粒子占据的像素数, 且m≥2; f(x_i, y_i)为第i个像素的灰度值。基于式(1)的灰度值加权质心计算法, 获取每张图片粒子的坐标, 并将其标注, 如图 5所示。同时, 对于二值化图像, 我们还得到了其中每个粒子的面积(pixel2), 这为粒子轨迹的求取提供了一组用以辨识的参数。通过文件批处理命令, 每张图片上粒子的特性矩阵D以*.mat格式输出, 便于下面MATLAB调用。

(3) 粒子标定与轨迹提取:粒子的标定是关键, 找到每两幅图中的对应粒子, 采用一组坐标值序列标定得到轨迹。这里, 采用面积和坐标作为特征量, 对粒子进行标定。

首先对第一第二幅图做粒子的初始标定, 从第一幅图片的粒子跟踪第二幅图片中的粒子, 并以此作为后续标定的初值。对于第一幅图中的某一粒子坐标为(x₀, y₀), 在其几何临域(这里采用扇形加矩形区域)内进行搜索, 下一幅图中所有处在该粒子邻域内的粒子均选为候选粒子; 对候选粒子逐个进行分析, 在候选粒子中标记变量不为1的且选择面积最相近的标定为该粒子下一时刻的状态, 并且置此两点标记变量flag为1以保证不会再次被标定; 由两组坐标分别相减得到粒子当前速度(pixel/f); 对第一张图中所有例子重复上述步骤。

初始时的大范围搜索图示, 此邻域判定条件的数学表达式:

$ \begin{align} &(\sqrt{{{\left( x-{{x}_{0}} \right)}^{2}}+{{\left( y-{{y}_{0}} \right)}^{2}}} <{{r}_{0}}\cap {\rm{arctan}}\left( \left| \frac{y-{{y}_{0}}}{x-{{x}_{0}}} \right| \right)<k)\cup \\ &\left| x-{{x}_{0}} \right|<{{x}_{l}}\cap {{X}_{i}}\cdot \left| y-{{y}_{0}} \right|<{{y}_{l}}(\sqrt{{{\left( x-{{x}_{0}} \right)}^{2}}+{{\left( y-{{y}_{0}} \right)}^{2}}}<{{r}_{0}}\cap \\ &{\rm{arctan}}\left( |\frac{y-{{y}_{0}}}{x-{{x}_{0}}}|<\mathit{k} \right)\cup \left| x-{{x}_{0}} \right|<xl\cap {{x}_{i}}\left| y-{{y}_{0}} \right|<yl \\ \end{align} $

(2)

然后用前面的结果作为条件对后续所有的图中的粒子进行标定, 算法如下:以当前粒子作为基准, 做出速度矢量, 矢量末端即为粒子下一时刻位置的预测; 把预测位置的邻域内所有粒子选为候选粒子。其圆形邻域数学表达为$\sqrt {{{\left({x + {u_0} -{x_0}} \right)}^2} + {{\left({y + {v_0} -{y_0}} \right)}^2}} < {r_1}$, 在候选粒子中选取面积最相近的粒子, 标定为该粒子下一状态; 重复进行预测, 标定, 直到下一时刻找不到候选粒子, 该粒子从图像序列中消失; 重复上述步骤, 遍历此时图片中每一点, 如图 5。

这里把图像进行倒序处理, 从后往前进行标定, 由于在粒子数密度较低的散射区域内, 粒子标定比较准确, 以此作为初值对高密度的火焰中心区域进行标定, 可以有效降低在高密度区内标定的偏差。

粒子标定之后得到了大量粒子坐标的时间序列, 把该坐标序在坐标系下按时序连接, 就可以得到粒子的运动轨迹。本研究通过MATLAB程序把这800张图中所有标记过的点的坐标、速度, 记录在矩阵里, 共记录了27313组数据。

3.4 建立模型 3.4.1 建模条件

(1) 假设所有燃烧粒子形状为球形。

(2) 火焰各区域温度变化较大, 但是燃烧粒子温度远高于其周围温度, 不仅燃烧粒子不断产生气体, 而且燃烧粒子呈现非刚性, 故假设忽略热泳力。

(3) 假设某取样时间内火焰中心单位时间发出的燃烧粒子数目恒定。

(4) 假设火焰中心发射的所有粒子, 在以火焰中心为顶点, 火焰中心轴线为轴线的一定圆锥角范围内, 且各个方向的单位时间内粒子发射数目相同。

(5) 假设经过一段时间, 系统为相对稳态。

(6) 假设流场气相流速与燃烧粒子速度差值为一常量。

3.4.2 计算结果

对燃烧粒子进行受力分析:主要有重力、浮力和曳力及其他影响。由牛顿第二定律可得下一方程:

$ \frac{{\rm{d}}{{u}_{p}}}{{\rm{d}}t}={{F}_{D}}\cdot \left( u-{{u}_{p}} \right)+\frac{{{g}_{x}}\cdot \left( {{\rho }_{p}}-\rho \right)}{{{\rho }_{p}}}+{{F}_{x}} $

(3)

式中，u_p为粒子速度, m·s^-1; t为时间, s; F_D, 曳力, N; u, 气相速度, m·s^-1; g_x, 重力加速度, m·s^-2; ρ_p, 燃烧粒子密度, kg·m^-3; ρ, 气相密度, kg·m^-3; F_x, 其它影响力, N。

上面方程右式第一项F_D·(u－u_p)为曳力, (${F_D} = \frac{{18u}}{{{\rho _p} \cdot d_p^2}} \cdot \frac{{{C_D} \cdot \mathit{Re}}}{{24}}, $${C_D} = {a_1} + \frac{{{a_2}}}{{\mathit{Re}}} + \frac{{{a_3}}}{{\mathit{R}{\mathit{e}^\mathit{2}}}}$), 右式第二项$\frac{{{g_x}\left({{\rho _p}-\rho } \right)}}{{{\rho _p}}}$为重力与浮力的合力, 右式第三项F_x为其他力的影响。

式中，d_p, 燃烧粒子直径, m; Re, 雷诺数; C_D, a₁, a₂, a₃为常数。

根据假设二, 燃烧粒子在流场中忽略升力, 布朗力, 热泳力对运动的影响, 故F_x=0。由于流场气相流速明显高于燃烧粒子, 所以可作假设, 流场气相流速与燃烧粒子速度差值为一常量, 记为C=u－u_p。

经过上述假设, 可得方程:

$ \frac{{\rm{d}}{{u}_{p}}}{{\rm{d}}t}={{F}_{D}}\cdot C+\frac{{{g}_{x}}\cdot \left( {{\rho }_{p}}\rho \right)}{{{\rho }_{p}}} $

(4)

根据上述受力分析, 结合牛顿第二定律, 得出下列常微分方程组:

$ \left\{ \begin{align} &\ddot{x}=\frac{{\rm{d}}{{v}_{x}}}{{\rm{d}}t}=-f\cdot {\rm{sin}}\theta \\ &\ddot{y}=\frac{{\rm{d}}{{v}_{y}}}{{\rm{d}}t}=-f\cdot {\rm{cos}}\theta-h \\ &{{v}_{x}}=\frac{{{\rm{d}}x}}{{\rm{d}}t} \\ &{{v}_{y}}=\frac{{\rm{d}}y}{{\rm{d}}t} \\ &{\rm{tan}}\theta =\frac{{{v_x}}}{{{v_y}}}=\frac{{\dot{x}}}{{\dot{y}}} \\ \end{align} \right. $

(5)

式中，${\ddot x}$, 水平方向加速度, m·s^-2; ${\ddot y}$, 垂直方向加速度, m·s^-2; ν_x, ${\dot x}$, 水平方向分速度, m·s^-1; ν_y, ${\dot y}$, 垂直方向分速度, m·s^-1; x, 横坐标, m; y, 纵坐标, m; θ, 速度方向与竖直方向夹角, rad。

根据初始条件t=0, v=v₀, R=0, 求解微分方程组, 可以求出x=x(t), y=y(t)的函数关系式。从而可得燃烧粒子的速度$v = \sqrt {{{\dot x}^2} + {{\dot y}^2}} = v\left(t \right)$。消去时间参数t, 可得v=v(x, y), 将燃烧粒子的坐标与速度建立了函数关系式, 得到燃烧粒子的速度区域分布。

3.4.3 问题求解

根据前文粒子的标定与轨迹的提取, 可以通过色彩、等高线、粒子数量等方式表述显示的粒子速度与密度分布, 如图 6所示。

从图 6a和图 6b可以看出, 粒子在火焰中心区域速度最大; 伴随着粒子的上升, 粒子的运动速度明显减小。火焰中正在燃烧粒子的速度分布有个梭型平衡区, 如图 6b内环区域内(每个等速等高线形成的近似图形), 而非预想粒子速度随着高度的增加, 速度迅速衰减。

从图 6c和图 6d可以看出, 粒子在火焰中心区域最集中, 数量所占比例最大; 随着燃烧粒子的上升, 区域粒子数量明显逐渐减少; 燃烧的粒子上升主要集中在一个以火焰为中心的锥角范围内, 在此范围外, 几乎没有粒子出现。

4 结论

基于HSC实验数据, 采用二值化图形处理方法, 以及粒子运动轨迹模型, 得出了烟火药燃烧粒子的粒子速度、粒子分布的计算模拟图, 并通过实验数据对数学模型进行了验证, 并得出如下结论:

(1) 采用灰度值加权质心计算法, 并通过倒序处理数据, 可以很好地解决烟火药粒子的燃烧轨迹;

(2) 喷口位置并不是粒子速度的最高区域。速度分布的等高线在内环区有个梭型平衡区, 证明正在燃烧粒子离开喷口后仍处在加速运动。

(3) 由于HSC可拍摄从点火到熄灭, 故模型可进行烟火药燃烧火焰的全过程分析。