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

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

0 引言

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

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

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

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

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

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

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

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

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

燃燒室出口和燃?xì)鉁u輪進(jìn)口燃?xì)馕锢砹髁科胶夥匠?/p>

燃?xì)鉁u輪出口和動力渦輪進(jìn)口燃?xì)馕锢砹髁科胶夥匠?/p>

動力渦輪出口和尾噴管進(jìn)口燃?xì)馕锢砹髁科胶夥匠?/p>

軸功率平衡方程

式中：Wg為燃?xì)馕锢砹髁?，下?biāo)數(shù)字為相應(yīng)截面編號；Ngt為燃?xì)鉁u輪輸出功率；Ncp為組合壓氣機(jī)消耗功率；Nex為附件提取功率；ηmgt為燃?xì)鉁u輪軸機(jī)械效率；εi（i=1,2,3,4）為平衡方程的殘差。

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

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

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

2.1 待修正參數(shù)

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

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

2.2 目標(biāo)函數(shù)

選擇燃油物理流量Wf、組合壓氣機(jī)增壓比πcp、發(fā)動機(jī)進(jìn)口空氣物理流量Wa2、動力渦輪進(jìn)口燃?xì)饪倻豑45以及輸出功率Ne作為目標(biāo)參數(shù)。

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

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

目標(biāo)函數(shù)定義為

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

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

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

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

速度表示為

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

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

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

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

所有粒子的速度和位置具有上下限，由實(shí)際問題決定。由式（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ā)動機(jī)部件特性。

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

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

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

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

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

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

3 部件特性圖修正方法

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

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

圖4 壓氣機(jī)流量-壓比特性

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖11 修正前、后燃?xì)鉁u輪流量特性對比

圖12 修正前、后燃?xì)鉁u輪效率特性對比

4 結(jié)論

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

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

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

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