基于热网络法的风电机组主轴承温度场分析

Temperature Field Analysis of Main Bearings in Wind Turbines Based on Thermal Network Method

摘要

针对结构型式为三排圆柱滚子的风电机组主轴承的发热温升问题，建立了一种基于热网络法的轴承温度场分析方法。通过对轴承滚子进行切片处理，利用变形协调条件建立轴承的受力平衡方程组，求解得到滚子切片载荷，并对滚子切片与滚道表面之间的滑动运动进行分析；然后，根据轴承零件之间的载荷作用和相对滑动运动，得到滚子、滚道、保持架和密封圈相互接触部位的摩擦热源发热功率；最后，根据轴承的摩擦发热部位及内部具体结构，采用热网络法设置轴承内部的网络节点，并建立轴承的热网络平衡方程组。针对具体的2.0 WM风电机组主轴承，求解分析了润滑黏度、轴向游隙、滚子数量及转速对轴承温度的影响。

Abstract

Based on thermal network method，a bearing temperature field analysis method was established to address the heating and temperature rise problems of the main bearings in wind turbine with a three-row cylindrical roller structure. By slicing the rollers of the bearings， a force balance equation system of the bearings was established by using deformation coordination conditions， the roller slice loads were obtained by solving the equation systems. The sliding motions between the roller slices and the raceway surfaces were also analyzed. Then， based on the load actions and relative sliding motions between the bearing parts， the heat generation powers of the friction heat source at the contact positions of the rollers， raceways， retainers and sealing rings were obtained. Finally， according to the friction heating positions and specific internal structures of the bearings， the thermal network method was used to set up the internal network nodes of the bearings， and the thermal network balance equation system of the bearings was established. For a specific the main bearing in a 2.0 MW wind turbine， the effects of lubrication viscosity， axial clearance， roller number and rotational speed on bearing temperature were solved and analyzed.

Graphical abstract

Author summay

李云峰^*（通信作者），男，1973年生，教授。研究方向为滚动轴承的设计理论、性能分析与应用。 E-mail： [e-mail]。

0 引言

能源短缺和环境污染问题日益突出，开发利用可再生能源来替代传统的化石能源成为应对这一问题的有效途径。近几年，随着能源结构转型的加速，风力发电已成为发展最快的可再生能源利用方式之一。根据全球风能协会（GWEC）的预测^［1］，在未来5年内全球将新增风电装机容量791 GW，年均新增风电装机容量超过158 GW，风电机组的装机空间巨大。

主轴承是风电机组中的核心部件，承担着风力载荷和重力载荷的共同作用。主轴承的运行工况条件复杂多变，具有故障率高、维护成本高和国产化率低的特点，是风电装备技术领域的关注热点^［2-8］。在现有的风电机组中，主轴承的结构类型有多种，其中，三排圆柱滚子轴承具有承载能力高、工艺性能优、结构紧凑的优势。现有针对三排圆柱滚子轴承的相关研究工作主要集中在滚子载荷分布^［9-10］、套圈应力^［11］、套圈变形^［12］、滚道寿命^［13］等的分析计算。三排圆柱滚子轴承的轴向推力滚子不满足纯滚动的几何条件，滚子与滚道之间存在滑动摩擦。在连续回转运动过程中，滚子与滚道之间的滑动接触摩擦将导致轴承的发热和温升。温升进一步带来轴承套圈的热胀，这将造成轴承配合状态和载荷分布的改变；过量的温升还会造成润滑条件的恶化，这些变化将导致轴承疲劳寿命的缩短和可靠度的降低。控制温升是提高轴承运转可靠性的途径之一。轴承温升的有效控制建立在对轴承的热生成和热传导机理进行系统分析的基础之上。贾晨鹤等^［14］采用局部生热法计算了结构型式为双列圆锥滚子的风电机组主轴承内部各部位的摩擦生热功率，通过有限元仿真分析得到了轴承的温度场分布。向东等^［15］构建了风电齿轮箱的热网络模型，推导出网络节点的稳态热力学方程组和瞬态热力学方程组，求解分析了实际工况下的齿轮箱油温、高速轴轴承温度和齿轮箱入口油温。宋海等^［16］分别采用热网络法和数值仿真法分析了风电齿轮箱高速输出轴轴承的温度分布，并通过热网络法分析了不同工况下的轴承温度及其对轴承变形和接触应力的影响。高文君等^［17］提出了一种改进的多节点热网络模型，在轴承内部滚动体与滚道接触部位的油膜中心增加一个网络节点，分析了航空发动机主轴承的温度分布。孔亚楠等^［18］建立了磁悬浮保护轴承中各零件之间的生热及热传递关系，运用热网络法和有限元法分析了轴承的温度分布规律。

本文针对结构型式为三排圆柱滚子的风电机组主轴承的发热温升问题，建立一种基于热网络法的轴承温度场分析方法，并深入分析轴承参数变化对温升的影响规律，可为轴承的参数设计提供依据。在对轴承进行力学建模分析的基础上，得到滚子、滚道、保持架和密封圈相互接触部位的摩擦热源发热功率；在轴承的轴向截面内设置合理的节点分布，构建节点之间的热平衡方程组，通过求解计算得到轴承的温度场分布。

1 发热模型

1.1 力学平衡分析

三排圆柱滚子风电机组主轴承由套圈、滚子、保持架和密封圈组成。如图1所示，上、下风向的轴向滚子承受轴向载荷和倾覆力矩，径向滚子承受径向载荷。轴承的内圈连接在机架上固定不动，外圈随风轮转动。

当外部径向载荷F_{\mathrm{r}}、轴向载荷F_{\mathrm{a}}、倾覆力矩M共同作用于轴承外圈时，外圈将产生相应方向的径向位移{\delta }_{\mathrm{r}}、轴向位移{\delta }_{\mathrm{a}}和倾角位移\theta。

为便于研究轴向滚子与滚道之间的相对滑动，对滚子进行切片处理。将各排滚道中的每个滚子沿轴向分为n_{\mathrm{s}}个切片，则每个滚子切片的厚度为w=L_{\mathrm{w}}/n_{\mathrm{s}}，其中L_{\mathrm{w}}为滚子长度。

