朱金輝,左應(yīng)紅,劉 利,牛勝利,商 鵬,李夏至,王學(xué)棟

(西北核技術(shù)研究所,西安 710024)

核輻射是核爆炸特有的殺傷因素。核爆炸的裂變核反應(yīng)和聚變核反應(yīng)過程會產(chǎn)生多種具有強輻射性質(zhì)的輻射源,包括瞬發(fā)核輻射粒子、緩發(fā)核輻射粒子及裂變產(chǎn)物衰變釋放的射線粒子[1-2]。核爆炸輻射源在爆心周圍大氣或介質(zhì)中傳輸,與大氣或介質(zhì)相互作用形成核輻射場,可采用具有一定時空分布的輻射環(huán)境參數(shù)進行描述,如輻射劑量、粒子注量、能譜、時間譜及角度譜等[3-6]。不同爆炸場景的核輻射環(huán)境特征差異較大,大氣層核爆炸與高空核爆炸是2類典型爆炸場景[7-9]。大氣層核爆炸場景中,按時間維度可將輻射環(huán)境分為持續(xù)時間較短的瞬發(fā)核輻射[10-11]、緩發(fā)核輻射及持續(xù)時間較長的放射性沾染[13-14]。高空核爆炸核輻射環(huán)境,可分為瞬時輻射環(huán)境與長期輻射環(huán)境,瞬時輻射環(huán)境包括中子、γ射線、X射線及核電磁脈沖[15]等,長期輻射環(huán)境包括碎片云、附加電離層和人工輻射帶等[16-17]。

核爆炸輻射環(huán)境場的特征取決于核輻射粒子在介質(zhì)中的輸運過程,即輻射粒子與物質(zhì)發(fā)生相互作用的過程。研究核爆炸輻射環(huán)境最常用的方法是模擬粒子輸運的蒙特卡羅方法[18-20]。蒙特卡羅方法,又稱隨機模擬方法,是通過直接跟蹤模擬大量粒子的運動軌跡,利用隨機數(shù)模擬粒子輸運中的隨機過程,抽取碰撞時的靶核、碰撞類型及出射角等參數(shù),統(tǒng)計粒子運動過程中對待求物理量的貢獻,獲得統(tǒng)計性結(jié)果。采用蒙特卡羅方法可實現(xiàn)核爆炸輻射環(huán)境的輸運計算。

本文作者團隊長期從事核爆炸輻射環(huán)境的理論與數(shù)值模擬研究,運用蒙特卡羅方法解決了核爆炸輻射環(huán)境模擬中遇到的難點技術(shù)問題,發(fā)展了一系列數(shù)值計算模型和方法,在大氣層和高空核爆炸輻射環(huán)境方面都取得了相應(yīng)進展。本文將詳細介紹蒙特卡羅方法在大氣層核爆炸復(fù)雜條件下的瞬發(fā)和緩發(fā)輻射環(huán)境高效高精度數(shù)值模擬方面,高空核爆炸X射線、附加電離層、人工輻射帶和核電磁脈沖等輻射環(huán)境數(shù)值模擬方面的應(yīng)用和發(fā)展情況。

1 蒙特卡羅方法在大氣層核輻射環(huán)境模擬中的應(yīng)用與發(fā)展

大氣層核輻射環(huán)境模擬計算需解決不同輻射源在大范圍近地面環(huán)境中的模擬問題,涉及復(fù)雜條件和深穿透等問題。在長期實踐中,本文作者團隊針對相關(guān)具體問題發(fā)展了一系列蒙特卡羅減方差方法和建模方法,在大氣層核輻射模擬中得到了較好的應(yīng)用。

1.1 用于大氣層核輻射深穿透問題的權(quán)窗逼近方法

核爆炸產(chǎn)生的早期核輻射環(huán)境在地面附近大氣中影響范圍可達3～5 km,而核爆炸產(chǎn)生的中子、γ射線在近地面大氣中的自由程為百米量級。大空間范圍內(nèi)的早期核輻射輸運問題屬于典型的射線深穿透問題[1],對于此類問題,采用蒙特卡羅直接模擬方法難以獲得可信的輻射環(huán)境參數(shù)。為在有限的時間內(nèi)獲得預(yù)定的計算精度,建立了柵元權(quán)窗結(jié)合密度逼近迭代的減方差方法。首先建立γ射線長距離輸運分層模型,并合理設(shè)置計算區(qū)域幾何網(wǎng)格。圖1為γ射線近地面長距離輸運示意圖。通過不同密度條件下的多次迭代計算,利用非深穿透輸運計算得到的權(quán)窗逐步生成更優(yōu)的適合深穿透輸運計算的權(quán)窗。利用權(quán)窗重要性函數(shù)引導(dǎo)粒子的分裂與輪盤賭,并采用環(huán)探測器記錄,實現(xiàn)了γ射線在低空稠密大氣中(含地面反照)長距離輸運的蒙特卡羅模擬。

圖1 γ射線近地面長距離輸運示意圖Fig.1 Diagram of long distance γ-ray transport near ground

