一種基于頻響因子的轉(zhuǎn)子系統(tǒng)基礎(chǔ)沖擊響應(yīng)計(jì)算方法

謝最偉，賀少華，吳新躍

(海軍工程大學(xué) 船舶與動(dòng)力學(xué)院機(jī)械工程系，武漢 430033)

對于動(dòng)力學(xué)系統(tǒng)響應(yīng)的求解即動(dòng)力學(xué)微分方程的求解，一般有三種解法，解析解法、近似解析解法和數(shù)值解法，只有極少數(shù)動(dòng)力學(xué)方程存在解析解，所以大部分動(dòng)力學(xué)方程只能采用近似解析解法和數(shù)值解法，近似解析方法主要針對弱非線性系統(tǒng)的穩(wěn)態(tài)響應(yīng)的求解，如攝動(dòng)法、諧波平衡法、多尺度方法、平均法等，而對于數(shù)值解法，則適應(yīng)于任何的系統(tǒng)響應(yīng)求解，主要包括顯式算法中的中心差分法、Runge-Kutta法和隱式算法中的 Wilson-θ法、Houbolt法、Newmark 法等［1－4]。

動(dòng)力學(xué)系統(tǒng)(結(jié)構(gòu))外界沖擊本文的定義:沖擊是一種振動(dòng)系統(tǒng)運(yùn)動(dòng)量的瞬時(shí)突然變化，其動(dòng)能傳輸是在一個(gè)與振動(dòng)系統(tǒng)固有周期相比的短時(shí)間內(nèi)進(jìn)行的，沖擊因其時(shí)間歷程比較復(fù)雜而不同于撞擊。以艦載設(shè)備的水下非接觸爆炸沖擊為例，沖擊理解為對艦艇結(jié)合部(外殼、甲板、艙壁等)的一種非接觸的相當(dāng)遠(yuǎn)距離的爆炸(當(dāng)水下爆炸時(shí)，氣泡不觸及外殼結(jié)構(gòu)或受壓物體結(jié)構(gòu))作用，其結(jié)果為在艦船正常運(yùn)行中不會(huì)出現(xiàn)的一種振動(dòng)系統(tǒng)“艦艇”運(yùn)動(dòng)量的巨大干擾。距離很近的水下爆炸或空中爆炸所引起的作用，另作特別處理，所以，結(jié)構(gòu)沖擊理解為:在承壓的艦艇結(jié)合部中和在設(shè)備基座中的機(jī)械干擾。沖擊載荷理解為由于爆炸產(chǎn)生的機(jī)械載荷。

旋轉(zhuǎn)機(jī)械的瞬態(tài)沖擊響應(yīng)分析區(qū)別于一般的靜態(tài)條件下結(jié)構(gòu)的相應(yīng)分析，主要體現(xiàn)在以下兩個(gè)方面:①軸承－轉(zhuǎn)子的交互效應(yīng);② 陀螺效應(yīng)。除了以上兩個(gè)因素，當(dāng)系統(tǒng)基礎(chǔ)遭受旋轉(zhuǎn)運(yùn)動(dòng)激勵(lì)時(shí)Coriolis效應(yīng)將不可忽略。在過去的幾十年中，國內(nèi)外學(xué)者對于旋轉(zhuǎn)機(jī)械(轉(zhuǎn)子－軸承系統(tǒng))在地震激勵(lì)下的響應(yīng)進(jìn)行了研究［5－9]，取得了一些成果，如:旋轉(zhuǎn)激勵(lì)成分對響應(yīng)的影響;線性理論的應(yīng)用局限性;數(shù)值算法的改進(jìn)等等。

