电机轴承动力学特性分析及疲劳寿命计算

Dynamic Characteristics Analysis and Fatigue Life Calculation of Motor Bearings

摘要: 深沟球轴承是电机中的关键部件，明确轴承的动力学特性对轴承疲劳寿命预测和改进设计具有重要意义。基于接触力学和流体力学，文章结合深沟球轴承内部各部件之间的几何变形与受力情况，建立了轴承动力学模型，根据滚珠、套圈和保持架之间的协调变形关系，分析了轴承的动力学特性。首先，以深沟球轴承为对象，结合轴承几何结构和内部相互作用关系，建立了轴承动力学模型；然后，以轴承工作中的典型工况作为轴承动力学模型的仿真输入数据；其次，采用四阶龙格–库塔法和牛顿–拉夫逊法对动力学模型进行求解；最后，基于L-P理论，对轴承疲劳寿命进行计算。结果表明，保持架打滑率的稳定速率与转速呈反相关、与载荷呈正相关，低转速高载荷工况下的打滑率能更快地趋于稳定；滚珠与保持架之间的碰撞力在启动阶段震荡明显，且与载荷呈正相关，高载荷低转速工况下的碰撞力能更早地进入稳定状态；轴承寿命与转速和载荷呈反相关关系。

Abstract: Deep groove ball bearings are key components in electric motors, and clarifying their dynamic characteristics is of great significance for bearing fatigue life prediction and design improvement. Based on contact mechanics and fluid mechanics, combined with the geometric deformation and force conditions among the internal components of the deep groove ball bearing, a bearing dynamic model was established. This model analyzes the dynamic characteristics of the bearing according to the coordinated deformation relationship among the balls, rings, and cage. First, taking the deep groove ball bearing as the subject and combining its geometric structure and internal interaction relationships, the bearing dynamic model was established. Then, typical working conditions during bearing operation were used as input data for the bearing dynamic model simulation. Next, the fourth-order Runge-Kutta method and the Newton-Raphson method were used to solve the dynamic model. Finally, the bearing fatigue life was calculated based on the L-P theory. The results show that the stabilization rate of the cage slip rate is inversely related to the rotation speed and positively related to the load, meaning the slip rate stabilizes faster under low-speed, high-load conditions. The collision force between the ball and the cage clearly oscillates during the starting phase and is positively related to the load, with the collision force entering a stable state earlier under high-load, low-speed conditions. The bearing life is inversely related to the rotation speed and the load.

1. 引言

电机是电驱动系统中的主要动力源，而轴承是电机中的重要传动部件，轴承的性能对电驱动系统的稳定性有直接影响[1] [2]，因此有必要对电机轴承展开动力学特性分析。

目前对轴承动力学特性的仿真分析主要包括轴承有限元软件仿真和代码仿真。于洁[3]基于滚动轴承动力学、摩擦学和传热学，建立了球轴承的完全拟动力学模型，利用该模型分析了轴承结构参数对轴承动态性能和温升的影响。陈世金[4]建立了方柱形、圆柱形和球形兜孔保持架球轴承动力学模型，分析了不同类型的兜孔对轴承动力学特性的影响。分析指出，球形兜孔保持架的轴向摆动范围在三种类型兜孔中最小，但磨损率略高。赵志强[5]建立了轴承力学和弹流润滑的耦合模型，补充了润滑油膜较厚时轴承处于油膜润滑和润滑油膜较薄时轴承处于干摩擦之间的摩擦力计算，进一步反应真实的润滑状态。于拂晓[6]建立了深沟球轴承的拟静力学模型，并使用牛顿–拉夫逊法进行求解，分析了转速和载荷对轴承内部载荷的影响规律。

刘合涛[7]利用ABAQUS软件对深沟球轴承三维模型进行了动力学仿真，分析了轴承滚珠在运转过程中的应力分布状态和滚珠与保持架之间的载荷大小，得到了应力主要集中在滚珠和保持架接触部分以及接触载荷随时间呈周期性变化的结论。王赵蕊佳[8]使用有限元仿真软件分析了转速、载荷及润滑油供油量对轴承温度的影响，结果表明转速和载荷与轴承温度呈正相关关系，供油量与轴承温度呈反相关。

本文基于接触力学和流体力学，建立了深沟球轴承的动力学模型，并采用四阶龙格–库塔法和牛顿–拉夫逊法进行求解，利用求解结果分析了轴承的动力学特性，并开展了疲劳寿命评估。

2. 轴承动力学模型建立

2.1. 动力学模型坐标系

