1 引言

航空发动机轮盘的疲劳失效通常发生在结构的榫槽、螺栓孔、销钉孔等应力集中的缺口部位，而轮盘的失效往往会对发动机的运行产生灾难性的影响^[1]，因此对应力集中部位疲劳寿命的准确预测对于轮盘的抗疲劳设计很有意义，针对应力集中部位疲劳寿命预测方法主要有名义应力法、局部应力应变法和场强法^[2]，其中局部应力应变法由于考虑了塑性应变的影响被广泛地应用到轮盘的寿命预测及疲劳可靠性分析中^[3~5]。但局部应力应变法将缺口根部最危险点的应力应变等效为光滑试件上作用的应力应变，这种假设仅在缺口根部严重进入塑性，疲劳区各点应力应变变化不大时才基本成立，而随着对航空发动机推重比要求的不断提高，势必要求发动机的轮盘更加轻质化，以满足其结构效率的要求。研究表明传统局部应力应变法由于其不能考虑缺口区域应力应变分布的影响，对疲劳寿命预测结果偏于保守且稳定性差^[6]，因此传统方法将不能满足航空发动机轮盘疲劳寿命及疲劳可靠性的精度要求。考虑到局部应力应变法局限性，姚卫星^{[7, 8]}提出的场强法认为结构的疲劳寿命除最大应力应变外还与应力峰值点周围的应力应变场有关，并提出场强因子很好地解释了结构的多种疲劳行为，被认为是另一种具有工程应用前景的疲劳寿命分析方法。许多研究已经分别从理论、试验和工程中验证了该方法的精度及稳定性，尚德广^[9]应用场强法对缺口试件的疲劳寿命进行预测，并与局部应力应变法预测结果进行对比发现场强法理论预测结果精度更高；张成成等^[10]则通过系列缺口试件疲劳试验从试验角度验证了理论分析得出的结论；李岩等^[11]应用场强法对某发动机压气机轮盘的疲劳寿命进行预测，并与试验结果吻合良好。在考虑参数不确定性情况下，轮盘的疲劳可靠性研究大多基于点应力应变理论的局部应力应变法^[12~14]，而基于场强法进行疲劳可靠性分析的研究较少。

本文以涡轮盘作为数值算例，利用应力应变场强法对轮盘进行疲劳寿命分析。在考虑各参数随机性的情况下，利用分布式协同响应面模型建立疲劳寿命与随机参数的代理模型，并进行疲劳寿命可靠性分析。

2 应力应变场强法

针对以应变作为控制参量的航空发动机轮盘疲劳寿命分析，定义应变场强因子

$ {\varepsilon _{{\rm{FI}}}} = \frac{1}{V}\int\limits_\mathit{\Omega } {f\left( {{\varepsilon _{ij}}} \right)\varphi \left( r \right){\rm{d}}v} $

(1)

式中ε_FI为应变场强度，Ω为缺口疲劳破坏区，V为Ω的体积，f (ε_ij)为破坏应变函数，φ(r)为权函数。其中Ω只与材料性能有关，为以应力峰值点为中心，数个晶粒大小为半径的球或椭球^[5]；f (ε_ij)对于各向同性的弹塑性金属材料为Von-Mises等效应变公式；φ(r)表示破坏区内r处的应变对峰值应变的贡献(图 1)，其定义如下

$ \varphi \left( {{r_\theta }} \right) = 1 - Cr\left( {1 + \sin \theta } \right) $

(2)

C为相对应变梯度

$ C = \left| {\frac{1}{{{\varepsilon _{\max }}}} \cdot \frac{{{\rm{d}}{\varepsilon _{\left( {r,\theta } \right)}}}}{{{\rm{d}}{r_\theta }}}} \right| $

(3)

3 应变—寿命方程

应变寿命方程用来描述材料应变与寿命的关系，其中总应变寿命方程的应用最为广泛，即Manson-Coffin公式，同时在航空发动机的复杂载荷下，平均应力对疲劳寿命的影响是不可避免的，因此本文选用Morrow平均应力修正公式，其表达式为

