离心泵内纤维悬浮流的PIV实验及CFD-DEM模拟研究

张启华 张为栋 闫召旭 谢之昂 康 顺

（江苏大学流体机械工程技术研究中心，江苏镇江，212013）

采用离心泵输送纤维悬浮液在制浆、造纸流程及废水处理过程中应用广泛。尤其近年来，随着纸张消费量的上升，作为纸浆循环的动力设备，纸浆泵所产生的能耗约占工艺流程总能耗的15%。同时，泵的稳定运行对造纸终端产品的质量也至关重要。因此，深入开展纸浆泵内部流动的研究，对了解长纤维悬浮液的流动特性及推进纸浆泵节能降耗均有重要的工程意义。

纤维悬浮液的运动规律以及纤维的取向分布状态造成的流变特征，是纤维悬浮液研究的两个主要方面。可通过引入黏度作为剪切速率的函数来表征纸浆纤维悬浮液的非牛顿行为。Cotas等［1］提出了一种拟均匀方法来研究管道中桉树纤维悬浮液的湍流。目前，大多数泵内两相流模拟采用基于Mixture的模型，即视两相均为具有连续速度场的流动介质。Liu等［2］基于Mixture模型模拟了粒径1～12 mm的颗粒在不同工况下的运动状态并总结了不同粒径粒子的分布规律。Zhang等［3］基于Mixture模型研究了颗粒浓度、体积分数及粒径对泵性能的影响。但Mixture模型忽略了颗粒间的碰撞以及颗粒与流体的相对运动，存在一定的模拟误差。且颗粒体积的变大或颗粒浓度的提高，颗粒的拓扑结构、自身运动以及颗粒间相互作用不可忽略。随着待解决问题的复杂程度越来越高，具有较高计算效率的离散元（DEM）法在泵两相流工程中得到了推广应用。基于丰富的颗粒表征、更复杂的颗粒初始化，以及颗粒-颗粒与颗粒-部件的接触模型，DEM法能更直观地展现各类复杂流程中颗粒的运动。Varga等［4］使用流体力学（CFD）-DEM方法模拟管内固液两相流动并与实验中不同材料颗粒对管道的磨损情况结果对比，发现材料属性对管壁磨损程度影响较大。Huang等［5］采用CFD-DEM方法具体分析泵内颗粒的运动情况，发现颗粒趋向于在叶片工作面聚集。以往泵工程的CFD-DEM研究主要以球颗粒为主，而对于长颗粒方面的研究仍鲜见报道。

目前，两相流可视化中颗粒相的追踪主要采用高速摄影［6］。近年来，粒子图像测速技术（PIV）以其非接触式测量方式，通过瞬态记录大量空间点上的速度分布信息，从而能够提供丰富的流场空间结构以及流动细节。通过PIV测试技术，科研人员已在离心泵的清水流场［7-8］以及基于玻璃球颗粒固液两相流实验测量方面［9-11］取得较大的进展。宋天雄［12］采用PIV及MATLAB软件后处理技术，针对纸浆纤维悬浮液在磁力搅拌器流场中的流变特性开展实验研究，但仅就悬浮液的涡量进行了分析。在纸浆泵方面主要集中在性能实验研究，进而探讨纸浆浓度及含气率对泵性能的影响［13-14］。Vaezi等［15］研究了输送麦秸和玉米秸颗粒对离心泵性能的影响，发现了颗粒独特的摩擦损失行为（即减阻效应）而引起泵效率提高的现象。叶道星等［16-18］对中浓浆泵及湍流发生器开展了一系列性能研究，这些研究都由于浆液不透明等因素，无法进一步通过PIV实验对内流场进行可视化研究。在纸浆流场研究中，由于纸浆到达一定浓度时会缠结成絮状及纸浆纤维膨胀的特性，会导致激光无法穿透，给PIV实验的追踪测量增加难度。随着时间推移，纸浆会产生发酵、变质等现象，不易储存。并且若要控制纸浆纤维尺寸相对均匀，需对其进行机械或化学处理，工作量大且会增加新的化学元素。自然材料中，头发纤维的性质很接近纸浆纤维，其具有较好的悬浮性且对激光不敏感，因此本研究采用头发纤维作为实验材料以便进行后续的PIV可视化研究。

