滑移边界下气膜厚度阶跃变化的阶形板间气体流动

张建勋， 伍 林

(上海应用技术大学机械工程学院，上海 201418)

微尺度下的气体流动中由于存在气流的非平衡效应，气体介质的连续型分布假设不再适用，通过Boltzmann方程求解微尺度下的气体流动[1-2]存在数学上的复杂性，应用边界速度滑移修正 Navier-Stokes(N-S)方程进行求解相对快速简捷且模型适应性较好。

基于Maxwell基本假设提出的滑移系数为固定值的一阶滑移[3]边界条件，适用于气膜间克努森数处于滑流区间时的微尺度流动。杨琴等[4-5]通过矩方法对Boltzmann方程求解并推导出了滑移系数随克努森数变化的扩展速度滑移边界条件[6]，在近过渡区内与Boltzmann方程解吻合性较好，但是在克努森数接近甚至大于1时偏差仍然较大。Wu[7]从分子动力学理论出发，细化了微尺度下气体分子与固体壁面间的碰撞运动并将其分为两类进行研究，推导出了滑移系数随气膜间克努森数变化且适用于任意克努森数的新滑移边界条件，在气膜厚度连续变化的楔形模型行中与F-K模型解呈现良好的吻合性。

基于阶形板的动压原理[8]衍生出了一系列沟槽类气体轴承，如广泛应用于干气密封等工业领域的平面螺旋槽气体止推轴承。Zhang等[9]、Wang等[10]应用不同的数值方法对平面螺旋槽模型进行了求解。但沟槽类气体轴承气膜厚度存在阶跃变化，槽区和台区克努森数跨区域分布，当台区克努森数进入过渡区时槽区可能还处于滑流区，基于无滑移和一阶滑移边界条件得到的数值结果将会与现实结果产生较大偏差。作为沟槽类气体轴承的基本构成单元，微尺度下阶形板间的气体流动特性的深入研究对沟槽类气体轴承的发展具有重要意义。

代入新滑移边界条件得到修正后的雷诺方程，应用有限体积法结合Newton-Raphson方法迭代求解了不同参数下阶形板模型中的气膜压力和克努森数分布曲线，沿气膜厚度方向的流速分布曲线，气膜承载力随速度参量的变化趋势等数值结果，对比了不同边界条件下所得计算结果之间的偏差并分析其产生原因，为相关微型沟槽类气体轴承的设计提供参考。

1 模型介绍

在流体润滑中为提高动压轴承的相关性能，转子或轴承上可以开有各种形式的槽，以提高轴承在工作过程中产生的动压效应。阶形板模型所产生的气膜压力分布曲线可以看作各种沟槽类动压气体轴承气膜压力分布的一个基本单元。图1为阶形板间的气体流动模型示意图，上板为台阶形平板，下板为平行板。槽区高度为h1，宽度为b1，台区高度为h2，宽度为b2，假定上板固定不动，下板以速度u0沿X轴正方向移动。

图1 阶形板间气体流动示意图Fig.1 Schematic diagram of gas flow between stepped plates

关于阶形板间的流体流动产生的动压效应可以用流量守恒原则解释是其物理原理，进气端与出气端的流量应该相等，槽区和台区的流量也应该相等。阶形板间气体流速在槽区的分布为斜边内凹的直角三角形，台区为斜边外凸的直角三角形。相应的槽区应存在沿X轴正方向的正压力梯度，台区应存在沿X轴正方向的负压力梯度，在槽台交界处取得压力极值。图1中气体流速分布示意图不考虑滑移边界条件，仅作为阶形板间动压原理分析阐述说明。

2 模型求解

2.1 滑移修正控制方程

微尺度下气膜间气体分子的非连续型运动程度增强，阶形板间存在上升和下降的气膜压力分布，应用滑移系数固定的一阶滑移边界条件统一描述各个位置的气体分子界面滑移状态显然是不符合现实物理规律的，新滑移边界条件不仅滑移系数随气膜间当地克努森数的变化而变化，而且结合分子自由程与气膜厚度的关系将气体分子的碰撞运动分为两类进行讨论，可以更好地描述微尺度下的气体流动特性。

