航空发动机限寿件疲劳可靠度计算新方法

2019-12-27游令非张建国周霜杜小松

航空学报 2019年12期

游令非，张建国，周霜，杜小松

1. 北京航空航天大学可靠性与系统工程学院, 北京 100083 2. 北京航空航天大学可靠性与环境工程技术国防重点实验室, 北京 100083

在航空发动机适航要求中，将原发失效能够引起发动机危害性影响的部件定义为发动机限寿件(Engine Life Limited Part,ELLP)，如旋转轮盘和大型旋转封严装置等，其可靠水平直接决定着发动机整体可靠性和工作性能，设计中主要通过降低ELLP的失效概率来提高整机的可靠性和安全性[1-2]。而对航空发动机造成危害性后果的事件有很大部分是由其机械结构的疲劳所导致。对压气机盘/涡轮盘进行基于疲劳可靠度计算的概率风险评估成为航空发动机型号适航取证过程以及保证其结构完整性的关键技术和重要实施步骤之一。

目前国内外针对以轮盘为代表的ELLP疲劳可靠性做了大量研究，美国军方[3]首先将概率统计理论引入到传统轮盘疲劳寿命预测方法中，形成了轮盘寿命可靠性设计方法。Melis等[4]对航空燃气轮机压气机盘的寿命进行了预测和可靠性分析。Gao等[5]应用分布式协同响应面法对航空发动机涡轮进行了疲劳损伤的可靠性分析。Zhu 等[6-7]分别建立了涡轮叶片轮盘疲劳可靠性评估的计算试验框架和概率分析框架。将超速试验与随机有限元分析相结合，量化了试验数据、材料性能和载荷的不确定性以及材料变化和荷载变化对可靠性结果的影响。Hu 等[8]针对合金涡轮盘的镍基高温合金，建立了基于外加机械功密度的概率模型，利用线性异方差函数对蠕变疲劳寿命中的非常数偏差进行评价。Wang等[9]探索了涡轮盘低周疲劳寿命预测的区域可靠性方法以及裂纹的闭合行为。李岩等[10]提出了一种基于Kriging和蒙特卡罗半径外自适应重要抽样(Monte Carlo Radius-Outside Adaptive Importance Sampling,MCROAIS)方法混合的结构概率风险评估方法，构建了压气机盘风险概率模型。陈志英等[11]结合应力应变场强法和响应面法建立了轮盘疲劳可靠性模型并进行了可靠度计算。除此以外，随机过程、神经网络和支持向量机等方法也广泛地应用到相关问题中[12-14]。

目前开展ELLP结构概率风险评估的难点在于非线性、小失效概率事件的定量要求验证问题。但是对于高度非线性极限状态函数以及小失效概率事件，以上方法均有一定的局限:① 阶矩法对于非线性较高的问题计算精度不高；② 响应面法与数值模拟相结合使计算步骤增加，特别是针对小概率失效事件计算效率低下。因此，一些基于方差抽样缩减技术的抽样方法[15-17]因具有较高的精度和鲁棒性成为解决该问题的选择。

