胡波，庞明军

（常州大学机械与轨道交通学院，江苏省绿色过程装备重点实验室，江苏常州213164）

非牛顿流体内气泡自由上浮运动广泛存在于自然界和工业生产过程中，如生物发酵、废水处理和生命维持系统等[1-3]。气泡的运动特性会直接影响液相的含气率，从而进一步影响气液相间的传热传质效率。因此，深入理解气泡在非牛顿流体中的运动特性，对于泡状设备的优化设计和气液相间传热传质效率的改善具有重要意义。

为了理解气泡在非牛顿流体内的运动特性（如运动轨迹、气泡形状、气泡上升速度和尾流等），国内外学者针对单气泡在静止非牛顿流体中的运动特性进行了大量的研究[4-10]。Premlata等[5]利用数值模拟和实验测试的方法研究了气泡在牛顿流体和剪切稀化流体中的上浮特性，发现气泡在牛顿流体中的上升轨迹为直线；而在剪切稀化流体中，气泡的运动轨迹与剪切稀化程度相关，剪切稀化程度越强，气泡的上升轨迹越不稳定，以之字形或螺旋形轨迹上浮。姜韶堃等[11]通过实验研究了气泡在羧甲基纤维素（carboxymethyl cellulose，CMC）溶液中的运动特性，分析了气泡振荡频率和黏度振荡频率，发现液相局部表观黏度的变化会对气泡振荡变形产生重要影响。Li等[1]利用粒子图像测速法和有限差分法研究了气泡在剪切稀化流体中周围液相表观黏度的分布特征，发现气泡周围黏度分布与流场直接相关，气泡尾部会形成类活塞流，活塞流边缘周围流体的剪切速率较大，而活塞流中心处流体的剪切速率较小，表观黏度呈现出与剪切速率相反的分布趋势。后来Li等[12]又利用level-set法数值研究了剪切稀化流体中气泡变形和流线分布特征，发现球形气泡尾部流线和内环流呈上下左右对称分布，而椭球形气泡尾部流线和内环流仅是左右对称、上下分布不对称，且气泡尾部的流线和内环流形状更细长些。当气泡形状变为球帽形时，其尾涡尺寸变得较大。Vahabi等[13]采用weakly-compressible smoothed particle hydrodynamics（WC-SPH）法研究了气泡在剪切稀化流体的上浮过程，发现随着流体剪切稀化程度的增强，气泡逐渐从球形变为下端带裙状的帽形气泡，在气泡周围出现了低黏度区域。张菊[14]通过实验也发现了同样的现象，她还指出剪切稀化特性使气泡上浮过程的阻力减小、速度增大。

Oshaghi等[15]采用VOF（volume of fluid）法研究了低毛细数（Ca）和低邦德数（Bo）下剪切稀化流体中气泡的生成和脱离过程并与牛顿流体进行了对比，发现剪切稀化流体和牛顿流体中气泡生成过程瞬时接触角的变化趋势相同；而气泡在剪切稀化流体中更易脱离，脱离时体积更小。Anjani等[16]利用level-set法研究了气泡在剪切稀化流体中的上浮过程，发现气泡尾部出现了高黏度区。除了对剪切稀化流体内气泡上浮运动特性的研究外，一些学者对剪切稠化和屈服应力流体中气泡的上浮运动特性也进行了相关研究。如Ohta等[17]采用level-set耦合VOF（CLSVOF）法，研究了气泡在剪切稠化流体中的上浮运动，发现随着流体剪切稠化程度的增强，气泡的变形和上浮速度减小。Sikorski等[18]实验发现屈服应力对小体积气泡的影响更强，更易出现较大的长径比。

