基于真實數(shù)據(jù)反演的風扇轉(zhuǎn)子本機平衡方法

陳立芳，孫亞冰，*，周書華，高強，喬保棟，李棟

1.北京化工大學 發(fā)動機健康監(jiān)控及網(wǎng)絡(luò)化教育部重點實驗室，北京 100029

2.北京化工大學 高端壓縮機及系統(tǒng)技術(shù)全國重點實驗室，北京 100029

3.中國航發(fā)沈陽發(fā)動機研究所，沈陽 110015

風扇轉(zhuǎn)子往往工藝平衡狀態(tài)良好，但工作狀態(tài)時不平衡振動依然較大，需在發(fā)動機試車臺上進行本機平衡［1-5］。當前，渦扇發(fā)動機風扇轉(zhuǎn)子本機平衡多采用三圓法［6-9］進行，通過3 次試重獲得不平衡量，至少5 次啟停車完成風扇轉(zhuǎn)子本機平衡，雖抑振效果良好，但啟停次數(shù)過多使得動平衡效率低且成本高。

此外，航空發(fā)動機研究院所也經(jīng)常采用影響系數(shù)法進行風扇轉(zhuǎn)子本機平衡。姜廣義等［10］采用葉片檢測系統(tǒng)測量定制磁鋼的鍵向信號，一次試重即獲得了配平大小和方位，實現(xiàn)了低速（3 000 r/min 左右）下的風扇轉(zhuǎn)子本機平衡。然而實際工作中，機匣振動經(jīng)復雜路徑傳遞會導致影響系數(shù)穩(wěn)定性差，往往需要增加試重次數(shù)。而且，發(fā)動機各機臺間動力學特性差異較大，測得的影響系數(shù)難以推廣應(yīng)用。

虛擬動平衡法通過仿真構(gòu)建與實際轉(zhuǎn)子和參數(shù)相符的動力學模型，在模型上添加虛擬激勵獲得響應(yīng)，計算影響系數(shù)進行動平衡的方法。目前，該方法基于鍵向信號可以實現(xiàn)無試重虛擬動平衡［11-16］。賓光富等［17］在汽輪機模擬試驗臺上進行過有效驗證，王維民等［18］在離心式氨氣壓縮機上也有成功應(yīng)用。結(jié)合粒子群尋優(yōu)算法，陳立芳等［19］提出一種基于無鍵相虛擬動平衡方法，在有限元模型基礎(chǔ)上利用粒子群算法尋優(yōu)轉(zhuǎn)子不平衡量，通過最多2 次試重完成動平衡，試驗臺測試抑振效果顯著。但以上方法均對模型的精度要求高，難以應(yīng)用在具有薄壁機匣復雜結(jié)構(gòu)的航空發(fā)動機風扇轉(zhuǎn)子上。

因此，針對航空發(fā)動機風扇轉(zhuǎn)子本機平衡的重大需求，本文結(jié)合反演理論和虛擬動平衡方法，提出一種啟停次數(shù)少，無需鍵相信號、能匹配不同機臺動力學特性的本機平衡方法。本方法融合轉(zhuǎn)子動力學和數(shù)據(jù)反演理論，通過分析風扇轉(zhuǎn)子振動傳遞路徑明確動剛度特征參數(shù)，基于模型和動平衡數(shù)據(jù)提取機臺的動力學特征參數(shù)信息，建立對應(yīng)機臺反演特征參數(shù)庫用于計算轉(zhuǎn)子不平衡量。本文的創(chuàng)新之處在于，有效利用發(fā)動機現(xiàn)場平衡數(shù)據(jù)，提高動平衡效率。

1 反演理論和風扇轉(zhuǎn)子反演模型

1.1 反演理論

反演模擬原理［20-21］是將一個物理系統(tǒng)抽象為一個數(shù)學上的參數(shù)化模型，利用可測參數(shù)的某些實測結(jié)果推斷模型參數(shù)的真值。系統(tǒng)被參數(shù)化后，模型可用一組參數(shù)來表達，模型參數(shù)所構(gòu)成的空間稱為模型空間M={mi}，可直接觀測數(shù)據(jù)構(gòu)成數(shù)據(jù)空間D={di}。在實際反演問題中，根據(jù)D中的觀測數(shù)據(jù)反推M中的待求參數(shù)［22］。如果系統(tǒng)為一線性系統(tǒng)，則相應(yīng)模型可用線性矩陣方程表達，可將其表示為

