非循環對稱力作用下離心葉輪的諧響應分析

肖 軍， 趙遠揚， 舒 悅

(1. 合肥通用機械研究院 壓縮機技術國家重點實驗室，合肥 230031； 2. 合肥通用環境控制技術有限責任公司，合肥 230031)

非循環對稱力作用下離心葉輪的諧響應分析

肖 軍1,2， 趙遠揚1,2， 舒 悅1,2

(1. 合肥通用機械研究院 壓縮機技術國家重點實驗室，合肥 230031； 2. 合肥通用環境控制技術有限責任公司，合肥 230031)

為研究流場中離心葉輪受氣流激振的受迫振動問題，利用葉盤結構的循環對稱性編制計算程序，實現了使用循環對稱模型進行非循環對稱力作用下離心葉輪的諧響應分析。針對空間周期性的來流條件，由各扇區所受周期激勵的相位差建立載荷向量，然后根據按扇區展開的葉輪各節徑模態的振型數據，建立模態坐標系下的諧響應運動方程，計算得到葉輪的諧響應位移;以一半開式徑向葉輪為例，進行了葉盤結構在氣動力非循環對稱分布情況下的諧響應分析。循環對稱模型與整體模型計算的位移響應結果符合較好，兩個共振點附近激振頻率下的振動響應均呈現出較高的幅值，表明了算法的正確性和預測氣流激勵導致共振的有效性。

氣流激振；非循環對稱力；諧響應；節徑；共振

現代工業離心壓縮機中，因周向進氣畸變、旋轉失速等原因誘發的氣流激振是導致葉片高周疲勞的重要因素，對由氣流激振引起的離心葉輪受迫響應問題展開深入研究對于抑制高周疲勞、提高葉輪的結構可靠性具有重要意義。

葉盤結構是一種典型的循環對稱結構，各葉片與其對應的輪盤部分組成葉盤結構的一個基本扇區。然而現代離心葉輪往往采用三元葉片，幾何形狀較為復雜，若對整個輪盤結構進行有限元建模，其網格規模通常很大，因而形成對計算資源的較高需求。為簡化計算分析工作，利用結構的循環對稱性，只對葉輪的一個扇區進行建模，進而求出整個葉輪的解，已成為目前工程上廣泛應用的理論分析方法[1-3]。目前多數商業有限元軟件雖然支持循環對稱結構的諧響應分析，但具體功能受限。如ANSYS軟件對壓力分布載荷只能指定零相位角，而離心葉輪各流道的氣動力往往并不同步變化，并且該軟件無法基于模態坐標展開循環對稱結構的諧響應分析，因此完全采用軟件實現氣流激勵下葉輪的諧響應計算受到諸多限制。

流場中的離心葉輪往往因上游尾流漩渦、周向進氣畸變及旋轉失速等原因形成氣流激振，目前關于氣流激勵引起葉片受迫響應的研究已有開展，但多數研究針對軸流壓氣機或汽輪機葉片。較近的如FRIDH等[4]研究了畸變插板導致部分進氣條件下軸流渦輪葉片的受迫振動問題，研究表明，由部分進氣引起的流動阻塞和抽吸直接導致了大量低階的振動響應，并且轉子數目與低階激勵的組合是誘發共振的重要原因。BLADH等[5]以一軸流壓氣機后部級為對象，研究了輪轂段軸向間隙的泄漏流動對下游轉子葉片的激振問題。結果表明，因間隙泄漏流動誘發的激振力遠高于上游尾跡引起的激振力。針對離心壓縮機葉輪的氣流激振問題亦有相關研究，DICKMANN等[6]針對一渦輪增壓器葉輪的受迫振動展開了實驗和數值研究，氣流激勵由葉輪與蝸殼和進氣彎管的干涉誘發，將流場計算得到的葉表非定常氣動力逐點映射到有限元網格上，求得的葉輪應變數據與實驗吻合很好。WHITE等[7]分析了某5級離心壓縮機葉輪的高周疲勞破壞問題，研究綜合歷史運行數據及數值分析結果，認為上游靜子葉片的尾跡流是葉輪疲勞斷裂的主要激振源，而處于堵塞區間的工況亦加劇了激振。然而，已有研究多使用商業軟件實現計算，且對葉盤結構往往采用整體建模，顯著增大了計算規模。亦有研究利用循環對稱性對一個扇區建模計算受迫響應，毛義軍等[8]采用諧響應方法計算了某葉輪在氣動力作用下的振動響應，其預測的葉片應力集中位置和實際斷裂位置一致，但該研究只考慮了頻率因素，忽略了激振階數及輪盤節徑對共振的影響機制?；谝陨显?，結合前述商業軟件的功能限制，有必要對離心葉輪的受迫響應問題展開深入的分析研究。

