航空發動機轉子系統建模方法和振動特性分析

張 力 洪 杰 馬艷紅

(中航工業沈陽發動機設計研究所，沈陽110015)(北京航空航天大學 能源與動力工程學院，北京100191)

對轉子系統進行動力學特性分析時，主要有傳遞矩陣法和有限元法.隨著計算機技術和有限元商用軟件的發展，有限元方法以其便捷高效的特點廣泛應用于復雜轉子系統的動力學計算.近些年，大量的轉子系統動力學研究基于有限元方法，建模時的單元選取主要分為梁單元和實體單元兩類.文獻［1－2］早年采用梁單元對轉子系統建模分析，考慮了轉動慣量和彎曲效應.文獻［3］的梁單元模型中，可以考慮轉子的陀螺力矩效應.之后，文獻［4－5］分別將瑞利梁理論和鐵木辛柯梁理論應用到臨界轉速分析中.文獻［6］采用梁單元對發動機整機建模并縮減自由度.文獻［7］采用實體單元對錐形軸段建模，較梁單元有更好的準確性.文獻［8］建立了轉子的實體單元循環對稱模型，同樣表明在計算不平衡響應時比梁單元具有提升.文獻［9］采用ANSYS軟件對轉子系統進行建模和動力學分析，證明程序可以考慮轉子系統的陀螺力矩效應.文獻［10］采用基于實體單元的整機動力特性分析方法，對復雜轉子系統建模時，考慮了陀螺力矩的影響.文獻［11］提出了兩種模化長葉片方法，對汽輪機轉子進行建模并求解臨界轉速，與試驗結果對比分析表明，該方法可以準確地模擬輪盤和長葉片的陀螺力矩效應.渦扇發動機轉子系統與地面旋轉機械不同，由于對重量的控制，轉子系統的剛度一般較弱，并含有大尺寸、大質量、大轉動慣量的葉盤結構，需要充分考慮其結構特征和動力學特征的影響.在建模方法上，主要問題集中在選用的單元和模型結構的簡化形式，將直接影響結果的準確性和計算經濟性.

本文對航空發動機典型轉子系統進行結構特征分析，對比不同單元類型和不同模型簡化方法對求解轉子系統動力學特性的影響，提出一種等效剛性圓環法對具有柔性轉子特征的系統進行有限元建模分析.

1 單元類型對轉子動力特性的影響

以某簡單盤軸轉子為例，分析離心預應力效應、陀螺力矩效應以及旋轉軟化效應對共振轉速的影響.模型為簡單的單盤轉子，如圖1所示.相關參數:L=305 mm，a=500 mm，d1=30 mm，d2=300 mm，B=25 mm，彈性模量 E=2.06×105MPa，泊松比 μ=0.3，密度 ρ=7 800 kg/m3，左支點為剛性固定鉸支，右支點為剛性可動鉸支.

圖1 單盤轉子模型示意圖

分別采用梁單元和實體單元對0轉速和500 rad/s轉速進行共振轉速求解，求得前三階共振轉速結果對比如表1、表2和表3所示，相對應的模態振型如圖2所示.

表1 共振轉速計算結果對比表(0轉速)

由表1～表3的計算結果對比，可得結論如下:

對于簡單盤軸轉子，在不考慮轉速影響時，采用梁單元和實體單元計算得到的各階固有頻率基本一致，相對誤差很小;轉速對轉子系統共振轉速的影響主要表現在離心載荷效應、陀螺力矩效應以及旋轉軟化效應上，其中離心載荷效應和陀螺力矩效應對轉子系統各階共振轉速影響較大，相應的影響程度和所對應振型有關，而旋轉軟化效應對轉子系統橫向振動的影響很小，因此可以忽略;梁單元模型由于求解程序的限制，不能考慮離心預應力效應和旋轉軟化效應的影響，因此在某些階共振轉速的求解中可能帶來較大誤差.

2 轉子動力特性相似與有限元建模

對于復雜的航空發動機轉子系統進行有限元模型時，必須對其結構進行必要地簡化，如何簡化主要取決于結構對轉子系統動力學特征的影響.主要表現在對軸的彎曲剛性，質量分布和轉動慣量的陀螺力矩效應等幾個方面.對于具體結構的的簡化主要分為:對盤軸的簡化，可采用三維實體單元或是一維梁單元模型;對葉片的簡化，可采用三維實體單元和質量單元模型兩種.

2.1 航空發動機典型柔性轉子系統

