颗粒床内固液过滤的三维CFD-DEM模拟

肖 桐， 王千红， 沈盈莺， 惠 卿， 陈 晗， 程振民

（华东理工大学化工学院，化学工程联合国家重点实验室，上海 200237）

颗粒床过滤是一种典型的深床过滤，与表面过滤方式不同，它主要将悬浮物拦截在过滤介质的内部。颗粒床过滤器结构简单、成本低廉、运行费用少、过滤效果好，在石油化工、生物制药、污水处理、农田灌溉等领域应用广泛[1-4]。为深入了解深床过滤机理和规律，更好地指导过滤器的工艺设计和操作运行，研究者们构建了不同的经验和理论模型，包括唯象模型[5]、迹线理论模型[6]、随机模型[7]、逾渗理论模型[8]和网络模型[9]等。但是这些经验或理论模型有许多缺陷：唯象模型是一种经验模型，有些参数并没有实际的物理意义；迹线理论模型不能反映过滤过程中的水头损失影响；随机模型不能预测滤出液的颗粒浓度；逾渗理论模型不适合非截留机理主导过滤的情况；网络模型假设的通道和实际的过滤床层孔道差别较大。因为这些数学模型的局限性，所以需要探索更为先进的研究方法。

随着计算机技术的发展和离散单元法的提出，基于计算流体力学（CFD）和离散单元法（DEM）的数值模拟研究方法被越来越多的学者所采用。Zhao等[10]用耦合CFD-DEM的方法模拟与采矿和岩土工程相关应用中的流体和颗粒相互作用的行为。方乙琳等[11]运用CFD-DEM耦合的手段，结合标准k-ε模型（k为湍流动能，ε为湍流耗散率）、Eulerian多相流模型，模拟沸腾床外排废催化剂的分选，得到催化剂分离状态及流场变化，探究分选柱中催化剂的分选机理和颗粒运动特性。谢晓翠等[12]利用CFD软件中的流体体积函数(VOF)模型捕捉气液相界面流动，利用离散相模型（DPM）考察纳米颗粒与水相的相互作用，进而模拟纳米颗粒、水和空气的流动。Tao等[13]采用三维多时间尺度CFD-DEM耦合方法研究不同类型纤维阵列对微米级颗粒的过滤，并比较了以平行和交错方式均匀排列的纤维的过滤性能。Schilling等[14]用Euler-Lagrange方法和CFD-DEM耦合数值模拟滤料纤维上颗粒的迁移和沉积行为，并分析了过滤纤维上由沉积颗粒构成的树枝状结构的重新分配。Qian等[15]在研究纤维层过滤气溶胶时，构建出基于随机算法的三维纤维床层，并将CFDDEM耦合模型应用到纤维床层过滤中，探究不同空隙率的纤维床层中杂质颗粒的沉积形态，以及过滤效率和床层压降随表观过滤速率、床层空隙率、杂质颗粒大小的变化。Wang等[16]基于CFD模拟，对两种典型填料颗粒过滤器的颗粒过滤特性进行了数值模拟，得到滤料颗粒表面上的流体流动、压降、颗粒沉积特性以及颗粒沿过滤床的沉积分布。

本文采用CFD-DEM耦合方法，对高20 mm、颗粒直径为1.0 mm的三维随机堆积颗粒床进行固液过滤模拟研究。这种双平台耦合计算的方法，可以避免刻画颗粒内部的复杂流道，大大减少前期准备的工作量，为颗粒床过滤的模拟研究提供新思路。

1 模型描述和数值计算方法

1.1 数值计算模型

1.1.1 流体相的控制方程 在Eulerian模型中，考虑到杂质颗粒沉积对滤料床层空隙率的影响，在Navier-Stokes方程中引入床层空隙率 ( φ) ，从而实现流体对颗粒及颗粒对流体的双向耦合分析。流体相的连续性方程和动量方程分别表述如下：

其中：μ是流体黏度，Pa·s； u 是流体速度，m/s； ρf是流体密度，kg/m3； g 是重力加速度，m2/s；S是作用在网格单元体积的流体阻力 FD的总和，N； Δ V 是网格单元体积，m3。

对于填料床内流体的流动状态的判断，景有海[17]依据前人的研究，提出Reynolds数计算公式（式（4）），只有当 Re 小于2时，流体流动状态才是层流。

其中：dp为颗粒直径； v 是滤层中水流平均速度； α 是滤料颗粒的面积形状系数，取值为1。据此计算的层流临界速率为0.006 8 m/s，而本实验中的流速均大于0.13 m/s，所以，应该采用湍流模型计算。

