风电调心滚子轴承的接触载荷与寿命可靠性分析

Analysis of contact load and lifespan reliability of wind power self-aligning roller bearing

摘要

以风力发电机主轴轴承为研究对象，建立了调心滚子轴承的静力学模型，评估了轴承的内部接触特性和疲劳寿命。同时，采用Romax模型验证了本文模型的准确性。此外，基于自适应Kriging结合Monte Carlo模拟法（AK-MCS）和L-P疲劳寿命理论提出了调心滚子轴承可靠性分析方法。最后，对轴承的寿命可靠性进行了定量评估。结果表明：在本文坐标系下，正向轴向力和调心角的增大使得调心滚子轴承两列滚子的偏载现象加强，寿命可靠度降低；随着滚子轮廓半径、滚子直径和初始接触角的增大，轴承的寿命可靠度升高；而随着内、外圈滚道轮廓半径和弹性模量的增大，轴承的寿命可靠度降低，其中滚子轮廓半径和内圈滚道轮廓半径对轴承的寿命可靠度影响最大。

Abstract

Taking the bearing for the main shaft of a wind turbine as the research object， this paper establishes a statics model of the self-aligning roller bearing， which is subsequently employed to evaluate the internal contact characteristics and fatigue life of the bearing. Furthermore， the accuracy of the proposed model is validated using the Romax model. Additionally， a reliability analysis methodology for the self-aligning roller bearing is proposed， grounded in the Adaptive Kriging Monte Carlo Simulation （AK-MCS） and the L-P fatigue life theory. Finally， a quantitative assessment of the life reliability of the bearing is conducted. The results indicate that， within the coordinate system adopted in this study， the increase in positive axial force and self-aligning angle intensifies the uneven load distribution among the two rows of rollers in the self-aligning roller bearing， leading to a decrease in life reliability. Conversely， an enlargement in the roller contour radius， roller diameter， and initial contact angle enhances the life reliability of the bearing. However， an increase in the contour radii of the inner and outer raceways and elastic modulus results in a reduction in the life reliability. Notably， the roller contour radius and the inner raceway contour radius exhibit the most significant influence on the life reliability of the bearing.

Graphical abstract

关键词

机械设计及理论 / 调心滚子轴承 / 接触载荷 / 自适应Kriging / 轴承寿命可靠性

Key words

mechanical design and theory / self-aligning roller bearing / contact load / adaptive Kriging / reliability of bearing life

0 引 言

近年来，我国持续推进风电产业的大力发展，风机装载容量呈现逐年增加的趋势^［1］。然而，风力发电机一般安装在环境恶劣的地区，且主机通常位于数十米以上的高空中，使得风机的运行保养和故障维修非常不便。作为风机主轴系统的核心部件之一，主轴调心滚子轴承在各种极端工况下的可靠性严重影响风机运行的稳定性。因此，急需开发一种风力发电机主轴调心滚子轴承可靠性分析方法。

Hertz接触理论研究了物体接触力与变形的关系，为轴承载荷分析提供了重要依据^［2］。Ohta等^［3］在赫兹理论的条件下，给出了球轴承在组合载荷下的载荷-位移关系式。Yang等^［4］通过改变轴承的转速、轴向和径向载荷探究了轴承位移、载荷及刚度的变化规律。Xu等^［5］提出了一种双列角接触球轴承的五自由度模型，分析得到了内圈错位将对轴承接触载荷与刚度产生显著影响的结论。Wu等^［6］提出了一种采用有限元分析预紧载荷对圆锥滚子轴承接触应力影响的方法，分析了不同预紧力下的接触应力分布和变形情况，并与实验结果进行了比较。

对于轴承的疲劳寿命，许多学者也进行了大量研究。Xu等^［7］根据Gumbel分布算法、雨流计数算法、线性累积损伤定律和S-N曲线理论得到了轴承的等效损伤载荷，采用有限元法计算了其疲劳寿命。Bogdan等^［8］分析了轴承内圈错位对圆柱滚子轴承疲劳寿命的影响。Mostafa等^［9］讨论了压力、温度与运行时间等因素对轴承滚动接触疲劳的影响。Warda等^［10］建立了一种考虑轴承几何参数、综合载荷和套圈不对中等因素的圆柱滚子轴承寿命预测模型。Shi等^［11］将滚动轴承准动力学、微观热分析、非高斯曲面模拟技术和有限元应力分析相结合，提出了一种高速重载球轴承相对疲劳寿命的计算方法。

然而，多数研究认为轴承参数是不变的，实际上由于加工和装配等因素，其结构材料参数存在随机不确定性。这种特性使得接触载荷模型产生一定误差，进而影响轴承的疲劳寿命可靠性分析。本文基于Hertz接触理论，建立了调心滚子轴承三自由度静力学模型。以某型号风力发电机主轴轴承为研究对象，分析计算了轴承的接触载荷和疲劳寿命。采用Romax软件进行仿真验证。然后，考虑部分结构参数为随机变量，运用AK-MCS方法进行调心滚子轴承的疲劳寿命可靠性分析。最后，探究了轴承的轴向力与调心角对其接触载荷分布及寿命可靠度的影响。本文提出的方法可以为提高风电调心滚子轴承的寿命可靠性提供一定的理论参考依据。

