水力發電機組軸系不確定性量化及參數敏感性分析

閆懂林，鄭 陽，劉宛瑩，陳啟卷

(武漢大學 動力與機械學院，武漢 430072)

水力發電機組作為水電站中能量轉換的樞紐，其故障或事故，會直接影響水電站乃至電網的安全穩定[1]。因此，為了避免一些重大安全事故地發生，對水力發電機組結構特性和動力學響應進行研究十分必要。

到目前為止，已經有大量學者致力于水力發電機組結構建模及動力學特性的分析研究。Xu等[2]綜合考慮了阻尼、油膜和磁拉力對水力發電機組的影響，以拉格朗日方程為基礎建立了一個機組軸系的分數階非線性動力學模型，并討論了參數變化對機組非線性特性的影響效應；Huang等[3]提出了包含不對中和質量偏心故障的水力發電機組轉子系統數學模型，并分析了故障耦合作用下的系統動力學特性；Zhang等[4]推導了不同荷載條件下的水電機組水力激勵的數學表達式，研究了機組突增負荷過程中軸系的振動響應特性；Yan等[5]建立了包含葉片的貫流式水力發電機組軸系統模型，討論了葉片對機組軸系振動特性的影響機理?？偟膩碚f，這些研究可以為水力發電機組的設計、優化、運行和維護提供一定的參考。但值得注意的是，當前的研究主要都是基于確定性模型和方法展開，很少提及參數不確定性的影響。實際上，對于水力發電機組來說，存在許多不確定性因素，如制造公差、材料特性的固有隨機性、磨損、應力變化、隨機載荷等，這些不確定性對機組的安全和可靠性會造成潛在的威脅[6]。而現有的確定性分析方法對水力發電機組的不確定性問題可能是局限的或無效的，因此，為了探究不確定參數對機組結構性能和動力學響應的影響，需要引入不確定性分析的方法和理論來展開研究。同時，考慮到結構的不確定性較難發現，其破壞性在實際工程中也最為顯著，對其研究也相對較少，故本研究將重點研究機組結構參數的不確定性。

大量研究表明，不確定性量化和參數敏感性分析是探索物理系統不確定問題的兩種有效方法[7]。一方面，不確定性量化是用來評估系統輸入到輸出過程中不確定性的傳播，通常采用輸出響應的均值、方差和置信區間作為量化指標[8]；另一方面，參數敏感性分析是為了量化不確定參數對系統輸出不確定性的貢獻，可以分為局部敏感性分析和全局敏感性分析。當模型的輸出不確定性比較突出，或者模型具有較強的非線性時，全局敏感性分析比局部敏感性分析更具優勢[9]，主流的全局敏感性分析方法有回歸分析法、Morris篩選法、Sobol法、傅里葉幅度檢驗法和擴展傅里葉幅度檢驗法等?？紤]到水力發電機組軸系統是一個復雜的非線性系統，這里將利用全局敏感性分析來討論機組結構參數不確定對系統輸出的影響效應。

對于高自由度的水力發電機組轉子結構系統，在計算不確定性參數的各階敏感性指標時，傳統的暴力Monte-Carlo模擬法需要大量的樣本，計算成本很高甚至難以接受。為了解決這一問題，本文將引入基于多項式混沌展開的代理模型技術[10]，建立描述機組不確定結構參數與輸出變量之間關系的代理模型用于不確定性量化研究。同時，由于多項式混沌展開的形式可以很容易地轉化為方差分解的形式，只需對展開系數進行簡單的后處理便可以得到參數的敏感性指標，避免了大量的樣本計算，故本文將采用這種方法來對機組結構參數不確定性的影響展開研究。

綜上所述，本文提出了一套以廣義多項式混沌展開為基礎的水力發電機組軸系建模、不確定性量化和參數敏感性分析的統一框架。主要貢獻如下：① 將水力發電機組軸系結構特性和動力響應的研究擴展到不確定性范圍，并建立了一個包含不確定性信息的機組軸系動力學模型；② 對于不確定參數下的機組結構固有特性和振動響應的不確定性量化，除了廣泛使用的均值、方差和各階統計矩外，還利用最大熵原理得到了隨機輸出的具體概率密度函數表達式；③ 基于廣義多項式混沌方法，快速獲得了不確定參數對固有頻率(靜態響應)和振動響應(動態響應)的全局靈敏性指標，指導水電機組的實際生產實踐。

