叶轮转子碰摩的非线性动力学响应

刘 昕， 张华彪， 孙小磊， 占传林， 赵庆军,3

(1．中国科学院 工程热物理研究所，北京 100190；2．中国科学院大学，北京 100190；3．中国科学院 轻型动力重点实验室，北京 100190)



叶轮转子碰摩的非线性动力学响应

刘 昕1,2， 张华彪1， 孙小磊1， 占传林1， 赵庆军1,3

(1．中国科学院 工程热物理研究所，北京 100190；2．中国科学院大学，北京 100190；3．中国科学院 轻型动力重点实验室，北京 100190)

叶轮转子的碰摩是旋转机械的常见故障，严重的碰摩可能导致转子失稳、叶片断裂等后果。针对叶轮转子-机匣的碰摩进行研究，将叶片简化为一端固支的悬臂梁，考虑其在旋转状态下的离心刚化效应，机匣的局部变形采用Hertz模型表征，推导了叶片-机匣接触力的表达式，分析了叶片长度、截面惯性矩、叶片截面积和转速对接触力的影响，结果表明当叶片长度越小、截面惯性矩越大、转速越高接触力越大。基于所得的接触关系，给出了叶轮转子碰摩力模型，对叶轮转子碰摩的稳态响应进行数值模拟，结果表明系统在碰摩时的响应频率主要包括转速激励频率、反向进动频率及其线性组合，叶片的扰动将导致主要响应频率附近出现边频，边频在不同转速情况下存在多种不同形式。当反向进动频率相对于转速近似线性变化时，边频和转速与叶片数目线性相关。

叶轮转子；碰摩；柔性机匣；边频

碰摩是叶轮机械常见故障之一，碰摩的发生可能导致机匣变形、轴承磨损、叶片断裂甚至转子失稳，最终造成严重的后果。

国内外学者针对转子系统的碰摩问题开展了广泛的研究。ZHANG等[1]分析了全周碰摩的精确解和解的稳定域。JIANG等[2]研究了同步全周/部分碰摩运动的稳定性。ZHANG等[3]采用了C-L方法对同步全周碰摩运动的分岔响应进行了分析。CHU等[4]利用Fourier分析和Floquet理论对跨中单转子的全局性分岔特性和稳定性进行研究，并通过数值模拟验证了通往混沌的三种路径。上述文献在研究中使用了圆盘转子模型，虽然便于分析，却不能很好的描述叶片和机匣的碰摩现象。PADOVAN等[5]将叶片简化为悬臂梁推导了叶片-刚性机匣碰摩力公式。PADOVAN等[6-8]针对不同情况的碰摩侵入量进行分析。JIANG等[9]和MA等[10]分别考虑转速和离心强化的影响，提出了不同的叶片机匣碰摩力模型，他们的研究更加侧重于叶片本身的响应。目前，针对工程中常见的叶轮转子-机匣碰摩问题的研究并不多见。

本文将叶片简化为悬臂欧拉梁，并考虑其在离心力作用下的刚度变化，利用Hertz模型描述叶片和机匣碰摩时的接触关系，构造了一种新的叶片-柔性机匣碰摩力的模型，并利用该模型对叶轮转子-柔性机匣碰摩的非线性动力学响应进行研究。

1 叶轮转子-机匣接触力模型

将叶片假设为一端固定于轮盘的等截面悬臂梁，图1中给出了叶片的受力情况。

图1 叶片简化模型Fig.1 Schematic diagram of rotating blade

在竖直方向存在离心力产生的分布载荷qω和法向接触力FN，水平方向存在摩擦力FT，用库伦模型描述，有

FT=μFN

(1)

式中：μ为摩擦因数。对分布载荷qω积分可得作用在距离叶根x位置处截面上的离心力，

Fω(x)=

(2)

式中：ρ,ω,A,Rr,lb分别为叶片密度，工作转速，叶片横截面积，轮盘半径和叶片展向长度。根据欧拉梁方程可得叶片的挠度方程，有

(3)

把离心力表达式(2)代入，展开，有

y″+(c+bx+ax2)y=hx

(4)

式中:

(5)

式(4)是一个非齐次非线性的微分方程，直接求解比较困难，这里采用幂级数法进行求解，设解为

(6)

考虑叶片叶根处挠度和转角为0，有y(0)=0,y′(0)=0。将式(6)代入式(4)，忽略高次项，令其余项系数为0，可解得

(7)

式中:n为解的阶次。下面对所得幂级数解进行分析，如未特殊声明，本文中的计算选取的参数由式(8)给出

ρ=8.1g/cm3,A=64mm2,Rr=73mm,μ=0.13

