基于Workbench的高压圆盘气体轴承三维流场分析*

郭良斌 吴永良

(1.武汉科技大学冶金装备及其控制教育部重点实验室 湖北武汉 430081；2.武汉科技大学机械传动与制造工程湖北省重点实验室 湖北武汉 430081)

气体轴承因其高精度、长寿命和低功耗等特点，已被广泛运用于纺织机械、仪表机床等现代工业领域[1]。传统的静压气体轴承刚度低、承载力小、供气压力低，只适用于轻载精密支撑[2]。为此，研究人员对气体轴承的发展与应用进行了深入研究[3-4]。徐凡[5]参考风洞实验中有关收缩段的设计理论，设计出气体轴承新型双对称收缩段壁面型线，使气体轴承内部流场得到有效改善，消除了供气压力的限制。刘剑和郭良斌[6]对高压圆盘气体轴承可压缩边界层形态进行了详细分析，研究了气膜厚度、轴承半径以及供气压力对气膜边界层发展特征的影响。郭良斌和刘剑[7]对不同供气总压时高压圆盘气体轴承边界层的对流换热系数进行了进一步研究，通过对气膜内的流场进行数值模拟，分析了速度边界层的特征，并预测了气膜内边界层的转捩点。然而，因为该高压圆盘气体轴承具有双对称性，上述研究均是在双对称条件下取模型的1/2或1/4进行简化，从而得到具有对称面和对称轴的二维简化模型及其流场特性。目前，国内外已有不少学者对气体轴承的三维流场特性进行了研究[8-14]，对未来研究气体轴承特性具有很好的参考价值。由于高压圆盘气体轴承模型的特殊性，以及气膜间隙内湍流的不确定性，对该模型进行三维流场分析显得尤为重要。

为了研究高压圆盘气体轴承完整模型的三维流场特性，本文作者采用Ansys17.0仿真软件内的Workbench工作平台，利用Solidworks专业建模软件对高压圆盘气体轴承进行三维建模，采用Workbench自带的网格划分工具Mesh对几何模型进行网格划分，并利用Fluent流体分析模块进行计算，最终得到了高压圆盘气体轴承的三维流场，并与二维简化模型流场计算结果进行了对比。

1 轴承结构原理与内流域抽取

1.1 轴承结构原理

由于文中主要研究的是高压圆盘气体轴承内流场状况，故对轴承本身建模作简化处理，忽略稳流管与工作圆盘之间的密封垫以及工作圆盘上、下端面的倒角。高压圆盘气体轴承的三维模型结构如图1所示。

由图1可知，完整的高压圆盘气体轴承由上、下圆盘组件构成；上圆盘组件包括上稳流管和上工作圆盘，二者通过螺纹连接紧固在一起，上圆盘组件与下圆盘组件完全对称。稳流管长72 mm，工作圆盘厚30 mm，两工作圆盘之间的工作气膜间隙为0.2 mm，空气由上、下稳流管的供气孔进入，经气膜间隙从两工作圆盘间的外边缘射流流出。由于采用Workbench中的Fill功能进行内流域抽取，该功能需保持入口和出口均为封闭状态。为了使整个计算域包含气膜出口缝隙射流区域，将图1的原始三维模型修改为如图2所示的结构。

图1 高压圆盘气体轴承三维结构

图2 修改后的高压圆盘气体轴承三维结构和剖面

即在轴承气膜出口添加轴向长度20.2 mm、径向长度15 mm的三维域作为气体出口的远场，用薄壳包裹，壳厚2 mm。

1.2 内流域抽取

将图2模型导入Workbench中的Fluid Flow(Fluent)分析系统，利用Workbench自带的Design Modeler建模模块的Fill功能对该轴承内流场进行流体域抽取，再通过压缩轴承实体得到轴承内流道的几何模型，如图3所示。

图3 轴承内流道几何模型

2 网格划分与边界条件设置

2.1 网格划分

如图4所示，Workbench中的Fluid Flow(Fluent)分析系统会自动保存Design Modeler模块中的建模，并将其导入Mesh网格划分模块。采用的网格划分方式为自动网格划分，通过全局网格参数设置，将网格精度relevance、网格细化程度relevance center以及跨度中心角span angle center调至最优，网格平滑度smoothing调至中等；选择尺寸函数类型为curvature，系统会根据曲率法向角度自动确定细化边与曲面处的网格单元尺寸大小，设置曲率法向角度curvature normal angle为30°，最小网格尺寸为0.05 mm，最大面部尺寸为1.2 mm，最大四面体尺寸为1.2 mm，网格增长率为1.05，其余保持默认设置。最终划分网格平均质量为0.83，平均Aspect ratio为1.9，网格节点数约30万个，网格数量约160万个。

