萬 周 劉璟澤 張大海 陳 強 唐振寰 費慶國

(1東南大學機械工程學院， 南京 211189)(2東南大學江蘇省空天機械裝備工程研究中心， 南京 211189)(3中國航空發(fā)動機集團有限公司湖南動力機械研究所， 株洲 412002)

航空發(fā)動機的高保真建模是其研制過程中的重要問題[1].支承參數(shù)是航空發(fā)動機動態(tài)特性的關鍵影響因素，對于航空發(fā)動機整機動力學建模至關重要，而支承結構中的非線性因素對支承參數(shù)識別存在一定影響，因此有必要開展考慮非線性的支承參數(shù)識別研究.

準確有效的航空發(fā)動機部件模型是精確識別支承參數(shù)的基礎.近年來，航空發(fā)動機部件建模技術不斷發(fā)展，基于梁單元的航空發(fā)動機建模方法利用非線性力模型表征支承結構中的非線性因素，目前已成為十分有效的建模方法[2-5].

獲得準確的轉靜子模型后，一些學者針對航空發(fā)動機支承參數(shù)識別問題開展了相關研究.屈美嬌等[6]運用支持向量機(SVM)與遺傳算法實現(xiàn)了航空發(fā)動機靜態(tài)支承剛度識別.費慶國等[7]提出了一種基于徑向基神經(jīng)網(wǎng)絡(RBFNN)的有限元模型修正方法.以上研究雖未能實現(xiàn)復雜非線性結構的參數(shù)識別，但其提出的機器學習參數(shù)識別方法為解決復雜非線性結構的動力學反問題提供了新的研究思路.

隨著機器學習的發(fā)展，性能強大的神經(jīng)網(wǎng)絡層出不窮[8-10].長短期記憶(LSTM)神經(jīng)網(wǎng)絡在處理時間序列時表現(xiàn)出優(yōu)秀的長期預測能力，使得求解復雜非線性結構的動力學反問題成為可能.

本文提出一種基于LSTM的航空發(fā)動機整機支承剛度識別方法.首先，基于整機模型獲得不同支承剛度對應的位移響應；然后，使用深度學習網(wǎng)絡構建支承剛度與位移響應之間的非線性關系；最后，利用訓練好的深度學習網(wǎng)絡實現(xiàn)支承剛度識別.本文通過與真實值的對比驗證了該方法的有效性，通過與基于RBFNN和SVM構建的識別模型進行對比，驗證了該方法的優(yōu)越性.

1 航空發(fā)動機整機建模方法

1.1 軸承力模型

由軸承外圈、軸承內(nèi)圈和數(shù)個滾珠組成的滾動軸承模型如圖1所示，其中，O1為軸承外圈圓心，O2為軸承內(nèi)圈圓心，r為軸承滾道內(nèi)徑，R為滾道外徑，Nb為滾珠個數(shù)，ωb為軸承保持架的轉速.

對于本文所研究的滾動軸承，其滾珠在滾道內(nèi)等距排列，在滾道內(nèi)做純滾動運動，軸承內(nèi)外圈與滾珠的接觸符合非線性赫茲接觸理論.設軸承游隙為δ0，對于第j個滾珠，只有當其與滾道的法向接觸變形量，δj=xcosθj+ysinθj-δ0>0時，該滾珠與滾道之間才會產(chǎn)生作用力Fj，基于非線性赫茲接觸理論，軸承產(chǎn)生的軸承力為[11]

圖1 滾動軸承示意圖

(1)

(2)

式中，F(xiàn)j=Cbδj3/2H(δj)，Cb為赫茲接觸剛度，H為亥維賽函數(shù)；θj為第j個滾珠轉過的角度，θj=ωbt+2π(j-1)/Nb，ωb=ωrr/(R+r)，ωr為轉子轉速.滾動軸承的變剛度振動(VC)頻率為轉子旋轉頻率的BN倍，即ωVC=ωrBN=ωrNbr/(R+r)[12].

1.2 油膜力模型

圖2為擠壓油膜阻尼器動力學模型.其中，Os為阻尼器內(nèi)環(huán)中心，xs與ys分別為Os的橫坐標和縱坐標，Oc為阻尼器外環(huán)中心，xc與yc分別為Oc的橫坐標和縱坐標，Rs與Rc分別為阻尼器內(nèi)外環(huán)半徑，F(xiàn)t與Fr分別對應油膜的切向與徑向作用力.

圖2 擠壓油膜阻尼器動力學模型

基于短軸承半油膜模型，可知切向油膜力與徑向油膜力分別為[13]

(3)

(4)

其中

(5)

(6)

