黏性颗粒在聚偏氟乙烯纤维膜内部沉积特性数值模拟

诸文旎，刘倩倩，JAMSHAID Hafsa，李 莉，祝国成

(1.浙江理工大学纺织科学与工程学院(国际丝绸学院)，杭州 310018；2.Department of Knitting, National Textile University, Faisalabad 37610, Pakistan；3.香港理工大学纺织与服装学院，香港 999077)

大气中的悬浮颗粒物对人体的健康有着很大的威胁，尤其是粒径小于5 μm的颗粒会进入人体肺部而不易排出，极大地增加了患呼吸道疾病的风险[1]。通过空气过滤纤维膜对颗粒物进行过滤与捕集是当前改善空气质量与人们生活质量最为有效的方法。

当前，空气过滤研究主要是针对无黏颗粒开展，关于黏性颗粒研究较少。然而，在潮湿环境下空气中飘浮的颗粒物会吸附水膜，使颗粒间产生黏性力，这使得其过滤机理有别于普通无黏颗粒。据国家统计局数据，2019年中国34个主要城市年平均湿度为65.12%[2]，这表明这些城市一年中大部分时间均处于非干燥环境，使得空气中的黏性颗粒含量高，因此研究黏性颗粒沉积特性与过滤机理对改进纤维膜的过滤性能有重要意义。

早期的纤维膜过滤机理主要围绕单根圆柱纤维展开研究[3]，与实际的纤维集合体形态差别较大，因而许多学者针对复杂结构的纤维膜过滤机理，通过计算流体力学(Computational fluid dynamics, CFD)进行了大量纤维过滤过程的研究，在构建纤维过滤材料三维模型上，主要有两种方式，其一是对实际纤维材料的测试照片进行技术处理再识别后转换为三维模型；其二是直接通过计算机建模。Sambaer等[4]提出的方法基于扫描电镜图识别顶部纤维的中心线，再旋转生成三维纤维，通过不断叠加纤维层组成具有一定厚度的三维模型；Faessel等[5]对木质纤维板进行X射线成像，通过程序分析、提取出纤维的长度、曲率、直径等信息再转换为三维模型；Tafreshi等[6]模拟了微观结构纤维材料的三维模型，计算了双模态介质的渗透率，并为双模态介质建立单模态等效直径，从而利用现有单模态材料的表达式进行过滤效率的预测；Wang等[7]模拟了不同规律排列方式的多层纤维过滤器的过滤过程，表明了交错模型的捕集性能优于平行模型；Tao 等[8]用离散单元法(Discrete element method, DEM)研究了纤维阵列排布方式对过滤性能的影响，发现不同的纤维阵列排布方式会影响压降和过滤效率；Hosseini等[9]通过程序建立相同直经的三维层状结构模型并模拟了纳米纤维膜在过滤气溶胶粒子时的压降和收集效率；钱付平等[10]采用随机算法生成在水平面相互交错的单层纤维，再将多个纤维层平行叠加合成三维模型，基于离散单元法耦合模拟了其气固流动情况，结果表明耦合模拟结果符合已有的经验模型，此外，纤维膜孔隙率、流速和颗粒大小对纤维膜过滤性能有重要影响。与实际的纤维膜相比，基于材料成像建模的方法能较好地还原单根纤维体的形貌，而计算机直接建模的方法从纤维的排列、空间分布等角度出发，能更好地调控模型整体，但都与实际的纤维膜存在差异。

为了使三维模型更贴近实际情况，实现局部与整体的平衡，本文通过观察、测量真实纤维膜材料，总结其形态结构特征后，通过Matlab编写程序，建立具有随机直径、随机空间分布特征的三维纤维模型。该模型不需要进行多次转换或叠加，可一次性生成，建模过程相对简便，其次该模型虽然直接由程序生成，但能够针对不同的滤材，自由调控纤维直径，灵活性大，通用性强。

