基于OpenMC的瞬发中子衰减常数计算模块开发与验证

李铸伦，谢金森,*，徐士坤，邓年彪，苑旭东，于 涛

(1.南华大学 核科学技术学院, 湖南 衡阳 421001； 2.南华大学 湖南省数字化核反应堆工程与技术研究中心, 湖南 衡阳 421001)

通过反应堆物理实验进行核设计验证是开发新型反应堆的必要环节，其中与反应性(或有效增殖因数keff)相关的物理量(如控制棒价值、停堆深度等)的实验测量及利用其对核设计程序进行校核是一项重要的工作内容。实验得到的反应性通常是以缓发中子有效份额βeff为单位的相对反应性(单位为$)，而理论计算得到的反应性是以keff换算得到的绝对反应性，两者之间的相互验证必须依赖中子动力学参数作为桥梁。随着加速器驱动次临界反应堆(ADS)向工程迈进[1]，以及核动力反应堆在次临界状态下实现控制棒刻度等实际需求出现，次临界(特别是较深次临界)绝对反应性的精确测量正逐步得到重视。

次临界反应性测量方法的一类经典方法(脉冲中子源(PNS)法、Rossi-α法等)是基于瞬发中子基波衰减常数α的。这类方法可通过式(1)得到绝对反应性ρabs，其中Λ为中子代时间。谢金森等[2-3]从中子时空动力学方程出发，指出用瞬发中子基波衰减常数α获得绝对反应性时，中子动力学参数计算中应采用瞬发α基波中子注量率。

(1)

蒙特卡罗方法具有良好的几何适应性，采用连续能量截面数据还可使其具备良好的能谱适应性，被认为是反应堆中子输运模拟的一类终极方法。利用蒙特卡罗方法进行瞬发α本征值及瞬发α基波中子注量率计算时，通常将α本征值问题转换为keff本征值问题间接求解，即k-α迭代策略。该策略通过在keff本征值问题中引入时间吸收项α/v，通过迭代寻找keff=1时的α值求解瞬发α本征值问题[4]。Cullen等[5]采用k-α迭代思想，利用TART程序对Godiva及其衍生基准题进行了模拟计算。经典的通用蒙特卡罗程序MCNP4C同样基于k-α迭代求解瞬发α本征值问题，但它对α的初始猜测值较为敏感，同时选取中子平均裂变寿命的倒数作为迭代参数，导致在中子速度较低或次临界度偏深时误差较大，甚至计算无法进行[6-7]。文献[8-10]为更好地进行反应性引入事故下的多群时空动力学计算，在OpenMC上开发了缓发α本征值计算模块，对次临界问题，该模块强制α本征值逼近负的最小缓发中子先驱核衰减常数，无法用于次临界瞬发α本征值问题计算。

针对MCNP4C因k-α迭代参数及初始α猜想值导致的次临界度较深时瞬发α本征值误差较大甚至无法计算的问题，以及文献[8-10]的工作无法实现次临界状态瞬发α本征值计算的问题，本文基于传统的k-α迭代方法，首先对迭代参数进行重新选择，其次在粒子输运过程中对瞬发、缓发中子分别考虑，基于开源蒙特卡罗程序OpenMC的二次开发，保留其计算缓发α本征值问题的功能，同时增加单独可选的瞬发α本征值问题计算功能。利用简单几何的Godiva衍生基准题及复杂几何和燃料的MUSE-4次临界实验装置对OpenMC-PA进行验证。

1 OpenMC介绍

OpenMC是由麻省理工学院(MIT)计算反应堆物理组于2011年开始研发、并在2012年12月发行第1版的并行蒙特卡罗计算程序，其主要目的是开发一种能根据研究需要易扩展的高性能、免费开源的模拟计算程序。OpenMC输入文件是一系列的XML格式文件，以便于数据在不同程序和接口之间的交换。具体包括几何结构、材料定义、计算参数设置、计数类型、绘图及加速收敛设置等输入文件[11]。

2 瞬发α本征值问题理论计算方法

2.1 理论模型

无外源增殖系统的非稳态中子输运问题可由式(2)描述：

(2)

式中：t为时间；v为中子速度(标量)；f(r,E′,Ω′→E,Ω)为散射函数，表示碰撞前中子的能量为E′，运动方向为Ω′，碰撞后中子能量变为E而运动方向为Ω的概率；φ′与φ分别为碰撞前后的中子注量率，φ′=φ′(r,E′,Ω′,t)、φ=φ(r,E,Ω,t)；r为中子群体的空间坐标；Σt为宏观总截面；Σs为宏观散射截面；Σf为宏观裂变截面；χ为中子裂变能谱；ν为每次裂变释放的平均中子数；λ为缓发中子先驱核有效衰减常数；Cj为缓发中子先驱核密度；下标j表示缓发中子分组号，p和d分别表示瞬发中子和缓发中子。

