基于流固耦合方法的某風扇超聲速非失速顫振研究

張 帥，陶 冶，高 磊，張霞妹

(中國飛行試驗研究院 發(fā)動機所，陜西 西安 710089)

0 引 言

航空發(fā)動機的風扇/壓氣機轉子工作在高壓比、高轉速、高負荷的級環(huán)境中，實際流動中存在風扇的長葉片與流體間強烈的耦合作用，這會引起葉片的氣動彈性不穩(wěn)定現(xiàn)象，嚴重時會導致發(fā)動機非包容性破壞[1]。因此，對風扇/壓氣機葉片的氣動彈性響應預測和顫振產生及發(fā)展研究具有重要意義。發(fā)動機風扇/壓氣機中典型的顫振邊界包括：亞/跨聲速失速顫振、超聲速失速顫振、超聲速非失速顫振、堵塞顫振、彎扭耦合顫振、A100型顫振等類型[2]。

超聲速非失速顫振通常發(fā)生在大型跨聲速或超聲速風扇/壓氣機上，典型特征就是發(fā)生顫振的葉片葉尖馬赫數(shù)為超聲速，是當今最引人關注的一種顫振現(xiàn)象[3-5]。這種顫振之所以引起人們的強烈關注，不僅僅是因為渦扇發(fā)動機的應用越來越廣泛，更是因為這種顫振類型自身極其危險。D. G. Halliwell描述葉輪機械超聲速非失速顫振兩個典型特征：顫振邊界隨轉速變化而快速變化、一般為突發(fā)型顫振[6]。

國內外針對葉輪機械中的葉片氣動彈性問題的研究大多是通過流固耦合的分析方法來實現(xiàn)的。徐敏等[7]通過流固耦合的方法來預測葉輪機械的氣動彈性穩(wěn)定性。王征等[7]以簡化后的某風扇轉子為研究對象，對該風扇轉子進行全通道的雙向流固耦合數(shù)值模擬，研究其耦合前后氣動及結構性能的變化情況，并通過模態(tài)耦合等方法求解轉子葉片顫振問題，最終獲取其顫振邊界。

隨著民用渦扇發(fā)動機的涵道比及風扇負荷的不斷增加，風扇葉片的超聲速非失速顫振問題日益突出。筆者將針對高轉速下的風扇轉子設計轉速下近效率最高處的超聲速非失速顫振問題展開研究。

1 葉輪機械顫振分析方法

能量法是目前葉輪機械葉片顫振分析方法中應用最廣泛的一種，通過計算振動葉片表面非定常氣動力在一個振動周期內對葉片所做的積累功，進而以氣動力或氣動阻尼系數(shù)作為系統(tǒng)顫振分析的衡量標準，這種方法一般能夠較精確地對葉片的顫振現(xiàn)象進行預判。而通過流固耦合方法獲得葉片表面的非定常氣動力，不僅考慮到氣動力對葉片的影響，還考慮到葉片變形對轉子流場的影響，是一種比較接近實際流動的計算方法。

針對葉輪機械的超/跨聲速非失速顫振問題，采用基于流固耦合方法的能量累積法進行求解。目前，對于壓氣機葉片顫振問題的研究多以二維葉片算例為主，研究三維氣動彈性問題研究較少，且大多數(shù)不直接涉及顫振邊界問題研究。由于發(fā)動機中壓氣機內部流動具有強烈的三維效應，傳統(tǒng)方法選取三維模型的典型特征截面進行研究，這導致研究模型與三維流動實際狀況的誤差大，而按照三維問題進行研究就會更接近壓氣機內部的真實流動狀況。

1.1 葉輪簡化數(shù)學模型

由于流體計算軟件中不存在移相周期性邊界條件，在對葉片排的顫振問題進行數(shù)值模擬時，需進行多通道甚至全通道計算，這對計算資源提出了很高要求。為了簡化葉片排的非定常計算，引入葉間相位角(Inter-Blade Phase Angle，簡稱IBPA)概念對葉片排的計算模型進行簡化，這種簡化方法可以大大減小三維問題非定常計算的計算量。

通過實驗經(jīng)驗發(fā)現(xiàn)，當葉片振動時，葉輪機械上相鄰葉片之間存在固定的相位差，而相位差的大小與葉輪機的顫振特性緊密相關。于是，Lane[8]提出了具有葉間相位角的顫振模型，該模型假定振動中所有葉片振動頻率與幅值都相同，且相鄰葉片間存在固定相位差。當整個葉輪機振動時，波動的葉片間貌似存在一種橫波在葉片之間傳播，這種橫波被稱為行波模型。因為葉輪機械的典型特征就是旋轉對稱性，可以通過行波模型簡化葉片排計算模型，替代非定常計算中的多通道或是全通道計算模型。

