动压轴承流固界面非一致滑移流动分析方法及分布特征

张镜洋，陈哲，吕元伟，孙义建，张靖周，陈卫东，3，罗欣洋

1.南京航空航天大学 航天学院，南京 211106

2.南京航空航天大学 能源与动力学院，南京 210016

3.南京航空航天大学 机械结构力学及控制国家重点实验室，南京 210016

4.南京机电液压工程研究中心 航空机电系统综合航空科技重点实验室，南京 211106

非接触式动压气体轴承具有轻质、低摩阻、无油污染等优势，作为转子支撑机构可大幅提高高速旋转机械的极限转速和功重比，在制冷透平、冲压涡轮、冷却风扇、小型涡轮发动机等航空航天旋转机械上具有突出的应用优势［1-4］。以大预紧力装配和高温涂层为特征的第四代多叶波箔动压气体轴承显著改善了承载性能和极端工况下的运行稳定性［5］。该类轴承利用转子与流体间粘性剪切力驱动气体周向流动，在变截面微间隙内形成动压压缩效应而提供承载。动压气体轴承转静间隙气膜流动为典型的“泰勒-库特”流动［6］。由于气膜流动间隙尺度小（10−6m 量级）、转子转速高（105r/min 量级）、剪切效应强，气膜流动达到滑移流范畴［7］。在稀薄效应和强剪切复合作用下，微尺度间隙内会发生近壁面气体界面滑移现象。由于界面滑移与环楔形间隙承载气膜压力分布强相关，轴承承载等核心性能对其极其敏感，动压气体轴承间隙内气固界面滑移机制及其对轴承性能的影响成为该领域当前的研究热点。

众多研究者相继利用滑移速度模型和极限剪应力模型等对动压气体轴承间隙气膜滑移流动开展研究。学者认为滑移速度模型下界面滑移速度与速度梯度成正比，逐步发展出1.0 阶［8］、2.0 阶［9］、1.5 阶滑移速度模型［10］，并基于上述模型推演出修正的雷诺方程，系统地研究了滑移对动压气体轴承间隙内压力分布和气膜承载力的影响。譬如，冯凯［11］和黄海［12］等分别基于1.0 阶和2.0 阶滑移速度模型，分析了轴承间隙气膜压力分布以及承载力变化规律。在此基础上，Wu［13］推演了基于分子动力学理论的滑移速度模型，解决了克努森数Kn>1 时气膜间隙流动处于分子流状态下的界面滑移模型问题，后续兰正义和伍林［14］、燕震雷和伍林［15］利用Wu 滑移模型分析了界面滑移对不同类型动压气体轴承性能的影响。极限剪应力模型则认为当流体界面剪应力达到流固界面切应力的极限值时才发生滑移［16］。Spikes和Granick［17］提出了一种极限剪应力为滑移判据复合滑移速度模型的轴承微间隙流动界面滑移模型。马国军［18］则通过计算和实验发现滑移速度模型适合低剪切界面滑移，极限剪应力模型更适合高剪切率的界面滑移。王丽丽等［19］以径向动压轴承为研究对象，基于极限剪应力推导了无滑移、转子侧滑移、静子侧滑移、两侧均滑移4 种不同滑移状态雷诺方程，其分析结果表明了轴承间隙流动处于4 种不同滑移状态时承载压力存在明显差异。李旺［5］、张镜洋等［20］以2 种类型的波箔型动压气体轴承为研究对象，基于极限剪应力模型分析了不同界面滑移状态对轴承静特性的影响，结果表明箔片侧滑移对静特性有利，轴面侧滑移对静特性不利。近年来，Bhattacharya［21］、Rao［22］等设计了动压气体轴承局部滑移速度梯度分布的气固界面，发现不同局部滑移分布形态能够显著影响轴承性能。实际上，为了形成提供承载的气膜压力，实现在负载作用下动压轴承转子的偏心转动，轴承楔形间隙周向流动尺度原本就存在差异，有研究表明其周向克努森数Kn相差可达一个数量级［23］，这会导致近壁面气体界面滑移状态存在非一致性。而前述研究中，所涉及的滑移速度模型多认为近壁面气体周向各局部同步产生滑移；所涉及的极限剪应力模型虽然考虑了转静两侧不一致滑移状态，但是在周向上多依据近壁平均剪应力关系来计算一致的滑移状态。滑移状态一致假设或人为设定滑移速度分布，与轴承间隙周向非一致的流动尺度及可压缩气体状态参数相矛盾。尤其在飞行器大过载间隙大偏心下这种矛盾会被加剧，而且第四代多叶波箔型动压轴承具有多楔型间隙特征、间隙内流动尺度分布更为复杂［5］。因此，工程应用中亟待发展一种适用的动压气体轴承间隙周向非一致滑移流动模型，并能掌握其复杂滑移状态的演化机制。

