滑梁式空气箔片轴承是一种具有发展潜力的轴承结构. 与波箔型空气箔片轴承相比，其制作工艺简单，成本较低. 鉴于现阶段对其力学及润滑性能的理论研究尚不够深入，研究结构参数对滑梁式空气箔片轴承润滑性能的影响规律十分必要.

目前对空气箔片轴承的研究主要集中于波箔型结构. Walowit等 ^[1]首次构建箔片弹性变形模型并求得弹性系数，结合一维Reynolds方程研究轴承的静特性. Heshmat等^[2]提出波箔的线性弹簧等效模型，利用有限差分法计算轴承在不同尺寸和运行参数下的静特性. Carpino等^[3-4]在考虑箔片弯曲、润滑膜和弹性支承共同作用的基础上，求解Reynolds方程，得到气膜压力分布情况，利用有限元法预测箔片结构变形量. 许怀锦等^[5-6]通过数值方法求解压力控制方程和气膜厚度方程，对比研究波箔型轴承和刚性轴承的润滑性能. Kim等^[7]借助直接隐式积分法求解箔片轴承的非线性动力模型，研究轴承在不同设计参数和运行条件下对质量不平衡激励的响应. Ku等^[8]用光学跟踪系统测量箔片轴承在各种运行条件下的波箔二维变形量，发现设计参数会影响波箔结构特性. Radil等^[9]试验研究半径间隙对波箔型轴承气膜承载力的影响及结构参数对承载能力和热学性能的影响. Li等^[10]通过试验研究非均匀波箔型轴承支承的空压机轴颈的动力学特性，发现减小半径间隙可以抑制轴颈振动. 许浩杰等^[11]综合考虑箔片结构变形、可压缩气体Reynolds方程、能量方程和表面粗糙度，研究空气箔片径向轴承的静动特性，同时分析箔片结构尺寸参数如何影响波箔轴承静动特性.

波箔的加工生产对模具的依赖度高，而轴承性能对波箔尺寸变化敏感，这使得模具加工的精度要求很高，导致生产成本居高不下. 随着激光切割技术的发展，滑梁式箔片轴承的优势逐渐显现. 由于滑梁箔片只需通过激光切割即可满足高精度加工要求，与常见的波箔型轴承相比，滑梁式箔片轴承具有显著的经济性优势. 然而，目前公开文献中对其性能的研究仍处于探索阶段. 冯凯等^[12]研究三瓣式弹性箔片轴承，用数值研究证明该结构可以提高轴承稳定性和承载力. 李长林等^[13-15]对三瓣结构的滑梁式轴承进行研究，计算滑梁形变并对该轴承的静动特性及其支撑的转子系统的性能进行分析. 吴洋等^[16]采用非线性薄板单元建立滑梁箔片力学模型，并对该轴承静态特性进行精确计算. 但这些研究均针对三瓣结构的滑梁式空气箔片轴承. 顾晨昀等^[17]对整周结构的滑梁式空气轴承的疲劳寿命进行理论研究，但未进一步对其承载能力和润滑性能进行分析.

在前人研究的基础上，以梯形滑梁结构的空气箔片轴承为研究对象，应用力学和流体润滑理论，构建流固耦合计算模型. 通过数值计算的方法，深入研究结构参数对其润滑性能的影响规律.

图1为滑梁式空气箔片轴承的基本结构. 其中，轴承的顶箔和底箔均为一片式的完整箔片，为整周结构. 底箔为滑梁箔片，其主要参数如图2所示.滑梁在底箔片上按照一定的规律排列，其形状由长度L₀、宽度d和梯形斜率γ决定；L_C为箔片周向长度，L_B为箔片轴向长度，m和n分别为周向和轴向滑梁个数. 当滑梁箔片弯曲时，箔片上的滑梁将保持直线形状，并在安装后与轴套内表面接触，起到支撑顶箔的作用，如图3所示.其中，R为轴颈半径，R_s为顶箔安装后的圆弧半径，O和O'分别为轴套和转子的中心，e₀为转子的偏心距，ψ₀为偏位角，ω为转子旋转的角速度. 由于收敛间隙的存在，轴颈在高速转动时产生气膜动压，此时，底箔上的滑梁弯曲，滑梁的两端会在轴套上滑动.

由于箔片的厚度通常不超过0.3 mm，可以采用Kirchhoff理论，使用薄板单元对箔片进行力学建模. 顶箔为规则的矩形结构，故采用图4中的矩形单元建立力学模型. 图中，x和y为单元内的局部坐标.

基于薄板理论，挠度w与转角θ_x、θ_y满足

(1)$ \left. {\begin{gathered} {{\theta _x} = {{\partial w} \mathord{\left/ {\vphantom {{\partial w} {\partial y}}} \right. } {\partial y}},} \\ {{\theta _y} = - {{\partial w} \mathord{\left/ {\vphantom {{\partial w} {\partial x}}} \right. } {\partial x}}.} \\\end{gathered}} \right\} $

