基于歐拉法的核爆煙塵顆粒群運動模擬

郭思禹，鄭 偉，張 偉，郭 俊

（1.西安交通大學，西安，710025；2.北京宇航系統工程研究所，北京，100076）

0 引言

核爆煙塵是核爆炸火球熄滅后形成的一種蘑菇狀、具有放射性的煙云和塵柱，持續時間長，覆蓋范圍廣，煙塵顆粒物形成的熱力學環境會對廣域空間中的人員以及設備造成損傷［1］。

煙塵顆粒物的產生與分布涉及等離子體物理、輻射流體力學、凝聚態物理等多個學科。目前學術界尚難精確模擬核爆炸爆源初場［2］，主要研究核爆后煙塵顆粒物運動規律，以及沉降到地面后產生的輻射效應［3］，極少關注煙塵顆粒群在大氣中的傳輸擴散。

劉朝輝等［4］模擬了煙塵中單個顆粒的大氣運動，建立了一個拉格朗日模型，將顆粒分為大顆粒與小顆粒，分別給出其在煙塵中的軌跡。但是，無論是大、小顆粒的劃分方法，還是顆粒初始條件的設置，都缺乏理論依據，文獻中也沒有涉及對顆粒物質量的模擬研究。卓俊等［5］建立了一個顆粒群拉格朗日模型，模擬顆粒速度、質量，但該模型以穩定煙塵為起始時刻，不能對煙塵上升階段進行表征。鄭毅等［6］采用中尺度氣象軟件RAMS模擬了上升階段的煙塵顆粒流動，但仿真中將固定粒徑顆粒的速度設為定值，且在模型中只考慮了大氣中的放射性活度濃度，沒有研究質量濃度。

本文在二維圓柱坐標系中基于歐拉網格方法建立了煙塵顆粒群運動模型，考慮了歷史試驗數據與真實大氣條件，對顆粒速度、質量進行時空表征，研究了核爆煙塵顆粒群的傳輸擴散規律。

1 理論建模

地面爆和近地面爆的物理模型如圖1所示，爆后蘑菇狀煙塵的發展由大氣中各種固體顆粒群的運動來體現，而地面沙粒、土壤、巖石等物質的成分是構成這些顆粒的主體。這些原本在地面的顆粒由于受到爆炸后氣體的抽吸帶動進入大氣中，在煙塵穩定后顆粒群又開始沉降，沉降過程中伴隨著擴散，其中，小粒徑顆粒相比大粒徑顆粒的遷移距離更遠。

圖1 地面和近地核爆炸顆粒運動物理模型Fig.1 Physical model of particle motion in ground and nearground nuclear explosions

同氣體相一樣，在氣固兩相流中，固體相同樣滿足質量守恒定律，固體相顆粒質量濃度c（單位為kg·m-3）滿足［7］：

式中u，v，w分別為x，y，z方向上的顆粒平均移動速度。

如果核煙塵中只有均勻的氣相流速，即假設層流條件下，顆粒群在氣流中的運動如圖2a 所示，事實上，大氣中存在著劇烈的湍流運動，使顆粒群與空氣之間強烈地混合和交換，使得顆粒不僅沿氣流流動方向傳輸，而且還會向各個方向擴散，如圖2b所示。

圖2 層流和湍流情況下顆粒群的擴散Fig.2 Diffusion of particle population in laminar and turbulent flow conditions

考慮由湍流引起的速度脈動和濃度漲落，可將速度和濃度寫成平均值與脈動值之和：

將式（2）代入式（1），并按流體力學中的雷諾平均法則取平均，經整理可得：

運用梯度輸送理論［7］，任意物理量的脈動輸送值與該特征量的平均值的梯度成比例關系。若該物理量為擴散物理的質量濃度c，則有：

式中Kx，Ky，Kz分別表示x，y，z方向的湍流擴散系數。

將式（4）代入式（3），有：

不考慮風速，認為核煙塵顆粒群在三維空間呈軸對稱分布，即可在二維圓柱坐標系中建立顆粒群運動模型。

將顆粒都近似成球形，顆粒群按粒徑劃分成多個區間，以某一中值粒徑下的顆粒濃度等效代替這一區間顆粒群的濃度，有：

式中Ci=Ci(z，r，t)為第i種粒徑顆粒時刻t時在坐標點(r，z)的質量濃度；w(z，r，t)為顆粒群層流速度；Kz(z，r，t)，Kr(z，r，t)分別為顆粒湍流擴散系數的垂直分量和水平分量。

2 數值計算方法

在明確煙塵顆粒群運動物理模型的基礎上，基于歐拉視角的數學求解方法按照如圖3所示的流程進行。

圖3 顆粒群運動模型算法流程Fig.3 Algorithm flow of particle population motion model

