多支承轉(zhuǎn)子系統(tǒng)辛空間傳遞矩陣法及應(yīng)用

張娟娟 ， 崔 升， 馮永新

(1. 復(fù)旦大學(xué) 航空航天系，上海 200433； 2.廣東電網(wǎng)公司 電力科學(xué)研究院，廣州 510600)

多支承轉(zhuǎn)子系統(tǒng)辛空間傳遞矩陣法及應(yīng)用

張娟娟1， 崔 升1， 馮永新2

(1. 復(fù)旦大學(xué) 航空航天系，上海 200433； 2.廣東電網(wǎng)公司 電力科學(xué)研究院，廣州 510600)

由變分原理出發(fā)，引入對(duì)偶變量，將多支承轉(zhuǎn)子系統(tǒng)導(dǎo)入辛體系，推導(dǎo)出轉(zhuǎn)子系統(tǒng)在辛空間中的傳遞辛矩陣格式,體現(xiàn)了保辛的性質(zhì)。建立多支承轉(zhuǎn)子系統(tǒng)模型，運(yùn)用MATLAB對(duì)模型進(jìn)行計(jì)算分析，求得固有頻率等模態(tài)參數(shù)，計(jì)算結(jié)果與有限元方法計(jì)算值相吻合。通過(guò)與傳統(tǒng)方法的計(jì)算結(jié)果比較證明了該計(jì)算方法的有效性。求解得到轉(zhuǎn)子系統(tǒng)各階臨界轉(zhuǎn)速，探討了支承剛度的改變對(duì)轉(zhuǎn)子各階臨界轉(zhuǎn)速的影響，為多支承轉(zhuǎn)子系統(tǒng)的優(yōu)化設(shè)計(jì)提供了依據(jù)參考。

轉(zhuǎn)子系統(tǒng)；辛空間；傳遞矩陣；應(yīng)用

轉(zhuǎn)子支承系統(tǒng)是旋轉(zhuǎn)機(jī)械的主要部件，國(guó)內(nèi)外許多學(xué)者已對(duì)其進(jìn)行了大量研究[1-9]。傳遞矩陣法自50年代中期被應(yīng)用于轉(zhuǎn)子動(dòng)力學(xué)分析，至今仍是轉(zhuǎn)子動(dòng)力學(xué)的主要數(shù)值分析方法之一[10]。袁惠群等[11]在Lagrange體系下運(yùn)用傳遞矩陣方法對(duì)轉(zhuǎn)子動(dòng)力學(xué)特性進(jìn)行了深入分析。在Lagrange體系下用傳遞矩陣法對(duì)轉(zhuǎn)子系統(tǒng)進(jìn)行動(dòng)力學(xué)分析時(shí)，首先要對(duì)轉(zhuǎn)子系統(tǒng)進(jìn)行離散，簡(jiǎn)化為集總質(zhì)量模型，把轉(zhuǎn)子系統(tǒng)分為圓盤(pán)、無(wú)質(zhì)量軸段和支承等單元，對(duì)各個(gè)單元進(jìn)行受力分析，將各個(gè)單元兩端狀態(tài)向量之間的傳遞關(guān)系用傳遞矩陣來(lái)表示，根據(jù)連接條件求得轉(zhuǎn)子系統(tǒng)的總傳遞矩陣，最后根據(jù)邊界條件，求得渦動(dòng)頻率等模態(tài)參數(shù)。

對(duì)轉(zhuǎn)子的動(dòng)力學(xué)分析同樣可以在辛體系下進(jìn)行。鐘萬(wàn)勰院士[12]將對(duì)偶變量體系引入到陀螺系統(tǒng)，從一個(gè)全新的角度研究了轉(zhuǎn)子動(dòng)力學(xué)問(wèn)題，證明了辛空間中分析陀螺系統(tǒng)的優(yōu)越性。隋永楓等[13-14]在辛空間中運(yùn)用有限元等計(jì)算方法對(duì)轉(zhuǎn)子進(jìn)行了研究，更好地解決了一系列轉(zhuǎn)子動(dòng)力學(xué)問(wèn)題。本論文將多支承轉(zhuǎn)子系統(tǒng)導(dǎo)入辛體系，結(jié)合轉(zhuǎn)子系統(tǒng)自身的特點(diǎn)，給出了在辛空間中的傳遞矩陣方法，體現(xiàn)了保辛的性質(zhì)及計(jì)算方法的優(yōu)越性，對(duì)運(yùn)用傳遞矩陣法分析多支承轉(zhuǎn)子動(dòng)力學(xué)特性提供了一種新的思路。

