气泡羽流摆动特性及附壁效应研究

2020-09-15程以炫董恒瑞闫睿奎

高校化学工程学报 2020年4期

程以炫, 董恒瑞, 闫睿奎, 危卫

(1. 武汉理工大学能源与动力工程学院, 湖北武汉 430063; 2. 中交广州航道局有限公司, 广东广州 510220)

1 前言

气泡羽流是一种常见的两相流动，广泛应用于化工和冶金工业等领域[1]。气液反应器在化工领域一直被重视和研究[2-3]，气泡柱发生器需要快速、有效地产生均匀大小气泡。曝气池内需让氧气均匀分布被充分吸收，从而更好地处理污水、改变水质[4]。膜生物反应器中需提高气含率使得气体被微生物吸收利用[5]。在气隙式膜蒸馏工艺中需气液两相混合均一，提高传质效率[6]。但都存在一个共同的问题，气流在上升中会倾斜发生附壁效应，这是工程实践需要尽量避免的。

已有文献表明，气泡羽流会因附壁效应使得气泡流偏向某侧壁面甚至贴着壁面流动[7～9]。曲径[10]通过改变曝气量、曝气位置和液体性质等条件，利用摄影技术观察气泡羽流的形态，得出附壁效应的存在减轻羽流震荡。马霞等[11]采用双欧拉法模拟矩形体积内的流动特性，发现当纵横比为1.0 时，气泡羽流流动平稳；当纵横比为2.25 和3.0 时，气泡羽流的稳定性变差。王海博等[12]采用流体体积函数(volume of fluid，VOF)方法模拟1 mm 曝气孔径下曝气强度对反应器内流场特性的影响，气泡脱离孔口呈现裙状形(skirted shape)上升。徐礼嘉等[13]采用湍流模型对鼓泡塔进行非稳态模拟，观测到气泡羽流的震荡行为，这是由于模型的低雷诺特性导致。王蒙等[14]发现，当流域形状的高度增加时，气泡羽流的震荡频率和运动分布会导致气泡的结构更加不稳定。VAN LUC NGUYEN 等[15]采用三维涡单元法来计算连续流动，通过运动方程对气泡进行跟踪。较小的气泡更容易被大尺度的涡结构所携带，并且传输距离与气泡直径成反比。LIU 等[16-17]对矩形气泡柱发生器进行了无量纲分析，得出震荡特性公式。

气泡流可能在附壁效应作用下偏移并大量聚集，从而对产生的气泡分布造成影响，所以厘清水域中气泡羽流的偏移规律很有必要。然而，定量研究气泡羽流上升方式及其偏移程度的文献较少。鉴此，本文采用VOF 方法[18]对气泡的上升运动行为及偏移特性进行研究，通过设定矩形水槽模型内不同工况，研究气泡羽流偏移距离和倾斜角度，得出气泡流摆动特性与其主要影响因素之间的定量关系，并分析其上升规律。

2 数学模型与理论基础

2.1 守恒方程及湍流模型

气液两相流满足连续性方程和动量方程。本文使用基于雷诺时均(Reynolds-averaged navier-storkes equation，RANS)方法的κ-ε 湍流模型(κ-ε turbulent model)来求解两相流动的湍流特征。其湍动能κ 的输送方程为

湍动能耗散率ε 的输送方程为

式中：μ 为湍流黏度；G 为湍动能生成项；ρ为密度；κ-ε 模型中各系数的取值分别为σκ= 1.0，σε= 1.3，1Cε= 1.44，2Cε= 1.92，Cμ= 0.09。

2.2 VOF 模型及表面张力

VOF 模型假定水和空气是连续且不可压缩的两相。当需要追踪空气和水的界面时，可以采用这种模型。在VOF 模型中，动量方程被细分为各自相的动量方程，且各相体积分数总和为1。表面张力模型采用连续表面力模型，VOF 计算中附加的表面张力导致了动量方程中的源项。跨过表面的压降依赖于表面张力系数σ 和通过两个半径的正交方向量度表面曲率R1和R2：