当外部载荷作用时，外圈的轴向位移{\delta }_{\mathrm{a}}和倾角位移\theta将产生轴向滚子与轴向滚道之间的接触变形。根据轴承的变形协调条件，上风向和下风向第j个滚子的第k个切片与滚道之间的接触变形量{\delta }_{j,k}^{\mathrm{u}}、{\delta }_{j,k}^{\mathrm{d}}分别为

{\delta }_{j,k}^{\mathrm{u}}={\delta }_{\mathrm{a}}+\{\mathrm{0.5}{d}_{\mathrm{m}}^{\mathrm{u}}-[\mathrm{0.5}{L}_{\mathrm{w}}^{\mathrm{u}}-w^{\mathrm{u}}(k-\mathrm{0.5})]\}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j^{\mathrm{u}}\mathrm{t}\mathrm{a}\mathrm{n}\theta -\mathrm{0.5}{u}_{\mathrm{a}}-c_{j,k}^{\mathrm{u}}（1）

{\delta }_{j,k}^{\mathrm{d}}=-{\delta }_{\mathrm{a}}-\{\mathrm{0.5}{d}_{\mathrm{m}}^{\mathrm{d}}-[\mathrm{0.5}{L}_{\mathrm{w}}^{\mathrm{d}}-w^{\mathrm{d}}(k-\mathrm{0.5})]\}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j^{\mathrm{d}}\mathrm{t}\mathrm{a}\mathrm{n}\theta -\mathrm{0.5}{u}_{\mathrm{a}}-c_{j,k}^{\mathrm{d}}（2）

式中：d_{\mathrm{m}}^{\mathrm{u}}、d_{\mathrm{m}}^{\mathrm{d}}分别为上风向和下风向滚子的分布圆直径；L_{\mathrm{w}}^{\mathrm{u}}、L_{\mathrm{w}}^{\mathrm{d}}分别为上风向和下风向滚子的长度；w^{\mathrm{u}}、w^{\mathrm{d}}分别为上风向和下风向滚子切片的厚度；{\psi }_j^{\mathrm{u}}、{\psi }_j^{\mathrm{d}}分别为上风向和下风向第j个滚子的圆周方位角；c_{j,k}^{\mathrm{u}}、c_{j,k}^{\mathrm{d}}分别为上风向和下风向第j个滚子的第k个切片处的凸度修形量；u_{\mathrm{a}}为轴承的轴向游隙。

同理，当外部载荷作用时，外圈的径向位移{\delta }_{\mathrm{r}}和倾角位移\theta将产生径向滚子与径向滚道之间的接触变形。径向第j个滚子的第k个切片与滚道之间的接触变形量{\delta }_{j,k}^{\mathrm{r}}可表示为

{\delta }_{j,k}^{\mathrm{r}}={\delta }_{\mathrm{r}}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j^{\mathrm{r}}+\frac{\mathrm{1}}{\mathrm{2}}{w}^{\mathrm{r}}(k-\frac{\mathrm{1}}{\mathrm{2}})\theta \mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j^{\mathrm{r}}-\frac{u_{\mathrm{r}}}{\mathrm{2}}-c_{j,k}^{\mathrm{r}}（3）

式中：w^{\mathrm{r}}为径向滚子切片的厚度；{\psi }_j^{\mathrm{r}}为径向第j个滚子的圆周方位角；c_{j,k}^{\mathrm{r}}为径向第j个滚子的第k个切片处的凸度修形量；u_{\mathrm{r}}为轴承的径向游隙。

根据接触力与接触变形之间的关系，可得上风向和下风向第j个滚子的第k个切片与滚道之间的法向接触载荷q_{j,k}^{\mathrm{u}}、q_{j,k}^{\mathrm{d}}分别为^［19］

q_{j,k}^{\mathrm{u}}=\frac{w^{\mathrm{u}}({\delta }_{j,k}^{\mathrm{u}}{)}^{\mathrm{1.11}}}{(f_{j,k}^{\mathrm{u}}{)}^{\mathrm{1.11}}{A}^{\mathrm{1.11}}(L_{\mathrm{w}}^{\mathrm{u}}{)}^{\mathrm{1.11}}}（4）

q_{j,k}^{\mathrm{d}}=\frac{w^{\mathrm{d}}({\delta }_{j,k}^{\mathrm{d}}{)}^{\mathrm{1.11}}}{(f_{j,k}^{\mathrm{d}}{)}^{\mathrm{1.11}}{A}^{\mathrm{1.11}}(L_{\mathrm{w}}^{\mathrm{d}}{)}^{\mathrm{1.11}}}（5）

式中：A=\mathrm{1.36}{\eta }^{\mathrm{0.9}}，\eta为两接触体的综合弹性常数；f_{j,k}^{\mathrm{u}}、f_{j,k}^{\mathrm{d}}分别为上风向和下风向第j个滚子的第k个切片与滚道之间接触变形量的修正系数。

径向滚子第j个滚子的第k个切片与滚道之间的法向接触载荷q_{j,k}^{\mathrm{r}}可表示为^［19］

q_{j,k}^{\mathrm{r}}=\frac{w^{\mathrm{r}}({\delta }_{j,k}^{\mathrm{r}}{)}^{\mathrm{1.11}}}{(f_{j,k}^{\mathrm{r}}{)}^{\mathrm{1.11}}{A}^{\mathrm{1.11}}(L_{\mathrm{w}}^{\mathrm{r}}{)}^{\mathrm{1.11}}}（6）

式中：f_{j,k}^{\mathrm{r}}为径向第j个滚子的第k个切片与滚道之间接触变形量的修正系数；L_{\mathrm{w}}^{\mathrm{r}}为径向滚子的长度。

根据轴承的力学平衡条件，轴承外圈在径向载荷F_{\mathrm{r}}、轴向载荷F_{\mathrm{a}}、倾覆力矩M和滚子切片载荷q_{j,k}^{\mathrm{u}}、q_{j,k}^{\mathrm{d}}、q_{j,k}^{\mathrm{r}}的共同作用下处于平衡状态，由此可得如下平衡方程组：

