趙洪利， 高晶東

(中國民航大學航空工程學院，天津 300300)

航空發(fā)動機部件特性圖(Map圖)，也被形象地稱為“萬有特性曲線”，即可以從該特性圖上讀取發(fā)動機各個單元體全包線范圍內(nèi)的換算流量、效率、壓比、換算轉(zhuǎn)速等性能參數(shù)[1-2]。通常由發(fā)動機制造廠商通過大量的試車測試繪制而成，成本高昂，一般不會向外界提供。另外，在進行發(fā)動機建模仿真過程中，發(fā)動機部件特性圖的準確與否直接影響所建立模型的精度[3-7]。Kong等[3-5]基于系統(tǒng)識別方法和遺傳算法完成發(fā)動機部件特性圖修正；Li等[6-7]使用遺傳算法結(jié)合最小二乘法等進行耦合因子函數(shù)系數(shù)項尋優(yōu)，修正發(fā)動機部件特性圖；劉盾[8]、白磊等[9]利用非線性最小二乘擬合結(jié)合模型辨識等方法進行了發(fā)動機部件特性圖耦合修正；肖洪等[10]、楊天南等[11]使用單純形和遺傳算法建立自適應模型進行了發(fā)動機部件特性模擬。但是研究內(nèi)容大都進行了發(fā)動機部件特性設計點匹配研究，對于非設計點匹配的論述較少。而且，使用的模型識別和遺傳算法等計算量大，迭代效率較低。

基于發(fā)動機試車臺數(shù)據(jù)，利用發(fā)動機通用部件特性進行發(fā)動機部件特性修正，可以獲取接近發(fā)動機實際工況的部件特性圖，進而準確地完成面向?qū)ο笮桶l(fā)動機的建模仿真，實現(xiàn)發(fā)動機各個不同工作狀態(tài)下的準確性能計算，對于實現(xiàn)準確的發(fā)動機狀態(tài)監(jiān)控與故障診斷意義重大。

1 發(fā)動機試車臺數(shù)據(jù)

以某型號渦扇發(fā)動機試車臺數(shù)據(jù)為支撐進行其部件特性圖修正，根據(jù)試車臺數(shù)據(jù)與發(fā)動機原理等相關(guān)知識，利用試車臺數(shù)據(jù)作為已知條件，基于MATLAB/Simulink自主建立發(fā)動機穩(wěn)態(tài)模型，進行逆向工程求取發(fā)動機設計點及其他各個穩(wěn)態(tài)工作點的性能參數(shù)[12-13]。

結(jié)合發(fā)動機原理相關(guān)知識與試車臺數(shù)據(jù)，對該型渦扇發(fā)動機進行站位劃分，具體站位示意圖如圖1所示，各個站位所代表的截面含義在表1中給出了相應說明。

表1 某型渦扇發(fā)動機站位說明

圖1 某型渦扇發(fā)動機站位示意圖Fig.1 Turbofan engine station indication

對發(fā)動機試車臺數(shù)據(jù)進行整理，如表2所示，作為發(fā)動機穩(wěn)態(tài)模型仿真與工作點計算的已知條件。

表2 部分該型渦扇發(fā)動機試車臺數(shù)據(jù)

2 牛頓迭代法設計

航空發(fā)動機模型是一種典型的多元非線性模型，因此，其共同工作方程求解則轉(zhuǎn)化為多元非線性方程的求解問題[14-15]。牛頓迭代法(N-R iteration)是一種有效的非線性方程求解方法。利用牛頓迭代法進行所建立發(fā)動機模型的迭代求解，在傳統(tǒng)牛頓迭代方法基礎上，增加變化迭代步長的設計，以增加方程迭代求解的速度，提高模型收斂性。

在進行發(fā)動機狀態(tài)參數(shù)計算過程中，未知參數(shù)為：x=[πfan,πl(wèi)pc,πhpc,πhpt,πl(wèi)pt,n1,n2]T，在進行牛頓迭代法求解過程中，依據(jù)發(fā)動機推力等級等已知狀態(tài)參數(shù)合理給出上面七個未知參數(shù)的初猜值，進行非線性共同工作方程求解：

φi(πfan,πl(wèi)pc,πhpc,πhpt,πl(wèi)pt,n1,n2)=ei

(1)

式(1)中：πfan、πl(wèi)pc、πhpc、πhpt、πl(wèi)pt分別表示發(fā)動機依次表示風扇、低壓壓氣機、高壓壓氣機、高壓渦輪、低壓渦輪的增壓比或落壓比；n1和n2分別為低壓轉(zhuǎn)子相對轉(zhuǎn)速和高壓轉(zhuǎn)子相對轉(zhuǎn)速，等于實測轉(zhuǎn)速值除以最大轉(zhuǎn)速值。