式中：p1和p2是2 种流体界面两侧的压力。

2.3 附壁效应

附壁效应[19]又称康达效应，在流体力学中具体指的是高速射流从狭窄的通道突然进入扩大的区域时，射流会进行附壁流动而发生变化。其速度大小会减弱，速度方向会偏向于空腔一侧的壁面。这是由于射流卷吸作用使得射流两侧存在低压区。在压力差的影响下，在低压区的核心区域形成涡旋。在实际扰动的存在下，一旦射流稍微偏向某侧，该侧的涡流加强，造成该侧的低压区压力更低。在这种正反馈循环的驱使下，射流最终会附着在该侧壁面上，形成附壁效应。附壁效应的影响因素主要是外部的流域结构以及内部的射流速度和特性长度。本文研究的是水槽内气泡上升过程的附壁效应及其摆动特性，如图1所示。

图1 气泡上升过程中的摆动特性示意图 Fig.1 Schematic diagram of bubble oscillation characteristics during rising

定义无量纲修正密度弗劳德数：

式中：v 为空气的入口速度，Δρ 为水和空气密度比，g 为重力加速度(取9.81 m·s-2)，d 为空气的孔口直径。

定义无量纲偏移距离：

式中：D 为气泡偏移的最大水平距离，L 为水槽半宽。显而易见，当η=1，发生附壁效应；当η=0，气泡未发生任何偏移。

模型采用最大倾斜偏移量所在位置点、孔口中心位置点组成的线段与孔口中心垂直线的夹角定义为倾斜角度θ。

定义水槽高宽比：

式中：H 为水槽高度。

3 仿真与实验

3.1 物理模型

本文研究水槽在静水中气泡生成上升直至自由液面(最上端边界)的过程，在模型底部中心设置直径为d 的气泡生成孔，并由此射入空气(常温常压)进入水槽；为保证更准确观察完整的气泡以及对周围流场的影响，在划分网格时对气泡孔口附近局部加密。物理参数：水的密度取998.2 kg·m-3，空气和水的黏度分别取1.789 4×10-5和0.001 Pa·s，水和空气的表面张力取0.07 N·m-1。随着空气速度的增大，其密度也会随着增大。查资料可得空气速度和密度的对应关系，在此不再累赘。

3.2 边界条件、初始化及求解方案

仿真部分使用商业软件ANSYS-Fluent。1 m×1 m 模型网格尺寸为2 mm(进口局部加密为0.2 mm)、网格结构为四边形网格，网格数量为251 161。其他模型网格尺寸和网格结构一样，网格数量随之变化。计算采用标准κ-ε 湍流模型。将底部的气泡孔设置为速度入口，自由液面设置为压力出口，壁面设置为无滑移壁面。求解拉普拉斯方程初始化流场的速度和压力。计算为瞬态类型，时间步长为0.01 s。采用有限体积法二阶迎风格式将控制方程在网格上进行离散。压力速度耦合采用SIMPLEC 算法，其中动量方程和湍流方程采用二阶迎风格式离散，压力项则采用PRESTO!算法。各算例的工作参数选取如表1 所示，共150 组算例。

表1 各算例的工作参数 Table 1 Boundary conditions and working conditions of the cases studied

3.3 模型验证

通过与DEEN 等[20]实验比较，在瞬时Y 轴速度v～曲线处，仿真结果与实验数据吻合较好。在轴向位置中心，VOF 模型对Y 轴速度的估计略高；而在中心与柱壁之间的区域，轴向速度略有低估。在近壁面区域，实验值的轴向速度比模拟值速度略高，这是造成壁面效应的原因之一。在气液两相同时存在的区域，除了两相之间没有速度滑移外，数据比较是不合理的。中心区域的数据偏差可能是由于不同的滑移条件造成的。验证了本文模拟计算方法的准确性，即用此模型可以对水槽内附壁效应分析。一般情况下，无论是模拟还是实验结果，轴向速度的变化趋势都是相同的，这表明该模型在大多数情况下都能较准确地预测水槽中的气泡上升。

