考虑流体惯性的球面螺旋槽轴承特性分析

Characteristics Analysis of Spherical Spiral Groove Bearings Considering Fluid Inertia

摘要

球面螺旋槽轴承作为飞轮储能装置中转子系统的重要承载和振动控制结构，其刚度和阻尼直接影响转子系统的振动行为，进而影响飞轮储能装置的性能与效率.为研究轴承的润滑特性，基于球坐标系下考虑流体惯性的雷诺方程，采用小扰动法推导扰动压力方程，利用有限体积法和有限差分法求解偏微分方程，获得油膜压力分布.通过数值计算分析考虑流体惯性和不考虑流体惯性的油膜压力、承载能力和动态特性.结果表明，当转速和偏心率较大时，流体惯性存在明显的影响，考虑流体惯性的油膜压力、承载力和阻尼更大，刚度更小.

Abstract

Spherical spiral groove bearings （SSGBs） serve as an important bearing and vibration control structure of the rotor system in a flywheel energy storage device， and their stiffness and damping directly affect the vibration behavior of the rotor system， thus controlling the performance and efficiency of the flywheel energy storage device. To study the lubrication characteristics of the bearings， based on the Reynolds equation considering fluid inertia in the spherical coordinate system， the perturbation pressure equation was derived by using the perturbation method. Then， the partial differential equations were solved by using the finite volume and the finite difference methods to obtain the pressure distribution of the oil film. Furthermore， the pressure， load-carrying capacity， and dynamic characteristics of the oil film considering the fluid inertia and disregarding fluid inertia were analyzed by numerical calculations. The results show that when the rotational speed and eccentricity are larger， the fluid inertia has an obvious influence； the oil film pressure， load-carrying capacity， and damping considering fluid inertia are larger， and the stiffness is smaller.

Graphical abstract

关键词

球面螺旋槽轴承 / 流体惯性 / 动态特性 / 油膜压力 / 小扰动法

Key words

spherical spiral groove bearing / fluid inertia / dynamic characteristic / oil film pressure / perturbation method

球面螺旋槽轴承（spherical spiral groove bearing， SSGB）是一种特殊设计的轴承，其内部具有沿球面螺旋状排列的凹槽，可以同时承受径向与轴向载荷，具有自润滑、高刚度和低摩擦等优点^［1］，在飞轮储能装置上有着广泛的应用.SSGB作为重要的支承结构，其油膜刚度与阻尼特性对于转子系统振动行为有着很大的影响.但是目前的研究主要集中在气体润滑上，缺少对于油润滑SSGB的分析，油和气体润滑存在着较大的差异.因此，有必要建立油润滑SSGB的润滑模型，探究其动态特性的变化规律.

学者对SSGB流体润滑模型进行了大量研究.Zhang等^［2］研究了柔性气体轴颈轴承转子的动力学特性；Duan等^［3］对多微孔气润滑静压推力轴承的静特性进行了仿真，分析了其承载能力.然而，上述研究均为气体润滑轴承，没有考虑油润滑.Xue等^［4］建立了球面混合陶瓷滑动轴承的瞬态载荷转子模型，研究了加速过程中轴承的稳定性；Sato等^［5］采用小扰动法分析了球面静压轴承的动态特性；Huang等^［6］分析了人字槽高速混合静压轴承的静态与动态特性.然而，上述工作均未考虑螺旋槽的影响.Yang等^［7］建立了螺旋微槽耦合推力轴承的摩擦动力学模型，研究了螺旋微槽对轴承性能的影响；Jiang等^［8］基于层流空化流模型，研究了水润滑织构螺旋槽止推轴承的静态特性；Xu等^［9］分析了水润滑流体动压人字形槽轴颈轴承低速运转时的动力学特性.然而，上述研究的对象均为推力轴承，并未涉及球面轴承.

本文以球面螺旋槽轴承为研究对象，建立了考虑惯性效应的雷诺方程，通过小扰动法得到稳态雷诺方程的扰动压力，并采用有限体积法与有限差分法相结合的方法，得到稳态压力与扰动压力分布.数值计算获得了SSGB的刚度和阻尼系数，并探究了轴承动态特性系数的变化规律，以期为SSGB的设计及轴承稳定性分析提供一定的理论基础.

1 SSGB润滑分析数学模型

为了探究SSGB的润滑特性，采用牛顿流体假设，忽略压力沿半径方向的变化，通过小扰动法建立了考虑流体惯性的油膜压力方程.

