周梦得， 李佳玉

(1.南京理工大学 能源与动力工程学院, 江苏 南京 210094; 2.电子设备热控制工业信息化部重点实验室， 江苏 南京 210094)

0 引言

随着激光制导和光电侦察技术的迅猛发展，使得战场透明度大大提高，目标隐蔽难度逐渐增大。为了争取防御方的战争主动权，光电对抗技术应运而生。烟幕就是一种简单且高效的光电对抗手段，并且由于烟幕成本低，效能高的特点而使其迅速成为现代高科技战场上最普遍的无源光电对抗手段之一[1-2]。

烟幕通过其中颗粒对入射激光的吸收和散射作用来达到降低激光武器效能的目的[3-4]。然而在烟幕云团的形成过程中，烟幕中颗粒由于相互碰撞发生凝聚现象，所以烟幕颗粒往往不是单颗粒而是以团聚体的形态存在的。由于团聚体是目前多种烟幕云团的基本组成单元，所以应先解决单个团聚体的散射及消光特性[5]。

为描述团聚体结构，国内外学者提出了许多物理模型，较常用的是分散限凝聚(DLA)模型和凝聚体- 凝聚体(CCA)模型。其中DLA模型以一个“种子”作为凝结核心，体系中颗粒随机运动并与之碰撞而发生团聚(包含颗粒与颗粒、颗粒与团聚体两种团聚过程)[6-7]；CCA模型则没有固定的凝结核，体系中颗粒随机进行布朗运动并发生碰撞团聚现象(包含颗粒与颗粒、颗粒与团聚体、团聚体与团聚体3种团聚过程)[8]。Li等基于DLA算法模拟生成了飞机羽流中粒径呈对数正态分布的氧化铝颗粒团聚体模型，并运用多球T矩阵(MSTM)方法计算得到了散射相函数、线性极化度随散射角的变化关系[9]。类成新基于CCA模型运用蒙特卡罗方法模拟计算了颗粒单体随机排布的烟尘团聚体的散射偏振特性，得到了团聚体的穆勒矩阵各元素随组成团聚体的颗粒数变化关系[10]。本文基于CCA理论建立了随机粒径(颗粒直径)的团聚体算法程序，此算法可以生成具有不同颗粒半径的团聚体模型，并基于此团聚体模型运用MSTM方法计算分析了烟幕中随机粒径团聚体颗粒的散射、衰减特性。

本文所计算的烟幕单体颗粒为表面光滑的球体，然而在实际烟幕中，团聚体的孔隙结构、表面粗糙度、表观密度等是客观存在的影响因素，这些因素对团聚体的散射特性具有明显的影响，采用有限元数值方法可以同步模拟计算团聚体的孔隙结构、表面粗糙等因素影响规律，但是目前缺乏上述影响因素的详细实验资料，因此可能会造成理论模拟结果较大偏离实际情况。本文所用MSTM理论是比数值方法更为精确的理论计算方法，可用于检验数值模型的精确性，它的局限之处就是目前只能处理单体颗粒为光滑球体的团聚体，如果用MSTM理论模拟实际团聚体，可根据实际清况，获得单体颗粒的等效光学参数(等效复折射率或等效介电常数)，然后再结合MSTM理论计算分析，由于本文篇幅有限，同时缺少烟幕颗粒的详细实验资料，因此计划先重点分析颗粒尺寸的影响规律，后期将分析其他因素对团聚体散射特性的影响程度，最终结合实际情况，同步考虑各类因素的综合影响规律。

