劉正先，王生玲

(天津大學機械工程學院，天津 300350)

離心葉輪旋轉失速和喘振造成了流場的不穩(wěn)定性，限制了設備穩(wěn)定工況范圍，制約了設備安全高效運行．目前，主要通過觀察性能曲線中峰值流量和靜壓脈動特征，運用時間傅里葉分析、行波能量法、相關系數(shù)法、小波分析和空間傅里葉分析對捕捉到的信號進行處理，從而判定葉輪是否進入失速狀態(tài)及失速團運動規(guī)律．

Camp等[1]和 Day[2]通過觀察性能曲線峰值點與失速點流量的關系，及機匣上靜壓是否發(fā)生突變判定失速，其準確度受經(jīng)驗和主觀因素制約．吳艷輝等[3]通過對非定常壓力頻譜分析，根據(jù)失速擾動頻率分析失速原因，但傅里葉變換結果是整個時間范圍內(nèi)的頻率特性，而不能反映信號在某一時刻的本質特征．程曉斌等[4-5]、Lin 等[6]、符嬈等[7]和 Zhang 等[8]等通過小波分析識別軸流葉輪內(nèi)失速信號，可以在時頻域內(nèi)分析失速流場，但受時域和頻域分辨率的制約和折中，對離心葉輪失速流場中復雜的壓力信號不能得到理想的結果．

空間傅里葉分析計算簡單，實時性好，不受采樣頻率影響，可以得到擾動信號空間不均勻程度隨時間的變化規(guī)律．Garnier等[9]通過分析軸流葉輪失速實驗數(shù)據(jù)的空間傅里葉結果，發(fā)現(xiàn)低速葉輪中振幅最大的為1階擾動，高速葉輪中為2階擾動，通過分析擾動相位變化，估算失速團周向運動速度．但受當時實驗和計算條件限制，僅在壁面周向布置 8個靜壓探針，只能得到 4階模態(tài)，無法指出階數(shù)的物理意義.Joshua等[10]也在軸流葉輪失速實驗中布置 8個周向靜壓探針，并通過分析1階模態(tài)的相位變化計算失速團周向運動速度．通過實驗研究雖然能夠得到更真實的流場信息，但存在較大的測量局限性，使空間傅里葉結果不能細致反映流場信息．

計算流體力學的發(fā)展使葉輪全三維非定常計算開始應用于失速研究中，鞠鵬飛[11]對一軸流葉輪失速過程進行了全通道非定常計算，在周向布置與葉片數(shù)一致的 36個監(jiān)視點，通過空間傅里葉分析，得出葉輪失速流場的周向不均勻度與失速團個數(shù)的關系，認為振幅最大的3階模態(tài)對應3個失速團．

本文通過對 Eckardt[12]離心葉輪全通道非定常計算，在不同葉高、不同流向位置的周向曲線上均勻布置 1024個監(jiān)視點，通過空間傅里葉結果的幅值和相位變化，研究失速流場分布和失速團運動規(guī)律，進行失速信號識別，為葉輪機械失速實時監(jiān)測、主動控制、優(yōu)化設計[13]和設備安全運行提供理論依據(jù)和指導．

1 數(shù)值計算

1.1 計算網(wǎng)格及方法

以 Eckardt離心葉輪為研究對象，選取入口段、葉輪和無葉擴壓器為計算域，幾何參數(shù)可參見Eckardt[12]實驗測量葉輪，葉頂間隙為前緣 0.5mm、尾緣 0.7mm，設計轉速為 14000r/min，設計流量5.31kg/s．通過圖1所示的3種數(shù)目的網(wǎng)格獨立性驗證，考慮計算成本并盡可能顯示失速渦分布，單通道網(wǎng)格數(shù)選擇 150×104：入口段 8×104，無葉擴壓器45×104，葉輪 97×104．其網(wǎng)格見圖 2．將單通道網(wǎng)格繞旋轉軸復制 20個形成完整葉輪網(wǎng)格，網(wǎng)格總數(shù)為3000×104，進行全通道非定常計算．

