0.   引　言

随着社会和经济的高速发展，能源的需求越来越强烈，传统的煤、石油等自然资源面临枯竭，同时环境将受到污染。因此，探索经济、绿色和可持续开发的清洁能源对经济的稳定发展和环境的持续保护具有重要意义。其中，核电作为一种绿色和可持续开发的清洁能源，其长期安全地发展不仅能很好地解决能源枯竭和环境污染问题，而且还能强有力地体现一个国家的雄厚实力。核电发展中，安全是核电发展的关键。因此，对核电安全理论和技术的研究具有重要的价值。在反应堆发生冷却剂丧失事故（LOCA）后，堆芯热量由于长时间得不到转移而使堆芯熔融，然后掉落到反应堆压力容器（RPV）下封头上，进而对RPV下封头产生烧蚀现象。对此，众多科研工作者对RPV在堆芯熔融 物条件下结构的稳定性和完整性进行了广泛地研究。例如，赵等^[1]采用数值方法分析了1000 MW级压水堆安全壳的压力和温度分布情况。Anderson等^[2]采用实验方法分析了AP600反应堆RPV在假定事故条件的热的转移问题。Zhang等^[3]研究了RPV在应急停堆事故中的失效准则。Theofanous等^[4]和Guan等^[5]基于2层熔池结构，分析了RPV下封头内壁面的热负荷、对流和辐射热流分布等。Fu等^[6]和Tsai等^[7]采用2层熔池结构分析了在堆芯熔融事故下熔池的形成过程和热分布问题。在有限元方法的框架下，张等^[8]基于生死单元法提出了相应的移动边界分析技术来模拟堆芯熔融物对RPV内壁的烧蚀过程。Zhan等^[9]基于ANSYS软件，采用2层结构熔池分析了RPV内壁厚度的变化、温度和热应力的分布等。总的来说，以上方法能较好地研究在堆芯熔融事故下RPV内壁的烧蚀过程、温度和热应力分布情况，然而却很少有研究同时考虑了在堆芯熔融事故下RPV内壁的烧蚀和RPV在内压作用下的断裂这一耦合过程。

在2000年，Silling^[10-11]提出了一种采用积分方程来作为控制方程的非局部方法，即近场动力学方法。近场动力学方法具有分析断裂等间断不连续问题的自然优势，且自其提出以来已经得到了深入发展和广泛应用。例如，Huang等^[12]利用近场动力学方法研究了脆性材料的变形和断裂问题；Bobaru ^[13]提出了瞬态传热问题的近场动力学模型；Li等^[14]发展了拉压不同模量结构静力问题的近场动力学模型，并成功用于桁架和张拉整体结构的大变形分析；Zhang等^[15]发展了饱和多孔介质液固耦合问题的近场动力学模型，并用其分析了饱和多孔介质固结和动力响应等。因此，基于近场动力学方法在瞬态传热问题和断裂问题中的良好应用，本文采用近场动力学方法研究在堆芯熔融事故下RPV内壁的烧蚀和RPV在内压作用下的断裂过程。另外，在堆芯熔融物对RPV内壁的烧蚀过程中，由于RPV厚度随着烧蚀的过程会逐步减小，即RPV内壁的几何形状会随时间发生变化，因此，本文还在近场动力学方法的框架下，提出了一种简便和高效的移动边界模型来分析这种变边界问题。

1.   近场动力学热力耦合问题控制方程

根据键型近场动力学基本理论^[11]，考虑结构在外力作用下的变形和热膨胀效应，任意物质点$ {{{\boldsymbol{x}}}^i} $的近场动力学运动方程表示为：

式中，$\rho $为材料的密度；t为时间；${\ddot {\boldsymbol{u}}}$为加速度；$ {{{\boldsymbol{x}}}^i} $和$ {{{\boldsymbol{x}}}^j} $分别为物质点$i$和$j$的初始位置；ξ为物质点在初始构型中的相对位置；${{\boldsymbol{\xi}} = }{{{\boldsymbol{x}}}^j}{ - }{{{\boldsymbol{x}}}^i}$；η为物质点在当前构型中的相对位置，${{\boldsymbol{\eta}} = }{{{\boldsymbol{x}}}^j}{ + }{{{\boldsymbol{u}}}^j}{ - }{{{\boldsymbol{x}}}^i}{ - }{{{\boldsymbol{u}}}^i}$；${{{\boldsymbol{u}}}^i}$和${{{\boldsymbol{u}}}^j}$分别为$i$和$j$的位移；${\mathcal{H}_{{{x}^i}}}$为$ {{{\boldsymbol{x}}}^i} $的邻域；${V_{{{{\boldsymbol{x}}}^j}}}$为$ {{{\boldsymbol{x}}}^j} $的体积；$ {{{\boldsymbol{b}}}^i} $为作用于$ {{{\boldsymbol{x}}}^i} $上的单位体积力；对势力$ {{\boldsymbol{f}}}\left( {{{\boldsymbol{\xi}} ,{\boldsymbol{\eta}} },\theta ,t} \right) $为$ {{{\boldsymbol{x}}}^j} $作用在$ {{{\boldsymbol{x}}}^i} $上的单位体积密度力。

