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

（1.中航工業(yè)沈陽(yáng)發(fā)動(dòng)機(jī)設(shè)計(jì)研究所，沈陽(yáng)110015；2.北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京100191）

基于能量法和特征值法的顫振預(yù)測(cè)數(shù)值方法研究

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

（1.中航工業(yè)沈陽(yáng)發(fā)動(dòng)機(jī)設(shè)計(jì)研究所，沈陽(yáng)110015；2.北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京100191）

為研究適用于工程設(shè)計(jì)的顫振數(shù)值預(yù)測(cè)方法，采用能量法和特征值法對(duì)出現(xiàn)顫振現(xiàn)象的真實(shí)葉片進(jìn)行預(yù)測(cè)。通過(guò)將2種數(shù)值方法的預(yù)測(cè)結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證數(shù)值預(yù)測(cè)方法的精度和工程適用性。結(jié)果表明：2種方法都較準(zhǔn)確的判斷了某典型風(fēng)扇轉(zhuǎn)子葉片的氣彈穩(wěn)定性，與試驗(yàn)結(jié)果基本吻合，具有較高的工程實(shí)用價(jià)值。

顫振預(yù)測(cè)；能量法；特征值法；非定常氣動(dòng)力；航空發(fā)動(dòng)機(jī)

0 引言

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

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

本文基于發(fā)展的能量法和特征值法，以某發(fā)生顫振的風(fēng)扇葉片為對(duì)象，采用數(shù)值模擬的方法研究了葉片顫振的機(jī)理。

1 數(shù)值方法和求解技術(shù)

1.1 能量法

能量法是從能量交換的角度來(lái)考察葉片是否顫振，在葉片的1個(gè)振動(dòng)周期內(nèi)，葉片振動(dòng)系統(tǒng)從外界獲得的能量與系統(tǒng)阻尼消耗的能量的正負(fù)決定了葉片是否顫振。葉輪機(jī)械結(jié)構(gòu)的機(jī)械阻尼遠(yuǎn)小于氣動(dòng)阻尼[5]，因此，葉片振動(dòng)系統(tǒng)從外界獲得的能量主要是氣動(dòng)力作功，通過(guò)判斷非定常氣動(dòng)力在1個(gè)周期內(nèi)對(duì)葉片作功的正負(fù)來(lái)判定是否發(fā)生顫振。

葉片在單位面積上1個(gè)周期內(nèi)的氣動(dòng)功為

阻尼的作用就是消耗能量，在振動(dòng)分析中經(jīng)常采用能量方法將復(fù)雜的、非線性的阻尼模型簡(jiǎn)化為線性黏性阻尼[6]。對(duì)于氣動(dòng)阻尼，當(dāng)流場(chǎng)繞流情況較好，沒(méi)有強(qiáng)的非線性如漩渦等特征時(shí)，可以用等效黏性阻尼來(lái)近似表示氣動(dòng)阻尼，等效的原則是：非定常氣動(dòng)功在1個(gè)周期內(nèi)對(duì)葉片作功等于氣動(dòng)阻尼在1個(gè)周期內(nèi)消耗的能量。

式中為各方向氣動(dòng)功的疊加；qcfd為正則坐標(biāo)系下的模態(tài)振幅。

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

非定常計(jì)算完成后，通過(guò)計(jì)算1個(gè)周期內(nèi)的氣動(dòng)功及等效模態(tài)氣動(dòng)阻尼比，作為判定顫振是否發(fā)作的依據(jù)。若W＜0，ζ＞0則氣動(dòng)力作功為負(fù)，系統(tǒng)穩(wěn)定；反之則系統(tǒng)不穩(wěn)定。

1.2 特征值法

采用影響系數(shù)法[9]建立整環(huán)葉柵模型，假設(shè)葉柵在某一階模態(tài)振蕩作用下，只有1片葉片振動(dòng)，其它葉片不動(dòng)，通過(guò)非定常計(jì)算得到葉柵所有葉片與參考葉片運(yùn)動(dòng)相關(guān)的非定常氣動(dòng)力。

