過流條件下泵站管道振動響應混沌特性分析

張建偉， 張翌娜， 程夢然， 王立彬

(1.華北水利水電大學 水利學院，鄭州 450046； 2.黃河水利職業技術學院，河南 開封 475004)

泵站管道長期運行過程中，易受自然及人為因素的影響而導致產生不同程度的振動[1]，長時間、無規律的管道振動會導致管系及其附屬系統結構產生松動，嚴重時會造成災難性的破壞[2]，因此，分析泵站管道的振動特性以避免其不利振動具有重要的研究價值。

混沌是強非線性結構在振動中表現出來的一種獨特的力學現象，大多學者認為輸液管道振動屬于弱非線性，主要集中于對水流等激振源所產生的混沌現象的研究。在非線性問題的研究上，對管道本身存在的混沌現象研究較為匱乏。研究發現輸液管道的混沌現象并不僅僅取決于結構非線性的強弱，對某些結構，即使是弱非線性或是線性結構，也會出現混沌運動[3]。Pa?Doussis等[4]研究了在非線性運動約束及定常內流作用下懸臂輸流管道的動力學問題，發現了該系統下的混沌運動。Tang等[5]以增加非線性力的方式得到輸流管道的混沌特性，發現混沌的出現與否受管道內流速影響最大。Sinir[6]研究了勻速輸流微曲管道的非線性振動問題，并在其橫向振動中觀察到周期性和混沌運動。Zhao等[7]分析了熱荷載作用下水流脈動激勵產生的管道振動中的混沌現象，得到頻率響應與水流脈動速度之間的關系。

管道系統混沌特性的研究大多基于油氣壓力管道及特定非線性約束條件的輸流管道數學模型，而對泵站管道系統的混沌特性分析較少，且以往研究成果中大部分僅對振動系統的混沌特性進行分析，并未進一步探究使其產生混沌特性的振動激勵源。本文以某灌區泵站管道工程實測振動響應為研究對象，利用飽和關聯維數、最大Lyapunov指數等混沌特性定量判別指標，分析不同工況下泵站管道各測點振動響應的混沌特性；采用改進的變分模態分解(improved variatronal mode decomposition, IVMD)方法分析典型工況下測點振動響應序列，研究各個本征模函數(intrinsic mode function, IMF)成分及混沌特性，從而得到引起泵站管道產生混沌特性的振動激勵源。

1 基本理論

1.1 混沌特性識別方法

目前，混沌特性識別方法大致分為定性分析方法和定量分析方法。其中定性分析包括相軌直接觀察法、Poincare截面法、功率譜分析法等，具有簡單易行的特點，但僅限于判定系統是否具有混沌特性，無法實現不同工況下的橫向比較；定量分析方法如應用廣泛的飽和關聯維數法、最大Lyapunov指數法等可以通過對參數值大小的比較，反應不同條件下響應振動復雜性及混沌特性的強弱[8]。為了提高結論的可靠性，本文同時選取飽和關聯維數法、最大Lyapunov指數法作為泵站管道振動響應混沌特性識別指標。

1.1.1 飽和關聯維數法

關聯維表征動力系統的緊湊程度，用以體現系統的復雜程度，當某一系統飽和關聯維數存在時，即可判定其具有混沌特征[9]。對于一個m維相空間，其關聯函數C(r)定義為

(1)

式中：C(r)為關聯函數；r為任一正常數；M為相點個數；‖Yi-Yj‖為相空間內i,j兩點間的距離；H(t)為Heaviside函數，定義為

當時間序列具有混沌特征時，對于正值r，關聯函數C(r)與r的關系為

C(r)∝αrD2

(2)

式中：α為常數；D2為關聯維數，可以通過log2C(r)～log2r圖的斜率求得，即

(3)

1.1.2 最大Lyapunov 指數法

Lyapunov指數是根據相軌跡有無擴散運動特征來判別系統的混沌特性。一般地，正的Lyapunov指數代表的方向對吸引子起支撐作用；而負的Lyapunov指數對應的收縮方向，在抵消膨方向的作用后，貢獻吸引子維數的分數部分，故正Lyapunov指數是混沌的一個顯著特征。設某一系統最大Lyapunov指數為λ1，則λ1取值為正即可判定系統存在混沌成分，其值大小反映混沌程度的強弱。

Rosenstein等[10]提出了一種計算λ1的小數據量法，其基本步驟如下：

