轉子葉片氣彈穩定性與強迫響應分析

韓 樂，王延榮，2

（1.北京航空航天大學能源與動力工程學院，北京100083；2.北京航空航天大學江西研究院，南昌330096）

0 引言

葉輪機械中的氣動彈性現象導致葉片高循環疲勞失效，其中氣動彈性失穩和強迫振動響應是2個主要方面。隨著航空發動機性能的大幅提升，葉片氣動負荷增大，結構更加緊湊，顫振和強迫振動響應問題愈加突出。

氣動彈性失穩（顫振）是彈性系統在氣流中的耦合自激振動，屬于氣動彈性力學中的動穩定性問題。當轉子葉片發生顫振時，流體激勵源于葉片振動本身，在激勵與振動相互作用下短時間內振幅急劇增大，進而造成葉片破壞。目前常用的氣動彈性穩定性分析方法包括能量法[1-2]和特征值法[3-4]。Erdos等[5]和He[6]將相位延遲（phase lag）法應用于氣彈穩定性預測，使能量法的計算速度大幅提高；通過結合Hanamur等[7]研究的影響系數法使特征值法的分析效率顯著提高；趙瑞勇等[8]和Fu等[9]利用上述方法對風扇葉片的氣動彈性穩定性也進行了分析。

轉子葉片強迫振動響應通常是由前排葉片尾跡或下游葉片的勢擾動引起的，激勵頻率與轉子轉速和引起激勵的結構特征直接相關，且不受葉片振動的影響?；诮怦罘椒▉矸治鋈~片強迫振動響應。Kielb等[10]將強迫響應問題的求解分為3個階段：識別來自上下游流動的不均勻性；計算葉片非定常氣動力；獲取葉片振動響應。文獻[11-12]中將最后1個階段又分為非線性結構建模和柔性葉片的氣動彈性響應求解。

此外，還可以采取耦合法分析葉片氣動彈性問題，包括強耦合法和弱耦合法。強耦合法是將流體力學控制方程和結構動力學方程統一求解，由于2個計算域的時間和空間尺寸相差偏大，此方法求解難度較大，且計算量也極大；弱耦合法則利用每個子系統原有的求解方法，通過固體域與流體域相關邊界條件的交替迭代更新數據，實現對耦合作用的模擬，雖然耗費很多計算資源，但部分學者仍采用該方法分析氣動彈性問題。Sadeghi等[13]采用時域法分析了NASA67轉子葉片扭轉模態下的氣動彈性，發現在某折合頻率下葉片會出現失穩現象；徐可寧[14]自主編寫了有限元程序，結合計算流體力學軟件分析了錯頻葉盤振動的局部化現象；Espinal[15]利用弱耦合方法分析了某軸流跨聲速葉片的氣動彈性問題，驗證了相位延遲法的可靠性。在多數情況下，由于考慮到計算能力，利用弱耦合法在多級環境中分析氣動彈性問題通常仍是難以開展的。

針對上述問題，本文分別以壓氣機轉子葉片和高壓渦輪轉子葉片模型為例，對氣動彈性穩定性和強迫振動響應進行了分析，給出了壓氣機轉子葉片氣動彈性穩定性特征和高壓渦輪轉子葉片強迫振動響應特征。

1 分析方法

1.1 氣動彈性穩定性分析方法

目前多采用能量法和特征值法進行氣動彈性穩定性預測，由于耦合法計算量過大，難以應用于工程，一般多用于校核所發展的新方法。下面主要介紹顫振預測的能量法和特征值法。

1.1.1 能量法

假設顫振以葉片某一階固有振型出現，通過計算1個振動周期內葉片與流場間的能量交換來預測顫振。如果1個振動周期內周圍氣流對葉片所作的非定常氣動功總和Waero為正，則葉片發生氣動彈性失穩，如果考慮機械阻尼，則機械阻尼耗散功與氣動功之和決定了系統的穩定性。

基于能量法預測顫振的主要流程如下：

（1）計算葉片的固有模態，分析中可考慮穩態氣動力和離心力的影響；

（2）將所關注的葉片固有模態位移映射到流場網格葉片表面上；

