张文龙，宁尚雷，靳海波，马磊，何广湘，杨索和，郭晓燕，张荣月

（1 北京石油化工学院新材料与化工学院，北京 102617；2燃料清洁化及高效催化减排技术北京市重点实验室，北京 102617）

引 言

鼓泡塔反应器因操作方便、结构简单、传质和传热性好等优点，被广泛地应用于石油化工、煤化工、环境工程和食品工程等领域[1]，如费-托合成、环己烷氧化、甲醇合成和对苯二甲酸合成等[2]。近年来，我国已经成为精对苯二甲酸（PTA）的生产大国，其中对二甲苯（PX）氧化反应器是PTA 生产过程中的关键部分，其决定了PTA 的产品质量和生产速率[3-4]。反应物性与操作条件是影响鼓泡塔内流体力学的主要因素，同时其气含率和气泡尺寸分布是设计和分析鼓泡塔反应器的重要参数，因此，研究有机体系下鼓泡塔流体力学特性对反应器的设计与放大具有重要意义。

目前，关于气液鼓泡塔数值模拟的研究仍处于发展阶段[5]。同时，双流体模型中破碎和聚并模型较为复杂、参数较多，使得鼓泡塔的数值模拟成为一个充满机遇与挑战的研究领域。研究者采用计算流体力学（CFD）方法模拟鼓泡塔反应器内的多尺度流动问题，模拟结果的准确性取决于模型对塔内多相湍流描述是否合理，如相间作用力的设置、气泡聚并和破碎的描述以及能量耗散机制等。Zhang等[6-7]考察了气泡聚并、聚并效率、破裂所需的临界能量、气泡大小等因素对CFD-BPM 耦合模型数值计算的影响；Yu等[8]采用CFD-PBM 模拟来预测分散相液滴尺寸分布，并研究脉冲塔内局部液-液流动行为；Yang 等[9-10]考察了压力对气泡破碎的影响，并在气泡破碎模型中加入了密度修正项；Gong 等[11]提出了一种改进的模拟湍流引起气泡聚并的理论模型用来预测临界聚并速度以及气泡大小分布；王铁锋等[12-14]的研究认为气泡群曳力系数与单个气泡及气泡群的运动速度有关，考虑了能量和压力约束条件，并通过湍流涡体动能和气泡表面能增量的大小来判断气泡破碎的原因；Prince 等[15]的聚并效率模型考虑了液膜排水模型的机理；Luo等[16]的聚并效率模型考虑了能量守恒；Lehr 等[17]的聚并效率模型考虑了临界速度；Luo 等[18]认为，当湍流涡旋的湍动能大于气泡表面能量增量时，气泡发生破碎，这是一种能量约束的现象；Shi 等[19]结合这些结果，通过引入等效直径进一步修正了气泡破碎模型；Hinze[20]研究发现破碎函数可以通过惯性力和表面张力相比来定义；Kolmogorov 在Weber 数的基础上提出的破碎函数，认为气泡破碎行为和Weber 数有关[21]。因此，上述研究推动了气泡破碎模型与聚并模型在气液鼓泡塔数学模拟的发展。

本文在课题组前期工作的基础上[22-26]，在空气-醋酸体系的鼓泡塔内对曳力模型和聚并模型进行修正，通过CFD-PBM 耦合模型进行了二维和三维数值模拟，研究了醋酸浓度对鼓泡塔内流体力学参数的影响，并与实验数据进行对比，验证该模型的可行性。

1 数学模型

1.1 双流体模型

气液系统采用欧拉-欧拉模型，控制方程如下：质量守恒方程

动量守恒方程

式中，i代表气相g 或液相l；α、ρ、u、t分别表示相含率、密度、速度矢量和时间；P表示压力；μeff表示流体有效黏度。

1.2 湍流方程

选择雷诺时均法（RANS）中的标准k-ε模型[27]，方程具体描述如下：

式中，Cε1=1.44,Cε2=1.92,Cμ=0.09,σk=1.0,σε=1.3。

1.3 相间作用力模型

欧拉-欧拉模型需要相间作用力来封闭，相间作用力的选择对于模拟结果的精度十分重要。本文模拟选择了曳力、升力、壁面润滑力和湍流扩散力，忽略了虚拟质量力的作用。具体表达式为：

