航空发动机限寿件概率失效风险评估的等效应力转化

丁水汀，郭祎玮，2，李 果，3，周 煜

(1.北京航空航天大学能源与动力工程学院航空发动机气动热力国家级重点实验室，北京100191；2.中国航发商用航空发动机有限责任公司，上海200241；3.多伦多大学航空航天研究所，多伦多M3H5T6；4.北京航空航天大学交通科学与工程学院，北京100191)

1 引言

通常，因失效可能引起航空发动机危害性事故的部件定义为航空发动机限寿件[1]，包括发动机的主要转子结构件和部分静子结构件，如各类旋转轮盘和承压机匣等。航空发动机适航规章FAR/CCAR-33.70条款及AC 33.14、AC 33.70-1要求，限寿件必须开展寿命期内的概率失效风险评估工作。开展概率风险评估可在限寿件设计生产阶段进行，预测其是否满足适航规章[2]对失效概率的要求。

有关限寿件概率失效风险的研究，国外从评估方法到基础数据库积累均取得了良好进展。美国联邦航空局(FAA)与西南研究院等科研机构和GE、普惠等航空发动机制造商合作，共同开发了成熟的概率失效风险分析方法。该方法起初由Leverant[3]、McClung[4]等提出，目标是解决钛合金压气机盘材料Hard α夹杂缺陷引起的失效问题。之后，Momin[5]和Wu[6]等进一步研究了分析流程，分别从断裂力学模型建立、随机变量分布及无损探伤检查等方面对失效风险分析方法进行了最终的完善，形成了达尔文软件[3-6]，从而建立了一整套完整的分析理论及流程。

受制于技术封锁和经验缺乏，国内对限寿件概率失效风险分析的研究目前尚处于起步阶段。丁水汀等[7]通过对适航条款及国内外相关研究的分析，对概率失效风险评估方法进行了研究；赵合阳等[8]利用GH4133材料低循环疲劳试验数据，对涡轮盘低循环疲劳失效进行了概率寿命分析；杜少辉等[9]在裂纹扩展安全寿命的基础上，对轮盘的失效概率进行了分析。

目前，概率失效风险评估主要通过断裂力学来分析部件的失效[10]，利用应力-强度干涉理论，并考虑材料参数、缺陷分布及载荷等变量的随机性，对失效概率进行定量分析[3，11]。其中，多场耦合的应力分布结果是失效风险分析的关键输入条件。通常，限寿件在发动机飞行循环中承受的应力是随时间变化的瞬态结果，但在失效风险分析中要求应力只能以一个载荷对的形式来执行断裂力学计算[2]。因此，对于限寿件应力分析，必需将飞行循环内的变化应力转化为断裂力学计算所需的一个载荷对。为此，断裂力学计算时，通常，选取区域内最大应力作为整个区域的等效应力[12]，即基于局部的等效应力转化方法。但这一方法在计算确定寿命期内的概率失效风险时，得到的结果通常偏保守。

本文以某型航空发动机离心压气机盘这一典型限寿件为例，提出一种将区域内所有节点在飞行循环内的应力转化为等效应力的方法，即基于全局的应力转化方法。该方法考虑了瞬态情况，且尽可能反映载荷在过渡过程中的变化，更贴近发动机真实运转情况，可获得更为准确的失效风险评估结果。

2 瞬态载荷计算

2.1 计算方法

采用单向流固耦合[13-15]方法对轮盘类限寿件进行温度及应力计算，以获得压气机盘随时间变化的应力载荷分布。即先采用CFX将瞬态边界条件代入分析流程计算瞬态温度载荷，随后将温度结果以载荷步的形式导入ANSYS模块进行静力分析，获得瞬态应力分布。

2.2 研究模型

以某型航空发动机离心压气机轮盘为研究对象(图1(a))，建立流体域模型(图1(b))。为减少分析中的计算量，取轮盘的1/10建立模型(图1(c))。

图1 研究模型Fig.1 Research model

2.3 有限元网格

图2为轮盘子午面网格节点分布。如图所示，概率失效风险评估建立在多场耦合获取的轮盘载荷分布的基础之上，因此选用轮盘子午面载荷分布来表征轮盘载荷分布结果。同时，分网时对有限元网格加以约束，即该截面上网格节点与流体网格具有完整的节点对应关系，以满足后续概率失效风险评估需要。轮盘子午面网格节点共1 300个。

