李 磊 李曉燕 黃 瑋 蔣樹斌 伍曉利 楊桂霞

（中國工程物理研究院 核物理與化學(xué)研究所 綿陽 621900）

提高單板鈷源劑量率MCNP計算效率的模擬研究

李 磊 李曉燕 黃 瑋 蔣樹斌 伍曉利 楊桂霞

（中國工程物理研究院 核物理與化學(xué)研究所 綿陽 621900）

蒙特卡羅方法是目前準(zhǔn)確的吸收劑量率計算方法，但其較長的模擬耗時阻礙了它在工業(yè)鈷源輻射加工和輻照實驗中的應(yīng)用。模擬耗時、模擬精度以及模擬值與實測值的相對偏差是表征蒙特卡羅計算效率的重要指標(biāo)。針對8.4PBq的單板鈷源輻照裝置，討論了并行線程數(shù)、記數(shù)方法、記數(shù)柵元尺寸、γ致電子的處理方式和截斷能5種參數(shù)對蒙特卡羅程序MCNP吸收劑量率計算效率的影響。利用實驗測量結(jié)合模擬試算的方法，給出了在保證一定精度和相對偏差前提下，使得模擬耗時最少的參數(shù)組合，提高了MCNP計算效率。結(jié)果如下：超線程模式下的并行計算、*F6記數(shù)方法、柵元邊長為1cm、γ輸運模式、γ截斷能為100keV。

單板鈷源輻照裝置，吸收劑量率，MCNP，計算效率

單板鈷源輻照裝置γ輻射場具有均勻區(qū)域大（幾十至幾百cm）、劑量率高（幾百至上千Gy.min?1）的特點，除輻射加工外，還適合于大尺寸部組件輻照考核以及高劑量率條件下組件/材料輻射效應(yīng)研究。了解輻射場吸收劑量率分布情況是準(zhǔn)確劑量率輻照預(yù)先定位及累積吸收劑量掌控的重要前提。蒙特卡羅方法是目前準(zhǔn)確的吸收劑量計算方法，但存在收斂速度慢、計算時間長的缺點，阻礙了其在單板鈷源輻照裝置吸收劑量率計算中的應(yīng)用。MCNP (Monte Carlo N Particle Transport Code)是功能強大的蒙特卡羅程序，可模擬電子、γ和中子在復(fù)雜幾何中的輸運過程，其模擬準(zhǔn)度已經(jīng)通過諸多實驗的驗證，在輻射防護與輻照加工等領(lǐng)域發(fā)揮了重要的作用。文獻[1–3]利用MCNP計算了單板鈷源輻射場在不同介質(zhì)中的吸收劑量率，其中文獻[1]報道的模擬耗時為40 h。

模擬耗時表征模擬速度快慢，模擬精度表征模擬值的統(tǒng)計誤差，是蒙特卡羅方法固有的偏差。模擬耗時、模擬精度以及模擬值與實測值的相對偏差（簡稱相對偏差）可用于描述模擬計算效率。除問題的固有復(fù)雜程度外，MCNP模擬計算效率還與模擬參數(shù)相關(guān)，包括記數(shù)方法、記數(shù)柵元尺寸、截斷能以及γ致電子的處理方式等。(1) *F6和*F8是MCNP 中計算吸收劑量的主要記數(shù)方法[4]，*F8通過直接記錄柵元中γ和電子損失的能量給出吸收劑量，是MCNP中最接近吸收劑量定義的記數(shù)方法；*F6通過估算柵元的γ體通量給出比釋動能，其它一些記數(shù)方法，比如面通量(F2)、體通量(F4、FMESH)等也可以通過注量-劑量轉(zhuǎn)換法來計算比釋動能，但是其本質(zhì)都和*F6一樣，即通過γ注量與響應(yīng)函數(shù)的卷積計算比釋動能，在帶電粒子平衡(Charged Particle Equilibrium)條件下，比釋動能與吸收劑量比較接近；一般而言，*F8的耗時較長、相對偏差小，而在模擬精度相同時，*F6耗時較少；(2) 較大尺寸的記數(shù)柵元可以增大有效事件的概率，以減小模擬精度和柵元空間分辨率為代價，來提高精度、減小模擬耗時；(3) 電子在介質(zhì)中的頻繁碰撞會耗費大量的模擬時間，取較低的截斷能，可以保證模擬精度和準(zhǔn)度，但模擬耗時會大為增加。此外，MCNP5支持并行計算，可在保證模擬精度和相對偏差不變的條件下，大幅減小模擬耗時。

