张宏博，赵 晨，彭星杰，赵文博，陈 长，李 庆，于颖锐，宫兆虎，曾 未，刘 琨，饶俊杰，王 博

(中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610213)

近年来，借助于高性能计算平台的长足进步，数字化反应堆的高保真、高分辨率建模与多物理模拟计算能力，对核能系统的设计改进与机理认知显现出重要且特殊的意义。其中，高保真中子学数值模拟技术作为数字化反应堆体系的核心组成之一，成为业内研究热点。近年来，仅在以特征线方法(MOC)为基础的非均匀输运理论体系框架下，就陆续涌现出了DeCART[1]、nTRACER[2]、MPACT[3]、NECP-X[4]和STREAM[5]等代表性软件。

SHARK程序是由中国核动力研究设计院正在研发的基于全堆芯确定论非均匀输运理论体系的高保真中子学计算软件。该软件采用构造实体几何(CSG)的建模方式，从多群数据库的截面与共振数据出发，采用改进子群方法刻画有效共振截面的复杂非均匀效应[6]，采用一步化的二维/一维(2D/1D)或准三维(Quasi-3D)MOC开展堆芯非均匀输运计算[7]。SHARK程序采用面向对象的设计理念，采用Python与C++混编开发，具有良好的可维护性和用户友好性。

本文介绍SHARK程序的研发情况，并采用一系列基准题对该程序进行测试验证。

1 方法模型

1.1 多群数据库

SHARK程序现阶段采用基于ENDF/B-Ⅵ.8核评价库的45群中子能群结构多群数据库[8]，这种能群结构既能较为充分地捕捉压水堆堆芯内的中子能谱特点，又不会为堆芯规模的高保真计算带来过重的计算负担，符合程序的应用堆型定位。SHARK程序数据库包含各类核素共计280余种，基本覆盖了程序目标应用对象的燃料、慢化剂/冷却剂、毒物、结构材料及相关燃耗链核素，考虑了锕系核素、裂变产物核素、可燃吸收体核素的裂变、俘获、(n，2n)/(n，3n)和衰变等主要反应类型。

1.2 全堆芯共振算法

在子群方法中，有效共振截面σx的表达形式为：

(1)

式中：wx，n和σx，n为x类型子群参数；下标a为吸收类型，n为子群编号；φn为子群通量。

在式(1)中，若要求取σx，必须求解具有输运方程形式的子群方程以获得子群通量。为提高堆芯层面共振计算效率，引入慢化方程近似解将式(1)改写为：

(2)

(3)

(4)

式中：σb，n为子群背景截面；λ为中间共振因子；Σp为宏观势散射截面；Σe，n为宏观等效截面；N为共振核素核子密度。

利用式(2)～(4)计算有效共振截面的关键是获得Σe，n，而直接求解同样是不经济的。由于Σe与吸收截面σa呈现相对平缓的单调关系，故可构造Σe(σa)函数关系线性对数插值表，以减少子群方程求解频次。

子群等效截面插值表Σe(σa)是与共振能群和共振核素相关的。根据核素共振峰分布及目标堆芯的特点，可采用共振核素分类及共振能区分段两种策略削减计算成本。

基于上述简化，以MOC求解形如式(5)的子群方程，获得子群通量并完成制表：

λΣp+(1-λ)Σpφm(r)m=1,…,M

(5)

式中：Ω为方向单位向量；r为中子向径；ψm为子群m的角通量。

当同一材料区域内含有多种共振核素时，采用迭代思想处理多核素共振干涉效应。

对于共振散射截面，可根据邦达连科迭代收敛后的有效共振吸收积分，利用二分法反向搜索背景截面，并据此计算有效共振散射截面。

对于非共振核素，采用移出截面修正因子方法相对准确刻画轻核的慢化效应，以考虑共振所致的能谱变化引发的包壳、慢化剂等材料的散射截面变化。

1.3 非均匀输运算法

SHARK程序包含2D/1D MOC和Quasi-3D MOC两种非均匀输运求解器。总体思想上，二者都是将三维输运问题转化为轴向若干层相互耦合的二维MOC问题进行求解的。目前程序主要采用2D/1D MOC。

