基于嵌入式谱随机有限元法的转子系统随机不平衡响应特性分析

周生通， 李鸿光， 张 龙， 周新建

(1.华东交通大学 机电与车辆工程学院，南昌 330013；2.上海交通大学 机械系统与振动国家重点实验室，上海 200240)



基于嵌入式谱随机有限元法的转子系统随机不平衡响应特性分析

周生通1， 李鸿光2， 张 龙1， 周新建1

(1.华东交通大学 机电与车辆工程学院，南昌 330013；2.上海交通大学 机械系统与振动国家重点实验室，上海 200240)

原始不平衡量在转轴中沿轴线呈任意空间分布，因而实际转子的不平衡量需用分布不平衡量和点不平衡量两种方式共同表达。考虑两种不平衡量的随机性，以随机场和随机变量的方式模拟转子中存在的两类不平衡量，并采用嵌入式谱随机有限元方法建立随机不平衡量下的转子系统不平衡响应随机分析模型。结果显示，提供的随机表征方法能较好地模拟两类转子不平衡量的随机性；不平衡量的随机性影响着转子不平衡响应的变异性，可以发现：① 共振频率附近的位移响应幅值标准差也较大；② 位移响应幅值的变异系数受低频段内共振频率的影响较小，但随着共振频率阶次的增高影响也逐渐增大；③ 位移响应相位角的标准差变化情况与位移幅值的变异系数变化情况是基本类似的。

质量不平衡；谱随机有限元方法；随机不平衡响应

一个转子在设计上一般都使它相对于旋转轴线是轴对称的。但是由于工艺上的一系列因素，最后装配完毕的转子总是不能做到动力上的完全轴对称，也就是存在一定的不平衡量，这种不平衡量通常称之为原始不平衡量[1]。造成原始不平衡量的因素主要有：转子材质的不均匀性，键槽不对称引起的不平衡，转子加工中总是产生一些圆度偏差和偏心等。对于汽轮机转子尚有各个叶片之间的差别，叶片锁口及末叶片的不对称等不平衡因素。对于压缩机尚有叶轮的不平衡量。所有这些因素造成的不平衡量都属于随机性质的[1]。

像叶片和叶轮位置处的不平衡量是转子系统不平衡量的主要贡献源。在理论分析时，一般将叶片和叶轮简化为刚性圆盘，相应的不平衡量则以力的形式施加在圆盘的几何中心处，用不平衡质量、偏心距和相位角三个量确定不平衡的大小和方向。不过，如果将连续转轴分割为微小厚度的圆片，那么正如单个圆盘一样，每一个圆片上都存在一定的不平衡量，且大小和方向各不相同。也就是说，转轴中的不平衡量分布函数实为一条任意的、随机的空间分布曲线[1-2]。虽然不平衡量的空间任意分布事实很早就被人们认识到，但直到1993年才被LEE等[3]在传递矩阵方法中予以考虑。后来，SHIH等[4]又提出了识别柔性转子不平衡量分布情况的方法。该方法将分段连续的质量偏心曲线用傅里叶级数展开，但这种方法不太适用于转子动力学的有限元分析程序，为此YANG等[5]提出采用多项式曲线模拟转子不平衡质量偏心曲线，以便借助有限元程序更容易地识别转子中存在的不平衡量分布。DEEPTHIKUMAR等[6]在此基础上进一步发展了同时存在分布不平衡和轴弯曲故障的柔性转子模态平衡方法。不过，上述这些研究都是从确定性的角度处理转子不平衡量的。

现有对转子不平衡量随机性的考虑多是针对圆盘不平衡量的[7-8]，即将圆盘的偏心距、质量或相位角作为随机变量处理。不过，STOCKI等[9]在转子振动离散性分析中则考虑了分布不平衡量的不确定性，其采用的思路是将分布不平衡量表示为转子前四阶振型的加权叠加，权值设为随机变量，并使用非嵌入式的谱方法开展了相应的随机分析工作。本文同样从随机角度考虑转子系统中的不平衡量，并以随机场和随机变量的形式将它们表达为两种形式，即集中不平衡量(Lumped Unbalance)和分布不平衡量(Distributed Unbalance)。最后，对某一50 MW汽轮机转子同时计及这两种不平衡量，并基于嵌入式谱随机有限元方法建立该转子的随机不确定性分析模型，实现其不平衡响应的随机分析。

1 转子不平衡量的随机表征