關(guān)于MCNP 模擬參數(shù)的選取沒有統(tǒng)一的標(biāo)準(zhǔn)，大多根據(jù)蒙特卡羅理論、實際需求和經(jīng)驗粗略選取參數(shù)，通過試算配合有限的實測值來調(diào)整參數(shù)，直至模擬值與實測值相符合后確定參數(shù)。邱有恒等[4]研究了柵元尺寸、δ電子、能量子步數(shù)和電子截斷能對計算效率的影響，建議使用*F6與*F8聯(lián)合記數(shù)，并提出了次級電子截斷能的自適應(yīng)方法來提高*F8的記數(shù)效率。林輝等[5]研究了頭部放療實例中MCNP截斷能、記數(shù)方法和γ致次級電子產(chǎn)生參數(shù)對*F8記數(shù)效率和速度的影響。目前國內(nèi)尚無關(guān)于提高鈷源輻照裝置吸收劑量率MCNP計算效率方面的報道。

本工作以某輻照站裝源量為8.4PBq的單板鈷源輻照裝置為對象，討論了MCNP中并行線程數(shù)、記數(shù)方法、記數(shù)柵元尺寸、截斷能和γ致次級電子處理方式對空場吸收劑量率模擬計算效率（模擬耗時、模擬精度以及模擬值與實測值的相對偏差）的影響，從而尋求在一定模擬精度和平均相對偏差前提下，使得模擬耗時最少的參數(shù)組合，為利用MCNP計算大、中型鈷源輻照裝置輻射場吸收劑量率提供參考數(shù)據(jù)。

1 輻照裝置與模擬過程

輻照裝置的設(shè)計裝源量為18.5PBq，最大可容納180根源棒。模擬時的裝源量為8.4PBq，根據(jù)劑量場均勻性的優(yōu)化設(shè)計，50根源棒[6]（活度分布見表1）分三層垂直排布在源架上。圖1給出了MCNP中幾何模型的示意圖，表2給出了主要結(jié)構(gòu)及其參數(shù)。模擬中的初始γ：(1) 位置，在單根源棒中均勻分布，不同源棒中的數(shù)量與活度成正比；(2) 能量，1.17MeV或1.33MeV；(3) 方向，4π各向同性。以矩形空氣塊作為記數(shù)柵元，分別采用*F8和*F6記數(shù)。整個模擬計算是在服務(wù)器（硬件：CPU，Intel Xeon X5675，3.07GHz；內(nèi)存：32GB；Windows Server 2003操作系統(tǒng)）上進行的，MCNP版本為MCNP5，并行計算程序為MPI version1.2.5。

表1 源棒活度分布Table 1 Activity distribution of source bars.

圖1 MCNP中幾何模型示意Fig.1 Scheme of geometry in MCNP.

表2 模擬中的主要結(jié)構(gòu)及參數(shù)Table 2 Major structures and their parameters used in simulation.

利用重鉻酸銀化學(xué)劑量計測量空間點的吸收劑量率[7]，為減小測量誤差對單個空間點連續(xù)進行3次測量后取平均值。本文的吸收劑量率均為空氣介質(zhì)吸收劑量率。實測值為水吸收劑量率，按文獻[8]方法換算為空氣吸收劑量率（模擬計算結(jié)果表明換算方法的誤差約為1.9%）。在距離源架30–115cm的常用輻照空間中，模擬和實測了10個空間點的吸收劑量率。

模擬精度σstat定義為多個空間點平均吸收劑量率X的統(tǒng)計誤差，計算式分別為：

