大型先进压水堆IVR策略评估和裂变产物分布计算

2022-07-29夏兆阳周志伟石雪垚

原子能科学技术 2022年7期

夏兆阳，周志伟，石雪垚，常愿

(1.清华大学核能与新能源技术研究院，北京 100084；2.中国核电工程有限公司，北京 100840)

反应堆发生严重事故时，堆芯燃料损坏熔化并同时释放大量裂变产物到安全壳中[1]。在无有效缓解措施投入情况下，熔融材料将熔穿下封头，掉落到反应堆堆腔，与堆腔的混凝土相互反应(MCCI)，进一步释放裂变产物到安全壳中。已有学者[2-5]对熔融物熔穿下封头的严重事故中，反应堆内裂变产物的释放和迁移行为进行了研究和分析，其研究成果已用于严重事故下反应堆的源项控制和管理。

为避免熔融物掉落到堆腔内发生MCCI、燃料冷却剂相互作用(FCI)等反应而威胁安全壳的完整性，美国圣芭芭拉大学[6]提出了熔融物堆内滞留(IVR)策略。该策略目前已用于三代先进压水堆设计中，如AP1000[7]、华龙一号[8]、APR1400[9]等。大型先进压水堆的堆腔注水系统(CIS)[8]通过向压力容器外表面和保温层之间的流道注水来冷却下封头，带走堆芯熔融物衰变热，避免下封头失效，以实现IVR。中国核动力研究设计院以概率论的方式论证了能动CIS的有效性[10]。中国核电工程有限公司评估了全厂断电(SBO)叠加小破口事故下非能动CIS的有效性[11]。美国圣芭芭拉大学[6]参考AP600压力容器几何搭建ULPU-2000实验设备，模拟了压力容器外部冷却(ERVC)条件下下封头外表面热流密度的变化，并根据该装置的实验Ⅱ和Ⅲ结果拟合了下封头外表面CHF(临界热流密度)实验关系式[12]。之后，为更精确地得到AP1000堆型压力容器外表面冷却效果，圣芭芭拉大学的Esmaili等[13]将实验Ⅲ的保温层流道优化为环形流道，得到了适用于AP1000的外表面CHF实验Ⅴ结果。美国爱达华国家实验室[14]使用SCDAP/RELAP5-3D软件对APR1400在失水事故下的ERVC缓解能力进行了计算。法国原子能委员会的Tarabelli等[15]使用ASTEC V1.2软件评估了VVER-440/213的IVR策略。捷克核研究所的Duspiva等[16]也使用MELCOR 1.86软件对VVER-1000/320在大破口事故下ERVC的缓解能力进行了计算，但未得到长期稳定熔融池结果。上海交通大学的金越等[17]利用MELCOR 1.86软件对1 700 MW大功率反应堆在SBO叠加大破口(LLOCA)严重事故下的ERVC进行了评估，计算通过更改临界热流传热公式实现反应堆IVR策略。韩国核安全研究所的Lim等[9]利用严重事故软件MELCOR 2.1对APR1400反应堆在IVR-ERVC策略下的熔融池传热进行了计算，重点关注了事故进程、堆芯退化行为以及熔融池行为。

本文利用严重事故分析软件评估大型先进压水堆CIS对一回路冷管段双端断裂严重事故的缓解能力。计算分析能动CIS投入下先进压水堆事故进程、熔融池传热行为、堆芯裂变产物的释放和迁移行为，为事故管理和放射性防护提供参考。同时将CIS下结果与无CIS事故后果进行比较和分析，综合评估IVR策略对事故和源项控制的能力。

1 几何模型建立

大破口严重事故下，能动CIS从内置换料水箱(IRWST)取水，通过CIS泵注入到压力容器外表面和保温层之间的流道，以900 m3/h流量对压力容器进行长期冷却[10]。冷却水携带熔融池热量从通道流出，进入到安全壳大空间，经过热量排出系统(PCS)和壁面结构的冷却回到IRWST，以保证水箱有足量水源支持CIS长期取水。反应堆CIS的冷却循环流程示意图如图1所示。

图1 反应堆CIS示意图Fig.1 Schematic of reactor CIS

