基于微元轨迹的密封动力特性系数理论识别方法

顾乾磊， 张万福， 张 尧， 陈璐琪， 李 春， 杨建刚

(1. 上海理工大学 能源与动力工程学院,上海 200093；2. 东南大学 火电机组振动国家工程研究中心,南京 210096)

密封是汽轮机、燃气轮机、压缩机等大型旋转机械中抑制流体泄漏的关键部件，直接影响机组的高效运行[1]。然而，随着机组容量及工质参数的不断提高，密封引起的气流激振问题越来越突出，对高参数大型透平机械安全与稳定运行形成极大挑战[2-5]。近年，对于密封的研究重点也因此由传统的静力特性(如泄漏性能)逐渐转向动力特性。与滑动轴承类似，通常用8个或12个动力特性系数(刚度、阻尼及质量惯性项)来描述密封的动力特性，但密封间隙通常是轴承间隙的2～10倍，密封内部流动为更加复杂的三维可压缩湍流。因此，密封动力特性系数的高精度识别及气流激振的准确定位是目前国内外研究的热点与难点。

目前对于密封动力特性系数的识别方法主要有理论、试验及计算流体力学(Computational Fluid Dynamics, CFD)。密封理论计算在早期的Bulk-Flow模型[6-7]之后，主要以控制体模型为主，学者先后提出单控制体模型[8-9]、双控制体模型[10-13]、三控制体模型[14-15]，取得了一定的预测效果，然而控制体模型在处理偏心、腔室周向速度、能量损耗系数、边界条件等方面做了大量假设，特别是随先进密封几何结构复杂性与工质参数的提高，以及阻塞流动带来的较差敛散性，限制了其应用的广度和深度。试验是密封动静特性研究的重要手段[16-17]，主要方法有动态压力测试法和机械阻抗法，后者由于具有较好的信噪比，在实际中被广泛应用。Texas A & M大学的Childs教授和Vance教授等在实验方面做了大量研究工作，从20世纪80年代开始至今一直从事密封(如梳齿、蜂窝、孔型、刷式、袋式等)动静特性的试验识别。1986年第一次给出了密封直接阻尼系数(最高压力8.25 bar、最高转速8 000 r/min的试验值[18]；1988年对试验台进行了改进(最高转速16 000 r/min)，并与控制体模型进行了对比，发现双控制体模型的预测结果要优于单控制体模型[19]；2002年Dawson等[20-21]在已有轴承试验台基础上对其改进，密封试验最高压力达17.2 bar，2003年试验参数进一步提高至最高压力84 bar、最高转速29 000 r/min。Tiwari[22]对密封试验的研究进展做了很好的综述。国内关于密封动力特性的试验研究相对较少，杨建刚等[23-24]提出密封动力特性系数识别的影响系数法和基于不平衡同频激励的密封动力特性系数识别方法。然而，开展试验研究，特别是不同类型、多种尺寸与工况下的密封性能测试需要大量时间和人力物力投入，代价较高，而且试验在流场细节等方面的测试较为困难。

随着计算机硬件的高速发展及其计算能力的不断提高，计算流体力学方法正在受到越来越多的应用。从早期基于有限差分、有限体积等方法的小规模编程计算[25-26]，到目前以大型商业计算软件(如FLUENT、CFX、STAR CD等)为主的CFD三维湍流模拟。密封动力特性系数的CFD识别方法主要有如下三种：

(1)有限小扰动法。该方法与滑动轴承动力特性系数的识别类似，通过在平衡点附近给转子正交方向上较小的位移扰动和速度扰动，进而得到刚度和阻尼系数，如式(1)所示

(1)

(2)旋转坐标系法[27-28]。如图1所示，该方法应用旋转坐标系将转子实际的非稳态涡动问题转变为准稳态模型，计算速度较快。

瞬态涡动法。该方法基于CFD非稳态求解及动网格方法，对密封内流场进行计算。西安交通大学李军教授团队做了大量的密封数值模拟工作，提出了多频涡动瞬态计算方法[29-31]，并针对蜂窝密封、孔型密封、袋式密封等开展了系统的数值仿真研究。

图1 准稳态模型示意图Fig.1 Schematic diagram of quasi-steady state model

有限小扰动法和旋转坐标系方法使用较为方便，然而在实际应用过程中局限性较多，如：① 转子实际涡动轨迹较为复杂，如存在偏心、椭圆涡动等；② 密封动力特性系数与偏心率、涡动频率、转速等都存在依赖关系；③ 复杂密封壁面(蜂窝、孔型、袋式等)的模拟较难实现。瞬态涡动法避免了上述问题，通过在不同方向施加激振计算密封气流力及密封动力特性系数，然而针对转子实际任意轨迹的计算还需要更深入地研究。

本文提出基于微元轨迹的密封动力特性系数理论识别方法，建立任一频率激励下流体动力学识别模型，应用CFD瞬态求解及动网格方法得到密封气流力特性及动力特性系数，并与试验值进行对比，验证模型与计算方法的准确性。

1 密封动力特性系数识别方法

图2 气缸-密封系统模型Fig.2 Cylinder-seal system model

图2中：1为转子涡动轨迹；2为涡动转子；3为静子(气缸)；4为静平衡位置的转子；Fη,Fξ为转子在η和ξ方向所受气流力；Fe,Fα为转子在e和α方向所受气流力。

转子做椭圆轨迹涡动，涡动转速为Ωi(i=1, 2, 3, …,n)，(X,Y)坐标系下的涡动轨迹为

(2)

式中:θ为椭圆轨迹倾斜角θ∈(0°,360°)；a和b分别为椭圆轨迹的长、短半轴长度；X0和Y0分别为涡动轨迹中心在(X,Y)坐标系下的横纵坐标。

在(e,α)坐标系中，椭圆轨迹方程为

(3)

相应的转子涡动速度为

(4)

以O1点为原点建立坐标系，忽略气体质量惯性力影响，转子小涡动轨迹下密封动力学模型可线性简化为

(5)

将式(3)与式(4)代入式(5)得

(6)

应用瞬态法得到任意涡动转速Ωi(i=1, 2, 3, 4, …,n)转子所受气流力，在t=0和t=T/4时刻，转子在e方向和α方向气流力为Fe(t=0, Ω=Ωi),Fα(t=0,Ω=Ωi),Fe(t=T/4,Ω=Ωi),Fα(t=T/4,Ω=Ωi)。

将t=0，t=T/4代入式(6)可得

ΔFe(t=0,Ω=Ωi)=-a·Kee-b·Ceα·Ωi

(7)

ΔFα(t=0,Ω=Ωi)=-a·Kαe-b·Cαα·Ωi

(8)

ΔFe(t=T/4,Ω=Ωi)=-b·Keα+a·Cee·Ωi

(9)

ΔFα(t=T/4,Ω=Ωi)=-b·Kαα+a·Cαe·Ωi

(10)

以涡动频率为横坐标，ΔFe(t=0,Ω=Ωi),ΔFα(t=0,Ω=Ωi),ΔFe(t=T/4,Ω=Ωi)和ΔFα(t=T/4,Ω=Ωi)为纵坐标，如图3所示绘制气流力之差随涡动频率变化图，分别记作曲线Ω-ΔFe(t=0,Ω=Ωi),Ω-ΔFα(t=0,Ω=Ωi),Ω-ΔFe(t=T/4,Ω=Ωi),Ω-ΔFα(t=T/4,Ω=Ωi)。

对于任意涡动频率Ωi，取微元轨迹(Ωi,Ωi+ΔΩ)，

图3 气流力之差随涡动频率变化Fig.3 Variation of seal force difference with whirling frequency

当ΔΩ→0时，微元段内动力特性系数不随涡动转速变化而发生改变，即当涡动转速为Ωi和Ωi+ΔΩ时密封动力特性系数Kee,Keα,Kαe,Kαα,Cee,Ceα,Cαe,Cαα不发生改变。将涡动转速为Ωi时对应的作用在密封段转子上的力ΔFe(t=0,Ω=Ωi)和涡动转速为Ωi+ΔΩ时对应的作用在密封段转子上的力ΔFe(t=0,Ω=Ωi+ΔΩ)代入式(7)得

ΔFe(t=0,Ω=Ωi)=-a·Kee-b·Ceα·Ωi

(11)

ΔFe(t=0,Ω=Ωi+ΔΩ)=-a·Kee-b·Ceα·(Ωi+ΔΩ)

(12)

对微元段ΔΩ取极值得

(13)

由式(13)得拟合曲线在横坐标为Ωi处切线的斜率。因此过曲线Ω-ΔFe(t=0,Ω=Ωi)上任意点作切线，切线斜率为式(7)中的系数-bCeα，解得Ceα即切点横坐标对应的涡动转速下的交叉阻尼系数，切线在纵坐标上的截距即为式(7)中的常数项-aKee，解得Kee即切点横坐标对应的涡动转速下的直接刚度系数。同理可得不同涡动转速下的动力特性系数Kαe,Cαα,Keα,Cee,Kαα,Cαe。

由于(e,α)坐标比(η,ξ)坐标超前一个θ角，两坐标系关系为

(14)

(15)

将式(15)用矩阵形式表示

(16)

由图2可知，用矩阵的形式表达ΔFe,ΔFα,ΔFη,ΔFξ之间的关系

(17)

易得

(18)

将式(18)用矩阵形式表示

(19)

将式(16)与式(19)合并得刚度系数为

(20)

同理，阻尼系数为

(21)

2 数值计算方法

2.1 几何模型

本文数值模拟的密封几何尺寸为Ertas试验数据。表1给出了密封的具体尺寸，密封齿处细节尺寸如图4所示。

表1 密封几何尺寸

(a)整体模型

(b)局部尺寸图4 密封几何模型Fig.4 Seal geometry

2.2 数值模型

网格划分是CFD计算中的重要环节，网格质量和布置的合理性对计算精确性和收敛速度具有很大影响。为了提高计算精度，采用了结构化网格，并对流动变化剧烈的齿顶处进行适当加密。本文采用壁面函数法将壁面的物理量与湍流核心区联系，并将y+值控制在30～150区间内，密封流场网格如图5所示，网格节点数共计114万个。

图5 密封流场网格划分Fig.5 Grid distribution

2.3 计算工况

表2为本文数值计算的参数细节。工质为理想空气，采用k-ε湍流模型，湍流强度5%，转子和静子壁面设置为绝热光滑无滑移壁面。在密封进口设置总温、总压，出口设置平均静压出口。瞬态计算的时间步长依据涡动频率进行调整，本文将转子涡动轨迹简化为圆轨迹，涡动中心为密封几何中心。在计算模拟时采用了动网格方法，运用CEL设定网格运动轨迹，网格的运动由转子的涡动方程决定。当转子上偏心和垂直于偏心方向受力呈现周期性变化、不同周期内对应点受力值相差小于0.01 N、进出口流量相差小于0.1%时，认为计算收敛。

表2 工况参数

3 计算结果与讨论

将数值模拟得到的直接刚度系数Kavg、直接阻尼系数Cavg、交叉刚度系数Kxy和有效阻尼系数Ceff与Ertas等的试验结果进行对比分析。其中Cavg,Kavg,Ceff的定义为

(22)

式中:Cxx和Cyy分别为x与y方向上的直接阻尼系数；Kxx和Kyy分别为x与y方向上的直接刚度系数。

3.1 密封动力特性系数

图6给出了直接刚度系数随涡动频率的变化趋势。模拟所得直接刚度系数Kavg在数值上小于试验值，特别是在低频段，符号呈相反趋势，即密封表现为负刚度，且符号不随涡动频率变化而改变，该结论与文献[32]一致。随涡动频率的增加，二者逐渐吻合，刚度系数绝对值越来越大。

图6 直接刚度Kavg随涡动频率变化Fig.6 Variation of direct stiffness with swirl frequency

图7给出了直接阻尼系数Cavg随涡动频率的变化趋势。可以看出，本文计算结果与试验基本吻合，约为600 N·s/m。

图7 直接阻尼Cavg随涡动频率变化Fig.7 Variation of direct damping with swirl frequency

图8给出了交叉刚度Kxy随涡动频率的变化趋势，计算与试验值随涡动频率增加，基本保持不变，约为50 kN/m。

图8 交叉刚度Kxy随涡动频率变化Fig.8 Variation of cross stiffness with swirl frequency

有效阻尼系数是评价密封性能的关键参数，图9给出了有效阻尼随涡动频率的变化趋势。可以看出，计算结果略大于试验值，当涡动频率大于100 Hz时，有效阻尼系数随涡动频率增加基本保持不变。

图9 有效阻尼Ceff随涡动频率变化Fig.9 Variation of effective damping with swirl frequency

3.2 密封静态稳定性分析

上述理论与试验结果都表明，该密封具有较好的动态稳定性，即具有正的有效阻尼系数。然而，本文与Ertas等的研究结果均表明直接刚度系数呈现负值，如图6所示。密封静态力是转子零转速或没有涡动情况下受到的静态激振力，该力对转子的影响主要表现为改变转子系统刚度，引起转子临界转速的变化，进而导致转子不稳定振动区域的改变；另一方面对试验气缸有负面影响，即会导致气缸与转子的碰摩，Picardo等[33-34]为了避免这种负面影响，专门在气缸上额外增加支撑结构(stiffener)。为此下文对该密封内静态力的形成原因开展进一步分析。

表3给出了静态时转子受力与泄漏情况。静态工况下，其在偏心状态下所受气流激振力均为正值，即直接刚度为负值，与图6识别结果相吻合。转子在偏心自转时会产生垂直于偏心方向的切向力，即产生正的交叉刚度系数。密封泄漏量在转子转动情况下会略微降低。

表3 密封静态特性

图10给出了静态下密封内小间隙与大间隙压差沿泄漏方向变化情况。除第一个密封腔，其余部位的小间隙与大间隙的压力差值小于0，即产生负刚度，加剧转子偏离中心程度。

图10 密封间隙压差沿泄漏路径分布Fig.10 Pressure difference distribution in the seal clearance along leakage path

图11给出了密封最大和最小间隙内气流速度沿泄漏方向分布情况。小间隙处速度均大于大间隙，惯性力占主导作用，小间隙压力小于大间隙，导致出现负刚度，即静态不稳定。

图11 密封内大小间隙处速度分布Fig.11 Velocity distribution in the seal clearance

4 结 论

(1)应用瞬态动网格方法，提出基于微元轨迹的密封动力特性系数理论识别方法，可满足任意转子椭圆涡动轨迹与任意涡动频率的密封动力特性系数求解。

(2)研究表明，本文理论计算与试验结果具有很高的吻合度，尤其对衡量密封系统稳定性的有效阻尼系数，具有较高的识别精确度。涡动频率对试验密封直接阻尼系数、有效阻尼系数和交叉刚度系数影响较小。

(3)试验密封具有静态不稳定性。计算表明密封最大间隙气流速度小于小间隙，惯性力占主导作用，大间隙处平均压力大于小间隙处，压力差加剧了转子偏心，造成密封系统静态不稳定。