馬 斌賴 歡呂 波侯 予

（1.中國(guó)空氣動(dòng)力學(xué)研究與發(fā)展中心；2.西安交通大學(xué)）

0　引言

隨著工業(yè)技術(shù)的發(fā)展，要求旋轉(zhuǎn)機(jī)械的工作轉(zhuǎn)速和效率不斷提高。轉(zhuǎn)子-支承系統(tǒng)的穩(wěn)定性是設(shè)計(jì)轉(zhuǎn)子系統(tǒng)的關(guān)鍵。當(dāng)系統(tǒng)工作轉(zhuǎn)速接近系統(tǒng)的橫向振動(dòng)頻率時(shí)，轉(zhuǎn)子會(huì)產(chǎn)生強(qiáng)烈的振動(dòng)和噪聲，嚴(yán)重時(shí)會(huì)造成系統(tǒng)損毀[1]。因此進(jìn)行轉(zhuǎn)子動(dòng)力學(xué)分析成為系統(tǒng)設(shè)計(jì)階段的重要內(nèi)容之一?，F(xiàn)代的轉(zhuǎn)子動(dòng)力學(xué)計(jì)算方法可以分為兩大類：傳遞矩陣法和有限元法。本文針對(duì)基于傳遞矩陣算法的商業(yè)軟件DYNAMICS進(jìn)行了介紹，對(duì)兩種不同結(jié)構(gòu)的透平膨脹機(jī)轉(zhuǎn)子進(jìn)行了分析，與通用有限元軟件ANSYS計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比研究。

1　轉(zhuǎn)子動(dòng)力學(xué)分析的基本理論

一個(gè)機(jī)械系統(tǒng)的運(yùn)動(dòng)微分方程可表達(dá)為

式中：M，C，K為系統(tǒng)的質(zhì)量、阻尼和剛度矩陣；Z為系統(tǒng)的廣義坐標(biāo)矢量；F是作用在系統(tǒng)上的廣義外力。對(duì)于轉(zhuǎn)子系統(tǒng)，還需要考慮轉(zhuǎn)動(dòng)帶來(lái)的陀螺效應(yīng)、轉(zhuǎn)動(dòng)剛度和阻尼。因此轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)微分方程可表達(dá)為

式中：C是阻尼矩陣，非對(duì)稱陣；G是陀螺矩陣，反對(duì)稱陣；K是剛度矩陣的對(duì)稱部分；S是它的不對(duì)稱部分。上述矩陣均是與轉(zhuǎn)速ω有關(guān)的函數(shù)。

傳遞矩陣法和有限元法均基于上述理論，求解式（2）的特征值來(lái)求解轉(zhuǎn)子的臨界轉(zhuǎn)速。傳遞矩陣法是對(duì)各個(gè)軸段的截面狀態(tài)向量進(jìn)行逐段推算，直到計(jì)算到轉(zhuǎn)軸的另一端，主要特點(diǎn)是：矩陣的階數(shù)不隨系統(tǒng)的自由度數(shù)增大而增加，因而編程簡(jiǎn)單，占內(nèi)存少，運(yùn)算速度快[2]。有限元法是將連續(xù)系統(tǒng)離散化以變分原理為基礎(chǔ)形成的一種數(shù)值近似解法[3]，隨著轉(zhuǎn)子有限元模型的不斷完善以及計(jì)算機(jī)的發(fā)展，使得有限元可以對(duì)大型復(fù)雜結(jié)構(gòu)進(jìn)行數(shù)值計(jì)算[4]。但系統(tǒng)復(fù)雜時(shí)會(huì)導(dǎo)致自由度數(shù)特別大，耗費(fèi)計(jì)算時(shí)間[5]。

2　DYNAMICS臨界轉(zhuǎn)速計(jì)算

