用于pin-by-pin输运计算的SP3解析基函数展开节块方法研究

彭良辉，汤春桃，毕光文，张宏博，杨 波

(上海核工程研究设计院有限公司，上海 200233)

三维全堆芯pin-by-pin输运计算是下一代堆芯物理分析方法(NGM)的重要研究内容。相较于传统的组件-堆芯“两步法”，该方法可获得更加精确的燃料棒层级的通量及功率分布。三阶简化球谐函数(SP3)方法是实现三维全堆芯pin-by-pin输运计算的重要途径之一。SP3方法的精度高于扩散方法，其计算量约为扩散方法的2倍，且远小于其他输运算法，如特征线方法(MOC)、离散纵标方法(SN)等[1]。选取SP3方法作为实现三维全堆芯pin-by-pin输运计算的途径，可较好地平衡计算精度与计算量之间的矛盾。此外，SP3方法的控制方程为两个扩散类型的方程，可沿用成熟的扩散节块方法对其进行求解。

本文采用解析基函数展开节块(AFEN)方法[2-3]求解SP3方法的控制方程，并对3种偏流表达形式的计算精度进行检验。

1 理论模型

1.1 SP3方法的控制方程及边界条件

多群SP3方法的控制方程[4]如下：

(1)

(2)

以上方程中：

(3)

Σr0=Σt-Σs,Σr2=Σt+0.8Σr0

(4)

(5)

(6)

基于Marshak边界条件，可得u方向上各阶偏流表达形式[5]为：

(7)

(8)

为了便于采用原有的扩散节块展开方法(NEM)求解SP3方法的控制方程，美国普度大学研发了PARCS程序对偏流表达形式进行截断[6]，即有：

(9)

(10)

此外，在日本NFI研发的SCOPE2程序中，为简化节块内偏流及平均通量响应矩阵，便于程序并行化，采用了与菲克定律类似的偏流表达形式[1，7]：

(11)

对于中子能群g，式(7)～(11)中的偏流表达形式可统一改写为如下形式：

(12)

表1 偏流表达形式的系数矩阵

获得偏流表达形式后，在两相邻节块k和节块k+1的交界面处，满足如下的偏流连续条件：

(13)

(14)

1.2 通量近似解

多群SP3方法的控制方程可改写为如下矩阵形式[3]：

(15)

设(λm,um),m=1,…,2G为矩阵[A′]的特征值与特征向量。

um=[u1,m…ui,m…u2G,m]T

m=1,…,2G

(16)

则式(15)存在如下通量解析解：

φ=[U]ψ

(17)

其中，[U]=[u1…um…u2G]。ψ中的第m分量为：

m=1,…,2G

(18)

其中：el为任意单位向量；r为位置向量；Aml和Bml为待定系数。

图1示出单位特征向量及定解条件示意图，其中S表示面。如图1所示，选择如下的3个单位特征向量，对式(18)中的辅助函数Ψm进行3阶截断(即取l=1,2,3)：

e1=ex,e2=ey,e3=ez

(19)

图1 单位特征向量及定解条件示意图

此时，辅助函数Ψm的近似解表达式为：

Ψm=em·cm

(20)

其中：

(21)

cm=[Am,1Am,2Am,3Bm,1Bm,2Bm,3]T

m=1,…,2G

(22)

将式(20)代入式(17)有：

φ=[U]Ψ=[U][E]C

(23)

其中：

(24)

(25)

式(23)即为通量近似表达式。其中：矩阵[U]为式(15)系数矩阵[A′]对应的特征向量矩阵，矩阵维数为2G×2G；矩阵[E]的矩阵元为含空间位置变量的指数函数，矩阵维数为2G×12G；C为待定系数向量，维数为12G，即每个节块中共有12G个待定系数。

1.3 节块响应矩阵

在某一节块中，为求出所有通量待定系数，需提供12G个定解条件。如图1所示，可选取节块6个边界处的各阶出射或入射偏流为定解条件，反解待定系数。为此首先需建立偏流与待定系数之间的关系式。令：

(26)

结合式(23)和(26)有：

(27)

(28)

其中：

(29)

其中，3种偏流表达形式下的系数矩阵[Mn]如表1所列。

jout=[P]C

(30)

同理，以节块6个边界处的入射偏流为定解条件，可得：

jin=[Q]C

(31)

其中：

(32)

(33)

由式(30)和(31)可得节块内偏流的响应矩阵为：

jout=[P][Q]-1jin

(34)

对式(23)进行体积平均，即对[E]的矩阵元进行体积平均，可得：

φave=[T]C

(35)

结合式(31)可得节块内平均通量的响应矩阵为：

φave=[T][Q]-1jin

(36)

式(34)及(36)即为节块出射偏流、节块平均通量相对于节块入射偏流的响应关系。

1.4 源迭代方法

