小流量系数离心压缩机级的泄漏及轮阻损失的数值研究*

2022-12-21于学智闻苏平胡小文王志恒杜俊杰王懿

风机技术 2022年5期

于学智闻苏平胡小文王志恒杜俊杰王懿

（1.西安交通大学能源与动力工程学院；2.沈阳鼓风机集团股份有限公司；3.美的集团中央研究院）

0 引言

小流量系数离心叶轮一般采用闭式叶轮，闭式叶轮离心压缩机主要损失包括轮阻损失、泄漏损失和流动损失。在小流量系数离心式压缩机级中轮阻损失和泄漏损失所占比例尤为突出，在流量系数低于0.02时甚至超过输入功的50%。

离心压缩机小流量系数级基于叶片宽度的雷诺数Re（Re=U2b2/ν）是重要的特征参数，Re数的变化对级性能有显著的影响。Wiesner人[3]总结了在各种级的几何结构、转速及各种介质条件的雷诺数对离心压缩机级性能的影响。指出雷诺数和相对粗糙度的变化对级性能有明显的影响，并提出了改进的经验公式。1951 年Batchelor[4]对旋转圆盘间的不可压缩流体的流动特性进行了研究，求解了无限大旋转圆盘的微分方程组，得出当两个圆盘转动一致时，圆盘之间的核心流体以固定的周向速度做刚体运动，同时在两盘附近存在两个边界层，通常这种流动被称作Batchelor 型流动。1953 年Stewartson[5]采用数值方法研究了高速转动圆盘系统圆盘间的流动。通过求解流动控制方程发现旋转圆盘间隙内的边界层之间的流体切向速度都为零，通常这种流动称为Stewartson型流动。1991年Morse等人[6-7]提出了修正低雷诺数k-ω湍流模型，并进行了旋转圆盘内部流动的研究，获得了满意的预测结果。Owen等人[8]用数值方法求解旋转圆盘间的流动特性。Owen 研究发现，有过流情况的旋转圆盘间隙内，在较小的半径区域存在Stewartson型流动，并称这一区域为源区域，当旋转盘对流体的泵送能力增强，在较大的半径区域形成Batchelor型流动，并且根据湍流参数可以确定Batchelor型流动产生的半径区域。1984 年Chew[9-10]第一个用低雷诺数k-ε 模型数值研究了带离心过流的转静盘系统流动问题，研究结果不是很理想。主要原因是简单的使用对数壁面函数规律连接粘性亚层引起，且由于其边界层为三维湍流边界层，因此近壁面完全湍流区的速度方向和壁面剪切应力矢量方向相同的假设必然会引起较大的误差。随后又使用了混合长度模型对旋转圆盘内的流动进行计算，并和实验进行了对比，发现混合长度模型的预测精度较为准确。

国内也陆续开展了离心压缩机小流量系数基本级的研究。张勇等人[11]数值研究了叶片扩压器在小流量系数基本级中的应用，并和无叶扩压器结果进行对比。冀春俊等人[12]对一离心压缩机小流量系数基本级扩压器内部流动进行数值分析。马新民等人[13]对一小流量系数离心叶轮进行数值研究，分析了叶轮内部二次流动的发展过程。徐立群等人[14]数值研究了迷宫密封泄漏对小流量系数离心叶轮气动性能的影响。谭佳健等人[15]对一离心压缩机模型级轮盘侧和轮盖侧间隙的压力分布和转子的轴向力进行了测量，发现在叶轮轮盘和轮盖两侧间隙内压力沿着径向成线性分布，叶轮的轴向力是由静态力和动态力组成，叶轮轴向力随叶轮转速的增大而增大。王维民[16]通过数值计算方法建立了叶轮和密封间隙的模型，分析了轮盘外侧间隙流动和压力分布的关系，并在产品上进行了验证。王志恒等[17]数值研究了一小流量系数离心叶轮的叶轮两侧间隙流动，发现小流量系数离心叶轮性能受到密封泄漏和轮阻损失影响较大，指出在模拟时需要予以考虑。闻苏平[18]进行了旋转圆盘间气体流动数值模拟，采用不同湍流模型数值结果差异，并和实验对比，指出在旋转圆盘研究中，采用SST 模型和k-ω 模型进行计算，结果更为准确。因此本文采用k-ω 湍流模型对离心式压缩机基本级叶轮外侧间隙空腔及密封结构内的泄漏气体流动特征进行数值研究。

