核材料炸藥化學(xué)爆炸條件下煙云源項的實驗與仿真

段中山，過惠平，羅永鋒，羅昆升，尚愛國，龔朋彬

(1.火箭軍工程大學(xué)，陜西 西安 710025；2.陸軍勤務(wù)學(xué)院，重慶 401331；3.火箭軍研究院，北京 100094)

在核材料與炸藥組成的體系中，一旦炸藥發(fā)生爆炸將導(dǎo)致核材料氣溶膠化并產(chǎn)生大量放射性煙云，對周圍人員及環(huán)境造成嚴(yán)重輻射危害，常見核材料炸藥爆炸事件主要包括臟彈恐怖襲擊、核武器化學(xué)爆炸以及粗糙核裝置爆炸。目前的研究大多運用高斯煙云模式、隨機(jī)游走算法等風(fēng)場擴(kuò)散模型開展外場核材料化學(xué)爆炸事故評估，以獲取事故后放射性污染強度和范圍[1-3]。不管采用哪種方法，都需先獲取核材料參數(shù)、炸藥參數(shù)、穩(wěn)定煙云高度H、穩(wěn)定煙云半徑R、空載氣溶膠釋放率(ARF)、可吸入比率(RF)和氣溶膠可吸入釋放率(RRF，RRF=ARF×RF)等源項參數(shù)。由于外場煙云穩(wěn)定狀態(tài)大多在高空位置，H、R、ARF、RF難以實驗測量，而大多只能采用經(jīng)驗估算。當(dāng)前外場煙云高度與半徑計算大多采用Church公式(H=76M0.25、R=13.13M0.25，M為當(dāng)量)[4-5]，而ARF大多取值地面氣化率，即假設(shè)爆轟生成的氣溶膠全部進(jìn)入到穩(wěn)定煙云中，RF取值于地面粒徑測量結(jié)果[6-8]。事實上Church公式主要由53 kg以上TNT爆炸及核爆炸實驗結(jié)果擬合得到，未討論其用于計算小尺度TNT爆炸煙云高度誤差，而從外場炸藥爆炸到形成穩(wěn)定煙云這段時間已有一定比例氣溶膠顆粒就地發(fā)生了驅(qū)散或沉降，并非所有地面氣溶膠全部進(jìn)入煙云并輻射下游地區(qū)，即ARF、RF值與地面氣化率、粒徑測量值都存在明顯差距，簡單且粗糙的源項估算方法導(dǎo)致事故評價結(jié)果誤差較大。

為解決以上問題，本文擬開展煙云擴(kuò)散抬升階段外場實驗與仿真，建立一種可用于煙云擴(kuò)散抬升階段的數(shù)值仿真方法，并用煙云實驗測量值驗證仿真方法的正確性，最后在仿真結(jié)果中提取實驗測量難以獲取的煙云參數(shù)及氣溶膠參數(shù)變化規(guī)律。

1 爆炸煙云氣溶膠擴(kuò)散理論與實驗

1.1 核化學(xué)爆炸煙云擴(kuò)散理論

圖1 SCM示意圖Fig.1 Schematic of SCM

源項模型(SCM)將炸藥爆轟過程分為3個階段[9]，具體如圖1所示。第1階段是爆轟熱力學(xué)階段。爆炸瞬間，高溫、高壓的氣體產(chǎn)物迅速膨脹，同時溫度、壓力不斷下降，當(dāng)火球膨脹至地面煙團(tuán)直徑最大時刻，高溫、低密度且包含放射性氣溶膠的初始煙團(tuán)形成[6]。第2階段是爆炸煙云主動抬升階段。爆炸膨脹完成后，煙團(tuán)在大氣浮力作用下將脫離地面開始擴(kuò)散上升，同時空氣不斷進(jìn)入、密度變大、溫度降低，期間大部分氣溶膠顆粒將隨煙云擴(kuò)散上升，大粒徑顆粒由于重力作用而沉降，小粒徑顆粒由于碰撞而凝并，沉降和凝并則會導(dǎo)致ARF和RF發(fā)生變化，該階段是一個復(fù)雜的流體力學(xué)過程，可采用計算流體動力學(xué)(CFD)方法開展數(shù)值仿真。第3階段是大氣被動擴(kuò)散階段。煙團(tuán)到達(dá)頂高后基本穩(wěn)定，煙團(tuán)和煙柱中的氣溶膠隨風(fēng)場逐步擴(kuò)散、沉降并輻射下游區(qū)域。