1 水力發電機組軸系建模

水力發電機組的軸系簡化結構如圖1所示，其動力響應主要由軸系、邊界條件和外部激勵決定，下面將分別建立各子單元的模型，并將其組成完整的軸系動力學模型。

圖1 機組簡化結構圖

1.1 軸系統

在軸系中，主要包括轉軸、發電機轉子和水輪機轉輪，這里將采用有限元法建立它們的模型。

對于轉軸，采用2節點8自由度的Timoshenko梁單元進行建模，其質量矩陣M(e)，剛度矩陣K(e)，和陀螺矩陣G(e)見文獻[11],下標(e)為單元矩陣。梁單元每個節點有兩個平移自由度(x,y)和兩個旋轉自由度(θx,θy)，對于第i個單元，其坐標可以寫成ue=[x(i),y(i),θx(i),θy(i),x(i+1),y(i+1),θx(i+1),θy(i+1)]。

發電機轉子和水輪機轉輪的質量矩陣分別為

Mrotor=diag{m1,m1,Jd1,Jd1}

(1)

Mrunner=diag{m2+mwater,m2+mwater,Jd2,Jd2}

(2)

式中：m1和m2分別為發電機轉子和水輪機轉輪的質量；Jd1和Jd2分別為它們的截面慣性矩；mwater為水體的附加質量。

發電機轉子和水輪機轉輪的陀螺矩陣可以表示為

(3)

(4)

式中，Jp1和Jp2分別為轉子和轉輪的極慣性矩。

根據這些子單元的矩陣，可以構造出機組軸系的完整質量、剛度和陀螺矩陣，構造過程如圖2所示。圖2中從左上角到右下角的6個節點分別為：上導軸承、發電機轉子、下導軸承、聯軸器、水導軸承及水輪機轉輪。M(e)iL和M(e)iR為單元質量矩陣M(e)的左上角和右下角子矩陣，G(e)iL和G(e)iR為單元陀螺矩陣G(e)的左上角和右下角子矩陣，n1～n5為疊加的單元矩陣個數。

圖2 軸系整體矩陣的裝配

最后，考慮到水力發電機組是一個多自由度系統，其阻尼采用瑞利阻尼理論進行假設，具體的計算方法和公式見文獻[12]。

1.2 邊界條件

機組軸系的邊界條件主要為各個導向軸承提供的維持軸系旋轉的支撐力，以油膜力[13]的形式表示

(5)

式中：Fx-oil和Fy-oil分別為x和y方向的油膜力；Fx0和Fy0為靜止工作點的油膜力；ΔFx-oil和ΔFy-oil為軸頸偏離靜止工作點引起的油膜力增量，其可以用式(6)計算

(6)

式中：kij(i,j=x,y)為剛度系數；dij(i,j=x,y)為阻尼系數；x,y為軸向位移。

1.3 外部激勵

外部激勵對水力發電機組的動態響應有明顯的影響，在本研究中，考慮了質量偏心和不平衡磁拉力兩種最典型的激勵。具體的模型表達式見Xu等和Zhang等的研究。

1.4 運動方程

基于1.1節～1.3節各個子模型，機組在確定性框架下的軸系動力學模型可以寫為

(7)

考慮到機組結構參數的不確定性，在不確定性框架下，水力發電機組的運動方程可表示為

(8)

2 廣義多項式混沌展開

2.1 展開式

假設Y=Y(ξ)為一個隨機物理模型，ξ為輸入隨機變量的集合,ξ= {ξ1,ξ2, …,ξd}，其中的ξi為第i個輸入參數，則該模型的輸出響應可以用多項式混沌方法重構[14]為

(9)

式中：Y為系統的輸出響應；d為隨機輸入的個數；b0,bi1,i2,…為展開系數；Φn(·)定義了n階多項式;ψ0為基函數；i1，i2,…為各參數展開項編號。

