接觸界面對(duì)拉桿組合柔性轉(zhuǎn)子軸承系統(tǒng)的非線性動(dòng)力特性影響

易 均，劉 恒，劉 意，于 明，王為民

(1.西安交通大學(xué) 機(jī)械工程學(xué)院，西安 710049;2.東方汽輪機(jī)有限公司，四川 德陽(yáng) 618000)

現(xiàn)代燃?xì)廨啓C(jī)、航空發(fā)動(dòng)機(jī)等大型動(dòng)力系統(tǒng)的核心部件均為典型的拉桿柔性組合轉(zhuǎn)子。此類轉(zhuǎn)子的多個(gè)獨(dú)立輪盤(pán)間依靠接觸界面?zhèn)鬟f各向作用力，并在拉桿預(yù)緊下組合為一體。重量輕、易于裝配且具有良好的冷卻效果是此類轉(zhuǎn)子的特點(diǎn)，針對(duì)高速化、高效化的工業(yè)發(fā)展，應(yīng)用前景非常廣泛。

目前，針對(duì)一般整體轉(zhuǎn)子穩(wěn)定性的討論已有很多，如剛性軸的穩(wěn)定性有文獻(xiàn)［1］等，簡(jiǎn)單Jeffcott柔性轉(zhuǎn)子軸承系統(tǒng)穩(wěn)定性有文獻(xiàn)［2］等，一般柔性轉(zhuǎn)子軸承系統(tǒng)則有文獻(xiàn)［3］等。然而拉桿組合轉(zhuǎn)子輪盤(pán)間非線性本質(zhì)的接觸效應(yīng)，導(dǎo)致此類組合轉(zhuǎn)子其結(jié)構(gòu)本質(zhì)上是不連續(xù)的，傳統(tǒng)的視作整體轉(zhuǎn)軸的轉(zhuǎn)子處理方式不適用于此類分布式特種轉(zhuǎn)子。針對(duì)拉桿組合轉(zhuǎn)子動(dòng)力特性的研究，饒柱石等［4-5］提出了拉桿組合轉(zhuǎn)子的一般力學(xué)模型，王艾倫等［6-8］則進(jìn)行了計(jì)及接觸界面的臨界轉(zhuǎn)速計(jì)算，但是此處非線性動(dòng)力特性的研究仍然十分有限。

運(yùn)用文獻(xiàn)［9］提出的粗糙機(jī)械結(jié)合面接觸剛度的研究方法，通過(guò)計(jì)算與真實(shí)接觸面等條件下的微元模型界面接觸剛度并將之進(jìn)行面積擴(kuò)展，得到輪盤(pán)界面接觸剛度;將接觸界面處理成無(wú)質(zhì)量均質(zhì)彈簧［5］，采用哈密頓原理完成了計(jì)及接觸界面的轉(zhuǎn)子建模，得到了拉桿轉(zhuǎn)子輪盤(pán)間界面產(chǎn)生的附加剛度矩陣;而后結(jié)合對(duì)整體轉(zhuǎn)軸應(yīng)用計(jì)及軸向力的鐵木辛格梁軸單元建立的有限元模型，得出了轉(zhuǎn)子軸承系統(tǒng)的動(dòng)力學(xué)方程。最后運(yùn)用打靶法結(jié)合Floquet穩(wěn)定性分叉理論，對(duì)比了計(jì)及接觸界面前后的系統(tǒng)穩(wěn)態(tài)周期解的穩(wěn)定性邊界和分叉形式，數(shù)值結(jié)果表明接觸界面對(duì)不平衡轉(zhuǎn)子動(dòng)力特性的影響不可忽略。

1 接觸界面建模

拉桿組合轉(zhuǎn)子的輪盤(pán)是通過(guò)拉桿預(yù)緊而組合在一起的。其輪盤(pán)間接觸界面在預(yù)緊力作用下發(fā)生變形產(chǎn)生作用力而傳遞軸向力，因此可以將接觸界面單獨(dú)視作一個(gè)無(wú)質(zhì)量的面彈簧單元，在對(duì)拉桿組合轉(zhuǎn)子建模時(shí)，將接觸界面作為一個(gè)附加彈簧單元的形式添加到系統(tǒng)，如圖1所示。

