許貝貝，陳帝伊，李歡歡，閆懂林

(1.西北農(nóng)林科技大學(xué) 水利水電科學(xué)研究院，西安 712100；2.西北農(nóng)林科技大學(xué) 中國旱區(qū)節(jié)水農(nóng)業(yè)研究院，西安 712100)

近年來，歐洲許多國家和俄羅斯等投入大量財力物力對服役多年水電站重建或改造以提高其服役年限和運(yùn)行性能[1]。顯然，隨時間流失，水電站各獨(dú)立子系統(tǒng)必然出現(xiàn)設(shè)備老化、機(jī)組振動及噪聲變大、絕緣性變差等現(xiàn)象。從非線性動力學(xué)角度講，可認(rèn)為由各獨(dú)立子系統(tǒng)結(jié)構(gòu)參數(shù)、性能參數(shù)等異變性引起水電站軸系振動擺度惡化，即可等價于：水電站系統(tǒng)非線性動力演變已成為水電站軸系振動擺度惡化主要原因。

沖擊式水電站由壓力引水管道、水輪發(fā)電機(jī)組及其輔助設(shè)備構(gòu)成的具有復(fù)雜水機(jī)電磁耦合關(guān)系的復(fù)雜非線性系統(tǒng)[2-6]，其建模理論與方法的發(fā)展經(jīng)歷了從單一領(lǐng)域獨(dú)立分塊建模到多領(lǐng)域統(tǒng)一建模的發(fā)展階段：①各組成模塊動力穩(wěn)定性研究迅速發(fā)展，形成較為成熟的單一領(lǐng)域分塊建模理論與方法[7-9]；在壓力引水管道方面，基于計算流體力學(xué)理論，重點(diǎn)研究引水管道內(nèi)部在恒定流狀態(tài)→非恒定流狀態(tài)→恒定流狀態(tài)過渡過程下壓力交替升降過程，即水擊過程。解析法、圖解法和數(shù)值解法是探求該過渡過程下壓力變化特征主要求解方法。目前，數(shù)值解法作為主流計算水擊方法，已被廣泛驗證其正確性和實用性。在水輪發(fā)電機(jī)組方面，針對水輪機(jī)過渡過程動態(tài)特性，可分為外特性法和內(nèi)特性法。外特性法，往往采用水力機(jī)組全特性曲線并聯(lián)立水輪發(fā)電機(jī)組運(yùn)動方程式，通過推求過渡過程線求解水力機(jī)組各工況瞬變規(guī)律。內(nèi)特性法，采用水輪機(jī)構(gòu)造參數(shù)推求暫態(tài)過程下水輪機(jī)動態(tài)力矩和水頭非線性方程，求解水力機(jī)組暫態(tài)工況下動態(tài)力矩和水頭變化特征。針對水輪發(fā)電機(jī)組軸系，重點(diǎn)研究機(jī)組在水機(jī)電磁耦合作用下軸心軌跡等參數(shù)瞬態(tài)變化規(guī)律。目前，上述計算方法已被國內(nèi)外大部分學(xué)者采用。在機(jī)組輔助設(shè)備方面，分為機(jī)械和和電氣兩種類型，其目的為有效提高電網(wǎng)安全可靠性，其模型結(jié)構(gòu)可見參考文獻(xiàn)[6]。目前，兩種類型調(diào)速器成為國內(nèi)外各類水電站重要組成部分。②國內(nèi)外水力發(fā)電系統(tǒng)多領(lǐng)域統(tǒng)一建模理論與方法近年來有明顯進(jìn)展[10-19]；水力發(fā)電系統(tǒng)受水力、機(jī)械和電磁耦合因素影響，其統(tǒng)一建模體現(xiàn)為多領(lǐng)域及各學(xué)科交叉融合趨勢。

