基于高阶CRAM的燃耗程序开发与验证

张彬航，杨森权，陈云龙，洪 锋,*，袁显宝，张永红，唐海波，冯虹瑛

(1.三峡大学 机械与动力学院，湖北 宜昌 443002;2.中核武汉核电运行技术股份有限公司，湖北 武汉 430074)

燃耗计算是反应堆设计及燃料管理的重要环节，也是反应堆燃耗信任制技术应用及乏燃料后处理的基础[1]。燃耗计算的核心是求解燃耗方程，由于堆芯燃料在中子辐照过程中产生大量的裂变产物，且不同核素间的半衰期差别极大，导致燃耗方程组具有大型、稀疏、刚性强的特点。目前求解的方法主要包括线性链解析方法、常微分方程组差分方法和矩阵指数法。线性链解析方法将堆内复杂的核素转换关系拆分成一系列线性反应链后进行求解，具有较高的计算精度，但为了平衡计算精度和效率，往往需要引入截断阈值舍弃重要性较低的反应链。同时对于深燃耗问题，已有研究表明线性链解析方法存在较大的计算误差[2]。常微分方程组差分方法则是通过采取特殊的数值计算技巧处理燃耗方程组的刚性问题，计算速度较快，但计算精度对时间步长的选取较为敏感，代表程序为FISPACT[3]。矩阵指数法主要包括泰勒展开与截断方法、切比雪夫有理近似方法(CRAM)等。ORIGEN程序采用泰勒展开与截断方法进行燃耗方程求解，具有速度快、步长包容性好的特点，但需对短寿命核素单独处理后近似求解[4]。CRAM具有速度快、精度高且不需单独对短寿命核素进行处理的特点，在很多燃耗程序中得到了广泛的发展与应用，包括Serpent[5]、RMC[6]和JMCT[7]等。但目前已有研究表明，CRAM对于在长时间步长下核素系统的衰变计算精度较差，难以满足高保真燃耗计算的发展需求[8-9]。

本文基于高阶CRAM研制点燃耗程序ICRAM，同时内耦合于蒙特卡罗输运程序OpenMC，形成燃耗计算分析程序OPICE，通过一系列基准题的计算，验证OPICE的计算精度和准确性。

1 高阶切比雪夫有理近似方法

在中子辐照下，堆内燃料中核素成分随时间的变化可用燃耗方程组描述，其矩阵形式如下：

(1)

式中：N(t)为t时刻核素的核子密度向量；A为系数矩阵，由系统中所有核素的反应率和衰变率组成。因此，可得矩阵方程(式(1))的解：

N(t)=eAtN(0)

(2)

(3)

(4)

N(t)≈α0N(0)-

(5)

式中，I为单位矩阵。

目前基于CRAM发展的大部分燃耗程序皆使用了PFD形式的有理近似展开，当阶数k取到14阶或16阶时就能获得较为理想的计算结果。但如果核素核子密度在某一燃耗步长内急剧下降，那么PFD形式的CRAM计算精度就会下降或导致计算错误[12]。因此，本文基于CRAM发展了不完全局部分解(IPF)形式的高阶有理近似展开方法来满足高保真燃耗计算的精度需求。

ml=1,2,…，μ

(6)

式中：ml为低阶有理函数的阶数；μ为整数。

(7)

(8)

N(t)=eAtN(0)≈

(9)

基于上述高阶CRAM研制了点燃耗程序ICRAM，其中阶数k的取值为16、32和48。由式(5)和(9)可看出，在相同阶数下，PFD和IPF形式的CRAM计算量相当，且在求解过程中涉及多次系数矩阵的求逆，对发展大规模高保真燃耗计算的效率会存在一定影响，因此对ICRAM的效率优化将是未来研究的重点。

2 蒙特卡罗输运-燃耗耦合计算方法

蒙特卡罗输运-燃耗耦合计算实质上是蒙特卡罗中子输运计算与点燃耗计算过程的相互耦合。按照耦合方式的不同，分为外耦合和内耦合两种。外耦合方法通过第三方接口程序和外部存储数据来实现蒙特卡罗输运程序和点燃耗程序之间的数据传递，导致计算效率低，计算规模受限。新近开发的蒙特卡罗燃耗程序多采用内耦合的方式，即在蒙特卡罗程序中嵌入点燃耗计算模块，能减少用户工作量，提高计算效率。本文在完成点燃耗程序ICRAM研制的基础上，进一步将OpenMC(V0.10.0)与ICRAM进行内耦合开发，形成了燃耗计算程序OPICE。

