周斌 于全芝 胡志良 陳亮 張雪熒 梁天驕

1) (散裂中子源科學(xué)中心, 東莞 523803)

2) (中國科學(xué)院物理研究所, 北京 100190)

3) (中國科學(xué)院近代物理研究所, 蘭州 730000)

高能質(zhì)子在散裂靶中的能量沉積是散裂靶中子學(xué)研究的重要內(nèi)容之一, 準確掌握高能質(zhì)子在散裂靶中引起的能量沉積分布與瞬態(tài)變化是開展散裂靶熱工流體設(shè)計的重要前提.本文采用MCNPX, PHITS 與FLUKA 三種蒙特卡羅模擬程序, 計算并比較了高能質(zhì)子入射重金屬鉛靶、鎢靶的能量沉積分布及不同粒子對總能量沉積的占比貢獻; 針對高能質(zhì)子入射金屬鎢靶的能量沉積實驗數(shù)據(jù)空白, 采用熱釋光探測器陣列測量了250 MeV 質(zhì)子束入射厚鎢靶的能量沉積分布, 實驗結(jié)果表明蒙特卡羅模擬程序在散裂靶中能量沉積的計算結(jié)果具有較高的可靠性.

1 引 言

近年來, 我國先后啟動了中國散裂中子源(China Spallation Neutron Source, CSNS)項目[1],以及未來先進核裂變能-加速器驅(qū)動次臨界系統(tǒng)(accelerator driven sub-critical system, ADS)項目[2].CSNS 項目采用1.6 GeV 的高能質(zhì)子束流轟擊鉭包鎢靶[3], 在靶上發(fā)生散裂反應(yīng)并產(chǎn)生大量中子, 為中子散射及相關(guān)領(lǐng)域提供科學(xué)研究與實驗平臺; 在ADS 項目中, 高能質(zhì)子轟擊重金屬靶, 為ADS 次臨界反應(yīng)堆提供持續(xù)中子外源.在散裂靶中, 高能質(zhì)子主要通過與靶原子的激發(fā)電離過程,以及與靶核發(fā)生散裂反應(yīng)過程損失能量.掌握高能質(zhì)子在散裂靶中能量沉積的產(chǎn)生機制, 分析熱量沉積的空間分布與瞬態(tài)變化, 是開展散裂靶熱工流體設(shè)計的重要前提.

目前, 高能質(zhì)子束在散裂靶中的能量沉積數(shù)據(jù)主要通過MCNPX[4], PHITS[5], FLUKA[6]等蒙特卡羅粒子輸運程序計算獲得, 模擬計算的準確度依賴于程序中物理模型對各個能量沉積途徑的處理及采用的計算方法.MCNPX 程序中的核子級聯(lián)模型包括BERTINI, ISABEL, CEM2K 和INCL4,各個模型的主要區(qū)別在于對級聯(lián)相互作用過程與對核密度等的近似處理; PHITS 采用JAM 模型作為強子級聯(lián)模型, 在低能端采用JENDL 核數(shù)據(jù);FLUKA 采用PEANUT 模型與ENDF 核數(shù)據(jù)進行相關(guān)計算.

受實驗條件如束流能量、束流時間等的限制,高能質(zhì)子在散裂靶中的能量沉積測量數(shù)據(jù)非常有限.國際上采用量熱法進行靶內(nèi)能量沉積測量[7?12],即通過測量靶內(nèi)各處由于能量沉積導(dǎo)致的溫度變化而獲得能量沉積分布, 但是只有束流較強時才能使靶體溫度上升到可準確測量的范圍.文獻[13]表明在較低束流強度條件下, 在散裂靶內(nèi)布置熱釋光探測器(thermoluminescence detector, TLD), 通過測量TLD 的劑量分布可以分析散裂靶內(nèi)的能量沉積.本文采用MCNPX, PHITS 與FLUKA 三種蒙特卡羅程序, 首先計算了不同能量高能質(zhì)子入射重金屬鉛靶、鎢靶的總能量沉積及其在入射深度的分布, 并通過與已有的測量數(shù)據(jù)進行比較, 評估不同級聯(lián)模型的計算結(jié)果; 其后計算了不同能量的質(zhì)子在散裂靶中產(chǎn)生的不同類型粒子, 并獲得各種粒子對總能量沉積的占比貢獻; 最后利用蘭州重離子加速器[14]提供的250 MeV 質(zhì)子束流, 采用TLD陣列測量250 MeV 質(zhì)子束入射鎢靶的能量沉積分布, 測量結(jié)果驗證了蒙特卡羅程序計算質(zhì)子入射散裂靶能量沉積的可靠性.

