计算流体润滑频变动力学特性的一种完整变量扰动偏导数法

毕春晓， 韩东江， 杨金福

(1. 中国科学院 工程热物理研究所，北京 100190；2. 中国科学院大学，北京 100049)

高速旋转机械应用广泛，在能源动力、国防军工、航空航天、车辆、舰船等具有高性能、极端参数、长寿命要求的领域发挥着重要作用[1]。高效率和结构紧凑的特点使旋转机械在保证可靠性的同时向高转速发展，因此高速轴承和密封的性能分析、优化设计及其轴系稳定性控制已成为关键问题[2]。

使用工质流体作为润滑介质的轴承既能彻底避免工质被油污染、降低密封技术难度，同时缩短转子以提高涡轮机轴系的动力学性能[3]，是国内外相关领域的重要发展方向之一[4]，比如自润滑轴承用于低温工程[5]和有机朗肯循环(organic rankine cycle, ORC)的工质泵[6]，液体火箭发动机的高速涡轮泵[7-8]。1970年美国的空间Brayton循环项目中的发电装置最终采用He-Xe工质润滑箔片轴承[9]。超临界二氧化碳(S-CO2)闭式Brayton循环在太阳能光热发电[10]、船舶和飞行器余热利用[11-12]等领域有着广阔前景。目前S-CO2自润滑轴承已用于多个高速涡轮机样机[13-14]，其结构形式的研究和设计已引起国内外越来越多学者的重视[15-17]。

特殊参数下的流体通常具有复杂的热物性，从而带来润滑膜流动的多种附加效应，如表1所列举。具体到数学模型即润滑雷诺方程中不只含有膜厚和压力两个显函数变量，还包含热物性变量以及表示附加效应的非线性函数。

表1 特种参数流体润滑的附加效应Tab.1 Additional effects of fluid lubrication with special parameters

上述文献对于轴承或密封的润滑流动附加效应研究集中于稳态特性，然而这些复杂效应更加导致动力学特性的求解难度明显增大。Lund[26]提出计算滑动轴承刚度阻尼系数的摄动法并结合“瓦块装配法”求解可倾瓦轴承同步动力学系数，即轴颈和瓦块的扰动频率等于轴颈旋转的圆频率。虽然该方法在推导过程中通过归并消掉了位移和速度扰动量，从而天然满足临界特性与线性失稳分析所要求的平衡位置邻域内扰动[27]。然而大量试验结果表明，无论是油润滑可倾瓦轴承[28]还是气体轴承[29-30]，在相同的旋转频率(转速)下，动力学刚度和阻尼系数随扰动频率变化。文献[31]通过含扰动频率的复指数描述小扰动下的膜厚和压力动态变化，从而能够处理空气轴承雷诺方程中的时滞压力项，并揭示了可压缩润滑轴承动力学特性的频率效应。而且频率扰动方法能够纳入轴瓦与弹性支承的运动和动态变形引起的频率效应，实现箔片轴承或可倾瓦轴承的动态耦合分析[32-33]。比如Yang等[34]将偏导数法用于可倾瓦空气轴承，求解包含瓦块和支点扰动的全局动力学系数，并在此基础上提出一种稳定性主动控制方法。李长林[35]将接触约束引入结构扰动模型，分析了多滑动梁空气箔片轴承的线性稳定性，其中的动力学系数反映了箔片结构之间库伦摩擦阻尼起到的作用。

然而对于真实气体，密度-压力呈非线性关系，沿用空气轴承的处理方式(将状态方程代入雷诺方程消去热物性变量)将会使方程变得非常复杂，用来求解频变动力学特性更加困难。而且对于热力学行为严重偏离理想气体的超临界流体，需要采用基于数值计算的模型或数据库才能得到较为准确的热物性结果[36]。温建全基于小扰动法分析了S-CO2箔片轴承的动特性，并与位移速度增量法进行对比，由于在推导扰动方程时只考虑了压力扰动，从而文中对两种方法结果差异的解释为：热物性和湍流系数也会随着转子扰动发生动态变化[37]。Bi等[38]对偏导数法做出改进，通过建立动态映射引入包含扰动密度、黏度、扰动雷诺数和湍流系数的完整变量扰动，首次给出任意扰动频率下考虑真实气体和湍流动态影响的超临界二氧化碳润滑轴承动力学刚度阻尼系数的求解方法。江锦波等[39]进一步拓展了Bi等的方法，将动态离心惯性效应引入完整变量扰动，用于超临界二氧化碳干气密封的动力学特性研究。

本文给出考虑完整变量扰动的偏导数法，推导可压缩湍流润滑的扰动雷诺方程，介绍动力学系数求解方法，适用于具有附加效应的轴承与密封的频变动力学特性计算。基于准静态处理的可压缩湍流润滑雷诺方程，通过对比偏导数法和有限增量法求解的S-CO2圆柱瓦轴承刚度阻尼系数，验证了本文的偏导数法可确保动力学系数包含附加效应动态变化的影响。同时澄清基于频率扰动的和有限位移速度增量的两类方法各自的特点，并以S-CO2和油润滑高速可倾瓦轴承为案例展现本文方法的应用。

