基于EMMS模型的气固鼓泡床的模拟及气泡特性的分析

吴迎亚，彭丽，高金森，蓝兴英

（中国石油大学（北京）重质油国家重点实验室，北京 102249）

基于EMMS模型的气固鼓泡床的模拟及气泡特性的分析

吴迎亚，彭丽，高金森，蓝兴英

（中国石油大学（北京）重质油国家重点实验室，北京 102249）

基于EMMS曳力模型，采用双流体的方法对气固鼓泡床内的气固流动特性进行模拟，建立基于图像处理气泡特性的分析方法，重点研究了不同表观气速下气泡在床层内分布特性，包括气泡平均当量直径、气泡速度和气泡球形度的轴向分布，以及气泡的生命周期。研究结果表明，小气泡多集中在床层底部和壁面区域，而大气泡多集中在床层中间区域。随着表观气速的增加，床层高度不断增加，气泡的球形度降低，气泡的大小、出现频率、上升速度以及生命周期均增加；然而，当表观气速增大到一定程度，继续增加气速对气泡的上升速度影响不大。

气固鼓泡床；多相流；计算流体力学；TFM；EMMS；气泡

引　言

气固鼓泡床具有良好的混合、传热及传质效果，在过程工业中有着广泛的应用［1］。计算流体力学方法（CFD）［2］被广泛应用于研究鼓泡床内气固流动行为的研究，其中欧拉-欧拉双流体模型（TFM）的应用最为广泛。双流体模型将气相和颗粒相均看成连续介质，满足质量、动量和能量守恒方程，其中颗粒相的压力、黏度等物性采用颗粒动理学（KTGF）描述，气固相间相互作用采用曳力模型描述［3-4］。而对于FCC这类Geladart A类颗粒而言，在流动过程中极易发生团聚现象，造成体系内较大颗粒聚团的生成，从而使得曳力降低。因此，采用Gidaspow等［5-8］经典曳力模型难以准确地模拟FCC颗粒在鼓泡床内的气固流动现象，高估了床层的膨胀率［9-12］。为此，很多研究者对曳力模型进行修正，其中应用最广泛的是能量最小化多尺度（EMMS）曳力模型［13］。研究者采用双流体模型耦合 EMMS曳力模型模拟气固循环流化床内气固流动特性［14-16］，研究结果发现，采用EMMS曳力模型能较好地模拟出气固循环流化床内的气固流动特性［17-21］。Shi等［22］进一步将应用于气固循环流化床模拟中的 EMMS曳力模型拓展并应用于气固鼓泡床的数值模拟，发展了新的基于气泡修正的EMMS曳力模型，并能较为准确地预测鼓泡床内的床层膨胀率等特性。此外，气固鼓泡床内气泡和颗粒之间存在着强烈的作用，加强了鼓泡床内气固的接触效率，保证了其良好的传热、传质以及反应特性。众多研究者通过实验和模拟研究发现［23-26］，鼓泡床内气泡的分布和颗粒的循环方式是相互影响的：气泡和颗粒的返混决定了气固的反应特性；气泡的大小、形态及运动规律决定了气固传热传质特性。因此，研究鼓泡床内气泡分布和气泡运动规律对进一步认识鼓泡床内的气固流动、传热及传质机制有着重要的意义。而前人的研究［23-26］工作中一般只关注了气泡的平均性质，而对气泡在床层内的瞬时分布特性以及气泡的动态变化过程缺乏深入研究。本文首先基于EMMS曳力模型分析了鼓泡床内的气固流动特性，并在此基础上建立了识别气泡和分析气泡的方法，分析了鼓泡床内气泡的平均性质、气泡在床层内的整体分布特性以及气泡的生命周期。

1　模拟对象

图1　模拟对象的几何结构［9］Fig. 1　Geometry structure of simulation domain［9］

表1　模拟条件Table 1　Simulation conditions

