醫院中子照射器I型堆熱中子束流孔道等效平面源的模擬計算

朱養妮，江新標，趙柱民，張 良，周永茂

(1.西北核技術研究所，西安 710024;2.中國核工業集團中原對外工程公司，北京 100191)

1 前言

硼中子俘獲治療(boron neutron capture therapy，BNCT)是一種針對腦部神經膠質瘤和黑色素瘤等惡性腫瘤較為理想的治療方法。該方法的優點在于，腫瘤細胞中引入含有10B的化合物后[1]，利用低能中子束照射頭部，形成低能中子對富集10B的腫瘤細胞的靶向照射，通過10B(n，α)7Li反應，放出α粒子、7Li以及足以殺死腫瘤細胞的2.79 MeV能量。從而達到最大限度地殺死腫瘤細胞、保護正常組織的目的。

醫院中子照射器I型堆(IHNI－1)B堆芯設計有BNCT熱中子束流孔道[2，3]，治療時利用從熱中子束流孔道引出的中子束照射于頭部病變部位，幾何示意圖見圖1。熱中子束流孔道的模擬計算屬于粒子深穿透問題，需要特殊的抽樣技巧，如果從堆芯開始跟蹤粒子、模擬計算人體頭顱等效模型內的劑量分布將耗費大量機時，且計算結果方差較大;而從BNCT熱中子束流孔道出口處跟蹤粒子、模擬計算人體頭顱等效模型內的劑量分布只需較少機時，其收斂速度較快、計算結果方差較小，因此，利用MCNP/4B程序在BNCT束流孔道出口處建立等效平面源的方案，為人體頭顱等效模型劑量分布的快速計算提供了一個較為可靠的中子、γ平面源。

圖1 IHNI－1堆B堆芯熱中子束流孔道幾何結構示意圖Fig.1 Transverse sectional chart of the thermal neutron duct of IHNI-1 reactor with B-core

2 IHNI－1堆熱中子束流孔道模擬計算

熱中子束流孔道的數值計算采用蒙特卡羅耦合抽樣技巧和分裂與輪盤賭技巧(或權窗游戲)，計算模式為中子和γ耦合輸運，中子和γ粒子抽樣重要性(IMP:N，P)在堆芯內置1，沿束流孔道方向依次增大。堆芯核功率為30 kW，堆芯單位時間內產生的中子數為每秒2.27948×1015。

堆芯臨界計算是中子束流孔道參數計算的前提，全堆的有效倍增系數keff不同，孔道出口的束流參數也將不同，因此，需要進行臨界搜索計算，得到keff=1的堆芯裝載模式，在此基礎上進行的束流孔道模擬計算結果才可靠。

2.1 堆芯臨界搜索計算

B堆芯由中心控制棒柵元、燃料柵元、水柵元組成，共計11圈，340個柵位。冷態臨界搜索計算時，根據臨界計算結果，調整最外圈燃料柵元和水柵元的數量。

IHNI－1堆B堆芯冷態臨界搜索計算采用了3種堆芯布置，計算結果見表1。由表1可知，堆芯裝載 319.16根燃料元件、235U裝載量為1176.6269 g時，堆芯達到冷態臨界。

表13種堆芯布置下有效倍增系數keff和冷態臨界235U裝料量蒙特卡羅模擬計算Table 1 The Monte Carlo calculation results of keffand load of235U for 3 core schemes

2.2 熱中子束流孔道參數計算

采用蒙特卡羅耦合抽樣法進行了熱中子束流孔道模擬計算，圖2給出了熱中子束流孔道軸線上不同位置處中子通量密度的變化曲線。由圖2可知，堆芯核功率為30 kW時，熱中子通量密度φth在鈹反射層內達到最大值 1.59999 ×1012cm－2·s－1，該點到堆芯中心點的距離為12.5 cm。同時，熱中子束流孔道出口中心處的熱中子通量密度達到2.37795 ×109cm－2·s－1。

3 中子、γ等效平面源計算模型

理論上，置于反應堆BNCT中子源外的人體頭顱等效模型內的通量密度分布可進行一次蒙特卡羅中子—γ耦合輸運計算得到，但該方法計算速度慢、計算精度較差，為了加快計算速度，提高計算精度，將上述一次計算過程劃分為兩步:a.采用蒙特卡羅耦合抽樣方法計算反應堆中子束流孔道出口處垂直于束流孔道軸線的整個平面上的中子、γ平面源;b.采用MCNP/4B程序分別計算中子、γ平面源在人體頭顱等效模型內所產生的中子、γ通量密度。

圖2 熱中子束流孔道中心軸線上中子通量密度衰減曲線Fig.2 Attenuation curve of neutron flux density along axial line of thermal neutron duct

上述中子、γ 平面源 φ(r，E，μ)是空間、能量和方向的函數，其中，r為距離，單位為cm;E為能量，單位為MeV;μ為散射角余弦，數值范圍－1～1。由于MCNP/4B程序在平面源問題的輸運計算時，只需考慮 μ＞0的源項 φ+(r，E，μ)，為了簡化平面源的計算模型，設正向平面源 φ+(r，E，μ)為

