【干货】气液鼓泡塔CFD模拟研究现状

鼓泡塔的性能受许多因素的影响，如：塔径、表观气速、分布器是设计、堆置液面的高度、液体的物性等。在这些因素中，分布器的设计对气含率、气泡直径分布以及混合效果有着重要的影响。因此，从底层机理了解分布器对鼓泡塔性能的影响，对于鼓泡塔的设计与放大至关重要。CFD技术的发展以及电脑计算速度的提高，促进人们采用CFD模拟三维复杂多相体系。正确的采用CFD技术有助于反应器放大研究以及了解鼓泡塔内流体动力学行为。目前，文献中已经证明CFD模拟可以准确预测实验室级甚至工业级鼓泡塔内多相时均流体力学行为。

CFD模拟被用来评估塔径对反应器性能的影响，并根据模拟结果拟合反应器放大的关联式，CFD模拟还用来研究分布器以及鼓泡塔内构件的影响。目前的CFD模拟中绝大多数假定分布器曝气均匀，此外也有部分学者研究分布器提供非均匀曝气的情况。Li等采用CFD研究研究不同分布器对两相流流体力学行为的影响，得出非对称开孔方法以及采用较小孔径的分布器有助于增加鼓泡塔的混合性能。Ranade和Tayalia研究单环或者双环分布器对鼓泡塔内流体力学行为的影响；Michelea和Hempela通过实验以及CFD模拟，研究了烧结板以及单孔分布器对流体力学行为影响；Akhtar等采用CFD模拟研究了不同开孔率对流体力学行为的影响。Bahadori和Rahimi通过CFD模拟研究了不同开孔数目对气含率以及液速的影响。

目前，文献主要采用双流体模型模拟实验室规模甚至工业规模的鼓泡塔，该方法主要是基于局部微元上求解的质量以及动量守恒方程。从理论上讲，该方法是采用平均化的方法描述局部微元的物理现象，由于这些微元之间的相互作用，最终导致鼓泡塔流型的转变。

在采用欧拉双流体甚至多流体模型模拟鼓泡塔时面临着以下三点问题：1采用哪种模型考虑湍流影响；2.采用哪种封闭变量描述气-液之间的相互作用；3.如何考虑局部的气泡直径分布，因为气泡直径影响湍流以及相间作用力。

Jokobsen等认为在双流体框架下RANS（Reynolds-Average Navier-Stokes）方程不适合描述鼓泡塔内湍流现象，这主要是因为k-ε湍流模型是基于各向同性的假设，而鼓泡塔内的湍流属于各项异性。然而，最近的一些研究表明k-ε湍流模型能够获得与LES或者Reynolds Stress模型近似的模拟结果。在Standard、RNG、Realizable形式的k-ε湍流模型中，Standard k-ε湍流模型使用次数最多。

Olmos等认为在鼓泡塔中，部分液相湍流是由于气泡诱导作用产生的，此时湍动能量通过气泡尾迹消耗。因此在模拟湍流时，需要在k和ε输运方程中添加源相来考虑气泡诱导产生的影响。

在绝大多数研究中并没有考虑气相内湍流的影响，这主要因为研究者认为与液相中气泡诱导湍流相比，气相中的气体湍动对鼓泡塔内流型的演化发展影响很小。但Laborde-Boutet等认为对于气-液鼓泡体系而言，当气体作为离散相且液体作为连续相时，该假设合理。但是在湍动鼓泡流型下，局部某些地方的气含率可能会高于0.5，此时气体为连续相，液体为离散相，上述假定未必正确。

气-液之间的相间作用力有很多种，主要包含曳力、升力、湍流扩散力以及虚拟质量力等，在这些力中，曳力对流场的影响占主导地位。曳力主要指气泡在液相中运动所受到的阻力，主要由形体曳力以及粘滞曳力两部分组成。形体曳力是由气泡周围不均匀的压力造成的，形体曳力大小与气泡尺寸以及形状有关。粘滞曳力与流体粘度有关，是由于粘性应力引起的，主要存在于边界层内。当气泡雷诺数增加到一定程度后，流体开始从气泡表面分离，在气泡后面形成涡，从而导致形体曳力的增加。可以认为当气泡雷诺数大于某一临界值时，粘滞曳力可以忽略，此时曳力主要是由于流体开始在气泡后面分离脱落造成的。在纯净或者在受污染的液体中流体的分离脱落程度不一样，流体是否分离也取决于气泡形状以及相界面的的流动。一般情况下，形体曳力以及粘滞曳力均采用曳力系数CD描述，影响曳力系数的因素有很多，如气泡雷诺数、Eotvos数、湍流程度以及气泡表面纯净度等。Kendoush等详细的讨论气泡雷诺数、Weber数、Eotvos数以及弗劳德数对气泡曳力系数的影响。Pan等认为当气泡的形变可以忽略时，曳力系数仅与流型以及液体物性有关。

