基于非線性支承參數(shù)的渦軸發(fā)動機燃發(fā)轉(zhuǎn)子不平衡量識別*

廖三豐 陳 成

（中國航發(fā)湖南動力機械研究所；航空發(fā)動機振動技術(shù)航空科技重點實驗室）

0 引言

渦軸發(fā)動機燃氣發(fā)生器(簡稱燃發(fā)轉(zhuǎn)子)是典型的高速兩支點轉(zhuǎn)子，兩支點指的是其支承布局為1-0-1，其工作轉(zhuǎn)速一般超過40 000RPM，往往在二階臨界轉(zhuǎn)速以上。轉(zhuǎn)子不平衡微弱變化將產(chǎn)生較大的不平衡激振力，所以前后彈性支承采用的是擠壓油膜阻尼結(jié)構(gòu)（SFD），轉(zhuǎn)子不平衡振動產(chǎn)生的動能通過SPD轉(zhuǎn)化為滑油的內(nèi)能被“帶走”，達到減振的目的。而SFD的阻尼剛度表現(xiàn)出明顯的非線性特征，在轉(zhuǎn)子不平衡量識別過程中如何處理支承參數(shù)的非線性問題一直是一個較大的難題。

上個世紀(jì)70 年代以來，發(fā)展了許多基于線性理論的無試重動平衡的方法[1～4]，無論是振型平衡法還是影響系數(shù)法都是建立在線性理論基礎(chǔ)之上的，此類方法都以降低轉(zhuǎn)子響應(yīng)為目標(biāo)，不關(guān)注轉(zhuǎn)子實際工作中不平衡量惡化規(guī)律，無法實現(xiàn)對轉(zhuǎn)子在全工作轉(zhuǎn)速內(nèi)的不平衡量識別。在旋轉(zhuǎn)機械故障診斷方面，也無法準(zhǔn)確考慮到支承結(jié)構(gòu)非線性因素的影響，轉(zhuǎn)子不平衡過大的情況只停留在定性的層面上[5-7]。國內(nèi)外已有大量關(guān)于轉(zhuǎn)子不平衡在線識別的研究，饒柱石[8]等基于參數(shù)識別、子結(jié)構(gòu)模態(tài)分析及阻抗匹配原理，通過建立轉(zhuǎn)子—軸承系統(tǒng)頻響函數(shù)實現(xiàn)對轉(zhuǎn)子不平衡量在線識別，該方法在識別過程中因計算規(guī)模小而具有快速的特征。畢士華等[9]提出了一種考慮了油膜力的轉(zhuǎn)子不平衡量時域方法，該方法通過求解線性矩陣方程組，無法給出不平衡量具體分布及相位關(guān)系。張召翠[10]等，基于遺傳算法開展非線性軸承參數(shù)及不平衡量參數(shù)識別。劉淑蓮等[11]采用全息譜技術(shù)提取轉(zhuǎn)子共頻分量，從而識別非線性系統(tǒng)的不平衡量。以上識別方法計算復(fù)雜、識別變量多，在渦軸發(fā)動機實際工程應(yīng)用中難以推廣。

目前，在發(fā)動機整機動力學(xué)仿真分析研究中和轉(zhuǎn)子不平衡惡化的在線監(jiān)控的工程應(yīng)用中，缺乏一套高效的不平衡量識別方法[12-13]。本文開展的基于非線性支承參數(shù)高速轉(zhuǎn)子不平衡量識別研究實現(xiàn)了全工作轉(zhuǎn)速下轉(zhuǎn)子不平衡的“可視化”。本文以轉(zhuǎn)子臺上轉(zhuǎn)子測點工頻信號為目標(biāo)，開展非線性支承參數(shù)下高速轉(zhuǎn)子不平衡量參數(shù)識別方法研究。

1 渦軸發(fā)動機燃發(fā)轉(zhuǎn)子動力學(xué)參數(shù)

1.1 燃發(fā)轉(zhuǎn)子不平衡量變化特征

渦軸發(fā)動機燃氣發(fā)生器轉(zhuǎn)子為典型的高速剛性轉(zhuǎn)子，見圖1。在裝配中，中心拉桿上的一、二級壓緊螺母給整個轉(zhuǎn)子施加軸向預(yù)緊力，但在工作過程中會產(chǎn)生損失，主要原因是：一發(fā)動機轉(zhuǎn)子在工作過程中，軸流壓氣機部分主要受到逆氣流方向的軸向力，渦輪部分受到氣流方向的軸向力，抵消了一部分預(yù)緊力；二高速轉(zhuǎn)子的超高工作轉(zhuǎn)速，巨大的離心力會引起轉(zhuǎn)子變形，由于Poisson效應(yīng)會造成零部件發(fā)生軸向收縮。變小的端齒壓緊力造成“壓不住”的情況發(fā)生，連接結(jié)構(gòu)不穩(wěn)定，轉(zhuǎn)子不平衡量隨之發(fā)生非線性變化。