1.1 轴承的结构

SSGB的结构如图1所示.在半球面上有均匀分布的螺旋槽，其中h₀为转子无偏心时的油膜厚度，h_g为槽深，b_r为台区宽度，b_g为槽区宽度，θ_g为螺旋槽起始角，β为螺旋角，R为轴承半径，ω为转子转速，槽台宽度比γ=b_g/b_r.

1.2 轴承润滑的数学模型

在考虑流体惯性的纳维-斯托克斯方程和连续性方程的基础上，假设黏度不变、忽略油膜压力沿膜厚方向的变化，建立球面坐标系下的雷诺方程^［5］如式（1）~式（3）所示.

\begin{array}{c}\frac{\partial }{\partial \phi }\left(H^{\mathrm{3}}\frac{\partial p_{\mathrm{0}}}{\partial \phi }\right)+\mathrm{s}\mathrm{i}\mathrm{n}\theta \frac{\partial }{\partial \theta }\left(\mathrm{s}\mathrm{i}\mathrm{n}\theta H^{\mathrm{3}}\frac{\partial p_{\mathrm{0}}}{\partial \theta }\right)=\\ \Lambda \mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta \frac{\partial H}{\partial \phi }+\mathrm{2}\Lambda \mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta \frac{\partial H}{\partial t}，\end{array}(1)

\begin{array}{c}\frac{\partial }{\partial \phi }\left(H^{\mathrm{3}}\frac{\partial p_{\mathrm{1}}}{\partial \phi }\right)+\mathrm{s}\mathrm{i}\mathrm{n}\theta \frac{\partial }{\partial \theta }\left(\mathrm{s}\mathrm{i}\mathrm{n}\theta H^{\mathrm{3}}\frac{\partial p_{\mathrm{1}}}{\partial \theta }\right)=\\ \Lambda \mathrm{s}\mathrm{i}\mathrm{n}\theta \left(\frac{\partial S_{\mathrm{1}}}{\partial \theta }+\frac{\partial S_{\mathrm{2}}}{\partial \phi }\right)，\end{array}(2)

p=p_{\mathrm{0}}+Re\times p_{\mathrm{1}}.(3)

其中：轴承数Λ=6μωR²/h₀²p_a；雷诺数Re=ρh₀²ω/μ；无量纲油膜厚度H=h/h₀；μ为润滑油的动力粘度；ρ为润滑油的密度；p_a为标准大气压；p₀为不考虑惯性的无量纲压力；p₁为惯性引起的压力变化；p为考虑惯性后的压力；θ为子午线夹角；φ为周向夹角；t为时间；S₁和S₂的表达式见文献［5］.

转子偏心示意图如图2所示，稳态下轴承台区和槽区任意一点（φ，θ）的无量纲油膜厚度H表达式如下：

H=\mathrm{1}+{\epsilon }_x\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\phi +{\epsilon }_y\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\phi +{\epsilon }_z\mathrm{s}\mathrm{i}\mathrm{n}\theta ,(4)

H=\mathrm{1}+\frac{h_{\mathrm{g}}}{h_{\mathrm{0}}}+{\epsilon }_x\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\phi +{\epsilon }_y\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\phi +{\epsilon }_z\mathrm{s}\mathrm{i}\mathrm{n}\theta .(5)

式中：偏心率ε_x =e_x/h₀，ε_y =e_y/h₀，ε_z =e_z/h₀，e_x，e_y，e_z 分别为轴承在x，y，z三个方向的偏心距.

2 轴承动态特性的求解

采用保角变换法简化雷诺方程，利用小扰动法推导扰动压力方程，根据润滑状态给出了连续性边界条件和恒压边界条件.

2.1 扰动压力方程

式（1）、式（2）所示的雷诺方程形式过于复杂，不易在球面坐标系下求解，为此引入保角变换

\alpha =-\mathrm{l}\mathrm{n}\left(\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}\right),\mathrm{d}\alpha =-\mathrm{d}\theta /\mathrm{s}\mathrm{i}\mathrm{n}\theta .(6)

式中，α为变换后的坐标.

代入式（1）、式（2）可得直角坐标下的雷诺方程：

