气泡群动力特性模拟分析

张淑君,吴锤结

(1.河海大学浅水湖泊综合治理与资源开发教育部重点实验室,江苏南京 210098;2.大连理工大学机械工程系,辽宁大连 116014)

气泡群运动的研究在环境工程、水利工程和航运工程等方面具有重要的理论和现实意义[1-3].如污水处理厂的曝气工艺,河道水环境治理中的人工受氧过程以及水源地水质改善等都与气泡群的运动密切相关.随着气泡数量的增加,气泡之间的影响和作用变得更加复杂,特别是当粒子浓度很大时对粒子的空间分布也有不同程度的影响.早期工业气泡流特性研究只集中于体特性,如平均压降、平均流量或平均剪切应力等,且仅限于一维[4],而更细致的气泡研究应着重于气泡本身和其引起的周围流场的三维运动特性,这要借助于各种界面模拟手段.Semereka等[5-6]采用边界积分方法模拟了大量气泡在势流中的运动,结果表明,气泡上升速度在经过一段时间后基本相同且在上升过程中气泡水平聚集,这一结论与试验结果不尽相同的原因是模拟过程中没有考虑气泡尾部漩涡的作用.Spelet等[7-8]模拟分析了湍流对气泡群平均上升速度的影响,发现气泡群在湍流中的上升速度只有其在静水中上升速度值的一半左右.由于文中把气泡假设成点颗粒,且没有考虑气泡间的相互作用,因此该结论只适用于体积分数很小的气泡运动.在考虑了更多因素如惯性、黏性等,Hu等[9-10]采用非结构网格模拟了100个固体颗粒的运动,但由于在每一时间步均须重构流场网格,与目前广泛采用的单流体模型相比有很大的不足.目前以Tryggvason[11-13]为代表的锋面跟踪法在模拟气泡群方面取得了一定进展,但多数假设气泡在上升过程中不发生碰撞、融合等现象.本文采用VOF(volume-of-fluid)中的PLIC(piecewise linear interface calculation)界面重构方法研究多个三维气泡的动力学特性及其对周围流场的扰动.

1 控制方程及其离散

1.1 控制方程

考虑表面张力的动量方程为

式中:σ——表面张力 ;p——压强;μ——动力黏度;κ——界面的曲率;δs——与界面有关的 Dirac分布;n——界面单位法向矢量;D——应力张量,满足

不可压缩流体连续方程

采用VOF方法模拟气泡界面的体积函数C满足

对于两相流,方程(1)中的 μ和ρ均由体积函数C决定

式中 ρ1,ρ2,μ1,μ2分别为2种流体的密度和动力黏度.

1.2 界面跟踪

本文采用Youngs的PLIC方法在单个网格内用直线段近似界面方法进行界面重构.

首先确定界面法向量.界面法向量m(不一定是单位矢量)定义在网格中心,步长为h,则

以二维为例,则网格中心法向量m可由4个角点处的m值确定

界面的位置由界面法向量和同样定义在网格中心的体积函数C确定,并采用Lagrangin方法跟踪随流动传播的界面.

物理量定义采用MAC交错网格,并采用经典的中心差分格式.另外,由于在界面处密度、黏度和压强突变,计算过程中采用多重网格技术,每个节点处的密度和黏度采用简单的网格体积平均进行计算,压强通过Possion方程进行求解.

为了验证模型,文献[14]模拟了气泡在水中的连续变形过程,并与Walters等[15]的试验结果进行了比较.

2 模拟结果分析

本文模拟了三维圆形气泡群在静止无限流场中重力的作用下上升和变形的过程.气泡初始速度为零,在浮力的作用下在静止的流场中开始上升,边界采用无滑移条件.模拟区域为长方体,计算网格为66×66×128,时间步长和空间步长分别为2.5×10-6和1.95×10-4,周围流体与气泡密度和黏度比分别为ρf/ρb=40,μf/μb=80,Eo=ρbgd2/σ描述气泡的尺寸,其中d为气泡直径.

2.1 气泡群的上升

气泡群在上升的过程中通常伴随着碰撞、变形和融合的过程,雷诺数及气泡数量对其亦有不同的影响.

