姬聯(lián)濤，袁浩強(qiáng)，荊岫巖，王 璞，劉 昊，楊雍杰，李超順

（1.中國(guó)電力科學(xué)研究院有限公司，江蘇 南京 210003；2.華中科技大學(xué) 土木與水利工程學(xué)院，湖北 武漢 430074；3.國(guó)家電網(wǎng)有限公司，北京 100031；4.國(guó)網(wǎng)湖南省電力有限公司，湖南 長(zhǎng)沙 410004）

0 引言

水電作為一種可再生能源，由于其諸多天然優(yōu)勢(shì)，在我國(guó)使用極為廣泛，各流域大中小型水電站建設(shè)數(shù)量眾多。隨著國(guó)內(nèi)更多中高水頭水電站的建設(shè)與增容改造，機(jī)組運(yùn)行工況更加惡劣、穩(wěn)定運(yùn)行要求更高，在實(shí)際運(yùn)行過(guò)程中，當(dāng)機(jī)組轉(zhuǎn)速上升到一定值時(shí)，激振力頻率與軸系固有頻率相近，且作用方向與振型方向一致時(shí)，軸系會(huì)產(chǎn)生強(qiáng)烈振動(dòng)，該轉(zhuǎn)速被稱為臨界轉(zhuǎn)速［1］。若機(jī)組在此轉(zhuǎn)速下長(zhǎng)時(shí)間運(yùn)行，可能造成結(jié)構(gòu)損壞、噪聲加劇，嚴(yán)重危及廠房及工作人員的生命財(cái)產(chǎn)安全，因而對(duì)水電機(jī)組軸系振動(dòng)特性尤其是臨界轉(zhuǎn)速的研究愈發(fā)重要。

目前針對(duì)此類轉(zhuǎn)子系統(tǒng)振動(dòng)特性分析，國(guó)內(nèi)外學(xué)者已經(jīng)做了大量的研究：榮吉利等［2］基于有限元法自主編制程序ROTOR，對(duì)紅石和隔河巖水電站軸系的自振特性進(jìn)行計(jì)算與分析，并與試驗(yàn)結(jié)果進(jìn)行對(duì)比；HUA 等［3］利用Muszynska 非線性密封力，構(gòu)建轉(zhuǎn)子與密封的動(dòng)力模型，并通過(guò)分岔圖和龐加萊圖分析了密封如何干擾系統(tǒng)的穩(wěn)定性；楊曉明等［4］采用有限元法建立水電機(jī)組軸系橫向振動(dòng)力學(xué)模型，并在考慮軸承簡(jiǎn)單耦合關(guān)系下求解了固有頻率與臨界轉(zhuǎn)速；LUNENO 等［5］針對(duì)組合軸承對(duì)水電機(jī)組轉(zhuǎn)子不同方向的耦合作用，分析其改變機(jī)組固有頻率與振動(dòng)模式的機(jī)理；翟黎明等［6］利用SAMCEF ROTOR 程序建立了抽水蓄能發(fā)電機(jī)組轉(zhuǎn)子-軸承-電磁系統(tǒng)三維有限元模型并分析了不平衡磁拉力與導(dǎo)軸承剛度變化對(duì)臨界轉(zhuǎn)速的影響；馬晨原等［7］采用有限元分析軟件探究了發(fā)電機(jī)質(zhì)量偏心及導(dǎo)軸承剛度變化對(duì)機(jī)組振動(dòng)特性的影響，并與現(xiàn)場(chǎng)測(cè)量數(shù)據(jù)進(jìn)行對(duì)比；徐瑞紅等［8］利用ANSYS 軟件對(duì)某水電站軸系進(jìn)行模態(tài)分析與臨界轉(zhuǎn)速計(jì)算，得出其發(fā)生共振可能性很小的結(jié)論。

以上研究為水電機(jī)組軸系的動(dòng)態(tài)特性分析與結(jié)構(gòu)優(yōu)化提供了借鑒，但多數(shù)模型過(guò)于簡(jiǎn)化，各類軸承與密封裝置的耦合作用未充分考慮，不能夠真實(shí)地反映出機(jī)組的邊界條件與工作狀態(tài)。鑒于此，本文針對(duì)某水電站的混流式機(jī)組，考慮上導(dǎo)、水導(dǎo)和推力軸承的耦合關(guān)系，分別構(gòu)建無(wú)機(jī)械故障、轉(zhuǎn)子或轉(zhuǎn)輪質(zhì)量偏心、聯(lián)軸法蘭不對(duì)中的有限元分析模型，進(jìn)行轉(zhuǎn)子動(dòng)力學(xué)仿真求解臨界轉(zhuǎn)速，從而分析軸系的動(dòng)態(tài)性能。

