基于非定常空化流动的离心泵涡旋结构数值分析

2023-02-01司乔瑞

西华大学学报（自然科学版） 2023年1期

伍杰，邱宁，朱涵，徐佩，司乔瑞

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

离心泵是一种重要的流体机械[1]，内部流动复杂且不稳定，时常伴随着流动分离、失速及回流等现象的发生[2−4]，同时泵内压力低于汽化压力还会造成空化的发生[5−6]。因此，泵内的复杂流动结构还需要进一步研究探讨[7]。

空化是水力机械中一种常见的现象。在离心泵当中，积聚的空泡能够堵塞叶轮流道，造成泵的性能下降。Medvitz 等[8]模拟了离心泵的空化流动，发现快速水头下降现象的物理机制是气泡爆炸性增长，堵塞了一半以上的叶轮流道。空泡的溃灭还会对叶片表面形成冲击，造成表面材料的剥蚀。Qiu 等[9]以水翼的研究载体，对空泡破灭时产生的能量冲击进行了模拟预测，发现最大能量冲击位于空腔闭合线附近。

在传统的流场研究中，研究人员以压力变化作为内部流场分析的重要参数。为了更加深入地了解流场特征，基于速度变化的涡动力学从提出之后便被广泛应用[10−11]。李等[12]将实验与模拟相结合，运用正则化螺旋度法，完成了混流泵启动时瞬时流场中涡核的提取，发现涡旋结构会随着转速的增加而出现正反交替现象，但在转速稳定之后流动也逐渐平稳。Lin 等[13]研究了高压涡轮的尾迹涡和非定常流动，对周期性尾迹涡和涡轮损失进行了捕捉与分析。考虑到流量大小对于流动状态的影响，Zhang 等[14]基于数值模拟方法研究了低比转速泵的复杂流动，将3 种不同流量工况方案对比发现，随着流量的降低，泵内流动变得紊乱，在小流量工况条件下已形成明显的涡流结构。敏等[15]也得出相似结论，流道内的涡量随着流量的增加而逐渐减小，并且随着介质流动，旋涡发生破裂并逐渐减小。Posa 等[16]对比研究了在设计工况与非设计工况下带有不同扩压器的离心泵性能，发现流动分离和回流现象是由流经叶轮的压力梯度所决定的。由于非定常流场比较复杂，为了揭示流动原理，Zhang 等[17]研究了离心泵内的非定常流场，将压力脉动与涡旋结构演化相结合，证明了叶片尾缘脱落涡与压力脉动存在明显的相关性。曹等[18]研究发现局部压力较低的吸力面脱落的旋涡逐渐向上游回流，并发展成为独立的分离涡。Sun 等[19]研究发现旋转方向不同的涡旋对会导致空腔周期性脱落，空化的发展能够增加涡旋强度，而增强的涡旋又反过来促进空腔的脱落。

随着计算流体力学的发展，有学者发现以往的涡流提取方法存在一些局限性。Liu 等[20−22]基于涡旋的物理特征，将其分为旋转与非旋转部分，提出了新Omega 涡识别方法。该方法能够很好地避免传统涡识别方法在阈值选取的不确定性。本文基于密度修正模型（DCM）对RNGk−ε湍流模型进行修正，模拟了额定流量、不同空化条件下的非定常流场结构；基于平面速度分量，分析二维涡旋结构的演变以及与空化的相互作用；采用新Omega涡识别方法对三维涡旋结构与空化进行耦合分析。

1 离心泵参数和实验装置

1.1 离心泵参数

本文以IS65-50-174 型单级单吸式离心泵为研究载体，该测试泵的计算域三维结构如图1 所示，设计参数如表1 所示。

表1 测试泵设计参数Tab.1 Design parameters of test pump

图1 离心泵计算域三维结构Fig.1 3D calculation domain of centrifugal pump

1.2 实验装置

本次实验在江苏大学流体机械工程技术研究中心进行，实验台结构如图2 所示，测试设备参数如表2 所示。测试系统主要包括水循环回路和数据采集系统两部分。其中，通过调节测试泵出口管路阀门来控制流量，并基于流量计进行流量数据采集，在进出口管路布置压力传感器以监测压力。

