基于OpenFOAM的铅基快堆流体通道内两相流研究

杨 晟，张有鹏

(复旦大学 现代物理研究所，上海 200433)

能源是人类社会赖以生存的根本，核能作为能源结构的重要组成部分，对人类的发展至关重要。第4代核能系统是目前先进核能技术的代表，相比于传统的核能利用设备，在燃料增殖、系统安全性等方面有了长足的提升和改进。铅基快堆是具有代表性的6种4代堆中发展潜力巨大的一种，具备小型化的潜质，在氢气制备、发电、燃料后处理等方面有巨大优势，适用于更多分布式小型发电的场景，目前各国都对其投入了巨大的研发力度。

在铅基快堆研究中，由于重金属的热工水力特性与传统介质差异巨大，铅/铅铋的导热性和热膨胀性好，且铅的辐射屏蔽性造成实验上的研究非常困难，不能采用X光等装置研究内部流场的情况。国内关于铅基快堆的热工水力研究大部分采用数值模拟，且大部分研究集中于堆内单相流，而关于堆内两相流的研究较少。当铅基快堆发生蒸汽发生器管道破裂事故(SGTR)时，二回路内的高压水会进入常压的一回路，从而在堆芯形成两相流动[1]，因此对堆内两相流的模拟十分必要。2017年Ge等[2]采用k-ε模型研究了单相条件下通道内的热工水力现象。Nishi等[3]和Suzuki等[4]对铅铋合金两相泡状流进行了研究，得到在低含气率下堆内流体的宏观流动特性。国内西安交通大学[5-6]和上海交通大学[7]对铅堆内的三角形通道分别用子通道程序和OpenFOAM程序模拟了单相流体的热工水力特性。中国核动力研究设计院[8]对以水为介质的两相流进行了CFD仿真，为后续以铅铋合金为介质的两相流仿真提供了基础。

开源的CFD软件OpenFOAM可通过自行编写求解器的办法实现流体力学方程的求解，因此本文拟采用该软件，应用基于VOF方法的数值模拟，构建铅基快堆中常见的三角形内通道模型，从而研究铅基快堆通道内的气液两相流热工水力特性。

1 物理模型建构

1.1 两相流模型

计算流体力学是通过利用质量、动量和能量守恒的控制方程组，在边界条件给定的情况下，使用现代高性能计算机模拟温度、压力、速度等相关物理量的方法。OpenFOAM是一款基于C++的开源CFD软件，是基于开源操作系统Linux设计开发的。相较于目前已成熟化的商用Fluent、CFX等软件，OpenFOAM的开源性使其求解控制方程时能更加灵活地根据用户需求自行编写求解器，让操作者有更加灵活的自定义功能。同时，它基于C++的编译环境，对于类的应用，如继承和派生方面，有了更多的交流性，能发挥较商用软件更好的交互性。在两相流模拟中，VOF方法是一种重要的数值模拟手段，它通过设定每个网格内的两相流含量，从而勾勒出两相流的界面。

质量方程可描述为：

(1)

动量方程可描述为：

(2)

能量方程可描述为：

(3)

式中：t为时间；U为流体速度；ρ为气液两相流体的混合密度；τ为流体表面黏性应力；g为重力加速度；T为流体温度；σ为表面张力系数；α为两相流中液相的占比；k为流体传热系数；κ为表面曲率；cv为比定容热容。

由于方程结构简单，可编写和可操作性强，VOF模型广泛应用于多相流的模拟仿真。

1.2 通道参数设置与网格划分

铅/铅铋合金具有沸点高、熔点低、导热性好，以及优越的中子经济性，是先进反应堆冷却剂的优秀材料。对于常规设计的液态铅/铅铋合金冷却剂先进反应堆，燃料棒的排列方式通常按照三角形方式，构成六边形组件，如图1所示，内部可分成冷却剂流动的内通道、边通道和角通道[9]。本文选取的研究通道尺寸如下：燃料棒直径6.55 mm、燃料棒间中心距9.17 mm、冷却剂进口流速0.21 m/s、进口温度583 K、轴向长度750 mm。网格划分是对所需计算的流体单元进行空间区域离散化的过程，根据通道参数划分网格，由于通道的几何参数具有对称性，因此在几何模型建构中采用对称性边界条件，网格在z方向的截面如图2所示，运用这种方法，大幅缩减了通道内的网格数量，提升了计算效率。最终的网格如图3所示，网格质量检查结果表明网格满足模拟条件。

图1 燃料组件子通道Fig.1 Fuel assembly subchannel

图3 网格划分Fig.3 Grid division

1.3 组件的网格和边界条件

