湍流模型对5×5格架棒束通道流动传热数值模拟影响分析

晁嫣萌，杨立新，张玉相，庞铮铮

(1.北京交通大学 机械与电子控制工程学院，北京 100044；2.中科华核电技术研究院有限公司，广东 深圳 518026)

定位格架是反应堆燃料组件的重要部件，影响着堆芯的热工水力性能。开发自主知识产权的燃料组件需深入了解定位格架对热工水力特性的影响，相比于周期长、费用高的试验研究，CFD模拟已成为安全、快速的重要设计辅助手段之一，应用CFD进行燃料组件内定位格架对流动传热特性的影响分析具有重要工程价值和意义。

1997年，Lee等[1]运用非线性k-ε湍流模型在未作任何模型参数调整情况下开展了数值仿真，得出该模型很大程度上低估了湍流强度，说明湍流模型对CFD模拟结果有着重要的影响。2002年，Smith等[2]采用RNGk-ε模型对带搅混翼格架燃料组件内流场进行了研究，结果与实验数据吻合较好。2003年，Yadigaroglu等[3]对棒束通道的数值计算结果表明，标准k-ε模型在窄缝空间不能很好地预测湍流流动。2012年，Liu等[4]采用Fluent软件，利用简化的燃料组件棒束通道模型，计算了6种不同湍流模型下Nu的沿程分布和典型位置的周向分布，其结果表明，Realizablek-ε湍流模型结合近壁面函数模拟结果与实验结果吻合较好。Holloway等[5-6]对不同类型的格架进行了单相流动换热实验研究，对棒束通道内换热系数沿程及周向分布进行了测量，为后续开展CFD模拟的标定提供了重要实验参考数据。

目前研究表明，采用不同湍流模型会使相同CFD模型的模拟结果不同。本文针对中国广核集团研发的格架结构，建立包含5×5典型格架结构和棒束通道CFD分析模型，采用ANSYS CFX软件，通过模拟获得不同湍流模型下的燃料组件棒束沿程及周向换热系数分布，与其他学者的相关实验数据进行对比分析，通过对3个典型搅混效果评价因子的分析，探讨子通道内流动与换热的内在关系，同时对比不同湍流模型对结果的影响，以得到最适宜本文所研究定位格架及棒束通道内流动传热特性计算的湍流模型。

1 几何模型

定位格架几何结构复杂，因弹簧、刚突与燃料棒间是线接触或小弧面接触，使局部网格划分困难，网格质量差，对CFD模拟结果影响较大。本文简化方法是在弹簧、刚突与燃料棒间留出0.1 mm的缝隙，如图1所示，同时删除原条带上的一些定位孔等不影响流动特性的微小结构。

图1 接触部位的简化

2 网格生成

采用非结构化四面体网格加拉伸网格的形式，应用ICEM CFD完成模型网格划分。格架部分由于结构复杂使用非结构化网格，流道部分几何规则且长径比较大，故采用拉伸网格，网格区域划分如图2所示。在燃料棒与刚突、弹簧缝隙处采用线控制和面控制进行网格细化，如图3所示。在格架条带和燃料棒束壁面设置附面层网格以提高网格质量和效率，如图4所示。

图2 网格类型区域的划分

图3 格架局部表面网格示意图

图4 附面层局部网格细节图

本文对格架四面体网格大小、拉伸网格数量以及附面层网格参数进行了敏感性分析，最终确定的计算网格模型参数列于表1。根据计算结果对燃料棒表面的Y+(第1层网格中心到壁面的无量纲距离)取平均值，如表2所列。不同湍流模型Y+值均在3～5范围内，说明计算模型壁面网格能满足各湍流模型要求。

表1 网格数量统计

表2 不同湍流模型Y+值

3 计算模型

3.1 边界条件

计算边界条件与Holloway等[5-6]的实验工况保持一致。实验工况为常压，冷却剂水入口温度为20 ℃，加热棒表面热流密度为10 kW/m2，Re分别为28 000和42 000，对应的CFD计算模型具体边界条件设置列于表3，表中未提到的边界参数均按软件默认设置。

表3 边界条件设置

计算采用不可压缩稳态流动传热模型，选取6种湍流模型进行了对比分析，水工质物性取CFX软件IAPSW数据库参数，计算求解采用高阶精度差分格式，收敛标准为各物理量残差小于5×10-5，同时出口流量和平均温度不再发生变化。

