謝卓健，閆懂林，肖 漢，麥先春，王峻峰，陳啟卷

（1.四川大唐國際甘孜水電開發(fā)有限公司，四川康定626001；2.武漢大學(xué)動力與機(jī)械學(xué)院，武漢430072）

0 引 言

在過去的二十年里，我國水電的開發(fā)和利用得到了快速的發(fā)展，在西南地區(qū)，各種大流量、高水頭、包含復(fù)雜水力系統(tǒng)的水電站不斷興起。到2019年底，我國的水電總裝機(jī)容量已經(jīng)達(dá)到3.56 億kW［1］。對于水電站來說，由于電網(wǎng)故障等原因，運(yùn)行中隨時可能會經(jīng)歷機(jī)組甩負(fù)荷過程，而一些水電站的重大事故，例如飛逸、壓力鋼管爆裂等，往往就是由于這種較為極端的工況引起［2］。如何有效地保證甩負(fù)荷過渡過程中機(jī)組的安全可靠，是實(shí)際水電工程中所面臨的重要問題。

到目前為止，針對水電機(jī)組甩負(fù)荷過渡過程進(jìn)行了大量研究。部分學(xué)者對實(shí)際電站發(fā)生的甩負(fù)荷事故進(jìn)行了分析和討論。黃菊珠［3］分析了山美水電站因110 kV 線路開關(guān)跳閘，造成3臺機(jī)組滿負(fù)荷高頻切機(jī)，引發(fā)機(jī)組過速停機(jī)的事故，并提出了相關(guān)的預(yù)防措施。曲東［4］分析了平班水電站在甩負(fù)荷試驗(yàn)時調(diào)速器調(diào)節(jié)失靈，導(dǎo)葉不能回關(guān)導(dǎo)致機(jī)組過速停機(jī)的事故，提出了相應(yīng)的處理方法，并對調(diào)速器重要參數(shù)進(jìn)行了優(yōu)化整定。張偉［5］針對西溪電站甩負(fù)荷試驗(yàn)時出現(xiàn)的機(jī)組過速停機(jī)及導(dǎo)葉快關(guān)時間不穩(wěn)定等問題進(jìn)行了分析，并提出了故障處理方法。除此之外，也有大量學(xué)者對水力發(fā)電機(jī)組甩負(fù)荷過渡過程進(jìn)行了數(shù)值模擬研究。付亮［6］針對有調(diào)壓井長引水隧洞水電站，利用數(shù)值模擬的方法研究了調(diào)壓井水位波動與流道壓力的變化過程。王威［7］等以伊朗卡侖水電站為例，分析了具有上游調(diào)壓室、一管兩機(jī)布置的水電站機(jī)組甩負(fù)荷過渡過程，并討論了導(dǎo)葉關(guān)閉規(guī)律對轉(zhuǎn)速上升和水擊壓力的影響，提出了合理的導(dǎo)葉關(guān)閉規(guī)律。葉文波等［8］針對南方某大型水電站，進(jìn)行了機(jī)組甩負(fù)荷大波動過渡過程的計(jì)算，分析了引水系統(tǒng)最高壓力升高和機(jī)組轉(zhuǎn)速的最大升高特性。陳勝等［9］研究了水電站中機(jī)組相繼甩負(fù)荷情況下調(diào)壓室涌浪的疊加特性。張洋［10］研究了設(shè)置氣墊式調(diào)壓室的超長引水隧洞水電站的甩負(fù)荷過渡過程，提出了氣墊式調(diào)壓室對系統(tǒng)水擊壓力的影響效應(yīng)。在這些分析和討論中，都是在確定性框架下研究水電站機(jī)組的甩負(fù)荷過渡過程，而并沒有考慮到參數(shù)不確定性對甩負(fù)荷過渡過程的影響。實(shí)際上，水電站過渡過程分析往往參數(shù)眾多且包含復(fù)雜的非線性動態(tài)過程，參數(shù)尺度多樣化及制造誤差等多種不確定性因素均會導(dǎo)致計(jì)算結(jié)果的不確定性。所以，在不確定框架下研究參數(shù)對甩負(fù)荷過渡過程的影響效應(yīng)就顯得更有必要。

