U-50Zr螺旋十字燃料热力耦合特性分析

蔡孟珂， 丛腾龙， 顾汉洋

(上海交通大学 机械与动力工程学院， 上海 200240)

为了提升现有压水堆的功率水平和安全裕量，研究者多从燃料材料和燃料几何2个角度对传统核燃料进行改进[1-5]。在材料方面，U-50Zr合金材料具有较高的热导率、改进的抗辐照性能、较好的包壳相容性以及在闭式循环中的优势，因此极具潜力[6]。在核燃料服役的过程中，气体裂变产物的堆积使得燃料晶胞内和晶界间形成微小气泡或孔隙，进而改变燃料的宏观热力性质。目前，针对传统UO2陶瓷燃料的孔隙效应的研究已较为完备；然而，针对U-50Zr合金材料的研究却几乎为空白。为了准确地预测U-50Zr合金核燃料的堆内辐照性能，需要对其材料孔隙效应进行研究。在燃料几何方面，螺旋十字燃料(helical cruciform fuel，HCF)是一种结合了棒状与板状燃料元件优点的革新型设计。HCF燃料的横截面为十字型，其翼片沿轴向方向扭转形成螺旋结构，可提高燃料的换热面积，并增强流体的搅浑换热性能。同时，燃料的翼根区离燃料中心的径向距离较短，减小芯体导热路径，降低燃料中心温度。另外，在轴向每扭转1/4周期平面内，相邻燃料间的翼尖区相互接触，形成“自定位”结构，故可省略定位格架，并在一定程度上弥补燃料格架带来的压降损失。

在新型核燃料的设计和优化中，燃料的热工水力特性以及热力响应反馈分析极为重要。目前，对于HCF的研究主要聚焦于热工水力领域，主要包括HCF的压降损失，搅浑特性以及传热能力。文献[7-8]通过实验和数值模拟的方法，对HCF的压降特性进行了研究。文献[9-14]通过视网探针和CFD模拟技术对HCF棒束的横向搅浑性能进行了研究，并预测了其在单相和两相情况下的传热特性。李惠悦等[15]采用SSTk-ω湍流模型对不同螺距的星形螺旋燃料棒束通道内流动换热特性进行数值模拟研究。Shirvan[16]基于两相欧拉方法研究了高压下HCF的沸腾危机现象。除了上述热工水力研究之外，燃料热力性能分析也对核燃料的设计和优化具有重要意义。然而，目前针对HCF的燃料性能分析的研究十分稀少，邓阳斌等[17]对单根HCF的热力性能分析进行了初步预测，但是，由于HCF复杂的三维几何结构，需要对多棒束下HCF的热力响应特性进行更加深入和全面的研究，以捕捉包壳翼尖处的接触应力特性。

本文从材料和燃料分析2方面出发，对U-50Zr合金材料HCF进行多物理场耦合分析。在材料方面，本文运用分子动力学软件LAMMPS，建立了带空隙U-50Zr合金材料的半经验修正模型。在燃料分析方面，本文考虑了包括U-50Zr材料孔隙退化模型在内的多物理场效应，对较深辐照下HCF的热力响应特性进行了耦合模拟和分析。

1 热力耦合数学物理模型

本文从材料和燃料分析2方面出发，对U-50Zr合金材料HCF进行多物理场耦合分析，如图1所示。将分别给出U-50Zr合金弹性模量、热导率、热力耦合计算模型。

1.1 杨氏模量的计算

为了获取不同温度下的U-50Zr合金材料的弹性模量，常用的方法是进行单轴拉伸实验模拟。对大多数材料而言，杨氏模量可看作是各向同性的。然而，U-50Zr合金是由δ相UZ2组成的，该δ相UZ2是六方点阵排布的[18]，因此U-50Zr合金的弹性模量是各向异性的。为了获取更准确的弹性模量预测值，单轴拉伸实验需在相互正交的3个方向上进行。在三维正交坐标系下，应力和应变的关系为：

(1)

式中：下标11、22、33、12、13和23分别表示x、y、z、xy、xz、yz方向；Cij为弹性模量张量。为了便于工程使用，往往将该6×6弹性模量矩阵转化成等效杨氏模量，其转化关系为：

(2)

(3)

(4)

(5)

(6)

(7)

式中：K是体积模量；G是剪切；E是杨氏模量。

