汽輪機轉子系統中共振組織的全局拓撲規律分析

徐 璐， 褚衍東， 楊 瓊， 李險峰

(1.蘭州交通大學 機電工程學院，蘭州 730070； 2.蘭州交通大學 數理學院，蘭州 730070)

《世界能源展望》中提到可再生能源的全球份額在發電領域迅猛增長,現有常規發電廠的目標之一是補償剩余負荷，即使在將來的能源情景中，汽輪機依舊起著重要作用。因此，獲取更多汽輪機轉子系統動力學信息，對其機理進一步地認識，最終進行非線性動力學設計，具有重要的現實意義和廣泛的工程背景。

研究大型轉子系統的穩定性與分岔理論具有重要意義。描述系統動力學問題的微分方程具有非線性特性，在一定條件下必然發生分岔和混沌現象[1-5]。劉小峰等[6]利用分岔圖、Poincaré截面圖等，揭示了交叉剛度對碰摩轉子系統運行狀態的影響。徐文標等[7]基于Paris裂紋擴展規律，在參數平面內分析了裂紋擴展對含拉桿裂紋的轉子系統動力學特性的影響，發現了多種分岔現象。張金劍等[8]利用單參數數值工具探討了質量偏心、勵磁電流和氣隙動靜偏心對水電機組碰摩轉子-軸承系統的影響。我國對于汽流激振的研究起步較晚，是從1980年開始的。柴山等[9-11]對汽輪機的汽流激振力做了較詳細的研究，包括葉片優化設計、汽流激振力計算公式的推導及計算程序的設計等。接著，文獻[12-17]充分利用分岔圖、Poincaré截面圖、軸心軌跡圖、頻譜圖等單參數數值工具討論了高壓或低壓汽輪機轉子系統的分岔和混沌行為。最近，曹麗華等[18-19]以耦合了非線性汽流激振力和油膜力的超超臨界汽輪機轉子系統為研究對象，基于分岔圖、最大Lyapunov指數圖等，揭示了非線性汽流激振力對系統的影響。翁雷等[20]對汽輪機轉子系統做了較全面的分析，分別探討了葉尖汽流激振對汽輪機轉子、基礎松動轉子[21]、裂紋轉子[22]、碰摩轉子[23]及裂紋-碰摩[24]轉子系統的影響。羅躍剛等[25]分析了汽流激振和偏心量對密封-轉子系統的影響。

綜上所述，盡管關于轉子系統非線性動力學行為的研究有豐碩的成果，但是都是通過單參數分岔分析實現，導致動力學信息有限且不完整。這表明我們對大型旋轉機械的雙參數結構的理解還很遠，更不用說更高維的參數空間。因此，針對大型旋轉機械系統的非線性動力學行為的研究，真正的挑戰是發現并總結二維參數平面內分岔行為的拓撲普遍性。

基于以上原因，本文以油膜支撐，在汽流激振力、非線性油膜力和不平衡離心力共同作用下的汽輪機轉子系統為研究對象，利用高清的雙參數平面周期穩定相圖、分岔圖、Lyapunov指數圖和時間序列圖，分析了轉子系統內周期吸引子的全局拓撲規律，目的是補充和完善已得到的研究成果，更全面地掌握汽輪機轉子系統的非線性動力學特性，揭示更多隱藏在系統中的動力學信息，為大型旋轉機械實際設計及故障預測提供有價值的依據。

1 系統力學模型及運動方程

1.1 非線性間隙汽流激振力計算公式

圖1為系統受力簡圖，其中fax，fay是汽輪機非線性汽流激振力Fa在x和y方向的分力。圖2為沿環向的兩個橫剖面展開圖。

圖1 系統受力圖Fig.1 Loading diagram of the system

圖2 汽流在靜、動葉片間流動示意圖Fig.2 The schematic diagram of air-exciting force moving between stator and rotor blades

則非線性汽流激振力的無量綱模型為：

Fa=A1δE+A3δ3E3

(1)

其中

1.2 油膜力計算公式

忽略扭轉振動和陀螺力矩，在圖3中：O1,O2分別為軸承內瓦和轉子幾何中心；O3為轉子質心；轉子兩端由半徑為R，長度為L的滑動軸承支撐；m1,c1和m2,c2分別為轉子在軸承處和圓盤處的等效集中質量和結構阻尼；Fx,Fy為非線性油膜分量；h為平均油膜厚度。

圖3 油膜支撐轉子系統示意圖Fig.3 A schematic diagram of the rotor system supported by oil film

圖3是油膜支撐轉子示意圖，利用短軸承模型[26]，則設x，y為軸承位移，則無量綱油膜力在x和y方向的分力為