參數(shù)敏感性分析是量化系統(tǒng)參數(shù)在設(shè)計(jì)或設(shè)定的可能范圍內(nèi)變化時對系統(tǒng)輸出響應(yīng)影響的一種方法［11］。通過敏感性分析，可以獲得對系統(tǒng)輸出影響較為顯著的輸入?yún)?shù)，無論在設(shè)計(jì)、運(yùn)行還是維護(hù)過程中，只需要重點(diǎn)關(guān)注這些參數(shù)而忽略一些不敏感的參數(shù)，就可以大大簡化工作量。通常的參數(shù)敏感性分析包括局部敏感性分析和全局敏感性分析，局部敏感性分析只能反映單一參數(shù)變化對系統(tǒng)響應(yīng)的影響，而全局敏感性分析則可以有效地反映多個參數(shù)同時變化對系統(tǒng)響應(yīng)的影響［12］。在全局敏感性分析方法中，主要包括：回歸分析法［13］、篩選法［14］，Sobol 指數(shù)法［15］和擴(kuò)展傅里葉幅度檢驗(yàn)法（extended Fouri?er amplitude sensitivity test，EFAST）［16］等。其中EFAST 方法可以獲得輸入?yún)?shù)之間的交互效應(yīng)，已經(jīng)在很多領(lǐng)域得到了廣泛應(yīng)用［17，18］。因此，本研究中將利用EFAST 全局敏感性分析方法來研究水力發(fā)電機(jī)組在甩負(fù)荷過渡過程中輸入?yún)?shù)對轉(zhuǎn)速、水擊壓力和調(diào)壓室涌浪的敏感性特性。

相比于之前的研究，本文的主要創(chuàng)新在于考慮參數(shù)的不確定性，在不確定框架下建立了水力發(fā)電機(jī)組的甩負(fù)荷過渡過程的模型，并以此模型為基礎(chǔ)進(jìn)行了數(shù)值研究，量化分析了輸入?yún)?shù)對機(jī)組甩負(fù)荷過程中機(jī)組轉(zhuǎn)速、水擊壓力和調(diào)壓室涌浪的敏感性特性。相關(guān)的理論方法和研究結(jié)論可以應(yīng)用于水電站的設(shè)計(jì)及維護(hù)工作。

1 水力發(fā)電機(jī)組甩負(fù)荷過渡過程建模

水電站的引水系統(tǒng)在空間分布上主要包括上游水庫，引水隧洞，上游引水調(diào)壓室，壓力鋼管，蝸殼，尾水管及尾水明渠段。本節(jié)將首先利用特征線方法建立該電站引水系統(tǒng)的數(shù)學(xué)模型，同時結(jié)合水輪機(jī)特性曲線的插值模型，調(diào)速器數(shù)學(xué)模型，發(fā)電機(jī)及負(fù)載數(shù)學(xué)模型，建立完整的水電機(jī)組甩負(fù)荷過渡過程的數(shù)學(xué)模型。

1.1 引水系統(tǒng)模型

1.1.1 有壓過流管道特征線方程

利用特征線法將有壓輸水系統(tǒng)非恒定流的運(yùn)動方程和連續(xù)性方程轉(zhuǎn)化為在特征線dx/dt=V±a上的常微分方程組［19］，可以得到：

式中：H表示測壓管水頭；a表示水擊波速；f為達(dá)西-維斯巴哈摩擦系數(shù)；D表示管道直徑；A為流道截面面積；Q為流道中的流量。在這里，V和H是與時間和空間相關(guān)的函數(shù)，其中t為時間，x表示管道長度。

在利用特征線法離散化求解微分方程時，需要對有壓管道按照等距離進(jìn)行劃分，得到長度為Δx的N段管，并用序號i進(jìn)行標(biāo)記，定義起始端為i=1，則末端為i=N+1。設(shè)定仿真步長為ΔT=Δx/a。利用特征線法進(jìn)行計(jì)算機(jī)數(shù)值求解的原理如圖1所示，橫坐標(biāo)表示管道的空間分布，縱坐標(biāo)表示采樣時刻分布。從中可以看出，根據(jù)上一時刻的節(jié)點(diǎn)i-1 和節(jié)點(diǎn)i+1 的值，可以計(jì)算出下一時刻節(jié)點(diǎn)i的水頭、流量值。

圖1 特征線法原理示意Fig.1 Schematic diagram of characteristics method

對公式（1）和（2）的特征線方程積分并離散化，可以得到t+ΔT時刻節(jié)點(diǎn)i處的水頭、流量關(guān)系方程式為：

