张竞宇 , 李富 , 孙玉良
清华大学 核能与新能源技术研究院, 北京 100084

收稿日期：2012-07-11
基金项目：国家自然科学基金资助项目（11375099, 11605058）；国家科技重大专项资助项目（ZX06901）
作者简介：张竞宇 (1984—), 男, 博士研究生
通信作者：李富, 研究员, E-mail：lifu@mail.tsinghua.edu.cn

摘要：球床高温气冷堆在初装堆芯建立过程中，会出现堆芯燃料球比例、空隙填充物质和顶空腔高度等不断变化的情况。基于确定论方法的设计程序VSOP在进行初装堆芯物理计算时需要采用多种特殊处理。为了明确这些特殊处理引入的近似程度，该文选用通用MonteCarlo程序更精细地模拟了初装堆芯。对比结果表明：对于初装堆芯，VSOP程序采用等效球体模型进行热谱计算、通过流修正方法对堆芯空隙填充物质反应性价值进行计算、以真空气氛近似氦气气氛、采用分方向扩散系数的方法对大尺寸顶空腔进行计算所引入的误差都很小。因此，VSOP程序对初装堆芯的这些近似处理是合理可行的。
关键词：球床高温气冷堆初装堆芯VSOP程序Monte Carlo程序验证
Neutronics calculation methods for the first core of the pebble bed high temperature gas cooled reactor
ZHANG Jingyu, LI Fu, SUN Yuliang
Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China


