寬弦風(fēng)扇葉片顫振預(yù)測(cè)的工程研究

章嘉麟，丁建國(guó)

(中國(guó)航發(fā)商用航空發(fā)動(dòng)機(jī)有限責(zé)任公司，上海201108)

1 引言

航空發(fā)動(dòng)機(jī)壓縮部件的顫振多發(fā)生于大展弦比、小剛度葉片中，一旦發(fā)作，可在極短時(shí)間內(nèi)導(dǎo)致葉片疲勞斷裂失效[1]。隨著民用航空發(fā)動(dòng)機(jī)持續(xù)向高效率、高可靠性、高推重比的方向發(fā)展[2]，發(fā)動(dòng)機(jī)風(fēng)扇葉片尺寸不斷增大，各種新結(jié)構(gòu)和新材料也在風(fēng)扇葉片上得以應(yīng)用，如鈦合金空心風(fēng)扇葉片、復(fù)合材料風(fēng)扇葉片等。這些技術(shù)在減輕葉片質(zhì)量、提高推重比的同時(shí)，也大大降低了葉片剛度，導(dǎo)致風(fēng)扇葉片顫振問(wèn)題更加突出。因此，防止風(fēng)扇葉片顫振已經(jīng)成為大涵道比民用航空發(fā)動(dòng)機(jī)研制過(guò)程中所面臨的一個(gè)重要問(wèn)題，發(fā)展適用于工程設(shè)計(jì)的風(fēng)扇葉片顫振預(yù)測(cè)方法，對(duì)于發(fā)動(dòng)機(jī)安全設(shè)計(jì)、縮短設(shè)計(jì)周期、節(jié)省研制經(jīng)費(fèi)[3-4]都具有深遠(yuǎn)意義。

早期的顫振預(yù)測(cè)主要采用經(jīng)驗(yàn)法和半經(jīng)驗(yàn)法[5]進(jìn)行。這類方法以相似理論為指導(dǎo)，將大量的葉片顫振統(tǒng)計(jì)數(shù)據(jù)整理成經(jīng)驗(yàn)性預(yù)測(cè)準(zhǔn)則或經(jīng)驗(yàn)型曲線，引入工程試驗(yàn)數(shù)據(jù)，用以預(yù)測(cè)新設(shè)計(jì)葉片的顫振穩(wěn)定性[6]。近年來(lái)，由于計(jì)算機(jī)性能的高速發(fā)展和非定常流場(chǎng)模擬技術(shù)的日漸成熟，基于CFD的方法逐漸成為顫振預(yù)報(bào)的主流[7-9]。葉輪機(jī)械顫振問(wèn)題的數(shù)值研究方法可分為經(jīng)典法和耦合法[10-11]兩種。經(jīng)典法通過(guò)對(duì)葉輪機(jī)械結(jié)構(gòu)和氣動(dòng)模型的適當(dāng)簡(jiǎn)化，忽略流體與結(jié)構(gòu)的耦合關(guān)系，將結(jié)構(gòu)與流體分開(kāi)求解。常用的經(jīng)典法有特征值法和能量法，Bendiksen等[12]采用特征值法對(duì)多種葉柵顫振問(wèn)題進(jìn)行了研究，Hall[13]和He[14]等則用能量法對(duì)三維葉片顫振進(jìn)行了分析。耦合法的流體和結(jié)構(gòu)求解是交互影響的，使得流體與結(jié)構(gòu)求解能夠考慮更多的非線性因素。Debrabandere等[15]采用耦合法對(duì)某壓氣機(jī)轉(zhuǎn)子葉片在設(shè)計(jì)工況下的氣動(dòng)彈性穩(wěn)定性進(jìn)行了研究。Schoenenborn等[16]采用三維有限元模型模擬葉片，利用Navier-Stokes方程求解流場(chǎng)，研究了壓氣機(jī)在喘振工況下的葉片顫振問(wèn)題。胡運(yùn)聰[17]和楊真青[18]等運(yùn)用耦合法對(duì)二維葉柵顫振問(wèn)題進(jìn)行了求解，指出振動(dòng)葉柵中的流動(dòng)具有較強(qiáng)的非線性特征，與非耦合情況存在較大差異。由于流動(dòng)的非定常性和求解的復(fù)雜性，完全耦合計(jì)算過(guò)程需要投入大量的資源，國(guó)內(nèi)尚未在工程設(shè)計(jì)中廣泛采用。

