易思宇，刘宙宇，尹 文

(西安交通大学 核科学与技术学院，陕西 西安 710049)

包覆燃料颗粒是一种广泛应用于高温气冷堆的燃料类型，具有良好的安全性。近年来，为提高压水堆中燃料的安全性，多个实验室开展了将这种燃料颗粒弥散到碳化硅基体中形成燃料密实体应用于先进压水堆的研究，这种燃料类型称为全陶瓷微胶囊封装(FCM)燃料。FCM燃料存在双重非均匀性，包括燃料密实体中颗粒和基质的微观非均匀性和燃料棒之间的宏观非均匀性，这给中子学计算带来了挑战。

对于确定论方法，双重非均匀性给共振和输运计算都带来了困难。在确定论方法中，空间被离散成多个平源区，且每个平源区都被视为均匀材料，但FCM燃料中所使用的TRISO颗粒尺寸过小，难以进一步离散。为解决此问题，人们最初提出了基于体积守恒的均匀化方法，将颗粒与基质打混处理，使双重非均匀性系统转化为单重非均匀性系统，但这样会由于忽略颗粒的空间自屏效应而引入巨大误差，为改善精度，一些基于确定论的双重非均匀性处理方法被提出。

Luo等[1]提出了反应率等价物理转换(RPT)方法，该方法将颗粒弥散在基质内的真实模型转化为内部为均匀燃料、外部包围一层基质的等效模型，并以此实现与真实问题相同的空间自屏效应。但这种方法需使用蒙特卡罗程序提前计算出真实问题的反应性，并根据反应性守恒来建立等效模型。

Sanchez等[2-5]和Pormraing[5]提出了基于统计学的方法，称为Sanchez-Pormraning方法。该方法广泛应用于中子学计算程序，如APPOLO、HELIOS、DeCART。其主要思想是使用概率论与数理统计中的数学期望的思想来建立通量的期望，然后引入弦长的概率分布来建立中子平衡方程，并将其称为更新方程。然后使用更新方程计算获得各平源区的均匀化截面和源项。该方法精度较高，但相比其他方法其计算效率较低。

其他研究者[6-8]提出了基于各种模型的等效均匀化方法，即缺陷因子方法。此方法易于实现且能有效处理双重非均匀性，在VOSP[9]、SCALE[6]和NECP-X[10]等程序中得到了使用。基于等效一维球模型的缺陷因子方法是其中常见的方法，该方法首先在一维或二维等效模型中求解慢化方程或碰撞概率以获得缺陷因子，然后通过缺陷因子将燃料棒中的颗粒与基质等效均匀化，最终实现双重非均匀性系统到单重非均匀性系统的转化。尽管研究人员基于缺陷因子方法进行了很多研究，但该方法仍存在一些问题。首先，计算缺陷因子的等效模型是通过体积权重守恒来建立的，这引入了颗粒均匀分布在无限基质中的假设，无法考虑宏观的非均匀性；其次，非共振能区的缺陷因子是通过在三维穿透模型中保证首次飞行碰撞概率守恒获得，但这种方法无法满足FCM燃料中高精度计算的要求。为在计算缺陷因子时考虑宏观非均匀性，并获得更准确的多群缺陷因子，本文在先前研究[10]的基础上提出一种改进的缺陷因子方法。

1 理论模型

1.1 概述

传统的压水堆燃料只需考虑燃料棒间的单重非均匀性，但FCM燃料还需考虑来自颗粒和基质间的非均匀性。早期的研究中，通过简单的体积权重守恒实现颗粒与基质的均匀化，但颗粒的光学厚度大于中子的平均自由程，且颗粒内部的通量会受空间自屏影响发生剧烈变化，简单的体积均匀化忽略了颗粒的空间自屏效应，会带来较大误差。为实现更准确的均匀化，提出了可衡量颗粒空间自屏效应的缺陷因子方法。缺陷因子方法的主要思想是通过缺陷因子来衡量颗粒与基质，以及颗粒内部的空间自屏效应，并实现等效均匀化。本文首先通过缺陷因子实现等效均匀化，然后使用Liu等[11]提出的全局-局部方法处理具有单重非均匀性的栅格系统中的共振自屏效应，最终获得适用于FCM燃料的多群有效自屏截面。为保证等效均匀化时的准确性，本文提出一种改进的缺陷因子计算方法。等效均匀化是缺陷因子方法的基础思想，其基于燃料密实体内的反应率守恒，如式(1)所示。

(1)

通过对式(1)简单转化，可得到式(2)：

(2)

