0.   引　言

与第二代核电相比，第三代核电具有更好的经济型与安全性。我国第三代核电“华龙一号”设置了堆腔注水冷却系统（CIS），可在严重事故情况下对压力容器外壁进行冷却，从而实现堆芯熔融物堆内滞留（IVR），防止放射性物质泄漏。在这种极端的情况下，反应堆堆芯熔融物可能坍塌到反应堆压力容器（RPV）下封头上，使RPV内壁承受上千摄氏度的高温，而下封头外部又存在CIS注入的大量冷却水。一方面，高温会使下封头材料的力学性能急速下降；另一方面，下封头内、外壁较高的温度差会导致极高的热应力，这2个原因的共同作用会引发下封头材料出现严重的蠕变损伤^[1]，给RPV的完整性带来极大挑战。为了防止RPV在严重事故情况下发生蠕变失效，保证一回路压力边界的完整性，在设计时必须要考虑其下封头的高温蠕变行为。

针对堆芯熔化的严重事故，美国、德国和瑞典等国家先后采用RPV缩比模型开展了一系列的实验研究^[2-3]，为了解严重事故的发生过程提供了大量的实验数据，也为研究事故后的应急措施提供了依据。为了准确分析严重事故工况下材料的高温蠕变行为，需要不断完善材料的高温性能数据。20世纪90年代起，德国、法国、美国等国家陆续开展了一系列材料的高温蠕变性能测试。德国卡尔斯鲁厄核能研究中心（KFK）的Müller和Kuhn对德国RPV用20MnMoNi55钢开展了高温蠕变试验，试验温度范围为700~950℃。美国爱达荷州国家实验室（INL）的Thinnes等针对美国RPV下封头用SA533B1钢进行了拉伸和蠕变试验^[4]，其温度范围为627~1200℃，得到了材料在高温下的屈服强度、抗拉强度和蠕变试验数据。法国原子能委员会的Sainte等对RCC-M标准中的RPV用16MND5钢的高温蠕变性能进行了测试^[5]，温度范围达到600~1300℃，获得了高温下16MND5材料的拉伸性能和蠕变性能。

而国内针对RPV材料蠕变行为的研究主要集中在蠕变本构理论和数值模拟分析上，谢志刚^[6]等对RPV用508-Ⅲ钢的蠕变损伤进行了试验研究，并提出了基于微观损伤修正的蠕变本构方程；罗娟^[7]等采用有限元方法对严重事故工况下RPV的蠕变行为进行了模拟，给出了RPV在严重事故工况下的主要失效模式。目前，国内针对RPV用钢的高温性能测试及其蠕变本构的研究尚处于起步阶段。

本文针对国产RPV用16MND5钢开展了高温蠕变性能试验测试，并基于应变强化基本蠕变本构和延性耗竭蠕变损伤模型建立16MND5材料的蠕变损伤本构模型，可为RPV蠕变行为的分析和模拟提供基础。

1.   国产16MND5钢蠕变性能测试

为了获得国产16MND5钢的蠕变特性，本文基于GB/T 2039—2012 金属材料单轴拉伸蠕变试验方法^[8]开展了试验测试，试验件为圆柱形试件，试件材料为基于RCC-M规范设计的RPV用16MND5钢，试验温度范围为600~900℃。由于该材料在600℃以上高温时屈服强度下降的很快，如表1所示，故不同温度下的试验应力水平应根据当前温度的屈服强度制定。

表  1  国产RPV用16MND5钢在600℃以上的屈服强度

Table  1.  Yield Strength of 16MND5 Steel for Domestic RPV above 600℃

温度/℃	600	700	800	900	1000
屈服强度/MPa	243.3	99.3	50.7	36.0	23.5

下载: 导出CSV 

| 显示表格

图1为不同温度下的蠕变试验结果，这些蠕变试验数据为下文的蠕变本构模型拟合提供了基础。

图  1  国产RPV用16MND5钢的高温蠕变试验结果

Figure  1.  Creep Test Results of 16MND5 Steel for Domestic RPV at High Temperature

下载: 全尺寸图片 幻灯片

2.   蠕变损伤本构模型

2.1   经典蠕变本构模型

蠕变效应与温度、应力及时间相关，在高温及恒定载荷条件下，材料的蠕变变形量随时间增加而增长。对于一般材料，其蠕变应变可以表示为应力、温度以及时间的函数：