1.3.1 曳力 本文基于Li 等[28]的能量最小多尺度（EMMS）模型，在不同质量分数的醋酸-空气体系中，通过表面张力函数f(σl/σ0)修正曳力模型：

简化并优化后的关系式为：

1.3.2 升力 气泡在液相剪切流中上升时，会产生与气泡上升方向垂直的一种径向力，称为升力。径向升力对于气含率的径向分布至关重要，本文采用文献[29]的升力，具体表达式为：

1.3.3 壁面润滑力 壁面润滑力是由靠近壁面处气液速度梯度引起的，是使气泡远离鼓泡塔壁面的一种力。本文采用文献[30]的壁面润滑力，具体表达式如下：

式中，CWD、CWC为无量纲常数；m=1.5～2。CW为壁面润滑力系数，是关于Eo的函数，具体表达式如下：

1.3.4 湍流扩散力 湍流扩散力是由液相湍流旋涡引起的，此力使得径向气含率分布更均匀。本文采用文献[31]的湍流扩散力，具体表达式为：

式中，kl为湍流强度；CTD为湍流扩散力系数，取值范围0.1～1，本文中取值为1。

1.4 群体平衡模型

群体平衡模型（PBM）主要是描述反应器中颗粒或气泡大小分布的一种方法。在鼓泡塔中，气泡主要发生破碎和聚并，因此该状况下群体平衡模型表示如下：

式中，V为母气泡体积；V′为子气泡体积；n(V′,t)为体积为V′的气泡数密度函数；c(V-V′,V′)为气泡聚并速率；β（V,V′）为体积V的气泡破裂成体积V′的子气泡分布函数。

1.4.1 聚并模型 气泡聚并速率可表示为：

气泡间碰撞频率为：

本文基于文献[16]的聚并模型，在空气-醋酸体系中，通过引入表观气速和表面张力修正项对聚并系数进行了修正，如表1所示。

表1 聚并系数的修正Table 1 Correction of the coalescence coefficient

修正的聚并系数Ce具体表达式如下：

式中，σ0为25℃下标准浓度醋酸的表面张力;σ为实际醋酸溶液的表面张力。

1.4.2 破碎模型 本文主要采用Luo 等[18]的气泡破碎模型进行数值模拟，该模型具体表达式为：

式中，β的取值为2.047；ξmin为湍流涡和母气泡最小尺寸比，即ξmin=λmin/d，其中λmin为Kolmogorov最小湍流涡尺寸的11.4倍；cf为气泡表面能增加量。

2 实验装置和模拟设置

2.1 实验部分

实验在玻璃鼓泡塔内进行，实验装置和测量方法详见文献[32]，塔内直径D=0.15 m，塔高H=2.2 m，采用差压变送器、光纤探针、ERT 技术手段测量实验数据。同样塔底部用以储存气体的气室以及多孔气体分布板均采用不锈钢材料制作完成，气体分布板上共有ϕ2 mm×19的开孔，开孔率约为0.338%。

2.2 物理特性

实验采用符合标准的工业冰醋酸（100%的醋酸），并根据其密度配制不同浓度的醋酸溶液，对所配溶液进行密度、表面张力、黏度的测量，不同浓度醋酸溶液的物性如图1所示。

图1 不同质量分数的醋酸物理性质Fig.1 Physical properties of acetic acid with different mass fractions

2.3 数值模拟条件的设置

数值模拟使用ANSYS FLUENT 软件作为计算平台，气相为常温常压下的空气，其作为离散相；液相为不同质量分数的醋酸溶液，其作为连续相。对于软件程序，设置速度为鼓泡塔塔底入口边界条件，压力为塔顶出口边界条件。鼓泡塔塔高为2.2 m，初始静止液面高度为1.1 m，实验值在塔高为0.86 m 处进行采集。模拟的时间步长固定为0.002 s，并认为在80 s内达到了准稳态。

3 空气-醋酸体系二维数值模拟

3.1 网格无关性验证

图2 描述了表观气速为0.094 m/s、系统压力为101325 Pa下，采用网格数为2200、3080和4590个三种网格验证径向气含率分布的计算精度。从结果看出，径向气含率随着网格数的增加基本不发生变化，因此，在综合考虑模拟的计算精度和计算时间的情况下，本模拟中最终选择网格数为3080 个的Grid网格。

图2 不同网格质量对径向气含率的影响Fig.2 Influence of different mesh quality on radial gas holdup

