, 周凯 清华大学 机械工程系, 北京 100084
收稿日期：2017-05-24
基金项目：赛尔网络下一代互联网技术创新项目（NGII20150803）
作者简介：李雅哲(1990-), 男, 博士研究生。E-mail:yz-li13@tsinghua.edu.cn
摘要：空气静压轴承以空气作为润滑剂，尤其适用于高速应用场合。该文基于Reynolds方程，在极坐标条件下推导出止推空气静压轴承的有限差分模型。提出改进型流量误差反馈迭代法并应用于止推空气静压轴承的迭代计算中，采用伴随迭代随时更新的流量相对误差作为反馈变量，并随时修正速度因子。研究改进型流量误差迭代法的收敛特性，具体分析速度因子和迭代初值对收敛速度的影响。结果表明：迭代初值对收敛速度影响不明显，改进型方法可以增加速度因子裕度。通过研究轴承节流参数对轴承性能的影响表明：配合良好的节流孔和气隙可以显著提升轴承刚度。
关键词：流量误差反馈迭代法止推空气静压轴承Reynolds方程有限差分法(FDM)
Flow-difference feedback iteration method for aerostatic bearings
LI Yazhe
, ZHOU KaiDepartment of Mechanical Engineering, Tsinghua University, Beijing 100084, China
Abstract: Aerostatic bearings using air as the lubricant are especially useful in high-speed situations. This study of aerostatic bearings used a finite difference model of an aerostatic thrust bearing in polar coordinates based on the Reynolds equation with an improved flow-difference feedback iteration method. The relative error in the flow equilibrium updated by the iteration is used as the feedback variable to rectify the convergence rate factors. The initial value does not strongly influence the convergence rate, while the improved flow-difference feedback iteration method increases the effectiveness of the convergence rate factors. An analysis of the influence of the throttle parameter on the bearing performance shows that the proper orifices and clearance significantly improve the bearing stiffness.
Key words: flow-difference feedback iteration methodaerostatic thrust bearingReynolds equationfinite difference method (FDM)
与机械轴承相比，空气静压轴承具有低摩擦、低损耗和寿命长等优点；空气静压轴承以压缩空气为润滑介质，相比液体静压轴承，具有低发热和环境友好等优势。但是空气静压轴承刚度往往较低，因此空气静压轴承节流参数的优化设计十分重要。准确、高效的计算方法是空气静压轴承设计的重要工具，也是提升轴承承载力和刚度的关键技术。
在空气静压轴承的发展早期，静压轴承的计算和设计主要依赖于工程经验[1]。伴随数值计算的发展，以Reynolds方程为基础的数值解法，成为空气静压轴承的主流计算手段[2]，国内外均有****对其进行完善和改进，并将其应用于静压轴承设计中。Fransse等[3]利用计算流体力学(computational fluid dynamics，CFD)研究了止推空气静压轴承的动力学特性；Nicoletti等[4]通过修正Reynolds方程，并应用于多孔材料的空气静压轴承的静态和动态性能研究中；Zhu等[5]利用大涡模拟方法研究了空气静压轴承的流场瞬态特性；也有****利用ANSYS和FLUENT等商用软件进行计算，例如Neves等[6]进行了针对径向空气静压轴承的流量系数的研究。
有限差分法(finite difference method，FDM)和有限单元法(finite element method，FEM)是Reynolds方程的2种常用数值解法。Gero等[7]在稳态黏性可压缩条件下，对以上2种计算方法的效率和精度进行了比较，FDM在对比实验中具有较小的计算误差和较高的计算效率；Tala-Ighil等[8]利用FDM研究了空气静压轴承的热效应。应用有限差分法进行静压轴承数值计算，关键在于寻找适用于压力比例系数环的迭代方法。Lo等[9]对连续分割法和比例分割法[10]进行了比较；Chen等[11]采用过度松弛法(successive over-relaxation，SOR)进行迭代计算，并探究了几何参数对空气静压轴承的性能影响。流量误差反馈迭代法(flow-difference feedback iteration method，FDFIM)[12]相比以上迭代方法具有明显的收敛速度优势；适用于迭代量数量多、迭代量之间相互耦合影响的复杂情况。在文[12]中，详述了该方法应用于径向空气静压轴承的详细建模和求解过程，并利用该方法研究了径向轴承节流参数对轴承性能的影响。FDFIM中的速度因子是该方法的关键迭代参数，与节流参数不匹配的速度因子会造成迭代效率低下，或是迭代发散的情况。但是现有方法中，速度因子选择有赖于经验，并存在速度因子裕度小的问题。
本文提出改进型流量误差反馈迭代法，以获得更好的速度因子裕度，并将其应用于止推空气静压轴承计算中。重点研究改进型流量误差迭代法的收敛特性，具体分析速度因子、迭代初值对迭代速率的影响，并利用该方法研究止推空气静压轴承节流参数对静压轴承承载力和刚度的影响。
1 止推空气静压轴承流场模型数值计算的基础是建立数学模型，本文针对止推空气静压轴承进行控制方程的无量纲化和离散化，并给定边界条件。
1.1 控制方程及其离散化空气静压轴承工作时，需向节流孔供给高压空气。对于止推空气静压轴承(见图 1)，当轴承有偏心量e时，上下部分气隙(润滑流场厚度)有所差异，位于小气隙处的气压大于大气隙处的气压，压差的累加形成了空气轴承的承载力[13]。
 |
