李佰霖, 方子鈴, 貴樹鑫, 王云峰

(1.三峽大學(xué) 電氣與新能源學(xué)院,湖北 宜昌 443002;2.三峽大學(xué) 梯級水電站運行與控制湖北省重點實驗室,湖北 宜昌 443002;3. 福建仙游抽水蓄能有限公司,福建 莆田 351267)

水力振動是抽水蓄能電站中一種常見的危害機組安全的現(xiàn)象,水力振動產(chǎn)生的壓力遠(yuǎn)超管道正常運行的安全范圍,會造成爆管等事故,嚴(yán)重威脅抽水蓄能電站的安全運行。索麗生等[1]曾通過實例對水電站水力振動的發(fā)生進行了分析與預(yù)測;周建旭等[2-3]基于改進的水體彈性模型揭示了壓力管道的水力振動特性;段煉等[4-6]將矩陣法應(yīng)用于水電站的水力振動研究,并通過實例分析了該方法的可靠性。除了以上方法,胡光華等[7]基于特征線法對不同波動的耦合效應(yīng)進行了分析,劉冬等[8]在此基礎(chǔ)上提出了改進的特征線法;樓勇等[9]則通過CFD(computational fluid dynamics)建立三維模型的方式得到了自激振動時的壓力與流量規(guī)律;程加堂等[10]將粒子群BP(back propagation)神經(jīng)網(wǎng)絡(luò)算法應(yīng)用于水電機組振動的故障診斷。

球閥在抽水蓄能電站中起到控制流量與穩(wěn)定壓力的作用,當(dāng)球閥密封劣化漏水會導(dǎo)致球閥變成柔性閥而發(fā)生自激振動[11]。張紹春[12]對柔性閥的自激振動進行分析得到在有分叉管時,自激振動的最大壓力可以達到靜水壓力的3倍;朱淵岳等[13-14]通過特征線法(method of characteristics,MOC)和阻抗法計算的仿真結(jié)果得出了自激振動的原因,并提出了消除自激振動的措施;李紅輝等[15]建立了自激振動的數(shù)值仿真模型[16],并通過經(jīng)驗小波變換對模擬自激振動信號和實際信號的頻譜進行分析,以驗證模型的可行性;斯靜等[17]設(shè)計了兩個對比方案,得出了管道與閥門漏水量對自激振動振幅和周期的影響。除此之外,有些學(xué)者則從機理角度對自激振動進行了分析, Bouzidi等[18]研究了彈簧加載閥的自激振動機制,發(fā)現(xiàn)振動幅度與管道長度和彈簧剛度呈正相關(guān)。Awad等[19]建立了一個簡單的球閥模型來研究環(huán)形密封的自激振動機理,并分析了各種因素對模型穩(wěn)定性的影響。

基于以上研究現(xiàn)狀,過往學(xué)者對自激振動原理進行了詳細(xì)研究缺乏與現(xiàn)場實際相結(jié)合,為此本文通過特征線法建立了自激振動系統(tǒng)的綜合仿真模型,研究框架如圖1所示。期望通過物理系統(tǒng)和數(shù)字模型的完全結(jié)合來分析球閥的自激振動特性,從而對現(xiàn)場工作進行指導(dǎo),這有利于減少自激振動對水力設(shè)備的破壞,維護抽水蓄能電站的正常運行。

圖1 球閥自激振動研究框架Fig.1 Research framework for self-excited vibration of ball valve

1 建 模

某抽水蓄能電站的結(jié)構(gòu)如圖2所示。本節(jié)選擇圖2中的虛線部分進行建模。該結(jié)構(gòu)主要包括上游水庫、主管道L1以及分叉管道L2,L3、閥門1和閥門2。在該模型中,自激振動時由于閥門2處密封漏水引起的,因此閥門2被視為柔性閥門。閥門1是關(guān)閉的,因此閥門1被視為封閉端。

圖2 某抽水蓄能電站的引水結(jié)構(gòu)Fig.2 The water diversion structure of a pumped storage power station

