重型燃气轮机可倾瓦轴承瓦块摆振特性研究

2021-11-17冯永志孟凡刚金英泽袁小阳陈予恕

振动与冲击 2021年21期

冯永志，闫岗，孟凡刚，金英泽，袁小阳，陈予恕

(1.哈尔滨工业大学航天学院，哈尔滨 150001；2.哈尔滨电气股份有限公司，哈尔滨 150028；3.西安交通大学现代设计及转子轴承系统教育部重点实验室，西安 710049)

重型燃气轮机广泛应用于能源发电和工业驱动等领域，是国家机械工业的重大核心装备，而轴承作为各类旋转机械中的关键功能保障部件，其稳定性直接影响燃气轮机的工作性能[1]。可倾瓦轴承具有相互独立的瓦块，这些瓦块能够随着载荷波动、转速变化而绕其背部支点做自适应摆动，减小了切向力，具有比其他轴承更高的稳定性，可以为燃气轮机的安全运行提供保障，因而在现代大型机组上得到了广泛应用[2]。

虽然可倾瓦轴承具有很高的稳定性，但并非完全稳定，虞烈[3]从理论上证明了在不考虑瓦块惯性、支点摩擦力及支点变形的情况下，可倾瓦轴承依然存在不稳定因素。在工程上，可倾瓦轴承同样出现了失稳现象。某西门子9F燃气轮机6号轴承在SSS离合器啮合或脱扣期间发生低频振动突变[4]；某燃气轮发电机与汽轮机啮合过程中经常出现幅值超过500 μm的突变振动，导致机组无法安全运行[5]；某1 000 MW汽轮机组在调试运行期间可倾瓦轴承上发生了非稳态振动[6]；某660 MW核电汽轮机3号轴承在机组调门试验、启停机以及改变机组功率时出现低频振动分量，使3号轴承振动幅值达到190 μm，接近机组停机限值[7]。这些理论结果与故障现象均证明了可倾瓦轴承并非天然稳定。

这些低频失稳故障的诱因有很多，其中油膜失稳有相对成熟的理论能够解释，可以通过减小轴承顶隙、增大轴承侧隙、抬高轴承标高等方法来消除不稳定振动[8-9]。由系统内其他部件引起的耦合振动可以通过振动成分对比找到诱因，进而解决振动问题，如由轴向振动引发的径向低频振动可通过减小推力轴承间隙值来消除。还有一些失稳现象与上述诱因导致的低频故障特征很相似，但难以用传统理论解释，虽然这种失稳现象也呈现低频特征但它却容易导致上瓦疲劳损伤[10]，这种失稳现象的工程表现为上瓦周期性的拍打轴颈。因此，业界将这种失稳现象定义为“拍瓦”故障，Yang等[11-12]指出这种故障现象与瓦块摆振有关，张楚等[13]对这种由瓦块摆振引发的低频振荡问题进行了研究并比较了瓦块颤振与油膜失稳之间的差别。但目前关于可倾瓦轴承的故障研究多集中于油膜失稳，对瓦块子结构的研究并不透彻，为了消除由瓦块引起的故障现象，有必要对瓦块的摆振特性进行研究。王永亮等[14-15]基于稳态雷诺方程分析了轴承参数对瓦块摆动特性的影响规律。但从瓦块动态响应着手进行的研究鲜有报导，鉴于瓦块动态时变非线性效应极强的工作状态，有必要针对此特性进行润滑及动力学建模并展开相关分析。

因此，本文基于考虑紊流效应的动态雷诺方程，建立了计入瓦块惯性的可倾瓦轴承瓦块摆振动力学分析模型，利用变量分离法推导了可倾瓦块油膜压力的近似解析表达式，采用欧拉法计算了瓦块的动态响应情况，并通过试验验证了分析模型和解析计算方法的正确性，研究了预负荷系数对瓦块摆振特性的影响规律，细化分析了阶跃激励下的瓦块响应，可为研究由瓦块摆振引起的可倾瓦轴承失稳故障提供参考。

1 可倾瓦块摆振分析模型

1.1 动力学方程

本文的主要研究对象为四瓦可倾瓦轴承，分析模型如图1所示，分析时将瓦块视为单自由度系统考虑瓦块绕支点摆动，不考虑瓦块变形，瓦块受重力和油膜力的复合作用。

图1 可倾瓦轴承动力学分析模型Fig.1 Dynamic analysis model of tilting pad bearing

轴颈和轴瓦的无量纲运动方程可描述为

(1)

惯量J由下式给出[16]

(2)

式中：ρP为瓦块密度；L为瓦块长度；α为瓦包角；R为轴颈半径；hP为瓦块厚度；ξ为支点系数。

1.2 油膜力计算

在等温、不可压缩和无惯性流动的条件下，第i块轴瓦的无量纲紊流雷诺方程为

(3)

式中：Pi为无量纲的液膜压力；Hi为无量纲的液膜厚度；η为长径比；τ为无量纲的时间；Φ为无量纲周向坐标；λ为无量纲轴向坐标；KΦ和Kλ为紊流系数。

