固定瓦-可倾瓦微气体轴承-转子系统的非线性运动分析

张永芳, 吕烨迪, 肖良君, 赵晶群, 刘 成

(1. 西安理工大学 印刷包装与数字媒体学院，西安 710048；2. 西安交通大学 机械结构强度与振动国家重点实验室，西安 710049)

固定瓦-可倾瓦微气体轴承-转子系统的非线性运动分析

张永芳1,2, 吕烨迪1,2, 肖良君1,2, 赵晶群1,2, 刘 成1,2

(1. 西安理工大学 印刷包装与数字媒体学院，西安 710048；2. 西安交通大学 机械结构强度与振动国家重点实验室，西安 710049)

针对固定瓦-可倾瓦微气体轴承，考虑气体稀薄效应，给出了具有Burgdorfer一阶滑移速度边界的Reynolds方程，运用微分变换法结合有限差分法求解方程，得到了各瓦块在单块瓦坐标系中的非线性气膜力，然后通过组装技术获得了固定瓦-可倾瓦微气体轴承的气膜力。针对固定瓦-可倾瓦微气体轴承支撑的刚性Jeffcott转子系统，运用转子中心轨迹图、Poincaré映射图和时间历程图分析了转子的不平衡响应，比较了努森数、转速等对转子系统非线性特性的影响。数值结果表明：稀薄效应对转子中心的运动轨迹有较大影响，转子系统的不平衡响应表现为周期一运动、周期三运动和周期四运动。

微气体轴承；固定瓦-可倾瓦；转子系统；非线性

随着旋转机械向着高速、微型、无污染的方向发展，微型气体轴承得到了越来越广泛的应用，尤其在MEMS、磁记录系统等领域[1-4]。微型气体轴承本质上是非线性的，目前已有许多对其非线性动力学特性及稳定性进行研究的报道[5-9]。张永芳等[10]针对高速气体轴承-转子系统在主动控制过程中无法避免的时滞现象，建立了含有双时滞控制的三轴向槽动压气体轴承-转子系统的动力学模型，运用微分变换法，得到了非线性气膜压力分布；基于精细积分法，构造了时滞非线性轴承-转子系统动力响应的求解方法；研究了系统的非线性响应，分析了时滞量和反馈控制增益对系统响应的影响。以上研究指出非线性气膜力对转子系统的稳定性有着重要的影响，揭示了微气体轴承-转子系统丰富的非线性动力学特性，但没有对气体的稀薄效应加以考虑。随着轴承尺寸的变小，空气平均分子自由程和轴承半径间隙的比值(即努森数)随之增大，气流在轴承表面产生速度滑移, 对轴承的性能产生较大的影响，因此在微型气体轴承-转子系统的运动分析中，应考虑气体的稀薄效应[11-13]。Zhang等[14]指出努森数在滑移流区和过渡流区的取值范围内，稀薄效应、可压缩性和随机粗糙度三者的耦合作用对气体轴承的分析与设计具有重要的影响。杨琴等[15]采用双向隐式差分法求解了修正的雷诺方程，结合转子运动方程，采用4阶龙格-库塔法计算研究了气体稀薄效应对气体轴承-转子系统不平衡响应的影响。以上研究表明考虑气体的稀薄效应对于微型气体轴承-转子系统的运动分析是非常必要的。

目前，对固定瓦-可倾瓦微气体轴承-转子系统进行非线性运动分析的报道较少。固定瓦-可倾瓦气体轴承[16]在轴承工作时，可倾瓦块可随转子运动状态的改变而绕支点自由摆动，从而使固定瓦-可倾瓦气体轴承具有良好的自对中能力和较高的稳定性。本文以固定瓦-可倾瓦微气体轴承-转子系统为研究对象，考虑了气体稀薄效应的影响，针对基于Burgdorfer一阶滑移速度边界条件的Reynolds方程，利用微分变换法，结合有限差分法，对修正的Reynolds方程进行求解，得到了单块瓦坐标系下的非线性气膜力，然后通过组装技术获得了固定瓦-可倾瓦微型气体轴承的气膜力。在此基础上，对固定瓦-可倾瓦微型气体轴承支撑的刚性Jeffcott转子系统的非线性动力学特性进行了研究，分析了稀薄效应、转速等对转子系统非线性特性的影响，运用转子中心轨迹图、Poincaré映射图和时间历程图分析了转子的不平衡响应。