建立的深沟球轴承模型由外圈、内圈、滚珠和保持架组成，轴承整体坐标系O-XYZ以O为原点；假设外圈固定，内圈随传动轴旋转；内圈坐标系O₂-X₂Y₂Z₂以O₂为原点，伴随内圈的质心位移产生X、Y和Z三个方向上的位移[x₂; y₂; z₂]和绕Y轴和Z轴的偏转角θ_y和θ_z；对每个滚珠设置独立的坐标系，滚珠坐标系O_b-X_bY_bZ_b以球心O_b为原点，使用滚珠方位角θ_j表示滚珠在运动过程中所处位置。轴承及坐标系如图1所示。

Figure 1. Schematic diagram of the bearing and coordinate systems

图1. 轴承及坐标系示意图

2.2. 套圈与滚珠相互作用模型

轴承转动时，外、内套圈与滚珠相互作用，对滚珠产生接触应力和拖动力。套圈与滚珠的相对位置如图2所示。

Figure 2. Schematic diagram of raceway and rolling element positions. (a) Schematic diagram of the relative position of the outer ring raceway and rolling element, (b) Schematic diagram of the relative position of the inner ring raceway and rolling element

图2. 滚道和滚珠相对位置。(a) 外圈滚道和滚珠相对位置，(b) 内圈滚道和滚珠相对位置

图2中，后缀j为滚珠编号；O_rj为套圈滚道中心；r_bj为轴承整体坐标系中心到滚珠中心的距离；r_rj为轴承整体坐标系中心到套圈滚道中心的距离。外圈滚道中心到外圈中心的距离为：

$r_{r1}=D_m/2+\left(f_1-0.5\right)\cdot D_w$ (1)

内圈滚道中心到内圈中心的距离为：

$r_{r2}=D_m/2-\left(f_2-0.5\right)\cdot D_w$ (2)

式中，f₁为外圈滚道的曲率系数，f₂为内圈滚道的曲率系数。

滚珠和滚道之间的接触变形量方向沿着矢量O_rjO_bj的方向，计算表达式为：

${\delta }_j=\left(f_{1\left(2\right)}-0.5\right)\cdot D_w-\left|O_{rj}{O}_{bj}\right|$ (3)

O_rjO_bj表示滚珠中心和套圈滚道中心之间的相对位置，可由滚珠和套圈的位移来表示：

$\{\begin{array}{l}{O}_{rj}{O}_{bj}\left(外圈\right)=O{O}_{bj}-O{O}_{rj}=\left\{\begin{array}{c}{x}_j\\ r_j\sin {\theta }_j\\ r_j\cos {\theta }_j\end{array}\right\}-\left\{\begin{array}{c}0\\ r_{r_j}\sin {\phi }_j\\ r_{r_j}\cos {\phi }_j\end{array}\right\}\\ O_{rj}{O}_{bj}\left(�圈\right)=O{O}_{bj}-O{O}_2-O_2{O}_{rj}=\left\{\begin{array}{c}{x}_j\\ r_j\sin {\theta }_j\\ r_j\cos {\theta }_j\end{array}\right\}-\left\{\begin{array}{c}{x}_2\\ y_2\\ z_2\end{array}\right\}-T_{ir}\left\{\begin{array}{c}0\\ r_{r_j}\sin {\phi }_j\\ r_{r_j}\cos {\phi }_j\end{array}\right\}\end{array}$ (4)

式中，x_j为滚珠沿X方向的轴向位移；r_j为滚珠沿OO_bj方向的径向位移；${\phi }_j$ 为第j个滚珠的方位角；T_ir为内圈坐标到整体坐标的转换矩阵：

$T_{ir}=\left[\begin{array}{ccc}\cos {\theta }_y\cos {\theta }_z& \cos {\theta }_y\sin {\theta }_z& \sin {\theta }_y\\ -\sin {\theta }_z& \cos {\theta }_z& 0\\ -\sin {\theta }_y\cos {\theta }_z& -\sin {\theta }_y\sin {\theta }_z& \cos {\theta }_y\end{array}\right]$ (5)

根据滚珠和套圈之间的接触变形量δ_j，结合赫兹接触理论[9]，求解接触点上的接触椭圆长半轴a、接触椭圆短半轴b、接触椭圆的椭圆率k、接触载荷Q_j以及最大接触应力q_max。

滚珠与套圈之间存在起到润滑作用的油膜[10]，根据Hamrock-Dowson油膜厚度公式计算得到油膜厚度：

$h_{\text{min}}=\frac{0.2207{\alpha }_0^{0.49}{D}_w^{0.466}{\left(1-\gamma \right)}^{1.15}{\left(1+\gamma \right)}^{0.68}\left(1-e^{-0.68k_2}\right)}{{\left({\eta }_0\cdot ni\cdot D_m\right)}^{-0.68}{E^{\prime }}^{0.117}{Q}_{\text{max}}{\left(2\right)}^{0.073}}$ (6)

ni为电机运行转速；γ为无量纲几何参数，$\gamma =D_w/D_m$ ；α₀为粘压系数；η₀为动力黏度；Q_max为表示滚珠所受最大载荷；$E^{\prime }$ 为轴承材料的等效弹性模量；D_m为轴承节圆直径。

将滚珠所受拖动力进行离散化处理：

$F_{br}={\displaystyle \sum _{k=-n}^n\left[\frac{3}{4}Q\left(1-{\overline{x}}_k^2\right)\Delta \overline{x}\cdot \left(\frac{{\kappa }_0}{p_k{h}_k}{\kappa }_k^{*}\Delta u_k\right)\right]}$ (7)

式中，n为单侧切片数量；Q为滚珠与滚道之间的接触载荷；${\overline{x}}_k^{}$ 为无量纲切片位置，${\overline{x}}_k^{}=a\cdot \frac{k}{n}$ ，a为

接触椭圆长半轴；${\kappa }_0$ 、${\kappa }_k^{*}$ 为润滑油流变参数[11]；p_k为单个切片上的赫兹接触压力；h_k为切片所处位置的油膜厚度，为了简化计算，此处取h_k = h_min；$\Delta u_k$ 为滚珠和滚道的相对滑动速度：

$\Delta u_k=v_b+\left(T_{ib}\cdot {\omega }_b\right)\times r_{cj}-v_2-{\omega }_2\times R_{cj}$ (8)

v_b为滚珠的质心位移线速度；w_b为滚珠的质心位移角速度；r_cj为滚珠中心到接触点的位置；v₂和w₂分别为内圈滚道的位移线速度和角速度，由于假设外圈固定，在计算外圈滚道接触的情况时取0；T_ib为滚珠和整体坐标系的转换矩阵：

$T_{ib}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\theta }_j& \sin {\theta }_j\\ 0& -\sin {\theta }_j& \cos {\theta }_j\end{array}\right]$ (9)