在纤维空气过滤膜材料中，聚偏氟乙烯(Poly-vinylidene fluoride, PVDF)材料耐腐蚀性好，抗拉伸强度高，在空气过滤领域有较大的发展潜力。因此本文在观察总结PVDF纤维膜的结构特征后，建立相应的三维模型；在离散单元法中引入(Johnson-Kendall-Roberts, JKR)接触模型，描述黏性颗粒在PVDF纤维膜表面的沉积行为，探究表面能、碰撞恢复系数和粒径对黏性颗粒沉积特性的影响，为深入研究黏性颗粒的纤维介质过滤机制提供理论参考。

1 基于实际PVDF纤维膜建立三维模型

为了获得纤维过滤介质的三维模型结构，本文首先对静电纺PVDF纤维膜扫描电镜成像，观察与测量对比并进行分析；根据该实际材料的结构特征，采用Matlab随机算法进行建模，这样建立的模型与实际情况具有很好的相似性。

1.1 PVDF纤维膜结构分析

图1为静电纺PVDF纤维膜的电镜照片。图1中纤维表面基本光顺，纤维排列与堆叠不规则，纤维轴向无序；将照片局部放大可以看出纤维排列较为紧密，轴向弯曲度小，膜的孔隙小。随机选取100个数据点测量纤维直径，对结果进行统计并计算平均值。纤维直径分布情况如图2所示，直径范围在 0～4.5 μm，纤维直径基本遵循正态分布，平均直径为1.875 μm。经测量，纤维膜填充率为17.13%±0.22%，厚度为(50.30±0.58)μm。

图1 PVDF纤维膜电镜

图2 纤维直径分布情况

1.2 纤维膜模型建立

针对实际PVDF膜的形貌与结构特征，在Matlab中编写相应的程序，生成不同纤维直径、厚度、方向的纤维结构。模型中将单根纤维形态简化为圆柱体。其次，由于静电纺丝的工艺受到许多因素的干扰，特别是纺丝液在射流阶段中溶剂挥发，纺丝液黏度增加，阻碍纤维的拉伸，导致纤维直径的不均匀，以及纤维固化、冷却后形态多样，从断面上看并不是所有纤维都处于或接近水平状态，而是存在不规律的角度或轴向弯曲，因此除了横向平面上的随机分布以外，设定纤维模型在纵向截面上随机分布，纤维与水平面存在无序的绕轴角度。

模型建立逻辑为：由random函数生成随机点坐标，任意两个非同面点坐标随机连线得到纤维轴中心线；生成符合正态分布的随机数矩阵作为纤维直径值；不同的纤维轴线和直径值随机组合生成各个纤维圆柱体，共同构成一个三维模型。

在可输入参数中，由模拟单元(设定长、宽(厚度)、高的模型外框架)来限制纤维长度；厚度、直径与纤维根数来调控固体体积分数(Solid volume fraction, SVF)的大小。其次，直径数值不得超过厚度值，以保证模型的完整性。纤维直径的改变不影响模型空间结构形态，模拟单元尺寸为80 μm×50 μm×80 μm，纤维直径为10～15 μm，SVF为20.15%，如图3所示。

图3 三维纤维模型

2 CFD-DEM数值计算

2.1 黏性颗粒受力模型

在气固两相流中，对于黏性颗粒体系，其受力主要有重力、曳力、接触力和黏性力，其中，黏性力可分为接触式黏性力与非接触式黏性力，接触式黏性力包括毛细力，非接触式黏性力包括范德华力与静电力。常见非干燥环境中，颗粒表面电荷数量少，易流失，产生的静电力可忽略不计；随着湿度增大，颗粒在空气中极容易吸水，表面凝聚水分子增多，毛细力将大大增强颗粒间团聚作用，这种情况下毛细力往往是范德华力的十几倍或几十倍[11]，因此颗粒黏性力仅考虑毛细力，忽略范德华力。

采用JKR法向接触模型，该模型将颗粒球体定义为软球，可以模拟颗粒间、颗粒与纤维膜的接触及弹性形变，并准确描述细微颗粒间的黏附作用[12-13]。为简化模型构建，Chokshi等[14]在JKR模型基础上建立了无量纲表达式来描述法向接触力与接触半径、法向重叠量之间的关系。

