基于二階攝動(dòng)的水電機(jī)組振動(dòng)傳導(dǎo)研究

職保平， 王 宇， 秦凈凈， 于 洋， 張宏戰(zhàn)

(1.小流域水利河南省高校工程技術(shù)研究中心， 河南 開封 475004；2.西藏自治區(qū)水利電力規(guī)劃勘測(cè)設(shè)計(jì)研究院，拉薩 850000； 3. 大連理工大學(xué)，遼寧 大連 116024)

在水電機(jī)組振動(dòng)分析中，振動(dòng)傳導(dǎo)的分析對(duì)電站安全穩(wěn)定性分析具有至關(guān)重要的意義，目前水電機(jī)組的振動(dòng)分析集中于振源的發(fā)生作用機(jī)理、響應(yīng)的作用機(jī)制的研究，在耦合振源及振動(dòng)體系的研究，但鮮見于振動(dòng)傳導(dǎo)過程以及放大作用機(jī)制的研究。

水電機(jī)組的振動(dòng)研究，由于結(jié)構(gòu)巨型化和復(fù)雜化產(chǎn)生的制造安裝誤差以及結(jié)構(gòu)材料參數(shù)的不確定性是客觀存在的[1]，整個(gè)水電站振動(dòng)系統(tǒng)屬于典型的隨機(jī)結(jié)構(gòu)系統(tǒng)問題。在處理隨機(jī)參數(shù)結(jié)構(gòu)振動(dòng)問題時(shí)，現(xiàn)有的分析方法大多局限于一階區(qū)間分析方法，如利用Kronecker代數(shù)[2]、Hardmarda積[3]進(jìn)行一階攝動(dòng)概率、區(qū)間分析結(jié)構(gòu)振動(dòng)[4]問題，利用非概率分析方法求解結(jié)構(gòu)動(dòng)力響應(yīng)的區(qū)間參數(shù)[5]，以及基于模態(tài)疊加法與一階攝動(dòng)法分析結(jié)構(gòu)不確定荷載的影響[6]等。一階分析方法中，假設(shè)隨機(jī)參數(shù)具有小擾動(dòng)，但并不適用于參數(shù)變化范圍較大的情況，采用一階區(qū)間分析方法對(duì)結(jié)構(gòu)動(dòng)力響應(yīng)范圍進(jìn)行估計(jì)可能會(huì)失效[7]。二階攝動(dòng)的研究中，目前從解動(dòng)力方程等的角度進(jìn)行分析，即將二階攝動(dòng)轉(zhuǎn)化為一階問題進(jìn)行分析，目前在該類方法上取得較多成果[8-9]。

由于水電機(jī)組中存在上導(dǎo)軸承、下導(dǎo)軸承、推力軸承等含有油箱的結(jié)構(gòu)，其油膜剛度的變化范圍較大，采用一階攝動(dòng)方法分析時(shí)，將存在隱患，因此，在分析一階攝動(dòng)法的基礎(chǔ)上，提出采用二階攝動(dòng)法進(jìn)行水電機(jī)組振動(dòng)傳導(dǎo)分析，明確二階攝動(dòng)，特別是含二階偏導(dǎo)數(shù)的攝動(dòng)方法的分析效果。本文針對(duì)簡(jiǎn)單模型、水輪機(jī)簡(jiǎn)化模型進(jìn)行Monte-Carlo、一階攝動(dòng)、二階攝動(dòng)方法的對(duì)比分析，明確二階擾動(dòng)項(xiàng)對(duì)分析結(jié)果的影響，確定二階攝動(dòng)法的效果，最終為解決參數(shù)含大擾動(dòng)參數(shù)時(shí)，方法的選擇提供依據(jù)。

1 基于二階泰勒展開的振動(dòng)特性分析

固有頻率(模態(tài)頻率)的計(jì)算相對(duì)簡(jiǎn)單，也能夠反映系統(tǒng)的振動(dòng)特性，因此以推導(dǎo)固有頻率為核心，ω(a,t)，a為結(jié)構(gòu)的質(zhì)量、剛度、阻尼及響應(yīng)組成的參數(shù)向量。當(dāng)結(jié)構(gòu)參數(shù)具有擾動(dòng)性時(shí)，假定參數(shù)服從正態(tài)分布，針對(duì)固有頻率可進(jìn)行如下分析：

