接觸界面對拉桿組合柔性轉子軸承系統的非線性動力特性影響

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

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

現代燃氣輪機、航空發動機等大型動力系統的核心部件均為典型的拉桿柔性組合轉子。此類轉子的多個獨立輪盤間依靠接觸界面傳遞各向作用力，并在拉桿預緊下組合為一體。重量輕、易于裝配且具有良好的冷卻效果是此類轉子的特點，針對高速化、高效化的工業發展，應用前景非常廣泛。

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

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

1 接觸界面建模

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

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

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

1.1 彈簧單元剛度確定

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

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

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

1.2 附加彈簧單元建立

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

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

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

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

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

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

Δx2取一階近似后可得:

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

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

在軸向z方向上，有:

可知:

其中:

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

2 系統方程

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

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

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

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

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

3 分析方法

假設作用于系統的載荷為周期為T的周期載荷:

則系統的穩態解為周期解、偽周期解及混沌解。對于軸承轉子系統而言，工頻周期解是其本征解，隨外參數轉速ω、不平衡量e等的變化，此解將發生失穩，而產生新的穩態解形式，如偽周期解及混沌解，或發生絕對失穩，因此求取系統的周期解及其穩定性分叉規律便成為研究此類系統非線性特性的主要內容。

3.1 周期解求解

應用打靶法求解式(13)動力系統周期解的問題可轉化為兩點邊值問題:

其中:

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

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

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

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

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

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

3.2 穩定性判別的FLOQUET理論

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

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

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

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

4 計算實例及結果對比

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

圖5 拉桿組合轉子軸承系統結構及單元劃分示意圖Fig.5 Structure and finite element of FRRBS

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

4.1 計入接觸界面對系統e-ω特性的影響:

圖6為界面計入前后系統e-ω轉遷圖中各類符號均為計及界面影響后的標記，無界面時見文獻［12］。

通過計算，得到了計入接觸界面前后此拉桿組合轉子軸承系統穩態解隨參數ω和e的分布變化規律，如圖6。

可以看出，計及接觸界面后系統穩態解的穩定性分叉規律與整體轉子相比較，同樣具有以下特征:

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

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

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

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

(3)在不同參數區域，界面的作用力大小不一樣，其對系統穩定性的影響不同。對高轉速區域失穩轉速的影響相對較大，低轉速區域失穩轉速的影響相對較小。

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

4.2 定參數下帶界面轉子系統典型軌跡

此處采用圖 6 中標記的 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)四個典型點的轉速和不平衡量參數，對比了界面計入前后系統的典型軌跡。結果如圖8所示，其中圖 8(a)代表整體轉軸的轉子軸承系統，圖8(b)代表計及接觸界面的轉子軸承系統。

圖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表明在高轉速區域小不平衡量時，整體轉軸系統的穩態周期解在接觸界面影響下變為偽周期解，且振動量明顯增大;圖9表明低轉速區域接觸界面使得整體轉軸系統由同步周期解變為倍周期解，且振動加劇;圖10表明在C點參數區域接觸界面使得倍周期運行狀況的系統進入到穩態同步周期解，一定程度上改善了系統穩定性;圖11表明在D點參數區域發生倍周期分叉的整體轉軸系統，在計入接觸界面以后更容易發生偽周期分叉而導致失穩。

5 結論

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

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

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

(3)在高轉速區較大不平衡量時，轉子系統的分叉形式可能受界面影響，由倍周期分叉失穩變為位周期分叉失穩。

［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］汪光明，饒柱石，夏松波.拉桿轉子力學模型的研究［J］.航空學報，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］饒柱石.拉桿組合式特種轉子力學特性及其接觸剛度的研究［D］.哈爾濱:哈爾濱工業大學，1992.

［7］王艾倫，駱 舟.拉桿轉子軸向振動的動力學模型［J］.機械工程學報，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］何 鵬，劉占生，張廣輝，等.分布式拉桿轉子動力學建模與分析［J］.汽輪機技術，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］.西安交通大學學報，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］劉 恒，陳 麗.軸向均布拉桿柔性組合轉子軸承系統的非線性動力特性［J］.機械工程學報，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.