(2)

其中

(3)

(4)

(5)

(6)

1.3 系統動力學方程

設x1，y1為轉子左端軸承處的徑向位移；x2，y2為轉盤處的徑向位移，則系統的無量綱運動方程為

(7)

式中，τ=ωt。式(7)中的參數意義及取值，如表1所示。

表1 汽輪機轉子系統參數意義及取值Tab.1 The parameter meanings and values of the system

2 計算方法

在本文，雙參數平面內的數值結果主要使用isoperiodic圖[27]呈現。Isoperiodic圖不僅顯示了一個周期內尖峰數量的積累，即波形的復雜性，以及在參數平面內尖峰數量積累的地方以及積累的快慢，而且也可以有效地用于處理試驗數據。對于大型旋轉機械系統，制作isoperiodic圖是一項非常耗時的任務。因此，采用GPU并行計算技術。

為了生成isoperiodic穩定相圖，借助GPU并行計算的優勢，將感興趣的參數窗口劃分為N×N的等距點網格，通?？蛇_2 000×2 000=4×106。對于每個點，通過對方程進行數值積分來獲得時間演化。對方程組式(7)使用標準四階Runge-Kutta算法。為了進行積分，從任意選擇的初始條件開始，從左到右水平掃描參數，通過“跟隨吸引子”向右進行。在所有情況下，都會放棄起初2×105步以避免瞬態響應。在對時間步長進行后續80×105積分后，可以獲得isoperiodic穩定相圖。每一張isoperiodic圖會使用19種顏色的調色板顯示每個周期的尖峰數量。通過循環19個基本顏色“模19”來繪制超過19個尖峰的振蕩，即通過為它們分配一個顏色索引。黑色用于表示“混沌或擬周期”，即缺乏數字可檢測的周期性。

同時，本文采用一種新的方法，即根據運動方程初始小的擾動副本(“克隆”)的時間演化，計算模型的李雅普諾夫指數譜，其中雅可比矩陣的復雜計算替換為線性化版本的運動方程的時間演化?！翱寺》ā崩钛牌罩Z夫指數的簡要并行計算代碼請參考文獻[28-29]。

3 具體數值結果分析

在圖4中偏心量e和旋轉速度ω同時變化，進汽速度v保持200不變，剛度k=2.5×107。通過右側的顏色棒很容易區分周期窗口的周期數、大小及形狀。將圖4(a)的周期序列劃分為兩組，這兩組周期結構在圖4(b)和圖4(c)中進一步放大詳述。這些周期區域具有自相似性。

圖4 汽輪機轉子系統在(ω,e)參數平面內的穩定周期相圖Fig.4 The isoperiodic diagrams of the steam turbine rotor system in (ω,e) parameter plane

在圖4(b)中，13條周期數不同的“細絲”嵌入在大片的黑色區域內。為了更加深入地理解這組周期序列的演化細節，沿著藍色直線做單參數分岔圖。在圖5(a)中，水平掃描e，隨著e的增大且旋轉速度ω固定在2 000，13個周期窗口依次鎖定在非周期區域內：P-3→P-17→P-14→P-11→P-8→P-13→P-18→P-5→P-17→P-12→P-7→P-16→P-9。圖5(b)呈現了通過“克隆”法求得的最大Lyapunov指數圖，可見堆積在其中的周期區域將擬周期區域瓜分，這是一個周期解與擬周期解頻繁轉換的過程。擬周期解隨參數變化在某些參數域鎖相到周期解，因此這些周期區域稱為共振組織或鎖模結構。

圖5(a)中A點處(e=0.060 8×10-3)的狀態傳遞矩陣的特征值為：λ=[-0.760 7±0.131 8i,0.359 6±0.658 1i, 1.012 1, 0.356 8]，主導Floquet乘子|λ|max=1.012 1>1，即|λ|max通過正實數穿過單位圓，因此在A點處發生了鞍結分岔(SN)。在B點處(e=0.070 5×10-3)的狀態傳遞矩陣的特征值為：λ=[-0.542 7±0.086 0i, 0.446 1±0.702 0i, 1.052 7, 0.211 8]，主導Floquet乘子|λ|max=1.052 7>1，即|λ|max通過正實數穿過單位圓，則在B點處發生了鞍結分岔(SN)。

由此可得：周期解通過鞍結分岔進入擬周期解，即兩根“細絲”之間的邊界為鞍結分岔曲線。

