大破口始发应急工况预测方法研究

王 宁，王任泽，杨亚鹏，冯宗洋，贾林胜，刘一宁，梁博宁

(中国辐射防护研究院，太原 030006)

福岛事故表明虽然反应堆严重事故发生概率极低，但是它仍然可能发生。福岛事故前，事故后果的预测和评价始终是核与辐射应急管理的一项重要内容，受到各国学术界的重视[1]。福岛事故后国际原子能机构(IAEA)提出，应直接根据堆芯损伤评价的结果指导防护行动[2]。在应急情况下依据几个关键的反应堆工况参数实时、快速地预测未来的反应堆工况与事故进程，可以为应急决策提供支持，为场内和场外防护行动的实施赢得更多的时间。

由于神经网络的方法具有计算速度快这一显著优点，所以国内外学者使用该方法进行了一些事故工况预测的研究。

韩国学者Seung Geun Kim等使用支持向量的计算器，利用核电站状态变量输入在短时间内的变化，预测严重事故中可能发生事件的发生时间[3]。韩国学者Soon Ho Park等使用模糊神经网络预测严重事故中的压力容器水位[4]，因为在严重事故中压力容器的水位无法测量。

韩国学者Dong Yeong Kim等使用人工智能的方法预测安全壳内氢气浓度的变化[5]，该研究将破口位置分为热管段、冷管段以及蒸汽发生器传热管三种。破口尺寸分为210个步长：对于冷却剂丧失事故(LOCA)，破口尺寸范围为双端断裂尺寸的1/10000～1；对于蒸汽发生器，传热管破裂数量从1根变化到210根。破口尺寸的误差为0.4%。预测安全壳氢气浓度所用的两个输入值为破口尺寸、停堆时长；使用MAAP4的模拟数据进行模型的开发和验证。

韩国学者Man Gyun Na,Sun Ho Shin等使用概率神经网络方法预测何时堆芯发生裸露、堆芯出口温度(CET)温度何时超过648.9 ℃(一般认为此时开始严重事故管理)、压力容器何时失效[6-7]。

清华大学核能与新能源技术研究院针对高温气冷堆示范电站，开发了一套核应急响应支持系统[8]。该系统有两个任务：一是基于测量仪器的异常进行故障诊断；二是基于测量仪器的读数和操作员的行动预测事故进程。为了更好地执行这两项任务，开发了动态贝叶斯网络，在该网络中使用贝叶斯机器学习(LBP)算法进行事故诊断和预测。

SESAME由IRSN开发，用于在LOCA或蒸汽发生器传热管破裂(SGTR)事故下诊断、预测核电站状态并分析源项，可以给出主要的事故特征参数，包括破口位置和尺寸、堆芯失水剩余时间、裂变产物释放到环境中的活度随时间的变化等。

日本学者Tsutomu Ishigami等为了实时预测压水堆安全壳失效时间，开发了一项分析技术[9]。预测安全壳失效时间的思路是：熔融堆芯熔穿压力容器后，冷却剂经压力容器下腔室泄漏进入安全壳使得安全壳压力升高，通过简单的方程，建立未来安全壳压力与当前安全壳压力的关系，从而预测未来的安全壳压力，得到安全壳失效时间。将得到的结果与THALES、STCP等成熟程序的模拟结果进行了对比，验证了该方法的准确性。

预测的关键在于超实时计算，上述研究通过神经网络、简化方程等方法实现了应急工况进程的预测。本文尝试使用经典公式预测应急工况进程，针对M310反应堆建立了大破口失水事故(LBLOCA)始发应急工况进程预测方法。根据事故特点将复杂的一回路进行合理简化，对质量守恒与能量守恒方程进行近似求解。将堆芯活性区划分成4个径向环和10个轴向层，即40个栅元。使用实验关联式计算堆芯换热，得到包壳温度，然后根据包壳温度可判断堆芯工况。基于本文建立的方法，可以开发大破口始发应急工况进程预测程序。

