鼓泡塔内空气-醋酸体系流体力学参数的CFD-PBM耦合模型数值模拟

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

2、数值模拟的研究仍处于发展阶段5。同时，双流体模型中破碎和聚并模型较为复杂、参数较多，使得鼓泡塔的数值模拟成为一个充满机遇与挑战的研究领域。研究者采用计算流体力学（CFD）方法模拟鼓泡塔反应器内的多尺度流动问题，模拟结果的准确性取决于模型对塔内多相湍流描述是否合理，如相间作用力的设置、气泡聚并和破碎的描述以及能量耗散机制等。Zhang等6-7考察了气泡聚并、聚并效率、破裂所需的临界能量、气泡大小等因素对CFD-BPM耦合模型数值计算的影响；Yu等8采用CFD-PBM模拟来预测分散相液滴尺寸分布，并研究脉冲塔内局部液-液流动行为；Yang等9-10考察了压力对气泡破碎的影响，并在气泡破碎模型中加入

3、了密度修正项；Gong等11提出了一种改进的模拟湍流引起气泡聚并的理论模型用来预测临界聚并速度以及气泡大小分布；王铁锋等12-14的研究认为气泡群曳力系数与单个气泡及气泡群的运动速度有关，考虑了能量和压力约束条件，并通过湍流涡体动能和气泡表面能增量的大小来判断气泡破碎的原因；Prince等15的聚并效率模型考虑了液膜排水模型的机理；Luo等16的聚并效率模型考虑了能量守恒；Lehr等17的聚并效率模型考虑了临界速度；Luo等18认为，当湍流涡旋的湍动能大于气泡表面能量增量时，气泡发生破碎，这是一种能量约束的现象；Shi等19结合这些结果，通过引入等效直径进一步修正了气泡破碎模型；Hinze20

4、研究发现破碎函数可以通过惯性力和表面张力相比来定义；Kolmogorov在Weber数的基础上提出的破碎函数，认为气泡破碎行为和Weber数有关21。因此，上述研究推动了气泡破碎模型与聚并模型在气液鼓泡塔数学模拟的发展。本文在课题组前期工作的基础上22-26，在空气-醋酸体系的鼓泡塔内对曳力模型和聚并模型进行修正，通过CFD-PBM耦合模型进行了二维和三维数值模拟，研究了醋酸浓度对鼓泡塔内流体力学参数的影响，并与实验数据进行对比，验证该模型的可行性。1 数学模型1.1 双流体模型气液系统采用欧拉-欧拉模型，控制方程如下：质量守恒方程tii+iiui=0,i=g,l(1)动量守恒方程tiiui+

5、iiuiui=-iP+ieffu+uT+Fi,j+iig(2)式中，i代表气相g或液相l；、u、t分别表示相含率、密度、速度矢量和时间；P表示压力；eff表示流体有效黏度。1.2 湍流方程选择雷诺时均法（RANS）中的标准k-模型27，方程具体描述如下：k方程kt+uk=+tkk+Gk-(3)方程t+u=+t+kC1Gk-C2(4)t=lCk2(5)式中，C1=1.44,C2=1.92,C=0.09,k=1.0,= 1.3。1.3 相间作用力模型欧拉-欧拉模型需要相间作用力来封闭，相间作用力的选择对于模拟结果的精度十分重要。本文模拟选择了曳力、升力、壁面润滑力和湍流扩散力，忽略了虚拟质量力的作

6、用。具体表达式为：Fi,j=FD+FL+FT+FW(6)1.3.1 曳力本文基于Li等28的能量最小多尺度（EMMS）模型，在不同质量分数的醋酸-空气体系中，通过表面张力函数f(l/0)修正曳力模型：CDdb=g,smalldsmallCD,smallUg,smallg,small2+g,largedlargeCD,largeUg.largeg,large2g,totalUg2(7)简化并优化后的关系式为：CDdb=fl0159-2166Ug+12359Ug2(8)fl0=7.1l0-5l02-1(9)1.3.2 升力气泡在液相剪切流中上升时，会产生与气泡上升方向垂直的一种径向力，称为升力。径