\frac{\partial }{\partial \phi }\left(H^{\mathrm{3}}\frac{\partial p_{\mathrm{0}}}{\partial \phi }\right)+\frac{\partial }{\partial \alpha }\left(H^{\mathrm{3}}\frac{\partial p_{\mathrm{0}}}{\partial \alpha }\right)=\Lambda {\lambda }^{\mathrm{2}}\frac{\partial H}{\partial \phi }+\mathrm{2}\Lambda \frac{\partial H}{\partial t},(7)

\frac{\partial }{\partial \phi }\left(H^{\mathrm{3}}\frac{\partial p_{\mathrm{1}}}{\partial \phi }\right)+\frac{\partial }{\partial \alpha }\left(H^{\mathrm{3}}\frac{\partial p_{\mathrm{1}}}{\partial \alpha }\right)=\Lambda \frac{\partial S_{\mathrm{1}}}{\partial \alpha }+\Lambda \lambda \frac{\partial S_{\mathrm{2}}}{\partial \theta }.(8)

其中，变换参数λ=sinθ=2e ^α /（1+e^2α ），α∈［α₁，α₂］.

当θ=0时，式（8）所示的保角变换不存在，为此可以引入一个很小的角度θ_a.由于在球窝中心流场压力的连续性，其驻点附近微小邻域的压力同驻点的压力非常接近，故近似用θ_a处的压力近似驻点的压力^［10］.

为了求解轴承的动态刚度和阻尼，给油膜3个方向的偏心率分别施加小扰动Δε_x，Δε_y，Δε_z，计算扰动引起的油膜厚度变化，从而建立扰动压力方程.

假设轴承稳态条件下的压力为p_s，无量纲油膜厚度为H_s，忽略高阶项，施加小扰动后无量纲压力和油膜厚度为

\begin{array}{c}p=p_{\mathrm{s}}+\mathrm{\Delta }p=\\ p_{\mathrm{s}}+\mathrm{\Sigma }{p}_{\xi }+\mathrm{\Sigma }{p}_{\dot{\xi }},\xi =x,y,z.\end{array}(9)

式中：\mathrm{\Delta }p为施加小扰动后无量纲压力的变化；p_{\xi }与p_{\dot{\xi }}分别为偏心率变化与偏心率一阶导数变化引起的扰动压力；\xi =x,y,z为不同方向的坐标；\dot{\xi }=\partial \xi /\partial t为坐标的一阶导数.

\begin{array}{c}H=H_{\mathrm{s}}+\mathrm{\Delta }H=H_{\mathrm{s}}+\mathrm{\Delta }{\epsilon }_x\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\phi +\\ \mathrm{\Delta }{\epsilon }_y\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\phi +\mathrm{\Delta }{\epsilon }_z\mathrm{s}\mathrm{i}\mathrm{n}\theta .\end{array}(10)

将式（9）、式（10）代入式（1）~式（3）得

\begin{array}{c}\frac{\partial }{\partial \phi }\left({\left(H_{\mathrm{s}}+\mathrm{\Delta }H\right)}^{\mathrm{3}}\frac{\partial \left(p_{\mathrm{s}\mathrm{0}}+\mathrm{\Delta }{p}_{\mathrm{s}\mathrm{0}}\right)}{\partial \phi }\right)+\\ \frac{\partial }{\partial \alpha }\left({\left(H_{\mathrm{s}}+\mathrm{\Delta }H\right)}^{\mathrm{3}}\frac{\partial \left(p_{\mathrm{s}\mathrm{0}}+\mathrm{\Delta }{p}_{\mathrm{s}\mathrm{0}}\right)}{\partial \alpha }\right)=\\ \Lambda {\lambda }^{\mathrm{2}}\frac{\partial \left(H_{\mathrm{s}}+\mathrm{\Delta }H\right)}{\partial \phi }+\mathrm{2}\Lambda {\lambda }^{\mathrm{2}}\frac{\partial \left(H_{\mathrm{s}}+\mathrm{\Delta }H\right)}{\partial t},\end{array}(11)

