双叶片泵内压力脉动的数值模拟

0 引 言

双叶片泵主要用于输送含有固体物质的液体，在冶金、矿山、化工等工业和河道疏浚、污水处理等部门中有着广泛的应用，深入研究其内部流动对于提高效率和改善其性能具有重要的意义。

国内学者主要针对清水离心泵及轴流泵内的压力脉动数值计算进行研究[1-5]。其中文献[1]采用“冻结转子法”处理离心泵内叶轮与蜗壳间动静耦合流动的参数传递和相互干扰问题问题，得出了离心泵内压力分布及速度分布,文献[2]对离心泵内叶轮与蜗壳动静耦合诱发压力脉动频谱进行了试验研究；模拟方面，文献[3]对5个叶片的普通清水离心泵泵内叶轮与蜗壳动静耦合进行了研究，分析了在设计流量点下，泵内蜗壳流道及叶轮流道内部监测点压力波动及压力脉动频谱，而文献[4]及文献[5]则分别对于高速泵及轴流泵内叶轮与导叶内部流场及其压力脉动进行了分析，得出压力脉动频谱趋势与试验值基本吻合。国外在泵内压力脉动的研究方面起步较早，已有较深入的研究[6-9]，1982年 Dring等人[6]提出导致泵内压力脉动 2种因素，分别为动叶轮与蜗壳之间的相互作用及泵内冲击。文献[7-9]对清水离心泵内压力脉动进行了深入研究，并得出一致的结论，即离心泵内压力脉动的主频等于叶轮转频与叶片数的乘积。

本考研究对象为双叶片特殊形式叶轮的离心泵，采用RNGκ-ε模型，在对蜗壳内压力脉动波动及其频谱数值模拟进行研究的同时，针对泵内叶轮进口及蜗壳出口受脉动的影响进行了分析，揭示了叶轮进口、叶轮与蜗壳耦合处及蜗壳出口压力脉动的主要影响因素及压力脉动产生的机理，对于掌握双叶片泵内压力脉动产生的原理及其特性具有重要的指导意义。

1 叶轮基本参数和计算域

叶轮的设计参数：流量Q=250m3/h，扬程H=25m，转速1 470 r/min，比转速126，设计模型为离心式双叶片泵。采用PRO/E软件生成三维计算区域模型，如图1所示，为使模拟结果更加稳定，对叶轮进口进行了适当的延伸，整个模型由1个动叶轮水体、1个静止蜗壳水体及进口延伸段水体组成。模型中未考虑容积损失及叶轮前后盖板与流体的摩擦机械损失。

图1 计算区域图
Fig.1 Calculation domain

2 网格的划分及监测点

2.1 网格划分

采用FLUENT中自带的网格划分软件GAMBIT进行网格划分，其中网格类型采用的是四面体非结构化网格，叶轮的网格数为581 014，蜗壳网格数为371 367，通过网格检查，叶轮和蜗壳计算区域的网格等角斜率和等尺寸斜率均不超过0.83，符合网格的等角斜率和等尺寸斜率不超过0.85的质量要求，三维网格计算区域图如图2所示。

图2 网格划分图
Fig.2 Sketchmap of mesh

2.2 监测点定位

图3为模拟中监测点示意图，图中除进口监测点7位于叶轮进口与进口延伸段的交界面中心位置上（参照图1），其他监测点位于叶轮的中间截面，点1为隔舌处监测点，点2～5为叶轮出口每隔90°取一监测点，蜗壳出口处取监测点6，其中监测点1完成叶轮进口压力脉动的监测，点1～5主要监测叶轮蜗壳耦合处压力脉动，监测点6可以实现泵出口流动受压力脉动的影响。上述监测点可完整地实现对泵内流动状态的监测。

图3 泵内各监测点示意图
Fig.3 Indicator points in the pump

3 数值求解方法

3.1 边界条件

3.1.1 进口边界条件

计算域进口边界条件采用的是速度边界条件，由质量守恒和进口无旋假设得进口速度即为轴向速度 win，切向速度及径向速度为0。其进口湍动能k可按以下公式计算：

式中：k——湍动能，J；win——轴向速度，m/s。

湍动能耗散率ε可按下式计算：

式中：l——进口混合长度，m；Cμ——经验常数，其值取0.09[10]。

3.1.2 出口边界条件

在计算域的出口边界上速度及压力未知，采用自然出流边界条件。

3.1.3 壁面条件

壁面用于限定fluid和solid区域，在本考模拟中，固体的壁面采用无滑移边界条件。根据壁面相邻的单元区的运动状态，定义壁面的运动状态，运动速度的大小采用相对相邻单元区速度（relative to adjacent cellzone）来定义。

3.2 数学模型

本考模拟为在定常计算基础上进行非定常计算，定常模拟采用多重坐标系下的 MRF求解，动量方程和连续性方程采用的是 SIMPLEC算法联立求解，计算求解模型采用的是分离求解器（segregated solver）。非定常计算时叶轮与蜗壳之间的耦合，也即动静耦合，采用了适用于非定常计算的滑移网格技术，动量方程和连续性方程采用的是PISO算法联立求解。