采用Mindlin-Deresiewicz切向接触模型[15]，当两个颗粒接触时发生滑动摩擦，两颗粒的接触点存在不为零的相对速度，滑动摩擦力总是小于最大静摩擦力；当该相对速度为零时，转变为滚动摩擦，滚动摩擦阻碍颗粒继续滚动，使其速度减小直至停止滚动，颗粒间摩擦转变为静摩擦，因此，颗粒间的滑动摩擦与滚动摩擦最终都会转变为静摩擦状态，因此切向接触力主要考虑为静摩擦力。

图4为颗粒接触模型，颗粒因碰撞发生弹性形变和相对滑移，产生法向接触力与切向接触力。图4 中，R为颗粒半径，m；a为接触半径，m；δn为法向重叠量，m；δt为切向位移量，m。

图4 颗粒接触模型

两物体间法向接触力Fn的计算基于材料表面能、相互接触参数以及接触半径得出，如式(1)：

(1)

当颗粒处于静态黏附的平衡状态时，此时黏性力Fc如式(2)—式(5)：

Fc=3πγR

(2)

(3)

(4)

(5)

式中：E为杨氏模量，Pa；γ为表面能，J/m2；a0为静态黏附时接触半径，m；δn0为静态黏附时法向重叠量，m。

切向接触力Ft如式(6)：

(6)

式中：μs为静摩擦系数；δt,max为切向位移量最大值，m。

2.2 数值模型的建立

CFD-DEM耦合采用欧拉-拉格朗日模型，其中欧拉方法描述流体相，拉格朗日方法描述颗粒相，该模型适用于颗粒相固体体积占比较小的情况。通过SIMPLE算法求解，流体连续性方程、动量守恒方程与本构方程分别如式(7)—式(9)：

(7)

(8)

(9)

式中：ρ为流体密度，kg/m3；v为流体速度，m/s；μ为流体动力黏度，Pa·s；p为流体压力，Pa；t为时间，s；σy为y方向应力，Pa。

忽略范德华力，颗粒受力与运动方程如式(10)和式(11)：

(10)

(11)

式中：mp为颗粒质量，kg；vp为颗粒运动速度，m/s；FA为颗粒所受曳力，N；FC为黏性力，N；Fn为颗粒间法向接触力，N；Ft为颗粒间切向接触力，N；Ip为颗粒转动惯量，kg·m2；ωp为角速度，r/s；Mi为颗粒与颗粒或其他固体壁面间的接触力矩，N·m。

模拟选用欧拉方法作为耦合方法，耦合设置中，为增加计算稳定性，减小了松弛系数，设定动量松弛系数为0.7，体积松弛系数为0.6。

CFD部分采用Fluent软件对气相流场进行瞬态计算求解，在完成一个时间步迭代后计算颗粒曳力，在这个结果的基础上，EDEM软件进行固相计算求解，更新颗粒运动状态并传递至Fluent进行下一次计算，由此循环，使气固两相流间质量、动量与能量进行交换，实现完全双向耦合。

2.3 边界条件及计算工况

计算域模拟流体为空气，流动的雷诺数Re计算式如式(12):

(12)

式中：l为特征长度，m。根据模拟设置参数，ρ=1.225 kg/m3，μ=1.789 4×10-5Pa·s，l采用流场总长度, 计算域前后流场长度分别为25 μm，过滤区域纤维厚度为50 μm，即l=100 μm。根据预设流场入口速度0.1 m/s，计算可得Re ≈ 0.68，而层流的判定条件是Re<2 300则流动为层流，因此流场内气流为层流状态，流动方式为不可压缩层流。

图5为模拟区域的边界条件设置。计算域入口速度为0.1 m/s，出口表压为0，四周是壁面，纤维介质表面采用无滑移边界条件。颗粒生成面也位于计算域入口，颗粒初始速度0.1 m/s，以1×105个/s的速度产生，注入到计算域后随气流运动，直至被拦截捕获或逃逸至出口。

图5 边界条件

