水輪發電機組碰摩轉子系統模型的周期解及數值分析

刁洋洋

（山東華邦建設集團有限公司，山東 濰坊 262500）

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

1 水輪發電機組碰摩轉子系統模型及運動方程

1.1 轉子-軸承模型

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

圖1 水輪發電機轉子-軸承系統示意圖

1.2 建立運動微分方程

1.2.1 模型基本假設

將發電機轉子簡化成質量為m1的圓盤結構，轉子的支撐結構為可滑動的軸承，該軸承為均勻的對稱結構，其上、下端質量均為m2。轉子位于上、下導軸承的中點，以線性方式處理軸承的阻尼和剛度。轉子的振動效應較復雜，存在多個方向的分量，該模型僅考慮橫向振動，不考慮陀螺力矩和扭轉效應的影響。

1.2.2 轉子-軸承碰摩運動微分方程

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

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

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

H函數式如公式（2）所示。

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

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

2 碰摩轉子系統模型的周期解及數值分析

2.1 周期解的近似表達式

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

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

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

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

2.1.2 周期解近似表達式的推導過程

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

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

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

2.2 數值分析

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

2.2.1 周期1的近似解

2.2.1.1 設置對照求解方法

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

2.2.1.2 數值分析結果

數值分析結果具體如下。

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

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 求解法在相同條件下得到的數值解如圖2（a）所示，二者的數據高度吻合，說明求解結果具有良好的可信度。上述2 種方法求解的轉子X方向時域圖如圖2（b）所示，其中τ為周期數。Runge-Kutta 法于大約第260 個周期達到穩態解，所用時間約為17.42s。與Runge-Kutta 法相比，HB-AFT 法更早達到穩態，所用時間約為7.72s。可見，HBAFT 法不僅具有較高的精確度，還能有效縮短求解時間。

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

2.2.2 周期運動的穩定性分析

2.2.2.1 改變軸徑間隙的周期解

在周期1 的運動分析中，將軸徑間隙設置為0.0002m，為了進一步驗證周期運行是否具有穩定規律，將軸徑間隙擴大至0.0014m，其他模擬條件保持不變。同樣利用HBAFT 法和Runge-Kutta 法進行迭代計算。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 穩定性分析

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

表2 不同軸徑間隙對應的最大Floquet 乘子及周期運動穩定性結果

3 結語

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