圖1 接觸界面彈簧替代示意圖Fig.1 Contact interface replaced by springs

附加彈簧單元?jiǎng)偠却笮〖礊榻缑娼佑|剛度，由預(yù)緊力大小和接觸界面屬性(材料，粗糙度，形貌度，波浪度等)決定。在彈性范圍內(nèi)，預(yù)緊力越大則其剛度越大;接觸界面粗糙度、形貌度和波浪度等都是影響界面接觸剛度大小的幾何因素，其中粗糙度為研究各種幾何因素的基礎(chǔ)。此處研究采用的模型其接觸界面僅僅計(jì)入了粗糙度的影響，對(duì)形貌度等其他幾何因素影響的研究，需綜合粗糙度的影響，并在此基礎(chǔ)上進(jìn)一步深入討論。實(shí)際情況中接觸界面形貌度和波浪度對(duì)界面剛度的影響均會(huì)大于單一的粗糙度的影響，二者導(dǎo)致的接觸界面剛度變化對(duì)動(dòng)力特性的影響也會(huì)比單一的粗糙度因素更大。

1.1 彈簧單元?jiǎng)偠却_定

直接在米級(jí)的全尺寸模型中研究微米級(jí)的粗糙度對(duì)界面剛度的影響，巨大的計(jì)算尺度跨越使得此處分析非常困難。本文采用等條件的微元有限元模型分析與宏觀尺寸相結(jié)合的跨尺度計(jì)算方法，從而準(zhǔn)確得到輪盤(pán)間接觸界面剛度，其步驟為:① 計(jì)算輪盤(pán)間如圖1所示A接觸界面面積s2，對(duì)輪盤(pán)接觸界面進(jìn)行受力分析，得到預(yù)緊后的真實(shí)界面壓力P;② 根據(jù)文獻(xiàn)［9］方法建立在壓力P作用下的面積為s1、與輪盤(pán)接觸面等條件的微觀有限元模型如圖2所示，對(duì)此模型進(jìn)行應(yīng)力應(yīng)變分析并提取該微元法向接觸剛度k;③ 由于線性面彈簧與彈簧面積成正比，輪盤(pán)間法向接觸剛度kf由式(1)得到，切向剛度kq根據(jù)式(2)［10］得到。

圖2 微觀尺度的有限元計(jì)算模型Fig.2 Micro-scale finite element model

式中:A=2，υ=0.3 為泊松比，取kq=0.82kf。

1.2 附加彈簧單元建立

輪盤(pán)間接觸界面處理成無(wú)質(zhì)量、法向剛度為kf、切向剛度為kq的附加彈簧單元，其在系統(tǒng)中的附加剛度矩陣可根據(jù)彈性勢(shì)能定律推導(dǎo)得到。附加彈簧單元?jiǎng)菽艽笮∮山佑|界面兩側(cè)附屬節(jié)點(diǎn)的動(dòng)態(tài)坐標(biāo)確定。對(duì)于接觸端面中心點(diǎn)來(lái)說(shuō)，其坐標(biāo)由固結(jié)連接于盤(pán)軸中心節(jié)點(diǎn)的動(dòng)態(tài)坐標(biāo)(x，y，ψ，φ)確定(輪盤(pán)與其鄰近軸段簡(jiǎn)稱盤(pán)軸單元)，因而其上距盤(pán)軸中心r的一點(diǎn)的坐標(biāo)可如下確定為:建立系統(tǒng)的絕對(duì)坐標(biāo)系oxyz，在此坐標(biāo)系中t時(shí)刻盤(pán)軸1中心的坐標(biāo)為(x1，y1，z1)，角坐標(biāo)為(ψ1，φ1，ωt)。然后通過(guò)移軸將坐標(biāo)系原點(diǎn)移至所考慮盤(pán)軸的中心，得到盤(pán)軸坐標(biāo)系o'x'y'z'，如圖3所示，再將圓盤(pán)繞y'軸旋轉(zhuǎn)φ1角到達(dá)o1x1y1z1，再繞x1軸反向旋轉(zhuǎn)ψ1角到達(dá)o2x2y2z2，最后再繞z2軸以角速度ω正向旋轉(zhuǎn)ωt角，最終到達(dá)圓盤(pán)的隨動(dòng)坐標(biāo)系o3x3y3z3。

