基于MANTA/SMART三维物理热工耦合的落棒事故分析

冯英杰，李昌莹，肖 红，胡友森

(1.中广核研究院有限公司，广东 深圳 518031；2.环境保护部 核与辐射安全中心，北京 100082)

事故分析是核电厂安全分析的重要组成部分，它研究核电厂在故障工况下的行为，是核电厂设计过程的重要步骤[1]。落棒事故是指由单一的电气故障或机械故障引起的某一控制棒子组中任意数目的控制棒或整个控制棒子组落入堆芯的事故，属于设计基准二类工况(DBC-2)[2-3]。

落棒事故是典型的二类反应性事故，畸变的堆芯功率和一、二回路控制系统的响应对事故的进程都有很大影响[4]。因此分析方法中需将这两方面的特征都刻画出来。CPR1000核电厂是中国广核集团在法国M310堆型的基础上，通过技术改进而提出的自主技术品牌[5]。传统的CPR分析方法采用保守的确定论分析方法和部分统计法，设计保守性较大。为提升机组经济性、挖掘热工裕量，需考虑重新开发新的落棒分析方法[6]。

本研究拟提出一套基于MANTA/SMART的三维物理热工耦合的落棒事故分析方法，实现探测阶段和瞬态分析阶段的三维模拟，更真实地模拟落棒事故后反应堆堆芯行为和一、二回路及控制系统的响应，从而实现为落棒事故释放更大的设计裕量的目的。

1 研究方法与假设

1.1 技术方案总体描述

MANTA程序是法国AREVA公司开发的用于模拟压水堆正常运行和非破口类事故瞬态的分析程序。该程序由4部分组成，包括热工水力模块，燃料热力分析模块，堆芯物理模块，以及核蒸汽供应系统、堆芯保护与控制和仪控模块。

SMART程序是一个三维两群堆芯扩散-燃耗计算程序，是商用核电设计软件SCIENCE程序包的一部分[7]，由法国AREVA公司开发。在APOLLO2-F数据库的支持下，使用SMART程序可进行稳态、瞬态工况下全堆芯三维双群中子扩散-燃耗计算，得到整个组件和组件内棒中子通量以及功率分布。

本研究基于MANTA/SMART程序提出的三维物理热工耦合的落棒事故分析方法，能更好地反映堆芯功率畸变和一、二回路控制系统的响应。图1为MANTA/SMART耦合的数据流示意图。MANTA程序承担热工水力计算以及控制、保护系统的模拟，SMART程序承担堆芯中子学的计算，二者在约定的时间步上交换物理热工参数、控制棒棒位、探测器信号等数据，从而实现计算的耦合。

图1 物理热工耦合的数据流示意图Fig.1 Data flow diagram of neutronics/thermal-hydraulics coupling

三维落棒事故分析方法主要分为3个步骤，总体分析流程如图2所示。

图2 总体分析流程Fig.2 General analysis process

1) 数据准备阶段。利用三维物理程序，计算后续分析需使用的核数据。

2) 探测阶段。通过三维物理热工耦合模拟堆外中子探测器的响应，根据负中子通量变化率确定未能触发停堆信号的落棒工况。为得到足够保守的探测结果，需构建1个堆芯模型，惩罚影响事故后果的关键参数，包括缓发中子份额、多普勒功率系数和慢化剂温度系数等，同时需对轴向功率偏移(AO)和初始棒位进行敏感性分析。

3) 瞬态分析阶段。对于未触发停堆的落棒组合，利用三维物理热工耦合模拟堆芯中子学参数及热工水力参数在瞬态过程中的变化。同样地，为得到保守的偏离泡核沸腾比(DNBR)结果，需构建堆芯模型惩罚关键参数，并在此基础上开展瞬态模拟。最后，论证在任意时刻反应堆状态均满足安全分析的设计准则要求。

传统CPR分析方法采用三维物理加一维热工的方法，步骤较多，计算工作量较大，分析中考虑了过多的保守性，虽然计算结果满足DNBR准则，但热工裕量较低[8]；法国EPR落棒事故分析方法与CPR分析方法较类似，本质上是三维中子学加零维热工的方法，未考虑物理和热工参数之间的耦合。美国AP1000[9]落棒事故分析方法的思路与上述两种方法的不同，该方法没有直接计算最小DNBR，而是通过验证事故下的径向功率焓升因子(FΔH)低于限值来保证事故下DNBR不超出限值，但本质上它也是三维中子加零维热工的分析方法，未考虑物理和热工参数之间的耦合。法国AREVA公司开发的先进落棒事故分析方法，采用了三维物理热工耦合的分析工具，在分析思路上与CPR分析方法有相似之处，也是通过一维计算画出探测曲线，并根据特征参数筛选后续分析的工况。

