賈忠湖，郭衛(wèi)剛，馬 力

航空發(fā)動(dòng)機(jī)部件特性修正技術(shù)研究

賈忠湖1，郭衛(wèi)剛1，馬 力2

（1.海軍航空工程學(xué)院飛行器工程系，山東煙臺(tái)264001；2.海軍駐沈陽(yáng)地區(qū)發(fā)動(dòng)機(jī)專業(yè)軍事代表室，沈陽(yáng)110043）

針對(duì)傳統(tǒng)部件特性修正方法未考慮發(fā)動(dòng)機(jī)多狀態(tài)導(dǎo)致修正精度不高的問(wèn)題，提出了1種基于粒子群優(yōu)化和滑動(dòng)最小二乘法的某型發(fā)動(dòng)機(jī)部件特性修正方法。該方法利用粒子群優(yōu)化算法分別求得在發(fā)動(dòng)機(jī)不同狀態(tài)下的修正系數(shù)，并以這些系數(shù)為基礎(chǔ)，采用滑動(dòng)最小二乘方法擬合修正系數(shù)曲面，從系數(shù)曲面上獲取原有部件特性圖上各點(diǎn)對(duì)應(yīng)的修正系數(shù)，從而得到修正特性。試驗(yàn)結(jié)果表明：該方法克服了傳統(tǒng)方法的不足，提高了特性修正精度，為開(kāi)展單機(jī)監(jiān)控和視情維修提供準(zhǔn)確的部件數(shù)據(jù)基礎(chǔ)。

部件特性修正;粒子群算法;滑動(dòng)最小二乘法;航空發(fā)動(dòng)機(jī)

0 引言

不同型號(hào)的發(fā)動(dòng)機(jī)有著不同的特性，甚至同一類型、同一批次的發(fā)動(dòng)機(jī)由于制造和安裝工藝的差異，其部件特性也會(huì)有所差異，使得每臺(tái)發(fā)動(dòng)機(jī)的性能不完全相同。因此，發(fā)動(dòng)機(jī)模型的仿真結(jié)果與其真實(shí)性能存在一定誤差。因此，有必要針對(duì)單臺(tái)發(fā)動(dòng)機(jī)的性能對(duì)模型進(jìn)行修正，使模型能夠與真實(shí)發(fā)動(dòng)機(jī)匹配[1-5]。

國(guó)內(nèi)外專家提出了一些針對(duì)發(fā)動(dòng)機(jī)部件特性的修正方法。文獻(xiàn)[6]提出了經(jīng)典的修正因子法，可根據(jù)試驗(yàn)數(shù)據(jù)，采用優(yōu)化算法求出；文獻(xiàn)[7]發(fā)展了1種預(yù)測(cè)發(fā)動(dòng)機(jī)部件特性的自適應(yīng)模型方法；文獻(xiàn)[8]分析了部件特性自適應(yīng)模型的不足，利用“部件特性刪除法”求解修正因子。

經(jīng)典的修正因子法最大的改進(jìn)是不再視修正因子為定值，而是隨狀態(tài)變化而變化?；谶@樣的思路，本文提出1種更加有效地綜合利用航空發(fā)動(dòng)機(jī)多個(gè)狀態(tài)試車數(shù)據(jù)的修正方法。

1 部件級(jí)穩(wěn)態(tài)發(fā)動(dòng)機(jī)模型的建立