采用上述柵元權(quán)窗結(jié)合密度逼近迭代減方差方法,計算了離地面400 m的γ射線在爆心投影點周圍5 km范圍內(nèi)輸運形成的地面γ射線輻射場,通過多次密度迭代,最終給出實際密度條件下的計算結(jié)果。柵元權(quán)窗減方差方法,γ注量Φγ及相對偏差σr隨測點距離d的變化關(guān)系如圖2所示。5 km測點的注量、相對偏差及FOM因子F隨模擬粒子數(shù)Np的變化關(guān)系如圖3所示。由圖2和圖3可見,該減方差方法模擬得到γ射線注量的相對偏差小于5%,反映計算效率的FOM因子隨模擬粒子數(shù)的增長而保持穩(wěn)定,表明該減方差方法可有效解決γ射線長距離輸運的模擬問題。

上述方法的核心是通過迭代計算得到合適的權(quán)窗參數(shù),對中子在大范圍空間內(nèi)輸運的深穿透問題,該方法同樣適用。對于介質(zhì)密度變化不大的模型,可采用網(wǎng)格權(quán)窗替代柵元權(quán)窗以降低建模的復(fù)雜性,同時應(yīng)用密度逼近迭代方法獲取優(yōu)化的網(wǎng)格權(quán)窗參數(shù),進而對中子、γ射線輸運進行引導(dǎo)。

圖2 柵元權(quán)窗減方差方法,γ注量及相對偏差隨測點距的變化關(guān)系Fig.2 Φγ and σr vs. d by using cell-based weight window variance reduction method

(a) Φγ vs. Np

(b) σr and F vs. Np

1.2 傳輸介質(zhì)變化條件下緩發(fā)γ輻射輸運自適應(yīng)建模方法

核爆炸早期核輻射中的緩發(fā)輻射源是隨時間變化的復(fù)雜體源,最主要成分為裂變碎片衰變釋放的γ輻射。緩發(fā)γ輻射在傳輸中會受沖擊波擴展過程的影響,沖擊波會使爆心周圍一定范圍內(nèi)的空氣密度分布發(fā)生變化,導(dǎo)致緩發(fā)γ輻射劑量出現(xiàn)流體力學(xué)增強效應(yīng)[5]。

為計算流體力學(xué)增強效應(yīng)對γ射線輸運的影響,建立了沖擊波導(dǎo)致傳輸介質(zhì)變化條件下的緩發(fā)γ輻射大氣輸運自適應(yīng)建模方法。在沖擊波在特定時刻形成的大氣密度場分布條件下,現(xiàn)有γ輻射大氣輸運幾何模型通常劃分為多個均勻步長網(wǎng)格,難以反映波陣面處大氣密度的劇烈變化。同時在大氣密度變化平緩區(qū)域劃分了過多的網(wǎng)格,降低了模型描述的精確性和增加了幾何柵元。針對上述問題,利用自適應(yīng)方法將均勻步長網(wǎng)格調(diào)整為非均勻步長網(wǎng)格,即針對大氣密度變化較大的區(qū)域增加了網(wǎng)格數(shù)目,避免了等步長網(wǎng)格中因步長設(shè)置過大而導(dǎo)致的幾何模型精確性降低;針對大氣密度變化不大的區(qū)域減少了網(wǎng)格數(shù)目,避免了等步長網(wǎng)格中過多的幾何柵元而導(dǎo)致計算量的增加。

針對10 kt當(dāng)量爆炸沖擊波在爆后0.2 s時的大氣密度分布函數(shù),取步長為64 m和步長4 m分別構(gòu)建均勻粗網(wǎng)格模型和細網(wǎng)格模型。同時以64 m為默認步長構(gòu)建自適應(yīng)網(wǎng)格模型,建立一個半徑為2 km的1維球?qū)ΨQ幾何模型,大氣為理想均勻大氣,密度為1.225 kg·m-3。1維球?qū)ΨQ幾何模型的3種網(wǎng)格劃分方法如圖4所示。

圖4 1維球?qū)ΨQ幾何模型的3種網(wǎng)格劃分方法Fig.4 Three grid partitioning methods for 1D spherically symmetric geometric models

采用蒙特卡羅模擬程序開展上述1維球?qū)ΨQ幾何模型的γ輻射輸運計算,并沿徑向在大氣密度變化較大的150～350 m范圍內(nèi)均勻布放21個探測點。模擬粒子數(shù)設(shè)置為1.0×107。建立的γ輻射1維球?qū)ΨQ幾何模型柵元數(shù)和對應(yīng)的模擬計算時間如表1所列。

表1 1維球?qū)ΨQ幾何模型柵元數(shù)和對應(yīng)的模擬計算時間Tab.1 Number of grid elements and simulation time for one-dimensional spherically symmetric geometric model

以步長為4 m的均勻細網(wǎng)格模型為基準,計算得到自適應(yīng)網(wǎng)格模型和均勻粗網(wǎng)格模型與基準值(均勻細網(wǎng)格模型)相對偏差隨距離的變化關(guān)系,如圖5所示。由圖5可見,在計算結(jié)果均收斂的情況下,自適應(yīng)網(wǎng)格模型模擬結(jié)果與基準結(jié)果基本重合,平均相對偏差為0.13%,而均勻粗網(wǎng)格模計算結(jié)果與基準結(jié)果則出現(xiàn)了約為5%的相對偏差,平均相對偏差為2.9%,但自適應(yīng)網(wǎng)格模型模擬計算時間較步長為4 m的均勻精細網(wǎng)格基準模型減少了52.0%。