1 理論基礎(chǔ)

1.1 連續(xù)梁處傳遞矩陣的推導(dǎo)

本文轉(zhuǎn)子系統(tǒng)模型的連接梁采用各向同性的鐵木辛柯梁[15]，需要計(jì)入其轉(zhuǎn)動(dòng)慣量和剪切變形。

圖1 單元梁示意圖Fig.1 Beam element diagram

(1)

變形能表達(dá)式為

(2)

則拉格朗日密度函數(shù)為

(3)

(4)

引入復(fù)變量，設(shè)R=h1+ih2，則式(4)可以寫(xiě)為

(5)

在時(shí)域中采用Fourier展開(kāi)法，即令R=q(x)e-iωt，可以得到

(6)

式(6)中空間部分的變分為

(7)

其中，

(8)

引入變量q的對(duì)偶變量

(9)

(10)

(11)

(12)

方程的解為

v=eHliv0

(13)

式中：li為此段梁的長(zhǎng)度；v0和v分別為梁左端和右端的狀態(tài)向量。其中，eHli為一個(gè)辛矩陣。

1.2 對(duì)支承及圓盤(pán)處連接條件的處理

對(duì)偶變量p的物理意義為

(14)

式中，Q跟M分別為梁的剪力和彎矩。因?yàn)樗芯康氖歉飨蛲缘霓D(zhuǎn)子系統(tǒng)，y方向與z方向具有相同的運(yùn)動(dòng)規(guī)律，所以取y方向進(jìn)行分析可知，其全狀態(tài)向量實(shí)質(zhì)是v=[yθy-QyMy]。

將轉(zhuǎn)子的彈性支承及圓盤(pán)看為一段梁左端的邊界條件進(jìn)行處理。取y方向進(jìn)行分析，如果左端為彈簧，彈簧的剛度為ky，耦合剛度為kyz，對(duì)于各向同性轉(zhuǎn)子系統(tǒng)，[yθy-QyMy]T=[Zθz-QzMz]T，則該段(第i段)左端的狀態(tài)向量為

(15)

設(shè)單位矩陣為

(16)

單位辛矩陣為

(17)

若左端為圓盤(pán)，圓盤(pán)的質(zhì)量為m，則該段(第i段)左端的狀態(tài)向量為

(18)

由辛矩陣的性質(zhì)知，辛矩陣的乘積仍為辛矩陣 ，所以該算法是保辛的。

1.3 固有頻率及臨界轉(zhuǎn)速的計(jì)算

設(shè)求得的總傳遞矩陣為

(19)

即

vR=TzvL

(20)

根據(jù)模型的邊界條件，可得

即可計(jì)算出固有頻率或臨界轉(zhuǎn)速。

2 模型的建立

各項(xiàng)同性轉(zhuǎn)子簡(jiǎn)化模型如圖2所示，該轉(zhuǎn)子系統(tǒng)由一根轉(zhuǎn)軸，三個(gè)支承，兩個(gè)圓盤(pán)組成。轉(zhuǎn)軸被支承及圓盤(pán)分為四段，每段的長(zhǎng)度為2 m， 轉(zhuǎn)軸的橫截面為圓形，轉(zhuǎn)軸的彈性模量為E=2×1011N/m，泊松比為0.3，密度為ρ=7 800 kg/m3，截面的慣性矩為J=4.952×10-4m4，三個(gè)支承處的等效支承剛度分別為k11=3.92×108N/m，k12=4.90×108N/m,k13=4.00×108N/m，兩個(gè)圓盤(pán)的質(zhì)量分別為m1=1.764×102kg,m2=2.067 8×

102kg, 轉(zhuǎn)動(dòng)慣量分別為jd1=10.78 kg·m2,jp1=20.58 kg·m2,jd2=29.498 kg·m2,jp2=58.212 kg·m2。

圖2 轉(zhuǎn)子模型圖Fig.2 Rotor model diagram

3 轉(zhuǎn)子系統(tǒng)計(jì)算分析

3.1 固有頻率及振型的計(jì)算

在轉(zhuǎn)子動(dòng)力學(xué)分析中，有限元方法是一種有效的計(jì)算方法，發(fā)展已十分完善，在工程中有著十分廣泛的應(yīng)用[16-17]。本文運(yùn)用MATLAB進(jìn)行編程計(jì)算，以有限元方法的計(jì)算結(jié)果為基準(zhǔn)，將本文計(jì)算結(jié)果與傳統(tǒng)方法的計(jì)算結(jié)果進(jìn)行比較。

