双吸离心泵叶轮内泥沙磨损非定常特性研究

张自超 李延频 陈德新

(1.华北水利水电大学电力学院，郑州 450045；2.华北水利水电大学河南省流体机械工程技术研究中心，郑州 450045；3.华北水利水电大学乌拉尔学院，郑州 450045)

0 引言

广泛应用于黄河沿岸灌溉泵站中的双吸离心泵，各部件存在明显的磨损问题。叶轮作为高速旋转的部件，其磨损问题更为严重，对性能影响更显著。因此，研究双吸离心泵叶轮的磨损特性具有重要的意义。

目前，数值模拟是研究水力机械泥沙磨损的主要手段。较多学者采用定常计算研究水力机械的磨损特性[1-8]。然而，泵内存在明显的动静干涉作用，必然导致固液两相流的内流场具有非定常特性。如压力脉动就是动静干涉引起的，这说明动静干涉对流场的影响不可忽视。

越来越多的学者关注了泵内固液两相流的非定常特性[9-13]。然而，水泵的磨损与固相速度、固相浓度和固相与壁面的碰撞角等流场参数有关，因此，水泵的磨损也必然具有非定常特性。

叶轮作为泵内受动静干涉作用影响最大的部件，其磨损的非定常特性必然显著。而现有研究关注较多的是扬程、功率、压力脉动等的非定常特性，并未对磨损的非定常特性进行过多研究，尤其是对双吸离心泵叶轮内的泥沙磨损非定常特性的相关研究较少。本文以双吸离心泵为研究对象，采用欧拉-欧拉方法对其内部固液两相流进行非定常计算，研究不同工况下，叶轮壁面上多个监测点的固相速度、固相浓度和磨损率的非定常特性，并与定常计算结果进行对比。

1 欧拉-欧拉方法

1.1 控制方程

欧拉-欧拉方法的控制方程为

(1)

(2)

其中

αl+αs=1

式中k——相类别(l为液相，s为固相)

α——体积分数μt——湍动粘度

v——时均速度p——压强

Fi——相间作用力gi——体积力

κ——泥沙扩散系数ρ——密度

t——时间μk——动力粘度

λ——应力xi、xj——坐标分量

vki、vkj——速度梯度

其中，相间阻力的计算体现在式(2)的相间作用力项Fi中，而含泥沙扩散系数项是根据梯度输移假定，通过泥沙扩散系数建立各相浓度、速度脉动关联项与浓度梯度的关系得到的，具体可参考文献[14]。

1.2 磨损模型

图尔萨大学的磨损模型是目前泥沙两相流磨损计算的经典模型[15]，考虑了颗粒冲击速度、冲击角和壁面材料的硬度对磨损的影响，即

(3)

其中

F(θ)=5.4θ-10.11θ2+10.93θ3-6.33θ4+1.42θ5

(4)

ER=WRmsand/Acell

(5)

式中WR——磨损失重率，壁面磨损失去的质量与颗粒质量的比值，kg/kg

BH——壁面材料的布氏硬度

Fs——颗粒形状系数

vs——颗粒入射速度，m/s

θ——颗粒冲击角，rad

a、C——经验常数，取2.41、2.17×10-7

ER——磨损率，kg/(m2·s)

msand——固相质量流量，kg/s

Acell——计算单元壁面网格面积

棱角颗粒Fs=1，半球形颗粒Fs=0.53，球形颗粒Fs=0.2，本文取Fs=0.2。

2 计算对象和设置

2.1 计算对象

本文的研究对象为山西尊村引黄泵站的双吸离心泵，主要几何参数：额定流量Qn为3 m3/s，额定转速n为490 r/min，设计扬程H为32 m，比转数ns为162，叶轮直径D为1.12 m，叶片数z为6，泵进口直径D1为1.1 m，泵出口直径D2为0.9 m。

计算含沙条件为：含沙质量浓度5 kg/m3、体积分数0.188%时，颗粒密度2 650 kg/m3；沙粒粒径ds为25 μm的悬移质典型粒径。

2.2 网格及边界条件

双吸离心泵的计算网格由Gambit生成，经过网格无关性检查，最终确定各部分网格单元数分别为：叶轮603 620；吸水室564 882；压水室400 447；叶轮前后腔水体64 288，口环间隙较小，计算复杂，而且对叶片磨损影响不大，因此本文不考虑口环间隙；进口延长段44 100；出口延长段56 202；总计网格数为1 733 539。网格无关性检查如图1所示，双吸离心泵计算域及网格划分情况如图2所示。

图1 网格无关性检查Fig.1 Grid-independent verification

