薛志恒，刘志宏，赵 晶，佘 顶

(清华大学 核能与新能源技术研究院，先进核能技术协同创新中心，先进反应堆工程与安全教育部重点实验室，北京 100084)

我国在役运行的第2代压水堆，其在线监测功能由核仪表系统[1]外核仪表连续监测反应堆功率变化及功率分布，监测反应堆径向功率倾斜和轴向功率偏差等。目前，使用堆外探测数据进行堆芯功率精细重构的方法主要有传输矩阵法、谐波综合法及神经网络法。

传输矩阵法认为堆外探测器信号与功率分布之间存在线性关系，堆芯功率可由探测器示数矩阵计算得到。这种做法可重构出堆芯轴向功率分布，但精确度尚有提升空间。清华大学核能与新能源技术研究院李富等[2]提出了谐波综合法。谐波综合法建立了各因素对功率分布变化影响的特征分布，但未考虑工况变化下截面的变化，可能导致无法准确模拟真实通量。夏虹等[3]提出利用RBF神经网络实现对堆芯功率三维分布的实时监测，其有统计数学上的意义，但忽略了扩散方程以及有效中子增殖因数等物理约束。

本文提出一种基于堆外探测实现堆芯功率重构的方法，并通过算例来验证该方法的可用性。本方法对引起反应堆堆芯功率波动的影响因素进行建模，建立影响因素的扰动与截面变化、堆外探测器示数之间的关系，利用特征统计算法(CSA)搜索扰动模型参数并通过扩散方程验证参数的准确性。

1 堆芯功率分布重构方法

贾健等[4]在研究基于堆内监测系统的堆芯功率重构方法时，提出了一种将功率偏差归结于组件等效截面参数偏差的功率分布重构方法。该方法的实质是将功率分布重构问题转化为以截面参数为自变量、堆芯功率分布为目标函数的优化搜索问题。本文将上述方法作进一步的改进，基于某供热堆堆型探讨利用堆外探测器提供的数据来进行堆芯功率分布的重构方法可行性。

该方法首先分析反应堆堆芯功率波动的原因，并对各影响因素进行建模；建立影响因素的扰动与截面变化、堆外探测器示数之间的关系。当堆芯功率波动时，调整影响因素模型参数，由上述关系得到堆外探测器计数，与堆外探测器的实际计数进行比较，如果不满足偏差要求，则重新调整影响因素的模型参数，直到满足偏差要求。

堆外测量系统需由堆内测量系统定期校正。本文假定堆外测量系统已由堆内测量系统定期校正过，在分析功率分布变化原因时，可忽略富集度偏差、设备加工精度等因素，只分析有限的几个因素，如温度场波动、氙振荡和控制棒移动等。本文仅分析氙毒模型和温度场波动模型。

氙浓度的变化会导致堆芯局部中子注量率的变化[5]。基于模型参数个数的要求和氙振荡的实际现象，给出下式：

ΔXe(x,y,z)=(a1x+a2y+a3z)·

exp(-a4x2-a5y2-a6z2)

(1)

式中：ΔXe(x,y,z)为氙毒在空间(x,y,z)点处的变化量；a1、a2、a3、a4、a5、a6为系数，是未知量。

堆芯冷却剂流经堆芯带走核燃料释放出的能量，温度场如果出现异常波动，可能的原因有燃耗变化、控制棒移动、环路流量失衡、堆芯入口温度变化及入口流量分配、包壳破裂等原因。借助谐波综合法的思想，使用特征分布的组合来表达温度场波动情况。本文为说明问题，用下式表达：

ΔT(x,y,z)=b1x+b2y+

b3z+b4xy+b5xz+b6yz

(2)

式中，b1～b6为温度扰动模型的系数。

进行扩散方程求解前需得到核截面数据。求解扰动变化量与截面变化之间的关系，常见的方法有线性插值法、切比雪夫多项式插值[6]等。本文采用线性插值法，公式如下：

ΔΣ1=c1ΔXe

(3)

ΔΣ2=c2ΔT

(4)

式中：ΔΣ为截面参数的实际变化量，cm-1；ΔXe为空间某位置氙毒变化百分比，100%；ΔT为空间某位置温度变化值，℃；c1、c2为常数系数，可通过计算得到。

