螺旋十字燃料棒束热力耦合特性研究

丛腾龙 刘峪洁 郭 辉 肖 瑶 顾汉洋

1（上海交通大学 核科学与工程学院 上海 200240）

2（哈尔滨工程大学 核科学与技术学院 哈尔滨 150009）

燃料组件是反应堆的核心释热部件，是决定反应堆安全性和经济性的关键要素之一，也是反应堆研发和设计的重点。螺旋十字燃料（Helical-Cruciform Fuel，HCF）是一种革新型的燃料组件，HCF元件具体结构如图1所示，十字形横截面结构沿轴向旋扭形成螺旋十字，通过紧密排布，可在高度方向上形成周期性自支撑结构，无须定位格架。同时，十字截面提高了堆芯比换热面积、降低导热距离，且可以通过螺旋结构增强通道间扰流交混，可在维持堆芯安全裕量的前提下提高堆芯功率密度，用于先进紧凑型堆芯的设计。目前俄罗斯已经将具有螺旋翼片的燃料组件用于高温气冷堆中［1］，美国也在积极推进HCF组件的工程应用，图2为美国Lightbridge公司在美国能源部的资助下提出的正方形栅格HCF组件设计方案，与传统17×17圆棒组件在尺寸上兼容，可以提供更高的单盒组件功率输出［2］。

图1 单根螺旋十字燃料及自支撑结构示意图(a) 单根燃料，(b) 自支撑结构Fig.1 Sketch for the single fuel rod and the self-support structure (a) Single fuel rod, (b) Self-support structure

图2 Lightbridge公司设计用于压水堆的HCF组件截面图Fig.2 Cross section of the HCF assembly for PWR designed by Lightbridge

现有针对HCF组件的研究多集中在流动传热方面，包括单相阻力［3-5］、交混［3，6］、对流传热［5］、沸腾传热［6］和临界热流密度［4，6］，证明了HCF组件相比传统带格架圆棒组件在流阻、交混、传热和临界热流密度方面具有优势。然而，HCF组件依靠自身接触支撑形成自定位结构，在反应堆运行过程中，燃料相对冷态会出现明显的热膨胀，棒束间点接触区域可能会出现应力集中，局部过高的应力会使包壳出现塑性变形乃至破裂，危害燃料组件的完整性，MIT的Deng等［7］对辐照条件下单根燃料的应力应变进行了分析，通过引入基体和包壳在辐照作用下的相关热、力学模型，得到金属基体的HCF单棒使役性能，但由于仅考虑单棒，对于多棒之间接触及其对应力应变影响规律的考虑仍不充分。鉴于此，本文建立典型HCF棒束单元的热力耦合分析模型，研究了正常稳瞬态和事故工况下HCF元件的应力应变分布，评价了HCF组件抵御事故的能力，为HCF组件设计提供了参考。

1 螺旋十字燃料热力耦合数学物理模型

相对于冷态停堆工况，反应堆热态运行阶段HCF温度会上升200~400 K，产生热膨胀和热应力，HCF棒束翼片顶部的接触位置出现应力集中，在热应力和局部集中应力共同作用下燃料包壳可能出现大幅度的弹性和塑性变形。此外，十字截面导致的燃料棒内部非均匀径向导热使得燃料棒内部温度差异很大，影响材料的热力学物性，进一步复杂化HCF元件内部应力应变场。本文分别建立燃料元件的温度场和结构场控制方程，给定热对流和约束边界条件，通过热-力耦合求解热力方程获得燃料元件的传热和力学响应。

1.1 基本控制方程

1.1.1 温度场方程

螺旋十字燃料组件的温度场可看作含有内热源的导热问题，包壳和芯块两个固体区域的温度分布可以由导热微分方程确定：

稳态工况下式（1）可以改写为：

式中：ρ为密度，kg‧m-3；c为比热容，J‧（kg·K）-1；T为温度，K；τ为时间，s；λ为导热系数，W‧（m·K）-1；φ为单位时间内单位体积中产生的热量，W‧m-3，对于包壳，该项为零。

1.1.2 结构力学场控制方程