综上所述，尽管国内外学者针对单气泡在非牛顿流体的运动特性做了一定的研究，但大都只研究了流体不同剪切稀化程度（即流变指数n）对气泡运动特性的影响。作为事实，特征时间λ是剪切稀化流体的一个重要物性参数，它对流体的表观黏度影响巨大，然而目前关于特征时间对气泡水动力学特性的影响缺乏系统深入的研究。从液相本构方程来看，特征时间对液相表观黏度的影响必然会对上浮气泡的形状、终端速度、气泡尾涡以及气泡周围液相表观黏度分布产生重要的影响。为此，本文运用VOF数值模拟方法详细研究了气泡在剪切稀化流体内的上浮运动特性，以便理解特征时间对气泡上浮运动特性的影响规律和影响机理。

1 物理模型和计算工况

1.1 几何模型

图1 计算区域

为了节约计算成本，将计算区域简化为一个二维矩形区域，如图1所示。让直径为d=6mm初始形状为圆形的气泡，在浮力的作用下从低部开始自由上浮。为了消除壁面对气泡上升运动的影响，同时保证气泡在竖直方向上有足够的空间能够达到稳定运动。计算区域的尺寸L和H设置为25d和40d。表1列出了文献中计算区域的尺寸。从表1可以看出，目前计算区域的宽度大于表1中所有文献的宽度；而在高度上，只有文献[8]的高度稍大于本文。这表明当前计算区域的尺寸能够满足计算的要求。为了消除底壁对气泡运动的影响，让气泡在距下壁面20d/3处开始释放。

表1 计算区域尺寸对比

1.2 控制方程

本次数值模拟是基于气液两相均为不可压缩流体，两相流均为层流且在恒温情况下进行的。气液两相流的控制方程如式(1)、式(2)所示。

式中，p为压强，Pa；u为速度矢量，m/s；Fs为表面张力所引起的体积力，N；μ(F)为流体局部平均黏度，Pa·s；ρ(F)为流体局部平均密度，kg/m3；F为计算域内的相函数；∇为哈密顿算子；D为应变率张量。应变率张量计算如式(3)所示。

流体局部平均密度ρ(F)和黏度μ(F)按式(4)、式(5)计算。

式中，下角标l和g分别表示液相和气相。

1.3 液相本构方程

采用Carreau模型来描述液相的黏度，如式(6)所示。Carreau模型的优点在于对低剪切速率和高剪切速率下的幂律行为都能给予准确描述。一旦计算区域中每个节点上的局部剪切速率已知，就可以得到流体的表观黏度。

式中，n、λ、μ∞、μ0为液相黏度特性参数，其大小取决于液相自身属性。n为流变指数，当n=1时，液相为牛顿流体；当n＜1时，为剪切稀化流体。根据文献[12,16]，λ为流体的特征时间。当保持稳态的液相突然受到扰动时，液相中的高分子链受到扰动、平衡态遭到破坏，而特征时间λ是描述液相从非平衡态恢复到平衡态所需要的时间。λ越大，表示流体从受到扰动到恢复平衡态所需的时间越长。当λ较大时，气泡变形非常不稳定，计算难以收敛，故目前计算最大λ设置为2s。当λ=0或n=1时，为牛顿流体。μ0和μ∞分别表示零剪切速率和无穷剪切速率下液相的表观黏度。

1.4 VOF方法

相比于其他界面捕捉方法，VOF法对于气液界面的捕捉具有计算精度高和计算速度快的优点[20]，为此运用VOF法捕捉气液界面。其相函数输运方程如式(7)所示。

式中，u为速度矢量；F为整个计算域的相函数，其值如式(8)所示。

计算出相函数之后，可利用式(9)计算出气相的局部体积分数。

式中，下角标s为网格编号；As为编号为s网格单元的面积。

基于上面的计算，利用分段线性法（PLIC）[21]进行界面重建。其原理是用直线线段近似表示单个网格内气液两相的界面，具有误差较小、精度较高[22]的优点。通过求解单个网格内气泡的局部体积分数梯度来获得界面的法向量。法向量计算如式(10)所示。

通过当前网格内气泡局部体积分数ψ和法向量n就可以获得气液界面的精确位置。在每个网格中，下一时刻的气泡局部体积分数是当前时刻网格内气泡局部体积分数与这段时间流入网格内气泡体积分数之和。通过更新下一时刻网格内的气泡局部体积分数，就可以对整个流场中的气泡界面进行重构[23]。

