冀锐敏，严 睿，康旭忠，杨 璞，于世和，邹 杨

(中国科学院 上海应用物理研究所，上海 201800)

中国科学院于2011年启动了战略性先导专项“未来先进核能系统——钍基熔盐堆(TMSR)核能系统”[1]，目标包括建造1座固态燃料钍基熔盐堆(TMSR-SF1)。TMSR-SF1是热功率为10 MW的氟盐冷却高温堆(FHR)，是21世纪初美国新提出的第4代先进反应堆，结合了高温气冷堆和熔盐堆的优点[2]。TMSR-SF1采用直径6 cm的球形燃料和2LiF-BeF2冷却剂。由于密度小于熔盐，燃料球从堆芯底部投入，靠浮力上升后逐步构筑堆芯。已开展的相关堆积实验[3-6]表明，TMSR-SF1堆积密度低于HTR-10，且不确定性较大。

TMSR-SF1堆芯结构是物理和热工设计的重要基础，也是其不确定性的主要来源。为解决该问题，本文采用蒙特卡罗程序(MCNP)[7]基于TMSR-SF1设计完成高保真物理分析。TMSR-SF1中子物理计算面临很大挑战：上万颗燃料球在堆芯内随机排布以及每颗燃料球的燃料区域内均有上万个包覆颗粒随机排布。蒙特卡罗程序通常采用3种方法完成球床堆[8-11]随机结构的模拟：1) 建立简立方(SC)、体心立方(BCC)等规则结构完成随机结构的近似(方法1)。HTR-10首次临界计算中曾采用了该方法，并获得了与实验较为吻合的结果[12-13]。2) 不提前精确定位，在计算中对其随机定位[14-17]，MCNP尚不具备该功能(方法2)。3) 用数值方法模拟随机结构(方法3)。通常采用离散单元法(DEM)和软球模型进行受力分析，获取各燃料球的精确位置[11,17]。此方法有3点不足：1) 计算模型导致球发生形变；2) 输入参数中有些依据经验输入，如恢复系数、摩擦系数等；3) 所需的计算资源多。综合以上因素，本文使用方法1和方法3模拟TMSR-SF1堆芯结构。方法3选择工业领域内矿团堆积中使用的RSA方法[18-19]，其基本思想是逐个加入新球，直至完成给定区域的填充。新球位置的计算参考Halo方法[20]，根据现有球床计算出所有可能的新球的位置。随后对上述位置进行判断，排除超出边界、发生重叠和不能获得稳定支撑的球体位置[21]，以获取稳定的堆积。为加速计算，对计算区域进行网格化，对各球体位置进行像素化。基于上述两种方法建模，分析不同燃料球床以及燃料球描述模型对TMSR-SF1主要中子物理参数的影响，如keff、堆芯能谱、控制棒价值和温度系数等。

1 计算方法

TMSR-SF1堆芯结构非常复杂，故临界计算选用MCNP(版本6-6.1.1)，截面数据库采用ENDF/B-Ⅵ.0。

1.1 随机球床

本文分别采用随机模型和规则模型来模拟随机球床。

1) 随机模型

计算机模拟已成为随机球床的主要研究方法[14-16]，分为物理和几何两类。物理方法如分子动力学和DEM方法重在球体的堆积过程，而几何方法重在堆积结果。几何方法主要分为序列添加和集合重排两种，具体的实现方法各不相同。本文采用随机序列添加(RSA)方法，其基本思想为每次迭代过程中生成1个新球到现有系统中。新增球需与已填充球或壁面有3个接触点。新球的产生采用Halo方法[20]：1) 与a球相切的所有球心可形成1个球面，以a球的球心为原点，2r(r为球半径)为半径，简称a球Halo(图1a)。2) a球和b球的Halo交集，即与a球和b球均相切的球心集合，可形成1个圆(图1b)，标识为圆(cc，rc)。向量cc表示圆心位置，rc为圆半径，可由式(1)计算得到。向量a和b为a球和b球的球心位置。3) 与球a、b

a——a球Halo；b——a球和b球Halo的交集；c——c球Halo与圆(cc,rc)的交集；d——球心集合{P1,P2}图1 Halo方法示意图Fig.1 Halo algorithm scheme

和c均相切的球心可通过投影计算。c球Halo在圆(cc，rc)所属平面的截面同样为1个圆(cP，rP)，如图1c所示。圆心位置cP和圆半径rP可由式(2)计算得到。向量c为c球球心位置，向量n为a球到b球方向的单位向量，λ为向量ccc在向量n上的投影。4) 式(3)中的判据L可判断圆(cc，rc)和圆(cP，rP)是否存在交集：若L>0，存在2个交点；若L=0，存在1个交点；若L<0，没有交点。5) 若存在交集，交点即为与a、b和c球同时相切的新球球心。球心位置P1,2可通过式(4)、(5)计算得到。cm为两个新球的中点，α和h为简化公式引入的系数。