本文以多叶波箔型动压气体径向轴承为研究对象，推演局部界面剪应力与极限剪应力计算局部滑移速度的方法，建立基于剪应力方程、滑移修正雷诺方程、气膜厚度方程耦合迭代的流固界面非一致滑移流动分析方法；通过与传统滑移模型流场分布特征以及实验值的分析比较验证模型有效性；进而数值研究旋转及几何参数变化下的轴承间隙滑移流动演化规律。以期为飞行器高速旋转机械非接触转子支撑结构设计提供技术支撑。

1 物理模型和数值方法

1.1 物理模型

多叶波箔动压气体轴承由轴承套、平箔和波箔组成，波箔作为平箔的支撑元件，一端固定在轴承套上，如图1（a）所示。其中，O为轴承套圆心，O′ 为轴颈圆心，E0为偏心距，偏心率e=E0/c，c为轴承间隙尺度，ω为旋转角速度，θ为转动角，P0为环境压力。在转子偏心和装配后的多楔形箔片结构下流道宽度即气膜厚度h沿周向呈现明显变化，并且周向压力P多峰非均匀分布，这是引起周向克努森数强不均的原因，如图1（b）所示。轴承结构参数如表1 所示。

表1 轴承模型结构参数Table 1 Structural parameters of bearing model

图1 轴承结构和间隙参数示意图Fig.1 Schematic diagram of bearing structure and pa‐rameters in clearance

1.2 滑移速度模型与极限剪应力模型

1）滑移速度模型

滑移速度模型由牛顿内摩擦定律和分子动力学方程导出，适用于描述微尺度下滑移流动边界条件，基本表达为

式中：us为滑移速度；u为周向线速度；y为径向位置；λ为平均分子自由程；α为动量协调系数；b为速度偏导二阶项系数。

式（1）中关联了近壁法向速度梯度，反映了剪应力对滑移的影响；平均分子自由程λ、动量调节系数α组成速度梯度项系数，反映微尺度下分子与壁面作用对于界面滑移的影响程度；利用速度偏导二阶项系数b修正了速度梯度高阶影响程度。为适应不同特征下的滑移流动模拟，学者通常将滑移速度模型分为1.0 阶、2.0 阶、1.5 阶，三者的区别主要在于二阶项系数b，分别为0、1/2、2/9。动压轴承领域学者也分别采用过上述3 种滑移速度模型［8-10］对其间隙流动界面滑移进行模拟，利用滑移速度修正原无滑移边界速度。轴面侧、箔片侧近壁气体线速度u1′、u2′的计算公式分别为

式中：u0为轴颈旋转线速度。

2）极限剪应力模型

式中：m、n分别为周向和轴向网格数；x为周向位置；Pa为近壁面气体的平均压力；τ0为初始剪应力；β为调节系数。

3）传统滑移模型的问题分析

动压气体轴承间隙尺度为μm 级，其间隙流动克努森数Kn范围为10−5～10−2［15，23］，间隙内局部处于滑移流范畴。从前述动压气体轴承界面滑移模型应用研究中可以看出，以往多认为近壁面气体周向滑移状态一致，即认为当流固界面出现滑移时则整面各局部均滑移，未细致考虑流固界面会出现局部滑移与局部非滑移区的混合状态。且滑移速度模型运用流固界面各处相同的动量协调系数α经验值，极限剪应力模型则通过比较面平均的界面剪应力、极限剪应力，判断其有无界面滑移或计算滑移速度大小。然而，动压轴承转子的偏心转动和间隙高度的周向非一致，会诱发周向强不匀且非线性的径向速度梯度分布（见图2），并且周向不均匀压力会使得局部可压缩气体密度等物性变化较大，表征局部流动尺度的克努森数Kn也会有较大差异。此情形下动压气体轴承间隙周向应具有不同的滑移状态（见图2）或尺度。但是，在传统滑移模型中这种影响滑移的局部非一致性因素并未被充分考虑。图2 中：u1、u2分别为无滑移时轴面侧和箔片侧线速度；us1、us2分别为轴面侧和箔片侧滑移速度；下标A、B、C 分别代表图1（a）中不同截面。

图2 不同截面界面滑移示意图Fig.2 Schematic diagram of interface slip in different sections

1.3 界面非一致滑移流动模型

