风电圆锥滚子主轴承接触轮廓修型分析

Contact Profile Modification Analysis of Tapered Roller Main Bearing in Wind Turbine

摘要: 针对风力发电机圆锥滚子主轴承，研究7 MW风电主轴承多工况的接触机理。通过雨流计数法对主轴承载荷谱进行等效力分析，根据静力学平衡方程计算轴承内部的滚子载荷分布并通过有限长赫兹接触理论计算滚子与滚道的接触应力。进一步研究不同滚子修型对滚子与滚道之间的接触应力的影响。研究表明：不同轮廓形貌的母线修型和同轮廓下的修型凸度对滚子轮廓应力有着明显影响。

Abstract: The contact mechanism of 7 MW wind turbine main bearing with multiple working conditions is investigated for wind turbine tapered roller main bearing. The load spectrum of the main bearing is analyzed by the rain-flow counting method, the roller load distribution inside the bearing is calcu-lated according to the hydrostatic equilibrium equation, and the contact stresses between the roller and raceway are calculated by the finite-length Hertz contact theory. The effect of different roller profiles on the contact stresses between the rollers and the raceways is further investigated. The study shows that the different profiles of busbar modification and the convexity of modification un-der the same profile have a significant effect on the roller profile stresses.

1. 引言

随着传统能源的枯竭和能源开发带来的环境问题日益恶化，风能作为一种天然、可靠、清洁的能源异军突起，在能源市场上需求量很大，在新旧能源过渡中占据着不可或缺的地位 [1] [2] 。主轴轴承是风电机组的关键机械受载支承部件，主轴轴承的要求寿命在20年以上，但实际主轴轴承在服役早期就形成一定程度磨损。风机主轴轴承的运行工况非常复杂，许多研究人员针对影响主轴轴承运行的因素做了一定的研究。

郑林征 [3] 建立双列圆锥滚子主轴承的力学分析模型，计算极限和等效工况下的接触载荷与接触应力，并研究了轴承游隙对轴承疲劳寿命的影响。褚景春等 [4] 对双列调心滚子主轴承的运行中偏载问题，设计了非对称滚道结构，通过理论计算与有限元分析表明了改进的设计对偏载的问题有缓解效果。杨剑飞等 [5] 建立双列球面滚子轴承在联合载荷作用下的受力与位移方程，结合有限长线接触理论分析滚子与滚道的接触应力，并探究了球面半径等因素对接触应力的影响。肖文 [6] 建立了风机主轴双列圆锥滚子轴承模型，研究了极限工况下的倾覆力矩对风电轴承接触载荷与变形分布的影响规律。Jiang [7] 建立了一种新的考虑轴承间隙和预紧力的五自由度圆锥滚子轴承模型，并应用在2个单列圆锥滚子轴承支撑的风电主轴承系统中，最后以有限元模型、Romax模型和Harris模型的仿真结果进行了比较，验证了所提模型的有效性。Li [8] 研究了双列圆锥滚子轴承在外加载荷和角不对中作用下的力矩、接触载荷和疲劳寿命结果表明了外载荷、轴向预变形、轴承转速和角不对中对双列圆锥滚子轴承接触和疲劳特性有显著的影响。Zheng [9] 研究了在振荡外载荷和转速作用下的风机主轴的双列圆锥滚子轴承力学模型，计算了滚子的内载荷分布，然后采用基于弹性接触模型的迭代法研究了滚轮的接触压力。

本文旨在针对风电机组主轴的圆锥滚子轴承组合，建立一种计算轴承内部滚子载荷分布的5自由度力学模型，借助于静力平衡方程计算滚子内部载荷分布；在此基础上，运用有限长线接触理论对滚子作切片处理求解滚子与滚道之间的接触应力分布。研究滚子修型对滚道接触应力的影响。

2. 圆锥滚子主轴承接触模型

2.1. 主轴承结构与载荷谱处理

风电机组服役条件恶劣，主轴承受载复杂，需同时承载叶片轮毂所受风载荷和主轴重量，在收到轴向力(F_x)和径向力(F_y、F_z)作用的同时，还收到倾覆力矩(M_y和M_z)作用，如图1(a)所示，此外，风速的大小、方向在实时变化，所以主轴承时序载荷谱呈现无规律性。通过雨流计数法处理主轴承时序载荷谱，得到等效载荷。本文是以7 MW风电机组双支撑主轴承组合结构，如图1(b)所示。表1中的轴承参数与载荷谱均由江苏京冶海上风电轴承制造有限公司提供。

(a) 主轴承系统坐标 (b) 主轴承结构

Figure 1. Wind turbine main bearing system coordinates and structure

图1. 风电主轴承系统坐标和主轴承结构

Table 1. Structural parameters of main bearing

表1. 主轴承结构尺寸

