基于等效几何理论的复杂几何燃料共振计算方法研究

曹 璐，刘宙宇，贺清明，曹良志

(1.西安交通大学 核科学与技术学院，陕西 西安 710049；2.西北核技术研究所，陕西 西安 710024)

传统共振计算方法包括等价理论[1]、子群方法[2]和超细群方法[3]，这些方法在处理复杂几何燃料的共振计算时均存在各自的问题。等价理论存在几何适应性问题，虽然通过三项有理近似可近似处理复杂几何的逃脱概率，但计算精度不够[4]；子群方法，理论上几何适应性取决于子群固定源方程求解器的几何处理能力，而随着计算机性能的提升，特征线方法(MOC)广泛应用于各确定论软件中，因此子群方法可借用特征线方法的几何适用性处理复杂几何燃料，但由于子群方法计算中需要利用共振积分表，采用均匀的积分表计算精度不够，而采用非均匀的共振积分表，需要提前针对几种固定几何的燃料进行共振积分表的计算，因此也存在几何适应性问题；直接超细群方法仅用于小规模问题的计算，如用于作共振积分表的均匀问题、一维燃料棒和一维燃料板，这是由于在应用超细群方法时，需要计算区域之间的碰撞概率，对于一维燃料棒或一维燃料板问题，其区域之间的碰撞概率可采用解析解获得，而对于复杂几何燃料，碰撞概率的计算存在困难。国际上有研究[5]将特征线方法应用于超细群慢化方程的求解，以解决超细群方法的几何适应性问题，但在应用该方法时，又需要提前作表用于加速，而作表又使得该方法只能适用于作表中采用的几何针对高保真数值反应堆技术对计算精度和计算效率的要求，国际上对全堆芯高效、高精度共振计算方法开展了大量的研究[6-7]。NECP-X[8]是由西安交通大学核工程计算物理实验室(NECP)自主开发的数值反应堆物理计算程序，已实现棒状、板状等几何形式燃料的全堆芯高精度共振计算[9-10]，但对于非规则几何燃料(如蜂窝煤状燃料)，由于采用特征线方法[11]求解伪核素子群方程，因此计算效率很低。本研究提出一种基于等效几何理论的共振计算方法用于处理各种几何形式的燃料，并将该方法应用于NECP-X，以解决复杂几何燃料的计算效率问题。

1 理论模型

1.1 基于等效几何理论的共振计算方法

共振计算的目的是获得共振能区的中子注量率，归并获得共振能量段的有效自屏截面，该过程如式(1)所示：

(1)

式中：下标x和g分别为反应道和能群编号；r、Ω和E分别为空间变量、角度变量和能量变量；ΔEg为能群的能量区间；σ和φ分别为微观截面和中子注量率。

共振计算的难点在于能群内中子通量密度的获得，对于多区系统，中子通量密度可通过慢化方程求解：

Σt,f(E)φf(E)Vf=Pf→f(E)Sf(E)Vf+

(2)

式中：f为目标燃料区编号；f′为其他燃料区编号；m为慢化剂区编号；E为能量点，eV；F为所有燃料区集合；M为所有慢化剂区集合；Σt,f为燃料区的宏观总截面，cm-1；φf为目标燃料区的中子通量密度，cm-2·s-1；V为区域体积，cm3；S为区域中子源强度，cm-3·s-1；Pf′→f为从其他燃料区产生的中子到目标燃料区第1次发生碰撞的概率；Pm→f为从慢化剂区产生的中子到目标燃料区第1次发生碰撞的概率。

基于互易公式，对慢化剂区源项采用窄共振近似，式(2)可写为：

Σt,f(E)φf(E)=[1-Pf→M(E)]Sf(E)+

Pf→M(E)Σt,f(E)/E

(3)

式中：Pf→M(E)为目标燃料区产生的中子到慢化剂区第1次发生碰撞的概率。

由式(3)可知，对于不同几何的燃料，如果能保证目标燃料区产生的中子在慢化剂中发生首次碰撞的概率Pf→M(E)相等，则能保证不同系统中的中子注量率相等，由式(1)可知，不同燃料的有效自屏截面将相同。

另外，根据丹可夫修正因子的定义可得：

(4)

