基于多重分形的離心壓縮機出口動態壓力非線性特征研究及其在喘振識別中的應用

劉 雁, 高 寬, 何 浩, 馬 凱

(西北工業大學 機電學院，西安 710072)

離心壓縮機在運行過程中，若流量減小到一定程度時，將發生非正常工況下的流動失穩，即喘振。喘振作為壓縮機的一種固有特性，在對壓縮機內部造成巨大沖擊的同時也使整個機組產生強烈振動，可使壓縮機轉子和靜子受交變力作用斷裂，從而引起葉片飛出，甚至發生機毀人亡的事故[1]。因此，尋找一種準確快速的喘振識別及預測方法是保證壓縮機高效運行的必要條件。

近年來，研究人員分別從喘振發生機理、分析預測及防喘振控制等方面進行了研究，旨在避免壓縮機工作時喘振的發生。在喘振機理研究方面，其中一項研究熱點就是壓縮機喘振的模型研究。目前，國內外較為廣泛用于研究的模型是Greitzer模型及其衍生的喘振模型，通過研究其結構參數的改變對系統壓比和流量變化的影響，分析系統的工作狀態，以期達到避免喘振的目的[2-5]。在喘振分析和預測方面，目前采用較多的方法是通過測試壓縮機的出口壓力和壓縮機喉部壓差判斷系統是否工作在安全工作區。

雖然離心壓縮機系統設置了安全工作區，可以使壓縮機遠離喘振區，達到離心壓縮機防喘振的效果，但由于防喘振線的設置需要考慮安全裕度，相應的就要犧牲壓縮機的工作效率。近年來，研究人員發現，當壓縮機發生喘振之前，由于管網中積蓄了高壓氣體，在壓縮機的葉道中會出現氣流的嚴重分離，即旋轉失速。而旋轉失速作為離心壓縮機產生喘振的前奏，可以在壓縮機的出口管道中的氣體動態壓力上體現出來[6]。隨著對離心壓縮機喘振機理研究的深入，研究人員發現在壓縮機的出、入口及葉輪處的壓力信號屬于非線性信號。因此，一些非線性動力學的分析方法如小波圖像處理技術、神經網絡和單重分形等被引入并應用于壓縮機的出、入口的氣體動態壓力測試，目的是實現壓縮機的喘振預測[7-9]。而在防喘振控制方面，目前也較多的采用如模糊控制、神經網絡可變控制以及非線性主動控制等非線性控制方法實現[10-11]。

分形理論作為分析非線性系統中信號的一類重要方法，廣泛地應用于湍流、地球科學、遺傳學、生命科學和金融領域[12]。而離心壓縮機出口動態壓力信號作為一種典型的非線性信號，可以采用分形理論對其進行分析。在此動機下，我們采用了單重分形理論對其進行分析，分析出系統在不同工作狀態下的動態壓力變化的自相關特性。隨著研究的深入，我們發現，當系統從穩態進入到喘振狀態，其多重分形特性表現得更加明顯。基于此背景，我們以800 kW離心壓縮機系統在不同工況下的出口動態壓力特征為分析對象，從非線性角度分析其多重分形特征，研究其精細結構的動力學特征，以期達到早期識別和預測喘振的目的，為壓縮機喘振的主動控制奠定基礎。

1 多重分形

簡單地說，多重分形即為分形對象的分解，分解后每部分均有自己的分形維數。通過多重分形可以描述分形體不同層次的特征，即其局部奇異特性。

1.1 多重分形及其譜參數

假定一單位區間具有單位質量，將其分成若干子區間，每個區間長度為δ，同時為第k個子區間分配一個非負的包含時間點t0的序列μk(k≥1)，則定義μk在時間點t0處的奇異性指數α(t0)為

(1)

奇異性指數α(t)也稱作Hǒlder指數，在分形理論中表征分形體某小區域的分維，又稱為局部分維，其值的大小可以反映小區域生長幾率的大小。

對于一個時間序列，如果式(1)中的極限不存在，則表示序列在t0處的奇異性指數沒有定義；如果α(t)為一常數，則說明該序列的奇異性可以僅用一個全局尺度指數來表征，這種分形特征就是單重分形特征；如果α(t)隨時間t的變化而變化，即其尺度行為特征與時間有關，則表示該序列具有多重分形特性。與單重分形相比，多重分形的概念推廣和拓展了人們對尺度行為的理解和認識，與時間相關的尺度特性能夠反映出局部時間范圍內的不規則現象[13]。

