基于气泡群相间作用力模型的加压鼓泡塔流体力学模拟

刘鑫，张煜，张丽，靳海波

（1北京石油化工学院化学工程学院，北京 102617；2中科合成油技术有限公司研发中心，北京 101407）

基于气泡群相间作用力模型的加压鼓泡塔流体力学模拟

刘鑫1,2，张煜2，张丽2，靳海波1

（1北京石油化工学院化学工程学院，北京 102617；2中科合成油技术有限公司研发中心，北京 101407）

目前，多数文献报道了冷态加压湍动鼓泡塔内流动特征，并且通过实验数据回归相关经验关联式。然而，此类关联式适用范围有限，难以直接外推到工业鼓泡塔反应器条件。因此，在FLUENT平台上建立了基于气泡群相间作用力的、动态二维加压鼓泡塔计算流体力学模型。通过数值模拟考察了操作压力为0.5～2.0 MPa，表观气速为0.20～0.31 m·s−1，内径0.3 m鼓泡塔内流场特性参数分布，并且与冷态实验数据进行比较。结果表明，采用修正后的气泡群曳力模型、径向力平衡模型以及壁面润滑力模型描述气泡群相间作用力，能够较为准确地反映平均气含率和气含率径向分布随操作压力和表观气速变化的规律。

鼓泡塔；气液两相流；计算流体力学；气泡群曳力模型；径向力平衡模型；壁面润滑力模型

引 言

煤间接液化技术作为煤炭清洁高效利用的重要手段之一，它是以煤炭为原料制取目标产品为液态燃料的技术。随着以费托合成工艺为核心的煤炭间接液化技术的重要性日益突显，作为费托合成反应的载体，费托合成反应器得到了快速发展。费托合成反应器先后经历了固定床、流化床和浆态床 3个不同阶段，其中采用浆态床反应器生产的柴油，具有高能量密度、低排放、品质好等优点，因此，浆态床已成为煤制油进行费托合成反应的主流反应器[1]。

工业鼓泡塔反应器通常带压操作，并且表观气速较高，处于湍流鼓泡状态。目前对于加压鼓泡塔的研究大多集中于冷态实验测量和流动规律的总结，根据实验数据回归经验关联式，缺乏能够详细描述加压鼓泡塔的流体力学模型。Krishna等[2]根据加压气泡界面的Kelvin-Helmholtz不稳定性分析和气泡曳力模型，将压力影响作为气相密度修正项ρG/ρG,1atm引入曳力系数计算公式。建立了双气泡二维轴对称CFD模型，用CFX软件模拟了塔径0.15 m加压湍动鼓泡塔充分发展段的液相轴向扩散系数及气含率，并与冷态实验数据进行比较。在高气速下，该模型计算的气含率明显高于实验测量值。Chen等[3]建立了结合气泡群体平衡方程（PBM）的 3D双流体和混合物模型，用 FLUENT模拟了塔径为0.162 m加压鼓泡塔流场。考虑到Krishna等[2]的模型在高气速下对气含率高估，认为造成偏差的原因可能是，高气速下两相湍流和尾涡的卷吸作用影响显著。为了改进模型的计算准确度，Chen等[3]引入(ρG/ρG,1atm)0.25修正气泡曳力项。该模型采用耦合PBM的3D模拟，为保证计算准确度，离散网格数量随鼓泡塔尺寸提高而快速增加，需要大量计算资源。以现有计算机运算能力，仍然难以采用此模型进行大型加压鼓泡塔反应器 3D模拟和放大设计。

Joshi[4]认为气泡所受的横向作用力决定其沿塔径向的运动模式，是湍动鼓泡塔形成气含率不均匀分布的重要影响因素；Tabid等[5]通过3D动态模拟鼓泡塔流场，比较了不同相间作用力的影响，认为曳力、升力和湍流扩散力对气含率和轴向液速分布影响显著。文献[6-8]通过引入气泡径向力平衡机制，确定气含率径向分布，一维及二维数值模拟结果与实验数据符合较好。综上所述，二维轴对称模型耦合适当的相间作用力模型，能够在可接受的计算代价下，得到较为准确的常压鼓泡塔流场模拟结果。本文拟建立动态二维轴对称加压湍动鼓泡塔模型，包括径向力平衡机制、气泡群曳力模型和壁面润滑力，根据内径0.3 m加压冷态鼓泡塔冷态实验数据，拟合气泡群曳力系数参数，使模型具备预测加压湍动鼓泡塔反应器内多相流场特性，为工业浆态床反应器放大和内构件优化提供指导。