某大涵道比風(fēng)扇葉片在設(shè)計(jì)中采用了寬弦和復(fù)合彎掠等技術(shù)，在性能考核試驗(yàn)過(guò)程中風(fēng)扇葉片發(fā)生了顫振現(xiàn)象。本文分別基于經(jīng)驗(yàn)法和數(shù)值模擬方法，以發(fā)生顫振的風(fēng)扇葉片為對(duì)象進(jìn)行顫振評(píng)估。以期通過(guò)對(duì)其的分析研究，掌握該類葉片的顫振機(jī)理；同時(shí)，將計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比，以驗(yàn)證顫振預(yù)測(cè)方法的準(zhǔn)確性。

2 風(fēng)扇葉片顫振預(yù)估方法

2.1 經(jīng)驗(yàn)法

經(jīng)驗(yàn)法[19]包括單參數(shù)法、雙參數(shù)法和三參數(shù)法等。對(duì)于傳統(tǒng)的三參數(shù)法，通常只計(jì)算單一特征截面的折合速度或折合頻率、進(jìn)口相對(duì)馬赫數(shù)以及氣動(dòng)攻角，并與試驗(yàn)得到的經(jīng)驗(yàn)曲線族進(jìn)行比較，判斷顫振是否發(fā)生。本文采用經(jīng)過(guò)改進(jìn)的三參數(shù)法，選擇兩個(gè)特征截面，對(duì)風(fēng)扇葉片進(jìn)行顫振分析。計(jì)算過(guò)程為：選擇75%和100%葉高處的兩個(gè)特征截面，計(jì)算1～3階的75%葉高處的折合速度和100%葉高處的模態(tài)振型，分別以折合速度和模態(tài)振型為橫、縱坐標(biāo)在經(jīng)驗(yàn)曲線圖上查找得出顫振臨界攻角。如果75%葉高處的實(shí)際氣動(dòng)攻角小于臨界攻角，則葉片穩(wěn)定；反之，則葉片發(fā)生顫振。與傳統(tǒng)意義上的三參數(shù)法相比，改進(jìn)后的三參數(shù)法增加了葉尖截面的模態(tài)振型計(jì)算，考慮了葉片尖部的扭轉(zhuǎn)和變形，有利于提高預(yù)測(cè)精度。

2.2 數(shù)值模擬方法

本文在對(duì)風(fēng)扇葉片顫振進(jìn)行數(shù)值模擬預(yù)測(cè)中，采用基于頻域方法的能量法，通過(guò)引進(jìn)合理假設(shè)，實(shí)現(xiàn)流體和固體的弱耦合[20]，計(jì)算量和計(jì)算精度都能滿足工程需求。假設(shè)所有葉片以相同的頻率和振幅做簡(jiǎn)諧振動(dòng)，相鄰葉片的振動(dòng)相差一個(gè)相同的葉片間相位角(IBPA)，計(jì)算一個(gè)周期內(nèi)非定常氣動(dòng)力對(duì)葉片所做的功。計(jì)算使用一方程S-A湍流模型；對(duì)流項(xiàng)空間離散采用中心差分+二階/四階人工粘性；虛擬時(shí)間項(xiàng)積分采用顯式(Runge-Kutta)和隱式(LU-SGS)相結(jié)合的方法；非定常流動(dòng)采用非線性諧波平衡方法進(jìn)行頻域模擬；并行計(jì)算采用OPENMPI。葉片顫振穩(wěn)定性使用能量法判據(jù)[21]，即在葉片的一個(gè)振動(dòng)周期內(nèi)，葉片振動(dòng)系統(tǒng)從外界獲得的能量與系統(tǒng)阻尼消耗的能量的正負(fù)決定了葉片是否顫振[22]。累積功W的計(jì)算公式為：

式中：AX為葉片振動(dòng)位移，Af1為葉片所受非定常力，φ1為非定常力相位角。

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

該大涵道比風(fēng)扇性能試驗(yàn)件，在85%轉(zhuǎn)速外涵逼喘過(guò)程中，接近喘振邊界時(shí)風(fēng)扇葉片振動(dòng)超限，圖1為監(jiān)控到的葉尖振幅曲線。可見(jiàn)11號(hào)葉片葉尖振幅最大值達(dá)到13.61 mm；葉片從32.875 s起振，在35.671 s葉尖振幅達(dá)到最大值，歷時(shí)2.796 s。以該發(fā)生顫振的風(fēng)扇葉片為研究對(duì)象，給定85%轉(zhuǎn)速近喘點(diǎn)的工況，分別基于經(jīng)驗(yàn)法和數(shù)值模擬方法，評(píng)估葉片的顫振穩(wěn)定性，并與試驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證。