多重分形可用于表示分形體中小區域的分維，如果小區域的數目很大，α(t)則將單重分形的分形指數(又稱為Hurst參數)擴充到多值，這就需要知道α(t)取不同值的概率，才能分析時間序列的特征。于是，得到了一個由不同α(t)組成的無窮序列構成的譜f(α)。f(α)的物理意義是表示相同α值子集的分形維數，即表示α(t)出現的概率，稱為多重分形頻譜(簡稱多重分形譜)。f(α)的取值應在[0，1]之間，通常是上凸形狀(∩)，也可以將多重分形譜理解為指數α(t)的概率分布(即密度)[14]。

1.2 多重分形譜物理量之間的關系

為了計算多重分形譜，通常采用盒子計數法，即把采樣區間內的時間序列作為一維點集，用“盒子”進行覆蓋。用度量為ε的“盒子”，將時間序列分為N個小區間。其中，ε趨近于0。將時間序列在該小區間的概率定義為Pi(ε)整個序列內奇異性指數相同的小區間個數定義為N(ε)，則Pi(ε)和N(ε)的關系式定義為

Pi(ε)∝εα

(2)

N(ε)∝ε-f(α)

(3)

其中，α可以反映時間序列在不同劃分尺寸ε下，物理量分布概率Pi(ε)隨ε變化的各區間的性質，α越大，子集的概率越小(因為ε<1)。

統計物理給出了規則和不規則多重分形譜的計算方法。首先，定義一個配分函數

(4)

式(4)中，q稱為權重因子，多重分形譜通過不同q的取值將分形體分成具有不同層次的區域來研究[14-15]。

若配分函數收斂，則等式

χq(ε)=ετ(q)

(5)

成立，即可以得到

τ(q)=lnχq(ε)/ln(ε)

(6)

其中，τ(q)稱為結構函數。

設N(Pi)是概率為Pi(ε)的小區間的數目，則有

(7)

將式(2)、(3)代入式(7)，并根據式(6)可以得出多重分形譜f(α)和奇異性指數α

f(α)=q·α-τ(q)

(8)

(9)

在分形結構中，單重分形可以看作是多重分形的特殊情況。當系統具有單重分形特征時，α(t)為一常數，用Hurst參數表示，f(α)為1。此時，式(8)可以變換為

1=q·H(q)-τ(q)

(10)

式(10)中，H(q)是Hurst參數，為一常數。此時，τ(q)為H的線性函數[16]。對于具有單重分形特征的時間序列，由于H(q)是常數，τ(q)和q滿足線性關系；而當序列具有多重分形特征時，τ(q)則表現出非線性特性。因此，結構函數的形狀呈現凸形是序列具有多重分形特性的典型標志。

1.3 多重分形譜與離心壓縮機出口動態壓力

奇異性指數α(t)的大小取決于離心壓縮機系統在其動態行為中的信息，不同的α(t)值表示了系統在不同層次的奇異性程度。其中，多重分形譜參數αmin和αmax分別是最大和最小概率分布Pi(ε)的奇異性指數。因此，多重分形譜f(α)可以度量離心壓縮機出口動態壓力分形結構的復雜程度。定義多重分形譜的寬度

Δα=αmax-αmin

(11)

Δα的大小反應出口動態壓力概率分布的均勻度，可以表征在標度不變的情況下，所分析時段內出口動態壓力分布的均勻程度，即出口動態壓力的幅度變化。譜寬度Δα越大表示出口動態壓力分布越不均勻，波動越劇烈。Δα=0則對應完全均勻分布的狀況，即當壓力無任何波動時，多重分形譜將退化成二維空間中的一個點，此時就是均勻分形。

f(αmin)和f(αmax)分別表示最大值、最小值概率子集Pi(ε)的分形維數，可以分別描述出口動態壓力最大值和最小值出現的次數。定義多重分形譜參數Δf(α)

Δf(α)=f(αmin)-f(αmax)

(12)

Δf(α)能夠描述最大、最小概率的數目之比[17]。若多重分形譜的左端低于右端，圖形向右鉤，即Δf(α)<0，表示出口動態壓力序列在這段時間內處于最小值的概率較大；若多重分形譜左端高于右端，圖形向左鉤，即Δf(α)>0，表示出口動態壓力序列在這段時間內的幅值處于最高值的概率較大。

2 出口動態壓力采集及頻譜分析

2.1 采集系統

圖1為用于分析離心壓縮機工作狀態的采集系統，其中離心壓縮機的葉輪葉片數為16片，葉輪的額定轉速為9 620 r/min，馬赫數為0.6，由一臺800 kW的直流電動機拖動。處于室溫狀態的大氣，首先經過空氣過濾器進入工作管路，在風室中等待壓縮，當壓縮機工作時進入離心壓縮機。壓縮后的氣體通過出口管道及直徑分別為250 mm和100 mm的兩個電動防喘振閥，進入排氣消聲器，最終排入大氣。