1 系统方程

图1给出了固定瓦-可倾瓦微气体轴承-对称刚性Jeffcott转子系统的模型示意图，其动力学方程为

(1)

式中：m为转子质量的二分之一；x，y为转子中心在x，y方向的位移分量；fx，fy为轴承在x，y负方向的非线性气膜力分量；g为重力加速度；ex，ey为转子不平衡质量偏心距在x，y方向的分量；fex，fey为不平衡力在x，y方向的分量；ω为转子的角速度。

图1 固定瓦-可倾瓦微气体轴承-对称刚性转子系统示意图

Fig.1 Schematic diagram of a rigid rotor system with micro fixed-tilting pad micro gas-lubricated bearings support

为了便于计算，利用转子角速度ω，轴承的半径间隙c，环境压力pa，轴颈半径R对转子系统的动力学方程进行无量纲化，为此引入以下无量纲变量

将以上无量纲变量代入式(1)中，可得转子系统的无量纲动力学方程为

(2)

式中：M为质量矩阵；X″为加速度向量；W为重量向量；Q为外部激励力向量；F为非线性气膜力向量。其表达式分别为

2 固定瓦-可倾瓦微气体轴承非线性气膜力

2.1 固定瓦-可倾瓦微气体轴承模型

2.2 修正的Reynolds方程及边界条件

考虑气体稀薄效应，由于稀薄气体动力学基本方程(Boltzmann方程)为一阶形式[17]，因此，基于一阶滑移边界的Reynolds方程的无量纲形式为

(3)

式中：φ为沿轴承圆周方向的无量纲坐标；λ为沿轴承轴向方向的无量纲坐标；P为无量纲气膜压力；H为无量纲气膜厚度；Λ为轴承数；Kn0为努森数。其表达式分别为

式中：h为有量纲气膜厚度；μ为气体动力黏度；λa为环境压力下气膜分子的平均自由程。

(a) 固定瓦-可倾瓦微气体轴承端视图

(b) 固定瓦-可倾瓦微气体轴承剖视图

固定瓦-可倾瓦微气体轴承的无量纲压力分布的边界条件可以表示为

式中：φi1、φi2为第i块瓦上气膜的起始角和终止角。

通过对每块瓦气膜区域上的压力积分，便可获得轴瓦径向和切向的无量纲承载力：

2.3 微分变换法求解修正的Reynolds方程

微分变换法具有收敛速度快和计算误差小等特点。为了减少迭代求解非线性气膜力的计算工作量和提高计算精度，本文采用微分变换法求解基于一阶速度滑移边界的Reynolds方程[18]。令Q=P2，I=H2，J=H3，并对式(3)进行微分变换，则式(3)可以表示为

(4)

根据卷积、一阶中差商和二阶中差商的定义，式(4)可以表示为

(5)

式中：i，j分别为轴承周向和轴向方向上的节点编码；k，l，m为微分变换的阶数(k≥l≥m)；Δφ，Δλ，Δτ分别为轴承周向差分步长、轴向差分步长和时间步长。

由于轴颈中心线与轴承的中心线是平行的，所以气膜厚度沿轴承轴向方向没有变化，即∂Hi,j/∂λ=0。当k=0，l=0，m=0时，式(5)可表达为

(6)

当k=0时,有：

(7)

将式(7)代入式(6)，有：

(8)