1 数学模型

1.1 守恒方程

本文所提出的双流体加压鼓泡塔模型，采用雷诺时均的控制方程[9-11]。与操作压力相比，加压鼓泡塔静液压相对较小，可近似认为气体密度基本不随轴向高度变化，冷态条件下可忽略两相间的传热和传质作用，将相间作用力以体积源项的方式，添加至动量守恒方程组中，控制方程可表示如下：

连续性方程为

文献[9, 11]给出了气液两相表达形式一致的动量守恒方程

文献[12-17]报道采用标准k-ε模型能够描述鼓泡塔内液相湍流场，本文也采用该模型计算液相湍流动能k和湍流动能耗散率ε的分布，以及湍流黏度μt，模型方程为

式中，参数Cμ=0.09，C1ε=1.44，C2ε=1.92，σk=1.0，σε=1.3。

1.2 相间作用力模型

鼓泡塔双流体模型中，主要通过体积平均的相间作用力源项描述气液相间动量传递[4-5,10-11]。本文所提出的加压鼓泡塔流体力学模型，主要考虑气泡曳力、升力、湍流扩散力和壁面润滑力。另外，虚拟质量力表现为颗粒在加速运动时，使周围流体也加速运动所需要额外的力，Khopkar等[18]认为虚拟质量力在气体分布器附近更重要，而在流动主体区域的作用远小于曳力，且本文通过实际工况模拟发现，添加虚拟质量力模型对模拟结果几乎无影响，因此本文忽略了虚拟质量力。

1.2.1 曳力 鼓泡塔内气泡曳力是相间作用力研究的重点之一，也是现有各种流体力学模型主要考虑的相间作用之一，甚至是许多模型中唯一考虑的相间作用力。对于鼓泡塔双流体模型，动量守恒方程中的曳力源项通常表述为单位体积内气泡群所产生的合曳力，表示为

在加压湍动鼓泡流状态下，气含率相对较高，气泡在浮升过程中会受到附近气泡的影响，需要引入气泡群曳力系数进行修正，文献中通常将气泡群曳力系数CD与单个气泡曳力系数CD0与气含率相关联进行修正。

文献中给出的单气泡曳力系数模型通常与气泡Reynolds数Re相关。Schiller & Naumann模型常用于气液两相流模拟，在鼓泡塔流场数值模拟过程中稳定性较好，模型方程如下

与 Schiller & Naumann模型相比，Morsi & Alexander模型更为完整，其曳力系数对气泡Reynolds数的划分更详细，曳力系数与Reynolds数的关系可表示为

式中，a1、a2、a3为模型参数，随Reynolds数取值区间的改变而对应不同的数值。但是，该模型的计算稳定性不及其他模型。

Lau等[19]和 Roghair等[20]采用直接数值模拟（front-tracking model，前缘追踪模型）单个气泡及气泡群浮升过程，通过模拟结果拟合得到气泡群曳力系数与气含率和Eo数的关联式

但是，Roghair等[20]未考虑操作压力对气泡群曳力系数的影响，因此，本文在式（9）的基础上，根据冷态实验[21]测量的全塔平均气含率（压差法测量）和气泡直径（电导探针测量），调整单气泡曳力模型中的曳力系数，使平均气含率的模拟结果符合实验测量值。根据调定后的曳力系数、气含率、气泡直径及新添加的量纲1压力项，拟合模型参数，得到适用于加压条件下改进的气泡群曳力系数计算公式

其中，CD,0表示单个气泡的曳力系数，本文采用Schiller & Naumann模型计算CD,0，湍动鼓泡塔中气泡的相对Reynolds数Re通常均大于1000，CD,0可取常数0.44。

