用于推進(jìn)軸系振動分析的改進(jìn)數(shù)值組裝法

巫 頔， 謝溪凌， 張志誼, 2

(1. 上海交通大學(xué) 機(jī)械系統(tǒng)與振動國家重點(diǎn)實(shí)驗(yàn)室， 上海 200240； 2. 上海交通大學(xué) 振動、沖擊、噪聲研究所， 上海 200240)

船舶推進(jìn)軸系振動校核是一項(xiàng)重要的設(shè)計(jì)內(nèi)容，高精度的動力學(xué)模型和計(jì)算方法是獲得較精確的振動特性和振動響應(yīng)的基礎(chǔ)。有限元法是目前廣泛使用的數(shù)值方法，但對于諸如梁的一維結(jié)構(gòu)，可以通過振動微分方程在梁的自然分段內(nèi)精確描述動力響應(yīng)，得到精度更高的計(jì)算結(jié)果。目前，已有多種變截面梁橫向振動精確求解方法，如傳遞矩陣法、譜元法、數(shù)值組裝法等，其中傳遞矩陣法在船舶軸系振動分析中應(yīng)用最為廣泛。傳遞矩陣法中一般取梁單元端點(diǎn)撓度、轉(zhuǎn)角、彎矩與剪力作為梁單元的狀態(tài)變量，通過梁單元內(nèi)部的傳遞關(guān)系和節(jié)點(diǎn)連續(xù)性條件，將各單元傳遞矩陣按序累乘[1]。傳遞矩陣法縮減了系統(tǒng)方程階數(shù)，計(jì)算速度快，但是系統(tǒng)矩陣的病態(tài)易導(dǎo)致特征方程的數(shù)值不穩(wěn)定，且梁單元內(nèi)部振動信息無法獲得。Wu等[2]從界面連續(xù)性條件出發(fā)，提出連續(xù)質(zhì)量傳遞矩陣法(continuous-mass transfer matrix method, CTMM)，并用于求解多跨變截面Timoshenko梁橫向振動特征值問題。CTMM方法采用梁單元振動微分方程的4個系數(shù)作為梁單元的狀態(tài)變量，便于連續(xù)表示梁單元內(nèi)部振動。Zhang等[3]在CTMM中考慮了梁的橫向與扭轉(zhuǎn)彈簧約束、附加彈簧-質(zhì)量系統(tǒng)、彈簧-質(zhì)量-彈簧約束以及聯(lián)軸器連接等各種常見邊界條件。徐穎蕾[4]推導(dǎo)了CTMM的模態(tài)正交性，并使用模態(tài)疊加法求解軸系的強(qiáng)迫振動響應(yīng)。Xie等[5]在界面連續(xù)性條件中加入了外載荷向量，并在頻域直接求解變截面Timoshenko梁的強(qiáng)迫振動響應(yīng)，避免了模態(tài)疊加過程中的截?cái)嗾`差。CTMM方法融合了傳統(tǒng)傳遞矩陣法的優(yōu)勢，但依然存在高頻數(shù)值發(fā)散問題，影響特征頻率與強(qiáng)迫振動響應(yīng)的求解。

為避免固有頻率求解中的數(shù)值發(fā)散問題，Riccati傳遞矩陣法[6]將經(jīng)典傳遞矩陣法中兩端狀態(tài)變量的對應(yīng)關(guān)系轉(zhuǎn)換為在每一截面上位移與載荷的對應(yīng)關(guān)系，減小了矩陣?yán)鄢苏`差，提高了系統(tǒng)行列式的數(shù)值穩(wěn)定性，但奇點(diǎn)干擾和多次求逆帶來的數(shù)值誤差仍會影響固有頻率的求解精度。王正[7]對Riccati傳遞矩陣法在求特征值時(shí)遇到的奇點(diǎn)問題進(jìn)行了討論，并提出了消除奇點(diǎn)干擾后的頻率方程式；黃娟等[8]在計(jì)算汽輪發(fā)電機(jī)軸系時(shí)使用消去奇點(diǎn)干擾的改進(jìn)Riccati方法。Bestle等[9]為提高Euler梁橫向振動的計(jì)算頻率范圍，提出了拆分梁段與變換局部坐標(biāo)系方向等方法，但對高頻數(shù)值發(fā)散的抑制效果有限，無法顯著提高求解頻率范圍。