lb=45mm,I=210mm4,E=2.1×1011N/m2

(8)

图2给出了不同阶次的幂级数解与通过数值积分所得的挠度曲线的对比，其中横坐标为沿展向长度x，纵坐标为挠度y，可以看到对于本文所研究的叶片，5阶幂级数解和数值积分的误差已经小于1%，因此本文下面的计算取n=5。

图2 叶片挠度曲线Fig.2 Deflection curve of the blade

由式(6)和(7)，利用近似关系计算叶尖的径向变形量为

(9)

同时，考察转子运动过程中的转子位移与机匣变形量之间的关系，如图3。

图3 叶片-机匣碰摩几何位置关系Fig.3 Geometric model of contact between blades and casing

2δ(lb+Rr)cos(π-T))1/2

(10)

(11)

叶片和机匣的接触采用Hertz模型，接触力为

(12)

式中:α为与机匣和叶片材料相关的系数，反映了机匣的局部变形程度，其表达式为：

(13)

式中:E1,E2,η1,η2分别为机匣和叶片的弹性模量和泊松比，本文中分别都取2.1×1011N/m和0.2。

((lb+Rr)2+δ2-2δ(lb+Rr)cos(π-T))1/2-

Rc=α(FN)2/3+δb(FN)

(14)

式中:δb(FN)是关于FN的高次多项式。注意到上述表达式中左边项实际上就是叶片和机匣变形量之和，记做δtotal，而右边为一个关于FN代数多项式，使用叶片与机匣间的静态间隙c对δtotal进行无量纲化，有

(15)

(16)

2 叶片参数对接触力的影响

根据上面求得的接触力公式，给定不同叶片参数代入计算，可以得到叶片参数和转速对接触力的影响。

图4～图7分别给出了叶片截面积、叶轮转速、叶片截面惯性矩和叶片长度对接触力的影响。可以看到不同叶片截面积和不同转速条件下相同变形量对应的接触力差异不大，而截面惯性矩和叶片长度对接触力的影响较大。接触力随着叶片截面积的增大而略有增大。当转速较高时，离心强化效应更强，叶片-机匣接触的等效刚度更高，接触力更大。对于一定的总体变形量，叶片长度越大，截面惯性矩越小，叶片就越容易弯曲变形，对应的机匣的变形量就越小，而根据Hertz接触模型，较小的变形量对应较小的接触力，因此较大的叶片长度或较小的弯曲刚度对应较小的接触力。

图4 叶片截面积对接触力的影响Fig.4 Contact force under different cross sections of blade

图5 转速对接触力的影响Fig.5 Contact force under different angular speeds

图6 叶片截面惯性矩对接触力的影响Fig.6 Contact force under different moments of inertia

图7 叶片长度对接触力的影响Fig.7 Contact force under different lengths of blade

3 叶轮转子运动方程

考虑跨中叶轮转子的碰摩，叶轮转子-机匣系统的动力学方程可写作

(17)

式中:转子质量为m，机匣质量为M，Ce,Ke,eu分别为阻尼系数、支承刚度和不平衡量，角标r和c分别代表转子和机匣。Fx和Fy分别为碰摩力在x和y方向分量，当无碰摩时为0，发生碰摩时对单个叶片有

(18)

Xc=x,Yc=y,τ=ωnt

(19)

(20)

为避免对全周叶片逐个检测判断是否发生碰摩，叶片/机匣碰摩事件的检测采用PADOVA给出的方法[5]。图8给出了发生碰摩时叶轮转子和机匣的几何关系。根据转子进动角φ与自转角sτ，确定可能发生接触的中间位置的叶片的角位置T，从而可以得到中间位置叶片与机匣的距离Dc以及相邻两叶片对应的距离Dr和Dl，根据Dc、Dr和Dl即可判断叶片与机匣是否发生碰摩。具体方法为：① 若Dc>0，说明没有发生接触;② 若Dc<0，说明中心位置的叶片与机匣发生了接触;③ 再分别将两侧叶片与机匣的距离Dr和Dl视为Dc，重复步骤①，直到两侧叶片与机匣距离都大于0为止，从而得到所有与机匣接触叶片的侵入量。按照前面所述的方法求出每个叶片的碰摩力FN，并按对应角度位置求出总的矢量和，从而得到碰摩力Fx和Fy。将碰摩力代入方程，利用变步长数值积分方法即可对式(17)求解。

图8 转子机匣几何位置关系Fig.8 Geometric model of rotor and casing

4 叶轮转子-机匣碰摩响应