1.2.2 升力 鼓泡塔内气泡浮升运动过程中，由于液相在其运动方向的两侧流动不对称，会造成气泡两侧压力不平衡，从而产生垂直于气泡运动方向的升力。大气泡在鼓泡塔内受到液相剪切流动、曳力以及涡旋的作用后，易发生变形，并且在气泡的后方产生偏斜尾涡，推动气泡横向运动。Drew等[22]给出了单位体积内，气泡群所受合升力的计算公式

在鼓泡塔充分发展段，如果气泡升力系数 CL所取的符号不同，气泡升力会指向塔中心或者塔壁，使气含率径向分布对应形成中心高的抛物线分布，或者形成壁面峰的鞍形分布，此外，升力系数 CL取值大小会显著影响气含率径向分布的梯度。

1.2.3 湍流扩散力 对于湍动鼓泡塔流场，湍流扩散力体现了气液两相随机运动所产生的扩散作用，在该力的作用下气含率径向分布趋于均匀。Lopez de Bertodano[23]将气泡的随机运动类比于分子热运动，推导出单位体积气泡受到液相涡旋所受的扩散力表达式

考虑到湍流扩散力主要在气含率较高的区域起作用，因此，引入极限函数 fTD,limiting对湍流扩散力作用区域进行调整，湍流扩散力可表示为

其中，εG,1和εG,2为湍流扩散力起作用的局部气含率下限及上限值，本文所述模型中将其分别取值为0.3和0.7。

文献[7-8]中认为鼓泡塔内径向力作用力对气含率的分布起着主导作用，当升力和湍流扩散力达到动态平衡，气含率的分布也随之趋于稳定。张煜等[7]通过常压鼓泡塔冷态实验数据和一维模型计算结果，确定了气泡升力系数 CL和湍流扩散力系数CTD符合如下函数关系

式（15）中的系数0.2是一个固定的值，不随实验条件而改变。本模型假设该径向力平衡机制同样适用于加压湍动鼓泡塔流场模拟，湍流扩散力系数 CTD取 1.0，液含率与升力系数对应情况如表 1所示。

表1 不同液含率下的升力系数Table 1 Lateral lift coefficient under different liquid holdup

1.2.4 壁面润滑力 当浮升气泡接近壁面时，气泡表面液体速度分布受到壁效应的影响，表现为气泡与壁面之间的液相排出速率下降，而在另一侧的排液速率增加，从而产生将气泡推离壁面的力，称之为壁面润滑力[24]。本文模型采用 Tomiyama[25]的壁面润滑力模型，单位体积的气泡所受的壁面润滑力可表述为

式中，壁面润滑力系数CWL可表示为

式中，yw为壁面距离，Cw可表示为气泡Eo数的分段函数

在进行加压鼓泡塔轴对称二维动态模拟时发现，当网格密度较高时，在充分发展段经常出现壁面气含率高中心低的情况，这与实验测量结果不符。当模型中考虑壁面润滑力后，能够有效避免出现气含率分布的非物理解。

2 模型设置及求解方法

本文参照文献[21]所述内径0.3 m冷模加压鼓泡塔装置生成二维计算网格，计算区域宽度为0.15 m、高度为6.6 m。根据表观气速和操作压力设置初始静液位高度H0，防止膨胀床层高度超过计算域高度，造成液体溢出，静液位设定高度为2.5～3.0 m。设定H0以下的液含率为1.0，H0以上的液含率为 0，气液相均处于静止状态。边界条件为：鼓泡塔底部的孔板分布器（r/R＜0.8）简化为气相速度进口边界条件；塔顶气相出口设置为压力出口；中心轴采用轴对称边界条件；壁面的气液相均设定为无滑移条件，湍流方程设定为标准壁面函数，如图 1所示。本文采用文献[21]的实验数据作为检验模型的依据，其中局部气含率、气泡上升速度及气泡Suater直径，是采用双电导探针，根据气液两相电导率差异较大的原理来测量的，对气泡上升速度的测量是通过考察电导针处，一段时间时间内统计所有气泡上升速度的数量平均获得的。

