趙洪利， 高晶東

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

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

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

1 發(fā)動(dòng)機(jī)試車臺(tái)數(shù)據(jù)

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

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

表1 某型渦扇發(fā)動(dòng)機(jī)站位說明

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

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

表2 部分該型渦扇發(fā)動(dòng)機(jī)試車臺(tái)數(shù)據(jù)

2 牛頓迭代法設(shè)計(jì)

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

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

φ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ā)動(dòng)機(jī)依次表示風(fēng)扇、低壓壓氣機(jī)、高壓壓氣機(jī)、高壓渦輪、低壓渦輪的增壓比或落壓比；n1和n2分別為低壓轉(zhuǎn)子相對(duì)轉(zhuǎn)速和高壓轉(zhuǎn)子相對(duì)轉(zhuǎn)速，等于實(shí)測(cè)轉(zhuǎn)速值除以最大轉(zhuǎn)速值。

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

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

(2)

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

(3)

選取歐幾里得范數(shù)作為解的誤差檢驗(yàn)?zāi)繕?biāo)函數(shù)，即：

(4)

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

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

3 部件特性圖修正

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

3.1 部件特性圖設(shè)計(jì)點(diǎn)耦合

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

(5)

(6)

(7)

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

以風(fēng)扇為例，將部件特性圖按照設(shè)計(jì)點(diǎn)進(jìn)行耦合因子修正，修正前后部件特性圖對(duì)比如圖3所示，求得耦合因子結(jié)果如表3所示。

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

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

3.2 部件特性圖非設(shè)計(jì)點(diǎn)耦合

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

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

表4 非設(shè)計(jì)點(diǎn)待優(yōu)化的耦合因子修正項(xiàng)

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

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

圖4 非設(shè)計(jì)點(diǎn)風(fēng)扇部件特性圖修正前、后對(duì)比Fig.4 Comparison of off-design point correction

最后，達(dá)到模型計(jì)算參數(shù)與目標(biāo)參數(shù)匹配的目的。在進(jìn)行穩(wěn)態(tài)模型迭代計(jì)算時(shí)，以式(8)作為誤差檢驗(yàn)?zāi)繕?biāo)函數(shù)。

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

(8)

式(8)中：j為非設(shè)計(jì)點(diǎn)待優(yōu)化的耦合因子修正項(xiàng)個(gè)數(shù)；x為目標(biāo)參數(shù)，下標(biāo)cal,test分別表示模型計(jì)算結(jié)果與實(shí)際試車臺(tái)數(shù)據(jù)。

4 修正結(jié)果分析

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

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

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

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

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

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

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

5 結(jié)論

基于試車臺(tái)數(shù)據(jù)，在耦合法的基礎(chǔ)上進(jìn)行了發(fā)動(dòng)機(jī)部件特性圖的多點(diǎn)匹配修正，得出如下結(jié)論。

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

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

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