南京航空航天大学学报  2017, Vol. 49 Issue (3): 389-395   PDF    
航空发动机轮盘概率风险评估方法
毕苏艺1, 孙有朝1, 李龙彪1, 陈健2, 杨坤2, 侯乃先2, 曾海军2    
1. 南京航空航天大学民航学院,南京, 211106;
2. 中国航空发动机集团商用航空发动机有限责任公司,上海, 201100
摘要: 建立了一种基于概率损伤容限的航空发动机压气机轮盘风险评估方法,考虑了缺陷发生率、缺陷分布、应力与缺陷检测概率对涡轮盘失效风险的影响。基于适航规章与咨询通告的要求,建立了考虑检查与不考虑检查两种不同的概率风险评估流程。结合有限元应力分析和断裂力学方法计算了压气机轮盘失效风险随循环变化曲线,与咨询通告AC 33.1 4-1给出的案例进行了对比,验证了本文提出方法的准确性与合理性。
关键词: 航空发动机     轮盘     裂纹扩展     蒙特卡洛仿真     风险评估    
Probabilistic Risk Assessment Method of Aeroengine Disk
BI Suyi1, SUN Youchao1, LI Longbiao1, CHEN Jian2, YANG Kun2, HOU Naixian2, ZENG Haijun2    
1. College of Civil Aviation, Nanjing University of Aeronautics & Astronautics, Nanjing, 211106, China;
2. Commercial Aircraft Engine Co., LTD, Aero Engine Corporation of China, Shanghai, 201100, China
Abstract: A risk assessment method of aeroengine compressor disk based on probabilistic damage tolerance is established. The effects of defect occurrence, defect distribution, stress and defect detecting probability on the risk of compressor disk failure are considered. The risk assessment processes with and without considering the inspection are established based on the requirements of the airworthiness regulations and advisory circular. By combining the finite element analysis and the fracture mechanics approach, the changing curve of the compressor disk failure risk versus cycle number is given. The result is compared with that in AC 33.14-1 to test and verify the accuracy and rationality of the developed method.
Key words: aeroengine     disk     crack propagation     Monte Carlo simulation     risk assessment    

随着全球航空事业的发展,航空发动机的安全性逐渐受到重视[1-2]。在航空发动机适航技术研究中,将其原发失效能够引起发动机危害性影响的部件定义为发动机寿命限制件(Engine life limited part,ELLP)[3],例如,旋转轮盘、大型旋转封严装置等,设计中主要通过降低发动机寿命限制件的失效概率来提高整机的安全性。风扇盘作为航空发动机的典型寿命限制件,在美国汽车工程师协会(Society of automotive engineers,SAE)关于发动机转子非包容性事故统计报告[4-6]以及发动机工业界记录[7]中,都有相应失效事故的记录。目前航空发动机风扇部件主要为含硬α缺陷的钛合金轮盘,此类轮盘在成型过程中存在的硬α相会导致轮盘在“安全寿命”达到之前发生破裂,因此在风扇盘设计和定寿过程中,必须评估含硬α缺陷轮盘在寿命周期内的失效概率。美国航空工业界提出采用基于概率风险评估(Probabilistic risk assessment,PRA)的部件寿命管理方法进一步降低发动机寿命限制件的失效概率,并且美国联邦航空局(Federal aviation administration,FAA)在适航规章中也提出了相关要求——在发动机的联合定义阶段通过系统安全性分析确定发动机寿命限制件后,必须通过风险评估表明发动机寿命限制件在预期使用寿命期内的失效概率风险低于10-8次/飞行小时,发动机才能获得最终的型号合格证。因此,对发动机寿命限制件在使用寿命期内的失效进行概率风险评估成为航空发动机适航取证过程的关键技术和实施步骤之一。

目前国际上关于概率风险评估研究无论在评估体系还是在基础数据获得方面都取得了很好的进展[8]。在国内,随着粉末盘在发动机上的应用,很多研究学者开展了粉末盘裂纹扩展和可靠性研究工作,徐凌志等[9]建立了考虑应力及断裂参数随机性的裂纹扩展失效功能函数,提出了亚表面裂纹扩展寿命计算方法。魏大盛等[10]在有限元分析结果的基础上,建立了针对低循环疲劳裂纹萌生和裂纹扩展失效模式的可靠度计算模型。目前国内系统研究概率风险评估的文献很少,这将极大地制约大型商用航空发动机的发展与适航取证。