图2 计算域及网格划分Fig.2 Computational domain and grid

应用Fluent进行数值模拟，运用欧拉-欧拉方法对双吸泵进行固液两相流数值计算。采用Phase Coupled SIMPLE算法求解二阶迎风格式的离散方程。湍流模型采用RNGk-ε模型，因为该模型可以更好地处理高应变率及流线弯曲程度较大的流动[16]；考虑到相间阻力和滑移速度的影响，湍流多相流模型采用Dispersed turbulence模型，该模型是目前应用较多的一种多相流湍流模型[17-18]；泥沙扩散系数模型采用Diffusion-in-VOF模型进行计算。壁面磨损计算采用图尔萨大学的磨损模型[15]。

计算域进口采用速度进口，分别给定固、液两相的速度，并给定固相体积分数；出口采用自由出流；过流部件内壁面，不考虑粗糙度，设为光滑壁面，对液相采用无滑移壁面边界条件，对固相采用自由滑移壁面边界条件，近壁区采用标准壁面函数。

以定常计算的结果为初始条件进行非定常计算，共计算10个叶轮旋转周期，每个旋转周期计算120个时间步，即每旋转3°为一个时间步。

2.3 阻力模型选择

相间阻力模型是固液两相流计算中重要的相间作用因素，合理地选择阻力模型可以提高固液两相流的计算精度。

常用的固液相间阻力模型有Syamlal-O’Brien模型[19]、Wen-Yu模型[20]、Gidaspow模型[21]和Huilin-Gidaspow模型[22]等。

Syamlal-O’Brien模型是基于流化床中颗粒的沉降得到的，因为模型中考虑了固体剪切粘度参数，随着固相浓度趋于零，固相剪切粘度也趋于零[23],并不适用于稀疏固液两相流的计算[24]。Wen-Yu模型被推荐用于固相体积分数小于20%的情况[17,25]，而Gidaspow模型和Huilin-Gidaspow模型是Wen-Yu模型和Ergun模型的组合，这两种模型可以用于更高固相浓度的两相流计算。但是，本文含沙体积分数远小于20%，此时，Wen-Yu模型、Gidaspow模型和Huilin-Gidaspow模型的本质均为Wen-Yu模型，尤其是在固相体积分数小于1%的情况下。因此，本文选用Wen-Yu模型作为相间阻力模型进行计算。

2.4 叶片监测点布置

为了监测叶轮内部固相参数的脉动情况，需在叶片壁面上设置监测点，各监测点位置如图3所示。

图3 叶片上监测点设置示意图Fig.3 Monitoring locations on blade

分别在叶片吸力面和压力面上各设置3个监测点，沿叶片吸力面中间流线从进口到出口分别设置3个监测点，分别为：S1位于叶片吸力面进口处，S2位于叶片吸力面中间位置，S3位于叶片吸力面出口处，同样对应压力面设置3个监测点P1、P2、P3。6个监测点位于叶片上并随叶片一起旋转。

各监测点监测的固相参数为固相速度、固相浓度和监测点处的壁面磨损率和冲击角，各监测值只记录参数的大小，不监测方向。

3 数值计算验证

3.1 外特性验证

为了验证数值计算的合理性，将数值计算得到的双吸离心泵扬程和效率与试验值作对比，并分析误差，如图4所示。计算工况分别为0.4Qn、0.6Qn、0.8Qn、Qn、1.2Qn。

图4 双吸离心泵外特性计算值和试验值对比Fig.4 Comparison of calculated value and experimental value for external characteristics of double suction centrifugal pump

由图4可知，计算得到的效率和扬程与试验值均较为接近，但有一定的误差。扬程的计算值高于试验值，最大误差为0.8Qn处的3.6%；而效率计算值较接近试验值，最大误差为Qn处的2.1%。因此，外特性计算值与试验值较吻合，计算误差在可接受的范围内，计算结果合理、可靠，可用于进一步的分析和研究。

3.2 固液两相流固相体积分数分布验证

为了进一步对固液两相流流场进行验证，选取文献[26]中的一个圆管为计算对象，如图5所示。采用与上述双吸离心泵完全相同的计算模型进行固液两相流计算，将计算得到的固相体积分数分布与试验值进行对比验证。

图5 计算圆管示意图Fig.5 Sketch of calculated circular pipe

圆管直径Dp为0.15 m，长度L为3 m。进口流速vl为3 m/s，进口平均含沙体积分数αvf为15%。颗粒粒径ds为90 μm。

圆管水平布置，坐标原点位于圆管底部，y轴竖直向上，线1为验证固相体积分数分布的位置，距离进口3L/4且竖直向上。

