核数据处理软件NECP-Atlas中的光子相关数据计算方法研究

祖铁军，徐 宁，尹 文，曹良志，吴宏春

(西安交通大学 核科学与技术学院，陕西 西安 710049)

核数据处理软件是联系评价核数据与反应堆核设计的纽带，可为反应堆核设计提供基础核数据库(应用核数据库)，是关键的反应堆核设计软件之一，其精度将直接影响最终核设计结果的精度。长期以来，国内外主要使用美国洛斯阿拉莫斯国家实验室的NJOY[1]、国际原子能机构的PREPRO[2]等程序制作应用核数据库，由此会产生两方面问题：一是这些软件主要针对中子输运计算的需求开发，能提供的数据类型有限，不能满足反应堆核设计各方面的需求；二是核数据处理理论模型、数据库制作方法研究基本停滞，并未随核设计方法的发展得到相应的发展，成为核反应堆设计精度提高的障碍。近年来，随着核反应堆高保真数值模拟技术的发展，国际上对应用核数据库制作方法开展了一定研究[3]，这些研究仅是在NJOY程序基础上对应用数据进行的针对性的改进。2015年国际原子能机构成立了核数据处理项目[4]，旨在推动核数据处理方法及软件的发展，国际上一批新的核数据处理软件参与了该项目，如中国核数据中心的Ruler、西安交通大学的NECP-Atlas、日本原子能机构的FRENDY、法国原子能委员会的GALIEE等。

NECP-Atlas是西安交通大学研发的核数据处理软件[5]，旨在建立核数据处理研究平台，开展高精度的核数据处理方法研究，实现核数据处理软件的自主化。该软件自2015年开始研发，目前已具备丰富的核数据处理能力，可为反应堆核设计提供不同类型的核数据，包括中子输运计算的反应截面数据库[6]、燃耗数据库[7-8]、屏蔽数据库[9]、反应截面的多群协方差数据[10]、裂变产额的协方差数据[7-8]、释热与辐照损伤数据[11]、热散射律数据[12]、活化及源项数据库[13]等；开展了先进的核数据处理方法研究，如共振弹性散射效应处理方法[14-15]、热能散射数据计算方法[16]等，提高了应用数据库的计算精度；基于NECP-Atlas软件研制了国产化的应用核数据库，将我国最新的评价核数据库CENDL-3.2应用于核反应堆物理设计及屏蔽设计软件[17-18]。

光子是核反应堆内重要的辐射粒子，是堆内释热率、材料辐照损伤等的重要来源。目前国际上现有的核数据处理软件仅可给出中子的非弹性散射、辐射俘获、裂变等核反应释放的瞬发光子，不能给出裂变产物衰变过程中释放的光子，这将影响反应堆内功率或释热率的计算精度[19]。对于辐照损伤截面数据，目前核数据处理软件仅考虑中子入射产生的辐照损伤，不能提供光子引起的辐照损伤截面数据，而Remec等[20]和Alexander等[21]的研究发现，光子的离位损伤是高通量核素生产堆、沸水堆压力容器材料辐照损伤的重要组成部分。针对以上问题，本文拟进行光子产生数据、光子辐照损伤数据计算方法研究，以完善NECP-Atlas的功能。

1 多群缓发光子产生矩阵计算功能

光子为中性粒子，与中子类似可在堆芯内输运一定距离，在输运过程中将其能量沉积在发生核反应的位置。核数据处理软件需要为光子输运计算提供光子产生、光子与原子反应相关的截面数据。堆芯内光子按产生方式可分为2类：第1类是在中子引发的非弹性散射、辐射俘获、裂变等反应发生的瞬间释放的光子，称为瞬发光子；第2类是裂变产物衰变过程中释放的光子，称为缓发光子。传统的核数据处理程序可根据中子反应核数据计算中子引发核反应产生瞬发光子的截面，但未考虑缓发光子的产生数据。对于燃料组件，光子对总释热率的贡献约占10%[19]，其中缓发光子贡献约为30%[22]，所以一般忽略缓发光子对释热的贡献。但钠冷快堆等反应堆中包含了钢反射层、控制组件等不含裂变材料的组件，这些组件内的释热率主要来自于光子释热，文献[19]中指出，快堆反射层内光子对释热率的贡献达90%，因此，需要精确的光子产生数据。NECP-Atlas根据评价数据库提供的裂变产额、衰变子库提供的数据，计算产生多群缓发光子产生矩阵，其计算理论如下。

多群缓发光子产生矩阵的计算公式可表示为：

(1)

其中：yi(h→g)为裂变核素i的第h群中子对第g群缓发光子的产额；Eh和Eh-1为中子能群边界；p为裂变产物编号；yi→p为裂变核素i产生裂变产物p的产额，由裂变产额子库给出；Ein为入射中子能量；Dg为多群衰变光子产额，根据衰变子库中核素p衰变光子能谱的形式分别采用式(2)、(3)进行计算，前者用于光子能谱为离散能级的情况，后者用于连续能谱情况。