1 数值计算

1.1 几何模型

本文所研究的模型为小流量系数离心式压缩机基本级，其子午结构如图1 所示，叶轮半径R2=225mm，轮盖外侧密封前端半径R1in=124mm，轮盘外侧密封半径为R2in=98mm。密封1 的齿间距为5mm，齿数为5，密封间隙为0.3mm。密封2 齿间距为4mm，密封间隙为0.25mm，密封齿数为5。

图1 模型级子午几何结构Fig.1 Meridian geometry structure of model stage

1.2 网格划分

对叶轮的轮盘、轮盖外侧间隙及叶轮轮盖、轮盘处的密封均进行建模，采用分块结构化的方法对级生成高质量网格。为减小交界面处的插值误差，划分网格时尽量使交界面处网格均匀分布，并使交界面两边的网格在各个方向上网格分布基本一致。为了准确获取固体壁面的压力，在壁面处设置了边界层网格，以保证在计算过程中所有靠近壁面的第一层网格均落在粘性底层，即Y+值在10 以下。边界层第一层网格厚度为0.01mm，边界层网格最小正交性15.5，最大长宽比2053，网格节点总数6253159 个。流场计算网格如图2和图3所示。

图2 模型级单通道网格模型Fig.2 Single-channel grid model of model stage

图3 模型级网格子午投影图Fig.3 Grid meridian projection of Model stage

1.3 数值方法及边界条件

湍流模型选用k-ω 方程。为减小数值扩散引起的误差，密度、动量、湍动能、湍流耗散率、扩散项和能量对应的输运方程的对流项均采用二阶中心差分格式进行离散。采用空气作为介质进行分析，气体状态方程采用理想气体状态方程。

轮盖外侧结构进口按基本级叶轮出口的CFD相应计算结果给定静压，出口给定总温总压条件。轮盘外侧密封结构出口给定静压，进口给定总温总压。壁面边界条件设置为固定壁面、绝热、无滑移。当全局残差达到10-6时认为计算收敛。

2 泄漏损失计算

2.1 泄漏损失模型

Japikse[19]曾给出通过叶轮前后空腔和迷宫密封泄漏量的简单计算模型。无论是轮盘还是轮盖侧，叶轮出口和密封位置的速度和静压都可以写成相对总压力守恒表达式。相对总压力是可压缩流体转子焓的不可压缩当量，以轮盖侧密封为例，可得到如下方程式：

式（1）必须给定叶轮外侧空腔中的相对速度，包括叶顶处和密封位置处的相对速度。通常相对速度很难给定，因此只能使用近似表达式。每个位置的相对速度可以近似表示为当地叶轮圆周速度的分数（W=FU），可以假定叶顶和密封位置的相对速度也为当地叶轮圆周速度的一个分数。那么可得到如下表达式：

则密封前端的压力可表示成：

经过密封结构的静压差（P1-）由式（3）将以上公式密度平均值上面的斜杠简化掉。静压能通常可以使用密封结构的动能分数表示，则可以得到如下方程：

则有：

泄漏流量表示为：

因此泄漏量可表示成空腔F2因子和损失系数CS的函数。早期F2因子通常采用0.25，而Stepanoff[20]在水泵中使用F2为0.75，因此推荐将实际运行情况下测量通过密封结构前后气体的静压，通过这个压差可以推导出不同测量数据下的F2因子。在实际工业应用中，F2因子的范围在0.09～0.81之间。使用Japikse的这个模型可以用来粗略地估算泄漏损失，他的模型中还有一个损失系数CS需要估算。

