2.2 旋转填充床内流体力学特性

2.2.1 液体流动模型

由于旋转填充床内的液体流动与传统填料塔内的液体流动有很大不同,因此,要想了解旋转填充床内液体的流动特性,就必须从理论分析入手,并做合理的简化,即通过建立描述床内液体流动形态的物理模型,建立能描述液体流体力学特性的数学模型[1,3,7]

电视摄像与高速频闪照相的实验结果表明,在低转速下(约小于60g),填料内的液体主要是以填料表面上的液膜(单面膜)和覆盖填料孔隙的液膜(双面膜)两种状态存在;而在高转速下(约大于100g),液体主要是以填料丝上的膜与空间的液滴两种形态存在,并伴有少量的液丝。由于旋转造成高的剪切力,使液体在填料上的膜很薄,仅有几十微米或更小,所以,雷诺数Re很低(通常小于30),液体在填料上的流动可按层流处理。虽然局部位置上填料丝的位向随机性很强,但其上的液体流动均可简化看作呈膜状沿填料的轴向与周向流动的合成。基于以上分析,可建立液体在填料内流动的物理模型。

2.2.1.1 液体流动的物理模型

液体流动物理模型的假设如下:

①填料内的液体分为填料丝上的液膜和空间的液滴两部分;

②填料上的液膜流动分为沿填料丝轴向的降膜与沿周向的绕流两种;

③填料上的液膜流动为层流;

④液滴通过下层填料丝时即被捕获重新形成新的液滴。

2.2.1.2 液体流动的数学模型

利用简化的Navier-Stokes方程,可得到液体轴向降膜与周向绕流流动(图2-9)的速度方程,分别为:

 (2-1)
 (2-2)

式中 uz——液体在填料丝上的轴向流速;

 uθ——液体在填料丝上的周向流速;

 ν——运动黏度;

 ω——角速度;

 Ri——第i层填料半径;

 δ——液膜厚度。

图2-9 液体轴向降膜与周向绕流流动的示意图

2.2.2 液膜厚度

Munjal[8]、竺洁松[9]等人基于对旋转圆盘和叶片上液膜厚度进行分析研究来估计填料上的液膜厚度。Munjal等人的结果为:

 (2-3)

式中 Qw——单位宽度表面上的液体流量;

 ν——动力黏度;

 R——转子半径;

 ω——角速度。

预计的液膜厚度为8×10-5 m(600r/min)和5×10-5 m(1200r/min)。郭锴[1]将干湿填料情况下的电视摄像结果,用图像分析仪进行分析对比,得到了不同条件下的填料上的液膜厚度的定量数据,结果见表2-1。

表2-1 丝网填料的丝上液膜厚度的测定和分析结果

郭奋[7]结合前面的式(2-1)、式(2-2)和表2-1中的实验结果,得到:

 (2-4)

式中 ν——动力黏度;

 af——填料比表面积;

 ω——角速度;

 R——转子半径。

利用上式,可计算丝网填料上的平均液膜厚度。

2.2.3 液滴直径

床层半径R处填料丝网空隙中的最大液滴直径可由液滴的受力分析得出(参见图2-10)。

图2-10 液滴受力分析

由于对丝网上的液体,离心力是主要的支配力,故忽略重力及气体曳力的影响。此时,液滴主要受两个力作用,一个是离心力:

 (2-5)

另一个是表面张力:

 (2-6)

能够维持不被离心力撕碎的最大液滴直径是由上述两个力平衡决定,即:

 (2-7)
 (2-8)

设液滴的平均直径为:

 (2-9)

式中 dmax——最大液滴直径;

 R——转子半径;

 ω——角速度;

 ρ——液体密度;

 σ——表面张力;

 B——常数,可由照相分析结果得出。

液滴直径的测量结果如表2-2液滴直径测试结果所示。

表2-2 液滴直径测试结果

由表2-2可以计算出,液滴平均直径为最大直径的1/3~1/4。

2.2.4 液体在填料中的平均径向速度

设液体在填料中的平均径向速度与Lω2R的关系为:

 (2-10)

