基于MATLAB的民航發(fā)動機渦輪葉片壽命建模

袁忠大，程秀全，王大偉

(1.廣州民航職業(yè)技術(shù)學(xué)院飛機維修工程學(xué)院，廣東廣州 510403；2.中國民航大學(xué)航空工程學(xué)院，天津 300300)

0 前言

在民用航空的維修工作中，發(fā)動機渦輪葉片故障導(dǎo)致的飛機臨時停場、換發(fā)現(xiàn)象層出不窮[1-2]，嚴重損害了各司的安全效益和經(jīng)濟效益。民航發(fā)動機渦輪葉片作為時壽件[3]，目前知曉其使用壽命完全基于發(fā)動機使用手冊，但手冊的編寫是基于民航發(fā)動機的材料及制造條件，而非基于用戶的使用時間、使用頻率及飛機運營當(dāng)?shù)氐臍夂驐l件。同時，發(fā)動機渦輪葉片的壽命預(yù)測直接關(guān)系到航空公司維修方案的制定及航材的儲備數(shù)量[4]，因此，各大航空公司的發(fā)動機管理中心(Engine Management Center，EMC)及可靠性管理部門都在積極對公司發(fā)動機渦輪葉片的使用數(shù)據(jù)進行統(tǒng)計分析，從而在發(fā)動機使用手冊的基礎(chǔ)上，期望依靠數(shù)理統(tǒng)計的方法切實保障公司的安全效益和經(jīng)濟效益。

國內(nèi)外學(xué)者基于數(shù)理統(tǒng)計的方法對機械電子系統(tǒng)的使用壽命進行了大量研究。蔡文斌等[5]基于三參數(shù)Weibull分布模型研究了超高強度抽油桿的概率疲勞壽命，在模型的求解過程中利用概率加權(quán)矩法估計三參數(shù)Weibull分布。韓威等人[6]研究了基于PCA和Weibull分布的滾動軸承剩余壽命預(yù)測方法，采用了極大似然估計法求解模型。王能歡等[7]提出了二、三參數(shù)Weibull分布的乘用車包修數(shù)據(jù)的分析方法。 劉建功等[8]進行了某型鋼板的疲勞試驗并得出三參數(shù)Weibull分布擬合結(jié)果相對更好的結(jié)論。WANG等[9]通過擬合軸承壽命數(shù)據(jù)，實現(xiàn)對軸承的健康管理。STRZELECKI[10]用三參數(shù)Weibull分布模型確定了不同應(yīng)力水平下的低失效概率疲勞壽命。TOASA CAIZA、 UMMENHOFER[11]基于Weibull分布進行了S355J2鋼的P-S-N曲線擬合。 ELMAHDY[12]提出了一種利用不同Weibull模型對具有失效模式的系統(tǒng)組件的壽命數(shù)據(jù)進行建模的方法。

本文作者運用三參數(shù)Weibull分布建立該型發(fā)動機渦輪葉片壽命數(shù)據(jù)的可靠性模型。為保證求解模型的計算精度，采用牛頓迭代法及三參數(shù)相關(guān)系數(shù)優(yōu)化法對該型渦輪葉片的可靠性模型進行計算[13]。同時，采用MATLAB軟件編寫計算程序并對計算結(jié)果進行K-S假設(shè)檢驗。

1 渦輪葉片壽命建模及求解

1.1 三參數(shù)Weibull分布模型

三參數(shù)Weibull分布是一種較為完善的分布，在擬合隨機數(shù)據(jù)時有很大的靈活性，對不同形狀的頻率分布有很強的適應(yīng)性，當(dāng)形狀參數(shù)取不同值時，它可以等效或接近于其他一些常用的分布。

航空產(chǎn)品的疲勞壽命和強度分布采用三參數(shù)Weibull分布可以很好地描述。但是三參數(shù)Weibull分布的參數(shù)估計比較復(fù)雜。三參數(shù)Weibull分布比二參數(shù)Weibull分布包含更多的統(tǒng)計信息，特別是對以損耗失效為特征的機械部件壽命評估中，采用三參數(shù)比采用二參數(shù)Weibull進行擬合及參數(shù)估計的精度更高，也更能反映產(chǎn)品可靠性的實際情況，因此應(yīng)用也更加廣泛。

三參數(shù)Weibull分布的概率密度函數(shù)、累積故障概率函數(shù)、可靠度函數(shù)、故障率函數(shù)、可靠壽命函數(shù)、平均故障間隔分別為

(1)