考查地震或水下非接觸爆炸沖擊下旋轉(zhuǎn)機(jī)械的動(dòng)態(tài)響應(yīng)特性，一般從研究轉(zhuǎn)子系統(tǒng)基礎(chǔ)沖擊響應(yīng)出發(fā)。對于轉(zhuǎn)子－軸承系統(tǒng)的沖擊響應(yīng)研究，主要方法為:首先對轉(zhuǎn)子－軸承系統(tǒng)進(jìn)行建模，用基礎(chǔ)運(yùn)動(dòng)(位移、速度和加速度)來模擬基礎(chǔ)沖擊，建立系統(tǒng)力學(xué)方程;然后通過數(shù)值計(jì)算求解方程得到響應(yīng)。通過響應(yīng)分析檢驗(yàn)以下幾個(gè)方面的內(nèi)容:①潤滑油膜在任何時(shí)候必須保持一個(gè)最小的油膜厚度以避免轉(zhuǎn)子和軸承不會(huì)發(fā)生碰摩;②軸承基礎(chǔ)能夠承受沖擊軸承反力;③轉(zhuǎn)子沖擊位移(變形)或應(yīng)力在安全的范圍之內(nèi)。由于陀螺效應(yīng)和軸承的交互效應(yīng)，轉(zhuǎn)子系統(tǒng)運(yùn)動(dòng)方程系數(shù)矩陣呈非對稱性，不能在模態(tài)坐標(biāo)下解耦，無法利用常規(guī)模態(tài)疊加法求解，所以已往的研究一般采用數(shù)值積分如Newmark法等進(jìn)行迭代求解，隨著有限元和計(jì)算機(jī)技術(shù)水平的飛速發(fā)展，目前，借助嵌入了這些數(shù)值積分算法的商業(yè)有限元軟件進(jìn)行沖擊響應(yīng)計(jì)算分析成為了國內(nèi)的主流。但是，必須指出的是，絕大部分有限元商業(yè)軟件的知識產(chǎn)權(quán)屬于國外，過分依賴這些商業(yè)軟件將使我國本已滯后的機(jī)械設(shè)備抗沖研究特別是理論研究更加落后國外。同時(shí)，數(shù)值積分法相對線性疊加法要耗費(fèi)更多的計(jì)算資源和計(jì)算時(shí)間。

轉(zhuǎn)子系統(tǒng)復(fù)模態(tài)分析法相對實(shí)模態(tài)分析法能夠更好地表達(dá)轉(zhuǎn)子的運(yùn)動(dòng)規(guī)律，因?yàn)閺?fù)模態(tài)包含了轉(zhuǎn)子的旋轉(zhuǎn)方向和振動(dòng)幅值雙重信息［10]，本文正是從這一思想出發(fā)，將提出并實(shí)現(xiàn)一種復(fù)數(shù)域內(nèi)轉(zhuǎn)子系統(tǒng)沖擊響應(yīng)計(jì)算方法，無需坐標(biāo)解耦但仍可以利用線性疊加法進(jìn)行響應(yīng)求解，從而克服常規(guī)模態(tài)疊加法和傳統(tǒng)數(shù)值積分算法的上述缺點(diǎn)，為轉(zhuǎn)子系統(tǒng)提供一種普適的沖擊瞬態(tài)響應(yīng)線性疊加計(jì)算方法。

1 理論推導(dǎo)

1.1 復(fù)數(shù)域內(nèi)轉(zhuǎn)子系統(tǒng)運(yùn)動(dòng)微分方程的推導(dǎo)

一個(gè)基礎(chǔ)受加速度運(yùn)動(dòng)激勵(lì)的n維線性轉(zhuǎn)子系統(tǒng)(軸線沿x方向)的運(yùn)動(dòng)方程可以寫為:

式中:Y、Z描述著基礎(chǔ)的運(yùn)動(dòng)，y、z為系統(tǒng)相對基礎(chǔ)的位移(變形)向量，它們均為n維。

令:

則:

其中:

將式(2)代入式(1)，等式兩邊左乘［T]－1:

進(jìn)一步簡寫成:

1.2 特征方程的演變和響應(yīng)求解

對于g(t)，有下面的傅里葉展開式

對于線性系統(tǒng)，同樣有解的形式為:

式中的Gf、Gb、Pf和Pb均為n維向量。由上式我們可以發(fā)現(xiàn)，轉(zhuǎn)子在YZ平面的任一頻率運(yùn)動(dòng)被分解成了相同頻率的一正一反兩個(gè)圓周運(yùn)動(dòng)，也就是通常所說的正向和反向渦動(dòng)(forward和backward rotating)。當(dāng)兩渦動(dòng)運(yùn)動(dòng)的幅值相等即時(shí)，合成運(yùn)動(dòng)為直線運(yùn)動(dòng)，這對應(yīng)于一般的非旋轉(zhuǎn)振動(dòng)系統(tǒng)的情況;當(dāng)兩渦動(dòng)運(yùn)動(dòng)的幅值不相等時(shí)，合成運(yùn)動(dòng)為橢圓運(yùn)動(dòng)。所以，可以將一般的非旋轉(zhuǎn)振動(dòng)系統(tǒng)的直線運(yùn)動(dòng)看成渦動(dòng)的一種特例。

將式(4)，式(5)代入式(3)，方程兩邊相等的唯一條件是方程兩邊相同頻率ejωkt項(xiàng)前系數(shù)相等，相同頻率e－jωkt前系數(shù)也相等。對于任一給定k，可以得到:

上面兩式表達(dá)的信息是相同的，所以只需考查其中的任一個(gè)。

但是，對于轉(zhuǎn)子系統(tǒng)，由于陀螺效應(yīng)和軸承－轉(zhuǎn)子的交互效應(yīng)，使得方程(3)的系數(shù)矩陣一般不再對稱，模態(tài)向量關(guān)于系數(shù)矩陣的正交性不再存在，無法將系統(tǒng)方程在模態(tài)坐標(biāo)下解耦，模態(tài)疊加法求解系統(tǒng)響應(yīng)失敗。

將式(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的形式為:

即為系統(tǒng)的模態(tài)向量。故式(7)可進(jìn)一步演變?yōu)?

當(dāng)為各向同性轉(zhuǎn)子時(shí)，A12=A21=0。在本文中，將A11、A12、A21和 A22稱為“頻響因子(frequency response coefficient)”。

如上所述，得到各頻率成分響應(yīng)后，將響應(yīng)合成即可得到系統(tǒng)總的響應(yīng)。在實(shí)際中，基于這樣的常識:當(dāng)結(jié)構(gòu)受沖擊作用時(shí)，在某種意義上結(jié)構(gòu)的響應(yīng)本質(zhì)上是一種振蕩，響應(yīng)時(shí)主要出現(xiàn)的頻率是沖擊的優(yōu)勢頻率和結(jié)構(gòu)的固有頻率。所以，只需考查結(jié)構(gòu)固有頻率和沖擊優(yōu)勢頻率處的響應(yīng)值。

1.3 響應(yīng)合成

在實(shí)際機(jī)械抗沖考核中，我們需要計(jì)算其在標(biāo)準(zhǔn)沖擊激勵(lì)下各部分的最大瞬態(tài)響應(yīng)，并不需要整個(gè)時(shí)間歷程的響應(yīng)。通過上小節(jié)求得各頻率處的沖擊響應(yīng)，需要對各頻響幅值進(jìn)行組合。組合方法可以仿照美國海軍 DDAM(dynamic design analysis method)［12]，使用下列三種方法之一:

(1)絕對值求和(ABS):

本文公式:

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

本文公式:

(3)美國海軍研究實(shí)驗(yàn)室求和(NRL)

本文公式:

式中:Pf(ωj)和Pb(ωj)為任一頻率正反渦動(dòng)幅值。

2 實(shí)例分析

2.1 某剛性轉(zhuǎn)子系統(tǒng)模型和參數(shù)

一柔性支撐剛性轉(zhuǎn)子系統(tǒng)如圖1所示，系統(tǒng)的詳細(xì)參數(shù)如表1所示，該系統(tǒng)的運(yùn)動(dòng)狀態(tài)完全由旋轉(zhuǎn)軸兩端與支撐連接點(diǎn)的運(yùn)動(dòng)狀態(tài)決定。