基于节块中出射偏流、平均通量与入射偏流间的响应关系，可构造如下的源迭代流程：1) 初始化系统keff和各节块边界处的入射偏流jin；2) 计算或更新各节块中的响应矩阵[P][Q]-1、[T][Q]-1；3) 根据式(34)，由入射偏流jin计算出射偏流jout；4) 根据式(36)，由入射偏流jin计算节块平均通量φave；5) 根据节块交界面连续条件和外边界条件，由出射偏流jout更新入射偏流jin；6) 由节块平均通量φave更新keff；7) 判断keff和节块平均通量φave是否收敛，若未收敛则返回第2步，若收敛则结束。

此外，为提高计算效率，可采用粗网有限差分加速方法对以上源迭代过程进行收敛加速[4]。

2 数值检验

采用式(7)～(11)对应的3种偏流表达形式，开发了用于pin-by-pin输运计算的SP3解析基函数展开节块程序AFEN_SP3，对该程序进行数值检验。

2.1 自定义一维基准例题

该例题由C5G7-MOX-2D基准例题[8]演化得到，其材料及几何布置如图2所示，从左到右依次为H2O、UO2、4.3%MOX、7.0%MOX和8.7%MOX。各材料区的宽度均为20 cm，左边界为真空边界，右边界为反射边界。各材料的7群反应截面均根据C5G7-MOX-2D基准例题由特征线输运程序PEACH[9]经通量-体积均匀化得到。

图2 自定义一维基准例题示意图

该例题的参考解由MCMG程序计算得到[4]。AFEN_SP3程序计算时，网格宽度为1 cm，计算结果列于表2。其中扩散程序为半解析扩散节块程序Stella[4]，网格划分方式与AFEN_SP3程序保持一致。特征值keff计算偏差为绝对偏差，栅元功率计算偏差为相对偏差，以下算例相同。可见，对于该一维例题，3种偏流表达形式下的SP3程序计算精度基本相同，且均优于扩散程序。

表2 自定义一维基准例题计算结果

2.2 自定义二维基准例题

该例题由C5G7-MOX-2D基准例题[8]演化得到，其材料及几何布置如图3所示，且几何尺寸及布置方式与C5G7-MOX-2D基准例题相同。各材料的7群反应截面均根据C5G7-MOX-2D基准例题[8]由特征线输运程序PEACH[9]经通量-体积均匀化得到。

该例题的参考解由特征线输运程序PEACH[9]计算得到。AFEN_SP3程序计算时，网格大小为1.26 cm×1.26 cm，计算结果列于表3。可见，对于该二维例题，3种偏流表达形式下的SP3程序计算精度基本相同，且均优于扩散程序。

图3 自定义二维基准例题示意图

表3 自定义二维基准例题计算结果

注：1) 各栅元功率相对偏差绝对值的平均值

2.3 三维KUCA基准例题

该例题为日本京都大学建立的KUCA(Kyoto University Critical Assembly)三维轻水堆基准例题，该基准例题的几何结构如图4所示，包含燃料区、反射层和控制棒，各区的截面参数列于表4[10]。

该例题的参考解由TEPFEM程序中的P3模型计算得到[11]。AFEN_SP3程序计算时，网格大小为1 cm×1 cm×1 cm，计算结果列于表5。由表5可见：SP3程序的计算精度明显优于扩散程序；3种偏流表达形式中，Marshak偏流表达形式计算精度最高，其余两种偏流表达形式的计算精度在可接受范围内。

如图4，取以线段L(x∈[0 cm，25 cm]，y=2.5 cm，z=0.5 cm)为中心的1排栅元，其裂变源归一化快群平均通量分布如图5所示。图5中，以Marshak偏流下通量计算结果为基准，进行相对计算偏差统计，用于比较各偏流形式下计算结果的差别。

图4 三维KUCA基准例题示意图

表5 三维KUCA基准例题计算结果

图5 三维KUCA基准例题典型区域快群通量分布

可见，堆芯内部(包括燃料区和控制棒区)各偏流表达形式下的计算结果差别小于1%，靠近堆芯外部计算相对偏差变大，最外围节块相对偏差最大，约为4%。结合式(13)及各偏流表达形式可见，内边界偏流连续条件均可由交界面处通量及通量一阶偏导数连续得到。根据式(14)，3种偏流表达形式间的差异主要体现在外边界真空条件之上。因此，当堆芯外围反射层区域较大时，如上述的一维及二维基准例题，3种偏流表达形式下的节块相对功率分布计算结果差别较小。

3 结论

基于SP3解析基函数展开节块方法，在3种不同偏流表达形式下开发了pin-by-pin输运程序AFEN_SP3。数值检验结果表明：AFEN_SP3程序的开发是正确的；3种偏流表达形式下的SP3程序计算精度均优于扩散程序；3种偏流表达形式中，Marshak偏流表达形式计算精度最高，其余两种偏流表达形式计算精度稍差，但仍在可接受范围内。简化后的偏流表达形式进一步简化了SP3方法，可采用成熟的扩散节块程序稍加改动进行求解。简化后的偏流表达形式有利于进一步提高SP3全堆pin-by-pin输运程序的计算效率，有利于程序并行化算法的设计。