设单元e内节点i的挠度及x和y向弯角分别为$ w_i^e $、$ \theta _{{x_i}}^e $和$ \theta _{{y_i}}^e $，利用平面Hermite插值，可以得到单元内挠度w^e的表达式：

(2)$ {w^e} = {{\boldsymbol{N}}^e}{{\boldsymbol{\alpha }}^e} = \sum\limits_{i = 1}^4 {{\boldsymbol{N}}_i^e{\boldsymbol{\alpha }}_i^e} . $

式中：${\boldsymbol{\alpha }}^e $为该单元的位移阵列，N^e为对应的插值形函数， $ {\boldsymbol{\alpha }}_i^e = {\left[ {\begin{array}{*{20}{c}} {w_i^e},&{\theta _{{x_i}}^e},&{\theta _{{y_i}}^e} \end{array}} \right]^{\text{T}}} $.

单元e内的弯曲应变κ^e可以表示为

(3)$ \begin{split} &{{\boldsymbol{\kappa }}^e} =\\& \left[ {\begin{array}{*{20}{c}} { - {\partial ^2}w^e{\text{/}}\partial {x^2}} \\ { - {\partial ^2}w^e{\text{/}}\partial {y^2}} \\ { - 2{\partial ^2}w^e{\text{/}}\partial x\partial y} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - {\partial ^2}{{\boldsymbol{N}}^e}{\text{/}}\partial {x^2}} \\ { - {\partial ^2}{{\boldsymbol{N}}^e}{\text{/}}\partial {y^2}} \\ { - 2{\partial ^2}{{\boldsymbol{N}}^e}{\text{/}}\partial x\partial y} \end{array}} \right]{{\boldsymbol{\alpha }}^e} = {{\boldsymbol{B}}^e}{{\boldsymbol{\alpha }}^e}. \end{split}$

式中：B^e为单元应变矩阵. 单元的总位能$ \varPi _{\mathrm{d}}^e $为

(4)$ \varPi _{\text{d}}^e = {\left( {{{\boldsymbol{\alpha }}^e}} \right)^{\text{T}}}\left\{ {\frac{1}{2}\left[\int {{\left( {{{\boldsymbol{B}}^e}} \right)}^{\text{T}}}{\boldsymbol{D}}{{\boldsymbol{B}}^e}{\mathrm{d}}x\right]{{\boldsymbol{\alpha }}^e} - {{\boldsymbol{P}}^e}} \right\}. $

式中： P^e为对应的单元载荷阵列；D为弹性矩形，且

$ {\boldsymbol{D}} = \dfrac{{E{t^3}}}{{12\left( {1 - {\mu ^2}} \right)}}\left[ {\begin{array}{*{20}{c}} 1&\mu &0 \\ \mu &1&0 \\ 0&0&\dfrac{1 - \mu }{2} \end{array}} \right],$

其中E、μ和t分别为箔片的弹性模量、泊松比和厚度. 基于最小位能原理，式(4)的变分为0，化简后得到单元e的平衡方程：

(5)$ {{\boldsymbol{K}}^e}{{\boldsymbol{\alpha }}^e} = \left[\int {{\left( {{{\boldsymbol{B}}^e}} \right)}^{\text{T}}}{\boldsymbol{D}}{{\boldsymbol{B}}^e}{\mathrm{d}}x\right]{{\boldsymbol{\alpha }}^e} = {{\boldsymbol{P}}^e} . $

式中：K^e为单元e的刚度矩阵.

计算顶箔中所有单元的${{\boldsymbol{\alpha }}^e} $、K^e和P^e并进行集成，得到顶箔整体的位移阵列${\boldsymbol{\alpha }}_{{\mathrm{top}}} $、整体刚度矩阵K_top和整体载荷阵列P_top，并满足

(6)$ {{\boldsymbol{ K}}_{\text{top}}}{{\boldsymbol{\alpha }}_{\text{top}}} = {{\boldsymbol{P}}_{\text{top}}}. $

考虑到滑梁箔片复杂的几何结构，在进行力学分析时需要进行适当的简化. 如图5所示，将滑梁箔片划分为3个部分，分别为滑梁、连接梁和边缘部分. 其中，边缘部分为规则的矩形结构，可以通过式(5)建立其单元刚度矩阵，并在集成后得到边缘部分的结构刚度矩阵K_b. 连接梁又分为周向部分和轴向部分，两者存在交点T. 滑梁和轴向连接梁相交于Q点.

如图6所示，在交点T处，周向梁与轴向梁有相同的挠度w. 此外，两者的扭角与弯角应当是耦合的，即周向梁扭角$ \varphi _{{T}}^{\text{C}} $与轴向梁弯角$ \theta _{{T}}^{\text{A}} $相等，周向梁弯角$ \theta _{{T}}^{\text{C}} $与轴向梁扭角$ \varphi _{{T}}^{\text{A}} $相等.因此，需要同时考虑连接梁的挠度w、弯角θ和扭角φ，其中，$ \theta = {{\partial w} \mathord{\left/ {\vphantom {{\partial w} {\partial l}}} \right. } {\partial l}} $，l为梁的长度. 将连接梁离散为两点形式的梁单元，并设定单元e 的节点位移为$ w_i^e $、$ \theta _i^e $和$ \varphi _i^e $（i=1, 2）. 分别通过Hermite插值和Lagrange插值得到连接梁单元e的位移w^e和φ^e的函数：

