刘 浩， 须 颖， 王本明， 姚建华， 陶文彬

(广东工业大学 机电工程学院， 广东 广州 510006)

引言

半导体检测设备在芯片制造领域有重要作用，这些设备对运动导轨的精度和速度有着很高的要求。普通机械刚性接触导轨很难实现高速高精度运动。气浮导轨以其非接触、低阻尼的优势，可以实现亚微米甚至纳米级的运动精度，因而成为实现高速高精密运动的主要导轨形式[1-4]。作为气悬浮导轨的核心零件，气浮轴承的合理设计能够显著提高运动台的承载力和精度。气浮轴承的静态特性与气膜厚度有着很大关系，气膜太薄对导轨及轴承面的加工要求太高，气膜太厚则其承载力和静刚度又会显著下降。有效预测气浮轴承在不同气膜厚度下的静态特性有重要工程价值。

气浮轴承的节流方式有小孔节流、表面节流及狭缝节流等。小孔节流气浮轴承具有加工简单、综合性能良好等优势；同时，小孔节流气浮轴承的研究也比较成熟，这类气浮轴承的设计准则已基本建立[5-7]。带腔小孔节流气浮轴承具有较高的承载力和刚度，但是容易发生振动[8-10]。带均压槽结构的气浮支承轴承通过设计均压槽形状可以改善气膜内压力分布[11-12]。狭缝节流气浮轴承具有更好的稳定性[13-14]。

多孔介质气浮轴承有别于以上几种节流方式，其特点在于采用多孔介质材料节流，具有更高的承载力、刚度和较好的稳定性[15-16]。其次，多孔材料表面上有成千上万个出气孔，因而多孔介质气浮轴承对气浮面刮伤不敏感。多孔介质气浮轴承的另一个重要研究对象是多孔介质材料孔隙结构。通常采用渗透率综合评价孔隙结构阻碍介质流动的程度[17-18]。而贯通孔隙度和渗透率有很大的联系，其他孔隙对气体流动没有影响。

但是，上述研究在分析多孔介质渗透率时，多通过二维成像手段获取多孔介质孔隙度，没有考虑三维孔隙连通性的影响。本研究引入micro-CT断层扫描方法获得多孔介质孔隙结构并分离出孤立孔隙和贯通孔隙。理论计算了多孔介质的黏性系数和惯性系数，并进行了试验验证。采用CFD方法进行仿真，研究了多孔介质止推轴承在不同形状、不同气膜厚度下的承载力和刚度，为设计多孔介质气浮轴承提供参考。

1 物性参数分析

1.1 渗透系数

多孔介质前后的压差与介质内部微通道特征有很大的关系，比如孔隙大小、孔隙形状、贯通孔隙度等。气体通过多孔介质流动物理模型示意图如图1所示，浅色为固相，深色为孔隙相。

图1 多孔介质流动示意图

气流从进口流入，沿着贯通孔隙流出，贯通区域壁面拖拽气流。当多孔介质内部气速较低时，流动为层流状态，多孔介质流动方向压差与速度成线性关系，方程为：

(1)

式中， Δp—— 压差

L—— 多孔介质流动方向的长度

μ—— 流动介质的黏度

K—— 渗透率

v—— 多孔介质前端表面气速

(2)

式中,Q—— 多孔介质前端工况体积流量

A—— 多孔介质截面积

当多孔介质压差与速度遵循式(1)时，流动状态为Darcy流；随着流速增大，压差与速度不再遵循式(1)，流动将脱离Darcy状态进入Forchheimer状态，压差与速度的关系可用式(3)描述：

(3)

式中，ρ—— 流动介质密度

C—— 惯性系数

(4)

假设流动是稳态充分发展的，流动主要沿着水平方向。从孔隙尺度上，内部的速度场、压力场比较复杂。Forchheimer方程的推导是基于管道流，因此在推导过程中采用了体积平均的方法。假设在微小控制容积长度dx内，压力和速度是都是均匀分布的。

平均速度可以表示为质量的函数，其中G为质量流量，A为多孔介质截面积。

(5)

考虑到空气的压缩性，引入理想状态方程，压缩空气的密度ρ可以表示为：

(6)

式中，p—— 压缩气体的绝对压力

R —— 特定气体常数

T—— 温度