烟幕的种类有很多，常见的有雾油烟幕、石墨烟幕、铜粉烟幕和赤磷烟幕等[11-14]。其中：雾油烟幕具有良好的空气悬浮稳定性，但对10.6 μm激光辐射干扰作用较小；石墨烟幕和铜粉烟幕虽对10.6 μm激光干扰作用明显，但沉降作用较为严重[11,15-16]；在多种烟幕中，赤磷烟幕不但对可见光，1～3 μm近红外光和8～12 μm远红外光具有良好的遮蔽效果，并且对不同波段的红外激光也有显著的干扰作用，且空气悬浮稳定性较强[17]。由之前的研究已知飞机羽流中的氧化铝随机粒径团聚体与等粒径团聚体在消光特性方面存在较大差异[9,18]，因此本文以粒径分散程度作为研究重点，并以赤磷团聚体为例深入探究粒径的随机分散程度对团聚体散射特性的具体影响程度。本文中颗粒一律指单个颗粒而非团聚体。

1 赤磷团聚体的物理模型

赤磷烟幕中的颗粒很多是以团聚体的形式存在的，并且由于颗粒自身的表面张力和剪切力作用，团聚体中单体数量一般在30～60个左右[18]。由王玄玉等[11]在进行赤磷烟幕对远红外激光衰减特性的测试实验给出的数据可知，赤磷颗粒的平均直径大概在0.5～13 μm之间，并且随时间增长，烟幕平均直径逐渐下降[17]。赤磷团聚体的粒径近似服从对数正态分布，其一般形式[19]为

(1)

(2)

(3)

本文采用CCA算法程序生成粒径呈对数正态分布的团聚体，图1通过流程图的方式简单介绍了算法的实现过程，图2是CCA算法生成的赤磷团聚体实例，且粒径都遵循均值为2 μm，标准差为0.8 μm的对数正态分布。

图1 CCA算法程序流程图Fig.1 Flow chart of CCA algorithm

图2 CCA算法生成的团聚体实例Fig.2 Aggregation models generated by CCA algorithm

2 MSTM方法

MSTM方法最早由Mishchenko等[20]提出，其在随机取向团聚体颗粒的散射特性计算方面具有极大的优势。MSTM方法的基本思想[21]是将散射体的入射场、散射场应用矢量球谐波函数展开，且由于麦克斯韦方程及其边界条件的线性化，可以运用一个传输矩阵来表示入射场展开系数与散射场展开系数之间的关系，这样的矩阵就是T矩阵。由于团聚体的T矩阵计算过程比较复杂，这里只做简要概述。团聚体的散射场由各个单体颗粒的散射场叠加而成，即

(4)

(5)

(6)

(7)

(8)

式中：Tj为团聚体中第j个单体的T矩阵；A(·)、B(·)表示分块矩阵的元素块，括号内容表示该块中所含元素；rl,j为rl-rj. 将(8)式反推可得

(9)

式中：Ti,j为将第i个单体的入射场展开系数转换为第j个单体的散射场展开系数的T矩阵。最终，可推得整个团聚体的T矩阵为

(10)

由此，可推团聚体的消光因子Qe、散射因子Qs、吸收因子Qa和散射相函数P，详细推导过程可见文献[22]。

3 模型验证与假设分析

此次模拟研究中用到了MSTM方法来计算团聚体的散射特性，并运用离散偶极子近似(DDA)方法来进行球形颗粒与非球形颗粒的的对比分析研究。因此为了验证计算模型的可靠性，现运用广义Mie(GMM)理论对MSTM方法和DDA方法进行对比验证。GMM是Mie理论在多球领域的拓展，具有极高的求解精度，因此作为验证模型可靠性的依据。

由于本文所采用的MSTM方法主要应用于计算随机排布的团聚体颗粒，DDA方法主要应用于计算单体颗粒，因此模型验证采用分开验证的方式。采用GMM分别验证两种计算模型在团聚体颗粒方面和单体颗粒方面的计算准确性。本次模型验证中，设定赤磷颗粒的复折射率为1.712+0.690i[18]、等粒径团聚体的单体数目为50、入射波波长为10.6 μm的远红外平面波。等粒径团聚体即指由相同粒径的颗粒组成的团聚体。