式中，i為空間點的序號；Xi為i點處的吸收劑量率；D(X)為樣本方差；NPS為抽樣數(shù)；E(X)為平均吸收劑量率X的期望；E(Xi)為i點處吸收劑量率的期望；σstat–i為i點處吸收劑量率計算值的統(tǒng)計誤差；N為空間點的總數(shù)，為10。

平均相對偏差σm定義為多個空間點吸收劑量率實測值與模擬值相對偏差的均值，其計算式為：

式中，j為測量序號；P為測量次數(shù)，為3；Mi?j為i點處吸收劑量率的第j次測量值；Si為i點處吸收劑量率的模擬值，即式(2)中的E(Xi)；σm?i為i點處模擬值相對實測值的偏差。

討論步驟為：(1) 柵元尺寸對*F8和*F6計算效率的影響；(2) 并行線程數(shù)對*F6模擬耗時的影響；(3) 電子截斷能、γ致次級電子處理方式對*F6計算效率的影響；(4) 光子截斷能對*F6計算效率的影響。前一步確定的參數(shù)，用于后續(xù)步驟中的模擬。模擬中采取接續(xù)運行的方式，直到吸收劑量率峰值的統(tǒng)計誤差小于0.5%、σstat小于1%。平均相對偏差σm的限值為2%，模擬的粒子數(shù)NPS在2×108–2×109。

2 結(jié)果與討論

2.1 柵元尺寸對計算效率的影響

本工作以矩形空氣塊作為記數(shù)柵元，用邊長表征柵元的尺寸。表3給出了利用單核計算所得邊長L對計算耗時t、模擬精度σstat和平均相對偏差σm的影響。圖2給出了吸收劑量率模擬值和計算值的對比。由表3和圖2可知：(1) 對于*F8記數(shù)方法，由于空氣密度較低，對γ和電子的阻止能力弱，使得有效事件的概率低，邊長5 cm的柵元、2×109個粒子數(shù)可使σstat=0.39%滿足要求，但此時σm=16.43%仍不滿足要求，這是因為模擬值是整個柵元介質(zhì)吸收劑量率的平均值，若柵元過大，模擬值將不能用于準(zhǔn)確描述柵元中任意點的吸收劑量率，因此*F8記數(shù)方法不能同時獲得滿意的模擬精度和平均相對偏差，不適用于本文所討論的問題；(2) 對*F6記數(shù)方法，近似為常數(shù)，該結(jié)論可用于模擬參數(shù)的粗選；(3) 對*F6記數(shù)方法，計算速度與*F8相近(約3000photon.s?1)，但收斂速度比*F8快得多，采用邊長為0.5–5cm的柵元均可獲得模擬精度和平均相對偏差較好的模擬結(jié)果，0.5cm情況所需的模擬耗時約為1cm情況的5倍；(4) 記數(shù)方法一定程度決定了模擬計算的可行性，合適的柵元尺寸可在保證模擬精度和平均相對偏差前提下大幅減小模擬耗時，因此記數(shù)方法對計算效率的影響較大。

表3 不同柵元邊長下的MCNP計算效率Table 3 MCNP computational efficiency under different cell size.

圖2 吸收劑量率模擬值和計算值的對比 (a) 吸收劑量率值，(b) 模擬值相對測量值的偏差*F8和*F6記數(shù)柵元的邊長分別為5cm和1cm，(x,y)=(0,0)Fig.2 Contrast of simulated and measured absorbed dose rate, where (x,y)=(0,0). The size of *F8 and *F6 were 5 cm and 1 cm, respectively. (a) Value, (b) Relative deviation, measurements as benchmark

文獻[1?2]中柵元邊長的最小值均為1cm，對于單板鈷源輻照裝置1cm的空間分辨率足以較好地描述劑量率場的分布。綜上，后續(xù)討論中選用*F6記數(shù)方法、柵元邊長1cm作為模擬參數(shù)，模擬的粒子數(shù)均為3.9×108。

2.2 并行線程數(shù)對模擬耗時的影響

加速比S是衡量并行計算性能的重要指標(biāo)，其計算公式[9]為：

式中，n為參與并行計算的線程個數(shù)；t1為單線程計算所用時間，min；tn為n個線程計算所用時間，min。