式中：C为丹可夫修正因子；Pe,iso(E)为孤立系统下目标燃料区的逃脱概率。

所以Pf→m(E)可表示为：

Pf→m(E)=[1-C(E)]Pe,iso(E)

(5)

因此，如果保证不同几何燃料中的(1-C(E))和Pe,iso(E)相等，则不同几何燃料可获得相同的有效自屏截面。

1.2 逃脱概率的计算

(6)

后续有研究采用贝尔因子aB对式(5)进行改进：

(7)

在轻水堆计算中，贝尔因子的选区范围为1.1～1.4。

为进一步提高计算精度，Carlvik等[13]针对棒状几何燃料提出了两项有理近似公式：

(8)

针对板状几何燃料，Roman等[14]同样给出了两项有理近似表达式：

(9)

但对于复杂几何的燃料，难以获得简单直接的逃脱概率计算公式，而是通过求解慢化方程获得。孤立系统中，根据中子守恒可得：

Σt,f(E)φf(E)Vf=Pf→f(E)Sf(E)Vf+

Pm→f(E)Sm(E)Vm

(10)

采用互易关系，式(10)可转化为：

Σt,f(E)φf(E)=Pf→f(E)Sf(E)+

(11)

对于孤立系统，Pf→m(E)=Pe,iso(E)，所以式(10)可写为：

(12)

通过上述推导可知，在孤立系统下可通过计算燃料区的中子注量率及式(12)获得燃料区的逃脱概率。所以对于复杂几何燃料可构造孤立系统问题，然后采用可处理复杂几何的MOC求解该孤立问题用以获得燃料区中子注量率进而获得该几何燃料的逃脱概率，再将其转化为逃脱概率相同的一维圆柱(或平板)燃料模型。

以蜂窝煤状燃料为例，基于逃脱概率守恒将蜂窝煤状燃料等效为棒状燃料的过程示于图1。对于蜂窝煤状燃料，采用MOC求解该燃料的孤立问题获得该燃料的逃脱概率Pe,iso,1(E)，然后进行不同半径下相同材料的棒状燃料进行孤立问题的计算，获得Pe,iso,2(E,r)，当Pe,iso,2(E,r)=Pe,iso,1(E)时所获得的燃料半径即为等效棒状燃料的半径。

图1 基于等效几何理论的蜂窝煤状燃料的等效过程

1.3 燃料外围结构材料的等效处理

燃料外围的结构材料会影响燃料区的能谱，因此需要考虑结构材料的影响。对于复杂几何燃料外围的结构材料，通过保证燃料中产生的中子到结构材料的首次碰撞的概率守恒，在等效一维圆柱燃料外围建立一层等效的结构材料。

在孤立系统中针对包壳区域建立中子平衡方程，可得：

Σt,clad(E)φclad(E)Vclad=Pf→clad(E)Sf(E)Vf+

(13)

式中：Σt,clad(E)为目标包壳区总截面；clad为包壳材料集。

基于互易关系(式(14))，式(13)可写为式(15)。

Σt,clad(E)Pclad→f(E)Vclad=Pf→clad(E)Σt,f(E)Vf

Σt,clad(E)Pclad→m(E)Vclad=Pm→clad(E)Σt,m(E)Vm

(14)