1.5 表面张力模型

计算时，采用Brackbill等[24]提出的连续表面张力模型来计算表面张力，并且将表面张力看作体积力加入到动量方程中的右边，考虑气泡加入对液相流动的影响。因此，动量方程中的源项Fs如式(11)所示。

式中，δ为狄拉克分布函数；k为界面曲率，可用单位法向量的散度表示［式(12)］。

1.6 数值方法

计算时，在如图1所示的AB边设置为压力出口边界条件，AC、BD和CD边均设置为无滑移边界条件。因为结构化网格生成速度快、质量好、易于实现计算区域的边界拟合且一般相对于非结构化网格来说，具有计算精度高、计算结果更容易收敛的优点。所以采用结构化网格划分计算区域。运用商业软件ANSYS18.0进行计算，在计算过程中，压力-速度场耦合使用带压力隐式的算子分裂（PISO）算法求解。压力和动量的离散格式分别为交错压力（PRESTO）和二阶迎风格式。时间步长经计算检验后选择为5×10-4s，当速度和动量残差小于10-6时，可视为计算取得收敛。计算时，首先初始化，使计算域充满静止液体，然后在其底部放入初始形状为圆形的气泡，进行两相流计算。

1.7 计算条件和计算工况

气泡在上浮过程中主要受到重力、黏性力、表面张力的影响，这三个力可用Ga（Gallilei）数和Eo（Eötvös）数来表征，如式(13)和式(14)所示。Ga数表示重力与黏性力的比值，Eo数表示重力与表面张力的比值。鉴于剪切稀化流体的表观黏度与特征时间λ紧密相关，为此研究了不同剪切稀化程度和表面张力下λ对气泡运动特性的影响。计算时，忽略重力水平的影响（即初始Ga数保持不变），在设计工况时，取Ga=3，λ取0.2、0.6、1.2和2，对应的n分别取0.2、0.4、0.6、0.8和1，Eo数分别取5、20、100和200。具体工况见表2。

式中，ρl为液相密度，kg/m3；g为重力加速度，m/s2；d为气泡直径，mm；μ0为液相初始黏度，Pa·s；σ为表面张力，N/m。

2 结果分析与讨论

为了对比分析，对所有计算工况的结果都进行了量纲为1处理。气泡不同时刻中心高度（Hg）用最大中心高度（Hg,max）进行量纲为1化；气泡不同时刻瞬时速度U用最大瞬时速度Umax进行量纲为1化；剪切速率γ˙用最大剪切速率γ˙max进行量纲为1化；液相表观黏度μl用液相零剪切黏度μ0进行量纲为1化，时间t用时间步长量纲为1化处理。

2.1 时间步长选取

数值计算时，时间步长的选取非常关键。若时间步长取得过大，会导致计算误差偏大、准确性较低，甚至会引起计算发散；若时间步长取得过小，则会大大增大计算量、耗时过长。为了选取合适的时间步长，选择在较密的网格尺寸下计算最苛刻的工况（网格独立性验证工况与此相同），分别对比了0.002s（Step1）、0.001s（Step2）、0.0005s（Step3）和0.00025s（Step4）4个时间步长气泡的瞬时高度和瞬时速度，计算结果如图2和图3所示。从图中可看出，随着时间步长的减小，气泡上浮过程中瞬时高度（和瞬时速度）不再发生变化。0.002s和0.001s与0.0005s气泡的瞬时高度（和瞬时速度）差距明显，而0.0005s和0.00025s气泡的瞬时高度（和瞬时速度）基本趋于一致。为了减小计算时长，又同时保证计算结果的准确性，时间步长选为0.0005s。

表2 计算工况

图2 不同时间步长下气泡的中心高度

图3 不同时间步长下的气泡瞬时速度

2.2 网格独立性验证