本文考虑近壁面气体周向局部滑移状态的非一致性，建立滑移修正雷诺方程，在此基础上推演局部剪应力计算局部滑移状态和滑移速度的方法。

无滑移时的无量纲雷诺方程［11-12］为

考虑了轴面与箔片两侧不同的滑移速度，引入无量纲近壁气体速度对式（9）进行修正，则得到滑移修正的无量纲雷诺方程［7-15］为

式中：±取舍原则为保证τ为正值。

根据牛顿内摩擦定律，轴面侧、箔片侧界面剪应力τ1、τ2分别为

借鉴极限剪应力经验关系［16］，获得关联局部气膜压力的极限剪应力方程

式中：τs为局部极限剪应力；τ0为局部初始剪应力；β为局部调节系数。因两侧材料及压力相同，故两侧具有相同的τs，根据文献［18］，取τ0=6 Pa，β=0.000 7。

根据前述界面剪应力和极限剪应力公式，结合式（17），可求出后文轴面侧和箔片侧局部滑移速度。

根据局部剪应力和极限剪应力方程［24］，可求出轴面侧和箔片侧局部滑移速度us1、us2，计算公式为

式（18）、式（19）说明若局部近壁气体流固界面剪应力大于极限剪应力，则判断局部发生滑移，反之则局部不发生滑移，即滑移速度为0。如此，则利用局部界面剪应力和极限剪应力关联了局部流动状态，来计算周向非一致的滑移速度。

1.4 耦合箔片变形的气膜厚度方程

无量纲气膜厚度的表达式为

式中：hf为箔片变形量。

箔片变形矩阵根据柔度和压力矩阵［5］求出

式中：U为柔度矩阵；P为气膜压力矩阵。

1.5 网格划分和边界条件

对轴承间隙内气体进行周向和轴向网格划分，将圆柱面沿周向展开，如图3 所示。将周向和轴向分别划分为m、n个单元，形成m+1、n+1个节点；将轴承进口和出口设置为大气压；周向展开时轴向两侧为循环边界，具体边界条件为

1.6 计算流程和工况

根据前述模型建立滑移修正雷诺方程、局部剪应力方程、弹流润滑气膜厚度方程耦合迭代的界面非一致滑移流动分析方法，计算流程如图4所示。将雷诺方程耦合气膜厚度方程，对其差分方程采用超松弛迭代法求解，直至气膜压力收敛，压力残差收敛公式为

图4 计算流程图Fig.4 Calculation flow chart

式中：i、j分别为周向和轴向网格点坐标，取值范围为［0，m］和［0，n］；k为雷诺方程迭代周期数；k′为滑移速度迭代周期数。然后通过两侧剪应力计算两侧局部滑移速度，并结合计算出的两侧局部滑移速度计算滑移速度残差，轴面侧和箔片侧滑移速度残差可分别表示为

判断轴面侧和箔片侧局部各节点滑移速度残差是否收敛，若不收敛则重新计算气膜厚度和气膜压力，经多次循环迭代后可获得计算域内所有节点的局部滑移速度、滑移状态、气膜压力、气膜厚度以及承载力，且承载力计算公式为

本文计算工况参数范围为：轴承数Λ∈[0.15，0.90]、偏心率e∈[0.1，0.8]。

1.7 计算方法验证

对本文提出的数值计算方法进行验证。首先，保持轴承参数与文献［25］中实验一致，并将计算结果与文献［25］实验结果进行对比，如图5（a）所示。结果表明：不同滑移模型计算的承载力随转速变化规律一致，转速较低时，不同滑移模型计算结果接近，随转速升高，不同滑移模型计算结果均高于实验值。与其他模型相比，非一致滑移模型的计算结果与实验值最为接近，最大误差为8.5%。此外，构建轴承间隙压力分布实验系统，实验装置如图5（b）所示。实验测试了轴承间隙c=600 μm、偏心率e=0.6 时轴承Z=L/2 处周向压力分布，将实验结果与计算和结果进行对比，如图5（c）所示。结果表明：n0=6×103，9×103r/min 工况下，数值计算结果与实验测得的气膜压力分布规律一致，且除高压区和低压区附近外结果接近，实验与数值模拟平均表面压力的误差在11%以内。与n0=6×103r/min工况相比，转速提升至n0=9×103r/min 工况，高转速下轴自身的动力学变形效应，使得转轴变得不稳定，气膜厚度处于变化之中，尤其在轴承气膜压力高压区与低压区位置附近处（位于最小和最大气膜厚度附近）对气膜厚度变化极度敏感，这使得理论值与实验值相差较大。这说明了本文数值计算方法的有效性。

