鼓泡塔反应器中两相流动CFD-PBM 耦合数值模拟

贾翔飞， 钱嘉澍， 吴幼青， 陈剑佩

（华东理工大学 1. 化学工程联合国家重点实验室；2. 煤气化及能源化工教育部重点实验室，上海 200237）

我国的自然资源呈现贫油、少气、富煤的现象，每年需要进口大量的石油。为了保障能源安全，我国对煤液化技术进行了技术储备并且成功将该技术工业化。煤液化技术的关键设备是液化反应器，其本质上属于鼓泡塔反应器。该反应器具有良好的混合性能并且生产能力大、结构简单、易于操作[1-2]。

煤液化反应器中涉及到氢气（气）、供氢溶剂（液）、煤粉颗粒（固）的流动问题，但反应器中的煤粉颗粒较小并且浓度较低，煤粉颗粒实际上是伴随着液相一起运动[3]，因此通常将煤粉颗粒与液相作为混合均匀的浆相处理，将气-液-固三相流动简化为气-浆两相流动。根据反应器操作的表观气速不同，反应器内的状态可以分为均匀鼓泡状态、过渡状态、湍动鼓泡状态。在均匀鼓泡状态下，鼓泡塔内的气泡尺寸小而均匀，气含率与表观气速呈现近似线性关系[4]。工业生产过程中的鼓泡塔反应器一般处于湍动鼓泡状态[5]，在该状态下，流动结构由于气泡大量的聚并和破碎行为而变得非常复杂。强烈的聚并和破碎行为使得鼓泡塔内的气泡尺寸分布范围较宽，并对鼓泡塔内与反应器设计和放大有关的气含率、液速、热量和质量传递速率等参数产生影响[1]。因此，反应器内的气泡破碎聚并以及基础流体力学行为对反应器的设计放大以及安全运行有着重要影响。

在对反应器进行设计和放大时，计算流体动力学（CFD）逐渐成为一种重要的工具[2,6-8]。在综合考虑计算量和计算精度的前提下，通常采用Eulerian-Eulerian 方法[9]对鼓泡塔进行模拟。使用该方法研究气液鼓泡塔的流体力学行为关键在于合理地描述气-液相间作用力和气泡的大小分布。相间作用力包括曳力、升力、壁面润滑力及湍流扩散力。曳力是由气液之间的摩擦产生的，曳力的大小决定了气泡在液相中的上升速度并最终影响鼓泡塔内的气含率的大小，而升力、壁面润滑力和湍流扩散力会影响气含率在径向上的分布。

由于气泡与气泡之间存在复杂的相互作用，鼓泡塔内的气泡存在局部气泡直径分布。早期的研究者通常将该气泡直径分布简化为一个恒定的当量气泡直径[10-12]。在均匀鼓泡状态下，气泡之间的相互作用较弱，气泡的尺寸分布较窄，气泡直径均一，通过调节曳力模型及气泡直径会获得与实验值较吻合的结果。但在湍动鼓泡状态下，气泡之间强烈的相互作用造成气泡较宽的尺寸分布，因而上述对气泡直径的简化无法获得与实验值吻合的结果。

由于气泡之间的相互作用对鼓泡塔内的流体力学参数有重要影响，因而群体平衡模型（PBM）逐渐被应用于描述鼓泡塔内的气泡尺寸的变化及分布，并与CFD耦合以对鼓泡塔进行模拟[13-16]。多年来，人们进行大量的工作来建立合理的聚并及破碎模型[16-19]，并希望将其应用于CFD-PBM耦合模型中。

本文将一个新颖的气泡破碎模型[19]植入PBM并与欧拉双流体模型进行耦合。首先在ANSYS Fluent 平台中将CFD模型所需要的相间作用力模型通过用户自定义函数（UDF）加入，然后将PBM 模型中聚并、破碎模型源项通过UDF植入并与CFD模型耦合建立CFD-PBM 模型。采用该模型对实验室尺寸的圆柱鼓泡塔进行二维轴对称模拟，并与文献中的实验数据进行比较。

1 模型方法

