李兆軍，張光政，劉福秀，丁 江

（廣西大學(xué)機(jī)械工程學(xué)院，廣西南寧 530004）

0 引言

水輪發(fā)電機(jī)組在運(yùn)行過(guò)程中的工況復(fù)雜多變，致使水流常常處于非穩(wěn)定狀態(tài)，從而使轉(zhuǎn)輪葉片在復(fù)雜的水流激勵(lì)作用下產(chǎn)生劇烈振動(dòng)，影響水輪機(jī)的安全穩(wěn)定運(yùn)行［1］，因而實(shí)現(xiàn)轉(zhuǎn)輪葉片振動(dòng)特性的有效調(diào)控使其振動(dòng)特性滿足水輪機(jī)正常工作是十分必要的。目前，學(xué)者們對(duì)轉(zhuǎn)輪葉片振動(dòng)特性控制方面的研究主要是集中在如何避免轉(zhuǎn)輪葉片產(chǎn)生水力共振等方面。例如，史廣泰等［2］采用ANSYS 平臺(tái)對(duì)水輪機(jī)轉(zhuǎn)輪進(jìn)行了模態(tài)計(jì)算和振動(dòng)分析，分析了機(jī)組產(chǎn)生共振的運(yùn)行條件；李海亮等［3］采用有限元法研究了混流式水輪機(jī)轉(zhuǎn)輪振動(dòng)特性以及各部件對(duì)轉(zhuǎn)輪固有頻率的影響規(guī)律；張新等［4］為了使轉(zhuǎn)輪的固有頻率避開(kāi)外界的激振頻率，利用CFD 軟件對(duì)水輪機(jī)轉(zhuǎn)輪進(jìn)行了有限元模態(tài)分析。但這些研究均未討論轉(zhuǎn)輪葉片振動(dòng)特性對(duì)工況參數(shù)的敏感程度。

動(dòng)態(tài)響應(yīng)靈敏度能夠反映設(shè)計(jì)變量或者參數(shù)對(duì)結(jié)構(gòu)振動(dòng)指標(biāo)的影響，許多學(xué)者開(kāi)展了相關(guān)方面的研究，例如，孫維光等［5］基于Sobol’法推導(dǎo)了動(dòng)力包機(jī)組振動(dòng)對(duì)懸掛隔振參數(shù)的全局靈敏度；王長(zhǎng)科等［6］對(duì)列車(chē)車(chē)體的扭轉(zhuǎn)振動(dòng)進(jìn)行了靈敏度分析，得到了對(duì)車(chē)體扭轉(zhuǎn)振動(dòng)頻率影響較大的設(shè)計(jì)參數(shù)；劉柏希等［7］利用偏導(dǎo)數(shù)法分析了RV 減速器系統(tǒng)固有頻率對(duì)轉(zhuǎn)動(dòng)慣量和剛度的影響。因而，積極開(kāi)展水輪機(jī)轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)靈敏度方面的研究對(duì)于有效控制轉(zhuǎn)輪葉片振動(dòng)特性具有十分重要的意義。然而，目前關(guān)于轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)靈敏度問(wèn)題的研究鮮見(jiàn)文獻(xiàn)報(bào)道。

本文以水輪機(jī)轉(zhuǎn)輪葉片為研究對(duì)象，根據(jù)轉(zhuǎn)輪葉片的動(dòng)力學(xué)模型，得到轉(zhuǎn)輪葉片的動(dòng)態(tài)響應(yīng)表達(dá)式，在此基礎(chǔ)上利用Sobol’法推導(dǎo)轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)工況參數(shù)的靈敏度，揭示工況參數(shù)對(duì)轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)的影響，并通過(guò)實(shí)例進(jìn)行分析。

1 轉(zhuǎn)輪葉片動(dòng)力學(xué)模型

考慮到轉(zhuǎn)輪葉片的幾何形狀以及在工作過(guò)程中的受力情況，其往往會(huì)發(fā)生彈性變形，故將轉(zhuǎn)輪葉片視為彈性體，采用矩形平面殼單元進(jìn)行轉(zhuǎn)輪葉片的模擬。利用有限元法，建立轉(zhuǎn)輪葉片的動(dòng)力學(xué)方程為［8］：