2 能量沉積模擬計算

2.1 計算方法

高能質(zhì)子在散裂靶中的能量沉積計算分別采用MCNPX, PHITS 和FLUKA 三種粒子輸運程序開展.在采用MCNPX 程序進行計算時, 分別采用了BERTINI 模型、ISABEL 模型、CEM2K 模型與INCL4 模型作為級聯(lián)模型.在采用PHITS 程序和FLUKA 程序進行計算時, 則分別采用了JAM級聯(lián)模型與PEANUT 模型.計算所選用的散裂靶模型為直徑為20 cm, 長度為50 cm 的圓柱體.鉛的密度為11.3 g/cm3, 同位素204Pb,206Pb,207Pb和208Pb 的質(zhì)量豐度分別為1.4%, 24.1%, 22.1%和52.4%; 鎢的密度為19.0 g/cm3, 同位素180W,182W,183W,184W 和186W 的質(zhì)量豐度分別為0.12%,26.5%, 14.31%, 30.64%和28.43%.為了對比于國際上已有的實驗數(shù)據(jù), 驗證模擬計算方法的可靠性, 質(zhì)子能量分別為800, 1000, 1200 MeV, 從圓柱體靶端面垂直入射, 質(zhì)子在入射橫截面上的抽樣概率服從二維高斯分布, 半高寬均為2.4 cm.

2.2 質(zhì)子入射鉛靶的能量沉積計算結(jié)果

圖1(a)給出了MCNPX, PHITS 與FLUKA 三種粒子輸運程序?qū)θ肷淠芰繛?00, 1000, 1200 MeV的質(zhì)子入射鉛靶的總能量沉積.同時, 圖1(a)還給出了Belyakov-Bodin 等[7]采用熱電偶法在鉛靶中測量到的能量沉積實驗值.可以看出, 三種粒子輸運程序的計算結(jié)果相對于實驗值整體偏大; 與MCNPX 相比, PHITS 與FLUKA 的 計算結(jié)果更大;即使采用MCNPX 程序進行計算, 不同的級聯(lián)模型所對應(yīng)的模擬結(jié)果也各不相同: 采用ISABEL模型的計算值較大, INCL4 與BERTINI 模型次之, CEM2K 模型的計算結(jié)果最小.對1000 MeV的質(zhì)子入射能量而言, PHITS 的模擬值比測量值高約16.1%, FLUKA 的模擬值比測量值高15.6%,MCNPX-CEM2K 的模擬值比測量值高約9%, 采用CEM2K 級聯(lián)模型獲得的能量沉積值最接近實驗測量值, 計算結(jié)果的統(tǒng)計誤差可忽略不計.

圖1 (a)不同模擬程序?qū)︺U靶總能量沉積計算的對比;(b) 鉛靶中能量沉積線性密度的軸向分布Fig.1.(a) Comparison of total energy deposition in lead target calculated by different Monte Carlo codes; (b) axial distribution of linear density of energy deposition in lead target.

散裂靶單位長度中的能量沉積值被稱為能量沉積線性密度.圖1(b)給出了采用MCNPX-CEM-2K 模型計算出的能量沉積線性密度在鉛靶中的分布.通過與文獻[7]中的測量數(shù)據(jù)進行比較, 可以看出, 高能質(zhì)子在鉛靶長度方向的能量沉積線性密度計算結(jié)果與測量結(jié)果在整體上能較好符合; 隨著質(zhì)子入射深度的增加, 能量沉積線性密度計算值的偏差略有增大.