F_{\mathrm{a}}={\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{u}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}{q}_{j,k}^{\mathrm{u}}-{\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{d}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}{q}_{j,k}^{\mathrm{d}}（7）

F_{\mathrm{r}}={\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{r}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}{q}_{j,k}^{\mathrm{r}}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j^{\mathrm{r}}（8）

M={\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{u}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}{q}_{j,k}^{\mathrm{u}}\{\frac{d_{\mathrm{m}}^{\mathrm{u}}}{\mathrm{2}}-[\frac{L_{\mathrm{w}}^{\mathrm{u}}}{\mathrm{2}}-w^{\mathrm{u}}(k-\frac{\mathrm{1}}{\mathrm{2}})]\}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j^{\mathrm{u}}-

{\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{d}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}{q}_{j,k}^{\mathrm{d}}\{\frac{d_{\mathrm{m}}^{\mathrm{d}}}{\mathrm{2}}-[\frac{L_{\mathrm{w}}^{\mathrm{d}}}{\mathrm{2}}-w^{\mathrm{d}}(k-\frac{\mathrm{1}}{\mathrm{2}})]\}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j^{\mathrm{d}}（9）

式中：Z^{\mathrm{u}}为上风向滚子的数量；Z^{\mathrm{d}}为下风向滚子的数量；Z^{\mathrm{r}}为径向滚子的数量。

上述方程组可通过Newton-Raphson迭代方法求解未知\mathrm{位}移量{\delta }_{\mathrm{r}}、{\delta }_{\mathrm{a}}和\theta，并进一步得到各排滚子的各切片所承受的载荷。

1.2 滚子运动分析

由于轴向滚子与滚道接触各点所处位置的回转半径不同，在保持架窗孔的约束作用下，滚子与滚道之间在接触部位将产生相对滑动。假设轴向滚子与滚道之间在滚道节圆处为纯滚动，而在接触线上其余各点存在相对滑动。如图2所示，O点为轴承的回转中心，{\omega }_{\mathrm{i}}、{\omega }_{\mathrm{e}}分别为轴承的内圈和外圈转动角速度，{\omega }_{\mathrm{c}}为保持架的公转角速度，滚子在保持架窗孔中随保持架一起公转。

根据滚道节圆处纯滚动运动的假设，可得到上风向和下风向滚子的自转角速度{\omega }_{\mathrm{s}}^{\mathrm{u}}、{\omega }_{\mathrm{s}}^{\mathrm{d}}分别为

{\omega }_{\mathrm{s}}^{\mathrm{u}}=\frac{d_{\mathrm{m}}^{\mathrm{u}}}{\mathrm{2}{D}_{\mathrm{w}}^{\mathrm{u}}}({\omega }_{\mathrm{i}}-{\omega }_{\mathrm{e}})（10）

{\omega }_{\mathrm{s}}^{\mathrm{d}}=\frac{d_{\mathrm{m}}^{\mathrm{d}}}{\mathrm{2}{D}_{\mathrm{w}}^{\mathrm{d}}}({\omega }_{\mathrm{i}}-{\omega }_{\mathrm{e}})（11）

式中：D_{\mathrm{w}}^{\mathrm{u}}、D_{\mathrm{w}}^{\mathrm{d}}分别为上风向和下方向滚子的直径。

上风向和下风向保持架的公转角速度{\omega }_{\mathrm{c}}^{\mathrm{u}}、{\omega }_{\mathrm{c}}^{\mathrm{d}}分别为

{\omega }_{\mathrm{c}}^{\mathrm{u}}=({\omega }_{\mathrm{i}}-{\omega }_{\mathrm{e}})/\mathrm{2}（12）

{\omega }_{\mathrm{c}}^{\mathrm{d}}=({\omega }_{\mathrm{i}}-{\omega }_{\mathrm{e}})/\mathrm{2}（13）

上风向第j个滚子上的第k个切片与内圈滚道和外圈滚道之间的滑动速度\mathrm{\Delta }{v}_{j,k}^{\mathrm{u}}、\mathrm{\Delta }{\mu }_{j,k}^{\mathrm{u}}分别为

\mathrm{\Delta }{v}_{j,k}^{\mathrm{u}}=({\omega }_{\mathrm{i}}-{\omega }_{\mathrm{c}}^{\mathrm{u}})y_{j,k}^{\mathrm{u}}-\mathrm{0.5}{D}_{\mathrm{w}}^{\mathrm{u}}{\omega }_{\mathrm{s}}^{\mathrm{u}}（14）

\mathrm{\Delta }{\mu }_{j,k}^{\mathrm{u}}=({\omega }_{\mathrm{e}}-{\omega }_{\mathrm{c}}^{\mathrm{u}})y_{j,k}^{\mathrm{u}}-\mathrm{0.5}{D}_{\mathrm{w}}^{\mathrm{u}}{\omega }_{\mathrm{s}}^{\mathrm{u}}（15）

式中：y_{j,k}^{\mathrm{u}}为上风向第j个滚子的第k个切片厚度中心到轴承中心的距离。

下风向第j个滚子上的第k个切片与内圈滚道和外圈滚道之间的滑动速度\mathrm{\Delta }{v}_{j,k}^{\mathrm{d}}、\mathrm{\Delta }{\mu }_{j,k}^{\mathrm{d}}分别为

\mathrm{\Delta }{v}_{j,k}^{\mathrm{d}}=({\omega }_{\mathrm{i}}-{\omega }_{\mathrm{c}}^{\mathrm{d}})y_{j,k}^{\mathrm{d}}-\mathrm{0.5}{D}_{\mathrm{w}}^{\mathrm{d}}{\omega }_{\mathrm{s}}^{\mathrm{d}}（16）

\mathrm{\Delta }{\mu }_{j,k}^{\mathrm{d}}=({\omega }_{\mathrm{e}}-{\omega }_{\mathrm{c}}^{\mathrm{d}})y_{j,k}^{\mathrm{d}}-\mathrm{0.5}{D}_{\mathrm{w}}^{\mathrm{d}}{\omega }_{\mathrm{s}}^{\mathrm{d}}（17）

