压力容器下封头大变形模型的开发及在FOREVER实验中的应用分析

翟润泽 张斌,2 杨皓 高鹏程 唐绍伟 单建强,2

1（西安交通大学 能源与动力工程学院 西安 710049）

2（西安交通大学 动力工程多相流国家重点实验室 西安 710049）

自20世纪80年代以来，尤其是三哩岛核电厂和福岛核电厂发生堆芯熔化严重事故后，对严重事故研究一直是核安全领域的重点课题。这两次事故提醒人们核电站存在发生严重事故的可能性，在核电站的设计以及应对措施上，需要考虑堆芯熔化以及熔融物重新定位到下封头的情况。将堆芯熔融物包容在压力容器内的策略（In-vessel Retention，IVR）是作为缓解严重事故的一项重要措施，该策略已成功应用于AP600、AP1000、华龙一号及CAP1400等先进堆型的严重事故管理中。近年来，随着对清洁能源的重视，各国政府和公众越来越关注核电厂的安全问题［1-2］。

IVR策略有效性的量化直接取决于两个参数：熔池加载于压力容器壁面的热流密度和压力容器外部冷却的排热能力［3］。在IVR策略中，熔融物作用于压力容器下封头，使其产生不可忽视的塑性变形［4］，下封头的变形会改变堆腔的冷却流道，这将直接影响压力容器外部冷却的排热能力和IVR策略的成功实施，建立能准确预测在IVR策略中下封头行为的模型具有重要意义。

国内外已经开展了许多堆芯熔化严重事故下反应堆压力容器（Reactor Pressure Vessel，RPV）热力学行为相关的实验研究和数值模拟计算。较为著名的下封头形变实验有：LHF（Lower Head Failure）实验、OLHF（OECD Lower Head Failure）实验和FOREVER实验。1997年，美国SNL（Sandia National Laboratories）实验室的Chu等［5］开展了LHF实验项目，研究了反应堆容器在严重事故条件下由于热负荷和压力载荷而造成的失效，研究重点是局部传热分布的峰值效应，以及贯穿件的影响。OLHF实验是LHF实验项目的扩展，深入研究了在低到中等冷却剂系统（Reactor Coolant System，RCS）压力（2～5 MPa）和下封头壁厚的典型温度梯度（200～400 K）下，RPV下封头失效的模式、时间和形变［6］。LHF实验和OLHF实验都是使用加热器加热容器内气体来模拟堆芯碎片向反应堆容器的能量传递，存在一定的失真。1999年，KTH（Kungliga Tekniska Högskolan）的Sehgal等［7］开展了FOREVER实验，实验容器的几何比例与原型LWR（Light Water Reactor）相比是1:10，实验是在1 473～1 573 K下浇注二元氧化物（CaO+B2O3）熔体来模拟熔池与下封头的换热过程，再现了原型的熔池对流过程和容器壁的温度场，能够模拟下封头所受重力载荷，是作为下封头大变形模型验证实验的很好选择。

在数值模拟计算方面，国际上广泛使用的一体化严重事故分析程序，例如MAAP（Multisensory Attention Assessment Protocol）［8］、MELCOR（Melanoma Clinical Outcomes Registry）［9］和ASTEC（Accident Source Term Evaluation Code）［10］，由于缺少计算RPV热力学行为的机理模型，不能分析事故过程中下封头的整体形变过程，无法准确预测失效时间、模式和位置等，模拟效果较差。例如，MELCOR程序简单假设RPV累积损伤份额超过18%时，则判定RPV发生破裂［10］。因此，开发一种准确度高的机理模型是很有必要的。

针对ISAA（Integrated Severe Accident Analysis）程序LHTCM（Lower Head Thermal Creep Module）的简化薄膜应力模型十分简单和缺乏计算变形模块的问题［11］，本文从机理出发，基于Timoshenko板壳理论、Nortron蠕变定律和大变形塑性理论开发了机理模型—下封头大变形模型，并将该模型集成到一体化严重事故分析程序ISAA中，由于OLHF实验有一定的失真，本文选择在FOREVER-EC2实验中进行应用研究，验证该模型在压水堆堆芯熔化严重事故中预测RPV整体形变以及失效时间和位置等情况的准确性。