式中：Hi(0)为上一时刻的气膜厚度；Pi,j(0)为上一时刻的压力分布；Hi(1)，Pi,j(1)分别为求解气膜压力分布的中间变量。Hi(1)的表达式为

(9)

当k=1；l=0, 1；m=1, 0时，式(5)可表示为

(10)

当k=1时,有：

(11)

将式(7)和式(11)代入式(10)，有：

(12)

式中：Hi(2)和Pi,j(2)为求解气膜压力分布的中间变量。Hi(2)的表达式为

(13)

将Hi(0)，Pi,j(0)，Hi(1)代入式(8)可求得Pi,j(1)，再将Hi(0)，Pi,j(0)，Hi(1) ，Pi,j(1) ，Hi(2)代入式(12)即可求得Pi,j(2)。则气膜压力Pi,j可以通过下式求得：

2.4 单块瓦非线性气膜力

为了求解固定瓦-可倾瓦微气体轴承非线性气膜力，首先需要在单瓦坐标系中求解修正的Reynolds方程，得到单瓦坐标系中的非线性气膜力。图3和图4给出了第i块瓦的计算坐标系及非线性气膜力分量的关系图，其中ε为轴颈中心相对于轴承中心的偏心率，εi为轴颈中心相对于第i块瓦的偏心率，εξi和εηi分别是εi在单瓦坐标系ξi和ηi方向上的分量,γi是第i块瓦在单瓦坐标系中的偏位角。

图3 第i块瓦计算坐标

图4 单瓦坐标系中气膜力分量的关系图

根据图3中的几何关系可以得到第i块轴瓦的偏心率εi为

εi=

(14)

式中，Ψ=c/R为轴承的间隙比。

根据轴瓦的平衡条件:Fξi=0,有:

即:

(15)

式中：Pi为第i块瓦上的压力分布；Fθi、Fεi为第i块瓦上的切向和径向气膜力分量，其表达式分别为

在第i块轴瓦坐标系中，气膜厚度为

(16)

在第i块轴瓦坐标系中,相对摆角αi/Ψ的计算采用以下形式:

因为εi<1,cosγi≤1,Ψ<<1,则上式变为

(17)

对于固定轴瓦，气膜力的计算较为简单，不需要修正摆角，直接调用气膜力计算程序就可以得到在轴承坐标系中的非线性气膜力。对于可倾轴瓦，还必须在轴承运行过程中不断修正摆角，以保证可倾轴瓦平衡[16]。摆角的修正是一个迭代过程，迭代步骤为：

(1) 给定预负荷δ、偏心率ε、偏位角θ、轴瓦的支点位置角φki，设定可倾瓦相对摆角的初始值为(αi/Ψ)0；

(2) 将δ，ε，θ，φki和αi/Ψ代入式(14)求得εi；

(3) 将εi，ε，θ，φki和(αi/Ψ)0代入式(17)求得γi0；

(4) 计算第i块瓦气膜的起始角和终止角；

(5) 根据式(15)计算γi；

(6) 根据式(17)修正相对摆角αi/Ψ；

2.5 固定瓦-可倾瓦微气体轴承气膜力的组装

(1) 给定初始参数:预负荷δ，偏位角θ，偏心率ε，可倾瓦的位置角φki和初始相对摆角(αi/Ψ)0；

(2) 若为可倾轴瓦，计算各瓦坐标系ξiOηi下的偏心率εi，偏位角γi。然后调用气膜力计算程序，求得瓦块坐标系下的非线性气膜力Fξi和Fηi，并利用下式将其转化到轴承坐标系xOby下的非线性气膜力FXi和FYi。最后判断瓦块是否平衡，如果不满足则利用牛顿迭代法修正相对摆角。

(3) 若为固定瓦块，则直接调用单块瓦的非线性气膜力计算程序，求得固定瓦在轴承坐标系xOby下的非线性气膜力Fx和Fy。