20世紀(jì)70年代前，分?jǐn)?shù)階微積分理論僅僅是一個數(shù)學(xué)領(lǐng)域純理論問題，其詳細(xì)發(fā)展歷史過程請參閱具體參考文獻(xiàn)[20-21]。1965年，美國耶魯大學(xué)Prof. Mandelbrot提出分形概念。自此以后，分?jǐn)?shù)階微積分理論被廣泛應(yīng)用于電氣工程、機(jī)械學(xué)、湍流、物理及控制理論、黏彈性動力學(xué)和轉(zhuǎn)子動力學(xué)等，尤其是分?jǐn)?shù)階阻尼系統(tǒng)非線性動力學(xué)研究，因為分?jǐn)?shù)階微積分強(qiáng)大記憶特性使得其描述復(fù)雜系統(tǒng)非線性動力演化過程更加真實、簡潔。Cao等[22]嘗試研究分?jǐn)?shù)階阻尼碰摩轉(zhuǎn)子軸承系統(tǒng)的非線性動力學(xué)特征。Yan等[23]通過改進(jìn)傳統(tǒng)整數(shù)階黏彈性模型，提出履帶式車輛分?jǐn)?shù)階模型并從實驗上驗證其有效性。劉鋮等[24]針對雙饋風(fēng)力發(fā)電系統(tǒng)設(shè)計分?jǐn)?shù)階自抗擾廣域阻尼控制器用以增強(qiáng)互聯(lián)網(wǎng)電網(wǎng)阻尼控制能力并驗證方法的有效性?？梢?，考慮分?jǐn)?shù)階阻尼是水電站系統(tǒng)非線性動力演化過程研究所必要的。

綜上所述，本文針對沖擊式水電站強(qiáng)非線性和多子系統(tǒng)耦合性，引入分?jǐn)?shù)階阻尼項，在考慮壓力引水管道、水輪發(fā)電機(jī)組和調(diào)速設(shè)備基礎(chǔ)上，建立沖擊式水電站分?jǐn)?shù)階動力學(xué)模型，研究分?jǐn)?shù)階階次、勵磁電流、上導(dǎo)軸承剛度、下導(dǎo)軸承剛度、水導(dǎo)軸承剛度和機(jī)組轉(zhuǎn)速參數(shù)值變化對機(jī)組軸心偏移、旋轉(zhuǎn)速度、噴針行程、水輪機(jī)流量和水頭的影響規(guī)律。

1 考慮分?jǐn)?shù)階阻尼沖擊式水電站模型

1.1 噴嘴流動模型

沖擊式水輪機(jī)是按動量定理工作的水力原動機(jī)，壓力鋼管內(nèi)有壓流經(jīng)噴嘴噴入大氣之中，形成射流，進(jìn)而推動水輪機(jī)轉(zhuǎn)輪旋轉(zhuǎn)。本文通過對噴嘴模型進(jìn)行簡化，獲得結(jié)構(gòu)簡圖如圖1所示。

圖1 沖擊式水輪機(jī)噴嘴結(jié)構(gòu)模型

噴嘴流動面積可寫為[8]

(1)

根據(jù)式(1)，沖擊式水輪機(jī)流量可表示為

(2)

式中：Q為噴嘴出射流量；V為射流速率；Sj為出射面積；κ為射流系數(shù)；Z為噴嘴數(shù)目；g為重力加速度；H為水輪機(jī)水頭。根據(jù)式(1)和(2)可得

(3)

式中：q為流量Q相對偏差值；h為水頭H相對偏差值。下角標(biāo)a和b分別表示噴針位置a和位置b。式(3)兩邊同時對時間t求導(dǎo)可得

(4)

1.2 壓力引水管道模型

引水管道水力損失為[7]

hq=h0-h

(5)

式中：h0為靜水頭；h為水輪機(jī)水頭。根據(jù)孔口出流原理，壓力引水管道出口處流量可表示為

(6)

式中：xr為額定負(fù)荷下主接力器位移標(biāo)幺值。

式(6)兩邊同時對時間t求導(dǎo)可得

(7)

根據(jù)式(5)和(6)，壓力引水管道模型可寫為

(8)

式中：Te為彈性時間常數(shù)；Zn為管道水力浪涌阻抗規(guī)格化值。

1.3 發(fā)電機(jī)負(fù)載與調(diào)速器模型

發(fā)電機(jī)負(fù)載動態(tài)特性可描述為[7]

(9)

式中:Ta為水輪機(jī)轉(zhuǎn)動部分慣性時間常數(shù)；Tb是負(fù)載轉(zhuǎn)動部分慣性時間常數(shù)，Tab=Ta+Tb；ω為機(jī)組轉(zhuǎn)動角速度；D為發(fā)電機(jī)阻尼系數(shù)。水輪機(jī)力矩mt采用IEEE Group提出的簡單非線性模型：mt=Ath(q-qnl)-Dtyω。

