一種分析具有局部非線性動力系統(tǒng)的二次降維方法

楊雄民 陳可軒 閆志勇

(1.杭州汽輪機股份有限公司工業(yè)透平研究院；2.浙江省工業(yè)汽輪機轉子動力學研究重點實驗室）

0 引言

對于轉子-軸承耦合系統(tǒng)以及其他一些具有非線性結合面的局部非線性動力系統(tǒng)，其系統(tǒng)維數較高,可以達到103以上，雖然非線性自由度只占系統(tǒng)自由度的很少一部分，但卻決定著系統(tǒng)的動力特性。直接利用數值方法求解上述系統(tǒng)固有值問題和非線性動力特性，不僅計算量巨大，而且由于舍入誤差、求解精度不協調等造成迭代不收斂、軸頸碰擦等迫使計算中斷。

圍繞系統(tǒng)的降維，眾多學者進行了大量的研究，其主要思想是利用模態(tài)綜合方法降低系統(tǒng)維數，進而通過求解降階方程來討論原系統(tǒng)的固有值問題和非線性動力特性。Neleson 等首先引入了模態(tài)綜合思想，用于復雜轉子系統(tǒng)的穩(wěn)定性分析和瞬態(tài)響應分析[1]。Rouch和Kao 利用Guyan 直接縮聚方法來降低轉子有限元模型的自由度規(guī)模[2]。進一步地，張家忠等提出局部非線性力之作用在系統(tǒng)的個別節(jié)點上，把系統(tǒng)的自由度分為線性自由度和非線性自由度，利用模態(tài)縮聚方法縮減線性自由度而把非線性自由度完全保留在物理空間中[3,4]。對于局部非線性動力系統(tǒng)，通常關心其靜平衡位置的擾動周期解及其分岔，然而上述模態(tài)變換方法模態(tài)變換矩陣僅僅與結構的質量矩陣、阻尼矩陣和剛度矩陣有關，而沒有考慮局部非線性力的影響。由于局部非線性力一般為系統(tǒng)在靜平衡位置上擾動位移和擾動速度的函數，與系統(tǒng)的各固有特性矩陣相耦合，共同決定著系統(tǒng)的運動模態(tài)形態(tài)。鄭鐵生等首先通過變分原理計算非線性油膜力及其Jacobin 矩陣，在形成系統(tǒng)的阻尼矩陣，剛度矩陣時考慮油膜力的影響，形成系統(tǒng)的模態(tài)變換矩陣，進而對原系統(tǒng)進行降階[5]?？紤]局部非線性力的影響，系統(tǒng)質量矩陣一般為對稱矩陣，而剛度矩陣和阻尼矩陣為非對稱大型帶狀稀疏矩陣。傳統(tǒng)的QR方法勢必破壞系統(tǒng)矩陣的上述特點，花費大量的機時和存儲空間。Meirovitch方法針對無阻尼回轉系統(tǒng)K，M 矩陣對稱，D 矩陣反對稱情況，無法推廣到一般系統(tǒng)[6]。非對稱型的Lanczos方法數值穩(wěn)定性差，而對稱性廣義Lanczos 方法勢必引入大量的復數運算[7]。鄭鐵生等提出了一種非對稱矩陣結構系統(tǒng)固有值分析的廣義逆迭代方法，不需要把系統(tǒng)改寫為一個2n階線性問題，直接在原n階規(guī)模上進行反迭代，從而把高維系統(tǒng)固有值問題降階為一個小型線性標準特征值問題[8]。鄭兆昌等把阻尼陀螺系統(tǒng)運動方程轉化為一階狀態(tài)方程，利用Arnoldi方法對系統(tǒng)進行降階，進而通過Schur分解把系數矩陣轉化為上三角矩陣來求解系統(tǒng)的特征值及動力響應[9]。王陶提出了針對局部非線性問題的混合坐標自由界面子結構模態(tài)綜合方法[10]。

由于截斷誤差以及降階模型階數選擇不合理等因素，使得降階后系統(tǒng)特征值出現偽根，在利用數值方法計算降維后，系統(tǒng)方程長期動力特性時發(fā)散。本文基于廣義反迭代思想，提出一種可以用于求解局部非對稱非線性系統(tǒng)動力特性的二次降維方法。通過二次降維，把原高維局部非線性動力系統(tǒng)降階為一個等價的低維系統(tǒng)。針對降維系統(tǒng)周期解的求解，給出其配套打靶法算法。最后以滑動軸承-轉子動力系統(tǒng)為例，運用本文算法求解其非線性動力響應。

2 局部非線性動力系統(tǒng)的二次降維

不失一般性，令局部非線性動力系統(tǒng)方程為

其中，M為系統(tǒng)的質量矩陣，G為系統(tǒng)的阻尼矩陣(陀螺矩陣)，K為系統(tǒng)的剛度矩陣，g(t)為線性激勵力，f為系統(tǒng)的非線性力項，只作用在系統(tǒng)的某些局部自由度上，可以表示為

考慮到系統(tǒng)的靜平衡方程

其中，g*為靜態(tài)線性力，系統(tǒng)運動方程可以簡化為

取變換矩陣P和Q滿足代入方程（5），可以得到，

存在上Hessenberg矩陣H∈?m×m滿足[8]

結合方程(6)和方程(7)可得

根據方程(11)可以發(fā)現，二次降維后系統(tǒng)自由度的維數為k，遠小于原系統(tǒng)(1)的自由度n，可以利用直接數值方法如Runge-Kutta法等來研究降維后系統(tǒng)的動力特性。

