文丘里管式微气泡发生器内单气泡碎化行为的数值模拟

丁国栋，陈家庆，蔡小垒，叶帆，李振林，姬宜朋，郭靖

（1 北京石油化工学院机械工程学院，北京102617；2 中国石油大学（北京）机械与储运工程学院，北京102249；3 中国石化西北油田分公司石油工程技术研究院，新疆乌鲁木齐830011）

管式紧凑型微气泡发生器作为一种新兴微气泡发生手段，正受到越来越多的关注[1-2]。其中文丘里管式微气泡发生器更是广泛应用于过程强化、水质净化[3-4]、矿物浮选[5]、水体增氧[6]等工程领域。尽管该类型微气泡发生器具有结构简单、携气量大等优点，但同时也存在所产生微气泡的平均粒径偏大、粒径分布不均匀等缺点。虽然国内外学者在文丘里管式微气泡发生器成泡粒径测量表征和应用效果等方面已开展了大量研究，但在气泡破碎过程及成泡机理方面的研究工作相对较少，具有代表性的仅有日本筑波大学Fujiwara[7]和Nomura[8]团队、哈尔滨工程大学阎昌琪教授团队[9-10]、四川大学孙立成教授团队[11-13]。更为重要的是，相关研究均采用可视化实验方法进行，但受限于所采用的实验方法和测量仪器，仅能观测表观气泡的破碎现象，无法全面、准确获取内部流场信息，进而无法准确分析气泡破碎原因；此外也不可否认，实验过程中的不可控因素较多，实验结果的偶然性和随机性较强。随着计算机软、硬件技术的迅速发展，计算流体动力学（CFD）数值模拟已经成为实验测量的有效补充，并应用于气液两相流的研究。例如：西安交通大学的Guo 等[15]运用VOF 模型对T形微流道内泰勒气泡的生成过程进行数值模拟研究，对比分析了液相黏度、表面张力和气液速度等因素对微流道内气液两相流流动特征的影响；河北工业大学的Chen等[16]对射流式微气泡发生器的气液两相射流流动进行数值模拟研究，分析了流场内不同横截面处的速度、压力和湍动能的变化情况；中国矿业大学的Wang 等[17]借助ANSYS Fluent 商用软件对微气泡发生器内的气液两相流动进行三维流场数值模拟，并对喷嘴与喉管的间距进行了优选。在气泡破碎方面，美国克拉克森大学的Qian 等[18]运用格子-玻尔兹曼方法对零重力条件下均匀湍流场中的气泡变形和破碎现象进行了数值模拟研究，分析了气泡破碎与韦伯数之间的关系；荷兰埃因霍芬理工大学的Lau等[19]采用Euler-Lagrange框架进行泡状流下气泡破碎特征数值模拟研究，并对湍流场中气泡粒径分布进行预测。在微气泡发生器中单气泡数值模拟方面，哈尔滨工程大学的唐文偲等[20]通过在文丘里式气泡发生器扩张段壁面处设置一个刚性小球模拟气泡，采取观察小球周围的速度和压力变化，分析小球受力情况的方式，间接探究气泡的破碎特征；大连海事大学Li 等[21]则借助Fluent 6.3 模拟二维文丘里流道中单气泡在不同进水压力下的变形及破碎过程，指出湍流波动和剪切应力是造成文丘里流道内气泡破碎的原因。不难看出，迄今研究主要针对宏观气液两相流来进行，气泡破碎过程的数值模拟研究也仅针对二维或微型文丘里流道。有鉴于此，本文将借助CFD 商业软件ANSYS Fluent 中的VOF 多相流模型，对文丘里管式微气泡发生器内单气泡的破碎过程进行研究，通过模拟初始气泡在三维流场中的变形和破碎行为，探究初始气泡位置、水流流速、扩张段张角等因素对气泡碎化特性的影响，讨论分析文丘里管式微气泡发生器生成气泡粒径不均匀的原因，提出内置安放带导流叶片轴向静止起旋元件的改进方案，为文丘里管式微气泡发生器的结构优化设计和发泡性能提升提供参考。

1 CFD数值模拟模型的构建

1.1 模型结构和碎化原理