对先进压水堆的安全壳隔间和主系统回路进行控制体建模，包括17个安全壳内系统隔间(包含3个蒸汽发生器隔间、稳压器隔间、安全壳上部空间等主要隔间)以及完整的3个冷却剂回路和压力容器结构，如图2所示。整个堆芯结构被认为是中心轴对称结构，堆芯几何划分为6个径向环和23个轴向层。为研究下封头内熔融池特性，下封头几何占有8个轴向层，如图3所示。下封头壁面为半球形，划分为9段，每段壁面单独计算与内部熔融池和外部冷却水的热量交换。

图2 反应堆安全壳控制体划分Fig.2 Control volume model of reactor containment

图3 反应堆堆芯和下封头节点划分Fig.3 Nodalization of reactor core and lower head

2 理论模型

2.1 下封头内熔融池模型

严重事故发生后，反应堆堆芯燃料棒和结构首先熔化，失去完整形状，形成碎片向下掉落，部分通过流道掉落到下封头内，部分被支撑结构阻塞；当堆芯支撑结构坍塌时，上方支撑的碎片和坍塌材料全部掉落到下封头内，形成熔融池。对于下封头内的熔融池，主要有两层熔融池模型和三层熔融池模型[13,18]。在严重事故分析软件中，熔融池模型为两层熔融池模型，上层为密度较轻的金属熔融池(MP2)，主要由未氧化包壳材料锆(Zr)和支撑结构材料不锈钢组成；下层为密度较大的氧化物熔融池(MP1)，主要由熔化的堆芯燃料(UO2)和氧化锆(ZrO2)组成。而熔融池外侧是由熔融池材料冷凝形成的碎片硬壳(PD)，如图4所示。

图4 下封头内熔融池的结构示意图Fig.4 Schematic of molten pool in lower head

相邻下封头单元内的熔融池会对流混合，具有均匀的材料、放射性核素和相同的温度组成。独立的熔融池会通过对流换热将热量传递给周围的支撑结构、相邻熔融池和上表面接触流体，还将通过辐射换热将热量传递给位于上表面的堆芯热构件。软件中两层熔融池的热平衡方程如下：

(1)

QMP2→MP1-QMP2→fluid-QMP2→rad

(2)

压力容器下封头壁面接受来自内部的热量(包括固体颗粒碎片和基底材料的传热)、熔融池的直接传热及壁面节点之间的传热。下封头的热量qlh,out由外部所在控制体中的流体带走，该热量定义为：

qlh,out=hatm(1-FPL)Alh(Tlh,out-Tatm)+

hpool,lhFPLAlh(Tlh,out-Tsat,pool)

(3)

其中：FPL为控制体中水池的比例；Alh为下封头外表面面积，m2；Tlh,out为下封头外表面温度，K；Tatm和Tsat,pool分别为大气温度和水池饱和温度，K，计算时忽略水池沸腾前的传热量；hatm和hpool,lh分别为下封头向大气和水池的对流传热系数，W/(m2·K)。hatm默认设定为10 W/(m2·K)；hpool,lh默认采用热构件沸腾传热关系式计算，本文则采用简单的池内泡核沸腾传热关系式，计算公式如下：

hpool,lh=hNB=34.5p1/4ΔT1.523

ΔT<23.4 K

(4)

hpool,lh=1.41×107p1/4ΔT-2.575

ΔT≥23.4 K

(5)

其中：p为水池的压力；ΔT=Tlh,out-Tsat,pool为表面过热度。假设控制体内压力不变，沸腾传热系数随表面过热度的变化如图5所示。

图5 沸腾传热系数随过热度的变化Fig.5 Boiling heat transfer coefficient vs superheat temperature

池内沸腾传热分为泡核沸腾、稳定膜态沸腾和处于两者间的过渡沸腾传热。泡核沸腾的临界热流密度(CHF)是下封头壁面角度的函数，且随水池的热力学性质变化，其表达式为：

qCHF(θ)=(0.034+0.003 7θ0.656)×

(6)

其中：θ为下封头壁面角度，不同壁面段需要对角度作平均；ρl和ρv分别为饱和温度下水和蒸汽的密度，kg/m3；hlv为水的汽化潜热，J/kg；g为重力加速度，m/s2；σ为水表面张力，N/m。

