NECP-Atlas中热中子散射律数据生成模块的开发与验证

汤勇强，祖铁军，易思宇，徐嘉隆，曹良志，吴宏春

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

热能区中子散射数据对反应堆设计、辐射屏蔽计算、临界安全分析以及冷中子源设计等的计算结果具有重要影响。在热能区，慢化剂靶核的化学键作用和热运动都会强烈影响中子的散射行为，由于热能区中子的能量与靶核相当，中子在与靶核的作用过程中可能从靶核获得能量或失去能量。与自由靶核原子相比，这些效应将影响热能区中子的散射截面，以及次级中子能量和角度分布。评价核数据库(ENDF)采用热散射律(TSL)数据描述以上效应，由于国际上缺少热散射实验，且热散射实验无法对每种慢化材料都覆盖整个热能区的入射中子能量范围和慢化材料的温度范围，因此，ENDF中的TSL数据需结合已有实验，通过物理模型进行评价。ENDF中含有TSL数据的材料极为有限，最新发布的ENDF/B-Ⅷ.0中，仅24类常用慢化材料以及实验冷中子源中的冷慢化剂材料含有TSL数据[1]，这些材料的TSL数据主要是采用美国洛斯阿拉莫斯国家实验室核数据处理程序NJOY的LEAPR模块计算获得的。

目前ENDF中的TSL数据主要存在如下问题：1) 先进反应堆中使用的慢化剂材料缺少TSL数据，如熔盐堆中的FLiBe；2) 材料的TSL数据与其晶体结构具有很强的相关性，如先进反应堆内的ZrH2材料，其中H核素含量不同使得其晶体结构不同[2]，数值结果显示，如果直接使用评价核数据库中给出的ZrH2的TSL数据，会对计算结果造成一定的误差。针对以上问题，需根据材料实际的晶体结构计算产生TSL数据。目前国际上的研究多是采用NJOY程序的LEAPR模块计算产生所需材料的TSL数据，但ENDF/B-Ⅷ.0中TSL的主要贡献者之一北卡罗来纳大学Hawari教授研究指出，LEAPR中使用的近似会引入较大的误差[3]，限制了其使用范围。

近年来，随着反应堆高保真计算方法的发展，核数据的精度成为限制反应堆设计精度提高的主要问题，其中TSL数据也是提高核数据精度的重要内容，因此，需要对TSL数据进行深入研究。Hawari教授针对以上提到的LEAPR问题开发了FLASSH程序[3]。然而，FLASSH未提及对非固体靶核TSL的处理能力。

为满足反应堆设计对高精度核数据的需求，西安交通大学NECP实验室自主研发了核数据处理软件NECP-Atlas[4]，本文在该程序的基础上，开发通用的TSL数据计算模块，去除传统方法中的晶体立方近似和原子位置近似，采用基于各向异性位移参数(ADPs)[5]的TSL数据计算方法，可得到任意固体晶体结构的TSL数据；在固体散射靶的模型基础上，加入扩散模型、离散谐振子模型以及舍尔德近似，使其可产生多原子分子固体材料、液体慢化材料的TSL数据。

1 理论方法

1.1 热中子散射公式

热能区中子在多原子的散射系统中，由于有多个原子参与散射，利用量子力学处理势散射的问题时，需要用多体薛定谔方程来描述该多体系统。因此，一般引入第一玻恩近似和费米赝势求解多体薛定谔方程，热能区中子散射的双微分截面[6]可表示为：

(σcohS(α,β,T)+σincSs(α,β,T))

(1)

式中：E为入射中子能量；E′为出射中子能量；μ为实验室坐标系下的散射角余弦；T为材料温度；kB为玻尔兹曼常数；σcoh为束缚核的相干散射截面；σinc为束缚核的非相干散射截面；α为中子的无量纲动量转移量；β为中子的无量纲能量转移量；S(α,β,T)为热散射函数；Ss(α,β,T)为自散射函数。α和β分别表征中子散射前后角度和能量的变化量，可分别表示为：

(2)

(3)

热散射函数S(α,β,T)可分成两部分：

S(α,β,T)=Ss(α,β,T)+Sd(α,β,T)

(4)

其中：自散射函数Ss(α,β,T)计入非相干散射效应；相互散射函数Sd(α,β,T)加上自散射函数，即S(α,β,T)计入相干散射效应。