圖3 剛性盤(pán)軸段坐標(biāo)變換示意圖Fig.3 Coordinate transformation of rigid plate

因此，在盤(pán)軸的隨動(dòng)坐標(biāo)系o3x3y3z3中，接觸界面的固接點(diǎn)A(xA，yA，zA)的坐標(biāo)可表示為式(3)，其中x1，y1，z1分別是盤(pán)1 的動(dòng)態(tài)位移，ψ1，φ1則為盤(pán)1 的動(dòng)態(tài)轉(zhuǎn)角，(為附加彈簧單元的端點(diǎn)在圓盤(pán)上的位置轉(zhuǎn)角。接觸界面另一端對(duì)應(yīng)的固結(jié)點(diǎn)B(xB，yB，zB)的坐標(biāo)可類似表出。

式中:a=rcos(ωt+γ)。

根據(jù)線性彈簧變形的勢(shì)能原理，接觸界面的彈性勢(shì)能表示為:

在徑向x方向上的勢(shì)能Ux，可根據(jù)彈簧此方向上變化量產(chǎn)生的勢(shì)能對(duì)接觸界面積分得出:

Δx2取一階近似后可得:

式(7)代入式(5)可知:

在徑向y方向上的勢(shì)能Uy，其彈簧勢(shì)能與x方向推導(dǎo)類似，有:

在軸向z方向上，有:

可知:

其中:

可推出接觸界面處理為附加彈簧單元后，其附加剛度矩陣ke為:

2 系統(tǒng)方程

對(duì)于預(yù)緊后的組合轉(zhuǎn)子整體轉(zhuǎn)軸，直接采用了計(jì)及軸向力的鐵木辛格梁軸單元進(jìn)行離散［11］，如圖4所示。在整體轉(zhuǎn)軸的動(dòng)力方程中，根據(jù)哈密頓原理加入輪盤(pán)間接觸界面、剛性圓盤(pán)和周向均布拉桿［12］的影響后，最終給出預(yù)緊后的拉桿組合轉(zhuǎn)子軸承系統(tǒng)動(dòng)力學(xué)方程:

式中:MS=Mr+Md，GS=Gr+Gd，

其中:Ms，Gs，Ks∈Rm×m分別是系統(tǒng)的質(zhì)量、阻尼和剛度矩陣;Qs，F(xiàn)s∈Rm×m分別是系統(tǒng)的重力和外力矢量;Ke為附加彈簧單元提供的附加剛度矩陣;Mr，Gr，Qr∈Rm×m分別為轉(zhuǎn)軸的質(zhì)量、阻尼、剛度矩陣;Md，Gd，Qd，fexd分別為剛性圓盤(pán)的質(zhì)量、陀螺矩陣及重力和不平衡力矢量;foil為軸承力;系統(tǒng)的位移矢量為:xS={x1y1ψ1φ1…xmymψmφm}T。

其中:xj，yj和ψj，φj(j=1，2，…，p)分別為第j個(gè)節(jié)點(diǎn)的橫向位移及轉(zhuǎn)動(dòng)自由度。

圖4 預(yù)緊后的鐵木辛格梁軸有限單元模型Fig.4 Timoshenko beam-shaft element after tightening

3 分析方法

假設(shè)作用于系統(tǒng)的載荷為周期為T(mén)的周期載荷:

則系統(tǒng)的穩(wěn)態(tài)解為周期解、偽周期解及混沌解。對(duì)于軸承轉(zhuǎn)子系統(tǒng)而言，工頻周期解是其本征解，隨外參數(shù)轉(zhuǎn)速ω、不平衡量e等的變化，此解將發(fā)生失穩(wěn)，而產(chǎn)生新的穩(wěn)態(tài)解形式，如偽周期解及混沌解，或發(fā)生絕對(duì)失穩(wěn)，因此求取系統(tǒng)的周期解及其穩(wěn)定性分叉規(guī)律便成為研究此類系統(tǒng)非線性特性的主要內(nèi)容。

3.1 周期解求解

應(yīng)用打靶法求解式(13)動(dòng)力系統(tǒng)周期解的問(wèn)題可轉(zhuǎn)化為兩點(diǎn)邊值問(wèn)題:

其中:

式中:μ為系統(tǒng)參數(shù)，如系統(tǒng)轉(zhuǎn)速ω、不平衡量e等等。將式(15)積分一個(gè)周期T，則如下式得到滿足時(shí)可求得一個(gè)T周期的周期解Xs。

對(duì)于一個(gè)給定的參數(shù)μ=μs，對(duì)應(yīng)的周期解Xs可通過(guò)對(duì)式(17)應(yīng)用牛頓迭代法求得。其雅克比矩陣為:

矩陣J可通過(guò)將下式與式(15)一起對(duì)于系統(tǒng)軌跡Xs(t0+t)積分一個(gè)周期求出:)

