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

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

泵站管道長期運(yùn)行過程中，易受自然及人為因素的影響而導(dǎo)致產(chǎn)生不同程度的振動[1]，長時間、無規(guī)律的管道振動會導(dǎo)致管系及其附屬系統(tǒng)結(jié)構(gòu)產(chǎn)生松動，嚴(yán)重時會造成災(zāi)難性的破壞[2]，因此，分析泵站管道的振動特性以避免其不利振動具有重要的研究價值。

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

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

1 基本理論

1.1 混沌特性識別方法

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

1.1.1 飽和關(guān)聯(lián)維數(shù)法

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

(1)

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

當(dāng)時間序列具有混沌特征時，對于正值r，關(guān)聯(lián)函數(shù)C(r)與r的關(guān)系為

C(r)∝αrD2

(2)

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

(3)

1.1.2 最大Lyapunov 指數(shù)法

Lyapunov指數(shù)是根據(jù)相軌跡有無擴(kuò)散運(yùn)動特征來判別系統(tǒng)的混沌特性。一般地，正的Lyapunov指數(shù)代表的方向?qū)ξ悠鹬巫饔?；而?fù)的Lyapunov指數(shù)對應(yīng)的收縮方向，在抵消膨方向的作用后，貢獻(xiàn)吸引子維數(shù)的分?jǐn)?shù)部分，故正Lyapunov指數(shù)是混沌的一個顯著特征。設(shè)某一系統(tǒng)最大Lyapunov指數(shù)為λ1，則λ1取值為正即可判定系統(tǒng)存在混沌成分，其值大小反映混沌程度的強(qiáng)弱。

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

(4)

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

定義該鄰點(diǎn)對j個離散時間步之后的距離為

(5)

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

(6)

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

1.2 改進(jìn)的變分模態(tài)分解

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

(7)

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

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

L(mk,wk,λ)=

(8)

式中：α為懲罰因子，確保信號重構(gòu)的準(zhǔn)確性；λ(t)為拉格朗日乘子，用來強(qiáng)化約束的嚴(yán)謹(jǐn)性；〈·〉為內(nèi)積運(yùn)算；δ(t)為脈沖函數(shù)。

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

(9)

(10)

(11)

VMD在對振動響應(yīng)序列進(jìn)行分解時，模態(tài)總數(shù)K的確定是至關(guān)重要的一步，模態(tài)參數(shù)K的選取極大影響結(jié)果的準(zhǔn)確性[15]。若K值大于信號分解得到有用成分的個數(shù)，則會產(chǎn)生信息疊加；若K值小于信號分解得到有用成分的個數(shù)，導(dǎo)致部分有限帶寬的固態(tài)模量不能被分解出來。針對K值選取問題提出基于互信息法的改進(jìn)變分模態(tài)分解方法。

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

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

(12)

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

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

(13)

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

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

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

步驟2K=K+1,執(zhí)行外層循環(huán)；

步驟4令n=n+1,執(zhí)行內(nèi)層循環(huán)；

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

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

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

2 管道振動響應(yīng)混沌特性分析

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

圖1 管道測點(diǎn)布置圖Fig.1 Layout of pipeline measuring points

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

表1 管道原型試驗(yàn)測試工況Tab.1 Test conditions of pipeline prototype

圖2 工況4測點(diǎn)1 z向振動Fig.2 Point 1 of condition 4 of z-axis

首先，對管道振動響應(yīng)時間序列進(jìn)行相空間重構(gòu)，即計算時間延遲τ和嵌入維數(shù)m。通過多向比較選取自相關(guān)函數(shù)法求取τ，選取CAO法求取m，限于篇幅，只列出工況4測點(diǎn)1z向振動響應(yīng)的τ和m計算過程。

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

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

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

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

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

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

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

為揭示關(guān)聯(lián)維數(shù)計算值分布規(guī)律，作不同工況下各方向測點(diǎn)振動響應(yīng)飽和關(guān)聯(lián)維數(shù)D2變化曲線，如圖6所示。

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

圖6 各方向測點(diǎn)關(guān)聯(lián)維數(shù)變化曲線Fig.6 Curves of correlation dimension of points in different directions

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

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

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

圖8 各方向測點(diǎn)最大Lyapunov指數(shù)變化曲線Fig.8 Curves of largest Lyapunov exponent of points in different directions

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

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

圖9 工況4測點(diǎn)1 z向振動響應(yīng)頻譜圖Fig.9 Spectrogram of point 1 z-axis vibration response under condition 4

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

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

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

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

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

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

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

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

4 結(jié) 論

(1)不同工況下泵站管道測點(diǎn)振動響應(yīng)飽和關(guān)聯(lián)維數(shù)D2均存在，說明泵站管道振動具有混沌特征。

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

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