基于EMMS方法的鼓泡塔反应器CFD及群平衡模拟

王珏，杨宁

（1中国科学院过程工程研究所多相复杂系统国家重点实验室，北京 100190；2中国科学院大学，北京 100049）

基于EMMS方法的鼓泡塔反应器CFD及群平衡模拟

王珏1,2，杨宁1

（1中国科学院过程工程研究所多相复杂系统国家重点实验室，北京 100190；2中国科学院大学，北京 100049）

能量最小多尺度 （energy-minimization multi-scale，EMMS） 方法已经被应用于气液体系中群平衡（population balance model，PBM）模型的改进。EMMS模型可计算气泡破碎聚并过程的能量，进而获得聚并速率的修正因子。应用这一模型对高气速鼓泡塔进行了模拟计算，并进一步对比了均一尺径模型、CFD-PBM 模型以及CFD-PBM-EMMS模型的模拟结果与实验数据。结果表明，在高表观气速条件下，基于EMMS方法的群平衡模型可以更加准确地预测鼓泡塔中不同高度的气泡尺径分布和轴向液速，同时提高了对整体气含率和局部气含率的模拟准确性。在表观气速为0.16 m·s-1和0.25 m·s-1时，CFD-PBM-EMMS模型对气泡尺径分布的预测精度更高，同时整体气含率模拟的相对误差下降为5%和15%，局部气含率模拟平均相对误差下降为8%和17%。

计算流体力学；群平衡模型；鼓泡塔；气含率；气泡尺径分布

引 言

鼓泡塔反应器具有结构简单、操作简便、良好的传热传质效率等优点，被广泛应用于化学工程、生物工程、石油工程等领域[1]。随着计算机技术的发展，计算流体力学（computational fluid dynamics，CFD）方法已经成为鼓泡塔反应器模拟和优化的重要工具。对于鼓泡塔中的气液两相流动常见的模拟方法有欧拉-欧拉方法和欧拉-拉格朗日方法。其中，欧拉-拉格朗日方法需要对每一个气泡进行跟踪，对计算资源要求高，适用于气含率较小的反应器模拟[2]；而欧拉-欧拉方法假定离散相与连续相之间可相互渗透，不直接追踪每一个气泡的运动，从而显著地降低了计算量，被广泛地应用于实验室规模和工业规模的鼓泡塔模拟之中。

欧拉-欧拉双流体模型中，气液相间作用力模型对模拟结果至关重要[3]。相间作用力包括曳力、升力、虚拟质量力、湍流扩散力等，其中曳力对气液系统的影响占主导地位[4-7]。在曳力模型中，气泡尺寸常被作为一个重要的输入参数，直接影响模拟结果。模型设置中，通常依据实验或经验将气泡直径设定为定值。然而实际系统中，由于气泡之间存在破碎和聚并的现象，导致气泡直径在反应器内呈现动态时空非均匀分布特征。针对此特征，群平衡模型（population balance model，PBM）提供了一个解决方案，对气泡数密度建立传输方程，考虑气泡聚并和破碎对气泡数密度变化的影响。因此，群平衡模型耦合计算流体力学方法（CFD-PBM）对鼓泡塔进行模拟已经成为近年来的热点之一[8]。

群平衡模型中，聚并核函数和破碎核函数直接决定了气泡粒径分布规律的模拟结果[9]。聚并过程通常简化为二元聚并的情况[10]，并基于物理现象推导和实验测量提出了不同的聚并核函数[11-12]。Chen等[13]发现现有的破碎聚并函数并不能使气泡群的聚并速率和破碎速率在鼓泡塔稳定鼓泡区达到动态平衡，这与实验现象不符。Wang等[14]提出，当气泡之间的距离比气泡湍动路径长时，气泡的实际碰撞频率要比Prince模型小，导致模型高估了气泡间的聚并效率。Bhole等[15]认为现有聚并模型在计算碰撞频率时使用对应尺度涡的速度代替气泡速度，而气泡与湍动涡之间是存在速度差的，所以聚并速率应该乘以一个小于1的系数。其他学者也提出，需要使用小于1的修正因子来修正聚并函数以减小聚并速率[16-18]。肖颀[19]和 Yang等[20]率先提出了基于 EMMS方法的群平衡（PBM）模型模拟鼓泡塔的CFD-PBM-EMMS模型，模拟结果表明了该模型可以合理地模拟实验室规模中低气速鼓泡塔中的气泡粒径分布。Qin等[21]采用了EMMS-PBM的方法对液-液体系的转子定子中的液滴尺径分布进行了修正，提高了原有模型的准确度。本文针对高表观气速的气液鼓泡塔，应用CFD-PBM-EMMS方法对高气速气液鼓泡塔进行了模拟，通过与实验数据对比验证了该方法的可行性。

