基于Timoshenko梁单元的径向波箔轴承箔片变形分析

蔡鹏程，马希直，原兆祥

(南京航空航天大学 直升机传动技术国防科技重点实验室，江苏 南京 210016)

0 引言

最近30年来，以波箔轴承为代表的弹性箔片轴承越来越被关注和研究，波箔轴承的复杂箔片结构形式和受力变形方式以及轴承运行过程中的气热弹耦合问题使其理论模型的建立变得十分困难，很难准确预测轴承的静动态特性[1]。

1975年，WALOWIT J A等人最早建立了波箔轴承箔片变形模型，考虑了波箔片与轴承套之间的摩擦，将箔片变形与Reynolds方程耦合求解[2]。1983年，HESHMAT H等人忽略摩擦，将波箔结构视为线性弹簧、波箔刚度分布平均，引入箔片变形系数，给出计算箔片变形的Heshmat公式[3]。1999年，LORDANOFF I 考虑了波箔片与轴承套之间的摩擦，给出了固定端和自由端波纹的刚度计算公式[4]。不足的是，这些计算公式都没有考虑平箔片在相邻波纹之间的变形，并且忽略了相邻波纹之间的相互作用。1992年，HESHMAT H在高倍放大镜下观察了箔片结构的受力变形情况，发现波箔受压会产生滑移现象。因此将摩擦因素和波纹相互作用考虑后对之前的理论进行了修正[5]。

2007年，LE LEZ S等人提出了一种新的波箔轴承箔片结构的模型，把每个波纹简化成相互连接的线性弹簧，如图1所示。箔片结构多自由度模型通过线性弹簧的刚度充分考虑了平箔片与波箔片之间的摩擦、波箔片与轴承套之间的摩擦预计相邻波纹之间的相互作用[6]。

图1 箔片结构多自由度(NDOF)模型[6]

2009年，SAN ANDRÉS L和KIM T H采用一维梁单元和二维板壳单元计算平箔片的变形，计算了波箔轴承的最小气膜厚度和动态刚度阻尼系数，通过与实验结果对比后，指出一维梁模型能更准确地预测波箔轴承气膜分布[7]。

2013年，徐方程研究了平箔片剪切变形对径向波箔轴承的静特性的影响，采用二维厚板单元对平箔片进行建模，研究结果表明平箔片的剪切变形对轴向气膜厚度的变化量有显著影响[8]。

2019年，ARGHIR M等人基于箔片结构的多自由度模型，指出该模型存在未考虑平箔片和波箔片、波箔片和轴承套之间接触间隙的不足，基于接触理论建立了考虑摩擦和接触间隙的箔片结构模型[9]。

本文采用Timoshenko梁单元对箔片结构进行建模，对平箔片和波箔片之间的摩擦、波箔片和轴承套之间的摩擦建立库仑摩擦模型，求解箔片变形并与线性弹簧模型和多自由度模型进行比较，基于求解出的变形分析径向波箔轴承的静特性。

1 箔片结构建模

采用Timoshenko梁单元对平箔片和波箔片进行离散，变形示意图如图2、图3所示。求解每个单元的刚度矩阵在整体坐标系中组装成整体刚度矩阵Kg，求解Kg·U=F即可得到各节点的位移，从而得到箔片的变形。

图2 平箔片变形示意图

图3 波箔片变形示意图

假设波箔片梁单元模型每个发生滑移单元的滑移方向都是从固定端向自由端滑移，波纹之间的平段为刚性，即：长度不变，平箔片和波箔片之间的摩擦系数为η。波箔片和轴承套之间的摩擦系数为μ。以两个波纹的波箔片为例，波箔片梁单元库仑摩擦模型受力情况如图4所示。

图4 波箔片库仑摩擦模型

2 波箔片变形分析

对有10个波纹的波箔片进行分析，其中第1个波纹与固定端连接，第10个波纹右端自由，有关的箔片参数见表1。

表1 箔片参数表

2.1 不同载荷下各波箔片模型的变形

对波箔片施加4种不同形式的压力载荷，作用在波纹顶端，每种载荷的单位面积等效载荷都为2×105Pa，得到的各个波纹的变形如图5-图8所示。图中水平轴为从固定端到自由端的波纹序号，竖直轴表示波纹顶端的竖直方向的位移，正值表示随压力向下变形。

图5 均布载荷下波箔片的变形

图6 递减载荷下波箔片的变形

图7 递增载荷下波箔片的变形

图8 对称递增-递减载荷下波箔片的变形

从图5-图8可以看出，对于线性弹簧模型，由于不考虑波纹之间的相互作用，大大低估了箔片结构的刚度，除了固定端和自由端的两个波纹，其余波纹的变形都是Timoshenko梁模型和NDOF模型的两倍以上。由于Iordanoff公式考虑了波箔片和轴承套之间的摩擦，变形相对于Heshmat公式小一些。

