崔 立 鄭建榮

1.上海第二工業(yè)大學(xué)，上海，201209 2.華東理工大學(xué)，上海，200237

0 引言

滾動軸承支承的轉(zhuǎn)子系統(tǒng)由于軸承非線性剛度而出現(xiàn)復(fù)雜的非線性動力學(xué)特性。近年來，剛性轉(zhuǎn)子滾動軸承系統(tǒng)的分叉及混沌特性研究得到了廣泛的重視。

Jang等［1］考慮滾動軸承的5個自由度以及表面波紋度建立了剛性轉(zhuǎn)子滾動軸承系統(tǒng)的非線性動力學(xué)方程，發(fā)現(xiàn)套圈的波紋會導(dǎo)致徑向位移和角位移響應(yīng)峰值附近出現(xiàn)邊頻帶。Harsha［2－3］考慮軸承間隙、不平衡力研究了轉(zhuǎn)子系統(tǒng)的動力學(xué)響應(yīng)，發(fā)現(xiàn)系統(tǒng)的動力學(xué)行為與滾動體的通過頻率有關(guān)，且當(dāng)通過頻率及其諧波與固有頻率相等、外圈波紋度階數(shù)與滾動體數(shù)目相等時振動幅值將變得極大。Bai等［4］研究了主軸－滾動軸承系統(tǒng)非線性動力學(xué)行為，發(fā)現(xiàn)軸承間隙的減小有利于提高主軸軸承系統(tǒng)的穩(wěn)定性。高尚晗等［5］研究了主軸－滾動軸承在負(fù)游隙情況下的機(jī)床主軸－滾動軸承系統(tǒng)的非線性動力學(xué)特性，揭示了主軸系統(tǒng)的混沌演化過程。崔立等［6］研究了圓柱滾子軸承剛性轉(zhuǎn)子系統(tǒng)周期運(yùn)動分岔特性，發(fā)現(xiàn)隨著徑向間隙、阻尼和力矩的變化，周期運(yùn)動將產(chǎn)生倍周期或Hopf分岔，分岔轉(zhuǎn)速隨參數(shù)變化而改變。

以上的研究對象均為剛性轉(zhuǎn)子，模型涉及轉(zhuǎn)軸的彎曲變形和陀螺力矩等參數(shù)。隨著旋轉(zhuǎn)機(jī)械轉(zhuǎn)速的提高，柔性轉(zhuǎn)子系統(tǒng)的設(shè)計與分析變得越來越重要。近年也有一些研究柔性轉(zhuǎn)子軸承系統(tǒng)的文章發(fā)表，如：Laha等［7］研究了油膜軸承支承的柔性轉(zhuǎn)子系統(tǒng)的分叉行為，分析了轉(zhuǎn)軸的材料、剛度和質(zhì)量等參數(shù)對轉(zhuǎn)子系統(tǒng)分叉行為的影響。Villa等［8］、Sinou［9］研究了球軸承支承的柔性轉(zhuǎn)子系統(tǒng)的非線性動力學(xué)行為，發(fā)現(xiàn)系統(tǒng)響應(yīng)中存在跳躍現(xiàn)象和超諧波，但其研究僅考慮了4個自由度的轉(zhuǎn)軸節(jié)點(diǎn)和2個軸承自由度，難以滿足實際工況的需求。

本文采用12自由度Euler－Bernoulli桿單元建立柔性轉(zhuǎn)子滾子軸承系統(tǒng)的非線性動力學(xué)模型，研究系統(tǒng)的混沌行為，分析軸承結(jié)構(gòu)參數(shù)與轉(zhuǎn)軸結(jié)構(gòu)參數(shù)對轉(zhuǎn)子系統(tǒng)混沌行為的影響規(guī)律。

1 計算模型

圖1所示為柔性轉(zhuǎn)子軸承系統(tǒng)模型，采用12自由度的Euler－Bernoulli桿單元，圓盤簡化為質(zhì)點(diǎn)并考慮其質(zhì)量與轉(zhuǎn)動慣量，轉(zhuǎn)子由兩個滾動軸承支承。