(4)

式中，p為時間序列的平均周期。

定義該鄰點對j個離散時間步之后的距離為

(5)

對每個j，求出所有i的lndi(j)平均

(6)

式中,q為非零lndi(j)的數目，采用最小二乘法做出回歸直線，則該直線斜率即λ1。

1.2 改進的變分模態分解

VMD是一種多分量自適應信號分解的新方法[11]，與傳統信號分解方法相比，能夠有效避免模態混疊和過分解等缺陷[12]。VMD利用變分約束將給定信號f分解為K個模態函數mk(t)，各個本征模函數IMF的帶寬特定有限，且每個模態函數都分布在中心脈動頻率周圍，變分約束模型表達式為[13]

(7)

式中:{mk}為分解得到的K個IMF分量，{mk}={m1,m2,…,mk}；σ(t)為脈沖函數；{wk}為各IMF分量的中心頻率，{wk}={w1,…,wk}。

為了完成輸入信號f的自適應分解，得到帶寬之和最小的IMFs，引入拓展的Lagrange表達式

L(mk,wk,λ)=

(8)

式中：α為懲罰因子，確保信號重構的準確性；λ(t)為拉格朗日乘子，用來強化約束的嚴謹性；〈·〉為內積運算；δ(t)為脈沖函數。

為解決以上變分約束問題，采用對偶分解和交替方向乘子算法[14]，一直更新mk,wk與λ(t)，使其循環迭代求取式(8)的鞍點，即為式(7)的最優解，模態分量函數mk和中心頻率wk，如式(9)和式(10)所示

(9)

(10)

(11)

VMD在對振動響應序列進行分解時，模態總數K的確定是至關重要的一步，模態參數K的選取極大影響結果的準確性[15]。若K值大于信號分解得到有用成分的個數，則會產生信息疊加；若K值小于信號分解得到有用成分的個數，導致部分有限帶寬的固態模量不能被分解出來。針對K值選取問題提出基于互信息法的改進變分模態分解方法。

互信息(mutual information，MI)反映兩個隨機變量間的關聯性，能更好地辨別相關程度[16]。互信息表示為

I(X,Y)=H(Y)-H(X|Y)

(12)

式中：H(Y)為Y的熵;H(Y|X)為X已知時Y的條件熵。當I(X,Y)=0時，X與Y相互獨立。

計算IVMD分解得到的各IMF分量與原始信號的互信息Ik，并利用式(13)進行歸一化處理，進而判斷各模態分量與原始信號的相關性，即原信號是否完全被分解。

(13)

式中,σi為求得各IMF分量的歸一化互信息值，i=1,2,…,k。當σi低于σ=0.02時，認為IMF分量中已經不含有效的特征信息，原信號已分解完畢，停止運算。

采用互信息法自適應確定K值的具體算法流程如下：

步驟1初始化n=n+1,令K=1；

步驟2K=K+1,執行外層循環；

步驟4令n=n+1,執行內層循環；

步驟5對一切w≥0，根據式(9)和式(10)分別更新mk和wk；

步驟6由式(11)更新λ；

步驟8循環步驟2～7直至設定閥值σ大于原信號f分解得到的各IMF與f的歸一化互信息值σi，即I(f-∑mk,f)<σ，結束循環。

2 管道振動響應混沌特性分析

以某泵站2號壓力管道為研究對象，在主管和兩個支管上共選取6個測點，每個測點的x,y,z三個方向各布置1個891-2型拾振器(檔2)，測點布置如圖1所示。

圖1 管道測點布置圖Fig.1 Layout of pipeline measuring points

原型試驗共選取了4種工況采集管道振動響應，各工況描述、采樣時間及采樣頻率，如表1所示，典型振動響應時程如圖2所示。

表1 管道原型試驗測試工況Tab.1 Test conditions of pipeline prototype

圖2 工況4測點1 z向振動Fig.2 Point 1 of condition 4 of z-axis

首先，對管道振動響應時間序列進行相空間重構，即計算時間延遲τ和嵌入維數m。通過多向比較選取自相關函數法求取τ，選取CAO法求取m，限于篇幅，只列出工況4測點1z向振動響應的τ和m計算過程。

在自相關函數法計算τ的過程中，當函數值下降到初始值的(1-1/e)時，對應的τ值即為所求時間延遲。自相關函數計算結果如圖3所示。