3.2 不同曳力模型对径向气含率的影响

在80%醋酸的条件下，将S-N、Tomiyama 单气泡曳力模型以及本文修正的曳力模型等模拟的径向气含率结果进行对比分析，如图3 所示。单气泡直径的设定采用光纤探针实验测得的平均值，在0.071 和0.094 m/s 的表观气速下，测得的气泡直径分别为7.0 和6.5 mm。从图中可以看出，径向气含率随表观气速的增加而增加，从塔中心到塔壁处以抛物线形式分布，越接近塔壁处，径向气含率越小。而且可以发现壁面处的径向气含率模拟结果明显低于光纤探针测得的实验值，这是因为当鼓泡塔的塔径较小时，具有较强的边壁效应，导致塔壁面处预测值偏低。

图3 不同曳力模型和表观气速下径向气含率对比Fig.3 Comparison of radial gas holdup under different drag force models and superficial gas velocities

通过对三种曳力模型的模拟结果进行对比，结果表明，修正后的曳力模型，直接以CD/db的形式计算曳力的大小，因此不会受到气泡尺寸大小的影响，与光纤探针测得的实验值吻合较好；而S-N单气泡曳力模型的预测结果偏低，Tomiyama 单气泡曳力模型的预测结果明显偏高，与王钰等[33]通过基于EMMS 方法的鼓泡塔反应器CFD 模拟结果相一致。

3.3 轴向高度对气含率的影响

图4 为醋酸浓度为100%，塔轴向高度为0.76 m和0.86 m 处的径向气含率分布。从图中可以看出，径向气含率和轴向液速呈现出塔中心大，边壁减小的趋势。在0.76 m和0.86 m处径向气含率分布基本重合，说明在鼓泡塔的这两个轴向高度处已经为充分发展阶段，与文献[34]中列管型鼓泡塔中流动发展规律相一致。

图4 不同轴向高度下的径向气含率分布Fig.4 Radial gas holdup and axial liquid velocity distribution at different column heights

3.4 空气-醋酸体系径向气含率

图5 为空气-醋酸体系下径向气含率分布。从图中可以看出，在塔高0.86 m 处，模型的模拟结果与光纤探针测得的实验值基本吻合。在相同浓度下，径向气含率随表观气速的增大而增大。且径向气含率从塔中心到塔壁处呈现减小趋势，在塔壁处径向气含率减小幅度较大。

图5 空气-醋酸体系下径向气含率分布Fig.5 Radial gas holdup distribution with air-acetic acid system

3.5 空气-醋酸体系平均气含率

图6为不同浓度和表观气速下醋酸平均气含率分布。从图中可以看出，实验值和模拟结果在10%的误差范围内基本吻合。醋酸浓度在50%～100%的范围内时，平均气含率先增加后减小，浓度在70%～80%之间时，平均气含率存在最大值。醋酸浓度在50%～80%范围内时，表面张力呈现减小的趋势，表面张力减小使得塔内气泡稳定性减弱，大气泡破碎成小气泡，气含率增加。醋酸浓度在80%～100%范围内时，表面张力减小，气泡破碎的概率增大。但同时发现，在此醋酸浓度范围内，黏度在1～3 mPa·s的低黏度范围内是逐渐减小的。Ruzicka等[35]研究认为黏度在0～3 mPa·s 时，气含率随黏度的增加而增加。因此，在此浓度范围内，黏度减小导致气含率减小。黄娟等[36]通过实验研究发现，气含率在醋酸浓度为60%～80%的范围内存在最大值，也验证了本文的模拟结果。

图6 不同醋酸浓度和表观气速下平均气含率分布Fig.6 Distribution of average gas holdup under different acetic acid concentrations and superficial gas velocities

3.6 空气-醋酸体系轴向液速

图7 为不同浓度和表观气速下轴向液速分布。由图可以看出，塔内轴向液速沿塔中心向上运动，且轴向液速从塔中心到塔壁处逐渐减小，在塔壁处轴向液速为负值，说明塔壁面处轴向液速沿壁面向下运动，而且塔中心轴向液速都在0.2～0.4 m/s 范围内。在相同的醋酸浓度下，轴向液速随着塔内表观气速的增加而增加，这是因为表观气速越大，向上运动的气体对塔中心的液体施加的轴向向上的作用力越大，因此轴向液速越大。在相同的表观气速下，轴向液速随着醋酸浓度的增大基本不变。Yan等[24]研究认为，在不同的表面张力下，塔内的轴向液速基本不变。Krishna 等[37]通过实验研究发现，液体黏度对塔内轴向液速的影响作用可以忽略。

