核爆炸地面感生放射性輻射環(huán)境的蒙特卡羅模擬

商 鵬，黃流興，牛勝利，左應(yīng)紅，朱金輝

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

核爆炸釋放的中子經(jīng)大氣輸運(yùn)后，一部分入射到地面的中子被地面介質(zhì)中的某些元素俘獲，發(fā)生活化反應(yīng)產(chǎn)生放射性核素。這些放射性核素衰變放出β射線(xiàn)和γ射線(xiàn)，稱(chēng)為感生放射性。土壤的感生放射性是地面放射性沾染的來(lái)源之一，在爆高較高的條件下，土壤感生放射性是地面放射性沾染的主要來(lái)源[1]，對(duì)周?chē)藛T和環(huán)境的影響不容忽視。

土壤的感生放射性主要由中子與土壤中的多種元素與中子發(fā)生俘獲反應(yīng)引起的。原子核俘獲中子后被活化，形成激發(fā)態(tài)的復(fù)合核，該復(fù)合核會(huì)在很短的時(shí)間內(nèi)(約10-14s)退激發(fā)向外輻射出一條或多條γ射線(xiàn)，稱(chēng)為中子俘獲反應(yīng)，即 (n,γ) 反應(yīng)[2]。土壤中能被活化的元素很多，一般而言，活化后的感生放射性主要來(lái)自28Al，56Mn，24Na，59Fe，對(duì)應(yīng)的靶核分別是27Al，55Mn，23Na，58Fe。

一般用距地面一定高度處的γ劑量率來(lái)表征感生放射性形成的輻射環(huán)境。文獻(xiàn)[1, 3-4]中采用點(diǎn)源近似法[5]，結(jié)合經(jīng)驗(yàn)累積系數(shù)修正推導(dǎo)了探測(cè)點(diǎn)處的γ放射性活度公式，并通過(guò)半無(wú)限空間積分，進(jìn)一步給出了土壤感生放射性在探測(cè)點(diǎn)處產(chǎn)生的γ劑量率表達(dá)式。上述方法在計(jì)算過(guò)程中進(jìn)行了較多的簡(jiǎn)化和近似，如中子注量采用能量小于0.5 eV的等效熱中子[6]近似及中子活化截面采用等效熱中子平均活化截面近似等，只能給出粗略的結(jié)果。

為合理而相對(duì)準(zhǔn)確地計(jì)算感生放射性，武禎等[7]針對(duì)電子直線(xiàn)加速器感生放射性的計(jì)算，對(duì)比分析了點(diǎn)源近似法、徑跡長(zhǎng)度法、功率吸收法和蒙特卡羅模擬法，研究表明，蒙特卡羅模擬法考慮反應(yīng)全面，可同時(shí)計(jì)算直接和間接感生放射性，計(jì)算步驟簡(jiǎn)單，結(jié)果準(zhǔn)確，是計(jì)算感生放射性的“最好”方法。本文細(xì)致分析了感生放射性的物理過(guò)程，基于蒙特卡羅方法建立了核爆炸地面感生放射性?xún)刹接?jì)算方法。

1 核爆炸地面感生放射性?xún)刹接?jì)算方法

蒙特卡羅模擬方法能較逼真地描述大量粒子在介質(zhì)中輸運(yùn)的物理過(guò)程，且收斂速度與問(wèn)題維度無(wú)關(guān)，程序結(jié)構(gòu)簡(jiǎn)單，易于實(shí)現(xiàn)，因此，蒙特卡羅方法也被喻為“最后的方法”。模塊化粒子輸運(yùn)程序包PHEN[8]由西北核技術(shù)研究所開(kāi)發(fā)，能模擬中子、光子、電子及質(zhì)子與物質(zhì)的相互作用過(guò)程，如中子的彈性散射、非彈性散射和輻射俘獲等，光子的光電吸收、康普頓散射和電子對(duì)產(chǎn)生等。它還可針對(duì)特定物理問(wèn)題進(jìn)行進(jìn)一步開(kāi)發(fā)，統(tǒng)計(jì)多種物理量。目前，PHEN程序并不具備直接模擬感生放射性的功能，需通過(guò)構(gòu)建中子輸運(yùn)分層統(tǒng)計(jì)模型和γ輻射體源模型，并采用模型參數(shù)轉(zhuǎn)換技術(shù)使二者結(jié)合，進(jìn)而開(kāi)展模擬計(jì)算，即地面感生放射性?xún)刹接?jì)算方法，計(jì)算流程如圖1所示。

1.1 中子輸運(yùn)分層統(tǒng)計(jì)模型

核爆炸后輸運(yùn)到地面的有直射和散射中子，本文建立了中子輸運(yùn)分層統(tǒng)計(jì)模型，如圖2所示。

