刁洋洋

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

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

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

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

在直角坐標(biāo)系o-xyz中，水輪發(fā)電機(jī)組的定子和轉(zhuǎn)子運(yùn)動(dòng)系統(tǒng)如圖1所示，其中橢圓實(shí)線部分為發(fā)電機(jī)定子，橢圓虛線部分表示發(fā)電機(jī)轉(zhuǎn)子，定子和轉(zhuǎn)子的幾何中心為o點(diǎn)，軸徑的初始中心為s點(diǎn)。當(dāng)系統(tǒng)處于靜止?fàn)顟B(tài)時(shí)，o點(diǎn)和s點(diǎn)重合；當(dāng)系統(tǒng)處于運(yùn)動(dòng)狀態(tài)時(shí)，受振動(dòng)效應(yīng)的影響，o點(diǎn)和s點(diǎn)有可能出現(xiàn)位置偏移。A點(diǎn)和B點(diǎn)分別為上導(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 建立運(yùn)動(dòng)微分方程

1.2.1 模型基本假設(shè)

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

1.2.2 轉(zhuǎn)子-軸承碰摩運(yùn)動(dòng)微分方程

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

表1 庫(kù)侖摩擦力的基本觀點(diǎn)

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

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

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

式中：m1為電機(jī)轉(zhuǎn)子圓盤(pá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'分別為對(duì)坐標(biāo)的二階微分、一階微分；Y''、Y'、Z''、Z'、W''、W'的含義與X''、X'類(lèi)似。

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

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

2.1.1 運(yùn)動(dòng)微分方程的解析方法

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

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

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

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

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

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

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

2.2 數(shù)值分析

以某型水輪發(fā)電機(jī)為數(shù)值分析對(duì)象，其對(duì)應(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)子和軸承的對(duì)應(yīng)長(zhǎng)度記為L(zhǎng)r和Lb，取值分別是0.15m、0.3m。定子和轉(zhuǎn)子的平均氣隙長(zhǎng)度為δ0=0.0045m，轉(zhuǎn)子的質(zhì)量偏心e0=0.0006m。勵(lì)磁電流Ij=4A，碰摩擦系數(shù)f取0.01。將空氣的磁導(dǎo)系數(shù)記為μ0，取值為4π×10-7H/m，μ為潤(rùn)滑油絕對(duì)黏度，取值為1.8×103Pa·s。

2.2.1 周期1的近似解

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

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

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

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

首先，周期解的求解結(jié)果。水輪機(jī)大軸轉(zhuǎn)速ω取值為13rad/s，軸徑間隙設(shè)置為0.0002m，時(shí)域離散點(diǎn)的數(shù)量設(shè)置為144 個(gè)，迭代精度設(shè)定為1×10-11，諧波次數(shù)設(shè)置為8。根據(jù)隱式諧波平衡法的實(shí)施原理，通過(guò)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）。

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

圖2 HB-AFT 求解法與Runge-Kutta 求解法性能對(duì)比

2.2.2 周期運(yùn)動(dòng)的穩(wěn)定性分析

2.2.2.1 改變軸徑間隙的周期解

在周期1 的運(yùn)動(dòng)分析中，將軸徑間隙設(shè)置為0.0002m，為了進(jìn)一步驗(yàn)證周期運(yùn)行是否具有穩(wěn)定規(guī)律，將軸徑間隙擴(kuò)大至0.0014m，其他模擬條件保持不變。同樣利用HBAFT 法和Runge-Kutta 法進(jìn)行迭代計(jì)算。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 理論可用于分析周期運(yùn)動(dòng)的穩(wěn)定性，通過(guò)特征乘子判斷運(yùn)動(dòng)過(guò)程是否穩(wěn)定。如果最大特征乘子小于1，表明周期運(yùn)動(dòng)是穩(wěn)定的；當(dāng)最大特征乘子大于1 時(shí)，表明周期運(yùn)動(dòng)不穩(wěn)定（出現(xiàn)分叉）；如果最大特征乘子等于1，說(shuō)明周期運(yùn)動(dòng)處于臨界穩(wěn)定狀態(tài)[6]。設(shè)置不同的軸徑間隙，計(jì)算最大Floquet 乘子，結(jié)果見(jiàn)表2。從表2 可知，當(dāng)可調(diào)參數(shù)軸徑間隙為0.0001m、0.0002m、0.0009m 和0.0014m時(shí)，系統(tǒng)呈現(xiàn)出穩(wěn)定的周期運(yùn)動(dòng)，其他軸徑間隙均出現(xiàn)了分叉。

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

3 結(jié)語(yǔ)

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