2019, Vol. 40
2. 东北大学 航空动力装备振动及控制教育部重点实验室, 辽宁 沈阳 110819
2. Key Laboratory of Vibration and Control of Aero-Propulsion System, Ministry of Education, Northeastern University, Shenyang 110819, China
转子系统的滞后突跳行为是指在非线性因素的激励下, 转子系统在水平或竖直方向上出现与Duffing软硬特性曲线类似的现象, 即共振峰值发生移动导致转子系统在不同于临界转速的位置发生大幅振动从而引发的突跳行为.滚动轴承-转子系统中作为主要支承构件的滚动轴承具有诸多非线性因素, 由此导致转子系统在运行过程中存在较强的非线性特性, 这对于转子系统运行的稳定性和安全性具有极大的影响.滚动轴承的非线性因素包括轴承赫兹接触刚度、轴承游隙等, 这些非线性因素引起转子系统诸多复杂的运动行为.因此, 国内外学者针对滚动轴承-转子系统开展了大量的研究工作.
Yamamoto等[1]通过建立球轴承-转子系统动力学模型, 针对滚动轴承间隙的非线性影响进行研究. Tiwari等[2]研究了轴承球体径向内部间隙对转子平衡动力学的影响.Villa等[3]针对复杂柔性转子轴承系统进行了稳定性分析.Peletan等[4]针对转子系统的稳态与动态特性提出准周期谐波平衡法(HBM)与伪弧长延拓法相结合的算法.Sun等[5]研究了具有摩擦冲击的双转子系统的稳态响应, 其研究结果表明了质量、偏心率、轴间刚度和转速比等参数对双转子系统动态特性的影响.Sundararajan等[6]用弧长延伸法求得了转子系统的动态特性.Ganesan[7]分析了转轴及轴承相关参数对转子系统振动特性的影响.Wang等[8]研究不同轴承配置对转子系统稳定性的影响, 并提出一种精确识别系统阻尼比的方法.张智勇[9]构建两自由度深沟球轴承模型, 以HB-AFT为手段对轴承转子系统进行追踪.白长青等[10]研究了轴承-转子系统在不同轴承间隙、不同转速下的稳定性, 其研究表明，轴承间隙是影响转子系统稳定性的一个重要因素.Ma等[11]建立了转子叶片系统的新动力学模型并利用有限元法和实验验证了该模型的正确性.Han等[12]研究了柔性转子-轴承系统在时间周期基角运动下的参数不稳定性并推导了时变基角运动下转子系统的精确有限元模型.
本文旨在利用同伦弧长延拓的HB-AFT法并结合Floquent稳定性理论, 计算转子系统滞后突跳区间, 分析滚动轴承的非线性特性对转子系统滞后区间的影响.
1 滚动轴承-转子模型为研究轴承非线性特性对转子系统振动特性的影响, 如图 1所示建立滚动轴承-转子系统动力学模型, 并作如下假设:①轴承外圈固结于轴承座, 且轴承座始终处于静止状态; ②转盘在转轴的正中位置, 两端轴承呈对称分布, 整个转子系统呈刚性; ③滚动轴承内圈与转轴整体建模, 滚动体等效为具有刚度和阻尼的弹簧, 并保留轴承的游隙、赫兹接触、变柔度这三个主要非线性因素.
|
图 1 滚动轴承-转子系统动力学模型 Fig.1 Dynamic model of the rolling bearing rotor system |
图 1中, 转盘几何中心O为坐标轴Oxyz原点；Oe为质心, e为偏心距；FL, x, FL, y, FR, x, FR, y为左、右轴承产生的非线性支反力.转子系统振动微分方程为
|
(1) |
式中: 
由于在建立滚动轴承-转子系统时将滚动轴承内圈与转轴进行整体建模, 且转子系统为刚性轴, 故系统振动响应同时反映了转子部分振动特性及滚动轴承振动特性.
2 计算转子系统的滞后突跳区间由于滚动轴承赫兹接触与轴承间隙的存在, 转子系统在运行过程中会出现有线性恢复力与无线性恢复力的情况, 因此, 转子系统的动力学微分方程属于分段非线性方程, 普通的谐波平衡法不适用于本文建立的转子系统非线性微分方程.针对此类问题, Yamauchi提出了一种半数值半解析的隐式谐波平衡法[13], 其主要思想是将系统响应及非线性部分分别设成周期响应解, 并对其进行离散化处理, 即将1个周期T分为N份, 每段时间Δt=T/N, 根据系统方程分别在每个离散点处进行频域时域转换, 建立谐波系数之间的关系, 并进行迭代求解, 从而得到相应条件下的系统响应.此方法称为频时域转换谐波平衡法(简称HB-AFT), 其优点在于避免了求解过程中复杂的积分处理, 对于复杂非线性系统求解具有较高的实用性.
2.1 HB-AFT的一般形式转子系统动力学方程写成矩阵形式
|
(2) |
式中：M表示转子系统质量矩阵；D表示转子系统阻尼矩阵；K表示系统刚度矩阵；X为动力学方程的解；
设动力学方程的解X与非线性作用力FNL为傅里叶级数形式:
|
(3) |
|
(4) |
式中：K为谐波项数；A0, C0, Ak, Bk, Ck, Sk为谐波系数.
对式(3)求取一阶、二阶导数可得
|
(5) |
将式(3)~式(5)代入式(2), 使等式两端对应的常数项及三角函数项系数分别相等, 常数部分为
|
(6) |
第k次谐波项结果为
|
(7) |
|
(8) |
将所有谐波项用对角矩阵表示, 可得(2K+1)r维矩阵:
|
(9) |
式(9)中
|
(10) |
式中r为转子系统自由度数目.
FNL为关于包含位移的非线性项的隐函数的傅里叶表达式, 因此, 式(9)可以写成如下形式:
|
(11) |
P, Q系数均为未知数时, 该系统方程含有(4K+2)r个未知数, 而式(9)中仅包含(2K+1)r个代数方程组, 无法直接对P, Q进行求解, 因此需要建立非线性力与位移的隐函数关系, 即令Q中的未知数用P中未知数表示.
|
(12) |
Q由C0, Ck, Sk组成, FNL与P中元素相关, 通过式(12)建立Q与P的关系, 得(2K+1)r个未知数和(2K+1)r个方程组.给定P一组初始值, 利用Newton-Raphson方法迭代得到P, 具体迭代格式如下:
|
(13) |
|
(14) |
式(14)中的∂H/∂P，∂H/∂Q可通过式(9)、式(10)得到, dQ/dP可通过式(12)求导得到.为避免积分过程, 采用IDFT方法将时域信息离散为以下形式[9]:
|
(15) |
式中, n=0, 1, 2, …, N－1.这里X(n, P)表示由位移谐波系数P表示的时域位移响应X(τ)在离散时间点n·ΔT处的值, 这里ΔT=Tω/N, Tω为ω对应的周期, N为时域离散点数.则非线性力FNL时域离散信息为
|
(16) |
经DFT变换，式(16)可表示为
|
(17) |
|
(18) |
|
(19) |
式中, k=1, 2, 3, …, K.通过式(17)~式(19), 建立了由时域信息和频域信息转换得来的Q与P的代数关系式, 得到式(14)中的dQ/dP, 进而可通过Newton-Raphson迭代过程求解未知向量P.
迭代求解不动点P的HB-ATF过程如下：
① 设置迭代初值P(0)及各项初始参数;
② 采用AFT过程给出迭代Jacobian矩阵J(0);
③ 根据式(14)中的Newton-Raphson方法进行迭代, 求得P(1);
④ 将P(1)值赋予P(0), 循环步骤②和③, 直到m次迭代后P(m)－P(m－1)的范数小于指定精度ε, 也可通过判定H范数小于ε来得到最优解.具体流程如图 2所示.
|
图 2 AFT计算方法流程图 Fig.2 Flow chart of AFT calculation method |
由以上求解过程可见, 非线性部分被设置为谐波解形式后, 只要非线性部分与系统位移响应函数关系确定且存在导数关系, 则非线性部分的复杂程度并不影响离散HB-AFT方法的使用.如式(16)所示, 复杂的非线性力都是由未知自变量的代数关系式确定的, 那么就可以根据式(13)~式(19)求得Newton-Raphson迭代计算所需的Jacobian矩阵J.因此, 即使对于包含分段函数、间隙和分数指数的非线性参激系统, HB-AFT方法求解依然有效.
2.2 利用HB-AFT结合同伦弧长延拓对系统解进行追踪为采用HB-AFT方法追踪系统的全局周期轨线, 进而研究系统的滞后突跳行为, 需要将HB-AFT过程中的Newton-Raphson迭代过程与弧长延拓法相结合.为此将HB-AFT方法求解变柔度振动周期响应的谐波平衡方程H(P, Q)=0的问题转化为
|
(20) |
采用HB-AFT法，结合弧长延拓法，并采用Newton-Raphpson迭代校正方法追踪系统的过程如下[9]:
① 利用HB-AFT方法得到相应转速对应的精确解P, 给定系统一组初始解(λ, P)0;
② 根据延拓方向ν0, 设定一个较小的迭代步长δ(本文初始迭代步长δ=0.5δ), 得到初始预估解(λ, P)1;
③ 对(λ, P)1通过Newton-Raphson进行迭代, 直至收敛至容许误差范围ε内, 得到新的精确解(λ, P)*;
④ 计算(λ, P)1与(λ, P)*形成的切向量角度aν, 若aν < 6°, 则取步长δ=2δ, 以(λ, P)*为初值，利用预估-校正方法继续进行下一步计算; 若aν>18°, 则放弃(λ, P)*, 取步长δ=0.5δ, 重新进行计算.
⑤ 重复步骤②, ③和④, 直至λ达到设定值, 得到解(λ, P)*, 输出P*.
对于HB-AFT结合同伦弧长延拓法并采用Newton-Raphson迭代校正方法追踪系统分支过程如图 3所示.
|
图 3 预估-校正追踪系统解的过程 Fig.3 Process of tracking system solution with predictor and corrector |
假设球轴承有j个滚动体, 第j个滚动体的位置角为θj, 则
|
(21) |
式中: 
第j个滚动体与滚道的接触变形量为
|
(22) |
式中:x, y分别为轴承在X, Y两个方向的位移量; δ0为轴承游隙.
滚道上产生的非线性接触变形载荷为
|
(23) |
式中“+”表示:当δj>0时, 滚道上的非线性接触变形为δj; 当δj≤0时, 滚道上的非线性接触变形为0, 即滚动体与环道不接触, 没有接触变形载荷.Kb为滚动体与滚道的接触变形系数, 则
|
(24) |
式中Ki, Ko分别为滚动体与内、外圈滚道之间接触点的刚度.
所以, 球轴承产生的轴承载荷为
|
(25) |
在HB-AFT计算过程中要求对谐波项数K进行选取, 如果项数太少, 将导致计算结果不准确；如果项数过多, 计算效率就会大幅降低.图 4为AFT算法与数值算法的对比示意图.
|
图 4 谐波项数与离散点的选取对比图 Fig.4 Comparison of harmonic frequency and discrete point selection (a)—K=9, N=150; (b)—K=9, N=30; (c)—K=18, N=80; (d)—K=18, N=150. |
如图 4所示, 虚线代表AFT算法, 实线代表数值算法.图 4d为谐波项数K=18时, AFT方法与数值计算方法结果的比较, 可以发现, 相对图 4a, 吻合度大幅提高.
如图 4b所示, 在选取离散点数N值时, 如果N取得过小, 在DFT过程中将会发生频率混叠现象, 计算结果将严重失真；选取过大将增大计算量.通过对比图 4a, 4b, 4c和4d, 本算例取谐波项数K=18, 时域离散点数N=150.
2.4.2 计算滞后突跳区间滚动轴承-转子系统相关参数如表 1所示.
|
表 1 滚动轴承-转子系统参数 Table 1 Parameters of the rolling bearing rotor system |
以保持架转速Ω为控制参数, 采用嵌入弧长延拓的HB-AFT方法追踪轴承-转子系统的幅频响应曲线, X方向代表转子系统的竖直方向, Y方向代表转子系统的水平方向, 计算结果如图 5所示.
|
图 5 幅频响应曲线 Fig.5 Frequency response curve |
图 5中, 虚线部分代表不稳定区间段转子系统在X, Y方向可能存在双稳态特性.系统在X方向的239.2~264.7 rad/s区间出现了向左偏移的共振峰值, 即表现为软的滞后特性.可能由于X方向受重力影响导致只在该方向上具有滞后特性.
系统在Y方向的198.8~210 rad/s区间出现了交叉结构, 即系统在Y方向上具有软、硬滞后特性.在X, Y方向378~422 rad/s区间表现为不稳定的幅频响应曲线.如图 5a所示，伴随着控制参数Ω的变化, 系统X方向上的幅频响应曲线从A点跳跃到A′上, 此时的振幅相较原来大幅提升, 称点A是系统的一个分岔点.随着控制参数Ω的减小, 在转接点B处产生的亚谐振动使C点向下跳到C′点, 此时振幅相较原来大幅下降.
综上, 系统在X方向上的参数激励接触共振频率区间产生了软的滞后突跳现象, 在Y方向上产生了软硬共存的滞后突跳现象.滞后分岔的双稳态乃至多稳态行为, 会给系统带来突跳冲击作用, 冲击作用是裂纹产生的主要因素, 而裂纹衍生是产生严重断裂故障的重要机制[14], 这对于转子系统来说是不希望发生的.
2.4.3 利用Floquent理论进行稳定性判定Floquent理论可用于分析非线性动力学周期响应的局部稳定性, 当系统设计参数集发生改变时, 根据Floquent特征乘子模可判断非线性动力系统是否失稳以及周期解的分岔形式.本文采用Hsu[15]方法求得Floquent特征矩阵的全部特征值, 并且取最大值作为Floquent乘子.利用Floquent乘子判断稳定性：当Floquent乘子模的最大值小于1时稳定; 大于1不稳定; 等于1临界稳定.可根据模最大的FIoquet乘子穿出复平面上单位圆时的情况将失稳分岔方式分类[14].
如表 2所示, 当转子系统的解逆时针绕过点A时, 发现特征矩阵的Floquent乘子λ在Ω=264.73 rad/s处从实部+1处穿出单位圆, 即转子系统发生鞍节型分岔进入了非稳定区间.
|
|
表 2 分岔点A附近的Floquent乘子λ Table 2 Floquent multipliers λ around the branch point A |
将上述算例中转子系统的质量设为20 kg, 阻尼系数设为200 N·s/m, 轴承游隙为10 μm, 计算滚动轴承赫兹接触刚度分别为7.5×109, 1×1010, 1.25×1010和1.5×1010N/m时整个转子系统的幅频响应曲线.为更直观地观察系统滞后特性曲线, 选取X, Y方向上滞后特性比较明显的170~250 rad/s响应区间进行研究.
如图 6所示, 当轴承赫兹接触刚度由7.5×109N/m变化到1×1010N/m时, X方向上的滞后突跳区间由208~223 rad/s变化到225~243 rad/s，且伴随着滞后区间的增大；当轴承赫兹接触刚度增加到1.25×1010N/m时, 滞后特性消失.当轴承赫兹接触刚度增加到1.5×1010N/m时同样没有发现滞后现象的发生.当轴承赫兹接触刚度为7.5×109N/m时Y方向上没有滞后突跳现象发生, 当轴承赫兹接触刚度增加到1×1010N/m时开始发生滞后突跳且滞后突跳区间逐渐右移, 振幅随之逐渐减小.研究整个变化过程得出结论:随着轴承赫兹接触刚度的增加, 转子系统的滞后区间呈现右移的趋势；在X方向上随刚度增大滞后区间逐渐减小直至消失, 在Y方向上刚度增加到一定程度会产生滞后突跳现象.
|
图 6 不同赫兹接触刚度的幅频响应曲线 Fig.6 Frequency response curves with different Herz contact stiffness |
计算阻尼系数分别为200, 800, 1 200，2 500 N·s/m时转子系统的幅频响应曲线.在研究轴承-转子系统的振动特性受阻尼系数的影响时, 保持其他参数不变.
如图 7所示, 阻尼系数在200 N·s/m经由800, 1 200 N·s/m变化到2 500 N·s/m的过程中, 阻尼系数主要影响转子系统滞后突跳区间的幅值, 转子系统在X, Y方向上的振幅随着阻尼系数的增加而减小, 且当阻尼系数大于某一数值时, 转子系统在X, Y方向的滞后效应将会消失.
|
图 7 不同阻尼系数的幅频响应曲线 Fig.7 Frequency response curves with different damping coeffcients |
计算轴承游隙分别为10, 20, 30, 50, 60 μm时转子系统的幅频响应曲线, 在研究轴承游隙对系统滞后特性的影响时, 保持其他参数不变.
图 8为轴承游隙分别为10, 20, 30, 50, 60 μm时转子系统在X, Y方向上的幅频响应曲线，其中断点部分为失稳区间，连续曲线部分为稳定运行区间.从Y方向的幅频响应曲线可观察到轴承游隙对整个转子系统滞后突跳区间的影响.转子系统Y方向的幅频响应曲线出现了交叉结构, 一般意味着系统此时伴随着软硬滞后特性并存的现象[16].同时, 随着轴承游隙的增加, 转子系统的滞后运动区间明显右移且滞后区间也逐渐变大, 并且滞后区间的幅值明显增大, 对系统的冲击作用更大.当游隙减小到一定程度时滞后特性将会消失, 转子系统的振幅将会减小.
|
图 8 不同轴承游隙的幅频响应曲线 Fig.8 Frequency response curves with different bearing clearance |
1) 受滚动轴承非线性因素的影响, 转子系统可能会呈现滞后突跳的特性.
2) 滚动轴承的游隙与等效刚度会影响转子系统滞后突跳区间的范围, 同时也会影响滞后突跳区间的振幅, 所以轴承游隙与轴承等效刚度对转子系统的滞后突跳的特性具有重要的影响.
3) 支承阻尼系数相对于轴承游隙、轴承等效刚度而言对转子系统滞后突跳行为的影响较小.在系统临界转速附近区域, 阻尼系数几乎不影响转子系统滞后突跳区间的范围, 但会降低滞后突跳区间的振幅, 故在转子系统的参数设计时, 适当增加阻尼系数对转子系统有一定的积极作用.
| [1] |
Yamamoto T, Ishida Y, Ikeda T, et al.
Nonstationary vibration of a rotating shaft with nonlinear spring characteristics during acceleration through a critical speed[J]. Nonlinear Dynamics, 1990, 1(5): 341–358.
DOI:10.1007/BF01893168 |
| [2] |
Tiwari M, Gupta K, Prakash O.
Effect of radial internal clearance of a ball bearing on the dynamics of a balanced horizontal rotor[J]. Journal of Sound & Vibration, 2000, 238(5): 723–756.
|
| [3] |
Villa C, Sinou J J, Thouverez F.
Stability and vibration analysis of a complex flexible rotor bearing system[J]. Communications in Nonlinear Science & Numerical Simulation, 2008, 13(4): 804–821.
|
| [4] |
Peletan L, Baguet S, Richardet G J, et al.Use and limitations of the harmonic balance method for rub-impact phenomena in rotor-stator dynamics[C]// ASME Turbo Expo 2012: Turbine Technical Conference and Exposition.[S.l.], 2012: 647-655.
https://www.researchgate.net/publication/267503726_Use_and_Limitations_of_the_Harmonic_Balance_Method_for_Rub-Impact_Phenomena_in_Rotor-Stator_Dynamics |
| [5] |
Sun C, Chen Y, Hou L.
Steady-state response characteristics of a dual-rotor system induced by rub-impact[J]. Nonlinear Dynamics, 2016, 86(1): 1–15.
|
| [6] |
Sundararajan P, Noah S T.
Dynamics of forced nonlinear systems using shooting/arc-length continuation method—application to rotor systems[J]. Journal of Vibration and Acoustics, 1997, 119(1): 9–20.
DOI:10.1115/1.2889694 |
| [7] |
Ganesan R.
Effects of bearing and shaft asymmetries on the instability of rotors operating at near-critical speeds[J]. Mechanism & Machine Theory, 2000, 35(5): 737–752.
|
| [8] |
Wang W M, Li Q H, Gao J J, et al.
An identification method for damping ratio in rotor systems[J]. Mechanical Systems & Signal Processing, 2016, 68/69: 536–554.
|
| [9] |
张智勇.球轴承-转子系统变柔度振动的分岔与滞后行为[D].哈尔滨: 哈尔滨工业大学, 2015.
( Zhang Zhi-yong.Bifurcations and hysteresis of varying compliance vibrations of a ball bearing-rotor system[D].Harbin: Harbin Institute of Technology, 2015. http://cdmd.cnki.com.cn/Article/CDMD-10213-1015957637.htm ) |
| [10] |
白长青, 许庆余, 张小龙.
考虑径向内间隙的滚动轴承平衡转子系统的非线性动力稳定性[J]. 应用数学和力学, 2006, 27(2): 159–169.
( Bai Chang-qing, Xu Qing-yu, Zhang Xiao-long. Nonlinear stability of balanced rotor due to the effect of ball bearing internal clearance[J]. Applied Mathematics and Mechanics, 2006, 27(2): 159–169. DOI:10.3321/j.issn:1000-0887.2006.02.005 ) |
| [11] |
Ma H, Lu Y, Wu Z, et al.
A new dynamic model of rotor-blade systems[J]. Journal of Sound & Vibration, 2015, 357: 168–194.
|
| [12] |
Han Q, Chu F.
Parametric instability of flexible rotor-bearing system under time-periodic base angular motions[J]. Applied Mathematical Modelling, 2015, 39(15): 4511–4522.
DOI:10.1016/j.apm.2014.10.064 |
| [13] |
Yamauchi S.
The nonlinear vibration of flexible rotors, 1st report, development of a new analysis technique[J]. Transactions of the Japan Society of Mechanical Engineers C, 1983, 49(446): 1862–1868.
DOI:10.1299/kikaic.49.1862 |
| [14] |
Nayfeh A H, Balachandran B.
Applied nonlinear dynamics[M]. Weinheim: Wiley-VCH, 2004: 187-207, 300-331, 449-450.
|
| [15] |
Hsu C S.
On approximating a general linear periodic system[J]. Journal of Mathematical Analysis & Applications, 1974, 45(1): 234–251.
|
| [16] |
Yao M H, Zhang W.
Using the extended Melnikov method to study multi-pulse chaotic motions of a rectangular thin plate[J]. International Journal of Dynamics and Control, 2014, 2(3): 365–385.
DOI:10.1007/s40435-013-0031-z |