宋漢強，錢仁軍，滕懷亮，閆思齊

（1.海軍研究院，上海200436；2.海軍航空大學，山東煙臺264001）

航空發(fā)動機數學模型是發(fā)動機研制過程中的重要組成部分，是發(fā)動機控制系統(tǒng)和性能評估系統(tǒng)的基礎。好的仿真模型能很大程度縮短發(fā)動機的研制周期，節(jié)約研制經費，而模型的精度直接決定了模型準確性及其實用價值[1]。與實際工作過程相比，發(fā)動機建模過程中進行了諸多簡化，加上制造工藝水平的影響，發(fā)動機模型所輸出的氣路性能參數與實際發(fā)動機試車數據往往存在較大差異。因此，對發(fā)動機數學模型進行修正研究使其能與真實數據匹配是必要的[2]。

國外在20世紀90年代就開始了航空發(fā)動機部件特性和模型參數修正方法的研究，其中，A Stamatis[3-4]最早提出將發(fā)動機實測數據和發(fā)動機模型輸出值的最小誤差作為目標函數，通過迭代部件特性修正系數來實現發(fā)動機模型修正的方法。C Kong[5]在研究中運用壓比的三次多項式表示壓氣機流量，用流量的三次多項式表示效率，最后運用遺傳算法和發(fā)動機實測數據求解出多項式系數，實現對壓氣機特性圖的修正。

國內對發(fā)動機部件級模型修正方法也開展了相關研究，段守付[6]最早提出加權函數自適應建模方法，通過迭代修正因子使得加權函數值最小，實現對發(fā)動機模型參數的修正。白磊[7]基于發(fā)動機試驗數據，運用變分加權最小二乘法對其進行模型辨識分析，從而實現對部件特性的修正。李冬[8]和潘鵬飛[9]都是基于遺傳算法對發(fā)動機部件模型參數進行的修正，朱正琛[10-11]采用微分進化算法對發(fā)動機部件特性進行修正，其均使修正后的模型輸出與實驗數據更一致，提高了模型精度。

本文基于發(fā)動機仿真模型，選擇合適部件特性作為修正因子，將發(fā)動機整機性能實測數據和發(fā)動機模型輸出值的最小誤差作為目標函數，并采用具有高效率和強全局搜索能力的粒子群算法（Particle Swarm Optimization，PSO）來優(yōu)化求解目標函數，從而實現對發(fā)動機模型的修正。

1 建立發(fā)動機模型

根據某型渦扇發(fā)動機的氣路結構和各截面之間的氣動熱力學關系，劃分計算截面如圖1所示，各截面定義如表1所示。

圖1 渦扇發(fā)動機截面圖Fig.1 Section diagram of turbofan engine

表1 發(fā)動機各截面定義Tab.1 Definition of engine sections

依據設計點參數建立的發(fā)動機熱力學模型，其輸入參數為大氣溫度、壓力、高度和馬赫數，模型輸出參數為各截面總溫和總壓，推力、燃油流量等性能參數[12]。本文求解發(fā)動機模型所選定的初猜值為低壓轉速、低壓壓比函數、高壓壓比函數、渦輪前溫度、高壓渦輪流量和低壓渦輪流量，求解過程對應的平衡方程分別為：

高壓渦輪流量平衡方程為

低壓渦輪流量平衡方程為

高壓渦輪與壓氣機的功率平衡方程為

低壓渦輪與風扇的功率平衡方程為

混合室入口內涵與外涵的靜壓平衡方程為

尾噴口的面積平衡方程為

式（1）～（6）中：Wg為燃氣流量計算值；W′g為特性數據插值結果；LC為風扇和壓氣機計算功率；LT為渦輪特性數據插值后計算功率；P55和P25分別為混合室入口內涵與外涵的靜壓。

當發(fā)動機處于穩(wěn)態(tài)工作時，模型求解問題轉化為6 個共同工作方程組成的非線性方程組求解問題[13]。本文利用牛頓-拉夫森法進行迭代求解，如殘差滿足要求，則認為結果收斂，模型求解完成[14]。

2 基于PSO的發(fā)動機模型修正方法

2.1 修正因子選擇

部件特性數據是發(fā)動機部件級模型的基礎，特性數據是否符合發(fā)動機真實情況，直接關系發(fā)動機模型的求解精度。因此，發(fā)動機模型修正通常都是通過修正部件特性數據實現精度的提高。

以壓氣機特性中的效率修正因子為例，其表達式為：

式（7）中：ηCH為發(fā)動機原特性數據中效率；ηact,CH為發(fā)動機真實效率；xη為修正因子。

本文以發(fā)動機中間狀態(tài)實測數據為基礎，結合發(fā)動機模型共同工作方程初猜值種類，選定6 個修正因子，分別為風扇流量、風扇效率、壓氣機流量、壓氣機效率、低壓渦輪和高壓渦輪效率，記為X=[ xi] ，i=1,2,…,6。