DYNAMICS由俄羅斯頂級(jí)透平機(jī)械工業(yè)領(lǐng)袖——Alfa-Tranzit公司研發(fā)，其核心技術(shù)團(tuán)隊(duì)來(lái)自國(guó)際久負(fù)盛名的莫斯科航空航天學(xué)院（AMI），合作單位涉及俄羅斯SU系列及米格系列戰(zhàn)機(jī)引擎的制造商。DYNAMICS是一款針對(duì)于旋轉(zhuǎn)機(jī)械設(shè)計(jì)與分析的轉(zhuǎn)子動(dòng)力學(xué)軟件，致力于高性能旋轉(zhuǎn)機(jī)械轉(zhuǎn)子系統(tǒng)的振動(dòng)、平衡和穩(wěn)定科學(xué)研究，可以實(shí)現(xiàn)轉(zhuǎn)子系統(tǒng)的結(jié)構(gòu)建模、仿真分析和故障診斷。它的分析設(shè)計(jì)對(duì)象從簡(jiǎn)單的Jeffcott轉(zhuǎn)子到多軸系轉(zhuǎn)子、齒輪系統(tǒng)，再到復(fù)雜的機(jī)匣-多轉(zhuǎn)子系統(tǒng)，能夠模擬簡(jiǎn)單的剛性支承、正交彈性支承和復(fù)雜的滾動(dòng)軸承、滑動(dòng)軸承、擠壓油膜阻尼器（SFD）支撐等，同時(shí)能夠提供長(zhǎng)短軸承理論、半膜假設(shè)、全膜假設(shè)等理論，此外還可以對(duì)外界阻尼、材料內(nèi)阻、間隙、碰撞等因素加以考慮[6]。圖1為DYNAMICS界面。

圖1　DYNAMICS界面圖Fig.1　Interface of DYNAMICS

2.1　前處理

DYNAMICS的數(shù)據(jù)結(jié)構(gòu)由子系統(tǒng)、轉(zhuǎn)配體和子模型三級(jí)組成，具有多級(jí)轉(zhuǎn)子結(jié)構(gòu)完整模型的文本與圖形并存的前處理能力；提供標(biāo)準(zhǔn)的建模單元庫(kù)，包括梁、殼、集中質(zhì)量、耦合、支撐、齒輪、超單元等；提供標(biāo)準(zhǔn)類型的材料庫(kù)；可以施加的載荷包括：不平衡力、一般形式的動(dòng)態(tài)力、葉片飛出、基礎(chǔ)運(yùn)動(dòng)等；支持轉(zhuǎn)子模型的2D和3D可視化。

2.2　求解計(jì)算

軟件使用傳遞矩陣法進(jìn)行轉(zhuǎn)子動(dòng)力學(xué)求解，計(jì)算速度快、精度能夠滿足工程需求?？梢赃x擇的求解模塊包括：自振頻率、臨界轉(zhuǎn)速、不平衡響應(yīng)、穩(wěn)定性、瞬態(tài)響應(yīng)等。

2.3　后處理

軟件可以輸出位移、速率、加速度、反作用力和力矩、能量等計(jì)算結(jié)果，允許對(duì)時(shí)間信號(hào)進(jìn)行充分的分析，可以得到平均值、軌跡、傅立葉變換、瀑布圖等計(jì)算結(jié)果。

3　計(jì)算結(jié)果

本文以箔片動(dòng)壓氣體軸承支撐的直徑Φ1=25mm雙盤轉(zhuǎn)子透平膨脹機(jī)和靜壓氣體軸承支撐的Φ2=17mm單盤轉(zhuǎn)子透平膨脹機(jī)為研究對(duì)象，轉(zhuǎn)子參數(shù)見表1，轉(zhuǎn)子結(jié)構(gòu)如圖2所示。動(dòng)壓軸承采用鈹青銅（QBe1.7）和不銹鋼（0Cr18Ni9）兩種材質(zhì)的雙層鼓泡箔片動(dòng)壓軸承，軸承長(zhǎng)度L1=36mm，軸承名義直徑D1=25.04mm，設(shè)計(jì)半徑間隙C1=20μm，剛度與阻尼根據(jù)文獻(xiàn)[7]計(jì)算結(jié)果選取，如圖3所示。靜壓軸承采用雙排環(huán)面切向小孔供氣形式，每排8個(gè)小孔均布，軸承長(zhǎng)度L2=27mm，軸承直徑D2=17.04mm，設(shè)計(jì)半徑間隙C2=20μm，節(jié)流孔直徑d=0.3mm，節(jié)流孔距軸承端面距離l=6.5mm，供氣壓力0.45MPa，剛度和阻尼如圖3所示。從圖3可以看出，切向供氣具有較小的交叉剛度和較大的主阻尼，因此具有較好的穩(wěn)定性[8]。