图4 轴承整体网格与局部网格

2.2 网格无关性验证

一般而言，模型划分网格数量越多，其计算结果精度越高，但过高的网格数量，会延长计算时间，降低工作效率，甚至普通的计算机难以进行正常工况计算。因此划分合适的网格数量显得尤为重要。

如图5所示为83万～207万网格数量下的轴承气膜对称线上马赫数的分布状况。可知，不同网格数量下在轴承入口和收缩段部分(72.39～132.39 mm)以及平行气膜区(42.39～72.39 mm、132.39～162.39 mm)其马赫数曲线基本重合，仅在收缩段出口处(72.39、132.39 mm)以及缝隙射流区(24～42.39 mm、162.39～178 mm)等处出现小幅变化。但马赫数总体变化趋势是一致的，且160万网格数量时在轴承气膜出口处能更好反映出马赫数最大值。考虑到计算机配置及程序运行等情况，文中采用160万网格数量的三维划分模型进行工况计算。

图5 网格无关性验证结果

2.3 边界条件设置

如图6所示，轴承圆盘直径为120.54 mm，稳流管气流入口直径为20 mm，供气孔直径为10 mm。将上、下稳流管入口面设置为压力入口IN，计算域与远场外壁的接触面设置为压力出口OUT，压力出口OUT轴向宽度20.2 mm，其余面WALL-WL、WALL-BEARING、WALL-ZZCM、WALL-FARFIELD均设为壁面。入口总压设置为2 MPa(20 atm)，初始表压设置为1.998 MPa(19.98 atm)。湍流项选择湍流强度和水力直径，其中湍流强度设为7%，入口水力直径为20 mm，入口总温为300.15 K。出口表压设置为0.1 MPa(1 atm)，湍流强度设为7%，出口水力直径为150.54 mm，总温为300.15 K。其中，完全发展的管流的核心湍流强度用以下经验公式计算：

图6 轴承各边界命名

(1)

式中：I为湍流强度；u′为湍流脉动速度；uavg为平均速度；ReDH为按水力直径计算得到的雷诺数，圆管水力直径即为其自身直径D：

DH=D

(2)

所有壁面均为固定无滑移恒温壁面，壁温300.15 K。

3 控制方程与湍流模型

3.1 控制方程

流体处于稳态流动时，其三维可压缩的、笛卡尔坐标系下的控制方程[7]为

(3)

(4)

(5)

(6)

3.2 湍流模型

完整的控制方程需要引入湍流模型将其封闭，建立雷诺应力方程组。文中数值计算使用Realizablek-ε模型。 Realizablek-ε是近期出现的带旋流修正的k-ε模型，对比标准的k-ε模型，其增加了一个与湍流黏性有关的公式，为湍流耗散率增加了新的传输方程，该模型直接的好处就是对于平板和圆柱射流的发散比率有着更精确的预测。Realizablek-ε模型的方程如下：

Gb-ρε-YM+Sk

(7)

(8)

4 数值计算与结果分析

4.1 计算参数设置

采用Workbench中的Fluent模块，对模型进行求解。选取基于压力的绝对稳态求解器，操作压力设为0，开启能量方程与湍流方程，设置流入气体为空气，并假设其为连续的理想气体，其密度采用理想气体状态方程，黏度采用萨瑟兰公式：

(9)

式中：T0=273.15 K；μ0=1.716×10-5Pa·s；Ts=110.56 K。空气导热系数、比定压热容以及分子质量设为默认常数，湍流模型使用标准壁面函数。求解方法采用SIMPLEC算法，压力、密度以及湍流动能等差分格式均保持默认设置，各项欠松弛因子、各项残差标准也保持默认设置。

4.2 结果及分析

4.2.1 三维气膜流场静压分析