假定葉排中的葉間相位角為σ，則葉排中每過2π/σ個葉片通道，葉片上下周期就會滿足周期性條件，因而整個葉片排計算模型就可以縮減到2nπ/σ個葉片通道，其中n為大于等于1的正整數(shù)。葉片排的簡化周期性邊界如圖1所示。

圖1 葉柵周期性邊界

通過移相周期性邊界條件將葉片排的計算域進一步簡化為單葉片通道。如圖2所示，假設行波方向與葉片運動方向相反，進行移相后，在周期性邊界條件上滿足條件：

(1)

從式(1)看出，任何相位差σ均可應用于移相周期邊界條件上。但在實際的葉輪中，葉片的數(shù)量是有限的，并非所有的IBPA都會出現(xiàn)。假定葉輪機械中有NBL個葉片，那么可能出現(xiàn)的葉間相位角為：

(2)

圖2 葉柵移相周期性邊界

1.2 能量法基本原理

能量法是屬于解耦方法的一種。假設葉片發(fā)生振動，并以某一固有模態(tài)振型的形式，在其對應的固有頻率下進行振動。然后通過計算葉片在一個振動周期內氣動力對葉片的做功情況來判斷能量的流動方向。若氣流對葉片做正功，葉片在振動過程中從氣流中源源不斷地吸取能量，那么振動將可能趨于發(fā)散；當氣流對葉片做負功時，氣動力在振動過程中起到氣動阻尼的作用，振動則趨于逐漸收斂。因此，可以通過定義氣動阻尼系數(shù)來判斷系統(tǒng)的穩(wěn)定性。由于能量法將氣動力的功作為判據(jù)，故將氣動阻尼定義為與其相對應的負值。

=Aζ0πsin(φ)

(3)

對所得到的模態(tài)氣動功進行傅里葉變換，得到一階諧波。并對氣動力系數(shù)無量綱化之后，定義氣動阻尼系數(shù)(Aerodynamic Damping)。

(4)

(5)

(6)

當阻尼系數(shù)Ξ<0時，系統(tǒng)趨于不穩(wěn)定；當阻尼系數(shù)Ξ>0時，系統(tǒng)趨于穩(wěn)定。

1.3 基于流固耦合方法的振動模態(tài)耦合求解

在對葉輪機械的葉片耦合問題進行研究時，需要對其結構動力學方程進行求解。認為葉片在振動過程中是屬于彈性形變的，因此采用振動模態(tài)對葉片結構的變形進行線性描述，假定其矢量位移為：

(7)

流固耦合問題的時域數(shù)值模擬過程中，流體域與固體域在交界面處的數(shù)據(jù)傳遞是非常重要的一步。由于在求解過程中，流體域網(wǎng)格與固體域網(wǎng)格在耦合交界面處完全不匹配，特別是流體域在葉片表面所生成的網(wǎng)格密度遠大于固體域所生成的網(wǎng)格，這樣就需選取一種數(shù)值插值方法來搭起流體域與固體域之間數(shù)值傳遞的橋梁。使用有限元方法計算葉片結構的模態(tài)振型之后，通過插值的方法將振型的變形后的坐標賦值于流體域的計算網(wǎng)格上。一般常用的數(shù)值插值方式有通用網(wǎng)格界面法(Generalized Grid Interface, GGI)與型函數(shù)法。

2 葉柵顫振特性分析與驗證

2.1 葉柵算例簡介

由于針對葉片結構特性測量的實驗數(shù)據(jù)很少，針對三維葉片結構特性的實驗數(shù)據(jù)更是無法獲取。因此，選擇二維葉柵Standard Test Configuration 4標準算例(簡稱STCF 4算例)作為驗證能量法預測葉片顫振特性可靠性的算例[9]。圖3是文獻中實驗所使用的試驗裝置實物圖，在渦輪靜子上共安裝了20個根梢比為0.8的渦輪葉片，葉片通過彈性梁與輪轂間進行連接。在顫振試驗中彈性梁發(fā)生彎曲，帶動葉片發(fā)生振動，而其振動方向與葉片弦向的夾角為60.4°。

圖3 STCF 4標準算例葉片實物圖

