氮化铀热中子截面的第一性原理计算

王立鹏 江新标 吴宏春 樊慧庆

1)(西安交通大学核科学与技术学院,西安 710049)

2)(西北核技术研究所,西安 710024)

3)(西北工业大学材料学院,西安 710072)

氮化铀(UN)因其较好的热物性和耐事故容错性成为先进动力堆的候选燃料,但目前热能区缺少可靠的UN热中子截面数据,这对于热中子反应堆物理计算是很不利的.本文基于量子力学的第一性原理,利用VASP/PHONON软件模拟计算了UN的声子态密度,以此为积分得到UN的定容比热容,并基于新制作的声子态密度,采用核截面处理程序NJOY/LEAPR,利用热中子散射理论,得到UN的S(α,β)数据,进而研究UN的热中子散射截面,并与传统压水堆的二氧化铀(UO2)进行对比.结果表明:优化的晶格参数与数据库符合较好,UN声子态密度的声子项和光子项较UO2的分隔更加明显,定容比热容计算结果与实验值一致,基于该声子态密度计算得到的UN中238U的非弹性散射和弹性散射截面比相同温度下UO2中238U小,UN中N仅考虑了非相干散射部分,随着温度升高,UN弹性散射截面变小,非弹性散射变大,并在高能段趋于自由核散射截面.本文的研究结果填补了UN热中子截面数据的缺失,为下一步系统研究UN燃料在轻水堆中的中子学性能奠定了基础.

1 引 言

相较于二氧化铀(UO2)燃料,氮化铀(UN)具有铀密度高、熔点高、热导率高、热膨胀系数低、辐照稳定性好、裂变气体释放率低等优点,成为先进动力堆的候选燃料,也是新型的耐事故容错燃料,具有较好的发展前景[1,2].热中子反应堆设计计算使用的核材料的热中子截面对反应堆临界安全特性、中子能谱等都会产生较大影响,需要较为可靠的热中子截面数据才能准确地计算出堆芯的物理参数.目前国际上已建立的多种常用慢化剂材料的热中子截面产生技术和处理方法大多数基于半经验提出的简化模型,再引入较多近似从而得到热中子截面数据.近年来,国外发展的先进模拟与仿真技术,依托高性能计算能力,以“第一性原理”的物理学为基础建立模型,取代经验公式,已经实现了许多复杂系统的预测性模拟,包括热中子截面的模拟计算[3−5].以往的热中子截面库仅给出了传统的UO2燃料的热中子截面数据,比如MCNP5自带的热中子截面ENDF70SAB数据文件,缺少UN燃料的热中子数据[6],这样就会给以UN为燃料的先进轻水堆或者其他类型热中子反应堆的物理计算带来较大误差.最新版的ENDF/B VIII.0中,美国北卡罗来纳州立大学的LEIP实验室利用“in-house code”制作了UN的TSL(thermal scattering library)文件库[7−9],相干弹性散射部分归并到U核内,不相干弹性散射归并到14N核内,并对NJOY程序[10,11]中LEAPR模块进行了大幅度修改,其中相干弹性散射部分采用了Cubic approximation和Exact Debye-Waller Matrix方法,仅依靠现有TSL库,没有相关程序,无法对其进行更加深入的研究.因此,本文基于UN热中子数据的迫切需求,研究了UN精确声子态密度的产生方法和验证,对NJOY现有模型稍加改动,使其适用于对UN的研究,以补充UN燃料的热中子截面数据.

热能区的中子散射截面不单纯与中子能量有关,还与散射介质的温度及物理、化学性质有关,它反映的是材料自身的声子态密度特征,热中子截面产生技术与处理方法是一项多学科相互交叉的难题.在热中子反应堆中,对于能量低于1 eV的中子,由于中子能量与散射核的热能是可比较的,不能简单认为靶核是静止的.由于热中子与运动的靶核发生碰撞时,能从靶核获得能量,所以,出射中子的能量可能大于入射中子的能量,这即是热能区中子的向上散射现象.在分子或固体中,散射核与邻近核之间存在着相互作用,原子核处于束缚状态,与中子发生碰撞时不能自由地反冲.由于较低能量中子的德布罗意波长可与分子或晶体内核的间距相比较,与不同核发生散射的中子间可能发生干涉效应[12,13].一般用热散射S(α,β)来表示热中子截面,ENDF格式评价库中有对热中子截面专门的TSL文件描述[14],特定的评价数据库仅给出了少数常用慢化剂材料的TSL文件,一般用户无法根据自己的需要处理评价库以外的材料,要制作UN的热中子散射截面数据就需要UN准确的声子态密度和可靠的热中子截面处理方法.