图3 运用MSTM方法和GMM方法计算的等粒径赤磷团聚体消光因子与粒径的关系Fig.3 Extinction factors and particle sizes of red phosphorus aggregates with the same particle size calculated by MSTM and GMM methods

图3为MSTM方法和GMM方法在计算随机排布的等粒径赤磷团聚体消光因子时的对比图。MSTM方法与GMM方法都是应用矢量球谐函数(VSWF)展开入射场、散射体内部的场及散射场来进行计算求解，因此误差较小，并且由图3可以看出，两条曲线吻合程度很好，最大误差为0.32%，由此可以证明MSTM方法计算团聚体的可靠性。

图4为DDA方法和GMM方法在计算赤磷球形单体颗粒消光因子时的对比图。由于DDA方法在颗粒尺度参数过大时会出现误差增大的情况，现讨论其在本文所涉及的尺度范围内的适用性。Draine等[23]根据粒子光学特性推导出的DDA方法适用条件：

|np|fd≤1,

(13)

式中：|np|fd设为A，np为颗粒的复折射率，f是入射波频率，d是偶极子棱长。本次计算偶极子数设定为60 000，颗粒的最大尺度参数为3.85，计算得到A最大为0.29，小于1，因此满足适用条件。观察图4可以看出两条曲线吻合程度很好，最大误差为0.84%，由此可以证明DDA方法计算赤磷单体颗粒的可靠性。

图4 运用DDA方法和GMM方法计算的赤磷单体颗粒消光因子与粒径的关系Fig.4 Extinction factors and particle sizes of red phosphorus monomer particle calculated by DDA and GMM

本文模拟研究的对象为赤磷烟幕团聚体，为计算模拟过程的可行性，现作如下假设：1)赤磷团聚体的单体颗粒为球形颗粒；2)赤磷颗粒团聚过程中复折射率不发生变化。

图5为利用DDA方法计算的具有相同体积的球体、椭球体和立方体赤磷单体颗粒的消光因子与等效直径的关系图，图6为椭球与球体、立方体与球体的消光因子差异随等效直径增加的曲线图，椭球体的3个半径之比为4∶4∶3. 由图5可以看出：立方体与球体差异较大，最大差异为10.74%；而椭球与球体的最大差异只有3.64%；但由于实际情况中方向的随机性以及自然界中的最小能量规则，故实际情况下颗粒多为椭球体与球体，因此将团聚体的单体假设为球形具有可行性。

图5 球体、椭球体和立方体赤磷单体颗粒消光因子与等效直径的关系Fig.5 Relation among extinction factors and equivalent diameters of red phosphorus monomer particles with three different shapes(sphere, ellipsoid (a∶b∶c=4∶4∶3), and cube)

图6 椭球体与球体、立方体与球体的消光因子差异与等效直径的关系Fig.6 Relations among extinction factor errors and equivalent diameters of ellipsoid and sphere, cube and sphere

赤磷烟幕中的团聚体在形成过程中，由于自身成分的变化，化学反应以及杂质等各种因素的影响，其复折射率会发生一定的变化，现就其对团聚体消光因子的影响进行分析。图7为复折射为1.712+0.690i的赤磷团聚体消光因子随折射率实部和虚部变化的曲线图。由图7可以看出，消光因子随复折射率变化呈线性增长趋势，但影响程度不大。赤磷团聚体复折射率实部与标准值差值为0.05时，误差达到最大为2.40%；赤磷团聚体复折射率虚部与标准值差值为0.05时,误差达到最大为1.23%. 由此可以看出，赤磷颗粒团聚过程中的复折射率的变化对团聚体的消光因子影响程度较小。因此假设烟幕颗粒团聚过程复折射率不发生变化以此来简化模型具有一定可行性。

图7 赤磷团聚体颗粒的消光因子与复折射率的关系Fig.7 Relation between extinction factor and complex refractive index of red phosphorus aggregates

4 结果与讨论

