第一作者毛文貴女，博士生，副教授，1975年1月生

通信作者韓旭男，博士,教授，博士生導(dǎo)師，1968年6月生

滑動(dòng)軸承-轉(zhuǎn)子系統(tǒng)Riccati-Newmark加速度傳遞矩陣法

毛文貴1,2，韓旭1，劉桂萍1

(1. 湖南大學(xué)機(jī)械與運(yùn)載工程學(xué)院,長(zhǎng)沙410082；2. 湖南工程學(xué)院機(jī)械工程學(xué)院,湖南湘潭411101)

摘要：為克服傳遞矩陣法數(shù)值不穩(wěn)定及非線性瞬態(tài)響應(yīng)分析中滑動(dòng)軸承矩陣建立困難問題，提出Riccati-Newmark加速度傳遞矩陣法。借助Newmark加速度法建立傳遞矩陣，采用Taylor級(jí)數(shù)預(yù)估滑動(dòng)軸承軸心下一時(shí)刻位移、速度建立滑動(dòng)軸承矩陣；據(jù)邊界條件用Riccati傳遞矩陣法求滑動(dòng)軸承-轉(zhuǎn)子系統(tǒng)非線性瞬態(tài)響應(yīng)，提高數(shù)值穩(wěn)定性。以單圓盤轉(zhuǎn)子系統(tǒng)為例，與傳統(tǒng)軸頸下一時(shí)刻位移、速度近似線性擾動(dòng)處理的瞬態(tài)響應(yīng)對(duì)比分析，驗(yàn)證此方法的有效性；討論不同轉(zhuǎn)速下線性、非線性油膜力的瞬態(tài)軌跡。

關(guān)鍵詞：非線性轉(zhuǎn)子系統(tǒng)；瞬態(tài)響應(yīng)；Newmark加速度法；Riccati傳遞矩陣法

基金項(xiàng)目：國(guó)家自然科學(xué)基金11202073;高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金(20130161130001);“高檔數(shù)控機(jī)床與基礎(chǔ)制造裝備”科技重大專項(xiàng)2012(2012ZX04002-091)

收稿日期：2014-08-07修改稿收到日期：2014-09-25

中圖分類號(hào):TH113 .3,1;O322

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2015.20.014

Abstract:In order to eliminate the numerical instability of transfer matrix method and to build a transfer matrix of nonlinear elements( bearings), the Riccati transfer matrix method was extended to the transient analysis of nonlinear rotor-bearing systems. In the method, the transfer matrix was obtained with the aid of the Newmark acceleration formulation, the deflections and velocities at the stations containing nonlinear element (bearings) were predicted by Taylor series, and then the deflections, velocities and accelerations at all stations were solved by using the Riccati transfer matrix method combined with the Newmark acceleration formulation integration according to the boundary conditions. An example of single disc rotor system was presented and the results were compared with those by a transient analysis considering linear perturbation and linear oil film force to verify the effectiveness of the proposed method.

Riccati transfer matrix method combined with Newmark acceleration formulation integration for analysing sliding bearings and rotor system

MAOWen-gui1,2,HANXu1,LIUGui-ping1(1. HNU College of Mechanical and Vehicle Engineering, Hunan University, Changsha 410082, China；2. College of Mechanical Engineering, Hunan Institute of Engineering, Xiangtan 411101, China)

Key words:nonlinear rotor system; transient response; Newmark acceleration method; Riccati transfer matrix method