Abstract: During the building of the first core of the pebble bed high-temperature gas-cooled reactor (PBR), the ratio between the fuel spheres and graphite spheres, the fill material in the gaps in the core and the cavity height above the pebble bed will change. Thus, the neutronics design of the first core using the deterministic VSOP code needs some special treatments. The approximations for the first are an evaluated using a more detailed model of the first core analyzed by a general Monte Carlo code. Comparisons of the predictions of the two codes show that the special treatments for the first core used with the VSOP code will not cause significant errors include using the equivalent sphere model to calculate the thermal spectrum, calculating the reactivity of the gap fill material using the streaming correction method, heating the helium atmosphere as a vacuum, and using anisotropic diffusion coefficients to model the large top cavity. Thus, these special treatments for the first core used by the VSOP code are acceptable.
Key words: pebble bed high-temperature gas-cooled reactorfirst coreVSOP codeMonte Carlo codeverification
2008年2月15日，国务院批准了高温气冷堆核电站示范工程重大专项总体实施方案，目标是建设世界上第一座具有第四代核能系统安全特性的20万千瓦级模块式高温气冷堆核电站 (high-temperature gas-cooled reactor-pebble-bed module, HTR-PM)^[1]示范工程。
HTR-PM属于球床模块式高温气冷堆，其物理设计主要采用德国Jülich研究中心开发的基于确定论方法的VSOP程序^[2]，具有燃料元件设计、截面加工、中子能谱计算、少群扩散计算、燃耗计算、燃料倒料模拟、热工水力学计算等稳态物理和热工设计能力。
为高效地完成工程设计，VSOP程序采用了很多近似处理，包括分离的共振积分计算、组件输运计算 (0维或1维) 和全堆扩散计算 (2维或3维)，包覆颗粒在燃料球中随机分布和球形元件在球床中随机分布的处理，非规则区域等效为规则区域，采用分方向扩散系数的方法计算顶空腔，侧反射层孔道与周围石墨的均匀化处理等。
VSOP程序广泛用于球床高温气冷堆的物理设计，并采用AVR (arbeitsgemeinschaft versuchs reactor)、THTR (thorium high-temperature reactor)、HTR-10(10 MW high-temperature gas-cooled reactor-test module)^[3]的运行数据进行过验证，证明这些近似处理是比较有效的。但是，球床高温气冷堆的初装堆芯情况更复杂，需要更多的近似处理，这些近似处理的效果还需要进一步验证。
根据当前的认识，Monte Carlo程序可以对复杂几何结构进行精细描述、对中子输运过程进行直接模拟，其截面库中所包含的核素也是比较完整的，可以作为很多程序和处理方法的验证基准。因此，本文选用通用Monte Carlo程序 (MCNP)^[4]，精细描述无燃耗的初装堆芯模型，并与VSOP程序的分析结果进行比较，可以对VSOP程序的处理方法的合理性进行进一步的验证。
1 初装堆芯的物理计算方法球床高温气冷堆的运行模式是：燃料球从堆芯活性区顶部加入，在裂变发热过程中逐步流动到活性区底部，再卸出堆外；卸出堆外的燃料球测量其燃耗，未达到规定燃耗的燃料球，要重新装入堆芯，达到规定燃耗的燃料球，存入乏燃料库；装料和卸料在反应堆运行时连续进行。因此，长期运行以后，堆芯活性区每处都是由不同燃耗水平的燃料球混合而成，平均燃耗从堆芯上部到下部逐渐增加，而且可以达到一个稳定不变的状态，即平衡堆芯。但是，在反应堆开始满功率运行时，出于安全性和经济性的考虑，堆芯活性区会装载零燃耗的新燃料球和部分石墨球，此时的状态称之为初装堆芯。HTR-PM的初装堆芯的一个方案是，由低富集度燃料球 (²³⁵U的质量富集度为4.1%) 和石墨球按照7：8的数量比例混合而成。
VSOP程序将堆芯活性区划分成若干物质区，每个物质区单独进行能谱计算。对于快中子能区和超热中子能区，中子平均自由程大于球形元件的尺寸，能谱计算比较简单，只需把整个物质区直接均匀化后求解多群中子慢化方程即可^[5]。对于热中子能区，中子平均自由程较小，不能忽略物质区的几何结构。VSOP程序热谱计算模块的几何描述采用一维球模型或柱模型^[6]，对于全部由燃料球组成的平衡堆芯，直接采用燃料球模型、碰撞几率法和白边界条件求解多群中子输运方程得到热谱。但是，对于由燃料球和石墨球按比例均匀混合而成的初装堆芯，需要把物质区内的燃料球、石墨球和球床空隙 (球与球之间的空隙) 等效成一维球体模型。具体做法是将燃料球置于等效球体模型的中心，外面包覆体积等效后的石墨球球壳和球床空隙球壳，随后采用碰撞几率法和白边界条件求解多群中子输运方程得到热谱，这种处理方法显著改变了石墨球和球床空隙的形状以及石墨球的数量，与真实情况差别较大。
中子能够穿过球床空隙离开堆芯，导致中子泄漏增大，称为“中子流效应”^[7]。VSOP程序将球床均匀化后，由于忽略了中子流效应，使得中子平均自由程被低估。在扩散计算中，为了考虑中子流效应造成的中子运动增强的影响，需要借助公式D_corr，g=C_str，gD_g来提高扩散系数，称为“流修正”，修正因子C_str，g的取值由真空状态下的实验数据结合理论推导得到。但是在初装堆芯建立过程中，需要经历物理临界、燃料装载、气氛切换、升温升压和功率提升阶段^[8]，堆芯空隙 (包括顶空腔、球床空隙和反射层空隙等) 会填充湿空气、氮气和氦气，对这种非真空状态需要研究流修正的必要性。与此同时，由于VSOP程序的多群常数库中缺少核素氦的截面数据，因此直接以真空气氛来近似反应堆正常运行时的氦气气氛，可能会引入误差。
球床高温气冷堆在堆芯球床和顶反射层之间存在着一个圆柱形的空腔，称之为“顶空腔”。空腔区通常不能直接使用扩散方法。一种解决手段是建立全堆输运模型，但是所需的计算量过大。另一种解决手段是在顶空腔内填充稀释的石墨，然后进行全堆扩散计算，石墨的填充浓度预先借助输运计算来标定。这种方法在为顶空腔引入扩散计算所需的扩散系数的同时，也引入了不应该存在的中子反应截面，因此计算精度不高，并且对于顶空腔高度变化的情况，需要重新计算石墨的填充浓度，操作过程繁琐。VSOP程序采用的解决手段是将顶空腔作为扩散区域处理，吸收截面、散射截面和裂变截面都设为0，通过优化计算产生沿径向的扩散系数D_r和沿轴向的扩散系数D_z，并代入扩散方程求解顶空腔的中子通量分布，称为“分方向扩散系数的方法”^[9]。该方法对于小尺寸顶空腔的计算结果是准确的 (KAHTER实验装置顶空腔高50 cm，有效增殖因子k_eff的偏差为-0.05%；PNP-3000实验堆顶空腔高100 cm，k_eff的偏差为-0.01%)。HTR-PM由物理临界过渡到初装堆芯的过程中，伴随着燃料装载，顶空腔高度会从947 cm逐渐减小到67 cm，对于这种顶空腔高度大幅变化的情况，需要研究分方向扩散系数的方法是否仍然适用。
2 MCNP对球床堆芯的模拟2.1 用重复结构近似随机分布HTR-PM初装堆芯含有约42万个球形元件，燃料球体积份额为7/15，剩下的部分是石墨球，每个燃料球又含有上万个包覆颗粒。燃料球中包覆颗粒的分布和球床中球形元件的分布均具有一定的随机性，对于这种数量庞大的随机分布，MCNP不能精确描述，只能近似为重复结构进行描述。具体方法是：确定一个基本的栅格；用该栅格按照一定的方式填充给定的体积。MCNP中用输入卡LATTICE定义重复结构，LAT=l意味着栅格是一个六面体，LAT=2意味着栅格是一个六棱柱，随后将栅格在目标区域内进行无限阵列填充^[10]。
2.2 球床模型大量实验表明随机堆积的球床的填充率为0.61，VSOP程序中就采用此值。用MCNP模拟球床时，最基本的要求是满足与随机分布具有相同的球床填充率，此外还需要满足燃料球的体积份额为设计值，以保证堆芯的石墨和重金属含量与实际情况一致。满足以上要求的球床模型有多种，本文仅建立球床的真实模型和等效模型。
2.2.1 真实模型燃料球由半径为2.5 cm的球形燃料区和内径为5 cm的石墨球壳 (壳厚度为0.5 cm) 构成，球形燃料区内含有石墨基体和包覆颗粒，包覆颗粒含有燃料核芯和4个包覆层，包覆颗粒外面是石墨基体^[11]。
研究表明，标准的燃料球按正方体栅格和正六棱柱栅格无限阵列排列时，在最紧密的排列方式下，其栅格填充率分别为0.523 598和0.604 045，均小于0.61。为此，需要构造出一种特殊形状的燃料球来满足栅格填充率等于0.61的要求，其在正方体栅格中的剖面图如图 1所示。