Wu等[10]提出的數(shù)值組裝法(numerical assembly method，NAM)以梁單元微分方程的4個系數(shù)作為狀態(tài)變量，將n個梁單元對應(yīng)的4n個方程一次性求解。由于矩陣方程只求解一次，避免了傳遞矩陣法多次求解矩陣方程帶來的誤差放大問題，因此計(jì)算頻率范圍相較之有一定提高。Wu等解決了梁上分布附加彈簧-質(zhì)量塊的均勻Euler-Bernoulli梁特征值問題。Lin等[11]利用NAM計(jì)算了附加集中質(zhì)量、集中轉(zhuǎn)動慣量、彈簧與扭簧等多種集中元素的多跨Timoshenko梁的固有頻率和模態(tài)振型。Yesilce等[12]使用NAM分析附加多個彈簧與質(zhì)量塊的多跨變截面梁，考慮了軸向作用力對系統(tǒng)橫向振動固有頻率的影響。Xu等[13]使用NAM計(jì)算Euler-Bernoulli階梯梁，通過局部坐標(biāo)系拓寬計(jì)算頻帶。NAM避免了誤差的累乘放大，因而能夠提高計(jì)算頻率范圍。但對于梁段單元參數(shù)差異過大的系統(tǒng)，通過NAM法建立的系統(tǒng)矩陣仍存在病態(tài)嚴(yán)重的問題，制約了軸系高頻段振動的求解。

本文首先分析了現(xiàn)有計(jì)算方法的誤差來源，然后使用局部坐標(biāo)系推導(dǎo)多跨Timoshenko階梯梁的橫向振動方程，并考慮推進(jìn)軸系上的質(zhì)量、慣量、彈簧及聯(lián)軸器等集中元素，進(jìn)而給出行歸一化措施，降低矩陣方程的條件數(shù)，提高推進(jìn)軸系高頻振動計(jì)算精度。

1 基于Timoshenko梁理論的改進(jìn)數(shù)值組裝法

1.1 Timoshenko梁橫向振動方程

考慮轉(zhuǎn)動慣量和剪切變形對橫向振動的影響，設(shè)第i梁段橫向振動位移為yi(x,t)，轉(zhuǎn)角為φi(x,t)，對應(yīng)的運(yùn)動方程可表示為

(1)

(2)

式中：ρi為材料密度；Ii為截面轉(zhuǎn)動慣量；Ei為彈性模量；Gi為剪切模量；Ai為截面面積；κi為Timoshenko梁剪切系數(shù)。設(shè)yi(x,t)=Yi(x)ejωt，φi(x,t)=Φi(x)ejωt，將轉(zhuǎn)角φi與撓度yi解耦可得

(3)

(4)

式(4)為4階齊次微分方程，撓度Yi(x)的解為振型函數(shù)，如式(5)所示

Yi(x)=C1icoshλi(x-Xi-1)+C2isinhλi(x-Xi-1)+

(5)

類似的，轉(zhuǎn)角的振型函數(shù)Φi(x)為

Φi(x)=C1iqisinhλi(x-Xi-1)+C2iqicoshλi(x-

(6)

其中

1.2 界面連續(xù)性條件

在軸系橫向振動時(shí)，梁段節(jié)點(diǎn)滿足撓度、轉(zhuǎn)角、彎矩與剪力的協(xié)調(diào)條件，如式(7)

(7)

式中：KRi為梁段交界處轉(zhuǎn)角彈簧剛度；KTi為線彈簧剛度。當(dāng)?shù)趇節(jié)點(diǎn)有集中質(zhì)量mi和集中轉(zhuǎn)動慣量Ji時(shí)，KRi=-ω2Ji，KTi=-ω2mi。在局部坐標(biāo)系下，第i梁段的坐標(biāo)原點(diǎn)定義為Xi-1,R，Xi,L為第i梁段的終點(diǎn)。將Xi,R=0，Xi,L=li代入式(6)，并將撓度、轉(zhuǎn)角、彎矩與剪力的表達(dá)式展開，可得矩陣方程

(8)

其中

方程式(8)可簡寫為AcriC(i)=Acli+1C(i+1)，Acri即為第i梁段右端參數(shù)矩陣，Acli+1為第i+1梁段左端參數(shù)矩陣，C(i)為第i梁段微分方程解的系數(shù)向量。若第i節(jié)點(diǎn)有聯(lián)軸器，則節(jié)點(diǎn)位移與轉(zhuǎn)角不連續(xù)，式(7)中的轉(zhuǎn)角與位移協(xié)調(diào)條件變?yōu)?/p>

Yi(Xi,L)-Yi+1(Xi,R)=

(9)

(10)

式中：KITi為聯(lián)軸器橫向剛度；KIRi為轉(zhuǎn)角剛度。相應(yīng)的，矩陣方程式(8)的形式變?yōu)?/p>

(11)

1.3 邊界條件