本文在非定常流計算得到葉表氣動力的基礎上編制程序，實現了利用循環對稱性對離心葉輪進行氣流激勵下諧響應分析的完整計算流程。對一半開式徑向葉輪展開周向畸變來流條件的諧響應分析，計算結果顯示了算法預測葉輪受氣流激勵導致共振的有效性。

1 數值方法

1.1 循環對稱結構的諧響應動力方程

對于任意葉盤結構，激振力作用下的無阻尼結構動力方程可表示為

(1)

式中:矩陣和向量均在各扇區局部坐標系下組裝，質量矩陣M對稱正定，剛度矩陣K對稱半正定，已包含離心力引起的預應力效應，F為外載荷向量，Δ為位移向量。定義坐標變換

(2)

V為范德蒙矩陣，定義為

(3)

式中:N表示扇區數，

(4)

(5)

將式(2)代入式(1)并左乘VH得到

(6)

式中:r表示節徑，且有

(7)

(8)

這樣結構動力方程就被正交分解為N個子空間的方程，復雜度大大降低。另根據式(9)，只需求解r= 0, 1, …,Nf對應的方程，N為偶數時，Nf=N/2，N為奇數時，Nf= (N-1)/2。針對循環對稱結構，建立n個自由度的單扇區有限元模型，求解各子空間的方程式(6)，再由式(2)即可得到位移解。

對于多載荷步的諧響應問題，需要在每個載荷步上求解n自由度方程式(6)，計算仍較為耗時，且得到的結果是所有結構點位移的頻域響應，數據量較大。工程應用中往往關注激振力頻率及階次是否導致結構某階振型的共振響應，為此考慮將結構動力方程轉換至模態空間求解。令qr=ψreiΩrt，在子空間中求解廣義特征值問題

(9)

即可得到n個廣義特征值Ωr及n組廣義特征向量ψr，且考慮到式(8)，有

(10)

將ψr對αr正交歸一化，有

(11)

以N為偶數為例，由式(2)得到按扇區展開的正交歸一實模態矩陣如下

(12)

該模態振型是在各扇區的局部坐標系中表達的，而ANSYS軟件計算的各扇區模態振型是在總體坐標系中表達的。為此本文先使用ANSYS建立基本扇區模型，導出其各節徑下的剛度矩陣及質量矩陣，之后編制程序調用大型稀疏矩陣解算器FEAST[9]求解廣義特征值問題式(9)，得到基本扇區各節徑的振型Ψr。FEAST每次求解得到的Ψr可相差一個常數eiδ，但不影響式(12)的正交歸一性。

引入變換

(13)

(14)

其中

(15)

f(t)=ΨHF(t),F(t)=[F1(t)F2(t) …FN(t)]T

(16)

計入阻尼后，式(14)成為

(17)

為進行諧響應計算，移除平均氣動力，將非定常氣動力表達為若干個諧波的疊加，即

(18)

(19)

(20)

1.2 周向畸變來流下的激振力分布