液壓隨動系統(tǒng)動態(tài)特性可表示為

(10)

式中:Ty為延遲時間常數(shù);kp、ki和kd分別為比例、積分和微分調(diào)節(jié)增益;s為指令信號;y為導(dǎo)葉開度。

1.4 水輪發(fā)電機(jī)組軸系模型

通過考慮水力不平衡力重新構(gòu)建水輪機(jī)動力矩與發(fā)電機(jī)角速度表達(dá)式，從而連接水輪機(jī)調(diào)節(jié)系統(tǒng)和水輪發(fā)電機(jī)組軸系統(tǒng)模型。發(fā)電機(jī)轉(zhuǎn)子在x方向和y方向受到分?jǐn)?shù)階阻尼力可描述為

(11)

式中:c為阻尼系數(shù);(x01,y01)為轉(zhuǎn)子軸心在x方向和y方向偏移量;α為分?jǐn)?shù)階階次。

機(jī)組軸系拉格朗日函數(shù)可表示為

L=T-U

(12)

式中:T為軸系動能;U為軸系勢能。

根據(jù)公式(12)，系統(tǒng)方程可表示為

(13)

式中：MgB為發(fā)電機(jī)額定轉(zhuǎn)矩；Fx-ump和Fy-ump為x方向和y方向不平衡磁拉力；Fx和Fy為非線性油膜力；Fxf和Fyf為分?jǐn)?shù)階阻尼力；Fx-rub和Fy-rub為x和y方向碰摩力。不平衡磁拉力和非線性油膜力解析表達(dá)式見參考文獻(xiàn)[11]，碰摩力見參考文獻(xiàn)[15]。根據(jù)式(13)，系統(tǒng)方程為

(14)

式中：變量x1、x2、x3和x4為中間變量參數(shù)，并無實際物理意義(詳見參考文獻(xiàn)[5, 25])；變量vx01和vy01分別為轉(zhuǎn)子軸心在x方向和y方向變化速率；變量m1和m2分別為發(fā)電機(jī)轉(zhuǎn)子和水輪機(jī)轉(zhuǎn)輪質(zhì)量；e1和e2分別為發(fā)電機(jī)轉(zhuǎn)子和水輪機(jī)轉(zhuǎn)輪質(zhì)量偏心距；r為發(fā)電機(jī)轉(zhuǎn)子和水輪機(jī)轉(zhuǎn)輪形心在水平方向距離；θ和φ分別為發(fā)電機(jī)轉(zhuǎn)子和水輪機(jī)轉(zhuǎn)輪轉(zhuǎn)角；變量k1、k2和k3分別為上導(dǎo)軸承剛度、下導(dǎo)軸承剛度和水導(dǎo)軸承剛度值。

論文中包含四個子系統(tǒng)：壓力引水管道子系統(tǒng)、水輪機(jī)力矩及發(fā)電機(jī)子系統(tǒng)、水輪發(fā)電機(jī)組軸系子系統(tǒng)和調(diào)速子系統(tǒng)。其中，壓力管道參數(shù)和噴嘴流量的耦合可用于水輪機(jī)入口處的流量q和水頭h，其詳細(xì)關(guān)系見公式(1)；壓力引水管道參數(shù)(公式(1))與軸系振動方程(公式(3))通過水輪機(jī)力矩及發(fā)電機(jī)參數(shù)(公式(2))建立耦合關(guān)系，耦聯(lián)參數(shù)包括流量q、水頭h、導(dǎo)葉開度y、力矩mt和角速度ω。調(diào)速子系統(tǒng)通過發(fā)電機(jī)角速度與壓力引水管道參數(shù)建立耦合關(guān)系。四個參數(shù)詳細(xì)耦合關(guān)系如圖2所示，四個子系統(tǒng)模型參數(shù)耦合見公式(1)～(4)。

圖2 引水管道子系統(tǒng)、水輪機(jī)力矩及發(fā)電機(jī)子系統(tǒng)、水輪發(fā)電機(jī)組軸系子系統(tǒng)和調(diào)速子系統(tǒng)耦合關(guān)系

Fig.2 Coupling relationship of pressure penstock, hydro-turbine torque, generator, shafting, and governor

