趙瑞勇，楊 慧，王延榮

（1.西安航天動力研究所，西安 710100；2.北京航空航天大學(xué)能源與動力工程學(xué)院，北京 100191）

0 引言

現(xiàn)代葉輪機(jī)械正朝著高性能和大容量方向發(fā)展，由流體誘發(fā)的葉片振動問題也越來越嚴(yán)重[1]。其中顫振就是1種典型的流體誘發(fā)自激振動[2]。嚴(yán)重的葉片顫振可能破壞發(fā)動機(jī)支撐結(jié)構(gòu)，甚至導(dǎo)致整架飛機(jī)失事。因此，在發(fā)動機(jī)研制過程中，必須避免發(fā)動機(jī)在整個飛行包線范圍內(nèi)發(fā)生顫振[3]。各類葉片顫振都屬于氣動彈性穩(wěn)定性問題。氣動彈性試驗破壞性極大，因此，通常應(yīng)用計算機(jī)和CFD技術(shù)對氣彈穩(wěn)定性進(jìn)行數(shù)值仿真研究。葉輪機(jī)自身十分復(fù)雜，采用直接流固耦合（FSI）計算分析耗時長，很難讓人接受；目前，廣泛采用能量法進(jìn)行葉輪機(jī)械顫振研究，通過計算以轉(zhuǎn)子振動的某一自然振型表現(xiàn)的顫振與流場之間的能量平衡來預(yù)測顫振。計算中通常忽略非定常氣動力對振動模態(tài)的影響。因此，該方法實質(zhì)上是不耦合，或精確地說是弱耦合的，也就不包括非定常氣動力對振型和頻率的反饋影響?；谀芰糠ǖ暮喕匠讨詰?yīng)用廣泛，主要是因為作用于葉片的非定常氣動力與力矩對金屬葉片的振動模態(tài)與頻率影響較小，可以獨立進(jìn)行結(jié)構(gòu)動力學(xué)分析與非定常氣動分析，極大簡化了分析工作。在葉片氣動彈性仿真計算中，由于結(jié)構(gòu)變形，每一計算時間步都需要生成自適應(yīng)計算網(wǎng)格。采用合理的網(wǎng)格更新方法[6]和尋求高效率、高精度的CSD/CFD數(shù)據(jù)交換方法是耦合技術(shù)應(yīng)用的關(guān)鍵。

本文針對剛性葉片由于變形較小，發(fā)展了3維線性插值理論來進(jìn)行CSD/CFD耦合界面數(shù)據(jù)交換。

1 葉片結(jié)構(gòu)動力學(xué)方程

本文使用有限元法計算葉片的固有振動特性，未考慮葉片間的相互影響，假設(shè)每個葉片的振動頻率都相同，模型不計機(jī)械阻尼的影響。葉片離心負(fù)載下振動的結(jié)構(gòu)動力學(xué)控制方程為

式中：[M]為單元質(zhì)量矩陣；[KL]為單元彈性剛度矩陣；[KG]為單元離心剛度矩陣；（σ0）為離心力引起的初始應(yīng)力向量；[MG]為單元離心質(zhì)量矩陣；{x¨}、{x}分別為單元的加速度向量和位移向量；Ω為葉片的轉(zhuǎn)動角速度；{FC（Ω2）}為單元節(jié)點平衡力向量，此處為離心分布的離散形式。

使用通用軟件ANSYS來計算式（1）。在結(jié)構(gòu)分析中，葉片固連在輪轂上，在葉根處為零線位移和零角位移約束，葉片其它位置為自由約束。給出葉片的各階振型和固有頻率，流場中葉片按照第i階振型和固有頻率進(jìn)行諧振動

式中：Z（x,y,z,t ）為葉片表面各點的位移矢量；Φi（x,y,z）為葉片的第i階振型矢量；q0為廣義坐標(biāo)振幅；ωi為第i階振型的固有頻率。

2 氣動控制方程

假設(shè)葉片在做簡諧振動的情況下引入運動網(wǎng)格，并對柱坐標(biāo)系下的雷諾平均N－S方程進(jìn)行求解來得到振動葉片所受到的非定常氣動力以及氣動力所做的功。

柱坐標(biāo)系（φ，r，z ）下帶有角速度為ω和有限體積表面運動速度Vg（ug,vg,wg）的3維積分型N－S方程為

式中：U為氣動參數(shù)；E、F、G為無黏通量；Sl為源項。具體表達(dá)式見參考文獻(xiàn)[7]。

湍流模型采用對分離模擬較好的SST模型[4],方程求解采用CFD軟件FLUENT。