式中：M為轉(zhuǎn)輪葉片整體質(zhì)量矩陣，且M=Mt+Mp；Mt為單元的質(zhì)量矩陣；Mp為附加質(zhì)量矩陣，與水力參數(shù)有關(guān)，因而質(zhì)量矩陣M不僅與轉(zhuǎn)輪葉片的結(jié)構(gòu)參數(shù)有關(guān)，而且還和水力參數(shù)有關(guān)，為流固耦合項(xiàng)；K為轉(zhuǎn)輪葉片整體剛度矩陣；C為轉(zhuǎn)輪葉片阻尼矩陣；、、U分別為系統(tǒng)廣義位移向量、廣義速度向量、廣義加速度向量；Mtr為轉(zhuǎn)輪葉片的自激慣性力項(xiàng)，它不僅與轉(zhuǎn)輪葉片的運(yùn)動(dòng)參數(shù)有關(guān)，而且還與轉(zhuǎn)輪葉片的結(jié)構(gòu)參數(shù)和水力參數(shù)有關(guān)，為流固耦合項(xiàng)；F1為卡門(mén)渦引起的壓力脈動(dòng)，可由文獻(xiàn)［9］確定，F(xiàn)2為尾水管渦引起的壓力脈動(dòng)，可由文獻(xiàn)［10］確定，F(xiàn)3為動(dòng)靜干涉引起的壓力脈動(dòng)，可由文獻(xiàn)［11］確定。

根據(jù)動(dòng)力學(xué)方程，轉(zhuǎn)輪葉片頻率方程可表示為：

根據(jù)式（2），即可求出轉(zhuǎn)輪葉片的固有頻率。進(jìn)一步，利用振型疊加法求解出轉(zhuǎn)輪葉片的振動(dòng)響應(yīng)：

式中：N為轉(zhuǎn)輪葉片自由度數(shù)；Φ為正則模態(tài)矩陣；η為相應(yīng)的模態(tài)坐標(biāo)列陣。

2 動(dòng)態(tài)響應(yīng)靈敏度分析

靈敏度分析可以反映輸入變量對(duì)輸出變量的影響程度，從數(shù)學(xué)上可表示為：若某一函數(shù)f（x）可導(dǎo)，則其一階靈敏度可表示為［12］：

式中：xi為所分析參數(shù)變量的集合；S又稱f（x）對(duì)xi的一階微分靈敏度。

這里采用Sobol’法來(lái)推導(dǎo)轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)工況參數(shù)的靈敏度，Sobol’法是一種基于方差的蒙特卡洛法，其大致過(guò)程如下所示［13］：

首先，將轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)表示為U=f(X)，其中X=(x1，x2，…，xk)為影響轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)的k個(gè)工況參數(shù)所組成的參數(shù)向量，進(jìn)而可建立輸入?yún)?shù)的空間域?yàn)镻k。

其次，將f（X）表示為2k個(gè)遞增項(xiàng)之和，即：

式中：f0為常數(shù)。其余子項(xiàng)對(duì)其所包含的變量的積分一定為零，即：

根據(jù)式（5）和（6），可以看出各項(xiàng)之間都是正交的，如果(i1，i2，…，is) ≠(j1，j2，…，jn)，有：

進(jìn)一步有：

對(duì)式（5）取多重積分，即可得到各子項(xiàng)的表達(dá)式：

式中：X-i表示xi之外的其他工況參數(shù)，同理X-(ij)表示xi和xj之外的其他工況參數(shù)。

基于上述推導(dǎo)，可得到函數(shù)的總方差和偏方差的表達(dá)式為：

最后，各工況參數(shù)變量對(duì)轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)的靈敏度可表示為：

式中：Di1，i2，…，is和D分別為f（x）的偏方差和總方差；Si1為轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)工況參數(shù)x1一階靈敏度，反映x1單獨(dú)變化對(duì)轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)的影響，同理可得其他轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)其他工況參數(shù)變量的一階靈敏度。

3 實(shí)例分析