缺陷因子是子区域的通量和平均通量的比值，但真实问题的通量分布是无法提前获得的，所以需构建等效模型，在等效模型中获得通量分布。

1.2 等效模型的建立

传统的缺陷因子方法通过体积权重守恒来建立等效一维球模型，即颗粒尺寸保持不变，基质外半径通过体积比计算获得，如式(3)所示。

(3)

其中：Rm为等效一维球模型中的基质外半径；Rt为TRISO颗粒的外半径；F为TRISO颗粒的填充率。在白边界条件下，体积权重等效模型可看作是由基质中的颗粒组成的无限阵列，如图1所示。该模型很明显无法考虑真实栅格系统中燃料棒之间的包壳和慢化剂(图2)。为考虑来自包壳和慢化剂的影响，本文提出一种基于颗粒丹可夫因子的等效一维球模型，以此来反映颗粒在真实栅格系统中的自屏效应。

图1 燃料颗粒组成的无限阵列Fig.1 Infinite array of fuel particels

图2 装载FCM燃料的真实栅格系统(a)和等效一维球模型(b)Fig.2 Real lattice system loaded with FCM fuel (a) and equivalent 1D sphere model (b)

1) 双重非均匀系统中颗粒丹可夫因子的计算

本文对颗粒丹可夫因子的计算参考Ji等[12]提出的方法，等效模型的推导都将引入颗粒包裹层和基质的均匀化，并将均匀化后的非燃料材料统称为基质。颗粒的丹可夫因子包含两部分：燃料棒间的丹可夫因子和燃料棒内的丹可夫因子。燃料棒内的丹可夫因子可通过式(4)计算获得。

(4)

其他变量可通过以下公式获得：

t=2(rp-rk)

(5)

(6)

(7)

(8)

(9)

其中：rp为颗粒的半径；rk为燃料核的半径；λ为基质中的平均自由程；〈l〉为无限基质中颗粒-燃料核之间的非燃料材料的平均弦长；frac′为燃料填充率；frac为颗粒填充率。

至此获得了所有棒内颗粒的丹可夫因子组成部分的表达式。燃料棒间颗粒的丹可夫因子可通过式(10)计算获得。

(10)

其中，Cb为燃料棒的丹可夫因子，可将燃料棒当作黑体，然后使用中子流方法获得。由此可通过式(4)计算获得燃料棒内颗粒的丹可夫因子，通过式(10)计算获得燃料棒间颗粒的丹可夫因子，将两者加和，即可获得双重非均匀性系统下颗粒的整体丹可夫因子，如式(11)所示。

CDH=Cintra+Cinter

(11)

2) 等效半径的推导

假设中子是均匀且各向同性地从燃料核外表面射出，白边界条件下等效一维球模型即可看作是由排布在基质中的燃料核组成的无限阵列[13](图1)，因此等效模型的丹可夫因子Ceq满足式(12)，即：

(12)

由式(7)、(12)可获得燃料的等效填充率frac′,eq，如式(13)所示。

(13)

而等效模型的外径Req和燃料填充率有如下关系：

(14)

因此，Req是颗粒丹可夫因子的函数，写作Req(Ceq)。一旦获得双重非均匀系统下颗粒的丹可夫因子，即可通过式(13)、(14)获得Req(Ceq)。由此即可建立与真实双重非均匀性系统具有相同丹可夫因子的等效一维球模型。

1.3 均匀化多群有效截面的计算

等效均匀化多群有效截面的计算分为共振能量段和非共振能量段两部分。这两个部分中，缺陷因子都是必需的参数，对于非共振能量段，本文提出了一种获得缺陷因子的新方法。

1) 共振能量段的均匀化多群有效截面

共振能量段的多群有效截面计算由等效均匀化和全局-局部共振计算两部分组成，具体流程如图3所示。本文使用全局-局部方法进行共振计算，其中局部计算使用超细群方法，因此需计算超细群缺陷因子来均匀化共振段的超细群截面。

图3 共振段等效多群有效截面的计算流程Fig.3 Calculation flowchart for equivalent multi-group effective cross sections in resonance energy range

(1) 等效均匀化

使用基于丹可夫因子的等效一维球模型，后续的计算与Yin等[10]提出的过程相同。通过碰撞概率法求解如式(15)所示的超细群慢化方程，以获得超细群通量。

(15)

其中：下标g为超细群能群，i和j为子区域编号；Qjg为体积平均慢化源项；Pjig为j区产生的中子在i区发生第1次碰撞的概率。

超细群缺陷因子可通过式(16)获得：

(16)

可通过式(17)获得超细群均匀化截面：

