徐路程， 郝雪穎， 肖凱濤， 宋偉偉， 陳春生

(1.軍事科學院 防化研究院, 北京 102205; 2.國民核生化災害防護國家重點實驗室, 北京 102205)

0 引言

煙幕是應用最廣泛的一種無源干擾方式，它通過在空中施放大量的氣溶膠微粒來改變光輻射的大氣傳播特性，從而掩蓋要保護的目標[1]。煙幕因較高的效費比以及在對抗光電制導武器方面具有顯著成效，受到各國軍隊的普遍重視[2]。按照成煙方式分類，發(fā)煙器材可以分為爆炸分散型、燃燒反應型和機械噴射型，其中爆炸分散型采用爆炸方式成煙。初始云團是指發(fā)煙劑在爆炸動能的作用下迅速膨脹，當在空氣阻力作用下達至平衡時形成的一個球冠狀或球形煙團，此后煙團在風、湍流及溫差作用下在大氣中傳播與擴散[3]。關于初始云團的形成過程，學者們已經開展了很多研究工作[4-7]。煙幕仿真對于煙幕戰(zhàn)術使用有較強的參考價值，張楠等[8]結合無風條件的高斯模型和紅外煙幕消光模型建立了紅外煙幕三維空間的透過率模型；花超等[9]、邱繼進等[10]運用拉赫特曼擴散理論模擬瞬時體源的擴散過程，分別對煙幕干擾紅外導引頭和紅外制導反艦導彈的效果進行了量化研究。相比于假定風場均勻穩(wěn)定的高斯模型和拉赫特曼擴散理論，基于計算流體力學(CFD)方法能夠得到關注區(qū)域中風場的更多結構特征[11]，從而提高煙幕擴散和消光模擬的可信度。此外，拉赫特曼擴散理論中只考慮了風場的擴散作用，忽略了源的熱效應如煙劑燃燒放熱的影響。徐路程等[12]采用CFD方法和離散相模型對連續(xù)點源煙幕進行了數值模擬。

為了模擬爆炸型煙幕彈擴散過程中的風場特征和熱力效應，本文根據試驗數據建立爆炸型煙幕彈初始云團的數學模型，分析煙劑燃燒放熱對初始云團的熱力效應；運用CFD方法求解模擬區(qū)域風場，并運用離散相模型進行擴散模擬；通過計算面密度對煙幕的有效遮蔽區(qū)域進行模擬；將模型計算結果與試驗數據進行對比，以檢驗模型的有效性。

1 材料、現場布置及試驗

試驗中使用爆炸型煙幕彈為圓柱型裝藥，裝藥結構如圖1所示，中心擴爆藥柱為6.3 g鈍化黑索今(RDX)炸藥，周圍發(fā)煙劑的主要成分為紅磷、氧化劑、鈍化包覆劑、黏合劑，質量104 g，其中紅磷84 g. 掛彈高度為1.5 m. 試驗時大氣穩(wěn)定度為中性層結，氣溫10 ℃，相對濕度40%，平均風速為2 m/s. 場地下墊面為平坦的裸露地面，布設如圖2所示，在掛彈點上風向和下風向分別豎立標桿，攝像方向與風向垂直，高速攝影采用美國DEL Imaging Systems公司生產的MotionXtra HG-LE高速攝像機，拍攝參數為分辨率1 128×752，幀速率1 000幀/s. 通過攝錄像及圖像分析法[3]對煙幕的長度和高度進行計算。圖3所示為試驗中爆炸生成的初始云團，經測量，煙幕有效長度為4.91 m，高度為1.21 m.

圖1 裝藥結構圖Fig.1 Diagram of charge structure

圖2 試驗場地布設示意圖Fig.2 Schematic diagram of test site layout

圖3 試驗中初始云團Fig.3 Initial smoke cloud in the test

2 數值模擬

為便于后續(xù)風場和擴散的建模和計算，建立坐標系如圖4所示：x軸正方向指向下風方向，y軸正方向與攝像方向一致，坐標系為笛卡爾坐標系，z軸正方向按照右手定則豎直向上，爆心在坐標系中的坐標為(100 m,100 m,1.5 m)，逆推可得坐標原點的位置。