以國(guó)內(nèi)某700 MW 大型混流式水輪機(jī)為例，主要計(jì)算參數(shù)如下：額定轉(zhuǎn)速n=107.14 r∕min；額定流量Q=554.52 m3∕s；額定負(fù)荷W=571 MW；轉(zhuǎn)輪最大外徑為D=8 051.8 mm；轉(zhuǎn)輪葉片數(shù)Z1=13，單個(gè)轉(zhuǎn)輪葉片的質(zhì)量mt=7.7 × 103kg，轉(zhuǎn)輪葉片的材料密度為ρ=7.85×103kg∕m3，轉(zhuǎn)輪葉片的彈性模量為E=210 GPa，泊松比μ=0.3；導(dǎo)葉數(shù)Z2=26；轉(zhuǎn)輪葉片出水邊厚度δ=0.026 m；轉(zhuǎn)輪葉片出水口流角β=20°；斯特魯哈數(shù)sh=0.26；轉(zhuǎn)輪葉片尾部形狀影響系數(shù)k=2.4。考慮到水輪機(jī)工況參數(shù)較多，這里選擇一個(gè)易于控制的工況參數(shù)（轉(zhuǎn)速）進(jìn)行分析，假設(shè)轉(zhuǎn)速的變化范圍為0.5n至1.2n，轉(zhuǎn)速每上升0.1n取一個(gè)轉(zhuǎn)速點(diǎn)進(jìn)行分析，共8 個(gè)轉(zhuǎn)速點(diǎn)，記為ni（i=1，2，…，8），為了簡(jiǎn)化計(jì)算，這里假設(shè)ni服從均勻分布，且取值范圍為［ni-0.02n，ni+0.02n］，如表1 所示。采用拉丁超立方采樣法對(duì)ni進(jìn)行采樣，采樣次數(shù)設(shè)置為2 000次，并基于Sobol’法計(jì)算轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)轉(zhuǎn)速的靈敏度，如圖1所示。

表1 轉(zhuǎn)速概率分布 r∕minTab.1 Probability distribution of speed

圖1 轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)轉(zhuǎn)速的靈敏度曲線Fig.1 Sensitivity curve of blade dynamic response to speed

從圖1 中可以以看出，在轉(zhuǎn)速?gòu)膎1變化至n8的過(guò)程中，轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)轉(zhuǎn)速的靈敏度呈先增大后減小再增大的趨勢(shì)，當(dāng)轉(zhuǎn)速變化至n6（0.6n）以及n8（1.2n）時(shí)，轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)轉(zhuǎn)速的靈敏度出現(xiàn)了兩個(gè)峰值，說(shuō)明在這兩個(gè)轉(zhuǎn)速附近轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)轉(zhuǎn)速的變化是較為敏感的。

從上述分析可以看出，轉(zhuǎn)速對(duì)轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)的影響是隨轉(zhuǎn)速的大小而不斷變化的，如果要想實(shí)現(xiàn)在水輪機(jī)運(yùn)行過(guò)程中對(duì)轉(zhuǎn)輪葉片振動(dòng)的控制，應(yīng)關(guān)注靈敏度較大的兩個(gè)位置，即轉(zhuǎn)速為n6和n8。

為了驗(yàn)證仿真結(jié)果的正確性，進(jìn)行了水輪機(jī)變轉(zhuǎn)速實(shí)驗(yàn)，首先進(jìn)行測(cè)點(diǎn)及裝置布置，即在轉(zhuǎn)輪葉片上布置應(yīng)變片，并將應(yīng)變信號(hào)采集儀安裝在主軸軸頂；然后讓水輪機(jī)分別在每個(gè)轉(zhuǎn)速下運(yùn)行并采集動(dòng)應(yīng)力信號(hào)；最后進(jìn)行數(shù)據(jù)處理并計(jì)算動(dòng)應(yīng)力通頻值，通過(guò)試驗(yàn)采集得到轉(zhuǎn)輪葉片出水邊某一測(cè)點(diǎn)的動(dòng)應(yīng)力峰峰值隨轉(zhuǎn)速的變化曲線，如圖2所示。