图1 基于U-50Zr合金材料下螺旋十字燃料的多物理场耦合Fig.1 The multi-physics coupling analysis for HCF based on U-50Zr alloy

1.2 热导率的计算

U-50Zr合金的热导率kt由声子导热率kp和电子导热率ke2部分组成：

kt=kp+ke

(8)

本文采用非平衡分子动力学方法预测材料的声子导热率，其原理由图2所示。在长平板的两端施加较冷的温度边界，在平板的中段施加较热的温度边界，通过测量平板材料的温度梯度和热流密度，进而由傅里叶导热定律求出声子导热率为：

(9)

式中：vh和vc是热端和冷端的原子平均速度；t是导热时间；S是平板端面截面积大小；dT/dx是温度梯度。

图2 非平衡分子动力学原理求声子热导率示意Fig.2 The sketch of the NEMD method to obtain phonon conductivity

同时，电子导热率跟电子阻抗值有关，二者的关系可通过Wiedemann-Franz定律描述：

(10)

式中：L是洛伦兹常数；ρe是电子阻抗。本文电子阻抗的计算借用了Zhou等[19]提出的经验模型：

(11)

式中：B是电子-声子散射时间；C是电子-电子散射时间；ρres是残余电阻率；ρsat是饱和电阻率。

1.3 热力耦合计算

热力耦合计算由瞬态热传导方程和力学本构方程组成。其中，热传导方程由时间项，扩散项和源项构成：

(12)

(13)

(14)

式中：σ是柯西应力；Fv是单位体积力；E是弹性模量；v是泊松比；ε是弹性应变。在非弹性变形阶段，总应变由弹性应变、热应变、塑性应变、肿胀应变和蠕变应变5部分组成。为了实现完全热力耦合，导热和变形方程同时求解，并同时考虑热对力、力对热的影响，故全局控制方程为：

(15)

式中：θ是节点温度；u是节点位移；Kθθ是温度刚度矩阵；Kuu是位移刚度矩阵；Kθu是位移对温度影响的刚度矩阵；Kuθ是温度对位移影响的刚度矩阵；Rθ和Ru分别是温度和位移方程的残差。

2 热力耦合结果分析

2.1 材料孔隙效应分析

本文预测了U-50Zr合金在温度范围为300～800 K的等效杨氏模量。结合裂变气泡在燃料中的尺寸大小，孔隙尺寸半径设置为0.5～1.2 nm，对应的孔隙率大小为0.37%～7.99%，详见表1。图3给出了不同温度、不同孔隙率下，杨氏模量的变化趋势图。由图可见，随着温度和孔隙率的升高，U-50Zr合金的杨氏模量逐渐减小。

表1 杨氏模量计算中孔隙半径和孔隙率的设置

图3 不同温度和孔隙率下U-50Zr合金的杨氏模量云图Fig.3 The contour of Young’s modulus of U-50Zr at different temperatures and porosities

同理，本文预测了U-50Zr合金在温度范围为300～800 K的热导率，孔隙尺寸和孔隙率的设置情况见表2。图4给出了不同温度、不同孔隙率下，热导率的变化趋势图。由图4可见，随着孔隙率的增大，U-50Zr合金的热导率逐渐降低；但是，随着温度的升高，U-50Zr合金的热导率逐渐增高，这是因为高温提高了电子导热率。

表2 热导率计算中孔隙半径和孔隙率的设置

图4 不同温度和孔隙率下U-50Zr合金的热导率云图Fig.4 The contour of thermal conductivity of U-50Zr at different temperatures and porosities

为了更加准确地表述杨氏模量、热导率和孔隙率的关系，本文建立了带空隙U-50Zr合金材料的半经验修正模型为：

Ep=(150.1-0.052 7T)(1-9.5P)

(16)

kp=(2.245 71-7.73×10-3T+

1.071 43×10-5T2)(1-2.15P)

(17)

该空隙效应模型可初步预测U-50Zr合金的杨氏模量和热导率，其适用于低空隙率(p<10%)和纳米级气泡(气泡半径为纳米量级)。

2.2 燃料热力响应分析