2 热中子散射理论

热中子相干散射和非相干散射包括弹性散射和非弹性散射部分,弹性散射不会带来系统能量的变化,但是由于热中子能量与晶格振动能态相当,认为热中子弹性散射是整个晶格的散射,这样靶核的有效质量将会很大,系统在中子散射过程中不会失去能量.相反,热中子的非弹性散射会带来能量的损失或者增加,这主要是由碰撞核处于激发态引起的,即非弹性散射会伴随原子的一个或多个声子的发射或吸收.虽然热中子的非弹性散射不会引起整个靶核本身的激发态,但可以引起分子(或晶体)中原子的一个或几个振动量子态的改变.

热中子截面通常可以分为3部分[15].

1)非弹性散射:包括相干和非相干非弹性散射,对所有物质都重要,用热散射律表示.非弹性散射的双微分散射截面是散射律S(α,β)的函数:

其中E和E′是实验系中入射和散射中子的能量,Ω是实验系下的散射角度,σcoh是材料的束缚核的相干散射截面,σinc是材料的束缚核的非相干散射截面,α为动量转移量,β为能量转移量,S(α,β)为热中子散射律.表达式分别为

其中A表示散射核的原子质量M和中子质量的比值;κ为散射矢量;kT是温度,单位是eV;~为普朗克常数;ε为能量变化量;SS(α,β)为不考虑相互作用的自散射律,Sd(α,β)为考虑相干效应的分立散射律.由(1)式可以看出,要求非弹性散射截面需要知道束缚核的截面以及相应的散射律.以往在热中子散射计算中通常引入“非相干近似”,比如NJOY程序中的LEAPR模块,即认为(4)式中的Sd(α,β)为零,这样(1)式变为

对于引入非相干近似后的 (5)式,认为SS(α,β)是高斯分布的函数,这时,热中子散射律SS(α,β)可以写成

其中的ρ(β)是以β为单位的声子态密度,并且是归一化的.

当α和|β|值较大时(α > αsw,|β|> βsw), 即在入射能量较高时,α和β已经超出了S(α,β)范围,这时就需要引入短时间碰撞近似(SCT)[10,11,15],从而热中子的非弹性散射表达变为

进而得到“短时间碰撞下”的双微分截面为

2)非相干弹性散射:对含氢的固态物质重要,如ZrHx、聚乙烯、固态轻水等.在含氢固体中,有一种由于零阶声子项产生的弹性散射(不造成能量损失),称为“非相干弹性散射”,对于非相干弹性散射,其双微分散射截面为

µi是第i个等概率区间的散射角余弦的上限,¯µi是该区间内平均散射角余弦值,N是等概率区间的个数,并且µ0=−1.

3)相干弹性散射:对晶体重要,如石墨、铍和UO2等.在包含相干散射的固体中,组成固体的晶体不同平面的原子会发生干涉散射,在ENDF库中,这一过程被称为“相干弹性散射”,因为没有能量的损失.

多晶材料的双微分相干弹性散射表达式如下:

其中σb是材料有效的束缚核相干散射截面,Ei是所谓的布拉格阈值,fi是对应的晶格学结构因子,不同晶格结构的fi不同.

3 计算方法与建模

氮化铀燃料的原子结构如图1所示,它与碱金属卤化物NaCl,KCl和MgO的结构类似,这种结构被称之为面心立方结构(FCC),它的原胞结构中只包含两个原子,每个原子被6个其他原子所包围.UN燃料结合了金属燃料和氧化物燃料的双重优点,既有像金属燃料一样的高热导率和高密度,又有像氧化物燃料一样的高熔点和较高的结构完整性.

图1 UN原子结构图Fig.1.Atom structure of UN.