为了保证网格对计算结果不产生影响，建立了Grid1、Grid2、Grid3和Grid4四套网格系统，对应的网格尺寸分别为0.40mm、0.30mm、0.20mm和0.15mm。在验证网格独立性的时候，选择最苛刻的工况进行验证（此时气泡的变形程度最大），即Ga=3，Eo=200，n=0.2和λ=2。气液相密度比ρr=ρg/ρl=0.001，黏度比为μr=μg/μ0=0.01。液相零剪切黏度μ0=0.49Pa·s，液相无穷剪切黏度μl=0.0098Pa·s。计算结果如图4和图5所示。

图4 不同网格尺寸下气泡的中心高度

图5 不同网格尺寸下气泡的瞬时速度

从图4和图5可以看出，随着网格的不断加密，计算结果不再发生变化。其中，格子Grid1与Grid2的计算结果存在一定的偏差，而格子Grid3与Grid4的结果就基本相同。为了减小计算量且能获得较为精确的计算结果，最终选取网格系统Grid3进行计算。

2.3 结果准确性检验

为了验证计算结果的准确性，首先将气泡的终端速度（UT,cal）与文献[25]中的实验值（UT.exp）进行了对比，计算时表面张力系数设置为0.072N/m，其余参数见表3。为了对比，表3中的计算工况与文献[25]中的实验工况一致。对比表明，目前计算的终端速度与文献[25]中的实验结果相差约10%。其次，按照文献[4]的实验工况设置了计算工况，对比了气泡的形状和尾涡结构，如图6所示。从图中可以看出气泡的形状和尾涡结构与文献[4]中的实验结果高度相似。可见，目前的计算方法能够满足计算要求，微小误差存在的原因可能是因为目前计算采用的是二维模型。

2.4 气泡形状分析

气泡的形状会直接影响气泡的水动力特性，为此首先分析了不同剪切稀化程度和表面张力下，特征时间λ对气泡形状的影响，如图7和图8所示。研究时假定重力保持不变，这样表面张力的影响可用Eo数来体现。

表3 气泡终端速度对比

图6 Ga=116和Eo=94.3时气泡形状和尾涡结构

从图7可以看出，当Eo数相同时（对于所有的Eo数），流体的剪切稀化程度不同（n不同），λ的变化对气泡变形的影响不同。当流体剪切稀化程度较弱（n=0.8）时，气泡的横向尺寸随λ的增大稍微增大，厚度（气泡顶部到底部的距离）基本不变。这是因为此工况下的流体与牛顿流体相似，稳定性较好，气泡自由上浮运动对周围流体的扰动较弱，所以λ的变化对气泡变形的影响较弱。但当流体的剪切稀化程度逐渐增强（n＜0.8）时，随着λ的增大，气泡表现出厚度明显减小、横向尺寸显著增大的变形趋势，且球帽形气泡逐渐出现裙状，显现出不断向下生长的趋势，发展成带裙状（裙状脱离和未脱落）的帽形气泡。当n相同时，对于不同的Eo数，λ的增大也会引起气泡明显的变形。对于较大的n（n≥0.6），当Eo数较小（如Eo=5）时，气泡的形状随着λ的增大变形成椭球形；当Eo数较大（如Eo=100、200）时，气泡的形状随λ的增大变形成球帽形或底部两端带裙状的帽形气泡。值得注意的是，当n较小时（n≤0.4），气泡的形状比较特殊。当Eo数（Eo=5）较小时，气泡的形状在小的λ（λ=0.2）下保持为椭球形；随着λ的增大，气泡的扁平程度逐渐增大，变成了顶部中央下凹、两端向上凸起的椭球形气泡；且λ越大，这种变形越明显。当Eo数增大至20时，在小的λ（λ=0.2）下，气泡变形成底部水平的帽形气泡；随λ的增大，气泡形状呈现出与Eo=5相反的趋势，气泡最终变成了顶部上凸、两端水平的帽形气泡。当Eo数进一步增大至100或200时，随着λ的逐渐增大，气泡裙状不断向下生长，形成底部两端带裙状的帽形气泡；当λ增大到一定程度（λ≥1.2）时，气泡裙状出现了与气泡脱离的现象；λ越大，这种现象越明显，气泡变形越剧烈。总之，在剪切稀化程度较弱（n=0.8）的流体中，λ的增大对气泡的变形影响不大，但当流体的剪切稀化程度增强（n＜0.8）时，λ对气泡形状的影响随之增大。小Eo数下，λ的增大降低了气泡的厚度，使气泡呈现为椭球形；在大Eo数下，λ的增大不仅降低了气泡的厚度，还出现了裙状，且裙状尺寸不断增大直至脱落。除此之外，当λ相同时，在较小的Eo（如Eo=5）数下，随着n的减小，气泡的厚度逐渐减小，扁平程度增大；而在较大的Eo数（如Eo=100、200）下，n的减小使得气泡尾部两端的裙状不断增大甚至出现脱落现象。这与文献[7-8,26]报道的结果相似，详细内容可参阅文献，不再赘述。