該方法能有效減少幾何柵元個數(shù),節(jié)約大量的內(nèi)存空間和計算時間。該方法采用自適應(yīng)網(wǎng)格模型,在保證計算精確性的前提下,能有效提升模型計算效率,是一種既能精確描述大氣密度分布情況,又盡量具有較少幾何柵元的精細高效建模方法。該方法還可通過在γ輻射輸運圓柱體模型的r-z方向進行2維自適應(yīng)建模,適用于沖擊波影響下的非均勻大氣-土壤模型。

圖5 自適應(yīng)網(wǎng)格模型和均勻粗網(wǎng)格模型與基準值(均勻細網(wǎng)格模型)相對偏差隨距離的變化關(guān)系Fig.5 Relative deviation of adaptive mesh model, and uniform coarse mesh model compared with the datum value (uniform fine mesh model) vs. distance

1.3 大規(guī)模復(fù)雜城市建模與基于網(wǎng)格的權(quán)窗減方差方法

當(dāng)核爆炸發(fā)生于城市場景,爆炸產(chǎn)生的中子、γ射線在城市中輸運時,除受空氣的散射外,還會受城市建筑物的散射、反射與衰減等復(fù)雜過程。準確計算城市大氣層核爆炸產(chǎn)生的瞬發(fā)核輻射環(huán)境參數(shù),主要面臨2個難題:一是如何將大規(guī)模城市建筑物復(fù)雜幾何批量轉(zhuǎn)換為蒙特卡羅粒子輸運計算模型;二是如何解決多棟建筑物厚屏蔽體及地表大氣中數(shù)千米遠距離輸運的深穿透模擬問題。針對這2個難題,發(fā)展了基于GIS數(shù)據(jù)的自動化輻射輸運城市幾何建模技術(shù)及基于網(wǎng)格權(quán)窗的減方差模擬技術(shù)。

城市復(fù)雜條件下的輻射建?？煞譃?個步驟:(1)采用ArcGIS、freeCAD等進行后處理開發(fā),將描述城市建筑物矢量輪廓的shapefile數(shù)據(jù),自動轉(zhuǎn)換為.step格式的3維幾何模型;(2)在3維CAD模型中,將城市建筑的墻體及樓層板設(shè)置為一定厚度,如墻體厚為30 cm,樓層板厚為10 cm,目前的模型中忽略了窗戶、門框及房間內(nèi)部隔墻等細節(jié);(3)采用SuperMC軟件將.step模型轉(zhuǎn)換為可供蒙特卡羅計算的輸入卡,該過程中需進行幾何交疊修復(fù)等模型預(yù)處理,并對地面混凝土幾何進行建模,設(shè)置建筑材料、重要性及源等必要的計算參數(shù);(4)手動完善蒙特卡羅計算的輸入卡,如設(shè)置空氣區(qū)域及空氣密度,減方差設(shè)置等。

由于城市幾何模型的尺度大、建筑物多、記錄范圍廣,會出現(xiàn)深穿透現(xiàn)象,直接蒙特卡羅模擬時,離源較遠區(qū)域位置的模擬結(jié)果相對偏差很大,因此基于網(wǎng)格的權(quán)窗,發(fā)展了針對城市核輻射粒子輸運的減方差技術(shù)。該減方差技術(shù)主要有4個步驟:(1)設(shè)置一個用于產(chǎn)生權(quán)窗的環(huán)形柵元能量沉積計數(shù):距離地面為1 m,高為5 m的靠近城市邊界的包圍整個2 km × 2 km模擬區(qū)的圓環(huán);(2)模擬區(qū)域用mesh劃分為同軸圓柱殼虛擬網(wǎng)格,柱面間距為40 m,圓柱中心軸垂直穿過爆心;(3)再進行直接蒙特卡羅模擬,依據(jù)各虛擬網(wǎng)格中粒子對環(huán)形柵元能量沉積計數(shù)結(jié)果的貢獻計算得到權(quán)窗,必要時對權(quán)窗進行迭代優(yōu)化;(4)通過權(quán)窗函數(shù)引導(dǎo)射線粒子輸運至離爆心較遠的區(qū)域,實現(xiàn)減方差模擬計算。圖6為建立的城市條件下輻射輸運幾何模型。

針對某城市模型,當(dāng)量為1 kt威力條件下,采用基于網(wǎng)格權(quán)窗減方差技術(shù),模擬得到圖7為典型城市大氣層核爆炸瞬發(fā)γ劑量場計算結(jié)果,如圖7所示。模擬粒子數(shù)為109,在未采用減方差技術(shù)時,相對偏差很大,采用減方差技術(shù)后,劑量計算結(jié)果的相對偏差明顯降低,滿足應(yīng)用要求精度。

(a) Model geometry transformed by SuperMC

(b) Model used in MCNP input

(a) No variance reduction

(b) Weight window based variance reduction

1.4 基于伴隨輸運的核輻射環(huán)境場模擬方法

