向文平，黃社華，莊克云

（1.國能大渡河瀑布溝水力發(fā)電總廠，四川漢源 625300；2.武漢大學(xué)水利水電學(xué)院，湖北武漢 430072）

0 引言

近年來，隨著制造業(yè)和能源工業(yè)的快速發(fā)展，水電站向著大容量的方向發(fā)展［1］，例如2015 年底，水電站裝機(jī)容量已經(jīng)達(dá)到320.03 GW。但是這些變化將增加水輪發(fā)電機(jī)組軸系統(tǒng)因水力、機(jī)械和電氣因素引起的不穩(wěn)定振動［2-12］，這不僅影響了機(jī)組正常運(yùn)行，有的還危及到機(jī)組的安全。而且水輪發(fā)電機(jī)組軸系統(tǒng)在運(yùn)行過程中，機(jī)械、水力和電氣因素相互作用。例如，當(dāng)水流作用于水輪發(fā)電機(jī)組旋轉(zhuǎn)部分時(shí)，會造成發(fā)電機(jī)轉(zhuǎn)子與定子間的氣隙不對稱，從而產(chǎn)生磁力，加劇振動。并且，當(dāng)旋轉(zhuǎn)部件運(yùn)動狀態(tài)變化時(shí)，也會影響水輪機(jī)流場和發(fā)電機(jī)磁場。因此，為保證水電站的安全穩(wěn)定運(yùn)行，對不穩(wěn)定因素下的水輪發(fā)電機(jī)組進(jìn)行非線性動力特性研究是十分必要的。

長期以來，為了研究機(jī)械和電氣不穩(wěn)定因素的影響，學(xué)者們提出了大量的非線性數(shù)學(xué)模型［13-18］。一方面，關(guān)于機(jī)械不穩(wěn)定因素，Yan等［13］將分?jǐn)?shù)階導(dǎo)數(shù)引入水輪發(fā)電機(jī)機(jī)組，研究了多故障下彎扭耦合發(fā)電機(jī)轉(zhuǎn)子軸系統(tǒng)的動力學(xué)行為。Huang等［14］建立了不對中和碰摩故障下，轉(zhuǎn)子系統(tǒng)的動力學(xué)模型，并采用數(shù)值積分法對系統(tǒng)的動力學(xué)行為進(jìn)行了分析。Ma等［15］系統(tǒng)地研究了具有平行不對中和角度不對中的懸架轉(zhuǎn)子系統(tǒng)在起降過程中的油膜失穩(wěn)規(guī)律。另一方面，關(guān)于電氣不穩(wěn)定因素，Perers等［16］報(bào)道了對靜態(tài)偏心量為20%的水輪發(fā)電機(jī)的不平衡磁拉力的飽和效應(yīng)的研究，并采用有限元法和一個(gè)簡單的解析模型并行計(jì)算了各種空載電壓和負(fù)載下的磁拉力?；贘effcott 轉(zhuǎn)子模型，Xiang等［17］分析了永磁同步電機(jī)轉(zhuǎn)子系統(tǒng)的剛度特性，研究了不平衡磁拉力影響下的非線性動態(tài)行為。Zarko等［18］利用有限元方法計(jì)算了凸極同步發(fā)電機(jī)轉(zhuǎn)子軸在有磁場存在下的偏心運(yùn)動所產(chǎn)生的不平衡磁力。

水力不平衡是引起水輪發(fā)電機(jī)組振動的一個(gè)復(fù)雜而不可預(yù)測的因素，其機(jī)理是流動沿圓周分布不均勻。可能的原因包括葉片流道不一致、葉片型線不一致、葉尖間隙不一致和葉片接力行程不一致［5］。特別是對于葉片流道不一致來說，主要是由于葉片出口安裝角、接力行程和導(dǎo)葉開度不一致造成的。目前為止，雖然對水力不平衡有了一些研究，但主要集中在水力不平衡故障的監(jiān)測與診斷方面［19-21］。這種方法不能很好地滿足水電機(jī)組的經(jīng)濟(jì)和安全需求。因此，為了研究水力不平衡對水電機(jī)組非線性動態(tài)特性的影響，本文以非對稱翼型NACA2415 為研究對象，在雷諾數(shù)為0.5×105到3.5×105的范圍內(nèi)，對NACA2415 翼型非恒定繞流進(jìn)行CFD 數(shù)值模擬，并根據(jù)儒可夫斯基定理，建立適用于非對稱翼型水力不平衡力的非線性動力學(xué)數(shù)學(xué)模型。進(jìn)而根據(jù)拉格朗日方程建立水機(jī)電耦合故障下水輪發(fā)電機(jī)組非線性動力學(xué)數(shù)學(xué)模型。然后通過MATLAB 對該模型進(jìn)行計(jì)算分析，研究不同雷諾數(shù)下水輪發(fā)電機(jī)組隨葉片出口安裝角偏差變化的動力學(xué)特征。

