虹吸管路内气泡动力学行为数值模拟

李兴雨 李琳 谭义海

摘 要：为了研究虹吸管中气泡运动行为与变形特征，对虹吸管中气泡运动过程进行数值模拟，结果表明：上行管内气泡纵横比越小，气泡速度越大，流量的变化对气泡纵横比影响很小，且气泡形状与静水中气泡相图相吻合，气泡在上升过程中会发生横向振荡，振荡幅度随着气泡直径的增大先增大后减小;下行管内流量不同，气泡运动方向也不同，气泡变形情况与上行管内一致;中行管内气泡为贴壁状态流动，随着气泡直径增大，气泡速度减小，当气泡直径增加到一定数值时，气泡会发生破裂，出现较小的气泡，气泡速度又出现增大的现象。

关键词：虹吸管;数值模拟;气泡速度;纵横比;振荡幅度

中图分类号：TV672.5 文献标志码：A

doi：10.3969/j.issn.1000-1379.2020.12.032

Abstract：The aim of this study was to detect the movement behavior and deformation characteristics of the air bubbles in siphon pipeline. It was essential to use numerical methods to figure out the movement of the air bubbles in the siphon pipeline. The results show that the bubble aspect ratio in the ascending tube is smaller and the bubble velocity is bigger. The change of flow has little effect on the aspect ratio and the bubble shape is consistent with its phase diagram in still water. During the progress of rising， the bubble oscillates laterally and the amplitude of oscillation first increases and then decreases with the increase of bubble diameter. The flow shows a wide difference in the downpipe， the movement direction of the bubble is also different and the bubble deformation is in agreement with the ascending tube. The bubbles in the middle pipe flow in the way with attaching the sidewall. Meanwhile， the bubble diameter increases with the decreasing of the bubble velocity. The bubble breaks at a certain diameter and smaller bubbles appear， then bubble velocity increases again. The results can provide a reference for the study of gas-liquid two-phase flow in siphon pipeline.

Key words： siphon; numerical methods; bubble velocity; aspect ratio; oscillation amplitude

水平管段距离长、真空度大的虹吸管输水管道是坎儿井式地下水库的重要组成部分。李琳等[1]、许史等[2-3]在新疆台兰河某地下水库长距离虹吸管道水力学模型试验中发现，安装高度小于7 m时，虹吸管内出现伪空化现象，流动介质由单一液相流转变为气泡流、过渡流和气团流。王梦婷[4]、张小莹等[5]的试验结果表明，虹吸管安装高度为6～8 m时，各水头运行时上行管均以气泡流为主，进入水平管路后，低水头运行时气泡在运动过程中不断聚合形成气团和长度为30～60 cm的大气囊，气囊体积随运行时间增加而增大，当其贯穿整个断面时管路断流。为了保证虹吸管正常运行，应尽量避免气泡在一定条件下聚合成气囊，而这一问题与气泡在动水和负压运行的管道内的动力学行为密切相关。

气泡以分散相的形式在水中的运动过程是一种复杂、不稳定和非线性的水动力学现象。Sanada等[6]、Wang等[7]通过试验研究和数值模拟对不同黏度液体中的气泡运动进行了研究，得到了不同直径下气泡的形状。程军明等[8]对气泡在静水中上升和破裂过程进行了数值模拟，获得了气泡上升速度和气泡形状之间的关系。陈如艳等[9]通过试验研究得到水平管内气液两相流中液相流速、气相流速和气泡尺寸的相互影响规律。Bhaga等[10]研究了浮力气泡在黏性液体中的上升运动过程，得到了气泡形状及终端上升速度与雷诺数Re、莫顿数Mo和埃奥特沃斯数Eo之间关系的气泡相图。Ellingsen等[11]使用高速摄像机对气泡在各种局部流动条件及不同水质情况下的动力学特性进行了研究，并得到了在不同条件下气泡的运动及形变过程。Ohta等[12-13]先后运用VOF数值模型，揭示了当流动条件为低Mo和Eo时，气泡上升行为与其初始形状密切相关。前人的研究主要集中于静水中气泡的运动特性，而对于动水和负压运行下虹吸管内的气泡运动特性尚未见报道。本文应用标准k-ε模型，同时结合VOF方法对虹吸管上、中、下行管路中单个气泡在不同流量条件下的运动行为进行模拟，为探明虹吸管内气泡聚合形成气囊和气泡独立运动产生条件奠定基础，同时为实际工程中合理布置虹吸管路排气装置、保證管路正常运行提供设计依据。