载荷谱中的F_x作雨流技术分析，可以看出原本的F_x走势是随着时间作不规则变化，如图2(a)所示，不便于后续对主轴承作静力学分析时的作用力加载。经过雨流计数法编译后的F_x呈现出均值和幅值的变化，可以确定F_x的等效力，如图2(b)所示。以此重复计算3种工况载荷谱内的F_x、F_y、F_z、M_y和M_z，并记录表2。

(a) 原始载荷谱 (b) 雨流化处理

Figure 2. Compilation of rain-flow counting method for F_x load spectra

图2. F_x载荷谱的雨流计数法编译

Table 2. Main bearing equivalent loads for 3 operating conditions

表2. 3种工况的主轴承等效载荷

2.2. 主轴承静力学模型

由于风机主轴承处于低速运转，忽略掉滚子自重和转到引起的离心力，只考虑内外滚道和挡边的接触力，如图3所示。

Figure 3. Main bearing force diagram

图3. 主轴承受力图

对于双列圆锥滚子主轴承整体而言，5个自由度受力分析如下。

$F_x={\displaystyle \underset{i=1}{\overset{Z}{\sum }}\left(Q_{1ei}-Q_{2ei}\right)}\times \sin a_e$ (1)

$F_y={\displaystyle \underset{i=1}{\overset{Z}{\sum }}\left(Q_{1ei}+Q_{2ei}\right)}\times \cos a_e\cos \varphi$ (2)

$F_z={\displaystyle \underset{i=1}{\overset{Z}{\sum }}\left(Q_{1ei}+Q_{2ei}\right)}\times \cos a_e\sin \varphi$ (3)

$M_y=\frac{d_m}{2}{\displaystyle \underset{i=1}{\overset{Z}{\sum }}\left(Q_{1ei}-Q_{2ei}\right)}\times \sin a_e\sin \varphi$ (4)

$M_z=\frac{d_m}{2}{\displaystyle \underset{i=1}{\overset{Z}{\sum }}\left(Q_{1ei}-Q_{2ei}\right)}\times \sin a_e\cos \varphi$ (5)

式中Q_i、Q_e、Q_f为接触载荷，下标i、e、f分别表示滚子–内滚道接触、滚子–外滚道接触和滚子–挡边接触，1、2则表示滚道列数，d_m是节圆直径。

其中再以单个滚子作受力分析，得到接触载荷Q_i、Q_e、Q_f的关系如下：

$Q_i=Q_e\frac{\sin \left({\alpha }_e+{\alpha }_f\right)}{\sin \left({\alpha }_i+{\alpha }_f\right)}$ (6)

$Q_f=Q_e\frac{\sin \left({\alpha }_e-{\alpha }_f\right)}{\sin \left({\alpha }_i+{\alpha }_f\right)}$ (7)

当求出外接触载荷Q_e时，可以通过三者之间的关系求解内接触载荷Q_i和挡边接触载荷Q_f的值。由于风机主轴承为重载接触，在倾覆力矩作用下，轴承内外圈将产生一定的倾角θ，使得滚子接触呈现一定的变形量。根据弹性接触理论，第i个滚子与外滚道之间的接触载荷与变形量之间的关系为：

$Q_{1ei}=k_{ne}{\delta }_{1ni}^{1.11}$ (8)

$Q_{2ei}=k_{ne}{\delta }_{2ni}^{1.11}$ (9)

其中 $k_{ne}$ 是外滚道接触处总刚度系数，根据LUNDBERG的研究结果：

$k_{ne}=8.06\times {10}^4{l}^{0.89}{\left\{1+{\left[\frac{\sin \left({\alpha }_e+{\alpha }_f\right)}{\sin \left({\alpha }_i+{\alpha }_f\right)}\right]}^{0.89}\cos \left({\alpha }_e-{\alpha }_i\right)\right\}}^{-1.11}$ (10)

其中，l为滚子的有效接触长度。

第i个滚子的总位移：

${\delta }_{ni}={\delta }_{ri}\cos {\alpha }_e+{\delta }_{ai}\sin {\alpha }_e$ (11)

其中 ${\delta }_{ri}$ 是滚子的径向位移， ${\delta }_{ai}$ 是滚子的轴向位移：

${\delta }_{ri}=\left({\delta }_r+0.5\theta C\right)\cos \phi$ (12)

${\delta }_{ai}={\delta }_a+0.5d_m\theta \cos \phi$ (13)

以上公式中， ${\delta }_r$ 、 ${\delta }_a$ 和 $\theta$ 为未知量，先由几何关系得到第i个滚子法向位移 ${\delta }_{ni}$ ，接着通过接触总刚度，得到第i个滚子的接触载荷，再通过力学平衡方程迭代求解出 ${\delta }_r$ 、 ${\delta }_a$ 和 $\theta$ ，最后反求接触载荷。结合表2的主轴承结构尺寸和表1的3种工况负载，求解出主轴承滚子分别与外滚道、内滚道、挡边的接触载荷。

