1. 引言

铅铋池式快堆是一种极具发展潜力的第四代核能系统，在燃料增殖和嬗变方面具有独特优势，具有良好的非能动安全特性和经济性，且有利于实现小型化，是目前国际核能领域研究的热点。池式快堆一回路关键设备浸泡在一个大型铅铋池中，这些设备在与容器外部相连接的过程中，穿过堆池的自由液面及其上的氩气空间，由于液体表面的波动，被覆盖的氩气可能以夹带气泡的形式混入液态铅铋中。这种携带的氩气会导致以下运行问题：a) 大量气体流经主循环泵会导致气蚀等恶劣流动现象出现，影响主循环泵的工作性能，降低循环泵的运行寿命。b) 热交换效率降低，进而影响堆芯冷却。c) 被夹带气体通过主循环泵吸入堆芯将严重威胁反应堆安全，其带来的正反应性可能引起反应堆功率震荡，更严重的会导致堆芯传热恶化、烧毁。

国内外学者通过实验和数值计算的方法针对各种类型的气体夹带进行了研究，相应为本研究中自由液面波动致气体夹带的影响因素和表征方法提供参考。K. Satpathy [1] 提出可能存在的四种类型的气体夹带机制：涡旋诱导夹带、液体下落诱导夹带、排水型涡旋形成的夹带和气液界面剪切导致的夹带。Eguchi [2] 等进行了小规模水力实验，发现流体黏度的增加可以抑制夹带，只有当流体黏度的弗劳德数大于临界值时才会发生夹带。巫英伟 [3] 等针对池式钠冷快堆中的液体下落型气体夹带展开了数值研究，分析了几何尺寸、初始液位高度和入口速度对气体夹带的影响。

总的来看，目前关于池式快堆气体夹带的研究多针对钠冷快堆展开，而本研究的对象为池式铅铋快堆，与以往研究的冷却剂物性不同所导致控制夹带过程的重要参数的变化也是待解决的关键问题之一。本研究基于OpenFOAM平台利用流体体积界面跟踪方法对液体下落型气体夹带特性展开研究，同时分析几何结构参数、热工水力参数对池中气液两相相互作用的影响规律。本研究阐明了四代铅铋池式快堆发生气体夹带的过程与机理，为池式快堆的安全运行提供一定的建议和指导。

2. 理论基础模型及验证

2.1. VOF模型控制方程

VOF模型适用于模拟相互之间有清晰边界的情况，能够准确地跟踪和描述这些相互作用的表面，并允许在相交界面上进行质量、动量和能量的交换。VOF模型的基本思想是通过跟踪每个单元中相的体积分数来表示不同相的存在 [4] ，在流场中定义了一个变量 $\alpha$ 来表示液相的占比。对于某一个网格单元的气液两相系统，如果此网格单元内充满了气体，则 $\alpha =0$ ；如果此网格单元内充满了流体，则 $\alpha =1$ ；如果 $\alpha$ 的值介于0和1之间，则此网格单元内为气液混合状态。

气液两相混合物密度 $\rho$ 和粘度 $\mu$ 可以由液相的体积分数计算得出，如式(1)和式(2)：

$\rho =\alpha {\rho }_l+\left(1-\alpha \right){\rho }_v$ (1)

$\mu =\alpha {\mu }_l+\left(1-\alpha \right){\mu }_v$ (2)

下标 $l$ 和 $v$ 分别代表液相和气相。

通过求解体积分数 $\alpha$ 的输运方程来模拟流场中不同相的运动，如式(3)，该方程描述了相对速度、质量传输以及相的变化，其求解允许在空间和时间上追踪相界面的演变。

$\frac{\partial \alpha }{\partial t}+\nabla \cdot \left(\alpha U\right)=0$ (3)

其中， $U$ 为混合物速度。

式(4)和式(5)为使用VOF模型计算的守恒方程，与单相流不同的是，VOF模型在动量方程中附加了随空间变化的粘度，以及附加表面张力项。

质量方程：

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho U\right)=0$ (4)

动量方程：