通过用户自定义方程UDF，将前文所述的各气泡界面力模型嵌入FLUENT计算程序，全塔设定单一气泡直径等于秦玉建[21]的冷态实验测量值。模型采用 FLUENT 15.0进行动态模拟，采用Phase-Coupled SIMPLE算法求解压力-速度耦合，离散空间梯度使用Least Squares Cell based方法，导数采用一阶迎风格式，时间导数采用一阶隐式格式，为防止发散，初始时间步长为0.001 s，待流场充分发展后，可以逐渐增加步长至0.004 s。模拟结果表明，当计算物理时间达到20 s后，加压鼓泡塔充分发展段的径向气含率分布基本不变，可近似认为流场充分发展，以表观气速0.31 m·s−1，操作压力为0.5 MPa，距气体分布器2.4 m处气含率径向分布为例，如图2所示。本文取模拟时间50 s后的流场进行时均，图3给出了不同采样时间对应的时均气含率分布，可以看出当流场基本稳定后，采样时间大于25 s后的气含率分布基本一致。因此，文中所述模拟均提取流场发展50 s后的数据，再时均50 s作为数值模拟结果。

图4给出了网格数量对充分发展段气含率径向分布模拟结果的影响，考察了网格数量为3960、5940、17160和47520个4种不同的网格，并且均在近壁区域加密了4层边界层网格。从图中可以看出，加密网格对中心气含率分布影响较小，主要使近壁面气含率分布变得更为平滑。综合考虑计算代价和数值模拟精度，采用网格数为5940个。

图1 加压鼓泡塔结构及边界条件Fig.1 Schematic diagram of pressurized bubble column and boundary conditions

图2 模拟的气含率径向分布随时间的变化Fig. 2 Simulated radial gas holdup profiles varied with time

图3 模拟时均气含率分布随采样时间的变化Fig. 3 Simulated radial gas holdup profiles varied with sampling time

图4 网格密度对气含率径向分布的影响Fig. 4 Influence of mesh density on radial gas holdup distribution (VG=0.27 m·s−1，P=0.5 MPa，H=2.6 m)

3 结果与讨论

3.1 流场分布

鼓泡塔流场的最主要特征是气含率在径向上的不均匀分布，即中心区域气含率高，边壁处气含率较低。这种气含率的不均匀分布形成密度差，从而驱动液相流动，形成液相塔内的大尺寸循环流动，呈现出中心区液相向上运动，近壁区液体向下回流的规律。图5是CFD模拟得到的，表观气速为0.27 m·s−1，操作压力为1.5 MPa，内径0.3 m加压鼓泡塔内时均气含率二维分布。模拟结果表明，在鼓泡塔底部，回流液体在壁面附近富集，使气泡向塔中心集中。经过约1倍塔径发展后，气含率径向分布呈现出中心高、边壁低的趋势，并且不随高度变化，形成充分充分发展段，直至接近塔上部液位处。在气体分布器上方2.6 m和3.0 m处，冷态测量得到的气含率径向分布基本相同，可认为处于充分发展段，模拟的气含率分布与实验结果符合较好，这表明，本文所提出的基于气泡群界面力的加压鼓泡塔二维模型，能够确定气含率分布。

图5 气含率二维分布模拟结果与实验数据的比较Fig.5 Comparison of simulated 2-D gas holdup distribution and experimental data(VG=0.27 m·s−1，P=1.5 MPa)

3.2 曳力模型比较分析

图6（a）为操作压力为0.5 MPa时，不同的曳力模型计算的充分发展段气含率径向分布与实验结果的比较，与Schiller & Naumann模型相比，改进的气泡群曳力模型和Morsi & Alexander曳力模型模拟结果更接近实验值。Schiller & Naumann模型所获得模拟结果更接近实验值；而当操作压力增至1.5 MPa时，Schiller & Naumann模型和Morsi & Alexander模型的气含率分布与低压时的模拟结果基本相同，不能反映压力增加后局部气含率增大的趋势。然而，改进的气泡群曳力模型能够描述压力对局部气含率的影响，并且模拟结果与实测值符合较好，如图6（b）所示。

