黎荣健,肖曙红,张岩

(广东工业大学 机电工程学院,广州 510006)

与传统的油润滑轴承和滚动轴承相比,箔片气体轴承因结构简单,使用寿命长,维护成本低以及能在高速、高温等极端环境下运行[1-2],成为高速涡轮机械的绝佳选择[3]。在已开发的箔片轴承类型中,径向箔片气体轴承使用最广泛,轴套内部放置一个柔性箔结构, 该结构由一个波箔和一个顶箔组成,这种弹性支承结构增强了轴承的承载能力和阻尼[4],克服了气体黏度低带来的性能缺陷。

快速、准确预测轴承性能对于箔片气体轴承的开发至关重要,然而,由于摩擦使得轴承系统具有高度非线性,仍然很难准确描述轴承的工作特性。文献[5-6]用一种弹簧模型模拟箔片,将各个波箔片等效为无相互作用的弹簧,但忽略了顶箔的刚度,箔片之间的摩擦和波箔之间的相互作用,导致结构刚度被低估。不过该模型因其结构简单,易于编程仍被广泛使用[7-9],文献[10]就利用此类模型在研究中首次考虑了箔片结构的阻尼效应,并预测了箔片轴承的刚度和阻尼系数。虽然上述模型已被广泛采用,但由于转轴与顶箔之间碰撞的相互作用被丢弃,因而具有局限性。为获得更精确的模拟结果,开始使用有限元法对轴承进行建模:文献[11]使用二维梁单元对顶箔和波箔建模,并使用黏滑算法[12]考虑摩擦,建立了考虑滞回特性的模型,预测表明一个静载荷对应多个平衡位置;文献[13]同样采用梁单元为基础的有限元结构模型,但不同之处在于采用非线性接触的数值程序处理各接触区域的接触力,并将该方法应用于三维模型[14]。另外还可使用壳单元模拟箔片[15-18],如文献[15]提出了一个考虑大位移的非线性模型,并将箔片建模为壳单元,文献[18]以惩罚法和正则化摩擦定律为基础对波箔进行建模。还有使用有限元商业代码[19-20]对箔片进行建模的方法,其考虑了箔片之间的相互作用,但求解较为复杂,在建立非线性摩擦模型时数值收敛性不好。

为简化计算量并保证计算精度,开发了近似有限元结构的桁架系统。文献[21]用一种结构简单的桁架系统模拟波箔片,每个波箔片被建模为3个基本弹簧,弹簧的刚度由卡斯提利亚诺理论获得,虽然结构简单,但与有限元仿真结果表现出更好的一致性,还证明了波箔相互作用的重要性。文献[22-23]采用上述桁架模型,在波箔与轴套的接触部位增加了垂直自由度,并考虑了顶箔的影响和轴承系统的接触/分离行为,箔片结构的法向和切向接触力分别采用增广拉格朗日乘子法和惩罚法处理,分析结果与文献[24-25]一致。

为准确描述摩擦对轴承系统特性的影响,本文提出一种考虑非线性接触的静态结构模型,其能捕捉柔性结构的非线性与时变性;主要讨论波箔与轴套之间以及箔片之间的点线接触状态;利用惩罚法离散每个接触区域的力学向量,通过正则化技术改进库伦摩擦模型,采用牛顿-拉弗森法求解平衡方程;最后通过转子推拉仿真得到箔片结构的力学行为与能量耗散等特性。

1 理论模型

1.1 波箔建模

以轴心为原点建立全局坐标系OXY,则径向箔片气体轴承简化图如图1所示,为简化计算量,波箔采用文献[21]提出的简化桁架结构进行建模。对波箔建立局部坐标系Oxy,2个波箔片组成的箔片数学模型如图2所示,可以扩展到任意数量波箔片,因为考虑了波箔、顶箔与轴套之间的松紧接触,在每个凸块底部都需要增加一个额外的自由度(图2中蓝色箭头),即垂直位移,并且需要重新求解波箔的刚度。

图1 径向箔片气体轴承简化图

图2 2个波箔片组成的箔片数学模型

(1)

进一步转化为全局坐标系的局部刚度矩阵kb,即

(2)

(3)

c=cosβ,

s=sinβ,

式中:T为坐标转换矩阵;β为两坐标系之间的夹角。

1.2 顶箔建模

(4)

式中:E为材料弹性模量,GPa;A为单元横截面积,m2;I为横截面积的惯性矩,m4;Le为单元长度,m。

