有限長軸承支撐的轉子系統(tǒng)非線性動力學分析

黑棣1,2呂延軍2張永芳2

1.陜西鐵路工程職業(yè)技術學院，渭南，7140002.西安理工大學，西安，710048

摘要：基于π油膜假設，運用變分約束原理和分離變量法求得了有限長動壓軸承非線性油膜力的近似解析表達式。通過計算比較可知，該方法和有限元法的計算結果非常接近，而且可以節(jié)約大量計算時間。在此基礎上，以轉子的角速度、圓盤處的偏心量及轉軸的剛度為控制參數(shù)，運用Newmark法分析了軸承-轉子系統(tǒng)的非線性動力學行為。數(shù)值結果揭示系統(tǒng)具有周期運動、倍周期運動、準周期運動、六周期運動、十周期運動等復雜豐富的非線性動力學特征。

關鍵詞：非線性動力學；軸承-轉子系統(tǒng)；有限長軸承；穩(wěn)定性；分岔

中圖分類號：TH133

收稿日期：2015-01-04

基金項目：國家自然科學基金資助項目(51375380)；陜西鐵路工程職業(yè)技術學院常規(guī)項目(2013-19)

作者簡介：黑棣，男，1983年生。陜西鐵路工程職業(yè)技術學院機電系講師，西安理工大學機械與精密儀器工程學院博士研究生。主要研究方向為非線性動力學及控制、潤滑理論及新型軸承技術。發(fā)表論文10余篇。呂延軍，男，1972年生。西安理工大學機械與精密儀器工程學院教授、博士研究生導師。張永芳，女，1975年生。西安理工大學印刷與包裝工程學院副教授。

Analysis of Nonlinear Dynamics of Rotor System Supported by Finite Long Bearing

Hei Di1,2Lü Yanjun2Zhang Yongfang2

1.Shaanxi Railway Institute, Weinan，Shaanxi，714000

2.Xi’an University of Technology，Xi’an，710048

Abstract:Based on the assumption of π oil film, an approximation expression of nonlinear oil film force of finite hydrodynamic bearing was obtained by variational constraint principal and the separation of variables. The results calculated by the proposed method are in good agreement with the oil film forces by the finite element method, and the computing cost is reduced greatly. Based on the nonlinear oil film force, the speed of rotor, eccentricity of disk and stiffness of shaft were taken as control parameters, nonlinear dynamic behaviors of bearing-rotor system were analyzed by Newmark method. Numerical results reveal periodic, double-periodic, quasi-periodic, 6-periodic, 10-periodic of rich and complex nonlinear behaviors of the system.

Key words: nonlinear dynamics；bearing-rotor system；finite long bearing；stability；bifurcation

0引言

大型旋轉機械廣泛用于燃氣輪機、航空發(fā)動機、工業(yè)壓縮機及各種電動機等機械裝置中，是國家基礎設施和基礎工業(yè)中最關鍵、最核心的設備之一，在電力、航空、冶金、石化、紡織、機械、動力、運輸?shù)刃袠I(yè)中有著舉足輕重的地位。軸承-轉子系統(tǒng)是旋轉機械的核心部件，且在高速、精密旋轉機械中的應用日益廣泛，所以對軸承-轉子系統(tǒng)非線性動力學特性的研究十分重要。

軸承-轉子系統(tǒng)的非線性動力學行為和穩(wěn)定性的研究主要是基于無限短軸承模型[1-6]和無限長軸承模型[7-10]。當軸承的長徑比非常小的時候(無限短軸承假設)，油膜力沿軸承周向的分布就可以忽略。相反，當軸承的長徑比非常大的時候(無限長軸承假設)，油膜力沿軸承軸向的分布也可以忽略。

無限短和無限長軸承模型是非常簡單的模型，并不能反映實際的情況。為了更準確地計算油膜力，學者們提出了許多計算油膜力的方法。Zheng等[11]提出了一個里茲模型來計算油膜力，這種方法可以節(jié)約大量的計算時間。Xiao等[12]基于變分不等式理論提出了一個快速有效的計算油膜力的算法，此算法將油膜力及其Jacobis的計算轉化為求解一組線性代數(shù)方程?；诙鄥?shù)法，Vignolo等[13]給出了無限長軸承假設下油膜力的近似解析解，并運用此方法求得了靜態(tài)油膜力。Hirani等[14]、Bastani等[15]運用不同的無限短和無限長軸承模型的組合來近似求解有限長軸承油膜力。Wang等[16]基于動態(tài)π油膜假設得到了有限長軸承油膜力的解析解。