本文所分析的數據由動態壓力采集模塊完成，數據為800 kW離心壓縮機系統在不同工況下的出口管道的動態壓力。動態數據采集系統采用美國Crystal Instruments公司的動態信號分析儀CoCo-80，采樣頻率為20.48 kHz。傳感器選用美國Kulite公司的動態壓力傳感器。

系統通過調節大防喘振閥和小防喘振閥的開度實現對出口壓力的控制，隨著防喘振閥的開度逐漸減小，系統逐漸由穩定狀態進入喘振狀態。整個實驗過程為635.5 s，本文選取其中180～330 s共150 s的出口動態壓力數據進行分析，將第180 s定義為時間起點。這段時間的數據包含了系統由穩態過渡到喘振的過程。

圖1 采集系統

2.2 離心壓縮機出口動態壓力的波形與頻譜

圖2為離心壓縮機的防喘振閥由打開到關閉時，離心壓縮機出口動態壓力的變化波形。可以看出，系統經歷了穩態、過渡過程，最終進入喘振狀態。在防喘振閥打開時，出口動態壓力處于穩態階段，如圖2中的0～85 s時的狀態，此時出口動態壓力幅值較小，波形表現出隨機特性。隨著防喘振閥開度的減小，壓縮機經過過渡過程進入喘振狀態，出口動態壓力幅值逐漸增大，如圖2中的85～96 s。隨著失穩狀態的加深，系統進入喘振狀態，出口波形開始表現出規律振動，波動周期明顯變長，如圖2中的96～150 s。

(a) 0～150 s

(b) 穩態時(0～10 s)

(c) 過渡過程時(85～95 s)

(d) 喘振時 (100～110 s)

圖3為這三種工作狀態時部分時段的頻譜圖，從圖3可以看出，壓縮機在穩態時出口動態壓力的頻譜較為分散，而隨著喘振的加深，頻譜越來越向低頻處集中。進一步，從圖3(c)可以看出，離心壓縮機喘振時的頻譜主要集中在0～0.2 Hz，能量高度集中，對系統的破壞力很強，極易發生事故。

(a) 穩態時(0～10 s)

(b) 過渡過程時(85～95 s)

(c) 喘振時(100～110 s)

3 離心壓縮機出口動態壓力的多重分形特征

3.1 離心壓縮機出口動態壓力的結構函數

圖4所示為離心壓縮機出口動態壓力的結構函數τ(q)的隨權重因子q的變化曲線。這里選擇離心壓縮機在穩態、過渡過程和喘振三種狀態的時段，每個時段選取5 s，以1 s為分析單元進行分析。從圖4可以看出，在穩態時，結構函數隨權重因子q的變化曲線十分平直，當系統進入非穩定狀態時，曲線出現曲度，而隨著系統經過過渡過程進入喘振狀態，曲線的曲度加深。說明系統在穩態時的單重分形特征明顯高于非穩態。系統從過渡過程進入喘振狀態，系統的多重分形特征相應地變得十分明顯。因此，我們進一步從多重分形的角度研究系統的非線性特征。

(a) 0～5 s

(b) 85～90 s

(c) 100～105 s

3.2 離心壓縮機出口動態壓力的多重分形譜

同樣選取穩態、過渡過程和喘振三種狀態，6個時段分析，每個時段選取5 s，以1 s為分析單元畫出多重分形譜如圖5所示。從圖5可以看出，幾乎所有時段的多重分形譜圖都表現出向左鉤的圖形特征，即Δf(α)>0，表示這些時段出現最大值的概率大于出現最小值的概率。在穩態時，Δα很小，Δf(α)也較小，表明在此階段離心壓縮機的出口動態壓力分布較均勻，系統發生喘振的幾率較小。而隨著系統由穩態經過過渡過程進入到喘振狀態，Δα和Δf(α)都在增大，說明系統動態壓力分布變得不均勻。而且在過渡過程時，Δα和Δf(α)的變化量都較大，說明此階段系統的出口動態壓力的分布十分不均勻，而且最大值、最小值出現的概率也不穩定。系統處于不穩定狀態，極易引發失穩現象，從而發生喘振。

(a) 0～5 s

(b) 80～85 s

(c) 85～90 s

(d) 90～95 s

(e) 105～110 s

(f) 110～115 s

3.3 權重因子的選擇與多重分形譜