反应堆运行阶段螺旋十字燃料组件的温度会升高，由于燃料间的自支撑结构限制了燃料的热膨胀，燃料内部产生热应力，均质各向同性体满足的基本方程为：

式中：σij为应力张量，Pa；εij为应变张量；ui为位移张量，m；ν为泊松比；E为弹性模量，Pa；δij为kronecker符号；α为热膨胀系数，K-1；ΔT为与参考温度的温差，K；εpl为塑性应变张量，在弹性阶段为零。本文采用von Mises屈服准则判断材料是否进入塑性阶段，有效应力σs可以表示为：

式中：σ1、σ2、σ3分别为第一、第二、第三主应力。当材料进入塑性阶段后，根据Levy-Mises增量理论，εpl可以写成：

1.2 燃料的材料模型

由于U-50wt%Zr合金的热导率高、裂变原子密度大且辐照肿胀率低，螺旋十字燃料采用U-50wt%Zr合金作为芯块材料，同时采用Zr-4合金作为包壳材料。包壳和芯块通过冶金结合的方式连接，中间无气隙。

1.2.1 芯块

在一些严重事故中芯块中心温度可能超过800 K，此时U-50wt%Zr合金会发生相变，导致材料的物性出现变化，但国内外对相变后U-50wt%Zr合金的研究较少，且本文研究工况范围内芯块温度均低于相变点温度（800 K），故不考虑相变对物性的影响。U-50wt%Zr合金的热膨胀系数αf为［8］：

芯块的密度［9］为 9.5×103kg‧m-3，导热系数和比热分别依据 Bauer［9］和 Fedorov［10］的研究结果计算。根据修正后的混合物定律［11］，芯块的弹性模量Ef和泊松比νf可以表示为：

式中：EU为588 K下纯铀的弹性模量，EU=1.6×1011Pa；P为燃料的孔隙率，对于金属芯体的燃料，在未经辐照条件下孔隙率为零；WZ为Zr的质量百分比，本文中WZ=0.5；TmU为纯铀的熔点，TmU=1 405 K；νU为588 K下纯铀的泊松比，νU=0.24。

1.2.2 包壳

当300 K＜T＜1 800 K时，包壳导热系数λc的表达式为［12］：

包 壳 的 密 度［9］为 6.55×103kg ‧m-3，采 用MATPRO［12］物性手册中的模型计算包壳的比热，包壳的热膨胀系数［13］与温度无关，为5.58×10-6K-1，弹性模量Ec和泊松比νc采用Fisher模型［13］：

进入塑性阶段后包壳的应力强化曲线可以拟合成［13］：

式中：σtrue和εtrue为真实应力和真实应变为真实应变率，如

1.3 有限元模型

图2所示的燃料组件包括多根燃料棒，对整盒组件进行力学分析所需的计算代价极大，由圣维南定律可知，组件外壳的约束仅影响边缘燃料肋片顶部区域的受力，中心燃料的力学行为主要与相邻燃料的约束有关。考虑到本文的目的是研究中心螺旋十字燃料棒在不同传热学边界条件下的热力响应，因此，仅以包含中心目标棒及周围的3×3棒束单元作为计算域进行热力耦合分析。此外，考虑到3×3棒束中4根角棒与中心棒之间无直接力学作用，分析中忽略4根角棒，仅建立图3所示的5根棒的计算域模型，单根燃料棒总长为0.8 m，其中燃料段长0.6 m，两端的渐变段和固定段各为0.05 m，通过周围4根棒为中间棒提供较为真实的约束条件。

图3 计算域模型 (a) 燃料组件的3D视图，(b) 扭转角90°、180°、270°截面图，(c) 单根燃料几何尺寸Fig.3 Computational domain of the rod bundle (a) 3D view of the rod bundle, (b) Cross section at the 90°, 180°, 270° plane,(c) Dimension of a single rod

1.3.1 网格划分