式中，μ為滑油黏度；Ro為油膜半徑，Ro=(Rs+Rc)/2；Lo為油膜長度；Co為油膜間隙，Co=Rc-Rs.

作用于水平、垂直方向的油膜力分量為

Fox=-Frcosφ+Ftsinφ

(7)

Foy=-Frsinφ-Ftcosφ

(8)

1.3 航空發(fā)動機整機模型

本文的研究對象為如圖3所示的航空發(fā)動機整機模型.圖中，klx、kly和clx、cly分別為左側鼠籠彈性支承的橫、縱向支承剛度和阻尼，krx、kry和crx、cry分別為右側鼠籠彈性支承的橫、縱向支承剛度和阻尼，cblx、cbly和cbrx、cbry分別為左右側滾動軸承橫、縱向阻尼，kflx、kfly和kfrx、kfry分別為前后安裝節(jié)橫、縱向支承剛度，cflx、cfly和cfrx、cfry分別為前后安裝節(jié)橫、縱向阻尼.Folx、Foly和Forx、Fory分別為左右側擠壓油膜阻尼器產(chǎn)生的橫、縱向油膜力，F(xiàn)blx、Fbly和Fbrx、Fbry分別為左右側滾動軸承產(chǎn)生的橫、縱向軸承力.

圖3 航空發(fā)動機整機模型

基于Timoshenko梁單元與剛性圓盤單元對轉子進行建模，將機匣視為不旋轉的軸，基于Timoshenko梁單元建立其有限元模型.各單元節(jié)點編號如圖3所示，則航空發(fā)動機整機動力學方程為

(9)

結構阻尼矩陣Cs采用瑞利阻尼的形式，即Cs=α1M+α2K，其中α1、α2為常數(shù)，可通過結構固有頻率與阻尼比求得[5].

2 航空發(fā)動機支承剛度識別方法

2.1 長短期記憶神經(jīng)網(wǎng)絡

長短期記憶神經(jīng)網(wǎng)絡是一種帶有“門”結構的循環(huán)神經(jīng)網(wǎng)絡.其在隱層中引入了能夠判斷舊有信息是否有用的單元(見圖4)，解決了處理長時間序列時易陷入梯度消失或梯度爆炸的問題，具有更強的長期預測能力[14].

圖4 LSTM的基本單元

長短期記憶神經(jīng)網(wǎng)絡基本單元中參量的計算公式如下[14]：

ft=sigmoid(Wfxxt+Wfhht-1+bf)

(10)

it=sigmoid(Wixxt+Wihht-1+bi)

(11)

gt=tanh(Wgxxt+Wghht-1+bg)

(12)

ot=sigmoid(Woxxt+Wohht-1+bo)

(13)

St=gt?it+St-1?ft

(14)

ht=τ(St)?ot

(15)

式中，ft、it、gt、ot、ht和St分別為遺忘門、輸入門、輸入節(jié)點、輸出門、中間輸出和狀態(tài)單元的狀態(tài)；Wfx、Wfh、Wix、Wih、Wgx、Wgh、Wox、Woh分別為相應門與輸入xt和中間輸出ht-1相乘的矩陣權重；bf、bi、bg、bo分別為相應門的偏置；?表示向量中元素按位相乘.

2.2 識別流程

本文利用航空發(fā)動機整機模型獲得各類樣本，然后以樣本中的位移響應為輸入、支承剛度為輸出訓練網(wǎng)絡模型，最后將實測位移響應輸入到事先訓練的網(wǎng)絡模型中，直接識別出支承剛度的實際值.該方法主要步驟如下：

①建立航空發(fā)動機整機模型，分析其結構特點，確定各處支承剛度的取值范圍.

②在目標轉速下分析支承剛度對輪盤橫向位移響應的影響，選擇影響顯著的支承剛度作為仿真變量.

③在仿真變量的取值范圍內(nèi)選取典型的剛度值，使用數(shù)值積分方法求解整機模型的不平衡響應，獲得目標轉速下對應于不同支承剛度的橫向位移響應.

④預處理橫向位移響應，將支承剛度值轉化為序列，以預處理后的位移響應為輸入，支承剛度序列為輸出構建訓練樣本與微調樣本.

⑤在支承剛度取值范圍內(nèi)隨機選取剛度值，計算對應的位移響應，并以步驟④中的方式將其構造為測試樣本.

⑥建立以LSTM為核心層的深度學習網(wǎng)絡模型，使用訓練樣本預訓練該網(wǎng)絡.

⑦設定預訓練最大容許相對誤差δpt，使用微調樣本測試已預訓練網(wǎng)絡模型，當識別相對誤差小于δpt時，使用微調樣本微調網(wǎng)絡模型，當識別相對誤差大于δpt時，重復步驟⑥.

