基于粒子群極值Kriging模型的壓氣機(jī)葉盤時變可靠性分析

王遠(yuǎn)卓, 韓治國, 楊子煜, 路成, 薛小鋒

1.北京航空航天大學(xué) 宇航學(xué)院, 北京100191; 2.西北工業(yè)大學(xué) 航天學(xué)院, 陜西 西安710072;3.復(fù)旦大學(xué) 航空航天系, 上海200433; 4.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安710072

低壓壓氣機(jī)轉(zhuǎn)子系統(tǒng)對于航空發(fā)動機(jī)的正常運轉(zhuǎn)起到尤為重要的作用，其功用是通過壓縮氣體來提升整個系統(tǒng)的熱力循環(huán)效率。低壓壓氣機(jī)葉盤工作環(huán)境嚴(yán)苛，承受著高壓、高轉(zhuǎn)速等交變載荷作用[1]，一旦發(fā)生功能失常，將對航空發(fā)動機(jī)產(chǎn)生二次損傷，造成災(zāi)難性后果。因此，為了保障航空發(fā)動機(jī)系統(tǒng)安全運行，有必要開展低壓壓氣機(jī)葉盤時變可靠性分析。

近年來，國內(nèi)外學(xué)者基于蒙特卡洛(Monte Carlo,MC)模擬與近似解析法開展了許多研究，對航空發(fā)動機(jī)壓氣機(jī)葉盤進(jìn)行了可靠性分析。趙洪利等[2]基于MC模擬，給出了用于評估航空發(fā)動機(jī)發(fā)生故障可能性的風(fēng)險預(yù)測方法；Narayanan等[3]提出了基于MC模擬的航空發(fā)動機(jī)結(jié)構(gòu)疲勞壽命估算方法；楊周等[4]根據(jù)模擬的隨機(jī)參數(shù)與隨機(jī)響應(yīng)的關(guān)系，通過一次二階矩法對發(fā)動機(jī)進(jìn)行可靠性分析；張飛[5]對一次二階矩法在均值點進(jìn)行線性化展開的基礎(chǔ)上進(jìn)行改進(jìn)，提出了通過在驗算點進(jìn)行線性化對發(fā)動機(jī)葉盤進(jìn)行可靠性分析的改進(jìn)一次二階矩方法；白冰等[6]針對強(qiáng)非線性問題高精度可靠性求解困難的問題，依據(jù)二次二階矩建立了可靠度指標(biāo)分析方法。上述研究工作為航空發(fā)動機(jī)壓氣機(jī)葉盤結(jié)構(gòu)可靠性分析提供了思路，但MC模擬需要通過大規(guī)模計算來實現(xiàn)可靠性分析，耗時較長且不能高效完成具有動態(tài)時變特性的可靠性評估；近似解析法適用于功能函數(shù)已知的可靠性分析，對于壓氣機(jī)葉盤這類復(fù)雜結(jié)構(gòu)可靠性分析，其功能函數(shù)通常是未知的。為合理規(guī)避MC模擬與近似解析法存在的缺點，許多專家學(xué)者通過構(gòu)建代理模型進(jìn)行結(jié)構(gòu)可靠性分析。高陽等[7]提出了徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)模型，對高壓渦輪盤疲勞可靠性分析進(jìn)行了研究；劉佳奇等[8]基于智能神經(jīng)網(wǎng)絡(luò)提出了一種航空發(fā)動機(jī)運行狀態(tài)下的時變安全性分析方法；馬小駿等[9]基于最小二乘支持向量機(jī)算法對航空發(fā)動機(jī)壽命進(jìn)行了分析。上述代理模型雖然一定程度上提高了可靠性分析的效率，但是對于具有時變特性的可靠性分析問題，需要建立大量的代理模型獲取結(jié)構(gòu)可靠性指標(biāo)，其分析流程過于復(fù)雜。因此，部分學(xué)者將極值思想引入代理模型來簡化可靠性計算的負(fù)擔(dān)。Zhang等[10-11]為有效地進(jìn)行渦輪葉盤疲勞可靠性及中心齒輪可靠性分析，基于響應(yīng)面方法，借鑒極值思想提出了多重極值響應(yīng)面策略； Bai等[12]針對航空發(fā)動機(jī)葉盤的非確定性因素提出了一種基于動態(tài)子結(jié)構(gòu)的極值響應(yīng)面法。然而，上述方法在構(gòu)建功能函數(shù)過程中采用的是基于二次多項式的形式，其可靠性分析精度與工程需求存在偏差。

