基于COMSOL的COPRA熔融池实验仿真研究

朱光昱 闵金坤 张 丽 元一单

1（中国核电工程有限公司核电安全研究中心 北京 100840）

2（生态环境部核与辐射安全中心 北京 100082）

3（清华大学工程物理系 北京 100084）

熔融物堆内滞留技术（In-vessel Retention，IVR）是目前广泛采用的核反应堆严重事故堆芯熔融物滞留方案，其有效性取决于反应堆压力容器下封头壁面是否可以有效导出熔融物的衰变热。为了研究反应堆严重事故后压力容器内熔融池与冷却壁面的换热特性，国内外设计了许多压力容器下封头二维切片 实 验 台 架 ，例 如 COPO［1］（COrium Pool avec formation de glace aux parois）、BALI［2］（BAin Liquide）、SIMECO［3］（Simulation of In-vessel MElt COolability）以及由中国核电工程有限公司和西安交通大学针对华龙一号堆型搭建的COPRA 实验台架［4−5］（COrium Pool Research Apparatus）。 其 中 ，COPRA实验采用了穆尔比为20%NaNO3-80%KNO3的二元混合物作为实验工质，该二元混合物与反应堆内真实熔融物UO2-ZrO2的相图存在相似性［4］，可以较为准确地模拟下封头熔融池内的熔融物相变行为和高瑞利数自然对流工况，为IVR 相关系统设计提供了必要的技术储备。由于熔融池内的物理变化过程十分复杂，限制了COPRA 实验中测量手段，需要建立仿真计算模型来弥补实验的不足。对此，本文基于COMSOL Multiphysics多物理场计算软件建立了一种熔融池计算模型，通过耦合湍流场和传热相变过程，实现了对COPRA 实验熔融池内的流动传热、相变过程的仿真模拟。

1 仿真模型

1.1 数学物理模型

采用COMSOL 软件中的非等温传热模块耦合实验中的流体传热、相变和湍流流动。穆尔比为20%NaNO3-80%KNO3的二元混合物具有较大固液相线温差，适宜采用流体传热模块中的相变材料模型进行模拟，其原理如图1所示，温度低于固相线区域采用固相物性计算，温度高于液相线区域采用液相物性计算，在固液相线之间采用固液两相按体积平均的物性参数进行计算。图1中θ为相体积分数；ΔT为固液相线温差；Tpc为COMSOL中的相变温度，这里设置为固相线温度和液相线温度的均值。参考Gaus-Liu［6］和 Janz［7］的 研 究 结 果 ，20%NaNO3-80%KNO3的二元混合物的固相线温度为497.15 K，液相线温度为557.15 K，由固相线到液相线的相变潜热为 161 960 J·kg−1，其他主要物性参数如表 1所示。

图1 相变模型原理图Fig.1 Schematic diagram of phase change model

湍流流动模块采用低雷诺数k-ε模型进行计算，参考相变材料模型设置思路，该模块中的粘度项设置为固液两项的体积平均值，可通过设置较大的粘度来模拟固相不存在流动性的情况。通过预计算对比，当固相的粘度系数高于103Pa·s后对计算的影响很小。

求解过程采用分离式求解器，时间步进选择向后差分，计算容差因子为0.5，计算加速性和稳定性设置为Anderson加速度。

1.2 计算域设置

图2 为参考COPRA 实验台架的切片结构建立的二维计算模型。台架整体呈半径为2.2 m的1/4圆形，加注熔盐后台架可分为底部的熔融池区域和顶部的辐射换热区域，熔盐加注的最大高度为1.9 m。其中1/4圆形半径上的边界（c、d和f）为绝热边界，圆弧边界（a和b）外侧为水冷流道，其平均温度约维持在323.15 K，可以认为是等温边界条件。

顶部的辐射换热区域可认为是4个表面组成的封闭腔辐射传热问题，其辐射换热网格图如图3 所示，其中绝热壁面c、d可认为是重辐射面，则熔融池上表面（e）向水冷壁面辐射热流Φeb可通过式（1）计算：

表1 20%NaNO3-80%KNO3混合物物性Table 1 Property of 20%NaNO3-80%KNO3

图2 计算域和网格划分示意Fig.2 Schematic diagram of calculation domain and mesh generation

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

