NPTS程序質(zhì)子輸運(yùn)蒙特卡羅模擬方法研究

韓立會(huì),張華陽,衷 斌,成 立,沈華韻

(1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京 100094;2.中國工程物理研究院研究生院,北京 100085)

蒙特卡羅方法已經(jīng)成為解決不同物理問題的有效方法,在復(fù)雜模型粒子輸運(yùn)問題模擬方面具有顯著優(yōu)勢(shì)[1]。基于蒙特卡羅方法開發(fā)的質(zhì)子輸運(yùn)程序具有廣闊的應(yīng)用前景,例如,在質(zhì)子束放射治療中,基于醫(yī)用輻射劑量體模模擬肝臟質(zhì)子治療,得到覆蓋整個(gè)腫瘤的最佳質(zhì)子能量區(qū)間,并且可以模擬計(jì)算出二次輻射傳輸?shù)闹凶印⒐庾訉?duì)患者器官造成的非期望劑量,有助于優(yōu)化治療劑量以外的輻射劑量[2];在質(zhì)子加速器的屏蔽設(shè)計(jì)方面,蒙特卡羅方法通常被認(rèn)為是復(fù)雜加速器設(shè)計(jì)分析的最精確方法[3],如應(yīng)用蒙特卡羅程序模擬優(yōu)化質(zhì)子治療室的屏蔽材料及方案[4]、中子靶物理計(jì)算分析[5]、改進(jìn)束流裝置優(yōu)化質(zhì)子束性能[6];在太空輻射防護(hù)方面,空間輻射環(huán)境中占比最高的質(zhì)子不僅危害飛行機(jī)組人員健康,而且會(huì)損壞衛(wèi)星和高空飛行器中的微電子設(shè)備,導(dǎo)致嚴(yán)重故障,蒙特卡羅程序可以用來模擬質(zhì)子在圖像傳感器中產(chǎn)生的位移損傷,對(duì)其在軌安全運(yùn)行的評(píng)估具有重要的研究意義[7]。

國外蒙特卡羅程序發(fā)展較早,具備質(zhì)子輸運(yùn)模擬能力的程序主要包括PHITS[8]、MCNPX[9]、GEANT4[10]、FLUKA[11]等。國內(nèi)SuperMC采用Bertini核內(nèi)級(jí)聯(lián)模型開發(fā)了質(zhì)子輸運(yùn)模塊,可模擬能量范圍至10 GeV[12]。胡志良等開發(fā)了一款粒子核內(nèi)級(jí)聯(lián)蒙特卡羅模擬程序CBIM,可模擬45 MeV以上質(zhì)子核內(nèi)級(jí)聯(lián)過程[13]。Zhao等[14]基于INCL核內(nèi)級(jí)聯(lián)模型與ABLA磨損燒蝕模型在OpenMC上開發(fā)了質(zhì)子中子輸運(yùn)程序IMPC-MC1.0。目前國內(nèi)具備復(fù)雜幾何體中質(zhì)子輸運(yùn)模擬能力的程序較少,且基于核內(nèi)級(jí)聯(lián)模型、蒸發(fā)裂變模型開發(fā)的質(zhì)子輸運(yùn)模擬程序適用的能量范圍通常較高,在幾十MeV以上。

NPTS程序(neutron-photon transport simulation program)是由北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所開發(fā)的一款定常中子-光子耦合蒙特卡羅輸運(yùn)程序。本文應(yīng)用濃縮歷史方法基于帶電粒子輸運(yùn)連續(xù)慢化理論、多重散射理論與能量岐離分布等關(guān)鍵物理模型與ENDF/B-Ⅶ.0核數(shù)據(jù)庫開發(fā)質(zhì)子輸運(yùn)模擬程序,支持1 keV～150 MeV范圍內(nèi)質(zhì)子的輸運(yùn)模擬計(jì)算,實(shí)現(xiàn)質(zhì)子、中子、光子、電子多粒子耦合輸運(yùn)功能,有效拓展NPTS程序在帶電粒子輸運(yùn)領(lǐng)域的應(yīng)用,提高NPTS的應(yīng)用價(jià)值,同時(shí)豐富國內(nèi)中低能質(zhì)子輸運(yùn)程序的研究方法。