由于铅铋合金具有更强的输热能力，可容许堆芯在更高的功率密度下运行，因此在相同功率水平下，铅铋合金冷却的反应堆堆芯较传统轻水反应堆更紧凑。在组件的物理模型建构中，按照目前的ADS设计参数[10]，燃料棒的直径为6.55 mm、燃料棒间的中心距为9.17 mm。为使计算结果收敛性更高，在网格划分中采用结构化网格，同时在燃料棒的壁面和六边形组件的外壁面设置边界层。无绕丝情况下组件z截面上的网格剖面示于图4。通过OpenFOAM中的“checkMesh”网格质量检查方法进行检查，结果表明，网格条件符合计算要求。

图4 组件的网格划分Fig.4 Grid division of assembly

图5 计算边界条件Fig.5 Calculation boundary condition

根据中国科学院[11]和中国科技大学[12]的选取方法确定边界条件，所图5所示。设置出口处为自由流出，压力为0 Pa；进口处为速度进口条件，冷却剂的进口流速为0.21 m/s，进口温度为583 K；外壁面为绝热无滑移壁面条件；燃料棒表面为均匀热流密度的无滑移壁面条件。

1.4 铅铋合金物性参数及湍流模型

在计算过程中，采用铅铋合金的热物性参数，由于铅铋的特殊热物性，在模拟中会有与常规冷却剂不同的现象，根据OECD/NEA[13]的推荐，选取如下。

密度为：

ρLBE=11 096-1.323 6T

(4)

动力黏度为：

μLBE=0.004 56-7.03×10-6T+

3.61×10-9T2

(5)

比定压热容为：

cp,LBE=159.0-0.027 2T+7.12×10-6T2

(6)

热导率为：

λLBE=3.61+1.517×10-2T-

1.741×10-6T2

(7)

在两相流的气相流动中，选取的气相为空气，使用理想气体模型的热物性参数进行模拟。在湍流模型的选择中，根据前人[14]的研究结果选择k-ε模型。该模型具有计算简单、在实验室研究和工程实践中应用广泛等优点，适于管内流动的数值模拟。根据OpenFOAM的推荐公式，k-ε模型可通过以下3个方程进行描述：

(8)

(9)

l=0.007D

(10)

采用以上步骤即可完成求解器和求解环境的构建，物理模型的网格划分，湍流模型、初始边界条件和铅铋热物性参数的设定。完善物理量的求解格式后，运行OpenFOAM求解器，即可求解出相应的物理量。

2 实例验证

2.1 内通道验证

为验证OpenFOAM软件在铅基快堆模拟中的准确性，采用进口为铅铋流体在堆芯内均匀加热的方式进行模拟验证，以相同条件下西安交通大学实验校核的SACOS-PB子通道程序的分析结果[6]进行对比，结果示于图6。从图6可看出，冷却剂温度随轴向距离的增加逐渐升高。两程序在相同的进出口条件下，冷却剂温度的上升趋势一致，在各轴向截面内的数据吻合良好[6,15]，因此可认为OpenFOAM模拟的结果准确、实验方案可行，可用于铅基快堆通道的两相流模拟分析。

图6 冷却剂温度的轴向变化Fig.6 Axial variation of coolant temperature

2.2 组件单相流验证

设置进口边界条件为流速不变的纯铅铋流动，即将VOF模型中的气相流动速度设置为0 m/s，模拟的出口温度分布示于图7。从图7可看出，在出口处，温度呈现中间高、周围低的分布，与文献[11]相吻合，从而验证了OpenFOAM对燃料组件模拟的准确性。同时从图7还可看出，铅铋在组件内流动时，组件盒的六边形角通道的温度最低，这是由于此区域的受热面积最小，出口冷却剂的温升也小。

图7 单相流出口温度分布Fig.7 Outlet temperature profile of single phase flow

3 不同条件下OpenFOAM的流动传热模拟

3.1 进口流速的影响

不同流速下，两相流内气体聚集现象的激烈程度也会不同，从而在两种介质中会存在不同的动量和能量交换过程。设置进口流速为0.21、0.25、0.30、0.35 m/s，并控制入口气体质量分数(入口含气率)为10%，考察通道出口处温度的变化，结果示于图8。从图8可看出，随着流速的增大，出口处温度明显下降。在建模过程中，确定的壁面边界条件为等功率边界条件，因此随着进口流速的增大，出口冷却剂的温度降低。虽然流速的增大促进了流体间动量和温度的交混，但两相流中气体的气泡在流动过程中位置更加随机，使传热面更不稳定，形态随时间不断变化。同时，新进入的冷流体与壁面接触的时间缩短也使得出口处流体温度下降。

图8 不同进口流速下的出口温度Fig.8 Outlet temperature under different inlet flow rates

3.2 入口含气率的影响

气体以不同入口含气率进入通道时，气体之间会发生聚集等现象，同时会与液体发生动量和能量交换，使得铅铋冷却剂在通道内的流动与单相流时不同。

图9 不同入口含气率下出口处气泡形态Fig.9 Outlet gas form of different inlet gas contents

设置入口含气率分别为1%、5%、10%、20%，固定进口流速为0.21 m/s，采用均匀进气的方式进行OpenFOAM模拟，所得通道出口处的气泡形态如图9所示。