中子源[6]各向同性發(fā)射中子，在均勻大氣中輸運(yùn)，最終入射到地面土壤中，與土壤中的靶核發(fā)生相互作用，產(chǎn)生放射性核素28Al，56Mn，24Na，59Fe，記錄不同深度處放射性核素的數(shù)密度。考慮到不同深度的土壤活化程度并不相同，傳統(tǒng)方法中采用的均勻化假設(shè)與實(shí)際情況相差較大，所以采用中子輸運(yùn)分層統(tǒng)計(jì)方法。中子與土壤介質(zhì)中的原子核發(fā)生相互作用，可能發(fā)生(n,γ)，(n,α)，(n,p)，(n,2n)，(n,3n)等多種反應(yīng)。其中，感生放射性γ劑量率計(jì)算主要關(guān)注中子在土壤中輸運(yùn)并發(fā)射γ射線(xiàn)的過(guò)程，即(n,γ)反應(yīng)過(guò)程。因此，需記錄中子輸運(yùn)過(guò)程中發(fā)生的(n,γ)反應(yīng)的次數(shù)。在PHEN程序包中，將中子在輸運(yùn)過(guò)程中與土壤某核素發(fā)生的所有(n,γ) 反應(yīng)進(jìn)行統(tǒng)計(jì)，得到土壤中該核素被活化后生成的放射性核素個(gè)數(shù)；對(duì)不同的核素重復(fù)上述過(guò)程，得到所有放射性核素的數(shù)密度。

上述計(jì)算過(guò)程可表示為

(1)

其中：Φ(E)為分能群中子注量，cm-2；σc(E)為(n,γ)反應(yīng)的截面，b；常數(shù)η為任意可用于歸一化處理的標(biāo)量，也稱(chēng)為歸一化因子。當(dāng)η的單位為b-1·cm-1時(shí)，式(1)的計(jì)算結(jié)果即為單位體積內(nèi)發(fā)生(n,γ)反應(yīng)的總數(shù)nc，它在數(shù)值上與放射性核素的數(shù)密度nrad相等。

設(shè)距爆心為R，距地面深度為x，能量為E的中子注量為Φ(R,x,E)，則距爆心R處，不同深度x，第i種的放射性核素的數(shù)密度nrad(R,x，i)可表示為

(2)

其中：ni為第i種靶核的原子數(shù)密度，cm-3，與土壤中該元素的含量ci及天然豐度αi有關(guān)，可表示為

(3)

其中：NA為阿伏加德羅常數(shù)，mol-1；ρg為模擬區(qū)域土壤介質(zhì)的密度，取1.4 g·cm-3；Mi為第i種靶核的摩爾質(zhì)量，g·mol-1，數(shù)值與相對(duì)原子質(zhì)量Ai相等。對(duì)比式(1)和式(2)可知，歸一化因子ηi=ni×10-24b-1·cm-1。

計(jì)算中使用的幾種元素的計(jì)算參數(shù)及由式(3)計(jì)算得到的歸一化因子如表1所列。

表1 27Al，55Mn，23Na，58Fe的相關(guān)計(jì)算參數(shù)Tab.1 Computational parameters of 27Al，55Mn，23Na and 58Fe

1.2 γ輻射分層體源模型

核爆炸發(fā)生后，土壤活化產(chǎn)生的放射性核素發(fā)生衰變，向外各向同性地放出γ射線(xiàn)，這個(gè)過(guò)程相對(duì)較慢，導(dǎo)致γ劑量率隨時(shí)間變化。而放射性核素單次衰變放出γ射線(xiàn)并被探測(cè)器探測(cè)的過(guò)程可視為瞬態(tài)過(guò)程。

28Al，56Mn，24Na，59Fe衰變放出的γ射線(xiàn)參數(shù)如表2所列。由表2可知，衰變?chǔ)镁哂须x散的能量值，且分支比rB也有所差別，這是由放射性核素自身性質(zhì)決定的。

表2 28Al，56Mn，24Na，59Fe的衰變?chǔ)蒙渚€(xiàn)參數(shù)[9]Tab.2 Parameters of decay γ rays for 28Al，56Mn，24Na and 59Fe

為模擬放射性核素出射γ射線(xiàn)對(duì)探測(cè)點(diǎn)探測(cè)劑量的貢獻(xiàn)，需將中子輸運(yùn)分層統(tǒng)計(jì)模型計(jì)算得到的放射性核素?cái)?shù)密度分布轉(zhuǎn)換為γ體源模型中不同深度處的γ射線(xiàn)的發(fā)射強(qiáng)度。同時(shí)，γ輻射源還應(yīng)當(dāng)體現(xiàn)表2列出的能量、分支比及產(chǎn)額等特性。放射性核素發(fā)射γ射線(xiàn)數(shù)與放射性核素?cái)?shù)密度及分層?xùn)旁w積相關(guān)，因此，構(gòu)建γ輻射分層體源模型，模型特征參數(shù)由中子輸運(yùn)分層統(tǒng)計(jì)模型進(jìn)行參數(shù)轉(zhuǎn)換得到，如圖3所示。