图2 含完整节点的轮盘子午面Fig.2 Meridian plane of disc with complete nodes

2.4 边界条件和材料属性

针对分析模型，由空气系统程序获得典型飞行循环瞬态工况的边界条件。由于本文的主要目的是给出等效应力转化方法，因此对边界条件进行了简化以模拟发动机的工作过程，见图3。以图3(a)中转速随时间的变化为例，随着发动机起动，压气机转速逐渐增加，经历一段时间波动后趋于稳定，最终工作一段时间后转速下降。整个过程历时10.0 s，其他边界条件与转速相匹配，基本模拟了发动机起动-巡航-减速-停机这4个主要工作状态。需指出，该方法可扩展至全部飞行循环。

压气机盘的材料为钛合金TC4，其各项属性见表1[16]。表中，a为热膨胀系数，E为弹性模量，ρ为密度，λ为导热系数，μ为泊松比。分析中为简化起见，不考虑轮盘材料各属性随温度的变化；流体域介质为理想可压缩气体，且边界无滑移，壁面光滑[17]。

2.5 瞬态温度场计算结果

利用CFX软件对已绘制好的有限元网格加入相应的瞬态边界条件进行计算。由于模型取自整盘的1/10，对流体和固体截面加入周期性对称边界条件。为得到较好的瞬态收敛结果，以稳定段边界计算出的稳态结果作为初值[18-19]。获得轮盘工作过程瞬态温度场后，将CFX计算结果按时间步导出节点温度数据。同时利用2.3节轮盘子午面的节点及单元和节点的对应关系，通过MATLAB编制程序，实现每一时间点对应温度数据的自动导入，并绘制各时刻子午面温度分布云图，见图4。图中,横纵坐标轴代表子午面上节点坐标，不同颜色表示对应节点温度值，下文类似。图中选取了初始(0 s)、应力最大(0.5 s)、应力波动(1.0 s)、稳态(4.0 s和 8.0s)及终止(10.0 s)6个最有代表性的时刻，来表征轮盘工作过程中载荷随时间的变化。

图3 边界条件随时间的变化Fig.3 Boundary conditions change as a function of time

表1 TC4材料属性Table 1 Material properties of TC4

图4 轮盘子午面温度分布云图Fig.4 Temperature distribution nephogram on meridian surface of compressor disk

2.6 瞬态应力场计算结果

应力场采用热固间接耦合分析方法[20]计算，其中间接法体现在计算时按照载荷步循环重复进行应力场计算。本研究中，每一时刻的应力场计算均是基于先期获得的温度场作为热边界条件施加在模型上，并通过MATLAB程序实现循环。图5给出了各时刻子午面应力分布的有限元云图。图6为在发动机一个飞行循环中，一级压气机盘最大应力点处应力随时间的变化。

图5 轮盘子午面应力分布云图Fig.5 Stress distribution nephogram on meridian surface of compressor disk

图6 一级压气机盘最大应力点处应力随时间的变化Fig.6 Stress changes with time at peak stress

3 全局等效应力转化

本文参考传统安全寿命理论[21-23]，考虑发动机整个飞行循环，在等损伤前提下，提出并采用了一种基于全局的应力转化方法获得概率失效风险评估的应力输入条件。其具体流程如图7所示。

图7 全局应力转化流程Fig.7 Global stress transformation process

3.1 雨流计数法及线性累计损伤理论

利用安全寿命理论进行低周疲劳应力变化分析时，通常需要利用循环计数法将轮盘不规则的应力随时间的变化转化为若干个应力循环对。目前工程中使用最广泛并为美国FAA认可的循环计数法是雨流计数法[24]。

利用雨流计数法获得与瞬态载荷等效的一组载荷对后，通过查寻对应材料的S-N曲线可得到每个载荷对对应的循环寿命Ni。由于计数后每个载荷对的应力比不同，需要利用Goodman方程进行统一：

式中：σa为任意实际应力循环对的应力幅值，σm为平均应力，σb为材料静强度。

根据线性累计损伤理论，若构件在恒幅交变载荷Si作用下对应循环次数ni，则疲劳损伤Di为：

对于每个转化后的等幅载荷，若破坏准则为：

可以认为，若每一个等幅载荷单独加载的疲劳寿命之和D1与一个载荷D经过相同多次加载后的疲劳寿命D2相同，则载荷可以等效。