2 非線性動力演化過程

分岔圖是系統(tǒng)狀態(tài)變量與異變參數(shù)構(gòu)成二維空間極限集隨參數(shù)變化圖形，反應(yīng)系統(tǒng)隨參數(shù)變化非線性動力演化情況，故可用來表征沖擊式水電站軸系振動擺度。本文采用冶勒水電站六噴嘴水輪發(fā)電機(jī)組參數(shù)：算例參數(shù)來自于冶勒水電站六噴嘴機(jī)組參數(shù)。其中，參數(shù)HB=580 m,Qb=23.5 m3/s,Z=6;κ=0.985;Zn=1.501 s;Te=0.515 5 s;g=9.81 m/s2;se=0.23 m;d1max=0.085 8 m來自于參考文獻(xiàn)[26]。參數(shù)m1=1.2×104kg;m2=2.1×105kg;e1=3×10-3m;Mgb=2.25×108N·m;J1=7.9×107;J2=3.5×107s;Tab=10 s來自于參考文獻(xiàn)[27]。根據(jù)參考文獻(xiàn)[28]，發(fā)電機(jī)阻尼系數(shù)D一般視為常數(shù)，取值在0～3范圍內(nèi)比較合理，故本文中取值D=0.5；s為參考輸入值，s=0。發(fā)電機(jī)轉(zhuǎn)子和水輪機(jī)轉(zhuǎn)輪形心在水平方向距離為估算值，r=0.000 02 m。PID三參數(shù)kp=1 s;ki=10 s;kd=3.5 s;均為估算值。轉(zhuǎn)輪質(zhì)量偏心也為估算值，e2=0.5×10-3m。變量k1、k2和k3分別為上導(dǎo)軸承剛度、下導(dǎo)軸承剛度和水導(dǎo)軸承剛度值。其值對軸系振動擺度值影響很大，故文獻(xiàn)中在模型耦合壓力引水管道參數(shù)后，著重研究其值對軸系振動擺度的影響規(guī)律。取值變化范圍(107，108)。在超過108值時，軸系振動擺度發(fā)散，系統(tǒng)已不可控。系統(tǒng)初值為[x1,x2,x3,x4,x,ω,h,q,x01,y01,vx,vy,φ]T=[0.001,0.001,0.001,0.000 4,0.004,-7×10-9,0.001 34,0.001 021,0.000 001,0.000 001,0,0,0]T。

2.1 勵磁電流影響

本節(jié)研究不同分?jǐn)?shù)階階次下勵磁電流值異變對系統(tǒng)動力演化過程影響，如圖3所示。

圖3 不同分?jǐn)?shù)階階次下勵磁電流值異變對軸系振動動力演化影響

Fig.3 The shaft vibration versus the excitation current with different level of fractional order

從圖3可總結(jié)以下四方面結(jié)果：①不同分?jǐn)?shù)階階次下，系統(tǒng)動力失穩(wěn)所對應(yīng)的閾值不同，但最終都會經(jīng)隨機(jī)振動態(tài)進(jìn)入到不可控態(tài)；②分?jǐn)?shù)階階次取值不影響系統(tǒng)動力演化方式，演化方式均是通過振動態(tài)→隨機(jī)振動態(tài)→不可控態(tài)；③勵磁電流對系統(tǒng)動態(tài)動力演化過程影響很大，但從電站運(yùn)行角度，在勵磁電流異變值不超過閾值情況下，水電站軸系振動擺度均滿足電站運(yùn)行要求；④隨勵磁電流逐漸增大，軸心偏移量增大，機(jī)組振動增強(qiáng)，且影響程度很大，極易引發(fā)碰摩事故。振動態(tài)進(jìn)入隨機(jī)振動態(tài)閾值(定義閾值A(chǔ))隨分?jǐn)?shù)階階次變化情況如圖4(a)所示。隨機(jī)振動態(tài)進(jìn)入不可調(diào)態(tài)閾值(定義閾值B)隨分?jǐn)?shù)階階次變化情況如圖4(b)所示。

觀察圖4，分?jǐn)?shù)階階次取值對閾值A(chǔ)影響不大，基本在740 A左右變化；分?jǐn)?shù)階階次取值對閾值B影響很大，整個變化范圍可延展至(750，2 750)A，且具有先突然增大后逐漸減小變化規(guī)律。

