基于逾渗理论的多孔过滤介质孔径分布估计方法

丁彬彬，李朝林

（哈尔滨工业大学深圳研究生院环境科学与工程研究中心，广东 深圳 518055）

引 言

孔径分布（pore size distribution,PSD）是多孔过滤介质重要性质之一，多孔材料过滤器的主要功能是截留流体中分散的固体颗粒，因此其孔径分布就决定了过滤精度和截留效率[1-2]，故其表征方法受到广泛关注[3]。

多孔材料的孔径是指多孔体中孔隙的等效直径，其测定方法也有很多，如断面直接观测法[4]、压汞法[5]、气体吸附法[6]、悬浮液过滤法[7]、X 射线小角度散射法[8]等。多孔材料的孔道形状复杂，影响孔径测量的因素也会很多，使用以上测定方法得出的结果也会有差异，故孔径的测定方法原理或过程最好与最终的使用情况相一致，对于过滤材料最好用过滤法。悬浮液过滤法测定多孔材料的孔径分布，其原理是对过滤前后悬浮液中粒子的粒度分布和浓度变化规律进行定量分析[9]，从而得出多孔体的孔径分布状况。

用于过滤法的数学模型主要有微管随机行走模型[10]、种群平衡理论[11-12]、各种Boltzmann 模 型[13-16]、平行管模型[7,17]。粒子的平均截留速率可用种群平衡理论[11,18]描述，它已成功应用于非完全堵塞过滤数据的处理[7,18]，其缺点是由实验数据估计的微管有效长度远大于平均孔径[7]。结果表明多孔介质中存在远大于孔径的特定相关长度。因此，需要新的随机理论来解释这些实验结果。

逾渗理论应用于多孔介质[19-21]总是和网格模型紧密联系在一起。网格模型主要关注孔尺度的物理性质，而逾渗理论则注重于随机性对宏观特性[22]的影响，如流体性质以及它们的相互作用等。网格模型结合逾渗理论也被用于过滤实验数据的处理[23]。

本文结合逾渗理论和网格模型，建立了过滤系数与孔径分布的联系方程，利用此方程结合过滤实验数据计算了过滤介质的孔径分布，并与其他方法获得PSD 参数进行比较和过滤过程的计算机模拟验证，探索该方法的可行性。

1 实验材料和方法

1.1 实验装置和材料

实验采用自制的深床过滤小试装置，装置图及流向如图1所示。其核心部分为深床过滤柱，其分解图如图2所示。整个柱子主要采用有机玻璃制成，除此之外，金属网为不锈钢材质，还有以橡胶做成的胶垫。每两个主要配件之间用胶垫相隔，并以玻璃胶粘住，以防止装置漏液，螺丝穿过所有配件将它们连接起来。深床过滤柱高5.8 cm，直径为4.0 cm，不锈钢金属网孔径为45 μm，填料球形玻璃珠粒径为 30～150 μm。过孔径为150 和58 μm 的筛后玻璃珠的粒径为58～150 μm（Medium58/150），过筛后的玻璃珠的粒径分布采用激光粒度分析仪Malvern Mastersizer 2000 来测量，检测结果见图3。羧基聚苯乙烯微球购自阿拉丁试剂有限公司，其表面覆盖有羧基基团，这些基团在碱性溶液（pH≈10） 中带负电，通过DLVO 理论[24]计算，粒子间相互作用力为斥力使得胶体溶液稳定而不会发生粒子聚集沉淀。研究中用的微球的尺寸见表1。

图1 深床过滤小试实验流程图及流向Fig.1 Scheme of experimental apparatus and flow direction

图2 深床过滤柱分解图Fig.2 Exploded diagram of column used for experiment

图3 Medium58/150 玻璃珠的粒径分布Fig.3 Grain size distribution of medium 58/150

表1 小试实验所用聚苯乙烯微球的尺寸及 对应的归一化出水浓度Table 1 Particle sizes and corresponding normalized effluent concentrations from experiment and reference

1.2 实验步骤

各种尺寸的单分散的羧基聚苯乙烯微球用NaOH 溶液配制成悬浮液，然后用超纯水稀释为 pH≈10 用于后续实验，微球悬浮液的浓度约为10 mg·L-1。调节泵速为10-5m·s-1以获得稳定的由下至上的液流，减小粒子由于密度差异造成的沉积和流动死区。低的粒子浓度(10 mg·L-1) 和短的注射时间可确保忽略由于粒子阻塞孔道引起的PSD的改变。对于出水的粒子浓度进行持续监测计数直至达到稳定状态。