用表2-1中数据拟合可以得到:

 (2-11)

式中 L——液体通量;

 ω——角速度;

 ab——系数;

 R几何平均半径,,R1为转子内半径,R2为转子外半径。

2.2.5 持液量

Basic和Dudukovic[10,11]用电导的方法对填料层的持液量进行了研究。通过建立简化的物理模型,得到了计算径向平均持液量的数学模型。郭锴[1]也用电导的方法对床层中的持液量进行了研究。结果表明,持液量随液量的增加而增大,随转速的升高而减小。

持液量与平均径向速度并不是独立的,二者有如下关系:

 (2-12)
 (2-13)

式中 εL——持液量,液体体积/填料体积;

 L——液体通量;

 u——平均径向速度;

 R——转子半径;

 ω——角速度。

杨宇成和陈建峰等[4]通过X-ray CT的技术定量测定了液体在泡沫镍填料(the nickel foam packing)和丝网填料(the wire mesh packing)内的持液量(εL)。图2-11表示转速(N)对填料区持液量的影响规律。由图可知,当转速从500r/min提高到1000r/min的过程中,持液量有快速的下降趋势,之后随着填料转速的增加,持液量下降速度变缓。图2-12表示了流量(QL)对填料内持液量的影响。持液量随着液量的增加而逐步升高;在相同转速下,随着液量的增加,持液量的增加速度比较缓慢,这意味着液量对持液量的影响并不明显。这一结果与前人的研究结论基本一致。此外,泡沫镍填料的持液量要大于金属丝网填料。

图2-11 转速对持液量的影响

图2-12 液体流量对持液量的影响

图2-13表示液体黏度(μ)对填料层内持液量的影响规律。由图可知,在高转速下,液体黏度对持液量的影响程度明显低于低转速下的。这可能是由于高转速下填料层内液体获得较大的离心力,提高了剪应力,使得液膜变薄,流速增快,从而降低了高转速下黏度对持液量的影响趋势;在同一转速下,金属丝网填料的持液量随黏度增加速度要小于泡沫镍填料。

根据实验数据,分别对两种填料建立相应的经验关联式:

 (2-14)

图2-13 液体黏度对持液量的影响

 (2-15)

式中 Ka——卡皮查准数=μ4g/σ3ρ

 Ga——伽利略准数=gd3p/ν2

2.2.6 液膜在填料丝网上流动的Re计算

由式(2-1)、式(2-2),可计算出轴向降膜和周向绕流流动的最大流速及平均流速。例如,在条件R1=0.032m,R2=0.082m,Ri=0.082m,H=0.1m,af=457m-1r0=0.00013m,L0=2m3/h,N=1500r/min,18℃,水,ν=1.066m2/s下,得到如下计算结果:

从计算结果可看出,由于液膜非常薄,填料丝网上的液膜流速很低,所以雷诺数Re很小,说明液体流动模型中关于丝网上的液膜流动为层流的假设合理。

2.2.7 液泛

Munjal[8]、朱慧铭[12]和王玉红[13]等对逆流旋转填充床中的液泛特性进行了研究,结果表明,液泛气速随液量的增加而减小,随转速的增加而增大。总体而言,逆流旋转填充床的液泛线速要比填料塔中的整砌拉西环的液泛线高40%左右。

2.2.8 气相压降

旋转填充床的压降主要由外腔(转子外缘和机器外壳之间)、转子和转子内腔(转子内缘到气体出口)三部分组成。压降主要是在转子的填料和内腔中。一般情况下,旋转填充床的压降在600~1500Pa左右。或者说,整体压降不高于传质效果与之相当的填料塔或筛板塔。

Keyvani[14]、Kumar[15]、朱慧铭[12]和沈浩[16]等人对逆流旋转填充床气相压降进行了实验研究和模型计算。总的趋势是气相压降随转速和气量的增加而增加,而在一定的液量范围内,随液量的增加而减小。特别值得一提的是,郑冲[17]等人对逆流旋转填充床中气相压降做了系统的研究工作,他们在实验中发现了三个重要的与传统填料塔完全不同的现象。第一,在实验范围内,湿床压降小于干床压降;第二,在固定填料外径的前提下,填料越厚,压降越小;第三,转子动力消耗随气体流量的增加而减小。基于简化的Navier-Stokes方程,对旋转填充床中的压降进行了数学模拟,得到了计算压降的半经验式,通过数学模拟,以上三个现象得到了初步的解释,如图2-14和图2-15所示。