圖1 剛性轉(zhuǎn)子系統(tǒng)的坐標(biāo)系和支撐Fig.1 Coordinate system and supports of a rigid rotor system

表1 系統(tǒng)參數(shù)Tab.1 System parameters

2.2 運(yùn)動(dòng)微分方程的建立

建立如式(1)的系統(tǒng)相對基礎(chǔ)的運(yùn)動(dòng)方程式:

式中:

{FY1FY2FZ1為支撐軸承等效的彈簧－阻尼單元支撐力。在上式中，與式(1)不同的地方在于它還考慮了基礎(chǔ)的速度激勵(lì)，由于地震沖擊或是艦載設(shè)備水下非接觸爆炸沖擊一般只采用加速度形式來表示，所以在這里忽略速度激勵(lì)項(xiàng)。

進(jìn)一步的，令:

則:

式中:

正如上文所述，方程中速度向量前的矩陣加進(jìn)了陀螺效應(yīng)和軸承等效彈簧－阻尼單元的貢獻(xiàn)，G為反對稱矩陣(skew-symmetric)(如果a=b)，而如果軸承的交叉阻尼項(xiàng)又不相等，這使得總體阻尼效應(yīng)矩陣變得十分復(fù)雜，不再正對稱，這與一般非旋轉(zhuǎn)振動(dòng)系統(tǒng)的阻尼效應(yīng)矩陣總被處理成對稱的(如比例粘性阻尼)是不同的，剛度效應(yīng)矩陣由于軸承等效彈簧－阻尼單元的貢獻(xiàn)一般也不再對稱(如果軸承交叉剛度項(xiàng)不相等)。

將方程(10)轉(zhuǎn)換成式(3)的形式:

式中:

后續(xù)公式的推導(dǎo)過程與上一小節(jié)相同，不再贅述。

2.3 計(jì)算

最終計(jì)算得到的系統(tǒng)特征值λ和模態(tài)向量{Pf1，i=1，2，3，4，如表 2 所示。

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

同時(shí)，對應(yīng)于式(8)的算術(shù)表達(dá)式為:

按照本文方法計(jì)算在圖3所示加速度沖擊激勵(lì)下系統(tǒng)的響應(yīng)，首先，按照上文所示方法將激勵(lì)信號經(jīng)傅里葉變換成圖3右所示頻率－幅值曲線，很顯然，在本例中，圖3的幅值(用A表示)與式(4)中Gf、Gb的關(guān)系為A=G=G。根據(jù)式(12)、圖2和圖3，在系統(tǒng)固

fb有頻率 λ1= λ2=123.2 rad/s，λ3=416.8 rad/s，λ4=456.9 rad/s和沖擊優(yōu)勢頻率 λ5=9.8 rad/s處，選取響應(yīng)合成方法為上文中的絕對值求和(ABS)，最終計(jì)算得到的最大位移響應(yīng)為4.30 m。

然后，采用經(jīng)典Newmark數(shù)值積分方法在時(shí)域直接求解系統(tǒng)響應(yīng)，由于Newmark數(shù)值積分方法無條件穩(wěn)定且具有二階精度，在沖擊和瞬態(tài)響應(yīng)計(jì)算中最為常用［13]。Newmark數(shù)值積分法計(jì)算得到的最大位移值為3.56 m。

2.4 誤差分析及說明

上小節(jié)計(jì)算結(jié)果表明，直接數(shù)值積分方法計(jì)算結(jié)果要明顯小于本文方法計(jì)算結(jié)果，誤差的主要原因在于本文方法類似于譜分析方法如DDAM，都是一種典型的保守計(jì)算方法，在響應(yīng)合成時(shí)沒有考慮各模態(tài)響應(yīng)的相位差，進(jìn)行的是簡單的“絕對值”相加。文獻(xiàn)［13] 指出，即使為非密集模態(tài)情況，兩種合成方法(NRL和SRSS)的結(jié)果與時(shí)域合成結(jié)果相比，在特定條件下也能產(chǎn)生40%的相對誤差，但這并不會(huì)妨礙其在工程中的應(yīng)用，仍具有工程應(yīng)用性，如DDAM，因?yàn)樵趯?shí)際中，設(shè)備的抗沖擊設(shè)計(jì)本身就需要一個(gè)安全裕值。本文提出的方法的誤差在理想的范圍之內(nèi)，所以可用于工程實(shí)際。