在大氣層核輻射環(huán)境模擬中,常需計算不同高度輻射源在地面附近產(chǎn)生的輻射場參數(shù)。采用常規(guī)正向輸運模擬方法,需進行多次反復(fù)模擬計算。通過引入伴隨輸運模擬方法[21],可將求解一系列“單源-多探測器”的正向輸運問題轉(zhuǎn)換為求解一個“單探測器-多源”伴隨輸運問題,從而實現(xiàn)一次計算給出不同輻射源項條件下地面附近不同測點的輻射場參數(shù)。核輻射正向與伴隨輸運模型如圖8所示。

(a) Forward transportation model

(b) Adjoint transportation model

伴隨輸運是正向輸運的逆過程。在伴隨輸運方式下,從探測位置處抽取的伴隨中子或γ首先被“賦予”某一權(quán)重,經(jīng)與大氣和地面相互作用的逆向輸運模擬后,得到不同爆高和不同相對距離處的中子或γ的伴隨注量。伴隨輸運給出的伴隨注量具有反映正向源位置處中子或γ重要性參數(shù)的物理意義,即一個源中子或γ對探測器響應(yīng)的貢獻[22]。

實際計算中,地面伴隨源項對不同高度H,不同距離S處的伴隨注量對應(yīng)于高度為H的源在距地面S處產(chǎn)生的輻射注量。

在伴隨模型中,所有屏蔽介質(zhì)的幾何及材料屬性均不發(fā)生變化,與正向模型保持一致,需改變的是源與探測器的相對位置及特性,在伴隨模型中,伴隨源位于正向模型中探測器的位置,即地面附近各向同性點源,伴隨源的能量按照響應(yīng)函數(shù)進行抽樣;而探測器則位于正向模型中源的位置,針對具體問題可根據(jù)模擬需要設(shè)置一系列探測器。圖9為一次伴隨輸運與多次正向輸運計算結(jié)果的對比。由圖9可見,二者吻合良好。從兩種輸運模式的計算時效性來看,無論是正向輸運還是伴隨輸運,對于千米量級的長距離輸運蒙特卡羅計算,都需使用減方差方法來加速問題收斂,且所需的計算時長也相當(dāng)。需要說明的是,伴隨計算目前只能采用多群方式進行模擬,計算精度與多群截面數(shù)據(jù)的精度相關(guān)。在上述算例中,伴隨計算與正向輸運的相對偏差大多小于15%,滿足工程應(yīng)用需求。由圖9對比結(jié)果可見,利用蒙特卡羅伴隨方法計算瞬發(fā)輻射環(huán)境參數(shù),工程應(yīng)用價值十分顯著。

(a) Different height with the same source

(b) Different source with the same height

1.5 大空間輻射場計算的全局減方差方法

在早期核輻射和近地面電磁脈沖源項計算中,往往需計算全空間輻射場。第1.1節(jié)中采用的局域減方差方法(local variance reduction, LVR)通常只針對特定目標探測器,難以給出理想的全空間輻射場。與之相比,全局減方差方法(global variance reduction, GVR) 通過全局權(quán)窗來控制全空間模擬粒子的分布,整體提高全空間計算效率,更適合全空間輻射場的求解計算。

在大氣層核輻射環(huán)境模擬直接采用現(xiàn)有的全局減方差方法面臨兩個問題:一是現(xiàn)有全局減方差方法適用于等體積的柵元/網(wǎng)格,當(dāng)計數(shù)柵元/網(wǎng)格體積差異較大時,體積大的柵元/網(wǎng)格內(nèi)模擬粒子數(shù)較多,體積小的柵元/網(wǎng)格內(nèi)模擬粒子數(shù)較少,使全局計算效率降低;二是計算邊界附近的非計數(shù)區(qū)占用的計算資源較多,輻射場計算模型一般可劃分為計數(shù)區(qū)和非計數(shù)區(qū),非計數(shù)區(qū)域內(nèi)模擬粒子可輸運至計數(shù)區(qū)并影響其計數(shù),將全局權(quán)窗參數(shù)統(tǒng)一應(yīng)用于計數(shù)區(qū)與非計數(shù)區(qū),非計數(shù)區(qū)內(nèi)(尤其是土壤等吸收截面較大的介質(zhì)內(nèi))模擬粒子過度分裂,占用較多計算資源。因此,本文提出基于體積和非計數(shù)區(qū)修正的全局減方差方法,能有效提高全空間輻射場計算效率。

基于粒子通量的GVR方法的權(quán)窗下限可表示為

(1)

其中:Φi為第i個柵元/網(wǎng)格內(nèi)粒子的平均注量;Max(Φi)為所有柵元/網(wǎng)格中最大粒子平均注量;β為柵元/網(wǎng)格權(quán)窗上限下限之比。

引入體積修正后的權(quán)窗下限可表示為

(2)

其中:Vi為第i個柵元/網(wǎng)格的體積;Vs為點源所在的柵元/網(wǎng)格的體積。

非計數(shù)區(qū)內(nèi)模擬粒子只有輸運至計數(shù)區(qū)才能對計數(shù)區(qū)粒子注量有貢獻。非計數(shù)區(qū)柵元/網(wǎng)格內(nèi)粒子輸運至計數(shù)區(qū)柵元/網(wǎng)格的概率為e-s/λ(s為非計數(shù)區(qū)柵元/網(wǎng)格到計數(shù)區(qū)柵元/網(wǎng)格的距離,λ為粒子的平均吸收自由程)。引入非計數(shù)區(qū)修正后的非計數(shù)區(qū)的權(quán)窗下限可表示為