另外, 我們采用MCNPX 程序中的CEM2K級聯(lián)模型計算得到了不同能量質(zhì)子入射鉛靶時, 不同粒子對總能量沉積值的占比貢獻, 如表1 所示.計算數(shù)據(jù)表明, 質(zhì)子作為入射粒子, 通過與靶核外電子的電離、激發(fā)過程及與靶核的散裂反應(yīng)過程損失能量, 對總能量沉積的貢獻最大; 鉛是原子序數(shù)較高的材料, 對光子的阻止本領(lǐng)很強, 這導(dǎo)致了質(zhì)子與散裂靶相互作用產(chǎn)生的光子在靶中的能量沉積較大; 隨著入射質(zhì)子能量的增大, 散裂反應(yīng)產(chǎn)生的次級帶電粒子增多, 預(yù)平衡與蒸發(fā)過程的剩余核會達到更高的自旋和激發(fā)態(tài), 這導(dǎo)致了質(zhì)子對總能量沉積的占比逐漸減小, 輕帶電粒子與光子的占比貢獻逐步增大.

表1 CEM2K 級聯(lián)模型計算質(zhì)子入射鉛靶產(chǎn)生的不同粒子對總能量沉積的占比貢獻Table 1.The calculated contribution of different particles to the total energy deposition in lead target by CEM2K-Cascade-Mode.

為了對比不同級聯(lián)模型的次級粒子計算結(jié)果,我們還采用MCNPX 中的BERTINI, ISABEL,INCL4 級聯(lián)模型分別對1000 MeV 質(zhì)子進行了模擬, 計算結(jié)果如表2.從表2 可以看出, 4 種級聯(lián)模型對初級質(zhì)子電離作用的能量沉積計算結(jié)果幾乎相等, 這是因為初級質(zhì)子的電離作用僅受到入射粒子和散裂靶的影響, 與級聯(lián)過程無關(guān); 相對其他兩種級聯(lián)模型, CEM2K 對次級質(zhì)子沉積能量計算值最小, 對光子、氘、氚、氦-3 的能量沉積計算值最大, 這是因為CEM2K 融合了DUBNA 級聯(lián)模型、Exciton 激子模型、GEM 蒸發(fā)模型、RAL 裂變模型、Fermi-breakup 模型及光核反應(yīng)模型等, 在計算剩余核、裂變產(chǎn)物、輕核產(chǎn)生及碰撞-反彈過程時更為準確[15].

2.3 質(zhì)子入射鎢靶的能量沉積計算結(jié)果

同樣, 對不同能量的質(zhì)子入射鎢靶產(chǎn)生的能量沉積進行計算, 圖2(a)為采用不同計算模型獲得的總能量沉積隨入射質(zhì)子能量的變化.與鉛靶的情況類似, PHITS 與FLUKA 的計算結(jié)果相對于MCNPX 各個級聯(lián)模型偏大.圖2(b)為采用MCNPXCEM2K 模型計算的能量沉積線性密度分布.在現(xiàn)階段, CSNS 項目采用了鉭包鎢作為散裂靶, ADS項目則采用鎢鎳合金靶, 然而, 由于國內(nèi)外尚無高能質(zhì)子入射鎢靶產(chǎn)生的能量沉積實驗數(shù)據(jù), CSNS與ADS 的熱工流體設(shè)計只能根據(jù)計算數(shù)據(jù)開展.對于中高能質(zhì)子入射鎢靶的能量沉積計算亟需實驗數(shù)據(jù)驗證和可靠程度評估.

3 質(zhì)子在鎢靶中能量沉積的實驗驗證

3.1 實 驗

目前, 國內(nèi)僅有蘭州重離子加速器國家實驗室的高能質(zhì)子實驗平臺可提供能量為250 MeV 的質(zhì)子束流, 流強約為107proton/s.由于束流強度偏低, 我們采用TLD[13,16]的測量方法, 通過測量鎢靶中TLD 的劑量分布, 從而獲得能量沉積數(shù)據(jù), 并與模擬計算的結(jié)果進行對比.