应用SST k-ω湍流模型计算颗粒床层内部的湍流流动。SST k-ω模型从k-ε和k-ω模型发展而来，可以在广泛的领域中独立于k-ε模型，使得在近壁自由流中k-ω模型有广泛的应用范围和精度。模型中的湍流动能k和比耗散率ω根据式（5）和式（6）计算：

其中： ui是流体速度的分量；Gk表示湍流的动能；Gω是由 ω 方程产生的项；Γk与 Γω分别代表k与 ω 的有效扩散项；Yk与Yω分别代表k与 ω 的发散项；Dω代表正交发散项；Sk与 Sω为自定义项。其他详细的SST k-ω湍流计算模型可以参见文献[18-19]。

1.1.2 固体相的离散模型 颗粒床过滤中包含两种离散相：一是滤料床层的颗粒固相，二是流体中杂质颗粒相。悬浮液中的杂质颗粒在流体中运动，单一颗粒所受外力来自颗粒与颗粒之间作用力以及流体对颗粒的作用力。运用Lagrange方法研究杂质颗粒在流体中的运动轨迹。由于模型中的杂质颗粒足够大，所以，本文只考虑流体对颗粒的曳力、颗粒自身的重力以及颗粒与颗粒的碰撞等因素，忽略范德华力、双电层作用力和布朗运动等因素对颗粒的影响。颗粒在流场中的运动方程为：

其中：mp是颗粒质量，kg； up是颗粒速度，m/s； Fg是颗粒的重力，N；fD(u-up) 是流体对颗粒的曳力，N。

fD可通过式（8）求得：

其 中非线性曳力系数 (CD) 为：

颗粒与颗粒之间的接触作用采用Hertz-Mindlin接触模型[20]进行计算，这种接触作用可用非线性的软球模型表示。颗粒与颗粒之间的作用力和力矩表示如下：

其中： IP是颗粒的转动惯量，kg·m2； wp是颗粒的转动速度，rad/s； Mi是颗粒碰撞扭矩，N·m。详细的颗粒与颗粒的作用力和力矩计算参见文献[21]。

1.1.3 耦合计算方法 CFD与DEM的耦合采用Eulerian耦合法，除考虑液相和固相之间的动量交换外，还考虑固体颗粒对流体相的影响。CFD平台采用Fluent软件，DEM平台采用EDEM软件，两个平台之间的耦合程序采用自编UDF（User-Defined Function，用户自定义函数）。如图1所示，首先在Fluent中进行液相流场初始化，在一次时间步长内，流场计算至收敛；然后在EDEM中确定颗粒位置。考虑到流体对颗粒的曳力、颗粒自身以及颗粒之间的作用等，计算颗粒碰撞方程，更新颗粒位置信息；将颗粒对流体的作用力、颗粒位置等信息，传递给Fluent进行下一时间步长内液相流场的计算。

图 1 CFD-DEM耦合计算流程Fig. 1 Flow chart of the CFD-DEM coupling

1.2 三维过滤模型

图2所示为随机堆积的颗粒床层。通过EDEM中圆柱形“颗粒工厂”快速生成随机位置的颗粒，在重力作用下于过滤床内堆积，进而形成随机堆积的三维过滤模型。图3所示为数值计算模型中的颗粒过滤床，直径（R）为10 mm，高（H）为20 mm，其中球形颗粒的直径（D）为1.0 mm，床层空隙率（φ）为0.43。额外的入口和出口延伸区域用以维持颗粒均匀的入口速度并防止其从出口回流。流体介质为常温水，密度约为1 000 kg/m3，黏度约为0.001 Pa·s。具体几何参数和表观过滤速率（uA）等模拟条件见表1，颗粒相关参数的设置见表2。网格单元采用结构性网格划分，网格数为2 261，最小体积为6.5×10-10m3，最大体积为3.7×10-9m3。Fluent中的时间步长设置为2×10-3s，为了准确捕获颗粒接触行为，EDEM的时间步长比Fluent的时间步长通常要小很多，这里两者比值设置为1∶100。所有的模拟计算均在Intel®Xeon®E5645 2.4 GHz（24核，48 GB内存）服务器上进行，计算的平均时长约为48 h。

图 2 随机堆积的颗粒床层Fig. 2 Randomly packed granular layer

2 数值模型的实验验证

在颗粒床过滤实验中，过滤效率 ( η ) 定义为：

其中： c1为母液质量浓度，mg/L； c2为滤出液质量浓度，mg/L。假设每个悬浮颗粒质量相等，采取计数的方法代替质量浓度的测量。因此，过滤效率等于母液和滤出液中颗粒数量之差与母液中颗粒数量的比值，初始过滤效率（η0）定义为：在零时刻，滤层处于干净状态下的过滤效率。

图 3 DEM计算区域和边界条件Fig. 3 Illustration of the DEM solution domain and boundary conditions