7、向升力对于气含率的径向分布至关重要，本文采用文献29的升力，具体表达式为：FL=-CLglug-ulul(10)CL=-0.2CTDl2l(11)1.3.3 壁面润滑力壁面润滑力是由靠近壁面处气液速度梯度引起的，是使气泡远离鼓泡塔壁面的一种力。本文采用文献30的壁面润滑力，具体表达式如下：FW=CWLlg(ug-ul)-(ug-ul)nr2nr(12)CWL=CWmax0,1CWD+1-yWCWCdbyWyWCWCdbm-1(13)式中，CWD、CWC为无量纲常数；m=1.52。CW为壁面润滑力系数，是关于Eo的函数，具体表达式如下：CW=0.47Eo1e-0.933Eo+0.1791Eo50

8、.00599Eo-0.01875Eo330.17933Eo(14)1.3.4 湍流扩散力湍流扩散力是由液相湍流旋涡引起的，此力使得径向气含率分布更均匀。本文采用文献31的湍流扩散力，具体表达式为：FT=CTDlklg(15)式中，kl为湍流强度；CTD为湍流扩散力系数，取值范围0.11，本文中取值为1。1.4 群体平衡模型群体平衡模型（PBM）主要是描述反应器中颗粒或气泡大小分布的一种方法。在鼓泡塔中，气泡主要发生破碎和聚并，因此该状况下群体平衡模型表示如下：nV,tt时间项+nV,tubV,t对流项+vdVbdtnV,t传质和压力变化的影响=120VnV-V,tnV,tcV-V,VdV-聚并

9、引起的源项0nV,tnV,tnV,VdV聚并引起的汇项+0V,VbVnV,tdV破碎引起的源项-bVnV,t破碎引起的汇项(16)式中，V为母气泡体积；V为子气泡体积；n(V,t)为体积为V的气泡数密度函数；c(V-V,V) 为气泡聚并速率；（V,V）为体积V的气泡破裂成体积V的子气泡分布函数。1.4.1 聚并模型气泡聚并速率可表示为：cdi,dj=cdi,djPcdi,dj(17)气泡间碰撞频率为：cdi,dj=241/3di+dj2di2/3+dj2/31/2(18)本文基于文献16的聚并模型，在空气-醋酸体系中，通过引入表观气速和表面张力修正项对聚并系数进行了修正，如表1所示。表1聚并系

10、数的修正Table 1Correction of the coalescence coefficientUg/(m/s)Concentration/%(mass)/0CeexpsimError/%0.047501.441.200.1200.113-6.190.071501.441.150.1520.1583.950.094501.441.100.1850.2022.650.047601.361.150.1390.130-6.920.071601.361.100.1740.1771.720.094601.361.050.2120.2214.250.047701.291.100.1430.141-

11、1.400.071701.291.050.1820.1872.750.094701.291.000.2200.2325.450.047801.201.050.1400.134-6.830.071801.201.000.1740.1824.600.094801.200.950.2070.2269.170.047901.111.000.1210.116-4.130.071901.110.950.1600.1653.130.094901.110.900.1950.2086.670.0471001.000.950.1170.108-7.700.0711001.000.900.1530.1551.310

12、.0941001.000.850.1930.1993.11新窗口打开|下载CSV修正的聚并系数Ce具体表达式如下：Ce=0.9Ug-0.0200.62(19)式中，0为25下标准浓度醋酸的表面张力;为实际醋酸溶液的表面张力。1.4.2 破碎模型本文主要采用Luo等18的气泡破碎模型进行数值模拟，该模型具体表达式为：bfv,d=0.92381-gd21/3min11+211/3exp-12cfl2/311/3d5/3d(20)式中，的取值为2.047；min为湍流涡和母气泡最小尺寸比，即min=min/d，其中min为Kolmogorov最小湍流涡尺寸的11.4倍；cf为气泡表面能增加量。2 实

13、验装置和模拟设置2.1 实验部分实验在玻璃鼓泡塔内进行，实验装置和测量方法详见文献32，塔内直径D=0.15 m，塔高H=2.2 m，采用差压变送器、光纤探针、ERT技术手段测量实验数据。同样塔底部用以储存气体的气室以及多孔气体分布板均采用不锈钢材料制作完成，气体分布板上共有2 mm19的开孔，开孔率约为0.338%。2.2 物理特性实验采用符合标准的工业冰醋酸（100%的醋酸），并根据其密度配制不同浓度的醋酸溶液，对所配溶液进行密度、表面张力、黏度的测量，不同浓度醋酸溶液的物性如图1所示。图1图1不同质量分数的醋酸物理性质Fig.1Physical properties of acetic

14、acid with different mass fractions2.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下，采用网格数为

