基于能量法和特征值法的顫振預測數值方法研究

劉一雄，劉廷毅，王德友，叢佩紅，王延榮

（1.中航工業沈陽發動機設計研究所，沈陽110015；2.北京航空航天大學能源與動力工程學院，北京100191）

基于能量法和特征值法的顫振預測數值方法研究

劉一雄1，劉廷毅1，王德友1，叢佩紅1，王延榮2

（1.中航工業沈陽發動機設計研究所，沈陽110015；2.北京航空航天大學能源與動力工程學院，北京100191）

為研究適用于工程設計的顫振數值預測方法，采用能量法和特征值法對出現顫振現象的真實葉片進行預測。通過將2種數值方法的預測結果與試驗結果進行對比，驗證數值預測方法的精度和工程適用性。結果表明：2種方法都較準確的判斷了某典型風扇轉子葉片的氣彈穩定性，與試驗結果基本吻合，具有較高的工程實用價值。

顫振預測；能量法；特征值法；非定常氣動力；航空發動機

0 引言

隨著設計要求的提高，航空發動機日益向高負荷、高效率、高可靠性、高推重比的方向發展。而隨著推重比的提高，級壓比增加，葉片變得更薄，非定常氣動載荷加大，再加上新型輕質材料的應用，使得葉片剛度大幅降低，導致發動機葉片振動尤其是顫振問題更加突出。葉輪機械葉片顫振預測方法從早期基于相似理論的經驗法發展到了當前的數值預測方法[1]。由于經驗法過分依賴于試驗數據的積累和經驗的判斷，不能準確預測一些新型壓氣機/風扇葉片。因此，發展適用于工程設計的顫振數值預測方法對發動機安全設計、縮短設計周期、節省經費等都有深遠意義。葉片顫振問題的數值研究方法主要分為經典方法和耦合方法[2]。經典方法通過對氣動及結構模型作一定的簡化，忽略流體與結構的耦合關系，將流體與結構分開求解，其主要包括能量法和特征值法。由Carta[3]提出的能量法基于以下假設：顫振以轉子振動的某一自然振型出現，通過計算這一振型與流場之間的能量平衡來預測；特征值法是以葉片的振動方程為基礎，將葉片的結構動力方程與非定常氣動特性結合起來，形成統一的氣動彈性方程，將穩定性問題轉化為特征值問題，根據特征值虛部的正、負來判斷葉片是否發生顫振。

雖然能量法與特征值法求解顫振問題并不能真實描述顫振的物理本質及流固耦合關系，但是隨著計算能力的提高，通過對模型的適當修改并采用合理的假設，是可以實現所謂的弱耦合[4]的，并且計算量小、計算周期短、精度較高，是比較不錯的工程選擇。

本文基于發展的能量法和特征值法，以某發生顫振的風扇葉片為對象，采用數值模擬的方法研究了葉片顫振的機理。

1 數值方法和求解技術

1.1 能量法

能量法是從能量交換的角度來考察葉片是否顫振，在葉片的1個振動周期內，葉片振動系統從外界獲得的能量與系統阻尼消耗的能量的正負決定了葉片是否顫振。葉輪機械結構的機械阻尼遠小于氣動阻尼[5]，因此，葉片振動系統從外界獲得的能量主要是氣動力作功，通過判斷非定常氣動力在1個周期內對葉片作功的正負來判定是否發生顫振。

葉片在單位面積上1個周期內的氣動功為

阻尼的作用就是消耗能量，在振動分析中經常采用能量方法將復雜的、非線性的阻尼模型簡化為線性黏性阻尼[6]。對于氣動阻尼，當流場繞流情況較好，沒有強的非線性如漩渦等特征時，可以用等效黏性阻尼來近似表示氣動阻尼，等效的原則是：非定常氣動功在1個周期內對葉片作功等于氣動阻尼在1個周期內消耗的能量。

式中為各方向氣動功的疊加；qcfd為正則坐標系下的模態振幅。

在進行非定常計算時，固體網格比流體網格稀疏很多，通過將固體節點上的位移插值到流體節點上實現數據交換[7]。由于葉片的振動會引起流場變化，且對全部流場進行網格點更新會浪費大量時間，因此,為了反映葉片振動對流場的作用同時為了保證效率，Slone[8]等考慮將流體域分為固定域和可動域，這樣只需要實時更新可動域的網格就可以較好地模擬流場與固體交互影響的振蕩過程。

非定常計算完成后，通過計算1個周期內的氣動功及等效模態氣動阻尼比，作為判定顫振是否發作的依據。若W＜0，ζ＞0則氣動力作功為負，系統穩定；反之則系統不穩定。

1.2 特征值法

采用影響系數法[9]建立整環葉柵模型，假設葉柵在某一階模態振蕩作用下，只有1片葉片振動，其它葉片不動，通過非定常計算得到葉柵所有葉片與參考葉片運動相關的非定常氣動力。

