寬弦風扇葉片顫振預測的工程研究

章嘉麟，丁建國

(中國航發商用航空發動機有限責任公司，上海201108)

1 引言

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

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

某大涵道比風扇葉片在設計中采用了寬弦和復合彎掠等技術，在性能考核試驗過程中風扇葉片發生了顫振現象。本文分別基于經驗法和數值模擬方法，以發生顫振的風扇葉片為對象進行顫振評估。以期通過對其的分析研究，掌握該類葉片的顫振機理；同時，將計算結果與試驗結果進行對比，以驗證顫振預測方法的準確性。

2 風扇葉片顫振預估方法

2.1 經驗法

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

2.2 數值模擬方法

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

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

3 計算結果及分析

該大涵道比風扇性能試驗件，在85%轉速外涵逼喘過程中，接近喘振邊界時風扇葉片振動超限，圖1為監控到的葉尖振幅曲線。可見11號葉片葉尖振幅最大值達到13.61 mm；葉片從32.875 s起振，在35.671 s葉尖振幅達到最大值，歷時2.796 s。以該發生顫振的風扇葉片為研究對象，給定85%轉速近喘點的工況，分別基于經驗法和數值模擬方法，評估葉片的顫振穩定性，并與試驗結果進行對比驗證。

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

3.1 經驗法計算過程及結果分析

3.1.1 計算過程

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

圖3和圖4示出了計算所得外涵特性線與試驗結果對比。可以看到：壓比特性的計算結果與試驗結果較接近；效率特性方面，由于設計壓比較低和受探針測量精度限制，外涵所錄取的溫升效率較計算結果偏高，但趨勢基本一致。取特性曲線中近喘點(圖4中藍色框所示)，將外涵OGV前的總壓比和總溫比徑向分布與級間探針測量結果進行對比，如圖5和6所示。考慮到溫度測量誤差，可認為近喘點的流場計算較為準確。

圖3 風扇外涵效率特性Fig.3 Adiabatic efficiency characteristic lines of fan bypass

圖4 風扇外涵壓比特性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

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

表1 風扇整體葉盤材料性能參數Table 1 Material properties of fan blisk

圖7 風扇轉子有限元模型Fig.7 Finite element model of fan rotor

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

表2 風扇葉片固有頻率Table 2 Natural frequency of fan blades

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

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

圖9 葉身應力分布云圖Fig.9 Stress distribution cloud diagram of blade

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

(3)以折合速度為橫坐標，模態振型為縱坐標，在經驗曲線圖上查出1～3階的顫振臨界攻角。

(4)計算風扇葉片75%葉高處的實際攻角i。i=相對氣流角-金屬角，相對氣流角來自三維計算結果，金屬角為葉片造型結果。

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

3.1.2 結果分析

表3示出了風扇葉片的經驗法顫振評估結果。可看出1階模態下，風扇葉片的實際攻角大于顫振臨界攻角，評估結果為葉片發生顫振，與試驗現象一致。

表3 經驗法顫振評估結果Table 3 Flutter evaluation results of empirical method

3.2 數值模擬方法計算過程及結果分析

3.2.1 計算過程

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

(1)使用網格轉換程序將Numeca生成的網格文件轉換為Turbo3D程序可讀的網格文件，設定輸入邊界條件(進口總壓、總溫、氣流角、轉速、摻混面模型、湍流模型、CFL數、人工粘性系數、壁面函數、多重網格層數等)、初場文件(各葉片排的壓力、軸向速度、切向速度、徑向速度、密度等)和交界面設置文件(各交界面位置的徑向網格對接點數，出口邊界條件等)，完成85%轉速近喘點定常計算。

(2)將強度計算得到的風扇模態變形量數據連續插值到CFD網格上，并設置葉間相位角。插值方法為二維線性插值，兼容強度有限元網格中包含的非結構化六面體網格。對于有限元網格上存在的三棱柱網格、四面體網格等非結構化六面體單元，均按照20節點輸出的方式，確定單元的外表面。插值原理是使用CFD網格表面的四邊形節點的三個點構成單獨的三角形面，選取離此三角形質心最近的有限元網格上的點賦值給此三角形的三個節點，以提高插值的準確性。

(3)以定常計算和模態插值結果為輸入，采用頻域諧波方法進行非定常計算。假設非定常頻率是葉片振動自然頻率的線性疊加，葉片按預先設定方式、振幅振動。葉片間相位角計算范圍為-160.0h～160.0h，對應的節徑為-8～8。

(4)計算1～3階模態下各葉間相位角的累積功，若各模態各葉間相位角的累積功都小于0則葉片穩定，若在某階模態某葉間相位角的累積功大于0則存在顫振風險。

3.2.2 結果分析

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

圖11為最不穩定點節徑-1時葉片表面的氣動功(氣動功沿表面積分即為累積功)分布，上圖為吸力面，下圖為壓力面。可以看到，在吸力面葉尖激波位置處出現氣動功峰值。通過減弱葉尖激波強度以減小當地的氣動功，從而減小1階的累積功，抑制顫振風險。

圖10 前3階振型對應不同節徑的累積功Fig.10 Worksum of different nodal diameter for the first three mode shapes

圖11 節徑-1時風扇葉片表面氣動功分布Fig.11 Fan blade surface worksum distribution with the nodal diameter equals to-1

4 結論

基于經驗法和數值模擬方法對發生顫振的風扇葉片進行了顫振預測計算分析，得到如下結論：

(1)三參數法和數值模擬方法的分析結果均表明葉片發生顫振，且與試驗結果一致。兩種方法的計算量均在可承受范圍內，工程上具備實用價值。

(2)根據三參數法的預測結果，可通過以下手段降低葉片顫振風險：減小近外涵喘點風扇葉片中上部葉高處的氣動攻角，增大葉片1階頻率，增大葉片75%葉高處弦長并減小100%葉高處弦長，減小葉片尖部扭轉角。

(3)數值模擬方法的分析結果表明，葉片中上部激波位置處氣動功出現峰值。通過減小喘點風扇葉片攻角，改善流動狀態，以及減弱葉尖激波，可以有效降低顫振風險。