通過廣義求逆可得模型不可測參數(shù)的解，將其表示為

式中：D為模型可測參數(shù)；M為模型不可測參數(shù)；G為相關(guān)系數(shù)矩陣；G-1為廣義逆矩陣。

慮及在風扇轉(zhuǎn)子動力學問題中，無法直接測量轉(zhuǎn)子的不平衡量U，而能直接觀測的是由不平衡量U產(chǎn)生的振動響應(yīng)X。因此，可根據(jù)反演理論，以多轉(zhuǎn)速下的轉(zhuǎn)子不平衡量U作為模型中的待求參數(shù)，振動響應(yīng)X作為可觀測數(shù)據(jù)，構(gòu)建風扇轉(zhuǎn)子振動系統(tǒng)的參數(shù)化模型，可為其他相似機臺的轉(zhuǎn)子不平衡量求解提供技術(shù)參數(shù)支撐。

1.2 構(gòu)建風扇轉(zhuǎn)子反演模型

由于風扇轉(zhuǎn)子結(jié)構(gòu)復雜，振動測點通常布置在進氣機匣表面。根據(jù)轉(zhuǎn)子動力學理論將風扇轉(zhuǎn)子振動傳遞系統(tǒng)描述為轉(zhuǎn)子-軸承-機匣系統(tǒng)，發(fā)生在風扇轉(zhuǎn)子上的不平衡振動經(jīng)軸承（支點1）傳遞到機匣，振動傳遞系統(tǒng)及路徑分解如圖1所示。

圖1 風扇轉(zhuǎn)子振動傳遞系統(tǒng)及路徑分解Fig.1 Fan rotor vibration transmission system and path decomposition

為獲得反演特征參數(shù)，根據(jù)激勵與響應(yīng)之間的相互關(guān)系，分解振動傳遞路徑，提出風扇轉(zhuǎn)子振動傳遞過程中的4 個動剛度參數(shù)：

式中：Kcs為機匣到軸承支點1 的動剛度；Fs為支點1 的力；Xc為機匣的上測點的振動位移；Ksd為支點1 到風扇轉(zhuǎn)子動平衡平面的動剛度；Fd為轉(zhuǎn)子上的不平衡力；Xs為支點1 的位移；Kss為支點1到支點1 的動剛度；Kcd為機匣到風扇轉(zhuǎn)子動平衡平面的動剛度。

根據(jù)“動平衡平面（d）-支點1（s）-機匣（c）”的振動傳遞關(guān)系可得：

上述動剛度參數(shù)結(jié)合數(shù)據(jù)反演原理可構(gòu)建風扇轉(zhuǎn)子振動傳遞系統(tǒng)的參數(shù)化模型如圖2所示。以多轉(zhuǎn)速下的不平衡量U={Ui}為待求模型參數(shù)，機匣振動位移X={Xi}為可觀測量。動剛度Kcd、Ksd、Kss和Kcs為可計算模型特征參數(shù)。

圖2 風扇轉(zhuǎn)子振動系統(tǒng)參數(shù)化模型Fig.2 Parametric model of fan rotor vibration system

2 真實數(shù)據(jù)反演的本機平衡方法

2.1 建立風扇轉(zhuǎn)子有限元模型

根據(jù)某型渦扇發(fā)動機風扇轉(zhuǎn)子真實結(jié)構(gòu)建立有限元模型。建模原則：一、二、三級葉片用等質(zhì)量和等轉(zhuǎn)動慣量的剛性盤單元表示；忽略倒角、螺紋孔等結(jié)構(gòu)特征；支點1 和支點2 用軸承單元表示，假設(shè)支承剛度各向同性，無交叉剛度和阻尼［23］；集中質(zhì)量作用處、軸徑變化處、軸承支承處均確保有單元節(jié)點，如圖3所示：軸單元238個，盤單元3 個，軸承單元2 個，共243 個單元。其中支點1 位置為8 號節(jié)點，支點2 位置為102 號節(jié)點，試重位置為11 號節(jié)點。

圖3 風扇轉(zhuǎn)子有限元模型Fig.3 Finite element model of fan rotor

2.2 基于實測臨界轉(zhuǎn)速優(yōu)化模型

