基于嵌入式譜隨機有限元法的轉子系統(tǒng)隨機不平衡響應特性分析

周生通， 李鴻光， 張 龍， 周新建

(1.華東交通大學 機電與車輛工程學院，南昌 330013；2.上海交通大學 機械系統(tǒng)與振動國家重點實驗室，上海 200240)

?

基于嵌入式譜隨機有限元法的轉子系統(tǒng)隨機不平衡響應特性分析

周生通1， 李鴻光2， 張 龍1， 周新建1

(1.華東交通大學 機電與車輛工程學院，南昌 330013；2.上海交通大學 機械系統(tǒng)與振動國家重點實驗室，上海 200240)

原始不平衡量在轉軸中沿軸線呈任意空間分布，因而實際轉子的不平衡量需用分布不平衡量和點不平衡量兩種方式共同表達?？紤]兩種不平衡量的隨機性，以隨機場和隨機變量的方式模擬轉子中存在的兩類不平衡量，并采用嵌入式譜隨機有限元方法建立隨機不平衡量下的轉子系統(tǒng)不平衡響應隨機分析模型。結果顯示，提供的隨機表征方法能較好地模擬兩類轉子不平衡量的隨機性；不平衡量的隨機性影響著轉子不平衡響應的變異性，可以發(fā)現：① 共振頻率附近的位移響應幅值標準差也較大；② 位移響應幅值的變異系數受低頻段內共振頻率的影響較小，但隨著共振頻率階次的增高影響也逐漸增大；③ 位移響應相位角的標準差變化情況與位移幅值的變異系數變化情況是基本類似的。

質量不平衡；譜隨機有限元方法；隨機不平衡響應

一個轉子在設計上一般都使它相對于旋轉軸線是軸對稱的。但是由于工藝上的一系列因素，最后裝配完畢的轉子總是不能做到動力上的完全軸對稱，也就是存在一定的不平衡量，這種不平衡量通常稱之為原始不平衡量[1]。造成原始不平衡量的因素主要有：轉子材質的不均勻性，鍵槽不對稱引起的不平衡，轉子加工中總是產生一些圓度偏差和偏心等。對于汽輪機轉子尚有各個葉片之間的差別，葉片鎖口及末葉片的不對稱等不平衡因素。對于壓縮機尚有葉輪的不平衡量。所有這些因素造成的不平衡量都屬于隨機性質的[1]。

像葉片和葉輪位置處的不平衡量是轉子系統(tǒng)不平衡量的主要貢獻源。在理論分析時，一般將葉片和葉輪簡化為剛性圓盤，相應的不平衡量則以力的形式施加在圓盤的幾何中心處，用不平衡質量、偏心距和相位角三個量確定不平衡的大小和方向。不過，如果將連續(xù)轉軸分割為微小厚度的圓片，那么正如單個圓盤一樣，每一個圓片上都存在一定的不平衡量，且大小和方向各不相同。也就是說，轉軸中的不平衡量分布函數實為一條任意的、隨機的空間分布曲線[1-2]。雖然不平衡量的空間任意分布事實很早就被人們認識到，但直到1993年才被LEE等[3]在傳遞矩陣方法中予以考慮。后來，SHIH等[4]又提出了識別柔性轉子不平衡量分布情況的方法。該方法將分段連續(xù)的質量偏心曲線用傅里葉級數展開，但這種方法不太適用于轉子動力學的有限元分析程序，為此YANG等[5]提出采用多項式曲線模擬轉子不平衡質量偏心曲線，以便借助有限元程序更容易地識別轉子中存在的不平衡量分布。DEEPTHIKUMAR等[6]在此基礎上進一步發(fā)展了同時存在分布不平衡和軸彎曲故障的柔性轉子模態(tài)平衡方法。不過，上述這些研究都是從確定性的角度處理轉子不平衡量的。

現有對轉子不平衡量隨機性的考慮多是針對圓盤不平衡量的[7-8]，即將圓盤的偏心距、質量或相位角作為隨機變量處理。不過，STOCKI等[9]在轉子振動離散性分析中則考慮了分布不平衡量的不確定性，其采用的思路是將分布不平衡量表示為轉子前四階振型的加權疊加，權值設為隨機變量，并使用非嵌入式的譜方法開展了相應的隨機分析工作。本文同樣從隨機角度考慮轉子系統(tǒng)中的不平衡量，并以隨機場和隨機變量的形式將它們表達為兩種形式，即集中不平衡量(Lumped Unbalance)和分布不平衡量(Distributed Unbalance)。最后，對某一50 MW汽輪機轉子同時計及這兩種不平衡量，并基于嵌入式譜隨機有限元方法建立該轉子的隨機不確定性分析模型，實現其不平衡響應的隨機分析。