2.2 上導(dǎo)軸承剛度影響

本節(jié)研究不同分?jǐn)?shù)階階次下上導(dǎo)軸承剛度異變對軸系振動動力演化過程影響，如圖5所示。

觀察圖5，隨上導(dǎo)軸承剛度減小，機(jī)組振動在維持一段時間等幅振蕩后，最終通過隨機(jī)振動態(tài)進(jìn)入不可控態(tài)。當(dāng)α=1.00，閾值A(chǔ)所對應(yīng)上導(dǎo)軸承剛度值最大；當(dāng)分?jǐn)?shù)階階次增加或減小，閾值A(chǔ)值均會減小。當(dāng)α取值小于0.9或大于1.1，系統(tǒng)動力演化方式也發(fā)生了變化，出現(xiàn)由穩(wěn)定直線過渡到連續(xù)波動狀態(tài)并逐漸進(jìn)入隨機(jī)振動態(tài)趨勢。

(a) 周期振動進(jìn)入點(diǎn)移動情況

(b) 隨機(jī)振動進(jìn)入點(diǎn)移動情況

圖4 不同分?jǐn)?shù)階階次下系統(tǒng)進(jìn)入振動周期態(tài)和隨機(jī)振動態(tài)所對應(yīng)勵磁電流值的移動情況

Fig.4 The bifurcation points versus the excitation current with different level of fractional order

圖5 不同分?jǐn)?shù)階階次下上導(dǎo)軸承剛度異變對系統(tǒng)動力演化過程影響

Fig.5 The shaft vibration versus the stiffness of upper guide bearing with different level of fractional order

2.3 下導(dǎo)軸承剛度影響

本節(jié)研究不同分?jǐn)?shù)階階次下，下導(dǎo)軸承剛度異變對系統(tǒng)動力演化過程影響，如圖6所示。

觀察圖6，隨下導(dǎo)軸承剛度增加，改變分?jǐn)?shù)階階次對系統(tǒng)動力演化過程影響很大，但均以Hopf分岔形式失去穩(wěn)定。具體來說，α分別取0.75、1.20和1.25時，系統(tǒng)失穩(wěn)后均直接進(jìn)入隨機(jī)振動態(tài)；α分別取0.80、0.85和1.15時，系統(tǒng)通過周期振蕩直接進(jìn)入大幅度隨機(jī)振動態(tài)。當(dāng)α分別取0.90、0.95、1.00、1.05和1.10時，系統(tǒng)通過周期振蕩先進(jìn)入小幅度隨機(jī)振動狀態(tài)，而后進(jìn)入大幅隨機(jī)振動態(tài)。

圖6 不同分?jǐn)?shù)階階次下下導(dǎo)軸承剛度異變對系統(tǒng)動力演化過程影響

Fig.6 The shaft vibration versus the stiffness of lower guide bearing with different level of fractional order

2.4 水導(dǎo)軸承剛度影響

本節(jié)研究不同分?jǐn)?shù)階階次下，水導(dǎo)軸承剛度異變對系統(tǒng)動力演化過程影響，如圖7所示。

圖7 不同分?jǐn)?shù)階階次下水導(dǎo)軸承剛度異變對系統(tǒng)動力演化過程影響

Fig.7 The shaft vibration versus the stiffness of water guide bearing with different level of fractional order

水導(dǎo)軸承剛度變化和下導(dǎo)軸承剛度對系統(tǒng)軸系振動擺度影響相似，這里不進(jìn)行詳細(xì)描述，主要介紹不同之處。在α=0.95時，隨水導(dǎo)軸承剛度增加，先出現(xiàn)Hopf分岔，接著進(jìn)入穩(wěn)定周期振蕩，進(jìn)而進(jìn)入一個小幅度隨機(jī)振動態(tài)；α分別取值0.95、1.00和1.05時，系統(tǒng)在小幅度隨機(jī)振蕩后并未直接過渡到大幅度隨機(jī)振蕩態(tài)，而是進(jìn)入周期振蕩，最后由周期振蕩直接進(jìn)入大幅度隨機(jī)振蕩態(tài)。

2.5 發(fā)電機(jī)角速度軸心動態(tài)演化過程

