不同風場下TNT 炸藥爆炸煙云的擴散模型及特性*

段中山，龔朋彬，袁 偉，過惠平，羅永鋒，羅昆升

（1. 火箭軍工程大學核工程學院，陜西 西安 710025；2. 中國人民解放軍陸軍勤務(wù)學院環(huán)境工程教研室，重慶 401331；3. 火箭軍研究院，北京 100094）

炸藥爆炸煙云抬升高度變化與爆炸當量密切相關(guān)，建立煙云擴散高度變化模型可用于反演未知爆炸事故類型、當量或推測戰(zhàn)場武器型號、威力及毀傷效應(yīng)，實現(xiàn)未知爆炸源實時快速偵測。而對于臟彈、核武器炸藥化學爆炸或其他“生化”爆炸，此類爆炸往往具有高度隱蔽性、煙雨危害性和應(yīng)急時效性，事故初期難以通過調(diào)查現(xiàn)場炸坑痕跡、玻璃受沖擊痕跡和爆炸震動記錄等傳統(tǒng)反演方法獲取其爆炸當量及煙云高度[1]。由于移動攝像、視頻監(jiān)控、地面遙感、衛(wèi)星觀測等影像獲取手段的使用與普及，使得利用爆炸視頻中煙云高度變化信息快速反演炸藥當量、預(yù)測污染煙云最終高度成為可能?？傊?，研究煙云擴散高度變化規(guī)律對戰(zhàn)場爆炸遠程評測、爆炸事故非現(xiàn)場偵測、有害物質(zhì)爆轟污染快速預(yù)警與評估等都具有重要價值。

為獲取復(fù)雜風場下爆炸煙云擴散高度變化過程及最終高度，本文擬采用仿真與實驗相結(jié)合的方法開展研究。首先采用CFD 方法計算煙云時空分布模型；然后利用煙云時空分布實驗結(jié)果驗證建模方法和參數(shù)設(shè)置的正確性，確定一套可表征煙云擴散時空分布的建模與參數(shù)設(shè)置方法；最后建立不同風速下煙云擴散高度變化模型和最終高度計算方法并討論風對煙云高度影響過程與機理。

1 爆炸煙云擴散理論

1.1 爆炸煙云擴散物理過程

炸藥爆炸煙云擴散過程可分為4 個階段[4]。階段1：爆轟階段，炸藥引爆后爆轟產(chǎn)物不斷膨脹直至地面煙團直徑達到最大值，炸藥參數(shù)、裝置結(jié)構(gòu)和地面環(huán)境決定了初始煙團形態(tài)、尺寸、密度、溫度及垂直動量等參數(shù)，該過程幾乎不受大氣條件影響。階段2：熱抬升階段，地面煙團在大氣浮力作用下將脫離地面開始熱抬升，此階段煙云平均溫度高于環(huán)境溫度，大多存在化學燃燒和地面反沖沖量，持續(xù)時間幾秒到十幾秒。階段3：持續(xù)上升階段，浮力和慣性作用將煙云以擴散與湍流的形式抬升至穩(wěn)定狀態(tài)，持續(xù)時間達數(shù)十秒至上百秒。階段2 和階段3 可統(tǒng)稱爆炸煙云主動抬升過程，過程中水平風和大氣穩(wěn)定度影響煙云與大氣熱交換速度、空氣混入速度等，進而改變煙云高度分布[11]。階段4：被動擴散階段，煙云穩(wěn)定后幾乎不再上升，隨著大氣風場被動擴散并輻射下游區(qū)域。

1.2 爆炸煙云擴散理論方程

爆炸煙云主動抬升過程中主要受重力、上升阻力和大氣密度差帶來的浮力等影響，其擴散過程可以用N-S（奈維-斯托克斯）方程描述[11]：

式中：ρ 為流體密度，煙云簡化為不可壓縮時可認為ρ 為常數(shù)；u 為流體速度；p 為壓強；μ為流體黏性系數(shù)；f 為外力項。

式(1)為流體質(zhì)量守恒方程；式(2)為流體動量守恒方程，右邊四項分別為外力項、壓強項、黏性項、對流項。

式中：a、b 為系數(shù)項，T 為煙云溫度，Tatm為環(huán)境溫度，y 垂直向上向量。