顯然當(dāng)以［Xs(t0)，I］作為初始值積分一個(gè)Poincare 映射周期T，δS(t0)=I時(shí)，δS(t0+t)=J。以上即為打靶法的基本思路。在此基礎(chǔ)上，采用周期解延續(xù)追蹤算法［13］當(dāng)外參數(shù)變化時(shí)對(duì)周期解進(jìn)行預(yù)估校正，有效地得到外參數(shù)變化時(shí)周期解的演化規(guī)律。當(dāng)已經(jīng)求得外參數(shù)μ=μn時(shí)的周期解Xn，則第n+1步解的迭代初值為:

而后在μ=μn+1處用打靶法進(jìn)行校準(zhǔn)即可。?H(X，μ)/?μ可以通過(guò)將式(20)對(duì)軌跡Xs(t0+t)積分求得:

當(dāng)δSμ(t0)=0 時(shí)，δSμ(t0+T)=?H(X，μ)/?μ。

3.2 穩(wěn)定性判別的FLOQUET理論

本文周期解判穩(wěn)采用Floquet理論［14］，通過(guò)求解周期解的狀態(tài)轉(zhuǎn)移矩陣即雅克比矩陣J的特征值(Floquet乘子)來(lái)進(jìn)行。當(dāng)所有Floquet乘子均位于復(fù)平面的單位圓內(nèi)時(shí)，周期解穩(wěn)定，而其隨外參數(shù)變化穿越單位圓的不同，周期解會(huì)發(fā)生不同的分叉形式:

(1)當(dāng)一個(gè)模最大的Floquet乘子由(1，0)穿出單位圓時(shí)，周期解失穩(wěn)分叉的可能方式主要有鞍結(jié)分叉、叉型分叉和對(duì)稱破損分叉等多種情況;

(2)當(dāng)一個(gè)模最大的Floquet乘子由(-1，0)穿出單位圓時(shí)，周期解將通過(guò)倍周期分叉而失穩(wěn);

(3)當(dāng)一對(duì)模最大的Floquet乘子以共軛復(fù)數(shù)方式(虛部不為0)穿出單位圓時(shí)，周期解經(jīng)Hopf型偽周期分叉產(chǎn)生偽周期解。

4 計(jì)算實(shí)例及結(jié)果對(duì)比

如圖5所示的一個(gè)典型拉桿組合轉(zhuǎn)子軸承系統(tǒng)，其轉(zhuǎn)子可以看作是由一根整體轉(zhuǎn)軸和固結(jié)于其上的4個(gè)剛性輪盤(pán)，4個(gè)輪盤(pán)通過(guò)3個(gè)接觸界面連接，并由固結(jié)預(yù)緊的12根周向均布拉桿組成，此處預(yù)緊量大小參照工程中重型燃機(jī)選取的拉桿長(zhǎng)度的千分之一。計(jì)入接觸界面以前，系統(tǒng)視作為一根整體轉(zhuǎn)子，無(wú)界面單元存在，此時(shí)系統(tǒng)離散為11段12個(gè)節(jié)點(diǎn);計(jì)入接觸界面以后，輪盤(pán)間界面處理為附加彈簧單元，此時(shí)系統(tǒng)離散為17段18個(gè)節(jié)點(diǎn)。模型材料為密度ρ=7 800 kg/m3，彈性模量E=2.1×1 011的鋼;接觸界面粗糙度0.4 μm(kf=1.72e12，kq=1.41e12)，界面摩擦系數(shù) 0.15，軸徑dshaft=0.08 m、長(zhǎng)度lshaft=1.1 m;界面效應(yīng)計(jì)入前后模型均采用文獻(xiàn)［12］中建立的Pinkus無(wú)限長(zhǎng)軸承模型，它是比較簡(jiǎn)單的軸承力形式之一。