$\frac{\partial \rho U}{\partial t}+\nabla \cdot \left(\rho UU\right)-\nabla \cdot \tau =-\nabla P_{rgh}-g\cdot h\nabla \rho +\sigma \kappa \nabla \alpha$ (5)

式中： $t$ 为时间； $U$ 为混合流体的速度； $\rho$ 为气液两相混合物密度； $\tau$ 为流体表面黏性应力； $g$ 为重力加速度； $T$ 为流体温度； $\sigma$ 为表面张力系数； $\alpha$ 为两相流中液相的占比； $\kappa$ 为表面曲率 [5] 。

在OpenFOAM中利用VOF模型求解两相流的求解器为interFoam，interFoam本身不具有温度场求解功能，这里植入温度输运方程，如式(6)。

能量方程：

$\frac{\partial \rho T}{\partial t}+\nabla \cdot \left(\rho UT\right)-\nabla \cdot \left(D_{eff}\nabla T\right)=0$ (6)

$D_{eff}$ 是混合相的导热系数：

$D_{eff}=\frac{\alpha k_l}{c_{vl}}+\frac{\left(1-\alpha \right)k_v}{c_{vv}}$ (7)

$k$ 为流体传热系数， $c_v$ 为比定容热容，下标 $l$ 和 $v$ 分别代表液相和气相。

2.2. VOF模型验证

溃坝问题在计算流体力学中由于其自由表面变化多样，常被用来验证数值模拟的可靠性。因此，利用溃坝问题来验证VOF模型在气体夹带研究中的适用性。如图1所示，该问题的实验装置是一个长度为4 m、高度为3 m的开放式水箱。实验前，左侧的水保持在长度为1 m、高度为2 m的位置范围内，其余部分为空气，模拟水柱坍塌晃动直至平稳的状态。

左墙水位高度随时间变化的曲线如图2所示。结果表明，仿真结果与以往的实验和仿真结果吻合较好 [1] [6] [7] ，验证了在这类问题中采用VOF模型的合理性。

3. 液体下落型气体夹带模拟

3.1. 计算模型及条件设置

使用VOF模型对气体夹带进行瞬态模拟，对于三维模型，需要很高的计算能力以及较大的计算资源。为了简化计算，针对热池局部区域建立二维模型，由于反应堆的对称性，只需建立整个模型的二分之一。二维模型及其简化示意图如图3所示，对于液体下落型气体夹带，由控制插头(CP)和主换热器(IHX)的两个轴组成的纵向截面被视为特征平面，如图3(a)所示。在X中红色虚线所包含的区域是要关注的主要区域，经过简化后如图3(b)所示，堆芯与控制插头之间的空隙冷却剂具备向热池流动的趋势，可视作局部区域的入口，冷却剂流入主换热器的入口可视作从出口流出局部区域。

如图3(c)所示，在局部的计算区域中，液态铅铋从左侧壁面底部的入口以 $V_{in}$ 的速度流入，从右侧壁面距离底部高度为 $H_O$ 的出口流出。整个过程模拟了冷却剂通过堆芯和控制插头之间的轴向间隙进入，进而流入热池，最终从主换热器(IHX)的入口离开热池的整个过程。局部区域入口和出口的高度分别表示为 $h_i$ 和 $h_o$ 。W为局部区域宽度，H为局部区域高度，H_l为液面初始高度。模拟工质为液态铅铋以及覆盖氩气，铅铋物性根据OECD出版的铅铋物性计算手册计算，工作温度根据工程设计经验建议选取 [8] [9] ，工质主要物理性质如表1所示。

湍流模型选用雷诺时均(RANS)模型，具体是由Launder与Spalding提出的标准k-epsilon模型，适用于高雷诺数流动计算。局部区域周围壁面设置为无滑移条件，入口和出口分别被施加有恒定速度边界条件和恒定压力边界条件。采用PIMPLE算法求解压力–速度耦合方程，压力和温度方程的残差设定为 ${10}^{-7}$ ，速度方程的残差设定为 ${10}^{-6}$ ，时间步长设定为0.001 s，计算过程中自动调节。