煙云密度場和溫度場受到煙云速度場的傳送作用，可用如下關(guān)系表示：

2）其他土地利用類型都有不同程度的減少，其中水域和耕地減少的最為明顯，1985－2000年和2000－2016年水域的減少速度分別為0.08%和0.19%，耕地的減少速度分別為0.12%和0.14%；

煙云擴散過程中存在熱量交換，遵循能量守恒傳遞方程：

式中：E 為流體微團總能量； hj′、 Jj′分別表示組分 j′的焓和擴散通量；keff為有效熱傳導(dǎo)系數(shù)；Sh為熱源項，包括了化學反應(yīng)熱及其他體積熱源項。式(5)右邊前3 項分別為熱傳導(dǎo)、組分擴散和黏性耗散能量輸運。

煙云擴散N-S 方程待求解變量為速度u、密度ρ、溫度T 和壓強p。CFD 計算方法一般通過劃分待求解區(qū)域網(wǎng)格來求解N-S 微分方程組，計算后可得各網(wǎng)格點上各時刻物理變量的分布場。

2 爆炸煙云擴散實驗與仿真

2.1 煙云擴散實驗設(shè)計與實施

實驗選取1、16、62 kg TNT 先后開展5 組外場煙云擴散實驗研究。實驗分組如表1 所示，其中1 kg TNT 無風實驗用于近距離觀察煙云擴散形態(tài)演化過程及其特征，16 和62 kg 實驗用于驗證無風與典型風場下的煙云擴散時空分布參數(shù)及模型。

選擇陰天傍晚時刻進行實驗以保證大氣處于穩(wěn)定狀態(tài)，對TNT 裸裝藥進行地面引爆。實驗布局如圖1 所示，主要過程可分為4 步：(1) 設(shè)計1、16、62 kg 的TNT 炸藥化學爆炸裝置，為盡可能長時間地記錄煙云，在實驗爆炸裝置外圍均勻布置少量煙云示蹤劑；(2) 對于1 kg TNT 爆炸，通過26 m 吊塔鋼絲線直接標定煙云高度，對于16 和62 kg TNT 爆炸，則通過釋放繩索牽引的高空氣球后，采用激光測距儀測量氣球高度，結(jié)合角度測量儀標定高空視野；(3) 根據(jù)測風儀測風結(jié)果選擇風速到達實驗要求時刻起爆裝置，同時記錄爆炸煙團初始形態(tài)及時空分布；(4) 根據(jù)視場標定結(jié)果，利用專業(yè)像素分析軟件對實驗煙云圖片信息開展計算，獲取爆炸煙云擴散時空分布數(shù)據(jù)。

表 1 實驗分組表Table 1 Experiment group table

圖 1 實驗場地布局Fig. 1 Experimental site layout

2.2 煙云擴散仿真思路與方法

煙云抬升演化過程主要指階段2 和3，階段1 可采用實驗經(jīng)驗公式結(jié)合有限元軟件AUTODYN 計算共同判定地面膨脹最大時刻非均勻的初始煙云形態(tài)、尺寸、密度、溫度及垂直動量，將初始煙云參數(shù)優(yōu)化后耦合到可編程軟件FLUENT 中構(gòu)建煙云擴散仿真模型[12-14]。

初始煙云參數(shù)：地面煙團最大膨脹半徑符合實驗經(jīng)驗公式r=1.93M0.32/ (T1/3 600)1/3，T1為爆溫（2 861 K）[9]；AUTODYN 爆轟計算模型顯示幾十毫秒內(nèi)初始煙云達到最大膨脹體積，由于地面影響初始煙團可簡化為半橢球形態(tài)耦合到CFD 模型中，1、16、62 kg TNT 的地面煙團直徑約為4.17、10.12、15.61 m，高約3.13,、7.59、11.71 m，火球膨脹后壓力與空氣壓力相當，整體符合ρ1=ρ0(T0/ T1)[9]，其中ρ1為煙團密度，ρ0與T0為空氣密度與溫度，ρ1均值取0.12 kg/m3，由于地面反沖效應(yīng)導(dǎo)致煙云存在約5 m/s 初始垂直動量，由于溫度迅速下降取整體煙團均值溫度約1 000 K，煙團組分按TNT 爆炸產(chǎn)物體積分數(shù)設(shè)定，假設(shè)初始煙團參數(shù)分布均勻且忽略顆粒攜帶、炸藥拋灑和后續(xù)化學燃燒[15-16]。