1 湍流数学模型

1.1 控制方程

1.2 标准k-ε模型和VOF法

在整个虹吸管内，除了弯头处以外，其他地方流线基本相互平行，因此选择计算量小、计算精度和收敛性均较好的标准k-ε模型进行计算。湍动能k和耗散率ε方程为式中：σk、σε分别为湍动能k和湍动能耗散率ε对应的Prandtl数，均为1.39;t为时间;Gk为湍动能k的生成项，Gk=μt（uixj+ujxj）uixj，其中μt=ρCμk2ε，Cμ=0.084 5;C1ε=1.44;C2ε=1.92;C3ε=0.09;Gb是因浮力影响引起的湍动能而产生的，Gb=βgiμtσtuixi，β=0.012，σt为紊动普朗特数，其他参数取值见文献[15];YM为有压缩流脉动膨胀对总耗散率的影响值，在本文中模拟的液相认为是不可压缩的，因此不考虑该项。

精确描述气液两相的运动界面是研究气泡行为的关键。VOF方法追踪的是网格中流体体积，具有容易实现、计算量小和精度高等优点，因此本文采用VOF法追踪气液交界面。当αw=0时，表示管道内没有水，被气充满;当水的体积分数αw=1时，表示管内被水充满;当0≤αw≤1时，表示管内被水和气充满。αw的控制方程为

2 计算区域的离散及边界条件

由于本文模拟的管道为圆形截面，气泡为球形，管道与气泡都具有轴对称性，因此建立二维模型进行数值模拟。模型由上行管、中行管和下行管三部分构成，其中上行管和下行管高度均为6 m，中行管长度为18 m，管径为2 cm。分别对单个气泡在上行管、中行管、下行管内运动过程进行数值模拟研究。虹吸管模型边界采用速度进口边界条件，液相初始速度分别为0.06、0.15、0.2 m/s且垂直于管路的横断面，气相的初始速度为0。管道出口为压力边界，相对压强为0，管道壁面设为无滑移固壁，采用标准壁面函数进行修正。当t=0时，气泡位于上行管距离管口50 cm中心处。计算区域网格划分采用四边形网格和三角形网格，由于整体划分网格数目较多，计算时间过长，因此对上行管、中行管和下行管网格局部加密。气泡在不同管段运动时，需运动一段距离后速度才能达到稳定。在上行管中，对气泡起点以上3 m范围内的管道进行加密;当气泡在中行管运动时，对距离气泡起点0.5 m处之后的3 m范围内的管道进行加密;下行管中，由于流量不同，气泡运动方向也不同，因此针对不同流量分别做了加密方案，流量较小时，对气泡起点以上3 m范围内的管道进行加密，流量较大时，对气泡起点以下3 m范围内的管道进行加密;网格单元最小尺寸为4×10-4m。控制方程的离散采用有限体积法，从稳定性、精确性、适用性方面考虑，控制方程中的对流项和扩散项均采用二阶迎风格式进行离散。离散后的线性代数方程组采用SIMPLEC算法迭代求解，计算时间步长为1×10-4 s。

使用标准的k-ε模型和VOF法，通过FLUENT软件进行模拟，设置表面张力系数σ=0.072 8 N/m，模拟虹吸管输水流量分别为0.018 8、0.047 1 L/s，气泡直径分别为1、4、5、6、7 mm的管路内气泡运动行为特征。

3 计算结果与分析

3.1 数学模型的验证

采用文献[16]和文献[17]的试验结果对本文的数学模型及经验参数进行验证。

3.1.1 物理模型试验装置

文献[16]试验装置见图1。试验装置高度为1 m，长、宽分别为20.32、7.56 cm，液面高度为0.8 m。使用空气注射器注入所需大小的气泡，通过水箱与管道的压力差使得气泡进入试验装置。气泡出口位置为半球形倒置的杯罩，手动旋转杯罩释放气泡，使用高速摄像机在测量区域进行测量。

3.1.2 数学模型的建立与对比