对迷宫密封结构的泄漏流量计算，Egli[21]也提供了一个计算公式：

式中，A为密封齿与叶轮形成的通流面积；CD为流动系数，通常取0.67～0.95；n为密封齿数；Pin为密封前静压；Pout为密封后静压。

泄漏损失和密封齿数、密封间隙大小和密封结构的压差直接相关。而且当叶轮工作流量越小时，压差越大，而且泄漏量增大迅速，泄漏损失相应也越大。

假定忽略气流节流时的损失及齿间隙前后的可压缩性时，提供了一种关于密封泄漏量的计算公式，当密封间隙中的气流速度小于声速时，漏气量采用：

考虑气体膨胀时，由于摩擦引起的损失以及气流在间隙中所获得的动能，并不能在随后的空腔中完全转化为热能，而有部分转化成压能的因素。D为间隙的平均直径，s为间隙大小，分别为密封结构前的比容和静压，n为密封齿数，P1为密封结构后的静压。显然这个公式中密封前后的静压，密封前的压力也是需要根据经验公式得到的，可应用Japikse[19]的方法。

当密封间隙中的气流速度出现声速时，一般声速最先出现在第一个密封间隙中，泄漏量采用如下公式：

计算密封泄漏时，首先应该判断密封结构中的气流是否已经达到音速，然后选择不同的公式计算漏气量。可用如下临界压力公式来判断：

密封前后压力比为上式倒数，当压力比超过上式的倒数时，表面气流已经达到声速。

显然利用泄漏损失模型进行简单的泄漏量估算是均需要对一些系数进行确定，最直接的方式可以通过实验来确定，然而在实际情况下，这种实验难度较大，因此，本节使用数值计算来研究基本级的密封泄漏结构的流动问题。

2.2 数值计算结果与分析

在机器马赫数为0.6 时，叶轮在三个不同工作流量下的密封结构的数值计算获得的压力场如图4所示，这三个流量分别对应0.72Qdes、1.0Qdes、1.5Qdes。其中Qdes为叶轮设计流量。从图4 可以看出，在轮盖密封结构空腔中的静压从叶轮外径往内半径方向逐渐降低，压力梯度随叶轮流量增大略减小，且叶轮外侧空腔中静压沿轴向几乎没有压力梯度，压力等值线近似为水平。尽管从流量3 到流量1 的密封压差相差3 倍，空腔中的压力差却变化很小。因此，类似旋转圆盘结构的空腔具有稳定的维持从叶轮出口到密封结构前端截面的压降能力。

图4 不同工作流量工况下叶轮轮盖密封结构压力分布图Fig.4 Pressure distribution diagrams of impeller shroud seal structure under different working flow conditions

台阶迷宫密封结构内部的温度场分布如图5所示，当叶轮工作在大流量时，叶轮出口压力不高，导致叶轮进出口压差小，泄漏量会降低，叶轮轮盖外侧面与气流摩擦做功大部分转化成气体内能，因此，叶轮在大流量工况下台阶密封结构内部的温度高于设计流量和小流量工况的温度。由图5还可以看出，台阶密封前后的温度变化不大，略微有些升高，可近似为一个等温过程。台阶密封结构的流线分布在图6 中表示出来。在不同的叶轮工作流量下流体通过密封齿间隙后突扩形成一个顺时针方向的大涡和一个由于台阶导致的逆时针方向的小涡。通过这些涡结构使得压能在每个密封齿后面的小空腔中损失掉而转化成热量。

图5 不同叶轮流量工况下的台阶密封结构温度场Fig.5 Temperature fields of step seal structure under different impeller flow rate conditions

图6 不同叶轮流量工况下的台阶密封结构流线Fig.6 Streamlines of step seal structure under different impeller flow rate conditions