图2-14 气体切向速度随半径的变化关系

图2-15 压降随半径的变化关系

李树华[18]用五孔探针测试了转子内腔流场的压力和速度,并利用氨气、氮气和二氧化碳三种气体,测试了密度对旋转填充床压降的影响。结果表明,旋转填充床内腔的速度和压力场是轴对称场;气体压力变化明显区域是气体流道突变区;减少气体压降应避免流道突然改变,从而减小气体压降和能量损失;旋转填充床压降与气体密度成正比。

经理论分析和实验数据拟合,发现转子内腔气体压力沿半径的分布符合下式:

 (2-16)

式中,abn是轴向位置的函数。因此,腔内气体压力分布是轴向和径向位置的函数。

实验测定的腔内流场情况见图2-16。图2-16中所示的2-2和3-3截面处,主流的外面有涡流区存在。流体在此构成环流,它的能量来自与主流的动量交换,这部分能量消耗在流体内部和涡流与器壁的摩擦上,最后变成热量,成为涡流损失。

图2-16 旋转填充床内腔气体流线示意图

2.2.9 旋转填充床内气液流动的CFD模拟

近年来,随着电子计算机及相关技术的高速发展,计算

流体力学(CFD)已逐渐成为流动、传递与反应器模拟研究和结构优化的有效手段。通过建立旋转填充床内气液流动的CFD方法,模拟研究旋转填充床内气相和气液两相流动过程,揭示填料层流体力学特征,对旋转填充床内结构进行优化,可为拓展旋转填充床的工业应用提供理论依据和技术支持。

2.2.9.1 旋转填充床气相流动的三维CFD模拟[19]

采用多孔介质模型模拟填料区,建立了旋转填充床气相流场的三维CFD模型,分别采用两种结构差别较大的旋转填充床压降实验数据进行模型验证。旋转填充床的几何结构如图217所示,使用四面体网格,按照T-grid方式对几何模型进行网格划分。网格间距设为5,孙润林等[20]使用的几何模型的网格数为3062309个,而Liu等[21]的几何模型的网格数为656555个。

图2-17 旋转填充床的三维几何模型

连续性方程如下:

 (2-17)

式中,ijk代表三维坐标的三个方向,Sm为连续性方程源项。对于动量守恒方程,在RPB的模拟中可以分成两个部分,对于不旋转区域的动量守恒方程如下:

 (2-18)

式中,p为静压力,SF是动量守恒方程的源项,可以添加其他体积力。采用旋转模型——

MRF方程去模拟旋转区,把动量守恒方程中的速度参数u进行了如下的修改:

 (2-19)

式中,用于代表移动旋转速度参数vr,可由下式求得:

 (2-20)

式中,vt代表移动区域内的平移速度,由于无平移,vt=0。ω代表角速度,用于模拟旋转速度。把修正后的旋转速度ur代入到动量守恒方程(2-18)中即可模拟该计算域的旋转效应。

采用多孔介质模型模拟填料区,将其添加在旋转区域动量守恒方程的源项中,方程如下:

 (2-21)

方程由两部分组成,等式右边第一部分为黏性阻力项,第二部分为惯性阻力项。这里u代表的是流动区域的速度,KpermKloss分别代表得是多孔介质的渗透性和能量损失系数,采用半经验方程——Ergun方程[22]进行求解,方程如下:

 (2-22)

动量守恒方程中的雷诺应力项采用Boussinesq的假设进行求解,该假设认为湍流中雷诺应力与时均速度成正比,该项可由下式得出:

 (2-23)

式中,δij为Kronecker符号,其数值为:

 (2-24)

式中,k为湍动能,可表示为:

 (2-25)

