基于隐式重启Arnoldi方法的中子扩散本征值问题求解及其降阶研究

向钊才 陈洽锋 赵鹏程 张庆航

（南华大学 核科学技术学院 衡阳 421001）

多群中子扩散方程求解是核工程领域中的一个重要问题，通过求解中子扩散方程，可获得反应堆内中子分布、功率分布以及临界状态等信息，对研究反应堆物理特性、优化反应堆设计和提高反应堆的安全性具有重要意义［1］。然而，实际反应堆工况复杂多变，传统源迭代方法求解收敛速度慢，计算耗时长，难以在短时间内准确预测反应堆内中子注量率分布，因此有必要发展高保真、高效率的模型以便于重构反应堆内中子注量率分布。

模型降阶技术可以将高维度方程转化为低维度方程，降低计算复杂度并提高计算效率。20世纪70年代首次出现在工程领域文献中，经过几十年的发展，模型降阶已成为提高数值计算效率的有效方法［2］。隐式重启Arnoldi方法（Implicitly Restarted Arnoldi Method，IRAM）和本征正交分解方法（Proper Orthogonal Decomposition，POD）都是一种正交投影的模型降阶方法。近年来，利用隐式重启Arnoldi方法求解中子扩散方程的高阶谐波已有一定研究。黄义超等［3］基于隐式重启Arnoldi方法对中子扩散方程高阶谐波进行求解，结果显示IRAM方法在求解速度上相较传统源修正法提升了10倍，同时保持了相同的精度；王常辉等［4］利用高阶谐波的完备性和正交性重构堆芯中子注量率分布，并对反应堆功率进行在线监测，监测误差在±3%以内；谢金森等［5］基于重启Arnoldi方法求解了瞬发α本征值问题，验证了IRAM在求解瞬发α本征值问题的正确性。

POD法最早由Pearson于1901年提出［6］，后由Hotelling于1933年进一步发展［7］，1987年，Sirovich引入Snapshots方法求解POD基，极大地简化了数据处理过程，促使POD降阶方法得以快速发展［8］。关于POD法在中子扩散方程求解中的应用有了一些进展。Buchan等［9］采用基于时间依赖特征函数解的实例构造快照矩阵（Snapshot）并求解中子扩散方程的本征值问题；Sartori等［10］比较了模态法和POD法在求解时空动力学问题中的保真度，结果表明POD法具有更高的保真度；羊俊合等［11］将POD法与Galerkin投影结合求解中子输运方程，实现了较高的保真度和较高的计算效率，并且可重构堆芯中子注量率分布。然而POD方法对样本数据的模拟参数较为敏感，当样本输入数据质量较差时，输出结果的准确性就会变差，因此需要高精度模型以提供样本数据［12］。

采用隐式重启Arnoldi方法求解本征值问题的中子扩散方程得到不同宏观截面状态下的高阶谐波数据构建快照矩阵，基于POD与Galerkin投影相结合的方法构建POD-Galerkin低阶模型，并重构稳态TWIGL基准题中子注量率分布进行验证。该方法将利用高保真、高效率的IRAM生成快照矩阵并与POD-Galerkin降阶技术相结合，以提高计算效率和预测准确性，为堆芯中子注量率重构提供一种有效方法。

1 理论模型

1.1 中子扩散方程

中子扩散方程是中子输运方程扩散近似的结果，描述了单位体积内中子数守恒过程。多群中子扩散方程组是对中子扩散方程能量变量的离散，即将中子能量分成若干群。非临界状态下稳态表达式如下：

式中：ϕg是第g群中子注量率，cm-2·s-1；Dg为第g群中子扩散系数，cm；∑t，g是第g群中子移出截面，cm-1；∑s，g'→g是第g'群中子散射入g群的散射截面，cm-1；∑f，g是第g群中子裂变截面，cm-1；v为相应的裂变产额；χg为中子能量对应于能群g的能谱；keff为有效增值因子。

多群中子扩散方程用矩阵表示如下：

式中：L表示泄漏和移出矩阵；S是散射源矩阵；F是裂变源矩阵；是本征值。

通过式（2）可知中子扩散方程是一个广义本征值问题，经过转化可变成标准本征值问题：

式中：k是本征值。对于含增值介质的问题，k有n个正的离散解k1>k2>…>kn>0，对应谐波ϕ1，ϕ2，…，ϕn，其中n为谐波阶数。最大本征值k1是基波本征值，它就是有效增值系数keff；ϕ1是k基波，表示系统临界时中子注量率分布。而当系统偏离临界较远时，ϕ1并不表示系统真实的中子注量率分布。

1.2 隐式重启Arnoldi方法（IRAM）

针对式（3）的特征值问题，Arnoldi方法选取一个初始单位向量x0构造m（m>n）维克雷洛夫（Krylov）子空间：

再通过修正的格拉姆-施密特正交化（Modified Gram-Schmidt）求得该子空间的一组正交基Vm=[v1，v2，…，vm]，其中：