图7 空气-醋酸体系下轴向液速分布Fig.7 Axial liquid velocity distribution with air-acetic acid system

3.7 空气-醋酸体系气泡直径分布

图8为不同浓度和表观气速下醋酸的气泡直径分布。由图8可看出，在相同的醋酸浓度下，随着表观气速的增加，径向气泡直径分布也随着增大。这是因为随着塔内表观气速增大，增加了塔内气泡的聚并，小气泡数量减少，大气泡数量增多，径向气泡直径分布增大。在相同的表观气速下，径向气泡直径分布受到液相醋酸黏度和表面张力的影响，黏度和表面张力的变化导致气泡发生聚并和破碎。醋酸浓度在50%～100%的范围内，随着黏度增加，径向气泡直径分布相应增加，当醋酸浓度达到80%时，黏度减小，径向气泡直径分布相应减小；而在此醋酸浓度范围内，表面张力在一直减小，径向气泡直径分布也相应减小。因此，醋酸浓度为100%时，其黏度和表面张力都为最小值，此时气泡直径分布最小。

图8 不同醋酸浓度和表观气速下径向气泡直径分布Fig.8 Radial bubble diameter distribution under different acetic acid concentrations and superficial gas velocities

3.8 空气-醋酸体系气泡数密度分布

图9 为不同浓度和表观气速下气泡数密度分布。由图9 可看出，气泡在直径为7 mm 左右所占有的比例最大，这也与图8 光纤探针测得的气泡直径值分布相符合。在相同表观气速下，醋酸浓度为50%～100%范围内，气泡数密度变化程度不大。在醋酸浓度为50%时，气泡尺寸分布较窄，而醋酸浓度为70%和80%时，气泡尺寸分布较宽，此时黏度较大。而在醋酸浓度为100%时，表面张力和黏度等物理性质达到最小值，此时气泡尺寸分布最窄。因为随着醋酸浓度的增加，液相醋酸黏度先增大后减小，表面张力一直减小。液相黏度增加，气泡发生聚并概率增加，小气泡聚并成大气泡，气泡尺寸分布变宽；表面张力减小，大气泡变得不稳定，破碎成小气泡，塔内小气泡数量增多增加了气泡在塔内的停留时间，气含率增加，气泡分布范围变窄。在相同的醋酸浓度下，表观气速的增加使得进入塔内的气体流量增加，气泡数密度增加。

图9 气泡数密度分布（z=0.86 m）Fig.9 Bubble number density distribution(z=0.86 m)

4 醋酸体系三维数值模拟

4.1 三维网格划分

本部分为不同浓度醋酸三维数值模拟，采用三维非结构性网格的划分，三维网格的网格数为93856 个。与二维数值模拟相比，三维数值模拟可以得到塔截面的气含率分布状况、塔内瞬时气含率分布以及塔内流型的变化情况。因此，进行三维数值模拟对了解塔内气液相的瞬时变化情况具有一定的指导意义。

4.2 三维径向截面气含率分布

通过ERT 实验测得的云图和三维数值模拟的云图进行了对比分析，图10 为50%～80%醋酸的平均气含率云图分布。从图中可以看出，三维数值模拟的云图和ERT 测得的云图变化规律基本一致。为了可视化不同浓度醋酸在不同表观气速下的变化规律，使用统一的图例（0～0.40）进行说明。图中中心处黄色区域代表较高的气含率，而边壁蓝色区域代表较低的气含率。随着表观气速的增加，图像中的黄色区域逐渐增加，蓝色区域逐渐减小。说明在相同醋酸浓度下，气含率随着表观气速的增加而增加。同时，从图中可以看出，当醋酸浓度为70%时，塔中心气泡量相对较大，与黄娟等[36]实验结果一致。

图10 不同表观气速下不同浓度醋酸横截面气含率云图分布（z=0.86 m）Fig.10 Cross-section contours of gas holdup distribution of different mass fractions of acetic acid under different superficial gas velocities（z=0.86 m）

4.3 三维柱体及轴向截面气含率分布