1 質(zhì)子輸運(yùn)蒙特卡羅模擬方法及物理模型

1.1 濃縮歷史方法

描述帶電粒子輸運(yùn)的玻爾茲曼方程為:

(1)

蒙特卡羅方法模擬帶電粒子輸運(yùn)計(jì)算,最大的問題在于帶電粒子的碰撞次數(shù)過于頻繁。為了減少計(jì)算量,Berger提出了濃縮歷史的方法[15]。Larsen[16]已證明濃縮歷史方法是足夠小步長條件下玻爾茲曼輸運(yùn)方程的近似解,本質(zhì)是人為的解耦了式(1)中的流、角散射和慢化的過程,將位置變化、角度偏轉(zhuǎn)和能量變化3個(gè)過程用一種相對(duì)簡單的蒙特卡羅方法來模擬。

對(duì)能量及步長的模擬,程序假設(shè)質(zhì)子每走一步,能量按對(duì)數(shù)減少:

(2)

式中,K=(1/2)1/m,表示質(zhì)子每走m步能量損失1/2,m為度量質(zhì)子游走步長的一個(gè)重要參數(shù),當(dāng)m為6時(shí),每個(gè)能量步對(duì)應(yīng)能量損失為10.9%;m為10時(shí),能量損失為6.7%。NPTS支持修改m,快速估算模擬結(jié)果時(shí)可將m設(shè)置稍小些,增加每個(gè)步長的能量損失,加快計(jì)算速度。

為提高模擬計(jì)算精度,在主步之間劃分子步,獲得更小的子步長,子步的能量間距為:

(3)

式中,ns為根據(jù)介質(zhì)材料劃分的子步數(shù)量,由經(jīng)驗(yàn)公式確定。

假設(shè)質(zhì)子沿軌跡的微分能量損失dE/dx是連續(xù)變化的函數(shù),對(duì)于給定能量損失ΔE的一個(gè)質(zhì)子所走的軌跡步長表示為:

(4)

NPTS應(yīng)用辛普森積分公式[17]得到與能量損失對(duì)應(yīng)的步長,并通過步長累加得到質(zhì)子的連續(xù)慢化射程。

1.2 帶電粒子的能量損失

具有一定動(dòng)能的質(zhì)子入射到靶物質(zhì)中會(huì)與原子核及核外電子發(fā)生庫侖相互作用,發(fā)生動(dòng)量轉(zhuǎn)移與電荷交換,通過連續(xù)小能量轉(zhuǎn)移,入射粒子逐漸損失能量,速度變慢。慢化過程中,電子阻止是造成質(zhì)子能量損失的主要原因,在能量較低(E/A<10 keV)時(shí)需要考慮核阻止。

1) 電子阻止

本文應(yīng)用帶有殼層效應(yīng)修正[18-21]及密度效應(yīng)修正[22]的Bethe公式計(jì)算入射帶電粒子的電子阻止,表達(dá)形式[23]如下:

(5)

式中:ρ為材料密度;re為電子經(jīng)典半徑;u為原子質(zhì)量單位,即碳原子重量的1/12;β為相對(duì)光速速度;fi為元素的原子組分;Ai為元素的原子質(zhì)量;Ii為元素的平均電離能;z為入射粒子質(zhì)量數(shù);m為入射粒子質(zhì)量;ci為殼層效應(yīng)修正系數(shù);δ為密度效應(yīng)校正系數(shù);εmax為質(zhì)量為m的粒子與一個(gè)質(zhì)量為me的軌道電子(可以認(rèn)為是自由電子)碰撞過程中可能轉(zhuǎn)移的最大能量:

(6)

由于Bethe公式推導(dǎo)過程中采用一階玻恩近似,在限制條件β>0.05z2/3范圍內(nèi)準(zhǔn)確性較高,對(duì)應(yīng)于質(zhì)子能量約為1.17 MeV,對(duì)于其他重帶電粒子如氘粒子、α粒子則分別對(duì)應(yīng)2.34、4.65 MeV[24]。