利用美国MaterialDesign公司研制的MedeA[16]材料计算平台,通过调用平台下的VASP和PHONON软件完成UN材料的声子态密度计算.本文首先计算了UN在基态时的能量与结构,交换关联函数使用Perdew在1991年提出的梯度密度修正近似GGA(Gradient Corrected Approximation)[17],U和N的截断能量选取为400 eV,计算采用周期性超晶格方法,UN立方相的布里渊区积分在6×6×6的Monkhost-Pack格子中进行.根据体系的周期性,移动每个原子的位置,计算出原胞中所有原子受力,利用第一性原理线性响应理论计算Hellmann-Feynman(HF)力常数,进而得到UN的声子谱,流程如图2所示.

图2 VASP/PHONON产生声子态密度的流程图Fig.2.Flowchart of phonon density of states generation in VASP/PHONON code.

热力学函数内能(E)、熵(S)、自由能(F)、恒容热容(Cv)与声子谱密切相关,一旦得到了UN的声子谱,这些热力学参数就可以在简谐近似的模型下确定,因为这些参数主要采用了声子谱作为它们积分的权重谱,其中Cv是比较重要的参数,因为它可以在实验中准确测量,通常用它来检验声子谱计算的准确性,在PHONON/VASP程序的计算中,比热容有如下表示:

其中r是晶胞内的自由度个数,kB是玻尔兹曼常数,T是温度,~是简化普朗克常数.

因为14N在热能区具有较高的热中子吸收截面,而且还会产生半衰期较长的14C,所以在UN的实际应用中,多采用高富集度的U15N.UN中U的热中子分析采用和UO2一致的方法,仅处理UN中U的238U核素,表1所列为UN中238U,14N和15N各核素原子量以及各自束缚核的相干和非相干散射截面[18],其中自由散射截面由下式得到:

从表1中数据可以看出,238U采用相干弹性散射截面和非相干非弹性散射,14N和15N采用非相干弹性散射截面是比较合理的.因此,在NJOY/LEAPR的计算中,仅需要对UN中238U的相干弹性散射截面的晶格学结构因子进行重新构置,其他参数按照表中数据和声子态密度进行输入,而对于非弹性散射截面部分,则采用程序自带的“非相干近似”忽略其中的相干部分.

表1 UN各核素束缚核的不同反应截面Table 1.Dif f erent cross sections of isotopes in UN.

4 计算结果与讨论

4.1 晶格参数的比较

UN的几何结构优化采用平面波赝势程序VASP进行,初始参数选取Pearson数据库中的实验值,对晶格参数和原子位置进行模拟,利用VASP程序包中的structure optimization选项将UN结构松弛到它的最低能量状态,电子能量收敛的标准选取为1×10−5eV,平面波截断能为400 eV,分别用广义梯度近似GGA和局域密度泛函LDA进行计算,结果如表2所列.可以看出,采用GGA的赝势优化得到的晶格参数更符合真实值.因此,接下来基于GGA赝势优化得到的结构进行声子谱的计算.

表2 晶格参数的对比Table 2.Comparison in lattice parameters.

4.2 声子谱的比较

UN声子色散关系在布里渊区不同方向的计算结果如图3所示,较低的分支为声子项,较高的分支为光子项.声子态密度的计算结果如图4所示,与声子色散关系结论一致,低频的声学支主要是体系的整体运动,光学支是原子间的相互运动,而且UN的声子项和光子项分隔较为明显.将UO2的分声子态密度与UN进行比较[19],如图5和图6所示.由图5可以看出两者的U分声子态密度相差不多,而由图6可以看出UO2中O的分声子态密度比UN中N的分声子态密度作用范围更广,但绝对值有所降低,同时UO2的声子项和光子项分隔不明显,说明声子态密度不仅与核素相关还与其晶体结构密切相关.

图3 UN声子色散关系Fig.3.Phonon dispersion of UN.

图4 UN声子态密度图Fig.4.Phonon density of states in UN.

图5 UN和UO2中U的分声子态密度Fig.5.Phonon spectrum of U in UN and UO2.

图6 UN中N和UO2中O的分声子态密度Fig.6.Phonon spectrum of N in UN and O in UO2.

4.3 比热容的比较

图7 所示为本文利用VASP+PHONON软件计算的UN比热容随温度的变化,同时与文献[20]中分子动力学(molecular dynamics,MD)的模拟结果以及文献[21]中的实验结果进行了比较.从图7可以看出,本文模拟的Cv较文献[20]中MD模拟的值更接近实验结果,同时,高温下三者都趋近于佩蒂特杜隆极限.通过比较可以看出,本文计算得到的UN声子态密度是较为准确的.