图2 实验台结构Fig.2 Experimental bench structure

表2 实验设备参数Tab.2 Parameters of experimental equipment

2 数值计算方法

2.1 网格划分

图3 所示为采用ANSYS-ICEM 软件对叶轮流域的网格划分。与四面体的非结构网格相比，六面体的结构网格具有更少的网格数量，可以提高计算效率。又由于结构网格更为规整，网格控制也更为容易，尤其是叶片壁面附近网格。因此，为了保证数值计算精度的同时还能提高计算效率，测试泵采用六面体结构化网格划分方法。

图3 叶轮网格Fig.3 Impeller grid and details

考虑到网格对计算的影响，现选取额定工况进行网格无关性验证，并以扬程的变化作为网格适用性的评判标准。网格无关性验证如图4 所示。根据网格验证计算结果，最终选择网格数为2.54×106的方案作为计算所用网格。无关性验证后的最优网格方案可以在数值计算时降低网格因素对模拟计算结果准确性的干扰，同时还能避免由过多网格数目造成计算资源的浪费。为能够更好地捕捉到叶片壁面附近的流场变化，需要保证第一层网格高度在适当范围之内。本文将第一层网格高度设为0.042 mm，网格生长比为1.5，叶片壁面Y+数值小于100，达到计算要求，其分布如图5 所示。

图4 网格无关性验证Fig.4 Grid-independent verification

图5 叶片壁面Y+分布Fig.5 Yplus distribution of blade wall

2.2 数值计算

在数值计算中，进口边界条件为压力进口，其值为101325 Pa，流量出口数值为50 m3/h。叶轮与蜗壳计算域的壁面粗糙度设为0.03 mm，其余壁面均设为无滑移壁面。叶轮设置为旋转域，转速2900 r/min，其他部分设为静止域。将蜗壳、进口延长段分别与叶轮相接的面设为动静交接面，其余交接面以静静交接的方式进行耦合。

由于非定常流动状态较为复杂，为获得准确的内部流道模拟结果，本文选用RNGk−ε湍流模型进行模拟计算。该湍流模型可以很好地捕捉高速流动介质的流动状态，提高了旋涡流动的模拟精度。同时，为更加精确地获得离心泵内的空化发展情况，考虑气相和液相混合的可压缩性，使用DCM 方法进行混合密度修正。湍流黏度定义如下：

式中，N根据文献推荐取10[23]。

3 结果分析

3.1 实验验证与定常分析

本文采用实验与模拟对比验证测试泵在额定转速下不同流量工况的扬程，如图6 所示。实验与模拟均呈现出随着流量的增大而扬程逐渐降低的趋势。在额定工况时，实验测得扬程为35.71 m，模拟计算扬程为36.39 m，相对误差为1.9%，计算结果较为可靠，可以进行进一步的空化计算。

图6 不同流量下的扬程曲线Fig.6 Flow rate-head curve

图7 所示为额定工况下离心泵模拟与实验的空化性能曲线对比。在实验中，保持流量不变，使用真空泵以抽真空的方式降低进口压力，从而获得不同进口压力下的汽蚀余量。从图7 中可以看出，在空化初始阶段，a点的扬程略有上升，但并不影响计算结果，此时也不会有空化现象的发生，这一点在传统意义上被定义为非空化条件。随着进口压力的降低，b点的扬程开始出现下降趋势，将该点称为初始空化条件。当进口压力降低到c点所对应的压力时，离心泵已经出现了较为明显的扬程下降，叶轮当中也会出现明显的空化。随着进口压力的进一步降低，扬程出现急剧下降的趋势，当扬程下降3%，即空化性能曲线上的d点，普遍认为此时已经处于临界空化状态，并将该点称为临界空化条件。当进口压力低于d点时，扬程陡降，并出现严重空化，此时也已失去工程意义。

图7 空化性能曲线Fig.7 Cavitation performance curve