在2D/1D MOC中，三维中子输运方程经轴向和径向积分，可得基本耦合方程为：

(6)

(7)

SHARK程序采用空间有限差分格式的离散纵标方法(SN)求解轴向一维方程(式(7))，以充分考虑轴向的非均匀输运效应。

对于2D/1D MOC数值计算的潜在不稳定性问题，采用各向同性泄漏近似和泄漏项分割方法进行弥补。

而在Quasi-3D MOC中，通过引入轴向差分格式近似，将三维中子输运方程转化为二维耦合方程：

(8)

(9)

由式(8)、(9)可见，Quasi-3D MOC不含有横向泄漏项，方程右项恒正，因此较2D/1D MOC具有更好的数值稳定性[7]。

1.4 加速与并行计算技术

对于大规模计算，加速与并行计算技术是十分必要的。加速技术方面，SHARK程序采用最优扩散粗网有限差分(odCMFD)加速方法[9]，该方法较好地改善了传统粗网有限差分(CMFD)方法的数值稳定性。

并行策略方面，SHARK程序采用基于空间区域分解的MPI+OpenMP并行方案，即采用MPI技术处理空间区域分解并行，各子区分别独立进行特征线追踪计算，并与内边界相邻子区通信(图1)；采用OpenMP技术对特征线产生和追踪过程中的部分循环体进行并行。这种混合并行策略可有效分散计算与存储负荷，具有良好的并行效率。

图1 SHARK程序空间区域分解并行策略

综合以上方法模型，图2给出程序的计算流程和基本逻辑关系。输运计算最终以特征值和平源区裂变率为收敛判定准则。

图2 SHARK程序基本计算流程

2 数值结果

SHARK程序目前已具备定态微观问题的中子学计算能力，各项测试验证均在有序推进。本文给出了部分基于国际著名基准题的计算结果。

2.1 宏观BEAVRS堆芯问题

该问题改编于装载有193盒UO2组件的大型商用压水堆堆芯[7]。堆芯中布置有1.6%、2.4%、3.1% 3种富集度的燃料组件。堆芯活性区高为380 cm，并在轴向两端各有40 cm的铁水反射层(图3)。

图3 BEAVRS堆芯几何示意图

该问题参考解来自蒙特卡罗程序RMC，以验证SHARK程序输运求解器。其中，SHARK程序采用的活性区轴向层高为9.5 cm；在每个象限，采用8个方位角和3极角Yamamoto最佳求积组；特征线宽度为0.02 cm；特征值和源项的收敛准则分别为1 pcm和10-4。BEAVRS堆芯问题计算结果列于表1，棒功率分布偏差如图4所示。由于高保真计算节点资源限制，以136核并行为基准，对比204、408及578个CPU核心数量，计算时间及并行效率列于表2。由计算结果可见，对于普通商用压水堆堆芯，非均匀输运具有较为良好的精度及效率。

表1 BEAVRS堆芯问题计算结果

表2 BEAVRS堆芯问题并行计算效率

图4 BEAVRS堆芯问题径向棒功率相对偏差分布