為了解決上述問題，本文基于Kriging模型，結(jié)合極值思想和粒子群算法(PSO)，提出了粒子群極值Kriging模型方法(PSOEKM)，用以實現(xiàn)航空發(fā)動機(jī)低壓壓氣機(jī)葉盤時變可靠性分析。其中，極值思想用以簡化動態(tài)分析流程，以輸出響應(yīng)的極值代替整個時變過程；PSO算法用以替代梯度下降，實現(xiàn)Kriging模型超參數(shù)尋優(yōu)；Kriging模型用于建立輸出響應(yīng)與輸入?yún)?shù)之間關(guān)系，進(jìn)而實現(xiàn)可靠性分析。此外，通過多種方法對比驗證了本文研究方法的有效性與可行性。

1 粒子群極值Kriging模型

1.1 分析流程

為了合理高效實現(xiàn)航空發(fā)動機(jī)低壓壓氣機(jī)葉盤時變可靠性分析，將極值思想和PSO算法引入Kriging模型，提出PSOEKM方法，其原理為：結(jié)合極值思想處理壓氣機(jī)葉盤輸出響應(yīng)的時變過程，運用輸出響應(yīng)極值將動態(tài)分析轉(zhuǎn)化為瞬態(tài)分析；采用PSO算法實現(xiàn)目標(biāo)函數(shù)求解，獲取Kriging模型的超參數(shù)；利用Kriging模型描述輸出響應(yīng)與輸入?yún)?shù)之間的關(guān)系，構(gòu)建功能函數(shù)并實現(xiàn)時變可靠性分析?；赑SOEKM方法的航空發(fā)動機(jī)低壓壓氣機(jī)葉盤時變可靠性分析流程如圖1所示。

圖1 基于PSOEKM方法的航空發(fā)動機(jī)低壓壓氣機(jī)葉盤時變可靠性分析流程

由圖1可知，基于PSOEKM方法的壓氣機(jī)葉盤動態(tài)時變可靠性分析包括確定性分析、樣本獲取、模型構(gòu)建、時變可靠性分析4個部分。具體分析流程如下：

1) 構(gòu)建航空發(fā)動機(jī)低壓壓氣機(jī)葉盤三維有限元模型，設(shè)置材料參數(shù)、邊界條件與動態(tài)載荷；執(zhí)行動態(tài)確定性分析；獲取輸出響應(yīng)隨時間變化的規(guī)律并結(jié)合極值思想確定輸出響應(yīng)極值與研究時刻點。

2) 基于輸入變量隨機(jī)分布特征，獲取足夠的輸入樣本，通過確定性分析計算各輸入樣本對應(yīng)的輸出響應(yīng)值；確定用以建模的訓(xùn)練樣本和用以模型驗證的測試樣本。

3) 結(jié)合訓(xùn)練樣本，運用PSO算法實現(xiàn)模型超參數(shù)尋優(yōu)，并求解PSOEKM的模型系數(shù)，構(gòu)建壓氣機(jī)葉盤輸出響應(yīng)與輸入?yún)?shù)之間的關(guān)系模型。其中，基于PSO算法的模型超參數(shù)求解步驟為：參數(shù)初始化，適應(yīng)度函數(shù)值計算，更新粒子的位置與速度、個體極值與群體極值，完成迭代獲取模型超參數(shù)。

4) 在建立的壓氣機(jī)葉盤輸出響應(yīng)模型基礎(chǔ)上，構(gòu)建其極限狀態(tài)函數(shù)，采用MC方法進(jìn)行大量模擬實現(xiàn)壓氣機(jī)葉盤時變可靠性分析并輸出結(jié)果。

