槽道湍流内气泡瞬态受力数值研究

庞明军，孙涛，代军荣（常州大学机械工程学院，江苏 常州 213164）



研究开发

槽道湍流内气泡瞬态受力数值研究

庞明军，孙涛，代军荣
（常州大学机械工程学院，江苏 常州 213164）

摘要：为了深入理解气泡在湍流场内的运动机制，使用欧拉-拉格朗日单向耦合数值方法，详细分析了气泡在低雷诺数槽道湍流场内的瞬态受力情况；其中液相湍流速度场采用直接数值模拟方法求解，气泡的瞬态受力由牛顿运动方程计算；计算时，考虑了相间阻力、剪切升力、压力梯度力、虚拟质量力、重力对气泡运动的影响。目前的计算结果表明：气泡所受的瞬态作用力分量同时取决于重力作用方向、液相流向以及气泡所处的法向位置；不同方向上、不同位置处影响气泡运动的主要作用力分量是不同的；相比较而言，与重力方向垂直的剪切升力分量在近壁面区域为影响气泡运动的主要作用力，压力梯度力的法向分量在近壁面区域之外为影响气泡运动的主要作用力，相间阻力分量在整个槽道区域内均为影响气泡运动的主要作用力，除了竖直槽道近壁面处之外、虚拟质量力也均为影响气泡运动的主要作用力。

关键词：欧拉-拉格朗日模型；单向耦合；直接数值模拟；气泡；受力分析

鉴于良好的物理特性，泡状流被广泛应用于金属冶炼、矿物浮选、泡状反应堆、微气泡减阻和流体输送等工程领域。对于泡状流而言，气泡在液相中的分布会直接影响液相的流动特性、压降分布以及传质传热等特性，为了理解气泡的分布模式和分布机理，国内外学者开展了大量的研究[1-4]。然而由于问题的复杂性，目前关于气泡的分布模式仍未取得一个普遍认可的机理解释。特别是从气泡受力角度出发，研究气泡运动和分布的研究者为数不多。陈凤等[5]研究了电场力作用下沸腾气泡的电场力情况，闫超星等[6]研究了摆动情况下气泡惯性力和浮力的变化情况。据本文作者所知，仅 GIUSTI等[7]给出了湍流场内气泡受力时均值的分布情况。事实上，气泡的瞬态运动与其瞬态受力有关，研究气泡在湍流场的瞬态受力情况，有助于了解处于湍流场不同位置的气泡，支配其运动的主要相间力有哪些，也有助于理解气泡是如何从湍流外层运动到湍流内层。鉴于实验测量的困难性，本文使用数值模拟方法对槽道湍流内气泡的瞬态受力进行了深入的分析和研究。

1 计算条件和数值方法

1.1 计算条件

因为相间作用力（如相间阻力和剪切升力）与气泡的滑移速度密切相关，而气泡的滑移速度又与气泡所受浮力有关，因此，为了理解重力在槽道内作用方向对气泡瞬态受力的影响，本文设计了两大类计算工况：一类是让气泡浮力方向与液相主流方向一致；另一类是让气泡浮力方向与液相主流方向垂直。同时为了研究气泡位置对气泡瞬态受力的影响，每一大类计算工况又分别包含了3种计算工况，即让气泡分别位于黏性底层（即气泡位置y+=3.75）、过渡层（即气泡位置 y+=10）和湍流外层（即气泡位置 y+=145），以便观察不同重力方向和不同气泡位置气泡瞬态受力情况。其中y+为量纲为1的壁面位置坐标，气液相的具体计算参数见表 1。目前使用的网格已在文献[8-9]进行了详细的精度验证，本文不再赘述，表明网格精度能够满足计算要求。

计算区域为一个二维槽道，几何形状和详细尺寸见图1，图中x、y、z方向分别对应液相的流向、法向和展向。如图所示，对于竖直槽道，重力方向与流向平行；当计算水平槽道时，将重力方向逆时针旋转90°与流向垂直。计算时，对于液相，在流向和展向设置周期边界条件，在壁面法向设置无滑移边界条件；对于气泡，同样在流向和展向设置周期边界条件，为了更好地观察气泡在不同位置时瞬态受力情况，气泡在法向的位置始终保持不变。

表1 气液相计算参数

图1 计算区域和尺寸

1.2 控制方程