图7 给出了机器马赫数分别为0.6，0.7 和0.8，在叶轮设计流量工况下，三个不同相对半径位置r*（r*=r/r2，r为径向半径，r2为叶轮出口半径，在台阶密封结构前端位置对应的r*=0.551）的相对周向速度Vt*和相对径向速度Vr*分布情况。与闭式旋转圆盘系统相似，周向速度Vt*分布也表现出三层结构的流动型式。但是由于泄漏流量的存在使得径向速度Vr*在间隙内不存在径向速度为0的核心区。在机器马赫数分别为0.6，0.7和0.8下，密封空腔中的周向相对速度Vt*和径向相对速度Vr*分布型态几乎完全相同，说明对同一个密封空腔结构，周向速度Vt*和径向速度Vr*基本不随叶轮转速的变化而变化。

图7 设计流量下不同机器马赫数周向相对速度和径向相对速度分布Fig.7 Circumferential relative velocity and radial relative velocity distributions of different Mach numbers at design flow rate

图8是不同马赫数下F2因子随叶轮的流量工况的关系。图中显示对于具有相同密封空腔结构的叶轮，机器马赫数的变化对F2因子的影响很小，且叶轮流量的变化对F2因子的影响可以忽略，仅在流量非常大时，F2因子随流量增大略有增加。除在大流量下，所研究的流量系数为0.01028 的基本级叶轮轮盖的空腔F2因子的数值范围为0.74～0.76，与Stepanoff[20]给出的水泵中的F2=0.75 非常接近。因此在估算小流量系数离心叶轮轮盖密封结构前的压力时，选取F2=0.75是合理的。

图8 轮盖密封空腔中的F2因子Fig.8 The impeller shroud seal F2 factor of the cavity

图9 所示的是在机器马赫数为0.6 时，叶轮设计流量Qdes下轮盘侧密封结构中的压力分布，可以看到和轮盖密封结构的压力分布不同，在轮盘密封结构中，压力分布呈凹形，也就是说叶轮出口和回流器压力都较高，而轮盘外侧空腔底部位置压力最低。尽管泄漏方向是从回流器到叶轮出口方向，但是叶轮出口压力是高于轮盘空腔内半径处的压力的，并靠轮盘的旋转作用保持这个压力梯度。和轮盖外侧空腔相比还有一个不同是轮盘外侧空腔呈现大半径处的大压力梯度，而轮盖外侧空腔中的压力梯度几乎一致。从密封齿位置的压力分布放大图可以看到，密封齿的节流作用使得密封结构前后保持大的压力梯度，密封的作用是通过节流的过程达到，而密封齿之间的各个小空腔内压力变化很小。

图9 轮盘密封压力分布Fig.9 Disc seal pressure distributions

定义相对泄漏量为密封结构，实际泄漏质量与通过叶轮实际流量的比值。数值计算得到的轮盖外侧和轮盖外侧的相对泄漏量曲线如图10 所示。图中显示，随着叶轮流量的增加，相对泄漏量逐渐降低，这主要是由于叶轮流量增加以及叶轮出口压力随叶轮流量增加而逐渐降低等因素造成的。

图10 轮盖和轮盘侧相对泄漏量随叶轮流量变化曲线Fig.10 The relationship between the relative leakage rate of the shroud and the disc side and the impeller flow rate

机器马赫数对叶轮轮盘和轮盖外侧的泄漏量影响不大，随机器马赫数的增大略增加。在叶轮设计流量工况下，轮盖密封产生的泄漏量约为叶轮流量的4.5%，轮盘密封产生的泄漏量约为叶轮流量的2.7%。显然轮盖密封泄漏量占泄漏损失的主要部分。而轮盘密封的泄漏量同样不可忽视，相对轮盖密封结构的泄漏量而言，轮盘外侧泄漏量约为轮盖密封泄漏量的一半还多。

3 轮盖侧轮阻损失计算与分析

叶轮轮阻损失通过计算叶轮外侧面所受转动扭矩求得。可以把叶轮看成一个旋转圆盘系统。而轮盖侧的密封结构可以分成两部分：从密封结构前端到叶轮进口和从叶轮出口到密封结构前端。通常把后一部分结构等效成转静盘系统进行研究，其转静盘系统的旋转盘扭矩占据轮盖外侧扭矩的主要组成部分。