从三维柱体图11可以看出，气含率随表观气速的增加而增加，且塔中心气含率较高，塔壁面处气含率较低。三维柱体图可以更全面地看出塔中心和塔壁面气含率沿塔轴向高度的变化趋势。也可以看出气泡在上升过程中的流动形态，即以S 形螺旋式上升。

图11 80%醋酸三维柱体气含率分布Fig.11 Three-dimensional cylinder gas holdup distribution map of 80%acetic acid solution

4.4 不同时间三维轴向截面气含率分布

在表观气速为0.094 m/s 的情况下，研究了1～30 s 时间轴向截面的气含率变化，图12 为浓度80%的醋酸在不同时间三维轴向截面气含率云图分布。从图中可以看出，随着模拟时间的增加，初始液面从1.1 m 上升到1 s 时的1.2 m，最后稳定在1.4 m。气泡流动从一开始沿塔高呈现出向上缓冲状态，到5 s时垂直贯穿整个液面高度，且呈现出中心对称分布，达到25 s 后呈现为S 形上升的趋势。而且从整体上来看，轴向瞬时气含率分布呈现出塔中心高、塔壁面处低的分布情况，这也与塔内实际的流动情况相符。

图12 80%醋酸轴向截面气含率云图分布Fig.12 Axial-section contours of gas holdup distribution of 80%acetic acid solution

4.5 不同轴向高度处三维径向截面气含率分布

图13 为浓度80%的醋酸三维轴向截面气含率充分发展状态下的云图分布。从图中可以看出，在表观气速为0.094 m/s 的条件下，不同轴向高度的气含率呈现中心高、边壁低的趋势，而且较高的气含率基本分布在塔中心，轴向高度的变化对径向截面气含率的大小分布影响不大。

图13 不同浓度醋酸不同轴向高度处径向截面气含率云图分布Fig.13 Radial-section contour distribution of gas holdup at different axial heights of different mass fractions of acetic acid solution

5 结 论

本文在冷态空气-醋酸体系中进行了CFDPBM 耦合模型数值模拟研究，通过二维和三维数值模拟方法探究了鼓泡塔内流体力学参数的变化规律，具体结论如下。

（1）通过冷态空气-醋酸体系下S-N 单气泡曳力模型、Tomiyama 单气泡曳力模型和修正的曳力模型等模拟结果的比较，S-N 单气泡曳力模型径向气含率预测值偏低，Tomiyama 单气泡曳力模型径向气含率预测值明显偏高，而修正的曳力模型具有较好的预测性。

（2）冷态空气-醋酸体系二维CFD-PBM 耦合模型数值模拟，探究了不同浓度醋酸对鼓泡塔流体力学参数的影响。通过模拟分析和实验测得的径向气含率与径向气泡直径等参数，发现模拟结果和实验值吻合较好，说明修正后的模型具有较好的预测性。通过液体黏度和表面张力对平均气含率的影响发现，在50%～100%醋酸浓度范围内，平均气含率先增加后减小，当醋酸浓度为70%～80%时，平均气含率达到最大值。

（3）冷态空气-醋酸体系三维CFD-PBM 耦合模型数值模拟结果表明，模拟测得的径向截面气含率云图分布和ERT 测得的云图分布结果基本相符，说明修正后的模型具有良好的预测性。不同时间的轴向截面气含率云图分布展示出了初始时间段塔内气含率的变化规律，而塔高的变化对径向截面气含率分布的影响作用很小，且气泡在塔内呈现S 形螺旋上升状态。

符 号 说 明

CD——曳力系数

Ce——聚并系数

CL——升力系数

CTD——湍流扩散力系数

CW——壁面润滑力系数

D——鼓泡塔内径，m

db——气泡直径，mm

Eo——Eötvös数

FD——曳力，N

FL——升力，N

FT——湍流扩散力，N

FWL——壁面润滑力，N

fv——气泡破碎比

H——鼓泡塔高度，m

H0——静液面高度，m

Pc(di,dj)——尺寸为di和dj的气泡间的聚并效率

r/R——径向位置

T——液体温度，℃

Ug——表观气速，m/s

ul——轴向液速，m/s

α——气含率

ɛ——湍流耗散率，m2/s3

ζ——气泡相对直径

ζmin——气泡最小相对直径

µl——液体黏度，Pa⋅s

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

ρl——液体密度，kg/m3

σl——液体表面张力，N/m

ϖc(di,dj)——尺寸为di和dj的气泡间的碰撞频率，m3/s