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

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

0 引言

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

為了模擬爆炸型煙幕彈擴(kuò)散過程中的風(fēng)場特征和熱力效應(yīng)，本文根據(jù)試驗數(shù)據(jù)建立爆炸型煙幕彈初始云團(tuán)的數(shù)學(xué)模型，分析煙劑燃燒放熱對初始云團(tuán)的熱力效應(yīng)；運用CFD方法求解模擬區(qū)域風(fēng)場，并運用離散相模型進(jìn)行擴(kuò)散模擬；通過計算面密度對煙幕的有效遮蔽區(qū)域進(jìn)行模擬；將模型計算結(jié)果與試驗數(shù)據(jù)進(jìn)行對比，以檢驗?zāi)Ｐ偷挠行浴?/p>

1 材料、現(xiàn)場布置及試驗

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

圖1 裝藥結(jié)構(gòu)圖Fig.1 Diagram of charge structure

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

圖3 試驗中初始云團(tuán)Fig.3 Initial smoke cloud in the test

2 數(shù)值模擬

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

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

2.1 初始云團(tuán)數(shù)學(xué)模型

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

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

2)初始云團(tuán)為橢球體；

3)初始云團(tuán)內(nèi)濃度均勻分布，即粒子在橢球體內(nèi)均勻分布。

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

(1)

可以運用MATLAB軟件實現(xiàn)上述構(gòu)造初始云團(tuán)的方法。具體算法如下：

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

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

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

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

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

圖5 根據(jù)試驗數(shù)據(jù)生成的初始云團(tuán)示意圖Fig.5 Sketch map of initial smoke cloud generated from test data

2.2 風(fēng)場模擬

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

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

1)連續(xù)性方程

(2)

2)動量方程

(3)

3)能量方程

(4)

4)湍流脈動動能方程

(5)

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

5)湍流耗散率方程

(6)

式中：σε為湍動能耗散率的普朗特數(shù)；c1、c2為經(jīng)驗常數(shù)。

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

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

(7)

(8)

(9)

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

(10)

式中：CS為經(jīng)驗常數(shù)。

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

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

圖6 網(wǎng)格無關(guān)性檢驗Fig.6 Grid independence test

2.3 初始云團(tuán)的熱力效應(yīng)

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

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

2.4 擴(kuò)散模擬

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

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

(11)

(12)

(13)

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

(14)

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

2.5 面密度計算

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

(15)

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

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

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

(16)

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

3 結(jié)果分析

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

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

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

圖7 不同粒子數(shù)的計算結(jié)果Fig.7 Calculated results using different particle numbers

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

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

圖9 模擬數(shù)據(jù)與試驗數(shù)據(jù)對比Fig.9 Comparison between simulated data and test data

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

4 結(jié)論

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