工况设置中，CFD时间步长为5×10-6s，DEM步长为5×10-8s，模拟时长为5 ms。如图6所示，网格无关性检验中，纤维一周的点数在22以后纤维模型压降较为稳定，综合考虑计算成本与计算精度的影响，选择纤维一周网格节点数量为22点，最终网格数约为 22 万。

图6 网格数量对压降的影响

碰撞恢复系数是反应颗粒非弹性碰撞程度的量，碰撞系数的改变会影响碰撞后颗粒振动的剧烈程度，从而影响颗粒的流动特性[16]；固体表面能γ影响表面润湿性，从而影响物体间接触黏性力，因此将这两个参数作为自变量。材料的物理属性与接触参数设置如表1与表2所示。颗粒属性参考粉煤灰，为常见的工业粉尘之一；纤维模型主要考虑其空间结构，其属性设置能够体现材料的差异，与仿真结果无联系，可按需进行相应设置，表2中接触参数的选取参考了前人数值模拟的经验[18]。

表1 材料物理属性参数

表2 接触参数

3 捕集效率的经验公式

颗粒的捕集是拦截、扩散、惯性作用等多种捕集机理共同作用的结果，对于单纤维的联合捕集效率，可以Kuwabara模型[17]计算，如式(13):

(13)

式中：α为纤维过滤区域的填充率，%；E∑为单纤维过滤效率；H为纤维介质厚度，m；df为纤维直径，m。忽略重力影响，E∑表达式[19]如式(14):

E∑=1-(1-ER)(1-EI)(1-ED)

(14)

式中：ER、EI、ED分别为拦截、惯性、扩散作用时的单纤维过滤效率。ER、EI、ED的表达式如表3所示[20-22]。

表3 捕集效率表达式

4 结果与分析

4.1 表面能对颗粒沉积的影响

图7是粒径dp为4 μm的颗粒分别在表面能为0.01、0.03、0.05 J/m2时在纤维介质表面沉积情况。可以看出，被捕获颗粒在纤维体表面的沉积并不是均匀的，主要集中在纤维体迎风面杂乱分布，后方的纤维体主要起到补充拦截作用，其次纤维表面有较多的树突结构颗粒沉积。0.01～0.05 J/m2的范围内，表面能数值减小，颗粒受到黏性力小，但颗粒沉积数量却有所增多，可能的原因是颗粒表面能较小时，颗粒内外受力不平衡的程度小，其受力相对稳定，最早进入流场的颗粒在被纤维体拦截捕获后黏附更多颗粒，形成稳定的树突结构沉积，从而对颗粒有较高的捕获效率。

图7 粒径为4 μm颗粒沉积情况

以表面能为0.01 J/m2时的工况为例，颗粒在纤维表面的树突沉积形态如图8所示，与图9由Schilling等[23]对纤维表面颗粒沉积与团聚的观测结果符合，说明本文对颗粒沉积情况的模拟是可靠的。

图8 表面能为0.01 J/m2时颗粒沉积形态

图9 纤维表面颗粒沉积与团聚的观测

图10是恢复系数为0.5，表面能为0.01 J/m2时，不同颗粒粒径和捕集效率的关系，由于粒径增大，惯性碰撞效应增强，因此颗粒捕集效率增大，将模拟结果和计算值比较，可以看到总体趋势是一致的，其中，2、3 μm颗粒的模拟值略小于计算值，原因可能是Kuwabara模型除了考虑拦截作用、惯性碰撞的影响，还考虑了布朗扩散作用，而由于2、3 μm颗粒的粒径小，布朗运动会更加剧烈，因此计算值较高。

图10 计算值与模拟值的对比

