(1.上海工程技术大学机械与汽车工程学院 上海 201620；2.上海交通大学机械系统与振动国家重点实验室，振动噪声研究所，船舰设备噪声与振动控制技术国防重点学科实验室 上海 200240)

液体动压轴承常用于高转速及高精度旋转机械，其可靠性会直接影响到整个机器的精度、使用寿命、经济性等。动压轴承一直是转子动力学和轴承润滑研究领域的热点问题，若轴承在设计、安装、维护等方面有失误，会出现振动噪声等,甚至造成重大损失。高转速下系统的稳定性已成为转子动力学分析的重要课题，其中的关键问题之一是轴承动态刚度的准确计算。

确定油膜的刚度就需要求解油膜力，而求解油膜力就需要求解油膜的压力分布，理论计算上通常采用不同形式的雷诺方程进行求解，很多学者专家在过去的几十年间进行了诸多研究。1965年，LUND[1]运用数值计算方法求解雷诺方程，得到了刚度系数，建立了 “轴承-转子系统”模型，从理论上解释了系统涡动失稳的形成原因。通常情况下，雷诺方程作为润滑方程用于计算动压轴承的刚度[3-5]，但对于雷诺方程求解不可避免地存在一系列的假设，与实际存在一定的偏差。雷诺方程很难准确反映转子高速转动过程中引起的周向惯性效应、动态挤压效应等的影响，油膜刚度计算精度受到影响，因此有必要直接以Navier-Stokes(N-S)方程为基础，精确研究主轴转速、位移扰动等对轴心轨迹、油膜刚度等的影响规律。文献[6]应用基于N-S方程的计算流体动力学软件FLUENT进行了动压轴承油膜力的计算研究工作，证实了用流体计算软件求解动压轴承油膜力的可行性；文献[7]进一步引入了动网格模型，计算了基于CFD 动网格模型的液体动静压轴承油膜刚度；文献[8]研究计算了动特性系数，但计算油膜力时没有考虑气穴效应；文献[9]应用FLUENT计算了气穴影响下的椭圆轴承油膜压力场，证明考虑气穴影响对椭圆轴承的模拟结果较不考虑气穴影响更为合理。本文作者以圆柱径向动轴承为例，在ADINA软件中建立流固耦合模型，考虑气穴效应，直接以Navier-Stokes 方程为基础，采用ALE法，有效避免了油膜的网格畸变的问题，计算了油膜瞬态压力场，进而得到轴心轨迹，动态油膜刚度，结合试验验证了该计算方法的有效性和可行性，为动态油膜刚度的计算提供了一种方法。

1 刚度计算方法

刚度计算方法是在ADINA软件基础上，对轴颈施加正弦位移小扰动，求解N-S方程得到每一个迭代步的轴颈位置变动前后流体的瞬态压力场和不同方向的瞬态油膜力；结合差分计算模型，计算轴承油膜动态刚度的一种方法。

1.1 刚度的差分计算模型

轴颈在平衡位置受外载荷作用，施加位移小扰动，油膜力会发生变化。在小扰动条件下，轴承油膜刚度可采用差分法求解[7]，如式(1)所示。

(1)

式中:ΔFdij为油膜力的变化；kxx和kyy是直接刚度；kxy和kyx为交叉刚度； Δx和Δy为扰动位移。

采用差分模型计算轴承的刚度需要考虑位移对不同方向油膜力的影响。由式(1)可知，计算油膜刚度系数，需要：①在轴的重力作用下，从初始时刻(轴与轴承同心)状态下，开始计算瞬态油膜压力场，该求解是基于气穴效应的瞬态流固耦合模型，获得轴心平衡位置，画出轴心轨迹；②在轴心平衡位置施加轴颈位移扰动，求解瞬态油膜力，得到每一迭代步位移变动前后瞬态油膜力。

1.2 瞬态油膜力控制方程

1.2.1 流体控制方程

在CFD软件中，油膜流动的基本动力润滑理论是N-S方程(Navier-Stokes Equations)，不考虑流体的黏温效应，动力润滑方程包括连续性方程和动量方程[10-12]:

(2)

(3)

式中：t是时间；ρf是流体密度；τf是流体区域的应力张量；v是速度矢量；ff是流体部分的体力矢量。

τf=(-p+λ·v)I+2μe

(4)

用任意的拉格朗日-欧拉(ALE)公式求解上述方程，从而求得流体区域的变形。

1.2.2 固体域控制方程

对固体域施加位移扰动时，固体域的控制方程可用式(5)描述。

(5)

1.2.3 流固耦合方程

在流固耦合界面上，流体和固体将同时满足位移和力平衡的条件。位移协调方程和力平衡方程分别如式(6)和(7)所示。

df=ds

(6)

n·τf=n·τs

(7)

式中：df和ds是流体和固体在流体-结构界面上的节点的位移;τf和τs是流体和固体在流体-结构界面上的节点的应力；n是方向向量。

以上参数仅在流固耦合面上。流固耦合方程采用直接法求解，耦合面上满足力和位移平衡条件，当力和位移的相对残差小于0.01时即认为该迭代步收敛。

1.2.4 气穴模型

不采用气穴模型时，软件模拟时油膜压力分布中存在较大负压，会导致油膜力计算偏差较大[9]。文中流体计算采用了相变模型。流体和气体相变通过求流体域压力函数的导数来实现。SUN和BREWE[13]曾说明这种相变模型能够真实反映薄油膜动态情况下气相的变化情况。计算中，润滑油看成不可压缩的流体，气相看成可以轻微压缩的流体。气相的体积分数f通过气体的体积与总体积的比来确定。m定义为任意情况下流体-气体的体积比，由下式得到：

m=ml+f(mv-ml)

(8)

