基于EMMS的循环流化床流域研究

陈恺成，田于杰，李飞,3，吴昊，王维

（1中国科学院过程工程研究所多相复杂系统国家重点实验室，北京100190；2中国科学院大学中丹学院，北京100049；3中国科学院绿色过程制造创新研究院，北京100190；4Department of Chemical and Biochemical Engineering，Technical University of Denmark，2800 Kongens Lyngby，Denmark）

引 言

气固流态化系统广泛应用于化工、冶金、能源、环境等领域，具有气固接触效率高、传递性能好、操作灵活等优点[1]。然而，气固流态化体系是一个典型的非平衡系统，气体和颗粒在微尺度上的多种非线性非平衡相互作用导致介尺度上气泡、聚团等非均匀结构的产生，并进一步影响宏尺度的流动行为。随着操作条件、物性和床层几何结构的不同，流化床呈现出各种复杂的流域，如鼓泡流化、湍动流化、快速流化、气力输送等[2-4]。鼓泡流化与湍动流化的操作气速较低，颗粒通量较小，又被称为传统流化床，常用于甲醇制烯烃等过程[5]；以快速流化为代表的循环流化床，常用于重质油催化裂化[6]、丁醇直接氧化制顺酐[7]等过程。

不同的流域不仅对应不同的操作条件，而且在气固接触效率、颗粒停留时间分布等方面存在明显差异[8]，因而直接影响着流化床反应器的设计与选型。同时，在流化床的逐级放大过程中，伴随反应装置与操作条件的改变，可能发生的流域转变也会影响传热、反应等性能[9]。因此，如何定性乃至定量地预测流域及其转变，对于指导工业生产及放大具有重要意义。

自20 世纪40 年代起，关于流域转变的理论建模[10-12]与实验研究[13-16]从未停止。床层压降测量[17]、颗粒成像测速[18]、正电子发射颗粒追踪[19]、X 射线断层扫描[20]、声波发射[21]等多种测量手段都曾被用于表征各流域的流动特征，压降波动的变化[17]、流动结构的转变[22-23]、操作变量间函数关系的改变[4,24]等被作为判定流域转变的标准，界定流域转变的经验关联式也有很多[16,25-29]。这催生了多种流域划分的结果。然而，一方面，由于不同测量手段得到的结果不同[30]、不同研究者关注的流动特征不同、不同物性系统与几何装置中的流动情况不同[31-32]，这些流域划分往往具有较强的经验性。流域划分的结果并不统一，甚至有所冲突，反而为流域转变的研究带来了更多的困惑。尤其是在对循环流化床（快速流化床）所涉及流域的界定上，不同研究者之间存在很大分歧。在Grace[2]总结过去研究而绘制的流域图上，传统床、循环床、输送床之间边界模糊、互相交叠，便是这一点的体现。另一方面，实验研究费时费力、成本较高、可以提取的流动信息有限，难以作为统一流域研究的工具。

随着数值模拟技术的发展与计算性能的提高，计算流体力学（CFD）模拟逐渐成为研究多相流动的有力手段。CFD 模拟可以展示更丰富的流动信息与流场细节，且相比实验研究更为经济、快捷，使得系统研究多种流域成为可能。在各种模拟方法中，双流体模型（two-fluid model,TFM）将气体与颗粒作为连续的两相处理，以其较低的计算成本而被工程计算所关注。但经典TFM 中使用的气固相间曳力基于均匀流动的假设，其计算结果与气固流化床中的真实情况存在较大偏差[33]。随着介尺度理论[34-35]的兴起，多尺度CFD 方法[36]考虑了流化床中存在的非均匀结构（如气泡、聚团），并被证明能够更准确地预测流化床中的流动情况。基于稳态的能量最小多尺度（EMMS）曳力模型考虑了非均匀结构对气固相间作用力的影响，已被证明适用于鼓泡床、湍动床、快速床等主要流域的模拟[37-41]，并能合理地捕捉到“噎塞”等流域转变现象[37,42]。

本文对前述流域研究工作进行总结，对循环流化床中不同的流域划分进行比较。以基于稳态EMMS的曳力模型对气固循环流化系统进行模拟研究，尝试对循环床流域转变所涉及的现象给出一定的解释。

1 流域划分的争议

