棒束内湍流流动特性的PIV与CFD研究

王 汉，王舜琦，陆道纲

(华北电力大学 核科学与工程学院，北京 102206)

在压水堆燃料组件内，冷却剂的流动和传热特性直接影响包壳和燃料芯块的温度分布，对反应堆的安全性和经济性至关重要。定位格架是燃料组件的重要组成部分，不仅能支撑和定位燃料棒，而且具有增强局部扰动、强化包壳与冷却剂对流传热的作用。国内外众多研究者[1-5]对定位格架作用下燃料组件的热工水力特性进行了实验研究，结果表明，冷却剂流经定位格架交混翼后产生强烈的横向流动，增强了相邻子通道间的交混。由于这种交混效应，带有交混翼的燃料组件临界热流密度可较无交混翼的燃料组件提高20%以上[5]。然而，近期的实验研究[6-8]表明，定位格架引入的横向速度和横向湍流强度沿流动方向以指数形式衰减，在4.5～10倍水力直径长度内几乎完全湮灭，流动趋于稳定。棒束稳定流动区内冷却剂的流动特性是限制燃料组件热工水力性能的瓶颈之一，有必要进行深入研究。

本文采用粒子图像测速(PIV)、大涡模拟(LES)以及SSG雷诺时均模拟相结合的方法，对3×3棒束内水的流动特性进行研究。

1 实验装置

1.1 实验回路

图1 实验回路系统Fig.1 Experimental loop system

图1为PIV实验测量系统示意图，主要包括水箱、离心泵、涡轮流量计、实验段以及管路和阀门等。水箱中的去离子水由变频离心泵驱动后分为两个支路，一路作为旁路直接流回水箱，另一路为主路，水流经涡轮流量计后进入实验段，经PIV测量流场后流回水箱。实验段内的工质流量由主路和旁路阀门以及离心泵的变频器配合调节。实验系统的管路、泵和阀门全部为不锈钢材质，水箱由亚克力玻璃制成。

1.2 实验段

图2 实验段结构Fig.2 Structure of test section

棒束实验段结构如图2所示，去离子水由底部进入实验段，在腔室内充分混合后向上流经3×3棒束区，随后在顶部腔室汇集后流出实验段。3×3棒束由直径为9.5 mm的不锈钢棒组成，相邻棒的中心距为12.6 mm，实验段水力当量直径为7.54 mm。中部为长度280 mm的透明区域，其中6根不锈钢棒由透明FEP管代替(如横截面A-A)，外壁由亚克力玻璃加工成型，通过法兰与上游和下游实验段连接。25 ℃时FEP的折射率为1.338，与水的折射率1.33非常接近，通过匹配折射率的方法能减小两种材质交界面附近示踪粒子的位移误差。此外，FEP管外径设计为9.5 mm、壁厚设计为0.5 mm，薄壁可保证测量区处于高透明状态。在实验棒束与顶部和底部腔室以及法兰连接的位置均安装有定位格架，用来支撑和定位不锈钢棒。本实验主要关注棒束稳定流动区的湍流特性，定位格架的结构不同于真实燃料组件中的格架，而是在实现定位作用的同时最大限度地增加流通面积，减小对流场的干扰。

1.3 PIV测量系统

PIV测量系统由激光源、相机、同步器、光学透镜组以及图像处理软件组成。实验采用Dantec Dynamics公司生产的Nd:YAG双脉冲激光器，波长为532 nm，发射频率设置为7.5 Hz，每束激光能量为100 mJ。激光经过光学透镜组后变成厚度约为1 mm的片光源，从侧面照射实验段并穿过两排FEP管中间的区域(图2横截面A-A)。采用平均粒径为10 μm的空心玻璃微珠作为示踪粒子，其密度为1.04 g/cm3，对水具有良好的跟随性。实验所用相机型号为TSI-630059 PowerView 4M Plus，分辨率为2 048×2 048 pixel，配备有Nikon可调焦距镜头。

采用Insight 4G软件进行互相关计算和图像后处理，判读区大小设置为32×32 pixel，重叠窗口为50%，图像分辨率为0.053 mm/pixel。由于粒子位移不能超过判读区域的1/4，根据流速的不同，两束激光发射时间间隔设置为60～180 μs。在每个工况下，共拍摄2 000组瞬时速度图像进行计算和后处理，并采用总体平均后得到平均速度。PIV系统的测量误差主要来源于标定误差、时间间隔误差、粒子拖尾、后处理算法误差等，经计算后，平均速度的合成不确定度为2.75%。

2 数值模拟

2.1 几何模型及边界条件

数值模拟所采用的几何模型与实验段相同，即直径为9.5 mm的不锈钢棒按3×3结构排列，布置在横截面为37.8 mm×37.8 mm的方腔内，模型长度为150 mm。在绝热条件下以水为工质进行计算，温度为20 ℃，压力为101 325 Pa。为加快求解速度，进出口采用周期性边界条件，设置恒定质量流量为0.789 kg/s，对应轴向平均速度为1 m/s。3×3棒束表面和方腔壁面均设置为无滑移边界条件。