引入简谐近似，即原子偏离平衡位置的位移较小的情况下，晶体中原子间的相互作用可近似正比于原子的位移。可通过对散射函数的展开将热中子散射律分解为弹性和非弹两部分。对于晶体，这种展开方式称为声子展开。散射函数用声子展开为：

(5)

(6)

其中，上标0表示0声子产生(或湮灭)，代表弹性散射，中子与靶核散射前后能量没有改变；带有其他更大数值上标的项表示1，2，3等声子项，表示有声子产生(或湮灭)，中子散射后获得能量或失去能量，代表非弹性散射项。热能区中子散射可分为4部分：相干弹性散射、相干非弹性散射、非相干弹性散射和非相干非弹性散射。

通用形式的双微分截面可进一步表示为：

(7)

1.2 固体热中子散射律数据计算模型

固体材料是由大量粒子(原子、离子等)所组成的具有强相互作用的集体。在一定温度下，晶体格点粒子在平衡位置附近小振幅振动，称为晶格振动。晶格振动是由不同频率的简正振动叠加而成的一种集体振动。对于给定的晶体，总的振动模式数目是一定的，按振动频率有一分布，用晶格振动模式密度(或频率分布函数)描述。从量子力学的观点看，表征原子集体振动的能量是量子化的，每个振动模式能量的最小单位称为声子，因此晶格振动模式密度又称声子态密度。根据声子态密度，便可对固体慢化材料的TSL数据进行研究[7]。

1) 固体非弹性散射

中子与所有慢化材料的散射过程中都存在非弹性散射，包括相干部分和非相干部分。一般来说，相互散射律部分较弱，引入非相干近似，即忽略相干散射项中的相互散射律部分，于是，声子展开的非弹性散射可表示为：

σine(E→E′,μ,T)=

(8)

(9)

其中，ρ(β)为声子态密度，是计算TSL数据的重要输入参数，可通过实验方式或第一性原理计算获得。

2) 固体非相干弹性散射

(10)

其中，w(T)为德拜-沃勒积分，表示评价数据库中描述非相干弹性散射的TSL数据。

(11)

3) 固体相干弹性散射

(12)

计算相干弹性散射截面的关键问题是得到S0的计算表达式，其介绍详见文献[6]。从第一性原理出发，推导得到了便于计算的相干弹性散射截面表达式：

(13)

式中：N为晶胞中的原子数；v0为晶胞的体积；κ为倒易晶格矢量；τ为中子散射过程中的动量转移矢量；dj和σcoh,j分别为晶胞中第j个原子的位置和束缚态相干散射截面；e-2Wj(τ,T)为第j个原子的德拜-沃勒因子；Wj(τ,T)为德拜-沃勒系数。

传统方法中，为简化德拜-沃勒因子的计算，采用了2个近似。首先是晶体立方近似，即假设在散射系统中每个原子的位移都是各向同性的，对于各原子可通过使用统一声子态密度来计算德拜-沃勒系数。于是，德拜-沃勒系数可用统一的表达式表示：

(14)

其中，Mj为晶胞中第j个原子的质量。接下来引入第2个近似，即原子位置近似。该近似假设晶胞中原子的位移与原子所处位置和原子类型无关。因此，德拜-沃勒系数可进一步简化为W(τ,T)。采用传统方法的相干弹性散射截面的最终形式为：

(15)

然而，对于非立方晶体，原子在不同位置不同方向所受的作用力是不相同的，这些力是相关的，原子在一个方向上的运动会导致在另一个方向上的力发生变化。原子的各向异性热运动可用ADPs均方位移矩阵U(T)表示，U(T)是一3×3的实对称矩阵，表示原子在平衡位置附近位移的均方振幅，具有长度平方量纲。严格的通用固体相干弹性TSL数据可由式(16)获得，称为各向异性位移参数法[8-9]。

(16)

2π2τTUj(T)τ

(17)

传统方法中，德拜-沃勒系数Wj(τ,T)仅与倒格子向量τ的长度和材料温度有关。采用ADPs方法的Wj(τ,T)不仅与倒格子向量τ的长度有关，还与其方向有关。至此，严格的相干弹性散射截面计算方法推导完成。由于ADPs方法是基于不同原子均方位移计算的，因此，该方法适用于任何固体材料的任何晶体结构。

1.3 液体或气体的散射律计算模型