臨界轉(zhuǎn)速是轉(zhuǎn)子動力學特性的重要表征參數(shù)，其中支承剛度對臨界轉(zhuǎn)速的影響最大［24-27］，利用機臺的實測臨界轉(zhuǎn)速對支承剛度進行優(yōu)化以提高風扇轉(zhuǎn)子有限元模型與實際機臺的匹配度。

首先，確定模型優(yōu)化的目標參數(shù)?；谀Ｐ瓦M行臨界轉(zhuǎn)速分析，探究支點1 和支點2 支承剛度對臨界轉(zhuǎn)速的影響規(guī)律?？紤]到臨界轉(zhuǎn)速裕度需大于20%，因此，以支點1支承剛度2.2×107N/m，支點2 支承剛度1.0×108N/m 為基準值，支承剛度變化范圍為正負基準值的1/2，每次固定一個支點剛度，展開前三階臨界轉(zhuǎn)速分析，如圖4所示。

圖4 兩支點支承剛度對臨界轉(zhuǎn)速的影響Fig.4 Influence of two support stiffness on critical speed

由圖 4（a）可見，支點1 支承剛度可有效調(diào)節(jié)一階臨界轉(zhuǎn)速，對二階臨界轉(zhuǎn)速影響不大，對三階臨界轉(zhuǎn)速有一定的調(diào)節(jié)效果；由圖 4（b）可見，支點2 支承剛度對一階臨界轉(zhuǎn)速影響不大，對二階臨界轉(zhuǎn)速幾乎無影響，可有效調(diào)節(jié)三階臨界轉(zhuǎn)速。對比可知，支點1 的支承剛度對一階臨界轉(zhuǎn)速影響較大，而支點2 對一階臨界轉(zhuǎn)速幾乎無影響。因此，可通過優(yōu)化支點1 支承剛度來匹配機臺的一階臨界轉(zhuǎn)速。

其次，基于粒子群算法進行支點1 剛度優(yōu)化。以某機臺實測數(shù)據(jù)為例，風扇轉(zhuǎn)子振動響應(yīng)隨轉(zhuǎn)速變化如圖5所示，可知該機臺的一階臨界轉(zhuǎn)速N1為6 200 r/min 左右，按照經(jīng)驗值固定支點2 的支承剛度為1.0×108N/m，通過粒子群尋優(yōu)算法得到支點1 的剛度。

圖5 實測某機臺振幅隨轉(zhuǎn)速變化曲線Fig.5 Variation curve of amplitude with speed for machine

利用自主編寫的基于Python 的可快速調(diào)用的動力學軟件和風扇轉(zhuǎn)子有限元模型迭代計算一階臨界轉(zhuǎn)速。根據(jù)粒子群算法原理，設(shè)計臨界轉(zhuǎn)速為f(K)，將轉(zhuǎn)子支點1 的支承剛度K作為粒子的位置，將實測與設(shè)計臨界轉(zhuǎn)速的差F(K)的絕對值作為粒子群的適應(yīng)度函數(shù)：

計算流程和結(jié)果如圖6 和圖 7所示。由圖 7可見，利用粒子群優(yōu)化算法最終得出K1為2.3×107N/m，將其添加到有限元模型中進行計算得一階臨界轉(zhuǎn)速為6 203 r/min，如圖8所示，與實測臨界轉(zhuǎn)速誤差僅為0.05%，表明有限元模型得到優(yōu)化，且匹配實際機臺特性。

圖6 粒子群算法優(yōu)化K1Fig.6 Particle swarm optimization K1

圖7 支點1 剛度優(yōu)化Fig.7 Stiffness optimization of fulcrum 1

圖8 模型優(yōu)化后臨界轉(zhuǎn)速Fig.8 Critical speed after model optimization

2.3 基于模型和動平衡數(shù)據(jù)反演特征參數(shù)

2.3.1 計算Kss和Ksd

基于風扇轉(zhuǎn)子優(yōu)化后的模型計算Kss和Ksd，流程如圖9所示。

圖9 Kss 和Ksd 計算流程Fig.9 Kss and Ksd calculation process

2.3.2 計算Kcs

在試驗臺架上采用錘擊法可以測試Kcs，但對某機臺進行動剛度Kcs測試時，發(fā)現(xiàn)前后2 次錘擊測試的動剛度Kcs誤差較大，如圖10所示，會導致不平衡量計算錯誤。