1 翼型NACA2415非恒定繞流

1.1 計(jì)算條件和網(wǎng)格劃分

模型建立在C 型的計(jì)算域中，計(jì)算域的左邊界距翼型尖端9c（c為翼型的特征長度），右邊界與翼型尾端相差14c，上下邊界到弦線各差10c。計(jì)算區(qū)域采用不同尺寸的四邊形結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分，如圖1所示。

圖1 翼型NACA2415計(jì)算網(wǎng)格區(qū)域Fig.1 Calculation grid region of the NACA2415 asymmetric airfoil blade

本文湍流模型采用Spalart-Allmaras 模型，邊界條件為為壓力遠(yuǎn)場邊界條件，離散方法為有限體積法，壓力和速度耦合采用SIMPLEC 算法；非定常計(jì)算時(shí)間步長取0.001 s，最大迭代為20次，迭代步數(shù)為2 000 步。

1.2 升阻力特性分析

升力是流體流經(jīng)翼型表面時(shí)翼型受到的垂直于流動方向的作用力，阻力則是翼型受到的與流向相反的作用力。為便于分析，基于無量綱的升力系數(shù)（CL），阻力系數(shù)（CD）進(jìn)行討論，其表達(dá)式為：

式中：L為翼型升力；D為翼型阻力；ρ為介質(zhì)密度；V∞為自由來流速度；c為特征尺度，取為翼型弦長。

圖2 給出了雷諾數(shù)為3.5×105時(shí)，升、阻力系數(shù)仿真數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)［22］的對比。

圖2 升、阻力系數(shù)仿真值與實(shí)驗(yàn)值［22］對比Fig.2 Comparison of lift and drag coefficients between the numerical calculation and experiment data［22］

由圖2可以看出，在研究攻角范圍內(nèi)，基于S-A湍流模型所得到的升阻力值隨攻角的變化曲線與實(shí)測值變化曲線基本一致，僅在臨界失速角附近誤差較大，說明本文采用S-A 湍流模型開展翼型非定常繞流研究是準(zhǔn)確可靠的。

基于此，通過MATLAB 進(jìn)行擬合，得到升力系數(shù)CL和升阻比λ近似計(jì)算公式如下：

改變馬赫數(shù)和來流速度，使雷諾數(shù)相應(yīng)變化為2.2×105和0.5×105，通過MATLAB 擬合得到的升力系數(shù)CL和阻力與升力比λ近似計(jì)算公式如下：

2 水力不平衡力非線性動力學(xué)數(shù)學(xué)模型

2.1 單個(gè)葉片受力分析

將轉(zhuǎn)輪內(nèi)水流絕對速度V沿葉片方向和轉(zhuǎn)輪圓周方向進(jìn)行分解，分別定義為相對速度W和圓周速度U。對于單個(gè)葉片，由儒可夫斯基定理得葉片的徑向力如下［21］：

式中：CL為轉(zhuǎn)輪葉片升力系數(shù)；F為水輪機(jī)葉柵中葉片最大投影面積；λ為升阻比；Wm為葉片前后相對流速矢量平均值；βm為平均相對速度的方向角；γ是轉(zhuǎn)輪葉片周圍液體的單位重量。

設(shè)葉片進(jìn)、出口水流相對速度分別為W1和W2，葉片進(jìn)、出口安裝角分別為β1和β2，將x軸設(shè)為圓周方向，y軸設(shè)為x軸垂直方向，則W1，W2和Wm可記為x軸和y軸方向速度分量的矢量之和：

假設(shè)葉片位置角為ζ，可得葉片徑向力在X和Y方向的分力為：

式中：ξ=ξ0+ωt，ξ0為所選葉片的初始位置角。

2.2 葉片水力不平衡力

在實(shí)際工程中，對于水輪機(jī)葉片，由于制造誤差、運(yùn)行中的變形和磨損，出口水流角總是存在差異，導(dǎo)致總徑向水推力不為零。這里假設(shè)轉(zhuǎn)輪所有葉片中只有一對葉片出口開口不均勻，即存在偏差。定義這一對葉片出口安裝角分別為β21，β22，且β22-β21=χ 為安裝角偏差。基于式（9）、（13）、（14）和（15），轉(zhuǎn)輪的徑向不對稱力即水力不平衡力如下：