式中：Pn為第n個(gè)節(jié)徑的非定常氣動(dòng)力；an為n節(jié)徑對(duì)應(yīng)的氣動(dòng)力影響系數(shù)；Qn為n節(jié)徑對(duì)應(yīng)的復(fù)數(shù)振幅。在不同葉片間相位角下，通過(guò)對(duì)Pn的疊加可以得到葉片的非定常氣動(dòng)力。

將物理坐標(biāo)系下的振動(dòng)方程變換到模態(tài)坐標(biāo)系下，并進(jìn)一步變換到行波坐標(biāo)系下，可以得到

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

從式（6）和（7)可知，特征值的實(shí)部對(duì)應(yīng)振動(dòng)方程的剛度，虛部對(duì)應(yīng)振動(dòng)方程的阻尼。因此，虛部為正時(shí)，非定常氣動(dòng)力對(duì)系統(tǒng)起到穩(wěn)定作用；虛部為負(fù)時(shí)，則非定常氣動(dòng)力對(duì)系統(tǒng)起到激勵(lì)作用，系統(tǒng)不穩(wěn)定。

2 計(jì)算結(jié)果及分析

某風(fēng)扇轉(zhuǎn)子葉片發(fā)生了典型的顫振現(xiàn)象，以該葉片為模型，給定80%轉(zhuǎn)速近喘點(diǎn)附近的工況，基于能量法和特征值法，研究了葉片在1階彎曲模態(tài)振蕩作用下的氣動(dòng)彈性穩(wěn)定性，并與試驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證。

2.1 能量法計(jì)算結(jié)果及分析

采用8節(jié)點(diǎn)6面體單元，在80%轉(zhuǎn)速下，建立有限元模型如圖1所示，通過(guò)模態(tài)分析得到了該葉片1階彎曲模態(tài)的振型及其動(dòng)頻，如圖2所示。

由于葉片間相位角對(duì)氣動(dòng)彈性穩(wěn)定性的影響較大[10]，因此,建立了整環(huán)葉柵模型（如圖3所示）。給定節(jié)徑數(shù)為n，整環(huán)葉片數(shù)N，則各相鄰葉片間相位角

圖1 有限元模型

圖21 階彎曲模態(tài)（289 Hz）

圖3 整環(huán)葉柵流場(chǎng)模型

為了加快非定常計(jì)算的收斂速度，首先，進(jìn)行定常分析，并以此作為非定常的初始條件。然后計(jì)算在1階彎曲模態(tài)振蕩作用下非定常氣動(dòng)力對(duì)整環(huán)葉片的所做氣動(dòng)功和振蕩穩(wěn)定后的模態(tài)氣動(dòng)阻尼比，模態(tài)氣動(dòng)阻尼比隨節(jié)徑的變化規(guī)律如圖4所示。

圖4 模態(tài)氣動(dòng)阻尼比隨節(jié)徑變化規(guī)律

從圖4中可見(jiàn)，模態(tài)氣動(dòng)阻尼比隨節(jié)徑數(shù)增加呈現(xiàn)正弦變化規(guī)律，且在2～8節(jié)徑處，非定常氣動(dòng)力對(duì)葉片作功為正，模態(tài)氣動(dòng)阻尼比為負(fù)，葉片吸收能量，系統(tǒng)不穩(wěn)定，這與試驗(yàn)中葉片發(fā)生顫振的現(xiàn)象一致。2.2 特征值法計(jì)算結(jié)果及分析

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

圖5 氣動(dòng)力影響系數(shù)幅值分布