圖3 工況4測點1 z向振動τ計算圖Fig.3 τ calculation of point 1 z-axis vibration under condition 4

τ值確定后，采用CAO法確定嵌入維數m。E1(m)用來確定最小的嵌入維數m，即當E1(m)隨m趨于平穩時對應的m值即為最小嵌入維數。參數E2(m)用于表明時間序列特性，即對任意m，若E2(m)都在 1 附近，則表明該時間序列是隨機序列。具有混沌特性時間序列的E2(m)值逐漸接近 1，可通過此特征定性直觀的判斷一個時間序列是否具有混沌特性。當E1(m)不再隨嵌入維數明顯變化，且E2(m)函數值趨于1時，對應的m值即為最優嵌入維數。由圖4可知，工況4測點1z向振動響應的最優嵌入維數m=11。

圖4 工況4測點1 z向振動m計算圖Fig.4 m calculation of point 1 z-axis vibration under condition 4

其次，選取飽和關聯維數法、G-P算法[17](Grassberger和Procaccia提出的一個從單變量時間序列計算吸引子飽和關聯維數的算法，簡稱G-P算法)及小數據量法分別計算飽和關聯維數和最大Lyapunov指數，并利用兩種指標進行混沌特性分析。

依據前文計算所得最佳時間延遲τ，嵌入維數取m=2,4,6,…,20，按式(3)中關聯函數關系分別作不同m的log2C(r)～log2r對數關系曲線圖，并選取其中的近直線段擬合其斜率，該斜率即為對應嵌入維數下的關聯維數，隨著嵌入維數的增加，關聯維數達到飽和值時即為飽和關聯維數D2。典型測點的飽和關聯維數計算圖，如圖5所示。順箭頭方向分別為m=2,4,6,…,20時的關系曲線。

圖5 工況4測點1 z向振動D2計算圖Fig.5 D2 calculation of point 1 z-axis vibration under condition 4

為揭示關聯維數計算值分布規律，作不同工況下各方向測點振動響應飽和關聯維數D2變化曲線，如圖6所示。

由圖6可知：①整體而言，關聯維數分布在1.156～5.283，說明管道各向振動響應具有混沌特性；②對比三個方向振動響應，主管道軸向(y向)測點關聯維數明顯小于其他兩方向測點，表明管道軸向振動呈現出較低維的混沌吸引子，需要較少的獨立控制變量即可描述其動力系統的變化規律；③同測點下，工況4(4機組、5機組穩定運行)各測點飽和關聯維數D2值大于其他工況，而工況3(4機組單獨關閉)相應D2值小于其他工況，說明機組的運行增加了管道振動的不確定性；④同工況下，靠近機組測點1、測點5及管道岔管處測點4的D2值相對較大，表明泵站管道振動復雜性受機組及流態影響較大。

圖6 各方向測點關聯維數變化曲線Fig.6 Curves of correlation dimension of points in different directions

為對比驗證上述結果，利用最大Lyapunov指數λ1對泵站管道的混沌特性做進一步分析。根據前文求得時間延遲τ和嵌入維數m，采用小數據量法計算最大Lyapunov指數λ1。圖7為典型測點λ1計算圖，可知分離因子y(i)函數值在近線性增長后趨于穩定，通過最小二乘法擬合增長段直線斜率，該值即為λ1。不同振動方向測點λ1值，如圖8所示。

圖7 工況4測點1z向λ1計算圖Fig.7 λ1 calculation of point 1 z-axis vibration under condition 4

由圖8可得： ①各測點不同工況下最大Lyapunov指數λ1取值在0.032 3～0.073 4，均大于零，說明管道實測振動響應具有明顯的混沌特征；②工況4(4機組、5機組穩定運行)各測點λ1相較其他工況更小，而工況3(4機組單獨關閉)λ1值相對偏大，可知隨著工況中兩機組的開啟，最大Lyapunov指數λ1取值隨之減小，表明隨著機組穩定運行，流態逐漸穩定，管道振動的混沌特性逐漸減弱；③同工況下，靠近機組測點1、測點5λ1值較小，而管道岔管處測點4λ1值較其他測點更大，說明岔管處流態使管道振動更具有混沌特性。上述分析與關聯維數D2計算結果互為補充，說明機組運行時流態對泵站管道混沌特性的影響較大。