1 数学模型和数值模拟

本文中的三维CFD模型以欧拉-欧拉模型为基础，采用重整化群（RNG）k-ε湍流模型和Schiller-Naumann曳力模型。同时，耦合群平衡模型（PBM）计算气泡尺径分布，其中描述气泡聚并和破碎过程分别采用 Luo聚并模型和 Luo &Svendsen破碎模型，并使用离散法（class method，CM）求解群平衡方程。

1.1 控制方程

气液系统模拟采用欧拉-欧拉模型，欧拉-欧拉模型又称双流体模型（two-fluid model，TFM）。模拟中把液体设置为连续相，气体设置为离散相，并且不考虑气体的压缩和相间质量、热量的传递。控制方程可以表示为以下形式。

质量守恒方程

动量守恒方程

式中，α、ρ、u、t分别表示相含率、密度、速度矢量和时间；P表示压力；μeff表示流体有效黏度；Mi指相间动量交换项。

式中，各项对应的有效黏度μeff是液相分子黏度μL和湍流黏度μt之和。其中湍流黏度为

式中，Cμ=0.0845，湍动能 k和湍动能耗散率ε可从湍流模型耦合计算得到。

1.2 湍流模型

湍流模型按照文献[22-23]中的选择，采用重整化群（RNG）k-ε混合相模型（Mixture），模拟湍动状态下的气液系统[13]。

式中，C1ε=1.42，C2ε=1.68。

在此 Mixture模型中，气液两相共用一套方程，方程中采用混合密度ρm和速度um，计算公式如下

湍流动能产生项Gk,m用以下公式计算得到

式中，η≡Sk/ε(S 是张量的积)，η0=4.38，β=0.012。

1.3 相间作用力

根据Sanyal等[4]和Van Baten等[5]的研究结果，本文的模拟计算只考虑曳力，忽略其他相间作用力，曳力系数选择 Schiller-Naumann[24]模型模拟气液体系相间作用力[11]。

其中，曳力系数为

式中，CD是气泡群曳力系数，CD,0为单气泡曳力系数，CD=CD,0(1-αg)p，本模拟中p=1。

1.4 群平衡模型

群平衡模型对气泡在时间和空间中进行的变化进行追踪，可以很好地预测气泡尺径的分布。在群平衡模型中气泡尺径的变化与聚并速率和破碎速率密切相关。而沸腾、冷凝还有气泡在边界上的生成和减小在本文中暂不考虑。

式中，n是气泡的数密度；c，b函数分别为聚并核函数和破碎核函数；β是母气泡破碎时子气泡的分布函数；式（13）右端的第1项和第2项分别表示由聚并引起的气泡产生和消亡；式（13）右端的第3项和第4项则分别表示由破碎引起的气泡产生和消亡。

1.4.1 破碎模型 本文选用Luo & Svendsen[25]的破碎模型

其中

1.4.2 聚并模型 Luo[26]的聚并模型为

式中，聚并速率c等于碰撞频率ϖc与聚并效率Pc的积。

其中

1.5 聚并速率修正因子

能量最小多尺度（EMMS）模型由Li[27]提出，现已广泛应用于气固流化床反应器中。Yang等[28]基于EMMS模型提出了适用于气-液体系的DBS模型，并且应用于鼓泡塔、搅拌槽等气液反应器。对于鼓泡塔体系，DBS模型将气相简化为“大气泡相”和“小气泡相”两相，通过大小气泡的直径(dL和dS)、气含率(αL和αS)和表观速度(Ug,L和 Ug,S)来描述气液系统的状态。DBS模型的主要模型方程[28]如下。

质量守恒方程

力平衡方程