圖3所示為典型高涵道比渦扇發動機高低壓轉子結構，其特征是低壓轉子穿過高壓轉子，軸徑細且跨度長，風扇多采用變截面彎扭復合葉片，通常有3～5級渦輪，葉片長且數目多，因而葉片質量、慣性力矩和離心力均會對轉子系統的動力特性產生復雜影響，錐殼結構較多.其一方面，轉子系統結構復雜，尤其是變截面彎扭葉片的大量應用，完全按照幾何特征建模，耗時費力，計算經濟性差;另一方面，考慮到質量和轉動慣量分布、離心力和陀螺力矩效應等均會對轉子系統的動力學特性產生重要影響，需要通過定量分析對轉子動力特性的具體影響，來選取這些力學特征參數和建立相應的模型.

圖3 典型高涵道比渦扇發動機低壓轉子

2.2 盤軸簡化建模

對于轉子系統主要組成結構盤－軸的有限元建模，通常采用有兩種方法:一是梁單元模型，忽略盤的幾何特征，使用質量單元將等效質量加入模型中;另一種是采用實體單元建立模型，這種建模方法可以最大限度的保留其幾何和力學特征.下面以一個大涵道比渦扇發動機低壓轉子為例，通過轉子系統的動力學特性計算結果對比，分析兩種建模方法的差異.

在發動機結構方案設計階段，通常采用梁－剛性輪盤-彈簧模型進行模擬，這種模型計算效率高，適用于有大量的設計參數需要調整的初步階段.將低壓轉子沿軸向劃分為150個梁單元，梁單元模型分段如圖4所示.在軸的梁單元建模中，通過梁單元直徑的不同可以調整軸彎曲剛度，滿足與實際軸的力學等效.對于風扇、增壓級和渦輪葉片和輪盤，則只考慮其質量作為集中質量單元(考慮質量和轉動慣量)和外力作用在其質心所在的軸向位置上.軸承支點位置建立彈簧單元，用以等效模擬支承剛度.

圖4 低壓轉子梁單元振型

對于實體單元建模，采用8節點6面體的實體單元對盤、軸和鼓筒等部件進行單元離散，建立實體單元模型如圖5所示.對于葉片的力學特征，則通過在相應每級輪轂外緣上沿周向均布多個集中質量單元，進行力學特征模擬和等效.在3個支承位置建立彈簧單元，模擬支承剛度對轉子系統的動力學影響.

圖5 低壓轉子實體單元模型

采用ANSYS程序，分別對0 r/min，1 000 r/min，3000 r/min和6000 r/min 4個轉速下的振動固有特性進行求解，計算中考慮陀螺力矩的影響，各階模態振型圖如圖6所示，分別求得兩種模型前4階共振轉速如表4和表5所示.

圖6 轉子各階共振轉速下振型(三維實體單元模型)

依據表 4、表 5數據，繪制兩種模型的Campbell圖如7所示.

由圖7可見，分別采用梁單元和實體單元模型，在零轉速時各階固有頻率相對誤差很小.但隨轉速的提高，有限元模型對各階動頻曲線的影響規律產生了很大的不同.①在零轉速時各階固有頻率相對誤差很小;②對于轉子系統的第1，2階模態動頻曲線，兩種模型的差異較小，但實體單元模型隨轉速的增長速度略高于梁單元模型;③對于第3，4階模態動頻曲線，2種模型的差異較大，梁單元模型隨轉速的增長速度遠遠高于實體單元模型.

表4 轉子系統前4階共振轉速及振型(一維梁單元模型)

表5 轉子系統前4階共振轉速及振型(實體單元模型)

圖7 兩種模型的Campbell圖對比

分析兩種模型其原因如下:

對于第1，2階模態振型，由于其振型是整體一彎和渦輪平動，陀螺力矩的影響較小，建模方法上的差異影響不大，但是在梁單元模型中，不能較好地模化錐殼(進氣錐、渦輪盤根部錐殼段)和直徑突變帶來的軸段彎曲剛度變化，因而存在一些差異.

第3，4階模態振型分別為壓氣機和渦輪的俯仰振動.陀螺力矩對軸的附加剛度影響很大，在梁單元模型中，由于對盤和鼓筒進行局部剛化，沒有考慮輪盤變形的影響，這相當于增強了陀螺力矩效應，在俯仰振型的橫向振動中體現得較為明顯.

2.3 葉片簡化

大涵道比渦扇發動機低壓轉子的另外一個典型結構特征是帶有大尺寸的風扇葉片，大風扇葉片的尺寸和質量效應及盤片耦合振動將對轉子系統的橫向振動帶來影響，如果對各級葉片不做簡化，采用精確三維實體模型固然可以獲得轉子系統橫向振動模態的準確解，但計算量過大，同時由于盤片耦合效應，使得計算結果的后處理工作過于繁重，因此本節針對簡單的葉-盤-軸轉子模型，研究帶大質量葉片的轉子系統實體建模方法.