为保证下封头传热不恶化，必须保证下封头外壁面处于沸腾传热，即热流密度低于CHF。Theofanous基于AP600反应堆压力容器的ULPU实验测试了下封头CHF特性，根据设备Ⅱ、Ⅲ实验结果拟合得到的CHF关系式[12]如下：

qCHF-AP600(θ)=(490+30.2θ-0.888θ2+

1.35×10-2θ3-6.65×10-5θ4)×103

(7)

针对AP1000反应堆下封头，西屋公司将设备Ⅴ实验的保温层流道改造为环形流道，增强了压力容器外表面传热。虽然该实验数据尚未发布，但西屋公司认为CHF关系式可在设备Ⅲ实验数据的基础上再扩大1.44倍[13]。本文选择流道冷却剂温度最低点的热力学参数，计算不同CHF关系式随角度的分布并与实验关系对比结果，如图6所示。由图6可见，软件中原始CHF关系式计算的临界热流密度qCHF远小于两次实验值，在模拟计算过程中无法满足下封头熔融物滞留的需求。为防止计算过程中下封头传热恶化，根据已知AP1000实验关系式，扩大软件关系式常数项，拟合得到了新的传热关系式，如式(8)所示。

图6 不同关系式计算的临界热流密度随下封头角度的变化Fig.6 Change of heat flux critical with lower head angle

qCHF-fix(θ)=(0.089 5+0.001 627θ0.970 7)×

(8)

2.2 核素释放和形态模型

在反应堆运行过程中，重核裂变过程中产生的裂变产物有200多种不同核素[1]。根据裂变产物相近的化学性质，软件将核素划分为16组[19]。反应堆安全主要关心具有高裂变产额、中等半衰期和放射性特征的裂变产物，其中主要有惰性气体、挥发性裂变产物和非挥发性裂变产物3大类，其核素组成列于表1。CsI类核素是由Cs类与I类核素反应后形成的新类别化合物。计算时假设当Cs类核素与I类核素从堆芯释放出来后，立即反应产生CsI核素，并作为新的类别处理。

表1 裂变产物核素分组Table 1 Radionuclide class composition

放射性核素存在于堆芯燃料和间隙中，其质量则根据堆芯功率分布分配。当燃料棒包壳完好时，从燃料中释放的裂变产物核素积存在于间隙中；当燃料棒包壳失效时，间隙中核素便泄漏到压力容器及安全壳环境中。本文选择CORSOR-M模型[19]计算核素释放速率：

(9)

其中：k0为系数，min-1；Q为指数系数，kcal/mol；R为一定值，等于1.987 cal/(mol·K)。Xe和Ba类等核素的k0和Q系数取值列于表2。

表2 CORSOR-M释放模型系数Table 2 Coefficient of CORSOR-M release model

堆芯裂变产物通过破损燃料和堆芯-混凝土相互作用释放到控制体中，核素以气溶胶或蒸汽的形态存在。软件计算中，核素状态由蒸汽压(p)和体积温度(T)的关系式决定[19]。对于惰性裂变产物Xe，有：

lgp=1/T

(10)

该蒸汽压低于1 mmHg，即惰性核素Xe始终为蒸汽。对于非挥发性裂变产物Ba，有：

lgp=-7 836/T+6.44T≥1 000 K

(11)

当蒸汽温度超过1 000 K时，蒸汽压开始增加，非挥发性核素Ba蒸汽的释放质量增加。对于挥发性裂变产物CsI，其蒸汽压变化分为3段，如式(12)～(14)所示。

lgp=-10 420/T+19.7-3.02lgT

600 K≤T<894 K

(12)

lgp=-9 678/T+20.3-3.52lgT

894 K≤T<1 553 K

(13)

lgp=-7 304/T+7.58T≥1 553 K

(14)

CsI的蒸汽压随着温度的升高而增加。当温度低于最低温度时，蒸汽压为0，即裂变产物为气溶胶。

(15)

其中：Mw为水蒸气的摩尔质量；Ti为沉积表面的温度。对于Xe、Ba、CsI，蒸汽的饱和浓度是表面温度的应变量。而裂变产物蒸汽冷凝到任意表面i的质量Mi定义为：

(16)

