基于蒙卡方法的高温气冷堆随机堆积球床模拟研究

Simulation Study of Randomly Packed Pebble Bed in High-Temperature Gas-Cooled Reactor Based on Monte Carlo Method

摘要: 对于采用球形元件的高温气冷堆，燃料球与石墨球经堆芯顶部加料管装入后，在堆芯中随机堆积成球床。为实现对球床堆物理特性的准确掌握，开展精确的物理仿真十分必要。本文以球床高温气冷堆为研究对象，首先采用离散元方法(DEM)和程序获取球床堆芯随机堆积下所有燃料球和石墨球的空间坐标，随后采用蒙卡程序构建出球床堆芯真实的随机堆积模型以进行有效增殖因子k_eff计算，此外，本文还基于“蒙卡重复结构规则排布 + 燃料球和石墨球随机分配”的方法构建了球床堆芯近似的规则排布模型。对于10兆瓦高温气冷实验堆(HTR-10)初次临界问题的计算结果表明：球床随机堆积模型和规则排布模型的k_eff计算结果都与实验结果较为接近(偏差小于0.4%)，说明本文建立的两种球床模型都具有较高的计算精度；当燃料球占比增加和球床高度增加时，两种球床模型的k_eff计算结果更为接近；相较于球床随机堆积模型，规则排布模型的蒙卡计算时间更短，并且易于满足蒙卡程序对于球形元件模拟规模要求，更推荐其作为工程计算模型。

Abstract: In high-temperature gas-cooled reactors using spherical fuel elements, fuel pebbles and graphite pebbles are loaded from the top of the core through the charging tube and randomly packed into a pebble bed. To accurately understand the physical characteristics of the pebble bed, conducting precise physical simulations is essential. This paper takes a pebble-bed high-temperature gas-cooled reactor as the research subject. First, the spatial coordinates of all fuel and graphite pebbles in the randomly packed core are obtained using the Discrete Element Method (DEM) and corresponding programs. Then, a Monte Carlo program is employed to construct an actual random packing model of the pebble bed core for calculating the effective multiplication factor (k_eff). Additionally, this paper develops an approximate regular arrangement model of the pebble bed core based on the method of “Monte Carlo repeated structures with regular arrangement + random assignment of fuel and graphite pebbles”. The calculation results for the initial criticality problem of the 10 MW high-temperature gas-cooled reactor (HTR-10) indicate that the k_eff results from both the random packing model and the regular arrangement model are relatively close to the experimental results (with deviations less than 0.4%), demonstrating that both pebble bed models established in this study have high computational accuracy. When the proportion of fuel pebbles increases or the height of the pebble bed rises, the k_eff results from the two models become even closer. Compared to the random packing model, the regular arrangement model requires shorter Monte Carlo computation time and is easier to meet the simulation scale requirements of Monte Carlo programs for spherical elements, making it more recommended as an engineering computational model.

1. 引言

球床模块式高温气冷堆是国际公认的具有优异安全特性的先进堆型。该堆型采用氦气作为冷却剂，堆芯使用由包覆颗粒燃料构成的球形燃料元件，并支持在不停堆条件下进行连续换料。为实现对球床堆物理特性的准确掌握，开展精确的物理仿真十分必要。然而，大量球形元件(燃料球和石墨球)在堆芯内呈随机堆积的空间分布，这给确定论方法的建模分析带来了很大困难。以往对于球床堆芯的物理建模，通常采用蒙特卡罗程序的重复结构功能来实现，即将燃料球和石墨球的随机堆积近似为重复结构规则排布，以控制蒙卡程序的栅元数量在其能力范围内(小于10万个栅元)并提高计算效率。相关研究主要包括：在清华大学的经荥清和杨永伟等人对10兆瓦高温气冷实验堆(HTR-10)初次临界装料预估研究中[1]，采用确定论程序VSOP基于球床堆芯近似处理模型计算了不同球形元件数量下堆芯有效增殖因子k_eff，最终预估HTR-10初次临界所需的混合球数量与实验值16,890个较为接近，误差小于1%；在清华大学的张竞宇和李富等人基于蒙特卡罗方法对VSOP模型验证研究中[2]，基于球床模块式高温气冷堆核电站(HTR-PM)堆芯比较了VSOP程序和蒙特卡罗程序的计算结果，结果表明VSOP模型中不同的近似处理方法会带来不同的偏差，但是与最精细的蒙特卡罗模型在有效增殖因子k_eff方面的计算结果差别不大，从而验证了VSOP模型的可靠性；在清华大学的常鸿和杨永伟等人对HTR-10初次临界物理计算的蒙特卡罗方法模型分析研究中[3]，采用两种蒙特卡罗程序基于四种不同的重复结构栅元建立了球床堆芯模型，四种模型计算的临界混合球数和初次临界试验误差范围为[−4.6%, +1.3%]。