取該葉片的50%葉高位置進行二維簡化計算。在此葉片截面上，葉片的柵距為0.056 5 m，葉片弦長0.074 4 m，葉片安裝角為56.65°，進口氣流角與來流攻角在圖中均已有表示。

2.2 STCF 4葉柵氣動性能分析

對STCF 4葉柵的Case627實驗狀態(tài)進行定常狀態(tài)下的氣動性能研究。其中，進口氣流角為-15.2°，進口總壓為160 900 Pa，進口氣流總溫為288.15 K，進氣攻角為41.45°，出口背壓為100 700 Pa。

圖4表示STCF 4葉柵在Case627實驗狀態(tài)下，葉片表面壓力的數(shù)值計算結果與實驗測量結果的對比圖。從對比中發(fā)現(xiàn)，在此流動狀態(tài)下，數(shù)值計算與實驗測量的葉片表面壓力值的分布不僅在變化趨勢上相同，而且在數(shù)值上的分布也一致。

圖4 STCF 4 Case627葉片表面壓力分布對比

2.3 葉柵顫振特性分析

由STCF 4葉柵氣動性能結果分析可知，葉柵Case627工況氣動性能的實驗值與數(shù)值計算值吻合較好。因此，基于Case627工況進行簡化假設及非定常計算研究。圖5所示為文獻[10]中給出的STCF 4中的行波假設模型。

對葉片的運動位移進行簡化，其運動方程為：

h=h0cos(ωt+nσ)

(8)

式中：h為葉片沿振動方向(與葉片弦向夾角60.4°)的位移值，h0=0.000 299 m為運動中的位移幅值；ω=936 rad/s為振動的圓頻率；n為葉片編號；σ為葉間相位角IBPA。

圖5 STCF 4中行波方向的定義

非定常氣動力系數(shù)及氣動阻尼系數(shù)定義為：

(9)

(10)

圖6為Case627工況下，葉片的穩(wěn)定區(qū)域(IBPA)數(shù)值計算與實驗測量值的對比圖。從對比圖中分析發(fā)現(xiàn)，Case627中氣動阻尼系數(shù)在分布趨勢上吻合，對于葉柵彈性失穩(wěn)區(qū)域的預測與實驗數(shù)據(jù)基本一致。

圖6 Case627工況下氣動阻尼系數(shù)對比圖

通過以上標準二維葉柵算例的計算分析，驗證了基于流固耦合方法的能量累積法對于超/跨聲速工況的顫振預測的正確性。對于下文復雜幾何的三維轉子葉片的顫振特性問題研究奠定了堅實基礎。

3 某風扇轉子超聲速非失速顫振研究

圖7所示為在進行風扇葉片顫振特性數(shù)值分析時所采取方法的流程示意圖。

圖7 葉片顫振分析流程圖

其過程是先對固體結構進行結構動力學分析，分析預應力作用下的葉片模態(tài)振動狀況，得到葉片的振型、自振頻率等結構參數(shù)；再基于該轉子以此階模態(tài)振動的假設來進行非定常的流固耦合數(shù)值模擬；最后根據(jù)能量法計算葉片表面累積功或氣動阻尼系數(shù)來判斷其氣動彈性穩(wěn)定性。

選取該風扇轉子在設計轉速下的近效率最高點處的工況進行研究，此工況進口總壓Pt=101,325 Pa、總溫288.15 K，出口平均背壓Pb=115 kPa。由葉輪機械轉子葉片實驗研究發(fā)現(xiàn)，轉子葉片顫振受其前三階模態(tài)影響較大，高階模態(tài)由于其本身對應的葉片振動頻率很高，而實際中葉片發(fā)生振動的頻率較小且振動時幅值較大，所以文章主要對該風扇轉子前三階模態(tài)進行分析。

3.1 風扇轉子計算模型簡化

在進行風扇轉子的顫振研究時，假定葉片以某階模態(tài)按照前文假設的行波模型的方式進行簡諧振動，則其關于位移的運動方程如下所示：

h=φζ=φζ0cos(ωt)

(11)

式中：h表示葉片表面位移矢量；φ表示振型的矢量；ζ=ζ0cos(ωt)為系統(tǒng)中模態(tài)的廣義坐標；ζ0為廣義坐標幅值；ω為振動圓頻率。