表1　轉(zhuǎn)子參數(shù)表Tab.1　Parameters of different rotors

圖2　轉(zhuǎn)子結(jié)構(gòu)圖（單位：mm）Fig.2　structure of different rotors（unit:mm）

圖3　軸承剛度和阻尼圖Fig.3　Variations of the stiffness coefficients and the damping coefficients versus rotational speed

圖4為DYNAMICS分析獲得不銹鋼雙層鼓泡箔片軸承支撐的Φ1=25mm轉(zhuǎn)子系統(tǒng)Campbell圖?？梢?，該系統(tǒng)呈現(xiàn)出高頻渦動(dòng)頻率與低頻渦動(dòng)頻率兩極分化的現(xiàn)象。低頻渦動(dòng)共4階，轉(zhuǎn)速均在1.5×104rpm（250Hz）以下，而高頻渦動(dòng)（本文只給出了其中的前2階）的轉(zhuǎn)速都在18×104rpm（3000Hz）以上。圖5為Φ1=25mm轉(zhuǎn)子的振型圖，8 499rpm，10 182rpm，186 466rpm對(duì)應(yīng)反向渦動(dòng)，11 698rpm，14 131rpm，208 858rpm對(duì)應(yīng)正向渦動(dòng)。前4階為剛體運(yùn)動(dòng)模態(tài)，5，6兩階為彎曲模態(tài)。造成以上結(jié)果的原因是轉(zhuǎn)子剛度遠(yuǎn)大于軸承剛度，使轉(zhuǎn)子彎曲模態(tài)的頻率遠(yuǎn)高于轉(zhuǎn)子剛體運(yùn)動(dòng)模態(tài)的頻率。低頻的前4階臨界轉(zhuǎn)速與高頻的5，6階臨界轉(zhuǎn)速之間較寬的區(qū)域意味著轉(zhuǎn)子系統(tǒng)具有較寬的安全工作轉(zhuǎn)速范圍?？梢?.5×104rpm到18×104rpm之間都是該轉(zhuǎn)子系統(tǒng)的安全工作轉(zhuǎn)速。18×104rpm這一轉(zhuǎn)速已經(jīng)遠(yuǎn)遠(yuǎn)高于透平實(shí)際工作轉(zhuǎn)速，因此在理論上透平只需要順利跨越前四階臨界轉(zhuǎn)速便能在較寬的轉(zhuǎn)速范圍內(nèi)穩(wěn)定工作。

圖4　不銹鋼氣浮軸承支撐的Φ1=25mm轉(zhuǎn)子系統(tǒng)Campbell圖Fig.4　Campbell diagram of the rotor system supported by stainless steel foil bearings(Φ1=25mm)

圖5　不銹鋼軸承支撐的Φ1=25mm轉(zhuǎn)子振型圖Fig.5　Mode sharp of the rotor system supported by stainless steel foil bearings(Φ1=25mm)

圖6為DYNAMICS分析獲得靜壓軸承支撐的Φ2=17mm轉(zhuǎn)子系統(tǒng)Campbell圖。該系統(tǒng)同樣呈現(xiàn)出高頻渦動(dòng)頻率與低頻渦動(dòng)頻率兩極分化的現(xiàn)象。低頻渦動(dòng)共4階，轉(zhuǎn)速均在2.5×104rpm（416.7Hz）以下，而高頻渦動(dòng)（本文給出了其中的前2階）的轉(zhuǎn)速都在15×104rpm（2 500Hz）以上。圖7為Φ2=17mm轉(zhuǎn)子的振型圖，其中17 343rpm，23 482rpm，156 145rpm對(duì)應(yīng)反向渦動(dòng)，18 053rpm，24 786rpm，173 205rpm對(duì)應(yīng)正向渦動(dòng)。前4階為剛體運(yùn)動(dòng)模態(tài)，5，6兩階為彎曲模態(tài)。造成以上結(jié)果的原因同樣是轉(zhuǎn)子剛度大于軸承剛度，使轉(zhuǎn)子彎曲模態(tài)的頻率高于轉(zhuǎn)子剛體運(yùn)動(dòng)模態(tài)的頻率。與Φ1=25mm轉(zhuǎn)子系統(tǒng)相比，靜壓軸承剛度和阻尼大于動(dòng)壓軸承，而其安全工作轉(zhuǎn)速在2.5×104rpm到15×104rpm之間，明顯小于Φ1=25mm轉(zhuǎn)子系統(tǒng)。這是由于Φ2=17mm轉(zhuǎn)子為單止推盤偏置，屬于細(xì)長(zhǎng)結(jié)構(gòu)，剛度較低。文獻(xiàn)[9、10]分析結(jié)果表明相同軸徑（Φ=25mm）條件下，這種轉(zhuǎn)子低階臨界轉(zhuǎn)速對(duì)軸承剛度和阻尼敏感度較大，穩(wěn)定性較差，彎曲模態(tài)的頻率較低，工作范圍較小。