对于实际上随机堆积的球床堆芯，基于规则几何或规则排布的近似方式虽然便于程序建模，但是会引入不同程度的计算偏差，因此需要发展一种高精度方法来真实模拟球床堆芯的随机堆积。

本文以HTR-10初次临界问题为研究对象，首先采用离散元方法(DEM)和程序获取球床堆芯随机堆积下所有燃料球和石墨球的空间坐标，随后采用蒙卡程序构建出球床堆芯真实的随机堆积模型以进行有效增殖因子k_eff计算，并通过与实验结果的对比来验证基于离散元方法模拟球床堆芯随机堆积的准确性。此外，考虑到蒙卡程序对于球形燃料元件模拟规模和效率的限制，本文还研究了便于工程应用的“蒙卡重复结构规则排布 + 燃料球和石墨球随机分配”的球床堆芯模拟方法，并分析了其k_eff计算精度。

2. HTR-10结构特点

2.1. 反应堆结构

Figure 1. Physical calculation zoning model of initial fuel loading for HTR-10

图1. HTR-10初装堆物理计算分区模型

Table 1. Main parameters of HTR-10 reactor model

表1. HTR-10反应堆模型主要参数

物理量/单位	参数值
燃料球床体积/m^3	5
堆芯等效直径/cm	180
堆芯等效高度/cm	196.5
顶反射层厚度/cm	90
锥形反射层厚度/cm	37.5
底反射层厚度/cm	121.2
侧反射层石墨等效厚度/cm	77.8
反射层石墨密度/(g/cm^3)	1.76
反射层石墨杂质总硼当量/(μg/cm^3)	8.52
顶含硼碳砖厚度/cm	40
底含硼碳砖厚度/cm	30
底不含碳硼砖厚度/cm	70
侧含硼碳砖等效厚度/cm	22.2
碳砖密度/(g/cm^3)	1.59

HTR-10初装堆物理计算分区模型如图1所示。其堆芯由球形燃料元件和石墨及碳砖构件组成，堆芯燃料区等效直径为180 cm，等效高度为196.5 cm [4]。反应堆进行装料时，混合球从顶部装料管中进入堆芯燃料区，并堆积形成一个圆锥形的锥顶。在堆芯下部有一个圆锥形的堆料区和圆柱形的卸料管，底部圆锥形的堆料区使得离中心较远的燃料球也能从上至下顺利流动至卸料管。在堆芯燃料区外侧的石墨侧反射层中布置了10个控制棒孔道、3个样品孔道、7个硼吸收球孔道和20个冷氦气流通孔道[5]。在启堆前下部堆料区和卸料管中装满石墨球，堆芯初始混合装载的燃料球和石墨球比例为0.57:0.43 [6]。HTR-10反应堆模型主要结构参数如表1所示。

2.2. 燃料球结构

HTR-10反应堆使用氦气作为冷却剂，石墨作为慢化剂。其堆芯燃料区采用“全陶瓷型”球形燃料元件，该元件整体直径为6 cm。其结构由内外的石墨层构成：外部是厚达1 cm的石墨壳层，起机械保护作用；内部是直径5 cm的燃料区，由包覆燃料颗粒均匀弥散于石墨基体中构成。每个这样的燃料区平均含有约8000个包覆颗粒燃料核心[7]。

包覆颗粒燃料核心是一种直径0.5 mm的微型球体，其内部为UO₂颗粒。为约束裂变产物和承受气体裂变产物产生的内压力，UO₂颗粒外部设计了四层包覆结构，依次为疏松热解碳层、内层高密度热解碳层、热解碳化硅层和外层高密度热解碳层。单个燃料球的重金属装载量为5 g，其中²³⁵U的富集度为17%。堆芯中混合装载的石墨球密度为1.84 g/cm³，直径与燃料球相同，均为6 cm。

3. 球床随机堆积模拟

3.1. 基于离散元法的球床随机堆积模型