1 轉(zhuǎn)子動(dòng)力學(xué)仿真原理

轉(zhuǎn)子動(dòng)力學(xué)分析主要包括穩(wěn)定性分析、臨界轉(zhuǎn)速計(jì)算和不平衡響應(yīng)三部分［9］，起初科研工作者們多數(shù)采用傳遞矩陣法，但隨著計(jì)算機(jī)運(yùn)算性能的提高，研究者們更加青睞于精度更高、直觀性更強(qiáng)的有限元法?；谟邢拊?jiǎng)恿W(xué)仿真推導(dǎo)出來(lái)的完整轉(zhuǎn)子動(dòng)力學(xué)微分方程為：

求解轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速，就是在給定轉(zhuǎn)子系統(tǒng)轉(zhuǎn)速Ω的前提下，求解系統(tǒng)渦動(dòng)轉(zhuǎn)速ω。當(dāng)所得渦動(dòng)轉(zhuǎn)速ω與給定轉(zhuǎn)速Ω相等時(shí)，便稱此轉(zhuǎn)速為臨界轉(zhuǎn)速，該渦動(dòng)轉(zhuǎn)速下的振型就是臨界轉(zhuǎn)速的振型。

2 物理模型

2.1 幾何參數(shù)

由水輪機(jī)與發(fā)電機(jī)圖冊(cè)可知該機(jī)組為二導(dǎo)懸式混流式水輪發(fā)電機(jī)組，即存在兩個(gè)導(dǎo)軸承與一個(gè)推力軸承，利用三維可視化建模軟件UG，采用傳統(tǒng)的自底向上建模方法，即先對(duì)機(jī)組的轉(zhuǎn)輪、水輪機(jī)主軸、發(fā)電機(jī)大軸、轉(zhuǎn)子機(jī)架、磁極、磁軛、推力頭等部件進(jìn)行精細(xì)化建模隨后裝配，最終得到軸系計(jì)算簡(jiǎn)圖與三維裝配模型如圖1所示。

圖1 軸系計(jì)算簡(jiǎn)圖與三維裝配模型Fig.1 Shafting calculation diagram and three-dimensional assembly model

2.2 材料參數(shù)

根據(jù)圖紙結(jié)合各部分的物理特性，軸系各部件材料參數(shù)統(tǒng)計(jì)如表1所示，在ANSYS中建立材料庫(kù)，并根據(jù)表1對(duì)機(jī)組各部分進(jìn)行材料賦值。

表1 軸系各部件材料參數(shù)統(tǒng)計(jì)Tab.1 Material parameter statistics of shafting components

3 有限元模型

3.1 網(wǎng)格劃分

為了節(jié)約計(jì)算資源，并保證仿真結(jié)果的精度，采用映射劃分、自由劃分與掃掠劃分相結(jié)合的方法，并在網(wǎng)格不連續(xù)處設(shè)置綁定接觸，單元類型為六面體20 節(jié)點(diǎn)的SOLID186 與四面體10 節(jié)點(diǎn)的SOLID187 單元，經(jīng)劃分的模型節(jié)點(diǎn)數(shù)為359 150，單元數(shù)為164 660，結(jié)果如圖2所示。

圖2 水輪發(fā)電機(jī)組網(wǎng)格劃分結(jié)果圖Fig.2 Grid division result of hydro-generator set

3.2 邊界條件

該二導(dǎo)懸式水輪發(fā)電機(jī)組的固有頻率很大程度上受到水導(dǎo)、上導(dǎo)、推力軸承等邊界條件的影響，故在保證高精度物理模型的基礎(chǔ)上合理簡(jiǎn)化邊界條件尤為重要。對(duì)于實(shí)際的水電機(jī)組，軸承座、軸瓦和機(jī)墩等都具有一定彈性，在使用彈簧阻尼器模擬各軸承時(shí)，其彈簧剛度可視為油膜剛度、軸承座與機(jī)墩剛度三者的串聯(lián)，由于本研究著重分析軸系的振動(dòng)特性，故假設(shè)軸承座與機(jī)墩為剛性，僅考慮油膜動(dòng)力特性參數(shù)。