图7(a)所示为气膜流场静压云图。可知，轴承气流入口压力最高，气流到达收缩段前压力下降并不明显；从收缩段开始，压力以稳流管轴线为中心，轴对称向四周逐渐降低，但整个轴承气膜内压力仍维持在0.5 MPa(4.93 atm)以上，仅在气膜出口处很窄的范围内压力快速下降，达到了良好的高压润滑效果；气流经气膜间隙出口流出后，气流逐渐降至略大于环境压力。利用Fluent软件的Plane Tool工具新建气膜对称面，再将轴承圆盘及稳流管的轴对称面与其相交，截得气膜对称线。图7(b)所示为气膜对称线上静压分布曲线。可知，静压分布曲线以R=102.39 mm的铅垂线为对称轴左右对称，以右半边曲线为例，气流压力在供气孔区域a(102.39～110 mm)下降缓慢，而在收缩段区域b(110～138.9 mm)气流压力快速下降至0.89 MPa(8.9 atm)，而后气流到达平行气膜区c(138.9～62.39 mm)，压力下降至0.5 MPa(4.93 atm)。在气膜出口处很窄范围内压力迅速降至0.11 MPa(1.1 atm)，最后在缝隙射流区d(162.39～178 mm)压力维持在略大于环境压力。

图7 气膜流场静压云图及气膜对称线上静压分布曲线

4.2.2 三维气膜流场马赫数分析

图8(a)所示为气膜流场马赫数云图。可知，稳流管气流入口处马赫数最低，气流到达收缩段前马赫数上升不明显，马赫数不超过；从收缩段开始，马赫数以稳流管轴线为中心，轴对称向四周逐渐升高，可以观察到，位于收缩段终点的马赫数已升至高亚音速，进而在平行气膜区有小幅度降低；而后在气膜出口处因气流压力远大于环境背压，使气流速度上升至音速，随后，在气膜出口外进一步形成超音速圆盘缝隙射流，最终因膨胀波和压缩波的交替作用马赫数逐渐趋于0，气流逐渐静止。图8(b)所示为气膜对称线上的马赫数分布曲线。可见，马赫数曲线以R=102.39 mm的铅垂线为对称轴左右对称，取右半边曲线为例，气流在供气孔区域a(102.39～110 mm)的马赫数为0.2，经过收缩段b(110～138.9 mm)后，马赫数上升至0.92，且在平行气膜区c(138.9～162.39 mm)的绝大部分区域，气流都维持在马赫数0.7～1的高亚音速状态。在气膜间隙出口处，由于流通截面的突然增大，气流因膨胀波的作用升至超音速马赫数1.36，经过缝隙射流区d(162.39～178 mm)后，马赫数逐渐趋于0。

图8 气膜流场马赫数云图及气膜对称线上马赫数分布曲线

4.2.3 三维模型与二维简化模型的计算结果对比

4.2.3.1 二维简化模型

二维模型及其边界条件如图9(a)所示，即OUT为压力出口，IN为压力入口，SYM为对称面，AXIS为对称轴，其余WALL-WL、WALL-BEARING、WALL-ZZCM、WALL-FARFIELD均为壁面。采用ICEM对二维模型进行网格划分，并对边界层进行局部加密，控制y+值在2.0以内，最终划分网格平均质量为0.96，平均Aspect ratio为1.98，网格数量约90万个，收缩段入口网格划分如图9(b)所示。利用Fluent求解二维流场。入口总压为2 MPa(20 atm)，初始表压为1.998 MPa(19.98 atm)，湍流项选择湍流强度和水力直径，其中湍流强度为7%，入口水力直径为20 mm，入口总温为300.15 K；流体出口设置为压力出口，出口表压为0.1 MPa(1 atm)，湍流强度为7%，出口水力直径为120 mm，总温为300.15 K；采用SIMPLE算法，其余设置同三维模型。

图9 二维模型边界条件及收缩段入口网格划分

4.2.3.2 二维模型网格无关性验证

如图10所示为网格数量在65万～97万时二维模型对称面上的马赫数分布状况。

图10 二维网格无关性验证结果

由图10可知，随着网格数量的增加，轴承气膜出口处(60～62 mm)最大马赫数会发生小幅变化，而不同网格数量时马赫数整体变化趋势是一样的，其马赫数曲线基本重合。考虑到计算机配置及程序运行等情况，文中采用90万网格数量的二维划分模型进行工况计算。

4.2.3.3 2种模型的计算结果比较