无量纲膜厚Hi由下式给出

(4)

式中：A为预负荷系数；βi为第i块瓦的支点位置角；ε为偏心率；θ为偏位角；δi为第i块瓦的无量纲的瓦块摆角；ψ=cp/R为间隙比，其中cp为瓦块径向间隙，R为轴颈半径。

Ng-Pan模型[17]紊流系数KΦ和Kλ由下式给出

(5)

式中,Re为雷诺数。

液膜流态有三种情况：层流状态、过渡状态和紊流状态，它们由两个临界雷诺数800和1 500来界定[18]。

Re=

(6)

式中：c1-4是通过确保Re及其在过渡带开始和结束处的一阶导数的连续性来确定的，c1=-23/3 430 000，c2=409/17 150，c3=-8 672/343，c4=2 880 000/343；Rem是平均雷诺数；Rel是局部雷诺数。

平均雷诺数Rem和局部雷诺数Rel可以表示为

(7)

式中：ρ是流体密度；ω是旋转角速度；μ是流体黏度；Hm是无量纲平均膜厚；cp为瓦块径向间隙。

Hm定义为

(8)

式中：Φin为进液边位置角；Φout为出液边位置角。

在动态Gümbel边界条件下，使用变量分离法求解方程式(3)，无量纲膜厚可以表示为

Hi=1+εicosφ(φ=φ-θi)

(9)

式中

(10)

Pi(φ,λ)=PL(φ)P*(λ)

(11)

式(11)中PL为无限长轴承压力，由下式给出

(12)

式中

(13)

P*由下式给出

(14)

(15)

瓦块i在X和Y方向上的无量纲承载力和无量纲力矩为

(16)

1.3 可倾瓦轴承响应计算方法

采用欧拉法求解轴颈和瓦块的运动方程，得到轴颈和瓦块的不平衡响应，可倾瓦轴承在下一时刻的状态变量可通过当前时刻的状态变量和时间步长计算出来，如式(17)

科学家们早已发现可以通过从大气层中提取二氧化碳来减缓二氧化碳带来的温室效应。但已有的解决方案都严重受限于实施难度，特别是经济困难。菱镁矿可以将二氧化碳固定在体内，从而减少大气中的二氧化碳含量。一吨自然形成的菱镁矿可以固定将近半吨的二氧化碳，但这个固化过程非常缓慢。现在，研究者们有史以来第一次描述了在低温下合成菱镁矿的技术，同时阐述了一项可以大幅加快菱镁矿结晶过程的技术。目前，我们认识到这还只是一个实验中的解决方案。在我们确定菱镁矿可以用于碳封存之前，还需要验证该过程能否被放大到相当的规模。这取决于几个因素，包括碳的价格和封存技术的改进，但我们现在至少知道是可能有可行性的。

(17)

式中：u为可倾瓦轴承的状态向量；Δt为计算时间步长。

以本文计算分析所用的四瓦可倾瓦轴承为例，上式可展开为

(18)

不平衡响应的求解流程为：(1)给出可倾瓦轴承结构和运行参数，选定轴颈初始位移、速度和各瓦块的初始摆角、摆动速度以及时间步长和计算终止时间；(2)根据解析法计算轴承的液膜力和各瓦块所受力矩；(3)根据欧拉法计算下一时刻的轴颈位移和瓦块摆角；(4)判断是否达到计算终止时间，若已经达到则输出计算结果，若未达到则重复步骤(2)和步骤(3)，直至满足计算终止时间。本文计算时取轴颈初始位移、速度和瓦块的初始摆角、摆动速度均为0，取时间步长为1×10-6s。

2 理论模型的试验验证

本文提出的解析方法计算效率可比有限差分法等其他方法高几个数量级，能够为可倾瓦轴承转子系统大规模非线性时变润滑动力学运算提供高效仿真工具，为验证该解析方法的正确性和精度，同时验证动力学分析模型的正确性，本章利用卧式转子-可倾瓦轴承试验台对理论模型及计算方法进行了试验验证。

2.1 试验装置

为验证分析模型及解析计算方法的正确性，在燃气轮机模化转子-轴承试验台上进行了可倾瓦轴承的瓦块摆振特性试验。试验转子为燃气轮机模化转子，由55 kW三相异步电机经齿轮箱驱动，试验轴承为四瓦可倾瓦轴承，实物图如图2所示，结构及运行参数如表1所示，图3为试验台实物图。

图2 试验用可倾瓦轴承Fig.2 Tilting pad bearing for text

图3 试验台实物图Fig.3 Physical map of test bench

表1 可倾瓦轴承结构运行参数

2.2 测试方案

本试验通过正对瓦背安装的电涡流传感器测量瓦块摆动时瓦背测点处的位移，受试验台结构所限，无法测量下瓦，因此仅对上瓦进行测量。电涡流传感器通过支撑轴承上所钻的螺纹孔固定，为防止干扰保证传感器灵敏度，传感器探头处的孔径要大一些，考虑到支撑轴承与轴承箱之间的装配精度问题，轴承箱上不加工螺纹孔，钻一个孔径略大于螺纹孔径的通孔，这样可保证即使存在装配精度问题传感器仍可顺利安装，具体如图4所示。