1.1 柔性閥門建模

抽水蓄能電站球閥根據(jù)漏水特性可分為剛性閥門和柔性閥門,通過分析活動環(huán)的受力情況可以判斷它能否正常與球閥表面貼合來確定閥門的漏水特性?；顒迎h(huán)受到水壓力的合力如下

F=P2S2+P3S3-P1S1

(1)

式中:F為活動環(huán)受到水壓力的合力;P1為輸入腔的水壓;P2為閥體內(nèi)部水壓;P3為輸出腔的水壓;S1為輸入腔的水壓對活動環(huán)的壓力面積;S2為閥體內(nèi)的水壓對活動環(huán)的壓力面積;S3為輸入腔的水壓對活動環(huán)的壓力面積;S1=S2+S3。

工作密封的結(jié)構(gòu)如圖3所示。剛性閥門漏水特性分析:正常情況下密封不會發(fā)生漏水。由于輸入腔的水是從球閥的上游側(cè)管道引入的,因此輸入腔的水壓等于閥體內(nèi)部水壓。輸出腔內(nèi)的水是常壓水,其壓力遠(yuǎn)小于輸入腔的水壓,并且輸入腔的水壓對活動環(huán)的壓力面積大于閥體內(nèi)部水壓對活動環(huán)的壓力面積?；顒迎h(huán)受到的水壓力將其向左推動,使其與球閥表面貼合。

圖3 工作密封結(jié)構(gòu)Fig.3 Structure of working seal

柔性閥門漏水特性分析:當(dāng)輸入腔與輸出腔之間或閥體與輸出腔之間的密封泄漏時,輸出腔的水壓也等于球閥的上游側(cè)管道壓力,活動環(huán)受到的水壓力處于平衡狀態(tài)?；顒迎h(huán)不會向閥體移動,也無法完成工作密封的投入。如果此時出現(xiàn)微小的壓力擾動,使閥門漏水面積減小,球閥上游側(cè)的管道壓力會隨之增大。閥體、輸入腔和出輸出腔中的水壓也會隨著上游水壓的增加而增加。由于輸出腔的水壓變化是由密封漏水引起的,因此它的水壓變化是一個緩慢的過程,具有一定的滯后性。因此,輸出腔的水壓小于輸入腔的水壓。此時,活動環(huán)受到水壓力向左,閥門的漏水面積進一步減小。球閥上游水壓的持續(xù)增加導(dǎo)致活動環(huán)繼續(xù)受到向左的力的影響。這形成了一個正反饋過程。

半個周期后,從上游水庫反射的水錘波導(dǎo)致閥門上游側(cè)壓力減小,閥體、輸入腔和輸出腔中的水壓也隨之減小。由于輸出腔的水壓變化具有一定的滯后性,此時輸出腔的水壓大于輸入腔的水壓?；顒迎h(huán)受到的水壓力向右,導(dǎo)致閥門漏水面積增大,這也是一個正反饋過程。水錘波在管道中來回傳播使這一過程不斷重復(fù),這會導(dǎo)致壓力振幅不斷增大最終發(fā)生自激振動。

柔性閥門的漏水面積變化如圖4所示。柔性閥的漏水面積與水壓之間的關(guān)系為

圖4 柔性閥門漏水面積Fig.4 Leakage area of flexible valve

(2)

式中:H0為初始壓力;A1為靜壓下的漏水面積,H1為密封處水壓。

柔性閥門與剛性閥門的漏水特性如圖5所示。剛性閥門的漏水量和閥門前后的壓差的關(guān)系式如下

圖5 剛性閥門和柔性閥門的漏水特性Fig.5 Leakage characteristics of rigid valve and flexible valve

(3)

式中:k為常數(shù);Q為閥門的漏水量;H為閥門前后的壓差。

柔性閥門的漏水量和閥門前后的壓差的關(guān)系式如下

(4)