采用有限體積法數(shù)值求解三維 RANS方程并采用標準k-e湍流模型[14]．邊界條件設置如下：入口軸向進氣，總溫 288K，總壓 101325Pa，定常計算時動靜交界面采用凍結轉子法，以大流量工況為起始點，通過減小出口流量不斷逼近失速工況．以定常計算結果作為初始條件，動靜交界面選擇瞬態(tài)轉靜子模型，給定流量出口，對設計工況點(D)、從實驗失速點開始計算的兩個失速工況點(S1、S2)進行全通道非定常計算，計算總時間為葉輪旋轉 5圈經(jīng)歷時間21.427ms，時間步長為葉片通過周期的 1/20，葉片通過周期(blade pass time，BPT)定義為

圖1 網(wǎng)格獨立性驗證Fig.1 Grid independence test

式中：n為葉輪轉速；Z為葉片數(shù)，對于設計轉速下Eckardt葉輪，1BPT＝0.21427ms．

1.2 氣動性能驗證

為了驗證數(shù)值計算結果的可靠性，下面進行氣動性能的驗證．表1給出了設計轉速下離心葉輪3個工況的氣動性能，與實驗數(shù)據(jù)相比符合較好，最大誤差僅2.21%．

表1 CFX和實驗性能對比Tab.1 Difference of performance between CFX and experiment

圖 3為定常計算和全通道非定常計算總壓特性曲線，計算結果與實驗值[15]整體趨勢符合很好，且非定常計算時均結果更接近實驗值，說明本文的數(shù)值計算方法可以準確反映葉輪流場的總體特征．

圖3 總壓特性曲線Fig.3 Characteristics of total pressure rise

2 流場特征

為了分析葉輪不同工況下非定常流場特征，沿流向等距離做10個截面，在流向截面上給出Q準則的渦系強度，其定義為

式中?和S分別為速度梯度 2階張量的對稱張量(應變率張量)和反對稱張量(渦張量)，其為正值時代表該區(qū)域旋轉部分起主導作用．

圖4為設計工況下流動渦空間分布，每個流道流動一致，在葉輪流向 70%下游開始出現(xiàn)由離心葉輪幾何折轉產(chǎn)生的二次渦(secondary vortex，SV)，且沿流向二次渦逐漸增大，從圖 4(b)的放大圖可以看出，在設計工況下葉輪下游吸力面?zhèn)?SV1、SV2)和壓力面?zhèn)?SV3)存在二次渦，且沿流向到葉輪出口流動渦范圍逐漸增大．

圖4 設計工況流動渦分布Fig.4 Flow vortex distribution under the design condition

圖 5表明在失速工況下，葉輪入口近葉頂區(qū)域80%葉高以上存在 4個周期分布的失速團(stall cells，SC)SC1～SC4，每個失速團影響范圍為 5個流道，因此將葉輪每 5個通道作為一個研究對象，用黃色葉片隔開，并依次沿葉輪旋轉方向(即θ正方向)編號為流道 1～5．從圖 5(c)的放大圖可以看出，受上游失速團的影響，葉輪流向 50%下游就開始出現(xiàn)二次渦，且每個流道二次渦強度不均勻，與設計工況(圖 4(b))相比，失速流道 3和流道 4流動渦影響范圍減小，主要在葉頂間隙區(qū)域，而上游沒有失速團的流道 1和流道 2，流動渦范圍與設計工況差別不大，僅位置發(fā)生改變．

圖5 失速工況流動渦分布Fig.5 Flow vortex distribution under the stall condition

3 空間傅里葉分析

3.1 流場數(shù)據(jù)監(jiān)視點布置

為了較系統(tǒng)全面地分析失速團在徑向、流向和周向的運動規(guī)律，以葉輪入口為基準，在流向(streamwise，st)定義葉輪入口為 0，出口為 100%，在徑向(span，sp)定義輪轂側為 0，輪蓋側為 100%．在90%葉高的 2%、40%和 80% 流向位置的 3條周向曲線，2%流向的 20%、40%、60%、80%、90%、95%和 98%葉高位置的 7條周向曲線上分別均勻布置1024個監(jiān)視點，提取流場數(shù)據(jù)進行空間傅里葉分析；在相對坐標系中的3個流向位置某一通道內(nèi)布置7個不同葉高的監(jiān)視點，如圖 6所示．其命名規(guī)則為流向位置用st表示，葉高位置用sp表示，即st80sp90表示80%流向90%葉高位置，以此類推．