图7 剪切稀化流体内不同工况气泡变形

特征时间λ是指液相从非平衡态恢复到平衡态所需要的时间。λ越大，液相受到扰动后恢复到平衡态所需的时间越长。从微观结构来看，λ的大小与液相内部溶解的分子链的柔顺性有关，λ越大，柔顺性越差，流体达到平衡态所需要的时间越长。对应的计算过程收敛所需要的时间越长。从液相的本构方程也可以看出，在给定的n、液相零剪切黏度（μ0）和无穷剪切黏度（μ∞）下，λ越大，液相的表观黏度越低，气泡的Ga越大。Ga越大，气泡形状越不稳定、容易发生变形。另外，根据高分子结构构象理论[27]，气泡周围液相分子链由于受气泡上浮运动剪切的影响，由初始的蜷曲状被迫沿着流动方向伸展取向，使液相的表观黏度减小，表现出剪切稀化效应[28]。λ越大，分子链在外力作用下达到稳定变形所需的时间越长，剪切稀化效应越强，所以气泡变形也越明显。

2.5 气泡尾涡特征

气泡形状的变化会对气泡局部水动力学特性产生重要的影响，特别是气泡的尾涡会随气泡形状不断变化，最终达到稳定状态。为了清楚地展示在不同剪切稀化程度和表面张力下特征时间λ对气泡尾涡的影响，图9给出了部分剪切稀化流体不同Eo数下λ对气泡尾涡的影响，图中红色轮廓线为气泡的形状。为了减小篇幅，仅给出λ=0.2和λ=2两种情况的尾涡图。

从图9可以看出，当液相的剪切稀化程度较强（n=0.2）时，对于每个Eo数，由于λ对气泡形状的影响导致其对气泡尾涡也产生了明显的影响。当Eo数等于5或20时，由于气泡小的变形，λ=0.2时，气泡尾部未出现尾涡；但当λ=2时，由于气泡大的变形，在其尾部出现了一对尺寸较大且形状对称的尾涡。当Eo等于100或200时，气泡发生了较复杂的变形，λ=0.2时气泡尾部出现了一对对称分布的尾涡；当λ=2时，气泡尾部出现了一小一大两对对称的尾涡，且尾涡结构变化较为复杂。当液相的剪切稀化程度较弱（n=0.8）时，小的Eo数下，对于所有的λ气泡尾部均未出现尾涡；大的Eo数下，对于所有的λ气泡尾部均出现一对对称分布的尾涡，且尾涡的尺寸随λ的增大而增大。而当λ相同时，随着Eo数的增大，在剪切稀化较弱的（n=0.8）流体内气泡尾涡结构的变化趋势与牛顿流体（n=1或λ=0）中气泡尾涡的变化趋势高度相似；但在剪切稀化程度较强（n=0.2）的流体中，气泡的变形加剧，其尾涡尺寸明显变大。此现象与Premlata等[8]和倪明玖[26]报道的结果相似。这说明当液相的剪切稀化较弱且Eo数较小时，λ对气泡形状和尾涡影响不大；当流体剪切稀化程度较强或Eo数较大时，λ对气泡形状和尾涡影响非常明显。

