趙　琪，尹韶平，王　中，郭　君

（1. 中國(guó)船舶重工集團(tuán)公司 第705研究所，陜西 西安，710075； 2. 水下信息與控制重點(diǎn)實(shí)驗(yàn)室，陜西 西安，710075）

基于MATLAB的魚雷推進(jìn)軸系彎曲振動(dòng)渦動(dòng)頻率計(jì)算

趙琪1，2，尹韶平1，2，王中1，2，郭君1

（1. 中國(guó)船舶重工集團(tuán)公司 第705研究所，陜西 西安，710075； 2. 水下信息與控制重點(diǎn)實(shí)驗(yàn)室，陜西 西安，710075）

在魚雷轉(zhuǎn)子系統(tǒng)中，為避免系統(tǒng)發(fā)生共振，常需要計(jì)算其渦動(dòng)頻率。鑒于商業(yè)有限元分析軟件的輸入?yún)?shù)有限，且無法滿足特殊情況下的計(jì)算需求，該文根據(jù)轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)微分方程及其求解過程，運(yùn)用MATLAB軟件編程來計(jì)算系統(tǒng)的渦動(dòng)頻率，然后與ANSYS軟件計(jì)算結(jié)果進(jìn)行對(duì)比。結(jié)果表明，所編寫的MATLAB程序合理，可對(duì)魚雷推進(jìn)軸系進(jìn)行計(jì)算，也可進(jìn)行后續(xù)開發(fā)。

魚雷推進(jìn)軸系； 彎曲振動(dòng)； 渦動(dòng)頻率

0　引言

一般工業(yè)中常見的機(jī)器都裝有旋轉(zhuǎn)部件即轉(zhuǎn)子。 轉(zhuǎn)子連同它的軸承和支座等統(tǒng)稱為轉(zhuǎn)子系統(tǒng)。轉(zhuǎn)子系統(tǒng)的振動(dòng)是多樣的，包括轉(zhuǎn)軸的縱向振動(dòng)、扭轉(zhuǎn)振動(dòng)和彎曲振動(dòng)等。其中轉(zhuǎn)軸的彎曲振動(dòng)較為復(fù)雜，牽扯的因素較多。該文就是以轉(zhuǎn)軸的彎曲振動(dòng)作為研究對(duì)象。

求解系統(tǒng)的渦動(dòng)頻率及臨界轉(zhuǎn)速，是轉(zhuǎn)子動(dòng)力學(xué)研究的重要內(nèi)容之一［1］。當(dāng)轉(zhuǎn)子達(dá)到臨界轉(zhuǎn)速時(shí)，往往產(chǎn)生劇烈振動(dòng)，為了避免這種現(xiàn)象，就必須計(jì)算轉(zhuǎn)子的渦動(dòng)速度。計(jì)算方法有傳遞矩陣法和有限元法等。采用有限元法建立的模型能考慮的因素比傳遞矩陣多，且數(shù)值穩(wěn)定，計(jì)算精度較高，因此得到了非常廣泛的應(yīng)用。現(xiàn)在應(yīng)用較多的有限元分析軟件，如ANSYS，ABAQUS等，雖然都可以進(jìn)行軸系振動(dòng)特性分析，但其輸入?yún)?shù)有限（如彈簧的剛度阻尼系數(shù)等），而實(shí)際情況又都較為復(fù)雜，存在彎、扭、縱的耦合振動(dòng)，需要考慮更加全面的邊界條件； 另一方面，當(dāng)商業(yè)軟件不能滿足特定情況的計(jì)算需求時(shí)，要對(duì)其進(jìn)行更多功能的開發(fā)較為困難。因此該文采用MATLAB編程計(jì)算和分析，實(shí)現(xiàn)對(duì)輸入?yún)?shù)的控制，也為后續(xù)其他功能的實(shí)現(xiàn)打好基礎(chǔ)。

1　有限元分析理論基礎(chǔ)

1.1渦動(dòng)頻率

通常轉(zhuǎn)軸的兩支點(diǎn)在同一水平線上，中心線是水平的。當(dāng)圓盤受橫向力作用，圓盤或轉(zhuǎn)軸的中心o′在互相垂直的2個(gè)方向作頻率同為wn的簡(jiǎn)諧運(yùn)動(dòng)。一般情況下，中心o′的軌跡為一橢圓，如圖1所示。o′的這種運(yùn)動(dòng)是一種“渦動(dòng)”或稱“進(jìn)動(dòng)”，自然頻率wn稱為渦動(dòng)角速度［2］。

