压水堆核电站堆芯物理/热工水力耦合特性研究

郑 勇，彭敏俊，夏庚磊，刘新凯

(哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001)

反应堆堆芯是一涉及反应堆物理、热工水力等诸多方面的复杂系统。一方面，堆芯内裂变功率的分布很大程度上决定了热工水力参数的变化。同时，热工水力参数的分布又会通过慢化剂温度效应和燃料多普勒效应来影响反应性的变化。RELAP5是研究反应堆热工水力现象的最佳估计系统程序，其早期版本的中子物理计算采用点堆模型，因此为了研究堆芯的某些不对称现象或局部具有强烈物理/热工水力反馈作用[1-3]时，国外学者将RELAP5程序作为热工水力程序通过耦合接口与其他物理程序进行耦合计算来研究反应堆堆芯多物理场耦合现象[4-6]。

RELAP5-3D是RELAP5系列程序的最新版本，取得的最主要进展是引入了三维物理计算程序Nestle和三维水力学部件，具备对堆芯进行三维物理/热工水力耦合计算的能力。RELAP5-3D程序经美国GSE公司改造，实现了在仿真平台SimExec上运行的目的，最终形成了RELAP5-HD程序。本文针对秦山核电二期工程堆芯部分进行仿真研究，基于两步法首先利用RELAP5-HD程序建立堆芯活性区的物理/热工水力耦合模型；其次由组件程序HELIOS计算出各类组件在大量工况下的少群截面参数，进行拟合处理后作为RELAP5-HD程序的输入值；最后进行稳态计算和事故工况仿真研究，以验证所建立的耦合模型的准确性。

1 RELAP5-HD模型建立

1.1 研究对象

秦山核电二期工程一号机组反应堆堆芯由121个燃料组件组成。每个燃料组件含有264根燃料棒、24个导向管和1个测量管。堆芯活性段高度为365.8 cm。为展平堆芯功率分布和实现分区换料，第1循环堆芯燃料按照235U富集度不同分3区装载。根据富集度和含有可燃毒物棒不同，本文将堆芯组件分为7组(表1)。

1.2 中子物理模型

RELAP5-HD的堆芯物理计算模型主要基于北卡罗来纳州州立大学开发的Nestle程序，并将原点堆模型中的衰变热计算模型进行了适当修改加入到其中，该程序的基本模型是三维带6组缓发中子的少群(2群或4群)扩散方程，利用成熟的节块展开法(NEM)和非线性迭代手段加速求解堆芯的三维功率分布，该程序可对堆芯稳态特征值问题和瞬态问题进行求解，且同时具有对正方形组件和六角形组件进行计算的能力[7]。

表1 需进行组件计算的燃料组件

RELAP5-HD中每个能群需输入的截面类型有5个(Σa、Σf、Σs、D、B2)，并提供5种截面反馈模型由用户选择。本文能群选择2群，截面反馈采用其中的GEN模型，在该模型中考虑了慢化剂温度、密度、硼浓度和燃料芯块温度对反应截面的影响。

(1)

(2)

C=min{1.0，max[0.0，hinsert/hl]}

(3)

根据选定的燃料组件类型和截面反馈模型，由HELIOS-1.11计算7种典型组件在平均燃耗点下，不同燃料温度、冷却剂密度、慢化剂温度和硼浓度时的反应截面，然后根据得到的多工况少群常数，选择其中1组工况作为基态截面利用Matlab7.1进行多元函数拟合，以获得7种典型组件的截面拟合系数。经检验，利用拟合公式求出的截面与原始截面的相对误差上限仅为0.3%。

1.3 物理/热工水力节块划分与耦合方式

根据节块法的求解需要，进行堆芯物理网格划分时，在径向上以1个组件为1个节块，而轴向上则考虑求解的稳定和快速收敛要求，尽可能以正方体为1个节块进行划分，最终将堆芯活性区划分为121×16=1 936个节块。堆芯热工水力节点在径向上以每个组件对应1个冷却剂通道，轴向网格划分考虑节点数量和计算特性的影响，将每个冷却剂通道沿高度划分成8个水力学控制体，外加顶部和底部各1个轴向反射层，反射层厚度取1个物理节块。如图1所示，在堆芯活性段每2个轴向物理节块映射1个热工水力节块，图中A～S为控制棒组件，围板/反射层组件的功率在稳态和瞬态时均设为定值，即1%满功率。在堆芯进出口使用时间相关控制体部件保证并联通道具有相同的压力边界，忽略环路对堆芯计算的影响，且不考虑相邻通道之间的相互交混。

a——径向节块划分；b——轴向节块划分

2 稳态结果分析

