基于子群方法的双重非均匀性共振计算方法研究

李 颂，张 乾,*，张志俭，赵 强，梁 亮，梁越超，张晋超，娄 磊，李满仓

(1.哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001；2.中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610213)

随着新型压水堆燃料的发展，双重非均匀性逐渐成为了堆芯计算的难点。传统燃料只包含燃料棒与慢化剂之间的一重非均匀性，而全陶瓷微密封(FCM)燃料中的TRISO以颗粒形式弥散分布在碳化硅基体中[1]，组成了第二重非均匀性。由于弥散颗粒具有随机排布的性质，传统共振方法无法进行双重非均匀性系统的共振计算。

为解决此问题，国际上主要提出了3种方法。第1种方法为丹可夫方法[2]，通过中子逃脱概率计算丹可夫因子来修正有效共振截面，主要用于高温气冷堆，不适用于FCM燃料球的多层结构。第2种方法为缺陷因子方法[3-6]，通过颗粒与基体的等效均匀化，将问题简化为传统一重非均匀系统。该方法简单高效，但无法给出颗粒内部的精细通量与截面分布。第3种方法为Sanchez-Pomraning方法[7]，可直接求解双重非均匀性系统特征线方法(MOC)输运方程，给子群共振计算提供了可行性。然而，子群方法难以处理双重非均匀性系统中的共振干涉，传统Bondarenko方法精度较低，干涉因子方法需建立复杂的双重非均匀性积分表，伪核素方法则需求解大量超细群方程，计算效率较低。为精确计算双重非均匀性系统的共振截面，本文提出耦合Sanchez-Pomraning方法的改进的子群方法(ISSP)。

1 理论模型

1.1 基于Sanchez-Pomraning方法的MOC模型

传统MOC无法直接对FCM燃料随机分布的弥散颗粒进行特征线扫描。为解决此问题，Sanchez-Pomraning方法提出了“更新方程”的概念。首先，通过迭代法计算基体和颗粒混合的等效总截面，如式(1)所示。基于等效总截面，平源区的出射通量如式(2)所示， 其中φas为渐进通量，计算公式如式(3)所示。平源区等效源项q=Σtφas。

(1)

(1-e-ΣtL)(φas-φin)

(2)

(3)

式(2)建立了传统MOC与双重非均匀性系统间的联系。对不含弥散颗粒的平源区，φas计算方法与传统MOC相同。在弥散燃料区域，基体通量φm可通过直接求解式(1)～(3)得到。基于基体通量，颗粒内部各层通量则通过式(4)求得。

(4)

1.2 改进的子群方法

与传统子群方法[8]相比，改进的子群方法采用精细能群结构来规避共振干涉处理。子群参数包括子群截面σi及子群权重ωi，可通过中间共振近似及帕德近似方法计算得到[9]。有效共振截面σx与子群参数间的关系为：

(5)

其中：g为当前能群；φi为子群通量，需求解子群固定源方程(式(6))得到。

(6)

其中：Ω和r分别为角度和位置变量；σt′,i=σt,i-(1-λ)σs,i，σt，i和σs，i分别为子群总截面和子群散射截面；Nr为共振核素r的核子密度；λn为核素n的中间共振因子；Σp，n为势散射截面。式(6)右侧为源项，左侧第2项为移出项。

在双重非均匀性条件下，子群移出项、源项与通量计算需要式(1)～(3)来进行修正。为减少固定源方程求解次数，本文结合HELIOS[10]和nTRACER[11]程序对子群固定源方程的优化，提出了等效单群固定源优化方法[9]。将所有共振能群等效为1个共振群，子群固定源方程根据等效共振群的子群截面进行计算，并通过插值得到实际子群通量。等效能群的子群截面计算方法由文献[9]提供，本文选取取值范围为10-24～10-20cm2的8个等效子群截面，可满足不同类型问题计算。

为保证输运计算效率，本文采用中子慢化方程来获得共振能群能谱。其移出项由细群总截面计算得到，源项来自上游能群的下散射。求出第1共振群的通量后，下游共振群的源项即可根据散射矩阵和上游能群通量逐群求得。第1共振群的散射源项来自于快群，其能谱由渐进通量1/E近似提供。在双重非均匀性条件下，慢化方程移出项、源项与通量计算同样需要式(1)～(3)来修正，并通过式(7)归并得到最终有效多群截面。

(7)

1.3 基于Sanchez-Pomraning方法的超细群方法

本文选取的超细群在共振能量段共有34 000群，中子与燃料核素碰撞后所能获得的最大对数能降远大于超细群的勒宽，因此在每个超细群能量段内中子最多只会碰撞1次。定义fk,n为中子飞越k个超细群后与核素n碰撞的概率，如式(8)所示，每个超细群的源项可由式(9)计算得到。