AC 33.14-1[11]给出了对轮盘进行风险评估需要考虑的参数,基本给出了实施概率风险评估所需的方法与流程。本文针对航空发动机压气机轮盘概率风险评估方法进行研究,考虑了缺陷发生率、缺陷分布、应力不确定性与缺陷检测概率对轮盘失效风险的影响,基于适航规章与咨询通告的要求,建立了考虑检查与不考虑检查两种不同的概率风险评估流程,结合有限元应力分析和断裂力学方法计算了压气机轮盘失效风险随循环变化曲线,与咨询通告AC 33.14-1给出的案例进行了对比。

1 概率损伤容限评估方法

轮盘风险分析主要的不确定性是缺陷发生的小概率或者随机尺寸,形状和位置的工具损伤。在设计中,通过有效检查和去除有缺陷的轮盘降低风险。检查的不确定性包括维修次数和与检测方法有关的检测概率(Probability of detection,POD)曲线。其他不确定性还包括材料属性,负载状况以及几何形状。缺陷应力强度因子随着这些随机变量以及缺陷尺寸的变化而变化,当应力强度因子大于某一给定值时,轮盘发生断裂。

本文将轮盘划分成有限个区域,假设在每个区域有可能存在一个裂纹(由于区域存在一个裂纹的可能性已经很小,不考虑一个区域存在两个或两个以上裂纹的情况),计算缺陷发生概率Pd,利用概率损伤容限法和蒙特卡洛方法计算轮盘每个区域条件失效概率Pfd,区域失效风险Pfi(i=1, 2, 3, …, M)为

$ {P_{{\rm{f}}i}} = {P_{{\rm{fd}}}} \times {P_{\rm{d}}} $ (1)

轮盘失效风险Pf_disk

$ {P_{{\rm{f\_disk}}}} = 1-\prod\limits_{i = 1}^M {\left( {1-{P_{{\rm{f}}i}}} \right)} $ (2)

Pfi较小时,式(2) 近似为

$ {P_{{\rm{f\_disk}}}} = \sum\limits_{i = 1}^M {{P_{{\rm{f}}i}}} $ (3)
1.1 缺陷发生概率

钛合金轮盘在成型过程中会产生硬α缺陷,在一定质量W的材料上,用超越曲线描述轮盘内部的缺陷,包括每个轮盘缺陷发生率和缺陷尺寸概率分布。超越曲线绘制于对数坐标纸上,y轴表示超过一个缺陷尺寸的缺陷数量,x轴表示缺陷长度或者面积,用A表示。数值Nd[Amin]表示质量为W的材料中超越最小缺陷尺寸数量。轮盘缺陷发生率可用式(4) 表示

$ \alpha = {N_{\rm{d}}}\left[{{A_{\min }}} \right] \cdot \frac{{{W_{\rm{d}}}}}{W} $ (4)

式中Wd为轮盘质量。

假设轮盘中缺陷在任何位置发生都是等可能的,在一块定义的区域中缺陷发生率azone即为α与区域体积比例的乘积,即

$ {\alpha _{{\rm{zon}}\;{\rm{e}}}} = {N_{\rm{d}}}\left[{{A_{\min }}} \right] \cdot \frac{{{W_{\rm{d}}}}}{W} \cdot \frac{{{V_{{\rm{zon}}\;{\rm{e}}}}}}{{{V_{{\rm{disk}}}}}} $ (5)

式中:Vdisk为轮盘体积,Vzone为区域体积。式(5) 可简化为

$ {\alpha _{{\rm{zon}}\;{\rm{e}}}} = {N_{\rm{d}}}\left[{{A_{\min }}} \right] \cdot \frac{{{W_{{\rm{zon}}\;{\rm{e}}}}}}{W} $ (6)