放射性核素以气溶胶形态在安全壳中扩散沉积时，气溶胶颗粒的沉积速率受到重力、布朗扩散、热泳和扩散泳的作用[20]。颗粒粒径越大，受到的重力作用越大；粒径越小，受到的布朗沉积影响越大。在安全壳大空间中，重力沉积起主要作用。

3 计算结果

本文对CIS投入下先进压水堆严重事故后果进行了计算。计算选取冷段双端断裂严重事故，得到了IVR有效下安全壳压力环境和熔融池热力学结果，并分析了CIS措施对源项分布的影响。事故计算前，计算得到的反应堆稳态运行结果列于表3。可见，计算结果与反应堆设计参数相对误差较小，说明了反应堆几何模型搭建和相关参数设置的合理性。

表3 反应堆稳态参数计算结果Table 3 Result of reactor steady state parameter

3.1 事故进程

根据事故假设条件，在0时刻先进压水堆一回路冷管段发生双端断裂大破口事故(LLOCA)，能动安注系统不可使用，CIS可正常投入使用，非能动安全壳热量导出系统(PCS)、消氢系统均可投入使用，其他安全措施假设失效。计算得到的CIS投入和未投入两种假设下的反应堆事故进程列于表4。CIS通过冷却压力容器外表面带走部分堆芯热量，使得支撑结构失效和压力容器裸露发生时间变晚，整体事故进程延后。

表4 大破口事故进程Table 4 Main event of LLOCA severe accident

大破口事故发生后，一回路冷却剂通过破口迅速向安全壳隔间喷放，液态水进入到低压安全壳内迅速闪蒸，安全壳压力突升，最大压力在35 s时达到0.36 MPa。之后，一回路压力和安全壳的压力达到平衡，破口流量变小，通过热构件和水池吸收热量，安全壳内整体压力回落。由于CIS的投入，压力容器外表面和保温层间冷却通道中的冷却水带走部分堆芯热能，释放到安全壳内的热量较无CIS投入时的少，压力回落更快更低，如图7所示。随着事故的进行，冷却通道内冷却水在下封头衰变热的加热下压力和温度持续升高。约25 000 s后，安全壳内压力和温度回升。当冷却水压力等于饱和压力(图8)时，冷却水在冷却通道内从单相流动变为两相流动，以水蒸气的形态流出通道进入安全壳隔间，造成安全壳内水蒸气分压迅速增加，整体压力上升。PCS在81 482 s投入使用，带走安全壳内多余热量，阻止压力和温度继续快速增长。事故后72 h，安全壳压力维持在0.35 MPa左右，远高于无CIS投入下大破口事故后期压力。

图7 安全壳压力变化Fig.7 Change of containment pressure

图8 CIS冷却通道压力变化Fig.8 Pressure in CIS coolant channel

CIS未投入时，熔融物掉落到堆腔，与混凝土反应的大部分氢气被氢气复合器消除，但后期安全壳中氧气体积浓度低于0.5%，复合器停止工作，氢气浓度直线上升，并直接影响安全壳后期压力和温度，如图9所示。由图9可知，虽然后期氧气浓度过低，安全壳整体处于惰化状态，但是局部过高的氢气浓度会增加氢气爆炸的风险。CIS的投入避免了MCCI，锆水反应产生的氢气量是MCCI产生的2%。整个事故过程中，安全壳氢气浓度低于燃烧下限(4%体积浓度)，大幅降低了氢气燃烧和燃爆的威胁。

图9 安全壳内氢气质量变化Fig.9 Change of hydrogen mass in containment

3.2 熔融池传热分析

当堆芯开始熔化并失去完整结构，堆芯材料将掉落到下腔室。13 500 s时下封头容纳所有堆芯熔融物。熔融物总质量为149 t，主要由UO2燃料、锆和不锈钢金属及金属氧化物构成，如图10所示。随着下封头内UO2的增加，下封头内衰变热逐渐增加，最大值约为23.6 MW。

图10 下封头熔融物质量Fig.10 Mass of molten core materials in lower head