式中：y_{j,k}^{\mathrm{d}}为下风向第j个滚子的第k个切片厚度中心到轴承中心的距离。

1.3 摩擦热源分析

轴承的发热功率在数值上等于摩擦力与速度的乘积或摩擦力矩与转速的乘积。三排圆柱滚子风电机组主轴承的摩擦来源主要包括滚子与滚道之间的滑动摩擦和滚动摩擦、保持架与套圈之间的滑动摩擦、密封圈与套圈之间的滑动摩擦、滚子与润滑脂之间的搅动摩擦。

1） 滚子与滚道之间存在相对滑动，由上风向和下风向滚子与滚道之间相对滑动摩擦产生的发热功率H_{\mathrm{s}}^{\mathrm{u}}、H_{\mathrm{s}}^{\mathrm{d}}分别为

H_{\mathrm{s}}^{\mathrm{u}}={\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{u}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}\frac{(\mathrm{\Delta }{v}_{j,k}^{\mathrm{u}}+\mathrm{\Delta }{\mu }_{j,k}^{\mathrm{u}})}{h^{\mathrm{*}}}{\int }_{-\mathrm{2}b}^{\mathrm{0}}{\eta }_{\mathrm{m}}\mathrm{d}x（18）

H_{\mathrm{s}}^{\mathrm{d}}={\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{d}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}\frac{(\mathrm{\Delta }{v}_{j,k}^{\mathrm{d}}+\mathrm{\Delta }{\mu }_{j,k}^{\mathrm{d}})}{h^{\mathrm{*}}}{\int }_{-\mathrm{2}b}^{\mathrm{0}}{\eta }_{\mathrm{m}}\mathrm{d}x（19）

式中：h^{*}为最小油膜厚度；b为滚子与滚道之间的Hertz接触半宽；{\eta }_{\mathrm{m}}为切片压力区内对应位置的润滑脂黏度。

2） 滚子与滚道之间的存在弹性接触变形，由于上风向和下风向滚道材料存在弹性滞后而产生的发热功率H_{\mathrm{t}}^{\mathrm{u}}、H_{\mathrm{t}}^{\mathrm{d}}分别为^［20］

H_{\mathrm{t}}^{\mathrm{u}}={\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{u}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}[\xi \sqrt[]{\frac{\mathrm{\pi }{q}_{j,k}^{\mathrm{u}}}{\mathrm{2}\eta D_{\mathrm{w}}^{\mathrm{u}}}}{\delta }_{j,k}^{\mathrm{u}}(|{\omega }_{\mathrm{i}}-{\omega }_{\mathrm{c}}^{\mathrm{u}}|+|{\omega }_{\mathrm{e}}-{\omega }_{\mathrm{c}}^{\mathrm{u}}|)y_{j,k}^{\mathrm{u}}]（20）

H_{\mathrm{t}}^{\mathrm{d}}={\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{d}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}[\xi \sqrt[]{\frac{\mathrm{\pi }{q}_{j,k}^{\mathrm{d}}}{\mathrm{2}\eta D_{\mathrm{w}}^{\mathrm{u}}}}{\delta }_{j,k}^{\mathrm{d}}(|{\omega }_{\mathrm{i}}-{\omega }_{\mathrm{c}}^{\mathrm{d}}|+|{\omega }_{\mathrm{e}}-{\omega }_{\mathrm{c}}^{\mathrm{d}}|)y_{j,k}^{\mathrm{d}}]（21）

由于径向滚道材料存在弹性滞后而产生的发热功率H_{\mathrm{t}}^{\mathrm{r}}为^［20］

H_{\mathrm{t}}^{\mathrm{r}}={\displaystyle \sum _{j=\mathrm{1}}^{Z^{\mathrm{r}}}}{\displaystyle \sum _{k=\mathrm{1}}^{n_{\mathrm{s}}}}[\xi \sqrt[]{\frac{\mathrm{\pi }{q}_{j,k}^{\mathrm{r}}}{\mathrm{2}\eta D_{\mathrm{w}}^{\mathrm{u}}}}{\delta }_{j,k}^{\mathrm{r}}(|{\omega }_{\mathrm{i}}-{\omega }_{\mathrm{c}}^{\mathrm{r}}|+|{\omega }_{\mathrm{e}}-{\omega }_{\mathrm{c}}^{\mathrm{r}}|)y_{j,k}^{\mathrm{r}}]（22）

式中：\xi为材料的弹性滞后系数；\eta为两接触体的综合弹性常数。

3） 轴承在工作运转时，套圈通过配合引导面对保持架进行回转定位，保持架引导面与套圈的接触滑动摩擦产生的发热功率H_{\mathrm{c}}为^［21］

H_{\mathrm{c}}=\frac{\mathrm{1}}{\mathrm{2}}{m}_{\mathrm{c}}g{f}_{\mathrm{c}}{d}_{\mathrm{c}}{\omega }_{\mathrm{c}}（23）

式中：m_{\mathrm{c}}为保持架质量；g为重力加速度；f_{\mathrm{c}}为套圈与引导面的摩擦因数；d_{\mathrm{c}}为保持架引导面直径。

4） 双唇密封圈与内圈挡边之间通过安装过盈量产生弹性接触力以达到密封效果，在内圈旋转过程中的滑动摩擦产生发热功率H_{\mathrm{s}\mathrm{l}}为

H_{\mathrm{s}\mathrm{l}}=\frac{\mathrm{3}\kappa EI}{\mathrm{2}{L}_{\mathrm{s}\mathrm{l}}^{\mathrm{3}}}\mathrm{\pi }{f}_{\mathrm{s}\mathrm{l}}{d}_{\mathrm{s}\mathrm{l}}^{\mathrm{2}}{\omega }_{\mathrm{i}}（24）

式中：\kappa为密封唇过盈量；E为密封唇材料的弹性模量；I为单位长度密封唇的截面惯性矩；L_{\mathrm{s}\mathrm{l}}为密封唇伸出部分的宽度；f_{\mathrm{s}\mathrm{l}}为密封唇接触面的摩擦因数；d_{\mathrm{s}\mathrm{l}}为密封唇接触面的回转直径。

