堵流工况下棒束子通道流场实验与模拟研究

陈登高，毕景良，黄彦平，袁德文，昝元峰，徐建军

(中国核动力研究设计院 核反应堆热工水力技术重点实验室，四川 成都 610213)

核电作为一种高效的清洁能源受到了能源行业的青睐，并逐渐成为一种重要的能源生产方式。压力容器中的堆芯是热能产生的部位，也是反应堆的关键部件之一，其热工水力特性直接关系到燃料元件的冷却效率及安全问题。在极端情况下，燃料组件子通道可能因为棒束肿胀、异物堆积等发生通道堵流现象，传热恶化，甚至超温、堆芯融化等危险工况，因此有必要开展相关研究。目前，针对通道堵流的研究主要集中在矩形通道的研究，如郭玉川[1]研究了板元件在堵塞情况下的流场分布特性，Xu等[2]及Yuan等[3]研究了板元件加热情况下堵流时流动阻力、传热等特性。此外，采用模拟方法研究子通道堵塞特性也是常用的方法，如李健全等[4]使用系统程序研究了堵塞情况下的反应堆整体热工水力变化；Lee等[5]、Guo等[6]、Salama等[7-9]则使用CFD分析软件计算了堵流情况下矩形通道的流动和传热特性。目前已有大量研究人员使用激光粒子测速(PIV)方法开展了子通道流场测量[12-15]，但仅有很少工作是堵流情况下棒束子通道流场的PIV测量[10-11]，且尚无堵塞情况下子通道流量的测量及分析。堵流下的子通道流量分布是建立系统程序的重要数据，因此有必要开展相关的研究。本文利用PIV方法测量5×5棒束堵塞情况下子通道的流场分布特性及子通道流量分布特性，并基于自研数值分析软件麒麟KILI V1.0和商用CFD软件进行堵流现象的数值模拟，结合实验结果分析棒束子通道在局部通道堵塞情况流场、子通道流动分布等特性。

1 实验装置

本文首先采用可视化实验方法对堵流情况下5×5棒束组件子通道流量分配进行了测量，实验回路如图1所示，主要包括实验本体、水箱、主泵、流量计、阀门和测控仪表等。主泵从水箱中取水加压后一路进入实验本体，一路回流至水箱，两路都设置调节阀以调节进入实验本体的水流量。通过文丘里流量计测量进入本体的水流量，通过压力变送器和热电偶测量水压和温度，所有测量仪表都进行了检定或校准以满足相应的测量精度要求。

图1 棒束组件内流场测量实验装置流程图Fig.1 Schematic diagram of rod bundle flow field measurement facility

图2a所示为实验本体结构示意图，实验本体主体为四面可视的矩形通道结构，进出口有稳流装置，内部放置5×5棒束组件。组件由支撑段和测量段(图2b)组成，支撑段为硬质材料，测量段为设计外径9.5 mm、中心距12.6 mm的薄壁FEP管，以最大程度接近水的折射率。堵流件安装在FEP管外，堵流区域为3×3棒束区域(图2c)，中间通道堵流比例最大，堵流中心通道(图2c中红色部分，C2-3、C2-4、C3-3、C3-4)堵流比例为72%，堵流边通道1(图2c中黄色部分，C1-3、C1-4、C2-2、C3-2、C4-3、C4-4)堵流比例约为20%，堵流通道2(图2c中绿色部分，C1-2、C4-2)堵流比例约为10%。堵流件安装位置距入口大于40D，距出口距离也在20D以上，以最大程度减弱入口和出口效应，测量区域约为20D。图2c中绿色虚线所示位置为拍摄竖直流场所在截面位置示意，具体为棒中心线和通道中心线。

a——本体结构；b——测量段；c——棒束截面结构图2 实验本体及棒束组件结构示意图Fig.2 Schematic diagram of test facility and rod bundle

PIV面激光以xz平面照射测量区域，高速相机则沿y方向拍摄。采用跨帧模式进行图像采集，然后使用专业软件进行流速场计算，即获得了测量区域的流场分布。相机拍摄帧率为15 Hz，每次针对同一工况拍摄200组以上照片，然后求取平均流场。本文还基于PIV测量结果采用数值积分方法获得了棒束子通道流量结果，具体测量和计算方法参考文献[16]。

2 数值模拟方法