多重分形是通過不同q的取值將分形體分成具有不同層次的區域進行研究的。理論上q值的范圍越大越好(-∞

圖6所示為系統在穩態、過渡過程和喘振三種狀態中各選取1 s畫出的權重因子q與多重分形譜的關系曲線。其中，40～41 s和140～141 s權重因子的范圍為-80

(13)

其中，fmax為1。由于對于同一個|q|有2個k值，在實際計算時可以選擇較大的值，然后通過設定k的閾值，如不超過0.1%，即可選取q的范圍。按照這種選擇方法，同時結合計算機運算工作量，在保證不發生溢出錯誤的情況下，0～94 s和106～150 s的壓力序列選取|q|max=80，而94～106 s的壓力序列選取|q|max=20。即可以保證全面分析動態壓力的概率分布。

(a) 40～41 s

(b) 94～95 s

(c) 140～141 s

(a) 40～41 s

(b) 94～95 s

(c) 140～141 s

(a) 40～41 s

(b) 94～95 s

(c) 140～141 s

3.4 離心壓縮機出口動態壓力的多重分形譜參數

圖9為離心壓縮機出口動態壓力的多重分形譜參數隨時間的變化曲線，從圖中可以看出離心壓縮機從穩態到喘振，多重分形譜的各參數都發生了明顯的變化。通過研究各參數的變化規律，得出以下結論：

(1) 在0～85 s時，Δα接近于零，αmin最大而αmax最小，表明在此階段系統出口動態壓力分布均勻，波動較為平穩，出現喘振的幾率較小。

(2) 在85～96 s時，αmax的值迅速增大，而αmin的值則迅速減小。這使得Δα迅速增加，系統出口動態壓力的分布變得不均勻。f(αmax)在這一階段波動的十分劇烈，表明壓力波動出現最小值的概率變化較大，而f(αmin)波動變化不明顯，表明壓力波動出現最大值的概率變化不大。Δf(α)在此階段出現了波動向上的變化趨勢，表明壓力波動出現最大值的概率與最小值的概率相比在增加，系統開始變得不穩定。

(3) 在96～150 s時，離心壓縮機處于喘振狀態，此階段Δα的寬度遠高于穩態時，表明壓縮機的出口動態壓力的分布不均勻，壓力波動較為劇烈。多重分形譜參數在喘振階段多處于較大的波動狀態，表明系統的工作狀態不穩定。與穩定狀態相比，f(αmax)在這一階段的值小于穩態時的值，表示此階段出現最小值的概率小于穩態時的概率。f(αmin)在喘振階段的波動和穩態時的波動量相比變化不大，即出現最大值的概率和在穩態時相當。因此，Δf(α)總體值也較大，表明出現最大值的概率與出現最小值的概率比值在增加。但是，f(αmin)在喘振初期出現了短暫的峰值平臺，表明此時段出現最大值的概率最大，而且Δf(α)的峰值也出現在此短暫的時段，即出現最大值的概率與最小值的概率比值的最大值，表明此時段為出口動態壓力波動最劇烈的時段。

進一步分析0～150 s的多重分形譜參數的變化規律。可以看出，所有多重分形譜參數在離心壓縮機的過渡狀態均出現了突變特征。其中，Δα、αmax和αmin在喘振時的波動量遠大于穩態時的波動量，系統壓力分布的不均勻程度較大。對于參數Δf(α)、f(αmax)和f(αmin)來說，在喘振時出現最小值的概率有所下降，而且在喘振初期，參數Δf(α)和f(αmin)出現了明顯的高峰平臺和峰值特征。因此，這幾個參數可以進一步用于離心壓縮機的初始喘振的預測。

(a) Δα

(b) αmin

(c) αmax

(d) Δf(α)

(e) f(αmin)

(f) f(αmax)

4 結 論

以800 kW離心壓縮機不同工作狀態下的出口動態壓力為研究對象，選取包含穩態、過渡過程和喘振狀態150 s的數據進行分析，從多重分形的角度對其分形特征進行研究。離心壓縮機出口動態壓力結構函數的曲線在穩態時表現出線性特征，而在非穩態時表現出多重分形特征。表明離心壓縮機出口動態壓力的多重分形譜與離心壓縮機的工作狀態存在較為明顯的關聯關系。多重分形譜特征的分析結果顯示，多重分形譜圖在穩態時最窄，而經過過渡過程進入喘振，寬度明顯增加，動態壓力分布不均勻，并且在過渡過程時譜寬度的變化量較大，極易發生喘振。出口動態壓力的多重分形譜參數在離心壓縮機的過渡狀態均出現了突變特征，尤其是Δα、αmin、αmax出現了明顯的波動，而且這些參數在喘振時的波動量遠大于穩態時的波動量。因此，可以采用這些參數判斷系統進入過渡過程狀態。而且，參數f(αmin)和Δf(α)在喘振初期也出現了明顯的局部特征，兩個參數可以用于判斷系統是否進入喘振。因此，可以采用出口動態壓力的多重分形譜參數實現離心壓縮機初期喘振的識別和預測，并為進一步完成離心壓縮機喘振的主動控制奠定基礎。