1.2 數(shù)學(xué)模型

參考Kriging模型,PSOEKM模型y(x)可描述為[13]

y(x)=GT(x)α+z(x)

(1)

式中:y(x)為輸出響應(yīng)值;α為待定系數(shù)向量;z(x)為誤差修正函數(shù);G(x)為基函數(shù)矩陣。

與傳統(tǒng)的響應(yīng)面法相比,Kriging模型的優(yōu)勢主要體現(xiàn)在誤差修正項z(x),該項的引入可進(jìn)一步提升了模型變量之間的關(guān)聯(lián)程度,進(jìn)而提升了代理模型的精度。通常,誤差修正函數(shù)z(x)服從均值為0、方差為σ2的正態(tài)分布[14],即

(2)

式中:E為函數(shù)均值;var為函數(shù)方差;cov為函數(shù)的協(xié)方差;R為相關(guān)函數(shù)矩陣,p,q=1, 2, …,m;m為訓(xùn)練樣本的數(shù)量;xp,xq為第p和q個訓(xùn)練樣本。

xp與xq的相關(guān)函數(shù)為

(3)

式中:i=1, 2, …,n,n為訓(xùn)練樣本的維度;θi為第i個模型的超參數(shù)。

相關(guān)函數(shù)常用的形式有指數(shù)型、高斯型及線性型等,本文選取高斯型函數(shù)作為相關(guān)函數(shù)的形式,原因在于該形式函數(shù)具有較好的擬合性能[15],高斯型相關(guān)函數(shù)(3)式可表達(dá)為

(4)

因此,基于樣本的相關(guān)函數(shù)矩陣R為

R=

(5)

進(jìn)而,模型中的待定系數(shù)向量α可由(6)式確定

α=(GT(x)R-1G(x))-1GT(x)R-1Ytrain

(6)

式中,Ytrain為訓(xùn)練樣本的輸出向量。

則誤差修正函數(shù)為

z(x)=rT(x)R-1(Y-GTα)

(7)

式中,rT(x)表示單個樣本與所有樣本相關(guān)性大小。