3.2.1 導(dǎo)軸承

為了準(zhǔn)確模擬上導(dǎo)和水導(dǎo)軸承的油膜剛度與阻尼，本文采用Lund 提出的八參數(shù)模型［10］，其對(duì)應(yīng)ANSYS 的彈簧阻尼單元COMBI214，包括線性化的4 個(gè)剛度系數(shù)與4 個(gè)阻尼系數(shù)，COMBI214單元幾何結(jié)構(gòu)如圖3所示。

圖3 COMBI214單元幾何結(jié)構(gòu)Fig.3 The geometry of COMBI214 unit

將支承視為各向同性，忽略交叉剛度與阻尼的影響，其主剛度與主阻尼可通過(guò)油膜力求得，但是油膜力表達(dá)式無(wú)法直接由流體潤(rùn)滑的Reynolds 方程導(dǎo)出，故選擇無(wú)限短軸承理論［11］對(duì)軸承做近似處理，得到導(dǎo)軸承無(wú)量綱剛度系數(shù)和阻尼系數(shù)的解析表達(dá)式為：

式中：K為無(wú)量綱剛度系數(shù)；C為無(wú)量綱阻尼系數(shù)；B為導(dǎo)軸承長(zhǎng)度；D為導(dǎo)軸承直徑；ε為軸頸中心的偏心率，等于偏心距離與軸承半徑間隙的比值。

導(dǎo)軸承剛度系數(shù)與阻尼系數(shù)的有量綱表達(dá)式為：

式中：k為有量綱剛度系數(shù)；c為有量綱阻尼系數(shù)；μ為潤(rùn)滑油度；n為軸頸轉(zhuǎn)速；L為軸承軸向長(zhǎng)度；φ為軸承間隙比，等于軸承半徑間隙與軸頸半徑的比值。

將涉及到的相關(guān)參數(shù)代入以上公式，軸頸中心偏心率取0.4，求得額定轉(zhuǎn)速下各導(dǎo)軸承剛度系數(shù)與阻尼系數(shù)：上導(dǎo)軸承剛度系數(shù)為8.55×106N/mm，阻尼系數(shù)為1.39×105（N·s）/mm；水導(dǎo)軸承剛度系數(shù)為7.09×107N/mm，阻尼系數(shù)為1.37×106（N·s）/mm。

3.2.2 推力軸承

考慮到推力軸承對(duì)軸系不僅有軸向的支撐作用，當(dāng)鏡板隨軸的彎曲振動(dòng)而傾斜時(shí)，油膜力還會(huì)提供恢復(fù)力矩，相當(dāng)于一個(gè)扭轉(zhuǎn)彈簧。因而將推力軸承約束視為伸縮彈簧與扭轉(zhuǎn)彈簧的共同作用，其油膜剛度系數(shù)與阻尼系數(shù)可通過(guò)文獻(xiàn)［12，13］所提供的解析表達(dá)式并結(jié)合自身參數(shù)求得，軸向剛度系數(shù)取為8.50×108N/mm，阻尼為1.01×106（N·s）/mm；扭轉(zhuǎn)剛度系數(shù)取為3.18×109（N·s）/rad，阻尼為7.64×107（N·s）/mm。

最終得到軸系三維有限元模型如圖4所示。

圖4 軸系三維有限元模型Fig.4 Three dimensional finite element model of shafting

4 計(jì)算結(jié)果與分析

4.1 無(wú)機(jī)械故障的軸系振動(dòng)特性分析

所謂無(wú)機(jī)械故障是指理想狀態(tài)下水電機(jī)組各部件的加工與設(shè)計(jì)階段完全一致，不存在材料的缺損、發(fā)電機(jī)轉(zhuǎn)子與水輪機(jī)轉(zhuǎn)輪的偏心與變形、聯(lián)軸法蘭不對(duì)中、軸系整體不對(duì)中等故障，且設(shè)備的安裝精度達(dá)到100%。利用軸系有限元模型考慮陀螺效應(yīng)在八個(gè)轉(zhuǎn)速點(diǎn)下進(jìn)行復(fù)模態(tài)分析，采用QR Damped 法提取各階模態(tài)，水電機(jī)組軸系前八階振型如圖5所示。

