水輪發(fā)電機(jī)組碰摩轉(zhuǎn)子系統(tǒng)模型的周期解及數(shù)值分析

刁洋洋

（山東華邦建設(shè)集團(tuán)有限公司，山東 濰坊 262500）

水輪發(fā)電機(jī)組的轉(zhuǎn)子和定子之間具有一定的軸徑間隙，受振蕩作用的影響，轉(zhuǎn)子在水平方向存在小幅度的偏移風(fēng)險，進(jìn)而導(dǎo)致定子和轉(zhuǎn)子間出現(xiàn)碰摩，探索碰摩轉(zhuǎn)子系統(tǒng)的周期解能夠建立科學(xué)的運動控制參數(shù)、降低碰摩對轉(zhuǎn)子和定子的破壞程度。由于運動過程屬于強(qiáng)非線性微分方程，難以直接求解，因此運用數(shù)值分析法進(jìn)行迭代計算。

1 水輪發(fā)電機(jī)組碰摩轉(zhuǎn)子系統(tǒng)模型及運動方程

1.1 轉(zhuǎn)子-軸承模型

在直角坐標(biāo)系o-xyz中，水輪發(fā)電機(jī)組的定子和轉(zhuǎn)子運動系統(tǒng)如圖1所示，其中橢圓實線部分為發(fā)電機(jī)定子，橢圓虛線部分表示發(fā)電機(jī)轉(zhuǎn)子，定子和轉(zhuǎn)子的幾何中心為o點，軸徑的初始中心為s點。當(dāng)系統(tǒng)處于靜止?fàn)顟B(tài)時，o點和s點重合；當(dāng)系統(tǒng)處于運動狀態(tài)時，受振動效應(yīng)的影響，o點和s點有可能出現(xiàn)位置偏移。A點和B點分別為上導(dǎo)軸承、下導(dǎo)軸承的幾何中心，ω為大軸轉(zhuǎn)速。將轉(zhuǎn)子的重心記為G，轉(zhuǎn)子的質(zhì)量偏心可為向量e0，大軸的旋轉(zhuǎn)偏心可為向量e。

圖1 水輪發(fā)電機(jī)轉(zhuǎn)子-軸承系統(tǒng)示意圖

1.2 建立運動微分方程

1.2.1 模型基本假設(shè)

將發(fā)電機(jī)轉(zhuǎn)子簡化成質(zhì)量為m1的圓盤結(jié)構(gòu)，轉(zhuǎn)子的支撐結(jié)構(gòu)為可滑動的軸承，該軸承為均勻的對稱結(jié)構(gòu)，其上、下端質(zhì)量均為m2。轉(zhuǎn)子位于上、下導(dǎo)軸承的中點，以線性方式處理軸承的阻尼和剛度。轉(zhuǎn)子的振動效應(yīng)較復(fù)雜，存在多個方向的分量，該模型僅考慮橫向振動，不考慮陀螺力矩和扭轉(zhuǎn)效應(yīng)的影響。

1.2.2 轉(zhuǎn)子-軸承碰摩運動微分方程

在僅考慮橫向振動的情況下，轉(zhuǎn)子在水平方向會產(chǎn)生一定的振動幅值，轉(zhuǎn)子和定子之間存在特定的距離，當(dāng)轉(zhuǎn)子的水平運動幅值超過該距離時，轉(zhuǎn)子和定子就會產(chǎn)生碰摩作用。從理論上講，碰摩過程必然會引發(fā)熱效應(yīng)，并且定子會發(fā)生不可預(yù)知的形變[1]。為了便于建立運動微分方程，可將定子形變視為線性形變，同時不考慮碰摩產(chǎn)生的熱效應(yīng)。綜合以上假設(shè)條件，定子和轉(zhuǎn)子之間的摩擦力符合庫侖摩擦力（見表1）的應(yīng)用條件。轉(zhuǎn)子處于運動狀態(tài)，當(dāng)產(chǎn)生碰摩作用時，以滑動摩擦力為主，將摩擦力分解至x軸和y軸，分解后的摩擦力如公式（1）所示。

