基于接触非线性的箔片结构力学特性分析

摘要：箔片轴承中顶箔与波箔之间的非线性接触力学问题具有非线性和时变性，建立考虑非线性接触的理论模型研究箔片结构的力学行为，采用增量迭代法和牛顿-拉弗森法求解系统平衡方程，通过该模型进行转子推拉仿真获得箔片结构变形量与响应力之间的关系以及箔片的能量耗散。结果表明:摩擦是箔片结构产生滞回行为的主要因素，随着摩擦因数的增大，箔片结构刚度增大，滞回曲线面积也随之增大；在一个加载、卸载回路中，加载过程比卸载过程具有更高的静载力；靠近固定端的波箔片比靠近自由端的波箔片具有更高的刚度。关键词：滑动轴承;箔片轴承;气体动压轴承;滞回曲线;能量耗散与传统的油润滑轴承和滚动轴承相比,箔片气体轴承因结构简单,使用寿命长,维护成本低以及能在高速、高温等极端环境下运行[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径向箔片气体轴承简化图Fig.1Simplifieddiagramofradialfoilgasbearing图22个波箔片组成的箔片数学模型Fig.2Amathematicalmodeloffoilcomposedoftwobumpfoils弹簧刚度与文献[21]中一致,每个弹簧均被视为基本单元。由于弹簧刚度ki已经量化,使用直接刚度法更容易求解单元的刚度矩阵,因此局部坐标系中单元的刚度矩阵可表示为(1)进一步转化为全局坐标系的局部刚度矩阵kb,即(2)(3)c=cos

β,s=sin

β,式中:T为坐标转换矩阵;β为两坐标系之间的夹角。1.2顶箔建模顶箔由欧拉-伯努利梁单元进行建模,每个波箔片对应的顶箔被离散为3个单元,每个单元有2个节点和3个自由度,展开的顶箔单元水平放置在桁架结构上,因此顶箔单元的全局坐标系与局部坐标系重合。在全局坐标系中,顶箔单元的刚度矩阵ktop表示为顶箔的局部单元刚度矩阵由伽辽金法获得,即(4)式中:E为材料弹性模量,GPa;A为单元横截面积,m2;I为横截面积的惯性矩,m4;Le为单元长度,m。箔片的整体刚度矩阵K=diag(Kb,Ktop)通过单元节点的编号组合单元刚度矩阵获得,其中Kb=∑kb,Ktop=∑ktop。全局结构的广义位移坐标向量其中ub为波箔单元节点的位移坐标向量,ut为顶箔单元节点的位移坐标向量。1.3摩擦接触建模顶箔和轴套视为主体,波箔视为从体。在外力作用下,箔片的变形增加了箔片之间的接触面积,理论上,每个波箔与顶箔之间应建立多个接触点对,由于与箔片的长度相比变形的尺寸较小,因此仅选择每个波箔上的最高节点作为接触节点。箔片结构及接触点示意图如图3所示,波箔上的绿色节点为接触节点,也被视为检测点,被检测节点位于红色单元中,蓝色节点是与接触对相关的节点。图3箔片结构及接触点示意图Fig.3Diagramoffoilstructureandcontactpoints顶箔与波箔之间的接触示意图如图4所示,Q为波箔的接触节点,P为Q在顶箔上的投影,P与Q之间的相对位移为(下角标t表示切向方向,n表示法向方向)图4顶箔与波箔之间的接触示意图Fig.4DiagramofcontactbetweentopfoilandbumpfoiluP-uQ=Ncuc,(5)Nc=[-N1(ξ)-N2(ξ)

I2×2],式中:uP为P点位移向量,需要通过节点P1和P2的位移向量插值得到;uP,uP分别为节点P1,P2的位移向量,为转角向量;uQ为Q点的位移向量;Nc为形状函数,利用插值法得到相应坐标;uc为形状位移;Ni(ξ)为梁单元的形函数;ξ为P点对应量纲一的坐标,ξ(P1)=-1,ξ(P2)=1;I为单位矩阵。上述向量在全局坐标系中定义,为方便计算,将其引入到接触条件中,在接触点对的局部坐标系中重新定义(5)式,即(6)L=-I2×2,(7)式中:uB,uA分别为P,Q在局部坐标系中的位移向量;L为坐标变换矩阵。将接触面离散化后,每个接触区域的分布接触力转化为离散形式,节点Q的法向接触力Fn表示为(8)式中:α为惩罚数,实际计算中不可能将α设置为无穷大以获得准确解,需要根据实际问题选择一个可接受的值,本文α取1.0×108

N/mm。利用库伦摩擦模型将节点Q的切向接触力Ft表示为(9)式中:分别为当前加载步骤与前一个加载步骤的切向相对位移,为摩擦因数。时,接触点处于黏着状态,时,接触点在大小和方向上存在不受限制的相对滑动。由于接触状态的不连续,数值计算难以收敛,因此使用符号函数来近似平滑不连续,即(10)式中:c1为控制正则化摩擦模型与库伦摩擦模型之间接近度的平滑参数,本文c1取100m。结合(9)—(10)式得到Q点切向力为(11)根据反作用力原理可得到相关节点对和接触节点对的等效节点接触力向量,即(12)FB=[Ft

Fn]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,N)T]T,(17)D=[(uc,1)T

(uc,2)T

…(uc,N)T]T,(18)式中:Δul为位移向量的增量;为第l步的全局切向刚度矩阵;Qc,D分别为所有接触点对与相关节点的接触力向量和位移向量。可以通过使用相应的选择矩阵从全局向量中提取。任一接触点对的等效节点接触力qc的导数为(19)最后得到线性化的系统方程为(20)(21)通过(20)式得到增量Δul,则新的近似解为ul+1=ul+Δul,通常该解满足不了原非线性方程,当残差变得足够小时,u是正确解。平衡方程的收敛准则为(22)2结果与分析转子不对中时,随着转轴在垂直方向移动,箔片结构在垂直和水平方向都发生变形,箔片的垂直变形决定了轴承的承载能力,而水平运动产生阻尼。箔片气体轴承结构参数见表1,利用Matlab将数学模型转化为程序语言进行仿真,转轴的最大径向位移设为79.5μm(即2.5C),程序模拟转轴在正负X,Y方向上逐步压向箔片结构,当转轴位移超过C时,转轴的表面和箔片结构开始相互作用,达到最大径向位移后,转轴逐渐被拉回初始位置。表1箔片气体轴承结构参数Tab.1Structuralparametersoffoilgasbearings2.1箔片结构力学行为不同摩擦因数下X,Y方向静载荷与转轴位移的关系如图5所示:除-X方向外仿真结果与文献[20]完全非线性模型的结果非常吻合,因为仿真程序阐明了加载和卸载,所以预期的滞回曲线清晰可见;静载荷为转轴挤压顶箔时受到顶箔的反作用力,方向为转轴径向移动方向,并与接触点处切线垂直,加载过程的静载荷比卸载过程的静载荷大,这是因为卸载开始时一部分力平衡了摩擦力,接触点不会立即滑动;载荷与偏转曲线在轻载荷区域不是线性的,而在重载荷区域几乎是直线,这是因为在轻载荷区域只有少数波箔发生变形,当变形波箔的数量不再增加时,载荷与偏转曲线的斜率几乎恒定;随着摩擦因数的增大,环路变得越来越明显,当没有摩擦时,加载过程曲线与卸载过程曲线完全重合;在-X方向上,仿真预测结果明显小于完全非线性模型预测结果,文献[22]指出,虽然固定端附近的几个波箔的刚度非常大,但自由端附近的波箔刚度非常小,因此完全非线性模型可能会高估-X方向上的结构刚度。

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