wth,nonc=wth,nes/λ

(3)

其中:wth,nonc為非計數(shù)區(qū)柵元/網(wǎng)格的權(quán)窗下限;wth,n為離非計數(shù)區(qū)最近的計數(shù)區(qū)柵元/網(wǎng)格的權(quán)窗下限。

圖10為全局減方差方法流程圖。具體步驟為:(1)首先建立輻射場輸運計算模型,采用蒙特卡羅方法跟蹤輻射場中模擬粒子的輸運過程,并統(tǒng)計計算所有柵元/網(wǎng)格內(nèi)粒子信息;(2)根據(jù)粒子信息計算全局減方差方法的權(quán)窗參數(shù),對權(quán)窗參數(shù)進行體積修正和非計數(shù)區(qū)修正;(3)采用修正后的權(quán)窗參數(shù),跟蹤模擬粒子在權(quán)窗引導(dǎo)下的輸運過程,并統(tǒng)計計算所有柵元/網(wǎng)格內(nèi)粒子信息;(4)判斷計算得到粒子信息的相對偏差是否收斂,如未收斂,則返回步驟(2)迭代,直到收斂。

采用直接模擬和修正后的GVR方法模擬計算了爆高為20 m,高度為0～1 km,水平距離為5 km范圍內(nèi)γ射線注量。兩次模擬的CPU總運行時間都設(shè)置為400 min。圖11直接模擬和修正后的GVR方法計算γ注量相對偏差的分布。由圖11可見,直接模擬得到大空間輻射場的平均相對偏差為20%,最大為100%;修正后的GVR方法模擬迭代3次得到的大空間輻射場平均相對偏差為1.2%,最大為4.2%,與直接模擬相比,平均相對偏差降低一個量級以上。

圖10 全局減方差流程圖Fig.10 Flow chart of the corrected GRV method

(a) Direct simulation

(b) GVR

2 蒙特卡羅方法在高空核爆炸輻射環(huán)境模擬中的應(yīng)用與發(fā)展

按照時間尺度,可將高空核爆炸產(chǎn)生的輻射環(huán)境,分為瞬時輻射環(huán)境與長期輻射環(huán)境,蒙特卡羅粒子輸運方法在高空核爆炸輻射環(huán)境的計算中也有廣泛運用。

2.1 裸核區(qū)效應(yīng)

X射線是一種光子,和源自核過程的γ射線不同,X射線主要來自于核外電子的躍遷[1],因此,能量低于γ射線,但又顯著高于可見光。一般而言,X射線的波長主要集中在0.01 ～10 nm之間[23]。對應(yīng)的單光子能量約為0.145～124 keV。X射線的能量遠小于電子對效應(yīng)的反應(yīng)閾值,因此X射線和物質(zhì)的主要相互作用為光電離效應(yīng)和康普頓效應(yīng)[24]。核爆炸中的X射線主要來自于高溫等離子體的輻射復(fù)合過程。核爆炸火球的初始溫度高達數(shù)千萬開甚至更高,且處于熱平衡狀態(tài)[25]?？梢曌饕粋€黑體[16]。因此,可用一個或多個黑體譜的組合來表征核爆炸X射線的能譜。

X射線的環(huán)境參數(shù)主要包括X射線的能注量和能譜。能注量是指通過空間中某一點單位面積上的能量?？梢宰C明,對于包含散射作用的物理過程,在經(jīng)過一定距離的輸運后,已發(fā)生散射的粒子的比例會上升并成為主要成分。對于X射線來說,當(dāng)單光子能量約為30 keV時,光電吸收截面和康普頓截面近似相等。因此,對于較硬的核爆炸X射線源,輸運時散射效應(yīng)非常明顯。故X射線的輸運過程常用蒙特卡羅方法進行模擬。X射線的模擬需要注意3個問題:一是長距離輸運時的方差控制;二是垂直大氣結(jié)構(gòu)的精確表征;三是非稀薄介質(zhì)中的近場裸核區(qū)效應(yīng)。

對于長距離輸運的方差控制問題,可使用權(quán)窗和源偏倚的方法削減方差隨距離增加的趨勢。通過多次迭代,反復(fù)優(yōu)化空間中的權(quán)窗分布和源的偏倚設(shè)置,可實現(xiàn)較長距離的輸運,同時保持較低的方差水平。

對于大氣垂直結(jié)構(gòu)表征問題,引入國際標準大氣模型[26],建立了一套可考慮不同經(jīng)緯度、海拔(0～1 000 km)、季節(jié)及太陽輻射強度的大氣組分計算程序,可在蒙特卡羅輸運中較準確地描述介質(zhì)屬性。大氣組分的垂直分布如圖12所示。