3.2 湍流模型

本文进行了涡粘和雷诺应力两类湍流模型计算，其中包括涡粘模型中的标准k-ε、RNGk-ε和SST 3种模型以及雷诺应力模型中LRR-IP、LRR-QI和SSG 3种模型。涡粘模型是假设雷诺应力与平均速度梯度呈正比，雷诺应力模型未应用涡粘假设，而是求解流体雷诺应力输运方程。

标准k-ε模型目前在科学研究及工程实际中得到最为广泛的检验和应用，但其用于强旋流动和弯曲壁面流动时仍会产生一定的失真，RNGk-ε对湍流黏度进行了修正，SST模型考虑了湍流剪切应力的输运，对在逆压力梯度下流动的分离起点和发展能预测得很准。与涡粘性模式相比，雷诺应力模型考虑了流动中的旋转及旋转流动的影响，包含以下因素：流线曲率、变形率突变、二次流及浮力，更适用于应变场复杂的各向异性流场，尤其是存在大流线曲率或涡旋的流场。3种雷诺应力模型由不同学者发展，有各自不同的模型常数。

6种湍流模型计算时，除SST外，均采用了默认的Scalable壁面函数，SST模型使用自动壁面函数。

4 结果分析

本文以Holloway实验结果为参考，对比分析了不同湍流模型CFD计算结果。Holloway实验格架为不带弹簧刚突的撕裂式搅混翼结构，本文CFD模型格架为带弹簧刚突的分开式搅混翼结构，图5为两种格架结构俯视图。通过数值结果与实验结果的对比，分析了格架不同搅混翼类型对流动的影响，同时在模型结构完全一致的棒束通道区域进一步验证了不同湍流模型的计算结果。下面从沿程压降、流动搅混性能评价因子以及Nu分布对比分析不同湍流模型的计算结果。

图5 不同搅混翼格架俯视图

4.1 沿程压降

图6为两种雷诺数工况下格架和棒束通道计算压降与实验结果对比。图6a示出了格架段的压降和棒束通道的压降测点位置，图6b、c分别为两种雷诺数工况下的各湍流模型计算压降与实验结果比较情况，图中两条水平直线位置分别代表实验的两段压降测量值。

由图可知，对于格架段压降，6种湍流模型结果均高于实验值，雷诺应力模型和实验值较吻合，涡粘模型结果偏高；对于棒束区压降，除SST模型外，其余数值预测结果均低于实验值，雷诺应力模型平均低20%。仅就棒束区域压降计算结果看，标准k-ε和SST模型与实验最接近。对比CFD模型和实验模型格架结构，可知格架段压降计算值偏高显然是由于弹簧/刚突和搅混翼结构不同而引起的，因格架段压降中也包含一完整棒束区压降，考虑模拟结果在棒束区的偏低量，可推算出由于格架结构不同所引起的压降约占总压降的15%。

4.2 流场搅混特性

大量文献[7-10]报道了格架搅混翼能在下游一定距离的横截面产生并保持一定强度的涡流和交叉流。存在于子通道的涡流引起的离心力能带走燃料棒壁面形成的汽泡，存在于子通道间的交叉流可均衡子通道间的焓升。不同形式的评价因子被定义来研究搅混翼所引起的流动特性，虽然其具体形式上略有差别，但所描述的流动本质相同，这些因子可概括为涡流搅混因子、交叉流搅混因子及湍流强度因子，这3个是影响燃料棒表面换热特性及临界热流密度的重要流动参数。

涡流搅混因子Svortex可用来衡量涡流对流体的搅混作用，定义为：

其中：R为子通道中心到燃料棒表面的距离；r为通道截面上流域内各点到中心的距离；vt为横截面上引起涡流的横向速度；ua为横截面上各点的轴向速度。

图6 沿程压降

交叉流搅混因子Fcross用来衡量子通道间交叉搅混的强度，定义为：

其中：s为燃料棒表面间距(栅距-棒直径)；Vcross为穿过棒间截面的交叉流速度分量；Ubulk为截面平均轴向流速。

湍流强度因子Tt为：

其中：k为湍动能；U为子通道内轴向平均流速。

取中心棒束左上角的子通道对上述3个因子进行积分计算，图7示出了6种湍流模型3个因子在两种雷诺数工况下的沿程变化曲线。图中横坐标采用了轴向距离与燃料棒直径比值的无量纲参数，Dh为燃料棒直径，零点位置为格架出口位置。由图7可知，两种Re工况下的搅混和交叉流因子在数量级及沿程分布趋势上均一致，湍流强度因子分布趋势一致，在高Re工况下该值略低。