圖6 數(shù)據(jù)監(jiān)視點布置Fig.6 Arrangement of the monitor points

3.2 空間傅里葉分析方法

空間快速傅里葉變換(spatial fast Frourier transform，SFFT)分析可以將周向分布的信號分解為空間坐標的簡諧運動，得出各階模態(tài)的幅值和相位．在葉輪周向曲線上布置N個監(jiān)視點，其觀測數(shù)據(jù)為，傅里葉級數(shù)的展開式表示為

其中各項傅里葉系數(shù)分別為

根據(jù)方程組唯一解條件可知經(jīng)傅里葉變換可以得到m≤N/2階模態(tài)．經(jīng)快速傅里葉變換將N個信號分解為N/2+1個空間模態(tài)對應的復數(shù)形式的傅里葉系數(shù)，即

其幅值為

每個幅值的帶寬為2/N，該信號的實際振幅

相位為

通過分析不同時刻的幅值譜，可以了解某一時刻脈動信號周向分布的不均勻度和周期性特征；通過分析幅值譜隨時間的變化規(guī)律，量化失速流場隨時間的變化規(guī)律；通過對主要模態(tài)的相位研究，計算得到非定常擾動傳播速度和方向．

3.3 空間傅里葉分析結果

采用式(3)～(7)計算得到空間傅里葉結果的振幅，分析設計工況和失速工況葉輪入口不同葉高靜壓空間分布，圖 7所示為 2%流向位置空間快速傅里葉變換結果．圖 7表明，設計工況下只存在 20階及其倍數(shù)階模態(tài)，而失速工況下 20階模態(tài)仍然存在，但明顯小于4階振幅，這與圖5中顯示的4個失速團一致．可見空間傅里葉分析主要模態(tài)的階數(shù)代表該幅值的脈動信號在周向出現(xiàn)的次數(shù)，且幅值越大，說明該階模態(tài)在周向周期出現(xiàn)的規(guī)律越顯著，以此確定葉輪流場是否存在失速擾動和失速團個數(shù)．上述分析說明在設計工況下流場只受葉片通過頻率影響，而失速工況下受4個失速團擾動頻率的影響更強烈，且集中在 80%葉高以上，與上文失速擾動主要在葉輪入口近葉頂區(qū)域一致，且 90%葉高振幅最大，說明該位置失速團擾動最強．

圖7 2%流向位置靜壓空間快速傅里葉變換結果Fig.7 SFFT results of static pressure at 2% streamwise

接下來分析不同流向位置 90%葉高的空間傅里葉結果．從圖 8可以看出，不同流向位置均存在振幅不同的4階模態(tài)和振幅基本相同的20階模態(tài)．在葉輪入口(流向 2%)，11倍葉片通過周期 BPT以前存在明顯的4階模態(tài)，即存在4個失速團，之后4階模態(tài)振幅逐漸減小，4階及其倍數(shù)階模態(tài)開始在3kPa以下波動，說明此時葉輪入口出現(xiàn)更多(8、12或 16個)幅值更小的失速擾動；在葉片折轉位置(流向 40%)，21BPT之前 4階模態(tài)出現(xiàn) 4個局部峰值，說明此時失速團擾動已經(jīng)影響到流向 40%位置，之后 4階擾動開始逐漸增強，在 73BPT時達到最大值，并開始小幅波動，說明此時在該位置存在 4個明顯的失速團；在葉輪下游區(qū)域(流向 80%)，17BPT之前，受上游失速團影響，4階擾動出現(xiàn)了3個局部峰值之后振幅逐漸增大，在 79BPT時達到最大峰值后開始小幅波動，且該位置受葉片通過影響的20階擾動更明顯.

圖8 90%葉高空間快速傅里葉變換振幅分析Fig.8 Amplitudes of the SFFT results at 90% span

綜上所述，可確定葉輪中主要存在 4個失速團，且失速首先發(fā)生在葉輪入口近葉頂區(qū)域，并隨著失速團的脫落、破碎逐漸向下游移動，在 14BPT之后失速團移動到葉輪流向 40%位置，在 17BPT之后失速團影響到葉輪流向 80%位置，60BPT之后葉輪流道內(nèi)4階模態(tài)振幅在6kPa以上小幅波動，此時葉輪進入深度失速工況．