1 调心滚子轴承接触载荷分析

调心滚子轴承的双列滚子形式可以承受较强的径向载荷和一定的轴向载荷。由于内外圈之间没有约束关系，调心滚子轴承不能承受力矩作用^［12］。

为了分析轴承在外部载荷下的滚子-滚道接触力及相应位移，分别以外圈中心和每个滚子中心为原点，建立轴承全局和局部坐标系，如图1所示。 F =［F_x，F_y，F_z ］为轴承内圈受到的外载荷， δ =［δ_x，δ_y，δ_z ］为内圈在各载荷作用下的对应位移。

调心滚子轴承的滚子与内外圈以球面滚道的形式接触，当内外圈存在不同轴心的角度误差时，两者的径向截面之间将产生一定的倾斜角，实现自动调心^［13］。本文假设轴承外圈固定，当内圈发生图示顺时针方向的偏斜时，调心角α_s为正，反之为负。

如图2所示，调心滚子轴承滚子的排列方式有并排和错排分布两种类型，即第1列滚子与第2列滚子之间可能存在交错角θ_1t2。

如图2规定滚子编号j_1/2，则第1、2列第j个滚子的位置角为：

{\theta }_{\mathrm{1}j/\mathrm{2}j}=\frac{\mathrm{2}\mathrm{\pi }}{Z}\left(j-\mathrm{1}\right)+{\theta }_{\mathrm{0}}(+{\theta }_{\mathrm{1}t\mathrm{2}})（1）

式中：交错角θ_1t2仅对第2列滚子的位置角有效。

根据几何关系，可得到内圈滚道第j个轴向截面沿Y_r1j/r2j 轴的位移为：

\left\{\begin{array}{l}{\delta }_{\mathrm{i}\mathrm{1}j}={\delta }_y\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{1}j}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}j}+{\delta }_z\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{1}j}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{1}j}-{\delta }_x\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{1}j}\\ {\delta }_{\mathrm{i}\mathrm{2}j}={\delta }_y\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{2}j}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{2}j}+{\delta }_z\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{2}j}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{2}j}+{\delta }_x\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{2}j}\end{array}\right.（2）

式中：α_1j 和α_2j 为第1、2列第j个滚子的接触角。

当轴承产生调心角α_s时，各滚子的接触角α_1j/2j 将随位置角θ_1j/2j 的不同而发生变化，具体为：

（3）

式中：上面和下面符号分别用于第1、2列滚子；{\alpha }_{\mathrm{i}\mathrm{n}\mathrm{i}}为轴承的初始接触角。

当轴承内圈在全局坐标系中产生位移时，滚子也会发生相应的运动。因此，第1、2列第j个滚子与内外圈滚道之间的相对位移分别为：

\left\{\begin{array}{l}{\delta }_{\mathrm{i}\mathrm{r}\mathrm{1}j/\mathrm{i}\mathrm{r}\mathrm{2}j}={{\delta }_{\mathrm{i}\mathrm{1}j/}}_{\mathrm{i}\mathrm{2}j}-{\delta }_{r\mathrm{1}j/r\mathrm{2}j}\\ {\delta }_{\mathrm{o}\mathrm{r}\mathrm{1}j/\mathrm{o}\mathrm{r}\mathrm{2}j}={\delta }_{\mathrm{r}\mathrm{1}j/\mathrm{r}\mathrm{2}j}\end{array}\right.（4）

式中：δ_r1j 、δ_r2j 分别为第1、2列第j个滚子沿Y_r1j 、Y_r2j 轴的位移，下角标i和o分别表示内外圈。

进一步可以得到第1、2列第j个滚子与内外圈滚道之间的法向接触变形分别为：

{\delta }_{\mathrm{i}\mathrm{1}jp/\mathrm{o}\mathrm{1}jp/\mathrm{i}\mathrm{2}jp/\mathrm{o}\mathrm{2}jp}={\delta }_{\mathrm{i}\mathrm{r}\mathrm{1}j/\mathrm{o}\mathrm{r}\mathrm{1}j/\mathrm{i}\mathrm{r}\mathrm{2}j/\mathrm{o}\mathrm{r}\mathrm{2}j}-\frac{P_{\mathrm{d}}}{\mathrm{4}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{1}j/\mathrm{2}j}}（5）

式中：P_d为调心滚子轴承的径向游隙。

