1. 引言
新的核燃料在制造完成后，燃料供应商会利用合适的运输手段将其运输到核电厂，入堆后进行使用[1]。另外，随着乏燃料池累积的乏燃料越来越多，后续可能会运输到别的燃料贮存地、燃料后处理厂或热室进行相关试验，此时就有乏燃料的运输需求[2]。在新燃料及乏燃料的运输过程中，往往会由于运输设备的启动、加速、制动以及运输道路的具体情况等，对燃料造成一定的振动与冲击，严重情况下会产生诸多不利影响[3]-[6]。为了防止燃料在运输途中可能造成的不良后果，不同国家或地区都制定了相关的核燃料运输法律规范[7]。但是，目前国内还没有一套完善的核燃料运输规范，而且当前我们并不明确运输过程中加速度可能造成的实际影响，尤其是燃料棒内芯块的状态。为了进一步了解运输载荷对燃料组件的影响，了解可能存在的易导致燃料组件损坏的工况，需要在运输加速度载荷下对燃料进行三维有限元分析，通过仿真分析的结果研究燃料在运输加速度载荷下的响应。
本文基于ABAQUS开展新燃料在运输载荷下的三维有限元分析，其结果可为燃料运输试验的设计和运输规范的制定提供一定参考。
2. 燃料棒参数及模型简化
2.1. 燃料棒主要参数及力学性能
表1给出了建模过程中考虑的燃料棒主要的几何参数[8]。
PNNL提出了一个模型，用于计算燃料棒包壳力学性能随各种因素的变化[9]。该模型描述了燃料棒包壳的应力–应变关系，弹性变形采用Hook定律，塑性变形采用幂指数形式，即：
Table 1. Main geometrical parameters of fuel rods
表1. 燃料棒主要几何参数
属性 |
参数值 |
芯块直径(mm) |
8.19 |
芯块高度(mm) |
13.5 |
碟形深度(mm) |
0.3 |
碟形球半径(mm) |
15.0 |
包壳厚度(mm) |
0.57 |
包壳外径(mm) |
9.5 |
(弹性变形阶段) (1)
(塑性变形阶段) (2)
其中，σ为应力(MPa)，ε为应变，E为杨氏模量(MPa)，K为强度系数(MPa)，
为应变速率(s−1)，n为应变硬化指数，m为应变速率指数。
将燃料棒芯块视为弹性材料，虽然在仿真中没有明确地模拟裂纹，但在估计其力学性能时考虑了裂纹的影响。将包壳视为弹塑性材料，其应变破坏极限在1%~4%之间，现选取1%作为保守值。燃料棒包壳和芯块的力学性能如表2所示，图1为参考表2生成的包壳应力–应变曲线，在转折点之前为弹性变形，转折点之后进入塑性变形阶段。
Table 2. Material parameters of fuel rods
表2. 燃料棒材料参数
包壳 |
Zr-4 |
新燃料芯块 |
UO2 |
密度(kg/m3) |
6590 |
密度(kg/m3) |
10970 |
弹性模量(GPa) |
77.64 |
弹性模量(GPa) |
201 |
强度系数(GPa) |
15.21 |
泊松比 |
0.32 |
应变强化指数 |
0.171 |
屈服强度(MPa) |
2146 |
应变速率指数 |
0.015 |
|
|
剪切模量(GPa) |
29.22 |
|
|
泊松比 |
0.328 |
|
|
屈服强度(MPa) |
764.92 |
|
|
极限拉伸强度(MPa) |
868.85 |
|
|
Figure 1. Stress-strain curve of fuel rod cladding
图1. 燃料棒包壳应力–应变曲线
2.2. 燃料棒模型的简化
模型简化的目的在于以较低的计算成本，准确预测在一定加速度载荷条件下燃料棒的应力分布。由于燃料棒属于截面尺寸远小于轴向长度的细长梁结构，对此类结构利用全尺寸三维网格直接计算会导致模型规模过大、计算成本过高，为此，首先建立了燃料棒的精细有限元模型，进行静力分析，得到燃料棒在静力载荷下的响应，并确定导致燃料棒失效的极限载荷。然后，在静力分析结果的基础上不断迭代、校准，得到一个等效梁模型。
在实际的燃料棒运输过程中，燃料棒承受着惯性力、支撑点施加的挤压载荷、与其他燃料棒碰撞产生的挤压载荷等复杂的载荷。对于细长型构件而言，在上述载荷条件下弯曲是使结构产生破坏的主要原因。因此，考虑施加弯矩对燃料棒段进行静力分析。在纯弯曲情况下，加载条件在燃料棒的整个长度上是一致的，由这种加载方式得出的结构性能与燃料棒段的长度无关。
以含有3个燃料芯块的燃料棒段作为模型简化的参考，采用六面体单元对包壳和芯块进行网格划分，如图2所示。其中，芯块之间、芯块与包壳之间通过定义无摩擦硬接触(Frictionless hard contact)来模拟其相互接触关系。模型中总共包含3120个单元和4404个节点。
Figure 2. Finite element model of fuel rod
图2. 燃料棒有限元模型
两端截面的节点分别利用MPC (多点耦合算法)参考点连接，并在这两个参考点上施加弯矩，通过逐渐提高弯矩的大小，确定可能使燃料棒产生破坏的极限载荷。
选择1%的塑性应变和2146 MPa的抗压强度分别作为包壳破坏和芯块的破坏判据。当施加的弯矩为35 N*m时，包壳的塑性应变为1%，芯块的Von-Mises应力仅为1600 MPa，此时包壳发生破坏。由此可以得到燃料棒段的临界弯矩为35 N*m，此时应力和应变云图如图3、图4所示。
为了使简化梁模型的力学响应与精细模型足够接近，应该仔细选择梁的材料属性和横截面参数。理想情况是，简化梁与精细模型的响应在达到临界载荷之前都完全一致的。然而，实际上很难找到这种绝对精确简化模型参数。根据前面的静力分析结果，燃料棒的破坏是由包壳达到1%塑性应变引起，由于塑性应变可以根据梁的曲率计算得到，而曲率可以根据沿燃料棒中心线的三个点的位移计算获得。选择弯矩–位移曲线上的三个点进行校准，如图5(a)所示。最左边的点选在弹性区域，最右边的点对应于产生1%塑性的临界力矩。
简化梁的直径与精细模型的包壳外径相同，设定为9.5 mm。对简化梁与精细模型施加同样大小的载荷和约束条件，以二者的变形误差为目标，经过不断的迭代校准，在极限载荷范围内寻求简化梁材料参数的最优解。当简化梁模型与精细模型在给定弯矩下具有相同的变形时，可以选择此时的参数作为简化梁的最终参数。最终得到精细模型与梁模型的弯矩–位移曲线对比如图5(b)所示。
Figure 3. Stress distribution of fuel rods under bending load
图3. 燃料棒弯曲载荷下的应力分布
Figure 4. Strain distribution of fuel rods under bending load
图4. 燃料棒弯曲载荷下的应变分布
Figure 5. Bending moment-displacement curve of fuel rod: (a) fine model and calibration point; (b) comparison between fine model and simplified beam model
图5. 燃料棒弯矩–位移曲线：(a)精细模型及校准点；(b)精细模型与简化梁模型对比弯曲载荷下的应变分布
3. 有限元模拟
3.1. 燃料棒动态响应分析
根据已得到的简化梁模型的等效材料参数，利用1D单元建立燃料棒的简化梁模型，同时采用Shell单元建立隔架有限元模型。根据实际运输过程燃料棒在空间隔架中的放置状态，燃料棒与隔架接触点(凸点和弹簧片)通过bush单元模拟隔架对燃料棒的支撑作用，bush单元的刚度由弹簧片和凸点决定，燃料组件有限元模型如图6所示。
Figure 6. Finite element model of fuel assembly
图6. 燃料组件有限元模型
动态响应分析所施加的载荷以真实的运输数据为参考，文献[10]记录了从内华达州的水星到新墨西哥州的阿尔布开克运输一个44,000磅和一个56,000磅核燃料运输容器的冲击和振动数据。由原始的道路数据转换成的加速度谱密度(ASD)如图7绿线和蓝线所示。由于原始数据的间隔过多且都是阶跃突变的载荷，难以在仿真模拟中直接使用，因此采用光滑包络载荷谱线来代替原始的载荷数据(图7红线)，作为本文分析中所施加的运输载荷，具体值如表3。考虑保守性，本文同时还分析了2倍运输载荷下的动态响应。
Figure 7. Original acceleration spectral density and envelope
图7. 原始加速度谱密度及包络线
Table 3. Main geometrical parameters of fuel rods
表3. 燃料棒主要几何参数
频率(Hz) |
ASD (G2/Hz) |
5 |
1.8e−3 |
20 |
1.8e−3 |
25 |
8.0e−4 |
125 |
8.0e−4 |
135 |
5.5e−4 |
265 |
5.5e−4 |
530 |
1.0e−4 |
1100 |
3.0e−6 |
2000 |
3.0e−6 |
3.2. 燃料芯块冲击分析
为了验证燃料芯块在运输过程中的完整性，选择最严酷的冲击条件进行数值模拟，即端部芯块固定，将其他芯块视为一个整体，当这一整体受到轴向脉冲加速度作用时，端部芯块受到的脉冲力可达到最大。
为了提高计算效率，将模型简化为局部细化的弹簧质量点模型。其中，燃料棒的包壳利用壳单元进行建模，端部的2.5个芯块进行精细化建模。此时可以将单个芯块的刚度等效为弹簧刚度，其他芯块的总质量等效为集中质量点。最终简化后的有限元模型如图8所示。
所施加的轴向脉冲加速度大小为4 g，持续时间分别为1 ms、3 ms、5 ms。
Figure 8. Simplified model of pellets impact analysis
图8. 芯块冲击分析简化模型
4. 结果分析
4.1. 燃料棒
对整个燃料组件施加表3中的加速度功率谱，分析结果显示米塞斯应力均方值(RMISES)最大的一根燃料棒位于组件中间偏上的位置，最大值为2.90 MPa，出现在该燃料棒端部，如图9所示。在2倍加速度功率谱激励下，最大RMISES出现的位置相同，其值为4.35 MPa，仍远小于包壳材料的屈服强度。因此在正常运输工况下，新燃料组件出现包壳变形的风险极小。
应力最大位置处的功率谱密度(PSD)如图10所示，从中可以看出PSD的主峰在62.5 Hz。
Figure 9. The fuel rod with maximum RMISES
图9. RMISES最大的燃料棒
Figure 10. Power spectral density at the location of maximum stress
图10. 应力最大位置处的功率谱密度
对应力最大的单根燃料棒进行模态分析，得到的前11阶模态如表4所示(其中频率相同的模态算做一个)。结合PSD图可以看出，PSD主峰对应第7阶模态。第7阶的模态图如图11所示，这也解释了最大应力出现的位置。
Table 4. First 10th order modes of the single fuel rod
表4. 单根燃料棒前10阶模态
模态阶数 |
频率(Hz) |
1 |
34.090 |
2 |
36.720 |
3 |
41.112 |
4 |
46.800 |
5 |
52.887 |
6 |
58.717 |
7 |
62.515 |
8 |
105.99 |
9 |
123.27 |
10 |
130.50 |
Figure 11. 7th order modal diagram for the single fuel rod
图11. 单根燃料棒第7阶模态图
4.2. 燃料芯块
燃料芯块在4 g脉冲加速度的作用下，应力情况如图12所示。可以看出，随着脉冲持续时间的增加，芯块的应力水平也随之增加。当脉冲持续时间达到5 ms时，芯块的最大应力为2140 MPa，已经非常接近芯块的屈服强度，此时芯块有一定的开裂风险。
需要说明的是，简化的质量点模型无法考虑芯块与包壳之间的摩擦力，轴向摩擦力有可能减弱芯块受到的轴向冲击力，这一点有待后续深入研究。
5. 总结
本文对新燃料运输加速度载荷下的燃料棒及燃料芯块的受力状态进行了三维有限元分析，主要成果
Figure 12. Stress distribution of pellets under impulse loading
图12. 芯块脉冲载荷下的应力分布
可以归纳为以下几点：
1) 对燃料棒进行了精细建模，并开展了弯曲载荷下的受力分析，在弯曲载荷下，包壳会先于芯块发生损伤，因此以包壳弯矩–位移曲线为基准，对燃料棒进行了等效梁模型的简化，大大提高计算效率；
2) 对燃料组件进行了动力学响应分析，计算结果表明，在正常的运输载荷下，燃料棒包壳最大应力值远未达到材料的屈服强度，发生形变的风险极小；
3) 开展了燃料芯块在轴向脉冲加速度作用下的受力分析，当燃料芯块收到幅值为4 g，持续时间为5 ms的脉冲加速度时，芯块最大应力为2140 MPa，存在一定的开裂风险。
基金项目
本研究由“压水堆运行核燃料组件零破损关键技术”研究课题资助。