本文的模拟工况为M cKeen等［9］建立的三维柱形气固鼓泡床冷态实验。模拟对象如图1所示，具体的实验条件及相关参数如表1所示。在室温下进行了4个表观气速（0.05、0.10、0.15、0.20 m·s-1）的实验。虽然三维模拟能更加准确地反映出气固体系的复杂流动特性，但是考虑到计算能力的限制，二维模拟作为一种较为准确和经济的简化方法，被众多研究者采用［21，27-28］。因此，本文在模拟计算过程中采用二维模拟，计算时间步长为0.001 s，气体流经整个床层需要的时间约为5 s，为了取得较稳定的床层，共模拟了鼓泡床流动时间为25 s。在下文后续分析中，除床层以及气泡的瞬时特性外，床层内颗粒体积分数分布和气泡的平均性质的分析均取5～25 s之间的平均数据。进行网格无关性研究后，发现轴向和径向上采用5 mm × 5 mm的均匀网格满足计算需求。因此，本模拟选择的网格尺寸为 5 mm × 5 mm的均匀网格。

表2　模型的控制方程Table 2　Governing equations

2　双流体模型

本文采用欧拉双流体方法（TFM）结合颗粒动理学理论（KTGF）对气固鼓泡床进行模拟研究，模型的控制方程组如表2所示，关于模型的本构方程及相关表达式见文献［29］。

表3　不同表观气速下的曳力模型Table 3　Correction factor for different superficial gas velocity

图2　气泡的表征方法Fig. 2　Method of bubble characterization

气固之间的曳力模型采用EMMS曳力模型［22］。EMMS曳力模型是在Wen-Yu曳力模型的基础上修正得到。其曳力系数表达式如下

其中f为修正因子［22］，不同气速下修正因子如表3所示。

图3　不同表观气速下床层内颗粒体积分数瞬时分布云图（10 s）Fig. 3　Instantaneous profiles of solids volume fraction at different superficial gas velocity at 10 s

3　结果与讨论

3.1气泡分析方法

为了研究鼓泡床内的气泡特征，本研究基于模拟得到的鼓泡床内气含率分布云图建立识别和分析气泡的方法，如图2所示。具体步骤如下：① 通过选择一个合适的气含率阈值识别和提取气泡［23-26］（一般定义鼓泡床中固含率小于 0.2的区域为气泡）；② 根据所选定的阈值将模拟得到的鼓泡床内气含率分布云图二值化；③ 通过边缘检测算法，确定气泡的边缘；④ 计算出气泡的中心位置和面积；⑤ 对相邻帧的气泡进行对比分析，求取气泡运动的速度。此外，为了进一步表征气泡的特征，对气泡进行球形近似，通过气泡的当量直径和球形度表征气泡的特性。其定义如下

其中C为球形度，A为气泡的面积，S为气泡的周长。C越接近于1，表示气泡接近完美球形。

3.2气固流动特性

图4　不同表观气速下时均颗粒体积分数的轴向分布（5～25 s）Fig. 4　Axial profiles of time-averaged solids volume fraction under different superficial gas velocity at 5—25 s

图3为模拟得到的不同表观气速下床层颗粒体积分数在模拟时间为10 s的瞬时分布云图。图4为不同表观气速下床层时均颗粒体积分数的轴向分布图。由图3和图4可知，随着表观气速的增大，床层高度不断增加；在床层高度为0～50 cm内，床层内的颗粒体积分数随表观气速的增加不断降低；同时沿床层高度方向，气泡的大小也明显变大。在较低气速下（0.05 m·s-1），床层界面清晰；而在较高气速下（0.20 m·s-1），床层界面有较大波动。表4对比了不同表观气速下床层高度的模拟结果和实验数据［9］。由表 4可知，模拟计算值略高于实验值，这可能是由于模拟计算中忽略了多孔型分布板和电容层析成像装置的影响，造成了模拟值偏高。