为模拟球床堆中混合球位置的随机分布特性，本文采用基于离散元方法(DEM)的EDEM仿真软件[8]，模拟混合球在重力作用下的随机下落过程，并据此确定其最终空间位置。EDEM软件基于牛顿运动定律，能够准确模拟燃料球的运动轨迹及球与球、球与壁面之间的碰撞行为。

在建立球床随机堆积模型时，首先构建一个圆柱形区域作为燃料堆积区，其高度为196.5 cm，半径为90 cm，侧壁与底面设置为封闭结构，以承接下落的球形燃料元件。在该圆柱区域顶部中心位置设置一根半径为25 cm的圆柱形装料管，确保生成的小球在管内运动而不超出边界。在装料管上方设置一个半径为25 cm的圆形入口面，模拟实际装料过程中燃料球与石墨球经由该面进入装料管的行为。该圆形面在模型中定义为“混合球工厂”，负责生成燃料球与石墨球。球形元件在重力作用下下落，期间经历多次碰撞，最终随机分布于燃料堆积区，形成随机堆积结构。

图2展示了所建立的几何结构及模拟过程：下部浅灰色圆柱体为燃料堆积区，上部黄色圆柱管为装料管，顶部橙色圆形区域为小球生成面。其中，紫色小球代表燃料球，黄色小球代表石墨球。

Figure 2. Simulation process of DEM method for pebble bed

图2. DEM方法的球床模拟过程

在DEM方法中，小球被视为弹性体，其相互之间及与管壁接触时会发生挤压变形并产生重叠。然而，这种几何重叠在后续蒙卡建模时会导致计算中断。为避免该问题，将小球模型半径设置为3.001 cm。经验证，此设置可有效避免几何曲面重叠，且半径仅增大0.001 cm对球床堆芯填充率的影响可忽略不计。因此，在基于所获球心坐标建立燃料球与石墨球模型时，其半径仍取为3 cm。

HTR-10初装堆采用燃料球与石墨球均匀混合方案，二者比例为0.57:0.43。在模拟中，两种小球随机于装料管内生成，并在重力作用下自由下落。过程中考虑小球之间及小球与管壁之间的摩擦与碰撞行为，相关摩擦系数与弹性常数参照文献[9]设定。当所有小球生成、下落并达到静止状态后，球床随机堆积模拟结束，最终状态如图3所示。

Figure 3. Final state of the randomly packed pebble bed

图3. 球床随机堆积最终状态

经EDEM软件计算，所得随机堆积球床的填充率为0.61，与实验值高度吻合。最终，通过EDEM软件导出两类小球的球心坐标，从而获得随机堆积状态下所有球形元件的空间分布数据。

为了研究球床顶锥对于蒙卡模拟的影响，本文的随机堆积方式包含两种构型：其一为混合球随机下落，在堆顶形成约20˚的锥角；其二为采用混合球进行随机下落模拟后，截取具有同样数量混合球的等效圆柱高度的球床区域，使顶部锥角为0˚，如图4所示，图中黄色与紫色小球为保留的混合球，顶部橙色小球为删除部分。

Figure 4. Model of cropped mixed balls at top cone

图4. 截取顶锥混合球模型

将EDEM软件计算得到的混合球位置坐标转换为蒙卡程序栅元卡格式，用于构建蒙卡程序的输入文件。在蒙卡程序中，每个燃料球或石墨球均被定义为一个独立的栅元，并配以相应的曲面卡描述其几何边界，所有球形栅元以外的区域均以空气介质填充，从而真实完整地描述了整个随机堆积球床的几何结构。

3.2. 球床随机堆积模型计算结果分析

根据文献的实验结果，HTR-10初次临界时(k_eff = 1)测得混合球数为16,890个。在EDEM软件中采用同样的混合球数16,890个，同时保证燃料球：石墨球 = 0.57:0.43以及堆芯填充率为0.61。蒙卡程序采用国际通用的MCNP5，搭配两个数据库分别进行计算，其一选用ENDF/B-VII的70c数据库，其二选用传统的基于ENDF/B-VI的60c数据库，中子代数和每代粒子数设置为500和10,000，k_eff计算结果如表2所示。

Table 2. Calculation results of k_eff for HTR-10 randomly packed pebble bed at the first critical height based on DEM method

表2. 基于DEM方法模拟HTR-10初次临界球床随机堆积的k_eff计算结果

	球床随机堆积构型	混合球总数	keff	标准方差
ENDF/B-VII 70c	锥顶	16,890	0.99921	0.00050
	平顶	16,890	0.99761	0.00046
