Massively Parallel Simulation and Optimization of Advanced Nuclear Reactors with Dispersed Particle Fuel
-
摘要:
颗粒燃料是将核燃料制成颗粒并弥散在基体中的一种新型燃料构型，广泛应用于高温气冷堆、空间堆、氟盐冷却高温堆等先进堆型中. 以高温气冷堆和空间堆为例，基于开源蒙特卡罗程序OpenMC研究了适用于颗粒燃料临界计算的虚拟网格模拟加速方法，并在山河超算平台开展了超10万核心的大规模并行测试. 结果表明，高温气冷堆模型的有效增殖因数计算结果与石岛湾核电站实验数据符合较好，验证了程序及模型的准确性. 在性能方面，虚拟网格方法与OpenMC此前的真实网格方法相比，在存储空间和计算速度上均有明显提升，高温气冷堆虚拟网格模型的内存和耗时分别为真实网格模型的0.2%和82%；此外，由于虚拟网格方法简化了模型几何，其间接实现了更好的负载均衡，使得程序拥有了更高的并行效率. 对于强可扩展性，在10752核规模的测试中，虚拟网格的并行效率为83.4%，而真实网格为63.6%；对于弱可扩展性，虚拟网格模型在131600核并行效率为83.1%，而真实网格为66.1%.
Abstract:Dispersed particle fuel is a new type of nuclear fuel that is in shape of small spheres and dispersed in a matrix. It has been widely used in advanced reactors such as high-temperature gas-cooled reactors (HTRs), space reactors and fluoride-salt-cooled high-temperature reactors. This study, taking an HTR and a space reactor as examples, develops a virtual lattice method based on the open-source Monte Carlo code OpenMC to speed up dispersed particle fuel criticality simulation. Parallel tests on the scale of 100000 cores are carried out on Shanhe supercomputing platform. The keff result of the HTR model agrees well with Shidao-Bay nuclear power plant experiment, indicating that the code is of high accuracy. As for the performance of the code, results show that the virtual lattice model is of less memory footprint and higher computational efficiency than the original physical lattice model. The memory-consumption and time-consumption of the HTR virtual lattice model are 0.2% and 82% that of the physical lattice model respectively. And thanks to the simplification of geometry, the virtual lattice model is of higher parallel efficiency. For strong scalability, the parallel efficiency of the virtual lattice model with 10752 cores is 83.4% while that of the physical lattice model is 63.6%. And for weak scalability, the parallel efficiency of the virtual lattice model with 131600 cores is 83.1% while that of the physical lattice model is 66.1%.
-
Keywords:
- nuclear reactor /
- dispersed particle fuel /
- Monte Carlo /
- massively parallel /
- lattice speed up
-
随着能源和环境问题日益严峻，社会亟需发展区别于传统化石能源的新能源. 其中，核能由于其稳定性好、能量密度大、碳排放量低的优点而成为解决能源与环境问题的有力手段. 此外，由于以上优点，核能还是航天器、潜艇、水面舰艇等特种设备的良好供能方式. 因此，核能存在着较大的发展机遇.
由于核能实验研究成本较高，故计算机模拟计算成为了核能领域的重要研究方式. 其中，以反应堆堆芯核反应率分布为主要研究对象的“反应堆物理”计算有助于获得反应堆有效增殖因数（keff）、功率分布等关键参数，并为堆芯热工计算、核电厂系统计算等其他模拟研究提供必要的参数，是核能模拟计算研究的重要组成部分. 蒙特卡罗（Monte Carlo，MC）方法[1]是使用计算机生成随机数来模拟中子等微观粒子在堆芯内的运动、散射、吸收等随机行为的计算方法，其数学模型较为贴近堆芯实际的物理过程，故计算结果较为精确. 但同时MC方法需要对堆芯的大量粒子分别进行模拟，因此也存在着计算量大的缺点. 而由于堆芯中各粒子的行为是相对独立的，故MC方法可以将不同粒子直接分配给不同线程（进程）进行计算，这样的并行模式天然具有非常好的并行可扩展性，与之相比，基于区域分解的并行模式的并行效率随并行规模的增加而显著降低[2]，因此是一种非常适合并行计算的方法. 目前，MC方法已经在反应堆物理计算领域得到了广泛地应用[3-6].
当前，随着核技术的发展，一种新型的核燃料：颗粒燃料逐渐得到了越来越广泛的应用. 诸如球床式高温气冷堆[7]、空间反应堆[8]、氟盐冷却高温堆[9]、耐事故燃料[10]等多种先进核能装置中均使用了这种燃料. 颗粒燃料是将核燃料活性成分，例如UO2制成小颗粒，然后弥散在基体中，而非传统棒状燃料使用整块均匀的UO2. 可见，颗粒燃料的使用必将导致堆芯几何中存在大量随机球体. 以一种高温气冷堆HTR-PM（high-temperature gas-cooled reactor pebble-bed module）为例，其堆芯存在约42万个随机堆积的球形燃料元件，每个燃料元件中又包含约1万个随机排列的燃料颗粒[7]，如图1所示，因此其全堆总共存在约40亿个随机排列的球体. 如此复杂的几何结构会为反应堆物理计算，尤其是MC计算带来极大的困难，大大增加其计算耗时. 例如经本文验证，在相同的MC计算规模（即模拟的中子数量，根据统计学大数定律，计算的中子数越多，结果越精确，标准差越小）以及并行规模下，如果不使用任何加速方法，一个HTR-PM的球形燃料元件的计算耗时将达到一个传统压水堆棒状燃料元件的400倍以上. 因此有必要开展适用于颗粒燃料几何的加速方法研究.
本文以其中几何较为复杂的球床式高温气冷堆HTR-PM[7]为例，在开源MC程序OpenMC[11]的基础上进行加速方法的研究.
HTR-PM是由清华大学核能与新能源技术研究院主导研发的球床式高温气冷堆，具有第4代先进核能系统[12]特征. 因其具有良好的固有安全性，HTR-PM受到了学界的广泛认可. 其示范堆山东石岛湾核电站属于我国16个国家科技重大专项之一，目前已经完成了建设并实现了并网发电. 因此，针对HTR-PM的研究将具有很高的工程应用价值.
OpenMC是由美国麻省理工学院主导开发的开源MC中子输运程序. 除麻省理工学院外，阿贡国家实验室、密歇根大学等多个研究机构也对其开发工作做出了贡献. OpenMC程序具有较高的几何建模和数据统计灵活性，并具备高性能并行能力，支持MPI（message passing interface）和OMP（open multi-processing）混合并行. 作为开源程序，OpenMC自2011年问世以来已经在世界范围内进行了大量的反应堆物理计算，其各项功能得到了学界的充分验证，具有较高的可靠性. 因此OpenMC适合用于开展本项研究.
除HTR-PM外，本文还建立了一种使用颗粒燃料的空间反应堆模型，以对优化后的程序进行辅助验证.
1. 模型建立
1.1 随机球体位置的生成
为建立颗粒燃料的高精度模型，必须准确设定每个燃料颗粒的位置. 对于HTR-PM中的数十亿燃料颗粒来说，通过实验方法测量其位置显然是不现实的，因此颗粒的位置也必须通过计算机模拟计算得出.
OpenMC程序本身具备生成随机颗粒位置的功能. 其使用的是RSP（random sequential packing）和CRP（close random packing）结合的方法[13]. RSP是在规定的空间内随机选取球体位置并去除重叠的球体的方法. 这一方法在填充率较低时拥有较高的计算效率，但其无法生成38%以上填充率的随机颗粒，且其生成的颗粒是均匀分布在空间内的，无法模拟重力作用下的球堆. CRP方法则可以生成填充率更高的随机颗粒，但其计算效率较低，同时也无法模拟重力的影响（由于本文不涉及CRP方法，故对其原理不再赘述）. OpenMC在生成随机颗粒时，会根据目标填充率进行判断，如果填充率较低，则使用RSP方法，否则将使用CRP方法. 除此之外，OpenMC还支持用户从外部输入颗粒的位置，为其与其他程序的耦合提供了可能.
HTR-PM堆芯中包含2重随机球体，首先是堆芯容器中堆积的燃料球，其次是每个燃料球中弥散的燃料颗粒，如图1所示. 在建立OpenMC模型时，可以首先建立单个燃料球及其内部的燃料颗粒的模型，然后使用“重复结构”功能将该燃料球重复地堆积在堆芯容器内. 如此一来，HTR-PM模型中的随机球体就分为了“堆芯容器中的燃料球”和“每个燃料球中的燃料颗粒”2部分. 因此建模时无需生成一个40亿量级的球堆，而只需生成一个40万量级和一个1万量级的球堆即可. 对于单个球内的燃料颗粒，其填充率约为7%，且均匀分布在燃料球基体中，因此非常适合使用OpenMC自带的RSP方法. 而燃料球则是在重力的作用下堆积在堆芯容器内，且填充率约为61%，因此不适合直接使用OpenMC进行计算，而需要单独进行处理.
本文使用了DEM（discrete element method）程序LAMMPS[14]来生成燃料球的位置信息. DEM方法是使用经典力学模拟每个颗粒运动过程的方法. 在LAMMPS模型中，首先设置好堆芯容器的几何，然后在高处生成燃料球，再令其在重力作用下自然下落并堆积在容器内即可. 由于这一物理模型较为接近HTR-PM的实际装料过程，因此可以得到较为准确的结果. LAMMPS模型的具体参数与文献[15]中的一致.
值得一提的是，HTR-PM模型并非包含42万个燃料球的平衡堆芯模型，而是包含约5万个燃料球和28万个石墨球的初装堆模型[16]，这是因为当前HTR-PM示范堆山东石岛湾核电站投产不久，并未达到平衡堆芯状态，而仅有初装堆的运行数据. 因此，建立初装堆模型可以与石岛湾核电站的实验数据进行对比，以更加有力地验证程序的准确性.
初装堆状态下，堆芯的燃料球和石墨球并非均匀混合，而是先在堆芯底部堆积约23万个石墨球，然后在其上方混合堆积约5万个燃料球和5万个石墨球. 由于LAMMPS程序只关注球的经典力学特征，不关注其核物理性质，而燃料球和石墨球的经典力学特征，如质量、摩擦系数、弹性系数等相近，因此这2种球在LAMMPS模型中不作区分. 在LAMMPS给出球的位置后，再人工确定其中哪些是石墨球、哪些是燃料球，最后将数据输入到OpenMC中.
LAMMPS程序得出的球床堆积率为61%，与文献[16]中的符合较好.
与HTR-PM不同，本文模拟的空间反应堆[17]堆芯为包含1519个圆形冷却剂孔道的六棱柱，其几何中仅存在1种随机球体，即弥散在堆芯固体区域内的大量燃料颗粒. 该堆芯中颗粒填充率约为30%，且均匀分布在基体中，因此使用RSP方法计算即可. 但是，由于全堆燃料颗粒数量极多（约3亿个），故直接使用OpenMC在整个堆芯内部随机生成颗粒位置是不实际的，此外OpenMC还不支持在带有冷却剂孔道的复杂几何中生成随机颗粒位置. 综上所述，本文在建立空间堆模型时同样使用了“重复结构”功能：将一个冷却剂孔道周围的六棱柱区域视为一个单元，使用RSP方法在此实心六棱柱内部随机生成燃料颗粒，并人工去除生成在冷却剂孔道中的颗粒位置，然后将1519个这样的单元重复填充在堆芯中，即可实现空间堆随机颗粒的建模.
1.2 HTR-PM与空间堆全堆模型
在1.1节所述的生成随机颗粒位置方法的基础上，本文根据石岛湾核电站初装堆相关参数[16]，建立了HTR-PM堆芯的OpenMC模型，如图2所示. 并根据文献[17]建立了空间堆堆芯的OpenMC模型，如图3所示.
2. 加速方法
2.1 网格加速的原理
前文已经提到，大量随机球体的几何结构会大大增加MC方法的计算量. 究其原因，是程序无法对随机排列的球体进行有效地检索. 例如在图4（a）所示的情形中，大量随机排列的燃料颗粒空间中存在一个中子，而程序需要确认中子具体位于哪个颗粒内. 为此，程序只能将中子的坐标一一带入颗粒的几何信息中进行验证，如式（1）所示：
|{\boldsymbol{r}}_{0}-{\boldsymbol{r}}_{i}|\le {R}_{i}^{} \text{，} (1) 其中 {\boldsymbol{r}}_{0} 是中子位置， {\boldsymbol{r}}_{i} 是第i个燃料颗粒的球心位置， {R}_{i}^{} 是第i个燃料颗粒的半径. 如果对于某一燃料颗粒i，式（1）得到满足，则中子位于该颗粒内或其表面，而如果中子不位于任何燃料颗粒内，则中子位于基体中.
显然，由于模型中颗粒的数量极多，以上“遍历所有颗粒”的过程必然耗费大量时间，因此有必要进行优化.
针对这一问题，“网格加速”方法是一种有效的手段. 如图4（b）所示，在空间中划分规则的网格，使得每个单元格的范围为
x\in ({x}_{0}+\left(i-1\right){P}_{x},{x}_{0}+i{P}_{x}) \text{，} y\in ({y}_{0}+\left(j-1\right){P}_{y},{y}_{0}+j{P}_{y}) \text{，} (2) z\in \left({z}_{0}+\left(k-1\right){P}_{z},{z}_{0}+k{P}_{z}\right) \text{，} 其中， ({x}_{0},{y}_{0},{z}_{0}) 为网格空间角点坐标； {P}_{x} ， {P}_{y} ， {P}_{z} 分别为x，y，z方向的网格栅距；i，j，k分别为x，y，z方向的网格编号.
完成网格划分后，即可建立单元格与燃料颗粒的对应关系，将与某单元格相交的所有燃料颗粒的几何信息与该单元格编号 (i,j,k) 关联（该过程仅在程序初始化时进行，计算过程中无需重复）. 如此，则当进行中子位置的判断时，程序即可首先根据式（2）求出中子所在的单元格 (i,j,k) ，然后仅遍历该单元格内的所有燃料颗粒即可. 例如图4（b）所示的情形中，只需对3个实线的燃料颗粒进行验证，而无需验证虚线的颗粒，这就大大节省了计算时间. 以上即为网格加速方法的基本原理.
2.2 虚拟网格与真实网格
当前OpenMC中已经支持网格加速方法，但其使用的是“真实网格”，而本文为OpenMC添加的是“虚拟网格”的方法，本节将介绍二者的不同.
在真实网格模型中，网格分界面是模型中真实存在的曲面，模型将被网格面划分为若干长方体. 如图5（a）所示，使用真实网格方法将模型划分为左右2个单元格，则最终的模型是由左右2个长方体拼接而成的.
而虚拟网格则是一套假想的网格，其模型几何中不存在网格面，不会对模型进行分割. 如2.1节所述，网格加速方法本质上是建立一种各球体与其所在区域的对应关系，以便根据粒子位置缩小所需遍历的球体数量. 因此，网格面的存在是没有必要的，只需按照网格位置对各球体进行编号即可. 例如图5（b）所示的情形下，使用虚拟网格划分模型，则只需对左侧4个球体添加编号“1”，代表其属于左侧单元格；对右侧4个球体添加编号“2” ，代表其属于右侧单元格；对中间球体同时添加编号“1，2”，代表其同时属于2个单元格. 模型的几何不存在变化.
虚拟网格方法与真实网格相比，具有4个优点：
1）虚拟网格对存储空间的占用更小. 这是因为真实网格的模型几何由大量的单元格几何体拼接而成，其复杂度远超虚拟网格模型的单一几何体. 以图2所示的HTR-PM堆芯模型为例，为达到良好的加速效果，需要做出较细的网格划分，例如将堆芯容器内堆积的燃料球以及石墨球划分为33×33×154的网格，同时将单个燃料球内的燃料颗粒划分为24×24×24的网格. 显然，如果使用真实网格，则如此大规模的单元格必将占据大量存储空间. 表1展示了该网格划分下虚拟网格与真实网格模型占用的存储空间对比，可见使用虚拟网格方法将大大节约存储空间.
表 1 HTR-PM模型存储空间占用Table 1. Storage Space Occupation of HTR-PM Model项目 虚拟网格 真实网格 输入文件大小/MB 86 407 预处理占用内存/GB 1.8 6.2 运行占用内存/GB 0.409 213 2）虚拟网格的计算速度更快. 具体表现在2方面. ①在程序初始化阶段，虽然虚拟网格模型需要额外对球体进行编号，但由于其模型几何大大简化，使得初始化的整体耗时远小于真实网格模型. ②在程序计算阶段，由于虚拟网格模型不存在网格面，故程序无需对中子穿过这些网格面的过程进行模拟，因此虚拟网格模型也将在计算阶段节约时间. 表2展示了相同计算规模和并行规模下，HTR-PM堆芯模型的计算时间对比. 可见虚拟网格模型的初始化时间和计算时间分别约为真实网格模型的13%和82%，拥有更高的计算效率.
表 2 HTR-PM模型计算耗时Table 2. Time-Consumption of HTR-PM Models 项目 虚拟网格 真实网格 初始化时间 19 150 计算时间 2545 3112 表3则展示了空间堆模型的计算时间对比，其虚拟网格的初始化时间和计算时间分别为真实网格的13%和46%，与HTR-PM模型相比效率的提升更为明显.
表 3 空间堆模型计算耗时Table 3. Time-Consumption of the Space Reactor Models 项目 虚拟网格 真实网格 初始化时间 52 390 计算时间 1600 3464 3）虚拟网格的数据统计更加便捷. 由于真实网格模型存在网格面对模型进行分割，其难免会切断模型中随机分布的球体，使得一个完整的球在模型中实际由2个，甚至多个“半球”拼接而成. 例如图5（a）中，中间的球体被网格面分割为了2部分. 因此如果用户需要以单个球为单位统计某参数（例如功率分布），其必须人工将多个“半球”指定为一个整体进行统计，这为用户操作带来了困难. 而虚拟网格模型中不存在网格面，其所有球均为一个完整的球体，因此用户可以便捷地进行数据统计.
4）使用虚拟网格加速的模型拥有更好的并行可扩展性. 这是因为虚拟网格模型几何更简单，使得单个中子输运过程的耗时更短、随机性更小. 由于在并行计算中，程序只有在每个进程各自完成其“负责”的全部中子的计算后，才会进入下一阶段，故更小的随机性就意味着进程间相互等待的时间更短，即程序拥有更好的负载均衡. 图6（a）展示了HTR-PM虚拟网格和真实网格的并行强可扩展性对比，其中每代中子数保持10752000不变，总代数为20，并行规模为112~10752核；图6（b）为弱可扩展性对比，其中每代中子数为核数的1000倍，总代数为20，并行规模为112~131600核. 本文针对图6的每个数据点均进行了10次重复测量，并标出了其平均值以及10次测量结果的变化范围. 其中某次测量的并行效率是用测量得到的计算时间与基点（112核）下10次测量的平均时间计算得出的. 由于数据的波动性，某些“并行效率”可能超过了100%.
![]()
图 6 虚拟网格与真实网格并行可扩展性对比Figure 6. The parallel scalability comparison of the virtual lattice and physical lattice可见，在强扩展问题中，当并行规模达到约3000核时，虚拟网格的并行效率显著高于真实网格；而对于弱扩展问题，当并行规模达到约3000核时虚拟网格的并行效率略高于真实网格；而达到约10000核时该趋势则更为显著. 因此可以认为虚拟网格拥有更好的并行可扩展性. 此外，真实网格测量结果的波动范围明显远大于虚拟网格，这可能是因为真实网格数据量更大、中子历史更复杂，其在缓存读写等方面具有更大的不确定性.
2.3 虚拟网格加速在OpenMC程序的实现
在OpenMC程序中，首先需要对初始化部分进行修改，以便生成虚拟网格. 具体的方法为：根据式（2）的规则建立单元格编号与其位置的对应关系，并建立对应的3维数组，使得数组编号与单元格编号一一对应. 然后对于每一个单元格，遍历所有随机球体，如果某一球体与单元格存在交集，则程序将该球的信息，包括其编号、球心位置、半径储存在单元格对应的数组元素内. 在对所有的单元格完成遍历后，就成功建立了燃料颗粒的几何信息与单元格编号的关联，即生成了虚拟网格.
除生成虚拟网格外，还涉及计算过程中对虚拟网格的利用问题，需要对OpenMC进行修改.
OpenMC计算中的几何处理过程包含3项主要任务：
1）已知中子的位置，求其位于哪个几何体内；
2）已知中子的位置及其运动方向和运动距离，求其即将与哪个曲面相撞；
3）已知中子正在穿过某一曲面，求其即将进入哪个几何体.
其中，任务1已经在2.1节中进行了讨论，此处不再赘述.
对于任务2，由于任意曲面都可以表示为 f(x, y,z)=0 ，因此在已知中子当前位置 {\boldsymbol{r}}_{0} 和其运动方向 {\boldsymbol{r}}_{\mathrm{d}\mathrm{i}\mathrm{r}} （单位向量）后，只需求解关于运动距离 l 的方程：
f\left({\boldsymbol{r}}_{0}+l\cdot {\boldsymbol{r}}_{\mathrm{d}\mathrm{i}\mathrm{r}}\right)=0\text{，} (3) 即可求得中子运动到该曲面上所需的距离，如果方程没有正实数解，则中子不会与该曲面相撞，此时程序会将碰撞距离 l 记为正无穷. 在模型空间的所有曲面中， l 值最小的曲面即为即将与中子相撞的曲面. 此处值得一提的是，假如最小相撞距离大于中子本次运动的距离，则中子不会与任何曲面碰撞，而是会在当前所处的几何体内发生核反应.
可见，这项任务的计算量同样会受到遍历的曲面数量的影响. 在大多数情形下，所需遍历的曲面数量不多，无需进行优化，而当“中子位于基体内，求其即将与哪个球体相撞”时，需要对所有球面进行遍历，其计算量过大，需要借助虚拟网格加以简化. 简化的方法与任务1类似，程序首先可以根据中子的位置找出中子所在的单元格，然后仅遍历该单元格内的所有球面即可，如果这些球面均不满足碰撞要求，则只需对下一个单元格内的球面进行遍历，如此循环直至找到即将碰撞的球面，这样就有效减少了需要检验的曲面数量.
任务3的逻辑则较为简单. 在中子穿过曲面离开某几何体时，程序只需遍历与该几何体相邻的所有几何体，判断中子是否在其内部即可. 但这一过程同样存在计算量问题：当中子穿过球体表面进入基体时，程序需要判断中子是否在基体空间内，这需要对组成基体空间的大量曲面（这些曲面包括基体内的所有球面）一一进行验证，极为耗时. 而实际上这一过程是没有必要的，中子穿过球体表面离开球体时必然会进入基体. 因此本文针对这种特殊情况进行了优化，取消了其判断过程，直接认定中子进入了基体，即可显著减小计算量.
表4展示了在56个进程的纯MPI并行下，对于高温气冷堆单个燃料球模型（中子历史为100代，每代1万个）使用不同优化水平时的计算耗时对比. 可见对3项任务的优化均可使计算效率得到明显提升. 值得一提的是，如果使用真实网格方法对该模型进行加速，其计算耗时为65 s，可见其计算耗时高于虚拟网格方法.
表 4 不同虚拟网格优化下HTR-PM单个燃料球模型计算耗时对比Table 4. Comparison of the Time-Consumption of Single Pebble Model in HTR-PM Under Different Virtual Lattice Optimizations优化水平 计算耗时/s 加速比/% 无优化（不使用任何网格） 5500 仅优化任务2 2900 90 优化任务2, 3 800 588 优化全部任务 50 10900 而对于HTR-PM全堆模型（中子历史为1000代，每代10万个），在相同并行规模下（10个MPI进程，每个进程56个OMP线程）与不进行优化相比，优化全部任务后其计算耗时可以从估算耗时约10700 h减小到539 s，计算效率大幅度提升，加速效果更为显著.
3. 结果分析
3.1 程序运行环境
本文所做的计算均在国家超级计算济南中心的山河超级计算平台开展，其采用Intel Xeon 6258R型号28核CPU处理器，主频2.7 GHz，每个计算节点包含2个CPU，其内存为192 GB. 使用的操作系统为Linux version 3.10.0-957.el7.x86_64，编译所用的软件为cmake/3.21.1，gcc/9.4.0，hdf5/1.8.13，mpich/3.4.2.
需要特别指出的是，表1所示真实网格算例的内存占用（213 GB）超过了山河超算平台单个节点的内存，该算例实际是在北京超算平台测试的. 而本文其他真实网格模型则略微缩小了网格规模，使其能够在山河平台运行.
3.2 计算结果分析
本节将介绍HTR-PM和空间堆模型的定量化计算结果，以验证优化后程序的性能及其正确性.
有效增殖因数（keff）可以描述核反应堆的临界程度，是衡量反应堆物理计算结果的重要指标. 山东石岛湾核电站于2021年成功临界，其临界时keff =1. 而本文按照其参数建立的HTR-PM虚拟网格模型的keff计算结果为0.99854±0.00010，二者相差146 pcm. 同时，如果使用OpenMC程序本身的真实网格方法计算，则其keff=0.99851±0.00010，与虚拟网格仅相差3 pcm. 可见，无论是与实验值对比，还是与OpenMC程序本身对比，HTR-PM虚拟网格模型的结果都符合得很好.
而本文模拟的空间堆还处于理论阶段，不存在相应实验装置. 同时，文献[17]并未对空间堆燃料颗粒进行显式建模，而是将材料混合后建立了均匀燃料模型. 研究表明，这样的处理方法会使计算结果存在很大偏差[18-20]，不具有参考价值. 综上所述，本文的空间堆虚拟网格模型缺少实验或计算结果作为验证基准，因此其只能与OpenMC自身的真实网格方法进行对比验证. 经计算得到：虚拟网格模型的计算结果为keff=1.10547±0.00007，真实网格的结果为keff=1.10551±0.00007，二者符合较好.
综上所述，可以认为使用虚拟网格优化后的程序具有良好的准确性. 图7展示了HTR-PM的堆芯功率分布，其中白球为石墨球，无功率. 得益于2.2节所述的虚拟网格数据统计方面的便利性，该功率分布结果可以精确到单个燃料球.
图8展示了空间堆模型的堆芯径向功率分布，结果可精确到全堆尺寸的1/8000.
3.3 并行性能分析
3.3.1 并行策略优化
OpenMC支持MPI和OMP混合并行，而在并行总核数相同的情况下，不同的MPI进程数和OMP线程数搭配得到的计算效率不同. 因此，有必要寻找一种最优的并行策略. 表5展示了同样使用单节点56核并行的情况下，HTR-PM虚拟网格模型在不同线程、进程数搭配时的计算时间，其中每种搭配分别测量了10次结果，计算了其平均值、标准差.
表 5 不同并行策略下HTR-PM模型计算耗时比较Table 5. Comparison of the Time-Consumption of HTR- PM Model with Different Parallel StrategiesMPI进程数 每个进程的OMP线程数 计算时间/s 1 56 277±3 2 28 215±4 4 14 204±2 8 7 204±3 14 4 205±2 28 2 214±1 56 1 226±2 由表5可见，当进程数由1变为2时，计算效率存在明显提升. 这可能和山河超算的硬件条件有关，其单个56核的节点实际包含了2个28核的CPU，因此相比单个进程的设置，2个进程的设置可以简化CPU之间的通信，进而节省计算时间. 此后随着进程数的增加，计算效率存在先升后降的趋势，这是因为MPI并行虽然计算效率更高，但其通信也更为耗时，因此过少或过多的MPI进程数设置都不利于获得更高的计算效率.
考虑到相比于“2个进程、28个线程”的并行策略，其他设置对效率的提升不显著，且更多的MPI进程设置会大量占用节点内存，不利于大规模数据统计，因此本文最终选取的并行策略是：每个节点2个进程、每个进程28个线程.
3.3.2 大规模并行效率分析
本文借助山河超算集群对HTR-PM模型进行了大规模并行计算，测试了其在10万核规模下的并行性能. 图9展示了HTR-PM模型的强、弱可扩展性，其计算规模与2.2节相同：对于强可扩展问题，每代中子数保持10752000不变，总代数为20，并行规模为112~10752核（2~192节点）；对于弱可扩展问题，每代中子数为核数的1000倍，总代数为20，并行规模为112~131600核（2~2350节点）. 图9中统计的“计算时间”是指排除程序初始化以及写输出所用时间后的运行时间. 对于每个数据点，本文均进行了10次测量，并选取了其平均值进行分析.
![]()
图 9 HTR-PM虚拟网格模型的并行可扩展性Figure 9. Parallel scalability of HTR-PM virtual lattice model当计算规模不随并行规模变化时，程序在10752核并行下与112核并行相比具有83. 4%的并行效率，其理想计算时间与实际计算时间分别为46 s和55 s. 而当计算规模与并行规模成正比时，程序在131600核并行下与112核并行相比具有83. 1%的并行效率，其理想计算时间与实际计算时间分别为46 s和56 s. 可见程序在大规模并行时仍然具有较高的并行效率，其并行可扩展性良好.
此外注意到在弱扩展问题中，当并行规模由14336核（256节点）增大到28672核（512节点）时，程序的并行效率突然下降，由95%降为84%，而此后并行效率随并行规模的下降反而不明显. 目前分析这一现象可能由MPI通信协议的变化引起. OpenMC程序的通信耗时主要花费在进程间的中子源同步（synchronizing fission bank）上，在弱可扩展问题中，虽然每个进程的中子数是固定的，但进程间所需传递的中子数量仍会随进程总数的增加而不断增大[21]. 而MPI在通信的信息量较大时，通信协议会随数据量有所调整，从而导致计算耗时的阶跃式增加[22]. 针对这一现象，在后续研究工作中将开展更加深入的分析和测试.
4. 总结与展望
本文在开源MC程序OpenMC的基础上，开发了适用于颗粒燃料的虚拟网格加速方法. 通过构建虚拟网格并对中子输运过程的3项主要任务进行优化，有效提高了程序的计算效率. 该方法的计算耗时不但远低于不使用网格的算例，且与OpenMC此前的真实网格方法相比，计算耗时也更小. 针对HTR-PM模型和空间堆模型，虚拟网格的计算时间分别减小到了真实网格的82%和46%. 同时虚拟网格相较真实网格还具有内存占用小、数据统计方便、并行可扩展性好的优点.
在定量化结果方面，HTR-PM模型的keff计算结果与石岛湾核电站实验值符合较好，因此可以认为程序具有良好的准确性. 除此之外，本文还计算了HTR-PM和空间堆的功率分布，并得益于虚拟网格方法数据统计的便利性，实现了HTR-PM精确到单个燃料球的堆芯功率分布计算.
最后，本文对优化后的程序进行了大规模并行测试，结果表明：对于强扩展问题，使用虚拟网格加速的HTR-PM模型在10752核并行下与112核并行相比，并行效率为83.4%；而对于弱扩展问题，程序在131600核并行下与112核并行相比并行效率为83.1%，因此可以认为程序具有良好的并行可扩展性.
本文所建模型及测试更侧重于对程序并行性能的探讨，模型本身的计算并不需要10万核规模的并行. 在此后的研究中，本文优化的程序将可以用于真正需要大规模并行的、更加复杂的模型中，例如带有燃耗的HTR模型.
同时，本文程序的并行性能仍存在优化空间. 例如进行燃耗计算时，OpenMC在数据读取和输出部分的并行可扩展性较差，导致其占用了大量时间. 这些问题均有待进一步的研究.
作者贡献声明：李睿涵负责算法优化的实现与部分模型的建立以及论文的撰写；侯叶凡负责部分模型的建立；李玉辉负责开展大规模并行计算；刘召远指导并行计算工作以及修改论文；张海红指导建模工作；梁金刚指导算法优化以及修改论文.
-
![]()
图 6 虚拟网格与真实网格并行可扩展性对比
Figure 6. The parallel scalability comparison of the virtual lattice and physical lattice
![]()
图 9 HTR-PM虚拟网格模型的并行可扩展性
Figure 9. Parallel scalability of HTR-PM virtual lattice model
表 1 HTR-PM模型存储空间占用
Table 1 Storage Space Occupation of HTR-PM Model
项目 虚拟网格 真实网格 输入文件大小/MB 86 407 预处理占用内存/GB 1.8 6.2 运行占用内存/GB 0.409 213 
下载: 导出CSV
表 2 HTR-PM模型计算耗时
Table 2 Time-Consumption of HTR-PM Model
s 项目 虚拟网格 真实网格 初始化时间 19 150 计算时间 2545 3112
下载: 导出CSV
表 3 空间堆模型计算耗时
Table 3 Time-Consumption of the Space Reactor Model
s 项目 虚拟网格 真实网格 初始化时间 52 390 计算时间 1600 3464
下载: 导出CSV
表 4 不同虚拟网格优化下HTR-PM单个燃料球模型计算耗时对比
Table 4 Comparison of the Time-Consumption of Single Pebble Model in HTR-PM Under Different Virtual Lattice Optimizations
优化水平 计算耗时/s 加速比/% 无优化（不使用任何网格） 5500 仅优化任务2 2900 90 优化任务2, 3 800 588 优化全部任务 50 10900
下载: 导出CSV
表 5 不同并行策略下HTR-PM模型计算耗时比较
Table 5 Comparison of the Time-Consumption of HTR- PM Model with Different Parallel Strategies
MPI进程数 每个进程的OMP线程数 计算时间/s 1 56 277±3 2 28 215±4 4 14 204±2 8 7 204±3 14 4 205±2 28 2 214±1 56 1 226±2
下载: 导出CSV
-
[1] Raeside D E. Monte-Carlo principles and applications[J]. Physics in Medicine and Biology, 1976, 21(2): 181−197 doi: 10.1088/0031-9155/21/2/001
[2] 郭红,李艳,安恒斌. 氧碘化学激光器数值模拟中的多块并行通信算法[J]. 计算机研究与发展,2016,53(5):1166−1172 doi: 10.7544/issn1000-1239.2016.20148444 Guo Hong, Li Yan, An Hengbin. A parallel communication algorithm in supersonic COIL’s calculations using multiblock mesh[J]. Journal of Computer Research and Development, 2016, 53(5): 1166−1172 (in Chinese) doi: 10.7544/issn1000-1239.2016.20148444
[3] Sood A, Forster R A, Archer B J, et al. Neutronics calculation advances at Los Alamos: Manhattan project to Monte Carlo[J]. Nuclear Technology, 2021, 207: S100−S133 doi: 10.1080/00295450.2021.1956255
[4] Mohsen M Y M, Hassan M S, Aziz M, et al. Investigating the neutronic, thermal-hydraulic, and solid mechanics analysis for AP-1000 nuclear reactor [J/OL]. Energy Sources Part A—Recovery Utilization and Environmental Effects, 2021[2022-07-23].https://www.tandfonline.com/doi/full/10.1080/15567036.2021.1912215
[5] Al Zain J, El Hajjaji O, El Bardouni T, et al. Neutronic and burn-up calculations of the (ThO2-UO2) pin cell benchmark using DRAGON5 and MCNP6.2 codes with ENDF/B-VIII. 0 nuclear data library[J]. International Journal of Energy Research, 2021, 45(8): 11538−11551 doi: 10.1002/er.6460
[6] Rao Junjie, Shang Xiaotong, Yu Ganglin, et al. Coupling RMC and CFD for simulation of transients in TREAT reactor[J]. Annals of Nuclear Energy, 2019, 132: 249−257 doi: 10.1016/j.anucene.2019.04.039
[7] Zhang Zuoyi, Dong Yujie, Li Fu, et al. The Shandong Shidao Bay 200 MWe high-temperature gas-cooled reactor pebble-bed module demonstration power plant: An engineering and technological innovation[J]. Engineering, 2016, 2(1): 112−118 doi: 10.1016/J.ENG.2016.01.020
[8] Li Zeguang, Sun Jun, Liu Malin, et al. Design of a hundred-kilowatt level integrated gas-cooled space nuclear reactor for deep space application [J/OL]. Nuclear Engineering and Design, 2020[2022-07-24].https://www.sciencedirect.com/science/article/pii/S0029549320300649?via%3Dihub
[9] Jiang Dianqiang, Zhang Dalin, Li Xinyu, et al. Fluoride-salt-cooled high-temperature reactors: Review of historical milestones, research status, challenges, and outlook [J/OL]. Renewable & Sustainable Energy Reviews, 2022[2022-07-24].https://www.sciencedirect.com/science/article/pii/S1364032122002581?via%3Dihub
[10] Awan MQ, Cao Liangzhi, Wu Hongchun, et al. Neutronic design study of a small modular IPWR loaded with accident tolerant-fully ceramic micro-encapsulated fuel[J]. Nuclear Engineering and Design, 2018, 335: 18−29 doi: 10.1016/j.nucengdes.2018.04.023
[11] Romano P K, Horelik N E, Herman B R, et al. OpenMC: A state-of-the-art Monte Carlo code for research and development[J]. Annals of Nuclear Energy, 2015, 82: 90−97 doi: 10.1016/j.anucene.2014.07.048
[12] Zohuri B. Nuclear Reactor Technology Development and Utilization [M]. Amsterdam: Elsevier, 2020
[13] Jodrey W S, Tory E M. Computer-simulation of close random packing of equal spheres[J]. Physical Review A, 1985, 32(4): 2347−2351 doi: 10.1103/PhysRevA.32.2347
[14] Thompson A P, Aktulga H M, Berger R, et al. LAMMPS — A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales [J/OL]. Computer Physics Communications, 2021[2022-05-20].https://www.sciencedirect.com/science/article/pii/S0010465521002836?via%3Dihub
[15] Rycroft C H, Dehbi A, Lind T, et al. Granular flow in pebble-bed nuclear reactors: Scaling, dust generation, and stress[J]. Nuclear Engineering and Design, 2013, 265: 69−84 doi: 10.1016/j.nucengdes.2013.07.010
[16] She Ding, Xia Bing, Guo Jiong, et al. Prediction calculations for the first criticality of the HTR-PM using the PANGU code [J/OL]. Nuclear Science and Techniques, 2021[2022-05-15].https://link.springer.com/article/10.1007/s41365-021-00936-5
[17] 杨谢. 兆瓦级空间棱柱堆概念设计与安全分析[D]. 北京:清华大学,2019 Yang Xie. Conceptual design and safety analysis of megawatts space prismatic reactors[D]. Beijing: Tsinghua University, 2019 (in Chinese)
[18] Auwerda G J, Kloosterman J L, Lathouwers D, et al. Effects of random pebble distribution on the multiplication factor in HTR pebble bed reactors[J]. Annals of Nuclear Energy, 2010, 37(8): 1056−1066 doi: 10.1016/j.anucene.2010.04.008
[19] Wang M J, Peir J J, Sheu R J, et al. Effects of geometry homogenization on the HTR-10 criticality calculations[J]. Nuclear Engineering and Design, 2014, 271: 356−360 doi: 10.1016/j.nucengdes.2013.11.062
[20] She Ding, Liu Zhihong, Shi Lei. An equivalent homogenization method for treating the stochastic media[J]. Nuclear Science and Engineering, 2017, 185(2): 351−360 doi: 10.1080/00295639.2016.1272363
[21] Romano P, Forget B. Parallel fission bank algorithms in Monte Carlo criticality calculations[J]. Nuclear Science and Engineering, 2012, 170(2): 125−135 doi: 10.13182/NSE10-98
[22] Gropp W, Lusk E, Doss N, et al. A high-performance, portable implementation of the MPI message passing interface standard[J]. Parallel Computing, 1996, 22(6): 789−828
计量
- 文章访问数: 173
- HTML全文浏览量: 16
- PDF下载量: 46