1.2 爆炸煙云時空分布實驗

為獲取煙云擴(kuò)散時空分布過程并檢驗數(shù)值仿真結(jié)果的正確性，本文實施了TNT爆炸煙云擴(kuò)散外場實驗。實驗設(shè)計了1 kg密度為1.63 g/cm3的TNT裸裝藥爆炸裝置，裝置上圍均勻包裹少量惰性金屬粉末用于煙云示蹤，大氣穩(wěn)定度、風(fēng)速、地面條件、炸藥殼體等都將影響煙云擴(kuò)散，實驗選擇在大氣穩(wěn)定(陰天、傍晚)、無風(fēng)(v<0.5 m/s)、硬質(zhì)水泥地面條件下開展。實驗簡易流程如圖2所示，首先根據(jù)經(jīng)驗公式對煙云的擴(kuò)散范圍進(jìn)行大致估算[4]，標(biāo)定爆炸煙云擴(kuò)散視野并調(diào)整照相機(jī)視場，使煙云在可視階段的飄散范圍全部包含在視場范圍內(nèi)；然后選取理想實驗條件時刻起爆裝置，同時用照相與攝像方式記錄爆炸煙團(tuán)初始形態(tài)及動態(tài)演變過程；最后根據(jù)視場標(biāo)定結(jié)果，使用像素計算軟件對實驗視頻及圖片進(jìn)行分析，獲取爆炸煙云擴(kuò)散時空分布參數(shù)。

圖2 煙云擴(kuò)散實驗簡易流程Fig.2 Simple process of cloud diffusion experiment

2 爆炸煙云氣溶膠擴(kuò)散的理論方程與仿真計算

2.1 煙云擴(kuò)散流體動力學(xué)基本方程

根據(jù)煙云浮力擴(kuò)散理論描述，爆炸煙云源項主要受重力、上升阻力和大氣密度差帶來的浮力等影響，煙云在浮力作用下的擴(kuò)散上升過程中遵循質(zhì)量、動量和能量守恒方程。

質(zhì)量守恒方程為：

(1)

其中：ux、uy、uz為x、y、z方向的速度分量；t為時間；ρ為密度。

煙云擴(kuò)散過程中滿足x、y、z3個方向的動量方程，以x方向為例：

(2)

能量守恒方程為：

(3)

其中：E=h-p/ρ+u2/2為流體微團(tuán)總能量，包含內(nèi)能、動能和勢能；h為焓，J/kg；hj為組分j的焓，J/kg；keff=k+kt為有效熱傳導(dǎo)系數(shù)，W/(m·K)，k為流體本征熱導(dǎo)率，kt為湍流熱傳導(dǎo)系數(shù)，根據(jù)所用的湍流模型來確定；Jj為組分j的擴(kuò)散通量；Sh為化學(xué)反應(yīng)熱及定義的其他體積熱源項[10]。

2.2 煙云擴(kuò)散仿真建模方法與計算