由表1及圖 3～圖5 可知，本文方法計(jì)算的固有頻率與有限元方法的計(jì)算結(jié)果非常接近，誤差在可接受的范圍內(nèi)，具有較高的精確度，同時(shí)運(yùn)用本文方法得到的振型曲線與有限元方法得到的振型曲線相吻合。與傳統(tǒng)傳遞矩陣法相比，該方法在計(jì)算高階渦動(dòng)頻率時(shí)避免了傳統(tǒng)方法容易出現(xiàn)數(shù)值不穩(wěn)定等缺點(diǎn)，體現(xiàn)出了本文計(jì)算方法的優(yōu)越性，說(shuō)明本文的方法是一種行之有效的計(jì)算方法。

表1 前四階固有頻率及相對(duì)誤差

分別運(yùn)用有限元方法及本文方法所求得的前三階振型圖如圖3～圖5。

圖3 一階振型曲線Fig.3 First order mode curve

圖4 二階振型曲線Fig.4 Second order mode curve

圖5 三階振型曲線Fig.5 Third order mode curve

3.2 臨界轉(zhuǎn)速的計(jì)算

在實(shí)際工程中，通常只需考慮同向渦動(dòng)時(shí)的臨界轉(zhuǎn)速。為了繪制坎貝爾圖，需要求出多個(gè)自轉(zhuǎn)頻率所對(duì)應(yīng)的各階渦動(dòng)頻率。在本文中，自轉(zhuǎn)角速度在0～1 000 rad/s內(nèi)以50 rad/s為步長(zhǎng)取了21個(gè)點(diǎn)，圖6中的交點(diǎn)所對(duì)應(yīng)的轉(zhuǎn)速即為轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速。

圖6 坎貝爾圖Fig.6 Campbell diagram

計(jì)算得到各階渦動(dòng)頻率所對(duì)應(yīng)的臨界轉(zhuǎn)速，如表2所示。

表2 各階臨界轉(zhuǎn)速

表2的計(jì)算結(jié)果與圖6相對(duì)應(yīng)，也進(jìn)一步說(shuō)明了本文方法的有效性。由圖6及表2可看出，在本文模型中，自轉(zhuǎn)角速度的變化對(duì)渦動(dòng)頻率的影響不明顯，第一階臨界轉(zhuǎn)速和第二階臨界轉(zhuǎn)速較為接近，第三階臨界轉(zhuǎn)速和第四階臨界轉(zhuǎn)速較為接近，實(shí)際工程中的轉(zhuǎn)子在設(shè)計(jì)時(shí)，應(yīng)當(dāng)避免其工作時(shí)的轉(zhuǎn)速在臨界轉(zhuǎn)速附近，以防止振動(dòng)幅值過(guò)大。

3.3 頻響分析

自轉(zhuǎn)角速度在0～1 000 rad/s內(nèi)以40 rad/s為步長(zhǎng)點(diǎn)，算出了轉(zhuǎn)子左端點(diǎn)在每個(gè)自轉(zhuǎn)角速度下所對(duì)應(yīng)的頻響函數(shù)曲線，繪制出瀑布圖，如圖7所示。

圖7 瀑布圖Fig.7 Waterfall plot

由圖7可知，在本文模型中，轉(zhuǎn)子的自轉(zhuǎn)角速度對(duì)其渦動(dòng)頻率的影響不明顯。圖7所示的計(jì)算結(jié)果跟圖6所示的計(jì)算結(jié)果相吻合，同時(shí)也驗(yàn)證了本文求解方法的有效性。

3.4 不同因素對(duì)臨界轉(zhuǎn)速的影響

在實(shí)際工程中，通常都是通過(guò)調(diào)節(jié)支承的剛度來(lái)對(duì)轉(zhuǎn)子系統(tǒng)進(jìn)行優(yōu)化，本論文將討論支承剛度的改變對(duì)多支承轉(zhuǎn)子系統(tǒng)各階臨界轉(zhuǎn)速的影響的幾種情況。在討論端部支承剛度的改變時(shí)，選取轉(zhuǎn)子系統(tǒng)的左端為例進(jìn)行討論。在討論同時(shí)改變兩個(gè)支承剛度對(duì)各階臨界轉(zhuǎn)速的影響時(shí)，通過(guò)比較發(fā)現(xiàn)同時(shí)改變中間跟端部支承剛度對(duì)臨界轉(zhuǎn)速的影響最大，故本文只給出了此種情況下的結(jié)果，同時(shí)改變兩端支承剛度時(shí)可用類似方法進(jìn)行分析。

(1)分別改變左端及中間彈性支承的剛度對(duì)各階臨界轉(zhuǎn)速的影響。