目前采用欧拉-拉格朗日单向耦合模型模拟气液两相流动，计算时不考虑气泡对液相流动的影响。为了准确求解液相的瞬态速度场，液相流场采用直接数值模拟方法求解，气泡的瞬态受力由牛顿运动方程进行计算。计算时，所有的控制方程进行了量纲归一化处理，与长度相关的量用半槽道高度（h）量纲为1化，与速度相关的量用摩擦速度（uτ）量纲为1化，与时间相关的量用h/uτ量纲为1化，具体过程详见文献[8]。

量纲归一化液相的连续性和动量方程见式(1)、式(2)。

量纲归一化的气泡运动方程见式(3)。

式(3)中，等号左边项为气泡惯性力项，等号右边项分别为压力梯度力、相间阻力、虚拟质量力、剪切升力和浮力项；CD、CV和CLF分别为相间阻力系数、虚拟质量力系数和升力系数。对于球形气泡，虚拟质量力系数为 0.5[10]，阻力系数和升力系数分别由式(4)[11]和式(5)[12]来计算。

在式(5)中

式中Srb、Reb分别为量纲归一化的剪切率和气泡雷诺数，ζ=SrbReb。

1.3 数值方法

目前的计算由自编程序完成，液相控制方程采用有限差分法进行离散，空间项采用二阶精度的差分格式；在流向和展向进行均匀网格划分，在法向进行非均匀网格划分，越靠壁面网格越密；为了避免压力场的不稳定性，采用交错网格系统，速度存储于单元格的边界上，而压力存储于网格中心；时间的推进，其压力项采用隐式格式，速度的时间推进采用二阶精度的Adams-Bashforth格式；速度和压力的耦合采用MAC方法；压力泊松方程由多重网格进行求解。采用二阶精度的 Crank-Nicholson格式对微气泡的加速度方程进行积分，得到微气泡的速度；利用同样的方法对速度方程积分得到气泡流向和展向上的位移。有关计算方法和计算模型的详细介绍和准确性验证详见文献[8]，该作者已采用该模型成功模拟了微气泡对液相湍流的调制情况。

2 结果分析

下面给出不同工况气泡瞬态受力随计算时间的变化情况，图 2、图 3中纵坐标为气泡瞬态受力情况，横坐标为量纲归一化的计算时间[t*=t/(h/uτ)]。图中负的瞬态力表示该力指向坐标轴负方向。

2.1 竖直槽道气泡瞬态受力分析

图2给出竖直槽道内气泡瞬态受力情况，即气泡所受浮力方向与液相流向一致情况下，气泡在流向和法向上的受力情况，有关相间力的具体计算公式详见式(3)。从图2可以看出，因湍流场的不稳定性，处于不同位置的气泡所受相间作用力随着计算时间的推移在作不稳定的波动，特别是处于过渡层和外层上的气泡，其法向分力波动幅度较大；而且气泡所受相间力的大小与其所处位置密切相关。在流向上，随着远离壁面，气泡所受的相间阻力（FD）和虚拟质量力（FVM）在增大；在法向上，随着远离壁面，气泡所受的剪切升力（FL）在减小，相间阻力（FD）是先增大后减小，而虚拟质量力（FVM）和压力梯度力（Fp）是在逐渐增大。

从图2中可以清楚得看到，气泡所处位置不同，不同方向上控制气泡运动的相间力也不同。当气泡位于黏性底层时（y+=3.75），在流向上控制气泡运动的力主要为浮力（FG）和相间阻力（FD），在法向上控制气泡运动的力主要为剪切升力（FL）和相间阻力（FD）。当气泡位于过渡层时（y+=15），在流向上控制气泡运动的力主要为浮力（FG）和相间阻力（FD），虚拟质量力（FVM）也不容忽视；在法向上控制气泡运动的力主要为压力梯度力（Fp）和虚拟质量力（FVM），但相间阻力（FD）也较为重要。气泡位于湍流外层与气泡处于过渡层相似，流向上控制气泡运动的力也主要为浮力（FG）和相间阻力（FD），同样虚拟质量力（FVM）也较为重要；但在法向上，不同于气泡位于过渡层时的情况，控制气泡运动的力主要为虚拟质量力（FVM）和压力梯度力（Fp）。

图2 竖直槽道内气泡瞬态受力情况

2.2 水平槽道气泡瞬态受力分析