OPICE的简要计算流程如图1所示。首先通过中子输运计算得到各燃耗区中子通量密度、单群反应率等中子学参数。然后经过加工处理后传递给燃耗计算模块并对各燃耗区执行点燃耗计算，得到步长末各燃耗区的核素种类及核子密度。最后更新相应燃耗区核素的核子密度并传递给输运计算进行下一时间步长的耦合计算，依次迭代直至燃耗寿期末。OPICE是基于OpenMC进行输运计算来获得不同燃耗区的中子通量密度和单群反应截面。而在蒙特卡罗输运计算中，中子通量的计数结果对源粒子进行了归一化处理，并不能真实反映燃耗区的中子通量密度。因此需根据当前燃耗步长内的热功率进行换算得到每个燃耗区的真实中子通量密度。假设堆芯共有N个燃耗区，每个燃耗区内有M个核素，热功率为P，计算公式如下：

(10)

式中：φreal为真实中子通量密度；φstatic为归一化中子通量密度；Vj为燃耗区j的体积；Ri和Qi分别为核素i的裂变反应率及可利用能。为进一步提高耦合计算精度，OPICE采用了预估-校正和子步法两种耦合策略。

图1 OPICE燃耗计算流程图Fig.1 Flow chart of burnup calculation of OPICE

同时OPICE还支持点燃耗计算功能，此时OPICE采用精细燃耗数据库进行点燃耗计算。该数据库整合了ORIGEN-S和ORIGEN-2的最新燃耗数据库，包含了1 487个核素、23种中子截面反应、11种衰变反应及30种重核的裂变产额。对于特定的24种锕系核素和32种裂变产物核素给出了精确的裂变能和中子俘获能，其余核素则分别认为裂变能近似为200 MeV，中子俘获能近似为5 MeV，以满足中子通量密度与功率的精确转换[15]。

3 数值验证与分析

3.1 高阶CRAM的数值精度

为验证分析高阶CRAM的数值精度，测试用例选取了富集度为3.4%的UO2燃料棒进行测试，计算该燃料棒辐照至燃耗深度为1 MW·d/kg(U)后，再冷却100 d的核素核子密度的变化情况。参考解选用了基于MATLAB符号工具箱和高精度计算工具箱研制的TTA算法的计算结果。在OPICE中，分别使用了16、32、48阶PFD和IPF两种算法进行计算。为保证计算效率，核素核子密度的截断值设为1.0×10-30cm-3，两种算法相对于TTA算法计算结果的相对偏差如图2所示。

如图2a所示，当阶数k同取16时，PFD形式的计算精度比IPF形式的计算精度稍差。随着阶数k提高至32，由图2b可见IPF形式的计算精度相比于PFD形式有了显著提高，大部分核素核子密度的相对偏差在10-5以内，其中核素核子密度大于10-6cm-3的所有核素的相对偏差均在10-13以内。当阶数k取至48时，IPF形式计算得到的核素核子密度大于10-20cm-3的所有核素的相对偏差在10-10以内，其中大部分核素核子密度的相对偏差在10-14以内，其精度优势进一步扩大，表明高阶CRAM的计算结果是足够精确的。值得注意的是，对于传统PFD形式，随阶数的提高，计算精度并不会显著提高。其中一个重要原因可能是在PFD形式中，相应的留数会随阶数的提高而增大，导致计算结果对舍入误差更敏感，计算精度下降[11]。

为评估高阶CRAM的步长包容性，将上述测试例题中UO2燃料棒的冷却时间步长分别取至10、103和106a。OPICE分别采用了16、32、48阶的PFD和IPF形式的CRAM进行计算，核素核子密度截断值取1.0×10-20cm-3，相对偏差及单步计算时间列于表1。由表1可见，在相同时间步长下，IPF形式的计算精度优于PFD形式。同时随时间步长的增加，高阶CRAM的计算结果充分体现了其高精度性和优异的步长包容性，单步时间步长可达到106a。注意到表1中对于48阶PFD形式的大步长计算，核素相对偏差小于10-10的核素百分比显著低于16阶和32阶的计算结果。其可能原因如前文所述，随阶数的提高，高阶PFD形式中留数的增大会引起其计算结果对舍入误差的敏感性增强，导致数值精度下降。同时随时间步长的增长，更多核素的核子密度低于截断值，对表1的统计结果没有贡献。后续工作还需进一步采用更多例题进行测试分析。在计算效率方面，单步计算时间随阶数的提高近似成倍增加。计算效率主要由两方面因素影响：一是随阶数的提高，计算过程中系数矩阵的多次求逆导致了耗时略长；二是为保证计算精度，ICRAM采用了多精度浮点数据库MPFR进行浮点数运算，使得计算时间增加[16]。因此计算效率优化将是下一步研究工作重点。

a——16阶；b——32阶；c——48阶图2 IPF和PFD形式计算结果的相对偏差 Fig.2 Relative deviation of calculation result of IPF and PFD forms

表1 高阶CRAM的步长包容性计算结果Table 1 Calculation result of step size tolerance of high-order CRAM

3.2 OECD/NEA压水堆栅元燃耗基准题