基于本文所建立的CCA模型并运用MSTM方法，计算了随机粒径且随机排布的赤磷烟幕团聚体在10.6 μm远红外平面波照射下的散射相函数、散射因子Qs和消光因子Qe，并对相函数在散射角上进行积分，得到了团聚体的总体散射能量Etot. 本次模拟中，设定赤磷团聚体的复折射率为1.712+0.690i[18]，团聚体单体数目为50，粒径分布函数为对数正态分布，其半径均值为1 μm，为能更细致地观察到粒径分散程度对散射特性的影响，将粒径标准差设定为3个不同的变化范围，分别是0.01～0.1 μm，0.1～1 μm，1～10 μm. 如图8所示。

图8 不同粒径标准差变化范围内相函数与散射角的关系Fig.8 Relation between phase function and scattering angle within the range of different standard particle size deviations

图8为标准差从0.01～0.1 μm、0.1～1 μm和1～10 μm变化时，团聚体的相函数随散射角的变化曲线，中间为相函数在散射角0°～15°范围内变化的局部放大图。由图8可以看出，当入射波长、团聚体中的颗粒数、复折射率等参量不变，仅改变粒径分布的标准差时，相函数的变化趋势保持不变，粒径的分散程度主要影响的是团聚体的前向散射，即散射角为0°～15°范围内的相函数。为能更直观体现随标准差的增加，赤磷团聚体散射特性的变化趋势，对相函数在散射角上进行积分，得出团聚体总体的散射能量Etot，如图9所示。

图9 总体散射能量与粒径标准差的关系Fig.9 Relation between total scattering energy and standard particle size deviation

图9是团聚体总体散射能量Etot随粒径标准差从0.01～10 μm的逐级变化曲线。由图9可以看出：标准差在0.01～0.1 μm范围内变化时，Etot没有明显变化，数据较为稳定；标准差在0.1～1 μm范围内变化时，Etot呈现出上升趋势，但曲线出现较大的波动；标准差在1～10 μm范围内变化时，Etot呈现出上升趋势，且曲线出现极大的波动。这表明团聚体粒径分布标准差以0.1 μm及以上量级增加时，团聚体的总体散射能量Etot受到较大影响，还表明粒径分布标准差变化的同时，出现了其他影响团聚体散射特性的因素，导致曲线出现波动。

图10为赤磷团聚体的散射因子Qs、消光因子Qe、吸收因子Qa随标准差从0.01～10 μm的逐级变化曲线。由图10可以看出，随标准差的逐级增加，消光因子Qe在均值3.83附近产生较大数据波动，最大浮动为25.5%，吸收因子Qa呈现下降趋势，散射因子Qs呈现上升趋势，且吸收因子Qa与散射因子Qs都出现较大的数据波动，针对此种情况将结合具体的模型结构进行详细分析。

图10 消光、散射、吸收因子与粒径标准差的关系Fig.10 Relation among extinction factor, scattering factor, absorption factor and standard particle size deviation

图11为运用本文所建立的CCA算法生成的赤磷团聚体模型结构图和粒径分布柱状图，其中图11(a)～图11(h)、图11(j)为A型团聚体，图11(i)、图11(k)～图11(o)为B型团聚体。为了研究随粒径分散程度的增加，团聚体总体散射能量Etot和消光因子Qe曲线出现巨大波动的具体因素，模拟生成了多组具有相同均值、标准差和单体数量的团聚体模型，设定入射光为10.6 μm的远红外平面波，复折射率为1.712+0.690i，单体数量为50，对数正态分布的半径均值为1 μm，标准差为0.8 μm，图11中，Etot为对相函数在散射角上积分得到的总散射能量。

