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

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

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

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

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

1 爆炸煙云氣溶膠擴散理論與實驗

1.1 核化學爆炸煙云擴散理論

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

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

1.2 爆炸煙云時空分布實驗

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

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

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

2.1 煙云擴散流體動力學基本方程

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

質量守恒方程為：

(1)

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

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

(2)

能量守恒方程為：

(3)

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

2.2 煙云擴散仿真建模方法與計算

根據SCM對煙云擴散過程與機理的描述，采用CFD方法對煙云擴散開展數值仿真[11]。通過在Fluent求解器中編寫UDF源項，結合火球現場測量與公式R=1.93M0.32/(T/3 600)1/3(T=2 861 K)確定煙團直徑和高度，1 kg TNT爆炸將形成半橢球形態(tài)煙團，水平直徑約4.17 m、高約3.13 m。火球膨脹到最大尺寸后壓力略低于或等于空氣壓力，由ρ=ρ0(T0/T)可得煙團整體密度范圍為0.08～0.12 kg/m3。煙團組分按TNT化學爆炸式理想氣體產物體積分數設定，由于氣溶膠顆粒整體體積分數小于1/10，采用離散相模型反復創(chuàng)建出空氣動力學直徑為1、10、30、50、70、100 μm的球形顆粒來描述氣溶膠顆粒。由于氣溶膠比例較小，在外場開放空間碰撞并凝并概率低，計算中暫不考慮其對ARF和RF的影響[12]。考慮到溫度下降迅速，將煙團混合物平均溫度設為約1 000 K，結合實驗現象設置煙團的初始翻滾抬升速度約為5 m/s。

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

3 實驗與仿真結果討論

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

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

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

圖4 氣溶膠煙云擴散仿真結果Fig.4 Simulated aerosol cloud diffusion

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

3.2 放射性煙云時空分布參數分析

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

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

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

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

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

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

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

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

圖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) 煙云中氣溶膠源項參數變化規(guī)律

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

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

4 結論

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

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

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

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