在二維圓柱坐標系中，將左邊界取為核煙塵的對稱軸，如圖4所示，則左邊界滿足第二類邊界條件：

圖4 顆粒群運動模型示意Fig.4 Schematic of particle population movement model

下邊界取為地面，地面與上方存在質量交換。上邊界與右邊界應取適當大的值，確保邊界處濃度為0即可。

被帶到大氣中的土壤總質量為［3］

式中kH=0.077 41；Y（kt，TNT）為核爆當量；Hˉ為爆炸比高，0 ≤Hˉ≤7.19 m/t1/3.4。

對核煙塵放射性沾染尺寸分布的測量［8］表明，煙塵粒子的數量按粒徑近似服從對數正態分布：

將粒子看成等密度球體，粒子表面積與粒徑成二次方，體積與粒徑成三次方，則粒子表面積與體積也均近似服從對數正態分布：

設定1 000 kt 內華達地面爆場景，煙塵顆粒平均密 度ρS=2 600 kg/m3，μN=ln0.407 μm，μS=2.94 μm，μV=5.61 μm，σ=ln 4 μm。

將煙塵顆粒群按粒徑分為16 個區間，各區間的初始顆粒數如表1所示。

表1 1 000 kt內華達地面爆顆粒群分譜Tab.1 1 000 kt Nevada ground burst particle population spectra

假設初始煙塵是一個以爆心為球心、與地面相接的球狀物，初始半徑為

假定煙塵顆粒群在初始煙塵中均勻分布，初始時間為

本研究采用均勻正方形進行空間的網格劃分，如圖4 所示。應用時域有限差分法對式（6）進行離散求解，時變項采用顯示格式，對流項采用一階迎風格式，擴散項采用二階中心差分格式［9］，如下：

簡化得到離散的時間步進公式：

離散的左邊界條件為

離散的初始質量為

在軸向上，顆粒會受到重力、浮力和流體的阻力等，當顆粒速度小于周圍流體速度時，阻力表現為曳力。

忽略浮力及其他力，將w取為重力與阻力達到平衡時的值［10］：

式中wg為核煙塵中氣體的發展速度與密度，根據典型核試驗數據取值［8］；CD=0.44為阻力系數；ρg為大氣背景密度［11］。

針對wg和w，取向上為正方向，經推導得到：

水平方向湍流擴散系數取為［12］

垂直方向湍流擴散系數取為

設置煙塵穩定前時間步長0.05 s，穩定后時間步長2 s，空間步長Δz=Δr=200 m，時間尺度1 h，空間尺度為z方向20.5 km、r方向12 km。

在物理上判斷解是否收斂的方法是在每次時間步進計算后，即每對式（14）進行一次求解后，判斷總顆粒質量是否不變，即質量是否守恒。

3 結果分析

基于上述數學物理方法，最終計算得到16 個粒徑顆粒群的時空分布，粒徑中值13.3 μm 顆粒群質量濃度在4個時刻的空間分布如圖5所示。

圖5 粒徑中值13.3μm顆粒群數量時空分布Fig.5 Spatial distribution of the number of particle population with a median particle size of 13.3μm

在圖5 中，標注了質量濃度大于1×10-3g/m3的顆粒群的軸向遷移速度，可以看出，由于上方大氣空氣稀薄，距地面更高處的速度大于低處的速度。

顆粒在升降過程中還會不斷擴散，在空間近似滿足高斯分布，由于在柱坐標系下建模，明顯可見顆粒群向r軸正方向遷移。

在爆后60 s，中值粒徑13.3 μm 顆粒群的質量濃度中心大概在地面上方約5 km處，爆后5 min升高到約13 km，之后緩慢降低，爆后1 h濃度中心為地面上方約8 km。此外，由于顆粒群的擴散，濃度中心處的質量濃度一直隨時間減少，爆后60 s 為5～6 g/m3，爆后1 h為0.15 g/m3。

將模擬出的爆后1～5 min的5個時刻下2.37 μm顆粒群在空間的近似高斯分布的3σ值作為煙塵高度與寬度，并將此結果與已有文獻中核試驗數據［8］作對比，結果如表2 所示。由表2 可知，本模型得到的顆粒群運動仿真結果與已有核試驗數據一致性良好。

表2 煙塵顆粒群模擬結果對比Tab.2 Comparison of simulation results of fallout particle population

4 結束語

本文基于歐拉方法構建了核爆煙塵顆粒運動模型，得到了與真實核爆炸運動過程相似的現象，模擬結果與已有核試驗數據基本一致。因此，采用該方法模擬煙塵顆粒群空間傳輸是可行的，且可節省大量時間和費用。但由于沒有采用計算流體力學方法，渦流現象難以模擬，也沒有考慮風速的影響，這兩個方面的研究有待于今后繼續開展。