三维耦合分析软件的落棒事故分析

贺青云，陈 俊，马忠英，彭思涛，厉井钢

(中广核研究院有限公司，广东 深圳 518026)

落棒事故是指由单一的电气故障或机械故障引起的单个控制棒或整个控制棒子组件落入堆芯的事故，属于二类事故工况[1-2]。1个或多个控制棒落入堆芯会引入负反应性，导致堆芯功率降低。若落棒引发堆外探测器量测的负中子通量变化率达到停堆信号的整定值，反应堆停堆。若没有触发停堆信号，堆芯功率下降及一、二回路功率之间的不平衡导致堆芯入口温度下降。落棒事故可能会导致堆芯功率畸变，同时影响一、二回路控制系统的响应和事故的进程，是反应堆安全分析中最重要的事故之一[3]。然而，传统的落棒事故分析方法步骤繁多，计算量大，分析中考虑了过多的保守性，虽然计算结果满足DNBR准则，但热工裕量较低[4]。近年来，随着数值模拟工具和硬件的发展，耦合分析软件逐渐发展成熟，逐步具备工程应用条件，三维耦合分析软件作为数值反应堆的核心组成部分，将打破反应堆工程研发模拟的流程瓶颈。

自2000年以来，国内外利用耦合软件对落棒分析开展了诸多研究，如Ramos等[5]利用耦合软件pCTF-PARCSv2.7开展了控制棒下落瞬态分析，并将模拟结果与真实试验测量结果进行了对比；Mesado等[6]利用TRACE/PARCS耦合软件对三维落棒事故开展模拟和对热工水力变量进行不确定性和敏感性分析；Marina等[7]开发了RELAP5/PARCSv2.7耦合软件，对压水堆中子和热工水力参数局部偏差的响应进行了分析，并与真实电厂数据和SIMULATE-3K模拟结果进行了对比验证；冯英杰等[8]基于MANTA/SMART软件建立了三维物理热工耦合的落棒事故分析方法，与CPR落棒分析方法进行了对比分析。本文采用三维耦合分析软件，针对华龙一号开展不同棒组的落棒事故瞬态分析，相关结果将为核电安全性和经济性提供参考，也为数值反应堆研发提供一定的研究基础和技术支撑。

1 三维耦合软件

1.1 耦合方法

由中广核研究院有限公司自主开发的耦合软件，已经过与国际基准题TMI-1蒸汽管道破裂事故的对比验证，其中的核设计分析软件[9]和系统瞬态分析软件[10]分别经过了大量的例题验证，并应用在商用压水堆的分析中。耦合软件GINKGO/COCO[11-12]能获得较为精细的堆芯物理和热工水力分布，获取落棒事故发生后堆芯内部的功率分布和温度分布。三维耦合分析软件的框架图如图1所示，一维系统分析软件作为主控制端，采用动态链接库DLL[13]耦合技术，控制三维物理核设计软件的计算、函数与参数传递及计算状态监控等。

图1 三维耦合分析软件的框架图

耦合计算是在统一的时间步长Δt内进行，系统分析软件将Tn时刻的热工参数通过动态链接函数传递给物理软件，等待三维物理软件的计算，直到整个内部迭代计算收敛后，将Tn+1时刻的三维功率分布传递给热工软件，随后通过热工软件计算整个系统热工参数，以此交替进行完成整个瞬态计算。在耦合计算中，考虑到热工软件计算速度较快，为减少整个耦合的计算时间，在1个耦合时间步长Δt内采用多步热工计算。对于耦合时间步的控制，采用的是物理计算的反应性对其进行自动调整，以保证瞬态计算的收敛且不消耗过多的计算资源。

不同工况的最佳应用分析耦合软件不同：GINKGO/POPLAR，模拟反应堆的运行工况(模拟时间尺度量级为天量级)，其分析工具为一维耦合软件，模拟的时间步通常为秒量级，主要用于设计和优化I&C；GINKGO/COCO，一定条件下的瞬态计算，可用于燃料至包壳的相互作用分析，隔离阀误开，超载荷运行和满功率下的无控制落棒，反应性引入事故等，其分析工具为三维耦合软件；COCO/LINDEN，限制于堆芯内部的瞬态工况，如弹棒事故的模拟；GINKGO/COCO/LINDEN，如蒸汽管道破裂事故，未能紧急停堆的预计瞬变工况。本文用到的耦合软件有一维耦合分析软件和三维耦合分析软件，前者为保守分析工具，仅作为对比，后者为本文中的三维落棒分析用耦合软件。

1.2 基准题验证结果