rT=[(R(θ,x,x1)R(θ,x,x2…R(θ,x,xm)]。

依據(jù)最小二乘理論,模型誤差的方差估計值為

(8)

模型超參數(shù)傳統(tǒng)上采用梯度下降法,通過最大化適應(yīng)度函數(shù)進(jìn)行求解獲取。然而,對于涉及多變量參數(shù)高非線性問題,梯度下降策略尋優(yōu)效率下降,甚至無法得到全局最優(yōu)解。因此,本文將最大化問題轉(zhuǎn)化為最小化問題[16]

(9)

式中,f(θ)為適應(yīng)度函數(shù)。

為了實現(xiàn)最小化問題求解,采用粒子群算法通過求解(9)式實現(xiàn)模型超參數(shù)θ的尋優(yōu)。粒子群算法通過粒子的隨機(jī)性、記憶性、不確定性,改善了基于梯度下降模型參數(shù)尋優(yōu)容易陷入局部最優(yōu)的問題。相比遺傳算法,PSO算法在更新過程中實時追蹤當(dāng)前最優(yōu)粒子,故收斂速度較快[17],可以較好、較快地滿足時變條件下的可靠性分析。粒子群算法粒子的速度和位置更新如(10)式所示

(10)

式中:k為迭代次數(shù);v(k)為當(dāng)前粒子速度;x(k)為當(dāng)前粒子位置;c1和c2為學(xué)習(xí)因子;β為隨機(jī)產(chǎn)生的隨機(jī)數(shù);xpbest為個體最優(yōu)值;xgbest為全局最優(yōu)值。

基于上述分析,結(jié)合訓(xùn)練樣本運用粒子群算法可實現(xiàn)模型超參數(shù)尋優(yōu),進(jìn)而結(jié)合最小二乘法實現(xiàn)模型系數(shù)求解,建立PSOEKM模型(功能函數(shù))。

1.3 時變可靠性分析

基于所構(gòu)建的功能函數(shù),結(jié)合分析目標(biāo)的許用值,可建立發(fā)動機(jī)壓氣機(jī)葉盤徑向變形的極限狀態(tài)函數(shù),其可表達(dá)為[18]

(11)

在此基礎(chǔ)上,運用MC法對極限狀態(tài)函數(shù)進(jìn)行大量抽樣,利用概率統(tǒng)計分析實現(xiàn)時變可靠性分析,其原理實質(zhì)是借鑒伯努利大數(shù)定律,即

(12)

式中:nB為n次試驗下事件B出現(xiàn)的次數(shù);P為事件B出現(xiàn)的概率;ε為任意很小的正數(shù)。

根據(jù)(12)式,當(dāng)試驗次數(shù)n趨近于無限大時,n次試驗中滿足可靠性閾值要求的次數(shù)將收斂到實際的可靠度,則可靠度計算原理為

(13)

式中:R(t)為可靠度函數(shù);pr(t)為可靠密度函數(shù);r(t)為0～t時間段內(nèi)累積的故障數(shù);N0為t=0時規(guī)定條件下進(jìn)行工作的產(chǎn)品數(shù)。

2 算例分析

2.1 動態(tài)確定性分析

壓氣機(jī)工作環(huán)境復(fù)雜,葉盤作為壓氣機(jī)重要的部件之一,其可靠性直接關(guān)乎航空發(fā)動機(jī)整個系統(tǒng)的安全性。選取某型航空發(fā)動機(jī)低壓壓氣機(jī)葉盤作為研究對象,考慮流體載荷和結(jié)構(gòu)載荷的作用,對其徑向變形進(jìn)行動態(tài)可靠性分析。為了實現(xiàn)低壓壓氣機(jī)葉盤徑向變形動態(tài)可靠性分析,首先建立研究對象的三維模型(如圖2所示),結(jié)合構(gòu)建的三維模型,進(jìn)而生成低壓壓氣機(jī)葉盤的有限元模型(如圖3所示),航空發(fā)動機(jī)低壓壓氣機(jī)葉盤由四面體單元組成,包括150 837個節(jié)點和89 424個單元。

圖2 低壓壓氣機(jī)葉盤三維模型

圖3 低壓壓氣機(jī)葉盤有限元模型

選取1Cr11Ni2W2MoV作為低壓壓氣機(jī)葉盤材料,其密度為7.8×103kg/m3、泊松比為0.3、彈性模量為2.06×1011Pa。此外,考慮進(jìn)口流速、出口壓力、密度和轉(zhuǎn)速等作用,其中進(jìn)口流速為124 m/s,出口壓力為5.07×105Pa,以及轉(zhuǎn)速為1 168 rad/s。需要說明的是:轉(zhuǎn)速是隨時間變化的,其變化特性參考文獻(xiàn)[16]確定。在有限元模型和載荷參數(shù)設(shè)置的基礎(chǔ)上,結(jié)合有限元法和有限元體積法對壓氣機(jī)葉盤徑向變形進(jìn)行動態(tài)確定性分析。在動態(tài)確定性分析過程中,分別針對流場和結(jié)構(gòu)場進(jìn)行設(shè)置:在流場中設(shè)置進(jìn)口流速和出口壓力載荷;在結(jié)構(gòu)場中,設(shè)置材料參數(shù)、轉(zhuǎn)速。最后對低壓壓氣機(jī)葉盤設(shè)置圓柱支撐約束?；诖_定性分析得到低壓壓氣機(jī)葉盤徑向變形隨時間變化曲線如圖4所示。

圖4 低壓壓氣機(jī)葉盤徑向變形隨時間的變化曲線

由圖4可知,低壓壓氣機(jī)葉盤徑向變形的最大值出現(xiàn)在分析時域[165 s, 200 s],任意選取其中某一時刻作為研究時刻點,本文選取t=172 s,則此時的徑向變形分布云圖如圖5所示。由圖5可以看出,低壓壓氣機(jī)葉盤徑向變形的最大位置位于葉尖尖端,同時將該部位作為時變可靠性分析研究的部位。

圖5 低壓壓氣機(jī)葉盤徑向變形分布云圖

2.2 樣本獲取

基于低壓壓氣機(jī)葉盤徑向變形動態(tài)確定性分析,結(jié)合其最大值出現(xiàn)的時刻和位置,選取進(jìn)口流速v、出口壓力pout、密度ρ和轉(zhuǎn)速w作為隨機(jī)輸入變量,假設(shè)這些變量參數(shù)均服從正態(tài)分布且相互獨立,其數(shù)值分布特征如表1所示[19-20]。

依據(jù)表1輸入變量的隨機(jī)分布特征,運用拉丁超立方抽樣獲取隨機(jī)輸入變量的樣本,并結(jié)合多次動態(tài)確定性分析計算輸入樣本對應(yīng)的輸出響應(yīng)極值(最大徑向變形),進(jìn)而形成含有輸入變量和輸出響應(yīng)的樣本。其中,隨機(jī)選取一部分樣本作為訓(xùn)練樣本,用以構(gòu)建低壓壓氣機(jī)葉盤徑向變形的PSOEKM模型(功能函數(shù));剩余的樣本用以驗證PSOEKM模型的精度。

表1 輸入變量的隨機(jī)分布特征

2.3 模型構(gòu)建

基于動態(tài)確定性分析獲取的訓(xùn)練樣本,結(jié)合1.2節(jié)數(shù)學(xué)模型構(gòu)建原理,運用粒子群算法實現(xiàn)PSOEKM模型超參數(shù)尋優(yōu)。其中,粒子群算法尋優(yōu)的初始化條件設(shè)置參數(shù)為:學(xué)習(xí)因子c1=0.56,c2=1.95、最大迭代次數(shù)為25、種群規(guī)模為20、進(jìn)化速度范圍為[-13,13]、進(jìn)化位置范圍為[0.1,20]。基于粒子群算法的適應(yīng)度函數(shù)進(jìn)化曲線如圖6所示。

圖6 基于粒子群算法的適應(yīng)度函數(shù)進(jìn)化曲線

從圖6可以看出,整個粒子適應(yīng)度函數(shù)收斂的過程中,開始收斂速度較快,隨著迭代次數(shù)的增加,收斂速度變慢直至穩(wěn)定。群體粒子適應(yīng)度通過7次迭代基本收斂,在第12次迭代時目標(biāo)值達(dá)到最小,穩(wěn)態(tài)值為2.9×10-4。此時,PSOEKM模型超參數(shù)θ=[0.716,13.500,0.521,0.223]。

在得到模型超參數(shù)的基礎(chǔ)上,為構(gòu)建低壓壓氣機(jī)葉盤徑向變形的功能函數(shù),采用最小二乘法實現(xiàn)PSOEKM模型相關(guān)函數(shù)矩陣R與待定系數(shù)α向量的求解

(14)

2.4 可靠性分析

基于建立的功能函數(shù),結(jié)合低壓壓氣機(jī)葉盤徑向變形的許用值,其極限狀態(tài)函數(shù)根據(jù)(11)式可表達(dá)為

(15)

利用MC法對(15)式進(jìn)行10 000次仿真抽樣,則低壓壓氣機(jī)葉盤徑向變形的仿真歷史與分布直方圖分別如圖7～8所示。由圖7和圖8可知,低壓壓氣機(jī)葉盤徑向變形服從正態(tài)分布。當(dāng)?shù)蛪簤簹鈾C(jī)葉盤徑向變形許用值為1.594×10-3m(該值基于3σ原則確定)時[16],可靠度為99.76%,滿足工程需求。

圖7 低壓壓氣機(jī)葉盤徑向變形仿真歷史

圖8 低壓壓氣機(jī)葉盤徑向變形分布直方圖

3 粒子群極值Kriging模型驗證

3.1 建模特性

為了說明PSOEKM在建模特性方面的優(yōu)勢,基于測試樣本,通過均方根誤差ERMS和平均相對誤差Ear對其建模精度進(jìn)行研究,其計算原理見(16)式。在PSOEKM建模特性過程中,通過極值響應(yīng)面法(extremum response surface method,ERSM)和極值Kriging模型(extremum Kriging model,EKM)進(jìn)行對比分析。需要說明的是:所有的計算在相同的計算環(huán)境下執(zhí)行;ERSM和EKM模型均是采用相同的訓(xùn)練樣本構(gòu)建;在建模精度分析過程中,以低壓壓氣機(jī)葉盤徑向變形動態(tài)確定性分析結(jié)果作為參考;在建模效率方面,則是以ERSM建模時間作為參考。建模精度與效率分析結(jié)果如表2所示。

(16)

表2 建模精度與效率分析結(jié)果

由表2可以看出,PSOEKM的建模時間相對于EKM和ERSM略有增長,其原因是PSO在模型超參數(shù)尋優(yōu)過程中存在多次迭代計算,但是其建模效率是可接受的。在建模精度方面,相較于EKM和ERSM,所提出的PSOEKM具有一定的優(yōu)勢,PSOEKM的均方根誤差和平均絕對誤差小于ERSM和EKM的均方根誤差和平均絕對誤差,并且PSOEKM的建模精度相對于ERSM和EKM分別提高了4.62%和4.38%。因此,在建模效率可接受范圍內(nèi),PSOEKM在建模精度方面具有明顯的優(yōu)勢。

3.2 仿真性能

PSOEKM方法的仿真性能主要從仿真效率和分析精度兩方面進(jìn)行研究,結(jié)合多種方法對比分析的手段加以說明,其中所用的方法有直接模擬法、ERSM和EKM?；诒?中輸入變量數(shù)值隨機(jī)的分布特征,分別將4種方法用于執(zhí)行不同頻次(102,103和104)的MC抽樣,實現(xiàn)低壓壓氣機(jī)葉盤徑向變形時變可靠性分析。其中,由于計算限制沒有執(zhí)行104次直接模擬。此外,在仿真性能分析過程中,以直接模擬的分析結(jié)果作為參考?；贓RSM、EKM和PSOEKM的仿真性能分析結(jié)果如表3～5所示。

表3 仿真效率分析結(jié)果 s

表4 仿真可靠性分析結(jié)果

表5 仿真精度分析結(jié)果 %

由表3可知,代理模型(包括PSOEKM、EKM和ERSM)的仿真時間遠(yuǎn)遠(yuǎn)小于直接模擬的仿真時間,由表4～5可以看出,在不同的仿真模擬次數(shù)下,PSOEKM的仿真結(jié)果優(yōu)于ERSM和EKM的仿真結(jié)果,并且其分析結(jié)果與直接模擬的結(jié)果最為接近;相較于ERSM和EKM,PSOEKM的仿真精度分別提升了0.61%和0.11%。因此,PSOEKM在仿真特性方面具有一定的優(yōu)勢。

4 結(jié) 論

基于Kriging模型,結(jié)合PSO與極值思想,提出了PSOEKM方法,在考慮流體載荷和結(jié)構(gòu)載荷交互作用的情況下,本文實現(xiàn)航空發(fā)動機(jī)低壓壓氣機(jī)葉盤徑向變形時變可靠性分析,并通過方法對比驗證了所提出方法的有效性和可行性。本文主要結(jié)論如下:

1) 考慮多種時變載荷的影響,開展了低壓壓氣機(jī)葉盤徑向變形動態(tài)確定性分析,獲取輸出響應(yīng)隨時間變化規(guī)律,并確定了輸出響應(yīng)極值與研究時刻點;

2) 基于PSOEKM實現(xiàn)了低壓壓氣機(jī)葉盤徑向變形建模,通過大量仿真模擬,當(dāng)徑向變形許用值為1.594×10-3m時,其可靠概率為99.76%;

3) PSOEKM在建模精度方面具有一定的優(yōu)勢,其精度相對于ERSM和EKM分別提高了4.62%和4.38%;

4) 在不同模擬次數(shù)下,在仿真效率滿足工程需求的前提下,PSOEKM仿真精度均優(yōu)于ERSM和EKM。