将中子注量率和缓发中子先驱核密度随时间的变化规律写成如下指数形式：

Cj(r,t)=Cα,j(r)eαt

(3)

将式(3)代入式(2)可得到广义α本征值方程[8]：

(4)

由于瞬发中子与缓发中子具有明显不同的时间特性，因此式(4)中的α本征值可分为两族：其中一族绝对值较小，即缓发α本征值；另一族绝对值较大，即瞬发α本征值。对于次临界系统，缓发与瞬发α本征值均小于0；对于临界与缓发超临界系统，缓发α本征值≥0，而瞬发α本征值≤0；对于瞬发超临界系统，缓发与瞬发α本征值均大于0。

对于次临界问题，若直接数值求解式(4)，由于占优本征值是缓发α本征值，因此最终迭代结果将返回αd，随着次临界度的加深，缓发α本征值将回归到负的最小缓发中子先驱核衰减常数，即-min{λj|j=1,2,…,J}，这符合物理上中子注量率的衰减速度不能超过寿命最长的缓发中子先驱核组的衰减常数[12]。

(5)

2.2 k-α迭代基本思想

瞬发α本征值的确定论方法迭代方案描述如下。

(6)

(7)

(8)

(9)

式中，l为迭代的代数索引。式(6)即k-α迭代的工作原理，裂变项滞后表示抽样的裂变中子被存储并用于下一次迭代时的固定源，归一化后如式(9)所示。标准的k本征值模拟计算程序原理如式(7)所示。式(8)所示的瞬发α本征值更新迭代过程是通过使用新的特征函数估计在整个相空间上对式(5)进行积分，使用归一化条件和本征值的定义来替代吸收和泄漏的积分[8]。

在式(6)中，对于任意的α都有唯一的k与之对应，即k=k(α)。将系统调整至“虚拟临界”的过程，即是寻找一个α值，可以使得k=k(α)=1，也即求解非线性方程k(α)-1=0。求解该非线性方程的方法一般是用迭代法，比较常用的迭代格式为：

αl+1=αl+cl(kl-1)

(10)

式中，cl为迭代参数，单位为μs-1，会随着迭代的进行而改变。迭代参数的选取不仅会影响收敛速度，更会影响误差估计[6]。在文献[15]中，Brockway等建议cl=αl+gn，gn为中子速度与宏观总截面乘积的期望值；而MCNP4C使用的迭代参数cl=αl+1/lf，lf为中子平均裂变寿命[16]。文献[17]对3种迭代参数(cl=αl+1/lf、cl=1/lf、cl=400)进行了比较，结果表明偏离临界越远，MCNP4C给出的统计误差越大，反映了其迭代参数的选取并不完全可靠。

在OpenMC-PA中选择的迭代参数计算为：

(11)

式中：d为中子连续两次反应间穿行的距离；s为积分变量，ds表示穿行距离上的微分量；DTl为中子在平均径迹长度上的时间吸收权重[18]，它反映了中子随机游走过程中在径迹长度上的衰减过程，更符合真实的物理过程。

3 基于OpenMC的程序实现

根据确定论方法的迭代思想，在中子输运方程中引入k作为参数计算瞬发α本征值，当k=1时，式(6)与式(5)完全一样，则k-α迭代的思想即是加入时间吸收截面[α/v]，将α本征值问题转化为k本征值问题的反问题，即通过调节α本征值对k本征值进行迭代计算，直至k本征值收敛于1，此时求得的α本征值即是系统瞬发中子衰减常数。

根据以上计算思想，OpenMC-PA总体通过两层迭代计算：内迭代完成k本征值的迭代计算，外迭代完成瞬发α本征值的迭代计算。其中，瞬发α本征值计算模块采用牛顿-拉夫逊迭代算法判定收敛条件。计算流程图如图1所示，其计算基本流程如下。

图1 基于OpenMC的α本征值计算模块流程Fig.1 Flow of α eigenvalue calculation module based on OpenMC

1) 开始内迭代，给出初始估计值kTl=1、初始值αeff=0、收敛判定误差ε=1×10-8、迭代的相对误差与标准残差均为1，确定粒子的初始能量、初始位置与初始运动方向，l=0时进行第1次迭代，总共计算L代。

2) 采用系统径迹长度估计模式统计粒子移动轨迹预估子代粒子总数Nl。

3) 确定下一个碰撞点位置并判断粒子是否穿出材料，穿出材料则粒子死亡停止追踪该粒子。

4) 根据截面确定碰撞对象。

5) 选取[0,1)范围内的随机数ξ对粒子抽样判断是否发生[α/v]反应。发生[α/v]反应且α>0则判定时间吸收，粒子死亡；发生[α/v]反应且α<0则判定δ散射，粒子能量与方向不变，权重变为原来的2倍；不发生[α/v]反应，则按照α=0进行碰撞抽样。