一個程序通常含有并行部分以及無法并行化的串行部分，因此使用n個線程并行計算能達到的加速比只能≤n。英特爾處理器的超線程技術(shù)通過采用特殊硬件指令把一個物理線程模擬兩個邏輯線程，邏輯線程也可以并行工作。由于邏輯線程共享物理線程的資源，超線程模式下的加速比通常不能達到非超線程模式的兩倍。默認(rèn)情況下計算機操作系統(tǒng)開啟超線程模擬，用戶可通過設(shè)置BIOS進入非超線程模式。

MCNP程序的并行計算性能與問題的類型、復(fù)雜程度及參數(shù)設(shè)置等因素有關(guān)[10]。表4給出了本文所述問題并行線程數(shù)n與加速比S的關(guān)系。由表4可知，(1) n1/n2=2時，S1/S2為1.20–1.28，即超線程模式下的并行計算加速比約為非超線程模式的1.20–1.28倍；(2) 超線程模式下23個線程并行計算可獲得最大加速比13.1；(3) 若繼續(xù)增加并行線程數(shù)，可獲得加速比的持續(xù)增大，其對計算效率的影響大于柵元尺寸。后續(xù)討論中，均采用超線程模式下23個線程并行計算。

表4 并行計算線程數(shù)目與加速比的關(guān)系Table 4 Relationship between parallel threads and speed-up ratio.

2.3 電子截斷能、γ致電子的處理方式對計算效率的影響

*F6統(tǒng)計柵元中的γ，除γ與介質(zhì)的碰撞過程外，γ致電子與介質(zhì)碰撞過程也會產(chǎn)生次級γ，這些γ會影響*F6的統(tǒng)計結(jié)果。MCNP 中γ致電子有三種處理方式[11]：(1) γ-電子耦合輸運模式(即MODE P E)，電子正常輸運；(2) 電子致次級γ近似處理模式(即MODE P)，此時電子不輸運，而采用厚靶韌致輻射模型(Thick Target Bremsstrahlung)模式做近似處理；(3) γ輸運模式(即PHYS卡IDES=1)，電子也不輸運，產(chǎn)生后僅在“原地”沉積能量，無次級γ的產(chǎn)生，該模式對*F6統(tǒng)計結(jié)果無貢獻。在方式(1)中，采用較大的電子截斷能Ee-cut，可獲得模擬速度的提升，但會影響模擬值的平均相對偏差。

表5給出了不同電子截斷能和γ致電子的處理方式對計算效率的影響：情況(1) γ-電子耦合輸運模式，Ee-cut為1keV?1MeV，1keV是MCNP默認(rèn)值；情況(2)電子致次級γ近似處理模式；情況(3) γ輸運模式。根據(jù)文獻[5]，若兩條峰值誤差小于2%的曲線相比，其90%的劑量點的相對誤差小于3%，就可以認(rèn)為兩條曲線近似相等。圖3以γ-電子耦合輸運模式Ee-cut=1keV為基準(zhǔn)，對情況(1)、(2)和(3)的計算值進行了對比。由表5和圖3可知，(1) 三種情況下的模擬精度σstat一致，平均相對偏差σm的變化不大，說明電子截斷能、γ致電子的處理方式對模擬結(jié)果的影響可忽略，這是由于空氣密度較低，γ與之發(fā)生作用的概率小，使得γ致電子的產(chǎn)額低，其對*F6記數(shù)的貢獻也就小，因此改變γ致電子的處理方式不足以引起*F6統(tǒng)計結(jié)果明顯改變，相較而言，柵元尺寸對計算效率的影響更大；(2) 由于沒有γ致電子及其次級γ的輸運，情況(3)模擬耗時最少，與默認(rèn)值相比耗時減小約4倍。

綜上，后續(xù)選用PHYS卡IDES=1（即γ輸運模式）作為模擬參數(shù)。

表5 不同電子截斷能和γ致次級電子處理方式下計算效率Table 5 Computational efficiency under different electron cut-off energy and methods dealing with γ-induced electrons.