式中,μt是湍流黏度,与湍流运动状态有关。采用k-ε两方程去求解湍流情况,在两方程的

假设中,湍流黏度由湍动能k和湍动能耗散率ε这两个方程去获得,求解方程如下:

 (2-26)

把式(2-26)的假设代入到公式(2-23)中,即可得到两方程假设下的雷诺应力求解方程:

 (2-27)

式中,Cμ为无量纲系数,设为0.09。采用标准k-ε模型,该模型因具有稳定、经济且合理

等优点,一直被用于模拟多种工况。模型中两个方程的表示分别如下:

 (2-28)
 (2-29)

式中,常数σk=1,C=1.44,C=1.92和σε=1.3。

气体入口为速度进口,方向与进口管截面垂直,按照真实管径大小设定入口处的水力直径,并且估算湍流强度的大小;气体出口为压力出口,其压力大小与大气压相同;壁面为无滑移的固壁边界,填料区与空腔区的交界面为“interface”。图2-18和图2-19分别描述了气量(G)和转速(N)对旋转填充床压降(ΔP)的影响规律。结果表明,压降随着气量或转速的增加而增加。由图2-18(a)和图2-19(a)可知,采用标准k-ε模型比realizable k-ε模型得到的模拟结果更接近实验值;随着转速增大,采用realizable k-ε模型的模拟值逐步与实验值误差缩小;在本模拟的操作参数范围内,标准k-ε模型更合适。

图2-18 气量对压降的影响

图2-19 转速对压降的影响

在工业应用中,旋转填充床常采用径向入口方式。这种结构可导致气体在填料外缘存在严重的分布不均现象,如图2-20(a)所示。为此,提出在填料外缘和气体入口间增加挡板解决此问题。模拟结果表明,该方法能有效地改善填料外缘气体分布不均状况[图2-20(b)]。为实现优化设计,我们采用下式定量分析气体分布情况:

 (2-30)

式中,为参与计算的网格数,为选取的角度范围(分析的区域为整个填料外表面,所以),h是填料的厚度,urθiz)是单个网格上的径向速度,而是整个填料外表面的平均径向速度。的值越大就代表气相分布越不均匀,相反值越小,就代表气相分布越为均匀。

图2-21和图2-22分别表示挡板宽度D和挡板离轴距离L对填料外缘气相分布均匀性的影响规律。由图2-21可知,最低点在D=200~250mm的范围内,表明较宽的挡板要优于较窄的挡板;由图2-22可知,最低点在L=250.5mm,表明此处气体分布最为均匀。

图2-20 径向速度云图

图2-21 挡板宽度D对的影响

图2-22 挡板离轴距离L对的影响

不同开孔率曲面挡板的结构如图2-23所示,其开孔率对的影响规律如图2-24。结果表明,开孔后的曲面挡板对气相分布的效果要优于不开孔的曲面挡板;三种不同开孔率的曲面挡板对气相分布的优化效果相近。这为气体入口管的优化设计提供了思路。

图2-23 具有不同开孔率的曲面挡板结构图

孔排布:(a)ζ=10%,5行11列;(b)ζ=20%,7行15列;(c)ζ=30%,7行23列

图2-24 曲面挡板上开孔率ζrMfB(h,)的影响

2.2.9.2 旋转填充床气液两相流动的三维CFD模拟[23]

旋转填充床内液体经分布器喷到转子内缘后,在填料区受到离心力和高速旋转填料的撞击切割,形成液膜、液线、液滴等流动形态,液体为分散相,气体为连续相,其相含率(持液量)在10%以内。通过在外空腔区外壁面添加液封装置,实现了气液逆流过程的三维模拟。旋转填充床的几何结构及局部放大如图2-25和图2-26所示。采用边长为1mm的正方体模拟填料,由空隙率为0.95,计算出正方体的数量。采用真实尺寸的旋转填充床网格数量太大(达上千万),因此模拟的旋转填充床在保证径向尺寸不变的基础上,按比例大幅度缩小了轴向高度,具体模型参数如表2-3所示。模型采用间距为0.5mm的六面体网格进行网格划分,网格总数为986504个。

