在过去几十年里, 膜过滤技术已逐渐成为水处理领域一项重要工艺。膜过滤工艺是一种压力驱动过程, 以膜两侧的压力差为驱动力, 以膜为过滤介质, 在一定的压力下, 实现了对原液的分离和浓缩的目的。膜工艺除了在工业生产中得到大量应用, 现已成为传统水处理工艺的一种附加或替代工艺^[1,2]。

膜污染是水体中杂质、杂质与膜表面一系列复杂的物理、化学和生物作用的结果^[3]。膜污染会导致过滤阻力 (TMP) 增大和通量下降, 化学清洗的频率和能耗增加, 进而缩短膜寿命及提高制水成本, 已阻碍膜分离技术的广泛应用^[4~6]。因此如何减少膜污染是膜技术进一步发展关键所在。

然而对膜污染的研究经常会受到实验条件和过滤装置的限制, 只能局限于定性或是半定量研究, 无法对膜污染这一问题深入到机理层面进行研究分析^[7,8]。计算流体力学 (CFD) 是通过计算机数值计算和图像显示, 对流体流动和热传导等物理现象做出的分析的技术, 所以能够很好克服单纯试验分析的缺点^[9]。

在早期的膜过滤CFD仿真研究中, 将膜表面视为不可渗透壁面或是多孔介质层, 对膜内流场进行研究分析。后来Pinho等^[10]将传质模型耦合到膜CFD模型中, 对膜边界进行定义, 得到膜通道内物质浓度分布和截留率, 并与试验结果相吻合;Ahmad等^[11], 采用薄膜理论, 利用CFD软件Fluent v6用户自定义函数 (UDF) 描述浓差极化现象, 成功模拟了狭通道内的层流状态下膜表面浓差极化现象, 预测出膜内浓差极化层厚度随雷诺数与跨膜压差的变化规律。但其仅是通过定义物质浓度分布关系实现了膜的“可渗透性”, 在模拟结果上具有一定局限性;Rajabzadeh等^[12]在研究电解酸化大豆提取液在错流过滤中的可逆与不可逆污染行为时, 利用达西定律, 耦合串联阻力模型 (resistance-in-series) 对过滤过程中的通量模型进行了更为细致的定义, 成功模拟出膜通量、阻力系数与边界物质浓度随时间变化趋势并与试验结果相吻合;宋卫臣等^[13,14]利用有限元软件Comsol Multiphsics耦合了多孔介质流场、自由开放流场, 并基于死端超滤建立起来的污染模型以及文献中的数据, 对吸附速率系数和解吸附速率系数进行了数值拟合, 建立了描述完整超滤错流污染过程的数学模型。

目前, 膜过滤模型的研究依然停留在对膜内流场状态的研究, 而对膜过滤过程中的膜污染动态变化过程研究分析鲜见报道。本文基于Marcos等^[15]建立的回流错流过滤污染模型的基础上, 通过Comsol Multiphsics有限元软件自由流动、Darcy定律以及稀物质传递模块结合优化处理的膜阻力公式, 模拟了无回流错流过滤的膜通量变化趋势, 并通过试验加以验证, 从而得到中空纤维超滤膜膜内污染以及浓差极化层厚度变化。

试验过滤过程如图1所示, 膜组件采用错流过滤, 恒压运行方式, 操作压力为0.10MPa, 在原料箱中使用自配2mg/L腐殖酸原水进行过滤试验。试验装置为某公司提供的多功能膜试系统, 如图2所示。过滤试验用超滤膜组件特征参数如表1所示。

用式 (1) 与式 (2) 、式 (3) 来描述流体在膜内的质量守恒与自由流动, 用式 (4) 来描述污染物质 (污染物质) 的分布情况。在确定控制方程的同时需要对模拟的流体进行以下假设:流体是不可压缩的、稳态的、等温的;流体具有恒定的密度与粘度系数^[16];溶剂的扩散系数为常数^[17];通道入口速度为完全发展流体;膜表面无滑移^[18];局部通量由阻力串联模型决定^[19]。

连续性方程:

N-S方程:

对流扩散方程:

式中v_X———沿X方向速度, m/s;

u_Y———沿Y方向速度, m/s;

ρ———溶液密度, kg/m³;

P———压强, Pa;

μ———动力粘度系数, Pa·s;

D———扩散系数, m²/s;

C———溶质浓度, mg/L。

模型通过耦合Comsol Multiphysics软件中的自由流动、Darcy定律以及稀物质传递模块, 模型如图3所示, 且边界设定如下:

入口:膜通道入口处假定流动是完全发展的, 且为抛物线流。入口处物质浓度即为原料液浓度。

出口:出口压强等于大气压。

膜壁面:u_Y=v_p

