箔片轴承转子系统的非线性动力学特性研究*

雷 群 宋宜强 杜建军

(1.广州市昊志机电股份有限公司 广东广州 511356;2.哈尔滨工业大学(深圳) 广东深圳 518055)

近年来，随着精密加工技术和超高速旋转机械的应用，人们对其中关键的部件——轴承的性能要求越来越高。从20世纪起，气体润滑技术发展迅速，气体轴承开始被广泛应用。气体轴承是用气体代替润滑油进行润滑、使用气膜提供支撑力的机械构件，特别适合于轻型、高速的涡轮机械，例如涡轮压缩机和空气循环机等[1-3]。气体箔片轴承是气体动压轴承的一种，由于其出色的性能，近年来越来越受到人们的重视。相比于传统的刚性气体动压轴承，气体箔片轴承具有寿命长、承载能力大、可靠性高和稳定性好的优点，并且可以在高温高速情况下工作。

1968年，LUND[4]在研究气体动压轴承时提出了摄动法，即用4个刚度系数和4个阻尼系数来计算气膜力。这种方法将油膜或气膜线性化表示，优点是建模简单、计算速度快，因此获得了广泛的应用[5-9]。但该方法仅限于小扰动及低速转动情况。轨迹法能够弥补摄动法这种不足，该方法通过联立求解雷诺方程、转子动力学方程和箔片变形方程，能够获得系统在各时刻的运动情况，因此有助于更好地研究系统的动力学特性，从而得到了广泛应用[10-13]。

在使用轨迹法时，一种常用的做法是将雷诺方程、转子动力学方程和箔片变形方程在每一个时间步内分开求解，由于各状态量不能同步更新，需要保持积分步长很小才能保证算法收敛，计算量很大。为了克服弱耦合解法的弱点，本文作者建立了箔片轴承转子系统的强耦合数学模型，实现雷诺方程、转子动力学方程和箔片变形方程的同步求解，从而使计算效率得到较大的提高。

1 非线性动力学耦合模型

图1所示是第一代气体箔片轴承的结构示意图。当转子发生偏移并高速旋转时，转轴与平箔片之间产生楔形间隙，从而为轴承提供承载力。气膜压力作用在平箔片上会使平箔片和波箔片产生弹性变形，文中将箔片结构视为普通弹性元件进行分析。在受到气膜力后，箔片结构的变形能够动态地调整气膜间隙。在求解系统的动态响应时，主要需要考虑3个方面：瞬态气膜压力的产生、转子受气膜力引起的涡动以及箔片结构受力变形。在建立相应的方程时，需要将一些相关的状态量表示出来。

图1 气体箔片轴承结构示意

1.1 气膜压力求解

在研究气体轴承转子系统的动力学特性时，通过求解可压缩时变雷诺方程得到气膜压力。在等温状态下，气体的时变雷诺方程为

(1)

式中：h为气膜厚度；p为气膜压力；z为轴向坐标；μ为气体黏度；t为时间；ω为转子转速；R为转子半径。

(2)

式中：θ为周向角度；pa为环境大气压；c为名义间隙；Λ为轴承数。

将式(1)变换为式(2)是为了方便后面建立状态方程。设箔片变形为q(θ)， 取向外为正，考虑箔片的变形，则气膜厚度为

h=c+xsinθ-ycosθ+q

(3)

式中：x为转子在x方向的位移；y为转子在y方向的位移；坐标轴方向在图1中显示。

(4)

使用有限差分法求解式(2)的雷诺方程，网格划分如图2所示。

图2 有限差分法的网格划分

计算域采用矩形网格离散为Nx×Nz个点，轴承的两端作为边界点，此处的气压等于外界大气压，不作为状态变量。对于气膜承载力，在计算中使用量纲一化形式，在考虑方向后积分得到计算公式：

(5)

(6)

式中:Fx和Fy分别为在x方向和y方向的气膜承载力;L为箔片轴承宽度。

1.2 转子动力学方程