圖1　轉(zhuǎn)子渦動(dòng)示意圖Fig. 1　Schematic of rotor whirling

1.2運(yùn)動(dòng)微分方程

轉(zhuǎn)子系統(tǒng)通常由剛性圓盤、彈性軸段和軸承座等部件組成，因此可沿軸線將轉(zhuǎn)子系統(tǒng)劃分為圓盤、軸段和軸承座等單元，各單元間彼此在結(jié)點(diǎn)處聯(lián)結(jié)。結(jié)點(diǎn)通常選在圓盤中心、軸頸中心以及軸的變截面等位置。如以軸承座中心線為s軸建立坐標(biāo)，則轉(zhuǎn)子系統(tǒng)的具體結(jié)構(gòu)見圖2［3］。

圖2　轉(zhuǎn)子系統(tǒng)結(jié)構(gòu)Fig. 2 Structure of a rotor system

基于上述模型，建立轉(zhuǎn)子系統(tǒng)運(yùn)動(dòng)微分方程

式中： ｛U｝，｛Q｝分別為系統(tǒng)位移向量及廣義力；［M］，［C］和［K］分別為系統(tǒng)質(zhì)量矩陣、回轉(zhuǎn)矩陣和剛度矩陣。系統(tǒng)矩陣是由單元矩陣根據(jù)對(duì)應(yīng)結(jié)點(diǎn)疊加而成的。

1.3結(jié)點(diǎn)為4自由度的單元矩陣形式

若考慮每個(gè)結(jié)點(diǎn)為4自由度狀態(tài)，則任一截面的位移可表示為

設(shè)剛性圓盤的質(zhì)量、直徑轉(zhuǎn)動(dòng)慣量和極轉(zhuǎn)動(dòng)慣量分別為m，Jd和Jp，則圓盤的質(zhì)量矩陣［Md］和回轉(zhuǎn)矩陣［Gd］為

設(shè)軸段單元長(zhǎng)度為l，半徑為r，單位長(zhǎng)度質(zhì)量為μ。則其移動(dòng)慣性矩陣［MsT］、轉(zhuǎn)動(dòng)慣性矩陣［MsR］和回轉(zhuǎn)矩陣［Gs］分別為［4-5］

可得軸段單元的質(zhì)量矩陣為

軸段單元?jiǎng)偠染仃嚍?/p>

根據(jù)各單元間的相互作用關(guān)系，將單元矩陣進(jìn)行集總，形成系統(tǒng)整體矩陣。

1.4結(jié)點(diǎn)為6自由度的單元矩陣形式

若考慮每個(gè)結(jié)點(diǎn)為6自由度狀態(tài)，圓盤單元只有一個(gè)結(jié)點(diǎn)，其位移向量可表示為

其質(zhì)量矩陣及回轉(zhuǎn)矩陣可表示為

其中： mx，my和mz為質(zhì)量單元3個(gè)方向的質(zhì)量； Jx為極轉(zhuǎn)動(dòng)慣量； Jy，Jz為直徑轉(zhuǎn)動(dòng)慣量。

軸段單元為兩結(jié)點(diǎn)單元，每個(gè)結(jié)點(diǎn)均有6個(gè)自由度，其單元的剛度矩陣為［6］

式中： A為橫截面面積； Ixx是橫截面對(duì)x軸的極慣性矩； Iyy和Izz是單元截面對(duì)y和z軸的主慣性矩。

軸段單元的一致質(zhì)量矩陣為［7-8］

回轉(zhuǎn)矩陣為

將上述單元矩陣疊加，即可得到系統(tǒng)矩陣。

2　運(yùn)動(dòng)微分方程的求解

系統(tǒng)運(yùn)動(dòng)微分方程的齊次式為

在結(jié)構(gòu)動(dòng)力學(xué)中，求解這種大型稀疏帶狀矩陣特征值問題的程序，通常要求在對(duì)稱矩陣情況下使用。而在轉(zhuǎn)子動(dòng)力學(xué)中，由于軸承動(dòng)力特性參數(shù)中包括交叉剛度和交叉阻尼，且二者往往并不相等，使得剛度矩陣［K］失去了對(duì)稱性，因此一般采用復(fù)模態(tài)理論求解轉(zhuǎn)子動(dòng)力學(xué)方程。

復(fù)模態(tài)理論是將運(yùn)動(dòng)方程的齊次式化為1階微分方程組來求解，即令

則運(yùn)動(dòng)方程的齊次式可改寫為