（3）指定流場中的葉片按照所得到固有模態振動，進行非定常流場分析；

（4）計算收斂后，提取葉片表面的力和位移，計算1個振動周期內的非定常氣動功，進而得到氣動阻尼；

（5）根據氣動阻尼的正負判斷系統氣彈穩定性。

在得到葉片表面節點集中力和節點位移后，非定常氣動功可表示為

式中：i為節點編號；nnode和ntstep分別為葉片表面流體節點總數和1個振動周期內的時間步數；j為時間步編號；F為節點集中力，N；D為節點相對于葉片平衡位置的位移，m。

本文以Moffatt等[16]提出的氣動模態阻尼比（Aerodynamic Modal Damping Ratio，AMDR）來表征氣動阻尼的大小，基于等效黏性阻尼的概念，氣動模態阻尼比ζaero可表示為

式中：q為指定模態的正則化幅值；ω為給定模態的固有振動圓頻率，rad/s。

基于等效黏性阻尼的形式，可直接與機械阻尼等代數相加得到總阻尼。

在相位延遲邊界條件下，通過存儲1個周期內邊界面上的所有信息實現對模型的簡化[5]。再利用時間范疇的傅里葉級數分解見式（3），進一步避免存儲周期邊界上的全部信息，僅需存儲傅里葉系數An就可在后處理中通過這些系數求出任意時間點的變量信息[6]。

當考慮非零葉間相位角時，需要指定相鄰葉片按照一定的時間相位差振動，相鄰2片葉片振動時間差為

式中：σ為葉間相位角。

1.1.2 特征值法

對于調諧葉盤結構系統，不考慮機械阻尼時，物理坐標系的運動方程為

式中：u為位移向量；M、K和f分別為質量矩陣、剛度矩陣和葉片振動所引起的非定常氣動力向量。

假設葉片按照簡諧形式振蕩u=Ueλt，并將位移幅值向量U寫成模態疊加的形式，式（5）可表示為

式中：Φ為模態矩陣；q為模態位移向量；Ω為對角矩陣。

當Φ中只包含某一階模態時，Ω的表達式為

式中：Ω中各元素為葉片各階次頻率的平方。

對于調諧系統，式（7）的各對角線元素相等，每個元素代表相應葉片的固有頻率，通過修改對角線元素便于引入錯頻。

當結合影響系數法時，式（6）中模態坐標系下的氣動力向量可表示為

式中：A為模態坐標系下的氣動力影響系數矩陣，代表了各葉片對參考葉片的影響，對調諧系統，A為循環矩陣，有如下形式

當僅考慮1階模態時，式（9）中的各分塊矩陣變為常數。將式（8）帶入式（6）中，可得到如下的特征值

通過求解式（10）所表示的特征值，便可得到復特征值λ，其中實部代表阻尼，虛部代表頻率。系統發生顫振（氣彈失穩）的依據為：存在特征值滿足

當Φ為質量歸一模態振型時，氣動阻尼表示為

1.2 強迫振動響應分析方法

在強迫振動響應分析中，當葉片發生振動時，其表面的壓力擾動主要源于上下游氣流的周期性激勵，而與葉片本身的變形關系不大。因此，先計算1個葉片掃掠周期內剛性葉片周圍的非定常流場，再將其作為外部激勵施加于葉片的有限元模型上。雖然這種假設對計算精度的影響不大，但可以極大地提高計算效率，適用于工程分析，其分析流程如圖1所示。

圖1 強迫振動響應分析方法及流程

主要步驟為：

（1）建立葉片流體域計算模型，進行定常模擬；

（2）建立葉片有限元模型，進行模態分析；

（3）通過葉片Campbell圖確定共振轉速，明確重點關注的模態階次；

（4）結合模態分析結果，找到所關注模態階次振動應力較大的位置作為瞬態計算的監測點；

（5）分析在轉靜干涉作用下的非定常流場特征，并提取1個周期內葉片表面非定常氣動力；

（6）計算在共振狀態下葉片強迫響應特征，當振動穩定后終止計算；

（7）給出各監測點沿應力最大方向的動應力響應，利用傅里葉變換分析動應力特征。

2 壓氣機轉子葉片氣動彈性穩定性分析