1.1 双流体模型的控制方程及湍流模型

文献[9]中提到，采用重整化群k-ε（RNGk-ε）湍流模型可以更好地描述鼓泡塔中的流动。本文采用RNGk-ε湍流模型计算液相中的湍动能和湍动耗散率。

1.2 相间作用力模型

1.2.1 曳力 曳力的大小会影响气泡在液相中的上升速度，并进一步影响鼓泡塔内气含率的大小，曳力通常采用式（3）计算[14]：

1.2.2 升力 液相中速度梯度的存在会使气泡受到一个横向升力，这会对流体力学参数的径向分布产生影响，升力系数的计算采用Tomiyama 模型[20]。

1.2.3 壁面润滑力 气泡受到的壁面润滑力采用Tomiyama 模型[20]。

1.2.4 湍动耗散力 湍动耗散力采用Lopez de Bertodano模型[21]。

1.3 群体平衡模型

1.3.1 群体平衡方程（PBE） 在鼓泡塔中，存在着数量非常庞大的气泡，这些气泡会发生聚并和破碎，最终使鼓泡塔中的气泡存在一定的尺寸分布。PBE（式（15））以及合理的聚并、破碎源项的引入使得鼓泡塔中的气泡尺寸分布得以描述。

其中：Si为第i组气泡中PBM 的源项[22]，在目前的研究工作中，成核及生长速率被忽略，只考虑聚并生成、聚并消失、破碎生成及破碎消失项（分别表示为BC,i，DC,i，BB,i，DB,i）。各源项的计算采用文献[22]的方法。

1.3.2 聚并核函数 鼓泡塔中引起气泡之间相互碰撞的机理通常有3种，分别是液相湍流引起的碰撞、气泡上升速度差异引起的碰撞以及大气泡尾涡引起的碰撞。但不是每次气泡碰撞都会使两个气泡发生聚并，所以通过引入气泡聚并效率的概念来表示两个气泡聚并的可能性。最终，可以通过式（17）来计算直径为di和dj的气泡的聚并速率。

鼓泡塔中总的聚并速率可以认为是3种机理聚并（湍流涡引起的聚并、气泡尾涡引起的聚并、气泡上升速度不同引起的聚并）速率的线性叠加。

（1）湍流涡引起的聚并

①碰撞频率

②聚并效率

根据文献[25]提出的液膜理论，气泡的相互作用时间和气泡聚并所需时间的相对大小决定着气泡聚并的概率。文献[25-26]给出了直径为di、dj的气泡的聚并效率。

（2）由气泡尾涡引起的聚并

①碰撞频率

鼓泡塔中气泡在尺寸较大时会由球形变成球帽形，并在气泡的尾部形成涡流。此时若其他气泡进入该大气泡的尾涡影响区域，就会被该气泡的尾涡夹带而加速上升。该过程会使大气泡和夹带气泡发生碰撞并聚并，文献[16]提出尾涡引起的气泡碰撞频率为：

由于小气泡的尾涡相对较小，因此这里只考虑大气泡尾涡引起的聚并，所以在碰撞频率中加入因子Θ。

②聚并效率

尾涡引起的聚并效率采用文献[27]提出的关联式。

（3）由于气泡上升速度不同引起的聚并

①碰撞频率

由气泡上升速度不同引起的聚并模型与液相湍流引起的聚并模型相似[23]，其气泡碰撞频率表达式为：

②聚并效率

由于气泡上升速度不同造成的聚并，文献[14, 24]将聚并效率设置为0.5，即

1.3.3 破碎核函数 对于气泡的破碎过程，本文只考虑湍流涡引起的破碎，并且破碎过程为二元破碎，即每个母气泡破碎成两个子气泡，故气泡破碎速率可采用式（26）进行计算[16,18]。

（1）气泡破碎碰撞频率

目前，一些研究者[17,28-30]认为大尺寸的湍流涡会把部分能量传递给气泡从而使气泡发生破碎。因此将λ≤d的涡和λ>d的涡与气泡的碰撞频率分别进行计算[19]。

