在现代战争中，随着精确制导武器的使用，对无源干扰的需求也与日俱增^[1-7]。在无源干扰中，烟幕占据重要位置，对烟幕作战效能的评估也成为研究热点。烟幕的作战效能与烟幕浓度及面密度紧密相关，计算烟幕浓度及面密度首先要知道爆炸云团的起始半径和高度，即烟幕云团初始参数。爆炸型烟源高度和半径的定义是：爆炸能量使所形成的烟幕云团膨胀扩展，与此同时能量逐渐散失，膨胀过程结束时烟团的最大高度称为初始云团高度，最大半径称为初始云团半径^[8]。关于烟幕初始云团参数的研究，朱晨光等^[9]建立了烟幕云团的膨胀模型，该模型假设烟幕云团膨胀过程始终受膨胀力和空气阻力作用；陈宁等^[10-11]建立了真空环境中烟幕云团形成阶段的膨胀模型，得到了烟幕云团在膨胀过程中体积及质量浓度与烟幕粒子运动速度运动时间的关系；本文中对烟幕膨胀过程理论模型作出改进：把云团的膨胀过程分为2个阶段，分别为等熵膨胀阶段和自由膨胀阶段，在此基础上建立烟幕云团膨胀的理论模型，该模型能够描述给定装置烟幕云团膨胀的基本规律，可将其用于爆炸发烟装置初始云团参数的计算。

1.   模型建立及爆炸过程分析

采用的模型为球形装药，配方是烟火药和轻质碳基干扰剂混合物。装药密度为1.1 g/cm³，其中碳基干扰剂单体(下文统称粒子微元)呈现多孔颗粒状，外形近似球体，半径为0.5 mm，密度为0.005 g/cm³。装药半径为13 mm，壳体材料为牛皮纸，壳体厚度为0.5 mm，采用中心点火方式，如图 1所示。

图  1  发烟装置模型截面图

Figure  1.  Model of smoke generator

下载: 全尺寸图片 幻灯片

发烟剂爆炸后，形成一个高温高压云团^[12]，其组分是气/固混合物。通常情况下，炸药的爆轰过程^[13-15]是非常短促的，因此，假定爆轰是瞬间完成的，即采用瞬时爆轰模型。基于瞬时爆轰假设，可使问题的研究有如下简化：(1)高温高压云团中气体为理想气体，第1阶段膨胀过程绝热等熵；(2)高温高压云团的膨胀过程视为一个不断扩大的球体，球体半径为r，质量为m；(3)假设有1个粒子微元始终处在云团边界，质量为dm，受产物膨胀力的作用面积为dS，粒子微元体积与云团的体积相比较足够小；(4)燃爆瞬间，t₀=0，初始云团半径r=r₀，第1个阶段的等熵膨胀完毕时t=t₁，云团的半径为r=r₁，当粒子微元速度变为零时t=t₂，云团的半径r=r₂。高温高压云团的膨胀过程分为2个阶段，分别为等熵膨胀阶段和自由膨胀阶段，如图 2所示。

图  2  烟幕云团及粒子微元受力分析示意图

Figure  2.  Schematic diagram of smoke cloud and force analysis on micro-unit

下载: 全尺寸图片 幻灯片

第1阶段为燃爆产物等熵膨胀阶段，在该阶段，粒子微元在炸药爆轰能量驱动下膨胀(由于爆轰能量驱动力远大于空气阻力和重力，此阶段忽略空气阻力、重力)，直至云团内部压力等于大气压时停止；

第2阶段为自由膨胀阶段，粒子微元只受重力和空气阻力作用(为了便于计算，暂时忽略重力)，直至在空气阻力作用下停止，此时形成的烟幕云团称为烟幕初始云团。

由粒子微元的受力分析得，其在第1阶段烟幕云团等熵膨胀时主要受到云团内部压力作用^[15]：

式中：p为云团压强，Pa。

根据上文假设，第1阶段为等熵过程，根据等熵过程理论有：

式中:p₀为高温高压云团初始压强，Pa；ρ₀为高温高压云团初始密度，kg/m³，κ为等熵指数。

将式(2)代入式(1), 得:

在第2阶段，粒子微元主要受到空气阻力的作用：

式中:C为空气阻力系数，ρ′为标准大气密度，kg/m³。式(3)~(4)分别为烟幕云团膨胀过程中第1、2阶段膨胀过程方程。

2.   基于龙格-库塔方法的模型计算

式(3)~(4)均为二阶非线性微分方程，一般说来不容易求出解析解，但可以通过数值方法求出其数值解^[16]。如龙格-库塔法^[17-19]，龙格-库塔法是一种间接采用泰勒级数展开而求解常微分方程初值问题的数值方法。其基本思想是利用在某点处值的线性组合构造公式，使其按泰勒展开后与初值问题的解的泰勒展开相比，有尽可能多的项完全相同，以确定其中的参数，从而保证算式有较高的精度。

以四阶龙格库塔为例，截断误差为R_h⁽⁴⁾=O(h⁵)，是关于步长h的无穷小量。下面给出最常用的四阶经典龙格-库塔公式：