(8)

(9)

其中：N为核素总数；K为中子与核素n碰撞后能飞越的最大能群数；Δug为第g群勒宽；u0、u、u′为能量变量；Σn，s为散射截面；α=(A-1)2/(A+1)2，A为相对原子质量。

由式(8)可得fk=e-Δuk-1fk-1，将此关系代入式(9)可得超细群源项，如式(10)所示。因此，在求出第1共振群的通量后，下游能群源项即可逐群求解固定源方程获得。第1共振群的源项由渐进通量近似求解，结合式(1)～(3)求得所有超细群通量后，再根据给定宽群能群结构通过式(7)进行并群计算。

(10)

2 数值验证

本文使用ENDF-Ⅶ版数据库，采取HELIOS-1.11程序的47群能群结构[10]，其共振能量范围为1.855～9 118 eV。本文多群数据库由开源程序NJOY-2016[13]制作得到，问题基准解由开源蒙特卡罗程序OpenMC[14]提供。蒙特卡罗建模采用显式随机模型，粒子总数设置为100 000个，粒子代数为1 600代，前300代不计入统计。

2.1 FCM单栅元问题

1) 不同填充率

图1 FCM燃料几何结构示意图Fig.1 Scheme of geometry structure of FCM fuel

单栅元FCM燃料几何结构如图1所示。本算例采用富集度为17.8%的碳化铀燃料，各材料及几何信息由文献[15]提供，系统温度为300 K。本文依次选取1%、5%、10%、20%、30%和40%的填充率进行计算，有效增殖因数计算结果列于表1，并与基于Sanchez-Pomraning方法的超细群方法(UFGSP)的计算结果进行对比。

从表1可见，ISSP与UFGSP计算结果与基准解一致，各填充率偏差均小于100 pcm，而ISSP的效率更高。同时，表1分别列出了ISSP子群固定源和慢化方程两个模块的计算时间以及共振和本征值计算的总时间。子群时间在各填充率下均低于总时间的4%，而慢化方程时间最高占比为28.2%，因此ISSP的并群效率仍有可提升空间。以填充率10%和40%为例，图2示出不同填充率下238U和235U吸收截面的计算精度。

表1 不同填充率下keff的计算偏差Table 1 Calculating deviation of keff for different packing fractions

图2 不同填充率下的238U(a)和235U(b)吸收截面的相对偏差Fig.2 Relative deviation of absorption cross-section of 238U (a) and 235U (b) under different packing fractions

ISSP对40%填充率的238U吸收截面的相对偏差均在1%以内，对10%填充率的最大相对偏差上升至1.9%。由于偏差较大能群的截面远小于其他共振群，因此对最终结果影响很小。UFGSP同样有上述偏差变化趋势，且所有能群相对偏差小于1%。整体而言，ISSP和UFGSP均可较高精度地提供各填充率下的共振截面，而ISSP具有更高的计算效率。

2) 毒物问题

本文采用如图3所示的两种毒物颗粒，第1种是在TRISO燃料颗粒内增加1层毒物材料(QUADRISO)，第2种是毒物以单独毒物球(BISO)的形式与TRISO颗粒一起弥散在基体中。本文采用B4C和Gd2O3两种毒物材料，其中单颗粒问题的QUADRISO填充率为30%，双颗粒问题中的TRISO颗粒填充率为33.8%，BISO为10.8%。系统温度为300 K，材料核子密度与几何尺寸来自文献[15-16]。

a——QUADRISO；b——BISO图3 可燃毒物颗粒的几何结构Fig.3 Geometry structure of burnable poison

B4C毒物在QUADRISO和BISO两种形式下的238U吸收截面的相对偏差如图4所示，ISSP和UFGSP均具有较高的计算精度，所有能群的相对偏差均小于1%。QUADRISO的相对偏差偏大，主要是因为与燃料内核毗邻的毒物层导致颗粒内部的通量波动更为剧烈。含Gd毒物中的155Gd和157Gd对反应率影响较大，其吸收截面的相对偏差如图5所示。

图4 含B4C毒物燃料QUADRISO和BISO颗粒238U吸收截面的相对偏差Fig.4 Relative deviation of 238U absorption cross-section in QUADRISO and BISO with B4C poison