图7 UN比热容随温度的变化Fig.7.Heat capacity changes of UN with temperature.

4.4 热中子散射律、Debye Waller积分和Teff

UN中N和U的S(α,β)随β的变化情况如图8所示,图中给出了α=0.5,α=1.33和α=10的计算结果.从图8(a)中可以看出,在α较小时,S(α,β)随β变化明显,表现出显著的振动,这种振动是声子谱导致的,随着α值的增大,β的变化范围慢慢增大,同时振动减弱,这就要求在较大动量变化情况下进行计算时,β的取值范围要尽可能大一些,这样能更好地反映振子情况.从图8(b)可以看出,UN中U的S(α,β)随β变化的振荡特性不明显,这主要是因为U原子质量比较大,散射中的最大能量损失很小,这种情况下就没有必要扩展β到比较高的能量.

图8 (a)UN中N的S(α,β)随β的变化情况;(b)UN中U的S(α,β)随β的变化情况Fig.8. (a)S(α,β)of N in UN changes with β;(b)S(α,β)of U in UN changes with β.

图9 (a)293.6 K温度下UO2中16O和UN中15N(14N)的非弹性散射和弹性散射截面;(b)293.6 K温度下UO2中U和UN中U的非弹性散射和弹性散射截面Fig.9.(a)Inelastic and elastic cross sections of16O in UO2and15N(14N)in UN at 293.6 K;(b)inelastic and elastic cross sections of U in UO2and UN at 293.6 K.

由(8)式—(10)式可以得到,较高的入射能量超出了S(α,β)范围,引入SCT后,热中子的非弹性散射由Teff决定,由 (11)式—(16)式可以得到,非相干弹性散射截面由Debye Waller积分γ(0)决定.表3给出为UN中238U和14N(15N)的Debye Waller积分和Teff数值,从表3可以看出有效温度要略高于实际温度.

表3 Debye Waller积分和TeffTable 3.Debye Waller integral and Teffparameters.

4.5 热中子截面

为了比较UO2和UN热中子截面的差异,分别对比了293.6 K温度下UO2中16O和UN中14N,15N以及UO2中U和UN中U的非弹性散射截面和弹性散射截面,如图9所示.同时给出了不同温度下UN中14N以及UN中U非弹性散射截面和弹性散射截面,如图10所示.由图9和图10可以看出,UO2中16O考虑了相干弹性散射部分,UN中14N和15N忽略了弹性散射的相干部分,并且14N的非弹性散射和弹性散射截面略高于15N,14N的非弹性散射截面与16O较为接近,相同温度下UN中238U非弹性散射和弹性散射截面比UO2中238U偏小.随着温度的升高,UN中N和U的非弹性散射截面升高,中子与UN作用更激发晶格态,从而获得声子;相反,弹性散射截面随温度的升高是降低的,当能量为1 eV左右时弹性散射截面基本为零,总截面主要是非弹性散射截面的贡献,等于自由核散射截面,与自由气体模型一致.

图10 非弹性散射和弹性散射截面随温度的变化 (a)UN中14N;(b)UN中UFig.10.Inelastic and elastic cross sections changes with temperatures:(a)N in UN;(b)U in UN.

5 结论与展望

本文采用第一性原理,基于热中子散射理论分析制作了氮化铀的热中子截面,提出合适的UN热中子截面的预测方法.从UN晶格参数出发,利用第一性原理晶格动力学的直接方法得到了UN声子态密度,补充了热能区多温度点的UN热中子散射截面,填补UN热能区数据的缺失;分别将影响较大的中子热散射律、热中子弹性和非弹性散射截面与传统UO2进行对比,相同温度下,UN中238U非弹性散射和弹性散射截面比UO2小,UN中N忽略了相干散射部分;随着温度升高,UN弹性散射截面变小,非弹性散射变大;弹性散射截面在高能段等于零,非弹性散射截面在低能区主要通过声子的吸收获得能量,截面的变化符合1/v规律,在中能区,通过与UN核碰撞产生或发射声子,随着能量的升高截面增加,在高能区符合自由核模型.本文的研究结论揭示了UN在热中子反应堆中的热中子散射特性,为以UN为燃料的热中子反应堆的研究奠定了基础.