为了进一步探究λ与气泡尾涡之间的关系，图10给出了Eo=200和n=0.6工况下λ对气泡尾涡特征的影响。从图中可以看出，随着λ的增大，气泡尾涡尺寸不断增大。当λ增大至1.2时，气泡尾涡由一对演变为两对；当λ继续增至2时，由于气泡底部两边裙状向下增长，使得气泡尾部两对对称尾涡的尺寸明显增大。可见，即使剪切稀化程度和表面张力相同，λ的改变也对气泡尾涡有着重要的影响。

2.6 气泡终端速度

图9 t=0.75s时气泡不同工况下的尾涡图

气泡的终端速度也是表征气泡局部水动力特性的一个重要参数，它会影响液相的含气率，而含气率的高低又会对两相间的传质传热效率产生重要影响。为了研究λ对气泡终端速度的影响，图11给出不同n和Eo数下气泡终端速度随λ的变化趋势。如图所示，对于相同的n和Eo数，气泡的终端速度随着λ的增大逐渐增大。在λ=0.2～0.6s区间内，气泡终端速度增幅较大；而在λ=0.6～2s区间内，气泡终端速度增幅较小。当n≤0.4或Eo数较小（Eo=5）时，气泡终端速度随着λ的增大，其增幅要比n≥0.6或Eo数（Eo＞5）较大的工况明显。这表明，对于剪切稀化程度较强或表面张力较大的流体，气泡的终端速度更容易受λ的影响。此外，在相同的λ下，对于给定的n，气泡的终端速度随着Eo数的增大不断减小；而对于给定的Eo数，随着n的减小，气泡的终端速度不断增大。这与文献[8,26]得到的不同剪切稀化程度和不同Eo数对气泡终端速度的影响相似。气泡终端速度与气泡上浮达到稳态运动的受力有关。当气泡运动达到稳态运动时，气泡主要受浮力和阻力的作用。当气泡的体积相同时，气泡所受的浮力是相同的。因此，气泡终端速度不同，主要是因其阻力不同而导致的。而气泡所受的阻力包括形状阻力和摩擦阻力。形状阻力的大小取决于气泡形状对尾涡的影响，尾涡强度和影响区越大，气泡受到的形状阻力越大。摩擦阻力与液相的黏度有关，当n越小和λ越大时，相同的剪切速率下，气泡周围液相的表观黏度越小，气泡上浮过程受到的摩擦阻力则越小。在两种阻力的联合作用影响下，导致了目前气泡终端速度的不同。为了理解气泡所受摩擦阻力的大小，下面分析了气泡周围剪切速率和表观黏度的分布情况。

图11 气泡终端速度

2.7 气泡周围液相剪切速率

对照图11，图12给出不同n和Eo=5、Eo=200时气泡周围液相剪切速率随λ的变化云图。从图中可以看出，气泡顶部和底部两翼均出现了高剪切速率区，这与Battistella等[29]的报道相似。从图12可以看出，当Eo数相同时，对于每个n，气泡周围液相高剪切速率区的面积都随λ的增大而增大，且气泡尾翼两侧高剪切速率区逐渐被拉长；而对于不同的λ，随着n的减小，气泡周围液相高剪切速率区也会增大。这表明，在相同的Eo数下，增大λ（n不变）与减小n（λ不变）对剪切速率的影响相似，都会增大气泡周围液相高剪切速率区的范围。当n相同时，λ对气泡周围液相高剪切速率区的影响与Eo数有关。当Eo=5时，对于相同的n，随λ的增大，气泡周围高剪切速率区的影响范围明显比Eo=200的大，气泡尾翼两侧条状高剪切速率区呈竖直分布。

图12 不同工况下气泡周围液相的剪切速率分布云图