式中:D1為閥門直徑;x0為閥門前后壓差為0時的密封間隙;Af為閥門密封蓋的面積;K為閥門密封的彈性系數(shù);ρ為水的密度。

由于密封的壓力水是從球閥上游側(cè)的管道引入的,因此H1=H。因此,柔性閥門漏水量可以表示為

(5)

1.2 特征線法

以水錘基本方程為基礎(chǔ),通過特征線法將水錘基本方程轉(zhuǎn)化為有限差分方程進行計算。特征線法的矩形計算網(wǎng)絡(luò)如圖6所示。將管道劃分為每段長為Δx,將時間被劃分為每個時間段Δt=Δx/a。C+為正特征線,C-為負(fù)特征線,i為管道節(jié)點。計算從t=0開始。除邊界點外的其他節(jié)點的參數(shù)可以通過特征線方程計算獲得,邊界點可以通過特征線方程結(jié)合邊界條件計算獲得。

圖6 特征線法的矩形計算網(wǎng)絡(luò)Fig.6 Rectangular calculation network of MOC

特征線方程

C+∶HPi=CP-BQPi

(6)

C-∶HPi=CM+BQPi

(7)

CP=Hi-1+BQi-1-RQi-1|Qi-1|

(8)

CM=Hi+1-BQi+1+RQi+1|Qi+1|

(9)

上游邊界條件

HP=Hu

(10)

式中,Hu為上游水庫水頭。

下游邊界條件

(11)

2 參數(shù)優(yōu)化

由于閥門漏水量無法準(zhǔn)確測量,因此很難獲得柔性閥門的特征參數(shù)k1和k2,無法精確模擬該抽水蓄能電站球閥自激振動過程。為了使模擬的自激振動過程符合實際情況,需要對柔性閥門的特性參數(shù)進行優(yōu)化。

2.1 適應(yīng)度函數(shù)

為了找到k1和k2的最優(yōu)解本文將仿真結(jié)果與監(jiān)測數(shù)據(jù)相結(jié)合得到了反映特征參數(shù)優(yōu)良程度的適應(yīng)度函數(shù)。由于監(jiān)測數(shù)據(jù)的采樣率較低,很難充分反映壓力隨時間的變化過程。如果通過比較每個時間點的數(shù)據(jù)來評估仿真結(jié)果和監(jiān)測數(shù)據(jù)之間的一致性就會產(chǎn)生很大的誤差。為了使比較結(jié)果更加準(zhǔn)確,本文對監(jiān)測數(shù)據(jù)與仿真結(jié)果的上下包絡(luò)線進行了比較,以確定實際和模擬自激振動發(fā)展過程的一致性。形成監(jiān)測數(shù)據(jù)包絡(luò)線方程的點如下

(12)

式中:HP(t) 為壓力數(shù)據(jù);Hm為上包絡(luò)線的點;Hn為下包絡(luò)線的點;t為自激振動時間。

為了使曲線更平滑,對曲線進行擬合。由于曲線的函數(shù)類型尚未確定,因此可以通過泰勒級數(shù)展開將曲線的表達式轉(zhuǎn)換為多項式函數(shù),最小二乘法通過計算最小誤差的平方和來找到擬合多項式函數(shù)的最佳函數(shù)匹配。擬合后上下包絡(luò)線的函數(shù)表達式為

(13)

式中:Hm1(t) 為監(jiān)測數(shù)據(jù)優(yōu)化后的上包絡(luò)線;Hn2(t) 為監(jiān)測數(shù)據(jù)優(yōu)化后的下包絡(luò)線。

仿真結(jié)果的包絡(luò)線是通過相同的方法求得,將監(jiān)測數(shù)據(jù)與仿真結(jié)果的包絡(luò)線進行了比較,并采用時間乘絕對誤差積分準(zhǔn)則作為性能指標(biāo),得到反映仿真結(jié)果和監(jiān)測數(shù)據(jù)一致性程度的函數(shù)如下

(14)

式中:Hm2(t) 為仿真結(jié)果的上包絡(luò)線;Hn2(t) 為仿真結(jié)果的下包絡(luò)線。