此外，Timoshenko梁模型和NDOF模型的变形折线比较吻合，验证了所建立的Timoshenko梁模型的合理性。从图中可以看出，各个波纹变形与载荷的关系是非线性的，即使相同的载荷，靠近固定端的波纹变形也是明显小于靠近自由端的波纹。从整体上看，Timoshenko梁模型的变形相对于NDOF模型要大一些。4种载荷中，对称递增-递减载荷最符合径向波箔轴承润滑气膜压力的周向变化。

2.2 摩擦系数对波箔片变形的影响

图9为均布载荷下、基于Timoshenko梁建立的波箔片模型在不同摩擦系数下的变形。从图中可以看出，当η=μ=0时，变形折线与图5中的线性弹簧模型相似，只是变形量要小一半。可见波纹之间的相互作用要靠摩擦来实现，不考虑摩擦就相当于将波箔片简化为一个个波纹组成的简单结构，波纹的变形是随载荷线性变化的；即使不考虑摩擦，Timoshenko梁模型的刚度也达到线性弹簧模型刚度的两倍左右。由此可见，使用线性弹簧模型计算的箔片变形是有很大误差的。从图中还可以看出，随着摩擦系数的增加，两端不发生滑移的波纹数目会增加，发生滑移后波纹的变形会突然增加。只要有波纹发生滑移，此后直到自由端的波纹变形由载荷的大小主导，摩擦系数大小对其影响十分微小。这一点在自由端的波纹变形上体现得最明显。

图9 均布载荷下波箔片的变形

2.3 剪切刚度对波箔片变形的影响

为了研究剪切刚度对于波箔片变形的影响，分别计算了4种摩擦系数下Timoshenko梁和Euler梁建立的波箔片模型的各个波纹的变形，如图10所示。从图中可以明显看出，不考虑剪切变形的Euler梁模型的变形比Timoshenko梁模型要小，所以剪切刚度会使波箔片变形增加，不考虑剪切刚度会对波箔片刚度产生较高估计。

图10 剪切刚度对波箔片变形的影响

3 Reynolds方程求解与静特性分析

定常理想气体润滑无量纲Reynolds方程为

(1)

采用Newton-Raphson迭代法结合有限差分法求解式(1)，采用Timoshenko梁箔片模型计算箔片变形wd，根据式(2)更新气膜厚度，不断迭代直到满足收敛条件式(3)，根据求得的气膜压力积分求解承载力和偏位角等静特性参数。

(2)

(3)

计算的轴承参数见表2，箔片参数与表1相同。

表2 轴承、气体参数

图11为转速为30 000 r/min时各个箔片模型的轴承承载能力随偏心率变化的折线图。从图中可以看出，在偏心率<0.4，也就是气膜压力较小时，各模型箔片变形差异较小，气膜压力基本相同，所以承载能力基本相同；当偏心率超过0.5，更容易变形的线性弹簧模型(Heshmat公式、Iordanoff公式)承载能力要远小于其他模型。

图11 各箔片模型承载力随偏心率变化图

箔片变形量接近的NDOF模型、Timoshenko梁模型和Euler梁模型承载能力相近，没有考虑平箔片在波纹之间“凹陷”变形的NDOF模型承载能力要比梁模型大。考虑了剪切刚度的Timoshenko梁模型虽然计算的变形比Euler梁模型稍大，但承载能力要比Euler梁模型稍强。为了研究其中的原因，计算了一定偏心率下轴承中截面的无量纲气膜厚度，如图12所示。

图12 轴承中截面无量纲气膜厚度图

从图12可见，一定的偏心率下，Euler梁模型的无量纲气膜厚度比Timoshenko梁模型大，也就是箔片变形量小一些，但是它们的最小气膜厚度十分接近。所以Euler梁模型波纹之间的平箔片变形要更大，波箔之间的“凹陷”会更深，更容易造成气体的端部泄漏，减弱流体的动压效果，造成计算的承载力要比变形量更大的Timoshenko梁模型更小。

4 结语

本文基于Timoshenko梁建立了箔片模型和波箔片库仑摩擦模型，分析了箔片变形和轴承承载能力，得出以下结论：

1)Timoshenko梁模型的变形与NDOF模型的变形比较吻合，远小于Heshmat公式和Iordanoff公式计算结果；

2)平箔片和波箔片之间的摩擦与波箔片和轴承套之间的摩擦会提高波箔片的刚度，摩擦系数越大，箔片变形越小；

3)Timoshenko梁模型的箔片变形比Euler梁模型更大，但Euler梁模型平箔片的“凹陷”更明显，减弱了动压效果，造成Timoshenko梁模型承载力更高一些。