气泡周围液相剪切速率分布的不同与气泡周围液相的速度梯度有关，速度梯度大的地方，液相剪切速率越高，反之亦然。液相速度梯度分布受气泡终端速度和尾涡的影响。气泡的终端速度越大，气泡周围液相高剪切速率区越大，如图11和图12所示。对于Eo=5和n=0.2的工况，λ从0.2增至2时，以及Eo=200和λ=0.6的工况，n从0.8减小至0.2时，气泡终端速度的增大导致气泡周围液相高剪切速率区的增大。也就是说，气泡终端速度越大，其引起周围流体的速度也越大，气泡周围大的液相速度与远离气泡处静止的液相产生了较大的速度梯度，在气泡周围形成大的高剪切速率区。另一方面，在气泡尾涡的两侧，流线越密、速度越大，也导致此部分与远离气泡处静止的流体产生了较大的速度梯度，因而在气泡尾翼形成两长条状高剪切速率区。在相同的n下，对于Eo=200的工况，随λ的增大，气泡尾涡横向尺寸逐渐变宽，使得气泡尾翼两侧长条状高剪切速率区分布比相应条件下Eo=5的工况中更宽。

2.8 液相表观黏度分析

对于剪切稀化流体，液相的速度场会影响表观黏度的分布，而液相的表观黏度又会影响气泡的水动力学特性，因此有必要分析液相表观黏度的分布情况。对照图12和图13可以看出，液相表观黏度的分布趋势与剪切速率的分布趋势恰好相反，即低剪切速率区对应高表观黏度区、高剪切速率区对应低表观黏度区。气泡周围液相低表观黏度区的增大会减小气泡上浮过程所受到的摩擦阻力，使气泡的终端速度增大。这也部分解释了图11气泡终端速度不同的原因。另外，从图13中也发现，当Eo=200时，在气泡尾部出现了黏度盲区，即局部高黏度区域，该区域内流体的表观黏度几乎未受到剪切稀化的影响。文献[13,25]等报道过类似的现象。必须指出的是，目前关于黏度盲区的报道基本上来源于数值模拟结果，尚未见过相关结果的实验报道。在Eo=200的工况中，随着λ的增大，黏度盲区逐渐脱离气泡尾部并发生破碎现象。而在Eo=5的工况中几乎未出现这种现象。由于黏度盲区的存在，导致Eo=200工况气泡周围液相低表观黏度区始终小于对应条件下Eo=5的工况。所以，对于所有Eo=200的工况，气泡上浮过程所受到的摩擦阻力始终大于对应条件下Eo=5的工况，导致前者气泡的终端速度始终小于后者的（图11）。

图13 不同工况下气泡周围液相表观黏度分布云图

3 结论

鉴于理解气泡在非牛顿流体运动过程对工农业生产的重要性，运用VOF数值模拟方法，详细分析了不同剪切稀化程度和表面张力下特征时间对气泡水动力特性的影响。目前发现：在相同的特征时间下，流变指数和表面张力对气泡水动力特性的影响与文献[6-7,11,14,26]报道相似。除此之外，关于特征时间对气泡水动力特性的影响，主要得出以下几点结论。

（1）气泡的水动力学特性与液相的特征时间密切相关，液相的剪切稀化程度越强（流变指数n越小）或液相的表面张力越小（Eo数越大），特征时间的变化对气泡变形和尾涡的影响则越大。

（2）在同一剪切稀化程度和同一表面张力水平下，随特征时间的增大，气泡终端速度的增幅不同；在0.2～0.6s区间内，增幅较大，而在0.6～2s区间内，增幅较小；而在剪切稀化程度较强或表面张力较大的液相中，气泡终端速度易受特征时间的影响。

（3）在同一剪切稀化程度和同一表面张力水平下，特征时间越大，气泡的终端速度越大，尾涡的强度越大、范围越宽，导致气泡周围高剪切速率区和低表观黏度区的范围也越宽。

（4）对于低表面张力的工况（Eo=200），在气泡尾部出现了黏度盲区，且随特征时间的增大，黏度盲区逐渐脱离气泡尾部并破碎；黏度盲区的出现，减小了气泡周围液相低表观黏度区的面积，增大气泡上浮过程的摩擦阻力，降低了气泡终端的速度。