$ \frac{{\Delta \varepsilon }}{2} = \frac{{{{\sigma '}_{\rm{f}}} - {\sigma _{\rm{m}}}}}{E}{\left( {2{N_{\rm{f}}}} \right)^b} + {{\varepsilon '}_{\rm{f}}}{\left( {2{N_{\rm{f}}}} \right)^c} $

(4)

式中Δε为总应变，σ_m为平均应力，E为弹性模量，σ′_f为疲劳强度系数，ε′_f为疲劳延性系数，b为疲劳强度指数，c为疲劳延性指数，N_f为循环数。在计算寿命时用Δε_FI作为Δε与σ_m一起带入式(5)即可求出结构的疲劳寿命。

4 疲劳可靠性分析模型

疲劳可靠性通常需要数万次的抽样以得到疲劳寿命的概率响应分布，通过蒙特卡洛抽样会花费大量的时间；而场强法涉及参数众多，直接进行响应面拟合将需要大量的样本点，这势必会因为弹塑性非线性有限元分析造成大量的计算时间，同时场强法由于需要计算场径并进行应变场拟合又大大影响了计算效率；基于以上考虑本文利用分布式协同响应面^[15~17]解决以上问题。

首先建立弹塑性有限元分析部分输出应变场强ε_FI和应力峰值σ_max的响应面函数

$ {\varepsilon _{{\rm{FI}}}}\left( {{\mathit{\boldsymbol{X}}_1}} \right) = {\mathit{\boldsymbol{A}}_1} + {\mathit{\boldsymbol{B}}_1}{\mathit{\boldsymbol{X}}_1} + \mathit{\boldsymbol{X}}_1^{\rm{T}}{\mathit{\boldsymbol{C}}_1}{\mathit{\boldsymbol{X}}_1} $

(5)

$ {\sigma _{\max }}\left( {{\mathit{\boldsymbol{X}}_2}} \right) = {\mathit{\boldsymbol{A}}_2} + {\mathit{\boldsymbol{B}}_2}{\mathit{\boldsymbol{X}}_2} + \mathit{\boldsymbol{X}}_2^{\rm{T}}{\mathit{\boldsymbol{C}}_2}{\mathit{\boldsymbol{X}}_2} $

(6)

式中A_i, B_i, C_i, X_i分别对应各自的常数项，一次项系数矩阵，二次项系数矩阵和随机变量。式(5)，(6)被称为分布式响应面函数。得到ε_FI和σ_max的响应之后，将输出响应作为疲劳寿命分析函数中的输入随机变量，则疲劳寿命分析响应面函数

$ {N_{\rm{f}}}\left( {\mathit{\boldsymbol{\hat X}}} \right) = \mathit{\boldsymbol{\hat A}} + \mathit{\boldsymbol{\hat B\hat X}} + {{\mathit{\boldsymbol{\hat A}}}^{\rm{T}}}\mathit{\boldsymbol{\hat C\hat A}} $

(7)

式中$\boldsymbol{\hat A}, \boldsymbol{\hat B}, \boldsymbol{\hat C}$为响应面函数系数，随机变量$\boldsymbol{\hat X} $包括分布式响应面输出ε_FI和σ_max以及其他随机变量，则式(7)被称为协同响应面函数。

通过以上分布式协同响应面的建立，将疲劳可靠性分析分为两个分布式响应面式(5)，(6)，并利用协同响应面式(7)得到结构的概率寿命响应。通过这样的分解与整合有效地降低了利用场强法进行可靠性分析的模型复杂程度，保证了可靠性模型精度和效率。

5 涡轮盘疲劳可靠性分析

综合以上理论，建立基于应力应变场强法轮盘疲劳寿命分析方法流程(如图 2)，并以某涡轮盘为数值算例验证了方法的有效性。

5.1 场径计算

不同材料的晶粒尺寸和微观结构不同，则不同材料的疲劳损伤区Ω也不同，宏观上认为损伤区是以应力峰值点为中心，数个晶粒大小为半径的球，球的半径称为应力场径。材料的场径一般通过试验确定。本文数值算例中的涡轮盘材料为镍基高温合金GH4169，材料300℃试件高周疲劳应力寿命曲线实验数据及拟合曲线如图 3，应力比0.1。光滑试件应力集中系数为1，缺口试件应力集中系数为3，由疲劳试验数据可知GH4169材料光滑试件的300℃高周疲劳强度为585MPa，缺口试件的300℃高周疲劳强度为235MPa，疲劳缺口减缩系数为2.49。

利用有限元软件对缺口试件进行有限元分析，由于场径大小一般为数个晶粒大小，单个金属晶粒直径大小为0.01mm量级，因此将缺口区域的网格加密，节点距离为0.02mm左右。左端施加固定约束，右端施加235MPa的拉力，温度为300℃，缺口试件的尺寸及有限元等效应力计算结果如图 4。

由结果可知试件最大应力为714.89MPa，应力集中系数714.89/235=3.04，模拟实验基本符合实际。同时由上图可知缺口区域的场强方向为箭头r方向，利用Ansys后处理中的路径功能，得到场强方向上的应力分布如图 5。

因为场径大小一般为数个晶粒大小，数个晶粒的尺度不会超过0.1mm，该尺度下应变随场径变化如图 6，可见在0.1mm范围内的应力应变沿场径变化可以近似为直线^[18]。

设直线方程为

$ f\left( {{\varepsilon _{ij}}} \right) = {\varepsilon _{\left( {r,\theta } \right)}} = {\varepsilon _{\max }} - \frac{{{\varepsilon _{\max }} - {\varepsilon _R}}}{R}r $

(8)

式中ε_max应变峰值，R为场径大小，ε_R为场径为R处的应变。结合式(2) ~ (4)得

$ \varphi \left( {{r_\theta }} \right) = 1 - \frac{{r\left( {1 + \sin \theta } \right)\left( {{\varepsilon _{\max }} - {\varepsilon _R}} \right)}}{{{\varepsilon _{\max }} \cdot R}} $

(9)

式(9)带入式(1)进行积分忽略高阶小量后，得简化表达式

$ {\varepsilon _{{\rm{FI}}}} = {\varepsilon _{\max }} - \frac{{24 + 3{\rm{ \mathsf{ π} }}}}{{32{\rm{ \mathsf{ π} }}}}\left( {{\varepsilon _{\max }} - {\varepsilon _R}} \right) $

(10)

同理地到应力场强表达式

$ {\sigma _{{\rm{FI}}}} = {\sigma _{\max }} - \frac{{24 + 3{\rm{ \mathsf{ π} }}}}{{32{\rm{ \mathsf{ π} }}}}\left( {{\sigma _{\max }} - {\sigma _R}} \right) $

(11)

根据式(11)得到不同场径下的应力场强，通过多项式拟合得到场径与应力场强的关系，如图 7。

根据场强法理论，若缺口根部的应力应变场强度的历程与光滑试件的应力应变场强度的历程相同，二者具有相同的寿命。因此当试件的应力场强等于光滑试件的应力场强585MPa时，即为材料的场径大小，通过计算材料的场径为81.3μm。

5.2 确定性涡轮盘有限元分析

本文选取某销钉连接涡轮转子作为算例，涡轮轴安装凸肩处开有销钉孔，通过径向销钉与盘轴过盈配合实现传扭和定心。选取1/16循环对称结构进行热弹塑性及接触非线性有限元分析^{[19, 20]}，并在接触区域进行网格加密以满足应变场强的计算。根据实际载荷谱^[21]选取该涡轮盘的主循环进行计算，最大转速为11425r/min，为0→11425r/min→0即R=0的脉动循环。结构几何模型及有限元模型见图 8。

考虑叶片、轮盘的离心力的同时考虑盘—销接触应力，在接触区域施加0.01mm的过盈量，采用经验公式模拟轮盘的温度场进行热固耦合计算

$ t = {t_0} + \frac{{\left( {{t_{\rm{a}}} - {t_0}} \right)\left( {{R^m} - R_0^m} \right)}}{{\left( {R_{\rm{a}}^m - R_{\rm{0}}^m} \right)}} $

(14)

式中t₀为盘心温度，t_a为盘缘温度，对于镍基高温合金m=4，其中轮缘温度取650℃，轮心和销钉连接处温度取300℃。强化准则采用随动强化准则，屈服准则采用Von-Mises屈服准则。最大转速下轮盘热弹塑性分析结果如图 9。

根据计算结果，涡轮盘最大应力点出现在涡轮轴安装凸肩的销钉孔P处，最大应变ε_max为0.51239%，根据上一部分计算的场径81.3μm，提取P为中心，图中r方向为场强方向上节点的应变，通过多项式拟合如图 10所示。

可见，轮盘的疲劳破坏区将内应变随场强变化仍近似为直线，图中r=81.3μm对应的应变值即为ε_R，将ε_R带入公式(6)，计算得P处的应变场强ε_FI为0.50609%，利用通用斜率法得到300℃的材料疲劳参数，对于脉动循环平均应力σ_m = σ_max/2，将σ_m和σ_FI带入式(5)得到涡轮盘疲劳寿命为13115.22循环数。

5.3 选取随机变量

随机变量的选取分为两部分，弹性塑性分析部分选取弹性模量E，销钉孔处温度T，转速ω和销钉装配过盈量f作为随机变量，疲劳可靠性分析部分将最大应力σ_max和应变场强ε_FI考虑成随机变量的同时选取材料疲劳参数作为随机变量进行疲劳可靠性分析。选取的随机变量及其分布特征见表 1。

5.4 建立分布式协同响应面

首先建立最大应力σ_max和应变场强ε_FI关于随机变量X =[E, T, ω, f]^T的响应面。基于Box-Behnken矩阵设计方法，由42组试验数据拟合得到σ_max和ε_FI关于随机变量X =[E, T, ω, f]^T的响应面函数，通过拉丁超立方抽样方法进行1.0×10⁵次抽样，得到σ_max和ε_FI的概率响应，如图 11，12所示。由抽样结果拟合，σ_max和ε_FI服从正态分布，其均值分别为1033.9MPa，0.50868%，标准差分别为97.634MPa，0.049149%。

将σ_max和ε_FI的输出响应作为输入随机变量，同时考虑$ \boldsymbol{\hat X} = {[E, \sigma {'_{\rm{f}}}, \varepsilon {'_{\rm{f}}}, b, c, {\varepsilon _{{\rm{FI}}}}, {\sigma _{\max }}]^{\rm{T}}}$对疲劳寿命Nf的影响，同样基于Box-Behnken矩阵试验设计方法，得到协同响应面函数，通过1.0×10⁴次抽样得到疲劳寿命N_f的概率响应，如图 13所示。通过拟合，寿命响应N_f服从对数正态分布。

5.5 灵敏度分析

对可靠性分析中的随机变量进行灵敏度分析，以获得疲劳寿命对每个随机变量的敏感程度，如图 14所示。其中疲劳延性指数c是影响疲劳寿命分散性的最主要因素，弹性模量E对疲劳寿命影响最小。根据灵敏度分析结果，可以为轮盘的疲劳可靠性设计提供方向。

5.6 疲劳可靠性分析

利用响应面模型得到疲劳寿命累计概率密度函数，如图 15所示，计算不同可靠度下的轮盘疲劳寿命，见表 2。

由表 2可知，局部应力应变法由于没有考虑应力应变场的分布，理论分析的确定性寿命与概率寿命与场强法相比偏于保守，结合文献[10]中试验数据，理论上本文方法分析结果精度更高。而可靠性分析与确定性分析对比可以发现，若涡轮盘的疲劳可靠性要求是疲劳寿命为6000循环数的可靠度不低于99.9%，确定性分析考虑2倍的安全系数后的安全疲劳寿命为6557循环。可靠性分析结果为9595循环，可见确定性分析确定的安全寿命相比可靠性分析确定的寿命更加保守。

6 结论

通过研究得到结论如下：

(1) 基于应力应变场强法对轮盘进行疲劳可靠性分析，99.90%可靠度下的疲劳寿命为9595循环，考虑2倍安全系数的安全寿命为6557循环，可靠性分析结果表明轮盘具有一定的寿命裕度，可对轮盘进一步优化设计，以提高结构效率。

(2) 利用分布式协同响应面方法建立涉及多变量、多学科的疲劳可靠性分析模型，可以解决可靠性分析中的效率与精度问题。

(3) 通过对比理论分析结果，当疲劳破坏区具有较大应力应变梯度时，传统方法结果趋于保守，而本文方法由于考虑缺口应力应变分布对疲劳寿命的影响，理论上能够提高疲劳可靠性分析结果的准确性。