155Gd共振峰主要集中在10 eV以下能量段，该范围能群内ISSP和UFGSP对QUADRISO和BISO的最大截面相对偏差均低于2%。在29.02～47.81 eV及130.1～2 034.7 eV之间，155Gd的最大相对偏差增至3.3%左右，但由于此能量段截面很小，因此对全局反应率的影响也较小。ISSP和UFGSP对157Gd的计算精度和155Gd类似，其最大相对偏差在2.6%左右。毒物问题的keff计算精度列于表2，其偏差均在150 pcm以内。总体来看，ISSP和UFGSP均可对毒物材料进行准确的共振计算。

图5 155Gd(a)和157Gd(b)吸收截面的相对偏差Fig.5 Relative deviation of absorption cross-section of 155Gd (a) and 157Gd (b)

表2 毒物问题keff的偏差Table 2 Deviation of keff for burnable poison problem

3) 燃耗问题

以填充率为30%的算例为初始条件计算FCM燃料的燃耗过程，其中每个燃耗步的核子密度由OpenMC提供，并由ISSP和UFGSP进行共振和输运计算。有效增殖因数的偏差如图6所示，ISSP和UFGSP在大多数燃耗步的偏差在100 pcm以内。两种方法在燃耗过程初始计算精度相近，随燃耗深度的增加，共振核素大量累积导致干涉效应复杂，计算的偏差会有所增大，在90 GW·d/tU之后偏差逐渐平稳。从整体上看，ISSP和UFGSP均可较好地计算FCM燃料的燃耗过程。

图6 燃耗过程中keff计算的偏差Fig.6 Deviation of keff during burn-up process

2.2 组件问题

1) 常规FCM燃料组件

本文选取了13×13二维燃料组件，其中包括160根FCM燃料棒和9根导向管，如图7所示。其中，FCM燃料选自2.1节填充率为30%的算例，导向管材料与几何组成参考文献[15]，系统温度为300 K。以图7所标燃料棒1和2为例，其238U与235U吸收截面的相对偏差如图8所示。

图7 13×13 FCM燃料组件几何结构Fig.7 Geometry structure of 13×13 FCM fuel lattice

由图8可见，238U吸收截面的相对偏差在大部分能群小于1%，最大相对偏差分别为-1.54%和-1.42%。235U的计算精度较高，ISSP的最大相对偏差为1.10%，UFGSP在所有能群均小于1%。有效增殖因数和计算时间列于表3，ISSP和UFGSP计算的偏差分别为16 pcm和-1 pcm，二者均可对FCM燃料组件进行精确计算，而ISSP的计算时间远小于UFGSP。

图8 13×13组件燃料棒1、2的238U(a)和235U(b)吸收截面相对偏差Fig.8 Relative deviation of absorption cross-section of 238U (a) and 235U (b) in fuel pin 1 and 2 of 13×13 lattice

表3 FCM燃料组件问题keff的偏差Table 3 Deviation of keff for FCM fuel lattice

2) 含毒物FCM燃料组件

本文构造了4×4含毒物FCM燃料组件，其中包括12根常规FCM燃料棒、2根导向管和2根毒物棒，如图9所示。本算例毒物棒选自2.1节含Gd毒物的QUADRISO燃料，其余材料及几何均与13×13算例相同，系统温度为300 K。图9所标燃料棒1和2的238U以及燃料棒1中155Gd吸收截面的相对偏差如图10所示。

相比于常规FCM燃料组件，ISSP和UFGSP对毒物组件238U吸收截面计算的相对偏差有所增加，其最大偏差为2.59%，对大多数能群仍保持1%左右，这是由于通量在毒物棒附近剧烈变化所导致。155Gd截面精度也具有相似的规律。该组件本征值计算结果列于表3，相比于常规组件，该组件有效增殖因数偏差有所增加，但仍在可接受范围内。同时ISSP具有更高的计算效率。

图9 含Gd毒物4×4 FCM燃料组件几何结构Fig.9 Geometry structure of FCM 4×4 fuel lattice with Gd poison

图10 4×4组件燃料棒1、2的238U(a)和棒1的155Gd(b)吸收截面相对偏差Fig.10 Relative deviation of absorption cross-section of 238Uin fuel pin 1 and 2 (a) and 155Gd in fuel pin 1 (b) of 4×4 lattice

3 结论

为解决双重非均匀性系统的共振计算问题，提出了改进的子群方法。该方法通过精细共振能谱来规避共振干涉处理，并采用等效单群固定源优化来提高计算效率。通过调用基于Sanchez-Pomraning方法的MOC输运模块分别计算子群固定源方程和慢化方程，得到FCM燃料颗粒和基体等材料的多群有效截面，最后再调用该输运模块进行本征值计算。数值结果表明，ISSP可充分考虑双重非均匀性的影响，能精确并高效地计算FCM燃料多种工况下的共振截面。