2.2 二维组件问题(VERA#2)

VERA#2为二维单组件基准题。该基准题较为全面地考虑了燃料富集度、温度、可燃毒物、控制棒等常见因素给商用压水堆组件带来的非均匀效应。

表3列出VERA#2基准题的计算统计结果，其中参考解来自蒙特卡罗程序KENO-Ⅵ[10]。SHARK程序在每个象限采用16个方位角和3极角Yamamoto最佳求积组；特征线宽度为0.03 cm(含IFBA问题特征线宽度为0.005 cm)；特征值与源项的收敛准则分别为1 pcm和10-5。对于特征值结果，SHARK程序的最大偏差为-220 pcm，大多数算例的偏差均在±150 pcm以内。对于棒功率分布结果，绝大多数算例的最大相对偏差(MAX)不足±0.30%，均方根偏差(RMS)不足0.10%；对于含有强吸收体的算例(典型为控制棒)，最大偏差水平有所增大。由此可见，对于二维单组件系列问题，程序的计算精度是合理和稳定的。

表3 VERA#2基准题计算结果

2.3 三维组件问题(VERA#3)

VERA#3基准题是一个压水堆三维单组件问题，包含无可燃毒物(3A@600 K)和含有16根Pyrex毒物(3B@565 K)两个算例。这两个算例改编自真实的商用压水堆组件。除径向布置，该基准题力求还原燃料组件在轴向上的端塞、压紧弹簧、管座、底座和定位格架等结构细节。因此该问题可较为充分地反映软件对压水堆轴向非均匀效应的模拟能力。

同样以KENO-Ⅵ为参考，并默认采用相同的轴向层划分，表4列出SHARK程序计算VERA#3基准题的宏观统计结果，图5示出VERA#3基准题径向棒功率偏差分布，图6示出VERA#3基准题轴向功率分布。由表4可见，特征值的计算偏差较小，径向棒功率偏差与二维组件问题的偏差水平相一致(图5)。

算例：a——3A；b——3B

表4 VERA#3基准题计算结果

对于轴向功率分布，SHARK程序的最大相对偏差仅有2.61%，且出现在顶部靠近反射层的低功率区域。同时，轴向功率平均偏差(AVG)和均方根偏差也说明了轴向功率分布的模拟结果是准确的。由图6可见，轴向功率形状与KENO-Ⅵ程序契合，且格架效应得到了准确体现。

算例：a——3A；b——3B

2.4 三维多组件问题(VERA#4)

VERA#4基准题是一个压水堆组件3×3排布问题，在径向和轴向上都有着比前述算例更强的非均匀效应。该系列问题包含3个二维算例和12个三维算例。本文对这些算例进行了精细建模与验证计算。

1)二维算例

根据中心位置RCCA组件插入情况不同，二维算例包括控制棒全提(ARO)算例4A、银铟镉(AIC)控制棒算例4B和碳化硼(B4C)控制棒算例4C。

表5列出VERA#4二维多组件基准题计算结果。由表5可见，SHARK程序能较准确地捕捉径向多组件的非均匀效应，特征值偏差不超过-43 pcm。图7示出棒功率相对偏差分布，少数临近控制棒位置的棒功率分布相对偏差最大达到2%左右(相对功率低至0.500以下)，绝大多数棒功率相对偏差不超过±1%。

表5 VERA#4二维多组件基准题计算结果

图7 VERA#4二维多组件基准题径向功率相对偏差分布

2)三维算例

与VERA#3基准题相同，VERA#4基准题同样考虑了轴向的精细结构。此外，每根控制棒在轴向上被分为两段，底部为AIC吸收体，上部则为B4C吸收体。该问题共包含12个分别处于不同控制棒位的算例，侧重于比较特征值、控制棒微分价值(DRW)、控制棒积分价值(IRW)和轴向功率偏移(AO)。这些参数能较为全面地刻画程序对功率分布的计算精度水平。

表6列出VERA#4三维多组件基准题的宏观物理参数计算结果。由表6可见，特征值偏差不超过100 pcm。对于棒价值，各棒位处DRW与IRW结果均吻合良好(图8)。同时，AO数据与参考解也十分契合。

表6 VERA#4三维多组件基准题计算结果

图8 VERA#4三维多组件基准题控制棒价值

3 总结

本文介绍了中国核动力研究设计院数字化反应堆高保真中子学程序SHARK的基本研发情况。SHARK程序采用匹配目标堆芯物理特性的多群数据库、共振计算与输运计算方法模型，同时应用了CSG建模、MPI+OpenMP并行等先进技术手段。数值结果显示，SHARK程序已具备定态微观中子学问题计算能力，对于商用压水堆相关基准题具有良好的计算精度和效率。根据项目规划，后续将对SHARK程序持续开发完善，陆续完成精细燃耗、多物理耦合、时空动力学等方法模型与代码的开发，并逐步推进更加深入和广泛的测试验证工作。