航空动力学报
JOURNAL OF AEROSPACE   POWER
1998年 第13卷  第2期科技期刊
转子系统瞬态响应计算的
精细迭代积分方法
西安交通大学　　姜培林　虞　烈
【摘要】　提出了计算转子系统瞬态响应的精细迭代积分方法。采用精细积分法计算状态转移矩阵，可以得到齐次方程的精确解；结合轴承所提供的油膜力的特点，提出一种隐式积分法，允许使用较大的步长，并对各矩阵及向量适当分块，使非线形方程组的维数大为减小。此方法无条件稳定，可以保证较高的精度与效率。
主题词：　转子　瞬态响应　积分法　迭代法
分类号：　V231.9　O317
1　引　言
       转子—轴承系统的非线性动力学研究近来引起了广泛的兴趣，系统瞬态响应的计算是其中一个重要环节。至今几乎所有的求解结构动力方程的积分方法都被应用于转子系统瞬态响应的计算中，如Wilson-θ法、Newmark-β法、Euler法、Houbolt法以及Runge-Kutta法等，但这些方法普遍存在计算稳定性与精度的问题。钟万勰提出了计算结构瞬态动力响应的精细积分法［1］，文献［2］曾指出这种方法无条件稳定，对于刚性方程采用精细积分法可以得到较为准确的解，因此对于齐次方程来说，其步长不受限制。精细积分法的优越性已为人们所逐渐认识。沈为平、林家浩研究了各种外载作用时的精细积分法［3，4］。已有的有关精细积分法的应用都是求解初值齐次方程或求与时间相关的载荷下的响应。由于轴承所提供的油膜力与轴颈的运动状态有关，其中并不显含时间，因此现有的这些方法不适合于转子—轴承系统瞬态响应的计算。本文充分利用精细积分法计算精度高、稳定性好的优点，提出了一种分块迭代的方法。
2　方法描述
　　由滑动轴承或挤压油膜阻尼器支承的转子系统，其运动方程可以表示为如下的形式：

(1)
初值x(t0)=x0，(t0)=0，其中M，G，K∈RN×N为质量、阻尼及刚度矩阵，一般为时不变矩阵；x∈RN为系统的广义坐标；r=rt+rn,rt∈RN为仅与时间有关的力矢量，如不平衡力等，rn∈RN(x，)为轴承提供的油膜力矢量。
　　根据文献［1］中提供的变换方法，将式(1)变形为：
                           

其中：T(τ)=exp(τH)　，　　τ=t-t0　　
   精细积分法的特点在于状态转移矩阵T(τ)的计算上，文献［1］中提出了 一种2w算法来计算T(τ)，它将时间分成2w等份，即△t=τ/2w，对exp(△tH)进行泰勒级数展开，而exp(Hτ)的计算则利用指数矩阵的特点对exp(△tH)进行2w次自乘得到。它引起的截尾与一般的计算机的截断误差处于同一数量级上，因此得到的T(τ)可以认为是精确解。