图1 文丘里流道的简化结构示意图

文丘里管式微气泡发生器的主体结构由注气管、进水管、收缩段、喉管、扩张段和出水管等组成。工作过程中，外界气源由注气管注入并在水流裹挟下生成初始气泡，初始气泡向下游运动并在水力作用下被剪切破碎成微米级气泡。本文主要开展气泡在文丘里流道内破碎过程的研究，因此省略注气成泡段，仅保留由进水管、收缩管、喉管、扩张段和出水管构成的文丘里流道，简化结构形式及流动区域特征尺寸如图1所示。

一般而言，初始气泡发生碎化是由于外部作用力的效果强于气泡内部维持稳定的作用力，所生成微气泡的粒径主要由紊流惯性力Fi和表面张力Fs共同决定，如式(1)、式(2)。紊流惯性力Fi由单位体积的能量耗散率ε决定[20]，适用条件为湍流中液体受到的平均力与施加在相同大小的流体元件上的力相同[22]。韦伯数We=Fi/Fs，并由此得到微气泡粒径表达式如式(3)[20]。

式中，ρ为流体密度；d为气泡直径；g为正比例常数；σ 为表面张力系数；ε 表示能量耗散率，能量耗散率越大，生成微气泡的粒径越小。

与直管流道相比，文丘里流道“收缩+扩张”的特殊结构形式使分散相气泡除受到湍流引起的作用力之外，也受到压力梯度力等其他作用力的剧烈影响[23-29]。日本筑波大学Fujiwara等[7]的实验研究发现，气泡在文丘里流道喉管处发生变形，进入扩张段后，由于气泡头部和尾部位置处的压力差使得气泡发生破碎，结果说明扩张段内静压回升是造成气泡发生碎化的原因；四川大学孙立成教授课题组[12]基于文丘里流道内单气泡受力的初步研究结果表明，在扩张段内激增的压力梯度力和附加质量力导致气泡运动速度减小，并与周围流体形成强烈相互作用，进而促进气泡发生变形和破碎。

1.2 CFD数值模拟方法

1.2.1 几何模型与网格划分

图2 文丘里流道三维模型及网格图

CFD数值模拟用几何模型的结构尺寸与实验用文丘里管式微气泡发生器的结构尺寸完全一致，同时为消除进、出口段对内部流场的影响，将进水管和出水管延长至100mm。以进水截面中心原点为参考零点，三维图形及网格划分如图2所示，计算域网格为O形六面体结构。

1.2.2 数值模拟模型选择与边界条件设置

为研究文丘里流道内气泡的变形及破碎过程，多相流模型选择VOF 模型。ANSYS Fluent 商业软件的湍流模型是基于Reynolds时均方程而建立，常用类型主要包括标准k-ε模型、RNG k-ε模型和雷诺应力模型等[30-31]。考虑到RNG k-ε 模型相比于标准k-ε模型在ε方程中添加了一个条件，有效改善了精度并考虑到湍流中的涡旋，拥有更高的可信度和精度；同时考虑到雷诺应力模型计算时需要更多的内存消耗，并且由于方程间的强耦合性而致使收敛性差，更适合复杂强湍流工况。因此，本文选择RNG k-ε湍流模型来模拟文丘里流道内的湍流流动过程。

参照相关学者数值模拟的设置经验[32-34]，同时考虑到本文重点关注雷诺数对气泡破碎的影响以及压降变化，因此设置入口边界条件设为速度入口（velocity inlet），出口边界条件设为压力出口（pressure outlet）；同时，为简化模型，去除其他无关因素对模拟结果的影响，设置壁面为无滑移壁面。湍流定义方法选择水力直径和湍流强度来进行，水力直径等于文丘里流道的入口直径，而湍流强度I的计算公式为式(4)。

工作介质为水（主相）和空气（次相），主相水的密度为998kg/m3，次相空气的密度为1.225 kg/m3，气液两相间的黏性系数为0.075n/m。采用控制体积法离散控制方程，压力和动量分别采用PRESTO！格式和二阶迎风格式离散。

1.3 网格无关性验证