圖8為葉片-盤的有限元網格圖.在圖中各節點的位移矢量可表示為

其中

質量矩陣和剛度矩陣為

圖8 葉片-盤有限元網格

如果將葉片簡化為集中質量加在輪緣上，這時其對應的葉片位移δk則被忽略.

質量矩陣和剛度矩陣則變為

式(7)、式(8)兩個矩陣的變化將對轉子系統的固有振動特性產生一定的影響.

此外，另一個重要的力學特征，葉盤的耦合振動被忽略，這將使帶大尺寸風扇盤結構的轉子系統的橫向振動特性分析產生大的誤差.

建立葉片-盤-軸轉子模型如圖9所示，用于分析結構間的耦合振動對轉子的橫向振動特性的影響.

圖9 葉盤軸轉子模型

葉-盤-軸轉子系統的前4階模態振型如圖10所示.

圖10 葉-盤-軸系統的前4階振型

在不考慮轉速影響時，求得精確有限元模型以及各種簡化模型的共振轉速計算結果見表6.

表6 各種簡化模型固有頻率對比表

為便于分析，表中以整體有限元模型為基準值，各階固有頻率取為1，各種簡化模型的固有頻率值為相對比值.

由表6可見，對于第3階彎曲模態，由于其振型為純粹的軸段彎曲，葉盤幾乎沒有位移，因此，對于該階振動模態，各種葉盤等效方法對計算結果影響不大;對于轉子系統的1階、2階彎曲模態，葉片質量和轉動慣量的分布對固有頻率的影響非常明顯，因此忽略葉片或者將葉片作為質量單元的等效方法將會帶來一定的誤差.采用第3種簡化方法，將葉片作為剛體處理，并采用實體模型建模，可以準確的考慮葉片質量和轉動慣量，又可以忽略葉片和盤局部耦合振動帶來的復雜模態特征.

3 三維實體模型的動力學修正

由上述對比計算分析，可看出對于傳統的有限元建模方法，在帶有大尺寸葉盤結構的轉子系統中，進行轉子動力學計算分析時需要考慮:離心載荷效應的剛度增強，輪盤的彈性變形，葉-盤-軸耦合振動效應對轉子橫向振動的影響，錐形結構的模化.同時，還要力求控制有限元模型的規模，計算結果便于分析.因此，為準確等效大尺寸葉片的力學特征，大尺寸葉片等效為剛性環盤的方法對轉子系統進行有限元建模.

3.1 等效剛性環盤建模

既要準確反映大涵道比發動機低壓轉子的結構特征，又要盡可能減小模型規模和便于分析，應當在動力學上修正原有模型，使得在模型中可以用較少的單元來描述其質量和陀螺力矩的影響.可以采用等效剛性環盤來模擬葉片的質量分布和轉動慣量，軸承部分仍然使用相互垂直的一對彈簧單元代替.這樣整個轉子和葉片的離心效應、陀螺力矩以及旋轉軟化效應都能夠考慮進去，而且振動耦合也會更加真實，同時還可以濃縮模型的單元數，保證計算高效率.

對葉片的等效剛性環盤的建立需要滿足下列模化原則:①模化前后質量m不變;②模化前后極轉動慣量Jρ不變;③模化前后葉片離心力F不變;④環盤填充了原型葉片間的間隙，增強了周向的剛性，應減弱模型結構變化對結果的影響;⑤對轉子橫向振動的結果分析，避免由于葉片和葉盤的耦合振動增加分析的難度.

應用上述等效原則，對圖5所示轉子的實體模型進行動力修正.等效環盤需要確定的參數有外徑 R、內徑 r、密度 ρ、厚度 h、彈性模量 E.等效盤的內徑r即為葉根處的半徑，為已知.R，ρ，h可通過方程式(9)～式(11)聯立解出，原型參數可通過軟件從三維模型提取.

此外，建立環盤有限元模型時，可選用各項異性材料，通過提高徑向彈性模量使之產生剛化效應，減低周向彈性模量使之減少周向的約束，更加準確地模擬大尺寸葉片的力學特征.建立有限元模型如圖11所示.

圖11 低壓轉子動力學修正模型

3.2 計算結果

求得動力學修正實體模型在各轉速下的共振轉速如表7所示.

表7 不同轉速下轉子系統共振轉速

將修正實體單元的計算結果與前節梁單元與傳統實體單元的計算結果進行對比，各動頻曲線對比如圖12所示.

圖12 帶大風扇柔性轉子系統振動特性對比分析

3.3 結果分析