1.1 分布不平衡量

图1 转子中的不平衡量Fig.1 Rotor unbalances

将分布不平衡量在单个转子单元内产生的动能代入拉格朗日方程，即可得到对应的不平衡力单元向量：

(1)

(2)

为了处理积分中的随机量，将一维随机场u(x,ω)和φ(x,ω)分别作Karhunen-Loeve展开并分别保留至L1阶和L2阶，有：

(3)

(4)

(5)

式(5)由于包含了连乘的指数函数，进一步处理会比较繁琐。这里将相位角随机场φ(x,ω)退化为正态随机变量，即取：

(6)

这样，式(5)就可简化为：

(7)

进一步，取指数函数的级数展开近似表达：

(8)

则由式(7)表示的转子空间随机不平衡量产生的单元力向量就可整理为如下形式：

(9)

式中:Qkj是维数为8×1的列矩阵，表达式为：

(10)

更特别地，若不平衡量分布密度u(x,ω)也同样退化为正态随机变量，即：

(11)

则有：

(12)

式中:

(13)

1.2 集中不平衡量

对于集中不平衡量，只需用随机变量就可完全表征不平衡量大小和相位角的随机性。在文献[8]中已经给出了当两者都是正态随机变量但相位角均值为零时的不平衡量表征公式。这里则直接将其扩展为不平衡量大小和相位分别是具有任意分布参数的正态随机变量的情况，对应的集中不平衡量表征公式为：

(14)

式中:Qkj是维数为4×1的列矩阵，表达式为：

(15)

2 不平衡响应的嵌入式谱随机有限元方程

谱随机有限元方程是依据多项式混沌基的正交性质建立的。从转子随机不平衡量的表征公式(9)、(12)和(14)中可以看到，作为右端项的不平衡力Fu是以多项式的形式表达的，无法直接利用多项式混沌基的正交性。需要将Fu的表达式改写成以多项式混沌展开的形式，即：

(16)

将原来Fu的表达式改写为式(16)的形式后，即可建立起转子系统不平衡响应的谱随机有限元方程：

(17)

式中:Z(ω)是转子系统的动刚度矩阵。在不考虑其他随机因素的情况下，式(17)最终可转化为如下P个线性方程组：

Z(ω)δi=Fi,i=0,1,…,P-1

(18)

解方程组即可求得不平衡响应的多项式混沌展开中的广义坐标值δi,i=0,…,P-1，从而表达出随机不平衡响应Δ，即为：

(19)

(20)

3 数值算例

求解某一50 MW汽轮机转子动平衡系统(如图2所示)的随机不平衡响应。设汽轮机转轴中的不平衡量分布密度的大小u(x,ω)服从高斯随机场且具有指数协方差函数，相应的相位角φ服从正态随机量；同时，取作用在圆盘10和圆盘19上的不平衡量U10和U19的大小和相位角亦均服从正态分布。具体的随机量设定值如表1所示。

图2 某50 MW汽轮机转子动平衡系统Fig.2 A 50 MW steam turbine rotor dynamic balancing system

元件属性名称分布类型均值标准差备注圆盘10不平衡量U10正态分布0.190.038相位角φ10正态分布0π/15圆盘19不平衡量U19正态分布0.30.06相位角φ19正态分布π/3π/15转轴相位角φ正态分布π/6π/60分布不平衡量大小u(x)一维高斯随机场0.010.001设转轴分布不平衡量随机场具有指数协方差函数,相关长度6.36m,区间为[0,6.36]

该模型中，单个集中不平衡量采用2维5阶精度描述，分布不平衡量采用5维5阶精度描述，故整个谱随机有限元模型将形成一个具有9维5阶的多项式混沌展开模型，共计P=2 002项多项式混沌基函数。进而，依据式(18)～(20)即可得到转子上任意节点位移的随机响应情况。限于篇幅，这里仅给出圆盘10位置处的位移响应概率统计信息，如图3～7所示。图3(a)和(b)所示分别为圆盘10处的水平和竖直位移频响函数幅值的均值、上下包络线和标准差曲线；图4给出了对应水平和竖直位移频响幅值的变异系数随激励频率的变化情况；图5和图6则类似地给出了水平和竖直位移相位角的统计信息；图7(a)和(b)给出了在特定转速频率(25.5 Hz和40.0 Hz)下水平位移的直方图和基于核密度估计的概率密度信息，其中频率25.5 Hz靠近转子系统一阶临界转速，频率40.0 Hz则远离共振频率。