作为在可靠度计算中广泛应用的数值模拟方法，重要抽样方法(Important Sampling, IS)[18]将抽样密度函数的抽样中心移到设计验算点(以下简称验算点)，可以使更多的样本点集中在失效域，提高了抽样的效率。但在未知验算点(或可靠度指标)的情况下，只能通过蒙特卡罗仿真(Monte Carlo Simulation,MCS)法或解析近似法求解。这极大地限制了它们的适用范围。为了克服重要抽样法对验算点和可靠度指标已知的依赖性，产生了蒙特卡罗自适应重要抽样(Monte Carlo Adaptive Important Sampling,MCAIS)法的思想。通过自动迭代不断改进MCS的重要抽样密度函数，使得抽样点越来越趋近于失效域，使抽样中心逐渐趋近于验算点。Au和Beck[19]采用基于Metropolis算法的马尔科夫链模拟方法生成失效域样本，而后运用自适应密度估计方法构造重要抽样密度函数。Sankaran和Prakasb[20]将自适应重要抽样法应用在大结构系统的可靠性分析中。黄毅等[21]采用蒙特卡罗方法对中心柱设计过程中的主要影响因素——低周疲劳进行可靠性评价。陈向前等[22]提出了一种基于样本概率密度加权的采样中心确定方法，该方法在自适应重要抽样的基础上增加了有效抽样中对失效概率贡献大的样本出现的概率。马纪明等[23]提出一种改进的自适应重要抽样方法，该方法通过失效样本估计IS密度函数的初始参数，提高了算法的收敛时间。戴鸿哲等[24]提出了一种自适应Metropolis算法和快速高斯变换技术的结构可靠性分析高效自适应重要抽样方法。但在一些情况下，自适应重要抽样法并不能得到理想的验算点近似结果，比如变量的联合概率密度函数的标准差比较大时，或者概率密度函数的梯度和极限状态函数的梯度在某一区域大致相同时等。将会使近似验算点和真实验算点间具有较大的误差。

为了克服以上问题，本文提出一种具有自更新机制的半径外自适应重要抽样(Auto Updating Monte Carlo Radius-Outside Adaptive Importance Sampling,AUMCROAIS)疲劳可靠性分析方法，本方法可以兼顾搜索效率和计算精度，在高效地寻找近似验算点的同时，能较为精确地计算可靠度。本文提出的方法总体思路如下：首先利用MCAIS快速逼近验算点附近，在MCAIS算法停止后，在以此时的近似验算点为中心进行极坐标抽样，并依次构造主动学习函数进行评判，即对极坐标抽样点先进行极限状态函数的主动学习函数的构造，根据判据自动选择那些离极限状态面近的抽样点，而后针对这些点进行最优半径的主动学习函数的构造并找到最小值点，此时，近似验算点得到了更新，而后以更新近似验算点为中心，利用半径外重要抽样计算可靠度，继续适当缩减方差并重复上述步骤，直到满足收敛条件，得到最终的可靠度。本方法通过不断的更新确定出最优抽样半径，加速失效概率计算的收敛，不受限于自适应抽样的收敛条件，保证了效率同时提高了精度。AUMCROAIS疲劳可靠性分析方法解决了疲劳可靠性问题中小概率失效和极限状态非线性较强的情况，支撑了适航规章要求(FAR33.75条款)[25]的实施，及对发动机失效概率事件分析的要求。最后在工程案例上利用本方法分别和蒙特卡罗方法、自适应抽样法和一阶可靠性方法进行对比，说明了本方法的适用性、精确性和高效性。

1 MCROAIS方法原理

图1 MCROAIS法流程图
Fig.1 Flow chart of MCROAIS method

(1)

MCROAIS虽然相比MCS提高了抽样效率，但是MCROAIS依然存在一些缺陷，表现在：① 小概率情形下且变量的标准差较大时，将会导致β在近验算点附近变化缓慢，将提前结束迭代过程；② 概率密度函数的梯度和极限状态函数的梯度在某一区域大致相同，导致β变化缓慢，将提前结束迭代过程；③ 若非线性较高时同时收敛偏离验算点，每次迭代时标准差的减小将不会使近似验算点的收敛路线得到修正，反而将加速迭代结束，得到不精确的结果；④ 若收敛准则严格(即提高精度要求)，则计算效率将大大增加。为此，在第2节中将讨论如何基于MCROAIS法加以改进并构造新的抽样方法的以同时兼顾抽样效率和计算精度。

2 基于MCROAIS法的极坐标抽样和主动学习函数