⑧設定微調最大容許相對誤差δft，使用測試樣本測試已微調網(wǎng)絡模型，當識別相對誤差小于δft時，認為該模型即為最終的支承剛度識別模型，當識別相對誤差大于δft時，重復步驟⑦.

⑨將實驗測得的不平衡位移響應輸入識別模型，輸出序列的平均值即為目標轉速下目標支承剛度的識別值.

3 支承剛度識別算例

3.1 航空發(fā)動機整機振動響應分析

3.1.1 模型參數(shù)

航空發(fā)動機整機模型的材料參數(shù)如表1所示.用D表示梁單元直徑，L表示梁單元長度，下標數(shù)字表示節(jié)點編號，其幾何尺寸如下：D1-2=D2-3=D4-5=10 mm，D3-4=12 mm，L1-2=L6-7=25 mm，L2-3=L7-8=165 mm，L3-4=L8-9=190 mm，L4-5=L9-10=180 mm，機匣內(nèi)徑為100 mm，外徑為110 mm，轉盤直徑為80 mm，轉盤厚度為10 mm.

整機模型支承參數(shù)如下：clx、cly、crx、cry、cblx、cbly、cbrx、cbry設為20 kN/(m·s-1)，cflx、cfly、cfrx、cfry設為2 kN/(m·s-1).滾動軸承參數(shù)如表2所示，擠壓油膜阻尼器參數(shù)如表3所示.

對于結構阻尼矩陣Cs，取第1階阻尼比為 0.02，第2階阻尼比為0.04，求得轉子結構阻尼矩陣的α1為5.62，α2為7.1×10-5，機匣結構阻尼矩陣的α1為64.48，α2為3.6×10-6.

表1 部件材料參數(shù)

表2 滾動軸承參數(shù)

表3 擠壓油膜阻尼器參數(shù)

3.1.2 模型驗證

軸承力和油膜力作為支承結構中的非線性因素，其準確性對整機振動特性影響較大，因而本文將分別對軸承力模型和油膜力模型進行驗證.為了減小轉子與機匣之間的耦合振動對轉子非線性響應的影響，在驗證過程中將klx、kly、krx、kry設為15 MN/m，kflx、kfrx設為5 GN/m，kfly、kfry設為7 GN/m.

本文通過對比分析滾動軸承的VC振動規(guī)律驗證軸承力模型的有效性[12].圖5(a)給出了轉速為100 r/min、軸承游隙為20 μm時，基于整機模型計算獲得的節(jié)點3橫向位移響應；圖5(b)為橫向位移響應對應的幅值頻譜.

(a) 時域波形

(b) 幅值頻譜

由圖5(a)可知，轉子運動呈周期性，并表現(xiàn)出滾珠在滾道中運動的過程；由圖5(b)可知，轉子振動表現(xiàn)出VC頻率(滾珠通過頻率)及其諧波.本文計算結果與文獻[12]的計算結果相符，因而驗證了本文軸承力模型的有效性.

本文通過驗證擠壓油膜阻尼器的減振作用來驗證油膜力模型的有效性[13].圖6為不同轉速下轉子軸頸處(節(jié)點11)的穩(wěn)態(tài)位移響應振幅.

圖6 振動幅值對比

由圖6可知，擠壓油膜阻尼器在臨界轉速附近具有明顯的減振作用，油膜力對轉子振動的抑制使得11號節(jié)點的振幅一直小于油膜間隙Co，擠壓油膜阻尼器的減振作用得以驗證，從而驗證了本文油膜力模型的有效性.

3.1.3 變量選擇

考慮到支承剛度與位移響應之間的非線性映射關系可能并不單調，本文采用文獻[6]中的變量選擇方法選取仿真變量.

根據(jù)工程經(jīng)驗，klx=kly，krx=kry，考慮其取值范圍為10～20 MN/m；kflx、kfrx取值范圍為45～55 MN/m；安裝節(jié)縱向剛度很強，不考慮其對位移響應的影響，將kfly、kfry設為70 MN/m.

本文的目標轉速為3 000 r/min，通過比較該轉速下各剛度對節(jié)點3橫向位移響應的影響，選擇影響顯著的支承剛度作為仿真變量.

分別設klx、kly、krx、kry為15 MN/m，kflx、kfrx為50 MN/m，單獨改變其中一個支承剛度，求解對應的位移響應.以一段位移響應的信號能量為比較對象，信號能量等于數(shù)字信號序列所有元素的平方和，分析各剛度對信號能量的影響，影響規(guī)律如圖7所示.

由圖7(c)和(d)可知，鼠籠彈性支承的剛度變化對信號能量的影響最為顯著.因此，將兩側鼠籠彈性支承的剛度klx、kly和krx、kry作為仿真變量，并將kflx、kfrx設為50 MN/m.