利用Fluent自带截面工具截取气膜对称线，并获取线上的马赫数与静压数据。因为二维简化模型流场的气膜对称线只有整个气膜的1/2，为了方便对比，使用了MatLab软件编程，对二维计算数据进行对称扩展，得到整个气膜的对称线上的数据，在origin中绘图进行对比，如图11所示。

由图11(a)可知，二维简化模型与三维模型的静压变化趋势基本一致，且变化的位置也基本一致，但在平行气膜区开始部分(130～140 mm)，二维简化模型静压曲线压力明显高于三维模型静压曲线的压力。其原因在于三维模型使用Workbench自带的mesh模块进行网格划分，其网格划分能力稍差，难以添加高质量边界层网格，对于高压圆盘气体轴承的微小气膜间隙(0.2 mm)难以划出较好的结构化网格，整体网格质量一般(0.83左右)；而二维简化模型使用ICEM专业网格划分软件，网格划分能力强，可添加边界层并合理控制y+值，网格质量较好(0.96以上)，因而造成了计算结果上的小幅偏差，以及出现了不一样的极值点。因此，二维模型网格质量比三维模型网格质量更高，其最大值显示更为精确。

图11 2种模型气膜对称线上静压和马赫数对比

由图11(b)可知，二维简化模型与三维模型的马赫数变化趋势基本一致，且变化的位置也基本一致，但在气膜间隙出口处，二维简化模型马赫数曲线的最大值明显高于三维模型马赫数曲线的最大值，其原因仍是二维网格质量更高使其最大值更为精确；二维模型缝隙射流区较窄(70 mm-58.9 mm=11.1 mm)，三维模型缝隙射流区较宽(178 mm-165 mm=13 mm)，且缝隙射流区(165～178 mm)三维模型马赫数曲线明显高于二维简化模型马赫数曲线，其原因在于三维模型中气体流出气膜间隙后，在缝隙射流区不仅有径向的流动，还存在周向发散的流动，这是二维简化模型难以预测到的现象。

4.2.4 不同湍流强度对高压圆盘气体轴承静压和马赫数的影响

虽然湍流强度理论上按公式(1)进行计算，但实际上通过管道供气，其湍流强度可能会发生改变。研究不同湍流强度下高压圆盘气体轴承静压和马赫数的变化情况，对后续实验有着重要指导作用。不同湍流强度下高压圆盘气体轴承静压和马赫数分布情况如图12所示。

图12 不同湍流强度下气体轴承静压和马赫数分布情况

由图12(a)可知，不同湍流强度时供气孔区域a、收缩段区域b、平行气膜区c以及缝隙射流区d的静压分布曲线基本重合，可见不同湍流强度对高压圆盘气体轴承内气膜压力无明显影响。其原因在于湍流强度和轴承气膜压力之间并无直接关联。

由图12(b)可知，不同湍流强度时供气孔区域a、收缩段区域b以及平行气膜区c的马赫数分布曲线基本重合，而在缝隙射流区d及气流远场区域(12.39～26 mm、178～192.39 mm)马赫数存在较明显的变化。其原因在于气流经平行气膜出口流出时加速为超音速，且湍流流动存在不确定性，湍流脉动速度在不断发生变化。

5 结论

(1)高压圆盘气体轴承内，气膜压力随着半径方向R的增大而逐渐减小。平行气膜区绝大部分区域气膜压力在0.5 MPa(4.93 atm)以上，仅在气膜出口处很窄的范围内压力快速下降。经超音速圆盘缝隙射流流出后，气流压力逐渐降至略高于环境压力。由于超音速流动时外部环境低压不影响上游压力分布，故文中的高压圆盘气体轴承可以实现良好的高压润滑作用。

(2)就整体变化趋势而言，二维简化模型与三维模型计算结果基本保持一致，出现小范围差别的主要原因在于，三维模型使用Workbench自带的mesh模块进行网格划分，其网格划分能力稍差，难以添加高质量边界层网格。但二者马赫数、静压曲线变化趋势基本一致，较好地反映了高压圆盘气体轴承内部流道的流动规律。

(3)不同湍流强度对高压圆盘气体轴承内流场的静压无明显影响，但对缝隙射流区及气流远场区域的马赫数有着小幅影响。因此后续实验中，测量气膜内流动参数时，可以排除供气管道布置方式改变造成入口湍流强度变化这一影响因素。