(18)

1) 有效宽度模型

有效宽度模型一般用于描述主散射体的扩散项，如H2O中的H，其散射律函数解析表达式为：

(19)

其中：K1为第2类修正贝塞尔函数；ωt为平动权重；c为无量纲扩散常数。

2) 自由气体模型

自由气体模型一般用于描述次散射体的扩散项，如H2O中的O，其表达式为：

(20)

1.4 多原子分子内振动模式

无论是固体还是液体散射靶，若是由多原子分子构成，此时分子内部的振动不可忽略，需引入离散谐振子来描述多原子分子的内部振动模式。如甲烷中每个甲烷分子包含4个H原子和1个C原子，其分子内部具有4个特征能量简正振动，分别为162、190、361、374 meV。

(21)

其中：ωi和βi为第i个谐振子的权重和特征能量参数；In(x)为第1类修正贝塞尔函数。通过与式(18)相似的卷积公式，即可将离散振动考虑到热中子散射律中。

1.5 液体分子间的相干效应

前文提到，对于非弹性散射一般可忽略其中相干散射中的相互散射律部分，然而对于束缚核相干散射截面σcoh占总束缚核散射截面σb(σb=σcoh+σinc)比例较大的核素，相干散射贡献部分无法满足非相干近似条件，需引入舍尔德修正因子s(κ)来考虑非弹性散射中的相互散射律部分。如Márquez-Damián提出的重水CAB模型[11]中，由于σD,coh/σD,b=5.592/7.64=0.731 94，在计算D2O中D的非弹性散射时，需考虑相干散射中相互散射律部分的影响。采用舍尔德近似[12]，相干TSLS(α,β,T)可由非相干TSLSs(α,β,T)来计算。

S(α,β,T)=Ss(α/s(κ),β,T)*s(κ)

(22)

(23)

非弹性散射截面为：

(24)

2 程序开发与数值结果

2.1 程序概述

基于以上热中子散射理论，在核数据处理程序NECP-Atlas中开发了sab_calc模块，用于产生慢化材料TSL数据，为反应堆设计提供高精度的TSL数据。sab_calc模块包含NJOY中LEAPR模块的所有功能，且在相干弹性散射中去除传统方法中的晶体立方近似和原子位置近似，采用ADPs方法，得到高精度的相干弹性TSL数据。利用sab_calc模块产生的TSL数据，借助NECP-Atlas程序中其他已有的模块，即可产生中子学计算所需的热散射截面数据。

产生的TSL数据包括描述非弹性散射的Stot(α,β,T)、非相干弹性散射的德拜-沃勒积分w(T)，以及相干弹性散射的Si(E,T)。目前，国际上的ENDF均采用ENDF-6格式[14]存储数据，为使所开发程序输出的数据与ENDF一致，程序根据慢化材料的类型，将计算得到的TSL数据按ENDF-6格式进行储存。该格式按照主散射体的非弹性散射类别、主散射体的弹性散射是否为相干散射、是否有次散射体分成表1所列类型。

2.2 非弹性散射截面的验证

液态重水的非弹性散射极为复杂，包括连续谱项、扩散项、离散振荡频谱和舍尔德近似等的处理，因此，本文以液态重水为例，将NECP-Atlas产生的非弹性散射TSL数据和热散射截面与NJOY2016进行对比，验证程序的正确性。

图1为NECP-Atlas中sab_calc计算的D2O中D在293.6 K下的非弹性TSL数据，每条曲线表示1个α值下Stot(α,β,T)随β的变化曲线。可看到，在α很小时，曲线上出现了几个峰值，证实了水分子的扩散运动和水分子内的振动运动。图2为NECP-Atlas采用D2O中D在293.6 K的TSL数据计算的非弹性散射概率密度及其与NJOY2016的相对偏差，图中每条曲线代表1个入射能量。从图2b可看出，NECP-Atlas与NJOY2016产生的散射概率密度符合良好，在小部分能区相对偏差略偏大，主要在准弹性散射峰附近。准弹性散射指中子散射前后能量的变化远小于中子本身的能量。相对偏差较大是因为在准弹性散射峰附近，散射概率变化剧烈，需要更密集的出射能量网格才能满足线性化收敛准则，在NJOY2016程序内存中对出射能量网格采用8位有效数字，而NJOY2016输出的文件中，有效位数最多仅为7位，输出过程中会产生偏差，由于输出数据有效位数的限制，在程序内部采取更高有效位数实际上是无意义的，NECP-Atlas在判断收敛时，考虑了输出数据有效位数的限制，避免了以上问题。以上不同的收敛策略导致NECP-Atlas与NJOY2016得到的散射概率相对偏差较大。D2O中D在293.6 K、入射能量0.025 3 eV的几个出射能量点及散射概率列于表2。可看到，在NJOY2016内存中有5组有效数据，而在输出文件中有效数据仅为(0.025 299 94，370.374 3)和(0.025 299 95，371.729 0)，证实在输出过程中产生了偏差。为进一步证明，图3示出了两者基于相同能量网格D2O中D在293.6 K、入射能量0.025 3 eV的非弹性散射概率的相对偏差，可见相对偏差均小于10-4，两者结果吻合很好，证明由不同收敛策略产生的不同能量网格会带来偏差。

