闫佳佳1,刘占生1,张广辉1,李豫1,徐亮2
(1. 哈尔滨工业大学 能源科学与工程学院 哈尔滨 150001; 2. 中航工业金城南京机电液压工程研究中心 南京 211102)
摘要:
为提高箔片轴承性能，提出一种新型的平箔片厚度轴向变化的波箔型箔片轴承.该新型箔片轴承的平箔片在轴承中间区域有一周向贯通的宽槽，凹槽深度在微米量级.考虑平箔片厚度的轴向变化，采用二维厚板有限元模型计算平箔片刚度.基于有限差分法和Newton-Raphson迭代法耦合求解压力控制方程及箔片变形方程，获得新型箔片轴承的气膜压力分布及承载力；基于小扰动法，计算轴承动力学特性系数.仿真研究平箔片宽槽深度对新型箔片轴承性能的影响.结果表明：与传统箔片轴承相比，平箔片厚度轴向变化的新型箔片轴承可有效减少润滑气体端泄，增大流体动压效果，提高轴承承载力并改善动力学特性；随着平箔片宽槽深度的增大，新型箔片轴承的承载力及主刚度进一步增大；平箔片宽槽深度为8 μm的新型箔片轴承，其承载力和主刚度分别增大约11.8 %、12.1 %.
关键词: 波箔型箔片轴承 平箔片 承载力 动力学特性 名义平均半径间隙 平箔片宽槽深度
DOI：10.11918/j.issn.0367-6234.201703024
分类号:V229；TH133.37
文献标识码:A
基金项目:国家自然科学基金(11572098)；航空科学基金(2015287008)
Performance of a novel foil bearing with top foil thickness variation in axial direction
YAN Jiajia1,LIU Zhansheng1,ZHANG Guanghui1,LI Yu1,XU Liang2
(1. School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, China; 2. Nanjing Engineering Institute of Aircraft Systems, Jincheng Aviation Industry Corporation of China, Nanjing 211102, China)
Abstract:
To improve the foil bearing performance, a novel foil journal bearing with top foil thickness variation in axial direction was proposed. A wide groove extending in circumferential direction, within micron magnitude, was manufactured in the mid-plane zone of top foil. To consider the top foil thickness variation, the 2D thick plate finite element model was employed to simulate the top foil deformation. Based on Newton-Raphson iterative method and finite difference method, the coupled pressure equation and foil deflection equation were solved, and the gas pressure distribution and load capacity of the novel foil bearing were obtained. The dynamic characteristics were solved using perturbation method and the influence of the wide groove depth on bearing performance was investigated theoretically. The results indicate that the novel foil bearing with top foil thickness variation axially can decrease the gas edge leakage and improve hydrodynamic pressure effectively, so the load capacity and dynamic characteristics are improved. Moreover, with the increment of the wide groove depth on top foil, the load capacity and direct stiffness are enhanced further. For the foil journal bearing with a wide groove depth of 8 μm, the load capacity and direct stiffness increase by about 11.8% and 12.1% respectively, compared with the traditional foil journal bearing.
Key words: bump-type foil bearing top foil load capacity dynamic characteristics nominal average radial clearance wide groove depth on top foil
闫佳佳, 刘占生, 张广辉, 李豫, 徐亮. 平箔片厚度轴向变化的新型箔片轴承性能[J]. 哈尔滨工业大学学报, 2018, 50(1): 59-67. DOI: 10.11918/j.issn.0367-6234.201703024.
YAN Jiajia, LIU Zhansheng, ZHANG Guanghui, LI Yu, XU Liang. Performance of a novel foil bearing with top foil thickness variation in axial direction[J]. Journal of Harbin Institute of Technology, 2018, 50(1): 59-67. DOI: 10.11918/j.issn.0367-6234.201703024.
基金项目 国家自然科学基金(11572098);航空科学基金(2015287008) 作者简介 闫佳佳(1987—)，男，博士研究生;
刘占生(1962—)，男，教授，博士生导师 通信作者 刘占生，lzs@hit.edu.cn 文章历史 收稿日期: 2017-03-06
Contents -->Abstract Full text Figures/Tables PDF
平箔片厚度轴向变化的新型箔片轴承性能
闫佳佳1, 刘占生1
, 张广辉1, 李豫1, 徐亮2 1. 哈尔滨工业大学 能源科学与工程学院 哈尔滨 150001;
2. 中航工业金城南京机电液压工程研究中心 南京 211102
收稿日期: 2017-03-06
基金项目: 国家自然科学基金(11572098);航空科学基金(2015287008)
作者简介: 闫佳佳(1987—)，男，博士研究生;
刘占生(1962—)，男，教授，博士生导师
通信作者: 刘占生，lzs@hit.edu.cn
摘要: 为提高箔片轴承性能，提出一种新型的平箔片厚度轴向变化的波箔型箔片轴承.该新型箔片轴承的平箔片在轴承中间区域有一周向贯通的宽槽，凹槽深度在微米量级.考虑平箔片厚度的轴向变化，采用二维厚板有限元模型计算平箔片刚度.基于有限差分法和Newton-Raphson迭代法耦合求解压力控制方程及箔片变形方程，获得新型箔片轴承的气膜压力分布及承载力；基于小扰动法，计算轴承动力学特性系数.仿真研究平箔片宽槽深度对新型箔片轴承性能的影响.结果表明：与传统箔片轴承相比，平箔片厚度轴向变化的新型箔片轴承可有效减少润滑气体端泄，增大流体动压效果，提高轴承承载力并改善动力学特性；随着平箔片宽槽深度的增大，新型箔片轴承的承载力及主刚度进一步增大；平箔片宽槽深度为8 μm的新型箔片轴承，其承载力和主刚度分别增大约11.8 %、12.1 %.
关键词: 波箔型箔片轴承 平箔片 承载力 动力学特性 名义平均半径间隙 平箔片宽槽深度
Performance of a novel foil bearing with top foil thickness variation in axial direction
YAN Jiajia1, LIU Zhansheng1
, ZHANG Guanghui1, LI Yu1, XU Liang2 1. School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, China;
2. Nanjing Engineering Institute of Aircraft Systems, Jincheng Aviation Industry Corporation of China, Nanjing 211102, China
Abstract: To improve the foil bearing performance, a novel foil journal bearing with top foil thickness variation in axial direction was proposed. A wide groove extending in circumferential direction, within micron magnitude, was manufactured in the mid-plane zone of top foil. To consider the top foil thickness variation, the 2D thick plate finite element model was employed to simulate the top foil deformation. Based on Newton-Raphson iterative method and finite difference method, the coupled pressure equation and foil deflection equation were solved, and the gas pressure distribution and load capacity of the novel foil bearing were obtained. The dynamic characteristics were solved using perturbation method and the influence of the wide groove depth on bearing performance was investigated theoretically. The results indicate that the novel foil bearing with top foil thickness variation axially can decrease the gas edge leakage and improve hydrodynamic pressure effectively, so the load capacity and dynamic characteristics are improved. Moreover, with the increment of the wide groove depth on top foil, the load capacity and direct stiffness are enhanced further. For the foil journal bearing with a wide groove depth of 8 μm, the load capacity and direct stiffness increase by about 11.8% and 12.1% respectively, compared with the traditional foil journal bearing.
Key words: bump-type foil bearing top foil load capacity dynamic characteristics nominal average radial clearance wide groove depth on top foil
弹性箔片动压气体轴承是以柔性表面为支承的自作用式动压挠性轴承.与传统油润滑轴承相比，气体箔片轴承无需供油系统，具有可靠性高、环境与温度适应性强、维修周期长、摩擦损失小、抗冲击、自适应好等诸多优点[1]. 20世纪70年代以来，随着航空航天、原子能、低温工程等技术的发展以及对箔片轴承研究的不断深入，气体箔片轴承逐渐应用到了飞机空气循环机、低温透平膨胀机以及压缩机等高速轻载旋转机械中.与润滑油相比，润滑气体的黏度低，导致箔片轴承的承载力和动力学特性系数偏小[2].自箔片轴承[3]的概念诞生以来，为提高其承载力和动力学稳定性，各国学者提出了多种结构形式的气体箔片轴承.其中，多叶型箔片轴承和波箔型箔片轴承受到了广泛关注.
20世纪60年代，美国的Garrett AiResearch公司(目前为Honeywell)最早对多叶型箔片轴承进行了研究[1].大量的理论及试验研究证明，多叶型箔片轴承具有充足的承载力和阻尼特性; 此外，弹簧支承型多叶箔片轴承的提出进一步增大了其承载力[4-5].
与多叶型箔片轴承相比，波箔型箔片轴承的承载力相当，但阻尼略小.它由平箔片和波箔片组成，波箔片产生弹性变形，平箔片提供光滑表面.为提高轴承性能，国内外学者通过调整波箔片形式，提出了多种复杂结构的波箔型箔片轴承，包括双层波箔箔片轴承[6]、波箔沿周向分隔的箔片轴承[7]、可变节距的3瓦片箔片轴承[8]、周向和轴向刚度均变化的波箔箔片轴承[9]等.从提高轴承预载的角度出发，又发明了多楔型波箔箔片轴承[10]、加薄垫片的波箔箔片轴承[11]等.通过调整楔的型线或垫片厚度，轴承性能得到了提高.
在多叶型和波箔型箔片轴承的基础上，通过对柔性支承结构进行改进，提出了多种其他结构形式的箔片轴承.以压缩弹簧代替波箔片的箔片轴承，可有效降低转子过临界时的振动峰值[12].金属橡胶箔片轴承的提出，改善了箔片轴承的阻尼特性和稳定性[13-14].西安交通大学的赖天伟等[15]提出的多层鼓泡弹性支承箔片轴承在抑制转子不稳定涡动方面取得了显著的效果.湖南大学的冯凯等[16]发明的混合金属橡胶-波箔型箔片轴承，集成了金属橡胶箔片轴承和波箔型箔片轴承的优点.
针对不同结构形式的箔片轴承，根据轴承柔性支承结构的刚度是否沿轴向、周向、径向的单一或多个方向变化，Dellacorte和Valco将箔片轴承划分为3代[17].随着箔片轴承由第1代向第3代发展，其性能有显著的提高.分析其原因，周向或径向刚度变化的柔性支承结构可优化箔片轴承的收敛楔，而轴向刚度的变化可减少润滑气体的端泄，有助于增大气体动压效果，从而提高轴承承载力及动力学特性等性能.
为提高箔片轴承性能，本文从减少润滑气体端泄的角度出发，提出了一种平箔片厚度轴向变化的波箔型动压气体箔片轴承，并仿真研究了轴承的静、动态性能.其中，通过有限差分法和Newton-Raphson迭代法耦合求解获得轴承的气膜压力分布、承载力，采用小扰动法获得动力学特性系数.为考虑平箔片厚度的轴向变化，平箔片刚度采用考虑剪切变形的二维厚板有限元模型计算.系统分析了平箔片厚度轴向变化剧烈程度对新型箔片轴承静、动态性能的影响规律.
1 平箔片厚度轴向变化的波箔型箔片轴承结构图 1为平箔片厚度轴向变化的波箔型箔片轴承结构示意，主要由波箔片、平箔片和轴承座组成.波箔片充当类似弹簧的作用，在气体动压作用下产生弹性变形.轴承座上开有轴向狭槽用来插入平箔片.平箔片和波箔片均采用焊接的方式固定到轴承座上.
Figure 1
图 1 平箔片厚度轴向变化的波箔型箔片轴承结构 Figure 1 Schematic of the gas foil bearing with top foil thickness variation in axial direction 该新型箔片轴承与传统箔片轴承的显著区别是平箔片厚度沿轴向有所变化.如图 1所示，平箔片中间区域有一宽槽，位于转子一侧，宽槽深度在微米量级.得益于表面微加工技术的发展，对该平箔片表面宽槽的加工可采用激光表面加工技术、化学侵蚀、机械磨削等实现.该新型箔片轴承保留了第一代波箔型箔片轴承的优点.同时，平箔片表面的宽槽使得中间区域气膜厚度增大，两端气膜厚度减小，减少了润滑气体端泄，有效提高气体动压，改善了轴承性能.
2 理论分析 2.1 平箔片二维厚板模型本文采用文献[18]的二维厚板模型来分析平箔片变形，其单元结构及自由度如图 2所示，图中u、υ、w为单元内任一点沿x、y、z方向的位移，θx、θy为绕x、y方向的转角.
Figure 2
图 2 二维厚板单元的结构 Figure 2 Geometry of a rectangular element for 2D thick plate 该模型考虑了平箔片的弯曲变形及剪切变形，其剪切应变γxz、γyz与位移的关系如下：
${\gamma _{xz}} = \frac{{\partial u}}{{\partial z}} + \frac{{\partial w}}{{\partial x}},\;\;\;{\gamma _{yz}} = \frac{{\partial v}}{{\partial z}} + \frac{{\partial w}}{{\partial y}}.$
其中：
$u\left( {x,y,z} \right) = z{\theta _y}\left( {x,y} \right),$
$v\left( {x,y,z} \right) = - z{\theta _x}\left( {x,y} \right).$
厚板单元的总应变能表示为弯曲变形能和剪切变形能的总和：
$\begin{array}{l}{U_{\rm{e}}} = \frac{1}{2}\int_A {\frac{{h_{\rm{p}}^3}}{{12}}{{\left\{ \mathit{\boldsymbol{\chi }} \right\}}^{\rm{T}}}\mathit{\boldsymbol{D}}\left\{ \mathit{\boldsymbol{\chi }} \right\}{\rm{d}}A} + \\\;\;\;\;\;\;\;\;\frac{1}{2}\int_A {\mathit{\boldsymbol{\kappa }}{h_{\rm{p}}}{{\left\{ \mathit{\boldsymbol{\gamma }} \right\}}^{\rm{T}}}{\mathit{\boldsymbol{D}}^{\rm{s}}}\left\{ \mathit{\boldsymbol{\gamma }} \right\}{\rm{d}}A} .\end{array}$ (1)
式中：hp为单元厚度，A为单元面积，κ为剪切系数.
$\left\{ \begin{array}{l}\mathit{\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}}{\frac{E}{{{{\left( {1 - v} \right)}^2}}}}&{\frac{{Ev}}{{{{\left( {1 - v} \right)}^2}}}}&0\\{\frac{{Ev}}{{{{\left( {1 - v} \right)}^2}}}}&{\frac{E}{{{{\left( {1 - v} \right)}^2}}}}&0\\0&0&{\frac{E}{{2\left( {1 - v} \right)}}}\end{array}} \right],\\{\mathit{\boldsymbol{D}}^{\rm{s}}} = \frac{E}{{2\left( {1 + v} \right)}}\left[ {\begin{array}{*{20}{c}}1&0\\0&1\end{array}} \right],\\\left\{ \mathit{\boldsymbol{\chi }} \right\} = \left[ {\begin{array}{*{20}{c}}{ - \partial {\theta _y}/\partial x}\\{\partial {\theta _x}/\partial y}\\{\partial {\theta _x}/\partial x - \partial {\theta _y}/\partial y}\end{array}} \right],\left\{ \mathit{\boldsymbol{\gamma }} \right\} = \left[ {\begin{array}{*{20}{c}}{{\theta _y} + \partial w/\partial x}\\{ - {\theta _x} + \partial w/\partial y}\end{array}} \right].\end{array} \right.$ (2)
式中：E、v分别为平箔片材料的弹性模量、泊松比.
单元内任一点的位移可由节点位移和形函数计算得到：
$\left[ {\begin{array}{*{20}{c}}w\\{{\theta _x}}\\{{\theta _y}}\end{array}} \right] = \mathit{\boldsymbol{N}}{\left\{ \mathit{\boldsymbol{w}} \right\}_{\rm{e}}}.$ (3)
式中：{w}e为节点位移向量，N为形函数矩阵，可表示为
$\left\{ \mathit{\boldsymbol{w}} \right\}_{\rm{e}}^{\rm{T}} = \left\{ {{w_1}\;{\theta _{x1}}\;{\theta _{y1}} \cdots {w_4}\;{\theta _{x4}}\;{\theta _{y4}}} \right\},$
$\mathit{\boldsymbol{N}} = \left[ {\begin{array}{*{20}{c}}{{N_1}}&0&0& \cdots &{{N_4}}&0&0\\0&{{N_1}}&0& \cdots &0&{{N_4}}&0\\0&0&{{N_1}}& \cdots &0&0&{{N_4}}\end{array}} \right].$
将节点j的量纲一的坐标记为(ξj, ηj)，则形函数Nj可表示为
${N_j} = \frac{1}{4}\left( {1 + {\xi _j}\xi } \right)\left( {1 + {\eta _j}\eta } \right).$
其中：
$\xi = x/a,\;\;\;\;\eta = y/b.$
将方程(3)代入方程(1)和(2)中，得
${U_{\rm{e}}} = \frac{1}{2}\left\{ \mathit{\boldsymbol{w}} \right\}_{\rm{e}}^{\rm{T}}{\mathit{\boldsymbol{k}}_{\rm{e}}}{\left\{ \mathit{\boldsymbol{w}} \right\}_{\rm{e}}}.$
式中：ke为平箔片厚板单元的刚度矩阵，可写为
${\mathit{\boldsymbol{k}}_{\rm{e}}} = {\mathit{\boldsymbol{k}}^{\rm{f}}} + {\mathit{\boldsymbol{k}}^{\rm{s}}}.$
其中，kf和ks分别对应弯曲刚度矩阵和剪切刚度矩阵.
当平箔片材料各向同性时，kf和ks表示为
${\mathit{\boldsymbol{k}}^{\rm{f}}} = \frac{{Eh_{\rm{p}}^3}}{{48ab{{\left( {1 - v} \right)}^2}}}\left[ {\begin{array}{*{20}{c}}{\mathit{\boldsymbol{k}}_{11}^{\rm{f}}}&{}&{}&{对称}\\{\mathit{\boldsymbol{k}}_{21}^{\rm{f}}}&{\mathit{\boldsymbol{k}}_{22}^{\rm{f}}}&{}&{}\\{\mathit{\boldsymbol{k}}_{31}^{\rm{f}}}&{\mathit{\boldsymbol{k}}_{32}^{\rm{f}}}&{\mathit{\boldsymbol{k}}_{33}^{\rm{f}}}&{}\\{\mathit{\boldsymbol{k}}_{41}^{\rm{f}}}&{\mathit{\boldsymbol{k}}_{42}^{\rm{f}}}&{\mathit{\boldsymbol{k}}_{43}^{\rm{f}}}&{\mathit{\boldsymbol{k}}_{44}^{\rm{f}}}\end{array}} \right],$
${\mathit{\boldsymbol{k}}^{\rm{s}}} = \frac{{Eh_{\rm{p}}^3}}{{48ab{\beta _{\rm{s}}}}}\left[ {\begin{array}{*{20}{c}}{\mathit{\boldsymbol{k}}_{11}^{\rm{s}}}&{}&{}&{对称}\\{\mathit{\boldsymbol{k}}_{21}^{\rm{s}}}&{\mathit{\boldsymbol{k}}_{22}^{\rm{s}}}&{}&{}\\{\mathit{\boldsymbol{k}}_{31}^{\rm{s}}}&{\mathit{\boldsymbol{k}}_{32}^{\rm{s}}}&{\mathit{\boldsymbol{k}}_{33}^{\rm{s}}}&{}\\{\mathit{\boldsymbol{k}}_{41}^{\rm{s}}}&{\mathit{\boldsymbol{k}}_{42}^{\rm{s}}}&{\mathit{\boldsymbol{k}}_{43}^{\rm{s}}}&{\mathit{\boldsymbol{k}}_{44}^{\rm{s}}}\end{array}} \right].$
式中：βs=Ehp2/12κGb2为剪切因子，G为剪切模量.矩阵中各元素的表达式见文献[18].
上述刚度矩阵中计入了平箔片单元厚度hp的影响.考虑平箔片中间区域宽槽对单元厚度的影响，将单元刚度矩阵集成，得到平箔片的总体刚度矩阵KT.
2.2 箔片结构变形方程波箔片的刚度矩阵KB采用Iordanoff公式[19]计算，针对一端固定一端自由和两端自由的波箔给出了不同的刚度计算表达式.综合波箔片刚度矩阵和平箔片刚度矩阵，最终得到箔片结构整体刚度矩阵KG为
${\mathit{\boldsymbol{K}}_{\rm{G}}} = {\mathit{\boldsymbol{K}}_{\rm{T}}} + {\mathit{\boldsymbol{K}}_{\rm{B}}}.$
图 3给出了箔片结构相互作用关系，分析时认为轴承座是完全刚性的.波箔片线性弹簧的一端固定到轴承座内表面上，另一端连接到平箔片对应网格节点.平箔片一端固定到轴承座内表面.对于该固定端，有位移边界条件w=θx=θy=0.平箔片其余节点均可自由横向变形及转动.
Figure 3
图 3 箔片结构及轴承座间的相互作用 Figure 3 Interaction between foil structures and bearing housing 结合上述边界条件，箔片结构的位移可由式(4)计算得到：
${\mathit{\boldsymbol{K}}_{\rm{G}}}{\left\{ \mathit{\boldsymbol{U}} \right\}_{\rm{G}}} = {\left\{ \mathit{\boldsymbol{F}} \right\}_{\rm{G}}}.$ (4)
式中：{U}G为箔片结构节点位移矢量，{F}G为作用在平箔片表面的气膜压力矢量，表达式分别为
$\left\{ \mathit{\boldsymbol{U}} \right\}_{\rm{G}}^{\rm{T}} = \left\{ {\begin{array}{*{20}{c}}{{w_1}}&{{\theta _{x1}}}&{{\theta _{y1}}}&{{w_2}}&{{\theta _{x2}}}&{{\theta _{y2}}}& \cdots \end{array}} \right\},$
$\left\{ \mathit{\boldsymbol{F}} \right\}_{\rm{G}}^{\rm{T}} = \left\{ {\begin{array}{*{20}{c}}{{p_1}}&0&0&{{p_2}}&0& \cdots \end{array}} \right\}.$
其中：p1, p2, …为各节点处气膜压力.
提取所有节点的横向位移w，得到箔片结构的径向变形为
${\left\{ {{\mathit{\boldsymbol{U}}_R}} \right\}^{\rm{T}}} = \left\{ {\begin{array}{*{20}{c}}{{w_1}}&{{w_2}}& \cdots \end{array}} \right\}.$
2.3 气膜厚度表达式图 4给出了该新型箔片轴承的纵剖图.从图 4中可以看出轴承气膜厚度的分布规律.
Figure 4
图 4 新型箔片轴承纵剖图 Figure 4 Longitudinal section of the novel gas foil bearing 定义轴承名义平均半径间隙为C，轴承中间区域平箔片宽槽深度为Td，则新型箔片轴承两端区域和中间区域的间隙分别为
${C_1} = C - {T_{\rm{d}}}/2,$
${C_2} = C - {T_{\rm{d}}}/2.$
考虑箔片结构径向变形UR及转子的偏心，得到如下气膜厚度表达式：
$h = \left\{ \begin{array}{l}{C_1} = e\cos \left( {\theta - \varphi } \right) + {U_R},\;\;\;两端区域;\\{C_2} = e\cos \left( {\theta - \varphi } \right) + {U_R},\;\;\;中间区域.\end{array} \right.$
式中：e为偏心距，φ为偏位角，θ为轴承周向坐标，如图 5所示.
Figure 5
图 5 新型箔片轴承坐标系 Figure 5 Coordinate system of the journal bearing 箔片轴承量纲一的气膜厚度表示如下：
$\bar h = \left\{ \begin{array}{l}{C_1}/C + \varepsilon \cos \left( {\theta - \varphi } \right) + {{\bar U}_R},\;\;\;两端区域;\\{C_2}/C + \varepsilon \cos \left( {\theta - \varphi } \right) + {{\bar U}_R},\;\;\;中间区域.\end{array} \right.$
其中：
$\varepsilon = e/C,\;\;\;\;{{\bar U}_R} = {U_R}/C.$
2.4 压力控制方程对于小尺寸箔片轴承，分析时可作如下假设：
1) 由于h/L ?1，则y方向的Navier-Stokes方程可简化为?p/?y=0.即压力沿气膜厚度方向的变化可以忽略.
2) 对于标准大气压下的小尺寸轴承，雷诺数通常＜0.05，由量纲分析可知气体惯性力可以忽略.
3) 气体流动遵循无滑移边界条件.
4) 气体流动近似为等温过程.
根据上述假设，可得到如下理想气体的雷诺方程：
$\begin{array}{l}\frac{1}{{{R^2}}}\frac{\partial }{{\partial \theta }}\left( {p{h^3}\frac{{\partial p}}{{\partial \theta }}} \right) + \frac{\partial }{{\partial z}}\left( {p{h^3}\frac{{\partial p}}{{\partial z}}} \right) = \\\;\;\;\;\;\;\;6\mu \omega \frac{\partial }{{\partial \theta }}\left( {ph} \right) + 12\mu \frac{\partial }{{\partial t}}\left( {ph} \right).\end{array}$
式中：p为气膜压力，Pa；R为轴颈半径，m；L为轴承轴向长度，m；μ为气体动力黏度，Pa·s；ω为轴颈转速，rad·s-1；z为轴承轴向坐标，m；t为时间，s.
定义量纲一的参数如下：
$\bar z = \frac{z}{{0.5L}},\bar h = \frac{h}{C},\bar p = \frac{{\bar p}}{{{p_{\rm{a}}}}},\bar t = \tau t.$
式中：z为量纲一轴向坐标；p为量纲一气膜压力；pa为环境压力，Pa；t为量纲一时间；τ为涡动转速，rad·s-1.
最终，可得到量纲一的雷诺方程为
$\begin{array}{l}\frac{\partial }{{\partial \theta }}\left( {\bar p{{\bar h}^3}\frac{{\partial \bar p}}{{\partial \theta }}} \right) + {\left( {\frac{{2R}}{L}} \right)^2}\frac{\partial }{{\partial \bar z}}\left( {\bar p{{\bar h}^3}\frac{{\partial \bar p}}{{\partial \bar z}}} \right) = \\\;\;\;\;\zeta \frac{\partial }{{\partial \theta }}\left( {\bar p\bar h} \right) + 2\gamma \zeta \frac{\partial }{{\partial \bar t}}\left( {\bar p\bar h} \right).\end{array}$ (5)
式中：ζ为轴承数，γ为涡动频率比，且
$\zeta = \frac{{6{R^2}\mu \omega }}{{{p_{\rm{a}}}{C^2}}},\;\;\;\;\;\gamma {\rm{ = }}\frac{\tau }{\omega }.$
3 数值求解箔片轴承的压力控制方程(5)为非线性方程，无法求其精确解.本文采用Newton-Raphson迭代法和有限差分法相结合的方法求解其静特性.首先，定义关于p的非线性函数为
$F\left( {\bar p} \right) = \frac{\partial }{{\partial \theta }}\left( {\bar p{{\bar h}^3}\frac{{\partial \bar p}}{{\partial \theta }}} \right) + {\left( {\frac{{2R}}{L}} \right)^2}\frac{\partial }{{\partial \bar z}}\left( {\bar p{{\bar h}^3}\frac{{\partial \bar p}}{{\partial \bar z}}} \right) - \zeta \frac{\partial }{{\partial \theta }}\left( {\bar p\bar h} \right).$
由Newton-Raphson迭代法可知
$F\left( {{{\bar p}^n}} \right) + F'\left( {{{\bar p}^n}} \right){\delta ^n} = 0\left( {{\delta ^n} = {{\bar p}^{n + 1}} - {{\bar p}^n},n = 0,1,2, \cdots } \right).$
对F(p+βδ)进行Taylor级数展开，进而对β求一阶导数.令β=0，得
$\frac{{\rm{d}}}{{{\rm{d}}\beta }}F\left( {\bar p + \beta \delta } \right)\left| {_{\beta = 0}} \right. = F'\left( {\bar p} \right)\delta = - F\left( {\bar p} \right).$ (6)
式(6)将关于p的非线性函数转化成了关于压力增量δ的线性偏微分方程，有效简化了求解.基于有限差分法，将各偏微分项表示为中心差分格式，最终可得到如下求解方程组：
${a_{i,j}}{\delta _{i - 1,j}} + {b_{i,j}}{\delta _{i + 1,j}} + {c_{i,j}}{\delta _{i,j}} + {d_{i,j}}{\delta _{i,j - 1}} + {e_{i,j}}{\delta _{i,j + 1}} = - {S_{i,j}}.$
式中：i、j分别表示周向和轴向的节点序号，ai, j、bi, j、ci, j、di, j、ei, j、Si, j为系数矩阵对应各元素.
给定初始气膜厚度及压力，通过迭代求解可获得压力和气膜厚度的数值解，进而计算得到该新型箔片轴承的承载力等静态性能：
${{\bar W}_x} = 2\int_0^1 {\int_0^{2{\rm{ \mathsf{ π} }}} {\left( {\bar p - 1} \right)\sin \theta {\rm{d}}\theta {\rm{d}}\bar z} } ,$
${{\bar W}_y} = 2\int_0^1 {\int_0^{2{\rm{ \mathsf{ π} }}} {\left( {\bar p - 1} \right)\cos \theta {\rm{d}}\theta {\rm{d}}\bar z} } ,$
$\bar W = \sqrt {{{\left( {{{\bar W}_x}} \right)}^2} + {{\left( {{{\bar W}_y}} \right)}^2}} ,W = 0.5LR{p_{\rm{a}}}\bar W,$
$\tan \varphi = - {{\bar W}_x}/{{\bar W}_y}.$
式中：Wx、Wy分别为x、y方向量纲一的承载力，W为轴承总的量纲一的承载力，W为轴承总承载力.
采用小扰动法求解新型箔片轴承的动力学特性系数[20].给定量纲一的小扰动位移(ΔX, ΔY)和小扰动速度(Δ?, Δ?)，并对压力、气膜厚度、箔片变形进行Taylor级数展开：
$\left\{ \begin{array}{l}\bar p = {{\bar p}_0} + {{\bar p}_x}\Delta X + {{\bar p'}_x}\Delta \dot X + {{\bar p}_y}\Delta Y + {{\bar p'}_y}\Delta \dot Y,\\\bar h = {{\bar h}_0} + {{\bar h}_x}\Delta X + {{\bar h'}_x}\Delta \dot X + {{\bar h}_y}\Delta Y + {{\bar h'}_y}\Delta \dot Y,\\{{\bar U}_R} = {{\bar U}_{R0}} + {{\bar U}_{Rx}}\Delta X + {{\bar U'}_{Rx}}\Delta \dot X + {{\bar U}_{Ry}}\Delta Y + {{\bar U'}_{Ry}}\Delta \dot Y.\end{array} \right.$ (7)
将方程(7)代入压力控制雷诺方程、箔片变形方程，联立求解各扰动方程，得到新型箔片轴承的刚度、阻尼系数：
$\left\{ {\begin{array}{*{20}{c}}{{{\bar K}_{xx}}}&{{{\bar K}_{xy}}}\\{{{\bar K}_{yx}}}&{{{\bar K}_{yy}}}\end{array}} \right\} = - 2\int_0^1 {\int_0^{2{\rm{ \mathsf{ π} }}} {\left\{ {\begin{array}{*{20}{c}}{{{\bar p}_x}\sin \theta }&{{{\bar p}_y}\sin \theta }\\{ - {{\bar p}_x}\cos \theta }&{ - {{\bar p}_y}\cos \theta }\end{array}} \right\}{\rm{d}}\theta {\rm{d}}\bar z} } ,$
$\left\{ {\begin{array}{*{20}{c}}{{K_{xx}}}&{{K_{xy}}}\\{{K_{yx}}}&{{K_{yy}}}\end{array}} \right\} = \frac{{0.5LR{p_{\rm{a}}}}}{C}\left\{ {\begin{array}{*{20}{c}}{{{\bar K}_{xx}}}&{{{\bar K}_{xy}}}\\{{{\bar K}_{yx}}}&{{{\bar K}_{yy}}}\end{array}} \right\},$
$\left\{ {\begin{array}{*{20}{c}}{{{\bar C}_{xx}}}&{{{\bar C}_{xy}}}\\{{{\bar C}_{yx}}}&{{{\bar C}_{yy}}}\end{array}} \right\} = - 2\int_0^1 {\int_0^{2{\rm{ \mathsf{ π} }}} {\left\{ {\begin{array}{*{20}{c}}{{{\bar p'}_x}\sin \theta }&{{{\bar p'}_y}\sin \theta }\\{ - {{\bar p'}_x}\cos \theta }&{ - {{\bar p'}_y}\cos \theta }\end{array}} \right\}{\rm{d}}\theta {\rm{d}}\bar z} } ,$
$\left\{ {\begin{array}{*{20}{c}}{{C_{xx}}}&{{C_{xy}}}\\{{C_{yx}}}&{{C_{yy}}}\end{array}} \right\} = \frac{{0.5LR{p_{\rm{a}}}}}{{C\tau }}\left\{ {\begin{array}{*{20}{c}}{{{\bar C}_{xx}}}&{{{\bar C}_{xy}}}\\{{{\bar C}_{yx}}}&{{{\bar C}_{yy}}}\end{array}} \right\}.$
式中：Kij、Cij分别为箔片轴承量纲一的刚度、阻尼系数，Kij、Cij分别为箔片轴承刚度、阻尼系数.
最终，得到新型箔片轴承静、动态性能的仿真分析流程如图 6所示.
Figure 6
图 6 数值求解流程 Figure 6 Flow chart of numerical computation 4 仿真结果及分析表 1列出了本文新型波箔型箔片轴承的结构参数.对于波箔型箔片轴承，波箔片的局部变形将使波箔片与平箔片的接触面积增大，进而导致平箔片的刚度也有所增大.文献[21]中通过引入平箔片沿周向方向的刚化系数来进行修正.令该系数为4，则平箔片沿周向方向的弹性模量修正为856 GPa.
表 1
表 1 新型箔片轴承结构参数 Table 1 Parameters of the novel foil journal bearing 轴承长度/ mm轴颈直径/
mm表面宽槽轴向长度/
mm名义平均半径间隙/
μm箔片材料箔片弹性模量/
GPa箔片泊松比
1616.96812Inconel 7182140.29
转速/
(r·min-1)波箔高度/
mm平箔、波箔厚度/
mm波箔片单元长度/
mm波箔片半波纹长度/
mm气体黏度/
(Pa·s)环境压力/
Pa
50 0000.50.101 64.21.751.932×10-5101 325
表 1 新型箔片轴承结构参数 Table 1 Parameters of the novel foil journal bearing
考虑到该新型箔片轴承关于轴承中截面对称，因而在求解其静、动态特性时只取轴承轴向的一半.同时，对平箔片二维厚板刚度和雷诺方程的求解采用相同的网格划分，周向和轴向的单元数为72×14.
4.1 承载力分析图 7给出了不同偏心距下新型箔片轴承承载力与平箔片中间区域宽槽深度Td的关系曲线.轴承工作转速为50 000 r·min-1.可以看出，随着平箔片宽槽深度Td增大，轴承的承载力有显著增大，且偏心率越大，轴承承载力的增大效果越显著.对于偏心距16 μm、平箔片宽槽深度8 μm的新型箔片轴承，其承载力比传统箔片轴承增大约11.8 %.
Figure 7
图 7 轴承承载力-平箔片宽槽深度曲线 Figure 7 Load capacity versus wide groove depth of top foil 为研究该新型箔片轴承承载力增大的机理，图 8~11对比分析了新型箔片轴承和传统箔片轴承气膜厚度分布及压力分布结果.其中，轴承偏心距为16 μm，新型箔片轴承的平箔片宽槽深度为8 μm，工作转速均为50 000 r·min-1.对于传统箔片轴承，其气膜厚度沿轴向方向基本为一定值.而新型箔片轴承，由于平箔片中间区域宽槽的影响，导致轴承中间区域的气膜厚度有所增大，而两端区域的气膜厚度减小.该结构特点有利于润滑气体在中间区域聚集，减少气体端泄，从而提高润滑气体压力，提高轴承承载力.由图 10可以看出，新型箔片轴承在中间区域的气膜压力均处于较高的水平.另一方面，该区域气膜压力的增大又会反作用于箔片结构，使得平箔片在相邻波箔之间“下凹”变形更加明显，导致局部的气体动压效果更为突出，从而增大润滑气体局部压力，提高轴承承载力水平.
Figure 8
图 8 平箔片厚度轴向变化的箔片轴承气膜厚度分布(Td=8 μm, e=16 μm) Figure 8 Film thickness distribution of the novel foil bearing with top foil thickness variation(Td=8 μm, e=16 μm) Figure 9
图 9 传统箔片轴承气膜厚度分布(Td=0 μm, e=16 μm) Figure 9 Film thickness distribution of the traditional gas foil bearing(Td=0 μm, e=16 μm) Figure 10
图 10 平箔片厚度轴向变化的箔片轴承压力分布(Td=8 μm, e=16 μm) Figure 10 Pressure distribution of the novel foil bearing with top foil thickness variation(Td=8 μm, e=16 μm) Figure 11
图 11 传统箔片轴承压力分布(Td=0 μm, e=16 μm) Figure 11 Pressure distribution of the traditional gas foil bearing(Td=0 μm, e=16 μm) 箔片轴承润滑气体的端泄流量可根据润滑气体沿轴向方向速度分布积分得到，即
$\begin{array}{l}{Q_{{\rm{leak}}}} = 2\int_0^{{\theta _{{\rm{cr}}}}} {\left( {R\int_0^h {{v_z}\left| {_{z = 0.5L}{\rm{d}}y} \right.} } \right){\rm{d}}\theta } = \\\;\;\;\;\;\;\;\;\;\; - 2\int_0^{{\theta _{{\rm{cr}}}}} {\left( {\frac{{R{h^3}}}{{12\mu }}\frac{{\partial p}}{{\partial z}}\left| {_{z = 0.5L}} \right.} \right){\rm{d}}\theta } .\end{array}$ (8)
式中θcr为润滑气体压力变为负压的临界值.
利用公式(8)可计算得到不同平箔片宽槽深度Td下新型箔片轴承的气体端泄流量，轴承偏心距为16 μm，结果见表 2.分析可知，虽然平箔片宽槽深度Td的增大会使轴承压力增大，压力沿轴向的梯度?p/?z也随之增大，但此时轴承两端边缘处润滑气膜厚度的减小起主导作用.因此，导致新型箔片轴承气体的端泄量随宽槽深度Td的增大逐渐减小，从而进一步验证了该新型箔片轴承承载力增大的机理.
表 2
表 2 新型箔片轴承气体端泄量(e=16 μm) Table 2 Side leakage flow of the novel foil bearing(e=16 μm) 宽槽深度/ μm气体端泄流量/ (10-4 m3·min-1)02.62
22.25
41.99
61.82
81.68
101.56
121.46
141.36
161.27
表 2 新型箔片轴承气体端泄量(e=16 μm) Table 2 Side leakage flow of the novel foil bearing(e=16 μm)
图 12给出了不同平箔片宽槽深度Td下轴承最大压力沿轴向方向的变化曲线，图中箭头方向表示平箔片中间区域宽槽深度Td在增大.对于传统箔片轴承，轴承最大压力在中间截面处达到最大，由中间向两端逐渐减小，并在轴承两端面减小到环境压力.随着平箔片宽槽深度Td的增大，轴承最大压力在平箔片的宽槽边界附近急剧增大.这主要是由于气膜厚度在该边界处的突变造成的，气膜厚度的变化越剧烈，轴承压力就增大得越为明显.
Figure 12
图 12 轴承最大压力沿轴向方向分布(e=16 μm) Figure 12 Pressure at line of maximum pressure along bearing length direction(e=16 μm) 4.2 动力学特性基于小扰动法，计算得到新型箔片轴承刚度系数随平箔片宽槽深度Td的关系如图 13所示，其中轴承静态载荷为16 N，转速50 000 r·min-1，涡动频率比取为1.可以看出，与主刚度相比，新型箔片轴承的交叉刚度约低1个数量级，且两个方向的交叉刚度正负号相反.同时，随着平箔片宽槽深度Td增大，轴承的气膜压力水平提高，导致其主刚度系数也逐渐增大，而交叉刚度的绝对值逐渐减小，这意味着轴承的稳定性有所提高.对于静态载荷为16 N、平箔片宽槽深度为8 μm的新型箔片轴承，其y方向的主刚度系数比传统箔片轴承增大约12.1 %.
Figure 13
图 13 刚度系数-平箔片宽槽深度曲线(轴承静载：16 N) Figure 13 Stiffness coefficients versus groove depth of top foil(Static load: 16 N) 轴承阻尼系数随平箔片宽槽深度变化见图 14.新型箔片轴承的交叉阻尼也小于主阻尼，且两个方向交叉阻尼正负号相反.随着平箔片宽槽深度Td增大，轴承x方向的主阻尼呈增大趋势，而y方向的主阻尼逐渐减小.当平箔片宽槽深度Td在0~16 μm变化时，两个方向交叉阻尼的变化均很小.对于静态载荷为16 N、平箔片宽槽深度为8 μm的新型箔片轴承，其主阻尼系数比传统箔片轴承增大约11.0 %.
Figure 14
图 14 阻尼系数-平箔片宽槽深度曲线(轴承静载：16 N) Figure 14 Damping coefficients versus groove depth of top foil(Static load: 16 N) 图 15给出了不同平箔片宽槽深度Td下新型箔片轴承各刚度系数随转速的变化曲线.轴承静态载荷取为16 N，涡动频率比为1.可以看出，随着转速增大，轴承的主刚度系数均显著增大，而交叉刚度系数明显减小.在0~120 000 r·min-1转速内，平箔片宽槽深度Td越大，轴承的主刚度系数也越大，同时随转速变化得也更加剧烈.
Figure 15
图 15 轴承刚度系数-转速曲线(轴承静载：16 N) Figure 15 Stiffness coefficients versus rotating speed(Static load: 16 N) 不同转速下的新型箔片轴承阻尼系数如图 16所示.轴承运行参数与图 15相同.随着转速增大，x方向的主阻尼呈现先增大后减小的趋势，而y方向的主阻尼逐渐减小.与轴承刚度系数的变化规律类似，平箔片宽槽深度Td越大，轴承的主阻尼系数也越大，同时随转速变化得也更加剧烈.
Figure 16
图 16 轴承阻尼系数-转速曲线(轴承静载：16 N) Figure 16 Damping coefficients versus rotating speed(Static load: 16 N) 5 结论本文提出了一种新型的平箔片厚度轴向变化的波箔型箔片轴承.为考虑平箔片厚度的轴向变化，采用二维厚板有限元模型计算平箔片刚度.仿真计算了新型箔片轴承承载力、动力学特性等性能，研究了平箔片中间区域宽槽深度对新型箔片轴承静、动态性能的影响规律.主要结论如下：
1) 与传统箔片轴承相比，平箔片厚度轴向变化的箔片轴承承载力有所增大，且平箔片宽槽深度越大，轴承承载力增加越显著.这是由于平箔片中间区域存在凹槽，导致新型箔片轴承的中间区域气膜厚度增大，而两端气膜厚度减小，有利于润滑气体的聚集，减小气体端泄，从而增大流体动压水平，提高承载力.
2) 在0~120 000 r/min的转速范围内，平箔片厚度轴向变化的箔片轴承主刚度和主阻尼比传统箔片轴承要大.同时，平箔片宽槽深度越大，主刚度和主阻尼增加越显著，从而箔片轴承的稳定性得以提高.
3) 在较低转速范围内增大转速时，箔片轴承的主刚度和主阻尼有所增大，交叉刚度逐渐减小，说明在一定的转速范围内，增大转速有利于转子稳定性的提高.
4) 对于平箔片宽槽深度为8 μm的新型箔片轴承，其承载力比传统箔片轴承增大约11.8 %，主刚度增大约12.1 %.
参考文献
[1] AGRAWAL G L. Foil air/gas bearing technology-an overview[C]// ASME 1997 International Gas Turbine and Aeroengine Congress and Exhibition. Orlando: International Gas Turbine Institute, 1997: V001T04A006. DOI: 10.1115/97-GT-347.
[2]DELLACORTE C. Oil-free shaft support system rotordynamics: past, present and future challenges and opportunities[J].Mechanical Systems and Signal Processing, 2012, 29: 67-76.DOI: 10.1016/j.ymssp.2011.07.024
[3]BLOK H, Van ROSSUM J J. The foil bearing-a new departure in hydrodynamic lubrication[J].ASLE Journal of Lubrication Engineering, 1953, 9(6): 316-320.
[4]HESHMAT C A, HESHMAT H. An analysis of gas-lubricated, multileaf foil journal bearings with backing springs[J].ASME Journal of Tribology, 1995, 117(7): 437-443.DOI: 10.1115/1.2831272
[5]ARAKERE N K. Analysis of foil journal bearings with backing springs[J].Tribology Transactions, 1996, 39(1): 208-214.DOI: 10.1080/10402009608983522
[6] HESHMAT H. Multi-stage support element for compliant hydrodynamic bearings: US 4300806[P]. 1981-11-17.
[7] GRAY S, BHUSHAN B. Support element for compliant hydrodyna-mic journal bearings: US 4274683[P]. 1981-06-23.
[8]HESHMAT H, SHAPIRO W, GRAY S. Development of foil journal bearings for high load capacity and high speed whirl stability[J].Journal of Lubrication Technology, 1982, 104(2): 149-156.DOI: 10.1115/1.3253173
[9] HESHMAT H. High load capacity compliant foil hydrodynamic journal bearing: US 5902049[P]. 1999-05-11.
[10]SIM K, LEE Y B, KIM T H. Rotordynamic analysis of an oil-free turbocharger supported on lobed gas foil bearings-predictions versus test data[J].Tribology Transactions, 2014, 57(6): 1086-1095.DOI: 10.1080/10402004.2014.937885
[11]SIM K, LEE J, LEE Y B, et al. Rotordynamic performance of shimmed gas foil bearings for oil-free turbochargers[J].ASME Journal of Tribology, 2012, 134(3): 31102.DOI: 10.1115/1.4005892
[12]SONG J, KIM D. Foil gas bearing with compression springs: analyses and experiments[J].ASME Journal of Tribology, 2007, 129(3): 628-639.DOI: 10.1115/1.2736455
[13]SAN ANDRES L, CHIRATHADAM T A. A metal mesh foil bearing and a bump-type foil bearing: comparison of performance for two similar size gas bearings[J].Journal of Engineering for Gas Turbines and Power, 2012, 134(10): 102501.DOI: 10.1115/1.4007061
[14]SAN ANDRES L, CHIRATHADAM T A. Performance characteristics of metal mesh foil bearings: predictions versus measurements[J].Journal of Engineering for Gas Turbines and Power, 2013, 135(12): 122503.DOI: 10.1115/1.4025146
[15]赖天伟, 郑越青, 陈双涛, 等. 多层鼓泡弹性支承箔片动压气体轴承的实验研究[J].东北大学学报(自然科学版), 2014, 35(3): 415-418.
LAI Tianwei, ZHENG Yueqing, CHEN Shuangtao, et al. Experimental study on the performance of multi-decked elastic support protuberant foil gas bearing[J].Journal of Northeastern University (Natural Science), 2014, 35(3): 415-418.
[16]FENG K, LIU Y, ZHAO X, et al. Experimental evaluation of the structure characterization of a novel hybrid bump-metal mesh foil bearing[J].ASME Journal of Tribology, 2016, 138(2): 21702.DOI: 10.1115/1.4031496
[17]DELLACORTE C, VALCO M J. Load capacity estimation of foil air journal bearings for oil-free turbomachinery applications[J].Tribology Transactions, 2000, 43(4): 795-801.DOI: 10.1080/10402000008982410
[18]刘占生, 徐方程, 张广辉, 等. 基于二维厚板模型的波箔片轴承静特性[J].航空动力学报, 2012, 27(6): 1405-1415.
LIU Zhansheng, XU Fangcheng, ZHANG Guanghui, et al. Static characterization of bump-type gas foil bearing: integration of top foil 2-D thick plate model[J].Journal of Aerospace Power, 2012, 27(6): 1405-1415.
[19]IORDANOFF I. Analysis of an aerodynamic compliant foil thrust bearing: Method for a rapid design[J].ASME Journal of Tribology, 1999, 121(4): 816-822.DOI: 10.1115/1.2834140
[20]PENG J P, CARPINO M. Calculation of stiffness and damping coefficients for elastically supported gas foil bearings[J].ASME Journal of Tribology, 1993, 115(1): 20-27.DOI: 10.1115/1.2920982
[21] SAN ANDRES L, KIM T H. Improvements to the analysis of gas foil bearings: integration of top foil 1D and 2D structural models[C]// ASME Turbo Expo 2007: Power for Land, Sea, and Air. Montreal: International Gas Turbine Institute, 2007: 779-789. DOI: 10.1115/GT2007-27249.