低能部分程序基于Lindhard理論計(jì)算電子阻止,電子阻止公式[25]表示為:

(7)

隨著入射質(zhì)子能量降低速度減小,會(huì)從介質(zhì)原子中獲得電子,從而導(dǎo)致入射粒子有效電荷發(fā)生改變,使用Barkas經(jīng)驗(yàn)公式計(jì)算有效電荷[27]:

z*=z(1-e-125βz-2/3)

(8)

當(dāng)質(zhì)子能量降至NPTS模擬的截止能量1 keV時(shí)對(duì)應(yīng)有效電荷為0.087 5個(gè)正電荷。

2) 核阻止

對(duì)于具有非常小比能量的入射粒子,特別是重帶電粒子入射到原子序數(shù)較大的靶材料中,由于粒子與靶原子核之間的彈性庫侖碰撞導(dǎo)致的粒子減速變得顯著。Lindhard等[28]給出了控制這一過程的通用曲線,Steward對(duì)該曲線進(jìn)行最小二乘法擬合得到核阻止計(jì)算公式[29]。

(9)

3) 阻止本領(lǐng)計(jì)算結(jié)果分析

SRIM是一個(gè)計(jì)算離子在物質(zhì)中能量損失與射程的程序包,自1985年推出以來被廣泛應(yīng)用[30]。分別使用SRIM程序和NPTS程序計(jì)算質(zhì)子在水、銅、鉛3種常見材料中的電子阻止本領(lǐng)與核阻止本領(lǐng),對(duì)比結(jié)果如圖1所示。能量高于1 MeV時(shí),不同材料與SRIM計(jì)算結(jié)果的相對(duì)偏差在2%以內(nèi)。當(dāng)NPTS與SRIM計(jì)算結(jié)果相對(duì)偏差在5%以內(nèi)時(shí),對(duì)應(yīng)于水、銅、鉛的質(zhì)子能量分別為0.7、0.225、0.250 MeV。

圖1 質(zhì)子在多種材料中的阻止本領(lǐng)Fig.1 Stopping power of proton in different materials

針對(duì)驗(yàn)證的3種材料,如圖1a所示,NPTS與SRIM在低能段的電子阻止本領(lǐng)偏差較大,特別是對(duì)于水材料,Vagena等[31]對(duì)GEANT4、MCNP、PHITS、SRIM等程序計(jì)算的阻止本領(lǐng)比結(jié)果顯示,不同程序計(jì)算1 MeV質(zhì)子在水材料中的阻止本領(lǐng)偏差為5%左右,NPTS與SRIM對(duì)于入射能量在1 MeV的質(zhì)子阻止本領(lǐng)偏差為1%。NPTS處理非單質(zhì)介質(zhì)的能量損失時(shí)基于布拉格相加法則,SRIM阻止本領(lǐng)計(jì)算更多考慮了分子結(jié)構(gòu),在布拉格法則與含有H、C、N、O等元素的化合物之間進(jìn)行了校正[30]。通過與SRIM結(jié)果對(duì)比,建議NPTS使用能量范圍大于1 MeV,尤其對(duì)化合物阻止本領(lǐng)的計(jì)算。圖1b為核阻止對(duì)比結(jié)果,因?yàn)楹俗柚怪辉诘湍芮闆r下對(duì)總的能量損失有貢獻(xiàn),所以本文僅對(duì)比1 MeV以下兩個(gè)程序的核阻止本領(lǐng)的計(jì)算結(jié)果,3種材料核阻止的計(jì)算結(jié)果相對(duì)偏差均在8%以內(nèi),對(duì)于鉛這種原子序數(shù)大的靶材料,核阻止相對(duì)偏差在2%以內(nèi)。

此外同樣對(duì)比了氘粒子與α粒子在3種材料中總的阻止本領(lǐng),如圖2所示,當(dāng)總的阻止本領(lǐng)相對(duì)偏差控制在5%以內(nèi)時(shí)對(duì)應(yīng)的氘粒子與α粒子能量分別為2.5、5.5 MeV。以上對(duì)比結(jié)果表明NPTS在較大能量范圍內(nèi)與SRIM吻合較好,可用于計(jì)算帶電粒子在介質(zhì)中的平均能量損失以及相應(yīng)的輸運(yùn)計(jì)算。