發(fā)動(dòng)機(jī)部件級(jí)解析模型建立的基本步驟和思路見(jiàn)文獻(xiàn)[9-10]，采用變比熱技術(shù)建立該型發(fā)動(dòng)機(jī)各部件的氣動(dòng)熱力模型。整臺(tái)發(fā)動(dòng)機(jī)在非設(shè)計(jì)點(diǎn)工作時(shí)各部件狀態(tài)是互相制約的，即需滿足共同工作條件。共同工作條件包括：流量平衡、功率平衡、靜壓平衡和轉(zhuǎn)速平衡。穩(wěn)態(tài)工作點(diǎn)求解的基本思路是對(duì)從進(jìn)氣道入口至尾噴管出口截面所建的部件模型進(jìn)行計(jì)算，遇到未知量就假定為已知，稱為獨(dú)立變量，按氣流流過(guò)的順序，完成所有部件氣動(dòng)熱力計(jì)算，從而得到表示發(fā)動(dòng)機(jī)工作的6個(gè)方程，稱為發(fā)動(dòng)機(jī)共同工作方程。根據(jù)該型發(fā)動(dòng)機(jī)的調(diào)節(jié)規(guī)律確定出的6個(gè)待求獨(dú)立變量為：風(fēng)扇工作點(diǎn)位置參數(shù)βcL、高壓轉(zhuǎn)速nH、壓氣機(jī)工作點(diǎn)位置參數(shù)βcH、燃燒室出口總溫、高壓渦輪單位焓降和低壓渦輪單位焓降。則共同工作方程組為

式（1）表示的共同工作方程是包含獨(dú)立變量為未知數(shù)的6個(gè)隱式非線性方程，穩(wěn)態(tài)工作點(diǎn)求解的數(shù)學(xué)問(wèn)題即為求解以獨(dú)立變量為未知數(shù)的非線性隱式方程組。要獲得共同工作方程需進(jìn)行1次完整的發(fā)動(dòng)機(jī)氣動(dòng)熱力計(jì)算，因此，該非線性方程組是非常復(fù)雜的。簡(jiǎn)而言之，求解發(fā)動(dòng)機(jī)部件級(jí)模型的數(shù)學(xué)問(wèn)題就是求解1個(gè)復(fù)雜的多維非線性隱式方程組。

2 基于粒子群優(yōu)化的方程求解

從上節(jié)分析可知，部件特性修正本質(zhì)上是1個(gè)非線性方程組的優(yōu)化過(guò)程。對(duì)這類問(wèn)題的優(yōu)化求解可采用傳統(tǒng)的迭代法，但存在受初值影響較大，不易收斂等不足[11]。并且隨著目標(biāo)方程數(shù)目增加，求解變得更加困難；粒子群算法是近年來(lái)提出的群智能優(yōu)化算法，適合解決復(fù)雜的非線性優(yōu)化問(wèn)題。因此，采用粒子群算法對(duì)共同工作方程組進(jìn)行求解。

粒子群優(yōu)化算法是1種模擬鳥(niǎo)群覓食過(guò)程中遷徙和聚集行為的進(jìn)化算法[12-13]。假設(shè)在1個(gè)n維的目標(biāo)搜索空間中，有N個(gè)粒子組成1個(gè)群體，其中第i個(gè)粒子表示 1 個(gè) n 維的向量 xi=（xi1，xi2，…，xin），i=1，2，…，N，則每個(gè)粒子的位置就是1個(gè)潛在的解。第i個(gè)粒子的飛行速度也是1個(gè)n維向量，記為Vi=（Vi1,Vi2，…，Vin），記第i個(gè)粒子迄今為止搜索到的最優(yōu)位置為 Pi=（Pi1,Pi2，…，Pin），整個(gè)粒子群迄今為止搜索到的最優(yōu)位置為 Pg=（Pg1,Pg2，…，Pgn），w 為慣性權(quán)重。在每一步迭代中，粒子可根據(jù)以下公式更新自己的速度和位置

