气体箔片轴承静态工作点求解

邓志凯，程文杰，曹广东，肖玲，李明

(西安科技大学 理学院，西安 710054)

气体箔片轴承(Gas Foil Bearings，GFB)作为一种无油轴承，具有摩擦功耗低，工作温度范围宽，耐冲击和装配对中要求低等优点。GFB已经成功应用于极端环境(dmn值超过3.0×106mm·r/min和538 ℃的服役环境)的气体透平机中， 并且在燃料电池、发电机、飞机推进器、气体处理领域有广泛的应用前景[1-3]。然而，GFB的设计需考虑多种工况，同时涉及多学科交叉的知识，其静动态性能及性能极限的计算往往基于刚性表面的假定，该假定限制了流体动力润滑理论在一定工况下的应用。

国内外诸多学者致力于建立更完善的模型来计算轴承的性能，以减少刚性表面轴承的假定对计算结果造成的误差。文献[4]在略去波箔之间的相互作用及库伦摩擦的基础上，采用柔度系数描述底层支承拱箔弹性变形对气膜厚度的影响，为降低气弹耦合计算的复杂性，对顶层箔片的变形未深入考虑。由于箔片结构在气膜压力作用下产生弯曲变形，轴承间隙随之改变，气膜厚度也相应增大，因此，气膜与箔片结构之间存在强耦合作用[5]。文献[6]采用二维Navier-Stokes方程求解流场效应，结合箔片结构的有限元法对流体-箔片-轴颈之间的相互作用进行了稳态、准瞬态和瞬态动力学分析，在分析中得到了耦合轴颈运动、流体、箔片结构的动态性能，并与刚性表面轴承的动态性能进行了比较。文献[7]通过在有限元公式中耦合气体流动和箔片的变形来分析轴承的静态性能，研究发现气膜出口区的膜厚变化是由于顶层箔片的脱离所引起。文献[8-9]采用有限差分法以更简单的方式研究了类似的流固耦合分析。诸多研究表明,对气体箔片轴承的完整分析需要考虑库伦摩擦的影响。文献[10-11]的研究指出拱箔与顶层箔片之间的微滑动产生了摩擦阻尼，该研究中提出的半经验模型分析了单个拱箔的摩擦阻尼及其与相邻拱箔的相互作用。

以上对气体箔片轴承的研究考虑了箔片结构与气膜之间的耦合作用，使轴承静动态性能更接近实际值。然而，关于轴承静态工作点的确定却鲜有文献提及，已有文献中给出的轴承性能多是基于给定偏心率(直接偏心法)来计算轴承的气体流量、摩擦损耗和承载力。为确定任意给定载荷下气体箔片轴承的静态工作点，本文提出了一种“二分法搜索+不动点迭代”的求解策略，其中，二分法用于搜索偏心率，不动点迭代用于寻找偏位角。弹流耦合求解中雷诺方程采用基于质量流量守恒的差分法离散进行超松弛迭代求解，顶层箔片采用Kirchhoff薄板模型进行有限元求解，以整周式GFB及三瓦插入式GFB为算例进行求解及分析。

1 气体箔片轴承模型及假定

1.1 气体箔片轴承模型

本文重点研究弹流耦合过程对轴承静态工作点的影响，由于顶层箔片变形会影响气膜厚度分布，为方便计算，采用单一的厚箔片替代完整的箔片结构。从轴承承载机理上分析，箔片结构的总刚度矩阵K=Kt+Kb，其中Kt为顶层箔片刚度矩阵，Kb为波箔刚度矩阵，当采用厚箔片时，顶层箔片刚度矩阵可近似替代箔片结构的总刚度矩阵。GFB的结构形式如图1所示，图中O，O′分别为轴承中心与轴颈中心，Ω为转子工作转速，对于三瓦插入式GFB，定义α为每块瓦的瓦张角，β为瓦位角(由轴承上方垂线到第1块轴瓦进气边的位置角)，ξ为槽宽包角。

1.2 计算假定

在确定轴承静态工作点的过程中，为减少耦合分析中每次迭代的时长，求解静态气膜压力场及顶层箔片变形时基于以下假定：

1)气体润滑过程是恒温的；

2)流场是定常的，轴承端部及开槽处的气膜压力为标准大气压；

3)不考虑顶层箔片的曲率效应。

2 静态气膜压力的求解

2.1 求解理论

