錳浴法測量中子源發射率中各修正因子的蒙特卡羅計算

李春娟，劉毅娜，張偉華，王志強

(中國原子能科學研究院 計量測試部，北京 102413)

錳浴法測量中子源強度用于標定各種便攜式中子源的發射率，是國際上普遍采用的絕對測量方法，有較高的準確度。它是將待測中子源置于體積很大的裝有MnSO4溶液的容器中心，中子在水中充分慢化后被俘獲，溶液中的55Mn吸收熱中子后，生成放射性核素56Mn，通過56Mn→56Fe+β—+γ發生衰變，半衰期為(2.579±0.003) h。當中子源在錳浴中放置足夠長時間后，錳池中的中子數將達平衡狀態，通過測量此時錳液中56Mn發射出的γ射線，即可得出中子源發射率。

根據錳浴法測量中子源發射率的原理，需對由中子源發出的未被MnSO4溶液中55Mn俘獲的中子組分進行修正，需考慮的修正項有：1) 錳池邊界泄漏中子修正L；2) 溶液中氧和硫通過核反應引起的快中子損失修正O，簡稱氧硫吸收；3) 中子源及源容器對從溶液中反射來的熱中子吸收修正S，簡稱源的自吸收；4) 被Mn俘獲的熱中子占被整個溶液俘獲熱中子的比例，即f。通常，只有不足50%的源中子會被溶液中的55Mn俘獲，因此，修正因子準確與否直接影響測量結果的可靠性，而由于傳統公式的局限性，不能詳細考慮輸運過程及各修正因子之間的相互影響，使得總修正因子的計算較困難，蒙特卡羅(MC)模擬方法在這些計算中發揮了至關重要的作用[1-4]，目前國際普遍采用MC模擬的方法[1]。

中國原子能科學研究院計量測試部(國防科技工業電離輻射一級計量站)建立了錳浴法測量中子源強度標準裝置，并參與了由國際計量局(BIPM)組織的中子源強度國際比對[5]，比對中子源為241Am-Be(α,n)源，其為ISO標準源。本工作利用大型通用MC粒子輸運程序MCNP[6]模擬計算測量過程中的各修正因子，對其不確定度進行詳細評定，以提供一種用于MC模擬計算結果不確定度評定的方法。

1 計算模型的建立

1.1 幾何模型的建立

錳浴法測量中子源強度標準裝置如圖1所示，系統主要由錳池、源容器、循環泵、流量計、測量室、兩個NaI(Tl)探測器及相應的電子學系統組成。測量時將中子源置于錳浴中心36 h，待錳池中的錳液達飽和放射性后，保持硫酸錳溶液在錳池和屏蔽樣品容器內連續循環，對56Mn溶液的飽和活化計數A進行測量。將已知活度(絕對測量得到)的56Mn溶液倒入錳池中，通過兩路NaI(Tl)探測器測量錳液中56Mn的γ計數，確定錳浴探測效率ε[7]。

圖1 中子源強度標準裝置示意圖

圖2 模擬計算的幾何模型

通過MC粒子輸運程序MCNP對測量過程中的各修正因子進行模擬計算。對修正因子的計算有直接影響的結構主要包括不銹鋼錳浴外殼、硫酸錳溶液、241Am-Be(α,n) ISO標準源(X3封裝外殼)、Al桶及不銹鋼容器，據此建立MC模擬計算幾何模型(圖2)。不銹鋼錳浴外殼直徑為100 cm、厚度為4 mm，通過烘干法測量硫酸錳溶液的濃度，計算得到溶液的密度為1.243 9 g/cm3，氫錳原子數比為58.08。241Am-Be(α,n)源為國際電離輻射咨詢委員會組織的國際比對CCRI(Ⅲ)-K9.AmBe的傳遞標準源，源活性區為圓柱形，外徑為17.6 mm，高度為16.3 mm，活性區內AmO2的質量分數為6.29%、Be的質量分數為93.71%，活性區密度為1.24 g/cm3。源采用雙層不銹鋼封裝，封裝外殼型號為X3，外觀直徑為22.4 mm，高度為31 mm。源的承托物包括鋁桶及不銹鋼容器，鋁桶外徑為33 mm，厚度為1.25 mm，不銹鋼容器外徑為38 mm，厚度為1 mm。

1.2 源的描述及中子反應截面數據庫的選取

實驗用241Am-Be(α,n)源的活性區內密度、發射α的氧化镅顆粒大小和封裝外殼均滿足ISO標準要求，故其封裝外殼外的中子能譜與ISO 8529[8]中給出241Am-Be(α,n)源中子能譜一致。國際推薦將比對源置于球形容器中，并對容器抽真空，然后將其置于錳浴中心進行測量。此情況下計算時可將源近似描述為從球形容器內表面發射的各向同性球面源，源中子能譜可描述為ISO 8529中給出的241Am-Be(α,n)源中子能譜。