(17)

(2) 全局-局部共振计算

将燃料棒等效均匀化后，带有FCM燃料的栅格系统可被处理为典型的压水堆栅格系统，即燃料棒中的燃料密实体为均匀材料。对于具有单重非均匀性的典型栅格系统，使用全局-局部共振方法进行自屏计算。

首先通过中子流方法获得燃料棒的丹可夫因子Cb，在黑体假设下，有：

(18)

其中：φ0为孤立棒中燃料区的通量；φ为栅格系统中燃料区的通量。φ可通过MOC方法求解单群固定源问题获得。在黑体假设下，φ0等于中子从燃料到慢化剂的逃脱概率，孤立棒问题的逃脱概率可通过Carlvik方法获得。

然后根据丹可夫因子守恒搜索慢化剂的外径，如式(19)所示。

(19)

其中：rM为慢化剂的外半径；Pe为中子从燃料棒中逃脱的概率；Pf-m(rM)为在燃料中产生的中子在慢化剂中发生第1次碰撞的概率。

至此即可建立针对每根燃料棒的一维栅元模型，最后在等效一维圆柱栅元模型中通过超细群方法获得共振能量段的多群有效截面。

2) 非共振能量段的均匀化多群有效截面

非共振能量段无需进行共振计算，所以该能量段的多群截面只需使用多群缺陷因子进行等效均匀化即能获得均匀化多群有效截面。具体流程如图4所示。

图4 非共振段等效多群有效截面的计算流程Fig.4 Calculation flowchart for equivalent multi-group effective cross section in non-resonance energy range

在之前的研究中，非共振能量段缺陷因子的计算参考She等[8]的方法，其主要思想是在一个三维穿透模型中保证首次飞行碰撞概率守恒，所以只需计算中子在各子区域的首次飞行碰撞概率和穿透概率。She等的方法有较高的计算效率，但在填充率较低或较高的情况下计算精度不足。

为获得更准确的缺陷因子，本文采用在等效一维球模型中进行特征值计算的方法获得通量分布。首先需获得进行特征值计算必需的多群有效截面，因为前文已获得了共振区的超细群通量，所以只需归并截面即可获得有效共振截面，而非共振区可直接使用典型能谱下的有效截面。然后是求解特征值问题，本文使用碰撞概率法进行求解。

2 数值结果

从单栅元问题、二维单组件问题、三维全堆问题等多方面对本文提出的方法进行验证。并在单栅元规模下，对不同颗粒尺寸、填充率、带可燃毒物等工况进行验证。参考解来自蒙特卡罗程序，为获得准确结果，蒙特卡罗计算时采用显示建模的方法。蒙特卡罗程序使用的连续能量库和NECP-X使用的截面库皆基于ENDF/Ⅶ.0。NECP-X的多群结构为WIMSD69，其中共振能群为13～45群。各区域的温度，如果没有提及，均为600 K。

2.1 单栅元问题

单栅元和TRISO颗粒的几何布置示于图5，其尺寸和材料列于表1。材料的核子密度列于表2。在蒙特卡罗计算中，使用30万个粒子，共500代进行模拟，其中前200代为非活跃代。由于TRISO颗粒在尺寸、填充率、富集度等方面都有多种选择，为方便讨论，设计了表3所列的一些常见工况的颗粒参数。

图5 FCM燃料栅元和TRISO颗粒Fig.5 FCM pin cell and TRISO particle

表1 FCM燃料栅元的尺寸和材料Table 1 Dimension and material of FCM cell

表2 材料的核子密度Table 2 Atomic number density of each material

表3 一些常见工况的颗粒设计参数Table 3 Particle design parameters for some common cases

1) 颗粒尺寸和填充率

为更全面地反映方法的几何适应性，在工况1的基础上通过改变燃料半径和疏松热解碳层厚度设计了表4所列的系列不同颗粒尺寸工况，其中包含了工况2～5。验证在15%、25%、35、45%的填充率下进行。

表4 用于几何适应性测试的TRISO颗粒尺寸Table 4 TRISO particle sizefor geometric adaptivity test

为说明方法的改进，分别采用新旧方案对问题进行计算。其中新方案为使用本文提出的基于丹可夫因子的等效一维球模型与求解特征值问题获得非共振区缺陷因子的计算过程，而旧方案为使用体积权重等效一维球模型和在三维穿透模型中计算非共振区缺陷因子的计算过程。两方案计算结果与蒙特卡罗结果的偏差列于表5。