滑動(dòng)軸承作為彈性支承可降低軸系結(jié)構(gòu)剛度，使軸系在較寬速度、載荷范圍內(nèi)無磨損工作，具有高回轉(zhuǎn)精度。軸系響應(yīng)作為研究轉(zhuǎn)子動(dòng)力學(xué)的主要內(nèi)容，其準(zhǔn)確性較大程度上取決于油膜力處理。由于油膜力的強(qiáng)非線性，滑動(dòng)軸承支承軸系為非線性系統(tǒng)。傳遞矩陣法[1-2]具有矩陣階數(shù)不隨軸系自由度增加而增加、編程簡(jiǎn)單等優(yōu)點(diǎn)，最適用轉(zhuǎn)子類鏈?zhǔn)较到y(tǒng)，多用其分析研究線性轉(zhuǎn)子系統(tǒng)動(dòng)力響應(yīng)問題，而對(duì)非線性轉(zhuǎn)子系統(tǒng)動(dòng)力響應(yīng)分析研究有限。對(duì)非線性轉(zhuǎn)子系統(tǒng)，滑動(dòng)軸承傳遞矩陣只有獲得非線性油膜力方能真正建立，非線性油膜力與軸心位移、速度有關(guān)，而非線性系統(tǒng)動(dòng)力響應(yīng)未知，因此獲得滑動(dòng)軸承傳遞矩陣較困難。加之滑動(dòng)軸承強(qiáng)非線性精確表達(dá)式不易獲得，從而阻礙滑動(dòng)軸承-轉(zhuǎn)子系統(tǒng)非線性研究進(jìn)展。轉(zhuǎn)子系統(tǒng)瞬態(tài)響應(yīng)分析中關(guān)于軸承油膜力處理有兩種方式：線性擾動(dòng)即對(duì)下一時(shí)刻軸頸位移用上一時(shí)刻位移代替；將線性化油膜力用4個(gè)剛度、阻尼系數(shù)線性化表示。李蘋等[3]采用Riccati傳遞矩陣與Newmark數(shù)值積分法研究滑動(dòng)軸承-轉(zhuǎn)子系統(tǒng)非線性瞬態(tài)響應(yīng)?？紤]油膜非線性，通過CFD軟件模擬[4-5]，但對(duì)軸心位移、速度用線性擾動(dòng)近似。線性油膜力、線性擾動(dòng)均不能真實(shí)反映軸系強(qiáng)非線性。用傳遞矩陣法研究軸系非線性瞬態(tài)響應(yīng)時(shí)傳遞矩陣中含材料、結(jié)構(gòu)及時(shí)間等參數(shù)，由于一般軸系中材料彈性模量加大而瞬態(tài)響應(yīng)分析中時(shí)間間隔較小，采用力、位移狀態(tài)變量表示傳遞矩陣易出現(xiàn)較0(1)量級(jí)大得多、小得多的數(shù)，會(huì)引起新病態(tài)矩陣[6]。

對(duì)此，本文提出非線性擾動(dòng)(Taylor級(jí)數(shù)預(yù)估)及非線性油膜力的Riccati-Newmark加速度傳遞矩陣法。用力、加速度狀態(tài)變量避免數(shù)值不穩(wěn)定性[7]。為獲得非線性系統(tǒng)準(zhǔn)確的油膜響應(yīng)，通過Taylor級(jí)數(shù)預(yù)估滑動(dòng)軸承非線性元件下一時(shí)刻位移、速度，借助Newmark加速度法建立傳遞矩陣，用Riccati傳遞矩陣法求滑動(dòng)軸承-轉(zhuǎn)子系統(tǒng)非線性瞬態(tài)響應(yīng)，獲得滑動(dòng)軸承下一時(shí)刻位移及速度,通過精度控制誤差獲得滑動(dòng)軸承位移、速度。以單圓盤彈性轉(zhuǎn)子系統(tǒng)為例，分別采用線性油膜力、非線性油膜力，線性擾動(dòng)、非線性擾動(dòng)對(duì)轉(zhuǎn)子系統(tǒng)瞬態(tài)響應(yīng)比較分析，驗(yàn)證本文方法的有效性。

1Newmark加速度法

求解非線性運(yùn)動(dòng)微分方程通常用數(shù)值解法，如Runge-Kutta法、Newmark法[8]等。前者收斂速度慢，易發(fā)散且積分步長(zhǎng)對(duì)結(jié)果精度及正確性影響較大。后者當(dāng)參數(shù)選擇合適時(shí)積分無條件穩(wěn)定，結(jié)果與步長(zhǎng)大小無關(guān)。本文用后者求解非線性瞬態(tài)響應(yīng)。

對(duì)已知上一步位移、速度、加速度及求解下一步位移、速度、加速度Newmark法進(jìn)行變換，獲得Newmark加速度法，即

(1)

2改進(jìn)的軸系各元件傳遞矩陣

傳遞矩陣法分析軸系結(jié)構(gòu)時(shí)將軸系離散為無質(zhì)量軸段與集總圓盤(F=0)、有支承集總圓盤(F≠0)等單元，見圖1。模塊中m，Jp，JD為集總圓盤節(jié)點(diǎn)處質(zhì)量、極轉(zhuǎn)動(dòng)慣量、直徑轉(zhuǎn)動(dòng)慣量；D，d，L為無質(zhì)量軸段外徑、內(nèi)徑及跨度。

