透平式能量回收一体机水润滑转子稳定性研究

宁延强，李 岩，张德胜，叶晓琰，丁竞飞

（江苏大学 国家水泵及系统工程技术研究中心，江苏镇江 212013）

0 引言

海水淡化高压泵与能量回收一体机［1-2］是膜法海水淡化系统节约能源、降低成本的核心部件。通过将泵叶轮与透平叶轮同轴连接，泵端提供超高压海水，为海水通过半透膜提供必要条件，未通过半透膜的高压浓盐水则流入透平端，推动透平叶轮旋转回收其压力能。一体机主轴采用水润滑轴承支撑，相比油润滑更为环保、经济，但由于水的黏度较低，因此对设计要求则更高。由于一体机要在高转速下运行［3］，而水润滑轴承-转子系统［4-6］作为一体机的核心部件，对保障一体机的平稳运行起着至关重要的作用。

一体机的工作转速高达20 000 r/min，且水润滑形成条件苛刻，因此一体机中水润滑转子系统的稳定性分析［7-9］不可忽视。而对一体机中水润滑转子系统的研究仅停留在由单个水润滑滑动轴承支撑，而未考虑口环间隙中水膜产生流体动压润滑效应对一体机内部转子稳定性及临界转速的影响。目前对转子的稳定性判断有多种方式，主要包括实验及数值计算。侯峰等［10］对双螺杆制冷压缩机转子轴心轨迹波动特性进行了试验研究；叶晓琰等［11］研究了不同半径间隙下水润滑轴承轴心轨迹的变化；ZHAO 等［12］对大长径比下的水润滑橡胶轴承进行实验研究并确定其动力特性［13］；GAO 等［14］用CFD［15］软件对水润滑径向滑动轴承流体动力润滑承载能力进行分析，并提出了一种新型过渡弧轴承轴瓦，有利于提高水动力承载能力；DU 等［16］研究了不可压缩层流静压轴承的动力学特性，将槽间流动连续方程从雷诺方程解耦，提高了计算效率；谢帆等［17］基于有限差分法计算了径向滑动轴承油膜压力分布；李超等［18］对考虑间隙影响的涡旋压缩机转子系统的动态特性进行了研究；LI 等［19］采用窄轴承理论求解滑动轴承液膜的非定常雷诺方程并计算了动力学特征系数；成金贵等［20］运用平均雷诺方程和固体接触理论建立了止推轴承接触面间摩擦润滑理论模型并进行了试验验证；LI 等［21］推导了考虑五自由度运动时判别转子稳定性的临界质量快速计算公式。

综上可知，现有研究均未考虑口环间隙水膜的动压润滑效应，且尚无证据表明该效应可以忽略。本文通过求解Reynolds 方程和转子动力学特征方程，考虑口环间隙水膜动压润滑效应的影响，研究了一体机转子系统的临界转速及系统的临界质量，为保障海水淡化高压泵与透平式能量回收一体机水润滑转子系统的稳定运行提供参考。

1 模型及参数

本文研究对象为海水淡化高压泵与透平式能量回收一体机中水润滑转子系统，其结构如图1所示。其转子结构左端为透平叶轮，右端为泵叶轮，中间部分为水润滑轴承，泵叶轮进口延伸段处的轴承为口环。一体机中转子部件的主要结构参数见表1。

表1 水润滑转子系统结构参数Tab.1 Structural parameters of water lubricated rotor system

图1 一体机水润滑转子系统结构示意Fig.1 Structure diagram of water lubricated rotor system of all-in-one machine

2 水润滑轴承动静特性分析

在一体机中，转子系统的动力学特性与轴承的支撑特性密切相关，水润滑滑动轴承及口环轴承的动力特性计算与分析是整个转子系统的核心。

2.1 水润滑轴承静特性计算

水润滑轴承原理如图2 所示，A点为沿y轴逆时针转动φ后的位置，ϕ为偏位角，θ=φ-ϕ，c为半径间隙，则水膜厚度h=c+ecosθ，e为偏心距。

图2 水润滑轴承原理Fig.2 Schematic diagram of water-lubricated bearings

图3 水膜网格划分Fig.3 Water film meshing

水润滑滑动轴承及口环轴承均属于圆柱形流体动压轴承，通过求解雷诺方程得到液膜压力分布。无量纲化后的雷诺方程为：

式中，L为轴承长度；P为水膜压力；R为轴颈半径。

