液体静压轴承中油膜除承受外力载荷、减轻摩擦、消除磨损外，从另一个角度讲，它也是液压主轴动力学系统中的一个关键环节。将液体静压主轴系统中的油膜看作线性化的弹簧和阻尼作用，通过4个动态刚度系数和4个阻尼系数来衡量系统的动力学特性。轴心轨迹是另外一种评判主轴系统稳定性的方法，用一种更加直观地方式反映了主轴瞬时运动情况及运动稳定性。轴承间隙中油膜厚度通常处于微尺度尺寸范围或更低，不同的油膜厚度直接关系到润滑效果和轴承性能^[1]。固-液界面上的速度滑移现象同时得到了理论和实验验证^[2-4]。速度滑移对轴承静特性必然会产生一定的影响^[5-7]，但对液压主轴中轴承动态特性的影响目前研究的还比较少。李凌祥^[8]采用单组周向均匀分布的4个弹簧-阻尼单元C0MBIN14模拟液体静压轴承，每个单元的刚度值为其理论值的1/3，分别对主轴-轴承系统进行模拟分析。杨小高等^[9]建立固定节流径向滑动轴承动态模型，进行动态特性分析。江桂云等^[10]建立了静压轴承系统的流动连续性方程，结合轴承-主轴系统动力学方程，推导出静压轴承系统动态特性的传递函数，分析了小孔节流静压轴承的动态特性。商远杰等^[11]应用弹簧阻尼单元模拟静压轴承，建立电主轴传动装置的有限元模型，分析其动态性能。Lee等^[12]运用实验手段并通过动力学8个特征系数对柔性转子系统的动力学特性进行了分析。谢沛霖等^[13]求解传统雷诺方程通过8个特征系数对毛细管节流的静压轴承的动力特性进行了分析，但没有进一步求解转子的运行状态。Rowe^[14]指出8个特征参数可以相对合理的用于轴承的动力学设计，并用来预测主轴旋转时的状态。轴心轨迹的计算方法通常有两种，一种称为汉氏法，另一种称为荷氏法。帅旗^[15]以荷氏法为理论基础推导提出三油叶轴承的动载轴心轨迹求解方法。Ramesh等^[16]以轴心轨迹作为轴承稳定性的评判依据，研究了轴承内表面粗糙度对轴承稳定性的影响。

作者把液体静压轴承的研究拓展到微尺度，在求解动力特征系数的同时，进一步将动力特征系数用于液压主轴轴心轨迹的建模，分析主轴轴心的运行状态。通过求解不定常工况下油膜压力分布的雷诺方程得非线性油膜力，求解出能够反映油膜动态特性8个特征系数，通过8个特征系数对主轴的动态特性进行分析，并将计入滑移情况下的动态特征系数与无滑移下的特征系数进行对比，分析微尺度速度滑移现象对主轴系统动态特性的影响规律。建立液压主轴存在质量偏心下的动力系统模型，并将刚度系数和阻尼系数用于液压主轴的轴心轨迹模型的建立，分析不同偏心质量距下主轴轴心的运行轨迹，对维持和提高液压主轴系统的动力学稳定性提供理论依据。

1 液体静压轴承-转子系统动力学模型 1.1 几何模型

液体静压轴承的结构如图 1所示，在一定的外载荷的作用下，转子轴心和轴承中心不再重合，而是产生了一定的偏差，转子在外载荷的作用下最终的平衡位置为O₁，两中心产生的位置偏差为e，角度偏差为θ，描述位置的角度为φ，设润滑油以恒定的油压P_s流入1、2、3、4腔。图 1坐中，x、y、z分别为轴承油膜的周向、径向、轴向坐标。

1.2 不定常工况下雷诺方程

液体静压轴承的动力学特性，反映的是转子在偏离原平衡位置时，并在此位置附近做微幅振动时，相应的油膜力的变化情况。因此，必须以不定常工况的雷诺方程作为分析、计算的基础。