其中，位移均值(实线)、标准差(实线)、上包络(点线)、下包络(虚线)和1 000个Monte Carlo样本(灰线)图3 圆盘10的位移幅值均值和标准差信息Fig.3 Mean and standard deviation of displacementamplitudes at disk 10

图4 圆盘10处水平和竖直位移幅值的变异系数曲线Fig.4 The coefficient of variation curves of horizontal and vertical displacement amplitudes at disk 10

从图3(a)和(b)中可以看到，受到转子随机不平衡量的影响，圆盘10处的水平和竖直位移频响幅值曲线表现出了离散特性，并围绕着均值曲线上下随机波动，而偏离均值的大小与激励频率有关，当激励频率与系统固有频率一致，即发生共振时，偏离量最大，即在共振处出现有较大的标准差峰值。但此处的标准差大小并不能准确反映位移频响幅值随机离散的程度，图4则采用变异系数表达离散程度以消除均值的影响。从图4中可以看到，水平和竖直位移频响函数幅值的变异系数在小于40 Hz的频段内基本没有变化，而这一频段内就包含有转子系统的第一阶临界转速频率，且此临界频率处的响应幅值标准差就比较大。随着激励频率的增加，当达到其他高阶临界转速时，变异系数就开始有较大的浮动，可以看到阶数越高变异系数峰值也就越大。这说明转子的不平衡随机因素对系统高阶共振频率处位移幅值的离散程度影响较大，而对低阶影响则较小。如在第一阶临界转速频率处的影响几乎可以忽略，因为它与附近非共振频率处的离散程度几乎一致。

其中，位移相位角均值(实线)、上包络(点线)、下包络(虚线)和1000个Monte Carlo样本(灰线)图5 圆盘10的位移相位角信息Fig.5 Phase angles of displacement at disk 10

图5(a)和(b)则分别给出了随机不平衡量影响下的水平和竖直位移相位角的离散情况，图6(a)和(b)中显示的是对应相位角的标准差和变异系数情况。可以看到，与位移幅值离散特性表达不同，位移相位角的离散程度用标准差表征更好，且此时的水平和竖直位移相位角标准差的变化情况与对应的水平和竖直位移幅值变异系数的变化情况类似，如图4所示。其中的原因主要是由于位移幅值均值及其标准差的数值大小在不同的频率范围内存在有量级上的较大差别，当表示离散特性时不同量级上的标准差之间不易反映整体离散程度的大小，此时就需要剔除均值的影响用变异系数来表示离散程度，但对于位移相位角来说，其标准差不存在量级差别问题，而均值在一定的范围内(不超过360°)变化，此时用标准差表示离散程度会更好，若采用变异系数则效果上不佳(如图6(a)和(b)虚线所示)。

图6 圆盘10的位移相位角标准差和变异系数频谱曲线Fig.6 The standard deviation and coefficient of variation curves of displacement phase angle at disk 10

图7 圆盘10的水平位移在不同转速频率下的均值、标准差、概率密度函数直方图和核密度估计曲线Fig.7 Means, standard deviations, histograms of probability density functions, kernel density estimation curves of the horizontal displacement at disk 10 under different rotating frequencies

进一步，从图7(a)和(b)中可以看到某一特定激励频率下水平位移频响幅值的概率分布比较接近正态分布，这与我们在例子中所采用的随机量类型均为正态型的有较大关系。其中，在第一阶临界转速频率附近(25.5 Hz处)时，其对应位移幅值的均值为1.214×10-3m、标准差为2.006×10-4m；在非共振区域(40 Hz处)时，均值为4.728×10-5m、标准差为7.960×10-6m。从量级上对比，可以看到两个频率处的标准差在数值上相差有两个数量级，但这也只能表明系统在共振频率处的响应离散幅值更大；而从幅值离散程度上对比，可以计算出两频率处的幅值变异系数分别为0.165 2和0.168 4，比较接近，说明实际中这两个频率处的随机离散程度差别并不大，与图4结果一致。

4 结 论

本文采用分布不平衡量和集中不平衡量两种概念描述转子中存在的不平衡量。在考虑不平衡量随机性的基础上，推导了分布不平衡量和集中不平衡量的力表达式，并利用嵌入式谱随机有限元方法建立了转子系统不平衡响应的随机分析模型，求解了某一50 MW汽轮机转子在动平衡摆架支撑系统上的随机不平衡响应概率信息。总结算例结果，可以看到：