式中：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 水輪機(jī)邊界

水輪機(jī)的瞬態(tài)工作水頭定義為水輪機(jī)蝸殼進(jìn)口和尾水管出口的水頭差，如圖2所示，當(dāng)水輪機(jī)上、下游節(jié)點(diǎn)為1 和2 時，則水輪機(jī)工作水頭為HT=H1-H2，瞬態(tài)的過機(jī)流量為QT=Q1=Q2。

圖2 水輪機(jī)簡化結(jié)構(gòu)圖Fig.2 Structure chart of hydro-turbine

根據(jù)水輪機(jī)上下游進(jìn)出口節(jié)點(diǎn)的水頭平衡方程，可以得到：

水庫邊界，上游調(diào)壓室邊界，蝸殼及尾水管當(dāng)量管的計(jì)算可以參考文獻(xiàn)［19］。

1.2 水輪機(jī)特性曲線插值模型

依據(jù)文獻(xiàn)［20］的方法，對原始的水輪機(jī)綜合特性曲線進(jìn)行擴(kuò)展，可以得到Q11=f1（n11，y）和M11=f2（n11，y）的關(guān)系，用于甩負(fù)荷過渡過程的計(jì)算。在這兩個關(guān)系式中，n11表示單位轉(zhuǎn)速，M11表示單位力矩，y表示導(dǎo)葉開度。以當(dāng)前研究的水電站為例，得到擴(kuò)展后的Q11-n11-y和M11-n11-y關(guān)系曲線如圖3所示。

圖3 單位流量和單位力矩圖Fig.3 Diagrams of unit discharge and unit moment

1.3 發(fā)電機(jī)及負(fù)載數(shù)學(xué)模型

水力發(fā)電機(jī)組的力矩平衡方程［19］可以表示為：

式中：J為機(jī)組轉(zhuǎn)動部分的轉(zhuǎn)動慣量，J=GD2/4，t·m2；GD2為飛輪力矩，t·m2；ω表示機(jī)組旋轉(zhuǎn)角速度，rad/s；M和Mg分別是水輪機(jī)的主動力矩和阻力矩，kN·m。

假定水輪機(jī)力矩在ΔT時刻里線性變化，則可以利用梯形積分公式對上式積分。同時，將旋轉(zhuǎn)角速度（ω）轉(zhuǎn)換為轉(zhuǎn)速（n），可以得到：

至此，水力發(fā)電機(jī)組甩負(fù)荷過渡過程的確定性模型被建立，在這個確定性模型中將確定性的參數(shù)用不確定的分布或者范圍進(jìn)行替代，即可以獲得不確定性的模型，以此為基礎(chǔ)，相關(guān)的參數(shù)敏感性可以利用方便的非侵入式EFAST方法進(jìn)行討論。

2 全局敏感性分析

2.1 全局敏感性分析方法

本文采用EFAST 方法來對水電機(jī)組甩負(fù)荷過渡過程的參數(shù)敏感性進(jìn)行分析。EFAST 法的原理是在輸入?yún)?shù)的多維空間中挑選合適的搜索曲線，并給所有輸入?yún)?shù)分配非線性相關(guān)的整數(shù)頻率值，從而引入搜索函數(shù)，使模型作為獨(dú)立參數(shù)的周期函數(shù)，則可以把多維積分問題轉(zhuǎn)換為一維積分問題。然后通過把目標(biāo)函數(shù)變?yōu)楦道锶~級數(shù)，計(jì)算出各頻率的傅里葉頻譜曲線，由頻譜曲線獲得由輸入?yún)?shù)xi引起的系統(tǒng)輸入方差和總方差，其比值則被定義為該參數(shù)的敏感度指標(biāo)。文獻(xiàn)［16-18］中給出了具體的算法推演流程，這里僅給出敏感度計(jì)算的公式。

總方差可以寫成各個子方差和的形式：

式中：V為模型總方差；Vi表示單一輸入變化引起的模型方差；Vij表示通過xi和xj交互作用引起的方差；Vij…k表示xi，xj，…，xk交互作用引起的耦合方差。

主敏感性（一階敏感性）指標(biāo)表示單個參數(shù)獨(dú)立作用時對系統(tǒng)輸出響應(yīng)的貢獻(xiàn)，可以表示為：

總敏感性指標(biāo)為：