膜壁面假设为无滑移壁面, x方向的速度v_X=0, y方向速度u_Y大小即为渗透通量v_p大小。在进行过滤试验中, 膜对配水浊度去除效果很好, 因此假设在过滤过程中膜表面对污染物质完全截留, 膜表面物质浓度变化如式 (5) 所示:

膜表面处的y方向速度u_Y即渗透通量j_v采用达西定律式 (6) 来定义:

式中μ₀———原液粘度系数;

R_T———总膜阻系数;

———坐标为x位置的膜表面t时刻的压强;

P_permeate———膜渗透侧的压强, 在模型中认为等于大气压强。

在过滤过程中, 污染物质在膜表面积累形成滤饼层或者吸附在膜表面以及浓差极化现象出现, 都会导致总阻力系数R_T增大, 所以将R_T以下几个部分:

式中R_cake———滤饼层产生的阻力系数;

R_a———物质在膜表面或膜孔之间吸附产生阻力系数;

R_m———膜纯水过滤时的阻力系数。

Song等^[20]研究表明, 错流过滤过程中通量会有3个变化阶段:初始快速显著下降期、逐渐下降期和后续期。其中快速显著下降期与物质吸附产生的阻力系数R_a相关, 而逐渐下降期与滤饼层阻力系数R_cake有关;宋航等^[21]在对吸附阻力研究分析得到:在过滤初期的吸附阻力系数R_a在过滤初期上升会快, 然后会趋于平缓, 最终会趋近于或维持在一个稳定阻力系数值, 如式 (9) 所示:

式中a, b———吸附阻力系数R_a的模型参数, a的取值影响吸附阻力系数趋于稳定所用时间, b的取值即为趋于稳定时的吸附阻力系数的倒数。

根据Tansel等^[22]人的研究:膜阻力系数分为随时间变化阻力系数和恒定阻力系数, 而随时间变化的阻力系数随时间的增加量正比于已经产生的阻力系数。在本文中, 由于吸附阻力系数R_a在短时间内即可达到饱和值, 所以只将滤饼阻力系数R_cake定义为随时间变化阻力系数。因此滤饼阻力增加量方程可用式 (10) 表示

式中k———比例因子;

R_cake———滤饼阻力系数。

由于浓差极化使边界层上的溶质达到饱和浓度时, 溶质在膜表面析出, 从而形成滤饼层并覆盖膜表面, 增加渗透阻力^[15,23]。本文定义浓差极化所产生阻力系数R_pol作为R_cake初始值, 则滤饼层阻力系数随时间变化如式 (11) 、式 (12) 所示

模型所需要的参数中膜丝自身阻力系数R_m由纯水试验确定。结合文献[21, 24]对吸附阻力R_a测定的方法, 通过水力清洗过滤原料液后的膜管, 去除膜表面所形成的滤饼层后, 在同样操作压力下, 利用渗透液进行过滤, 最后得到吸附在膜上的阻力系数R_a。因此, 依次测定出膜管过滤不同时长原料液所产生的吸附阻力系数R_a值, 建立1/R_a~1/t的线性回归方程即:, 得到参数a、b, 如图4所示。

参数:a=5.0×10^-12, b=0.74×10^-12, 相关系数R²=0.98。

由式 (8) 、式 (9) 、式 (12) 可得:

式中R_T、R_m以及R_cake0可由试验得到, R_a由上述拟合方程得到, 建立ln (ΔR_cake) ~t的线性回归方程, 如图5所示。

由由回归方程取ln (A) 、斜率k, 得到参数取值A=0.187×10¹², k=37.0。

表2为模型最终参数设置。

计算机采用自配高性能电脑, 操作系统为Win 7, CPU采用intel酷睿i7-4790K, 运行内存32G, 已确保模拟高效进行。为了更好将N-S描述的自由流场、溶质对流扩散场以及膜阻力模型耦合, 求解软件采用具有强大多场耦合功能的有限元软件COMSOL Multiphysics 4.4, 膜内为自由流动与稀物质传递模块的耦合, 膜表面采用Darcy定律模块进行设定, 从而避免了编写复杂的计算机语言。模型为二维旋转对称模型, 采用均匀分布的网格, 网格形状为矩形。在轴长方向上设置单元数为1 000, 垂直于膜表面方向上设置单元数20。网格总数为20 000。

图6为TMP=0.10 MPa时, 模型预测通量以及试验通量随时间变化趋势图。由图可知, 按图将过滤过程可分为3个阶段:显著下降期、逐渐下降期和后续下降期, 且模拟结果与试验结果较吻合。与Song等^[20]的研究结果不同的是, 在进行验证试验中未呈现出其所提到的后续稳定期, 这可能是由于实验滤饼层未能完全覆盖膜表面。图6中每个时间点上的通量J, 都为膜丝局部通量的平均值, 通过根据模型局部渗透速率v, 按照等式而来:

