IVR策略下反应堆压力容器变形数值模拟

朱光昱， 靖剑平， 石兴伟， 左嘉旭， 刘宇生，4， 温爽*

(1.生态环境部核与辐射安全中心， 北京 100082； 2.国家环境保护核与辐射安全审评模拟分析与验证重点实验室， 北京 102488; 3.中国核电工程有限公司， 北京 100084; 4.哈尔滨工程大学核科学与技术学院， 哈尔滨 150001)

核电厂严重事故中，堆芯构件由于失去冷却而熔化形成高温熔融池[1]。在中国的第三代先进压水堆发展过程中，倾向于采用堆内熔融物滞留技术(in-vessel retention, IVR)作为核电厂严重事故的缓解措施[2]，该策略通过对反应堆压力容器(reactor pressure vessel, RPV)外壁面进行冷却带走堆内熔融池产生的衰变热，从而保持 RPV 的结构完整性。在高温条件下，RPV下封头会发生明显的热膨胀，从而导致其外部冷却流道结构与初始设计存在偏差[3]。如果上述变化导致局部传热恶化，熔融池导出热量大于局部CHF值，则可能发生RPV熔穿进而导致IVR策略失效，因此在IVR策略设计过程中，需要对严重事故后RPV的形变进行研究。

在高温熔融池的作用下，RPV壁面会发生熔化进而降低RPV的壁面厚度，最终加剧了RPV 壁面的变形。因此在对RPV下封头的热膨胀变形影响进行评估时，也必须考虑钢结构熔化带来的影响。在熔融池的演化过程中，由于其中的轻金属与氧化物不相溶，最终会形成双层熔融池结构[4]。顶部较轻的金属层导热性能良好，产生的热聚集效应会使与金属层接触的RPV下封头内壁面的熔化可能更加明显，在研究过程中同时需要考虑熔融池分层带来的影响。

在工程设计上，一般通过数值模拟的方式对严重事故后RPV下封头形变进行分析。温爽等[5]和田欣鹭等[6]基于ANSYS研究了热膨胀和高温蠕变对下封头失效的影响。罗娟等[7]通过弹塑性模型分析了下封头自重、熔融池内压和外部冷却水压力对下封头变形的影响。刘兆阳等[8]采用ABAQUS软件对示范快堆熔融物收集装置进行了安全分析。高永建等[9]在完成温度场分析的基础上，通过Larson-Miller参数积累损伤理论对RPV蠕变进行了研究。上述工作中均采用了将熔融池内部与下封头解耦的方式完成力学分析，可能引入一定的不确定性。对此，现采用多物理场耦合仿真平台COMSOL，建立一个热工水力和固体力学耦合计算模型，对单层和双层熔融池结构下RPV下封头的形变进行研究，为严重事故下RPV外壁面的换热情况研究提供必要的输入参数。

1 计算模型

1.1 计算域设置

图1所示为参考HPR1000堆型设计参数建立的二维轴对称双层熔融池计算模型，熔融物区域设置为流体域，其中氧化物层高度为1.58 m，金属层高度为0.56 m。下封头区域设置为固体域，下封头内径为4.4 m，下封头壁面厚度为0.15 m。图1中r表示轴对称模型的半径，g为重力加速度，模型中取9.8 m/s2。单层熔融池计算模型与图1类似，但不存在分层结构，总高度为1.8 m。为了减少边缘效应对力学计算的影响，下封头顶部区域延伸了长度为3 m的不锈钢壁面，图1中并未完全示出。

当前COMSOL计算模型中，熔融池和下封头的温度场采用固体和流体传热模块计算，熔融池湍流场采用低雷诺数k-ε模型计算，下封头的热膨胀过程采用固体力学模块计算。熔融物和下封头的相变过程采用相变模型计算。计算域中，金属层顶部和下封头延伸部分的内壁面均设置为表面对环境辐射条件，发射率分别设置为0.45和0.8[4]。下封头外侧设置为400 K的恒温壁面[5]。对于双层熔融池结构，可以认为衰变热仅存在于氧化物层中，计算公式[4]为

Q=0.125fdP0t-0.28

(1)

式(1)中：Q为堆芯衰变热功率，MW；fd为衰变热修正系数，设置为0.75[12]；P0可认为是反应堆正常运行时的热功率，本文取3 050 MW；t为停堆时间。

图1 计算域和网格划分示意Fig.1 Calculation domain and mesh generation

1.2 熔融池物性设置

熔融池内通过引入相变模型模拟凝固过程，材料密度采用Boussinesq近似处理。参考VULCANO熔融物铺展实验相关数值模拟经验[10]，熔融池中糊状区的黏度为

μ=μle2.5Cθs

(2)

式(2)中：μl为液相黏度；系数C一般取3.5～8[11]；θs为固相体积分数，可直接在计算过程中实时调用参数。

单层和双层熔融池内其他计算所需物性参照文献[12-14]设置，膨胀系数参考液态金属和陶瓷氧化物设置，具体如表1和表2所示。

表1 单层熔融池物性参数表

表2 双层熔融池物性参数表

1.3 压力容器物性设置

计算过程中依然采用相变材料模型计算下封头熔化过程，并为此设置了10 K的相变区间，RPV的具体物性设置如表3所示。

为了引入固体力学模型，下封头需要设置为固体计算域，这导致无法直接计算RPV壁面熔化后混入熔融物的过程，因此对下封头物性进行了如下设置。通过预计算判断，当前熔融池参数设置下，仅下封头与金属层接触区域会发生熔化，而熔化后的钢液相对于金属层总质量而言很小，因此可以将下封头熔化后区域的材料热工物性参数设置为与金属层一致。

此外，对于下封头熔化区域，在材料本身的导热系数之外将湍流导热添加至熔化后区域的导热系数当中，湍流导热的计算公式为