圖1 燃氣發(fā)生器轉(zhuǎn)子結(jié)構(gòu)簡圖Fig.1 Schematic diagram of gas generator rotor structure

1.2 燃發(fā)轉(zhuǎn)子動力學(xué)參數(shù)的等效

圖2 渦軸發(fā)動機燃發(fā)轉(zhuǎn)子模型，由徑向葉輪盤（紅色）、離心葉輪盤（綠色）和燃氣渦輪葉片盤（藍色）組成，工作轉(zhuǎn)子前三階振型分別為平動、擺動和彎曲振型。圖2所列出的為工作中的轉(zhuǎn)子振型，由于轉(zhuǎn)子完整臨界轉(zhuǎn)速較高，工作狀態(tài)下轉(zhuǎn)子以第二階振型為主。

圖2 燃發(fā)轉(zhuǎn)子示意圖Fig.2 Schematic of gas generator rotor structure

實際中的轉(zhuǎn)子客觀存在的不平衡量在轉(zhuǎn)子工作時會產(chǎn)生較大的不平衡激振力，這是轉(zhuǎn)子振動的主要的激振源。燃發(fā)轉(zhuǎn)子可簡化為圖3，圖中ε 為工作中的轉(zhuǎn)子偏心量，k表示支承剛度，c表示支承阻尼。

圖3 燃發(fā)轉(zhuǎn)子不平衡響應(yīng)計算用圖Fig.3 Diagram for calculation of gas generator rotor unbalanced response

1.3 燃發(fā)轉(zhuǎn)子支承參數(shù)的非線性

燃發(fā)轉(zhuǎn)子前后支承一般是帶有鼠籠式彈性支承的定心式SFD，其中結(jié)構(gòu)形式見圖4。

圖4 SFD結(jié)構(gòu)示意圖Fig.4 Schetch of SFD

滑油通過供油孔進入供油槽，滑油在彈性支承與阻尼器之間的間隙中形成一圈油膜，轉(zhuǎn)子運轉(zhuǎn)產(chǎn)生偏心量，受擾動油膜產(chǎn)生的油膜力阻礙轉(zhuǎn)子發(fā)生偏心，流動的滑油把它“吸收”的轉(zhuǎn)子振動動能“帶”出去。SFD剛度阻尼參數(shù)理論推導(dǎo)方法常采用“短”“長”軸承假設(shè)，阻尼器長度與直徑之比小于0.25，兩端不封嚴(yán)時，采用短軸承假設(shè)；密封較好則采用長軸承假設(shè)，具體推導(dǎo)過程及計算公式在相關(guān)設(shè)計手冊已有推導(dǎo)，本文不再贅述。

基于轉(zhuǎn)子試驗臺中振動數(shù)據(jù)的不平衡量識別過程中的k采用彈性支承剛度，而基于整機試車中振動數(shù)據(jù)的轉(zhuǎn)子不平衡量識別中的k應(yīng)包含機匣的剛度，是一個隨工頻變化的動剛度。剛度、阻尼的非線性會直接影響發(fā)動機轉(zhuǎn)子不平衡量識別的準(zhǔn)確性。

1.4 不平衡量識別中參數(shù)描述

考慮動剛度及非線性油膜剛度阻尼的轉(zhuǎn)子運動方程可以寫成

式中，q 為轉(zhuǎn)子廣義坐標(biāo)系矩陣；m 為高速轉(zhuǎn)子系統(tǒng)質(zhì)量矩陣；d 為包含系統(tǒng)內(nèi)阻尼和陀螺力矩矩陣；se為包含系統(tǒng)動剛度的矩陣。由于篇幅所限，以下參數(shù)介紹不帶下標(biāo)，α 為待識別不平衡量位置矩陣（一般選取轉(zhuǎn)子動平衡去材料界面位置）；Fi為不平衡量大??；ω 為轉(zhuǎn)子工作轉(zhuǎn)速；φi表示不平衡量角度；β 為SFD位置矩陣；f 為油膜力。SFD 剛度阻尼通過試驗實測或理論公式計算得到。本文轉(zhuǎn)子不平衡識別主要用于整機動力學(xué)研究的變量輸入，不涉及動平衡試驗，則各位置等效不平衡量周向角度關(guān)系以相對角度表示。待識別不平衡參數(shù)為Fi與φi，共計2n-1個。

