膜生物反应器内气液与污泥三相流动的数值模拟

李会强, 刘 峤, 王 青, 徐 农, 3, 董 强, 3, 丁爱琴, 范 茏*

(1. 合肥学院, 能源材料与化工学院, 合肥 230601; 2. 中国科学院过程工程研究所, 生化工程国家重点实验室, 北京 100190 3. 南京工业大学, 材料化学工程国家重点实验室, 南京 211816)

活性污泥会沉积在膜生物反应器(MBR)的超滤膜表面或进入膜孔内部造成堵塞,导致膜运行通量下降[1].为了缓解这一现象,通常向膜池内曝气[2].如何在控制膜污染的前提下,提高曝气效率,降低运行能耗,是目前MBR应用过程中亟待解决的问题[3].

计算流体力学(CFD)是研究流体力学特性的有效方法,在MBR中有着广泛的应用[4].通过对MBR内部流动的CFD模拟,得到优化条件,从而减少膜污染,为实际工程设计提供理论依据和数据支持.如Ndinisa等[5]模拟了浸没式膜单元内的两相流,通过整体优化膜单元来减少污垢的形成.林进等[6]利用流体体积模型模拟气升式陶瓷膜中的气液两相流,初步提出气含率的变化规律.Brannock等[7-9]模拟了平板膜和中空纤维膜的MBR的混合特性.Wu等[10]分步模拟了MBR内的气液两相流,从层流到湍流,再到非牛顿流体、多孔介质模型逐步耦合完成.Wang等[11]对平板膜生物反应器进行了气液两相流模拟,得到了优化的通道间隙和曝气速度.Cao等[12]对全尺寸的商用MBR进行了气液两相流模拟,发现挡板的存在以及挡板与最外层膜板之间的距离对全尺寸MBR的面积加权剪切力影响较小,曝气量的增加对面积加权剪切应力有显著影响.曹迎晨等[13]对单个气泡和气泡群在MBR平板膜内的流体力学特性模拟,得到单个气泡上升速度与直径的关系,及气泡群时,气液两相剪切力的变化.Liu等[14]通过模拟全尺寸MBR内的气液两相流动,对膜间距、曝气设计及操作参数进行了优化设计.以上这些文献[5-14]的研究中,均将真实的气-水-污泥三相体系简化为气-液两相,这一定程度上能够描述MBR内的流动,但无法得到污泥的流体力学特性,及其对膜污染的影响,无法深入地理解MBR内复杂的流动现象.MBR文献中关于气-液-固三相的模拟较少,仅Bayat等[15]采用欧拉-欧拉模型对气升式MBR进行了两相和三相模拟,随着空气流速和MLSS的增加,膜污染会降低.Kang等[16]将欧拉多相流模型和多孔介质相结合,模拟比较了中试单个膜组件和工厂实际运行中多个膜组件内气、液、生物质的流体力学特性.Amini等[17]用欧拉-欧拉多相流模型对MBR进行了两相和三相模拟,将混合液作为牛顿流体,而膜作为固体壁面处理.

本文将污泥作为单独的一相处理,对MBR内的气-水-污泥三相体系进行模拟,考察了污泥本身的流动特性,探讨了气泡、曝气位置等对污泥沉降的影响,以期优化结构,减少膜污染的发生.由于缺乏三相的数据,文中还将气液两相流模拟结果与文献比较验证.

1 数学模型和模拟方法

1.1 物理模型

本文采用与Wu等[10]相同的MBR膜单元,考虑到膜结构的对称性,选取结构的1/4进行模拟(如图1).采用ANSYS WORKBENCH建立模型,整体模型为1 525×2 350×4 000 mm3的长方体,其中膜单元为625×1 350×2 000 mm3的长方体.

1.膜单元;2.气体出口;3.对称面;4.平面x=225 mm;5.曝气管图1 流体计算域的1/4Fig.1 One quarter of the MBR

膜单元区域采用六面体网格,膜单元以外区域采用四面体网格.曝气面的尺寸非常小,因此对曝气面进行网格加密,以保证计算精度.网格计算域如图2所示.

图2 计算域网格(黑框为膜单元区域,其余的为膜单元外区域)Fig.2 Computational domain grid

