熔盐冷却球床实验堆堆内热源分布计算与分析

杨 璞，严 睿，邹 杨，于世和，朱贵凤，周 波, 马玉雯,*

(1.中国科学院 上海应用物理研究所，上海 201800；2.中国科学院大学，北京 100049)

TMSR-SF1是热功率为10 MW的熔盐冷却球床实验堆，燃料元件采用球形燃料，冷却剂采用2LiF-BeF2(7Li丰度>99.99%)，相比于传统水堆，TMSR-SF1在堆芯结构、冷却剂、燃料组件等方面有很大的不同，而准确计算堆内释热率分布对于这种新型反应堆的热工水力设计、瞬态分析、结构力学设计等都有重要意义。本文使用蒙特卡罗计算程序MCNP(版本5-1.51)对TMSR-SF1寿期初(BOL)及寿期末(EOL)堆内部件(包括燃料球、冷却剂、控制棒、石墨反射层及合金结构材料等)的释热率进行详细计算。MCNP程序可通过F6及F7能量计数卡来统计裂变碎片、中子、瞬发γ和俘获γ的能量沉积，但该程序无法直接计算缓发β及缓发γ的能量沉积，因此一般通过统一乘1个系数来计入这两者对堆内部件释热率的贡献，但这种方法忽略了缓发β与缓发γ在堆内能量沉积分布的不同。本文通过使用光子产生偏倚卡(pikmt)，经过3次MCNP输运计算，分别统计裂变碎片、中子、瞬发γ、俘获γ、缓发β及缓发γ在TMSR-SF1堆内不同部件内的能量沉积，并最终计算出部件内的释热率。

1 释热率计算方法

TMSR-SF1堆内释热率分布计算需要分别进行3次MCNP输运计算，即：1) 栅元体积计算；2) 使用F6:n、F6:p和F7:n统计卡的常规MCNP计算；3) 使用光子产生偏倚卡及F6:p统计卡的MCNP计算。计算结果处理需要两步，首先需通过常规MCNP计算F6:n、F6:p和F7:n卡给出的统计结果来求得平均单次裂变放出的可利用能量(Q)，然后根据Q值来计算堆内不同部件的总释热率及体积释热率。

1.1 Q值计算方法

根据反应堆总热功率和平均单次裂变放出的可利用能量(即Q)可计算反应堆单位时间的裂变次数，因此对于反应堆释热率计算，Q值的计算是非常重要的一步。反应堆Q值由以下几部分组成[9]：裂变碎片的能量贡献(Qfp)、中子的能量贡献(Qn)、瞬发γ的能量贡献(Qγp)、俘获γ的能量贡献(Qγc)、裂变产物衰变放出的缓发β粒子能量贡献(Qβ)和缓发γ的能量贡献(Qγd)。MCNP程序中的F6和F7卡可用来统计裂变碎片、中子、缓发和俘获γ的能量。其中F6:n用于统计裂变碎片和中子的能量，F6:p用于统计瞬发和俘获γ的能量，F7用于统计裂变碎片、中子、瞬发γ的能量且与F6卡统计原理不同，F7卡假设所有能量就地沉积。因此，可用下式计算Q值[10]：

(1)

(2)

(3)

由于MCNP本身无法直接计算缓发γ的能量贡献Qγd和缓发β的能量贡献Qβ，因此这两部分直接使用ENDF/B-Ⅶ.0数据库中235U的缓发γ和缓发β的能量数据，分别取6.33 MeV和6.50 MeV[11]。一般反应堆内的缓发γ主要来源于235U裂变产物的衰变，因此使用235U的缓发γ的能量作为整个反应堆的Qγd值是一个合理的假设[10,12-13]。β粒子穿透力很低，可认为β粒子的全部能量就地沉积，因此堆内缓发β的能量沉积近似和反应堆的裂变率呈正比，对于使用U燃料的反应堆，这部分能量可直接使用235U的缓发β的能量数据[10]。

1.2 释热率计算方法

反应堆堆内部件的释热率主要来源于裂变碎片、中子、缓发β、瞬发γ、俘获γ和缓发γ的能量沉积。MCNP通过栅元来对反应堆内的实际部件进行建模，某个栅元i内裂变碎片、中子、瞬发和俘获γ的能量沉积可使用下式计算：

(4)

(5)

采用SPSS15.0统计学软件,计量资料采用()的形式表示,用t检验;计数资料采用x2检验,P<0.05为差异具有统计学意义。

(6)

栅元i的释热率为上述所有粒子在栅元i内能量沉积的总和(式(7))。栅元i的释热率除以其栅元体积便得到栅元i的体积释热率(式(8))。反应堆的总释热率(即总功率)为所有栅元释热率的总和(式(9))。

(7)

(8)

(9)

图1 释热率计算流程Fig.1 Computational step for heat rate

1.3 计算流程