\begin{array}{c}\frac{\partial }{\partial \phi }\left({\left(H_{\mathrm{s}}+\mathrm{\Delta }H\right)}^{\mathrm{3}}\frac{\partial \left(p_{\mathrm{s}\mathrm{1}}+\mathrm{\Delta }{p}_{\mathrm{s}\mathrm{1}}\right)}{\partial \phi }\right)+\\ \frac{\partial }{\partial \alpha }\left({\left(H_{\mathrm{s}}+\mathrm{\Delta }H\right)}^{\mathrm{3}}\frac{\partial \left(p_{\mathrm{s}\mathrm{1}}+\mathrm{\Delta }{p}_{\mathrm{s}\mathrm{1}}\right)}{\partial \alpha }\right)=\\ \Lambda \frac{\partial (S_{\mathrm{1}}+\mathrm{\Delta }{S}_{\mathrm{1}})}{\partial \alpha }+\Lambda \lambda \frac{\partial (S_{\mathrm{2}}+\mathrm{\Delta }{S}_{\mathrm{2}})}{\partial \phi }.\end{array}(12)

\left.\begin{array}{c}{p}_{\mathrm{s}}=p_{\mathrm{s}\mathrm{0}}+Re\times p_{\mathrm{s}\mathrm{1}}；\\ \mathrm{\Delta }p=\mathrm{\Delta }{p}_{\mathrm{0}}+Re\times \mathrm{\Delta }{p}_{\mathrm{1}}.\end{array}\right\}(13)

其中：\mathrm{\Delta }H为施加小扰动后的油膜厚度变化；p_s0为不考虑惯性的稳态压力；p_s1为惯性引起的压力变化；\mathrm{\Delta }{p}_{\mathrm{0}}，\mathrm{\Delta }{p}_{\mathrm{1}}分别为小抗动引起的不考虑惯性的压力变化与惯性压力变化.

将式（11）、式（12）展开，并对比Δε_x，Δε_y，Δε_z 前的系数，可得扰动压力方程组为

\left.\begin{array}{l}\frac{\partial }{\partial \phi }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\mathrm{s}\mathrm{0}}}{\partial \phi }\right)+\frac{\partial }{\partial \alpha }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\mathrm{s}\mathrm{0}}}{\partial \alpha }\right)=\Lambda {\lambda }^{\mathrm{2}}\frac{\partial H_{\mathrm{s}}}{\partial \phi };\\ \frac{\partial }{\partial \phi }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\xi \mathrm{0}}}{\partial \phi }\right)+\frac{\partial }{\partial \alpha }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\xi \mathrm{0}}}{\partial \alpha }\right)=\Lambda {\lambda }^{\mathrm{2}}\frac{\partial H_{\xi }}{\partial \phi }-\frac{\partial }{\partial \phi }\left(\mathrm{3}{H}_{\mathrm{s}}^{\mathrm{2}}{H}_{\xi }\frac{\partial p_{\mathrm{s}\mathrm{0}}}{\partial \phi }\right)-\frac{\partial }{\partial \alpha }\left(\mathrm{3}{H}_{\mathrm{s}}^{\mathrm{3}}{H}_{\xi }\frac{\partial p_{\mathrm{s}\mathrm{0}}}{\partial \alpha }\right),\xi =x,y,z;\\ \frac{\partial }{\partial \phi }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\dot{\xi }\mathrm{0}}}{\partial \phi }\right)+\frac{\partial }{\partial \alpha }\left(H_s^{\mathrm{3}}\frac{\partial p_{\dot{\xi }\mathrm{0}}}{\partial \alpha }\right)=\mathrm{2}\Lambda {\lambda }^{\mathrm{2}}{H}_{\xi },\xi =x,y,z;\\ \frac{\partial }{\partial \phi }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\mathrm{s}\mathrm{1}}}{\partial \phi }\right)+\frac{\partial }{\partial \alpha }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\mathrm{s}\mathrm{1}}}{\partial \alpha }\right)=\Lambda \frac{\partial S_{\mathrm{1}}}{\partial \alpha }+\Lambda \lambda \frac{\partial S_{\mathrm{2}}}{\partial \phi };\\ \frac{\partial }{\partial \phi }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\xi \mathrm{1}}}{\partial \phi }\right)+\frac{\partial }{\partial \alpha }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\xi \mathrm{1}}}{\partial \alpha }\right)=\Lambda \frac{\partial S_{\mathrm{1}\xi }}{\partial \alpha }+\Lambda \lambda \frac{\partial S_{\mathrm{2}\xi }}{\partial \phi },\xi =x,y,z,\dot{x},\dot{y},\dot{z};\\ S_{\mathrm{1}\xi }=\frac{\partial S_{\mathrm{1}}}{\partial \mathrm{\Delta }{\epsilon }_{\xi }},S_{\mathrm{2}\xi }=\frac{\partial S_{\mathrm{2}}}{\partial \mathrm{\Delta }{\epsilon }_{\xi }},\xi =x,y,z,\dot{x},\dot{y},\dot{z};\\ p_{\xi }=p_{\xi \mathrm{0}}+Re\times p_{\xi \mathrm{1}},\xi =x,y,z,\dot{x},\dot{y},\dot{z};\\ H_x=\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\phi ,H_y=\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\phi ,H_z=\mathrm{s}\mathrm{i}\mathrm{n}\theta .\end{array}\right\}(14)