目前计算圆盘扭矩有两种办法，第一种通过分析转静盘流动，通过计算带动系数K并利用Haaser等人[22]的公式得到；第二种采用Daily 等人[23]的实验数据总结的经验公式得到。

转静盘系统的流动受间隙比（G=s/R）和旋转雷诺数Re=ωr2/v这两个几何参数以及流量、进口预旋等流动参数控制。假设G＜＜1且Re＞＞1，对应于旋转圆盘分区流动中的分离边界层湍流区Ⅳ（Batchelor 型流动）。则可以将轮盖外侧以及腔体分别对应于旋转圆盘的转盘与静盘，则摩擦扭矩可表示为：

切向摩擦应力定义为：

摩擦系数采用Haaser 等人[2-6]在文献中提出的通用表达式：

对转静盘系统有：

显然，如果驱动系数K为已知，那么轮盖外侧面的摩擦扭矩就可以计算出来，问题的关键在于带动系数K的准确计算。根据前人的研究，驱动系数K不仅是叶轮半径的函数，而且和叶轮流量和预旋值有关。

Poncet 等人[24]在文献中将流量Qt和带动系数K联系起来。根据文献[25]假设边界层的速度分布满足经典的1/7指数规律。那么摩擦系数Cf可以根据文献[26]Dean 方程给出：Cfmean=，其中Remean为基于Umean的雷诺数，此时Ekman 边界层的平均周向速度为Umean=(1+K)ωr/2，特征长度为Ekman边界层厚度δE，则摩擦应力为：

通过标识壁面上的剪应力形成的角度为α，剪应力的周向分量可表示为：

切应力的径向分量可以通过边界层控制容积的离心力和切应力平衡获得：

假设流线的角度α沿着半径基本保持不变，则对于Ekman边界层有：

然后假设Bödewadt边界层径向摩擦由径向流量QB决定，并假设在核心区的径向速度为0（这个假设对于Batchelor型流动是满足的），由连续性方程有：

式中，Qt为旋转圆盘的总流量；VE为Ekman 边界层的平均速度，和边界层的最大速度成正比（VE=βωr）。

使用以上方程可以得到：

定义一个新参数Cqr=，其中Rer=ωr2/ν为当地雷诺数，那么K的最终表达式为：

其中，a和b为实验确定的常数，b可以通过让系统的流量为0时的K值计算求得。然而a和b的值随进口预旋的不同而改变。

除了以上先通过计算带动系数K来计算圆盘扭矩的办法外，对有过流量的旋转圆盘系统Daily 等人[23]做过研究，实验拟合给出了在间隙比G分别为0.0276、0.069和0.124，湍流参数λturb为0～0.06之间时测量扭矩系数结果的经验公式：

其中，λturb=CwRe-0.8为湍流流动参数；K*为无过流时的带动系数；Cw=mleak/(μr2)为无量纲流量。本文的λturb为0.06，间隙比G为0.0178，然而Daily 等的公式未考虑流动预旋的影响，因此本文对考虑进口预旋时的计算扭矩系数进行了研究，并提出计算扭矩系数的公式。

4 预旋对扭矩系数的分析

轮盖外侧密封空腔流动时，必须考虑进口流体预旋对轮盖空腔内流动的影响。预旋Kp定义为Kp=Vt2/ωr2，其中Vt2为叶轮出口周向速度。图11 给出了机器马赫数为0.6时，不同进口预旋对轮盖外侧空腔内的流动结构的影响，为显示方便将径向尺寸进行压缩。可以看到所有进入轮盖外侧空腔的流体很快进入静盘侧（左边）边界层，随后部分流体离开静盘侧边界层，穿过核心区进入转盘侧边界层，部分流体则直接进入齿形密封结构。在空腔内部则形成旋涡，随着进口预旋的增大，空腔内旋涡中心渐渐远离进口区。当无进口预旋时，在空腔内形成一个旋涡，而当进口预旋增加到1时，空腔顶部又出现一个小旋涡。