2 优化问题的求解

本文选用特征统计算法[7]作为全局优化算法。该算法的基本流程如下：1) 在定义域内随机给出若干组解N0，并计算各自的目标函数值；2) 将得到的目标函数值排序，选取其中最好的N1个目标函数值，并选出其对应的N1组解；3) 描绘出每个特征量在定义域上的分布曲线；4) 根据分布曲线继续选取N0组解，计算出各自的目标函数值后，在N0+N1中选取最好的N1个目标函数值及其对应的解；5) 重复第3步及第4步，直至循环结束。

在Qt平台上使用C++语言对特征统计算法进行编程。本文涉及到的目标函数值是由扰动多项式的系数作为特征量组成的，目标函数为：

(5)

式中：Ni为第i组解；n为探测器能提供的数据个数；Pk0为第k个探测器提供的数据；Pki为通过计算得到的数据。F(Ni)值越小，表示重构的精确度越高。

3 计算测试模型

将堆芯各处的裂变中子对堆外探测器读数的贡献率作为堆外探测器响应函数，能准确描述堆内功率分布和堆外探测器读数的关系。实际的反应堆响应矩阵计算非常复杂。本文为了说明问题，进行了简化，忽略慢化剂温度变化、控制棒移动等对响应矩阵的影响，将其定为常数矩阵。程序中设定每个探测器仅受同一象限的堆芯影响，该象限对其他象限的探测器没有贡献。

本文采用某供热堆为模型，全堆共有148个组件，结构如图1、2所示。在扩散方程计算时，在堆芯x轴(数字方向)、y轴(字母方向)上按组件划分了14个节块；在z轴上从下到上(冷却剂流动方向)平均划分了16个节块。堆外共有4组功率量程探测器，分别位于4个象限内，每组由6个敏感段构成。

探测器响应矩阵为四维矩阵，每个元素表示堆芯某象限空间节块对该象限探测器某敏感段的响应。现将响应矩阵分为两个二维矩阵表示，第1个矩阵为堆芯轴向各层节块对6个探测器敏感段的响应矩阵(表1)，第2个矩阵为同一层内各节块对该层的贡献矩阵(表2)。

图1 某供热堆堆外探测器布置Fig.1 Out-core detector location of nuclear heating reactor

图2 某供热堆堆外探测器结构Fig.2 Out-core detector structure of nuclear heating reactor

4 计算测试结果

用C++语言在Qt平台上将中子扩散程序与特征统计算法耦合。特征统计算法提供特征值及计算特征量的收敛情况，中子扩散程序计算目标函数值，两者耦合可重构出堆芯功率分布。在扰动因素产生之前，各组件功率积分归一化后的功率分布示于图3。

4.1 温度场波动

表1 各层节块对探测器的响应Table 1 Response of each layer of node to detector

表2 各节块对该层的贡献Table 2 Contribution of each node to this layer

图3 某供热堆稳态功率分布Fig.3 Power distribution of nuclear heating reactor

在耦合程序中，限定特征量的搜索范围为[-1.0，1.0]，同时在计算截面变化时，设定截面必须为正，同时截面相对变化绝对值不超过100%。温度场波动模型如下：

ΔT(x,y,z)=0.1x+0.2y+0.3z+

0.1xy+0.1xz+0.1yz

(6)

经耦合程序计算后，得到最优的重构结果列于表3，对应的扰动模型多项式系数列于表4。

表3 温度场波动重构结果对比Table 3 Comparison of reconstruction results of temperature field fluctuation

表4 温度场波动模型系数对比Table 4 Comparison of coefficients of temperature field fluctuation model

表3中，目标函数值F可通过堆外探测器提供数据由式(5)计算得到，表征整体重构结果的精确程度；ΔPremax表示各节块中最大的相对功率偏差，其可表示搜索方法精细功率重构的可靠程度。表4中，Set表示温度扰动模型的设定系数值，Cal表示重构后计算得到的扰动模型系数值，Cal系数的精确程度由系数偏差向量的2-范数‖ΔX‖2即式(7)表示，该值代表搜索方法找出功率扰动原因的可靠程度。重构计算后，F为7.97×10-5，ΔPremax由6.72%降至0.067%，‖ΔX‖2为8.75×10-3，可说明重构结果准确且可靠。