旋轉數不斷變化，結果形成了“魔鬼樓梯”[32](如圖5(c))。隨著e的不斷增加，r/s不是線性增長，而是一系列的跳躍和不同大小的臺階，不同的臺階代表著不同的鎖模頻率，鎖模的大小則由臺階的寬度來反映。同時，這些臺階具有分形的自相似性。沿著臺階，整個“魔鬼樓梯”平穩地單調遞增，表明周期運動被鎖定在關于e的一個不斷增加且有限范圍內的某個方向上。

圖5 沿著圖4(b)中藍色直線的分岔圖、指數圖及“魔鬼樓梯”Fig.5 The bifurcation diagram,Lyapunov exponent diagram and “Devil’s staircase” along the blue line of Fig.4(b)

圖4(c)放大了第二組共振組織。正如圖6所示，有7個周期窗口嵌入在擬周期區域內：P-11→P-18→P-7→P-17→P-10→P-13→P-16。通過計算在C點處(ω=1 812)的狀態傳遞矩陣的特征值：λ=[0.530 3±0.721 0i,-0.382 4±0.508 4i,0.130 5±0.495 2i]，由3對共軛復數組成，主導Floquet乘子|λ|max=1.002 5>1，即|λ|max通過其中一對共軛復數穿過單位圓，因此在C點處發生了Hopf分岔。由于Hopf分岔的存在導致了共振組織的出現。

圖6 沿著圖4(c)中白色直線的分岔圖及“魔鬼樓梯”Fig.6 The bifurcation diagram and “Devil’s staircase” along the white line of Fig4.(c)

這組共振組織的旋轉數的分布也滿足法里樹序列的拓撲規律，且形成了單調遞增的“魔鬼樓梯”。

圖8(c)～圖(i)顯示了圖7(a)中被淡藍色的A,B,...,G點所標記的參數組合下的強度振蕩的時間演化過程，并通過波形變形說明了新尖峰的產生，揭示了一種按有規律的方式所組織的周期振蕩。綠色垂直線標記振蕩周期T1。圖8(c)～圖8(i)中的時間演化表明了在轉子模式的穩定復雜性中隱藏的規律：每一族似乎都是一些準相同模式的固定組合的串聯，尖峰越多，波形就越多。

圖7 k×v參數平面內穩定振蕩和混沌振蕩的曲折分布Fig.7 The distribution of periodic and chaotic oscillations in k×v parameter plane

圖8 沿著圖7(a)中黃色直線的分岔圖、指數圖及時間歷程圖Fig.8 The bifurcation diagram,Lyapunov exponent diagram and time series diagrams along the yellow line of Fig7.(a)

在圖7(b)中，隨著v的增大，順著水平方向掃描圖像，有8個周期窗口堆積在黑色的擬周期區域內。通過僅僅調整一個參數而得到的單參數分岔圖9(a)，相應地，圖9(b)中有8個負指數被鎖在虛直線以下及圖9(c)中單調遞減的“魔鬼樓梯”。

圖9 沿著圖7(b)中藍色直線的分岔圖、指數圖及“魔鬼樓梯”Fig.9 The bifurcation diagram,Lyapunov exponent diagram and “Devil’s staircase” along the blue line of Fig7.(b)

圖10呈現了參數空間k×ω，k×e，v×e內共振組織的全局拓撲規律。很容易識別出這些不同的雙參數平面具有相同的拓撲規律，即法里樹序列。因此，滿足法里和的共振組織分布是該系統的共性。在旋轉機械中，轉子的質量偏心會產生振動頻率為ω的周期性激振力。若轉子彈性力有非線性就會出現許多非線性共振現象。

圖10 k×ω，k×e，v×e平面內共振組織的全局拓撲機制Fig.10 The global topological mechanism of periodic oscillations of k×ω，k×e，v×e parameter planes

4 結 論

針對大型旋轉機械系統，對其中的動力學現象做詳細的分類是一個費時費力的工作，由于GPU并行計算方法的卓越性能和算法的不斷改進，本文從二維參數角度揭示了汽輪機轉子系統中共振組織排列的拓撲規律，可以得到：滿足法里和序列的共振組織分布是該系統的共性。從非線性動力學角度來看，二維參數平面相圖顯示了波形的復雜性，以及周期吸引子在控制參數空間中積累的速度和位置。因此，這一共性有利于預測在更寬的參數空間動力學行為的演變；從參數匹配的角度來看，周期振蕩的展開會表明系統對某些參數的敏感程度，因此實際試驗也可以從以上的圖表中受益，并提供需要什么樣的精度才能區分相圖中呈現的周期和混沌窗口的復雜交替。從實際應用角度來看，拓撲規律的提取是將復雜的機械現象轉化為簡單的數學問題，為優化設計、參數匹配及故障預測提供了可視化、可比較的圖表，有利于對轉子系統非線性動力學特性更進一步的理解。