图11 不同预旋对轮盖密封空腔内流线的影响Fig.11 The influence of different prewhirl on the streamline in the shroud seal cavity

由预旋引起的这种流动形式的改变必将引起空腔内部的带动系数K发生变化。图12给出了机器马赫数为0.6时，在叶轮设计流量下的不同轮盖密封进口气流预旋对轮盖空腔中的带动系数K随Cqr的变化曲线。图中可以看到进口预旋对带动系数K的影响，在接近轮缘处（Cqr=0.004），带动系数随进口预旋增大急剧地增大。且随半径减小，带动系数K随进口预旋量的增加而增加的速度越来越接近。在Ponect 等人[24]的研究中给出了他们的实验关系式为K=-1。

本文也采用和Ponect 等人[24]相似的关系式K=-1 来拟合图12 中进口预旋为0 时的情况，得到下式：

拟合的结果为：a=5.32，b=0.63。拟合的b值和Ponect 等人[24]的实验结果完全一致，拟合的a偏小。产生差异的主要原因是不同的几何结构以及气体的可压缩性造成的。注意到本文拟合的公式（23）在Cqr=0 也就是说无系统过流量时，可以得到K=0.438，Owen[8]总结层流解的K=0.382，Daily[23]的实验数据为0.42～0.475。

尽管在预旋为0的情况下，K值可以以Cqr的5/7次方的规律表达出来，但是从图12中可以看到，随着预旋的增加，在大半径处（对应于小的Cqr）越来越偏离以上拟合方程（23）。因此在大预旋的情况下，K值用式（23）计算则误差较大，必然对计算轮盖摩擦扭矩带来较大的误差，而实际叶轮轮盘侧密封进口预旋并不为0。

图12 不同进口预旋时K与Cqr的变化关系Fig.12 The relationship between K and Cqr with different inlet prewhirl

为评价进口预旋对轮盖外侧摩擦扭矩及扭矩系数的影响，图13 给出了在叶轮设计流量下三个机器马赫数Mu的摩擦扭矩及扭矩系数随预旋的变化关系。扭矩系数定义为：

图13 轮盖外侧扭矩与扭矩系数随进口预旋的变化Fig.13 The torque and torque coefficient of the outer side of the impeller shroud vary with the Prewhirl of the inlet

式中，A为轮盖外侧面面积。

从图13可以看出预旋对轮盖扭矩以及扭矩系数CM产生显著的影响，在进口预旋为0.25附近，在所有机器马赫数下的扭矩系数CM均出现最大值，随着进口预旋的进一步增加叶轮轮盖侧扭矩系数CM逐渐降低。轮盖扭矩随转速的增大而增大，扭矩系数随转速的增大（因为旋转雷诺数Re也增大）而减小。

Daily[23]公式考虑了过流量、间隙比和雷诺数变化对圆盘扭矩系数CM的影响，但并未考虑进口预旋对扭矩系数的影响。这里根据扭矩系数的计算结果来拟合扭矩和进口预旋的关系，并修正经验公式得到计算轮盖外侧轮阻扭矩系数的计算公式如下：

其中，Kp为轮盖密封进口的预旋值系数。

图14出了式（25）得出的扭矩系数与进口预旋的关系曲线，点数据为数值计算的扭矩系数CM，曲线为式（25）的结果。通过计算叶轮外侧面的扭矩并利用公式（24）计算求得。在三个机器马赫数下，与数值计算的比较可以看出，Daily 的公式考虑加入进口预旋后的结果与数值计算结果符合较好，误差最大出现在预旋为0.2附近，当进口预旋大于0.7时，误差很小，而对实际叶轮轮盖的进口预旋通常接近0.7。因此，采用本文修正的方程可以较准确地计算出轮盖外侧的扭矩系数CM。