3.2. 网格无关性分析

针对简化后的局部区域模型，利用ANYSY ICEM建立非结构化四边形网格。对于表2中的工况，网格尺寸取0.2~0.5 mm进行网格独立性分析，计算结果与网格数之间的相关性如图4所示，当网格尺寸为0.29 mm时，自由液面平均速度和液面隆起高度分别为0.113 m/s、0.023 m；网格尺寸为0.2 mm时，自由液面平均速度和液面隆起高度分别为0.114 m/s、0.023 m，自由液面平均速度相比较变化仅为0.877%，因此，可以认为数值计算结果与网格大小无关，从而很好地验证了网格的无关性。考虑到节省计算资源和计算成本，选择0.29 mm的网格尺寸进行数值计算，对应计算网格数量为129,688个，同时夹带产生的气泡直径通常在1毫米到几毫米的范围内，因此选择0.29 mm的网格尺寸是合理的。当网格尺寸为0.29 mm时，所有网格质量均高于0.95。请注意，这里为使表征参数特征明显，选取的工况设定较大的入口速度。

4. 液体下落型气体夹带特性与参数影响研究

4.1. 气体夹带评价模型研究

液体下落型气体夹带是由液体流撞击自由表面并在液面上方形成局部隆起并再次回落引起。如图5所示，液体下落型气体夹带发展的整个过程按照时间推进可分为初始阶段，液面波动阶段，夹带阶段和结束阶段。

初始阶段：该阶段是流体静止时的初始状态。

液面波动阶段：如图5(b)所示，从入口流入区域的流体导致左侧液面上升，靠近出口的部分流体流出导致右侧液面略微下降。如图5(c)所示，进入的流体首先沿底部向右侧流动，然后在右壁的冲击下逐渐变为沿右壁向上流动，导致右侧液面隆起。自由表面的速度也随着动能转化为势能而逐渐降低。如图5(d)所示，随着速度的持续增加，自由表面涌浪从右上角向左下角流动，由于流入率和流出率之间不平衡，左侧液体积累导致左侧液面膨胀。

夹带阶段：图5(e)为夹带开始的阶段，气泡进入铅铋的点被称为再淹没点，这时液体从高处坠落，再淹没点附近的自由表面剧烈扭曲，卷吸夹带覆盖氩气进入液态铅铋中，通常，气体夹带是指至少形成一个气泡，气泡四周都是液体。图5(f)为气泡持续夹带进入铅铋的过程阶段。图5(g)为夹带完成的阶段，此时存在大量随局部流场向下运动的气泡，自由表面不再发生气泡夹带卷吸。

结束阶段：随着时间推进，流出和流入之间达到了动态平衡，液面波动趋于平缓，如图5(h)所示。同时在某些情况下，特别是在入口速度较低的情况下，不会发生气体夹带。

计算结果可表示为气液两相份额分布图，定性判断发生气体夹带的方法为观察到液相中至少有一个气泡从两相界面夹带卷吸进入到液相中，并满足气泡四周都是液体的情况，否则不发生夹带。由于定性判断方法易导致偏差，因此采用定量判断方法。如式(8)和式(9)所示，根据夹带进气相所占面积 $A_g$ 以及夹带进气相份额平均值 ${\alpha }_{average}$ ，定量判断是否发生气体夹带现象。

$A_g={\displaystyle \underset{i=1}{\overset{n}{\sum }}{A}_i{\alpha }_i}$ (8)

${\alpha }_{average}=\frac{{\displaystyle {\sum }_{i=1}^n{\alpha }_i}}{n}$ (9)

其中， $\alpha$ 为气相份额，i为单元网格序号，n为网格数量。

利用再淹没速度 $V_{rs}$ 、再淹没角 ${\theta }_{rs}$ 以及液面隆起高度 $h_s$ 三个参数表征自由表面气体夹带特性。如图5(e)所示，再淹没速度 $V_{rs}$ 为该时间刻最大的自由表面速度。再淹没角 ${\theta }_{rs}$ 为再淹没速度矢量方向与水平线的夹角。液面隆起高度 $h_s$ 为自由表面发生形变的最高点与再淹没速度所在位置点的差值。