圖3 不同γ致電子處理方式下的相對偏差Fig.3 Relative deviation under different methods dealing with γ-induced electrons.

2.4 γ截斷能對計算效率的影響

表6給出了不同γ截斷能Eγ-cut對t、σstat和σm的影響。由表6可知，(1) 直到Eγ-cut=300keV，σstat和σm仍然滿足要求；(2) Eγ-cut=1keV時模擬值較實測值偏大，隨著γ截斷能增大，*F6計數(shù)率減小，σm出現(xiàn)先減小后增大的現(xiàn)象；(3) 光子截斷能的提高，對模擬速度的影響并不顯著，相較γ致電子的處理方式和電子截斷能而言，γ截斷能對計算效率的影響更小。

圖4以Eγ-cut=1keV (MCNP默認(rèn)值)為基準(zhǔn)，討論了不同γ截斷能對計算效率的影響。由圖4可知，(1) 隨著γ截斷能的增大，*F6模擬值逐漸減小，這是由于能量低于Eγ-cut的γ被“截斷”后對*F6沒有貢獻，增大Eγ-cut使得進入記數(shù)柵元γ的數(shù)量減??；(2) Eγ-cut相同時，隨著距離增大相對偏差逐漸變大，這是由于墻面散射γ的影響，距離越大越靠近墻面，散射γ比重越大；(3) 直到Eγ-cut=0.5MeV，90%劑量點的相對偏差都小于3%。

出于保守考慮，建議采用Eγ-cut=0.1MeV，與Eγ-cut默認(rèn)值相比，模擬耗時減小約1.1倍。

表6 不同γ截斷能下的計算效率Table 6 Computational efficiency under different effect of γ cut-off energy.

圖4 不同γ截斷能下的相對偏差Fig.4 Relative deviation under different effect of γ cut-off energy.

3 結(jié)語

對本工作，單板鈷源輻照裝置空場吸收劑量率MCNP模擬計算的總結(jié)：(1) 宜采用超線程模式下的并行計算，23個線程并行計算可將模擬耗時減小約13.1倍；(2) *F6記數(shù)方法適用于計算單板鈷源輻照裝置空場吸收劑量率，而*F8不適用；(3) 對于矩形空氣記數(shù)柵元，柵元邊長×度的值近似為常數(shù)(≈570)，邊長1cm的柵元，可在保證空間分辨率條件下獲得模擬精度和平均相對偏差滿意的模擬結(jié)果；(4) 相比γ截斷能而言，電子截斷能和γ致電子處理方式對計算效率的影響較大，采用γ輸運模式（即模擬參數(shù)PHYS卡IDES=1）、100keV的γ截斷能，可將模擬耗時進一步減小約4.4倍；(5) 不同因素對計算效率影響大小為：記數(shù)方法＞并行線程數(shù)＞柵元尺寸＞電子截斷能和γ致電子的處理方式＞γ截斷能。

綜上，本工作在保證模擬精度(1%)和相對偏差(2%)的條件下將MCNP模擬耗時減小約57倍，獲得了計算效率的提高。

1 郭平穩(wěn). Co-60輻照裝置劑量場分布的計算與測量研究[D]. 長沙: 國防科學(xué)技術(shù)大學(xué), 2004

GUO Pingwen. A study on computation an measurement of dose field distribution of Co-60 irradiation facility[D]. Changsha: National Defense Science and Technology University, 2004

2 劉江平.60Co輻照加工產(chǎn)品劑量分布的蒙特卡羅模擬[J]. 輻射研究與輻射工藝學(xué)報, 2010, 28(1): 57–61

LIU Jiangping. Monte Carlo simulation for the dose mapping of products irradiated by60Co γ-ray[J]. Journal of Radiation Research and Radiation Processing, 2010, 28(1): 57–61

3 Oliveira C, Ferreira L M, Gon-calves I F, et al. Monte Carlo studies of the irradiator geometry of the Portuguese gamma irradiation facility[J]. Radiation Physics and Chemistry, 2002, 65: 293–295