圖1　轉(zhuǎn)子軸承系統(tǒng)模塊化圖 Fig.1 bearings-rotor systems modular figure

每種單元均有傳遞矩陣。Riccati傳遞矩陣法為提高傳遞矩陣數(shù)值穩(wěn)定性，將狀態(tài)向量中元素分為兩部分，即一部分由對(duì)應(yīng)初始截面狀態(tài)向量中具有零值元素組成，另部份由其余互補(bǔ)元素組成。狀態(tài)向量為

q=(My,Qy,Mz,Qz?θy，Y,θz,Z)T=(f?e)T

式中：Mz，My為彎矩；Qz，Qy為剪力；Z，Y為位移；θz，θy為轉(zhuǎn)角。

用力、位移狀態(tài)變量傳遞矩陣法研究軸系非線性瞬態(tài)響應(yīng)時(shí)矩陣中含材料、結(jié)構(gòu)及時(shí)間參數(shù)，由于彈性模量較大、時(shí)間間隔較小，矩陣中出現(xiàn)較0(1)量級(jí)大得多、小得多的數(shù)會(huì)引起新病態(tài)矩陣。本文用式(1)改造力、位移狀態(tài)變量表示的傳遞矩陣，獲得用力、加速度狀態(tài)變量表示的傳遞矩陣為

2.1集總圓盤點(diǎn)傳遞矩陣

由達(dá)郞伯爾原理可導(dǎo)出集中圓盤運(yùn)動(dòng)微分方程為

(2)

式中：ω為自轉(zhuǎn)角速度；CD為集總圓盤處外阻尼系數(shù)；g為重力加速度；eμ為偏心矩。

用式(1)改造式(2)，得用力、加速度狀態(tài)變量表示的集總圓盤點(diǎn)傳遞矩陣為

(3)

式中：D1，D2為集總圓盤點(diǎn)傳遞矩陣，即

式中：Δt為時(shí)間間隔。

2.2無質(zhì)量軸段場(chǎng)傳遞矩陣

由文獻(xiàn)[9]可知無質(zhì)量軸段y方向左右兩端狀態(tài)變量關(guān)系為

(4)

動(dòng)力學(xué)方程為

(5)

式中：K為軸段彈性剛度。

無質(zhì)量軸段z方向左右兩端狀態(tài)變量關(guān)系與動(dòng)力學(xué)方程類似，僅g=零。用式(1)、(5)改造式(4)，獲得用力、加速度狀態(tài)變量表示的無質(zhì)量軸段場(chǎng)傳遞矩陣為

(6)

式中：B1，B2為無質(zhì)量軸段場(chǎng)傳遞矩陣，即

式中：

V=6EJ/aGAL2為考慮軸段的剪切影響，a為截面系數(shù)(實(shí)心軸0.886，空心軸2/3)，G為剪切彈性模量，A為截面積

2.3滑動(dòng)軸承點(diǎn)傳遞矩陣

滑動(dòng)軸承類似有支承的集總圓盤(F≠0)，動(dòng)力學(xué)方程為

(7)

由于式(7)中含非線性項(xiàng)，因此將滑動(dòng)軸承位移、速度及非線性項(xiàng)作為已知量，加速度作為傳遞變量，得

(8)

式中：C1，C2為滑動(dòng)軸承點(diǎn)傳遞矩陣，即

3軸系響應(yīng)Riccati法

對(duì)每一元件(集總圓盤、無質(zhì)量軸段及滑動(dòng)軸承)由式(3)、(6)、(8)獲得相應(yīng)傳遞矩陣，集總圓盤及滑動(dòng)軸承作為節(jié)點(diǎn)處理，與其右邊軸段組成典型單元。該單元傳遞矩陣為

(9)

引入Riccati變換[10]，建立力狀態(tài)變量及加速度狀態(tài)變量之關(guān)系為

(10)

由式(9)、(10)變換得

(11)

(12)

(13)

(P)j+1=(U11P+Ff)j-(S)j+1(U21P+Fe)j

(14)

(Q)j=-(T)j(U21P+F1e)j。

4Riccati-Newmark加速度傳遞矩陣法