文献[18]中指出计算区域的宽度大于4倍气泡直径时，可忽略边界对气泡运动特性的影响。为了验证拟采用的数学模型，先建立高度为1 m、直径为50 mm的一端封闭一端开口的竖直管道模型，出口边界设为压力出口（相对压强为0），网格尺寸分别为6×10-4、5×10-4、4×10-4、3×10-4 m，分别模拟了直径为4、6、8 mm的气泡在静水中的上升和变形过程，气泡直径为6 mm时不同网格尺寸下气泡速度随时间变化情况见图2。由图2可知，网格尺寸越小，气泡速度越快，当网格尺寸小于4×10-4 m时，气泡速度不再随着网格尺寸的变化而发生变化。网格尺寸越小，需要的网格数量越大，计算耗时就越长，为了减少计算时间，选择网格尺寸为4×10-4m进行计算。由图2还可得到，网格尺寸为4×10-4m时模拟的气泡速度与实测的结果吻合较好。当t=0.15 s时，气泡速度模拟结果与试验结果相差最小，相对误差为1.7%;当t=0.25 s时，气泡速度模拟结果与试验结果相差最大，相对误差为4.9%。

将气泡直径为6 mm时的计算结果与文献[17]中实测的气泡相图及文献[16]的试验结果进行对比，结果如表1和图2所示。表1为不同直径气泡在上升过程中不同时刻所对应的Re和Eo，并将气泡形状模拟结果与文献[17]结果进行对比。从表1中可以看出3種直径的气泡运动过程与文献[17]中气泡相图的结果一致。文献[19]中提到，气泡在流体中上升时，气泡形状主要取决于气泡雷诺数Re、埃奥特沃斯数Eo和莫顿数Mo三个无因次准数，它们的表达式分别为

3.2 气泡竖直运动和水平运动过程及形状变化

在负压条件下，不同流量时虹吸管上、中、下行管内直径d分别为5、6、7 mm的气泡运动过程及形状变化见图3～图5。气泡形变过程用不同时刻的纵横比（即气泡短轴与长轴的比值）变化来表示，t=0时，纵横比为1，气泡纵横比越接近1表明气泡变形量越小。

图3为上行管内气泡运动变形过程，t=0时气泡为球形。当流量为0.018 8 L/s时，随着气泡上升，气泡变成椭球状，随着直径的增大，变形越来越明显，如d=5 mm的气泡经过4 s后，纵横比由t=0时的1变为0.64，而d=7 mm的气泡经过相同时间纵横比由t=0时的1变为0.55。当流量由0.018 8 L/s增大至0.047 1 L/s时，气泡仍然保持椭球状上升，纵横比随气泡直径增大而减小，但是流量的增大对气泡形状的变化影响很小，如d=5 mm和d=7 mm的气泡经过6 s后纵横比分别由1变为0.63和0.56，与流量为0.018 8 L/s时的纵横比基本相同。

图4为球形气泡在下行管随水流运动时气泡变形过程。气泡在下行或上行过程中主要经历了球形和椭球形的变形过程。计算不同时刻的Re、Eo和Mo，与文献[17]的气泡相图进行对比，结果表明：动水条件下虹吸管上、下行管内气泡在不同时刻的形状变化过程模拟结果与文献[17]的静水条件下气泡运动变形过程一致。当流量为0.018 8 L/s时，因流速较小，流体曳力（物体在流体中有相对运动时，会受到流体的阻力，阻力大小由相对速度差决定）较小，在浮力作用下d=5 mm的气泡在下行管内向上运动，气泡运动速度小于上行管中气泡运动速度，此时下行管气泡纵横比较大，气泡运动过程中变形量略小于上行管内的。如t=4 s时，d=5 mm的气泡在下行管内气泡纵横比比上行管内增加2%。当流量为0.047 1 L/s时，流体曳力大于浮力，气泡随水流一起向下运动，随着气泡速度的增大，相对速度差减小，流体曳力减小，下行管气泡运动速度小于上行管内的，因此气泡变形量仍小于上行管内的气泡变形量，经过4 s后，d=7 mm的气泡纵横比比上行管中增大了9%。

图5为气泡在中行管运动变形过程，不同流量下气泡变形过程基本一致。当t=0时，气泡为球形。当t>0时，在浮力和壁面黏附力作用下气泡变为半椭球状，自管道中心向上管壁运动，保持贴壁流运动。当Q=0.018 8 L/s时，经过相同时间，气泡纵横比随气泡直径的增大而减小。如d=5 mm和7 mm的气泡经过4 s后纵横比分别由1变为0.24和0.17，对同一直径气泡而言，经过相同时间，气泡纵横比随流量增大而减小。如Q=0.018 8 L/s时，d=5 mm的气泡经过4 s后纵横比较Q=0.047 1 L/s时增大了11%。