OECD/NEA压水堆栅元燃耗基准题来自于欧洲经济合作与发展组织核能署(OECD/NEA)发布的燃耗信任制系列基准题之一，其计算目标是压水堆组件中的单个栅元，并比较不同程序系统对燃料中同位素成分进行燃耗计算和分析的能力。OECD/NEA发布了16个机构共计21套程序的计算结果，包括对乏燃料有效增殖因数有着重要影响的核素质量等结果。该基准题包含4个完整的运行循环，各循环的燃耗时间、冷却时间及硼浓度列于表2。根据功率和最终燃耗深度的不同，分为工况A(27.35 GW·d/tHM)、工况B(37.12 GW·d/tHM)和工况C(44.34 GW·d/tHM)。具体的几何与材料参数详见文献[17]。

表2 运行循环参数Table 2 Parameter of burnup operation

OPICE采用ENDF/B-Ⅶ数据库进行计算，表3列出工况B下核素质量的计算结果。由表3可见，OPICE对大部分核素质量计算结果的相对偏差在5%以内，吻合良好。另外部分核素如237Np、149Sm与实验值的相对偏差较大，但通过对比可发现，这2个核素的计算结果落在其他程序的结果范围内。由表3还可看出，21套程序的计算结果存在一定范围，主要原因是不同程序进行输运计算时所使用的方法不同，对中子通量密度和反应率的计算会造成一定的误差。同时参考程序所使用的数据库主要来自ENDF/B-Ⅴ、ENDF/B-Ⅵ等，在进行燃耗计算时，大部分程序直接使用少群截面库或ORIGEN-2的数据库，对最终的计算结果也会造成一定的误差。总体来看，OPICE具有良好的燃耗计算精度。

表3 工况B下OPICE核素质量的计算结果Table 3 Calculation result of isotopic mass of OPICE under condition B

3.3 快堆燃耗基准题

该基准题是OECD/NEA发布的一个二氧化钚燃料的快堆基准题，其计算目标是1个快中子增殖堆的等效模型。该快堆R-z几何结构如图3所示，包括堆芯内部区域、堆芯外部区域、控制棒孔道区域和反射层区域，其中不同富集度的MOX燃料组件被布置在堆芯内、外两个区域。该堆的热功率为1 500 MW，负载因子为0.8，具体的几何、材料和功率参数详见文献[18]。

图3 快堆基准题R-z几何结构图Fig.3 R-z geometry sketch for fast reactor benchmark

不同研究机构发布的计算结果包括了有效增殖因数、核素核子密度等物理量随燃耗的变化。此外从文献[19]中找到了OMCB的计算结果。OPICE采用ENDF/B-Ⅶ数据库进行计算，在燃耗过程中设置了3个有效满功率天(EFPD)作为重要时间点，分别为0、375和625 EFPD。表4列出不同机构反应率亏损率的计算结果。由表4可见，OPICE的计算结果落在了其他机构发布的计算结果之间。若以发布结果的均值作为参考值，则在375 EFPD和625 EFPD时，OPICE对应的差异为-2.28%和-1.48%，最大差异来自PNC(J3.2)的-2.81%和-3.10%。表5进一步列出堆芯燃耗区中重要核素在燃耗周期内核素质量的变化情况。以各机构的平均值作为参考值，可看出OPICE的计算结果同样在各机构的范围内，且具有较小的相对偏差，证明了OPICE对快谱燃耗系统具有较好的计算精度。

表4 快堆基准题中初始keff和反应率亏损率的对比Table 4 Comparison of initial keff and reactivity loss ratio for fast reactor benchmark

表5 快堆基准题中核素质量的变化对比 Table 5 Comparison of isotopic mass variation for fast reactor benchmark

4 结论

本文基于高阶CRAM研制了点燃耗程序ICRAM，并与蒙特卡罗输运程序OpenMC内耦合形成了一套燃耗计算分析程序OPICE。在点燃耗算法的研究中，高阶IPF形式的CRAM相比于传统PFD形式的CRAM，具有数值稳定性好、计算精度高、步长包容性更好等特点，更符合高保真燃耗计算的需求。通过数值验证与分析可发现，高阶的CRAM单步时长可达到106a，并仍可保持很高的计算精度。在耦合方式选择中，OPICE采用内耦合的方法避免了传统外耦合导致计算效率低、计算规模受限等问题，并采用了预估-校正和子步法两种耦合策略保证计算精度。选用了OECD/NEA压水堆栅元燃耗基准题和快堆燃耗基准题对OPICE进行验证，程序计算结果与实验值及各程序的参考值吻合良好，初步证明了程序的正确性和有效性。后续工作中将重点对高阶CRAM中的系数矩阵求逆效率、数据存储结构等方面进行优化，同时开展大规模并行燃耗计算方法研究，以保证OPICE在大规模高保真燃耗计算中的求解效率，使其具有更加普适的工程应用价值。