OECD/NEA已发布了一系列的基准题用于反应堆时空动力学与热工耦合问题的计算准确性验证[14-15]。其中MSLB基准题使用三维的中子动力学堆芯模型进一步验证耦合软件的计算能力，适合于分析耦合堆芯和系统相互作用下的复杂瞬态工况和测试热工耦合问题。该基准题还具备由于环路不对称的冷却剂引起堆芯显著的时空效应，且瞬态事故中假设卡棒。三维物理热工耦合软件已经过国际基准题TMI-1蒸汽发生器事故算例的对比验证[11]。堆芯稳态平均轴向功率和瞬态堆芯平均裂变功率分别与国际基准题参与者平均值的对比结果如图2所示，可看出，三维耦合软件的计算结果与基准题平均值较为接近，详细误差数据可见文献[16]。

a——堆芯平均轴向功率分布；b——堆芯瞬态平均裂变功率

2 落棒事故与计算结果

2.1 落棒事故特征

当发生落棒事故后，反应堆自动保护是由负中子通量变化率保护系统提供，该系统特征如下：1)在落棒事故发生后的最初几秒中，核功率迅速降低；2)保护系统为2/4逻辑，假设单一故障，当通道的第三高负信号达到负中子通量变化率整定值时便会触发停堆；3)堆外探测器的响应依赖于落棒引入的反应性价值大小和落棒棒组在堆芯所处的位置；4)更多的保护通过其他不同的可预知途径来提醒操纵员反应堆运行处于不正常状态。

若堆芯保护不充分，对于不会触发中子通量变化率停堆信号的落棒瞬态，可能因为瞬态进程中的功率畸变和高功率而导致在某些工况下发生偏离泡核沸腾。落棒事故属于二类事故工况，必须保证其后果不能超过燃料设计限制，即瞬态中最小DNBR必须在限值之上。落棒事故分析分为两个阶段：1)探测分析阶段，为了确定在哪些落棒组合和条件情况下将触发反应堆跳堆信号；2)热工水力瞬态分析阶段，针对前一阶段中确定的未能触发跳堆的落棒情况进行落棒后的热工水力模拟分析并计算瞬态全过程中最小的DNBR，以验证燃料设计限值总能得到满足。

2.2 落棒算例

同一子组中一组棒、两组棒、三组棒或四组棒的下落均可构成落棒事故。落棒事故的分类需要考虑同一子组控制棒组的对称性，堆芯布置和落棒组合分类示意图如图3所示，堆芯布置的最外围为反射层，内部为燃料区域控制棒棒组，分为停堆棒(SA，SB，SC，SD)、温度调节棒(R)、功率调节棒(G1，G2，N1，N2)。其中G棒为灰棒，N和S棒为黑棒。

a——堆芯布置；b——落棒组合

本文选取3种不同的落棒方式，且考虑落棒后不停堆工况进行计算，3个算例分别为：Case1，一维落棒事故工况计算，初始工况反应堆处于首循环，燃耗初期，初始AO为+9%，控制棒落棒价值为307 pcm；Case2，两束棒组落棒事故工况计算，落棒编号为L05+E05，具体位置如图3a所示，落棒为停堆棒SD棒，初始工况反应堆处于首循环，燃耗初期，初始AO为+9%，落棒价值为334 pcm，R棒的初始棒位为46步；Case3，两束棒组事故工况落棒，落棒编号为K06+F06，具体位置如图3a所示，落棒为温度调节棒R棒。除落棒价值为269 pcm不同外，其余反应堆初始状态同Case2。Case1计算采用的棒组无法与图3b中的所列棒组合一一对应，所选算例落棒价值为一假定值，计算方式为一维耦合分析软件，Case2和Case3两个算例采用三维耦合软件GINKGO/COCO计算。

2.3 落棒事故计算

1)Case1算例计算

Case1采用一维耦合软件进行分析，选取了控制棒引入不同的反应性等多个状态点，并对系统主要参数做了保守性假设。耦合分析软件计算的堆芯核功率和二次侧汽轮机功率随时间变化的对应关系如图4a所示。由于核功率与汽轮机功率不匹配会导致R棒的运动，进而会导致核功率在60 s内与汽轮机功率不能达到一个相对的稳定值。图4b为一维落棒分析3个回路热管段的温度对比，热管道1和热管道3曲线基本重合，30 s后3个环路热管道的温度基本一致。

图4 一维落棒分析核功率与汽轮机功率(a)和环路温度(b)

利用CHF关系式W3公式[17]计算的DNBR变化如图5所示，结合图4a功率匹配曲线，汽轮机功率与核功率的不匹配，落棒会导致堆芯核功率的总下降，R棒的自动调节会使得功率回升，DNBR的变化呈现与功率相反的趋势。