本文運用變分原理和分離變量法求得了有限長軸承油膜力的近似解析解，即在求解雷諾方程的時候，將雷諾方程中的壓力分布函數(shù)看作是2個獨立函數(shù)的乘積(一個函數(shù)僅僅是軸承周向坐標的函數(shù)，另一個函數(shù)僅僅是軸承軸向坐標的函數(shù))。在求得有限長軸承油膜力的近似解析解之后，針對柔性對稱轉子，分析了此軸承-轉子系統(tǒng)的非線性動力學行為。

1轉子系統(tǒng)的動力學方程

有限長滑動軸承支撐的對稱柔性轉子模型如圖1所示。

圖1　有限長滑動軸承—柔性轉子系統(tǒng)示意圖

有限長滑動軸承支承的對稱柔性轉子系統(tǒng)動力學方程可寫為

(1)

2ma+mb=m

式中，M、K分別為轉子的質量矩陣和剛度矩陣；f為軸承的非線性油膜力；Q為施加在轉子上的周期激勵力(與角速度同步的不平衡激勵或氣動力激勵)；x為轉子的位移；m為轉子總質量；ma、mb分別為集中到2個軸頸及圓盤處的質量；E、I分別為轉子的彈性模量及轉子截面的赤道慣性矩；l為兩個軸承的中心距；ω為轉子角速度；g為重力加速度；eji(i=x，y；j=a，b)為軸承及圓盤處不平衡質量偏心距在x、y方向的分量；fxa、fya分別為2個軸承在x軸和y軸負方向的非線性油膜力分量。

對于軸承-轉子系統(tǒng)，為了理論推導和計算方便，引入如下量綱一變量：

式中，X、Y分別為軸頸中心x、y方向的位移分量；e為偏心量；c為軸承半徑間隙；m′為轉子質量的一半；σ為Sommerfeld數(shù)(Sommerfeld數(shù)是表征系統(tǒng)載荷的參數(shù))；B為軸承寬度；η為潤滑油動力黏度；R為軸承半徑。

為了研究方便，引入系統(tǒng)參數(shù)：

該參數(shù)對于具體的軸承-轉子系統(tǒng)，是一個常量。

2有限長軸承的油膜力計算

2.1油膜力的計算方法

圖2給出了有限長動壓軸承示意圖及其計算坐標。圖2中，Ob為軸承中心，Oj為軸頸中心，Φ為從y軸負方向順時針轉至油膜位置的角度，θ為偏位角，φ為從OjOb延長線順時針方向轉至油膜位置的角度，fr、ft分別為作用在軸頸上的非線性油膜力徑向和切向分量，fx、fy分別為作用在軸頸上的非線性油膜力在x軸和y軸負方向的分量，h為油膜厚度，r為軸頸半徑，W為軸承載荷。

圖2　有限長軸承圖

對于徑向滑動軸承，將潤滑油膜視為不可壓縮的流體，則其一般的Reynolds方程為

(2)

h=c+αcosφ

式中，Gφ、Gz為紊流因子；p為油膜的壓力；γ為軸承周向坐標；z為軸承軸向坐標；U為軸頸的周向速度；α為偏心距。

引入以下量綱一變量：

λ=2z/Bε=α/cτ=ωtP=p/p0

p0=2ωη/ψ2H=1+εcosφ

式中，λ為量綱一軸向坐標；τ為量綱一時間；H為量綱一油膜厚度；ε為偏心率；ψ為間隙比；P為量綱一的油膜壓力。

將式(2)代入式(1)可得

(3)

本文基于π油膜假設，運用變分原理和分離變量法對式(3)進行求解，求得非線性油膜壓力分布的近似解析解。求解過程中，所用到的邊界條件：油膜起始邊處的油膜壓力為0，即φ=0，P=0；油膜終止邊處的油膜壓力為0，即φ=π，P=0；油膜在終止邊處連續(xù)，即φ=π，?P/?φ=0。

