趙瑞勇，楊 慧，王延榮

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

0 引言

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

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

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

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

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

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

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

2 氣動(dòng)控制方程

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

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

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

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

3 能量法-氣動(dòng)功

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

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

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

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

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

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

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

圖1 CSD網(wǎng)格

圖2 CFD網(wǎng)格

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

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

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

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

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

用Fortran編制插值程序，對(duì)NASA R67葉片前6階模態(tài)進(jìn)行了驗(yà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 動(dòng)網(wǎng)格模型

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

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

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

6 算例與結(jié)果

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

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

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

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

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

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

6.2.1 NASAR67葉片有限元計(jì)算

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

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

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

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

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

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

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

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

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

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

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

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

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

7 結(jié)論

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

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

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

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

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

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

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

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

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

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

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

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

[9]張瀟,王延榮.基于能量法的葉片顫振邊界預(yù)測(cè)方法[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è)計(jì)轉(zhuǎn)速下內(nèi)部流場(chǎng)的數(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.