从表5可见，旧方案燃料核直径为0.08 cm、填充率为0.45的工况中，kinf的绝对偏差超过300 pcm，即对于有较大燃料核直径、较大填充率和较小疏松热解碳层厚度的工况，旧方案无法获得准确的计算结果。而新方案中这几个工况都得到了极大的改善。对其中的工况1、3、4进行深入分析，具体计算结果列于表6。

表5 新旧方案中不同颗粒尺寸和填充率下kinf的绝对偏差Table 5 Absolute bias of kinf in different particle sizes and packing fractions with new and old schemes

表6 部分工况的kinf和等效半径Table 6 kinf and equivalent radius in part of cases

在新方案中，基于丹可夫因子的等效一维球模型能考虑燃料棒外慢化剂的影响，因此相比旧方案有更大的基质外径来反映更加真实的屏蔽效应。将新旧方案中的238U和235U有效共振自屏截面与蒙特卡罗结果进行比较，结果如图6所示。

图6 238U和235U共振能量段有效吸收截面相对偏差比较Fig.6 Comparison of relative bias of effective absorption cross section of 238U and 235U

由图6可见，相比旧方案，新方案在20～27群的238U吸收截面明显减小且更接近参考值，整个共振能量段最大相对偏差小于2.5%。235U的吸收截面有类似情况，20～27群的截面得到了改善，虽然最后几群的截面更加偏离参考值，但变化很小，且最大相对偏差都小于1%。新方案中的等效模型更大的基质外径使得其共振段吸收截面减小，增强了空间自屏效应，使燃料核中的缺陷因子减小，而缺陷因子的减小会反映到吸收截面上，最终吸收截面减小使得kinf增大。这也能解释在旧方案中当燃料核直径增大时会出现kinf往负方向偏离参考值的趋势，且这种趋势会随着颗粒填充率的增大、疏松热解碳层变薄更加明显。这是因为增大燃料核直径和填充率、减小疏松热解碳层厚度都会使燃料占比增大，这时体积权重模型导致的忽略慢化剂所带来的影响就会被放大。

除共振能量段的改进，新方案中热能区的宏观截面同样也得到了改善，如图7所示，燃料密实体中宏观裂变截面和宏观总截面相比旧方案都有所减小，且更接近参考值，这进一步说明了新方案中获得多群缺陷因子方法的准确性。

图7 燃料密实体宏观裂变截面和总截面相对偏差比较 Fig.7 Comparison of relative bias of macro fission and total cross sections in fuel lump

2) 燃料温度和富集度

工况1、6、7是关于温度和富集度的敏感性测试，kinf的计算结果列于表7。可见，在不同温度和富集度下kinf与参考值符合得较好。不同工况下235U和238U的有效截面相对偏差比较示于图8。可见，238U有效截面的最大相对偏差小于2%，235U有效截面的最大偏差小于1%，均与蒙特卡罗计算所得的参考解符合较好。

图8 新方案中不同温度及富集度下238U和235U共振能量段有效吸收截面相对偏差比较Fig.8 Comparison of relative bias of effective absorption cross section of 238U and 235U in different termperatures and enrichments in new scheme

表7 不同温度和富集度下kinf和绝对偏差Table 7 kinf and absolute bias in different temperatures and enrichments

3) 燃料和毒物燃耗

为验证本文提出的方法在燃耗计算时的准确性，基于工况4进行燃耗计算，并在工况4的基础上设计了带毒物的工况。在毒物工况中，使用QUADRISO颗粒替代TRISO颗粒。QUADRISO颗粒即在TRISO颗粒的燃料核外额外包裹1层毒物，本文使用10 μm厚的Er2O3。相关的颗粒尺寸列于表8。

表8 QUADRISO颗粒尺寸Table 8 Size of QUARISO particle

整个燃耗过程从0 GW·d/tU到72.5 GW·d/tU。整个燃耗过程中kinf及其绝对偏差的变化示于图9。对于没有毒物的工况，kinf的绝对偏差从最开始的-110 pcm到20 GW·d/tU的10 pcm，随后缓慢增长到72.5 GW·d/tU的30 pcm。其中最大绝对偏差为-120 pcm。对于有毒物的工况，kinf的绝对偏差从最开始的-80 pcm增长到最后的-230 pcm。相比没有毒物的工况，毒物的出现使得偏差略有增加，但仍在可接受范围内。

图9 无毒物和有毒物燃料栅元在燃耗情况下的kinf及其绝对偏差Fig.9 kinf and its absolute bias in depletion for case without and with poison

2.2 二维组件问题