1 事故情景与假设

本文分析的事故情景为大破口失水叠加全厂断电(SBO)始发严重事故。建模使用的M310反应堆部分参数列于表1。

基于事故情景做如下假设：

(1)事故前反应堆处于100%功率稳定运行状态；

(2)大破口和全厂断电于零时刻同时发生；

(3)由于大破口发生后0.1 s内一回路压力即可降到冷却剂饱和压力，瞬间产生大量蒸汽，空泡效应引入的负反应性可使反应堆停堆，故假设0 s时反应堆停堆；

(4)高压安注、低压安注、主辅给水失效，非能动中压安注能够正常运行；

(5)主泵惰转时间为30 s；

(6)安全壳泄露速率为0.1%/d；

(7)仿真计算以下封头失效为终止事件，安全壳的响应不在本文的研究范围内；

(8)破口发生在稳压器所在环路。

2 应急工况预测方法

2.1 守恒方程与整体思路

流体的质量守恒和能量守恒方程是模型的基础，守恒方程的形式取决于所采用的两相流模型。本文使用基于两相分相流模型的守恒方程，如式(1)～(4)所示[10]。

表1 建模所用的部分参数Tab.1 Part of parameters of the model

①气相质量守恒(连续性方程)

(1)

②液相质量守恒方程

(2)

式中，ρl为液相密度，kg/m3；Ul为液相流速，m/s。

③气相能量守恒方程

(3)

④液相能量守恒方程

1.2.2 疼痛护理。疱疹可伴随不同程度或不同类型的疼痛感，如刀割样疼痛、烧灼样等，尤其是眼周疱疹，疼痛感极易影响患者身心健康。本报告中患者为刀割样疼痛，痛感强烈并对护理内容实施造成影响。因此，护理人员有必要给予疼痛护理。采用通俗易懂的语言将疱疹相关信息对其进行讲解，告知患者疼痛感是由于病毒侵犯神经导致，使患者对疾病的认知增强，提高对疼痛的耐受力。如痛感较强，可给予止痛药处理，有效减少痛感对患者造成的影响。

(4)

求解守恒方程需要初始条件、边界条件以及结构方程，对守恒方程离散成非线性方程组，进而迭代求解。如果根据事故特点将复杂的一回路进行合理简化，对守恒方程进行近似求解，使用经典公式近似估算应急工况进程，计算速度将大大加快，可以实时评价严重事故应急工况，也可以预测事故工况的发展。

通过泄漏流量、安注流量计算堆芯水位，然后进行堆芯换热计算，最后计算堆芯包壳温度，通过包壳温度判断应急工况。

2.2 一回路的简化

大破口失水事故具有如图1所示的事故特点[11]，本文根据这些特点，将大破口失水事故分为喷放、再灌水/再淹没、“长期冷却”三个阶段，针对每个阶段采用不同的计算方法与模型。

图1 大破口始发应急工况特点Fig.1 Characters of LBLCOA initiated emergency condition

在大破口失水事故期间，一回路迅速泄压，稳压器很快失效，故可忽略稳压器的作用；在喷放阶段，蒸汽发生器起到最终热阱的作用，但是喷放阶段持续时间很短，喷放阶段结束后，冷却剂流经蒸汽发生器U型管的流量很小，可以忽略蒸汽发生器的影响。为满足预测的计算速度要求，可将稳压器和蒸汽发生器忽略，反应堆一回路可等效简化为一个圆柱体，该圆柱体的体积等于一回路冷却剂所占体积；在该圆柱体冷管段所对应的高度上，冷却剂通过破口发生泄漏。模型示意图如图2所示。

图2 一回路简化示意图Fig.2 Simplification of the primary loop

2.3 堆芯换热计算

堆芯换热是极为重要的结构方程。有研究表明[12]，严重事故期间辐射换热量对包壳温度以及氢气产生量的贡献不大(如图3、图4所示)，所以在堆芯换热计算中忽略辐射换热。本文依据包壳温度判断堆芯工况，而应急工况下包壳温度和燃料芯块温度基本一样，所以可忽略燃料棒径向温度分布，不计算径向导热量。