对气液体系的能量关系的描述如下。

总能量

介尺度能耗

微尺度能耗

式中，CD,p为固体颗粒曳力系数；CD,b为气泡曳力系数；ϖ为涡与气泡碰撞频率；αb为气含率；Pb为气泡破碎概率；fv为破碎前后子气泡与母气泡体积比；cf为气泡破碎时表面积增加量，等于σ为表面张力。

在鼓泡塔的气液体系中，存在气相控制和液相控制两种机制。当气相控制为主导时，液相的湍动能趋于最小（Nturb→ min），而液相控制时，气泡在湍流中表面振荡能耗最小（Nsurf→ min）。在这两种控制机制都发挥作用的复杂流型之下由两种机制协调控制，DBS模型稳定条件为

DBS模型在当前的气液体系研究中的新进展见文献[29-30]。Yang等[31]在DBS的基础上提出了气液体系DBS曳力模型，在CFD模拟鼓泡塔的气含率方面优势明显，并推广至搅拌槽中的气含率预测[32]。最近，肖颀[19]和Yang等[20]将气液体系内介尺度能量的描述和 PBM 模型相结合，提出了对聚并模型的修正方法，如图1所示。气液体系在稳定状态下，微尺度能量趋于最小，同时不同尺径之间的气泡破碎和聚并处于动态平衡。此时在聚并相前加入抑制聚并的修正因子 C，在使用 PBM 求解的介尺度能耗 Nbreak_PBM和 DBS求解的介尺度能耗Nbreak_DBS相等时，就得到由DBS模型进行修正的PBM修正因子。

破碎核函数：Luo & Svendsen[25]；聚并核函数：Luo 对应修正因子关系式[26]。

图1 聚并核函数修正因子求解流程[20]Fig. 1 Flow chart for calculation of coalescence corrector[20]

2 建模对象和方法

2.1 结构模型设置和网格无关性验证

本文用 GAMBIT建立鼓泡塔的三维物理模型并且划分网格，在进气口采用面进气的方式，对分布器上部鼓泡区域模拟计算。鼓泡塔塔径390 mm，塔高2000 mm，液面高为1000 mm。为了获得合适尺寸的网格，本文对粗（coarse）、中（medium）、细（fine）3种不同尺寸的网格尺寸进行了分析。图2表示了三维网格划分方式，并给出了中等尺寸网格的俯视图和侧视图。表1为分布器不同的划分尺寸。模型分布器与McClure等[33]的保持一致，中间部分为分布器管道壁面，两边为进气区域。进气区域采取面进气的模型设置，各个区域网格结构如图2所示。

图2 网格结构Fig. 2 Mesh configuration

表1 网格无关性验证模型网格尺寸及数目Table 1 Mesh size and cells adopted for grid independence test

图3给出了3种不同尺寸网格模拟得到的不同高度截面的平均气含率分布。3种不同尺寸的网格均能大致反映出鼓泡塔中气含率的轴向变化趋势。粗网格的网格数目为17万个，因为网格数目偏低，导致计算时产生了一定的误差。而在网格数目大于36万个之后，网格数目对模拟结果的影响较小。因此，考虑到计算速度和模拟结果准确性，本文后续模拟工作均采用中等尺寸的网格。

图3 不同网格尺寸模拟所得局部气含率分布Fig.3 Local gas holdup at different column height of different mesh size

2.2 数学模型求解设置：

模拟计算基于Ansys Fluent 15.0软件，采用非稳态模型方程。鼓泡塔入口处采用gas-velocity边界条件，出口处使用pressure-outlet边界条件。壁面使用无滑移边界。为了保持计算的稳定性，计算初始时使用时间步长Δt=0.0001 s，计算1000步，增大时间步长至Δt=0.001 s，计算1000步至1.1 s，最终使用Δt=0.005 s稳定计算至120 s。压力速度耦合方式采用SIMPLE格式离散，动量方程、体积分数以及湍流方程均采用一阶迎风格式离散，松弛因子使用0.75的默认设置。

