某艦用燃機高壓渦輪壓氣機轉子傳遞矩陣建模

馬 超，湯華濤，賀培文

(1.中國人民解放軍91343部隊，山東 威海 264200;2.海軍工程大學機械工程系，湖北 武漢 430033;3.中國人民解放軍91329部隊，山東 威海 264200)

0 前言

燃氣輪機是艦船、航空動力和工業生產中常用的一種發動機，而燃氣輪機的轉子是燃機的重要組成部分，國內外學者對其展開了廣泛的研究。目前，學者們采用的建模方法大體可以歸類為多體系統動力學、傳遞矩陣法和有限單元法3種。但隨著計算量的加大和計算精度的不斷提高，這3種方法的缺點日益突出，為了解決這些問題，必須尋找一種新的適合燃氣輪機轉子的計算方法。

20世紀末，芮筱亭等將多體系統動力學和傳遞矩陣法相結合，成功解決了復雜多體發射系統特征值問題，并首次提出多體系統傳遞矩陣法［1］。經過10多年的不斷完善和工程應用以及廣泛的國際學術交流，最終形成了較為系統的多體系統傳遞矩陣法。它建模靈活，簡潔，程式化程度高，計算量小，已得到廣泛應用。由于該方法的計算特點，它特別適合解決鏈式結構系統問題，如連續梁，渦輪軸等。因此，本文將選擇多體系統傳遞矩陣法建立某燃氣輪機高壓渦輪壓氣機轉子計算模型。

由于該轉子的工作特點，其結構比較復雜，轉子內部有多處凸起，如果采用傳統的偏微分方程來推導其傳遞矩陣，難度太大。本文以某燃機高壓渦輪壓氣機轉子為例，利用有限元法推導轉子各段的剛度矩陣，再通過矩陣變換求出其傳遞矩陣，最后將各段傳遞矩陣組合為整體傳遞矩陣，建立高壓轉子的多體系統傳遞矩陣模型。

1 凸起部分傳遞矩陣

圖1為某燃機高壓轉子結構示意圖，轉子有多處凸起部分，應該先單獨考慮凸起部分，分別推導其傳遞矩陣，然后將各段傳遞矩陣組合為整體傳遞矩陣。

可得轉子的凸起部分簡化為圖2中結構，并以凸起中線為軸，將凸起分為左右2個部分。單獨考慮凸起左半部分，將其劃分為2個梁單元，編號為Ⅰ、Ⅱ，如圖3所示。若只考慮梁單元的橫向振動，可假設2個梁單元的剛度矩陣分別為KⅠ和KⅡ，2個剛度矩陣都是四階矩陣［2－4］，其表達式如下:

將2個梁單元的剛度矩陣集成為整體剛度矩陣［5］，則可得到整體剛度方程

圖1 某燃氣輪機高壓渦輪壓氣機轉子Fig.1 High-pressure compressor rotor of gas turbine

圖2 凸起部分結構圖Fig.2 Structure of rotor with convexity

圖3 凸起梁單元模型Fig.3 Model of beam element with convexity

式中:Fi和Mi(i=1，2，3)為節點i所受的力和力矩;ui和θi(i=1，2，3)為節點i的位移和轉角。

假設此凸起的節點1不受力，即 F1=M1=0，則

矩陣A即整個軸結構中節點2到節點3的剛度矩陣，剛度矩陣A考慮了梁單元1對整個軸結構的影響，且A為四階矩陣，可將A寫成如下形式:

定義狀態矢量 Z2={u2θ2M2F2}T，Z3={u3θ3M3F3}T，將式(7)展開，用 Z2表示Z3，則其傳遞方程［1］為

用同樣的方法，可以推導出凸起右半部分的傳遞矩陣。

2 轉子整體傳遞矩陣模型

轉子無凸起部分的傳遞矩陣文獻 ［1］中已給出，且前文已經利用有限元法推導出轉子凸起部分的傳遞矩陣，但由于在推導過程中使用的是無質量梁單元模型，所得到的傳遞矩陣也沒有質量信息，無法計算其振動，因此在轉子整體傳遞矩陣模型中必須考慮其質量。本文的處理方法是在轉子各段中點處添加集中質量點，其質量為各段總質量。將圖1中的轉子模型分段，并依次編號為0～46，如圖4所示，其中在中心線以上的編號表示集中質量點。