1.3 其他PSD 估计方法

由于本文中的多孔介质是由球形玻璃珠堆积形成的，因此可以采用堆积法获得PSD，即通过计算机模拟符合图3分布的3D 球形颗粒的堆积可以得到每个球形颗粒的坐标，采用Delaunay 四面体将 堆积的小球进行分组，四面体的4 个顶点为小球的球心，通过求解与组成每个三角面的球形颗粒外切的最小外切圆来获得孔道分布函数PSD。

此外由于本文同时监测了进出水的粒子浓度，为便于与其他过滤模型比较，也采用平行管模型求解PSD，即通过向多孔介质中注入单分散的胶体粒子，通过监测穿透粒子的浓度并根据模型公式回归计算获得[7,17]。这种方法不需要介质颗粒的粒径分布。平行管模型公式[17]如下

式中，Ce(rs)和C0(rs)分别为粒径为rs的粒子进出水粒子浓度；rs为悬浮颗粒半径；L为填充柱的长度；l为两个孔室间的平均距离；fns为三角形孔室中粒子不能通过部分的流量比例[17]

式中，rp为孔道的半径；k1（rp）为孔道半径为rp的水力传导率，k1=πrp4/8；H(rp,x,t) 为随时间和空间变化的孔径分布密度函数。

2 理论分析

2.1 模型基本假设

图4 4×4 对角正方网格的流向示意图Fig.4 Schematic representation of 4×4 diagonal square lattice with periodical boundaries

本文中实验和网格模拟中粒子穿透和逾渗都只在一个方向发生，这种情况下，采用2D 和3D 的网格只是在网格规模和网格配位数上有差别，如果网格规模相同，网格配位数也相同，且采用PSD 参数也一致，那么粒子的移动模式基本相同。为了模拟计算程序简便，采用2D 网格表示多孔介质，而且也有较多的文献[2,25]表明2D 网格模拟的结果与3D 过滤实验数据也是相符的。采用配位数为4 对角 正方形网格，是由于堆积形成的孔道的配位数一般为4，文献[25]的研究中也表明配位数为4 时，实验和模拟结果基本一致。本文以2D 对角正方格来代表多孔介质，圆柱形孔道（键）半径rp由已知的孔道分布函数f(rp) 确定，键长度为常数。网格采用周期边界条件以避免表面效应。图4显示2D 周期正方格的各节点连接情况，图中的箭头表示流向。在这个模型中，假设粒子在网格中的行走是相互独立的，并且粒子捕捉主要机制是截留机理。由DLVO理论计算可知，高pH 和低盐度使得粒子介质表面均带有大量的负电荷，它们之间主要的相互作用力主要体现为斥力。在这种状况下，粒子与介质间的黏附作用可以忽略，截留成为粒子捕捉的主要机制。同时粒子在网格上移动遵循流量有偏[2]的原则。

2.2 统计参数

假设孔道的PSD 服从对数正态分布，其孔道半径rp的分布密度函数可以写为

在深床过滤中，由于本文研究的过滤机理只考虑截留，也就是说只有比悬浮颗粒粒径大的孔才可以让颗粒通过。那么，比较重要的一个参数就是“较大孔”（比悬浮颗粒大的孔）的比例，如式(4)

式中，rs为悬浮颗粒半径，m；rp为孔的半径，m；f(rp)为孔径分布函数。

假设单个孔道的流阻服从Poiseuille 定律，那么悬浮颗粒可以通过较大孔的平均概率为

本文中孔道的PSD 服从对数正态分布，代入其密度函数f(rp) 的表达式，式(5)可以化简如下

在逾渗阈值以下，颗粒通过介质的概率为0，在键逾渗的逾渗阈值处，f∗l的值为

式中，f∗c是基于流量有偏计算的逾渗阈值；rsc是可通过粒子半径的阈值。

2.3 基于逾渗理论的过滤过程分析

设网格中大于rs的键为活性键，当粒子的半径大于网格中所有的键时，所有的键都是非活性键，粒子不能通过网格介质，随着粒子半径的减小，网格中的活性键增加并相互连接形成集团，根据逾渗理论，当rs达到逾渗阈值时，由活性键组成的集团即粒子可行通路将贯穿整个网络，粒子将通过网格介质。从图5中可知，在此模型中，只有行走在集团主干上的粒子才能通过网格介质，而进入集团分支的粒子将被网格介质截留。