式中，$ {\varepsilon _{{\text{cr}}}} $为蠕变应变；$ \sigma $为应力；$ T $为温度；$ t $为时间；$ {f_1}\left( \sigma \right) $、$ {f_2}\left( T \right) $及$ {f_3}\left( t \right) $分别为应力、温度及时间的函数，与材料特性密切相关。

通过选择合适的应力函数、温度函数以及时间函数，可将材料的蠕变应变表达为：

式中，$ B $、$ m $及$ n $均为材料参数，其中$ n $也被称为稳态蠕变指数。式（2）就是经典的蠕变本构模型——Norton模型。

材料的蠕变失效过程可分为3个阶段：第一阶段（减速蠕变阶段）、第二阶段（稳态蠕变阶段）及第三阶段（加速蠕变阶段），在加载瞬间产生初始弹性应变之后，3个阶段会随着时间增加逐渐发生。典型的蠕变应变随时间的变化曲线如图2所示。由于经典蠕变本构模型如Norton模型，以及在此基础上提出的应变强化模型、时间强化模型等，均未考虑蠕变损伤的情况，往往只能较为精确地描述前2个阶段的蠕变变形，无法准确描述蠕变的第三阶段。因此，需要在经典蠕变本构模型的基础上合理地耦合蠕变损伤模型，实现对全部3个阶段蠕变过程的合理描述。

图  2  典型的蠕变应变-时间曲线

$ {t_{\text{f}}} $—发生蠕变断裂的时间；$ {\varepsilon _{\text{f}}} $—发生蠕变断裂时的应变，即蠕变延性　　　　　　　

Figure  2.  Typical Creep Strain-Time Curve

下载: 全尺寸图片 幻灯片

2.2   耦合损伤的蠕变本构模型

蠕变损伤力学是介于微观力学、连续介质力学和材料科学之间的一门交叉学科，其建立了高温材料性能劣化和材料宏观破坏之间的联系。基于蠕变损伤力学，可以建立描述材料全寿命蠕变行为的耦合蠕变损伤的蠕变本构模型。

对于一般材料，在式（1）的基础上，考虑损伤的蠕变应变率$\left( \dfrac{{{\text{d}}{\varepsilon _{{\text{cr}}}}}}{{{\text{d}}t}}\right) $可以表示为$ \sigma $、$ t $、$ T $以及反映材料劣化程度的一系列损伤变量（$ {\omega _1} $，$ {\omega _2}, \cdots, {\omega _n} $）的函数，其数学关系式为：

式中，$ {f_4}\left( {{\omega _1},{\omega _2},\cdots,{\omega _n}} \right) $表示蠕变应变率与损伤变量之间的函数关系式。

通过对不同材料的蠕变-时间曲线演化形式的观察，结合不同函数的曲线形状与功能，最后通过选取合适的应力函数、温度函数、时间函数以及蠕变损伤函数，就可以得到耦合损伤的蠕变本构模型。

2.3   本文蠕变损伤模型的主控方程

16MND5钢的蠕变实验数据表明，该材料具有明显的蠕变三阶段特性。因此，本文采用应变强化的蠕变本构模型描述材料的第一、第二阶段蠕变变形行为，通过耦合基于延性耗竭理论的蠕变损伤演化方程来描述材料的蠕变损伤，同时加入多轴应力影响因子，以建立16MND5钢的多轴蠕变损伤本构模型。具体方法如下：

应变强化蠕变本构的基本表达式为：

式中，$ {\dot \varepsilon _{{\text{cr}}}} $为蠕变应变率；A为材料参数；$ {\sigma _{{\text{eq}}}} $为等效应力；$ \sigma _{{\text{eq}}}^{n - 1} $为与n有关的等效应力；S为偏应力。

在式（4）的基础上耦合蠕变损伤，则相应的蠕变损伤本构模型中的${\dot \varepsilon _{{\rm{cr}}}}$可表示为：

式中，$ \omega $为材料蠕变损伤变量。

通过对16MND5钢的蠕变曲线分析发现，其蠕变变形具有很强的应力水平依赖性，其材料参数也应该具有相应的应力水平依赖性，因此本文对基于应变强化的蠕变本构模型做了进一步地修正，依据经验和试算提出了具有应力水平依赖性的n，其表达式为：