本文提出的落棒事故分析方法，采用纯三维的分析思路，能释放出关键参数的保守性，进而挖掘出安全裕量。另外，在探测阶段对所有落棒组合都进行了三维模拟，分析逻辑更严密，有利于后续批量化自动化计算的实现。

1.2 分析假设

三维落棒分析方法所使用的程序有SMART、MANTA和FLICA Ⅲ-F。其中，SMART用于中子学数据计算，在三维耦合中承担物理计算部分；MANTA用于热工水力计算，在三维耦合中承担热工水力计算部分；FLICA Ⅲ-F用于子通道分析[10-11]。

在分析中，假设初始的功率、冷却剂平均温度以及稳压器压力皆取名义值，这些参数的不确定性在DNBR限值中考虑；另外，堆芯流量取热工水力设计流量(3台主泵运行工况)；同时，考虑不同的轴向功率分布以及温度控制棒组的位置。

反应堆通过堆外负中子通量变化率提供保护。若落棒引入的反应性足够大，4台堆外中子探测器中有2台的信号超过高负中子通量变化率阈值，则会触发反应堆停堆。考虑堆外探测器的单一故障准则，即4台堆外中子探测器中有3台的信号超过高负中子通量变化率阈值，才认为该落棒工况触发反应堆停堆。

除平均温度控制系统外，其余的反应堆控制系统假设不适用[12-13]。落棒引入负反应性导致堆芯功率下降，而温度控制棒组的抽出会加剧落棒后堆芯功率的增加[14]。

事故分析考虑的参数不确定性或设计裕量列于表1。对于探测阶段，原则是使落棒组合更难被探测到。对于瞬态分析阶段，原则是使堆芯热工水力状态更恶劣，瞬态过程中的DNBR更小。

表1 事故分析考虑的参数不确定性或设计裕量Table 1 Parameter uncertainty or design margin considered in accident analysis

2 计算结果与分析

以某大型压水堆核电厂为分析对象，采用本文提出的基于MANTA/SMART三维物理热工耦合的方法论，进行落棒事故分析。

该堆型反应堆堆芯有177盒燃料组件，其中控制棒燃料组件68盒，由于堆芯布置旋转对称，根据其特性，控制棒可分为17个子组，每组4束。

堆芯物理热工水力模型的节点划分方案如图3所示。对于中子学模型，在径向方向，每个组件划分为4个节点；在轴向方向，堆芯活性区划分为16层。对于热工水力模型，在径向方向，每个组件作为1个热工水力通道；轴向方向划分为16层，与中子学模型相对应。

图3 堆芯物理热工水力模型节点划分Fig.3 Nodalization of netronics/thermal-hydraulics model

2.1 探测阶段

对寿期初氙平衡(BLX)、寿期中(MOL)和寿期末(EOL)所有落棒组合进行模拟，判断每个工况是否触发负中子变化率高停堆信号，筛选出不可探测工况。根据分析结果，BLX未探测工况有34个，MOL未探测工况有49个，EOL未探测工况有54个。以K02+B06+P10+F14子组(BLX)为例，该子组中各落棒组合的堆外探测器功率量程中子通量变化率变化曲线示于图4。工况K02表示始发事件为控制棒组K02落入堆芯，工况K02B06表示始发事件为控制棒组K02和B06同时落入堆芯，以此类推。由图4可看出，工况K02、K02B06和K02F14的TILT3(绝对值第3大的负中子通量变化率)在落棒过程中均小于负中子通量变化率停堆阈值，所以没有触发高负中子通量变化率停堆信号；而工况K02B06P10和K02B06P10F14则超出停堆阈值，因此触发停堆。

传统CPR分析方法使用三维SMART计算所有工况的径向再分布因子，然后采用一维程序计算几个典型工况下TILT3的临界值，画出探测曲线，从而筛选为触发停堆的工况[15]。

