挠性支承可倾瓦轴承动力特性研究

杨期江， 李伟光， 赵学智， 郭明军

(1. 广州航海学院 轮机工程学院，广州 510725; 2. 华南理工大学 机械与汽车工程学院，广州 510640)

传统机械支点可倾瓦轴承通过轴瓦摇动或滑动，实现轴承刚度与阻尼系数低交叉耦合。挠性支承可倾瓦轴承(Flexure Pivot Tilt Pad Bearing，FPB)通过挠性支承(如图1所示)，实现轴瓦支承点的弯扭耦合，达到相同的低交叉耦合和高稳定性，同时消除了传统机械支点的磨损和高接触应力。

图1 传统机械支点与挠性支点Fig.1 Traditional mechanical pivot and flexible pivot

轴承间隙是滑动轴承设计的重要参数，影响轴承转子稳定性，临界转速和动态响应。挠性支承可倾瓦轴承采用电火花技术加工轴瓦，瓦块与壳体之间通过挠性梁连接，使其成为一个整体。这种设计消除了瓦块与壳体的装配误差(如图2所示)，支点疲劳，瓦块卸载的颤振问题等传统可倾瓦轴承的缺陷。

图2 传统机械支承与挠性支承可倾瓦轴承公差带Fig.2 Traditional mechanical pivot and flexible pivottilting pad bearing tolerance

国外Armentrout等[1]推出了全新的一体式可倾瓦轴承的设计，即挠性支点可倾瓦轴承(FPB)，相对与传统可倾瓦滑动轴承(Tilting Pad Journal Bearing, TPJB)，它在转子动力学方面更具有优势。利用电火花加工(Electrical Discharge Machining, EDM)工艺整体线切割出 FPB挠性支点、轴瓦基体和瓦座。Kepple等[2]报道第一个FPB已安装于1994年，讨论了可倾瓦轴承的缺点，如支点磨损，设计，生产，安装和维护的复杂性，并将传统可倾瓦轴承与新的轴承设计“挠性支承可倾瓦轴承(FPB)”进行比较，提供了FPB应用在高速涡轮压缩机上所测得的轴承-转子稳定性试验数据，发现挠性支承与传统可倾瓦轴承在稳定性方面基本一致。Zeidan[3]讨论了存在于所有的旋转机械中油膜轴承交叉耦合特性及其对转子系统稳定性的影响。接着介绍了之前所广泛使用可倾瓦轴承的一些问题和缺点，并讨论用FPB来解决这些问题。Chen等[4]在高速一体齿轮压缩机将其与传统的轴承进行比较，并阐述了优化设计之后FPB的效果，即降低了摩擦损失并具有更宽的工作温度范围。De Choudhury等[5]研究了在高速离心压缩机上，具有一个新的径向和推力轴承挠性支承5瓦可倾瓦轴承和推力轴承系统中的效果。他们的测试显示出较低的油温上升，较少的摩擦功率损失。目前该领域发展的方向：①与挤压油膜阻尼器串联组成全新的阻尼减振滑动轴承[6-7]；②将该类型结构的轴承应用于无油透平机械-如微型燃气轮机，改用空气作为润滑剂，大大拓宽了其工作转速与温度[8]；③引入流体静压，发展成为动静压混合润滑轴承[9-10]。

国内Feng等[11]提出了一种挠性支撑可倾瓦与金属橡胶串联的空气轴承，建立了理论计算模型，并对其阻尼特性进行了研究。但目前国内在这一领域处于起步时期，本文从挠性支承可倾瓦轴承的瓦块与轴颈的动力学特性出发，研究挠性支承可倾瓦轴承频率缩减动力学模型，提出一种基于PDE工具箱快速求解轴承静态以及频率扰动压力Reynolds方程，考虑可倾瓦轴承的瓦块自由度，采用Newton-Raphson 迭代法可计算得到给定载荷和转速工况下的轴承的静平衡位置的数值计算方法。并与国外的最新试验数据进行对比分析，验证本文的动力学模型及仿真计算方法。为该类型轴承在透平涡轮上的应用做好理论工作准备。

1 流体动力润滑的基本方程

润滑油膜的惯性力，以及热效应影响油膜动力及减振特性。本文引入带惯性项的不可压缩非定常工况Reynolds 方程[12]

(1)

滑动轴承油膜二维偏微分能量方程为[13]

(2)

式中：x为轴承周向坐标；y为径向坐标(膜厚方向)；z为轴承间隙轴向坐标；h为油膜厚度；μ为滑油动力黏度， Pa·s；P为油膜压力，Pa；t为时间，s；Θ为周向位置角。

2 挠性支承可倾瓦轴承动力学模型