式中：S_1ξ，S_2ξ 为S₁，S₂对Δε_ξ 的偏导数；Δε_ξ 为不同坐标下偏心率的变化；p_ξ0为不考虑惯性时的扰动压力；p_ξ1为惯性引起的扰动压力变化；H_x，H_y，H_z 分别为无量纲油膜厚度H对Δε_x，Δε_y，Δε_z 的偏导数.

求解式（14）的边界条件如下：

\left.\begin{array}{l}{p}_{\mathrm{s}\mathrm{0}}\left(\phi ,{\alpha }_{\mathrm{2}}\right)=p_{\mathrm{s}\mathrm{1}}\left(\phi ,{\alpha }_{\mathrm{2}}\right)=\mathrm{0};\\ p_{\xi \mathrm{0}}\left(\phi ,{\alpha }_{\mathrm{2}}\right)=p_{\xi \mathrm{1}}\left(\phi ,{\alpha }_{\mathrm{2}}\right)=\mathrm{0};\\ p_{\mathrm{s}\mathrm{0}}\left(\phi ,\alpha \right)=p_{\mathrm{s}\mathrm{0}}\left(\phi +\mathrm{2}\pi ,\alpha \right);\\ p_{\xi \mathrm{0}}\left(\phi ,\alpha \right)=p_{\xi \mathrm{0}}\left(\phi +\mathrm{2}\pi ,\alpha \right);\\ p_{\mathrm{s}\mathrm{1}}\left(\phi ,\alpha \right)=p_{\mathrm{s}\mathrm{1}}\left(\phi +\mathrm{2}\pi ,\alpha \right);\\ p_{\xi \mathrm{1}}\left(\phi ,\alpha \right)=p_{\xi \mathrm{1}}\left(\phi +\mathrm{2}\pi ,\alpha \right);\\ \frac{\partial p_{\mathrm{s}\mathrm{0}}}{\partial \alpha }{|}_{\alpha ={\alpha }_{\mathrm{1}}}=\frac{\partial p_{\mathrm{s}\mathrm{1}}}{\partial \alpha }{|}_{\alpha ={\alpha }_{\mathrm{1}}}=\mathrm{0};\\ \frac{\partial p_{\xi \mathrm{0}}}{\partial \alpha }{|}_{\alpha ={\alpha }_{\mathrm{1}}}=\frac{\partial p_{\xi \mathrm{1}}}{\partial \alpha }{|}_{\alpha ={\alpha }_{\mathrm{1}}}=\mathrm{0}.\end{array}\right\}(15)

式中：{\alpha }_{\mathrm{1}}、{\alpha }_{\mathrm{2}}为求解边界处的坐标，根据式（6）可以得到{\alpha }_{\mathrm{1}}=\mathrm{0},{\alpha }_{\mathrm{2}}=-\mathrm{l}\mathrm{n}\left(\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\theta }_{\mathrm{a}}}{\mathrm{2}}\right).

2.2 有限差分法求解

油膜厚度的不连续性导致膜厚偏导数项在槽台交界处是无意义的，故采用九点差分格式的有限体积法进行离散处理^［11］，网格划分示意图如图3所示.图中i，j为网络格点的编号，L_ij 为线积分路径.

在控制体积Ω_ij 上对雷诺方程进行积分

\begin{array}{c}{∯}_{{\Omega }_{ij}}^{}\frac{\partial }{\partial \phi }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\mathrm{s}\mathrm{0}}}{\partial \phi }\right)+\frac{\partial }{\partial \alpha }\left(H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\mathrm{s}\mathrm{0}}}{\partial \alpha }\right)=\\ {∯}_{{\Omega }_{ij}}^{}\Lambda {\lambda }^{\mathrm{2}}\frac{\partial H_{\mathrm{s}}}{\partial \phi }.\end{array}(16)

由格林公式可得