帶入式（17）可得

由此可求解特征值v及特征向量｛V0｝。特征值v即為系統(tǒng)在給定自轉(zhuǎn)角速度下的渦動(dòng)頻率，且根據(jù)特征向量｛V0｝可得相應(yīng)的模態(tài)振型。

3　MATLAB編程

3.1編程流程

根據(jù)上述對(duì)系統(tǒng)運(yùn)動(dòng)微分方程及其求解過程的分析，編寫MATLAB程序。

首先對(duì)轉(zhuǎn)子系統(tǒng)進(jìn)行結(jié)點(diǎn)單元的劃分，以建立有限元模型。其次計(jì)算結(jié)點(diǎn)處的單元矩陣，并將單元矩陣按照結(jié)點(diǎn)位置進(jìn)行疊加，形成系統(tǒng)整體矩陣，得到微分方程中的矩陣參數(shù)。最后根據(jù)支承情況，對(duì)系統(tǒng)剛度矩陣進(jìn)行修正處理。在求解過程中，根據(jù)復(fù)模態(tài)理論對(duì)系統(tǒng)矩陣進(jìn)行擴(kuò)充，將運(yùn)動(dòng)方程轉(zhuǎn)化成1階微分方程組，然后可直接利用MATLAB中eig程序（求特征值和特征向量）進(jìn)行求解。

程序具體流程如圖3所示。

圖3　編程流程圖Fig. 3　Flowchart of programming

按照結(jié)點(diǎn)處為4自由度和6自由度狀態(tài)分別編程。4自由度和6自由度程序的單元疊加方式一致，但在對(duì)剛度矩陣進(jìn)行修正時(shí)會(huì)有所差異。

3.24自由度剛度矩陣修正

對(duì)于一般的軸承單元，其剛度、阻尼系數(shù)矩陣分別為

若假設(shè)系統(tǒng)具有N個(gè)結(jié)點(diǎn)，其間用N-1個(gè)軸段連接而成，則系統(tǒng)的位移向量由式（2）改寫為

考慮支承處的軸承作用，則此時(shí)系統(tǒng)運(yùn)動(dòng)方程為

假設(shè)第j個(gè)結(jié)點(diǎn)為軸承單元，則在2N×2N階的矩陣［c11］，［c21］，［c12］，［c22］和［k11］，［k12］，［k21］，［k22］中，除了第2s（j）-1行和2s（j）-1列中有元素cxx，cxy，cyx，cyy，kxx，kxy，kyx，kyy外，其余元素都為零。

則運(yùn)動(dòng)微分方程的各個(gè)參數(shù)為如下形式

3.36自由度剛度矩陣修正

考慮軸承單元6個(gè)方向的連接剛度，則其剛度系數(shù)矩陣形式為

在對(duì)系統(tǒng)剛度矩陣進(jìn)行修正時(shí)，可將軸承的剛度系數(shù)矩陣按照其結(jié)點(diǎn)位置直接疊加到系統(tǒng)矩陣中進(jìn)行計(jì)算。

4　算例分析

ANSYS幫助文檔中VM254算例［9］的模型為典型的轉(zhuǎn)子動(dòng)力學(xué)模型。它由剛性圓盤、變截面的彈性軸段及滑動(dòng)軸承組成，具體結(jié)構(gòu)見圖4。在變截面處、支承位置和圓盤位置設(shè)置結(jié)點(diǎn)，則可沿軸線自左至右將軸段劃分為18段，總計(jì)19個(gè)結(jié)點(diǎn)。模型具體參數(shù)可參考ANSYS幫助文檔。

ANSYS運(yùn)行結(jié)果的坎貝爾圖如圖5所示。圖中交點(diǎn)的橫坐標(biāo)即為臨界轉(zhuǎn)速值。應(yīng)用所編寫的MATLAB程序?qū)ι鲜瞿Ｐ瓦M(jìn)行計(jì)算。按4自由度狀態(tài)得出的坎貝爾圖見圖6。

按6自由度狀態(tài)得出的坎貝爾圖見圖7。

圖4　轉(zhuǎn)子模型Fig. 4　Rotor model

圖5　ANSYS得出的坎貝爾圖Fig. 5　Campbell diagram from ANSYS

圖7　6自由度狀態(tài)的坎貝爾圖Fig. 7　Campbell diagram of six-degree-of-freedom state

在輸入轉(zhuǎn)速為105r/min時(shí)的渦動(dòng)頻率計(jì)算結(jié)果如表1所示。