在本文后面的计算中，如未特殊声明，取阻尼系数Cre=2 000 Ns/m,Cce=6 000 Ns/m,，不平衡量eu=7×10-6m，叶尖与机匣的静态间隙值为c=6×10-5m，支撑刚度Kre=2.5×108N/m,Kce=5×108N/m，转子与机匣质量取m=M=10 kg。

作为对照，本文计算了圆盘转子与机匣的碰摩响应，接触力分别采用线性当量刚度模型和Hertz模型进行对比。由于本文中涉及到的叶片/机匣侵入量最大不超过间隙值c的15%，在此范围内对Hertz模型给出的侵入量-接触力曲线进行线性拟合，从而得到一个恒定刚度Kc近似代替Hertz模型。根据本文中的参数有Kc=1.5×109N/m，碰摩力表达式为

FN=Kcδtotal

(21)

式中:δtotal为叶片/机匣侵入量。图9给出了两种模型条件下接触力的对比。

图9 线性当量刚度模型与Hertz模型的碰摩力关系Fig.9 Contact force under linear model and Hertz model

图10和图11分别给出了当量刚度模型和Hertz模型条件下转子发生碰摩时稳态响应的频谱云图，横坐标为无量纲频率，s>0.95时转子和机匣发生碰摩，位移响应频谱云图上出现了非协调进动的频率成分。图12给出了圆盘转子碰摩的位移响应频谱和轴心轨迹。

图10 当量刚度模型圆盘转子碰摩响应频谱云图Fig.10 Spectrum contour of a rubbing 0-blade rotor with constant casing equivalent stiffness

图11 Hertz模型圆盘转子碰摩响应频谱云图(无叶片)Fig.11 Spectrum contour of a rubbing 0-blade rotor with Hertz model

图12 当量刚度模型圆盘转子的碰摩响应(s=1.04)Fig.12 Rubbing response of the no-blade rotor with constant casing equivalent stiffness (s=1.04)

可以看到图10中除了存在转速激励频率fω外，还存在转子反向进动频率fN及2fω-fN和2fN-fω等组合响应频率。将这些响应记作fO，则fO可统一写作如下形式

fO=n(fω-fN)+fN,

n=…,-2,-1,0,1,2,…

(22)

式中:n为整数。另外，在图10中fO的两侧还存在一对幅值更小的等间隔分布的频率fS。对比图11中Hertz模型机匣的转子系统的碰摩响应，可以看出两者的频谱基本一致。

图13给出了采用当量刚度模型的叶轮转子碰摩响应的频谱云图，可以看到叶轮转子的碰摩响应频谱和圆盘转子的频谱相比有较大区别。当s<0.97时，频谱和圆盘转子类似只存在fO与fS，而当s>0.97，fO与fS不再随着转速线性变化，且在这些频率附近出现不同形式的边频fside。图14给出了对应的叶轮转子碰摩的响应频谱和轴心轨迹。

图13 当量刚度模型叶轮转子碰摩响应频谱云图(叶片数目17)Fig.13 Spectrum contour of a rubbing 17-blade rotor with constant casing equivalent stiffness

图14 当量刚度模型叶轮转子的碰摩响应(叶片数目17，s=1.1)Fig.14 Rubbing response of the blade-rotor with constant casing equivalent stiffness (17 blades, s=1.1)

图15 当量刚度模型叶轮转子碰摩响应频谱云图局部，叶片数目17Fig.15 Details of the spectrum contour of a rubbing 17-blade rotor with constant casing equivalent stiffness

(23)

式中:n为整数,Nblade为叶片数目。图16给出了Hertz模型机匣的碰摩叶轮转子的频谱云图，注意到不同形式的边频依然存在，说明边频并不依赖与特定的碰摩力模型，而是叶轮碰摩转子的独特性质。

图16 Hertz模型叶轮转子碰摩响应频谱云图，叶片数目17Fig.16 Spectrum contour of a rubbing 17-blade rotor with Hertz model

5 结 论

本文提出了一种新的叶轮转子叶尖-机匣碰摩力模型，此模型考虑了高速旋转叶片的离心强化作用，使用Hertz模型表征机匣局部变形情况下叶片和机匣的接触。对该模型的分析表明，在叶片对机匣的侵入量相同的情况下，叶片展向长度越短、截面惯性矩越大、工作转速越高，单个叶片上承受载荷越大。

基于叶轮转子叶尖-机匣碰摩力模型，建立了叶轮转子-机匣碰摩的运动方程并进行了数值模拟，计算结果表明在叶片与机匣发生碰摩时，位移响应频谱上存在反向进动频率和转速激励频率的线性组合，对于叶轮转子还在主要响应频率附近还存在较小的边频响应。边频在不同转速情况下存在多种不同形式。当反向进动频率相对于转速近似线性变化时，边频与主要响应频率之差相对于转速线性变化，在位移响应频谱云图上的边频的斜率和主要响应频率的斜率之差与叶片数目之比为有理数，体现了叶轮转子碰摩的独特性质。