實際運行中離心葉輪可能經歷如上游尾渦脫落、旋轉失速等復雜工況，使葉輪流場呈現出遠離設計狀態的非定常變化特征，此時的進氣參數往往沿周向不均勻分布，并呈現空間周期性的變化(見圖1)，θ為周向角，U為氣流參數。以沿周向運動的失速團為例，動葉每掃過一個失速團，在葉輪坐標系中觀察，該葉片通道的流場即變化一個周期，相鄰動葉通道的流場變化只相差一個固定時間，體現在激振力上，即為相鄰扇區的一個固定相差?；儦饬餮刂芟蜻\動，設T為動葉掃過一個畸變區的時間，M為畸變區數目，則

ΔΦj=2πjM/N

(21)

圖1 周期畸變來流條件

(22)

1.3 氣動力插值及等效結點載荷

流場計算與諧響應分析均用各自的網格體系，在流體網格上獲取的氣動力數據需向結構網格插值，并轉換為等效結點載荷，才能進行諧響應方程的求解。首先針對每個結構表面點，根據流體表面網格中心的坐標搜索距其最近的流體表面網格，然后用反距離加權法插值得到該結點的氣動壓力值。設結構點到距其最近的流體表面網格頂點的距離為di，該表面網格頂點的壓力為pi，則結構點壓力為

(23)

氣動壓力作用于單元面，需將其轉換為等效結點載荷后再施加于結構方程右端。根據虛功原理，等效結點載荷應與單元面分布載荷在任何虛位移上所做的功相等。對于結構表面點i有

(24)

式中:fi為等效結點載荷;Ni為形函數，n為單元面法向量。

2 算例分析

本文針對離心葉輪在周向畸變來流下的受迫振動展開諧響應分析，研究對象為一半開式徑向葉輪，出口直徑280 mm，轉速14 250 r/min，葉片數16，葉輪上游畸變區數目為13。對流場的非定常計算完成后，獲得以傅里葉系數形式存儲在流體網格上的氣動壓力?？烧J為氣動壓力的各階諧波項僅與來流參數的空間畸變有關，而與畸變區相對轉子的周向運動速度無關。提取前幾階諧波項，用“1.3”節所述方法將其插值到結構表面網格。圖2所示為分別由氣動網格和結構網格表達的一階諧波的葉表非定常壓力幅值分布，顯見相同色溫設置下，幾乎所有區域的數據都符合很好，在葉片前緣附近，因氣動網格較密，兩圖略有差異。

圖2 葉表非定常壓力幅值的分布(一階諧波)

將包含幅值和相位信息的壓力諧波項轉換為等效結點力的形式后，即可進行諧響應方程的求解。計算采用10節點四面體單元，葉輪材料的彈性模量為7.2×1010Pa，泊松比0.33，密度2 700 kg/m3，對軸孔表面施加徑向及周向位移約束，對靠軸肩的結點施加軸向約束。材料阻尼比取0.002，氣動阻尼比按照文獻[10]的方法估計。

(a) 循環對稱模型 (b) 整體模型

分別采用循環對稱模型和整體模型進行葉輪的諧響應計算，兩套結構網格如圖3所示，網格數據在表1中列出。由表可以看出，整體模型的數據規模約為循環對稱模型的16倍。對于本文基于模態坐標的諧響應分析，整體模型最主要的問題在于各扇區的網格劃分不能完全做到循環對稱的分布，因此無法準確呈現各節徑模態的振型，進而影響到諧響應計算的結果無法準確捕捉激振力對某節徑模態的共振激勵效果。表2所列為采用循環對稱模型和整體模型計算得到的振型頻率數據，顯見對于循環對稱模型同一節徑和頻率的兩個振型，整體模型計算得到的頻率略有差異，可以判斷整體模型計算的這兩個振型亦不嚴格正交。振型計算的差異將對諧響應分析的結果產生影響。

表1 有限元模型規模比較

表2 振型頻率比較