最后，反查S-N曲线即可得到等效的应力载荷对[21-23]。至此，通过线性累计损伤理论即可把雨流计数后的一组载荷对转化为可用于断裂力学计算的一个载荷。

3.2 转化过程

对压气机盘进行瞬态分析时，需要对盘-盘腔模型每个节点在整个发动机飞行循环中的应力变化进行雨流计数处理。为提高效率，通过C语言编制“四点法雨流计数”程序[25-26]完成转换工作。该程序主要由两部分组成：①简化数据，剔除不会对材料造成疲劳损伤的、相邻的、应力值相同的点，提取瞬态应力的波峰和波谷；②依照四点法编写雨流计数程序，其计算结果为每个点所受的与瞬态应力等效的恒幅载荷，如表2所示。表中，恒幅载荷最大应力幅SI和恒幅载荷最小应力幅SI′构成了Si。根据材料手册确定应力转化中所使用的S-N曲线，且材料屈服极限为973 MPa，应力比为-1。S-N曲线可转化为如图8中便于程序计算的标准幂函数形式[27]。图9为全局应力转化程序框图。首先将雨流计数法得到的一组载荷对利用Goodman方程统一为应力比为-1后，再导入S-N曲线得到对应的寿命，进而利用线性累计损伤理论计算全局等效应力，最后通过程序得到整个飞行循环轮盘的全局等效应力分布。

表2 每个点雨流计数程序运行结果Table 2 Rain-flow program operation results

图8 幂函数拟合S-N曲线Fig.8 Power function fittingS-Ncurve

图9 全局应力转化程序框图Fig.9 Global stress transformation program diagram

3.3 转化结果

图10给出了经转化的典型飞行循环下压气机轮盘子午面的全局和局部等效应力的分布云图。可看出，与最大应力时刻的局部等效应力相比，全局等效应力的应力水平显著下降，但分布规律相似。

4 概率失效风险分析计算

4.1 计算方法

图10 两种等效应力分布云图Fig.10 Two types of equivalent stress distribution nephograms

根据文献[28]中的概率失效风险计算方法，采用分区评估方式，并结合应力-强度干涉原理、线弹性断裂力学及概率统计建立分析模型；应用Monte Carlo法、Runge-Kutta法和Newman法，确定模型的数值解法[29-31]。在此基础上，可以将概率失效风险评估的流程，编制为失效风险的计算程序。其中，限寿件的应力分布将影响断裂力学计算中裂纹扩展驱动力的大小从而影响轮盘失效风险，限寿件的温度分布变化会影响材料的裂纹扩展速率和断裂韧性，如图11所示。因此，当失效概率在某一应力水平(图中为500 MPa)之下时，可认为不发生失效；当应力超过这一水平时，失效概率增长率随应力的增加而增加，直到完全失效为止。

图11 失效概率随应力和温度分布的变化Fig.11 The change of failure probability with the stress and temperature distribution

4.2 计算结果

将全局等效应力分布与基于局部的应力分布分别代入概率失效风险分析流程，得到两种不同应力转化方法下的压气机盘失效概率，如图12所示。结果表明，使用循环数达到20 000时，基于局部应力转化方法获得的失效概率为2.2 10-11次/飞行循环，而基于全局应力转化的为4.8 10-12次/飞行循环。这表明，与基于局部应力转化的方法相比，采用FAA所推荐的基于全局应力转化方法具有更小的失效概率。尽管这一过程需要后期的验证性试验数据支持，但上述典型示例表明转换方法不同对失效概率有着显著且重要的影响。

图12 不同应力转化方法下的压气机盘失效概率Fig.12 Failure probability of compressor disk under two different stress transformation methods

5 结论

通过对概率失效风险评估输入条件的两类转化方法进行对比，在相同飞行循环(20 000)内，采用全局应力转化方法获得的等效应力分布作为输入条件代入概率失效风险评估流程时获得的失效概率更小(全局 4.8 10-12次/飞行循环，局部 2.2 10-11次/飞行循环)，表明不同等效应力转化方法对概率失效风险有显著影响。同时，也揭示出基于局部应力转化方法获得的失效概率是否过于保守，是我国民机型号未来适航取证过程中需仔细考虑的问题。

本文的研究成果对概率失效风险评估中关键输入参数应力分布数据的确定给出了定量参考，后续可以通过协调行业内相关发动机研究所，以真实压气机盘为分析对象，获取工作状态下轮盘各项边界条件的变化数据，以失效风险分析的数值计算结果对试验结果进行评估。