式中：Pn為第n個節徑的非定常氣動力；an為n節徑對應的氣動力影響系數；Qn為n節徑對應的復數振幅。在不同葉片間相位角下，通過對Pn的疊加可以得到葉片的非定常氣動力。

將物理坐標系下的振動方程變換到模態坐標系下，并進一步變換到行波坐標系下，可以得到

將式（5）代入式（4)可以得到

從式（6）和（7)可知，特征值的實部對應振動方程的剛度，虛部對應振動方程的阻尼。因此，虛部為正時，非定常氣動力對系統起到穩定作用；虛部為負時，則非定常氣動力對系統起到激勵作用，系統不穩定。

2 計算結果及分析

某風扇轉子葉片發生了典型的顫振現象，以該葉片為模型，給定80%轉速近喘點附近的工況，基于能量法和特征值法，研究了葉片在1階彎曲模態振蕩作用下的氣動彈性穩定性，并與試驗結果進行對比驗證。

2.1 能量法計算結果及分析

采用8節點6面體單元，在80%轉速下，建立有限元模型如圖1所示，通過模態分析得到了該葉片1階彎曲模態的振型及其動頻，如圖2所示。

由于葉片間相位角對氣動彈性穩定性的影響較大[10]，因此,建立了整環葉柵模型（如圖3所示）。給定節徑數為n，整環葉片數N，則各相鄰葉片間相位角

圖1 有限元模型

圖21 階彎曲模態（289 Hz）

圖3 整環葉柵流場模型

為了加快非定常計算的收斂速度，首先，進行定常分析，并以此作為非定常的初始條件。然后計算在1階彎曲模態振蕩作用下非定常氣動力對整環葉片的所做氣動功和振蕩穩定后的模態氣動阻尼比，模態氣動阻尼比隨節徑的變化規律如圖4所示。

圖4 模態氣動阻尼比隨節徑變化規律

從圖4中可見，模態氣動阻尼比隨節徑數增加呈現正弦變化規律，且在2～8節徑處，非定常氣動力對葉片作功為正，模態氣動阻尼比為負，葉片吸收能量，系統不穩定，這與試驗中葉片發生顫振的現象一致。2.2 特征值法計算結果及分析

特征值法定常計算采用能量法的結果，并以此作為非定常的初始條件。進行非定常計算分析時，假定參考葉片以1階彎曲模態振動，其他葉片不動，進行非定常計算分析。待振蕩穩定后，提取物理坐標系下各葉片表面的非定常氣動力及其振動位移，并將之轉換到模態坐標系下，得到各葉片的非定常氣動力影響系數的幅值分布，如圖5所示。從圖5中可見，參考葉片的振動對其自身的影響最大，對其相鄰2個葉片的振動影響次之，對其它葉片的影響較小，基本可以忽略。因此，氣動力影響系數矩陣采用5對角矩陣即可滿足計算精度[11]。

圖5 氣動力影響系數幅值分布

基于特征值法得到復數形式的特征值，特征值在復平面的分布如圖6所示。復數特征值的實部表示振動系統的頻率，虛部表示振動系統的氣動阻尼，虛部的正負決定了系統的氣動彈性穩定性。從圖6中可見，特征值的虛部有負值出現，表明系統發生了氣動彈性失穩，有顫振現象發生，這與試驗結果一致。

圖6 特征值在復平面的分布

在行波坐標系下，得到了氣動力影響系數的共軛轉置矩陣，其虛部為負表示葉柵氣動彈性失穩。氣動力影響系數共軛轉置矩陣的虛部隨節徑數的分布如圖7所示。

圖7 行波坐標系下氣動力影響系數的虛部分布

對比圖4可以發現采用能量法與特征值法計算的結果相似，都在2～8節徑處，判據為負，系統不穩定，發生顫振。

3 結論

基于能量法和特征值法對某發生典型顫振的風扇轉子葉片進行了計算分析，結果表明：采用能量法計算得到了整環葉片在2～8節徑處有正的氣動功和負的模態氣動阻尼比，認為葉片發生了顫振，其葉片間相位角隨節徑變化規律與采用特征值法求得的行波坐標系下氣動力影響系數隨節徑數變化的規律相似。這2種方法分別從能量交換和振動方程2方面探討了葉片顫振發作的機理，均能準確預測顫振，對葉片安全設計有指導意義，具有一定的工程實用價值。參考文獻：

[1]張明明，周盛.葉輪機械葉片顫振研究的進展與評述[J].力學進展，2011，41（1）：26-38. ZHANG Mingming,ZHOU Sheng.A review of the research of blade flutter in turbomachinery[J].Advances in Machines，2011,41（1）：26-38.(in Chinese)

