一种基于频响因子的转子系统基础冲击响应计算方法

谢最伟，贺少华，吴新跃

(海军工程大学 船舶与动力学院机械工程系，武汉 430033)

对于动力学系统响应的求解即动力学微分方程的求解，一般有三种解法，解析解法、近似解析解法和数值解法，只有极少数动力学方程存在解析解，所以大部分动力学方程只能采用近似解析解法和数值解法，近似解析方法主要针对弱非线性系统的稳态响应的求解，如摄动法、谐波平衡法、多尺度方法、平均法等，而对于数值解法，则适应于任何的系统响应求解，主要包括显式算法中的中心差分法、Runge-Kutta法和隐式算法中的 Wilson-θ法、Houbolt法、Newmark 法等［1－4]。

动力学系统(结构)外界冲击本文的定义:冲击是一种振动系统运动量的瞬时突然变化，其动能传输是在一个与振动系统固有周期相比的短时间内进行的，冲击因其时间历程比较复杂而不同于撞击。以舰载设备的水下非接触爆炸冲击为例，冲击理解为对舰艇结合部(外壳、甲板、舱壁等)的一种非接触的相当远距离的爆炸(当水下爆炸时，气泡不触及外壳结构或受压物体结构)作用，其结果为在舰船正常运行中不会出现的一种振动系统“舰艇”运动量的巨大干扰。距离很近的水下爆炸或空中爆炸所引起的作用，另作特别处理，所以，结构冲击理解为:在承压的舰艇结合部中和在设备基座中的机械干扰。冲击载荷理解为由于爆炸产生的机械载荷。

旋转机械的瞬态冲击响应分析区别于一般的静态条件下结构的相应分析，主要体现在以下两个方面:①轴承－转子的交互效应;② 陀螺效应。除了以上两个因素，当系统基础遭受旋转运动激励时Coriolis效应将不可忽略。在过去的几十年中，国内外学者对于旋转机械(转子－轴承系统)在地震激励下的响应进行了研究［5－9]，取得了一些成果，如:旋转激励成分对响应的影响;线性理论的应用局限性;数值算法的改进等等。

考查地震或水下非接触爆炸冲击下旋转机械的动态响应特性，一般从研究转子系统基础冲击响应出发。对于转子－轴承系统的冲击响应研究，主要方法为:首先对转子－轴承系统进行建模，用基础运动(位移、速度和加速度)来模拟基础冲击，建立系统力学方程;然后通过数值计算求解方程得到响应。通过响应分析检验以下几个方面的内容:①润滑油膜在任何时候必须保持一个最小的油膜厚度以避免转子和轴承不会发生碰摩;②轴承基础能够承受冲击轴承反力;③转子冲击位移(变形)或应力在安全的范围之内。由于陀螺效应和轴承的交互效应，转子系统运动方程系数矩阵呈非对称性，不能在模态坐标下解耦，无法利用常规模态叠加法求解，所以已往的研究一般采用数值积分如Newmark法等进行迭代求解，随着有限元和计算机技术水平的飞速发展，目前，借助嵌入了这些数值积分算法的商业有限元软件进行冲击响应计算分析成为了国内的主流。但是，必须指出的是，绝大部分有限元商业软件的知识产权属于国外，过分依赖这些商业软件将使我国本已滞后的机械设备抗冲研究特别是理论研究更加落后国外。同时，数值积分法相对线性叠加法要耗费更多的计算资源和计算时间。

转子系统复模态分析法相对实模态分析法能够更好地表达转子的运动规律，因为复模态包含了转子的旋转方向和振动幅值双重信息［10]，本文正是从这一思想出发，将提出并实现一种复数域内转子系统冲击响应计算方法，无需坐标解耦但仍可以利用线性叠加法进行响应求解，从而克服常规模态叠加法和传统数值积分算法的上述缺点，为转子系统提供一种普适的冲击瞬态响应线性叠加计算方法。

1 理论推导

1.1 复数域内转子系统运动微分方程的推导

一个基础受加速度运动激励的n维线性转子系统(轴线沿x方向)的运动方程可以写为:

式中:Y、Z描述着基础的运动，y、z为系统相对基础的位移(变形)向量，它们均为n维。

令:

则:

其中:

将式(2)代入式(1)，等式两边左乘［T]－1:

进一步简写成:

1.2 特征方程的演变和响应求解

对于g(t)，有下面的傅里叶展开式

对于线性系统，同样有解的形式为:

式中的Gf、Gb、Pf和Pb均为n维向量。由上式我们可以发现，转子在YZ平面的任一频率运动被分解成了相同频率的一正一反两个圆周运动，也就是通常所说的正向和反向涡动(forward和backward rotating)。当两涡动运动的幅值相等即时，合成运动为直线运动，这对应于一般的非旋转振动系统的情况;当两涡动运动的幅值不相等时，合成运动为椭圆运动。所以，可以将一般的非旋转振动系统的直线运动看成涡动的一种特例。