本研究基于有限元软件Abaqus的Standard隐式求解器，综合考虑基于U-50Zr的HCF多物理场现象：热膨胀、辐照肿胀、蠕变、弹塑性变形等，并植入了式(16)和式(17)中的孔隙效应模型，对深辐照下HCF的热力响应特性进行分析。HCF在线功率水平18.2 kW/m下辐照了2.4 a。若以发生裂变的原子数占总原子数的百分比(per fissions of initial mental atoms, FIMA)为燃耗单位，其燃耗深度达到了14.1% FIMA。图5给出了典型HCF燃料的几何结构。另外，由于HCF具有翼片结构，因此翼根区域的导热路径小于翼尖区域，因此大部分热量从翼根区域出流出，造成翼根处的热流密度远大于翼尖处。

图5 螺旋十字燃料3维几何构型Fig.5 The 3D geometry of helical cruciform fuel

在传热方面，图6和图7分别给出了未辐照情况下，HCF的温度和热流密度分布云图。由于U-50Zr体系下的HCF具有较高的热导率和换热性能，因此HCF的整体温度偏低，在轴向250 mm处的活性区顶端，最高温度为674.4 K，远低于常规陶瓷棒状燃料。

图6 未辐照下螺旋十字燃料的温度云图Fig.6 Temperature contour of unirradiated HCF

图7 未辐照下螺旋十字燃料的热流密度云图Fig.7 Heat flux contour of unirradiated HCF

图8为包壳外壁面的温度和热流密度的周向分布图。显然，在45°、135°、225°、315°(即翼根)处的温度和热流密度均大于翼尖处。其中，温度在翼根处达到最大值583 K，在翼尖处达到最小值564 K；热流密度在翼根处达到最大值1 150 kW/m2,在翼尖处达到最小值281 kW/m2。图9为燃料中心温度随燃耗变化的影响。随着燃耗的加深，燃料产生中的裂变气体不断积累，进而使燃料产生孔隙，减小燃料的热导率。但是由图9看出，燃料中心温度随燃耗的变化不明显，从未经辐照到辐照末期，温度提升不超过0.3 K，这充分体现了U-50Zr体系下HCF的耐辐照稳定性。

图8 未辐照下螺旋十字燃料温度和热流密度周向分布图Fig.8 Circumferential distribution of temperature and heat flux of unirradiated HCF

在受力方面，图10 (a)给出了未辐照下，轴向扭转90°平面内HCF燃料的Mise应力分布情况，应力最大值位于产生接触的包壳翼尖区域。由图10 (b)所示，当燃耗达到14.1% FIMA时，应力最大值转移到包壳翼尖区域两侧，且应力集中现象较未辐照时有所缓解，应力峰值也相应降低。图11给出了轴向扭转90°平面内HCF燃料翼尖区域应力峰值随燃耗的变化情况。在燃耗从0～1.1% FIMA内，由于蠕变效应的影响，应力峰值从326.7 MPa下降到307.3 MPa。随着燃耗进一步的深入，辐照肿胀加大了燃料向外膨胀的趋势，因此使得接触处的应力峰值重新开始增大，当燃耗达到3.3% FIMA 时，应力峰值增大到322.1 MPa。此后，燃料应力在蠕变的影响下开始回落。最终，当燃耗达到14.1% FIMA 时，翼尖处应力峰值达到了313.3 MPa。在相同情况下，Geelhood等[20]提出的经验模型计算得到锆合金包壳的极限断裂强度值为500 MPa左右。因此，在较深燃耗下，U-50Zr体系下HCF的包壳应力最大值低于材料失效强度，表明其拥有较好的力学安全性能。

图9 螺旋十字燃料中心温度随燃耗变化Fig.9 Plot of central temperature of HCF

图10 螺旋十字燃料的Mises应力云图Fig.10 Mises stress contour of HCF

图11 螺旋十字燃料应力峰值随燃耗变化图Fig.11 Plot of peak stress of HCF under different burnups

3 结论

1)在低孔隙率、低温条件下，U-50Zr合金材料的杨氏模量随孔隙率、温度的升高而降低；U-50Zr合金材料的热导率随孔隙率的升高而降低，随温度的升高而升高。

2)螺旋十字燃料存在周向温度、热流密度分布的不均匀现象。

3)螺旋十字燃料的燃料温度较低，在整个计算寿期内，燃料中心温度较为稳定，且不超过675 K；另外，燃耗深度对温度的影响不显著。

4)U-50Zr材料的物性孔隙模型的准确性仍需要未来的辐照实验数据的验证和支撑。