為了表達方便，將式(3)改寫為

E(P)=f

(4)

式(4)定義在平面(φ，λ)的有界區(qū)域Ω上，邊界為Γ，設H1(Ω)是Soblev空間，B(u，v)：H1(Ω)×H1(Ω)→R是強制對稱連續(xù)雙線性泛函，其中

(5)

u,v∈H1(Ω)

f(v)是H1(Ω)上的線性連續(xù)泛函：

f(v)=?ΩfvdΩ

(6)

定義二次泛函：

J(v)=B(u，v)/2-f(v)

(7)

根據(jù)潤滑力學第一變分原理，如果P是式(4)的解，那么P也同時滿足如下泛函的極值問題：

(8)

其中，檢驗函數(shù)集K是H1(Ω)中的非空閉凸集。

油膜壓力只存在于油膜收斂區(qū)，則式(8)只在0≤φ≤π和-1≤λ≤1范圍內(nèi)成立。假設非線性油膜壓力試驗函數(shù)是2個獨立函數(shù)的乘積即P(φ，λ)=P*(φ)ζ(λ)，其中，P*(φ)是關于φ的未知函數(shù)，ζ(λ)是λ的未知函數(shù)。要想求得油膜壓力試驗函數(shù)P(φ，λ)，則需分別求解2個獨立函數(shù)P*(φ)和ζ(λ)。

將試驗函數(shù)P(φ，λ)=P*(φ)ζ(λ)代入式(7)可得

(9)

對式(9)取一階變分，由極值條件得

(10)

由于P*(φ)符合一階變分條件，故有W=-K，則式(10)可寫為

(11)

由式(11)即可求解出

(12)

通過以上方法可以求出獨立函數(shù)ζ(λ)的表達式，由于非線性油膜壓力試驗函數(shù)P(φ，λ)=P*(φ)ζ(λ)，所以要想得到油膜壓力試驗函數(shù)P(φ，λ)的表達式，則還需求出獨立函數(shù)P*(φ)的表達式?；跓o限長軸承假設，P*(φ)的表達式可以通過以下的方法求出。

無限長軸承是基于軸承寬度B遠大于軸承直徑d這個假設，即B?d。當B?d時，油膜力沿軸向的變化率較其沿周向的變化率可以忽略不計(?p/?z??p/?φ)，即λ方向的壓力梯度可忽略，也就是說無限長軸承的油膜壓力函數(shù)僅僅是關于軸承周向坐標φ的函數(shù)。由于P*(φ)也僅僅是關于軸承周向坐標φ的函數(shù)，所以可以將無限長軸承油膜壓力分布函數(shù)作為P*(φ)。

基于無限長軸承假設，式(3)可以寫為

(13)

對式(13)積分兩次，即可求得

(14)

通過對P(φ，λ)積分，我們可以求出徑向、切向的量綱一油膜力Fr和Ft：

(15)

由Fr、Ft可以求出x、y方向的油膜力分量Fx和Fy：

(16)

2.2油膜力的計算結果

2.1節(jié)提出了求解有限長軸承油膜力近似解析解的方法，為了驗證此方法的可行性，將本文提出方法的計算結果、無限長軸承模型的計算結果和FEM計算結果進行比較。對比參數(shù)如下：偏心率ε=0.5，x、y方向速度分別為vx=0.01，vy=0，有限長滑動軸承寬徑比B/d=1.0。圖3給出了3種方法的計算結果。

(a)三種方法計算的F x

(b)三種方法計算的F y 圖3　三種方法計算的油膜力的比較

從圖3可以看出，本文提出方法的計算結果與FEM的計算結果吻合得非常好，無限長軸承模型計算的結果與FEM的計算結果相差較大。

圖4給出了本文提出方法的計算速度和FEM的計算速度。從圖4可以看出，本文提出方法的計算速度遠高于FEM的計算速度，即本文提出的方法可以節(jié)約大量的計算時間。

圖4　本文方法與FEM方法的計算速度

3數(shù)值算例

3.1以量綱一角速度為控制參數(shù)

(a)軸承處

(b)圓盤處 圖5　不同轉速下轉子軸承處和圓盤處y方向的分岔圖