5） 滚子运转时搅动润滑脂所产生的发热功率H_{\mathrm{\upsilon }}为^［22］

H_{\mathrm{\upsilon }}={\mathrm{10}}^{-\mathrm{7}}{f}_{\mathrm{0}}({\upsilon }_{\mathrm{0}}{n)}^{\frac{\mathrm{2}}{\mathrm{3}}}{d}_{\mathrm{m}}^{\mathrm{3}}{\omega }_{\mathrm{c}}（25）

式中：f_{\mathrm{0}}为轴承结构类型系数；{\upsilon }_{\mathrm{0}}为润滑油黏度；n为轴承的转速；d_{\mathrm{m}}为轴承的节圆直径。

2 热网络模型

2.1 热网络法

热网络法借助电学上的基尔霍夫节点电流定律，将热学问题比拟为电学问题，根据研究对象的结构及热源，将需要获得温度的部位设置为节点，节点之间的热量传递（导热、对流、辐射）均通过传热热阻进行连接，从而形成热网络。当整个系统达到稳态温度状态时，任意节点上输入的热量等于输出的热量，以此建立节点的热平衡方程组。通过求解节点的线性方程组来获得各个节点的稳态温度。

根据稳态时任意节点v流入热量等于流出热量，该节点的净热流量等于0，可得

{\displaystyle \sum _{\mu =\mathrm{1}}^{N_v}}\frac{T_{\mu }-T_v}{R_{\mu ,v}}+Q_v=\mathrm{0}　 v=\mathrm{1,2},\cdot \cdot \cdot ,N（26）

式中：T_{\mu }为第\mu个节点的温度；T_v为第v个节点的温度；R_{\mu ,v}为第\mu个节点与第v个节点之间的热阻；Q_v为节点v的内热源；N_v为与节点v存在热量交换关系的节点个数；N为整个系统的节点总数。

假设无传热关系的节点之间热阻无限大，则整个系统的热平衡方程组可用矩阵形式表示为

\mathit{R}\mathit{T}=\mathit{Q}（27）

式中：\mathit{R}为热阻矩阵，由各个节点之间的传热热阻组成；\mathit{T}为待求解的节点温度矩阵；\mathit{Q}为热源矩阵。

式（27）中，\mathit{T}和\mathit{Q}均为N维列向量，对于整个轴承系统，除去有摩擦发热的部位外，其余节点处的发热量均为0。通过求解该矩阵即可求得所有节点温度。轴承热网络模型的建立主要包括如下三个步骤：

1）节点设置。根据轴承的结构及热源分布在轴承系统中设置不同的网络节点，以便将系统的热传递行为进行离散化建模。

2）构建闭合回路。通过传热热阻连接这些网络节点，构建一个闭合的热流回路，以确保热量在系统内部能够连续传递。

3）建立热平衡方程组。基于传热学和基尔霍夫定律，为每个网络节点建立热平衡方程，构建整个系统网络节点的热平衡方程组，求解方程组得到整个系统的稳态温度分布。

2.2 热阻计算

基于轴承内部各部位节点位置分布的特点，考虑各节点之间的热量传递时，也考虑热传导和热对流的作用，忽略热辐射的影响。根据傅里叶导热定律和牛顿冷却公式，在一维稳态传热过程中传导热热阻计算公式为

R_{\mathrm{c}}=L/(\lambda A)（28）

式中：L为导热物体特征长度；\lambda为材料的热导率；A为换热面积。

对流换热热阻的计算公式为

R_{\mathrm{t}}=\mathrm{1}/(hA)（29）

式中：h为对流传热系数。

将轴承内部各部分均简化为薄圆环、平壁或圆柱体模型，再根据热阻定义及材料和尺寸参数，求得各节点间的导热热阻。对于对流换热，其传热系数可表示为

h=Nu{\lambda }_{\mathrm{a}}/L（30）

式中：Nu为努塞尔数，取决于换热方式；{\lambda }_{\mathrm{a}}为空气的热导率。

由于轴承的内圈固定、外圈旋转，因此，将外圈与空气之间视为强迫对流换热，而内圈与空气之间视为自然对流换热。强迫对流换热时努塞尔数的计算公式为^［23］

N{u}_{\mathrm{f}}=\mathrm{0.133}R{e}^{\mathrm{2}/\mathrm{3}}P{r}^{\mathrm{1}/\mathrm{3}}（31）

式中：Re为空气的雷诺数；Pr为空气的普朗特数。

自然对流换热时努塞尔数的计算公式为^［23］

N{u}_{\mathrm{n}}=C{(GrPr)}^{\tau }（32）

式中：Gr为格拉晓夫常数；C、\tau为热力学常数，根据Gr的计算值确定。

2.3 温度节点划分

在对三排圆柱滚子风电机组主轴承系统结构及热源进行分析的基础上，根据结构各组成部分之间的传热关系，在轴承的轴向截面内设置热网络的47个关键温度节点。温度节点分布如图3所示，其中空心点表示不含热源的温度节点，实心点表示含热源的温度节点。设置环境温度T_{\mathrm{0}}和润滑脂温度T_{\mathrm{o}\mathrm{i}\mathrm{l}}作为已知边界条件。

轴承零件之间的摩擦所产生的发热量按照1∶1的比例分配到相邻的接触部位，并将上述划分的网络节点按照传热关系进行连接。根据设定的热流网络通过热阻进行热传递，可以得到主轴承的热网络模型如图4所示。图4中，长方形表示两个相邻节点之间的热阻，空心点表示不含热源的节点，实心点表示热源节点，其中，节点9、13、15、18、19、20、24、25、26为滚子-滚道热源节点，节点12、14、22、28、32、34为保持架热源节点，节点31、35为密封圈热源节点。根据建立的热网络模型，由基尔霍夫定律可知每个节点的热流量守恒。