大涡模拟是介于直接数值模拟和湍流模式之间的一种直接数值模拟方法，目前已有部分泵内部流动研究采用大涡模拟的方法，在分析泵内压力脉动方面大涡模拟也已得到运用[11-13]，将随着计算机的发展成为数值模拟的一种发展趋势，但同时相比于大涡模拟，雷诺时均方法更加成熟，对网格、计算机要求低，同时计算时间短，在诸多文献[3-5]中已证实可以模拟泵内的压力脉动，其中RNGκ-ε模型在形式上类似于标准κ-ε模型，但是在计算功能上强于标准 κ-ε模型[14]，综合考虑本考计算采用雷诺时均方法，湍流模型为RNGκ-ε模型。

4 结果与分析

4.1 外特征预测

图4为根据模拟计算结果分析所得预测值[15]及试验得到测量值的扬程-流量对比图，模拟中扬程为根据模拟中的瞬时扬程计算出的有效扬程，虚线代表模拟值，实线代表试验值。由图4可以看出，预测值与试验结果值在趋势上完全一致，与试验值相比在设计流量及大流量时模拟值吻合的较好，小流量时模拟值小于试验值，但不超过 4%。由此可见，模拟中采用的数学模型比较准确地预测泵的外特性，保证了以下进一步分析的正确性。

图4 流量扬程曲线
Fig.4 Curves of head and flow rate

4.2 隔舌处绝对速度分布

图5为3种不同流量下在定常计算时隔舌附近的绝对速度矢量分布图。

蜗壳内产生压力脉动，其主要原因是叶轮与蜗壳存在的动静干涉，特别是叶轮与隔舌所导致流动的不对称性[3]。因此，在分析非定常前对定常计算的不同流量隔舌附近的流动状态对压力脉动产生有重要意义。如图5所示，在小流量时，在蜗壳扩散管区域的速度较小，而在隔舌靠近叶轮区域速度大于其他流量下的速度；大流量时，在隔舌下端存在明显的低速区同时有一部分流体流向了叶轮内，而靠近蜗壳第1断面处一部分流体流向隔舌，这2股流体在隔舌处相互碰撞，在蜗壳第1断面附近形成低速漩涡；设计流量下，流动较为稳定，速度分布均匀。从图5中同时能看出，叶轮出口速度在遇到隔舌时其速度矢量方向产生突变，这种突变必然在隔舌处产生一组方向相反的正负压力波，其中一压力波直接传向蜗壳出口，另一压力波传向蜗壳内；而叶轮出口速度分布并不均匀，各流道中水流流动与蜗壳的相对位置的不同而变化[1]，因此这种压力波的频率为叶轮的转频与叶片数的乘积，可以预测蜗壳与叶轮动静耦合压力脉动主频与隔舌处产生的压力波频率有关，在叶轮转动方向上这种压力波将逐渐衰落。

图5 不同流量下隔舌附近速度分布
Fig.5 Absolute velocities close to tongueat different flow rate

注：Q——设计流量，下同

4.3 叶轮蜗壳耦合处压力脉动分析

图6为叶轮与蜗壳耦合处1～5监测点在流量Q、0.5Q及1.5Q下的时域图，图7为叶轮与蜗壳耦合处1～5监测点在流量Q、0.5Q及1.5Q下的频谱图。

非定常的计算在定常计算的基础上进行，图6中不同流量下各监测点的时域图为非定常计算收敛后再进行四周期的计算所得到的数据值。其中横坐标为时间，纵坐标为静压值。图7中的频谱图由收敛后四周期内监测点的监测值运用 Origin8.0软件通过快速傅里叶变换所得，其中横坐标为频率，纵坐标为幅度值，幅值的大小代表受其对应的频率影响的程度，幅值越大影响越大，图中 1～5代表各监测点。

图6 不同流量下各检测点时域图
Fig.6 Time-domain spectrums of indicator pointsat different flow rate

图7 不同流量下各检测点频谱图
Fig.7 Frequency spectrums of indicator pointsat different flow rate

参照图 6中各监测点在不同流量下的压力脉动时域图，从图中可以看出各监测点压力波形都具有明显的周期性，在小流量时各监测点的压力平均值最小，大流量时压力平均值最大，这与离心泵流量与扬程关系相符。不同流量下隔舌处监测点1的变化幅度要明显大于其他监测点。图中各监测点波形存在相位差，其中监测点2与监测点4的相位一致，监测点3与监测点5的相位一致，监测点1相位与监测点2相差不多，导致相位差主要由非定常计算所导致，在非定常计算中周期由式（3）计算。

式中：T——周期，s；z——叶片数；n——泵转速,r/min。每一周期叶轮旋转的角度为180°，分100时间步长计算，即每一时间步长为周期的1%。

