基于OpenFOAM的5×5棒束流动数值模拟

黄瑀琦，蔡杰进

(华南理工大学 电力学院，广东 广州 510640)

反应堆燃料组件是堆芯的主要部件，是堆芯热工水力性能的重要影响因素之一。与研究周期长、费用高的试验研究相比，CFD模拟已成为安全、快速研究手段之一，因此利用CFD进行燃料组件子通道的流动分析研究具有重要意义。燃料组件棒束子通道内湍流流动模拟问题是核热工水力领域长期存在的问题之一，而决定模拟质量的核心便是其内部的算法代码。

MATIS-H基准题是韩国原子能中心主导启动的第二个检验CFD代码的国际基准题试验计划(IBE-2)，其设计目的是检验CFD软件代码在三维子通道单相湍流中的模拟质量。MATIS-H基准题的格架是标准压水堆燃料组件的5×5格架，试验使用了两种搅混翼，分别为分离式和旋流式。

Chang等[1]总结了MATIS-H的试验数据，至2010年12月，共有25款商用和开源的CFD软件代码进行盲测，参与模拟测试的网格数从330万至1.44亿不等，近避面y+值从10-2到150不等。根据测试结果得出一个针对单相流动模拟的最佳预测指导(BPG)[2]，本文中的模拟设置参考了该最佳预测指导。Cinosi等[3]利用商业CFD软件STAR-CCM+计算MATIS-H基准题，指出基于RANS模型求解的速度和压力在峰值的预测上有偏差，低估了格架下游的湍流强度。Mikuž等[4]在2012年新欧洲核能会议上提出OpenFOAM在反应堆热工水力领域，尤其是子通道CFD模型代码研究上具有广泛的应用前景。其团队针对不同工况下棒束子通道，在小尺度几何模型下，同时进行WALE-LES模型和DNS模型的模拟。利用DNS模型数据验证植入的WALE-LES模型，证明OpenFOAM具有模型植入能力[5]。商用CFD软件为便于工程使用以及保护自身知识产权，其绝大多数模型的核心算法代码均被封装。用户并不能阅读和修改模型代码。当商用软件对特定问题进行模拟时(如单相棒束流动)，对部分物理量(如速度、压力)的预测值较为满意，但部分物理量(如湍流强度)预测值与试验值偏差较大，且这种偏差并非来自于网格和收敛精度的影响。Chen等[6]使用商业软件CFX v14.5研究了基于RANS方法的4种流动模型对子通道棒束流动的CFD模拟，在速度和压力上能得到与试验一致性较好的结果，但在湍流强度上却和试验数据相差2～5倍。用户只能通过更换模型才能解决这种问题，但更换模型时需要重新考虑新模型的适用范围和重新设置各种参数，甚至新模型在计算资源占用和耗时上是旧模型的几倍甚至几十倍，因此开源CFD软件在近年来备受研究者的青睐。本文借助开源软件OpenFOAM v8.0，对MATIS-H基准题的几何模型进行全尺寸建模，模拟得到速度分布和压力分布，并研究OpenFOAM对子通道CFD模拟的有效性。

1 5×5燃料组件建模

1.1 几何建模

以MATIS-H基准题为标准建立几何模型，搅混翼为分离式叶片，各围板尺寸均为宽33 mm、长63.40 mm、厚度0.4 mm。定位格架的几何结构复杂，因弹簧、刚突与燃料棒间是线接触或弧面接触，导致局部网格划分难度大，网格质量差，目前的研究大多对定位格架进行简化。本文的简化方法是将弹簧、刚突与燃料棒间改为面接触，并保留定位孔等结构。简化后的格架如图1所示。

图1 简化后的格架Fig.1 Simplified grid

1.2 网格划分

通过OpenFOAM中的blockMesh工具实现外流域的划分，snappyHexMesh工具进行细化网格。利用blockMesh工具生成1个65 mm×65 mm×520 mm长方体外流域。在前处理环节上，OpenFOAM允许用户使用其他商用软件进行网格划分等前处理操作，本文为验证OpenFOAM自身的模拟有效性，利用OpenFOAM具有的snappyHexMesh工具，仅在格架处设置细化区域，并将划分精度等级设置为1级、2级、3级、4级，以便进行网格敏感性分析。其中，3级精度网格划分如图2所示。将不同精度网格所计算出的轴向平均速度与MATIS-H的试验数据进行对比。

网格敏感性分析结果列于表1，可见，4级精度网格数是3级精度网格的2.1倍，但计算精度仅提升了1.33%，计算时长是其3.6倍，3级精度网格已能满足计算要求，误差小于5%。本文将采用3级精度网格进行模拟。

图2 3级精度网格划分Fig.2 Level 3 mesh

表1 网格敏感性分析Table 1 Sensitivity analysis of mesh

