解析函数展开节块法的数值稳定性分析

孙 伟，李 庆，王 侃

(1.清华大学 工程物理系，北京 100084；

2.中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610041)

方形几何节块法中，横向积分是其中的关键技术，但在六角形几何中直接采用横向积分技术时横向泄漏会出现奇异项，奇异项的出现给节块中子注量率的求解带来困难[1]。针对此问题，韩国KAIST的Cho等[2]提出了解析函数展开节块(AFEN)法。AFEN法是在六角形几何内不分离变量而直接利用严格满足中子扩散方程的解析基函数，把节块内的各群中子注量率近似展开。然后根据节块间耦合条件，将节块面平均净流表示为节块体平均中子注量率和面平均中子注量率的函数关系，代入节块中子平衡方程形成全堆关于节块体平均中子注量率和面平均中子注量率的方程组，分别用加速的高斯塞德尔方法和源迭代方法求解中子注量率和堆芯有效增殖因数keff。

研究发现利用AFEN法数值求解一些特殊问题时，如keff接近于某一节块的无限介质增殖因数kinf(节块在堆芯中实际净泄漏接近于零)时，面平均净流对中子注量率的关系系数出现震荡，使得关于中子注量率方程组的系数矩阵出现病态，进而方程组的求解不收敛，无法得出中子注量率的解。方形几何中，基于横向积分技术的非线性迭代解析节块法的这类不稳定性问题，Joo等[3]做了相关研究，对病态节块采用线性近似或节块展开法(NEM)代替解析节块法求解；而六角形几何中，针对中子注量率直接展开的AFEN法的这类不稳定性问题，尚未见相关研究，本文提出两种解决方法来解决上述问题，即截断近似方法和泰勒展开方法。

1 AFEN法数值不稳定性分析

1.1 AFEN法公式推导

对于反应堆内一个均匀的六角形节块，其三维多群稳态中子扩散方程[1]可写成如下形式：

(1)

其中：r为空间位置变量；V为节块体积；Φ(r)=[Φ1(r)，Φ2(r)，…，ΦG(r)]T为G群中子注量率构成的向量。矩阵Λ(keff)的元素Λgg′为：

νΣfg′)/Dg

(2)

式(1)的解析解主要取决于矩阵Λ(keff)的特征值λm及相应特征向量um的数值类型。对于一般多群问题，λm既可能为实数也可能为复数，本文仅介绍λm均为实数的情况。

利用Λ(keff)的λm及um可将每群中子注量率展开为：

(3)

其中：U为矩阵Λ(keff)所有特征向量构成的矩阵；第m列对应于特征向量um；μgm为矩阵U的元素；Aml、Bml、Am0为展开系数；el为节块展开时所选坐标轴上的单位向量，l取4时其具体表达式参见文献[1]。

SN、CS为双曲或三角函数：

(4)

利用净流和中子注量率的关系可得节块内面平均净流与面平均中子注量率、体平均中子注量率的关系：

(5)

将式(5)代入式(6)可得全堆芯所有节块关于面平均中子注量率、体平均中子注量率的方程组，可简写为式(7)，对式(7)的求解采用源迭代法，具体求解过程参见文献[4]。

,2,…,G

(6)

MX=N

(7)

式中：系数矩阵M与群截面、堆芯几何参数、组件不连续因子有关；X为堆芯所有面平均中子注量率构成的向量；N与体平均中子注量率有关。

1.2 AFEN法数值不稳定性原因探究

式(5)中系数矩阵F的每个元素都是双曲或三角函数的表达式，与组件半宽度H、λm有关。表1列出矩阵元素F(G，1)(具体表达式详见式(8)和(9)，其他元素的表达式在形式上与F(G，1)的相同)随λm的变化情况(H取7.35)。可看出，当λm>0时，随λm的减小，F(G，1)变化平缓，但λm减小到某一值时，F(G，1)出现畸变，在λm接近于0时F(G，1)变为无穷大；当λm<0时，在接近于0时F(G，1)也为无穷大，通过观察F(G，1)的求解式可发现，当λm接近于0时，F(G，1)求解式各项分子、分母均接近于0，当两个接近于0的数运算时会引入舍入误差[5]。矩阵F元素的畸变会引起式(7)中系数矩阵M的元素也相差几个数量级甚至为无穷大，导致M变为刚性矩阵，式(7)内迭代求解不收敛。

表1 F(G，1)值随特征值的变化情况

当λm>0时，有：

(8)

当λm<0时，有：

F(G,1)=(H|λm|(-1+cos(H|λm|0.5))-

H2|λm|1.5sin(H|λm|0.5))/6(-2+(2+

H2|λm|)cos(H|λm|0.5))+|λm|0.5×

(-1+cos(H|λm|0.5))/6(-H|λm|0.5×

cos(H|λm|0.5)+sin(H|λm|0.5))+

H|λm|sin(H|λm|0.5)/3(-H|λm|0.5×

cos(H|λm|0.5)+sin(H|λm|0.5))+

|λm|0.5(H|λm|0.5cos(H|λm|0.5)×

(-1+H|λm|0.5cos(0.5H|λm|0.5))-

sin(H|λm|0.5))/2(-2+H2|λm|+

(2+H2|λm|)cos(H|λm|0.5)-

H|λm|0.5sin(H|λm|0.5))

(9)

系数矩阵F的元素值出现畸变，是由于矩阵Λ(keff)出现接近于0的特征值，本文仅给出2群截面时的公式推导解释特征值接近于0的原因，2群时矩阵Λ(keff)的具体形式及截面间的转换关系如下：

Λ(keff)=

(10)