R_cj为内圈滚道中心到接触点的位置，通过滚珠、内圈滚道和接触点之间的相对位置表示为：

$R_{cj}=r_{cj}+O{O}_{bj}-O_2{O}_{rj}$ (10)

其中，r_cj通过对切片坐标系下的滚珠和接触点之间的距离进行坐标系转换，表示为：

$r_{cj}=T_{bcn}\cdot r_{cj}^{cn}=T_{bcn}\cdot \left[\begin{array}{c}\overline{x}a\\ \overline{y}b\\ \sqrt{{\overline{R}}^2-{\left(\overline{x}\cdot a\right)}^2}-\sqrt{{\overline{R}}^2-a^2}+\sqrt{D_w^2/4-a^2}\end{array}\right]$ (11)

$\overline{x}$ 和$\overline{y}$ 为无量纲参数，用于描述单个切片在整体接触椭圆上的位置，取值范围为[−1, 1]；$\overline{R}$ 为滚珠和套圈滚道的等效半径(m)，$\overline{R}=\frac{2f_{1\left(2\right)}\cdot D_w}{2f_{1\left(2\right)}+1}$ ，$f_{1\left(2\right)}$ 为外圈(内圈)滚道的曲率系数；T_bcn为滚珠和接触点坐标系的转换矩阵：

$T_{bcn}=\left[\begin{array}{ccc}\cos {\gamma }_1& 0& -\sin {\gamma }_1\\ \sin {\gamma }_1\cdot \sin {\gamma }_2& \cos {\gamma }_2& \cos {\gamma }_1\cdot \sin {\gamma }_2\\ \sin {\gamma }_1\cdot \cos {\gamma }_2& -\sin {\gamma }_2& \cos {\gamma }_1\cdot \cos {\gamma }_2\end{array}\right]$ (12)

式中，γ₁和γ₂分别为滚珠和滚道在X-Z平面和Y-Z平面上的接触角，${\gamma }_1=\arctan \frac{e_{brx}^b}{e_{brz}^b}$ ， ${\gamma }_2=\arctan \frac{-e_{bry}^b}{\sqrt{{\left(e_{brx}^b\right)}^2+{\left(e_{brz}^b\right)}^2}}$ ；$e_{br}^b$ 表示滚珠中心和套圈滚道中心的相对位置。

滚珠所受拖动力矩通过拖动力和转动半径叉乘得到：

$M_{br}=F_{br}\times r_{cj}$ (13)

2.3. 保持架与滚珠相互作用模型

滚珠在兜孔内运动时会与兜孔的不同部位产生接触，为了区分接触区域，将椭圆形的兜孔简化为矩形和两个半圆，如图3所示。

图中，S_p为兜孔矩形部分的长度；R_p为兜孔圆形部分的半径；O_pj为兜孔中心、O_prj为兜孔圆形部分中心、O_bj为滚珠中心，后缀j为滚珠编号；${\delta }_{bh}$ 是滚珠和保持架兜孔之间的间隙，根据接触部位区分为：