图6转轴位移为2.5C时各波箔的承载力Fig.6Loadcapacityofeachbumpfoilwhenshaftdisplacementis2.5C2.2固定端对箔片刚度的影响选择+X方向分析固定端对箔片刚度的影响,如图7所示,随着转轴在+X方向上被推拉(最大径向位移为2.5C),第10到第18个波箔与转轴表面接触,波箔14的垂直位移最大,波箔14左右两侧各选一个波箔(波箔12与16)进行对比分析。图7波箔与转轴表面接触Fig.7Contactsofbumpfoilswithshaft转轴在+X方向上被推拉至径向位移为2.5C的过程中,波箔12,14,16的滞回曲线如图8所示:无论加载还是卸载过程,靠近固定端的波箔具有更大的刚度;在整个回路中,波箔垂直位移与静载荷之间呈线性关系,符合实际情况。转轴最大径向位移分别为60,70,80μm时波箔14的滞回曲线如图9所示:在卸载过程中曲线不重合,这是固定端施加反作用力的结果。通过上述仿真结果可以看出,固定端对箔片结构力学性能的影响不能忽略。图8波箔12,14,16的滞回曲线Fig.8Hystereticcurvesofbumpfoils12,14,and16图9转轴最大径向位移分别为60,70,80μm时波箔14的滞回曲线Fig.9Hystereticcurvesofbumpfoil14whenmaximumradialdisplacementofshaftis60,70and80μmrespectively转轴最大径向位移分别为60,70,80μm时波箔12与波箔16的滞回曲线如图10所示:波箔12更接近于固定端,在大回路中表现出更高的刚度。这是因为仿真时转轴中心左侧的波箔受固定端响应力的影响较大;在卸载过程中,固定端的响应力有助于箔片克服摩擦力,促进接触节点的滑动,预载荷较小时该现象更为明显;随着预载荷逐渐增大,固定端的响应力对波箔的影响变小,这是因为转轴位移增加到一定程度时反作用力比摩擦力小很多,从而导致预载荷较小的转轴在卸载过程中更快进入滑动状态。

图10波箔12,16的垂直位移与静载荷的关系Fig.10Relationshipbetweenverticaldisplacementandstaticloadofbumpfoils12and16转轴径向位移为2.5C,摩擦因数为0.1,不同方向上波箔的耗散能量如图11所示,在图11的基础上继续研究不同摩擦因数下箔片的总耗散能量,结果如图12所示:影响耗散能量的因素为静载荷和接触点的水平位移,能够看出各波箔耗散能量与图6中的静载荷分布