采用六面体网格对HCF组件进行网格划分。单根燃料在包壳外表面沿周向存在184个节点，其中肋片顶部区域的网格较密，可以提高接触计算的精度。从燃料中心到肋片顶部径向共有13个节点，在肋片根部燃料的热流密度较高，温度梯度较大，为准确计算肋片根部的温度场，从燃料中心到包壳外表面沿径向共划分18个节点。各处包壳中沿径向均有4个节点。HCF组件燃料段的轴向网格节点数为180，单根HCF棒束网格总数为317 088，共计374 000个节点，具体如图4所示。

图4 单根燃料的网格划分(a) 单根燃料网格，(b) 轴向截面网格Fig.4 Mesh partitioning of a single fuel(a) Mesh of the single rod, (b) Mesh of the axial section

1.3.2 对流边界条件

由于本文重点研究燃料的力学行为，为简化计算忽略燃料周围的流体区域，在燃料组件的包壳外表面上给定第三类边界条件，即对流传热系数和主流温度。第三类边界条件可表示为如下形式：

式中：n为壁面外法线；Tf为主流温度，K；Tw为壁面温度，K；h为壁面对流传热系数，W‧（m2·K）-1，模拟单相液体对流或两相沸腾传热工况时，分别采用Dittus-Boelter公式［14］和Chen公式［15］计算换热系数，稳态工况中流动传热边界条件见表1，在两相工况的计算达到稳定后引入破口失水事故和反应性引入事故以分析燃料在更高的运行参数下抵御事故的能力，瞬态事故工况的事故序列见表2。

表1 稳态工况中的传热边界Table 1 Heat transfer boundary under steady conditions

表2 瞬态事故工况的事故序列Table 2 Time sequence of the events for transient accidents

1.3.3 结构约束条件

对HCF组件进行力学模拟的外部约束条件如下：

1）燃料组件的四周设置一个无应变的刚性外壳（图3中6处），外壳的内壁冷态时与边缘棒发生点接触，在接触面的法向边缘棒的位移等于零，边缘棒与外壳和中心棒之间均存在摩擦，摩擦系数为0.1；

2）对燃料底部的固定段及圆形端面（图3中5处）施加径向、周向和轴向的固定约束，使该固定段的位移为零，模拟下管座的支撑和固定作用；

3）对燃料顶部的固定段（图3中1处）施加径向、周向的固定约束，在圆形端面设置轴向的弹簧约束，弹簧的刚度k为1.45×1010N‧m-1，使该固定段在径向和周向的位移等于零，轴向位移x受到弹簧反力F限制，具体方程为F+kx=0，模拟上管座的定位和压紧作用。

1.4 耦合及数值求解方法

本文使用ANSYS Workbench进行热力耦合分析，由于燃料棒变形幅度较小，不考虑燃料棒变形对其内部导热的影响。先通过热单元steady-state thermal（transient thermal）求解燃料的稳态（瞬态）温度场，再将数据单向传递至结构单元static structural进行稳态（瞬态）力学分析。采用Newton-Raphson迭代法求解燃料在接触区域和进入塑性阶段后的非线性方程，位移的容差因子为0.5%。

2 结果分析

本文对图3所示的计算域开展数值模拟研究，获得中心棒及周围四根棒的热力响应，并重点对中心棒的热力行为展开分析。螺旋十字燃料组件沿轴线发生扭转，扭转的角度与轴向高度成正比，为表明燃料组件的力学行为在轴向的周期性变化规律并对轴向高度进行无量纲化，可定义扭转角δ（rad）代替轴向高度z（0≤z≤0.6 m）表示燃料段的轴向位置，δ的表达式如下：

相邻燃料在扭转角为0°、90°、180°、270°和360°位置发生接触。

2.1 单相液体对流传热工况

