楊雄民 陳可軒 閆志勇

(1.杭州汽輪機(jī)股份有限公司工業(yè)透平研究院；2.浙江省工業(yè)汽輪機(jī)轉(zhuǎn)子動(dòng)力學(xué)研究重點(diǎn)實(shí)驗(yàn)室）

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

3 非線性動(dòng)力特性求解及算法

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

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

利用Newton-Raphson方法構(gòu)造如下的迭代格式

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

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

4 數(shù)值算例

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

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

表1 葉輪質(zhì)量及轉(zhuǎn)動(dòng)慣量Tab.1 Impeller mass and moment

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

表2 不同降階模型計(jì)算響應(yīng)的誤差Tab.2 Calculated response errors of different reduced order model

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

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

圖2 B1點(diǎn)分岔圖Fig.2 Bifurcation of Point B1

圖3 準(zhǔn)周期運(yùn)動(dòng)Fig.3 Quasi-periodic motion

圖4 1/7亞諧運(yùn)動(dòng)Fig.4 1/7 sub harmonic motion

圖5 混沌運(yùn)動(dòng)Fig.5 Chaotic motion

5 結(jié)論

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