基于歐拉法的核爆煙塵顆粒群運(yùn)動模擬

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

（1.西安交通大學(xué)，西安，710025；2.北京宇航系統(tǒng)工程研究所，北京，100076）

0 引言

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

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

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

本文在二維圓柱坐標(biāo)系中基于歐拉網(wǎng)格方法建立了煙塵顆粒群運(yùn)動模型，考慮了歷史試驗(yàn)數(shù)據(jù)與真實(shí)大氣條件，對顆粒速度、質(zhì)量進(jìn)行時空表征，研究了核爆煙塵顆粒群的傳輸擴(kuò)散規(guī)律。

1 理論建模

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

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

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

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

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

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

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

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

運(yùn)用梯度輸送理論［7］，任意物理量的脈動輸送值與該特征量的平均值的梯度成比例關(guān)系。若該物理量為擴(kuò)散物理的質(zhì)量濃度c，則有：

式中Kx，Ky，Kz分別表示x，y，z方向的湍流擴(kuò)散系數(shù)。

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

不考慮風(fēng)速，認(rèn)為核煙塵顆粒群在三維空間呈軸對稱分布，即可在二維圓柱坐標(biāo)系中建立顆粒群運(yùn)動模型。

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

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

2 數(shù)值計算方法

在明確煙塵顆粒群運(yùn)動物理模型的基礎(chǔ)上，基于歐拉視角的數(shù)學(xué)求解方法按照如圖3所示的流程進(jìn)行。

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

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

圖4 顆粒群運(yùn)動模型示意Fig.4 Schematic of particle population movement model

下邊界取為地面，地面與上方存在質(zhì)量交換。上邊界與右邊界應(yīng)取適當(dāng)大的值，確保邊界處濃度為0即可。

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

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

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

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

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

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

表1 1 000 kt內(nèi)華達(dá)地面爆顆粒群分譜Tab.1 1 000 kt Nevada ground burst particle population spectra

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

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

本研究采用均勻正方形進(jìn)行空間的網(wǎng)格劃分，如圖4 所示。應(yīng)用時域有限差分法對式（6）進(jìn)行離散求解，時變項(xiàng)采用顯示格式，對流項(xiàng)采用一階迎風(fēng)格式，擴(kuò)散項(xiàng)采用二階中心差分格式［9］，如下：

簡化得到離散的時間步進(jìn)公式：

離散的左邊界條件為

離散的初始質(zhì)量為

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

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

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

針對wg和w，取向上為正方向，經(jīng)推導(dǎo)得到：

水平方向湍流擴(kuò)散系數(shù)取為［12］

垂直方向湍流擴(kuò)散系數(shù)取為

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

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

3 結(jié)果分析

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

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

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

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

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

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

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

4 結(jié)束語

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