图5　不同表观气速下时均颗粒体积分数的径向分布（5～25 s）Fig. 5　Radial profiles of time-averaged solids volume fraction under different superficial gas velocity at 5—25 s

表4　不同表观气速下床层高度的模拟结果和实验结果对比Table 4　Height in simulation and experiment under different superficial gas velocity

图6　不同表观气速下时均颗粒速度的径向分布（5～25 s）Fig. 6　Radial profiles of time-averaged solids velocity under different superficial gas velocity at 5—25 s

图5和图6分别为不同表观气速下床层时均颗粒体积分数和颗粒速度在床高H = 12、28和42 cm处的径向分布。从图5可以看出，在床层的不同高度截面处，颗粒体积分数沿径向均呈现中间低边壁高的不均匀分布规律，由中心到边壁，分布曲线变陡，颗粒体积分数变大，在壁面处达到最大值。在较低气速下，颗粒体积分数较高，而且在径向上分布较为均匀；在较高气速下，中心处呈现较低的波谷，且气速越大波谷越低。同一气速下，在床层的不同高度截面处，颗粒体积分数分布情况变化不大。由图6可以看出，在床层的不同高度截面处，均呈现出床层中心区域颗粒速度大，边壁区域颗粒速度小的分布趋势，且随着气速的增大，床层中心区域的颗粒速度不断增加。

3.3气泡特征

气泡在鼓泡床中有极其重要的作用，气泡的生成、聚并以及破碎决定了气固鼓泡床的传热和反应效率。在准确模拟鼓泡床内气固流动特性的基础上，采用基于图像处理气泡的分析方法（见3.1节），进一步分析了气固鼓泡床内的气泡特性。

图7为不同表观气速下，鼓泡床内气泡平均当量直径的轴向分布。由图7可知，随着床层高度的增加，气泡逐渐增大。在较低的气速下（0.05 m·s-1和 0.10 m·s-1），当床层高度增加到一定程度，气泡大小趋于稳定，停止生长，直至在气固相界面发生破碎。在较高的气速下（0.15 m·s-1和 0.20 m·s-1），由于气体表观速度较大，气固之间的湍动激烈，随着床层高度的增加，气泡的大小总体趋于增大的趋势，但是由于激烈的湍动，造成一定的波动。图8是不同表观气速下，气泡平均上升速度的轴向分布。总体上来说，随着床层高度的增加，气泡上升速度增加。然而，当表观气速增大到一定程度，继续增大表观气速对气泡上升速度的影响不大。图9是不同表观气速下，鼓泡床内气泡的球形度的轴向分布。由图9可以发现，在较低气速下，气泡较小，球形度较大；在较高气速下，气泡较大，球形度较小。这是由于气固之间强烈的湍动容易造成气泡聚并和破碎，导致气泡的球形度减小。在床层底部区域，气泡处于生成和生长期，气泡球形度较大，几乎近似于球形；而在床层的中上部，由于气泡的聚并和破碎，引起气泡拉长扭曲，最终导致气泡的球形度减小。

图7　不同表观气速下气泡平均当量直径的轴向分布（5～25 s）Fig. 7　Axial profiles of time-averaged bubble equivalent diameter under different superficial gas velocity at 5—25 s

图8　不同表观气速下气泡平均速度的轴向分布（5～25 s）Fig. 8　Axial profiles of time-averaged bubble velocity under different superficial gas velocity at 5—25 s

图9　不同表观气速下气泡球形度的轴向分布（5～25 s）Fig. 9　Axial profiles of bubble sphere coefficient under different superficial gas velocity at 5—25 s

