作为航空发动机中极其重要的部件, 滚动轴承的寿命问题已经成为提高我国航空发动机推重比的障碍之一.轴承游隙是轴承设计与使用中的焦点问题, 对轴承寿命的影响十分显著.近年来, 有学者对轴承游隙做出了一些研究^[1], 并逐渐考虑了轴承结构尺寸、材料性能以及工况的随机性.随着计算机辅助工程(computer aided engineering, CAE)技术与大型商业有限元仿真软件的不断发展, 越来越多的工程结构可靠性问题的极限状态函数都是建立在有限元仿真的基础之上, 即是隐式的.目前较为流行的解决隐式结构可靠性问题的手段是通过构建一个简单的等价模型来替代原来的复杂极限状态函数, 进而对等价模型进行可靠性分析^[2].多项式响应面法^[3]是研究最早的一种代理模型近似技术, 实际上, 由于许多工程力学模型并不是都具有很强的非线性, 而响应面法由于其简单有效, 目前在可靠性分析中仍被广泛地使用.自从响应面法被Faravelli首次引入到结构可靠性理论后, 就有许多学者对其做出了大量的改进.而在众多的改进方案中脱颖而出的是具有较高拟合精度和计算效率的自适应响应面法^[4].在这之后, 研究人员发现确定多项式响应面的试验点对结果有很大影响, DOE选取不当会使结果完全错误.为了使试验点能够距离设计点更近, Kim等将星形DOE向极限状态面上映射, 提出了梯度投影法^[5].Das等则通过累积使用已有的试验点来改进自适应响应面法^[6].类似地, Nguyen根据随机变量对极限状态函数的灵敏度以及初始试验点的规模在半星形的DOE内不断增加样本点^[7].本文在建立响应面时, 根据误差预测标准以及交叉验证方法对多项式响应面中最重要的项进行筛选, 并提出一种在多维情况下更合理的实验设计方案.通过这两方面对经典的自适应响应面法做出改进, 使其更加适合多维结构可靠性分析问题.最后结合有限元法解决了滚动轴承游隙的可靠性分析问题.

1 稀疏响应面

假设基于有限元的滚动轴承工作游隙极限状态函数已经确定为y=g(x), 式中y为工作游隙极限状态函数的真实响应, x为影响工作游隙的随机变量.若利用多项式响应面来近似g(x), 则有

(1)

式中: 为响应面函数; ε为误差. 的表达式一般为

(2)

式中:z_i可以为线性项、平方项或交叉项; a_i为各项系数; P为项数.回归方程的矩阵形式可表示为

(3)

式中:y=(y₁, y₂，…，y_N)^T为学习样本点的真实响应向量; Z为N×P回归因子矩阵; a=(a₁, a₂，…，a_P)^T为各项的系数矩阵; ε=(ε₁, ε₂，…，ε_P)^T为误差矩阵.根据回归理论, 利用最小二乘法可估计出系数矩阵, 即

(4)

由于滚动轴承游隙可靠性问题的随机变量多达14个, 若采用全模型响应面, 将导致过拟合现象^[8].这一缺陷能够通过减少响应面中的项来改善.

较为简单的衡量响应面模型质量的方式通常有经验均方误差ε_emp或决定系数R², 其表达式为

(5)

(6)

式中y为经验均值.MSE的表达式为

(7)

本文在估计预测误差期望时采用两类方法:

1) 罚函数法.C_p统计量常被用来估计MSE.若根据N个样本点拟合出具有ξ项的稀疏响应面, 且误差平方和, 则C_p的值为

(8)

式中，为全模型响应面中σ²的无偏估计量.C_p统计量能够反映出项数的增加对预测误差的作用.具有最小C_p值的响应面可视为最优模型.其他两种罚函数法分别为AIC准则和BIC准则, 计算公式分别为

(9)

(10)

2) 交叉验证法.K折交叉验证将N个学习样本点分割成K个子样本, 一个单独的子样本被保留作为验证模型的数据, 而其他K-1个子样本用来训练得到响应面.平均K次误差结果得到一个单一估测值ε_CV:

(11)

根据式(11)衍生出预测系数Q²:

(12)

罚函数法计算方便, 计算量小.交叉验证法更准确但计算量较大.在对响应面中的项进行选择时应充分利用这两类方法的优势.

2 基于极变换的实验设计

实验设计对响应面的拟合有较大影响, 如果每一次构建响应面时能够增加对临界样本点的拟合, 则最终的响应面能够对真实的极限状态函数做很好的近似.

2.1 极变换

极变换由Hurtado首次提出^[9].假定在标准正态随机变量空间u中, 首先根据FORM(一次可靠性方法)法得到其设计点为u^*.若空间中的一个样本向量为u_i=(u₁, u₂，…，u_d), d为随机变量空间的维度, 设计点单位向量w则被定义为

(13)

定义样本的极特征分别为样本向量的长度r和样本向量与设计点单位向量之间夹角φ的余弦, 即

(14)

(15)

文献[9]已经讨论了在小概率情况下, 从标准正态随机变量空间映射到平面内的失效类数据点具有显著的聚集特性.

2.2 基于极变换的实验设计方案

本文利用极限状态方程的二次近似在极特征平面内的映射来对临界样本点进行选择.根据极特征的定义有

(16)

且

(17)

将式(16)与式(17)代入g(u)=0, 得到极限状态方程的二次近似在极特征平面内的映射, 即

(18)

若给式(18)中的λ赋予不同的值, 则其对应的不同曲线位置如图 1所示, 可以看出, λ值越大, 曲线的位置越靠上.对于极特征平面内的每一个数据点v, 其对应的λ值可由式(18)推导而来, 即

(19)

本文通过调整两条具有不同λ值的曲线的位置, 将与两条曲线中间数据点对应的标准正态空间中的点作为最接近临界状态的样本点, 具体步骤如下.

1) 将由响应面法拟合出的显式极限状态函数g(x)转换到标准正态随机变量空间, 产生N个Monte Carlo样本点并将其映射到极特征平面内;

2) 对于极特征平面内的每一个数据点, 计算与其对应的极限状态函数g(x)的响应值, 并根据式(19)计算与其对应的λ值;

3) 确定上端曲线的λ值λ_max.首先寻找与具有最小|λ|值的曲线对应的λ值λ₀, 令λ_max=λ₀.若λ值为λ_max的曲线上方有安全类数据点, 则将λ_max更改为所有比λ_max大的λ中的最小值, 直到曲线上方全部为失效类数据点为止, 即确定了λ_max;

4) 确定下端曲线的λ值λ_min.首先令λ_min=λ₀.若λ值为λ_min的曲线下方有失效类数据点, 则将λ_min更改为所有比λ_min小的λ中的最大值, 直到曲线下方全部为安全类数据点为止, 即确定了λ_min;

5) 将与λ值为λ_max和λ_min的两条曲线间的数据点对应的样本点作为下次响应面拟合的候选实验点, 如图 2所示.

3 滚动轴承游隙可靠性分析

为了提高滚动轴承工作游隙极限状态函数的拟合精度并进行可靠性分析, 本文在稀疏响应面的基础上, 结合基于极变换的抽样方法, 提出一种混合可靠性算法, 其计算流程如图 3所示.

为了避免传动轴发生抱死事故, 本文将工作游隙的安全阈值定为0 μm, 则工作游隙的极限状态函数为g(x)=μ_w(x)-0, 其中μ_w(·)为工作游隙函数.可以看出, 建立工作游隙的极限状态函数实际是要确定工作游隙的计算函数, 而工作游隙计算方式的不同则对应着极限状态函数的不同形式.本文在工作游隙的可靠性模型中同时考虑滚动轴承的尺寸、材料属性以及过盈配合、离心力和温升等随机因素的影响.则滚动轴承工作游隙即为原始游隙与过盈配合、离心力和温升引起的轴承游隙变化的差值, 即

(20)