(a)周期軌跡(b)Poincaré映射點列 圖6　轉子軸承處的周期軌跡和 Poincaré映射點列( =0.95)

(a)準周期軌跡

(b)Poincaré映射點列 圖7　轉子軸承處的準周期軌跡和 Poincaré映射點列( =1.05)

(a)十周期軌跡

(b)Poincaré映射點列 圖8　轉子軸承處的十周期軌跡和 Poincaré映射點列( =1.10)

3.2以轉軸量綱一剛度為控制參數(shù)

圖9為軸承處和圓盤處轉子中心y方向的位移隨轉軸量綱一剛度變化的分岔圖。由圖9可以看出，轉子的運動非常復雜，當轉軸的量綱一剛度從2變化到6時，轉子所經(jīng)歷的運動分別為周期運動、準周期運動、周期運動、倍周期運動。由此可以看出，剛度對轉子系統(tǒng)的運動影響較大。圖10所示為轉軸的量綱一剛度k=2.5時，轉子軸承處轉子的運動軌跡和Poincaré映射點列。隨著轉軸剛度的增加，系統(tǒng)將發(fā)生擬周期分岔。圖11為k=3.7時，轉子軸承處準周運動軌跡、Poincaré映射圖。隨著轉軸剛度繼續(xù)增加，系統(tǒng)將發(fā)生倒分岔，即由準周期運動變?yōu)橹芷谶\動。當轉軸剛度再繼續(xù)增加，系統(tǒng)將會發(fā)生倍周期分岔，即由原來的周期運動轉化為倍周期運動。圖12為k=5.2時，軸承處轉子的倍周期運動軌跡、Poincaré映射圖。從圖12b可以看出，倍周期軌跡在穿越Pioncaré截面時為2個不動點。轉軸的剛度繼續(xù)增加，系統(tǒng)將再次發(fā)生倒分岔，即由原來的倍周期運動變?yōu)橹芷谶\動。

(a)軸承處

(b)圓盤處 圖9　不同剛度下轉子軸承處和圓盤處y方向的分岔圖

(a)周期軌跡(b)Poincaré映射點列 圖10　轉子軸承處的周期軌跡和 Poincaré映射點列(k=2.5)

(a)準周期軌跡

(b)Poincaré映射點列 圖11　轉子軸承處的準周期軌跡和 Poincaré映射點列(k=3.7)

(a)二周期軌跡

(b)Poincaré映射點列 圖12　轉子軸承處的二周期軌跡和 Poincaré映射點列(k=5.2)

3.3以圓盤處量綱一偏心量為控制參數(shù)

(a)軸承處

(b)圓盤處 圖13　轉子軸承處和圓盤處y方向的分岔圖

(a)準周期軌跡

(b)Poincaré映射點列 圖14　轉子軸承處的準周期軌跡和 Poincaré映射點列( =0.01)

(a)六周期軌跡

(b)Poincaré映射點列 圖15　轉子軸承處的六周期軌跡和 Poincaré映射點列( =0.27)

4結語

首先，本文基于有限長軸承模型假設，運用變分原理和分離變量法求得了有限長軸承的近似解析解。在求解雷諾方程的時候，將油膜壓力分布函數(shù)P(φ，λ)看成是兩個獨立函數(shù)P*(φ)、ζ(λ)的乘積，其中，P*(φ)是關于φ的周向函數(shù)，ζ(λ)是關于λ的軸向函數(shù)。周向函數(shù)P*(φ)用無限長軸承假設下的油膜壓力分布來代替，軸向函數(shù)ζ(λ)通過變分原理求得。在求得P*(φ)和ζ(λ)之后便可以求出有限長軸承油膜力的近似解。將本文提出方法的計算結果與有限元和無限長軸承模型的計算結果進行比較可知，本文提出的方法與有限元的計算結果非常接近，而且計算速度要比有限元方法的計算速度快很多，即本文提出的油膜力計算方法不僅能夠保證計算的精確度，而且還能夠提高計算效率，節(jié)約大量計算時間。