對(duì)線性軸系可據(jù)邊界及初始位移、速度條件求出各節(jié)點(diǎn)初始加速度，再通過軸系響應(yīng)求解方法計(jì)算下一時(shí)刻各節(jié)點(diǎn)加速度，由式(1)獲得各節(jié)點(diǎn)位移及速度。而對(duì)非線性軸系含非線性元件節(jié)點(diǎn)的傳遞矩(如滑動(dòng)軸承，式(8))，只在非線性油膜力已知后其節(jié)點(diǎn)傳遞矩陣才真正建立。因此用軸系響應(yīng)Riccati法計(jì)算加速度前須已知非線性項(xiàng)，否則無法計(jì)算。對(duì)此，本文采用Taylor級(jí)數(shù)預(yù)估非線性元件的下一時(shí)刻t+Δt的位移、速度，計(jì)算式為

(15)

本文所提Riccati-Newmark加速度傳遞矩陣法流程見圖2。流程中按式(16)修正t+Δt時(shí)刻位移及速度，即

(16)

式中：ε為迭代收斂需滿足的精度。

采用圖2流程獲得t+Δt時(shí)刻各節(jié)點(diǎn)位移、速度。重復(fù)可得t+2Δt，t+3Δt，….等各時(shí)刻各節(jié)點(diǎn)位移、速度及加速度響應(yīng)。

圖2　Riccati-Newmark加速度傳遞矩陣法 Fig.2 Riccati-Newmark transfer matri acceleration formulation integration method

5數(shù)值例子

采用單圓盤轉(zhuǎn)子系統(tǒng)模型：軸承間隙h0=100 um，軸承寬Lc=28.7mm，轉(zhuǎn)軸跨度L=2 000 mm，直徑d=57.4 mm；集總圓盤直徑D=100 mm, 質(zhì)量為5 kg；圓盤處偏心距eu=1 mm，軸段質(zhì)量密度7.80 g/cm3，彈性模量E=211 GPa，泊松比0.3。油膜動(dòng)力粘度0.016 8 Pa.s。軸段劃分72段，73個(gè)節(jié)點(diǎn)，兩軸承分別作用于第6、第66節(jié)點(diǎn)，圓盤作用于第36節(jié)點(diǎn)；滑動(dòng)軸承用capone圓軸承模型。采用本文提出的方法對(duì)滑動(dòng)軸承-轉(zhuǎn)子系統(tǒng)進(jìn)行非線性瞬態(tài)響應(yīng)分析。

在低速1 200 r/min下用本文方法及Runge-Kutta法(ode45位移法)獲得軸承瞬態(tài)、穩(wěn)態(tài)軌跡見圖3。由圖3看出，兩方法瞬態(tài)解完全吻合，穩(wěn)態(tài)軌跡基本吻合。由此驗(yàn)證本文方法的正解性。Runge-Kutta法為對(duì)下一時(shí)刻軸承處軸頸位移、速度采用線性擾動(dòng)方式，線性擾動(dòng)最后穩(wěn)態(tài)下位移較本文方法偏大，用線性擾動(dòng)處理軸心下一時(shí)刻位移相對(duì)安全。

針對(duì)Taylor級(jí)數(shù)預(yù)估軸心下一時(shí)刻位移的非線擾動(dòng)方法，本文分別對(duì)比分析不同轉(zhuǎn)速下線性、非線性油膜力作用的軸心軌跡，結(jié)果見圖4。由圖4看出，較低轉(zhuǎn)速下，開始線性油膜力Y向(重力作用方向)響應(yīng)較大，穩(wěn)定后響應(yīng)范圍也較大。因此小擾動(dòng)下工程中一般用8個(gè)系數(shù)處理非線性油膜力；但高轉(zhuǎn)速時(shí)非線性油膜力響應(yīng)幅值較線性大。

分析圖3、圖4結(jié)果可知，小擾動(dòng)下采用線性擾動(dòng)響應(yīng)較非線性擾動(dòng)響應(yīng)大。用線性擾動(dòng)模擬下一時(shí)刻軸頸處位移、速度瞬態(tài)分析設(shè)計(jì)偏于安全；但線性擾動(dòng)化處理對(duì)高精密結(jié)構(gòu)設(shè)計(jì)較粗糙。高速下油膜力非線性特性對(duì)軸系瞬態(tài)響應(yīng)有一定影響，采用非線性油膜力及非線性擾動(dòng)的瞬態(tài)響應(yīng)幅值更大。即轉(zhuǎn)速達(dá)一定范圍時(shí)不能視油膜力為線性。