以图6中时域图监测点2与监测点3为例，监测点2与监测点 3空间上相差90°（见图3），由非定常中计算周期式（3）计算得叶轮旋转180°时间为0.020408 s，则监测点2与监测点3相差90°需要时间差为0.010204 s，等于时域图中相位差。

频谱图具有2个重要意义，一个是从时域的波形中分解出各个频率所含有震动成分大小，从中可以获得波形的卓越震动频率；另一个重要的意义是将时域的波形转换到频域，或将频域的信号返回到时域。因此对于本考可以从时域图中分析出频谱图，利用频谱图分析出造成泵内压力波动的卓越频率，从而找出压力波动主要影响因数，卓越频率为幅值最大点对应的频率。参照图 7中各监测点在不同流量下压力脉动频谱图，在不同流量下各监测点主要受低频的影响，而受高频区影响很小与文献[7-8]结果一致；各监测点的卓越频率即主频都在50 Hz左右，这刚好等于叶轮转频与叶片数的乘积（叶轮转频为叶轮每秒旋转的次数，单位 Hz），与 4.2节中的预测相符。对比图7中不同流量下的频谱图，同一监测点在小流量时主频的幅值要大于大流量与设计流量下的幅值，而在设计流量Q时具有最小的幅值，说明在设计流量点下泵内产生的压力脉动最小。同一流量下，从隔舌附近监测点1到监测点5主频的幅值呈递减的趋势，反映了在隔舌处产生的压力脉动波随着叶轮转动方向在逐渐减小，在监测点5时达到最小值，以图7 a所示主频的幅值从点1的62 000Pa2/Hz减小到点5时的30 000Pa2/Hz，同样与 4.2节预测一致，压力脉动在隔舌处开始沿着叶轮旋转方向上逐渐衰落。

4.4 进出口压力脉动分析

图8为叶轮进口在不同流量下的压力脉动时域图及频谱图，图9为蜗壳出口在不同流量下的压力脉动时域图及频谱图。

叶轮进口流动在泵内产生压力脉动之前，但其流动状态受叶轮蜗壳动静耦合作用的影响，通过数值模拟发现叶轮进口也存在压力脉动，但与叶轮与蜗壳耦合处压力脉动存在差异。如图8所示，从时域图中可以看出，进口监测点压力波的周期为其他监测点的2倍，从频谱图中可以看出，在不同流量下监测点的主频都为25 Hz左右，等于叶轮的转频，且其主频的幅值（1 200 Pa2/Hz）相对于叶轮与蜗壳耦合处脉动主频的幅值（35 000 Pa2/Hz）很小，因此叶轮进口流动相对较为稳定，且其压力波动主要受叶轮的转频脉动的影响，叶轮转频造成的脉动主要原因为叶轮流道壁厚的不均，这同时也反映进口主频幅值小的原因。图8中压力脉动时域图中看出不同流量下进口静压变化较大，由于不同流量下的模拟考察压力相同，而大流量时叶轮进口流速相对较大，因此其静压值小。

图9为蜗壳出口监测点在不同流量下的时域图及频谱图。从图中时域图可知，在大流量时静压值要小于小流量时的静压值，不同流量下具有相同周期的正弦波；同时由于蜗壳出口受泵内叶轮与蜗壳动静耦合作用特别在隔舌附近对流体进入蜗壳扩散段的影响，在频谱图上表现为与叶轮与蜗壳耦合处具有相同的主频（50 Hz），同样在小流量时主频的幅值最大，而在设计流量点时具有最小的幅值，这同样表明在蜗壳出口压力脉动为泵内产生压力脉动的延续，具有相同的主频，但其强度已减弱。

图8 不同流量进口监测点时域图及频谱图
Fig.8 Time-domain spectrums and frequency spectrums of indicator points of inletat different flow rate

图9 不同流量蜗壳出口监测点时域图及频谱图
Fig.9 Time-domain spectrums and frequency spectrums of indicator points of volute casing outletat different flow rate

5 结论

通过对双叶片泵进行非定常计算，分析不同监测点在不同工况下的压力时域图，压力波呈正弦周期性变化，进口监测点压力波的周期为其他监测点的2倍；同时对不同监测点不同工况下频谱图分析，泵内蜗壳内压力脉动产生的主要原因为叶轮与蜗壳的动静耦合，特别与叶轮出口速度在通过隔舌时矢量方向突变产生的压力波有关，且这种压力波在叶轮旋转方向上逐渐衰退导致压力脉动沿蜗壳圆周方向上的变化逐渐减弱；在不同工况下同一监测点的主频不变，但在小流量及大流量时主频的幅值增大；叶轮进口处压力脉动主频等于叶轮转频，叶轮与蜗壳耦合处及蜗壳出口压力脉动主频相等，为叶轮转频与叶片数的乘积，即叶片的通过频率，对于掌握双叶片泵内压力脉动产生的原理及其特性具有重要的指导意义。