式中Wzone为区域质量。

1.2 缺陷尺寸分布

缺陷的累积分布函数(Cumulative distribution function,CDF)为[12]

$ {\rm{CDF = 1}}-\frac{{{N_{\rm{d}}}\left( A \right)-{N_{\rm{d}}}\left( {{A_{\max }}} \right)}}{{{N_{\rm{d}}}\left( {{A_{\min }}} \right)-{N_{\rm{d}}}\left( {{A_{\max }}} \right)}} $ (7)

式中:A为缺陷尺寸(Amin < A < Amax);Nd(Amax)为超越最大尺寸Amax的缺陷数量,Nd(Amin)为超越最小尺寸Amin的缺陷数量;Nd(A)为超越尺寸A的缺陷数量。AC 33.14-1[11]给出了缺陷尺寸与缺陷数量的曲线,称为超越曲线,如图 1所示,图中横坐标为缺陷尺寸,纵坐标为超越缺陷尺寸的缺陷数量。从图中可以看出缺陷尺寸越大,超越缺陷尺寸的缺陷数量越少。在缺陷面积较小时,超越缺陷尺寸的缺陷数量随着尺寸的增加急剧减小;当缺陷尺寸超过一定数值时,超越缺陷尺寸的缺陷数量减小的趋势趋于平缓。将缺陷尺寸与超越缺陷尺寸的数量取对数后,缺陷尺寸A与超越缺陷尺寸的数量N满足关系

图 1 超越曲线 Figure 1 Exceedance curve

$ \lg N = m\lg A + n $ (8)

式中:mn为缺陷分布参数。

通过蒙特卡洛方法随机生成累计分布函数CDF值,利用式(7,8) 得到裂纹初始长度a0

$ {a_0} = \sqrt {\frac{{10\frac{{\lg \left\{ {{N_{\rm{d}}}\left( {{A_{\max }}} \right) + \left( {1- {\rm{CDF}}} \right)\left[{{N_{\rm{d}}}\left( {{A_{\min }}} \right)-{N_{\rm{d}}}\left( {{A_{\max }}} \right)} \right]} \right\} -n}}{m}}}{{\rm{\pi }}}} $ (9)
1.3 应力不确定性

为了模拟轮盘应力不确定性,假设区域应力服从均值为σ0,方差为λ2的正态分布,σ0为轮盘区域有限元分析获得的应力均值,其概率密度函数为

$ f\left( \sigma \right) = \frac{1}{{\sqrt {2{\rm{\pi }}} \lambda }}\exp \left( {-\frac{{{{\left( {\sigma-{\sigma _0}} \right)}^2}}}{{2{\lambda ^2}}}} \right) $ (10)
1.4 缺陷检测概率

缺陷检出概率与检测手段和缺陷尺寸相关。当检测手段越先进,缺陷尺寸越大,缺陷被检出的概率就越高。图 2给出了硬α缺陷检测概率曲线[11],其中横坐标表示缺陷尺寸(长度或面积),纵坐标表示该尺寸缺陷对应的检出概率。从图中可以看出,缺陷面积越大,缺陷检测概率越高,缺陷面积超过16 mm2时,检出概率几乎为1。

图 2 缺陷检测概率曲线 Figure 2 Probability of detection curve

2 概率风险评估流程

采用蒙特卡洛方法进行轮盘风险评估,针对不考虑检查与考虑检查建立两种不同的蒙特卡洛仿真模型,计算轮盘在20 000个飞行循环发生失效风险。

2.1 不考虑检查

该仿真将轮盘划分为有限个区域,每个区域指定100 000个样本。假设在每个区域只可能存在一个裂纹,结合断裂力学计算在20 000个循环发生断裂的裂纹数量,由此得到区域条件失效概率,计算每个区域的缺陷发生概率得到区域失效概率,利用式(2) 得到轮盘失效概率。

具体模拟步骤如下:

步骤1  根据有限元应力分析结果,划分轮盘区域,同一区域内应力值接近。

步骤2  定义蒙特卡洛模拟次数Ns(Ns=100 000),初始缺陷设定在区域最大应力节点位置,初始失效样本数Nf=0。

