滑动轴承轴心轨迹研究

Study of the Centre Motion Trajectory of Journal Bearings

摘要: 滑动轴承轴心轨迹可以反映出轴承的运行状态。通过对轴心轨迹的研究可以分析轴承设计于使用工况的合理性。本文基于纳维–斯托克斯方程(N-S方程)和流量连续性方程推导了滑动轴承油膜润滑雷诺方程，使用有限体积法求解雷诺方程，通过逐次超松弛迭代法(SOR)求解线性方程组，获得油膜压力分布，构建了轴心轨迹计算模型。基于该模型分析了主轴转速、配合间隙和主轴质量对轴心轨迹的影响。

Abstract: The centre motion trajectory of the journal bearing can reflect the operating state of the bearing. Through the study of the spindle trajectory, the rationality of the bearing design can be analyzed. In this paper, the Reynolds equation for oil film lubrication of journal bearing is derived based on the Navier-Stokes equation (N-S equation) and the flow continuity equation, and the Reynolds equation is solved by the finite volume method, and the linear equation system is solved by the successive over relaxation iterative method (SOR) to obtain the oil film pressure distribution. The calculation model of centre motion trajectory was constructed. Based on this model, the effects of spindle speed, mating clearance and spindle quality on the spindle trajectory were analyzed.

1. 引言

轴承是旋转机械中必不可少的零件，其能够在保证旋转精度的前提下降低旋转机械的旋转摩擦力。目前轴承可分为滑动轴承和滚动轴，滑动轴承由于其承载能力强、噪声小、结构简单和体积小等特点被广泛应用[1]，尤其是结构紧凑的旋转机械中[2]。

在轴心轨迹方面，不同的轴心轨迹可以反映出不同的故障信息，它可以直观反映出系统的工作状态，对轴承故障诊断具有重要意义。在轻微的不对中的工况下轴心轨迹则呈椭圆形；在不对中方向上施加一个中等负载后，轴心轨迹则变为香蕉形；严重不对中的情况使转子轴心轨迹呈现外“8”字形，油膜涡动引起的轴心轨迹为内“8”字形[3]。轴心轨迹不仅可以用于轴承故障诊断，还可用于异形截面加工，在外界人为干预的情况下，轴心轨迹可以按照预期轨迹运行。沈庆崇等人提出电磁铁–流体润滑滑动轴承系统，通过主动控制电磁铁输入电流来控制滑动轴承运动轨迹[4]。

在计算流体力学方面，由于雷诺方程是一个非线性的偏微分方程，很难获得方程的解析解，因此需要通过数值方法获取数值解。目前针对雷诺方程的数值计算方法主要包括有限差分法、有限体积法和有限元法[5]。Qiu等人分别采用有限差分法和无限小扰动法计算轴承油膜的非线性动系数，结果表明使用这两种方法计算结果差距很小[6]。敏政等人使用有限差分法求解了油膜压力分布的数值解[7]。张伟等人分析了有限体积法和有限元法之间的优缺点，有限元法对椭圆形问题有更好的适应性，但计算量偏大。有限体积法相较于有限元法，其计算过程占用内存更小[8]。

液压领域中，旋转机械主要包括液压泵和液压马达。在设计过程中考虑到更紧凑的结构，轴向柱塞泵部分位置并未使用滚动轴承而是使用滑动轴承。为了防止主轴出现较大范围跳动进而影响柱塞泵寿命，本文基于某型号轴向柱塞泵后盖中滑动轴承开展轴心轨迹分析。

2. 滑动轴承油膜方程建立

两个相互倾斜的平板之间会形成一个楔形间隙。当其中一个平板沿间隙收敛的方向移动时，填充在两平板间的粘性流体会被带入该楔形间隙。随着流体的不断流入，压力逐渐增大，最终形成一个稳定的动压油膜，从而避免两个平板直接接触，实现了流体润滑。

2.1. 油膜压力方程

纳维–斯托克斯方程(N-S方程)是一组用于描述粘性流体流动的偏微分方程组。在直角坐标系中其表达形式如式(1)所示。