對式(9)中各項進行重新排序，并進行有限項截斷，則原模型可以近似表示為

(10)

式中：r為截斷維度；Nd為多維索引的集合；α= {α1, …,αd|αi≥0}為Nd中包含d個元素的多維指標；Ad,r為一個根據r和d定義的區域；βα為重排后的多項式系數；ψα(ξ)為正交基函數，對于有多個隨機輸入的問題，基函數可通過張量積[15]進行構造。

截斷后的多項式混沌展開項數由截斷維度和隨機變量個數決定，用式(11)計算

(11)

當多項式混沌展開項數確定后，式(10)可以寫成一個更緊湊的形式

(12)

式中：β為展開系數的向量，β= {β0, …,βN-1}T；ψ(ξ)為基函數的向量，ψ(ξ)= {ψ0(ξ), …,ψN-1(ξ)}；βj為整理后的展開系數；ψj(ξ)為對應的基函數。

2.2 展開系數

假定Yn= {Y1,Y2, …,Yj,…,Yn}T為系統在n個樣本ξ= {ξ1,ξ2, …,ξn}下的輸出矢量，各實際值與估算值之差ε越小，估計越準確，所以對多項式混沌展開系數的總體估計可以轉化為如下的最小化問題

(13)

式中：Yj為第j個樣本下的系統輸出響應值；ξj為第j個輸入樣本。

式(13)的最優解即多項式混沌展開系數可以用式(14)計算

β=(HTH)-1HTYn

(14)

式中：Yn為n個樣本下的系統輸出矢量；H為輔助計算矩陣，其可以表示為

(15)

3 不確定性量化和全局敏感性分析

3.1 不確定性量化

最大熵原理表明：為了獲得具有部分已知信息事物的概率分布形式，在約束條件下，選擇熵最大的概率分布是最優的。在本研究中，對于原點矩確定的隨機輸出響應，有多個不同的概率密度函數與之對應，其中最優概率密度函數是熵值最大的一個。換句話說，對不確定系統隨機輸出的概率密度函數的求解可以轉換為如下的優化問題

(16)

式中：H為熵；f(x)為隨機變量x的概率密度函數；μs為s階原點矩；下標s為原點矩的階次；Ps(x)為已知約束條件下的基函數。

利用拉格朗日乘子法求解式(16)，可以得到最優的概率密度函數表達式為

(17)

式中：ηj為待定系數；Pj(x)為第j階的函數;S為展開的最大階次。

Abramov的研究表明，冪函數Ps(x)=xs可以作為基函數用于求解概率密度函數，且四階原點矩已經可以保證足夠的精度，故輸出響應的概率分布函數可以寫成

(18)

3.2 全局敏感性分析

定義一個在Ad,r中的多維指標Li1,…,ie為

Li1,…,ie=

(19)

基于式(19)，截斷的多項式混沌展開式式(10)可以被重新組織為標準的方差分解形式

(20)

所以，基于方差的各階敏感性指標可由式(21)直接計算得到

(21)

式中：V(Y)為輸出的總方差；Sik為主效應指標；Si1,…,ie為交互效應指標。

而總的敏感性指標為

(22)

4 數值研究

包含不確定參數的水力發電機組軸系固有特性及動力學響應不確定性量化及全局敏感性分析的具體流程，如圖3所示。依據文獻[16]，軸系不確定性參數和必要的確定性參數分別如表1和表2所示。這里選擇d= 7和r= 3可以滿足多項式混沌的估計效果，對應的樣本數設置為240，采用拉丁超立方抽樣確定樣本。

表1 不確定結構參數及其分布

表2 確定性參數

圖3 數值研究流程圖

4.1 驗 證

隨機選擇120個樣本，對比原始模型和代理模型下機組的振動響應來驗證多項式混沌展開代理模型的準確性，對比結果如圖4所示。從圖4可知，多項式混沌展開模型可以很好地反映原始模型中不確定參數與隨機輸出的關系。

(a)轉子

4.2 對固有頻率的不確定性量化和參數敏感性分析