图4为鼓泡塔中不同高度下的Sauter平均尺径（d32）和7 mm的气泡尺径体积分数随着时间变化的模拟结果。由图4（a）可以看出，在40 s之后，不同高度截面的 d32基本相同，认为此时鼓泡塔已经进入较为稳定的状态。由图中可以观察到，气泡群的平均尺径沿着塔高逐渐增加。由图4（b）可以看出，在40 s之后，7 mm气泡的体积分数呈现出周期性的变化趋势，流场数据均取40～120 s的数据进行时间平均。

本文中群平衡模型的求解方法采用离散法，离散法具有概念简单、容易理解的特点，能够给出直观的颗粒粒级分布。

离散法处理各个气泡粒级

图4 气泡尺寸随时间的变化趋势Fig. 4 Instantaneous variation of bubble size

同时气泡个数和气含率之间存在

式中，Ni、αi、Vi分别表示第i个气泡粒级的气泡数密度函数、气含率、体积。

TFM 中用于描述气泡变化的气含率αi项和PBM 方程中描述不同气泡浓度的 Ni的关系即可获得。结合McClure等[33]实验数据，将气泡尺径范围设定为1，3，5，……，39 mm 20个粒级，入口的初始颗粒粒径使用Bhole等[15]使用的5 mm。

3 模拟结果与讨论

本文使用3种模型进行模拟，如表2所示。case 1采用均一气泡尺径的双流体（TFM）模型；case 2采用群平衡（PBM）模型耦合双流体模型，群平衡（PBM）核函数采用 Luo聚并函数和 Luo &Svendsen破碎模型；case 3也采用群平衡（PBM）模型耦合双流体（TFM）模型，不同的是加入了EMMS方法对聚并模型修正。经计算，表观气速为0.16 m·s-1时修正因子为 0.1230，表观气速为 0.25 m·s-1时计算出的修正因子为0.1051。3种模型中均使用Schiller-Naumann曳力模型，入口处气泡的尺径均为5 mm。

表2 不同模型算例的设置Table 2 Different model settings

3.1 气含率

图5为不同模型鼓泡塔整体气含率模拟和实验结果的对比。由图可知高气速下传统的双流体模型明显高估气含率，0.16 m·s-1和0.25 m·s-1表观气速下相对误差分别为40%和80%。可能原因在于，高气速条件下，气泡受到液相湍动的影响较为明显，气泡在不同高度段有着不同的尺径分布，气泡的实际尺径不能使用入口气泡的尺径近似处理。Schiller-Naumann曳力模型在计算气含率方面，也受到尺径的影响。Xu等[34]发现，在表观气速为0.1 m·s-1时，气泡尺径设置为2.5、4.5、6.5 mm均一气泡尺径时，采用Schiller-Naumann曳力模型的条件下发现，气含率会随着气泡尺径的增加而降低。case 1假定了均一的气泡尺径，没有描述气泡经过聚并破碎后平均尺径沿塔高增大的现象，导致模拟结果中的气含率偏高。相比较而言 PBM 模型描述了气泡上升过程中的破碎与聚并现象，得到了气泡的尺径分布。未经修正的群平衡模型略低估了气含率，在0.16 m·s-1和0.25 m·s-1的表观气速下整体气含率模拟数据和实验结果之间的相对误差为20%和15%，而经过修正之后的群平衡模型因提供了更为准确的气泡尺径信息，从而预测的气含率模拟结果和实验数据之间的相对误差为5%和15%。

图5 整体气含率模拟数据和实验结果[33]的对比Fig. 5 Comparison between different models and experiments[33]

图6 不同模型整体气含率云图Fig. 6 Contour plots for different models

图6为不同模型的时均气含率分布云图。所有case均能反映出鼓泡塔中气含率分布的结构特征。进气口处受分布器受非进气区域影响，塔中心气含率偏低。在鼓泡塔中上高度处，case 1的气含率分布较为均匀，这与使用均一气泡尺径相关。在Usup=0.25 m·s-1的情况下，较高的气含率导致模拟液面高于鼓泡塔顶部。同时在case 2、case 3中，可以明显地看到鼓泡塔中的气含率在轴向和径向不同位置的不均匀分布；其中case 2由于气含率偏低，导致鼓泡液位高度与未鼓泡前的静液位高度相比差距并不明显。

图7 局部气含率模拟数据和实验结果[33]的对比Fig.7 Comparison of local gas holdup between different models and experiments[33]