圖10 錘擊測試所得動剛度KcsFig.10 Dynamic stiffness Kcs from hammer test

當前，發(fā)動機機臺具有大量動平衡數(shù)據(jù)，而動平衡數(shù)據(jù)中包含著轉(zhuǎn)子動力學信息，未能有效利用。因此，本文提出一種基于動平衡數(shù)據(jù)反演參數(shù)Kcs的方法。由圖1 的振動傳遞路徑分析可知，機匣到動平衡平面的動剛度關(guān)系如圖11所示。通過反演理論提取動平衡數(shù)據(jù)中的相關(guān)動力學信息，由于轉(zhuǎn)子的不平衡力主要由不平衡量產(chǎn)生的離心力所致，因此可得到式（6），進而得到Kcd；結(jié)合式（4）即可得到Kcs，如式（7）所示。

圖11 振動傳遞路徑中的動剛度Fig.11 Dynamic stiffness in vibration transmission path

通過該方法可優(yōu)化出不同轉(zhuǎn)速下的Kcs，圖12所示是同一機臺通過2 種動平衡方法所得Kcs，對比錘擊測試所得Kcs可知，經(jīng)過2 種動平衡數(shù)據(jù)優(yōu)化得到的動剛度較為接近，具有較高的精確度，證實了利用動平衡數(shù)據(jù)優(yōu)化出的動剛度Kcs更加準確，可用于風扇轉(zhuǎn)子本機平衡。

圖12 不同方法計算的Kcs 參數(shù)Fig.12 Kcs parameters calculated by different methods

2.4 基于模型和特征參數(shù)反演不平衡量

確定風扇轉(zhuǎn)子反演模型的特征參數(shù)后，由式（6）和式（7）可建立多轉(zhuǎn)速下的轉(zhuǎn)子不平衡量、模型特征參數(shù)和機匣響應(yīng)之間的線性矩陣方程：

通過對式（8）進行廣義求逆可得：

根據(jù)式（8）和式（9）只能求得不平衡量的大小，并不能確定不平衡量的位置。為此選擇圖1所示風扇轉(zhuǎn)子前支點軸承內(nèi)部改裝后的分油盤作為試重平面，建立試重與原始不平衡量大小和位置的關(guān)系，如圖13所示。圖13 中，U0為轉(zhuǎn)子原始不平衡量大??；U1為試重大??；U01為添加試重后合成不平衡量大小；α為試重角度；β為不平衡量的位置；γ為試重與不平衡量的夾角。U0和U01均可通過式（9）算出，則夾角γ為

圖13 不平衡量位置識別Fig.13 Unbalance location identification

由于反余弦函數(shù)的對稱特性，在[-π，π]范圍內(nèi)的同一函數(shù)值會對應(yīng)2 個自變量，且兩者互為相反數(shù)，因此可得到不平衡量的2 個位置為

根據(jù)式（10）和式（11）可以得出不平衡量的2 個位置，若再進行一次試重，則可以得到不平衡量的唯一位置。至此，基于真實數(shù)據(jù)反演的風扇轉(zhuǎn)子本機平衡方法原理如圖14所示。

圖14 基于真實數(shù)據(jù)反演的風扇轉(zhuǎn)子本機平衡方法Fig.14 Principle of fan rotor local balancing method based on real data inversion

3 本機平衡現(xiàn)場試驗

3.1 建立機臺特征參數(shù)庫

不同機臺發(fā)動機動力學特性差異較大，需分別建立各自的特征參數(shù)庫。以具有接近臨界轉(zhuǎn)速（6 200 r/min 左右）的2 機臺為例，基于2.1 節(jié)～2.3節(jié)的方法建立特征參數(shù)庫，如表1和表2所示。

表1 機臺1 動剛度KcsTable 1 Machine 1 dynamic stiffness Kcs

表2 機臺2 的動剛度KcsTable 2 Machine 2 dynamic stiffness Kcs

以機臺2 的轉(zhuǎn)速為參考，機臺1 內(nèi)對應(yīng)參考轉(zhuǎn)速的Kcs通過插值法得出，對比2 機臺的動剛度Kcs，如圖15所示。由圖15 可知，在臨界轉(zhuǎn)速6 200 r/min 前后，機匣測點至支點1 的動剛度Kcs呈先減小后增大的趨勢；2機臺的動剛度Kcs在一階臨界轉(zhuǎn)速之前差異不大，臨界轉(zhuǎn)速之后差異增大。因此在對不同的機臺進行動平衡時，需建立各自的反演模型參數(shù)庫。