基于上述分析，Yuan 等[23]提出了两条幂律描述过滤系数与网格小孔流量比间的关系，这里简要推导粒子尺寸在远离逾渗阈值时（rs≪rsc）的幂律关系，rs在逾渗阈值附近时，粒子出水浓度非常低，不好检测，没法验证，这里不做分析。设悬浮液带着悬浮颗粒流入无限集团主干的概率为B，那么1-B就表示悬浮液流入分支的概率。根据逾渗理论[26-27]，假设悬浮液流入无限集团的比例为，则B可以表示为

图5 过滤网格的入口、出口、集团主干和分支示意图Fig.5 Schematic representation of a network backbone,branch and dead end with line inlet and outlet nodes

由于颗粒只有处于主干才能通过整个介质，那么流入分支的颗粒只能被捕捉，即1-B近似表示颗粒被捕捉的概率。在rs≪rsc时，集团主干占的比例远大于分支，悬浮液流入无限集团的比例为1-B可以表示为

集团分支的密度随着粒子尺寸的减小急剧减少，以lw表示网格上两分支的平均距离，根据逾渗理论lw可能符合如下关系[23]

而过滤系数λ表示的是单位长度多孔过滤介质捕捉概率

定义指数β=1-γ，K为比例常数，则有

3 结果与讨论

3.1 实验结果和PSD 参数的获得

实验测得的归一化出水粒子浓度如表1所示，随着悬浮粒子半径的增大，出水浓度急剧降低。根据经典过滤理论[28]过滤系数的计算公式如下

式中，λ为深床过滤实验对应的过滤系数；L0为填充柱子的长度。由逾渗理论可知过滤系数λ满足式(12)，由此可推出

考虑到n次实验用的悬浮粒子半径为rs1，rs2，···，rsn，每次实验的归一化出水浓度Ce(rsi)/C0(rsi)由实验测得，设多孔过滤介质孔径分布满足对数正态分布，每次实验数据满足如下方程

式(15)中K、μ、σ、β为未知参数，参数的获得可通过代入实验数据利用非线性最小二乘法求解如下的优化问题

求解的结果见图6和表2，如果以μm 为单位，求解获得PSD 的统计参数σ=0.535，μ=1.799，由此可知平均孔径为6.98 μm，孔径方差为 4.01 μm，指数β=0.80，比例系数K为-15.846（L0的单位为mm），为了对上面的算法进行评估并和其他方法 比较，接着采用其他方法估算了过滤介质的PSD 参数。

图6 归一化出水粒子浓度与粒子半径的关系Fig.6 Predicted and experimental normalized breakthrough concentration versus particles size

表2 3 种方法估计的堆积玻璃球（Medium58/150）形成多孔介质的孔径分布参数Table 2 Pore size distributions evaluated from grain sizedistributions and experiment with three methods

3.2 其他方法获得PSD 参数

采用堆积方法和平行管模型求解获得的PSD参数见表2，由表可知采用平行管模型求解获得PSD 参数和孔径与逾渗模型求得结果非常接近，平行管模型求解获得的PSD 统计参数σ=0.580，μ=1.780，同时代入式(1)求解归一化出水浓度Ce/C0的预测值，如图6所示，由逾渗模型和平行管模型计算得到Ce/C0也与实验值一致。采用堆积法获得的PSD 参数与前面两种方法计算的结果有比较大的差别，平均孔径偏大为9.10 μm，孔径方差偏小为3.30 μm，这可能是由于采用堆积体系偏小的原因，由于计算机计算能力的限制，用于堆积法的粒子数目只有5 万个。

3.3 其他文献数据的验证

为了检验模型的适用性，本文还对文献[17,23]中实验数据做了计算验证，见表1中Medium 30/125的数据，计算结果见表3，其结果与本文计算结果类似，平行管模型求解获得PSD 参数和孔径与逾渗模型求得结果非常接近，文献中PSD 参数采用的是基于笛卡儿理论的Monte Carlo 拉丁超立方抽样模拟[7]，该方法假设介质颗粒是紧密接触密堆积，与实际的随机自然堆积不相符，因此计算的孔径方差明显偏小，而且文献[23]中采用此PSD 用于模拟得到结果与实验结果并不相符。