将式(4)，式(5)代入式(3)，方程两边相等的唯一条件是方程两边相同频率ejωkt项前系数相等，相同频率e－jωkt前系数也相等。对于任一给定k，可以得到:

上面两式表达的信息是相同的，所以只需考查其中的任一个。

但是，对于转子系统，由于陀螺效应和轴承－转子的交互效应，使得方程(3)的系数矩阵一般不再对称，模态向量关于系数矩阵的正交性不再存在，无法将系统方程在模态坐标下解耦，模态叠加法求解系统响应失败。

将式(6)写成一般矩阵束的本征方程形式(形如(sA －B)x=0)，令 s=jωk:

或者:

可以证明上式中的矩阵束 s［A] －［B] 为简单(simple)矩阵束［11]，于是有:

式中，Li和 Ri分别为 λi［A] －［B] 的左右特征向量即:

λi为特征值，由=0决定，且特征向量满足条件:LTAR=I，LTBR=Λ，Λ为由特征值组成的对角矩阵。在本例中Ri和Li的形式为:

即为系统的模态向量。故式(7)可进一步演变为:

当为各向同性转子时，A12=A21=0。在本文中，将A11、A12、A21和 A22称为“频响因子(frequency response coefficient)”。

如上所述，得到各频率成分响应后，将响应合成即可得到系统总的响应。在实际中，基于这样的常识:当结构受冲击作用时，在某种意义上结构的响应本质上是一种振荡，响应时主要出现的频率是冲击的优势频率和结构的固有频率。所以，只需考查结构固有频率和冲击优势频率处的响应值。

1.3 响应合成

在实际机械抗冲考核中，我们需要计算其在标准冲击激励下各部分的最大瞬态响应，并不需要整个时间历程的响应。通过上小节求得各频率处的冲击响应，需要对各频响幅值进行组合。组合方法可以仿照美国海军 DDAM(dynamic design analysis method)［12]，使用下列三种方法之一:

(1)绝对值求和(ABS):

本文公式:

(2)平方和的均方根(SRSS)

本文公式:

(3)美国海军研究实验室求和(NRL)

本文公式:

式中:Pf(ωj)和Pb(ωj)为任一频率正反涡动幅值。

2 实例分析

2.1 某刚性转子系统模型和参数

一柔性支撑刚性转子系统如图1所示，系统的详细参数如表1所示，该系统的运动状态完全由旋转轴两端与支撑连接点的运动状态决定。

图1 刚性转子系统的坐标系和支撑Fig.1 Coordinate system and supports of a rigid rotor system

表1 系统参数Tab.1 System parameters

2.2 运动微分方程的建立

建立如式(1)的系统相对基础的运动方程式:

式中:

{FY1FY2FZ1为支撑轴承等效的弹簧－阻尼单元支撑力。在上式中，与式(1)不同的地方在于它还考虑了基础的速度激励，由于地震冲击或是舰载设备水下非接触爆炸冲击一般只采用加速度形式来表示，所以在这里忽略速度激励项。

进一步的，令:

则:

式中:

正如上文所述，方程中速度向量前的矩阵加进了陀螺效应和轴承等效弹簧－阻尼单元的贡献，G为反对称矩阵(skew-symmetric)(如果a=b)，而如果轴承的交叉阻尼项又不相等，这使得总体阻尼效应矩阵变得十分复杂，不再正对称，这与一般非旋转振动系统的阻尼效应矩阵总被处理成对称的(如比例粘性阻尼)是不同的，刚度效应矩阵由于轴承等效弹簧－阻尼单元的贡献一般也不再对称(如果轴承交叉刚度项不相等)。

将方程(10)转换成式(3)的形式:

式中:

后续公式的推导过程与上一小节相同，不再赘述。

2.3 计算

最终计算得到的系统特征值λ和模态向量{Pf1，i=1，2，3，4，如表 2 所示。

表2 特征值和特征向量Tab.2 Eigenvalues and eigenvectors

同时，对应于式(8)的算术表达式为:

按照本文方法计算在图3所示加速度冲击激励下系统的响应，首先，按照上文所示方法将激励信号经傅里叶变换成图3右所示频率－幅值曲线，很显然，在本例中，图3的幅值(用A表示)与式(4)中Gf、Gb的关系为A=G=G。根据式(12)、图2和图3，在系统固

fb有频率 λ1= λ2=123.2 rad/s，λ3=416.8 rad/s，λ4=456.9 rad/s和冲击优势频率 λ5=9.8 rad/s处，选取响应合成方法为上文中的绝对值求和(ABS)，最终计算得到的最大位移响应为4.30 m。

然后，采用经典Newmark数值积分方法在时域直接求解系统响应，由于Newmark数值积分方法无条件稳定且具有二阶精度，在冲击和瞬态响应计算中最为常用［13]。Newmark数值积分法计算得到的最大位移值为3.56 m。

2.4 误差分析及说明

上小节计算结果表明，直接数值积分方法计算结果要明显小于本文方法计算结果，误差的主要原因在于本文方法类似于谱分析方法如DDAM，都是一种典型的保守计算方法，在响应合成时没有考虑各模态响应的相位差，进行的是简单的“绝对值”相加。文献［13] 指出，即使为非密集模态情况，两种合成方法(NRL和SRSS)的结果与时域合成结果相比，在特定条件下也能产生40%的相对误差，但这并不会妨碍其在工程中的应用，仍具有工程应用性，如DDAM，因为在实际中，设备的抗冲击设计本身就需要一个安全裕值。本文提出的方法的误差在理想的范围之内，所以可用于工程实际。