(3)

式(3)中：kT为湍流流动产生的导热系数；PrT为湍流普朗特数，取0.9；μ为湍流动力黏度；Cp为比热容，两项均参照熔化后的金属层参数设置。

PRV的杨氏模量通过插值函数设置[7]，具体如表4所示。对于固体力学计算，将熔化区域的力学参数直接设置为0会导致计算发散。对此，对表3中熔化后的泊松比和表4中液相线以上的杨氏模量设置了极小的数值，以达到忽略熔化区域力学影响的目的。

表3 压力容器物性参数表

表4 不同温度下的杨氏模量

1.4 网格敏感性分析

采用瞬态方式求解，通过设置较长的计算时长获得各物理场达到稳定状态时的参数。考虑停堆3 h后已经形成双层熔融池结构，此时氧化物层的衰变热为21.23 MW。图2所示为不同网格数计算得到的RPV下封头外壁面应力分布。当网格加密至30 000以上后，网格变化对计算得到的应力分布变化的影响降低至 1 %以内，双层熔融池后续计算采用的网格数为32 298。对于单层熔融池计算模型，采用相同的方法进行网格敏感性分析，最终采用的网格数为21 279。

图2 不同网格数下封头外壁面应力Fig.2 Stress distribution on the outer surface of lower head under under different mesh number

2 数值模拟结果分析

2.1 熔融池温度场分析

图3 单层熔融池和双层熔融池的温度分布 Fig.3 Temperature distribution of ingle-layer corium pool and two-layer corium pool

计算过程中，熔融池以参考温度作为初始化温度，下封头的初始温度设置为450 K。图3为衰变热为23.78 MW时的单层熔融池和衰变热为21.23 MW时的双层熔融池温度分布情况。本文中，相位角指的是与计算域中对称轴的夹角。对于单层熔融池结构，下封头中相位角大于25°的区域均不同程度的超过了不锈钢的熔点。双层熔融池中，与氧化物层接触的下封头区域温度较低不会超过熔点，对于与金属层区域接触的区域，在热聚集效应作用下，这部分下封头内壁温度会显著升高并超过不锈钢的熔点。

图4 单层熔融池和双层熔融池的下封头von Mises应力分布 Fig.4 Distribution of von Mises stress in lower head of single-layer corium pool and two-layer corium pool

2.2 RPV热膨胀形变分析

图4为单层熔融池和双层熔融池状态下封头的von Mises应力分布情况，高温作用下不锈钢壁面的熔化情况也可在图中看出。可见下封头热膨胀产生的应力主要集中在压力容器外壁面处。对于发生熔化的区域，在较高的温度和壁面厚度由于熔化削弱的双重作用下，下封头外侧von Mises应力可达其他区域的两倍以上。

图5为单层熔融池和双层熔融池状态下RPV的热膨胀形变量分布，为了便于观察，图5所示的膨胀量显示为实际值的5倍。计算域中，RPV延伸区域顶部3 m处的边界条件为纵向固定约束，导致计算域内所有的纵向热膨胀形变会逐渐向下累积。同时，由于对称轴相当于横向位移约束，所以横向形变的则会逐渐向下封头的高相位角积累。在各处的局部热膨胀形变和上述两种积累效应的共同作用下，单层熔融池的最大形变出现在相位角40°左右的区域，最大形变量达到19.7 mm。双层熔融池的最大形变量出现在与金属层接触的区域，其中最大的膨胀位移达到了16.9 mm。对于下封头的底部区域，虽然两种熔融池状态下的局部的热膨胀形变均很小，但在纵向的积累效应作用下依然存在约15 mm的膨胀位移。

2.3 RPV内压影响分析

除衰变热和下封头壁面熔化外，RPV内压也是其形变的重要影响因素[5-7]，在泄压系统的作用下，可以认为熔融池形成时，RPV的内压一般不会超过1 MPa。本文选取双层熔融池结构对内压影响进行分析，通过施加边界载荷改变内压。图6所示为不同衰变热和内压条件下的下封头外壁面膨胀位移随相位角的变化。图6中衰变热21.23 MW和14.39 MW分别对应停堆后约3 h和12 h的时刻。衰变热为21.23 MW时，1 MPa和2 MPa内压分别使得最大热膨胀位移增加了约0.9 mm和1.4 mm，这将导致RPV对冷却流道的挤压更加明显，使外流道传热性能发生变化，有可能造成冷却恶化和IVR策略的失效。随着衰变热降低，下封头各相位角的径向膨胀位移呈现降低的趋势，然而当内压较高时这种收缩较为缓慢，从而导致RPV冷却流道变形时间增长。

图5 单层熔融池和双层熔融池的RPV的热膨胀形变量 Fig.5 RPV thermal expansion deformation of single-layer corium pool and two-layer corium pool

图6 下封头外侧径向热膨胀位移Fig.6 Thermal expansion deformation displacement of lower head outboard

3 结论

根据HPR1000堆型设计参数，采用COMSOL软件搭建了一个耦合计算模型，对严重事故后单层和双层熔融池结构下的RPV热膨胀形变进行了分析。得出以下结论。

(1)严重事故后RPV内的熔融池将导致IVR下封头内壁面局部区域发生熔化，受局部高温和不锈钢壁面变薄的共同影响，RPV外冷却流道受到明显挤压。对于单层熔融池结构，最大形变出现在相位角40°左右的区域，最大形变量达到19.7 mm。对于双层熔融池结构，最大形变出现在相位角80°左右的区域，最大形变量达到16.9 mm。

(2)RPV膨胀形变量随着内压增高而明显增加，由此可能导致外流道传热性能发生变化，进而造成冷却恶化和IVR策略的失效，在工程设计过程中应予以关注。