2.1 計算模型

以高壓壓氣機模型中間級的轉子葉片為研究對象，相關設計和材料參數見表1。建立孤立轉子葉片模型并計算在級條件下的氣動阻尼，各參數包括進口總溫、總壓和氣流角、出口靜壓，見表2，其子午面流道如圖2所示。

表1 轉子葉片的無量綱設計參數和材料參數

表2 級條件下進出口平均參數

圖2 轉子葉片及上下游流道

2.2 定常流場分析

建立轉子葉片的流場計算模型，根據葉片在多級環境中的位置截取計算域，如圖3所示。計算網格采用O4H拓撲結構，沿流向布置69層網格；周向布置57層網格（包括主流區33層和葉片周圍O網格13層）；徑向布置53層（間隙內13層）；單通道網格約為26.1萬。湍流模型選擇標準k-ε雙方程模型，第1層網格距壁面距離為1×10-5m。進口給定總溫、總壓和速度方向，出口給定靜壓條件，相關條件由多級計算結果提取。計算結果見表3。從表中可見，轉子的氣動效率較高。在不同條件下葉尖流場如圖4所示。從圖中可見，此時為在亞聲速條件下葉片全展向工作。

圖3 轉子葉片定常流場分析計算域

表3 定常流場計算結果

圖4 典型葉高下的流場

進一步給出了多級環境和孤立轉子模擬（級條件下）50%葉高截面的壓力分布，如圖5所示。從圖中可見，壓力面吻合較好，吸力面存在一定差異。其原因是在多級環境下，流動偏離設計狀態并逐級放大所致。

圖5 典型葉高葉片表面壓力分布

2.3 模態分析

建立轉子葉片有限元模型及冷熱態葉型，如圖6所示。從圖中可見，采用8節點6面體單元，分別沿徑向和弦向各布置41層和29層網格節點，厚度方向布置3層網格節點，節點和單元總數分別為3567和2240。在設計轉速下的靜力分析結果如圖7所示。由離心力作用引起的最大靜變形約為0.55 mm。隨后進行含預應力的模態分析，重點考查葉尖前緣附近應力較大的振型(第4和7~9階模態)，其振型和模態應力分布如圖8所示。

圖6 轉子葉片有限元模型和冷熱態葉型對比

圖7 葉片在離心力作用下的靜變形

圖8 轉子葉片振型（左）和模態應力（右）

2.4 振蕩流場和氣動阻尼分析

分別采用能量法和特征值法對轉子葉片的氣動阻尼進行分析。非定常流場的多通道計算域如圖9所示。在計算過程中指定葉片按照不同模態振動，分析各節徑下的氣動阻尼。在設置葉片振幅時，按照50 MPa振動應力對應的振幅進行折算，得到各階振幅和AF值，見表4。AF=2πfA，代表了葉片的振動速度。以定常計算結果作為非定常計算的初始流場，在1個振動周期內劃分60個非定常時間步，當流場參數達到穩態振蕩時停止計算。分別通過式（2）和式（12）計算氣動阻尼。

圖9 非定常流場計算域模型

表4 在設計轉速下葉片的動頻、振幅及AF值

以第8階模態為例，葉片表面非定常氣動力幅值分布如圖10所示。從圖中可見，對高階模態而言，參考葉片振動引起的非定常氣動力主要集中在中間2~3個通道內。

圖10 第8階模態下葉片表面非定常氣動力幅值

在第8階模態下氣動阻尼隨節徑的變化如圖11所示。從圖中可見，通過2種方法得到的氣動阻尼較為一致，特別是最小氣動阻尼相互吻合。

圖11 第8階模態下氣動阻尼

在各階模態下最小氣動阻尼及其對應的節徑見表5。從表中可見，高階模態氣動阻尼反而更小，若附近存在激勵成分，可能引起較顯著的振動應力。結合圖8中模態分析結果，一旦激起第7、8和9階模態，在葉尖近前緣和近尾緣區域可能出現較大的振動應力。在第4階模態下，氣動阻尼隨模態頻率增加而減小。因此對該轉子葉片而言，需要重點關注第2階弦向彎曲（第7階模態）和第8、9階模態的氣動阻尼。

