田 益，馬向超，蔣佳麗

（西安電子科技大學，陜西西安 710071）

0 引 言

在現(xiàn)代軍事戰(zhàn)爭中，柴油燃燒所形成的煙塵是非常普遍的。例如：坦克、裝甲車等重型武器運行過程中，柴油經(jīng)過柴油機燃燒后會釋放出大量的煙塵；加油站和油井等能源設(shè)施也經(jīng)常被作為軍事打擊的目標；在某些特殊情況下，為了隱蔽目標，有時也會故意點燃油井、油桶等形成大量燃燒煙塵，以減少目標被命中的幾率，如海灣戰(zhàn)爭中伊拉克故意點燃油井躲避美軍紅外武器的攻擊。

隨著紅外成像導(dǎo)引頭技術(shù)［1］的快速發(fā)展，雖然采用第三代凝視紅外焦平面探測器的導(dǎo)引頭系統(tǒng)綜合性能大幅提高，但是燃燒煙塵等擴展干擾物依然是制約紅外制導(dǎo)導(dǎo)彈作戰(zhàn)效果的主要措施之一。當導(dǎo)引頭和目標之間存在燃燒煙塵時，由于其消光遮蔽作用，可破壞導(dǎo)引頭對目標的發(fā)現(xiàn)、識別和跟蹤，使得處于搜索狀態(tài)的導(dǎo)彈不能轉(zhuǎn)入跟蹤狀態(tài)；處于跟蹤狀態(tài)的導(dǎo)彈則會丟失目標，由跟蹤轉(zhuǎn)入搜索，或者導(dǎo)致其跟蹤誤差增大。因此，柴油燃燒所形成的煙塵不僅是己方裝備對抗紅外成像制導(dǎo)武器的潛在手段之一，也是己方裝備的導(dǎo)彈紅外導(dǎo)引頭所面臨的潛在干擾之一。

目前對柴油燃燒煙塵的研究主要在于對煙塵的成分進行物化分析，如柴油完全燃燒時的主要產(chǎn)物是二氧化碳和水，不完全燃燒時還產(chǎn)生一氧化碳甚至碳微小顆粒。

由于柴油的燃燒煙塵干擾普遍存在且容易釋放，本文擬通過對這種煙塵的物理化學特性、時空運動特性及其紅外輻射的空間、時間和光譜傳輸特性等進行研究，構(gòu)建燃燒煙塵的擴散運動和紅外輻射傳輸特性模型，為紅外制導(dǎo)武器抗燃燒煙塵干擾能力的評估提供基礎(chǔ)。

1 柴油燃燒煙塵的收集和微觀形貌分析

柴油燃燒煙塵的組成成分較為復(fù)雜，受燃燒條件的影響很大。我們使用煤油燈對柴油進行直接可控燃燒，然后采集其燃燒煙塵。

直接采集煙塵所使用的方法如圖1 所示，其目的是為了獲得用于微觀形貌分析的顆粒物。首先，將導(dǎo)電膠帶粘貼到電鏡樣品臺上，靠近燃燒煙塵排放通道，使煙塵顆粒直接流經(jīng)導(dǎo)電膠帶被沾附于樣品臺，從而獲得燃燒煙塵樣品。

圖1 燃燒煙塵采集示意圖Fig.1 Schematic of smoke collection from diesel combustion

圖2為所收集的燃燒煙塵樣品的掃描電子顯微鏡（scanning electron microscope，SEM）照片，顯示了500 nm 標尺下的表面形貌。從圖2 可以看出，燃燒煙塵主要是由小尺寸的顆粒組成，小顆粒物的形狀近似為球形，尺寸在20～40 nm之間。

圖2 柴油燃燒煙塵的SEM圖（標尺500 nm）Fig.2 SEM image of smoke from diesel combustion（scale ruler 500 nm）

2 紅外輻射傳輸模型

2.1 單個粒子消光特性