圖1 85%轉(zhuǎn)速外涵逼喘過(guò)程葉片葉尖振幅曲線Fig.1 Blade tip vibration amplitude curve during surge at 85%speed

3.1 經(jīng)驗(yàn)法計(jì)算過(guò)程及結(jié)果分析

3.1.1 計(jì)算過(guò)程

根據(jù)85%轉(zhuǎn)速近喘點(diǎn)工況，首先完成定常氣動(dòng)計(jì)算。三維流場(chǎng)計(jì)算采用商業(yè)軟件NUMECA的FINE/TURBO求解。計(jì)算設(shè)置如下：采用S-A模型；計(jì)算域進(jìn)口按照試驗(yàn)時(shí)的真實(shí)情況給定總溫、總壓和進(jìn)口氣流角；風(fēng)扇內(nèi)外涵分別給定出口平均靜壓；固壁為絕熱、無(wú)滑移邊界條件；轉(zhuǎn)靜交界面選擇為周向守恒型[23-25]。計(jì)算網(wǎng)格如圖2所示，風(fēng)扇和外涵出口導(dǎo)流葉片(OGV)的網(wǎng)格拓?fù)浣Y(jié)構(gòu)為默認(rèn)的O4H，內(nèi)涵通道中葉片的網(wǎng)格拓?fù)浣Y(jié)構(gòu)為HOH。離開(kāi)葉片表面第一層網(wǎng)格的距離為5 10-6m，y+值保持在10以下。總網(wǎng)格量約300萬(wàn)。

圖3和圖4示出了計(jì)算所得外涵特性線與試驗(yàn)結(jié)果對(duì)比?？梢钥吹剑簤罕忍匦缘挠?jì)算結(jié)果與試驗(yàn)結(jié)果較接近；效率特性方面，由于設(shè)計(jì)壓比較低和受探針測(cè)量精度限制，外涵所錄取的溫升效率較計(jì)算結(jié)果偏高，但趨勢(shì)基本一致。取特性曲線中近喘點(diǎn)(圖4中藍(lán)色框所示)，將外涵OGV前的總壓比和總溫比徑向分布與級(jí)間探針測(cè)量結(jié)果進(jìn)行對(duì)比，如圖5和6所示。考慮到溫度測(cè)量誤差，可認(rèn)為近喘點(diǎn)的流場(chǎng)計(jì)算較為準(zhǔn)確。

圖3 風(fēng)扇外涵效率特性Fig.3 Adiabatic efficiency characteristic lines of fan bypass

圖4 風(fēng)扇外涵壓比特性Fig.4 Total pressure ratio characteristic lines of fan bypass

圖5 外涵OGV前總壓比Fig.5 Inlet total pressure ratio of bypass OGV

圖6 外涵OGV前總溫比Fig.6 Inlet total temperature ratio of bypass OGV

其次，完成風(fēng)扇葉片強(qiáng)度計(jì)算。風(fēng)扇轉(zhuǎn)子為整體葉盤(pán)結(jié)構(gòu)，共有18片葉片，其材料為YZ-TC4鍛件，表1給出了材料性能參數(shù)。葉盤(pán)前緣與帽罩用花鍵連接，葉盤(pán)安裝邊與風(fēng)扇軸用螺栓連接，輪盤(pán)后端面為自由狀態(tài)。網(wǎng)格單元采用六面體，單元類型為二階solid186，單元數(shù)為114 336，節(jié)點(diǎn)數(shù)為459 790，建立的有限元模型如圖7所示。計(jì)算時(shí)，施加的邊界條件及載荷為：①對(duì)風(fēng)扇盤(pán)軸連接法蘭(面3)施加軸向和周向位移約束，風(fēng)扇盤(pán)周期對(duì)稱面(面1、面2)施加周期對(duì)稱約束；②對(duì)風(fēng)扇盤(pán)腔表面施加均布的氣動(dòng)載荷，以模擬腔壓對(duì)風(fēng)扇盤(pán)的影響；③對(duì)風(fēng)扇整體葉盤(pán)施加離心力載荷，對(duì)葉身和葉盤(pán)流道表面施加氣動(dòng)力載荷，對(duì)葉片表面進(jìn)行氣動(dòng)力插值；④設(shè)置分析類型為大變形靜態(tài)分析進(jìn)行計(jì)算，并在靜強(qiáng)度分析基礎(chǔ)上進(jìn)行考慮應(yīng)力剛化和旋轉(zhuǎn)軟化效應(yīng)的模態(tài)分析；⑤在笛卡爾坐標(biāo)系下，提取計(jì)算工況下葉尖前、后緣節(jié)點(diǎn)前3階模態(tài)的位移量、比例系數(shù)和頻率，以及葉尖偏移和扭轉(zhuǎn)變形計(jì)算結(jié)果。