为了更加详细地分析气泡在鼓泡床内分布特性，统计了20～21 s时间内不同表观气速下床层内的所有气泡分布情况，如图10所示。由图10可知，在较低气速下，床层内的气泡数目较少，气泡较小，大部分气泡趋于在床层中上部生成。随着气速的增大，床层内气泡的数目不断增加，气泡的大小也在逐渐增加。总体来看，小气泡多集中在床层底部区域，而大气泡多集中在床层的中上部；从径向分布来看，大气泡多集中在床层中间区域，小气泡多集中在床层壁面区域。且随着表观气速的增加，气泡的运动路径呈弯曲上升趋势，中心区域的气泡在上升过程中，不断发生着气泡的增长和聚并，气泡的大小不断增大。在鼓泡床内，生命周期较长的气泡对床层内颗粒的混合、传热、传质和反应性能有着较大的影响。进一步统计了床层10～25 s内生命周期大于0.5 s的气泡出现的次数，并统计出床层气泡热区分布，如图11所示。由图11可知，在较低的气速下，床层内仅存在少量的长周期气泡。随着气速的增大，长周期气泡的数目不断增加，且多集中在床层中部。

图10　20～21 s内不同表观气速下床层内的气泡分布Fig. 10　Distributions of bubble under different superficial gas velocity at 20—21 s

图11　不同表观气速床层内气泡生命周期大于0.5 s的气泡热区图Fig. 11　Bubble hot spots at different superficial gas velocity for bubble life cycle of over 0.5 s

4　结　论

本文首先采用欧拉双流体方法结合 EMMS曳力模型对气固鼓泡床内的气固流动特性进行模拟研究，并建立了识别和分析鼓泡床内气泡特性的方法，分析了不同表观气速下气泡沿床层轴向方向上的分布特性，包括气泡平均当量直径、气泡速度以及气泡球形度的轴向分布，进一步详细地分析了不同气速下气泡在鼓泡床内的分布特性，得到了以下结论：随着表观气速的增大，床层高度不断增加，气泡的球形度降低；在床层高度为0～50 cm内，床层内的固含率随表观气速的增大不断降低；同时沿床层高度方向，气泡的大小也明显变大。在鼓泡床底部，随着表观气速的增大，气泡产生的频率显著增加。在较低表观气速下（0.05 m·s-1），床层界面清晰；而在较高表观气速下（0.20 m·s-1），床层界面出现较大波动；在较高的表观气速下（0.15 m·s-1和0.20 m·s-1），由于气固之间存在强烈的湍动作用，随着床层高度的增加，气泡上升速度增加。但当表观气速增大一定程度，继续增加气速对气泡上升速度的影响不大；小气泡多集中在床层底部和壁面区域，而大气泡多集中在床层中间区域。且随着表观气速的增大，气泡的运动路径呈弯曲上升趋势，中心区域的气泡不断生长和聚并，导致气泡不断增大；在较低的表观气速下，床层内仅存在少量的长周期气泡。随着表观气速的增大，长周期气泡的数目不断增加，且多集中在床层中部。本文建立了识别和分析鼓泡床内气泡特性的方法，仍缺乏对气泡的聚并和破碎过程的追踪，后期的工作仍需进一步完善。

符号说明

A ——气泡直径，m

C ——气泡球形度

CD——曳力系数，kg·m-3·s-1

f ——修正因子

I ——单位张量

kΘs——颗粒脉动动能扩散系数

S ——气泡周长，m

u ——速度/表观气速，m·s-1

V ——速度矢量，m·s-1

β ——曳力，kg·m2·s-1

γΘs——碰撞耗散能，kg·m-1·s-3

ε ——空隙率

Θs——颗粒温度，m2·s-2

ρ——密度，kg·m-3

τ ——应力

下角标

g——气相

p——颗粒

s——固相/颗粒相

References

［1］ 金涌， 祝京旭， 俞芷青. 流态化工程原理［M］. 北京: 清华大学出版社， 2001.

JIN Y， ZHU J X， YU Z Q. Fluidization Engineering Principles［M］. Beijing: Tsinghua University Press， 2001.