如图11所示，直径为4 μm颗粒在不同表面能下的颗粒捕集情况。可以看出，表面能不变，捕集量随时间逐步增加；颗粒表面能小，捕集数量更多。初期进入流场中的颗粒随气流运动一段距离后与纤维发生碰撞和黏附，因此模拟前期主要是纤维拦截导致颗粒沉积，而后期进入流场的颗粒减少了与纤维体发生碰撞黏附的几率，大量颗粒是因沉积层的黏附作用而被捕获，并且黏附在纤维表面形成树突结构沉积，进一步提高了过滤效率。而表面能的增大会引起颗粒内外受力不平衡的程度增大，使颗粒的受力不稳定，导致其不易黏附或是已形成的树突结构不稳定，因此颗粒容易因曳力或其他颗粒碰撞等作用脱离黏附，继续随气流运动，最终导致颗粒捕集数量有所减少。

图11 不同表面能下颗粒捕集数量

图12为不同表面能工况下，直径4 μm颗粒间接触矢量随时间变化的情况。对于表面能为 0.03 J/m2以下的颗粒，模拟过程中颗粒受力较稳定，互相接触颗粒的沉积团聚状态稳定，其接触矢量也逐渐平稳，整体趋于稳定状态耗时较短，沉积颗粒的挤压程度小。当表面能数值在0.04 J/m2及以上时，模拟初期，颗粒流刚进入流场，只有少量颗粒被拦截与黏附，接触数量少，因此颗粒间接触矢量小，随着更多颗粒被拦截、黏附与沉积，颗粒处于受挤压状态，接触矢量逐渐增大，此外较高的能量导致颗粒受力不平衡，模拟过程中颗粒受到的接触力波动较大，整体趋于稳定状态的耗时更长，沉积团聚的颗粒堆积挤压程度也较大。

图12 不同表面能下颗粒间接触矢量

4.2 恢复系数对颗粒沉积的影响

图13给出了颗粒粒径与恢复系数的最终捕集效率关系，工况中颗粒表面能均为0.01 J/m2。从最终捕集效率来看，颗粒恢复系数的变化对直径为2、3 μm的颗粒影响较小，随着恢复系数增大，其捕集效率稍有增加；对粒径为4 μm及以上的颗粒，在恢复系数为0.1至0.5的范围内增幅较大，特别是粒径为5 μm的颗粒，增量超过50%，后有些许波动。综合来看，颗粒恢复系数为0.5时颗粒捕集效率的吻合度最好，总体趋势是先增加后趋于稳定。恢复系数较小时，颗粒间碰撞的剧烈程度低，颗粒形变恢复慢，接触作用时间长，黏附作用强，模拟中被纤维拦截黏附的颗粒容易因其他运动颗粒的碰撞-黏附作用而脱离沉积状态，继续随气流运动，纤维体上不易产生颗粒沉积层，因此总体对颗粒的捕集能力有较大影响，当恢复系数在0.5及以上时，颗粒间发生接触碰撞相对剧烈，颗粒恢复形变快，黏性力的作用时间短，已沉积颗粒不易被带走，滤饼过滤作用加强，因此捕集效率增加，综合来看，颗粒碰撞恢复系数为0.5时捕集效果最好。

图13 不同恢复系数的捕集效率

图14为模拟过程中，粒径为4 μm颗粒所受到的颗粒间平均接触力。从图14中能够看出，接触力均大于0，说明过程中颗粒间的接触力均为斥力[18]。颗粒恢复系数较小时，沉积颗粒数量少，并且陆续回到运动状态，极大减少了颗粒间互相接触的几率，同时两两接触的颗粒中，一方处于已沉积或即将沉积的低速运动状态，另一方从颗粒入射面出发进入计算域，处于高速运动状态，两者发生碰撞时速度差大，又因恢复系数小，颗粒间趋于非弹性碰撞，两者趋于黏合在一起运动，表面黏附力大，为保持受力稳定就要互相给予较大的斥力；恢复系数较大时，两颗粒趋于完全弹性碰撞，颗粒接触时间短，形变小，颗粒间黏附力较小，因此接触力也较小，同时由于黏性力小，已沉积颗粒不容易被高速运动颗粒带走，能使纤维表面的沉积速率保持稳定，使捕集效率明显提高。

图14 不同恢复系数下颗粒间法向接触力

4.3 粒径对颗粒沉积的影响