1 下封头大变形模型

1.1 概述

本文基于大变形塑性理论构建的下封头大变形失效模型的计算流程如图1所示，首先确定当前时刻的6个待定参数，按角度对下封头进行分段，利用大变形理论计算出下封头当前形状，通过受力分析分段计算出此时的应力和应变，当下封头节点计算完毕，利用虚功原理计算下一时刻的6个待定系数，最后根据失效准则判别下封头是否失效，若失效则结束计算，若不失效将进行下一时刻计算。

图1 下封头大变形失效模型的计算流程图Fig.1 Calculation flowchart for the large-deformation model of lower head

1.2 力学分析

下封头的角度以及半径等几何关系如图2所示，角度方向将下封头划分为N个单元，每个单元再沿径向或厚度方向划分j个节点进行计算，下封头大变形模型的具体步骤如下。

图2 下封头模型图Fig.2 Schematic of the lower head

1.2.1 参考Timoshenko板壳理论确定下封头节点受力情况

在RPV单元的中心点P处建立三维坐标轴，如图3所示，半径定义为RPV单元中心点到原点的距离。坐标轴的三个方向分别代表：壳体半径的切线、壳体环向的切线和壳体表面的法线。根据Timoshenko板壳理论，下封头壁面满足如下三个力学平衡方程：

图3 下封头节点受力分析Fig.3 Stress analysis of the lower-head node

其中：式（1）对应与球面剖面线相切方向上力的平衡；式（2）对应于与法线方向上力的平衡；式（3）定义了力矩的平衡。N1、N2分别为单位长度上的环向和切向应力，N·m-1；Q为横向剪切力，N·m-1；M为弯矩；x为RPV球壳当前点距离Z轴的距离，m；r1为曲率半径，m；ψ为表面切线与水平线的倾角，rad；Pn为半径方向压力，Pa；Pt为切线方向压力，Pa；e为当前RPV的厚度，m。

对于法线方向上，O'P以下部分承受自身重量WO'以及下腔室压力p，在壁厚中心点受力分析如图4所示。

图4 下封头壁厚中心点受力分析Fig.4 Stress analysis of the center node of the wall thickness of the lower head

在壁厚中心点处沿Z轴方向建立力学平衡方程，即：

下封头O'P以下部分的自重表示为：

通过式（4）和式（5），可得到N2表达式为：

将式（6）代入式（1）得到N1的表达式为：

在横向剪切力Q确定之后所有应力均可求解。Rabotnov公式给出了横向剪切力Q的近似解析表达式：

式中：R0、e0分别为初始的平均半径和壁厚，m；φ为半球壳的初始极角，rad；ν为泊松比；λ为修正项。

1.2.2 计算下封头节点应力和应变

根据第四强度理论假设即可得出：

不锈钢在高温环境下会产生蠕变应变，蠕变是与高温相联系的，持续的高温蠕变导致容器壁发生大变形。本文采用Norton蠕变定律计算高温蠕变，表达式如下：

上述介绍的等效应变率以及等效应力都是以每个RPV单元的中心位置确定的，实际上由于蠕变对温度十分敏感，下封头内壁面的应力要低于外壁面的应力，需要计算沿壁厚方向上的应力变化。由于下封头形变过程中没有发生明显的倾斜弯曲，与拉伸应变相比壳体经历的弯曲很小，因此，可以近似认为等效应变率沿壁厚保持不变，可以得到各个沿壁厚各个节点的应力：

其中：ξi表示下封头沿壁厚方向划分的节点，通过沿壁厚等效应变不变假设可以计算出壁厚上各个节点的等效应力值。

1.3 大变形塑性理论

下封头的初始结构与其形变后形状之间的对应关系如图2所示，建立的壳体形变位移函数表达式为：