图7是不同高度处的局部气含率沿径向的分布。McClure等[33]选取了鼓泡塔中段和接近相界面处的两个截面记录了实验数据。case 1、case 2和case 3的模拟数据同实验结果相比，气含率的径向分布结果均呈现出中间高，边壁低的分布规律。分布器非均匀分布结构使得轴心处的气含率略有偏低，3个case也都捕捉到了这一现象。与整体气含率相同的是，局部气含率的结果也表现出了case 1明显高估了气含率，case 2、case 3在误差允许的范围内都可以较为可靠地反映鼓泡塔z=550 mm处的局部气含率径向分布。case 2的整体气含率低于case 3。

3.2 气泡大小分布

3.2.1 修正因子对气泡尺径分布的影响 图8给出了 case 2、case 3 在 0.16 m·s-1和 0.25 m·s-1表观气速下鼓泡塔z=550 mm和z=1050 mm处气泡尺径分布的模拟值和实验值的对比，其中case 2中的气泡峰值位于39 mm的粒级。小于39 mm的其他粒级的气泡，分布较为平缓。分析原因：① case 2 中的气泡粒级受到1～39 mm粒级区间设置的限制和聚并核函数对聚并速率的高估，使粒径的峰值为 39 mm，而且从37 mm到39 mm的粒级间陡然上升，由于原始聚并模型高估了聚并速率[11]，模拟得到的粒径分布显示出大气泡较多的趋势；② 从曲线整体分析，1～37 mm的气泡尺径分布较为平缓，与实验结果相比，说明了模拟中气泡的聚并被高估了。case 3尝试使用修正因子抑制聚并函数，修正气泡聚并速率。修正后的模型中，不同尺径气泡的体积分布范围集中于5～15 mm，9 mm的气泡具有最高的体积分数。两种不同的模型对比，修正后的群平衡模型预测的粒径分布更加符合实际。实验中体积分数最高的气泡粒级为 5 mm，与模型预测结果仍有一定差距。McClure等[33]分析实验结果认为，在H/d=2、H/d=3截面上，气泡尺径分布在实验允许的误差范围内是相同的，而带有修正因子的模拟结果也反映了这一点。因此，在H/d≥2截面上，气泡的聚并和破碎达到了动态平衡。

图8 气泡尺径分布模拟数据和实验结果[33]的对比Fig.8 Comparison of bubble size distribution between different models and experiments[33]

3.2.2 Sauter平均直径 除模拟出高气速下气泡尺径在不同位置处的动态变化外，群平衡模型可进一步获得气泡的Sauter平均尺径（d32）信息，此平均直径将用于双流体模型中的耦合计算。图9给出了case 1中的双流体模型采用均一气泡的设定，气泡直径保持了5 mm的初始值不变。未经修正的PBM和 CFD耦合模型计算出的气泡在不同径向位置上的d32在8～18 mm范围变化，而修正后的PBM模型由于聚并受到了抑制，平均尺径减小，H/d=2到H/d=3径向位置上平均粒径在7～8 mm范围变化。

3.3 轴向液速

图10 中可以看出，不同气泡尺径的条件下，3个不同case模拟的轴向液速相差不大，说明气泡尺径对轴向液速的模拟的影响不大。分布器结构导致鼓泡塔轴心处液速较为平缓，与塔中心处的气含率较低（图7）情况相同。在z=1050 mm处，模拟的轴向液速都较实际值偏低，与徐琰等[35]对于离子液体的气液两相流的模拟类似，其模拟结果也出现接近塔顶处液速减小的趋势。z=550 mm处，不同模型之间差别很小，对于液速的预测大体能符合实验中获得的真实流速。

4 结 论

准确地模拟气泡尺径分布，描述鼓泡塔中不同空间位置的气泡大小，对高气速下气液两相流体系有着更为重要的意义。本文采用CFD-PBM模型对高气速下的气液鼓泡塔体系进行了模拟，应用EMMS/DBS模型计算出不同气速下的聚并速率修正因子，研究了不同群平衡模型下的气液两相流的流场信息，获得以下结论。

