基于非線性諧波法評估壓氣機轉靜干涉誘發的強迫響應

孫海，王培屹，鄭寧，楊琳，陳育志

（1.中航工業沈陽發動機設計研究所，沈陽110015；2.北京航空航天大學能源與動力工程學院，北京100083）

0 引言

在葉輪機械工程應用中，流體誘發振動導致葉片斷裂是1種主要的結構破壞形式。其引發的故障在各類發動機故障中占相當高的比例，且危害較大，是航空燃氣渦輪發動機結構完整性和可靠性的主要威脅。而轉靜干涉是葉輪機械流體誘發強迫振動的重要因素，是由轉、靜子葉片間的相互轉動而產生的非定常流動現象。國內外雖然對由轉靜干涉引起的尾流激振或者由勢流擾動而引起的強迫響應開展了一些研究，但在氣動層面非定常脈動壓力的計算精度、效率和強迫響應阻尼比的取法以及工程適用性等方面還存在不足。

本文采用非線性諧波法計算非定常流動產生的諧波壓力，進而計算該壓力在葉片上的強迫響應，計算精度和效率都有一定提升，建立了適合工程應用的、能夠評估由于轉靜干涉誘發的強迫響應的方法。

1 研究背景

一般來說，葉輪機械中的流體誘發振動問題可以分為氣動彈性穩定性問題和響應問題2種基本類型。氣動彈性穩定性問題對應葉輪機械的自激振動，而響應問題對應強迫振動。與葉片的自激振動相比，在航空發動機的發展歷程中，由流體激勵誘發的強迫振動所造成的葉片振動破壞故障的次數要遠超過由葉片顫振造成的故障。

當前對轉靜干涉研究較多的是上游葉片對下游葉片的尾流激振和下游葉片的勢流擾動對上游葉片非定常表面力的周期性影響。

由于尾流激振力或勢流擾動很難用具體的解析表達式表示，因而數值模擬計算尾流激振力是工程中最主要的方法。A.J.Sander和Fleeter[1]研究了1種可用于預測軸向和離心式葉輪機械響應的模型，用于分別計算葉片排在上下游勢流和在上游尾跡作用下的強迫響應。結果證明，無論是在軸流還是在離心式壓氣機中，葉片的強迫響應幅值為同一量級。因此，對于葉片排間距較近的壓氣機，流勢與葉片尾跡具有同樣的重要性。Y·T·Lee和J·Z·Feng[2]等研究了不可壓流在多級葉排的流動情況，分別研究了轉子/靜子、靜子/轉子2種模型。研究指出，在軸向間距為10%弦長時，尾跡效應提供了氣動載荷峰值的10%。針對轉子數目的影響，研究了轉子1/靜子/轉子2模型，并且轉子1和靜子、轉子2和靜子之間的軸向間距可調，發現不同軸向間距最小值到最大值的非定常氣動力變化較大。Bjorn Laumert和Hans.Martensson[3]通過3維數值模擬，得出在亞聲速工況下，轉子葉片壓力面上壓力值主要受靜葉尾流的影響，而在吸力面，壓力值主要受靜葉勢流和尾跡共同影響。由于3維轉子葉片沿徑向形狀不同，使得葉片表面壓力分布沿徑向也有所不同。葉片表面激振力的大小沿葉根向葉尖幅值逐漸減小。Stuart M[4]提出1種松散的耦合方式對葉片在氣流激振力下強迫響應進行預估計算。假定氣動力對葉片的模態沒有影響，先在ANSYS中計算葉片某一模態的振型和頻率，然后將其帶入CFD進行頻域求解，計算氣動力和阻尼力，再帶回有限元模型計算該模態下的振動響應。王梅等[5]通過研究均勻葉柵的尾跡引起的高頻振動問題，把已有的氣動計算方法和強度計算方法整個串起來，為工程應用初步建立起了1個尾流激振情況下葉片振動應力預估的半經驗方法。具體是通過計算轉子通道流場，得到轉子葉片表面的非定常壓力幅值沿葉片表面的分布，然后在結構計算中將葉片表面的壓力幅值等效地轉化為有限元節點上的節點力，再通過有限元方法計算葉片的振動應力。但是，非定常流場計算采用的參數多項式方法精度不足，阻尼的取法也使響應計算的準確性受到影響。孟越等[6]針對轉子葉片在前排靜子葉片尾流激振[7]情況下的位移和應力進行預估計算，采用3維非定常求解，將3維非定常氣動力引入轉子葉片有限元結構計算中，對尾流激振情況下葉片強迫響應問題[8]進行瞬態分析。通過分析可知，葉片響應大小與葉片表面所受氣動激振力大小、相位和分布有關。寧方飛[9]以1.5級高負荷風扇為研究對象，利用諧波法[10-11]捕獲了轉子激波向上游的傳播，以及轉子尾跡對下游靜子的影響，結果與常規非定常方法吻合，但并未涉及強迫響應計算。