图 1 特殊形状的燃料球的剖面图

图选项

特殊形状的燃料球的构造方法如下：燃料球的球形燃料区的结构保持不变，石墨球壳外径扩大后被正方体栅格剪切，石墨球壳外径和正方体栅格边长的选取需要同时满足石墨球壳的体积不变且栅格填充率等于0.61这2个要求。
石墨球的实际半径为3 cm，为满足栅格填充率等于0.61的要求，在MCNP中建模时也需要按照上述方法处理。
球床采用层重复结构描述，即构造2个基本层，其燃料球和石墨球在径向上按7：8的数量比例均匀排布，但是排布方式并不相同，交替重复2个基本层填充球床，使得在轴向上燃料球和石墨球也是均匀交错排布，从而避免产生燃料球柱或石墨球柱。
这种球床模型只是略微改变了燃料球、石墨球和球床空隙的形状，没有改变石墨球的数量，与实际情况最为接近，称为球床的真实模型。后续将在此模型的基础上模拟初装堆芯建立过程中堆芯空隙填充物质和顶空腔高度变化的情况。
2.2.2 等效模型为了模拟VSOP程序所采用的等效球体模型，在MCNP中建立同样的等效立方体模型，方法如下：燃料球外层的石墨球球壳外径扩大后被正方体栅格剪切以满足栅格填充率等于0.61的要求，栅格的剩余部分填充真空。通过重复这样的等效立方体模型即可产生球床的等效模型。
2.3 反射层模型在VSOP程序中，由于几何处理能力的限制，需要将侧反射层中的控制棒孔道、吸收球孔道、冷氦气孔道与周围石墨均匀化，将侧反射层中的控制棒与周围石墨均匀化后调整硼质量浓度来保证控制棒价值不变，将球床的顶锥和底锥等效成相同体积的圆柱，如图 2所示。

图 2 反射层物理计算分区模型

图选项

出于程序验证的需要，在MCNP中建立初装堆芯模型时，反射层结构与VSOP程序的结构完全一致，核素质量浓度 (特别是石墨中硼杂质质量浓度) 的分布也完全相同。
3 计算结果分析3.1 初装堆芯热谱通过在MCNP中建立初装堆芯球床的真实模型和等效模型，可以评估VSOP程序采用等效球体模型进行初装堆芯热谱计算的准确性。
MCNP的计算环境为真空，温度293.6 K，连续能量点截面库66C是在评价核数据库ENDF/B-VI Release 8的基础上加工得到，并采用了热化截面库GRPH.60T修正石墨的热散射效应，每次循环包含10000个粒子，共1000次有效循环。在0~1.9641 eV的热中子能区，采用了与VSOP程序相同的能群划分，初装堆芯热谱的计算结果如图 3所示。