可倾瓦轴承几何及坐标系如图3所示，x方向为轴瓦圆周方向，z方向为轴线方向。 全局坐标系(X,Y,Z)以轴承中心为原点，eX,eY为轴心在X,Y方向上的位移。 瓦块位置角θ=x/R。

图3 挠性支承瓦块坐标关系及支点运动位移Fig.3 Flexure pivot tilting pad geometry and pivot displacement

图3中Θl为轴瓦起始边角度，Θt为轴瓦终止边角度，Θp为轴瓦支点位置角； 挠性支承可倾瓦轴承支承位置具有柔性，轴颈在旋转过程中，瓦块在支点处具有摆角(δk)， 径向(ξk)以及切向(ηk)位移。

2.1 瓦块运动方程

所有瓦块油膜力的矢量和必须与轴颈外部载荷(WX,WY)平衡：

(3)

第kth块瓦运动方程为

(4)

(5)

2.2 油膜力与力矩

油膜力与力矩可分为静态与动态部分

(6)

其中静态平衡位置(0阶)油膜力与力矩为

(7)

式中： {Zαβ}为油膜力动刚度系数， 油膜力动刚度系数是可以根据轴心平衡位置(eX0,eY0)得到

(8)

式中：hx=cosθ，hy=sinθ。

2.3 扰动分析

当轴颈中心在静态平衡(eX0,eY0)处做频率为ω， 幅值为(ΔeX,ΔeY)的小扰动。瓦块也做同频小扰动。因此轴颈、瓦块的扰动方程可表示为：

eX(t)=eX0+ΔeXeiωt，eY(t)=eY0+ΔeYeiωt，

式中：k表示瓦块数， 即为：k=1, 2, …,Npad轴瓦油膜厚度与油膜压力分布都由静态域(0阶)与扰动域(1阶)组成，即为

(9)

式中： Δhk=ΔeXcosθ+ΔeYsinθ+Δξkcos(θ-Θp)+(Δηk-RΔδk)sin(θ-Θ), ΔPk=PXΔeX+PYΔeY+PξΔξk+PηΔηk+PδΔδk。

2.4 缩减动力系数

可倾瓦轴承动力系数取决于轴颈与轴瓦平衡位置，以及特定的激振频率(ω)，通常认为同步激振频率(ω=Ω)，或者不同步激振频率(ω<Ω)。轴承的总自由度为2+3Npads即轴颈中心位移沿{X,Y}轴线，为2自由度，另一方面，每块瓦有1个转动自由度和2个径向及切向位移自由度。由此可知，可倾瓦轴承包含了轴颈及轴瓦的耦合运动。

假设轴瓦与轴颈做频率为(ω)的简谐运动，将式(6)代入式(3)、式(4)中得到频率缩减的动刚度系数。

(10)

3 支点刚度

对于挠性支承轴瓦，将挠性支承轴瓦视作一端固定，一端自由的悬臂梁，而轴瓦本身被视为集中惯量[14]，具体如图4所示。

图4 挠性支点简化悬臂梁Fig.4 Cantilever bearm model of a tilting pad with flexural web

(11)

(12)

式中：mp是瓦块质量;ro,p和rr,p是内、外瓦块半径；ρ为轴承材料的密度；l为瓦块轴向长度；Larc,p为垫的弧长。

忽略挠性支承的圆角，同理可得挠性支点的压缩与切向刚度为[14]

(13)

式中：E为材料弹性模量；w为挠性支承宽度；L为轴瓦宽度；h为挠性支承高度。

4 平衡位置迭代算法

可倾瓦轴承动态特性的传统差分法仿真思路多是在给定偏心率的条件下进行编程计算，需经过插值方法，曲线拟合出实际轴承工作转速与载荷下的油膜力以及平衡位置，再调用扰动压力Reynolds方程计算程序，进而积分得到实际轴承工作转速与载荷下的动力系数。但此方法过于繁琐。本文提出采用Newton-Raphson 迭代法可直接计算得到给定载荷和转速工况下的轴承的静平衡位置。联合Matlab以及Comsol PDE工具包，设计程序结构和循环迭代算法。最终可获得平衡状态下静态特性、缩减动力学系数。其静态平衡位置迭代计算的基本流程如图5所示，具体过程如下：

图5 挠性支承可倾瓦轴承的瓦块及轴颈中心静平衡位置迭代算法Fig.5 Flow chart of an algorithm to find the equilibriumposition of the the pads and journal in the flexible-pivottilting pad bearing

(14)

瓦块摆角平衡位置计算迭代公式为

(15)

(2) 采用的Newton-Raphson方法对作用于轴心位置的加载力进行平衡迭代计算，即

F0-W0

(16)

轴心位置的迭代计算公式为