3.2 非定常空化性能分析

若流体的流动状态不随时间发生变化，这种流动状态称为定常流动。随着时间不断变化的非定常流动更贴近自然界中流体的流动状态。在非定常计算时，将计算总时长设为0.4138 s，即叶轮旋转运动20 个周期所用时间，时间步长设为2.298 9×10−4s，即叶轮每旋转4°进行一次计算，每个时间步长内迭代1～10 次，尽可能达到收敛精度，同时叶轮每旋转20°保存一个结果文件。将数值计算的收敛残差精度设为1×10−4，以此保证计算结果的可靠性。

图8 所示为典型时刻3 个不同空化条件下叶轮流道中截面液相速度流线–气相体积分数分布情况。在NPSHa=3.59 m 条件下，流道截面上的空腔覆盖面积较小，液相速度流线分布更为均匀合理，而随着汽蚀余量的逐渐降低，空腔覆盖面积逐渐增大，流道内的流线也逐渐紊乱起来。尤其是，当NPSHa=2.125 m（临界空化条件）时，由于流道内的空腔尾部迅速加厚卷曲，挤压流道空间，造成流道堵塞，从而在卷曲的空腔之后区域形成一个由于液体介质回流而产生的局部旋涡，表现出液相速度流线逐渐汇聚，形成一个或多个流线团。对比不同空化条件下的截面流线，可以看出空腔在较大程度上影响着流道内的介质流动状态。发生较为严重的空化时，叶轮流道下游区域形成较大程度的回流，直接造成泵的性能下降。为进一步分析NPSHa=2.125 m 时的空腔分布情况，图9 所示为典型时刻不同流道截面气相体积分数分布。可以发现，更为靠近前盖板的Span=0.8 截面上空腔覆盖面积最大，说明空化最为严重的区域位于前盖板附近，且空腔覆盖面积向后盖板方向逐渐减小。在Span=0.2 截面，空腔主要分布于叶片吸力面前端，很好地附着于叶片表面，而在吸力面中后部区域还存在已完全卷曲但体积分数较小的空腔。相较于Span=0.2 截面，在Span=0.5 截面上叶片表面空腔附着区域沿着流向方向有较大程度的扩大，在叶片中上游区域空腔很好地附着于吸力面表面，而空腔尾部迅速扩展变厚，并在叶轮旋转效应作用之下空腔尾部卷起脱离叶片表面，造成流道堵塞，严重影响液体介质流动，造成泵的性能下降。在Span=0.8 截面，后端得到充分发展的空腔开始附着于邻近叶片压力面，至此，同时连接着吸力面与压力面的空腔已完全将流道堵住，流动性能完全恶化。

图8 不同汽蚀余量叶轮流道中截面液相速度流线–气相体积分数分布Fig.8 Distribution of liquid velocity streamline-vapor volume fraction in the middle section of impeller passage with different NPSHa

图9 不同叶轮流道截面气相体积分数分布（NPSHa=2.125 m）Fig.9 Distribution of vapor volume fraction in different impeller passage sections（NPSHa=2.125 m）

3.3 二维涡旋结构与演变分析

在高速旋转的离心泵中，内部复杂流态监测捕捉较为困难，通常依靠数值模拟的手段完成内部流场分析。涡量作为一个显著的湍流特征，分析涡量的分布与演变可以更好地理解泵内湍流的非定常流场。式(3)给出了y方向的涡量定义[14]。

式中，vx与vz分别为x、z方向的速度分量。由于本文分析的泵模型绕y轴旋转，故采用x、z方向的速度分量来进行涡量的定义。