1.3 模拟条件设置

入口速度设为1 m/s，密度设为994 kg/m3，与MATIS-H试验设置一致。出口选择压力法向梯度为0，壁面采用无滑移边界处理，平均静压为0 MPa，残差控制收敛精度设置为10-6，雷诺数约为1.10×105，参数设置列于表2。

求解器选用simpleFOAM。simpleFOAM是一种基于SIMPLE算法的求解器，SIMPLE算法由Patankar等[7]于1972年提出，并迅速成为计算不可压缩流体的主流方法。

表2 参数设置Table 2 Parameter setting

2 结果分析

2.1 SST k-ω模型试算

在OpenFOAM下，设置如图3所示4个截面和截面上的通道A作为数据获取源。采用SSTk-ω模型在通道A上进行轴向速度场试算，将试算结果与MATIS-H基准题中的试验数据进行对比。

截面3、4的通道A轴向速度如图4所示，其中V为入口初速度，w为实际速度。截面3的通道A处于搅混翼下游，出现较为明显的速度波动，OpenFOAM中实现的SSTk-ω模型低估了轴向速度波动幅度。截面4位于远离搅混翼的下游，SSTk-ω模型的轴向速度预测效果较好，总体试算结果较为吻合，达到了商业CFD分析软件的应用水平[8]。

2.2 SST k-ω模型模拟结果与讨论

截面1位于格架下方15 mm处，截面2位于格架中部位置，距离格架入口15 mm处，截面3位于格架上方搅混翼下游出口，截面4位于远离格架出口60 mm位置。各截面轴向速度v示于图5。

图3 横截面选取和通道A示意图Fig.3 Selection of plane and channel A

图4 通道A轴向速度对比Fig.4 Comparison of axial speed of channel A

图5 各截面轴向速度Fig.5 Axial velocity of each plane

由图5可见，位于格架中部的截面2，由于单弹簧和双弹簧实体的阻碍作用，弹簧附近的轴向流动速度较小，形成角通道流动速度大、棒束边界流动速度小，受实体挤压导致局部流速达到1.80 m/s。截面3位于搅混翼下游，轴向速度在搅混翼的作用下出现规则对称的扰流，如图6所示。在搅混翼下游出现对称的轴向速度分布，中心通道的流速较角通道流速大，从流线图可看出，在搅混翼下游出现涡流，这种涡流结构之间相互干涉，在搅混翼分布密集的中心通道，涡流之间干涉最为明显。

各截面横向流动速度u示于图7。因位于搅混翼下游，受搅混翼的扰动影响，横向流动速度开始增强，并规则地出现在搅混翼附近，形成1组大小相等、方向相反的横向流动速度，且在搅混翼片两侧的流动速度最大，最大横向流动速度达到入口速度的53%。位于远离格架的截面4上，搅混翼的搅混效果减弱，最大横向流动速度减小到入口速度的20%，但由搅混翼形成的横向流动区域不断扩大，在格架中心通道的横向流场范围最广，这是由于中心通道的搅混翼片数目密度最大所致。

图6 截面3的速度矢量图和流线图Fig.6 Vector and streamlinediagrams of plane 3

图7 各截面横向流动速度Fig.7 Cross flow velocity of each plane

各截面压力分布云图示于图8。在格架中部压力呈现中心通道高、四周通道较低，这是由于格架产生的束流效应使冷却剂向中心流道靠拢，压力升高。位于搅混翼下游的截面3上，经过搅混翼的搅混作用后，横截面上的流体压力分布发生很大变化。在靠近搅混翼处，流体压力较小，在远离搅混翼处，流体压力较大。这是由于搅混翼加强冷却剂扰流作用，使高压区和低压区的分布明显。位于格架下游的横截面4，经过一段规则的流动段，搅混翼对流体的扰动逐渐减弱，高压区和低压区相互扩散，右上角压力与左下角压力呈中心对称，略高于中心通道压力。

图8 各截面压力分布云图Fig.8 Pressure distribution of each plane

3 结论

本文以MATIS-H试验和文献[8-9]中的数值模拟为事实基础，利用OpenFOAM进行数值模拟。在搅混翼的作用下，搅混翼下游截面3上轴向流速最大达到1.28 m/s，相比于入口速度提升了28%，最大横向流动速度达到入口速度的53%，中心通道冷却剂搅混效果明显。经数据对比得出，速度和压力的模拟结果相比试验有10%～20%的波动，属于模型所产生的正常误差范围。压力和温度模拟趋势和物理现象与试验吻合度较高。OpenFOAM具备一定的单相棒束流动数值模拟能力，并能在标准求解器进行控制方程的植入，具有一定的研究前景。

本文仅是初步利用OpenFOAM进行的单相棒束流动模拟，未来还需进行多相流动的棒束流动研究，深入细致流场结构的研究和分析。