根據(jù)SCM對煙云擴(kuò)散過程與機(jī)理的描述，采用CFD方法對煙云擴(kuò)散開展數(shù)值仿真[11]。通過在Fluent求解器中編寫UDF源項，結(jié)合火球現(xiàn)場測量與公式R=1.93M0.32/(T/3 600)1/3(T=2 861 K)確定煙團(tuán)直徑和高度，1 kg TNT爆炸將形成半橢球形態(tài)煙團(tuán)，水平直徑約4.17 m、高約3.13 m?；鹎蚺蛎浀阶畲蟪叽绾髩毫β缘陀诨虻扔诳諝鈮毫Γ搔?ρ0(T0/T)可得煙團(tuán)整體密度范圍為0.08～0.12 kg/m3。煙團(tuán)組分按TNT化學(xué)爆炸式理想氣體產(chǎn)物體積分?jǐn)?shù)設(shè)定，由于氣溶膠顆粒整體體積分?jǐn)?shù)小于1/10，采用離散相模型反復(fù)創(chuàng)建出空氣動力學(xué)直徑為1、10、30、50、70、100 μm的球形顆粒來描述氣溶膠顆粒。由于氣溶膠比例較小，在外場開放空間碰撞并凝并概率低，計算中暫不考慮其對ARF和RF的影響[12]?？紤]到溫度下降迅速，將煙團(tuán)混合物平均溫度設(shè)為約1 000 K，結(jié)合實驗現(xiàn)象設(shè)置煙團(tuán)的初始翻滾抬升速度約為5 m/s。

計算采用60 m×60 m二維空間。網(wǎng)格從中間到兩邊、從下到上以1.02倍步進(jìn)遞增方法劃分，該方法劃分的網(wǎng)格在煙云擴(kuò)散區(qū)域密集度高、計算更準(zhǔn)確，下邊界為墻面，其他邊界為自由界面。選取非定常、2階方程計算方式及標(biāo)準(zhǔn)k-ε湍流方程。對x方向和y方向速度、能量、各物相分?jǐn)?shù)及k值等進(jìn)行收斂設(shè)定和殘差監(jiān)視，計算中殘差線下行收斂并達(dá)到計算要求精度[10]。

3 實驗與仿真結(jié)果討論

3.1 放射性煙云演變形態(tài)分析

圖3a～h較完整地記錄了不同時刻實驗煙云形態(tài)，圖中煙云抬升速度由快到慢，濃度不斷降低直至穩(wěn)定狀態(tài)。圖3a～d顯示，煙云擴(kuò)散前期由于示蹤劑驅(qū)散導(dǎo)致形態(tài)不夠規(guī)整。圖3e～h顯示，后期煙云在浮力主導(dǎo)下其擴(kuò)散的整體呈現(xiàn)橢球形態(tài)，煙團(tuán)上端分布均勻、下端煙柱明顯。

圖3 實驗煙云擴(kuò)散過程Fig.3 Experimental cloud diffusion process

圖4 氣溶膠煙云擴(kuò)散仿真結(jié)果Fig.4 Simulated aerosol cloud diffusion

圖4為煙云擴(kuò)散到約25 m高度的氣溶膠分布計算結(jié)果(圖3f時刻)，圖片清晰顯示出此時氣溶膠隨煙云擴(kuò)散的整體分布與形態(tài)，由于煙云密度梯度差產(chǎn)生的浮力將氣團(tuán)由傳統(tǒng)的高斯形態(tài)濃度分布轉(zhuǎn)變成由兩個反旋渦環(huán)流組成的渦流分布，橫截面上展現(xiàn)出由兩個反旋渦環(huán)流所組成的“腎形狀”模式，即“蘑菇云”現(xiàn)象。渦環(huán)圍繞其渦核由內(nèi)向外翻卷，直徑逐漸變大且上升過程非常穩(wěn)定，仿真計算結(jié)果與實驗觀察吻合，表明氣溶膠上升驅(qū)動力主要來自于低密度高溫?zé)焾F(tuán)浮力并非完全依賴初始爆轟驅(qū)散。

3.2 放射性煙云時空分布參數(shù)分析

煙云擴(kuò)散高度、直徑(寬度)和寬高比隨時間變化的實驗和仿真結(jié)果示于圖5。