图5 DNBR随时间的变化

2)Case2算例计算

耦合软件GINKGO/COCO计算的核功率和汽轮机功率对比如图6a所示。由于控制棒的下落，堆芯功率会明显下降，在约3.0 s时控制棒会落到底部。R棒控制系统反馈的来源有3列信号控制，第1列信号是堆芯的平均温度信号反馈，第2列信号是汽轮机的功率反馈，第3列信号是堆芯的核功率。R棒的动作如图6b所示，由于控制棒下落，堆芯功率下降，使得功率与汽轮机功率不一致，会使得R棒有动作，但R棒的动作由于机械控制的速度限制，不会立刻引起堆芯的功率调整。由于汽轮机功率与核功率的不匹配使控制棒会往上提升，结合图6a中两功率匹配曲线，当汽轮机功率与核功率达到匹配1∶1时，R棒不再动作。

图6 L05+E05两组落棒耦合计算功率(a)和R棒棒位(b)结果

反应堆堆芯反应性随时间的变化如图7a所示，落棒会对堆芯引入负反应性，而R棒的自动调节会使堆芯反应性上升，在20 s左右堆芯进入超临界状态，此时会由于反馈机制使R棒的动作进入微调状态，即R棒的调节速度变慢，存在R棒停止运动的状态，最终进入核功率和汽轮机功率相对稳定的状态。图7b为利用W3公式计算DNBR随时间的变化，随功率的下降，DNBR会上升，最终由于功率回升，DNBR逐步下降。

图7 L05+E05两组落棒的反应性(a)和DNBR(b)随时间的变化

因为非中心对称的控制棒组落棒，使得周围的组件径向功率峰下降，非落棒区域的功率峰上升，功率中心非对称图如图8a所示。L05+E05两组控制棒的下落会导致功率峰往落棒中心对称区域侧偏移，呈现出V字形。最大功率峰落在距离落棒最远的中间线上，最高值为1.472。这个位置靠近下侧的R棒位置。因为R棒的提棒动作，会在堆芯区域引入正的反应性，从而导致堆芯功率正向偏移平均功率。后续可通过子通道分析更准确地分析最小DNBR的位置。图8b为L05+E05两组落棒工况下3个冷管段环路进入堆芯下腔室入口的温度随时间的变化。可看出，60 s后3个环路的温度不同。环路1的温度最低，环路3最高，这与一维落棒计算环路温度变化有显著差异。

图8 L05+E05两组落棒事故径向功率分布(a)和环路温度(b)的变化

3)Case3算例计算

Case3算例计算的堆芯2D径向功率分布如图9a所示，整体趋势与Case2算例落棒工况相同，最大功率峰与平均功率的比值最大值为1.564 2，高于Case2计算的1.427 2。Case3算例计算的3个环路进入堆芯下腔室入口的温度随时间的变化如图9b所示。3个环路的温度变化趋势一致，但值有差异，环路1的温度最低，环路3最高，这与Case2算例结果趋势一致，但环路温度有所不同，60 s后Case3计算的3个环路温度分别比Case2计算的高。

图9 K06+F06两组落棒事故径向功率分布(a)和环路温度(b)的变化

4)Case2与Case3算例结果对比

图10a为两组落棒的核功率随时间的变化对比。Case3(K06+F06)的落棒价值低于Case2(L05+E05)的落棒棒组，可看出，Case3算例计算的功率回升值高于Case2，前者的调节速度稍晚于后者，最终两者都会达到相对稳定的状态。R棒的调节范围对比示于图10b，可看出，Case3棒组的调节范围大于Case2的计算。

图10 L05+E05和K06+F06两组落棒事故核功率(a)和R棒棒位(b)对比

3 结论

基于开发的GINKGO/COCO耦合软件，开展了停堆棒组落棒和R棒棒组两组落棒计算对堆芯功率分布的影响研究，对比分析了两组落棒在R棒自动调节后的稳定功率。通过3组不同落棒分析，可观察到非中心对称的棒组落棒事故会导致堆芯功率出现中心不对称，中心不对称的落棒棒组会导致堆芯功率出现中心不对称，从而使得堆芯出口环路温度不同。R棒棒组落棒价值虽小，但调整幅度大，落棒价值越大，调节后的稳态功率回升相比初始稳态差异越大。落棒事故算例W3公式计算的DNBR变化趋势与功率呈现相反规律，后续可结合更加精细的子通道耦合分析开展局部精细DNBR计算。