圖2 轉(zhuǎn)輪葉片動(dòng)應(yīng)力隨轉(zhuǎn)速的變化曲線Fig.2 Variation curve of dynamic stress with rotating speed

從圖2 中可以看出，隨著轉(zhuǎn)速∕額定轉(zhuǎn)速的增大，轉(zhuǎn)輪葉片的動(dòng)應(yīng)力響應(yīng)也隨之增大，當(dāng)轉(zhuǎn)速比在［0.8，1］這一區(qū)間內(nèi)，轉(zhuǎn)輪葉片動(dòng)應(yīng)力響應(yīng)上升較為緩慢，處于一個(gè)較為穩(wěn)定的狀態(tài)，總的來(lái)說(shuō)，轉(zhuǎn)輪葉片動(dòng)應(yīng)力響應(yīng)隨機(jī)組的轉(zhuǎn)速的增加而單調(diào)上升。

根據(jù)圖2 所示試驗(yàn)所測(cè)得的轉(zhuǎn)輪葉片動(dòng)應(yīng)力隨轉(zhuǎn)速的變化曲線，計(jì)算其每個(gè)轉(zhuǎn)速比處的斜率，并進(jìn)行歸一化處理，得到轉(zhuǎn)輪葉片動(dòng)應(yīng)力對(duì)轉(zhuǎn)速的靈敏度曲線，并將仿真得到曲線與其進(jìn)行對(duì)比，如圖3所示。

從圖3 中可以看出，通過(guò)實(shí)測(cè)數(shù)據(jù)繪制成的靈敏度曲線和仿真所得到的靈敏度曲線基本一致，其峰值都出現(xiàn)在轉(zhuǎn)速為n6和n8處，只是靈敏度的大小有所差異。此外，該靈敏度曲線也與圖2 中轉(zhuǎn)輪葉片動(dòng)應(yīng)力隨轉(zhuǎn)速的變化趨勢(shì)較為符合，故可以說(shuō)明本文通過(guò)仿真計(jì)算得到的結(jié)果是較為準(zhǔn)確的。根據(jù)靈敏度分析的結(jié)果，可得到對(duì)轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)變化最敏感的轉(zhuǎn)速范圍，從而為轉(zhuǎn)輪葉片的振動(dòng)控制提供依據(jù)。

圖3 轉(zhuǎn)輪葉片動(dòng)應(yīng)力對(duì)轉(zhuǎn)速的靈敏度曲線Fig.3 Sensitivity curve of blade dynamic stress to speed

4 結(jié)論

Sobol’法是一種基于方差的蒙特卡洛法，采用Sobol’法能夠有效地分析和計(jì)算轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)工況參數(shù)的靈敏度。首先，建立轉(zhuǎn)輪葉片的動(dòng)力學(xué)模型，推導(dǎo)轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)與各相關(guān)參數(shù)的內(nèi)在關(guān)系式，并構(gòu)建相關(guān)輸入?yún)?shù)的空間域；然后，將轉(zhuǎn)輪葉片的動(dòng)態(tài)響應(yīng)函數(shù)式表示為2k個(gè) 遞增項(xiàng)之和，并且使各項(xiàng)之間都是相互正交的；最后，基于轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)函數(shù)的總方差和偏方差表達(dá)式，得到轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)靈敏度計(jì)算公式。實(shí)例分析結(jié)果表明，轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)轉(zhuǎn)速的靈敏度呈先增大后減小再增大的趨勢(shì)，在轉(zhuǎn)速∕額定轉(zhuǎn)速為0.6 以及1.2 時(shí)，轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)轉(zhuǎn)速的靈敏度出現(xiàn)了兩個(gè)峰值，說(shuō)明在這兩個(gè)轉(zhuǎn)速附近轉(zhuǎn)輪葉片動(dòng)態(tài)響應(yīng)對(duì)轉(zhuǎn)速的變化是較為敏感的。本文的研究為進(jìn)一步深入分析轉(zhuǎn)輪葉片的振動(dòng)特性與工況參數(shù)之間的內(nèi)在關(guān)系提供一個(gè)可行的方法。