[2]劉文閣，馮毓誠.對經驗法預測葉片顫振的分析及葉片顫振數據庫的建立[J].北京航空學院學報,1986（4）:113-122. LIU Wenge,FENG Yucheng.Analysis of empirical flutter predicition method of blades and establishment of a flutter data bank[J].Journal of Beijing Institute of Aeronautics and Astronautics,1986（4）:113-122.(in Chinese)

[3]周盛.葉輪機氣動彈性力學引論 [M].北京:國防工業出版社,1986：241. ZHOU Sheng.Introduction of aeroelastic in turbomachinery[M]. Beijing:National Defence Industry Press,1986:241.(in Chinese)

[4]Marshall J G,Imregum M A.Review of aeroelasticity methods with emphasis on turbomachinery applications[J].Journal of Fluids and Structure,1996,10（3）:237-267.

[5]Catra F O.Coupled blade-disk-shroud flutter instabilities in turbojet engine rotors[J].Journal of Engineering for Power, 1967,89（3）:419-426

[6]Kaza K R V,Kielb R E.Effects of mistuning on bendingtorsion flutter and response of a cascade in imcompressible flow[R]. AIAA-1981-0602.

[7]Stuart M,LI He.Blade forced response prediction for industrial gas turbines,part I:methodologies[C]//Proceeding of ASME Turbo Expo:Power for Land,Sea,and Air,2003:38640.

[8]Kielb R E,Chiang H D.Recent advancements in turbomachinery forced response analysis[R].AIAA-1992-0012.

[9]胡海巖，孫久厚，陳懷海.機械振動與沖擊[M].北京:航空工業出版社，1998:285. HU Haiyan,SUN Jiuhou,CHEN Huaihai.Mechanical vibration and shock[M].Beijing:Aviation Industry Press,1998：285.(in Chinese)

[10]徐敏，陳士櫓.CFD/CSD耦合計算研究[J].應用力學學報，2003，21（2）:33-36. XU Min,CHEN Shilu.Study of data exchange method for coupling computational CFD/CSD[J].Chinese Journal of Applied Mechanics,2003,21（2）:33-36.(in Chinese)

[11]Slone A K,Pericleous K,Bailey C,et al.A finite volume unstructured mesh approach to dynamic fluidstructure interaction:an assessment of the challenge of predicting the onset of flutter[J].Applied Mathematical Modelling,2004,28:211-239. [12]Hanamura Y,Tanaka H,Yamaguchi K.A simplified method to measure unsteady forces acting on the vibrating blades in cascade[J].Bulletin of the JSME,1980,23:880-887.

[13]徐可寧.葉輪機械葉盤結構氣動彈性算法研究及程序實現[D].北京：北京航空航天大學,2011. XU Kening.On the algorithm of aeroelasticity and the program implementation of bladed disk assemblies in turbomachines[D].Beijing:Beihang University,2011.(in Chinese)

[14]張小偉,王延榮.葉間相位角對葉片顫振的影響[J].航空動力學報，2010，25（2）:412-416. ZHANG Xiaowei,WANG Yanrong.Influence of interblade phase angle on the flutter of rotor blades[J].Journal of Aerospace Power,2010,25（2）:412-416.(in Chinese)

[15]張小偉.葉輪機械葉片氣彈穩定性數值預測方法研究[D].北京：北京航空航天大學,2011. ZHANG Xiaowei.Numerical prediction method for aeroelastic stability of blades in turbomachines[D].Beijing：Beihang University,2011.(in Chinese)

Numerical Study of Flutter Prediction Based on Energy Method and Eigenvalue Method

LIU Yi-xiong1,LIU Ting-yi1,WANG De-you1,CONG Pei-hong1,WANG Yan-rong2
（1.AVIC Shenyang Engine Design and Research Institute,Shenyang 110015,China;2.School of Jet Propulsion,Beihang University, Beijing 100083,China）

In order to study the numerical method of flutter prediction applied on the engineering design,the flutter phenomenon of a fan blade was predicted by energy method and eigenvalue method.The precision and engineering applicability of the numerical methods on flutter prediction were validated by comparing the predicted results and test results.The results show that the energy method and eigenvalue method are accuracy to diagnosis the aeroelastic stability of a fan blade,which in good agreement with the test results and are of good engineering practical value.

flutter prediction;energy method;eigenvalue method;unsteady aerodynamic forces;aeroengine

V 211.1+5

A

10.13477/j.cnki.aeroengine.2014.06.009

2013-07-06

劉一雄（1988），男，在讀碩士研究生，研究方向為航空發動機結構強度振動及氣動彈性穩定性；E-mail:364835973@qq.com。

劉一雄，劉廷毅，王德友，等.基于能量法和特征值法的顫振預測數值方法研究[J].航空發動機，2014,40（6）:43-46.LIUYixiong,LIUTingyi, WANGDeyou,et al.Numerical studyofflutter prediction based on energymethod and engevalue method[J].Aeroengine，2014,40（6）:43-46.