式中，$ {n_1} $、$ {n_2} $、$ {n_3} $为材料参数，可通过描述不同载荷水平的蠕变曲线拟合获得。

16MND5钢的蠕变损伤变量通过延性耗竭理论获得，该理论假设当蠕变裂纹尖端的局部蠕变应变累积达到临界蠕变延性时，裂纹即往前扩展，此时，损伤变量达到临界值。因此，蠕变损伤变量的演化模型为：

式中，$ \dot \omega $为蠕变损伤率；$ \varepsilon _{\text{f}}^{\text{*}} $为多轴蠕变延性，即发生蠕变断裂时的多轴蠕变应变。

由于裂尖区域应力状态的复杂性，通常多轴蠕变延性与单轴蠕变延性有所不同，两者的比值$ {{\varepsilon _{\text{f}}^*} / {{\varepsilon _{\text{f}}}}} $称为多轴延性因子。多轴延性因子是架构起单轴和多轴蠕变延性非常重要的桥梁，也是实现以延性耗竭理论为基础的单轴蠕变损伤模型向多轴拓展的重要方法。

为了描述多轴延性因子，在孔洞长大理论、孔洞形核理论及试验数据拟合的基础上，已发展了多种多轴蠕变延性因子模型（MCDF）。基于幂率控制蠕变孔洞长大理论，Wen和Tu等^[9]以指数方程的形式定义了MCDF，其表达式为：

式中，$ {\sigma _{\text{h}}} $为静水压力，其与等效应力的比值$ {{{\sigma _{\text{h}}}}/{{\sigma _{{\text{eq}}}}}} $即为应力三轴度。

该MCDF模型很好地解决了蠕变寿命计算值过于保守的缺陷，得到了与实验更加接近的结果。基于此，本文的损伤模型选取Wen和Tu提出的MCDF模型中的多轴延性因子。

当恒定温度时，${\varepsilon _{\rm{f}}}$与${\dot \varepsilon _{{\text{cr}}}}$的关系为：

式中，$ d $为蠕变延性乘子；$\beta$为温度相关的材料参数；$\dot \varepsilon _{{\rm{cr}}}^\beta$为恒定温度下的蠕变应变率。d和β可通过实验数据拟合获得；

对于本文所研究的16MND5钢，与单轴蠕变延性相关的参数$ d $也具有一定的载荷水平依赖性，因此采用与$ n $的修正相类似的方法，对参数$ d $也进行了修正：

式中，$ {d_1} $、$ {d_2} $、$ {d_3} $均为材料参数，可通过曲线拟合获得。

联立式（5）~式（10），即可构成本文蠕变损伤本构模型的主控方程。

3.   16MND5钢的蠕变参数确定及验证

3.1   蠕变损伤本构模型的参数确定

本文的蠕变损伤本构模型涉及多个参数，在进行验证之前，需要针对16MND5钢的每个蠕变参数进行拟合确定。参数确定的过程如下：

（1）针对某一温度水平下的蠕变曲线，由于蠕变第一阶段损伤很小，因此假定此时损伤为零，用式（5）拟合第一阶段蠕变曲线，获得材料参数$ A $、$ m $和$ n $。在此过程中，$ A $、$ m $为确定值，不随应力水平变化，而允许$ n $随应力水平变化。

以800℃下16MND5钢的蠕变试验曲线（图1c）为例，在20、15、10 MPa载荷水平下第一阶段的蠕变曲线拟合结果如图3所示。根据图1可获得材料参数$ A = 5.73 \times {10^{{{ - }}22}} $、$ m = {{ - }}0.865 $以及与应力水平相关的$ n = 12.81 $（$ {\sigma _{{\text{eq}}}} $=10 MPa）、$ n = 11.26 $（$ {\sigma _{{\text{eq}}}} $=15 MPa）、$ n = 10.54 $（$ {\sigma _{{\text{eq}}}} $=20 MPa）。

图  3  800℃下蠕变试验第一阶段曲线拟合

Figure  3.  Curve Fitting of the First Stage of Creep Test at 800℃　　　　　　

下载: 全尺寸图片 幻灯片

（2）根据第一步的结果，绘制$n \text{-} {\sigma _{{\rm{eq}}}}$曲线，然后采用式（6）拟合曲线即可得到材料参数n₁、n₂、n₃。