为了便于计算机进行运算，对方程传统的雷诺方程进行无量纲化，取特征压强为p₀，特征滑移长度b和油膜特征厚度为轴承半径间隙h₀，油膜特征轴向长度为轴承长度的一半L/2，油膜特征周向长度为轴承半径R，令p=p₀P，其中, ${p_0} = \frac{{6w\eta {R^2}}}{{h_0^2}}$ ，b=Bh₀，h=Hh₀， $z = \lambda \frac{L}{2}\left( { - 1 \le \lambda \le 1} \right)$ ，x=Rφ(0≤φ≤2π)， $\left( {X,Y} \right) = \frac{{\left( {x,y} \right)}}{{{h_0}}}$ ， $\left( {\dot X,\dot Y} \right) = \frac{{\left( {\dot x,\dot y} \right)}}{{{h_0}{w_0}}}$ ，τ=w₀t，将其代入传统雷诺方程式得到Reynolds方程的无量纲形式为：

$ \begin{array}{l} \left( {1 + 6B} \right)\frac{\partial }{{\partial \varphi }}\left( {{H^3}\frac{{\partial P}}{{\partial \varphi }}} \right) + \\ \left( {1 + 6B} \right)\left( {\frac{{2R}}{L}} \right)\frac{\partial }{{\partial \lambda }}\left( {{H^3}\frac{{\partial P}}{{\partial \lambda }}} \right) = \frac{{\partial H}}{{\partial \varphi }} + 2\frac{{\partial P}}{{\partial \tau }} \end{array} $

(1)

1.3 油膜的动刚度和阻尼

轴承油膜的动刚度是主轴抵抗位移挠动的能力。轴承油膜刚度定义为:轴心有一微小挠动时，油膜力增加量对微小挠动位移的比值，用k_ij表示，第1个i表示载荷增量的方向，第2个下标j表示位移增量的方向，即：

$ {k_{ij}} = {\left( {\frac{{\partial {F_i}}}{{\partial j}}} \right)_0} $

(2)

式中, 0表示转子平衡位置。

轴承油膜阻尼用来表征油膜抵抗转子偏离原平衡位置时产生的速度扰动的能力。这样4个油膜阻尼系数用c_ij表示，第1个i表示载荷增量的方向，第2个下标j表示速度增量的方向，即：

$ {c_{ij}} = {\left( {\frac{{\partial {F_i}}}{{\partial j}}} \right)_0} $

(3)

这4个动刚度和4个阻尼系数即为油膜的8个动力特征系数。

所以油膜力可以表示为:

$ {F_x} = {({F_x})_0} + {k_{xx}}\Delta x + {k_{xy}}\Delta y + {c_{xx}}\dot x + {c_{xy}}\dot y $

(4)

$ {F_y} = {({F_y})_0} + {k_{yx}}\Delta x + {k_{yy}}\Delta y + {c_{yx}}\dot x + {c_{yy}}\dot y $

(5)

动刚度当量 ${K_0} = 6\eta L{w_0}\left( {\frac{R}{{{h_0}}}} \right)$ ，阻尼当量 ${C_0} = 6\eta L{\left( {\frac{R}{{{h_0}}}} \right)^3}$ ，用参数K₀和C₀动刚度系数和阻尼系数进行无量纲化。式(2)和(3)中的下标ij分别代表xx、xy、yx和yy。

2 液体静压轴承动力学特征求解方法 2.1 载荷增量法计算油膜动刚度 2.1.1 动刚度系数K_xx和K_yx计算

如图 2所示，O为轴承的中心位置，O^′为转子的平衡位置，转子沿x方向分别产生微小的位移扰动Δx后的位置为O₁和O₂。

转子在平衡位置O^′沿x方向向左产生微小扰动后，轴心位置移动到O₁处，此时，x和y方向的油膜力分别为F_x1和F_y1；转子在平衡位置O^′沿x方向向右产生微小扰动后，轴心位置移动到O₂处，此时，x和y方向的油膜力分别为F_x2和F_y2。根据油膜刚度的计算公式，此时油膜动刚度可以表示为：

$ {K_{xx}} = \frac{{{F_{x1}} - {F_{x2}}}}{{2\Delta x}} $

(6)

