王 乐， 苏军伟， 郑西朋， 杨顺生

(1.西南交通大学土木工程学院，四川 成都610031;2.西安交通大学人居环境与建筑工程学院，陕西 西安710049)

化工领域的鼓泡塔作为一种常见的曝气反应器，由于其结构简单、传质效率高等优点，也被用于环境工程领域［1］，并被称为塔式曝气池(后文简称曝气池)［2］.对曝气池进行的研究，有利于优化曝气方式、提高曝气池气含率、改进曝气池结构.数值模拟方法由于良好的经济性，并且可以获得曝气池内流场形态、气相分布、湍动能、气含率等更多的细节信息，被广泛地应用于研究曝气池气泡羽流特性［3］、液相速度［4］、气泡直径［5］等动力学行为.

曝气池内，气泡由于受到表面张力、升力、曳力等的影响，在上浮的过程受气泡之间的不同速度差、湍流以及气泡尾涡等因素发生聚并;当湍流涡体的惯性力大于气泡表面张力引起的附加力时，气泡不断变形直至破碎成更小粒径气泡［6］.传统的CFD(computational fluid dynamics)方法将气泡视为单一粒径，未考虑气泡的聚并和破碎效应，忽视了气泡大小分布对相间作用力以及液相湍动的影响，无法准确计算气液体系的相间面积及预测气液体系的传质行为［7］.众多学者采用群体平衡模型考虑了气泡的聚并和破碎效应，并与欧拉双流体模型耦合对气液两相流动进行研究［8-10］.Wang 等［8］采用CFD-PBM(population balance model)对鼓泡塔进行模拟，结果显示该模型对于预测塔式曝气池的气泡尺寸分布、界面面积、气液传质速率方面准确有效.Gupta等［9］采用CFD-PBM方法对矩形鼓泡塔进行数值模拟，分析不同的聚并破碎求解方法以及升力、虚拟质量力对曝气池内流场及气含率的影响.徐琰等［10］针对咪唑类离子液体介质，基于CFDPBM方法，研究不同表观气速下鼓泡塔内气液两相速度场分布、气含率以及气泡尺寸等流体动力学性质，取得了满意的效果.Li等［11］对柱形鼓泡塔8种不同入口分布及形状进行了模拟，显示增加通气管能够有效改善塔内气含率.

然而以上的研究多针对气液相间曳力、液体性质、气泡聚并破碎过程的求解算法等方面，对气泡的聚并和破碎效应下气体入口位置不同导致曝气池内气相分布、液相流场、气泡数密度等方面仍鲜有研究.本文采用FLUENT流体力学计算软件，基于欧拉双流体模型耦合群体平衡模型，在对计算域网格及气相分布与文献［12］的曝气试验对比验证的基础上，对不同入口位置下的实验室尺度曝气池内气液两相流进行研究，分析曝气池内气含率、气泡数密度、流场等流体动力学行为.研究结果为曝气池的优化设计及入口位置选择提供更多依据.

1 数学模型

采用欧拉双流体模型用于模拟气液两相流动，并耦合群体平衡模型模拟气泡的聚并与破碎效应，模拟过程中忽略气液两相间的质量和热量传递，且认为气相和液相为不可压缩流体.

1.1 欧拉双流体模型

质量守恒方程

式中:α为体积分数;ρ为密度，kg/m3;u为速度;τ为剪切应力，Pa;p为压强，Pa;g为重力加速度，9.8 m/s2;下标q为相区分，为g时表示气相，为 l时表示液相;F为两相的相间力，N;本文考虑了升力FL及曳力FD.如方程(3)～(5)所示.

动量守恒方程

式中:d为气泡直径，m;CD为曳力系数，采用适用于不同形状气泡的 GRACE模型［13］，具体如方程(6)～(14)所示.

式中:Re为雷诺数;

其中:M为Morton数，

J为分段函数，

其中，E为Eotvos数，

曝气池内气液两相湍流求解采用混合RNG(re-normalization group)k-ε模型，该模型被广泛用于曝气池内湍流的模拟［14-15］，具体表述为