表1 庫侖摩擦力的基本觀點

式中：Fx_rub、Fy_rub為撞摩擦力在x、y軸的分量；e1為發(fā)電機(jī)轉(zhuǎn)子的徑向位移；f為摩擦力系數(shù)；kr為定子徑向剛度；δ0為定子和轉(zhuǎn)子的平均氣隙長度；H為函數(shù)。

H函數(shù)式如公式（2）所示。

根據(jù)現(xiàn)有研究成果，水輪機(jī)轉(zhuǎn)子-軸承系統(tǒng)在運行過程中主要存在3 種作用力，分別為不平衡磁拉力（Unbalanced Magnetic Pull，UMP）、碰摩力和非線性油膜力。將軸承徑處的坐標(biāo)記為（Z，W），則定子-轉(zhuǎn)子系統(tǒng)的運動微分方程組如公式（3）所示。

式中：m1為電機(jī)轉(zhuǎn)子圓盤的質(zhì)量；m2為軸承上端或者下端所集中的質(zhì)量；轉(zhuǎn)子質(zhì)量偏心為e0；c1為轉(zhuǎn)子處的阻尼；c2為軸承處的阻尼；Fx_ump、Fy_ump是不平衡磁拉力在x軸、y軸方向的分量；Ke為大軸的剛度；fx、fy為非線性油膜力在x軸、y軸的分量；將轉(zhuǎn)子外圓幾何中心的坐標(biāo)記為（X，Y，0）；X''、X'分別為對坐標(biāo)的二階微分、一階微分；Y''、Y'、Z''、Z'、W''、W'的含義與X''、X'類似。

2 碰摩轉(zhuǎn)子系統(tǒng)模型的周期解及數(shù)值分析

2.1 周期解的近似表達(dá)式

2.1.1 運動微分方程的解析方法

諧波平衡法（Harmonic Balance，HB）是處理非線性問題的常用方法，但公式（3）屬于強(qiáng)非線性微分方程，HB 方法難以對該方程進(jìn)行解析[2]。隱式諧波平衡法（HB-AFT）由HB 方法和時域頻域轉(zhuǎn)換技術(shù)（Alternating Frequency/time，AFT）綜合發(fā)展而來，能夠解析強(qiáng)非線性微分方程，其原理如下。

假設(shè)存在某個非線性系統(tǒng)，其二階微分方程如公式（4）所示。

為了得到公式（4）的周期解，假設(shè)x與t的關(guān)系式為x（t），并且x（t）為周期性函數(shù)，其周期為T，于是有x（t）=x（t+T）。此時函數(shù)f（·）與x（t）建立了相同的正交基，可分別對其實施傅里葉展開，進(jìn)而獲得二者的隱式非線性代數(shù)關(guān)系式。

2.1.2 周期解近似表達(dá)式的推導(dǎo)過程

根據(jù)HB-AFT 方法的應(yīng)用原理，應(yīng)該先對公式（3）進(jìn)行非線性函數(shù)諧波平衡化，再從中獲取主方程，具體過程如下：對方程在X、Y、Z、W這4 個坐標(biāo)點的周期解實施傅里葉展開，得到相應(yīng)的傅里葉級數(shù)形式和線性外力的傅里葉級數(shù)[3]。例如非線性外力如公式（5）所示。式中：cx0、cy0、cz0以及cw0均為常數(shù)；cxk、cyk、czk以及cwk為余弦項相關(guān)的系數(shù)；dxk、dzk、dwk以及dyk為正弦項相關(guān)系數(shù)。

