严重事故熔融池实验热工水力特性仿真分析

朱光昱，张宁娜，王昆鹏，石兴伟，左嘉旭,*，刘福东

（1.生态环境部核与辐射安全中心，北京 100082；2.国家环境保护核与辐射安全审评模拟分析与验证重点实验室，北京 102488；3.中国核电工程有限公司核电安全研究中心，北京 100840）

压水堆核电厂发生严重事故后，堆芯由于失去冷却而熔化，由此形成的高温熔融物会聚集在压力容器下封头内形成熔融池[1]。对此，部分电厂设计了堆腔注水冷却系统，通过在压力容器外部注水的方式来冷却下封头，以避免压力容器在严重事故后失效。该系统的有效性取决于熔融池施加到下封头上的热流密度是否低于冷却水侧的临界热流密度。为了给系统设计提供必要的依据，国内外设计了多个实验台架并使用多种熔融物模拟物对下封头内熔融池换热进行了研究，其中包括切片形式的COPO[2]（COrium Pool avec formation de glace aux parois）、BALI[3]（BAin Liquide）、SIMECO[4]（Simulation of In-vessel MElt COolability）和国内的 COPRA 等实验台架[5]（COrium Pool Research Apparatus）以及下封头三维缩比形式的 LIVE 实验台架（Late In-Vessel Phase Experiments）[6]。

为了准确模拟下封头熔融池内高瑞利数自然对流工况以及熔融物相变过程，LIVE 实验台架和 COPRA 实验台架采用了与原型熔融物UO2-ZrO2的相图存在相似性的20 mol%NaNO3-80 mol%KNO3混合物作为实验工质[5,6]。由于熔融池内的物理变化过程十分复杂，LIVE 实验仅测量了熔融池中的部分区域的温度值，出于弥补实验中可采集数据局限性的目的，Zhang等[7]采用基于源项的SIMPLE 算法建立仿真计算模型，通过仿真手段研究了LIVE 实验的熔融池流场和温度分层情况。上述方法需要通过自编译程序进行计算，为了获得简便的熔融池数值模拟方法，本文基于成熟的商用软件COMSOL Multiphysics 建立了一种熔融池计算模型，对LIVE 实验熔融池内的温度场和速度场进行了仿真模拟，测试了COMSOL 在熔融池仿真方面的适用性。

1 仿真模型

1.1 相变模型

20 mol%NaNO3-80 mol%KNO3的二元混合物的物性可根据Gaus-Liu[6]的研究结果设置，固相线温度为497.15 K，液相线温度为557.15 K，由固相线到液相线的相变潜热为161 960 J·kg-1，其他主要物性参数如表1 所示。结合以往仿真经验等[5,7]研究结果，由于 20 mol%NaNO3-80 mol% KNO3的二元混合物具有较大固液相线温差，采用相变模型可以更好地对流场进行模拟。COMSOL 中相变模型原理如图1 所示，图中θ为相体积分数，ΔT为固液相线温差，固相线温度和液相线温度的均值Tpc即为COMSOL 中相变温度，对于在固液相线之间的混合物，在计算过程中采用固液两相按体积平均的物性参数进行计算。

图1 相变模型原理图Fig.1 The schematic of phase change

表1 20 mol% NaNO3-80 mol% KNO3 混合物物性Table 1 The property of 20 mol% NaNO3-80 mol% KNO3

根据Zhang等[5]研究，可采用低雷诺数k-ε模型对熔融池湍流场进行计算，混合物的粘度设置为固液两相的平均值，本文通过设置较大的粘度来模拟固相不存在流动性的情况，当固相的粘度系数高于103Pa·s 后对计算的影响很小。

1.2 计算域设置

本次模拟工作针对LIVE-L5 实验[8]进行。图2 为LIVE 实验台架示意图，其整体为典型压水堆压力容器下封头的1/5，内径为1 m。实验过程中通过熔盐中的电加热器模拟熔融物衰变热，下封头外部设有与 IVR（In-Vessel Retention，IVR）技术类似的冷却流道。实验中加注的熔盐为 210 L，形成的熔融池高度为0.43 m。熔融池上部为熔盐液面、水冷壁面和顶部盖板组成的封闭腔体。腔体内辐射换热网格图如图3 所示，其中顶部壁面保温性能很好可认为是重辐射面，则熔融池上表面向水冷壁面辐射热流Φ可通过式（1）计算：

图2 LIVE 实验系统示意图Fig.2 The schematic of the LIVE test vessel

图3 辐射传热等效网络图Fig.3 The radiation heat transfer network

图3 和式（1）中：Eb1、Eb2和Eb3以及J1、J2和J3分别为熔盐液面、水冷壁面和顶部盖板的黑体辐射热流密度和有效辐射热流密度；ε1和ε2为熔盐液面和台架内不锈钢壁面的发射率，参考以往仿真经验分别取0.44 和0.8[5]；A1、A2和A3为分别为熔盐液面、水冷壁面和顶部盖板的面积；X1,2、X1,3和X2,3为各表面间的角系数；Rt为图3 中所示的各辐射热阻的总等效热阻，可以根据各热阻串并联关系计算，具体计算方法在文献［9］中有详细介绍。

式中：σ——黑体辐射常数；

T1——熔盐液面的温度；