总结文献中关于流域研究的工作，发现影响流域转变的因素大致包括：（1）操作条件；（2）体系的物性参数；（3）反应器的几何结构。基于这些影响因素的不同组合，流域图大致可分为以下几类：第一类：仅考虑操作条件的流域图，如Yerushalmi等[4]、Bi等[22]绘制的流域图；第二类：考虑物性参数的流域图，如Geldart[31]对流化床颗粒的分类、Grace[2]以无量纲气速与Archimedes 数为坐标绘制的流域图；第三类：考虑床层几何结构的流域图，如Shaul 等[32]考虑反应器高径比对流域转变影响的流域图。各类流域图使用的基本变量及其所划分的流域汇总见表1。

表观气速小于输送气速(Utr)时的流态化通常被称为传统流化床，流域划分较为清晰。这一阶段典型的流域有鼓泡流化和湍动流化。一般认为，床层压降波动最大时对应的气速为鼓泡床向湍动床的转变气速[30]。对于Geldart A 类颗粒，在鼓泡床出现之前存在床层均匀膨胀的流域；而对于床层直径较小的鼓泡床，当气泡直径大于反应器直径的60%时会出现“节涌流化”。

表1 文献中的流域图总结Table 1 Regime diagrams in literature

表观气速大于输送气速时的流域被统称为循环流化床，其流域转变更为复杂。对于循环床向相邻流域的转变，一些研究者认为这种转变仅与气速有关，如Yerushalmi 等[4]、Avidan 等[24]、Horio[3]。尽管这些研究者也采取了两种变量绘制流域图，然而其流域划分的标准却是“一维”的：当气速达到输送气速时，系统由湍动床过渡到快速床；当气速达到输运气速(Utp)时，系统由快速床过渡到稀相输送。按这种方式进行流域划分，转变气速通常仅取决于系统的物性，虽然简单易求，但却难以解释颗粒浓度或颗粒通量对流域转变的影响。

认识到颗粒浓度或通量对流域转变的影响之后，后续研究者进行的流域划分更多是“二维”的：流域转变的边界不仅取决于气速，同时也受颗粒通量的影响。持这种观点的有Bai 等[13]、Bi 等[22]、Sun等[43]。图1汇总展示了这类流域图，表2具体罗列了这些研究者所划分的流域及其特征。

Bai 等[13]定义的快速床具有明显的轴/径向颗粒浓度分布，在“局部与整体上都有很强的非均匀性”，其操作下限（UFD,即由输送床向快速床的过渡）对应上稀下浓结构的出现，操作上限(UTF)伴随着上部稀相区大量聚团破碎、垮塌。Bi 等[22]定义的快速床具有颗粒循环通量大、轴径向颗粒浓度分布明显等特点，其操作下限为出现聚团时的气速(VCA)，操作上限为出现节涌时的气速(VCC)。Sun 等[43]定义的快速床同样具有明显的径向非均匀结构，其操作下限(Utp)来自于对文献中经验关联式的平均，而操作上限(Utr)对应床层崩塌、通量锐减等现象。此外，该研究还将快速床流域细分为两个子流域：固含率小于10%的低密度循环流化床(LDCFB)与固含率大于10%的高密度循环流化床(HDCFB)，二者具有不同的轴/径向流动结构，分别适用于循环床燃烧与催化裂化这两种主要的循环床工业应用。不同研究者基于大量实验数据拟合了快速床的边界，这些边界及其适用的物性和床层几何尺寸范围见表3。可以看出，各研究者对于快速床流动特征的描述相近，然而划分快速床与毗邻流域的经验关联式却并不统一，流域图的构成形式也不相同。

表2 循环床的流域细分及特征描述Table 2 Regime classifications and characteristics in CFBs

图1 文献中不同循环床流域划分的比较Fig.1 Comparison of regime classifications of CFBs in literatures

在Wang 等[42]之前的工作中，曾基于经典EMMS模型给出了循环床本征流域图的划分，将流域划分为稀相输送、密相流动及二者之间的噎塞型转变，如图2 所示。其中，噎塞曲线上的流态具有明显的径向环核结构和轴向上稀下浓颗粒分布。然而，在噎塞曲线的下方，表2 中所列的一些非均匀特征就已出现，具体表现为浓相和稀相的空隙率不相等[46]；而在噎塞曲线的上方，实际操作的通量显然难以无限增大。因此，与完全不考虑几何结构影响的本征流域图不同，如何确定现实操作图中的快速床上下限以及流动结构的特征，值得进一步关注。

图2 基于EMMS的循环床本征流域图[42]Fig.2 EMMS-based intrinsic regime diagram of CFBs [42]

综上所述，流域划分的争议主要集中在循环流化床之中，而文献中对快速床流域如何界定未有定论。本文将以基于稳态EMMS的曳力模型模拟不同操作条件下的算例，致力于确定快速床的操作范围。

