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

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

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

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

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

1 接觸界面建模

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

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

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

1.1 彈簧單元剛度確定

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

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

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

1.2 附加彈簧單元建立

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

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

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

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

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

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

Δx2取一階近似后可得:

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

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

在軸向z方向上，有:

可知:

其中:

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

2 系統(tǒng)方程

對于預(yù)緊后的組合轉(zhuǎn)子整體轉(zhuǎn)軸，直接采用了計及軸向力的鐵木辛格梁軸單元進(jìn)行離散［11］，如圖4所示。在整體轉(zhuǎn)軸的動力方程中，根據(jù)哈密頓原理加入輪盤間接觸界面、剛性圓盤和周向均布拉桿［12］的影響后，最終給出預(yù)緊后的拉桿組合轉(zhuǎn)子軸承系統(tǒ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分別為剛性圓盤的質(zhì)量、陀螺矩陣及重力和不平衡力矢量;foil為軸承力;系統(tǒng)的位移矢量為:xS={x1y1ψ1φ1…xmymψmφm}T。

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

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

3 分析方法

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

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

3.1 周期解求解

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

其中:

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

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

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

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

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

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

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

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

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

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

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

4 計算實例及結(jié)果對比

如圖5所示的一個典型拉桿組合轉(zhuǎn)子軸承系統(tǒng)，其轉(zhuǎn)子可以看作是由一根整體轉(zhuǎn)軸和固結(jié)于其上的4個剛性輪盤，4個輪盤通過3個接觸界面連接，并由固結(jié)預(yù)緊的12根周向均布拉桿組成，此處預(yù)緊量大小參照工程中重型燃機(jī)選取的拉桿長度的千分之一。計入接觸界面以前，系統(tǒng)視作為一根整體轉(zhuǎn)子，無界面單元存在，此時系統(tǒng)離散為11段12個節(jié)點;計入接觸界面以后，輪盤間界面處理為附加彈簧單元，此時系統(tǒng)離散為17段18個節(jié)點。模型材料為密度ρ=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、長度lshaft=1.1 m;界面效應(yīng)計入前后模型均采用文獻(xiàn)［12］中建立的Pinkus無限長軸承模型，它是比較簡單的軸承力形式之一。

圖5 拉桿組合轉(zhuǎn)子軸承系統(tǒng)結(jié)構(gòu)及單元劃分示意圖Fig.5 Structure and finite element of FRRBS

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

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

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

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

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

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

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

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

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

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

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

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)四個典型點的轉(zhuǎn)速和不平衡量參數(shù)，對比了界面計入前后系統(tǒng)的典型軌跡。結(jié)果如圖8所示，其中圖 8(a)代表整體轉(zhuǎn)軸的轉(zhuǎn)子軸承系統(tǒng)，圖8(b)代表計及接觸界面的轉(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ū)域小不平衡量時，整體轉(zhuǎn)軸系統(tǒng)的穩(wěn)態(tài)周期解在接觸界面影響下變?yōu)閭沃芷诮?，且振動量明顯增大;圖9表明低轉(zhuǎn)速區(qū)域接觸界面使得整體轉(zhuǎn)軸系統(tǒng)由同步周期解變?yōu)楸吨芷诮?，且振動加?圖10表明在C點參數(shù)區(qū)域接觸界面使得倍周期運(yùn)行狀況的系統(tǒng)進(jìn)入到穩(wěn)態(tài)同步周期解，一定程度上改善了系統(tǒng)穩(wěn)定性;圖11表明在D點參數(shù)區(qū)域發(fā)生倍周期分叉的整體轉(zhuǎn)軸系統(tǒng)，在計入接觸界面以后更容易發(fā)生偽周期分叉而導(dǎo)致失穩(wěn)。

5 結(jié)論

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

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

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

(3)在高轉(zhuǎn)速區(qū)較大不平衡量時，轉(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é)報，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)子軸向振動的動力學(xué)模型［J］.機(jī)械工程學(xué)報，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)子動力學(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é)報，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)的非線性動力特性［J］.機(jī)械工程學(xué)報，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.