表1　渦動(dòng)頻率計(jì)算結(jié)果Table 1　Calculation results of whirling frequency

對(duì)比三者間的坎貝爾圖可以看出，6自由度程序得出的結(jié)果與ANSYS比較吻合。通過分析在轉(zhuǎn)速為105r/min時(shí)的計(jì)算結(jié)果可以看出，6自由度程序所得結(jié)果較4自由度程序更接近ANSYS的計(jì)算結(jié)果。所以，可以認(rèn)為6自由度程序是合理的，可以用來計(jì)算魚雷推進(jìn)軸系的渦動(dòng)頻率。

5　6自由度程序在魚雷推進(jìn)軸系的應(yīng)用

5.1魚雷推進(jìn)軸系結(jié)構(gòu)及其有限元模型

魚雷推進(jìn)軸系結(jié)構(gòu)如圖8所示。動(dòng)力裝置與花鍵軸之間、花鍵軸與尾軸之間通過花鍵相連，推進(jìn)器轉(zhuǎn)子安裝在尾軸上?；ㄦI軸及尾軸主要用于將動(dòng)力裝置的輸出功率傳遞給推進(jìn)器，而推進(jìn)器工作時(shí)產(chǎn)生的推力通過推力軸承傳遞至雷尾殼體，推動(dòng)魚雷航行。在花鍵軸與尾軸相連的一端裝有彈性花鍵套（圖中標(biāo)記2）。尾軸的支撐是靠隔板的滾珠軸承和尾蓋處的滑動(dòng)軸承，在滾珠軸承及滑動(dòng)軸承外側(cè)都安裝了金屬橡膠隔振器（圖中標(biāo)記5和9）。

魚雷推進(jìn)軸系簡(jiǎn)化模型（將推進(jìn)器作為質(zhì)量單元處理）如圖9所示。

建立魚雷推進(jìn)軸系的有限元模型，可以沿軸線將花鍵軸和尾軸劃分為若干個(gè)結(jié)點(diǎn)，如圖10所示。對(duì)于彈性花鍵套，可以采用等效軸段法模擬其橫向剛度和轉(zhuǎn)角剛度。圖中，結(jié)點(diǎn)1為花鍵軸前端，與動(dòng)力裝置相連； 結(jié)點(diǎn)2為花鍵套簡(jiǎn)化位置； 結(jié)點(diǎn)8，9和10為彈性花鍵套簡(jiǎn)化位置； 結(jié)點(diǎn)11和19為軸承支撐位置； 結(jié)點(diǎn)21為推進(jìn)器安裝位置。

圖8　魚雷推進(jìn)軸系結(jié)構(gòu)Fig. 8　Structure of torpedo propulsion shafting

圖9　推進(jìn)軸系簡(jiǎn)化模型Fig. 9　Simplified model of torpedo propulsion shafting

圖10　花鍵軸和尾軸的結(jié)點(diǎn)劃分Fig. 10　Elements dividing of spline-shaft and tail-shaft

5.2計(jì)算結(jié)果與分析

應(yīng)用6自由度MATLAB程序?qū)︳~雷推進(jìn)軸系進(jìn)行渦動(dòng)頻率計(jì)算，得出的坎貝爾圖見圖11。

選取轉(zhuǎn)速1150 r/min，1500 r/min和2050 r/min時(shí)的前4階渦動(dòng)頻率如表2所示。利用ANSYS軟件進(jìn)行計(jì)算的結(jié)果見表3。對(duì)比推進(jìn)軸系不同轉(zhuǎn)速下的前4階渦動(dòng)頻率可得，MATLAB程序計(jì)算結(jié)果與ANSYS計(jì)算結(jié)果基本一致，證明該文所編程序是合理的。

6　結(jié)束語

根據(jù)轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)微分方程，從2個(gè)方面進(jìn)行編程，分別計(jì)算了結(jié)點(diǎn)為4自由度狀態(tài)和6自由度狀態(tài)時(shí)的結(jié)果。6自由度狀態(tài)考慮了結(jié)點(diǎn)在3個(gè)方向的位移和轉(zhuǎn)角，并在進(jìn)行系統(tǒng)剛度矩陣修正時(shí)，對(duì)軸承剛度和阻尼系數(shù)的處理方式也比較完整，更接近于實(shí)際情況。

基于此程序，可以對(duì)轉(zhuǎn)子系統(tǒng)在其他實(shí)際狀態(tài)下進(jìn)行渦動(dòng)頻率及臨界轉(zhuǎn)速等因素的分析，也可進(jìn)行后續(xù)的開發(fā)。