如图3所示，轴承滚子轮廓半径略小于内外圈滚道轮廓半径，两者将形成接触椭圆，2a和2b分别为椭圆的长半轴和短半轴，L_w为滚子有效长度。Palmgen和Lundberg在滚子与滚道的线接触判断条件中指出，当2a<L_w时，属于点接触；当2a>L_w时，属于线接触。在中等载荷条件下，调心滚子轴承滚子⁃滚道接触类型属于Hertz点接触^［14］。

根据Hertz点接触理论，轴承第1、2列第j个滚子与内外圈滚道之间的法向接触载荷分别为：

\left\{\begin{array}{l}{Q}_{\mathrm{i}\mathrm{1}j/\mathrm{i}\mathrm{2}j}=K_{\mathrm{i}}{\delta }_{\mathrm{i}\mathrm{1}jp/\mathrm{i}\mathrm{2}jp}^{\frac{\mathrm{3}}{\mathrm{2}}}\\ Q_{\mathrm{o}\mathrm{1}j/\mathrm{o}\mathrm{2}j}=K_{\mathrm{o}}{\delta }_{\mathrm{o}\mathrm{1}jp/\mathrm{o}\mathrm{2}jp}^{\frac{\mathrm{3}}{\mathrm{2}}}\end{array}\right.（6）

式中：K_i、K_o分别为内外圈滚道和滚子的接触刚度系数。若δ_i1jp/i2jp 、δ_o1jp/o2jp 的计算值为负数，表明滚子未受载，应令相应的Q_i1j/i2j 、Q_o1j/o2j 等于零。

其中，K_i和K_o可表达为^［15］：

K_{\mathrm{i}/\mathrm{o}}=\frac{\mathrm{2.0678}\mathrm{\pi }}{\mathrm{3}\left(\frac{\mathrm{1}-{{\nu }_{\mathrm{i}/\mathrm{o}}}^{\mathrm{2}}}{E_{\mathrm{i}/\mathrm{o}}}+\frac{\mathrm{1}-{{\nu }_{\mathrm{r}}}^{\mathrm{2}}}{E_{\mathrm{r}}}\right)}\times

\sqrt[]{\frac{\left[\mathrm{2.000}\mathrm{6}+\mathrm{1.193}\mathrm{6}\left(\frac{R_{\mathrm{i}\mathrm{r}x/\mathrm{o}\mathrm{r}x}}{R_{\mathrm{i}\mathrm{r}y/\mathrm{o}\mathrm{r}y}}\right)\right]{\left(\frac{R_{\mathrm{i}\mathrm{r}y/\mathrm{o}\mathrm{r}y}}{R_{\mathrm{i}\mathrm{r}x/\mathrm{o}\mathrm{r}x}}\right)}^{\mathrm{1.272}}}{{\left[\mathrm{1.527}\mathrm{7}+\mathrm{0.602}\mathrm{3}\mathrm{l}\mathrm{n}\left(\frac{R_{\mathrm{i}\mathrm{r}y/\mathrm{o}\mathrm{r}y}}{R_{\mathrm{i}\mathrm{r}x/\mathrm{o}\mathrm{r}x}}\right)\right]}^{\mathrm{3}}\sum {\rho }_{\mathrm{i}\mathrm{r}/\mathrm{o}\mathrm{r}}}}（7）

式中：ν_i/o、ν_r和E_i/o、E_r分别为内外圈和滚子的泊松比及弹性模量；R_irx/orx 和R_iry/ory 分别为内外圈与滚子的接触椭圆在长半轴和短半轴方向上的等效半径；∑ρ_ir/or分别为内外圈与滚子接触处的主曲率和。

内外圈滚道和滚子接触处的等效半径分别为：

\left\{\begin{array}{l}\frac{\mathrm{1}}{R_{\mathrm{i}\mathrm{r}x/\mathrm{o}\mathrm{r}x}}=\frac{\pm \mathrm{2}\gamma }{D(\mathrm{1}\mp \gamma )}+\frac{\mathrm{2}}{D}\\ \frac{\mathrm{1}}{R_{\mathrm{i}\mathrm{r}y/\mathrm{o}\mathrm{r}y}}=-\frac{\mathrm{1}}{r_{\mathrm{i}/\mathrm{o}}}+\frac{\mathrm{1}}{R}\end{array}\right.（8）

式中：上面和下面符号分别用于内、外圈；D为滚子直径；R、r_i和r_o分别为滚子、内圈和外圈滚道在轴向截面内的轮廓半径。

\gamma =\frac{D\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{i}\mathrm{n}\mathrm{i}}}{d_{\mathrm{m}}}（9）

式中：d_m为节圆直径。

\sum {\rho }_{\mathrm{i}\mathrm{r}/\mathrm{o}\mathrm{r}}=\frac{\pm \mathrm{2}\gamma }{D(\mathrm{1}\mp \gamma )}-\frac{\mathrm{1}}{r_{\mathrm{i}/\mathrm{o}}}+\frac{\mathrm{2}}{D}+\frac{\mathrm{1}}{R}（10）

式中：上面和下面符号分别用于内、外圈。

本文调心滚子轴承用于风力发电机主轴，内圈转速低，圆周运动产生的离心力可忽略不计，即内外圈滚道与滚子之间的法向接触载荷近似相等，可列得第1、2列第j个滚子的平衡方程分别为：

\left\{\begin{array}{l}{Q}_{\mathrm{i}\mathrm{1}j}-Q_{\mathrm{o}\mathrm{1}j}=\mathrm{0}\\ Q_{\mathrm{i}\mathrm{2}j}-Q_{\mathrm{o}\mathrm{2}j}=\mathrm{0}\end{array}\right.（11）

同时，滚子对内圈的法向接触载荷应与外载荷F_x 、F_y 和F_z 平衡，可列得轴承整体平衡方程为：

\left\{\begin{array}{l}{F}_x+{\displaystyle \sum _{j=\mathrm{1}}^Z}{Q}_{\mathrm{i}\mathrm{1}j}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{1}j}-{\displaystyle \sum _{j=\mathrm{1}}^Z}{Q}_{\mathrm{i}\mathrm{2}j}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{2}j}=\mathrm{0}\\ F_y-{\displaystyle \sum _{j=\mathrm{1}}^Z}{Q}_{\mathrm{i}\mathrm{1}j}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{1}j}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}j}-{\displaystyle \sum _{j=\mathrm{1}}^Z}{Q}_{\mathrm{i}\mathrm{2}j}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{2}j}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{2}j}=\mathrm{0}\\ F_z-{\displaystyle \sum _{j=\mathrm{1}}^Z}{Q}_{\mathrm{i}\mathrm{1}j}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{1}j}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{1}j}-{\displaystyle \sum _{j=\mathrm{1}}^Z}{Q}_{\mathrm{i}\mathrm{2}j}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{2}j}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{2}j}=\mathrm{0}\end{array}\right.（12）