在非稀薄介質(zhì)中的近場輸運中,裸核區(qū)現(xiàn)象對X射線環(huán)境參數(shù)的影響較大。對于這一問題,基于獨立原子近似[27]和殼層電離模型[28]開發(fā)了裸核區(qū)數(shù)值計算程序[29],能給出不同爆炸場景下的裸核區(qū)參數(shù)[30],包括裸核區(qū)大小、沉積在裸核區(qū)中的能量份額及裸核區(qū)邊界能譜。通過以上方法,實現(xiàn)了大氣層內(nèi)核爆炸X射線環(huán)境的模擬計算功能。爆高不同時,高空核爆炸X射線能注量ΨX模擬計算結(jié)果如圖13所示。

圖12 大氣組分的垂直分布Fig.12 Vertical distribution of atmosphere components

(a) Downward direction

(b) Upward direction

通過大量計算,得到了大氣層內(nèi)核爆炸X射線環(huán)境的一些基本規(guī)律:對于爆炸高度在50 km以上的場景,X射線在向下傳輸時存在一個截止高度,對于不同能量的X射線,這個值為20～50 km,X射線在這一高度層幾乎全部沉積,無法向下穿透;對于爆炸高度在50 km以下的場景,X射線注量隨著爆心距的增加快速下降,其本身的作用效果(這里不包含由X射線引發(fā)的電磁脈沖等其他間接效應(yīng))范圍較為有限,一般在千米量級;對于150 km及以上高度的核爆炸,往往可以使用簡單的平方衰減率近似估算其能注量。

2.2 附加電離源建模及其在非均勻大氣中的輸運

高空核爆炸產(chǎn)生的瞬發(fā)輻射(中子、γ和X射線)和緩發(fā)輻射(γ射線和β粒子)會電離高空大氣,造成電離層電子數(shù)密度劇增,形成所謂的“附加電離層”,從而破壞無線電通信、雷達等系統(tǒng)的通信鏈路。附加電離層模擬的難點技術(shù)問題之一就是電離源建模及其在非均勻大氣中輸運的蒙特卡羅模擬。

瞬發(fā)中子、γ和X射線在爆炸瞬間釋放出來,可視為各項同性的點源。緩發(fā)γ射線和β粒子則由裂變產(chǎn)物形成的碎片云持續(xù)釋放,輻射源是一個時間空間持續(xù)變化的體源。當(dāng)爆炸高度和當(dāng)量不同時,碎片云形狀可由球形演化為扁橢球形或倒梨形等,可采用流體力學(xué)模型描述碎片云的時空演化過程[31-32],計算得到碎片云運動參數(shù)。針對碎片云形狀隨時間不斷演化特性,采用輻射源分層抽樣的方法得到緩發(fā)γ射線的初始位置。根據(jù)不同時刻的碎片云高度和幾何形狀,將碎片云等份額劃分為多層,每層等效為一定高度上均勻分布的薄圓柱體,倒梨形碎片云形成的緩發(fā)γ射線源如圖14所示。通過分層均勻抽樣實現(xiàn)對不規(guī)則體源的抽樣。

圖14 倒梨形碎片云形成的緩發(fā)γ射線源Fig.14 The delayed γ-ray source from the inverted pear-shaped debris

由于高空大氣密度隨高度增加而指數(shù)衰減,傳統(tǒng)蒙特卡羅方法要實現(xiàn)步長抽樣,只能將高空大氣分層處理,劃分為內(nèi)部密度均勻的密集薄層。每層大氣厚度越小,計算結(jié)果越精確,同時計算量也越大。本文建立了基于質(zhì)量厚度抽樣的蒙特卡羅方法,在粒子輸運過程中,對兩次碰撞點之間的空氣質(zhì)量厚度進行抽樣,再將質(zhì)量厚度轉(zhuǎn)換為每一步的步長。質(zhì)量厚度抽樣代替步長抽樣,無須對高空大氣進行分層處理,可簡化大氣幾何模型,提高蒙特卡羅方法計算效率。采用該方法建立了中子、光子在高空大氣中輸運的專用蒙特卡羅程序MCATNP[33]。圖15為MCNP和MCATNP程序計算給出的X、γ和中子能注量隨距離的變化關(guān)系。其中,計算場景爆高取為100 km,當(dāng)量取為1 Mt TNT,X射線能譜采用黑體溫度為1.5 keV的普朗克譜,γ射線和中子能譜分別采用氫彈出殼γ射線和中子的能譜[8]。由圖15可見,X、γ和中子能注量均隨距離的增大成指數(shù)衰減關(guān)系,且MCATNP模擬結(jié)果與MCNP結(jié)果吻合得較好,表明質(zhì)量厚度抽樣方法可用于計算高空核爆炸產(chǎn)生的X、γ和中子能注量。

圖15 MCNP和MCATNP程序計算給出的X、γ和中子能注量隨距離的變化關(guān)系Fig.15 The X, γ and neutron energy fluences vs. distance computed by MCNP and MCATNP

基于碎片云運動參數(shù)和裂變產(chǎn)物設(shè)置電離源,采用基于質(zhì)量厚度抽樣的蒙特卡羅方法跟蹤模擬瞬發(fā)輻射和緩發(fā)γ輻射電離大氣分布,結(jié)合緩發(fā)β粒子電離大氣分布數(shù)值求解電離連續(xù)性方程組,進而得到高空核爆炸電離大氣產(chǎn)生的附加電離層環(huán)境。爆高為50 km,當(dāng)量為50 kt TNT,爆后5 min電離層電子數(shù)密度的分布如圖16所示。由圖16可見,最高處比背景電離層高出3～4個量級。

