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

袁忠大，程秀全，王大偉

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

0 前言

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

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

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

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

1.1 三參數Weibull分布模型

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

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

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

(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 三參數相關系數優化法及牛頓迭代法求解模型

改寫三參數Weibull分布的分布函數[14]為

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

(7)

取雙對數后可得：

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)較復雜，可借助數值分析計算方法，應用計算機進行求解[15-16]。

令

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

(30)

(1)計算E(0)

(2)計算E(tmid)

2 渦輪葉片可靠性模型實際應用

2.1 渦輪葉片壽命分布形式

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

2.2 渦輪葉片可靠性計算

渦輪葉片故障數據如表1所示。

表1 渦輪葉片故障數據

首先用MATLAB軟件[20]中的函數weibplot()對數據進行分析。若壽命數據來源于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個數據分布在圖中的直線附近，故渦輪葉片的壽命符合Weibull分布。

設目標函數為

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

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

先編寫牛頓迭代法的M函數文件：

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”為文件名創建MATLAB應用文件，最后編寫MATLAB調用程序：

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

將它應用于MATLAB的主窗口，運行后輸出結果如圖2所示。

圖2 位置參數輸出結果

圖3 位置參數最終結果

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

局部收斂性定理內容如下：

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

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

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

(1)假如在fzero中直接應用字符串表示被解函數(目標函數)，容易出錯，因此先構造內聯函數：

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

函數零點分布觀察圖源程序：

>>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 位置參數近似初始解

(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)輸出圖形和數據如圖5所示。

圖5 圖4的部分放大圖

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

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

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

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

圖6 渦輪葉片可靠性指標

2.3 K-S假設檢驗

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

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

表2 分布參數估計及K-S檢驗結果

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

3 結論

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

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