图 3 初装堆芯热谱

图选项

由图 3可见，MCNP采用等效模型计算得到的初装堆芯热谱与真实模型的热谱完全一致，说明热谱主要取决于宏观物质组成，与石墨球和球床空隙的形状关系不大，因此VSOP程序采用等效球体模型计算初装堆芯热谱，尽管改变了石墨球和球床空隙的形状以及石墨球的数量，所得计算结果仍是准确的。
3.2 堆芯空隙填充物质改变在初装堆芯建立过程中，堆芯空隙会涉及到填充物质的改变，包括真空、0.1 MPa湿空气 (水蒸气密度为1.72066×10^-5 g/cm³)、0.1 MPa氮气和7 MPa氦气，在MCNP中分别模拟以上气氛，反应性价值的计算结果见表 1。
表 1 MCNP堆芯空隙填充物质的反应性价值

填充物质	keff	标准偏差	反应性价值
真空	1.17842	0.00029	无
0.1 MPa湿空气	1.16033	0.00029	-1.323%
0.1 MPa氮气	1.15591	0.00029	-1.653%
7 MPa氦气	1.18206	0.00030	?0.261%

表选项






由于空气和氮气中含有大量¹⁴N，会吸收中子，所以即使压力很低，相比较真空环境，k_eff仍然会有明显下降。氦原子吸收中子的能力较弱，作为轻核，中子慢化能力较强^[12]，并且由于反应堆正常运行时氦气压力比较高 (氦原子密度比较大)，使得堆芯能谱热化，从而引起k_eff提高，但是幅度非常小，意味着VSOP程序直接以真空气氛近似氦气气氛是合理的。以上效应对于初装堆芯建立过程中反应性的控制尤为重要。
VSOP程序的计算环境为温度293.15 K，多群常数库是在评价核数据库JEF-1、ENDF/B-IV和ENDF/B-V的基础上加工得到。VSOP程序的多群常数库中缺少核素氦的截面数据，无法计算氦气的反应性价值，对湿空气和氮气的反应性价值的计算结果见表 2。
表 2 VSOP堆芯空隙填充物质的反应性价值

填充物质	keff		反应性价值
	流修正	无流修正	流修正	无流修正
真空	1.155 92	1.169 88	无	无
0.1 MPa湿空气	1.138 43	1.151 95	-1.329%	-0.298%
0.1 MPa氮气	1.133 35	1.146 82	-1.723%	-0.686%

表选项






比较表 1和表 2的数据可见，如果不加入流修正，VSOP程序的计算结果与MCNP程序的差别较大，这是因为即使堆芯空隙填充了气体物质，但是由于压力比较低 (密度比较小)，对中子的散射作用不明显，部分中子仍然会通过球床空隙泄漏出去。加入流修正以后，VSOP程序计算的0.1 MPa空气和0.1 MPa氮气的反应性价值与MCNP程序计算的结果符合较好。
3.3 顶空腔高度变化分别在MCNP和VSOP程序中，将球床的高度从220 cm逐步提升至1 100 cm，研究Δk_eff=(k_eff，VSOP-k_eff，MCNP)/k_eff，MCNP的变化规律，计算结果如图 4所示。

图 4 Δkeff随球床高度逐渐增加的变化规律

图选项

由图 4可见，VSOP程序采用分方向扩散系数的方法对于小尺寸顶空腔的计算结果是非常准确的，对于大尺寸顶空腔的计算结果存在一定偏差，但是在核设计要求的1%允许范围内，因此可以用VSOP程序模拟HTR-PM初装堆芯建立过程。
4 结??论本文研究了VSOP程序进行初装堆芯热谱计算的方法，在MCNP中采用重复结构建立了初装堆芯球床的真实模型和等效模型，并将2种模型计算得到的热谱进行比较；研究了VSOP程序所采用的流修正和分方向扩散系数的方法，在MCNP中，以球床的真实模型为平台，模拟了初装堆芯建立过程中堆芯空隙填充物质和顶空腔高度变化的情况。
结果表明：对于初装堆芯，VSOP程序采用等效球体模型进行热谱计算是准确的，加入流修正后对堆芯空隙填充物质反应性价值的计算误差较小，直接以真空气氛近似氦气气氛是可行的，分方向扩散系数的方法对于大尺寸顶空腔也是适用的。
由于采用了大量的近似处理，VSOP程序计算效率高，适合于球床高温气冷堆的工程设计。对于全部由燃料球组成的平衡堆芯，VSOP程序已经经过广泛的验证，但对于由低富集度燃料球和石墨球均匀混合而成的初装堆芯，VSOP程序所采用的近似处理方法的适用性需要进行分析和验证。