表5 不同模態下的最小氣動阻尼及其節徑數

同樣以第8階模態為例，最小氣動阻尼所在節徑下葉片表面的氣動功密度分布如圖12所示。從圖中可見，葉身表面大部分區域氣動功幅值接近0，只是在葉尖近前緣區域有1個明顯的負功區。

圖12 第8階模態60節徑下氣動功密度

3 強迫振動響應分析

3.1 計算模型

對渦輪模型的高壓轉子葉片進行分析，其流場計算域如圖13所示。從圖中可見，模型包括高低壓渦輪轉子、上游進口導葉和下游支板。模型中的設計參數和設計工況下的進出口邊界條件見表6。4排葉片不滿足約化條件，所以采用全環模型進行非定常計算。在計算域中各排葉片的網格見表7。每排葉片各通道的網格拓撲和網格數保持一致。在非定常計算中給定時間步，確保高、低壓渦輪葉片均經過整數時間步后通過彼此1個葉片柵距。

圖13 渦輪葉片CFD計算模型

表6 渦輪轉子葉片無量綱設計參數和進出口氣動參數

表7 各排葉片網格數 萬

3.2 定常和非定常流場分析

在葉片壓力面和吸力面的90%、70%和50%葉高大約10%、50%和90%弦長處分別設置壓力監測點，如圖14所示。為保證信號的周期性，提取壓力信號完整周期進行傅里葉分析。以弦向S1位置的3個監測點為例，壓力隨時間的變化曲線如圖15所示，此時流場中存在多個擾動成分。

圖14 葉片表面監測點

圖15 監測點壓力隨時間的變化

以吸力面S1和S3監測位置為例，對葉片表面監測點的壓力數據進行時頻分析，得到頻域信號如圖16所示。從圖中可見，主要激勵為上游進口導葉（12192 Hz）及其倍頻成份和在尾緣附近存在低壓渦輪引起的激勵（37060 Hz）。

圖16 監測點非定常氣動力頻譜分析

流場90%和50%葉高截面的相對馬赫數分布結果如圖17所示。從圖中可見，在高壓渦輪轉子葉片尾緣存在明顯的激波結構，導致低壓葉片帶來的勢擾動難以向上游傳播，因此高壓渦輪轉子葉片主要感受來自進口導葉帶來的擾動，幅值可達60~70 kPa。

圖17 90%和50%葉高截面相對馬赫數分布

3.3 模態分析

基于流場模型建立有限元模型，如圖18所示。通過葉片冷熱態轉換獲得冷態葉型[17]。采用8節點6面體單元，其中節點數為3450，單元數為2552，K447A為葉片材料。

圖18 高壓渦輪葉片有限元模型

對高壓渦輪葉片進行靜力分析，在溫度場和離心載荷（100%轉速）作用下高壓渦輪的靜變形如圖19所示。從圖中可見，溫度場形成較大的徑向變形，離心力和氣動力則使葉片更加“展開”，形成較大的周向變形。在幾種載荷影響下，葉片最大變形均出現在葉尖尾緣區域。

圖19 不同載荷下高壓渦輪葉片靜變形

進行含預應力的模態分析，結合前排導葉和低壓渦輪的激勵得到高壓渦輪轉子葉片的Campbell圖，如圖20所示。在設計轉速下，進口導葉通過頻率位于第4、5階模態之間，低壓渦輪帶來的激勵靠近第16階模態。因此，重點關注高壓渦輪葉片在上述幾階模態和激起振動所需能量較低的前3階模態，模態頻率見表8。從表中可見，離心力使葉片剛性增強，頻率加大；溫度使葉片固有頻率降低。

圖20 高壓渦輪轉子葉片Campbell圖

表8 高壓渦輪葉片主要階模態固有頻率（Hz）和阻尼

根據模態應力較大點確定強迫振動響應分析的監測點，如圖21所示。各階模態應力和振型分布如圖22所示。模擬貼應變片的方式并結合模態分析結果，根據式（13）確定高壓渦輪監測點的監測方向，見表9。

圖21 高壓渦輪葉片振動應力監測點位置