图4 堆外探测器功率量程中子通量变化率变化曲线Fig.4 Neutron flux rate of excore detector

而在三维分析方法中，则针对每一落棒工况，计算4个堆外探测器负中子变化率随时间的变化，从而判断在落棒过程是否达到停堆阈值。这使得探测阶段的分析去除了一定的保守性。

2.2 瞬态分析阶段

利用MANTA/SMART程序，针对未探测到的落棒工况进行瞬态模拟，并使用FLICA Ⅲ-F进行子通道分析，得到每个落棒组合的最小DNBR。通过比较各工况的计算结果，最终得到的各燃耗步的最小DNBR列于表2。该循环下最小DNBR为1.79，设计裕量为28.4%。传统分析方法下，落棒事故裕量一般在15%左右，某些换料设计中甚至会到10%以内。因此，采用三维分析方法能有效提高落棒事故设计裕量。值得注意的是，该方法下可得到具体落棒组合对应的DNBR，而传统CPR方法则是用包络的方式确定最恶劣工况并对其进行瞬态模拟。各燃耗步发生最小DNBR时刻的主要物理热工参数列于表2。

图5为L03E13工况整个事故过程中FΔH和AO随时间的变化曲线。控制棒的下落导致堆芯功率迅速发生畸变，因此FΔH急剧增大；堆芯功率减小和冷却剂温度降低，使得温度控制棒组往外抽出，对径向功率分布畸变有一定程度的缓和作用，FΔH减小并趋于稳定。同时，温度控制棒组的抽出也使得AO逐渐增大，轴向功率分布变得恶劣。

对于FΔH，在CPR分析方法中，先计算燃料管理方案中每个工况落棒前后的FΔH，从而得到所有落棒组合的FΔH变化率，进而描点画出包络线，在子通道分析中使用的是一个定值；而三维分析方法可精细刻画出事故进程中径向功率分布的变化，在子通道分析中FΔH使用的是动态值。从图5可看出，虽然落棒后FΔH较大，但此时功率较小，最小DNBR一般出现在事故后半程，因此从FΔH上能挖掘到裕量。

表2 最小DNBR及主要物理热工参数Table 2 Minimum DNBR and principle neutronics/thermal-hydraulics parameter

图5 L03E13工况FΔH和AO随时间的变化Fig.5 Curve of FΔH and AO for case L03E13

对于AO，由于CPR分析方法采用的是一维模型，子通道分析中轴向功率分布取堆芯平均值；而三维分析方法则取热组件对应的轴向功率分布进行计算，因此，从图5能看出三维方法的AO明显小于CPR分析方法的AO。

图6为最小DNBR时的堆芯功率分布，其中L03和E13是落棒位置，此处功率份额显著减小，约为0.43；而远离落棒位置的组件的功率份额则明显增大，约为1.45。可见，落棒造成堆芯径向功率分布发生明显畸变。在轴向上，能明显看出堆芯功率向上偏移，这是由于温度控制棒抽出所导致的。图7为使用两种分析方法计算得到的核功率变化趋势对比，CPR分析方法计算得到的功率回调高于三维分析方法。对CPR传统分析方法和MANTA/SMART三维分析方法进行归纳比较，结果列于表3。

图6 最小DNBR时的堆芯功率分布Fig.6 Power distribution at moment of minimum DNBR occurrence

图7 核功率变化趋势对比Fig.7 Comparison of nuclear power

表3 落棒事故分析方法对比Table 3 Comparison of control rod drop accident analysis methods

3 结论

本研究提出了一套基于MANTA/SMART三维物理热工耦合的落棒事故分析方法，并论述了该方法的原理、工具、流程、假设及基于该方法的分析结果。分析表明，该三维落棒分析方法能更真实地体现堆外探测器的探测逻辑和反映落棒事故瞬态过程中主要参数的变化，事件序列更贴近实际情况，变化符合预期，具有可行性。在三维模型中，落棒及平均温度控制棒动作所造成的径向和轴向功率分布畸变，较一维模型更真实。一方面，三维分析方法能模拟出事故过程中FΔH随时间的变化，用以代替CPR方法中的固定保守值；另一方面，采用热组件AO代替堆芯平均AO进行子通道分析。三维分析方法对事故的关键参数刻画得较精细，为落棒事故分析释放出一定的设计裕量。