3 結(jié)論

由于陀螺效應(yīng)和轉(zhuǎn)子－軸承的交互效應(yīng)，旋轉(zhuǎn)機(jī)械的沖擊瞬態(tài)響應(yīng)計(jì)算比一般系統(tǒng)要復(fù)雜得多，以往采用的數(shù)值積分算法雖能解決任意復(fù)雜系統(tǒng)的響應(yīng)求解問題，但這些算法往往需要耗費(fèi)大量的計(jì)算資源和計(jì)算時(shí)間，為突破這一瓶頸，本文提出了一種復(fù)數(shù)域內(nèi)的轉(zhuǎn)子系統(tǒng)基礎(chǔ)沖擊響應(yīng)計(jì)算方法，由于采用了線性疊加原理，該方法較好地克服了傳統(tǒng)數(shù)值積分方法的上述缺點(diǎn)。同時(shí)，由于轉(zhuǎn)子系統(tǒng)運(yùn)動(dòng)方程系數(shù)矩陣呈非對稱性，不能在模態(tài)坐標(biāo)下解耦，無法利用常規(guī)模態(tài)疊加法求解，所以，本文提出的方法可以視為一種改進(jìn)的瞬態(tài)響應(yīng)線性疊加計(jì)算方法。首先將激勵(lì)和響應(yīng)傅里葉展開成復(fù)數(shù)形式，包括正向旋轉(zhuǎn)項(xiàng)和反向旋轉(zhuǎn)項(xiàng)，根據(jù)方程左右兩邊相同頻率前系數(shù)相等的事實(shí)得到特征方程，將特征方程寫成簡單矩陣束的本征方程形式，求得矩陣束的本征值和本征向量，將本征向量正規(guī)化，進(jìn)一步得到矩陣束的逆陣，逆陣元素即為“頻響因子”，將逆陣與激勵(lì)相乘即可得到頻率響應(yīng)幅值，將所有頻率響應(yīng)成分疊加即可得到系統(tǒng)響應(yīng)。通過一個(gè)工程實(shí)例，比較了所提方法與Newmark數(shù)值積分方法的計(jì)算結(jié)果，誤差在合理的范圍內(nèi)，從而驗(yàn)證了其在工程應(yīng)用中的合法性。本文所提方法可以作為轉(zhuǎn)子系統(tǒng)基礎(chǔ)沖擊響應(yīng)和瞬態(tài)響應(yīng)計(jì)算的一種普適方法。“頻響因子”在本文中為首次提出。

［1] 汪 玉，趙建華，杜儉業(yè)，等.基于多體動(dòng)力學(xué)有限元計(jì)算的一種艦用柴油機(jī)抗沖擊性能仿真［J] .振動(dòng)與沖擊，2009，28(11):87 －90，129.

［2] 高浩鵬，黃映云，王雪山.船用柴油機(jī)運(yùn)動(dòng)機(jī)構(gòu)動(dòng)力學(xué)動(dòng)態(tài)分析方法研究［J] .振動(dòng)與沖擊，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] 汪 玉，華宏星.艦船現(xiàn)代沖擊理論及應(yīng)用［M] .北京:科學(xué)出版社，2005.

［13] 劉建湖.艦船非接觸水下爆炸動(dòng)力學(xué)的理論與應(yīng)用［D] .無錫:中國艦船科學(xué)研究中心，2002.