(7)$ {{\boldsymbol{w}}^e_l} = {\boldsymbol{N}}_l^e{{\boldsymbol{\beta }}^e} = \sum\limits_{i = 1}^2 {{\boldsymbol{N}}_{{l_i}}^e{\boldsymbol{\beta }}_i^e} . $

式中：${\boldsymbol{\beta }}^e $为连接梁单元的位移阵列；${\boldsymbol{N}}_l^e $为对应的形函数矩阵；$ {{\boldsymbol{w}}^e_l} = {\left[ {{w^e}},\;{{\varphi ^e}} \right]^{\text{T}}} $；$ {\boldsymbol{\beta }}_i^e = {\left[ {w_i^e},\;{\theta _i^e},\;{\varphi _i^e} \right]^{\text{T}}} $；$ {\boldsymbol{N}}_{l_i}^e = \left[ {{{{\boldsymbol{N}}}}_{{\theta _i}}^e},\;{N_{{\varphi _i}}^e} \right] $，其中$ {{{\boldsymbol{N}}}}_{{\theta _i}}^e $与$ N_{{\varphi _i}}^e $分别为弯曲与扭转的插值形函数.

$ {{\boldsymbol{w}}^e_l} $随l的变化率为

(8)$ {{\partial {{\boldsymbol{w}}^e_l}} \mathord{\left/ {\vphantom {{\partial {{\boldsymbol{w}}^e}} {\partial l}}} \right. } {\partial l}} = \left[ {\begin{array}{*{20}{c}} {{{\partial {w^e}} \mathord{\left/ {\vphantom {{\partial {w^e}} {\partial l}}} \right. } {\partial l}}} \\ {{{\partial {\varphi ^e}} \mathord{\left/ {\vphantom {{\partial {\varphi ^e}} {\partial l}}} \right. } {\partial l}}} \end{array}} \right] = {\boldsymbol{B}}_l^e{{\boldsymbol{\beta }}^e}. $

式中：$ {\boldsymbol{B}}_l^e $为广义单元应变矩阵. 结合梁的弯曲与扭转变形公式，可以得到连接梁单元的总位能$ \varPi _l^e $为

(9)$ \varPi _l^e = {\left( {{{\boldsymbol{\beta }}^e}} \right)^{\text{T}}}\left\{ \frac{1}{2}\left[{{\int {\left( {{\boldsymbol{B}}_l^e} \right)} }^{\text{T}}}{{\boldsymbol{D}}_l}{\boldsymbol{B}}_l^e{\mathrm{d}}l\right] {{\boldsymbol{\beta }}^e} - {\boldsymbol{P}}_l^e\right\}. $

式中：广义弹性矩阵$ {{\boldsymbol{D}}_l} = {\mathrm{diag}}\left( {E{I_{\text{z}}},G{I_{\text{p}}}} \right) $ ，其中I_z和I_p分别为连接梁截面的中性轴惯性矩和极惯性矩，G为箔片的切变模量；$ {\boldsymbol{P}}_l^e $为连接梁单元的载荷阵列. 基于最小位能原理，式(9)的变分为0，化简后得到：

(10)$ {\boldsymbol{K}}_l^e{{\boldsymbol{\beta }}^e} = \left[\int \left( {{\boldsymbol{B}}_l^e} \right) ^{\text{T}}{{\boldsymbol{D}}_l}{\boldsymbol{B}}_l^e{\mathrm{d}}l \right] {\boldsymbol{\beta }}^e = {\boldsymbol{P}}_l^e. $

式中：$ {\boldsymbol{K}}_l^e $为连接梁单元的刚度矩阵.

建立所有连接梁单元的$ {\boldsymbol{K}}_l^e $，结合周向梁与轴向梁处的耦合条件，在集成后得到连接梁整体的刚度矩阵K_l.

滑梁在结构中起到支撑的作用，其两端与轴套接触，如图7所示. 图中，Q点为滑梁与轴向连接梁的连接点，也是滑梁对顶箔的支撑点. 因此，可以通过分析滑梁在Q点的挠度与载荷关系，将序号为L的滑梁等效为刚度为k_L的支撑弹簧.当滑梁弯曲时，其两端与轴套接触，滑梁与轴套间存在法向接触力F_N和摩擦力F_f，滑动摩擦系数为f，取f=0.3. 经过坐标转换后，滑梁的端部受到轴向力F_ξ和横向力F_η的作用，并且滑梁的水平位移u和竖直挠度w是耦合的.

若采用线形梁单元进行力学分析，则无法考虑F_ξ和u对挠度w的影响，这显然是不合理的. 因此，结合梁的大变形特点对滑梁的挠度进行数值研究，并对滑梁形变做以下假设：

1）在弹性范围内；

2）为纯弯曲，忽略切应力和拉伸或压缩导致的形变.