将式(5)、式(6)带入到式(4)中，然后沿着x方向进行积分，边界条件为：x=0，p=p1；x=L，p=p2，最终得到如下可压缩Forcheimer方程：

f=aG+bG2

(7)

其中，f表示压力的平方在流动方向上的变化率。

(8)

(9)

(10)

实验采用的多孔介质是多孔石墨，材料实验参数如表1所示。

表1 多孔介质参数

实验装置及多孔介质夹具如图2所示，多孔介质密封在金属套管中，用速干胶水T1401密封周边，防止漏气。空气压缩机提供气源，并经过了一系列干燥去油处理；减压阀调节多孔介质前端供气压力；气罐稳定多孔介质前端供气压力；差压表检测多孔介质前后的压差；多孔介质前后放置2个温度传感器用来检测温度变化；流量计为差压式质量流量计，其内部的多层栅格将气流状态调整为层流；为了避免杂质颗粒堵塞栅格，在流量计前端放置了2个杂质过滤器，过滤能力为5 μm。

图2 实验示意图

图3是试验采集到的数据。空气在多孔介质内部的流动主要受到贯通孔隙壁面黏性力和惯性力的作用，这两种作用可以通过黏性渗透系数和惯性渗透系数来表征。多孔介质空气轴承使用的材料渗透率一般介于3.13e-15～8.44e-14 m2，在这个数量级下，多孔介质内部的雷诺数很小，黏性力占主导地位。

图3 质量流量及压差的关系曲线

图4为f和质量流量的关系曲线。可以看出，f和质量流量几乎接近线性的关系，这说明多孔介质内部的气体流动主要是黏性力主导，也就是渗透率系数起到了主要的作用，而惯性渗透率系数可以忽略。采用最小二乘法对实验数据进行线性拟合得到方程为：

图4 f与质量流量的关系曲线

f=9.936e19G2+1.64e17G

(11)

方差为0.999，得到多孔介质的黏性渗透系数为4.7447e-15 m2，惯性渗透系数为312.2390 m-1。

1.2 孔隙度

孔隙度的定义为多孔介质单位体积中孔隙所占的比例。孔隙度进一步细化可以分为贯通孔隙度和孤立孔隙度，如图5所示。贯通孔隙在流动方向上是贯通的；孤立孔隙指的是多孔介质中闭死的孔隙，与外部大气没有任何连通通道。

图5 孔隙结构细分示意图

文献[20-21]计算孔隙度多采用SEM图像法、压汞/压水法、煮沸法等。SEM图像分析对材料表面二维图像进行阈值分割出孔隙，进而得到面孔隙度，将面孔隙度视为三维空间上的孔隙度。这种方法没有考虑孔隙在三维结构上的连通性。

因此，本研究采用micro-CT计算机断层扫描方法分离多孔介质的贯通孔隙度，其扫描原理如图6所示。射线源发射的X光穿透样品，部分射线被样品所吸收，剩余量会打在探测板上被感应元件探测到。样品旋转360°，探测板共采集数百至数千帧数据，之后数据被重构，从而还原出样品的三维结构。试验数据采集使用的设备为天津某公司制造的nanoVoxel 3000型X射线三维显微镜。

图6 micro-CT扫描原理示意图

三维多孔介质图像能够显示孔隙在三维方向上的连通关系，图像后处理进行中值滤波降低噪声，改善图片质量。采用分水岭算法分割孔隙，通过逻辑运算将孔隙结构细分为贯通孔隙、孤立孔隙。图7a为贯通孔隙，图7b为孤立孔隙。

图7 孔隙结构图

多孔介质孔隙度分析的关键之一在于找到其最小重复单元(Representative Elementary Volume,REV)。REV参考BEAR[22]提出的定义：不断增大截取的子样本体积并计算对应的孔隙度，当孔隙度不再随着子样本体积V增大而变化时，认为这个体积为REV。图8显示了REV寻找的曲线图，当体积为0.25 mm3时孔隙度φ基本保持恒定。

图8 样本最小重复单元及孔隙度的关系曲线

2 静态特性分析

多孔介质的阻尼作用以沉积项的形式附加在标准控制方程。因此，在计算流场时有如下假设：

(1) 多孔介质内部的堵塞效应不会被考虑；

(2) 多孔介质内部的流动假设为层流；

