馬 超，湯華濤，賀培文

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

0 前言

燃?xì)廨啓C(jī)是艦船、航空動(dòng)力和工業(yè)生產(chǎn)中常用的一種發(fā)動(dòng)機(jī)，而燃?xì)廨啓C(jī)的轉(zhuǎn)子是燃機(jī)的重要組成部分，國內(nèi)外學(xué)者對其展開了廣泛的研究。目前，學(xué)者們采用的建模方法大體可以歸類為多體系統(tǒng)動(dòng)力學(xué)、傳遞矩陣法和有限單元法3種。但隨著計(jì)算量的加大和計(jì)算精度的不斷提高，這3種方法的缺點(diǎn)日益突出，為了解決這些問題，必須尋找一種新的適合燃?xì)廨啓C(jī)轉(zhuǎn)子的計(jì)算方法。

20世紀(jì)末，芮筱亭等將多體系統(tǒng)動(dòng)力學(xué)和傳遞矩陣法相結(jié)合，成功解決了復(fù)雜多體發(fā)射系統(tǒng)特征值問題，并首次提出多體系統(tǒng)傳遞矩陣法［1］。經(jīng)過10多年的不斷完善和工程應(yīng)用以及廣泛的國際學(xué)術(shù)交流，最終形成了較為系統(tǒng)的多體系統(tǒng)傳遞矩陣法。它建模靈活，簡潔，程式化程度高，計(jì)算量小，已得到廣泛應(yīng)用。由于該方法的計(jì)算特點(diǎn)，它特別適合解決鏈?zhǔn)浇Y(jié)構(gòu)系統(tǒng)問題，如連續(xù)梁，渦輪軸等。因此，本文將選擇多體系統(tǒng)傳遞矩陣法建立某燃?xì)廨啓C(jī)高壓渦輪壓氣機(jī)轉(zhuǎn)子計(jì)算模型。

由于該轉(zhuǎn)子的工作特點(diǎn)，其結(jié)構(gòu)比較復(fù)雜，轉(zhuǎn)子內(nèi)部有多處凸起，如果采用傳統(tǒng)的偏微分方程來推導(dǎo)其傳遞矩陣，難度太大。本文以某燃機(jī)高壓渦輪壓氣機(jī)轉(zhuǎn)子為例，利用有限元法推導(dǎo)轉(zhuǎn)子各段的剛度矩陣，再通過矩陣變換求出其傳遞矩陣，最后將各段傳遞矩陣組合為整體傳遞矩陣，建立高壓轉(zhuǎn)子的多體系統(tǒng)傳遞矩陣模型。

1 凸起部分傳遞矩陣

圖1為某燃機(jī)高壓轉(zhuǎn)子結(jié)構(gòu)示意圖，轉(zhuǎn)子有多處凸起部分，應(yīng)該先單獨(dú)考慮凸起部分，分別推導(dǎo)其傳遞矩陣，然后將各段傳遞矩陣組合為整體傳遞矩陣。

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

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

圖1 某燃?xì)廨啓C(jī)高壓渦輪壓氣機(jī)轉(zhuǎn)子Fig.1 High-pressure compressor rotor of gas turbine

圖2 凸起部分結(jié)構(gòu)圖Fig.2 Structure of rotor with convexity

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

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

假設(shè)此凸起的節(jié)點(diǎn)1不受力，即 F1=M1=0，則

矩陣A即整個(gè)軸結(jié)構(gòu)中節(jié)點(diǎn)2到節(jié)點(diǎn)3的剛度矩陣，剛度矩陣A考慮了梁單元1對整個(gè)軸結(jié)構(gòu)的影響，且A為四階矩陣，可將A寫成如下形式:

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

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

2 轉(zhuǎn)子整體傳遞矩陣模型

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

圖4 轉(zhuǎn)子多體系統(tǒng)傳遞矩陣模型Fig.4 Transfer matrix method of multibody system model of rotor

定義狀態(tài)矢量Z0，1， Z2，1， Z2，3， …， Z44，45，Z46，45的形式均為 Z=［u，θ，M，F(xiàn)］T，定義各段傳遞矩陣U1，U2，U3，…，U45，U46，其中轉(zhuǎn)子各段的傳遞矩陣已得到。根據(jù)多體系統(tǒng)傳遞矩陣法的知識，集中質(zhì)量點(diǎn)的傳遞矩陣為［5］

式中:m為集中質(zhì)量點(diǎn)的質(zhì)量;ω為系統(tǒng)固有頻率。最終可得到轉(zhuǎn)子的總傳遞方程為

已知U為四階矩陣，設(shè)其表達(dá)式為

在此模型中，轉(zhuǎn)子兩端固定約束，故有:

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

3 算例分析

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

表1 計(jì)算結(jié)果對比Tab.1 Comparison of results

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

本文方法還繼承了Matlab自編程序靈活多變的特點(diǎn)，可以改進(jìn)程序，控制誤差，提高計(jì)算精度，還可以優(yōu)化算法以減少計(jì)算時(shí)間。

4 結(jié)語

本文以某燃?xì)廨啓C(jī)高壓渦輪壓氣機(jī)轉(zhuǎn)子為例，利用有限元法推導(dǎo)轉(zhuǎn)子各段的剛度矩陣，再通過矩陣變換求出其傳遞矩陣，最后將各段傳遞矩陣組合為整體傳遞矩陣，建立高壓渦輪壓氣機(jī)轉(zhuǎn)子的多體系統(tǒng)傳遞矩陣模型。通過算例證明，該方法計(jì)算準(zhǔn)確，易于實(shí)現(xiàn)，為用有限元法和多體系統(tǒng)傳遞矩陣法求解燃?xì)廨啓C(jī)轉(zhuǎn)子動(dòng)力學(xué)問題提供了新思路。

［1］芮筱亭，贠來峰，陸毓琪，等.多體系統(tǒng)傳遞矩陣法及其應(yīng)用［M］.北京:科學(xué)出版社，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］李開禧，趙廣坡，熊曉莉.薄壁梁的單元?jiǎng)偠染仃嚰捌鋺?yīng)用［J］.重慶建筑大學(xué)學(xué)報(bào)，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］王勖成，邵敏.有限單元法基本原理和數(shù)值方法［M］.北京:清華大學(xué)出版社，1997.

［5］芮筱亭，劉亞飛.梁系統(tǒng)振動(dòng)的傳遞矩陣法［J］.力學(xué)與實(shí)踐，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］關(guān)琦，金鶴，新力.某型燃?xì)廨啓C(jī)低壓渦輪壓氣機(jī)轉(zhuǎn)子動(dòng)力學(xué)分析［J］.艦船科學(xué)技術(shù)，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.