(2)

(3)

(4)

(5)

其中，〈Ek〉为第k种衰变能谱的平均衰变光子能量。以JEFF-2.2评价核数据库为例，部分核素的Ak值列于表1。由表1可看出，Ak值基本在1.0附近，表明了上述对衰变光子数据处理的准确性。

表1 部分核素的Ak值Table1 Ak values of some nuclides

在评价核数据库的458反应道(MT=458)给出了核素裂变后释放的衰变光子能量，为保证以上通过采用衰变数据计算的光子能量与458反应道给出的总能量一致，将式(1)计算的多群缓发光子产生矩阵做进一步处理：

Yi(h→g)=f(g)·yi(h→g)

(6)

(7)

其中，Qd为评价核数据库中458反应道给出的平均裂变缓发光子能量。对于未给出MT=458反应道的评价核数据库，缓发光子产生矩阵未进行以上修正。

多群缓发光子产生矩阵计算流程如图1所示。首先利用评价核数据库衰变子库给出的衰变光子能谱，基于式(2)、(3)计算多群衰变光子产额Dg；然后利用裂变产额子库给出的独立裂变产额，基于式(1)计算多群缓发光子产生矩阵yi(h→g)；最后为保证缓发光子释热计算的准确性，利用中子子库中给出的平均裂变缓发光子能量Qd对计算所得多群缓发光子产生矩阵进行修正，得到最终的多群缓发光子产生矩阵Yi(h→g)。若评价核数据库中未给出平均裂变缓发光子能量，则不对yi(h→g)进行修正。

图1 多群缓发光子产生矩阵计算流程Fig.1 Workflow to generate delayed photon library

2 光子离位辐照损伤截面计算功能

光子在输运过程中与堆内材料原子发生反应，包括康普顿散射、电子对效应、光电效应等。这些反应将产生次级电子，电子与原子核发生碰撞后会造成原子离位损伤。因此，计算光子的离位损伤截面首先需建立电子离位辐照损伤截面计算模型，然后在此基础上建立光子离位辐照损伤截面计算模型。需要指出的是，此处建立的电子离位辐照损伤截面计算方法，可直接单独用于电子入射时产生的离位辐照损伤截面的计算。

2.1 电子离位损伤截面计算方法

电子与原子核发生碰撞的离位辐照损伤截面的计算公式为：

(8)

其中：υ(T)为材料的离位损伤函数，可采用K-P模型、NRT模型或ARC模型[11]计算获得；T为反冲核的能量；Tmax为电子能产生的反冲核能量的最大值。

(9)

其中：E为电子能量；M为材料原子核质量；m0为电子静止质量；c为光速。

式(8)中，dσ/dT为电子的微分散射截面，可表示为Mott截面形式：

(10)

σMott=σRR(θ,Z,E)

(11)

其中，σR为卢瑟福微分散射截面，可表示为：

(12)

对于R(θ,Z,E)，文献[23]针对1 keV～900 MeV的电子与1≤Z≤118的元素反应，提出一种插值多项式形式，如式(13)所示：

(13)

(14)

(15)

参照文献[23]针对正负电子分别拟合获得系数bk,j(Z)。最终，通过以上方法可计算不同入射电子能量E、不同散射角度θ下的电子微分散射截面。

2.2 光子离位损伤截面计算方法

光子辐照损伤截面可通过对3种光子-原子反应道产生的离位损伤截面得到：

σt(Eγ)=σCS(Eγ)+σPE(Eγ)+σPP(Eγ)

(16)

其中：σCS(Eγ)、σPE(Eγ)、σPP(Eγ)分别为康普顿散射、光电效应和电子对效应对总光子离位损伤截面的贡献；Eγ为入射光子能量。为获得光子的离位损伤截面，需分别计算这3个分反应道的离位损伤截面。

1) 康普顿散射

对于康普顿散射，其离位损伤截面可表示为：

(17)

(18)

式(17)中的积分上限可表示为：

(19)

式(17)中电子离位损伤函数可通过下式积分得到：

(20)

其中：NV为材料的原子核密度；σD(T)为电子离位损伤截面，通过式(8)计算；S(T)为电子阻止本领，本文采用ICRU1984电子组织本领数据库[25]计算，对于该数据库中未给出正电子组织本领的材料采用文献[26]中的方法近似计算。

2) 光电效应

对于光电效应，电子的动能可通过以下公式计算：

E=Eγ-Be

(21)

其中，Be为电子结合能，一般是几百eV的量级，相比于光子能量这个数值很小，可忽略不计。因此，对于一给定的入射光子能量，出射电子能量是一固定值，光电效应的光子辐照损伤截面可表示为：

σPE(Eγ)=σPE(E)n(E)

(22)