反应堆正常运行时，绝大部分冷却剂处于单相液态，因此本文首先开展单相液体对流边界条件下的螺旋十字燃料组件热力行为研究，冷却剂入口温度为553.85 K，入口速度为1 m‧s-1，燃料棒的体积释热率为740 MW‧m-3，在冷却剂入口温度、速度和棒功率的基础上，通过热平衡计算冷却剂温度。利用Dittus-Bolter公式基于通道平均热工参数计算第三类边界条件所需的对流换热系数，换热系数沿棒周向为常数。通过计算获得中心棒的von Mises应力分布如图5所示，可以看出相邻燃料发生接触的自支撑平面上肋片顶部区域的von Mises应力较高，相对其他区域更可能发生塑性变形甚至破裂，因此本文主要对接触平面上燃料的受力和变形进行分析，其中扭转角180°平面上的中心棒温度、von Mises应力和塑性应变分布（同高度截面冷却剂平均温度为559.4 K）如图6所示。燃料温度最大值677 K处于芯块中心处，周围的等温线形状为圆形。而靠近包壳的燃料温度由于受冷却剂对流传热影响更大，等温线形状与燃料的十字形轮廓相似。等温线形状的变化使得燃料内部径向的温度梯度沿周向出现较大差异，因温度梯度产生的热应力在周向上也呈现不均匀分布（如图6（b）所示）。中心棒与边缘棒接触的区域温度梯度较小，但由于接触约束的存在，受应力集中作用，该处von Mises应力出现极大值，超过了屈服极限，发生明显的塑性变形（图6（c））。从图6（d）和（e）可以看出，燃料芯块和包壳中弹性应变和热膨胀应变的最大值仅为塑性应变最大值的十分之一，因此本文仅对塑性应变进行分析。

图5 中心棒von Mises应力分布云图Fig.5 Distribution of von Mises stress on the center rod

图6 180°平面中心棒的温度、von Mises应力和应变分布(a) 温度分布，(b) von Mises应力分布，(c) 塑性应变分布，(d) 弹性应变分布，(e) 热膨胀应变分布Fig.6 Temperature, von Misesstress and strain at the 180° plane (a) Temperature distribution, (b) von Mises stress distribution,(c) Plastic strain distribution, (d) Elastic strain distribution, (e) Thermal strain distribution

为进一步量化相邻棒约束作用和燃料内部温度梯度对中心棒应力分量的影响，图7给出了中心棒扭转180°平面包壳外表面处三个应力分量的分布曲线（平面内0°位置见图3（b））。图中0°、90°、180°和270°的位置区域（以下称翼尖处）受接触约束影响更大，由于接触约束主要限制燃料的径向热膨胀，对周向和轴向的热膨胀影响较小，因此翼尖处的径向应力为-700 MPa左右的压应力，周向应力和轴向应力仅分别为径向应力的55%和70%。而45°、135°、225°和315°的位置区域（以下称翼根处）离约束位置较远，受接触约束影响较小，从图6（a）可以看出，翼根处的温度梯度较大，说明该处的热应力与温度梯度相关，温度梯度使得包壳外侧的温度较低，热膨胀量比内侧小，限制了内侧的热膨胀，导致包壳外侧在周向和轴向受到反作用的拉应力，数值分别可以达到250 MPa和130 MPa，由于文章中忽略了流体区域，不考虑流体压力对燃料包壳的影响，径向应力接近于零。

图7 180°平面包壳外表面处的应力分量分布（加粗刻度线为零应力线）Fig.7 Distributions of stress components on the outer surface of cladding at the 180° plane(the bold scale line represents the zero stress scale)

螺旋十字燃料棒在扭转角为90°、180°和270°平面发生约束定位，不同扭转角平面冷却剂的平均温度不同，图8展示了不同扭转角平面包壳外表面处von Mises应力和塑性应变的分布情况。可以看出，三个扭转角平面包壳外表面处的von Mises应力和塑性应变曲线均以90°为周期变化，von Mises应力在翼尖和翼根处较大，分别为300 MPa和216 MPa左右，塑性应变在翼尖处达到最大值0.04，说明发生约束定位的平面包壳外表面处的von Mises应力和塑性应变的变化趋势基本相同。但从局部放大图中可以看出，在扭转角为90°和270°平面上，von Mises应力和塑性应变曲线分别向逆时针方向和顺时针方向偏转，这是燃料轴向沿螺旋线向两个固定端膨胀的结果。由于燃料顶部固定端的弹簧刚度较大，燃料所受约束在轴向几乎关于扭转角为180°平面对称，因此这一平面上偏转现象不明显。此外区别于点接触时的应力应变尖峰，单相工况下翼尖处出现了一段跨度约2°，峰值为300 MPa的von Mises应力平台和数值为0.02~0.04的塑性应变平台，说明翼尖处已经从冷态时的点接触转变为面接触。