数值模拟基于麒麟KILI V1.0进行，该软件采用浸入界面(IST)方法进行网格划分，基于NS方程求解流场、温度场、压力场等。该软件支持常见CAD软件三维几何模型的导入，然后进行IST方法的笛卡尔网格划分，而后进行边界条件设置、模型设置、求解设置等，最后进行迭代求解和后处理。为了保证透光性，实验所用FEP管壁厚仅约0.5 mm，因此易变形导致其定位与设计不同。为了保证数值模拟所建模型最大程度地还原实验棒束，基于PIV所拍摄的棒束图进行图像处理并测量了棒束的真实尺寸，测量原理示意图如图3所示，基于PIV拍摄通道截面图，测量棒束通道参考尺寸LRef和待测量尺寸LB对应的像素数量PRef和PB，然后进行长度计算：

(1)

图3 棒束安装后实际尺寸测量原理示意图Fig.3 Schematic diagram of size measurement of installed rod bundle

图4为采用浸入界面方法时网格划分示意图。图4a为浸入界面方法的基本原理，即在网格划分时不区分流体与固体域，使用流固界面符号距离函数φs进行数值地区分流体与固体域。在求解NS方程时，通过引入Heaviside函数H(φs)自动地区别流体与固体域。该函数在流体内等于1，在固体内等于0，在界面上等于0.5。带堵流件的棒束几何由三维CAD软件根据实验装置尺寸测量结果建立，然后导入麒麟KILI V1.0中进行网格划分，图4b的固体域即是绘制的棒束三维几何。划分网格时直接根据实验装置矩形流道大小绘制矩形流体域，即图4c中的背景网格，这样三维几何就浸没在背景网格中(流体域)。然后基于重叠网格技术在堵流件区域添加加密的子网格对堵流块进行加密，以提高堵流件区域流场解析精度。本次模拟主要针对堵流件上下游约9D的范围，为了避免入口效应，也绘制了堵流件上游10D的网格。堵流件对子通道流场有很强的影响，局部网格精细程度对堵流件下游临近流场模拟精细程度影响很大，因此网格数量对流场模拟结果影响很大，因此本文也探索了不同网格数量对模拟结果的影响。

图4 采用浸入界面方法时网格划分示意图Fig.4 Grid generation of simulated rod bundle using IST method

本次模拟针对实验工况，所涉及的模型为冷态流动，模拟为稳态模拟，主要求解x、y、z方向的速度以及压力。模拟边界条件为流量入口，流量、压力、温度根据实验测量结果设置，模拟中迭代计算到残差小于10-5或不再减小为止。

3 结果及分析

3.1 流速测量结果

本实验中进口流量稳定在约14 t/h，对应的水平截面平均质量流速约1 566 kg/(m2·s)，平均流速约1.6 m/s。实验中测量了堵流件上下游约7D范围内棒束中心竖直截面、子通道竖直中心截面的流场。测量得到的典型竖直截面流速使用截面平均流速进行了归一化处理，即实际流速/截面平均流速，结果如图5、6所示。图5、6所示结果皆为100组PIV测量结果的平均值，对应总时长约6.5 s，横坐标使用棒束间距P进行归一化处理，纵坐标使用棒束外径D进行归一化处理。可看出，在没有堵流件的正常状态下，堵流件所在位置上下游流场基本都以本子通道的竖直流动为主，流场未发生偏斜，不同子通道流速分布相似。在堵流情况下，C1和C2都出现了流场偏斜和流速不均的情况。对照图2c可知C1子通道最大堵流比例约20%，C1-3和C1-4子通道下游约3D范围内流速小于截面平均流速的50%，可见仅堵块下游少量流场受影响。对于C2-3和C2-4，子通道堵塞比例为72%，堵流件下游2D范围内流速整体小于平均流速的20%，超过4D后才恢复到约50%。C2-3和C2-4子通道在堵流件下游也出现了明显的局部回流和漩涡，局部流速小于截面平均流速的10%，堵塞现象十分明显。

a——C1通道；b——C2通道图5 非堵流状态下截面归一化流速测量结果Fig.5 Typical normalized velocity from experiment measuring under non-blocking condition

a——C1通道；b——C2通道；c——R2棒束中心截面；d——R3棒束中心截面图6 堵流工况下归一化流速测量结果Fig.6 Typical normalized velocity from experiment measuring under blocking condition

为了直观地分析堵流对子通道流速的影响，使用流速不均匀指数进行量化。该指数等于指定竖直截面高度z处流速绝对值与截面平均流速的均方根，具体计算公式为：

(2)