图1给出了三维球形气泡群在上升过程中的变化情况,气泡数量为144个,气泡沿z上升方向排成4层,每层36个,且3个方向气泡间距相同,气泡在上升过程中保持球形形状.从图中可以看出,在无量纲时间τ=45时,气泡由最初的整齐排列开始变化,占据中心位置的气泡上升速度要高于周围气泡,这使得每层由平面变为向上升方向突出的曲面,靠近底部的气泡由于受到前部气泡的阻挡只得朝着与上升方向垂直的方向运动,整个群体呈现扇形体的形状.随着时间的推移,外围的气泡开始发生碰撞并且开始融合,从τ=75时刻对应的图中可以观察到第1个融合后体积变大的气泡随着群体继续上升,随后大量气泡的融合过程开始,其中包括左右相邻和前后排列的气泡之间的融合.在τ=111时刻气泡的第2次融合过程亦已发生,此时气泡的体积相差显著,气泡之间的上升速度差别亦增大,因此气泡群体在上升方向距离加大,并在随机的位置上继续上升、碰撞、融合.这一现象与气泡数量减至Nb=64个时变化过程类似,但由于气泡数量少的缘故,气泡融合时间开始得稍晚一些.

随着雷诺数的增加,例如Re=74,Eo=6.9,气泡变形显著,为椭圆形,且其上升速度与小雷诺数时相比高出很多,变化过程也大不相同.气泡在初期也表现出向两侧扩展的扇形体,但程度不太明显.值得注意的是,与表面相邻的第2层气泡几乎整体与其发生碰撞、融合形成大气泡,留下3排体积不等的气泡继续上升,此时气泡间速度差别开始加大,但接下来的第2次融合过程与第1次类似,仍为整体融合,气泡数量再次大量减少.

泡群运动引起的速度分布见图2.从图2可以清楚地看到与少量气泡相比,多个气泡的同时出现使得气泡间的相互作用、相互影响范围扩大,因此对周围流体的扰动也大大增强.但当Re=74时,流体黏性降低,气泡上升速度加快,在其影响下,外围流体速度相对Re=0.8要大得多,气泡对流场的扰动也进一步增强.

图1 气泡群的上升过程(Nb=144,Eo=0.36,Re=0.8)Fig.1 Rising processes of bubble clusters(Nb=144,Eo=0.36,Re=0.8)

图2 气泡群运动引起的流场速度分布(Nb=144)Fig.2 Distribution of velocity fields induced by motion of bubble clusters(Nb=144)

图3 气泡群平均上升速度(Nb=144)Fig.3 Average rising velocities of bubble clusters(Nb=144)

与上升过程相对应可以分析气泡群整体上升速度的变化.图3给出了气泡群的平均上升速度,其中u*为速度.当Re=0.8,Eo=0.36时,初期速度由零迅速增大至0.8左右,开始振荡上升,此时由于前排气泡的阻挡,气泡群整体形状向外扩展成扇状体,上升速度趋缓.随着气泡之间距离的减小,少量气泡间开始发生碰撞与融合.从融合开始到第1个融合后大气泡的生成时段,气泡群速度经历了第1个融合后的速度降低与反弹[10],例如 τ=75时速度为1.15.随着大量体积变大气泡的出现及碰撞、融合的频繁发生,一部分气泡体积变大,上升速度加快,进而带动整体上升速度,表现为斜率变大、振幅增强和多个峰值的出现.图3(b)考察了气泡形状为椭圆形时平均速度的变化过程,显然其上升速度要远远高于小Re时的对应值.总体来说,气泡群速度在上升到9.0左右以后变化不大,这是由于气泡群体扩散现象不明显的原因.随后在第1次2排气泡的融合过程中,速度出现了明显的先降后升现象,并在第2次的大量融合时段出现第2个整体速度的降低与反弹现象.

2.2 气泡群引起的速度脉动

图4给出了在融合发生之前气泡表面的速度脉动,所跟踪气泡分别位于各层的中心,采样无量纲时间间隔为0.15,采样点数为400个.可以看出在气泡体积增大和上升速度同时增大的 τ=0～30时段,气泡表面x,y,z3个方向的速度u′,v′,w′脉动振幅都很大,这一现象在气泡数量较多Nb=144时更加明显.这是由于气泡数量的增加使得相互之间作用和影响增强的结果.但在数量较少时不同层面之间的气泡表面脉动程度相差较大,例如Nb=64时第2层速度脉动与位于表面第1层气泡间的差别较明显,并在最大速度到来之前出现2个脉动峰值.