{\oint }_{L_{i,j}}^{}{H}_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\mathrm{s}\mathrm{0}}}{\partial \phi }\mathrm{d}\alpha -H_{\mathrm{s}}^{\mathrm{3}}\frac{\partial p_{\mathrm{s}\mathrm{0}}}{\partial \alpha }\mathrm{d}\phi ={\oint }_{L_{i,j}}^{}\Lambda {\lambda }^{\mathrm{2}}{H}_{\mathrm{s}}\mathrm{d}\alpha .(17)

根据中矩形积分公式得到有限差分格式

\begin{array}{c}\left[H_{{\mathrm{s}}_{i+\mathrm{1}/\mathrm{2},j}}^{\mathrm{3}}\left(P_{\mathrm{s}{\mathrm{0}}_{i+\mathrm{1},j}}-P_{\mathrm{s}{\mathrm{0}}_{i,j}}\right)-H_{{\mathrm{s}}_{i-\mathrm{1}/\mathrm{2},j}}^{\mathrm{3}}\left(P_{\mathrm{s}{\mathrm{0}}_{i,j}}-P_{\mathrm{s}{\mathrm{0}}_{i-\mathrm{1},j}}\right)\right]\frac{\mathrm{\Delta }\alpha }{\mathrm{\Delta }\phi }+\\ \left[H_{{\mathrm{s}}_{i,j+\mathrm{1}/\mathrm{2}}}^{\mathrm{3}}\left(P_{\mathrm{s}{\mathrm{0}}_{i,j+\mathrm{1}}}-P_{\mathrm{s}{\mathrm{0}}_{i,j}}\right)-H_{{\mathrm{s}}_{i,j-\mathrm{1}/\mathrm{2}}}^{\mathrm{3}}\left(P_{\mathrm{s}{\mathrm{0}}_{i,j}}-P_{\mathrm{s}{\mathrm{0}}_{i,j-\mathrm{1}}}\right)\right]\frac{\mathrm{\Delta }\phi }{\mathrm{\Delta }\alpha }=\\ \Lambda {\lambda }^{\mathrm{2}}\left(H_{{\mathrm{s}}_{i+\mathrm{1}/\mathrm{2},j}}-H_{{\mathrm{s}}_{i-\mathrm{1}/\mathrm{2},j}}\right)\mathrm{\Delta }\alpha .\end{array}(18)

式中，下标对应图3中的网络格点.

式（14）中的其他方程也可建立类似的有限差分格式，采用超松弛（successive over-relaxation，SOR）迭代法求解即可得到扰动压力.

2.3 油膜刚度与阻尼的计算

油膜不同方向的承载能力可通过压力积分计算获得