7) 子代粒子数与父代粒子数相除得到径迹长度估计有效增殖因数kTl。

8) 根据瞬发α的初始估计值，通过式(6)进行外迭代计算新的瞬发α本征值。

9) 通过牛顿-拉夫逊迭代算法收敛计算瞬发α本征值，更新迭代的相对误差与标准残差，并作为下一代的α初始估计值代入下一轮计算。对中子注量率进行更新。

10) 将初始αeff的值替换更新为迭代计算的α本征值，输出当前第l代的α本征值(αeff)。

11) 若当前已完成预设最大迭代代数L，则程序结束。

4 基准题问题检验

选取Godiva衍生基准题[5]以及MUSE-4基准题[19]分别对OpenMC-PA进行超临界系统以及次临界系统的模拟计算。MUSE-4基准题模型的材料及参数参见文献[19]。

4.1 Godiva衍生基准题

Godiva1模型是一个高富集度的纯铀裸球系统，Godiva2与Godiva4模型是在Godiva1的基础上更改的均匀超临界裸球系统。Godiva2是密度翻倍的快中子超临界裸球系统，Godiva4是原子密度之比为100∶1的水与铀直接混合的热中子超临界裸球系统，模型及参数列于表1。Godiva3模型是在Godiva1的基础上更改的非均匀快中子超临界系统，在燃料的外侧增加了30 cm的水作为反射层。Godiva3模型及参数列于表2。

表1 Godiva1、Godiva2、Godiva4参数Table 1 Parameter of Godiva1, Godiva2 and Godiva4

表2 Godiva3参数Table 2 Parameter of Godiva3

分别以50代非活跃中子、300代活跃中子、每代106个中子进行计算，结果列于表3。由表3可见，OpenMC计算结果与TART的计算结果吻合较好，最大相对误差小于0.5%，以简单模型验证了超临界问题中本文程序计算结果的准确性。

以Godiva1基准题模型为参考，通过持续降低质量密度使其从临界状态转为次临界状态，随着次临界度的加深，分别以50代非活跃中子、300代活跃中子、每代104个中子，在OpenMC与MCNP4C中计算该次临界问题。计算结果列于表4。

表3 Godiva衍生基准题计算结果对比Table 3 Comparison of calculation result of Godiva derivative benchmark

表4 Godiva基准题次临界计算结果对比Table 4 Comparison of subcritical calculation result of Godiva benchmark

4.2 MUSE-4次临界实验装置

MASURCA反应堆是CEA在法国卡达拉赫研究中心负责运营的零功率装置。MUSE-4是基于MASURCA改造的次临界实验装置，堆芯模型如图2所示，不同类型MUSE-4组件的轴向视图如图3所示。该装置在976组燃料组件装载时的有效增殖因数keff=0.97[19]，缓发中子有效份额βeff=0.003 4，中子代时间Λ=0.74×10-6s。MUSE-4实验装置参数列于表5。

图2 MUSE-4的栅元装载形式Fig.2 Cell loading configuration of MUSE-4

图3 MUSE-4组件1～8和11的轴向视图Fig.3 Axial view of MUSE-4 tube 1-8 and 11

对MUSE-4实验装置模型在每代10 000粒子、100代非活跃代、1 000代活跃代下进行模拟计算，结果列于表6。

表5 MUSE-4实验装置参数Table 5 Parameter of MUSE-4 experimental device

表6 MUSE-4的模拟结果与基准题对比验证Table 6 Comparison and verification of MUSE-4 simulation result and benchmark result

由表6可看出，OpenMC得到的计算结果与基准题相比相对误差均在0.5%以内，验证了OpenMC在该模型下的计算精度满足工程设计需求。同时，由OpenMC-PA得到的α，利用式(1)回归计算得到的keff为0.970 1±0.000 1，与实验值吻合良好，相对误差小于0.02%，从另一个角度验证了OpenMC-PA计算瞬发α本征值的正确性。

5 结论

本文采用k-α迭代方法，选取径迹长度上裂变中子权重的平均修正参数DTl作为迭代参数，在粒子输运过程中对瞬发、缓发中子分别考虑，在开源蒙特卡罗程序OpenMC上开发了OpenMC-PA模块，增加了单独可选的瞬发α本征值问题计算功能。基于Godiva衍生基准题以及MUSE-4次临界实验装置，分别对超临界及次临界状态下的模拟计算结果进行了验证。结果表明，在OpenMC上模拟计算得到的瞬发α本征值结果与基准题结果吻合较好，在一定次临界度下可保持稳定运行，瞬发α本征值计算结果的相对误差稳定在0.5%以内，计算范围与计算时间均优于MCNP4C，验证了开发的OpenMC-PA模块的可靠性和准确性。后续可基于OpenMC继续开展改进的k-α迭代等方法优化瞬发α本征值的计算。

感谢密歇根大学Ilham Variansyah博士对本论文研究提供的帮助。