根据式（26）构建节点的热平衡方程，将滚子与套圈产生的热量以1∶1的方式分配到相邻网络节点，得到轴承系统的热平衡方程组如下：

\left.\begin{array}{l}\frac{T_{\mathrm{1}}-T_{\mathrm{0}}}{R_{\mathrm{1,0}}}+\frac{T_{\mathrm{1}}-T_{\mathrm{4}}}{R_{\mathrm{1,4}}}=\mathrm{0}\\ \frac{T_{\mathrm{2}}-T_{\mathrm{7}}}{R_{\mathrm{2,7}}}+\frac{T_{\mathrm{2}}-T_{\mathrm{3}}}{R_{\mathrm{2,3}}}+\frac{T_{\mathrm{2}}-T_{\mathrm{0}}}{R_{\mathrm{2,0}}}=\mathrm{0}\\ ⋮\\ \frac{T_{\mathrm{46}}-T_{\mathrm{40}}}{R_{\mathrm{46,40}}}+\frac{T_{\mathrm{46}}-T_{\mathrm{45}}}{R_{\mathrm{46,45}}}+\frac{T_{\mathrm{46}}-T_{\mathrm{47}}}{R_{\mathrm{46,47}}}+\frac{T_{\mathrm{46}}-T_{\mathrm{0}}}{R_{\mathrm{46,0}}}=\mathrm{0}\\ \frac{T_{\mathrm{47}}-T_{\mathrm{41}}}{R_{\mathrm{47,41}}}+\frac{T_{\mathrm{47}}-T_{\mathrm{46}}}{R_{\mathrm{47,46}}}+\frac{T_{\mathrm{47}}-T_{\mathrm{0}}}{R_{\mathrm{47,0}}}=\mathrm{0}\end{array}\right\}（33）

其中，T为温度，R为热阻。依据上述建立的47个节点平衡方程，将轴承结构参数及各部分产生的摩擦热功率代入方程组（式（33）），通过对方程组进行求解即可得到各温度节点的温度值。

3 实例计算

本文针对2.0 MW风电机组某型号三排圆柱滚子风电机组主轴承进行温度场分析，主轴承的结构参数如表1所示。

主轴承的转速为\mathrm{10.2}\mathrm{r}/\mathrm{m}\mathrm{i}\mathrm{n}，润滑脂的基础油为PAO40，运动黏度为396 \mathrm{m}{\mathrm{m}}^{\mathrm{2}}/\mathrm{s}，动力黏度为\mathrm{0.386}\mathrm{P}\mathrm{a}\bullet \mathrm{s}。轴承承受的当量载荷为：径向载荷F_{\mathrm{r}}=\mathrm{245}\mathrm{k}\mathrm{N}，轴向载荷F_{\mathrm{a}}=\mathrm{632}\mathrm{k}\mathrm{N}，倾覆力矩M=\mathrm{3278}\mathrm{k}\mathrm{N}\bullet \mathrm{m}。

3.1 热网络法

基于上述发热模型和热网络模型，通过编写程序来实现轴承温度场的数值求解，温度场的求解结果如图5所示。

从图5中轴承各节点的稳态温度分布数据中可以发现，轴承系统中各节点的温度明显高于环境常温，其内圈的平均温度为76.4 ℃，外圈的平均温度为71.0 ℃，内圈平均温度高于外圈平均温度（高5.4 ℃），原因在于内圈与空气之间为自然对流换热；滚子与滚道产生的摩擦发热功率要远大于其余部位产生的摩擦发热功率，且散热条件差，这导致滚子的最高温度可达到79.8 ℃，明显高于其余部位的温度；轴承的轴向滚子的最高温度高于径向滚子的最高温度（高4.1 ℃），轴向滚子的最高温度出现在上风向滚子处，这是因为上风向滚子承受着较大的轴向风力载荷，所产生的摩擦功率也最大。

3.2 有限元分析法

为验证上述热网络模型的合理性，通过商业有限元软件进行轴承的温度场分析。考虑到轴承结构沿圆周方向完全一致的特点，针对1个滚子直径厚度的轴承切片进行建模分析。在SolidWorks软件中建立该轴承切片的三维几何模型，将该三维模型导入到ANSYS Workbench中的Steady-State Thermal模块进行温度场分析。先完成轴承内外圈及滚子材料热导率参数的设置，在完成材料设置后对模型进行网格划分，并对滚子进行网格细化，划分后的网格模型如图6所示。

完成网格划分后对模型的边界条件进行设置。根据式（10）～式（14）计算出轴承滚道各部位的摩擦发热功率，在滚子、滚道、内圈挡边的摩擦接触发热部位添加热流载荷，并在滚子表面及轴承套圈外表面设置对应的对流传热系数，然后进行轴承的温度场分布求解，求解结果如图7所示。

由轴承热网络法与有限元法求解结果得到的温度场数据对比如表2所示。可以发现两种方法的计算结果十分接近，这验证了本文所建立的热网络模型的合理性和正确性。

图8为实际运行风电机组的主轴承温度监控曲线，时间跨度为27 d，每隔10 min采样一次，测温部位为轴承外圈。由于受到风速和风向随机变化的影响，风电机组主轴的转速具有波动性，且在有些时段低于额定转速，导致轴承的温度也表现为起伏变化。图9为主轴承温度数据的频数分布直方图，温度分布数据主要集中在60～70 ℃范围内。该数据与本文图5中热网络法的温度场计算结果一致，进一步验证了本文所建立的热网络模型的合理性和正确性。

3.3 温度的影响因素分析

由三排圆柱滚子风电机组主轴承的结构及各节点的温度分布可知，轴承中温度最高且最易发生温度升高的节点位于轴向滚子处，因此，在分析不同因素对轴承温度的影响时，将轴向滚子处的温度作为轴承的温度指标。

3.3.1 润滑油黏度对温度的影响