图7为不同操作压力下，上述3种曳力模型计算的轴向气速径向分布与实验数据的比较。冷态实验数据显示，当表观气速为0.27 m·s−1时，随着操作压力的升高，局部气相速度略有下降。从图7中可以看出，3种曳力模型对塔中心区域（r/R＜0.3）的气相速度较为接近，越靠近塔壁模拟计算值与实验值差别越大。当压力为0.5 MPa时，改进的气泡群曳力模型的计算结果比另两种模型更接近实验值。而当压力增加到1.5 MPa时，Schiller & Naumann模型和Morsi & Alexander模型计算得到气相速度与低压下的计算值相比，并没有如实验数据表现出下降趋势，反而略有增加，然而，改进的气泡群压力模型能够定性反映轴向气速随操作压力增大而减小的趋势。由于本文给出的气泡群压力系数表达式是基于平均气含率和气泡直径的冷态数据回归得到的，尚且没有考虑对气相上升速度的修正，这可能是导致中心区域以外的轴向气速模拟结果偏差较大的一个原因。在冷态实验中，气相是以气泡离散相存在的，测量的气相速度即是气泡的上升速度，而模拟出来的气相上升速度是将气相作为连续相考虑的，即作为连续介质的气相与液相充分混合，因此所计算出来的气相上升速度相较于实验测量气泡上升速度值偏小，这可能也是气相上升速度存在偏差的另一个原因。

图6 不同曳力模型对气含率径向分布的影响Fig.6 Effect of various drag model on radial gas holdup distribution(VG=0.27 m·s−1, H=2.6 m)

图7 不同曳力模型对气相上升速度径向分布的影响Fig.7 Effect of various drag model on gas rising velocity radial distribution(VG=0.27 m·s−1,H=2.6 m)

综上，改进的气泡群曳力模型能够用于不同操作压力下的鼓泡塔流场，尽管对气相上升速度分布模拟精度不够理想，但是定性能够体现压力对轴向气相速度的影响，并且高压条件下模拟得到的气含率径向分布与实验值符合较好。因此，该模型与Schiller & Naumann和Morsi & Alexander模型相比模拟准确度有所提高，基本能够描述湍动加压鼓泡塔内，操作压力对充分发展段气含率和轴向气速分布的影响趋势。

3.3 表观气速对流场的影响

表观气速是影响加压鼓泡塔内气液两相流动的关键参数，冷态实验结果表明，全塔平均气含率随表观气速提高而增大。表2列出了操作压力为0.5 MPa，表观气速为0.20～0.31 m·s−1时，全塔平均气含率模拟值与实验值的比较结果，可以看出模型对平均气含率的预报较为准确，最大误差绝对值不超过 2%，能够定量描述加压鼓泡塔平均气含率随表观气速变化的规律。

图8为上述表观气速下，模拟得到的局部气含率与实验结果的比较。冷态实验数据显示，随着表观气速的增加，局部气含率提高，气含率沿径向分布的陡峭程度基本不变。从图8中可以看出，模拟得到的气含率分布虽然略高于实验测量结果，尤其是表观气速为0.20 m·s−1时偏差相对较大，但是模拟的气含率分布趋势与实验结果较为接近。模拟和实验结果定量上的差异，可主要归结于两个方面：一方面，改进的气泡群曳力系数关联式[式（10）]中的参数是基于压差法测量得到的平均气含率的数据回归而来；另一方面，由于气泡在浮升过程中受到液相湍流涡旋的扰动，会发生一定随机摆动，造成气泡与电导探针接触点的气泡表面法向与探针角度过大，气泡易从探针尖端脱离，使得探针测量的局部气含率偏小。秦玉建[21]的冷态实验数据显示，通过电导针测量的局部气含率积分得到的截面平均气含率，通常略低于差压法的测量值，尤其当表观气速为0.20 m·s−1时，两者相差较为明显。

表2 不同表观气速下平均气含率模拟值与实验值比较Table 2 Comparison between simulation results and experiment values for different superficial gas velocity (P=0.5 MPa)

图9为不同表观气速下，轴向气相速度模拟与实验结果的比较，从图中可以看出，随着表观气速由0.20 m·s−1增加到0.31 m·s−1，实验测量的气相速度略有提高，而径向分布趋势基本相同；而模拟得到的气相速度随着表观气速增加而增大，且分布变得更为陡峭，与实测结果相差较多。

图9 表观气速对轴向气速分布的影响Fig.9 Effect of superficial gas velocity on axial gas rising velocity along radius (P=0.5 MPa, H=2.6 m)

3.4 操作压力对流场的影响