在堆芯向下坍塌阶段，堆芯释放大量裂变能和衰变能，并夹杂少量金属氧化反应产生的热量。此时，堆芯主要依靠内部冷却剂流失带走堆芯热量，其次依靠外部水池传热和热结构辐射(图11)。由于堆芯产热高于冷却剂带走的热量，堆芯热量不断增加，最高热能约为2.72×1011J(31 200 s时)。热量在熔融池内不断累积，熔融池内掉落的固体碎片不断被融化，形成稳定高度的熔融池结构，熔融池高度为段8，角度为80°，如图12所示。此后外部冷却水带走的热量超过堆芯衰变热释放，堆芯的总内能缓慢下降。

图11 堆芯能量变化Fig.11 Cumulative energy in reactor core

下封头外壁面热流密度随时间和角度的分布如图13所示。由图13可见，熔融池高热流密度主要集中在段4～6(30°～56°)之间，为氧化熔融池的边界范围。在下封头底部，碎片到下封头的传热系数和热流密度偏低。随着熔融池边界碎片的熔化(图12a、b)，氧化熔融池直接与内壁面接触，内壁面温度上升，外壁面热流密度增加(图13a)。氧化池上方边界(40°～56°)出现热流密度峰值(图13b)。随着外表面的冷却，下封头内壁面温度和外壁面热流密度整体下降，边界熔融池内出现碎片(106 000 s后)，碎片比例随着时间的增加而增加，导致段5热流密度迅速下降(图13a)。由于固体碎片具有较低的温度和传热系数，其存在会降低熔融池整体向下封头壁面的传热量。但软件计算低估了金属池向下封头壁面的传热，未形成金属层的热聚集效应[13]。且氧化熔融池和金属熔融池之间并未形成水平交界面，影响熔融池纵向和径向的传热。

图12 不同时刻下封头熔融池状态Fig.12 Molten pool configuration in lower head at different time

图13 下封头外壁面热流密度随时间和角度的变化Fig.13 Heat flux of lower head outer surface with time and angle

随着熔融物热量的变化，下封头外壁面的热流密度整体呈先增加后减小的趋势。根据计算，下封头外表面最大热流密度为1.02 MW/m2，出现在氧化池上部接触壁面，远小于实验拟合临界热流密度。拟合临界热流密度较计算得到的最大临界热流密度有40%的裕度。下封头各处的热流密度均小于临界热流密度(图13b)，保证了大破口事故下压力容器下封头的完整性；计算的熔融池内最大衰变热23.6 MW和最大热流密度1.02 MW/m2与中国核电工程有限公司自主开发程序的计算结果(25.3 MW和1.35 MW/m2)[11]接近，证明本模拟评估的下封头IVR结果有效、合理。

3.3 源项分布

由于IVR策略有效，裂变产物全部在压力容器内释放，无堆腔释放过程。裂变产物随着燃料棒包壳的失效，从燃料和燃料棒间隙释放到压力容器内。本文源项分析选取惰性气体、非挥发性和挥发性3类裂变产物。

惰性气体Xe的最终释放质量份额为100%，90%的惰性气体随着堆芯燃料的失效坍塌迅速释放，后期从熔融池中缓慢释放，到43 700 s释放停止(图14)。惰性气体以蒸汽的形态从堆芯大量释放，并通过破口从主系统迅速迁移到安全壳中。由于惰性气体不具备沉积特性，始终存在于大气中(图15)。在事故期间，累计有0.81%的惰性气体泄漏到环境中。而无CIS下惰性气体的环境泄漏质量份额为0.83%，CIS的投入对惰性气体迁移的影响小于3%。

图14 Xe在反应堆不同位置的释放份额Fig.14 Release of Xe in reactor

图15 Xe在安全壳内的形态Fig.15 Form of Xe in containment

随着事故的进行不断有非挥发性裂变产物Ba从堆芯燃料中释放，其释放份额示于图16。由图16a可见，到事故后72 h，有56.2%的非挥发性裂变产物从燃料中释放。虽然熔融池中Ba核素的释放速率低于MCCI，但是其释放时间较长，导致IVR下更多的Ba核素释放。在压力容器内释放的Ba核素有53.8%滞留在主系统内，仅向安全壳迁移2.3%。由图16b可见，无CIS投入下，Ba核素主要在堆腔内释放，大量Ba核素被直接释放到安全壳内。与无CIS投入的事故后果相比，IVR下Ba核素滞留在主系统的质量增加200倍，安全壳内质量降低91%。IVR下Ba核素泄漏到环境的比例仅0.001%，是无CIS工况下泄漏份额的26%。放射性核素Ba向安全壳以及环境的泄漏质量降低。