3.3 虹吸管内不同直径气泡速度与纵横比

图6为Q=0.018 8 L/s时不同直径气泡的速度、纵横比随时间的变化情况。其中：v、v1、v2分别为气泡直径为5、6、7 mm时的气泡运动速度，e、e1、e2分别为气泡直径为5、6、7 mm时的气泡纵横比。

气泡直径一定时，纵横比越大，气泡速度越小，这一规律与文献[16]中提到的静水条件下气泡运动过程中纵横比与其速度成反比的规律一致。流量增大时，气泡纵横比与速度的关系和小流量时一致（限于篇幅，文中未给出其他流量时的图形）。从图6（a）可以看出，气泡直径越大，速度越小，但不同直径气泡的速度相差不大，如d=5 mm和d=7 mm时，最大差为t=1.0 s时的12%，最小差为t=2 s时的2%。其原因是气泡直径增大，所受浮力增大，在浮力作用下，气泡与液体之间的相对速度差增大，导致流体曳力增大，抑制气泡上升，随着气泡直径的增大，管壁黏附力将会阻碍气泡运动，使得大、小气泡的运动速度差别不明显。

从图6（b）可以看出，不同直径的气泡在下行管内的流速均小于上行管内的，且同一时刻气泡纵横比越大，速度越小。Q=0.018 8 L/s时，在下行管中，气泡向上运动，由于水流方向与上行管内水流方向相反，因此速度比上行管内小。气泡速度随着气泡直径的增大略有减小，如气泡直径由5 mm增加到7 mm时，气泡速度由0.16 m/s减小19%，其原因是较大的气泡所受管壁黏附作用影响较大，使得阻力增大。当流量增加到0.047 1 L/s时，流体曳力增大，气泡向下运动，气泡的纵横比与速度的关系也与上述一致;气泡直径对气泡速度的影响很小。

从图6（c）可以看出，在中行管内不同直径的气泡先进行短暂的减速运动后达到终速度，与上行管相比速度有所减小，其原因是在中行管内气泡先向上运动，直到与管壁贴合，壁面黏附力增大，浮力不再做功，气泡速度减小。当直径d=5～6 mm时，气泡速度会随着气泡直径的增大而减小。当气泡直径进一步增大，如d=7 mm时，气泡速度增大，其原因是气泡直径过大，运动过程中发生了破裂，形成了较小的气泡，使得速度又出现增大的现象。由于中行管内是贴壁流，因此气泡纵横比与上行管相比相差很大，如Q=0.018 8 L/s时，d=5 mm的气泡经过4 s后，中行管内气泡纵横比较上行管内减小65%。

3.4 气泡运动轨迹

气泡在上升过程中会发生横向位移，其原因是气泡在上升的过程中会发生形变使得受力不平衡，导致了横向振荡。将气泡中心点的位置坐标（0，0.6）作为t=0时刻气泡的运动位置，流量为0.047 1 L/s时不同直径的气泡在上升过程的运动轨迹见图7，图中x/d表示气泡的中心点位置沿管径方向的相对位置，y/hs表示气泡的中心点位置与虹吸管高度的相对距离。由图7可知，气泡在上行管中运动时，随着气泡直径的增大，振荡幅度发生变化。当上升气泡直径较小时，如d=1 mm时，气泡上升中的振荡幅度较小，基本沿直線上升。随着气泡直径的增大，如d=4 mm时，气泡呈“之”字形上升。当气泡直径增大到7 mm时，气泡的振荡幅度却开始减小，其原因是直径较大的气泡在上升过程中受到壁面黏滞作用的影响增大，对气泡的振荡起阻碍作用，使得振荡幅度变小。在不同流量的情况下，气泡的运动轨迹呈现相似的变化。

4 结 论

本研究模拟了单个气泡在静水中的上升过程，并与文献[17]的气泡相图结果进行比较，验证了模拟结果的正确性，在此基础上对虹吸管道水流进行数值模拟研究，并对模拟结果进行分析。