根据N-S方程，针对滑动轴承中流体流动情况做出以下假设：

(1) 不考虑流体的压缩性，流体密度为常数。

(2) 轴承与轴之间的流体流动为层流。

(3) 惯性力远小于粘性力。

(4) 流体厚度极小，油膜压力在厚度方向上的变化忽略不计，即$p=p\left(x,y\right)$ 。

(5) 流体流速在厚度方向上的梯度远大于另外两个方向上的梯度，即$\frac{\partial u,v}{\partial z}\gg \frac{\partial u,v}{\partial x}$ ，$\frac{\partial u,v}{\partial z}\gg \frac{\partial u,v}{\partial y}$ 。

$\{\begin{array}{l}\rho \left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}\right)=\rho f_x-\frac{\partial p}{\partial x}+\mu \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)\\ \rho \left(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}\right)=\rho f_y-\frac{\partial p}{\partial y}+\mu \left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}+\frac{{\partial }^{2}v}{\partial z^2}\right)\\ \rho \left(\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}\right)=\rho f_z-\frac{\partial p}{\partial z}+\mu \left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}+\frac{{\partial }^{2}w}{\partial z^2}\right)\end{array}$ (1)

引入稳态下连续性方程：

$\frac{\partial v_{\text{x}}}{\partial x}+\frac{\partial v_y}{\partial y}+\frac{\partial v_{\text{z}}}{\partial z}=0$ (2)

得到滑动轴承油膜雷诺方程：

$\frac{\partial }{\partial x}\left(\frac{h^3}{\mu }\frac{\partial p}{\partial x}\right)+\frac{\partial }{\partial y}\left(\frac{h^3}{\mu }\frac{\partial p}{\partial y}\right)=6\omega r\frac{\partial h}{\partial x}+12\frac{\partial h}{\partial t}$ (3)

式中，$\mu$ ——流体动力粘度；

$\omega$ ——主轴转动角速度；

$r$ ——轴半径。

2.2. 油膜厚度方程

油膜压力与油膜的厚度密切相关，在进行油膜压力的计算之前，需要先确定油膜的厚度分布情况(图1)。

Figure 1. Schematic diagram of oil film thickness

图1. 油膜厚度示意图

在不考虑多场耦合的前提下，油膜厚度分布情况仅与轴颈的姿态变化有关。假设轴颈仅在$xoz$ 平面中平动，因此通过设定$e_x$ 和$e_z$ 来表述轴颈在滑动轴承中的姿态。其中，$e_x$ 和$e_z$ 分别为轴颈在$x$ 和$z$ 方向上的偏移值。整理后可得到简化后的油膜厚度方程。

$h=\sqrt{{\left(R\cos {\theta }_1-e_x\right)}^2+{\left(R\cos {\theta }_1-e_z\right)}^2}-r$ (4)

式中，$R$ ——滑动轴承内圈半径。

3. 雷诺方程求解

3.1. 有限体积法求解

本文选用有限体积法对雷诺方程进行求解。具体流程如下：首先将流体域展开成矩形，并基于上述简化条件将三维流动问题转化为二维流动问题。其次使用网格交错法分别划分流体域和油膜厚度场网格(如图2所示)，并计算每个点位的油膜厚度值。最终根据给定的初始条件，使用SOR迭代法对每一个节点处的压力值进行迭代。

采用笛卡尔直角坐标系，将滑动轴承流体域展开，流体域由三维简化为二维。其中，$x$ 方向为圆周方向，$y$ 方向为轴向方向。将流体域在$x$ 和$y$ 方向上分别划分成$pnx$ 和$pny$ 份，则在$x$ 和$y$ 方向上分别获得$pnx+1$ 和$pny+1$ 个节点(正方形节点)，所有正方形节点用于存储该点位上的压力值。基于网格交错法，在任意两个压力节点之间插入一个圆形节点，用于存储该点位上的油膜厚度值。压力节点既存储压力值也存储油膜厚度值。因此在压力场中网格节点数为$\left(pnx+1\right)\left(pny+1\right)$ ，油膜厚度场中网格节点数为$\left(2pnx+1\right)\left(2pny+1\right)$ 。