因子c1和c2分別用來(lái)調(diào)節(jié)該粒子向自身已經(jīng)歷過(guò)的最優(yōu)位置和同伴已尋找到的最優(yōu)位置方向飛行的最大步長(zhǎng)。增大c1和c2可使粒子向全局及個(gè)體最優(yōu)位置飛行的速度加快，從而加快收斂速度，但過(guò)大的c1和c2可能導(dǎo)致由于飛行速度過(guò)快而飛過(guò)目標(biāo)區(qū)域，選擇合適的c1和c2可以加快收斂速度又能避免陷入局部最優(yōu)，通常情況下取c1=c2=2。隨機(jī)數(shù)r1和r2使得粒子飛行時(shí)的速度具有了一定的隨機(jī)性，這種隨機(jī)性提高了粒子在解空間的搜索能力。較小的慣性權(quán)重w可加強(qiáng)粒子局部搜索能力，而較大的w可加快收斂速度，所以可通過(guò)調(diào)節(jié)w來(lái)兼顧收斂速度和收斂精度，w常在[0.4，0.9]之間變化，隨著迭代步數(shù)的增加而減小，使得早期具有較快的收斂速度，而后期又有較強(qiáng)的局部搜索能力。根據(jù)具體問(wèn)題，當(dāng)飛行速度足夠小或達(dá)到預(yù)設(shè)的迭代步數(shù)時(shí)，停止迭代，輸出最優(yōu)解。

式（1）表示的共同工作方程組的求解過(guò)程可轉(zhuǎn)化為如下優(yōu)化問(wèn)題進(jìn)行尋優(yōu)求解

對(duì)于上述優(yōu)化問(wèn)題，利用粒子群算法進(jìn)行求解，粒子 xi=（xi1，xi2，…，xi12），各維度分別對(duì)應(yīng) βcL，nH，βcH，適應(yīng)度函數(shù)選取為

3 基于滑動(dòng)最小二乘法的部件特性修正

利用該型發(fā)動(dòng)機(jī)多個(gè)工作狀態(tài)下的試車數(shù)據(jù)，以風(fēng)扇部件為例進(jìn)行多狀態(tài)部件特性修正。粒子群算法求解部件特性修正系數(shù)得到在已知7個(gè)狀態(tài)下的風(fēng)扇增壓比修正系數(shù)，比較發(fā)現(xiàn)在不同工作范圍得到的修正系數(shù)也不同。因此,傳統(tǒng)方法對(duì)于整個(gè)范圍的部件特性取相同的修正系數(shù)值是不準(zhǔn)確的。通過(guò)計(jì)算可以得到風(fēng)扇部件在較大范圍內(nèi)多個(gè)狀態(tài)的增壓比修正系數(shù)，但受試驗(yàn)所限，不可能對(duì)所有狀態(tài)進(jìn)行試驗(yàn)。因此,對(duì)于已知狀態(tài)之間的任意狀態(tài)，其修正系數(shù)可以根據(jù)該狀態(tài)周圍的已知狀態(tài)特性求取。如果利用所有已知狀態(tài)點(diǎn)的修正系數(shù)，能夠擬合出1個(gè)系數(shù)曲面，則任意狀態(tài)的修正系數(shù)可從該曲面獲得。

近年來(lái)滑動(dòng)最小二乘（MLS，Moving Least Square）在曲面擬合方面取得較好的效果[14-15]，本文即采用滑動(dòng)最小二乘法進(jìn)行曲面擬合。滑動(dòng)最小二乘法是1種基于點(diǎn)的近似方法，具有擬合精度高，通用性強(qiáng)的特點(diǎn)。其擬合函數(shù)采用隨坐標(biāo)變化的系數(shù)向量來(lái)代替?zhèn)鹘y(tǒng)最小二乘法的固定系數(shù)，因此，不存在對(duì)復(fù)雜、離散數(shù)據(jù)量的分段（分塊）擬合和平滑化等問(wèn)題，解決了傳統(tǒng)曲線曲面擬合過(guò)程中存在的困難。

假設(shè)計(jì)算區(qū)域內(nèi)擬合函數(shù)為

式中：x，y，z為計(jì)算區(qū)域的坐標(biāo)變量。

設(shè)在n個(gè)結(jié)點(diǎn)zi=（xi，y）i有已知的實(shí)測(cè)值，分別為σi=σ（xi，y）i，i=1，2，…，n。用矢量形式可表示為