圖5 拉桿組合轉(zhuǎn)子軸承系統(tǒng)結(jié)構(gòu)及單元?jiǎng)澐质疽鈭DFig.5 Structure and finite element of FRRBS

對(duì)于此系統(tǒng)，首先采用前述介紹的方法，建立此系統(tǒng)的動(dòng)力學(xué)模型，輪盤(pán)間接觸界面的附加剛度矩陣通過(guò)相應(yīng)的附加彈簧單元添加到系統(tǒng)中;轉(zhuǎn)子不平衡量的添加方式為在4個(gè)剛性輪盤(pán)上加同相位的無(wú)量綱質(zhì)量偏心距e;徑向軸承力添加到對(duì)應(yīng)軸承節(jié)點(diǎn)上，其主要研究?jī)?nèi)容和結(jié)果整理如下:

4.1 計(jì)入接觸界面對(duì)系統(tǒng)e-ω特性的影響:

圖6為界面計(jì)入前后系統(tǒng)e-ω轉(zhuǎn)遷圖中各類符號(hào)均為計(jì)及界面影響后的標(biāo)記，無(wú)界面時(shí)見(jiàn)文獻(xiàn)［12］。

通過(guò)計(jì)算，得到了計(jì)入接觸界面前后此拉桿組合轉(zhuǎn)子軸承系統(tǒng)穩(wěn)態(tài)解隨參數(shù)ω和e的分布變化規(guī)律，如圖6。

可以看出，計(jì)及接觸界面后系統(tǒng)穩(wěn)態(tài)解的穩(wěn)定性分叉規(guī)律與整體轉(zhuǎn)子相比較，同樣具有以下特征:

同步周期解的Hopf型偽周期分叉集Ⅱ與倍周期分叉集Ⅰ及e=ecr將所研究的參數(shù)域劃分為如圖6的同步周期解①、偽周期解②及倍周期解③;同步周期解的分叉在e較小時(shí)表現(xiàn)為偽周期分叉，而在e較大時(shí)發(fā)生倍周期分叉;不平衡量e較小時(shí)偽周期失穩(wěn)轉(zhuǎn)速隨e增大而略微增大，較大時(shí)倍周期失穩(wěn)轉(zhuǎn)速隨e增大而減小;系統(tǒng)響應(yīng)峰處倍周期分叉集呈鞍形變化。

然而接觸界面的計(jì)入，為系統(tǒng)帶來(lái)的影響也是顯著的:

(1)輪盤(pán)間接觸界面使得整體轉(zhuǎn)子的剛度降低，導(dǎo)致其第一臨界轉(zhuǎn)速明顯下降，如圖7所示波峰所在位置明顯左移，系統(tǒng)穩(wěn)定工作的e-ω參數(shù)區(qū)域整體左移。故在工程實(shí)際中此類組合轉(zhuǎn)子盡量提高輪盤(pán)間接觸界面加工質(zhì)量，以保證足夠的接觸剛度，可以有效降低接觸界面對(duì)臨界轉(zhuǎn)速的影響。

(2)接觸界面使倍周期分叉集Ⅱ和偽周期分叉集Ⅰ的臨界不平量ecr增大，說(shuō)明在第一臨界轉(zhuǎn)速以上工作的轉(zhuǎn)子，較大的不平衡量將使得系統(tǒng)更容易發(fā)生Hopf型偽周期分叉而失穩(wěn)。同時(shí)，當(dāng)系統(tǒng)工作在第一臨界轉(zhuǎn)速后鞍形區(qū)域內(nèi)時(shí)，定轉(zhuǎn)速下發(fā)生倍周期分叉的臨界不平衡量有所增大。此處界面影響一定程度上提高了系統(tǒng)在該e-ω參數(shù)區(qū)域穩(wěn)定運(yùn)行的能力。