2 非定常脈動力的研究對象、計算方法和結果

2.1 研究對象

針對某風扇試驗件進口可變彎度導向葉片開展數值模擬和試驗研究。風扇進口可變彎度導向葉片分為固定部分和可調部分，可調部分可繞轉軸旋轉調節以滿足不同轉速轉子預旋的需求，由于導葉可調部分只在轉軸處固支，與常規靜子的2段全部固定相比，更容易受轉靜干涉影響而被誘發振動。因此，針對后排轉子葉片的位勢作用對上游可調葉片的強迫響應進行分析。

2.2 非線性諧波法

非線性諧波法（nonlinear harmonic method）能夠在一定的計算量范圍之內，提供1種介于定常和非定常之間的計算方法，其計算精度較高能夠滿足計算要求。該方法將流體非定常擾動傅里葉分解為不同階次諧波疊加?；谇蠼鈺r間平均和不同諧波階次的定常輸運方程，能夠有效地在頻域內求解非定常問題，尤其是對于轉靜干涉問題有較高的計算效率和精度。

轉靜干涉是葉輪機械流動中最常見的非定常因素，對轉、靜子通道內流動的影響具有周期性。采用諧波函數逼近，經過周期性擾動近似以及線化假設，可以將N-S方程分解為時均方程和擾動方程，擾動方程的個數由諧波的階次確定，各階擾動方程之間相互獨立，在諧波的頻率給定之后，可以將擾動方程由時域轉化到頻域。最后在頻域內對時均方程和各階擾動方程耦合求解，就能得到流場各物理量的時均值和各階擾動值，將其疊加起來就得到近似解。可以通過選擇諧波的階數來控制求解的精度，諧波階次越高，求解精度也就越高。

引入周期性擾動假設后，在非定常流動中，守恒型變量可以分解為時均值與周期性擾動之和，而每個周期性的擾動可以用N階諧波來逼近，即

采用有限體積法的守恒型非定常流動方程可表示為

式中：Ωi為網格單元的體積；Fc和Fv分別為離散的對流項和黏性項；S為面積；Q為源項。

具體展開式見文獻[12]。展開式中包含的未知數數目是原N-S方程的2N+1倍，求解諧波方程僅相當于求解2N+1個時刻的定常方程。理論上，計算的諧波階次趨于無窮，其計算結果應當趨近于時間推進的非定常求解結果。一般選取前若干階（小于5）即可滿足工程需求。而且，一般認為，第1階非定常擾動在幅值上占絕對主導地位（對于尾跡或勢流造成的非定?，F象），高階擾動因其頻次過高，其距風扇葉片低階共振頻率過遠，因此影響較小。

2.3 計算設置及結果分析

網格使用AUTOGRID5生成。Vilmin認為若要成功運用N L H法模擬葉排間的流動狀況，上下游葉片周向網格數需要達到一定要求，具體見文獻[13]。根據網格數目限定，做出風扇網格，采用NUMECA的Fine進行計算。

以該多級風扇試驗件作為流場計算對象。計算時給定進口總溫28815K，進口總壓101325Pa，給定出口靠近輪轂位置靜壓值，且靜子出口靜壓滿足簡單徑向平衡方程。采用Harmonic Method進行非定常計算。非定常流場求解以相同進、出口邊界條件下的定常結果為初場。為提高計算效率，每個擾動的諧波數設置為2，每個葉排最大的擾動數設置為2。

風扇進口可調導葉第1階諧波壓力實部、虛部及幅值如圖1所示。第1階諧波壓力反映了后排轉子的非定常勢流影響，主要分布在可調導葉表面上半部，由此判斷，諧波壓力幅值較大區域可能會引起較大的強迫響應。