為了進一步分析失速團周向傳播速度，通過失速團擾動相位變化計算其角速度．

葉輪旋轉的角速度計算式為

可計算出葉輪旋轉角速度為1465.33rad/s．

失速團運動角速度為

可得失速團相對葉輪旋轉速度的比值為

圖 9給出了空間傅里葉結果 4階擾動的相位變化，圖中圓點大小代表振幅．

圖9 4階擾動相位Fig.9 Phases of the 4th mode

通過式(8)～(11)可得失速團周向運動速度，在葉輪入口時斜率更大，即運動速度更大，為0.73倍葉輪轉速；30BPT時40%流向位置出現(xiàn)速度為0.62倍葉輪轉速的失速團；40BPT時失速團到達流向 80%位置，60BPT之后，葉輪進入深度失速工況，流道內(nèi)不同流向位置失速團運動速度基本一致，約 0.62～0.64倍葉輪轉速．

3.4 時間傅里葉分析驗證

圖10 靜壓FFT結果Fig.10 FFT results for static pressure

為驗證第 3.3節(jié)空間傅里葉方法分析失速團運動規(guī)律的可行性，對3個流向位置不同葉高監(jiān)視點的靜壓時間信號進行 FFT分析，并用其振幅最大值歸一化，最終結果見圖 10．從圖中可以得出，不同流向位置葉輪失速擾動信號主要集中在80%葉高以上，并且由于監(jiān)視點位于相對坐標系中，沒有出現(xiàn)振幅較大的葉輪旋轉頻率 233.33Hz和葉片通過頻率4666.66Hz．在葉輪入口，主要存在頻率為 1980Hz和2333Hz的非定常擾動，且存在振幅為0.5左右的頻率分別為175Hz、757～1160Hz的擾動，說明入口失速擾動在整個數(shù)值計算時間范圍內(nèi)周期性不明顯，存在速度大小不等的失速擾動．在葉輪折轉位置和葉輪尾緣區(qū)域的失速擾動頻率一致，最大振幅頻率均為784Hz，且存在一個頻率175Hz的擾動．

綜上所述，整個葉輪中均存在175Hz的頻率，應為葉輪旋轉1周后失速團再次經(jīng)過該監(jiān)視點引起，可估算出該失速團速度約為 0.75倍葉輪轉速；葉輪入口振幅最大的頻率應為失速團移動一個通道引起，可估算出該失速團速度約為0.7倍葉輪轉速，且存在很多其他頻率的擾動，筆者推測是由于葉輪入口僅在11BPT之前存在強度較大的 4個失速團，而之后失速團脫落破碎為更多強度較小失速團，使各失速團運動速度不相等．在葉輪流道內(nèi)，失速團擾動頻率大約為 0.66倍葉輪轉速，與空間傅里葉相位分析結果基本一致，說明空間傅里葉分析不僅可以識別失速團個數(shù)，還可以給出失速團運動規(guī)律．

4 結 論

通過對 Eckardt離心葉輪全通道非定常不同葉高、不同流向位置周向靜壓數(shù)據(jù)的空間傅里葉分析，以及與失速渦分布及時間傅里葉結果的對比，得出以下結論．

(1) Eckardt離心葉輪在失速工況下，周向存在4個失速團，主要集中在 80%葉高以上．葉輪失速首先發(fā)生在葉輪入口近葉頂區(qū)域，并逐漸向通道下游移動，當下游存在 4個強度穩(wěn)定的失速團時，表明此時葉輪進入深度失速工況．

(2) 通過對不同流向位置周向空間傅里葉相位分析，得到失速團在周向的運動規(guī)律，失速團首先在葉輪入口以 0.73倍葉輪轉速周向運動，隨著失速團向下游移動，速度逐漸減小為0.63倍葉輪轉速．

(3) 在離心葉輪失速流場中，失速團運動規(guī)律復雜，存在周向、流向和徑向的三維運動．

綜上所述，空間傅里葉分析能準確識別葉輪失速流場中的不均勻性及失速信號周向運動規(guī)律隨時間的變化，可用于實時監(jiān)測和失速預警，對防止失速故障、進行失速主動控制和確保設備安全運行具有重要工程指導價值．

