文章快速检索  
  高级检索
航空发动机限寿件疲劳可靠度计算新方法
游令非1,2, 张建国1,2, 周霜1,2, 杜小松1,2     
1. 北京航空航天大学 可靠性与系统工程学院, 北京 100083;
2. 北京航空航天大学 可靠性与环境工程技术国防重点实验室, 北京 100083
摘要: 针对目前的航空发动机限寿件(ELLP)疲劳可靠性分析中的小失效概率事件以及其极限状态函数具有较强非线性的特点,提出了一种具有自更新机制的半径外自适应重要抽样(AUMCROAIS)疲劳可靠性分析方法。该方法首先利用蒙特卡罗自适应重要抽样(MCAIS)快速逼近真实设计验算点(MPP)附近,随后以近似设计验算点为中心进行极坐标抽样,并依次构造主动学习函数,对近极限状态函数和抽样半径进行最优选取,从而实现最优抽样半径的更新,通过不断的更新确定出最优抽样半径,加速失效概率计算的收敛。本方法提高了设计验算点的收敛速度同时保证了计算精度,解决了小失效概率事件以及强非线性极限状态函数可靠度计算难题,最后以某型发动机压气机轮盘为对象应用本方法,并与传统的蒙特卡罗仿真(MCS)方法、蒙特卡罗半径外自适应重要抽样法(MCROAIS)和一阶可靠性方法(FORM)进行了对比,验证了本方法的高效率、鲁棒性和仿真精度。
关键词: 疲劳可靠性     航空发动机     自适应重要抽样     自更新     主动学习    
A new method for fatigue reliability calculation of aero-engine limited life parts
YOU Lingfei1,2, ZHANG Jianguo1,2, ZHOU Shuang1,2, DU Xiaosong1,2     
1. School of Reliability and System Engineering, Beihang University, Beijing 100083, China;
2. Science and Technology on Reliability and Engineering Laboratory, Beihang University, Beijing 100083, China
Abstract: Aiming at the small failure probability events and the strong non-linearity of limit state function in the fatigue reliability analysis of aero-engine Limited Life Parts (ELLP), a fatigue reliability analysis method based on the Auto Updating Monte Carlo Radius-Outside Adaptive Importance Sampling (AUMCROAIS) is proposed. In this method, the Monte Carlo Adaptive Importance Sampling (MCAIS) is used to approach the Most Probable Point (MPP) efficiently, then the polar coordinate sampling is employed by taking the approximate MPP as the sampling center. An active learning function is constructed to optimize the near-limit state function and sampling radius, so that the optimal sampling radius can be updated. The optimal sampling radius can be determined through continuous updating, accelerating the convergence of failure probability. This method improves the convergence speed of MPP points and ensures the accuracy of calculation. It solves the problem of reliability calculation of small failure probability events and strong non-linear limit state function. Finally, taking a compressor disk of an engine as an application, and the efficiency, robustness and simulation accuracy of the proposed method are verified by comparing with the traditional Monte Carlo Simulation (MCS) method, the Monte Carlo Radius-Outside Adaptive Importance Sampling (MCROAIS), and First-Order Reliability Method (FORM).
Keywords: fatigue reliability     aero-engine     adaptive importance sampling     automatic updating     active learning    