1 轉子不平衡量的隨機表征

1.1 分布不平衡量

圖1 轉子中的不平衡量Fig.1 Rotor unbalances

將分布不平衡量在單個轉子單元內產生的動能代入拉格朗日方程，即可得到對應的不平衡力單元向量：

(1)

(2)

為了處理積分中的隨機量，將一維隨機場u(x,ω)和φ(x,ω)分別作Karhunen-Loeve展開并分別保留至L1階和L2階，有：

(3)

(4)

(5)

式(5)由于包含了連乘的指數函數，進一步處理會比較繁瑣。這里將相位角隨機場φ(x,ω)退化為正態(tài)隨機變量，即?。?/p>

(6)

這樣，式(5)就可簡化為：

(7)

進一步，取指數函數的級數展開近似表達：

(8)

則由式(7)表示的轉子空間隨機不平衡量產生的單元力向量就可整理為如下形式：

(9)

式中:Qkj是維數為8×1的列矩陣，表達式為：

(10)

更特別地，若不平衡量分布密度u(x,ω)也同樣退化為正態(tài)隨機變量，即：

(11)

則有：

(12)

式中:

(13)

1.2 集中不平衡量

對于集中不平衡量，只需用隨機變量就可完全表征不平衡量大小和相位角的隨機性。在文獻[8]中已經給出了當兩者都是正態(tài)隨機變量但相位角均值為零時的不平衡量表征公式。這里則直接將其擴展為不平衡量大小和相位分別是具有任意分布參數的正態(tài)隨機變量的情況，對應的集中不平衡量表征公式為：

(14)

式中:Qkj是維數為4×1的列矩陣，表達式為：

(15)

2 不平衡響應的嵌入式譜隨機有限元方程

譜隨機有限元方程是依據多項式混沌基的正交性質建立的。從轉子隨機不平衡量的表征公式(9)、(12)和(14)中可以看到，作為右端項的不平衡力Fu是以多項式的形式表達的，無法直接利用多項式混沌基的正交性。需要將Fu的表達式改寫成以多項式混沌展開的形式，即：

(16)

將原來Fu的表達式改寫為式(16)的形式后，即可建立起轉子系統(tǒng)不平衡響應的譜隨機有限元方程：

(17)

式中:Z(ω)是轉子系統(tǒng)的動剛度矩陣。在不考慮其他隨機因素的情況下，式(17)最終可轉化為如下P個線性方程組：

Z(ω)δi=Fi,i=0,1,…,P-1

(18)

解方程組即可求得不平衡響應的多項式混沌展開中的廣義坐標值δi,i=0,…,P-1，從而表達出隨機不平衡響應Δ，即為：

(19)

(20)

3 數值算例

求解某一50 MW汽輪機轉子動平衡系統(tǒng)(如圖2所示)的隨機不平衡響應。設汽輪機轉軸中的不平衡量分布密度的大小u(x,ω)服從高斯隨機場且具有指數協方差函數，相應的相位角φ服從正態(tài)隨機量；同時，取作用在圓盤10和圓盤19上的不平衡量U10和U19的大小和相位角亦均服從正態(tài)分布。具體的隨機量設定值如表1所示。

圖2 某50 MW汽輪機轉子動平衡系統(tǒng)Fig.2 A 50 MW steam turbine rotor dynamic balancing system

元件屬性名稱分布類型均值標準差備注圓盤10不平衡量U10正態(tài)分布0.190.038相位角φ10正態(tài)分布0π/15圓盤19不平衡量U19正態(tài)分布0.30.06相位角φ19正態(tài)分布π/3π/15轉軸相位角φ正態(tài)分布π/6π/60分布不平衡量大小u(x)一維高斯隨機場0.010.001設轉軸分布不平衡量隨機場具有指數協方差函數,相關長度6.36m,區(qū)間為[0,6.36]

該模型中，單個集中不平衡量采用2維5階精度描述，分布不平衡量采用5維5階精度描述，故整個譜隨機有限元模型將形成一個具有9維5階的多項式混沌展開模型，共計P=2 002項多項式混沌基函數。進而，依據式(18)～(20)即可得到轉子上任意節(jié)點位移的隨機響應情況。限于篇幅，這里僅給出圓盤10位置處的位移響應概率統(tǒng)計信息，如圖3～7所示。圖3(a)和(b)所示分別為圓盤10處的水平和豎直位移頻響函數幅值的均值、上下包絡線和標準差曲線；圖4給出了對應水平和豎直位移頻響幅值的變異系數隨激勵頻率的變化情況；圖5和圖6則類似地給出了水平和豎直位移相位角的統(tǒng)計信息；圖7(a)和(b)給出了在特定轉速頻率(25.5 Hz和40.0 Hz)下水平位移的直方圖和基于核密度估計的概率密度信息，其中頻率25.5 Hz靠近轉子系統(tǒng)一階臨界轉速，頻率40.0 Hz則遠離共振頻率。

