基于三维有限元波箔片模型的气体箔片轴承承载性能研究

梁 波，罗欣洋，宁家强，程定斌，徐科繁，张广辉，刘占生

(1.中国航空工业集团公司 南京机电液压工程研究中心航空机电系统综合航空科技重点实验室,江苏 南京 211106;2.哈尔滨工业大学 能源科学与工程学院,黑龙江 哈尔滨 150001)

箔片空气动压轴承(Air Foil Bearing)是通过弹性箔片将转轴包裹，由转子高速旋转产生的动压效应将箔片撑开，进而形成气膜使转子浮起的一种轴承形式。箔片空气动压轴承主要有悬臂型、缠绕型、波箔型和外楔型等形式[1]。图1为典型的波箔型径向空气箔片轴承。

图1 典型的空气箔片轴承结构形式[2]

箔片空气动压轴承在小型高速旋转机械中具有广阔的应用前景，尤其是在飞机空气循环机中，此时机载环境恶劣，转速变化剧烈，而箔片空气动压轴承既能适应高速运转条件下的稳定运行，又具有可靠性高、长寿命、无润滑油污染、温度适应性宽等优点[3]，因此已经取得广泛的应用[4]。

戚社苗等[5]借助MATLAB PDE工具箱对雷诺方程进行了高精度计算。徐方程等[6]考虑了箔片间的摩擦效应，结果表明摩擦效应增大了波箔的刚度。李洁[7]基于曲壳模型对箔片轴承的承载特性进行了分析，通过商业软件ANSYS模拟以及实验验证了其准确性。许浩杰等[8]研究了箔片结构参数对箔片轴承承载性能的影响，结果表明当偏心率一定时，轴承承载力与波箔厚度成正相关。Gu等[9]采用包含分离和摩擦行为的通用接触模型对平箔和波箔的相互作用进行了建模，结果表明该模型允许平箔与波箔分离，符合工程实际。

综上所述，箔片轴承的承载性能受箔片弹性变形影响较大，而后者与箔片结构参数密切相关，这一点有待进一步研究。本文以波箔型空气动压轴承为研究对象，结合有限元法建立了两种箔片变形模型，通过耦合求解雷诺方程以及箔片变形模型，获得了不同波箔片弹性和不同平箔片厚度对轴承承载力的影响，为后续波箔型空气动压轴承的设计制造提高理论参考。

1 雷诺方程

箔片空气动压轴承的控制方程是雷诺方程，其无量纲形式参见式(1)。雷诺方程是非线性的二维椭圆型方程，可通过有限差分法[10-16]或有限元法[17]求解。本文选择前者。并结合五点差分格式求解。此时轴向为周期边界条件，周向为大气边界条件

(1)

为提高程序在大偏心率时的收敛性，本文采用亚松弛法(SUR)[13]获取每一次迭代的初值，见式(2)，此时亚松弛系数为0.5

(2)

2 箔片变形模型

2.1 基于三维有限元法的箔片变形模型

文献中动压轴承的平箔片变形模型主要有以下三种：(1)忽略平箔片并假定波箔片与弹簧等价[8]；(2)考虑平箔片沿着周向的刚度[18]；(3)考虑平箔片与波箔片的摩擦力[11]。本文采用前两种模型计算箔片变形。

针对忽略平箔片并假定波箔片与弹簧等价的模型，此时气膜压力直接作用在弹簧上，参见图2(a)。

针对考虑平箔片沿着周向的刚度的模型，此时平箔片将简化为一系列固结在弹簧上的三节点梁单元，箔片即简化为弹簧与梁单元的组合体，参见图2(b)。

图2 两种箔片变形模型

2.2 箔片变形模型1的建模

针对模型(1)，每个节点上气膜压力产生的变形量按式(3)计算

δij=fij/K

(3)

式中δij——节点的变形量；

K——波箔片的弹性系数。