（1）高气速条件下的鼓泡塔中发生强烈的湍动，气泡的破碎和聚并导致入口处气泡尺径和塔中稳定鼓泡段的气泡尺径有着很大的区别，此时采用均一气泡尺径模型不能准确模拟塔内的所有流场信息。基于EMMS方法修正群平衡模型，考虑了气泡的聚并和破碎，获得了具有稳定尺径分布的气泡群信息，捕捉到气泡群尺径由入口处的均一粒径(5 mm)上升至平均粒径7 mm的动态变化。

图9 Sauter平均直径在不同表观气速下的径向分布Fig.9 Radial profiles of Sauter mean diameter at different superficial velocities

（2）传统的 PBM聚并模型需要经过修正才能准确地描述反应器中的实际气泡分布。一般情况下传统的聚并核函数模型均高估了聚并速率。使用EMMS-PBM 耦合方法计算出的修正因子有效地抑制了聚并，得到单峰分布的气泡尺径曲线。与未修正前的曲线相比，模拟预测更加接近实验结果。

（3）CFD-PBM-EMMS群平衡模型可以模拟气泡粒径从入口到出口处不同高度上的平均直径的变化，并且在稳定鼓泡区给出了较为准确的气泡平均尺径。在表观气速0.16 m·s-1和0.25 m·s-1时，使用CFD-PBM-EMMS模型模拟结果中，整体气含率的相对误差为 5%和 15%，局部气含率的平均相对误差为8%和17%；未经修正的CFD-PBM模型，整体气含率的相对误差为20%和15%，局部气含率的相对误差为23%和20%。CFD-PBM-EMMS方法优化了对高气速下气液流场的模拟。

（4）模拟结果表明，EMMS-PBM 耦合方法作为一种改进群平衡方程的思路是可靠的，此方法仍有进一步的改进空间。

符 号 说 明

b——破碎核函数，s-1

C——聚并修正因子

c——聚并核函数，m3·s-1

cf——破碎前后气泡表面积增加率

D——气泡直径，m

dL, dS——分别为大、小气泡直径，m

d32——Sauter平均直径，m

FD——气泡群曳力，N·m-3

k ——湍动能，m2·s-2

M ——相间动量交换项，kg·m-2·s-2

N——数密度函数

Nbreak——破碎聚并能耗，m2·s-3

图10 轴向液速模拟结果[33]与实验数据的对比Fig.10 Comparison of upward liquid velocity between different models and experiment[33]

Nsurf——气泡表面能，m2·s-3

NT——悬浮输送能，W·kg-1

Nturb——液相黏性耗散能，m2·s-3

P ——压强，Pa

Re ——Reynolds数

Ug,L, Ug,S——分别为大、小气泡表观气速，m·s-1

Usup——表观气速，m·s-1

V ——气泡体积，m3

We ——Weber数

α,αg——分别为相含率、气含率

αL,αS——分别为大、小气泡气含率

ε ——湍动能耗散率，m2·s-3

λ ——湍流涡尺度，m

μ ——动力学黏度，Pa·s

ξ ——涡尺寸和气泡直径比值

ξij——两聚并气泡直径之比

ρ ——密度，kg·m-3

σ——表面张力，N·m-1

ω——涡与气泡碰撞频率，s-1

ωc——两气泡碰撞频率，m3·s-1

[1] KANTARCI N, BORAK F, ULGEN K O. Bubble column reactors[J].Process Biochem., 2005, 40(7)： 2263-2283.

[2] SIMCIK M, MOTA A, RUZICKA M C, et al. CFD simulation and experimental measurement of gas holdup and liquid interstitial velocity in internal loop airlift reactor[J]. Chem. Eng. Sci., 2011,66(14)： 3268-3279.

[3] HIBIKI T, ISHII M. Lift force in bubbly flow systems[J]. Chem. Eng.Sci., 2007, 62(22)： 6457-6474.

[4] SANYAL J, VASQUEZ S, ROY S, et al. Numerical simulation of gas-liquid dynamics in cylindrical bubble column reactors[J]. Chem.Eng. Sci., 1999, 54(21)： 5071-5083.

[5] VAN BATEN J M, KRISHNA R. CFD simulations of a bubble column operating in the homogeneous and heterogeneous flow regimes[J]. Chem. Eng. Technol., 2002, 25(11)： 1081-1086.