式中:k为湍动能，m2/s2;ε为湍动能耗散率，m2/s3;C1、C2为常数;Gk，q为湍动能产生项;Πk，q和 Πε，q为气相对液相影响源项.

1.2 群体平衡模型

Hagesaether等［16-17］将气液体系宏观尺度群体平衡模型表述如下:

式中:n(v，t)为数量密度函数;Ba、Da、Bb、Db分别为气泡聚并生成项、气泡聚并消亡项、气泡破碎生成项和气泡破碎消亡项;β(v，v')为概率分布密度函数;g(v')为破碎频率;a(v，v')为聚并频率.

采用离散方法求解方程(17)，模拟破碎过程采用 Luo等［18］提出的机理模型;聚并过程采用Luo［19］提出的机理模型.离散过程中不同粒径的气泡分为10组，具体的气泡分组见表1所示［12］.

表1 气泡直径分组Tab.1 Bubble size group discretization

2 模型设置及计算方法

2.1 模型及网格

Ali等［4］研究发现塔式曝气池二维和三维模型在速度预测等方面较为一致，二维模型能够反映曝气池中的气液两相作用.因此建立了与罗玮［20］、肖浩飞［21］以及 Díaz 等［12］类似或一致的矩形曝气池模型作为二维简化物理模型，如图1(a)所示.模型的高度为 0.450 m，模型宽度为 0.200 m，设置了4个不同的入口位置算例，入口位置中心分别距离原点 0.100、0.110、0.125、0.146 m，分别称为算例1、算例 2、算例 3以及算例 4，入口的长度为0.018 m.

本文对5种不同单元数量的均一的结构化网格(5 763、7 353、10 050、14 400、22 500)进行模拟，模拟结果与实验对比后(见3.1节)，选择单元数量为10 050的网格划分，其中x方向网格数量为67，y方向网格数量为150，如图1(b)所示.

2.2 计算方法

采用有限体积法对第1章节中的方程离散，其中体积分数项采用QUICK格式进行离散，其他项的离散采用二阶迎风格式，速度和压力采用耦合的SIMPLE算法，入口采用速度进口边界条件，空气进口速度值为 0.026 67 m/s［22］，入口气泡粒径为5.95 mm，出口采用脱气边界条件，其他物理边界为固壁边界条件.在稳态下计算1 000步作为非稳态计算的初始场，并继续在非稳态下计算100 s，残差设置为 10－4，时间步长为 0.01 s，最大迭代步为 30 步.

图1 模型及网格划分(单位:m)Fig.1 Physical model configurations and mesh(unit:m)

3 验证

3.1 网格验证

对2.1节算例1的物理模型在5种不同的网格划分下，采用如上所述的数学模型和计算方法，对文献［12］中的曝气过程进行模拟，模拟过程考虑了最大迭代步及时间步长的影响，如表2所示.

表2 网格及计算参数验证Tab.2 Verification of mesh and calculated parameters

从表2可见网格数量增多达到22 725时，模拟结果反而与实验结果吻合较差，这与Díaz等［12］和Buwa等［23］的研究结果一致，其可能原因与湍动谱有关.时间步长及最大迭代步也会影响模拟结果，因此选择了网格数量为10 050的网格，时间步长为0.010 0 s，最大迭代步为30步作为计算参数进行模拟.此外，进一步验证了模拟得到的气相分布并与实验进行对比，见3.2节.

3.2 气相分布验证

采用上文所述的物理模型、网格划分、计算方法等对文献［12］中的曝气过程进行模拟，气相分布的模拟结果与实验对比如图2所示，发现模拟所得气含率分布呈现“之”字形，入口附近气含率最高，出口附近气含率较低，由底部至顶部气体分布范围逐渐变大与实验所得的气含率分布结果基本一致.

图2 瞬态气含率模拟与实验对比Fig.2 Comparison between experimental and computational results of instantaneous gas hold-up

图 3 为监测点(0.100，0.225)m 液相水平速度瞬时变化，选择与文献［12，22］相同的二维坐标监测点，显示速度随时间呈周期性震荡，这也意味着气泡羽流在曝气池中呈周期摆动现象，与实验观察相符.为进一步验证模拟结果的正确性，对模拟得到的总的气含率及气泡羽流的震荡周期与实验对比，同时为验证欧拉双流体模型耦合PBM的优越性，一并给出在相同的计算条件和数学模型下仅采用欧拉双流体模型单一粒径(5.95 mm)下的模拟结果，如表3所示.

图3 监测点瞬态水平方向液相速度Fig.3 Instantaneous horizontal liquid velocity at observation points

表3 气含率及气泡羽流摆动周期模拟结果与试验结果对比Tab.3 Comparison between experimental and calculated results of gas hold-up and bubble plume oscillation period

其中采用欧拉双流体模型耦合PBM模拟得到的总气含率与实验极为接近且优于欧拉双流体模型与实验2.60%的误差，此外采用欧拉双流体模型耦合PBM模拟得到的气泡羽流周期与实验误差为1.79%优于欧拉双流体模型的3.16%，进一步验证了本文方法的正确性及优越性.

4 结果与讨论

4.1 液相流动特性分析

图4为曝气池内液相速度矢量图，液相在气相驱动下向上流动经塔顶后沿两侧塔壁回流至塔底，在此过程中受入口位置不同的影响形成了不同形状的旋涡结构.当入口位置位于曝气池底部中心时，在曝气池两侧分别形成两个对称的漩涡，漩涡中心位于曝气池偏下部，这与Díaz等［12］模拟得到的结果一致.随着入口位置逐渐远离中心，如图4(a)～(d)所示，曝气池左侧上方、中部右侧以及底部区域形成了3个明显的漩涡，这主要由于气相贴近曝气池右侧上浮并驱动液相流动，在顶部区域受边界条件影响，部分液相从曝气池左侧回流形成左侧上方旋涡，剩余液相则从右侧回流形成中部右侧旋涡，在底部区域液相流动受两侧壁面及上部液相回流影响形成底部的旋涡.随着入口位置逐渐偏离曝气池底部中心，漩涡的强度逐渐增强.

图4 液相时均速度Fig.4 Time-aver a ged liquid velocity field

图5 (a)～(c)为监测点的水平方向液相速度，结合图3(算例1)的模拟结果，发现不同算例中均呈现水平液相速度周期性震荡，这意味着瞬时下曝气池内气泡羽流和流场结构也呈现着周期性变化.

图5 监测点瞬态水平方向液相速度Fig.5 Calculated instantaneous horizontal liquid velocity at observation point

此外，液相水平速度的振幅随入口位置逐渐偏离曝气池中心而逐渐降低，速度的震荡周期也逐渐加长.不同于算例1、2、3，在算例4中，一个完整的震荡周期内液相水平速度的极大值为负值，结合图4可以发现，此时底部区域的漩涡向上方膨胀并压缩上方漩涡导致监测点速度只有大小变化没有方向改变.

4.2 气相时均垂向速度分布及羽流摆动周期

图6为气相垂直方向速度分布，不同算例垂向速度分布均呈现在曝气池上部中间区域速度较高，意味着气泡受浮力等作用力加速上浮.随着入口位置逐渐偏离曝气池中心，曝气池底部区域气相垂向速度的峰值逐渐增大，尤其在算例3和4中，气相速度在曝气池底部靠近右侧壁面区域有明显的高值分布，其原因是由于气泡浮射流在边壁附近区域由于附壁效应(也称为Coanda效应)会向边壁倾斜甚至贴着边壁流动［24］，附壁效应的存在导致曝气池底部区域气泡上升路径变窄，相同曝气量下气速峰值不断增大.在不同算例曝气池靠近壁面以及底部区域，存在气体垂向速度分布极低的区域，这主要是由于这些区域气相分布极少造成的，同时入口位置的不同也会导致这些区域面积大小及分布存在差异，当入口位置距离底部中心较远时，如算例3、4，气相垂向速度的极低区分布在曝气池底部左侧，当入口位置处于底部中心或距离中心较近时，垂向速度极低区分布在曝气池壁面两侧附近.

图6 时间平均垂直方向气相速度分布Fig.6 Time-averaged distribution of the vertical gas velocity

气泡在液相中上升时，由于受到周边压力不平衡的影响，会产生与气泡上升方向垂直的剪切诱导力，同时升力的作用也导致气泡沿着剪切速度较小的方向移动，这些因素导致曝气池内旋涡结构改变以及气泡羽流不断摆动.图7为不同算例下的气泡羽流摆动周期.随着入口位置距离曝气池中心越来越远，羽流摆动周期逐渐增大，算例1的摆动周期最小，算例3的摆动周期最大，但在算例4中羽流摆动周期反而减小，这可能是由于入口位置距离壁面较近，羽流摆动受右侧壁面边界的影响所造成的.

图7 气泡羽流摆动周期Fig.7 Bubble plume oscillation period

4.3 气相分布及气含率

如图8所示，为100 s时间内的气相时间平均分布，在不同入口位置算例中，气相总体呈底部区域分布面积小，顶部区域分布面积广，当入口位置偏离曝气池底部中心时，底部气相分布面积逐渐减小气含率值逐渐增大，在图8(d)算例4中底部区域气相的分布面积最小.在不同算例中均存在着不同面积的气含率极低区域，其主要分布在曝气池底部及入口位置对侧的壁面区域.图8(a)算例1中，由于入口位置及模型的对称性以及受到液相旋涡的影响(见图4)，气相分布形状近似柱形呈x=0.1 m对称;在图8(c)算例3及图8(d)算例4中，受附壁效应的影响气相在底部区域贴近壁面分布，在曝气池中部受瞬时针旋涡的剪切使得气相分布逐渐摆脱右侧壁面且分布面积进一步扩大，在曝气池顶部受逆时针旋涡影响，气相分布在顶部偏右侧区域，最终整个气相分布呈类“之”字形.通过以上分析可以得出，针对矩形塔式曝气池的设计过程中，应避免入口位置距离壁面较近的设计方式，这会导致底部曝气效果较差且曝气区域的实际面积较小等缺陷;采用入口位于塔式曝气池中心的设计，能够保证曝气区域面积较大且不同高度区域曝气更为均匀.

图8 曝气池时间平均气含率分布Fig.8 Time-averaged gas hold-up distribution for various cases

气含率作为影响曝气池传质的重要因素之一其大小对于曝气设计具有重要意义.较高的气含率能够促进好氧微生物代谢，加速分解污水中有机物.如图9所示为不同入口位置下曝气池内的气含率，发现算例1、2、3的气含率逐渐增大但增长幅度较小，气含率的最大值出现在算例3，而在算例4中，气含率最小.由于影响气含率的因素众多，包括气体停留时间、气相分布面积、气体速度、气相掺混能力等，其中算例1至算例3气含率增大的可能原因与旋涡强度较强湍动能较大有关，导致气液掺混能力较强气含率最大.而算例4中，虽然旋涡强度及湍动能较大，但底部气速过快导致气体停留时间较少，此外气相分布面积也较小，这些综合因素导致了气含率的降低.

图9 气含率对比Fig.9 Comparison of gas hold-up

4.4 气泡数密度

在以往对曝气池的模拟过程中，着重于对液相速度、气相分布以及气含率等指标进行描述，而对于每立方米气泡数量的多少(即气泡数密度)的分析较少.在同样的气含率下，气泡的直径越小，则气泡数量越多气液传质比表面积越大，从而微生物生长代谢能够利用的氧气量也越多.图10为不同入口位置下曝气池一个周期内的气泡数密度的时均统计，发现入口位置对于气泡数密度的大小及分布没有十分显著的影响，在小于5.95 mm的不同粒径分组中，气泡数密度普遍高出大于5.95 mm的不同粒径分组中气泡数密度1～2个量级，气泡数密度最大值出现在直径为5.95 mm的分组.导致小于5.95 mm粒径组气泡数密度远高于大于5.95 mm粒径组气泡数密度的原因是由于注入的气泡粒径为5.95 mm，并且低速下气泡聚并尤其是大气泡的聚并较难发生，而气泡的破碎是由于湍流涡体的惯性力大于气泡表面张力引起的附加力导致的，气泡的破碎效应较强，Buwa等［23］在实验中也观察到了在低速下气泡不易聚并这一现象.在算例3和4中，对比气泡粒径大于5.95 mm的气泡分组发现气泡数密度均多于相同粒径组中算例1和算例2，其主要是由于在算例3和算例4曝气池底部气相分布有明显的高值区且分布范围更窄，有利于气泡的碰撞聚并.

图10 不同分组下气泡数密度Fig.10 Bubble number densities in different size groups

5 结论

本文采用欧拉双流体模型耦合PBM模型对4种不同入口位置下塔式曝气池内气液动力学行为进行数值模拟，探讨入口位置对气含率、气泡数密度、气相速度分布等的影响，得到以下结论:

(1)欧拉双流体模型耦合群体平衡模型能够考虑真实气泡破碎与聚并效应，相比于单一粒径的欧拉双流体模型能够更好的模拟二维条件下曝气池气液两相规律，模拟结果与实验相符.

(2)随着入口位置逐渐偏离曝气池底部中心，曝气池内由两个对称分布的旋涡发展为分布在曝气池左侧上部区域、中部右侧区域以及底部区域的3个不同大小旋涡，且旋涡强度逐渐增强;瞬时液相水平速度呈周期性震荡，振幅随入口位置远离曝气池中心而不断缩小.

(3)曝气池内气含率及其分布、气相垂向速度峰值、气泡羽流摆动周期均与入口位置有关，当入口位置偏离曝气池底部中心时，气相垂向速度峰值增大，气相的底部分布区域面积减小，而气含率及气泡羽流摆动周期均先增大后减小;在入口位置相对一侧壁面附近区域及底部区域分布着气含率极低区.

(4)入口位置的改变对曝气池内气泡数密度的大小及分布没有十分显著的影响，曝气池内气泡破碎形成的小气泡数量个数明显多于聚并形成的大气泡个数.

致谢:论文模拟过程中华东理工大学陈彩霞教授、江苏大学詹水清老师、上海交通大学的梁晓飞硕士对数值模拟提供了很好的建议及帮助，在此一并感谢!

［1］ 李孟，李向阳，王宏智，等.鼓泡塔气液两相流不同曳力模型的数值模拟［J］.过程工程学报，2015，15(2):181-189.LI Meng， LI Xiangyang， WANG Hongzhi， et al.Numerical simulation of gas-liquid two-phase flow in a bubble column with various drag models［J］.Chinese Journal of Process Engineering，2015，15(2):181-189.

［2］ 金丹.塔式曝气池内气液两相流动数值模拟及参数影响的研究［D］.哈尔滨:哈尔滨工业大学，2006.

［3］ 魏文礼，赵小军，刘玉玲.气泡浮力羽流动力特性三维数值模拟研究［J］.西北农林科技大学学报自然科学版，2014，42(5):229-234.WEIWenli， ZHAO Xiaojun， LIU Yulin.Threedimensional numerical simulation of dynamic characteristicsofair bubble buoyancy plume［J］.Journal of Northwest A ＆ F University:Nat.Sci.Ed.，2014，42(5):229-234.

［4］ ALI B A，KUMAR C S，PUSHPAVANAM S.Analysis of liquid circulation in a rectangular tank with a gas source at a corner［J］.Chemical Engineering Journal，2008，144(3):442-452.

［5］ 肖柏青，张法星，戎贵文.气泡尺寸对曝气池内气液两相流数值模拟的影响［J］.中国环境科学，2012，32(11):2006-2010.XIAO Baiqing， ZHANG Faxing， RONG Guiwen.Influence of the bubble size on numerical simulation of the gas- liquid flow in aeration tanks［J］. China Environmental Science，2012，32(11):2006-2010.

［6］ LEHR F，MILLIES M，MEWES D.Bubble-size distributions and flow fields in bubble columns［J］.AIChE Journal，2002，48(11):2426-2443.

［7］ 邢楚填.鼓泡床反应器实验研究及CFD-PBM耦合模型数值模拟［D］.北京:清华大学，2014.

［8］ WANG T，WANG J.Numerical simulations of gasliquid mass transfer in bubble columns with a CFDPBM coupled model［J］. Chemical Engineering Science，2007，62(24):7107-7118.

［9］ GUPTA A，ROY S.Euler-Euler simulation of bubbly flow in a rectangular bubble column:experimental validation with radioactive particle tracking［J］.Chemical Engineering Journal，2013，225(1/2/3/4/5/6):818-836.

［10］ 徐琰，董海峰，田肖，等.鼓泡塔中离子液体-空气两相流的 CFD-PBM耦合模拟［J］.化工学报，2011，62(10):2699-2706.XU Yan ，DONG Haifeng，TIAN Xiao，et al.CFDPBM coupled simulation of ionic liquid-air two-phase flow in bubble column［J］.CIESC Journal，2011，62(10):2699-2706.

［11］ LI G，YANG X，DAI G.CFD simulation of effects of the configuration of gas distributors on gas-liquid flow and mixing in a bubble column［J］. Chemical Engineering Science，2009，64(24):5104-5116.

［12］ D AZ M E， IRANZOA， CUADRA D， etal.Numerical simulation of the gas-liquid flow in a laboratory scale bubble column:Influence of bubble size distribution and non-drag forces［J］.Chemical Engineering Journal，2008，139(2):363-379.

［13］ CLIFT R，GRACE J R，WEBER M E.Bubbles，drops，and particles［M］.London:Academic Press，1978:97-146.

［14］ LIANG X F，PAN H，SU Y H，et al.CFD-PBM approach with modified drag model for the gas-liquid flow in a bubble column［J］.Chemical Engineering Research＆ Design，2016，112:88-102.

［15］ LABORDE-BOUTET C，LARACHI F，DROMARD N，et al.CFD simulation of bubble column flows:investigations on turbulence models in RANS approach［J］.Chemical Engineering Science，2009，64(21):4399-4413.

［16］ HAGESAETHER L，JAKOBSEN H A，SVENDSEN H F.A model for turbulent binary breakup of dispersed fluid particles［J］.ChemicalEngineeringScience，2002，57(16):3251-3267.

［17］ HAGESAETHER L，JAKOBSEN H A，KAI H，et al.A coalescence and breakup module for implementation in CFD-codes［J］. Computer Aided Chemical Engineering，2000，8(0):367-372.

［18］ LUO H，SVENDSEN H F.Theoretical model for drop and bubble breakup in turbulentdispersions［J］.Chemical Engineering Science，1996，66(5):766-776.

［19］ LUO H.Coalescence，breakup and liquid circulation in bubble column reactors［D］.Trondheim:Norwegian Institute of Technology，1993.

［20］ 罗玮.曝气池中气液两相流PIV实验研究及数值模拟［D］.西安:西安理工大学，2006.

［21］ 肖浩飞.曝气池内气液两相流CFD数值模拟［D］.上海:东华大学，2010.

［22］ 徐礼嘉，陈彩霞，夏梓洪，等.鼓泡塔内气泡羽流周期性摆动的数值模拟［J］.化学工程，2012，40(9):48-51.XU Lijia， CHEN Caixia， XIA Zihong， etal.Numericalsimulation of bubble plume periodic oscillation in bubble column［J］. Chemical Engineering，2012，40(9):48-51.

［23］ BUWA V V，RANADE V V.Dynamics of gas-liquid flow in a rectangular bubble column:experiments and single/multi-group CFD simulations［J］.Chemical Engineering Science，2002，57(22/23):4715-4736.

［24］ 肖柏青，张法星，刘春艳，等.曝气池内气泡羽流附壁效应的试验研究［J］.水力发电学报，2012，31(4):104-107.XIAO Baiqing，ZHANG Faxing，LIU Chunyan，et al.Coanda effect of bubble plume in aeration tanks［J］.Journal of Hydroelectric Engineering，2012，31(4):104-107.