1 考虑完整变量动态变化的偏导数法

1.1 瞬态雷诺方程的频率扰动分析

为了使偏位角迭代易于处理，本文采用始于竖直位置的润滑膜圆周方向角度θ，如图1所示。一般形式的可压缩湍流润滑雷诺方程为

(1)

式中：R为轴承半径；θ为环向角度坐标；z为轴向坐标；ρ为密度；μ为动力黏度；h为润滑膜厚；p为润滑膜压力；ω为旋转圆频率；t为时间；kθ和kz分别为圆周方向和轴向湍流润滑系数，可通过局部雷诺数Re的非线性函数描述湍流效应。

kθ=12+α1Reβ1=12+0.013 6Re0.9,
kz=12+α2Reβ2=12+0.004 3Re0.96

(2)

式中： 常数α1，β1，α2，β2的值由Ng-Pan湍流润滑模型确定，该模型在来流马赫数小于5的情形对可压缩和不可压缩润滑均适用[40]；层流时上述所有常数均等于零，由于润滑膜同时存在压差流和剪切流，以Re<1 000为层流判据。

式(1)不依赖密度和黏度模型，因而不局限于特定的润滑流体，基于该方程可发展出通用的润滑特性分析方法。

(3)

式中：ψ=C/R为间隙比；C为名义半径间隙；Λ=μaω/2paψ2为轴承数，是可压缩润滑问题特有的无量纲参数；μa为环境黏度；ρa为环境密度；pa为环境压力。

(4)

在图1中，ε和φ为轴颈在扰动状态下某一瞬时所在位置的偏心率和偏位角； 对应的ε0和φ0为轴颈稳态平衡位置的偏心率和偏位角；Ob为轴承中心；Oj和Oj0为瞬时位置和稳态平衡位置的轴颈中心。

在扰动状态下，可压缩湍流雷诺方程中的密度、黏度和湍流系数同样都发生与压力扰动和膜厚扰动耦合的动态变化。本文给出雷诺方程中完整变量的扰动形式并建立与压力和膜厚扰动的动态映射关系，进而可将基于数据库和数值方法的热物性模型纳入该求解体系。

(5)

(6)

局部雷诺数Re的动态变化反映了湍流对动力学特性的影响，引入由名义半径间隙、环境密度和环境黏度定义的平均雷诺数Rea=ρaωRC/μa，可得由无量纲变量表示的动态局部雷诺数。

(7)

由幂函数的泰勒公式，结合式(6)有

(8)

将式(4)、式(6)和式(8)代入式(7)，略去二阶及以上的扰动项，得到动态局部雷诺数Re的扰动形式。

(9)

式中，Re0为对应于轴颈平衡位置的稳态雷诺数，扰动雷诺数Red是扰动密度、扰动黏度和扰动膜厚的联合响应。

利用幂函数的泰勒公式不难推出Reβ。

(10)

(11)

1.2 稳态和扰动雷诺方程

(12)

(13)

1.3 动力学系数的求取方法

(14)

式中：L为轴承长；D为轴承直径。

式(15)可将动力学系数从耦合坐标系变换到x-y直角坐标系。

(15)

式中， [A]为旋转矩阵，其表达式为

(16)

数值积分公式和旋转矩阵的具体表达式与坐标轴指向和轴颈旋转方向有关，式(14)和式(16)对应图1的坐标系。对于可倾瓦轴承，扰动膜厚表达式中还包含瓦块摆动和支点运动等相关的扰动变量，并需要借助瓦块-支点动力学方程推导折合刚度和阻尼系数，具体过程同Yang等的研究。

2 完整变量扰动偏导数法的验证

2.1 有限增量法(仅用于验证)

(17)

(18)

求解各方向位移和速度增量对应的五个类似稳态雷诺方程形式的式(17)，通过增量压力分布与稳态压力分布的差值可提取刚度阻尼系数。

(19)

其中，

(20)

动力学系数从ε-φ坐标系转换到图1的x-y坐标系公式如式(21)所示，其中[A]与式(16)相同。

(21)

这类方法的命名并不统一，表达形式也不尽相同，文献[41]以载荷增量方式实现刚度阻尼系数的提取，方法名称中文直译为数值微分法；国内大部分资料[42]称为有限差分法，是通过在不同位置和速度下积分得到油膜力的差值实现的。温建全提出一种新的实现方式，并根据方法的特点称为位移速度增量法。本质上都是通过给定有限的轴颈位移和速度实现刚度阻尼系数求解，因而为了避免与数值方法中的差分或微分的概念造成混乱，本文称之为有限增量法，制造增量的方式可借鉴温建全的研究。

对于可压缩润滑来讲，该方法得到的是静态刚度和阻尼系数，即扰动频率无限趋近于零Ω→0+的极端情况。但是能够在准静态处理(去掉可压缩时滞项∂ρ/∂t)的前提下，自动包含密度、黏度和湍流系数的动态变化效应，因此可用来验证考虑完整变量动态变化的偏导数法。