圖4 幾何模型及邊界條件示意圖Fig.4 Schematic diagram of geometric model and boundary condition

2.1 初始云團數學模型

當初始煙幕的體積相對于煙幕的擴散體積不可忽略時[3]，應該使用體源模擬發(fā)煙源。本文研究的爆炸型煙幕彈符合這種特征。在初始云團的相關理論研究中[4,13]，一般假設初始云團內部濃度均勻分布。擴爆后形成的初始云團近似為橢球狀，并且一般采用起始半徑和高度兩個參數來表示尺寸[3,14-15]。初始云團的數學模型假設如下：

1)初始云團使用拉格朗日方法進行描述，即用粒子分布代表初始云團中的濃度分布；

2)初始云團為橢球體；

3)初始云團內濃度均勻分布，即粒子在橢球體內均勻分布。

令初始云團橫向半徑為r、高度為2h，爆心位置為(xc,yc,zc)，則構成初始云團的所有粒子坐標(x,y,z)滿足如下橢球方程：

(1)

可以運用MATLAB軟件實現上述構造初始云團的方法。具體算法如下：

1)通過MATLAB軟件隨機生成分別在[xc-r,xc+r]、[yc-r,yc+r]、[zc-h,zc+h]上滿足均勻分布的隨機數a、b、c；

2)如果(a,b,c)滿足上述橢球方程，則表明(a,b,c)是一個隨機生成的、橢球面以內的點，否則重新生成；

3)將(a,b,c)加入點序列(xi,yi,zi)；

4)重復上述步驟，直至點序列長度達到需要的總生成點數N.

取爆心坐標(xc,yc,zc)為(100 m,100 m,1.5 m)，r=2.455 m，h=0.605 m，運用以上方法生成N=10 000的初始云團示意圖如圖5所示。需要說明的是，在后續(xù)擴散模擬時，為匹配高分辨率的網格、保證濃度統(tǒng)計的精度，經無關性檢驗，實際仿真中使用的粒子數為N=400 000.

圖5 根據試驗數據生成的初始云團示意圖Fig.5 Sketch map of initial smoke cloud generated from test data

2.2 風場模擬

準確的風場模擬是氣溶膠擴散模擬的先決條件，本文采用CFD方法對風場進行數值模擬。氣溶膠在空間中擴散是一個三維過程，并且擴散過程一般都發(fā)生在大氣邊界層以內。因此本文中考慮垂直方向的計算高度為50 m，三維的計算區(qū)域為400 m(x軸方向)×200 m(y軸方向)×50 m(z軸方向)。

由于大氣邊界層中流速遠小于聲速，且初始云團形成后爆炸沖擊波對擴散區(qū)域的影響基本消失，可以假定大氣為不可壓縮理想氣體。描述不可壓縮黏性流體的控制方程采用Navier-Stokes(N-S)方程和能量守恒方程，湍流模型采用雷諾平均(RANS)方法中的標準k-ε模型(k為湍流脈動動能，ε為湍流耗散率)[16]。其中：

1)連續(xù)性方程

(2)

2)動量方程

(3)

3)能量方程

(4)

4)湍流脈動動能方程

(5)

式中：μt為湍流黏性系數；σk為湍流脈動動能的普朗特數。

5)湍流耗散率方程

(6)

式中：σε為湍動能耗散率的普朗特數；c1、c2為經驗常數。

以上方程組為非線性方程組，需要通過數值方法進行求解?？臻g上，對計算區(qū)域進行網格劃分，在爆心附近和地面附近等物理量梯度較大區(qū)域，局部進行網格加密。方程的求解使用Ansys Fluent軟件完成，動量、能量、湍動能和湍動能耗散率方程使用精度較高的2階迎風格式進行離散，壓力速度耦合采用穩(wěn)定性高的壓力耦合方程的半隱式(SIMPLE)方法。

對于入口邊界條件，Richards等[17]建議速度、湍動能、湍動能耗散率使用如下函數廓線表示：

(7)

(8)

(9)