ENDF/B-VII 60c	锥顶	16,890	1.00115	0.00019
	平顶	16,890	0.99916	0.00045

从表2可以看出，本文基于DEM方法模拟HTR-10初次临界球床随机堆积的k_eff计算结果与实验结果非常接近，说明本文建立的球床随机堆积模型的计算精度较高。此外，还可以看出，球床是否有锥顶对k_eff影响较小，差别小于0.2%，数据库选用基于ENDF/B-VII的70c数据库或者基于ENDF/B-VI的60c数据库对结果的影响也较小，差别小于0.2%。球床随机堆积过程因每次燃料球和石墨球下落的轨迹和最后位置都是随机的，导致每次模拟计算得到的k_eff存在波动，波动范围约为[−0.2%, +0.2%]。后续采用蒙卡程序模拟球床随机堆积时，模型选用锥顶构型，数据库选用基于ENDF/B-VI的60c数据库。

4. 球床规则排布模拟

基于离散元法真实模拟球床随机堆积，模拟精度较高，但是模拟时间较长。对于HTR-10初次临界问题(混合球数量16,890)，基于离散元法的EDEM仿真软件运行时间约为3天，因此对于更大规模的球床堆，比如高温气冷堆核电站示范工程(HTR-PM)堆芯共含有约42万个燃料球，采用离散元法模拟消耗的时间会变得不可接受。此外，蒙卡程序MCNP对于栅元数量也有限制(小于10万个栅元)，因此在实际工程计算时，通常基于重复结构功能将球床的随机堆积近似为规则排布，这样计算效率也更高。

4.1. 基于燃料球和石墨球随机分配的球床规则排布模型

采用蒙卡程序重复结构功能对球床堆芯进行规则排布近似时，先将球形元件在每一层规则排布，再逐层向上堆叠形成整个球床。实际球床堆的填充率为0.61，若直接以半径3 cm的球体、边长6 cm的立方体或边长3 cm的六棱柱栅格进行建模，所得填充率均低于该实际值。为此，需构造一种特殊几何形状的元件球，以满足填充率等于0.61的要求[10]，其结构如图5所示。经计算，对应立方体栅格边长为5.7022 cm，燃料元件球等效半径为3.013 cm。

Figure 5. Sectional view of a spherical element with a special geometric shape in a cubic grid

图5. 特殊几何形状元件球在正方体栅格中的剖面图

球床堆在初装堆和过渡过程阶段，会采用一部分石墨球来参与控制堆芯反应性。从初装堆到平衡堆芯过程中，球床高度以及燃料球和石墨球比例会持续变化，这就给蒙卡程序的模拟带来困难。在过去的研究工作中，采用规则排布近似球床堆芯时，通常通过人工排布来控制燃料球和石墨球比例，非常繁琐耗时。本文基于MATLAB平台开发了燃料球和石墨球随机分配程序，用于保证燃料球和石墨球比例为设定值时，实现燃料球和石墨球在球床栅格阵列中的位置随机分配，极大提高了建模效率。燃料球和石墨球随机分配程序的原理如下：将球床划分为按层排列，层高为特殊几何元件球所对应的正方体栅格边长5.7022 cm；将每层进一步划分为5.7022 cm × 5.7022 cm × 5.7022 cm的正方体阵列，按燃料球和石墨球的比例对每层的正方体阵列进行随机填充；将所有层进行堆叠得到三维正方体阵列结构，再由圆柱堆芯边界进行切割，位于圆柱以外的部分被舍弃，进而得到燃料球和石墨球在指定比例下随机分配的规则排布球床堆芯。

4.2. 球床规则排布模型计算结果分析

球床规则排布是按层级叠加，HTR-10球床初次临界高度在规则排布的两个层级之间，需要进行插值计算。HTR-10球床初次临界混合球总数为16,890个，相邻的两个层级高度为119.7462 cm和125.4484 cm，对应的混合球数量为16,435个和17,218个。此外，为了研究燃料球与石墨球比例变化带来的影响，人为设定了两种比例，即全部为燃料球、燃料球与石墨球数量相等。k_eff计算结果如表3所示。

Table 3. Calculation results of k_eff for HTR-10 pebble bed arranged according to the rules

表3. HTR-10球床规则排布的k_eff计算结果

燃料球:石墨球	混合球总数	规则排布keff	标准方差
全部为燃料球	16,435	1.06137	0.00035
全部为燃料球	17,218	1.07675	0.00035
57:43	16,435	0.98943	0.00036
57:43	17,218	1.00713	0.00033
1:1	16,435	0.96706	0.00034
1:1	17,218	0.98517	0.00034