(1)

λ=(c-cc)n,n=(a-b)/|a-b|

(2)

L=rP+rc-|cc-cP|

(3)

cm=(1-α)cc+αcP

(4)

P1,2=cm±h(n×cccm)/|n×cccm|

(5)

计算出的新球还需进行有效性判断：不得超出堆芯活性区；不能与已有燃料球发生重叠；落位应稳定[20-21]。在上述过程中需多次检索给定球的邻球。为加速检索，初始将整个堆芯区域进行网格化，预分为若干立方体区域，立方体的边长取燃料球的半径r。根据燃料球位置产生与之对应的像素。假设堆芯原点位置为O=[Ox，Oy，Oz]，位置为P=[Px，Py，Pz]的燃料球可定位为(i，j，k)，如式(6)所示。

(6)

图2 球床计算流程Fig.2 Flowchart for pebble bed generation

随机球床的计算流程如图2所示。简述如下：输入燃料球半径和堆芯活性区信息，将堆芯网格化，建立已填充燃料球序列和待填充燃料球序列，用于存储燃料球编号和位置信息；随机产生放置于堆芯活性区顶部的1个球；利用Halo方法结合已有燃料球计算所有可能的新球位置；对其进行有效性判断；将有效的燃料球放置于待填充球序列；依据填充的密集程度，对待填充燃料球进行排序；取出排名第一的球至已填充序列；根据新填充的燃料球重复上述过程，直至待填充序列为空集，即堆芯活性区已填满燃料球。

2) 规则模型

HTR-10在首次临界计算中曾采用一些规则的模型来模拟随机堆积的燃料球床，如BCC模型、FCC模型、CHPOP模型等[12-13]，并获得了与实验较为一致的结果。本文使用BCC模型进行等效处理，球床模型细节示于图3。

a——周期性规则填充的栅格；b——燃料球堆积水平方向；c——燃料球堆积竖直方向图3 燃料球堆积的规则等效模型Fig.3 Equivalent regular model of fuel pebble packing

1.2 燃料球

图4 燃料球模型与URAN扰动原理示意图Fig.4 Schematic of fuel pebble and URAN vibration principle

受限于MCNP程序中对栅元、曲面等编号数量的限制，并未对燃料球中的包覆颗粒进行随机位置的建模，而将其进行规则的建模，然后通过MCNP中的URAN卡对其位置进行随机扰动[7]。URAN卡实现扰动的基本原理图示于图4。详细步骤如下：1) 对燃料球中包覆颗粒详细建模；2) 将包覆颗粒填充至燃料填充区域中划分好的简单立方格中，简单立方格的尺寸由燃料球中包覆颗粒的堆积密度确定；3) 设定URAN卡，指定扰动的栅元编号以及扰动范围；4) 在计算过程中，待抽样到该栅元时，抽样执行随机扰动，将包覆颗粒随机移动(式(7))。其中，[x，y，z]为扰动后的位置；[x0，y0，z0]为扰动前的位置；ξ1、ξ2和ξ3为0～1的随机数；δx、δy和δz为指定的扰动范围。值得指出的是，该扰动受扰动范围的限定，通常限制于划分好的简单立方网格内。即实现了燃料球内包覆颗粒有限的随机分布，与真实的分布仍有一定的差异。单燃料球中包覆颗粒分析表明，相比于完全的随机，这种简单的处理方式会导致边界处包覆颗粒的切割，引起无限增殖因数低估[10]。

(7)

2 TMSR-SF1堆芯

TMSR-SF1使用全陶瓷型包覆颗粒制成的球形燃料元件、氟锂铍熔盐(FLiBe)作为冷却剂、石墨作为慢化剂和堆芯的结构材料，热功率为10 MW，其堆芯结构示意图如图5a、b所示，主要由堆芯活性区和石墨反射层组成。堆芯活性区由反射层的石墨构件围成，形成上下圆台和中间圆柱区域。圆柱直径为135 cm，高为180 cm；上下圆台最小直径为30 cm；圆台斜面与水平面的夹角为30°；反射层外围形状为圆柱体，直径为285 cm，高为306 cm，反射层外的堆芯围筒厚度为2 cm。堆芯活性区内燃料球随机堆积，燃料球之间的空隙形成不规则的流道，供冷却剂自下向上流动，带走裂变产生的热量。靠近活性区的石墨侧反射层中布置16个孔道，其中间隔布置的12根控制棒作为反应性控制系统，排空熔盐作为第2套停堆系统。燃料球由中心燃料填充区域和外部石墨壳组成，如图5c所示。TMSR-SF1设计中所用燃料球的基本参数列于表1。燃料填充区域内由包覆燃料颗粒(TRISO)和石墨基体混合而成。TRISO由UO2核心和4层包覆层组成。