由圖8～圖10可知，無(wú)論改變中間彈性支承還是端部彈性支承，隨著支承剛度的增加，各階臨界轉(zhuǎn)速都是先有一個(gè)快速上升的階段，然后趨于平穩(wěn)。支承剛度的改變對(duì)低階臨界轉(zhuǎn)速的影響較小，對(duì)高階臨界轉(zhuǎn)速影響較大。在本論文所討論的轉(zhuǎn)子系統(tǒng)中，通過(guò)對(duì)改變中間支承剛度跟端部支承剛度所得到的結(jié)果進(jìn)行對(duì)比可知，中間支承剛度對(duì)各階臨界轉(zhuǎn)速有較大的影響。

圖8 對(duì)一階臨界轉(zhuǎn)速的影響Fig.8 Influence on the first order critical speed

圖9 對(duì)二階臨界轉(zhuǎn)速的影響Fig.9 Influence on the second order critical speed

圖10 對(duì)三階臨界轉(zhuǎn)速的影響曲線Fig.10 Influence on the third order critical speed

(2)同時(shí)改變左端跟中間支承剛度對(duì)各階臨界轉(zhuǎn)速的影響。

由圖11～圖13可知，隨著兩個(gè)支承剛度的增加，各階臨界轉(zhuǎn)速都先有一個(gè)快速上升的階段，然后趨于平穩(wěn)。由圖11～圖13同樣可以看出本文算例中，同端部支承相比，中間支承剛度對(duì)各階臨界轉(zhuǎn)速的影響較大。

圖11 對(duì)一階臨界轉(zhuǎn)速的影響Fig.11 Influence on the first order critical speed

圖12 對(duì)二階臨界轉(zhuǎn)速的影響Fig.12 Influence on the second order critical speed

圖13 對(duì)三階臨界轉(zhuǎn)速的影響Fig.13 Influence on the third order critical speed

4 結(jié) 論

本文將多支承轉(zhuǎn)子的動(dòng)力學(xué)分析引入到辛體系中，運(yùn)用辛方法推導(dǎo)出了狀態(tài)向量之間的傳遞矩陣，該方法具有方便用MATLAB進(jìn)行編程，計(jì)算速度快的優(yōu)點(diǎn)。與傳統(tǒng)的傳遞矩陣法相比較，該方法避免了求解梁的聚集模型等一系列復(fù)雜操作，便于理解，在計(jì)算高階渦動(dòng)頻率時(shí)也可以得到比較精確的計(jì)算結(jié)果，體現(xiàn)出本文計(jì)算方法的有效性及優(yōu)越性。同時(shí)本論文對(duì)臨界轉(zhuǎn)速及支承剛度的改變對(duì)臨界轉(zhuǎn)速的影響進(jìn)行了分析，為轉(zhuǎn)子系統(tǒng)的優(yōu)化設(shè)計(jì)提供了依據(jù)。

[ 1 ] 孟光. 轉(zhuǎn)子動(dòng)力學(xué)研究的回顧與展望[J]. 振動(dòng)工程學(xué)報(bào), 2002, 15(1): 5-13. MENG Guang. Retrospect and prospect to the research on rotor dynamics[J]. Journal of Vibration Engineering, 2002, 15(1): 5-13.

[ 2 ] 楊永鋒, 任興民, 徐斌. 國(guó)外轉(zhuǎn)子動(dòng)力學(xué)研究綜述[J]. 機(jī)械科學(xué)與技術(shù), 2011, 30(10): 1775-1780. YANG Yongfeng, REN Xingmin, XU Bin. Review of international researches on rotor dynamics[J]. Mechanical Science and Technology for Aerospace Engineering, 2011, 30(10): 1775-1780.

