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

朱養(yǎng)妮，江新標(biāo)，趙柱民，張 良，周永茂

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

1 前言

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

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

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

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

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

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

2.1 堆芯臨界搜索計(jì)算

B堆芯由中心控制棒柵元、燃料柵元、水柵元組成，共計(jì)11圈，340個柵位。冷態(tài)臨界搜索計(jì)算時，根據(jù)臨界計(jì)算結(jié)果，調(diào)整最外圈燃料柵元和水柵元的數(shù)量。

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

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

2.2 熱中子束流孔道參數(shù)計(jì)算

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

3 中子、γ等效平面源計(jì)算模型

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

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

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

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

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

針對束流孔徑為10 cm和14 cm的孔道設(shè)計(jì)方案，將平面源的空間變量r分別劃分為5個網(wǎng)格，網(wǎng)格外徑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個網(wǎng)格，網(wǎng)格劃分見表4。

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

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

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

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

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

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

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

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

圖3 平面源等效轉(zhuǎn)換Fig.3 Equivalent surface source mutation

4 中子、γ等效平面源計(jì)算結(jié)果

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

中子束流孔徑分別采用φ10 cm和φ14 cm。當(dāng)中子束流孔徑為φ10 cm時，等效中子源的歸一常數(shù)為2.28224647×1011n/s，等效γ源的歸一常數(shù)為3.334001535×1010γ/s;當(dāng)中子束流孔徑為φ14 cm時，等效中子源的歸一常數(shù)為3.4613195×1011n/s，等效γ源的歸一常數(shù)為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 結(jié)語

采用蒙特卡羅程序MCNP對醫(yī)院中子照射器I型堆(IHNI－1)B堆芯進(jìn)行了臨界搜索計(jì)算，在此基礎(chǔ)上進(jìn)行了熱中子束流孔道模擬，計(jì)算了束流參數(shù)，采用等效平面源計(jì)算模型，建立了兩種出口孔徑的等效中子、γ平面源。BNCT熱中子束流孔道出口處平面源的建立，為人體頭顱等效模型吸收劑量分布的快速計(jì)算提供了一個較為可靠的中子、γ平面源[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] 江新標(biāo)，朱養(yǎng)妮.醫(yī)院中子照射器I型堆(IHNI－1)BNCT中子束的校核計(jì)算報(bào)告[R] .NINT09－IHNI1－007.西安:西北核技術(shù)研究所，2010.

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

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