在常規(guī)結(jié)構(gòu)中，系統(tǒng)、材料參數(shù)的隨機(jī)性是客觀存在的，將質(zhì)量、剛度、阻尼描述為一個(gè)隨機(jī)向量an×1，n為參數(shù)個(gè)數(shù)，假定這些參數(shù)的概率特性是已知的。若隨機(jī)參數(shù)服從正態(tài)分布，其線性變換也服從正態(tài)分布，乘法一般也服從正態(tài)分布，因此，若an×1服從正態(tài)分布，則一般a的函數(shù)，F(xiàn)服從正態(tài)分布，即固有頻率ω也服從正態(tài)分布，將隨機(jī)變量分解為確定部分和擾動(dòng)部分，即

a=ad+ap

(1)

(2)

(3)

即

(4)

類推可得

(5)

采用泰勒展開，忽略二階以上擾動(dòng)量，梯度函數(shù)gi可表示為

(6)

即

(7)

其數(shù)學(xué)期望為

(8)

(9)

可以看到的是，在計(jì)算期望時(shí)，一階泰勒項(xiàng)消除，但仍有二階項(xiàng)，即二階擾動(dòng)項(xiàng)其修正作用，忽略二階項(xiàng)將影響精度。

(10)

其中，定義

(11)

(12)

(13)

(14)

在計(jì)算中，通過常規(guī)結(jié)構(gòu)分析方法、有限元分析、矩陣分解可得到解析函數(shù)表達(dá)式，二階偏導(dǎo)的計(jì)算較為復(fù)雜，但可通過分解巨型方程的相關(guān)算法進(jìn)行計(jì)算。

上述方法給出了考慮隨機(jī)變量的傳遞率，以及傳遞率的期望、方差和傳遞系數(shù)的推導(dǎo)方法，在求解傳遞率的過程中只涉及了隨機(jī)參數(shù)的一階靈敏度，并且在整個(gè)運(yùn)算過程中只利用隨機(jī)變量的期望和方差，而隨機(jī)變量本身不參與運(yùn)算，大大簡(jiǎn)化了計(jì)算量，給解決工程實(shí)際問題帶來了方便。

當(dāng)采用攝動(dòng)法對(duì)于帶有不確定參數(shù)的結(jié)構(gòu)進(jìn)行動(dòng)力分析時(shí)，二階攝動(dòng)法能避免一階攝動(dòng)法中久期項(xiàng)問題，即在攝動(dòng)計(jì)算中引入的不應(yīng)存在的數(shù)值共振。文獻(xiàn)[10]研究表明當(dāng)結(jié)構(gòu)阻尼比與結(jié)構(gòu)固有頻率較大時(shí)，久期項(xiàng)效應(yīng)對(duì)一階(二階)攝動(dòng)響應(yīng)影響較小。本文所采用的二階攝動(dòng)展開方法適用于以上情況。當(dāng)結(jié)構(gòu)阻尼比較小，就會(huì)有久期項(xiàng)現(xiàn)象發(fā)生，此時(shí)可以用阻帶濾波的方式來抑制久期項(xiàng)影響[11]。

2 簡(jiǎn)單結(jié)構(gòu)驗(yàn)證計(jì)算

在驗(yàn)證方法正確性時(shí)，盡量采用簡(jiǎn)單結(jié)構(gòu)，不考慮阻尼的影響，因此，計(jì)算采用三結(jié)點(diǎn)兩單元的桿結(jié)構(gòu)進(jìn)行驗(yàn)證，如圖1所示。不考慮結(jié)構(gòu)的假定參數(shù)包含服從正態(tài)分布的擾動(dòng)性，各單元參數(shù)、隨機(jī)量均相同，加性噪聲a的均值為0，加性噪聲方差用N(ad，σ2)表示：?jiǎn)卧獎(jiǎng)偠菿為N(108N/m，1×1014); 單元質(zhì)量M為N(104kg，2.5×105)。

圖1 三結(jié)點(diǎn)兩自由度桿單元模型示意圖Fig.1 The schematic diagram about pole element model with three nodes and two free units

表1給出了將考慮參數(shù)擾動(dòng)性的結(jié)構(gòu)，按照Monte-Calo模擬，一階攝動(dòng)，二階攝動(dòng)三種形式進(jìn)行求解。由計(jì)算機(jī)生成參數(shù)的樣本空間，每個(gè)樣本空間含10 000個(gè)樣本，但由于計(jì)算機(jī)生成均值不絕對(duì)為設(shè)定值，其結(jié)果與解析解存在一定的差異，因此，算例中以生成樣本的均值為依據(jù)，計(jì)算所得的固有頻率并不嚴(yán)格等于解析解2.93×103(rad/s)2和1.71×104(rad/s)2；一階攝動(dòng)方法是忽略二階擾動(dòng)項(xiàng)所計(jì)算的結(jié)果，關(guān)于一階攝動(dòng)的計(jì)算，已發(fā)表文獻(xiàn)[12-13]有詳細(xì)的推導(dǎo)過程及驗(yàn)證方法；二階攝動(dòng)方法是根據(jù)式(9)、式(10)計(jì)算的結(jié)果。表1中Error(%)是一階攝動(dòng)、二階攝動(dòng)方法所計(jì)算的均值、方差與Monte-Calo模擬之間的誤差，計(jì)算時(shí)間由于工作環(huán)境不同而不同，但在同一環(huán)境下各種方法之間計(jì)算耗時(shí)可以一定程度上的反映方法的計(jì)算效率。