3 能量法-氣動功

考察在1個振動周期內(nèi)，流體與葉片發(fā)生能量交換，如果非定常氣動力對葉片作正功，則葉片從氣流中吸收能量，發(fā)生顫振[8]。

整個葉片的氣動功為葉片單位面積上1個周期內(nèi)氣動功在葉片表面上的積分

定義等效模態(tài)氣動阻尼[10]

式中：xcfd為葉片實際振幅；φ為模態(tài)振幅。

按照上述的能量法，當(dāng)非定常氣動力總功為正值時，模態(tài)氣動阻尼為負(fù)，葉片振動發(fā)散。

4 CSD/CFD數(shù)據(jù)插值

由于氣動計算網(wǎng)格比結(jié)構(gòu)計算網(wǎng)格密，如圖1、2所示，要進(jìn)行氣動與結(jié)構(gòu)耦合計算，快捷、高質(zhì)量的插值方法對計算時效和計算精度至關(guān)重要。本文改進(jìn)了文獻(xiàn)[9]中的3維線性插值算法，在程序中引入面積控制因子，減小了插值誤差。該算法易懂，比復(fù)雜插值方法效率高。非常適用于非結(jié)構(gòu)網(wǎng)格的計算。

圖1 CSD網(wǎng)格

圖2 CFD網(wǎng)格

該算法每個CFD格點選距離最近的10個FE網(wǎng)格點，引入面積控制因子選出其中不小于其值的最近的4個點，4個FE網(wǎng)格點每3個就可以組成1個平面三角形，求得每個平面三角形的單元中心，選定單元中心離CFD網(wǎng)格點最近的三角形平面作為插值平面△A1B1C1，如圖3所示。

圖3 FE網(wǎng)格點可選擇插值面

每個CFD網(wǎng)格點沿選定平面法向量投影如圖4所示。

圖4 CFD網(wǎng)格點投影到插值面

在點 A1、B1、C1上 分別作△A1B1C1的法線，使得 A1A2、B1B2、C1C2分別為點 A1、B1、C1的位移值。過點F做△A1B1C1的法線，在△A1B1C1和△A2B2C2上的交點分別為F1、F2，則F1F2為F點的位移值。如圖5所示。

用Fortran編制插值程序，對NASA R67葉片前6階模態(tài)進(jìn)行了驗證，都得到了很好的結(jié)果。以下僅以第1階模態(tài)X方向的位移為例。插值前FE網(wǎng)格第1階模態(tài)X方向位移如圖6所示，插值后CFD網(wǎng)格模態(tài)位移等值線如圖7所示。從圖中可見，插值前、后數(shù)值一致。

圖5 3維線性插值

圖6 固體域FE網(wǎng)格第1階X方向模態(tài)位移

圖7 流體域CFD網(wǎng)格插值結(jié)果

5 動網(wǎng)格模型

在葉片氣動彈性仿真計算中，葉片表面各節(jié)點每個時間步以某階模態(tài)振型按固有模態(tài)頻率做正弦運動，這樣由于結(jié)構(gòu)變形，每一計算時間步，周圍流場需要生成自適應(yīng)的計算網(wǎng)格。本文編制UDF程序用DEFINE_GRID_MOTION（）宏命令指定葉片表面各節(jié)點的運動形式，實現(xiàn)葉片的某階模態(tài)振動。

由于計算模型單元數(shù)較多，非定常計算量大，為保證黏性計算精度和葉片表面附近網(wǎng)格隨著葉片振動畸變小，采用彈簧光順動網(wǎng)格模型，更新葉片周圍流場網(wǎng)格。彈簧光順網(wǎng)格模型更新過程如圖8所示。

圖8 彈簧光順網(wǎng)格模型更新

6 算例與結(jié)果

6.1 NASA R67風(fēng)扇葉片定常場計算與結(jié)果分析

本例計算模型選取1個扇區(qū)通道，對稱面加周期邊界條件，如圖9所示。采用6面體非結(jié)構(gòu)網(wǎng)格，共317833個節(jié)點，299509個單元；葉片表面10198個節(jié)點，10072個單元。不考慮葉尖間隙的影響，壁面設(shè)定為絕熱條件。為了加速收斂，計算采用多重網(wǎng)格技術(shù)，V循環(huán)。在定常狀態(tài)下，入口總壓為1個標(biāo)準(zhǔn)大氣壓，風(fēng)扇葉片轉(zhuǎn)速按照設(shè)計轉(zhuǎn)速給定，通過改變流量得到風(fēng)扇性能曲線。試驗得到的堵塞流量為34.96 kg/s，計算得到的堵塞流量為34.77 kg/s。不同流量下總壓比計算值與試驗值[12]的對比如圖10所示，試驗值和計算值分別按照各自的堵塞流量進(jìn)行規(guī)一。從圖中可見，二者吻合較好，但是總體上計算壓比值比試驗值稍低。近最高效率點下的計算值和設(shè)計值見表1。從表中可見，二者誤差較小。