Σr1=Σt1-Σ1→1;Σr2=Σt2-Σ2→2;

χ1=1;χ2=0;Σ2→1=0

(11)

其中，Σr1、Σr2为1、2群移除截面。由式(10)、(11)可知矩阵Λ(keff)行列式为：

(12)

当Λ(keff)出现接近于0的特征值时，Λ(keff)行列式接近于0，得：

·

(13)

又根据文献[3]可知，2群截面时节块的kinf为：

·

(14)

由式(13)、(14)可知，当keff接近于某一节块的kinf时，矩阵Λ(keff)会出现接近于0的特征值，引起系数矩阵M的病态，出现不收敛情况。而kinf接近于keff的含义是节块在堆芯中净流为0，这多出现在堆芯中部，且堆芯含有多种富集度燃料，因为一般堆芯都有泄漏，净流很难为0。

1.3 AFEN法的数值稳定性方法介绍

出现不收敛的直接原因是矩阵Λ(keff)的特征值很小时，M变为刚性矩阵，导致迭代不收敛，为解决这一问题，提出2种方法。

1) 截断近似方法。从表1可知，在畸变前当特征值很小时，元素值基本不变，此时可将出现畸变的值截断，直接用畸变前一刻的值代替。具体方法是当0<λm<10-6时，λm近似取为10-6，当-10-8<λm<0时，λm近似取为-10-8(由于λm在(-0.01，0.01)范围内时元素值变化不大，且堆芯内出现特征值接近于0的节块极少，因此截断点的选择基本不影响堆芯计算结果的精度，本文截断点是参考F(G，1)变化趋势给出的)。截断近似的本质是对病态节块采取近似常解析基函数展开。

2) 泰勒展开方法。当-10-8<λm<10-6时，对矩阵F元素求解式中的双曲函数或三角函数进行泰勒展开，泰勒展开的目的是对求解式继续化简，避免出现分子、分母接近于0的项。表1列出5阶和7阶泰勒展开时元素F(G，1)的计算结果，可看出7阶展开求出的F元素值更接近于真实值。泰勒展开的核心思想是高阶多项式展开，与Joo等提出的对病态节块直接采取NEM方法相比，7阶泰勒展开保留原AFEN法的求解格式，在保证精度的同时更容易实现。

2 数值验证

HANDF-D[4]是基于解析函数展开节块法编制的可用于三维多群六角形几何计算的程序，其正确性已得到多个基准题的验证。为检验上述两种稳定性方法的有效性，将其添加到HANDF-D程序中。验证对象为改造后的三维VVER440基准题[1]，改造方法是减小VVER440材料1的两群吸收截面(表2)，使keff接近于材料2的kinf(1.082 55)。将TRIVAC[6]的计算结果作为基准解，TRIVAC程序是由加拿大开发的可计算三维多群六角形几何堆芯的中子学程序，采用了多种方法，包括有限差分方法、有限元方法等。

表2 VVER440材料1截面变化

在使用原HANDF-D程序计算改造后VVER440基准题(轴向分12层，每25 cm为1层)时，迭代中出现keff为1.082 56，接近于材料2的kinf(1.082 55)，矩阵Λ(keff)出现极小特征值1.429 2×10-7，迭代不收敛。表3列出采用两种稳定性方法时keff的结果，可看出两种方法计算均收敛且keff与TRIVAC的相比相对偏差仅为-0.02%。图1示出功率分布的比较结果，功率分布最大相对偏差为2.03%，出现在燃料与反射层交界处，这与文献[4]中指出的HANDF-D程序本身固有2.5%左右的功率偏差相当。从keff和功率分布结果可看出，两种稳定性方法计算结果精确可靠，均能很好地解决解析函数展开节块法计算不稳定性问题。

表3 改造后VVER440基准题keff计算结果

图1 改造后VVER440基准题功率分布比较

3 结论

通过计算分析发现利用AFEN法求解六角形堆芯中子注量率，当keff接近于某一节块kinf时，会使得此节块的矩阵Λ(keff)出现极小特征值。利用此极小特征值求解全堆芯关于中子注量率方程组的系数矩阵时，由于引入舍入误差使得系数矩阵变成刚性矩阵，导致迭代不收敛。针对此类问题，提出2种解决方法：截断近似方法和泰勒展开方法。对改造后的VVER440基准题验证表明，两种方法均可有效解决不收敛问题，计算结果和参考解符合很好这一事实也充分说明这两种方法有较高的精度。

参考文献：

[1] 夏榜样. 三维多群六角形节块法及其应用研究[R]. 西安：西安交通大学，2006.

[2] CHO N Z, NOH J M. Analytic function expansion nodal method for hexagonal geometry[J]. Nucl Sci Eng, 1995, 121(3): 245-253.

[3] JOO H G, JIANG Guobing, DOWNAR T J. Stabilization techniques for the nonlinear analytic nodal method[J]. Nucl Sci Eng, 1998, 130(1): 47-59.

[4] 孙伟，倪东洋，李庆，等. 三维多群六角形几何中子扩散程序开发[J]. 原子能科学技术，2013，47(10)：1 707-1 712.

SUN Wei, NI Dongyang, LI Qing, et al. Development of 3D multi-groups neutron diffusion code for hexagonal geometry[J]. Atomic Energy Science and Technology, 2013, 47(10): 1 707-1 712(in Chinese).

[5] 周铁，徐树方，张平文，等. 计算方法[M]. 北京：清华大学出版社，2006.

[6] HÉBERT A, SEKKI D. A user guide for TRIVAC version4[R]. Montreal: Insitut de Genie Nucleaire Department Degenie Mecanique-Ecole Polytechnique de Montreal, 2009.