(4) 气膜力的组装，求得轴承坐标系xOby下的非线性气膜力FX和FY。

3 数值算例

3.1 稀薄效应对转子非线性特性的影响

为了研究稀薄效应对转子非线性特性的影响，选用图1所示的Jeffcott刚性对称转子和图2所示的固定瓦-可倾瓦微气体轴承，对努森数为0(不考虑气体的稀薄效应)和0.03时的转子中心运动轨迹进行了分析。转子和轴承的参数选取为：转子长度为l=0.1 m，圆盘宽度8.25 mm，轴承直径为D=0.006 m，轴承宽度为B=0.006 m，半径间隙为c=3.116E-006 m，支点比为0.5，预负荷系数δ=0.2，转子速度ω=1 000 rad/s。

图5所示为努森数为0时转子中心的周期四运动，图6为努森数为0.03时转子中心的周期一运动。可以看出，考虑稀薄效应和不考虑稀薄效应对转子的非线性特性有较大影响，表明了考虑稀薄效应的必要性。

(a) ω=1 000 rad/s, Kn0=0时，轴颈中心周期四运动轨迹

(b)ω=1 000 rad/s, Kn0=0时，轴颈中心运动轨迹的Poincaé映射在X-Y平面上的投影

图5ω=1 000 rad/s,Kn0=0时，轴颈中心的周期四运动轨迹及其Poincaé映射

Fig.5 The period-4 orbit and Poincaré map of the center of the rotor forω=1 000 rad/s andKn0=0

(a) ω=1 000 rad/s, Kn0=0.03时，轴颈中心周期一运动轨迹

(b) ω=1 000 rad/s, Kn0=0.03时，轴颈中心运动轨迹的Poincaré映射在X-Y平面上的投影

图6ω=1 000 rad/s,Kn0=0.03时，轴颈中心的周期一运动轨迹及其Poincaé映射

Fig.6 The period-1 orbit and Poincaré map of the center of the rotor forω=1 000 rad/s andKn0=0.03

3.2 固定瓦-可倾瓦微气体轴承-刚性转子系统动力学特性计算

对图1所示的固定瓦-可倾瓦微气体轴承-刚性转子系统进行分析，转子和轴承的参数选取如下：努森数Kn0=0.01，转子为对称刚性转子，转子长度l=0.05 m，轴承直径为D=0.006 m，轴承宽度B=0.006 m，半径间隙为c=3.116E-006 m，支点比为0.5，预负荷系数δ=0.2，质量密度为7 800 kg/m3，质量偏心ex=ey=0.000 5 mm，可倾瓦的张角为75°，可倾瓦支点位置角分别为 45°和 315°，固定瓦的张角为170°，气体动力黏度为μ=1.8E-005 Pa·s。

当转速较低时，转子的不平衡响应是稳定的周期运动。图7给出了转子转速为ω=800 rad/s时转子中心的周期运动轨迹、转子中心运动轨迹的Poincaré映射在X-Y平面上的投影、转子中心轨迹在Y方向的时间历程。随着转子转速的增加，系统的响应为周期三运动。图8给出了转子转速为1 000 rad/s时，转子中心的周期三运动轨迹、转子中心运动轨迹的Poincaré映射在X-Y平面上的投影，转子中心运动轨迹在Y方向的时间历程。

(a) ω=800 rad/s, Kn0=0.01时，轴颈中心周期一运动轨迹

(b) ω=800 rad/s, Kn0=0.01时，轴颈中心运动轨迹的Poincaré映射在X-Y平面上的投影

(c) ω=800 rad/s，Kn0=0.01时，转子中心Y方向上的时间历程

图7ω=800 rad/s,Kn0=0.01时，轴颈中心的周期一运动轨迹、Poincaé映射和时间历程

Fig.7 The period-1 orbit，Poincaré map and time series of the center of the rotor forω=800 rad/s andKn0=0.01

4 结 论