图7a、b表明，涡流搅混因子受搅混翼影响在刚出格架时值最大，在0～3Dh区间迅速衰减，3Dh后又开始增加并在10Dh左右达到第2个峰值，这区间受搅混翼的影响显著，10Dh～25Dh区间搅混翼的影响逐渐减弱，其值波动衰减，25Dh后该因子已降至较低水平，并继续沿程线性下降，搅混翼的影响基本消失。各湍流模型结果趋势一致，除SST模型值略小外，其余曲线基本重合。

图7c、d表明，交叉流搅混因子刚出格架时较小，在搅混翼作用下迅速上升，并在3Dh处达到最大值，3Dh～10Dh间波动下降，该区间交叉搅混因子受搅混翼的影响显著，10Dh～25Dh区间搅混翼的影响逐渐减弱，该因子呈波动衰减，25Dh后该因子也降至较低水平，并沿程线性下降。雷诺应力3个湍流模型结果基本重合，涡粘模型中SST模型与其他结果偏差略大。结合图7a～d，可看到在涡流搅混因子较高的位置交叉流搅混因子较低(如格架出口及10Dh处)，同时当涡流搅混因子下降时交叉流搅混因子则有升高的趋势(10Dh～25Dh区间)，可见搅混翼引起的绕子通道中心涡流强度与子通道间交叉流强度相互耦合。

图7e、f表明，湍流强度因子出格架后迅速下降，在3Dh～10Dh间略有波动，5Dh～25Dh区间基本维持不变，在下游后段又呈现上升趋势。SST和标准k-ε模型湍流强度因子明显大于其他湍流模型的，这也间接说明了这两个湍流模型计算压降较大的原因。Re大的计算工况湍流因子量级略有减小，原因是Re增大导致轴向速度U增大，而湍动能k并不完全与U的二次方呈正比。

图7 不同湍流模型3种评价因子沿程变化

4.3 换热分析

1) 沿程平均Nu

取5×5格架中间燃料棒束表面，沿流动方向截取系列截面，对每个截面棒束圆周的Nu进行算术平均，绘制沿程平均Nu沿程变化曲线，如图8所示。图中分别示出了两种雷诺数工况条件下6种湍流模型计算结果与Holloway实验结果的对比曲线。

由图8可知，两种雷诺数工况下，随Re的增大棒束表面平均Nu也增大，但沿程变化趋势相同。6种湍流模型中，SST模型的计算结果明显高于其他模型和实验结果，这是由于SST模型采用的自动壁面函数对壁面网格更敏感导致的。在棒束通道区，涡粘类模型的Nu均高于雷诺应力模型的结果，在高雷诺数情况下该趋势更明显。雷诺应力类的3种湍流模型Nu基本重合，且在受格架影响较弱的棒束通道内其值与实验结果较为接近。

由图8还可知，6种湍流模型计算得到的平均Nu均是在格架出口位置最大，然后迅速减小，在3Dh～10Dh区间出现剧烈波动，在10Dh～25Dh区间呈现波动下降趋势，25Dh后Nu趋于稳定值。沿程Nu变化区间的特征位置与上节中分析的3个流动因子变化区间位置是相对应的，说明搅混翼引起的搅混流动直接影响着棒束表面的换热特性。Holloway实验中Nu曲线在离开格架位置后的0～5Dh区间迅速降低，然后有小幅上升，并在25Dh后Nu分别稳定在210和280。对比可知，在25Dh位置前，计算结果与实验值有较大差异，这是由于本文CFD模型采用了分开式搅混翼格架，搅混翼长度、宽度等形状参数与实验所用的撕裂式搅混翼格架不同(图5)，计算结果表明分开式搅混翼会引起更强的涡流和交叉流，从而提高了格架下游子通道的换热能力。经过25Dh后，流体进入棒束通道，此时格架和搅混翼对流动的影响已大幅减弱，而棒束通道实验段和CFD模型完全相同，所以这一区间计算与实验Nu的对比能用来更直接地验证模型的准确性。综合分析不同Re下计算工况的结果，雷诺应力类的3种湍流模型计算的Nu与实验值吻合得更好。

2) 棒束周向Nu