(3)在不同參數(shù)區(qū)域，界面的作用力大小不一樣，其對(duì)系統(tǒng)穩(wěn)定性的影響不同。對(duì)高轉(zhuǎn)速區(qū)域失穩(wěn)轉(zhuǎn)速的影響相對(duì)較大，低轉(zhuǎn)速區(qū)域失穩(wěn)轉(zhuǎn)速的影響相對(duì)較小。

(4)由于接觸界面導(dǎo)致整體轉(zhuǎn)子剛度的降低，轉(zhuǎn)子的振動(dòng)量較之整體轉(zhuǎn)子時(shí)有所增大。如圖7，其最大振動(dòng)量在y和x方向上均有顯著增大，對(duì)轉(zhuǎn)子的振動(dòng)特性影響明顯。

4.2 定參數(shù)下帶界面轉(zhuǎn)子系統(tǒng)典型軌跡

此處采用圖 6 中標(biāo)記的 A(11 200 r/min，2 μm)，B(6 500 r/min，6 μm)，C(10000 r/min，11 μm)和D(12 450 r/min，11 μm)四個(gè)典型點(diǎn)的轉(zhuǎn)速和不平衡量參數(shù)，對(duì)比了界面計(jì)入前后系統(tǒng)的典型軌跡。結(jié)果如圖8所示，其中圖 8(a)代表整體轉(zhuǎn)軸的轉(zhuǎn)子軸承系統(tǒng)，圖8(b)代表計(jì)及接觸界面的轉(zhuǎn)子軸承系統(tǒng)。

圖8 ω =11 200 r/min，e=2 μmFig.8ω =11 200 r/min，e=2 μm

圖9 ω =6 500 r/min，e=6 μmFig.9ω =6 500 r/min，e=6 μm

圖10 ω =1 000 r/min，e=11 μmFig.10ω =1 000 r/min，e=11 μm

圖11 ω =12 450 r/min，e=11 μmFig.11ω =12 450 r/min，e=11 μm

圖8表明在高轉(zhuǎn)速區(qū)域小不平衡量時(shí)，整體轉(zhuǎn)軸系統(tǒng)的穩(wěn)態(tài)周期解在接觸界面影響下變?yōu)閭沃芷诮?，且振?dòng)量明顯增大;圖9表明低轉(zhuǎn)速區(qū)域接觸界面使得整體轉(zhuǎn)軸系統(tǒng)由同步周期解變?yōu)楸吨芷诮猓艺駝?dòng)加劇;圖10表明在C點(diǎn)參數(shù)區(qū)域接觸界面使得倍周期運(yùn)行狀況的系統(tǒng)進(jìn)入到穩(wěn)態(tài)同步周期解，一定程度上改善了系統(tǒng)穩(wěn)定性;圖11表明在D點(diǎn)參數(shù)區(qū)域發(fā)生倍周期分叉的整體轉(zhuǎn)軸系統(tǒng)，在計(jì)入接觸界面以后更容易發(fā)生偽周期分叉而導(dǎo)致失穩(wěn)。

5 結(jié)論

本文針對(duì)拉桿組合柔性轉(zhuǎn)子軸承系統(tǒng)，完成了計(jì)及輪盤(pán)間接觸界面影響的動(dòng)力學(xué)系統(tǒng)建模，并對(duì)接觸界面對(duì)其非線性動(dòng)力學(xué)特性的影響進(jìn)行了研究，所得出主要結(jié)論如下:

(1)接觸界面使得整體轉(zhuǎn)子剛度降低，導(dǎo)致第一臨界轉(zhuǎn)速減小，并且使得系統(tǒng)振動(dòng)加劇，最大振動(dòng)幅值明顯增大;

(2)系統(tǒng)穩(wěn)定工作的同步周期解e-ω區(qū)域受接觸界面的影響整體左移，且在高轉(zhuǎn)速區(qū)域的影響量較之低轉(zhuǎn)速區(qū)域大;

