六角形特征线方法及区域分解并行

吴文斌， 王志强， 郑竞超， 谢择怿， 赵晨， 彭星杰

(1.中山大学 中法核工程与技术学院，广东 珠海 519082； 2.中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610213)

特征线方法是一种几何适应能力强、精度高、并行潜力大的中子输运计算方法，它已成为组件程序及高保真数值堆所采用的主流方法，目前已经有一批针对栅元呈正方形排列的特征线输运计算程序[1]。然而，针对六角形几何堆芯布置，相对来说开发滞后，目前，DeCART程序[2]、SONG输运计算程序[3]、NECP-X程序[4]和ANT-MOC程序[5]等具备针对六角形几何的特征线输运求解的能力。六角形燃料组件因其排列紧凑，可以最大化冷却剂流速及其传热效率，在先进多用途反应堆的堆芯设计中广泛使用[6]。因此，开发具备六角形几何处理能力的中子输运求解器具有重要的工程实用价值。

OpenMOC程序[7]是MIT开发和维护的MOC开源程序，可用于方形堆芯和组件的特征线输运计算。本文在OpenMOC程序的基础上，开发了六角形堆芯中子输运计算程序OpenMOC-HEX,实现六角形模块化射线布置和基于区域分解的消息传递接口(message passing interface，MPI)并行计算。并采用多个基准例题对程序进行数值验证和并行效率分析。

1 六角形特征线理论模型

1.1 六角形模块化射线布置

模块化射线布置中，特征线仅在小模块(栅元或组件)内生成；小模块的短特征线通过合理的布置能够精确连接，构成长特征线。模块化射线布置方法在矩形几何的MOC程序中广泛使用[8-11]。本文参考文献[2]将模块化射线布置方法应用到六角形几何中。模块化射线布置必须满足2个条件：1)模块内短特征线在模块间精确连接；2)模块反射边界上射线与反射线精确连接。而，射线角度、射线间距的选取只有满足特定的条件，才能使上述2个条件成立。本文射线布置方式为：

2)在(0,π/6)角度范围内方位角数目N/12，全部方位角分为N/12组；

3)(0,π/6)角度范围内，用户预期的方位角平均分布为：

(1)

(2)

式中：a为六角形棱边长度；ceil为向上取整。

图1 模块化射线布置示意Fig.1 Schematic diagram of modular ray arrangement

5)模块化修正后的每组方位角的第1个角度φs与射线间距δs为：

(3)

6)每组内6个方位角的射线间距相等，每组内的6个方位角δs为：

(4)

其中，极角求积组选取Tabuchi-Yamamoto求积组[12]。

1.2 六角形栅格几何处理

为了灵活快速定义六角形组件，需要建立六角形栅格(Lattice)。如图2所示，每个栅元单元可填充不同的燃料棒。在特征线方法中，需要获取特征线段所在网格的材料参数。因此，若给定任意点在笛卡尔坐标系下的坐标(x,y)，需求出在该点所属栅元单元的编号，进而得到具体的栅元材料信息(燃料棒信息)。六角形栅格几何处理中，编号求解方法如图2所示。

图2 六角形栅格几何处理Fig.2 Hexagonal lattice geometry processing

(5)

根据计算得到(i1,i2)编号坐标，坐标点可能位于(i1,i2)、(i1+1,i2)、(i1,i2+1)、(i1+1,i2+1)4个六角形栅元单元中，需要进行进一步的判断。具体判断方法为：计算坐标点(x,y)与4个六角形栅元单元中心点的距离，距离最短者则为坐标点(x,y)所在的栅元单元。

1.3 区域分解并行实现

如图3所示，每个MPI进程都负责处理一个六角形组件，不同组件间通过角通量通信将计算耦合在一起，实现MPI并行。程序使用非阻塞模式通信，以减少通信的等待时间，避免通信死锁。为了能将角通量数据准确发送给目标进程，需找出相邻组件的对应关系。

图3 4圈六角形组件构成的堆芯MPI进程编号Fig.3 MPI ranks diagram for a core composed of four-ring hexagonal assemblies

本文对于y方向的六角形组件的进程编号(rank)从最上方六角形组件开始，沿着顺时针方向，从外向里递增。例如，4圈六角形y方向组件构成的堆芯中，MPI进程编号(rank)如图3所示。

不同圈数下，六角形的邻边对应关系不相同。为了找出邻边六角形组件的进程编号(rank)，先根据进程编号(rank)对应的(i1,i2)坐标，利用(i1,i2)坐标找到对应边的邻居六角形(i1,i2)坐标。如图4所示，再根据邻边对应六角形组件的(i1,i2)坐标找出对应的进程编号(rank)。

图4 (i1,i2)坐标法寻找邻边区域Fig.4 (i1, i2) coordinate method to find neighboring regions

2 特征线程序数值验证

2.1 C5G7 六角形堆芯基准题