15、2200、3080和4590个三种网格验证径向气含率分布的计算精度。从结果看出，径向气含率随着网格数的增加基本不发生变化，因此，在综合考虑模拟的计算精度和计算时间的情况下，本模拟中最终选择网格数为3080个的Grid网格。图2图2不同网格质量对径向气含率的影响Fig.2Influence of different mesh quality on radial gas holdup3.2 不同曳力模型对径向气含率的影响在80%醋酸的条件下，将S-N、Tomiyama单气泡曳力模型以及本文修正的曳力模型等模拟的径向气含率结果进行对比分析，如图3所示。单气泡直径的设定采用光纤探针实验测得的平均值，在

16、0.071和0.094 m/s的表观气速下，测得的气泡直径分别为7.0和6.5 mm。从图中可以看出，径向气含率随表观气速的增加而增加，从塔中心到塔壁处以抛物线形式分布，越接近塔壁处，径向气含率越小。而且可以发现壁面处的径向气含率模拟结果明显低于光纤探针测得的实验值，这是因为当鼓泡塔的塔径较小时，具有较强的边壁效应，导致塔壁面处预测值偏低。图3图3不同曳力模型和表观气速下径向气含率对比Fig.3Comparison of radial gas holdup under different drag force models and superficial gas velocities通过对三种

17、曳力模型的模拟结果进行对比，结果表明，修正后的曳力模型，直接以CD/db的形式计算曳力的大小，因此不会受到气泡尺寸大小的影响，与光纤探针测得的实验值吻合较好；而S-N单气泡曳力模型的预测结果偏低，Tomiyama单气泡曳力模型的预测结果明显偏高，与王钰等33通过基于EMMS方法的鼓泡塔反应器CFD模拟结果相一致。3.3 轴向高度对气含率的影响图4为醋酸浓度为100%，塔轴向高度为0.76 m和0.86 m处的径向气含率分布。从图中可以看出，径向气含率和轴向液速呈现出塔中心大，边壁减小的趋势。在0.76 m和0.86 m处径向气含率分布基本重合，说明在鼓泡塔的这两个轴向高度处已经为充分发展阶段，

18、与文献34中列管型鼓泡塔中流动发展规律相一致。图4图4不同轴向高度下的径向气含率分布Fig.4Radial gas holdup and axial liquid velocity distribution at different column heights3.4 空气-醋酸体系径向气含率图5为空气-醋酸体系下径向气含率分布。从图中可以看出，在塔高0.86 m处，模型的模拟结果与光纤探针测得的实验值基本吻合。在相同浓度下，径向气含率随表观气速的增大而增大。且径向气含率从塔中心到塔壁处呈现减小趋势，在塔壁处径向气含率减小幅度较大。图5图5空气-醋酸体系下径向气含率分布Fig.5Radial

19、gas holdup distribution with air-acetic acid system3.5 空气-醋酸体系平均气含率图6为不同浓度和表观气速下醋酸平均气含率分布。从图中可以看出，实验值和模拟结果在10%的误差范围内基本吻合。醋酸浓度在50%100%的范围内时，平均气含率先增加后减小，浓度在70%80%之间时，平均气含率存在最大值。醋酸浓度在50%80%范围内时，表面张力呈现减小的趋势，表面张力减小使得塔内气泡稳定性减弱，大气泡破碎成小气泡，气含率增加。醋酸浓度在80%100%范围内时，表面张力减小，气泡破碎的概率增大。但同时发现，在此醋酸浓度范围内，黏度在13 mPas的低黏

20、度范围内是逐渐减小的。Ruzicka等35研究认为黏度在03 mPas时，气含率随黏度的增加而增加。因此，在此浓度范围内，黏度减小导致气含率减小。黄娟等36通过实验研究发现，气含率在醋酸浓度为60%80%的范围内存在最大值，也验证了本文的模拟结果。图6图6不同醋酸浓度和表观气速下平均气含率分布Fig.6Distribution of average gas holdup under different acetic acid concentrations and superficial gas velocities3.6 空气-醋酸体系轴向液速图7为不同浓度和表观气速下轴向液速分布。由图可以看

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

22、uid velocity distribution with air-acetic acid system3.7 空气-醋酸体系气泡直径分布图8为不同浓度和表观气速下醋酸的气泡直径分布。由图8可看出，在相同的醋酸浓度下，随着表观气速的增加，径向气泡直径分布也随着增大。这是因为随着塔内表观气速增大，增加了塔内气泡的聚并，小气泡数量减少，大气泡数量增多，径向气泡直径分布增大。在相同的表观气速下，径向气泡直径分布受到液相醋酸黏度和表面张力的影响，黏度和表面张力的变化导致气泡发生聚并和破碎。醋酸浓度在50%100%的范围内，随着黏度增加，径向气泡直径分布相应增加，当醋酸浓度达到80%时，黏度减小，径向