式(1)和式(2)为根据简化的N-S方程推导出的流速u和流量方程Qx：

(1)

(2)

式中:u为气体流速；Qx为气体流量；p为气膜压力；ρ为气体密度；h为气膜厚度；λ为分子自由程;η为气体黏度系数；M、N为滑移系数，如表1所示。表1中α为分子切向动量协调系数，即在边界层处发生漫反射的气体分子所占的比例，克努森数Kn为分子自由程与气膜厚度的比值，f=min[1/Kn,1]。

表1 滑移系数Table 1 Slip coefficients

(3)

根据式(3)定义可得到无量纲化的雷诺方程如式(4)所示。

(4)

式(4)中：P为气膜压力；h为气膜厚度；x为模型特征长度；Λx为定义轴承数。

2.2 限体积法网格划分与方程的离散

因阶形板模型在槽台分界处存在气膜厚度阶跃变化的特性，若要对边界压力极大值进行求解，需应用有限体积法将计算参量转移到有限体积边界，解决气膜厚度存在的奇异变化，具体的求解区域网格点划分形式如图2所示。

W为相邻左侧单元网格中心线；E为右侧单元网格中心线；C为有限体积单元中心线；w、e分别为有限体积单元的左右边界图2 网格单元划分Fig.2 Grid cell division

根据阶形板的几何特性共分为槽区(G)，台区(L)和槽台交界处(G-L)三类有限体积单元。

根据流量守恒原理，求解域内为无源流场故有限体积单元内部流量变化为零，任意有限体积单元左侧边界气体流入量Qw与右侧边界气体流出量Qe相等，如式(5)所示。式(6)、式(7)分别为新滑移边界条件下式(4)根据上述网格点划分形式得到的控制方程离散格式，可以看出因为根据分子自由程与气膜间隙之间的关系将气体分子碰撞运动进行分类讨论，进一步联系气膜间克努森数变化对滑移边界的影响，过渡区内气膜间克努森数大于1的方程离散格式参数相对复杂。

(5)

式(5)中：下角标w、e′表示该参量在网格单元中的位置。

当Kn≤1时，有:

(6)

当Kn>1时，有:

(7)

2.3 边界条件与计算过程

阶形板气体流入和流出端的环境压力值设为已知边界条件，沿气体流动方向依次分布的任意有限单元体可构成以Pw、Pc、Pe为未知量的非线性方程，结合两端已知边界条件在求解区域内可写成非线性三对角矩阵方程组的形式，可用Newton-Raphson方法可对其进行数值迭代求解，如式(8)所示。图3为算法流程图。

图3 计算流程图Fig.3 Algorithm flow chart

(8)

式(8)中：n为迭代次数；i表示流场模型中任意网格单元中心。

3 计算结果分析

3.1 模型参数与程序验证

微尺度下的气体流动特性与高空稀薄气体效应存在共性，微尺度效应是气膜间隙的减小导致克努森数的增大，稀薄气体效应是气膜间分子平均自由程的增大导致的克努森数增大。在研究微尺度气体流动特性时可以采用模拟增大环境分子平均自由程直接增大气膜间克努森数的方式进行相似理论研究。直接应用环境克努森数的改变可描述不同模型参数下的微尺度气体流动。阶形板几何模型计算参数取文献[3]，具体数值如表2所示。

表2 模型计算参数Table 2 Model calculation parameters

与文献[3]取相同计算参数，图4(a)中一阶滑移和无滑移边界条件下的阶形板气膜压力分布曲线为应用本文方法所得的数值结果，与文献[3]中的给出的解析解数值基本一致，可作为本文方法可靠性的依据。

图4 气膜压力与克努森数分布Fig.4 Distribution of gas film pressure and Knudsen number

3.2 气膜压力与克努森数