式中：βm1，βm2分別為這對葉片的平均相對速度的方向角，且表達(dá)式如下：

其中，雷諾數(shù)為3.5×105，2.2×105和0.5×105時(shí)，升力系數(shù)CL和升阻比λ分別按式（3）～（8）計(jì)算。

2.3 水輪發(fā)電機(jī)組非線性動力學(xué)數(shù)學(xué)模型

發(fā)電機(jī)轉(zhuǎn)子（o1）與水輪機(jī)轉(zhuǎn)輪（o2）的坐標(biāo)關(guān)系［23］如圖3所示。

圖3 水輪發(fā)電機(jī)軸系統(tǒng)的坐標(biāo)關(guān)系Fig.3 Coordinate relationships of the HGU shaft system

根據(jù)圖3，發(fā)電機(jī)轉(zhuǎn)子（o1）與水輪機(jī)轉(zhuǎn)輪（o2）的坐標(biāo)關(guān)系如下：

式中：（xo2，yo2）為水輪機(jī)轉(zhuǎn)輪軸坐標(biāo)；θ0為水輪發(fā)電機(jī)組的初始角度，θ=ωt+θ0；Φ0為發(fā)電機(jī)轉(zhuǎn)子的初始角，Φ=ωt+Φ0；r是o1和o2之間的距離；e1和e2分別為發(fā)電機(jī)轉(zhuǎn)子和水輪機(jī)轉(zhuǎn)輪的質(zhì)量偏心量。

水輪發(fā)電機(jī)組的拉格朗日函數(shù)定義為動能和勢能之差，可表示為：

進(jìn)而可得水輪發(fā)電機(jī)組軸系統(tǒng)非線性動力學(xué)數(shù)學(xué)模型如下：

式中：Fx-ump，F(xiàn)y-ump為不平衡磁拉力［23］；Fx-oil，F(xiàn)y-oil為油磨力［23］；Fx-f，F(xiàn)y-f為阻尼力［23］；Fx-rub，F(xiàn)y-rub為碰磨力［25］。

3 計(jì)算與分析

采用龍格-庫塔法進(jìn)行數(shù)值模擬，分析了雷諾數(shù)為0.5×105，2.2×105和3.5×105三種情況下，葉片出口安裝角偏差（χ）對水輪發(fā)電機(jī)組非線性動態(tài)特性的影響。步長為0.01，每次模擬迭代步長為5 000 步，計(jì)算指定的初始值［23］為［xo1，vx，yo1，vy］=［0.000 1，0.000 1，0.000 1，0.000 1］。此外，根據(jù)參考文獻(xiàn)［23-26］，計(jì)算中涉及的參數(shù)取值如下：m1=1.5×104kg，m2=1.1×104

kg，c=6.5×104N·s/m，k1=8.5×107N/m，k2=6.5×107N/m，e1=0.000 5 m，e2=0.000 5 m，i=800 A，m3=1.0×103kg，ω=3.925 rad/s，Δd=0.000 106 8 m，l=5.0 m，Δl=0.000 2 m，μ0=4π × 10-7H/m，δ=0.008 m，f=0.012，Kc=3×107N/m，Q=42.86 m3/s，D1=2.0 m，b0=0.5 m，Φ=0 rad，θ0=1 rad，φ0=0.8 rad，α=0.136 rad。

3.1 模型動態(tài)特性比較

為了研究水力不平衡力對HGSS 振動性能的影響，并驗(yàn)證本文模型的合理性，將本文模型與文獻(xiàn)［23］模型的結(jié)果進(jìn)行了對比，如圖4所示。文獻(xiàn)［23］中沒有考慮水力不平衡力因素，本文模型的軸心坐標(biāo)（xo1和yo1）的振動幅值遠(yuǎn)遠(yuǎn)大于文獻(xiàn)［23］中模型的振動幅值，說明水力不平衡力會加劇HGSS 的振動。事實(shí)上，在實(shí)際水電站中已經(jīng)檢測到水力不平衡力的影響，其典型特征正是HGSS的振動幅值被放大。

圖4 軸心位移動態(tài)特性比較Fig.4 Comparison of dynamic characteristics for axis coordinates（xo1 and yo1）

3.2 葉片出口安裝角偏差（χ）對機(jī)組振動的影響

在雷諾數(shù)為0.5×105，2.2×105和3.5×105三種工況下，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）隨葉片出口安裝角偏差（χ）變化的分岔圖如圖5、6和圖7所示，其中圖5（a）、圖6（a）和圖7（a）為軸心橫坐標(biāo)（xo1），圖5（b）、圖6（b）和圖7（b）為軸心縱坐標(biāo)（yo1）。