圖15 機臺1 和機臺2 的動剛度對比Fig.15 Comparison of dynamic stiffness of machine 1 and machine 2

3.2 現(xiàn)場平衡結(jié)果分析

應(yīng)用本方法對上述2 機臺進行現(xiàn)場動平衡，均獲得了較好的抑振效果，在全轉(zhuǎn)速范圍內(nèi)實現(xiàn)了100%的振動抑制。其中機臺1 降振達40%以上，機臺2 降振達43%以上。現(xiàn)以機臺2 為例闡述動平衡實施過程。首先，在試車臺架上拆下發(fā)動機進氣帽罩等相關(guān)組件，選擇圖1所示風扇轉(zhuǎn)子前支點軸承內(nèi)部改裝過的分油盤作為試重平面，配平時再折算到配重平面。選擇試重為1 123 g·mm∠-5°，實測機臺2 的初始振動及試重后振動如表3所示。

表3 初始及試重振動幅值Table 3 Initial and trial recovery amplitude

根據(jù)式（9）～式（11）計算可得各轉(zhuǎn)速下的不平衡量的大小及位置如表4所示。因現(xiàn)場實際記錄缺少部分振動響應(yīng)值，故未能計算出相應(yīng)配重位置。

表4 不平衡量的大小及位置Table 4 Size and location of unbalance

用表4計算結(jié)果進行動平衡，從表中可以看出，6 600 r/min 之后的不平衡量均相差不大，但綜合考慮到不平衡量與配重角度2方面因素，不難看出6 600 與7 600 r/min 2 個轉(zhuǎn)速下配重位置最為接近，因此取兩者平均值1 654.5 g·mm 作為配重大?。粎⒖棘F(xiàn)場試重角度210°的抑振效果，選取配重角度237.5°（排除112.5°），即：1 654.5 g·mm ∠237.5°。由于受平衡螺釘重量及位置限制，最終實際添加的矯正不平衡量為1 600 g·mm∠235°，配平前后的振幅變化如圖16所示。

圖16 平衡前后振幅變化Fig.16 Amplitude changes before and after equilibrium

根據(jù)配重結(jié)果可知在轉(zhuǎn)速范圍內(nèi)全部實現(xiàn)了降振，降振達43%～68%，如圖17所示。

圖17 振動降幅百分比Fig.17 Percentage of vibration decline

4 結(jié) 論

針對發(fā)動機風扇轉(zhuǎn)子現(xiàn)場動平衡效率低，且大量動平衡數(shù)據(jù)未能有效利用的問題，本文提出一種基于模型和現(xiàn)場真實動平衡數(shù)據(jù)反演的風扇轉(zhuǎn)子本機平衡方法。通過動平衡數(shù)據(jù)提取機臺的動力學特征參數(shù)信息并建立機臺的反演特征參數(shù)庫，用于下一次動平衡計算。通過機臺現(xiàn)場動平衡驗證了該方法的有效性。

1） 風扇轉(zhuǎn)子有限元模型的精確度是此動平衡方法的關(guān)鍵，所用模型的動力學特性必須與真實機臺保持高度的一致。

2） 利用錘擊法所測動剛度Kcs誤差較大，通過機臺的真實歷史動平衡數(shù)據(jù)優(yōu)化動剛度Kcs，可提高其精確度，可為風扇轉(zhuǎn)子本機平衡提供較高匹配度的機臺參數(shù)。

3） 本方法無需鍵相，最少試重1 次，至多試重2 次，減少了啟停機次數(shù)，且對風扇轉(zhuǎn)子能起到較好的抑振效果。

4） 風扇轉(zhuǎn)子作為發(fā)動機低壓轉(zhuǎn)子的前端，降低風扇轉(zhuǎn)子的振動對低壓渦輪轉(zhuǎn)子的振動也具有一定的改善。

本方法通過建立機臺反演特征參數(shù)庫，提高動平衡效率，與機臺匹配度越高動平衡效果越好。因此需要大量已有機臺的動平衡數(shù)據(jù)優(yōu)化風扇轉(zhuǎn)子的動剛度參數(shù)，這也造成本方法應(yīng)用的局限性。