不同环境克努森数Kn0下的阶形板间无量纲气膜压力分布和克努森数分布如图4所示。从图4可以看出，气膜压力在槽台交界处出现最大值，克努森数在槽台交界处出现阶跃变化，阶形板槽区和台区气膜间克努森数跨区域分布。不考虑边界速度滑移时的气膜压力数值明显大于考虑边界速度滑移时的计算数值，一阶滑移边界条件相对于新滑移边界条件明显过高地估计了阶形板间的气膜压力分布数值，过低地估计了克努森数分布数值。在其他模型计算参数保持不变的情况下，随着环境克努森数的增大特别是进入过渡区时，气膜气体可压缩性下降导致气膜压力分布和克努森数分布数值在阶形板台区和槽区的线性化增强，边界速度滑移效应增强且气体分子碰撞频率增大，气膜中气体分子运动矢能降低，气膜压力数值减小的同时气膜间克努森数数值增大。

3.3 气体流速分布

不同环境克努森数下阶形板间槽区和台区中点处沿高度方向的气体流速分布如图5所示，计算结果中的流速分布曲线外形与理论分析一致。从图5可以看出，相较于无滑移边界条件下的流速分布，考虑边界速度滑移时的计算结果在靠近下板处的气体流速滞后且小于下板运动速度，在靠近上板处的气体流速超前且大于静止的上板。随着环境克努森数的增大，特别是进入过渡区后可以看出，边界速度滑移效应的影响更加明显，阶形板槽区和台区靠近上板和下板边界处的气体流速超前和滞后的程度相较于滑流区明显增大。槽区靠近下板处的气体流速受边界速度滑移影响的程度要大于靠近上板处，台区靠近上板处的气体流速受边界速度滑移的影响程度要大于下板处。随环境克努森数的增大不同滑移边界条件下流速分布的线性化程度增强，相同条件参数下新滑移边界条件下的计算结果所呈现出的边界速度滑移程度要大于一阶滑移。

Ux为沿x轴方向的速度图5 气体流速分布Fig.5 Distribution of gas velocity

3.4 阶形板承载力

不同环境克努森数下阶形板无量纲承载力随下板速度参量的变化曲线如图6所示。从图6可以看出，环境克努森数的增大导致阶形板达到最大承载力的临界速度不断升高，相同环境克努森数下新滑移边界条件的临界速度大于一阶滑移。阶形板的最大承载力数值趋于稳定值，未达到临界速度之前边界速度滑移效应明显降低了阶形板承载力， 随着环境克努森数的增大，气膜间更多的分子动能减少表现为更低的承载力数值，一阶滑移相较于新滑移边界条件过高地估计了阶形板的气膜承载力数值。

图6 阶形板承载力Fig.6 Bearing capacity of step plate

4 结论

针对气膜厚度阶跃变化导致克努森数跨区间分布的阶形板结构，引入适用于任意克努森数的新滑边界条件修正雷诺方程，数值模拟了不同环境克努森数下阶形板间的气体流动特性，得到了气膜压力和克努森数分布，气体流速分布和气膜承载力等静态特性的数值结果。得到如下结论。

(1)边界速度滑移效应对阶形板间气体流动特性的影响程度随着Kn0的增大而增大，微观上体现在气膜压力数值的减小和气膜间Kn的增大，靠近上板和下板边界处的气体流速超前与滞后于固体板边界速度，宏观上体现在气膜承载力的减小与临界速度的增大。

(2)Kn0分别为0.05和0.5时，临界速度分别在6 000 mm/s和25 000 mm/s时得到接近的稳定承载力数值，增大下板移动速度可减小滑移边界对阶形板间气膜承载力的影响程度。

(3)气膜间Kn进入过渡区后一阶滑移与新滑移之间产生明显偏差，滑移系数固定的一阶滑移边界条件过低地估计了边界速度滑移效应的影响程度，应用滑移系数随Kn改变的新滑移边界条件可以更加准确地描述阶形板间气体流动。