Numerical Simulation of RPV with Complex Mechanical Behaviors under Core-melting Accident
-
摘要: 在反应堆冷却剂丧失事故(LOCA)后,堆芯熔融物将对反应堆压力容器(RPV)下封头产生烧蚀作用,为此模拟和分析在堆芯熔融事故下RPV内壁的烧蚀及复杂力学行为对于RPV的设计、事故预防和缓解都具有重要意义。本文首先给出了近场动力学方法分析含裂纹扩展的热力耦合问题的积分列式。其次,在近场动力学方法框架下,提出了一种简便和高效的移动边界模型。其中,通过引一个标量场来直接表征物质点的烧蚀状态,使得在计算过程中不需要不断地更新计算域,从而提高计算效率。最后,采用该方法模拟了堆芯熔融物对RPV内壁的动态烧蚀和RPV在内压作用下的裂纹扩展,计算结果表明在堆芯熔融事故下RPV中同时存在弹性变形和塑性变形,而且还存在损伤断裂等复杂力学行为。Abstract: After the severe loss of coolant accident (LOCA), the lower head of reactor pressure vessel (RPV) will be ablated by the core melt. Therefore, it is of great significance to simulate and analyze the ablation and complex mechanical behavior of the pressure vessel wall under the condition of core melt for the design, accident prevention and mitigation of RPV. This paper firstly presents the peridynamic formulation of coupling thermo-mechanical problem with crack propagation. Then, based on the peridyanmic framework, a simple and efficient moving boundary model is proposed. The abated status of the material points is directly characterized by introducing a scalar field, which makes it unnecessary to constantly update the calculation domain in the calculation process, thus improving the calculation efficiency. Finally, the present method is used to simulate the dynamic ablation of the core melt on the pressure vessel wall and the crack propagation of the pressure vessel under the action of internal pressure. The calculation results show that there are both elastic deformation and plastic deformation in RPV under the core melt accident, and there are also complex mechanical behaviors such as damage and fracture.
-
Key words:
- Peridynamics /
- Coupling thermo-mechanics /
- Ablation /
- Moving boundary /
- Numerical simulation
-
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 {\ddot {\boldsymbol{u}}}\left( {{{{\boldsymbol{x}}}^i},t} \right) = \int_{{\mathcal{H}_{{{{\boldsymbol{x}}}^i}}}} {{{\boldsymbol{f}}}\left( {{{\boldsymbol{\xi}} ,{\boldsymbol{\eta}} },\theta ,t} \right){\text{d}}{V_{{{\boldsymbol{x}}^j}}}} + {{{\boldsymbol{b}}}^i} $$ (1) 式中,$\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} $上的单位体积密度力。
$$ {{\boldsymbol{f}}}\left( {{{\boldsymbol{\xi}} },{{\boldsymbol{\eta}} },\theta ,t} \right) = c\omega \left( {{\boldsymbol{\xi}} } \right)\left( {s - {\alpha _\Theta }\tilde T} \right)\mu \left( {{{\boldsymbol{\xi}} },t} \right)\frac{{{{\boldsymbol{\eta}} - {\boldsymbol{\xi}} }}}{{\left\| {{{\boldsymbol{\eta}} - {\boldsymbol{\xi}} }} \right\|}} $$ (2) $$ {\boldsymbol{s}} = \frac{{\left\| {{\boldsymbol{\eta}} + {\boldsymbol{\xi }}} \right\| - \left\| {\boldsymbol{\xi }} \right\|}}{{\left\| {{\boldsymbol{\eta }} + {\boldsymbol{\xi }}} \right\|}} $$ (3) $$ \tilde T = \frac{{T\left( {{{{\boldsymbol{x}}}^i},t} \right) + T\left( {{{{\boldsymbol{x}}}^j},t} \right)}}{2} $$ (4) 式中,$ 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)$如下:
$$ \mu \left( {{{\boldsymbol{\xi}} },t} \right) = \left\{ {\begin{array}{*{20}{c}} {1}&{s \lt {s_0}} \\ {0}&{s \geqslant {s_0}} \end{array}} \right. $$ (5) 式中,$ {s_0} $为键断裂的临界伸长率。
考虑结构中热传导问题,假设裂纹两侧物质点之间的热传导系数为$ 0 $,则任意物质点$ {{{\boldsymbol{x}}}^i} $的近场动力学传热方程为:
$$ \begin{split} & \rho c_{\Theta}\frac{\partial T\left(\boldsymbol{x}^i,t\right)}{\partial t}=\int_{\mathcal{H}_{\boldsymbol{x}^i}}^{ }K\left[\delta\right]\mu\left(\boldsymbol{\xi},t\right)\omega\left(\boldsymbol{\xi}\right)\times \\ & \qquad\frac{\left[T\left(\boldsymbol{x}^j,t\right)-T\left(\boldsymbol{x}^i,t\right)\right]}{\boldsymbol{\xi}^{ij^2}}\text{d}V_{\boldsymbol{x}^j}+Q^i \end{split} $$ (6) 式中,$ {c_\Theta } $为材料的热容;$ K\left[ \delta \right] $为材料的微观热传导系数;$ \delta $为物质点的邻域半径大小;$ {Q^i} $为热源。值得注意的是,当$ {{{\boldsymbol{x}}}^i} $的邻域趋近于0时,式(6)可以演化到经典的连续介质模型下为热传导问题的微分方程。因此,在这种条件下,可以推导得到式(6)中材料的微观热传导系数[15],即
$$ K\left[ \delta \right] = \left\{ {\begin{array}{*{20}{c}} {\dfrac{k}{\delta },}&{\omega \left( {{\boldsymbol{\xi}} } \right) = 1,}&{1{\text{D}}} \\ {\dfrac{{4k}}{{{\text{π}}{\delta ^2}}}}, &{\omega \left( {{\boldsymbol{\xi}} } \right) = 1,}&{2{\text{D}}} \end{array}} \right. \qquad $$ (7) $$ K\left[ \delta \right] = \left\{ {\begin{array}{*{20}{c}} {\dfrac{k}{\delta },}&{\omega \left( {{\boldsymbol{\xi}} } \right) = 1 - \dfrac{{{\boldsymbol{\xi}} }}{\delta },}&{1{\text{D}}} \\ {\dfrac{{12k}}{{{\text{π}}{\delta ^2}}},}&{\omega \left( {{\boldsymbol{\xi}} } \right) = 1 - \dfrac{{{\boldsymbol{\xi}} }}{\delta },}&{2{\text{D}}} \end{array}} \right. $$ (8) 式中,$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)$,表示为:
$$ \varphi \left( {{{{\boldsymbol{x}}}^i},t} \right) = \left\{ {\begin{array}{*{20}{c}} {1}&{T\left( {{{{\boldsymbol{x}}}^i},t} \right) \lt \hat T} \\ {0}&{T\left( {{{{\boldsymbol{x}}}^i},t} \right) \geqslant \hat T} \end{array}} \right. $$ (9) 式中,$\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)可以改写为:
$$ \begin{split} & \rho\ddot{\boldsymbol{u}}\left(\boldsymbol{x}^i,t\right)=\int_{\mathcal{H}_{x^i}}^{ }c\omega\left(\boldsymbol{\xi}\right)\left(s-\alpha_{\Theta}\tilde{T}\right)\mu\left(\boldsymbol{\xi},t\right)\times \\ & \qquad\varphi\left(\boldsymbol{x}^i,\boldsymbol{x}^j,t\right)\frac{\left(\eta-\xi\right)}{ \left\| \boldsymbol{\eta}-\xi \right\| }\mathrm{d}V_{\boldsymbol{x}^j}+\boldsymbol{b}^i \end{split} $$ (10) $$ \begin{split} & \rho c_{\Theta}\frac{\partial T\left(\boldsymbol{x}^i,t\right)}{\partial t}=\int_{\mathcal{H}_{\boldsymbol{x}^i}}^{ }K\left[\delta\right]\mu\left(\boldsymbol{\xi},t\right)\varphi\left(\boldsymbol{x}^i,\boldsymbol{x}^j,t\right)\omega\left(\boldsymbol{\xi}\right)\times \\ & \qquad\frac{\left[T\left(\boldsymbol{x}^j,t\right)-T\left(\boldsymbol{x}^i,t\right)\right]}{ \left\| \boldsymbol{\xi} \right\| ^2}\text{d}V_{\boldsymbol{x}^j}+Q^i \\[-48pt] \end{split} $$ (11) $$ \varphi \left( {{{{\boldsymbol{x}}}^i},{{{\boldsymbol{x}}}^j},t} \right) = {\text{min}}\left\{ {\varphi \left( {{{{\boldsymbol{x}}}^i},t} \right),\varphi \left( {{{{\boldsymbol{x}}}^j},t} \right)} \right\} $$ (12) 或者:
$$ \varphi \left( {{{{\boldsymbol{x}}}^i},{{{\boldsymbol{x}}}^j},t} \right) = \varphi \left( {{{{\boldsymbol{x}}}^i},t} \right) \cdot \varphi \left( {{{{\boldsymbol{x}}}^j},t} \right) $$ (13) 式中,$\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$时刻的值可以由下面的式子计算得到,即:
$$ \rho {\ddot {\boldsymbol{u}}}_{t + {{\Delta }}t}^i = \mathop \sum \limits_{j = 1}^{{N_{{{{\boldsymbol{x}}}^i}}}} c\varpi \left( {{\boldsymbol{\xi}} } \right)\left( {s - {\alpha _\Theta }\tilde T_t^{ij}} \right)\mu _t^{ij}\varphi _t^{ij}{V_{{{x}^j}}}\frac{{{{\boldsymbol{\eta}} - {\boldsymbol{\xi}} }}}{{\left\| {{{\boldsymbol{\eta}} - {\boldsymbol{\xi}} }} \right\|}} + {{\boldsymbol{b}}}_{t + {{\Delta }}t}^i $$ (14) $$ \rho {c_\Theta }\dot T_{t + {{\Delta }}t}^i = \mathop \sum \limits_{j = 1}^{{N_{{{{\boldsymbol{x}}}^i}}}} K\left[ \delta \right]\omega \left( {{\boldsymbol{\xi}} } \right)\mu _t^{ij}\varphi _t^{ij}{V_{{{{\boldsymbol{x}}}^j}}}\frac{{T_t^j - T_t^i}}{{{{\left\| {{\boldsymbol{\xi}} } \right\|}^2}}} + Q_{t + {{\Delta }}t}^i $$ (15) 式中,${N_{{{{\boldsymbol{x}}}^i}}}$为${{{\boldsymbol{x}}}^i}$邻域中物质点的个数。然后,采用差分方法对位移和温度进行求解,物质点${{{\boldsymbol{x}}}^i}$在$t + \Delta t$时刻的速度、位移和温度可以表示为:
$$ {\dot {\boldsymbol{u}}}_{t + {{\Delta }}t}^i = {\dot {\boldsymbol{u}}}_t^i + {{\Delta }}t{\ddot {\boldsymbol{u}}}_{t + {{\Delta }}t}^i $$ (16) $$ {{\boldsymbol{u}}}_{t + {{\Delta }}t}^i = {{\boldsymbol{u}}}_t^i + {{\Delta }}t{\dot {\boldsymbol{u}}}_{t + {{\Delta }}t}^i $$ (17) $$ {{\boldsymbol{T}}}_{t + {{\Delta }}t}^i = {{\boldsymbol{T}}}_t^i + {{\Delta }}t{\dot {\boldsymbol{T}}}_{t + {{\Delta }}t}^i $$ (18) 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−3、热容为740 J/( kg·K)和传热系数为25.5 W/( m·K)。在瞬态分析中,RPV的初始温度为373 K,在RPV内壁上施加的随几何角度变化的热流值如图3所示。RPV外壁为沸腾水环境,外壁的温度保持为$373{\text{ K}}$不变。从图4可以看出,随着时间的增加,RPV的最薄区域的厚度越来越小。图5所示的RPV的烧蚀形状与文献[9]的结果一致,从而验证了本文方法模拟烧蚀过程的有效性和正确性。
图6给出了压力容器下封头烧蚀区沿壁厚方向的等效Mises应力分布情况。有无内压情况下曲线无明显误差。如图所示路径的应力分布特点及其与相应位置材料屈服极限的比较,可以看出整个壁厚方向上除了中间小部分区域外,其余区域处的等效Mises应力水平都远远高于材料在相应温度下的屈服强度。这必定会造成压力容器下封头这些区域的材料屈服。
另外,假定在$t = 2000{\text{ s}}$时,RPV外部流体的散热量和内部堆芯熔融产生的热量达到平衡条件(即满足临界热通量边界条件),此时RPV将不再产生烧蚀现象,且各部分的温度达到稳态分布。在此基础上,进一步考虑被烧蚀后的RPV在内压作用下的断裂破坏过程。图7给出了裂纹扩展路径云图,图7中损伤值表示物质点的邻域内断裂的键的数量与总的键的数量的比值。从图7中可以看出,在内压作用下,结构将从被烧蚀后的最薄处开始断裂。
5. 结 论
本文基于近场动力学基本理论,首先给出了适用于模拟RPV烧蚀及断裂过程的积分列式。其次,为了表征烧蚀过程中计算模型几何边界的变化情况,在近场动力学框架下,通过引入一个标量场来直接描述物质点的烧蚀状态,并提出了一种简便和高效的移动边界模型。基于这个移动边界模型,在计算的过程中不需要重复地更新计算域,也不需要像有限元方法需要不断地更新刚度矩阵,从而使得计算效率大大提高。另外,本文给出了热力耦合问题的显式的弱耦合求解策略,并给出了简要的计算流程。最后,本文通过堆芯熔融事故下RPV烧蚀及断裂的算例来验证了提出的方法的有效性和正确性。同时,计算结果中路径的应力分布可以看出整个壁厚方向上除了中间小部分区域外,其余区域处的等效应力水平都远远高于材料在相应温度下的屈服强度,因此在堆芯熔融事故下RPV还存在塑性变形、高温蠕变变形及损伤失效等复杂力学行为。
-
-
[1] 赵静,余红星,李锋. 1000MW级压水堆安全壳压力温度计算分析[J]. 核动力工程,2003, 24(5): 409-411,425. doi: 10.3969/j.issn.0258-0926.2003.05.003 [2] ANDERSON M H, HERRANZ L E, CORRADINI M L. Experimental analysis of heat transfer within the AP600 containment under postulated accident conditions[J]. Nuclear Engineering and Design, 1998, 185(2-3): 153-172. doi: 10.1016/S0029-5493(98)00232-5 [3] 张龙飞,房保国,李凤宇. 全厂断电引发的严重事故中反应堆压力容器失效机理研究[J]. 原子能科学技术,2012, 46: 305-308. doi: 10.7538/yzk.2012.46.03.0305 [4] THEOFANOUS T G, LIU C, ADDITION S, ANGELINI S, et al. In-vessel coolability and retention of a core melt[J]. Nuclear Engineering and Design, 1997, 169(1-3): 1-48. doi: 10.1016/S0029-5493(97)00009-5 [5] 关仲华,余红星,江光明. 堆芯熔融物在下腔室内冷却模型研究及缓解集热效应的对策[J]. 核动力工程,2008, 29(5): 72-76. [6] 傅孝良,杨燕华,周卫华,杨晓. CAR 1000的IVR有效性评价堆芯融化及熔池形成过程分析[J]. 核动力工程,2010, 31(5): 102-107. [7] TSAI F J, LEE M, LIU H C. Simulation of the in-vessel retention device heat-removal capability of AP-1000 during a core meltdown accident[J]. Annals of Nuclear Energy, 2017, 99: 455-463. doi: 10.1016/j.anucene.2016.09.052 [8] 张小英,姚婷婷,李志威,等. 堆芯熔融物对压力容器壁面烧蚀过程的数值模拟[J]. 核技术,2015, 38: 020606. doi: 10.11889/j.0253-3219.2015.hjs.38.020606 [9] ZHAN D K, LIU F Y, ZHANG X Y, et al. Ablation and thermal stress analysis of RPV vessel under heating by core melt[J]. Nuclear Engineering and Design, 2018, 330: 550-558. doi: 10.1016/j.nucengdes.2018.02.008 [10] SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces[J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175-209. doi: 10.1016/S0022-5096(99)00029-0 [11] SILLING S A. Linearized theory of peridynamic states[J]. Journal of Elasticity, 2010, 99(1): 85-111. doi: 10.1007/s10659-009-9234-0 [12] HUANG D, LU G D, QIAO P Z. An improved peridynamic approach for quasi-static elastic deformation and brittle fracture analysis[J]. International Journal of Mechanical Sciences, 2015, 94-95: 111-122. doi: 10.1016/j.ijmecsci.2015.02.018 [13] BOBARU F, DUANGPANYA M. The peridynamic formulation for transient heat conduction[J]. International Journal of Heat and Mass Transfer, 2010, 53(19-20): 4047-4059. doi: 10.1016/j.ijheatmasstransfer.2010.05.024 [14] LI H, ZHANG H W, ZHENG Y G, et al. A peridynamic model for the nonlinear static analysis of truss and tensegrity structures[J]. Computational Mechanics, 2016, 57(5): 843-858. doi: 10.1007/s00466-016-1264-4 [15] ZHANG H W, LI H, YE H F, et al. A coupling peridynamic approach for the consolidation and dynamic analysis of saturated porous media[J]. Computational Mechanics, 2019, 64(4): 1097-1113. doi: 10.1007/s00466-019-01695-2 -