圖3 兩種方法軸承軸心軌跡Fig.3Orbitofshaftcenterbytwomethods圖4 軸承軸心軌跡Fig.4Orbitofshaftcenter

6結(jié)論

(1)本文所提非線性滑動(dòng)軸承-轉(zhuǎn)子系統(tǒng)瞬態(tài)響應(yīng)新方法綜合Riccati傳遞矩陣法與Newmark加速法優(yōu)點(diǎn)，改進(jìn)傳遞矩陣的穩(wěn)定性。通過與傳統(tǒng)瞬態(tài)響應(yīng)分析的軸心下一時(shí)刻位移近似擾動(dòng)所得結(jié)果對(duì)比，驗(yàn)證該方法的準(zhǔn)確性。

(2)通過引入Taylor級(jí)數(shù)預(yù)估，解決非線性元件滑動(dòng)軸承傳遞矩陣建立；通過對(duì)比不同轉(zhuǎn)速非線性擾動(dòng)下線性油膜力、非線性油膜力產(chǎn)生的瞬態(tài)響應(yīng)知，高速下油膜力線性處理不能詮釋軸系非線性現(xiàn)象。

參考文獻(xiàn)

[1]何成兵,顧煜炯.增量傳遞矩陣法及其在軸系彎扭耦合振動(dòng)中的應(yīng)用[J]. 振動(dòng)工程學(xué)報(bào),2006,19(2):119-226.

HE Cheng-bing,GU Yu-jiong.Increment transfer matrix method and its application in the coupled flexural and torsional vibrations analysis[J]. Journal of Vibration Engineering, 2006, 19(2):119-226.

[2]高浩鵬,黃映云,劉鵬.引入傳遞矩陣法的復(fù)雜多體系統(tǒng)連接件建模方法研究[J].振動(dòng)與沖擊,2012,31(6):52-55.

GAO Hao-peng,HUANG Ying-yun,LIU Peng. Modeling of linkers in complex multibody system with transfer matrix method[J].Journal of Vibration and Shock,2012,31(6):52-55.

[3]李蘋,竇海波,王正.水輪發(fā)電機(jī)組主軸系統(tǒng)的建模及非線性瞬態(tài)響應(yīng)[J].清華大學(xué)學(xué)報(bào):自然科學(xué)版,1998,38(6)：123-128.

LI Ping,DOU Hai-bo,WANG Zheng. Modeling for rotor-bearing system of hydroelectric machines and its non-linear transient response[J].Journal of Tsinghua University: Sci. &Tech.,1998,38(6):123-128.

[4]Liu H P, Xu H, Ellison P J,et al.Application of computa-tional fluid dynamica and fluid-structure interaction method to the lubrication study of a rotor-bearing system[J].Tribol Letters,2012,38(3):325-336.

[5]Li qiang,Yu gui-chang,Liu shu-lian. Application of computational fluid dynamics and fluid structure interaction techniques for calculating the 3D transient flow of journal bearings coupled with rotor systems[J].Chinese Journal of Mechanical Engineering,2013,25(5):926-931.

[6]顧致平,劉永壽.非線性轉(zhuǎn)子系統(tǒng)中的傳遞矩陣技術(shù)[M].北京:科學(xué)出版社, 2010.

[7]王正. Riccati傳遞矩陣法的奇點(diǎn)及其消除方法[J]. 振動(dòng)與沖擊, 1987,6(2): 74-78.

WANG Zheng. Singular point and elimination method of riccati transfer matrix method[J]. Journal of Vibration and Shock, 1987,6 (2):74-78.

[8]徐斌. Matlab有限元結(jié)構(gòu)動(dòng)力學(xué)分析與工程應(yīng)用[M].北京:清華大學(xué)出版社, 2010.

[9]聞邦椿,顧家柳,夏松波,等. 高等轉(zhuǎn)子動(dòng)力學(xué)[M]. 北京：機(jī)械工業(yè)出版社, 2000.

[10]毛文貴,劉桂萍,韓旭,等.一種高效的電主軸系統(tǒng)復(fù)頻率計(jì)算方法[J].振動(dòng)與沖擊,2014,33(19):164-168.

MAO Wen-gui,LIU Gui-ping,HAN Xu,et al.An efficent method for calculating complex frequency of a motorized spindle[J].Journal of Vibration and Shock, 2014, 33(19): 164-168.