网格数目对数值计算结果的准确性影响很大，只有当网格达到一定数量，关键计算结果随网格数目的增加变化很小时，数值模拟计算的结果才具有意义。分别构建网格数为350370、470524、610598、771792、949784 的五组模型进行网格无关性验证。扩张段部分的湍流耗散率是本文重点关注的参量，截取扩张段中间位置（x=164mm）横截面，得到各组模型的平均湍流耗散率和相对偏差，计算结果如表1所示。由表1可知，随着网格数目的增加，相对偏差值逐渐减小，当网格数目由771792 增加至949784 时，相对偏差为0.91%，可见网格数为771792即可满足网格无关性要求。

表1 网格无关性验证结果

1.4 实体模型测试校验

为实体测试校验所建立CFD 数值模拟模型的准确性，搭建如图3所示的实验装置，对比分析实验结果与CFD 数值模拟结果。实验装置主体由水箱、离心泵、文丘里管、流量计、压力表和阀门等组成，实验过程中通过改变离心泵转速和侧线阀门开度，调节进水流量大小，观察并记录文丘里流道进、出口压力值，同时为降低偶然误差的影响，重复3次实验并取平均值，实验结果与CFD数值模拟结果如表2所示。由表可知，当选择RNG k-ε湍流模型时，进出口截面的平均压降与实验值最大偏差仅为1.14%，因此可见所建立的CFD数值模拟模型满足计算精度要求。

图3 实验流程示意图

表2 实体模型测试验证结果

为进一步证明CFD 数值模拟结果的可信性，参照文献[11]中的实验过程，设置入口雷诺数Re=20000，待单相流场稳定后，在离喉管入口x=125mm 处注入（Patch）进一个粒径为4mm 的初始气泡，观察其在喉管及破碎段的运动过程，CFD数值模拟的后处理结果如图4所示。由图可得，注入的气泡在喉管内基本保持形状不变，进入扩张段后气泡发生较大程度的变形及碎化现象。在文丘里流道内，气泡碎化过程的CFD 数值模拟结果与实验结果总体一致，但气泡在扩张段的碎化程度与实验值有所差异。分析原因主要是由于实验过程中，由注气管注入并产生的初始气泡具有一定的初速度，而CFD 数值模拟过程中Patch 进的初始气泡速度为零，此外当前采用的数值模拟手段无法完全模拟出流道内的湍流场，因此导致扩张段内的气泡模拟形态与实验观测结果有所差异[35]。

图4 气泡碎化过程云图

2 模拟结果及分析

2.1 初始气泡位置对气泡碎化特性的影响

2.1.1 轴向位置

为研究初始气泡轴向位置对气泡碎化特性的影响，设置进水流速为1m/s，待主流流场运行基本稳定后，分别在空间坐标(10,0,0)和(125,0,0)位置处Patch 进粒径为φ6mm 的单一气泡，分别代表单气泡由进水管注入和由喉管处注入。调整计算步长为0.00005s，每5步自动保存.cas和.dat文件，设置重力（gravity）方向为X 轴负方向，后处理得到如图5所示的气泡在文丘里流道内变形及破碎瞬态图。

由图5(a)可知，气泡在文丘里流道内经历微小形变、拉伸变形和破裂碎化3 个过程。由图5(b)可知，单个气泡由喉管处注入后，在喉管内并未发生明显变形，进入扩张段后，由于水流速度降低，静压回升，气泡开始发生较大变形，但并未发生气泡碎化现象。具体而言，在进水管段，气泡的运动较为平稳，仅发生微小形变；进入收缩段后，气泡开始发生拉伸变形，气泡逐渐向子弹外形演变；气泡进入喉管后，气泡尾部开始发生凹陷；进入扩张段后，气泡尾部冲破气泡头部，气泡发生破裂并生成子气泡。在进水管处Patch 进的气泡在后续流场中发生碎化的程度明显强于由喉管处Patch进的气泡，分析原因是由于，前者气泡经历更长时间的湍流碰撞和能量耗散，气泡表面的不稳定性增强，进入扩张段后更易发生碎化。

图5 不同轴向初始位置注入时的气泡变形及破碎瞬态图

图6 扩张段径向速度分布图