圖6　靜壓軸承支撐Φ2=17mm轉(zhuǎn)子系統(tǒng)Campbell圖Fig.6　Campbell diagram of the rotor system supported by hydrostatics gas bearings(Φ2=17mm)

圖7　Φ2=17mm轉(zhuǎn)子振型圖Fig.7　Mode sharp of the rotor system supported by hydrostatics gas bearings(Φ2=27mm)

4　理論計(jì)算結(jié)果對(duì)比

表2和表3分別給出了不同軸承支撐的Φ1=25mm轉(zhuǎn)子前6階臨界轉(zhuǎn)速ANSYS與DYNAMICS的計(jì)算結(jié)果，可以看出前4階臨界轉(zhuǎn)速的計(jì)算結(jié)果差異較大，臨界轉(zhuǎn)速提高差異變小，第5、6階臨界轉(zhuǎn)速的計(jì)算結(jié)果差異小于1%。ANSYS計(jì)算時(shí)間大于100s，DYNAMICS計(jì)算時(shí)間小于1s。造成低階臨界轉(zhuǎn)速計(jì)算結(jié)果差異較大的原因是由于ANSYS與DYNAMICS的插值方式不同，當(dāng)轉(zhuǎn)速小于40 000rpm時(shí)，ANSYS采用40 000rpm時(shí)的剛度和阻尼，當(dāng)轉(zhuǎn)速大于75 000rpm時(shí)，ANSYS采用75 000rpm時(shí)的剛度和阻尼；當(dāng)轉(zhuǎn)速超出給定范圍時(shí)，DYNAMICS采用內(nèi)部擬合值。因此低轉(zhuǎn)速時(shí)ANSYS所選剛度和阻尼較大，高轉(zhuǎn)速時(shí)ANSYS所選剛度和阻尼較小。因?yàn)檗D(zhuǎn)子剛體運(yùn)動(dòng)模態(tài)頻率取決于支撐的剛度，軸承剛度增加，轉(zhuǎn)子的臨界轉(zhuǎn)速提高[11,12]，而轉(zhuǎn)子彎曲模態(tài)頻率依賴轉(zhuǎn)子剛度，所以出現(xiàn)上述差異。同理，由于不銹鋼軸承剛度大于鈹青銅軸承剛度，因此不銹鋼軸承支撐轉(zhuǎn)子的前4階臨界轉(zhuǎn)速大于相對(duì)應(yīng)鈹青銅軸承支撐轉(zhuǎn)子的前4階臨界轉(zhuǎn)速，而兩者的第5、6階臨界轉(zhuǎn)速基本相同。

表2　鈹青銅軸承支撐Φ1=25mm轉(zhuǎn)子計(jì)算結(jié)果Tab.2　Calculation result of the rotor system supported by beryllium bronze foil bearings

表3　不銹鋼軸承支撐Φ1=25mm轉(zhuǎn)子計(jì)算結(jié)果Tab.3　Calculation result of the rotor system supported by stainless steel foil bearings