將經(jīng)過傅里葉展開的周期解和非線性外力代入公式（3）中，即可推導(dǎo)出有關(guān)常數(shù)項的代數(shù)方程，將該方程記為G。將X、Y、Z、W方向的周期解經(jīng)過傅里葉展開后同樣可得到常數(shù)項ax0、ay0、az0、aw0，相應(yīng)的余弦項系數(shù)為awk、azk、ayk、axk，相應(yīng)的正弦項系數(shù)包括bwk、bzk、byk以及bxk。

難以直接對隱式非線性方程G進(jìn)行求解，于是利用牛頓迭代法來求解，其實施步驟包括確定需要迭代的變量、構(gòu)造迭代關(guān)系式以及控制迭代過程[4]。利用X、Y、Z、W方向周期解傅里葉展開級數(shù)中的常數(shù)項、余弦系數(shù)項、正弦系數(shù)項，構(gòu)造出迭代關(guān)系式P，則有P=[ax0ay0az0aw0ax1ay1az1aw1bx1by1bz1bw1...axkaykazkawkbxkbykbzkbwk]。再將非線性外力經(jīng)過傅里葉展開所形成的常數(shù)項、正弦系數(shù)項、余弦系數(shù)項構(gòu)建立為迭代式Q，則有Q=[cx0cy0cz0c w0cx1cy1cz1cw1dx1dy1dz1dw1...cxkcykczkcwkdxkdykdzkdwk]。在求周期解的過程中，將迭代關(guān)系式Q作為已知因素，P為求解目標(biāo)，以不動點迭代法求解P，但計算過程涉及大量數(shù)學(xué)推導(dǎo)且較抽象，難以直觀展示，下文將通過數(shù)值分析展示求解過程。

2.2 數(shù)值分析

以某型水輪發(fā)電機(jī)為數(shù)值分析對象，其對應(yīng)的m1、m2質(zhì)量分別為60kg、25kg，轉(zhuǎn)子阻尼和軸承阻尼分別為c1=4000N·s/m、c2=1200N·s/m，轉(zhuǎn)子和軸承的半徑分別為Rr=0.06m、Rb=0.5m，將轉(zhuǎn)子和軸承的對應(yīng)長度記為Lr和Lb，取值分別是0.15m、0.3m。定子和轉(zhuǎn)子的平均氣隙長度為δ0=0.0045m，轉(zhuǎn)子的質(zhì)量偏心e0=0.0006m。勵磁電流Ij=4A，碰摩擦系數(shù)f取0.01。將空氣的磁導(dǎo)系數(shù)記為μ0，取值為4π×10-7H/m，μ為潤滑油絕對黏度，取值為1.8×103Pa·s。

2.2.1 周期1的近似解

2.2.1.1 設(shè)置對照求解方法

為了評價HB-AFT 方法在碰摩轉(zhuǎn)子系統(tǒng)周期解數(shù)值分析中的效果，將Runge-Kutta 方法作為對照組。Runge-Kutta方法是一種隱式或顯式迭代法，在非線性微分方程的求解中具有重要應(yīng)用，其優(yōu)點為求解精度高，缺點為實現(xiàn)原理較復(fù)雜，不易操作[5]。

2.2.1.2 數(shù)值分析結(jié)果

數(shù)值分析結(jié)果具體如下。

首先，周期解的求解結(jié)果。水輪機(jī)大軸轉(zhuǎn)速ω取值為13rad/s，軸徑間隙設(shè)置為0.0002m，時域離散點的數(shù)量設(shè)置為144 個，迭代精度設(shè)定為1×10-11，諧波次數(shù)設(shè)置為8。根據(jù)隱式諧波平衡法的實施原理，通過MATLAB 軟件求解出水輪機(jī)碰摩轉(zhuǎn)子系統(tǒng)周期1 的近似解析表達(dá)式，將X方向的周期解記為X（t），Y方向的周期解記為Y（t），結(jié)果如下。

X（t）=-1.6432231×10-11+6.13443091×10-3sin（ωt）-1.1020326×10-11cos（2ωt）+1.61331×10-11cos（3ωt）+1.12521×10-10sin（3ωt）；Y（t）=0+6.1343138×10-4cos（ωt）+1.10925×10-10cos（3ωt）。