本文通过有限元软件ABAQUS进行仿真计算以获取弹性系数。物理模型参见图3，其厚度为0.1 mm，弹性模量为E=2.0×1011Pa，泊松比υ=0.3。Abaqus中波箔片结构模型和计算结果参见图4。

图3 典型波箔片的结构参数/mm

图4 ABAQUS仿真结果

根据仿真结果，波箔片的弹性系数为K=4.167 GN/m3。文献中利用经验公式获取的波箔片弹性系数依次为4.7 GN/m3和5.0 GN/m3，显然本文的计算结果可信。

为进一步评估不同载荷作用下波箔片的弹性系数变化，本文绘制了载荷与波箔片最大位移的关系，如图5所示。可以看到，波箔片的弹性系数在一个较大范围内保持不变，因此，模型1可以假定波箔片与弹簧等价。

图5 不同载荷下的波箔片的最大位移

在此基础上，分别改变波宽和波高，然后计算波箔片的弹性系数，相关结果如图6所示。可以看到，波箔片的弹性系数随着波宽的增加逐渐降低，而波高在一定范围内影响较小。显然，波宽是一个重要设计参数，可以通过修改波宽改变波箔片的弹性系数。

图6 波箔片弹性系数与波宽、波高的关系

2.3 箔片变形模型2的建模

针对模型(2)，由于此时梁单元仅有三个节点，而二次分布有三个未知数，因此可以假设二次分布能逼近箔片单元所受的压力分布qi,j(ζ)，故有

qi,j(ζ)=(ai,jζ2+bi,jζ+ci,j)Δz

(4)

结合图2(b)有

(5)

对公式(5)求解有

(6)

在此基础上，由材料力学中梁单元变形的微分方程[19]有

(7)

式中ωi,j——箔片单元的变形量分布。

(8)

(9)

结合式(4)、(6)、(8)、(9)可获取平箔片的变形量分布，参见式(10)

(10)

需要注意的是，式(10)的前提是假设平箔片单元的两端是固支边界条件。当同时考虑波箔片的弹性效应时，平箔片单元的两端与轴瓦间是弹性连接。因此，需要结合式(3)获取平箔片单元的两端点的变形量δij。此时fi,j通过式(11)求解

(11)

以si,j表示平箔片总的变形量，显然其等于两端点变形量的线性插值结果与平箔片的变形量之和，参见式(12)，进而有(13)

(12)

(13)

3 数值计算模型

3.1 耦合计算流程

数值计算过程中，若箔片的刚度较小，气膜压力会使箔片产生较大的变形，最终容易造成计算发散，这个问题在大偏心率情况下尤其明显。为提高程序收敛性和稳定性，本文提出了小柔度步数值方法，如图7所示。

图7 气膜与弹性箔片的耦合计算流程图

小柔度步数值方法将计算过程分为多个柔度步。先采用小柔度计算使得气膜压力和气膜间隙满足收敛条件，然后逐渐加大柔度值，最终得到全柔度下气膜压力和气膜间隙的求解结果。经过验证，小柔度步数值方法可以显著提高程序的稳定性。

3.2 程序验证

为验证程序的准确性，本文按照文献[15]中的参数对动压轴承进行模拟，结果如图8所示。图9给

图8 按照文献[11]参数计算的气膜压力分布

图9 气膜压力与气膜间隙对比

出了文献[15]和本文计算结果的对比图。可以看到，中截面上气膜压力气膜压力偏差较小，误差小于3%，这说明本文程序可以准确模拟箔片空气动压轴承。

4 箔片弹性对承载力的影响分析

4.1 箔片变形模型1：波箔片弹性的影响

结合箔片变形模型(1)，本文按照3.1给定的计算流程，采用表1中的计算参数，对动压轴承进行了模拟，结果如图10所示。

表1 弹性箔片轴承计算参数

图10 箔片轴承气膜压力分布和气膜间隙分布

由图10可知，气膜压力较大处气膜间隙较小，同时由于大气边界条件的影响，边界点附近箔片变形量小，气膜间隙更小。