同样以16MND5钢在800℃下的蠕变曲线为例，根据步骤（1）拟合得到的3个应力水平下的$ n $，绘制$n \text{-} {\sigma _{{\rm{eq}}}}$曲线如图4所示，再采用式（6）拟合曲线，可以得到$ {n_1} = 9.917 $、$ {n_2} = 13.487 $、$ {n_3} = $$ {{ - }}0.154 $。

图  4  根据不同应力水平下的n值绘制的n-σ_eq曲线

Figure  4.  n-σ_eq Curve Drawn by n Values at Different Stress Levels　　　　　　

下载: 全尺寸图片 幻灯片

（3）在前两步拟合结果的基础上，根据金属材料特性及经验，取相同的$ \beta $，再针对不同应力水平拟合对应的$ d $即可绘制$d \text{-} {\sigma _{{\rm{eq}}}}$曲线，从而得到$ {d_1} $、$ {d_2} $、$ {d_3} $。

同样以16MND5钢在800℃下的蠕变曲线为例，取金属材料的一般值$ \beta = 0.45 $，可以得到不同应力水平下的$ d $值：$ d = 29972 $（$ {\sigma _{{\text{eq}}}} $=10 MPa）、$ d = $$ 14985 $（$ {\sigma _{{\text{eq}}}} $=15 MPa）、$ d = 9793.3 $（$ {\sigma _{{\text{eq}}}} $=20 MPa），绘制d-σ_eq曲线如图5所示，再用式（10）拟合曲线，可得$ {d_1} = 7040.8 $、$ {d_2} = 191043.5 $、$ {d_3} = {{ - }}0.212 $。

图  5  根据不同应力水平下的d值绘制的d-σ_eq曲线

Figure  5.  d-σ_eq Curve Drawn by d Values at Different Stress Levels　　　　　

下载: 全尺寸图片 幻灯片

基于第1节中16MND5钢在不同温度水平下的实验数据，重复上述一系列步骤，拟合得到了适用于16MND5钢的蠕变损伤本构模型参数，600、700、800、900℃ 4个温度水平下的蠕变损伤本构模型参数值见表2。

表  2  不同温度下16MND5钢的蠕变损伤本构模型参数

Table  2.  Creep Damage Constitutive Model Parameters of 16MND5 Steel at Different Temperatures

温度/℃	蠕变参数
	A	n1	n2	n3	m	d1	d2	d3	β
600	5.73×10−22	6.622	9.362	−0.041	−0.550	−11430.075	44826.022	−0.009	0.36
700	5.73×10−22	8.115	7.833	−0.070	−0.865	655.602	84209.310	−0.081	0.45
800	5.73×10−22	9.917	13.487	−0.154	−0.865	7040.816	191043.500	−0.212	0.45
900	5.73×10−22	10.001	11.464	−0.126	−0.865	158.151	36234.250	−0.183	0.45

下载: 导出CSV 

| 显示表格

3.2   蠕变损伤本构模型的验证

根据表2中参数，对16MND5钢在600、700、800、900℃ 4个温度水平下的蠕变演化曲线进行了数值模拟。模拟采用二次开发后的ANSYS程序，将蠕变损伤本构模型写入蠕变用户子程序。选用SOLID185单元对16MND5材料在不同温度、不同载荷水平下的蠕变曲线进行了模拟，结果如图6所示。

图  6  蠕变损伤本构模型的计算值与试验值的对比

Figure  6.  Comparison between Calculated and Experimental Values of Creep Damage Constitutive Model

下载: 全尺寸图片 幻灯片

从图6中可以看出，本文所发展的蠕变损伤本构模型可以合理描述16MND5钢在不同温度、不同应力水平下的3个阶段的蠕变变形行为。

4.   结　论

本文基于应变强化的基本蠕变本构模型与基于延性耗竭理论的蠕变损伤模型，建立了适用于16MND5钢的蠕变损伤本构模型。同时，立足于16MND5钢蠕变变形的特殊性，引入了反映载荷水平依赖性的稳态蠕变指数演化方程和单轴蠕变延性演化方程，并加入MCDF，发展了新的适用于16MND5钢的蠕变损伤本构模型。最后，对新发展的蠕变损伤进行了ANSYS大型有限元软件的数值实现，通过对实验工况的模拟，验证了本构模型的合理性。