由圖5 可知，前兩階振型表現(xiàn)為以兩個(gè)導(dǎo)軸承為支點(diǎn)的擺動(dòng)，轉(zhuǎn)輪下環(huán)處位移最大、推力頭次之；第三階為轉(zhuǎn)輪的扭轉(zhuǎn)變形，下環(huán)處位移最大且沿周向均勻分布；第四、五階振型的危險(xiǎn)點(diǎn)位于推力頭上端，其次是轉(zhuǎn)輪，其余部位幾乎不發(fā)生振動(dòng)；第六、七階主要表現(xiàn)為推力頭、轉(zhuǎn)子、轉(zhuǎn)輪的振擺，其中位移最大處仍位于推力頭上端與轉(zhuǎn)輪下沿；第八階振型下水輪機(jī)很安全，主要表現(xiàn)為發(fā)電機(jī)整體的擺動(dòng)，推力頭處位移最大。隨后利用復(fù)模態(tài)分析結(jié)果繪制該轉(zhuǎn)速范圍內(nèi)固有頻率與轉(zhuǎn)速相關(guān)激振頻率的關(guān)系曲線，即Campbell圖，結(jié)果如圖6所示。由于本研究主要針對(duì)軸系的徑向振動(dòng)特性，故排除扭轉(zhuǎn)振型對(duì)應(yīng)的模態(tài)分析結(jié)果，前四階固有頻率隨機(jī)組轉(zhuǎn)速的變化情況如表2所示。

表2 前四階固有頻率隨機(jī)組轉(zhuǎn)速的變化情況Tab.2 Variation of the first four order natural frequencies with unit speed

圖6 無(wú)機(jī)械故障下軸系Campbell圖Fig.6 Campbell diagram of shafting without mechanical failure

結(jié)合圖6 與表2 分析可知，隨著機(jī)組轉(zhuǎn)速增大，正向渦動(dòng)固有頻率上升，反向渦動(dòng)固有頻率下降，前四階臨界轉(zhuǎn)速分別為：74.35、105.64、124.55、169.38 rad/s，但第一、三階為反向渦動(dòng)，在實(shí)際運(yùn)行過(guò)程中該共振并不會(huì)被激發(fā)，因而其并不能稱作真正意義上的臨界轉(zhuǎn)速。對(duì)于水輪發(fā)電機(jī)組，一般要求第一階臨界轉(zhuǎn)速為飛逸轉(zhuǎn)速的1.4 倍以上［14］，由資料可知該機(jī)組的飛逸轉(zhuǎn)速為300 r/min，即31.42 rad/s，遠(yuǎn)小于正向渦動(dòng)的第一階臨界轉(zhuǎn)速，故滿足設(shè)計(jì)要求，機(jī)組在正常運(yùn)行時(shí)不會(huì)發(fā)生共振現(xiàn)象。

4.2 質(zhì)量偏心對(duì)軸系振動(dòng)特性的影響

由于水電機(jī)組在安裝階段的誤差或者長(zhǎng)期運(yùn)行過(guò)程中的磨損變形，在軸系的主要部位，即發(fā)電機(jī)轉(zhuǎn)子和轉(zhuǎn)輪往往會(huì)出現(xiàn)質(zhì)量偏心問(wèn)題，為探究其偏心率其對(duì)軸系振動(dòng)的影響，首先計(jì)算出理想工況下兩者的質(zhì)量、體積、質(zhì)心位置、轉(zhuǎn)動(dòng)慣量等物理參數(shù)，再以該發(fā)電機(jī)定轉(zhuǎn)子間的平均氣隙長(zhǎng)度，即27 mm 為基準(zhǔn)，由平行軸定理［15］計(jì)算出轉(zhuǎn)子與轉(zhuǎn)輪分別偏心20%、40%、60%后的參數(shù)變化，結(jié)果如表3所示。

表3 轉(zhuǎn)子與轉(zhuǎn)輪物理參數(shù)隨偏心率的變化情況Tab.3 Variation of rotor and runner physical parameters with eccentricity