本研究以头发纤维悬浮液作为研究对象，针对离心泵内细长纤维悬浮液内流场，首次采用CFD-DEM耦合计算及PIV可视化实验相结合的方法，研究了三叶片开式叶轮离心泵内不同工况、不同轴截面以及不同质量浓度的长纤维悬浮液的内流场流动特性。通过流场的对比分析，探讨长纤维对流态畸变的作用；同时，也对泵外特性进行了测量和分析。

1 耦合仿真分析

CFD-DEM数值模型的耦合计算包括两个步骤，首先，通过EDEM软件最大程度地反映颗粒特征，模拟工程中纤维颗粒真实形态，并生成研究所需的颗粒群；然后，通过与FLUENT耦合对泵内流场进行CFD模拟，建立纤维颗粒与流体间的信息传递。其中，EDEM软件负责纤维颗粒的运动计算，流场通过FLU⁃ENT软件进行耦合求解。

1.1 泵基本参数

本研究采用的是一台扭曲叶片开式叶轮纸浆泵，叶轮叶片数为3，模型泵的主要几何参数为：叶轮进口直径D1=70 mm，出口直径D2=184 mm，叶轮出口宽度b2=17.4 mm，叶片出口角β2=43°，蜗壳基圆直径D3=191 mm，蜗壳出口直径D4=70 mm，其三维模型如图1所示。本研究主要探究离心泵叶轮流道内纤维悬浮液的运动情况，从而在模拟中忽略间隙的影响。

1.2 颗粒模型的建立

颗粒形状与物料属性可通过EDEM软件建模了解。为尽可能模拟头发纤维，采用多个圆球颗粒相连成链状结构，如图2所示。单根头发纤维的长度为6 mm，质量为4.60×10-7kg。

图2 物料颗粒示意图Fig.2 Schematic diagramof solid material particles

1.3 模型控制方程

在本研究的流体模拟计算中，所有计算区域的流体相均在绝对坐标系内。流体的控制方程如下：

式中，ρf、u、t、P、g、μf分别表示流体的密度、流体速度、时间、压力、重力加速度和动力黏度，α表示流体相的体积分数，F表示其他附加力的合力。

本研究中，EDEM软件模拟采用基于欧拉-拉格朗日法追踪颗粒的受力与运动情况，颗粒体积分数均小于10%。颗粒受力情况求解的表达式如下：

1.4 参数设置

在耦合模拟前，需对EDEM软件和FLUENT软件做相关参数设置。本研究利用FLUENT求解流体运动，采用欧拉-拉格朗日法对泵内固液两相流进行数值模拟。考虑到流体与颗粒间的相互作用，采用双向耦合的瞬态模拟。在离心泵入口处定义进口速度为2.03 m s以模拟额定工况下的流体运动情况，入口处湍流强度和湍流黏性比分别为5%和10。出口设置为压力出口，其静压值设置为一个标准大气压。在叶轮入口前和蜗壳出口之后各延伸一段长度为3倍管径的直管以保证颗粒与流体充分混合运动。计算中假设壁面为无滑移壁面，不考虑粗糙度对颗粒的影响。考虑到涡流及旋转的因素，流体计算模型选择RNG k-ε湍流模型，选取SIMPLEC算法来耦合泵内流场的速度和压力，而动量、压力、湍流动能和耗散率的离散化求解均采用二阶迎风格式。

本研究利用EDEM软件求解颗粒相的受力情况，颗粒间的接触选用Hertz-Mindlin（无滑移）连接模型。在EDEM中，泵表面均被定义为有机玻璃，其密度为1200 kg/m3，泊松比为0.4066，剪切模量为1.41 GPa，杨氏模量为3.97 GPa，用来模拟头发纤维的颗粒密度为680 kg/m3，泊松比为0.3。叶轮转速设置为1450 r/min，与FLUENT设置保持一致。在直管进口处设置颗粒工厂，用于计算的颗粒均匀且随机地从直管入口处跟随流体一起运动。定义每秒入射20239、30344及40478个颗粒以表示颗粒质量浓度分别为0.1%、0.15%和0.2%。