圖1 柔性轉(zhuǎn)子滾動軸承系統(tǒng)模型

1.1 柔性轉(zhuǎn)子軸承系統(tǒng)動力學(xué)模型

柔性轉(zhuǎn)子系統(tǒng)的動力學(xué)方程為

式中，M為包含轉(zhuǎn)軸、圓盤和軸承的總質(zhì)量矩陣；C為總阻尼矩陣；G為總陀螺矩陣；K為總剛度矩陣；X為轉(zhuǎn)子各節(jié)點(diǎn)的位移向量；f（X，t）為包括轉(zhuǎn)子各節(jié)點(diǎn)軸承力、重力、轉(zhuǎn)子不平衡力和外載荷的矩陣。

先對轉(zhuǎn)子系統(tǒng)進(jìn)行軸段與節(jié)點(diǎn)劃分，然后進(jìn)行轉(zhuǎn)子系統(tǒng)的質(zhì)量矩陣、剛度矩陣、阻尼矩陣、陀螺矩陣、載荷矩陣求解，并按照節(jié)點(diǎn)的順序?qū)Ω骶仃囘M(jìn)行組裝，具體過程如下：首先計算桿單元的質(zhì)量矩陣、陀螺矩陣、剛度矩陣、載荷矩陣并組裝；然后將剛性圓盤的質(zhì)量矩陣、陀螺矩陣、不平衡力矩陣疊加到所在節(jié)點(diǎn)的相應(yīng)矩陣中；之后將支承軸承的質(zhì)量矩陣、非線性軸承力疊加到所在節(jié)點(diǎn)的相應(yīng)矩陣中；最后計算系統(tǒng)的結(jié)構(gòu)阻尼、軸承阻尼矩陣并組裝。

1.2 桿單元

圖2所示為針對空間桿單元分析獲得的空間桿單元兩端節(jié)點(diǎn)位移。

圖2 二節(jié)點(diǎn)空間桿單元

針對圖2中12自由度Euler－Bernoulli桿單元，若每個節(jié)點(diǎn)考慮6個自由度，則共有6個廣義位移和6個廣義力，其表達(dá)式為

式中，xi、yi、zi分別為節(jié)點(diǎn)i沿x、y、z方向的線位移；θxi、θyi、θzi分別為節(jié)點(diǎn)i處截面繞x、y、z軸的轉(zhuǎn)角位移；Fxi、Fyi、Fzi分別為節(jié)點(diǎn)i在x、y、z方向受到的軸向力和剪切力；Mxi、Myi、Mzi分別為節(jié)點(diǎn)i在x、y、z方向所受的扭矩和彎矩。

假設(shè)已知桿單元橫截面面積、截面慣性矩、單元的扭轉(zhuǎn)慣性矩、長度、材料彈性模量和剪切模量，則可根據(jù)有限元理論［10］求出桿單元的質(zhì)量矩陣Ms、剛度矩陣Ks、阻尼矩陣Cs、陀螺矩陣Gs。

1.3 剛性圓盤模型

當(dāng)軸上安裝有圓盤時，將其視為剛性圓盤，并將其質(zhì)量矩陣、陀螺矩陣、不平衡力矩陣疊加到所在節(jié)點(diǎn)的相應(yīng)矩陣中。

假設(shè)已知圓盤的質(zhì)量、半徑、不平衡質(zhì)徑積，則可建立剛性圓盤的質(zhì)量矩陣Md、陀螺矩陣Gd、不平衡力矩陣Fd，其中Fd為

式中，me為圓盤的不平衡質(zhì)徑積；ω為轉(zhuǎn)速。

1.4 滾子軸承模型

考慮普遍受載的圓柱滾子軸承，其模型如圖3所示，圖中，Dr為滾子直徑，D1、D2分別為外圈外徑和內(nèi)圈內(nèi)徑，δ0為徑向間隙，le為帶凸度滾子的長度，ls為帶凸度滾子直線部分的長度。