F(t)=1-e-[(t-t0)/η]m

(2)

R(t)=1-F(t)=e-[(t-t0)/η]m

(3)

(4)

tR=t0+η[-lnR(t)]1/m

(5)

θ=t0+ηΓ(1+1/m)

(6)

1.2 三參數(shù)相關(guān)系數(shù)優(yōu)化法及牛頓迭代法求解模型

改寫三參數(shù)Weibull分布的分布函數(shù)[14]為

1-F(t)=exp[-((t-t0)/η)m]

(7)

取雙對數(shù)后可得：

ln[-ln(1-F(t))]=mln(t-t0)-mlnη

(8)

作如下變換：

y=ln[-ln(1-F(t))]

(9)

x=ln(t-t0)

(10)

A=-mlnη

(11)

B=m

(12)

可得線性回歸方程

y=Bx+A=mx+A

(13)

(14)

(15)

(16)

式中：

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

因為

(25)

所以

(26)

令

(27)

(28)

將式(27)(28)代入式(26)中，得

lx0/lxx-ly0/lxy=0

(29)

式(29)較復(fù)雜，可借助數(shù)值分析計算方法，應(yīng)用計算機進行求解[15-16]。

令

E(t0)=lx0/lxx-ly0/lxy

(30)

(1)計算E(0)

(2)計算E(tmid)

2 渦輪葉片可靠性模型實際應(yīng)用

2.1 渦輪葉片壽命分布形式

以廣義Weibull分布為基礎(chǔ)，確定產(chǎn)品可靠度為R*時的壽命tR。如美國GE公司和P&W公司在計算渦輪葉片可靠性時常取R*=99.9%，允許零件的不可靠值為F=1/1 000，其含義是每1 000個零件中只允許1個零件報廢或到壽，相當(dāng)于正態(tài)分布中的-3.09σ設(shè)計標準，由此可見其安全性較高[18-19]。

2.2 渦輪葉片可靠性計算

渦輪葉片故障數(shù)據(jù)如表1所示。

表1 渦輪葉片故障數(shù)據(jù)

首先用MATLAB軟件[20]中的函數(shù)weibplot()對數(shù)據(jù)進行分析。若壽命數(shù)據(jù)來源于Weibull分布，則繪制的是直線。

采用下述命令可得到圖1所示的Weibull概率分布：

圖1 Weibull概率曲線

>>t=[5 805 5 835 5 865 5 895 5 925 5 955 5 985 6 015 6 045 6 075]；

>>weibplot(t)；

>>xlabel(‘Data{itx}’)；

>>ylabel(‘Probiblity{itP}’)；

>>title(‘Weibull Probiblity Plot’)；

由圖1可知：9個數(shù)據(jù)分布在圖中的直線附近，故渦輪葉片的壽命符合Weibull分布。

設(shè)目標函數(shù)為

E(t0)=lx0/lxx-ly0/lxy

為求出三參數(shù)Weibull分布中的位置參數(shù)，文中采用數(shù)值分析中經(jīng)典的牛頓迭代法并采用計算機分析[22]。

先編寫牛頓迭代法的M函數(shù)文件：

function x=Newt_g(f_name,x0,xmin,xmax,n_points)

clf,hold off

wid_x=xmax-xmin;dx=(xmax-xmin)/n_points;

xp=xmin:dx:xmax;

yp=feval(f_name,xp);

plot(xp,yp,′r′);

xlabel(′x′);ylabel(′f(x)′);title(′Newton Iteration′),hold on

ymin=min(yp);ymax=max(yp);wid_y=ymax-ymin;

yp=0.0*xp;

plot(xp,yp)

x=x0;

xb=x+999;

n=0;

del_x=0.001;

while abs(x-xb)>0.001

if n>300

break;

end

y=feval(f_name,x);

plot([x,x],[y,0],′black′);

plot(x,0,′o′)

fprintf(′n=%3.0f,x=%12.5ey= %12.5e ′,n ,x,y)；

if n<4

text(x,-wid_y/10,[num2str(n)]),

end

y_driv=(feval(f_name,x+del_x)-y)/del_x；

xb=x；

x=xb-y/y_driv；

n=n+1；

plot([xb，x]，[y，0]，′g′)

end

plot([x,x],[0.05*wid_y,0.2*wid_y],′r′);text(x,0.25*wid_y,′finalsolution′)

plot([x,(x-wid_x*0.004)],[0.01*wid_y,0.09*wid_y],′r′)