通过建立耦合接口传递颗粒与流体的计算数据，FLUENT计算时间步长设为5×10-5s，即叶轮旋转0.435°计算1次，最大迭代步数60，总步数20000。而EDEM中时间步长应小于FLUENT中时间步长且与FLUENT中时间步长成倍数关系，因此，EDEM中时间步长设置为1×10-6s。叶轮每旋转87°，数据保存1次，总时间0.8 s，即叶轮旋转19圈。

1.5 无关性验证

计算网格通过ICEM划分工具创建，考虑到非结构网格对于耦合计算适应性较好等因素，本研究采用四面体非结构网格对计算区域进行网格划分。为验证网格数对计算结果的影响，采用3组不同网格数方案进行计算，结果如表1所示。通过比较，且考虑到计算精度和计算资源的情况下，选用网格B进行后续数值计算。

表1 网格无关性模拟结果Table 1 Simulation results of grid independence

2 PIV实验研究

2.1 实验模型

为验证上述仿真模型的正确性并进一步分析流场，开展了相应的PIV实验研究。本研究采用的实验泵是一台典型的开式叶轮离心式纸浆泵，其主要部件为叶轮、蜗壳、耐磨挡板、轴和密封系统等。为便于进行PIV拍摄，采用半螺旋吸水室结构，而电机轴穿过吸水室与叶轮直连。蜗壳与半螺旋吸入室之间通过4个夹具固定。实验泵整体结构如图3所示，泵主要参数见1.1节。

考虑到有机玻璃透光性好，抛光处理后粗糙度较低等特点，模型泵中蜗壳和叶轮均由有机玻璃制造而成。前挡板使用不锈钢制造，为减小其对PIV激光造成的反射，将前挡板进行涂黑处理。实物图如图4所示。

2.2 实验装置及实验介质

为开展相关测试工作，搭建了一座开式实验台，主要由模型泵、进出口阀门、罐体和管道等组成。采用敞口塑料制罐，罐内采用搅拌器连续搅拌，以便于颗粒均匀混合。管路、连接弯头等全部采用DN80不锈钢管制作。

图3 模型泵结构图及零部件Fig.3 Structural drawingand part drawingof themodel pump

图4 模型泵蜗壳、叶轮和前挡板Fig.4 The volute,impeller，and front plate of the model pump

水泵实验台的测试精度通常用水泵效率测量的不确定度来评估。泵效率测量不确定度（eη）计算公式如下：

式中，eQ为流量测量的不确定度，eH为扬程测量的不确定度，eP为驱动机输入功率测量的不确定度。

实验中使用的DN80流量计经过法定计量部门的检定，符合0.2级基本精度要求，流量测量的不确定度为：

实验中使用的差压变送器经过法定计量部门的检定，则扬程H测量的不确定度为：

实验中使用的电流互感器精度等级为0.2级，即线性误差E=0.2%，则驱动机输入功率测量的不确定度为：

因此，泵效率测量的不确定度为：

综上，该实验台精度能够满足GB/T 3216—2016标准规定的2级精度要求。

PIV实验装置由PIV测试系统及外触发同步控制系统两部分构成。实验中采用美国TSI公司的粒子图像测速系统。该测试系统主要包括：YAG200-NWL型脉冲激光器、610015-SOL型光臂及其片光源透镜系统、610035型同步器、630059POWERVIEW型CCD相机。采用Insight 3G和Tecplot软件对PIV实验图像进行采集处理。

PIV实验采用了具有较好跟随性和散射型的二氧化硅空心玻璃球作为示踪粒子，该粒子直径为20μm，材质密度为1.05 g/cm3。PIV实验采取相平均的方法进行，即在每个相位点拍摄多组图片，使用Insight 3G软件将所拍摄的图片进行平均处理，得到对应的流场信息。本研究通过轴编码器和外触发同步控制系统控制相机在叶轮的叶片掠过蜗壳隔舌位置时进行拍摄。在头发纤维悬浮液PIV实验中，示踪粒子有一定几率被头发颗粒遮挡，继而影响后处理中示踪粒子的识别，因此实验中每个相位点拍摄200张图片以确保实验结果的准确性。