图8 不同扭转角平面包壳外表面处的von Mises应力和塑性应变分布(a) von Mises应力（加粗刻度线为零应力线），(b) 塑性应变Fig.8 Distributions of von Mises stress and plastic strain on the cladding outer surface with different helical angles(a) von Misesstress (the bold scale line means the zero stress scale), (b) Plastic strain

2.2 沸腾传热工况

针对堆芯过热通道中冷却剂可能出现的流动沸腾，本文研究了冷却剂处于15.5 MPa下的饱和温度617.9 K时体积释热率为740 MW‧m-3的棒束的热力响应，由于计算区域加热长度较短，传热量较少，并且本文重点关注棒中间远离两端影响的高度，即扭转角为180°的位置，因此假设冷却剂的质量含气率保持不变，为0.1。图9给出了包壳内外表面的温度分布，可以看出由于两相工况中冷却剂的温度较高，翼尖处和翼根处的温度分别比单相工况中高54 K和38 K。但棒束的线功率密度保持不变，单相和两相工况中翼尖处和翼根处包壳内外表面的温差相同，分别为5 K和20 K左右。

图9 不同工况中180°平面包壳内外表面的温度分布Fig.9 Temperature distribution on the cladding inner and outer surfaces at 180° plane under different operating conditions

为研究温度的差异对包壳力学行为的影响，图10比较了单相和两相工况扭转180°平面包壳外表面处的von Mises应力和塑性应变分布。由于翼根处的热应力与包壳内外表面的温差有关，而两种工况中温差基本相同，翼根处的von Mises应力在两种工况中同为215 MPa。从图10（b）可以看出，两相工况翼尖处的塑性应变较单相工况偏大15%~25%，说明翼尖处因热膨胀产生的塑性应变随包壳温度的上升而增加。但图10（a）中两相工况翼尖处的von Mises应力却较单相工况偏低6%~10%，这是由于温度影响了包壳的力学性能。图11展示了不同温度下包壳在塑性阶段的应力强化曲线，可以看出温度越高，包壳的屈服极限越低，说明相较于单相工况，两相工况在较低的应力水平下进入塑性阶段，塑性阶段应力随应变的增加速度较弹性阶段放缓；图11中曲线的斜率还随着温度的上升而减小，即同处于塑性阶段时两相工况中包壳的应力随塑性应变的增加速度比单相工况中慢。这两点原因共同导致两相工况中翼尖处von Mises应力比单相工况小。

图10 不同工况中180°平面包壳外表面的von Mises应力和塑性应变分布(a) von Mises应力（加粗刻度线为零应力线），(b) 塑性应变Fig.10 Distributions of von Mises stress and plastic strain on the cladding inner and outer surfaces at 180° plane under different operating conditions (a) von Mises stress (the bold scale line means the zero stress scale), (b) Plastic strain

图11 包壳材料在不同温度下应力强化曲线Fig.11 Stress intensification curves for cladding material under different temperature

2.3 反应性引入事故

本文研究的反应性引入事故为正常运行时控制棒突然提出堆芯，导致堆芯出现较大的正反应性的事故，即弹棒事故。由图10可知，单相和两相工况中扭转角180°平面包壳的应力应变最大值均出现在翼尖处和翼根处，出于安全分析的保守性考虑，本文对事故过程中包壳翼根处和翼尖处的力学行为进行研究并于包壳的内外表面设置了4个检测点，如图12（a）所示。反应性引入事故于1 s时在处于额定功率正常运行的反应堆中发生，1.1 s时控制棒完全弹出堆芯，导致燃料功率在1.24 s达到峰值，为额定功率的1.83倍，包壳的温度也随之上升，监测点中4处的温度最高，最大值出现在1.6 s，为662 K，可以将事故以该时刻为界分为前期和后期两个阶段。由于反应堆存在固有安全性，1.24 s后燃料功率不会继续增加。与此同时收到1.1 s发出的停堆信号后，控制棒在2 s时开始下落并在5 s时完全插入堆芯，燃料功率迅速下降，5.5 s时仅为额定功率的20%，包壳温度也降低至630 K以下，图12（b）为反应性引入事故中燃料功率［16］与各监测点温度的变化情况。