(1) 受到转子随机不平衡量的影响，系统不平衡响应(位移幅值频响函数)表现出了较大的随机离散特性，尤其是在共振频率附近，其对应的响应幅值标准差要明显高于非共振频率区；

(2) 但从位移频响幅值的变异系数上分析，在低频段内各频率上的变异系数基本保持不变(即使是在第一阶临界转速频率处)，而随着激振频率的增大，变异系数的浮动也开始增大，尤其是在高阶次的共振频率处，变异系数出现了较大峰值，这说明了转子不平衡随机因素对更高阶共振频率处的位移幅值离散程度影响更为明显；

(3) 对于位移相位角来说，其也表现出了较大的随机离散特性，但其离散程度更适合用相位角标准差来表征。

[1] 周仁睦. 转子动平衡——原理、方法和标准[M]. 北京: 化学工业出版社, 1992.

[2] 王悦武，田杜平，徐锡林. 高速平衡技术装备[M]. 上海: 上海交通大学出版社, 2013.[3] LEE A, SHIH Y, KANG Y. The analysis of linear rotor-bearing systems: a general transfer matrix method[J]. Journal of Vibration and Acoustics,1993, 115(4): 490-497.

[4] SHIH Y, LEE A. Identification of the unbalance distribution in flexible rotors[J]. International Journal of Mechanical Sciences,1997, 39(7): 841-857.

[5] YANG T, LIN C. Estimation of distributed unbalance of rotors[J]. Journal of Engineering for Gas Turbines and Power,2002,124(4): 976-983.

[6] DEEPTHIKUMAR M B, SEKHAR A S, SRIKANTHAN M R. Modal balancing of flexible rotors with bow and distributed unbalance[J]. Journal of Sound and Vibration,2013,332(24):6216-6233.

[7] ZHANG Y, WEN B, LEUNG A Y T. Reliability analysis for rotor rubbing[J]. Journal of Vibration and Acoustics,2001, 124(1): 58-62.

[8] DIDIER J, FAVERJON B, SINOU J. Analysing the dynamic response of a rotor system under uncertain parameters by polynomial chaos expansion[J]. Journal of Vibration and Control,2012, 18(5): 712-732.

[9] STOCKI R, LASOTA R, TAUZOWSKI P, et al. Scatter assessment of rotating system vibrations due to uncertain residual unbalances and bearing properties[J]. Computer Assisted Methods in Engineering and Science,2012, 19(2): 95-120.

Stochastic unbalance response characteristics of rotor systems based on intrusive spectral stochastic finite element method

ZHOU Shengtong1, LI Hongguang2, ZHANG Long1, ZHOU XinJian1

(1. School of Mechatronics & Vehicle Engineering, East China Jiaotong University, Nanchang 330013, China;2. State Key Lab of Mechanical Systems & Vibration, Shanghai Jiao Tong University, Shanghai 200240, China)

Mass unbalance in a flexible shaft is distributed arbitrarily along its axis, so distributed mass unbalance and lumped mass unbalance should be considered and presented, respectively for mass unbalance simulation of real rotor systems. Taking the random nature of these two unbalances in a flexible shaft into account, they are represented as random fields and random variables. Here, stochastic analysis of unbalance response of rotor systems with random mass unbalances were performed using the intrusive spectral stochastic finite element method. Numerical results illustrated that the presented stochastic representation methods of distributed and lumped mass unblances are suitable to simulate the random nature of mass unbalances. Results showed that the variations of unblance responses are largely affected by the randomness of unbalances, for example, ① the standard devitations of displacement amplitudes near the resonance frequencies become larger; ② the effects of resonance frequences with in a lower frequency band on the variation coefficient of displacement amplitudes are less, but the effects increase with increase in resonance frequency order; ③ the change curve of the standard deviation of displacement phase angle is similar to that of the variation coefficient of displacement amplitudes.

mass unbalance; spectral stochastic finite element method; stochastic unbalance response

国家自然科学基金(51505146)；江西省自然科学基金(20161BAB216135;20122BAB206027)

2015-07-02 修改稿收到日期：2015-09-27

周生通 男，博士，讲师，1984年生

李鸿光 男，博导，教授，1972年生

E-mail: hgli@sjtu.edu.cn;zhoust@ecjtu.edu.cn

TH113

A

10.13465/j.cnki.jvs.2016.19.008