表3 界定快速床操作上下限的经验关联式Table 3 Empirical correlations to outline the fast bed region

2 模型方法与模拟设置

本节介绍模拟算例的设置及使用的模型方法。相较之前的工作，模拟了多个气速、更大浓度范围内的算例。

模拟采用双流体模型，对气固两相的运动均采用连续介质模型进行描述。其中，颗粒相的应力采用颗粒动理论进行封闭，气-固相间曳力采用基于稳态EMMS 的曳力模型，模型的推导及验证工作见文献[41]，非均匀结构因子Hd的关联式见附录1。双流体模型的控制方程、本构关系以及用于封闭模型的关系式见附录2。

由于前述流域图采用的物性体系、反应器几何构造都不相同，针对每一系统进行模拟不仅会耗费大量时间，而且得到的结果不具有可比性。因此，本文采用Yerushalmi 等[4]的装置与物性系统，对不同气速、不同初始颗粒浓度下的算例进行模拟，着重研究操作条件对于流域转变的影响。模拟的几何构体、物性参数和操作条件分别见图3与表4。

模拟使用图3 所示的二维构体，在商业软件ANSYS Fluent®中进行。由于基于稳态EMMS 的曳力模型具有较好的网格无关性[41]，本模拟采用较粗的网格(15 mm×85 mm)以节省计算成本。所用网格在Gambit®2.4生成。模拟开始前，颗粒被设置以一定的初始浓度(εs0)均匀分散于整个反应器中，空气从反应器入口均匀进入，带动反应器内的颗粒向上流动，直到从位于反应器上部的出口流出。反应器出口压力为大气压，出口处的颗粒被全部返料进入底部入口进行物料循环。每个算例的模拟总时长为15 s，最后5 s 用于时均统计。详细的模拟设置见表5。

图3 模拟采用的几何构体Fig.3 Schematic drawing of riser in simulations

表4 模拟系统的物性参数与操作条件Table 4 Material properties and operating conditions in simulations

表5 模拟设置Table 5 Simulation settings

3 结果与讨论

3.1 快速床的操作下限

图4 展示了不同操作条件下床层高度为0.5 m处的固含率与竖直方向上颗粒速度和气体速度的时均径向分布。可以看出，当εs0=0.01 时，各气速下的固含率径向分布较为均匀。随初始颗粒浓度的增大，固含率逐渐呈现出中间稀、两边浓的径向分布，且气速越低、初始颗粒浓度越高，这种径向分布越明显。竖直方向的颗粒速度upy和气体速度ugy呈现出中心高、两边低的径向分布，随初始颗粒浓度的增大，气体和颗粒的速度分布都更加明显。当εs0=0.01 时，各气速下的upy在不同径向位置处均为正值，随着初始颗粒浓度增加，当εs0=0.03，气速为2.2、3.0、3.4 m/s 时，提升管壁面附近处出现了颗粒向下运动的现象，且气速越小，颗粒向下运动的速度越大；而当气速为4.0 m/s 时，在εs0=0.05 时才出现颗粒向下运动的现象。

颗粒浓度较低时所展现的近似均匀的流化状态常被认为是稀相输送的特征；随颗粒浓度的增加，颗粒开始在壁面附近形成聚团，由此而导致的径向非均匀结构、颗粒回流等现象则被认为是典型的快速床特征。若以壁面处颗粒速度由正转负作为快速床操作开始的标志，可将不同气速下该点对应的颗粒通量求出，得到快速床的操作下限，绘于图5中，并与Bai等[13]、Bi等[22]、Sun等[43]提出的经验关联式比较。可以看出，划定快速床操作下限的各种经验关联式之间并不吻合。这可能是由于：（1）物性与反应器几何结构的不同对流域转变造成了影响，因此对在不同装置中开展的、使用不同物性体系的实验结果进行直接拟合得到的关联式，在对某一特定体系进行预测时会产生误差；（2）不同研究者判定流域转变的依据并不完全相同，Bai等[13]的判断依据为轴向非均匀性的出现，Bi 等[22]的判断依据为出现颗粒聚团，Sun等[43]虽然没有提出相应的判断标准，但是对快速流化特征的描述与上文两个研究者相似，即具有明显的轴径向颗粒浓度分布。理论上讲出现颗粒聚团后，轴径向非均匀结构也会随之出现，但类似“S 形”轴向分布、“环-核型”径向分布等特征的出现难以精确定量判断，具有很强的经验性。本文将边壁处颗粒回流作为判定快速床操作下限的标准，更易观察测量也更为明确，据此得到的流域转变分界线落于Bi等[22]、Sun等[43]所提出的经验关联式之间，与Bai 等[13]定义的边界接近。在Gs-Ug流域图上，这条线的右下方表示稀相输送、左上方表示快速床和密相流态化。