图3给出了水平槽道内气泡瞬态受力情况，即液相主流方向与浮力（重力）方向相垂直情况下，气泡的瞬态受力。类似于竖直槽道内气泡的受力情况，水平槽道内气泡的瞬态受力也随着计算时间和气泡位置在变化。当气泡位于黏性底层时，在流向上控制气泡运动的力主要有阻力（FD）、剪切升力（FL）和虚拟质量力（FVM）；在法向上控制气泡运动的力主要为浮力（FG）和相间阻力（FD）。当气泡位于过渡层和湍流外层时，流向上控制气泡运动的力主要为相间阻力（FD）和虚拟质量力（FVM）；而法向上控制气泡运动的力主要为浮力（FG）、相间阻力（FD）、虚拟质量力（FVM）和压力梯度力（FP），而且阻力的方向基本与浮力的方向相反。

2.3 槽道内气泡受力分析总结

为了清楚展示气泡在槽道内的受力情况，根据图2和图3给出了微气泡受力分析示意图，如图4所示。可以看出，控制气泡运动主要的力与重力（浮力）方向密切相关。当气泡位于黏性底层时，在重力方向上，控制气泡运动主要的力为浮力（FG）和相间阻力（FD）；在与重力相垂直的方向上，控制气泡运动的力主要为剪切升力（FL）和相间阻力（FD），对于水平槽道，虚拟质量力（FVM）不容忽视。当气泡位于过渡层和湍流外层时，在重力方向上，控制气泡运动的力主要有浮力（FG）、相间阻力（FD）和虚拟质量力（FVM），对于水平槽道，压力梯度力（Fp）也很重要；在与重力相垂直的方向上，重力的影响比较复杂，对于竖直槽道而言，控制气泡运动的力主要为压力梯度力（Fp）和虚拟质量力（FVM），对于水平槽道而言，主要的力为阻力（FD）和虚拟质量力（FVM）。

图3 水平槽道内气泡瞬态受力情况

作为一个事实，重力）（浮力）的存在会使气泡在重力方向上比其周围的液体运动得快，即存在一个负的滑移速度，导致在重力方向上，气泡所受相间阻力的方向始终与浮力方向相反，且二者的数值也很大；另外，滑移速度的存在，再加上局部液相的速度梯度，必然会在与重力相垂直的方向上产生一个剪切升力，因液相的局部速度梯度在壁面附近最大，所以气泡在黏性底层时其所受的剪切升力较大，对其运动起主导作用。随着远离壁面，液相的速度梯度急剧下降，剪切升力变得不再重要。虚拟质量力是因气泡运动拖动局部液体一起运动而产生的一种力，所以对气泡而言是一种附加阻力，因此方向也始终与浮力方向相反。

图4 气泡主要受力示意图

对于充分发展的槽道湍流，当流动达到统计稳态时，液相在流向上的速度梯度基本为零，对应的短距离内的压力梯度也为零，因此在流向上气泡的压力梯度力可忽略不计；而液相在法向上的速度梯度较大，所以在法向上，气泡所受压力梯度力较大，且随着流场的变动在不断调整和变化。

从目前的分析可以看出，不论是水平槽道还是竖直槽道，气泡从槽道中心区域向壁面附近区域的运动与压力梯度力和虚拟质量力的法向分量有直接的关系。

3 结 论

利用欧拉-拉格朗日单向耦合数值方法详细调查了竖直和水平槽道湍流场内、不同法向位置处气泡的瞬态受力情况，可以总结以下几点结论。

（1）气泡瞬态受力随湍流场不断变化，其大小和方向同时取决于重力方向、液相流向以及气泡在槽道内壁面法向上的位置；不同方向上、不同位置处，影响气泡运动的主要作用力分量是不同的。

（2）与重力相垂直的剪切升力分量在近壁面区域处为控制气泡运动的主要作用力，在近壁面之外区域、压力梯度力法向分量为控制气泡运动的主要作用力，不同方向上的相间阻力和虚拟质量力分量均为控制气泡运动的主要作用力。

（3）相同工况下，法向上的虚拟质量力分量与法向上的压力梯度力分量有关；重力方向上的相间阻力分量取决于气泡所受浮力，非重力方向上的相间阻力在远离壁面区域处较小且随湍流场波动。

（4）气泡从槽道中心区域向壁面附近区域的运动与压力梯度力和虚拟质量力的法向分量有直接的关系。

