烟幕弹可以在短时间内形成半径数十米的烟幕，造成目视和红外屏蔽，有效对付热成像仪和激光测距仪的探测，在国内外均受到广泛关注和研究[1-4]。目前关于通过添加红外活化物质提高赤磷烟雾抗红外遮蔽性能的研究有很多[5-6]，但关于烟雾云团形成、扩散规律的研究相对较少。
在以烟幕弹为对象的相关研究中，M.E.Milham等[7]提出了磷衍生烟雾颗粒形成的描述模型，描述了磷衍生的烟雾消光对相对湿度的依赖性。朱晨光[8]针对赤磷烟幕的抗红外效能、对红外的遮蔽率展开了深入研究，同时，对烟雾云团的形成与扩散及微粒的运动性质进行了研究，提出了测试烟幕遮蔽性能的相对测量方法以及小波分析法，能够综合地评估红外烟幕的性能。潘功配等[5-6]研究了可膨胀石墨对发烟剂性能的影响，杜雪峰[9]等研究了真空中炭黑型发烟剂的红外消光能力,为发烟剂的改进提供了参考。王启超等[10]还针对碳类烟幕材料对太赫兹波的衰减特性进行了研究，结果表明在烟幕中加入碳黑可对入射太赫兹波进行遮蔽。
这里拟结合赤磷烟幕弹的工作过程，采用计算流体力学方法和数据处理方法对烟雾生成-扩散状态进行研究分析，建立烟雾流场的计算分析模型。计算模型从烟雾形成的两个阶段出发，将其分为初始烟雾场和发烟片飞散烟雾场分别进行计算，进而通过组合叠加方法，构建完整的烟雾场模型。对于发烟片飞散烟雾场，由于发烟片数量众多，同样采用叠加组合方式建立模型，即先建立单个发烟片的飞散轨迹模型和烟雾流场生成模型，再通过叠加组合的方式，构建完整的烟幕弹烟雾场模型。此外，本研究还对装药方式的影响进行了考察分析。
1 烟雾流场生成过程与方法建模
1.1 流场建立过程
一般情况下，烟幕弹被发射到指定空域后，点火装置在弹丸内部引燃扩爆药，当产生燃气压强超过弹壳的极限强度时，弹丸破裂并使得未经充分燃烧的发烟片四处发散，如图1所示。
图1 烟幕弹爆炸过程示意图
从图1可以看出，烟雾生成主要分布在两个阶段。第一阶段为弹壳破裂后，扩爆药爆炸产生初始的烟雾场；第二阶段为发烟片在飞散过程中，持续燃烧产生烟雾，形成垂帘烟幕。两个阶段产生的烟雾场叠加即为烟幕弹生成的最终烟雾流场。
1.2 建模分析方法
结合赤磷烟幕弹的工作过程，建模分析的主要内容分为3部分：
1) 结合扩爆药燃烧产物组成和弹丸破裂条件，利用非定常可压缩多组分流动扩散模型建立装药爆炸的初始流场模型，模拟初始烟雾场的形成和膨胀过程，并以达到初步稳定状态的流场状态作为初始流场状态。
2) 结合单个发烟片的飞散速度和燃烧释放产物，建立单个发烟片的运动和烟雾流场计算模型。其中发烟片运动采用动网格技术进行模拟；发烟片烟雾流场生成过程采用有限速率化学反应模型进行模拟，明确沿发烟片飞散轨迹的烟雾流场状态。
3) 假定每个发烟片沿其运动轨迹生成烟雾的状态相似，并且按一定的方式分散运动，如图2所示，采用密度叠加方式进行叠加组合，建立完整的烟雾流场结构。
图2 两种分布方式
分层均布式分散即将发烟片均分为多层，每层的发烟片再按一定角度均布。球体坐标系下均布即发烟片以烟幕弹为中心向空间分散，分散位置用立体角与平面角进行控制。
2 化学反应模型与发烟片轨迹模拟
2.1 含烟雾生成的流体计算模型
采用欧拉坐标系下的流动控制方程作为模型控制方程，控制方程组如下：
质量守恒方程：
(1)
动量守恒方程：
(2)
其中: ρ为密度；V为速度向量；v为速度在某一方向上分量的大小；τ为切应力；g和T分别代表作用在微元体上的重力体积力和其他外部体积力。
能量守恒方程：
(3)
式中，
代表流体微团的总能。
发烟片飞散过程中向外释放Cu和P2O5等燃烧产物，其中P2O5与大气中的水蒸汽产生化学反应生成的磷酸或偏磷酸，是发烟片烟雾流场的重要组成部分。针对这类化学反应，采用有限速率化学反应模型进行建模分析。
对于有限速率化学反应，混合物中各组分满足组分守恒方程：
(4)
式中： F为单位卷吸质量燃烧产热量；Vj为组分体积；Dim为质量扩散系数；
为质量变化率，可表示为：
(5)
式中：kf、kb为正、逆反应速率常数；Wj为组分j的分子量；
分别为反应物和生成物的化学当量系数[11]。
采用雷诺平均方法(RANS)对流动控制方程进行求解计算，并利用基于涡粘性假设的湍流模型对流动脉动项产生的雷诺应力和输运项进行封闭。计算中采用目前应用广泛、同时适用于可压缩和不可压缩流动的Realizable k-ε模型，求解带湍动能及其耗散率的输运方程。
2.2 发烟片运动及组分释放过程模拟
针对发烟片飞散的运动状态，可采用两种方法进行模拟。一种是采用自由来流方式模拟运动速度对烟雾流场的影响；另一种是采用动网格技术直接模拟发烟片的运动过程和运动轨迹。前者计算量相对较小，主要适用于定常气动模拟状态，后者适用于考察运动轨迹及流动的非定常发展过程，因此这里采用动网格技术模拟单个发烟片的飞散运动状态，并选择动态分层技术进行网格更新。
动态分层技术在边界上假定一个优化的网格层高度，在边界发生移动或变形时，如果临近边界的一层网格的高度同优化高度相比大到一定程度时，则在边界与相邻网格之间增加一层网格。相反，如果边界向计算域内运动，临近边界的一层网格被压缩到一定程度时，临近边界一层的网格将被删除。动态分层技术根据边界的移动量动态增加或减少边界上网格层，因此适用于结构化网格。
发烟片的运动如图3所示，即指定发烟片的初始速度，假定其在飞散过程中受到恒定阻力的作用，直至飞散速度为零。在此过程中，发烟片燃烧产生烟雾并向周围扩散。
图3 发烟片运动轨迹模拟示意图
3 烟雾粒子分布计算分析
3.1 实例模型与参数
1) 初始烟幕云团计算模型
烟幕弹在预定区域进行爆炸，由于其具有显著的轴对称结构，爆炸产生的初始流场将以弹丸为对称轴分布。因此采用轴对称模型建立初始烟幕云团计算模型如图4所示。
图4 初始烟幕流场计算模型示意图
由于最后的烟幕弹粒子分布状态由初始流场和单条发烟片的运动轨迹及粒子分布数据组合的方式来表征，所以初始流场与发烟片轨迹的模拟的准确程度相当重要。因此在实际模拟中，所有的网格都进行了必要的加密，以保证计算精度；以弹丸破裂时的状态作为初始条件，即初始压强设置为2 MPa，初始温度设置为1 500 K；在计算域的外边界，设置为标准大气条件。结合弹丸壳体厚度和相关研究，设置扩爆药燃烧生成的组分主要为H2、N2、CO和CO2(忽略含量小于5%的组分)，这些组分的体积百分比为10%、50%、10%和30%。
2) 单个发烟片分散运动与烟雾流场模型
研究所考察的烟幕弹发烟片尺寸为1 cm×1 cm，发烟片分散运动距离可达80 m。假定发烟片飞散运动2 s时间内，发烟片燃烧释放产物形成的烟雾流场径向扩散范围小于1.5 m，因此设置完整计算域为80 m×1.5 m，如图5所示。
图5 单个发烟片流场计算模型示意图
在计算模型中，发烟片以8 g/s的质量流量向外释放燃烧产物Cu、P2O5和NH3，释放产物的质量分数分别为0.3、0.6和0.1。模型中主要考虑P2O5与H2O之间的如下化学反应：
P2O5+3H2O=2H3PO4
(6)
P2O5+H2O=2HPO3
(7)
模型中假定Cu粒子直径5 μm，H3PO4粒子和HPO3粒子直径为1 μm。
3) 流场数据的叠加组合
在叠加模型中，将总计2 400块发烟片按照图2中的两种分散形式分布在空间里。叠加过程使用Matlab编写程序完成。该程序先读取完整的烟雾流场网格坐标，并将扩展数据(粒子数，温度等)设置为0。然后读取爆炸初始数据，将数据赋值给与烟雾流场网格坐标对应的节点。最后读取发烟片数据，将二维的单个发烟片烟雾颗粒的数据转换成三维坐标系下的数据，并赋值给完整流场中与网格坐标对应的节点。
对于分层均布式分散，坐标转换如下:
r=x2+y2
(8)
X=rsin(α)
(9)
Y=rcos(α)
(10)
Z=(position_start,position_end)
(11)
其中：x、y为发烟片计算结果中的数据点坐标；X、Y、Z为坐标转换后的坐标。position_start与position_end为烟雾流场轴线方向的起始与结束的位置。α为平面角，变化范围为[-π,π]。
球体坐标系下均布分散的坐标转换如下：
r=x2+y2
(12)
X=rsin(θ)cos(β)
(13)
Z=rsin(θ)sin(β)
(14)
Y=rcosθ-0.5·g·r2/v2
(15)
式中：g为重力加速度；v为发烟片初始速度； β为平面角； θ为立体角。
3.2 典型状态分析
利用上述方法和模型，对某赤磷烟幕弹烟雾流场粒子随时间变化的状态进行计算分析，获得其典型状态如下。表1为与扩散时间所对应的烟雾半径。
表1 烟雾半径变化
扩散时间/s0.080.881.78烟雾扩散半径/m6.6745.1173.33
1) 分层均布式叠加计算结果
图6给出了点火后，Cu粒子在0.08 s、0.88 s、1.8 s时的粒子数密度分布。在点火工作初期，发烟片便释放了部分Cu粒子，这些粒子主要集中在弹丸所在中心位置附近。随着发烟片的运动，Cu粒子向外飞散，在0.88 s时，粒子覆盖范围达到半径50 m左右，1.8 s时达到半径80 m左右。
图6 Cu粒子数密度
图7给出了0.18 s时,H3PO4粒子数密度的分布状态。在0.08 s时，部分P2O5与大气中的水蒸汽相互作用生成了磷酸和偏磷酸粒子。随着发烟片的运动，H3PO4粒子数密度增加。0.18 s时，整个烟幕基本形成。
图7 H3PO4粒子数密度
2) 球体坐标系下均布叠加计算结果
图8给出了球体坐标系下0.08 s、0.18 s时均布叠加的计算结果。Cu粒子的数密度分布规律与使用分层均布式叠加类似。
图8 Cu粒子数密度
4 结论
1) 不同时段的计算结果较好地表现了烟雾生成的过程，以及烟幕云团中烟雾微粒的粒度分布，为烟雾云团的遮蔽测试提供了参考。
2) 计算结果清晰地表现了两种分散方式产生的不同形状的烟幕，为烟幕弹装药结构的设计和改进提供了参考。
[1] Ernst-Christian Koch.Pyrotechnic smoke screen units for producing an aerosol impenetrable in the visible,infrared and millimetric wave range[P].United States:US 6,578,492 B1,Jun.17,2003.
[2] Josef Schneider,Ernst-Christian Koch.Method of producing a screening smoke with one-way transparency in the infrared spectrum[P].United States:US 6,484,640 B1,Nov.26,2002.
[3] 高志扬,祝利,韩书键.烟幕弹对毫米波末制导雷达的有效遮蔽区域研究[J].弹箭与制导学报,2015,35(6):62-64,69.
[4] 方良,郝建滨,朱璐,等.基于对抗反舰导弹的箔条质心干扰建模与仿真研究[J].兵器装备工程学报,2019,40(2):59-61.
[5] 潘功配,关华,朱晨光,等.可膨胀石墨用作抗红外/毫米波双模发烟剂的研究[J].含能材料,2007(1):70-72.
[6] 关华,潘功配,周遵宁,等.可膨胀石墨发烟剂对毫米波衰减性能的实验研究[J].火工品,2004(2):1-3.
[7] MILHAM M E,ANDERSON D H,FRICKEL R H.Infrared optical properties of phosphorus-derived smoke[J].Applied Optics,1982,21(14):2501-7.
[8] 朱晨光.提高赤磷烟幕抗红外效能研究[D].南京:南京理工大学,2005.
[9] 杜雪峰,潘功配,曲家惠,等.炭黑型发烟剂在真空中的红外消光性能研究[J].真空,2016,53(6):63-65.
[10] 王启超,汪家春,赵大鹏,等.碳类烟幕材料对太赫兹波的衰减特性[J].红外与激光工程,2017,46(5):201-205.
[11] 牛青林,杨霄,陈彪,等.高速滑翔目标点源红外辐射特征模拟及可探测性分析[J].红外与激光工程,2018,47(11):160-167.