显然，式（11）（12）的联立方程组中存在2Z+3个未知数，数量较多，可采用Newton-Raphson迭代法完成上述方程组的求解。调心滚子轴承接触载荷的总体计算流程如图4所示。

2 调心滚子轴承疲劳寿命计算

20世纪四十年代，G.Lundberg和A.Palmgrem经过理论分析和大量实验研究后，提出了滚动轴承疲劳寿命的计算方法，即著名的L⁃P理论。随着轴承材料和加工工艺的进步，L⁃P理论得到了修订和补充。经过一定简化后得到国际标准ISO 281，然而该标准忽略了各滚动体载荷不同的影响，只能进行近似的轴承寿命估算。因此，当能够获得准确的滚动体载荷分布时，应该优先采用L⁃P理论进行轴承疲劳寿命的计算。

对于具有点接触特征的双列调心滚子轴承，其疲劳寿命表达式为：

L={\left({L_{\mathrm{1}}}^{-\frac{\mathrm{10}}{\mathrm{9}}}+{L_{\mathrm{2}}}^{-\frac{\mathrm{10}}{\mathrm{9}}}\right)}^{-\frac{\mathrm{9}}{\mathrm{10}}}（13）

式中：L₁和L₂分别为第1、2列的滚子-滚道疲劳寿命，单位一般为百万转。

调心滚子轴承每一列的滚子-滚道寿命可按照单列点接触轴承的疲劳寿命公式计算，即：

L_{\mathrm{1}/\mathrm{2}}={\left({L_{\mathrm{1}\mu /\mathrm{2}\mu }}^{-\frac{\mathrm{10}}{\mathrm{9}}}+{L_{\mathrm{1}\nu /\mathrm{2}\nu }}^{-\frac{\mathrm{10}}{\mathrm{9}}}\right)}^{-\frac{\mathrm{9}}{\mathrm{10}}}（14）

式中：L_1μ 、L_2μ 和L_1ν 、L_2ν 分别为第1、2列滚子旋转和非旋转套圈的疲劳寿命。本文调心滚子轴承的旋转套圈是内圈滚道，非旋转套圈是外圈滚道。

内圈滚道和外圈滚道的疲劳寿命分别为：