球床规则排布模型在每一层级燃料球和石墨球的排布位置是随机的，导致每次模拟计算得到的k_eff存在波动，波动范围约为[−0.05%, +0.05%]。对表3中的球床规则排布k_eff进行插值计算，得到HTR-10球床初次临界高度的k_eff计算结果如表4所示。

Table 4. Calculation results of k_eff at the first critical height of HTR-10 pebble bed

表4. HTR-10球床初次临界高度的k_eff计算结果

燃料球：石墨球	混合球总数	规则排布keff	随机堆积keff	相对偏差
全部为燃料球	16,890	1.07031	1.07180	0.00139
57:43	16,890	0.99972	1.00115	0.00143
1:1	16,890	0.97758	0.98093	0.00342

从表4可以看出，当燃料球：石墨球 = 0.57:0.43时，本文基于“蒙卡重复结构规则排布 + 燃料球和石墨球随机分配”方法模拟HTR-10初次临界球床随机堆积的k_eff计算结果与实验结果(k_eff = 1)非常接近，说明本文建立的球床规则排布模型的计算精度较高。此外，规则排布模型的k_eff计算结果略小于随机堆积模型的，两者总体较为接近，当燃料球数量增多时，两者的差别进一步减小。

为了研究球床高度变化带来的影响，基于57%燃料球数量比例，逐渐提升球床高度至满功率额定高度。

Table 5. Calculation results of k_eff at different heights of HTR-10 pebble bed

表5. HTR-10不同球床高度的k_eff计算结果

混合球总数	规则排布高度(cm)	规则排布keff	随机堆积keff	相对偏差
12,522	91.23520	0.87478	0.87907	0.00488
14,087	102.6396	0.92540	0.92948	0.00439
15,653	114.0440	0.96916	0.97246	0.00339
17,218	125.4484	1.00713	1.00775	0.00062
18,783	136.8528	1.04043	1.03992	0.00049
20,348	148.2572	1.06868	1.06820	0.00045
21,914	159.6616	1.09350	1.09302	0.00044
23,479	171.0660	1.11671	1.11625	0.00041
25,044	182.4704	1.13713	1.13667	0.00040

从表5可以看出，随着球床高度的升高，规则排布模型的k_eff计算结果与随机堆积模型的差别逐渐减小，在初次临界高度附近两者的误差范围为[−0.00392, 0.00104]，在满功率额定高度附近两者的误差范围为[−0.00209, 0.00291]。

当中子代数和每代粒子数设置为500和10,000时，MCNP程序采用球床随机堆积模型的计算时间约为24小时，采用球床规则排布模型的计算时间约为3小时，意味着规则排布模型具有较高的计算效率，此外，研究发现当混合球数量进一步增加后，随机堆积模型的计算时间将大幅增加。

5. 总结

本研究基于HTR-10堆芯结构，分别构建了球床随机堆积与规则排布两种计算模型。随机堆积模型的构建过程如下：首先建立反应堆堆芯几何模型，随后采用离散元方法(DEM)模拟获得球床中混合球的随机空间坐标，最后通过蒙卡程序计算得到k_eff。规则排布模型则借助蒙卡程序的重复结构功能构造球床栅格阵列，并基于MATLAB平台开发了燃料球和石墨球随机分配程序，用于保证燃料球和石墨球比例为设定值时，实现燃料球和石墨球在球床栅格阵列中的位置随机分配。

对于HTR-10初次临界状态(室温，空气气氛，混合球装载数量16,890个，k_eff = 1)，两种球床模型的k_eff计算结果都与实验结果较为接近(偏差小于0.4%)，说明本文建立的球床随机堆积模型和规则排布模型都具有较高的计算精度。当燃料球占比增加和球床高度增加时，两种球床模型的k_eff计算结果的偏差进一步减小，更为接近。此外，相较于球床随机堆积模型，规则排布模型尽管在几何结构上做了一些近似，但是也实现了燃料球和石墨球比例为设定值时在球床栅格阵列中的位置随机分配，并且更易于满足蒙卡程序对于球形元件模拟规模和效率的要求，因而更推荐其作为工程计算模型。

基金项目

中央高校基本科研业务费专项(2024MS050)，国家磁约束核聚变能发展研究专项(2019YFE03110003)。

NOTES

^*第一作者。

^#通讯作者。