基于特征值法得到復(fù)數(shù)形式的特征值，特征值在復(fù)平面的分布如圖6所示。復(fù)數(shù)特征值的實(shí)部表示振動(dòng)系統(tǒng)的頻率，虛部表示振動(dòng)系統(tǒng)的氣動(dòng)阻尼，虛部的正負(fù)決定了系統(tǒng)的氣動(dòng)彈性穩(wěn)定性。從圖6中可見(jiàn)，特征值的虛部有負(fù)值出現(xiàn)，表明系統(tǒng)發(fā)生了氣動(dòng)彈性失穩(wěn)，有顫振現(xiàn)象發(fā)生，這與試驗(yàn)結(jié)果一致。

圖6 特征值在復(fù)平面的分布

在行波坐標(biāo)系下，得到了氣動(dòng)力影響系數(shù)的共軛轉(zhuǎn)置矩陣，其虛部為負(fù)表示葉柵氣動(dòng)彈性失穩(wěn)。氣動(dòng)力影響系數(shù)共軛轉(zhuǎn)置矩陣的虛部隨節(jié)徑數(shù)的分布如圖7所示。

圖7 行波坐標(biāo)系下氣動(dòng)力影響系數(shù)的虛部分布

對(duì)比圖4可以發(fā)現(xiàn)采用能量法與特征值法計(jì)算的結(jié)果相似，都在2～8節(jié)徑處，判據(jù)為負(fù)，系統(tǒng)不穩(wěn)定，發(fā)生顫振。

3 結(jié)論

基于能量法和特征值法對(duì)某發(fā)生典型顫振的風(fēng)扇轉(zhuǎn)子葉片進(jìn)行了計(jì)算分析，結(jié)果表明：采用能量法計(jì)算得到了整環(huán)葉片在2～8節(jié)徑處有正的氣動(dòng)功和負(fù)的模態(tài)氣動(dòng)阻尼比，認(rèn)為葉片發(fā)生了顫振，其葉片間相位角隨節(jié)徑變化規(guī)律與采用特征值法求得的行波坐標(biāo)系下氣動(dòng)力影響系數(shù)隨節(jié)徑數(shù)變化的規(guī)律相似。這2種方法分別從能量交換和振動(dòng)方程2方面探討了葉片顫振發(fā)作的機(jī)理，均能準(zhǔn)確預(yù)測(cè)顫振，對(duì)葉片安全設(shè)計(jì)有指導(dǎo)意義，具有一定的工程實(shí)用價(jià)值。參考文獻(xiàn)：

[1]張明明，周盛.葉輪機(jī)械葉片顫振研究的進(jìn)展與評(píng)述[J].力學(xué)進(jìn)展，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]劉文閣，馮毓誠(chéng).對(duì)經(jīng)驗(yàn)法預(yù)測(cè)葉片顫振的分析及葉片顫振數(shù)據(jù)庫(kù)的建立[J].北京航空學(xué)院學(xué)報(bào),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]周盛.葉輪機(jī)氣動(dòng)彈性力學(xué)引論 [M].北京:國(guó)防工業(yè)出版社,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]胡海巖，孫久厚，陳懷海.機(jī)械振動(dòng)與沖擊[M].北京:航空工業(yè)出版社，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ì)算研究[J].應(yīng)用力學(xué)學(xué)報(bào)，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]徐可寧.葉輪機(jī)械葉盤結(jié)構(gòu)氣動(dòng)彈性算法研究及程序?qū)崿F(xiàn)[D].北京：北京航空航天大學(xué),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]張小偉,王延榮.葉間相位角對(duì)葉片顫振的影響[J].航空動(dòng)力學(xué)報(bào)，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]張小偉.葉輪機(jī)械葉片氣彈穩(wěn)定性數(shù)值預(yù)測(cè)方法研究[D].北京：北京航空航天大學(xué),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），男，在讀碩士研究生，研究方向?yàn)楹娇瞻l(fā)動(dòng)機(jī)結(jié)構(gòu)強(qiáng)度振動(dòng)及氣動(dòng)彈性穩(wěn)定性；E-mail:364835973@qq.com。

劉一雄，劉廷毅，王德友，等.基于能量法和特征值法的顫振預(yù)測(cè)數(shù)值方法研究[J].航空發(fā)動(dòng)機(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.