4 邱有恒, 應(yīng)陽君, 王敏, 等. 光子能量沉積計算的一種新方法[J]. 原子核物理評論, 2011, 28(1): 97–102

QIU Youheng, YING Yangjun, WANG Min, et al. Several techniques on increasing computation efficiency of photon energy deposition[J]. Nuclear Physics Review, 2011, 28(1): 97–102

5 林輝, 吳宜燦, 陳義學(xué). 基于臨床實例的影響蒙特卡羅模擬程序MCNP計算精度和速度的若干參數(shù)模擬研究[J]. 原子核物理評論, 2006, 23(2): 237–241

LIN Hui, WU Yican, CHEN Yixue. Investigation of parameters influencing on simulation precision and speed of Monte Carlo MCNP based on a clinic case[J]. Nuclear Physics Review, 2006, 23(2): 237–241

6 GB 7465-2009, 高活度鈷60密封放射源[S]. 2009 GB 7465-2009, High activity cobalt-60 sealed radioactive sources[S]. 2009

7 JJG 1028-1991, 使用重鉻酸銀劑量計測量γ射線水吸收劑量標(biāo)準(zhǔn)方法[S]. 1991

JJG 1028-1991, Standard method to measure γ-rays absorbed doses of water by sliver dichromate dosimeter[S]. 1991

8 GB/T15447-2008, X, γ射線和電子束輻照不同材料吸收劑量的換算方法[S]. 2008

GB/T15447-2008, Conversion method of absorbed doses in different material irradiated by X, γ rays and electron beams[S]. 2008

9 鄧力, 張文勇, 徐涵, 等. 蒙特卡羅程序MCNPIIMCNP-5并行效率比較[J]. 計算機工程與科學(xué), 2009, 31(A1): 185–187

DENG Li, ZHANG Wenyong, XU Han, et al. The comparison of parallel efficiency between Monte Carlo code MCNP-II and MCNP5[J]. Computer Engineering & Science, 2009, 31(A1): 185–187

10 王磊, 王侃, 余綱林. MCNP 程序并行計算性能分析[J]. 核科學(xué)與工程, 2006, 26(4): 301–306

WANG Lei, WANG Kan, YU Ganglin. Analysis of parallel computing performance of the code MCNP[J]. Chinese Journal of Nuclear Science and Engineering, 2006, 26(4): 301–306

11 X-5 Monte Carlo Team. MCNP-a general Monte Carlo n-particle transport code[R]. Version 5, LA-UR-03-1987, 2003

CLC TL929

Study on improving MCNP computational efficiency of dose rate for a single-rack cobalt irradiation facility

LI Lei LI Xiaoyan HUANG Wei JIANG Shubin WU Xiaoli YANG Guixia
(Institute of Nuclear Physics and Chemistry, Chinese Academy of Engineering Physics, Mianyang 621900, China)

Background: Monte Carlo method is regarded as the most accurate method for dose rate calculation at present, whereas its time-consuming defect increases complication for its application in industrial irradiation processing and experiment. Purpose: For a single-rack cobalt irradiation facility with loaded activity of 8.4PBq, a study of optimizing MCNP (Monte Carlo N Particle Transport Code) parameters was done to improve computational efficiency of dose rate. Methods: The effects of five parameters on MCNP computational efficiency, including the thread numbers of parallel computation, tally method, cell size, cutoff energy, and methods dealing with γ-induced electrons were discussed. Results: With certain precision and relative error, the least time-consuming situation was achieved with a set of five parameters, including parallel computation with hyper-thread mode, *F6 tally, cell size of 1cm, γ transport mode and γ-cutoff energy of 100keV. Conclusion: In this work, the cost of time was reduced about 57 times, and the improvement of MCNP computational efficiency was realized.

Single-rack cobalt irradiation facility, Absorbed dose rate, MCNP (Monte Carlo N Particle Transport Code), Computational efficiency

TL929

10.11889/j.0253-3219.2015.hjs.38.030201

李磊，男，1986年出生，2011年于四川大學(xué)獲理學(xué)碩士學(xué)位，助理研究員，從事抗輻射加固及輻照工藝研究

2014-03-10，

2014-07-29