（1）气泡在上行管上升过程中，由球状变为椭球状，纵横比逐渐减小，气泡直径越大，纵横比变化越剧烈，而流量的变化对气泡纵横比的影响很小。气泡速度与纵横比有关，纵横比越大，气泡速度越小。气泡直径较小时，气泡沿直线上升，随着气泡直径增大，气泡发生振荡，呈“之”字形上升。当气泡直径增大到一定程度时，气泡受到壁面的黏滞作用影响增大，对气泡的振荡起阻碍作用，使得振荡幅度变小。

（2）气泡在下行管运动过程中，气泡变形情况与上行管相似，但气泡受到浮力和曳力作用，下行管内气泡速度和变形程度小于上行管内的。

（3）气泡在中行管运动过程中，气泡经过短暂时间运动至管道顶部，在此之后为贴壁运动，气泡由球状变为半椭球状。随着气泡直径的增大，气泡纵横比逐渐减小;纵横比大小与流量有关，流量越大，纵横比越小。气泡速度随气泡直径增大而逐渐减小，气泡发生破裂后的气泡速度大于原来气泡的速度。

参考文献：

[1] 李琳，邱秀云，许史，等.长距离虹吸管道输水水力学模型试验研究[J].南水北调与水利科技，2010，8（3）：106-109.

[2] 许史，李琳，邱秀云，等.长距离虹吸管输水试验研究初探[J].中国农村水利水电，2010（3）：70-72.

[3] 许史.长距离虹吸管输水试验研究[D].乌鲁木齐：新疆农业大学，2010：13-24.

[4] 王梦婷.真空管道气液两相流动水力特性试验研究[D].乌鲁木齐：新疆农业大学，2015：19-69.

[5] 张小莹，李琳，谭义海，等.虹吸管道坡度对气液两相流动特性影响的试验研究[J].农业工程学报，2017，33（14）：122-129.

[6] SANADA T，SUGIHARA K，SHIROTA M，et al. Motion and Drag of a Single Bubble in Super-Purified Water [J]. Fluid Dynamics Research，2008，40（7/8）：534-545.

[7] WANG Huanran，LI Yang，YANG Dong，et al. On the Shape Feature of a Single Bubble Rising in Viscous Liquids [J]. Engineering Thermophysics，2009，30（9）：1492-1494.

[8] 程军明，吴伟烽，聂娟，等.气泡在静水中上升破裂产生射流特性的数值模拟[J].江苏大学学报（自然科学版），2014，35（5）：513-517.

[9] 陈如艳，仇性启，王晓琦.水平管内气液两相流研究[J].石油化工设備，2007，36（6）：37-38.

[10] BHAGA D，WEBER M E. Bubbles in Viscous Liquids： Shapes，Wake and Velocities [J]. Fluid Mech， 1981，105：61-85.

[11] ELLINGSEN K，RISSO F. On the Rise of an Ellipsoidal Bubble in Water：Oscillatory Paths and Liquid Induced Velocity [J]. Computer Phys， 2008，227：3358-3382.

[12] OHTA M，IMURA T， YOSHIDA Y， et al. A Computational Study of the Effect of Initial Bubble Conditions on the Motion of a Gas Bubble Rising in Viscous Liquids [J]. Multiphase Flow，2005，31（2）：223-237.

[13] BONOMETTI T，MAGNAUDET J. Transition from Spherical Cap to Toroidal Bubbles [J]. Phys Fluids，2006，18（5）：1-12.

[14] BRACKBILL J，KOTHE D. A Continuum Method for Modeling Surface Tension [J]. Journal of Computational Phvsics，1992，100：335-345.

[15] ISSA R I，OLIVEIRA. Numerical Prediction of Phase Separation in TwoPhase Flow Through TJunctions [J]. Computers Fluids，1994，23（2）：347-372.

[16] 刘柳.垂直上升管中气泡动力学特性实验研究[D].长沙：中南大学，2013：13-17.

[17] 李仲春，宋小明，姜胜耀，等.静水中较大气泡运动特性实验研究[J].核动力工程，2015，36（1）：161-164.

[18] 何丹，李彦鹏，刘艳艳.初始形状对浮升气泡动力特性的影响[J].西安交通大学学报，2011，45（1）：43-47.

[19] 刘静如.非牛顿流体中多气泡相互作用、聚并与破裂过程的数值模拟[D].天津：天津大学，2014：1-10.

【责任编辑 许立新】