$\{\begin{array}{ll}{\delta }_{bh\left(�形部分\right)}=\left(R_p-\frac{D_w}{2}\right)-\sqrt{{\left(O_{prj}{O}_{bj}\right)}_x^2+{\left(O_{prj}{O}_{bj}\right)}_z^2},\hfill & {\left(O_{pj}{O}_{bj}\right)}_x-\frac{S_p}{2}\ge 0\hfill \\ {\delta }_{bh\left(矩形部分\right)}=\left(\frac{S_p}{2}-\frac{D_w}{2}\right)-\left|{\left(O_{pj}{O}_{bj}\right)}_z\right|,\hfill & {\left(O_{pj}{O}_{bj}\right)}_x-\frac{S_p}{2}<0\hfill \end{array}$ (14)

Figure 3. Schematic diagram of pocket contact parts. (a) Rectangular part contact, (b) Circular part contact

图3. 兜孔接触部位示意图。(a) 矩形部分接触，(b) 圆形部分接触

通过对比间隙${\delta }_{bh}$ 的初始值和几何计算值，当初始值大于计算值时，滚珠和保持架接触并产生接触变形，接触变形量${\delta }_j={\delta }_{bh}$ 。

O_prjO_bj为保持架兜孔圆形部分中心指向滚珠中心的矢量，计算公式如下：

$O_{prj}{O}_{bj}=O_{pj}{O}_{bj}-T_{ih}{T}_{hp}{\left[\begin{array}{ccc}0& \pm \frac{S_p}{2}& 0\end{array}\right]}^{\text{T}}$ (15)

式(15)中的正负值与矢量O_pjO_bj的方向相关。

O_pjO_bj为保持架兜孔中心和滚珠中心之间的相对位置：

$O_{pj}{O}_{bj}=\left\{\begin{array}{c}{x}_j\\ r_j\sin {\theta }_j\\ r_j\cos {\theta }_j\end{array}\right\}-\left\{\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right\}-T_{ih}\left\{\begin{array}{c}0\\ \frac{D_c}{2}\sin {\theta }_{pj}\\ \frac{D_c}{2}\cos {\theta }_{pj}\end{array}\right\}$ (16)

式中，x_c、y_c、z_c为保持架在三个方向上的质心位移；D_c为保持架的直径；θ_pj为保持架兜孔的方位角；T_ih为保持架坐标系与整体坐标系的转换矩阵，通过保持架的偏转角θ_hy和θ_hz表示为：

$T_{ih}=\left[\begin{array}{ccc}\cos {\theta }_{hy}\cos {\theta }_{hz}& \cos {\theta }_{hy}\sin {\theta }_{hz}& \sin {\theta }_{hy}\\ -\sin {\theta }_{hz}& \cos {\theta }_{hz}& 0\\ -\sin {\theta }_{hy}\cos {\theta }_{hz}& -\sin {\theta }_{hy}\sin {\theta }_{hz}& \cos {\theta }_{hy}\end{array}\right]$ (17)

T_hp为保持架兜孔坐标系与保持架坐标系的转换矩阵：

$T_{hp}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\theta }_{pj}& \sin {\theta }_{pj}\\ 0& -\sin {\theta }_{pj}& \cos {\theta }_{pj}\end{array}\right]$ (18)

滚珠和保持架兜孔之间存在一层微米级的润滑油膜，二者之间的间隙大小直接影响其接触状态，判断依据如下：

