基于CFD方法的铅铋冷却燃料棒束的热工水力特性分析

孙 畅，焦守华，柴 翔，杨 云

(上海交通大学 核科学与工程学院，上海 200240)

液态铅铋合金具有良好的中子学性能和较高的中子经济性，使反应堆可在常压下运行，同时实现高功率密度。这不仅大幅降低了由高压控制引起的失水事故的发生概率，同时也大幅降低了一回路设备对承压性能的要求，是我国未来铅基反应堆的重要材料之一[1]。铅铋冷却快堆燃料组件为六边形密集棒束结构，其采用了具有强化换热和减小机械震荡效果[2]的金属绕丝进行固定。

由于液态铅铋合金的高密度、强腐蚀及实验温度要求较高等因素，开展液态铅铋合金在燃料组件内的热工水力实验难度较大。此外由于液态铅铋合金的普朗特数Pr很低(约0.025)，若使用相似准则法在常规介质中进行模化实验也存在一定的挑战。综合实验技术、设备及经费等原因，目前国外较大规模地开展液态铅铋合金在燃料组件内的热工水力实验的科研机构仅有德国的KIT和意大利的ENEA实验室。德国KIT的KALLA实验室[3]开展了液态铅铋合金在燃料棒束间的热工水力特性研究实验，以研究冷却剂的流动阻力特性和流动传热特性。Litfin等[4]开展了液态铅铋合金和水冷却剂在19棒束结构内的流动阻力特性实验，并采用欧拉相似准则进行了关联比较。Pacio等[5]开展了使用格架定位的19棒束燃料组件内的液态铅铋合金流动传热行为实验，通过实验测得的压降和流动速度估算铅铋合金在组件内的流动特性，最终得出“液态铅铋合金与水等常规流体介质具有相似的流动特性”这一结论。意大利的Brasimone研究中心[6-7]对液态重金属作为反应堆冷却剂的可能性展开了研究，针对液态重金属堆型相关设备材料耐高温、抗腐蚀性测量，堆芯组件内的热工水力等开展了一系列实验。

针对Pacio等[8-9]设计的带绕丝的19棒束铅铋组件，本文使用商用计算流体力学(CFD)软件STAR-CCM+对其建立全尺寸模型，划分网格后选取数值模型进行模拟计算，将该组件中子通道内的流动传热与实验值进行对比验证，之后对组件的进口质量流量和总热功率进行敏感性分析。

1 正常工况模拟

本文以Pacio等[8-9]设计的组件为研究对象，组件由截面为正六边形的外套管与其中呈正三角形排列的19根被金属绕丝缠绕的加热棒组成。燃料组件几何参数列于表1。

表1 燃料组件几何参数Table 1 Geometric parameter of fuel assembly

1.1 几何模型

图1 燃料棒和子通道编号示意图Fig.1 Diagram of number of rod and sub-channel

1.2 网格划分

网格划分使用STAR-CCM+中自动生成非结构化网格，该工具可自动根据局部区域的几何特征修改附近网格尺寸，根据壁面网格质量选取不同的y+值。为尽可能提高模拟精度，开展了网格敏感性分析。网格划分时选取3种网格数量进行敏感性分析，在不同网格数量下计算得到1号加热棒平均温度在轴向的发展趋势，如图2所示。计算发现，3种网格数量在相同的计算条件下温度发展趋势相同，温度模拟值非常相近。为保证计算机可负荷的情况下尽可能提高网格质量，采用4 000万网格尺寸进行计算。各计算区域的网格细节如图3所示。

图2 网格敏感性分析Fig.2 Mesh sensitivity analysis

1.3 湍流数值模型

在各湍流数值模型中，k-ε模型能较好地模拟充分发展的湍流流动，k-ω模型则能较好地模拟各压力梯度下近壁面处的流动情况。Menter[11]提出SSTk-ω模型，不仅在近壁面处保留原始的k-ω模型，也对充分发展的湍流流动状况运用k-ω模型，k方程和ω方程为：

(1)

(2)

(3)

(4)

(5)

其中：ρ为密度；Gk为湍流动能k的产生项；Gω为比耗散率ω的产生项；t为时间；x为三维坐标系；μ为动力黏度；μt=ρkT，T为湍流时间尺度；σk、σω2、β*和β可通过混合函数求解得到；γ为混合系数；S为流体平均应力张量的模；U为流体速度；y为距壁面的距离；ν为流体的运动黏度；CDkw为与流体的横向扩散项相关的项。