单气泡的曳力系数与气泡大小、形状以及气泡周围流场有关，Ishii-Zuber根据气泡的形状拟合气泡曳力系数，将气泡分为球形，椭圆形以及球帽形，每种类型气泡具有不同的曳力表达式。Tomiyama根据溶液中表面活性剂的含量拟合不同的曳力模型关联式。

表格 1-3列举出文献中报道的常见的曳力系数经验关联式。绝大多数情况，曳力系数的经验关联式是通过实验测量单独孤立的气泡在静止的液体运动，并结合经验拟合而成，此时单气泡的曳力系数表示为。许多研究者指出，静止流场中测量单气泡曳力系数不适用于泡群体系。在多气泡体系中，气泡受周围气泡的影响，曳力随之发生变化，但是气泡的群体行为对单泡的曳力影响程度尚未等到严格的论证。一般情况下CFD模拟气泡群时，需要对单气泡曳力系数经验关联式进行修正。修正方法有许多种，最常用的修正表达式为：

式中p指的是修正因子，不同工况下选择不一样。一般情况下p值取1。

除了针对单气泡建立曳力关联式之外，有些研究者针对气泡群建立相应的曳力模型。如Krishna通过关联气泡群终端速度获得曳力系数，关联式如表格1所示。通过床层塌落的实验方法测量气泡群的上升速度，并将测量的上升速度代入关联式中获得气泡群对应的曳力系数。

除了基于实验测量以及经验关联获得之外，曳力模型也可以通过理论分析获得。Yang等基于EMMS的思想针对气-液体系建立了DBS模型，通过变分形式的稳定性条件封闭输运方程，通过求解DBS模型获得气-液体系的结构变量，然后将结构变量代入力平衡方程中获得有效曳力系数与气泡直径的比值（CD/db），并将其与CFD模拟进行耦合。与传统曳力模型相比，DBS曳力模型无需调节修正因子便可以准确预测鼓泡塔内局部气含率分布。但为了简化处理，在Yang等工作中CD/db值仅与全局表观有关。严格意义上讲，曳力系数的大小应该与局部流动结构参数有关而不是全局参数。因此，本文的另一个主要目的就是将原来的基于全局表观气速的DBS曳力模型进行改进，将CD/db拓展成与微元流动相关的曳力模型，并与CFD进行耦合。

表格1 常见曳力系数关联式

当气泡在剪切或者旋转的流场中运动，会受到一种侧向力即为升力。根据产生原因分类，升力可以分为两种：Magnus力以及Saffman力。Magnus力主要是由于气泡自身的旋转造成周围流动速度分布不均，从而产生压力差造成的。Saffman力是由于流场本身具有的速度梯度造成的。Zhang和Bogy认为除了气泡的旋转速度比剪切速度大很多的情况，一般情况下Saffman力会比Magnus力大一个数量级以上。在目前的CFD模拟鼓泡塔时，主要考虑由于流场本身具有速度梯度产生的升力，升力的计算表达式如下：

学者们研究了升力系数对CFD模拟结果的影响，Drew和Larey认为在非粘性流体中气泡通常为球形，升力系数取0.5，在粘性流体中，升力系数通常取值较小，如在Lahey工作中升力系数值为0.01，在Dhotre和Smith的工作中CL值为0.1，在Zun的工作中CL取值在0.25~0.30之间， Sundaresan和Sankaranarayanan指出在绝大多数情况下，CL应该取正值。此外也有许多文献报道升力系数为负值，升力系数CL取值大都在-0.1到-3之间。Kulkarni文章中指出，在模拟鼓泡塔时，升力系数CL不能假定为常数，升力系数应该随着径向位置的变化而变化，在鼓泡塔中心处升力系数最小且值为负值，随着径向位置的增加而增加，升力系数逐渐由负数变成正数，在无因次半径大约为0.8时，升力系数才开始为正；在靠近墙壁处，升力系数大小约为0.5。

Tomiyama认为升力系数不仅应该与径向位置有关，而且还与气泡直径有关。Lucas等认为升力能够影响流场的稳定，是影响流型转变重要因素之一。具有负的升力系数的大气泡会破坏流场的稳定，容易导致气含率径向分布不均匀，而具有正的升力系数的小气泡能促进流场的稳定。此外，湍流离散力也有促进流场稳定的功能，例如能一定程度上弥补由负升力系数对流场稳定程度的破坏。尽管学者们对升力进行了大量研究，但升力系数的取值仍然是一个值得争论的问题，且在CFD模拟时，考虑升力后容易造成数值发散。因此在本文的工作中暂不考虑升力的影响。