但本錳浴系統源的承托物為鋁桶及不銹鋼管容器，而源則置于承托物的底部，此結構會改變源中子入射的邊界條件。為使計算更接近真實情況，確保計算結果的可靠性，計算中將源項描述為從源活性區發射的均勻分布的圓柱形體源。封裝外殼內從源活性區發射的中子能譜是未知的，利用MCNP通過多次迭代的計算方法得到封裝外殼內的裸體芯區中子能譜。具體計算過程如下。

首先利用ISO 8529給出的能譜作為裸體芯區體源的初始能譜，并采用相同的能群結構記錄封裝外殼外的能譜，計算出其與ISO 8529能譜的差別，將其疊加在初始能譜上作為裸體芯區體源源中子能譜進行重新計算，再次記錄封裝外殼外的能譜并與ISO 8529能譜進行比較，如此迭代計算，直到包殼外的能譜與ISO 8529能譜非常接近，此時的源能譜即可作為最終的源活性區體源中子能譜。修正過程中各次計算輸入能譜的比較如圖3所示，各次輸出能譜與ISO 8529能譜的比較如圖4所示。

圖3 各次計算輸入能譜的比較

圖4 各次輸出能譜與ISO 8529能譜的比較

另外，由于中子與源活性區材料及其不銹鋼封裝外殼的(n,2n)反應，計算得到源材料發射1個中子，封裝外殼外的中子數約為1.016 6個，故將源中子的權重設為0.983 64，以使計算結果對源封裝外殼外的1個源中子進行歸一。

MCNP對中子輸運過程的MC模擬基于中子與物質相互作用的反應截面，本計算中使用的中子反應截面數據均選自于ENDF/B-7[9]評價庫，使用程序NJOY99[9]將其作為MCNP所需數據庫的格式。

2 修正因子計算方法

為得到各修正因子并進行自洽性驗證，模擬計算主要記錄以下各量：錳池邊界泄漏中子份額L，源活性區俘獲中子份額S1，不銹鋼封裝外殼俘獲中子份額S2，不銹鋼源容器俘獲中子份額S3，Al桶俘獲中子份額S4，硫酸錳溶液俘獲中子份額Nabs，硫酸錳溶液中55Mn通過56Mn(n,γ)反應所俘獲的熱中子份額f1，硫酸錳溶液中S通過S(n, γ)反應所俘獲的熱中子份額f2，硫酸錳溶液中S通過S(n, p)、S(n, α)反應所俘獲的快中子份額O1、O2，硫酸錳溶液中O通過O(n,γ)反應所俘獲的熱中子份額f3，硫酸錳溶液中O通過O(n, α)反應所俘獲的快中子份額O3，硫酸錳溶液中H通過H(n,γ)反應所俘獲的熱中子份額f4。其中，利用MCNP中的F1計算錳池邊界泄漏中子，利用柵元注量計數器F4及乘子卡FM計算得到上述各核素通過核反應俘獲的中子。硫酸錳溶液俘獲中子份額計算方法如下。

設硫酸錳溶液中55Mn、32S、16O、1H的原子核數分別為n1、n2、n3、n4，它們與能量為E的中子的吸收截面分別為σ1abs(E)、σ2abs(E)、σ3abs(E)、σ4abs(E)，中子注量率譜為φ(E)，則硫酸錳溶液俘獲中子總數Nabs即可通過式(1)得到。利用MCNP中的柵元注量計數器F4可記錄硫酸錳溶液中的平均中子注量率譜，乘子卡FM可通過引用各核素的原子核數、中子與各核素相互作用的反應截面及由F4記錄得到的中子注量率譜，得到硫酸錳溶液俘獲中子總數。由于該計算對源封裝外殼外的1個源中子進行歸一，故模擬計算得到的硫酸錳溶液俘獲中子總數等于硫酸錳溶液俘獲中子份額。

n3σ3abs(E)+n4σ4abs(E))φ(E)dE

(1)

其他各量的計算方法與硫酸錳溶液俘獲中子份額計算方法類似。通過以上各量即可得到各修正因子，其中錳池邊界泄漏中子修正因子L可直接計算得到，源的自吸收修正因子S=S1+S2+S3+S4，55Mn俘獲的熱中子占被整個溶液俘獲的熱中子的比例f=f1/(f1+f2+f3+f4)，氧硫吸收修正因子O=O1+O2+O3。得到各修正因子后，通過式(2)，即可得到中子源發射率Q。