質(zhì)子在鎢靶中的能量沉積的驗證實驗示意圖如圖3 所示.整個鎢靶厚度為4 cm, 由厚度分別為1.0, 1.0, 0.5, 0.5, 0.5, 0.5 cm 的幾個靶塊組成.實驗使用的TLD 是由防化研究院生產(chǎn)的GR100系列, 其主要成分是LiF(Mg, Ti).在鎢靶的質(zhì)子入射端面放置了一層TLD, 用于分析質(zhì)子束流強度與分布.在鎢靶靶塊之間布置了若干TLD, 用于分析鎢靶不同深度處的能量沉積分布.為了固定TLD的位置, 加工制作了1 mm 厚的硬鋁支架, 2 個TLD 的間距是3 mm.置于鎢靶內(nèi)部的TLD 大多數(shù)采用“十”字分布.為了降低加速器輸出質(zhì)子脈沖穩(wěn)定性和重復(fù)性影響, 實驗僅接受了1 個質(zhì)子脈沖束團的輻照, 這也能確保所有TLD 都能工作在劑量線性響應(yīng)區(qū)域.

表2 BERTINI, ISABEL, CEM2K 與INCL4 級聯(lián)模型計算1000 MeV 質(zhì)子入射鉛靶產(chǎn)生的不同粒子對總能量沉積值的占比貢獻Table 2.The calculated contribution of different particles to the total energy deposition in lead target by 1000 MeV protons with BERTINI, ISABEL, CEM2K, and INCL4 cascade mode.

圖2 (a)不同模擬程序?qū)︽u靶總能量沉積計算的對比; (b)鎢靶中能量沉積線性密度分布Fig.2.(a) Comparison of total energy deposition in tungsten target calculated by different Monte Carlo code; (b) axial distribution of linear density of energy deposition in tungsten target.

圖3 質(zhì)子在鎢靶中的能量沉積測量示意圖Fig.3.Schematic of the energy deposition measurement in a tungsten target incident by protons.

經(jīng)過輻照的TLD 采用RGD6 型熱釋光儀[17]進行數(shù)據(jù)讀出.RGD6 型熱釋光儀通過測量TLD的積分光產(chǎn)額并與標準60Co 伽瑪源劑量的標定光產(chǎn)額進行對比, 其輸出量的單位是μGy.第一層TLD 的劑量讀出結(jié)果如圖4 所示.可以看出, 在鎢靶幾何中心水平面偏左下方區(qū)域的TLD 劑量讀數(shù)值較大, 這個區(qū)域?qū)?yīng)著加速器輸出質(zhì)子束流的直接照射位置.根據(jù)第一層TLD 測得的質(zhì)子直射區(qū)域,結(jié)合圖3 的實驗布置, 可以看出僅有 (0, –1) 位置的TLD 在質(zhì)子入射深度的各層對應(yīng)位置有測量點.我們采用MCNPX 對單個能量為250 MeV的質(zhì)子直射TLD 時的總能量沉積開展了模擬計算.具有不同線性能量轉(zhuǎn)移值(linear energy transfer, LET)的帶電粒子對TLD 照射相同劑量后, 相對于等同劑量的60Co 伽瑪源射線照射, TLD 的積分發(fā)光量會發(fā)生變化, 這被稱為TLD 的相對發(fā)光效率[18].第一層TLD 被能量為250 MeV 的質(zhì)子直接入射, 對應(yīng)相對發(fā)光效率值[19]約為1.02.TLD測量值的單位是μGy, 計算值的單位是MeV/g, 通過單位轉(zhuǎn)換與發(fā)光效率修正, 對比得到在本次實驗中單個質(zhì)子脈沖輻照鎢靶時, 到達(0, –1)位置TLD 有效面積上的質(zhì)子數(shù)目約為4.74 × 106個.

圖4 第一層TLD 的劑量讀出值Fig.4.TLD dose readouts at the first layer.

3.2 實驗結(jié)果分析