通过第1节分析可知，MCROAIS法并不能在任何极限状态函数和变量分布，特别是小概率事件下得到理想的结果，本节将通过极坐标抽样和主动学习函数的构造，使其在满足自适应抽样收敛条件下，继续构造和识别离验算点近的样本点，从而达到自更新的目的。极坐标抽样即以每次自适应迭代的结果得出的近似验算点为中心继续构造抽样，进而用主动学习函数依次对贴近极限状态和最优半径进行筛选，更新出最优点，用以构造最优半径或者下一次迭代的起始点。极坐标抽样和主动学习函数相结合的方法提高了抽样效率，同时也考虑了计算精度。

2.1 极坐标抽样

以单次自适应抽样的近似验算点为中心进行极坐标抽样，极坐标抽样基于布朗运动的轨迹概率分布，其特点是分别基于极角和极径分布进行抽样，且样本点数控制在很小的范围内，大大地提高了抽样效率同时也使样本点更多地接近近似验算点一侧。下面以二维布朗运动加以简要说明，高维同理。设x方向和y方向的位移X和Y都是独立的随机变量，符合相同的正态分布N(0,σ2)。X和Y的联合分布是fX和fY的乘积：

(2)

式中：σ2为方差。

位移R2=X2+Y2是随机变量的函数，其分布函数为

(3)

式中：r为极径；θ为极角；D是x2+y2≤r2的区域，通过变成极坐标求得结果。对式(3)求导即得到概率密度分布:

(4)

(5)

(6)

则极径r和极角θ的联合分布为

pRΘ(r,θ)=

(7)

由于无规行走的极径和极角是相互独立的随机量，故2个边缘分布可以直接求积分，则在近似验算点x*对应的近似最优半径β*下极径和极角的分布为

(8)

(9)

图2 极坐标抽样
Fig.2 Polar coordinate sampling

2.2 主动学习函数的构造

(10)

3 AUMCROAIS可靠度计算方法

AUMCROAIS法通过对最优样本的不断更新，达到对小概率事件和强非线性极限状态函数的真实验算点的近似，从而得到最优半径，继而用MCROAIS法进行可靠度求解。其包括内循环和外循环2部分。

3.1 自更新机制

外循环包括以下几个部分：以每次改进的自适应抽样的近似验算点为中心进行的极坐标抽样(pR(β*),pΘ(θ))，而后依次针对“近极限状态面”和“最优半径”进行主动学习函数的构造和判别，首先建立近极限状态面主动学习函数LG(xi)并判别找出那些近极限状态面的点xi(i=1,2,…,N)，再建立最优半径主动学习函数Lβ*(xi)并根据判别，找出近极限状态面的点中半径最小的点，即最优验算点x*，并基于更新的验算点进行半径外重要抽样可靠度计算R。最终进行可靠度收敛判别，若不满足要求，则进行下一次外循环直至满足外循环的可靠度收敛判别。基于以上步骤可同时兼顾效率和精度，同时构建双主动学习函数，给极限状态建立了统一的评判标准，有助于快速高效地识别出最优验算点。

内循环即根据本算法的特点和需要改进的自适应抽样，旨在快速使近似验算点x*收敛于真实验算点x*附近。在迭代过程中，为使近似验算点快速收敛，内循环中自适应抽样的收敛准则可适当放大，同时，由于第1次外循环已经到达真实验算点附近，所以第2次以后的循环中(如果需要)，近似验算点的收敛步长将适当减小，则：在第1次循环中σhV(v)i=σhV(v)i-1×1，hV(v)为重要抽样密度函数，第2次及以后的循环中σhV(v)i=σhV(v)i-1×0.9。同时，为加速内循环的收敛，将γ设为1.02，以此来提高内循环的计算精度。

3.2 算法流程

步骤4若没有样本点落在失效域内，令βi=β**并把hV(v)移到Y**，否则比较β*>γβi-1，若不等式成立，令βi=β**，σhV(v)i=σhV(v)i-1×0.9(k≠1时)。i=i+1，当|βi-βi-1|<0.01时，转到步骤2，直到满足要求，记录此时的Y=x*，σhV(v)k=σhV(v)i×0.9。