Pf→clad(E)=(Σt,clad(E)φclad(E)Vclad-

(15)

从式(15)可看出，若要保证不同几何下Pf→clad(E)守恒，则需要保证不同几何下φclad(E)、Vclad/Vf、Pclad→m(E)Sm(E)/Σt,m、Pclad1→clad(E)·Sclad1(E)以及Sf(E)相等。对于原来的复杂几何燃料的孤立问题，Pf→clad(E)可通过特征线方法获得。

基于燃料区到包壳材料区碰撞概率守恒实现等效为棒状燃料(蜂窝煤状燃料)周围包壳区的过程示于图2。首先对图2a目标燃料周围的包壳采用MOC求解该燃料的孤立问题获得该燃料到包壳的碰撞概率，然后在1.2节中已获得的燃料棒半径的基础上，基于燃料到周围结构材料碰撞概率守恒采用二分法进行等效结构材料半径的搜索，获得的等效结构材料几何如图2b所示。

图2 基于等效几何理论的蜂窝煤状燃料的包壳区等效过程

1.4 丹可夫修正因子的处理

从1.1节可知，等效一维几何燃料的丹可夫修正因子也应等于复杂几何燃料在堆芯中的丹可夫因子。首先采用中子流方法获得原复杂几何燃料的丹可夫修正因子C(E)，然后搜索等效一维几何燃料外围的慢化剂尺寸，使C(pitch,E)=C(E)，C(pitch,E)为等效后几何在慢化剂尺寸为pitch时的丹可夫修正因子。由上文可知，丹可夫修正因子是能量相关的，即灰体丹可夫修正因子，有研究[15]表明，采用黑体丹可夫修正因子C(Black)和灰体丹可夫修正因子C(E)时，1-C(Black)较1-C(E)偏差小，对最终有效自屏截面的计算影响很小，所以本文采用黑体丹可夫修正因子进行模拟计算。

基于燃料区丹可夫修正因子守恒实现等效为棒状燃料(蜂窝煤状燃料)外围慢化剂的过程示于图3。根据图3a目标燃料全局计算获得的丹克夫修正因子，在1.2、1.3节中获得的燃料棒与周围结构材料半径的基础上，基于丹可夫修正守恒采用二分法进行等效慢化剂半径的搜索，获得的等效结构材料几何如图3b所示。

图3 基于等效几何理论的蜂窝煤状燃料的慢化剂等效过程

2 数值验证

本文在NECP-X程序中实现基于等效几何理论的复杂几何燃料共振计算方法。为分析计算精度，对多种非棒状栅元、组件问题进行计算分析。NECP-X采用基于ENDF/B-Ⅶ.0的拓展共振能群范围的WIMS格式69群核数据库，利用蒙特卡罗程序获得参考解。NECP-X程序的计算条件为：1/4卦限内有8个方位角和3个极角；特征线间距为0.01 cm；特征值收敛限为10-6；裂变率收敛限为10-5。为验证本文方法的计算精度和计算效率，将本文方法计算结果与可处理复杂几何燃料的传统子群方法以及基于MOC的伪核素子群方法[10]的计算结果进行比较。

2.1 非棒状燃料栅元问题

1)环板燃料栅元问题

环板燃料栅元问题燃料几何取自于ATR反应堆，该环板燃料栅元问题的几何以及平源区网格示于图4，计算中省略包壳的建模，燃料采用3.1%二氧化铀燃料，慢化剂采用轻水，环板燃料内外径分别为7.7 cm和7.750 8 cm，栅元长为6 cm，宽为1.6 cm，所有材料温度均为600 K，采用全反射边界条件，蒙特卡罗程序统计特征值结果为0.404 19±5.0×10-5。不同方法计算得到的235U和238U微观吸收截面及其偏差示于图5，不同共振计算方法的计算时间及特征值计算结果列于表1。

图4 环板燃料栅元问题的几何以及平源区网格

由图5可见，基于MOC的伪核素子群方法计算的235U微观吸收截面(σ)最大偏差发生在第27群，为1.69%，其他能群偏差均小于0.54%，238U微观吸收截面最大偏差也发生在第27群，为1.56%，基于等效几何理论的共振计算方法计算的235U微观吸收截面偏差均小于0.9%，但相对于基于MOC的伪核素子群方法，238U微观吸收截面偏差更大，最大偏差接近3.0%，传统子群方法计算的235U微观吸收截面最大偏差为5.67%，238U微观吸收截面最大偏差为-1.75%，可看出，3种方法中基于MOC的伪核素子群方法精度最高，其次为基于等效几何理论的共振计算方法，最后为传统子群方法，但从表1的特征值可看出，3种方法的特征值计算结果相近，特征值偏差分别为84 pcm、62 pcm、52 pcm，上述截面及特征值结果表明，基于等效几何理论的共振计算方法和基于MOC的伪核素子群方法在处理环板燃料共振计算中具有较好的计算精度，而从表1的共振计算时间可看出，采用基于等效几何理论的共振计算方法的计算速度较基于MOC的伪核素子群方法快近80倍，较传统子群方法快140.83倍，表明基于等效几何理论的共振计算方法具有较高的计算效率。

图5 环板燃料栅元问题中235U和238U吸收截面及偏差的计算结果