图10 所示为典型时刻不同流道截面气相体积分数–涡量分布。同图9 进行对比可知，图10 白色标注的区域1 为空泡覆盖区域，而在其后的蓝色标注区域2 为涡量分布区域。从图中能明显地看到，更为靠近后盖板的Span=0.2 截面中，在空腔卷起之后区域存在着许多高强度的正值涡量分布。这是因为该截面中，空腔覆盖区域较小，流道内的扰动较小，流态相对稳定，所以截面内的高强度涡量更为集中连贯。在Span=0.5 截面中，由于空腔卷曲更为严重，流道内的扰动大，流态不够稳定，集中连贯的高强度涡量破裂成为一些面积较小的分散涡量片。在Span=0.8 截面中，卷曲的空腔已附着于相邻叶片的压力面上，截面内的流道已被完全堵住，空腔之后的区域流动极不稳定，高强度涡量进一步分散。附着于相邻叶片压力面的空腔之后区域，即黑色标注3 所在区域也分布着面积可观的高强度涡量，该涡量还有向上游区域进一步发展的趋势。

图10 不同叶轮流道截面气相体积分数-涡量分布（NPSHa=2.125 m）Fig.10 Distribution of vapor volume fraction-vorticity in different impeller passage sections（NPSHa=2.125 m）

图11 所示为非定常流动状态下，单枚叶片旋转1/6 个周期，NPSHa=2.125 m 时离心泵叶轮与蜗壳4 个不同时刻的涡量分布。以图11（a）中灰白色叶片所在位置t=0 时刻，且叶轮旋转20°为一个时刻进行截面分析。从图中可以明显地看到，叶轮流道内部分布着许多速度为正的正向旋转涡量，而叶片尾缘为反向旋转的负值涡量。

图11 NPSHa=2.125 m 叶轮-蜗壳截面涡量分布及演变Fig.11 Distribution and evolution of vorticity in the impellervolute section of NPSHa=2.125 m

在不同时刻，同蜗壳流道相交接的叶片尾缘附着有尺寸不一的尾迹涡。由于叶片尾缘与蜗壳隔舌的动静干涉作用，叶片尾迹涡在隔舌附近变化较大。当t=0 时，灰色叶片正好旋转到隔舌前端，可以很明显地看到，在动静干涉作用之下，该枚叶片尾缘处的大尺寸负值尾迹涡开始发生脱落。在t=1/18T时刻，灰色叶片已完全扫过蜗壳隔舌，尾迹涡脱落完成，并附着于隔舌前端的蜗壳壁面。当t=2/18T时，灰色叶片尾迹涡重新附着生长，而隔舌前端涡量也基本消散，在其之后的叶片也携带着大量尾迹涡向隔舌靠近。在t=3/18T时刻，灰色叶片之后的叶片已旋转到隔舌前端，并开始重复t=0时刻在动静干涉作用下而发生的尾迹涡脱落现象。

对于4 个不同时刻，叶片吸力面中上游区域存在一个紧贴叶片壁面的负值涡量带。其后区域还有着正值涡量的分布，且该正值涡量带很快便发生卷曲而脱离叶片表面，并随着叶轮的旋转开始出现涡带脱落，然后在流道内逐渐消散。

观察图12、13 可以发现，从演变起始时刻t=0开始，流道内的涡量随着叶轮的旋转而呈现出逐渐减少的趋势，这与NPSHa=2.125 m 的涡量演变相对应。同时，有一个比较有意思的现象是尾迹涡附近还存在有些许正值高强度涡，并随着叶片的旋转而逐渐生长发展，最后脱落附着在蜗壳壁面。将3 个不同的空化条件对比可知，空化程度较小的工况，流道内的高强度涡量分布更少，流态更为简单合理。

图12 NPSHa=2.624 m 叶轮-蜗壳截面涡量分布及演变Fig.12 Distribution and evolution of vorticity in the impellervolute section of NPSHa=2.624 m

图13 NPSHa=3.59 m 叶轮-蜗壳截面涡量分布及演变Fig.13 Distribution and evolution of vorticity in the impellervolute section of NPSHa=3.59 m

图14 给出了NPSHa=2.125 m 时，一个旋转周期内典型叶片不同展长基于 Ωy的涡量分布。从图中可以看到，在叶片起始部分均存在一个涡量骤变的过程，并且不同展长位置的起始涡量还存在一定差异。这是由于来流首先冲击叶片前缘，从而在叶片前缘区域形成局部涡旋，而span=0.8 相对于span=0.2 与span=0.5 受到的冲击更大，所以在span=0.8上形成最大起始涡量。叶片吸力面前半段有高强度涡量分布，说明在流道上游区域流动变化较大，并且流动情况复杂，同时空化也容易在该区域内生长并发展。