以整周式GFB为例，无论是静态气膜压力还是动态气膜压力的求解，气膜厚度h都是气膜压力的决定性参数，其表达式为

h=C0+ecos(φ-θ)，

(1)

式中：C0为轴承间隙；e为偏心距；φ为轴承展开角；θ为偏位角。

在此基础上，静态气膜压力的求解可归结为定常二维雷诺方程[12]的求解

(2)

式中：p为气膜压力；x,y,z分别为轴承水平方向、竖直方向、轴向方向的坐标分量；μ为气体动力黏度；U为转子表面沿x方向上的速度分量。椭圆型偏微分方程(2)式无法进行解析求解，大多数学者直接从雷诺方程出发，采用有限差分法构造每一项的差分格式后代入原方程迭代求解。本文基于质量流量守恒原理[2,12]建立包含所有节点压力在内的差分方程，该方法在处理雷诺方程并不成立的特殊边界节点时依然具有统一的表达形式。

将轴承沿周向展开并对周向和轴向进行网格划分，如图2所示，任意节点(i,j)区域所满足的守恒关系为

图2 任意节点处的质量流量守恒关系Fig.2 Conservation of mass flow at any node

Qa+Qb-Qc-Qd-Qe-Qf+Qg+Qh=0，

(3)

式中：Q为质量流量。

根据文献[12]的推导，单位时间通过任一单位宽度周向和轴向截面的质量流量为

(4)

式中：P为量纲一的气膜压力；H为量纲一的气膜厚度；Λ为轴承数；ψ，λ分别为轴承量纲一的周向长度和轴向长度；Patm为标准大气压；R为轴承半径。

b～h点的推导与a点类似, 以a点为例，其质量流量、气膜厚度、气膜压力以及偏导数的半步长差分格式为

(5)

式中：Δψ，Δλ分别为轴承周向、轴向微段长度；i，j分别为轴承周向、轴向网格节点；P0为气膜压力的迭代初始预设值；Pa0为节点a处的气膜压力初始值。

将(5)式代入(4)式，即可得到任意节点(i,j)处的气膜压力差分方程。由于轴承端部及轴瓦进出口压力取环境压力，(4)式的边界条件为

P(ψ,0)=P(ψ,L)=1，

(6)

式中：L为轴承长度。

对已构造差分格式的所有节点，给出边界条件后可采用超松弛迭代法进行差分方程组的求解。

气膜力分量Fx，Fy及气膜承载力W的计算公式为

(7)

以上关于整周式GFB的求解过程同样适用于三瓦插入式GFB静态气膜压力的计算，需要注意的是，后者在除轴承端部之外的开槽处也取环境压力，其边界条件在满足(6)式的基础上还要满足

P(ψm,λ)=1;m=1,2,3，

(8)

此外，三瓦插入式GFB的W满足叠加原则,即

(9)

2.2 GFB气膜压力求解算例

为验证气膜压力计算程序的可靠性，以整周式刚性表面轴承为算例，首先给出本文计算与文献[2]计算的气膜压力分布对比(图3)，其中，轴承结构参数及其他性能参数见表1。仿真结果显示，本文计算得到的量纲一的承载力为3.73，与文献[2]的计算结果差值在5%以内。

图3 整周式GFB静态气膜压力分布对比Fig.3 Comparison of static gas film pressure distribution of single-pad GFB

表1 整周式GFB结构参数及其他性能参数Tab.1 Structure parameters of single-pad GFB and other performance parameters

为比较三瓦插入式GFB与整周式GFB的气膜压力分布特点，刚性表面假设下整周式GFB及三瓦插入式GFB的静态气膜压力分布分别如图4、图5所示，其中，三瓦插入式GFB的瓦张角α=114°，瓦位角β=38°，槽宽包角ξ=18°，其余参数与表1相同。与整周式GFB相比，三瓦插入式GFB在轴向方向上的气膜压力分布更均匀。2种GFB在轴承中截面(λ=1.5)上的气膜厚度和气膜压力分布分别如图6、图7所示。不难发现，三瓦插入式GFB中轴向槽的引入使气膜厚度沿轴承周向分布不再是连续函数，同时，在轴承周向方向上出现3个气膜压力峰，与整周式GFB相比，三瓦插入式GFB气膜压力沿轴承周向分布较为均匀，但气膜承载力有所下降。

图4 整周式GFB静态气膜压力分布Fig.4 Static gas film pressure distribution of single-pad GFB