\left\{\begin{array}{l}{L}_{\mathrm{1}\mu /\mathrm{2}\mu }={\left(\frac{Q_{\mathrm{1}c\mu /\mathrm{2}c\mu }}{Q_{\mathrm{1}e\mu /\mathrm{2}e\mu }}\right)}^{\mathrm{3}}\\ L_{\mathrm{1}\nu /\mathrm{2}\nu }={\left(\frac{Q_{\mathrm{1}c\nu /\mathrm{2}c\nu }}{Q_{\mathrm{1}e\nu /\mathrm{2}e\nu }}\right)}^{\mathrm{3}}\end{array}\right.（15）

式中：Q_1cμ、Q_2cμ 和Q_1cν、Q_2cν 分别为第1、2列滚子内圈和外圈滚道的额定动载荷；Q_1eμ、Q_2eμ 和Q_1eν、Q_2eν 分别为第1、2列滚子内圈和外圈滚道的当量载荷。

根据L⁃P理论，轴承内外圈滚道的额定动载荷由轴承的结构参数确定，故有Q_1cμ=Q_2cμ=Q_cμ，Q_1cν=Q_2cν=Q_cν，具体计算式为：

Q_{c\mu /c\nu }=\mathrm{98.1}{\left(\frac{\mathrm{2}R}{D}\frac{r_{\mathrm{i}/\mathrm{o}}}{r_{\mathrm{i}/\mathrm{o}}-R}\right)}^{\mathrm{0.41}}\frac{{(\mathrm{1}\mp \gamma )}^{\mathrm{1.39}}}{{(\mathrm{1}\pm \gamma )}^{\mathrm{1}/\mathrm{3}}}\times

{\left(\frac{\gamma }{\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{i}\mathrm{n}\mathrm{i}}}\right)}^{\mathrm{0.3}}{D}^{\mathrm{1.8}}{Z}^{-\frac{\mathrm{1}}{\mathrm{3}}}（16）

式中：上面和下面符号分别用于内、外圈。

对于与滚子实际载荷分布相关的当量载荷，其计算式为：