利用牛頓迭代法求解非線性方程組[式(1)]，迭代式為

x(k+1)=xk-J-1φ(xk), k=0,1,2，…

(2)

式(2)中：xk為上一步初猜值的解向量；x(k+1)為迭代一次后的修正解向量；J為非線性方程組的雅可比矩陣，其表達式為

(3)

選取歐幾里得范數(shù)作為解的誤差檢驗目標函數(shù)，即：

(4)

發(fā)動機穩(wěn)態(tài)模型迭代計算過程如圖2所示，模型初始化時，使用賦予模型的初猜值進行模型計算，然后計算模型殘差。當模型殘差ei達到要求的精度eps1時，變換擾動步長，加快迭代收斂進程；當模型殘差ei達到要求的精度eps2時，輸出最優(yōu)解向量。

圖2 模型穩(wěn)態(tài)工作點迭代過程Fig.2 Iteration process of engine steady working point

3 部件特性圖修正

通常，發(fā)動機部件特性圖的獲取方法分為兩種：一種方法是擁有大量的某型號發(fā)動機全包線范圍內(nèi)的試車數(shù)據(jù)，使用實際的試車數(shù)據(jù)擬合出各轉(zhuǎn)速下多個匹配的工作點，并通過多項式擬合得到各個工作狀態(tài)下轉(zhuǎn)速線；另一種是擁有的試車數(shù)據(jù)較少，通過選取合適的發(fā)動機通用部件特性圖，利用耦合法(scaling method)進行發(fā)動機設計點處參數(shù)修正，在此基礎上，將其他非設計點處的試車臺數(shù)據(jù)與發(fā)動機穩(wěn)態(tài)模型計算所獲得的參數(shù)進行多點匹配，通過逐段修正的方法獲取修正后的發(fā)動機部件特性圖[16-17]。由于獲取大量全包線工作范圍內(nèi)的試車數(shù)據(jù)可操作性差，不確定因素較多，且成本高昂，因此，采用部件特性圖的多點匹配的方法，基于多穩(wěn)態(tài)工作點的試車臺數(shù)據(jù)進行部件特性圖逐段修正，以此提高部件特性圖的準確性，進而建立準確性較高的發(fā)動機模型。

3.1 部件特性圖設計點耦合

首先，基于耦合法，進行發(fā)動機設計點處的參數(shù)匹配，使得發(fā)動機設計點處模型殘差ei達到預定的精度要求，進行發(fā)動機設計點匹配時，修正部件特性圖的耦合因子求解公式如下所示：

(5)

(6)

(7)

式中：f表示耦合因子；π表示壓比；Wa表示流量；η表示效率；含有下標D的參數(shù)表示目標發(fā)動機參考點數(shù)據(jù)；含有下標M、D的參數(shù)表示原特性圖的發(fā)動機參考點數(shù)據(jù)。

以風扇為例，將部件特性圖按照設計點進行耦合因子修正，修正前后部件特性圖對比如圖3所示，求得耦合因子結(jié)果如表3所示。

圖3 耦合法部件特性圖修正前后對比Fig.3 Original and scaled map comparison

表3 耦合因子與縮放結(jié)果示例

3.2 部件特性圖非設計點耦合

在進行部件特性圖非設計點修正時，設計點已經(jīng)完成了與試車數(shù)據(jù)的匹配，滿足精度要求，需要保證已經(jīng)完成修正的設計點精度，因此，后續(xù)修正工作以設計點為基準，不再對設計點所在轉(zhuǎn)速線的數(shù)據(jù)進行修正變動。

根據(jù)已獲取的發(fā)動機試車臺數(shù)據(jù)分析，所測試的發(fā)動機工作點均處于高功率穩(wěn)態(tài)下的工作狀態(tài)，因此，高壓渦輪部件基本處于臨界或超臨界狀態(tài)，渦輪部件特性圖的性能參數(shù)插值結(jié)果變化較小。因此，在進行非設計點匹配修正時，選取的待優(yōu)化參數(shù)如表4所示。

表4 非設計點待優(yōu)化的耦合因子修正項

以風扇、壓氣機等部件耦合因子修正項作為待優(yōu)化參數(shù)，如表4所示；以發(fā)動機試車臺數(shù)據(jù)中各個穩(wěn)態(tài)工作點測量參數(shù)為目標參數(shù)，如表2所示，包括風扇外涵出口壓力P13、低壓壓氣機出口壓力P25、低壓壓氣機出口溫度T25、高壓壓氣機出口溫度T3、低壓渦輪出口壓力P5、低壓渦輪出口溫度T5、低壓轉(zhuǎn)子轉(zhuǎn)速N1、高壓轉(zhuǎn)子轉(zhuǎn)速N2、發(fā)動機推力FN等，進行發(fā)動機穩(wěn)態(tài)模型迭代計算，求取各部件特性的耦合因子修正項。

