粒徑對煙幕中團聚體顆粒散射特性的影響分析

周夢得， 李佳玉

(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)保持團聚體的分散程度相同并結合物理模型觀察發現，赤磷團聚體中的少量極大粒徑單體顆粒會對煙幕團聚體的消光性能產生較大影響，這也是造成數據波動的主要原因。

綜上所述，赤磷團聚體中單體粒徑的分散程度以及是否存在極大粒徑的單體顆粒對于煙幕的消光特性具有較大影響，在新型煙幕開發研制過程中，合理控制煙幕團聚體的粒徑分布，盡量避免生成具有極大粒徑單體顆粒的團聚體 對于增強煙幕的消光性能具有重大意義。