在建立煙塵的紅外輻射傳輸模型［2］前，需要知道單個煙塵粒子的消光特性參量，如散射截面σsca、消光截面σext、不對稱因子g等。為了簡化，根據(jù)實驗分析結(jié)果，我們將煙塵粒子視為球形粒子。對于球形粒子，這些參量可由Mie散射理論求解得到，也可采用離散偶極近似（discrete dipole approximation，DDA）算法計算?？紤]到通用性，本文采用DDA 算法來計算單個煙塵粒子在特定波長下的消光特性。

當波長不同時，粒子的復(fù)折射率會發(fā)生改變，因此粒子對輻射的消光特性會隨著波長變化。Chang等［3］對煙塵顆粒提出了一組對數(shù)擬合函數(shù)，復(fù)折射率的表達式為n′＝n＋ik，n和k的擬合函數(shù)如下所示。

式中：λ為入射光波長，單位μm。

根據(jù)對燃燒煙塵微觀形貌的分析，球形粒子半徑取20 nm，依據(jù)Chang 等提出的復(fù)折射率模型，在紅外波段上等間隔選取m個采樣點來計算對應(yīng)波長下單個煙塵顆粒的散射截面、消光截面以及不對稱因子，最終得出煙塵粒子在采樣點處的光學參量。

2.2 粒子系消光特性

燃燒煙塵是由多個粒子組成的彌散粒子系。根據(jù)粒子復(fù)折射率及尺度參數(shù)計算出單個粒子消光截面及不對稱因子后，再結(jié)合粒子系的物理參數(shù)（如粒子濃度、形狀及粒徑）可求出粒子系的消光及散射特性。

為了簡化模型，本文將粒子系中的粒子視為相同粒徑的球體，則單個粒子在某波長處消光截面σext乘以網(wǎng)格區(qū)域單位體積內(nèi)粒子數(shù)N的結(jié)果即為粒子系在該波長下的消光系數(shù)βsys：

由于在單個顆粒的消光中，有一部分是前向散射，依然會進入視場，因此，需要考慮前向散射的貢獻。Henyey 和Greenstein 根據(jù)多年研究，提出了一個與已知實驗數(shù)據(jù)較為符合的散射相函數(shù)經(jīng)驗公式，簡稱為HG相函數(shù)［4］，如式（4）所示。

式中：g是一個無量綱的控制散射相函數(shù)分布的常數(shù)，為單個煙塵顆粒的不對稱因子。當g＜1時，散射是各項異性的，g越大，分布越集中，前向散射占據(jù)的比率越大；當g＝1時，散射方向完全集中在前向。

單個網(wǎng)格內(nèi)粒子系在各角度散射的能量比率為

式中：∑σsca是流場計算中在某波長下單個網(wǎng)格內(nèi)粒子的累加散射截面；Δs是垂直入射方向的單個網(wǎng)格橫截面積；PHG（Θ）為該波長處單個煙塵顆粒HG 的散射相函數(shù)。當Θ 取0 時，F(xiàn)（Θ）即為單個網(wǎng)格內(nèi)粒子系的前向散射比率。

2.3 粒子系光譜消光特性

燃燒煙塵顆粒的光學特性具有強烈的光譜依賴性，粒子系在某個波段內(nèi)的消光特性通常是對光譜平均獲得的，其中常用的平均方法是普朗克平均法，即對之前計算的粒子系在該特定波段內(nèi)的光譜消光特性結(jié)果進行普朗克（Planck）平均，普朗克平均消光系數(shù)為

式中：f1為在紅外波段λ1～λm上分別等間隔選取m個采樣點后計算采樣點對應(yīng)的消光系數(shù)βsys，并將βsys數(shù)值進行擬合而得的4 階多項式函數(shù)；λ1和λm分別為波段積分下限與上限；Ebλ為特定波長范圍內(nèi)等間距所取的各個波長對應(yīng)的黑體的光譜輻射照度。

同理，普朗克平均前向散射比率可表達為

式中：f2為利用在紅外波段λ1～λm上等間隔選取m個采樣點后計算采樣點對應(yīng)的前向散射比率F（Θ），并將F（Θ）數(shù)值進行擬合而得的4階多項式函數(shù)。

上述計算過程可用圖3所示流程圖表達。