式中：u*為摩擦速度；κ=0.42為卡門常數；zc=0.01 m為地面粗糙長度；cμ為經驗常數。對于下邊界，采用無滑移壁面邊界條件，需要特別注意的是，Ansys Fluent軟件中粗糙度用等效砂粒粗糙高度表示，它與上述空氣動力學地面粗糙長度之間的關系[17]為

(10)

式中：CS為經驗常數。

對于上、側邊界條件，Blocken等[18]建議采用與入口邊界保持一致的廓線，從而能夠最小地引入流向上梯度，保證流向上的均勻性。出口邊界采用自由流邊界。對于能量方程，除出口邊界以外，各邊界溫度均為283.15 K.

理論上，更密集的網格能夠得到分辨率更高的計算結果，本文采用3種網格劃分方式，分別將三維計算域劃分為10萬、15萬和20萬網格，并對風場進行穩(wěn)態(tài)計算。如圖6所示，不同網格數量情況下，在爆心位置處垂線的風場分布隨著網格加密，可以認為計算已經收斂，表明20萬網格已經能夠合理模擬風場特征，以下擴散模擬在此網格的風場中完成。

圖6 網格無關性檢驗Fig.6 Grid independence test

2.3 初始云團的熱力效應

爆炸燃燒成煙的煙幕云團因為燃燒釋放的熱量作用而溫度較高。因此當煙幕云團與周圍空氣達到壓力平衡時，煙幕云團的密度與周圍空氣相比較小，煙幕云團受浮力作用表現出上升運動[13]。此外，熱力作用還會增大湍流強度，使擴散更加強烈。

紅磷在空氣中與氧氣充分燃燒后生成十氧化四磷(P4O10)，與空氣中的水分相遇后生成具有不同結晶水的正磷酸。由于初始云團的形成過程時間較短，假定此過程中的放熱完全由紅磷與氧氣完全反應提供，此過程的反應熱可通過計算P4O10和紅磷的標準生成焓差值[19]計算得到，經計算為735.35 kJ/mol. 朱晨光[13]通過試驗發(fā)現，紅磷發(fā)煙劑燃燒產物分為兩個部分，一部分變成煙幕微粒，形成煙幕云團，另一部分燃燒產物結塊落在地面，落在地面的燃燒產物質量占50%. 根據第1節(jié)中發(fā)煙彈的成分組成，可以計算得到實際形成云團的紅磷燃燒放出的全部熱量為996.28 kJ，釋放出的熱量用于加熱初始云團。而發(fā)煙彈完成近場拋撒后一般在1～3 s內煙團溫度與大氣溫度逐漸達到一致[15]。于是，可以估算初始云團的平均溫度為297.45 K.

2.4 擴散模擬

為與拉格朗日方法描述的初始云團相匹配，煙幕擴散采用屬于拉格朗日方法的離散相模型[20]進行模擬。

根據牛頓第二定律，拉格朗日坐標下顆粒的動力學方程(以x軸方向為例說明)為

(11)

(12)

(13)

除以上力的作用外，顆粒還受到湍流擴散作用的影響，本文采用隨機游走(DRW)模型確定粒子的瞬時速度：

(14)

經煙箱中使用激光粒度儀測量，煙幕的平均粒徑dav=1.6 μm. 磷酸的密度取ρ=1 334 kg/m3[21].

2.5 面密度計算

對于網格統(tǒng)計得到的濃度的離散分布，面密度的計算可采用數值積分的方法實現(從y軸負方向進行觀察)，表示為

(15)

式中：Mik表示面密度矩陣M的分量；cijk表示煙幕濃度；Δy表示濃度cijk對應的光程。

需要說明的是，在計算風場時使用的網格尺度為米量級。從網格無關性檢驗結果來看，這個量級能夠滿足描述風場的需要。但在濃度統(tǒng)計及面密度計算時，因為源及擴散尺度都在1～10 m的量級，所以需要分辨率更高的網格來進行描述。本文中采用0.1 m的立方體網格進行濃度統(tǒng)計和面密度的計算。

對于一般的情況，網格3個方向的邊長分別為Δx、Δy、Δz，則網格濃度可表示為網格內粒子代表的質量和與網格體積之比為

(16)