参考文献

[1]	Journal of Central South University(Science and Technology), 41(2):649-654.-->ZHANG Zuoyi, WU Zongxin, WANG Dazhong, et al. Current status and technical description of Chinese 2×250 MWth HTR-PM demonstration plant[J]. Nuclear Engineering and Design, 2009, 239(7): 1212–1219. DOI:10.1016/j.nucengdes.2009.02.023
[2]	Journal of Central South University(Science and Technology), 41(2):649-654.--> Rütten H J, Haas K A, Brockmann H, et al. VSOP (99/05) Computer Code System for Reactor Physics and Fuel Cycle Simulation, Jül-4189[R]. Jülich, Germany: Forschungszentrum Jülich, 2005.
[3]	Journal of Central South University(Science and Technology), 41(2):649-654.--> International Atomic Energy Agency. Evaluation of High Temperature Gas Cooled Reactor Performance: Benchmark Analysis Related to Initial Testing of the HTTR and HTR-10, IAEA-TECDOC-1382[R]. Vienna, Austria: IAEA, 2003.
[4]	Journal of Central South University(Science and Technology), 41(2):649-654.--> Briesmeister J F. MCNP-A General Monte Carlo N-particle Transport Code Version 4C, LA-13709-M[R]. Los Alamos, USA: Los Alamos National Laboratory, 2000.
[5]	Journal of Central South University(Science and Technology), 41(2):649-654.--> Joanou C D, Dudek J S. GAM-I-A Consistant P1 Multigroup Code for the Calculation of Fast Neutron Spectra and Multigroup Constants, GA-1850[R]. San Diego, USA: General Atomic, 1961.
[6]	Journal of Central South University(Science and Technology), 41(2):649-654.--> Honek H C. THERMOS-A Thermalization Transport Theory Code for Reactor Lattice Calculation, BNL-5826[R]. New York, USA: Brookhaven National Laboratory, 1961.
[7]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Lieberoth J, Stojadinovic A. Neutron streaming in pebble beds[J]. Nuclear Science and Engineering, 1980, 76(3): 336–344.
[8]	Journal of Central South University(Science and Technology), 41(2):649-654.-->经荥清, 杨永伟, 古玉祥, 等. HTR-10初次临界装料预估[J]. 清华大学学报 (自然科学版), 2001, 41(4): 116–119.JING Xingqing, YANG Yongwei, GU Yuxiang, et al. Prediction calculation of HTR-10 fuel loading for the first criticality[J]. Journal of Tsinghua University (Science and Technology), 2001, 41(4): 116–119. (in Chinese)
[9]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Gerwin H, Scherer W. Treatment of the upper cavity in a pebble-bed high-temperature gas-cooled reactor by diffusion theory[J]. Nuclear Science and Engineering, 1987, 97(1): 9–19.
[10]	Journal of Central South University(Science and Technology), 41(2):649-654.-->常鸿, 杨永伟, 经荥清. 球床式高温气冷堆初次临界物理计算的蒙特卡罗方法模型分析[J]. 核动力工程, 2005, 26(5): 419–424.CHANG Hong, YANG Yongwei, JING Xingqing. Model analysis of Monte Carlo method for first criticality physics calculation in pebble bed high temperature gas-cooled reactor[J]. Nuclear Power Engineering, 2005, 26(5): 419–424. (in Chinese)
[11]	Journal of Central South University(Science and Technology), 41(2):649-654.-->经荥清, 杨永伟, 许云林. 蒙特卡罗方法用于HTR-10首次临界燃料装料预估的校算[J]. 核动力工程, 2005, 26(1): 28–34.JING Xingqing, YANG Yongwei, XU Yunlin. Application of Monte Carlo method for verification calculation in fuel loading prediction for first criticality of HTR-10[J]. Nuclear Power Engineering, 2005, 26(1): 28–34. (in Chinese)
[12]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Seker V, Colak V. HTR-10 full core first criticality analysis with MCNP[J]. Nuclear Engineering and Design, 2003, 222(2-3): 270–263.