$\{\begin{array}{l}{\delta }_{bh}<\sqrt{{\sigma }_b^2+{\sigma }_h^2},赫兹接触状态\\ {\delta }_{bh}>\sqrt{{\sigma }_b^2+{\sigma }_h^2},混合润滑状态\end{array}$ (19)

式中，${\sigma }_b$ 为滚珠的表面粗糙度，本文取0.02 μm；${\sigma }_h$ 为保持架的表面粗糙度，本文取0.4 μm。

当滚珠和保持架处于赫兹接触状态时，二者之间的拖动力为：

$F_{bp}=-F_e\cdot e_{pb}-\frac{F_{\mu e}}{\left|u_{pb}\right|}\cdot u_{pb}$ (20)

式中，F_e为滚珠和保持架兜孔之间的接触载荷，通过赫兹接触理论公式对接触变形量δ_j进行计算得到；$F_{\mu e}$ 为滚珠和保持架兜孔之间的摩擦力，$F_{\mu e}=F_e\cdot f_{\mu e}$ ，$f_{\mu e}$ 为摩擦系数；e_pb为滚珠中心和兜孔接触区域中心的相对位置单位矢量，计算公式如下：

$\{\begin{array}{l}{e}_{pb}{}_{\left(�形部分\right)}=\frac{{\left[\begin{array}{ccc}{O}_{prj}{O}_{bj}{}_x& O_{prj}{O}_{bj}{}_y& 0\end{array}\right]}^{\text{T}}}{\sqrt{{\left(O_{prj}{O}_{bjx}\right)}^2+{\left(O_{prj}{O}_{bjy}\right)}^2}}\\ e_{pb}{}_{\left(矩形部分\right)}=\pm 1\end{array}$ (21)

与矩形部分接触时，滚珠中心与接触点始终处于接触区域的法线上，其正负值取决于滚珠在接触点的受力方向。

u_pb为滚珠和保持架兜孔接触点处的相对运动速度，根据接触区域的不同区分为：

$\{\begin{array}{l}{u}_{pb}{}_{\left(�形部分\right)}=T_{ib}^{\text{T}}\cdot \left(v_b-v_p\right)+w_b\times \frac{D_w}{2\cdot e_{pb}}\\ u_{pb}{}_{\left(矩形部分\right)}=T_{ib}^{\text{T}}\cdot \left(v_b-v_p\right)+w_b\cdot T_{ib}^{\text{T}}\cdot T_{ip}\cdot \frac{D_w}{2}\cdot {\left[\begin{array}{ccc}\mathrm{sgn}\left(O_{pj}{O}_{bj}{}_x\right)& 0& 0\end{array}\right]}^{\text{T}}\end{array}$ (22)

v_b为滚珠质心的线速度(m/s)；v_p为保持架兜孔中心线速度(m/s)；w_b为滚珠的自转角速度(rad/s)，使用龙格库塔方法对平衡方程求解得到。

当滚珠和保持架的接触区域处于混合润滑状态时，滚珠表面受到润滑油膜压力造成的法向压力以及润滑油粘性剪切力，阻碍了滚珠和保持架的运动，在轴承内部产生拖动力：

$F_{bp}^{}=-\left(F_{ex}+F_{ez}\right)\cdot e_{pb}^{}-\frac{F_{\mu x}}{\sqrt{{\left(u_{pbx}\right)}^2+{\left(u_{pby}\right)}^2}}\cdot \left[\begin{array}{ccc}{u}_{pbx}& u_{pby}& 0\end{array}\right]-F_{\mu z}\cdot \frac{u_{pbz}}{\left|u_{pbz}\right|}$ (23)

式中，F_ex为滚珠在保持架兜孔沿滚珠自身的X-Y平面滚动时，由油膜压力造成的法向压力；F_ez为滚珠在保持架兜孔内绕Z轴自旋时，受油膜压力影响产生的法向压力；F_μx为滚珠在保持架兜孔内沿着X-Y平面滚动/滑动时，因润滑油粘性而产生的粘性剪切力；F_μz为滚珠在保持架兜孔内沿着Z轴自旋时产生的

润滑油粘性剪切力。对参数进行无量纲化处理：油膜在x方向上的无量纲分布$\overline{\epsilon }$ ，$\overline{\epsilon }=\frac{2x}{l}$ ；油膜在y方向上的无量纲分布$\overline{\zeta }$ ，$\overline{\zeta }=\frac{2y}{l}\sqrt{\frac{R_x}{R_y}}$ ；无量纲几何参数$\overline{L}$ ，$\overline{L}=\frac{l}{\sqrt{2R_x{h}_m}}$ 。l表示保持架的厚度(m)；R_x，R_y

表示滚珠和保持架兜孔接触表面在x方向和y方向上的等效半径；h_m表示滚珠和保持架之间的油膜厚度，并将油膜分布简化为一个半球面。

分别对雷诺方程压力计算公式和润滑油剪切力按照油膜压力的分布范围进行积分，得到法向压力和剪切力的分布，计算公式如下：

(24)

式(24)中，${\overline{\epsilon }}_0$ 对应油膜在x方向上压力为0的位置，${\overline{\zeta }}_0$ 对应油膜在y方向上压力为0的位置，油膜压力集中在${\overline{\epsilon }}_0\times {\overline{\zeta }}_0$ 的矩形区域中。由于保持架兜孔在y方向上为椭圆形的封闭曲线，假设其值在Y轴正方向与兜孔的交点处，本文取1。x方向上压力为0的位置通过以雷诺方程压力公式$\overline{P}\left(\overline{\zeta },{\overline{\epsilon }}_0\right)=0$ 作为边界条件进行确定，即：

$\overline{P}\left(\overline{\zeta },{\overline{\epsilon }}_0\right)={\displaystyle {\int }_{-1}^{\overline{\epsilon }}\frac{x^2-{\overline{\epsilon }}_0^2}{{\left(\frac{1}{{\overline{L}}^2}+x^2+{\overline{\zeta }}^2\right)}^3}\text{d}x}=0$(25)

最后，通过滚珠中心到接触点的向量与滚珠和保持架兜孔之间的拖动力叉乘计算，得到二者之间的拖动力矩：

$\{\begin{array}{l}{M}_{bp}^{}{}_{\left(�形部分\right)}=\left(\frac{D_w}{2\cdot e_{pb}^{}{}_{\left(�形部分\right)}}\right)\times F_{bp}^{}{}_{\left(�形部分\right)}\\ M_{bp}^{}{}_{\left(矩形部分\right)}=\left(T_{ib}^{\text{T}}\cdot T_{ip}\cdot \frac{D_w}{2}\cdot \left[\begin{array}{c}\mathrm{sgn}\left(O_{pj}O{}_{bjx}\right)\\ 0\\ 0\end{array}\right]\right)\times F_{bp}^{}{}_{\left(矩形部分\right)}\end{array}$ (26)

2.4. 保持架与套圈相互作用模型

在电机中，轴承处于高速且相对稳定的工作状态，因此本文采用外圈引导的方式，套圈与保持架之间的相互作用如图所示[12]。

图4是采用外引导方式下，保持架的受力情况示意，w_c为保持架的转速；${\phi }_d$ 为保持架的偏心角，${\phi }_d=\arctan \frac{y_c}{z_c}$ ；F_M为套圈和保持架之间的油膜压力对保持架产生的径向压力；$F_{\varphi }$ 为油膜对保持架的切向拖动力；$W_{\rho }$ 为油膜对保持架的阻碍力矩。计算公式如下：

Figure 4. Schematic diagram of the interaction between the cage and the rings

图4. 保持架与套圈相互作用示意图

$\{\begin{array}{l}{F}_M=-\frac{\eta R_g{B}_g^3}{C_g^2}\cdot \frac{{\epsilon }^2}{{\left(1-{\epsilon }^2\right)}^2}{\omega }_c\\ F_{\varphi }=-\frac{\pi \eta R_g{B}_g^3}{4C_g^2}\cdot \frac{\epsilon }{{\left(1-{\epsilon }^2\right)}^{3/2}}{\omega }_c\\ W_{\rho }=-\frac{2\pi \eta R_g^3{B}_g}{C_g}\cdot \frac{1}{\sqrt{1-{\epsilon }^2}}{\omega }_c\end{array}$ (27)

式中，$\eta$ 为润滑油粘度；R_g为引导面半径；B_g为引导面宽度；C_g为轴承引导间隙；$\epsilon$ 为保持架的偏心率，$\epsilon =\frac{r_c}{C_g}$ ；r_c为保持架质心与外圈中心之间的距离，$r_c=\sqrt{y_c^2+z_c^2}$ 。

通过转换矩阵，将保持架坐标系下的受力转换到轴承整体坐标系下：

$\left\{\begin{array}{c}{F}_{ly}\\ F_{lz}\\ M_l\end{array}\right\}=\left[\begin{array}{ccc}\cos {\phi }_d& -\sin {\phi }_d& 0\\ \sin {\phi }_d& \cos {\phi }_d& 0\\ 0& 0& 1\end{array}\right]\left\{\begin{array}{c}{F}_M\\ F_{\varphi }\\ W_{\rho }\end{array}\right\}$ (28)

2.5. 轴承平衡方程

轴承内部除了各部件之间的作用力，还需要考虑润滑油对部件产生的阻滞力和阻滞力矩。根据Schlichting理论[13]，滚珠所受润滑油阻力为：

$F_D=\frac{1}{2}{C}_D{\rho }_{oil}{\left({\omega }_j{D}_m/2\right)}^{2}A$ (29)

式中，${\rho }_{oil}$ 为润滑油密度；w_j为滚珠的公转速度；A为滚珠在润滑油中的投影面积，$A=\frac{1}{4}\pi D_w^2-l\cdot D_w$ ；C_D为阻力系数。

滚珠所受润滑油阻力矩为：

$M_e=\frac{1}{2}{\rho }_{oil}{\omega }_j{\left(D_w/2\right)}^5{C}_n$ (30)

式中，Re为雷诺系数，C_n为流体系数。

滚珠和轴承内圈的瞬时平衡位置由代数方程形式的平衡方程组确定：

$\{\begin{array}{l}-Q_{1j}\sin {\alpha }_{1j}+Q_{2j}\sin {\alpha }_{2j}+F_{bpjx}=0\hfill \\ -Q_{1j}\cos {\alpha }_{1j}+Q_{2j}\cos {\alpha }_{2j}+F_{bpjz}+m_b{\omega }_j^2{r}_j=0\hfill \end{array}$ (31)

对于含N个滚珠的轴承内圈，其平衡方程为：

$\{\begin{array}{l}{F}_X-{\displaystyle \sum _{j=1}^N\left(Q_{2j}\sin {\alpha }_{2j}-F_{br}\left(1\right)\right)}=0\hfill \\ F_Y-{\displaystyle \sum _{j=1}^N\left[\left(Q_{2j}\cos {\alpha }_{2j}-F_{br}\left(3\right)\right)\sin {\theta }_j-F_{br}\left(2\right)\cos {\theta }_j\right]}=0\hfill \\ F_Z-{\displaystyle \sum _{j=1}^N\left[\left(Q_{2j}\cos {\alpha }_{2j}-F_{br}\left(3\right)\right)\cos {\theta }_j-F_{br}\left(2\right)\sin {\theta }_j\right]}=0\hfill \\ M_Y-{\displaystyle \sum _{j=1}^N\left[r_2\left(Q_{2j}\sin {\alpha }_{2j}+F_{br}\left(1\right)\right)-f_2{D}_w{F}_{br}\left(1\right)\right]}\cos {\theta }_j=0\hfill \\ M_Z-{\displaystyle \sum _{j=1}^N\left[r_2\left(Q_{2j}\sin {\alpha }_{2j}+F_{br}\left(1\right)\right)-f_2{D}_w{F}_{br}\left(1\right)\right]}\sin {\theta }_j=0\hfill \end{array}$ (32)

式中，F_X、F_Y、F_Z、M_Y和M_Z为轴承所受外部载荷；${\alpha }_{1j}$ 和${\alpha }_{2j}$ 分别为滚珠和外圈以及内圈的接触角[14]。

滚珠的运动遵循牛顿第二定律，通过微分运动方程组进行求解：

$\{\begin{array}{l}{F}_{brj}\left(2\right)-F_D+F_{bpj}\left(2\right)=m_b{r}_{bj}{\ddot{\theta }}_j\hfill \\ M_{brj}\left(1\right)-M_{ej}\left(1\right)+M_{bpj}\left(1\right)=I_b{\dot{\omega }}_j\left(1\right)\hfill \\ M_{brj}\left(2\right)-M_{ej}\left(2\right)+M_{bpj}\left(2\right)=I_b{\dot{\omega }}_j\left(2\right)\hfill \\ M_{brj}\left(3\right)-M_{ej}\left(3\right)+M_{bpj}\left(3\right)=I_b{\dot{\omega }}_j\left(3\right)\hfill \end{array}$ (33)

保持架的运动微分方程为：

$\{\begin{array}{l}{\displaystyle \sum _{j=1}^N\left(-F_{bpj}\left(1\right)\right)}=m_c{\ddot{x}}_c\hfill \\ F_{ly}+{\displaystyle \sum _{j=1}^N\left(-F_{bpj}\left(2\right)\cdot \cos {\theta }_j-F_{bpj}\left(3\right)\cdot \sin {\theta }_j\right)}=m_c{\ddot{y}}_c\hfill \\ F_{lz}+{\displaystyle \sum _{j=1}^N\left(-F_{bpj}\left(2\right)\cdot \sin {\theta }_j-F_{bpj}\left(3\right)\cdot \cos {\theta }_j\right)}=m_c{\ddot{z}}_c\hfill \\ {\displaystyle \sum _{j=1}^N\left(-F_{bpj}\left(2\right)\cdot \frac{D_c}{2}\right)}+M_l=I_c{\ddot{\theta }}_c\hfill \end{array}$ (34)

式中，m_b为滚珠质量；m_c为保持架质量；D_c为保持架直径；I_b为滚珠的转动惯量；I_c为保持架的转动惯量[15]。

3. 动力学模型求解与验证

3.1. 动力学模型求解流程

轴承动力学模型求解流程如图5所示：

根据轴承型号和所处工况，设置结构参数、材料参数、工况参数，并利用经验公式求解轴承动力学模型的启动值，加快模型迭代的收敛速度。基于启动值，使用牛顿–拉夫逊法对平衡方程(32)和(33)进行求解，得到每个滚珠的轴向位移x_j、径向位移r_j，轴承内圈的位移x₂、y₂、z₂以及角位移θ_y、θ_z；使用四阶龙格–库塔法对平衡方程(34)和(35)进行求解，得到轴承的动力学特性参数。

Figure 5. Solution process of the bearing dynamic model

图5. 轴承动力学模型求解流程

3.2. 动力学模型验证

通过与Yuqing Liu [16]所建立的柔性保持架轴承模型进行对比，验证模型的可行性。参考文献的轴承结构参数和工况参数如表1所示。

Table 1. Reference parameter table

表1. 参考文献参数表

参数符号	参数名称	参数值
N	滚珠数量	7
Dm	节圆直径(mm)	36
Rb	滚珠半径(mm)	4.76
R1	外圈滚道半径(mm)	22.77
R2	内圈滚道半径(mm)	13.24
μ0	干摩擦系数	0.02
FX	轴承所受载荷(N)	500
aw2	轴承内圈角加速度(rad/s^2)	500

为了与参考模型的输入工况匹配，需要将固定的轴承内圈转速修改为随加速度变化的变转速，即$w_2=t\times a_{w2}$ ，t为转速达到预期的时间。

采用相同的结构参数与工况参数，与参考文献所建立模型仿真计算得到的滚珠与内圈接触载荷值进行对比，结果如图6所示。

图6中，黑色实线为参考文献值，红色虚线为本文仿真值。在内圈加速初期，滚珠与内圈之间碰撞频繁，因此降低第0 s到第0.2 s之间的采样频率，使得二者在加速初期的偏差较大，为13.3%。滚珠与内圈碰撞频率相对稳定之后的接触载荷平均偏差为3.7%，认为所建立的轴承动力学模型可以有效仿真轴承的动力学特性。

Figure 6. Contact load between the ball and the inner ring

图6. 滚珠与内圈的接触载荷

4. 轴承动力学特性分析与寿命计算

4.1. 多工况仿真分析

根据电机轴承的运行数据，采用K-Means算法提取典型工况，将工况划分为工况一(1730 RPM, 938 N)：低转速高载荷；工况二(3703 RPM, 903 N)：高转速高载荷；工况三(1785 RPM, 266 N)：低转速低载荷；工况四(3998 RPM, 206 N)：高转速低载荷四种工况，如表2所示[17]。

Table 2. Typical working conditions of motor bearings

表2. 电机轴承典型工况

工况	转速(RPM)	载荷(N)
工况一	1730	938
工况二	3703	903
工况三	1785	266
工况四	3998	206

基于典型工况，对轴承进行仿真，得到轴承的动力学特性。

由图7可知，随着轴承转速的降低和所受载荷的提高，轴承保持架的打滑率降低到稳定状态所需时长随之降低，说明保持架打滑率稳定速率与轴承转速反相关，与轴承所受载荷正相关。高载荷增强了滚珠的同步驱动能力，使得保持架能够更快地被滚珠加速到稳定的公转速度，从而更快地消除初始的打滑现象，使得打滑率更快趋于稳定。

Figure 7. Cage slip rate under typical working conditions

图7. 典型工况下的轴承保持架打滑率

4.2. 疲劳寿命计算

根据L-P轴承疲劳寿命理论，计算轴承整体的疲劳寿命：

$L_{10}^{-\frac{10}{9}}=\frac{1}{Z}\left({\displaystyle \sum _{j=1}^Z{L}_{1j}^{-\frac{10}{9}}}+{\displaystyle \sum _{j=1}^Z{L}_{2j}^{-\frac{10}{9}}}\right)+{{\displaystyle \sum _{j=1}^Z\left(\frac{L_{b1j}}{N_j}\right)}}^{-\frac{10}{9}}+{{\displaystyle \sum _{j=1}^Z\left(\frac{L_{b2j}}{N_j}\right)}}^{-\frac{10}{9}}$ (35)

式中，L_1(2)j为轴承套圈的等效疲劳寿命；L_b1(2)j为滚珠的等效疲劳寿命；N_j为轴承内圈旋转一周时滚珠的转数，具体的计算方法参考文献[18]。基于典型工况，计算得到不同工况下的轴承疲劳寿命，如图8所示。

Figure 8. Bearing fatigue life under typical working conditions

图8. 典型工况下的轴承疲劳寿命

由图可知，在高载荷的工况下，如工况一(1730 RPM, 938 N)和工况二(3703 RPM, 903 N)轴承的寿命较短；在低载荷的工况下，如工况三(1785 RPM, 266 N)和工况四(3998 RPM, 206 N)轴承的寿命显著增长；同时，随着转速提高，轴承的寿命相应减少，因此判断轴承的寿命与轴承转速和所受载荷呈反相关关系。

5. 结论

1) 与文献模型对比验证表明，滚珠与内圈接触载荷稳定后的平均偏差为3.7%，模型可有效表征轴承动态运行特性。

2) 通过对典型工况下的轴承动力学特性分析，揭示了轴承各部件之间的作用规律。保持架打滑率的稳定速率与转速呈反相关、与载荷呈正相关，低转速高载荷工况下的打滑率能更快地趋于稳定；滚珠与保持架之间的碰撞力在启动阶段震荡明显，且与载荷呈正相关，高载荷低转速工况下的碰撞力能更早地进入稳定状态。

3) 基于L-P理论的疲劳寿命计算表明，轴承寿命与内圈转速和所受载荷呈反相关关系，高载荷工况(工况一、工况二)下的轴承寿命较短，低载荷工况(工况三、工况四)下的轴承寿命显著提升。

建立的动力学模型及研究方法为深沟球轴承的动态特性分析、疲劳寿命预测提供了有效理论支撑，可为轴承结构优化设计和提高轴承耐久性提供参考依据。