圖3 粒子系在不同波段的消光特性及前向散射計算流程圖Fig.3 Flow chart for calculating the extinction properties and forward scattering of particles in different wavebands

2.4 出射強度

沿厚度方向?qū)焿m平均分割為Ω份，每份的厚度為Δx，如圖4所示。

圖4 煙塵平均分割為Ω份的示意圖Fig.4 Schematic of dividing the smoke into N parts

出射強度為

其中：I0為入射紅外輻射強度。

3 時空擴散模型

3.1 物理模型

燃燒煙塵的擴散運動符合流體的運動規(guī)律，因此可以基于流體力學理論對煙塵的時空擴散運動進行建模?；谖锢磉^程的煙塵擴散運動仿真的理論基礎(chǔ)是求解其對應(yīng)的Navier-Stokes（NS）方程，可使用半拉格朗日法與壓力泊松方程來高效地求解方程，并使用渦量約束法［5］來降低由半拉格朗日法引起的數(shù)值耗散問題，從而準確地模擬出煙塵流場的湍流現(xiàn)象。

一般而言，煙塵的黏性可忽略不計。當煙塵的擴散速度遠小于聲速時，可壓縮性的影響也可以忽略不計，因此視煙塵為無黏性、不可壓縮流體［6］。求解煙塵流場的方程可表達為［6］

式中：ρ為密度；p為壓強；a為體積力（如重力，慣性力）的加速度；u為速度。

流場中煙塵顆粒受到的浮力與溫度有關(guān)，溫度較高的煙塵顆粒受到浮力的作用上升，通常用如下經(jīng)驗公式描述［7］。

式中：α表示溫度對浮力的影響系數(shù)；T為煙塵溫度；Tair為環(huán)境溫度；y表示浮力方向上的單位向量。

溫度在空間分布中不是均勻的，導(dǎo)致浮力的空間分布也不均勻，所以該浮力對速度場的影響不會被壓力項完全消除，會形成對流。

通過對溫度的計算可進而求出煙塵顆粒所受的浮力大小與煙塵的自發(fā)輻射強度，為此需要引入一個溫度標量場。溫度場不僅自身具有擴散效應(yīng)，速度場也會對其產(chǎn)生對流影響。溫度場可表示為

式中：cT為溫度冷卻系數(shù)；Tmax為流場最高溫度；Tair為環(huán)境溫度。

風與大氣條件和海拔高度有關(guān)，流場中用于描述風場的公式為

式中：u0是處于海拔為z0時測量出的風速；z是地面海拔；r用于描述大氣中風速隨海拔高度的變化情況。在實際仿真計算中取z0為10 m，并設(shè)r的默認值為0.3。

3.2 基于網(wǎng)格的煙塵流體仿真方法

不可壓縮流體控制方程可化為拉格朗日形式和歐拉形式，兩種形式在理論上是等價的。在CFD 中，由于計算精度有限，且兩種形式的離散化方式不同，會形成兩種截然不同的仿真方法：基于粒子的方法、基于網(wǎng)格的方法。

本文采用的流體計算方法是歐拉法，將流場計算區(qū)域劃分為離散網(wǎng)格。一個三維空間的規(guī)則網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)就是一個三維數(shù)組，只需要記錄數(shù)組的大小，通過下標就能完成對數(shù)組元素的訪問。本文的流場計算采用規(guī)則的靜態(tài)結(jié)構(gòu)化網(wǎng)格來描述。依據(jù)流場維度將計算區(qū)域劃分為X×Y×Z的三維立方體網(wǎng)格。

網(wǎng)格將空間均勻地劃分，仿真過程中速度、溫度、壓強、密度這些物理量都記錄在網(wǎng)格單元中心。為了簡化流場邊界處理，在實際仿真過程中需要在流場周圍額外加上一層網(wǎng)格［8］。因此在程序中，對每個物理量分配的數(shù)組大小為（X＋2）×（Y＋2）×（Z＋2）。

3.3 紅外煙塵渲染