3.2 快速床的操作上限

图6 展示了气速为4 m/s，初始颗粒浓度分别为0.25、0.28 与0.30 时出口处颗粒通量随时间的变化。可以看出，在模拟刚开始时，三组结果都出现了极大的颗粒通量，这可能是由于模拟设置的初始状态为颗粒均匀分布于提升管内，因而导致当气体到达提升管出口附近时，大量悬浮颗粒被带出。在一段时间（约10 s）之后，εs0=0.25、εs0=0.28 两组的颗粒通量逐步稳定，在某一值附近上下波动，而εs0=0.30 的一组则呈现出完全不同的流动状态：颗粒通量始终剧烈波动，通量的最高值可达1500 kg/(m2·s)及以上，而最低值则接近于0。对应床层浓度的瞬时快照可以看出，此时在提升管中出现了类似“节涌”的现象：颗粒团与气体呈栓状向上运动，当颗粒栓经过出口时，监控得到的颗粒通量极大；而气栓经过时，颗粒通量则极小。节涌发生时，虽然平均颗粒循环通量很大，但因床层压降大幅涨落、气固混合程度小，操作极其不稳定。节涌的出现可以视作快速床操作的终止[47]。

图4 不同操作条件下的径向颗粒浓度与颗粒速度、气体速度分布（床高0.5 m处）Fig.4 Radial distribution of solids holdups,particle velocities and gas velocities(axial position:0.5 m)

图5 快速床的下边界Fig.5 Lower boundary of fast bed

图6 较大颗粒浓度下通量随时间的变化Fig.6 Variation of solids flux with time at higher solids holdup

其他气速下的模拟结果呈现出相似的变化规律：当初始颗粒浓度较低时，颗粒通量随模拟时间逐步稳定，而当颗粒浓度大于某一值后会出现节涌。对于某一气速，选取节涌尚未出现时初始颗粒浓度最大的算例，将其颗粒通量的时均值作为该气速下快速床操作的最大通量，即快速床操作的上边界，如图7 中实线所示。实线的右下方表示快速床和稀相输送，左上方由于节涌的存在无法稳定操作。与Bai 等[13]、Bi 等[22]、Sun 等[43]提出的经验关联式比较，发现模拟得到的边界位于不同的经验关联式之间。同样，这可能是因为物性条件与床层几何结构对流域转变有影响，因而对不同体系、不同装置中的实验结果进行平均，所得的关联式在对某一特定体系进行预测时会产生偏差。

至此，确定了快速床的操作范围，如图8中红色区域所示。相较于之前的经验关联式，该区域的划定基于CFD 模拟中的流动特征，原则上在对不同物性系统与几何结构的循环床模拟中仍能出现类似结果。

图7 快速床的上边界Fig.7 Upper boundary of fast bed

图8 快速床的操作范围Fig.8 Operating range of fast bed

3.3 快速床的流域细分

图9展示了不同气速下时均颗粒通量随初始颗粒浓度的变化。可以看到，在颗粒浓度较低与较高时，颗粒通量均随颗粒浓度的增加而增加。而在二者之间，随颗粒浓度增加，通量基本不变，在Gs-εs0图中呈“平台型”区域，该区域被称为“噎塞型转变”[42]。

图10 展示了气速为4 m/s 时，快速床流域内不同初始颗粒浓度下的轴/径向颗粒浓度分布。从10(a)中可以看出，在这个操作区间内，提升管内均存在着明显的“上稀下浓”分布，不同的是，随颗粒浓度的增加，底部的密相区域不断上移，逐步靠近提升管的出口。从图10(b)中可以看出，当初始颗粒浓度较低时，径向颗粒浓度呈“环-核型”分布，而随着颗粒浓度增加，径向颗粒浓度逐步变为类似凹形抛物线的分布。这可能是由于随着颗粒浓度增加，边壁处累积的颗粒聚团逐渐变浓、变大，并向管中心处发展，当聚团逐步拓展至提升管中心时，便形成了类似抛物线的径向浓度分布。

图9 不同气速下颗粒通量随初始颗粒浓度的变化Fig.9 Variation of solids flux with mean solids holdup at different gas velocities