本文设计了2个二维的13×13组件，第1个组件由带TRISO颗粒的燃料棒和导向管组成，第2个组件将第1个组件中部分燃料棒中的TRISO颗粒替换为QUADRISO颗粒，几何布置如图10所示。TRISO颗粒尺寸同工况4，QUATRISO颗粒尺寸如表8所列。材料组成和几何参数列于表9，材料的核子密度如表2所列。

图10 无毒物组件和有毒物组件的几何布置Fig.10 Geometry configuration of fuel lattice without and with poison

表9 导向管尺寸和材料Table 9 Dimension and material of guide tube

将kinf和功率分布皆与蒙特卡罗结果进行比较，结果见表10和图11。可看到，两个工况下kinf皆与参考值符合较好。对于不含毒物的组件，其栅元功率最大相对偏差为0.31%；在含毒物的组件中，栅元功率最大相对偏差为0.32%。

图11 有无毒物组件的功率分布及相对偏差Fig.11 Power distribution and relative bias for lattice with and without poison

表10 2D组件问题的kinf及其绝对偏差 Table 10 kinf and absolute bias of FCM lattice

2.3 三维全堆问题

根据Awan等[14]设计的装载FCM燃料的小型模块化堆，本文建立了1个有4类共37个组件的三维全堆问题，其中2类只装载了TRISO颗粒，其燃料分别为14.3%和10%的UC。另外2类中分别将前两者中部分燃料棒的TRISO颗粒替换为含有5 μm和10 μm厚毒物层的QUADRISO颗粒。具体组合列于表11。组件由13×13栅元构成，具体几何布置如图10所示。颗粒和栅元的尺寸基于2.1节中的工况1，所用材料如表2所列。

表11 4类组件的配置Table 11 Configuration of 4 kinds of assemblies

为简化计算，将巨大的环形反射层简化为一圈20 cm厚的反射层，轴向反射层被简化为上下各20 cm厚的反射层，如图12所示。至此建立起了一个由160 cm高活性区、40 cm厚反射层组成的简化全堆。在NECP-X计算中，将模型分为20层，每层10 cm。

图12 全堆和简化全堆的俯视图和侧视图Fig.12 Plan view and side view of whole and simplified FCM fuel core

使用本文提出的改进缺陷因子方法进行计算，以蒙特卡罗结果为参考。在特征值上，采用改进缺陷因子方法所得的keff为1.227 81，蒙特卡罗结果为1.226 67，两者间的绝对偏差为114 pcm，符合得很好。径向和轴向功率分布与参考解的比较如图13、14所示。轴向功率的最大相对偏差出现在最顶部，为2.28%，但两端的功率较低，所以偏差是可接受的。组件径向功率的相对偏差最大为-1.13%，出现在最中心的组件，整体符合较好。图15为有最大棒功率相对偏差的组件中的功率分布，其中棒功率最大相对偏差为1.88%，整体与参考解符合得较好。

图13 归一化轴向功率分布和相对偏差Fig.13 Normalized axial power distribution and relative bias

图14 全堆径向组件功率分布和相对偏差Fig.14 Radial assembly power distribution and relative bias in core

图15 组件内径向栅元功率分布和相对偏差Fig.15 Radial pin power distribution and relative bias in assembly

3 结论

本文提出了一种改进的缺陷因子方法用于获得FCM燃料中的有效多群截面。首先提出了基于丹可夫因子的等效一维球模型来代替体积权重模型，传统的体积权重模型忽略了燃料棒之间的非均匀性，因此当燃料密实体中燃料体积占比较大时会引入较大误差，而较大的燃料体积占比在FCM燃料设计中很常见，新的等效模型能考虑燃料棒间的非均匀性，可获得更精确的计算结果。其次在非共振段处理中，采用在等效模型中进行特征值计算的方法来获得多群缺陷因子，相比穿透模型的碰撞概率计算，该方法可获得更精确的结果。

该方法已在NECP-X中实现，并在不同工况和不同规模的问题中进行了验证。在单栅元问题中，计算了不同颗粒尺寸、填充率等工况，kinf偏差均在100 pcm以内，有效多群截面除少数能群相对偏差接近2.5%，大多数能群都在±1%以内。在二维组件问题中，kinf偏差在103 pcm以内，栅元功率最大相对偏差为0.32%。在三维全堆问题中，keff偏差为114 pcm，除两端功率较低外，其余活性区层轴向功率相对偏差均在1.5%以内，径向组件功率最大相对偏差为1.14%，径向积分栅元功率相对偏差为1.88%。所有计算结果与参考解符合较好，证明本方法能有效处理FCM燃料中的双重非均匀性，并获得了精确的有效多群截面。