李青松，张珂，燕震雷
(上海应用技术大学 机械工程学院，上海 201418)
以气体代替油脂作润滑的动压气体轴承依靠其收敛间隙内气体流动产生的动压效应来承载[1]。可倾瓦动压气体轴承由多个瓦块构成，各瓦块可以自由倾斜摆动进行自我调节，具有转速高，稳定性好，能耗低，摩擦小等优良性能， 在低温、高速旋转等机械中得到应用[1-5]；
但其结构复杂，制造困难[5-6]。
可倾瓦气体轴承在国外广泛研究运用而国内研究并不成熟：文献[7]提出了气体轴承的理念；
文献[1]和文献[8]首次研制出了以空气作润滑的气体轴承；
文献[9-10]在等温、考虑可压缩性的条件下计算分析了无限长气体轴承的性能；
文献[11]以可倾瓦动压气体轴承为研究对象计算分析了其静态性能；
文献[12]基于矢量叠加原理求解分析气体轴承静态特性参数；
文献[13-14]提出以装配法计算分析轴承的静态性能；
文献[15]采用线性方法对可倾瓦气体轴承的静、动态特性进行了分析研究；
文献[16-17]利用间隙耦合法分析了可倾瓦气体轴承静、动态性能。上述研究大部分是先计算单一瓦块承载力，再进行矢量合成的数值计算方法分析可倾瓦动压气体轴承的性能参数，但对气体轴承性能分析大部分不考虑滑移边界，仅利用连续介质模型进行分析。
轴承高速旋转时，转、静子两者较小的气膜间隙可导致两者与气膜润滑作用面产生稀薄气体滑移效应[18-19]。文献[20]最早利用线性玻尔兹曼方程推导出克努森数Kn变化下的润滑方程，证实滑移边界影响动压气体轴承的性能。滑移边界模型是通过引入速度滑移边界条件修正雷诺方程的方法分析滑移效应下气体轴承的各项性能，最常用的是1阶滑移模型。文献[21]考虑气体滑移效应提出2阶滑移修正边界模型，文献[22]建立了1.5 阶滑移速度边界模型。当0.001
因此，文献[27]从分子运动理论出发详细分析了微小间隙下分子与分子壁面之间的碰撞运动，将碰撞分子划分为临界角θ和小于临界角θ碰撞壁面，建立了基于任意克努森数Kn下的新滑移修正模型，该模型与传统FK模型的计算结果相近；
改进的滑移模型不仅便于应用且在高克努森数下所提供的预测结果优于其他模型；
文献[28] 对微型动压气体轴承和三瓦可倾瓦动压气体轴承研究分析，比较了连续介质模型、1阶滑移修正模型和WU新滑移修正模型下轴承的性能参数，WU新滑移模型下计算的轴承承载力数值更接近玻尔兹曼方程解。
四瓦可倾瓦轴承主要安装于第四代核电机组高温气冷堆氦风机中，其安全稳定性能对冷却堆芯系统正常运行，防止核电站发生事故极为重要。本文以四瓦可倾瓦动压气体轴承为研究对象，工作介质为氦气，考虑滑移影响，基于牛顿-拉弗森迭代法，利用有限差分离散方程求解连续介质、1阶滑移、WU新滑移模型修正后的主控雷诺方程；
分析不同模型下轴承承载力变化以及滑移边界对轴承性能的影响。
四瓦可倾瓦动压气体轴承结构示意图如图1所示，主要由1个轴颈和4个弧形瓦块组成，各瓦块在圆周方向均匀布置且可绕其支点自由摆动以适应工作位置。可倾瓦动压气体轴承的瓦块支点位置角和放置位置因瓦块数量的不同而不同；
支点位置角度和摆放位置影响轴承气膜间隙和承载能力[30-31]。图中：R为轴颈半径，θ为轴承偏位角，γ为支点位置角，φ为周向角坐标，e为轴颈中心相对支点圆圆心偏移距离，β为瓦块的支点的位置角，z为轴向坐标，h为气膜厚度。
图1 四瓦可倾瓦动压气体轴承结构示意图
2.1 修正雷诺方程
微观尺度上存在着一些不同于宏观尺度的物理现象，此时，宏观尺度的分析理论将不再适用；
微观尺度上发生的现象因其在宏观尺度上的影响较小而被忽略，但微观尺度下，这些现象成为影响力学性能的重要因素[32-33]。
可倾瓦动压气体轴承支点位置角、摆放位置因瓦块数目不同而有差异，并对轴承气膜间隙、承载力造成影响[16,29]。在轴承的分析计算中，很多参数通常以量纲一化的形式存在；
通过参数的量纲一化，可以将物理问题转化为数学问题且避免计算过程中不必要的因素，这样可有效提高计算精度。
(1)
式中：P为量纲一的压力；
p为压力；
p0为环境压力；
H为量纲一的气膜厚度；
Cb为瓦块内圆弧半径与轴颈半径之差；
φ，λ为量纲一的周向坐标。
滑移速度为气体流动速度方向与壁面切向速度之差，其大小正比于分子平均自由程。克努森数Kn是描述微观尺度物理现象时常用的重要参数[20]，恒温下：Kn=C1/(ph)，即分子自由程σ=C1/p与h的比值，其中C1取6.5×10-4。若气膜间隙较小时，Kn可能增大至气流过渡区域(0.1 表1 3种模型的修正系数 连续介质下不考虑滑移效应，此时无边界滑移的气体速度为 (2) 式中：δ为气膜边界。 笛卡尔坐标系下气体润滑轴承的N-S方程为 (3) (4) (5) 式中：u，v，w分别为x，y，z方向的气体速度； 简化(3)，(4)和(5)式可得 (6) 根据(2)，(6)式得连续介质(无滑移)、1阶滑移、WU新滑移模型的修正气体速度为 (7) (8) (9) 式中：uR为速度项系数。 从气体速度的表达式可以看出，不同模型的速度滑移方向相同。与连续介质模型相比，滑移模型引入边界速度条件使计算结果向真实结果接近[33-34]，其中WU新滑移模型系数的值随着Kn的变化而发生改变，比1阶滑移模型能更好的适应任意Kn的真实情况，也更好地与实际物理模型吻合。 等温条件下主控雷诺方程[35-36]量纲一的形式为 (10) 联立(7)—(10)式并化简整理，得到经滑移模型修正的可压缩气体的定常量纲一的雷诺方程为 Λ=6μω/(p0Ψ2)， Ψ=Cb/R， (11) 式中：Λ为量纲一的轴承数； 对轴承而言，各瓦块可以绕其支点自由摆动，摆动角度也会影响瓦块内表面与轴颈之间的间隙[16,29]。因此，气膜厚度和偏心距与瓦块自由运动有关[36-37]。单瓦块量纲一的气膜厚度为 H=1-(m+u0)cos(β-φ)+εcos(φ-θ)+ χ0sin(β-φ)， (12) χ0=χ/Ψ，m=1-Cp/Cb， 式中：ε为偏心率(轴颈中心距支点圆中心的距离)； 轴承运行时，瓦块边界与气体相通[38]，气膜压力在瓦块内表面与轴颈表面间的量纲一化边界条件为 P(φ=φ0,λ) =1,P(φ=φ1,λ) =1， P(φ,λ=l/(2R)) =1,P(φ,λ=-l/(2R)) =1， (13) 式中：l为瓦块长度； 数值算法的选择影响求解精度和计算速度[39]。进行数值模拟计算前对单一瓦块展开并进行离散网格，因为气膜厚度的方向数量级非常大，相对来说其他方向上的数量级可以忽略不记，因此即使瓦块展开成平面，也不会对最后的计算结果产生明显影响。基于以上考虑，忽略瓦块表面曲率对计算结果的影响，将各瓦块近似展开成矩形平面，并离散成40×40的网格。采用牛顿迭代法和有限差分法求解非线性雷诺方程。构造量纲一的气膜压力P的函数F(P)为 (14) 由牛顿迭代法得 F′(Pn)δ=-F(Pn)； δ=Pn+1-Pn， (15) 式中：δ为2次迭代值之差； 使用牛顿-莱布尼茨迭代法将P的非线性方程转化为压力增量ΔP的线性方程 (16) 基于有限差分法，利用中心差分格式，将(16)式展开(i=1～41，j=1～41)为 (17) 利用(17)式将(16)式整理为有规律的系数矩阵方程组，即 Ai,jPi,j+Bi,jPi+1,j+Ci,jPi-1,j+Di,jPi,j+1+ Ei,jPi,j-1=-Hi,j。 (18) 有限差分法离散得到各系数项为 根据已知条件迭代求解矩阵方程组以获得单个瓦块的量纲一的压力分布。迭代过程中，为保证计算过程的稳定收敛和加快收敛速度，一般采用松弛迭代法计算 (19) 为了确保迭代过程可终止，必须保证每次迭代计算的压力分布满足计算精度。按照收敛准则，当计算结果达到精度要求后，程序将跳出循环，从而终止迭代计算，允许相对误差eeps=1×10-4，收敛准则为 (20) 3.3.1 计算瓦块摆角 瓦块可以自由摆动，计算时需要判断摆角的位置。当瓦块处于静平衡位置时，瓦块不再摆动，此时的角度就是计算轴承静态特性时的摆角位置，瓦块的摆动角度与瓦块的支点力矩M、气膜压力分布和量纲一的力矩满足平衡关系，即 (21) 当计算结果不满足瓦块支点力矩平衡时，使用超大范围牛顿公式迭代产生新的瓦块摆角，迭代更新瓦块摆角直至支点力矩M为0，迭代公式和收敛准则为 (22) (23) 3.3.2 计算气膜合力 所有瓦块的合力满足静态平衡条件时，可得各单个瓦块的气膜压力分布，通过矢量叠加的方法计算轴承的气膜合力。轴承垂直和水平方向的量纲一的气膜合力Fx，Fy及量纲一的承载力FCL计算公式为 (24) 3.3.3 计算偏位角 压力分布计算时需确定偏位角θ，偏位角迭代格式和收敛准则为 θn+1=θn-ζarctan(Fx/Fy) ， (25) Fx/Fy≤η。 为验证本文数值计算模型的准确性，以文献[10]和文献[40]的计算结果与本文3种模型的计算结果进行对比，数值计算流程如图2所示，后文若无特殊说明，所得结果皆为量纲一的值。 图2 数值计算流程 文献[10,40]与本文3种模型计算结果的对比如图3所示：将文献[10]中的可倾瓦气体轴承参数代入数值计算流程，利用连续介质模型计算不同偏心率下的轴承承载力，并与文献[10]中的计算结果相比较，可以看出不同偏心率下的误差值均小于5.8%； 图3 本文3种模型与文献的计算结果对比 给定预载系数m=0.2，偏心率ε=0.5下的轴承参数见表2。 表2 四瓦可倾瓦动压气体轴承参数[41-42] 4.2.1 气膜厚度计算结果 m=0.2，ε=0.5时，瓦块的气膜厚度如图4所示：气膜厚度沿轴颈旋转方向逐渐减小。瓦块气膜厚度随周向节点的变化如图5所示：任意位置瓦块2，3周向节点方向气膜厚度小于瓦块1，4； (a)瓦块1 图5 瓦块量纲一的气膜厚度随周向节点的变化 4.2.2 气膜压力和克努森数的计算结果 以瓦块1为研究对象，考虑滑移效应，取m=0.2，ε=0.5，3种模型下瓦块1气膜压力分布(左图)和等高线(右图)如图6所示，连续介质模型下瓦块1的克努森数Kn分布如图7所示。滑移效应使瓦块1的气膜压力明显降低，且WU新滑移模型计算所得气膜压力分布小于1阶滑移模型； 4.2.3 偏心率对量纲一的气膜厚度、克努森数和承载力的影响 以瓦块1为研究对象，在m=0.2，ε取值为0，0.2，0.4，0.5，0.6，0.8时，探讨偏心率的变化对气膜厚度、克努森数、不同修正模型下轴承承载力的影响。瓦块1气膜厚度随偏心率的变化如图8所示：随偏心率增大，瓦块气膜厚度呈增大趋势。瓦块1的Kn值沿轴向节点随偏心率的变化如图9所示：偏心率越大，沿周向节点Kn值越大。不同模型下轴承承载力随偏心率的变化如图10所示：随偏心率增大，承载力FCL增大； 图8 瓦块1量纲一的气膜厚度随偏心率的变化规律 图9 瓦块1克努森数随偏心率的变化 图10 轴承承载力随偏心率的变化 因瓦块关于支点对称布置，所以瓦块1，2承载力相同，瓦块3，4承载力相同，计算不同模型下单个瓦块承载力随偏心率的变化规律如图11所示：对于单个瓦块，随偏心率增大，瓦块1，2承载力减小，瓦块3，4承载力增大。考虑滑移修正下计算所得承载力显著下降，且WU新滑移模型计算所得轴承承载力小于1阶滑移模型。 (a)瓦块1,2 4.2.4 预载系数和量纲一的轴承数对轴承承载 力的影响 以瓦块1为研究对象，探讨预载系数m和量纲一的轴承数对不同模型下轴承承载力的影响。ε=0.5时，不同模型轴承承载力随预载系数m的变化如图12所示：对同一个模型，预载系数m增大，轴承承载力随之增大； 图12 轴承承载力随预载系数的变化 m=0.2，ε=0.5时，3种模型下承载力随量纲一的轴承数的变化如图13所示：随轴承数增加，3个模型的轴承承载力均增加； 图13 轴承承载力随量纲一的轴承数的变化 m=0.2，ε=0.5时，3种模型单个瓦块承载力随量纲一的轴承数的变化如图14所示：随轴承数增大，瓦块1，2的承载力呈先增大后减小； (a)瓦块1，2 以四瓦可倾瓦动压气体轴承为研究对象，考虑滑移效应的影响并分析滑移模型对量纲一的气膜压力分布和承载力的影响。从气体润滑基本方程出发，利用1阶滑移模型和WU新滑移模型的气体速度推导出不同模型的定常量纲一化静态雷诺方程，讨论连续模型、1阶滑移模型、WU新滑移模型在不同轴承气膜厚度、压力等条件下的承载力。结论如下： 1)单个瓦块量纲一的气膜厚度均沿轴颈旋转方向逐渐减小； 2)可倾瓦动压气体轴承的承载力与偏心率、预载系数和量纲一的轴承数都成正比； 3)WU新滑移模型计算的承载力、气膜压力结果最小，1阶滑移模型次之，连续模型最大；
μ为气体动力黏度系数；
ρ为气体密度；
t为时间。
ω为轴颈旋转角速度。2.2 气膜厚度方程
χ0为量纲一的瓦块摆角；
m为预载系数；
Cp为支点间隙；
u0为量纲一的静态径向位移，取u0=0。2.3 边界条件
φ0为瓦块前缘进气端的角度；
φ1为瓦块后缘出气端的角度。3.1 牛顿迭代法转化雷诺方程
n=0,1,2…，
n为迭代次数。3.2 基于有限差分法离散方程
3.3 求解流程
4.1 程序验证
将文献[40]中的可倾瓦气体轴承参数代入数值计算流程，分别利用连续介质模型、1阶滑移模型和WU新滑移模型计算不同偏心率下的轴承承载力，并与文献[40]中的计算结果相比较，可以看出连续介质模型与文献[40]结果之间的误差最小，3种模型的误差值均小于7.2%，考虑滑移边界条件下所计算出的可倾瓦气体轴承的承载力相比文献[40]和连续模型显著降低。由上述结果可知，图2程序计算流程具有一定的准确性。
4.2 计算结果
瓦块1气膜厚度最大，瓦块3气膜厚度最小。
在瓦块1的边界处，气膜压力等于大气压力，3种模型在边界处计算的气膜压力都相同，接近出口位置时的气膜压力达到最大值；
除了边界处，连续介质模型计算出的气膜压力大于滑移模型，WU新滑移模型计算的气膜压力值小于1阶滑移模型。另外，气膜压力沿轴向方向先增大后减小，分析原因：瓦块1摆角在气体冲击下沿轴承旋转方向摆动，瓦块内表面与轴颈表面间隙逐渐减小，气膜压力逐渐增大，并在瓦块后段接近后缘部分的最小间隙处达到最大；
瓦块后缘出气端与气体环境相同，后缘部分压力急剧减小至和外界压力相同；
由于可倾瓦块沿轴向两端都与大气环境相同，所以气膜压力沿轴向表现为先增大后减小。图7中一部分润滑气膜处于气流过渡区，另一部分润滑气膜则处于滑移区；
WU新滑移模型在处理气流过渡区的计算性能优于连续介质模型和1阶滑移模型。
WU新滑移修正模型计算所得承载力结果最小，1阶滑移模型次之；
随偏心率增大，3种模型计算结果数值偏差也增大。
对相同预载系数，WU新滑移模型计算所得的承载力最小。
相同轴承数时，考虑滑移效应的模型轴承承载力较小；
随着轴承数增加，3种模型的轴承承载力差值减小。这是因为量纲一的轴承数只与角速度ω有关，因此随着角速度速的增加，轴承数增加，瓦块内表面与轴颈表面间隙产生的动压效果增强，滑移效应对轴承的影响增大，气膜压力增大，轴承承载力随之增大；
量纲一的轴承数一定时，考虑滑移效应的影响，瓦块内表面与轴颈表面间的气膜压力减小，导致承载力减小。
瓦块3，4的承载力随轴承数的增大而增大。考虑滑移效应的修正模型所计算承载力结果明显较小，即WU新滑移模型在承载力计算上有优势。
任意位置瓦块2，3周向节点量纲一的气膜厚度小于瓦块1，4。
滑移修正使轴承承载力显著下降，且WU新滑移模型计算的承载力小于1阶滑移模型；
随着偏心率、预载系数的增大，滑移效应对承载力结果的影响也增大。
基于滑移效应客观存在的事实，使用WU新滑移模型对微小间隙下可倾瓦动压气体轴承的静态性能计算更为准确。