對于n段變截面梁，共有n+1個節(jié)點(diǎn)，其中有n-1個內(nèi)節(jié)點(diǎn)，每個節(jié)點(diǎn)存在4個連續(xù)性條件構(gòu)成的界面連續(xù)性條件方程組；兩個邊界節(jié)點(diǎn)，各有兩個邊界條件方程。不同邊界條件對應(yīng)的方程組形式不同，在無位移約束，兩端有集中質(zhì)量、集中轉(zhuǎn)動慣量或彈簧支承的條件下，梁段邊界的剪力與彎矩平衡條件表示為

左端

(12)

右端

(13)

式中：KTL，KTR分別為左右兩端的橫向彈簧；KRL，KRR為左右兩端的扭轉(zhuǎn)彈簧；ML，MR，JL，JR為各端的集中質(zhì)量與集中轉(zhuǎn)動慣量。利用式(12)、式(13)的邊界條件方程組可以得到邊界節(jié)點(diǎn)的內(nèi)力表達(dá)式

BLC1=0,BRCn=0

(14)

式中：BL與BR為2×4矩陣，為自由-自由條件下邊界節(jié)點(diǎn)處的剪力和彎矩為零。矩陣BL與BR見附錄A所示。在無外力作用的條件下，將各梁段方程依次組裝并求解，即為NAM的矩陣方程

(15)

記方程式(15)為T(ω)C(ω)=0，通過求解T(ω)的特征方程T(ω)=0可求得系統(tǒng)固有振動特性；求解強(qiáng)迫振動響應(yīng)問題時(shí)，外力向量在方程右側(cè)，矩陣方程可簡寫為T(ω)C(ω)=F。

1.4 行歸一化加權(quán)陣

隨著頻率升高，NAM的特征方程(T(ω))曲線會在高頻出現(xiàn)數(shù)值發(fā)散現(xiàn)象。觀察方程式(14)可以發(fā)現(xiàn)，矩陣T(ω)的每一行與系數(shù)C的乘積都有明顯的物理意義，各子塊的四行與系數(shù)向量的乘積分別對應(yīng)梁段的撓度、轉(zhuǎn)角、彎矩與剪力。隨著頻率的升高，撓度與轉(zhuǎn)角所對應(yīng)的行元素會與剪力、彎矩所對應(yīng)的行元素存在數(shù)量級的差距，導(dǎo)致矩陣條件數(shù)升高，影響方程對舍入誤差的敏感度。當(dāng)矩陣病態(tài)嚴(yán)重，系統(tǒng)特征方程曲線就會發(fā)散。因此，采用適當(dāng)?shù)膶羌訖?quán)陣W(ω)左乘T(ω)，以改善矩陣的病態(tài)問題。加權(quán)矩陣W(ω)可取為

(16)

式中，m為矩陣T的階數(shù)。對T(ω)左乘加權(quán)矩陣，可平衡各行的數(shù)量級差距，有效降低矩陣的條件數(shù)，減弱方程式(8)對數(shù)值誤差的敏感度，在保證求解數(shù)值穩(wěn)定性的同時(shí)，提高了計(jì)算頻率范圍。因此，本方法實(shí)際上為改進(jìn)的數(shù)值組裝法(modified NAM, m-NAM)。左乘W(ω)后，特征方程式變?yōu)?/p>

Δ=W(ω)T(ω)=0

(17)

為避免數(shù)值過大對頻率方程式求根帶來的不便，參考Wu等的研究中對CTMM系統(tǒng)特征方程的變換，對Δ做如下變換

(18)

對變換后的特征方程δ(ω)求根，即可求得系統(tǒng)固有頻率。

2 數(shù)值算例驗(yàn)證

2.1 算例1(等截面解析算例)

對于均質(zhì)等截面Timoshenko梁，在自由-自由邊界條件下，系統(tǒng)矩陣行列式的展開式為

(19)

式(19)即為系統(tǒng)特征方程，對其求根可得系統(tǒng)固有頻率的解析解。以一段長度10 m、直徑0.06 m的圓截面鋼梁為例，分別使用CTMM、NAM、m-NAM與解析表達(dá)式計(jì)算系統(tǒng)的特征方程曲線。將梁分為3段，長度分別為1 m，8 m，1 m，計(jì)算結(jié)果如圖1所示。

圖1 系統(tǒng)特征方程曲線

CTMM也是針對梁段微分方程的系數(shù)進(jìn)行求解，將式(8)改寫為

(20)

通過遞推式，得到第一段系數(shù)C1與第n段系數(shù)Cn的關(guān)系，并利用邊界條件構(gòu)建系統(tǒng)矩陣

(21)

方程式(21)即為CTMM法的系統(tǒng)特征方程，其中包含各梁段子塊累乘的部分(第三、第四行)。在數(shù)值仿真中發(fā)現(xiàn)，分段數(shù)與分段長度對CTMM系統(tǒng)矩陣數(shù)值大小沒有影響，而各段的舍入誤差不斷累積并放大，且系統(tǒng)矩陣Z病態(tài)嚴(yán)重，方程的解易在微小數(shù)值擾動下發(fā)散。