a——整体网格；b——局部网格图3 流体区域和固体区域的网格Fig.3 Mesh around fluid region and solid region

由于铅铋合金普朗特数较低、热导率较高，与水冷却剂相比，铅铋合金在包壳表面形成的热边界层薄于流动边界场，故该模拟计算中采用SSTk-ω湍流模型。

湍流普朗特数Prt表现了湍流流动的一种属性，其数值为湍流动量扩散率与湍流热扩散率的比值，本文选取了不同学者关于Prt的研究结果[12-14]。

根据不同Prt关系式计算得到的1号加热棒平均温度T在轴向上的发展状况，如图4所示。计算发现，在相同的计算条件下的温度发展趋势相同，温度模拟值差异不大。由于Cheng等[14]关于Prt的研究更适用于恒定热流密度情况下的铅铋冷却剂流动计算，故选择其研究的Prt关系式。

图4 湍流普朗特数敏感性分析Fig.4 Turbulent Prandtl number sensitivity analysis

1.4 边界条件

在计算过程中，燃料组件的总功率为197 kW。

经历过17世纪“黄金时代”的阿姆斯特丹，如今的倚重科技创新，重塑城市生活方式，兼顾城市发展与居民生活质量，成为欧洲乃至全球市场中最具竞争力绿色城市之一，是欧洲最安全的城市。

入口边界条件为：质量流量为19 kg/s，冷却剂进口温度为473.15 K。组件外围壁面和上、下两个端面设为绝热边界条件以忽略铅铋合金对壁面间的热传递，另外将其余面设为交界面以实现不同计算区域间的热传递。压力设为出口边界条件，在加热棒内壁面处给定热流密度，热流密度qw由下式计算：

(6)

其中：N为加热棒数量；Q为总功率。

2 模型验证

在Pacio等实验中，若干数量的热电偶布置在沿轴向发展的3个不同位置处，其对应高度分别为绕丝螺距H的1/6、11/6和15/6倍。将各温度测量点的模拟值与实验结果进行对比，以验证几何建模与数值模型的适用性，结果如图5所示。图5中，Rod1表示1号加热棒，SCH-1表示1号子通道。

由图5可见，随着冷却剂充分发展，温度模拟值与实验结果的差距减小。在众多模拟值与实验结果的对比中，子通道温度的最大相对误差为4.12%，加热棒温度的最大相对误差为2.84%。由于增加考虑固体域内的导热，本文结果与Chai等[15]的结果相比，尤其在子通道处，模拟值与实验结果更为接近。

时，Pe与Nu的关系图5 不同高度处加热棒温度和子通道温度的实验结果与模拟值对比Fig.5 Comparison of temperature at different heights between simulation value and experimental data

3 燃料棒束热工水力特性分析

3.1 分析算例

本文旨在研究绕丝棒束关于进口质量流量和功率的敏感性，拟定的质量流量和功率工况列于表2。当功率为197 kW时，调整组件入口的质量流量，分别为Pacio等实验参数的0.25、0.5、1和2倍。当组件入口质量流量为19 kg/s时，将组件功率设为Pacio等实验参数的0.6、0.8、1和1.2倍。

3.2 结果分析

垂直于轴向、距进口2.25倍螺距(z=738 mm)和2.5倍螺距(z=820 mm)时，工况3下截面横向流速度分布如图6所示。由图6可见，随绕丝沿轴向高度的发展，冷却剂在子通道内发生交混，使得截面上的速度峰值从左上区域向右上区域转移。由于边角子通道的当量直径较大，因此摩擦阻力较小，该区域的冷却剂速度则相对较大。绕丝的存在使得子通道间隙宽度减小，铅铋合金的流通面积减小，则其流通速度增大。

表2 拟定的质量流量和功率工况Table 2 Selected case of mass flow rate and power

图7为与图6对应高度的温度截面，由于速度的不均匀性分布使得横截面上的温度分布呈不均匀状态。在边角通道周围加热功率较少，且铅铋合金速度较大，铅铋合金得不到充分加热而呈现出较低的温度分布。而在速度较低的子通道区域，由于较低的铅铋合金速度分布，铅铋合金可得到充分加热，故图7中中心区域的温度较高。

a——z=738 mm；b——z=820 mm图6 工况3下的截面横向流速度分布Fig.6 Velocity distribution of cross section transverse flow in case 3

a——z=738 mm；b——z=820 mm图7 工况3下的截面横向流温度分布Fig.7 Temperature distribution of cross section transverse flow in case 3