如图6所示为扩张段不同轴向位置（x=140mm、145mm、150mm）处的径向速度分布。由图可知，随着轴向位置由140mm 增大至150mm，液相平均速度由5.96m/s 降低至4.07m/s。液相速度的降低将导致气泡在扩张段运动过程中，气泡头部的运动速度低于气泡尾部。在惯性力作用下，气泡头部则极易切入气泡内部造成气泡破碎，形成如图5所示的破碎形态。

2.1.2 径向位置

在文丘里管式微气泡发生器内，初始气泡的生成位置并不局限于中心区域，在径向边壁区域同样有所分布。为研究径向边壁处初始气泡的碎化特性，待主相流场运行稳定后，在空间位置(10,6,0)处Patch 进粒径为φ6mm 的单气泡。如图7 所示为气泡在文丘里流道内的变形及破碎瞬时图，与图5(a)所示的气泡运动及破碎过程相比，初始气泡由边壁位置进入并在扩张段边壁处的碎化程度明显更强。由此可推测，在扩张段区域的湍流程度并非均匀分布，边壁区域远大于中心区域，这将直接导致扩张段区域径向位置气泡碎化程度不一致，生成微细气泡粒径分布不均。

图7 初始气泡位于边壁区域时的变形及破碎瞬态图

图8 扩张段横截面处的湍流参量

为定量表征扩张段区域径向方向的湍流流场，在扩张段中间位置（x=164mm）截面处得到径向方向的表观速度和湍流参数，绘制如图8所示的湍流流场分布图。由图可知，在中心区域（－0.004m≤y≤0.004m）径向范围内，表观速度值最大，湍动能、湍流强度和湍能耗散率值最小；由中心区域向边壁推移，速度值逐渐减小，湍动能和湍流强度先增加后减小。据此分析，在中间区域，水流速度不随径向位置发生改变，致使剪切作用较弱，并且由于层间脉动较小，导致湍动作用较弱。提取图8(a)所示边壁区域（-0.0085m＜y＜-0.0041m）内的速度参量并进行线性拟合，得到函数关系：u=989.09y+8.11，拟合度R2=0.97。由此可得，位于中心区域的气泡碎化程度较低，而在边壁区域速度梯度较大，剪切应力及湍动强度均较强，导致气泡碎化程度更高。由此可证明，由于扩张段区域中间位置的湍流参数远低于边壁区域，使得单气泡在扩张段中间区域的破碎程度远弱于边壁区域，中间区域生成子气泡的粒径也将大于边壁区域，进而造成文丘里管式微气泡发生器成泡粒径分布不均。

日本筑波大学Uesawa 和Nomura[36-37]在实验中发现，水流流速对湍流场内气泡破碎具有重要影响，水流速度越大则气泡碎化程度越强，产生微气泡的平均粒径越小。橡树岭国家实验室Kress[22]的研究结果表明，当液相雷诺数增加时，能量耗散率增加。图9 所示为进水流速分别为1m/s、1.5m/s、3m/s 和6m/s 时，在文丘里流道收缩段中间位置（x=164mm）处，湍流耗散率沿径向方向分布折线图。从图中可以看出，随着进水流速增加，平均湍流耗散率随着水流流速增加，但耗散率沿径向方向的分布形式并未改变，中间区域的湍流耗散率依旧很小。由式(3)可知，湍流场中的成泡粒径与湍流耗散率成反比，边壁区域内产生子气泡的粒径将降低，而在中心区域的子气泡的粒径并未发生变化，生成子气泡粒径分布不均匀的现象并不随着水流流速的增大而改善。

2.2 结构参数对子气泡粒径分布的影响

2.2.1 扩张段张角

图9 不同水流流速时湍能耗散率分布

上海交通大学Li等[38]在进行文丘里管式微气泡发生器成泡实验时发现，扩张段张角对生成微气泡的粒径分布影响较大，扩张段张角越大，则生成微气泡索特平均直径越小。分析原因是由于，扩张段张角的增大引起较大的压力梯度和速度梯度，湍流强度增强；同时扩张段张角增大将引起边界层分离，助力气泡的破碎细化。湍流耗散率（ε）与生成微气泡粒径正相关，在扩张段径向方向湍流耗散率的分布越均匀，则认为生成子气泡的粒径分布越均匀。定义标准差（σ）表征径向湍流耗散率的离散程度，标准差越小，则表明数据点与平均值差异越小，均匀度越高，如式(5)。