本文针对文献[8]的第1个燃料循环中燃耗为1 467 MW·d/tU的满功率工况点进行稳态计算。该工况点下主调节棒棒位为189步。由于硼在堆芯的分布较均匀，计算时通过调节堆芯进口硼浓度使堆芯达到临界。为便于与电厂的测量数据进行比较，给出的功率分布均进行了归一化处理。

图2为100%满功率工况下径向功率分布的比较。从图中相对误差(Num3)可看出，仿真结果和电厂参考值符合较好，只有极个别相对偏差达6%，误差来源主要包括实际测量的不确定性、反应截面计算误差、模型本身和入口条件的近似处理。组件计算中采用一维模型来产生围板/反射层等效均匀化常数(AXS-ADF)。由图2可见，径向反射层的存在，提高了边缘组件的中子利用率，对功率展平起到改善作用。

图2 100%FP工况下径向功率分布

图3为100%满功率工况下各组件通道的流量分配因子，文献[9]所进行的水力模拟实验结果也示于图3，该实验结果仅针对堆芯下腔室和流量分配孔板的几何结构对流量分配的影响进行研究，而未考虑堆芯活性区不同通道的积分功率水平不一致对流量分配的影响。本文在基于所有通道的几何结构一致的假设下，研究了堆芯功率分布对通道间流量分配的影响。由图3可看出，在不同的加热功率下，堆芯活性区的流量呈现明显的中心高、边缘低的分布特点。这主要是由于堆芯中心区域中子泄漏少，功率密度较边缘组件更高，而较高的燃料释热率使得燃料和冷却剂温度升高。一方面由于负温度系数的影响，较高的燃料和冷却剂温度会引入负反应性；另一方面，冷却剂温度升高而密度减小，因此流动阻力减少使该组件的流量增大。增大的流量会降低燃料和冷却剂的温度，对负反应性的引入起到抑制作用。堆芯的流量分配在这两方面因素的相互作用下达到平衡状态，呈现出高功率区大流量、低功率区小流量的特性。

图3 径向流量分配结果

表2为部分参数的测量值与仿真结果的比较。在本文计算工况下，临界硼浓度仿真计算值与电厂测量值的偏差极小，远小于一般的偏差允许值50 ppm。对于燃料芯块中心最高温度，文献[10]给出了COBRA-Ⅳ程序和COMBAT程序的计算结果，分别为1 937 ℃和1 895 ℃，与之相比，本文所计算的结果偏低，原因可能是所选燃耗工况点的功率峰与文献中不一致。除此之外，其他热工水力参数、径向功率峰因子和轴向功率峰因子均与实际测量值吻合得非常好。

表2 稳态计算结果与测量值的比较

表3为部分冷却剂通道出口温度计算结果。可看出，计算值和电厂测量值的绝对偏差较小，最大偏差为3.671 K(N09号通道)。这些数据可认为是图2、3结果的一个综合反映，表明流量和功率在堆芯的分布是合理的，与实际电厂测量值接近。

表3 冷却剂通道出口温度比较

3 掉棒事故工况结果

本文在100%功率稳态运行工况下，于290 s引入掉棒事故。假设掉棒事故发生后，主调节棒D的某束控制棒(图1中的L08)瞬间插入堆芯，并认为控制棒的插入不会对堆芯的几何结构和进口参数产生影响。

图4 掉棒事故下堆芯总功率的变化

图4为发生单束控制棒掉棒事故后反应堆堆芯总功率的变化。从图4可看出，由于掉落的控制棒瞬间引入了较大的负反应性，堆芯总功率急剧下降。较低的功率使得冷却剂和燃料温度相应下降，由于慢化剂温度效应和燃料多普勒效应，堆芯功率随后会呈现缓慢上升的趋势，最终稳定在一新的功率水平，整个功率波动过程持续35 s左右。

图5为掉棒事故后堆芯轴向和径向功率峰因子的变化。从图5可看出，掉棒事故会造成径向功率峰因子显著上升，而轴向功率峰因子轻微降低。这是由于掉棒事故后，控制棒的空间效应使径向功率分布发生畸变，径向归一化功率在掉棒的组件位置达到最小值0.392 1，而功率峰因子则达到1.387 5(图6)，在事故掉棒过程中甚至能达到1.43。在轴向上，由于单束控制棒掉落引起的空间效应并不显著，因此轴向功率峰因子变化较小。

图5 掉棒事故后功率峰因子的变化

图6为掉棒事故发生200 s后，堆芯各参数已基本稳定后的流量和功率分布情况。

图7为掉棒事故发生后部分热工水力参数的变化。对比图7可发现，掉棒事故对冷却剂流量分配因子的影响较小。在整个掉棒事故进程中，最大和最小流量分配因子的变化量均小于0.1。由于控制棒插入引起的功率突降使得最小流量因子出现在掉棒组件位置处，而多余的流量则分配到其他通道，所以引起的最大流量因子变化很小。