[ 3 ] 徐小明, 鐘萬(wàn)勰. 轉(zhuǎn)子動(dòng)力學(xué)的非線性數(shù)值求解[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2015,36(7): 677-685. XU Xiaoming, ZHONG Wanxie. Nonlinear numerical simulation of rotor dynamics[J]. Applied Mathematics and Mechanics, 2015,36(7): 677-685.

[ 4 ] 李強(qiáng), 馬龍, 許偉偉, 等. 橢圓軸承-轉(zhuǎn)子耦合系統(tǒng)動(dòng)力學(xué)特性研究[J]. 振動(dòng)與沖擊, 2016,35(11): 174-179. LI Qiang, MA Long, XU Weiwei, et al. Dynamic characteristics analysis of an elliptical bearing-rotor coupled system[J]. Journal of Vibration and Shock, 2016,35(11): 174-179.

[ 5 ] KIM J S, LEE S H. The stability of active balancing control using influence coefficients for a variable rotor system[J]. International Journal of Advanced Manufacturing Technology, 2003, 22(7/8): 562-567.

[ 6 ] 張文. 轉(zhuǎn)子動(dòng)力學(xué)理論基礎(chǔ)[M]. 北京: 科學(xué)出版社, 1990.

[ 7 ] 何洪軍, 荊建平. 非線性轉(zhuǎn)子-密封系統(tǒng)動(dòng)力學(xué)模型研究[J]. 振動(dòng)與沖擊, 2014, 33(10): 73-76. HE Hongjun, JING Jianping. Dynamic model of a nonlinear rotor-seal system[J]. Journal of Vibration and Shock, 2014, 33(10): 73-76.

[ 8 ] KIM T, NA S. New automatic ball balancer design to reduce transient-response in rotor system[J]. Mechanical Systems and Signal Processing, 2013, 37(1/2): 265-275.

[ 9 ] BEHZAD M, ALVANDI M, MBA D, et al. A finite element-based algorithm for rubbing induced vibration prediction in rotors[J]. Journal of Sound and Vibration, 2013, 332(21): 5523-5542.

[10] 鐘一諤,何衍宗,王正,等. 轉(zhuǎn)子動(dòng)力學(xué)[M]. 北京: 清華大學(xué)出版社, 1987.

[11] 袁惠群. 復(fù)雜轉(zhuǎn)子系統(tǒng)的矩陣分析方法[M]. 沈陽(yáng): 遼寧科學(xué)技術(shù)出版社, 2014.

[12] 鐘萬(wàn)勰. 經(jīng)典力學(xué)辛講[M]. 大連: 大連理工大學(xué)出版社, 2013.

[13] 隋永楓. 轉(zhuǎn)子動(dòng)力學(xué)的求解辛體系及其數(shù)值計(jì)算方法[D]. 大連： 大連理工大學(xué), 2006.

[14] 隋永楓, 高強(qiáng), 鐘萬(wàn)勰. 陀螺系統(tǒng)時(shí)間有限元方法[J]. 振動(dòng)與沖擊, 2012, 31(13): 95-98. SUI Yongfeng, GAO Qiang, ZHONG Wanxie. Time domain finite element method for gyroscopic systems[J]. Journal of Vibration and Shock, 2012, 31(13): 95-98.

[15] TIMOSHENKO S. 材料力學(xué)[M]. 北京: 科學(xué)出版社, 1978.

[16] 陳錫棟, 楊婕, 趙曉棟, 等. 有限元法的發(fā)展現(xiàn)狀及應(yīng)用[J]. 中國(guó)制造業(yè)信息化, 2010, 39(11): 6-8. CHEN Xidong, YANG Jie, ZHAO Xiaodong, et al. The status and development of finite element method[J]. Machine Design and Manufacturing Engineering, 2010, 39(11): 6-8.

[17] 孫海霞, 戴京濤, 姜偉, 等. 有限元法在機(jī)械工程中的應(yīng)用與發(fā)展[J]. 科技創(chuàng)新導(dǎo)報(bào), 2011(3): 84. SUN Haixia, DAI Jingtao, JIANG Wei, et al. Application and development of finite element method in mechanical engineering[J]. Science and Technology Innovation Herald, 2011(3): 84.

A transfer matrix method for a multi-support rotor systemin the symplectic space and its application

ZHANG Juanjuan1, CUI Sheng1, FENG Yongxin2

(1. Department of Aeronautics and Astronautics，F(xiàn)udan University, Shanghai 200433, China; 2. Guangdong Power Test & Research Institute, Guangzhou 510600, China )

Based on the variational principle, dual variables were introduced, and the multi-support rotor system was imported into symplectic space. The symplectic transfer matrix in the symplectic space was presented for the rotor system, which reflects the symplectic conservation property. A model of the rotor system was built. MATLAB was used to analyze the model. The natural frequency and other modal parameters were obtained. The calculated values are consistent with the values obtained by the finite element method. The effectiveness of the method was confirmed by comparing the results with traditional methods. After large amount of calculations the critical speeds were obtained. The relationship between the critical speeds and the support stiffness was studied. This work provides reference for the optimization of rotor systems.

rotor system; symplectic space; transfer matrix; application

廣東電網(wǎng)公司電力科學(xué)研究院科研基金

2016-10-12 修改稿收到日期： 2016-11-23

張娟娟 女，碩士，1990年1月生

崔升 男，博士，副教授，1967年9月生

O347.6

A

10.13465/j.cnki.jvs.2017.16.005