由于消光因子Qe是衡量烟幕消光能力的主要指标，图11中，为了便于观察粒径分布和团聚体结构对于烟幕消光特性的影响，赤磷团聚体结构图和粒径分布图皆以消光因子Qe从大到小顺序排列。为了叙述方便，这里设定半径0～2 μm为小粒径颗粒，半径2～5 μm为中、大粒径颗粒，半径5 μm以上为极大粒径颗粒。图11中，所有的团聚体模型大致可分为两种：一种为由小粒径颗粒和少数中、大粒径颗粒组成的团聚体，这里设为A型团聚体，如图11(a)～图11(h)、图11(j)；另一种为存在少数极大粒径颗粒单体的团聚体，这里设为B型团聚体，如图11(i)、图11(k)～图11(o)。观察图11中(a)～图11(o)可知，当团聚体粒径分布的均值和标准差相同时，A型团聚体的消光因子Qe普遍高于B型团聚体，其差值最高达到了1.40，占最大值比值的30.9%. 团聚体的总体散射能量Etot也是衡量烟幕消光能力的重要指标，因此接着对比各个团聚体的总体散射能量Etot. 观察图11可以发现，A型团聚体的散射能量Etot普遍低于B型团聚体，其差值最高达到了222.83，占最大值比值的32.5%. A型团聚体与B型团聚体的主要差别在于B型团聚体中存在少数极大粒径单体颗粒，由此可以推断，团聚体中的极大粒径单体会对其消光特性产生严重影响，这也是曲线产生巨大波动的主要原因。

当团聚体中出现大粒径单体颗粒时，其消光因子会下降的原因分析如下：由图5烟幕单体颗粒消光因子与粒径的关系图所示，烟幕单体颗粒的消光因子随粒径增加在粒径5～6 μm处达到极大值，之后逐渐下降，但总体降幅较小，因此粒径大于5 μm的烟幕颗粒单体具有较强的消光性能，而粒径小于5 μm的烟幕颗粒单体普遍消光性能较弱。本文设定团聚体中的烟幕单体颗粒粒径均值为2 μm，标准差为0.8 μm，因此团聚体中出现粒径超过10 μm的极大粒径单体颗粒时，会导致团聚体中的单体颗粒粒径分布更多地集中在2 μm处，而粒径大于5 μm的烟幕单体颗粒总量减少。由于团聚体的总体消光因子由各个单体颗粒消光因子叠加而成，处于消光因子峰值区颗粒减少会直接导致团聚体的总体消光因子下降。因此在烟幕制备过程中，尽量避免生成具有极大粒径单体颗粒的团聚体，可以有效增强烟幕的消光性能。

5 结论

本文基于CCA理论建立了随机粒径的团聚体算法程序，使其可以生成具有随机粒径、随机排布方式的团聚体颗粒模型，并将其应用于计算烟幕中赤磷团聚体的消光和散射特性。着重研究了粒径的分散程度对于团聚体消光散射系数、散射相函数以及总散射能量的影响，并针对研究过程中出现的数据波动现象，结合具体的团聚体物理模型，进行了深入理论分析与讨论。得出主要结论如下；

1)随着团聚体颗粒粒径分布标准差的逐级增加，赤磷团聚体的消光因子Qe在均值3.83附近产生较大数据波动，最大浮动为25.5%.

2)赤磷团聚体的总体散射能量随着团聚体粒径分散程度的逐级增加，呈现出了逐渐上升的趋势，并且同消光系数一样，团聚体的总体散射能量也出现了较大的数据波动。

3)保持团聚体的分散程度相同并结合物理模型观察发现，赤磷团聚体中的少量极大粒径单体颗粒会对烟幕团聚体的消光性能产生较大影响，这也是造成数据波动的主要原因。

综上所述，赤磷团聚体中单体粒径的分散程度以及是否存在极大粒径的单体颗粒对于烟幕的消光特性具有较大影响，在新型烟幕开发研制过程中，合理控制烟幕团聚体的粒径分布，尽量避免生成具有极大粒径单体颗粒的团聚体 对于增强烟幕的消光性能具有重大意义。