[1] ZHANG G F, XU W N, XU B, et al. Analytical study of nonlinear synchronous full annular rub motion of flexible rotor-stator system and its dynamic stability[J]. Nonlinear Dynamics, 2009, 57(4): 579-592.

[2] JIANG J, ULBRICH H. Stability analysis of sliding whirl in a nonlinear Jeffcott rotor with cross-coupling stiffness coefficients[J]. Nonlinear Dynamics, 2001, 24(3): 269-283.[3] ZHANG H, CHEN Y. Bifurcation analysis on full annular rub of a nonlinear rotor system[J]. Science China Technological Sciences, 2011, 54(8): 1977-1985.

[4] CHU F, ZHANG Z. Bifurcation and chaos in a rub-impact Jeffcott rotor system[J]. Journal of Sound and Vibration, 1998, 210(1): 1-18.

[5] PADOVAN J, CHOY F. Nonlinear dynamics of rotor/blade/casing rub interactions[J]. Journal of Turbomachinery, 1987, 109(4): 527-534.

[6] PADOVA C, BARTON J, DUNN M G, et al. Experimental results from controlled blade tip/shroud rubs at engine speed[J]. Journal of Turbomachinery, 2007, 129(4): 713-723.

[7] PADOVA C, BARTON J, DUNN M G, et al. Development of an experimental capability to produce controlled blade tip/shroud rubs at engine speed[J]. Journal of Turbomachinery, 2005, 127(4): 726-735.

[8] PADOVA C, DUNN M, BARTON J, et al. Controlled fan blade tip/shroud rubs at engine conditions[C]// ASME 2011 Turbo Expo: Turbine Technical Conference and Exposition. Vancouver,Canada, 2011: 951-961.

[9] JIANG J, AHRENS J, ULBRICH H, et al. A contact model of a rotating rubbing blade[C]// Proceedings of the 5th International Conference on Rotor Dynamics, echnische Universitat Darmstadt. Darmstadt, Germany, 1998: 478-489.

[10] MA H, TAI X, HAN Q, et al. A revised model for rubbing between rotating blade and elastic casing[J]. Journal of Sound and Vibration, 2015,337:301-320.

[11] JIANG J. The analytical solution and the existence condition of dry friction backward whirl in rotor-to-stator contact systems[J]. Journal of Vibration and Acoustics, 2007, 129(2): 260-264.

Nonlinear dynamics on rubbing of a blade-rotor system

LIU Xin1, 2, ZHANG Huabiao1, SUN Xiaolei1, ZHAN Chuanlin1, ZHAO Qingjun1,3

(1. Institute of Engineering Thermophysics, Chinese Academy of Sciences, Beijing100190, China;2. University of Chinese Academy of Sciences, Beijing100190, China;3. Key Laboratory of Light-duty Gas-turbine, Chinese Academy of Sciences, Beijing100190, China)

As a type of common failure of the rotating machines, rubbing can lead to catastrophic consequences such as instability of the rotor or blade damages. Rubbing between blades and casing was investigated in this paper. The blades were simplified as cantilever beams with centrifugal stiffening effect, while the blade/casing contact was modeled by Hertz contact force and an expression of blade/casing contact force has been derived. Based on the derived expression, influence of various parameters on the contact force was analyzed, such as the rotating speed, cross section of blade, length of blade, and bending stiffness of blade. The results show that blades with shorter length and larger moment of inertia or under higher rotating speed will lead to greater contact force on the casing. Meanwhile, a blade-rotor/casing rubbing model was derived with the given expression of contact force and numeric simulations with this model were also conducted. The results show that the responses on the spectrum consist of synchronous frequency, reverse precession frequencies, and their linear combinations. It was also observed that the blades could lead to side frequencies near the main frequencies. Different kinds of side frequencies were observed under different angular speeds. When the reverse precession frequencies are nearly linear with the angular speed, the side frequencies are also linear with the angular speed as well as the number of blades.

blade-rotor; rubbing, flexible casing; side frequencies

国家自然科学基金(11302223)；国家科技支撑计划资助项目(2012BAA11B01)

2015-06-11 修改稿收到日期：2015-09-24

刘昕 男，硕士生，1991年生

张华彪 男，博士，助理研究员，1984年生

TH113.1

A

10.13465/j.cnki.jvs.2016.20.028