邊界與網(wǎng)格：在二維空間中采用從中間到兩邊、從下到上遞增的步進網(wǎng)格劃分方式，這樣中間和下部網(wǎng)格密集利于關(guān)注煙云擴散；計算中發(fā)現(xiàn)邊界條件設(shè)置對煙云分布有影響，綜合分析與試算后將文獻[12]中邊界設(shè)置變?yōu)樽筮吔鐬镮nflow、上邊界與右邊界均為Pressure Outlet、底部邊界Wall 計算更科學；為使風場更真實，采用UDF 編程風速隨高度變化規(guī)律來設(shè)置風場，風速u=U10(z/10)a，U10為地面10 m 處風速，大小取0～6 m/s（0～4 級），z 為高度，a 為地面粗糙系數(shù)（取0.3），在空氣入口處加載UDF velocity[17]，5 級及以上風速屬較極端氣象條件暫不做討論。

算法與參數(shù)：采用混合物耦合計算、湍流標準k-ε 模型、隨時間變化的非定常流計算[17-18]；計算中對x 向和y 向速度、能量、各物相體積分數(shù)及k 值等進行收斂設(shè)定和殘差監(jiān)視。

3 結(jié)果分析

3.1 煙云擴散仿真與實驗結(jié)果分析

3.1.1 煙云形態(tài)時空分布結(jié)果分析

圖2 為1 kg TNT 爆炸后不同時刻煙云分布實驗與仿真結(jié)果，其中實驗吊塔高度26 m，仿真空間高度為40 m。實驗較完整記錄了煙云時空分布及形態(tài)變化過程，前期煙云由于地面反沖動量及浮力作用加速上升、寬高比例較小，后期在浮力渦環(huán)主導(dǎo)下煙云擴散整體呈現(xiàn)“球形態(tài)”。煙團上端邊界分明，下端煙柱明顯，由于大氣不確定性、示蹤劑非均勻拋灑等影響導(dǎo)致后期煙云左側(cè)出現(xiàn)少許不規(guī)則形態(tài)，30 s 左右煙云基本穩(wěn)定后幾乎不再上升。計算結(jié)果則顯示由于煙云密度梯度差產(chǎn)生的浮力將煙團由高斯形態(tài)濃度分布轉(zhuǎn)變成由兩個反漩渦環(huán)流組成的渦流分布，橫截面上展現(xiàn)出由兩個反漩渦環(huán)流所組成的“腎形狀”模式，也就是常見的“蘑菇云”現(xiàn)象。實驗和仿真結(jié)果可見，實驗煙云擴散過程中具體形態(tài)不及仿真煙云對稱和規(guī)整，但兩者的形態(tài)演變過程和時空分布參數(shù)基本一致，驗證了爆炸煙云低密度浮力煙團擴散機理和雙渦環(huán)反漩渦環(huán)流擴散方式的正確性。

圖 2 實驗與仿真煙云時空形態(tài)對比Fig. 2 Comparison of time and space patterns between experimental and simulated clouds

3.1.2 煙云擴散高度分布實驗與仿真結(jié)果

實驗煙云早期能清楚顯示上邊界，由于空氣稀釋中后期煙云上邊界模糊，特別是接近頂高時更是無法表征，除1 kg TNT 能測量實驗煙云最終高度外，16 和62 kg 煙云后期均難以準確測量，只能通過仿真表征其高度。圖3 結(jié)果顯示仿真與實驗煙云高度分布基本一致，實驗結(jié)果由于風場不確定性和視野誤差出現(xiàn)少量拐點，仿真曲線平滑性和規(guī)律性更好。實驗煙云前期（爆后10 s）高度略高于其仿真值，考慮為實驗煙云擴散前期存在地面反沖和燃燒反應(yīng)，中后期實驗與仿真規(guī)律幾乎一致，平均相對誤差約5%說明仿真結(jié)果可信度高。實驗和仿真均顯示，典型風場條件下煙云高度會降低，風對煙云第2 階段影響較小，16 kg TNT 爆炸前15 s 和62 kg TNT 爆炸前25 s 無風和典型風下煙云高度差距15%內(nèi)，考慮為前期快速加速和逐步減速過程中，盡管風加速煙云與空氣混合，但加速和減速過程總路程（高度）差距較小。風場加速空氣與煙云混合導(dǎo)致第3 階段煙云高度差距擴大，煙云快速瓦解、后期上升動力不足，典型風場煙云上升時間和最終高度明顯小于無風煙云，風場煙云將較快趨于穩(wěn)定。