步骤3  在第i个区域随机生成初始裂纹。

步骤4  计算裂纹长度达到临界值时的循环次数Nc。对比Nc与飞行寿命Ncycle(Ncycle=20 000),如果Nc >Ncycle,则不发生失效;否则发生失效,Nf=Nf+1,直到模拟次数达到Ns

步骤5  重复步骤2~4,获得各区域内条件失效概率Pfd=Nf/Ns

步骤6  计算缺陷出现概率Pd

$ {P_{\rm{d}}} = \frac{{{N_{\rm{d}}}\left( {{a_{\min }}} \right)}}{{{{10}^6}}} \times \rho \times V $ (11)

式中:ρ为材料密度,V为区域体积,106为钛合金材料质量,Nd(amin)为106 kg钛合金超越最小尺寸缺陷数量。

步骤7  计算各区域失效概率Pf

步骤8  计算整个轮盘失效概率Pf_disk图 3给出了不考虑检查的蒙特卡洛仿真模拟流程图。

图 3 不考虑检查蒙特卡洛方法模拟转子轮盘风险评估流程图 Figure 3 Flow chart of disk risk assessment based on Monte Carlo simulation without inspection

2.2 考虑检查

引入检查间隔和检测概率曲线后,蒙特卡洛模拟仿真具体流程如下:

步骤1  根据有限元应力分析结果,划分轮盘区域,同一区域内应力值接近。

步骤2  定义模拟次数Ns(Ns=100 000),检查间隔T,初始失效样本数Nf=0,k=0。

步骤3  随机生成累积分布函数值CDF,得到裂纹初始长度a0i=1。

步骤4  随机生成检测概率POD。

步骤5  当检查间隔小于轮盘飞行寿命Ncycle时,NN=i*T;否则NN=Ncycle

步骤6  由裂纹初始长度a0,区域峰值应力σmax,判断是否失效,如果失效,Nf=Nf+1,转至步骤8;否则,计算循环到NN时的裂纹长度a,计算POD(a)。

步骤7  如果POD(a)≥POD,则返回步骤3,否则

(1) NN==Ncycle,则转至步骤8;

(2) NNNcyclei=i+1,返回步骤4。

步骤8  对下一个样本重复步骤3~7,直至Ns个样本全部模拟完毕。

步骤9  计算获得各个区域的条件失效概率Pfd=Nf/Ns

步骤10  计算缺陷出现概率Pd

步骤11  计算区域风险Pf

步骤12  计算整个轮盘失效概率Pf_disk

图 4为考虑检查间隔和检测概率曲线,运用蒙特卡洛仿真模拟转子轮盘风险评估流程图。

图 4 考虑检查蒙特卡洛方法模拟转子轮盘风险评估流程图 Figure 4 Flow chart of disk risk assessment based on Monte Carlo simulation with inspection

3 案例分析

AC 33.14-1[11]提供了概率风险评估方法的算例,描述了发动机轮盘进行了20 000次简单循环载荷,其最大转速为6 800 r/min,在外径处模拟外部叶片负荷为50 MPa。分别计算轮盘无在役检查和10 000循环在役检查的断裂概率(见图 5)。

图 5 参数设置 Figure 5 Parameter setting of the case

选取与AC 33.14-1案例同样的输入参数计算轮盘风险,将结果与AC 33.14-1进行比较,验证本文提出方法的合理性。钛合金轮盘密度为4 450 kg/m3,安全寿命为20 000飞行循环,裂纹扩展模型计算公式选取Paris公式,即

$ \frac{{{\rm{d}}a}}{{{\rm{d}}N}} = C{\left( {\Delta K} \right)^{{n_1}}} $ (12)

式中:C=9.25×10-13; n1=3.87[11]

图 1给出的缺陷分布与缺陷数量的对数坐标关系图(已将AC中关系图单位化为国际单位),采用最小二乘法拟合获得缺陷分布参数m=-1.075 5,n=-0.363 5。