本文针对固定瓦-可倾瓦微气体轴承支撑的刚性Jeffcott转子系统，研究了气体稀薄效应和转速对系统非线性动力学行为的影响，结论如下：

(1) 运用微分变换法求解固定瓦-可倾瓦微气体轴承的非线性气膜压力分布，在提高计算精度的同时，减少了迭代求解气体轴承-转子系统非线性动力响应的计算工作量。

(a) ω=1 000 rad/s，Kn0=0.01时，轴颈中心周期三运动轨迹

(b)ω=1 000 rad/s，Kn0=0.01时，轴颈中心运动轨迹的Poincaré映射在X-Y平面上的投影

(c)ω=1 000 rad/s，Kn0=0.01时，转子中心Y方向上的时间历程

图8ω=1 000 rad/s,Kn0=0.01时，轴颈中心的周期三运动轨迹、Poincaé映射和时间历程

Fig.8 The period-3 orbit，Poincaré map and time series of the center of the rotor forω=1 000 rad/s andKn0=0.01

(2) 气体稀薄效应对转子中心运动轨迹有较大影响，在微气体轴承的研究中，气体稀薄效应的作用不可忽视。

(3) 运用转子中心轨迹图、Poincaré映射图和时间历程图分析了转子的非线性动力学行为，当转速较低时，转子的不平衡响应为稳定的周期运动，但转速增加时，转子的不平衡响应为周期三运动，周期三运动是通向混沌的典型道路。

本文的研究工作可为固定瓦-可倾瓦微气体轴承-转子系统的动力学设计提供理论技术参考，具有积极的指导意义。

[1] DESSORNES O，LANDAIS S，VALLE R，et al. Advances in the development of a microturbine engine[J]. ASME Journal of Engineering for Gas Turbines and Power，2014，136(7)：071201-1-071201-9.

[2] DEMIRI S，BOEDO S，HOLSEN L S. Wear characteristics of large aspect ratio silicon microbearing systems[J]. Wear，2014，312(1/2)：58-69.

[3] ISE T，ARITA N，ASAMI T，et al. Development of externally pressurized small-size conical-shaped gas bearings for micro rotary machines[J]. Precision Engineering，2014，38(3)：506-511.

[4] DAI X Y，LI H，SHEN S N，et al. Numerical simulation of bearing force over bit-patterned media using 3-D DSMC method[J]. IEEE Transactions on Magnetics，2015，51(11)：7209804.

[5] WANG W Z，LIU Y Z，JIANG P N. Numerical investigation on influence of real gas properties on nonlinear behavior of labyrinth seal-rotor system[J]. Applied Mathematics and Computation，2015，263：12-24.

[6] ZHANG Y F，HEI D，LU Y J，et al. Bifurcation and chaos analysis of nonlinear rotor system with axial-grooved gas-lubricated journal bearing support[J]. Chinese Journal of Mechanical Engineering，2014，27(2)：358-368.

[7] MEYBODI R R，MOHAMMADI A K，BAKHTIARI-NEJAD F. Numerical analysis of a rigid rotor supported by aerodynamic four-lobe journal bearing system with mass unbalance[J]. Communications in Nonlinear Science and Numerical Simulation，2012，17(1)：454-471.

[8] WANG C C. Bifurcation and nonlinear analysis of a flexible rotor supported by a relative short spherical gas bearing system[J]. Communications in Nonlinear Science and Numerical Simulation，2010，15(9)：2659-2671.

[9] HASSINI M A，ARGHIR M. A simplified and consistent nonlinear transient analysis method for gas bearing：extension to flexible rotors[J]. ASME Journal of Engineering for Gas Turbines and Power，2015，137(9)：092502-1-092511-10.

[10] ZHANG Y F，ZHANG S，LIU F X，et al. Motion analysis of a rotor supported by self-acting axial groove gas bearing system with double time delays[J]. Proceedings of the Institution of Mechanical Engineers Part C: Journal of Mechanical Engineering Science，2014，228(16)：2888-2899.