目前，并没有专门的设备来加工纤维颗粒以控制其长度完全一致，所以本实验通过人工筛选处理以确保头发纤维长度严格控制在2～8 mm，制取难度及工作量较大，因此，本研究对质量浓度分别为0.1%、0.15%、0.2%及0.3%的头发纤维悬浮液进行研究，暂不研究更高质量浓度的头发纤维悬浮液。

2.3 实验方案与测量工况

实验时，启动并调节变频器使得离心泵转速稳定在1450 r min。为验证实验台和实验方法的可靠性，需对实验台进行重复外特性实验。通过泵参数测量仪测量参数并求得离心泵的扬程及效率。

使用电子秤测量头发纤维质量，逐步将对应头发纤维加入罐中并开启搅拌器，使头发纤维在罐中分布均匀，进行离心泵叶轮内部流场的PIV测量，分别对0.6Qd、0.8Qd、1.0Qd、1.2Qd、1.4Qd5个工况点进行PIV实验，其中，Qd为实验泵额定工况，大小为33.5 m3/h。图5为PIV实验激光截面效果图。由图5可知，PIV拍摄效果较好，头发纤维的存在并没有遮挡示踪粒子，因此不需要对固相颗粒图像和液相示踪粒子图像进行分离。

实验泵采用了空间扭曲叶片，叶轮不同轴向位置将对叶轮流道的流场产生影响。为了研究扭曲叶片对叶轮流道内流场的影响，PIV实验选取叶轮区域内3个截面进行拍摄（见图6），这3个截面都垂直于泵轴，分别记为截面Z1、Z2、Z3，其中截面Z2位于叶轮出口宽度中间位置，截面Z1、Z3各距Z2截面6 mm。为表述方便，将隔舌位置处叶片记为叶片1，沿顺时针方向的第二、第三个叶片记为叶片2、叶片3。记叶片1左边流道为L1，流道L2、L3沿顺时针方向类推。

图5 PIV实验激光截面效果图Fig.5 Effect drawing of PIV test laser section

图6 PIV拍摄截面图Fig.6 Cross-section schematic of PIVmeasurement

3 结果与讨论

3.1 离心泵的性能

为得到准确的离心泵外特性结果，对实验泵的扬程、流量进行了重复性实验，头发纤维质量浓度为0.1%时，离心泵的外特性曲线如图7所示。从图7可以看出，3组外特性实验数据变化较小，曲线基本重合，最大误差在4%以内，证明了实验台及实验结果的可靠性。

图8为清水、0.1%、0.15%及0.2%头发纤维质量浓度条件下离心泵性能曲线对比图。从图8可以看出，不同头发纤维质量浓度下，泵扬程变化不明显，但都高于输送清水时的扬程。且随着泵内头发纤维质量浓度的提高，泵效率逐渐增大，纤维质量浓度达到0.1%及0.15%时已初具减阻趋势，纤维质量浓度达到0.2%时已完全显示出泵内因纤维颗粒存在而引起的减阻效应。

图8 不同纤维质量浓度离心泵性能曲线对比图Fig.8 Comparison of pumping performance curvesof different fiber concentrations

3.2 泵内流场验证分析

为了便于分析流体相对于叶轮的运动规律，采用自带的PIV速度处理软件将处理后的绝对速度矢量分解，从而得到叶轮流道内流体的相对速度。在进行PIV实验时，从泵轴左侧打入激光，由于泵轴会遮去部分激光，所以泵轴右侧没有激光，因此在后处理时略去这一区域。

