楊錕,張鑫,屠秋野,陳飛陽,鄭康
(1.西北工業(yè)大學 動力與能源學院,西安 710129) (2.中國航發(fā)湖南動力機械研究所 發(fā)動機總體研究部,株洲 412002)
航空發(fā)動機健康監(jiān)測的目標是掌握發(fā)動機各部件的性能退化程度,即獲得相應部件在特定工作狀態(tài)下健康參數(shù)的量值,它們通常與旋轉部件的流通能力、做功效率,燃燒室反應效率、燃油泵延遲,高/低壓軸承的摩擦損失,氣路的漏氣比例以及變幾何機構的控制參數(shù)等的變化相關。上述健康參數(shù)通常無法通過傳感器直接測量,因此需要采用氣路分析(Gas Path Analysis)技術并利用發(fā)動機氣路傳感器測得的各旋轉部件轉速以及部分截面的壓力和溫度等對相應的健康參數(shù)進行估計。
自從L.A.Urban[1]將基于線性模型的氣路分析方法引入航空發(fā)動機狀態(tài)分析以來,絕大多數(shù)的氣路分析研究均基于發(fā)動機的穩(wěn)態(tài)過程[2-7]。然而從工程實踐的角度來說,穩(wěn)態(tài)數(shù)據(jù)經(jīng)常難以獲得[8],且在實際某個航段獲得多個穩(wěn)態(tài)工作點的數(shù)據(jù)更不現(xiàn)實。在民航客機特別是國內(nèi)航線上運行的發(fā)動機,存在著相當比例的過渡工作狀態(tài),而這一狀態(tài)在軍用發(fā)動機的工作過程中占比甚至達到了70%[9]。此外,由于發(fā)動機的很多故障特征僅在過渡工作狀態(tài)下才能表現(xiàn)出來,使用相對更容易獲得的過渡態(tài)數(shù)據(jù)能夠更全面的對發(fā)動機的健康狀態(tài)進行評估。
目前過渡態(tài)氣路分析中比較有代表性的工作有:T.Groenstedt[10]發(fā)展了基于最小二乘的擴展常微分代數(shù)方程算法;Li Yiguang[11]提出了以積累誤差構建目標函數(shù)的理論;J.C.Simmons等[12]和J.R.Mc Cusker等[13]分別對過渡態(tài)輸出信號進行了小波分析,前者獲得了故障特征并以此為基礎開發(fā)出一套適合機載系統(tǒng)使用的快速診斷方法,后者實現(xiàn)了對測量參數(shù)的選擇并對性能參數(shù)進行了非線性最小二乘估計。但是上述研究均未對影響發(fā)動機動態(tài)性能的關鍵參數(shù)進行分析,也無法解決在測量傳感器數(shù)目少于待求健康參數(shù)數(shù)目情況下的大偏差發(fā)動機氣路分析問題。
本文首次提出基于過渡工作過程的航空發(fā)動機大偏差欠定氣路分析方法,是為在測量傳感器數(shù)量有限的條件下,實現(xiàn)對包含表征發(fā)動機機動性能在內(nèi)的大量關鍵氣路部件健康參數(shù)精確估計的系統(tǒng)方法。該方法基于發(fā)動機非線性模型,采用序列工作點分析方法,并應用間接遞歸牛頓-拉夫遜法強化非支配分類差分進化算法(Indirectly Recursive Newton-Raphson Method Enhanced Non-dominated Sorting Differential Evolution Algorithm)執(zhí)行健康參數(shù)估計的數(shù)值計算。其優(yōu)點是能夠區(qū)分在同一過渡工作過程中不同功率條件下退化程度的差異,且在提升發(fā)動機大偏差性能退化條件下全局尋優(yōu)能力的同時,平衡演化算法的高計算時間成本。通過對多組序列工作點算例的應用以及對算法性能的分析,驗證本文方法的有效性。
采用基于非線性模型的方法對發(fā)動機的健康狀態(tài)進行估計,研究對象為某型大涵道比分開排氣渦輪風扇發(fā)動機,其布局如圖1所示。
圖1 分排渦扇發(fā)動機示意圖Fig.1 Schematic diagram of separated flow turbofan engine
發(fā)動機模型的數(shù)學描述為一組非線性方程組:
(1)
z(t)=g[(x(t),u(t)]
(2)
式中:x為發(fā)動機性能參數(shù)向量;u為發(fā)動機輸入向量(包括飛行高度、馬赫數(shù)和燃油流量等);z為發(fā)動機測量參數(shù)向量。式(1)左邊部分分量為零。
將上述非線性模型進行擴展,嵌入部件健康參數(shù)向量θ,構建發(fā)動機非線性狀態(tài)分析模型:
(3)
z(t)=g[(x(t),u(t),θ(t)]
(4)
式中:θ各分量由對應部件性能參數(shù)的相對值定義,其形式為
(5)
式中:X為特定的部件性能參數(shù);下標d、c分別表示實際退化發(fā)動機與標準發(fā)動機;p為健康參數(shù)向量維度。
狀態(tài)分析模型中健康參數(shù)偏差的定義為
Δθi,d(t)=θi,d(t)-1 (i=1,…,p)
(6)
在對健康參數(shù)進行估計的過程中,需對非線性模型進行線性逼近并實施迭代計算,其過程如下:
(7)
(8)
(9)
式中:q為測量傳感器數(shù)目。
將式(7)寫為差值形式:
(10)
(11)
(12)
根據(jù)式(10)可以給出健康參數(shù)估計的最小二乘迭代格式:
(13)
式中:k為迭代步數(shù)。
利用式(13)進行迭代計算的原理如圖2所示。
圖2 迭代計算原理示意圖Fig.2 Schematic diagram of iterative computation
考慮式(13)的構造形式,Jacobi矩陣J的行數(shù)應該大于或等于列數(shù),而在測量傳感器數(shù)目較少的情況下,這一條件并不滿足。因此,本文采用基于多工作點分析(Multiple Operating Points Analysis,簡稱MOPA)[14]思想的序列工作點分析技術(Sequential Operating Points Analysis,簡稱SOPA)[15]對Jacobi矩陣進行擴維。SOPA方法在發(fā)動機的某段過渡工作過程中選取特定傳感器布局下一組采樣頻率為10 Hz(氣路傳感器信號的帶寬為5~20 Hz[16])的連續(xù)測量數(shù)據(jù)進行氣路分析。序列工作點選取在發(fā)動機某加速過程工作線上的分布如圖3所示。
圖3 序列工作點選取示意圖Fig.3 Schematic diagram of sequential operating points
在不考慮傳感器故障的前提下,序列工作點的數(shù)量l需滿足如下條件:
p≤l×q
(14)
(15)
(16)
式中:T為t時刻(如圖3所示)對應的采樣點。
基于序列工作點的迭代更新格式為
(17)
本文提出間接遞歸牛頓-拉夫遜法強化非支配分類差分進化算法,其流程如圖4所示,通過結合進化算法的全局搜索能力與確定性算法的局部逼近能力,為大偏差氣路分析問題提供一種解決思路。
在不考慮測量噪聲的情況下,在各序列工作點構造子目標函數(shù):
(18)
多目標進化算法能夠將散布在整個搜索空間的非支配個體逐步向局部最優(yōu)解或全局最優(yōu)解聚集,但是由于演化過程的隨機性,準確判斷非支配個體進入全局最優(yōu)解凸鄰域范圍的時機并不現(xiàn)實,故本文算法采用一種特殊的后驗偏好銜接(A Posteriori Preference Articulation)[17]方式:漸進后驗偏好銜接(Progressively A Posteriori Preference Articulation),即在種群的進化過程中,按照一定的代數(shù)間隔選擇非支配個體執(zhí)行最優(yōu)解決策,如果達到?jīng)Q策目標,則算法結束,否則繼續(xù)演化進程。其中,非支配個體的演化過程采用非支配分類差分進化算法(Non-dominated Sorting Differential Evolution,簡稱NSDE),通過將NSGA-Ⅱ[18]中的遺傳算法(GA)替換為基本差分進化算法(DE/rand/1)構建。
圖4 間接遞歸牛頓-拉夫遜法強化非支配 分類差分進化算法流程Fig.4 Indirectly recursive Newton-Raphson method enhanced nondominated sorting differential evolution flowchart
最優(yōu)解決策的過程采用間接遞歸牛頓—拉夫遜方法(Indirectly Recursive Newton-Raphson Method,簡稱IRNRM),其構建思路如下:
以某型大涵道比雙軸分排渦扇發(fā)動機作為研究對象進行方法驗證,其布局如圖1所示。給定發(fā)動機工作高度和馬赫數(shù),燃油計劃如圖5所示。
圖5 目標發(fā)動機燃油計劃Fig.5 Fuel schedule of target engine
待分析部件性能參數(shù)為風扇、增壓級、壓氣機、高壓渦輪、低壓渦輪的換算流量Γ和效率η,高壓軸承、低壓軸承的摩擦損失μ,壓氣機導葉調(diào)節(jié)角度β,燃燒室燃燒效率η以及燃油泵延遲τ,即p=15,如表1所示(第一列)。將測量傳感器布置給定為(如圖1所示):增壓級出口總壓p25、壓氣機出口總溫T3、壓氣機出口總壓p3,即q=3。在持續(xù)時間為5 s的加速工作過程中,于t分別為1.0、1.5和3.0 s處各選取一組序列工作點(如圖6所示)的測量數(shù)據(jù)分別采用如圖4所示算法進行氣路狀態(tài)分析。其中每組序列工作點包含5個連續(xù)的發(fā)動機狀態(tài)點以滿足式(14)條件,即l=5。
(a) 算例1
(b) 算例2
(c) 算例3圖6 在過渡工作過程中包含t=1.0 s,t=1.5 s 以及t=3.0 s時刻的序列工作點選取示意圖Fig.6 Schematic diagram of sequential operating points including t=1.0 s,t=1.5 s, and t=3.0 s
本文以預先植入性能偏差的發(fā)動機非線性模型仿真數(shù)據(jù)代替真實發(fā)動機的測量數(shù)據(jù),待估計性能參數(shù)偏差的植入值、搜索范圍,以及在各算例中使IRNRM決策計算過程收斂的非支配個體取值、健康參數(shù)估計結果和計算所需機時如表1所示。假設傳感器測量準確且不受噪聲影響,表中所列部件的健康參數(shù)偏差同時存在。另需說明的是,雖然本文在整個過渡工作過程中均采用了表1所列的植入健康參數(shù)偏差,但這并不代表在實際發(fā)動機工作過程中健康參數(shù)的偏差值是恒定的。
表1 發(fā)動機性能參數(shù)植入偏差、搜索范圍以及估計結果
從表1可以看出:采用本文所提出的方法,僅利用有限的測量傳感器布置即能實現(xiàn)大偏差范圍內(nèi)多于測量傳感器數(shù)目的大量健康參數(shù)的精確估計,而較低的計算時間成本也體現(xiàn)了本文優(yōu)化方法的計算效率。另需說明的是,本文所采用的傳感器布置(p25、T3、p3)在該過渡工作過程不同位置的序列工作點上應用時,其參數(shù)辨識能力存在一定差異,從而會影響算法的收斂性及計算效率,但是由于該部分內(nèi)容不屬于本文的討論范圍,在此不再贅述。
以圖6中算例1為例,進一步論證本文的參數(shù)估計方法:NSDE的種群大小設置為NP=100,差分變異操作中縮放因子F=0.5,交叉概率Cr=0.5。算例1在非支配解搜索過程中起始代(第0代)、中間代(第10代)以及終了代(第20代)非支配個體的分布情況如圖7所示。由于可視化呈現(xiàn)的維度限制,圖中僅選擇健康參數(shù)中壓氣機流量偏差、效率偏差兩個分量進行展示,縱坐標OF取值為式(18)中各子目標函數(shù)平方和的正二次方根。
(a) 起始代
(b) 中間代
(c) 終了代圖7 典型代數(shù)非支配解相對真實解的分布Fig.7 Distribution of nondominated solutions for representative generations relative to the true solution
從圖7可以看出:伴隨演化過程的進行,非支配解的數(shù)量逐漸增多,且在空間中向真實解聚集;根據(jù)圖4所列流程,圖7中所列非支配解均需經(jīng)歷決策過程;當某代種群非支配解集中某個體足夠接近真實解(即進入了真實解的凸鄰域范圍內(nèi))時,以其作為初值的IRNRM決策過程計算收斂并可最終獲得真實解;而圖7(c)所示終了代編號為13的非支配個體(個體取值如表1陰影列所示)即滿足了上述條件,以其作為初值進行IRNRM最優(yōu)決策計算收斂。
由決策計算中各步迭代結果的均方根誤差(RMSE)所表示的收斂歷史如圖8所示。
圖8 算例1終了代13#非支配個體IRNRM 決策過程收斂歷史Fig.8 Convergence history of No.13 nondominated individual in the 20th Generation of Case 1
(1) 本文提出的基于過渡工作過程的氣路分析方法,可以實現(xiàn)對壓氣機變導葉角度、燃油泵控制延遲等發(fā)動機關鍵機動性能參數(shù)的退化估計。
(2) SOPA方法既能解決工程實際中待求健康參數(shù)較多而測量參數(shù)不足所導致的欠定問題,又可以在最大程度上降低由多工作點分析方法的平均效應所引入的參數(shù)估計系統(tǒng)誤差。具體到本文結果,該方法能夠在僅利用3個低溫截面?zhèn)鞲衅鳒y量信息的條件下,確定發(fā)動機在不同功率條件下(針對不同算例)由15個健康參數(shù)描述的真實健康狀態(tài)。
(3) 本文提出的間接遞歸牛頓-拉夫遜法強化非支配分類差分進化算法,既具備全局搜索能力,又能夠實現(xiàn)高效高精度的局部逼近。該算法能在配置一般的個人電腦上于數(shù)小時內(nèi)完成大偏差參數(shù)空間內(nèi)健康參數(shù)的精確求解。
綜上所述,本文所提出的系統(tǒng)性方法能夠有效地解決涵蓋發(fā)動機機動性能退化的大偏差欠定氣路分析問題。