图 6给出了AC33.14-1中轮盘区域划分,读取AC33.14-1有限元模型,划分轮盘区域,如图 7所示。不同颜色表示不同大小的应力,单位为MPa,红色区域应力最大,蓝色区域应力最小。箭头表示相应裂纹的应力梯度,圆点表示裂纹所处位置。区域1 ~ 4为角裂纹,区域5 ~ 16,17和20为表面裂纹,剩余区域为内埋裂纹。3种裂纹的应力强度因子计算采用的是经验公式,裂纹扩展采用的是椭圆模型(双参数),相关计算公式参见应力强度因子手册[13] (注:图 6图 7盘方位坐标系不同,将图 6右侧应力体积定义示意图顺时针旋转90°得到图 7)。

图 6 轮盘区域划分 Figure 6 Zone dividing of disk

图 7 轮盘区域划分与裂纹设置 Figure 7 Zone dividing and crack setting of disk

利用蒙特卡洛方法进行模拟,36个区域风险,轮盘总风险,单次循环风险与每个区域风险占总风险比例图如图 8~11所示。假设单次循环风险为Ppercycle,则Ppercycle可由式(13) 得到。从图 8~10可看出,随着循环数增加,区域风险、轮盘总风险和单次循环风险将增加。从图 811可看出,区域20,31,19和30等风险占轮盘总风险比例较大。轮盘不同区域风险的大小与裂纹形式及区域应力等相关,高风险裂纹(例如,角裂纹等)以及高应力导致区域风险值较大。从图 10可看出,轮盘单次循环风险在20 000循环时为1.5×10-9,适航规章要求限寿件失效发生概率低于10-7,由于一般民用涡扇航空发动机限寿件为数十个部件,因此对单个限寿件部件而言,其设定的设计目标风险值设定为10-9较为合理,本算例结果为1.5×10-9,可以认为满足限寿件失效风险的适航要求。

图 8 轮盘多区域风险 Figure 8 Multiple zone risk of disk

图 9 轮盘总风险 Figure 9 Total risk of disk

图 10 轮盘单次循环风险 Figure 10 Risk per cycle of disk

图 11 轮盘区域风险占总风险比例 Figure 11 Percentage of zone risk of disk

$ {P_{{\rm{percycle}}}} = {P_{{\rm{f\_disk}}}}/{N_{{\rm{cycle}}}} $ (13)

引入检查间隔与检测概率曲线,检查间隔为10 000飞行循环,检测概率曲线如图 12所示。将上述输入条件代入仿真模型,36个区域风险、轮盘总风险、单次循环风险以及各区域风险占总风险比例如图 13~16所示。从图 13~15可看出,随着循环增加,区域风险,轮盘总风险以及单次循环风险逐渐增加,在10 000次循环时采取检查措施后,区域风险、轮盘总风险及单次循环风险先减小,然后随循环增加,采取检查措施有效降低了轮盘风险。结合图 1316可看出,区域29,20,31,21和30风险占轮盘总风险比例较大,20 000次循环时轮盘单次循环风险为1.44×10-9,满足适航规章要求。

图 12 检测概率曲线 Figure 12 Probability of detection curve

图 13 轮盘多区域风险 Figure 13 Multiple zone risk of disk

图 14 轮盘总风险 Figure 14 Total risk of disk

图 15 轮盘单次循环风险 Figure 15 Risk per cycle of disk

图 16 轮盘多区域风险所占比例 Figure 16 Percentage of zone risk of disk

AC33.14-1以轮盘单次循环风险作为验证风险评估流程合理性的指标,多次进行蒙特卡洛仿真,计算所得轮盘风险与AC33.14-1案例计算值进行比较如表 1所示。由表 1可知,本文计算结果在AC33.14-1允许值范围内,完成了风险评估流程正确性的验证工作。

表 1 计算结果与AC33.14-1对比 Table 1 Comparison result between results and AC 33.14-1

4 结论

本文针对航空发动机限寿件概率风险评估方法开展了研究,根据发动机限寿件的适航条款要求,建立了概率风险评估方法和流程,综合有限元应力分析、断裂力学分析、无损检测和概率分析,评估了在役检查中发动机转子轮盘发生断裂的风险,建立了与AC 33.14-1案例相同的轮盘模型,利用本文建立的风险评估流程计算了轮盘风险,研究发现:

(1) 轮盘服役过程中未采取检查措施总风险约为3×10-5,单次循环风险范围约为1.54×10-9 ~ 1.6×10-9,满足适航规章要求。

(2) 采取检查措施后轮盘风险降低,轮盘风险约为2.9×10-5,单次循环风险约为1.26×10-9 ~ 1.3×10-9,满足适航规章要求。

(3) 在10 000次循环实施了检查,有效降低了10 000次循环后每次循环风险值。10 000次循环后,单次循环风险先随循环减少,然后随循环增加,在20 000次循环时其单次循环风险低于10 000次循环时的单次循环风险。

计算结果满足AC 33.14-1要求,轮盘区域风险、轮盘总风险、单次循环风险变化趋势符合定性分析规律,本文建立的方法通过了考题验证,可以为航空发动机轮盘概率风险评估提供参考。

参考文献
[1] 金捷, 刘邓欢. 航空发动机燃烧室湍流两相燃烧模型发展现状[J]. 南京航空航天大学学报, 2016, 48(3): 303–309.
JIN Jie, LIU Denghuan. Recent advances in turbulent two-phase combustion models[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2016, 48(3): 303–309.
[2] 曾海军, 孙有朝, 李龙彪. 航空发动机风扇叶片振动适航符合性设计方法[J]. 南京航空航天大学学报, 2015, 47(6): 884–889.
ZENG Haijun, SUN Youchao, LI Longbiao. Approach for vibration airworthiness compliance of aeroengine fan blade[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2015, 47(6): 884–889.
[3] U.S. Department of Transportation, Federal Aviation Administration. Advisory circular 33.70-1: Guidance material for aircraft engine life-limited parts requirements [S].Washington, DC: FAA, AC 33.70-1, 2009.
[4] SAE International Group. Report on aircraft engine containment[R]. SAE AIR 1537, 1977.
[5] SAE International Group. Report on aircraft engine containment[R]. SAE AIR 4003, 1991.
[6] SAE International Group. Report on aircraft engine containment[R]. SAE AIR 1537A, 1996.
[7] CORRAN R, GORELIK M, LEHMANN D, et al. The development of anomaly distributions for machined holes in aircraft engine rotors[C]// ASME Turbo Expo 2006: Power for Land, Sea, and Air.[S.l.]:ASME, 2006:941-950.
[8] KAPPAS J. Review of risk and reliability methods for aircraft gas turbine engines[R]. Fishermans Bend, Victoria3207 Australia: DSTO Aeronautical and Maritime Research Laboratory, DSTO-TR-1306, 2002.
[9] 徐凌志, 吕文林. 粉末冶金涡轮盘裂纹扩展失效概率分析[J]. 机械科学与技术, 2000, 19(2): 210–212.
XU Lingzhi, LV Wenlin. Crack growth failure probability of PM turbine disk[J]. Mechanical Science and Technology, 2000, 19(2): 210–212.
[10] 魏大盛, 杨晓光, 王延荣. 基于缺陷分布形式的粉末冶金涡轮盘可靠度计算模型[J]. 机械工程学报, 2008, 44(11): 132–137.
WEI Dasheng, YANG Xiaoguang, WANG Yanrong. Model for calculating the reliability of powder metallurgy turbine disk based on distribution of defects[J]. Chinese Journal of Mechanical Engineering, 2008, 44(11): 132–137.
[11] U.S. Department of Transportation, Federal Aviation Administration. Advisory circular 33.14-1: damage tolerance for high energy turbine engine rotors[S]. Washington, DC: FAA, AC 33.14-1, 2001.
[12] ENRIGHT M P, HUDAK S J, MCCLUNG R C, et al. Application of probabilistic fracture mechanics to prognosis of aircraft engine components[J]. AIAA Journal, 2006, 44(2): 311–316. DOI:10.2514/1.13142
[13] 中国航空研究院. 应力强度因子手册[M]. 北京: 北京科学出版社, 1993.