表3 3 种方法估计的堆积玻璃球（Medium30/125）形成多孔介质的孔径分布参数Table 3 Pore size distributions evaluated from grain size distributions and experiment with three methods

3.4 采用不同PSD 参数进行过滤模拟的结果

为了检验获得PSD 的参数准确性，将获得的PSD 参数用于过滤过程模拟，即为网格上每个键随机设置半径rp，其分布函数f(rp)由PSD 参数确定，模型采用网格随机行走模型，为获得网格模型的出水粒子浓度，在每个节点需要运用质量守恒计算每条键的通量。连接节点i和相邻节点j的键的水力传导系数gij，可由Poiseuille 公式计算

式中，η为液体的动态黏度；l为微管长度；为微管半径。

质量守恒要求每个节点相连的微管流入和流出的流量相等，净流量和为0，公式描述如下

图7 不同PSD 参数下的模拟的Ce/C0 与粒子尺寸的关系Fig.7 Simulated and experimental normalized effluent concentration versus particles size with applied PSD parameter estimated by different methods

式中，∆pij为节点i和j之间的压力差。在网格的入口和出口出的压力设为常数，网格采用周期边界条件以避免表面效应。求解由上述方程组成的线性系统，可获得各个节点的压力值，进而利用Darcy 定律获得各个微管的流量和流场的细节。

由于模拟的条件为稀溶液和短时间注射，故由孔阻塞而引起的渗透率的降低可以忽略，粒子在网格上的行走是独立的，粒子是逐个地注入到100×100 对角正方格中，每次网格上只有一个粒子在移动。颗粒依据流量有偏[2]的原则选择即将行走的孔道。粒子行走的独立性保证在模拟中可以采用并行机制，节省模拟时间，充分利用计算资源。本文模拟采用两颗六核Xeon E5-2620 处理器，主频2.0 GHz。每次实验注入粒子数目为105确保统计分析的精度。

过滤模拟的结果见图7，由图可知采用平行管模型和逾渗模型求得求解获得PSD 参数用于模拟得到Ce/C0与实验值基本一致，说明本文中提出的基于逾渗理论的过滤模型用于PSD 的估计是比较准确和有效的，而采用基于笛卡尔理论的Monte Carlo 拉丁超立方抽样法和堆积方法获得PSD 用于过滤模拟，其结果与实验并不相符。

逾渗模型和平行管模型求解多孔介质PSD 参数都需要知道测试粒子的尺寸以及粒子的进出水浓度，不需要知道介质颗粒的粒径分布，对于非堆积形成的多孔介质也是适用的，如果没有合适粒子测试多孔介质，或者进出水粒子浓度变化不显著就不能采用逾渗模型和平行管模型。堆积方法是通过计算机模拟实际堆积，通过几何的方法求解PSD，需要知道介质颗粒的粒径分布，只适合球形或近似球形颗粒堆积形成的多孔介质，而且需要相当大规模的模拟才能得到较准确的结果。对于形状不规则介质颗粒堆积的多孔介质或者不是堆积形成的多孔介质并不适用。

4 结 论

本文结合逾渗理论和网格模型，建立了过滤系数与PSD 的联系方程，利用此方程结合过滤实验数据计算了过滤介质的孔径分布，并与其他方法获得PSD 参数进行比较和过滤过程的计算机模拟验证，逾渗模型估算的PSD 用于过滤模拟的结果和实验数据符合较好，表明基于逾渗理论的过滤模型用于PSD 参数的估计是比较准确和有效的。

[1]Xu S,Gao B,Saiers J E.Straining of colloidal particles in saturated porous media [J].Water Resources Research,2006,42 (12):W12S-W16S

[2]Rege S D,Fogler H S.Network model for straining dominated particle entrapment in porous media [J].Chemical Engineering Science,1987,42 (7):1553-1564

[3]Gibson L J,Ashby M F.Cellular Solids:Structure and Properties [M].2nd ed.Cambridge:Cambridge University Press,1999:532

[4]Langlois S,Coeuret F.Flow-through and flow-by porous electrodes of nickel foam (Ⅰ):Material characterization [J].Journal of Applied Electrochemistry,1989,19 (1):43-50

[5]Hoefner M L,Fogler H S.Pore evolution and channel formation during flow and reaction in porous media [J].AIChE Journal,1988,34 (1):45-54