設(shè)擬合函數(shù)可近似表示為

式中：b1（z）為m個(gè)線性無(wú)關(guān)的基函數(shù)，且b1=1；bT（z）為由這些基函數(shù)組成的m維函數(shù)矢量；λ（z）為m維系數(shù)矢量，待求。

與經(jīng)典的擬合方法不同，這里的待定系數(shù)λ（z）是隨著域內(nèi)變量z變化的。首先固定1點(diǎn)z^，稱為估值點(diǎn)。研究在這點(diǎn)鄰域內(nèi)的局部近似

所要求的近似函數(shù)與真實(shí)函數(shù)在已知的n個(gè)實(shí)測(cè)點(diǎn)上的加權(quán)平方范數(shù)最小，即

其中

由于基函數(shù)bi（z）已知，則測(cè)值點(diǎn)場(chǎng)坐標(biāo)zi已知，因此，矩陣B是由確定數(shù)值元素構(gòu)成的。權(quán)矩陣W是對(duì)角矩陣，其元素Wi（z^-zi）是非負(fù)的可導(dǎo)函數(shù)。為減少計(jì)算量，結(jié)點(diǎn)權(quán)函數(shù)只在其影響區(qū)域B（zi）內(nèi)非零，即

可解出

將式（14）代到式（9），得到在估值點(diǎn)（z^）附近的近似函數(shù)

對(duì)每個(gè)估值點(diǎn)z^必須確定1個(gè)相應(yīng)的系數(shù)矢量λ（z^），這是因?yàn)闄?quán)矩陣在數(shù)值上隨估值點(diǎn)的不同而變化。現(xiàn)在使量z^成為計(jì)算區(qū)域內(nèi)的任意1點(diǎn)，即z^可用z代替，則得到全局滑動(dòng)最小二乘法的近似函數(shù)

如果W（z）是常數(shù)矩陣，式（17）是1個(gè)經(jīng)典的、非活動(dòng)的、加權(quán)最小二乘法的近似函數(shù)。由于權(quán)函數(shù)是非線性函數(shù)，該式不能表述為1個(gè)解析函數(shù)，只能用數(shù)值方法求解。

權(quán)函數(shù)采取以下形式

式中：m為影響區(qū)域內(nèi)包含的結(jié)點(diǎn)個(gè)數(shù)；Si為與估值點(diǎn)距離最近的結(jié)點(diǎn)距離；Sj≠i為影響區(qū)域內(nèi)其它結(jié)點(diǎn)與估值點(diǎn)的距離。

4 試驗(yàn)分析

利用MLS和已知的修正系數(shù)點(diǎn)擬合計(jì)算獲得風(fēng)扇的增壓比修正系數(shù)曲面，如圖1所示（圖中曲面只畫出風(fēng)扇特性覆蓋的范圍）。用同樣方法獲得風(fēng)扇的流量特性修正系數(shù)曲面，如圖2所示。

利用得到的相應(yīng)修正系數(shù)曲面，可以獲得精確的修正特性圖。以風(fēng)扇增壓比特性為例，對(duì)于任何1個(gè)轉(zhuǎn)速線上的點(diǎn)（每個(gè)轉(zhuǎn)速線上有10個(gè)點(diǎn)），則依據(jù)其在原來(lái)特性圖上的坐標(biāo)，對(duì)應(yīng)到2個(gè)修正系數(shù)曲面上，找到平面上相應(yīng)位置對(duì)應(yīng)到特性曲面的z軸值即為該點(diǎn)的修正系數(shù)。利用所求的原特性圖上所有點(diǎn)的修正系數(shù)計(jì)算出修正后的風(fēng)扇增壓比和效率。修正前、后的風(fēng)扇增壓比和效率特性分別如圖3、4所示。