如1次衰變的γ射線(xiàn)產(chǎn)額為Y，能量離散分布，分布概率即為對(duì)應(yīng)能量的γ射線(xiàn)分支比。將源權(quán)重賦值為Y，則模擬結(jié)果的物理含義為體源中發(fā)生1次衰變放出的Y個(gè)具有不同能量的γ射線(xiàn)對(duì)探測(cè)點(diǎn)探測(cè)劑量的貢獻(xiàn)。

由上述分析可知，γ輻射分層體源模型中的輻射源為一個(gè)能量離散分布、位置具有特定分布及各向同性的復(fù)雜體源。離散能量抽樣和方向各向同性抽樣算法可參考文獻(xiàn)[13]。γ射線(xiàn)在分層體源模型中的位置抽樣采取封閉圓柱體包含多個(gè)體源柵元的舍選方法(rejection method,RM)，如圖3所示。由于地面不同深度處的放射性核素?cái)?shù)密度不同，衰變發(fā)射γ射線(xiàn)的強(qiáng)度也不同，即每個(gè)柵元發(fā)射γ射線(xiàn)的概率與柵元內(nèi)的放射性核素個(gè)數(shù)成正比。設(shè)柵元j經(jīng)歸一化后的γ射線(xiàn)發(fā)射概率為Pj，在模擬程序中對(duì)γ射線(xiàn)發(fā)射位置進(jìn)行抽樣時(shí)，設(shè)置一個(gè)虛擬圓柱體將柵元體源完全包圍，并在虛擬圓柱體內(nèi)進(jìn)行均勻抽樣；若抽樣得到的粒子位置位于柵元j中，則以Pj的概率接受該位置抽樣，否則舍去本次抽樣，重新抽取位置信息；重復(fù)上述舍選過(guò)程，直到一次抽樣的位置信息被接受。虛擬圓柱體需完全包圍發(fā)射γ射線(xiàn)的多個(gè)柵元，但不應(yīng)設(shè)置得過(guò)大，否則將降低位置抽樣效率。

利用PHEN程序包分別計(jì)算28Al，56Mn，24Na，59Fe單次放射性核素衰變放射的γ射線(xiàn)對(duì)探測(cè)點(diǎn)的平均劑量貢獻(xiàn)Di。對(duì)于整個(gè)多層γ輻射體源而言，柵元j中第i種放射性核素的數(shù)密度ni,j與體積Vj的乘積之和即為第i種放射性核素總數(shù)Ni，表示為

(4)

如果放射性核素全部衰變，則衰變產(chǎn)生的總劑量D(t)為所有放射性核素的劑量貢獻(xiàn)之和，表示為

(5)

(6)

其中：λi為第i種放射性核素的衰變常量；t為衰變時(shí)間。由式(6)可知，γ衰變劑量率受到放射性活度的影響，同樣具備指數(shù)衰減規(guī)律。

2 兩步計(jì)算方法應(yīng)用實(shí)例

以爆高為200 m，典型氫彈(歸一化到1 kt當(dāng)量)的核爆炸場(chǎng)景為例，利用第1節(jié)所述地面感生放射性?xún)刹接?jì)算方法，研究了γ輻射分層體源的半徑和厚度對(duì)感生放射性的影響，選取了一組較為合理的模型參數(shù)，并采用優(yōu)化參數(shù)計(jì)算了地面上方1 m處的γ總劑量率及放射性核數(shù)劑量貢獻(xiàn)比。

2.1 土壤層厚度及半徑參數(shù)

文獻(xiàn)[1]指出，地面上方1 m處的γ劑量率主要來(lái)自半徑為10～20 m和深度為10～20 cm土壤薄層。本文通過(guò)建立中子輸運(yùn)分層統(tǒng)計(jì)模型，模擬計(jì)算了半徑為20 m，土壤層厚度為2 cm和4 cm時(shí)，放射性核素28Al的數(shù)密度隨深度的變化關(guān)系，如圖4所示。

由圖4可見(jiàn)，放射性核素?cái)?shù)密度最大值均出現(xiàn)在10～20 cm深度處，與文獻(xiàn)[1]所述一致?？捎脠D5所示27Al的中子輻射俘獲截面隨能量的變化關(guān)系來(lái)解釋?zhuān)喝肷渲凶又心芰枯^高的快中子成分在土壤中輸運(yùn)并在被慢化過(guò)程中產(chǎn)生更多的熱中子，使土壤一定深度處的熱中子水平高于地面，而熱中子的輻射俘獲截面比快中子高的多，導(dǎo)致由輻射俘獲反應(yīng)產(chǎn)生的放射性核素的數(shù)密度隨深度呈先上升后下降的趨勢(shì)。放射性核素的數(shù)密度隨深度的變化關(guān)系如圖6所示。

