在航空发动机适航要求中，将原发失效能够引起发动机危害性影响的部件定义为发动机限寿件(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可以通过不断调整可靠度指标半径(β球)，使其逐渐接近可靠度指标。其通过搜寻随机变量概率密度f_X(·)值最大点而不断修正。设极限状态函数为g(X)，其中X=[x₁  x₂  …  x_n]，MCROAIS估计失效概率${\widehat P_{\rm{f}}}$的计算流程如图 1所示，图中：k为迭代次数；β为可靠度指标；γ为迭代点到均值点的距离；r₁和r₂为半径。

其中更新机制如下：若x_min存在，则β_k+1=β(x_min)；若x_max存在，则比较β(x_max)>γβ_k, γ=1.02，若不等式成立，β_k+1=β(x_max)；若不成立，进入下一次循环。估计失效概率${\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为总抽样次数；n_k为第k次迭代步骤里重要抽样所需样本点数；v_i为MCROIS中的抽样点，i为序列数；χ_n²(·)为自由度为n的卡方分布；I[·]为示性函数。

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

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

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

2.1 极坐标抽样

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

$ {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）

式中：σ²为方差。

位移R²=X²+Y²是随机变量的函数，其分布函数为

$ \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为极径；θ为极角；D是x²+y²≤r²的区域，通过变成极坐标求得结果。对式(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，方差都是σ² (各向同性)时，这时其极坐标极径r和极角θ的分布分别为一个Rayleigh分布和一个均匀分布。考虑位于极坐标中点(μ_r, μ_θ)处的无规行走，即其直角坐标x和y的期望值为μ_X和μ_Y，方差都是σ² (各向同性)。于是

$ {\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）

式中：I₀(x)是零阶第一类Bessel函数，I₀(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.2 主动学习函数的构造

为使得样本点更多落在抽样半径确定的球区域附近，从而加速失效概率计算的收敛，本文构造主动学习函数针对极坐标抽样的结果来识别最优样本点，使得目标函数$\widehat G$(x_i)尽可能接近极限状态函数的阈值a附近，范围为[a－ζ, a+ζ]，ζ为误差。当$\widehat G$(x_i)>a时，表示产品处于失效状态；当$\widehat G$(x_i)＜a时，表示产品处于安全状态；当$\widehat G$(x_i)=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$(x_i)和σ_{$\widehat G$(xi)}²分别为极坐标抽样样本点的极限状态函数或概率密度函数的预测值和方差；Φ(·)为标准正态累积分布函数；φ为标准正态概率密度函数；误差ζ可取max($\widehat G$(x_i))×0.1(i=1, 2, …, N)；x_i为第i个样本点。找出所有样本点主动学习函数值的最大值, 即max(L(x_i))(i=1, 2, …, N)。

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

3 AUMCROAIS可靠度计算方法

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

3.1 自更新机制

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

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

3.2 算法流程

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

步骤2  若为步骤10跳回，重置i=1；按重要抽样密度函数h_V(v)抽样获得n_k个样本点Y_j=[y_1j  y_2j  …  y_nj](j=1, 2, …, n_k)，其中若为第1次迭代(k=1)，σ_hV(v)=σ_Y。

步骤3  对于在失效域的样本点，计算f_X(_kY_j)，记录max(f_X(_kY_j))，将其样本设为Y^*，β^*=║Y^*║₂；对于不在失效域的样本点，计算g(Y_j)，记录min(g(Y_j))，将其样本设为Y^**，β^**=║Y^**║₂。

步骤4  若没有样本点落在失效域内，令β_i=β^**并把h_V(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^*为中心，进行极坐标抽样(p_R(r), p_Θ(θ))，样本点记为x_i(i=1, 2, …, N)。

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

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

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

步骤9  抽取落在抽样区域的随机样本Y′_j。随机样本Y′_j=[y′_{j, 1}  y′_{j, 2}  …  y′_{j, n}]，令a_j= $\frac{{\mathit{\boldsymbol{Y}}{'_j}}}{{\sqrt {\sum {{{(y{'_{i, j}})}^2}} } }}$，则Y′_j=R_ja_j，利用随机样本点Y′_j计算极限状态函数预测值$\widehat G$ (Y′_j)并得到I[$\widehat G$ (Y′_j)]的值：如果$\widehat G$ (Y′_j)＜0，则I[$\widehat G$ (Y′_j)]=1；否则，I[$\widehat G$ (Y′_j)]=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所示。

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。

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

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

压气机轮盘应力比R_L=0的S-N曲线(即压气机轮盘在0~29 189.6 N的脉动载荷下的疲劳强度与疲劳寿命曲线)可由N_life=Kσ_a^－m得到，则

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

式中：N_life为疲劳载荷循环次数；σ_a为循环疲劳应力幅值；S_max为应力场强；m和K为材料参数；压气机轮盘材料为TC11，由《飞机结构金属材料力学性能手册》^[26]得到其材料参数K=6.278 4×10¹⁵，m=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  S_max]；D_f为疲劳退化强度；D_f(t)为疲劳载荷作用下的退化量；n_max为疲劳应力S_max对应的循环载荷次数。考虑单个疲劳应力水平下的模型

$ \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材料参数m和K的均值、变异系数分别为μ_m=4.628、COV_m=0.01，μ_K=6.135 26×10¹⁵，COV_K=0.015，服从正态分布。

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

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

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

1) 仿真效率：通过MCROAIS收敛速度对比，可以验证本方法的仿真效率。如图 11所示，抽取12 000次、24 000次2种情况作为对比其可靠度指标β的收敛过程，可以看出本方法在抽样效率上的优越性，由于结合了极坐标抽样和主动学习函数，使近似验算点能迅速接近真实验算点x^*，不断更新最优抽样半径；另外，在抽样次数上与MCS和MCROAIS进行对比，在每一级循环下，MCS需10⁸次抽样；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点附近。综上所述，本方法在仿真效率上得到了很大的提升。

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左右。