表1 NASA R67定常計算結(jié)果

近最高效率點90%葉高處馬赫數(shù)如圖11所示。從圖中可見，90%葉高處葉片前緣斜激波和葉道中的垂直波呈“λ”形狀，同時還可以觀察到緊靠垂直激波后面有1道弱激波。

在設(shè)計轉(zhuǎn)速下，葉片前緣壓力值比后緣的低，吸力面上的激波波腳呈明顯的λ形，如圖12所示。

6.2 NASA R67葉片氣彈穩(wěn)定性計算與結(jié)果

6.2.1 NASAR67葉片有限元計算

葉片模型單元類型為solid45，X、Y、Z 3個方向共2×20×22個單元，共1449個節(jié)點?？紤]離心力，轉(zhuǎn)速為 16043r/min，泊松比 μ=0.3，彈性模 E=1.12E+11，選取鈦合金材料密度ρ=4440 kg/m3。約束方式為葉根固支全約束，葉尖自由。葉片模態(tài)和頻率如圖13所示。

6.2.2 第1階模態(tài)下葉片周期氣動功分布

品牌農(nóng)業(yè)區(qū)域效應(yīng)逐步顯現(xiàn)。臥龍區(qū)通過發(fā)展龍頭企業(yè)、培育農(nóng)業(yè)品牌等舉措，充分發(fā)揮品牌農(nóng)業(yè)企業(yè)的產(chǎn)業(yè)鏈優(yōu)勢，提高標(biāo)準(zhǔn)化、區(qū)域化、產(chǎn)業(yè)化經(jīng)營水平，提高農(nóng)業(yè)經(jīng)濟(jì)整體效益，有力地推進(jìn)區(qū)域經(jīng)濟(jì)的增長。如石橋的月季，現(xiàn)已發(fā)展到臥龍區(qū)的各個鄉(xiāng)鎮(zhèn)；石橋的老姜歷史悠久，姜質(zhì)優(yōu)良，在石橋的各個村都有不同規(guī)模的種植；謝莊龔河的玫瑰花、董營的紅薯已發(fā)展到周邊村種植；蒲山的桃、潦河坡的石榴、潦河的葡萄等因特色不同，形成了不同的旅游觀光區(qū)，而且由于產(chǎn)品品質(zhì)優(yōu)良，在南陽乃至武漢、上海等大中城市備受歡迎，品質(zhì)優(yōu)良價位高，品牌效應(yīng)有所顯現(xiàn)。

本文計算了前3階模態(tài)氣動功分布，限于篇幅僅給出第1階彎曲模態(tài)下耦合計算結(jié)果。

將有限元結(jié)果插值到如圖2所示的CFD網(wǎng)格節(jié)點上，得到葉片表面位移數(shù)據(jù)，調(diào)用FLUENT動網(wǎng)格模型和UDF程序，并對模態(tài)位移進(jìn)行歸一化，指定流場中葉片實際最大振幅為1 mm進(jìn)行葉片第1階模態(tài)振動下的非定常計算。葉片前3階模態(tài)振動下的總氣動功和氣動阻尼見表2。

表2 前3階模態(tài)氣動功和氣動阻尼比

根據(jù)能量法，由表2計算結(jié)果可知：前3階模態(tài)振動非定常氣動力對葉片作功為負(fù)，起到阻尼作用，系統(tǒng)能量耗散，葉片不會發(fā)生氣彈失穩(wěn)。當(dāng)振幅為1 mm時,在第1階模態(tài)（第1階彎曲）和第3階模態(tài)（第1階扭轉(zhuǎn)）下的氣動阻尼較大，在第2階模態(tài)（以彎為主的彎扭耦合振動）下的氣動阻尼較小。

在第1階模態(tài)下，葉片監(jiān)視點位移與非定常氣動力隨時間變化如圖14所示，由圖14（b、c、d）中可見，非定常氣動力較振蕩葉片位移變化滯后一定相位角，該相位角是葉片表面非定常氣動功變化因素之一。為了便于比較，放大了氣動力幅值倍數(shù)。