2.2 对比验证结果

表2 S-CO2 圆柱瓦轴承的结构和运行参数Tab.2 Parameters of S-CO2 cylindrical bearing

3 案例计算结果与讨论

3.1 案例1—超临界二氧化碳高速可倾瓦轴承

表3 S-CO2高速可倾瓦轴承的参数Tab.3 Parameters of S-CO2 tilting pad bearing

需要说明的是，本章两节内容中可倾瓦轴承的偏心率指的是以装配半径间隙Cb为参照的装配偏心率，记作ε′0，以区分润滑膜厚度公式中的偏心率ε0，后者以名义半径间隙(瓦块半径间隙)Cp为参照。两个偏心率的关系ε0=ε′0(1-m)为，其中m=1-Cb/Cp为可倾瓦轴承的预负荷系数。本章两个案例的结果与讨论中所提及的“偏心率”指装配偏心率。

图5显示了偏心率0.1，转速55 000 r/min的可倾瓦轴承动力学刚度和阻尼系数随着无量纲扰动频率Ω的变化。在Ω较低时，准静态假设的刚度系数结果接近完整变量扰动，而只考虑压力扰动导致刚度结果明显偏大。随着Ω增加，刚度系数明显增大，在Ω<1.5范围内只有压力扰动的结果均偏大。

随着的增加，阻尼系数先缓慢减小，大约从Ω>0.5开始明显减小，只考虑压力扰动或准静态假设都导致阻尼系数结果明显偏大。

图6为偏心率0.8，转速20 000 r/min的结果，刚度系数的变化趋势与图5高转速小偏心情形基本一致，只是增长较缓慢。阻尼系数随着Ω的增加而缓慢增大。完整变量扰动对阻尼系数结果的影响规律与图5相同。

对于可倾瓦轴承，可压缩时滞效应增大了刚度系数而减小了阻尼系数。忽略热物性和湍流的动态变化效应将导致刚度和阻尼系数明显偏大。

3.2 案例2—高速可倾瓦油润滑轴承

由于不可压缩润滑的雷诺方程中不存在压力或密度的时滞项，油润滑可倾瓦轴承的频率效应与润滑膜自身特性无关。因而也有一类方法采用Lund摄动法处理油膜扰动并与引入频率的瓦块和支点动力学模型结合，可求解不同频率的可倾瓦轴承动力学系数[43]。本文的偏导数法能够统一处理润滑膜和结构的频率扰动。首先将式(3)各项中表征可压缩性的密度变量去掉，然后用偏导数法计算高速可倾瓦油润滑轴承的折合刚度和阻尼系数，求解方程时采用雷诺边界条件处理油膜破裂。轴承的结构形式和坐标系选自手册《Journal Bearing Data Book》[44]第二部分No.49数据的五瓦瓦上加载可倾瓦动压轴承。Someya采用Lund摄动法计算可倾瓦轴承同步动力学系数，因此在本节的偏导数法计算中，无量纲扰动频率Ω=1。轴承的结构和运行参数列于表4。

表4 高速可倾瓦油润滑轴承的参数Tab.4 Parameters of tilting pad oil bearing

图8显示了可倾瓦轴承的无量纲折合刚度和阻尼系数随偏心率的变化，图8中纵坐标为无量纲动力学系数与无量纲承载力的比值，与手册《Journal Bearing Data Book》的无量纲方式一致。当偏心率大于0.6时，对于x方向的刚度阻尼，偏导数法的层流结果比手册中数据偏大，这是由于本节的计算没有采用质量守恒边界条件处理油膜破裂。

对比考虑湍流系数扰动和只有压力扰动的湍流结果，由于扰动项不影响静特性因而两种求解条件下无量纲承载力相同，可以看出忽略湍流系数的动态变化效应将使刚度系数结果偏大，偏心率越小差异越明显；而湍流系数扰动对阻尼系数几乎没有影响。

4 结 论

本文给出考虑完整变量扰动的偏导数法，从一般形式的可压缩湍流润滑雷诺方程出发，推导扰动变量的复合动态映射表达式以及扰动雷诺方程，并求解轴承动力学刚度与阻尼系数。该方法适用于具有真实气体、湍流等附加效应的轴承与密封流体膜频变动力学特性计算。对方法进行验证并给出应用案例，计算分析的结果可得如下结论。

(1) 在准静态假设下(可压缩湍流润滑雷诺方程去掉密度时滞项)用偏导数法求解S-CO2圆柱瓦轴承的刚度与阻尼系数，与有限增量法结果对比，验证了本文所提出的偏导数法求解动力学系数的准确性与有效性。

(2) 本文的偏导数法能够统一处理结构扰动和润滑膜可压缩性共同作用产生的动力学系数的频率效应。S-CO2高速可倾瓦轴承的分析结果表明，忽略压力和膜厚以外的扰动变量将对刚度与阻尼系数结果带来很大偏差。

(3) 高速可倾瓦油润滑轴承在润滑膜湍流效应明显时，若忽略湍流系数的动态变化效应将导致小偏心率下的刚度系数结果偏大。