1.2 多相流模型和湍流模型

采用软件FLUENT 2020R2对MBR的流动进行非稳态模拟,对气-液两相流或气-液-污泥三相流均使用欧拉模型.水为连续相,污泥和气泡均为分散相,污泥密度为1 100 kg/m3,体积浓度为1%,污泥颗粒假设为刚性小球,粒径40 μm.

采用张鹏等[18]提出的式(1)计算污泥黏度,温度25 ℃,获得污泥黏度为83.18 mPa·s.

(1)

式中:μ为黏度,Pa·s;T为绝对温度,K.采用标准 湍流模型处理湍流,如式(2)、式(3)所示.

(2)

(3)

式中:Cε1,Cε2,σk均为常数,分别取值1.44,1.92,1.3.rα为α相体积分数;ρα为α相密度,kg/m3;kα为α相湍动能,J;t为时间,s;Vα为α相速度,m/s;μα为α相黏度,Pa·s;μtα为α相涡黏性系数;Pkα为α相湍流能量生成项;εα为α相湍流耗散率.

1.3 非牛顿流体模型

为了与文献中的两相模拟结果进行比较验证,将MBR内的气-水-污泥三相体系简化为气-液两相流动处理,其中液体为水和活性污泥组成的混合物.活性污泥的存在使得混合物表现为非牛顿流体,因此在模拟气-液两相流时,将液体作非牛顿流体处理,液体密度取997 kg/m3.

Liu等[19]研究发现,Ostwald-de Waele模型式(4)适合描述污泥流动特性,其剪切力与污泥浓度、剪切速率之间的关系可以用式(5)和式(6)描述.