其中，位移均值(實線)、標準差(實線)、上包絡(點線)、下包絡(虛線)和1 000個Monte Carlo樣本(灰線)圖3 圓盤10的位移幅值均值和標準差信息Fig.3 Mean and standard deviation of displacementamplitudes at disk 10

圖4 圓盤10處水平和豎直位移幅值的變異系數曲線Fig.4 The coefficient of variation curves of horizontal and vertical displacement amplitudes at disk 10

從圖3(a)和(b)中可以看到，受到轉子隨機不平衡量的影響，圓盤10處的水平和豎直位移頻響幅值曲線表現出了離散特性，并圍繞著均值曲線上下隨機波動，而偏離均值的大小與激勵頻率有關，當激勵頻率與系統(tǒng)固有頻率一致，即發(fā)生共振時，偏離量最大，即在共振處出現有較大的標準差峰值。但此處的標準差大小并不能準確反映位移頻響幅值隨機離散的程度，圖4則采用變異系數表達離散程度以消除均值的影響。從圖4中可以看到，水平和豎直位移頻響函數幅值的變異系數在小于40 Hz的頻段內基本沒有變化，而這一頻段內就包含有轉子系統(tǒng)的第一階臨界轉速頻率，且此臨界頻率處的響應幅值標準差就比較大。隨著激勵頻率的增加，當達到其他高階臨界轉速時，變異系數就開始有較大的浮動，可以看到階數越高變異系數峰值也就越大。這說明轉子的不平衡隨機因素對系統(tǒng)高階共振頻率處位移幅值的離散程度影響較大，而對低階影響則較小。如在第一階臨界轉速頻率處的影響幾乎可以忽略，因為它與附近非共振頻率處的離散程度幾乎一致。

其中，位移相位角均值(實線)、上包絡(點線)、下包絡(虛線)和1000個Monte Carlo樣本(灰線)圖5 圓盤10的位移相位角信息Fig.5 Phase angles of displacement at disk 10

圖5(a)和(b)則分別給出了隨機不平衡量影響下的水平和豎直位移相位角的離散情況，圖6(a)和(b)中顯示的是對應相位角的標準差和變異系數情況?？梢钥吹?，與位移幅值離散特性表達不同，位移相位角的離散程度用標準差表征更好，且此時的水平和豎直位移相位角標準差的變化情況與對應的水平和豎直位移幅值變異系數的變化情況類似，如圖4所示。其中的原因主要是由于位移幅值均值及其標準差的數值大小在不同的頻率范圍內存在有量級上的較大差別，當表示離散特性時不同量級上的標準差之間不易反映整體離散程度的大小，此時就需要剔除均值的影響用變異系數來表示離散程度，但對于位移相位角來說，其標準差不存在量級差別問題，而均值在一定的范圍內(不超過360°)變化，此時用標準差表示離散程度會更好，若采用變異系數則效果上不佳(如圖6(a)和(b)虛線所示)。

圖6 圓盤10的位移相位角標準差和變異系數頻譜曲線Fig.6 The standard deviation and coefficient of variation curves of displacement phase angle at disk 10

圖7 圓盤10的水平位移在不同轉速頻率下的均值、標準差、概率密度函數直方圖和核密度估計曲線Fig.7 Means, standard deviations, histograms of probability density functions, kernel density estimation curves of the horizontal displacement at disk 10 under different rotating frequencies

進一步，從圖7(a)和(b)中可以看到某一特定激勵頻率下水平位移頻響幅值的概率分布比較接近正態(tài)分布，這與我們在例子中所采用的隨機量類型均為正態(tài)型的有較大關系。其中，在第一階臨界轉速頻率附近(25.5 Hz處)時，其對應位移幅值的均值為1.214×10-3m、標準差為2.006×10-4m；在非共振區(qū)域(40 Hz處)時，均值為4.728×10-5m、標準差為7.960×10-6m。從量級上對比，可以看到兩個頻率處的標準差在數值上相差有兩個數量級，但這也只能表明系統(tǒng)在共振頻率處的響應離散幅值更大；而從幅值離散程度上對比，可以計算出兩頻率處的幅值變異系數分別為0.165 2和0.168 4，比較接近，說明實際中這兩個頻率處的隨機離散程度差別并不大，與圖4結果一致。

4 結 論