式中：下标l表示流体，v表示气体。

当流体压力p大于定义的蒸发压力plv时，气相分数f=0；当流体压力p小于定义的冷凝压力pvl时，气相分数f=1。用如下函数近似定义气相体积分数的变化(式中α=0.1)：

(9)

(10)

(11)

冷凝压力(pvl)的取值通过敏感性分析后确定。假设轴为刚体，计算冷凝压力为-20、-60、-80 kPa时轴心位置；在取该压力值时，轴心位置的偏差小于1%。文中计算，蒸发压力(plv)和冷凝压力(pvl)分别取0和-60 kPa。

2 典型轴承的油膜刚度计算

2.1 轴承结构和建模

按试验台中径向圆柱动压轴承(如图1所示)尺寸建模，主轴轴颈直径是40 mm，轴承的直径是40.08 mm，轴承轴瓦的长度是32 mm。该轴承承受主轴重量为150 N，在ADINA中完成轴和油膜的流固耦合模型[14]，采用六面体网格划分模型，油膜周向划分320份,轴向划分30份，厚度方向分5层，共48 000个单元；轴的周向划分320份,轴向划分10份，径向划分12份。轴和油膜网格模型如图2所示。假设流体不可压缩，等温、黏度不变，层流。采用式(8)—(11)所表示的气穴模型，油膜与轴接触面施加切向速度，油膜与轴瓦接触的表面定义为固定壁面，轴上施加重力150 N，轴与流体的接触面定义为流固耦合面，采用瞬态分析。

图1 径向轴承结构

图2 三维模型网格图

2.2 轴心轨迹

求解时，初始位置轴与轴衬同心，油膜各个方向厚度相同，然后在轴的重力作用和轴颈转速共同作用下移动，油膜厚度发生改变，油膜产生油膜力，轴心位置不断改变。在Δλ/ΔY≤2%时(如图3所示)，认为轴到达平衡位置。图3中，ΔY为任意周期内点A、B坐标的平均值的绝对值(水平横线)。Δλ为点A或点B坐标到该横线距离中的较小值。转速为2 500 r/min时轴心轨迹如图4所示，1.2 s时达到平衡位置。

图3 轴心位置判断方法

图4 轴心轨迹

2.3 瞬态压力场

油膜的压力分布随轴心轨迹的变化而变化，初始时刻，油膜厚度各个方向相同，无油膜压力，在0.001 s时(图4中点C)油膜压力场如图5(a)所示, 油膜压力分布最大值为0.3 MPa，随着油膜压力场的变化，轴心上浮慢慢稳定下来，油膜压力分布最大值略有降低，迭代至1.2 s时(图4中点D)，油膜压力分布最大值为0.27 MPa，如图5(b)所示。

图5 不同时刻油膜压力分布

2.4 动态刚度

通过计算轴心轨迹，确定轴心平衡位置，然后再对轴施加正弦扰动，得到各个迭代步的瞬态油膜力，通过式(1)，来确定动态油膜刚度。以轴颈转速为2 500 r/min为例，计算油膜4个动态刚度kxx、kxy、kyx、kyy，其数值的数量级多为6、7、8，如图6—9所示。计算时每间隔4.88×10-4s读取对应位移与力的数值，总共读取2 048对数据，以对应试验时的2 048个采样点，各刚度值的数量级个数统计如图10所示。

图6 动态油膜刚度kxx

图7 动态油膜刚度kyx

图8 动态油膜刚度kxy

图9 动态油膜刚度kyy

图10 刚度的数量级个数

3 试验研究

试验台示意图、结构图如图11、12所示，试验和数值计算主要是研究轴承2的油膜刚度，传感器安装于轴承上，传感器输出的信号相对于轴瓦是离散数值，对离散点进行拟合得到轴瓦上动态油膜压力分布情况[15]。

图11 试验台示意图

图12 试验系统布置图

试验测得下半瓦传感器的一阶频谱幅值比上半瓦的一阶频谱幅值要大7倍，也就是下半瓦上测得的压力值要比上半瓦的大，这与数值模拟的结果相符。测量刚度时，采用电涡流位移传感器在主轴运转中记录下振动位移情况，第一时刻，记录下油膜力和振动位移；第二时刻，记录下油膜力和振动位移，这些参数就可以计算动态油膜刚度。转速为2 500 r/min，信号是在时间1 s内进行了个2 048个点的采样，计算得到动态刚度系数如图13所示。kxx、kxy、kyx、kyy数值的数量级多为6、7、8，总体与数值计算结果相符。由于数值计算时假设润滑油不可压缩、油温不变，忽略了主轴不平衡等因素，而试验时可能存在转子质量动不平衡、轴承加工误差等问题，会影响油膜动刚度系数的准确性，因此仿真结果和试验之间局部各数据点并非完全一致，但仿真结果与以上实测动态刚度数值相比较发现油膜动态刚度的变化趋势是一致的，两者基本吻合，说明基于气穴模型的数值仿真计算方法是可行和有效的，为动态油膜刚度的计算提供了一种方法。

图13 动态油膜刚度系数

4 结论

(1)采用瞬态分析，基于流固耦合方法和气穴模型，轴在重力和油膜力的作用下从轴与轴衬同心开始计算，逐渐稳定在平衡位置，再现了轴心动态平衡过程，计算过程更加符合工程实际。

(2)数值仿真油膜动态刚度与实测动态刚度数值的变化趋势是一致的，两者基本吻合，说明基于ADINA软件的数值仿真计算方法是可行和有效的，为动态油膜刚度的计算提供了一种方法。

(3)计算得到的油膜力和油膜厚度的动态发展过程、轴心轨迹和动态刚度在轴承设计中考虑瞬态阶段的安全性方面提供了理论上的依据。