1. 引言
随着21世纪航天技术的大力发展，大推力载运火箭技术得到广泛应用。大推力载运火箭多采用液氢液氧作为推进剂，因此，液氢液氧等低温液体在运输过程中的稳定性越来越被人们所重视，并成为传输工程成败的关键。低温液体在运输过程中会产生间歇泉现象[1] ，是指当低温液体在垂直管路运输过程中，由于外部容器漏热作用发生气化现象并产生气泡，当气泡不断增加会产生管路阻塞，并将液体挤出管路从而发生类似喷泉的现象，与此同时，管路上方的液体由于重力作用重新冲入管路，并且与管内的气体发生冷凝作用，并产生类似水锤的压力波动，这种压力波动会损害阀门、管路等部件 [2] [3] 。
自上世纪50年代起，间歇泉现象在多相流领域被很多学者所重视，Murphy [4] 针对长管路中的间歇泉展开了实验研究，给出了间歇泉和非间歇泉区别的曲线及经验公式。Kuncoro [5] 采用水和R113作为工质，测量出了管内发生间歇泉时的温度和压力分布数值，并认为间歇泉主要是由于液体过热所引起的，间歇泉产生的原因和管路的几何结构及温度分布有直接影响。Takashi [6] 采用水作为工质，研究了间歇泉产生的频率和强度，并用声波来控制间歇泉的产生。Khazaee [7] 以家村作为工质，研究了管壁温度与间歇泉产生周期的关系，并得出了两者的经验公式。Morioka [8] 以R113作为工质，研究了从底部和顶部向储存罐中添加R113时的间歇泉现象，得出了漏热量与流体速度的经验公式，并指出了间歇泉发生的主要区域。Angelo [9] 在NASA的报告中为了避免低温液体在运输工程中间歇泉的危害，研究了几种不同的储存罐对低温液体运输的影响。Tong [10] 通过实验手段研究了间歇泉产生过程中气泡的产生规律，提出间歇泉气泡的产生分为两个阶段，分别为产生阶段和生长阶段，并根据动力学方程列出了气泡脱离的经验公式。在此基础上，Mikic [11] 得到了统一的两阶段增长经验公式。
综上所述，间歇泉的流动具有复杂性和不稳定性等特点，特别是在低温液体工质中的间歇泉研究需要深入研究，本文在前人研究的基础上使用实验手段和数值模拟手段，对低温下垂直管路的间歇泉现象进行了实验和数值模拟研究，研究了气泡进入上部储存罐时气泡的破碎、聚合现象及底部的压力变化，为进一步研究低温液体运输安全提供可靠依据。
2. 数学模型及计算方法
本文选用VOF方法进行数值模拟，在每个控制体积单元里，所有相的体积分数之和等于1。假设用
表示第q相的体积分数，则在某控制体积单元，
的值有下列三种情况：
、
或
，分别表示该单元没有第q相、该单元只有第q相或该单元是有第q相和其它相的自由表面单元，且
满足方程
(1)
式中：n表示总相数。
每个控制体内混合流体的密度可表示为：
；每个控制体内混合流体的粘性可表示为:
。
各相之间界面的跟踪可通过求解体积分数的连续方程来完成：
(2)
对于基本相，其体积分数不是体积分数方程(2)求得的，而是由方程(1)得到。
VOF法在计算多相流时，只求解单一的动量方程，得到的速度场为各相共享，而各相的体积分数对动量方程的影响是通过属性
和
反映到动量方程中的，如式(3)所示：
(3)
3. 间歇泉数值模拟
本文使用FLUENT软件对垂直管路中的间歇泉现象进行数值模拟，并假定气泡进入存储罐时气泡长度不发生变化，气相和液相的温度相同，两相之间不发生传热和传质现象。几何模型与图2中几何模型一致，采用四边形网格，经过网格无关性测试后，采用网格数34,026的网格进行分析。两相流模型采用几何重构VOF模型，每个相的体积分数之和为1。采用PISO法对压力和速度进行耦合，表面张力采用Brackbill [12] 表面张力模型，其Courant数为0.25。边界条件和初始条件如图1所示，Ls为气泡长度，b为液膜厚度，Δy为气泡间距，壁面采用无滑移边界条件，存储罐出口为标准大气压。液相工质为液氮，其密度为807.779 kg/m3 (0 kPa)，粘度为0.164 × 10-3 Pa∙s，气相工质为氮蒸汽和空气，空气密度为1.185 kg∙m−3，粘度为1.79 × 10−5 Pa∙s，氮蒸汽密度为4.572 kg∙m−3，粘度为5.2557 × 10−6 Pa∙s。
4. 结果分析
图2为管路直径为0.018m，长度为5.5 mm单弹状气泡运动图像和瞬态压力图，a、b、c、d的图像与压力图中的a、b、c、d相对应。可以看出，弹状气泡进入存储罐后向上部液面移动，在移动过程中形成团状气泡，a、b、c处的压力波动较小，d处的压力波动值最大，幅度为0.093 kPa，此时无间歇泉现象发生。
图3为管路直径为0.018 m，长度为7.6 mm单弹状气泡运动图像和瞬态压力图，a、b、c、d的图像与压力图中的a、b、c、d相对应。可以看出，气泡的运动轨迹与图5中一致，最大压力波动幅度为1.098 kPa，此时无间歇泉现象发生。
图4为管路直径为0.018 m，长度为10 mm单弹状气泡运动图像和瞬态压力图，a、b、c、d的图像与压力图中的a、b、c、d相对应。可以看出，弹状气泡进入存储罐后向上部液面移动，在移动过程中形成団状气泡，a、b、c相应压力的波动幅度较小，d、e、f、g相应压力的压力波动很大，g点之后压力的波动幅度又变小。最大压力波动幅度为1.2 kPa，说明已有间歇泉现象发生。
图5为管路直径为0.018 m的连续多个弹状气泡底部压力的实验结果与数值模拟对比图，可以看出，
Figure 1. The boundary conditions of computational domain
图1. 计算区域的边界条件
单个气泡时，压力波动在0.2 s处有波动值，其峰值为10.46 kPa，当管内弹状气泡增多时，其整体压力幅值增大，压力波动的频率增大，数值模拟结果与实验结果相比，趋势相吻合，这是因为实验过程中无法确定弹状气泡的个数，因此只能通过压力幅值的趋势进行比较。
图6为初始气泡长度为10 mm时，不同的系统压力下管路底部与系统压力之间的关系，可以看出，
随着系统压力的增大，压力波动频率降低，压力信号的波动幅度减小，这表明大的系统压力能抑制间歇泉的发生。
5. 结论
本文采用数值模拟相结合的方法对垂直管路内的间歇泉现象进行了可视化动态特性研究，得到以下结论：
1) 弹状气泡由垂直管内进入上部存储罐时，破碎成为团状的大气泡，并沿着储罐中心位置向液面移动，在移动过程中，团状气泡发生破碎，并产生摇摆，最终溢出液面，液面伴随有液滴溅出，且液面高度发生周期性变化。
2) 长径比较小、热流密度很大工况下不会产生间歇泉。
3) 不同长度的弹状气泡对间歇泉有直接影响，当初始弹状气泡长度较小时，弹状气泡进入上部贮罐破碎成一个大的球团状气泡，引起管路底部压力波动较小，即间歇泉的强度小；初始弹状气泡长度较大时，弹状气泡进入上部贮罐破碎成多个大的云团状气泡，引起管路底部压力波动较大，即间歇泉的强度大。
4) 当垂直管内存在多个弹状气泡时，随着管内弹状气泡的增多，压力波动的频率增大，压力值整体
(a) 单个弹状气泡 (b) 2个弹状气泡
(c) 4个弹状气泡 (d) 实验结果 [13]
Figure 5. Compare between multiple consecutive slug bubbles simulation results and experimental data
图5. 多个连续弹状气泡模拟结果与实验数据的对比
Figure 6. Dynamic pressure signals at the bottom of the tube with different system pressure
图6. 不同系统压力下管路底部的动态压力信号
增大
基金项目
国家自然科学基金资助项目(50476015)。