2. 东华理工大学 机械与电子工程学院 南昌 330013
2. School of Mechanical and Electronic Engineering, East China University of Technology, Nanchang 330013, China
反应堆堆芯温度变化时，堆内各种材料的中子截面均随温度发生变化，从而引起反应性变化。因此，蒙特卡罗（以下简称蒙卡）程序在进行高保真反应堆模拟计算时，必须考虑温度变化对其造成的影响[1]。对于有共振行为的核素，当入射中子的能量处于不可分辨共振能区内时，共振峰会变得非常密集，由于实验测量仪器的分辨率限制，评价核数据无法给出单个共振峰的共振参数。目前蒙卡程序中使用比较广泛的方法是概率表修正法，该方法不直接使用不可分辨共振能区的截面，而是提供该能点处反应截面的概率分布。受截面的温度效应影响，对于给定的概率概率分布，其不可分辨共振截面亦因温度而发生变化。
消除上述影响的传统做法为：预先采用核截面处理软件（如NJOY的PURR模块）处理得到所需温度的不可分辨能区核截面[2-3]。但该模块处理时间较长，尤其是重核核素，其处理时间均达2 min以上，无法满足物理-热工耦合问题等需实时生成大量核素的温度相关截面的要求。为解决上述问题，随着计算机性能的不断发展，“在线”(On-the-fly)多温度核截面生成方法成为近年来研究的热点[4]。所谓“在线”，即在中子输运模拟过程中根据需要，将某一温度截面动态展宽至目标温度。目前，Walsh等[5]提出一种基于PURR模块中采用的共振序列抽样(Ladder Sampling)法在线生成不可分辨共振能区截面。该方法精度较高，且可以估计出使用概率表方法引起的不确定性，但由于实现过程类似于将NJOY的PURR模块耦合至蒙卡程序中，其计算效率较差，暂时无法满足不可分辨共振能区在线核截面生成实际需求。MCNP5 (Monte Carlo N Particle Transport Code) 1.5及以上版本附带的makxsf辅助小程序[6]采用线性插值方式生成温度相关热化截面及概率表数据，该程序采用的插值方法简单有效，但该程序只能用于截面的预处理操作，且由于输入文件存在较多限制，使得该程序使用存在诸多不便。
本文在充分调研国内外蒙卡中子输运中不可分辨共振能区的在线核截面生成方法现状基础上，根据不可分辨共振能区截面的特点，基于MCNP软件开展了蒙卡中子输运不可分辨共振区截面生成方法研究，提出使用基于内维尔(Neville)插值的不可分辨共振能区多温度截面生成方法，提高不可分辨共振能区截面生成的综合计算效率。
1 内维尔插值方法内维尔(Neville)插值法认为每个高次插值结果都是由低次的插值递推而来[7]。该算法的优点是可以直接计算出某一函数在某点的近似值，而无需函数的具体表达式，其函数在某点的高次插值结果是由前两个低次插值的结果再经线性插值得到[8]。
该算法可描述为：给定一个包含n+1点的序列(xi, yi)，并设各点的值均不相同，此时令该序列经n次插值得到的多项式P表示为：
| $ P\left( {{x_i}} \right) = {y_i}, \;\;i-0, 1, 2 \ldots, \mathit{n} $ | (1) |
引入j-i次插值多项式Pi, j(x)，使该多项式通过点(xk, yk)，其中k=i, i+1, …, j。显然，当j=i时，Pi, j(x)=P(x)。
当j-i=1时，对于任意
| $ \begin{array}{l} {\mathit{P}_{i, j}}\left( x \right) = \frac{{\left( {{x_{i + 1}}-x} \right)}}{{{x_{i + 1}}}}{y_i} + \frac{{\left( {x-{x_i}} \right)}}{{{x_{i + 1}}-{x_i}}}{y_{i + 1}}\\ \;\;\;\;\;\;\;\;\;\;\;\; = \frac{{\left( {{x_j} - x} \right)}}{{{x_j} - {x_i}}}P\left( {{x_i}} \right) + \frac{{\left( {x - {x_i}} \right)}}{{{x_j} - {x_i}}}P\left( {{x_j}} \right) \end{array} $ | (2) |
式中：插值多项式可由前一次插值多项式迭代表示。由此可得，当j-i > 1时：
| $ {\mathit{P}_{i, j}}\left( x \right) = \frac{{\left( {{x_j}-x} \right)}}{{{x_j}-{x_i}}}{P_{i, j-1}} + \frac{{\left( {x - {x_i}} \right)}}{{{x_j} - {x_i}}}{P_{i + 1}}\left( x \right) $ | (3) |
式(3)说明，对于之前每次的插值计算，增加节点时均可得到利用，而不必进行重复计算，即内维尔方法具有自适应性。
为进一步阐明该方法的自适应性，图 1表示的是当j-i=4时内维尔插值方法的迭代过程。
|
图 1 内维尔插值流程示意图 Figure 1 Diagram of the process in Neville interpolation |
不可分辨共振能区核截面插值的基本原理是：通过对事先准备好的参考温度概率表数据，在蒙卡输运过程中，通过插值得到目标温度的概率表，并由此得到相应的反应截面数值。图 2显示的是235U核素在不同温度（300 K、600 K、900 K、1200 K、1500 K与1800 K）条件下，某一概率处不同能量对应的截面，不同温度的不可分辨共振截面数据来源于JEFF3.2[9]。从图 2可以看出，不可分辨共振能区的截面在不同温度下呈现一定的线性关系，验证了插值方法的可行性。
|
图 2 不同温度条件下不可分辨共振能区截面 Figure 2 The cross sections of unresolved resonance region underin the condition of different temperatures |
由于不可辨共振能区概率表的能量框架内的能点较少，因此可按一定的温度间隔利用NJOY加工生成不同温度间隔的截面，并将概率表数据按核素名称全部提取出来，供蒙卡程序调用，如某一温度截面不在温度列表中时，采用内维尔插值方式得到目标温度截面。
由于概率表中的能量点相对较少，生成多个温度点的概率表参考库并不会占用太多存储空间。从表 1可以看出，包含265种核素、6个温度点的概率表数据大小仅为29.3 MB，处在可接受范围内。
|
表 1 不同温度概率表存储容量信息 Table 1 Probability table storage information with different temperatures |
为验证基于内维尔插值的不可分辨共振截面生成方法的准确性，图 3给出了在不同温度下，使用Neville插值方法(300 K、600 K、900 K)生成235U与238U总截面与NJOY分别进行加工所得500 K温度总截面的对比结果。结果显示：基于内维尔插值在不可分辨共振能区多温度截面与NJOY加工方法生成截面基本一致，最大偏差小于0.3%，验证了本文方法的准确性。
|
图 3 概率p=0.17413时235U核素(a)和p=0.15223时238U核素(b)不同方法生成总截面对比 Figure 3 Comparison of total cross section generationed by different methods of 235U at probability p=0.17413 (a) and 238U at probability p=0.15223 (b) |
为了宏观验证生成的概率表的可应用性，将基于内维尔插值的不可分辨共振能区的多温度概率表生成方法集成至MCNP5（以下简称MCNP_NevX，其中X表示使用温度点的数目），分别使用MCNP（含NJOY PURR加工所得目标温度概率表）与MCNP_NevX（不含NJOY加工所得目标温度概率表）进行了简单临界模型的中子输运计算，并对处理时间、有效增殖因数和能谱进行了对比与分析。对于可分辨共振能区截面，本文均假设已使用NJOY予以展宽，文中MCNP5版本为1.20，NJOY版本为99.396。
3.1 一维快堆堆芯模型Bigten1测试为验证概率表插值方法的有效性，选取了Bigten1即一维Bigten简单快堆堆芯模型，该模型由10%高浓铀组成的裸球堆芯和由贫铀组成的球壳反射层两部分组成，内外半径分别为30.48 cm和45.72 cm，堆芯布置如图 4所示。该模型总共涉及234U、235U和238U常用锕系核素[10]。
|
图 4 Bigten1球形反应堆模型 Figure 4 Spherical reactor model of Bigten1 |
计算条件为每代20000个中子，循环500代并舍去前50个非活跃代。本文对三种情况进行了计算，分别使用NJOY处理好的600 K温度截面（记为MCNP+NJOY）、使用400 K、700 K、900 K进行内维尔插值得到的600 K温度截面（记为MCNP_Nev3）与使用300 K、400 K、700 K、900 K进行内维尔插值得到的600 K温度截面（记为MCNP_Nev4）。计算处理器为Intel Core i7-6700，最终得到的计算结果如表 2所示。
|
|
表 2 Bigten1快堆堆芯模型计算结果及偏差 Table 2 Calculation result and deviations of fast reactor core of Bigten1 |
结果显示：在效率方面，使用内维尔插值方法相比与传统处理方法效率提高两个数量级；在精度方面，使用300 K、400 K、700 K、900 K温度点进行内维尔插值的结果与预处理值最为接近，相差仅为7×10-5，这也说明了原始节点越多，内维尔插值方法精度越高，但带来的影响便是计算效率相对于三个温度点内维尔插值方法有所降低。
为进一步比较概率表数据在线插值方法的精度，本文统计并比较了该模型的中子能谱，如图 5所示。可以看出，NJOY方法与在线插值方法的能谱的一致性也非常高，证明了基于内维尔插值的不可分辨共振截面生成方法的准确性，在实际反应堆分析计算中的精度与NJOY的PURR模块是基本一致的。
|
图 5 Bigten1模型能谱对比 Figure 5 Energy spectruma comparison of Bigten1 |
为进一步验证概率表数据在堆芯的中子输运计算的精度，选取了Bigten2即二维Bigten快堆堆芯模型[11]。该模型由10%高浓铀组成的圆柱体堆芯和由贫铀组成的圆柱反射层两部分组成，内外圆柱半径分别为26.67 cm和41.91 cm，内外圆柱高度分别是55.88 cm与96.52 cm，堆芯布置如图 6所示。该模型总共涉及234U、235U和238U常用锕系核素。
|
图 6 Bigten2圆柱形反应堆模型 Figure 6 Cylindrical reactor model of Bigten2 |
计算条件为每代20000个中子，循环500代并舍去前50个非活跃代。本文对三种情况进行了计算，分别使用NJOY处理好的600 K温度截面（记为MCNP+NJOY）、使用400 K、500 K、700 K进行内维尔插值得到的600 K温度截面（记为MCNP_Nev3）与使用400 K、500 K、700 K、900 K进行内维尔插值得到的600 K截面（记为MCNP_Nev4），最终得到的计算结果如表 3所示。
|
|
表 3 Bigten2快堆堆芯模型有效增殖因数计算结果及偏差对比 Table 3 Calculation results and deviations of comparison of effective multiplication factor for fast reactor core model of Bigten2 |
结果显示：效率方面，使用内维尔插值方法相比与传统处理方法效率提高两个数量级；精度方面，使用4个温度点进行内维尔插值得到的600 K概率表的结果与预处理值最为接近，仅为6×10-5，再次证明了基于内维尔插值的概率表在线生成方法的准确性。同时，本文统计并比较了该模型的中子能谱，如图 7所示。
|
图 7 Bigten2圆柱形反应堆能谱对比 Figure 7 Energy spectra comparison of Bigten2 cylindrical reactor |
可以看出，不可分辨共振区概率表在线插值方法与NJOY加工处理生成的概率表所产生能谱的一致性也非常高，其计算结果进一步证明了不可分辨共振能区概率表在线插值方法的准确性及有效性。
为了测试方法的稳定性，本文还对国际临界安全基准题库(International Criticality Safety Benchmark Evaluation Project, ICSBEP)中的一些快堆例题进行了测试，分别使用NJOY处理好的600 K温度截面（记为MCNP+NJOY）与使用400 K、700K、900 K进行内维尔插值得到的600 K温度截面计算有效增殖因数（记为MCNP_Nev3）进行对比。计算条件为每代20000个中子，循环500代并舍去前50个非活跃代。其计算结果如表 4所示。
|
|
表 4 ICSBEP测试例题keff计算结果及偏差对比 Table 4 Calculation result and deviation comparison of ICSBEP test example keff |
结果也证明，通过采用本文实现的不可分辨共振能区核截面生成方法，能够准确且高效地进行蒙特卡罗中子输运计算，证明了本文方法的可行性和有效性。
4 结语本文发展了基于插值方法的不可分辨共振能区多温度截面生成方法，提出使用内维尔插值方法在保证截面精度前提下提高截面生成效率。为了验证本文发展的方法，进行了典型核素与快堆堆芯基准例题的对比测试。计算结果表明：本文方法能够快速并准确地生成各种温度下的中子截面，在保证精度的前提下，典型核素的在线生成效率平均提高了两个数量级以上，初步验证了本文方法的正确性以及高效性，可用于反应堆多物理耦合计算。
| [1] |
谢仲生, 吴宏春, 张少泓. 核反应堆物理分析(修订本)[M]. 西安: 西安交通大学出版社, 2004, 22-23. XIE Zhongsheng, WU Hongchun, ZHANG Shaohong. Physical analysis of nuclear reactor (revised edition)[M]. Xi'an: Xi'an Jiaotong University Press, 2004, 22-23. |
| [2] |
MacFarlane R, Kahler A. Method for processing ENDF/BVⅡ with NJOY[J]. Nuclear Data Sheets, 2011, 111: 2739-2890. |
| [3] |
Macfarlane R, Muir D, Boicourt R, et al. The NJOY nuclear data processing system[R]. Version 2012. LA-UR-12-27079. Los Alamos National Laboratory, 2012.
|
| [4] |
Yesilyurt G, Martin W, Brown F. On-the-fly doppler broadening for Monte Carlo codes[J]. Nuclear Science and Engineering, 2012, 171(3): 239-257. DOI:10.13182/NSE11-67 |
| [5] |
Walsh J, Forget B, Smith K, et al. On-the-fly doppler broadening of unresolved resonance region cross sections[R]. LLNL-JRNL-704125-DRAFT. Lawrence Livermore National Laboratory, 2016.
|
| [6] |
Brown F. The makxsf code with doppler broadening[R]. LA-UR-06-7002. Los Alamos National Laboratory, 2006.
|
| [7] |
Neville E. Iterative interpolation[J]. The Journal of the Indian Mathematical Society, 1933, 20: 87-120. |
| [8] |
Press W, Teukolsky S, Vetterling W, et al. Numerical recipes[M]. 3rd ed. New York: Cambridge University Press, 2007, 118-124.
|
| [9] |
NEA Data Bank. JEFF-3.2 evaluated data library-Neutron data[EB/OL]. 2016. http://www.oecd-nea.org/dbforms/data/eva/evatapes/jeff_32/.
|
| [10] |
Frankle S. A suite of criticality benchmarks for validating nuclear data[R]. LA-13594. Los Alamos National Laboratory, 1999.
|
| [11] |
Kahler A, MacFarlane R, Mosteller R, et al. ENDF/B-VⅡ.1 Neutron cross section data testing with critical assembly benchmarks and reactor experiments[J]. Nuclear Data Sheets, 2011, 112(12): 2997-3036. DOI:10.1016/j.nds.2011.11.003 |