结合工程实际，分析所针对的润滑脂基础油动力黏度变化范围为0.\mathrm{2}~\mathrm{0.8}\mathrm{P}\mathrm{a}\bullet \mathrm{s}。图10为轴承工况不变时，轴向滚子温度随润滑油黏度的变化曲线，上风向滚子的温升变化幅度达到24.9 ℃，下风向滚子的温升变化幅度达到24.1 ℃。由图10可知，随着润滑油黏度的不断增大，轴承温度不断升高，这是因为滚子滚道的滑动摩擦会受到润滑状态的影响。润滑油脂黏度的增大导致油膜的剪切应力增大，滑动摩擦功率增大，以至于轴承温度随之升高。由此可知，在选择润滑油时，需要在温度范围内保证轴承处于良好润滑状态下，且尽可能地选择黏度变化小的润滑油以保证轴承处于较为稳定的运转状态。

3.3.2 轴向游隙对温度的影响

轴向游隙对温度的影响如图11所示，上风向滚子的温升变化幅度达到30.3 ℃，下风向滚子的温升变化幅度达到32.0 ℃。可以看出，轴向游隙的逐渐增大使得轴承温度呈现上升的趋势，当轴向游隙处于-0.03 mm左右时，轴承温度最低。随着轴承游隙的增大，起初轴承处于较大负游隙时，轴承内圈滚子与滚道因为过分压紧而导致滚子载荷较大，摩擦功率也同样较大，温度较高。当轴向游隙增大至-0.03 mm左右时，承载区滚子数量不变而单个滚子载荷减小，此时摩擦功率最小，温度最低，伴随着轴向游隙的进一步增大，承载区滚子数量开始减少，单个滚子的摩擦功率增大，温度开始攀升。

3.3.3 滚子数量对温度的影响

通过改变轴向滚子数量可以得到轴承温度变化，如图12所示，上风向滚子的温升变化幅度达到15.7 ℃，下风向滚子的温升变化幅度达到14.5 ℃。可以看出，随着滚子数量的增加，轴承的温度均呈现先快速下降、后缓慢下降的趋势。这是因为随着滚子数量的增加，单个滚子承受的载荷迅速减小，产生的摩擦功率也在迅速减小，而当滚子数量增加到一定数量后，承载区的滚子承受载荷变化不再明显，摩擦功率的减小幅度变小，温度下降趋势变缓。

3.3.4 转速对温度的影响

轴承转速对温度的影响如图13所示，上风向滚子的温升变化幅度达到50.4 ℃，下风向滚子的温升变化幅度达到50.3 ℃。可以看出，转速的增大使得轴承温度呈现快速上升的趋势。这是因为转速的增大使得上下风向的滚子与滚道之间的摩擦功率有明显的增大；虽然套圈与周围空气的对流传热系数也有所增大，从表面流失的热量有所增加，但套圈表面的散热功率远小于滚道的摩擦发热功率，轴承温度随着转速增大呈现上升趋势。

4 结论

1）轴承的温度分布不均匀，内圈平均温度比外圈平均温度高5.4 ℃，上风向滚子的最高温度达到79.8 ℃，明显高于其余部位的温度。轴承转速的增大会导致轴承温度显著升高，当转速从\mathrm{8}~\mathrm{20}\mathrm{r}/\mathrm{m}\mathrm{i}\mathrm{n}变化时，温升变化幅度高达50.4 ℃，这种温升现象在上风向滚子处表现得更为明显。

2）轴向游隙的增大使得轴承温度出现先下降后上升的趋势，当轴向游隙处于-0.03 mm左右时，轴承温度最低，且比游隙为0时的温度低11.7 ℃。轴向游隙通过影响滚子载荷和载荷分布范围来影响轴承的摩擦力矩。因此，存在一个最佳负游隙范围使得轴承的温度最低。

3）滚子数量的增加会使得轴承温度出现先明显降低，而当滚子数量增加到一定数量后，承载区的滚子承受载荷变化不再明显，摩擦发热功率的减小幅度变小，温度下降趋势变缓。当滚子数量由110增加至130时，轴承温度降低4.5 ℃。

参考文献

[1] Global Wind Energy Council. Global Wind Report 2024［EB/OL］.（2024-04-16）［2024-05-23］.

[2] LIU Zepeng， ZHANG Long. A Review of Failure Modes， Condition Monitoring and Fault Diagnosis Methods for Large-scale Wind Turbine Bearings［J］. Measurement， 2020， 149： 107002.

[3] 刘杰， 蒋树旗. 基于LSTNet网络的风电机组发电机主轴承温度预测研究［J］. 传感技术学报， 2022， 35（7）： 945-951.

[4] LIU Jie， JIANG Shuqi. Temperature Prediction of Main Bearing of Wind Turbine Generator Based on LSTNet［J］. Chinese Journal of Sensors and Actuators， 2022， 35（7）： 945-951.

[5] 尹诗， 侯国莲， 胡晓东， 等. 基于AC-GAN数据重构的风电机组主轴承温度监测方法［J］. 智能系统学报， 2021， 16（6）： 1106-1116.

[6] YIN Shi， HOU Guolian， HU Xiaodong， et al. Temperature Monitoring Method of the Main Bearing of Wind Turbine Based on AC-GAN Data Reconstruction［J］. CAAI Transactions on Intelligent Systems， 2021， 16（6）： 1106-1116.

[7] 向健平， 凌永志， 詹俊， 等. 基于SCADA系统的风电机组主轴承故障预警方法［J］. 电力科学与技术学报， 2019， 34（3）： 223-228.

[8] XIANG Jianping， LING Yongzhi， ZHAN Jun， et al. Fault Early Warning Method of Wind Turbine Main Bearing Based on SCADA System［J］. Journal of Electric Power Science and Technology， 2019， 34（3）： 223-228.

[9] 王洪斌， 王红， 何群， 等. 基于深度信念网络的风机主轴承状态监测方法［J］. 中国机械工程， 2018， 29（8）： 948-953.

[10] WANG Hongbin， WANG Hong， HE Qun， et al. Condition Monitoring Method for Wind Turbine Main Bearings Based on DBN［J］. China Mechanical Engineering， 2018， 29（8）： 948-953.

[11] 何群， 王红， 江国乾， 等. 基于相关主成分分析和极限学习机的风电机组主轴承状态监测研究［J］. 计量学报， 2018， 39（1）： 89-93.