$ {K_{yx}} = \frac{{{F_{y1}} - {F_{y2}}}}{{2\Delta x}} $

(7)

2.1.2 动刚度系数K_xy和K_yy计算

如图 2所示，O为轴承的中心位置，O^′为转子的平衡位置，转子沿y方向分别产生微小的位移扰动Δy后的位置为O₃和O₄。

转子在平衡位置O^′沿y方向向下产生微小扰动后，轴心位置移动到O₃处，此时，x方向和y方向的油膜力分别为F_x3和F_y3；转子在平衡位置O^′沿x方向向上产生微小扰动后，轴心位置移动到O₄处，此时，x方向和y方向的油膜力分别为F_x4和F_y4。根据油膜刚度的计算公式，此时油膜动刚度可以表示为：

$ {K_{xy}} = \frac{{{F_{x3}} - {F_{x4}}}}{{2\Delta y}} $

(8)

$ {K_{yy}} = \frac{{{F_{y3}} - {F_{y4}}}}{{2\Delta y}} $

(9)

2.2 扰动压力法计算油膜阻尼

如图 3所示，转子在外载荷的作用下平衡于O^′位置，在(Δe, Δθ)的扰动下，轴心位置移动到O₁处，此时油膜力沿Δe的方向和垂直于Δe的方向分别为F_e和F_θ，

$ {F_e} = \int_{ - 1/2L}^{1/2L} {\int_0^{2\pi } {p{\rm{cos}}\;\varphi R{\rm{d}}\varphi {\rm{d}}z} } $

(10)

$ {F_\theta } = \int_{ - 1/2L}^{1/2L} {\int_0^{2\pi } {p{\rm{sin}}\;\varphi R{\rm{d}}\varphi {\rm{d}}z} } $

(11)

结合式(3)，将式(10)和(11)分别对 ${\dot e}$ 和 ${\dot \theta }$ 求导，并结合式(4)可得：

$ {c_{ij}} = - \int_{ - L/2}^{L/2} {\int_0^{2\pi } {{p_j}{\rm{cos}}\;\varphi R{\rm{d}}\varphi {\rm{d}}z} } $

(12)

式(5)和(12)中i、j分别代表ee、eθ、θe和θθ。

将沿Δe的方向和垂直于Δe的方向的阻尼经过坐标变换，变换到(x, y)方向：

$ \begin{array}{l} \left( {\begin{array}{*{20}{c}} {{c_{xx}}}&{{c_{xy}}}\\ {{c_{yx}}}&{{c_{yy}}} \end{array}} \right) = \\ \left( {\begin{array}{*{20}{c}} {{\rm{sin}}\;\theta }&{{\rm{cos}}\;\theta }\\ {{\rm{cos}}\;\theta }&{{\rm{ - sin}}\;\theta } \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{c_{ee}}}&{{c_{e\theta }}}\\ {{c_{\theta e}}}&{{c_{\theta \theta }}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{\rm{sin}}\;\theta }&{{\rm{cos}}\;\theta }\\ {{\rm{cos}}\;\theta }&{ - {\rm{sin}}\;\theta } \end{array}} \right) \end{array} $

(13)

2.3 液体静压主轴系统动力学模型

图 4为动载荷下转子受力分析。将计算所得的油膜压力用于转子运动行为的描述，实现液体静压主轴系统的动力学分析。

假设转子受到单一的偏心质量动载荷影响，e_g为转子的质量偏心距，所以作用在转子上的不平衡载荷为：

$ {Q_x} = M{e_{\rm{g}}}w_0^2{\rm{sin}}\;({w_0}t) $

(14)

$ {Q_y} = M{e_{\rm{g}}}w_0^2{\rm{cos}}\;({w_0}t) $

(15)

根据式(6)和(7)将油膜力F_x(w₀t)和F_y(w₀t)用油膜的动刚度和阻尼系数，所以转子的运动方程的无量纲形式为：

$ {M_0}\ddot X + {K_{xx}}X + {K_{xy}}Y + {C_{xx}}\dot X + {C_{xy}}\dot Y = {M_0}{\varepsilon _{\rm{g}}}{\rm{sin}}\;\tau $