設計算的頻率范圍為1 500～3 500 Hz，分析了非定常壓力一階諧波作用下的諧響應結果，考慮到基于模態坐標的諧響應計算速度較快，每隔1 Hz即設置一個載荷步。圖4比較了采用循環對稱模型和整體模型計算的基本扇區葉片前緣葉頂結點的諧位移幅值，顯見兩種模型計算的x,y,z向位移響應相差均很小，顯示了基于循環對稱模型算法的正確性，結果的差異是由不同結構網格計算的振型數據存在差別所致。阻尼的存在主要對諧響應的相位產生影響，并限制共振頻率下響應幅值的增長。從圖4可知，激振力作用下位移響應在圖示頻率范圍內出現兩個峰值，相應頻率分別為2 016 Hz和3 184 Hz，這正與第3節徑的前兩個振型頻率重合。由表2可知,節徑3～8節徑的前一兩個模態頻率非常接近，但位移響應峰值不會出現在圖示頻率范圍非節徑3的頻率上。實際上，對于本文畸變區數目為13的情況，只有屬于節徑3的振型才較為危險，其它振型并不會被激發。對于循環對稱的葉盤結構，只有滿足[4]的條件才可能導致共振。

jM±r=lN

(25)

式中:j為激振力階次，r為節徑，l為任意整數。

綜合激振頻率和上述節徑條件，離心葉輪與畸變來流在相對轉速9 306 r/min(9 306/60×13 ≈ 2016)下的共振預測如圖5所示，實心黑表示不同節徑下葉輪的固有頻率，黑實線條線表示滿足式(25)的激振頻率，黑圓圈表示可能引起共振的區域。對于本文算例，畸變區和葉輪數目分別為13和16，一階諧波的激振力只能激發節徑3的振型，該節徑前兩個正交模態2 016.4 Hz的頻率剛好與9 306 r/min相對轉速下的激振頻率高度吻合；二階諧波的激振力可以激發節徑6的振型，但非定常壓力的二階諧波項較小，未計算其諧響應。

圖4 基本扇區葉片前緣葉頂結點的諧位移幅值比較

圖5 相對轉速9 306 r/min下的共振預測

圖6給出了不同結構模型計算的節徑4模態1的廣義位移響應，對整體模型還給出了節徑3模態1的結果作為比較。由圖6可知，循環對稱模型計算的諧響應位移幅值在共振區內的數量級僅為10-19，從數值角度看相當于零位移，而整體模型計算的位移幅值在共振區內的數量級為10-6，雖不是零位移，但比節徑3模態1的位移幅值已大大減小。當對葉輪整體建模的結構網格能完全呈現葉盤結構的循環對稱性時，位移幅值的數量級會進一步減小。其它節徑模態的廣義位移響應亦呈現類似的情況，如圖7所示節徑0模態3的諧位移幅值同樣極小，由激振條件(25)知該模態不會被激發。以上結果顯示了基于循環對稱模型的算法預測氣流激勵導致葉輪共振的有效性。

圖6 不同結構模型計算的廣義位移響應(節徑4模態1)

Fig.6 Generalized displacement responses calculated by two structural models (Mode 1 of nodal diameter 4)

圖7 節徑0模態3的諧位移幅值(循環對稱模型)

Fig.7 Harmonic displacement amplitude of mode 3 of nodal diameter 0 (Cyclic model)

3 結 論

利用離心葉輪的循環對稱結構，提出一種計算非循環對稱氣流激勵下葉輪諧響應的算法。實現了網格界面氣動力插值、等效節點載荷轉換、振型解算及擴展等關鍵技術。針對一氣流激勵下的半開式徑向葉輪，計算得到各節徑模態的諧響應，結果顯示算法能有效預測該葉盤結構的共振問題。算法利用循環對稱性建模，使數據量較葉輪整體建模大大減小，計算表明兩種建模方式得到位移響應結果符合較好。

本文研究僅針對諧調葉盤結構，而有關失諧葉輪的諧響應計算，需要進一步展開深入的研究。

[ 1 ] DICKENS J M. Numerical methods for dynamic substructure analysis [D]. Berkeley: University of California, 1980.

[ 2 ] PETROV E P. A method for use of cyclic symmetry properties in analysis of nonlinear multi-harmonic vibrations of bladed disks[J]. Journal of Turbomachinery,2004,126(1): 175-183.