3.1.4 樣本構造

在鼠籠彈性支承的取值范圍10～20 MN/m內(nèi)每隔1 MN/m取一個剛度值，則每個支承剛度都有11種取值，將序號為奇數(shù)的剛度取值進行組合，獲得36種不同的剛度組合以構造訓練樣本.將序號為偶數(shù)的5種取值組合成另外25種剛度組合，用來構造微調樣本.為了驗證網(wǎng)絡模型的準確性，在剛度取值范圍內(nèi)隨機生成5種剛度值，用來構造25個測試樣本.

(a) 前安裝節(jié)橫向剛度kflx

(c) 左側鼠籠彈性支承剛度klx、kly

基于Newmark-β法獲得這些剛度組合對應的位移響應，以節(jié)點1與節(jié)點3的橫向位移響應為原始數(shù)據(jù)，基于z-score方法將其標準化.同時將剛度值轉化為與位移響應長度相同的等值序列，并將其縮小107倍以統(tǒng)一輸入與輸出的量級.最后以位移響應為輸入數(shù)據(jù)，對應的剛度值序列為輸出數(shù)據(jù)構造各個樣本.

3.2 深度學習網(wǎng)絡建模

建立如圖8所示的深度學習網(wǎng)絡模型.其中，輸入層節(jié)點數(shù)為2，長短期記憶網(wǎng)絡層隱含層單元數(shù)為256，全連接層節(jié)點數(shù)分別為128和2，輸出層節(jié)點數(shù)為2，丟棄層歸零概率為50%，網(wǎng)絡模型的損失函數(shù)為均方誤差函數(shù).

圖8 深度學習網(wǎng)絡模型框架

取δpt=10%，設初始學習率為0.01，基于Adam優(yōu)化算法預訓練網(wǎng)絡，數(shù)次預訓練后，微調樣本的支承剛度識別相對誤差分別降為5.28%與 3.11%，開始對網(wǎng)絡模型進行微調.

取δft=1%，設初始學習率為1.0×10-4，使用微調樣本對網(wǎng)絡模型進行全局權重微調.微調網(wǎng)絡權重后，25個測試樣本的識別相對誤差絕對值的均值已小于1%，將此時的深度學習網(wǎng)絡作為支承剛度識別模型.

3.3 方法驗證

本文為位移響應附加信噪比為5%的高斯白噪聲模擬試驗數(shù)據(jù)，為了減小噪聲隨機性的影響，基于25個仿真樣本構造了500個含噪聲的新樣本.為了驗證本文方法的優(yōu)越性，在其他條件不變的情況下，分別基于RBFNN與SVM構建支承剛度識別模型，并對含噪聲樣本進行識別，3種網(wǎng)絡的識別結果如表4所示.

表4 3種網(wǎng)絡識別相對誤差絕對值的平均值 %

由表4可知：無噪聲時，3種網(wǎng)絡的識別精度較高，LSTM具有最高的識別精度；5%噪聲時，3種網(wǎng)絡均受到一定影響，但LSTM受到的影響最小，仍保持較高的識別精度.

以本文方法識別相對誤差最大的樣本為分析對象，其剛度識別值對應的位移響應與剛度真實值對應的位移響應如圖9所示，圖中實、虛線分別表示剛度真實值、剛度識別值對應的3號節(jié)點橫向位移響應.由圖可知，剛度識別值對應的位移響應與真實值對應的位移響應高度一致，利用剛度識別后的航空發(fā)動機整機模型能夠獲得更加準確的位移響應，進一步驗證了本文方法的有效性.

圖9 橫向位移響應對比

4 結論

1) 建立了包含擠壓油膜阻尼器的航空發(fā)動機轉子-滾動軸承-支承-機匣耦合動力學模型，通過仿真研究了滾動軸承的VC振動規(guī)律與擠壓油膜阻尼器的減振作用，并與現(xiàn)有文獻進行對比，驗證了本文模型的有效性.

2) 以航空發(fā)動機整機模型為識別對象，基于LSTM構建了支承剛度識別模型，該模型在無噪聲情況下的識別誤差小于1%，在5%噪聲情況下的識別誤差小于2%.

3) 相比于RBFNN與SVM，本文構建的深度學習網(wǎng)絡具有更高的識別精度，表明了本文方法在解決支承剛度識別問題時的優(yōu)越性.

4) 本文方法利用神經(jīng)網(wǎng)絡的泛化特性對支承剛度進行識別，無需迭代計算，避開了動力學反問題中復雜的尋優(yōu)過程，大大減少了單次運算的計算量，降低了計算機的運算負荷.