操作压力对流场的影响显著，实验表明随着操作压力的增加，气泡尺寸略有减小，气泡上升速率降低，使得气泡在流场内停留时间增加，使气含率明显提高。表3给出了在表观气速为0.27 m·s−1，压力为0.5 MPa～2.0 MPa下，平均气含率模拟结果与实验值的比较，模拟结果略低于实测值，最大相对误差的绝对值为3.3%。

表3 不同操作压力下平均气含率模拟值与实验值比较Table 3 Comparison between simulation results and experimental values under various operating pressure(VG=0.27 m·s−1)

从图10中可以看出，压力为0.5～1.5 MPa时，模拟的气含率分布略高于实验值，而当压力达到2.0 MPa时，塔中心处实验结果高于模拟值。数值模拟与实测的气含率分布曲线变化趋势基本一致，随着操作压力的增加两者的气含率均逐渐增加。但是，在部分操作压力下，模拟的局部气含率与实验结果相比还存在一定偏差。

图 11表明模拟的气相速度随着压力升高而下降，与实验结果所呈现的趋势基本一致。然而，模拟的轴向气相速度分布比实测分布更为陡峭，并且模拟结果在总体上也略低于实验测量结果。由于模型假设全塔气泡尺寸不变，而实际上湍动鼓泡塔内气泡直径存在一定空间分布，而气泡运动速度通常与其尺寸相关，使得模型虽然能够定性描述轴向气速分布随操作压力变化的基本趋势，但是在定量上还存在明显偏差。在后续工作中将引入数量衡算方程（population balance equation），模拟得到加压鼓泡塔内气泡尺度分布，并以此校正气相速度分布，提高模型的预测精度。

图10 压力对气含率径向分布的影响Fig.10 Effect of operating pressure on radial gas holdup distribution（VG=0.27 m·s−1, H=2.6 m）

图11 压力对轴向气速分布的影响Fig.11 Effect of operating pressure on gas rising velocity radial distribution（VG=0.27 m·s−1）

4 结 论

（1）本文采用双流体模型对内径0.3 m加压湍动鼓泡塔内的气液两相流动，进行了动态二维轴对称数值模拟。针对加压条件下，采用改进的气泡群曳力系数关联式，同时考虑气泡群径向力平衡机制，以及壁面润滑力建立的加压鼓泡塔流体力学模型，数值计算的稳定性和效率较高，模拟得到的气含率分布与实验结果符合较好，能够描述加压鼓泡塔内气含率和轴向气速分布随操作条件变化的趋势。

（2）与其他采用单气泡曳力系数的模型相比，改进的气泡群曳力模型能够较为准确地模拟加压鼓泡塔总体气含率和气含率分布，正确描述轴向气速分布随表观气速和压力的变化趋势。对处于加压湍动鼓泡状态的工业浆态床反应器内流场特性参数分布的预测和变化规律的认识具有重要意义。

符 号 说 明

CD——曳力系数

CD,0——单气泡曳力系数

CWL——壁面润滑力系数

CTD——湍流扩散力系数

D——塔径，m

dB——气泡直径，m

Fex——相间作用力，N·m−3

FL——气泡升力，N·m−3

FTD——气泡湍流扩散力，N·m−3

FWL——气泡壁面润滑力，N·m−3

g——重力加速度，m·s−2

H——轴向高度，m

k, kL——湍流动能, m2·s−2

Nc——计算网格数量

nW——指向壁面外的单位法向量

P——压力，Pa

P0——标况下大气压，Pa

R——塔半径，m

r——径向位置，m

u——气相或液相速度，m·s−1

uG——气相速度，m·s−1

uL——液相速度，m·s−1

VG——表观气速，m·s−1

yw——壁面距离，m

ε——湍流耗散率，m2·s−3

εG——局部气含率

εL——局部液含率

ε−G——平均气含率

ε−L——平均液含率

μt——湍流黏度，Pa·s

ρG——气体密度，kg·m−3

ρL——液体密度，kg·m−3

σ——液体表面张力，N·m−3

下角标

G——气相

L——液相

i——液相(1)或气相(2)

[1]许世峰, 王斯民, 李彩霞, 等. FT合成浆态床反应器的研究进展[J]. 化工进展, 2013, 32(S1): 4-8.XU S F, WAMG S M, LI C X, et al. Recent advances of the slurry-bed reactor FT synthesis[J]. Chemical Industry and Engineering Progress, 2013, 32(S1): 4-8.