使用六角形C5G7 2D基准题对程序进行验证，如图5所示，该堆芯由MOX燃料组件，UO2燃料组件和最外圈的反射组件组成，外边界为真空边界条件。

图5 六角形C5G7 2D基准题示意Fig.5 Diagram of the hexagonal C5G7 2D benchmark

图6是UO2燃料组件和MOX燃料组件的具体布置。每个六角形组件边长为11.9 cm，均采用如图7所示的栅元规格，由边长0.78 cm的六角形组成，栅元包含外径为0.54 cm的棒，用于填充燃料、裂变腔和导向管。反射组件全部由水组成，组件边长同为11.9 cm。

图6 C5G7 2D基准题燃料组件Fig.6 Hexagonal C5G7 2D benchmark fuel assemblies

图7 燃料、裂变腔和导向管栅元示意Fig.7 Schematic diagram of fuel, fission chamber and guide tube cell

特征线建模参数为60个均匀分布方位角，特征线径向间距为0.03 cm；对于燃料栅元和慢化剂栅元，网格划分参数均为径向2圈、周向4扇区；求解采用线性源近似；特征值收敛精度1×10-5，裂变率收敛精度10-5。测试平台采用AMD EPYC 7T83 64-Core Processor 处理器，主频最大2.45 GHz，内存容量512 GB，双路CPU共128计算核心256线程。

计算值keff为1.162 32，参考值为1.162 44，参考解来自蒙特卡罗程序McCARD[2]，其中功率最大棒的相对误差为-1.02%；功率最小棒的相对误差为0.17%；相对误差绝对值的最大值为1.56；相对误差绝对值的平均值为0.30；相对误差均方根为0.36；功率加权相对误差为0.26。

六角形C5G7 2D基准题棒功率分布相对误差如图8所示，其中误差为

(6)

式中P为棒功率值。

图8 六角形C5G7 2D基准题棒功率分布相对误差Fig.8 Relative error of power distribution of hexagonal C5G7 2D benchmark

由上述验证结果可知，keff计算结果与参考解相差12×10-5，燃料棒功率分布最大相对偏差为1.56%(出现在堆芯外部功率较低处)，表明程序计算精度良好。

2.2 多圈六角形组件并行效率基准题

该基准题由水、UO2和MOX燃料组成，材料的宏观截面由C5G7基准题提供。2圈组件堆芯结构如图9所示，堆芯由2圈共7个相同的组件构成。每个组件的中心和外层燃料由大栅元(图10(a))组成，材料填充为MOX-4.3%，中间层燃料由小栅元(图10(b))组成，材料填充为UO2，剩余部分用水填充。

图9 2圈六角形组件并行效率基准题示意Fig.9 Schematic diagram of parallel efficiency benchmarking of two-turn hexagonal components

图10 燃料栅元几何尺寸Fig.10 Fuel geometry

使用该基准题组件进行弱可扩展性分析。逐渐增加六角形的组件圈数进行测试，原本为2圈7个组件，增加一圈后为3圈19个组件，其中组件数Z与MPI进程数一致，对应关系为：

Z=(r-1)×r×3+1

(7)

式中r为圈数。

依次递增圈数进行弱可扩展性分析测试。测试不启用OpenMP并行，测试平台CPU为Intel Xeon Platinum 9242 CPU@2.30 GHz (单节点共96核)，使用ICC 17.0.5编译器和Intel 2017 MPI并行库。测试均为100次迭代。每种圈数计算4次，计算时间取后3次结果平均值。不同的圈数的测试结果如图11所示。

图11 弱可扩展性分析运行时间随MPI进程数变化Fig.11 Weak scalability analysis of runtime variation with the number of MPI processes

理想情况下，组件和MPI进程数增加，计算时间应不变。实际情况中，7～91进程在单节点内计算，随着进程数的增加计算时间较明显地增加；可能原因是单节点内使用核数的增加使得每个核平均可使用的缓存更少甚至主频降低，故将91进程的运行时间作为并行效率计算基准更为合理。91到1 027进程，计算时间只增加约2 s，并行效率为87%(91进程计算时间/1 027进程计算时间)，表明程序具有良好的并行扩展性和并行效率。

3 结论

1)本文基于开源特征线程序OpenMOC实现了六角形模块化射线布置和计算功能，并基于区域分解策略实现MPI多进程并行，开发了OpenMOC-HEX程序，使程序能对六角形堆芯进行计算。对基准题的测试表明：程序具有良好的计算精度，六角形C5G7 2D基准题计算解与参考解keff与参考解相差12×10-5，燃料棒功率分布最大偏差为1.56%；程序具有良好的MPI并行效率，满足计算精度和扩展性上的要求。

2)OpenMOC-HEX程序尚不具备CMFD加速功能，对于C5G7 六角形堆芯基准题外迭代次数约500次，计算代价较大。

程序暂未实现射线轴向的布置，无法计算三维六角形问题。下一步将实现基于空间区域分解的六角形CMFD加速功能和三维六角形全堆芯计算功能。