式中，$ c $为材料的微观弹性模量；$\left\| {\;} \right\|$为向量的模；$ \omega \left( {{\boldsymbol{\xi}} } \right) $为权函数；$ {\alpha _\Theta } $为热膨胀系数；$ s $为键的伸长率；$ \tilde T $为$ {{{\boldsymbol{x}}}^i} $和$ {{{\boldsymbol{x}}}^j} $的平均温度。

另外，为了表征键的断裂情况，在$ {{\boldsymbol{f}}}\left( {{{\boldsymbol{\xi}} ,{\boldsymbol{\eta}} },\theta ,t} \right) $中引入为键断裂状态的标量场$\mu \left( {{{\boldsymbol{\xi}} },t} \right)$如下：

式中，$ {s_0} $为键断裂的临界伸长率。

考虑结构中热传导问题，假设裂纹两侧物质点之间的热传导系数为$ 0 $，则任意物质点$ {{{\boldsymbol{x}}}^i} $的近场动力学传热方程为：

式中，$ {c_\Theta } $为材料的热容；$ K\left[ \delta \right] $为材料的微观热传导系数；$ \delta $为物质点的邻域半径大小；$ {Q^i} $为热源。值得注意的是，当$ {{{\boldsymbol{x}}}^i} $的邻域趋近于0时，式（6）可以演化到经典的连续介质模型下为热传导问题的微分方程。因此，在这种条件下，可以推导得到式(6)中材料的微观热传导系数^[15]，即

式中，$k$为材料在连续介质模型中的热传导系数；1D为一维；2D为二维。式（7）和式（8）分别为考虑权函数为常数和线性分布时得到的一维和二维问题材料的微观热传导系数。

2.   移动边界模型

在高温条件下的烧蚀现象是金属材料和结构破坏失效的重要原因之一。金属材料和结构在施加外部热源的作用下，当局部任意一个位置的温度达到材料的熔点时，材料和结构将发生熔融，即烧蚀现象。可以发现，金属材料和结构在高温烧蚀的过程中，其外貌和几何形状都随时间发生变化。因此，在采用数值方法模拟金属材料和结构烧蚀过程时，需要采用合理和高效的移动边界分析技术。接下来，在近场动力学方法框架下，将提出一种简便和高效的移动边界模型。图1所示为近场动力学的移动边界模型。在该移动边界模型中，首先将初始结构离散成若干大小一致的物质点。其次，如图1所示，将所有的物质点划分成3类，即：蓝色的内部物质点（IMP），紫色的边界上物质点（BMP）和红色的烧蚀物质点（AMP）。在初始的时候，结构中只存在IMP和BMP 2种物质点。在模拟烧蚀的过程中，随着外部热源的热量流入结构与结构中热的传导，当结构中任意物质点的温度达到材料熔点时，相应的物质点标记为AMP。同时，根据需要更新物质点IMP和BMP的类别。另外，为了表征任意${{{\boldsymbol{x}}}^i}$的烧蚀状态，引入一个标量场$\varphi \left( {{{{\boldsymbol{x}}}^i},t} \right)$，表示为：

图  1  近场动力学移动边界模型

Figure  1.  Moving Boundary Model of Peridynamics

下载: 全尺寸图片 幻灯片

式中，$\hat T$为材料融化时的临界温度；$\varphi \left( {{{{\boldsymbol{x}}}^i},t} \right) = 1$表示${{{\boldsymbol{x}}}^i}$未被烧蚀；$\varphi \left( {{{{\boldsymbol{x}}}^i},t} \right) = 0$表示${{{\boldsymbol{x}}}^i}$已经被烧蚀。

进一步地，将表征烧蚀状态的$\varphi $引入到${{{\boldsymbol{x}}}^i}$的运动方程和传热方程中，因此，式（1）和式（6）可以改写为：

或者：

式中，$\varphi \left( {{{{\boldsymbol{x}}}^i},{{{\boldsymbol{x}}}^j},t} \right)$为${{{\boldsymbol{x}}}^i}$和${{{\boldsymbol{x}}}^j}$相互作用阈值。

因此，式（10）和式（11）为考虑烧蚀和断裂过程热力耦合问题的控制方程。通过引入表征键断裂状态的标量场$\mu $和表征物质点烧蚀状态的$\varphi $，使得在计算过程中不需要不断地更新计算域，从而提高计算效率。同时，可以有效地结合近场动力学方法处理断裂问题和移动边界问题的优势，并将其应用到RPV下封头构件在高温和压力条件下烧蚀和断裂过程的数值模拟中。

3.   算法流程

为了提高计算效率，采用显式的弱耦合求解策略，即结构变形和热传导交替求解的方法。首先，将整个求解域在空间上离散成$N$个大小一致的物质点，同时将总的分析时间${t_{{\text{total}}}}$划分为${N_{{\text{step}}}} = \dfrac{{{t_{{\text{total}}}}}}{{\Delta t}}$个时间步（$\Delta t$为时间步长）。假设整个求解域内所有物质点在$ t $时刻的位移、速度和温度等信息是已知的，则${{{\boldsymbol{x}}}^i}$在$t + \Delta t$时刻的值可以由下面的式子计算得到，即：