值得注意的是，在发生夹带的情况下，数据是在夹带即将发生之前读取的，不发生夹带的情况下，数据在自由表面发生最大形变时读取。为提高数据测量精度并节省用于数据监测的资源，气液两相交界处自由表面在后处理软件ParaView中的数据监测标准为 $0.001<\alpha <0.999$ ， $\alpha$ 为气相份额。

4.2. 气体夹带参数影响分析

改变局部区域的入口高度 $h_i$ 、水平宽度W以及液面初始高度 $H_l$ ，设置不同工况条件进行数值计算，如表3所示，并对各参数的影响进行了分析。

对比发生气体夹带和不发生夹带流体的发展情况，图6为数据读取时刻下局部区域速度分布图，图7为自由表面数据读取点的速度随时间的变化，均以工况1a中入口速度 $V_{in}$ 为0.33 m/s和0.37 m/s为例。根据表3的工况进行数值计算，能够得到规律，入口速度较小时，不会引起气体夹带的发生，随着入口速度的增加，气体夹带的趋势增强。如图6和图7所示，入口速度增加，导致流场整体速度增加，自由液面的波动更加剧烈，从而夹带进入覆盖气的概率越大。

纵向比较1a~1c工况，其他几何条件均相同，仅改变液面的初始高度，可发现，随着液面的初始高度增加，发生气体夹带的最小速度增大，可见液面初始高度的增加阻碍气体夹带的趋势。如图8所示，逆时针流场顶端与自由液面的距离增加，自由液面波动趋势减弱，气液两相交界面流体速度方向与水平面的夹角趋于平缓，比较3a~3c工况同样可以印证规律。值得注意的是2a~2c工况与4a~4c工况，其中2a与4a工况，具有较大局部区域水平宽度(W)和较高的入口湍流强度(即较小入口速度V_in或较大入口高度h_i)，这种情况由于液体在两侧多次相遇而引起自由表面波动，与常规情况不同，称为偏差情况，因此水平宽度不能作为单一影响因素分析。

横向逐一比较1a~1c工况和3a~3c工况，其他几何条件均相同，仅改变入口高度。1a工况下发生气体夹带的最小速度为0.35 m/s，3a工况下发生气体夹带的最小速度为0.45 m/s；1b工况下发生气体夹带的最小速度为0.39 m/s，3b工况下发生气体夹带的最小速度为0.47 m/s；1c工况下发生气体夹带的最小速度为0.45 m/s，3c工况下发生气体夹带的最小速度为0.49 m/s，可见随着入口高度 $h_i$ 的增加，发生气体夹带的最小速度增大，发生夹带的趋势就越弱。如图9所示，入口高度越大，从入口进入区域的流体更多从出口直接流出，逆时针流场区域变小，导致自由液面波动趋势减弱，横向比较2a~2c工况和4a~4c工况亦然。

5. 结论

本文基于OpenFOAM使用VOF方法对铅铋池式快堆局部区域液体下落型气体夹带问题进行了数值计算，研究了几何结构参数、热工水力参数对气体夹带特性的影响，得到如下结论。

1) 液体下落型气体夹带发展的整个过程按照时间推进可分为初始阶段，液面波动阶段，夹带阶段和结束阶段。

2) 入口速度较小时，不会引起气体夹带的发生，随着入口速度的增加，气体夹带的趋势增强。

3) 随着液面的初始高度增加，发生气体夹带的最小速度增大，液面初始高度的增加阻碍气体夹带的趋势。

4) 入口高度越大，从入口进入区域的流体更多从出口直接流出，逆时针流场区域变小，导致自由液面波动趋势减弱。

本研究为铅铋池式快堆区域结构优化以及运行工况参数设定提供建议。

基金项目

南京航空航天大学2022研究生科研与实践创新计划项目(xcxjh20220610)。

参考文献

NOTES

^*通讯作者。

参考文献