2.2 网格划分

LES的思想是使用滤波函数在空间滤波，将湍流分为大尺度涡和小尺度涡，其中大尺度涡直接求解，而小尺度涡则采用亚格子尺度模型进行模化，过滤尺度主要取决于所划分的网格。LES的横截面网格划分如图3所示，近壁面第1个网格距离设置为0.003 mm，径向网格伸展率为1.2，沿流动方向和每根棒周向分别设置151和180个网格节点，总网格数为1 600万。计算过程中最大y+约为0.2，沿展向和流动方向的空间分辨率分别为9.5和76.2，满足大涡模拟的分辨率要求[9]。雷诺时均模拟的横截面网格划分与图3类似，经网格无关性验证后总网格数选取为300万。

图3 横截面网格划分Fig.3 Cross-sectional view of structural mesh

2.3 计算模型与求解设置

采用STAR-CCM+软件进行数值计算，其中LES使用WALE亚格子应力模型，对流项采用二阶中心差分离散，非稳态项采用二阶隐式方法离散。时间步长为5×10-5s，保证最大库朗特数小于0.5。计算物理时间设置为3 s，湍流充分发展后共采集46 944个瞬时速度做时间平均。内循环迭代次数设置为20次，每个时间步长内残差至少降低3个数量级。雷诺时均模拟采用SSG雷诺应力模型进行稳态计算，残差收敛标准为1×10-6。

3 结果分析与讨论

3.1 LES计算结果有效性验证

前已述及，为加快计算速度，LES几何模型的长度缩减为150 mm，同时进出口采用周期性边界条件，通过空间两点速度关联可判断所选模型的长度是否合理。图4为子通道中心处沿流动方向的无量纲速度关联Ri(y)。由图4可知，横向速度u和w沿流动方向的关联较弱，Ru(y)和Rw(y)在距离入口15 mm处即下降到0，但轴向速度v的关联性较强，Rv(y)需要30 mm的长度才将两点速度关联分离。通过对速度关联函数沿流动方向进行积分可获得大尺度涡的积分长度尺度。图4表明，所选取的几何模型长度是足够的，进出口周期性边界条件并未引入附加的周期性效应。

图4 沿流动方向两点无量纲速度关联Fig.4 Normalized two-point velocity correlation along streamwise direction

图5 横向及轴向速度的功率谱密度Fig.5 Power spectrum density of lateral and streamwise velocities

湍流能谱可粗糙地分为3个区域[10]：含能尺度区、惯性子区以及耗散区，LES一般要求将截断尺度设在惯性子区。对LES计算的子通道中心横向速度u和轴向速度v进行快速傅里叶变换，即得到图5所示的功率谱密度分布。由图5可看出，湍流能量主要保持在10 Hz以下的平均流内，湍流耗散开始于200 Hz。此外，惯性子区的斜率与理论上的Kolmogorov -5/3斜率重合，表明当前的网格能求解到惯性子区，即LES计算具有足够的网格分辨率。

3.2 实验与数值模拟结果

1) PIV测量结果

对PIV测得的2 000组瞬时速度矢量进行平均可得轴向平均速度V分布，如图6所示，实验过程中原点设置在测量区域的左上角。在名义速度为1 m/s时，轴向平均速度在0.7～1.4 m/s范围内，中心子通道的流速明显高于棒间隙区，产生这种现象的主要原因是流通截面的不均匀性。中心子通道流通截面积大、流动阻力小，流量再分配后导致该区域流速增大。棒间隙处流通面积小、流动阻力大，轴向速度减小。

图6 轴向平均速度矢量图Fig.6 Vector plot of streamwise mean velocity

实验测量了1、1.5、2 m/s 3个轴向名义速度下的流场，对应横截面雷诺数为7 540、11 310和15 080。图7为横截面轴向平均速度V与名义速度Vm的比值随雷诺数的变化规律。由图7可知，子通道中心速度最高，较棒间隙区高30%左右。此外，不同雷诺数下无量纲轴向平均速度曲线几乎重合，只是在子通道的中心附近，无量纲轴向平均速度随雷诺数的增大而略有减小。由图7还可看出，横截面速度分布并不严格对称，左侧Rod1间隙的流速较右侧Rod3低4%左右。可能的原因是FEP透明管Rod1发生了微小形变，或是该位置激光平面未严格穿过棒与棒的间隙中心，从而导致了一定的测量误差。

图7 无量纲轴向平均速度Fig.7 Normalized streamwise mean velocity

2) 数值模拟与实验结果对比