圖1 風(fēng)扇的增壓比修正系數(shù)曲面

圖2 風(fēng)扇的流量修正系數(shù)曲面

圖3 修正前、后的風(fēng)扇增壓比特性

為了驗(yàn)證本文方法得到修正特性的精度，分別利用經(jīng)典方法、曲線擬合方法和本文方法得到修正部件特性，進(jìn)行發(fā)動(dòng)機(jī)穩(wěn)態(tài)性能仿真計(jì)算。3種方法的發(fā)動(dòng)機(jī)模型非線性方程組均采用粒子群算法求解，算法均迭代運(yùn)算100次，并取相同的參數(shù)設(shè)置。以某渦扇發(fā)動(dòng)機(jī)2個(gè)不同的狀態(tài)I、II來(lái)驗(yàn)證由不同方法獲得部件特性的精度，2個(gè)狀態(tài)的對(duì)應(yīng)的轉(zhuǎn)速為：狀態(tài)I：Nh=94.35%,N1=86.56%；狀態(tài)II：Nh=92.59%,N1=83.54%。

圖4 修正前、后的風(fēng)扇效率特性

對(duì)于這2個(gè)狀態(tài)，分別基于不同方法得到的部件修正特性，利用粒子群算法進(jìn)行穩(wěn)態(tài)建模仿真。每種情況分別仿真計(jì)算20次，將運(yùn)算結(jié)果對(duì)應(yīng)的適應(yīng)度函數(shù)平均值顯示在表1中?？梢钥闯鲈诓煌瑺顟B(tài)下的計(jì)算結(jié)果有較大不同。在狀態(tài)I下，用3種方法的修正特性計(jì)算所得結(jié)果相差不大，本文算法的適應(yīng)度略小于其他2種方法的，表明本文方法略優(yōu)。然而對(duì)應(yīng)狀態(tài)II，經(jīng)典方法和曲線擬合方法的計(jì)算結(jié)果變差，而采用本文方法計(jì)算取得最好的結(jié)果（適應(yīng)度最小，表明計(jì)算誤差最小）。

表1 不同狀態(tài)下3種方法得到的適應(yīng)度函數(shù)平均值

5 結(jié)論

（1）研究了1種更加有效的綜合利用多個(gè)狀態(tài)試車數(shù)據(jù)的發(fā)動(dòng)機(jī)部件特性修正方法，該方法利用粒子群優(yōu)化算法分別求得在發(fā)動(dòng)機(jī)不同狀態(tài)下的修正系數(shù)，采用滑動(dòng)最小二乘方法擬合修正系數(shù)曲面，從而得到修正特性。

（2）該方法能夠提高修正精度，優(yōu)于其他傳統(tǒng)方法，為開(kāi)展基于單機(jī)監(jiān)控和視情維修提供了準(zhǔn)確的部件數(shù)據(jù)。

[1]Kamboukos P，Mathioudakis K.Comparison of linear and nonlinear gas turbine performance diagnostics[J].Journal of Engineering for Gas Turbines and Power，2005，127（1）：49-56.

[2]Zwebek A，Pilidis P.Degradation effects on combined cycle power plant performance-Part I:gas turbine cycle component degradation effects[J].Journal of Engineering for Gas Turbines and Power，2003，125（7）：651-657.

[3]Li Y G，Nilkitsaranont P.Gas turbine performance prognostic for condition based maintenance[J].Applied Energy，2009，86（10）：2152-2161.

[4]Volponi A J，DePold H，Ganguli R.The use of Kalman filter and neural network methodologies in gas turbine performance diagnostics:a comparative study[J].Journal of Engineering for Gas Turbines and Power，2003，125（7）：917-924.

[5]Fasching W A.CF6 jet engine performance improvement summary report[R].NASA-CR-2165612，1982.

[6]Stamatis A，Mathioudakis K，Papailiou K D.Adaptive simulation of gas turbine performance[J].Journal of Engineering for Gas Turbines and Power，1990，112（2）：168-175.