首先，确定初始条件。根据理论模型，在REAL软件(各物质的物化参数在REAL软件的数据库中有存储)中进行计算，瞬时爆轰后，爆轰产物的温度T=1607.29 K，p′=12.89 MPa，气体质量m₁=6.60 g，固体质量m₂=3.52 g。因此，高温高压云团的初始参数为：云团压力p₀=p′=12.89 MPa，爆炸瞬间高温高压云团半径r₀=r′=13 mm。

然后，编写MATLAB程序，得出云团半径随时间变化结果如图 3所示。由于第1阶段膨胀时间极短，为了区别2个阶段云团膨胀规律，图 3(a)所示的第1阶段膨胀时间为0~7 μs，图 3(b)所示的第2阶段膨胀时间为0~1 s。由图 3(a)可以看出等熵膨胀阶段为变加速运动，在高温高压云团初始膨胀的第1阶段结束时，云团半径近似为42.3 mm，约为初始半径13 mm的4倍，这是因为在第1阶段中，粒子微元在爆轰产生能量的驱动下，粒子微元的加速度、速度迅速增加，导致云团半径的迅速增加。由图 3(b)可以看出，在第2阶段，云团半径仍持续增加。在其后由于粒子微元仅受到空气阻力的作用，粒子微元的速度变化逐渐变缓，云团的膨胀速度也逐渐变慢，直至约1 s时终止在100 mm附近，膨胀结束。

图  3  云团半径随时间的变化时程曲线

Figure  3.  Histories of smoke cloud radius

下载: 全尺寸图片 幻灯片

3.   烟幕云团参数实验研究

3.1   实验原理及方法

根据前文中的理论模型，加工烟幕发生装置，并将其吊装在固定架上，在室内条件下进行实验。采用"摄像法"测试云团的膨胀过程及初始云团参数，系统示意图如图 4所示。具体原理如下：通过摄像机记录烟幕成形过程，测距仪、测角仪测得距离角度参数，然后通过图像分析软件去除背景、确定烟幕边界阈值并二值化、去除图像上的“噪声”将被测对象提取出来。图像二值化就是将图像上的像素点的灰度值设置为0或255，也就是将整个图像呈现明显的黑白效果，这样做方便提取图像特征，有利于对图片做进一步处理。用Matlab中的bwarea工具获取二值图像的面积，然后求解云团半径^[8]:

图  4  测试系统示意图

Figure  4.  Schematic diagram of testing system

下载: 全尺寸图片 幻灯片

式中:$B=\frac{\alpha l}{{{57.3}^{{}^\circ }}}, \alpha =\left( \arctan \frac{b}{2f} \right)/{{b}^{\prime }}$；其中α为显示屏张角, b为成像面宽度，mm；f为摄像机镜头焦距，mm；b′为显示窗口半宽度，B为距离放大倍数, l为摄像点至源点距离，m; s为二值化后图像中云团面积，m²。

3.2   实验结果与分析

采用高速摄影机为SONY880E，其距离放大倍数为15。截取视频中0~0.35 s烟幕云团图像, 如图 5所示。以图 5中最后一幅图为例说明利用MATLAB对结果进行处理计算的步骤和方法：(1)对图像进行二值化，如图 6所示；(2)利用图像处理软件，去掉图像噪声，如图 7所示；(3)在MATLAB中应用bearea函数计算燃爆产物的面积，并求解此面积下的等效半径。

图  5  云团图像

Figure  5.  Picture of smoke cloud

下载: 全尺寸图片 幻灯片

图  6  图像二值化处理

Figure  6.  Image binarization processing

下载: 全尺寸图片 幻灯片

图  7  图像去除噪声处理

Figure  7.  Image interference removal processing

下载: 全尺寸图片 幻灯片

依据上述方法，对测得的图像进行处理，然后根据式(6)计算云团半径，并与理论计算曲线进行比较分析，如图 8所示。从图 8中可以看出，无论在云团膨胀的初期，还是在自由膨胀阶段，云团的半径变化实验测试值要比理论计算值小。主要原因有如下方面：第1阶段持续时间极短，为微妙级别，高速摄影机来不及捕捉烟幕膨胀图像；发烟剂未完全反应，放出的能量小于理论计算值。故烟幕云团半径变化实验测试值要比理论计算值小。可根据实验值对理论模型进行修正，使理论计算更加符合实际情况。

图  8  云团半径变化时程曲线

Figure  8.  Histories of smoke cloud radius

下载: 全尺寸图片 幻灯片

4.   总结与展望

本文中基于一种发烟装置，通过理论假设、建模分析、理论计算等方法描述了该装置烟幕云团的膨胀过程。通过实验结果分析可知，该方法能够描述该装置烟幕云团扩散规律。要进一步提高初始云团参数的计算精度，需考虑壳体破碎因素，如果能准确计算壳体破碎时高温高压云团的压强温度等参数，准确性将进一步提高。但该模型仅对发烟装置缩比模型进行研究，实际发烟装置尺寸比本文中模型尺寸要大，形状多是圆柱体。要把该理论运用于发烟装置烟幕初始云团参数的计算，还需考虑缩比效应、解决圆柱体爆炸与球体爆炸等效问题，这将在未来的工作中做进一步研究。