改變?chǔ)幂椛潴w源模型土壤層的厚度，計(jì)算得到不同厚度γ輻射分層體源模型中28Al 的總數(shù)及衰變?chǔ)每倓┝?，如?所列。由表3可知，采用蒙特卡羅方法模擬計(jì)算地面γ劑量，模型厚度小于60 cm時(shí)，有一部分衰變?chǔ)么┩竿寥缹?，被終止輸運(yùn)，導(dǎo)致測(cè)點(diǎn)有效記數(shù)偏低；模型厚度大于60 cm后，記數(shù)基本穩(wěn)定。因此，后續(xù)計(jì)算中，土壤層厚度取60 cm。

確定土壤層厚度后，改變?chǔ)幂椛潴w源模型的半徑，得到衰變?chǔ)每倓┝侩S半徑的變化關(guān)系，如圖7所示。由圖7可見(jiàn)，模型半徑為40 m與模型半徑為100 m時(shí)，衰變?chǔ)每倓┝肯嗖畈怀^(guò)5%。為節(jié)省粒子輸運(yùn)時(shí)間提高計(jì)算效率，取模型半徑為40 m。

表3 不同厚度γ輻射分層體源模型中28Al 的總數(shù)及衰變?chǔ)每倓┝縏ab.3 N of 28Al and D for γ layered-radiator source model with different thickness

2.2 地面衰變?chǔ)脛┝柯视?jì)算

2.1節(jié)研究表明，采用蒙特卡羅方法計(jì)算土壤中放射性核素出射γ對(duì)地面上方探測(cè)點(diǎn)探測(cè)劑量的貢獻(xiàn)時(shí)，取模型半徑為40 m，厚度為60 cm可滿(mǎn)足計(jì)算要求。

按照兩步計(jì)算法計(jì)算流程，利用中子輸運(yùn)分層統(tǒng)計(jì)模型計(jì)算給出感生放射性核素的數(shù)密度分布，利用放射性核素?cái)?shù)密度分布及放射性核素衰變特性構(gòu)建γ輻射分層體源模型，模擬樣本數(shù)為1×107個(gè)光子，分別計(jì)算28Al，56Mn，24Na，59Fe單次衰變的平均劑量貢獻(xiàn)與放射性核素總數(shù)，如表4所列。

表4 28Al， 56Mn，24Na，59Fe的單次衰變平均劑量貢獻(xiàn)與放射性核素總數(shù)Tab.4 Average dose contributions of single decay and total radionuclide amounts for 28Al， 56Mn，24Na and 59Fe

將放射性核素單次衰變平均劑量貢獻(xiàn)與總數(shù)代入式(6)，可得土壤感生放射性總劑量率隨時(shí)間的變化關(guān)系，如圖8所示。由圖8可見(jiàn)，總劑量率整體呈下降趨勢(shì)，但下降速率發(fā)生了幾次較大的變化，使總劑量率變化曲線(xiàn)表現(xiàn)出“臺(tái)階”狀變化的形式。

這一時(shí)間特性同樣可從圖9所示放射性核素的劑量貢獻(xiàn)比隨時(shí)間的變化關(guān)系分析得到。由圖9還可見(jiàn)，4種放射性核素分別在不同的時(shí)間區(qū)間上發(fā)揮主要作用，具體而言，在爆炸后的數(shù)分鐘內(nèi)，地面感生放射性主要來(lái)自28Al，強(qiáng)度比其他放射性核素高2個(gè)量級(jí)左右；在數(shù)小時(shí)范圍區(qū)間，主要來(lái)自56Mn和24Na；而在200 h之后，基本只來(lái)自59Fe。

3 結(jié)論及展望

本文基于蒙特卡羅方法建立了核爆炸地面感生放射性?xún)刹接?jì)算方法。通過(guò)敏感性分析給出了構(gòu)建中子輸運(yùn)分層統(tǒng)計(jì)模型和γ輻射分層體源模型的典型參數(shù)，修正了傳統(tǒng)半經(jīng)驗(yàn)公式方法中的等效熱中子近似和等效熱中子平均活化截面近似，考慮了土壤活化后放射性核素?cái)?shù)密度隨深度的變化及γ射線(xiàn)在空氣中輸運(yùn)時(shí)的散射和衰減，并更細(xì)致地描述了感生放射性輻射環(huán)境形成的物理過(guò)程。本文方法還能應(yīng)用于人工放射性輻射源活化其他材料后的γ放射性輻射環(huán)境計(jì)算，為高精度輻射環(huán)境計(jì)算提供了有效技術(shù)手段。