以最小流化速度Umf和它与颗粒密度的乘积(ρsUmf)对图8 进行无量纲化，并将图9 中平台型区域对应的颗粒通量绘于该图中，可以得到图11中代表噎塞型转变的蓝色实线。在这条曲线的下部，颗粒通量与浓度均较小、提升管底部密相区的高度较低、颗粒浓度呈环-核型径向分布；在这条曲线的上部，颗粒浓度大，且其径向分布类似凹形抛物线。对照表2可以看出，噎塞曲线下方与上方的区域，分别与Sun 等[43]描述的低密度循环流化床(LDCFB)与高密度循环流化床(HDCFB)类似。值得一提的是，Sun 等[43]以εs=0.1 为界区分LDCFB 与HDCFB，然而图11 中的噎塞转变曲线却对应通量不变时的一系列浓度范围。如按文献[43]的方式，将图中蓝色和红色的区域分别定义为LDCFB 与HDCFB，则灰色区域为稀相输送；蓝色和红色的区域共同构成快速床流域；而空白区域代表无法稳定操作的区域。图11 中代表稀相输送向低密度循环流化床的转变边界拟合关系式为

图10 快速床流域内的轴径向颗粒浓度分布（Ug=4 m·s-1）Fig.10 Axial and radial distribution of solids holdups in fast bed regime（Ug=4 m·s-1）

低密度循环床向高密度循环床的转变边界

快速床的操作上限

从图2、图8 到图11，流域转变的边界由一条线、两条线拓展为三条线；稀相输送、快速床、低密度循环流化床、高密度循环流化床等流域被表示于流域图中。

图11 无量纲化的循环床流域图Fig.11 Dimensionless regime diagram of CFB

4 结 论

本文总结分析了气固流化系统，尤其是循环流化床流域划分的研究。并耦合双流体模型与基于稳态EMMS 的曳力模型，模拟了不同操作条件下的流动情况，结论如下。

（1）受操作条件、物性体系与床层几何结构共同影响，气固流化系统在循环床区域的流域划分存在较大分歧，各实验研究对快速床流域的界定不同。

（2）在颗粒浓度较低时，气固流化系统处于较为均匀的流动状态，随着颗粒浓度的增加，颗粒浓度的径向分布逐渐明显，提升管边壁处出现颗粒下滑。将颗粒下滑作为判定快速床出现的依据，可以得到稀相输送向快速床的转变边界，即快速床的操作下限。

（3）在颗粒浓度很大时，不同气速下均出现了节涌现象，代表着快速床操作的终止，即快速床的操作上限。

（4）以上基于流动结构得到的快速床上下限与噎塞型转变线将流域图分为了四部分，模拟在不同区域捕捉到了稀相输送、低密度循环床、高密度循环床等流域的流动特征。

符 号 说 明

Ar——Archimedes数

Cd0——单颗粒曳力系数

D——反应器直径，m

dp——颗粒直径，m

e——碰撞恢复系数

Gs——颗粒循环通量，kg/(m2·s)

g——重力加速度，m/s2

g0——径向分布函数

H——反应器高度，m

Hd——非均匀结构因子

p——压力，Pa

q——脉动能通量，J/(m2·s)

R——反应器半径，m

r——径向位置坐标，m

Re——Reynolds数

Ret——终端沉降Reynolds数

基于实施乡村振兴战略，促进城乡融合发展的需要，宅基地对外流转不仅包括农民、城市居民等自然人之间的流转，还包括农民与企业之间、农民与集体之间、农民与政府之间的流转，对于不同类型的流转，应分别制定相应的流转政策与办法，分类指导、分类规范管理。

U——表观速度，m/s

UFD——快速床的操作下限[13]，m/s

——无量纲表观气速[2]

Umf——最小流化速度，m/s

Uslip——表观滑移速度，m/s

UTF——快速床的操作上限[13]，m/s

Utr——湍动床向循环床的转变气速，m/s

Utp——循环床向输送床的转变气速，m/s

ugy——竖直方向上的真实气体速度，m/s

upy——竖直方向上的真实颗粒速度，m/s

VCA——快速床的操作下限[22]，m/s

VCC——快速床的操作下限[22]，m/s

β——气固相间曳力系数，kg/(m3·s)

γ——碰撞耗散能，J/(m3·s)

ε——体积分数

——颗粒温度，m2/s2

κ——脉动能传导率，kg/(m·s)

λ——体积黏度，Pa·s

μ——动力黏度系数，Pa·s

ρ——密度，kg/m3

τ——黏性应力张量，Pa

Ω——Beranek数

下角标

g——气相

p——颗粒

s——固相