图2 模型验证 Fig.2 Model verification

4 结果与分析

气泡孔产生的大量气泡在上升中，气泡受壁面影响，上升方式各不相同。计算时间设置为气泡到达自由液面时刻，此时水槽计算区域内的流态已发展充分。以H = 6 m、2L = 2 m、d = 2 mm 和v = 25 m·s-1模型为例，分析随时间变化的速度矢量图，结果如图3(a)所示。从图中看出，气泡从气孔上升，沿着水槽的中心线直线上升，气泡流周围因卷吸而产生漩涡，由于气泡群的不均匀性导致在气泡流两侧形成的漩涡并不相同，从而气泡流两侧产生压力差，左右会蛇形摆动。随着时间变化，气泡流摆幅加剧，漩涡垂向流动，压力差增大，最终可能贴到壁上。取10 s 时刻，体积分率为0.5 的等值线图，即表示气泡的轮廓，见图3(b)。

图3 气泡时变特性速度矢量图及10 s 时刻气泡轮廓图 Fig.3 Time-varying profiles of bubble speed vector and bubble contour map at 10 s

气泡羽流的倾斜程度直接反映其受到附壁效应的强弱，倾斜程度可以用θ、η 作为指标共同说明。下面对影响倾斜程度的因素进行分析。

4.1 弗劳德数的影响

通入空气速度v 和孔口直径d 等参数对Fr 均有影响，取实验台尺寸2L=1 m、H=1 m 模型，对不同Fr 值下的倾斜程度分析。

4.1.1 速度场分析

图4 是不同Fr 下的速度流线图(水平和垂直分别为X 和Y 方向)。选取气泡流基本垂直上升、单侧偏移、幅度较小的蛇形摆动、幅度较大的蛇形摆动等4 组典型算例。从图4(a)中看出，当Fr 较小时，流域内气泡沿着孔口中心线基本垂直上升至自由液面，垂直中心线两侧速度基本对称且产生相反方向漩涡，附壁效应不明显；从图4(b)中看出，随着Fr 的增大，气泡流两侧形成高压和低压区。此时，气泡会向某一侧倾斜，气泡两侧漩涡强度不对称，气泡流偏移侧(即低压区侧)速度变化较剧烈；从图4(c)中看出，气泡发生蛇形摆动。由于摆动幅度小，流域内也仅有2 个漩涡，水槽内除摆动部分以外的速度分布较均匀；从图4(d)中看出，气泡摆动幅度较大，气泡偏移程度变大，形成多个强度不一的漩涡，水槽内速度分布紊乱，压力损耗变大。这是由于Fr 较大时，气泡流对原压力场冲击能力更强，造成的游涡数量更多，影响范围更广，强度更为剧烈。

图4 不同Fr 下流域内速度流线图 Fig.4 Velocity flow diagram within the field under different Fr

4.1.2 摆动特性分析

图5 所示为不同Fr 值下的偏移距离和倾斜角度。从图中可以看出，随着Fr 增大，偏移距离和倾斜角度分别呈现指数衰减、先增长后摆动趋势。这是由于Fr 较小时，气泡基本垂直上升，两者均很小。伴随Fr 增大，气泡蛇形摆动，漩涡强度及数量增加，偏移距离也随之增大。同时，水槽内摆动紊乱继续加大，偏移角度也变大。Fr 增至一定数值后，紊乱程度达到某种剧烈，而水槽空间有限，所以气泡运动剧烈使得提前达到偏移距离最大值并保持在固定值附近，偏移角度大小会随漩涡位置变化而变化。

图5 不同Fr 下的η 和θ 值(ε = 1) Fig.5 η and θ under different Fr (ε = 1)

4.2 水槽结构尺寸的影响

水槽模型中高度H 和宽度2L 等因素合称为结构参数，它们对附壁效应均有影响，构建ɛ 值具体分析水槽内附壁效应。取v = 15 m·s-1、d = 0.8 mm 工况，对不同ɛ 值下倾斜程度分析。