τ=Kγn(0

(4)

τ=0.043 1MLSS0.89γ0.68(3≤MLSS≤8)

(5)

τ=0.041 2MLSS1.64γ0.45(8≤MLSS≤16)

(6)

因此,文中采用Ostwald-de Waele模型来处理液体的非牛顿流体特性,K为稠度系数;n为流动特性指数;τ为剪切力,N/m2;γ为剪切速率,s-1;MLSS为污泥浓度,g/L;污泥浓度取10.2 g/L.

1.4 多孔介质模型[20]

MBR中膜单元一般多采用的是中空纤维膜,膜单元由几十个膜组件组成,每个膜组件又包含上千根膜纤维,因此不可能对每根纤维内的流动进行模拟.紧密排列的膜纤维与多孔介质具有相似的特征,因此文中将膜单元简化为多孔介质处理.

由于膜单元类似于管束,黏性损失项可忽略不计,因此动量方程中应添加源项[10],即式(7)～(9):

(7)

(8)

(9)

1.5 边界条件及模拟方法

图1所示的模型结构中,计算域有两个对称面,膜单元设置为多孔介质区域,多孔介质区域有两面与对称面重叠,作为对称面处理,剩余两面代表膜单元的壁面.膜单元下方有6根曝气管,简化为6个宽度为2 mm,距离膜单元底面350 mm的曝气面,作为速度进口,曝气速度为3 m/s,气泡直径为6 mm,假设气泡为刚性小球,不会破裂或合并.

在气-液两相模拟时,顶部设置为“排气”,即只能气体从顶部流过,液体不会流出.在气-液-污泥三相模拟时,无法使用“排气”条件,因此在顶部左侧设置一小口,为压力出口.在进行网格划分时,将出口网格进行局部加密,尺寸设置为5 mm.

模拟时,将多相流模型、湍流模型、非牛顿流体模型、多孔介质模型等一次性加入到本文模型中进行求解.求解时,控制方程的压力-速度耦合相采用“Phase Coupled Simple”算法,动量、体积分数、湍流动能、湍流耗散率均设置为一阶迎风格式.

1.6 网格无关性验证

计算之前首先对网格的无关性进行研究,以确定合适的网格,保证计算结果的准确性.对模型进行了4种不同的网格划分,网格数量分别为584 109、747 797、917 811、1 025 309,考察平面z=-480 mm上y=750 mm处的液体流动速度.

如图3所示,网格数从584 109增加到917 811的过程中,液体流速改变较大,而网格数从917 811增加到1 025 309时,液体流速平均误差约为1.1%.因此,网格数为917 811时,计算结果已经足够准确,故选取917 811的网格进行计算.

图3 不同网格数下平面z=-480 mm上y=750 mm处液体流速变化图Fig.3 Variation of liquid velocity with grids at y=750 mm in the plane z=-480 mm

计算时发现,MBR内的流动在50 s后基本达到稳态,为了确保结果准确,下文均选择100 s时的结果进行分析.

2 结果与讨论

2.1 两相流模拟的对比与验证

目前国内外关于MBR内气-液-污泥三相流模拟的文献非常少,已有的文献[15-17]由于结构参数、操作条件等的差异以及数据的欠缺,无法具体比较.因此,将MBR内流动简化为气液两相来模拟,与文献的气液两相流结果比较验证.本文选取Wu等[10]的MBR进行模拟.Wu等[10]在进行气液两相流模拟时,分四步进行,以层流模拟为起点,在层流模拟结果的基础上,加入湍流模型,再依次加入非牛顿流体和多孔介质,逐步模拟,这样操作耗费了大量的时间.本文直接将湍流模型、非牛顿流体和多孔介质模型一次性纳入模拟中,相比于Wu等[10],操作简单,时间大幅度减少.

本文结果与Wu等[10]的比较如图4和图5所示.图4为MBR截面4(图1中所示)的液体速度矢量图,本文模拟结果与Wu等[10]的模拟结果基本一致.图5为不同结构时,膜单元内液体平均速度的比较,平均误差为25%,模拟精度在其它文献记载的范围之内(8%～30%)[22],证明本研究所建立的模型及采用的模拟方法是可靠的.

图4 液体速度矢量图对比Fig.4 Comparison of liquid velocity vectors

图5 膜单元内液体平均速度的比较Fig.5 Comparison of average liquid velocity in membrane unit

2.2 气-液-污泥三相流的模拟结果

2.2.1三相流与两相流的流动对比

三相流时平面4(图1所示)处的液体速度矢量如图6(a)所示,在膜单元下方,水随着曝气向上流动,部分液体到达膜单元底部时会绕过单元向上流动,这与两相流模拟时液体混合物随着曝气流入膜单元[图4(a)]有所不同,三相条件下模拟更为真实,污泥会堆积在膜单元底部,而不是进入膜单元,因此成为水流的阻力,使部分水流改变方向,绕过膜单元向上运动.

图6 三相流在平面4处的速度矢量图Fig.6 Velocity vectors at plane 4 for three-phase flow

图6(b)为污泥的速度矢量图,在MBR内大部分区域污泥的流动与水基本一致,速度差别不大.但由于污泥无法进入膜单元,因此在膜单元内没有固体流动.在膜单元的右上方,固体最不易达到,同样没有固体.

选取膜单元本身和膜单元外部这两个区域进行探究,将三相流时两区域的液体平均速度与两相流时进行比较,结果如图7所示.三相流动时,膜单元内外区域液体的平均速度差异较大,而两相流时差异较小.膜单元外部区域,三相流动时的液体速度比两相流时大,而膜单元内部,则远小于两相流.产生这种现象的主要原因是对污泥的不同处理方式造成的.三相流时,水和污泥为单独的两相,膜单元外部,水和污泥同时流动,膜单元内部仅有水通过,内外区域的差别不仅体现在膜单元的影响,还有污泥和水之间的相互作用.而两相流时,液体是水和污泥的混合物,同时在膜单元内外流动,因此三相流动的模拟更接近于真实的MBR,也说明对于MBR的模拟而言,本研究将污泥单独处理优于传统的两相流.

图7 膜生物反应器内各区域液体平均速度对比Fig.7 Comparison of average liquid velocity in MBR

图8给出污泥浓度分布,部分污泥会沉积在膜单元底部,这与Amini等[17]的三相模拟结果相类似,他们也发现固体颗粒更易于聚集在MBR底部,靠近膜组件壁面的位置.另一部分污泥会随水一起绕过膜单元向上流动,膜单元底部污泥的浓度明显大于其它区域,也说明膜单元底部更容易发生沉积性膜污染,这也是将MBR中气-水-污泥三相体系简化为气-液两相流无法得到的结果.

图8 平面4处污泥体积浓度分布Fig.8 Distribution of sludge volume fraction at plane 4

2.2.2气泡大小的影响

曝气可以增强膜表面剪切力[23],对膜表面起到冲刷作用[24-25],有效地控制污泥浓度分布,进而改善膜污染.本研究中,与曝气相关的因素主要包括气泡大小、曝气面积及其位置.因此首先考虑气泡大小的影响,将气泡直径由前述的6 mm分别改为2、4和8 mm,进行模拟.

图9为不同气泡直径时MBR中的污泥分布,气泡直径的变化对污泥浓度分布几乎没有影响,仅在膜单元右下方污泥沉积有所增加.图10是膜单元底面的液体平均速度,随气泡直径的增大,水进入膜单元时的速度减小,速度的降低使污泥沉积的可能性增大,这与图9中直径增大后污泥在膜单元右下方沉积增多相互吻合,也说明气泡直径的增大一定程度上增加了膜污染的可能性.

图9 污泥在不同气泡直径下的分布Fig.9 Distribution of sludge with different bubble diameters

图10 不同气泡直径下膜单元底面水的平均流速Fig.10 Average water velocity at the bottom of the membrane unit with different bubble diameters

Wu等[10]将膜单元内液体的流速按照数值由小到大排列,认为排列在末尾5%内的区域为可能发生膜污染的区域,即“低速区域”,并将排列在5%处的速度值定义为“参考流速”,用于评估膜污染潜能.

参考流速越高,膜单元内可能发生污染的区域则越少.

参考此法,我们同样对膜单元内液体的平均流速进行分析.图11所示为不同气泡直径下膜单元内部液体平均流速的频率分布,膜单元内部液体平均流速在0～0.150 m/s之间,各种气泡直径下,0.030～0.075 m/s之间的速度所占比例均较大,在60%以上.当气泡直径为2和8 mm时,速度为0～0.015 m/s的低速区域频率均接近5.0%,气泡直径为4和6 mm时,低速区域的频率则低于2.5%.这说明,曝气时气泡直径为4～6 mm,出现膜污染的可能性比较低.

2.2.3曝气面积的影响

将曝气面宽度由最初的2 mm分别调整为1、3和4 mm.如图12所示,随着曝气面积的增加,曝气量逐渐增大,膜单元内的气含率明显增大,大量气体集中在膜单元内部和顶部,造成流场的混乱.图13所示为不同曝气面宽度下污泥体积浓度分布.由图13可知,曝气面积对污泥浓度影响不大,随面积增大,膜单元底部沉积的污泥略有增加.

图12 不同曝气面宽度下平面4内气含率分布图Fig.12 Volume fraction of gas in plane 4 with different aeration area

图13 不同曝气面宽度下平面4内污泥体积浓度分布图Fig.13 Volume fraction of sludge in plane 4 with different aeration area

图14为膜单元内液体流速范围的频率分布.曝气面宽度为1 mm时,大部分速度集中在0～0.045 m/s,0～0.015 m/s的低速区域占15%以上,大大增加了膜污染的可能.随着曝气面积的增加,低速区域的占比明显降低,曝气面宽度≥2 mm时,0～0.015 m/s的速度出现频率仅为2.5%左右.综合气含率分布等结果分析,当曝气面宽度为2 mm时,曝气效率较高,且出现膜污染的可能较低.

图14 不同曝气面宽度下膜单元内部液体流速频率分布Fig.14 Frequency distribution of liquid velocity in membrane unit with different aeration surface widths

2.2.4曝气面位置的影响

将曝气面由初始的距离膜单元底面350 mm分别上移100、200 mm和下移100、200 mm进行模拟.图15为5种不同曝气位置的液体速度矢量图,随着曝气面与膜单元距离的增加,膜单元上部速度明显增大,膜单元底部液体的流动方向亦发生改变,由图15(a)中流向膜单元,变为图15(b)～15(e)中绕过膜单元向上流动,且随着与膜单元距离的增大,这部分液体的速度和流量都增加.

图15 曝气面在不同位置时平面4的液体速度矢量图Fig.15 Velocity vectors of water in plane 4 with different aeration positions

图16的污泥浓度分布中,距膜单元较近时,污泥易被水夹带流向膜单元,膜单元底部出现较多的沉积污泥,污泥浓度较高.随着曝气面与膜单元距离的增大,部分液体绕过膜单元向上流动,部分污泥也会随之绕过膜单元流走,与膜单元距离越大,液体携带远离的污泥越多,沉积在膜单元底部的污泥则越少,膜污染的可能性减少.

图16 曝气面在不同位置时平面4的污泥体积分布图Fig.16 Sludge distribution in plane 4 with different aeration positions

另外,随着曝气面位置的改变,膜单元内部低速区域的相对频率也发生改变,如图17所示.在上移和下移后,膜单元内部速度在低速区0～0.015 m/s之间的频率均有所增加,特别是下移200 mm后,低速区域频率接近15%,较大程度增加了膜污染的可能性;而在初始位置时,低速区域频率仅为2.5%[见图11(c)].综合考虑膜单元内液体流速及底部的污泥沉积,在初始位置附近时,MBR内部污染的可能性最小.

图17 不同曝气面位置时膜单元内部液体流速频率分布Fig.17 Frequency distribution of liquid velocity in membrane unit at different aeration surface positions

结合气泡大小、曝气面积以及曝气面位置的研究结果,为了减少膜污染,应将曝气管安装在接近膜单元底面与MBR底面中间的位置,并控制曝气孔的大小,使气泡尺寸维持在4～6 mm之间.

2.2.5固相浓度的影响

将固相浓度由1%调整为0.5%和1.5%,分别进行模拟.污泥分布如图18所示,随着固相浓度的增加,膜单元底部堆积的污泥量增加,导致液体流动受阻,如图19所示,在浓度由0.5%增加到1.5%的过程中,流入膜单元的液体减少,膜单元上方的液体更分散,不利于液体在膜单元内部的循环.同时,随着固相浓度增加,污泥在膜单元底部堆积,膜单元内气体量减少,如图20所示,在固相浓度为0.5%及1.0%时,气体分布较均匀,当固相浓度再增至1.5%,膜单元内部气体分布不均,右侧甚至没有气体,也不利于膜单元内部的循环流动.

图18 不同固相浓度下污泥体积分布图Fig.18 Sludge volume distribution with different solid concentrations

图19 不同固相浓度下液体速度矢量图Fig.19 Velocity vector of liquid with different solid concentrations

图20 不同固相浓度下气含率分布Fig.20 Volume fraction of gas with different solid concentrations

图21所示为不同固相浓度下膜单元内液体流速频率分布,膜单元内部液体流速大致分布在0.02～0.08 m/s范围内,固相浓度为0.5%和1.0%时,0～0.015 m/s的低速区域的频率约为2.5%.固相浓度为1.5%时,低速区域频率超过5%,这与上文所分析的污泥在膜单元底部堆积阻碍液体和气体流动相吻合.

图21 不同固相浓度时膜单元内部液体流速频率分布Fig.21 Frequency distribution of liquid velocity in membrane unit with different solid concentrations

3 结论

1) 污泥相的引入,使三相模拟中液体的流动与气-液两相模拟中的液体流动明显不同:三相模拟中,部分液体绕过膜单元向上流走,而不是流入膜单元;膜单元内外区域液体速度差异明显,而两相流计算的内外速度差别较小.本文获得这些流动特性有利于对MBR流场的深入解析.

2) 曝气气泡直径增大,液体进入膜单元的流速减小,发生污泥沉积的可能性增大.本文研究范围内,当曝气面距离膜单元底面350 mm,曝气面宽度2 mm,气泡直径4～6 mm时,发生污泥沉积产生膜污染的风险最小.MBR设计中,应将曝气管安装在膜单元与MBR底面接近中间的位置,并通过调节曝气孔,控制气泡尺寸在4～6 mm之间,以达到减小膜污染的目的.

3) MBR内固相浓度不断增加,会导致污泥在膜单元底部的堆积,进而阻碍液体、气体围绕膜单元产生的循环流动,一定程度上加剧膜污染的风险.