從圖5a可見，實驗煙云前期較仿真結(jié)果略高，考慮為實驗過程存在示蹤劑爆轟拋撒及地面反彈，而仿真模型則是理想煙團(tuán)浮力模型，中后期煙云變化規(guī)律基本一致，實驗和仿真結(jié)果相互驗證、數(shù)據(jù)可信度高。煙云高度隨時間變化整體呈現(xiàn)前快后慢的冪函數(shù)擴(kuò)散規(guī)律，與文獻(xiàn)[13]的結(jié)果一致。采用Church模型計算的1 kg TNT爆炸高度高達(dá)76 m，而仿真和實驗結(jié)果均顯示煙云約在30 s時刻幾乎不再上升，此時煙云穩(wěn)定高度約38 m，證明Church模型在計算小當(dāng)量爆高時誤差過大，小當(dāng)量爆炸煙云頂高可采用CFD仿真方法直接計算得到或采用小當(dāng)量實驗確定。

煙云寬度是計算煙云體積的重要依據(jù)，Church模型給出了煙云寬度約為高度0.4倍的結(jié)論。圖5b顯示，實驗和仿真煙團(tuán)寬度變化曲線規(guī)律較為一致，煙云寬度隨時間幾乎呈線性增長，但由于實驗存在拋撒導(dǎo)致測量寬度略大，最終煙云寬度約為22 m。圖5c為煙云寬度與高度數(shù)據(jù)對比，斜率顯示其寬高比先小后大，最終寬高比可達(dá)0.58，這一結(jié)果與德國爆炸實驗結(jié)果[14]吻合。

3.3 煙云中氣溶膠擴(kuò)散參數(shù)變化規(guī)律

1) 不同粒徑氣溶膠在煙云中的分布規(guī)律

由圖4可見，大部分氣溶膠顆粒隨向上運動的雙渦環(huán)集聚于爆炸煙云的頂部和莖部形成蘑菇云及煙柱，少量顆粒物將隨四周橫向運動的氣流發(fā)生驅(qū)散或由于重力作用沉降于地面。為更好地觀察煙云中顆粒物的分布規(guī)律，截取了煙云中1、10、50、100 μm顆粒物在煙云穩(wěn)定時刻的分布(圖6)。從圖6可見，此時渦環(huán)高度和寬度已達(dá)最大，組成煙云頂端的主要成分是1 μm和10 μm的顆粒物，50 μm顆粒物主要分布在煙云莖部或少量到達(dá)頂端，由于大粒徑顆粒物穩(wěn)定上升至頂端需要流場向心力最大位置才能完成，故卷入頂端的少量顆粒物主要位于外環(huán)，實驗后期也顯示煙云外環(huán)位置能見度低于中心位置，也證明外環(huán)處顆粒物濃度高或較多大粒徑顆粒物匯聚在此。100 μm顆粒物幾乎不能到達(dá)渦環(huán)頂高，大多集聚在煙柱中或沉積于地面。

2) 不同粒徑氣溶膠的擴(kuò)散沉降規(guī)律

煙云在擴(kuò)散過程中顆粒物將發(fā)生水平驅(qū)散和垂直沉降[15]，根據(jù)氣溶膠在100 s時間內(nèi)變化的統(tǒng)計結(jié)果(圖7)，可將其運動過程分為3個階段：第1階段發(fā)生在前10 s，煙云抬升時少量顆粒物隨周圍小渦環(huán)發(fā)生水平驅(qū)散或未卷入渦環(huán)，這些顆粒物在重力作用下迅速沉降于地面并造成煙云中顆粒物數(shù)目快速下降，顆粒物粒徑越大，沉降越快，比例越高，如100 μm顆粒物4 s時減少約40%，而可吸入粒徑顆粒則90%以上隨煙云上升；第2階段為顆粒物停留階段，顆粒物隨煙云上升過程中將經(jīng)歷浮力加速、重力與阻力減速、逐步從煙頂分開經(jīng)煙柱在空氣中懸浮、沉降過程，該階段顆粒物數(shù)明顯進(jìn)入一個停留期，粒徑越大停留期越短，如100 μm顆粒物在煙云中停留約40 s之后開始沉積于地面，而70 μm顆粒物在80 s時刻開始沉積于地面，其他粒徑顆粒物停留期均大于100 s；第3階段為懸浮、沉積階段，經(jīng)過不同停留期后將呈現(xiàn)大顆粒沉積、小顆粒懸浮現(xiàn)象，該規(guī)律與文獻(xiàn)[9]經(jīng)驗計算結(jié)果基本吻合。