由图9可看出，随入口含气率的增大，出口处气泡的体积增大，数量增多。在入口含气率不断增大的过程中，出口处会出现更多的气泡团块，主要集中在通道的中心部分。从不同入口含气率的出口气体形态可看出，出口气体基本在通道内部流动，粘附壁面的少。

3.3 空泡的轴向分布

设置入口含气率为10%，在轴向高度分别为0.15、0.35、0.65 m处截取3个z方向截面，截面的气泡形态如图10所示。图10表明，在入口为均匀混合两相流时，随着冷却剂的流动，在轴向上，两相流开始出现气泡的聚合现象。

从图10可看出，冷却剂在通道内流动的过程中，两相流在轴向的均匀性逐渐降低，出现单相聚集现象，随着轴向高度的增加，在通道的中心处会出现气泡聚合。轴向高度为0.15 m时，气体和铅铋冷却剂的边界还不明显，只是出现了VOF模型中气体体积占比较低的区域粘连，出现模糊的气液分隔边界。轴向高度为0.35 m时，通道的中心区域开始出现气体聚集区，图中蓝色区域即为小气泡的形态。当轴向高度为0.65 m，即接近出口处时，气泡已形成近似圆形的形态完整的气泡，其边界清晰。

3.4 网格含气率和速度与轴向高度的关系

气液两相流中气体的存在会使得流动形态更加紊乱，同时气相和液相不同的流动速度也会导致通道内混合流体的流动速度发生变化。

入口含气率为20%时，网格内的含气率和两相流体的速度与轴向长度的关系示于图11。从图11可看出，轴向长度增加时，网格内的含气率变得不均匀，从原先的20%变得时而增大，时而减小，这是气泡在轴向高度上聚集和破裂造成的。同时可看出，在含气率高的网格内，网格单元内的两相流体速度也加快，说明在两相流流动过程中，含气量大的网格两相流动速度也大。气体在通道内的流动速度较铅铋高，使得当网格内两相流的含气率上升时，两相流体速度明显加快。

图11 网格内两相流含气率和速度与轴向长度的关系Fig.11 Relationship between gas content and velocity with axial length

4 不同入口含气率下燃料组件内的流动模拟

图12 不同入口含气率下组件出口处的气体形态Fig.12 Assembly outlet gas form under different inlet gas contents

在不同入口含气率下，燃料组件内的流动和传热会发生改变，同时气泡间的聚合和分离使得组件内的流动情况十分复杂。为进一步研究燃料组件内边通道、角通道和内通道中气泡的分布，采用组件模型，设置入口流速为0.21 m/s，入口含气率分别为1%、5%、10%和20%进行模拟，结果示于图12。从图12可看出，随着入口含气率的增加，出口处的浅色斑块逐渐增多，表明出口气泡在增多，出口处含气量也在增加。入口含气率为1%时，出口的气泡形态并不明显，没有形成非常大的气泡，由于VOF模型描述的是气体与液体在混合状态下的百分比，所以在入口含气率不高时，液体与气体的交界面并不明显。当入口含气率上升后，图中的白点逐渐增多，即气泡的形态开始明显。综合图中组件出口气泡形态可知，气泡在六边形组件的角通道处易聚集，由于气体的导热性差以及比热容小，将导致在角通道处的温度发生畸变，可能出现燃料棒的温度骤升，造成烧毁。因此，如果模拟计算中发现角通道是组件中气体的聚合处，需对其深入研究。

5 结论

本文运用基于开源的CFD软件OpenFOAM及VOF模型研究了铅基快堆中铅铋冷却剂在燃料组件内的两相流动，通过划分模型网格，在OpenFOAM中植入铅铋冷却剂的物性参数，采用k-ε湍流模型，得到了铅铋冷却剂的出口条件，并采用子通道程序和单相流实验对其进行了验证。通过改变入口边界条件的办法，研究了内通道的两相流动，并分析了轴向高度上气泡的聚集情况。最后研究了燃料组件盒的两相流动特性。得到以下结论。

1) 冷却剂出口温度会随流速的增大而降低，而气泡的出口形态基本保持不变，但流体与加热壁面的交界面的高温区域变薄，这是气泡与液体在动量和能量交混过程中，在出口处流动速度增大、加热时间缩短造成的。

2) 随着入口含气率的增大，出口处气泡会增多，且气泡出现大的团块聚集。气液两相流在内通道的流动过程中，气相基本在通道内部流动，不与加热壁面粘附。

3) 冷却剂在内通道的流动过程中，两相流在轴向的分布均匀性逐渐降低，出现单相聚集的现象，随着轴向高度的增加，在通道的中心处会出现气泡聚合。

4) 在燃料组件内，组件出口处气泡在六边形组件的角通道处易聚集，由于气体的导热性差以及比热容小，将导致在角通道处的温度发生畸变，可能出现燃料棒的温度骤升，造成烧毁。