表1 兩種形式求解各階固有頻率的計(jì)算值

表1的結(jié)果中：①以一階攝動(dòng)的方法所計(jì)算的均值誤差在0.251%大于二階攝動(dòng)計(jì)算所產(chǎn)生的誤差0.245%的誤差，降低幅度為2.39%；②方差項(xiàng)的誤差由0.959%將為0.309%，降低幅度67.8%；③二階攝動(dòng)方法所計(jì)算的時(shí)間略大于一階方法，但幅度不大。

以上結(jié)果表明，二階攝動(dòng)方法在計(jì)算簡(jiǎn)單結(jié)構(gòu)時(shí)，的確產(chǎn)生了優(yōu)化作用，其中方差優(yōu)化更為明顯，顯然計(jì)算方法是正確。

3 水電機(jī)組簡(jiǎn)化模型的振動(dòng)特性分析

以水輪機(jī)豎向振動(dòng)傳遞至廠房結(jié)構(gòu)的路徑中的兩條：①轉(zhuǎn)輪～軸系～軸承～固定部件(機(jī)架、頂蓋)～廠房；②轉(zhuǎn)輪～轉(zhuǎn)輪負(fù)壓區(qū)～頂蓋～廠房為分析對(duì)象，建立簡(jiǎn)化模型，利用軸系-廠房振動(dòng)微分方程來分析其傳導(dǎo)率：假定系統(tǒng)為線性系統(tǒng)，應(yīng)用拉格朗日方程建立振動(dòng)微分方程

(15)

將軸系、轉(zhuǎn)子、下機(jī)架、頂蓋、機(jī)墩平衡方程歸并到一起，擴(kuò)展后得到總剛度陣為

(16)

采用集中質(zhì)量陣，最后形成總體質(zhì)量陣M=diag{m1,m2,m3,m4,m5,m6,m7,m8}，假設(shè)整個(gè)機(jī)組廠房結(jié)構(gòu)只在水輪機(jī)轉(zhuǎn)輪處受豎向簡(jiǎn)諧激勵(lì)，且初始相位在整個(gè)過程中不發(fā)生改變，可得

(-ω2M+iωC+K)U=F(t)

(17)

通過解析計(jì)算得到結(jié)構(gòu)各點(diǎn)的響應(yīng)向量U，從而根據(jù)各路徑經(jīng)過的結(jié)構(gòu)及其附在結(jié)構(gòu)，可計(jì)算水輪機(jī)通過軸系統(tǒng)和頂蓋系統(tǒng)傳遞至機(jī)墩的力。定義傳遞率為傳遞力的幅值與振源的幅值之比：βt=|ft/f0|，其中，f0為激振力，ft為傳遞力，根據(jù)隨機(jī)變量概率分析的代數(shù)綜合法(Algebra Synthesis Method)，傳遞率的期望和方差為以及傳遞系數(shù)

(18)

以某巨型水電站傘式機(jī)組振動(dòng)傳遞路徑，如圖2所示。不考慮蝸殼及其下部結(jié)構(gòu)影響，假設(shè)激勵(lì)為簡(jiǎn)諧激勵(lì)，隨機(jī)參數(shù)均值由水電站施工設(shè)計(jì)圖計(jì)算得到，各參數(shù)均值如下：m1=8.28×104，m2=1.04×106，m3=3.29×105，m4=9.00×105，m5=1.20×105，m6=1.15×105，m7=1.39×105，m8=8.92×105質(zhì)量m的單位為kg；k1=7.26×1010，k3=5.72×1010，k4=2.32×1010，k51=2.20×1012，k52=9.41×109，k6=7.70×109，k7=4.26×108，k81=1.73×108，k82=1.73×1010，剛度k單位N/m；c1=5.48×106，c3=4.11×106，c4=1.02×107，c51=2.57×107，c52=7.51×104，c6=1.64×106，c7=9.74×105，c81=2.23×105，c82=9.99×104，阻尼c的單位為N·s/m。