圖8 各方向測點最大Lyapunov指數變化曲線Fig.8 Curves of largest Lyapunov exponent of points in different directions

3 基于IVMD的多尺度混沌特性分析

泵站管道的振動特點具有不同于一般輸流管道的特殊性，這主要體現在泵站機組對連接管道振動的影響方面，其振源組成主要為管道輸水流動產生的低頻水流脈動以及機組運行產生的葉頻、轉頻及倍頻等。以工況4測點1z向振動響應為例，對其進行頻譜分析如圖9所示。參考文獻[18]，20 Hz,40 Hz及60 Hz分別為機組旋轉倍頻及轉輪葉片振動頻率，而0.5 Hz為泵站管道輸水湍流產生的水流脈動激勵頻率。可知水流脈動產生的管道振動激勵頻帶(0.5 Hz)相對較寬，寬峰功率譜屬于混沌系統的典型特征；機組運行產生的管道振動激勵(20 Hz,40 Hz,60 Hz)對應尖峰功率譜，周期性強，故推測泵站管道的混沌特性主要是輸水管道內湍流脈動引起的，而機組振動掩蓋了泵站管道的混沌特性。

圖9 工況4測點1 z向振動響應頻譜圖Fig.9 Spectrogram of point 1 z-axis vibration response under condition 4

為明晰使管道振動響應具有混沌特性的振動激勵，采用IVMD方法識別引起管道混沌特性的振動激勵源。

利用互信息法確定IVMD的模態參數K=4，管道工況4測點1z向振動響應經IVMD分解獲得4個IMF分量。分解所得IMF分量時程如圖10所示，各IMF分量與原振動序列的歸一化互信息值，如表2所示。

圖10 振動響應IVMD分解IMF時程圖Fig.10 Time histories of IMFs decomposed by IVMD

表2 振動響應各IMF歸一化互信息值Tab.2 IMF normalized mutual information value of vibration response

由表2可知，4個IMF分量歸一化互信息值均大于閥值0.02，滿足分解要求。由圖10可知，原振動響應經IVMD分解后可得4個頻率依次增加的IMF分量，經頻譜分析可得IMF1～IMF4的頻率分別對應原響應頻譜圖中4個主要頻率帶：0.5 Hz，20 Hz，40 Hz，60 Hz。利用飽和關聯維數和最大Lyapunov指數對分解所得各IMF進行混沌特性分析。

典型IMF混沌特征值計算過程如圖11所示，由圖11(a)、圖11(b)可以看出，IMF1飽和關聯維數D2=1.115，最大Lyapunov指數λ1=0.077 4，具有明顯混沌特性；IMF2～IMF4均無法找到D2對數曲線中近線性區域，無混沌特性，限于篇幅僅給出IMF2對數曲線斜率如圖11(c)所示，順箭頭方向分別為m=2,4,6,…，20時的關系曲線。

圖11 典型IMF分量混沌特征值計算圖Fig.11 Calculation of chaotic eigenvalues of typical IMFs

對比分析可得：①代表湍流脈動激勵的IMF1飽和關聯維數1.115為分數維，且最大Lyapunov指數0.077 4大于零，具有明顯混沌特性；代表機組運行振動激勵的IMF2～IMF4無混沌特征，表明機組的運行并不能引起泵站管道振動的混沌特性；②剔除無混沌特性的機組運行振動激勵成分(IMF2～IMF4)后，管道振動響應關聯維數D2由4.985降低至1.115，而最大Lyapunov指數λ1由0.033 5增大至0.077 4，即管道振動的復雜性下降，其混沌特性更為明顯，表明泵站管道振動時，輸水湍流激勵使其振動具有混沌特性，而機組運行產生的振動激勵掩蓋了泵站管道的混沌特性，增加了管道振動的不確定性。

4 結 論

(1)不同工況下泵站管道測點振動響應飽和關聯維數D2均存在，說明泵站管道振動具有混沌特征。

(2)不同工況下各測點最大Lyapunov指數λ1在0.051 3～0.077 4，隨著工況中兩機組的開啟，最大Lyapunov指數λ1取值隨之減小，表明隨著機組穩定運行，流態逐漸穩定，管道振動的混沌特性逐漸減弱。

(3)湍流脈動是引起泵站管道產生混沌特性的主要原因，而機組振動掩蓋了泵站管道的混沌特性，增加了管道振動的不確定性。