为简化上述计算过程中数据传输和消除容易出错的手动数据传输过程，本文使用Python语言编写用于输入输出处理的程序。图1为整个计算流程。计算过程包括3次MCNP运行：第1次运行采用固定源模式并配合材料无效卡(void)，用于计算栅元的体积；第2次运行则采用临界模式(使用kcode卡)配合F6:n、F6:p和F7:n计数卡，用来计算裂变碎片、中子、缓发β、瞬发γ和俘获γ的能量沉积；第3次运行在临界模式使用pikmt卡且配合F6:p计数卡，用以计算缓发γ的能量沉积。最后，采用自编程序自动读取输出文件中的统计结果，通过数据处理得到Q值及各栅元的释热率及平均体积释热率，并将这些计算结果整理和输出到excel文件中。

2 TMSR-SF1堆芯描述

TMSR-SF1堆芯主要由堆芯活性区、石墨反射层和哈氏合金围筒组成。堆芯活性区为燃料球和石墨球的随机堆积区，燃料球及石墨球间空隙形成熔盐随机不规则流道，供其自下往上流动带走热量，活性区由中心圆柱体和上下圆台组成，其中圆柱体直径为135 cm、高度为180 cm，上下圆台直径为30 cm、高度为30 cm。活性区外为石墨反射层，其内布置16个控制棒通道、2个燃料球注入通道、1个熔盐装卸通道、1个中子源通道、2个物理启动用测量通道、7个中子通量与能谱测量通道、6个温度测量通道、3个随堆监督通道及2个备用通道。每个通道内还包含套管，其中温度测量通道、中子通量与能谱测量通道、燃料球注入通道及中子源通道内为哈氏合金套管，其他通道内则为碳-碳复合材料(C-C)套管。反射层外为堆芯围筒，其厚度为3 cm。整个堆芯直径为292.6 cm、高度为308.4 cm。堆芯具体结构如图2所示。

图2 TMSR-SF1堆芯结构Fig.2 Schematic of TMSR-SF1 core

TMSR-SF1活性区装载的燃料球与石墨球直径为6 cm，其中燃料球具体结构示于图3，主要参数列于表1。

图3 燃料球及TRISO包覆颗粒结构Fig.3 Fuel pebble and TRISO coated particle structure

参数数值直径,cm6燃料区域直径,cm5石墨基体和石墨球壳密度,g/cm31.73235U富集度,%17.0单个燃料球U装量(纯铀),g7.00燃料颗粒直径,mm0.5UO2核芯密度,g/cm310.4包覆层材料(从内到外)PrC/PyC/SiC/PyC包覆层厚度,μm95/40/35/40包覆层密度,g/cm31.1/1.90/3.18/1.90

3 TMSR-SF1的MCNP模型

考虑到TRISO包覆颗粒在燃料球内的填充比一定的情况下，其随机分布与规则分布计算的燃料球无限介质增殖因数并无太大差别[15]，因此为简化建模过程，本文采用简单立方建模方式来模拟TRISO颗粒在燃料球内的排布，为保证单颗燃料球铀装量为7.00 g，简单立方的边长被设为0.178 cm。TRISO及燃料球模型结构如图4所示。

TMSR-SF1堆内燃料球采用随机堆积方式，燃料球在堆芯内的占空比为60%，但在MCNP建模过程中由于随机堆积模型存在较大的困难和不确定性[7]，因此，本文采用图5所示的体心立方排布方式(占空比60%)对堆芯内的燃料球进行等效建模，体心立方边长为7.224 cm。

图4 TRISO颗粒及燃料球模型Fig.4 Geometry model of TRISO and fuel pebble

图5 随机球床堆的燃料球堆积等效Fig.5 Equivalent of random packing of fuel pebble

图6为最终的MCNP模型。模型主要包括堆芯活性区、石墨反射层、堆内通道、通道套管、控制棒及堆芯哈氏合金围筒等结构。

图6 TMSR-SF1的MCNP模型Fig.6 MCNP model of TMSR-SF1

4 计算结果及分析

表2为TMSR-SF1在寿期初和寿期末时的Q值计算结果。从表2可见，寿期末的Q值相较于寿期初略有增加，这主要是由于寿期末有一定量的239Pu产生，而239Pu的单次裂变能大于235U的。

表2 TMSR-SF1 Q值计算结果Table 2 Fission Q value of TMSR-SF1 model