在機(jī)組發(fā)生故障后，很容易發(fā)生機(jī)組失速故障，機(jī)組轉(zhuǎn)速升高后會出現(xiàn)劇烈振動行為。機(jī)組軸心軌跡圖可用于分析機(jī)組失速后動力演化過程。圖8描述機(jī)組轉(zhuǎn)動角速度從0～850 rad/s變化軸心軌跡動力演化過程。

(a) 0<ω<200

(b) 200<ω<700

(c) 700<ω<850

圖8 機(jī)組轉(zhuǎn)動角速度從0～850 rad/s變化軸心軌跡動態(tài)演化過程

Fig.8 Dynamic evaluation of the center of the generator rotor when the generator speed changes in the interval (0, 850) rad/s

觀察圖8(a)，在0<ω<7 rad/s時，機(jī)組振動幅值基本為0 m；在7<ω<44 rad/s時，機(jī)組振動幅度隨轉(zhuǎn)速線性增加；在44<ω<55 rad/s時，振動幅度基本保持不變，且在ω=55 rad/s轉(zhuǎn)速發(fā)生了跳躍現(xiàn)象；在55<ω<189 rad/s時，振動幅度隨著轉(zhuǎn)速增加迅速增加，并在ω=189 rad/s達(dá)到該階段峰值；在189<ω<194 rad/s時，振動幅度迅速衰減，振動幅度減小，在194<ω<250 rad/s時，振動幅度基本不變。觀察圖8(b)，在250<ω<269 rad/s時，振動幅度基本不變；在269<ω<299 rad/s時，振動幅度略有增加，且在ω=299 rad/s時出現(xiàn)跳躍式下降；在299<ω<343 rad/s時，振動幅度先減小后趨于穩(wěn)定，且在ω=343 rad/s跳躍式上升，之后隨轉(zhuǎn)速升高振動幅度先上升后下降；當(dāng)401<ω<603 rad/s時，隨轉(zhuǎn)動角速度緩慢上升，振動幅度緩慢上升趨勢；在517<ω<530 rad/s，振動幅度先是跳躍式下降，緊接著出現(xiàn)跳躍式上升，且幅度大于之前幅度；在603<ω<700 rad/s時，振動幅度也隨轉(zhuǎn)速增加緩慢增大。觀察圖8(c)，當(dāng)700<ω<850 rad/s時，機(jī)組振動幅度剛開始緩慢增加，最后振幅迅速升高并達(dá)到最大值。

3 結(jié) 論

本文采用非線性動力學(xué)方法探究水輪發(fā)電機(jī)組結(jié)構(gòu)參數(shù)異變性對沖擊式水電站軸系振動擺度影響規(guī)律，得出如下結(jié)論：

(1) 分?jǐn)?shù)階階次值通過阻尼力響應(yīng)影響系統(tǒng)演化過程，對Hopf分岔影響不是很明顯，但對隨機(jī)振動幅值有明顯影響。

(2) 上導(dǎo)軸承、下導(dǎo)軸承和水導(dǎo)軸承剛度對系統(tǒng)穩(wěn)定性有著明顯影響；不同分?jǐn)?shù)階次下系統(tǒng)進(jìn)入隨機(jī)振動方式也不相同：①從穩(wěn)定周期振蕩直接進(jìn)入隨機(jī)振蕩；②從Hopf分岔過渡至隨機(jī)振蕩；③系統(tǒng)先進(jìn)入小幅度隨機(jī)振蕩，而后迅速過渡到大幅度隨機(jī)振蕩狀態(tài)。

(3) 通過對系統(tǒng)失速后振動幅度分析發(fā)現(xiàn)，系統(tǒng)振動幅度會隨轉(zhuǎn)速升高出現(xiàn)明顯跳躍現(xiàn)象。

(1)

(2)

(3)

(4)

聯(lián)立公式(1)～(4)，得

(5)

(6)

(7)

(8)

(9)

聯(lián)立公式(6)～(9)，可以得到

(10)

(11)

(12)

(13)

(14)

Mt-Me=(mt-me)MgB

(15)

聯(lián)立公式(5),(10)～(13)和(14)，得到

(16)

取

(17)

(18)

(19)

(20)

對公式(16)進(jìn)行化簡，得到

(21)

化簡可得發(fā)電機(jī)角速度為

(22)