(17)

(18)

径向位移迭代公式为

(19)

5 算例分析

挠性支承可倾瓦轴承动力特性仿真计算流程具体如图6所示。首先根据挠性支承可倾瓦的动力学模型，采用PDE数值计算非定常工况Reynolds方程及能量方程，根据图5所示的平衡位置计算及本文第5部分的缩减动刚度系数计算方法，获得给定工况下的轴心及瓦块平衡位置与频率缩减动刚度系数，再根据最小二乘法曲线拟合得到刚度Kij、阻尼Cij和附加质量Mij。

图6 挠性支承可倾瓦轴承动力特性仿真计算流程Fig.6 The flexuret ilting pad bearing support dynamiccharacteristics simulation process

5.1 挠性支承可倾瓦轴承参数

本章采用文献[16]中的四瓦挠性支承可倾瓦轴承进行仿真计算对比，具体轴承工况条件如表1所示。

表1 轴承工况条件

其中四瓦挠性支承可倾瓦轴承的参数具体如表2所示。

表2 轴承参数表

挠性支承可倾瓦轴承润滑油牌号为ISO VG32-透平油，其温黏方程为

μ(T)=0.045 3e-0.030 07(T-21.11)

(20)

5.2 挠性支承可倾瓦轴承动力参数识别

轴承动力参数实验采用倒置式试验台， 轴承安装在轴承座内部。图7描绘了轴承油膜八个线性化刚度和阻尼矩阵模型。

图7 线性化刚度与阻尼系数Fig.7 Linearized bearing coefficients

轴承座(质量为Ms)运动方程可以写为

(21)

(22)

式(22)是轴承油膜线性力-位移模型， 动力系数包括刚度Kij、 阻尼Cij和附加质量Mij。 Δx、 Δy表示轴颈和轴承座之间的相对运动位移。

将式(21)代入式(22)中，并进行傅里叶变换可得

(23)

式(23)左边为已知测量量， 式中Fx,Fy为输入激励力频域值；Ax、Ay为轴承座加速度频域值。因此4个轴承油膜动刚度未知量Hxx、Hxy、Hyx、Hyy可通过式(23)求出。

而轴承油膜动刚度与动力参数之间的关系为：

Hij=(Kij-Ω2Mij)+j(ΩCij)

Re(Hij)=Kij-Ω2Mij, Im(Hij)=ΩCij

式中： 下标i,j表示x,y方向，Ω是激振频率。因此动力参数(刚度Kij、 阻尼Cij和附加质量Mij)可采用最小二乘法对轴承油膜动刚度的实部与虚部进行曲线拟合得出[17]。

5.3 油膜动刚度系数

根据图6的挠性支承可倾瓦轴承动力特性仿真计算流程，采用考虑流体惯性的Reynolds方程，通过建立挠性支承可倾瓦轴承的频率-缩减动力学模型，得到在不同转速，载荷的工况条件下动刚度系数实部与虚部的动刚度系数计算结果，具体如图8所示。

转速为10 000 r/min，轻载(348.8 kPa)与重载(1 038.2 kPa)情况下，四瓦挠性支承可倾瓦轴承动刚度系数理论计算与试验对比分析如图7所示。随着激振频率的增加，动刚度实部直接项(Re(Hxx)与Re(Hyy))呈二次方下降趋势，轻载时曲线下降斜率比重载时要小，说明重载时油膜刚性大。随着激振频率的增加，动刚度实部交叉项(Re(Hxy)与Re(Hyx))逐渐降低，但降幅较小。随着激振频率的增加，动刚度虚部直接项(Im(Hxx)与Im(Hyy))呈一次方上升趋势。轻载时曲线上升斜率比重载时要大，说明重载时油膜阻尼大。随着激振频率的增加，动刚度虚部交叉项(Im(Hxy)与Im(Hyx))先上升后逐渐降低，但升降幅度十分微小。由此可知动刚度实部直接项(Re(Hxx)与Re(Hyy))理论计算与试验分析比较接近，但动刚度实部交叉项(Re(Hxy)与Re(Hyx))理论计算与试验分析结果有一定的差异。且动刚度虚部直接项(Im(Hxy)与Im(Hyx))理论计算与试验分析理论计算与试验分析结果的较吻合，但动刚度虚部交叉项(Im(Hxy)与Im(Hyx))理论计算与试验分析结果的差异比较明显。总体上10 000 r/min重载时，四瓦挠性支承可倾瓦轴承动刚度实部直接项(Re(Hxx)， Re(Hyy))与虚部直接项(Im(Hxx)， Im(Hyy))理论计算与试验分析结果更加吻合，主要原因是重载时油膜厚度比轻载时油膜厚度小，油膜惯性力影响减小，因此理论计算更接近于试验结果。