圖5 煙云高度、直徑和寬高比對比Fig.5 Comparison of cloud height, diameter and diameter and height ratio

圖6 穩(wěn)定時刻不同粒徑的顆粒物分布Fig.6 Particle distribution of different sizes at steady time

圖7 顆粒物在煙云中停留時間Fig.7 Residence time of particulate matter in cloud

3) 煙云中氣溶膠源項參數(shù)變化規(guī)律

圖8 煙云抬升過程中顆粒物的留存比Fig.8 Retention ratio of particles during rising of cloud

圖8為煙云抬升過程中顆粒物在空氣中的留存比(RR)。煙云主動擴(kuò)散到穩(wěn)定狀態(tài)，由于氣溶膠運動造成其ARF和RF與地面值相比發(fā)生了變化，1 μm和10 μm顆粒物在空氣中的留存比約90%，大部分分布在煙云上部渦環(huán)中；30 μm和50 μm顆粒物分別有63%和52%在空氣中，大多位于頂部煙云的煙柱或外圍；70 μm和100 μm顆粒物驅(qū)散或沉降較多，分別只有46%和36%留在空氣中。以O(shè)peration Roller Coaster(ORC)實驗質(zhì)量歸一化分布為例，地面實驗測量ARF為100%、RF為20%、RRF為20%，粒徑區(qū)間為<10、10～30、30～50、50～70、70～100、>100 μm時的質(zhì)量比(MR)分別為20%、25%、15%、10%、10%、20%[5]。煙云中ARF計算公式為ARF=∑RR×MR，RR值取粒徑留存比區(qū)間期望值，計算得到ARF約58%，由于ARF變化導(dǎo)致煙云中可吸入比例升高，RF由20%增至31%，RRF由20%變?yōu)?8%。這一結(jié)果說明煙云抬升階段顆粒物驅(qū)散與沉降會導(dǎo)致爆炸點附近輻射強度增加，進(jìn)入大氣擴(kuò)散模式的氣溶膠數(shù)量整體減少，但可吸入部分氣溶膠數(shù)量幾乎保持不變。

4 結(jié)論

本文通過對核材料化學(xué)爆炸煙云的實驗與仿真分析，得到以下結(jié)論。

1) CFD方法仿真煙云擴(kuò)散的過程、形態(tài)和時空分布參數(shù)與實驗吻合。該方法可用于核材料化學(xué)爆炸煙云擴(kuò)散實驗設(shè)計、源項計算、近場污染評估等。

2) 煙云流場由兩個反旋渦環(huán)流組成，橫截面展現(xiàn)出一個“腎形狀”模型，煙云中可吸入粒徑氣溶膠大部分分布在煙云渦環(huán)中且顆粒粒徑越小離渦環(huán)中心越近，粒徑在50 μm以上顆粒在重力作用下沉降于地面或集聚在煙柱中。1 kg TNT爆炸煙云高度約38 m、寬度約22 m，最終寬高比可達(dá)0.58。

3) 由于驅(qū)散與沉降效應(yīng)煙云中氣溶膠大多經(jīng)過初期減少、顆粒物停留、逐步沉降3個階段，沉降導(dǎo)致穩(wěn)定煙云ARF與RF均發(fā)生一定變化，以O(shè)RC實驗質(zhì)量歸一化分布為例，其煙云穩(wěn)定時刻ARF約為地面氣化率的58%，RF由20%增至31%，RRF約為18%。