在ANSYS 中采用mass21 單元分別模擬發(fā)電機(jī)轉(zhuǎn)子或水輪機(jī)轉(zhuǎn)輪偏心對(duì)軸系的影響，隨后以同樣方式繪制Campbell 圖并計(jì)算各階臨界轉(zhuǎn)速，得到臨界轉(zhuǎn)速隨轉(zhuǎn)子或轉(zhuǎn)輪偏心率的變化情況如表4所示。

表4 臨界轉(zhuǎn)速隨轉(zhuǎn)子或轉(zhuǎn)輪偏心率的變化情況Tab.4 The variation of critical speed with the eccentricity of the rotor or runner

由表4 可知，隨著發(fā)電機(jī)轉(zhuǎn)子或者水輪機(jī)轉(zhuǎn)輪偏心率的增大，前四階臨界轉(zhuǎn)速的變化微乎其微，變化幅度不超過(guò)0.01%，因而可以得出結(jié)論：至少在該偏心范圍內(nèi)，轉(zhuǎn)子與轉(zhuǎn)輪的質(zhì)量偏心不會(huì)對(duì)軸系臨界轉(zhuǎn)速產(chǎn)生影響，即臨界轉(zhuǎn)速與單純的質(zhì)量偏心無(wú)關(guān)。

4.3 聯(lián)軸法蘭不對(duì)中對(duì)軸系振動(dòng)特性的影響

由電站資料可知，水輪機(jī)軸與發(fā)電機(jī)軸通過(guò)法蘭結(jié)構(gòu)連接，聯(lián)軸螺栓為20 條鉸制孔銷釘螺栓，其規(guī)格為M115×6 mm，長(zhǎng)度為480 mm，在機(jī)組吊裝過(guò)程中不可避免地會(huì)造成一定的誤差，因而聯(lián)軸法蘭不對(duì)中現(xiàn)象普遍存在，為了探究偏心距離對(duì)臨界轉(zhuǎn)速的影響，本文分別建立了聯(lián)軸法蘭偏心10、50、100 mm的有限元模型并進(jìn)行轉(zhuǎn)子動(dòng)力學(xué)分析，得到臨界轉(zhuǎn)速隨聯(lián)軸法蘭偏心距離的變化情況如表5所示。

表5 臨界轉(zhuǎn)速隨聯(lián)軸法蘭偏心距離的變化情況Tab.5 Variation of critical speed with eccentric distance of coupling flange

由表5可知，隨著聯(lián)軸法蘭偏心距離的增大，低階臨界轉(zhuǎn)速呈減小態(tài)勢(shì)，高階臨界轉(zhuǎn)速略微增大，當(dāng)聯(lián)軸法蘭偏心距離為100 mm時(shí)，前三階臨界轉(zhuǎn)速減小，第四階臨界轉(zhuǎn)速增大，與對(duì)中情況相比，增減幅度最大不超過(guò)0.33%，且實(shí)際運(yùn)行過(guò)程中聯(lián)軸法蘭偏心距離遠(yuǎn)達(dá)不到該數(shù)值，因而可以得出結(jié)論：水電機(jī)組臨界轉(zhuǎn)速受到聯(lián)軸法蘭不對(duì)中的影響不大。

5 結(jié)論

基于限元法建立了水輪發(fā)電機(jī)組軸系高精度動(dòng)力學(xué)模型，考慮了上導(dǎo)、水導(dǎo)和推力軸承的耦合關(guān)系并計(jì)算出各向剛度與阻尼值，在無(wú)機(jī)械故障、轉(zhuǎn)子或轉(zhuǎn)輪質(zhì)量偏心、聯(lián)軸法蘭不對(duì)中3種工況下進(jìn)行了轉(zhuǎn)子動(dòng)力學(xué)仿真，得到各階模態(tài)固有頻率、振型以及臨界轉(zhuǎn)速。通過(guò)分析可知：該水電機(jī)組軸系在上述3 種工況下的臨界轉(zhuǎn)速均滿足要求，且受質(zhì)量偏心與聯(lián)軸法蘭不對(duì)中的影響很小，該研究對(duì)水輪發(fā)電機(jī)組結(jié)構(gòu)優(yōu)化設(shè)計(jì)具有一定的意義。