图2-25 三维旋转填充床几何结构

1—内空腔区;2—填料区;3—外空腔区

图2-26 局部放大图

表2-3 旋转填充床的几何尺寸

对于气液两相流动,连续性方程在单相的基础上添加了相含率α,来分别描述各相的连续性问题,具体方程如下:

 (2-31)

式中,下标q可以被写为g和l,分别代表了气相和液相,由于没有涉及传质反应,所以,源项值为0。采用欧拉模型,气液两相的动量方程是分别计算的,对于每一相态g或l,其在不旋转区域和旋转区域的方程表示如下:

 (2-32)
 (2-33)

式(2-33)通过添加来表示旋转区域的速度,而相与相之间信息,除了压力外,其他主要依靠源项SF来实现,如相间作用力:

 (2-34)

式中,Kgl(=Klg)是相间动能转换系数,对于流体间的模拟,可由公式(2-35)求出:

 (2-35)

式中,Ai表示相间面积,它在欧拉模型中的计算如下:

 (2-36)

式中,dp代表离散相的直径;αp代表p相所占的体积分率。而公式(2-35)中,τp代表离散相的松弛时间,被定义为:

 (2-37)

公式(2-35)中的f代表的是曳力函数,本部分使用常用的Schiller and Naumann提出的模型[24]进行求解:

 (2-38)

式中,CD为曳力系数,根据相间相对雷诺数Re的不同,可以分别求解:

 (2-39)

而雷诺数近似用下式计算:

 (2-40)

源项SF中还包含了表面张力项FV,采用由Brackbill提出的Continuum Surface Model(CSF)[25]进行求解:

 (2-41)

其中,γ是表面张力系数,设置为常数73.5dyn/cm(1dyn=10-5 N),为某相态相含率的变化梯度,K是表面曲率值,它可以由公式(2-42)计算得出:

 (2-42)

当液体出现在填料表面时,由于此时会出现气液固三相,需要对表面曲率值K进行修正,方法是引入接触角,根据杨式公式进行修正。

 (2-43)

式中,θw为接触角,设置为常数53°,分别代表了相对于壁面的切向和垂直方向的单位矢量。最后,动量方程式(2-32)和式(2-33)中脉冲速度继续沿用气相模拟的方法,采用两方程中标准k-ε模型进行求解。

气体和液体的入口都设置为速度进口,气液入口速度采用保证实验与模拟的流通通量相等。气体和液体的出口边界条件被设置为压力出口,出口表压设置为0Pa。整个模拟过程中操作压力被设置为标准大气压101325Pa。壁面对气相和液相都设置为无滑移的固壁边界,接触角设置为53°。压力-速度耦合问题通过Phase Coupled SIMPLE方法进行求解,网格梯度采用Least Squares Cell Based方式进行计算,采用一阶迎风方法进行离散。每个工况的时间步长设置为5×10-4s,计算残差值在1×10-3以下。模型模拟结果与实验结果的比较如图2-27~图2-30所示。

图2-27和图2-28分别表示转速和液量对旋转填充床湿床压降的影响规律。由图可知,压降随着转速和液量的升高而升高;实验值与模拟值在湿床压降的对比中,误差值要明显大于±15%。这可能是由于模拟中填料的结构尺寸与真实丝网填料相差太远导致的。为了减少计算量,模拟中填料层轴向厚度只有4mm,气相流动受到液相和壁面的影响更明显,立方体填料对液体的剪切能力弱于真实丝网填料,导致填料层内液体微元的尺寸增加。

图2-27 转速N对湿床压降ΔP的影响

图2-28 液体流量Ql对湿床压降ΔP的影响

图2-29 液体流量Ql对持液量εl的影响

图2-30 转速N对持液量εl的影响

图2-29和图2-30分别表示液量和转速对持液量的影响规律。结果表明,填料区持液量随着液量的增加而升高,随转速的降低而降低。持液量实验关联式计算值[4]与模拟值误差在±25%以内,表明本CFD模型相对合理。