表1 風(fēng)扇整體葉盤(pán)材料性能參數(shù)Table 1 Material properties of fan blisk

圖7 風(fēng)扇轉(zhuǎn)子有限元模型Fig.7 Finite element model of fan rotor

表2示出了風(fēng)扇葉片前3階模態(tài)下的固有頻率。圖8為風(fēng)扇整體葉盤(pán)葉身在85%轉(zhuǎn)速外涵喘點(diǎn)狀態(tài)下的徑向變形、軸向變形和總變形分布。由圖可知，最大徑向變形、軸向變形和最大總變形均出現(xiàn)在葉尖，分別為1.03 mm、-6.64 mm(向后)和8.71 mm。圖9為風(fēng)扇整體葉盤(pán)葉身的徑向應(yīng)力、第一主應(yīng)力和等效應(yīng)力分布云圖。可看出，最大第一主應(yīng)力出現(xiàn)在葉盆中部區(qū)域，為425.19 MPa。

表2 風(fēng)扇葉片固有頻率Table 2 Natural frequency of fan blades

最后，比較顫振臨界攻角和實(shí)際攻角。具體步驟為：

圖8 葉身變形云圖Fig.8 Deformation cloud diagram of blade

圖9 葉身應(yīng)力分布云圖Fig.9 Stress distribution cloud diagram of blade

(2)計(jì)算1～3階模態(tài)振型。模態(tài)振型的定義為Flap/(φ·B′)。其中：Flap為葉尖彎曲分量，即變形前后葉尖周向偏移量；φ為葉尖扭轉(zhuǎn)角，即變形前后弦線夾角；B′為風(fēng)扇葉片葉尖弦長(zhǎng)。

(3)以折合速度為橫坐標(biāo)，模態(tài)振型為縱坐標(biāo)，在經(jīng)驗(yàn)曲線圖上查出1～3階的顫振臨界攻角。

(4)計(jì)算風(fēng)扇葉片75%葉高處的實(shí)際攻角i。i=相對(duì)氣流角-金屬角，相對(duì)氣流角來(lái)自三維計(jì)算結(jié)果，金屬角為葉片造型結(jié)果。

(5)比較實(shí)際攻角與1～3階臨界攻角大小，若任意1階的實(shí)際攻角大于臨界攻角，則葉片發(fā)生顫振，反之則穩(wěn)定。

3.1.2 結(jié)果分析

表3示出了風(fēng)扇葉片的經(jīng)驗(yàn)法顫振評(píng)估結(jié)果?？煽闯?階模態(tài)下，風(fēng)扇葉片的實(shí)際攻角大于顫振臨界攻角，評(píng)估結(jié)果為葉片發(fā)生顫振，與試驗(yàn)現(xiàn)象一致。

表3 經(jīng)驗(yàn)法顫振評(píng)估結(jié)果Table 3 Flutter evaluation results of empirical method

3.2 數(shù)值模擬方法計(jì)算過(guò)程及結(jié)果分析

3.2.1 計(jì)算過(guò)程

顫振數(shù)值模擬采用Turbo3D軟件。該軟件使用頻域法進(jìn)行顫振分析，是一種流固解耦的能量法分析方法。主要計(jì)算過(guò)程如下：

(1)使用網(wǎng)格轉(zhuǎn)換程序?qū)umeca生成的網(wǎng)格文件轉(zhuǎn)換為T(mén)urbo3D程序可讀的網(wǎng)格文件，設(shè)定輸入邊界條件(進(jìn)口總壓、總溫、氣流角、轉(zhuǎn)速、摻混面模型、湍流模型、CFL數(shù)、人工粘性系數(shù)、壁面函數(shù)、多重網(wǎng)格層數(shù)等)、初場(chǎng)文件(各葉片排的壓力、軸向速度、切向速度、徑向速度、密度等)和交界面設(shè)置文件(各交界面位置的徑向網(wǎng)格對(duì)接點(diǎn)數(shù)，出口邊界條件等)，完成85%轉(zhuǎn)速近喘點(diǎn)定常計(jì)算。