\left.\begin{array}{l}{F}_x={\int }_{{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{2}}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{p}_a{R}^{\mathrm{2}}{p}_s{\lambda }^{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\phi \mathrm{d}\phi \mathrm{d}\alpha ;\\ F_y={\int }_{{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{2}}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{p}_a{R}^{\mathrm{2}}{p}_s{\lambda }^{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\phi \mathrm{d}\phi \mathrm{d}\alpha ;\\ F_{\mathrm{z}}={\int }_{{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{2}}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{p}_a{R}^{\mathrm{2}}{p}_s{\lambda }^{\mathrm{2}}\sqrt[]{\mathrm{1}-{\lambda }^{\mathrm{2}}}\mathrm{d}\phi \mathrm{d}\alpha .\end{array}\right\}(19)

对油膜力求导可得到刚度矩阵 K 和阻尼矩阵 C ：

\mathit{K}=\left[\begin{array}{ccc}{K}_{xx}& K_{xy}& K_{xz}\\ K_{yx}& K_{yy}& K_{yz}\\ K_{zx}& K_{zy}& K_{zz}\end{array}\right]={\int }_{{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{2}}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}\frac{p_{\mathrm{a}}{R}^{\mathrm{2}}}{h_{\mathrm{0}}}\left[\begin{array}{ccc}{p}_x{\lambda }^{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\phi & p_y{\lambda }^{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\phi & p_z{\lambda }^{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\phi \\ p_x{\lambda }^{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\phi & p_y{\lambda }^{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\phi & p_z{\lambda }^{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\phi \\ p_x{\lambda }^{\mathrm{2}}\sqrt[]{\mathrm{1}-{\lambda }^{\mathrm{2}}}& p_y{\lambda }^{\mathrm{2}}\sqrt[]{\mathrm{1}-{\lambda }^{\mathrm{2}}}& p_z{\lambda }^{\mathrm{2}}\sqrt[]{\mathrm{1}-{\lambda }^{\mathrm{2}}}\end{array}\right]\mathrm{d}\phi \mathrm{d}\alpha ,(20)

\mathit{C}=\left[\begin{array}{ccc}{C}_{xx}& C_{xy}& C_{xz}\\ C_{yx}& C_{yy}& C_{yz}\\ C_{zx}& C_{zy}& C_{zz}\end{array}\right]={\int }_{{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{2}}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}\frac{p_{\mathrm{a}}{R}^{\mathrm{2}}}{h_{\mathrm{0}}\omega }\left[\begin{array}{ccc}{p}_{\dot{x}}{\lambda }^{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\phi & p_{\dot{y}}{\lambda }^{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\phi & p_{\dot{z}}{\lambda }^{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\phi \\ p_{\dot{x}}{\lambda }^{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\phi & p_{\dot{y}}{\lambda }^{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\phi & p_{\dot{z}}{\lambda }^{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\phi \\ p_{\dot{x}}{\lambda }^{\mathrm{2}}\sqrt[]{\mathrm{1}-{\lambda }^{\mathrm{2}}}& p_{\dot{y}}{\lambda }^{\mathrm{2}}\sqrt[]{\mathrm{1}-{\lambda }^{\mathrm{2}}}& p_{\dot{z}}{\lambda }^{\mathrm{2}}\sqrt[]{\mathrm{1}-{\lambda }^{\mathrm{2}}}\end{array}\right]\mathrm{d}\phi \mathrm{d}\alpha .(21)

编程求解轴承特性时，计算流程如图4所示.图中p_{i,j}^{(k)}表示第k次迭代中p在（i，j）节点的值.

3 轴承静态特性的仿真分析

根据超松弛迭代法编程计算求解稳态压力方程，使用1/3辛普森法则计算分析惯性对压力及油膜力的影响.

3.1 油膜压力分布

根据表1的轴承参数计算，轴承在转速ω=200 r/min时的油膜压力分布如图5所示，在ω=8 000 r/min时油膜压力分布如图6所示，图7为\epsilon =\mathrm{0.2}时的油膜厚度.

通过有限元分析验证了模型的正确性.图5a为有限元分析结果，与图5b的有限差分模型结果对比可以看出，两者得到的油膜压力分布相似，两者的最大与最小压力均相差不大，说明了有限差分模型的正确性.

由从图5~图7可见，油膜厚度越小，压力越大，压力在台区明显增大，在槽区急剧下降，产生负压，油膜压力在每个槽区和台区呈现相同的变化规律.由于轴承偏心的存在，靠近偏心方向的台区油膜压力最大，从最大压力开始，沿周向台区的最大压力逐渐减小.此外，在低转速下惯性效应对于压力的影响并不明显；在高转速下雷诺数增大，流体惯性使油膜压力有明显的提高.

3.2 油膜力分析

假设轴承偏心率ε=0.2，油膜力与转速的关系如图8所示.可以看出，油膜力随着转速的增大而增大，油膜力和转速近似呈线性关系.此外，惯性会使油膜力增大，随着转速上升，惯性影响更加明显，转速直接影响了油膜惯性.

图9为转速ω=200 r/min时，不同偏心率下油膜力的变化曲线.可以看出，油膜力随偏心率增大呈非线性增大趋势，偏心率越大，油膜力变化越快.随着偏心率增加，考虑惯性与否对油膜力的影响更加显著.

4 轴承动态特性的仿真分析

根据得到的扰动压力方程，讨论了不同转速和偏心率下SSGB的动态特性.

4.1 转速对动态特性的影响

图10为ε=0.2时不同转速下油膜的直接刚度.油膜惯性使直接刚度略微下降，但影响很小，可以忽略不计.x方向的直接刚度最大，z方向的刚度较小.刚度随转速的增加近似线性增大.

图11为ε=0.2时不同转速下油膜的直接阻尼.可以看出，油膜惯性使直接阻尼略微提高.其中，z方向的直接阻尼最大，x，y方向的阻尼较小，转速对直接阻尼的影响很小.

rotational speeds

（a）—不考虑惯性； （b）—考虑惯性.

图12为交叉刚度随转速的变化曲线，受螺旋槽与球面轴承的影响，油膜存在较大的交叉刚度，K_xy 和K_zx 较大，K_yz 较小，这与螺旋槽的旋向相关，且刚度随转速的增大而近似线性增大.

4.2 轴承偏心率对动态特性的影响

图13为ω=200 r/min时，直接刚度与偏心率的关系.从整体上来看，刚度随着偏心率的增大呈现非线性增大的趋势，z方向的刚度最小，在ε=0.8附近时，y方向的刚度超过x方向.油膜惯性对刚度的影响较小，且随着偏心率的增加而增加.

图14为直接阻尼随偏心率的变化趋势，直接阻尼随着偏心率的增大呈现非线性增大的趋势，z方向的阻尼最大，x，y方向的阻尼相对较小.

eccentricities

（a）—不考虑惯性； （b）—考虑惯性.

图15为交叉刚度随偏心率的变化趋势，从图中可以看出，随着偏心率的变化，K_zx，K_yz 基本没有改变，而K_xy 呈现非线性增加.

5 结 论

1） 当转速较小时，流体惯性对压力的影响较小，可以忽略不计；而在转速较大时，由于雷诺数的增大，油膜压力有明显的提升.

2） 轴承的承载能力随转速的提升近似线性增长，随偏心率增大呈现非线性增加的趋势.流体惯性会提高油膜力的大小，在较大转速和偏心率的条件下，惯性对承载能力的影响不可忽略.

3） 转速的增加可提高轴承的刚度，转速对轴承刚度有较大影响，而阻尼基本不随转速变化.惯性效应对动态特性的影响较小.

参考文献

[1] Muijderman E A. Spiral Groove Bearings［M］. Boston： Springer-Verlag， 2013： 2-5.

[2] Zhang C B， Ao H R， Cheng Y J， et al. Theoretical and experimental studies on nonlinear dynamic performances of gas journal bearing with three pads-rotor system［J］. Nonlinear Dynamics， 2025， 113（9）： 9331-9349.

[3] Duan L Y， Xu M M. Study on static performance of gas-lubricated thrust bearing based on multi-microporous stainless steel plate［J］. Journal of the Brazilian Society of Mechanical Sciences and Engineering， 2021， 43（5）： 1-14.

[4] Xue W Q， Shen J F， Fan Y W. Study on dynamic characteristics of acceleration process of spherical hybrid sliding bearings rotor system［J］. Industrial Lubrication and Tribology， 2024， 76（3）： 365-373.

[5] Sato Y， Tamura A. Fluid inertia effects on the characteristics of spherical hydrostatic bearings： part 1 bearing forces［J］. Bulletin of JSME， 1983， 26（216）： 1059-1065.

[6] Huang X， Lin S Y， Liang J H. Modelling approach of hybrid air herringbone grooved journal bearing considering the angular misalignment and centrifugal expansion［J］. Tribology International， 2023， 189： 108928.

[7] Yang T Y， Xiao K， Xiang G， et al. Effects of rotating microgroove on tribo-dynamic performance of coupled bearings［J］. International Journal of Mechanical Sciences， 2024， 279： 109629.

[8] Jiang S Y， Liu P F， Lin X H. Study on static characteristics of water-lubricated textured spiral groove thrust bearing using laminar cavitating flow lubrication model［J］. Journal of Tribology， 2022， 144（4）： 041803.

[9] Xu G， Jiang S Y， Zhang C B， et al. Dynamic characteristics of hybrid water-lubricated herringbone groove journal bearing［J］. Industrial Lubrication and Tribology， 2025， 77（2）： 219-230.

[10] 郭良斌，王祖温，包钢，等. 新型环面节流静压气体球轴承压力分布的有限元计算［J］. 摩擦学学报， 2004， 24（6）： 531-535.

[11] Guo Liang-bin， Wang Zu-wen， Bao Gang， et al. Finite element analysis of the pressure distribution of externally pressurized spherical gas bearings with inherent compensation［J］. Journal of Tribology， 2004， 24（6）： 531-535.

[12] Chen S K， Chou H C， Kang Y. Stability analysis of hydrodynamic bearing with herringbone grooved sleeve［J］. Tribology International， 2012， 55： 15-28.

基金资助

国家自然科学基金资助项目(12272089)

粒子输运与富集技术全国重点实验室资助项目(WZKF-2024-6)

辽宁省“揭榜挂帅”科技攻关计划项目(2023JH1/10400068)

中央高校基本科研业务费专项资金资助项目(N2303002)