1.3 摩擦接触建模

顶箔和轴套视为主体,波箔视为从体。在外力作用下,箔片的变形增加了箔片之间的接触面积,理论上,每个波箔与顶箔之间应建立多个接触点对,由于与箔片的长度相比变形的尺寸较小,因此仅选择每个波箔上的最高节点作为接触节点。箔片结构及接触点示意图如图3所示,波箔上的绿色节点为接触节点,也被视为检测点,被检测节点位于红色单元中,蓝色节点是与接触对相关的节点。

图3 箔片结构及接触点示意图

顶箔与波箔之间的接触示意图如图4所示,Q为波箔的接触节点,P为Q在顶箔上的投影,P与Q之间的相对位移为(下角标t表示切向方向,n表示法向方向)

图4 顶箔与波箔之间的接触示意图

uP-uQ=Ncuc,

(5)

Nc=[-N1(ξ) -N2(ξ)I2×2],

上述向量在全局坐标系中定义,为方便计算,将其引入到接触条件中,在接触点对的局部坐标系中重新定义(5)式,即

(6)

L=-I2×2,

(7)

式中:uB,uA分别为P,Q在局部坐标系中的位移向量;L为坐标变换矩阵。

将接触面离散化后,每个接触区域的分布接触力转化为离散形式,节点Q的法向接触力Fn表示为

(8)

式中:α为惩罚数,实际计算中不可能将α设置为无穷大以获得准确解,需要根据实际问题选择一个可接受的值,本文α取1.0×108N/mm。

利用库伦摩擦模型将节点Q的切向接触力Ft表示为

(9)

(10)

式中:c1为控制正则化摩擦模型与库伦摩擦模型之间接近度的平滑参数,本文c1取100 m。

结合(9)—(10)式得到Q点切向力为

(11)

根据反作用力原理可得到相关节点对和接触节点对的等效节点接触力向量,即

(12)

FB=[FtFn]T,

式中:FB为B点的接触力向量;qc等效节点接触力向量。

由于考虑了摩擦接触,非线性行为被引入模型,假设将所有单元组合后形成的系统的全局非线性方程为

K(u)=FL=Ku-Qc,

(13)

式中:K(u)为u的非线性函数;FL为外力向量。

牛顿-拉弗森法是求解非线性方程的基本方法之一,其从一个假设解开始,通过不断的迭代逼近真实解,直到满足收敛准则。假设第l步的近似解已知,记作ul,则下一步的近似解ul+1近似为一阶泰勒展开式,即

(14)

进一步写为

(15)

(16)

Qc=[(qc,1)T(qc,2)T… (qc,Ncont)T]T,

(17)

D=[(uc,1)T(uc,2)T… (uc,Ncont)T]T,

(18)

(19)

最后得到线性化的系统方程为

(20)

(21)

通过(20)式得到增量Δul,则新的近似解为ul+1=ul+Δul,通常该解满足不了原非线性方程,当残差变得足够小时,u是正确解。平衡方程的收敛准则为

(22)

2 结果与分析

转子不对中时,随着转轴在垂直方向移动,箔片结构在垂直和水平方向都发生变形,箔片的垂直变形决定了轴承的承载能力,而水平运动产生阻尼。箔片气体轴承结构参数见表1,利用Matlab将数学模型转化为程序语言进行仿真,转轴的最大径向位移设为79.5 μm(即2.5C),程序模拟转轴在正负X,Y方向上逐步压向箔片结构,当转轴位移超过C时,转轴的表面和箔片结构开始相互作用,达到最大径向位移后,转轴逐渐被拉回初始位置。

表1 箔片气体轴承结构参数

2.1 箔片结构力学行为

不同摩擦因数下X,Y方向静载荷与转轴位移的关系如图5所示:除-X方向外仿真结果与文献[20]完全非线性模型的结果非常吻合,因为仿真程序阐明了加载和卸载,所以预期的滞回曲线清晰可见;静载荷为转轴挤压顶箔时受到顶箔的反作用力,方向为转轴径向移动方向,并与接触点处切线垂直,加载过程的静载荷比卸载过程的静载荷大,这是因为卸载开始时一部分力平衡了摩擦力,接触点不会立即滑动;载荷与偏转曲线在轻载荷区域不是线性的,而在重载荷区域几乎是直线,这是因为在轻载荷区域只有少数波箔发生变形, 当变形波箔的数量不再增加时,载荷与偏转曲线的斜率几乎恒定; 随着摩擦因数的增大,环路变得越来越明显,当没有摩擦时,加载过程曲线与卸载过程曲线完全重合;在-X方向上,仿真预测结果明显小于完全非线性模型预测结果,文献[22]指出,虽然固定端附近的几个波箔的刚度非常大,但自由端附近的波箔刚度非常小,因此完全非线性模型可能会高估-X方向上的结构刚度。