利用这组正交基将矩阵A约化为m维的上海森堡矩阵（Heisenberg）即：

通过求解Hm的本征值k和本征向量y即可得到A的本征值k和本征向量（谐波）ϕ：

Arnoldi方法将大型矩阵本征值问题投影到低维度Krylov子空间，通过求解低维度矩阵的本征值问题反演求得高维度矩阵的本征值和本征向量。然而Arnoldi方法需要储存Krylov子空间所有基向量，而且由于计算机舍入误差的积累可能导致基向量正交性的丧失，导致该方法计算效率和稳定性降低。隐式重启Arnoldi方法由Sorensen等提出［13］，该方法将Arnoldi过程与隐式位移迭代QR分解算法相结合，通过周期性地重启迭代过程来减小Krylov子空间的规模，并不断更新初始向量x0，从而保证基向量的正交性并降低计算成本。隐式重启Arnoldi方法的具体过程如下：

1）进行m步Arnoldi迭代，得到AVm=VmHm+；

2）求得Hm的前q个较小的特征值作为选取位移u1，u2，…，up；

3）QR迭代是一种数值计算方法，它通过迭代地进行正交三角分解（QR分解），将原始矩阵逐渐变换为对角矩阵或块对角矩阵的形式。对Hm进行q次隐式QR位移，得到正交矩阵Q，过程如下：

①令H1=Hm

②对i=1，…，q执行以下循环：

③得到最终正交矩阵：Q=Q1Q2…Qp

5）利用重启向量和截短形式回到第一步，重启整个过程直到收敛。

由于隐式重启Arnoldi方法在Arnoldi迭代过程中进行简单的自动处理，因此开发相关软件十分容易，目前已有许多程序包可以使用［14］。

1.3 POD-Galerkin法构建低阶模型

POD方法可以提取数据集中的主要模态或主成分，将高维数据进行降维处理的同时保留数据主要特征，通过POD方法，可以将中子扩散方程的高维度解降低到低维度解，从而提高计算效率和减少计算资源的需求。

通过IRAM可以求得不同宏观截面下的特征值ki和其对应谐波ϕi，利用这组谐波便可以构造一个快照（snapshot）矩阵M：

式中：L为快照数；ϕ1～ϕL为不同宏观截面下谐波。

利用式（9）的快照矩阵构造相关矩阵S：

式中：S为L阶方阵。

为获得快照矩阵中数据的特征信息，需求得相关矩阵S的特征值和特性向量。相关矩阵的特征值代表了数据在特定方向上的重要性，特征值越大，对数据贡献越大；特征向量代表了数据在特定方向上的分布情况，其方向表示数据在该方向上的主要特征，是一个相互正交向量，可构成POD正交基，用于描述数据的主要特征。求解公式如下：

定向越野运动不但考验学生的身体素质，对其综合能力也是一种考验。如，在开展定向越野运动时，首先，教师可多设定几条不同路线，并在不同路线上设置多个不同的目标，路线的选择与最终的成绩相关联，在这个思考、分析和选择的过程中，学生的综合能力也会得到相应的锻炼与提高。其次，教师可在地图上设置一些陷阱目标，从而锻炼学生的分析能力、选择能力和决断能力，这对学生未来的发展具有积极意义。

式中：λi和ψi为相关矩阵S的特征值和特征向量；φi为POD基。可以验证POD基满足正交性，即：

为提高计算效率，只需选取占比最重的前r个POD正交基，所取的POD基的个数即为模型的阶数，选取准则如下：

式中：r为POD正交基的个数，一般γ≥99，W为模态矩阵。

经过上述过程求得的POD基矩阵用最少的模态数描述快照矩阵的主要特征，反应数据特征信息的同时，极大提高了计算效率。中子扩散方程谐波可用这组正交基近似表示：

式中：bi为POD基φi(i=1，2，…，r)所对应的系数。只需求系数矩阵B便可求得中子扩散方程的谐波。

将式（16）代入式（3），再运用Galerkin法便可构建出中子扩散方程低阶模型：

从式（18）可知，POD-Galerkin法将高维度矩阵A本征值问题转化为低维度矩阵WTAW的本征值问题。通过求解低阶系统（式（18））的本征值和本征向量即可反演出全阶系统（式（3））的本征值和本征向量。

为验证低阶模型的准确度，需要构建谐波计算精度验证模型，采用王常辉等［4］提到的表征本征值和谐波计算精度的误差量：

采用有限差分法离散中子扩散方程，利用Matlab实现隐式重启Arnoldi方法求解本征值问题中子扩散方程获得谐波数据，并利用所得数据编制POD-Galerkin低阶模型代码。

2 二维稳态TWIGL基准题验证

2.1 TWIGL基准题介绍

TWIGL基准题是一个简化的矩形反应堆堆芯结构，包含三个区域，其中1区和2区是增值区，3区是包裹区，瞬态工况下1区的材料的热吸收截面发生变化来驱动反应堆功率变化，几何布置和边界条件如图1所示，各区域的物性参数如表1所示。