[12] HE Qun， WANG Hong， JIANG Guoqian， et al. Research of Wind Turbine Main Bearing Condition Monitoring Based on Correlation PCA and ELM［J］. Acta Metrologica Sinica， 2018， 39（1）： 89-93.

[13] 赵洪山，刘辉海. 基于深度学习网络的风电机组主轴承故障检测［J］. 太阳能学报， 2018， 39（3）： 588-595.

[14] ZHAO Hongshan， LIU Huihai. Fault Detection of Wind Turbine Main Bear Based on Deep Learning Network［J］. Acta Energiae Solaris Sinica， 2018， 39（3）： 588-595.

[15] GÖNCZ P， POTOFINIK R， GLODEŽ S. Computational Model for Determination of Dynamic Load Capacity of Large Three-row Roller Slewing Bearings［J］. Engineering Failure Analysis， 2013， 32： 44-53.

[16] ZHAO Haiming， ZHANG Tao， XIONG Zhihong， et al. Influence Analysis of Key Design Parameters on the Roller Safety Factor of Three-row and Four-column Roller Slewing Bearings Based on Analytical Method and FEM［J］. Journal of the Brazilian Society of Mechanical Sciences and Engineering， 2022， 44： 488.

[17] LI Yunfeng， JIANG Di. Strength Check of a Three-row Roller Slewing Bearing Based on a Mixed Finite Element Model［J］. Journal of Mechanical Engineering Science， 2017， 231（18）： 3393-3400.

[18] HE Peiyu， WANG Yun， WANG Hua. An Analysis Method of Carrying Capacity Accuracy of Three-row Roller Slewing Bearing［J］. Mechanika， 2021， 27（5）： 360-367.

[19] 李云峰， 高元安， 贾磊. 表面硬化滚道风电主轴承的疲劳寿命分析［J］. 太阳能学报， 2022， 43（6）： 176-182.

[20] LI Yunfeng， GAO Yuanan， JIA Lei. Analysis for Fatigue Life of Wind Turbine Main Bearing with Surface Hardened Raceways［J］. Acta Energiae Solaris Sinica， 2022， 43（6）： 176-182.

[21] 贾晨鹤， 陈观慈， 王燕霜， 等. 风电机组主轴用双列 圆锥滚子轴承的热分析［J］. 轴承， 2023（6）： 43-50.

[22] JIA Chenhe， CHEN Guanci， WANG Yanshuang， et al. Thermal Analysis of Double Row Tapered Roller Main Bearings for Wind Turbines［J］. Bearing， 2023（6）： 43-50.

[23] 向东， 韦尧中， 沈银华， 等. 面向风电齿轮箱油温超 限故障的热网络模型及其节点温度计算方法研究［J］. 机械工程学报， 2022， 58（9）： 119-135.

[24] XIANG Dong， WEI Yaozhong， SHEN Yinhua， et al. Research on Thermal Network Modeling and Temperature Calculation Method for Wind Turbine Gearbox Lubrication Oil Temperature Overrun Fault［J］. Journal of Mechanical Engineering， 2022， 58（9）： 119-135.

[25] 宋海， 穆塔里夫·阿赫迈德， 尼加提·玉素甫， 等. 基于热网络与数值仿真模拟的风电轴承瞬态温度研究［J］. 轴承， 2023（6）： 51-58.

[26] SONG Hai， MUTELLIP Ahmat， NIJAT Yusup， et al. Transient Temperature Analysis of Wind Turbine Bearings Based on Thermal Network and Numerical Simulation［J］. Bearing， 2023（6）： 51-58.

[27] 高文君， 刘振侠， 朱鹏飞， 等. 基于多节点热网络法 的航空发动机主轴承温度场分析［J］. 推进技术， 2019（2）： 382-388.

[28] GAO Wenjun， LIU Zhenxia， ZHU Pengfei， et al. Steady Thermal Analysis of Main-shaft Roller Bearing for Aero-engine Based on Multi-nodes Thermal Network Methods［J］. Journal of Propulsion Technology， 2019（2）： 382-388.

[29] 孔亚楠， 薛玉君， 叶军， 等. 热网络法和有限元法保 护轴承热分析研究［J］. 机械设计与制造， 2023（10）： 31-35.

[30] KONG Yanan， XUE Yujun， YE Jun， et al. Research on Thermal Analysis of Auxiliary Bearings Based on Thermal Network Method and Finite Element Method［J］. Machinery Design and Manufacturing， 2023（10）： 31-35.

[31] 邓四二， 李猛， 卢羽佳， 等. 推力滚针轴承摩擦力矩特性研究［J］. 兵工学报， 2015， 36（7）：1347-1355.

[32] DENG Sier， LI Meng， LU Yujia， et al. Study of Friction Torque of Needle Roller Thrust Bearing［J］. Acta Armamentarii， 2015， 36（7）： 1347-1355.

[33] 张占立， 周鹏举， 李文博， 等. YRT转台轴承摩擦力矩特性研究［J］. 兵工学报， 2019， 40（7）： 1495-1502.

[34] ZHANG Zhanli， ZHOU Pengju， LI Wenbo， et al. Study of Friction Torque Characteristics of YRT Rotary Table Bearing［J］. Acta Armamentarii， 2019， 40（7）： 1495-1502.

[35] 张占立， 张开哲， 王恒迪， 等. 球柱联合转台轴承摩擦力矩特性［J］. 航空动力学报， 2020， 35（12）： 2664-2672.

[36] ZHANG Zhanli， ZHANG Kaizhe， WANG Hengdi， et al. Friction Torque Characteristic of Ball-roller Combined Turntable Bearing［J］. Journal of Aerospace Power， 2020， 35（12）： 2664-2672.

[37] HARRIS T A. Rolling Bearing Analysis［M］. 4th ed. New York： John Wiley & Sons， Inc.， 2001.

[38] 杨世铭，陶文铨. 传热学［M］. 4版.北京： 高等教育出版社， 2007.

[39] YANG Shiming， TAO Wenquan.Heat Transfer ［M］.4th ed.Beijing： Higher Education Press， 2007.