使用式5计算两相邻节点之间的轴向导热量[13]:

(5)

式中，Φ为沿轴向通过两相邻栅元之间的热流量，W；λ为堆芯栅元的导热系数，W/(m·K)；A为栅元的横截面积，m2；Δt为相邻栅元的温度差，℃；Δl为相邻栅元中心点距离，m。

图3 氢气产生量对比[12]Fig.3 Comparison for hydrogen production[12]

图4 堆芯温度对比[12]Fig.4 Comparison for core temperature[12]

在喷放、再灌水以及再淹没阶段，堆芯冷却剂与燃料棒的对流方式属于强迫对流，使用式(6)～式(8)计算层流强迫对流换热系数，使用式(9)计算紊流强迫对流换热系数，取两者最大值作为强迫对流换热系数[14]:

Nu=C(n)gdev

(6)

(7)

(8)

Nu=0.023Re0.8Pr0.4

(9)

式中，C(n)为常数，对于棒束和单管均取4.36；(z-z0)为计算点距离流道入口的长度，m；Dh为流道的水力学直径，m；Re为雷诺数；Pr为普朗特数。

在事故的其他阶段，对流方式为自然对流，使用式(10)～(11)计算对流换热系数，同样取两式计算得到的最大值作为自然对流换热系数:

Nu=0.18Ra1/4(L/Dh)-1/9

(10)

Nu=0.065Ra1/3(L/Dh)-1/9

(11)

式中，L为流道长度，m；Ra为瑞利数。

2.4 堆芯节点划分

将堆芯划分为4个径向环和14个轴向层(其中10个轴向层位于堆芯活性区)，如图5所示。轴向层1代表二次支撑组件和堆芯支撑板，轴向层2代表流量分配孔板，轴向层3代表堆芯下栅格板，轴向层4～13表示堆芯活性区，轴向层14表示堆芯上栅格板。

图5 堆芯节点划分Fig.5 Node division in the core

2.5 应急工况判断依据

根据包壳温度可判断应急工况，判断依据列于表2。

表2 应急工况判断依据 [15]Tab.2 Judgment basis for emergency condition [15]

2.6 临界流

事故初期一回路内外压差很大，冷却剂泄漏出一回路管道后瞬间汽化，流动为两相临界流，物理现象比较复杂。使用Moody临界流计算泄漏流量。两相临界流动的质量流速是空泡份额的线性函数，如式(12)所示[14]:

(12)

式中，α为空泡份额；ρg为气相密度，kg/m3；ρl为液相密度，kg/m3；ρm为两相流密度，kg/m3，ρm=αρg+(1-α)ρl；GC(α)为两相临界流质量流速，kg/(m2·s)；GC(1)为气相临界流质量流速，kg/(m2·s)；GC(0)为液相临界流质量流速，kg/(m2·s)。

3 结论

本文针对M310反应堆，建立了大破口失水始发应急工况进程预测方法。控制方程基于两相分相流模型。求解守恒方程需要初始条件、边界条件以及结构方程，对守恒方程进行离散成非线性方程组，进而迭代求解。根据事故特点将复杂的一回路进行合理简化，对守恒方程进行近似求解。使用经典公式近似估算应急工况进程，计算速度大大加快，可以实时评价严重事故应急工况，也可以预测事故工况的发展。

通过泄漏流量、安注流量计算堆芯水位，然后进行堆芯换热计算，最后计算堆芯包壳温度，通过包壳温度判断应急工况。堆芯划分为4个径向环和14个轴向层。使用包壳温度判断应急工况。堆芯换热计算是极为重要的结构方程，使用实验关联式进行计算。

基于本文建立的应急工况预测方法，可以开发响应的大破口始发应急工况预测程序。针对其他应急工况的预测方法正在研究中。