参考文献

[1] SHEN X Z，MISHIMA K，NAKAMURA H.Two-phase phase distribution in a vertical large diameter pipe[J].Int.J.Heat Mass Transfer，2005，48（1）：211-225.

[2] OKAWA T，KATAOKA I，MORI M.Numerical simulation of lateral phase distribution in turbulent upward bubbly two-phase flow[J].Nucl.Eng.Des.，2002，213（2/3）：183-197.

[3] 刘国强，孙立成，田道贵，等.竖直圆管内空泡份额径向分布特性形成机制[J].核动力工程，2014，35（1）：42-45.

[4] 庞明军，徐一丹，魏进家.管道泡状流相分布模式和分布机理研究进展[J].化工进展，2014，33（11）：2829-2841.

[5] 陈凤，宋耀祖，陈民.电场作用下的气泡受力分析[J].工程热物理学报，2005，26（6）：146-148.

[6] 闫超星，阎昌琪，孙立成.周期性力场作用下气泡径向附加力分析[J].原子能科学技术，2013，47（12）：2220-2224.

[7] GIUSTI A，LUCCI F，SOLDATI A.Influence of the lift force in direct numerical simulation of upward/downward turbulent channel flow laden with surfactant contaminated microbubbles[J].Chem.Eng.Sci.，2005，60（22）：6176-6187.

[8] PANG M J，WEI J J，YU B.Numerical study of bubbly upflows in a vertical channel using the Euler–Lagrange two-way model[J].Chem.Eng.Sci.，2010，65（23）：6215-6228.

[9] YU B，KAWAGUCHI Y.Direct numerical simulation of visco-elastic drag-reducing flow: a faithful finite difference method [J].Journal of Non-Newtonian Fluid Mechanics， 2004，116（2/3）：431-466.

[10] TABIB M V， ROY S A，JOSHI J B.CFD simulation of bubble column——an analysis of interphase forces and turbulence models [J].Chemical Engineering Journal，2008，139（8）：589-614.

[11] LAIN S， BRODER D， SOMMERFELD M， et al.Modelling hydrodynamics and turbulence in a bubble column using the Euler-Lagrange procedure [J].Int.J.Multiphase Flow，2002，28（8）：1381-1407.

[12] LEGENDER D， MAGNAUDET J.The lift force on a spherical bubble in a viscous linear shear flow [J].J.Fluid Mech.， 1998，369：81-126.

第一作者及联系人：庞明军（1976—），男，博士，副教授，研究方向为气液两相流与湍流减阻。E-mail pangmj@cczu.edu.cn。

中图分类号：TK 124

文献标志码：A

文章编号：1000-6613（2016）07-1980-06

DOI：10.16085/j.issn.1000-6613.2016.07.007

收稿日期：2015-11-18；修改稿日期：2015-12-14。

基金项目：国家自然科学基金（51376026）及江苏省大学生创新创业训练计划重点项目（201510292016Z）。

Numerical investigation on instantaneous forces exerting on bubbles in channel turbulent flows

PANG Mingjun，SUN Tao，DAI Junrong
（School of Mechanical Engineering，Changzhou University，Changzhou 213164，Jiangsu，China）

Abstract：In order to deeply understand the motion mechanism of bubbles in turbulent flows，instantaneous forces exerting on bubbles in a channel turbulence flow with a low Reynolds number were analyzed in details with the Euler-Lagrange one-way method.The liquid-phase turbulence field was computed with direct numerical simulations（DNS） in Euler reference frame，and instantaneous forces on bubbles were solved with Newtonian motion equations including drag，shear lift，pressure gradient，virtual mass，gravity（i.e.，buoyant） and inertia forces in Lagrange reference frame.The present investigations showed that forces acting on bubbles depended simultaneously on the gravity direction，the flow direction of liquid and the bubble place in the wall-normal direction.The instantaneous force components controlling bubble motions changed with bubbles location and motion direction.Comparatively speaking，the component of the shear lift force perpendicular to gravity had an important influence on the bubble motion in the region near walls.The wall-normal component of the pressure gradient force was dominant for bubble motions in the region away from walls.The drag force components was very important for bubble motions in the whole channel region，and the virtual mass force played a dominant role in bubble motions in all regions except one near the walls for the vertical channel.

Key words：Euler-Lagrange model；one-way；direct numerical simulations；bubbles；force analysis