23、气泡直径分布相应减小；而在此醋酸浓度范围内，表面张力在一直减小，径向气泡直径分布也相应减小。因此，醋酸浓度为100%时，其黏度和表面张力都为最小值，此时气泡直径分布最小。图8图8不同醋酸浓度和表观气速下径向气泡直径分布Fig.8Radial bubble diameter distribution under different acetic acid concentrations and superficial gas velocities3.8 空气-醋酸体系气泡数密度分布图9为不同浓度和表观气速下气泡数密度分布。由图9可看出，气泡在直径为7 mm左右所占有的比例最大，这也与图8光纤探针测

24、得的气泡直径值分布相符合。在相同表观气速下，醋酸浓度为50%100%范围内，气泡数密度变化程度不大。在醋酸浓度为50%时，气泡尺寸分布较窄，而醋酸浓度为70%和80%时，气泡尺寸分布较宽，此时黏度较大。而在醋酸浓度为100%时，表面张力和黏度等物理性质达到最小值，此时气泡尺寸分布最窄。因为随着醋酸浓度的增加，液相醋酸黏度先增大后减小，表面张力一直减小。液相黏度增加，气泡发生聚并概率增加，小气泡聚并成大气泡，气泡尺寸分布变宽；表面张力减小，大气泡变得不稳定，破碎成小气泡，塔内小气泡数量增多增加了气泡在塔内的停留时间，气含率增加，气泡分布范围变窄。在相同的醋酸浓度下，表观气速的增加使得进入塔内的气

25、体流量增加，气泡数密度增加。图9图9气泡数密度分布（z=0.86 m）Fig.9Bubble number density distribution (z=0.86 m)4 醋酸体系三维数值模拟4.1 三维网格划分本部分为不同浓度醋酸三维数值模拟，采用三维非结构性网格的划分，三维网格的网格数为93856个。与二维数值模拟相比，三维数值模拟可以得到塔截面的气含率分布状况、塔内瞬时气含率分布以及塔内流型的变化情况。因此，进行三维数值模拟对了解塔内气液相的瞬时变化情况具有一定的指导意义。4.2 三维径向截面气含率分布通过ERT实验测得的云图和三维数值模拟的云图进行了对比分析，图10为50%80%醋酸

26、的平均气含率云图分布。从图中可以看出，三维数值模拟的云图和ERT测得的云图变化规律基本一致。为了可视化不同浓度醋酸在不同表观气速下的变化规律，使用统一的图例（00.40）进行说明。图中中心处黄色区域代表较高的气含率，而边壁蓝色区域代表较低的气含率。随着表观气速的增加，图像中的黄色区域逐渐增加，蓝色区域逐渐减小。说明在相同醋酸浓度下，气含率随着表观气速的增加而增加。同时，从图中可以看出，当醋酸浓度为70时，塔中心气泡量相对较大，与黄娟等36实验结果一致。图10图10不同表观气速下不同浓度醋酸横截面气含率云图分布（z= 0.86 m）Fig.10Cross-section contours of

27、gas holdup distribution of different mass fractions of acetic acid under different superficial gas velocities（z= 0.86 m）4.3 三维柱体及轴向截面气含率分布从三维柱体图11可以看出，气含率随表观气速的增加而增加，且塔中心气含率较高，塔壁面处气含率较低。三维柱体图可以更全面地看出塔中心和塔壁面气含率沿塔轴向高度的变化趋势。也可以看出气泡在上升过程中的流动形态，即以S形螺旋式上升。图11图1180%醋酸三维柱体气含率分布Fig.11Three-dimensional cylind

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

29、分布情况，这也与塔内实际的流动情况相符。图12图1280%醋酸轴向截面气含率云图分布Fig.12Axial-section contours of gas holdup distribution of 80% acetic acid solution4.5 不同轴向高度处三维径向截面气含率分布图13为浓度80%的醋酸三维轴向截面气含率充分发展状态下的云图分布。从图中可以看出，在表观气速为0.094 m/s的条件下，不同轴向高度的气含率呈现中心高、边壁低的趋势，而且较高的气含率基本分布在塔中心，轴向高度的变化对径向截面气含率的大小分布影响不大。图13图13不同浓度醋酸不同轴向高度处径向截面气含率云图分布Fig