由于不同的柔性閥門特性參數(shù)會導(dǎo)致不同的閥門漏水量,需要限定漏水量的范圍,以排除不符合要求的特征參數(shù)。從式(2)中可以得出,初始漏水面積為A1,漏水面積在0和2A1之間來回變化。因此,0

(15)

式中,Q0為閥門的初始漏水量。

對于不滿足式(15)的柔性閥特征參數(shù),在性能指標(biāo)中加入懲罰函數(shù)。懲罰函數(shù)如下

(16)

獲得適應(yīng)度函數(shù)如下

y(j)=F1(t)+ε(t)

(17)

2.2 粒子群優(yōu)化

采用粒子群算法進行優(yōu)化,粒子的速度和位置更新公式如下

(18)

式中:xj為粒子位置;vj為粒子的更新速度;w為慣性因子;c1和c2為加速度常數(shù);pbestj為局部最優(yōu)位置;gbestj為全局最優(yōu)位置。

將基于MOC的自激振動模型代入PSO進行參數(shù)優(yōu)化的流程如圖7所示。主要步驟如下:

圖7 參數(shù)優(yōu)化流程圖Fig.7 Parameters optimization flow chart

步驟1通過設(shè)置適當(dāng)?shù)南到y(tǒng)參數(shù)并限定粒子的范圍隨機生成粒子。

步驟2將每組粒子代入MOC程序得到仿真結(jié)果,通過現(xiàn)場狀態(tài)監(jiān)測裝置獲得監(jiān)測數(shù)據(jù),通過適應(yīng)度函數(shù)計算每組粒子的適應(yīng)度。

步驟3將初始適應(yīng)度作為每個粒子的局部最優(yōu)解,并找到全局最優(yōu)解。

步驟4通過速度和位置更新公式,更新粒子的速度和位置。超出范圍的粒子將根據(jù)邊界條件進行處理。

步驟5計算每個新粒子的適應(yīng)度與局部最優(yōu)解比較,找到新的局部最優(yōu)解,并從所有新的局部最優(yōu)解中找到全局最優(yōu)解。

步驟6重復(fù)步驟4和步驟5,直到找到最優(yōu)解。

抽水蓄能電站初始水頭為516 m,具體管道參數(shù)如表1所示。監(jiān)測數(shù)據(jù)的擬合曲線函數(shù)為4階,參數(shù)為a1=52.38,a2=-3.504×10-3,a3=6.581×10-5,a4=-2.323×10-7,a5=3.776×10-10,b1=50.9,b2=4.508×10-3,b3=-6.935×10-5,b4=2.321×10-9,b5=-3.7×10-10。粒子群優(yōu)化算法設(shè)置50個粒子并迭代30次;w=0.8;c1=1.5;c2=1.5。

表1 管道參數(shù)Tab.1 Pipeline Parameters

代入數(shù)據(jù)計算后,適應(yīng)度函數(shù)曲線如圖8所示。在20次迭代之后,曲線趨于平緩,結(jié)果達到穩(wěn)定。求得柔性閥的特征參數(shù)為k1=0.039 8,k2=3.19×10-5。

圖8 適應(yīng)度函數(shù)曲線Fig.8 Fitness function curve

優(yōu)化前后仿真結(jié)果與監(jiān)測數(shù)據(jù)的對比如圖9所示?？梢钥吹絻?yōu)化前的壓力數(shù)據(jù)與實際情況存在較大差距;優(yōu)化后仿真結(jié)果在400 s前與實際自激振動過程存在較小誤差,在400 s后的發(fā)展過程基本相同。通過粒子群算法優(yōu)化的柔性閥門特征參數(shù)符合實際情況。

圖9 優(yōu)化前后球閥上游側(cè)壓力對比Fig.9 Comparison of upstream side pressure of ball valve before and after optimization

3 分 析