[ 3 ] 李松濤, 許慶余, 張小龍. 透平機械離心壓縮機和離心泵葉輪動力分析的三維旋轉循環對稱CN群算法 [J]. 應用力學學報, 2005, 22 (2): 169-174. LI Songtao, XU Qingyu, ZHANG Xiaolong. Application of CNgroup cyclic symmetry method to three-dimensional dynamic analysis of centrifugal compressor and pump impeller [J]. Chinese Journal of Applied Mechanics, 2005, 22 (2): 169-174.

[ 4 ] FRIDH J, LAUMERT B, FRANSSON T. Forced response in axial turbines under the influence of partial admission [R]. ASME Paper GT2012-68303, 2012.

[ 5 ] BLADH R, ZHUANG Qingyuan, HU Jiasen, et al. Leakage-induced compressor blade excitation due to inter-segment gaps [R]. ASME Paper GT2012-70040, 2012.

[ 6 ] DICKMANN H P, WIMMEL T S, SZWEDOWICZ J, et al. Unsteady flow in a turbocharger centrifugal compressor: 3D-CFD-simulation and numerical and experimental analysis of impeller blade vibration [R]. ASME Paper GT2005-68235, 2005.

[ 7 ] WHITE N, LANEY S, ZORZI C. RCFA for recurring impeller failures in a 4.7 Mtpa LNG train propane compressor [C]// Proceedings of the Fortieth Turbomachinery Symposium, Houston, USA, 2011.

[ 8 ] 毛義軍, 祁大同, 許慶余. 離心壓縮機葉輪葉片疲勞斷裂故障的數值分析 [J]. 西安交通大學學報,2008,42(11): 1336-1339. MAO Yijun, QI Datong, XU Qingyu. Numerical study on high cycle fatigue failure of centrifugal compressor impeller blades [J]. Journal of Xi’an Jiaotong University, 2008, 42(11): 1336-1339.

[ 9 ] POLIZZI E. A high-performance numerical library for solving eigenvalue problems [R]. FEAST Solver v2.0 User’s Guide, 2012.

[10] ZEMP A, ABHARI R S, RIBI B. Experimental investigation of forced response impeller blade vibration in a centrifugal compressor with variable inlet guide vanes-Part 1: blade damping [C]// Proceedings of ASME Turbo Expo 2011, Vancouber, Canada, 2011.

Harmonic response analysis for a centrifugal impeller under non-cyclic symmetric exciting forces

XIAO Jun1,2, ZHAO Yuanyang1,2, SHU Yue1,2

(1. State Key Laboratory of Compressor Technology, Hefei General Machinery Research Institute, Hefei 230031, China；2. Hefei General Environmental Control Technology Co., Ltd., Hefei 230031, China)

In order to investigate forced vibration problems of centrifugal impellers subjected to gas excitation in flow fields, an algorithm to calculate harmonic responses of a centrifugal impeller under non-cyclic symmetric exciting forces was realized with a cyclie and symmetric model considering the cyclic symmetry of impeller structures. For the spatially periodic inflow conditions, firstly the load vector was constructed according to the phase differences of periodic excitations of sectors, and then utilizing the modal shape data of different nodal diameters of the impeller its harmonic response dynamic equations were established under the modal coordinate system to calculate its displacement responses. Finally, the harmonic response analysis for a radial impeller under non-cyclic symmetric aerodynamic exciting forces was performed. It was shown that the displacement response results of the radial impeller obtained with the cyclic and symmetric model agree well with those with the whole body model; the responses under the exciting frequencies near two resonance points have relatively higher amplitudes; the correctness and effectiveness of the proposed algorithm to predict the resonance of a centrifugal impeller caused by gas excitation are verified.

gas excitation; non-cyclic symmetric force; harmonic response; nodal diameter; resonance

國家青年科學基金 (2014011155);國家973計劃項目(2012CB026000)

2015-06-16 修改稿收到日期：2015-11-11

肖軍 男，博士，副研究員，1979年生

TH113.1