步骤5以x*为中心，进行极坐标抽样(pR(r),pΘ(θ))，样本点记为xi(i=1,2,…,N)。

步骤10根据下式确定第k次外循环下ELLP的失效概率:

4 应用实例

将本文提出的方法应用于某型发动机低压压气机轮盘，通过与MCS、MCROAIS法和FORM法进行对比，验证本文方法的高效率、鲁棒性和仿真精度。

针对影响发动机安全性的关键件压气机轮盘的主要失效模式—低循环疲劳断裂，本文选取某型发动机低压三级轮盘为典型关键件，开展发动机轮盘疲劳可靠性计算。

4.1 压气机轮盘疲劳试验及仿真验证

采用立式旋转试验系统开展疲劳寿命试验，试验在3 mmHg真空度下进行[10]，如图4所示。立式旋转系统控制轮盘的转速由低速500 r/min至高速9 550 r/min之间呈三角波式周期性变化，周期为30 s，如图5所示。

通过该试验可确定其疲劳危险位置，试验转速控制在500-9 550-500 r/min，共完成55 300次试验器循环，约461 h，完成114 265次试验器循环后，对该盘进行无损探伤检查，在销钉孔6点钟方向处出现了穿透性裂纹并延伸至盘体，分析可知低压三级轮盘销钉孔裂纹属于疲劳断裂，起源于销钉孔的内端面。

同时，利用有限元方法对转速为9 550 r/min的试验结果进行仿真验证。涡轮盘的可靠性主要取决于结构危险点的应力应变分布。由于结构上的形状和负载完全对称，因此本例中考虑轮盘的1/37，整体模型见图6，1/37模型见图7。为还原真实约束，仿真时添加空心销轴，轴材料为3Cr13，轮盘材料为TC11。结合结构循环对称的几何特性，对有限元模型施加以下位移约束：

1) 以发动机轮盘的轴线为z轴建立圆柱坐标系，在该圆柱坐标系内，如图8所示，约束A、B两面上全部节点沿旋转方向θ的位移。

2) 在总体笛卡尔坐标系下，约束图8中销钉两端面——C、D面，以及轮盘E、F面上全部节点沿x方向——即沿销钉轴向的位移。

叶片载荷施加在销轴上，载荷合力为29 189.6 N(这一载荷为航空发动机研究所相关部门给出的叶片产生的等效载荷)，方向为笛卡尔坐标系y方向。参数详情见表1。

图3 AUMCROAIS方法流程
Fig.3 Flowchart of AUMCROAIS method

图4 低压压气机轮盘试验组装
Fig.4 Test assembly of low pressure compressor disk

图5 500～9 550 r/min转速的变化示意图
Fig.5 Schematic diagram of change of rotation speed between 500 and 9 550 r/min

图6 低压压气机轮盘模型
Fig.6 Low pressure compressor disk model

图7 1/37轮盘和空心销模型
Fig.7 Model of 1/37 disk and hollow pin

图8 各面约束
Fig.8 All-sided constraints

选取表1中的数值作为仿真参数，利用ANSYS18.1进行仿真。模拟结果如图9所示。由结果可知应力应变水平最高处位于轮盘与销轴连接处(位置在销钉孔下部内端面，节点号1535)，此位置即结构的破坏点，最大应力为791.64 MPa。从而验证了试验的正确性和准确性，下面将用试验结果进行发动机轮盘第一级疲劳应力作用下(9 550 r/min)的疲劳可靠性计算。

表1 涡轮盘参数Table 1 Turbine disk parameters

图9 9 550 r/min转速下轮盘的应力云图
Fig.9 Stress cloud chart of disk at 9 550 r/min rotate speed

4.2 压气机轮盘疲劳可靠性分析

lgNlife=lgK-mlg(Smax-703.84)

(11)

式中：Nlife为疲劳载荷循环次数；σa为循环疲劳应力幅值；Smax为应力场强；m和K为材料参数；压气机轮盘材料为TC11，由《飞机结构金属材料力学性能手册》[26]得到其材料参数K=6.278 4×1015，m=4.736，利用应力-强度干涉理论构建模型为