计算结果如图6所示。图中y*=y/Dp，αs为圆管中线1上各点的固相体积分数，αvf为圆管进口平均固相体积分数。

图6 固相体积分数计算值和试验值对比Fig.6 Calculated solid concentration distribution compared with experimental data

由图6可知，本文采用的两相流计算模型得到的固相体积分数的计算值与试验值分布趋势一致，结果较为吻合，相差不大，最大误差为9.5%，计算误差在可接受的范围内。由此证明，本文采用的固液两相流计算模型合理、可靠，可用于流场的进一步计算和分析。

4 结果分析

4.1 叶轮磨损特性

4.1.1叶片上磨损率脉动特性分析

图7给出了各监测点不同工况时磨损率非定常和定常计算结果对比，图中，T为叶轮旋转周期。

由图7可知，叶轮表面各监测点的磨损率具有周期性，其波动周期等于叶轮旋转周期，这与动静干涉引起的叶轮内固液两相流流场具有脉动特性有关，这说明叶轮壁面的磨损具有脉动特性。同时，最大磨损率与最小磨损率之间相差较大，脉动特性明显。随着流量的增大，磨损率逐渐增大，脉动更加剧烈。在叶片吸力面，叶片进口处的磨损率较小，出口处的磨损率较大；而在叶片压力面，叶片中部的磨损率较小，最大磨损率在叶片出口处。叶片压力面出口处为叶轮磨损最大位置，而压力面中部的磨损率最小，这可能与该处流态变化有关。

图7 各监测点不同工况时磨损率非定常和定常结果对比Fig.7 Comparisons for unsteady and steady results of erosion rate under different operation conditions at each monitoring point

由图7可以看出，定常计算得到的磨损率远小于非定常的结果，为了更清楚地对比二者的计算结果，表1给出了各监测点不同工况时定常计算得到磨损率和非定常计算得到的最大磨损率的对比。由表1可以看出，定常计算得到的磨损率远小于非定常的结果，这说明采用定常方法计算双吸离心泵的壁面磨损率时，误差较大。而随着流量的增大，定常和非定常计算得到磨损率的误差增大。这说明定常计算并不能准确反映叶轮壁面上磨损的脉动特性。

表1 定常和非定常计算磨损率对比Tab.1 Comparison of unsteady and steady results of erosion rate kg/(m2·s)

4.1.2叶片上磨损分布特性分析

为了分析叶片表面上的磨损分布情况，图8给出了不同工况下，叶片吸力面和压力面的磨损云图。

为了对计算得到的磨损分布进行验证，图9给出了离心泵叶片表面的实际磨损情况，其中图9a为文献[27-28]得到的离心泵叶片表面磨损形貌试验结果示意图，图9b为双吸离心泵叶片现场磨损图。由于本文研究的双吸离心泵缺少磨损试验，图9b为同类型其它的双吸离心泵的实际磨损情况，根据磨损规律，其磨损情况可以定性地验证本文的数值计算得到的磨损分布情况。

由图8可知，叶片表面的磨损率随着流量的增大而增大，磨损面积也逐渐增大。各工况下，磨损最大的部位为头部和尾部，这与图7中非定常计算结果一致。吸力面中部的磨损面积大于压力面，而在头部和尾部，压力面的磨损大于吸力面，这与图9a的试验结果一致。同时根据图7得到的结论和图8的分布情况，可以看出，压力面尾部磨损最为严重，这与图9b中的双吸离心泵叶片现场磨损情况一致。同时也说明本文的计算结果可靠。

图8 不同工况下叶片表面的磨损分布云图Fig.8 Calculated erosion distributions of blade surface under different operation conditions

图9 离心泵叶片实际磨损情况Fig.9 Actual erosion conditions of centrifugal pump blade

4.2 固液两相流场特性分析

由式(3)～(5)可知，对于确定的计算对象和条件，壁面材料硬度、颗粒形状系数等均为定值，影响磨损率的主要流场参数为冲击角、固相体积分数和固相速度。为了进一步分析叶轮壁面磨损特性及两相流场特性，对各监测点在不同工况时的冲击角、固相体积分数和固相速度的脉动特性进行分析。

4.2.1冲击角及冲击角函数特性

冲击角对磨损率的影响，通过冲击角函数F(θ)来反映，冲击角函数F(θ)的具体形式见式(4)。图10为壁面冲击角的示意图，冲击角为固相速度与壁面切线的夹角，其大小为0°～90°，由图可知，固相颗粒在撞击壁面之后反射射出。

图10 壁面冲击角示意图Fig.10 Diagram of wall impact angle