表3为TMSR-SF1在寿期初和寿期末时堆内不同部件的释热率计算结果，其中合金构件包括堆芯围桶及通道合金套管。由表3可见，TMSR-SF1的堆内释热主要集中在燃料球内，其释热率占堆内总释热率的94%以上，熔盐和反射层由于在堆内体积较大，因此释热率也较高，分别占到总释热率的3%和1%以上，其他堆内部件占总释热率的比例都小于1%。与寿期初相比，寿期末燃料球释热率略有减小，这主要归因于燃料球的释热率主要来源于裂变碎片和裂变中子的能量沉积，且寿期末单位时间的裂变次数小于寿期初的。此外，控制棒与石墨球的释热率也有所减少，前者主要是寿期末控制棒移动到堆顶，导致能入射到控制棒的中子及γ减少造成的。而石墨球的释热率减少，主要原因是堆芯活性区中子和γ轴向通量密度分布的峰值上移，而石墨球位于堆芯活性底部，所以入射到石墨球的中子及γ也有所减少。而由于寿期末控制棒上移，堆内产生的裂变中子被控制棒吸收的概率减小，中子及中子被吸收后产生的俘获γ在堆芯活性区外的部件内沉积的能量增加，导致反射层、冷却剂熔盐、孔道C-C套管及堆内合金构件的释热率增加。

表3 TMSR-SF1堆内释热率分布计算结果Table 3 Heat rate distribution results of TMSR-SF1

图7为寿期初和寿期末堆内主要结构件(非燃料球)的释热来源。由图7可见，结构件的释热率主要来源于γ的能量沉积，其中包含合金材料的控制棒及合金构件的释热主要是来源于俘获γ的能量沉积，这是因为合金材料会俘获中子发生(n，γ)反应并放出大量的俘获γ，而俘获γ的大部分能量会就地在合金材料内沉积。

表4为不同部件内的平均和最大体积释热率。图8为寿期初和寿期末全堆体积释热率的精细分布。可见，燃料球的体积释热率远大于其他堆内部件，寿期初的最大体积释热率为22.14 W/cm3，寿期末的最大体积释热率为20.26 W/cm3。反射层靠近活性区的部分体积释热率较大，最大值达到约0.7 W/cm3，而远离活性区的部分体积释热率仅约10-5W/cm3量级。此外，相比于寿期初，寿期末反射层周向的体积释热率分布更加均匀。控制棒的最大体积释热率位于控制棒底端的配重合金棒头，寿期初和寿期末的最大体积释热率分别可达1.408 W/cm3和0.9502 W/cm3。堆内合金构件的最大体积释热率位于中子源通道的合金套管，寿期初为0.792 1 W/cm3，寿期末增加到0.822 4 W/cm3。

图7 TMSR-SF1堆内结构件释热来源Fig.7 Heat source of component in TMSR-SF1

堆内构件平均体积释热率/(W·cm-3)最大体积释热率/(W·cm-3)BOLEOLBOLEOL燃料球6.3636.34922.14020.260石墨球0.0470.0380.2230.196反射层0.0090.0110.7130.751冷却剂熔盐0.1290.1281.5341.464孔道C-C套管0.0260.0280.1550.125控制棒0.0740.0541.4080.950合金构件0.0300.0400.7920.822

根据经验，相对误差小于5%的结果通常认为是可靠的[16]。表5为上述计算结果的相对误差。体积计算投入粒子数为5×1010，体积计算结果的相对误差不超过0.02%。能量沉积计算中单代投入粒子数为100 000，有效迭代次数为4 550，所有F6、F7计数卡结果的相对误差都小于1%。

图8 TMSR-SF1体积释热率分布Fig.8 Volumetric heat rate distribution of TMSR-SF1

堆内构件相对误差/%体积计算BOLF6:n、F6:p、F7:n卡统计使用pikmt卡的F6:p统计体积计算EOLF6:n、F6:p、F7:n卡统计使用pikmt卡的F6:p统计燃料球0～0.010.01～0.020.01～0.020～0.010.01～0.020.01～0.02石墨球0～0.020.01～0.020.02～0.270～0.020.01～0.020.01～0.14反射层0～0.010.02～0.820.04～0.890～0.010.02～0.780.04～0.81冷却剂熔盐0～0.010.03～0.540.01～0.650～0.010.03～0.480.01～0.62孔道C-C套管0～0.020.02～0.460.02～0.570～0.020.02～0.450.02～0.55控制棒0～0.020.04～0.800.05～0.910～0.020.04～0.850.05～0.99合金构件0～0.020.05～0.760.04～0.81 0～0.020.05～0.720.03～0.78

5 结论

通过使用光子产生偏倚卡，经过3次MCNP输运计算，对TMSR-SF1寿期初和寿期末堆内的释热率进行了详细计算研究，得到了堆内不同部件的总释热率、体积释热率分布和最大体积释热率，并得到以下结论。

1) 燃料球释热率占堆内总释热率的94%以上，熔盐和反射层由于在堆内体积较大，因此释热率也较高，分别占总释热率的3%和1%以上，其他堆内部件释热率占总释热率的比例都小于1%。

2) 与寿期初相比，寿期末燃料球、控制棒与石墨球的释热率有所减少，而其他构件包括反射层、冷却剂熔盐、孔道C-C套管及合金构件的释热率都有所增加。

3) 包含合金材料的控制棒及堆内合金构件的释热率主要来源是俘获γ的能量沉积。控制棒的最大体积释热率位于控制棒底端的配重合金棒头，合金构件的最大体积释热率位于中子源通道的合金套管。