圖16 爆高為50 km,當(dāng)量為50 kt TNT,爆后5 min電離層電子數(shù)密度的分布Fig.16 The electron number density distribution at 5 min with the detonation height of 50 km and the yield of 50 kt TNT

2.3 基于蒙特卡羅與帶電粒子推進相結(jié)合方法的人工輻射帶數(shù)值模擬

高空核爆炸裂變碎片在運動過程中衰變產(chǎn)生的β電子被地磁場捕獲注入到自然輻射帶中,形成人工輻射帶。人工輻射帶能覆蓋全球、電子注量率高、持續(xù)時間長,會對衛(wèi)星等空間目標產(chǎn)生總劑量和充放電等多種破壞效應(yīng)。核爆炸釋放帶電粒子注入輻射帶的過程可采用基于蒙特卡羅與帶電粒子推進相結(jié)合[16]的方法,通過跟蹤核爆炸主要裂變衰變鏈碎片離子在地磁場和大氣作用下的運動,統(tǒng)計獲取碎片離子和衰變β電子的捕獲效率及爆后初始電子注量率隨經(jīng)度、磁殼參數(shù)、能量及赤道傾角余弦的分布。

核爆炸裂變產(chǎn)物有近千種放射性同位素,不同中子能量誘發(fā)裂變材料裂變產(chǎn)額不同。為簡化計算,選取總產(chǎn)額較高的質(zhì)量鏈并整理構(gòu)建裂變產(chǎn)物數(shù)據(jù)庫,包括所選取質(zhì)量鏈的總產(chǎn)額及各質(zhì)量鏈中考慮核素的質(zhì)量數(shù)、質(zhì)子數(shù)、獨立產(chǎn)額及半衰期等參數(shù)。

在構(gòu)建高空核爆炸后初始碎片離子的空間、速度及方向分布模型的基礎(chǔ)上,采用隨機蒙特卡羅抽樣方法,確定碎片離子的初始位置;采用Boris推進方法跟蹤碎片離子在地磁場和大氣作用下的運動[34];采用壽命概率模型結(jié)合裂變產(chǎn)物數(shù)據(jù)庫進行隨機抽樣,確定放射性碎片發(fā)生β衰變的時刻和位置,并根據(jù)裂變衰變β電子能譜抽取電子能量。模擬過程中如裂變產(chǎn)物發(fā)生β衰變,根據(jù)裂變衰變質(zhì)量鏈,裂變碎片采用新的子體,其速度和位置仍用母體的速度和位置,繼續(xù)跟蹤裂變碎片直到其子體為穩(wěn)態(tài)核或能量小于截斷能量。

根據(jù)上述方法通過跟蹤大量碎片離子的運動,得到裂變衰變β電子的位置、能量和方向,進而計算爆后初期赤道面電子注量率分布。在此基礎(chǔ)上,結(jié)合輻射帶電子長期擴散損失模擬方法,實現(xiàn)了高空核爆炸人工輻射帶環(huán)境參數(shù)的數(shù)值模擬。爆高為400 km,當(dāng)量1.4 Mt核爆炸的初始捕獲電子數(shù)Ne隨磁殼參數(shù)L的分布和電子全向注量率φe隨時間t的變化關(guān)系如圖17所示。由圖17(a)可見,初始捕獲電子隨磁殼參數(shù)的分布并不是以爆心為中心,而是在高L值上的分布占多,這是由于碎片離子的徑向擴散使其向更高方向運動;由圖17(b)可見,電子注量率隨時間近似成指數(shù)形式衰減關(guān)系,衰減到本底需2～3 a。

(a) Ne vs. L

(b) φe vs. t

2.4 基于蒙特卡羅模擬的高空核電磁脈沖源項計算方法

高空核電磁脈沖按照產(chǎn)生機理分為早期(E1)、中期(E2)和晚期(E3)電磁脈沖。直達γ、散射γ和中子次級γ射線是產(chǎn)生E1和E2的主要因素。在以往的E1模擬中,往往直接采用解析公式計算直達γ注量作為輻射源項,忽略了散射γ射線的影響,也忽略了γ射線時間譜、能譜和角度譜特性。E2的主要輻射源項為中子次級γ射線,目前尚無解析公式可計算其注量。因此,本文提出了基于蒙特卡羅模擬的高空核電磁脈沖源項計算方法,通過跟蹤模擬核爆炸γ射線和中子在高空非均勻大氣中的輸運過程,獲得E1和E2的源項參數(shù)。

向下傳播的γ射線在距地面20～40 km高度形成沉積區(qū),激勵產(chǎn)生電磁脈沖。核爆炸γ射線和中子由爆點輸運至沉積區(qū)的蒙特卡羅模擬面臨深穿透問題。圖18為γ射線和中子在高空大氣中輸運的幾何模型示意圖。由于計算空間尺度長達數(shù)百千米,且高空大氣不斷與γ射線和中子發(fā)生吸收和散射作用,直接蒙特卡羅模擬時,到達記錄點(尤其是高度為20 km的記錄點)附近的模擬粒子極少,使蒙特卡羅模擬的統(tǒng)計漲落較大,模擬結(jié)果不可信。必須采用源偏倚、幾何分裂和時間分裂等綜合減方差技術(shù)[35],提高計算效率與計算精度。爆高為80 km、當(dāng)量為1 Mt TNT核爆炸情況下,高度為20 km處的不同組成和方向的γ注量率φγ隨時間的變化關(guān)系如圖19所示。由圖19(a)可見,1 μs之前,直達γ射線占主導(dǎo);1～100 μs,散射γ射線占主導(dǎo);100 μs后,中子次級γ射線占主導(dǎo)。