图15为颗粒表面能0.01 J/m2，恢复系数0.5情况下，不同粒径黏性颗粒在纤维体上的过滤沉积情况，可以看出，粒径对颗粒沉积影响较大。直径小的颗粒质量小，在模拟过程中多数颗粒从纤维模型空隙中穿过，同时容易因曳力作用做绕流运动，不易被捕获，因此颗粒沉积较少；随着dp增大，质量、体积显著增大，颗粒受曳力影响小，纤维与沉积层拦截捕获作用显著，沉积与团聚的颗粒数量也更多。

图15 不同粒径颗粒沉积情况

图16是在5 ms时刻，颗粒表面能0.01 J/m2，恢复系数为0.5的工况下，流场内平均颗粒间法向接触力、颗粒与纤维间法向接触力情况，在经过了 5 ms 时长的运动，一部分颗粒已经穿过纤维逃逸，一部分颗粒已经沉积在纤维表面，还有部分颗粒刚从入口发射，正处于流场前端，无论是运动的还是已经沉积的颗粒，都可能与其他颗粒或纤维存在接触，这种接触主要是沉积颗粒间、沉积颗粒与纤维间的接触，并且力的大小与颗粒粒径有关。

图16 不同粒径下法向接触力情况

随着粒径增大，颗粒间接触力先增大，后逐渐减小，颗粒与纤维间的接触力随粒径的增加而增大。在粒径为2 μm时，颗粒较多作绕流运动，碰撞情况很少，颗粒间接触力、颗粒与纤维体接触力很小；可以看到粒径为3 μm时颗粒接触力的均值和离散程度都是最大的，相比2 μm粒径，3 μm粒径的颗粒质量、体积都有所增加，受曳力影响稍有减小，颗粒间接触与碰撞次数更多，模拟过程中，该粒径颗粒在纤维表面有少量沉积，但部分已沉积颗粒会受曳力影响而脱离纤维表面逃逸，而在这种运动-沉积-运动的过程中一直伴随着颗粒的接触与黏附，导致颗粒间接触与碰撞非常不稳定，所以其受到的力的离散程度最高；从4 μm粒径起，随着粒径增大，接触半径增大，而接触半径增大的程度是低于粒径增大的程度，因此颗粒间接触力减小，离散程度也有所减弱。对纤维体而言，颗粒与纤维间接触力随粒径非线性增大，这是由于颗粒质量增加，惯性作用更明显，颗粒对纤维体的碰撞接触作用增强，接触力增大。

5 结 论

通过扫描电镜成像获得的PVDF膜形态特征为基础，编写程序以生成直径随机、分布随机的纤维过滤介质三维模型，通过气固两相流耦合方法，对黏性颗粒的运动与沉积进行模拟，研究了不同表面能、颗粒间碰撞恢复系数和粒径对颗粒间、颗粒与纤维间接触与黏附作用、颗粒沉积及捕集情况的影响，得出以下结论：

a)颗粒的表面能是影响其沉积形态的因素之一，在0.01～0.05 J/m2的范围，相同粒径下，表面能为0.01 J/m2时，颗粒沉积数量最多，颗粒捕集效率高，同时沉积颗粒间的挤压程度小，颗粒受到的接触力波动小，接触矢量大，整体趋于稳定状态的耗时短。

b)颗粒恢复系数的变化对直径为2、3 μm的颗粒影响较小，对粒径为4 μm及以上的颗粒，恢复系数增大，捕集效率先增加再减少后大致趋于稳定状态，恢复系数为0.5时综合捕集效果最好，7 μm颗粒过滤效率为95.4%，此时颗粒处于不完全弹性碰撞，沉积团聚形态稳定，黏附效率高。

c)黏性颗粒的捕获受粒径影响较大，粒径为 2 μm 时颗粒沉积量小，颗粒法向接触力小；粒径在3 μm及以上时，随着粒径增大，拦截捕获作用明显，颗粒更容易堆积堵塞纤维模型空隙使拦截过滤作用显著，因此黏附效率明显提升；同时粒径增大，接触半径减小，颗粒间接触力减小，而惯性作用增强，与纤维间接触力随粒径增大呈非线性增大。