(a) 滚子–外滚道 (b) 滚子–内滚道 (c) 滚子–挡边

Figure 4. Contact load distribution of main bearing

图4. 主轴承接触载荷分布

根据建立的双列圆锥滚子轴承5自由度模型，求解得到主轴承内部载荷分布，如图4所示。以工况一为例，从计算结果中可以看出，滚子–外滚道接触中，左列最大载荷为296 KN，右列最大载荷为211 KN；滚子–内滚道接触中，左列最大载荷为296 KN，右列最大载荷为211 KN；滚子–挡边接触中，左列最大载荷为5195 N，右列最大载荷为3605 N。当主轴承受力矩作用时，第一列和第二列同位置的滚子存在较大的载荷差异。外滚道载荷Q_e、内滚道载荷Q_i、挡边载荷Q_f存在比例关系，所以三者只存在大小上的差别，在整圈分布中有着相同的分布规律。

3. 主轴承滚子接触应力

3.1. 有限长线接触赫兹应力分析

圆锥滚子主轴承滚子与滚道之间的接触属于线接触，在受重载的情况下，滚子接触面发生微变形而形成细长矩形接触痕迹，如图5(a)所示。

(a) 滚子接触示意图 (b) 有限元切片

Figure 5. Schematic diagram of roller contact stresses

图5. 滚子接触应力示意图

在滚子与滚道的接触截面中，应力轮廓为椭圆形，压力分布函数遵守椭圆轮廓，将接触区域沿长度方向划分为微小单元条n_ij，其长度与宽度2h × 2b，每个单元条上的接触压力分布为：

$q\left(x,y\right)=q_m\left(x\right)\sqrt{1-{\left(\frac{y}{b}\right)}^2},\left|y\right|\le b$ (14)

在微小单元条接触区域n_ij中心点(x_ij, y_ij)上，接触变形的计算公式可以表示为：

$\omega \left(x_{ij},y_{ij}\right)=\frac{1-{\nu }^2}{\pi E}{\displaystyle \underset{\Xi }{\iint }\frac{q_m\left(\tilde{x}\right)B\left(\tilde{x},\tilde{y}\right)}{\sqrt{{\left(x_{ij}-\tilde{x}\right)}^2+{\left(y_{ij}-\tilde{y}\right)}^2}}}\text{d}\tilde{x}\text{d}\tilde{y}$ (15)

将式(15)代入两个物体接触变形协调条件后，得到：

${\delta }_0={\omega }_1\left(x_{ij},y_{ij}\right)+{\omega }_2\left(x_{ij},y_{ij}\right)+z_1\left(x_{ij},y_{ij}\right)+z_2\left(x_{ij},y_{ij}\right)$ (16)

式中 $ij=1,2,3,\cdots ,MN$ ，MN为计算点数量。其中式(15)有MN个， ${\delta }_0$ 也未知，需要增加半椭圆分布下的接触平衡方程：

$Q={\displaystyle \underset{n_{ij}}{\sum }{\displaystyle \underset{{\Xi }_n}{\iint }{q}_m\left(\tilde{x}\right)\sqrt{1-{\left[\tilde{y}/b\left(\tilde{x}\right)\right]}^2}\text{d}\tilde{x}\text{d}\tilde{y}}}={\displaystyle \underset{n_{ij}}{\sum }{q}_m\left(x_{ij}\right)\frac{\pi 2h{b}_{ij}}{2}}$ (17)

(a) 滚子轮廓应力数值解 (b) 滚子轮廓应力仿真

Figure 6. Stress distribution on roller contact profile

图6. 滚子接触轮廓应力分布

如上图6(a)所示，将滚子与外滚道、挡边的接触载荷代入到上面各式中，迭代得到直线型滚子的接触应力分布。可以看出直线型滚子的端面应力急速上升，存在应力集中现象，滚子中部应力分布平缓。对滚子进行Workbench应力仿真分析，发现同样存在端部应力集中现象。应力集中使得滚子与滚道和挡边之间的滚动磨损变大，疲劳寿命降低。故此对主轴承滚子轮廓进行修型处理，使滚子轮廓应力分布更平缓。

3.2. 滚子修型分析

常见的滚子轮廓修型有部分圆弧直线曲线、对数型曲线和大圆弧曲线，每种修型适应不同的轴承应用。如图7所示，当滚子轮廓采用圆弧直线修型时，圆弧曲线与直线光滑，接触其母线方程为：