式中：φ、θ分别为变形前后的角度，rad；R0和r为下封头初始半径和形变后的半径，m。

根据虚功原理，不管允许的虚位移如何，内力做功和外力做功始终相等。因此，得到待定系数的显式递归求解如下：

式中：矩阵｛F｝和矩阵［Γ］是通过沿下封头曲线对模型提供的解析表达式进行数值积分而建立的，6个系数确定后即可确定下封头壁面的形变。矩阵［Γ］的最终表达式由下式给出：

1.4 失效准则

RPV由于高温蠕变造成的损伤遵循Kachanov提出的蠕变损伤准则，损伤份额如下式：

式中：D为损伤份额，当D达到1.0时，结构发生失效；Ak为温度相关的材料参数；k、B与材料和温度相关，具体取值详见参考文献［13］。

2 FOREVER-EC2实验的模拟

2.1 FOREVER-EC2实验简介

FOREVER实验于1999-2002年期间在KTH开展，反应堆压力容器下封头通过1∶10几何缩比模拟［14］。FOREVER-EC2实验容器由一个焊接在半球（材料为法国16MND5钢）上的筒体组成。容器内半径为188 mm，壁厚约15 mm。FOREVER-EC2实验设施由感应炉、容器、加热器、电源等组成。感应炉可制备高达15 L的熔体，可远程倾斜以将熔体添加到容器中。实验是在1 473～1 573 K下浇注二元氧化物（CaO+B2O3）熔体进行的。然后用加热器加热熔体，使其保持在1 373～1 473 K，并用氩气加压。

2.2 数值建模

下封头大变形模型是一个下封头行为模型，本文将下封头大变形模型集成到自主开发的一体化严重事故分析代码ISAA［15］中，取代现有的简单下封头行为模型。FOREVER实验的数值建模基于ISAA程序，建模数据主要来自于文献［14］，建立的数值模型如图5所示，控制体用CV-XXX表示，连接管用JP-XXX表示，换热壁用HTW-XXXXX表示。其中CV001是标准的外部大气环境，其内部参数恒定与时间无关，温度和压力始终为300 K和0.1 MPa。CV002是半圆形下封头控制体，内径为0.188 m，壁厚为0.015 m，该控制体内部是高温熔池和氩气。CV003是压力容器圆柱段控制体，内径为0.188 m，高度为0.4 m，该控制体内部充满氩气。容器壁面由热构件模拟，HTW00201表示下封头半球形容器壁，HTW00301和HTW00401分别表示容器垂直圆柱壁面和顶部法兰。

图5 FOREVER实验的数值模拟Fig.5 Numerical simulation of the FOREVER experiment

下封头节点划分如图5所示，下封头的底部为0°，水平位置为90°，壁面沿角度方向划分为9个径向段，沿壁厚方向划分为11个节点，细节图以第8环为例，展示了沿壁厚节点与温度的对应关系，TLH801表示第8环外壁温，TLH811表示第8环内壁温。

2.3 实验过程和模拟

FOREVER实验是在压力为0.1 MPa的环境下进行，主要包括4个阶段：初始浇注高温熔体阶段、初始加热阶段、升温和恒温保持阶段。详细实验过程如下：

1）初始浇注高温熔体阶段：实验是在1 473～1 573 K下向容器内部浇注高温熔体（CaO+B2O3），使下封头内壁温度高达900～1 450 K，容器的初始压力为0.1 MPa，由于容器内部温度的升高，压力增加到2.4 MPa；

2）初始加热阶段：0～215 min期间因为电力供应问题实验没有按照计划进行，加热器功率被限制在20 kW而不是预定的38 kW，由于加热功率较低，熔池温度波动性下降；

3）升温和恒温阶段：在215 min之后电力恢复，加热器功率加并保持在预定的38 kW，熔池温度迅速升温达到并保持在预定的1 373～1 473 K，控制体内部压力也保持在预定的2.5 MPa，容器发生蠕变，最终在402 min时容器发生失效。