［2］ ZHU H P， ZHOU Z Y， YANG R Y， et al. Discrete particle simulation of particulate systems: a review of major applications and findings ［J］. Chem ical Engineering Science， 2008， 63（23）: 5728-5770.

［3］ ZOU L M， GUO Y C， CHAN C K. Cluster-based drag coefficient model for simulating gas-solid flow in a fast-fluidized bed［J］. Chem ical Engineering Science， 2008， 63（4）: 1052-1061.

［4］ MAKKAW I Y， WRIGHT P， OCONE R. The effect of friction and inter-particle cohesive forces on the hydrodynam ics of gas-solid flow: a comparative analysis of theoretical predictions and experiments［J］. Powder Technology， 2006， 163: 69-79.

［5］ GIDASPOW D. Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions［M］. Boston: American Press， 1994.

［6］ WEN C， YU Y. Mechanics of fluidization［J］. Chem ical Engineering Progress Symposium Series， 1966， 62: 100-111.

［7］ ERGUN S. Fluid flow through packed columns［J］. Chem ical Engineering Progress， 1952， 48: 89-94.

［8］ TURTON R， LEVENSPIEL O. A short note on the drag correlation for spheres［J］. Powder Technology， 1986， 47: 83-86.

［9］ MCKEEN T， PUGSLEY T. Simulation and experimental validation of a freely bubbling bed of FCC catalyst［J］. Powder Technology，2003， 129（1）: 139-152.

［10］ FERSCHNEIDER G， MEGE P. Eulerian simulation of dense phase fluidized beds［J］. Rev. Inst. Fr. Pe’t.， 1996， 51（2）: 301-307.

［11］ KRISHNA R， VAN BATEN J M. Using CFD for scaling up gas-solid bubbling fluidized bed reactors w ith Geldart A powders［J］. Chem ical Engineering Journal， 2001， 82（2）: 247-257.

［12］ BAYLE J， MEGE P， GAUTHIER T. Dispersion of bubble flow properties in a turbulent FCC fluidized bed［M］// Fluidization X，Engineering Foundation. New York， 2001: 125-132.

［13］ WANG W， LI J H. Simulation of gas-solid two-phase flow by a multi-scale CFD approach-extension of the EMMS model to the sub-grid level［J］. Chem ical Engineering Science， 2007， 62（1/2）: 208-231.

［14］ LI F. Investigations on the turbulent gas-solid two-phase interactions in fluidized desulfurization process［D］. Beijing: Tsinghua University，2009.

［15］ LU B. EMMS-based meso-scale model and its application in simulating gas-solid two-phase flows［D］. Beijing: Chinese Academy of Sciences， 2009.

［16］ ZHOU Q， WANG J， LI J. Three-dimensional simulation of dense suspension upflow regime in high-density CFB risers w ith EMMS-based two-fluid model［J］. Chem ical Engineering Science，2014， 107: 206-217.

［17］ JIRDILOK V， GIDSSPOW D， DAMRONGLERD S. et al. Kinetic theory based CFD simulation of turbulent fluidization of FCC particles in a riser［J］. Chemical Engineering Science， 2006， 61（17）: 5544-5559.

［18］ QI H， LI F， XI B， et al. Modeling of drag with the Eulerian approach and EMMS theory for heterogeneous dense gas-solid two-phase flow［J］. Chemical Engineering Science， 2007， 62（6）: 1670-1681.

［19］ NIKOLOPOULOS A， ATSONIOS K， NIKOLOPOULOS N， et al. An advanced EMMS scheme for the prediction of drag coefficient under a 1.2 MWth CFBC isothermal flow（Ⅱ）: Numerical implementation［J］. Chemical Engineering Science， 2010， 65（13）: 4089-4099.

［20］ BENYAHIA S， SUNDARESAN S. Do we need sub-grid scale corrections for both continuum and discrete gas-particles flow models［J］. Powder Technology， 2011， 220: 2-6.