[7]吳虎，肖洪，蔣建軍.渦扇發(fā)動(dòng)機(jī)部件特性自適應(yīng)模擬[J].推進(jìn)技術(shù)，2005，26（5）：430-433.WU Hu，XIAO Hong，JIANG Jianjun.Adaptive simulation of turbofan engine component performance[J].Journal of Propulsion Technology，2005，26（5）：430-433.（in Chinese）

[8]陳玉春，黃興，徐思遠(yuǎn)，等.渦輪發(fā)動(dòng)機(jī)部件特性自適應(yīng)模型的確定方法[J].推進(jìn)技術(shù)，2008，29（2）：214-218.CHEN Yuchun，HUANG Xing，XU Siyuan，et al.A method of definition for the adaptive simulation model of gas turbine performance[J].Journal of Propulsion Technology，2008，29 （2）：214-218.（in Chinese）.

[9]廉小純，吳虎.航空燃?xì)廨啓C(jī)原理 [M].北京：國(guó)防工業(yè)出版社，2001：45-74.LIAN Xiaochun，WU Hu.Aviation gas turbine principle[M].Beijing：Defense Industry Press，2001：45-74.（in Chinese）

[10]Kong C，Ki J.A new scaling method for component maps of gas turbine using system identification[J].Journal of Engineering for Gas Turbines and Power，2003，125（1）：958-979.

[11]Li Y G，Pilidis P.GA-based design-point performance adaptation and its comparison with ICM-based approach [J].Applied Energy，2010，87（1）：340-348.

[12]Kennedy J，Eberhart R C.Particle swarm optimization[C]//Proceedings of IEEE International Conference on Neural Networks.Perth，Australia，1995：1942-1948.

[13]Eberhart R C，Kennedy J.A new optimizer using particle swarm theory[C]//Proceedings of the Sixth International Symposium on Micro Machine and Human Science，Nagoya:Japan，1995：39-43.

[14]Levin D.The approximation power of moving least-squares[J].Mathematics of Computation，1998，67（10）：1517-1531.

[15]Most T，Bucher C.A moving least squares weighting function for the element-free Galerkin method which almost fulfills essential boundary conditions[J].Structural Engineering and Mechanics，2005，21（3）：315-321.

Research on Correction of Aeroengine Component Characteristics

JIA Zhong-hu1,GUO Wei-gang1,MA Li2
（1.Department of Aircraft Engineering,Naval Aeronautical Engineering Institute,Yantai Shandong 264001,China;2.Navy Customer Representative Office of Engine in Shenyang,Shenyang 110043,China）

Aiming at the problem that the conventional correction methods doesn't consider the multi-condition of engine,the correction of component characteristic based on particle swarm and Moving Least Squares(MLS)was presented.The particle swarm method was used to solve the correction coefficients under different conditions.The MLSmethod was used to fit the curve consisted of correction coefficients,and got the correction coefficients from the orginal characteristic.The experiment results show the proposed method outperforms the conventional methods and improves the accuracy of component characteristic correction,which provide the accurate data for single aircraft monitoring and maintenance.

component characteristic correction;particle swarm algorithms;Moving Least Square method;aeroengine

V231.91

A

10.13477/j.cnki.aeroengine.2015.01.019

2014-01-21

賈忠湖（1964），男，博士，副教授，研究方向?yàn)轱w行器動(dòng)力學(xué)；E-mail：jiazh586@sina.com。

賈忠湖，郭衛(wèi)剛，馬力.某型發(fā)動(dòng)機(jī)部件特性修正技術(shù)研究[J].航空發(fā)動(dòng)機(jī)，2015，41（1）:94-98.JIA Zhonghu，GUOWeigang，MA Li.Research on correction of aeroenginecomponentcharacteristics[J].Aeroengine，2015，41（1）：94-98.

（編輯：張寶玲）