图4 传感器安装图Fig.4 Sensor installation diagram

瓦块的摆角时域响应与试验测得的位移曲线存在对应关系，可通过一定的数学变换将其统一到同一维度中进行对比，以3#瓦块为例进行说明，计算时采用忽略瓦块厚度的简化支点假设，几何关系如图5所示。

图5 几何关系示意图Fig.5 Geometric relation diagram

根据余弦定理可知，瓦块摆角为δ3时测点对应的位移s3为

(19)

式中,L为支点到测点之间的距离。

2.3 试验结果

试验时首先根据轴振识别不平衡量，然后计算在此激振源下的瓦块摆角响应，根据式(19)所示的转换关系将计算得到的3#瓦块的摆角时域响应转换为测点的位移响应，并与试验结果进行对比，如图6所示。可以看出，本文理论模型计算出的瓦块运动状态与试验得到的瓦块运动状态在时域波形上的变化趋势及摆动幅值是相同的，仅在瓦块由正变负和由负变正的极值点附近存在差异，这是因为实际工况中瓦块摆角换向时受到多种因素影响其非线性效应更强，因而不会像理论计算结果一样平滑，但在正常摆动过程中，理论模型可以很好的反映瓦块的运动状态。

(a) 试验结果的时域波形图

(b) 理论计算结果的时域波形图图6 试验结果与理论结果的时域波形对比图

图7对比了试验结果与理论结果的频谱图，可以看出，试验结果与理论计算结果均显示瓦块的摆动频率为转子的转动频率(30 Hz)，但理论计算结果的频谱图中存在倍频，这是因为理论计算时所取的时间步长为1×10-6s，且计算结果的时域波形变化较为平滑，而试验中的采样频率较低，数据点较少，时域波形存在突变及尖峰，因此造成了此项差异。

综上，虽然理论计算结果与试验结果存在细微差异，但这些差异不是分析模型的固有缺陷而是由数值计算带来的，计算模型的精度对分析规律来说是足够的，因此理论分析模型的正确性得以验证。

3 案例分析与讨论

3.1 分析参数

间隙是影响可倾瓦轴承瓦块摆振特性的关键参数，间隙有瓦块径向间隙和轴承径向间隙两种，预负荷系数可以表征两种间隙之间的关系问题，是轴承设计和计算的重要参数之一。

(a) 试验结果的频谱图

(b) 理论计算结果的频谱图图7 试验结果与理论结果的频谱对比图

假设轴承瓦块支点为刚性且位于同一个支点圆上，预负荷系数与瓦块支点的径向位置和瓦面的曲率半径有关，如图8所示。可倾瓦轴承预负荷的大小可用预负荷系数A表示

(20)

式中：cb为轴承径向间隙，cb=R1-R，R1为支点圆半径，R为轴颈半径；cp为瓦块径向间隙，cp=R2-R，R2为瓦面的曲率半径。

图8 可倾瓦轴承预负荷系数的定义Fig.8 Definition of preload factor of tilting pad bearings

本文计算采用的轴承参数如表2所示。

3.2 时域响应分析

图9比较了预负荷系数分别为0.4、0.5、0.6、0.7、0.8时，四块瓦的摆角时域响应。由图9可知，1#、3#、4#瓦块的摆角随预负荷系数的增加呈现先增大后减小的趋势，2#瓦块的摆角随预负荷系数的增大而减小。这是因为2#瓦为主承载瓦，1#、3#、4#瓦为非主承载瓦，所以呈现出了特殊的摆动特性。预负荷系数从0.4增加到0.6的过程中，1#瓦块和3#瓦块的摆角随预负荷系数的增加而增大，预负荷系数从0.6增加到0.8的过程中，1#瓦块和3#瓦块的摆角随预负荷系数的增加而减小，由数据分析可知，在0.5～0.7之间存在一个会使1#瓦块和3#瓦块摆角最大的预负荷系数。预负荷系数从0.4增加到0.7的过程中，4#瓦块的摆角随预负荷系数的增加而增大，预负荷系数从0.7增加到0.8的过程中，4#瓦块的摆角随预负荷系数的增加而减小，这说明在0.6～0.8之间存在一个会使4#瓦块摆角最大的预负荷系数。

表2 可倾瓦轴承计算参数

3.3 频谱特性分析

图10给出了预负荷系数分别为0.4、0.5、0.6、0.7、0.8时，四块瓦的摆角频谱图。由图10可知，各瓦块在不同预负荷系数下的摆动频率基本就是转动频率50 Hz，预负荷系数较大时，各瓦块均无明显倍频，预负荷系数较小时，1#瓦和2#瓦有微小倍频，3#瓦和4#瓦有明显倍频。各瓦块的振动幅值均随预负荷系数的增大而减小，减小的程度十分明显，均在90%以上。可见预负荷系数对瓦块摆振特性有着重要影响，且非承载瓦在预负荷系数较小时有特殊的摆振特性。