［21］ LU B， WANG W， LI J. Eulerian simulation of gas-solid flows with particles of Geldart groups A， B and D using EMMS-based meso-scale model［J］. Chemical Engineering Science， 2011， 66（20）: 4624-4635.

［22］ SHI Z， WANG W， LI J. A bubble-based EMMS model for gas-solid bubbling fluidization［J］. Chemical Engineering Science， 2011， 66（22）: 5541-5555.

［23］ VERMA V， PADDING J T， DEEN N G. Effect of bed size on hydrodynamics in 3D gas-solid fluidized beds［J］. AIChE Journal，2015， 61（5）: 1492-1506.

［24］ EVGENIDIS S P， KARAPANTSIOS T D. Effect of bubble size on void fraction fluctuations in dispersed bubble flows［J］. International Journal of Multiphase Flow， 2015， 75: 163-173.

［25］ LU Y， HUANG J， ZHENG P. Flow structure and bubble dynamics in supercritical water fluidized bed and gas fluidized bed: a comparative study［J］. International Journal of Multiphase Flow， 2015， 73: 130-141.

［26］ ASEGEHEGN T W， SCHREIBER M， KRAUTZ H J. Investigation of bubble behavior in fluidized beds with and without immersed horizontal tubes using a digital image analysis technique［J］. Powder Technology， 2011， 210（3）: 248-260.

［27］ LV X， LI H， ZHU Q. Simulation of gas-solid flow in 2D/3D bubbling fluidized beds by combining the two-fluid model with structure-based drag model［J］. Chemical Engineering Journal， 2014， 236: 149-157.

［28］ CLOETE S， JOHANSEN S T， AMINI S. Investigation into the effect of simulating a 3D cylindrical fluidized bed reactor on a 2D plane［J］. Powder Technology， 2013， 239: 21-35.

［29］ GIDASPOW D， JUNG J， SINGH R K. Hydrodynamics of fluidization using kinetic theory: an emerging paradigm: 2002 Flour-Daniel lecture［J］. Powder Technology， 2004， 148: 123-141.

Numerical simulation of gas-solid bubbling bed and bubble characteristics based on EMMS drag model

WU Yingya， PENG Li， GAO Jinsen， LAN Xingying
（State Key Laboratory of Heavy Oil Processing， China University of Petroleum， Beijing 102249， China）

The gas-solid flow in gas-solid bubbling fluidized bed was simulated by a combined approach of computational fluid dynam ics （CFD） and the two fluid model （TFM）， based on the modified EMMS drag model. A new image processing method derived from image calibration was applied to study bubble characteristics w ith focus on overall and radial distributions of bubbles as well as bubble profiles in terms of averaged equivalent diameter， rising velocity， degree of sphericity， and bubble lifetime at various superficial gas velocity. The results showed that smaller bubbles most located at the bottom of beds and wall region whereas larger bubbles most located at the central region of beds. With increase of superficial gas velocity， the bed height were gradually increased and bubble profiles of the averaged equivalent diameter， the frequency of occurrence， the rising velocity and the lifetime were all increased except that the degree of sphericity was decreased. Nevertheless， once the superficial gas velocity increased to a certain level， its increase would no longer have any significant impact on the rising velocity.

gas-solid bubbling fluidized bed； multiphase flow； computational fluid dynam ics ； TFM； EMMS；bubble

date: 2016-04-07.

LAN Xingying， lanxy@cup.edu.cn

supported by the National Basic Research Program of China （2012CB215003） and the New Century Excellent Talents Program in University（NCET-13-1027）.

TQ 021.1

A

0438—1157（2016）08—3259—09

10.11949/j.issn.0438-1157.20160441

2016-04-07收到初稿，2016-06-08收到修改稿。

联系人：蓝兴英。第一作者：吴迎亚（1990—），男，博士研究生。

国家重点基础研究发展计划项目（2012CB215003）；教育部新世纪优秀人才支持计划项目（NCET-13-1027）。