在自然界和工农业生产领域中,经常会看到气泡在液相中的上浮运动.鉴于气泡的运动特性会影响液相的流变特性以及气液相间的质量和热量交换,所以有必要对液相内气泡的运动规律开展深入的研究.现有的研究表明,气泡的运动特性除了受气泡大小、气液相热物理参数(黏度与表面张力系数)、重力水平和气泡界面洁净程度影响之外,还受气泡运动空间大小的影响[1-20].
通过对不同黏度Newton流体内气泡上浮运动特性的研究发现,气泡中心距其运动空间内壁距离的大小(即壁面效应)对气泡的变形、终端速度、尾涡结构和所受作用力都具有重要的影响[1-6].文献[1]指出,气泡的终端速度与壁面效应和表面张力系数密切相关;当气泡运动空间壁面法向尺寸大于10倍(优先15倍)气泡的当量直径时,壁面效应对气泡终端速度的影响可以忽略.文献[2]指出,当气泡运动空间尺寸大于8倍气泡的当量直径时,壁面效应即可以忽略.文献[3]指出,气泡的变形和上浮速度同时受气泡Reynolds数和Eötvös数、液相Morton数以及壁面效应控制;当气泡Reynolds数小于0.1时,气泡变形仅与壁面效应有关.文献[5]指出,当忽略重力时,壁面效应对气泡的热毛细运动也有重要的影响;当空间法向尺寸与气泡直径比大于等于3.75时,壁面效应对气泡的热毛细迁移速度的影响可以忽略.
目前的文献仅研究了壁面效应对Newton流体内气泡上浮运动特性的影响,而有关壁面效应对非Newton流体内气泡上浮运动特性的影响鲜有报道.一些学者仅研究了大空间环境下流变特性和气液相热物性参数对非Newton流体内气泡上浮运动特性的影响[7-12].为了洞察壁面效应对非Newton流体内气泡上浮运动特性的影响,采用数值方法深入研究了壁面效应对剪切稀化流体内气泡运动特性的影响.通过调节气泡直径与计算域宽度的比值,研究壁面效应对气泡形状、终端速度等参数的影响规律.
气泡Reynolds数对气泡的运动特性有重要的影响.当其大于200时,气泡的形状和尾流会呈现出三维特性;当其小于200时,仅呈现出二维特性[13].为了简化研究,只考虑气泡Reynolds数小于200的工况.本文计算气泡的Reynolds数约为20,所以几何模型可以使用二维轴对称模型,见图1.在图中,r和z分别表示坐标轴的径向和轴向;边界AB设置为压力出口,边界AD是对称轴(即轴对称边界条件),边界BC和CD是无滑移壁面.重力方向沿着z轴的负方向.BC的轴向高度H= 40d,以确保气泡的运动充分发展.初始时刻,气泡中心距计算域底部h=6d.d为气泡的初始直径, mm.
图1 几何模型
Fig.1 Geometric model
气泡上浮运动主要受重力、黏性力和表面张力的影响,这三个力对气泡运动的影响可用无量纲参数Gallilei数Ga和Eötvös数Eo来表征,其定义见式(1)和(2).Ga表示重力与黏性力的比值,Eo表示重力与表面张力的比值.壁面效应用计算域宽度与气泡直径的比值L/d表征,剪切稀化程度用流变指数n表征.文献[10]指出,当Ga较小或者Eo较大时,改变流变指数,气泡形状和终端速度的变化比较明显.所以选取较小的Ga和较大的Eo进行计算,即Ga=3和Eo=200.为了研究不同剪切稀化流体内壁面效应对气泡运动的影响,取n=1,0.8,0.6,0.4,0.2,L/d=0.75,1,1.5,2.5,3.5,5,6.5,8,10,12.5.
(1)
(2)
式中,ρl为液相密度,kg/m3;g为重力加速度,m/s2;μ0为液相初始黏度,Pa·s,σ为表面张力系数,N/m.
目前的计算是基于忽略流体的可压缩性和温度的变化、视流动为层流的假设开展的.为此,气液两相的控制方程可表示为
·u=0,
(3)
p+·[2μ(f)D]+Fs+ρ(f)g,
(4)
式中,p是压力,Pa;Fs是表面张力产生的体积力,N;ρ(f)是局部平均密度,kg/m3;u是速度,m/s;是Hamilton算子;μ(f)是局部平均黏度,Pa·s;D是应变率张量,即
(5)
流体局部平均密度ρ(f)和黏度μ(f)由式(6)和(7)计算:
ρ(f)=ρl(f)+ρg(1-f),
(6)
μ(f)=μl(f)+μg(1-f),
(7)
式中,下标l和g分别表示液相和气相,下文符号与此处相同.
因Carreau本构模型适用范围广,便于描述非Newton流体零剪切速率和无限剪切速率下的黏度.为此,选用Carreau模型作为剪切稀化流体的本构方程,即
(8)
在式(8)中,λ为时间常数;当λ=0时,流体为Newton流体.目前主要研究中低剪切速率下流体的特性,所以λ不宜过低,取λ=1.2.n是流变指数;当n=1时,流体为Newton流体.μ∞和μ0分别为无穷大剪切速率和零剪切速率下液相的黏度.计算表明,对于n=0.2的工况,由于剪切稀化效应较强,如果无穷大剪切速率黏度过低,会导致计算难于收敛,所以取μ0/μ∞=50.
要理解气泡在剪切稀化流体内的上浮运动特性,应准确捕捉气液相间的运动界面.鉴于VOF法具有计算精度高和计算速度较快的优点[10],故用来捕捉气相间的运动界面.VOF法的核心是相函数输运方程,即
f=0,
(9)
式中,f为整个计算域的相函数,其值定义为
(10)
则可利用式(11)计算出气泡的局部体积分数:
(11)
式中,下标s是网格的编号,As是编号为s网格单元的面积.目前的VOF法使用分段线性(PLIC)法来重构界面[14].该方法使用直线段(或平面)来近似某个网格内气液相的界面.界面法向量通过求解网格单元内气泡局部体积分数的梯度获得[14],法向量计算公式如下:
(12)
使用Brackbill等[15]所提出的连续表面张力模型来计算表面张力.在计算过程中,表面张力作为源项体积力Fs出现在动量方程中,即
(13)
式中,κ是相界面曲率,计算式为
κ=-(
(14)
式中,为单位法向量,
鉴于结构网格能够很好地拟合当前计算区域的边界,而且有利于相界面的重构,所以采用结构化网格生成网格系统.为了获得准确的计算结果和加快收敛速度,使用了自适应网格技术(AMR).计算过程中,采用h-adaptive自适应网格技术对网格进行重构.对于当前的研究,气泡仅沿z轴正向上升,在径向上不发生偏移.因此,采用区域自适应网格法对气泡运动区进行细化,网格系统如图2所示.当前的网格系统,有三个网格尺寸级别.首先,对整个计算域进行网格划分,网格尺寸为Δ1.接着用区域自适应方法对宽度为气泡直径3倍、高度与计算区域相同区域内的网格进行细化,网格尺寸为Δ2.为了准确捕捉气泡尾部的流动特征,再次采用相同的方法对宽度为气泡直径1.5倍、高度与计算区域相同区域内的网格进行细化,网格尺寸为Δ3.三个级别网格尺寸的关系为Δ1=2Δ2=4Δ3.
图2 网格系统
Fig.2 The mesh system
计算时先初始化,使计算域内充满静止的液体;然后放入初始形状为圆形的气泡,进行两相流动的计算.控制方程采用与时间相关的层流模型求解,速度和压力场的耦合采用PISO算法.采用一阶迎风格式和体积加权法对动量和压力项进行离散.动量、压力和体积力的亚松弛因子分别设为0.3、0.5和0.7,以保证较好的收敛性和稳定性;时间步长设置为5×10-5 s.当速度和动量残差低于1×10-6时认为计算已经收敛.
为了便于计算结果的分析和对比,以及消除初始参数对结果的影响,下面给出的结果都用相关参数进行了无量纲处理.无量纲化后的参数为:气泡终端速度U=Uf/Umax,气泡瞬时速度V=Vf/Vmax,剪切速率液相表观黏度μde=μf/μ0.其中,Uf为气泡终端速度,Umax为气泡最大终端速度,Vf为气泡瞬时速度,Vmax为气泡最大瞬时速度,为剪切速率,为最大剪切速率,μf为液相表观黏度.
为了避免网格系统对计算结果的影响,同时尽可能减小计算量,设计了四套不同的网格系统(即grid 1、grid 2、grid 3和grid 4),检查了网格尺度对计算结果的影响,其最小网格尺寸分别为d/40,d/60,d/80和d/100.为了保证计算精度,以最苛刻的工况,即壁面效应和流变指数最小的工况来验证网格无关性.计算参数为Ga=3,Eo=200,n=0.2,L/d=12.5,μr=0.01和ρr=0.001.图3为四套网格下气泡瞬时速度随计算时间的变化.
图3 四套网格下气泡瞬时速度
Fig.3 Instantaneous bubble velocities for 4 grid systems
气泡瞬时速度由下式计算:
V=(zt1-zt2)/Δt,
(15)
式中,zt1和zt2是t1和t2时刻下气泡中心位置,Δt为时间间隔.从图3可以看出,由网格系统grid 3和grid 4得出的计算结果高度吻合,可见继续加密网格不会影响计算结果.为了同时保证计算精度和减少计算量,选用最小网格尺寸为d/80的网格系统(grid 3)进行计算.
为了验证计算结果的准确性,将本文的计算结果与文献的数值[10,16]和实验[17-18]结果进行了对比,分别对比了气泡的形状、尾涡以及阻力系数,见图4~6.图6纵坐标为气泡的阻力系数,横坐标为表征流体黏弹性的无量纲参数,其计算方法见文献[16].可见,目前的数值结果和文献的实验和数值结果吻合得较好,所以目前的数值方法是可信的.
(a)文献[18]的实验结果 (b)文献[17](实线)和当前计算结果(虚线)
(a)The experimental result of ref.[18] (b)The numerical results of ref.[17](solid line)and the present computation(dash line)
图4 气泡形状对比(Ga=6.55,Eo=116和n=1)
Fig.4 Comparison of bubble shapes for Ga=6.55, Eo=116 and n=1
(a)左边是文献[10]的数值结果,右边是文献[18]的实验结果 (b)当前计算结果
(a)The numerical result of ref.[10]on the left side and (b)The present numerical resultthe experimental result of ref.[18]on the right side
图5 气泡形状和尾涡对比(Ga=93.9,Eo=44和n=1)
Fig.5 Comparisons of streamline patterns and bubble shapes for Ga=93.9, Eo=44 and n=1
图6 文献[16]实验与目前计算结果的阻力系数对比
Fig.6 Comparison of drag coefficients between the present numerical result and the experimental result of ref.[16]
要理解壁面效应对气泡上浮运动特性的影响,首先需明确其对气泡形状的影响.图7给出了不同流变指数下,壁面效应对气泡形状的影响.从图7可看出,对于所有的流变指数,随着壁面效应的增强(L/d的减小),气泡的纵向尺寸变大,而裙边长度减小.对于剪切稀化效应较强的流体,当L/d减小到0.75时,气泡的裙边几乎消失.对于剪切稀化较弱的流体,如n=0.8和1时,因气泡本身裙边很短,所以壁面效应对气泡裙边的影响不大.当L/d=2.5时,对于所有的流变指数,气泡形状都已基本稳定,继续增大L/d,气泡的形状不再发生明显的变化.对于所有的流变指数,气泡取得稳定变形所需的L/d近似相同.对于相同的L/d,剪切稀化效应对气泡变形的影响却很大.总之,强的壁面效应严重限制了气泡的横向变形,在相同的运动空间中,剪切稀化效应越强,气泡横向变形越大.
图7 壁面效应对气泡形状的影响
Fig.7 The wall effects on the bubble shapes
壁面效应对气泡形状的影响源于其对气泡变形过程的影响.图8给出了流变指数n=0.2时,壁面效应对气泡变形过程的影响.从图中可看出,相同时刻下,不同壁面效应所对应的气泡变形程度明显不同.当L/d=0.75时,受强壁面效应的影响,气泡在各个时刻下的总体变形均较小.随着壁面效应的减弱,不同时刻下气泡的变形明显增大,特别是在t=0.05时刻气泡变形显著增加.目前的结果与文献[19]给出的结果基本一致.可见,不同壁面效应下,气泡上升时的变形过程不同,导致稳态下气泡的形状不同.
气泡周围的流场与气泡的形状密切相关,而且它会对气泡的形状产生重要影响.为此,在研究气泡形状的基础上,下面分析了气泡周围液相的流场特征.图9给出了壁面效应对气泡周围流场的影响.图中黑色区域为气泡所占区域,左半部分为液相表观黏度分布,右半部分为液相剪切速率和流线分布.由图7可知,当L/d≥2.5时,气泡的形状基本稳定.为此,图9仅给出L/d≤3.5的结果.
图8 壁面效应对气泡变形过程的影响(n=0.2)
Fig.8 The wall effects on the evolution of bubble shapes for n=0.2
图9 壁面效应对气泡周围液相表观黏度、剪切速率和流线分布的影响
Fig.9 The wall effects on the apparent viscosity, the shear rate and the streamline pattern of liquid around the bubble
注 为了解释图中的颜色,读者可以参考本文的电子网页版本,后同.
从图9可以看出,流变指数(剪切稀化)和壁面效应对气泡周围的流场均有显著影响.对于剪切稀化较强(n=0.2和0.4)的流体,由于壁面效应对气泡变形的影响十分显著,导致气泡周围流场也受到明显的影响.随着壁面效应的减弱,气泡横向变形加剧、裙边变长,导致气泡尾涡逐渐形成、且尾涡尺寸逐渐增大.对于剪切稀化较弱(n>0.4)的流体,因气泡横向变形较小,气泡尾部虽无尾涡形成,但壁面效应对气泡周围流场的影响仍非常显著.从气泡周围的流线分布可看出,壁面效应对气泡周围液相流线的疏密有重要影响,也就是说气泡周围液相的速度大小受到重要影响.
流变指数(剪切稀化)和壁面效应对液相周围流场的影响,影响了气泡周围液相剪切速率的分布,进而对液相表观黏度分布也产生了重要影响.从图9中还可以看出,强的壁面效应导致流线被压缩,在气泡与壁面之间形成高剪切速率区.受气泡尾流的影响,在气泡尾部也出现了高剪切速率区,导致气泡周围液相的表观黏度出现了大范围的下降.除此之外,受气泡周围复杂流场的影响,在径向气泡边缘出现一个中等黏度区域,在气泡尾部还出现了高黏度区域.随着流变指数的增大(剪切稀化效应的减小),气泡尾部高黏度区域逐渐消失,但是径向气泡边缘中等黏度区域仍然存在.
随着壁面效应的减弱(L/d的增大,即L/d增大到1时),由于气泡变形程度的增大,导致气泡尾部高黏度区域随之增大.同时,由于液相剪切速率分布的改变,导致径向气泡边缘中等黏度区域的面积略微增大、并向外扩展.从剪切速率分布也可以看出,壁面附近没有出现高剪切速率区域,这说明此时壁面效应对气泡周围流场的影响已经减小.随着壁面效应的进一步减小,气泡周围液相表观黏度和剪切速率的分布不再发生变化.
总而言之,对于气泡周围液相的剪切速率分布而言,对于n=1的Newton流体,L/d=1.5时,剪切速率分布基本稳定;对于n=0.8,0.6,0.4和0.2的剪切稀化流体,L/d=2.5时,剪切速率分布基本稳定. 对于表观黏度分布而言,n=1时,液相表观黏度为定值;n=0.8,0.6,0.4和0.2时,当L/d=2.5时,液相表观黏度分布基本稳定.可见,对于n=1的Newton流体,L/d=1.5时,液相流场基本稳定;对于n=0.8,0.6,0.4和0.2的剪切稀化流体,L/d=2.5时,液相流场基本稳定.
图10 壁面效应对气泡终端速度的影响
Fig.10 The wall effects on the bubble terminal velocity
稳定后气泡的上浮速度是表征气泡运动特征的另一个重要参数.壁面效应对气泡形状及其周围流场的影响势必对气泡的终端速度也会产生重要影响.之前的分析表明,当L/d=3.5时,气泡的形状和其周围的流场基本稳定,不再随L/d的增大而变化.为了检验壁面效应对气泡终端速度的影响是否也符合这一规律,图10给出了不同流变指数下,壁面效应对气泡终端速度的影响.
图11 不同壁面效应下的压力分布
Fig.11 Pressure distributions under different wall effects
从图10可以看出,对于所有的流变指数,随着壁面效应的减弱(即L/d的增加),气泡的终端速度逐渐增大.可见,小的壁面距离极大地阻碍了气泡的上浮运动.对于不同的流变指数,壁面效应对气泡终端速度的影响也不尽相同.对于流变指数n=1和0.8的流体,当L/d=6.5时,气泡的终端速度趋于稳定.对于n=0.6,0.4和0.2的流体,当L/d=10时,气泡的终端速度才趋于稳定.可见,对于剪切稀化程度强的流体,气泡的上浮运动更易受壁面效应的影响.
上面的分析表明,对于所有的流变指数,气泡终端速度均随壁面效应的减弱而增大.除了壁面对气泡上浮运动的束缚,气泡上下两侧的压差或许对其也有重要的影响.为此,图11给出了气泡周围的压力分布.从图中可以看出,对于所有的流变指数,随壁面效应的减弱,气泡上下两侧的压差在减小,从而导致气泡终端速度增大.对于相同的壁面效应,随流体剪切稀化效应的减弱(即流变指数的增大),气泡上下两侧的压差增大,从而导致气泡终端速度减小.这或许是剪切稀化程度弱的流体,气泡终端速度低的原因之一.
另外,气泡终端速度的大小也受流体自身性质的影响.对于剪切稀化程度弱的流体,由于气泡周围液相的黏度变化不大(如对于n=1的流体,气泡周围液相的黏度始终等于液相的初始黏度),所以气泡所受的黏性阻力始终很大,因此其终端速度也最小.而对于剪切稀化程度强的流体,如n=0.2时,由于流动的剪切作用,气泡周围液相黏度大幅减小,所以气泡所受的黏性阻力也大幅减小,因此其终端速度也最大.
总的来说,气泡终端速度不仅与气泡形状、其周围流场分布有关,而且还受流体性质的影响.其可以综合反映各个参数对气泡上浮运动特性造成的影响.与气泡形状和其周围流场分布相比而言,壁面效应对气泡终端速度的影响最大,其能更直观地反映壁面效应对气泡上浮运动特性的影响. 所以,要消除壁面距离对气泡运动特性的影响,应以气泡终端速度不变作为判断几何模型大小的依据.
本文通过分析不同流变指数下,壁面效应对气泡运动特性的影响,可得出以下几点结论:
1)壁面效应对气泡形状的影响贯穿在气泡的整个变形过程中,壁面效应越强或流体的剪切稀化程度越弱,气泡的横向变形则越小,气泡整体变形程度也越小.
2)壁面效应越强或流体剪切稀化程度越弱,气泡尾涡越不易形成,气泡所受的压差阻力越大,所以气泡终端速度越小.
3)壁面效应和流体的剪切稀化程度会影响壁面与气泡间液相剪切速率的大小,强的壁面效应和强的剪切稀化程度导致高剪切速率区出现在壁面附近,引起该区域液相的表观黏度大幅度的下降.
4)壁面效应对气泡上浮运动特性的影响,与流体的性质有关;对于剪切稀化程度强的流体,气泡的上浮运动似乎更易受到壁面效应的影响.
5)在所分析的物理参数中,气泡终端速度对壁面效应的依赖性最大;如果为了消除壁面效应对气泡运动特性的影响,应以终端速度作为建模的判断依据.
[1] UNO S, KINTNER R C.Effect of wall proximity on the rate of rise of single air bubbles in a quiescent liquid[J].AICHE Journal, 1956, 2(3): 420-425.
[2] KRISHNA R, URSEANU M I, VAN BATEN J M,et al.Wall effects on the rise of single gas bubbles in liquids[J].International Communications in Heat and Mass Transfer, 1999, 26(6): 781-790.
[3] TSUGE H, HAMAMOTO S, HIBINO S.Wall effect on the behavior of single bubbles rising in highly viscous-liquids[J].Journal of Chemical Engineering of Japan, 1984, 17(6): 619-623.
[4] MUKUNDAKRISHNAN K, QUAN S P, ECKMANN D M, et al.Numerical study of wall effects on buoyant gas-bubble rise in a liquid-filled finite cylinder[J].Physical Review E, 2007, 76(3): 036308.DOI: 10.1103/PhysRevE.76.036308.
[5] ALHENDAL Y, TURAN A, KALENDAR A.Wall effects on the thermocapillary migration of single gas bubbles in stagnant liquids[J].Heat and Mass Transfer, 2017, 53(4): 1315-1326.
[6] LEE J, PARK H.Wake structures behind an oscillating bubble rising close to a vertical wall[J].International Journal of Multiphase Flow, 2017, 91: 225-242.
[7] 易妍妍, 王智慧, 杨超, 等.静止非牛顿流体中气泡生成过程的传质[J].化工学报, 2015, 66(11): 4335-4341.(YI Yanyan, WANG Zhihui, YANG Chao, et al.Mass transfer during single bubble growing in static non-Newtonian fluid[J].CIESC Journal, 2015, 66(11): 4335-4341.(in Chinese))
[8] OHTA M, KOBAYASHI N, SHIGEKANE Y, et al.The dynamic motion of single bubbles with unique shapes rising freely in hydrophobically modified alkali-soluble emulsion polymer solutions[J].Journal of Rheology, 2015, 59(2): 303-316.
[9] MIYAHARA T, YAMANAKA S.Mechanics of motion and deformation of a single bubble rising through quiescent highly viscous Newtonian and non-Newtonian media[J].Journal of Chemical Engineering of Japan, 1993, 26(3): 297-302.
[10] PREMLATA A R, TRIPATHI M K, KARRI B, et al.Dynamics of an air bubble rising in a non-Newtonian liquid in the axisymmetric regime[J].Journal of Non-Newtonian Fluid Mechanics, 2017, 239: 53-61.
[11] PORYLES R, VIDAL V.Rising bubble instabilities and fragmentation in a confined polymer solution[J].Journal of Non-Newtonian Fluid Mechanics, 2017, 241: 26-33.
[12] XU X F, ZHANG J, LIU F X, et al.Rising behavior of single bubble in infinite stagnant non-Newtonian liquids[J].International Journal of Multiphase Flow, 2017, 95: 84-90.
[13] KISHORE N, NALAJALA V S, CHHABRA R P.Effects of contamination and shear-thinning fluid viscosity on drag behavior of spherical bubbles[J].Industrial & Engineering Chemical Research, 2013, 52(17): 6049-6056.
[14] 黄萌.电场作用下油水乳液中水滴的运动和形变特性研究[D].博士学位论文.西安: 西安交通大学, 2015.(HUANG Meng.Investigation on motion and deformation of water drops in oil-water emulsion under electrostatic field[D].PhD Thesis.Xi’an: Xi’an Jiaotong University, 2015.(in Chinese))
[15] BRACKBILL J U, KOTHE D B, ZEMACH C.A continuum method for modeling surface tension[J].Journal of Computational Physics, 1992, 100(2): 335-354.
[16] DEKEE D, CARREAU P J.Friction factors and bubble dynamics in polymer solutions[J].The Canadian Journal of Chemical Engineering, 1993, 71(2): 183-188.
[17] DIMAKOPOULOS Y, PAVLIDIS M, TSAMOPOULOS J.Steady bubble rise in Herschel-Bulkley fluids and comparison of predictions via the augmented Lagrangian method with those via the Papanastasiou model[J].Journal of Non-Newtonian Fluid Mechanics, 2013, 200: 34-51.
[18] BHAGA D, WEBER M E.Bubbles in viscous liquids: shapes, wakes and velocities[J].Journal of Fluid Mechanics, 1981, 105: 61-85.
[19] KUMAR P, VANKA S P.Effects of confinement on bubble dynamics in a square duct[J].International Journal of Multiphase Flow, 2015, 77: 32-47.
[20] 徐丞君, 徐胜利, 刘庆源.修正压力梯度粒子近似SPH方法计算大密度比界面流动[J].应用数学和力学, 2019, 40(1): 20-35.(XU Chengjun, XU Shengli, LIU Qingyuan.Modified particle approximation to pressure gradients in the SPH algorithm for interfacial flows with high density rations[J].Applied Mathematics and Mechanics, 2019, 40(1): 20-35.(in Chinese))