假設(shè)滾子數(shù)目為N，使用切片法將滾子分成nr個圓片。對第j個滾子進(jìn)行受力分析，假設(shè)第j個滾子方位角為φj，根據(jù)赫茲接觸理論，并使用擬動力學(xué)方法建立滾子軸承的非線性平衡方程組，使用Newton－Raphson法可求出滾子與套圈的接觸力［11］。將各滾子的接觸力分解，即可求出圓柱滾子軸承的非線性軸承力矩陣：

圖3 圓柱滾子軸承結(jié)構(gòu)簡圖

式中，F(xiàn)2jk為第j個滾子的第k個切片與滾子軸承內(nèi)圈的作用力。

1.5 結(jié)構(gòu)阻尼求解

Rayleigh提出的結(jié)構(gòu)阻尼模型計算表達(dá)式為

式中，ε1、ε2為兩個振型的阻尼比，根據(jù)經(jīng)驗取ε1＝0.005，ε2＝0.01；ω1、ω2為轉(zhuǎn)子系統(tǒng)的二階固有頻率。

式（7）中轉(zhuǎn)子系統(tǒng)的固有頻率ω1、ω2可根據(jù)計算得到的質(zhì)量矩陣、剛度矩陣解｜K－Mω2｜＝0得到。

2 計算方法

采用Runge－Kutta法、Newton－Raphson法進(jìn)行非線性動力學(xué)方程組求解，根據(jù)FPA修正法確定求解周期，然后求解最大Lyapunov指數(shù)，判斷系統(tǒng)的動力學(xué)行為。

2.1 求解周期的確定

轉(zhuǎn)子軸承系統(tǒng)中存在軸承變剛度激勵，還可能存在不平衡力激勵。不平衡力產(chǎn)生的激勵周期為軸轉(zhuǎn)動周期的整數(shù)倍，但軸承的變剛度激勵周期往往不是軸轉(zhuǎn)動周期的整數(shù)倍，所以在判斷和求解時，采用修正的FPA法建立統(tǒng)一的求解周期［12］，其表達(dá)式定義為

式中，T為求解周期；Td為不平衡力激勵周期；TVC為軸承變剛度激勵周期；ε為常數(shù)，取ε＝0.01；K為比例系數(shù)，K ＝1，2，…，nk。

根據(jù)式（9）進(jìn)行循環(huán)計算，直至找到滿足其要求的K值，代入式（8）即得求解周期。

2.2 Lyapunov指數(shù)計算

Lyapunov指數(shù)表示在相平面中2條相鄰軌線間的距離隨時間的平均指數(shù)發(fā)散率，它明確地區(qū)分了確定性運(yùn)動和混沌運(yùn)動［13］。

對于連續(xù)系統(tǒng)有

設(shè)在τ0時刻‖δX（τ0）‖充分小，于是1維的Lyapunov指數(shù)可定義為

在n維連續(xù)系統(tǒng)中，δX（τ）在每個基底上有分量，每一個分量均可按上式求出一個λ，因此共存在n個Lyapunov指數(shù)λi，稱為Lyapunov指數(shù)譜。當(dāng)任意選取矩陣δX（τ）時，Lyapunov指數(shù)以概率1可能取得最大值，如果其中最大的Lyapunov指數(shù)λmax＞0，則該系統(tǒng)一定存在混沌運(yùn)動。因此，只要計算出系統(tǒng)的最大Lyapunov指數(shù)，就可以判斷系統(tǒng)是否處于混沌狀態(tài)。

3 計算與分析

圖4為柔性轉(zhuǎn)子系統(tǒng)簡圖，轉(zhuǎn)軸由兩個型號相同的滾子軸承支承，轉(zhuǎn)軸中有剛性圓盤，該圓盤可施加不平衡力。

圖4 柔性轉(zhuǎn)子系統(tǒng)簡圖

表1所示為滾子軸承的結(jié)構(gòu)參數(shù)，其中，軸承的彈性模量為204GPa，泊松比為0.3，阻尼為200N·s／m，軸承載荷為｛0，2000N，2000N，0，0｝。

