有限差分法求动静压轴颈轴承性能研究

王学敏，白红宇，庄明

(中国科学院等离子体物理研究所，安徽合肥230031)

由于气体润滑介质具有摩擦功率小、无污染、运动精度高等优点，特别是近年来随着核能和航空航天技术的发展，其应用前景广阔，通过求解Reynolds方程准确预测轴承性能是气体轴承应用的关键.目前，已经有多种数值方法应用于对Reynolds方程的求解［1］，其中由于有限差分法编程序灵活、成熟，其应用最为广泛.国内外虽然有多篇关于利用有限差分法求气体轴承性能的文章［2-4］，但给出的大都是未考虑偏位角影响的静态算例，在求解过程中对一些影响计算精度的因素也并未加以讨论，本文对用有限差分求解轴颈轴承性能过程中若干影响计算精度的问题进行了细致讨论，并得出了若干可以提高计算结果精度的结论.

1 Reynolds方程无量纲化及边界条件

很多文献中都有详细的Reynolds方程推导过程［5-6］，完整的非定常可压缩雷诺方程如下:

以上代入式(1)，得到稳态无量纲化Reynolds方程式如下:

式中:H=1+εcos θ，ε 为偏心率，式(2)的无量纲边界条件为:

2)气膜被剪开的两边界上:P(θ=0)=P(θ =2π)，即 P1，j=PΠ，j，其中 Π 为最大网格节点数，dP/dx(θ=0)=dP/dx(θ=2π)，即(P2，j-P1，j)/ΔX=(PΠ，j-PΠ-1，j)/ΔX.

2 Reynolds方程离散化

由于气膜厚度相对于其他几何尺寸小103～104倍，因此可以略去圆柱表面率的影响，为了便于程序编写，一般都是把圆柱轴承的气膜展开成如图1的一个平面［7-8］，将轴承面按两供气孔的中心线划分为N个区域，其中N为供气孔数(本文以双排8个供气孔为例).

图1 气膜划分区域示意图Fig.1 Computational domain of journal gas bearings

采用具有二阶精度的中心差分格式的超松弛迭代法(SOR)［9］对式(2)进行离散，得到离散化Reynolds方程式为

单个区域流出的质量流量包括:

轴承承载能力计算为

式中，Wn和Wt分别表示气膜内压力在偏心方向上对面积积分和与其成直角方向上对面积积分.

3 有限差分法求解过程若干问题

为了比较算法的有效性，本文利用具有解析解的ε=0时轴颈轴承模型，其模型如图2，具体参数如下:轴承直径D=14 mm，长度L=33 mm，两排环面节流供气小孔，每排8个，小孔直径d=0.4 mm，小孔至端面距离L1=9.5 mm，轴承平均半径间隙Cr=0.022 5 mm，偏心率ε=0，气体工质为氦气.

图2 ε=0时理想模型Fig.2 Ideal gas bearing model with ε =0

理想情况下，当满足收敛条件时，气膜内压力分布应该左右上下完全对称，承载能力=0.但实际情况会有以下原因造成载荷不为0，即存在一个小的载荷误差Wσ，正常情况下，载荷误差Wσ相比有偏心情况下的载荷为小值，不会对结果产生本质影响.

1)边界条件2)假设被剪开两条边上的压力其与相邻的两条边上压力相等，当网格数不足够多时，相邻两条边总会有一定压力差，会导致压力分布并不完全对称.

2)本程序利用辛朴森公式求压力对面积的积分得到承载能力也会一定误差.

3.1 前端角γ对结果影响

在求轴颈轴承静态特性时，将圆柱形气膜展开成如图1长方形计算区域过程中，从什么位置剪开，即气膜圆柱面剪开线与其垂直轴线组成的平面和按竖直轴线剪开平面的夹角γ(前端角)，对计算结果是否会产生影响进行探讨.本文按照 10°、22.5°和30°方向分别展开，并且都以两节流小孔的中心线将计算域划分为8个区域.图3为3种情况下的计算结果，从图中可以看出，载荷误差Wσ都在一个相对比较小的范围，并且没有递增或递减趋势，3种情况下得到的气体消耗量Q基本上都相等，因此说明前端角γ对求轴颈轴承静态性能并没有影响，即可以从任意位置将圆柱轴承展开成平面计算区域.

此规律为划分网格区域带来了便利，以8个供气孔为例，选择22.5°位置将气膜展开成长方形计算区域，可以得到大小相等的8个区域，使得在节流小孔位置布网格点更方便，另外还可以避免被剪开的两条边上碰到有节流小孔，使程序处理变得复杂.

图3 前端角γ对计算结果影响Fig.3 Load error and gas consumption versus pressure plot

3.2 网格数与收敛精度δ对结果影响

图4 网格数和收敛精度δ对计算结果影响Fig.4 Load error and gas consumption versus pressure plot

一般数值计算问题中，网格数与计算精度成正比，即网格数越多计算结果精度越高.但本文发现此规律在超松弛迭代法(SOR)中并不适用，网格数与计算结果精度并不是成正比，而是与收敛精度δ有关，当收敛精度8选择不当，网格数越多计算结果越不精确，有的甚至得到错误的结果，本文对网格数与收敛精度δ对计算结果的影响进行了研究.

图4表明，收敛精度取δ=10-4时，当长宽的网格数m=n=128和256时，虽然程序也满足了收敛条件，但载荷误差Wσ很大，特别是当m=n=256时，载荷误差δ能达到约20 N，很显然计算结果是不可信的.相比而言，网格数m=n=32和64时，载荷误差Wσ在一个比较小的范围，得到的结果更接近实际.当网格数m=n=128不变，提高收敛精度到δ=10-5，载荷误差Wδ在一个比较小的范围内变化，气体消耗量与网格数m=n=32和64时很接近，继续提高收敛精度到δ=5×10-6，载荷误差Wσ并没有继续减小，而也是在相同小的范围内变化.