紅外煙塵，特別是軍事上用于遮擋目標而使用的紅外煙幕，發(fā)煙器都會具有比環(huán)境溫度高的初始溫度，故需要考慮煙塵的紅外自發(fā)輻射。渲染紅外煙塵主要需要考慮的因素是透過率和對輻射的散射。煙塵粒子對光子的相互作用概率決定了其透過率，從視線方向穿透煙塵體的總透過率就是最終被渲染出來的α通道，α通道即在OpenGL 渲染引擎內(nèi)用于存儲透明度值的通道。煙塵的消光系數(shù)與消光截面和密度相關(guān)，輻射一部分會被外圍的煙塵吸收和散射，另一部分會透射出來并進入相機［9］。本文只考慮在視線積分路徑上的前向散射。

由于煙塵運動模型采用網(wǎng)格方法構(gòu)建，這種網(wǎng)格無法通過常規(guī)的光柵化方法進行渲染，故采用的是體繪制算法［10］，該算法能夠高效地處理煙塵體在視線方向上的透過率和前向散射，以達到實時性的需求。

本文采用了體繪制中最常用的光線投射算法，如圖5 所示，展示了從相機視口發(fā)出的射線穿過體紋理并進行采樣的過程。首先需要計算煙塵占據(jù)區(qū)域，當從相機引出的射線命中煙塵前表面時開始進行體繪制，光線傳輸方向與射線方向相反：射線不斷地在煙塵體積內(nèi)步進，不斷累計亮度和透過率，當射線離開長方體后表面，體繪制結(jié)束，其過程中累計的亮度和透過率就是煙塵的亮度和透過率。

將圖5 射線穿越體紋理的過程作為采樣合成過程，則合成公式可表達為

圖5 射線穿過體紋理［11］Fig.5 Rays passing through the volume texture

式中：Ci和Ai分別是在體紋理上采樣所得到的顏色值和不透明度，其實也就是體素中蘊含的數(shù)據(jù)；和表示累加的顏色值和不透明度。

基于本文所建立的紅外輻射傳輸模型，對體紋理的顏色值Ci與不透明度Ai有

式中：Iself為網(wǎng)格內(nèi)煙塵的自發(fā)輻射強度；Itrans為透射進來的輻射強度；對于煙塵自身紅外輻射，將其簡化為灰體輻射模型，且使其輻射強度與濃度成正比，該網(wǎng)格單元發(fā)出的自身紅外輻射強度為

式中：Ibλ為T溫度下λ波長的黑體輻射強度；λ1～λ2為所計算輻射的波段；ε為由煙塵粒子的等效發(fā)射率；N為網(wǎng)格內(nèi)的粒子數(shù)。

3.4 具體實現(xiàn)方法

針對煙塵的擴散運動模型，基于NS 流體力學控制方程，采用歐拉網(wǎng)格方法求解流場的擴散運動，主要包括：求解壓力項、計算對流、降低流場渦量和其他內(nèi)外力項的解算。將流場計算區(qū)域分為40×40×40個網(wǎng)格，對每個網(wǎng)格內(nèi)的煙塵運動狀態(tài)進行計算迭代，最終可得出網(wǎng)格內(nèi)的煙塵的速度、濃度、溫度等物理量。

針對煙塵的自發(fā)紅外輻射，首先建立煙塵的等效發(fā)射率模型，然后基于灰體模型建立煙塵粒子系的自身紅外輻射模型。針對煙塵紅外輻射的時間、空間、波段紅外輻射傳輸特性，采用DDA 算法，且分別在3～5 μm 以及8～12 μm 紅外波段等間隔選取40 個采樣點，計算單個煙塵粒子在不同波長處的光譜和波段光學參量，如散射與吸收截面、散射相函數(shù)等；繼而可計算粒子系在各波長采樣點處的消光系數(shù)和前向散射比例；然后，基于普朗克平均方法獲得煙塵粒子系在3～5 μm、8～12 μm 紅外波段的消光和散射特性；最后，采用對煙塵仿真區(qū)域進行等間隔采樣的光線投射算法，結(jié)合單個網(wǎng)格內(nèi)粒子系的消光、散射和自發(fā)輻射模型，從成像平面的每個像元引出一條射線，計算途經(jīng)每一條射線到達視點處的紅外輻射，如圖6所示。