這些參數(shù)中，油膜、水封剛度、阻尼非線性強(qiáng)，難以獲取化，本文根據(jù)其他電站的實(shí)測(cè)數(shù)據(jù)，利用相似原理獲取參數(shù)的均值，取值誤差較大，因此，取油膜和水封的剛度、阻尼值k51，k81，c51，c81的方差為0.15倍均值；頂蓋上承載的控制部件等附加部件較多，致使剛度、阻尼k82，c82不確定較強(qiáng)，其均值由計(jì)算采用折減系數(shù)得到，取其方差為0.1倍均值；其余參數(shù)方差為0.05倍均值。

圖3～圖6給出了Monte-Calo模擬，一階攝動(dòng)，二階攝動(dòng)的計(jì)算結(jié)果。其中，各方法與簡(jiǎn)單算例中計(jì)算一致，此時(shí)方程具有26個(gè)隨機(jī)參數(shù)。

圖2 傘式機(jī)組簡(jiǎn)化模型Fig.2 Simplified model of the umbrella unit

圖3為期望的計(jì)算結(jié)果，圖4為期望的誤差分析圖，其中，一階攝動(dòng)方法與Monte-Calo模擬計(jì)算的結(jié)果非常接近，兩條曲線的線型變化趨勢(shì)非常接近，一階攝動(dòng)法計(jì)算具有一定的波動(dòng)，從圖4可知，整體幅度不大；二階攝動(dòng)計(jì)算結(jié)果波動(dòng)性變化較大，特別是在固有頻率處波動(dòng)更為顯著，表明在期望計(jì)算時(shí)，二階項(xiàng)對(duì)結(jié)果產(chǎn)生一定負(fù)面效果。

圖5為方差的計(jì)算結(jié)果，圖6為期望的誤差分析圖，其中二階攝動(dòng)的計(jì)算，整體對(duì)一階項(xiàng)造成了擴(kuò)大作用，特別是在固有頻率處，放大作用更為明顯，并沒有產(chǎn)生相應(yīng)的優(yōu)化作用，反而起到負(fù)面效果。

圖3 不同方法對(duì)路徑傳遞率期望的頻域特性曲線Fig.3 Path transmissibility of the excitation frequency curve

圖4 一階、二階方法計(jì)算傳遞率期望的誤差分析圖Fig.4 Error graph of transmissibility excitation by first order and two order methods

圖5 不同方法對(duì)路徑傳遞率方差的頻域特性曲線Fig.5 Variance of the path transmissibility of the excitation frequency characteristic curve

圖6 一階、二階方法計(jì)算傳遞率方差的誤差分析圖Fig. 6 Error graph of transmissibility variance by first order and two order methods

4 結(jié) 論

以往關(guān)于二階項(xiàng)的推導(dǎo)是利用力學(xué)方程等關(guān)系式，將其轉(zhuǎn)化為一階攝動(dòng)方法，并未涉及二階項(xiàng)，特別是二階偏導(dǎo)數(shù)的計(jì)算，通過本文對(duì)二階攝動(dòng)方法的推導(dǎo)，明確了高階項(xiàng)的計(jì)算方式，并利用簡(jiǎn)單模型、復(fù)雜模型驗(yàn)證了方法的正確性，得到了如下結(jié)論：

(1)二階攝動(dòng)法在計(jì)算簡(jiǎn)單結(jié)構(gòu)時(shí)，能夠?qū)σ浑A攝動(dòng)法計(jì)算中所遇到的長久項(xiàng)進(jìn)行優(yōu)化，特別是方差項(xiàng)的優(yōu)化更為顯著。

(2)在水電機(jī)組傳導(dǎo)分析時(shí)，二階攝動(dòng)方法存在誤差，特別是在固有頻率處誤差更為嚴(yán)重，因此在復(fù)雜結(jié)構(gòu)體系中，選用二階攝動(dòng)方法需慎重考慮。

通過本文的研究，能夠明確以二階偏導(dǎo)數(shù)為基礎(chǔ)的二階攝動(dòng)法在計(jì)算水電機(jī)組振動(dòng)傳導(dǎo)分析中，能夠分析含較大擾動(dòng)項(xiàng)，如軸承油膜、不平衡磁拉力等參數(shù)的影響，但需要注意的是，在結(jié)構(gòu)的固有頻率處，采用二階偏導(dǎo)數(shù)的二階攝動(dòng)方法不再適用，可采用轉(zhuǎn)化為一階攝動(dòng)的二階攝動(dòng)方法進(jìn)行計(jì)算。本文的結(jié)論為接下來研究復(fù)雜結(jié)構(gòu)振動(dòng)傳導(dǎo)問題提供基礎(chǔ)