表1 二维稳态TWIGL基准题物性参数Table 1 Physical parameters of two-dimensional steadystate TWIGL benchmark problem

图1 TWIGL基准题区域划分、几何尺寸以及边界条件Fig.1 TWIGL benchmark problem area division, geometric dimensions, and boundary conditions

2.2 利用IRAM获取快照矩阵

由于TWIGL基准瞬态工况中，区域1中热中子吸收截∑a，2面发生变化，且变化范围在0.146～0.15 cm-1，因此，选择这一范围内的状态变化构造快照矩阵。利用IRAM方法分别求得∑a，2取值为0.148 cm-1和0.149 cm-1的前6阶本征值和前6阶谐波。IRAM本征值求解结果及误差如表2所示，前3阶谐波云图如图2所示。

表2 不同截面下前6阶本征值计算结果及误差Table 2 Calculation results and relative deviations of the first six eigenvalues for different cross-sections

图2 IRAM求解1区吸收截面为0.148时前3阶谐波快、热群分布（彩图见网络版）Fig.2 Distribution of the first three harmonics for fast and thermal groups when the absorption cross-section of Zone 1 is 0.148, solved by IRAM (color online)

从表2可知，IRAM求解εi误差数量级在10-14，具有很高的精确度，利用IRAM所求得的12组谐波构建快照矩阵M，通过式（10）～（15）求得一组POD基，该组POD基对应特征值与能量占比如表3所示。

表3 POD基对应的特征值及能量占比Table 3 Eigenvalues and energy ratios corresponding to POD basis

从表3可知，前6个POD基总能量占比高达99.99%，且每个POD基能量占比相近，因此选择前6个POD基构成模态矩阵W，它们用最少的模态数表示了快照矩阵的数值特征，这极大地减少了计算资源的消耗。

2.3 中子注量率分布重构

通过IRAM求解1区热中子吸收截面分别为1.48 cm-1和1.49 cm-1时的12组谐波，利用这组谐波构建快照矩阵并提取POD基，基于POD基的总能量占比构建模态矩阵W，最后再运用IRAM求解中子扩散低阶模型即可重构1区热中子吸收截面为1.50 cm-1时的本征值和谐波。所得结果如表4所示。

表4 TWIGL基准题堆芯前6阶本征值计算结果Table 4 Calculation results of the first six eigenvalues for the TWIGL benchmark problem core

从表4可知，采用POD-Galerkin方法构建的低阶模型在重构TWIGL基准题本征值和谐波方面，其误差εi随着本征值阶数而增加，但对于基波本征值和基波的预测该方法仍具有较好的精度。而实际反应堆往往在临界附近，这时基波本征值即为有效增值系数keff数值为0.913127，与参考解0.913214的误差仅为8.7×10-5。基波即为反应堆即为中子注量率分布，因此可重构堆芯中子注量率分布。重构中子注量率分布云图如图3所示，对角线上归一化中子注量分布与参考解［15］对比如图4所示。

图3 TWIGL基准题堆芯中子注量率分布 (a) 快群，(b) 热群（彩图见网络版）Fig.3 TWIGL benchmark problem core neutron flux distribution for fast group (a) and thermal group (b) (color online)

图4 对角线中子注量率分布与参考解对比 (a) 快群，(b) 热群Fig.4 Diagonal neutron dose rate distribution and comparison with reference solution for fast group (a) and thermal group (b)

从图3、4可以看出，降阶模型重构堆芯中子注量率分布与参考解基本吻合，对角线上快群和热群中子注量率最大相对误差为2.56%，在计算时间上，低阶模型计算一次的时间为0.039928 s，全阶模型计算一次的时间0.392965 s，低阶模型计算时间仅为全阶模型的10.18%，因此，低阶模型不仅可以准确重构堆芯中子注量率，而且可以有效减小计算用时，降低计算资源的消耗。

3 结语

采用隐式重启Arnoldi方法求解本征值问题中子扩散方程获得谐波数据，利用所得数据构建PODGalerkin低阶模型，并重构二维稳态TWIGL基准题中子注量率分布，所得结果如下：隐式重启Arnoldi方法在求解中子扩散方程的高阶本征值和高阶谐波时表现出较高的精度；基于POD-Galerkin低阶模型重构稳态TWIGL基准题中子注量率分布，在计算精度上，低阶模型预测对角线上快群和热群中子注量率最大相对误差为2.56%；在计算时间上，低阶模型计算时间仅为全阶模型的10.18%，因此低阶模型能在保持精度的同时有效降低计算资源的消耗。

本研究为堆芯中子注量率求解提供了一种可靠且高效的方法，该方法不仅可用于重构稳态时堆芯中子注量率分布，还具有在瞬态情况下预测中子注量率的潜力，有望在未来的应用中进一步拓展。

作者贡献声明向钊才负责起草文章，编写代码，分析/解释数据；陈洽锋负责采集数据，统计分析；赵鹏程负责对文章的知识性内容作批评性审阅，获取研究经费，指导；张庆航负责做表画图。