定義風扇轉子非定常模型的行波方向與葉片轉動方向相反。為確定風扇轉子葉片的顫振邊界，先對轉子葉片的前三階模態(tài)振動的氣動阻尼進行求解，來判斷轉子在各階模態(tài)下是否出現(xiàn)負阻尼；然后針對出現(xiàn)負阻尼的不同模態(tài)狀況，計算對應模態(tài)頻率下葉片在各個IBPA的氣動阻尼曲線；最后針對出現(xiàn)負阻尼對應的葉間相位角IBPA，改變葉片的振動頻率，計算該模態(tài)及IBPA下的阻尼系數(shù)，通過插值方法得到氣動阻尼零點，該點即為風扇轉子對應模態(tài)下的超聲速非失速顫振臨界點。

3.2 風扇轉子顫振特性分析

基于能量法理論，風扇轉子對應的非定常氣動力系數(shù)定義為：

(12)

(13)

模態(tài)分析得到該轉子葉片前三階模態(tài)頻率分別為189.34 Hz、599.32 Hz、967.69 Hz。圖8所示為風扇轉子葉片一階振型表面非定常積累功率分布，氣動力做功主要集中于葉片的上半部分。分析發(fā)現(xiàn)，當葉間相位角IBPA=120°時，葉片上半部分表面有較大面積的氣流正功和負功區(qū)域，負功所占面積及強度較大，葉片發(fā)生顫振的可能性較??；當葉間相位角IBPA=334.3°時，葉片表面正功區(qū)面積及強度較大，其發(fā)生顫振可能性增加。將風扇葉片的一階模態(tài)對應的氣動阻尼系數(shù)隨葉間相位角變化的關系繪制如圖9，葉片在此階模態(tài)振動時，氣動阻尼系數(shù)在IBPA=334.3°附近出現(xiàn)負值，即出現(xiàn)顫振不穩(wěn)定區(qū)域，這與葉片表面的能量分布分析相一致。

圖8 轉子一階振型葉片表面非定常積累功率分布

圖9 f=189.34 Hz一階模態(tài)下葉片氣動阻尼系數(shù)

該風扇轉子葉片的第二階、三階模態(tài)對應的氣動阻尼系數(shù)隨IBPA變化的分布如圖10、圖11所示。分析發(fā)現(xiàn)，轉子葉片的第二階、三階模態(tài)對應的氣動阻尼系數(shù)均為正值，說明在此工況下，該轉子葉片不會發(fā)生基于第二階、三階振動模態(tài)的顫振現(xiàn)象。

圖10 f=599.32 Hz二階模態(tài)下葉片氣動阻尼系數(shù)

圖11 f=967.69 Hz三階模態(tài)下葉片氣動阻尼系數(shù)

通過數(shù)值計算及分析發(fā)現(xiàn)，對于該風扇轉子，一階模態(tài)為其最危險的模態(tài)，并且其最危險的IBPA在334.3°附近。控制葉間相位角IBPA=334.3°不變，改變風扇轉子葉片的振動頻率，得到其不同葉片振動頻率對應的氣動阻尼系數(shù)，通過進一步處理得到轉子的顫振臨界點。

圖12為一階模態(tài)氣動阻尼系數(shù)在臨界點附近隨葉片振動頻率的變化圖。依據(jù)線性插值的方法對氣動阻尼系數(shù)零點進行插值，當氣動阻尼系數(shù)為零時，其對應的葉片振動頻率為f=201.68 Hz，即為一階模態(tài)對應的超聲速非失速顫振邊界。對于該風扇轉子而言，第一階模態(tài)是其前三階模態(tài)中最不穩(wěn)定的模態(tài)，所以認定一階模態(tài)顫振邊界即為該風扇轉子的顫振邊界。從以上分析可知，當葉片的一階模態(tài)頻率低于f=201.68 Hz時，該風扇轉子葉片會以一階模態(tài)、IBPA=334.3°的形式發(fā)生超聲速非失速顫振。

圖12 IBPA=334.3°時一階模態(tài)氣動阻尼系數(shù)

4 結 論

采用基于流固耦合方法的能量累積法對某典型大尺度風扇轉子的超聲速非失速顫振特性進行了研究，得到以下主要結論：

(1) 風扇轉子葉片在近效率最大點可能發(fā)生基于一階模態(tài)、IBPA=334.3°的超聲速非失速顫振，其邊界的一階葉片結構固有頻率為f=201.68 Hz。

(2) 改變系統(tǒng)結構阻尼和材料阻尼等可提高葉片固有頻率可達到避免發(fā)生超聲速非失速顫振的目的。