图8为LES、SSG雷诺时均模拟计算所得的轴向平均速度V和轴向均方根速度Vrms与PIV测量数据的对比，数值计算模型的原点位于中心棒的圆心。由图8a可知，SSG模型对中心子通道的速度预测明显高于实验数据，但对棒间隙(x=0 mm)附近轴向速度的预测较为准确。LES在子通道中心和近壁区较SSG模型有明显的改善，对棒间隙区的预测不理想，但总体上与实验测量值更为接近。由图8b可知，SSG雷诺时均模拟和LES均能定性地预测均方根速度的变化趋势，但在定量上预测偏低。相比而言，LES结果更为准确。

图8 数值计算结果与实验数据对比Fig.8 Comparison of numerical results and experimental data

图9 大涡模拟得到的雷诺应力分布Fig.9 Distribution of Reynolds stresses obtained from LES

LES得到的平均雷诺应力分布如图9所示，在3个雷诺正应力中，轴向分量v′v′较其他两个方向的分量更大，尤其在近壁区和棒间隙区。雷诺正应力u′u′的分布与v′v′相似，但数值上有所降低，而w′w′几乎保持水平。雷诺切应力u′v′沿x方向由负值逐渐转变为正值，呈波浪形变化。值得指出的是，图9所示的雷诺应力分布由LES直接求解的脉动速度得到，并未包含亚格子模型模化的部分。

3.3 大尺度流量脉动现象

图10为PIV测得的轴向瞬时速度云图，虽然经过2 000组瞬时速度平均后轴向速度的分布十分规则(图6、7)，但由图10可知，瞬时速度沿流动方向并未表现出平直分布的特征，而是呈S形变化，这一现象在文献[11]中称为大尺度流量脉动。产生这一现象的主要原因是轴向速度从子通道中心向棒间隙区逐渐减小，导致非常大的速度梯度，这一速度梯度与湍流剪应力共同作用，产生大尺度涡旋。由于大涡并未受到壁面的束缚，其在被主流输运的同时会穿过棒间隙并沿相反的方向旋转。这种大尺度的流量脉动一方面增强了相邻子通道间的湍流交混，强化了冷却剂与包壳的对流换热，另一方面也可能导致流致振动，对燃料棒束造成破坏。

由图10还可看出，大尺度脉动沿轴向的波长随雷诺数的增大而略有增加。在7 540、11 310和15 080 3个雷诺数下，大尺度脉动的半波长分别为23、25、29 mm。图11为LES和SSG模拟得到的速度云图，可看出LES能捕捉类似PIV测量的速度脉动，其脉动的半波长较PIV测量值偏高40%左右，而稳态SSG模拟得到的是横截面对称的速度分布，并不能预测棒束中的瞬态流动现象。

图10 PIV测得的轴向瞬时速度Fig.10 PIV measurement result of instantaneous streamwise velocity

图11 数值计算得到的轴向瞬时速度Fig.11 Numerical result of instantaneous streamwise velocity

3.4 湍流交混

湍流交混促进了子通道间的动量和能量传递，进而强化包壳与冷却剂之间的对流换热，是反应堆燃料组件设计时关注的一个重要热工水力问题。湍流交混系数表征相邻子通道间交混的强弱，其定义为横向有效交混流速与子通道轴向平均速度的比值[12]。本文选取了5个常用的交混系数(β)计算公式(表1)，其预测值与实验测量值的对比如图12所示。由图12可知，β实验值与预测值的变化趋势一致，均随雷诺数的增加而减小。大部分公式都过低地预测了湍流交混系数，实验值与Castellana公式[14]的计算值较吻合，定量上仅较预测值高10%左右。Castellana公式是压水堆燃料组件热工设计时广泛采用的经验公式，具有较高的准确性。实验值高于预测值的可能原因是棒束内存在大尺度的流量脉动现象，如图13子通道间的速度分布所示，主流方向的大尺度涡在被流体输运的同时，会周期性地穿过棒间隙，其横向流动在一定程度上增强了交混效应。

表1 湍流交混系数计算公式Table 1 Correlations of the turbulent mixing coefficient

图12 湍流交混系数实验值与经验公式对比Fig.12 Comparison of measured and predicted turbulent mixing coefficient

图13 相邻子通道间的速度分布Fig.13 Velocity distribution between adjacent subchannels

4 结论

本文采用实验测量和数值模拟相结合的方法，对3×3棒束内水的湍流流动特性进行了研究，得到以下结论：

1) 子通道中心的轴向速度较棒间隙区高30%左右，主要原因是流通截面的不均匀性导致了流量再分配；

2) 与SSG雷诺时均模拟相比，LES所得的轴向平均速度和均方根速度更接近实验数据，对湍流特性的预测较好；

3) PIV测量发现了沿流动方向的大尺度脉动现象，脉动波长随雷诺数的增加略有增大，LES能较好地捕捉这一过程，而稳态SSG雷诺时均模拟无法预测棒束内的流动特征；

4) 湍流交混系数随雷诺数的增加而减小，PIV测量值较压水堆采用的Castellana公式预测值偏高10%左右。