表1 滾子軸承的結(jié)構(gòu)參數(shù)

柔性轉(zhuǎn)子系統(tǒng)的結(jié)構(gòu)參數(shù)如表2所示，轉(zhuǎn)軸被劃分為4個軸段，轉(zhuǎn)子系統(tǒng)劃分為5個節(jié)點(diǎn)。轉(zhuǎn)軸的彈性模量為204GPa，泊松比為0.3。

表2 柔性轉(zhuǎn)子系統(tǒng)結(jié)構(gòu)參數(shù)

對圖4所示的柔性轉(zhuǎn)子系統(tǒng)進(jìn)行計算，判斷圓盤節(jié)點(diǎn)處的混沌行為，并分析軸承徑向間隙、圓盤不平衡力、轉(zhuǎn)軸剛度比等參數(shù)對系統(tǒng)混沌行為的影響規(guī)律。

3.1 徑向間隙影響分析

假設(shè)系統(tǒng)不平衡力為0，則軸承徑向間隙對系統(tǒng)混沌特性的影響如圖5所示。分別取40μm、60μm和80μm的徑向間隙計算系統(tǒng)的最大Lyapunov指數(shù)。

圖5 不同節(jié)點(diǎn)處的最大Lyapunov指數(shù)

系統(tǒng)最大Lyapunov指數(shù)小于0，表明系統(tǒng)運(yùn)動是穩(wěn)定的；當(dāng)最大Lyapunov指數(shù)等于0時，系統(tǒng)為倍周期或擬周期分叉運(yùn)動；當(dāng)最大Lyapunov指數(shù)大于0時，系統(tǒng)為混沌運(yùn)動。

為了驗證計算的準(zhǔn)確性，使用Runge－Kutta法進(jìn)行非線性動力學(xué)方程組求解，得到轉(zhuǎn)子系統(tǒng)響應(yīng)的頻譜圖和Poincaré截面。從圖5可知，當(dāng)徑向間隙為80μm、轉(zhuǎn)速為4200r／min時，系統(tǒng)的最大Lyapunov指數(shù)等于0（圖中fVC為滾子軸承的變剛度振動頻率）。圖6所示為圖5工況下圓盤節(jié)點(diǎn)處的頻譜圖和Poincaré截面，圖形表明系統(tǒng)為擬周期振動。

圖6 轉(zhuǎn)速為4200r／min、徑向間隙為80μm時的響應(yīng)

圖5表明，徑向間隙為80μm、轉(zhuǎn)速為8000 r／min時，系統(tǒng)的最大Lyapunov指數(shù)大于0，圖7所示為對應(yīng)該工況的圓盤節(jié)點(diǎn)處的頻譜圖和Poincaré截面，其表明系統(tǒng)為混沌運(yùn)動。

圖5還表明，當(dāng)徑向間隙為40μm時，系統(tǒng)的最大Lyapunov指數(shù)小于0，未出現(xiàn)混沌行為；當(dāng)徑向間隙為60μm 時，在0～2500r／min、6600～7200r／min轉(zhuǎn)速下的最大Lyapunov指數(shù)大于0，系統(tǒng)出現(xiàn)混沌行為；當(dāng)徑向間隙為80μm時，在0～3100r／min、5000～5400r／min、7500～8500 r／min轉(zhuǎn)速下出現(xiàn)混沌行為。

3.2 不平衡力影響分析

圖8所示為在其他條件不變前提下的不平衡力對系統(tǒng)混沌特性的影響。分別取圓盤的不平衡質(zhì)徑積為5×10－6kg·m、1×10－5kg·m、2×10－5kg·m，計算系統(tǒng)的最大Lyapunov指數(shù)。

從圖8可知，當(dāng)不平衡質(zhì)徑積為5×10－6kg·m時，在0～2000r／min，系統(tǒng)出現(xiàn)混沌行為；當(dāng)不平衡質(zhì)徑積為1×10－5kg·m時，在0～3700r／min，系統(tǒng)出現(xiàn)混沌行為；當(dāng)不平衡質(zhì)徑積為2×10－5kg·m時，在0～4400r／min，系統(tǒng)出現(xiàn)混沌行為?？梢钥闯?，不平衡力的存在會改變系統(tǒng)的混沌行為，系統(tǒng)出現(xiàn)混沌的轉(zhuǎn)速及范圍隨著不平衡力的增大逐漸增大。