3 非線性動力特性求解及算法

研究系統(tǒng)平衡位置的周期解及其穩(wěn)定性具有重要的工程和理論意義。直接采用基于Poincaré點映射的數值積分方法，由于求解局部非線性力及其Jacobin 矩陣如轉子-軸承系統(tǒng)需要求解Reynolds 方程，往往要花費大量的機時。構建合適的算法來求解降維系統(tǒng)的周期解，則成為本文算法能否用于工程實際的關鍵點。二次降維系統(tǒng)非線性動力方程簡記為

一般滿足非自治系統(tǒng)G(η,t)=G(η,t+T0)。在不引起混淆的前提下定義P為Poincaré 映射，設ηj為系統(tǒng)周期l解在Poincaré 截面上的不動點，則滿足

利用Newton-Raphson方法構造如下的迭代格式

對于迭代格式(12)中DηP(ηj)的計算，根據Poincaré映射定義可以按照式（15）進行計算

在t=T0時，Φ(T0)=DηP(ηj)即為Poincaré 映射在ηj處的Jacobin矩陣，可以根據Φ(T0)來判斷周期解的穩(wěn)定性。二次降維系統(tǒng)方程（11）和方程（15）同時迭代計算，進而可以根據式（14）打靶得到系統(tǒng)的周期解。

4 數值算例

轉子-軸承系統(tǒng)為一典型的局部非線性動力系統(tǒng)。為了驗證本文算法，對某一工程實際轉子進行非線性動力特性分析。圖1某工程軸承-轉子系統(tǒng)動力學模型，轉子各段長度及其外徑如表1 所示，由兩個相同的橢圓瓦軸承支承（B1 和B2，軸頸直徑140mm，瓦塊包角為150°，長徑比為1.0，間隙比為0.0025，瓦塊預負載系數為0.3，潤滑油粘度為0.022N·Sec/m2），附加葉輪（D1-D5）質量及轉動慣量如表2 所示。不平衡量施加在圓盤D2 和D4 處，相位分別為0°和180°，大小為3000g·mm。

圖1 某工程軸承-轉子系統(tǒng)動力學模型Fig.1 Dynamic model of rotor-bearing system

表1 葉輪質量及轉動慣量Tab.1 Impeller mass and moment

取轉速ω=3919rpm 和零初始條件計算其周期解殘差，表2 給出了部分降階模型誤差結果對比，可以發(fā)現本文算法通過二次降維，避免了因降維階數m設置不合理而出現偽根，具有較高的計算精度。表3 為ω=3000rpm時利用打靶法求解系統(tǒng)周期解的迭代過程，其中第一次打靶初始條件為η=00，可以發(fā)現本文周期解打靶收斂速度極快。

表2 不同降階模型計算響應的誤差Tab.2 Calculated response errors of different reduced order model

表3 打靶法求解周期解誤差Tab.3 Periodic relative errors calculated of shooting method

圖2為以轉速為參數系統(tǒng)分岔圖，其中橫軸轉速從2000rpm 到20000rpm，縱軸為B1 點x方向位移。從圖2可以看出，轉速小于3919rpm時，為周期1解；到轉速大于3919rpm 時，Floquent 乘子在-1處通過復平面上的單位圓，周期1 解失穩(wěn)而發(fā)生倍周期分岔，分岔后系統(tǒng)由不穩(wěn)定的周期解變?yōu)榉€(wěn)定的倍周期解。當轉速增至4400rpm 時，倍周期解的一對共軛Floquent 乘子通過復平面上的單位圓，發(fā)生準周期分岔。隨著轉速的增加，系統(tǒng)基本呈現為準周期解。從圖2可以看出，在某些轉速下仍有一些周期窗口出現，即出現各種亞諧振動（模態(tài)鎖定現象）。圖3為轉速6869rpm時系統(tǒng)B1點的軸心軌跡和Poincaré 圖。由圖3 可以看出，Poincaré 截面上的不動點已突變?yōu)榄h(huán)面集，系統(tǒng)為準周期運動。圖4為轉速18809rpm 系統(tǒng)的1/7 亞諧運動。圖6 為轉速14000rpm 系統(tǒng)B1 點的軸心軌跡、時間歷程曲線、Poincaré 圖和功率譜圖，可以發(fā)現其吸引子已呈現明顯的分數維特征，響應曲線包含了豐富的低頻特征，系統(tǒng)已進入混沌運動。

圖2 B1點分岔圖Fig.2 Bifurcation of Point B1

圖3 準周期運動Fig.3 Quasi-periodic motion

圖4 1/7亞諧運動Fig.4 1/7 sub harmonic motion

圖5 混沌運動Fig.5 Chaotic motion

5 結論

本文針對局部非線性動力系統(tǒng)質量、阻尼和剛度矩陣大型稀疏非對稱特點，充分考慮局部非線性力、截斷誤差等因素，提出了一種二次降維方法及其配套的周期解求解方法，把原高維系統(tǒng)降階為等價的低維系統(tǒng)，使得系統(tǒng)的維數得到大大降低。整個降階過程僅在n維系統(tǒng)上進行，第一次降階不涉及復數運算，第二次降維只需要求解一個低維特征值問題。數值試驗表明，本文算法具有比較高的數值精度和收斂精度，有利于認識透平機械非正常工況亞諧和超諧非線性振動。