NAM中梁段的誤差沒有經(jīng)過累乘放大，且由于矩陣病態(tài)主要來源于矩陣元素的數(shù)量級差距。在梁段子塊中，隨著頻率升高，雙曲余弦、雙曲正弦元素與正弦、余弦元素的數(shù)值(一、二列與三、四列)量級差距非常大，在運(yùn)算過程中舍入誤差會對矩陣方程的求解造成擾動。從圖1中可以看出，CTMM法的特征方程曲線在200 Hz以上完全失去信息，無法判斷固有頻率，NAM的計(jì)算頻率范圍也止于350 Hz，而m-NAM能夠很好地保持?jǐn)?shù)值穩(wěn)定性。從圖2中的對比可以看出，m-NAM的前150階固有頻率與對式(19)求根得到的固有頻率解析解吻合良好，固有頻率的相對誤差在0.5/10 000以內(nèi)。

圖2 模態(tài)頻率計(jì)算結(jié)果及相對誤差

2.2 算例2(推進(jìn)軸系算例)

為進(jìn)一步驗(yàn)證m-NAM在軸系動力學(xué)特性計(jì)算上的準(zhǔn)確性， 使用有限元方法建立算例軸系模型。軸系總長15 m，采用實(shí)心圓截面，材料為鋼，楊氏模量E=2.1×1011Pa，密度ρ=7 850 kg/m3，泊松比ν=0.3。

表1 算例系統(tǒng)軸段長度與截面參數(shù)

軸系共有4個支承，分別軸承k1，軸承k2，軸承k3以及軸承k4。軸系由支承位置與截面變化共分為11段，共12個節(jié)點(diǎn)，4個支承布置于節(jié)點(diǎn)2，5，8，11處，如圖3所示。k1=5×108N/m，k2=k3=k4=1×109N/m，各支承置于剛性基礎(chǔ)上。各軸段對應(yīng)長度與截面半徑如表1所示?？紤]到較高的計(jì)算頻率范圍，設(shè)置單元尺寸10 mm，共計(jì)15 000個單元。

圖3 軸系模型

分別使用NAM與m-NAM繪制系統(tǒng)矩陣特征方程隨頻率變化的曲線，并使用二分法對其求根。從圖4中可以看出，對于算例軸系，NAM的頻率特征方程在1 300 Hz左右開始出現(xiàn)數(shù)值發(fā)散，而m-NAM可以保持很好的數(shù)值穩(wěn)定性。

圖4 軸系特征方程曲線

對比圖5中m-NAM與FEM的計(jì)算結(jié)果可知，使用m-NAM得到的前40階固有頻率與有限元計(jì)算結(jié)果的相對誤差小于1/1 000，表明m-NAM具有良好的計(jì)算精度與數(shù)值穩(wěn)定性。

要計(jì)算軸系強(qiáng)迫振動響應(yīng)，將第i節(jié)點(diǎn)的外力作用在方程式(8)的等號右側(cè)，則界面連續(xù)性條件可簡寫為

AcliC(i)=AcriC(i+1)+Fi

(22)

圖5 系統(tǒng)前40階橫向振動固有頻率及相對誤差

(23)

式中，F(xiàn)L與FR為邊界作用力，為2×1向量。將邊界作用力與內(nèi)部節(jié)點(diǎn)的外力組裝成長4n×1的向量，置于式(15 )等式右側(cè)，則NAM計(jì)算強(qiáng)迫振動響應(yīng)的矩陣方程為

T(ω)C(ω)=F(ω)

(24)

求解方程式(24)的系數(shù)向量C，即可重構(gòu)外激勵下各梁段的位移、轉(zhuǎn)角、彎矩與剪力信息，得到結(jié)構(gòu)強(qiáng)迫振動響應(yīng)。

3 結(jié) 論

本文給出改進(jìn)數(shù)值組裝方法對推進(jìn)軸系進(jìn)行動力學(xué)建模，利用Timoshenko梁理論和節(jié)點(diǎn)連續(xù)性條件構(gòu)建階梯軸的控制方程，考慮質(zhì)量、慣量、彈性元件等集中元素。數(shù)值算例表明，m-NAM避免了CTMM與NAM中高頻數(shù)值發(fā)散問題，計(jì)算頻率范圍更寬。

此外，m-NAM可以容易地推廣到軸系縱向、扭轉(zhuǎn)振動的固有特性和強(qiáng)迫響應(yīng)計(jì)算。

附錄A

左端邊界矩陣BL

右端邊界矩陣BR

其中