Figure 2. Schematic diagram of meshing

图2. 网格划分示意图

在展开后的流体域选取一控制体用于表述有限体积法对雷诺方程的离散过程，图中阴影区域为一个控制体，对式(3)在控制体内积分得到：

${\displaystyle \underset{s}{\overset{n}{\int }}{\displaystyle \underset{w}{\overset{e}{\int }}\frac{\partial }{\partial x}}}\left(\frac{h^3}{\mu }\frac{\partial p}{\partial x}\right)\text{d}x\text{d}y+{\displaystyle \underset{s}{\overset{n}{\int }}{\displaystyle \underset{w}{\overset{e}{\int }}\frac{\partial }{\partial y}}}\left(\frac{h^3}{\mu }\frac{\partial p}{\partial y}\right)\text{d}x\text{d}y={\displaystyle ∯6\left(\omega r\frac{\partial h}{\partial x}+2\frac{\partial h}{\partial t}\right)\text{d}x\text{d}y}$ (5)

在控制体内，式(5)可以改写为：

$\frac{h_e}{\mu }\frac{p_E-p_P}{\delta x_{EP}}\Delta y-\frac{h_w}{\mu }\frac{p_P-p_W}{\delta x_{PW}}\Delta y+\frac{h_n}{\mu }\frac{p_N-p_P}{\delta y_{NP}}\Delta x-\frac{h_s}{\mu }\frac{p_P-p_S}{\delta y_{PS}}\Delta x=S\Delta x\Delta y$ (6)

在控制体四个方向上有：

${\left(\frac{h^3}{\mu }\frac{\partial p}{\partial x}\right)}_e=\frac{h_e^3}{\mu }\frac{p_E-p_P}{\delta x_{EP}}$ (7)

${\left(\frac{h^3}{\mu }\frac{\partial p}{\partial x}\right)}_w=\frac{h_w^3}{\mu }\frac{p_P-p_W}{\delta x_{PW}}$ (8)

${\left(\frac{h^3}{\mu }\frac{\partial p}{\partial y}\right)}_n=\frac{h_n^3}{\mu }\frac{p_N-p_P}{\delta y_{NP}}$ (9)

${\left(\frac{h^3}{\mu }\frac{\partial p}{\partial y}\right)}_s=\frac{h_s^3}{\mu }\frac{p_P-p_S}{\delta y_{PS}}$ (10)

将式(7)~(10)代入到式(6)中并整理可以得到：

${\alpha }_P{p}_P-{\alpha }_E{p}_E-{\alpha }_W{p}_W-{\alpha }_N{p}_N-{\alpha }_S{p}_S=S$ (11)

其中：

${\alpha }_E=\frac{h_E^3}{\mu }\frac{\delta y}{\delta x}$ (12)

${\alpha }_W=\frac{h_W^3}{\mu }\frac{\delta y}{\delta x}$ (13)

${\alpha }_N=\frac{h_N^3}{\mu }\frac{\delta x}{\delta y}$ (14)

${\alpha }_S=\frac{h_S^3}{\mu }\frac{\delta x}{\delta y}$ (15)

${\alpha }_{\text{P}}={\alpha }_{\text{E}}+{\alpha }_{\text{W}}+{\alpha }_{\text{N}}+{\alpha }_{\text{S}}$ (16)

$S=-6\omega r\left(h_{\text{e}}-h_{\text{w}}\right)\delta y-12{\displaystyle \underset{w}{\overset{e}{\int }}{\displaystyle \underset{s}{\overset{n}{\int }}\frac{\partial h}{\partial t}\text{d}x\text{d}y}}$ (17)

通过有限体积法的离散，该控制体内压力值可由$ax=b$ 形式的线性方程表示，其中$a$ 是以$a_{i,j}$ 构成的系数行向量，$x$ 是以压力值构成的压力列向量。遍历流体域内所有压力节点即可构建形如$AX=B$ 的线性方程组。