2 轉(zhuǎn)子系統(tǒng)參數(shù)識別

本文采用建立的高速轉(zhuǎn)子ANSYS 有限元模型，開展支承參數(shù)識別，并以ISIGHT 為平臺，調(diào)用ANSYS 軟件，開展發(fā)動機轉(zhuǎn)子不平衡量參數(shù)識別。

2.1 轉(zhuǎn)子的響應(yīng)

當(dāng)SFD結(jié)構(gòu)確定后，其阻尼、剛度是供油壓力、滑油粘度、轉(zhuǎn)子工作轉(zhuǎn)速和偏心率函數(shù)。一般情況下，供油壓力和滑油粘度是確定的，工作狀態(tài)下SFD的阻尼與剛度僅是轉(zhuǎn)子工作轉(zhuǎn)速及偏心率的函數(shù)。非線性支承參數(shù)轉(zhuǎn)子不平衡響應(yīng)計算流程見圖5，通過不斷修正支承剛度、阻尼使轉(zhuǎn)子不平衡響應(yīng)達到一個穩(wěn)定值。由于試驗數(shù)據(jù)的離散性，識別過程中基于實測SFD的剛度與阻尼數(shù)據(jù)采用線性插值獲取。每次迭代得到的SFD 的剛度與阻尼SFD 帶入到動力學(xué)有限元模型中即可計算得到支承位置的油膜力，帶入到式（1）即可實現(xiàn)對動力學(xué)方程的求解。

圖5 支承參數(shù)識別流程Fig.5 Support parameter identification process

2.2 不平衡量識別

轉(zhuǎn)子不平衡量識別過程可歸結(jié)為優(yōu)化問題。獲取轉(zhuǎn)子測點響應(yīng)有效值。參數(shù)識別過程的數(shù)學(xué)模型為:

數(shù)學(xué)模型中：Ex 表示試驗中測點響應(yīng)有效值向量；Re 表示仿真獲取的對應(yīng)試驗測點響應(yīng)有效值向量；V 表示仿真值向量與計算值向量之差的2范數(shù)；F 表示待識別的不平衡量；f表示變量上下邊界；φ 為不平衡量相對相角；φ 表示相對相角上下邊界，一般情況下取0～2π。

由于轉(zhuǎn)子支承參數(shù)的高度非線性，本文選用Pointer優(yōu)化器，該求解器能夠針對高度非線性和非連續(xù)性設(shè)計空間問題，通過自動捕捉設(shè)計空間的信息，靈活地組合算法形成一個參數(shù)識別策略。

3 算例研究

以某渦軸發(fā)動機高速轉(zhuǎn)子為不平衡量識別對象，轉(zhuǎn)子總質(zhì)量為20.303 0kg，軸慣性矩Ixc=62 848.30kg?mm2，直徑慣性矩Iyc=429375.08kg?mm2，Izc=429375.08kg?mm2。基于ANSYS建立該轉(zhuǎn)子有限元模型，見圖6。

圖6 燃發(fā)轉(zhuǎn)子有限元模型Fig.6 Finite element model of gas generator rotor

3.1 非線性支承參數(shù)轉(zhuǎn)子的響應(yīng)

基于建立的轉(zhuǎn)子有限元模型，結(jié)合圖5中非線性支承參數(shù)下轉(zhuǎn)子不平衡響應(yīng)計算方法，對比仿真給定不平量下定支承參數(shù)與變支承參數(shù)的轉(zhuǎn)子響應(yīng)。

本算例中前等效不平衡量位置為轉(zhuǎn)子一級壓氣機盤位置，大小為6g?mm，后等效不平衡量為渦輪盤位置，大小為8g?mm，前后不平衡量相對相位差為0 rad。為模擬轉(zhuǎn)子實際工作狀態(tài)，把該轉(zhuǎn)子模型置于整機環(huán)境中，建立靜子支承系統(tǒng)超單元有限元模型，整機中的燃發(fā)轉(zhuǎn)子見圖7。

圖7 整機中的燃發(fā)轉(zhuǎn)子Fig.7 Gas generator rotor in turboshaft engine

SFD參數(shù)見表1，其中潤滑油密度ρ=904kg/m3，運動粘度υ=4.2×10-6m2/s。

表1 轉(zhuǎn)子前后支承參數(shù)Tab.1 Parmeters of rotor supporting

