医院中子照射器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.