图9（a）显示了通过CFD-DEM耦合计算，相同质量浓度下叶轮转到PIV拍摄相同位置时，Z2截面处流体的绝对速度及相对速度分布；图9（b）给出了通过PIV拍摄得到头发纤维质量浓度为0.1%时，叶轮区域Z2截面处流体的绝对速度及相对速度分布。从图9可以看出，靠近叶片工作面处流体的绝对速度高于背面，从叶轮进口到出口，Z2截面处流体的绝对速度逐渐增大，叶轮出口附近Z2截面处流体的绝对速度有所下降。而Z2截面处流体的相对速度会从叶轮进口至出口处逐渐增大，且在叶轮流道中段靠近压力面附近可观察到低速区域，在此区域内流体的流速低于周围区域内的流速，但差距并不明显。本研究选取相同位置叶片附近流体速度进行分析，并将CFD-DEM耦合计算结果与PIV实验结果进行对比发现，使用CFD-DEM方法得到的流体速度分布及变化趋势与PIV实验结果吻合较好，说明CFD-DEM方法能够较好地模拟泵内复杂流场下的固液两相流运动，验证了采用CFD-DEM方法定性分析泵内流场的可行性，同时反映了PIV实验的准确性。

由图9可以看出，由于遮光、反光等因素影响，L1流道的PIV拍摄效果明显优于L2、L3，所以后续采用L1流道的PIV实验结果对泵内流场进行分析。关于离心泵内清水流场PIV实验，这方面已有许多的研究结论，本文不再赘述。

3.3 测试结果及分析

3.3.1 不同流量下叶轮流道内流体的相对速度分布

为更具体地分析叶轮流道内流体的相对速度随流量的变化情况，分别提取出头发纤维质量浓度为0.1%时，0.6Qd、1.0Qd和1.4Qd工况下Z2截面L1流道内叶片压力面、中间流线及吸力面上流体的相对速度分布并绘制曲线图，结果如图10～图12所示。其中，横坐标L表示距离叶轮进出口的位置，记r0为当前位置半径，r2为叶轮出口半径，0＜r0/r2＜0.6为叶轮轮毂区域，0.6＜r0/r2＜1为叶轮流道区域，因此，L=0代表在r0/r2=0.6位置处，L=1代表在r0/r2=1位置处。

由图10可以看出，在0.6Qd和1.0Qd工况下，叶片压力面附近流体的相对速度随着L的增大呈先减小后增大的趋势，即在流道中段压力面附近会产生低速区，且随着流量逐渐增大，低速区强度和面积逐渐减小，流量增加到1.4Qd时，低速区基本消失。由图11可知，在0.6Qd和1.0Qd工况下，中间流线流体的速度变化规律与压力面一致，说明叶轮流道中段亦存在低速区，但低速区强度和面积小于叶片压力面。由图12可知，在叶片吸力面附近，只有在小流量工况下会形成低速区。叶轮流道内，流体的相对速度总体随着流量的增大而增大。综上，当流量较小时，叶轮流道压力面及中间流道附近流体会形成低速区，随着流量增大，低速区的面积和强度会逐渐降低，中间流道内流体的低速区会逐渐消失。

3.3.2 额定流量下不同轴截面流体相对速度分布

本研究采用扭曲叶片叶轮，对前文提到的Z1、Z2和Z3截面进行PIV实验分析。由于不同轴向位置叶片的扭曲程度不同，因此对额定工况下L1流道内头发纤维质量浓度0.1%时，r0/r2=0.7、0.85及1的3个不同位置处从叶片吸力面到压力面流体的相对速度分布进行分析，结果如图13～图15所示。

图9 Z2截面处流体的速度分布Fig.9 Velocity distributions of fluids in Z2 cross-section

图10 不同流量下压力面流体速度分布Fig.10 Pressuresurfacevelocity distributionsof fluids at different flowrates

图11 不同流量下中间流线流体速度分布Fig.11 Intermediate flow channel velocity distri⁃butionsof fluidsat different flowrates

图12 不同流量下吸力面流体速度分布Fig.12 Suction surface velocity distributions of fluidsat different flowrates