[6] KRISHNAR, VAN BATEN J M, URSEANU M I. Three-phase Eulerian simulations of bubble column reactors operating in the churn-turbulent regime： a scale up strategy[J]. Chem. Eng. Sci., 2000,55(16)： 3275-3286.

[7] SCARGIALI F, D’ORAZIO A, GRISAFI F, et al. Modelling and simulation of gas-liquid hydrodynamics in mechanically stirred tanks[J]. Chem. Eng. Res. Des., 2007, 85(A5)： 637-646.

[8] WANG T F, WANG J F, JIN Y. A novel theoretical breakup kernel function for bubbles/droplets in a turbulent flow[J]. Chem. Eng. Sci.,2003, 58(20)： 4629-4637.

[9] 覃成鹏, 杨宁. 多相分散体系中气泡/液滴聚并和破碎的群平衡模拟[J]. 化学进展, 2016, (8)： 1207-1223.QIN C P, YANG N. Population balance modeling of breakage and coalescence of dispersed bubbles of droplets in multiphase systems[J].Progress in Chemistry, 2016, (8)： 1207-1233.

[10] RAMKRISHNA D. Population Balances： Theory and Applications to Particulate Systems in Engineering[M]. San Diego： Academic Press,2000.

[11] CHEN P, SANYAL J, DUDUKOVIĆ M P. Numerical simulation of bubble columns flows： effect of different breakup and coalescence closures[J]. Chem. Eng. Sci., 2005, 60(4)： 1085-1101.

[12] RIBEIRO C P, MEWES D. On the effect of liquid temperature upon bubble coalescence[J]. Chem. Eng. Sci., 2006, 61(17)： 5704-5716.

[13] CHEN P, DUDUKOVIC M P, SANYAL J. Three-dimensional simulation of bubble column flows with bubble coalescence and breakup[J]. AIChE J., 2005, 51(3)： 696-712.

[14] WANG T F, WANG J F, JIN Y. Theoretical prediction of flow regime transition in bubble columns by the population balance model[J].Chem. Eng. Sci., 2005, 60(22)： 6199-6209.

[15] BHOLE M R, JOSHI J B, RAMKRISHNA D. CFD simulation of bubble columns incorporating population balance modeling[J]. Chem.Eng. Sci., 2008, 63(8)： 2267-2282.

[16] YAO W, MOREL C. Volumetric interfacial area prediction in upward bubbly two-phase flow[J]. Int. J. Heat Mass Trans., 2004, 47(2)：307-328.

[17] MITRE J F, TAKAHASHI R S M, RIBEIRO C P, et al. Analysis of breakage and coalescence models for bubble columns[J]. Chem. Eng.Sci., 2010, 65(23)： 6089-6100.

[18] MUKIN R V. Modeling of bubble coalescence and break-up in turbulent bubbly flow[J]. Int. J. Multiphas Flow, 2014, 62： 52-66.

[19] 肖颀. 搅拌釜及鼓泡塔内气液两相流的多尺度模型[D]. 西安： 西安交通大学, 2015.XIAO Q. Multi-scale modeling of the gas-liquid flow in mixing tanks and bubble columns[D]. Xi’an： Xi’an Jiaotong University, 2015.

[20] YANG N, XIAO Q. A mesoscale approach for population balance modeling of bubble size distribution in bubble column reactors[J].Chem. Eng. Sci., 2017, DOI： 10.1016/j.ces.2017. 01.026.

[21] QIN C, CHEN C, XIAO Q, et al. CFD-PBM simulation of droplets size distribution in rotor-stator mixing devices[J]. Chem. Eng. Sci.,2016, 155： 16-26.

[22] 许婷婷. DBS曳力模型与湍流模型对气液CFD模拟的影响[D]. 北京： 中国科学院大学, 2014.XU T T. Influence of DBS drag model and turbulence models on CFD simulation of gas-liquid flow[D]. Beijing： University of Chinese Academy of Sciences, 2014.

[23] XU T, JIANG X, YANG N, et al. CFD simulation of internal-loop airlift reactor using EMMS drag model[J]. Particuology, 2015, 19：124-132.

[24] SCHILLER L, NAUMANN Z. A drag coefficient correlation[J]. Z.Ver. Deutsch. Ing., 1935, 77(1)： 318-320.