气泡破碎时必须同时满足能量约束和压力约束[13,15-16,31]。当λ>d时[28]，湍流涡中的能量只能部分作用于气泡并导致气泡的破碎，因此在量纲为一湍流涡能量中引入修正因子Ψ。

对气泡颈处的受力进行分析最终得到压力破碎因子为：

能量约束条件的推导可见文献[18]，能量约束因子为：

1.4 模型方程求解

基于文献[33]设置相关实验并进行二维轴对称模拟,采用空气-水作为实验介质，鼓泡塔的相关几何模型见图1（a）。鼓泡塔网格采用ANSYSICEM CFD划分并对进口处以及壁面处进行加密处理。最终经过网格无关化检验后确定四边形网格数量为9 581，计算过程中的时间步长取0.002 s。气体均匀地从底部r/R<0.9的中心区域网格进入（图1（b））。实验的气速范围为0.02～0.10 m/s，覆盖了均匀鼓泡流和湍动鼓泡流。所有计算均在ANSYSFluent 平台上进行，压力-速度方程使用SIMPLE算法处理，动量方程采用QUICK 格式，体积分数方程等均采用一阶迎风格式离散，松弛因子均为默认值。模拟中所用到的相间作用力模型均采用UDF加入Fluent 求解平台。

图1 鼓泡塔反应器几何模型（a）及网格截面（b）Fig.1 Geometry (a)and mesh (b)of bubble column reactor

鼓泡塔反应器的入口及出口边界条件为速度进口和压力出口，壁面边界条件为无滑移边界。PBE方程的求解采用离散法，同时破碎和聚并模型均采用UDF实现。本文气泡尺寸为0 .5≤db≤40.3 mm，分为20组，各组气泡大小的计算如式（36）所示。

2 结果与讨论

2.1 概述

采用CFD-PBM 耦合模型对鼓泡塔反应器内的气含率及轴向液速分布情况进行模拟，气泡聚并及破碎模型通过UDF接口植入ANSYSFluent 中。其中鼓泡塔反应器入口表观气速分别为0.02 m/s和0.10 m/s，并与相应表观气速下的实验数据[33]进行对比。

2.2 破碎模型验证

其中：L和T分别为长度尺度和时间尺度[34]，可分别通过式（39）、（40）进行计算。

图2 量纲为一破碎速率预测值与实验值对比Fig.2 Predicted and experimental values of dimensionless breakup rate

图2中实线代表模型预测的量纲为一破碎速率，数据点为文献[35-37]中的实验测量值。从图2中可以看出，在量纲为一气泡直径较大时模型预测值略小于实验值，在量纲为一气泡直径较小时模型预测值与实验值基本吻合。

图3比较了子气泡分布的模型预测值与实验测定值，其中实验数据来自文献[32]。实线为模型预测的子气泡分布概率密度曲线，灰色直方图为概率密度曲线在不同区间上的积分值，红色直方图为相应区间上的实验测量值。从图中可以看出由于模型采用了压力约束条件，在气泡破碎体积分数fv趋于0时子气泡分布的概率密度β（fv,d）曲线也趋于0。同时，从图3中也可以看出，模型预测的结果为气泡更倾向于不均匀破碎，这也与实验测定的结果相吻合[32]。

图3 子气泡分布预测值与实验值Fig.3 Predicted and experimental valuesof daughter bubble size distribution

2.3 升力的影响

升力是由于液相速度梯度差使得气泡在径向方向上受到的力，其决定了气泡在径向的运动方向。升力的加入会使大气泡向反应器中心移动，因此会对气含率和轴向液速在径向上的分布产生影响。图4和图5分别示出了添加升力对气含率和轴向液速在径向分布的影响，图4、5中实线分别为未添加升力时的气含率、轴向液速径向分布。从图4可以看出，当不加入升力时，气含率在径向上的分布更加平坦，在鼓泡塔中心处模拟值小于实验测量值，在靠近壁面处模拟值远大于实验测量值。由图5可以看出，升力的加入使得径向速度梯度变大，模拟值与实验值更加吻合。总体来说，升力对气泡在鼓泡塔内流体力学参数的径向分布十分重要，升力的加入使径向气含率的分布与实验值更加吻合。