圖2 氘粒子、α粒子在多種材料中總的阻止本領(lǐng)Fig.2 Stopping power of D and Alpha particles in different materials

1.3 多重散射

當(dāng)帶電粒子入射到靶物質(zhì)中時(shí),與原子核及核外電子發(fā)生大量的碰撞從而改變粒子的方向,應(yīng)用多次散射理論得到帶電粒子在軌跡末端產(chǎn)生的總偏轉(zhuǎn)角分布。NPTS應(yīng)用Striganov修正的Molière理論,分布形式如下:

F(θ,t)=

(10)

式中:z為入射粒子質(zhì)子數(shù);p為動(dòng)量;N為阿伏伽德羅常數(shù);A為靶核原子質(zhì)量;參數(shù)B、χmax的定義見文獻(xiàn)[32]。

Striganov模型和Molière的原始理論之間采用的微分散射截面計(jì)算不同,在Molière的理論中,微分散射截面由盧瑟福公式推導(dǎo)得出,帶有屏蔽參數(shù):

(11)

(12)

在Striganov模型中,微分散射截面結(jié)合粒子與原子核有限尺寸進(jìn)一步修正得到:

(13)

1.4 能量岐離

電子阻止與核阻止反映了質(zhì)子在單位距離內(nèi)的平均能量損失,而每一個(gè)能量步代表多次獨(dú)立的隨機(jī)碰撞的累積效應(yīng)。質(zhì)子在通過同樣長度的路程內(nèi)所損失的能量是不同的,存在一個(gè)概率分布,對(duì)于質(zhì)子、氘粒子等重帶電粒子應(yīng)用Vavilov分布描述,形式如下:

(14)

由于Vavilov解析解的復(fù)雜性,計(jì)算耗時(shí)嚴(yán)重,特別是在蒙特卡羅方法模擬帶電粒子輸運(yùn)過程,程序中實(shí)際應(yīng)用Rotondi和Montagna開發(fā)的近似模型以便快速計(jì)算Vavilov分布[34]。該模型根據(jù)κ選擇適當(dāng)?shù)臄M合函數(shù)對(duì)Vavilov進(jìn)行分段擬合。當(dāng)κ低于0.01時(shí),Vavilov分布近似使用朗道分布描述,對(duì)于高于12的κ,近似應(yīng)用高斯分布。

1.5 質(zhì)子核數(shù)據(jù)庫

質(zhì)子核反應(yīng)采用洛斯阿拉莫斯實(shí)驗(yàn)室開發(fā)的質(zhì)子核評(píng)價(jià)數(shù)據(jù)庫LA150H。LA150H數(shù)據(jù)評(píng)價(jià)使用核反應(yīng)數(shù)據(jù)計(jì)算程序GNASH進(jìn)行理論分析,并以實(shí)驗(yàn)數(shù)據(jù)為基準(zhǔn),數(shù)據(jù)為ACE(a compact ENDF)格式[35]。該數(shù)據(jù)庫對(duì)重要靶材料和屏蔽同位素的評(píng)估延伸至150 MeV,包含41種核素的評(píng)價(jià)數(shù)據(jù)。ENDF/B-Ⅶ.0中新增了3H(p,n)、3He(p,n)、6Li(p,n)、7Li(p,n)、9Be(p,n)、10B(p,n)、13C(p,n)的反應(yīng)數(shù)據(jù)[36],ENDF/B-Ⅷ.0增加了4He(p,n)反應(yīng)數(shù)據(jù),目前ENDF質(zhì)子子庫包含49個(gè)核素的LANL評(píng)估數(shù)據(jù)。質(zhì)子核數(shù)據(jù)庫支持連續(xù)能量中子、光子及帶電粒子耦合計(jì)算,包含次級(jí)輕粒子以及重反沖粒子的產(chǎn)生截面和能量角度分布。NPTS程序基于核數(shù)據(jù)庫精細(xì)描述帶電粒子核反應(yīng)及次級(jí)中子、伽馬及重帶電離子的產(chǎn)生過程。