圖6 軸心坐標(biāo)分岔圖（Re=2.2×105）Fig.6 Bifurcation charts of axis coordinates（Re=2.2×105）

由圖5、6和圖7可得，隨著葉片出口安裝角偏差（χ）從0 rad增大到0.091 rad，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）一開始保持幅值較小的兩周期振蕩，然后進(jìn)入短暫的且幅值較小的混沌1狀態(tài)，緊接著又從混沌狀態(tài)1進(jìn)入周期振蕩，且振動幅值越來越大，并且還出現(xiàn)了HOPF分岔現(xiàn)象，最后從分岔進(jìn)入永久的混沌狀態(tài)2，且振動幅值較大。

圖5 軸心坐標(biāo)分岔圖（Re=0.5×105）Fig.5 Bifurcation charts of axis coordinates（Re=0.5×105）

除此之外，比較圖5、圖6和圖7可以發(fā)現(xiàn)，Re=0.5×105時(shí)，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）在0.008 754 rd<χ<0.015 87 rad時(shí)進(jìn)入混沌1 狀態(tài)，在χ>0.084 25 rad 時(shí)進(jìn)入混沌2 狀態(tài)；Re=2.2×105時(shí)，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）在0.008 206 rd<χ<0.014 04 rad 時(shí)進(jìn)入混沌1 狀態(tài)，在χ>0.075 86 rad 時(shí)進(jìn)入混沌2狀態(tài)；Re=3.5×105，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）在0.006 383 rd<χ<0.011 31 rad 時(shí)進(jìn)入混沌1 狀態(tài)，在χ>0.055 07 rad 時(shí)進(jìn)入混沌2 狀態(tài)?；谝陨戏治?，隨著雷諾數(shù)的增大，一方面我們可以得到水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）的混沌1和混沌2 現(xiàn)象都前移了；另一方面，混沌1 的區(qū)間范圍變化不大，但是混沌2 的區(qū)間范圍明顯增大了，意味著隨著雷諾數(shù)增大，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）在葉片開口不均勻下更容易出現(xiàn)運(yùn)行不穩(wěn)定現(xiàn)象，也就是說隨著雷諾數(shù)的增大，很小的葉片出口不均勻的誤差就能對機(jī)組的穩(wěn)定運(yùn)行構(gòu)成威脅。

圖7 軸心坐標(biāo)標(biāo)分岔圖（Re=3.5×105）Fig.7 Bifurcation charts of axis coordinates（Re=3.5×105）

4 結(jié)論

本文首先提出一種采用SA 單方程紊流模型計(jì)算NACA2415 翼型升、阻力系數(shù)的方法，然后基于儒可夫斯基定理，建立了非對稱翼型葉片的水力不平衡力數(shù)學(xué)模型。結(jié)合不平衡磁拉力和碰磨力，建立了一個(gè)包含攻角和雷諾數(shù)的水輪發(fā)電機(jī)組水機(jī)電耦合故障非線性數(shù)學(xué)模型，并通過計(jì)算分析，研究了雷諾數(shù)和葉片出口安裝角偏差（χ）對水電機(jī)組振動特性的影響，且得到了以下結(jié)論：

（1）隨著葉片出口安裝角偏差（χ）從0 rad增大到0.091 rad，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）都經(jīng)歷兩周期振蕩-混沌1 狀態(tài)-兩周期振蕩-多周期振蕩-HOPF分岔-混沌2狀態(tài)。

（2）Re=0.5×105時(shí)，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）穩(wěn)定范圍為0 rad<χ<0.008 754 rad，0.015 87 rad<χ<0.084 25 rad；Re=2.2×105時(shí)，穩(wěn)定范圍為0 rad<χ<0.008 206 rad，0.014 04<χ<0.075 86 rad；Re=0.5×105時(shí)，穩(wěn)定范圍為0 rad<χ<0.006 383 rad，0.011 31 rad<χ<0.055 07 rad。

（3）隨著雷諾數(shù)的增大，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）的混沌現(xiàn)象前移，即在高雷諾數(shù)情況下，水輪發(fā)電機(jī)組軸心坐標(biāo)（xo1和yo1）在葉片開口不均勻下更容易出現(xiàn)運(yùn)行不穩(wěn)定現(xiàn)象。本文研究成果將有利于水輪機(jī)的優(yōu)化設(shè)計(jì)，并對水輪發(fā)電機(jī)組安全運(yùn)行提供參考依據(jù)。