T2——水冷壁面的温度。

由于熔盐的液位很高，此时水冷壁面在高度方向上的弯曲程度很小，因此可近似认为封闭腔体是扁平的圆柱体，并据此计算各辐射面之间的角系数。通过将黑体辐射公式带入式（1）可导出式（2），此时将式（2）的总等效热阻的倒数作为一个等效的表面发射率，水冷壁面的温度作为外界环境温度，便可将辐射换热区转化为熔融池上表面对外界环境的辐射换热过程，这样熔盐液面可设置为COMSOL 软件中的表面对环境辐射换热边界条件，从而不需要对辐射换热的区域进行建模。此外，模型还考虑了顶部通入氮气带走的热量，根据Zhang等[10]估算这部分提供的冷却效率最大不超过11.7 W。

综上所述，最终建立的轴对称计算模型如图4 所示。由于LIVE 实验相关文献中并未给出冷却流道的具体尺寸，因此只能估测水冷壁面对流冷却换热系数。对此，本文根据实验中的各稳态工况中[8]下封头壁面导出热量占总加热功率的比例来校核对流换热边界的热通量。

图4 计算域和网格划分Fig.4 The calculation domain and mesh generation

1.3 边界层和网格敏感性

参考以往的研究结果[7,10]，仿真得到熔融池内流速均很低，因此边界层变化对速度场中最大速度和平均速度影响很小。然而，通过预计算发现，边界层的设置对水冷却壁面附近相变计算的影响较为明显，经过测试，采用COMSOL 默认的边界层设置，在水冷壁面内设置8 层边界层网格即可满足一般的计算需要。图5 所示为加热功率为18 kW，初始化温度为615 K 时，不同网格数下熔融池平均温度随时间的变化。随着网格加密，计算得到的平均温度逐渐降低，在网格数加密至4 121 后，进一步加密的计算结果变化率在1%以内，因此最终计算采用的网格数为4 121。

图5 不同网格数下的熔融池平均温度Fig.5 The average corium pool in different mesh number

2 计算结果分析

图6 和图7 所示为热功率为18 kW 时熔融池的计算温度分布和速度场，在稳态工况下，冷却壁面附近熔盐受冷后形成了沿壁面向下的自然对流，其他不同高度区域存在多个涡流。受涡流搅拌的熔融池中、上部温度分布较为均匀，而底部则的呈现出了明显的热分层结构。图8 所示为热功率为5 kW 和18 kW 时稳态情况下，距离中心轴365 mm 处的熔融池温度随高度的变化。其中底部区域的计算值与实验值[7]相比偏差较大，这可能是因为LIVE 实验的下封头中底部布置了大量的加热和测温设备，实际实验条件下该区域的涡流搅拌作用可能并不明显，因此温度场仿真结果并不准确。对于熔融池中部和顶部区域，计算得到的温度分布与实验值十分接近，误差不超过4%。

图6 热功率18 kW 时熔融池的热分层结构Fig.6 Thermal stratification of the corium pool under the heating power of 18 kW

图7 热功率18 kW 时熔融池内速度场Fig.7 The velocity field of the corium pool under the heating power of 18 kW

图8 热功率为5 kW 和18 kW 时熔融池内温度稳态分布Fig.8 The steady corium pool temperature profile under the heating power of 5 kW and 18 kW

图9 和图10 示出了热功率为18 kW 时熔融池内的结壳厚度分布以及结壳厚度的计算值和实验值[7]。受重力和沿冷却壁面向下的自然对流影响，凝固的熔盐会向下封头底部迁移，因此结壳厚度总体呈现随着熔融池圆弧壁面径向角度增加而减小的趋势。在严重事故缓解过程中，堆内熔融物结壳可以形成一层热阻，从而有效阻止下封头的不锈钢壁面发生熔化。根据当前的实验和仿真得到得结壳分布情况，即使仅考虑单层熔融池，下封头的高相位角区域也属于薄弱环节。当前熔融物堆内滞留技术（Invessel Retention，IVR）普遍采用了从下封头底部向堆腔注水的方式，由于冷却流道的截面积随着相位角增加而增大导致在高相位角区域冷却水流速降低，同时冷却水在底部受热导致高相位角区域水温相对较高，因此高相位角区域的冷却效果相对底部较差，更容易发生熔穿等后果，在后续IVR 技术工程设计过程中，应考虑在下封头外冷却流道的高相位角区域增加强化换热设计。

图9 18 kW 时固相体积分数Fig.9 The soild volume fraction under the heating power of 18 kW

图10 18 kW 时冷却壁面上结壳厚度Fig.10 The crust thickness along the vessel wall under the heating power of 18 kW

3 结论

基于COMSOL Multiphysics 软件建立了核电厂严重事故后下封头熔融池的计算模型，采用LIVE 实验结果对模型进行了验证。仿真结果表明熔融池中除了沿冷却壁面向下的自然对流外，熔融池的上、中部还同时存在大量涡流，受其搅拌的区域分布较为均匀，而熔融池底部黏度较大无法产生涡流因此呈现出了明显的热分层结构。在稳态下，仿真得到的中部和顶部熔融池温度分布以及冷却壁面附近的结壳厚度与实验结果十分接近，说明当前COMSOL 模型适用于熔融池的仿真计算。