表 1 几何参数和模拟条件Table 1 Geometric parameters and the simulated conditions

表 2 模拟计算参数Table 2 Parameters in the simulated calculation

为了验证三维模型的CFD-DEM耦合计算方法的可行性，设计了小型过滤床实验，如图4所示。过滤床内的滤料是直径为1.0 mm的氧化铝颗粒，过滤床高为20 mm，杂质粒径分别选取100 μm和200 μm。不同表观过滤速率下过滤床对杂质颗粒的过滤效率的实验值与计算值比较如图5所示。由图5可以看出模型计算值和实验测量值基本吻合。

3 结果与讨论

3.1 初始过滤效率和床层压降

图6所示为双因素（杂质颗粒大小、表观过滤速率）下初始过滤效率的预测值和实验值的比较，以及由计算值拟合得出的多项式曲面过滤效率的预测。由图6可以看到，过滤床层对较大杂质颗粒的拦截效率高，表观过滤速率的增大会减小过滤效率。由于过滤效率与杂质颗粒直径（dp）、表观过滤速率有关，初始过滤效率可以写成：

图 4 过滤床示意图Fig. 4 Schematic diagram of the filter bed

图 5 过滤效率的实验值与模拟计算值比较Fig. 5 Comparison of the filtration efficiencies between experimental and simulated data

图 6 双因素下的初始过滤效率： (a) 拟合曲面；(b) 预测值与实验值Fig. 6 Initial filtration efficiency under two factors: (a) Fitting surface; (b) Predicted value and experimental value

利用Matlab的多项式函数进行非线性拟合，在130 mm/s＜uA＜330 mm/s，100 μm ＜ dp＜ 200 μm范围内得到初始过滤效率的计算式，如式（14）所示。需要注意 的是，式（14）中 uA、 dp的单位分别是mm/s、mm。

绘制出计算值的拟合曲面，并与100 μm和200 μm时的实验值进行对比，结果如图6（a）所示。

由回归分析得到的拟合曲面和多项式函数关系，进行初始过滤效率的预测。如图6（b）所示，100 μm和200 μm杂质颗粒的过滤效率，其预测值和实验值的偏差均在10%以内。Ergun方程[22]可以用来计算填料床层的压降值，如式（15）所示。

其中： ψ 为滤料颗粒的球形度，取值为1； de为滤料颗粒的体积等效直径，与 dp相等。

图7所示为过滤床的初始压降。由图7可得干净床层的初始压降的模型计算值与Ergun方程所得压降值的比较。Ergun方程的计算误差在±25%以内有效，虽然模型计算值相比Ergun方程计算值偏小，但误差都在25%以内。

3.2 过滤速率对颗粒沉积的影响

杂质颗粒在过滤中的一个重要行为是聚集。当杂质颗粒聚集在一起形成不同的簇（聚集体）时，则发生聚集现象。新形成的簇中的每个单独的颗粒将失去动力学独立性。较大的颗粒（聚集体）在过滤期间更可能沉积。当已经沉积的颗粒受到流体的冲刷、颗粒的碰撞等外力作用时，有可能再次剥离沉积表面，这些现象显然与过滤速率有关。

图 7 过滤床的初始压降Fig. 7 Initial pressure drop of the filter bed

图8示出了过滤一段时间（10 s）时，不同表观过滤速率下同一纵截面的颗粒沉积情况（蓝色代表颗粒速度接近零，红色代表颗粒速度最大，绿色代表颗粒介于两者之间，全文同），流体走向为从上至下。可以看到，当表观过滤速率为130 mm/s时，悬浮颗粒在滤料间隙成团的比例较大；当表观过滤速率增大到230 mm/s时，“团簇”的颗粒数量减少；当表观过滤速度增大到330 mm/s时，几乎看不到成团，且颗粒沉积量明显减少。

图 8 不同表观过滤速率下同一纵截面的颗粒沉积Fig. 8 Particle deposition in the same longitudinal section at different superficial velocities

在EDEM的后处理界面，使用“Clipping”功能将过滤床模型等分为20层，并逐个统计每层中杂质颗粒数随时间的变化关系，然后单独统计某一时刻不同层的颗粒数量。图9所示为相同时刻（10 s）、不同表观过滤速率下过滤床不同深度处的杂质颗粒数量。可以看出，随着表观过滤速率的增大，每层的杂质颗粒数量有减少的趋势，总的沉积颗粒数量相应减少。所以，表观过滤速率会影响杂质颗粒在滤层中的成团趋势，会使得已经沉积的颗粒再次返回流体中，进而影响滤层的过滤效率。