表1 环板燃料栅元问题中不同共振计算方法的计算时间及特征值计算结果

2)蜂窝煤状燃料栅元问题

本文根据文献[16]设计蜂窝煤状燃料栅元问题，燃料采用3.1%二氧化铀燃料，慢化剂采用轻水，所设计的蜂窝煤状燃料栅元问题的几何以及计算时采用的平源区网格示于图6，蜂窝煤状燃料内部的每个水洞尺寸均为0.05 cm，共分布3圈19个水洞，3圈水洞圆心所在半径分别为0.0、0.15、0.3 cm，蜂窝煤状燃料外径为0.4 cm，栅距为1.2 cm，所有材料温度均为600 K，计算时采用全反射边界条件。采用不同共振计算方法及蒙特卡罗程序计算235U和238U微观吸收截面及其偏差，结果如图7所示。由图7可见，基于等效几何理论的共振计算方法计算的235U的微观吸收截面最大偏差为-0.48%，发生在第27群，238U的微观吸收截面最大偏差为1.61%，发生在第24群；基于MOC的伪核素子群方法计算的235U的微观吸收截面最大偏差为0.65%，发生在第27群，238U的微观吸收截面最大偏差为1.88%，发生在第24群；传统子群方法计算的235U的微观吸收截面最大偏差为7.32%，发生在第25群，238U的微观吸收截面最大偏差为1.31%，发生在第24群。特征值计算结果列于表2，蒙特卡罗程序统计结果为1.373 45±4.0×10-5，基于等效几何理论的共振计算方法与基于MOC的伪核素子群方法计算的特征值偏差均小于50 pcm，而传统子群方法计算的特征值偏差为-164 pcm；基于等效几何理论的共振计算方法计算速度更快，相较于基于MOC的伪核素子群方法加速比为72.6，相较于传统子群方法加速比为152.87。上述结果表明，基于等效几何理论的共振计算方法在计算蜂窝煤状燃料时有较好的计算精度和较高的计算效率。

图6 蜂窝煤状燃料栅元几何及平源区网格

图7 蜂窝煤状燃料栅元问题中235U和238U吸收截面及偏差的计算结果

表2 蜂窝煤状燃料栅元问题中不同共振计算方法的计算时间及特征值计算结果

2.2 非棒状燃料组件问题

1)环板燃料组件问题

根据2.1节栅元问题建立环板燃料组件问题，共采用7个环板栅元，栅元几何如图8所示，所有材料的温度均为600 K，计算时采用全反射边界条件。采用蒙特卡罗程序统计特征值和图8中1号环板燃料的微观吸收截面，NECP-X计算条件与环板燃料栅元问题一致，采用不同方法计算的1号环板燃料235U和238U微观吸收截面及其偏差示于图9。由图9可见，基于等效几何理论的共振计算方法最大偏差发生在238U第23群，为-3.2%，基于MOC的伪核素子群方法最大偏差也在该能群，为-1.79%，而传统子群方法最大偏差发生在235U第25群，为5.57%。基于MOC的伪核素子群方法精度最高，基于等效几何理论次之，传统子群方法最差。特征值计算结果列于表3，蒙特卡罗程序统计结果为0.404 48±2.0×10-5，3种方法的特征值计算结果相差不大，而基于等效几何理论的共振计算方法的计算速度为基于MOC的伪核素子群方法的203.6倍，为传统子群方法的319.77倍。上述微观截面及特征值计算结果表明，基于等效几何理论的共振计算方法在环板燃料的共振计算中具有良好的计算精度，较基于MOC的伪核素子群方法具有非常明显的加速效果。

表3 环板燃料组件问题中不同共振计算方法的计算时间及特征值计算结果

图8 环板燃料组件几何

图9 环板燃料组件问题中235U和238U吸收截面及偏差的计算结果

2)蜂窝煤状燃料组件问题