圖 3 實驗與仿真煙云高度對比Fig. 3 Comparison of heights between experimental and simulated clouds

3.2 不同風速下煙云擴散高度模型

3.2.1 煙云擴散過程中高度變化模型

不同風速下的煙云高度分布數(shù)據(jù)如圖4～5 所示，可以看出，不同風速下煙云高度變化呈前快后慢的冪函數(shù)增長規(guī)律，遵循浮力加速與重力、大氣阻力減速物理模型。采用Matlab 冪函數(shù)數(shù)學模型H(t)=atb進行高度變化規(guī)律擬合（H 單位為m，t 單位為s），同公斤級下初始煙云H(1)高度相同[12]，a=H(1)=(6.3±1)M(0.29±0.03)，計算所得b 值如圖6 所示，發(fā)現(xiàn)b 值與風速呈線性關(guān)系，擬合后得到b=(0.5±0.003)?(0.024±0.001)v，擬合度R 為99%，故煙云高度隨時間變化模型H(t)=(6.3±1)M(0.29±0.03)t(0.5±0.003)?(0.024±0.001)v。該擴散規(guī)律可用于非現(xiàn)場方式反演爆炸當量并推測污染煙云高度，直接實現(xiàn)爆炸危害快速偵測與預(yù)警，也可結(jié)合其他現(xiàn)場痕跡計算方法精確反演爆炸當量[1]。

圖 4 16 kg TNT 爆炸的煙云高度變化Fig. 4 Change of cloud heights for the explosion of 16 kg TNT

圖 5 62 TNT 爆炸的煙云高度變化Fig. 5 Change of cloud heights for the explosion of 62 kg TNT

圖 6 高度變化參數(shù)(b)擬合直線Fig. 6 Fitting line of heights variation parameter b

3.2.2 煙云擴散最終高度計算模型

不同風速下的爆炸煙云最終高度數(shù)據(jù)如圖7～8 所示，可以看出最終高度與風速大小呈線性下降趨勢，高度與當量整體符合Hmax=cMd模型（Hmax單位為m，M 單位為kg）[12]，由實驗可得當M=1 時，c 值區(qū)間為33～36 m，高度模型可轉(zhuǎn)變?yōu)镠max(M,v)=(34.5±1.5)Md(v)，取33、34.5、36 分別代入16 和62 kg 最大高度下求解d 值，所得d 值如圖9 所示，可以看出d 值隨風速呈線性下降，擬合得到d(v)=(0.47±0.01)?(0.038±0.002)v，擬合度R 為92.6%，Hmax(M,v)=(34.5±1.5)M(0.47±0.01)?(0.038±0.002)v。該爆高計算模型與Church 公式相比，當v=0 時該模型在計算小于50 kg 當量爆高時準確度更高，在50 kg 以上與Church 計算結(jié)果差距不大；當v≠0 時，Church 公式無法準確計算爆高，而該模型則可用于不同風速下爆高計算。復(fù)雜擴散條件下煙云爆高計算模型的建立可使爆高計算誤差更小，進而提高煙云污染評估結(jié)果準確度。

圖 7 16 kg TNT 爆炸的煙云最終高度Fig. 7 Final cloud height for 16 kg TNT explosion

圖 8 62 kg TNT 煙云最終高度Fig. 8 Final cloud height for 62 kg TNT explosion

3.2.3 最終高度計算模型使用誤差分析