图12 反应性引入事故中燃料的功率和各监测点温度的变化 (a) 监测点位置，(b) 功率和温度的变化Fig.12 Variations of fuel power and temperature at different monitoring points(a) Locations of the monitoring points, (b) Power and temperature variations

反应性引入事故中温度的起伏变化会极大影响包壳的受力和变形，图13展示了翼尖处包壳内外表面的应力随时间的变化情况。可以发现，由于接触约束限制了包壳径向的热膨胀，翼尖处径向应力的变化幅度明显高于周向应力和轴向应力，其中包壳内表面的径向应力在1 s和7.5 s时的差值分别是周向应力和轴向应力的数百倍和4.6倍，包壳外表面则是3倍和1.7倍，说明翼尖处包壳的受力变化主要由径向应力引起。从von Mises应力的变化趋势中可以看出，2 s前翼尖处包壳内外表面的von Mises应力均基本保持不变，2 s后包壳内表面的von Mises应力始终随着燃料温度的降低而降低，7.5 s时不足100 MPa；而包壳外表面的von Mises应力则是先下降再上升，最小值52 MPa出现在2.36 s，7.5 s时von Mises应力上升至超过该温度下的屈服极限281 MPa。出现这一差异主要是由于包壳经历的力学过程不同：在反应性引入事故期间，燃料温度先上升后下降，包壳从受热膨胀变为冷却收缩，发生先加载后卸载的过程，但由于塑性变形不会在卸载过程中减小，燃料温度下降后塑性应变无法恢复到事故前的水平，燃料冷却收缩过程中翼尖处包壳外表面会出现反向加载的现象，径向应力在2 s后先逐渐降低至零然后反向增大，如图13（b）所示。而包壳内表面则始终处于卸载状态，径向应力在2 s后逐渐降低，如图13（a）所示。因此当燃料温度降低时包壳内外表面的von Mises应力呈现不同的变化趋势。

图13 反应性引入事故中翼尖处包壳内外表面应力的变化 (a) 包壳内表面，(b) 包壳外表面Fig.13 Stress at the cladding inner and outer surfaces of blade tip under reactivity insertion accident(a) Cladding inner surface, (b) Cladding outer surface

图14进一步给出了翼尖处包壳内外表面塑性应变随时间的变化规律。事故前期翼尖处包壳内表面的塑性应变随着温度的增加上升了15%左右；包壳外表面的塑性应变由于接触约束限制基本保持不变。事故后期包壳内表面的塑性应变保持不变，为0.036 9，说明包壳内表面仅发生卸载过程；而3 s后包壳外表面的塑性应变则从0.033下降到0.027，说明3 s后包壳外表面处发生反向加载过程并产生反向的塑性变形，抵消了一部分事故前产生的塑性应变，使得7.5 s时的塑性应变相对于事故前降低。

图14 反应性引入事故中翼尖处包壳内外表面塑性应变的变化 (a) 包壳内表面，(b) 包壳外表面Fig.14 Plastic strain at the cladding inner and outer surfaces of blade tip under reactivity insertion accident(a) Cladding inner surface, (b) Cladding outer surface

与接触约束影响翼尖处的径向应力不同，翼根处的受力与包壳内外表面的温差有关，图15展示了反应性引入事故中翼根处包壳内外表面应力的变化情况。由于翼根处包壳内外表面的温差从1 s时的21 K降低到7.5 s时的3.7 K，受温度梯度影响的周向应力和轴向应力也出现大幅下降，数值分别可以达到100~170 MPa和60~80 MPa，而径向应力的变化幅度则小于10 MPa，即径向应力的变化幅度远小于周向应力和轴向应力的变化幅度，von Mises应力的变化主要由周向应力和轴向应力的变化引起。从图15还可以看出，翼根处包壳内外表面的周向应力与轴向应力2 s前上升幅度不超过20 MPa，2 s后均呈现下降的趋势，没有出现方向变化，并且从图16可以看出，事故后期包壳翼根处的塑性应变保持不变，说明翼根处在事故后期仅发生卸载过程，没有出现反向加载。虽然翼根处的塑性应变随着温度的增加在1.6 s时达到事故前的数倍，但由于翼根处的塑性应变仅为翼尖处的1/10，因此从防止包壳破裂的角度考虑更需要关注翼尖处塑性应变的变化。