(7)

4.2 氙振荡

氙振荡常发生在堆芯中间部位，该处的功率波动较大，探测器对此处的氙振荡不够敏感。为提高计算精度，在程序中合理减小特征统计算法的搜索范围，核截面数据将根据实际变化率来设定核截面的变化范围。随机设定氙毒波动的模型如下：

ΔXe(x,y,z)=(-0.2x+0.3y+0.15z)·

exp(-0.03x2-0.05y2-0.07z2)

(8)

得到的重构结果列于表5、6。

表5 氙振荡重构结果对比Table 5 Comparison of reconstruction results of xenon oscillation

表5中，F由4.27×10-3降至7.28×10-6，ΔPremax从1.401%降至0.053%，重构结果精确。表6中重构后的系数向量偏差2-范数‖ΔX‖2为4.09×10-2，重构结果可靠。

表6 氙振荡模型系数对比Table 6 Comparison of coefficients of xenon oscillation model

4.3 氙毒-温度场波动

随机设定双影响因素的模型如下：

ΔT(x,y,z)=0.1x-0.2y+

0.3z+0.1xy+0.2xz-0.15yz

(9)

ΔXe(x,y,z)=(-0.2x+0.3y+0.15z)·

exp(-0.03x2-0.05y2-0.07z2)

(10)

确定每个参数的定义域后，进行重构计算。重构结果列于表7、8。表7中，F由2.65×10-2降至8.42×10-5，ΔPremax由6.78%降至0.08%，重构效果较好。表8中，重构后的系数向量偏差2-范数‖ΔX‖2为0.118，偏差略大，这是因为温度扰动与氙振荡同时发生时，两个因素对功率变化同时起作用，程序无法区分两个变量对功率的影响，这是所有基于堆外探测数据进行功率重构的方法共同面对的问题。

表7 氙毒-温度场波动重构结果对比Table 7 Comparison of reconstruction results of xenon and temperature field fluctuation

沿着堆芯中轴线的切面、垂直于堆芯中轴线从下到上第1、5、8、12、16层节块的重构前、后功率相对偏差数据分别示于图4、5。

表8 氙毒-温度场波动模型系数对比Table 8 Comparison of coefficients of xenon and temperature field fluctuation model

图4 堆芯中心轴向切面功率重构前、后相对偏差分布Fig.4 z-axis distributions of power relative deviation at core center before and after reconstruction

图5 z轴第1、5、8、12、16层节块的功率相对偏差分布Fig.5 Distribution of power relative deviations in sections 1, 5, 8, 12 and 16 of z-axis

5 结论与展望

本文研究了基于堆外测量数据进行功率重构的扰动搜索方法，通过构造氙振荡及温度场波动单因素算例及双因素算例，初步验证了方法的可行性。

单因素算例可重构出堆芯功率分布，其最大节块功率偏差可控制在0.1%以下，同时可精确搜索出影响因素模型的系数，该功能可用来提供堆芯功率扰动的具体原因分析；双因素算例的堆芯功率最大节块偏差降到0.1%以下，但扰动模型系数重构没有单因素算例准确，是因为双因素相互干扰，这是该类问题的固有性质。

本文为说明问题进行了大量简化。该方法如在实际工程中进行应用，还要进行以下几方面的研究与开发：

1) 扰动模型的合理性与可区分性，如果扰动之间有因果关系，可进行合并分析；

2) 截面数据与扰动之间的关系不仅是常系数线性关系，可继续研究两者之间的敏感性；

3) 堆内功率分布与堆外探测器之间的响应矩阵不是常数矩阵，该矩阵的影响因素可展开分析；

4) 在扰动模型系数较多时，堆外测量可能存在多解，特征统计算法对多解的甄别能力需进行研究；

5) 目前在役压水堆的堆外测量系统仅用来实现堆芯轴向功率分布的测量，如将该方法进行大型压水堆堆芯功率重构，还需进一步研究，结合堆芯测量系统的堆芯出口温度测量系统进行堆芯功率的三维立体重构。