图7为操作压力为0.10 MPa时, 在0min、60min和120min沿膜长局部渗透速率v分布图。由图可知, 渗透速率v沿着x轴方向线性下降, 这是由于膜内压强是沿x轴方向而减小的, 所以局部渗透速率也随着x轴方向减小, 在TMP为0.10 MPa下, 0min的直线斜率绝对值要大于60 min和120min, 这是由于随着过滤进行, 膜丝入口段和出口段受污染程度将逐渐趋近, 导致渗透速率沿X轴变化减小。同时, 随着过滤进行, 膜沿X轴相同位置的膜渗透速率u_Y都在逐渐减小。

膜污染会引起阻力系数增大导致过滤过程中膜通量下降, TMP为0.10 MPa时阻力系数变化, 如图8所示。

由图8可知, 试验测得的总阻力系数与模拟所得的R_T吻合。与此同时, 总阻力系数由上升期、平缓期和后续上升期3个阶段构成。过滤初期 (0~20min) , R_T随着物质吸附产生的阻力系数R_a增大而逐渐增大, 导致通量下降;平缓期 (20~60 min) , 由于物质吸附达到饱和, 此时滤饼层尚未形成, 污染物质在膜表面吸附和扩散达到平衡, 即R_T和膜通量均无明显变化;过滤后期 (60~120min) , 滤饼层形成, 导致其阻力系数R_cake迅速增大, 从而导致R_T快速增大和通量急剧下降。

浓差极化现象的产生是由污染物质与膜发生一系列作用产生的, 图9与图10分别在TMP=0.10MPa时污染物质沿x轴方向 (膜壁面) 和y轴方向 (选取出口x=200 mm) 的浓度分布图。由图9可知, 污染物质浓度C沿着x轴方向不断上升, 在出口处达到最大。这是由于出口溶液与膜面的切应力较小, 使得污染物质容易积累, 从而使其浓度C达到最大值。同时, 随着过滤的进行, 相同位置的污染物质浓度C也在逐渐增大。值得注意的是, 120min时污染物质沿x轴分布比60min更均匀, 由此证明了随着过滤的进行入口段和出口段受污染程度逐渐接近。图10反映的是出口处污染物质浓度C沿y轴方向的分布规律, 随着过滤的持续进行, 壁面污染物质浓度C也逐渐增大, 从0min的2 mg/L增加到120min的10mg/L, 且对于同一个时刻, 越靠近膜壁面, 污染物质浓度C也越高。根据文献^[25]对浓差极化层厚度的定义, 当膜内某处物质浓度时, 此区域正好发生浓差极化现象。模型中物质初始浓度C₀=2.0mg/L, 因此根据上述公式可得C>2.002mg/L即在膜出口截面上, 物质浓度大于2.002mg/L的区域即为浓差极化区域, 如图10所示。过滤从0 min进行到60 min, 浓差极化区域宽度从0增加到了0.5 mm, 60 min进行到120min时, 浓差极化区域仅增加了0.14 mm。由此说明, 在过滤前期和中期 (0~60min) , 污染物质在膜表面的累积是导致膜通量下降的主要因素;过滤进行到60~120 min, 快速形成的滤饼层导致膜阻力系数增加, 从而成为膜通量变化主要因素。

(1) 通过耦合自由流场、浓度场, 引入通量模型, 建立了中空纤维膜错流过滤二维CFD瞬态模型, 模型预测数据很好的拟合了试验数据。

(2) 通过定义总膜阻系数R_T=R_cake+R_a+R_m, 建立了瞬态的膜阻力串联模型。试验与模拟结果表明, 膜污染分为3个阶段, 首先是由物质吸附R_a引起R_T增大, 使得通量下降;然后由于污染物质在膜表面的吸附饱和, 滤饼层污染尚未形成时, R_T基本不变, 因此, 通量保持不变;最后由于滤饼层阻力系数R_cake快速增加, 导致R_T增大, 从而通量明显下降。

(3) 渗透速率u_Y沿着x轴方向线性下降, 这是由于局部渗透压强也随着x方向减小, 并且在同样的TMP下, 0min的直线斜率绝对值要大于60min和120min, 这是由于随着过滤进行, 膜丝入口段和出口段受污染程度将逐渐趋近, 导致渗透速率沿x变化减小。

(4) 污染物质浓度沿着x轴方向不断上升, 在出口处达到最大, 沿y轴方向, 越靠近膜壁面浓度越高;同一位置, 随着过滤进行, 污染物质浓度不断升高。通过污染物质浓度随时间变化图可以发现, 膜污染前期造成通量变化, 主要是物质吸附产生的阻力造成的, 后期通量变化主要是由于滤饼层的形成。