图7为根据实验测量结果得到的通道C1、C2、C3、C4在不同高度处的流速不均匀指数，其中虚线所在位置为堵流件，轴向位置0为堵流件轴向中心位置。可看出，通道C1和C4的流速不均匀指数变化较小，仅在虚线所示堵流件区域出现了较小的增加，并很快恢复正常，表明其受堵流件的影响较小，冷却剂可正常流动。而通道C2和C3流速不均匀指数在堵流件处出现了较大的增加，而后缓慢减小；且流速不均匀指数在堵流件上游1D处已开始增加，表明其对流场的影响在上游1D处已开始。虽然在堵流件下游约3D后流速不均匀指数下降较缓，但仍未恢复到堵流前水平。

图7 流速不均匀指数随z的变化Fig.7 Velocity uneven factor at different z values

3.2 子通道流量结果

本研究使用PIV测量获得了子通道流量结果，具体测量和计算方法参见文献[16]，本次流量测量范围为堵流件中心上游约2D至下游约5D范围。图8为获得的堵流件中心下游不同距离处截面竖直方向质量流速的分布，其中质量流速以平均质量流速1 500 kg/(m2·s)进行了归一化处理。从图8可看出，堵流件上游1D处水平截面质量流速分布仍较均匀，而在堵流件下游0.5D处，堵流中心通道(C2-3、C2-4、C3-3、C3-4)质量流速出现了较大降低，上述子通道质量流速普遍在平均值的30%以下。堵流区域临近的非堵流区域则出现了明显的质量流速增加，其最大质量流速超过平均值的2倍，表明堵流导致的子通道质量流速偏移十分明显。质量流速在堵流中心下游4D后恢复较明显，堵流中心质量流速已恢复到平均值的80%左右。图8a中部分棒束周围出现了局部较低的流速测量结果，通过分析原始PIV图片发现，导致这一现象的主要原因是部分棒束在加工制造中出现了难以察觉的刮伤，在激光照射下出现局部亮度过高的情况，进而导致这些局部区域的示踪粒子被掩盖，PIV相机拍摄不到清晰的移动粒子，进而计算为低流速区域。在实验中进行了多次更换新棒束，仍在局部出现了上述现象，但对真实流速分布没有影响，也不影响整体流场结果的判断。图8b中C1-3和C1-4子通道质量流速明显高于C4-3和C4-4子通道的现象，原因是安装后的带堵流件棒束在R1-X侧与流道壁面的间隙(测量值约4.07 mm)大于R5-X侧(测量值约2.17 mm)，导致在堵流截面处流体经C1-3和C1-4侧绕流更多，进而导致质量流速更多。该现象随着远离堵流区域而明显减小，这一差异对堵流件对冷却剂流动的改变趋势判断影响不大。

a——z=-1D；b——z=0.5D；c——z=1D；d——z=4D图8 堵流件下游不同距离处截面质量流速分布Fig.8 Mass velocity distribution of blockage downstream

为了进一步量化比较堵流对子通道流量的影响，本研究基于测量所得子通道流量计算了5×5棒束的4×4子通道截面流量不均匀指数n，其计算公式为：

(3)

图9为堵流件下游不同距离处子通道流量分布系数，其中折线为不同子通道归一化流量，柱状图为子通道截面流量不均匀指数，0所在位置为堵流件高度中心。可看出，在堵流区域上游1D处子通道截面流量不均匀指数明显增加，表明子通道流量分配已受到影响，堵流件下游0.5D范围内流量不均匀指数出现了突增，在1D后不均匀程度有了明显恢复。图9中曲线为堵流中心通道(C2-3、C2-4、C3-3、C3-4)、堵流边通道(C2-2、C3-2)及部分未堵流通道(C2-1、C3-1)的子通道归一化流量。与流量不均匀指数类似，3种类型通道都是在堵流件上游1D距离开始出现了流量变化，其中堵流中心通道流量降低，堵流边通道和未堵流通道则上升。在堵流件下游0.5D范围内，堵流中心通道流量出现陡降，最小流量仅为理论流量的25%左右，而堵流边通道流量则突增，最大流量约为理论流量的1.7倍，未堵流通道最大则约为1.3倍。由图9可知，中心堵流对所在子通道有强烈的限流作用，对邻近的子通道则有强烈的增加流量作用，对邻近未堵流通道有增流作用。3种子通道的流量变化都在堵流下游1D处出现了明显的缓解。

图9 堵流件下游子通道截面流量不均匀指数Fig.9 Sub-channel flow rate uneven index of blockage downstream

3.3 模拟结果及分析