煙云最終高度模型誤差主要有三個來源：實驗估計煙云最終高度誤差，線性擬合R 值過大和煙云擴散計算模型適應(yīng)條件。單純實驗估計最終測量高度定量誤差在10%～15%[4]，由于仿真計算能清晰顯示煙云頂高邊界，極大降低了這一誤差。模型線性擬合精度通過擬合參數(shù)R 評估，高度變化模型擬合度R＞0.99，最終高度模型擬合度R＞0.92，最終高度模型誤差在10%內(nèi)，該精度用于煙云污染評估可接受。最終高度模型忽略了較小地面因素及爆轟差異影響，計算模型適應(yīng)條件為穩(wěn)定大氣，相比不穩(wěn)定大氣下模型煙云高度整體偏低[4]，該模型計算結(jié)果可認為是誤差較小的爆炸煙云擴散高度下限值，煙云越低其計算所得危害結(jié)果卻代表近場最大危害結(jié)果[13]，在煙云污染應(yīng)急評估中采用高度下限值劃分的應(yīng)急區(qū)域是面積較大的保守結(jié)果，符合事故初期源項不確定情況下安全評價與應(yīng)急救援基本原則。

圖 9 最終高度參數(shù)d 的擬合直線Fig. 9 Fitting line of final height parameter (d)

3.3 不同風速下煙云擴散參數(shù)變化規(guī)律

針對前期不同風速下出現(xiàn)的高度變化差異，圖10 和11 提取了煙云擴散參數(shù)云圖中的最高溫度、最小密度和最大速度進行數(shù)值分析。煙云中最高溫度出現(xiàn)在左右渦環(huán)中心處，數(shù)據(jù)顯示溫度隨時間呈反比例衰減，熱釋放過程大多在5～10 s，當量越大下降至常溫的速度越慢，風速越高下降至常溫越快。煙云密度分布呈現(xiàn)從渦環(huán)中心到外環(huán)逐步增大的規(guī)律，實驗中也觀測到煙云濃度外圍大于中心處現(xiàn)象，煙云最小密度出現(xiàn)在渦環(huán)中心處，變化呈現(xiàn)前期快后期慢的指數(shù)函數(shù)增長規(guī)律，風速越大密度越快趨近空氣密度，無風條件下16 和62 kg 煙云最小密度增至1 kg/m3需2 和30 s，而風場條件下僅10 和12 s，這是導(dǎo)致煙云高度差距的主要原因。煙云上升過程中最大速度并非上升速度而是渦環(huán)中心處的漩渦速度，3 秒左右煙云渦環(huán)翻滾至最大速度之后逐步衰減，風速越大衰減速度越快，因其先加速后衰減的特性導(dǎo)致煙云前期盡管與空氣快速混合但高度差距較小，中后期速度和密度已接近空氣導(dǎo)致其提前進入穩(wěn)定狀態(tài)，煙云最終高度也因此出現(xiàn)較大差距。

圖 10 16 kg TNT 爆炸煙云擴散參數(shù)變化Fig. 10 Change of cloud diffusion parameters for 16 kg TNT explosions

圖 11 62 kg TNT 爆炸煙云擴散參數(shù)變化Fig. 11 Change of cloud diffusion parameters for 62 kg TNT explosions

4 結(jié) 論

(1) 爆炸煙云仿真與實驗結(jié)果顯示煙云時空分布過程與參數(shù)幾乎一致，說明基于CFD 建模的煙云擴散仿真方法可較好表征煙云擴散過程、形態(tài)及物理參數(shù)。

(2) 大氣穩(wěn)定不同風速條件下煙云可見階段高度分布擬合函數(shù)為H(t)=(6.3±1.0)M(0.29±0.03)t(0.5±0.003)?(0.024±0.001)v，該規(guī)律可初步用于實現(xiàn)非現(xiàn)場方式快速反演爆炸當量并推測煙云污染最終高度。風場下煙云擴散最終高度理論模型為Hmax(M,v)=(34.5±1.5)M(0.47±0.01)?(0.038±0.002)v，該模型可較準確計算復(fù)雜風場條件下事故污染煙云穩(wěn)定高度值。

(3) 爆炸煙云仿真結(jié)果顯示煙云流場由左右兩個反漩渦環(huán)流組成，溫度、速度最大值和密度最小值均出現(xiàn)在左右渦環(huán)中心處。實驗和仿真結(jié)果均顯示風將加快煙云與空氣混合速度，風速越大煙云溫度和速度衰減越快，密度加速趨于空氣。風場下煙云高度在前期熱抬升階段差距較小，但風對煙云抬升中后期瓦解作用明顯導(dǎo)致高度差距明顯，風速越大煙云抬升至穩(wěn)定的時間越短、煙云最終高度越低。