　　对于转子系统来说，极有可能存在刚性位移(如对称转子)，因此H-1不存在，此时不可以由文献［1］中给出的解析式来计算νk。本文作如下的处理。
　　由于油膜力与轴颈的位移与速度有关，而二者都是待求量，因此必须引入必要的假设。本文设f在区间［tk-1，tk］上呈线性变化，即f(t)=fnk-1+(fnk-fnk-1)(t-tk-1)/τ，由式(3)有：
νk=T(tk-tk-1)νk-1+∫tktk-1T(tk-s)［fnk-1+(fnk-fnk-1)(t-tk-1)/τ］ds
+∫tktk-1T(tk-s)ft(s)ds
=Tνk-1+Ffnk-1+Sfnk+btk=b′k
(4)
其中：b′k-1=Tνk-1+Ffnk-1+btk，F=∫tktk-1［1-(s-tk-1)/τ］T(tk-s)ds
      S=∫tktk-1(tk-s)/τT(tk-s)ds，btk＝∫tktk-1T(tk-s)ft(s）ds　　
   上式中各积分采用高斯型数值积分计算，因此可以避免因H阵奇异带来的困难。T(tj)的计算仍采用2w算法。一般tk-tk-1为常数，则F，S及T为常数，只需在开始时计算一次。由于fnk与νk有关，式(4)为隐式，因此必须迭代求解νk与fnk，对于大型转子—轴承系统来说，由于转子分段数多，自由度较大，因此表示的非线性方程组的维数是很高的。
　　实际上，轴承的油膜力仅作用于很少的几个位置上，仅与该处的运动状态有关，式(4)中的线性自由度远多于非线性自由度，因此将式(4)中的各矩阵及向量写成分块形式：

其中：表示仅与轴承力有关的坐标，x′1∈RN1表示其它坐标。

因此，有
(5)
(A′22-A′21A′-111A′12)x′2k＝b′2k-1+S22fn2k-A′21A′-111(b′k-1+S12fn2k)
(6)
　这是一低维的非线性方程组，采用迭代方法求解时，每一次迭代中，只需重新计算fn2k与x′2k，因此可以大幅度地减小计算量，尤其对于转子单元较多时。由式(10)可以得到x′1k
x′1k=A′-111(b′1k-1+S12fn2k-A′12x′12)
(7)
3 算　例　
　采用图1所示单圆盘对称转子的瞬态响应。用有限元法对转子进行离散建立系统的运动


图1　转子—轴承系统示意图
方程，方程中包含转子转角自由度，共有6个离散节点。轴承提供的油膜力采用文献［5］中提供的动态油膜力数据库方法获得。计算时采用如下的参数：转子D=4600 mm，d=400 mm，l=200 mm，t=200 mm，L=2600 mm；轴承长径比l/d=0.5，径向间隙C=0.4 mm，轴承有均匀分布的3块固定瓦。瓦张角φ0=100°，铅垂线通过第一瓦的中点；物性参数为：转子的密度ρ=7750 kg/m3，弹性模量E=2.1×1011Pa，润滑油的动力粘度μ=0.015 Pa·s；轴的转速ω=3000 r/min。


图2　τ=T/100时转子系统的瞬态响应(a)左轴承的轴心轨迹(b)左轴承轴心垂直位移的时间历程(c)左轴承轴心水平位移的时间历程(d)左轴承处轴颈 水平平面内的转角的时间历程
    转子由轴承的中心处自由释放，图2为系统的响应。由图可知，转子的轴心轨迹及轴的转角很快收敛到其静态平衡位置上。从图3可以看出：即使采用较大的步长，计算结果也比较可靠。同时与不分块的计算相比，分块以后计算速度提高了近10倍。


图3　τ=T/50时左轴承的轴心轨迹　图4　突加不平衡量时左轴承的轴心轨迹
    工程中常发生一种由叶片断裂等引起的突加不平衡量的故障。本文采用精细积分法对突加不平衡量时机组的响应进行了模拟。开始时 圆盘处有偏心距为0.04 mm的不平衡量，突加以偏心距为0.4 mm的不平衡量，左轴承的轴心轨迹如图4所示。此时机组的转速为1000 r/min，积分步长为T/500。既使采用较大的积分步长，也可以得到较为稳定的解。
4　结　论
(1)本文采用数值积分法计算状态转移矩阵，可以避免因刚体位移的存在而无法计算H-1的困难。(2)由于采用隐式，可以在较大的积分步长下保持较好的数值稳定性。(3)采取分块的方法将非线性代数方程组分成两部分：一部分为需迭代求解的非线性方程组，它的维数仅与轴承运动自由度数有关；另一部分为线性代数方程组，其维数虽然有时很高，但不需要迭代求解。对于大型转子系统来说，支承所占的自由度仅占很少的一部分，所以这种分块方法的优越性更显著。可以说本方法在计算稳定性、计算精度以及效率各方面都比其它方法优越。这种方法可以应有于外载荷与运动状态有关而不显含时间的非线性动力方程求解。
男　28　博士生　西安交通大学润滑理论及轴承研究所　
                                          参　考　文　献
1　Zhong W X,Williams F W.A Precise Time Step Integration Method.Proc.IMechE,1994,208(C6):427-430
2　Lin J H,Shen W P,Williams F W.A High Precision Direct Integration Scheme for Structures Subject to Transient Dynamic Loading.Computers and Structures,1995,56(1):113-120
3　Shen W P,Lin J H,Williams F W.Parallel Computing for the High Precision Direct Ingegration Method.Comput.Methods Appl.Mech Engg.,1995,126:315-331
4　沈为平，宋华茂．任意荷载下结构动力响应的并行算法．振动工程学报，1996，9(4)：333-340
5　Wang Wen ,Zhang Zhiming.Calculation of Journal Dynamic Locus Aided by Database of Non-Stationary Oil Film Force of Single Bush Segment.Asia-Pactific Vibration Conference 93,Japan:365-369
（责任编辑　王震华）