K(Smax-703.84)-m-nmax

(12)

K(Smaxi-703.84)-m-nmaxi

(13)

同时，由TC11进行了材料级试验[8]，统计得到TC11材料参数m和K的均值、变异系数分别为μm=4.628、COVm=0.01，μK=6.135 26×1015，COVK=0.015，服从正态分布。

第1级应力场强的均值和变异系数为μSmax1=791.64，COVSmax1=0.1，服从极值分布。利用本文中提出的AUMCROAIS疲劳可靠性分析方法，对第一级循环应力作用下轮盘疲劳失效概率与应力循环次数的关系进行求解，同时用MCS, FORM法与MCROAIS法进行比较。其中由于MCS抽样规模较大，认为其结果是最接近真实结果的并以MCS的结果作为精确度对比。

在发动机轮盘第1级疲劳应力作用下得到疲劳失效概率Pr随载荷循环次数的关系如图10所示，载荷循环从10 000～55 000，4种方法的可靠度指标计算结果见表2，可以看出随着疲劳载荷循环次数的增加，失效概率逐渐升高。

通过将本文方法的仿真效率、鲁棒性以及仿真精度与MCS方法进行比较，从而验证本文提出方法的适用性。

1) 仿真效率：通过MCROAIS收敛速度对比，可以验证本方法的仿真效率。如图11所示，抽取12 000次、24 000次2种情况作为对比其可靠度指标β的收敛过程，可以看出本方法在抽样效率上的优越性，由于结合了极坐标抽样和主动学习函数，使近似验算点能迅速接近真实验算点x*，不断更新最优抽样半径；另外，在抽样次数上与MCS和MCROAIS进行对比，在每一级循环下，MCS需108次抽样；MCROAIS分别需32次和25次迭代，每次1 000个样本点，则MCROAIS共需32 000和25 000个样本点；而本方法则分别需要2次外循环，11次内循环； 2次外循环，8次内循环；每次内循环需抽样1 000个样本点，外循环36个样本点，则本方法分别需要11 072和8 072个样本点。同时, 在与FORM法收敛速度对比上可以看出，AUMCROAIS的收敛速度与FORM基本相同，且在一定程度上要先于FORM到达真实MPP点附近。综上所述，本方法在仿真效率上得到了很大的提升。

图10 第一级疲劳应力作用下轮盘失效概率与疲劳应力循环次数的关系
Fig.10 Relationship between failure probability of disk and number of fatigue stress cycles under the first-order fatigue stress

表2 可靠度指标结果对比Table 2 Comparisons of reliability index results

图11 可靠度指标收敛速度对比
Fig.11 Comparison of convergence rate of reliability index

2) 鲁棒性：针对高非线性的风险概率模型的失效概率仿真，验证本方法的鲁棒性。如图10所示，围绕精确解(MCS结果)，AUMCROAIS和MCROAIS均有所波动，但总体来说本文方法波动较小，更具参考性。同时，FORM法的结果随着应力循环次数的增加越来越偏离真实值，而本方法依然在小范围内波动，说明本方法将在因高非线性而导致FORM偏离真实值或者不适用的情况下依然能给出近似解。

3) 仿真精度：通过仿真结果误差分析，来验证本文方法的仿真精度。仿真精度见表2，可以看出可靠度指标总体误差不大于0.052 84，总体结果和FORM法基本相近，相比MCROAIS精度得到了很大提升。

从计算结果对比中可以看出，在不降低仿真精度和仿真效率得情况下，AUMCROAIS得到了较满意的结果，解决了在小概率失效和非线性情况下MCROAIS并不能得到稳定的精确解的问题，同时由于极坐标抽样和主动学习函数的高效的抽样和判别能力，使本方法的仿真效率大大提高，因此通过上述分析验证了本文所提方法的较高的计算效率、鲁棒性和计算精度。