图5 三瓦插入式GFB静态气膜压力分布Fig.5 Static gas film pressure distribution of three-pad insertion GFB

图6 GFB气膜厚度分布Fig.6 Gas film thickness distribution of GFB

图7 GFB气膜压力分布Fig.7 Gas film pressure distribution of GFB

3 顶层箔片的弯曲变形

3.1 Kirchhoff板理论

在GFB的设计中，顶层箔片的结构形状和宏观变形是箔片轴承静动态性能的主导性因素，就无预紧箔片轴承而言，顶层箔片所具有的间隙接近于名义气隙的平均值[2]，因此，本文针对无预紧轴承进行研究。由于顶层箔片的厚度相比其长、宽尺寸小了几个数量级，顶层箔片的建模可以归结为横向复杂载荷下的薄板弯曲。

采用有限差分法计算薄板的弯曲变形存在收敛条件较高，复杂载荷下难以收敛的问题，相反，采用有限单元法则极大地缩减了求解时间，这对耦合求解气膜压力场及箔片变形来说无疑加快了每一步的迭代过程。因此，本文采用有限元法自主编程求解顶层箔片的弯曲变形。

在已有关于弹性板的弯曲理论中，Kirchhoff板理论更适用于顶层箔片的力学建模，由于忽略了板的横向剪切变形，板内所有的力学量都可用板中面的挠度ω(x,y)表示。顶层箔片简化后的Kirchhoff薄板模型如图8所示，该薄板模型中左右边界分别采用固支和简支。板中每个节点有3个广义位移分量(挠度ω，绕x轴的转角θx和绕y轴的转角θy)。在此基础上，由四节点矩形单元位移模式下的形函数得到薄板的单元刚度矩阵

图8 Kirchhoff薄板模型Fig.8 Model of thin plate in Kirchhoff

(10)

式中：B为几何矩阵；D为弹性矩阵；N为形函数，N=[N1N2N3N4]，其子矩阵N1～N4推导较为复杂，具体表达式可参考文献[13]。

将各单元的刚度矩阵进行组装，可建立薄板结构的整体刚度矩阵。若薄板单元受横向气膜力的作用，则等效结点力列阵为

(11)

式中：Fzi，Mθxi，Mθyi分别为作用在节点i上的节点力、绕x轴和绕y轴的节点弯矩。各节点的位移列阵为

(12)

式中：ωi为薄板节点的挠度；θxi，θyi分别为节点绕x轴和绕y轴的转角。

3.2 箔片变形求解算例

图9—图11给出了轴承数Λ=1，长径比L/D=1.5，偏心率ε=0.7，偏位角θ=38.08°下2种GFB的顶层箔片变形。其中，三瓦插入式瓦张角α=114°,瓦位角β=38°，箔片结构参数见表2。由图9、图10可知，尽管整周式GFB的箔片厚度是三瓦插入式GFB的2倍还多，但其径向变形依然比后者大出一个数量级,原因为三瓦插入式GFB的多槽结构使气膜沿周向分布相对均匀且总承载力有所降低,每块瓦(顶层箔片)的长宽比更接近于1，若将图8所示薄板模型简化为二维欧拉梁模型，在横截面不变的情况下，顶层箔片变形将与其板长呈三次非线性关系。因此，在箔片参数相同的情况下，相对整周式GFB，三瓦插入式GFB每块瓦的周向长度更短，顶层箔片变形更小，其变形后的承载轮廓面更接近于刚性轴承。由于气膜压力分布的均匀化，三瓦插入式GFB更有利于转子的旋转稳定性。在三瓦插入式GFB中，当某块箔片位于膜厚最小区域时，其径向变形约是其余箔片的5倍。

图9 整周式GFB气膜压力下顶层箔片的变形Fig.9 Deformation of top foil under gas film pressure for single-pad GFB

图10 三瓦插入式GFB顶层箔片的径向变形Fig.10 Radial deformation of top foil under gas film pressure for three-pad insertion GFB

图11 三瓦插入式GFB轴向方向膜厚最小处顶层箔片的径向变形Fig.11 Radial deformation of top foil at minimum gas film thickness in axial direction of three-pad insertion GFB

表2 顶层箔片的结构参数Tab.2 Structural parameters of top foil

由图11可知，位于轴承端部的顶层箔片径向变形较小，与文献[2,14]得出的结果相似，这种变形特征减少了气体在轴承端部的泄漏。