在215 min之后即升温和恒温阶段为本文模拟分析的重点，在该阶段中容器内外壁温达到并保持在实验预定的温度，直到容器发生破裂后停止实验，实验容器破裂发生在70°～80°区域，容器70º和80º内壁面温度随时间的变化分别如图6（a）和（b）所示，图中黑色散点为实验值，红色实线为本文模拟出的内壁面节点温度，可以看出，模拟出的内壁温与实验数据基本一致。由于缺少容器70º和80º外壁面温度随时间变化的实验值，本文将对比在300 min达到相对稳态时外壁面随角度变化的温度分布，从图7可以看出，模拟出的内壁温变化过程与实验过程基本一致，内外壁温模拟曲线与实验符合较好。

图6 下封头70º (a)和80º (b)内壁温实验数据与模拟结果（彩图见网络版）Fig.6 Experimental data and simulation results for the 70° (a) and 80º (b) inner wall temperature of the lower head (color online)

图7 300 min时外壁温实验数据与模拟结果Fig.7 Experimental data and simulation results for the external wall temperature at 300 min

压力变化曲线如图8所示，图中黑色散点为实验值，实验中容器内部压力通过气体系统控制，在215 min之后稳定在实验预定的2.5 MPa，在本文模拟中使用ISAA程序的控制函数来模拟出该控制过程，红色实线为数值模型中下封头控制体的压力，可以看出，模拟结果与实验中压力变化符合较好。

图8 FOREVER压力实验数据与模拟结果（彩图见网络版）Fig.8 FOREVER pressure experimental data and simulation results (color online)

3 计算结果与分析

3.1 下封头应力

下封头失效时壁厚中间节点应力沿角度分布如图9所示，下封头大变形模型从机理出发，在变形基础上对下封头壁面进行实时受力分析，得到等效应力即壁厚中间节点应力，该等效应力受内压、熔融物重量等因素影响。

图9 失效时壁厚中间节点应力分布Fig.9 Stress distribution in the middle node of the wall thickness at failure

在下封头失效节点沿壁厚的等效应力如图10所示，第1节点代表外壁面等效应力，第11节点代表内壁面等效应力。从图10可以看出，在下封头失效节点大变形模型计算出的等效应力沿壁厚由外到内逐渐减少，外壁面应力大于内壁面应力这是由于下封头受向外作用的剪力的作用，在下封头经线方向产生较大的经向弯曲作用，在下封头外壁面形成压缩应力，其内壁面引起拉伸应力，由泊松效应的作用，在下封头外壁面产生周向压缩应力，内壁面引起周向拉伸应力。在下封头外壁面上由弯曲引起的周向压缩应力与周向薄膜压缩应力相叠加，构成了下封头的最大应力。

图10 下封头失效处沿壁厚的等效应力Fig.10 Equivalent stress along the wall thickness at the failure position of the lower head

目前，很多模型不能计算沿壁厚上的等效应力分布，简单地将下封头壁厚中间节点应力作为沿壁厚所有节点的应力，比如Larson-Miller模型和原ISAA程序的LHTCM模型，这与下封头真实受力情况并不相符。

3.2 下封头变形

图11为大变形模型预计的下封头整体形变情况，图中黑色实线为实验容器的初始形状，红色实线为大变形模型预计的形变结果，最终呈现经典的鸡蛋状形态，从下封头的形状变化可以发现下封头整体向下发生形变，最大形变位移发生在下封头底部位置。原ISAA程序的LHTCM模型由于缺乏计算变形的模块，不能预测出下封头整体变形情况。

图11 下封头变形前后形状（彩图见网络版）Fig.11 Shape of the lower head before and after deformation(color online)

底部垂直位移如图12所示，底部最大垂直位移实验结果为2.2 cm。SLM（Surface Loads Mapping）模型和VLM（Volume Loads Mapping）模型［16］从194.5 min开始模拟下封头形变，并在194.5 min增加了初始位移3.55 mm，这两种模型预测的最大总变形为1.8 cm。

图12 下封头底部伸长量实验数据与模拟结果Fig.12 Experimental data and simulation results for elongation at the bottom of the lower head