2 程序開發(fā)

程序輸入文件采用XML語言格式,輸入文件以卡片形式書寫,主要分柵元卡、曲面卡、材料卡和物理卡四大類,其中物理卡含控制卡、計(jì)數(shù)卡、加速卡和源卡等描述模型物理相關(guān)參數(shù)的卡片。開始模擬質(zhì)子輸運(yùn)前,程序已基于控制卡中模擬質(zhì)子的能量信息計(jì)算并儲(chǔ)存了能量網(wǎng)格、平均能量損失、子步長、射程等信息。

單個(gè)質(zhì)子輸運(yùn)過程如圖3所示,NPTS模擬步驟如下。

圖3 質(zhì)子輸運(yùn)程序流程圖Fig.3 Flow chart of proton transport program

1) 首先進(jìn)入柵元循環(huán),得到質(zhì)子所在柵元填充的材料信息,根據(jù)經(jīng)驗(yàn)公式確定內(nèi)循環(huán)子步數(shù)。

2) 進(jìn)入能量循環(huán),根據(jù)質(zhì)子能量與靶材料信息查找步長信息得到子步長ds。

3) 進(jìn)入最內(nèi)層子步循環(huán)中,通過電子阻止與核阻止計(jì)算能量損失,根據(jù)損失能量前后的能量平均值查找已處理的截面信息并計(jì)算反應(yīng)距離dr,計(jì)算質(zhì)子到柵元邊界的距離db,確認(rèn)ds、dr、db中最小值為質(zhì)子輸運(yùn)距離,更新質(zhì)子位置。根據(jù)多重散射理論角分布計(jì)算質(zhì)子偏轉(zhuǎn)角度更新質(zhì)子方向。根據(jù)質(zhì)子輸運(yùn)距離判斷模擬過程,如果db最小,根據(jù)邊界條件及多重散射偏轉(zhuǎn)角度確認(rèn)質(zhì)子是否穿出當(dāng)前柵元,歷史結(jié)束或者進(jìn)入下一個(gè)柵元循環(huán);如果dr最小,則基于核數(shù)據(jù)庫進(jìn)行核反應(yīng)及次級(jí)粒子的模擬;如果ds最小,質(zhì)子能量網(wǎng)格發(fā)生變化且未達(dá)到截止能量則更新能量網(wǎng)格索引后重新進(jìn)入能量循環(huán),反之進(jìn)入下一次子步循環(huán)。

3 結(jié)果分析

3.1 中子產(chǎn)額分析計(jì)算

應(yīng)用NPTS程序模擬了不同能量(5、6、7、8、9和10 MeV)質(zhì)子入射9Be厚靶模型產(chǎn)生的次級(jí)中子情況(模擬入射質(zhì)子數(shù)為108,NPTS程序的漲落在0.5%以內(nèi)),與文獻(xiàn)[37]中的GEANT4、MCNP6模擬數(shù)據(jù)及Atta的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比。中子產(chǎn)額隨入射質(zhì)子能量的變化如圖4a所示,不同蒙特卡羅程序模擬得到的中子產(chǎn)額與Atta實(shí)驗(yàn)值的相對(duì)偏差列于表1。5 MeV時(shí)與實(shí)驗(yàn)值偏差較大,此時(shí)蒙特卡羅程序模擬結(jié)果相對(duì)偏差均在10%以上,NPTS與GEANT4、MCNP6模擬結(jié)果較為一致。圖4b為產(chǎn)生的次級(jí)中子的平均能量,出射中子的平均能量近似與入射質(zhì)子能量呈正比,NPTS與GEANT4模擬結(jié)果的相對(duì)偏差如表1所示,均小于2%。

表1 質(zhì)子轟擊9Be厚靶中子產(chǎn)額Table 1 Thick target neutron yields for 9Be metal

圖4 5～10 MeV質(zhì)子轟擊厚9Be靶產(chǎn)生的次級(jí)中子信息Fig.4 Neutron information of 5-10 MeV proton bombarding thick 9Be target