(2)

式中：A為測量得到的錳液飽和放射性活化計數；ε為錳浴系統的探測效率。

另外，若不考慮由中子源發出的未被MnSO4溶液中55Mn俘獲的中子組分，而直接考慮由中子源發出的最終被硫酸錳溶液中55Mn通過56Mn(n, γ)反應俘獲的中子份額F，也可直接由式(3)計算得到中子源發射率。由式(3)和式(1)計算得到的結果應該是一致的，F可通過MC模擬計算得到，即上述記錄量中的f1。這種計算方法較簡單，也可簡化不確定度評價過程。但為對測量過程中各影響量有較詳細的了解，并驗證計算模型的可靠性，有必要對由中子源發出的未被MnSO4溶液中55Mn俘獲的各中子組分進行計算。

Q=A/Fε

(3)

3 計算結果及實驗驗證

對10 000 000個粒子進行跟蹤模擬，計算結果列于表1、2。從表1、2可看出：計算結果是自洽的，即直接計算得到的硫酸錳溶液俘獲中子份額與分別計算得到的溶液中Mn、S、O、H所俘獲的中子份額總和在其統計不確定度范圍內一致，各材料的俘獲中子份額及泄漏中子份額總和近似為1，F與f(1-L-O-S)的計算結果在0.6%內一致；錳池邊界泄漏中子份額不足1.5%，可使源中子得到充分慢化并被錳浴吸收。

表1 各材料對中子吸收份額的MC模擬計算結果

表2 各修正因子的MC模擬計算結果

根據式(1)，將修正因子應用于溶液飽和放射性活度實驗測量結果，得到待測241Am-Be(α,n)中子源的發射率，與參加比對的各機構所提供的結果進行比較，結果列于表3。比對結果表明，除法國LEN-LNHB外，本結果與各機構的實驗室的結果在其不確定度范圍內一致，驗證了本計算模型的可靠性。

表3 各機構241Am-Be(α,n)中子源發射率對比[5]

4 模擬計算結果不確定度分析

不確定度評定是中子計量工作的重要組成部分，其表征了測量結果的分散性及可靠程度。根據式(1)及不確定度傳播規律，中子源發射率Q測量結果的不確定度來源主要有56Mn溶液的飽和放射性活化計數A的測量不確定度μ(A)、錳浴系統探測效率的不確定度μ(ε)及各修正量的不確定度。這里主要對各修正量計算結果的不確定度及其對中子源發射率Q測量結果的不確定度的貢獻進行分析。

MCNP程序計算各修正量的不確定度即錳池邊界泄漏中子修正的不確定度μ(L)、溶液中氧和硫通過帶電粒子反應引起的快中子損失修正的不確定度μ(O)、中子源及源容器對從溶液中反射回來的熱中子吸收修正的不確定度μ(S)以及由錳所俘獲的熱中子份額f的不確定度μ(f)。根據不確定度的傳播規律，各修正量的不確定度對中子源發射率不確定度的貢獻可分別用式(4)～(7)計算得到。

(4)

(5)

μf(Q)=μ(f)

(6)

(7)

各修正量的不確定度的來源主要有統計不確定度、反應截面的不確定度、MnSO4溶液密度的不確定度、源和源的承托物材料密度的不確定度及錳浴半徑的不確定度5個方面。首先評估得到各不確定來源引起的各修正量的標準偏差，認為其均滿足矩型分布，得到各不確定度分量，然后通過不確定度的合成法則即可得到各修正量的不確定度，評估結果列于表4。各不確定來源引起的各修正量的標準偏差的評估方法分別如下。

表4 MCNP程序計算各修正量的不確定度

4.1 統計不確定度

統計抽樣的好壞直接影響到計算結果的準確度，MCNP程序對其計算結果進行統計檢查并給出相應統計標準偏差，結果在輸出文件中給出?？赏ㄟ^增加計算粒子數和使用方差減小技術的方法減少統計標準偏差。這里取粒子數為10 000 000進行統計抽樣。

4.2 反應截面的不確定度

MCNP程序計算是基于各反應截面庫的，所用截面數據的不確定度直接影響計算結果。這里所使用的截面數據均選用美國的評價核數據庫ENDF/B-7[9]的評價數據。通過改變不同截面數據庫的辦法評定反應截面的改變對各修正量計算結果的不確定度帶來的貢獻。計算發現，截面數據庫的改變引起的L、O、S、f的最大改變分別為4.68%、21.8%、3.12%、0.04%。由于這里選用的截面數據是較新的數據，最終選定由截面帶來的各修正量的標準偏差分別為4.0%、16%、3.0%、0.04%。

4.3 MnSO4溶液密度的不確定度