湍流扩散力主要是由液相速度的湍流脉动造成的，Moraga等认为湍流扩散力影响着垂直管流中气含率径向分布曲线的形状。当采用负升力系数模拟垂直管流时，气体主要集中在管中心处，从而造成模拟结果与实验不符。此时若考虑湍流扩散力，能一定程度扼制气泡向管中心聚集。Lucas等指出由Lahey提出的湍流扩散力模型计算湍流扩散力偏小，无法扼制大气泡在管中心聚集，因此引入了一个与Eotvos有关的湍流扩散力。需要指出，目前文献中报道的湍流扩散力模型基本都是经验的模型，在CFD模拟时，通常与升力一起使用，用来部分抵消添加升力后产生的不合理的结果。在本文的工作中没有考虑升力，因此也忽略了湍流扩散力。

虚拟质量力又被称为附加质量力，指气泡穿过静止的液体时，气泡推动部分液体一起运动，这部分具有一定质量的液体应该获得与气泡相同的速度，就好像虚拟的增加了气泡的质量一样。Lahey等认为虚拟质量力对于气-液两相的加速运动模拟有很大影响。但是对于鼓泡塔而言，气体通过分布器进入鼓泡塔，在很短的时间内便能达到终端速度，因此气泡加速段很小。Oey等通过模拟分析发现曳力在所有相间作用力中占主导地位，决定着CFD模拟能否准确的捕捉鼓泡塔内羽流形态，而虚拟质量力仅起到微调的作用，不会影响整体的分布。

气泡尺寸是影响鼓泡塔性能的重要参数之一，决定着气泡的上升速度以及停留时间，并且影响鼓泡塔内气含率，气-液接触面积甚至传质速率。影响气泡直径分布的因素有很多，如流型、分布器、物性、温度以及压力等。气-液鼓泡塔内流型大致可分为三种：均匀鼓泡流型，过渡流型以及湍动鼓泡流型。均匀鼓泡流型一般在表观气速较低的情况下才会出现，并要求分布器的开孔数量多，孔的排布相对均匀，孔径较小，此时从分布器产生的气泡大小相同，形状几乎为球形。这些气泡安静有序的通过液相到达塔顶，在该流型下气泡几乎不会发生聚并破碎。随着表观气速增大，鼓泡塔由均匀鼓泡流型转化成过渡流型，气泡开始相互碰撞并且发生聚并，从而形成直径相对较大的气泡，直径较大的气泡受升力的作用往塔中心运动，直径较小的气泡开始往边壁运动，此时径向气含率开始成抛物线形态分布。随着表观气速进一步增大，鼓泡塔内流型转变成湍动鼓泡流型，在该流型下气泡聚并破碎作用明显，鼓泡塔内气泡直径分布较广。Ruzicka等指出分布器采用较大的孔径时，不论在哪种气速下鼓泡塔内流型都属于湍动鼓泡流型，此时鼓泡塔内存在明显的气泡直径分布。

对于空气-水体系而言，一般认为初始气泡直径随着压力的增大而减小。Luo等针对氮气-Paratherm NF体系也发现相同的规律。此外，Lin等还发现随着温度的增大气泡直径也会减小。但是有学者对此结论提出异议，Pohorecki等认为温度压力对气泡直径没有影响，Pohorecki等的实验主要针对氮气-水体系，发现无论怎么调节温度以及压力，平均气泡直径基本保持在6.67mm(±12%)。Schafer等认为文献中关于温度、压力对初始气泡直径影响分歧主要是由测量技术、鼓泡塔以及分布器结构不同造成，并通过实验分析得出降低表面张力、增加气体密度或者减小液体粘度，均能导致气泡直径变小，增加压力或者温度能够降低气泡的稳定直径。在水中或者有机溶剂中添加杂质会抑制气泡聚并，从而影响气泡直径，此时鼓泡塔内的气泡直径分布主要由初始气泡直径决定，而初始气泡直径主要受分布器的影响。