其中：n(E)为电子离位损伤函数，与式(8)中υ(T)的计算方法相同；σPE(E)为光电效应反应截面，可用Hall公式[27]表示：

(23)

其中：α=Z/137；γ为Lorentzian系数，可表示为：

(24)

3) 电子对效应

对于电子对效应，离位损伤截面可通过下式计算：

(25)

(26)

其中：σ∞=5.18×10-28cm2；变量F(s)通过下式计算：

(27)

其中：l=2；n=8。其他变量的计算公式如下：

s=E/(Eγ-2m0c2)

(28)

u=ln(Eγ/m0c2)

(29)

g(u)=-0.183 5u3+1.653u2-

2.154 3u+0.761 4

(30)

h(u)=0.219 3u+0.182 5

(31)

3 数值验证

3.1 多群缓发光子产生矩阵验证

本文针对美国阿贡国家实验室的EBR-Ⅱ钠冷快堆，采用先进反应堆计算程序SARAX基于欧洲快堆计算程序ERANOS的数据库进行光子释热率计算。EBR-Ⅱ堆芯采用六边形栅元布置，由16圈共637个组件构成，其中组件的对边距为5.892 9 cm。根据组件的分布，可将整个堆芯分成3部分，由内到外分别为燃料区、反射层区和外增殖区。在燃料区内除燃料组件外，还布置有8个控制棒组件和2个安全棒组件。ERANOS数据库主要基于JEFF-2.2[28]研制，包含多群瞬发与缓发光子产生矩阵。本文采用相同的评价核数据库制作多群瞬发和缓发光子数据库，并与基于ERANOS数据库的计算结果进行对比。采用本文方法获得的考虑缓发光子产生矩阵的数据库计算的堆内不同组件的光子功率及ERANOS的数据库的结果及比较示于图2。由图2可见，两者的最大相对偏差为0.93%，表明了本文方法的正确性。为体现缓发光子对光子功率计算结果的影响，计算了不考虑缓发光子和考虑缓发光子计算得到的光子功率的相对偏差，结果示于图3，其中，燃料组件的最大偏差为32.91%，控制棒组件和反射层组件最大相对偏差为32.58%、20.41%。

图2 组件光子功率计算结果Fig.2 Calculated result of photon heat for assembly

图3 考虑与不考虑缓发光子对组件光子功率计算结果的影响Fig.3 Photon heat for assemby with and ignoring delayed photon contribution

为考虑缓发光子对商用压水堆光子通量的影响，利用压水堆燃料管理程序NECP-Bamboo[29]对VERA-2D基准题的光子通量进行计算，中子能群结构选用WIMS-69群能群结构，光子能群结构选用美国洛斯阿拉莫斯国家实验室给定的48群能群结构[1]。不考虑缓发光子和考虑缓发光子计算得到的组件光子总通量示于图4。从图4可见，两者最大相对偏差为11.68%。

图4 VERA-2D组件光子总通量计算结果Fig.4 Photon flux calculation result for VERA-2D

3.2 电子离位辐照损伤截面的验证

电子离位辐照损伤截面的计算是光子离位辐照损伤截面计算的基础，因此，首先对其计算精度进行验证。选取美国橡树岭国家实验室Oen等[30]的正、负电子离位损伤截面作为参考解，该实验室采用K-P模型[5]用于材料原子离位损伤函数计算，NECP-Atlas计算中采用相同的模型。负电子、正电子引起的Fe、Au的离位辐照损伤截面示于图5。由图5可见，NECP-Atlas获得的结果与Oen的结果符合良好。

3.3 光子离位辐照损伤截面验证

选取Fukuya等[26]的结果作为参考解验证NECP-Atlas计算的光子离位损伤截面的准确性，结果示于图6。由图6可看出，NECP-Atlas计算得到的光子离位损伤截面与参考解吻合良好，证明了程序的准确性。

a、b——Fe和Au的负电子离位损伤；c、d——Fe和Au的正电子离位损伤图5 Fe和Au的负电子、正电子离位损伤截面Fig.5 Displacement cross section of Fe and Au caused by electron and positron

图6 Fe的光子离位损伤截面Fig.6 Gamma displacement damage cross section of Fe

4 结论

本文基于核数据处理程序NECP-Atlas开展了光子相关数据计算方法研究，在程序中建立了完善的光子数据计算方法，除传统核数据处理程序可产生的中子核反应释放的瞬发光子产生截面、光子与原子的反应截面，针对NECP-Atlas程序新开发了裂变产物衰变释放的缓发光子多群产生矩阵、光子离位辐照损伤截面计算功能。数值验证结果表明，钠冷快堆中缓发光子对非燃料组件的释热将产生显著影响；针对光子离位辐照损伤截面的数值验证证明了程序的正确性。此外，由于光子离位辐照损伤截面的计算需以电子离位辐照损伤截面的计算为基础，因此，形成了电子离位辐照损伤截面的计算功能，并对该功能进行了验证。