根据所得线性方程组，本文采用逐次超松弛迭代法(SOR迭代法)对其进行求解。逐次超松弛迭代法实际上是对高斯–塞德尔迭代算法(GS算法)的一种改进和优化。与传统的高斯–塞德尔方法不同，逐次超松弛迭代法在每一次迭代过程中，不仅利用当前的计算结果，还会结合每一次迭代计算的误差或变化量，并对这些结果进行加权处理，再将加权后的结果作为下一次迭代的输入。因此逐次超松弛迭代法能够有效地降低收敛步数，从而显著提高了收敛速度。对于$AX=B$ 形式的线性方程组，其迭代格式表示如下：

$x_i^{(k)}=x_i^{(k-1)}+\frac{\omega }{a_{ii}}\left(b_i-{\displaystyle \sum _{j=1}^{i-1}{a}_{ij}}{x}_j^{(k)}-{\displaystyle \sum _{j=i}^n{a}_{ij}}{x}_j^{(k-1)}\right)$ (18)

式中，$\omega$ ——松弛因子。

为了获取更快的计算速度，油膜压力和厚度数据采用列向量的存储方式。在二维矩阵中，可用$a_{(i,j)}$ 来表示任意节点下存储的数值($i$ 为行编号，$j$ 为列编号)，该方式可以直观体现出每个节点之间的位置关系，某个节点的周围节点非常容易表示。采用列向量的方式存储同一节点下的压力值可用$a_{(k)}$ 表示，但采用该方式存储数据，节点与节点之间的位置关系表述不直观。因此需要构建二维矩阵与列向量之间的关系，而$i$ ，$j$ 与$k$ 之间的关系依赖遍历网格节点的顺序。本文使用列扫描的方式，将二维矩阵转化为列向量。转化关系如式(19) (20)所示(图3)。

$a_{(i,j)}=a_{(k)}$ (19)

$k=\left(j-1\right)m+i$ (20)

式中，$m$ ——二维矩阵中行数量。

Figure 3. Pressure field scan sequence

图3. 压力场扫描顺序

根据式(18)和上述扫描方式，对于流体域内部节点(四条边上以外的节点)来说，其$p_E$ 和$p_N$ 位于$p_P$ 后，尚未计算出新结果所以采用第$k-1$ 次迭代值；$p_W$ 和$p_S$ 位于$p_P$ 前，已经存在新计算结果所以采用第$k$ 次迭代值。因此在流体域内部可以获得油膜压力迭代公式。

$p_P^{(k)}=p_P^{(k-1)}+\omega \left(\frac{S+a_W{p}_W^{(k)}+a_S{p}_S^{(k)}+a_N{p}_N^{(k-1)}+a_E{p}_E^{(k-1)}}{a_P}-p_P^{(k-1)}\right)$ (21)

在原模型中，矩形流体域的左右两边原本相连，因此在边界条件设置中需将两边设置为对称边界(上下两端点除外)，如式(22)所示。由于左右边界仅为对称边界没有确定的数值约束，因此左右边界上的压力值同样需要迭代计算获得，计算获得左边界压力数值分布后可根据式(22)确定右边界压力数值分布。对于左边界上节点$p_{(i,1)}$ 来说，其$p_W$ 值与位于右边界上同行节点$p_{(i,pnx+1)}$ 的$p_W$ 值相同，但位于$p_P$ 值之后，尚未计算获得，因此$p_W$ 应选取第$k-1$ 次迭值。而$p_E$ ，$p_N$ 和$p_S$ 与式(21)中取值一致。可以构建出如式(23)所示的左边界节点压力迭代公式。至此完成了对称边界条件设置。上下边两边则设置为压力边界，有确定数值约束，如式(24) (25)所示。

$p_{(i,1)}=p_{(i,pnx+1)}$ (22)