2.2 目標函數構建

發(fā)動機穩(wěn)態(tài)模型修正方法是通過選擇合適的修正因子X ，修正模型計算過程中的原有特性數據，從而使得模型的輸出參數與實測性能參數誤差減小。因此，通過選擇發(fā)動機實測性能參數構建模型修正的目標函數：

式（8）中：Ycal為模型參數計算值；Yact為實測參數測量值；m 為選擇實測性能參數的個數。

本文選擇推力F 、燃油流量WFT、渦輪后溫度T55、渦輪后壓力P55、壓氣機出口溫度T3、壓氣機出口壓力P3、風扇出口壓力P22作為目標參數。

通過目標函數可知，發(fā)動機模型修正方法就是通過合適的優(yōu)化算法，求解min FC( X )。

2.3 PSO優(yōu)化算法

PSO 算法的基本思想就是由M 個粒子組成的群體，以一定速度飛行，尋找最優(yōu)位置[15]。其中，M 個粒子代表修正因子的M 個候選解，目標函數計算值就是粒子的適應度，目標函數的解空間就是粒子運動空間。搜索過程中，根據單個粒子歷史最優(yōu)點Pbest和群體粒子歷史最優(yōu)點Gbest更新飛行速度，從而實現位置（修正因子的解）更新[16]。

本文將式中學習因子c1和c2均取值為2，其慣性權重w 的初始取值設為0.9，隨著迭代次數增加線性減小到0.4，方法如下：

式（12）中：k 為當前代數；Itermax代表最大進化代數。

結合發(fā)動機模型計算求解過程和實測試車數據，運用PSO算法修正發(fā)動機模型的流程圖如圖2所示。

圖2 運用PSO算法修正發(fā)動機模型流程圖Fig.2 Flow chart of modification of engine model using PSO algorithms

3 仿真結果與分析

3.1 修正前模型計算結果

本文以某發(fā)動機為例，根據其設計點參數建立該發(fā)動機通用模型，在標準大氣條件下計算得到其具體性能參數，并與發(fā)動機整機實測性能進行對比，結果如表2 所示。從表2 中可以看出，模型計算值與發(fā)動機實測數值存在明顯偏差。其中，推力的計算誤差為最大誤差，達到4.85%。

表2 修正前模型計算值與實測數值對比Tab.2 Comparisons between calculatedand measured values of the model before correction

3.2 修正后模型計算結果

采用前文所述PSO 算法對發(fā)動機修正模型的目標函數進行優(yōu)化求解，設定粒子種群數量為30，最大迭代次數為500，計算求出優(yōu)化后的修正因子為X=（0.933 5，1.03，0.999 6，0.999 6，0.947 2，0.911 1），粒子群適應度變化過程見圖3，最小值為1.223 1×10-4。

圖3 PSO算法迭代收斂圖Fig.3 Iterative convergence graph of PSO algorithms

將PSO 優(yōu)化求解得出的修正因子帶入發(fā)動機模型進行計算，得到修正后發(fā)動機模型中間狀態(tài)性能參數的計算值，將其與發(fā)動機實測數值進行對比，結果如表3 所示。從表中可以看出，模型計算值與發(fā)動機實測數值偏差較小。其中，最大誤差只有0.97%。

表3 模型修正后計算值與實測數值對比Tab.3 Comparisons between calculated and measured values after model modification

將發(fā)動機模型修正前與修正后中間狀態(tài)性能計算值與實測值的誤差進行對比，結果如圖4 所示。從圖中可以看出，相比修正前的發(fā)動機模型，修正后的模型性能參數計算值精度明顯提高。其中，后端渦輪等部件的精度明顯比風扇和壓氣機等前端部件提高的多。

圖4 模型修正前與修正后的計算值與實測值誤差對比Fig.4 Error comparison between calculated and measured values before and after model modification

4 結論

本文通過建立某型渦扇發(fā)動機氣動熱力學仿真模型，提出基于粒子群算法（PSO）的發(fā)動機模型修正方法，運用修正因子實現了模型精度的提高，通過發(fā)動機模型計算結果與發(fā)動機整機性能實測數據進行對比，得出以下結論：

1）運用PSO算法對發(fā)動機模型進行修正，能夠顯著提高模型的精度。修正前，發(fā)動機模型的性能計算值與實測數據對比，最大誤差到達4.85%。修正后，模型精度極大提高，最大誤差只有0.97%，修正效果良好。

2）從模型修正效果看，后端部件比前端部件精度提高更為明顯。修正前后結果表明，風扇和壓氣機壓比誤差變化不大，渦輪后溫度和壓力誤差減小很多，推力計算精度極大提高。