动压箔片轴承气膜厚度分布特点

王文升，郑越青,2，徐刚，徐庆元

(1.中国工程物理研究院 机械制造工艺研究所，四川 绵阳 621900；2.西安交通大学 能源与动力工程学院，西安 710049)

动压箔片轴承是一种由平箔和波箔组合的空气动压轴承，它以平箔作为轴承表面，以波箔作为弹性支承[1]，采用空气润滑。由于在运行过程中采用气体润滑，箔片轴承具有摩擦力低、损耗小及寿命长等特点；箔片的弹性可以平衡转子的不平衡性，具有很强的环境适应性；选择性能优良的轴承材料和涂层材料的箔片轴承可以耐受高温。箔片轴承的这些特点使其在诸多性能上超越了传统支承[1-3]，可以作为高速和超高速转子系统的支承。

箔片轴承特性仿真涉及流体力学与弹性力学，并且存在二者的耦合。目前最常用的方法是简单弹性基础模型[2-3]，该模型将平箔与波箔的组合简化为弹性均匀的柔性支承体，并且认为柔性支承体的刚度完全取决于波箔。该模型假设支承体表面任一位置的变形仅与支承刚度及该处压强有关，而与周围结构变形无关，即柔性支承体表面各点是相对独立的。该假设显然不符合平箔是连续介质这一基本事实。

为解决该模型存在的缺陷，在简单弹性基础模型基础上，细化了对平箔与波箔的数学描述，采用Euler梁单元描述平箔，简单弹簧单元描述波箔。以此模型为基础，对箔片轴承的气膜厚度与压力分布进行了分析。

1 数学模型

箔片轴承结构示意图如图1所示，其由构成轴承表面的平箔和提供弹性的波箔组成。平箔的一端与波箔固定，另一端自由。正常工作时，转子沿自由端至固定端的方向旋转。

1.1 压力场控制方程

箔片轴承气隙中的气体压力场采用Reynolds方程描述，为减小舍入误差，Reynolds方程通常采用无量纲形式，即

(1)

式中：X为无量纲周向坐标；Y为无量纲轴向坐标；P,H分别为无量纲气膜压力和无量纲气膜厚度；Λ为轴承数。无量纲数与轴承数的定义可参考文献[2]。

(1)式为非线性偏微分方程，需采用Newton-Raphson方法进行求解。将连续的Reynolds方程以中心差分形式展开，得到各差分节点处的离散Reynolds方程(角标规则如图2所示)

图2 网格与节点角标规则

(2)

构造除边界点外的各节点上的离散Reynolds方程组，并采用Newton-Raphson迭代法进行求解，获得各节点压力值。(2)式中P值在首次计算时需给定一组初始值，然后以迭代过程中获得的P分布不断更新。

1.2 气膜厚度数学模型

气膜厚度等于初始气膜厚度与平箔变形量的叠加[4]

h0=c+ecos(θ-θ0)+u，

(3)

式中：c为轴承名义间隙；e为偏心量；θ为周向角度；θ0为轴承偏位角；u为平箔变形量。

当转子旋转时，空气在黏性力作用下随转子在气隙中运动，当进入收敛区时气体被压缩而使压力提高。气体压力作用于平箔表面，一方面波箔发生弹性变形，另一方面平箔本身也会发生变形。为描述以上现象，对平箔与波箔建立不同的有限元模型，平箔采用Euler梁单元描述，而波箔采用简单弹簧描述(图3)。

图3 波箔与平箔的有限元模型

平箔的Euler梁单元刚度矩阵为

(4)

式中：E为弹性模量；Iz为平箔截面矩；l为单元长度。

波箔的简单弹簧单元刚度矩阵为

(5)

式中：k为波箔的刚度系数[4]；b为轴承宽度；l0为波拱宽度的一半；ν为泊松比；t为波箔片厚度。

2 计算方法

2.1 边界条件

计算过程采用有限元与有限差分2种计算方法，因此计算过程中使用两类边界条件：一类是压力边界条件；另一类是位移边界条件。

气膜间隙区域的4个边界，即轴承两侧以及箔片固定端与自由端，均与大气接触，因此环境压力边界条件为

(6)

简单弹簧单元的固定端施加位移边界为

u(n+1：n+imax)=0。

(7)

2.2 收敛标准

求解量包括气膜压力P、气膜厚度H与偏位角θ0，它们需要满足的收敛标准为

(8)

2.3 迭代过程

整个计算由3层迭代过程组成，最底层是Newton-Raphson迭代，用于求解Reynolds方程；中间层是气弹耦合迭代，用于获得满足耦合条件的气膜压力分布与气膜厚度；最上层为偏位角修正迭代，用于校正转子姿态，以满足受力平衡要求。详细计算过程如图4所示。

图4 程序流程框图

给定初始的轴承名义间隙与偏心量，假定偏位角为0，且平箔没有变形，根据(3)式计算出初始气膜厚度；然后以该初始气膜厚度求解Reynolds方程，计算压力分布；由压力分布获得平箔受力情况，将其带入有限元模型计算平箔变形，从而获得新的气膜厚度分布；再以该气膜厚度计算压力分布，如此迭代，直至收敛。

获得收敛解后，得到稳态下的气膜压力，然后通过对x和y方向进行积分，得到轴承承载力以及偏位角。承载力与偏位角计算方法见文献[4]。由于初始假定偏位角为0，该解并不满足受力平衡要求，因此以计算获得的偏位角替代之前的偏位角，重复上述过程。持续迭代，直至2次求解获得的偏位角误差小于收敛标准。

3 结果与分析

箔片轴承的气膜压力与气膜厚度分布分别如图5和图6所示，其转速为30 000 r/min，偏心率为0.9，其他轴承和气体参数见表1。由图5和图6可知，随着气膜厚度的减小，气体被压缩使压力迅速增大。由于重力作用气膜厚度减小区域位于轴承下方，因此压力增大区域也位于轴承下方，从而形成向上的承载力。通过气膜厚度减小段后，气体进入气膜厚度增大段，气压逐渐恢复到大气压，在气膜厚度较大区域(轴承上部分)气膜压力基本等于大气压[5]。

表1 轴承和气体参数

图5 无量纲气膜厚度分布

图6 无量纲气膜压力分布

以上结果与简单弹性基础模型的计算结果基本一致[6-7]，不同之处在于该模型的高压区出现小尺度空间上的压力波动，这一现象是由波箔下陷变形引起的。从图5和图6可见，在气膜厚度较小处，也就是压力较大区域，平箔出现局部下陷，该分析结果与文献[8]中的试验结果基本一致。

试验后的平箔表面存在摩擦形成的类似斑马线状磨痕(图7)，这是由运转过程中转子与定子碰撞摩擦所致。由于平箔的下陷效应，下陷区域没有发生接触摩擦，从而使磨痕呈现出斑马线状的条纹。

图7 试验后箔片的斑马线磨痕

4 结束语

建立了考虑平箔片下垂效应的动压箔片轴承气弹耦合模型，通过迭代方法耦合求解了箔片轴承的气膜压力分布与气膜厚度分布。基于仿真结果，对箔片轴承的气膜厚度分布与压力分布特点进行了分析，为平箔下陷这一试验现象提供了合理解释。