BNCT治療應(yīng)用的中子能量范圍通常在10 keV以下,而出射中子的平均能量與入射質(zhì)子的能量近似呈正比,在低能情況下鋰靶的中子產(chǎn)額、中子出射角度、中子能量方面相比鈹靶更適合作為BNCT的靶材料。使用NPTS對(duì)能量在7Li(p,n)7Be反應(yīng)閾值附近的質(zhì)子入射7Li厚靶模型進(jìn)行計(jì)算(模擬粒子數(shù)為109,此時(shí)對(duì)于不同算例NPTS程序的漲落在0.16%以內(nèi))。以Lee等[38]的研究數(shù)據(jù)為參考,對(duì)比NPTS的計(jì)算結(jié)果,如表2所列。中子產(chǎn)額NPTS的計(jì)算結(jié)果較Lee等的計(jì)算值在1.96 MeV左右最為接近,相對(duì)偏差整體在4%以內(nèi),出射中子的平均能量與Lee等的數(shù)據(jù)吻合較好,相對(duì)偏差小于1.5%,出射中子平均角度余弦,相比Lee等計(jì)算值偏差稍大,但基本在5%以內(nèi)。

表2 質(zhì)子轟擊7Li厚靶中子產(chǎn)額Table 2 Thick target neutron yields for 7Li metal

應(yīng)用NPTS程序模擬不同能量質(zhì)子入射鋰、鈹厚靶模型,通過與GEANT4、MCNP6及已有研究數(shù)據(jù)的對(duì)比,結(jié)果表明NPTS與主流蒙特卡羅模擬程序及實(shí)驗(yàn)值具有較好的一致性,驗(yàn)證了NPTS程序在模擬厚靶中子產(chǎn)額方面的正確性。

質(zhì)子加速器靶的設(shè)計(jì)過程中需要明確靶的厚度及幾何形狀,以獲得最優(yōu)的中子束流,應(yīng)用NPTS測(cè)試靶厚對(duì)中子產(chǎn)額的影響,計(jì)算了不同厚度7Li、9Be靶的中子產(chǎn)額,如圖5a、b所示,在一定質(zhì)子能量下,隨著厚度的增加,中子產(chǎn)額曲線會(huì)出現(xiàn)飽和,達(dá)到穩(wěn)定,達(dá)到穩(wěn)定的厚度近似等于慢化到產(chǎn)生中子反應(yīng)閾能的射程,對(duì)于鋰材料反應(yīng)閾能為1.88 MeV,鈹材料的反應(yīng)閾能為2.06 MeV。圖6為NPTS計(jì)算得到的能量在100 MeV以下的質(zhì)子在7Li、9Be靶當(dāng)中能量降低到7Li(p,n),9Be(p,n)反應(yīng)閾能時(shí)的射程,根據(jù)圖6b得到的最佳靶厚與圖5a、b程序模擬計(jì)算的中子產(chǎn)額最大值時(shí)的靶厚非常接近。

圖5 中子產(chǎn)額隨靶厚度的變化Fig.5 Neutron yields of different thicknesses of target

圖6 質(zhì)子連續(xù)慢化射程曲線Fig.6 Proton continuous slowing down approximation range

100 MeV以下質(zhì)子在7Li、9Be靶中慢化到反應(yīng)閾能以下的質(zhì)子射程如圖6a所示,根據(jù)NPTS給出的圖6曲線得到30 MeV質(zhì)子入射鈹靶對(duì)應(yīng)中子產(chǎn)額的最佳厚度為0.6 cm,這與Dorostkar等[39]測(cè)試得到的結(jié)果相一致。根據(jù)圖6給出的能量射程曲線可以快速估算出中子產(chǎn)額最大時(shí)對(duì)應(yīng)的靶厚,為BNCT靶設(shè)計(jì)提供有利參考。此外在應(yīng)用蒙特卡羅程序模擬計(jì)算中子產(chǎn)額問題過程中,建議靈活設(shè)置質(zhì)子截止能量,因?yàn)榈陀谥凶臃磻?yīng)閾能的質(zhì)子對(duì)中子產(chǎn)額沒有貢獻(xiàn),但計(jì)算低能質(zhì)子輸運(yùn)耗時(shí)相當(dāng)嚴(yán)重,合理設(shè)置質(zhì)子截止能量可以有效地加快計(jì)算速度。