[1] Camp T R，Day I J.A study of spike and modal stall phenomena in a low-speed axial compressor[J].Journal of Turbomachinery，1998，120(3)：393-401.

[2] Day I J.Stall inception in axial flow compressors[J].Journal of Turbomachinery，1993，115：1-9.

[3] 吳艷輝，安光耀，陳智洋，等.跨聲速壓氣機轉子近失速工況非定常流動及相關機理研究[J].推進技術，2016，37(10)：1847-1854.

Wu Yanhui，An Guangyao，Chen Zhiyang，et al.Numerical investigation into unsteady flow and its associated flow mechanism in a transonic compressor rotor at near stall conditions[J].Journal of Propulsion Technology，2016，37(10)：1847-1854(in Chinese).

[4] 程曉斌，聶超群，陳靜宜.軸流壓氣機旋轉失速先兆過程中的頻率階躍現(xiàn)象[J].工程熱物理學報，2000，21(1)：29-33.

Cheng Xiaobin，Nie Chaoqun，Chen Jingyi.The frequency step-up of rotating stall precursors in two axial compressor[J].Journal of Engineering Thermophysics，2000，21(1)：29-33(in Chinese).

[5] 程曉斌，陳靜宜.基于小波變換的離心壓氣機旋轉失速先兆時頻分析[J].工程熱物理學報，2000，21(3)：289-293.

Cheng Xiaobin，Chen Jingyi.The time-frequency analysis of rotating stall inception in the centrifugal compressor based on wavelet[J].Journal of Engineering Thermophysics，2000，21(3)：289-293(in Chinese).

[6] Lin F，Chen J，Li M.Wavelet analysis of rotor-tip disturbances in an axial-flow compressor[J].Journal of Propulsion and Power，2004，20(2)：319-334.

[7] 符 嬈，雷 勇.基于小波分析的壓氣機失速初始擾動信號處理[J].西北工業(yè)大學學報，2010，28(3)：421-424.

Fu Rao，Lei Yong.Applying wavelet analysis to processing precursor signals in compressor to obtain early warning of stall[J].Journal of Northwestern Polytechnical University，2010，28(3)：421-424(in Chinese).

[8] Zhang H D，Yu X J，Liu B J，et al.Using wavelets to study spike-type compressor rotating stall inception[J].Aerospace Science and Technology，2016，58：467-479.

[9] Garnier V H，Epstein A H，Greitzer E M.Rotating waves as a stall inception indication in axial compres-sors[J].Journal of Turbomachinery，1991，113(2)：290-301.

[10] Joshua D C，Scott C M.Analysis of axial compressor stall inception using unsteady casing pressure measurements[J].Journal of Turbomachinery，2013，135(2)：21036-1-21036-12..

[11] 鞠鵬飛.壓氣機旋轉失速數(shù)值模擬研究[D].北京：北京航空航天大學，2001.

Ju Pengfei.Numerical Investigation of the Compressor Rotating Stall[D].Beijing：Beihang University，2001(in Chinese).

[12] Eckardt D.Detailed flow investigations within a high speed centrifugal compressor impeller[J].Journal of Engineering for Gas Turbines & Power，1976，98(98)：391-402.

[13] 劉正先，韓 博，魯業(yè)明.目標優(yōu)化算法在葉片參數(shù)化設計中的應用[J].天津大學學報：自然科學與工程技術版，2017，50(1)：19-27.

Liu Zhengxian，Han Bo，Lu Yeming.Application of the objective optimization algorithm in parametric design of impeller blade[J].Journal of Tianjin University：Science and Technology，2017，50(1)：19-27(in Chinese).

[14] 劉正先，苗永淼.有曲率影響的湍流流場中幾種k-z模型的數(shù)值模擬[J].天津大學學報，2000，33(1)：98-101.

Liu Zhengxian，Miao Yongmiao.Numerical simulation of the effect of streamline curvature on turbulent flow field with threek-zmodels[J].Journal of Tianjin University，2000，33(1)：98-101(in Chinese).

[15] Eckardt D.Instantaneous measurements in the jet-wake discharge flow of a centrifugal compressor impeller[J].Journal of Engineering for Gas Turbines & Power，1975，97(97)：337-345.