图 9 不同表观过滤速率下不同床层深度处的杂质颗粒数量Fig. 9 Numbers of particles at different bed depths and superficial velocities

3.3 过滤时间对颗粒沉积的影响

图10示出了过滤床中不同深度的杂质颗粒随时间的沉积变化情况。将同一时刻（分别选取1，5、10 s）、不同颗粒层内的杂质颗粒数统计到一起，构成杂质颗粒数随过滤床深度的变化图，如图11所示。图10和图11均在表观过滤速率为130 mm/s，杂质颗粒直径为100 μm的条件下进行，其他条件参见表1。其中图10是分别选取距滤层顶部0、10 mm和19 mm截面处的杂质颗粒沉积情况。

从图10和图11可以看到，沿着过滤床从上到下，杂质颗粒数量呈现减少趋势；靠近入口处的杂质颗粒成团趋势明显；在靠近出口处，速度较大的杂质颗粒数量比例较大。再对比不同时间的颗粒沉积情况，取前3层作为表层，后3层作为底层。图11（a）中柱状分布较“尖锐”，表层颗粒数超过45%，而底层颗粒数仅为4%；图11（b）、（c）中柱状分布趋向平缓，表层颗粒数分别为29%、25%，底层分别为14%、17%，表层和底层颗粒数的差值进一步减少。从而说明，在过滤前期，靠近入口处的杂质颗粒数量变化较大；当过滤逐渐加深，靠近入口处的杂质颗粒数量变化减弱，较深的床层处也逐渐积累更多的杂质颗粒，最后床层积累的杂质颗粒数量越来越均匀。

图 11 过滤床内不同时刻不同深度处的杂质颗粒数占所有沉积颗粒数的比例Fig. 11 Ratio of impurity particles to all deposited particles at different depths of the filter bed at different time

3.4 颗粒沉积对过滤效率和压降的影响

随着过滤时间的进行，杂质颗粒的沉积会影响过滤床内的流体流动、过滤效率和床层压降。杂质颗粒随流体进入过滤床层，受到滤料颗粒的拦截作用，一部分杂质颗粒沉积在滤料颗粒表面，另一部分受前面颗粒架桥作用的影响而被拦截，但仍有被流体冲刷而脱离的可能性。其他的颗粒随流体被带出过滤床层。沉积在床层的颗粒会改变后来流体和颗粒的流道，影响其流动速率和被拦截的概率，进而影响过滤床层的过滤效率和压降。

图12（a）、12（b）分别为在表观过滤速率为230 mm/s、连续过滤条件下，过滤床的过滤效率和床层压降随单位体积沉积量增长的变化情况。可以看到，对于不同大小的杂质颗粒，过滤效率均有所提高，但是床层压降也随之增大。杂质颗粒的沉积会改变滤料床层的流体通道大小，使得通道尺寸小于杂质颗粒尺寸的可能性更大，因此床层拦截的杂质颗粒更多，过滤效率更高。小颗粒之间的缝隙比大颗粒之间的缝隙小，流体流过小颗粒之间的缝隙时的速率更大，并且床层空隙率也随颗粒沉积的增多而变小。床层压降与流速和空隙率是正相关关系，所以小颗粒杂质沉积产生的压降就较高。

图 12 过滤效率和床层压降随杂质颗粒单位体积沉积量的变化Fig. 12 Filtration efficiencies and pressure drop variations with the particle deposition mass per unit volume

4 结 论

采用CFD-DEM耦合方法，对高20 mm、颗粒直径为1.0 mm的三维随机堆积颗粒床层进行固液过滤研究。主要结论如下：

（1）随着床层的加深，压降值逐渐增大；流体因为颗粒之间的缝隙减小而加快速度。与Ergun方程相比，不同表观过滤速率下的三维模型压降值的计算偏差在25%以内。

（2）过滤床层对较大杂质颗粒的拦截效率高，表观过滤速率的增大会减小过滤效率。通过回归分析得到拟合函数，并进行初始过滤效率的预测，对于100 μm和200 μm的杂质颗粒的过滤效率，其预测值和实验值的偏差均在10%以内。

（3）沿着过滤床从上到下，杂质颗粒数量呈现减少趋势；靠近入口处的杂质颗粒成团趋势明显；在靠近出口处，速度较大的杂质颗粒数量较多。在过滤前期，靠近入口处的杂质颗粒数量变化较大；当过滤逐渐加深，靠近入口处的杂质颗粒数量变化减弱，床层积累的杂质颗粒数量越来越均匀。

（4）随着过滤时间的增加，杂质颗粒的沉积会影响过滤床内的流体流动，进而影响过滤效率和床层压降。对于不同大小的杂质颗粒，随着时间进行，其过滤效率均有所提高，但是引起的床层压降也随之增大。