渦軸發(fā)動機部件特性修正及更新方法

董 楨，周文祥，潘慕絢 ，黃開明

（1.南京航空航天大學江蘇省航空動力系統(tǒng)重點實驗室，南京210016；2.中國航發(fā)湖南動力機械研究所，湖南株洲412002）

0 引言

非線性部件級模型是航空發(fā)動機及其控制系統(tǒng)研制過程中不可或缺的仿真工具。模型的精度是建模過程中的1項關(guān)鍵技術(shù)指標，直接決定該模型的應(yīng)用價值和適用范圍。借助高精度的航空發(fā)動機性能仿真模型，通過計算機仿真及優(yōu)化，可以大大縮短發(fā)動機數(shù)控系統(tǒng)的設(shè)計及測試周期，降低研發(fā)成本，同時還可避免不必要的試車風險[1-2]。

影響發(fā)動機部件級模型精度的1個重要因素是發(fā)動機部件特性。對發(fā)動機而言，零部件制造和裝配過程必然存在著公差[3]，并且在發(fā)動機壽命期內(nèi)自然磨損、疲勞、腐蝕或積垢等原因也會引起相關(guān)部件的性能退化[4]，從而導致設(shè)計出的發(fā)動機期望部件特性與實際裝機的部件特性之間產(chǎn)生差異。若以期望部件特性為依據(jù)建立的數(shù)學模型表現(xiàn)出的整機性能，與實際發(fā)動機整機性能之間也必然存在一定差異。因此，模型修正技術(shù)對于提高發(fā)動機部件級模型精度非常重要。

國外學者早在20世紀90年代初就開始進行修正發(fā)動機部件特性和模型參數(shù)的方法研究。Stamatis A等[5-6]提出1種基于發(fā)動機實際測量數(shù)據(jù)的方法，將發(fā)動機模型的仿真輸出和試車數(shù)據(jù)誤差的加權(quán)平方和作為目標函數(shù)，采用非線性廣義最小殘差法求解共同工作方程；同樣在試車數(shù)據(jù)的基礎(chǔ)上，Joachim Kurzk[7-8]著眼于大氣條件對發(fā)動機各截面總溫、總壓及流量等氣路參數(shù)的影響，給出詳細的修正方法。之后，優(yōu)化算法逐漸引入部件特性修正中。Y.G.Li[9]采用牛頓-拉夫森算法求解不可測參數(shù)的修正因子，得到準確的不可測參數(shù)的真實值；Changduk Kong[10-12]從部件特性圖著手，發(fā)現(xiàn)壓氣機流量可以近似表示為壓比的3次多項式，效率可以近似表示為流量的3次多項式，其中多項式系數(shù)通過遺傳算法優(yōu)化得到。國內(nèi)學者在發(fā)動機部件級模型修正方面也進行了深入研究。段守付等[13]在模型中引入加權(quán)函數(shù)和一系列修正因子，改變修正因子的取值使加權(quán)函數(shù)最小，實現(xiàn)對發(fā)動機部件特性和模型的最優(yōu)修正；陳策等[14]建立BP神經(jīng)網(wǎng)絡(luò)對部件特性數(shù)據(jù)進行識別學習，實現(xiàn)特性數(shù)據(jù)的精確插值和對未知特性數(shù)據(jù)的推測；吳虎等[15]提出1種預測發(fā)動機部件特性的自適應(yīng)模型方法，以通用特性為基礎(chǔ)，運用單純形優(yōu)化方法預測出不同飛行條件下渦扇發(fā)動機的部件特性。上述修正方法在原理上具有相似性，但在求解非線性方程組或優(yōu)化算法方面存在一定差異。

本文利用具有高效率和全局搜索能力的粒子群優(yōu)化算法（Particle Swarm Optimization，PSO）來修正渦軸發(fā)動機的部件特性。在基本不改變部件級模型核心計算代碼的情況下，使用PSO算法計算出部件特性圖和其他參數(shù)的修正因子。

1 渦軸發(fā)動機建模

擬建模的渦軸發(fā)動機由進氣道（包括粒子分離器）、組合壓氣機（包括軸流壓氣機和離心壓氣機）、燃燒室、燃氣渦輪、動力渦輪和尾噴管組成，如圖1所示，各截面定義見表1。

圖1 渦軸發(fā)動機結(jié)構(gòu)

表1 渦軸發(fā)動機各截面定義