表4給出了Φ2=17mm轉(zhuǎn)子前6階臨界轉(zhuǎn)速ANSYS與DYNAMICS的計(jì)算結(jié)果，可以看出各階臨界轉(zhuǎn)速差異均小于2%。ANSYS計(jì)算時(shí)間為70s，DYNAMICS計(jì)算時(shí)間小于1s。對(duì)比表2～4可知，ANSYS與DYNAMICS計(jì)算結(jié)果，動(dòng)壓軸承支撐轉(zhuǎn)子系統(tǒng)前4階臨界轉(zhuǎn)速差異較大，靜壓軸承支撐轉(zhuǎn)子系統(tǒng)前4階臨界轉(zhuǎn)速差異較小。從圖3（a）可以看出，動(dòng)壓軸承主剛度Kxx，Kyy隨轉(zhuǎn)速變化較大，靜壓軸承主剛度Kxx，Kyy隨隨轉(zhuǎn)速變化較小，這是因?yàn)閯?dòng)壓軸承剛度主要依賴轉(zhuǎn)子轉(zhuǎn)速，靜壓軸承剛度主要依賴于軸承供氣壓力。如前文所述，ANSYS與DYNAMICS的插值方式不同，低轉(zhuǎn)速時(shí)，靜壓軸承剛度插值差異小，因此靜壓軸承支撐轉(zhuǎn)子系統(tǒng)低階臨界轉(zhuǎn)速計(jì)算結(jié)果差異較小。

表4　靜壓軸承支撐Φ2=17mm轉(zhuǎn)子計(jì)算結(jié)果Tab.4　Calculation result of the rotor supported by hydrostatics gas bearings

5　理論計(jì)算與實(shí)驗(yàn)結(jié)果對(duì)比

文獻(xiàn)[7]對(duì)于鈹青銅軸承和不銹鋼軸承支撐的Φ1=25mm轉(zhuǎn)子系統(tǒng)的自然頻率預(yù)測(cè)值分別為9546rpm（159.1Hz）和11 388rpm（189.8Hz），鈹青銅軸承支撐轉(zhuǎn)子非同步渦動(dòng)振動(dòng)頻率實(shí)驗(yàn)值分布在150Hz（9 000rpm）～200Hz（12 000rpm）之間，不銹鋼軸承支撐轉(zhuǎn)子非同步渦動(dòng)振動(dòng)頻率實(shí)驗(yàn)值分布在190Hz（11 400rpm）～210Hz（12 600rpm）之間。本文鈹青銅軸承和不銹鋼軸承支撐轉(zhuǎn)子系統(tǒng)3，4階臨界轉(zhuǎn)速計(jì)算結(jié)果9 394rpm、11 289rpm和11 698rpm、14 131rpm與非同步渦動(dòng)振動(dòng)頻率實(shí)驗(yàn)結(jié)果基本吻合。1，2兩階臨界轉(zhuǎn)速缺乏實(shí)驗(yàn)數(shù)據(jù)，高階臨界轉(zhuǎn)速已超出轉(zhuǎn)子工作范圍，無(wú)法驗(yàn)證。

圖8為Φ2=17mm轉(zhuǎn)子軸心軌跡變化圖，可以看出隨著轉(zhuǎn)速提高，轉(zhuǎn)子軸心軌跡變化明顯。轉(zhuǎn)子轉(zhuǎn)動(dòng)方向?yàn)槟鏁r(shí)針，轉(zhuǎn)速為20 660rpm時(shí)，轉(zhuǎn)子渦動(dòng)方向?yàn)轫槙r(shí)針（反進(jìn)），轉(zhuǎn)速為21 124rpm時(shí)，轉(zhuǎn)子軸心軌跡出現(xiàn)扭曲，是因?yàn)檗D(zhuǎn)子渦動(dòng)發(fā)生轉(zhuǎn)向。其他轉(zhuǎn)速時(shí)轉(zhuǎn)子渦動(dòng)方向?yàn)槟鏁r(shí)針（正進(jìn)）。轉(zhuǎn)速超過(guò)21 949rpm后，轉(zhuǎn)子軸心軌跡形狀保持穩(wěn)定，表明轉(zhuǎn)子跨過(guò)臨界轉(zhuǎn)速進(jìn)入穩(wěn)定運(yùn)行。