葉片壓力面和吸力面1個振動周期累積氣動功分布如圖15所示。在不同徑向位置（10%、30%、50%、70%、90%葉高），氣動功沿葉片弦向分布如圖16～20所示。文中選取定常場計算70%和90%葉高處相對馬赫數(shù)沿弦向分布作為激波參考位置，如圖21、22所示。

對比圖 15（a）、（b）可知，即使第1階彎曲振動葉片表面總的氣動功為負(fù)值，壓力面和吸力面局部仍存在氣動功為正值的區(qū)域（約70%～90%葉高，弦向0.6～0.9區(qū)域），對此應(yīng)予以關(guān)注。

從圖16、17中可見，在第1階模態(tài)振動、50%葉高以下，由于激波的減弱或消失，壓力面、吸力面氣動功沿弦向受激波影響小，基本呈對稱分布。

從圖19、20中可見，在70%、90%葉高處，壓力面、吸力面沿著弦向出現(xiàn)了絕對值較大的正功區(qū)。對比圖19、21可知，在70%葉高激波位置處（約弦向0.35），氣動功出現(xiàn)了波動；同樣，對比圖20、22可知，在90%葉高激波位置也出現(xiàn)了氣動功的波動。

由此可見，在跨聲速流動中激波的位置對葉片氣彈穩(wěn)定性有較大影響。

7 結(jié)論

（1）在前3階模態(tài)振動下，非定常氣動力對葉片作負(fù)功，氣動阻尼為正值，葉片不會發(fā)生顫振。

（2）基于3維線性插值理論發(fā)展的CSD/CFD網(wǎng)格數(shù)據(jù)交換程序，可以實現(xiàn)流固耦合計算，該程序?qū)W(wǎng)格限制性小，經(jīng)對比發(fā)現(xiàn)插值結(jié)果誤差小，從而減小了數(shù)據(jù)傳遞帶來的結(jié)果誤差。

（3）使用Fluent動網(wǎng)格模型對振蕩葉片擾流進(jìn)行了仿真計算，給出了葉片表面非定常周期氣動功的分布情況，可為預(yù)測顫振提供依據(jù)。

（4）通過分析不同徑向位置葉片沿弦線方向的氣動功分布得出：激波位置和位移、氣動力間相位角是影響非定常氣動功分布的重要因素。

[1]袁新，金琰，暢國勇.葉輪機(jī)械中流體激振問題的流固耦合研究[C]//第八屆全國空氣彈性學(xué)術(shù)交流會論文集，北京：中國空氣動力學(xué)會，2003：274-280.

[2]宋兆泓.葉片顫振[J].水利水電機(jī)械，2006（5）：45-52.

[3]蔣福慶.風(fēng)扇-壓氣機(jī)亞、跨音失速顫振的發(fā)作機(jī)理及預(yù)測技術(shù)[J].航空發(fā)動機(jī)，1996（3）：9-15.

[4]張陳安，史愛明，劉鋒，等.基于SST湍流模型的三維葉片氣動彈性問題研究 [C]//第十屆全國空氣彈性學(xué)術(shù)交流會論文集，北京：中國力學(xué)學(xué)會，2007：345-351.

[5]王延榮,郝燕平,宋兆泓.渦輪機(jī)械葉片顫振及其研究進(jìn)展[C]//中國航空學(xué)會第九屆航空發(fā)動機(jī)結(jié)構(gòu)強(qiáng)度振動學(xué)術(shù)交流會論文集，北京：中國航空學(xué)會，1998:265-271.

[6]徐敏，陳士櫓.CFD/CSD耦合計算研究 [J].應(yīng)用力學(xué)學(xué)報，2004，21（2）：33-36.

[7]胡運聰,周新海.振動葉柵非定常流動數(shù)值模擬與葉片顫振分析[D].西安:西北工業(yè)大學(xué),2003.

[8]趙瑞勇,楊慧，王延榮.葉片顫振數(shù)值模擬方法研究[D].北京:北京航空航天大學(xué),2010.

[9]張瀟,王延榮.基于能量法的葉片顫振邊界預(yù)測方法[D].北京:北京航空航天大學(xué),2008.

[10]Moffatt S，Li He. Blade forced response prediction for industrial gas turbines,Part 1: methodologies[R].ASME 2003- GT- 38640.

[11]楊策,老大中,蔣滋康.風(fēng)扇轉(zhuǎn)子在設(shè)計轉(zhuǎn)速下內(nèi)部流場的數(shù)值研究[J].自然科學(xué)進(jìn)展,1999,9（增刊1）：1312-1317.

[12]Anthony J S, Jerry R W, Michael D H，et al. Laser anemometer measurements in a transonic axial- flow fan rotor [R].NASA- TP- 1989- 2879.