a——堆芯三维图；b——堆芯俯视图；c——燃料球示意图图5 TMSR-SF1堆芯结构示意图Fig.5 Schematic diagram of TMSR-SF1

表1 TMSR-SF1球形燃料元件参数Table 1 TMSR-SF1 fuel element characteristic parameter

3 结果和讨论

3.1 球床结构

球床堆积模拟中的一个关键问题是球体堆积结构的描述。除堆芯平均堆积密度外，还分别使用投影以及轴向和径向堆积密度分布[22]两种方法给出两种不同建模方法所获球床的结构细节，如图6所示。图6a、b为RSA方法获得的随机球床在底面和侧面的投影，图6c、d为规则堆积方法获得的规则球床在底面和侧面的投影，图6e、f为与图6a、b对应的径向和轴向堆积密度分布，图6g、h为与图6c、d对应的径向和轴向堆积密度分布。可见，随机球床在靠近活性区边界处的区域以及球床顶部区域的堆积密度较高，轴向堆积密度随球床高度的降低而有所降低；规则球床在靠近活性区边界处的区域堆积密度较低，轴向堆积密度则为规律的锯齿状。对比可见，尽管两组球床模型的堆积密度均接近60%，但球床的细节结构有很大不同。

3.2 中子物理参数

使用3种建模方法计算满装载时的冷态反应性：1) RSA方法模拟随机堆积球床，依据各燃料球位置利用MCNP实现详细建模(方法1)；2) 等效规则球床，MCNP直接建模，通过不同的堆积密度拟合出给定堆积密度的情况(方法2)；3) 规则球床等效随机球床，在相同堆积密度的前提下，BCC结构填充原点有差异，以统计活性区边界对燃料球的不同切割效应(方法3)。为计算包覆颗粒随机分布的影响，在上述3种输入的前提下，通过MCNP程序URAN卡实现所有燃料球中包覆颗粒的随机移动，分别计算球随机分布以及包覆颗粒随机分布对计算结果的影响。

1)keff

MCNP计算中，每代粒子数100 000，循环总数1 000代，有效循环500代，计算得到有效增殖因数keff的标准偏差小于0.000 10。上述3种方法的计算结果列于表2。由表2可见，方法1的计算结果最大，方法2、3的计算结果非常接近。平均堆积密度为59.58%时，随机堆积模型结果较规则堆积模型结果约大500 pcm，接近0.5%。3种计算方法中由包覆颗粒随机分布导致的有效增殖因数差异分别为+31、+26、-21 pcm。由此可见，包覆颗粒的随机分布对满装载冷态有效增殖因数的影响较小，相对于燃料球随机分布其影响完全可忽略。

图6 球床堆积密度分布对比Fig.6 Comparison of packing density distribution of pebble bed

表2 TMSR-SF1满装载冷态有效增殖因数计算结果Table 2 keff for TMSR-SF1 fully-loaded cold state

2) 中子能谱

对于随机球床，距活性区边界约3个燃料球直径范围内的堆积密度出现较大振荡，因此将堆芯活性区(R<67.5 cm，R为模型径向距离)、活性区中心区域(R<50 cm)、活性区边界区域(50 cm

3) 控制棒价值和温度系数

为进一步量化不同模型的影响，基于第10组算例和前2种方法计算控制棒总价值和温度系数，结果列于表3。由表3可见，相对于随机模型，规则模型会高估控制棒价值约5.78%，而温度系数的结果较为一致。控制棒价值高估与堆芯有效增殖因数低估的原因一致，燃料球排布的简化增加了堆芯中子的泄漏。

图7 堆芯各区域中子通量密度分布对比Fig.7 Comparison of neutron flux density in different regions

表3 TMSR-SF1满装载冷态控制棒价值和温度系数Table 3 Control rod worth and temperature coefficient for TMSR-SF1 fully-loaded cold state

4 结论

基于TMSR-SF1设计模型，利用RSA方法(方法1)实现了随机球床的模拟，与BCC结构的规则等效球床进行了球床堆积(方法2)情况对比；基于MCNP程序中URAN卡完成了燃料球内包覆颗粒随机分布的模拟；针对上述计算模型进行了中子物理关键参数分析，得到如下主要结论。

1) 尽管随机球床和规则球床的平均堆积密度均接近60%，但球床结构的细节有很大不同。随机球床在靠近活性区边界约3个燃料球直径区域内形成较为规则的分布，径向堆积密度在该区域内有一定的振荡，在堆芯区域则较为平稳。随机球床的轴向堆积密度从顶部到底部逐步降低。规则球床的径向和轴向堆积密度均呈现锯齿状分布。

2) 由于球床结构的不同，规则球床模型会造成更多的中子泄漏，导致堆芯有效增殖因数的低估和布置于反射层的控制棒价值的高估，但对温度系数的影响可忽略。

3) 燃料球中包覆颗粒随机分布的影响远小于球床随机分布的影响。