利用有限差分法计算水膜压力，需对水膜进行正交网格划分，将水膜划分成m×n的网格，如图（i，j）处节点的压力，且ΓΔPλ=2/n，Δθ=（θH2-θ1）/m，3 所示，每一个节点有对应的坐标，Pi，j表示坐标为θ1为水膜起始角，θ2为水膜终止角。计算各节点的压力求得压力分布，将压力分布以水膜表面进行积分求得水膜承载力。

为了提高计算精度，采用半步长差分法，对水膜端面采用第一类边界条件，即λ=±1，P=0；水膜周向上采用雷诺边界条件，即起始边固定，则有θ=0，P=0；水膜终止边自由浮动，并取P=0，∂P/∂θ=0。采用超松弛迭代法计算各个节点的压力值，初始化各个节点压力值为P=0，代入雷诺边界条件，求解得到压力分布。

为求解水润滑轴承动力特性系数，需要确定轴承偏心率及偏位角并求得水膜承载力。将水膜承载力分解为x方向和y方向，分别用Fx，Fy表示，则有：

式中，Fx，Fy为无量纲水膜承载力；θ为水膜旋转角，以x轴负方向为水膜角初始位置。

2.2 水润滑轴承动特性计算

计算水润滑轴承刚度阻尼动特性系数采用瞬态雷诺方程求解，对其无量纲化后可表示为［22-23］：

式中，ε'，εϕ'为瞬时速度。

将式（4）依次对ε，ϕ，ε'，εϕ'求偏导得到扰动微分方程为：

由静平衡位置及水膜压力分布结合式（5）～（8）解得水膜扰动压力，进而求解水膜的刚度阻尼系数：

由上式解得滑动轴承及口环轴承动特性系数。

将式（9）（10）化为直角坐标形式：

3 水润滑转子系统稳定性分析

在一体机中，当转子受到扰动时能否保证系统平稳运行与转子的稳定性密切相关。水润滑转子系统的稳定性可由求解转子动力学特征方程，通过其特征值实部判别，并通过求解转子系统临界质量进行稳定性分析。

将一体机中转子简化为一维模型，由受力平衡、力矩平衡方程求得滑动轴承及口环轴承处所受到的压力，受力分析如图4 所示。其中G为轴质量，G1，G2分别为透平叶轮质量、泵叶轮质量，F1，F2分别为滑动轴承、口环轴承所提供的承载力。

图4 转子受力分析Fig.4 Rotor stress analysis

列出平衡方程为：

求解得F1，F2后，进而转化为水膜承载力，分别解得滑动轴承刚度、阻尼矩阵K1，C1及口环轴承的刚度、阻尼矩阵K2，C2：

一体机转子系统的动力学方程为：

式中，X为转子的位置向量；M为质量矩阵；C为阻尼矩阵；K为刚度矩阵；F为不平衡力。取F=0，则式（16）变为转子自由振动模型，改写为：

因此可得到转子系统的特征方程为：

式中，ξ为特征方程的根，为转子质量。

将式（18）展开可得：

式（21）的解即为系统的特征值，求得的特征值为两对共轭的复根，即ξ=αi+im βi，其中αi（i=1，2，3，4）为特征值实部，im为虚数单位，β i为虚部。当α i＜0 时，系统为稳定运行状态；当α i=0 时，系统处于稳定运行的临界状态；当α i>0 时，系统处于不稳定状态，并求得系统临界质量。

综上所述，求解水润滑转子系统的稳定性流程如图5 所示。

图5 水润滑转子系统稳定性程序框图Fig.5 Water lubricated rotor system stability program block diagram

4 验证

结合文献［24］中模型参数，选取与文献相同的计算模型参数，轴承直径为30 cm，间隙比为0.002，宽径比为0.8，转速为3 000 r/min，载荷为65 000 N。图6 示出MATLAB 计算结果与文献中计算结果对比，其中主刚度系数Kxx，Kyy误差为4%～4.4%，交叉刚度系数Kxy，Kyx误差为0.4%～4.3%，主阻尼系数Cxx，Cyy误差为3.3%～4.6%，交叉阻尼系数Cxy，Cyx误差为4.1%～4.2%。考虑到计算过程中迭代算法和收敛准则等差异，会导致存在一定的结果偏差，计算结果与文献中的结果基本吻合，验证了本文计算结果的可靠性。

图6 MATLAB 计算结果与文献结果对比Fig.6 Comparison of MATLAB calculation results and literature results

5 结果分析

考虑口环轴承前、后，轴承水膜压力分布如图7 所示，对比可知滑动轴承处水膜压力约减小了40%。当转子系统考虑口环间隙水膜所产生的流体动压润滑效应后，口环轴承会承担一部分支撑作用。由于在一定的载荷下，相应的水膜提供的承载力也是一定的，因此滑动轴承处水膜压力的减小量近似于口环轴承水膜压力的增量。