根据2.1节中栅元问题建立蜂窝煤状燃料组件问题，其组件几何如图10所示，所有材料温度均为600 K，计算时采用全反射边界条件，NECP-X计算条件与蜂窝煤状燃料栅元问题一致，除统计组件内所有栅元的功率外，还统计了图10中数字2处燃料栅元(记为栅元2)的235U和238U的微观吸收截面。采用不同方法计算的燃料栅元2的235U和238U微观吸收截面及其偏差示于图11，各方法的计算时间和特征值计算结果列于表4。由图11可知，基于等效几何理论的共振计算方法截面计算结果最大偏差为-1.72%，基于MOC的伪核素子群方法计算最大偏差发生在238U第25群，为1.69%，传统子群方法截面计算结果最大偏差发生在235U第25群，为7.47%。从表4的特征值计算结果可看出，基于等效几何理论的共振计算方法的特征值偏差为0 pcm，基于MOC的伪核素子群方法的特征值偏差为-35 pcm，传统子群方法的特征值偏差为-171 pcm，从计算时间来看，基于等效几何理论的共振计算方法相较于基于MOC的伪核素子群方法加速比为380.05，相较于传统子群方法加速比为211.47，基于MOC的伪核素子群方法较传统子群方法时间更长的原因在于，基于MOC的伪核素子群方法进行等效计算时需要多次采用MOC搜索二维等效问题，当燃料栅元增多时，基于MOC的伪核素子群方法的计算时间将会进一步增加。3种方法计算的功率分布基本一致，其中基于等效几何理论的共振计算方法的计算功率与参考解的偏差如图12所示，燃料栅元功率在靠近导向管栅元及仪表管栅元处较高，栅元功率最大偏差为-0.128%，发生在组件边缘处，可看出，由于蒙特卡罗程序统计的偏差，栅元功率存在不对称的情况，栅元功率RMS为0.043%。截面计算结果、特征值结果和功率分布结果均表明，基于等效几何理论的共振计算方法在处理蜂窝煤状燃料组件问题时精度较高，而传统子群方法截面计算结果较差，与基于MOC的伪核素子群方法相比，基于等效几何理论的共振计算方法在保证较高的微观吸收截面计算精度的情况下，计算时间大幅减少，与传统子群方法相比，基于等效几何理论的共振计算方法具有更高的计算精度和计算效率。

图10 蜂窝煤状燃料组件几何

图11 蜂窝煤状燃料组件问题中235U和238U吸收截面及偏差的计算结果

表4 蜂窝煤状燃料组件问题中不同共振计算方法的计算时间及特征值计算结果

图12 蜂窝煤状燃料组件问题栅元功率及偏差分布计算结果

2.3 JRR-3全堆芯问题

采用本文提出的基于等效几何理论的共振计算方法对JRR-3全堆芯问题进行模拟计算，几何与材料参数参见文献[17]，采用二维控制棒全提问题进行模拟，模拟中简化外部轻水反射层的建模，仅建模到重水反射层。

特征值参考解为1.374 20±3.6×10-5，NECP-X特征值计算结果为1.377 58，特征值偏差为338 pcm，相比上述问题，该问题的各项异性更强。JRR-3组件功率计算结果示于图13，网格功率及偏差分布计算结果示于图14，其中的标准燃料组件燃料板长度方向等距划分20个网格，跟随燃料组件等距划分16个网格。

由图13可见，最大组件功率偏差为-1.182%，组件功率分布呈中间低两边高的趋势。由图14a可见，网格高功率发生在组件角点附近，由图14b可见，虽然有些网格的功率偏差超过5%，但大部分功率偏差均小于±3%，网格功率RMS为1.088%。考虑到板状燃料研究堆的强各向异性，特征值与功率计算偏差在可接受范围内。

图13 JRR-3组件功率及偏差分布计算结果

图14 JRR-3网格功率及偏差分布计算结果

3 结论

针对全局-局部耦合策略中非棒状燃料的计算问题，提出了基于等效几何理论的复杂燃料共振计算方法，并将其应用于NECP-X程序中。该方法通过基于复杂几何燃料逃脱概率、复杂几何燃料到周围结构材料的碰撞概率以及丹可夫修正因子的守恒，将复杂几何燃料等效为规则的一维棒状燃料，在等效一维棒状燃料上利用伪核素子群计算燃料的有效自屏截面。数值结果表明：针对环板燃料栅元、蜂窝煤状燃料栅元、环板燃料组件、蜂窝煤状燃料组件和JRR-3全堆芯问题，该方法具有较高的计算精度与计算效率。