在航空发动机适航要求中,将原发失效能够引起发动机危害性影响的部件定义为发动机限寿件(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方法原理

蒙特卡罗半径外自适应重要抽样法(MCROAIS)是基于MCAIS而来。MCROAIS可以通过不断调整可靠度指标半径(β球),使其逐渐接近可靠度指标。其通过搜寻随机变量概率密度fX(·)值最大点而不断修正。设极限状态函数为g(X),其中X=[x1  x2  …  xn],MCROAIS估计失效概率${\widehat P_{\rm{f}}}$的计算流程如图 1所示,图中:k为迭代次数;β为可靠度指标;γ为迭代点到均值点的距离;r1r2为半径。

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

其中更新机制如下:若xmin存在,则βk+1=β(xmin);若xmax存在,则比较β(xmax)>γβk, γ=1.02,若不等式成立,βk+1=β(xmax);若不成立,进入下一次循环。估计失效概率${\widehat P_{\rm{f}}}$的计算公式为

$ {{\hat P}_{\rm{f}}} = \frac{1}{N}\sum\limits_{k = 1}^n {\sum\limits_{i = {n_{k - 1}} - 1}^{{n_k}} {I\left[ {g\left( {_k{v_i}} \right)} \right]\left( {1 - \chi _n^2\left( {{\beta _k}} \right)} \right)} } $ (1)
 

式中:N为总抽样次数;nk为第k次迭代步骤里重要抽样所需样本点数;vi为MCROIS中的抽样点,i为序列数;χn2(·)为自由度为n的卡方分布;I[·]为示性函数。

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

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

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

2.1 极坐标抽样

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

$ {f_{XY}}\left( {x,y} \right) = \frac{1}{{2{\rm{ \mathsf{ π} }}{\sigma ^2}}}\exp \left( { - \frac{{{x^2} + {y^2}}}{{2{\sigma ^2}}}} \right) $ (2)
 

式中:σ2为方差。

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

$ \begin{gathered} {F_R}\left( r \right) = \iint\limits_D {{f_{XY}}\left( {x,y} \right){\text{d}}x{\text{d}}y} = \hfill \\ \;\;\;\;\;\;\;\;\;\int_0^{2{{\rm{\pi }}}} {{\text{d}}\theta \int_0^R {\frac{1}{{2{{\rm{\pi }}}{\sigma ^2}}}\exp \left( { - \frac{{{r^2}}}{{2{\sigma ^2}}}} \right)r{\text{d}}r} } = \hfill \\ \;\;\;\;\;\;\;\;\;1 - \exp \left( { - \frac{{{R^2}}}{{2{\sigma ^2}}}} \right) \hfill \\ \end{gathered} $ (3)
 

式中:r为极径;θ为极角;Dx2+y2r2的区域,通过变成极坐标求得结果。对式(3)求导即得到概率密度分布:

$ {f_R}\left( r \right) = \frac{r}{{{\sigma ^2}}}\exp \left( { - \frac{{{r^2}}}{{2{\sigma ^2}}}} \right) $ (4)
 

R的期望是$\sqrt {\frac{{({\rm{ \mathsf{ π} }}\sigma )}}{2}} $>0。对于高维的具有类似的推导过程。当期望μX=μY=0,方差都是σ2 (各向同性)时,这时其极坐标极径r和极角θ的分布分别为一个Rayleigh分布和一个均匀分布。考虑位于极坐标中点(μr, μθ)处的无规行走,即其直角坐标xy的期望值为μXμY,方差都是σ2 (各向同性)。于是

$ {\mu _r} = \sqrt {\mu _X^2 + \mu _Y^2} $ (5)
 
$ {\mu _\theta } = \left\{ {\begin{array}{*{20}{l}} {\arctan \frac{{{\mu _Y}}}{{{\mu _X}}}}&{x > 0}\\ {\arctan \frac{{{\mu _Y}}}{{{\mu _X}}} + {\rm{ \mathsf{ π} }}}&{x < 0} \end{array}} \right. $ (6)
 

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

$ \begin{array}{l} {p_{R\mathit{\Theta }}}\left( {r,\theta } \right) = \\ \;\;\;\;\;\;\;\frac{r}{{2{\rm{ \mathsf{ π} }}{\sigma ^2}}}\exp \left( { - \frac{{{r^2} + \mu _r^2 - 2r{\mu _r}\cos \left( {\theta - {\mu _\theta }} \right)}}{{2{\sigma ^2}}}} \right) \end{array} $ (7)
 

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

$ \begin{array}{l} {p_R}\left( {{\beta ^*}} \right) = \int_0^{2{\rm{ \mathsf{ π} }}} {{p_{R\mathit{\Theta }}}} {\rm{d}}\theta = \\ \;\;\;\;\;\;\;\;\;\;\;\;\frac{{{\beta ^*}}}{{{\sigma ^2}}}{I_0}\left( {\frac{{{\beta ^*}{\mu _{{\beta ^*}}}}}{{{\sigma ^2}}}} \right)\exp \left( { - \frac{{{\beta ^{ * 2}} + \mu _{{\beta ^ * }}^2}}{{2{\sigma ^2}}}} \right) \end{array} $ (8)
 
$ \begin{array}{l} {p_\mathit{\Theta }}\left( \theta \right) = \int_0^\infty {{p_{R\mathit{\Theta }}}} {\rm{d}}{\beta ^ * } = \frac{1}{{\sqrt {8{\rm{ \mathsf{ π} }}{\sigma ^2}} }}{\mu _{{\beta ^ * }}}\cos \left( {\theta - {\mu _\theta }} \right) \cdot \\ \;\;\;\;\;\;\;\;\exp \left( { - \frac{{{\mu _{{\beta ^ * }}}{{\sin }^2}\left( {\theta - {\mu _\theta }} \right)}}{{2{\sigma ^2}}}} \right) \cdot \\ \;\;\;\;\;\;\;\;\left[ { - {\rm{erf}}\left( {\frac{{{\mu _{{\beta ^ * }}}\cos \left( {\theta - {\mu _\theta }} \right)}}{{\sqrt {2{\sigma ^2}} }}} \right) + 1} \right] + \\ \;\;\;\;\;\;\;\;\frac{1}{{2{\rm{ \mathsf{ π} }}}}\exp \left( { - \frac{{{\mu _{{\beta ^ * }}}^2}}{{2{\sigma ^2}}}} \right) \end{array} $ (9)
 

式中:I0(x)是零阶第一类Bessel函数,I0(x)= $\;\frac{1}{{\rm{ \mathsf{ π} }}}\int_0^{\rm{ \mathsf{ π} }} {\exp (x\cos \theta )d} \theta $;erf(x)是误差函数,erf(x)= $\frac{2}{{\sqrt {\rm{ \mathsf{ π} }} }}\int_0^{\rm{ \mathsf{ π} }} {{e^{ - {t^2}}}{\rm{d}}t} $。二维极坐标抽样示例如图 2所示,图中实心五角星(极限状态面上)为实际验算点,实心五角星(失效域内)为本次自适应迭代得出的近似验算点,空心五角星为按极坐标抽样得出的新的样本点,可以看出极坐标抽样不需要过多的样本点即可覆盖以β*为极径的区域,且越接近极限状态面越密,提高了抽样效率。

图 2 极坐标抽样 Fig. 2 Polar coordinate sampling
2.2 主动学习函数的构造

为使得样本点更多落在抽样半径确定的球区域附近,从而加速失效概率计算的收敛,本文构造主动学习函数针对极坐标抽样的结果来识别最优样本点,使得目标函数$\widehat G$(xi)尽可能接近极限状态函数的阈值a附近,范围为[aζ, a+ζ],ζ为误差。当$\widehat G$(xi)>a时,表示产品处于失效状态;当$\widehat G$(xi)<a时,表示产品处于安全状态;当$\widehat G$(xi)=a时,表示极限状态,a通常取值为0。本文构造了双主动学习函数,即先对“近极限状态函数面”进行选择,若不满足,则基于上一次近似验算点重新进行自适应抽样;若选择出近极限状态面的点,则再在这些点中,对“最优半径β* ”进行选择。主动学习函数为

$ \begin{array}{l} L\left( {{\mathit{\boldsymbol{x}}_i}} \right) = \left( {\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right) - a} \right) \times \\ \;\;\;\;\;\;\;\;\left[ {2\mathit{\Phi }\left( {\frac{{a - \hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}{{{\sigma _{\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}}}} \right) - \mathit{\Phi }\left( {\frac{{a - \zeta - \hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}{{{\sigma _{\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}}}} \right) - } \right.\\ \;\;\;\;\;\;\;\;\left. {\mathit{\Phi }\left( {\frac{{a + \zeta - \hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}{{{\sigma _{\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}}}} \right)} \right] - {\sigma _{\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}\left[ {2\varphi \left( {\frac{{a - \hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}{{{\sigma _{\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}}}} \right) - } \right.\\ \;\;\;\;\;\;\;\;\left. {\varphi \left( {\frac{{a - \zeta - \hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}{{{\sigma _{\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}}}} \right) - \varphi \left( {\frac{{a + \zeta - \hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}{{{\sigma _{\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}}}} \right)} \right] + \\ \;\;\;\;\;\;\;\;\left[ {\mathit{\Phi }\left( {\frac{{a - \hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}{{{\sigma _{\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}}}} \right) - \mathit{\Phi }\left( {\frac{{a - \zeta - \hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}{{{\sigma _{\hat G\left( {{\mathit{\boldsymbol{x}}_i}} \right)}}}}} \right)} \right] \end{array} $ (10)
 

式中:$\widehat G$(xi)和σ$\widehat G$(xi)2分别为极坐标抽样样本点的极限状态函数或概率密度函数的预测值和方差;Φ(·)为标准正态累积分布函数;φ为标准正态概率密度函数;误差ζ可取max($\widehat G$(xi))×0.1(i=1, 2, …, N);xi为第i个样本点。找出所有样本点主动学习函数值的最大值, 即max(L(xi))(i=1, 2, …, N)。

构造学习函数后,建立判别准则,即max(L(xi))<ε,ε为计算误差。本文中,先进行“近极限状态”判别,$\widehat G$(xi)和σ$\widehat G$(xi)2在这里分别为极坐标抽样样本点的极限状态函数的预测值和方差,均选取max($\widehat G$(xi))×0.1(i=1, 2, …, N)并找出满足L(xi)>0.000 1的样本点,当不满足判别条件时,说明没有近极限状态面,则需更新样本点重新进行极坐标抽样,继续迭代;当满足判别条件时,当前所得所对应的样本点xi(i=1, 2, …, M) (M为近极限状态样本个数),继续进行二次主动学习函数构造和判别,主动学习函数中的$\widehat G$σ2在这里分别为极坐标抽样样本点对应的最优半径β*的预测值和方差,分别记为$\widehat \beta $*(xi)和σ$\widehat \beta $*(xi)2,均选取max($\widehat \beta $*(xi))×0.1(i=1, 2, …, M),找出max(L(xi))(i=1, 2, …, M)即是最优样本点x*,用于构造抽样半径。

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 算法流程

步骤1  设定初始球半径βk,步骤k=1。将原随机变量X转化为标准正态空间的随机变量Y=[y1  y2  …  yn],相应的状态函数转换为g(Y)。选取初始迭代点,一般为随机变量均值点μY,初始化β=0, γ=1.02。

步骤2  若为步骤10跳回,重置i=1;按重要抽样密度函数hV(v)抽样获得nk个样本点Yj=[y1j  y2j  …  ynj](j=1, 2, …, nk),其中若为第1次迭代(k=1),σhV(v)=σY

步骤3  对于在失效域的样本点,计算fX(kYj),记录max(fX(kYj)),将其样本设为Y*β*=║Y*2;对于不在失效域的样本点,计算g(Yj),记录min(g(Yj)),将其样本设为Y**β**=║Y**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)i=σhV(v)i-1×0.9。

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

步骤6  基于xi(i=1, 2, …, N),建立近极限状态面主动学习函数LG(xi),进行近极限状态判据,取σ$\widehat G$ (xi)ζ均为max($\widehat G$ (xi))×0.1(i=1, 2, …, N),找出(L(xi))>0.000 1的样本点,当不满足判别条件时,则需更新样本点,重新进行极坐标抽样,继续迭代;当满足判别条件时,记样本点为xi(i=1, 2, …, M)。

步骤7  当前所得所对应的样本点构造最优半径主动学习函数Lβ*(xi),取σ$\widehat \beta $ *(xi)2ζ$\widehat \beta $*均为max($\widehat \beta $(xi))×0.1(i=1, 2, …, M),找出max (L(xi))(i=1, 2, …, M)即更新最优样本点x*

步骤8  抽取半径随机样本。根据设计验算点x*得到最优抽样半径r* = $\sqrt {\sum {{{(\mathit{\boldsymbol{x}}_j^*)}^2}} } $,确定抽样区域(r*, r*+ξ),其中,ξ≤3,抽取服从卡方分布χ2(n)的半径随机样本点R,并判断所产生的随机样本是否落入抽样区域,若是,则将该样本点记为Rj(j=1, 2, …, NNMCROIS),NNMCROIS为要取得半径样本点的总个数。

步骤9  抽取落在抽样区域的随机样本Yj。随机样本Yj=[yj, 1  yj, 2  …  yj, n],令aj= $\frac{{\mathit{\boldsymbol{Y}}{'_j}}}{{\sqrt {\sum {{{(y{'_{i, j}})}^2}} } }}$,则Yj=Rjaj,利用随机样本点Yj计算极限状态函数预测值$\widehat G$ (Yj)并得到I[$\widehat G$ (Yj)]的值:如果$\widehat G$ (Yj)<0,则I[$\widehat G$ (Yj)]=1;否则,I[$\widehat G$ (Yj)]=0。

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

$ {{\hat P}_{{\rm{f}},k}} = \left( {1 - \chi _n^2\left( {\beta _k^2} \right)} \right)\frac{1}{{{N_{{\rm{NMCROIS}}}}}}\sum\limits_{j = 1}^{{N_{{\rm{NMCROIS}}}}} {I\left[ {\hat G\left( {{\mathit{\boldsymbol{Y}}_j}} \right)} \right]} $

式中:βk=║r*║,当${\widehat P_{{\rm{f}}, n}} - {\widehat P_{{\rm{f}}, n - 1}}$ >0.001时,k=k+1,σhV(v)k=σhV(v)k-1×0.9,回到步骤2,直到满足要求。算法流程图如图 3所示。

图 3 AUMCROAIS方法流程 Fig. 3 Flowchart of AUMCROAIS method
4 应用实例

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

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

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

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

图 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

通过该试验可确定其疲劳危险位置,试验转速控制在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。结合结构循环对称的几何特性,对有限元模型施加以下位移约束:

图 6 低压压气机轮盘模型 Fig. 6 Low pressure compressor disk model
图 7 1/37轮盘和空心销模型 Fig. 7 Model of 1/37 disk and hollow pin

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

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

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

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

表 1 涡轮盘参数 Table 1 Turbine disk parameters
变量 物理意义 数值
E1/GPa 弹性模量(轮盘) 123
ε1 泊松比(轮盘) 0.33
ρ1/(g·cm-3) 密度(轮盘) 4.48
E2/GPa 弹性模量(销轴) 219
ε2 泊松比(销轴) 0.3
ρ2/(g·cm-3) 密度(销轴) 7.76
ω/(r·min-1) 转速 9 550
F/N 全部载荷合力 29 189.6

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

图 9 9 550 r/min转速下轮盘的应力云图 Fig. 9 Stress cloud chart of disk at 9 550 r/min rotate speed
4.2 压气机轮盘疲劳可靠性分析

压气机轮盘应力比RL=0的S-N曲线(即压气机轮盘在0~29 189.6 N的脉动载荷下的疲劳强度与疲劳寿命曲线)可由Nlife=am得到,则

$ \lg {N_{{\rm{life}}}} = \lg K - m\lg \left( {{S_{\max }} - 703.84} \right) $ (11)
 

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

$ \begin{array}{l} G\left( {n,\mathit{\boldsymbol{Z}}} \right)m = {{\bar D}_{\rm{f}}} - {{\bar D}_{\rm{f}}}\left( t \right) = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;K{\left( {{S_{\max }} - 703.84} \right)^{ - m}} - {n_{\max }} \end{array} $ (12)
 

式中:Z=[K  m  Smax];Df为疲劳退化强度;Df(t)为疲劳载荷作用下的退化量;nmax为疲劳应力Smax对应的循环载荷次数。考虑单个疲劳应力水平下的模型

$ \begin{array}{l} G\left( {n,\mathit{\boldsymbol{Z}}} \right)m = {{\bar D}_{\rm{f}}} - {{\bar D}_{\rm{f}}}\left( t \right) = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;K{\left( {{S_{\max i}} - 703.84} \right)^{ - m}} - {n_{\max i}} \end{array} $ (13)
 

同时,由TC11进行了材料级试验[8],统计得到TC11材料参数mK的均值、变异系数分别为μ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,可以看出随着疲劳载荷循环次数的增加,失效概率逐渐升高。

图 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
转速/(r·min-1) MCS AUMCROAIS AUMCROAIS误差 MCROAIS MCROAIS误差 FORM FORM误差
10 000 3.027 95 3.044 4 0.016 45 3.160 1 0.132 15 3.011 20 0.016 75
14 000 2.872 03 2.883 9 0.011 87 2.895 4 0.023 37 2.855 34 0.016 69
18 000 2.752 81 2.732 1 0.020 71 2.612 7 0.140 11 2.738 80 0.014 01
22 000 2.658 63 2.679 8 0.021 17 2.630 0 0.028 63 2.645 68 0.012 95
26 000 2.578 73 2.597 0 0.018 27 2.605 0 0.026 27 2.568 10 0.010 63
30 000 2.516 62 2.512 2 0.004 42 2.448 1 0.068 52 2.501 62 0.015 00
34 000 2.456 75 2.473 9 0.017 15 2.524 6 0.067 85 2.443 44 0.013 31
38 000 2.405 16 2.458 0 0.052 84 2.358 4 0.046 76 2.391 71 0.013 45
42 000 2.358 39 2.373 2 0.014 81 2.274 7 0.083 69 2.345 15 0.013 24
46 000 2.315 63 2.352 0 0.036 37 2.289 3 0.026 33 2.302 82 0.012 81
50 000 2.276 58 2.280 1 0.003 52 2.248 8 0.027 78 2.264 01 0.012 57
55 000 2.233 98 2.211 0 0.022 98 2.198 3 0.035 68 2.219 63 0.014 35

通过将本文方法的仿真效率、鲁棒性以及仿真精度与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点附近。综上所述,本方法在仿真效率上得到了很大的提升。

图 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并不能得到稳定的精确解的问题,同时由于极坐标抽样和主动学习函数的高效的抽样和判别能力,使本方法的仿真效率大大提高,因此通过上述分析验证了本文所提方法的较高的计算效率、鲁棒性和计算精度。

5 结论

本文针对小概率失效事件和强非线性可靠性分析问题,综合考虑了抽样效率和抽样精度,针对航空发动机结构疲劳可靠性分析问题,提出了AUMCROAIS方法,经具体案例验证表明:

1) 本文提出的AUMCROAIS方法实施方便,计算误差较小,对类似的小概率失效事件和强非线性可靠性分析问题具有一定的指导意义,例如本文中航空发动机结构疲劳可靠性分析案例中本方法的计算误差最大不超过0.052 84,和MCS方法计算结果贴合度较高。

2) 本文提出的AUMCROAIS方法可以结合仿真,在完善随机参数信息的情况下,可得到较精确的疲劳可靠度计算结果,同时计算效率大大提高,总算次数约为MCS方法的1/3左右。

参考文献
[1] U.S. Department of Transportation, Federal Aviation Administration. Guidance material for aircraft engine life limit parts requirement: Advisory Circular 33.70-1[R]. Washington, D.C.: FAA, 2009: 30-50.
[2] VITTALS S, HAJELA P, JOSHI A. Review of approaches to gas turbine life management: AIAA-2004-4372[R]. Reston, VA: AIAA, 2004.
[3] BEACHKOFSKI B K, GRANDHI R V. Probabilistic system reliability for a turbine engine airfoil[C]//ASME Turbo Expo 2004: Power for Land, Sea, and Air. New York: ASME, 2004: 171-179.
[4] MELIS M E, ZARETSKY E V, AUGUST R. Probabilistic analysis of aircraft gas turbine disk life and reliability[J]. Journal of Propulsion and Power, 1999, 15(5): 658-666.
Click to display the text
[5] GAO H, FEI C, BAI G. Reliability-based low-cycle fatigue damage analysis for turbine blade with thermo-structural interaction[J]. Aerospace Science & Technology, 2016, 49: 289-300.
Click to display the text
[6] ZHU S P, LIU Q, PENG W, et al. Computational-experimental approaches for fatigue reliability assessment of turbine bladed disks[J]. International Journal of Mechanical Sciences, 2018, 142: 502-517.
Click to display the text
[7] ZHU S P, LIU Q, ZHOU J. Fatigue reliability assessment of turbine discs under multi-source uncertainties[J]. Fatigue & Fracture of Engineering Materials & Structures, 2018, 4: 1291-1305.
Click to display the text
[8] HU D Y, MA Q H, SHANG L H, et al. Creep-fatigue behavior of turbine disc of superalloy GH720Li at 650℃ and probabilistic creep-fatigue modeling[J]. Materials Science & Engineering A, 2016, 670: 17-25.
Click to display the text
[9] WANG R Q, LIU X, HU D Y, et al. Zone-based reliability analysis on fatigue life of GH720Li turbine disk concerning uncertainty quantification[J]. Aerospace Science and Technology, 2017, 70: 300-309.
Click to display the text
[10] 李岩, 张曙光, 宫綦. 一种改进的航空发动机结构概率风险评估方法[J]. 航空学报, 2016, 37(2): 597-608.
LI Y, ZHANG S G, GONG Q. An improved probabilistic risk assessment method of structural parts for aeroengine[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(2): 597-608. (in Chinese)
Cited By in Cnki (1) | Click to display the text
[11] 陈志英, 谷裕, 周平, 等. 基于应力应变场强法的轮盘疲劳可靠性分析方法[J]. 推进技术, 2018, 39(2): 433-439.
CHEN Z Y, GU Y, ZHOU P, et al. Fatigue reliability analysis of disk based on stress-strain field intensity method[J]. Journal of Propulsion Technology, 2018, 39(2): 433-439. (in Chinese)
Cited By in Cnki | Click to display the text
[12] 刘君强, 谢吉伟, 左洪福, 等. 基于随机Wiener过程的航空发动机剩余寿命预测[J]. 航空学报, 2015, 36(2): 564-574.
LIU J Q, XIE J W, ZUO H F, et al. Residual lifetime prediction for aeroengines based on Wiener process with random effects[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(2): 564-574. (in Chinese)
Cited By in Cnki (7) | Click to display the text
[13] 姚伟, 白广忱. 基于Fourier正交基神经网络的涡轮盘低循环疲劳可靠性分析[J]. 装备制造技术, 2014(10): 132-134.
YAO W, BAI G C. Reliability analysis of low cycle fatigue of turbine disk based on fourier orthogonal neural network[J]. Equipment Manufacturing Technology, 2014(10): 132-134. (in Chinese)
Cited By in Cnki (2) | Click to display the text
[14] 皮骏, 马圣, 贺嘉诚, 等. 遗传算法优化的SVM在航空发动机磨损故障诊断中的应用[J]. 润滑与密封, 2018, 43(10): 95-103.
PI J, MA S, HE J C, et al. Application of genetic algorithm optimized SVM in aeroengine wear fault diagnosis[J]. Lubrication Engineering, 2018, 43(10): 95-103. (in Chinese)
Cited By in Cnki | Click to display the text
[15] TOMASSON E, SODER L. Improved importance sampling for reliability evaluation of composite power systems[J]. IEEE Transactions on Power Systems, 2017, 32(3): 2426-2434.
Click to display the text
[16] LI L, LU Z Z. Interval optimization based line sampling method for fuzzy and random reliability analysis[J]. Journal of Applied Mathematics, 2014, 38(13): 3124-3135.
Click to display the text
[17] XIONG B, TAN H F. A robust and efficient structural reliability method combining radial-based importance sampling and Kriging[J]. Science China Technological Sciences, 2018, 61(5): 724-734.
Click to display the text
[18] MELCHERS R E. Search-based importance sampling[J]. Structural Safety, 1990, 9: 117-128.
Click to display the text
[19] AU S K, BECK J L. A new adaptive importance sampling scheme for reliability calculations[J]. Structural Safety, 1999, 21: 135-158.
Click to display the text
[20] SANKARAN M, PRAKASB L. Adaptive simulation for system reliability analysis of large structures[J]. Computers and Stmctures, 2000, 77: 725-734.
Click to display the text
[21] 黄毅, 王明政, 颜寒, 等. 快堆中心柱低周疲劳可靠性评价[J]. 原子能科学技术, 2018, 52(1): 118-125.
HUANG Y, WANG M Z, YAN H, et al. Low cycle fatigue reliability evaluation of central column in fast reactor[J]. Atomic Energy Science and Technology, 2018, 52(1): 118-125. (in Chinese)
Cited By in Cnki | Click to display the text
[22] 陈向前, 董聪, 闫阳. 自适应重要抽样方法的改进算法[J]. 工程力学, 2012, 29(11): 123-128.
CHEN X Q, DONG C, YAN Y. Improved adaptive important sampling algorithm[J]. Engineering Mechanics, 2012, 29(11): 123-128. (in Chinese)
Cited By in Cnki (4) | Click to display the text
[23] 马纪明, 詹晓燕, 曾声奎. 基于自适应重要抽样的可靠性分析方法[J]. 北京航空航天大学学报, 2011, 37(9): 1142-1150.
MA J M, ZHAN X Y, ZENG S K. Reliability analysis method based on adaptive importance sampling[J]. Journal of Beijing University of Aeronautics and Astronautics, 2011, 37(9): 1142-1150. (in Chinese)
Cited By in Cnki (4) | Click to display the text
[24] 戴鸿哲, 赵威, 王伟. 结构可靠性高效自适应重要抽样方法[J]. 力学学报, 2011, 4(6): 1133-1140.
DAI H Z, ZHAO W, WANG W. An efficient adaptive important sampling method for structural reliability analysis[J]. Chinese Journal of Theoretical and Applied Mechanics, 2011, 4(6): 1133-1140. (in Chinese)
Cited By in Cnki (2) | Click to display the text
[25] U. S. Department of Transportation, Federal Aviation Administration Airworthiness Standards. Aircraft engines: CFR 14 Part 33[S]. Washington, D.C.: FAA, 2009: 10-30.
[26] 吴学仁. 飞机结构金属材料力学性能手册:静强度疲劳-耐久性[M]. 北京: 航空工业出版社, 1997: 30-35.
WU X R. Handbook of mechanical prosperities of aircraft structural metals:Static strength/durability[M]. Beijing: Aviation Industry Press, 1997: 30-35. (in Chinese)
http://dx.doi.org/10.7527/S1000-6893.2019.23228
中国航空学会和北京航空航天大学主办。
0

文章信息

游令非, 张建国, 周霜, 杜小松
YOU Lingfei, ZHANG Jianguo, ZHOU Shuang, DU Xiaosong
航空发动机限寿件疲劳可靠度计算新方法
A new method for fatigue reliability calculation of aero-engine limited life parts
航空学报, 2019, 40(12): 223228.
Acta Aeronautica et Astronautica Sinica, 2019, 40(12): 223228.
http://dx.doi.org/10.7527/S1000-6893.2019.23228

文章历史

收稿日期: 2019-06-19
退修日期: 2019-07-23
录用日期: 2019-08-07
网络出版时间: 2019-08-23 16:55

相关文章

工作空间