图15 反应性引入事故中翼根处包壳内外表面应力的变化 (a) 包壳内表面，(b) 包壳外表面Fig.15 Stress at the cladding inner and outer surfaces of blade elbow under reactivity insertion accident(a) Cladding inner surface, (b) Cladding outer surface

图16 反应性引入事故中翼根处包壳内外表面塑性应变的变化 (a) 包壳内表面，(b) 包壳外表面Fig.16 Temperature at the cladding inner and outer surfaces of blade elbow under reactivity insertion accident(a) Cladding inner surface, (b) Cladding outer surface

2.4 一回路压力边界破口失水事故

在处于满功率正常运行的反应堆中1 s时引入一回路压力边界破口失水事故，事故发生后冷却剂从破口中大量流出，堆芯压力迅速下降，1.6 s时降至12.6 MPa并触发低压停堆信号，忽略停堆响应时间和控制棒下落时间，不考虑冷却剂空泡效应的影响，反应堆功率在1.6 s前保持不变并在1.7 s降低至事故前的十分之一，破口失水事故中堆芯压力［17］和燃料功率［18］的变化曲线如图17所示。事故中假设一回路冷却剂始终处于饱和沸腾阶段，饱和温度随堆芯压力降低而降低。

图17 破口失水事故中堆芯压力和功率的变化Fig.17 Variations of core pressure and power under the LOCA condition

为分析一回路压力边界破口失水事故中包壳翼尖处和翼根处的温度、变形和受力的变化情况，仍选取§3.3中的4个监测点作为研究对象（图12（a）），图18和图19分别给出了破口失水事故中翼尖处各监测点应力和塑性应变的变化曲线。可以看出，从事故发生到停堆包壳内外表面的温差基本不变，但由于冷却剂的温度随堆芯压力下降而降低且燃料棒功率保持不变，包壳内外表面的温度下降了10 K，结合图19中停堆前包壳内外表面的塑性应变保持不变，说明停堆前翼尖处发生卸载过程，包壳内外表面的von Mises应力分别降低37 MPa和20 MPa。停堆后0.2 s内，包壳内外表面温度下降的平均值同为10 K，但温差则由于堆芯功率的剧烈下降而降低至停堆前的40%，包壳内外表面von Mises应力的变化幅度则可以达到停堆前的3.5倍和10倍，说明翼尖处包壳的受力不仅与限制热膨胀的接触约束有关，同时也受包壳内外表面温差的影响。从图18中还可以看出，停堆后短时间内包壳外表面的von Mises应力就已超过屈服极限，而包壳内表面的von Mises应力在9.5 s时仍小于200 MPa，图19中包壳外表面的塑性应变由0.034 2下降至0.018 4而包壳内表面的塑性应变保持0.034不变，说明包壳外表面先发生卸载后出现反向加载而包壳内表面仅存在卸载过程。

图18 破口失水事故中翼尖处包壳内外表面应力的变化 (a) 包壳内表面，(b) 包壳外表面Fig.18 Stress at the cladding inner and outer surfaces of blade tip under the LOCA condition(a) Cladding inner surface, (b) Cladding outer surface

图19 破口失水事故中翼尖处包壳内外表面塑性应变的变化 (a) 包壳内表面，(b) 包壳外表面Fig.19 Plastic strain at the cladding inner and outer surfaces of blade tip under the LOCA condition(a) Cladding inner surface, (b) Cladding outer surface

