謝卓健,閆懂林,肖 漢,麥先春,王峻峰,陳啟卷
(1.四川大唐國際甘孜水電開發(fā)有限公司,四川康定626001;2.武漢大學動力與機械學院,武漢430072)
在過去的二十年里,我國水電的開發(fā)和利用得到了快速的發(fā)展,在西南地區(qū),各種大流量、高水頭、包含復雜水力系統(tǒng)的水電站不斷興起。到2019年底,我國的水電總裝機容量已經(jīng)達到3.56 億kW[1]。對于水電站來說,由于電網(wǎng)故障等原因,運行中隨時可能會經(jīng)歷機組甩負荷過程,而一些水電站的重大事故,例如飛逸、壓力鋼管爆裂等,往往就是由于這種較為極端的工況引起[2]。如何有效地保證甩負荷過渡過程中機組的安全可靠,是實際水電工程中所面臨的重要問題。
到目前為止,針對水電機組甩負荷過渡過程進行了大量研究。部分學者對實際電站發(fā)生的甩負荷事故進行了分析和討論。黃菊珠[3]分析了山美水電站因110 kV 線路開關跳閘,造成3臺機組滿負荷高頻切機,引發(fā)機組過速停機的事故,并提出了相關的預防措施。曲東[4]分析了平班水電站在甩負荷試驗時調(diào)速器調(diào)節(jié)失靈,導葉不能回關導致機組過速停機的事故,提出了相應的處理方法,并對調(diào)速器重要參數(shù)進行了優(yōu)化整定。張偉[5]針對西溪電站甩負荷試驗時出現(xiàn)的機組過速停機及導葉快關時間不穩(wěn)定等問題進行了分析,并提出了故障處理方法。除此之外,也有大量學者對水力發(fā)電機組甩負荷過渡過程進行了數(shù)值模擬研究。付亮[6]針對有調(diào)壓井長引水隧洞水電站,利用數(shù)值模擬的方法研究了調(diào)壓井水位波動與流道壓力的變化過程。王威[7]等以伊朗卡侖水電站為例,分析了具有上游調(diào)壓室、一管兩機布置的水電站機組甩負荷過渡過程,并討論了導葉關閉規(guī)律對轉(zhuǎn)速上升和水擊壓力的影響,提出了合理的導葉關閉規(guī)律。葉文波等[8]針對南方某大型水電站,進行了機組甩負荷大波動過渡過程的計算,分析了引水系統(tǒng)最高壓力升高和機組轉(zhuǎn)速的最大升高特性。陳勝等[9]研究了水電站中機組相繼甩負荷情況下調(diào)壓室涌浪的疊加特性。張洋[10]研究了設置氣墊式調(diào)壓室的超長引水隧洞水電站的甩負荷過渡過程,提出了氣墊式調(diào)壓室對系統(tǒng)水擊壓力的影響效應。在這些分析和討論中,都是在確定性框架下研究水電站機組的甩負荷過渡過程,而并沒有考慮到參數(shù)不確定性對甩負荷過渡過程的影響。實際上,水電站過渡過程分析往往參數(shù)眾多且包含復雜的非線性動態(tài)過程,參數(shù)尺度多樣化及制造誤差等多種不確定性因素均會導致計算結(jié)果的不確定性。所以,在不確定框架下研究參數(shù)對甩負荷過渡過程的影響效應就顯得更有必要。
參數(shù)敏感性分析是量化系統(tǒng)參數(shù)在設計或設定的可能范圍內(nèi)變化時對系統(tǒng)輸出響應影響的一種方法[11]。通過敏感性分析,可以獲得對系統(tǒng)輸出影響較為顯著的輸入?yún)?shù),無論在設計、運行還是維護過程中,只需要重點關注這些參數(shù)而忽略一些不敏感的參數(shù),就可以大大簡化工作量。通常的參數(shù)敏感性分析包括局部敏感性分析和全局敏感性分析,局部敏感性分析只能反映單一參數(shù)變化對系統(tǒng)響應的影響,而全局敏感性分析則可以有效地反映多個參數(shù)同時變化對系統(tǒng)響應的影響[12]。在全局敏感性分析方法中,主要包括:回歸分析法[13]、篩選法[14],Sobol 指數(shù)法[15]和擴展傅里葉幅度檢驗法(extended Fouri?er amplitude sensitivity test,EFAST)[16]等。其中EFAST 方法可以獲得輸入?yún)?shù)之間的交互效應,已經(jīng)在很多領域得到了廣泛應用[17,18]。因此,本研究中將利用EFAST 全局敏感性分析方法來研究水力發(fā)電機組在甩負荷過渡過程中輸入?yún)?shù)對轉(zhuǎn)速、水擊壓力和調(diào)壓室涌浪的敏感性特性。
相比于之前的研究,本文的主要創(chuàng)新在于考慮參數(shù)的不確定性,在不確定框架下建立了水力發(fā)電機組的甩負荷過渡過程的模型,并以此模型為基礎進行了數(shù)值研究,量化分析了輸入?yún)?shù)對機組甩負荷過程中機組轉(zhuǎn)速、水擊壓力和調(diào)壓室涌浪的敏感性特性。相關的理論方法和研究結(jié)論可以應用于水電站的設計及維護工作。
水電站的引水系統(tǒng)在空間分布上主要包括上游水庫,引水隧洞,上游引水調(diào)壓室,壓力鋼管,蝸殼,尾水管及尾水明渠段。本節(jié)將首先利用特征線方法建立該電站引水系統(tǒng)的數(shù)學模型,同時結(jié)合水輪機特性曲線的插值模型,調(diào)速器數(shù)學模型,發(fā)電機及負載數(shù)學模型,建立完整的水電機組甩負荷過渡過程的數(shù)學模型。
1.1.1 有壓過流管道特征線方程
利用特征線法將有壓輸水系統(tǒng)非恒定流的運動方程和連續(xù)性方程轉(zhuǎn)化為在特征線dx/dt=V±a上的常微分方程組[19],可以得到:
式中:H表示測壓管水頭;a表示水擊波速;f為達西-維斯巴哈摩擦系數(shù);D表示管道直徑;A為流道截面面積;Q為流道中的流量。在這里,V和H是與時間和空間相關的函數(shù),其中t為時間,x表示管道長度。
在利用特征線法離散化求解微分方程時,需要對有壓管道按照等距離進行劃分,得到長度為Δx的N段管,并用序號i進行標記,定義起始端為i=1,則末端為i=N+1。設定仿真步長為ΔT=Δx/a。利用特征線法進行計算機數(shù)值求解的原理如圖1所示,橫坐標表示管道的空間分布,縱坐標表示采樣時刻分布。從中可以看出,根據(jù)上一時刻的節(jié)點i-1 和節(jié)點i+1 的值,可以計算出下一時刻節(jié)點i的水頭、流量值。
圖1 特征線法原理示意Fig.1 Schematic diagram of characteristics method
對公式(1)和(2)的特征線方程積分并離散化,可以得到t+ΔT時刻節(jié)點i處的水頭、流量關系方程式為:
式中:B=a/Ag;R=fΔx/(2gDA2);CPi,t+ΔT=Hi-1,t+BQi-1,t;BPi,t+ΔT=B+R|Qi-1,t|;CMi,t+ΔT=Hi+1,t-BQi+1,t;BMi,t+ΔT=B+R|Qi+1,t|。
1.1.2 水輪機邊界
水輪機的瞬態(tài)工作水頭定義為水輪機蝸殼進口和尾水管出口的水頭差,如圖2所示,當水輪機上、下游節(jié)點為1 和2 時,則水輪機工作水頭為HT=H1-H2,瞬態(tài)的過機流量為QT=Q1=Q2。
圖2 水輪機簡化結(jié)構(gòu)圖Fig.2 Structure chart of hydro-turbine
根據(jù)水輪機上下游進出口節(jié)點的水頭平衡方程,可以得到:
水庫邊界,上游調(diào)壓室邊界,蝸殼及尾水管當量管的計算可以參考文獻[19]。
依據(jù)文獻[20]的方法,對原始的水輪機綜合特性曲線進行擴展,可以得到Q11=f1(n11,y)和M11=f2(n11,y)的關系,用于甩負荷過渡過程的計算。在這兩個關系式中,n11表示單位轉(zhuǎn)速,M11表示單位力矩,y表示導葉開度。以當前研究的水電站為例,得到擴展后的Q11-n11-y和M11-n11-y關系曲線如圖3所示。
圖3 單位流量和單位力矩圖Fig.3 Diagrams of unit discharge and unit moment
水力發(fā)電機組的力矩平衡方程[19]可以表示為:
式中:J為機組轉(zhuǎn)動部分的轉(zhuǎn)動慣量,J=GD2/4,t·m2;GD2為飛輪力矩,t·m2;ω表示機組旋轉(zhuǎn)角速度,rad/s;M和Mg分別是水輪機的主動力矩和阻力矩,kN·m。
假定水輪機力矩在ΔT時刻里線性變化,則可以利用梯形積分公式對上式積分。同時,將旋轉(zhuǎn)角速度(ω)轉(zhuǎn)換為轉(zhuǎn)速(n),可以得到:
至此,水力發(fā)電機組甩負荷過渡過程的確定性模型被建立,在這個確定性模型中將確定性的參數(shù)用不確定的分布或者范圍進行替代,即可以獲得不確定性的模型,以此為基礎,相關的參數(shù)敏感性可以利用方便的非侵入式EFAST方法進行討論。
本文采用EFAST 方法來對水電機組甩負荷過渡過程的參數(shù)敏感性進行分析。EFAST 法的原理是在輸入?yún)?shù)的多維空間中挑選合適的搜索曲線,并給所有輸入?yún)?shù)分配非線性相關的整數(shù)頻率值,從而引入搜索函數(shù),使模型作為獨立參數(shù)的周期函數(shù),則可以把多維積分問題轉(zhuǎn)換為一維積分問題。然后通過把目標函數(shù)變?yōu)楦道锶~級數(shù),計算出各頻率的傅里葉頻譜曲線,由頻譜曲線獲得由輸入?yún)?shù)xi引起的系統(tǒng)輸入方差和總方差,其比值則被定義為該參數(shù)的敏感度指標。文獻[16-18]中給出了具體的算法推演流程,這里僅給出敏感度計算的公式。
總方差可以寫成各個子方差和的形式:
式中:V為模型總方差;Vi表示單一輸入變化引起的模型方差;Vij表示通過xi和xj交互作用引起的方差;Vij…k表示xi,xj,…,xk交互作用引起的耦合方差。
主敏感性(一階敏感性)指標表示單個參數(shù)獨立作用時對系統(tǒng)輸出響應的貢獻,可以表示為:
總敏感性指標為:
式中:V-i表示不包括參數(shù)xi的所有參數(shù)方差和。
以國內(nèi)某水電站為例,進行機組甩100%負荷時過渡過程計算的參數(shù)敏感性分析,涉及的主要參數(shù)如表1所示。敏感性分析的輸入?yún)?shù)均值選取設計值,方差分別取變異系數(shù)(=標準差/均值)為0.01和0.05兩種情況式的值,代表不同的不確定度。選擇Monte Carlo 方法進行參數(shù)的隨機采樣,采樣樣本為5 000,滿足EFAST 法對樣本個數(shù)的要求。另外,由于其他水工建筑物的影響,該電站不同高程下的調(diào)壓室面積會有所有差別,在高程Hs0<1 458.7 m 時,面積As2為801.975 m2;當1 458.7 m
表1 仿真參數(shù)Tab.1 Simulation parameters
在機組的甩負荷過程中,最高轉(zhuǎn)速、最大水擊壓力、調(diào)壓室最大涌浪和最小涌浪需要被重點關注,也是評估機組甩負荷過程安全性的重要指標,所以,在本節(jié)中,參數(shù)相對于這幾個輸出變量的敏感性將被重點分析。
機組甩負荷過程中最高轉(zhuǎn)速的全局敏感性分析結(jié)果如圖4所示。
圖4 最高轉(zhuǎn)速的全局敏感性Fig.4 Global sensitivity of maximum speed
以甩負荷過程中的機組最高轉(zhuǎn)速作為輸出時,無論在變異系數(shù)為0.01 還是0.05 的情況下,機組旋轉(zhuǎn)部分的飛輪力矩都是對其最為敏感的參數(shù),其他參數(shù)的影響微乎其微。另外,對比飛輪力矩的主效應和總效應指標,可以發(fā)現(xiàn)在不同的不確定度下其差值都很小,即飛輪力矩對甩負荷過程中轉(zhuǎn)速升高的影響也主要是通過自身的作用實現(xiàn),與其他參數(shù)的交互作用對轉(zhuǎn)速升高的影響可以忽略。最后,對比在變異系數(shù)為0.01 和0.05 情況下對最高轉(zhuǎn)速的參數(shù)敏感性特性,發(fā)現(xiàn)變異系數(shù)為0.05 時上游引水隧洞面積對最高轉(zhuǎn)速的敏感性要高于變異系數(shù)為0.01時的敏感性。其主要原因是,在系統(tǒng)變異系數(shù)從0.01 同尺度放大到0.05 時,引水隧洞面積對最高轉(zhuǎn)速影響的放大尺度要高于其他參數(shù)對最高轉(zhuǎn)速影響的放大尺度,從而導致不同不確定度下會出現(xiàn)這種差異。當然從整體的角度來說,飛輪力矩的絕對主導地位仍然不會受到威脅??傮w來說,飛輪力矩對于甩負荷過程中的機組轉(zhuǎn)速升高的影響最為顯著。
除了轉(zhuǎn)速的升高,水擊壓力也是甩負荷過程中影響機組安全的關鍵指標。如圖5所示,給出了機組甩負荷過程中最高水擊壓力的全局敏感性分析結(jié)果。
圖5 最高水擊壓力的全局敏感性Fig.5 Global sensitivity of maximum water hammer pressure
當變異系數(shù)為0.01 時,輸入?yún)?shù)對水擊壓力敏感性的主效應指標排序為:A2>a1>As4>r1>As3>a2>a4>a3>A5>f1>f2>As1>A4>A3>GD2>f5>f3>A1>a5>As2>f4>r2;輸入?yún)?shù)對水擊壓力敏感性的總效應指標排序為:A2>a3>As4>As3>A5>a1>f1>a2>f5>r1>a4>f2>a5>GD2r2>As1>A3>A1>As2>f4>A4>f3。無論是主效應還是總效應,壓力鋼管面積對水擊壓力的敏感性最為顯著,且從中可以看出A2獨立的影響占46.6%,與其他參數(shù)的交互效應占53.4%,也就是說A2對水擊壓力的獨立作用和與其他參數(shù)的交互作用效果基本持平。此外,蝸殼水擊波速a3,調(diào)壓室面積As4和As3的總效應指標都超過了0.4,它們對最高水擊壓力的影響也僅次于A2。與A2的作用方式不同的是,a3、As4和As3都是主要通過與其他參數(shù)的交互作用來實現(xiàn)對水擊壓力的影響。
當變異系數(shù)為0.05 時,輸入?yún)?shù)對水擊壓力敏感性的主效應指標排序為:a3>A2>r1>A4>As3>A1>A3>a4>f1>GD2>A5>As1>As4>As2>a5>a2>a1>f2>f3>f5>f4>r2;輸入?yún)?shù)對水擊壓力敏感性的總效應指標排序為:a3>A2>As3>r1>A4>f5>A5>f1>As4>GD2>f2>a1>f4>As2>a5>a4>a2>A3>r2>As1>A1>f3。在這種情況下,蝸殼水擊波速a3會超越壓力鋼管面積A2成為最敏感的參數(shù)。同時,a3獨立作用對水擊壓力影響的占比也從變異系數(shù)為0.01 時的11.51%增加到49.54%,這也就意味著隨著系統(tǒng)不確定度的增加,a3獨立的影響會被放大。其主要原因是系統(tǒng)不確定度的增加會使得參數(shù)的變化范圍也相應地擴大,在這個過程中,a3在偏離標準值越遠時其單參數(shù)項對水擊壓力影響的變化程度更加顯著,從而使得其獨立影響的占比在不確定度增加后提升而交互效應的占比下降。此外,相比于變異系數(shù)為0.01 時的敏感性指標,As4和As3的影響在變異系數(shù)為0.05時會被削弱。值得注意的是,在變異系數(shù)為0.05時,調(diào)壓室流入系數(shù)r1和蝸殼截面積A4對最大水擊壓力的影響會顯著提升。
圖6描述了機組甩負荷過程中調(diào)壓室最高涌浪的全局敏感性分析結(jié)果。在變異系數(shù)為0.01 時,調(diào)壓室截面積(As2、As3和As4),調(diào)壓室流入系數(shù)(r1)以及引水隧洞的面積(A1)對最高涌浪較為敏感,它們之間的總敏感性排序為:As3>A1>As4>r1>As2>其他。同時,可以注意到參數(shù)As3,A1,As4和r1主要是通過它們獨立的作用直接影響最高涌浪高度,而參數(shù)As2則是通過與其他參數(shù)的交互過程形成對調(diào)壓室最高涌浪高度的影響。剩余的參數(shù)對調(diào)壓室最高涌浪的影響則不明顯。在變異系數(shù)為0.05 時,輸入?yún)?shù)對調(diào)壓室最高涌浪的總敏感性排序仍然為:As3>A1>As4>As2>r1>其他。另外,與變異系數(shù)為0.01 時的參數(shù)敏感性指標進行對比,可以發(fā)現(xiàn)它們之間的差異并不明顯。也就是說,參數(shù)不確定性的變化對調(diào)壓室最高涌浪的影響比較微弱,可以忽略不計。
圖6 調(diào)壓室最高涌浪的全局敏感性Fig.6 Global sensitivity of highest surge
圖7給出了機組甩負荷過程中調(diào)壓室最低涌浪的全局敏感性分析結(jié)果。無論變異系數(shù)為0.01 還是0.05,調(diào)壓室面積As3和壓力隧洞面積A1是對調(diào)壓室最低涌浪影響最為顯著的兩個參數(shù),而其他參數(shù)的影響則相對較小。與調(diào)壓室最高涌浪的敏感性結(jié)果對比可以發(fā)現(xiàn),參數(shù)A1和As4的影響會明顯減小,但As3對最高涌浪和最低涌浪的影響無論在哪個不確定度下都十分顯著。
圖7 調(diào)壓室最低涌浪的全局敏感性Fig.7 Global sensitivity of minimum surge
本文以國內(nèi)某水電站為例,建立了水力發(fā)電機組甩負荷過渡過程的數(shù)學模型,并考慮到參數(shù)的不確定性,引入EFAST 方法對輸入?yún)?shù)在不同不確定度下進行了敏感性分析,可以得到以下結(jié)論。
(1)對于機組甩負荷過程中轉(zhuǎn)速的最大上升值,機組的飛輪力矩對其影響占據(jù)了絕對的主導地位,并且這種影響效應是飛輪力矩獨立作用的效果而不是與其他參數(shù)交互影響的結(jié)果,且參數(shù)的不確定度對參數(shù)敏感性的影響不明顯。
(2)在較小的不確定度下,壓力鋼管面積A2,蝸殼波速a3,調(diào)壓室面積As4和As3,以及尾水管面積A5是對水擊壓力影響最為顯著的5 個參數(shù)。而當參數(shù)的不確定度增加后,蝸殼波速a3對水擊壓力的影響會提升并成為影響最顯著的參數(shù),而調(diào)壓室面積As4和As3的影響在不確定度增大后會被削弱。
(3)對于調(diào)壓室涌浪,調(diào)壓室截面積As3及引水隧洞的面積A1對調(diào)壓室最高最低涌浪都較為敏感,而調(diào)壓室面積As2和As4以及調(diào)壓室流入系數(shù)r1僅對最高涌浪的影響較為明顯。同時,參數(shù)的不確定度并不會影響對調(diào)壓室涌浪的敏感性特性。
考慮到本文的主題是討論參數(shù)不確定框架下系統(tǒng)參數(shù)對輸出響應的敏感性特性,以及考慮到當前關于參數(shù)集引入數(shù)學模型展開敏感性分析的實際困難,這里沒有對特性曲線引起的模型不確定性進行分析,僅針對參數(shù)的不確定性進行了研究,在后續(xù)的工作中可以繼續(xù)考慮改進分析方法來展開模型不確定性的相關研究。□