图11为冲击角函数随冲击角变化规律，由图可知，随着冲击角的增大，冲击角函数先增大后减小，其最大值出现在50°附近，这说明存在一个使磨损最大的冲击角，此时的切削作用最显著。

图11 冲击角函数变化曲线Fig.11 Law of impact angle function

当冲击角为20°～80°时，冲击角函数大于1，说明冲击角增强了磨损损失；当冲击角小于20°或大于80°时，冲击角函数小于1，说明冲击角减弱了磨损损失；当冲击角达到50°附近时，冲击角函数最大，冲击角对磨损损失的影响最大。

图12给出了各监测点不同工况时，冲击角及其函数的波动规律。由图12可知，冲击角和冲击角函数随叶轮旋转呈周期性变化，其脉动曲线与图7中磨损率的脉动曲线较相似，这说明冲击角对磨损的影响至关重要，其决定了磨损率的脉动特性，以及最大磨损率和最小磨损率。冲击角具有增强或减弱磨损的作用。同一监测点处，不同工况下，冲击角大小有所不同，这与工况改变后，流场的流态有所改变，导致冲击角发生变化有关。

图12 各监测点不同工况时冲击角及其函数波动规律对比Fig.12 Comparisons for fluctuation of impact angle and its function under different operation conditions at each monitoring point

由此可知，冲击角对磨损的影响为增强或减弱磨损效果的作用，其决定了磨损率的脉动特性及最大磨损率、最小磨损率。

4.2.2固相体积分数特性

图13给出了各监测点不同工况时，固相体积分数的变化规律。由图13可知，固相体积分数随着叶轮旋转也具有周期性。同一监测点处，随着流量的增大，固相体积分数有所增加；而随着流量减小，固相体积分数随时间的波动更为显著，其原因可能是小流量时，偏离最优工况，流态变差，二次流和漩涡增多，流场波动明显。相同工况下，吸力面和压力面对应监测点的固相体积分数分布相差不大。吸力面和压力面从进口到出口，固相体积分数分布波动逐渐减小，这是因为叶片进口处的流态较为复杂，容易出现二次流和漩涡，在叶片出口处，流态基本较为稳定。

图13 各监测点不同工况时固相体积分数波动规律对比Fig.13 Comparison for solid volume fraction under different operation conditions at each monitoring point

通过与图7对比，发现固相体积分数的脉动曲线与磨损率脉动曲线差别较大，这说明固相体积分数并不能决定磨损率的脉动特性，只是对其大小有所影响。

4.2.3固相速度特性

图14给出了各监测点不同工况时，固相速度的变化规律。由图可知，随着叶轮旋转，固相速度的变化具有周期性。同一监测点处，随着流量增大，固相速度逐渐增大。通过对比图7可以发现，最大固相速度出现在叶片压力面出口处(P3)，最小速度位于压力面中间(P2)，这与图7中磨损率最大和最小的位置完全相同，这说明图7中压力面中间(P2)处出现磨损率最小主要是由该处的固相速度最小引起的，而叶片压力面出口(P3)处磨损率最大主要是由该处固相速度最大引起的。

图14 各监测点不同工况时固相速度波动规律对比Fig.14 Comparisons for solid velocity under different operation conditions at each monitoring point

对比图14和图7可以看出，磨损率与固相速度的脉动规律非常吻合，这说明固相速度对磨损大小的影响更为显著。

5 结论

(1)叶轮表面磨损率分布具有周期性，其脉动周期为叶轮旋转周期，其脉动特性明显。随着流量增大，磨损率增大，脉动特性更明显。叶片压力面出口处的磨损率最大。

(2)定常计算的磨损率远小于非定常结果，随着流量增大，二者误差增大。因此，定常计算结果不能准确预测叶轮壁面上的磨损损失。

(3)磨损最大的部位为头部和尾部。吸力面中部的磨损面积大于压力面，而在头部和尾部，压力面的磨损大于吸力面。

(4)叶轮壁面冲击角、冲击角函数、固相体积分数和固相速度分布均具有周期性。冲击角函数的脉动曲线与磨损率脉动曲线较相似，冲击角对磨损损失具有增强或减弱的作用，存在一个使磨损率最大的冲击角，其决定了磨损率的脉动特性。

(5)在小流量时，固相体积分数脉动明显。固相体积分数的脉动曲线与磨损率脉动曲线差别较大，其对磨损率脉动特性影响较小，对磨损大小有所影响。固相速度脉动规律及固相速度最大、最小值所处位置均与磨损率脉动特性相似，固相速度对磨损率影响显著。