根据几何关系和本构关系，滑梁的曲率半径ρ与弯曲内力M满足

(11)$ \frac{1}{\rho } = \frac{{{\text{d}}\theta }}{{{\text{d}}s}} = \frac{M}{{E{I_{\text{z}}}}}. $

式中：s为滑梁的长度. 如图8所示，形变后滑梁的位移u、w、θ与s之间满足几何关系：

(12)$ \left. \begin{gathered} {\mathrm{d}}w = {\mathrm{d}}s \cdot \sin\; \theta , \\ {\mathrm{d}}u = {\mathrm{d}}s \cdot \left( {\cos\; \theta - 1} \right). \\ \end{gathered} \right\} $

由式(12)得到：$ \theta = \arcsin \left(\dfrac{{{\mathrm{d}}w}}{{{\mathrm{d}}s}}\right) = \arcsin w' $，代入式(11)得到：

(13)$ \frac{{w''}}{{\sqrt {1 - w{'^2}} }}{\text{ = }}\frac{M}{{E{I_{\text{z}}}}} = \frac{{{\mathrm{d}}\theta }}{{{\mathrm{d}}s}}. $

将滑梁离散为两点形式的梁单元，结合式(7)，以节点位移$ w_i^e $、$ \theta _i^e $得到单元内挠度w的插值函数：

(14)$ w = \sum\limits_{i = 1}^2 {{\boldsymbol{N}}_{{\theta _i}}^e{{\left[ {\begin{array}{*{20}{c}} {{w_i}}，{{\theta _i}} \end{array}} \right]}^{\text{T}}}} = {\boldsymbol{N}}_s^e{\boldsymbol{w}}_s^e. $

式中：$ {\boldsymbol{w}}_s^e $ 和$ {\boldsymbol{N}}_s^e $分别为滑梁单元的节点位移阵列和对应的形函数. 令$ A = \dfrac{1}{{1 - w{'^2}}} $，计算滑梁单元的应变能$ \varPi _{{s_0}}^e $：

(15)$ \varPi _{{s_0}}^e = \frac{1}{2}\int M {\mathrm{d}}\theta = \frac{1}{2}\int {E{I_{\text{z}}}} Aw'{'^2}{\mathrm{d}}s. $

由于滑梁为变截面梁，式(15)中I_z可以采用线性插值形式. 滑梁单元的总位能$ \varPi _s^e $为

(16)$ \begin{split} \varPi _s^e = &\varPi _{{s_0}}^e - \int {\left( {\overline {F_\eta ^e} {\mathrm{d}}w+\overline {{M^e}} {\mathrm{d}}\theta } \right)} - \sum\limits_{i = 1}^2 {\left( {u_i^e - u_1^e} \right)F_{{\xi _i}}^e}= \\& \varPi _{{s_0}}^e - {\left( {{\boldsymbol{w}}_s^e} \right)^{\text{T}}}{\boldsymbol{P}}_s^e - F_{{\xi _2}}^e\int {{\mathrm{d}}u} . \end{split} $

式中：$ \overline {F_\eta ^e} $ 和$ \overline {{M^e}} $分别为滑梁单元上横向力和弯矩的分布函数，与对应的形函数相乘并积分后得到载荷阵列$ {\boldsymbol{P}}_s^e $； $ u_i^e $和$ F_{{\xi _i}}^e $分别为节点i的轴向位移和载荷. 对式(16)求变分得到：

(17)$ \text{δ} \varPi _s^e = \text{δ} \varPi _{{s_0}}^e - {\left( {\text{δ} {\boldsymbol{w}}_s^e} \right)^{\text{T}}}{\boldsymbol{P}}_s^e - F_{{\xi _2}}^e\int {{\mathrm{d}}\left( {\text{δ} u} \right)} . $

并且满足

(18)$ \begin{split} \text{δ} {\varPi }_{{s}_{0}}^{e}=&{\displaystyle \int E{I}_{\text{z}}\left(Aw''\cdot \text{δ} w''+{A}^{2}w{''}^{2}w'\cdot \text{δ} w'\right)}{\mathrm{d}}s\text=\\ &{\left(\text{δ} {{\boldsymbol{w}}}_{s}^{e}\right)}^{\text{T}}{{\boldsymbol{K}}}_{{s}_{0}}^{e}{{\boldsymbol{w}}}_{s}^{e},\end{split} $

(19)$ \begin{split} F_{{\xi _2}}^e\int {{\mathrm{d}}\left( {\text{δ} u} \right)} =& - F_{{\xi _2}}^e\int {\sin \theta \,\,\text{δ} \theta } \,\,{\mathrm{d}}s = - F_{{\xi _2}}^e\int {{A^{{1}/{2}}}} \\& w'\text{δ} w' {\mathrm{d}}s= - {\left( {\text{δ} {\boldsymbol{w}}_s^e} \right)^{\text{T}}}{\boldsymbol{K}}_{{s_F}}^e{\boldsymbol{w}}_s^e. \end{split} $

式中：$ {\boldsymbol{K}}_{{s_0}}^e = {\boldsymbol{K}}_{{s_1}}^e + {\boldsymbol{K}}_{{s_2}}^e ， $$ {\boldsymbol{K}}_{{s_1}}^e = \displaystyle\int {E{I_{\text{z}}}A{{\boldsymbol{C}}_2}} {\mathrm{d}}s, $$ {{\boldsymbol{K}}}_{{s}_{2}}^{e} = \displaystyle \int E{I}_{\text{z}}w{''}^{2} \times {A}^{2}{{\boldsymbol{C}}}_{1}{\mathrm{d}}s,\; $$ {\boldsymbol{K}}_{{s_F}}^e = \displaystyle\int {F_{{\xi _2}}^e{A^{{1}/{2}}}{{\boldsymbol{C}}_1}{\mathrm{d}}s} ,$ C₁=$ {\left( {{\boldsymbol{N}}_s^{e \prime}} \right)^{\text{T}}}{\boldsymbol{N}}_s^{e \prime}, $ C₂=$ {\left({{\boldsymbol{N}}}_{s}^{e \prime \prime}\right)}^{\text{T}} {{\boldsymbol{N}}}_{s}^{e \prime \prime} $.

将式(18)和式(19)代入式(17)，根据最小位能原理，梁单元总位能的变分为0，化简后得到：

(20)$ {\boldsymbol{K}}_s^e{\boldsymbol{w}}_s^e = \left( {{\boldsymbol{K}}_{{s_1}}^e+{\boldsymbol{K}}_{{s_2}}^e+{\boldsymbol{K}}_{{s_F}}^e} \right){\boldsymbol{w}}_s^e = {\boldsymbol{P}}_s^e. $

式中：$ {\boldsymbol{K}}_s^e $为滑梁单元的刚度矩阵. 由式(18)和式(19)可知，式(20)是关于${\boldsymbol{w}}_s^e $的非线性方程，$ {\boldsymbol{K}}_s^e $的值与挠度${\boldsymbol{w}}_s^e $及梁轴向力F_ξ的分布有关. 图7中滑梁只有两端承受轴向力F_ξ，因此滑梁的轴向内力应当处处相等.

对式(20)求微分得到：

(21)$ {\boldsymbol{K}}_{{\text{T}}s}^e\Delta {\boldsymbol{w}}_s^e = \Delta {\boldsymbol{P}}_s^e. $

式中：$ {\boldsymbol{K}}_{{\text{T}}s}^e $为单元切线刚度矩阵，且

$ {\boldsymbol{K}}_{{\text{T}}s}^e = {\boldsymbol{K}}_{{s}}^e+{\boldsymbol{K}}_{{\text{T}}{s_1}}^e+{\boldsymbol{K}}_{{\text{T}}{s_2}}^e+{\boldsymbol{K}}_{{\text{T}}{s_3}}^e+{\boldsymbol{K}}_{{\text{T}}{s_4}}^e $，其中

$ {\boldsymbol{K}}_{{\text{T}}{s_1}}^e = \int {2E{I_{\text{z}}}w'{A^2}{{\boldsymbol{C}}_2}{\boldsymbol{w}}_s^e{\boldsymbol{N}}_s^{e \prime}{\mathrm{d}}s} $，

$ {\boldsymbol{K}}_{{\text{T}}{s_2}}^e = \int {4E{I_{\text{z}}}w'w''{A^3}{{\boldsymbol{C}}_1}{\boldsymbol{w}}_s^e{\boldsymbol{N}}_s^{e \prime}{\mathrm{d}}s} $，

$ {\boldsymbol{K}}_{{\text{T}}{s_3}}^e = \int {2E{I_{\text{z}}}w''{A^2}{{\boldsymbol{C}}_1}{\boldsymbol{w}}_s^e{\boldsymbol{N}}_s^{e \prime \prime}{\mathrm{d}}s} $，

$ {\boldsymbol{K}}_{{\text{T}}{s_4}}^e = \int {F_{{\xi _2}}^ew'A{{\boldsymbol{C}}_1}{\boldsymbol{w}}_s^e{\boldsymbol{N}}_s^{e \prime}{\mathrm{d}}s} $.

对每根滑梁中所有的$ {\boldsymbol{w}}_s^e $、$ {\boldsymbol{K}}_s^e $和$ {\boldsymbol{P}}_s^e $以及$ \Delta {\boldsymbol{w}}_s^e $、$ \Delta {\boldsymbol{K}}_s^e $和$ \Delta {\boldsymbol{P}}_s^e $进行集成后，得到滑梁L的位移和载荷阵列及其微分$ {\boldsymbol{w}}_s^L $、$ {\boldsymbol{P}}_s^L $、$ \Delta {\boldsymbol{w}}_s^L $和$ \Delta {\boldsymbol{P}}_s^L $，以及滑梁L的刚度矩阵$ {\boldsymbol{K}}_s^L $和切线刚度矩阵$ {\boldsymbol{K}}_{{\text{T}}s}^L $，并满足

(22)$ {\boldsymbol{K}}_s^L{\boldsymbol{w}}_s^L = {\boldsymbol{P}}_s^L, $

(23)$ {\boldsymbol{K}}_{{\text{T}}s}^L\Delta {\boldsymbol{w}}_s^L = \Delta {\boldsymbol{P}}_s^L. $

在$ \Delta {\boldsymbol{w}}_s^L $与$ \Delta {\boldsymbol{P}}_s^L $中提取Q点对应的挠度项$ \Delta w_{{Q}}^L $与载荷项$ \Delta F_{{Q\text{z}}}^L $，则式(23)分解为

(24)$ \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{K}}_{{{{\mathrm{T}}Q}}}^L}&{{\boldsymbol{K}}_{{\text{TN}}{{{Q}}_1}}^L} \\ {{\boldsymbol{K}}_{{\text{TN}}{{{Q}}_{\text{3}}}}^L}&{{\boldsymbol{K}}_{{\text{TN}}{{{Q}}_{\text{4}}}}^L} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\Delta w_{{Q}}^L} \\ {\Delta {\boldsymbol{w}}_{{\text{N}Q}}^L} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\Delta F_{{Q \text{z}}}^L} \\ {\Delta {\boldsymbol{P}}_{{\text{N}Q}}^L} \end{array}} \right]. $

滑梁L在Q点处的支撑刚度k_L为

(25)$ {k_L} = {{\partial F_{{Q\text{z}}}^L} \mathord{\left/ {\vphantom {{\partial F_{{Q\text{z}}}^L} {\partial w_{{Q}}^L}}} \right. } {\partial w_{{Q}}^L}} = K_{{\text{T}Q}}^L. $

计算所有滑梁对Q点的支撑刚度并进行集成，得到滑梁的支撑刚度矩阵$ {{\boldsymbol{K}}_{\text{h}}} = {\mathrm{diag}}\left( {{k_1},{k_2}, \cdots ,}\right. $$\left. {k_L} ,\cdots \right) $.

对K_b、K_l与K_h进行集成，得到底箔的刚度矩阵K_B. K_B与所有滑梁在Q点的挠度阵列W有关. 需要先设定滑梁在Q点的初始挠度阵列W⁰，再结合图7中轴套对滑梁端部的接触约束，采用直接迭代法对式(22)确定滑梁端部载荷F_ξ，然后由式(23)和式(25)确定每个滑梁中Q点的支撑刚度k_L，并计算$ {\boldsymbol{K}}_{\text{B}}^1 $. 由$ {\boldsymbol{K}}_{\text{B}}^1 $计算出下一轮滑梁Q点的挠度阵列W¹，进行迭代. 迭代的收敛判据为

(26)$ {\varepsilon }_{w}^{t}={\Vert {{\boldsymbol{W}}}^{t}-{{\boldsymbol{W}}}^{t-1}\Vert }_{\infty }/{\Vert {{\boldsymbol{W}}}^{t}\Vert }_{\infty }\leqslant {10}^{-6}. $

式中：t为迭代层数，$ \varepsilon _w^t $为挠度的相对误差，在收敛后获得精确的K_B. 当给定气膜压力时，结合顶箔刚度矩阵K_top，采用拉格朗日法处理顶箔与底箔间的接触，得到所有箔片在特定气膜压力下的形变.

对于稳定运转的滑梁式空气轴承，可以看作等温状态，采用雷诺润滑方程：

(27)$ \frac{\partial }{{\partial \psi }}\left( {p{h^3}\frac{{\partial p}}{{\partial \psi }}} \right)+\frac{\partial }{{\partial \lambda }}\left( {p{h^3}\frac{{\partial p}}{{\partial \lambda }}} \right) = 6{\eta _{\text{A}}}\omega R\frac{{\partial \left( {ph} \right)}}{{\partial \psi }}. $

式中：ψ和λ分别为轴承的周向角度和轴向坐标，h和p分别为气膜的厚度和压力分布函数，η_A为空气的动力黏度.

当已知顶箔的挠度分布函数w时，气膜厚度的表达式为

(28)$ h = {h_0}+{e_0}\cos \left( {\psi - {\psi _0}} \right)+w. $

式中：取燕尾槽中心处的ψ=0；h₀为标准气膜厚度，即e₀=0时的半径间隙.

式(27)为椭圆方程，采用有限差分法求解. 在圆周方向上采用连续性条件；在轴向端面处，气膜压力满足边界条件：$ p = {P_0} $ （P₀为环境压强）.

设定转子偏心距e₀和偏位角$\psi _0$，给定初始的顶箔位移阵列$ {\boldsymbol{\alpha }}^0_{ \text{top}} $，并提取顶箔挠度w⁰，由式(28)得到气膜厚度h₀，再由式(27)求解气膜压力分布p⁰；通过p⁰生成顶箔的载荷阵列$ {\boldsymbol{P}}^0_{{\mathrm{top}} }$，并求解出新的顶箔挠度分布$ {\boldsymbol{\alpha }}^1_{{\mathrm{top}} } $，进行迭代，收敛判据为

(29)$ {\varepsilon }_{{\boldsymbol{\alpha}} }^{r}={\Vert {{\boldsymbol{\alpha}} }^r_{{\mathrm{top}}}-{{\boldsymbol{\alpha}} }^{r-1 }_{{\mathrm{top}}}\Vert }_{\infty }/{\Vert {{\boldsymbol{\alpha}} }^r_{{\mathrm{top}}}\Vert }_{\infty }\leqslant {10}^{-6}. $

式中：r为迭代层数，$ \varepsilon _{\boldsymbol{\alpha}} ^r $为迭代的相对误差.

迭代收敛后，得到给定$e_0$和$\psi _0$下的气膜压力分布p^r. 对p^r进行周向和轴向积分后，得到气膜压力的合力，其大小以及与竖直方向的夹角分别为$ F_e^0 $和$ \psi _e^0 $. 根据文献[18]，在得到气膜压力与厚度的分布后，可以进一步计算气膜内的速度场，并推导出轴承摩擦力矩M_f和端泄量Q_R的计算公式：

(30)$ {M_{\text{f}}} = \iint {\left( {\frac{{R\omega \eta _A }}{h}+\frac{h}{2}\frac{{\partial {p^r}}}{{R\partial \psi }}} \right)}{\mathrm{d}}\psi {\mathrm{d}}\lambda , $

(31)$ {Q_{\text{R}}} = \int {{\rho _{\text{A}}}\left( {\frac{{R\lambda h}}{2} - \frac{{{h^3}}}{{12\eta }}\frac{{\partial {p^r}}}{{R\partial \psi }}} \right)} {\mathrm{d}}\lambda . $

式中：M_f对整个气膜区域进行积分，Q_R对整个圆周方向进行积分，ρ_A为空气密度.

在实际试验中，载荷的大小是给定的，且一般沿竖直方向. 因此，须通过迭代修正的方式计算给定载荷$ {F_{\text{C}}} $下的偏心距和偏位角，收敛判据为

(32)$ \left.\begin{split} &{\varepsilon }_{F}^{q}={| {F}_{e}^{q}-{F}_{\text{C}}| }\leqslant {10}^{-3}\;\text{N}，\\ &{\varepsilon }_{\psi }^{q}={| {\psi }_{e}^{q}| }\leqslant {10}^{-3}\;\text{rad}.\end{split}\right\} $

式中：q为迭代层数，$ \varepsilon _F^q $和$ \varepsilon _\psi ^q $分别为$ F_e^q $和$ \psi _e^q $的绝对误差.

综上，当给定结构参数及载荷F_C时，箔片轴承的运行参数及润滑性能的数值求解方法如图9所示，包括箔片位移求解模块、流固耦合模块和偏心距与偏位角求解模块.

通过MATLAB编程实现图9中的计算流程. 为了检验所建模型的准确性，搭建试验台进行验证.

试验轴承、试验台及试验原理如图10所示，试验轴承的参数如表1所示. 表中，t₁为底箔厚度，t₀为顶箔厚度.

转子与伺服电机的轴通过联轴器直接相连. 箔片轴承套装在转子上，轴承上端通过细线与砝码相连. 砝码对轴承施加向上的拉力，根据力的相互作用，等价于转子对轴承施加向下的力. 轴套与力臂杆相连，力臂杆的一端以自由的方式压在压力传感器上. 当转子转速提高到一定程度时，形成流体动压润滑气膜. 随着转子的旋转，气膜受到剪切作用，使轴承套产生旋转的力矩，带动力臂杆给传感器施加压力. 当转子旋转时，通过传感器测得压力值F_SEN，可以进一步计算轴承的摩擦力矩M_EP：

(33)$ {M_{{\text{EP}}}} = {F_{{\text{SEN}}}}\left( {{l_{{\text{pole}}}}+{R_{{\text{SL}}}}} \right). $

式中：l_pole为力臂杆的长度，R_SL为轴套外表面半径.

当转速较低时，轴承未起飞，转子与箔片之间存在接触摩擦. 取转速n_r＞1×10⁴ r/min的试验数据与数值结果对比，如图11所示. 可以看出理论模拟结果与试验结果在数值及变化趋势方面均具有良好的一致性，两者吻合度较高，说明本研究建立的理论模型具有良好的可靠性.

对轴承进行数值研究，轴承结构参数如表1所示. 已知大气压强P₀=101 325 Pa，流体动力学黏度η_A=$ 1.932 \times {10^{ - 5}}{\text{ Pa}} \cdot {\text{s}} $，弹性模量E=210 GPa，泊松比ν=0.3. 当转速n_r=$ 4 \times 10^4{\text{ r}}/\min $、载荷F_C=8 N时，无量纲气膜压力$ \overline P $（$ \overline P = {p \mathord{\left/ {\vphantom {p {{P_0}}}} \right. } {{P_0}}} $）和气膜厚度h的分布如图12所示.

图12(b)中，轴向中面处（λ=0）的气膜厚度大于两端（$ \lambda = \pm 12.5{\text{ mm}} $），表明顶箔在中面上的形变较大，这与图12(a)中气膜压力在轴向的分布规律相对应. 气膜厚度的分布并不均匀，呈现局部不平整现象，这是因为滑梁在底箔上呈间断分布，导致其对顶箔的支撑不连续，从而使顶箔在形变后出现凹凸不平的情况. 在图12(a)中，局部不平整的气膜厚度使气膜压力的分布呈现出局部不平整现象，并且在压力越大的区域这种不平整现象越明显.

调整标准气膜厚度h₀，保持其他参数不变，计算结果如图13所示. 图中13 (a)、(b)为λ=0处的气膜压力及气膜厚度. 随着h₀的增大，气膜厚度显著增大，气膜的最大压力上升，动压范围减小. 图13(c)中，当 h₀由10 μm增大到25 μm时，摩擦力矩降低了约1 N·mm，这是因为气膜厚度的增大使得气膜中流场的速度梯度减小，从而降低了空气内摩擦力. 由于h₀反映轴承与转子间的初始间隙，影响气流的出口大小，h₀越大，轴承的端泄量就越大，这与图中的结果一致.

滑梁的尺寸参数受到滑梁数目与梯形形状的限制，长、宽过大将会导致结构之间的贯穿，斜率过大将不能维持滑梁的梯形形状. 研究滑梁尺寸对轴承气膜压力与厚度分布的影响，采用控制变量法，分别取$ {L_0} \in \left[ {1.50{\text{ , 2}}{\text{.40}}} \right]{\text{mm}} $ 、$ d \in \left[ {1.20{\text{ , 1}}{\text{.60}}} \right]{\text{mm}} $和$ \gamma \in [ 0.01, {\text{0.19}} ] $，计算结果如图14所示. 随着滑梁长度的增大，转子的偏心距逐渐增大，而偏位角逐渐减小. 这是因为滑梁长度的增大会降低滑梁的支撑刚度，从而降低整个轴承对转子的支撑刚度. 在相同的载荷下，转子在竖直方向会产生更大的位移. 图14(d)中，气膜厚度的增大主要是因为滑梁刚度的降低导致顶箔发生了更大的形变.

滑梁宽度的减小和斜率的增大也会降低滑梁的支撑刚度. 对比图14(d)~(f)可知，滑梁长度对气膜厚度的影响最大，而滑梁斜率的影响最小. 图14(a)~(c)中，滑梁长度的增大使气膜的最大压力降低，局部不平整现象减弱，而滑梁宽度和斜率的变化对气膜压力分布的影响很小.

相比于h₀，滑梁尺寸的变化对轴承摩擦力矩和端泄量的影响较小，因为与图13(b)相比，图14(d)~(f)中气膜厚度的变化总是在较小范围内，无法对气体的内摩擦力产生更大的影响.

分别调整顶箔和底箔的厚度t₀和t₁，保持其他参数不变，计算结果如图15所示. 其中，底箔厚度的变化同样会改变滑梁的支撑刚度，因此，图15(b)和15(d)中，随着底箔厚度增大，气膜的最大压力升高，气膜的厚度减小. 顶箔厚度的变化对气膜压力和气膜厚度的整体分布趋势的影响很小，但是随着顶箔厚度的降低，气膜压力和气膜厚度的局部不平整性更加明显. 当顶箔厚度从0.20 mm降低到0.10 mm时，在最小气膜厚度附近，气膜厚度的局部波动值从0.1 μm增大到 0.4 μm.

（1） 基于Kirchhoff理论建立顶箔的有限元模型，并对底箔结构进行简化处理. 在模型中，考虑底箔连接梁之间的位移耦合，根据梁的大变形理论计算滑梁的等效刚度. 进一步建立底箔的非线性刚度矩阵，并计算箔片的形变. 采用等温雷诺方程计算气膜压力分布，通过流固耦合方法构建滑梁式空气箔片轴承的力学模型，计算在特定载荷下轴承的性能参数. 该模型的数值结果与试验结果良好地吻合.

（2） 结合具体的算例，研究箔片的结构参数对轴承性能的影响. 结果表明，随着轴承与转子初始间隙的增大，气膜动压区域减小，气膜的最大压力上升，摩擦力矩减小，端泄量增大；滑梁尺寸和底箔厚度对气膜的影响是通过改变滑梁支撑刚度实现的，当滑梁对顶箔的支撑刚度增大时，气膜厚度减小，气膜的最大压力增大. 由于滑梁对顶箔的支撑不连续，气膜厚度与气膜压力呈现局部不平整现象，并且随着压强的增大与顶箔厚度的减小，局部不平整现象加剧.

（3） 通过数值计算方法研究滑梁结构参数对滑梁箔片轴承润滑性能的影响，得出一系列规律和结论. 然而，影响箔片轴承润滑性能的因素还有很多，例如表面粗糙度、温度、振动、形位公差等. 在今后的研究中还需要进一步地把这些因素纳入考虑范围，进行更加全面和深入的研究，以更准确地评估和优化滑梁箔片轴承的润滑性能.