$p_P^{(k)}=p_P^{(k-1)}+\omega \left(\frac{S+a_S{p}_S^{(k)}+a_W{p}_W^{(k-1)}+a_N{p}_N^{(k-1)}+a_E{p}_E^{(k-1)}}{a_P}-p_P^{(k-1)}\right)$ (23)

$p_{(1,j)}=p_a$ (24)

$p_{(pny+1,j)}=p_a$ (25)

式中，$p_a$ ——环境压力。

3.2. 轴心轨迹

在主轴运转过程中，主轴会受到多种载荷的作用。这些载荷包括回转部件由于残余不平衡所引起的周期性不平衡载荷、碰撞和冲击等瞬变载荷。由于这些载荷的影响，轴颈的轴心并非固定在某一位置，而是会随着时间的推移发生变化。为了更好地了解轴心的运动轨迹和行为，需要计算轴承中轴心轨迹的变化情况。本文将分析主轴在油膜压力和自身重力共同作用下的轴心轨迹情况。

在完成滑动轴承油膜计算后可利用式(26)计算获得油膜压力作用在轴颈上的作用力。

$F_p={\displaystyle \sum _{i=1}^{N_p}{p}_i{A}_{pi}}$ (26)

该作用力在三个坐标轴上的分量分别为：

$\{\begin{array}{l}{F}_{px}=-F_p\cos \theta ={\displaystyle \sum _{i=1}^{N_p}{p}_i{A}_{pi}\cos {\theta }_i}\\ F_{py}=0\\ F_{pz}=-F_p\sin \theta ={\displaystyle \sum _{i=1}^{N_p}{p}_i{A}_{pi}\sin {\theta }_i}\end{array}$ (27)

在油膜压力的作用下，轴颈会在$x$ 和$z$ 方向上产生加速度，使轴颈产生偏移。该偏移改变了轴颈原有姿态。进而导致原有油膜厚度分布情况随之变化。又由于油膜压力与油膜厚度有关，因此偏移将会导致油膜压力改变，进而会引起轴颈再次发生偏移。在滑动轴承使用过程中上述流程将不断循环。在一定条件下，该循环会达到一种动态平衡状态，轴颈偏移最终将收敛至某一位置。

4. 仿真求解

4.1. 仿真参数设定

基于上述构建的滑动轴承动力学模型及其求解算法，通过编写Matlab程序进行数值仿真。首先对稳态下结果进行仿真，由于未考虑动力学特性，油膜厚度会保持不变，因此需要人为设定偏心距初始值才能获得油膜压力分布，稳态仿真参数如表1所示。其次基于滑动轴承的动力学模型，引用文献[9]中的仿真参数进行仿真，验证该计算方法的准确性。最后基于已建立的模型开展对轴心轨迹影响因素的研究。

Table 1. Simulation parameters

表1. 仿真参数

名称	值	单位
轴承内径	30	mm
轴颈	29.9	mm
轴承长度	40	mm
x方向初始偏心距	0	mm
z方向初始偏心距	−0.025	mm
主轴转速	3000	rev/min
油液粘度	0.0277	Pa·s
环境压力	100,000	Pa

4.2. 稳态仿真结果

对上述模型进在轴颈29.9 mm，主轴转速3000 rev/min工况进行稳态仿真，图4、图5中分别展示了滑动轴承长度中点处截面(y = 20)上的压力分布情况和整个流体域中油膜压力分布情况。从图中可以看出，在0˚~270˚范围内，油膜压力逐渐增大到最大压力值，油膜最大压力为2.2 MPa，随后压力值急剧下降，在270˚~360˚范围内迅速下降至环境压力。这与经典滑动轴承的润滑理论一致。

Figure 4. Midpoint section pressure distribution

图4. 中点截面压力分布

Figure 5. Oil film pressure distribution

图5. 油膜压力分布

4.3. 动态仿真结果

从图6、图7中可以看出，文献中轴心轨迹最终收敛在偏心率为0.35，角度值为340˚位置处。本文仿真结果最终收敛于偏心率为0.37，角度值为337˚位置处。该模型的偏心率误差为5.7%，角度误差为0.8%。这一结果表明，该模型具有较高的精度，能够满足后续研究中的精度要求。

Figure 6. Simulation results

图6. 仿真结果

Figure 7. Simulation results (literature)

图7. 仿真结果(文献)

在半径间隙0.1 mm，主轴质量10 kg工况下，分析转速对轴心轨迹的影响。从图8、图9中可以看到，转速在3000 rev/min内，轴心轨迹逐渐趋于收敛，轴颈受到自重和油膜力的作用产生偏移，轴心由初始位置原点出发，不断调整，最终在点(30.36, −8.15)，(25.91, −4.86)，(22.47, −2.90)处平衡，主轴偏心率分别为0.31，0.26和0.23。此时油膜力与主轴重力、主轴外加载荷力保持平衡。轴颈中心在平衡点处做微小幅度的偏移，其位置基本稳定不变。并且随着转速逐渐升高，轴心收敛位置逐渐靠近原点。在2000~3000 rev/min工况下，轴心均在0.13 s时收敛，受转速影响较小。

Figure 8. Trajectory distribution at different speeds

图8. 不同转速下的轨迹分布情况

Figure 9. Trajectory convergence at different speeds

图9. 不同转速下的轨迹收敛情况

在主轴转速2000 rev/min，主轴质量10 kg工况下，分析配合间隙(半径间隙)对轴心轨迹的影响。从图10，图11中可以看到，在不同的配合间隙下，轴心收敛位置分别为(4.54, 0.34)，(30.36, −8.15)，(69.09, −42.34)。主轴偏心率分别为0.09，0.31和0.54。并且随着配合间隙逐渐升高，轴心收敛位置逐渐远离原点。由于配合间隙较小，轴颈的微小偏移会导致油膜压力产生较大变化，因此在小配合间隙情况下，轴心收敛过程会出现比较剧烈的抖动，收敛时间也随之加长，约为0.25 s。

Figure 10. Trajectory distribution at different mating clearance

图10. 不同配合间隙下的轨迹分布情况

Figure 11. Trajectory convergence at different mating clearance

图11. 不同配合间隙下的轨迹收敛情况

在半径间隙0.1 mm，主轴转速2500 rev/min工况下，分析主轴质量对轴心轨迹的影响。从图12，图13中可以看到，在不同的主轴质量下，轴心收敛位置分别为(25.91, −4.86)，(34.08, −11.62)，(39.82, −18.43)，主轴偏心率分别为0.26，0.36和0.44。随着主轴质量逐渐升高，轴心收敛位置逐渐远离原点，并且轴心收敛时间逐渐由0.14 s降低至0.11 s，变化幅度较小。

Figure 12. Trajectory distribution at different spindle quality

图12. 不同主轴质量下的轨迹分布情况

Figure 13. Trajectory convergence at different spindle quality

图13. 不同主轴质量下的轨迹分布情况

5. 结论

本文利用有限体积法求解油膜雷诺方程，通过编写Matlab程序实现对雷诺方程的求解。仿真结果表明，该计算方法与经典润滑理论结果吻合。通过与文献对比，该计算方法具有较高精度。基于该数值计算方法，分析了转速、主轴质量和配合间隙对轴心收敛位置和收敛速度的影响，仿真结论如下：

(1) 在本文仿真工况下，主轴轴心轨迹均达到收敛状态。

(2) 轴心收敛速度受主轴转速和主轴质量影响较小，受到配合间隙的影响较明显。具体表现为配合间隙越小，收敛时间越久。

(3) 三个参数均会对主轴最终稳定时的偏心率产生影响，可从增大柱塞泵使用转速，减小配合间隙和降低主轴质量三个方面降低主轴偏心率。

综上所述，为了实现主轴轴心的快速收敛，同时避免出现明显偏心现象。可通过调整主轴转速和主轴质量实现。然而，仅通过调整配合间隙，难以使上述两项指标同时达到最优值，因此未来可开展最佳配合间隙的研究。