plot([x,(x+wid_x*0.004)],[0.01*wid_y,0.09*wid_y],′r′)

再以“ex_5x2.m”為文件名創(chuàng)建MATLAB應(yīng)用文件，最后編寫MATLAB調(diào)用程序：

x=newt_g(′ex_5x2′,5 773.8,3 000,5 804,10 000)

將它應(yīng)用于MATLAB的主窗口，運行后輸出結(jié)果如圖2所示。

圖2 位置參數(shù)輸出結(jié)果

圖3 位置參數(shù)最終結(jié)果

在應(yīng)用牛頓迭代法求解時，函數(shù)文件中有一個x0值，它為根的近似迭代初值。其選取非常關(guān)鍵，選取得不合適，將導(dǎo)致程序不能正常運行[23]。原因為此處用的是局部收斂性定理。

局部收斂性定理內(nèi)容如下：

設(shè)s是方程f(x)=0的根，在包含s的某個開區(qū)間內(nèi)f″(x)連續(xù)且f′(x)≠0，則存在δ>0，當(dāng)x0∈[s-δ，s+δ]時，有Newton法產(chǎn)生的序列{xi}收斂于s；若f″(s)≠0且x0≠s，則序列{xi}是平方收斂的。

為選取合適的x0值，選用MATLAB中的近似估計函數(shù)零點的命令格式。

[z,f-z,exitflag]=fzero(fun,x0,options,P1,P2,…)

(1)假如在fzero中直接應(yīng)用字符串表示被解函數(shù)(目標函數(shù))，容易出錯，因此先構(gòu)造內(nèi)聯(lián)函數(shù)：

(2)作圖法觀察函數(shù)零點分布

函數(shù)零點分布觀察圖源程序：

>>x=0:10:5 804;

>>y_char=vectorize(y);

>>Y=feval(y_char,x);

>>clf,plot(x,Y,′r′);holdon,plot(x,zeros (size(x)),′k′);

>>xlabel(′x′);ylabel(′y(x)′),hold off

回車輸出圖形如圖4所示。

圖4 位置參數(shù)近似初始解

(3)利用zoom和ginput指令獲得零點的初始近似值。其源程序如下：

>>zoom on%在指令窗口中運行，獲局部放大圖

>>[t u]=ginput(n)；zoom off %在指令窗口中運行，用鼠標獲得n個零點的猜測值

>>t %顯示所得零點初始猜測值

(4)求靠近t(i)的零點(即為Newton法的初始值x0)

[ti,yi,exitflag]=fzero(y,t(i),[])

通過上述步驟(1)—(4)輸出圖形和數(shù)據(jù)如圖5所示。

圖5 圖4的部分放大圖

考慮到計算中取到小數(shù)點后三位，故選取圖5所示的結(jié)果較合理。

t=5.773 8×103，u1=1.275 0×10-16，exitflag=1，exitflag>0，則表明找到近似零點退出。

例如，給定R*=99.9%，則有

渦輪葉片可靠性指標如圖6所示。

圖6 渦輪葉片可靠性指標

2.3 K-S假設(shè)檢驗

渦輪葉片壽命數(shù)據(jù)的分布類型中假設(shè)H0，即假設(shè)壽命數(shù)據(jù)符合三參數(shù)Weibull分布模型。

利用K-S檢驗法求出三參數(shù)Weibull分布下的統(tǒng)計量觀測值Dn，給定顯著性水平α=0.1，查表可得到統(tǒng)計量的臨界值D30，0.1，若Dn

表2 分布參數(shù)估計及K-S檢驗結(jié)果

由表2可知：三參數(shù)Weibull分布K-S檢驗的統(tǒng)計量觀測值小于臨界值，因此接受原假設(shè)。

3 結(jié)論

(1)運用三參數(shù)Weibull分布建立了發(fā)動機渦輪葉片可靠性壽命的模型。在數(shù)據(jù)計算過程中采用了牛頓迭代法及三參數(shù)相關(guān)系數(shù)優(yōu)化法，有效提高了模型的計算精度。

(2)對發(fā)動機渦輪葉片壽命數(shù)據(jù)的計算中，應(yīng)用MATLAB軟件實現(xiàn)了計算機處理和人工處理相結(jié)合及發(fā)動機渦輪葉片壽命數(shù)據(jù)的分析結(jié)果圖形化，極大地方便了發(fā)動機性能工程師對發(fā)動機渦輪葉片的狀態(tài)評估與監(jiān)管，提高了發(fā)動機評估預(yù)測的效率。