图5 不同摩擦因数下静载荷与转轴位移的关系

摩擦因数分别为0,0.10,转轴达到最大径向位移(2.5C)时4个方向波箔上的静载荷仿真结果如图6所示:波箔从固定端开始编号,摩擦因数对各波箔之间刚度的相对大小有影响;对于-X方向,靠近固定端的第1和第2个波箔比其他波箔的静载荷大得多,自由端附近波箔的静载荷很小。

图6 转轴位移为2.5C时各波箔的承载力

2.2 固定端对箔片刚度的影响

选择+X方向分析固定端对箔片刚度的影响,如图7所示,随着转轴在+X方向上被推拉(最大径向位移为2.5C),第10到第18个波箔与转轴表面接触,波箔14的垂直位移最大,波箔14左右两侧各选一个波箔(波箔12与16)进行对比分析。

图7 波箔与转轴表面接触

转轴在+X方向上被推拉至径向位移为2.5C的过程中,波箔12,14,16的滞回曲线如图8所示:无论加载还是卸载过程,靠近固定端的波箔具有更大的刚度;在整个回路中,波箔垂直位移与静载荷之间呈线性关系,符合实际情况。转轴最大径向位移分别为60,70,80 μm时波箔14的滞回曲线如图9所示:在卸载过程中曲线不重合,这是固定端施加反作用力的结果。通过上述仿真结果可以看出,固定端对箔片结构力学性能的影响不能忽略。

图8 波箔12,14,16的滞回曲线

图9 转轴最大径向位移分别为60,70,80 μm时波箔14的滞回曲线

转轴最大径向位移分别为60,70,80 μm时波箔12与波箔16的滞回曲线如图10所示:波箔12更接近于固定端,在大回路中表现出更高的刚度。这是因为仿真时转轴中心左侧的波箔受固定端响应力的影响较大;在卸载过程中,固定端的响应力有助于箔片克服摩擦力,促进接触节点的滑动,预载荷较小时该现象更为明显;随着预载荷逐渐增大,固定端的响应力对波箔的影响变小,这是因为转轴位移增加到一定程度时反作用力比摩擦力小很多,从而导致预载荷较小的转轴在卸载过程中更快进入滑动状态。

转轴径向位移为2.5C,摩擦因数为0.1,不同方向上波箔的耗散能量如图11所示,在图11的基础上继续研究不同摩擦因数下箔片的总耗散能量,结果如图12所示:影响耗散能量的因素为静载荷和接触点的水平位移,能够看出各波箔耗散能量与图6中的静载荷分布略有不同;随着摩擦因数的增大,能量耗散越来越多;+Y方向的耗散能量最高,-X方向上靠近固定端的2个波箔能量耗散较高,但是靠近自由端的波箔刚度很低,所以总耗散能量偏低。

图12 不同摩擦因数下箔片的总耗散能量

3 结论

介绍了一种考虑非线性接触的箔片轴承分析模型,波箔使用桁架结构建模,顶箔使用梁单元建模;顶箔、波箔与轴套之间的接触被描述为点线接触,并通过基于惩罚法的接触算法进行说明;为提高数值计算的稳定性,采用正则化摩擦模型代替库伦摩擦模型;牛顿-拉夫森方法用于求解每个增量步骤中的系统方程,得到主要结论如下:

1)利用该模型研究了箔片结构的力学行为,通过转轴推拉仿真得到预期的滞回曲线,证明了摩擦是箔片产生滞后现象的主要因素。

2)分析转轴推拉仿真达到最大径向位移时+X,-X,+Y,-Y方向上各个波箔的静载荷分布,发现摩擦能够提升箔片的承载能力,摩擦因数对各个波箔之间刚度的相对大小有一定的影响。

3)观察单个箔片的滞回曲线,发现靠近固定端的箔片刚度更高,受固定端反作用力的影响,不同最大径向位移箔片的卸载曲线不重合。不同摩擦因数下各个方向箔片结构的耗散能量表明,相同外加条件下+Y方向的耗散能量最大。