優(yōu)化求得的部件特性圖修正因子對非設計點轉(zhuǎn)速線參數(shù)進行修正，獲得修正后的轉(zhuǎn)速線，如圖4所示，并重新進行部件特性圖插值計算得到修正后的匹配工作點。圖4中，B點為非設計點轉(zhuǎn)速線修正前匹配到的穩(wěn)態(tài)工作點，B′點為非設計點轉(zhuǎn)速線修正后重新進行部件特性圖匹配計算得到的穩(wěn)態(tài)工作點。

圖4 非設計點風扇部件特性圖修正前、后對比Fig.4 Comparison of off-design point correction

最后，達到模型計算參數(shù)與目標參數(shù)匹配的目的。在進行穩(wěn)態(tài)模型迭代計算時，以式(8)作為誤差檢驗目標函數(shù)。

min(error)=min[f1(e)+f2(e′)]=

(8)

式(8)中：j為非設計點待優(yōu)化的耦合因子修正項個數(shù)；x為目標參數(shù)，下標cal,test分別表示模型計算結(jié)果與實際試車臺數(shù)據(jù)。

4 修正結(jié)果分析

利用耦合法進行發(fā)動機設計點參數(shù)匹配修正后，發(fā)動機部件特性圖已經(jīng)完成了一次整體調(diào)整，更接近實際的發(fā)動機部件特性，但是進行發(fā)動機非設計點參數(shù)匹配時仍存在較大誤差。因此，在保持設計點處修正參數(shù)不變的基礎上，通過牛頓迭代法迭代計算獲取發(fā)動機部件特性圖耦合因子修正項，進行部件特性圖非設計點的匹配，修正非設計狀態(tài)下各單元體的部件特性圖。模型修正前后，各站位主要測量參數(shù)與模型計算匹配結(jié)果的相對誤差如圖5所示。

從圖5中可以看出，采用多點匹配方法完成模型修正后，發(fā)動機各站位參數(shù)匹配精度均有所提高，平均建模誤差由5.163 7%降低到多點匹配后的0.529 7%，模型精度得到很大程度上的提升。

圖5 模型修正前后主要測量參數(shù)誤差對比Fig.5 Difference comparison of engine testing points

模型修正前后發(fā)動機推力值誤差的變化趨勢如圖6所示。

表5 穩(wěn)態(tài)模型計算結(jié)果與試車臺數(shù)據(jù)對比

圖6 各工作狀態(tài)下推力值在模型修正前、后誤差對比Fig.6 Thrust difference comparison of multiple testing points

從模型計算結(jié)果與試車臺數(shù)據(jù)對比可以看出：在完成基于耦合法的設計點參數(shù)匹配后，發(fā)動機部件特性圖的整體精度有所提升，但是僅在發(fā)動機設計點附近參數(shù)匹配精度誤差較小，且隨著工作點逐漸遠離設計點，模型匹配誤差逐漸增大，最大誤差達到7.6%。在進行發(fā)動機部件特性圖的多點匹配修正后，模型匹配推力最大誤差由7.6%降低到1.86%，發(fā)動機各站位參數(shù)均滿足穩(wěn)態(tài)模型誤差小于3%的模型精度要求。

5 結(jié)論

基于試車臺數(shù)據(jù)，在耦合法的基礎上進行了發(fā)動機部件特性圖的多點匹配修正，得出如下結(jié)論。

(1)采用變步長牛頓迭代法可以有效加快模型迭代收斂進程，提高模型迭代收斂速度，以獲取模型最優(yōu)解。

(2)基于發(fā)動機試車臺數(shù)據(jù)，在耦合法的基礎上進行了發(fā)動機部件特性的多點匹配修正，穩(wěn)態(tài)模型平均誤差由5.163 7%降低到多點匹配后的0.529 7%，滿足穩(wěn)態(tài)模型誤差小于3%的模型精度，滿足工程應用要求。

(3)利用試車臺數(shù)據(jù)可以繪制出發(fā)動機高功率穩(wěn)態(tài)下的共同工作線，可以利用修正后的部件特性圖及發(fā)動機穩(wěn)態(tài)模型計算結(jié)果實現(xiàn)發(fā)動機狀態(tài)監(jiān)控與性能分析，為發(fā)動機下發(fā)計劃的制定提供參考。