式中:μ_w为工作游隙; μ_o为原始游隙; Δμ_p, Δμ_n和Δμ_T分别为过盈配合、离心力和温升引起的轴承游隙变化量.假设各随机变量均服从正态分布且相互独立.若借助商业有限元软件求解游隙的变化量, 其极限状态函数则是隐式的.轴承内圈与轴配合的有限元模型中, 轴承内圈的内表面是目标面, 而轴的外表面是接触面.在轴承外圈与轴承座配合中, 轴承外圈的外表面是目标面, 而轴承座的内表面是接触面.由于皆是轴对称模型, 因此只建立四分之一的实体以减少运算量.实体单元采用SOLID186单元.在均值点处仿真得到轴承内圈外径位移量Δμ_i=15.933 μm, 轴承外圈内径位移量Δμ_o=8.018 μm, 如图 4所示.

因此由于过盈配合引起的轴承游隙减少量应为Δμ_p=Δμ_i+Δμ_o=23.951 μm.轴承内圈以n=15 ×10³r/min旋转时, 因离心力而产生的位移量如图 5所示, 即轴承径向游隙的减少量为Δμ_n=12.976 μm.

轴承运转达到热平衡后, 根据轴承热变形仿真分析得到的结果如图 6所示.

从轴承内外圈的径向热变形图可知, 内圈外表面的热变形最大, Δu_Ti=29.727 μm, 外圈内表面的热变形最大, Δu_To=11.257 μm, 轴承内圈的膨胀量大于轴承外圈的膨胀量.滚动体热变形2Δu_Tb=16.844 μm.因此在本例中温升使轴承游隙减小.隐式极限状态函数计算得到的均值处的游隙变化量Δμ=Δμ_p+Δμ_n+Δμ_T =23.951+12.976+35.314=72.241 μm.则可得g(x)=100-72.241=27.759 μm.对于不同的随机变量向量x, 根据以上步骤均能返回极限状态值g(x), 即建立了滚动轴承工作游隙的隐式极限状态函数.

通过以上研究, 调用335次有限元程序得到滚动轴承工作游隙的极限状态函数为g(x)=μ_o-0.018D₁²+0.007D₁α+0.004D₁t₂+6.269D₁+0.032D_b²+0.001D_bα-0.023D_bt₃-1.616D_b-0.652δ_i-0.702δ_o+0.022E²+0.412Ed₂+4.763Ek+2.863Eρ-46.109E-0.015αd₂-0.122αt₁+0.170αt₂-0.048αt₃+0.622α-2.932d₂²-0.497d₂k+0.001d₂n-0.121d₂ρ-0.003d₂t₁+700.035d₂-6.272k²-0.003kn+1.095kρ+108.809k-0.001nρ-0.167n+0.018ρ²+14.592ρ+0.408t₁-0.699t₂+0.559t₃-42 238.896.

得到的设计点x^*为[0.042, -0.014, 0.045, -1.351×10^-8, -0.001, -0.029, 0.928, -0.842 4, 1.064, 1.870, -1.307, 0.740, 0.009, 0.333], 可靠度指标β为2.927 7.经过本文的稀疏响应面拟合, 再经过FORM(一次可靠性方法)以及IS(重要抽样法)计算的结果与直接Monte Carlo法进行对比, 结果如表 1所示.

由于基于有限元的极限状态函数涉及的随机变量个数较多, 而应用直接优化方法求解设计点的过程难以收敛, 因此不能用直接FORM法求解其可靠度, 本文方法避免了求解隐式函数的梯度, 并具有较高的精度.

4 结论

1) 应用稀疏响应面对滚动轴承工作游隙极限状态函数进行拟合的策略更加合理, 避免了用全响应面模型拟合由于随机变量个数较多而带来的过拟合现象.

2) 利用极限状态方程的二次近似在极特征平面内的映射来对临界样本点进行选择, 利用了极变换后安全与失效类样本点在平面内的聚集特性, 在多维情况下能更方便直观地选择处于临界状态的样本点, 提高了响应面的拟合效率.