$z\left(x\right)=\{\begin{array}{l}0,\left|x\right|\le l_0/2\\ R_C-\sqrt{R_C{}^2-{\left(x-l_0/2\right)}^2},l_0/2\text{<}\left|x\right|\le l_e/2\end{array}$ (18)

式中l₀为滚子直母线长度，圆弧半径 $R_C=\left[{\left(l_w-2c\right)}^2-{\left(l_w-2t\right)}^2\right]/\left(8{\delta }_0\right)$ ，c为滚子端面倒角尺寸，t为凸度修正圆弧部分长度， ${\delta }_0$ 为凸度值。

Figure 7. Circular linear roller profile

图7. 圆弧直线型滚子轮廓

如图8所示，当滚子轮廓为对数型时，为使母线为光滑完整的对数曲线，只在滚子的端部区域采用对数曲线修型，而在滚子中部仍采用直线，对数曲线通常要根据滚子的凸度值来设计。这时滚子母线函数形式为：

$z\left(x\right)=\{\begin{array}{l}0,\left|x\right|\le l_0/2\\ \xi {\delta }_0\ln \frac{1}{1-\eta {\left[2\left(x-l_0/2\right)l_e\right]}^2},l_0/2<\left|x\right|\le l_e/2\end{array}$ (19)

Figure 8. Logarithmic roller profile

图8. 对数型滚子轮廓

其中，l₀为滚子直母线长度， $\xi$ ， $\eta$ 为系数。当 $x=l_0/2$ 时， $z=0$ ；当 $x=l_e/2$ 时， $z={\delta }_0$ 。利用这些条件可以确定函数中的系数 $\xi$ ， $\eta$ 的值。对数曲线型母线滚子如图：

如图9所示，当滚子轮廓为大圆弧型时，其母线方程为：

$z\left(x\right)=R_W-\sqrt{R_W{}^2-x^2},\left|x\right|\le l_e/2$ (20)

Figure 9. Large arc roller profile

图9. 大圆弧型滚子轮廓

(a) 圆弧直线型(b) 对数型 (c) 大圆弧型

Figure 10. Numerical solution for contact stresses of trimmed rollers

图10. 修型滚子接触应力数值解

(a) 圆弧直线型(b) 对数型 (c) 大圆弧型

Figure 11. Simulation of contact stresses on trimmed rollers

图11. 修型滚子接触应力仿真

对主轴承轮廓修型滚子进行接触应力数值求解，结果如图10所示，圆弧直线修型滚子应力分布呈“凹”字型，即在圆弧与直线交接处存在着较大应力，中间段应力分布较为平缓；对数型修型滚子应力分布平缓且应力值大小适中，有效缓解了圆弧直线型滚子轮廓交界处的应力峰值；但随着滚子整体凸度的逐渐增大，大圆弧型滚子中间段应力显著提高，两端应力下降，整体呈“凸”字型。同时在Workbench中对修型滚子进行仿真分析，滚子应力云图如图11所示，各修型滚子的应力云图分布与应力数值解无较大差异，进一步证明了风力发电机圆锥滚子主轴承滚子修型轮廓的差异性影响了整个轴承滚子的接触性能。

4. 结论

建立7MW双支承风机双列圆锥滚子主轴承接触应力分析模型，研究多种载荷工况对主轴承滚子接触载荷和滚子轮廓应力的影响，进一步探究不同滚子修型对滚子轮廓应力影响；

1) 为解决双支承风机双列圆锥滚子主轴承在5个自由度受力联合作用下的内部滚子载荷分布问题，建立双列圆锥滚子轴承的5自由度计算模型。从计算结果可以看出第一、二列滚子在不同工况作用下的接触载荷存在明显差异；外滚道、内滚道和挡边与滚子的接触载荷相同的分布规律和相似的比例关系。

2) 对于双支承风机双列圆锥滚子主轴承而言，滚子的母线修型显著影响滚子–滚道的轮廓应力分布。无修型滚子两端有明显的应力集中，圆弧直线型存在局部较大应力；大圆弧修型会导致滚子中端应力过大，而对数型滚子应力分布平缓，大小适中，结合风电主轴承的应用情况，推荐主轴承滚子采用对数修型。

3) 本文研究了滚子修型对轮廓应力的影响，后期将细分每个修型凸度，研究主轴承滚子轮廓应力随凸度增加的变化关系；同时对主轴承疲劳寿命展开深入研究。

致谢

首先要感谢我的导师杨超君老师，在论文撰写的过程中，我遇到了许多挑战和困难，杨老师都以严谨的态度和专业的知识给予了我无私的指导和帮助。其次感谢江苏京冶海上风电轴承制造有限公司的相关人员，他们不仅提供本篇论文的轴承参数和载荷谱，更为我解决风电主轴承的相关难题提供了宝贵的经验指导。论文已经完稿，但惭愧地说，风电轴承的实际工程问题与学术研究有着较大的区别，这让我意识到自己的不足，同时也让我对日后的学术研究和工作充满敬畏！

参考文献