式中：V-i表示不包括參數(shù)xi的所有參數(shù)方差和。

2.2 仿真設(shè)計(jì)與參數(shù)選擇

以國內(nèi)某水電站為例，進(jìn)行機(jī)組甩100%負(fù)荷時過渡過程計(jì)算的參數(shù)敏感性分析，涉及的主要參數(shù)如表1所示。敏感性分析的輸入?yún)?shù)均值選取設(shè)計(jì)值，方差分別取變異系數(shù)（=標(biāo)準(zhǔn)差/均值）為0.01和0.05兩種情況式的值，代表不同的不確定度。選擇Monte Carlo 方法進(jìn)行參數(shù)的隨機(jī)采樣，采樣樣本為5 000，滿足EFAST 法對樣本個數(shù)的要求。另外，由于其他水工建筑物的影響，該電站不同高程下的調(diào)壓室面積會有所有差別，在高程Hs0<1 458.7 m 時，面積As2為801.975 m2；當(dāng)1 458.7 m1 476.5 m 時，面積As4為2 562 m2。

表1 仿真參數(shù)Tab.1 Simulation parameters

3 分析與討論

在機(jī)組的甩負(fù)荷過程中，最高轉(zhuǎn)速、最大水擊壓力、調(diào)壓室最大涌浪和最小涌浪需要被重點(diǎn)關(guān)注，也是評估機(jī)組甩負(fù)荷過程安全性的重要指標(biāo)，所以，在本節(jié)中，參數(shù)相對于這幾個輸出變量的敏感性將被重點(diǎn)分析。

3.1 轉(zhuǎn)速升高

機(jī)組甩負(fù)荷過程中最高轉(zhuǎn)速的全局敏感性分析結(jié)果如圖4所示。

圖4 最高轉(zhuǎn)速的全局敏感性Fig.4 Global sensitivity of maximum speed

以甩負(fù)荷過程中的機(jī)組最高轉(zhuǎn)速作為輸出時，無論在變異系數(shù)為0.01 還是0.05 的情況下，機(jī)組旋轉(zhuǎn)部分的飛輪力矩都是對其最為敏感的參數(shù)，其他參數(shù)的影響微乎其微。另外，對比飛輪力矩的主效應(yīng)和總效應(yīng)指標(biāo)，可以發(fā)現(xiàn)在不同的不確定度下其差值都很小，即飛輪力矩對甩負(fù)荷過程中轉(zhuǎn)速升高的影響也主要是通過自身的作用實(shí)現(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)速影響的放大尺度，從而導(dǎo)致不同不確定度下會出現(xiàn)這種差異。當(dāng)然從整體的角度來說，飛輪力矩的絕對主導(dǎo)地位仍然不會受到威脅?？傮w來說，飛輪力矩對于甩負(fù)荷過程中的機(jī)組轉(zhuǎn)速升高的影響最為顯著。

3.2 水擊壓力

除了轉(zhuǎn)速的升高，水擊壓力也是甩負(fù)荷過程中影響機(jī)組安全的關(guān)鍵指標(biāo)。如圖5所示，給出了機(jī)組甩負(fù)荷過程中最高水擊壓力的全局敏感性分析結(jié)果。

圖5 最高水擊壓力的全局敏感性Fig.5 Global sensitivity of maximum water hammer pressure

當(dāng)變異系數(shù)為0.01 時，輸入?yún)?shù)對水擊壓力敏感性的主效應(yīng)指標(biāo)排序?yàn)椋篈2>a1>As4>r1>As3>a2>a4>a3>A5>f1>f2>As1>A4>A3>GD2>f5>f3>A1>a5>As2>f4>r2；輸入?yún)?shù)對水擊壓力敏感性的總效應(yīng)指標(biāo)排序?yàn)椋篈2>a3>As4>As3>A5>a1>f1>a2>f5>r1>a4>f2>a5>GD2r2>As1>A3>A1>As2>f4>A4>f3。無論是主效應(yīng)還是總效應(yīng)，壓力鋼管面積對水擊壓力的敏感性最為顯著，且從中可以看出A2獨(dú)立的影響占46.6%，與其他參數(shù)的交互效應(yīng)占53.4%，也就是說A2對水擊壓力的獨(dú)立作用和與其他參數(shù)的交互作用效果基本持平。此外，蝸殼水擊波速a3，調(diào)壓室面積As4和As3的總效應(yīng)指標(biāo)都超過了0.4，它們對最高水擊壓力的影響也僅次于A2。與A2的作用方式不同的是，a3、As4和As3都是主要通過與其他參數(shù)的交互作用來實(shí)現(xiàn)對水擊壓力的影響。