3 结论

由于陀螺效应和转子－轴承的交互效应，旋转机械的冲击瞬态响应计算比一般系统要复杂得多，以往采用的数值积分算法虽能解决任意复杂系统的响应求解问题，但这些算法往往需要耗费大量的计算资源和计算时间，为突破这一瓶颈，本文提出了一种复数域内的转子系统基础冲击响应计算方法，由于采用了线性叠加原理，该方法较好地克服了传统数值积分方法的上述缺点。同时，由于转子系统运动方程系数矩阵呈非对称性，不能在模态坐标下解耦，无法利用常规模态叠加法求解，所以，本文提出的方法可以视为一种改进的瞬态响应线性叠加计算方法。首先将激励和响应傅里叶展开成复数形式，包括正向旋转项和反向旋转项，根据方程左右两边相同频率前系数相等的事实得到特征方程，将特征方程写成简单矩阵束的本征方程形式，求得矩阵束的本征值和本征向量，将本征向量正规化，进一步得到矩阵束的逆阵，逆阵元素即为“频响因子”，将逆阵与激励相乘即可得到频率响应幅值，将所有频率响应成分叠加即可得到系统响应。通过一个工程实例，比较了所提方法与Newmark数值积分方法的计算结果，误差在合理的范围内，从而验证了其在工程应用中的合法性。本文所提方法可以作为转子系统基础冲击响应和瞬态响应计算的一种普适方法。“频响因子”在本文中为首次提出。

［1] 汪 玉，赵建华，杜俭业，等.基于多体动力学有限元计算的一种舰用柴油机抗冲击性能仿真［J] .振动与冲击，2009，28(11):87 －90，129.

［2] 高浩鹏，黄映云，王雪山.船用柴油机运动机构动力学动态分析方法研究［J] .振动与冲击，2009，28(8):84 －87.

［3] Kim SJ，Cho J Y，Kim W D.From the trapezoidal rule to higher-order accurate and unconditionally stable timeintegration method for structural dynamics［J] .Computer Methods in Applied Mechanics and Engineering，1997(149):73－88.

［4] Fung T C.Unconditionally stable higher-order accurate collocation time-step integration algorithms for first-order equations［J] .Computer Methods in Applied Mechanics and Engineering，2000(190):1651 －1662.

［5] Lee A S，Kim B O，Kim Y C.A finite element transient response analysis method of a rotor-bearing system to base shock excitations using the state-space Newmark scheme and comparisons with experiments［J] .Journal of Sound and Vibration，2006(197):595－615.

［6] Hori Y，Kato T.Earthquake-induced instability of a rotor supported by oil film bearings［J] .Journal of Vibration and Acoustics， Transactions of the ASME， 1990， 112:160－165.

［7] Singh M P，Chang T S，Suarez L E.A responses spectrum method for seismic design evaluation of rotating machines［J] .Journal of Vibration and Acoustics，Transactions of the ASME，1992，114:454 －460.

［8] Gaganis B J，Zisimopoulos A K，Nikolakopoulos PG，et al.Modal analysis of rotor on piecewise linear journal bearing under seismic excitation［J] .Journal of Vibration and Acoustics， Transactions of the ASME， 1999， 121:190－196.

［9] Kostarev V V，Petrenko A V，Vasilyev P S.An advanced seismic analysis of an NPP powerful turbogenerator on an isolation pedestal［J] .Nuclear Engineering and Design，2007，237(12 －13):1315 －1324.

［10] Kessler C L.Complex modal analysis of rotating machinery［D] .UMI Microform 9933653，1999.

［11] Lee CW.Vibration analysis of rotors［M] .Kluwer Academic Publishers，1993.

［12] 汪 玉，华宏星.舰船现代冲击理论及应用［M] .北京:科学出版社，2005.

［13] 刘建湖.舰船非接触水下爆炸动力学的理论与应用［D] .无锡:中国舰船科学研究中心，2002.