[25] LUO H, SVENDSEN H F. Theoretical model for drop and bubble breakup in turbulent dispersions[J]. AIChE J., 1996, 66(5)：766-776.

[26] LUO H. Coalescence, breakup and liquid circulation in bubble column reactors[D]. Trondheim, Norway： the Norwegian Institute of Technology, 1993.

[27] LI J. Particle-fluid Two-Phase Flow： the Energy-Minimization Multi-Scale Method[M]. Beijing： Metallurgical Industry Press, 1994.

[28] YANG N, CHEN J H, ZHAO H, et al. Explorations on the multi-scale flow structure and stability condition in bubble columns[J]. Chem.Eng. Sci., 2007, 62(24)： 6978-6991.

[29] YANG N. Chapter five-mesoscale transport phenomena and mechanisms in gas-liquid reaction systems[J]. Advances in Chemical Engineering, 2015, 46： 245-280.

[30] XIAO Q, YANG N, LI J. Stability-constrained multi-fluid CFD models for gas-liquid flow in bubble columns[J]. Chem. Eng. Sci.,2013, 100： 279-292.

[31] YANG N, WU Z, CHEN J, et al. Multi-scale analysis of gas-liquid interaction and CFD simulation of gas-liquid flow in bubble columns[J]. Chem. Eng. Sci., 2011, 66(14)： 3212-3222.

[32] 肖颀, 杨宁. 基于EMMS模型的搅拌釜内气液两相流数值模拟[J].化工学报, 2016, 67(7)： 2732-2739.XIAO Q, YANG N. Numerical simulation of gas-liquid flow in stirred tanks based on EMMS model[J]. CIESC Journal, 2016, 67(7)：2732-2739.

[33] MCCLURE D D, NORRIS H, KAVANAGH J M, et al. Validation of a computationally efficient computational fluid dynamics (CFD)model for industrial bubble column bioreactors[J]. Ind. Eng. Chem.Res., 2014, 53(37)： 14526-14543.

[34] XU L, YUAN B, NI H, et al. Numerical simulation of bubble column flows in churn-turbulent regime： comparison of bubble size models[J].Ind. Eng. Chem. Res., 2013, 52(20)： 6794-6802.

[35] 徐琰, 董海峰, 田肖, 等. 鼓泡塔中离子液体-空气两相流的CFD-PBM耦合模拟[J]. 化工学报, 2011, 62(10)： 2699-2706.XU Y, DONG H F, TIAN X, et al. CFD-PBM coupled simulation of ionic liquid-air two-phase flow in bubble column[J]. CIESC Journal,2011, 62(10)： 2699-2706.

CFD-PBM simulation with EMMS correctors for bubble column reactors

WANG Jue1,2, YANG Ning1
(1State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences,Beijing 100190, China;2University of Chinese Academy of Sciences, Beijing 100049, China)

The energy-minimization multi-scale (EMMS) model has been introduced to improve the population balance modeling (PBM) of gas-liquid flows. The energy for bubble breakup and coalescence can be obtained from the EMMS model and then used to derive a correction factor for the coalescence rate. This new model is applied in this study to simulate the bubble columns of high flow rates. Simulations using the three different models, namely, the constant-bubble-size model, the CFD-PBM model and the CFD-PBM-EMMS model, are compared with experimental data. The simulation of CFD-PBM-EMMS gives better prediction for bubble size distribution and liquid axial velocity at different heights as well as the overall and local gas holdup. The relative error of global gas holdup reduces to 5% or 15%, and the mean relative error of local gas holdup reduces to 8% or 17% for 0.16 m·s-1or 0.25 m·s-1of superficial gas velocity.

CFD; population balance modeling; bubble column; gas holdup; bubble size distribution

date:2017-01-08.

Prof. YANG Ning, nyang@ipe.ac.cn

supported by the National Natural Science Foundation of China (91434121,91634203).

TQ 021.1

A

0438—1157（2017）07—2667—11

10.11949/j.issn.0438-1157.20170025

2017-01-08收到初稿，2017-03-13收到修改稿。

联系人：杨宁。

王珏（1991—），男，硕士研究生。

国家自然科学基金项目（91434121, 91634203）。