气泡直径分布影响气-液之间相间作用力大小，Tabib等指出升力系数应该与局部气泡直径相关，一般认为小气泡的升力系数为正，大气泡的升力系数为负。气泡直径同样影响着曳力大小，然而在传统的双流体CFD模拟中通常是基于单气泡模型框架，在单气泡模型中使用平均气泡直径作为一个集总参数，从而实现与CFD耦合。Xu等基于Schiller-Naumann曳力模型，通过调节平均直径的大小，发现随着气泡直径增大，预测的气含率减小。气泡直径影响着气泡的上升速度以及曳力大小，在过渡流型或者湍动鼓泡流型下气泡直径分布较广，此时若仍采用单一气泡直径有可能导致模拟结果与实验不符。因此Krishna等以及Yang等分别提出双气泡曳力模型。Krishna等将气相分为大、小两类气泡，大、小气泡具有不同的曳力系数，将气-液两相双欧拉模拟转化成大气泡-小气泡-液相三欧拉模拟，导致计算量的增加，并且在三欧拉模拟中忽略大小气泡之间的相互作用。Yang等人基于EMMS思想针对气-液体系提出双气泡理论流体力学模型即DBS模型，通过求解DBS模型获得大、小气泡的流体力学结构参数，并将这些结构参数代入气相力平衡方程中计算出适用于气-液两相的曳力系数，因此DBS曳力模型仍基于双流体模型框架。传统曳力模型仅仅计算曳力系数，在与CFD耦合时需要额外给定气泡直径，然而DBS曳力模型直接给出有效曳力系数与气泡直径的比值（CD/db），无需额外输入气泡直径。

在鼓泡塔内，气泡与气泡发生碰撞后有可能形成更大的气泡，大气泡与湍流涡发生相互作用有可能破碎成两个或者两个以上的小气泡。气泡的聚并破碎影响着鼓泡塔内气泡直径分布以及气-液相界面积，从而影响传质速率。因此，了解气泡聚并破碎机理对于鼓泡塔的设计与放大至关重要。气泡的聚并破碎有可能会引起鼓泡内能流型的变化，如下图所示，对于空气-水体系而言，将小气泡（直径小于5.5mm）从底部注入气-液并流的反应器内，此时反应器底部为均匀鼓泡流。小气泡受到的升力的作用往边壁运动，从而导致边壁处气泡数密度增加，容易聚并成直径较大的气泡（直径大于5.5mm），大气泡受升力的影响往塔中心运动，与其它气泡发生相互作用有可能形成更大直径的气泡。此时反应器顶部气泡直径分布较宽，该部分流型转化成过渡流型甚至柱塞流。

安静鼓泡流型到过渡流型机理分解图

学者们通过群平衡模型（Population Balance Model, PBM）研究气-液两相流中气泡的群体行为。PBM模型有着很久远的发展史，早在18世纪Ludwig Boltzmann提出的Boltzmann方程被认为是首个群平衡方程，在该方程中采用统计学分布的形式描述混乱状态中分子或者颗粒的分布。但是通用的群平衡概念模型开创于19世纪中期，Hulburt和Katz以及Pandolph和Larso分别基于统计学以及连续介质力学概念提出群平衡概念来求解由成核、增长以及凝聚导致的颗粒尺寸的变化。随后研究者基于上述基础进行一系列的改进。尽管如此，在求解实际问题时PBM的普遍适用性仍未得到实验有效的验证。

由于颗粒动理论的复杂性，只有将聚并以及破碎的核函数简化后，群平衡方程才有少量的解析解。为此，研究者们提出了许多数值解法来求解群平衡方程，最常用的方法有蒙特卡罗法（Monte Carlo Methods）、矩量法（Moments Methods）以及类方法（Class Methods）。蒙特卡罗法采用总体统计的方法求解群平衡方程，与其它方法相比，在追踪多维体系中颗粒变化方面更加灵活和准确，但是蒙特卡罗法的准确度与计算的颗粒数目成正比，因此计算量大，致使很难与CFD耦合。目前与CFD耦合主要采用的是矩方法以及类方法。

虽然目前已有大量的文献进行群平衡模拟研究，但总体上讲，目前仍存在以下几点问题限制群平衡的应用以及推广方面：

1.相间作用力模型的限制。目前绝大多数相间作用力模型主要基于单颗粒建立或者校正的，因此很难适用于颗粒紧密堆积的情况，因为在颗粒运动时会受到周围颗粒的影响；2.聚并破碎核函数的限制。聚并破碎核函数主要基于气泡为球形的假设，该假设决定了核函数只能适用于均匀鼓泡流型。在非均匀鼓泡流型下会存在直径较大的气泡，大气泡会扭曲变形并在后面产生尾迹，尾迹会对气泡聚并会产生明显影响，然而现有模型无法考虑这点；3.计算量的限制。为了降低计算成本提高计算精度，学者们也提出了各种算法将PBM与CFD耦合，如正交矩量法（Quadrature Method of Moments, QMOM）、直接正交矩量法等。尽管随着科技的发展，目前计算能力有了显著的提高，但是将PBM与CFD耦合时计算成本仍然巨大。

综上所述，基于欧拉框架下想要准确模拟鼓泡塔时，必须选用合适的湍流模型、相间作用力模型、以及采用适当的方法考虑局部气泡直径分布。了解这些模型特性以及所使用的范围，根据具体体系选择合适的模型。