風扇后排轉子50%、90%葉高相對馬赫數如圖2所示。后排轉子馬赫數圖中的空白部分為超聲區?？烧{導葉表面依舊是諧波壓力幅值。在轉子90%葉高吸力面有較強激波，而在50%葉高處沒有激波。由于可調導葉和轉子存在相對轉動，轉子上半部分存在的激波通過勢流作用不斷掃掠前排可調導葉，影響可調導葉上半部，尤其影響靠近尾緣的后半部分邊界層狀況，造成較大的諧波壓力脈動。

圖1 風扇進口可調導葉第1階諧波壓力

3 強迫響應分析

本文建立的強迫諧響應的分析模型，應用MATLAB和ANSYS的2次開發語言APDL編制流場與結構的接口程序，將流場計算的3維非定常激振力完整加載到結構中，并對結構進行強迫響應瞬態分析，具體計算流程如圖3所示。

選取可調導葉部分（略去固定支板）進行有限元建模，葉片在厚度方向分為2層，對前緣和尾緣進行了適當加密，整體網格數約為1.6萬，如圖4所示。模擬可調導葉真實的約束情況，給定葉尖前緣部分全約束，葉根前緣部分周向和軸向約束，可繞徑向軸轉動。

圖2 風扇后排轉子50%、90%葉高相對馬赫數和第1階諧波壓力幅值

圖3 強迫諧響應分析系統

圖4 有限元模型及施加約束

利用SURF154單元進行脈動壓力的傳遞，利用F E M模型表面節點以及CFD模型靠近葉片的第1層網格節點插值，如圖5所示。圖中紅色節點為F E M模型可調導葉部分，藍色為CFD模型（包括固定支板）。程序自動尋找較近點進行插值，因此多余的固定支板不影響插值結果。

對葉片進行模態分析，找尋到與勢流激勵相近的模態階數。在氣動載荷的激勵下，進行諧響應分析，模態阻尼比取5e-4von Misesstress如圖6所示。除去轉軸處固支帶來的應力集中，葉身振動應力幅值約為20MPa，與試驗應力水平相當。

不同模態阻尼比振動響應分析的動應力結果見表1。從表中可見，模態阻尼比是振動響應分析的重要影響因素。如果試驗能夠獲得較為準確的阻尼數據，則采用該方法能夠較好地評估葉片由于尾流激振或勢流擾動造成的強迫響應。然而，在工程設計階段，沒有試驗件的生產，無法通過試驗獲得阻尼數據，因此很難準確確定模態阻尼比，導致諧響應分析存在較大誤差。這就要分析在不需要阻尼比的情況下，評估這種強迫響應造成的高周疲勞水平。

圖5 CFD與FEM模型插值節點

圖6 阻尼比為5e-4時諧響應von Missesstress

表1 不同模態阻尼比下響應分析的動應力結果

4 模態激勵因子評估強迫響應

如前所述，在葉輪機械中，由轉靜干涉引起的葉片強迫響應振動主要是以轉速為基頻的周期性激勵。尤其是在某些特定階激勵頻率（可引起危險共振的，如第n階激勵頻率）作用下葉片的強迫響應。

與一般激勵不同，這里的激勵為非同步激勵，即除了頻率以外，在任意時刻，流場中的任意1點（包括葉片表面上）的作用力由2個參數確定，即：幅值和相位。這2個參數均是點的空間位置坐標的函數。

尾流激勵或者是勢流擾動所具有的非同步激勵場的特性，決定了葉片的強迫響應既取決于激振力場又取決于模態位移場。模態激勵因子可以包含這2方面性質，能夠對某階共振的相對危險性做出一定評價。

激振力、模態位移都以場的形式出現，其中激振力場又是場的復數形式的疊加。這就使得激勵-響應之間的關系可以從場的角度描述?？梢砸霘W氏空間的相關概念，從數學上討論激勵場與模態位移場的空間相對位置對葉片響應的影響，從理論上給出模態激振因子的解釋。

定義模態激勵因子[14]

式中：｛F｝為激勵的幅值向量；[Λ（ejθ）]為以激勵相位為對角線元素的對角矩陣；｝為對系統質量矩陣歸一的第s階模態向量。模態激勵因子的幾何意義為尾流激勵的空間分布與葉片振動模態的夾角。

激勵因子的幅值反映了尾流激勵或勢流擾動在某階模態投影的大?。幌嘟欠从沉思ふ窳雠c模態位移場間的夾角，夾角越小，作用越強。

以某風扇導向器葉片與轉子軸向間距改變對尾跡激振的影響為對象，對非同步激勵場作用下葉片高激振力場階共振的危險性評定方法進行校核，結果如圖7所示。

圖7 激勵因子及其和相角隨軸向間距變化

從圖中可見，隨著軸向間距的增加，激勵因子的模呈下降趨勢，這與實際物理現象吻合；相位變化，反映了上游尾跡進入下游轉子通道的時間先后發生了變化。可見，激勵因子能夠直觀定量地反映出結構、氣動參數變化對結構共振響應造成的影響，能夠真實地反映物理實際。

對可調導葉，在計算過程中提取載荷向量和模態向量，計算出激勵因子，如圖8所示。按經驗，模態激勵因子與之前計算的振動應力水平相當。本文雖然對轉靜干涉中的后排對前排的勢流擾動進行計算，其方法也同樣適用于尾流激振。

圖8 可調導葉勢流激勵狀態下的模態激勵因子

5 結論

利用非線性諧波法以及插值程序，把已有的氣動計算方法和強度計算方法結合起來，為轉靜干涉中的尾流激振和勢流擾動產生的強迫響應建立了適用于工程設計的、合理的流固耦合計算模型及研究流程。

（1）非線性諧波法可以有效計算轉靜干涉導致的勢流擾動和尾流激振脈動壓力，精度和計算效率能達到工程要求，比定常和非定常計算有有一定的優勢。

（2）對于強迫響應分析，結合強迫響應計算的應力和模態激勵因子，可以對這種勢流擾動或尾流激振誘發的強迫響應作出相應評估。

（3）隨著經驗的積累，模態激勵因子可以單獨作為評估參數，并可以在完成氣動設計后對這種轉靜干涉進行評估，迅速地回饋氣動設計，減少設計流程和時間，而且可以避免引入阻尼比的誤差，提高工程上評估這種高周疲勞的精度和效率。

[1] Sander A J,Fleeter S.Blading response to potential field interactions in axial and radial flow turbomachinery[J].Journal of Propulsion and Power，1998，14（2）：199-207.

[2] Lee Y T,Feng J Z.Potential and viscous interactions for a multiblade row compressor[C]//ASME Turbo Expo 2003，Collocated with the 2003 International Joint Power Generation Conference,USA:American Society of Mechanical Engineers，2003.

[3] Bjorn Laumert，HansMartensson.Investigation of unsteady aerodynamic blade excitation mechanisms in a transonic turbine stage,Part I:phenomenological identification and classification[J].Journal of Turbomachinery，2002，124（3）：410-418.

[4] Wei Ning，StuartM.Blade forced response prediction for industrial gas turbines—Part2:verification and application[C]//ASME Turbo Expo 2003,Collocated with the 2003 International Joint Power Generation Conference，USA:American Society ofMechanical Engineers，2003.

[5] 王梅，江和甫，呂文林.在尾流激振情況下葉片振動應力預估技術[J].航空動力學報，2007，22（4）：608-613.WANG Mei，JIANG Hefu，LV Wenglin.Method to predict the blade vibration stress induced by wake flow[J].Journal of Aerospace Power，2007，22（4）：608-613.(in Chinese)

[6] 孟越，李琳，李其漢.尾流激振情況下葉片強迫響應瞬態分析方法[J].北京航空航天大學學報，2006，32（6）：671-674.MENG Yue，LI Lin，LI Qihan.Transient analytical method of vane forcing response under stator-rotorwake influence[J].Journal of Beijing University of Aeronautics and Astronautics，2006，32（6）：671-674.(in Chinese)

[7] 孟越，李琳，李其漢.不對稱靜子尾跡流場激振力分析及計算方法[J].北京航空航天大學學報，2007，33（9）：1005-1008.MENG Yue，LI Lin，LIQihan.Investigation of force under asymmetry statorwake[J].Journal of Beijing University of Aeronautics and Astronautics，2007，33（9）：1005-1008.(in Chinese)

[8] 李潤澤，李琳.跨聲速工況下流體誘發葉片振動研究[J].航空動力學報，2008，23（4）：747-753.LIRunze，LI Lin.An investigation of flow induced blade vibrations at the transonic operating condition[J].Journal of Aerospace Power，2008，23（4）：747-753.(in Chinese)

[9] Du Pengcheng，Ning Fangfei.Application of the harmonic balance method in simulating almost periodic turbomachinery flows[R].ASME 2014-GT-25457.

[10] Ning F,Xu L.Numerical investigation of transonic compressor rotor flow using an implicit 3D flow solver with one-equation spalart-all maras turbulence model[R].ASME 2001-GT-0359.

[11] Ning F.MAP：a CFD package for turbomachinery flow simulation and aerodynamic design optimization[R].ASME 2014-GT-26515.

[12] He L.Harmonic solution of unsteady flow around blades with separation[J].AIAA Journal，2008，46（6）：1299-1307.

[13] Vilmin S，Hirsh C，Lorrain E.Unsteady flow modeling across therotor/stator interface using the nonlinear harmonic method[R].ASME-2006-GT-90210.

[14] Li Lin，Wang Peiyi.Evaluation of high-order resonance of blade underwake excitation[R].ASME 2010-GT-23148.