圖22 主要模態階次振型（上）和模態應力分布（下）

表9 高壓渦輪監測點的應力方向

式中：nx、ny和nz為最大應力方向余弦。

由應力分量確定最大振動應力為

3.4 強迫振動響應分析

3.4.1 模態阻尼在瞬態響應計算中阻尼以瑞利阻尼形式給出

式中：系數α和β由任意2階模態確定。假設該2階模態的阻尼比分別為ζi和ζj，頻率值分別為ωi和ωj，則

確定系數α和β后，得到任意階模態的阻尼比

由于高壓渦輪為整體葉盤結構，其結構阻尼較小，參考文獻[14]，取葉片第1、2階模態阻尼比為0.03%，此時阻尼系數α和β分別為9.55356×10-9和8.79657×10-9，由式（17）得到主要階模態阻尼比（表8）。

3.4.2 監測點振動應力

在瞬態動力學分析中，為保證得到穩態響應解，在計算過程中首先通過設置較大的阻尼系數使瞬態響應部分在短時間內快速衰減。在計算收斂后，某監測點時域結果如圖23所示。從圖中可見，隨著時間推進，葉片響應逐漸趨于較為穩定的周期振動，提取最后1個周期的結果用于振動應力分析。

圖23 瞬態響應計算監測點周向振動應力的收斂情況

以節點2879和2303為例，提取高壓渦輪轉子葉片監測點穩態振動應力并進行快速傅里葉變換，結果如圖24所示。對高壓渦輪葉片而言，振動主要由上游進口導葉引起，在低階模態上振動較為顯著。雖然激勵頻率更接近第5階模態固有頻率，但此時尚未達到共振狀態，最大振動應力出現在第3階模態頻率下，約為12 MPa（阻尼為0.033%）。

圖24 高壓渦輪轉子葉片監測點的振動應力

3.4.3 危險模態共振狀態分析

進一步分析共振狀態下的響應特征。由于流場計算過于耗時，且第5階模態共振轉速較為靠近計算轉速，工況偏差較小，故利用頻率縮放技術對共振狀態進行快速評估。由模態分析結果可知，第5階模態危險點為圖22中節點2225。對轉速進行縮放，使激勵與固有頻率接近，得到104%轉速下模態頻率、激勵頻率和阻尼，見表10。

表10 104%轉速下模態頻率、激勵頻率和阻尼

對激勵頻率進行縮放，研究高壓渦輪轉子葉片在上述轉速下的響應特征，近葉尖處監測點時域結果如圖25所示。從圖中可見，隨著時間的推進，振動表現為明顯的放大特征，表明此時激勵頻率與葉片固有頻率十分接近。提取最后1個周期結果用于振動應力分析。

圖25 104%轉速下節點2225周向振動應力時域結果

對監測點的振動響應進行傅里葉分析，高壓渦輪轉子葉片監測點的動態應力頻譜如圖26所示。在104%轉速下發生了第5階模態共振，此時危險點振動應力顯著增大，在0.14%阻尼下振動應力達到134 MPa。此外，對比響應和激勵的頻率（表10）可見，104%轉速下尚未達到最大共振狀態。

圖26 104%轉速下節點2225振動應力

4 結論

（1）給出了壓氣機轉子葉片第4（第1階弦向彎曲）、7（第2階弦向彎曲）、8（復合模態）和9階（復合模態）模態下的氣動阻尼，采用能量法和特征值法得到的計算結果有較好吻合性。高階模態（第7~9階模態）對應的氣動阻尼相對更小，葉身表面大部分區域非定常氣動功幅值接近0，只是在葉尖前尾緣區域存在較明顯的負功區，提供正阻尼。

（2）結合振型和應力分布分析，如果在激起葉片第7、8和9階模態振動，或者在這3階模態下發生氣彈失穩，葉尖前尾緣區域振動應力較大，易引起疲勞開裂。

（3）利用強迫振動響應分析方法對高壓渦輪葉片進行了分析，主要受上游進口導葉的影響，其帶來的非定常氣動力可達60~70 kPa，可能激起葉片第5階模態。當發生第5階模態共振時，振動應力可達134 MPa（阻尼0.14%）。