本文采用分布不平衡量和集中不平衡量兩種概念描述轉子中存在的不平衡量。在考慮不平衡量隨機性的基礎上，推導了分布不平衡量和集中不平衡量的力表達式，并利用嵌入式譜隨機有限元方法建立了轉子系統(tǒng)不平衡響應的隨機分析模型，求解了某一50 MW汽輪機轉子在動平衡擺架支撐系統(tǒng)上的隨機不平衡響應概率信息。總結算例結果，可以看到：

(1) 受到轉子隨機不平衡量的影響，系統(tǒng)不平衡響應(位移幅值頻響函數)表現出了較大的隨機離散特性，尤其是在共振頻率附近，其對應的響應幅值標準差要明顯高于非共振頻率區(qū)；

(2) 但從位移頻響幅值的變異系數上分析，在低頻段內各頻率上的變異系數基本保持不變(即使是在第一階臨界轉速頻率處)，而隨著激振頻率的增大，變異系數的浮動也開始增大，尤其是在高階次的共振頻率處，變異系數出現了較大峰值，這說明了轉子不平衡隨機因素對更高階共振頻率處的位移幅值離散程度影響更為明顯；

(3) 對于位移相位角來說，其也表現出了較大的隨機離散特性，但其離散程度更適合用相位角標準差來表征。

[1] 周仁睦. 轉子動平衡——原理、方法和標準[M]. 北京: 化學工業(yè)出版社, 1992.

[2] 王悅武，田杜平，徐錫林. 高速平衡技術裝備[M]. 上海: 上海交通大學出版社, 2013.[3] LEE A, SHIH Y, KANG Y. The analysis of linear rotor-bearing systems: a general transfer matrix method[J]. Journal of Vibration and Acoustics,1993, 115(4): 490-497.

[4] SHIH Y, LEE A. Identification of the unbalance distribution in flexible rotors[J]. International Journal of Mechanical Sciences,1997, 39(7): 841-857.

[5] YANG T, LIN C. Estimation of distributed unbalance of rotors[J]. Journal of Engineering for Gas Turbines and Power,2002,124(4): 976-983.

[6] DEEPTHIKUMAR M B, SEKHAR A S, SRIKANTHAN M R. Modal balancing of flexible rotors with bow and distributed unbalance[J]. Journal of Sound and Vibration,2013,332(24):6216-6233.

[7] ZHANG Y, WEN B, LEUNG A Y T. Reliability analysis for rotor rubbing[J]. Journal of Vibration and Acoustics,2001, 124(1): 58-62.

[8] DIDIER J, FAVERJON B, SINOU J. Analysing the dynamic response of a rotor system under uncertain parameters by polynomial chaos expansion[J]. Journal of Vibration and Control,2012, 18(5): 712-732.

[9] STOCKI R, LASOTA R, TAUZOWSKI P, et al. Scatter assessment of rotating system vibrations due to uncertain residual unbalances and bearing properties[J]. Computer Assisted Methods in Engineering and Science,2012, 19(2): 95-120.

Stochastic unbalance response characteristics of rotor systems based on intrusive spectral stochastic finite element method

ZHOU Shengtong1, LI Hongguang2, ZHANG Long1, ZHOU XinJian1

(1. School of Mechatronics & Vehicle Engineering, East China Jiaotong University, Nanchang 330013, China;2. State Key Lab of Mechanical Systems & Vibration, Shanghai Jiao Tong University, Shanghai 200240, China)

Mass unbalance in a flexible shaft is distributed arbitrarily along its axis, so distributed mass unbalance and lumped mass unbalance should be considered and presented, respectively for mass unbalance simulation of real rotor systems. Taking the random nature of these two unbalances in a flexible shaft into account, they are represented as random fields and random variables. Here, stochastic analysis of unbalance response of rotor systems with random mass unbalances were performed using the intrusive spectral stochastic finite element method. Numerical results illustrated that the presented stochastic representation methods of distributed and lumped mass unblances are suitable to simulate the random nature of mass unbalances. Results showed that the variations of unblance responses are largely affected by the randomness of unbalances, for example, ① the standard devitations of displacement amplitudes near the resonance frequencies become larger; ② the effects of resonance frequences with in a lower frequency band on the variation coefficient of displacement amplitudes are less, but the effects increase with increase in resonance frequency order; ③ the change curve of the standard deviation of displacement phase angle is similar to that of the variation coefficient of displacement amplitudes.

mass unbalance; spectral stochastic finite element method; stochastic unbalance response

國家自然科學基金(51505146)；江西省自然科學基金(20161BAB216135;20122BAB206027)

2015-07-02 修改稿收到日期：2015-09-27

周生通 男，博士，講師，1984年生

李鴻光 男，博導，教授，1972年生

E-mail: hgli@sjtu.edu.cn;zhoust@ecjtu.edu.cn

TH113

A

10.13465/j.cnki.jvs.2016.19.008