4 静态工作点求解

4.1 二分法搜索确定偏心率

对于任意给定的轴承参数及外载荷，偏心率对气膜承载力的影响远大于偏位角，因此，先采用二分法确定偏心率，在此之前，需预设初始偏位角θ0。

对于给定的偏心率搜索区间(ε1,ε2)，为判断预设区间的合理性，需要分别计算(ε1,θ0)，(ε2,θ0)下的气膜承载力与外载荷的差值W1，W2。若W1W2>0，则所求的偏心率不在搜索区间范围内，反之，则令εx=(ε1+ε2)/2，开始在所给区间内由(εx,θ0)计算气膜承载力与外载荷的差值Wx。为保证解的精度，需给出偏心率的允许误差Δ，总迭代次数Nt要满足

(13)

在搜索过程中，当W1Wx>0时，则令εx=ε1，反之，令εx=ε2，直至迭代次数n>Nt时迭代终止，此时，ε即为要找的偏心率。

4.2 不动点迭代寻找偏位角

对于任意给定的外载荷角αf，静态偏位角与αf差值δ不能过大，需要满足收敛条件

(14)

式中：Fx，Fy分别为(ε,θn)下的气膜力在水平与竖直方向的分量；γ为偏位角的允许误差，一般取γ<1.0×10-4。预设偏位角后，新的偏位角可按照如下迭代格式进行寻找

(15)

当δ<γ时，迭代终止。至此，偏心率、偏位角都已确定。

4.3 求解策略

在应用“二分法搜索+不动点迭代”求解轴承静态工作点的过程中，气膜压力采用有限差分求解，薄板变形运用有限元法求解，其中网格大小相同，单元及结点也完全相同，于是，相同结点的气膜力及薄板挠度的传递在程序上容易实现。此外，在迭代过程中，需要不断用箔片变形更新气膜厚度，具体求解步骤如下：

1)预设初始偏位角，采用二分法迭代寻找偏心率，给定偏心率的迭代区间；

2)由当前迭代步中的偏心率、偏位角计算静态气膜压力场中的气膜厚度h1，对压力场积分获得气膜合力F1，将气膜压力场代入箔片变形场中，计算得到箔片变形ω1；

3)将第2步中得到的箔片变形ω1代入静态气膜压力场中，修正气膜厚度得到新膜厚h2，重新计算气膜压力场并积分获得气膜合力F2，此时将气膜压力场带入箔片变形场得到箔片变形ω2；

4)重复第2步和第3步直到气膜合力Fn收敛；

5)判断气膜合力Fn与外载荷是否满足偏心率的收敛条件(满足二分法根的允许误差的最大迭代次数)，若满足，进入偏位角的迭代计算，若不满足，重复第2步到第5步；

6)在确定偏心率的基础上，采用不动点迭代寻找偏位角，由当前迭代步中的偏位角和已确定的偏心率计算气膜压力场及顶层箔片变形，迭代过程与第2步至第4步一致；

7)当外载荷角αf与偏位角θn满足偏位角的收敛条件(|αf-θn|≤10-4)时迭代终止，若不满足，重复上一步的迭代。

4.4 GFB静态工作点的求解算例

由于三瓦插入式GFB的求解流程与整周式GFB完全相同，因此，本文以具有代表性的整周式GFB为算例进行静态工作点的求解。表3给出了整周式刚性表面GFB在给定量纲一外载荷Wf下的静态偏心率和偏位角，其中Λ=0.6，L/D=1.0。

由表3可知：本文静态工作点求解方法得到的静态偏心率与文献[2]的相对误差最大为6.60%，随着偏心率的增大，相对误差显著减小，且小偏心率的相对误差并不会造成过大的轴承承载力误差；偏位角与文献[2]的相对误差都小于3.80%。上述计算结果表明本文所采用方法的可靠性。需要说明的是，文献[2]的承载力是给出转子偏心率与偏位角后计算得到的，本文则是根据文献计算的承载力来反向搜索轴承的静态偏心率与偏位角。实际上，更关心某个给定质量的转子运行稳定时的工作位置，按本文的算法，任意给定一个转子及外载荷，都可以很快确定它的静态工作点，这也是本文求解静态工作点的目的。

表3 整周式刚性表面GFB静态偏心率与偏位角Tab.3 Static eccentricity and attitude angle of single-pad rigid surface GFB