渦軸發(fā)動機部件級模型的輸入?yún)?shù)是大氣壓力、環(huán)境溫度、高度、馬赫數(shù)、燃氣渦輪轉(zhuǎn)速和動力渦輪轉(zhuǎn)速。輸出參數(shù)是燃油流量、各截面氣動參數(shù)（總溫、總壓等）和性能參數(shù)（動力渦輪功率、耗油率等）。在建立所有部件的氣動熱力學模型后，選定4個仿真初猜值（壓氣機、燃氣渦輪和動力渦輪壓比及燃氣渦輪轉(zhuǎn)速），其對應(yīng)的4個共同工作方程為

燃燒室出口和燃氣渦輪進口燃氣物理流量平衡方程

燃氣渦輪出口和動力渦輪進口燃氣物理流量平衡方程

動力渦輪出口和尾噴管進口燃氣物理流量平衡方程

軸功率平衡方程

式中：Wg為燃氣物理流量，下標數(shù)字為相應(yīng)截面編號；Ngt為燃氣渦輪輸出功率；Ncp為組合壓氣機消耗功率；Nex為附件提取功率；ηmgt為燃氣渦輪軸機械效率；εi（i=1,2,3,4）為平衡方程的殘差。

利用牛頓-拉夫森法迭代求解共同工作方程組就可以得到穩(wěn)態(tài)工作點。由上述4個共同工作方程結(jié)合各部件氣動熱力計算方程，構(gòu)成發(fā)動機的穩(wěn)態(tài)模型。在發(fā)動機穩(wěn)態(tài)工作時，若共同工作方程的殘差|εi|＜ε0，i=1,2,3,4（其中 ε0為控制精度，本文取 10-3），則認為方程收斂，停止迭代。

2 粒子群優(yōu)化算法

采用粒子群優(yōu)化算法修正發(fā)動機部件級模型，其原理是利用PSO優(yōu)化發(fā)動機模型的修正因子，使目標函數(shù)最小化或使適應(yīng)度函數(shù)最大化。

2.1 待修正參數(shù)

對渦軸發(fā)動機部件級模型而言，待修正的參數(shù)包括旋轉(zhuǎn)部件的部件特性（包括壓比、換算質(zhì)量流量和效率）、典型流道部件的總壓恢復系數(shù)、壓氣機引、放氣比例以及燃燒室效率。

修正設(shè)計點時，對上述參數(shù)都進行修正；修正非設(shè)計點時，本文只修正旋轉(zhuǎn)部件的部件特性，其他待修正參數(shù)采用設(shè)計點的修正結(jié)果。

2.2 目標函數(shù)

選擇燃油物理流量Wf、組合壓氣機增壓比πcp、發(fā)動機進口空氣物理流量Wa2、動力渦輪進口燃氣總溫T45以及輸出功率Ne作為目標參數(shù)。

發(fā)動機測量參數(shù)的誤差定義為

式中：ycal為計算參數(shù)；yact為基準參數(shù)。

目標函數(shù)定義為

式中：ai為各目標參數(shù)的權(quán)重系數(shù)，本文在優(yōu)化過程中暫取1。

2.3 粒子群算法優(yōu)化原理

粒子群算法是模擬鳥群覓食的過程[16-17]。每個優(yōu)化問題的解都是搜索空間的1只鳥，稱為“粒子”。每個粒子具有由優(yōu)化函數(shù)確定的適應(yīng)度值。每個粒子也具有1個速度，決定搜索的方向和距離。所有粒子追隨當前最優(yōu)粒子在解空間中搜索。

假設(shè)D維空間中的粒子群中的第i個粒子的位置表示為

速度表示為

當前粒子群中的最優(yōu)個體粒子表示為

歷代全局最優(yōu)個體粒子表示為

粒子群中的所有粒子根據(jù)式（11）和（12）更新速度和位置

式中：r1和r2取（0，1）中的隨機數(shù)；c1和c2為學習因子，通常取2；ω為慣性權(quán)重，取較大的值意味著粒子群具有較強的全局搜索能力，取較小的值意味著粒子群局部搜索能力達到最優(yōu)。

所有粒子的速度和位置具有上下限，由實際問題決定。由式（11）和（12）組成的粒子群算法稱為基本粒子群算法。

在此基礎(chǔ)上，如果慣性權(quán)重的值根據(jù)迭代次數(shù)而減少，則算法被稱為慣性權(quán)重線性遞減粒子群算法[18]。一般來說，ω的初始值設(shè)為0.9，然后根據(jù)迭代次數(shù)將其線性減小到0.4。粒子群優(yōu)化算法的慣性權(quán)重線性遞減是粒子群算法中比較常見的算法模型。本文采用慣性權(quán)重線性遞減粒子群優(yōu)化算法優(yōu)化渦軸發(fā)動機部件特性。