图14 NPSHa=2.125 m 一个旋转周期内叶片吸力面上涡量分布Fig.14 Distribution of vortex on the suction surface during one rotation cycle with NPSHa=2.125 m

随着时间的演变，叶片前半段的涡量出现正负交替，形成一个周期性的变化过程，而后半段的涡量数值在小范围内变化。由于空腔在叶片中段位置卷起，而涡旋随着卷曲的空腔脱离叶片表面，在叶片中段位置呈现出骤变的过程。在L=0.3～0.6 区间内，不同展长的空腔体积存在较大差异，所以在此区域的涡量出现较为明显的波动。

3.4 三维空化– 涡结构分析

涡是一种流体在旋转过程中产生的物理现象，是一个具有方向与大小的矢量，但以往常用的识别方法得到的涡是一个只有大小而没有方向的标量，不能完整地表现出涡的具体参数与物理特征。Liu 等[20]提出了Omega 涡识别方法，主要是将流场中的涡分解为旋转部分与非旋转部分来进行涡特征的体现。

式中：R代表旋转部分涡量；S代表非旋转部分涡量，在一般情况下R与S具有不同方向。

为了更好地表征具有旋转效应的涡量场，将表征旋转涡量与流场总涡量之比定义为 Ω，而该定义参数 Ω在0～1 范围内取值。可以看出，随着数值的增大，所表征的旋转涡具有更高的旋转效应，计算公式如下：

为更好地识别流场中的涡旋结构，选取反对称张量B大于对称张量A，即 Ω>0.5 时作为流场涡识别的判据条件。同时，Liu 等[20]在文献中给出了Ω=0.52的参考值来进行识别涡旋结构边界。与该Omega 涡识别方法相比，Q准则等一些常用的涡识别方法具有更大的取值范围，合适的涡识别参数选取较为不易。

图15 所示为NPSHa=2.125 m 时叶轮流道内气相体积分数为0.1 的空腔与Ω=0.52的涡旋结构分布。从图中可以看到，在流道上游区域，涡旋结构较为稳定，且被空腔覆盖区域很好地包裹起来。由于临界空化状态的空腔结构比较稳定，空腔内部流动扰动较小，从而空腔内部形成的涡旋与空腔有着相似的稳定结构。在叶轮旋转效应的作用下，叶片中部位置的空腔与涡旋结构均卷曲脱离叶片表面。在叶轮流道下游区域，由于没有了稳定的空腔作为涡旋的外部覆盖结构，涡旋结构变得不再稳定，从而形成许多随着介质流动而不断变化的小体积涡。

图15 NPSHa=2.125 m 空腔-涡旋结构分布Fig.15 Distribution of cavity-vortex structure with NPSHa=2.125 m

三维空腔与涡旋结构的提取，能够更好地对上述二维涡量的演变做出进一步的解释。在空化的作用影响之下，叶片吸力面的涡量表现出同空腔相似的变化规律，沿着叶片弦长方向，吸力面的涡量大小随着表面附着空腔的增厚而逐渐变大。在叶片中段位置，由于空腔卷曲引起涡旋结构脱离叶片表面，从而产生叶片表面涡量骤变的现象。由于流道下游区域的涡旋结构零散分布于其中，并没有完整地附着于叶片表面，因而在吸力面后半段的涡量没有太大波动，数值较小。

图16—17 分别对应NPSHa 为2.624 m 及3.59 m时的空腔– 涡旋结构。不同于临界空化工况，其余空化条件下涡旋结构较为稳定集中，在流道下游区域没有出现零散分布的现象。随着NPSHa 的增加，空腔体积减小，没有占据过多的叶轮流道空间，扰动减小，流态更为稳定且合理，流道内的涡旋结构也随之逐渐减少。