[6]Liu Peisheng (刘培生).Determining methods for aperture and aperture distribution of porous materials [J].Titanium IndustryProgress(钛工业进展),2006,23 (2):29-34

[7]Chalk P,Gooding N,Hutten S,You Z,Bedrikovetsky P.Pore size distribution from challenge coreflood testing by colloidal flow [J].Chemical Engineering Research and Design,2012,90 (1):63-77

[8]Banhart J.Manufacture,characterisation and application of cellular metals and metal foams [J].Progress in Materials Science,2001,46 (6):559-632

[9]Fan Yunge (范云鸽),Li Yanhong (李燕鸿),Ma Jianbiao (马建标).Characterization of the nanoscaled pores in porous polydivinylbenzen adsorbents [J].Acta Polymerica Sinica(高分子学报),2002 (2):173-179

[10]Yuan H,Shapiro A A.A mathematical model for non-monotonic deposition profiles in deep bed filtration systems [J].Chemical Engineering Journal,2011,166:105-115

[11]Santos A,Bedrikovetsky P.A stochastic model for particulate suspension flow in porous media [J].Transport in Porous Media,2006,62 (1):23-53

[12]Sharma M M,Yortsos Y C.A network model for deep bed filtration processes [J].AIChE Journal,1987,33 (10):1644-1653

[13]Tang G H,Ye P X,Tao W Q.Pressure-driven and electroosmotic non-Newtonian flows through microporous mediavialattice Boltzmann method [J].Journal of Non-Newtonian Fluid Mechanics,2010,165 (21/22):1536-1542

[14]Shapiro A A,Wesselingh J A.Gas transport in tight porous media:gas kinetic approach [J].Chemical Engineering Journal,2008,142 (1):14-22

[15]Zhu Weibing (朱卫兵),Wang Meng (王猛),Chen Hong (陈宏),Han Ding (韩丁),Liu Jianwen (刘建文).Lattice Boltzmann simulations for fluid flow through porous media [J].CIESC Journal(化工学报),2013,64 (S1):33

[16]Sun Meiyu (孙梅玉),Ji Zhongli (姬忠礼).Parallel computing of gas flow through ceramic filter by using lattice Boltzmann method [J].CIESC Journal(化工学报),2010,61 (6):1423

[17]You Z,Badalyan A,Bedrikovetsky P.Size-exclusion colloidal transport in porous media—stochastic modeling and experimental study [J].SPE Journal,2013 (8):620-633

[18]Shapiro A A,Bedrikovetsky P G,Santos A,Medvedev O O.A stochastic model for filtration of particulate suspensions with incomplete pore plugging [J].Transport in Porous Media,2007,67 (1):135-164

[19]Perrier E M A,Bird N R A,Rieutord T B.Percolation properties of 3-D multiscale pore networks:how connectivity controls soil filtration processes [J].Biogeosciences,2010,7 (10):3177-3186

[20]Berkowitz B,Ewing R P.Percolation theory and network modeling applications in soil physics [J].Surveys in Geophysics,1998,19 (1):23-72

[21]Selyakov V I,Kadet V.Percolation Models for Transport in Porous Media:with Applications to Reservoir Engineering [M].Springer,1997

[22]Berkowitz B,Balberg I.Percolation approach to the problem of hydraulic conductivity in porous media [J].Transport in Porous Media,1992,9 (3):275-286

[23]Yuan H,Shapiro A,You Z,Badalyan A.Estimating filtration coefficients for straining from percolation and random walk theories [J].Chemical Engineering Journal,2012,210:63-73

[24]Jia X,Gregory J,Williams R.Particle Deposition & Aggregation:Measurement,Modelling and Simulation [M].Elsevier Science,1998

[25]Ding B,Li C,Zhang M,Ji F,Dong X.Effects of pore size distribution and coordination number on the prediction of filtration coefficients for straining from percolation theory [J].Chemical Engineering Science,2015,127:40-51

[26]Reyes S,Jensen K F.Estimation of effective transport coefficients in porous solids based on percolation concepts [J].Chemical Engineering Science,1985,40 (9):1723-1734

[27]Larson R G,Davis H T.Conducting backbone in percolating Bethe lattices [J].Journal of Physics C:Solid State Physics,1982,15 (11):2327-2331

[28]Santos A,Barros I P H.Multiple particle retention mechanisms during filtration in porous media [J].Environmental Science & Technology,2010,44 (7):2515-2521