(3) 多孔介质的孔隙度是各向同性的；

(4) 不考虑多孔介质与激波的相互作用。

单相空气流动的体积平均质量和动量守恒方程如下所示：

(12)

(13)

热平衡方程为：

(14)

式中，U—— 速度向量

I—— 单位矩阵

S—— 多孔介质源项

γ—— 贯通孔隙度

Ef—— 总流体能量

Es—— 总固体能量

ρf—— 流体密度

ρs—— 固体密度

keff—— 有效热传导率

2.1 网格无关性验证

采用全结构六面体网格对多孔介质空气轴承流场特性进行数值计算。设计了3种相同的承载面积和介质厚度不同物理形状的多孔介质空气轴承，网格模型如图9所示。在气膜厚度为10 μm时，进行网格无关性验证，验证结果如图10所示，其中横坐标表示网格数量n，纵坐标表示承载力F。可以看出，多孔介质空气轴承总网格数在8×104左右，继续增加网格数量对结果几乎没有影响。

图9 不同形状多孔介质空气轴承网格

图10 网格数量及承载力的关系曲线

2.2 静态特性分析

静态特性包含承载力和静刚度。静刚度指在外部载荷变化过程中，不考虑气膜的动态变化，单位气膜厚度变化所需要的力。通常，静刚度越高则运动部分抗载荷变化的能力越强，保持高精度的能力越强。

流体控制方程为复杂的偏微分方程组，获得方程组的分析解是很困难的。数值解能够在一定的精度下为工程应用提供有效参考价值。数值模拟通常需要模型网格化、方程组离散化等一系列操作。其中，网格的质量对迭代计算是否收敛影响很大，各个形状网格的质量如表2所示。

表2 网格质量

基于正交模拟的原则，除了结构形状外，其他物理边界条件保持相同，详细数值模拟参数如表3所示。

表3 CFD模拟条件

由图11可以看出，随着气膜厚度l的减小多孔介质空气轴承的承载力F一直增加。这是由于随着气膜厚度的减小，气膜腔的表面积体积值之比增大，阻尼增大，气膜内的平均压力升高，空气的弹性模量增大，进一步压缩气膜需要更大的力。

图11 气膜厚度及承载力的关系曲线

由图12可以看出，随着气膜厚度的减小，静刚度会出现1峰值，之后气膜厚度的再次减小导致静刚度下降。不同形状的多孔介质空气轴承静刚度峰值及对应的气膜厚度也有略有差异。其中，圆形多孔介质气膜厚度为6 μm时静刚度峰值为82 N/μm，方形多孔介质在气膜厚度为5 μm时静刚度峰值为82 N/μm，三角形多孔介质气膜厚度为5 μm时静刚度峰值为88 N/μm。相同的截面积情况下，三角形的静态特性值略高一些。合理设计多孔介质截面的大小以及气膜的厚度能够很好均衡承载能力和静刚度。

图12 气膜厚度及刚度的关系曲线

为了进一步理解相同面积不同形状多孔介质气浮轴承的承载差异，采用有限体积法计算了气膜内的压力场，CFD模拟参数如表3所示，气膜厚度选定为5 μm，压力云图如图13所示。

图13 气膜内压力场力

沿着图13中3条对称线绘制压力值曲线，如图14所示,横坐标表示对称线长度x，纵坐标表示对应点静压力p。可以看出，圆形的高压区压力高且保持线最长，正方形高压区压力保持较高但保持线短，三角形的高压区压力相对不大。由此可见，高压区压力大小和面积决定了承载力的大小。

图14 对称线上压力分布

3 结论

完整提出了多孔介质静态特性预测方法，从多孔介质的物性参数探测、贯通孔隙度计算到基于有限体积的静态特性数值模拟，得到以下结论：

(1) 采用可压缩Forchheimer方程获得多孔介质黏性系数和惯性系数；

(2) 采用micro-CT方法得到多孔介质的最小重复单元，体素在1.48 μm，样本REV尺寸为0.25 mm3；

(3) 多孔介质空气轴承的承载能力随着气膜的减小而增大，刚度有1个峰值，圆形多孔介质峰值刚度为82 N/μm在6 μm处，方形为82 N/μm在5 μm处，三角形为88 N/μm在5 μm处，相同的截面积情况下，三角形的刚度较高。