圖7 轉(zhuǎn)速為8000r／min、徑向間隙為80μm時的響應(yīng)

圖8 不同不平衡力時的系統(tǒng)最大Lyapunov指數(shù)

3.3 剛度比影響分析

定義柔性轉(zhuǎn)子系統(tǒng)的剛度比為

式中，Kshaft為柔性轉(zhuǎn)軸的剛度；Kbearing為支承軸承剛度。

在其他參數(shù)不變的情況下，分別取長度為500mm、400mm、300mm的轉(zhuǎn)軸進(jìn)行剛度比計算，其對應(yīng)所得的剛度比分別為7.38×10－3、1.43×10－2、3.31×10－2，據(jù)此，可得出剛度對于柔性轉(zhuǎn)子系統(tǒng)動力學(xué)特性的影響。

圖9所示為不同剛度比時圓盤節(jié)點(diǎn)處的振幅隨轉(zhuǎn)速變化曲線，可以看出，隨著剛度比的增大，振幅峰值對應(yīng)的轉(zhuǎn)速逐漸增大，即臨界轉(zhuǎn)速增大；振幅峰值也隨著剛度比增大而增大，且臨界轉(zhuǎn)速附近的峰值逐漸變多，可以看出隨著轉(zhuǎn)軸剛度增大，軸承非線性振動對轉(zhuǎn)子系統(tǒng)的影響增大，導(dǎo)致系統(tǒng)響應(yīng)復(fù)雜。

圖9 不同剛度比時的圓盤節(jié)點(diǎn)振幅

圖10所示為不同剛度比時最大Lyapunov指數(shù)隨轉(zhuǎn)速的變化曲線。當(dāng)剛度比為7.38×10－3時，在0～800r／min出現(xiàn)混沌行為。當(dāng)剛度比為1.43×10－2時，在0～1200r／min、5200～5400r／min出現(xiàn)混沌行為。當(dāng)剛度比為3.31×10－2時，在0～3200r／min、4800～5900r／min出現(xiàn)混沌行為。

圖10 不同轉(zhuǎn)軸剛度比時的系統(tǒng)最大Lyapunov指數(shù)

可見，隨著轉(zhuǎn)軸剛度增大即剛度比增大，混沌運(yùn)動的區(qū)間發(fā)生改變。隨著剛度比的增大，系統(tǒng)的混沌區(qū)間增大，軸承引起的轉(zhuǎn)子系統(tǒng)非線性行為明顯。對比柔性轉(zhuǎn)子和剛性轉(zhuǎn)子的混沌特性，發(fā)現(xiàn)軸承的非線性接觸力對剛性轉(zhuǎn)子系統(tǒng)的混沌特性影響大于柔性轉(zhuǎn)子。

4 實驗

為驗證本文方法對柔性轉(zhuǎn)子系統(tǒng)動力學(xué)行為預(yù)測的準(zhǔn)確性，使用圖11所示的高速滾子軸承柔性轉(zhuǎn)子實驗器進(jìn)行測試。動力裝置的高速電主軸，最高轉(zhuǎn)速可達(dá)24 000r／min，實驗滾子軸承參數(shù)如表1所示，實驗過程中滾子軸承承受的徑向載荷為2000N。采用非接觸式的電渦流傳感器測量滾子軸承及圓盤處的徑向振動位移。

圖12所示為轉(zhuǎn)子圓盤處振動位移幅值計算結(jié)果與實驗結(jié)果對比，可以看出，計算得到的一階臨界轉(zhuǎn)速為4900r／min，實驗得到的一階臨界轉(zhuǎn)速為4700r／min，與計算結(jié)果較為接近；實驗測試振幅與計算結(jié)果也較為接近，證明了本文方法的準(zhǔn)確性。