在求得有限長軸承油膜力近似解的基礎上，分別以轉軸的角速度、轉軸的剛度和圓盤的偏心量為控制參數(shù)，運用軌跡圖、Poincaré映射圖、時間歷程和分岔圖分析了對稱柔性轉子-軸承系統(tǒng)的非線性動力學行為。從計算的結果可以看出，轉子的非線性運動行為是非常豐富復雜的，呈現(xiàn)出周期、倍周期、準周期、六周期、十周期等非線性動力學現(xiàn)象。

參考文獻：

[1]Wang Ligang，Cao Dengqin，Huang Wenhu. Nonlinear Coupled Dynamics of Flexible Blade-rotor-bearing Systems[J]. Tribology International，2010，43(4)：759-778.

[2]Castroa H F，Cavalca K L，Nordmann R. Whirl and Whip Instabilities in Rotor-bearing System Considering a Nnonlinear Force Model[J]. Journal of Sound and Vibration，2008，317(1/2)：273-293.

[3]Xie Wenhui, Tang Yougang, Chen Yushu. Analysis of Motion Stability of the Flexible Rotor-bearing System with Two Unbalanced Disks[J]. Journal of Sound and Vibration, 2008, 310(1/2): 381-393.

[4]Jing Jianping, Meng Guang, Sun Yi, et al. On the Non-linear Dynamic Behavior of a Rotor-bearing System[J]. Journal of Sound and Vibration, 2004, 274(3/5): 1031-1044.

[5]Wang Junguo, Zhou Jianzhong, Dong Dawei, et al. Nonlinear Dynamic Analysis of a Rub-impact Rotor Supported by Oil Film Bearings[J]. Archive of Applied Mechanics, 2013, 83(3): 413-430.

[6]Wang Jianke, Khonsari M M. Bifurcation Analysis of a Flexible Rotor Supported by Two Fluid-film Journal Bearings[J]. ASME Journal of Tribology, 2006, 128(3): 594-603.

[7]Li Dexin, Xu Jianxue. A Method to Determine the Periodic Solution of the Non-linear Dynamics System[J]. Journal of Sound and Vibration, 2004, 275(1/2): 1-16.

[8]Chang Jian，CaiWan. Non-linear Dynamic Analysis of Dual Flexible Rotors Supported by Long Journal Bearings[J]. Mechanism and Machine Theory, 2010, 45(6): 844-866.

[9]Wang Jianke, Khonsari M M. Effects of Oil Inlet Pressure and Inlet Position of Axially Grooved Infinitely Long Journal Bearings, Part I: Analytical Solutions and Static Performance[J]. Tribology International, 2008, 41(2): 119-131.

[10]Wang Jianke, Khonsari M M. Effects of Oil Inlet Pressure and Inlet Position of Axially Grooved Infinitely Long Journal Bearings, Part II: Nonlinear Instability Analysis[J]. Tribology International, 2008, 41(2): 132-140.

[11]Zheng Tiesheng, Yang Shuhua, Xiao Zhonghui, et al. A Ritz Model of Unsteady Oil-film Forces for Nonlinear Dynamic Rotor-bearing System[J]. ASME Journal of Applied Mechanics, 2004, 71(2): 219-224.

[12]Xiao Zhonghui, Wang Liping, Zheng Tiesheng. An Efficient Algorithm for Fluid Force and Its Jacobian Matrix in Journal Bearing[J]. ASME Journal of Triboogy, 2006, 128(2): 291-295.

[13]Vignolo G G, Barilá D O, Quinzani L M. Approximate Analytical Solution to Reynolds Equation for Finite Length Journal Bearings[J]. Tribology International, 2011, 44(10): 1089-1099.

[14]Hirani H, Athre K, Biswas S. Dynamically Loaded Finite Length Journal Bearings: Analytical Method of Solution[J]. ASME Journal of Tribology, 1999, 121(4): 844-852.

[15]Bastani Y, Queiroz M D. A New Analytic Approximation for the Hydrodynamic Forces in Finite-length Journal Bearings[J]. ASME Journal of Tribology, 2010, 132(1): 014502-1-9.

[16]Wang Yongliang, Liu Zhansheng, Kang Weijia, et al. Approximate Analytical Model for Fluid Film Force of Finite Length Plain Journal Bearing[J]. IMechE, Part C: Journal of Mechanical Engineering Science, 2012, 226(5): 1345-1355.

(編輯張洋)