基于暴力Monte-Carlo模擬和最大熵原理計算的機組軸系固有頻率的概率密度分布如圖5所示。從圖5可知，兩種方法下的概率分布是高度一致的。但其中暴力Monte-Carlo模擬是基于10 000個樣本，而最大熵原理僅僅基于240個樣本。換句話說，多項式混沌展開驅動的最大熵原理方法可以在大大縮減計算量的情況下有效地評估不確定參數影響下水力發電機組固有頻率的概率分布特性，且能給出連續的概率密度函數表達式，具體的表達式系數如表3所示。

(a)一階固有頻率

表3 固有頻率的概率密度函數系數

基于固有頻率的概率密度函數，可以得到如下結果：

(1)比較確定性框架下的固有頻率響應值(determinate value, DV)和不確定性框架下的均值(mean value, MV)，可以發現它們在一階固有頻率下的響應值非常接近。對比三階固有頻率，它們之間的差異被放大。而對比五階固有頻率，這種差異比一階和三階固有頻率的差異更明顯。對于高階固有頻率，其值在不確定系統中更容易偏離確定性系統的固有頻率值。

(2)根據計算的概率密度函數，可以得出一階、三階和五階固有頻率在95%置信區間下的輸出響應范圍分別為[14.733 8,18.368 9]，[45.584 3,55.673 3]和[88.123 7,107.056 8]。由于充分考慮了結構的不確定性，這些范圍對機組的設計優化比確定的值更有意義。

(3)根據各固有頻率的多項式混沌表達式，經過簡單后處理，得到各固有頻率的方差值分別為0.850 3，6.625 1和23.196 3，隨著固有頻率階數地增加而增大。換句話說，結構不確定性對高階固有頻率的影響更為顯著。

下面將利用全局敏感性分析來研究每個結構參數對固有頻率不確定性的具體貢獻。各參數對轉軸固有頻率敏感性的主效應和總效應指標在廣義多項式混沌方法和傳統暴力Monte-Carlo模擬下被分別提出，對比發現當基于廣義多項式混沌方法計算敏感性指標時在更小的計算成本下得到的精度結果與使用大量樣本的暴力Monte-Carlo模擬相同，如圖6所示。

基于圖6，一些關于固有頻率參數敏感性的結論可以被得到：首先，一階固有頻率對參數E和Jd1是最為敏感的，也就是說，一階固有頻率的不確定性主要來源于這兩個參數;同時，m2對一階固有頻率也有輕微的影響，而其他不確定參數的影響可以忽略不計。各參數對一階固有頻率的總效應指標與主效應指標非常接近，即各參數與其他參數的交互效應對一階固有頻率的影響不顯著。其次，對于三階固有頻率,參數E仍然是最主要的不確定性來源，參數Jd2，ρs和Gs次之。而參數m1，m2和Jd1對三階固有頻率不確定性的貢獻接近于零。此外，對于五階固有頻率，有3個顯著的靈敏度參數：ρs，E和Gs，其中參數ρs對五階固有頻率的影響最為顯著，而參數m1，m2和Jd2對五階固有頻率不敏感。最后，從整體的角度來說，無論是一階、三階還是五階固有頻率，參數E對它們的不確定性的貢獻都是比較明顯的。

(a)一階固有頻率

4.3 對振動響應的不確定性量化和參數敏感性分析

水力發電機組作為一種旋轉機械，其振動響應也是實際工程中評價其運行狀態的重要指標。因此，本節將繼續基于不確定性量化和參數全局敏感性分析過程研究結構參數不確定對振動響應的影響。各個導軸承和聯軸器處隨機振動響應的概率分布，如圖7所示?；谧畲箪卦淼母怕拭芏群瘮当磉_式系數，如表4所示。通過對比可以發現在暴力Monte-Carlo模擬和最大熵原理下的振動響應概率分布是一致的。

表4 振動響應的概率密度函數系數