为保证第1道放射性屏障的完整性，须对燃料芯块温度给予特别关注。从图7b可知，事故掉棒发生后，堆芯最热组件的最高芯块中心温度会经历先下降后上升的过程，这种变化趋势主要是由该组件的功率变化引起的。而对于发生掉棒事故的组件，一方面由于燃料元件具有较大比热容不会使中心温度发生急剧下降；另一方面由于冷却该组件的冷却剂流量在掉棒后会降低，缓解了对燃料元件的冷却效果，因此其中心最高温度会持续缓慢下降，约在掉棒事故发生200 s后达到稳定。综上分析可预测，发生单束控制棒掉棒事故除了使堆芯参数产生波动和分布畸变外，不会使堆芯燃料组件的完整性受到破坏。

图6 掉棒事故后的稳态功率和流量分布

图7 掉棒事故后部分热工水力参数的变化

4 结论

为准确描述反应堆堆芯发生的物理/热工水力耦合现象，本文使用RELAP5-HD程序，对秦山核电二期工程反应堆堆芯在平均燃耗为1 467 MW·d/tU时的稳态和事故工况进行了研究，并将稳态计算结果与电厂实测值进行了比较，通过本文的研究可得出以下结论：

1) 耦合计算具有较高的计算精度，除个别组件外，稳态功率分布相对误差可达5%以下。

2) 堆芯冷却剂通道轴向积分功率不一致使冷却剂流量在堆芯中呈现中间高、边缘低的分布，且与堆芯下腔室几何结构相比，功率分布不均匀才是导致流量分布呈现上述特点的主要原因。

3) 发生单束控制棒掉棒事故会使堆芯功率发生畸变，特别是径向功率峰值会发生剧烈波动，并产生流量的再分配，但对流量分配因子峰值的影响不到0.1，远小于对功率分布的作用。

4) 通过进行单束控制棒掉棒事故工况分析表明，该事故虽然会造成功率和热工水力参数的波动，但不会对放射性屏障的完整性构成威胁，获得的事故曲线趋势合理。

参考文献：

[1] ANTONELLA L, WALTER A, ALESSANDRO P, et al. Analyses of pressure perturbation events in boiling water reactor[J]. Annals of Nuclear Energy, 2008, 35: 1 199-1 215.

[2] TEWFIK H, ANIS B, E’KHIDER S, et al. Application of coupled code technique to a safety analysis of a standard MTR research reactor[J]. Nucl Eng Des, 2009, 239: 2 104-2 118.

[3] D’AURIA F, SOLOVIEV S, MALOFEEV V, et al. The three-dimensional neutron kinetics coupled with thermal-hydraulic in RBMK accident analysis[J].Nucl Eng Des, 2008, 238: 1 002-1 025.

[4] 陈玉清，蔡琦，于雷. 三维瞬态物理分析模型与RELAP5/mod3.2程序的耦合研究[J]. 核动力工程，2010，31(3)：9-13.

CHEN Yuqing, CAI Qi, YU Lei. Study on coupling of NLSANMT and RELAP5/mod3.2[J]. Nucl Power Eng, 2010, 31(3): 9-13(in Chinese).

[5] 周铃岚，张虹，臧希年，等. 耦合核反馈并联通道异相振荡研究[J]. 核动力工程，2011，32(6)：66-70.

ZHOU Linglan, ZHANG Hong, ZANG Xinian, et al. Research on out-of-phase oscillation in a nuclear-coupled parallel double-channel boiling system[J]. Nucl Power Eng, 2011, 32(6): 66-70(in Chinese).

[6] 王军. RELAP5/MOD3与THAS-PC4程序耦合方法研究[D]. 北京：中国原子能科学研究院，2005.

[7] The RELAP5-3D Code Development Team. RELAP5-3D code manual, Volume 1: Code structure, system models, and solution methods[R]. Idaho, America: Idaho National Laboratory, 2005.

[8] 赵强. 核电厂反应堆堆芯物理在线仿真系统研究[D]. 哈尔滨：哈尔滨工程大学，2006.

[9] 杨来生，宗桂芳，胡俊. 秦山核电二期工程反应堆水力模拟实验研究[J]. 核动力工程，2003，24(增刊)：208-226.

YANG Laisheng, ZONG Guifang, HU Jun. Reactor hydraulic simulation test study of Qinshan Phase Ⅱ NPP project[J]. Nucl Power Eng, 2003, 24(Suppl.): 208-226(in Chinese).

[10] 李经纬. 秦山核电二期工程反应堆热工水力设计[J]. 核动力工程，1999，20(4)：308-312.

LI Jingwei. Reactor thermal-hydraulic design for Qinshan Phase Ⅱ Nuclear Power project[J]. Nucl Power Eng, 1999, 20(4): 308-312(in Chinese).