图3和式（1）中：Ebb、Ebc、Ebd、Ebe分别为各表面上的黑体辐射热流密度；Jb、Jc、Jd、Je分别为各表面上的有效辐射热流密度；ε1和ε2为熔融池上表面和台架内不锈钢壁面的发射率，参考张卢腾［4］的成果分别取0.44 和 0.8；Ab、Ac、Ad、Ae分别为各表面的面积，对于二维模型，可通过将切片台架的厚度进行归一化计算各面的面积；Xe，b、Xe，d、Xe，c、Xb，d、Xd，c、Xb，c为表面 e和b、e和d、e和c、b和d、d和c、b和c间的角系数，可将封闭腔简化为长方形结构通过查表［8］获得；Rt为图3中所示的各辐射热阻，总等效热阻需根据各热阻串并联关系计算，具体计算方法在文献［8］中有详细介绍。

通过引入黑体辐射公式可将式（1）化为式（2）的形式，此时可将总等效热阻的倒数作为一个等效表面发射率考虑，则辐射换热区域可简化为熔融池上表面对外界环境的辐射换热过程，因此仅对由a、e和f 三个边界围成的熔融池区域的建立计算域即可。其中，将a 设置为边界恒温条件，e 设置为COMSOL软件中的表面对环境辐射换热边界条件，f设置为绝热边界条件。

式中：σ为黑体辐射常数；Te为熔融池上表面的温度；Tb为水冷壁面的温度，这里设置为“表面对环境辐射换热”边界条件中的外界环境温度。

1.3 网格敏感性分析

根据网格敏感性分析，当网格数大于10 924后，网格变化对计算得到的温度场分布和熔融物结壳厚度变化率降低至1%以内，所以本文仅选取流场中的最大计算流速作为分析依据。图4为加热功率为10 kW时的网格敏感性分析结果，随着网格加密，计算得到的流场中最大流速逐渐增加，在网格数加密至15 462后网格宽度对计算得到的最大速度变化率约为2%，因此最终计算采用的网格数为15 462。计算采用的初始时间步长为1 ms，由于COMSOL软件会自动根据局部误差对计算时间步长进行调整，因此未对时间步长的敏感性进行分析。

图4 不同网格数下的最大计算速度Fig.4 Maximum simulation velocity under different mesh number

2 计算结果分析

COPRA 台架通过在熔融池中的电加热器模拟熔融物衰变热，对于实验中的T1 工况［4］，加热功率按照A 阶段（功率15 kW，持续时间 23 700 s）、B 阶段（功率10 kW，持续时间40 578 s）、C 阶段（功率15 kW，持续至稳态）变化。计算过程中，参照文献［5］中的实验台架设计信息，将不同的加热功率折合成体积热源并采用阶梯函数控制，从而实现对实验过程的模拟。图5为A阶段熔融池的计算温度稳态分布，与实验结果类似，在熔融池内的温度总体呈现随高度增加逐渐升高的趋势，且不同高度出现了明显的热分层现象，在底部结壳的区域的温度梯度较大，中部和顶部的温度梯度较小。

图5 A阶段熔融池的热分层结构Fig.5 Thermal stratification of corium pool in stage A

图6为各阶段稳态下的熔融池温度分布的实验值和计算值。其中，底部区域的计算偏差相对较大，这可能是由于当前计算模型没有引入在熔融物凝固时会在硬壳和下封头间形成气体间隙导致的［4］。对于熔融池高度0.5 m以上区域，计算得到的温度分布与实验值十分接近，最大计算误差不超过1.5%。

图6 各阶段熔融池内温度稳态分布Fig.6 Steady temperature profile in different stages of corium pool

图7（a）为B 阶段的流场，熔融池内形成了由冷却壁面向下流动，在中部受热后在浮力作用下缓慢向熔融池顶部和冷却壁面流动的自然对流趋势。受自然对流和顶部辐射换热的影响，在熔融池内存在大量的涡流，这些涡流搅混的区域温度分布趋于均匀，因此易于在熔融池内形成热分层结构。图7（b）为B阶段的固态体积分数，在冷却壁面附近，受高速向下主流和重力的共同作用，凝固形成的糊状物会更容易向熔融池底部转移并最终结壳，从而使熔融物结壳的厚度呈现沿冷却壁面由熔融池底部到顶部逐渐降低的趋势。

图7 B阶段熔融池内速度场(a)和固相体积分数(b)Fig.7 Velocity field(a)and soild volume fraction(b)of corium pool in stage B

3 结语

采用COMSOL Multiphysics软件的非等温流动模块、相变材料模型建立的熔融池计算模型对COPRA实验进行了仿真模拟，结果表明：

1）稳态熔融池内自然对流形式为：由冷却壁面向下流动，在中部受热后在浮力作用下缓慢向熔融池顶部和冷却壁面流动。

2）在自然对流和顶部辐射换热的共同影响下，熔融池内会形成了大量的涡流。

3）在冷却壁面向下主流和重力作用下，熔融物凝固硬壳呈现沿冷却壁面由熔融池底部到顶部逐渐变薄的趋势。