转子动力学方程由牛顿第二定律获得。系统运行时，转子除了受到气膜力外，还需要考虑由于转子质心与形心不重合导致的自转惯性力。因此转子动力学方程表示为

(7)

(8)

式中:m为转子质量;Fub为不平衡力;S为重力;下标x和y分别表示在x和y方向。

对式(7)和(8)进行量纲一化，得到：

(9)

(10)

1.3 箔片变形方程

对于第一代波箔轴承，箔片沿轴向的变形可忽略不计[14]，箔片按照简单弹簧振动模型进行分析。根据文献[15]，箔片变形方程表示为

(11)

(12)

箔片的刚度计算通常使用HESHMAT等[16]在1983年提出的公式：

(13)

式中:E为弹性模量；tb为箔片厚度；ν为泊松比；l0为箔片波纹长度的1/2。

式(13)计算得到的箔片刚度是线刚度，进行量纲一化，式(11)对应的量纲一箔片面刚度为

(14)

式中：s为波箔节距。

1.4 系统的状态方程与耦合求解

为了建立整个气体箔片轴承转子系统的状态方程，实现强耦合求解，需要将时变雷诺方程、转子动力学方程和箔片变形方程以状态方程的形式表示。首先对时变雷诺方程进行处理，式(2)变换为

(15)

式中：∘为哈达玛积，表示向量的各项逐项相乘。

式(9)和(10)的转子动力学方程以向量形式表示为

(16)

式(11)的箔片变形方程同样以向量形式表示为

(17)

将式(15)、(16)和(17)联立起来，便得到箔片轴承转子系统的总体方程，以状态方程的形式表示为

(18)

其中，状态变量y为

(19)

求解式(18)的系统状态方程，便能同时得到瞬态气膜压力、转子位移和箔片变形量，实现强耦合求解，计算效率相比弱耦合解法有了较大提升。

2 仿真结果与分析

如无特别指出，文中仿真所选用的箔片轴承和转子的基本参数如表1所示。

表1 箔片轴承和转子的基本参数

2.1 模型验证

设转子初始位置的中心线与轴承的轴心重合，初始平动速度为0，转速为10 000 r/min，计算从量纲一时间为0～10π的转子运动轨迹与箔片变形，并与文献[17]中的数据比较来验证模型的正确性，如图3所示。

图3 转速为10 000 r/min时仿真结果与文献结果比较

可以看到，文中仿真得到的转子轨迹、箔片变形与文献[17]中的结果基本一致，说明文中建立的箔片轴承转子系统动力学模型是可行的。

2.2 转子转速对系统动力学特性的影响

转子转速对系统的稳定性和其他动力学特性产生重要影响。以转速为研究对象，使用不同转速，计算系统的运动状态。分别取转速为9 000、10 500和12 000 r/min，转子不平衡偏心量为1 μm，绘制系统进入稳定状态后，转子的运动轨迹和x方向的频谱图，如图4—6所示。

图4(a)—6(a)所示是转子的运动轨迹图，其横坐标为转子轴心的水平偏心率，纵坐标为轴心的竖直偏心率。可以看出，当转子转速为9 000 r/min时，转子的运动轨迹是一个规则的椭圆形；而在转速提高为10 500 r/min时，转子轨迹略微紊乱，但仍有迹可循；但当转速提升到12 000 r/min时，转子轨迹图变得复杂混乱，难以对其走向进行估计。

图4(b)—6(b)所示是对转子轴心在x方向的运动做快速傅里叶变换得到的频谱图。可以看出，当转速为9 000 r/min时，转子涡动与转速同频涡动；当提高转速到10 500 r/min时，除了1倍频外，还出现了0.3倍频和0.7倍频；当转速为12 000 r/min时，同样出现了多种频率，频率分量包含0.35倍频、0.65倍频和1倍频，且0.65倍频分量的振幅接近于1倍频分量的振幅，此时系统接近失稳。据频谱图可以知道，在转速较低时，转子做稳定的周期运动，但是当转速提高后，频谱图会出现其他频率分量，转子可能做拟周期运动。