Holloway实验给出了2.2Dh、6.5Dh和35.7Dh位置的Nu周向分布，取计算模型的中间燃料棒束对应位置截取3个截面，图9左侧示出了这3个位置棒束表面不同湍流模型计算Nu的分布，右侧示出了对应位置截面上的横向流速度分布矢量图。

由图9左侧的Nu曲线对比可知，除SST模型外，其他湍流模型预测值均与实验值较接近。2.2Dh和6.5Dh位置的Nu周向分布实验值与计算值在分布的波形上相位略有区别，这是由于计算模型弹簧刚突的存在以及搅混翼结构的不同所引起的。Nu分布明显受到搅混翼所引起的横向速度分布的影响，对照右侧的横向速度分布矢量图可发现：Nu较大的位置横向速度较大(如棒束圆周105°、345°等位置)，而Nu处于波谷位置的横向速度较小或局部有涡流出现(如棒束圆周45°、145°等位置)。在35.7Dh位置5种湍流模型的预测值与实验值已基本重合，搅混翼引起的横向流动影响也基本消失。

图8 不同湍流模型Nu沿程分布与实验值

图9 不同湍流模型Nu周向分布

5 结论

1) 涡粘模型中标准k-ε和SST模型对压降预测较好，而对换热系数的预测结果与其他湍流模型和实验值偏差较大，其中SST模型对Nu的过高预测与CFX中该模型采用的自动壁面函数有关；

2) 格架及搅混翼结构给流动子通道内带来强烈的涡流和交叉流，雷诺应力类湍流模型更适用于这种各向异性的流场模拟，3种雷诺应力模型结果基本一致，Nu与实验结果吻合较好，是后续进行这类格架详细性能分析及设计优化的首选湍流模型；

3) 搅混翼引起的沿程涡流搅混因子和交叉流搅混因子具有强烈的耦合关系，并共同影响着沿程平均Nu分布，这种影响趋势在格架下游10Dh后逐渐减弱，在25Dh后基本消失。

参考文献：

[1] LEE K B, JANG H G. A numerical prediction on the turbulent flow in closely spaced bare rod arrays by a nonlineark-εmodel[J]. Nuclear Engineering and Design, 1997, 172: 351-357.

[2] SMITH L D, CONNER M E, LIU B, et al. Benchmarking computational fluid dynamics for application to PWR fuel[C]∥Proceedings of the 10th International Conference on Nuclear Engineering. USA: ASME, 2002: 823-830.

[3] YADIGAROGLU G, ANDERANI M, DREIER J, et al. Trends and needs in experimentation and numerical simulation for LWR safety[J]. Nuclear Engineering and Design, 2003, 221: 205-223.

[4] LIU C C, FERNG Y M, SHIH C K. CFD evaluation of turbulence models for flow simulation of the fuel rod bundle with a spacer assembly[J]. Applied Thermal Engineering, 2012(40): 389-396.

[5] HOLLOWAY M V, McCLUSKY H L, BEASLEY D E, et al. The effect of support grid features on local, single-phase heat transfer measurements in rod bundles[J]. ASME Journal of Heat Transfer, 2004, 126: 43-53.

[6] HOLLOWAY M V, CONOVER T A, McCLUSKY H L, et al. The effect of support grid design on azimuthal variation in heat transfer coefficient for rod bundles[J]. ASME Journal of Heat Transfer, 2005, 127: 598-605.

[7] CUI X Z, KIM K Y. Three-dimensional analysis of turbulent heat transfer and flow through mixing vane in a subchannel of nuclear reactor[J]. Journal of Nuclear Science and Technology, 2003, 40(10): 719-724.

[8] LEE C M, CHOI Y D. Comparison of thermo-hydraulic performances of large scale vortex flow (LSVF) and small scale vortex flow (SSVF) mixing vanes in 17×17 nuclear rod bundle[J]. Nuclear Engineering and Design, 2007, 237(24): 2 322-2 331.

[9] IN W K, CHUN T H, OH D S, et al. CFD analysis of turbulent flows in rod bundles for nuclear fuel spacer design[C]∥Transactions of the 15th International Conference on Structural Mechanics in Reactor Technology. Seoul: SMiRT, 1999: 365-372.

[10] IN W K, CHUN T H, SHIN C H, et al. Numerical computation of heat transfer enhancement of a PWR rod bundle with mixing vane spacers[J]. Nuclear Technology, 2008, 161(1): 69-79.