式中，${N_{{{{\boldsymbol{x}}}^i}}}$为${{{\boldsymbol{x}}}^i}$邻域中物质点的个数。然后，采用差分方法对位移和温度进行求解，物质点${{{\boldsymbol{x}}}^i}$在$t + \Delta t$时刻的速度、位移和温度可以表示为：

4.   RPV下封头高温烧蚀与断裂过程数值模拟

采用上述近场动力学方法模拟AP600的RPV下封头在堆芯熔融事故下的烧蚀和断裂破坏过程。堆芯熔融物在RPV下封头形成的2层熔池简化结构见图2。其中，RPV内壁的半径$R = 2{\text{ m}}$，壁厚${d_{\text{V}}} = 0.2{\text{ m}}$。铀氧化物熔池体积${V_{\text{U}}} = 8.6{\text{ }}{{\text{m}}^3}$，金属熔池体积${V_{\text{M}}} = 8.6{\text{ }}{{\text{m}}^3}$。RPV材料的密度为6890 kg·m⁻³、热容为740 J/( kg·K)和传热系数为25.5 W/( m·K)。在瞬态分析中，RPV的初始温度为373 K，在RPV内壁上施加的随几何角度变化的热流值如图3所示。RPV外壁为沸腾水环境，外壁的温度保持为$373{\text{ K}}$不变。从图4可以看出，随着时间的增加，RPV的最薄区域的厚度越来越小。图5所示的RPV的烧蚀形状与文献[9]的结果一致，从而验证了本文方法模拟烧蚀过程的有效性和正确性。

图  2  RPV下封头高温烧蚀计算模型

θ—与中轴线夹角；H—熔融物厚度；下标M—金属熔池；下标U—铀氧化物熔池；下标V—RPV

Figure  2.  Computational Model of RPV Lower Head for Thermal Ablation

下载: 全尺寸图片 幻灯片

图  3  施加在RPV内壁上的热流密度随几何角度的变化

Figure  3.  Variation of Heat Flux Imposed on the RPV Inner Wall with the Angles

下载: 全尺寸图片 幻灯片

图  4  热烧蚀作用下RPV移动的内壁轮廓

Figure  4.  Profile of the Moving Inner Wall of RPV Caused by Thermal Ablation

下载: 全尺寸图片 幻灯片

图  5  不同时刻RPV温度和烧蚀形状分布图

Figure  5.  Temperature Contours and Ablated Shapes of RPV at Different Times

下载: 全尺寸图片 幻灯片

图6给出了压力容器下封头烧蚀区沿壁厚方向的等效Mises应力分布情况。有无内压情况下曲线无明显误差。如图所示路径的应力分布特点及其与相应位置材料屈服极限的比较，可以看出整个壁厚方向上除了中间小部分区域外，其余区域处的等效Mises应力水平都远远高于材料在相应温度下的屈服强度。这必定会造成压力容器下封头这些区域的材料屈服。

图  6  弹性分析条件下截面路径上等效Mises应力和温度分布

Figure  6.  Curves of Mises Stresses and Temperature Distribution for the Given Path under Elastic Condition

下载: 全尺寸图片 幻灯片

另外，假定在$t = 2000{\text{ s}}$时，RPV外部流体的散热量和内部堆芯熔融产生的热量达到平衡条件（即满足临界热通量边界条件），此时RPV将不再产生烧蚀现象，且各部分的温度达到稳态分布。在此基础上，进一步考虑被烧蚀后的RPV在内压作用下的断裂破坏过程。图7给出了裂纹扩展路径云图，图7中损伤值表示物质点的邻域内断裂的键的数量与总的键的数量的比值。从图7中可以看出，在内压作用下，结构将从被烧蚀后的最薄处开始断裂。

图  7  烧蚀后裂纹扩展路径云图

Figure  7.  Contour of Crack Paths After Ablation

下载: 全尺寸图片 幻灯片

5.   结　论

本文基于近场动力学基本理论，首先给出了适用于模拟RPV烧蚀及断裂过程的积分列式。其次，为了表征烧蚀过程中计算模型几何边界的变化情况，在近场动力学框架下，通过引入一个标量场来直接描述物质点的烧蚀状态，并提出了一种简便和高效的移动边界模型。基于这个移动边界模型，在计算的过程中不需要重复地更新计算域，也不需要像有限元方法需要不断地更新刚度矩阵，从而使得计算效率大大提高。另外，本文给出了热力耦合问题的显式的弱耦合求解策略，并给出了简要的计算流程。最后，本文通过堆芯熔融事故下RPV烧蚀及断裂的算例来验证了提出的方法的有效性和正确性。同时，计算结果中路径的应力分布可以看出整个壁厚方向上除了中间小部分区域外，其余区域处的等效应力水平都远远高于材料在相应温度下的屈服强度，因此在堆芯熔融事故下RPV还存在塑性变形、高温蠕变变形及损伤失效等复杂力学行为。