3.2 典型質(zhì)子治療室輻射場分布

質(zhì)子治療通常使用能量70～230 MeV之間的質(zhì)子束,這個(gè)能量范圍的質(zhì)子會(huì)產(chǎn)生高能中子、光子和其他帶電粒子,其中產(chǎn)生的次級(jí)中子在劑量方面占主導(dǎo)地位[40]。典型質(zhì)子治療室結(jié)構(gòu)如圖7所示[41],內(nèi)部空間尺寸為6.5 m×6.5 m,高度為3 m,周圍的混凝土墻壁厚度為1.5 m,從入口到出口的長度約9.5 m。1束150 MeV的單能質(zhì)子束從中心沿y軸正方向發(fā)射,撞擊位于治療室中心的圓柱形銅靶(直徑7 cm,長度7 cm,軸線為y軸),在該區(qū)域創(chuàng)建一個(gè)輻射場。GEANT4與NPTS計(jì)算的質(zhì)子治療室中子三維輻射場分布如圖8所示,圖中,橫縱坐標(biāo)為模擬計(jì)算劃分的網(wǎng)格數(shù)量,對(duì)應(yīng)1個(gè)網(wǎng)格xy方向尺寸為5 cm×5 cm。模擬粒子數(shù)均為1千萬,NPTS治療室內(nèi)部通量的平均漲落為0.5%,由于屏蔽墻體計(jì)數(shù)率低,漲落較大,模型整體平均漲落為7.1%。圖8顯示兩程序銅靶外中子通量分布吻合較好,屏蔽墻體外側(cè)通量計(jì)數(shù)較低,但從分布形式來看較為一致,驗(yàn)證了NPTS程序在計(jì)算復(fù)雜幾何體中質(zhì)子、中子輻射場分布的準(zhǔn)確性。

圖7 典型質(zhì)子治療室布置Fig.7 Layout of typical proton treatment room

圖8 質(zhì)子治療室內(nèi)中子輻射場分布Fig.8 Neutron radiation field distribution in proton treatment room

4 結(jié)論

本文基于NPTS程序,應(yīng)用帶電粒子輸運(yùn)的關(guān)鍵物理模型及ACE格式的核數(shù)據(jù)庫,實(shí)現(xiàn)了帶電粒子的輸運(yùn)與核反應(yīng)模擬,開發(fā)了適用范圍在1 keV～150 MeV的質(zhì)子輸運(yùn)蒙特卡羅模擬程序,并建立了重帶電粒子蒙特卡羅模擬方法與框架。應(yīng)用NPTS程序?qū)υ诔Ｒ姴牧现匈|(zhì)子、氘核、α粒子的阻止本領(lǐng)與SRIM進(jìn)行了對(duì)比分析,證明了NPTS在計(jì)算帶電粒子阻止本領(lǐng)方面的正確性。針對(duì)厚靶模型得到的中子產(chǎn)額、平均能量、角分布等數(shù)值結(jié)果與GEANT4、MCNP6及實(shí)驗(yàn)結(jié)果吻合較好,典型質(zhì)子治療室三維輻射場分布與GEANT4一致,驗(yàn)證了NPTS程序質(zhì)子、中子耦合輸運(yùn)功能的準(zhǔn)確性,有效擴(kuò)展了NPTS程序在質(zhì)子治療、加速器設(shè)計(jì)、屏蔽設(shè)計(jì)、輻射防護(hù)等方面的應(yīng)用。此外,NPTS程序模擬質(zhì)子輸運(yùn)應(yīng)用的物理模型同樣適用于其他重帶電粒子的輸運(yùn)計(jì)算,為進(jìn)一步研究其他重帶電粒子的輸運(yùn)模擬提供了程序基礎(chǔ)。