在此基础上，为更直观地比较箔片轴承与刚性轴承，图11给出了中截面上轴承气膜压力和气膜间隙分布。同时考虑到决定动压轴承承载能力的主要因素之一是最小气膜间隙，故此时选定刚性轴承偏心率为0.55，箔片轴承偏心率为0.8。

图11 相同最小气膜间隙下刚性轴承与箔片轴承的中截面气膜压力分布与气膜间隙分布

可以看到，当最小气膜间隙相同时，箔片轴承的承载区压力更大，整个承载区更宽，同时气膜间隙形状更陡，显然能提供更大的承载能力。

表2给出了两种轴承的承载能力[14]对比数据。可以看到，当最小气膜间隙相同时，箔片轴承的承载能力约为刚性轴承的1.5倍，同时其偏位角也较小，一定程序上具有更好的气膜稳定性。

表2 相同最小间隙下弹性箔片轴承与刚性轴承无量纲承载能力对比

为进一步分析箔片弹性对轴承承载力的影响，图12给出了不同箔片刚度下，箔片轴承的气膜压力分布和气膜间隙分布。

图12 不同刚度下箔片轴承中截面的气膜压力分布和气膜间隙分布

随着箔片刚度逐渐减小，轴承承载区宽度逐渐增大，气膜的最大压力逐渐增加。但如果将箔片刚度减小到一定数值时，如K=2.582 GN/m3，由于箔片刚度过小，轴承无法形成稳定的气膜，气膜压力会呈波动状态，参见图12，此时程序未收敛。

通过上述分析，显然，在动压轴承设计时，应该合理选择箔片刚度，以保证形成较高承载能力的气膜，但又不至于产生气膜震荡。

4.2 箔片变形模型2：平箔片厚度的影响

结合箔片变形模型(2)和计算参数表1，图13给出了相关的计算结果，此时平箔片厚度为0.1 mm。

图13 考虑平箔片变形的气膜压力分布与间隙分布

图14给出了两种箔片变形模型对应的数值模拟结果。可以看到，考虑平箔片变形后，在压力较大处，由于梁单元的变形，整个间隙分布会有微小的波动。另外，轴承的气膜间隙分布会沿圆周方向钝化，这是因为箔片变形模型2考虑了沿圆周方向的箔片刚度。

图14 两种模型下的中截面气膜压力分布和气膜间隙分布的比较

由节2.2、2.3可知，在箔片变形模型1中，节点变形量仅与该节点的压力有关；在箔片变形模型2中。节点变形量由相邻的沿圆周方向的几个节点共同决定，显然，后者更加准确。

图15给出了不同平箔片厚度条件下轴承气膜厚度和气膜间隙的分布情况。

图15 不同平箔片厚度下中截面气膜压力分布和气膜间隙分布的比较

可以看到，当平箔片厚度不同时，轴承气膜压力和气膜间隙分布基本一致。但当平箔片厚度降低到0.08 mm时，气膜间隙在高承载区有更多的波动，这是由于平箔片厚度降低时，平箔片的截面惯性矩会以三次方的速度降低，进而导致更大的变形量。另外，若将平箔片厚度增大到0.15 mm，此时波动情况会得到弱化。

5 结论

本文通过对雷诺方程以及箔片变形模型的耦合求解，获取了箔片弹性对空气动压轴承承载力的影响，同时数值结果与文献数据误差小于3%。相关结论如下：

(1)本文提出了小柔度步的数值方法，显著提高了程序的收敛性和稳定性。

(2)针对箔片变形模型1，相较于刚性轴承，波箔片的弹性会提高箔片轴承的承载区宽度，在同样的最小气膜间隙下，轴承承载能力会显著增强。

(3)波箔片的刚度应该合理设计，以形成承载能力较大的气膜形状并不引起气膜震荡。

(4)针对箔片变形模型2，平箔片会增加箔片的弹性效果，相较于箔片变形模型1，箔片轴承的间隙分布会在圆周方向钝化。

(5)平箔片厚度一般选择0.1 mm，过小的平箔片厚度会引起气膜震荡。