图5 有效性验证Fig.5 Validity verification

2 计算结果与分析

2.1 不同滑移模型对压力和滑移速度分布影响

图6 为偏心率e=0.6 时不同工况下局部克努森数Kn分布，从图中可以看出，轴承间隙内局部克努森数Kn沿周向分布变化明显，最大可达1.8×10−3、最小为4×10−4，相差几近一个数量级。尤其在周向θ=0.8π～1.5π 位置最小气膜厚度附近，局部克努森数Kn的范围为10−3～10−1，处于滑移流区。

图6 不同工况下局部克努森数Kn 分布（e=0.6）Fig.6 Kn distribution under different working condi‐tions（e=0.6）

图7、图8 分别为Λ=0.45，0.90 工况下滑移速度分布。可以看出，在两工况下，一阶滑移速度模型下轴面侧和箔片侧流固界面均发生滑移，局部滑移速度数值差异较小。极限剪应力模型下Λ=0.45 工况两侧流固界面均不发生滑移；而Λ=0.90 工况轴面侧流固界面均发生滑移，箔片侧流固界面均无滑移。而两工况中，非一致滑移模型下两侧局部发生滑移，滑移速度呈现多峰值分布，与气膜厚度分布状态相似，局部克努森数Kn较大区域滑移速度较大、较小区域无滑移或者滑移速度较小。本文数值方法获得的滑移速度分布与局部克努森数Kn分布更符合。

图7 Λ=0.45 时滑移速度分布Fig.7 Slip velocity distribution at Λ=0.45

图8 Λ=0.90 时滑移速度分布Fig.8 Slip velocity distribution at Λ=0.90

图9 为不同滑移模型下轴承Z=L/2 处压力分布。与无滑移模型相比，一阶滑移速度模型下Λ=0.45 工况的高压区压力峰值降低、低压区压力无明显变化，最大影响小于1.5%。Λ=0.90 工况的高压区压力峰值降低、低压区压力峰值升高，但影响较小，最大影响小于5%；极限剪应力模型下Λ=0.45 工况的压力分布与无滑移模型完全相同。Λ=0.90 工况的高压区和低压区压力峰值均降低，对高压区影响较大，可达20%，对低压区影响较小，最大为6%；非一致滑移模型下Λ=0.45 工况的高压区压力峰值升高、低压区压力峰值降低，最大影响小于2.3%。Λ=0.90 工况的高压区压力峰值降低、低压区压力峰值升高，最大影响分别为10.6%、5.0%。不同滑移模型与无滑移相比，压力分布趋势相同，但对决定轴承承载能力的压力峰值影响较大，尤其是轴承数越大影响程度越大，相较于无滑移模型分析结果，非一致滑移模型下压力峰值变化幅度位于滑移速度模型和极限剪应力模型之间。产生上述现象的具体原因：滑移速度模型和极限剪应力模型认为轴承间隙各处流动尺度相当，采用间隙各处相同的滑移模型参数，导致过高或过低估计边界滑移速度的影响；局部流动状态与动量协调系数α及面平均的剪应力-τ、-τs强相关，而传统滑移模型中并未考虑滑移速度与局部流动状态的耦合影响，在迭代收敛过程中仅进行了单次赋值，因此均出现了滑移速度分布不合理的情况。而本文提出的非一致滑移模型，考虑了偏心间隙尺度及局部压力分布对剪应力的影响，并将剪应力计算与雷诺方程进行耦合迭代，关联了局部流动状态和滑移速度的相互关系。从图10 中可以看出，经多次迭代后流固界面剪应力变化幅度较大，这证明了剪应力计算与雷诺方程耦合迭代的必要性。经耦合迭代后局部界面剪应力收敛至小于等于局部极限剪应力的合理范围。因此，关联滑移与局部流动状态进行动压气体轴承间隙滑移流动分析，可以有效解决周向非一致滑移速度的准确预计问题。

图9 不同滑移模型下气膜压力分布Fig.9 Pressure distribution at different slip models

图10 压力和剪应力随迭代次数变化Fig.10 Variation of pressure and shear stress distribu‐tion with iteration times

2.2 不同轴承参数对界面非一致滑移特征影响

2.2.1 轴承数