其次，求解方法性能評價。HB-AFT 求解法和Runge-Kutta 求解法在相同條件下得到的數(shù)值解如圖2（a）所示，二者的數(shù)據(jù)高度吻合，說明求解結(jié)果具有良好的可信度。上述2 種方法求解的轉(zhuǎn)子X方向時域圖如圖2（b）所示，其中τ為周期數(shù)。Runge-Kutta 法于大約第260 個周期達(dá)到穩(wěn)態(tài)解，所用時間約為17.42s。與Runge-Kutta 法相比，HB-AFT 法更早達(dá)到穩(wěn)態(tài)，所用時間約為7.72s。可見，HBAFT 法不僅具有較高的精確度，還能有效縮短求解時間。

圖2 HB-AFT 求解法與Runge-Kutta 求解法性能對比

2.2.2 周期運動的穩(wěn)定性分析

2.2.2.1 改變軸徑間隙的周期解

在周期1 的運動分析中，將軸徑間隙設(shè)置為0.0002m，為了進(jìn)一步驗證周期運行是否具有穩(wěn)定規(guī)律，將軸徑間隙擴(kuò)大至0.0014m，其他模擬條件保持不變。同樣利用HBAFT 法和Runge-Kutta 法進(jìn)行迭代計算。HB-AFT 法在新條件下的X方向和Y方向周期解X1（t）和Y1（t）如下。

X1（t）=3.2755×10-8+2.82821×10-8cos（0.25ωt）-2.6488×10-1sin（0.25ωt）-2.4338×10-8cos（0.5ωt）-4.490×10-8sin（0.5ωt）-1.2406×10-8cos（0.75ωt）-1.6481×10-8sin（0.75ωt）-1.7878×10-7cos（ωt）。

Y1（t）=0+2.6144×10-1cos（0.25ωt）-2.8433×10-8sin（0.25ωt）+3.1124×10-8sin（0.5ωt）-1.2481×10-8cos（0.75ωt）+1.7214×10-8sin（0.75ωt）+1.7654×10-7sin（ωt）。

2.2.2.2 穩(wěn)定性分析

Floquet 理論可用于分析周期運動的穩(wěn)定性，通過特征乘子判斷運動過程是否穩(wěn)定。如果最大特征乘子小于1，表明周期運動是穩(wěn)定的；當(dāng)最大特征乘子大于1 時，表明周期運動不穩(wěn)定（出現(xiàn)分叉）；如果最大特征乘子等于1，說明周期運動處于臨界穩(wěn)定狀態(tài)[6]。設(shè)置不同的軸徑間隙，計算最大Floquet 乘子，結(jié)果見表2。從表2 可知，當(dāng)可調(diào)參數(shù)軸徑間隙為0.0001m、0.0002m、0.0009m 和0.0014m時，系統(tǒng)呈現(xiàn)出穩(wěn)定的周期運動，其他軸徑間隙均出現(xiàn)了分叉。

表2 不同軸徑間隙對應(yīng)的最大Floquet 乘子及周期運動穩(wěn)定性結(jié)果

3 結(jié)語

對以上研究內(nèi)容進(jìn)行總結(jié)，得到以下3 個結(jié)論：第一，通過HB-AFT 法對碰摩轉(zhuǎn)子系統(tǒng)的運動方程進(jìn)行迭代結(jié)算，可有效求出X方向和Y方向的周期解。第二，與Runge-Kutta 法相比，HB-AFT 法均能達(dá)到基本相同的計算精度，計算耗時也有明顯降低。第三，當(dāng)軸徑間隙為0.0001m～0.0014m 時，僅有0.0001m、0.0002m、0.0009m、0.0014m 共4 種取值能夠使碰摩轉(zhuǎn)子系統(tǒng)達(dá)到穩(wěn)定運行狀態(tài)。