\left\{\begin{array}{l}{Q}_{\mathrm{1}e\mu /\mathrm{2}e\mu }={\left(\frac{\mathrm{1}}{Z}{\displaystyle \sum _{j=\mathrm{1}}^Z}{Q_{\mathrm{i}\mathrm{1}j/\mathrm{i}\mathrm{2}j}}^{\mathrm{3}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\\ Q_{\mathrm{1}e\nu /\mathrm{2}e\nu }={\left(\frac{\mathrm{1}}{Z}{\displaystyle \sum _{j=\mathrm{1}}^Z}{Q_{\mathrm{o}\mathrm{1}j/\mathrm{o}\mathrm{2}j}}^{\frac{\mathrm{10}}{\mathrm{3}}}\right)}^{\mathrm{0.3}}\end{array}\right.（17）

式中：Q_i1j 、Q_i2j 和Q_o1j 、Q_o2j 分别为第1、2列第j个滚子与内圈和外圈滚道的法向接触载荷。

由式（13）~（17），代入前文求解的滚子-滚道接触载荷值，可以得到疲劳寿命L。由于轴承的额定动载荷与其结构参数相关，需引入b_m修正系数来反映轴承材料和加工工艺的进步。对于调心滚子轴承取b_m=1.15，修正后的轴承疲劳寿命为：

L_{\mathrm{m}}={b_{\mathrm{m}}}^{\mathrm{3}}L（18）

3 基于AK-MCS的轴承疲劳寿命可靠性分析

由于加工和装配等误差，调心滚子轴承的几何结构和材料参数可能呈现一定的随机性，进而影响滚子接触载荷分析和疲劳寿命计算的准确性。以轴承的疲劳寿命是否达到一定的阈值为判别条件，建立极限状态函数为：

\mathit{G}=g(\mathit{X})=R(\mathit{X})-r^{\mathrm{*}}（19）

式中：R（ X ）为轴承在随机变量 X 影响下的寿命结果，本文确定随机变量为滚子轮廓半径R、滚子直径D、内圈滚道轮廓半径r_i、外圈滚道轮廓半径r_o、初始接触角α_ini和弹性模量E，即 X =［R，D，r_i，r_o，α_ini，E］^T；r^* 为寿命阈值，根据GL风机行业设计标准^［16］，风力发电机主轴轴承的运行寿命应达到20年，按主轴转速为10 r/min，换算结果为r^* =105.12×10⁶转。

当 G <0时，轴承疲劳寿命低于规定阈值，轴承处于失效域F=｛ x：g（ x ）≤0｝。失效概率为：

（20）

可靠性灵敏度定义为输入变量 X 分布参数的变化引起失效概率P_f 变化的比率，可衡量分布参数在给定点对失效概率的影响程度^［17］。具体表达式为：

\frac{\partial P_f}{\partial {\theta }_{{\mathit{X}}_i}^{\left(k\right)}}=\frac{\partial {\int }_F{f}_{\mathit{X}}\left(\mathit{x}\right)\mathrm{d}\mathit{x}}{\partial {\theta }_{{\mathit{X}}_i}^{\left(k\right)}}={\int }_F\frac{\partial f_{\mathit{X}}\left(\mathit{x}\right)}{\partial {\theta }_{{\mathit{X}}_i}^{\left(k\right)}}\mathrm{d}x（21）

式中：{\theta }_{{\mathit{X}}_{{}_i}}^{\left(k\right)}为 X_i 的第k个分布参数。

AK-MCS方法是自适应Kriging代理模型与Monte Carlo数字模拟法的结合，调用模型的次数和计算量更少，能得到更加准确的预测结果，有效提升可靠性分析的效率^{［17，18］}。主要步骤为：

（1）根据输入变量的概率密度函数，通过Monte Carlo模拟法产生MCS样本池 S_MC，该样本池由N_MC个样本组成。

（2）通过随机抽样法，从样本池 S_MC中选择N₁个样本，代入轴承疲劳寿命计算公式后得到式（20）极限状态函数值，构建初始训练集 T_MC。

（3）根据初始训练集 T_MC的数据建立Kriging代理模型g_K （ X ）。

（4）计算样本池 S_MC中每个样本的U学习函数值，其中U学习函数定义为：

U(\mathit{x})=\left|\frac{{\mu }_{g_K}(\mathit{x})}{{\sigma }_{g_K}(\mathit{x})}\right|（22）

当\underset{\mathit{x}\in {\mathit{S}}_{\mathrm{M}\mathrm{C}}}{\mathrm{m}\mathrm{i}\mathrm{n}}U(\mathit{x})\ge \mathrm{2}时，停止自适应学习过程，执行第（5）步；否则，更新样本点 x^u 为：

{\mathit{x}}^u=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{x}\in {\mathit{S}}_{\mathrm{M}\mathrm{C}}}{\mathrm{m}\mathrm{i}\mathrm{n}}U(\mathit{x})（23）

将｛ x^u，g（ x^u ）｝加入到 T_MC中，并更新第（3）步的g_K （ X ）。

（5）由当前Kriging代理模型g_K （ X ）进行可靠性分析，失效概率估计值及其变异系数分别为：

{\widehat{P}}_f=\frac{N_{g_K\le \mathrm{0}}}{N_{\mathrm{M}\mathrm{C}}}（24）

\mathrm{C}\mathrm{o}\mathrm{v}({\widehat{P}}_f)=\sqrt[]{\frac{\mathrm{1}-{\widehat{P}}_f}{(N_{\mathrm{M}\mathrm{C}}-\mathrm{1}){\widehat{P}}_f}}（25）

式中：N_gK≤0表示g_K （ x^（j））（j=1，2，…，N_MC）小于等于0的样本个数。当Cov（{\widehat{P}}_f）<5%时，认为{\widehat{P}}_f满足要求；否则，扩充样本池 S_MC后执行第（4）步。

（6）根据当前Kriging代理模型g_K （ X ）计算 S_MC每个样本 x^（j）（j=1，2，…，N_MC）的失效域指示函数值为：

I_F({\mathit{x}}^{(j)})=\left\{\begin{array}{l}\mathrm{0},g_K({\mathit{x}}^{(j)})>\mathrm{0}\\ \mathrm{1},g_K({\mathit{x}}^{(j)})\le \mathrm{0}\end{array}\right.({\mathit{x}}^{(j)}\in {\mathit{S}}_{\mathrm{M}\mathrm{C}})（26）

（7）计算可靠性灵敏度的估计值\frac{\partial {\widehat{P}}_f}{\partial {\theta }_{{\mathit{X}}_i}^{\left(k\right)}}及其变异系数\mathrm{C}\mathrm{o}\mathrm{v}\left(\frac{\partial {\widehat{P}}_f}{\partial {\theta }_{{\mathit{X}}_i}^{\left(k\right)}}\right)。

（8）当\mathrm{C}\mathrm{o}\mathrm{v}\left(\frac{\partial {\widehat{P}}_f}{\partial {\theta }_{{\mathit{X}}_i}^{\left(k\right)}}\right)<5%时，结束AK-MCS求解过程；否则，扩充样本池 S_MC后执行第（4）步。

4 数值算例

本文对某种2 MW风力发电机进行研究，主轴安装的调心滚子轴承型号为240/600 ECA/W33-SKF Explorer，如图5所示。具体结构参数见表1。

4.1　轴承接触载荷模型验证

兆瓦级风机一般在环境恶劣的地区；同时，将承受部件自身的重力、阵风的轴向冲击等，载荷复杂多变。本文按照风机的额定工况进行分析，主轴转速为10 r/min，轴承的轴向载荷为F_x =100 kN，径向载荷分别为F_y =700 kN和F_z =200 kN。

利用本文提出的理论模型与行业内的仿真软件Romax对滚子-滚道法向接触载荷进行计算，两者计算结果如图6所示。分析可知，第2列滚子的法向接触载荷整体上显著高于第1列，第1、2列的最大接触载荷分别为36.19 kN和128.58 kN。理论计算值与Romax仿真值之间的平均相对误差不超过0.1%，拟合效果较好，证明本文提出的调心滚子轴承滚子⁃滚道接触载荷模型具有较高的准确性。

由于风力载荷复杂多变，设轴向力F_x 变化范围为90~140 kN，滚子⁃滚道接触载荷计算结果如图7所示。随着轴向力F_x 的增大，第2列滚子的法向接触载荷整体上始终大于第1列，第1列最大接触载荷从41.94 kN不断减小至8.71 kN，第2列最大接触载荷从123.69 kN不断增加至141.11 kN。轴承内圈在轴向力的作用下造成两列滚子与滚道的接触状态呈现出明显不同，导致单列滚子承受偏载。

在风机实际运行过程中，主轴轴承的内外圈之间可能会因为轴的变形和安装误差等原因产生调心角α_s。对于调心滚子轴承，α_s取值为-2.5~2.5°。在额定工况下，当调心角α_s分别取不同值时，轴承滚子⁃滚道法向接触载荷的变化情况如图8所示。

从图8可知，当调心角为正值时，随着α_s的增大，第1列滚子的法向接触载荷逐渐降低，第2列滚子的载荷不断升高，滚子偏载的现象进一步增强；当调心角为负时，随着α_s绝对值的增大，第1列滚子的载荷呈现升高趋势，第2列滚子的载荷逐渐降低，滚子偏载的现象得到缓解。

4.2　疲劳寿命可靠性及灵敏度分析

本文研究的风电主轴调心滚子轴承的结构材料参数变量 X =［R，D，r_i，r_o，α_ini，E］^T相互独立，且符合正态分布，各变量的分布参数如表2所示。

在额定工况下，采用AK⁃MCS方法开展轴承的疲劳寿命可靠性及灵敏度分析。首先，根据随机参数的概率密度函数，通过Monte Carlo数字模拟法生成容量为4×10⁶组的样本池；之后，采用拉丁超立方抽样法，随机选取30组的变量样本，计算疲劳寿命L_m和极限状态函数值，以此构建Kriging代理模型并进行模型的自适应学习与迭代。

另外，随机抽取50组变量为测试样本，分别代入轴承疲劳寿命理论公式和自适应Kriging代理模型中计算，如图9所示。由图可知，自适应Kriging代理模型的预测值与理论值之间的相对误差在±0.15%内，证明本文建立的自适应Kriging代理模型具有较高的拟合精度。

为了研究轴向力的大小对轴承疲劳寿命可靠度的影响，利用AK-MCS方法进行可靠性分析，得到如图10所示的寿命可靠度曲线。当轴向力F_x =104 kN时，寿命可靠度为0.996 1，轴承开始发生失效，风机主轴轴承有较小的概率在20年内损坏；当F_x 超过110 kN时，轴承寿命可靠度随着F_x 的增大而急剧下降；当F_x 增大至133 kN时，寿命可靠度为0.001 1，轴承几乎完全失效，风机主轴轴承的寿命年限有相当大的概率不能达标。

图10显示轴承内圈受到的轴向力会降低轴承的寿命可靠度，而根据4.1节可知，轴向力和调心角的作用将影响两列滚子受力不均衡的程度，由此推测，轴承的寿命可靠度可能与两列滚子的偏载情况有关。在F_x =120 kN，F_y =700 kN和F_z =200 kN的条件下，设α_s取值为-2.5°~2.5°，利用AK⁃MCS方法计算得到如图11所示的轴承疲劳寿命可靠度曲线。

由图11可知，当轴承内圈不发生倾斜，即α_s =0时，轴承寿命可靠度为0.236 64，风机主轴轴承在寿命年限内失效的概率较大；当α_s 为负且绝对值增大时，寿命可靠度不断升高，直至α_s =-2°时，可靠度达到0.999 8，基本能满足运行寿命要求；而当α_s 以正值开始增大时，轴承寿命可靠度逐渐降低，直至α_s =1.5°时，可靠度降到0.000 5，主轴轴承几乎完全失效，不能正常工作。由此可知，在正向轴向力和径向力的作用下，内圈发生顺时针方向的倾斜将降低轴承的疲劳寿命可靠度，而逆时针方向的倾斜将有利于轴承寿命可靠度的提高。

为研究随机变量 X 各参数对轴承疲劳寿命失效概率的影响程度，利用AK-MCS进行可靠性灵敏度分析，如图12所示。分析可知，滚子轮廓半径R、滚子直径D和初始接触角α_ini与轴承的疲劳寿命失效概率成负相关，即轴承寿命可靠度随R、D和α_ini的增大而升高；内圈滚道轮廓半径r_i、外圈滚道轮廓半径r_o和弹性模量E与轴承的疲劳寿命失效概率成正相关，即轴承寿命可靠度随r_i、r_o和E的增大而降低。其中，滚子轮廓半径R对轴承的寿命失效概率的影响最大，占比达到59.62%；其次是内圈滚道轮廓半径r_i，占比为26.86%。因此，在调心滚子轴承的生产制造过程中，应注意严格控制滚子和内圈滚道轮廓半径的加工误差。

5 结束语

本文建立了风电调心滚子轴承接触载荷及疲劳寿命的计算模型。运用Romax软件进行仿真，验证了提出的模型具有较高的准确性。之后，采用AK-MCS方法对轴承的疲劳寿命可靠性进行了分析。最后，探究了部分外部因素及随机结构参数变量对轴承的接触载荷及寿命可靠性的影响。轴承在正向径向载荷的作用下，正向轴向力的增大和顺时针方向的调心角将加剧两列滚子的偏载现象，并且降低其疲劳寿命可靠度。轴承的寿命可靠度将随着滚子轮廓半径、滚子直径和初始接触角的增大而升高。随着内外圈滚道轮廓半径的增大而降低，其中滚子轮廓半径和内圈滚道轮廓半径对轴承寿命可靠度的影响较大。

参考文献

[1] 王同光, 秦海岩. 2014-2023年中国风电产业发展成就及未来展望[J]. 太阳能, 2024(7): 31-39.

[2] Wang Tong-guang, Qin Hai-yan. Development achievement and future prospect of wind power industry in China from 2014 to 2023[J]. Solar Energy, 2024(7): 31-39.

[3] Sinha R, Sahoo V, Paswan M. Radial load distribution by balls in a ball bearing with variable clearance[J]. Mechanics Based Design of Structures and Machines, 2023, 51(6): 3538-3563.

[4] Ohta H, Sakaguchi T, Uchiumi M. Load–displacement relationship of a ball bearing with axial, radial, and angular displacements for both the inner and outer rings[J]. Journal of Tribology, 2017, 139(1):No. 011103.

[5] Yang L, Xu T. Characteristics of high speed tapered roller bearing considering combined loads[C]∥Prognostics and System Health Management Conference, Harbin, China,2017: 1-7.

[6] Xu T, Yang L, Wu W, et al. Effect of angular misalignment of inner ring on the contact characteristics and stiffness coefficients of duplex angular contact ball bearings[J]. Mechanism and Machine Theory, 2021, 157: No.104178.

[7] Wu H, Zhou J, Li S. Theoretical and experimental analysis of contact stress in tapered roller bearings under preload[J]. Journal of Mechanical Science and Technology, 2016, 30(8): 3753-3762.

[8] Xu J W, Benson S, Wetenhall B. Comparative analysis of fatigue life of a wind turbine yaw bearing with different support foundations[J]. Ocean Engineering 2021, 235: 132-146.

[9] Bogdan W, Agnieszka C. Effect of ring misalignment on the fatigue life of the radial cylindrical roller bearing[J]. International Journal of Mechanical Sciences, 2016, 111: 464-474.

[10] Mostafa E L, Ling W, Terry J H, et al. Further understanding of rolling contact fatigue in rolling eleent bearings—a review[J]. Tribology International, 2019, 140: 644-691.

[11] Warda B, Chudzik A. Fatigue life prediction of the radial roller bearing with the correction of roller generators[J]. International Journal of Mechanical Sciences, 2014, 89: 299-310.

[12] Shi X, Wang L, Qin F. Relative fatigue life prediction of high-speed and heavy-load ball bearing based on surface texture[J]. Tribology International, 2016, 101: 364-374.

[13] 刘学婧, 刘磊, 徐华, 等. 滚子分布方式对球面滚子轴承服役性能的影响[J]. 西安交通大学学报, 2018, 52(12): 9-17.

[14] Liu Xue-jing, Liu Lei, Xu Hua, et al. Effect of rollers arrangement form on service performance of spherical roller bearing[J]. Journal of Xi'an Jiaotong University, 2018, 52(12): 9-17.

[15] 马方波. 脂润滑球面滚子轴承性能及润滑脂寿命预测方法研究[D]. 上海: 华东理工大学机械与动力工程学院, 2016.

[16] Ma Fang-bo. Study on the performance and grease life prediction method for grease-lubricated spherical roller bearings[D]. Shanghai: School of Mechanical and Power Engineering, East China University of Science and Technology, 2016.

[17] 罗继伟, 罗天宇. 滚动轴承分析计算与应用[M]. 北京: 机械工业出版社, 2009.

[18] Harris A T, Kotzalas N M. Essential Concepts of Bearing Technology(5th ed)[M]. Taylor and Francis:CRC Press, 2006.

[19] Germanischer L. Guideline for the certification of wind turbines [S]. Hamburg, Germany: Germanischer Lloyd, 2010: 221-223.

[20] 黄贤振, 孙楷铂, 栾晓刚, 等. 螺栓预紧连接可靠性灵敏度分析[J]. 吉林大学学报:工学版, 2023,53(8): 2219-2226.

[21] Huang Xian-zhen, Sun Kai-bo, Luan Xiao-gang, et al. Reliability sensitivity analysis of bolt pre-tightening connection[J]. Journal of Jilin University(Engineering and Technology Edition), 2023, 53(8): 2219-2226.

[22] Pei J Z, Ming W C, Zhou H Z, et al. A new active learning method based on the learning function U of the AK-MCS reliability analysis method[J]. Engineering Structures, 2017, 148: 185-194.

基金资助

国家自然科学基金项目(U23B2098)

国家自然科学基金项目(U22B2087)

辽宁省应用基础研究计划项目(2023JH2/101300160)