圖6 煙塵輻射傳輸計算示意圖［12］Fig.6 Schematic for calculating the radiation transfer of smoke

依據(jù)以上一系列理論和方法，我們使用OpenGL圖形API 與CUDA 并行程序搭建整個煙塵的計算和渲染的框架，實現(xiàn)了基于歐拉網(wǎng)格方法的煙塵擴散運動和紅外輻射傳輸特性的聯(lián)合解算。

4 結(jié)果與分析

煙團的運動過程主要由煙塵源（煙塵發(fā)生器）類型和邊界條件決定，運動導(dǎo)致煙塵的溫度和密度隨時間變化，所以在某一時刻的溫度和密度的分布決定了煙塵的亮度和透過率。下面將結(jié)合仿真結(jié)果分析擴散時間、風場邊界條件以及波段對煙塵運動的影響，并檢驗其對恒定溫度黑體的遮擋效果。

4.1 風場對仿真結(jié)果的影響

圖7 顯示的是在5 m×5 m×5 m 的空間內(nèi)，無風條件下不同時刻煙的擴散圖。發(fā)射源持續(xù)釋放初始速度為1 m/s 且速度方向垂直向上的煙塵，煙塵初始濃度為0.1 kg/m3。發(fā)射源溫度恒定為350 K，環(huán)境溫度設(shè)為280 K，并且在流場區(qū)域放置一個溫度與發(fā)射源相同的黑體，在每個時刻對煙塵的擴散情況進行計算仿真。

圖7 無風情況下不同時刻煙的擴散圖Fig.7 Spread of smoke at different times without wind

將煙塵源位置設(shè)置為側(cè)面，煙塵初始釋放速度設(shè)為3 m/s 且方向設(shè)為水平方向，其它初始條件保持不變。在有風情況下煙塵的時空運動如圖8所示。

圖8 有風情況下不同時刻煙的擴散圖Fig.8 Spread of smoke at different times in the wind

由仿真結(jié)果可以看出，在風力作用下，煙塵的整體形態(tài)發(fā)生了變化，向著風場方向偏移。一般而言，風速越大，偏移越明顯。經(jīng)與無風情況進行對比還可看出，在有風與無風的條件下，煙塵對350 K 方形黑體遮蔽能力不同：在有風的情況下，煙塵迅速擴散，這導(dǎo)致對黑體的遮蔽時間和效果都明顯下降。

4.2 不同波段下的仿真結(jié)果

設(shè)煙塵初始釋放速度為1 m/s，釋放濃度為0.1 kg/m3。量化范圍上限設(shè)置為350 K所對應(yīng)的黑體輻射亮度，下限設(shè)置為280 K。波段分別取3～5 μm 和8～12 μm，并在場景中加入與煙塵初始溫度（350 K）相同的黑體方塊，用于觀察煙塵的遮蔽特性，得到仿真結(jié)果如圖9所示，左側(cè)是8～12 μm波段，右側(cè)是3～5 μm。

圖9 不同波段紅外煙塵Fig.9 Smoke in different infrared bands

將3～5 μm 和8～12 μm 波段中煙塵對黑體的遮蔽情況進行對比，可知：在8～12 μm 波段內(nèi)輻射透過率較高，煙塵的遮蔽能力較弱。紅外圖像中煙塵底端比較暗，這是因為底部煙塵濃度較高，輻射衰減很快；而在紅外煙塵頂端，煙塵濃度與溫度較低，輻射衰減較慢。

5 結(jié)束語

由于柴油燃燒煙塵在戰(zhàn)場中對目標的遮蔽特性，本文建立了柴油燃燒煙塵的紅外仿真模型。通過獲取煙塵形貌特征、建立紅外輻射傳輸模型、建立煙塵時空擴散模型等，完成了對紅外煙塵的實時仿真渲染，研究了不同風場下的柴油燃燒煙塵干擾的紅外輻射特性。研究結(jié)果可為紅外制導(dǎo)武器抗燃燒煙塵干擾能力的評估提供基礎(chǔ)。