提取圖6 中1#測點響應(yīng)，仿真結(jié)果對比圖見圖8，變阻尼的轉(zhuǎn)子在過臨界時測點的不平衡響應(yīng)明顯低于定阻尼情況，說明SFD可有效抑制轉(zhuǎn)子的振動。

圖8 1#測點不平衡響應(yīng)曲線Fig.8 Unbalance response curve of 1#measuring point

3.2 不平衡量識別

轉(zhuǎn)子臺上前后彈支剛度為1.72×107N/m 和2.15×107N/m。前后平衡截面等效不平衡力分別為：

以該轉(zhuǎn)子在轉(zhuǎn)子試驗臺上試驗數(shù)據(jù)為目標(biāo)值，提取轉(zhuǎn)速32 900RPM、36 000RPM、38 000RPM、41 600 RPM和42 500RPM轉(zhuǎn)速下轉(zhuǎn)子1#測點與2#測點位移響應(yīng)有效值，見表2。

表2 測點位移響應(yīng)幅值Tab.2 The displacement response amplitude of the measuring point

本算例采用的仿真計算機為Intel i5 四核處理器，主頻3.4GHz，8GB內(nèi)存容量。

在識別程序中輸入轉(zhuǎn)速、前后軸承位移響應(yīng)有效值，以F1,F2和φ 為待識別參量，開展不平衡量識別研究。為縮減篇幅，本文只列出轉(zhuǎn)速32 900RPM下轉(zhuǎn)子不平衡量識別過程，在參數(shù)識別程序中分別輸入1#、2#測點位移響應(yīng)，識別過程見圖9 和圖10，每個迭代步中嵌套有圖5中支承參數(shù)識別流程，不同平衡量下識別的油膜阻尼、剛度也存在差異性，整個不平衡量及支承參數(shù)識別過程約3分鐘，其它轉(zhuǎn)速下識別時間為2～4分鐘。

圖9 1#測點位移響應(yīng)收斂曲線Fig.9 Displacement response convergence curve of 1#measuring point

圖10 2#測點位移響應(yīng)收斂曲線Fig.10 Displacement response convergence curve of 2#measuring point

輸入其他轉(zhuǎn)速及對應(yīng)的目標(biāo)值，識別結(jié)果見表3，繪制工作轉(zhuǎn)速內(nèi)轉(zhuǎn)子不平衡量變化曲線圖。圖11 表明：轉(zhuǎn)子不平衡量隨高速轉(zhuǎn)子轉(zhuǎn)速上升呈現(xiàn)放大的趨勢，高速轉(zhuǎn)子前后不平衡量相位差約180度。

表3 不平衡量識別結(jié)果Tab.3 The identification result of the unbalance

圖11 工作轉(zhuǎn)速內(nèi)轉(zhuǎn)子不平衡量Fig.11 Unbalance of the rotor in the working speed

該類發(fā)動機轉(zhuǎn)子裝機前不平衡量會控制在10g?mm以內(nèi)，工作中的轉(zhuǎn)子不平衡量往往會放大數(shù)十倍。通過本算例的研究發(fā)現(xiàn)，工作中的轉(zhuǎn)子不平衡量比裝機前需用不平衡量放大了10倍到20倍。

4 結(jié)論

本文建立起渦軸發(fā)動機燃發(fā)轉(zhuǎn)子有限元動力學(xué)模型，選取轉(zhuǎn)子動平衡位置處的等效不平衡量及相對相角為待識別變量，面向轉(zhuǎn)子臺中的轉(zhuǎn)子上傳感器測點工頻信號，基于ANSY編寫包含SFD非線性的轉(zhuǎn)子不平衡響應(yīng)仿真程序，基于ISIGHT 軟件采用Pointer 優(yōu)化器完成渦軸發(fā)動機高速兩支點轉(zhuǎn)子不平衡量的識別，研究發(fā)現(xiàn)：

1)本文所提出的不平衡量識別方法的識別效率較高，在文中的算例中，單個工作轉(zhuǎn)速下不平衡量識別的時間為2～4min，該方法在工程設(shè)計中具有較強的可操作性，識別結(jié)果較為準(zhǔn)確地反映實際轉(zhuǎn)子不平衡量的變化趨勢；

2)本文的算例研究發(fā)現(xiàn)，包含SFD非線性的轉(zhuǎn)子不平衡響應(yīng)仿真程序較為準(zhǔn)確地模擬轉(zhuǎn)子過臨界時響應(yīng)特征；

3)本文所識別得到的結(jié)果可作為整機振動仿真分析中輸入，為發(fā)動機整機不平衡響應(yīng)研究中的不平衡量的施加提供參考。