(2)將強(qiáng)度計(jì)算得到的風(fēng)扇模態(tài)變形量數(shù)據(jù)連續(xù)插值到CFD網(wǎng)格上，并設(shè)置葉間相位角。插值方法為二維線性插值，兼容強(qiáng)度有限元網(wǎng)格中包含的非結(jié)構(gòu)化六面體網(wǎng)格。對(duì)于有限元網(wǎng)格上存在的三棱柱網(wǎng)格、四面體網(wǎng)格等非結(jié)構(gòu)化六面體單元，均按照20節(jié)點(diǎn)輸出的方式，確定單元的外表面。插值原理是使用CFD網(wǎng)格表面的四邊形節(jié)點(diǎn)的三個(gè)點(diǎn)構(gòu)成單獨(dú)的三角形面，選取離此三角形質(zhì)心最近的有限元網(wǎng)格上的點(diǎn)賦值給此三角形的三個(gè)節(jié)點(diǎn)，以提高插值的準(zhǔn)確性。

(3)以定常計(jì)算和模態(tài)插值結(jié)果為輸入，采用頻域諧波方法進(jìn)行非定常計(jì)算。假設(shè)非定常頻率是葉片振動(dòng)自然頻率的線性疊加，葉片按預(yù)先設(shè)定方式、振幅振動(dòng)。葉片間相位角計(jì)算范圍為-160.0h～160.0h，對(duì)應(yīng)的節(jié)徑為-8～8。

(4)計(jì)算1～3階模態(tài)下各葉間相位角的累積功，若各模態(tài)各葉間相位角的累積功都小于0則葉片穩(wěn)定，若在某階模態(tài)某葉間相位角的累積功大于0則存在顫振風(fēng)險(xiǎn)。

3.2.2 結(jié)果分析

圖10為85%轉(zhuǎn)速外涵近喘點(diǎn)前3階振型對(duì)應(yīng)不同葉片間相位角的累積功隨節(jié)徑的變化特性。如圖中紅圈所示，對(duì)應(yīng)1階振型-1、-2節(jié)徑的累積功大于0，其中最大累積功出現(xiàn)在-1節(jié)徑，對(duì)應(yīng)為-20.0h葉片間相位角。根據(jù)能量法判據(jù)，該風(fēng)扇葉片存在顫振風(fēng)險(xiǎn)，與試驗(yàn)現(xiàn)象一致。

圖11為最不穩(wěn)定點(diǎn)節(jié)徑-1時(shí)葉片表面的氣動(dòng)功(氣動(dòng)功沿表面積分即為累積功)分布，上圖為吸力面，下圖為壓力面?？梢钥吹?，在吸力面葉尖激波位置處出現(xiàn)氣動(dòng)功峰值。通過(guò)減弱葉尖激波強(qiáng)度以減小當(dāng)?shù)氐臍鈩?dòng)功，從而減小1階的累積功，抑制顫振風(fēng)險(xiǎn)。

圖10 前3階振型對(duì)應(yīng)不同節(jié)徑的累積功Fig.10 Worksum of different nodal diameter for the first three mode shapes

圖11 節(jié)徑-1時(shí)風(fēng)扇葉片表面氣動(dòng)功分布Fig.11 Fan blade surface worksum distribution with the nodal diameter equals to-1

4 結(jié)論

基于經(jīng)驗(yàn)法和數(shù)值模擬方法對(duì)發(fā)生顫振的風(fēng)扇葉片進(jìn)行了顫振預(yù)測(cè)計(jì)算分析，得到如下結(jié)論：

(1)三參數(shù)法和數(shù)值模擬方法的分析結(jié)果均表明葉片發(fā)生顫振，且與試驗(yàn)結(jié)果一致。兩種方法的計(jì)算量均在可承受范圍內(nèi)，工程上具備實(shí)用價(jià)值。

(2)根據(jù)三參數(shù)法的預(yù)測(cè)結(jié)果，可通過(guò)以下手段降低葉片顫振風(fēng)險(xiǎn)：減小近外涵喘點(diǎn)風(fēng)扇葉片中上部葉高處的氣動(dòng)攻角，增大葉片1階頻率，增大葉片75%葉高處弦長(zhǎng)并減小100%葉高處弦長(zhǎng)，減小葉片尖部扭轉(zhuǎn)角。

(3)數(shù)值模擬方法的分析結(jié)果表明，葉片中上部激波位置處氣動(dòng)功出現(xiàn)峰值。通過(guò)減小喘點(diǎn)風(fēng)扇葉片攻角，改善流動(dòng)狀態(tài)，以及減弱葉尖激波，可以有效降低顫振風(fēng)險(xiǎn)。