式(1)中，φ+(r)為粒子正向流密度的空間分布;φ(r，E)為與空間有關的粒子正向流密度的能譜分布;φ+(r，μ)為與空間有關的粒子正向流密度的角分布;C為歸一化系數。

為了得到平面源的空間分布、能量分布和角分布，采用 MCNP/4B程序中相應的分段計數卡(FSn)、分段除數卡(SDn)、計數能量卡(En)和計數余弦卡(Cn)對平面源的空間、能量和角度進行網格劃分。

針對束流孔徑為10 cm和14 cm的孔道設計方案，將平面源的空間變量r分別劃分為5個網格，網格外徑ri(i=1，…，5)分別為 5 cm、10 cm、20 cm、30 cm、40 cm 和7 cm、10 cm、20 cm、30 cm、40 cm。

等效中子平面源劃分為37個能群，能群劃分見表2;等效γ平面源劃分為18個能群，能群劃分見表3;在平面源的角度上劃分為17個網格，網格劃分見表4。

表2 等效中子平面源能群劃分Table 2 Energy groups of equivalent neutron surface source

表3 等效γ平面源能群劃分Table 3 Energy groups of equivalent γ surface source

表4 等效平面源角分布網格劃分Table 4 Grids of angular distribution for equivalent surface source

采用蒙特卡羅耦合抽樣方法[4]計算了各空間網格內的多群中子、γ正向流密度的能譜分布φ(ri，Eg)，并計算了穿過平面源表面上各空間網格內的正向流密度的角分布φ+(r，μ)。

對φ+(ri，μj)進行方向歸并，得到各空間網格上的正向流密度φ+(ri)

式(2)即為平面源的空間分布。根據式(2)，可得平面源的總源強

式(3)中，Ai為各空間網格的面積，cm2;S為人體頭顱劑量分布計算時的源強歸一化常數。

由于MCNP/4B程序的輸入卡片不能模擬上述一維空間幾何各網格內的平面源，因此還需將此一維幾何平面源轉換為孔道設計的二維x－y等效幾何平面源，其平面源轉換的等效幾何見圖3，圖3中一維同心圓平面和二維x－y平面對應網格的面積相等。

圖3 平面源等效轉換Fig.3 Equivalent surface source mutation

4 中子、γ等效平面源計算結果

采用等效平面源計算模型對B堆芯熱中子束流孔道出口處的中子、γ等效平面源進行了模擬，計算了中子、γ正向流密度角分布、能譜分布等束流參數。

中子束流孔徑分別采用φ10 cm和φ14 cm。當中子束流孔徑為φ10 cm時，等效中子源的歸一常數為2.28224647×1011n/s，等效γ源的歸一常數為3.334001535×1010γ/s;當中子束流孔徑為φ14 cm時，等效中子源的歸一常數為3.4613195×1011n/s，等效γ源的歸一常數為3.27395314×1010γ/s。

圖4和圖5分別給出了中子束流孔徑為10 cm時，熱中子束流孔道出口處中子流密度的角分布和中子正向流密度的能譜分布。

圖4 熱中子束流孔道出口處群中子流密度的角分布Fig.4 Angular distribution of neutron flow density at exit of the thermal neutron duct

圖5 熱中子束流孔道出口處正向中子流密度的能譜分布Fig.5 Energy spectra distribution of neutron flow density at exit of the thermal neutron duct

5 結語

采用蒙特卡羅程序MCNP對醫院中子照射器I型堆(IHNI－1)B堆芯進行了臨界搜索計算，在此基礎上進行了熱中子束流孔道模擬，計算了束流參數，采用等效平面源計算模型，建立了兩種出口孔徑的等效中子、γ平面源。BNCT熱中子束流孔道出口處平面源的建立，為人體頭顱等效模型吸收劑量分布的快速計算提供了一個較為可靠的中子、γ平面源[5]。

[1] Choi J R，Clement S D，Harling O K，et al.Neutron capture therapy beams at the MIT research reactor[C] //Harling O K.Neutron Beam Design，Development，and Performance for Neutron Capture Therapy.New York:Plenum Press，1990:201 －218.

[2] Jiang Xinbiao，Zhu Yangni，Gao Jijin，et al.The conceptual calculation for the neutron beam device at＜In Hospital Neutron Irradiator＞ Mark 1[C] //Advances in Neutron Capture Therapy 2006:Proceedings of 12th International Congress on Neutron Capture Therapy，Takamatsu，Kagawa，Japan，2006.

[3] 江新標，朱養妮.醫院中子照射器I型堆(IHNI－1)BNCT中子束的校核計算報告[R] .NINT09－IHNI1－007.西安:西北核技術研究所，2010.

[4] 江新標，陳 達，謝仲生，等.反應堆孔道屏蔽計算的蒙卡方法[J] .計算物理，2001，18(3):285 －289.

[5] 江新標，張新軍，朱養妮，等.醫院中子照射器 I型堆(IHNI－1)人體頭顱等效模型劑量分布計算[R] .NINT09－IHNI1－002.西安:西北核技術研究所，2003.