相同载荷下刚性表面轴承与柔性表面轴承的静态偏心率和偏位角分别如图12—图14所示，其中轴承的结构参数见表4。

图12 刚性表面与柔性表面下轴承偏心率对比Fig.12 Comparison of eccentricity between rigid surface bearing and flexible surface bearing

图14 极坐标系下刚性、柔性表面轴承偏位角、偏心率的对比Fig.14 Comparison of attitude angle and eccentricity in polar coordinate system between rigid surface bearing and flexible surface bearing

表4 轴承的结构参数Tab.4 Structural parameters of bearing

由图12可知，随着外载荷的增加，柔性表面轴承相对刚性表面轴承的静态偏心率不断增大，在载荷最大处其偏心率差值约为0.05，在重载情况下，即使偏心率差值较小，也会造成较大的承载力误差。实际上，若箔片表现出更软的特性，其变形越大,刚性表面轴承的静态偏心率与实际差值也将越大。比较图12、图13可知，无论是刚性表面轴承还是柔性表面轴承，相同载荷条件下，偏位角随着偏心率的增加都会逐渐减小，不同的是，随着偏心率的增加，柔性表面轴承的偏位角递减的更快一些。

图13 刚性表面与柔性表面下轴承偏位角对比Fig.13 Comparison of attitude angle between rigid surface bearing and flexible surface bearing

图14给出了更直观的静态工作点“轨迹图”，不难发现，柔性表面轴承的静态工作点轨迹表现出了更明显的下沉，该现象与箔片变形产生的偏转有关。需要注意的是，(2)式等式右端的项表示由气膜旋转效应所贡献的气膜承载力，本文采用的定常流体润滑假设忽略了气膜挤压效应下的气膜承载力，若计入这种效应的影响，箔片将会产生更大的变形，此时，刚性表面轴承与柔性表面轴承的偏位角差值将会更大。因此，在GFB的动态性能计算中，更完善的模型应考虑雷诺方程中关于时间的微分项。

图15—图17分别给出了Wf=0.75,ε=0.62，θ=61.90°时轴承间隙最小处的量纲一的膜厚、量纲一的轴承承载力、箔片变形随迭代次数的变化曲线。其中，首次迭代计算出的气膜厚度为刚性表面轴承的膜厚h1；二次迭代时则引入首次计算的气膜合力F1下的箔片变形ω1，由图可知，最小膜厚经过箔片变形的修正会显著变大为h2(h2>h1)，此时气膜合力减小为F2(F2h3>h1)；重复以上过程，膜厚和承载力最终收敛于定值。实际上，上述流体-箔片的耦合分析中，也可理解为气膜和箔片同时经历了振动的瞬态和稳态过程。

图15 轴承量纲一的气膜厚度随迭代次数的变化Fig.15 Variation of dimensionless gas film thickness for bearing with iterations

图16 轴承量纲一的承载力随迭代次数的变化Fig.16 Variation of dimensionless load capacity of bearing with iterations

图17 顶层箔片量纲一的变形随迭代次数的变化Fig.17 Variation of dimensionless deformation for top foil with iterations

另外，以上计算结果显示，气体箔片轴承采用刚性表面假设会带来不小的误差，尤其是在大偏心率的情况下计算得到的结果很可能失真。由于本文算例中所采用的箔片尺寸相对较厚，箔片变形很小，在工程实际中箔片的厚度只有零点几毫米，往往会出现偏心率大于1的情况，因此，计入箔片的弯曲变形对气体箔片轴承的静、动态性能计算非常必要。

5 结论

通过“二分法搜索+不动点迭代”的求解策略寻找气体箔片轴承的静态工作点，得出以下结论：

1)采用本文设计的算法可确定任意给定外载荷下轴承的静态工作点，刚性表面轴承的偏心率与偏位角的计算结果与文献[2]的计算结果吻合较好(偏心率相对误差不超过6.60%，偏位角相对误差小于3.80%)。

2)无论是刚性表面轴承还是柔性表面轴承，相同载荷条件下，偏位角随着偏心率的增加都会逐渐减小，不同的是，相对于刚性表面轴承，柔性表面轴承偏位角递减的更快一些。

3)在计算气体箔片轴承时，采用刚性表面假设会带来不小的误差，尤其是在大偏心率的情况下计算得到的结果很可能失真。考虑到工程中往往会出现偏心率大于1的情形，计入箔片的弯曲变形来计算轴承的静动态性能更加合理。