圖9為Φ2=17mm轉(zhuǎn)子振幅、相位隨轉(zhuǎn)速變化圖，從圖9（a）中可以看出轉(zhuǎn)子X，Y方向的振幅和相位從10 000rpm到30 000rpm時(shí)變化劇烈，轉(zhuǎn)速超過(guò)30 000rpm后振幅和相位保持穩(wěn)定。從圖9（b）中可以看出在轉(zhuǎn)速19 561rpm、22 948rpm時(shí)轉(zhuǎn)子X方向振幅有極大值，轉(zhuǎn)速15 263rpm、20 660rpm和21 656rpm時(shí)轉(zhuǎn)子Y方向振幅有極大值。X方向振動(dòng)相位在19561rpm～20660rpm之間跨過(guò)90°。Y方向振動(dòng)相位分別在15 263rpm～17 770rpm、19 561rpm～20 660rpm、21 124～21 656rpm 之 間 跨 過(guò)180°。在轉(zhuǎn)速為20660rpm時(shí)，X方向振動(dòng)相位大于Y方向振動(dòng)的相位，因此出現(xiàn)圖7中渦動(dòng)與轉(zhuǎn)子轉(zhuǎn)動(dòng)反向的現(xiàn)象。由文獻(xiàn)[13-15]可知，發(fā)生相位偏轉(zhuǎn)的轉(zhuǎn)速為轉(zhuǎn)子的臨界轉(zhuǎn)速，但本文中轉(zhuǎn)子為渦輪驅(qū)動(dòng)，啟動(dòng)階段轉(zhuǎn)速很難穩(wěn)定在特定轉(zhuǎn)速，因此無(wú)法確定相位偏轉(zhuǎn)速度。綜合振幅與相位判斷轉(zhuǎn)子的臨界轉(zhuǎn)速為：17 000rpm、19 500rpm、20 500rpm、21 500rpm。計(jì)算結(jié)果分別為 17 343rpm、18 053rpm、23 483rpm、24 786rpm，計(jì)算結(jié)果略大于實(shí)驗(yàn)結(jié)果，造成差異的原因是計(jì)算時(shí)低轉(zhuǎn)速的剛度和阻尼由插值得到，與軸承實(shí)際剛度和阻尼有差距，想要得到準(zhǔn)確的臨界轉(zhuǎn)速需要根據(jù)實(shí)際的軸承剛度和阻尼計(jì)算。出現(xiàn)X，Y方向振動(dòng)幅值與相位差異的原因一方面是軸承本身為各向異性，另一方面是由于Y方向有重力作用，造成軸承在X，Y方向的受力以及剛度和阻尼差異較大，從而造成X，Y方向的振動(dòng)特性有較大區(qū)別。

圖8　Φ2=17mm轉(zhuǎn)子低轉(zhuǎn)速軸心軌跡變化圖Fig.8　Orbit of shaft center of the rotor at low speed(Φ2=17mm)

圖9　Φ2=17mm轉(zhuǎn)子振幅、相位圖Fig.9　Amplitude and phase position of the rotor supported by hydrostatics gas bearings(Φ2=17mm)

6　結(jié)論

本文通過(guò)比較ANSYS和DYNAMICS計(jì)算結(jié)果，并與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比分析，得出以下結(jié)論：

1）ANSYS作為通用有限元商業(yè)軟件，軸承剛度和阻尼超出給定范圍后為定值，在轉(zhuǎn)子動(dòng)力學(xué)分析上具有一定的局限性，計(jì)算耗時(shí)較長(zhǎng)；DYNAMICS作為專業(yè)的轉(zhuǎn)子動(dòng)力學(xué)分析軟件，軸承剛度和阻尼超出給定范圍后可線性插值，操作簡(jiǎn)單、運(yùn)算速度快、適用范圍廣。

2）轉(zhuǎn)子的低階臨界轉(zhuǎn)速隨軸承剛度增大而提高，高階臨界轉(zhuǎn)速依賴于轉(zhuǎn)子剛度，基本保持不變。

3）在轉(zhuǎn)子動(dòng)力學(xué)分析中，為了得到準(zhǔn)確的計(jì)算結(jié)果，需提供不同轉(zhuǎn)速下軸承的剛度和阻尼，尤其是臨界轉(zhuǎn)速附近轉(zhuǎn)速時(shí)軸承的剛度和阻尼。