(16)

$ {M_0}\ddot Y + {K_{yx}}X + {K_{yy}}Y + {C_{yx}}\dot X + \sqrt 2 {C_{yy}}\dot Y = {M_0}{\varepsilon _{\rm{g}}}{\rm{cos}}\;\tau + {M_{\rm{g}}} $

(17)

式中，当量质量 ${M_0} = \frac{{M{w_0}}}{{3\eta L}}{\left( {\frac{{{h_0}}}{R}} \right)^3}$ ，无量纲偏心质量距 ${\varepsilon _{\rm{g}}} = \frac{{{e_{\rm{g}}}}}{{{h_0}}}$ ，转子当量重量 ${{\bar M}_{\rm{g}}} = \frac{{{M_{\rm{g}}}}}{{3\eta LR{w_0}}}{\left( {\frac{{{h_0}}}{R}} \right)^2}$ 。

3 扰动压力分布求解 3.1 扰动压力求解程序流程

无论是油膜动刚度和阻尼的求解，还是轴心轨迹的求解都是以扰动压力分布的求解为基础的，而扰动压力的获得就是对方程(1)的求解，有限差分法是求解雷诺方程常用方法之一，唐倩等^[18]为研究滑动轴承合金层应力分布对雷诺方程采用了有限差分方法求解。采用有限差分法对扰动压力方程(1)进行求解，方程求解的程序流程：首先，开始输入初始参数，计算间隙及相关参数，设定方程扰动位移。其次，进行扰动压力场迭代并判断检验精度δ₁，若精度符合要去则进入下一步，否则返回继续迭代。然后，检验流量精度δ₂，若符合要求进入下一步，否则返回上一步继续迭代。最后，计算动力学相关性能，结束。

3.2 油垫网格划分

对第k油腔及其封油边进行网格划分(如图 5所示)，令封油边的边界均在网格的节点上，以利于边界条件的确定。用各节点的压力值构成各阶差商，近似取代Reynolds方程中的导数。

静压轴节点位置油膜压力值用P_{i, j}表示。图 5中, i为节点在周向方向的行数，j为节点在轴向方向的列数。设在周向方向上有m个格子，因此每个格子的宽度为 $\Delta \varphi = \frac{{4{\theta _2} - 2{\theta _1}}}{m}$ ，其中, θ₁、θ₂分别为油腔半角和油垫半角，在轴向方向上有n个格子，因此每个格子的长度为 $\Delta \lambda = \frac{2}{n}$ 。

3.3 雷诺方程离散化

利用网格划分各节点周围节点的差商来表示节点(i, j)上油膜压力P的1阶偏导数:

$ \left\{ \begin{array}{l} {\left( {\frac{{\partial P}}{{\partial \varphi }}} \right)_{i,j}} = \frac{{{P_{i + 1,j}} - {P_{i - 1,j}}}}{{2\Delta \varphi }},\\ {\left( {\frac{{\partial P}}{{\partial \lambda }}} \right)_{i,j}} = \frac{{{P_{i,j + 1}} - {P_{i,j - 1}}}}{{2\Delta \lambda }} \end{array} \right. $

(18)

任意节点(i, j)上无量纲油膜厚度的差分形式如下：

$ {\left( {\frac{{\partial H}}{{\partial \varphi }}} \right)_{i,j}} = \frac{{{H_{i + 1,j}} - {H_{i - 1,j}}}}{{2\Delta \varphi }} $

(19)

而无量纲油膜厚度H对τ的偏导数与转子沿x和y方向的运动速度有关，即：

$ \frac{{\partial H}}{{\partial \tau }} = \dot X{\rm{sin}}\;({\varphi _{i,j}}) + \dot Y{\rm{cos}}\;({\varphi _{i,j}}) $

(20)

将式(18)、(19)和(20)代入式(1)得到无量纲Reynolds方程的差分形式为：

$ \begin{array}{l} {P_{i,j}} = \frac{1}{{2({A_{i,j}} + {C_{i,j}})}}[({A_{i,j}} + {B_{i,j}}){P_{i + 1,j}} + \\ \;\;\;\;\;\;\left( {{A_{i,j}} - {B_{i,j}}} \right){P_{i - 1,j}} + {C_{i,j}}{P_{i,j + 1}} + {C_{i,j}}{P_{i,j - 1}} - {F_{i,j}}] \end{array} $

(21)

式中， ${A_{i,j}} = \frac{{\left( {{H_{i,j}} + 6B} \right)H_{i,j}^2}}{{\Delta {\varphi ^2}}}$ , ${B_{i,j}} = \frac{{3\left( {1 + 4B} \right){H_{i,j}}\left( {{H_{i + 1,j}} - {H_{i - 1,j}}} \right)}}{{4\Delta {\varphi ^2}}}$ , ${C_{i,j}} = \frac{{{\beta ^2}{H_{i,j}}^2\left( {1 + 6B} \right)}}{{\Delta {\lambda ^2}}}$ , ${D_{i,j}} = \frac{{{H_{i + 1,j}} - {H_{i - 1,j}}}}{{2\Delta \varphi }}$ , ${E_{i,j}} = 2(\dot X{\rm{sin}}\;({\varphi _{i,j}}) + \dot Y{\rm{cos}}\;({\varphi _{i,j}}))$ , ${F_{i,j}} = {D_{i,j}} + {E_{i,j}},\beta = \frac{{2R}}{L}$ 。

3.4 边界条件

对于静压轴承，需要通过流量平衡方程确定出节流后油腔压力，并以此作为求解Reynolds方程的一个边界条件。流入油腔的流量即为通过节流器流入油腔的流量，采用小孔节流器作为节流元件。根据流体力学的公式，计算流入节流器流量。(Q_i)_k为通过节流器流入k油腔的流量(k=1, 2, 3, 4)，流出油腔的流量(Q_o)_k可以通过对流速的积分求得。流入油腔的流量和流出油腔的流量相等(Q_i)_k=(Q_o)_k建立静压轴承的流量平衡方程。通过流量平衡方程能够计算出油腔压力(P_r)_k(k=1, 2, 3, 4)。同时，在计算封油边上的压力时，以油腔外的有关节点压力P_{i, j}作为边界条件，如下式：

$ {P_{i,j}} = {({P_{\rm{r}}})_k}。$

式中，i=i₁, i₁+1, …，i₂; j=j₁, j₁+1, …，j₂。

3.5 压力收敛准则

采用松弛法对压力场进行迭代计算，在每次迭代完成之后都要判断压力迭代结果是否已经达到足够的精度，从而决定是否终止迭代过程，可以根据收敛准则来判别。

$ \frac{{\sum\limits_{i = 2}^{m - 1} {\sum\limits_{j = 2}^{n - 1} {\left| {P_{i,j}^{\left( {k + 1} \right)} - P_{i,j}^{\left( k \right)}} \right|} } }}{{\sum\limits_{i = 2}^{m - 1} {\sum\limits_{j = 2}^{n - 1} {\left| {P_{i,j}^{\left( {k + 1} \right)}} \right|} } }} \le {\delta _1} $

(22)

式中，误差δ₁是一个非常小的量，其数量级一般为10^-3[20]时可以获得较高的计算精度。

3.6 流量收敛准则

在边界条件中计算流量平衡方程时，计算得到的流入油腔的流量(Q_i)_k与通过积分方法计算得到的流出油腔的流量(Q_o)_k不可能完全相等，当两者之间的误差达到一定程度时，即可认为流量是平衡的。

该误差通过δ₂来衡量，取δ₂=10^-3。当程序循环过程中流量平衡方程不能满足时需对每个油腔内压力(P_r)_k^(k+1)进行修正：

$ \left( {{P_{\rm{r}}}} \right)_k^{\left( {k + 1} \right)} = \left( {1 - {\omega _1}} \right)\left( {{P_{\rm{r}}}} \right)_k^{\left( {k + 1} \right)} + {\omega _1}\left( {{P_{\rm{s}}} - \frac{{8{{\left( {{Q_{\rm{o}}}} \right)}_k}}}{{{\pi ^2}d_0^2{\alpha ^2}}}} \right) $

(23)

式中，ω₁为流量修正因子，k=1, 2, 3, 4。

4 油膜动态特性数值仿真 4.1 油膜动刚度仿真

表 1为油膜4个动刚度系数在发生滑移和无滑移的情况下随偏心率的变化关系，b为滑移长度，K_xx、K_xy、K_yx和K_yy为油膜的4个无量纲刚度系数。

从表 1可以看出，无论滑移是否发生，4个刚度系数都随偏心率的增大而减小。速度滑移没有对油膜的4个动刚度系数与偏心率变化关系产生显著的影响，但在一定程度上降低了油膜的动刚度。动刚度系数K_xx、K_xy、K_yx和K_yy降低的最大比率分别为4.85%、4.85%、7.20%、7.16%，并且都发生在偏心率为0.8时。

表 2为在滑移长度为0和7 μm时，4个油膜阻尼系数和转子偏心率之间的变化关系。b为滑移长度，C_xx、C_xy、C_yx和C_yy为从油膜的4个无量纲刚度系数。从表 2中可以看出，轴承内固-液界面速度滑移并没有对4个阻尼系数与转子偏心率之间的变化关系产生显著的影响，速度滑移效应对阻尼系数C_yx、C_xy和C_yy的影响规律比较明显和单一；受速度滑移的影响，在不同的偏心率下，油膜的阻尼系数都有降低，速度滑移程度越大，油膜阻尼的减小程度越明显。阻尼系数C_xx受速度滑移影响比较复杂，在不同的偏心率下，表现为不同的影响规律。

4.2 不平衡载荷下轴心轨迹仿真

图 6为不同转子偏心质量距时转子旋转时轴心的行驶轨迹。

转子偏心质量距分别取为0.1、0.2、0.3、0.4，从图 6中可以看出:随着转子偏心质量距的不断增加，转子轴心的行驶轨迹范围不断扩大; 当偏心质量距为0.4时转子轴心轨迹与轴承间隙圆发生接触，此时转子与轴瓦发生接触，这将严重影响转子旋转的稳定性，该情况也是不允许出现的。

5 轴心轨迹实验检测

将标准球安放在液压回转轴系的工作台上，测微仪安放在轴系底座上，调整标准球与回转轴系的偏心，达到电感仪测量范围之内，稳定后开始检测，如图 7所示。

采样点512，采样间隔10 ms，得到主轴回转误差曲线，如图 8所示。

从图 8中可以看出，最大回转误差为2.6 μm，这说明轴心轨迹坐标小于图 6中偏心率为0.1时坐标，也即主轴系统的回转精度基本上没有受到不平衡量的影响。

6 结论

1) 应用有限差分法对扰动雷诺方程进行求解，解出油膜的扰动压力分布，通过载荷增量法和扰动压力法对油膜的4个动刚度系数和阻尼系数求解。在无速度滑移条件下4个动刚度系数都随偏心率的增加而减小；无滑移时，阻尼系数C_xx随偏心率的增加而减小，而阻尼系数C_xy、C_yx和C_yy会随偏心率的增加而增大。

2) 将固-液界面上的速度滑移考虑到液体静压轴承油膜性能的研究中，速度滑移并没有对4个动刚度系数和4个阻尼系数与偏心率的变化关系产生显著的影响，但在一定程度上降低了油膜的动刚度，动刚度系数K_xx、K_yx、K_xy和K_yy降低的最大比率分别为4.85%、4.85%、7.20%、7.16%，并且都发生在偏心率为0.8时；受速度滑移影响，在不同的偏心率下，油膜阻尼系数C_yx、C_xy和C_yy都有所降低。阻尼系数C_xx受滑移影响比较复杂，在不同的偏心率下，表现为不同的影响规律。

3) 将4个动刚度系数和4个阻尼系数用于建立转子的轴心轨迹模型，分析不同偏心质量距下转子轴心的运行轨迹。随着转子偏心质量距的不断增加，转子轴心的行驶轨迹范围也在不断扩大，当偏心质量距为0.4时转子轴心位移最大，此时系统也最不稳定。