4.2.1 速度场分析

图6 为不同高宽比ɛ 下流域内速度流线图。随着ɛ 值不断增大，气泡流从基本垂直上升、单侧摆动、蛇形摆动到附壁。从图中看出，当ɛ = 0.8 时，速度关于垂直中心线基本对称；当ɛ = 0.86 时，流域内存在明显的高低压区，气泡两侧漩涡明显；当ɛ = 1.5 时，流域内存在多个漩涡，气泡流周围流场紊乱；当ɛ = 5 时，气泡在上升中发生附壁效应，在附壁区周围漩涡较大。

图6 不同高宽比下流域内速度流线图 Fig.6 Velocity flow diagram within the field under different aspect ratios

4.2.2 摆动特性分析

图7 为不同ɛ 下的η 和θ。随着ɛ 的不断增大，气泡会随着水槽区域水平距离减小而在流域内扩散，继而导致气泡倾斜距离变大。当ɛ 大到一定值，气泡羽流发生附壁效应，并开始贴着壁面上升。而倾斜角度θ 随着ɛ 增大先微量减小后增大最后减小，这是由于ɛ 不大时，气泡上升保持微弱单侧偏移；接着达到贴壁上升的ɛ 值时，气泡会在上升一小段高度之后开始贴壁，继而导致θ 急剧增大；ɛ 继续增大，气泡流在有限水槽内形成不了足够的高低压区，所以最后θ 会急剧减小。

图7 不同ɛ 下的η 和θ 值(Fr = 67.15) Fig.7 η and θ under different ɛ (Fr = 67.15)

4.3 无量纲数据处理

从康达效应产生的情况来看，气泡羽流的摆动特性与水域内流场的变化有关，构建无量纲的η、θ 和操作参数Fr、水槽结构参数ɛ 的关系式：

图8 Fr 与η 的对应关系 Fig.8 Relationship between Fr and η

将Fr 与η、θ一一对应，通过不同ɛ 构建回归方程，其变量数值如表2。图8 为指定ɛ 值时最大偏移量与 Fr 的对应关系，得到的拟合方程表达式为：

表2 不同ɛ 下的Fr-η 拟合方程各变量数值 Table 2 Values of variables of the Fr-η fitting equation under different ɛ

此直线的拟合程度分别为R2=0.860、0.926、0.829、0.845 和0.917。由此说明，η 和Fr 拟合度较好。随着ɛ 增大，使得η 等于1 时的Fr 拐点值减小。

图9 为指定ɛ 值时最大偏移量对应的θ 与Fr的对应关系，其各系数值如表3 所示，得到的拟合方程表达式为：

图9 Fr 与θ 的对应关系 Fig.9 Relationship between Fr and θ

此拟合直线的拟合程度R2=0.817、0.812、0.827、0.799、0.807。由此说明，θ 和Fr 拟合度较好。

表3 不同ɛ 下的Fr-θ 拟合方程各变量数值 Table 3 Values of variables of the Fr-θ fitting equation under different ɛ

图10 不同Fr 和ɛ 情况下气泡上升规律 Fig.10 Bubble rising characteristics under different Fr and ɛ

规定在流域内气泡流没有出现拐角且气泡偏移很微弱的上升模式为基本垂直上升，代号为1；没有出现拐角且气泡偏移明显或仅出现一个拐角的模式为单侧摆动，代号为2；出现两个及以上拐角的模式为蛇形摆动，代号为3；有部分区域气泡贴着壁面上升模式为附壁上升，代号为4。所有算例的上升方式统计如图10 所示。当ɛ=0.3 时，上升方式主要以2(单侧摆动)为主；当ɛ=1 时，上升方式以2(单侧摆动)、3(蛇形摆动)为主，1(基本垂直上升)为辅；当ɛ=1.5、2.5、5 时，上升方式以3(蛇形摆动)为主，4(附壁上升)为辅。