图4 气泡表面的速度脉动(Eo=0.36,Re=0.8)Fig.4 Fluctuating velocities on bubble surface(Eo=0.36,Re=0.8)

3 结 语

本文采用VOF(volume-of-fFluid)中的PLIC(piecewise linear interface calculation)界面重构方法模拟研究三维气泡群的动力学特性及其对周围流场的扰动.结果表明,在气泡群上升过程中,伴随着随机的碰撞与融合等现象,气泡分布经历了由规则到扇形体再到随机的过程;在大雷诺数时融合过程表现为整体行为.与此相对应,气泡群平均上升速度在小雷诺数时随着时间的推移逐渐加快,并在气泡发生融合过程中出现震荡.气泡数量的增多对气泡本身和外围流体的速度扰动都有所增强.

[1]郑源,索丽生,屈波,等.有压输水管道系统含气水锤研究[J].河海大学学报:自然科学版,2005,33(3):277-281.(ZHENG Yuan,SUO Li-sheng,QU Bo,et al.Study onwater hammer with gas in pressuriedwater delivery system[J].Journal of Hohai University:Natural Science,2005,33(3):277-281.(in Chinese))

[2]FERRANTEA,ELGHOBASHI S E.On the accuracy of the two-fluid formulation in DNS of bubble-laden turbulent boundary layers[J].Physics of Fluids,2007,19:1-8.

[3]张淑君,吴锤结,王惠民.单个三维气泡运动的数值模拟[J].河海大学学报:自然科学版,2005,33(5):534-537.(ZHANG Shu-jun,WU Chui-jie,WANG Hui-min.Numerical simulation of the rise and deformation of a three dimenssionalgas bubble[J].Journal of Hohai University:Natural Science,2005,33(5):534-537.(in Chinese))

[4]HETSRONI G.Handbook of multiphase systems[M].Washington:Hemisphere-McGraw Hill,1982:550-556.

[5]SMEREKA S.On the motion of bubbles in a periodic box[J].J Fluid Mech,1993,254:79-112.

[6]SMEREKA S.A vlasov equations for soundwave propagation in bubbly fluids[J].J Fluid Mech,2002,454:287-325.

[7]SANGANI A S,DIDWANIA A K.Dynamic simulationsof flow of bubbly liquids at large Reynolds numbers[J].JFluidMech,1993,250:307-337.

[8]PELET P D,BIESHEUVEI A.On the motion of gas bubbles in homogeneous isotropic turbulence[J].J FluidMech,1997,336:221-244.

[9]HU H H.Direct numerical simulations of flows of solid-liquid mixtures[J].Phys Fluids,1996,22:335-352.

[10]JOHNSON A A,TEZDUYAR T E.3D simulation of fluid-particale interactions with the number of particles reaching 100[J].Comput Meth Appl Mech Eng,1997,145:301-321.

[11]RADL S,TRYGGVASON G,KHINAST J.Flow and mass transfer of fully resolved bubbles in non-Newtonian fluids[J].AIChE Journal,2007,53:1861-1878.

[12]LU J,TRYGGVASON G.Numerical study of turbulent bubbly downflows in a vertical channel[J].Physics of Fluids,2006,18(10):103-302.

[13]LU J,BISWAS S,TRYGGVASON G.Study of laminar bubbly flows in a vertical channel[J].International Journal of Multiphase Flow,2006,32:643-660.

[14]张淑君,吴锤结.气泡之间相互作用的数值模拟[J].水动力研究与进展:A辑,2008,23(6):681-686.(ZHANG Shu-jun,WU Chui-jie.Numerical simulation of the intera between two three-dimensional deformable bubbles[J].Chinese Journal of Hydrodynamics,2008,23(6):681-686.(in Chinese))

[15]WALTERS J K,DAVIDSION J F.The initialmotion of a gas bubble formed in an inviscid liquid,part 1:the two dimensional bubble[J].J Fluid Mech,1962,12:408-416.