根據圖7中機組不同位置振動響應的隨機分布可知，無論在什么位置，不確定框架下的振動響應均值都非常接近于確定系統的響應值。即確定性模型可以有效地描述不確定性框架下振動響應的均值信息。同時，在最大概率下，聯軸器和水導軸承處的振動響應值會明顯偏離均值，而在上導軸承和下導軸承處，這一現象并不明顯。其次，考慮95%的置信區間，可以得到軸系不同位置處隨機振動響應的范圍，分別為[4.562×10-5,7.820×10-5],[2.736×10-5,4.668×10-5],[3.054×10-5, 6.286×10-5]和[3.765×10-5,8.047×10-5]對應于上導軸承、下導軸承、耦合器和水導軸承，覆蓋長度分別為3.257×10-5,1.932×10-5,3.231×10-5,4.281×10-5。另外，不同位置處隨機振動響應的方差值分別為6.947 3×10-11,2.448 3×10-11,6.754 5×10-11和1.182 0×10-10。結合隨機振動響應的覆蓋長度和方差值，可以得到不同位置處的不確定度排序為：水導軸承>上導軸承>聯軸器>下導軸承。

(a)上導軸承

在完成振動響應的不確定性量化分析后，下面將通過全局敏感性分析確定各不確定參數對振動響應不確定性的貢獻。不同位置處振動響應的全局靈敏度指標，如圖8所示，其中包括主效應指數和總效應指數。對比兩種方法的結果，基于多項式混沌方法計算振動響應敏感性指標的可靠性同樣得到了驗證。

根據圖8中關于振動響應的全局靈敏度分析，可以得到：

(1)對于振動響應而言，無論在什么位置，各不確定結構參數的主效應指標都非常接近于總效應指標，即不確定結構參數對振動響應的交互影響不明顯，這與固有頻率的敏感性特性是相似的。

(2)不同位置下不確定參數的敏感性指標存在明顯差異。具體來說，對于上導軸承處的振動響應，最敏感的參數是m1，其次是E和Gs，最小的參數是ρs，m1，Jd1和Jd2，其中Jd1和Jd2的影響可以忽略不計。與上導軸承處的敏感性特性相比，參數m1對下導軸承處振動不確定性的貢獻明顯減小，而參數E和Gs的影響明顯增大。此外，聯軸器與水導軸承處振動響應的靈敏度指標相似，可以排序為E>m2>Gs>ρs>其他。

(3)在4.2節中，得到參數Jd1和Jd2對固有頻率有明顯的貢獻。但從圖8可以看出，參數Jd1和Jd2在任何位置對振動響應的影響都不顯著或可以忽略。這表明結構參數對固有頻率(靜態指標)的敏感性與對振動響應(動態指標)的敏感性并不直接相關。

(a)上導軸承

5 結 論

本文提出了一種基于廣義多項式混沌方法對水力發電機組軸系進行不確定性量化和全局敏感性分析的統一框架，并通過與暴力Monte-Carlo模擬結果的對比分析，驗證了其可靠性。以此為基礎，考慮結構參數的不確定性，對機組固有頻率和振動響應的不確定性和參數全局敏感性特性進行了分析，得到了如下的主要結論：

(1)得到了結構參數不確定下機組軸系固有頻率和振動響應概率密度函數的具體表達式、均值、方差和置信區間。

(2)結構參數不確定性對不同輸出的傳播特性表現出明顯差異。在相同的不確定輸入下，五階固有頻率比一階和三階固有頻率具有更大的不確定性，而不同位置下機組振動響應的不確定度表現為——水導軸承>上導軸承>聯軸器>下導軸承。

(3)參數對固有頻率(靜態指標)和振動響應(動態指標)的靈敏度特征不直接相關。例如，參數Jd1和Jd2只影響固有頻率，而不影響振動響應。同時，各階固有頻率對參數的靈敏度也沒有直接關系。

(4)參數E無論對各階固有頻率還是不同位置下的振動響應都保持較強的靈敏度，在實際工程應用中需要重點關注。

由于考慮了結構參數的不確定性，上述結果比基于確定性模型的結論對機組的設計、優化、運行和維護更有指導意義。同時，所提出的不確定性量化和全局敏感性分析框架也可以推廣到水力發電機組其他不確定參數的分析研究當中。當然，本文主要是針對結構參數的不確定性展開研究，而沒有考慮其他參數不確定性的影響，在后續的工作中可以進一步討論其他類型不確定性參數對機組響應的影響效應。