2.4 气含率

图6示出了不同表现气速下鼓泡塔反应器内气含率及轴向液速分布云图。由图6中气含率分布云图和图4中气含率径向分布曲线可以看出，鼓泡塔内气含率沿径向方向逐渐降低。由于鼓泡塔内存在液相的循环，液体在中心处上升而在壁面处下降，这种液相的循环使得气含率在鼓泡塔的中心较高，而在壁面处较低。从气含率云图与径向分布图可以看出，在低表观气速（0.02 m/s）下，气含率的径向分布比较均匀，但在高表观气速（0.10 m/s）下，鼓泡塔中心与壁面处气含率的径向分布差异较大。

由图4可以看出，在高表观气速下，鼓泡塔壁面处气含率的模拟值大于实验值；在低表观气速（0.02 m/s）下，气含率模拟值在径向略小于实验值。从整体上看，在低表观气速下，目前模型的计算值低估了气含率。

图4 不同表观气速下气含率径向分布Fig.4 Radial profiles of gas holdup at different superficial gas velocities

图5 不同表观气速下轴向液速径向分布Fig.5 Radial profiles of axial liquid velocity at different superficial gasvelocities

2.5 轴向气、液速

由图6中轴向液速分布云图及图5中轴向液速径向分布图可以看出，鼓泡塔内的轴向液速沿径向方向逐渐降低。由图6可以看出，在径向方向上，中心处的轴向液速最大，沿径向方向，轴向液速逐渐降低并变为负值。此时的负值表示液体的流动方向变为与原来流动方向相反。造成这种现象的原因是鼓泡塔内的气相在浮力的驱动下沿鼓泡塔上升，进而造成鼓泡塔内的液相产生这种轴向循环。

图5中将现有模型的模拟值与实验值进行了比较。在高表观气速下，即处于湍动鼓泡状态时[9]，模拟值与实验值吻合良好，但在中心处略微低估了轴向液速。在低表观气速下，即均匀鼓泡状态时，模拟值与实验值大致吻合，但在鼓泡塔中心处，模拟值要高于实验值。

图6 不同表观气速下鼓泡塔反应器内气含率及轴向液速分布云图Fig.6 Gas holdup and axial liquid velocity contours at different superficial gasvelocitiesin bubblecolumn reactor

图7 不同表观气速下气泡轴向速度径向分布Fig.7 Radial profiles of bubble axial velocity at different superficial gasvelocities

图7所示为鼓泡塔内不同表观气速下气泡轴向速度的径向分布，可以看出，随表观气速增加，气泡的轴向速度增加，并且气泡轴向速度径向分布的不均匀性也在增加。结合图5和图7也可以看出，气液之间的相对速度变大。这是因为随着表观气速的增加，反应器状态由均匀鼓泡状态变为湍动鼓泡状态[4]，气泡大小分布变得更宽，反应器内大气泡数量增加[16]，大气泡尾涡产生的加速效应使得气液之间的相对速度变大。

2.6 湍动耗散率

气泡的聚并与破碎速率均与湍动耗散率有关，并且气泡的聚并与破碎速率会对局部气泡直径产生影响，进而影响气液相间作用力的大小，因此合理地描述鼓泡塔内的液相湍动耗散率对于流场的准确模拟非常重要。

图8所示为不同表观气速下液相的湍动耗散率的径向分布。可以看出，湍动耗散率与表观气速呈正比关系。同时，由于壁面的存在以及流体黏性的作用使得壁面处液相的速度梯度较大，因此壁面处的湍动耗散率较大。

图8 不同表观气速下湍动耗散率径向分布Fig.8 Radial profiles of turbulent dissipation rate at different superficial gasvelocities

2.7 气泡尺寸分布

图9所示为鼓泡塔H/D=3截面处不同表观气速下气泡尺寸分布，该图横坐标为气泡直径，纵坐标为基于气泡体积处理得到的气泡尺寸分布概率密度函数。从图中可以看出，随着表观气速增加，H/D=3截面处平均气泡尺寸分布曲线变得更加低矮，同时大气泡的体积分数增加。