圖18 γ射線和中子在高空大氣中輸運的幾何模型示意圖Fig.18 The sketch of the geometry model of γ-rays and neutron transport in the upper atmosphere

粒子的角分布是高空核電磁脈沖源項的關(guān)鍵參數(shù)之一。蒙特卡羅模擬中常見的基于面計數(shù)的角分布統(tǒng)計方法僅適用于面流量記錄方式,難以給出滿足電磁脈沖計算精度的粒子角分布。本文提出了基于點通量記錄的粒子角分布計算方法。在點通量記錄中,設(shè)置參考方向,計算虛粒子的方向與參考方向的夾角,并根據(jù)夾角統(tǒng)計計算不同角度分區(qū)內(nèi)的粒子注量。由圖19(b)可見,在1 μs以前,直達γ射線占主導(dǎo)地位,γ射線方向幾乎均沿r方向;1 μs之后,由于散射γ和中子次級γ方向分布的復(fù)雜性,γ射線沿各方向分量均有分布。

(a) Composition

(b) Direction

HEMP模擬計算時所需求的最直接的參數(shù)就是康普頓電子信息。前人的研究中往往根據(jù)γ射線注量率、能量來計算康普頓電子參數(shù),忽略了隨時間變化的γ射線能量分布和角分布。本文通過在點通量記錄中構(gòu)造基于康普頓散射理論的探測器響應(yīng)函數(shù),直接統(tǒng)計計算康普頓電子的產(chǎn)生率,同時利用康普頓散射理論計算康普頓電子的數(shù)量、能量和方向等信息,最后根據(jù)康普頓電子產(chǎn)生率和能量方向信息計算得到初始康普頓電流密度。該方法直接由蒙特卡羅模擬給出初始康普頓電子信息,可考慮γ射線和康普頓電子的時間、能量和方向等影響因素。

根據(jù)蒙特卡羅模擬計算的康普頓電子信息,計算得到初始康普頓電流和傳導(dǎo)電流,代入麥克斯韋方程組,求解得到電磁脈沖波形與電場強度[36]。圖20為不同輻射源項產(chǎn)生的電場波形。由圖20可見,電場波形明顯具有直達γ、散射γ和中子次級γ射線的特征;直達γ射線引起早期E1快速上升的電場強度峰值,最大峰值約為50 kV·m-1;散射γ射線在1 ～100 μs占主導(dǎo)地位,對應(yīng)峰值電場強度可達100 V·m-1;中子次級γ射線在10-3～10-2s引起第二個峰,峰值電場強度約為10 V·m-1。

3 結(jié)論與展望

本文作者團隊經(jīng)長期持續(xù)研究,采用蒙特卡羅模擬方法,建立了大氣層核爆炸與高空核爆炸諸多核輻射環(huán)境參數(shù)的數(shù)值模擬計算模型與方法。在大氣層核爆炸方面,構(gòu)建了沖擊波影響下緩發(fā)γ輻射輸運模型和大規(guī)模復(fù)雜城市條件下粒子輸運模型,建立了權(quán)窗結(jié)合密度逼近迭代的減方差方法和GVR方法等一系列減方差方法,解決了核輻射在近地面大氣遠距離輸運模擬的深穿透問題。在高空核爆炸方面,建立了考慮裸核區(qū)的X射線輸運模型、附加電離源建模及其在非均勻大氣中的輸運模型、基于蒙特卡羅與帶電粒子推進相結(jié)合的人工輻射帶注入模型和基于蒙特卡羅模擬的高空核電磁脈沖源項計算方法,實現(xiàn)了高空核爆炸X射線、附加電離、人工輻射帶和核電磁脈沖等輻射環(huán)境數(shù)值模擬。

核爆炸輻射環(huán)境涉及因素多,機制機理復(fù)雜,空間尺度大,時間跨度大,在數(shù)值模擬方面存在較多需要進一步研究解決的問題。在核爆炸輻射環(huán)境計算中,某些計算模型過于簡單,氣象條件、地形地貌和建筑物類型等實際復(fù)雜因素考慮不充分不全面。在蒙特卡羅模擬中,為解決遠距離輸運帶來的深穿透問題,建立了相應(yīng)的減方差方法提高計算效率與精度,但一旦變換場景條件相關(guān)的參數(shù),需重新迭代計算權(quán)窗等參數(shù),耗費大量時間。

下一步將梳理總結(jié)大氣層和高空核爆炸核輻射環(huán)境模擬中的建模方法和蒙特卡羅減方差方法,建立一系列核輻射環(huán)境參數(shù)計算的標準規(guī)范,開發(fā)自主專用的核輻射環(huán)境數(shù)值模擬程序,切實提高核輻射環(huán)境數(shù)值模擬能力。