[2]KRISHNA R,VAN BATEN J M. Eulerian simulations of bubble columns operating at elevated pressures in the churn turbulent flow regime[J]. Chemical Engineering Science, 2001, 56(21): 6249-6258.

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

[4]JOSHI J B. Computational flow modeling and design of bubble column reactors[J]. Chemical Engineering Science, 2001, 56: 5893-5933.

[5]TABID M V , ROY S A, JOSHI J B. CFD simulation of bubble column—an analysis of interphase forces and turbulence models[J]. Chemical Engineering Journal, 2007, 139(3): 589-614.

[6]LUCAS D, KREPPER E, PRASSER H M. Prediction of radial gas profiles in vertical pipe flow on the basis of buble size distribution[J]. International Journal of Thermal Sciences, 2001, 40(3): 217-226.

[7]张煜, 李红波, 李兆奇, 等. 湍动浆态床流体力学研究(Ⅳ): 带垂直列管束的浆态床流体力学模型与模拟[J]. 化工学报, 2011, 62(12): 3373-3380. ZHANG Y, LI H B, LI Z Q, et al. Hydrodynamics of turbulent slurry bubble column (Ⅳ): Modeling and simulation of bubble column with internal pipe bundles[J]. CIESC Journal, 2011, 62(12): 3373-3380.

[8]李兆奇, 王丽军, 管小平, 等. 基于径向力平衡的鼓泡塔二维流体力学模型[J]. 化工学报, 2014, 65(11):4221-4230. LI Z Q, WANG L J, GUAN X P, et al. 2-D hydrodynamic model of bubble column based on lateral-force balance[J]. CIESC Journal, 2014, 65(11): 4221-4230.

[9]DREW D A. Mathematical modeling of two-phase fow[J]. Annual Review of Fluid Mechanics, 1983, 15(1): 261-291.

[10]RAFIQUE M, CHEN P, DUDUKOVIĆ M P. Computational modeling of gas-liquid flow in bubble column[J]. Reviews in Chemical Engineering, 2004, 20(3): 225-375.

[11]DREW D A, PASSMAN S L. Applied Mathematical Sciences[M]. New York: Springer, 1999.

[12]PFLEGER D, BECKER S. Modeling and simulation of the dynamic flow behaviour in a bubble column[J]. Chemical Engineering Science, 2001, 56(4): 1737-1747.

[13]MUDDE R F, SIMONIN O. Two and three dimensional simulation of a bubble plume using a two fluid model[J]. Chemical Engineering Science, 1999, 54(21): 5061-5069.

[14]BUWA V V, RANADE V V. Dynamics of gas iquid flows in rectangular bubble columns: experiments and single/multi groups imulations[J]. Chemical Engineering Science, 2002, 57(22): 4715-4736.

[15]PFLEGER D, GOMES S, GILBERT N, et al. Hydrodynamic simulation of laboratory scale bubble columns, fundamental studies of Eulerian-Eulerian modeling approach[J]. Chemical Engineering Science, 1999, 54(21): 5091-5099.

[16]李光, 杨晓钢, 戴干策. 鼓泡塔反应器气液两相流 CFD数值模拟[J]. 化工学报, 2008, 59(8): 1954-1965. LI G, YANG X G, DAI G C. CFD simulation of gas-liquid flow in bubble column[J]. Journal of Chemical Industry and Engineering (China), 2008, 59(8): 1954-1965.

[17]李光. 鼓泡塔反应器气液两相流数值模拟模型及应用[D]. 上海:华东理工大学, 2010. LI G. Two-phase flow dynamical simulations and modelling of bubble column reactors[D]. Shanghai: East China University of Science and Technology, 2010.

[18]KHOPKAR A R, RAMMOHAN A R, RANADE V V, et al. Gas-liquid flow generated by a Rushton turbine in stirred vessel: CARPT/CT measurements and CFD simulation[J]. Chemical Engineering Science, 2005, 60: 2215-2229.

[19]LAU Y M, ROGHAIR I, DEEN N G, et al. Numerical investigation of the drag closure for bubbles in bubble swarms[J]. Chemical Engineering Science, 2011, 66(14): 3309-3316.