为研究一回路压力边界破口失水事故中包壳翼根处应力和塑性应变的变化情况，本文绘制了翼根处各监测点应力和塑性应变随温度和包壳内外表面温差的变化曲线如图20和图21所示。从图20可以看出，停堆前翼根处包壳内外表面的温度下降了6 K而温差基本不变，von Mises应力也基本不变。停堆后0.05 s内，翼根处的温度同样下降了6 K左右，包壳内外表面的温差降低到停堆前的80%，von Mises应力出现下降，说明翼根处离接触位置较远，包壳整体温度的下降不会使该位置发生卸载，应力受包壳内外表面温差影响更大。停堆前由于温差存在0.2 K的小幅度上升，图21中可以观察到包壳内外表面的塑性应变增加了1%~3%，但仍远小于翼尖处。停堆后包壳的von Mises应力不断下降且未超过屈服极限，翼根处仅发生卸载过程，塑性应变保持不变。

图20 破口失水事故中翼根处包壳内外表面应力的变化 (a) 包壳内表面，(b) 包壳外表面Fig.20 Stress at the cladding inner and outer surfaces of blade elbow under the LOCA condition(a) Cladding inner surface, (b) Cladding outer surface

图21 破口失水事故中翼根处包壳内外表面塑性应变的变化 (a) 包壳内表面，(b) 包壳外表面Fig.21 Plastic strain at the cladding inner and outer surfaces of blade elbow under the LOCA condition(a) Cladding inner surface, (b) Cladding outer surface

3 结语

本文对典型螺旋十字燃料HCF棒束单元进行热力耦合特性研究，模拟了包括单相对流、饱和沸腾、反应性引入事故和失水事故在内的多种工况，重点分析中心棒在自定位平面上的传热状况和力学行为，评估了HCF组件的安全性能。根据计算结果，可以得到如下结论：

1） 所有工况下HCF组件中心棒包壳外表面处的von Mises应力和塑性应变最大值总是出现在相邻燃料的接触区域，其次是燃料翼片的翼根处。不同轴向位置发生约束定位的平面上包壳外表面在周向的受力与变形的变化趋势相同，但相位存在偏差。翼尖处的应力主要受接触约束影响，同时也和包壳内外表面温差有关；翼根处的应力则与包壳内外表面的温差有关。

2） 反应堆正常运行冷却剂处于单相液态时燃料包壳翼尖处会出现0.02~0.04的塑性应变平台和峰值为300 MPa的von Mises应力平台，说明翼尖处已经从点接触变为面接触。

3） 相较于单相液体冷却工况，饱和沸腾工况下包壳的平均温度高35~60 K，但包壳内外表面的温差相同，导致翼根处von Mises应力相等，翼尖处塑性应变偏大15%~25%，但由于温度越高包壳的屈服极限越低、塑性阶段应力随温度的升高对塑性应变的增加越不敏感，翼尖处von Mises应力反而偏小6%~10%。

4） 反应性引入和失水事故工况下，包壳的von Mises应力和塑性应变分别小于350 MPa和0.04，且包壳外表面温度低于锆水反应温度；停堆后，翼尖处包壳外表面发生反向加载并再次进入塑性阶段，导致该处塑性应变降低，von Mises应力再次上升并超过屈服极限，需要注意这一反向应力对燃料的二次损伤。

本文对HCF组件热力耦合特性进行了研究，得到了燃料组件的力学响应，但由于HCF组件的热流密度沿周向并非均匀分布，冷却剂在同一轴向平面上的温度也不尽相同，因此在热力耦合的基础上还需要对冷却剂的流动传热进行模拟。此外，本文的热力分析中仅考虑了寿期初的燃料组件，随着反应堆的运行，燃料燃耗深度的增加会改变芯块和包壳材料的物性，辐照对物性的改变以及辐照肿胀现象本文也未分析，可以增加冷却剂计算模块、建立燃料的辐照和燃耗模型使HCF组件的力学计算更加准确。

作者贡献声明丛腾龙负责立题和设计计算工况，数据分析，论文修改；刘峪洁负责数值模拟，数据分析，论文初稿；郭辉数值模拟；肖瑶技术支持；顾汉洋获取研究经费。