圖11 滾動軸承轉(zhuǎn)子系統(tǒng)實驗器

圖12 幅值計算結(jié)果與實驗結(jié)果對比

5 結(jié)論

（1）軸承徑向間隙是影響轉(zhuǎn)子系統(tǒng)非線性振動特性的重要參數(shù)。隨著軸承徑向間隙的增大，系統(tǒng)的混沌區(qū)間逐漸增大、變多。

（2）不平衡力的存在對系統(tǒng)混沌行為也有較大的影響。系統(tǒng)出現(xiàn)混沌的轉(zhuǎn)速及范圍隨著不平衡力的增大逐漸增大。

（3）隨著轉(zhuǎn)軸剛度比增大，轉(zhuǎn)子振幅峰值以及對應(yīng)的轉(zhuǎn)速均逐漸增大。共振轉(zhuǎn)速附近產(chǎn)生混沌運(yùn)動對應(yīng)的轉(zhuǎn)速增大時，系統(tǒng)混沌區(qū)間也增多，軸承非線性對系統(tǒng)響應(yīng)影響明顯。

［1］Jang G H，Jeong S W.Nonlinear Excitation Model of Ball Bearing Waviness in a Rigid Rotor Supported by Two or More Ball Bearings Considering Five Degrees of Freedom［J］.Transactions of the ASME，2002，124：82－90.

［2］Harsha S P.Non－linear Dynamic Analysis of an Unbalanced Rotor Supported by Roller Bearing［J］.Chaos，Solitons and Fractals，2005，26：47－66.

［3］Harsha S P.Non－linear Dynamic Response of a Balanced Rotor Supported by Rolling Element Bearings Due to Radial Internal Clearance Effect［J］.Mechanism and Machine Theory，2006，41：688－706.

［4］Bai Changqing，Xu Qingyu.Dynamic Model of Ball Bearings with Internal Clearance and Waviness［J］.Journal of Sound and Vibration，2006，294：23－48.

［5］高尚晗，龍新華，孟光.主軸－滾動軸承系統(tǒng)三種分岔形式［J］.振動與沖擊，2009，28（4）：59－64.Gao Shanghan，Long Xinhua，Meng Guang.Three Bifurcation Types in the Spindle－ball Bearing System［J］.Journal of Vibration and Shock，2009，28（4）：59－64.

［6］崔立，劉長利，鄭建榮.圓柱滾子軸承剛性轉(zhuǎn)子系統(tǒng)周期運(yùn)動分岔研究［J］.振動、測試與診斷，2011，31（5）：637－641.Cui Li，Liu Changli，Zheng Jianrong.Periodic Motion Bifurcation of Rigid Rotor System in Cylinder Roller Bearing［J］.Journal of Vibration，Measurement ＆ Diagnosis，2011，31（5）：637－641.

［7］Laha S K，Kakoty S K.Non－linear Dynamic Analysis of a Flexible Rotor Supported on Porous Oil Journal Bearings［J］.Communications in Nonlinear Science and Numerical Simulation，2011，16：1617－1631.

［8］Villa C，Sinou J J，Thouverez F.Stability and Vibration Analysis of a Complex Flexible Rotor Bearing System［J］.Communications in Non－linear Science and Numerical Simulation，2008，13：804－821.

［9］Sinou J J.Non－linear Dynamics and Contacts of an Unbalanced Flexible Rotor Supported on Ball Bearings［J］.Mechanism and Machine Theory，2009，44：1713－1732.

［10］Wang Maocheng.Finite Element Method［M］.Beijing：Tsinghua University Press，2008.

［11］Harris T A.Rolling Bearing Analysis［M］.New York，USA：John Wiley ＆Sons，2001.

［12］Choi S K，Noah S T.Response and Stability Analysis of Piecewise Linear Oscillations under Multi－forcing Frequencies［J］.Nonlinear Dynamics，1992，3：105－121

［13］Wolf A，Swift J B，Swinney H L.Determining Lyapunov Exponents from a Time Series［J］.Physica D.，1985，16：285－317.