图7 水膜压力分布Fig.7 Water film pressure distribution

当一体机转子系统考虑口环轴承后，滑动轴承及口环轴承的刚度系数随着转速变化规律如图8 所示。两轴承的主刚度系数各自近似相等，即K1xx=K1yy，K2xx=K2yy，随着转速增加，滑动轴承主刚度呈明显线性递增的趋势，而口环轴承仅小幅增加。滑动轴承与口环轴承的交叉刚度系数均呈x轴对称分布，即K1xy=-K1yx，K2xy=-K2yx，且滑动轴承的交叉刚度约为口环轴承的5 倍。随着转速增加，两轴承的交叉刚度系数均增大。

图8 刚度系数与转速的关系Fig.8 The relationship between stiffness coefficient and rotational speed

当一体机转子系统考虑口环轴承后，滑动轴承及口环轴承的阻尼系数变化曲线如图9 所示。随着转速的增加，滑动轴承及口环轴承主阻尼系数无明显变化。对于交叉刚度系数而言，当转速＜15 000 r/min 时，系数变化较为明显且呈减小趋势；当转速＞15 000 r/min 时，系数趋于稳定不再变化，且C1xy=C1yx，C2xy=C2yx。

图9 阻尼系数与转速的关系Fig.9 The relation between damping coefficient and rotational speed

考虑口环轴承前、后的动力特性系数分别为A，B，则动力特性系数变化率表示为：（B-A）/A×100%，则其变化率如图10 所示。考虑口环轴承后，系统动力特性系数均有所增加，系统交叉刚度系数及主阻尼系数增加量超过20%，主刚度系数和交叉阻尼系数增加5%左右。因此一体机转子系统中口环轴承的存在有利于水润滑轴承承载特性的提升。

图10 动特性系数变化率Fig.10 Rate of change of dynamic characteristic coefficient

添加口环轴承后特征值变化曲线如图11 所示，临界质量变化曲线如图12 所示。

图11 特征值变化曲线Fig.11 Eigenvalue change curve

特征方程的解为2 个共轭负根，且有2 个实部，随着转速的升高实部1 呈现先减小后增大的趋势，而实部2 单调递减且恒为负值。考虑口环轴承后，实部1 有所增加但上升趋势减缓，故而临界稳定状态时的转速增加。随着转速增加，系统的临界质量逐渐减小，在一体机工作转速下，考虑口环轴承支撑后系统的临界质量将提升8.4%左右。因此未考虑口环轴承将导致系统临界转速和临界质量的计算产生误差，从而影响水润滑转子系统稳定性设计的准确性。

不同长径比下水润滑转子系统临界质量变化如图13 所示，长径比的增加使得系统的临界质量增加，且在低转速下临界质量受长径比的影响较为明显，这是由于轴承的刚度增加，承载性能有所提升，增大轴承的长径比有利于提高系统稳定性。

图13 不同长径比的临界质量Fig.13 Critical mass of different length diameter ratio

不同半径间隙下的水润滑转子系统临界质量变化如图14 所示，随着轴承半径间隙的增加，系统的临界质量逐渐降低，轴承的承载能力下降导致系统在高转速下运行不稳定，因此采用减小半径间隙的方法来提升系统的临界质量，且临界质量随半径间隙的降低呈现非线性增加，半径间隙越小，临界质量的增量越大。

图14 不同半径间隙的临界质量Fig.14 Critical mass of gap with different radius

6 结论

（1）当未考虑口环间隙水膜时，滑动轴承水膜将会承受更高的压力；而考虑口环轴承则降低了滑动轴承的偏心率，有利于提高系统的稳定性。

（2）当未考虑口环轴承时，水润滑转子系统动特性系数的求解准确性将会降低；考虑口环轴承后，其动特性系数均有所增加，其中交叉刚度系数及主阻尼系数增幅超过20%，主刚度系数与交叉阻尼系数增加约为5%，且随转速的增加趋于稳定。

（3）当未考虑口环轴承时系统的临界转速和临界质量的求解结果存在一定偏差，均比考虑口环轴承时有所减小，其中系统的临界转速减小7%，临界质量减小8.4%。因此，考虑口环间隙水膜的动压润滑效应可提高一体机水润滑转子系统的稳定性。

（4）适当增加轴承的长径比、减小半径间隙可以有效提高其承载性能和临界质量，有助于增强其转子系统的稳定性。