由圖12可見，前兩階振型渦輪盤和風扇盤的俯仰較小，陀螺力矩對共振轉速的影響較小，相對誤差較小.第3，4階振型頻率的相對誤差較大，由于這兩階振型全都帶有較大的葉盤俯仰擺動，葉盤陀螺力矩在這種振型下的作用體現的明顯.在傳統實體單元建模方法中，葉片等效為集中質量單元，在采用ANSYS程序求解時無法考慮非軸對稱模型的陀螺力矩效應，因而弱化了陀螺力矩的影響.結果表明，采用等效環盤法對非連續的風扇葉片結構進行簡化和等效建模，可以修正傳統實體有限元模型中陀螺力矩效應對固有頻率的影響，模型可以較好地用于后續的轉子動力學計算分析.

4 結論

通過對轉子系統建模方法的研究，指出了不同單元類型和不同簡化方法在轉子動力學分析中的影響因素.①梁單元模型不能考慮離心載荷效應對結構剛度增強作用，而對于某些橫向彎曲振型，離心載荷效應影響較大，因此可能帶來較大誤差;②采用梁單元和集中質量單元建模時，忽略了輪盤的柔性變形對大尺寸葉盤結構的陀螺力矩的影響，使陀螺效應增強;③梁單元模型忽略了葉片-盤-軸耦合振動效應對轉子系統橫向振動的影響，從而可能導致梁單元模型的陀螺力矩影響偏大;④對于錐、殼結構，采用梁單元進行有限元離散時，將引入結構簡化的誤差;⑤對葉片采用實體單元離散，由于模型為非軸對稱，在ANSYS程序中無法考慮其陀螺力矩效應，同時由于葉片-盤局部耦合振動，使得篩選轉子橫向振動的后處理工作過于繁重，很難獲得轉子橫向振動的各階模態結果.

在綜合考慮各種有限元模型特點的基礎上，提出了適用于帶有大尺寸葉片-盤結構轉子系統動力學分析的修正實體模型建立方法，使用等效剛性環盤來模擬葉片質量和轉動慣量分布及其對轉子動力特性的計算影響.通過對某型大涵道比渦扇發動機低壓轉子系統共振轉速的計算結果對比，驗證了葉片力學特征等效剛度環盤可以有效地模擬葉片的陀螺力矩影響，模型規模易于控制，可用于葉片-盤-軸轉子系統的橫向振動特性分析.

References)

［1］Ruhl R L.Dynamics of distributed parameter rotor systems:transfer matrix and finite element techniques［D］.New York:Cornell University，1970

［2］Ruhl R L，Booker J F.A finite element model for distributed parameter turborotor systems［J］.Journal of Engineering for Industry，1972，94:128－132

［3］Thorkildsen T.Solution of a distributed mass and unbalanced rotor system using a consistent mass matrix approach［D］.Arizona:Arizona State University，1974

［4］Polk S R.Finite element formulation and solution of flexible rotor-rigid disc systems for natural frequencies and critical whirl speed［D］.Arizona:Arizona State University，1974

［5］Nelson H D，Mcvaugh J M.The dynamics of rotor bearing systems using finite elements［J］.Journal of Engineering for Industry，1976，98:593－600

［6］James T，Charles L，Jonathan L.Reduced-order modeling and wavelet analysis of turbofan engine structural response due to foreign object damage(FOD)events［J］.Journal of Engineering for Gas Turbines and Power，2007，129(3):814 – 826

［7］Nandi A，Neogy S.Modelling of rotors with three-dimensional solid finite elements［J］.Journal of Strain Analysis for Engineering Design，2001，36(4):359 – 371

［8］Lazarus A，Prabel B，Combescure D.A 3D finite element model for the vibration analysis of asymmetric rotating machines［J］.Journal of Sound and Vibration，2010，329(18):3780－3797

［9］Kumar S M.Rotor dynamic analysis using ANSYS［C］//Gupta K.IUTAM Symposium on Emerging Trends in Rotor Dynamics.New York:Springer Science Business Media，2011:153－162

［10］陳萌，馬艷紅，劉書國，等.航空發動機整機有限元模型轉子動力學分析［J］.北京航空航天大學學報，2007，33(9):1014–1016 Chen Meng，Ma Yanhong，Liu Shuguo，et al.Rotordynamic analysis of whole aero-engine models based on finite element method［J］.Journal of Beijing University of Aeronautics and Astronautics，2007，33(9):1014 – 1016(in Chinese)

［11］王超，王延榮，徐星仲，等.應用三維有限元法計算汽輪機轉子臨界轉速和模化長葉片［J］.動力工程，2007，27(6):840－844 Wang Chao，Wang Yanrong，Xu Xingzhong，et al.Simulation of long blades and calculation of critical speeds of the steam turbine rotor by 3D finite element methods［J］.Journal of Power Engineering，2007，27(6):840－844(in Chinese)