因此，在用有限差分解Reynolds方程过程中，网格数需要与收敛精度δ相匹配，以本文为例，本程序采用具有二阶精度的中心差分格式离散Reynolds方程，数值解精度可以达到O(Δx2)，步长Δx由网格点数(m，n)和轴承特征长度(或宽度)决定，收敛精度δ值至少要小于Δx2数量级才能保证超松弛迭代法(SOR)也具有二阶精度，本文中如果是m=n=64，δ值取1 ×10-4合适，如果 m=n=128，δ值至少需取5×10-5合适，依此类推.网格数增多需要同时提高收敛精度δ，网格数一定，当收敛精度δ值小于某临界值以后，再小对提高数值解精度并没有太大帮助，并且还增加了计算成本，因此一般情况下δ值取到1×10-6即可.

图5 程序流程图Fig.5 Flow chart for calculation program

3.3 偏位角α影响

偏位角α由负载和转速唯一确定，当转子的转速不为零时，轴颈的偏位角也不为零［10］，此时的气膜厚度表达式应该按实际偏位角方向来计算，而大多数求解过程都是忽略了偏位角α对气膜厚度的影响.

本程序先假设一个小的偏位角度，在偏位角方向上将气膜展开，此时的气膜厚度表达式即为实际的气膜厚度，当满足收敛条件以后，求得实际的偏位角α'=arctan(Wt/Wn)，判断α'和α是否相等，不相等则改变初设的α重新返回计算，直至α'和α的差值与α的比小于ζ，此时求得的α'即为实际情况下的偏位角.

为了节约篇幅，本文以2个节流孔来说明其迭代的流程图如图5，节流孔数增多就是再添加如区域I，Π的质量守恒循环，流程图中虚线部分表示有转速时需要加入的判断偏位角是否合适的嵌套循环(如图中A位置)，在求轴颈轴承静态性能时按实线流程直接计算轴承流量.

4 算例

文献［11］利用Marc软件中的轴承分析模块求解了如下参数的气体轴承静态性能，并且和试验进行了对比，具体参数如下:轴承直径D=45 mm，长度L=35 mm，共有两排环面供气孔，每排8个，节流孔直径d=0.3 mm，节流孔至断面距离l=10 mm，轴承平均间隙Cr=0.02 mm，供气压力ps=6.5×105Pa.

气体消耗量和承载能力随偏心率变化的对比关系分别如图6、7.

与文献中计算结果相比，误差在5%以下，产生误差的原因除了程序本身以外，还与各种方法对文献中未提供的其它参数取值差异有关，另外收敛精度和网格数的选择对计算结果也有一定影响.

图6 气体消耗量计算结果比较Fig.6 Comparison of the results for gas consumption

图7 承载能力计算结果比较Fig.7 Comparison of the results for load capacity

5 结论

1)本文利用有限差分法编写了一套可以求静压轴颈轴承动静态性能的程序，引入了零偏心情况下的理想模型，讨论了前端角γ对求轴颈轴承静态性能计算结果的影响，结果证明前端角γ对求轴颈轴承静态性能计算结果没有影响;

2)对网格数与收敛精度δ进行研究发现，网格数需要与收敛精度δ相匹配，否则对计算结果会造成很大影响;

3)引入了偏位角α的嵌套循环，考虑了实际气膜厚度表达式对轴承性能的影响，计算结果与有实验验证的算例比较，两者误差在5%以下，该方法已成功用于对俄罗斯氦透平膨胀机静压气体轴承端的改造上，效果良好.

［1］吴华鹏，孔宪梅.GMRES算法在雷诺方程数值解法中的应用［J］.润滑与密封，2000(1):2-4.

WU Huapeng，KONG Xianmei.The application of GMRES algorithm in calculation Reynolds equation［J］.Lubrication Engineering，2000(1):2-4.

［2］彭万欢，蒲如平.气体静压轴颈轴承的静态特性分析及其试验研究［J］.润滑与密封，2006(12):204-206.

PENG Wanhuan，PU Ruping.A static performance analysis and experimental study of aerostatic journal bearing［J］.Lubrication Engineering，2006(12):204-206.

［3］LO Chengying，WANG Chengchi，LEE Yuhan.Performance analysis of high-speed spindle aerostatic bearings［J］.Tribology International，2005，38:5-14.

［4］张鸣，朱煜，段广洪.基于FEMLAB的气浮轴承静态性能求解方法［J］.润滑与密封，2006(2):63-65.

ZHANG Ming ，ZHU Yu ，DUAN Guanghong.A new method based FEMLAB to predictstatic performance ofair bearing［J］.Lubrication Engineering，2006(2):63-65.

［5］十合晋一.气体轴承设计、制作与应用［M］.韩涣臣，译.北京:宇航出版社，1988:12-20.

［6］POWELL J W.Design of aerostatic bearings［M］.London:Machinery Publishing Co.Ltd，1970:43-47.

［7］王云飞.气体润滑理论与气体轴承设计［M］.北京:机械工业出版社，1999:19-30.

［8］刘暾，刘育华，陈世杰.静压气体润滑［M］.哈尔滨:哈尔滨工业大学出版社，1990:40-60.

［9］黄云清.数值计算方法［M］.北京:科学出版社，2009:70-85.

［10］池长青.气体动静压轴承的动力学和热力学［M］.北京:北京航空航天大学，2008:60-74.

［11］张恩龙.高速全支撑气体静压电主轴的承载特性研究［D］.广州:广州工业大学，2005:65-86.