对比图13～图15可知，在r0/r2=0.7处，叶片扭曲程度相对较大，因此各个截面上流体的相对速度相差较大且靠近压力面流体的速度变化幅度较大。而随着位置向出口方向移动，叶片的扭曲程度逐渐降低，各截面上流体的相对速度差值逐渐减小，到r0/r2=1即叶片出口位置处，叶片形状可近似于圆柱叶片，各截面流体的速度变化相对平滑，基本归于一致。

3.3.3 额定流量下不同头发纤维质量浓度Z2截面处流体的相对速度分布

为进一步研究纤维质量浓度对于纸浆泵内液相流场的影响，限于篇幅影响，本研究分别提取了Z2截面上L1流道内质量浓度0.1%和0.2%头发纤维悬浮液在0.6Qd、1.0Qd及1.4Qd工况和r0/r2=0.7及r0/r2=1情况下，沿吸力面至压力面的相对速度分布图（见图16和图17）。

图13 不同截面上r0/r2=0.7处流体的相对速度分布Fig.13 Relativevelocity distributionsof fluidsat r0/r2=0.7 and different cross-sections

图14 不同截面上r0/r2=0.85处流体的相对速度分布Fig.14 Relativevelocity distribution of fluidsatr0/r2=0.85 and different cross-sections

图15 不同截面上r0/r2=1处流体相对速度分布Fig.15 Relative velocity distributions of fluids atr0/r2=1 and different cross-sections

图16 r0/r2=0.7处不同流量时流体的相对速度分布Fig.16 Relative velocity distributions of fluids at r0/r2=0.75 and different flowretes

由图16及图17可知，随着流量逐渐增大，高质量浓度纤维悬浮液Z2截面内的相对速度会逐渐超过低质量浓度纤维悬浮液、各位置相对速度分布趋势相似度增大，且随着位置靠近叶轮出口区域，上述变化趋势加快。大流量工况下，各质量浓度纤维悬浮液在叶轮流道内从吸力面至压力面的相对速度总体呈下降趋势。且随着质量浓度增大，叶轮出口截面处流场的均匀性在各流量工况均有所改善，叶轮流道内各位置处流体的相对速度变化更趋平滑，有利于降低叶轮出口滑移率，提高叶轮效率。

图17 r0/r2=1处不同流量时流体的相对速度分布Fig.17 Relative velocity distributions of fluidsat r0/r2=1 and different flowrates

4 结论

本研究使用CFD-DEM耦合方法对离心泵内纤维颗粒悬浮液流场进行模拟和分析，发现其与PIV实验所得叶轮流道内速度场的分布趋势基本吻合，从而验证了CFD-DEM方法初步模拟研究纤维悬浮液内流场的可行性；进一步采用PIV实验结果对泵内流场进行具体分析，主要结论如下。

（1）随着泵内纤维质量浓度逐渐增大，泵效率随之提高，纤维质量浓度达到0.2%时已明显呈现出减阻效应。

（2）当输送相同质量浓度纤维时，大流量工况下叶轮流道内流体的流动稳定，未出现漩涡和低速区。而在小流量工况下，沿叶轮流道中段压力面附近流体会产生低速区，且低速区随着流量的减小逐渐向中间流道延伸。

（3）在垂直于泵轴的不同截面叶轮流道内，受叶片扭曲程度的影响，各截面处流体的相对速度存在差异。靠近轮毂侧，叶片扭曲程度较大时，各位置处流体相对速度分布相差较大。而随着位置向出口方向移动，叶片的扭曲程度逐渐降低，各截面上流体的相对速度差值逐渐减小，在靠近叶轮出口处各截面处流体相对速度变化相对平滑且基本归于一致。

（4）头发纤维质量浓度的变化会影响流道内流体的相对速度分布，当流量逐渐增大时，较高质量浓度纤维悬浮液内液相相对速度大小会逐渐超过低质量浓度纤维悬浮液，在大流量工况下，各质量浓度纤维悬浮液在叶轮流道内从吸力面至压力面的相对速度总体呈下降趋势。而随着纤维悬浮液质量浓度增大，叶轮出口截面流场的均匀性在各流量工况均有所改善，叶轮流道内各位置处流体的相对速度变化更趋平滑。上述结果可对理解长纤维悬浮液减阻效应提供参考。