图8示出工况3下截面横向流速度矢量。由图8a可看出，边缘子通道左上部区域中冷却剂的流动强度明显低于右下部区域，图8b中边缘子通道左上部区域中冷却剂的流动强度明显高于右上区域，这是因为绕丝在加热棒上沿轴向方向呈逆时针向上旋转，不同高度的截面上横向流动受到方位角的影响，使速度的方向随绕丝的转动而发生变化，最大速度位置也发生周期性旋转。在同一平面中，铅铋合金横向流速度在中心通道的分布较为平缓，而在边角位置的变化较剧烈。图9示出冷却剂流动过程中的流线图。

a——z=738 mm；b——z=820 mm图8 工况3下的截面横向流速度矢量Fig.8 Velocity vector of cross section transverse flow under case 3

图9 工况3下的速度流线图Fig.9 Streamline plot in case 3

考虑到棒束的对称性，如图10所示，选取棒1、4、11、12研究温度随轴向高度的变化。对于处在棒束中心区域的燃料棒，如棒1、4，其表面平均温度沿轴向的发展呈稳定的周期性震荡，这是因为绕丝的存在使中心通道间的铅铋合金呈稳态交叉流动，从而使温度呈现稳定的周期性震荡。对于处在棒束边缘区域的加热棒，如棒11、12，边角子通道内的二次流形态较为复杂，故未观察到加热棒温度沿轴向的周期性变化。

图10 工况3下不同加热棒的温度沿轴向的变化Fig.10 Temperature of different heating rods vs. axial position in case 3

将沿加热棒周向的平均温度作为包壳温度，研究质量流量对加热棒温度分布的影响。选取棒1、4、11和12为研究对象，采用式(7)对每根加热棒的包壳温度进行无量纲化：

(7)

其中：T′为无量纲的加热棒包壳温度；T为加热棒包壳温度；Tc、Th分别为加热棒包壳温度的初始值和最大值。

在总功率相同的条件下，不同质量流量下的温度分布如图11所示。由图11可见，质量流量的增加使得进口位置处加热棒的温度变化更加剧烈，随铅铋合金在轴向的发展，质量流量对加热棒中后段的温度变化幅度影响较小。

为研究燃料组件内冷却剂的流动特性，将子通道i与子通道j交界面中心位置处的横向速度定义为ωi，j(i

由图12可见，由于绝热外壁的存在和同一截面中绕丝相对交界面的方向变化，ω2，3和ω27，28分别在轴向呈现出不同的周期性变化。在划分网格过程中，由于绕丝向加热棒圆心处移动0.1 mm，铅铋合金流经该区域时，其流通面积急速降低，冷却剂横向交混速度则迅速增加，形成异常峰值。另外，随铅铋合金进口质量流量的增大，其交混程度变得更剧烈。

a——棒1；b——棒4；c——棒11；d——棒12图11 加热棒温度的质量流量敏感性分析Fig.11 Sensitivity analysis of mass flow rate for heating rod temperature

图12 ωi，j沿轴向高度随质量流量的变化Fig.12 ωi，j change along axial height with mass flow rate

在入口质量流量均为19 kg/s的条件下，研究加热棒温度的功率敏感性，如图13所示。图13结果与图10的相同，功率大小只影响加热棒温度的大小，但功率变化并不影响加热棒平均温度随轴向的变化。图14示出ωi，j沿轴向高度随功率的变化。由图14可见，ω2，3和ω27，28在不同功率下均呈现出相同的周期性变化，由此可见功率大小对冷却剂的横向交混速度影响较小。

a——棒1；b——棒11图13 加热棒温度的功率敏感性分析Fig.13 Power sensitivity analysis of heating rod temperature

图14 ωi，j沿轴向高度随功率的变化Fig.14 ωi，j change along axial height with power

4 结论

本文基于商用计算流体力学软件STAR-CCM+，使用流固耦合的方法对带有绕丝结构的19棒束铅铋组件建立几何模型，选取适当的网格数量，确定合适的湍流模型。模拟计算结果表明，由于边角子通道与中心通道间的摩擦阻力存在差异，使得铅铋合金在横截面上的速度呈不均匀性分布，进而导致温度分布的不均匀性。在进口质量流量为19 kg/s和总功率为197 kW工况下，中心区域和边缘区域的加热棒温度分别呈现稳态周期性震荡和非稳态震荡。另外，随质量流量的增加，子通道间横向交混增大。功率变化对通道间的横向交混速度的影响较小，对冷却剂温度的横向分布也无明显差异。