基于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软件实现了计算机处理和人工处理相结合及发动机涡轮叶片寿命数据的分析结果图形化，极大地方便了发动机性能工程师对发动机涡轮叶片的状态评估与监管，提高了发动机评估预测的效率。