當(dāng)變異系數(shù)為0.05 時，輸入?yún)?shù)對水擊壓力敏感性的主效應(yīng)指標(biāo)排序?yàn)椋篴3>A2>r1>A4>As3>A1>A3>a4>f1>GD2>A5>As1>As4>As2>a5>a2>a1>f2>f3>f5>f4>r2；輸入?yún)?shù)對水擊壓力敏感性的總效應(yīng)指標(biāo)排序?yàn)椋篴3>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獨(dú)立作用對水擊壓力影響的占比也從變異系數(shù)為0.01 時的11.51%增加到49.54%，這也就意味著隨著系統(tǒng)不確定度的增加，a3獨(dú)立的影響會被放大。其主要原因是系統(tǒng)不確定度的增加會使得參數(shù)的變化范圍也相應(yīng)地?cái)U(kuò)大，在這個過程中，a3在偏離標(biāo)準(zhǔn)值越遠(yuǎn)時其單參數(shù)項(xiàng)對水擊壓力影響的變化程度更加顯著，從而使得其獨(dú)立影響的占比在不確定度增加后提升而交互效應(yīng)的占比下降。此外，相比于變異系數(shù)為0.01 時的敏感性指標(biāo)，As4和As3的影響在變異系數(shù)為0.05時會被削弱。值得注意的是，在變異系數(shù)為0.05時，調(diào)壓室流入系數(shù)r1和蝸殼截面積A4對最大水擊壓力的影響會顯著提升。

3.3 調(diào)壓室涌浪

圖6描述了機(jī)組甩負(fù)荷過程中調(diào)壓室最高涌浪的全局敏感性分析結(jié)果。在變異系數(shù)為0.01 時，調(diào)壓室截面積（As2、As3和As4），調(diào)壓室流入系數(shù)（r1）以及引水隧洞的面積（A1）對最高涌浪較為敏感，它們之間的總敏感性排序?yàn)椋篈s3>A1>As4>r1>As2>其他。同時，可以注意到參數(shù)As3，A1，As4和r1主要是通過它們獨(dú)立的作用直接影響最高涌浪高度，而參數(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ù)敏感性指標(biāo)進(jìn)行對比，可以發(fā)現(xiàn)它們之間的差異并不明顯。也就是說，參數(shù)不確定性的變化對調(diào)壓室最高涌浪的影響比較微弱，可以忽略不計(jì)。

圖6 調(diào)壓室最高涌浪的全局敏感性Fig.6 Global sensitivity of highest surge

圖7給出了機(jī)組甩負(fù)荷過程中調(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

4 結(jié) 論

本文以國內(nèi)某水電站為例，建立了水力發(fā)電機(jī)組甩負(fù)荷過渡過程的數(shù)學(xué)模型，并考慮到參數(shù)的不確定性，引入EFAST 方法對輸入?yún)?shù)在不同不確定度下進(jìn)行了敏感性分析，可以得到以下結(jié)論。

（1）對于機(jī)組甩負(fù)荷過程中轉(zhuǎn)速的最大上升值，機(jī)組的飛輪力矩對其影響占據(jù)了絕對的主導(dǎo)地位，并且這種影響效應(yīng)是飛輪力矩獨(dú)立作用的效果而不是與其他參數(shù)交互影響的結(jié)果，且參數(shù)的不確定度對參數(shù)敏感性的影響不明顯。

（2）在較小的不確定度下，壓力鋼管面積A2，蝸殼波速a3，調(diào)壓室面積As4和As3，以及尾水管面積A5是對水擊壓力影響最為顯著的5 個參數(shù)。而當(dāng)參數(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ù)對輸出響應(yīng)的敏感性特性，以及考慮到當(dāng)前關(guān)于參數(shù)集引入數(shù)學(xué)模型展開敏感性分析的實(shí)際困難，這里沒有對特性曲線引起的模型不確定性進(jìn)行分析，僅針對參數(shù)的不確定性進(jìn)行了研究，在后續(xù)的工作中可以繼續(xù)考慮改進(jìn)分析方法來展開模型不確定性的相關(guān)研究?！?/p>