慣性權(quán)重定義為

式中：i為當前迭代次數(shù)；n為總迭代次數(shù)。

本文粒子群總優(yōu)化代數(shù)設(shè)置為300，所有待修正因子優(yōu)化范圍統(tǒng)一設(shè)置為[0.98，1.02]。適應(yīng)度在初始100代以內(nèi)急劇增加，然后減慢，在150代左右達到最大，如圖2所示。

部件特性修正流程如圖3所示。

需要特別說明的是，粒子群優(yōu)化迭代過程中更新的是部件特性圖，而不是部件流量或效率特性的修正因子，這樣做的好處是一旦優(yōu)化結(jié)束，可以直接輸出修正后的部件特性圖和其他有關(guān)模型參數(shù)，如總壓損失、引放氣比例等。

圖2 粒子群適應(yīng)度變化

3 部件特性圖修正方法

壓氣機流量-壓比特性如圖4所示。圖中描繪出待修正的一些穩(wěn)態(tài)工作點，A為設(shè)計點，B、C、D為非設(shè)計點。對于設(shè)計點A，部件特性圖將根據(jù)模型計算的結(jié)果與發(fā)動機設(shè)計點參數(shù)進行整體縮放。非設(shè)計點則是根據(jù)相對換算轉(zhuǎn)速從高到低的順序依次進行修正。例如，當修正工作點B時，保持L1不變并縮放L2，L4，…，Ln，以免影響設(shè)計點A的仿真精度。然后固定L1和L2，繼續(xù)修正L4及更低的等換算轉(zhuǎn)速線來修正工作點C。

上述修正方法假定在壓氣機2條等換算轉(zhuǎn)速線之間只存在1個待修正工作點，這一假設(shè)條件在模型修正過程中有時是不滿足的。圖4顯示在壓氣機L2和L42條等換算轉(zhuǎn)速線之間存在C和D 2個工作點，針對這種情況，可以先選擇距離設(shè)計點A較遠的工作點D作為待修正點，對L4及其以下的等換算轉(zhuǎn)速線進行整體修正，然后用工作點C來檢驗修正好的部件特性是否達到仿真精度要求。如若不滿足要求，可以在點C和點D之間補充1條等換算轉(zhuǎn)速線L3。這時問題轉(zhuǎn)化成2條等換算轉(zhuǎn)速線之間只有1個待修正點的情況，可根據(jù)前述方法對補作的特性線L3及其以下特性線進行修正，修正完畢后，繼續(xù)根據(jù)點D的仿真結(jié)果修正L4及其以下轉(zhuǎn)速特性線，如圖5所示。

圖4 壓氣機流量-壓比特性

圖5 補作等換算轉(zhuǎn)速線后壓氣機流量-壓比特性

為了驗證提出的部件特性修正方案是否可行，本文從發(fā)動機性能分析軟件GasTurb提取渦軸發(fā)動機通用部件特性圖，以GasTurb軟件計算結(jié)果作為模型修正基準值。在地面標況下（ISA H=0 km，Ma=0）選取燃氣發(fā)生器轉(zhuǎn)子相對物理轉(zhuǎn)速分別為1、0.975、0.941、0.907的 4個工作點 A、B、C、D 來修正通用部件特性圖。設(shè)計點A修正前、后各目標參數(shù)的穩(wěn)態(tài)誤差見表2。

表2 設(shè)計點A修正前、后穩(wěn)態(tài)誤差對比（ISA H=0 km，Ma=0，n g=100%，n p=100%）

由表中可知，修正設(shè)計點A后，模型精度大大提高，設(shè)計點處模型平均穩(wěn)態(tài)仿真誤差由4.7%下降至0.2%。設(shè)計點A修正前、后壓氣機流量特性變化如圖6所示。

圖6 設(shè)計點A修正前、后壓氣機流量特性對比

接著修正穩(wěn)態(tài)工作點B。穩(wěn)態(tài)工作點B修正前、后各目標參數(shù)的穩(wěn)態(tài)誤差見表3。

表3 穩(wěn)態(tài)工作點B修正前、后穩(wěn)態(tài)誤差對比（ISA H=0 km，Ma=0，n g=97.5%，n p=100%）

穩(wěn)態(tài)工作點B修正前、后壓氣機流量特性變化如圖7所示。