非挥发裂变产物主要以蒸汽的形式从熔融池释放到主系统中，并大量沉积在主系统的壁面上(图17a)。Ba核素以气溶胶的形态随着冷却剂流失从破口迁移到安全壳中，后期安全壳内Ba核素的质量不再明显增加(图16a和图18a)。对于无CIS事故，Ba核素从堆腔释放后，部分以蒸汽的形式沉积在堆腔壁面，扩散到安全壳其他隔间的Ba核素以气溶胶的形式存在，其气溶胶峰值质量为25.5 kg，是IVR下气溶胶峰值的32倍。

图16 Ba在反应堆不同位置的释放份额Fig.16 Release of Ba in reactor

图17 Ba在主系统中的存在形态Fig.17 Form of Ba in main system

图18 Ba在安全壳中的存在形态Fig.18 Form of Ba in containment

图19 无CIS下Ba蒸汽在安全壳不同隔间壁面的沉积质量Fig.19 Settlement mass of Ba vapor at containment wall without CIS

在重力作用下，悬浮气溶胶沉积到安全壳壁面和水池中。事故后72 h，无CIS下安全壳内Ba气溶胶有44%分布在水池中，56%沉积在壁面上，安全壳内悬浮气溶胶质量为0.03 kg。而IVR下，由于安全壳(事故后14 h)处于高温高湿环境，壁面上沉积的核素随着壁面冷凝水的冲刷，转移到安全壳底部水池中。最终安全壳中Ba气溶胶有96%分布在安全壳水池中，悬浮气溶胶质量降低到0.01 kg以下。

CIS下，挥发性裂变产物CsI共有95%释放。与惰性气体相同，挥发性裂变产物在堆芯坍塌过程中大量释放，如图20a所示。从堆芯释放的挥发性裂变产物有32.56%残留在主系统内，而62.39%随着喷放扩散到安全壳内。到事故后72 h，累计有0.004%的裂变产物泄漏到环境中。与无CIS事故下CsI核素分布结果(图20b)相比，残留在主系统回路的CsI核素比例增加1倍，迁移到安全壳中的质量减少21%，而泄漏到环境的CsI核素质量减少38%，说明IVR策略下CsI核素的迁移和泄漏得到有效控制。

图20 CsI在反应堆不同位置的释放份额Fig.20 Release of CsI in reactor

挥发性裂变产物以蒸汽和气溶胶的形式从堆芯释放并在主系统中扩散，如图21所示，蒸汽形态占比较大。由图21a可见，由于大破口事故，主系统内冷却水从破口快速流失，核素大部分沉积在壁面上。由于CIS的投入，压力容器壁面被冷却，沉积在主系统壁面上的CsI核素质量不变，沉积蒸汽并未像无CIS投入事故(图21b)中一样被再蒸发。在IVR事故下，放射性CsI能有效滞留在主系统中，而不发生核素向安全壳的二次转移。

图21 CsI在主系统中的存在形态Fig.21 Form of CsI in main system

CsI核素的第1次转移发生在事故后核素快速释放阶段，释放核素通过破口大量迁移到安全壳中，如图22所示。由图22a可见，CsI以气溶胶的形态在安全壳内扩散，悬浮气溶胶最高质量可达3.73 kg。随着事故的进行，安全壳内悬浮气溶胶在重力等作用下沉积在壁面和水池中。由于CIS的投入，安全壳后期大气相对湿度为1，大量水蒸气冷凝到壁面上，冲刷壁面上沉积的核素。最终CsI核素随着壁面水膜流入到安全壳底部水池。

图22 CsI在安全壳中的存在形态Fig.22 Form of CsI in containment

通过与无CIS事故下Ba和CsI气溶胶分布结果(图18b和图22b)对比，发现IVR下迁移到安全壳的Ba和CsI气溶胶总质量减少，放射性核素能被较好地控制在压力容器内，泄漏到环境的份额减小。安全壳内放射性气溶胶主要沉积在安全壳底部水池中，安全壳壁面上沉积核素的质量大幅减小，方便事故后放射性处理。总体来说，IVR下，安全壳内放射性核素能得到更好的控制和处理。