使用前述数值模拟方法进行数值迭代求解可得到5×5棒束在堵流情况下的模拟结果，本研究也使用商业CFD软件(ANSYS FLUENT 15.0)进行了模拟，网格绘制时尽量保证二者网格总数和最小网格尺寸相同，此外二者湍流模型(标准k-ε模型)、边界条件、求解算法等相同，以对比麒麟KILI V1.0与商业CFD软件。图10为模拟所得C1、C2子通道中心竖直截面流场的典型模拟结果，其中流速为使用截面平均流速归一化处理的流速，其中白色框部分为堵流件所在位置。C2子通道位于堵流中心通道，堵塞比例较高，模拟结果中堵流件下游形成了低速区，这与实验测量结果吻合。比较麒麟KILI V1.0和商用CFD软件的模拟结果可发现二者流场分布十分接近，表明麒麟KILI V1.0可很好地模拟流场分布。

a——C1，KILI V1.0；b——C2，KILI V1.0；c——C1，CFD软件；d——C2，CFD软件图10 模拟所得堵流情况下归一化流速云图Fig.10 Typical normalized velocity contour from simulation under blocking condition

为了更加直观地对比麒麟KILI V1.0和商用CFD软件差别，提取了堵流件下游不同距离处实验与模拟的结果，如图11所示(C2通道中心截面)。可看出，在两种距离处模拟结果可较好地反映实验结果，其中1D处堵流中心通道流速较低，实验测量中有强烈回流；麒麟KILI V1.0与商用CFD软件在细节流速模拟上还原实验测量值准确程度相当。

a——堵流件下游1D；b——堵流件下游4D图11 不同距离处流速的实验与模拟结果对比Fig.11 Comparison of experimental and simulated velocities at different locations

图12为模拟结果的均方根误差(RMSE)，其计算公式为：

(4)

式中：vi,Exp为流速测量值；vi,Simu为对应位置流速的模拟值。

a——z=1D；b——z=4D图12 模拟值均方根误差Fig.12 RMSE for different simulation methods

RMSE越小，表明模拟值与实验测量值偏差越小，模拟不确定度越小。麒麟KILI V1.0和商用CFD软件尽量保证网格总数和最小网格尺寸相同。可看出，整体上1D处的模拟值不确定度大于4D处。1D处回流较强，对湍流模拟准确度要求较高，因此采用RANS模型难免对局部细节模拟准确度有限。模拟不确定度整体随网格总数增加而降低，由于网格总数增加主要是对堵流件附近区域网格进行加密，因此1D处不确定度下降较快。在1D处RMSE随着网格总数由30万逐渐增加到220万不断降低，在4D处网格总数达到140万后RMSE下降很少，表明在1D处需要更精细的局部网格才能实现网格无关，在4D处因为湍流强度降低，网格大小降低到一定程度后即实现了网格无关。总体而言，麒麟KILI V1.0在1D处与商业数值分析软件不确定度相当，在4D处部分网格总数下麒麟KILI V1.0优于商业软件。网格总数较小的20万和40万，商业CFD软件模拟准确度绝对值上更低，主要原因是商业CFD软件采用贴体网格和局部边界层加密，更有利于湍流模拟。麒麟KILI V1.0采用IST技术的网格，网格划分不区分流固区域，网格不贴体，且采用局部整体网格加密，因此其湍流模拟准确度会降低。在网格加密到一定程度后，网格尺寸对模拟准确度影响较小，所以麒麟KILI V1.0在网格总数增加后与商业CFD软件差别减小，甚至更优。

图13为堵流件所在区域流场特征模拟结果。流线图表明堵流件下游区域产生明显的回流区，导致流速突降，这可能导致棒束冷却不足；矢量图用x方向流速值进行标志颜色，红色表示正向，蓝色表示负向，可看出，堵流件的阻塞导致所在子通道冷却剂向周围子通道流动，这是堵流中心子通道流量减少，堵流边通道流量增加的主要原因。

4 结论

本文基于棒束子通道流场可视化实验装置和5×5棒束子通道堵流实验本体获得了子通道堵流72%比例下的流场分布及子通道流量结果，并使用麒麟KILI V1.0对堵流情况下的子通道进行了数值模拟分析。结果表明：1) 在子通道堵流约72%的情况下，堵塞子通道上游约1D范围流场开始受到明显影响，堵流通道下游0.5D距离内堵流影响最强烈，堵塞中心子通道流速仅为平均流速的约30%，子通道流量仅约为理论流量的25%，堵流影响在1D后出现明显缓解，堵塞中心子通道流量恢复到理论流量的60%以上；2) 基于浸入界面方法的麒麟KILI V1.0能对堵流情况下的棒束子通道流场进行高准度的模拟，其模拟结果与商用CFD软件相当。模拟结果表明，堵流件下游的局部回流是导致流速降低的主要原因，堵流通道阻力增大，使得冷却剂向堵塞通道临近通道分流，导致堵塞通道流量降低，堵流中心通道相邻子通道流量增大。