由于穩(wěn)態(tài)工作點C、D同時在2條等換算轉(zhuǎn)速線之間，先修正穩(wěn)態(tài)工作點D。穩(wěn)態(tài)工作點D修正前、后穩(wěn)態(tài)誤差見表4。

修正后 A、B、D 3 個穩(wěn)態(tài)工作點目標參數(shù)相對誤差都在1%以內(nèi)，滿足精度要求。然后檢驗修正后的部件特性是否能夠使穩(wěn)態(tài)工作點C也滿足指標要求。修正后穩(wěn)態(tài)工作點C穩(wěn)態(tài)誤差見表5。

圖7 穩(wěn)態(tài)工作點B修正前、后壓氣機流量特性對比

表4 穩(wěn)態(tài)工作點D修正前、后穩(wěn)態(tài)誤差對比（ISA H=0 km，Ma=0，n g=90.7%，n p=100%）

表5 修正后穩(wěn)態(tài)工作點C穩(wěn)態(tài)誤差（ISA H=0 km，Ma=0，n g=94.1%，n p=100%）

從表中可見，Wf和T45的相對誤差仍高于1%。因為這2個參數(shù)是反映渦軸發(fā)動機經(jīng)濟性及整機性能指標非常重要的參數(shù)，很有必要進一步提升其仿真精度，所以此處考慮在點C、D之間補作1條等換算轉(zhuǎn)速線。補作的等換算轉(zhuǎn)速線需要根據(jù)實際情況來確定，這里取相鄰等換算轉(zhuǎn)速的中值。

補作等換算轉(zhuǎn)速線后穩(wěn)態(tài)工作點C、D修正前、后穩(wěn)態(tài)誤差見表6、7。修正前、后壓氣機流量特性對比如圖8、9所示，其中，L3為補作的等換算轉(zhuǎn)速線。

表6 補作等換算轉(zhuǎn)速線后工作點C修正前、后誤差對比（ISA H=0 km，Ma=0，n g=94.1%，n p=100%）

表7 補作等換算轉(zhuǎn)速線后工作點D修正前、后誤差對比（ISA H=0 km，Ma=0，n g=90.7%，n p=100%）

以上仿真結(jié)果表明，修正后點A、B、C、D目標參數(shù)穩(wěn)態(tài)仿真誤差均在1%以內(nèi)，滿足精度要求并且結(jié)果與預期相同。

圖8 穩(wěn)態(tài)工作點C修正前、后壓氣機流量特性對比

圖9 穩(wěn)態(tài)工作點D修正前、后壓氣機流量特性對比

壓氣機效率特性的修正與流量特性的修正方法一致，修正前、后壓氣機效率特性對比如圖10所示。

修正渦輪部件特性時，2條等換算轉(zhuǎn)速線之間存在B、C、D 3個穩(wěn)態(tài)工作點。采用與壓氣機部件特性相同的修正方法，由于等換算轉(zhuǎn)速線分布比較密集，修正穩(wěn)態(tài)工作點D之后的特性能夠使這3個穩(wěn)態(tài)工作點同時滿足精度要求。修正前、后渦輪流量特性對比如圖11所示，渦輪效率特性對比如圖12所示。

圖10 修正前、后壓氣機效率特性對比

圖11 修正前、后燃氣渦輪流量特性對比

圖12 修正前、后燃氣渦輪效率特性對比

4 結(jié)論

本文提出基于粒子群優(yōu)化算法（PSO）修正渦軸發(fā)動機部件特性的方法，通過模型計算結(jié)果與Gas-Turb仿真結(jié)果的對比，得到以下結(jié)論：

（1）用粒子群算法優(yōu)化并更新已有的部件特性圖可以有效減少部件級模型的穩(wěn)態(tài)仿真誤差。

（2）在模型修正過程中，不更新特性圖只更新修正因子的方法對修正因子缺乏約束，求解出的數(shù)學解可能存在對應(yīng)的發(fā)動機旋轉(zhuǎn)部件效率大于1的情況，這與實際情況不相符。而本文提出的方法在每次優(yōu)化修正結(jié)束后直接輸出更新的部件特性圖，可以直觀地查看更新后的部件特性，有效避免了這一問題。

本文所提出的部件特性修正方法還有不足之處，例如在修正發(fā)動機非設(shè)計點的過程中，只是對發(fā)動機部件特性中的等換算轉(zhuǎn)速線進行平移或者縮放處理，并沒有改變其形狀或趨勢，有待引入更多的高空穩(wěn)態(tài)工作點，進一步開展改變等換算轉(zhuǎn)速線局部形狀或趨勢的修正方法研究。