對于紅磷煙幕，其遮蔽質量與空氣的相對濕度密切相關。根據文獻[21]，在相對濕度為40%時，紅磷煙幕的遮蔽質量為Mb=0.36 g/m2. 于是在求得面密度矩陣M的基礎上，通過計算面密度等值線，并且繪制遮蔽質量對應等值線以內的范圍，即可得到煙幕的有效遮蔽區(qū)域。

3 結果分析

下面運用第2節(jié)方法，對試驗中彈爆后10 s內的煙幕效能進行仿真計算。

通過大量的人眼和心理測試統(tǒng)計可知，只有當透過率小于1.25%時，目標才會被遮蔽，這個透過率與上述遮蔽質量是對應的[1,22]。因此，可以將計算得到的有效遮蔽區(qū)域與試驗中的煙幕遮蔽區(qū)域進行對比驗證。第1節(jié)試驗中煙幕對可見光有效遮蔽區(qū)域邊界的判讀參照文獻[3,23]中的方法：將某一時刻煙幕圖像與煙幕完全消散后的背景圖像對齊，并借助圖像分析軟件及鼠標以肉眼判讀來判讀煙幕圖像的遮蔽邊界，進而確定煙幕有效遮蔽區(qū)域的形狀、有效煙幕長度和高度。

與風場計算類似，理論上，更多的粒子數會提高計算精度，但同時也會帶來更大的計算量。分別采用N為100 000、200 000、400 000共3種粒子數進行仿真計算，煙幕的長度和高度的計算結果分別如圖7所示。由圖7可見，隨著粒子數增長，各時刻的計算結果均已收斂。因此，選用N=400 000的計算結果進行后續(xù)的結果分析。

圖7 不同粒子數的計算結果Fig.7 Calculated results using different particle numbers

2 s、5 s、10 s 3個時刻試驗圖像和模擬圖像對比如圖8所示，在試驗圖像中用紅色實線對煙幕遮蔽邊界進行了標識，對應的模擬圖像為同一時刻計算得到的有效遮蔽區(qū)域。不同時刻有效煙幕長度和高度模擬數據和試驗數據的對比如圖9所示。

圖8 不同時刻的試驗圖像與模擬圖像對比圖Fig.8 Comparison between test images and simulated images at different time

圖9 模擬數據與試驗數據對比Fig.9 Comparison between simulated data and test data

由圖8、圖9可見，模型能夠較好地反映煙幕形狀及長度、高度的變化規(guī)律。爆炸型煙幕彈在彈爆后由于受到自身的熱力作用，垂直方向抬升、擴散明顯；當抬升到一定高度后，由于能量散失，溫度下降，抬升作用減弱。同時，隨著離地高度的增加，風速明顯增加，動力效應作用凸顯，橫向擴散更加明顯。相比于拉赫特曼模式，本文模型的計算結果與試驗結果更為接近，表明本文模型能夠較好地模擬爆炸型煙幕彈前期擴散階段的作用效能。但隨著時間的增加，兩種模型的計算結果與試驗結果的相對誤差均逐漸增大，原因可能在于風向的波動沒有在定常模型中合理地反映出來[14]，也有文獻[24]指出即使根據多次示蹤試驗數據，同一地點同一類型天氣給出的擴散參數，其大小可相差多倍，其不確定度是非常顯著的，這個問題有待后續(xù)研究逐步完善。

4 結論

本文建立了一個爆炸型煙幕彈遮蔽效能仿真模型。運用拉格朗日方法，通過在橢球體內均勻生成隨機粒子的方法構建初始云團，通過分析初始云團的熱力效應結合計算流體力學方法和離散相模型，對粒子的擴散進行模擬。最后通過濃度統(tǒng)計和計算面密度得到了煙幕對可見光的有效遮蔽區(qū)域。與試驗結果和拉赫特曼模式結果對比，證明了本文所提模型的有效性，表明該模型能夠在爆炸型煙幕彈的前期擴散階段較好地反映煙幕形狀及長度、高度的變化規(guī)律，能夠為爆炸型煙幕彈的效能評估提供有效技術手段，為復雜環(huán)境下的戰(zhàn)術使用提供參考。