Development of Computing Code for Full Spectrum Assembly
-
摘要: 为解决先进组件设计中存在的多样几何设计问题和使用慢化材料引入的复杂能谱问题，以进一步提升组件程序设计能力，本文基于非结构网格，设计了基于2164群结构的子群共振计算与特征线输运计算的组合策略，并完成了程序开发。程序应用了高效多核素共振干涉方法、散射源移位算法和千群级别的图形处理器（GPU）特征线并行方案以保证计算效率。对不同能谱、不同几何特征的先进组件设计的验证结果表明：与蒙特卡洛基准解相比，对于快谱组件问题，特征值偏差均低于72pcm（1pcm=10−5）；对于含慢化材料的快谱组件问题，特征值偏差均低于132pcm。因此，本文设计的计算方案能够处理具有复杂几何和复杂能谱的组件问题。
-
关键词:
- 广谱 /
- 组件计算 /
- 共振计算 /
- 特征线方法（MOC） /
- 图形处理器（GPU）
Abstract: In order to address the issues caused by the diverse geometric designs and complex energy spectrum problems introduced by the use of moderator materials in advanced assembly designs, and to further enhance the design capability of the assembly code, this paper designs a combination strategy and develops a corresponding code based on unstructured mesh. The strategy combines the subgroup method, which is based on a group structure including 2164 groups, with the method of characteristic transport. To ensure computational efficiency, the code employs the efficient multi-nuclide resonance interference method, offset algorithm for scattering source and parallel scheme of graphics processing unit (GPU) characteristics at thousand-group level. Advanced assembly designs with different energy spectra and geometric features are selected to validate the proposed method. The results indicate that, compared to the Monte Carlo reference, the deviations in eigenvalue are within 72pcm (1pcm=10−5) for problems with fast spectra, and within 132pcm for problems utilizing moderators. In conclusion, the calculation scheme proposed in this paper can handle assembly problems characterized by complex geometry and complex energy spectrum. -
0. 引 言
新型反应堆燃料组件设计由于同时使用高富集度燃料和慢化材料，组件能谱呈现快热耦合，即广谱特征，这对组件计算提出了新的要求，主要包括共振计算与输运计算两方面。在共振计算方面，除需考虑传统共振核素，如238U在传统共振能段的共振自屏效应外，还需考虑以结构材料铁为代表在快中子能段出现的共振散射效应。在输运计算方面，由于需使用较为致密的能群结构以捕获组件中的中子注量率（后文简称为注量率）变化情况，所以计算负担较重。此外，先进燃料组件采用的非对称设计，使得以往快谱组件计算中使用的一维等效模型不再适用[1-2]，需开发二维输运模型，这进一步增加了输运计算的负担。
目前用于广谱共振计算的方法主要有基于等价理论的Tone方法[3-4]和子群方法[5]。Tone方法通过迭代求解超细群输运方程与超精细群固定源方程，实现共振自屏截面与特征值、注量率分布的同时更新。子群方法与之类似，以迭代求解超细群输运方程与超细群子群方程的方式同时实现共振自屏截面与特征值、注量率分布的收敛。由于两种方法需要迭代求解输运方程与固定源方程，所以计算效率较低。在输运计算方面，基于图形处理器（GPU）的硬件加速方法已应用于压水堆组件计算程序[6]，且取得了较为满意的效果，但该加速方法在千群级别下的计算效果尚未明确。
针对上述问题，本文基于在细群结构下，群内注量率变化幅度较小这一假设，使用窄共振近似源项替代精细计算源项，将共振计算与输运计算解耦，达到提升计算效率的目的；其次，为处理先进组件设计中出现的非均匀共振自屏效应，建立二维共振模型，并耦合特征线方法（MOC）进行计算；最后，将GPU硬件加速应用于共振固定源计算和输运计算当中，实现计算加速。
1. 理论模型
1.1 共振计算方法
本文选择子群方法作为共振计算方法。在子群方法中，共振自屏截面通过使用预制的子群参数表与求解子群固定源方程获得的子群注量率计算得到。对于目标能群，可通过式（1）计算共振自屏截面$ {\sigma _{x,i}} $。
$$ {\sigma _{x,i}} = \frac{{\displaystyle\sum\limits_{k = 1}^N {{w_{k,i}}{\sigma _{x,k,i}}{\phi _{k,i}}} }}{{\displaystyle\sum\limits_{k = 1}^N {{w_{k,i}}{\phi _{k,i}}} }} $$ (1) 式中，i为当前共振核素；N为子群数量；x为共振截面类型；k为当前子群数；w、σ分别为子群权重与子群截面；$\phi $为子群标注量率。结合窄共振近似，子群角注量率可由求解子群固定源方程（2）获得，通过数值积分可得到子群标注量率。
$$ {\boldsymbol{\varOmega}} \cdot \nabla {\psi _{k,i}}({\boldsymbol{r}},{\boldsymbol{\varOmega}} ) + {\varSigma _{{\text{t}},k,i}}({\boldsymbol{r}}){\psi _{k,i}}({\boldsymbol{r}},{\boldsymbol{\varOmega}} ) = \frac{1}{{4\pi }}{\varSigma _{\text{p}}}({\boldsymbol{r}}) $$ (2) 式中，r、Ω分别为位置变量与角度变量；ψk,i为子群角注量率；Σt,k,i为子群总截面；Σp为势散射截面。其中，子群总截面和势散射截面可分别由式（3）和式（4）计算得到：
$$ {\varSigma _{{\text{t}},k,i}}({\boldsymbol{r}}) = {N_i}({\boldsymbol{r}}){\sigma _{{\text{t}},k,i}}({\boldsymbol{r}}) + \sum\limits_{m \ne i}^M {{N_m}({\boldsymbol{r}}){\sigma _{{\text{t}},m}}({\boldsymbol{r}})} $$ (3) $$ {\varSigma _{\text{p}}}({\boldsymbol{r}}) = \sum\limits_{m = 1}^M {{N_m}({\boldsymbol{r}}){\sigma _{{\text{p}},m}}({\boldsymbol{r}})} $$ (4) 式中，Ni为共振核素i 的核子数密度；Nm为核素m的核子数密度；M为核素种类；σt,k,i为共振核素i第k子群的微观总截面；σp,m为核素m的微观势散射截面。
程序采用邦达连科迭代方法处理不同核素间的共振自屏效应。为提升迭代更新效率，采用Gauss-Seidel更新方案，即在一次迭代计算中，将计算得到的核素共振截面直接用于后续核素的子群注量率计算。
1.2 共振数据库制作
子群方法中，子群数据库的制作与共振计算中对源项采用的近似方法密切相关。本文使用窄共振近似来近似子群源项，因此也基于窄共振近似制作子群参数库。在共振截面表方面，考虑到在先进堆芯设计中，燃料的具体成分与栅元几何无法确定，故选择基于均匀共振截面表制作子群参数。
最终以核数据处理程序NJOY[7]制作的均匀共振截面表为基础，结合窄共振近似与帕德近似，通过求解线性方程组的方式获得子群参数库。子群参数制作的伪代码如图1所示。在能群循环中，首先基于吸收截面变化情况对该能群是否进行共振做出判断。然后，通过选取背景截面-共振截面来构造线性方程组并求解。若循环所有背景截面-共振截面组合都未满足收敛准则，则保存最优子群参数。
1.3 截面修正
在快谱问题中，(n, 2n)、(n, 3n)反应及阈能反应作为中子的主要慢化方式而不可忽略。但是确定论方法由于没有反应道的概念，无法单独表示这些反应。因此，需对散射截面进行修正，以考虑这些反应的贡献。本文选择式（5）的修正方案，对吸收截面σa、散射截面σs、散射矩阵和总截面σt进行修正。
$$ \left\{ \begin{gathered} \sigma _{\text{a}}^{^{{\prime}}} = {\sigma _{\text{a}}} - {\sigma _{({\mathrm{n}},2{\mathrm{n}})}} - 2{\sigma _{{({\mathrm{n}},3{\mathrm{n}})}}} \\ \sigma _{\text{s}}^{^{{\prime}}} = {\sigma _{{\text{es}}}} + {\sigma _{{\text{level}}}} + 2{\sigma _{({\mathrm{n}},2{\mathrm{n}})}} + 3{\sigma _{{({\mathrm{n}},3{\mathrm{n}})}}} \\ {\sigma _{{\text{s}},g - > g^{^{{\prime}}}}} = {\sigma _{{\text{es}},g - > g^{^{{\prime}}}}} + {\sigma _{{\text{level}},g - > g^{^{{\prime}}}}} + 2{\sigma _{{({\mathrm{n}},2{\mathrm{n}})},g - > g^{^{{\prime}}}}} + 3{\sigma _{{({\mathrm{n}},3{\mathrm{n}})},g - > g^{^{{\prime}}}}} \\ {\sigma _{\text{t}}} = \sigma _{\text{a}}^{^{{\prime}}} + \sigma _{\text{s}}^{^{{\prime}}} \\ \end{gathered} \right. $$ (5) 式中，σ(n, 2n)为(n, 2n)反应截面；σ(n, 3n)为(n, 3n)反应截面；σlevel为阈能反应截面；σes为弹性散射截面；${\sigma _{{\text{s}},g - > g^{^{{\prime}}}}} $为第g群至第$g^{^{{\prime}}} $群的散射反应截面。
1.4 裂变谱计算
不同的可裂变核素具有不同的裂变谱，如图2所示。因此，在实际问题中应采用问题相关裂变谱。此外，由于先进燃料组件设计中采用复杂燃料设计方案，燃料成分复杂度的增加致使以往用于压水堆制作裂变谱的调制方法不再适用。针对这一问题，本文采取求解均匀问题计算裂变谱的方案。在共振计算结束后，针对每个可裂变平源区，求解均匀输运问题，并基于式（6）归并计算裂变谱χ。
$$ {\chi _g} = \frac{{\displaystyle\sum\limits_{m = 1}^M {{\chi _{g,m}}\sum\limits_{g^{^{{\prime}}} = 1}^G {{\nu_{g^{{\prime}} ,m}}{\varSigma _{{\mathrm{f}},g^{^{{\prime}}} ,m}}{\phi _{g^{^{{\prime}}}}}} } }}{{\displaystyle\sum\limits_{m = 1}^M {\sum\limits_{g^{^{{\prime}}} = 1}^G {{\nu_{g^{{\prime}} ,m}}{\varSigma _{{\mathrm{f}},g^{^{{\prime}}} ,m}}{\phi _{g^{^{{\prime}}}}}} } }} $$ (6) 式中，νg',m为核素m在$g^{^{{\prime}}} $群每次裂变反应释放的中子数量；$ {\varSigma _{{\mathrm{f}},g^{^{{\prime}}} ,m}} $为核素m在$g^{^{{\prime}}} $群的宏观裂变截面；G为能群总数。
1.5 能群结构选择
有研究表明[3]，在快谱问题中，应用的能群数量在超过2000群时，计算结果才会收敛。考虑到先进燃料组件中的快-热谱耦合问题，本文选取2164群结构用于共振-输运计算。2164群结构由从20 MeV至0.414 eV的等对数能降的2122群和低能段来自ECCO-1968群结构[8]的42群构成。2164群结构下，对数能降随能量的变化情况如图3所示。
1.6 GPU硬件加速
在共振计算和输运计算中采用超细群能群结构增加了计算负担，需进一步提升计算效率。与单块中央处理器（CPU）相比，单块GPU具有更多计算核心，每个计算核心都可独立执行指令。这一特性与MOC的算法特点十分契合。
因此，本文选择在固定源方程求解与输运计算中使用GPU硬件加速方案。GPU特征线并行计算的伪代码如图4所示。首先，基于已知注量率计算源项或直接对源项进行赋值处理；然后，对整体计算按照辐角、长特征线、特征线线段、能群、极角进行多次分解，为每一个计算核心指定计算内容；最后，完成平源区注量率归并和特征值更新。
在输运计算中，若能群数为G，散射矩阵为满阵，则计算散射源需进行的运算次数为G2。这直接导致在千群条件下，存储散射矩阵需要的内存和散射源项计算时间占比都大幅度增加。
为提升访存效率与降低显存消耗，在存储散射矩阵时采取行压缩格式，即存储散射矩阵中每一群中第一个非0元素位置、最后一个非0元素位置和行散射矩阵数据。为解决散射源项计算时间长的问题，程序前期以原子加技术替代逻辑判断操作，这在一定程度上提升了计算效率。但原子加技术需频繁访问全局显存，这直接降低了散射源计算的并行性。因此，本文提出数据移位技术，基于散射截面所处位置进行移位操作，根据最终移位结果计算散射源项。本文将在2.3节中对散射源优化效果展开分析。
程序整体计算流程如图5所示。在共振计算完成后，读入弹性散射概率矩阵、非弹性散射矩阵和裂变谱信息，构建散射矩阵和裂变矩阵。进而通过求解系列均匀问题，获得不同区域的裂变谱信息，用于组件输运计算。
2. 数值结果
本文使用ENDF-VII.0评价库作为数据库，问题基准由蒙特卡洛程序OpenMC[9]提供，基准解粒子总数设置为106个，粒子代数为550代，其中前50代不计入统计。为验证本文共振计算方法与程序开发的正确性与通用性，选择自拟快谱组件问题与含慢化材料的快谱组件问题进行验证，选择特征值偏差与注量率分布偏差作为评估指标。
2.1 快谱组件问题
自拟快谱组件问题网格如图6所示，其中红色区域为燃料，白色区域为SS316不锈钢，蓝色区域为液态钾冷却剂。图6a为常规组件设计，图6b参考半价值组件设计[10]，将常规组件中部分燃料更换为不锈钢棒获得的新型组件设计方案。MOC计算中特征线线宽为0.03 cm，辐角离散数为48，极角离散数为3，极角求积组为T（Tabuchi-Yamamoto）求积组，采用P0阶散射近似。
表1中给出了特征值对比结果。从表1可见，本文采用的耦合MOC的子群方法具有较高计算精度，与参考解的偏差均小于72pcm（1pcm=10−5）。从图7给出的注量率分布对比结果可以看出，本文采用的共振-输运计算方法获得的注量率分布精度较高。但同时注意到，虽然在半价值组件中引入了更多的结构材料，但能谱形状未发生明显改变，具有快谱特征，需对本文计算方法的通用性进行进一步验证。
表 1 快谱组件问题特征值对比Table 1. Eigenvalue Comparison for Fast Spectrum Problems问题 OpenMC 子群方法 特征值偏差/pcm 常规组件 1.26270 1.26215 −55 半价值组件 1.04132 1.04060 −72 2.2 含慢化材料的快谱组件问题
含慢化材料的快谱组件问题来自文献[11]，对文献中的4种组件设计进行测试，网格划分如图8所示，选用20%富集度的UO2作为燃料。图中灰色区域为燃料棒，亮青色区域为CO2冷却剂，绿色区域为慢化棒，最外围为不锈钢组件框。燃料棒内部分为3圈12扇区，慢化棒分为4圈12扇区，冷却剂区域分为2圈12扇区。MOC计算中特征线线宽为0.03 cm，辐角离散数为60，极角离散数为3，极角求积组为T-Y求积组。组件输运计算在不同阶散射近似下进行，其中P0阶散射未考虑输运修正。
特征值对比结果和每种组件设计对应的造成裂变的中子平均对数能降（EALF）在表2中给出。可以看出，慢化材料的引入使得组件EALF出现明显下降，从0.1 MeV降低至1.5 keV附近。在慢化棒数量相同的条件下，慢化棒聚集会强化慢化效果，进一步降低组件EALF。慢化棒的聚集同时也强化了慢化棒附近存在的角度各向异性效应。反映到特征值偏差上，表现为P0阶散射近似条件下，紧凑组件设计特征值偏差最大，环形组件设计偏差居中，优化组件设计偏差最小。在引入高阶散射近似后，计算精度有明显提升，精度提升在散射阶数为P2时达到收敛。
表 2 含慢化材料的快谱组件问题特征值和EALF对比Table 2. Comparison of Eigenvalue and EALF for Fast Spectrum Assembly Problems Using Moderator几何 OpenMC计算特征值 EALF/MeV 子群方法计算特征值 特征值偏差/pcm P0 P1 P2 P3 P0 P1 P2 P3 常规组件 1.45504 1.191×10−1 1.45528 1.45521 1.45526 1.45523 24 17 22 19 紧凑组件 1.45181 1.016×10−3 1.46121 1.44953 1.45057 1.45049 940 −229 −124 −132 环形组件 1.42411 1.433×10−3 1.42728 1.42215 1.42283 1.42283 317 −196 −128 −128 优化组件 1.40462 1.858×10−3 1.40527 1.40312 1.40346 1.40344 65 −150 −116 −118 图9给出组件注量率对比结果。可以看到，散射近似阶数的选取对含慢化材料的组件问题影响较大。对于含慢化材料的组件，在高阶散射近似下进行输运计算，可明显提升注量率计算精度。对比三种含慢化材料的组件设计结果，可以发现，随着散射阶数提升，注量率偏差在大多能群中有降低趋势。在0.1 eV以下的热能区域，由于未存在共振现象，截面信息从数据库中直接读取得到，无法考虑真实的环境效应，致使注量率偏差出现上升趋势。整体看来，本文选用方法对含慢化材料且EALF高于1 keV的快谱组件问题有效。
2.3 计算效率分析
本节针对GPU计算效率及散射源优化效果进行分析。测试问题以2.1节中的常规快堆组件为基础，通过调整燃料圈数的方式达到调整问题规模的目的。最终选定的测试问题及测试结果如表3所示。可以看出，相较于频繁访问全局显存的原子加技术，移位技术的应用直接提升了散射源计算的并行性，使得计算效率在不同规模问题下都有明显提升，并且随着规模的提升，计算效率也在进一步上升。在平源区数量接近10000的情况下，移位技术带来接近30%的计算效率提升。同时，在4090服务器、单卡计算条件下，本节所有测试问题的共振-输运总计算时间均低于24 s，可见采用GPU作为计算平台具有较高计算效率。
表 3 不同散射源计算方法下GPU计算时间对比Table 3. Comparison of GPU Computation Time under Different Scattering Source Calculating Methods栅元
圈数平源
区数显存/
GB4090服务器 3090服务器 原子加
技术耗时/s移位技
术耗时/s原子加
技术耗时/s移位技
术耗时/s3 648 1.95 4.2 3.4 18.3 7.2 4 1272 3.02 6.9 5.3 32.9 12.7 5 2112 4.65 9.5 7.0 50.5 18.1 6 3168 6.70 13.2 9.6 72.4 24.8 7 4440 9.13 17.2 12.4 98.2 32.8 8 5928 12.00 22.0 15.6 128.3 42.3 9 7632 15.20 27.5 19.4 162.6 52.9 10 9552 18.84 33.6 23.6 201.9 64.9 3. 结 论
本文将基于窄共振近似的子群方法和GPU硬件加速方法应用于千群广谱问题，并基于二维模型进行共振自屏计算和输运计算，对自拟快谱组件问题、半价值组件问题和含慢化材料的快谱组件问题进行了测试。数值结果显示，本文采用的方法可适用于复杂几何问题和多种能谱问题。对于常规快谱组件问题和半价值组件问题，程序计算特征值偏差普遍低于72pcm。对于存在局部慢化的快谱组件设计，高阶散射近似的重要性与慢化材料的聚集程度具有正相关性。对于测试问题，散射阶数为P2时达到收敛，收敛时特征值偏差普遍低于132pcm。在效率方面，相较于原子加技术，采用移位技术可进一步提升散射源项计算效率。在平源区数为9552时，采用移位技术可提升约30%的计算效率。此外，通过观察程序在不同规模下的表现，可以得到GPU计算平台具有较高计算效率的结论，对于所有测试问题总耗时均低于24 s。
-
表 1 快谱组件问题特征值对比
Table 1. Eigenvalue Comparison for Fast Spectrum Problems
问题 OpenMC 子群方法 特征值偏差/pcm 常规组件 1.26270 1.26215 −55 半价值组件 1.04132 1.04060 −72 
下载: 导出CSV
表 2 含慢化材料的快谱组件问题特征值和EALF对比
Table 2. Comparison of Eigenvalue and EALF for Fast Spectrum Assembly Problems Using Moderator
几何 OpenMC计算特征值 EALF/MeV 子群方法计算特征值 特征值偏差/pcm P0 P1 P2 P3 P0 P1 P2 P3 常规组件 1.45504 1.191×10−1 1.45528 1.45521 1.45526 1.45523 24 17 22 19 紧凑组件 1.45181 1.016×10−3 1.46121 1.44953 1.45057 1.45049 940 −229 −124 −132 环形组件 1.42411 1.433×10−3 1.42728 1.42215 1.42283 1.42283 317 −196 −128 −128 优化组件 1.40462 1.858×10−3 1.40527 1.40312 1.40346 1.40344 65 −150 −116 −118
下载: 导出CSV
表 3 不同散射源计算方法下GPU计算时间对比
Table 3. Comparison of GPU Computation Time under Different Scattering Source Calculating Methods
栅元
圈数平源
区数显存/
GB4090服务器 3090服务器 原子加
技术耗时/s移位技
术耗时/s原子加
技术耗时/s移位技
术耗时/s3 648 1.95 4.2 3.4 18.3 7.2 4 1272 3.02 6.9 5.3 32.9 12.7 5 2112 4.65 9.5 7.0 50.5 18.1 6 3168 6.70 13.2 9.6 72.4 24.8 7 4440 9.13 17.2 12.4 98.2 32.8 8 5928 12.00 22.0 15.6 128.3 42.3 9 7632 15.20 27.5 19.4 162.6 52.9 10 9552 18.84 33.6 23.6 201.9 64.9
下载: 导出CSV
-
[1] WEI L F, ZHENG Y Q, DU X N, et al. Development of SARAX code system for full-range spectrum adaptability in advanced reactor analysis[J]. Annals of Nuclear Energy, 2022, 165: 108664. doi: 10.1016/j.anucene.2021.108664 [2] LEE C, YANG W S. An improved resonance self-shielding method for heterogeneous fast reactor assembly and core calculations[C]//Proceedings of the 2013 International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering. Sun Valley: American Nuclear Society, 2013. [3] LEE C, YANG W S. MC2-3: multigroup cross section generation code for fast reactor analysis[J]. Nuclear Science and Engineering, 2017, 187(3): 268-290. doi: 10.1080/00295639.2017.1320893 [4] JEON B K, YANG W S, JUNG Y S, et al. Extension of MC2-3 for generation of multigroup cross sections in thermal energy range[C]//International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering. Jeju, Korea: American Nuclear Society, 2017. [5] MAO L, ZMIJAREVIC I, SANCHEZ R. Resonance self-shielding methods for fast reactor calculations—Comparison of a new tone’s method with the subgroup method in APOLLO3®[J]. Nuclear Science and Engineering, 2017, 188(1): 15-32. doi: 10.1080/00295639.2017.1332890 [6] CHOI N, KANG J, JOO H G. Preliminary performance assessment of GPU acceleration module in nTRACER[C]//Transactions of the Korean Nuclear Society Autumn Meeting. Yeosu, Korea: Korean Nuclear Society, 2018. [7] MACFARLANE R E, MUIR D W, BOICOURT R M, et al. The NJOY nuclear data processing system, version 2016: LA-UR-17-20093[R]. Los Alamos: Los Alamos National Laboratory, 2017. [8] RUGGIERI J M, TOMMASI J, LEBRAT J F, et al. ERANOS 2.1: international code system for GEN IV fast reactor analysis[C]//Proceedings of the International Congress on Advances in Nuclear Power Plants. Reno: American Nuclear Society, 2006. [9] 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 [10] FEI T, MOHAMED A, KIM T K. Neutronics benchmark specifications for EBR-II shutdown heat removal test SHRT-45R - revision 1: ANL-ARC-228 Rev. 1[R]. Argonne: Argonne National Laboratory, 2013. [11] WANG L J, LU D, YAO L, et al. Study on temperature feedback effect of supercritical CO2–cooled reactor[J]. Frontiers in Energy Research, 2021, 9: 764906. doi: 10.3389/fenrg.2021.764906 -
下载:
计量
- 文章访问数: 644
- HTML全文浏览量: 59
- PDF下载量: 26
- 被引次数: 0