在獲得與實際情況一致的模型后,本章通過控制變量法分析了初始擾動,管道長度,管道直徑,雙柔性閥門系統(tǒng)對自激振動的影響,確定自激振動特性為抽水蓄能電站的優(yōu)化設(shè)計提供了依據(jù)。

3.1 初始壓力擾動對自激振動的影響

僅改變閥門1處受到的初始壓力擾動,不改變其他參數(shù),設(shè)計對比方案1,其球閥上游側(cè)壓力如圖10所示,統(tǒng)計不同初始擾動其振幅時間如表2所示。

表2 不同壓力擾動的振幅時間Tab.2 Amplitude time of different pressure disturbances

圖10 方案1對比結(jié)果Fig.10 Comparison results of option 1

如表2所示:從1 m水頭到2.5 m水頭在兩種不同初始擾動情況下所需時間相同,均為105 s;從2.5 m水頭到400 m水頭在3種不同初始擾動情況下所需時間相同,均為600 s。在不同的初始壓力擾動下,自激振動從一個振幅到另一個振幅的增長速度是相同的。對比3種不同初始擾動的發(fā)展過程,可以發(fā)現(xiàn)自激振動在初始振幅較小的階段緩慢增加,在后期振幅較大時增加迅速,因此初始擾動越大,自激振動增長速度也就越快。

3.2 管道長度對自激振動的影響

為探究自激振動與不同管道長度的關(guān)系,設(shè)計對比方案2,僅改變主管道長度與分叉管道長度,其他條件不變,共設(shè)計了9種管道長度組合,其球閥上游側(cè)壓力如圖11所示,統(tǒng)計不同管道長度組合的振幅與頻率如表3所示。

表3 不同管道長度組合的振幅與頻率Tab.3 Amplitude and frequency of different pipeline length combinations

圖11 方案2對比結(jié)果Fig.11 Comparison results of option 2

在前6種管道長度組合中僅改變主管道長度,球閥上游側(cè)壓力信號的主頻率仍為1.67 Hz,而后3種管道長度組合中僅改變分叉管道長度,自激振動壓力信號的主頻率隨著分叉管道長度的增加而減小,這是由于水錘波在異徑點處發(fā)生了反射,因此自激振動頻率主要與分叉管道長度有關(guān)。

僅增大主管道長度,自激振動幅值并沒有隨之增加,由幅值結(jié)果可以判斷自激振動與不同管道長度的比例有關(guān)。在復(fù)雜管道的情況下,假設(shè)每條管道的水錘波速相等,當(dāng)主管道長度是分叉管道長度的偶數(shù)倍時的自激振動的增長速度遠(yuǎn)大于當(dāng)主管道長度不是分叉管道長度偶數(shù)倍時的增長速度。這是由于水錘波在管道中來回流動時,會在不同管道的異徑點處發(fā)生透射與反射。當(dāng)主管道長度是分叉管道長度的偶數(shù)倍時,不同頻率的水錘波之間的疊加效應(yīng)更明顯;當(dāng)主管道長度不是分叉管道長度的偶數(shù)倍時,不同頻率的水錘波之間的抵消效應(yīng)更明顯。僅增大分叉管道長度,自激振動的頻率減小,這就導(dǎo)致了自激振動的增長速度也隨之減小。

3.3 管道直徑對自激振動的影響

為探究管道直徑對自激振動的影響,設(shè)計對比方案3,僅改變主管道直徑與分叉管道直徑,其他條件不變,其球閥上游側(cè)壓力如圖12所示,不同管道直徑的自激振動幅值如表4所示。

表4 不同管道直徑的自激振動幅值Tab.4 Self excited vibration amplitudes of different pipeline diameters

圖12 方案3對比結(jié)果Fig.12 Comparison results of option 3

當(dāng)僅改變主管道直徑時,與初始狀態(tài)相比自激振動幅值基本未發(fā)生改變,因此主管道直徑對自激振動的增長速度影響不大;當(dāng)僅改變主管道直徑時,分叉管道與閥門直接相連,分叉管道直徑變化會導(dǎo)致閥門附近的水錘波變化從而直接影響自激振動的發(fā)展過程,因此從仿真結(jié)果可以看到自激振動的增長速度隨著分叉管道直徑的增大而減小;同時減小或增大主管道和分叉管道直徑的自激振動結(jié)果與以上分析結(jié)果相同。