[11] WHITE J. Combined effects of surface roughness and rarefaction in the region between high wave number-limited and high bearing number-limited lubricant flows[J]. ASME Journal of Tribology, 2015，137(1)：012001-1-012010-10.

[12] ZHANG X Q，WANG X L, LIU R, et al. Modeling and analysis of micro hybrid gas spiral-grooved thrust bearing for microengine[J]. ASME Journal of Engineering for Gas Turbines and Power, 2013，135(12)：122508-1-122515-8.

[13] ZHANG X Q，WANG X L，ZHANG Y Y. Non-linear dynamic analysis of the ultra-short micro gas journal bearing-rotor systems considering viscous friction effects[J]. Nonlinear Dynamics，2013，73(1)：751-765.

[14] ZHANG W M，MENG G，PENG Z K. Gaseous slip flow in micro-bearings with random rough surface[J]. International Journal of Mechanical Sciences，2013，68：105-113.

[15] 杨琴，刘宇陆，张海军，等. 气体稀薄效应对微机电系统(MEMS)气体轴承-转子系统不平衡响应的影响[J]. 西安交通大学学报，2015，49(7)：134-139.

YANG Qin，LIU Yulu，ZHANG Haijun，et al. Influence of gas rarefaction effect on unbalance response of micro-electro-mechanical (MEMS) gas bearing-rotor system[J]. Journal of Xi’an Jiaotong University，2015，49(7)：134-139.

[16] LU Y J，ZHANG Y F，SHI X L，et al. Nonlinear dynamic analysis of a rotor system with fixed-tilting-pad self-acting gas-lubricated bearings support[J]. Nonlinear Dynamics，2012，69(3)：877-890.

[17] SHEN S，CHEN G，CRONE R M, et al. A kinetic-theory based first order slip boundary condition for gas flow[J]. Physics of Fluids，2007，19(8)：086101-1-086106-6.

[18] 赵家奎. 微分变换及其在电路中的应用[M]. 武汉：华中理工大学出版社，1988.

Nonlinear dynamic analysis for a rotor system supported with fixed-tilting pad micro gas-lubricated bearings

ZHANG Yongfang1,2, LÜ Yedi1,2, XIAO Liangjun1,2, ZHAO Jingqun1,2, LIU Cheng1,2

(1.School of Printing, Packaging Engineering and Digital Media Technology, Xi’an University of Technology, Xi’an 710048, China；2.State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University, Xi’an 710049, China)

Aiming at fixed-tilting pad micro gas-lubricated bearings, considering gas rarefaction effect, Reynolds equation with Burgdorfer 1st-order slip velocity boundary was derived. This Reynolds equation was solved with the differential transformation method combined with the finite difference one, and then the nonlinear gas film force of each pad was calculated. Adopting the assembly technique, the gas film force of the fixed-tilting pad micro gas-lubricated bearing was obtained. For a rigid Jeffcott rotor system supported with fixed-tilting pad micro gas-lubricated bearings, its unbalance responses were analyzed using rotor center’s orbit diagrams, Poincaré map diagrams and time history diagrams. In addition, influences of Knudsen number and rotating speed on the nonlinear dynamic characteristics of the rotor were compared. The results showed that the gas rarefaction effect has a larger influence on orbits of the rotor’s centers; the unbalance responses of the rotor system are characterized as period-1, period-3, and period-4 motions.

micro gas-lubricated bearing; fixed-tilting pad; rotor system; nonlinear

国家自然科学基金(51505375)；机械结构强度与振动国家重点实验室开放课题(SV2016-KF-10)；陕西省教育厅科学研究计划项目(15JK1549,15JS068)

2016-01-30 修改稿收到日期：2016-04-16

张永芳 女，博士，副教授，1975年生

TH113

A

10.13465/j.cnki.jvs.2017.11.010