1. 引言
控制棒组件(Rod Cluster Control Assembly, RCCA)是一种用于控制反应堆内反应性的装置，通过将中子吸收材料插入或抽出反应堆堆芯来实现其功能。在核反应堆正常运行时，RCCA主要用于控制反应堆的反应性和功率分布，以确保反应堆安全稳定运行。在核反应堆发生事故时，RCCA通过迅速插入到堆芯中使反应堆紧急停堆，以防止核燃料的过热以及事故的进一步发展，保证核反应堆的安全性。
控制棒因自身具有细长的结构，当不稳定湍流状态的冷却剂施加作用力在控制棒组件上时，控制棒会受到冷却剂流动冲击的影响，与导向机构发生碰撞/切向滑动，产生微振磨损，进而影响控制棒可靠性。
目前，国内外对于控制棒包壳磨损的研究方法包括试验和数值模拟。试验研究有两种主要形式，一种是基于单根控制棒及其导向结构设计的试验台架研究控制棒流致振动及磨损问题。B. Reynier等[1]在试验条件接近压水堆主系统环境下，针对无水流流动情况，通过高压磨损釜试验台架对控制棒不锈钢包壳(316L)进行了两次局部冲击滑动磨损试验，研究了激励频率和运动信号类型对磨损过程和损伤的影响。A. Van Herpen等[2]针对控制棒组件因流致振动导致的磨损问题，采用CANDUSE磨损试验台架在室温下对控制棒304L不锈钢包壳进行了冲击滑动磨损试验，试验表明测试时间对材料的行为没有影响，室温下的主要磨损机制是表层的氧化，其次是因为冲击滑动运动导致的氧化物脱落。S. N. Kim[3]等为了确定控制棒组件振动的主要机制，对典型控制棒组件的微动磨损现象进行了试验，并分析了韩国Ulchin反应堆1、2号机组涡流探伤结果及其规律。另一种是基于简化的局部控制棒导向筒模型，研究控制棒包壳局部的磨损的机理和规律。A. Laure Kaiser等[4]在类似一回路的环境(温度、压力和水化学)下，采用摩擦磨损测试装置开展控制棒包壳相对于导向装置的滑动冲击运动试验研究。研究表明，即使在低接触压力条件下，机械效应在控制棒包壳的磨损中起着主导作用。E. Marc等[5][6]针对控制棒组件出现的磨损问题，采用了一种能独立控制切向滑动和法向力波动的微动磨损系统，设计了在不同环境下进行往复微动的试验装置，对316L不锈钢圆管和304不锈钢平板开展了微动磨损试验，研究了环境温度下不同环境(空气和溶液)中的微动滑动磨损行为，并提出了一种新的能量磨损损伤定律。
对于数值模拟研究，针对典型微动磨损现象，部分学者利用有限元程序建立了微动磨损的有限元分析模型和方法。霍永忠等[7]基于Archard模型以及多层节点更新方法，建立二维平面/柱面微动接触模型，通过有限元方法模拟了锆合金的微动磨损行为，结果表明部分滑移与整体滑移状态下的微动磨损行为存在差别。李玲等[8]基于修正的能量磨损模型，建立了二维平面/柱面的动态微动磨损模型，分析了不同加载条件下柱面/平面微动磨损规律。G. Kermouche等[9]通过有限元模拟和试验观察的方法，根据数值计算结果和试验观察结果，提出了一种磨损机理，认为带有冲击过程的微动滑动磨损行为可以用带有刚性球形粗糙体的双层材料的刮擦过程来模拟。对于控制棒组件的磨损现象和行为，鉴于控制棒在堆内环境和运行工况的复杂性，目前的数值模拟分析聚焦于控制棒的流致振动特性研究。张惠民等[10]使用ANSYS软件对上腔室控制棒导向筒的内外流场进行数值模拟，获得了流场的流速分布，探究了悬挂状态下控制棒在流场中的流致振动特性，通过自编程，着重模拟了进出口即横向流和纵向流流量改变条件下的控制棒振动特性。Zhang等[11]基于CFD方法使用Code_Saturne以及ANSYS FLUENT对控制棒导向筒内的流场和压力进行了详细的数值模拟及对比验证。C. Phalippou等[12]对典型AECL台架的冲击力和磨损率进行了试验和数值研究，测量了控制棒的模态振型和振动频率，并与有限元计算结果进行了比较，计算了磨损试样之间的碰撞和滑动运动，结果表明，冲击力、平均磨损速率以及试样的相对运动与试验数据较为接近。
由于控制棒组件结构的复杂性以及堆芯上腔室流场分布的影响，试验研究仍是研究控制棒组件磨损的主要手段。随着计算机技术的发展，基于数值模拟方法研究磨损行为日趋成熟。本研究将能量磨损模型和有限元分析方法成功应用到控制棒包壳的磨损研究中，并针对控制棒包壳和导向机构在不同工况条件下的微动磨损现象和行为开展研究，揭示不同参数条件对磨损深度和形貌计算的影响规律，为后续相关研究提供参考。
2. 模型和建模
2.1. 磨损模型
在微动磨损仿真中，选择合适的磨损模型至关重要。目前常用的经典磨损模型有Archard模型和能量磨损模型。
2.1.1. Archard模型
Archard模型[13]是一种很常用的磨损模型。在有限元分析金属滑动磨损性能时，Archard模型被普遍应用与仿真分析中也验证了其对磨损预测有较好的普适性。因此可以利用此模型分析控制棒包壳微动磨损性能。
Archard理论模型指出磨损量与滑动距离、接触面间接触压力成正比，同材料的硬度成反比，即：
(1)
式中，Vw为材料表面体积磨损量，mm3；s代表相对滑移距离，mm；H为材料的硬度；k是磨损系数；FN为接触压力，N。
2.1.2. 能量磨损模型
能量磨损模型[14]基于能量耗散概念，其假设磨损是由接触表面上的能量耗散导致的，其磨损体积是接触面上的摩擦力和滑移距离的乘积：
(2)
式中，V是磨损体积，mm3；Q是接触剪切力，N；S是滑移距离，mm；K是磨损系数。
总的磨损体积V是所有微元(Vi,i表示磨损内表面的位置)磨损体积之和，微元的磨损体积为微元磨损面积与磨损深度之积，因此，针对第i个微元，其磨损体积为：
(3)
式中，Ai为i位置处的微元磨损面积，mm2；Di为i位置处的磨损深度，mm，Si是i位置处的滑移距离，mm，Qi是i位置处的接触剪切力，N。从方程可以得到，i位置处的磨损深度为：
(4)
式中，
为i位置处的剪应力。
Archard磨损模型，以其直观的形式和简便的使用方法，在工程实践中被广泛采纳。然而，该模型也存在一些实质性的不足。磨损系数k，作为一个比例常数，对于给定的材料其值并非恒定不变，而是受载荷、滑动速度、温度、磨损类型及其它摩擦学条件的影响而变化。这意味着，即便对于相同的材料组合，k值也可能因测试条件、设计参数和接触几何的不同而有所差异。
能量磨损模型基于能量耗散原理，提供了一种更为深入的磨损机理解释。与Archard模型相比，能量磨损模型的磨损系数K不显著依赖于材料特性、位移幅值、接触几何等参数，从而降低了对特定测试条件的依赖性。尽管能量磨损模型在实施上需要大量的剪切摩擦力与滑移量数据，并且需要较高的采样密度以积分获得能耗值，增加了实现的复杂度，但其磨损系数不受摩擦系数变化影响的特点，具有更高的准确性。
因此，基于上述考虑，本研究选择采用能量磨损模型对控制棒包壳管的磨损问题进行计算分析。
2.2. 有限元建模
2.2.1. 磨损分析方法
随着磨损的进行，磨损界面的几何形状会随之发生变化。由于磨损导致的网格畸变可能会显著降低计算结果的精确度，甚至导致计算过程无法收敛。因此，有必要对磨损后的网格进行实时更新。本研究采用了ALE自适应网格技术，在分析对象出现材料的损失的情况下，对磨损发生位置的网格进行自适应调整，通过对网格进行约束并定义网格平滑算法，重新划分高质量的网格。
本文的磨损有限元分析流程如图1所示。采用有限元软件ABAQUS，利用Fortran语言编写用户子程序UMESHMOTION，开展控制棒包壳/导向机构的界面磨损情况模拟，具体步骤包括：确定磨损区域节点及磨损方向；开展有限元分析，求解获得磨损区域节点的接触应力以及相对滑移距离；在每一次增量步结束之后，调用子程序UMESHMOTION，根据磨损模型，计算磨损深度和宽度；根据ALE自适应网格设定的更新频率，更新节点i处坐标
，进而生成新的网格轮廓；如此反复上述流程。
在进行控制棒包壳管的微动磨损分析时，考虑到微动磨损循环次数通常极为庞大，可达数十万至数百万次，每次循环引起的磨损量极为微小，表面形貌的变化可能并不显著。鉴于此，在实施磨损的有限元分析过程中，广泛采纳了循环加速磨损模型的策略。该策略在预设的循环次数内，忽略由磨损引起的物体表面轮廓的渐进变化，而是在累积至一定循环次数之后，再对物体的接触表面轮廓进行更新。为了实现这一加速计算，本研究在磨损模型框架内引入了ΔN次循环载荷的概念，其中每个载荷步长代表ΔN次循环载荷的累积效应，以此提高计算效率，同时确保磨损模拟的准确性。
Figure 1.Flowchart for finite element analysis of wear
图1.磨损的有限元分析流程
2.2.2. 有限元建模
控制棒和导向机构的如图2所示，本研究的控制棒包壳切向振动幅度在100 μm左右，其和导向机构内弧面之间微动距离较小时，可采用包壳管和平板接触来近似包壳管和圆弧面的接触情况。本研究基于E. Marc等[5][6]圆管和平板磨损试验工况开展模拟，试验中采用AISI 316L不锈钢管与AISI 304L不锈钢平板，圆管样品来自真实的反应堆控制棒包壳管，外径为9.7 mm。试验如图3所示，圆管和平板之间相互接触，对圆管施加法向载荷P，同时圆管和平板之间发生位移幅值为δ、频率为f的切向滑动。试验过程考虑了两种工况：室温下的干燥环境；硼锂溶液的湿环境。
Figure2.Control rod guide structure
图2.控制棒导向结构[15]
Figure3.Tube-plate wear test
图3.圆管–平板磨损试验
建立如图4所示的包壳管–平板磨损模型，包壳外径9.7 mm、壁厚0.45 mm，包壳和平板的材料均为不锈钢管，分析过程中弹性模量和泊松比均取200 GPa和0.3。运动过程中，包壳管和平板始终接触，其中，下部平板保持固定，上部包壳管在顶部均匀载荷和切向力的作用下，进行切向往复运动。两者发生相对滑动，产生微动磨损。
(a) 三维模型 (b) 二维模型
Figure4.Finite element analysis model
图4.有限元分析模型
对于包壳管和平板接触，采用双部件表面–表面接触设置，对两个接触对进行双向接触设置。采用有限滑移算法，法向设置为硬接触，切向约束通过用户子程序FRIC_COEF定义时变摩擦系数，切向摩擦系数根据E. Marc试验提供的结果进行拟合，拟合结果如图5所示。
(a) 干环境下的摩擦系数 (b) 湿环境下的摩擦系数
Figure5.Friction coefficient obtained from test
图5.试验摩擦系数
分析中能量磨损模型的磨损系数基于E. Marc试验结果来考虑。对于干燥环境和湿环境，能耗法磨损系数为：
干燥环境下平板：
；
干燥环境下圆管：
；
湿环境下平板：
；
湿环境下圆管：
。
3. 有限元分析
3.1. 二维分析结果
磨损工况共有干燥环境、湿环境两种工况，切向微动位移幅值干环境为40、160 μm，湿环境为40、100 μm。磨损范围均大于切向位移幅值但受接触状态影响，在相同磨损次数下湿环境工况磨损范围小于干燥环境工况。最大磨损量出现在部件横向中心位置处，随微动位移幅值增大而增大。二维模型的磨损分析结果见图6。
(a) 湿环境下包壳管磨损 (b) 湿环境下平板磨损
(c) 干环境下包壳管磨损 (d) 干环境下平板磨损
Figure6.Wear analysis results under two-dimensional model
图6.二维模型磨损分析结果
3.2. 三维分析结果
三维有限元模型平板的磨损分析结果如图7所示，三维平板磨损在空间上呈对称分布，最大磨损深度出现在模型磨损区域的几何中心，在轴向截面上磨损形貌呈现U型。
(a) 平板磨损三维模型结果
(b) 横向磨损深度 (c) 轴向磨损深度
Figure7.Wear analysis results of plate under three-dimensional model in wet environment
图7.湿环境下三维模型平板磨损分析结果
三维有限元模型包壳管的磨损分析结果如图8所示，同样的，三维包壳管磨损在空间上呈对称分布，最大磨损深度出现在模型磨损区域的几何中心。受法向载荷与微动幅值影响，包壳管存在部分滑移、混合滑移与完全滑移状态。
(a) 圆管磨损三维模型结果
(b) 横向磨损深度 (c) 轴向磨损深度
Figure8.Wear analysis results of cladding under three-dimensional model in wet environment
图8.湿环境下三维模型包壳管磨损分析结果
3.3. 分析结果讨论
a) 由于粘着作用，当接触中心由位移产生的剪切应力无法克服剪切摩擦力时，磨损处于部分滑移状态。如图9所示，部分滑移状态下最大磨损深度出现在接触中心点的左右两侧且呈“M”型，完全滑移状态最大磨损深度增大且接触区域中心区域呈倒“U”型。完全滑移状态的磨损量远大于处于部分滑移状态的磨损量。
(a) 包壳管磨损 (b) 平板磨损
Figure9.Comparison of partially slip condition and gross slip condition
图9.部分滑移状态和完全滑移状态对比
b) 湿工况状态下的介质润滑作用导致了较小摩擦系数与磨损系数，在相同载荷及位移幅值条件下，湿工况状态下的磨损远小于干工况，两种环境的结果图10所示。
c) 相同工况下三维模型和二维模型的最大磨损量对比如图11所示，计算结果表明，对于包壳管，三维模型和二维模型最大磨损深度相差约为3.4%；对于平板，最大差别约为12.8%。三维模型计算结果更精细，但大幅增加了计算时间成本与不收敛风险。因此，在不同参数影响规律研究中采用二维模型开展分析。
(a) 包壳管磨损 (b) 平板磨损
Figure 10.Comparison of dry environment and wet environment
图10.干燥环境和湿润环境对比
(a) 包壳管磨损 (b) 平板磨损
Figure 11.Comparison of wear analysis results between 3D and 2D models
图11.三维模型和二维模型磨损分析结果对比
4. 参数影响规律
4.1. 加速因子影响
对干燥环境下位移幅值40 μm工况加速因子分别选取100次、2000次、4000次进行磨损计算，其结果如图12所示。加速因子100次与2000次时，平板与包壳管累计磨损量计算结果相近，4000次加速因子计算磨损量计算有较为明显的偏差。对于本研究分析工况，100次与2000次工况的最大磨损深度和最大磨损范围误差在3%以内。因此，磨损计算中选取2000次加速因子以提高计算速度并保证计算精度。
(a) 包壳管磨损 (b) 平板磨损
Figure 12.Effects of acceleration factors
图12.加速因子影响
4.2. 摩擦系数影响
图5所示试验测量的摩擦系数表明，随着磨损的进行，摩擦系数趋向稳定。采用稳定后的摩擦系数开展分析，并与随时间变化的摩擦系数分析结果进行对比，图13所示结果表明，恒定摩擦系数和时变摩擦系数的分析结果差异较小。
(a) 包壳管磨损 (b) 平板磨损
Figure 13.Effect of coefficient of friction
图13.摩擦系数影响
4.3. 分析步影响
对于状态位移幅值40 μm工况，分析增量步时长分别设置0.001 s、0.0025 s，研究不同增量步时长的影响，计算结果见图14。当增量步时长为0.0025 s，可以有效提高计算效率，但模型收敛性降低导致最终计算结果产生较大误差。因此，磨损分析时增量步时长不宜过大。
(a) 包壳管磨损 (b) 平板磨损
Figure 14.Impact of the analysis step size
图14.分析步的影响
4.4. 长期磨损规律
长期磨损下的规律如图15所示，随微动循环次数的增大，平板与包壳管磨损区域深度和宽度均有明显的增加，并且平板的最大磨损深度随循环次数的增长远快于包壳管。随着磨损过程的推进，磨损的深度上持续增长，但是磨损宽度上初期变化较快，后期变化有所减缓。对于包壳管，随着磨蚀范围和深度的扩大，接触位置和接触状态变化导致局部可能出现部分滑移现象，磨损深度曲线出现了明显的转折。
5. 结论
对干燥环境下位移幅值40 μm工况加速因子分别选取100次、2000次、4000次进行磨损计算，加速因子100次与2000次时，平板与包壳管累计磨损量计算结果相近，4000次加速因子计算磨损量计算有较为明显的偏差。对于本研究分析工况，100次与2000次工况的最大磨损深度和最大磨损范围误差在3%以内。因此，磨损计算中选取2000次加速因子，以提高计算速度并保证计算精度。
(a) 包壳管磨损 (b) 平板磨损
(c) 磨损深度变化
Figure 15.Long-term wear evolution
图15.长期磨损规律
本研究通过有限元分析软件ABAQUS，对控制棒包壳的微动磨损进行了深入研究。通过建立能量磨损模型，并采用UMESHMOTION子程序对磨损表面进行动态更新，本研究成功模拟了控制棒包壳在不同工况下的磨损行为。研究结果表明：
1) 三维模型与二维模型的对比分析显示，三维模型能够更准确地考虑结构轴向的影响，但是计算成本和不收敛风险有所增加；
2) 润滑条件对控制棒包壳的磨损有显著影响。在湿环境下，由于介质的润滑作用，摩擦系数和磨损系数均较小，导致磨损量远小于干燥环境；
3) 通过对加速因子、摩擦系数以及分析步长的参数影响规律研究，确定了合理的计算参数，为后续的磨损计算提供了依据。加速因子的选择对计算精度和速度有重要影响，而摩擦系数的恒定与时变对磨损结果影响较小。分析步长的选择则需要在计算时间和模型收敛性之间取得平衡；
4) 长期磨损规律的研究揭示了随着微动循环次数的增加，平板与包壳管的磨损区域宽度显著增加，且平板的最大磨损深度增长速度明显快于包壳管。
综上所述，本研究为控制棒包壳的磨损机制提供了深入的理解，并为后续控制棒包壳的寿命预测和优化设计提供了有价值的参考。未来的工作将进一步优化计算模型，探索不同材料和工况下的磨损行为，以及开发更高效的磨损预测模型。