| 图 1 止推空气静压轴承结构 |
| 图选项 |
空气静压轴承的润滑流场分析是黏性可压缩流体问题[14]，假设气压在润滑流场厚度方向上为常数，其压力分布服从Reynolds方程[15]，止推空气静压轴承润滑流场的数值计算在极坐标系下进行。忽略动压效应(低速条件)，在稳态等温过程条件下，Reynolds方程表达式为
| $\begin{array}{l}\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r{h^3}\frac{{\partial \left( {{p^2}} \right)}}{{\partial r}}} \right) + \\\frac{1}{{{r^2}}}\frac{\partial }{{\partial \theta }}\left( {{h^3}\frac{{\partial \left( {{p^2}} \right)}}{{\partial \theta }}} \right) = 0, \end{array}$ | (1) |
| $\frac{\partial }{{\partial \bar r}}\left( {\bar r\frac{{\partial \bar f}}{{\partial \bar r}}} \right) + \frac{1}{{\bar r}}\frac{\partial }{{\partial \theta }}\left( {\frac{{\partial \bar f}}{{\partial \theta }}} \right) = 0.$ | (2) |
| $\left\{ \begin{array}{l}\bar p = \frac{p}{{{p_{\rm{a}}}}}, \bar r = \frac{r}{{{R_0}}};\\\bar f = {{\bar p}^2}, \varepsilon = e/{h_0}.\end{array} \right.$ |
 |
| 图 2 极坐标系及有限差分网格 |
| 图选项 |
有限差分网格(见图 2)采用极坐标下的结构化网格，并将节流孔中心布置于网格节点。利用中心差分格式[16](5点差分)，可得式(2) 的离散形式为
| $\begin{array}{l}\left[{\frac{{{{\bar r}_{i, j}}}}{{{{\left( {\mathit{\Delta }\bar r} \right)}^2}}} + \frac{1}{{2\Delta \bar r}}, \frac{{{{\bar r}_{i, j}}}}{{{{\left( {\Delta \bar r} \right)}^2}}}-\frac{1}{{2\Delta \bar r}}, -\frac{{2{{\bar r}_{i, j}}}}{{{{\left( {\Delta \bar r} \right)}^2}}}-\frac{2}{{{{\bar r}_{i, j}}{{\left( {\Delta \theta } \right)}^2}}}, } \right.\\\;\;\;\;\;\;\;\;\;\;\;\;\left. {\frac{1}{{{{\bar r}_{i, j}}{{\left( {\Delta \theta } \right)}^2}}}, \frac{1}{{{{\bar r}_{i, j}}{{\left( {\Delta \theta } \right)}^2}}}} \right]\left\{ \begin{array}{l}{{\bar f}_{i + 1, j}}\\{{\bar f}_{i -1, j}}\\{{\bar f}_{i, j}}\\{{\bar f}_{i, j -1}}\\{{\bar f}_{i, j + 1}}\end{array} \right\} = 0.\end{array}$ | (3) |
| ${\left[{{{\bar f}_{i + 1, j}}\;\;{{\bar f}_{i-1, j}}\;\;{{\bar f}_{i, j}}\;\;{{\bar f}_{i, M}}\;\;{{\bar f}_{i, j + 1}}} \right]^{\rm{T}}}.$ | (4) |
| ${\left[{{{\bar f}_{i + 1, j}}\;\;{{\bar f}_{i-1, j}}\;\;{{\bar f}_{i, j}}\;\;{{\bar f}_{i, j-1}}\;\;{{\bar f}_{i, 1}}} \right]^{\rm{T}}}.$ | (5) |
| ${{\bar f}_{{\rm{air}}}} = 1.$ | (6) |
| ${{\bar f}_{{\rm{orifice}}}} = \frac{{{\beta ^2}p_0^2}}{{p_a^2}}.$ | (7) |
2 改进型FDFIM由于β在离散模型中仍然未知，节流孔处的边界条件尚不明确。需引入流体连续性方程描述的流量平衡关系，使数学模型完备，本文提出改进型FDFIM对β迭代。
假设空气流入节流孔是绝热过程，流入节流孔的质量流量[17]为
| ${{\dot m}_{{\rm{in}}}} = A\zeta \sqrt {2{\rho _0}{p_0}} \cdot \psi .$ | (8) |
考虑润滑层流体具有Newton黏性特征，则在上下流体边界r方向流速为0，利用简化的Navier-Stokes方程[18]
| $\left\{ \begin{array}{l}\frac{{\partial p}}{{\partial r}} = \mu \frac{{{\partial ^2}u}}{{\partial {z^2}}}, \\\frac{{\partial p}}{{\partial z}} = 0, \\\frac{{\partial p}}{{\partial \theta }} = r\mu \frac{{{\partial ^2}v}}{{\partial {z^2}}}, \end{array} \right.$ | (9) |
| $\begin{array}{l}\;\;\;\;{{\dot m}_r} = \int_0^h {\rho {\mu _r}\left( {{r_i}\Delta \theta } \right){\rm{d}}z = } \\\int_0^h {\rho \frac{{-1}}{{2\mu }}\frac{{\partial p}}{{\partial r}}z\left( {h-z} \right)\left( {{r_i}\Delta \theta } \right){\rm{d}}z} = \\\;\;\;\;\;\;\;\;\;\;\frac{{-\rho {h^3}{r_i}\Delta \theta }}{{12\mu }}\frac{{\partial p}}{{\partial r}}.\end{array}$ | (10) |
| $\frac{p}{\rho } = \frac{{{p_{\rm{a}}}}}{{{\rho _{\rm{a}}}}}.$ | (11) |
 |
| 图 3 r方向质量流量 |
| 图选项 |
| $\begin{array}{l}{{\dot m}_{\rm{r}}} = \frac{{{p_{\rm{P}}}{\rho _{\rm{a}}}{h^3}\left( {{p_{{\rm{in}}}}- {p_{{\rm{out}}}}} \right)\left[{N{R_2} + \left( {{R_1}-{R_2}} \right)\left( {i-1} \right)} \right]\mathit{\Delta }\theta }}{{12{p_{\rm{a}}}\mu N\mathit{\Delta }r}}, \\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\{ \begin{array}{l}{p_{\rm{P}}} = \frac{{{p_{i, j}} + {p_{i, j + 1}}}}{2}, \\{p_{{\rm{in}}}} = \frac{{{p_{i -1, j}} + {p_{i -1, j + 1}}}}{2}, \\{p_{{\rm{out}}}} = \frac{{{p_{i + 1, j}} + {p_{i + 1, j + 1}}}}{2}.\end{array} \right.\end{array}$ | (12) |
由流体连续性方程，在空气静压轴承任意一节流孔周围区域(如图 2阴影所示)，流体质量流量应守恒。定义流量平衡相对误差(简称为相对误差)为
| $E = \frac{{{{\dot m}_{{\rm{in}}}} + \int {\dot m_{\rm{r}}^{\left( {i = 1} \right)}{\rm{d}}\theta-\dot m_{\rm{r}}^{\left( {i = N} \right)}{\rm{d}}\theta } }}{{{{\dot m}_{{\rm{in}}}} + \int {\dot m_r^{\left( {i = 1} \right)}{\rm{d}}\theta } }}.$ | (13) |
改进型FDFIM为
| ${\beta _{\left( k \right)}} = \delta \cdot \sum {E_{\left( k \right)}^2.} $ | (14) |
 |
| 图 4 改进型FDFIM |
| 图选项 |
3 改进型FDFIM的收敛特性收敛速度和稳定性取决于计算参数的选择。迭代初值是常见的数值计算参数，而δ是改进型FDFIM的特有参数。
用于计算实验的止推空气静压轴承的几何参数如下：R1为75 mm，R2为50 mm，R0为62.5 mm，单排节流孔布置12个节流孔，节流孔直径d为0.2 mm，h0为20 μm，p0为0.5 MPa，转速ω为0。收敛要求E≤10-3。
设定迭代初值β(0)，选取120×10规格的网格，研究δ对收敛速度的影响。偏心率ε为0.1，δ为0.6、0.9、1.2、1.5时，分别历经419、236、207、167次迭代计算收敛，图 5是相应的(前40次)收敛曲线。而当δ>1.5时，收敛曲线发散。以上计算结果表明：较大的δ有较高的计算效率(即较少的迭代次数)。
 |
| 图 5 不同速度因子下迭代方法的收敛特性 |
| 图选项 |
设定δ为0.8，选取120×10规格的网格，研究迭代初值对收敛速度的影响。ε为0.1，迭代初值β(0)为0.4、0.5、0.6、0.7时，分别历经315、309、300、309次迭代计算收敛，图 6是相应的收敛曲线。越接近终点值的初值选择，所需的迭代次数越少。应用改进型FDFIM，初值条件不同导致的收敛速度差异不明显。同时，在不同的初值条件下，流量误差反馈法都表现出较好的收敛稳定性。
 |
| 图 6 不同迭代初值下迭代方法的收敛特性 |
| 图选项 |
设定迭代初值β(0)为0.5，选取120×10规格的网格，将改进型行FDFIM的δ裕度(保证收敛条件下)与原方法进行对比，结果见表 1。对于2种方法而言，当δ小于表中所示数值时，均能保证收敛，而表中数值越大则代表有越大的δ选择范围，该方法的δ裕度越大。对比实验结果表明：相较原方法改进型FDFIM在不同偏心率下均有更大的速度因子裕度。
表 1 流量误差反馈迭代法的δ裕度
| ε | δ裕度 | |
| FDFIM | 改进型FDFIM | |
| 0.1 | 0.4 | 1.5 |
| 0.2 | 0.4 | 1.0 |
| 0.3 | 0.3 | 0.8 |
| 0.4 | 0.3 | 0.7 |
表选项
4 止推空气静压轴承结构参数研究应用改进型FDFIM研究节流参数对止推空气静压轴承的承载力W和刚度Kw的影响。选取止推空气静压轴承几何参数：R1为75 mm，R2为50 mm，R0为62.5 mm，单排节流孔布置12个节流孔；选取64×10网格，收敛要求E≤10-6，建议δ为0.5。
设置p0为0.5 MPa，考察d为0.2、0.3、0.4 mm时，承载力和刚度随h0的变化规律。图 7a是止推轴承在ε为0.1处的承载力随气隙的变化规律，直径越小的节流孔在相对较小的平均间隙处取得越大承载力，而不同节流孔径下能取得的最大承载力大小相似。图 7b是径向轴承在ε为0处的刚度变化规律，越小节流孔对应的最大刚度处气隙越小；而越小节流孔下的最大刚度越大。
 |
| 图 7 节流孔直径对静压气体轴承承载力和刚度的影响 |
| 图选项 |
设置d为0.2 mm，考察p0为0.4、0.5、0.6 MPa时，承载力和刚度随轴承h0的变化规律。图 8a是止推轴承在ε为0.1处的承载力随p0的变化规律。承载力随着供气压力的增加而增大，并且不同压力下承载力极值点的平均气隙相近。图 8b是止推轴承在偏心率ε为0处的刚度变化规律，刚度随着供气压力的增大而增大，并且不同压力下刚度极值点的平均气隙相近。
 |
| 图 8 供气压力对静压气体轴承承载力和刚度的影响 |
| 图选项 |
设置p0为0.5 MPa和d为0.2 mm，考察止推空气静压气体轴承在不同ε下的承载力和刚度随h0的变化规律。图 9a是径向轴承承载力在ε分别为0、0.1、0.2、0.3时承载力的变化情况，承载力随偏心率增加而增大，不同偏心状态下对应的承载力极值点的平均气隙略有不同。图 9b是径向轴承承载力在ε分别为0、0.1、0.2、0.3时刚度的变化情况，ε为0时能取得的最大刚度值最大，不同偏心状态下对应的最大刚度处的平均间隙略有不同。
 |
| 图 9 静压轴承在不同偏心率下的承载力和刚度 |
| 图选项 |
5 结论本文提出一种改进型流量误差反馈迭代法，并将其应用于止推空气静压轴承流场的有限差分解法中。研究了改进型流量误差迭代法的收敛特性，结果表明：迭代初值对收敛速度影响不大，该迭代方法在不同初值条件均表现出较好的收敛行为；而选定适合的速度因子可以大幅提升收敛速度，改进型流量误差反馈迭代法改善了速度因子裕度。通过研究节流参数对轴承性能的影响表明：配合良好的节流孔径和气隙可以显著提高轴承刚度，而较小的节流孔可能获得最大的刚度；而在较大的供气压力和较小的偏心率下轴承可获得更大的承载刚度。
参考文献
| [1] | Powell J W. Design of Aerostatic Bearings[M]. London: Machinery Publishing, 1970. |
| [2] | Szeri A Z. Fluid Film Lubrication:Theory and Design[M]. Cambridge: Cambridge University Press, 2005. |
| [3] | Franssen R M, Potze W, de Jong P, et al. Large amplitude dynamic behavior of thrust air bearings:Modeling and experiments[J]. Tribology International, 2017, 109: 460–466. DOI:10.1016/j.triboint.2016.12.024 |
| [4] | Nicoletti R, Purquerio B D M, Silveira Z D C. The effect of permeability distribution on the numerical analysis of aerostatic ceramic porous bearings[J]. Lubrication Science, 2013, 25(2): 185–194. DOI:10.1002/ls.v25.2 |
| [5] | Zhu J, Chen H, Chen X. Large eddy simulation of vortex shedding and pressure fluctuation in aerostatic bearings[J]. Journal of Fluids & Structures, 2013, 40(7): 42–51. |
| [6] | Neves M T, Schwarz V A, Menon G J. Discharge coefficient influence on the performance of aerostatic journal bearings[J]. Tribology International, 2010, 43(4): 746–751. DOI:10.1016/j.triboint.2009.11.001 |
| [7] | Gero L R, Ettles C M M. An evaluation of finite difference and finite element methods for the solution of the Reynolds equation[J]. Tribology Transactions, 1986, 29(2): 166–172. |
| [8] | Tala-Ighil N, Fillon M. A numerical investigation of both thermal and texturing surface effects on the journal bearings static characteristics[J]. Tribology International, 2015, 90: 228–239. DOI:10.1016/j.triboint.2015.02.032 |
| [9] | Lo C Y, Wang C C, Lee Y H. Performance analysis of high-speed spindle aerostatic bearings[J]. Tribology International, 2005, 38(1): 5–14. DOI:10.1016/j.triboint.2004.04.008 |
| [10] | 刘暾, 刘育华, 陈世杰. 静压气体润滑[M]. 哈尔滨: 哈尔滨工业大学出版社, 1990.LIU Dun, LIU Yuhua, CHEN Shijie. Aerostatic Lubrication[M]. Harbin: Harbin Institute of Technology Press, 1990. (in Chinese) |
| [11] | Chen Y S, Chiu C C, Cheng Y D. Influences of operational conditions and geometric parameters on the stiffness of aerostatic journal bearings[J]. Precision Engineering, 2010, 34(4): 722–734. DOI:10.1016/j.precisioneng.2010.04.001 |
| [12] | Li Y, Zhou K, Zhang Z. A flow-difference feedback iteration method and its application to high-speed aerostatic journal bearings[J]. Tribology International, 2014, 84: 132–141. |
| [13] | 党根茂. 气体润滑技术[M]. 南京: 东南大学出版社, 1990.DANG Genmao. Technology of Gas Lubrication[M]. Nanjing: Southeast University Press, 1990. (in Chinese) |
| [14] | Hori Y. Hydrodynamic Lubrication[M]. Berlin: Springer Science & Business Media, 2006. |
| [15] | 王云飞. 气体润滑理论与气体轴承设计[M]. 北京: 中国机械工业出版社, 1999.WANG Yunfei. Gas Lubricated Theory and Design Manual of Gas Bearings[M]. Beijing: China Machinery Industry Press, 1999. (in Chinese) |
| [16] | Stoer J, Bulirsch R. Introduction to Numerical Analysis[M]. Berlin: Springer Science & Business Media, 2013. |
| [17] | 张君安. 高刚度空气静压轴承研究[D]. 西安: 西北工业大学, 2006. ZHANG Junan. Research on High Stiffness Aerostatic Bearing[D]. Xi'an:Northwestern Polytechnic University, 2006. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10699-2007214186.htm |
| [18] | Zheng S F, Jiang S Y. Improved finite difference method for pressure distribution of aerostatic bearing[J]. Journal of Southeast University, 2009, 25(4): 501–505. |