式中，σ为标准差；X为变量；μ为总体均值；N为总体样数。

因此，为研究扩张段张角对生成微气泡粒径分布的影响，构建扩张段张角分别为6°、8°、10°、12°和16°共5组模型进行对比研究。图10所示为不同扩张段张角时湍流耗散率的整体平均值（εa）、中间区域平均值（εm）和标准差（σ）。从图中可以看出，随着扩张段张角增大，湍流耗散率平均值（εa）增大，进而产生子气泡的平均粒径减小，与上海交通大学所得的实验结果一致。中心区域的平均湍流耗散率（εm）随着张角的增大而增大，并且在12°时达到最大；标准差（σ）则随着张角的增大而增大，均匀度降低。这说明，尽管随着扩张段张角的增大，生成的微气泡平均粒径减小，但同时也导致微气泡粒径分布均匀度降低的问题。

图10 不同张角时的湍流耗散率和标准差

2.2.2 内置静止起旋元件

图11 不同导流叶片出口角时的湍流耗散率

为提高文丘里管式微气泡发生器的成泡质量，本文作者课题组[39-40]在2017年左右研究发现，在文丘里流道前端内置安装带导流叶片的轴向静止起旋元件，可显著降低生成子气泡的平均粒径，并设计了一种轴向旋流式微气泡发生器。静止起旋元件上导流叶片的出口张角（α）对生成气泡粒径具有重要影响，为研究导流叶片出口角度对生成子气泡粒径分布的影响，构建了出口角度分别为15°、20°、25°、20°和无叶片共5 组模型进行数值模拟计算。图11 所示为不同叶片出口角时，扩张段径向平均湍流耗散率和中心区域平均湍流耗散率柱状图。从图中可以看出，在文丘里流道进水管内轴向安装带导流叶片的轴向静止起旋元件可大幅提高平均湍流耗散率，进而降低生成微气泡的平均粒径，与实验所得结论一致。同时，随着导流叶片出口角度增加，平均湍流耗散率（εa）先增加后降低，并在25°时最大，而中间区域平均湍流耗散率（εm）则随着出口角度的增大而减小，标准差（σ）也随之减小，扩张段径向方向的湍能耗散率分布更加均匀。

4 结论

本文借助商业CFD 软件ANSYSY Fluent 中的VOF 模型，对单气泡在文丘里管式微气泡发生器内的碎化行为进行数值模拟研究，进行了网格无关性验证且数值模拟结果与实验结论基本相符，在此基础进一步分析研究并得到如下结论。

（1）单气泡碎化发生在文丘里管式微气泡发生器扩张段内，湍能耗散率越大，则气泡碎化生成微气泡的粒径越小。

（2）扩张段中心区域的湍能耗散率远小于边壁区域，湍能耗散率径向位置分布差异将直接导致生成微气泡粒径分布不均。在轴向方向，由进水管注入的单气泡碎化程度强于由喉管处注入；在径向方向，由进水管偏心位置处注入的单气泡碎化程度强于中心位置处注入。

（3）进水管内安放带导流叶片的轴向静止起旋元件，不仅可提高扩张段的平均湍流耗散率，降低生成微气泡的平均粒径，而且使得径向湍流耗散率的标准差降低，增强扩张段径向湍流耗散率均匀度，进而提高生成微气泡的粒径分布均匀程度。

符号说明

d—— 气泡粒径，mm

F—— 相间作用力，kg/(m·s2)

g—— 正比例常数

I—— 湍流强度

N—— 总体样数

Re—— 雷诺数

v—— 进水流速，m/s

We—— 韦伯数

X—— 变量

α—— 导流叶片出口张角，（°）

ε—— 能量耗散率，m2/s3

θ—— 表面张力系数

μ—— 总体均值

ρ—— 流体密度，kg/m3

σ—— 标准差

下角标

a—— 整体平均值

i—— 紊流惯性力

m—— 中间区域平均值

s—— 表面张力