圖11　推進(jìn)軸系坎貝爾圖Fig. 11　Campbell diagram of propulsion shafting

表2　MATLAB計(jì)算結(jié)果Table 2　Calculation results of MATLAB

表3　ANSYS計(jì)算結(jié)果Table 3　Calculation results of ANSYS

［1］阮小麗. 基于傳遞矩陣法和有限元法的轉(zhuǎn)子動(dòng)力學(xué)分析［J］. 機(jī)電工程技術(shù)，2011，40（3）： 71-73.

Ruan Xiao-li. Based on Transfer Matrix Method and the Finite Element Method Analysis of Rotor Dynamics［J］. Mechanical & Electrical Engineering Technology. 2011，40（3）： 71-73.

［2］鐘一諤. 轉(zhuǎn)子動(dòng)力學(xué)［M］. 北京： 清華大學(xué)出版社，1984.

［3］王軍峰，孫康. 基于有限元法的轉(zhuǎn)子臨界轉(zhuǎn)速計(jì)算［J］.機(jī)械設(shè)計(jì). 2012，29（12）： 10-13.

Wang Jun-feng，Sun Kang. Rotor Critical Speed Calculation Based on Finite Element Method［J］. Journal of Machine Design.2012，29（12）： 10-13.

［4］汪夢(mèng)甫，王朝暉. 兩種質(zhì)量矩陣在梁模態(tài)分析中差異的比較［J］. 地震工程與工程振動(dòng). 2006，26（6）： 84-86.

Wang Meng-fu，Wang Zhao-hui. Analysis of Differences Between Two Types of Mass Matrixes in Beam Modal Analysis［J］. Earthquake Engineering and Engineering Vibration. 2006，26（6）： 84-86.

［5］徐榮橋. 結(jié)構(gòu)分析的有限元法與MATLAB程序設(shè)計(jì)［M］. 北京： 人民交通出版社，2005.

［6］車樹汶，陳權(quán)，樓松慶. 質(zhì)量矩陣模式對(duì)橋梁自振頻率的影響［J］. 蘭州交通大學(xué)學(xué)報(bào)，2003，22（6）： 80-83.

Che Shu-wen，Chen Quan，Lou Song-qing. Influence of Mass Matrix Model on Natural Vibration Frequency of Bridges［J］. Journal of Lanzhou Jiaotong University，2003，22（6）： 80-83.

［7］董軍，劉旭紅，姚順忠，等. 基于虛功原理表達(dá)的空間梁?jiǎn)卧獎(jiǎng)偠染仃嚪治觯跩］. 西南林學(xué)院學(xué)報(bào)，2002，22（4）：59-62.

Dong Jun，Liu Xu-hong，Yao Shun-zhong，et al. Analysis of Element Stiffness Matrix About Space Beam-Element Based on Principle of Virtual Work［J］. Journal of Southwest Forestry College，2002，22（4）： 59-62.

［8］王勖成. 有限單元法［M］. 北京： 清華大學(xué)出版社，2003.

［9］Vaugh N M. The Dynamics of Rotor-Bearing Systems Using Finite Element［J］. Journal of Engineering for Industry，1976，98（2）：593-600.

（責(zé)任編輯： 陳曦）

Calculation of Whirling Frequency in Flexural Vibration of Torpedo Propulsion Shafting Based on MATLAB

ZHAO Qi1，2，YIN Shao-ping1，2，WANG Zhong1，2，GUO Jun1
（1. The 705 Research Institute，China Shipbuilding Industry Corporation，Xi′an 710075，China； 2. Science and Technology on Underwater Information and Control Laboratory，Xi′an 710075，China）

To avoid resonance，it is necessary to calculate whirling frequency of a rotor system. The input parameters of commercial finite element analysis（FEA）softwares are limited，and these softwares can not complete the analysis of special cases. In this paper，according to the differential equation of motion for a rotor system and its solution process，a program is coded to calculate the whirling frequency by the software MATLAB. Comparison between the results of MATLAB and ANSYS shows that this MATLAB program is reasonable，and it can be used to calculate whirling frequency in flexural vibration of torpedo propulsion shafting or to perform subsequent development.

torpedo propulsion shafting； flexural vibration； whirling frequency

TJ630.1

A

1673-1948（2015）01-0007-07

2014-05-14；

2014-07-10.

趙琪（1989-），男，在讀碩士，研究方向?yàn)轸~雷總體技術(shù).