[20]ROGHAIR I VAN SINT ANNALAND M, KUIPERS H J A M. Drag force and clustering in bubble swarms [J]. AIChE Journal, 2013, 59(5): 1791-1800.

[21]秦玉建. 加压气液鼓泡塔流体参数的测量与CFD数值模拟[D]. 北京: 北京化工大学, 2012. QIN Y J. The measurement of hydrodynamic parameters and CFD simulation of gas-liquid flow behaviors in a pressurized bubble column[D]. Beijing: Beijing University of Chemical Technology, 2012.

[22]DREW D A, LAHEY R T. In Particulate Two-phase Flow[M]. Boston: Butterworth Geinemann, 1993: 509-566.

[23]LOPEZ DE BERTODANO M. Turbulent bubbly two-phase flow in a triangular duct [D]. New York: Rensselaer Rolytechnic Institute, 1992.

[24]ANTAL S P, LAHEY JR R T, FLAHERTY J E. Analysis of phase distribution in fully developed laminar bubbly two-phase flow[J]. Int. J. Multiphase Flow, 1991, 17(5): 635-652.

[25]TOMIYAMA A. Struggle with computational bubble dynamics[C]// ICMF'98, 3rd International Conference on Multiphase Flow. Lyon, France, 1998: 1-18.

Hydrodynamics simulation of pressurized bubble column based on bubble swarm interphase force models

LIU Xin1,2, ZHANG Yu2, ZHANG Li2, JIN Haibo1
(1School of Chemical Engineering, Beijing Institute of Petrochemical Technology, Beijing 102617, China;2Research and Development Center, Synfuels China, Beijing 101407, China)

The commercial bubble column reactors are normally operated under the conditions of pressurized state and churn turbulent bubbly flow which the superficial gas velocity always exceeded 0.1 m·s−1. However, most of the published literatures are focus on the cold model experimental investigation of the flow field characteristics in pressurized churn-turbulent bubble columns, and the empirical or semi-empirical correlations which were on the basis of the cold model experimental data to quantitatively estimate the characteristic parameters. Therefore, it is lack of the systematic study on the flow characteristics in the pressurized turbulent bubble columns by the hydrodynamic models with interactions between bubbles and liquid phase. In this work, a steady 2-D axisymmetric hydrodynamic two-fluid model of pressurized bubble column was developed on the platform of FLUENT 15.0 with the UDF sub-models which described the bubble swarm interface forces. The numerical investigations upon the flow field characteristics of pressurized churn turbulent bubble column with the operation pressure varying from 0.5 MPa to 2.0 MPa, and superficial gas velocity from 0.2 m·s−1to 0.31 m·s−1were performed. Moreover, the comparisons between the simulation results and cold model experimental data were implemented. The comparisons indicated that, the model with the bubble swarm interface forces including improved bubble swarm drag force model, lateral force balance model, and wall lubrication force model, is able to reflect the influences of operation pressure and superficial gas velocity upon flow pattern, that is, the elevatingoperation pressure leads to increase of time averaged gas holdup and decrease of time averaged axial gas phase velocity, whereas for the ascending superficial gas velocity, the gas holdup and gas phase velocity are both raising. Consequently this mathematical model could correctly predict the basic gas-liquid flow pattern in the turbulent pressurized bubble columns, and simulation results are in good agreement with the experimental data.

bubble column; gas-liquid flow; computational fluid dynamics; bubble swarm drag force model; lateral- force balance model; wall lubrication force model

Prof. JIN Haibo, jinhaibo@bipt.edu.cn

TQ 021.1

：A

：0438—1157（2017）01—0087—10

10.11949/j.issn.0438-1157.20160762

2016-06-01收到初稿，2016-10-20收到修改稿。

联系人：靳海波。

：刘鑫（1989—），男，硕士研究生。

国家自然科学基金重大研究计划项目（91634101）；北京市属高等学校高层次人才引进与培养计划项目（CIT&TCD20130325）。

Received date: 2016-06-01.

Foundation item: supported by the National Natural Science Foundation of China (91634101) and the Importation and Development of High-Caliber Talents Project of Beijing Municipal Institutions(CIT&TCD20130325).