圖4 轉子多體系統傳遞矩陣模型Fig.4 Transfer matrix method of multibody system model of rotor

定義狀態矢量Z0，1， Z2，1， Z2，3， …， Z44，45，Z46，45的形式均為 Z=［u，θ，M，F］T，定義各段傳遞矩陣U1，U2，U3，…，U45，U46，其中轉子各段的傳遞矩陣已得到。根據多體系統傳遞矩陣法的知識，集中質量點的傳遞矩陣為［5］

式中:m為集中質量點的質量;ω為系統固有頻率。最終可得到轉子的總傳遞方程為

已知U為四階矩陣，設其表達式為

在此模型中，轉子兩端固定約束，故有:

求解式(14)，即可得到系統的固有頻率［6］。

3 算例分析

以圖1中轉子為計算模型，用多體系統傳遞矩陣法計算其橫向振動固有頻率。為了驗證該方法的正確性，可與有限元軟件Ansys計算結果做比較［7］。在有限元前處理時，為了保證計算精度，將轉子的幾何模型導入有限元軟件HyperMesh中，采用六面體劃分網格，單元類型選用Solid185，最后得到的轉子有限元模型共有162214個單元，206196個節點。然后將有限元模型導入Ansys中，在轉子的兩端施加位移約束，計算其約束模態。本文只考慮轉子的橫向振動，在Ansys中計算轉子的前20階模態，并取出其振型為橫向振動的模態，與本文方法計算結果相比較 (見表1)。

表1 計算結果對比Tab.1 Comparison of results

從表1可看出，用2種方法計算得到的轉子前兩階模態頻率差距很小，驗證了本文方法的正確性。用Ansys軟件計算時，計算過程耗時3 h，加上前處理時間，共用時10 h。本文方法從編程到計算共用時5 h，其缺點是編程和程序調試過程較繁瑣，但是當方法成熟之后計算轉子的其他問題就比較方便，所以本文方法比Ansys軟件計算簡潔方便。

本文方法還繼承了Matlab自編程序靈活多變的特點，可以改進程序，控制誤差，提高計算精度，還可以優化算法以減少計算時間。

4 結語

本文以某燃氣輪機高壓渦輪壓氣機轉子為例，利用有限元法推導轉子各段的剛度矩陣，再通過矩陣變換求出其傳遞矩陣，最后將各段傳遞矩陣組合為整體傳遞矩陣，建立高壓渦輪壓氣機轉子的多體系統傳遞矩陣模型。通過算例證明，該方法計算準確，易于實現，為用有限元法和多體系統傳遞矩陣法求解燃氣輪機轉子動力學問題提供了新思路。

［1］芮筱亭，贠來峰，陸毓琪，等.多體系統傳遞矩陣法及其應用［M］.北京:科學出版社，2008.

［2］HOWSON R.Exact dynamic stiffness matrix for a thinwalled beam-column of doubly asymmetric cross-section［J］.Structural Engineering and Mechanics，2011，38(2):195－210.

［3］李開禧，趙廣坡，熊曉莉.薄壁梁的單元剛度矩陣及其應用［J］.重慶建筑大學學報，2004，26(5):35 －38.LI Kai-xi，ZHAO Guang-po，XIONG Xiao-li.Element stiffness matrix of thin-walled beam and its application［J］.Journal of Chongqing Jianzhu University，2004，26(5):35 －38.

［4］王勖成，邵敏.有限單元法基本原理和數值方法［M］.北京:清華大學出版社，1997.

［5］芮筱亭，劉亞飛.梁系統振動的傳遞矩陣法［J］.力學與實踐，1993(5):66－67.RUI Xiao-ting，LIU Ya-fei.Transfer matrix method on vibraton of beam［J］.Mechanics and Engineering，1993(5):66 －67.

［6］LOEWY R G，BHNTANI N.Combined finite elementtransfer matrix method［J］.Journal of Sound and Vibration，1999，226(5):1048 －1052.

［7］關琦，金鶴，新力.某型燃氣輪機低壓渦輪壓氣機轉子動力學分析［J］.艦船科學技術，2010，32(8):127 －132.GUAN Qi，JIN He，XIN Li.Analysis on rotor dynamic of the low tuibocompressor［J］.Ship Science and Technology，2010，32(8):127－132.