图9 H/D=3处不同表观气速下气泡尺寸分布Fig.9 Bubble size distribution at H/D=3 at different superficial gasvelocities

图10所示为H/D=3处，r/R=0和r/R=0.5处不同表观气速下的气泡尺寸分布。从图10（a）可以明显看出，随着表观气速增加，鼓泡塔由均匀鼓泡状态转变为湍动鼓泡状态，鼓泡塔中心位置的气泡尺寸分布变得更宽，大气泡的体积分数增加。这主要是由于湍动鼓泡状态下，较大的表观气速造成鼓泡塔内的气泡数目增加。单位体积内气泡数量的增加使得气泡之间的相互作用更加频繁，在液相湍流的诱导下气泡聚并产生更多的大气泡，而大气泡数量的增加也使得更多的气泡因为气泡尾涡效应和气泡上升速度差效应而发生聚并。另一方面，鼓泡塔内聚并产生的大气泡拥有更大的碰撞截面积，因此会有更多的湍流涡体与其发生相互作用。与此同时，高表观气速下湍动耗散率增加使得湍流涡体的湍动能增加，最终表现为气泡的破碎速率增加。这些因素的综合作用使得气泡尺寸的分布范围变宽。同时，在升力的作用下，大气泡趋向于向鼓泡塔中心运动，故可以看到高表观气速下，r/R=0处气泡尺寸分布较宽且大气泡数量较多。从图10（b）可以看出，沿径向方向远离鼓泡塔中心处，高表观气速与低表观气速的气泡尺寸分布差异在减小。但高表观气速下气泡尺寸分布范围依然大于低表观气速下的范围。

图10 不同表观气速下r/R=0（a）和r/R=0.5（b）的气泡尺寸分布Fig.10 Bubblesizedistribution at r/R=0(a),and r/R=0.5(b)at different superficial gasvelocities

图11所示为表观气速为0.10 m/s时不同径向位置处的气泡尺寸分布，从图中可以非常清晰地看出，沿径向靠近鼓泡塔中心，大气泡的体积分数在增加，小气泡的体积分数在减少。

图11 径向位置处气泡尺寸分布Fig.11 Bubblesizedistribution at different radial positions

3 结 论

将一个新型的破碎模型通过UDF植入ANSYS Fluent 软件中，并将聚并和破碎模型与双流体模型进行双向耦合建立一个通用的CFD-PBM 耦合模型，在考虑了不同机理对气泡破碎和聚并的影响以及气液相间作用力后，对处于均匀鼓泡状态与湍动鼓泡状态下的鼓泡塔流体力学参数进行研究。通过数值模拟得出以下结论：

（1）在对鼓泡塔进行模拟时，升力是否加入对于气含率和轴向液速等流体力学参数的轴向分布有重要影响，升力的加入使模拟数据在径向上的分布与实验值吻合得更好。

（2）在对破碎模型的气泡破碎速率和子气泡分布预测值进行验证之后，将模型在不同表观气速下的气含率、轴向液速模拟值与实验值进行了对比。结果显示CFD-PBM 耦合模型的计算值与实验值吻合较好，这也证明了该模型可以用于湍动鼓泡状态下鼓泡塔的流体力学模拟。

（3）鼓泡塔内液相湍动耗散率随表观气速的增加而增加。同时，壁面的存在以及流体黏性的作用使得壁面处液相的速度梯度较大，因此壁面处的湍动耗散率较大。

（4）随表观气速增加，鼓泡塔内的气泡数目增加，单位体积内气泡数量的增加使得气泡之间的相互作用更加频繁，气泡聚并速率增加。另一方面，鼓泡塔内聚并产生的大气泡拥有更大的碰撞截面积，因此会有更多湍流涡体与其发生相互作用。与此同时，高表观气速下湍动耗散率增加使得湍流涡体的湍动能增加，最终表现为气泡的破碎速率增加。这些因素的综合作用使得气泡尺寸的分布范围变宽。