MnSO4溶液的濃度的不確定度為0.5%，可計算得到其密度的不確定度，通過MCNP程序的微擾計算功能(PERT卡)得到MnSO4溶液的密度在其不確定度范圍內的改變對各修正量計算結果的影響，從而估計由MnSO4溶液密度的不確定度而引起的各修正量的標準偏差。

4.4 源和源的承托物材料密度的不確定度

源和源的承托物包括源活性區材料、不銹鋼封裝外殼、不銹鋼源容器以及Al桶，估計各種材料的不確定度，通過MCNP程序計算得到源和源的承托物的密度在其不確定度范圍內的改變對各修正量計算結果的影響，從而估計由源和源的承托物材料密度的不確定度引起的各修正量的標準偏差。

由于計算中源中子是從源封裝外殼內發射出的體源，故改變源活性區材料及不銹鋼封裝外殼的密度均會改變封裝外殼外的中子能譜，使其偏離ISO 8529所給出的能譜，故需在源活性區材料及不銹鋼封裝外殼的密度改變的情況下，重新對源中子能譜進行修正，而后得到各修正量，估計由源和源的承托物密度的不確定度而引起的各修正量的標準偏差。計算結果表明，源和源的承托物密度的改變主要對源的自吸收因子影響較大，對其他修正量的影響可忽略。

4.5 錳浴半徑的不確定度

將錳浴的半徑由50 cm改為50.5 cm，計算其對各修正量計算結果的影響，從而估計由錳浴半徑的不確定度而引起的各修正量的標準偏差。

5 結論

利用MCNP程序對錳浴法測量中子源發射率中各修正因子進行計算是可行的，并能彌補傳統計算方法的不足。本工作建立的計算方法和模型是合理的，通過多次迭代計算的方法得到合理的源活性區中子能譜有效地解決了計算中的一個難題，可將此計算方法推廣應用于其他相關問題的計算。另外，本文也為MC模擬計算不確定度的評定提供了一較合理的方法。

參考文獻：

[1] 許淑艷. 蒙特卡羅方法在實驗核物理中的應用[M]. 北京：原子能出版社，1996.

[2] 武禎，李君利. 用于探測器校正因子計算的Monte Carlo方法[J]. 清華大學學報：自然科學版，2006，46(9)：1 585-1 588.

WU Zhen, LI Junli. Monte Carlo method for calculating particle radiation detector correction factors[J]. J Tsinghua Univ: Sci & Tech, 2006, 46(9): 1 585-1 588(in Chinese).

[3] 李春娟，陳軍，吳海成，等.238U裂變電離室修正因子及探測響應的蒙卡模擬計算[J]. 原子能科學技術，2009，43(3)：208-214.

LI Chunjuan, CHEN Jun, WU Haicheng, et al. Monte Carlo simulation for the correction factors and detection response of238U fission chamber[J]. Atomic Energy Science and Technology, 2009, 43(3): 208-214(in Chinese).

[4] KHABAZ R. Appraisement of the correction factors for neutron reaction in the manganese bath using Monte Carlo calculation[J]. Radioanal Nucl Chem, 2012, 293: 455-462.

[5] ROBERTS N J, JONES L N, WANG Z, et al. International key comparison of measurements of neutron source emission rate (1999—2005): CCRI(Ⅲ)-K9.AmBe[J]. Metrologia, 2011, 48(Suppl.): 06018.

[6] BRIESMEISTER E. MCNP: A general Monte CarloN-particle transport code, Version 4C[R]. US: Los Alamos National Laboratory, 2000.

[7] 劉毅娜，王志強，李春娟，等. 錳浴法絕對測量中子源發射率[J]. 原子能科學技術，2013，47(6)：1 044-1 047.

LIU Yina, WANG Zhiqiang, LI Chunjuan, et al. Absolute measurement of neutron source emission rate with manganese bath method[J]. Atomic Energy Science and Technology, 2013, 47(6): 1 044-1 047(in Chinese).

[8] ISO 8529 Neutron reference radiations for calibrating neutron-measuring devices used for radiation protection purposes and for determining their response as a function of neutron energy[S]. Geneva, Switzerland: International Organization for Standardization, 1989.

[9] CHADWICKM B, OBLOZINSKY P, HERMAN M, et al. ENDF/B-Ⅶ.0: Next generation evaluated nuclear data library for nuclear science and technology[J]. Nuclear Data Sheets, 2006, 107(12): 2 931-3 060.

[10] MacFARLANE R E, MUIR D W, BOICOURT R M, et al. The NJOY nuclear data processing system, Volume Ⅰ: User’s manual, LA-9303-M[R]. US: Los Alamos National Laboratory, 1982.