轴承数反映了转静子相对运动速度和间隙尺度的平方之比的大小，对轴承动压压缩效应及周向滑移速度有着显著的影响。图11 为偏心率e=0.6 时非一致滑移模型下滑移速度随轴承数变化规律。随轴承数增大，箔片侧与轴面侧依次出现滑移。Λ=0.30 工况下，仅箔片侧发生滑移；Λ=0.45 工况下，轴面侧也发生滑移。随轴承数增大，轴面侧和箔片侧滑移区域沿周向从轴承最小气膜厚度处逐渐向两侧扩张且滑移速度不断增大，箔片侧局部滑移区域小于轴面侧，轴面侧滑移区域面积最高约为箔片侧的2 倍。图12为不同工况下Z=L/2 处剪应力分布，从图中可以看出，两侧剪应力峰值随轴承数增大而增大，易发生滑移。如图12（a）所示，Λ=0.30 工况下，仅箔片侧部分区域的气体界面剪应力达到极限剪应力；如图12（b）所示，Λ=0.60 工况下，轴面侧部分区域的气体界面剪应力也达到极限剪应力，故箔片侧和轴面侧依次发生滑移。从图12（c）中可以看出，随轴承数的进一步增大，两侧气体界面剪应力达到极限剪应力的区域沿周向从最小气膜厚度处向两侧逐渐扩张，且剪应力大处滑移速度较大。其中轴面侧剪应力更大，故轴面侧滑移区域扩张范围更广，这是两侧滑移速度分布随轴承数变化规律的产生内因。

图11 滑移速度分布随轴承数变化规律Fig.11 Variation of slip velocity distribution with bearing number

图12 不同工况下Z=L/2 处剪应力分布Fig.12 Shear stress distribution under different working conditions at Z=L/2

2.2.2 偏心率

偏心率反映了轴承间隙周向尺度变化幅度大小，同样对周向滑移速度分布有显著影响。图13 为轴承数Λ=0.60 时非一致滑移模型下滑移速度随偏心率变化规律。与轴承数增大时相同，随偏心率增大，箔片侧与轴面侧依次出现滑移。e=0.2 工况下，仅箔片侧发生滑移；e=0.5 工况下，轴面侧也发生滑移。随偏心率增大，轴面侧和箔片侧滑移区域沿周向从最小气膜厚度处逐渐向两侧扩张且滑移速度均不断增大，箔片侧滑移区域小于轴面侧。图14为不同工况下Z=L/2 处剪应力分布，从图中可以看出，两侧剪应力峰值随偏心率增大而增大，易发生滑移。如图14（a）所示，e=0.2 工况下，仅箔片侧部分区域的气体界面剪应力达到极限剪应力；如图14（b）所示，e=0.5 工况下，轴面侧部分区域的气体界面剪应力也达到极限剪应力。与轴承数增大时不同的是，偏心率增大时两侧滑移区域扩张范围较小，沿周向远离最小气膜厚度处的界面剪应力逐渐降低。其内因是增大偏心率使得间隙周向通道高度变化更为剧烈，偏心率增大仅使得最小气膜厚度附近径向速度梯度增大，而沿周向从最小气膜厚度处向两侧界面剪应力迅速降低，限制了滑移区域的扩张，这不同于轴承数提高整体界面剪应力的情况。

图13 滑移速度分布随偏心率变化规律Fig.13 Variation of slip velocity distribution with eccentricity

图14 不同工况下Z=L/2 处剪应力分布Fig.14 Shear stress distribution under different working conditions at Z=L/2

3 结论

为获得动压气体轴承变截面微间隙局部滑移特征，开展了界面非一致滑移分析方法及滑移速度演化规律研究，得到如下主要结论：

1）针对传统滑移模型未考虑变截面间隙流固界面局部流动尺度差异问题，推导了局部剪应力与滑移速度的关联方程，建立了局部剪应力方程、滑移修正雷诺方程、弹流润滑气膜厚度方程耦合迭代的流固界面非一致滑移流动分析方法，有效解决了周向非一致滑移速度的预计问题。

2）在本文研究参数范围内，与传统滑移模型相比，非一致滑移流动分析方法获得的滑移速度分布与局部克努森数Kn分布一致性更好。与无滑移相比，不同滑移模型压力分布趋势相同，但对决定轴承承载能力的压力峰值影响较大。一阶滑移速度模型和极限剪应力模型对压力幅值最大影响分别为5%、20%，而非一致滑移模型介于两者之间，最大影响为10.6%。

3）非一致滑移模型下，随轴承数和偏心率的增大，箔片侧与轴面侧依次出现滑移，且两侧滑移区域沿周向从轴承最小气膜厚度处逐渐向两侧扩张，轴面侧滑移区域面积最大为箔片侧的2倍。但轴承数增大时滑移区域面积迅速增大，而偏心率增大时引起局部流动尺度差异变大，仅最小气膜厚度处的滑移面积增大。