為了獲得不同深度位置TLD 的相對發(fā)光效率值, 利用SRIM[20]程序計算了質(zhì)子在穿透不同厚度的鎢靶后的剩余能量, 并轉(zhuǎn)換為等效水LET 值如圖5 所示.根據(jù)文獻[19], 當入射粒子的等效水LET值在0.4—1.2 keV/μm 范圍時, LiF(Mg, Ti)相對發(fā)光效率約為1.02.

由于質(zhì)子是重金屬鎢靶中能量沉積的主要貢獻者, 結(jié)合TLD 有效面積接受的入射質(zhì)子數(shù)目與相對發(fā)光效率, 可以將TLD 劑量測量值μGy 進行單位轉(zhuǎn)換, 得到單個質(zhì)子在TLD 中引起的能量沉積值, 即MeV/(g·proton).

我們選用MCNPX-CEM2K 級聯(lián)模型對 (0, –1)位置不同質(zhì)子入射深度的TLD 進行了模擬計算.首先, 根據(jù)圖3 中測量實驗布置的幾何條件建立計算模型; 然后, 利用程序中的能量沉積記錄卡對分布在鎢靶內(nèi)部不同質(zhì)子入射深度的TLD 中總能量沉積進行了記錄.TLD 能量沉積模擬計算值與實驗測量值, 以及其二者的比較如圖6 所示.

圖5 不同深度的鎢靶中質(zhì)子的平均能量與水等效LET 值Fig.5.Average proton energy in the tungsten target and equivalent LET in water.

圖6 鎢靶中TLD 的能量沉積測量值與計算值Fig.6.Energy deposition comparison between measurement and calculation of TLD in tungsten target.

可以看出, 總的說來, MCNPX-CEM2K 模型能夠較好地模擬鎢靶中熱量沉積的結(jié)果.通過細致地對比可以發(fā)現(xiàn), 在鎢靶較淺位置, 實驗測量值比模擬計算值略小; 隨著鎢靶深度的增加, 測量值與模擬值的偏差逐步增大; 在質(zhì)子射程末端位置, 二者相差最大.出現(xiàn)較大差別的原因有以下幾點: 首先, TLD 的有效面積較小, 入射質(zhì)子束細微的非垂直入射會對較深鎢靶位置的測量值帶來較大影響;第二, 隨著鎢靶深度的增加, 較高LET 的次級帶電粒子對總能量沉積的占比貢獻增大, 質(zhì)子對總能量沉積的貢獻份額減小, 這就需要對TLD 的相對發(fā)光效率進行更細致的修正.實驗測量數(shù)據(jù)的誤差為10%, 主要來源于質(zhì)子束斑的不規(guī)則、同批次TLD 出廠性能的差異、劑量讀出設(shè)備的系統(tǒng)誤差等.

4 總 結(jié)

本文分別采用MCNPX, PHITS 和FLUKA三種蒙特卡羅粒子輸運程序計算了高能質(zhì)子在金屬散裂靶中的能量沉積分布, 通過與已有實驗測量結(jié)果的比較, 表明這三種蒙特卡羅程序能夠?qū)︺U靶的能量沉積進行較準確的模擬計算, 相對于實驗測量值, 模擬計算的結(jié)果高約10%—15%, 即蒙特卡羅程序的模擬計算結(jié)果能夠比較保守地估計鉛靶中的熱量沉積.對3 個不同模擬程序及同一模擬程序的不同級聯(lián)模型的計算結(jié)果進行比較, 表明MCNPX 程序中的CEM2K 級聯(lián)模型的計算結(jié)果最接近實驗測量值.由于迄今國內(nèi)外尚無高能質(zhì)子在鎢靶中的能量沉積測量數(shù)據(jù), 我們利用熱釋光探測器LiF(Mg, Ti), 測量了能量為250 MeV 的質(zhì)子束在鎢靶中的能量沉積分布, 結(jié)果表明除了在質(zhì)子射程末端有較大的差異外, 能量沉積的實驗測量結(jié)果與MCNPX 的計算值有較好的吻合.通過與實驗數(shù)據(jù)的比較, 驗證了蒙特卡羅模擬程序?qū)δ芰砍练e計算具有較高的可靠性.