图8 10 000 r/min轻载与重载情况下理论动刚度与试验动刚度随频率变化曲线Fig.8 Theoretical and Experimental dynamic stiffnesses versus the excitation frequency at 10 000 r/min and light, high unit load

5.4 油膜刚度、阻尼与附加质量系数

应用本文算法(如图6所示)，采用最小二乘法拟合得到在不同转速，载荷的工况条件下油膜刚度、阻尼与附加质量系数。对不同转速下刚度、阻尼、附加质量系数随载荷变化的理论计算结果与文献[16]中试验数据曲线进行对比分析具体如图9、图10、图11所示。

(1) 油膜刚度系数

直接刚度项Kxx与Kyy理论计算结果是相等的，随着载荷(>3.1 kN)增大， 直接刚度项Kxx与Kyy理论计算与试验结果都呈线性增大趋势；当高转速(10 000 r/min与12 000 r/min)与低载荷(<6.16 kN)下理论计算数值比试验结果要低，高转速(10 000 r/min与12 000 r/min)与重载荷(≥6.16 kN)下理论计算数值与试验结果更加接近，说明高速重载下的直接刚度项Kxx、Kyy本文理论计算数值与试验结果较吻合； 交叉刚度项Kxy与Kyx在上述工况载荷下理论计算数值与试验分析结果也都较吻合。

图9 不同转速下刚度系数随载荷变化Fig.9 Stiffness coefficients vs. static loads fordifferent rotor speeds

(2) 油膜阻尼系数

直接阻尼项Cxx与Cyy理论计算结果是相等的。随着载荷(>3.1 kN)增大， 直接阻尼项Cxx与Cyy理论计算呈线性缓慢增大趋势。高转速(10 000 r/min与12 000 r/min)与低载荷下(<6.16 kN)理论计算数值比试验结果要低，高转速(10 000 r/min与12 000 r/min)与重载荷下(≥6.16 kN)理论计算数值与试验结果更接近，如说明高速时， 重载下的直接阻尼项Cxx与Cyy本文理论计算数值与试验结果较为吻合。 交叉刚度项Cxy与Cyx理论计算结果变化不明显，且在上述工况条件下理论计算数值与试验分析结果也都较吻合。

(3) 油膜附加质量系数

直接附加质量项Mxx与Myy理论计算结果是相等的。 直接附加质量项Mxx与Myy理论值随着载荷的增大呈线性增大趋势，而试验结果在高转速(10 000 r/min与12 000 r/min)与低载荷(<6.16 kN)下理论计算数值比试验结果要低，高转速(10 000 r/min与12 000 r/min)与重载荷下(≥6.16 kN)理论计算数值与试验结果更接近，如工况(转速12 000 r/min，轻载3.1 kN)下， 直接质量Mxx理论值为21.24 kg， 试验值为25.5 kg， 误差为16.7%； 工况(转速12 000 r/min，重载荷9.24 kN)下， 直接阻尼Mxx理论值为29.01 kN·s/m，试验值为27.94 kN·s/m，误差为3.82% ，说明高速时，重载下的直接质量项Mxx与Myy本文理论计算数值与试验结果较为吻合。 交叉质量项Mxy与Myx理论计算结果变化不明显，且在上述工况条件下理论计算与试验分析结果也都较吻合。

图10 不同转速下阻尼系数随载荷变化Fig.10 Damping coefficients vs. static loads in fordifferent rotor speeds

图11 不同转速下质量系数随载荷变化Fig.11 Added-mass coefficients vs. static loads in fordifferent rotor speeds

6 结 论

(1) 引入带惯性项的Reynolds方程，考虑动压润滑油膜温黏效应；首先在轴瓦与轴颈作同频振动的假设的基础上，根据瓦块支点受力的平衡方程，静平衡位置瓦块运动的微分方程等，建立了柔性阻尼支承可倾瓦轴承频率-缩减动力学模型。

(2) 提出一种挠性支承可倾瓦轴承的静平衡位置迭代计算方法，即采用PDE工具箱有限元数值方法可快速求解轴承非定常工况Reynolds 方程及二维能量方程，将轴承工作载荷和转速作为初始条件，并考虑可倾瓦轴承的瓦块摆角自由度，通过Newton-Raphson 迭代法可直接计算得到给定载荷和转速工况下的轴承的静平衡位置。最终可获得平衡状态动力学系数。

(3) 采用本文提出的动力学模型及平衡位置迭代计算方法对挠性支承可倾瓦轴承进行了动态特性参数分析研究，并与相关文献的理论及试验仿真结果进行了对比，验证了本文提出的动力学模型，平衡位置迭代算法。