3.4 雙柔性閥門系統(tǒng)對自激振動的影響

當(dāng)兩處閥門均發(fā)生密封漏水變?yōu)槿嵝蚤y門時,設(shè)計對比方案4,將閥門1由剛性閥門變?yōu)槿嵝蚤y門,其他管道等參數(shù)不變,設(shè)計了3種情況下的自激振動仿真,閥門2上游側(cè)壓力如圖13所示,3種情況下的自激振動幅值如表5所示。

表5 雙柔性閥門系統(tǒng)不同情況下的自激振動幅值Tab.5 Amplitude of self-excited vibration of double flexible valve system under different conditions

圖13 方案4對比結(jié)果Fig.13 Comparison results of option 4

當(dāng)系統(tǒng)中的兩個閥門均為柔性閥門時,其他條件均不發(fā)生改變,其仿真結(jié)果如情況1所示,300 s時的自激振動幅值明顯增大,這是由于有兩個閥門同時發(fā)生振動加速了自激振動的發(fā)展。在柔性閥門1處增加1個與柔性閥門2處相同的擾動,其仿真結(jié)果如情況2所示,300 s時的自激振動幅值僅有14 m,這是由于兩處閥門上的擾動產(chǎn)生的水錘波在管道中傳播時相向運動發(fā)生了抵消導(dǎo)致自激振動增長緩慢。在柔性閥門1處增加一個與柔性閥門2處相反的擾動,其仿真結(jié)果如情況3所示,300 s時的自激振動幅值與情況1相比進一步增大,這是由于兩處閥門擾動產(chǎn)生的水錘波在管道中傳播時同向運動發(fā)生了疊加導(dǎo)致管道中的擾動進一步增大,使得自激振動幅值也快速增大。當(dāng)柔性閥門擾動發(fā)生的時間不同時,其仿真結(jié)果如情況4、情況5所示,由于擾動有時差,水錘波在傳播時出現(xiàn)抵消或者疊加,最終導(dǎo)致自激振動結(jié)果不同。

4 結(jié) 論

為研究抽水蓄能電站球閥自激振動的發(fā)展過程本文提出了一種球閥自激振動的研究框架,通過物理系統(tǒng)的監(jiān)測數(shù)據(jù)來優(yōu)化數(shù)字模型的建立,并分析數(shù)字模型的自激振動特性,以指導(dǎo)物理系統(tǒng)的工作。在優(yōu)化過程中,本文基于MOC和PSO對模型參數(shù)進行優(yōu)化。通過將優(yōu)化后仿真結(jié)果與監(jiān)測數(shù)據(jù)進行比較,驗證了優(yōu)化參數(shù)的可靠性。最后通過兩個對比方案發(fā)現(xiàn)了自激振動特性,結(jié)論如下:

(1)初始擾動越大,自激振動增長速度也就越快。

(2)由于水錘波在異徑點處的反射,分叉管道長度會影響自激振動頻率,同時又會影響自激振動的發(fā)展過程。

(3)不同管道的水錘波速相同時,由于水錘波在異徑點處的透射與反射,當(dāng)主管道長度是分叉管道長度的偶數(shù)倍時自激振動的增長速度遠(yuǎn)大于當(dāng)主管道長度不是分叉管道長度偶數(shù)倍時的增長速度。

(4)分叉管道直徑會影響閥門附近的水錘波,從而導(dǎo)致自激振動的增長速度隨著分叉管道直徑變化而變化。

(5)當(dāng)系統(tǒng)中的兩個閥門均為柔性閥門時,兩處閥門處的擾動產(chǎn)生的水錘波會在管道中發(fā)生抵消或疊加進而導(dǎo)致自激振動增長速度的增大或減小。