高能質子在散裂靶中的能量沉積計算與實驗驗證*

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

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

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

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

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

1 引 言

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

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

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

2 能量沉積模擬計算

2.1 計算方法

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

2.2 質子入射鉛靶的能量沉積計算結果

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

圖1 (a)不同模擬程序對鉛靶總能量沉積計算的對比;(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]中的測量數據進行比較, 可以看出, 高能質子在鉛靶長度方向的能量沉積線性密度計算結果與測量結果在整體上能較好符合; 隨著質子入射深度的增加, 能量沉積線性密度計算值的偏差略有增大.

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

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

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

2.3 質子入射鎢靶的能量沉積計算結果

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

3 質子在鎢靶中能量沉積的實驗驗證

3.1 實 驗

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

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

表2 BERTINI, ISABEL, CEM2K 與INCL4 級聯模型計算1000 MeV 質子入射鉛靶產生的不同粒子對總能量沉積值的占比貢獻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)不同模擬程序對鎢靶總能量沉積計算的對比; (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 質子在鎢靶中的能量沉積測量示意圖Fig.3.Schematic of the energy deposition measurement in a tungsten target incident by protons.

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

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

3.2 實驗結果分析

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

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

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

圖5 不同深度的鎢靶中質子的平均能量與水等效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 模型能夠較好地模擬鎢靶中熱量沉積的結果.通過細致地對比可以發現, 在鎢靶較淺位置, 實驗測量值比模擬計算值略小; 隨著鎢靶深度的增加, 測量值與模擬值的偏差逐步增大; 在質子射程末端位置, 二者相差最大.出現較大差別的原因有以下幾點: 首先, TLD 的有效面積較小, 入射質子束細微的非垂直入射會對較深鎢靶位置的測量值帶來較大影響;第二, 隨著鎢靶深度的增加, 較高LET 的次級帶電粒子對總能量沉積的占比貢獻增大, 質子對總能量沉積的貢獻份額減小, 這就需要對TLD 的相對發光效率進行更細致的修正.實驗測量數據的誤差為10%, 主要來源于質子束斑的不規則、同批次TLD 出廠性能的差異、劑量讀出設備的系統誤差等.

4 總 結

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