图4 转速为9 000 r/min时转子的运动状态

图5 转速为10 500 r/min时转子的运动状态

图6 转速为12 000 r/min时转子的运动状态

图7所示是以转速为分岔参数绘制的转子在x方向和y方向的分岔图。分岔图是描绘状态变量随分岔参数变化的关系。可以看出，当转速达到9 500 r/min时，转子的运动出现分叉，说明此时转子进入拟周期运动状态。

图7 以转速为分岔参数的分岔图

2.3 转子质量对系统动力学特性的影响

改变转子质量进行数值仿真计算。取转子转速为10 000 r/min，转子质量分别为2、4和6 kg，绘制其运动轨迹图和频谱图进行分析，如图8—10所示。

图8 转子质量为2 kg时转子的运动状态

图9 转子质量为4 kg时转子的运动状态

图10 转子质量为6 kg时转子的运动状态

图8(a)—10(a)显示了转子的轨迹。可知，当转子质量为2 kg时，转子轨迹为规则的椭圆形；而在转子质量为4 kg时，轨迹变得比较复杂，且转子的振幅增大；当转子质量增大到6 kg时，转子轨迹进一步紊乱。这个规律与增加转子转速情况相似。

图8(b)—10(b)所示是对转子轴心在x方向的运动做快速傅里叶变换得到的频谱图。可以看出，当转子质量为2 kg时，转子涡动与转速同频涡动；在转子质量为4 kg时，除了基频，频谱图上还出现0.35倍频和0.65倍频；当转子质量为6 kg时，频谱图上出现0.45倍频、0.55倍频和1倍频。据此规律分析，当转子质量进一步增加时，将会出现0.5倍频的涡动，即为半速涡动。

图11所示是以转子质量为分岔参数绘制的分岔图。可知，当转子质量小于等于3 kg时，转子做周期运动；但是当转子质量大于3 kg时，开始出现其他频率的振动；随着转子质量的进一步增加，系统最终将失稳。

图11 以转子质量为分岔参数的分岔图

2.4 名义间隙对系统动力学特性的影响

以名义间隙为研究对象，分别取名义间隙为25、35和50 μm进行仿真，计算结果如图12—14所示。

图12(a)—14(a)显示了转子的轨迹。可知，名义间隙为25 μm时，系统不稳定，转子运动情况较混乱；而在名义间隙为35 μm时，转子的运动状态得到改善，运动轨迹偏移范围较小；当名义间隙为50 μm时，运动振幅进一步减小。

图12(b)—14(b)所示是对转子轴心在x方向的运动做快速傅里叶变换得到的频谱图。可以看出，当名义间隙为25 μm时，频谱图上除同频分量外，还出现0.2倍频和0.8倍频；在名义间隙为35 μm时，频谱图上主要频率为基频分量，其他另有0.7倍频和0.75倍频，但仅占很小部分；当名义间隙为50 μm时，主要频率成分仍为基频，还出现0.2倍频和0.8倍频。可以知道，当名义间隙增大到某一值时，再增加名义间隙对系统稳定性影响不大。

图12 名义间隙为25 μm时转子的运动状态

图13 名义间隙为35 μm时转子的运动状态

图14 名义间隙为50 μm时转子的运动状态

3 结论

(1)针对气体箔片轴承转子系统的动态响应进行研究，建立同步求解雷诺方程、转子动力学方程和箔片变形方程的强耦合数学模型，相比传统的弱耦合求解方法，能够进一步提升计算速度和提高计算精度，为箔片轴承的设计和性能预测提供了研究工具。

(2)探究不同设计参数对轴承转子系统性能的影响。通过改变转速、转子质量和名义间隙等参数进行数值仿真，并绘制转子运动轨迹图及频谱图，分析系统的稳定性。结果显示，当转子的质量增加或转速提高时，系统的稳定性呈下降趋势；反之，系统稳定性增强。另外，当名义间隙过小时，系统容易出现失稳状况，当名义间隙增大到某一值时，再增加名义间隙对系统稳定性影响不大。