(3)在高轉(zhuǎn)速區(qū)較大不平衡量時(shí)，轉(zhuǎn)子系統(tǒng)的分叉形式可能受界面影響，由倍周期分叉失穩(wěn)變?yōu)槲恢芷诜植媸Х€(wěn)。

［1］Lund J W，Nelson H B.Instability threshold of an unbalanced rigid rotor in short journal bearings［C］//IMechE，Proceeding IMechE Conference on Vibration in Rotating Machinery，1980，Cambridge，UK.IMechE，1980:91-95.

［2］ Kim Y B，Noahst S T.Bifurcation analysis for a modified jeffcott rotor with bearing clearance［J］.Nonlinear Dynamics，1990，1(3):221-241.

［3］Sundararajan P，Noah S T.Dynamics of forced nonlinear systems using shooting/arclength continuation methodapplication to rotor system［J］.ASME Journal of Vibration and Acoustics，1997，119(1):10-20.

［4］ Zheng T，Hasebe N.Nonlinear dynamic behaviors of a complex rotor-bearing system［J］.ASME Journal of Applied Mechanics，2000，67(3):485-495.

［5］汪光明，饒柱石，夏松波.拉桿轉(zhuǎn)子力學(xué)模型的研究［J］.航空學(xué)報(bào)，1993，14(8):419-423.

WANG Guang-ming， RAO Zhu-shi， XIA Song-bo. The analysis of mechanical model of rod fastening rotor［J］.Acta Aeronautica et Astronautica Sinica，1993，14(8):419-423.

［6］饒柱石.拉桿組合式特種轉(zhuǎn)子力學(xué)特性及其接觸剛度的研究［D］.哈爾濱:哈爾濱工業(yè)大學(xué)，1992.

［7］王艾倫，駱 舟.拉桿轉(zhuǎn)子軸向振動(dòng)的動(dòng)力學(xué)模型［J］.機(jī)械工程學(xué)報(bào)，2009，20(13):1524-1527.

WANG Ai-lun，LUO Zhou.Research on rod fastening rotor dynamicsaxialvibration［J］. JournalofMechanical Engineering，2009，20(13):1524-1527.

［8］何 鵬，劉占生，張廣輝，等.分布式拉桿轉(zhuǎn)子動(dòng)力學(xué)建模與分析［J］.汽輪機(jī)技術(shù)，2010，52(1):4-9.

HE Peng， LIU Zhan-sheng， ZHANG Guang-hui， et al.Dynamic modeling and analysis of distributed rod fastening rotor［J］.Turbine Technology，2010，52(1):4-9.

［9］李輝光，劉 恒，虞 烈.粗糙機(jī)械結(jié)合面的接觸剛度研究［J］.西安交通大學(xué)學(xué)報(bào)，2011，45(6):1-6.

LI Hui-guang，LIU Heng，YU Lie.Study on contact stiffness of rough mechanical joint surface［J］.Journal of xi'an Jiaotong University，2011，45(6):1-6.

［10］Valadez G M，Baltazar A，Joyce D R S.Study of interfacial stiffness ratio of a rough in contact using a spring model［J］.Wear，2010，268(3-4):373-379.

［11］ Nelson H D.A finite rotating shaft element using timoshenko beam theory［J］.ASME J.Mech.Des.，1980，102:793-803.

［12］劉 恒，陳 麗.軸向均布拉桿柔性組合轉(zhuǎn)子軸承系統(tǒng)的非線性動(dòng)力特性［J］.機(jī)械工程學(xué)報(bào)，2010，46(19):53-62.

LIU Heng， CHEN Li. Nonlinear dynamic analysis of a flexible rod fastening rotor bearing system［J］.Journal of Mechanical Engineering，2010，46(19):52-62.

［13］Fey R H B，Van Campen D H，De Kraker A.Long term structuraldynamics ofmechanical systems with local nonlineariyies［J］.ASME Journal of Vibration and Acoustics，1996，118(4):147-153.

［14］Ho Y S，Liu H，Yu L.Effect of thrust magnetic bearing on stability and bifurcation of a flexible rotor active magnetic bearing system［J］. ASME JournalofVibration and Acoustics，2003，125(3):307-316.