本文下封头大变形模型是从实验初始时间开始模拟下封头变形，变化趋势与实验基本一致，实验初期下封头由于受到高温熔体的冲击发生形变。在初始加热阶段由于加热功率较小，下封头温度波动性下降，不足以让下封头发生较为明显的形变。在加热功率恢复后，下封头开始发生明显形变，下封头变形随着时间而增加，直至失效。下封头大变形模拟结果为2.73 cm，略高于实验值，原因可能是高温蠕变受压力和温度影响较大［17］，大变形模型对高温较为敏感，对低温敏感度较低，在初始浇注熔体阶段，下封头壁温较低，大变形模型模拟出的底部垂直位移相对实验数据偏小。在最终的升温和恒温阶段，下封头壁温升至较高温度时，模拟出的底部垂直位移偏大且位移速率偏大。但整体而言，大变形模型模拟出的形变情况与实验结果基本符合。

3.3 失效时间及位置

RPV下封头的损伤主要是由高温蠕变造成的，遵循Kachanov提出的蠕变损伤准则，在大变形模型中选用Kachanov蠕变损伤准则，并以失效分数的形式来表示，当失效分数达到1（即100%）时，认为下封头该位置失效。下封头各个位置的失效分数变化如图13所示，第8环的失效分数最先达到1（即100%），其他环均未达到失效标准，即下封头破口位置发生在75°～85°。

图13 下封头各个位置的失效分数Fig.13 Failure fraction at each position of the lower head

本文将下封头大变形模型模拟结果与实验、SLM模型、VLM模型进行结果比较，如表1所示，实验测量的容器失效时间为402 min，破口位置大约在θ=75°的位置。大变形模型预计的失效时间为394.33 min，误差为1.9%，与实验结果吻合较好，大变形模型预测的破口位置都是75°～85°之间与实验结果一致。SLM模型预测失效时间相对较早，VLM模型预测失效时间与实验较为相近，但如果考虑到容器失效时间非常长（约402 min）和各种不确定性，VLM模型与大变形模型预测的失效时间之间的差异可以忽略不计，都能准确预测下封头失效时间。

表1 结果对比Table 1 Comparison of results

4 结语

本文为了解决ISAA程序LHTCM模型简化薄膜应力模型十分简单、需要更加机理的应力模型、缺乏计算变形模块等问题，基于Timoshenko板壳理论、Nortron蠕变定律和大变形塑性理论开发了机理模型—下封头大变形模型，并将该模型集成到一体化严重事故分析程序ISAA中，并在FOREVER-EC2实验中进行应用研究，验证了该模型在压水堆堆芯熔化严重事故中能准确预测RPV整体形变以及失效时间和位置等情况。主要研究结果总结如下：

1）下封头大变形模型从机理出发，在变形基础上对下封头壁面进行实时受力分析，得到等效应力即壁厚中间节点应力，并且相对于Larson-Miller模型和原ISAA程序的LHTCM模型，能计算沿壁厚的应力分布，外壁面应力大于内壁面应力。

2）当发生严重事故执行IVR策略时大变形模型能准确地预测出RPV下封头的失效时间、破口位置和下封头的变形情况。

3）相比较于缺乏计算变形模块的LHTCM模型，下封头大变形能准确预测出下封头整体变形情况。大变形模型对高温较为敏感，对低温敏感性较低，整体而言，大变形模型模拟出的形变结果与实验结果基本符合，变化趋势基本一致。

当发生堆芯熔化的严重事故时在执行IVR期间，下封头大变形模型从机理出发，能够对下封头所受应力、失效时间、位置和整体形变等关键信息进行有效评估，解决了LHTCM模型的简化薄膜模型十分简单、缺乏更加机理的应力模型、不能预测下封头整体变形等问题。

作者贡献声明翟润泽负责软件开发、稿件撰写；张斌负责调研、方法论；杨皓负责软件开发；高鹏程负责调研、修改文章；唐绍伟负责技术支持；单建强负责修改文章。