表1 ENDF-6格式中热散射律数据类型Table 1 Type of thermal scattering law data in ENDF-6 format

图1 D2O中的D在293.6 K下的散射律Fig.1 TSL for D in D2O at 293.6 K

图2 NECP-Atlas计算的D2O中D在293.6 K的非弹性散射概率密度(a)以及与NJOY2016的相对偏差(b)Fig.2 Thermal inelastic scattering probability for D in D2O at 293.6 K calculated by NECP-Atlas (a) and comparison with NJOY2016 (b)

表2 NJOY2016程序内存和输出文件数据对比Table 2 Comparison of values calculated by NJOY2016 in memory and output file

2.3 非相干弹性散射截面的验证

表3为本文计算的LiH中H的德拜-沃勒积分和有效温度与文献值的对比，可看出，本文计算的数据与文献值符合良好。

2.4 相干弹性散射截面的验证

采用传统的晶体立方近似方法和ADPs方法分别计算了金属Be的相干弹性TSL数据(图4a)，可看到，sab_calc采用传统方法与LEAPR结果符合很好，采用ADPs方法会导致TSL数据产生最大10%的偏差。这种偏差会导致由TSL数据计算的相干弹性散射截面产生最大10%的偏差(图4b)。

图3 基于相同能量网格非弹性散射概率的相对偏差Fig.3 Comparison of calculated thermal inelastic scattering probability based on the same secondary energy grids

表3 LiH中H的德拜-沃勒积分和有效温度Table 3 Debye-Waller integral and effective temperature for H in LiH

图4 296 K下金属Be的相干弹性TSL数据(a)和相干弹性散射截面(b)Fig.4 Coherent elastic TSL data (a) and coherent elastic scattering data (b) for metal-Be

为检验以上相干弹性散射数据的偏差对反应堆中子学计算结果造成的影响，使用NECP-Atlas程序基于ENDF/B-Ⅷ.0中子库和计算的TSL数据制作蒙特卡罗程序使用的ACE格式数据库，并采用MCNP程序计算了27个以金属Be为慢化剂的国际临界基准题(ICSBEP)[17]HEU-MET-THERM-026(HMT026)。图5示出了对金属Be分别基于LEAPR、晶体立方近似方法的sab_calc和ADPs方法的sab_calc制作的热散射库，计算得到的基准题keff，所有的统计偏差为24 pcm。结果表明，晶体立方近似方法的sab_calc结果和基于相同假设的LEAPR结果吻合，证明了sab_calc程序开发正确；其次，相对于LEAPR，sab_calc的ADPs方法对27个基准题均有改进，与实验基准结果的偏差平均降低60 pcm，最大偏差为HMT026_19基准题的135 pcm。

图5 金属Be采用ADPs方法的临界基准题keffFig.5 Calculated keff of critical benchmarks with metal beryllium using ADPs method

3 结论

基于热中子散射理论，本文在核数据处理程序NECP-Atlas中开发了TSL数据产生模块sab_calc，该模块可产生任意晶体结构的固体材料、部分液体材料的TSL数据，通过数值结果证明了程序的正确性。该程序采用ADPs方法，提高了相干弹性TSL数据的精度，基于ICSBEP基准题的计算结果表明，采用ADPs方法获得的金属Be的TSL数据，与传统晶体立方近似相比，会使keff的计算结果与实验基准值的偏差平均降低约60 pcm。