于 翔，周 松

(沈陽(yáng)航空航天大學(xué) 機(jī)電工程學(xué)院，遼寧 沈陽(yáng) 110136)

基于MATLAB多自由度單向串聯(lián)振動(dòng)系統(tǒng)的計(jì)算分析

于 翔，周 松

(沈陽(yáng)航空航天大學(xué) 機(jī)電工程學(xué)院，遼寧 沈陽(yáng) 110136)

在工程振動(dòng)問(wèn)題中，固有特性(固有頻率和主振型)是描述振動(dòng)系統(tǒng)的主要參數(shù)。根據(jù)矩陣迭代法，利用MATLAB/GUI界面設(shè)計(jì)平臺(tái)，提出了針對(duì)多自由度單向串聯(lián)振動(dòng)系統(tǒng)固有特性的計(jì)算方法，并通過(guò)算例求解和誤差分析，完成了合理性驗(yàn)證，為復(fù)雜振動(dòng)系統(tǒng)響應(yīng)的初步解耦奠定了基礎(chǔ)。

MATLAB；多自由度；振動(dòng)系統(tǒng)；矩陣迭代法；固有特性

振動(dòng)力學(xué)不僅是近代應(yīng)用力學(xué)的一個(gè)重要分支，而且是現(xiàn)今工業(yè)生產(chǎn)中分析和解決機(jī)械產(chǎn)品運(yùn)轉(zhuǎn)穩(wěn)定和壽命持久問(wèn)題的重要手段。隨著科學(xué)技術(shù)的高速發(fā)展，振動(dòng)力學(xué)在機(jī)械加工、航空航天以及交通運(yùn)輸?shù)裙I(yè)技術(shù)領(lǐng)域中占有愈來(lái)愈重要的地位[1]。其中，固有特性(固有頻率和主振型)作為描述振動(dòng)系統(tǒng)的主要參數(shù)[2-5]，在響應(yīng)模擬和仿真預(yù)測(cè)中發(fā)揮著重要的作用。因此，通過(guò)研究振動(dòng)系統(tǒng)固有特性的求解方法，可以為防止系統(tǒng)共振提供理論依據(jù)[6-8]，也為進(jìn)行系統(tǒng)響應(yīng)的初步解耦和機(jī)械結(jié)構(gòu)的強(qiáng)度分析奠定基礎(chǔ)[9]。然而，對(duì)于工程振動(dòng)問(wèn)題，無(wú)論是復(fù)雜的運(yùn)動(dòng)體還是離散的彈性體都可以通過(guò)有限元法化解為理想節(jié)點(diǎn)系統(tǒng)的振動(dòng)模型，但是由于系統(tǒng)自由度不唯一，并且以多自由度單向串聯(lián)振動(dòng)系統(tǒng)最為常見(jiàn)，因此如何準(zhǔn)確地求得振動(dòng)系統(tǒng)任意階固有特性的有效解，一直都是研究振動(dòng)力學(xué)問(wèn)題的關(guān)鍵。

1 固有特性的求解方法

多自由度振動(dòng)系統(tǒng)固有特性的算法主要有直接計(jì)算和近似計(jì)算2種。其中，對(duì)低自由度振動(dòng)系統(tǒng)，直接計(jì)算作為首選方法，既簡(jiǎn)便快捷又準(zhǔn)確有效，但是在實(shí)際工程中，大多數(shù)系統(tǒng)的振動(dòng)形式非常復(fù)雜且自由度很高，用頻率方程(特征方程)很難進(jìn)行直接求解。然而，通過(guò)相關(guān)研究發(fā)現(xiàn)，近似計(jì)算中的矩陣迭代法因?yàn)槠渚_度不依賴于假設(shè)振型，假設(shè)初值的好壞只影響迭代次數(shù)，與所求固有特性的精度無(wú)關(guān)，即迭代運(yùn)算總是收斂于振動(dòng)系統(tǒng)的最低階固有特性[1,10]。因此，對(duì)于實(shí)際振動(dòng)問(wèn)題，往往采用矩陣迭代的方法求解系統(tǒng)的固有特性。

在多自由度正定系統(tǒng)的自由振動(dòng)中，系統(tǒng)的主振型方程可分別表示為：

ω2{A}=[M]-1[K]{A}

(1)

(2)

式中：ω為固有頻率；{A}為主振型；[M]為質(zhì)量矩陣；[K]為剛度矩陣；[δ]為柔度矩陣。矩陣迭代法是從假設(shè)主振型出發(fā)，對(duì)以上兩式進(jìn)行矩陣迭代運(yùn)算。因?yàn)楣こ躺蠈?duì)系統(tǒng)的最低階或較低階固有特性比較重視[2-3]，所以通常只對(duì)式(2)進(jìn)行迭代運(yùn)算。引進(jìn)動(dòng)力矩陣：

[D]=[δ][M]

(3)

通過(guò)矩陣迭代，在規(guī)定的有效位數(shù)內(nèi)，當(dāng)發(fā)現(xiàn){A}k≈{A}k-1時(shí)，運(yùn)算結(jié)束。此時(shí)，{A}k或{A}k-1即為系統(tǒng)第一階主振型{A(1)}的近似值；系數(shù)ak即為系統(tǒng)第一階固有頻率平方倒數(shù)的近似值，即：

{A(1)}≈{A}k-1={A}k

(4)

(5)

利用式(2)不僅可以求出最低階固有特性，而且還能依次求得較低階固有特性。其方法是在求系統(tǒng)第(s+1)階固有特性時(shí)，每次迭代前都乘以清型變換后的動(dòng)力矩陣[D]s。根據(jù)式(2)，則有：

[D]s=[D][Q]s=

(6)

式中：ωj為第j階固有頻率；{A(j)}為第j階主振型；Mj為系統(tǒng)中第j個(gè)物體的質(zhì)量。通過(guò)上述方法得到清型變換后的動(dòng)力矩陣[D]s。然后，根據(jù)第一階固有特性的算法求出第(s+1)階固有特性，并以此類推求解系統(tǒng)的任意階固有特性，其對(duì)應(yīng)的程序流程圖如圖1所示。

圖1 矩陣迭代法求解固有特性的程序流程圖

2 程序設(shè)計(jì)

MATLAB 作為一款數(shù)據(jù)分析軟件，對(duì)矩陣迭代運(yùn)算尤其合適[3，11]。通過(guò)GUI界面設(shè)計(jì)，能為用戶提供人性化的操作界面[12]，便于針對(duì)多自由度振動(dòng)系統(tǒng)的計(jì)算分析。

2.1GUI界面操作

根據(jù)矩陣迭代法，基于MATLAB軟件設(shè)計(jì)的計(jì)算多自由度單向串聯(lián)振動(dòng)系統(tǒng)固有特性的GUI操作界面如圖2所示。通過(guò)輸入各單元的質(zhì)量和剛度以及振動(dòng)系統(tǒng)的初始假設(shè)振型，就可以自動(dòng)生成系統(tǒng)本身固有的質(zhì)量矩陣和剛度矩陣。同時(shí)，可以根據(jù)用戶需要，求出系統(tǒng)的任意階固有特性。

圖2 軟件操作界面

2.2MATLAB程序

根據(jù)矩陣迭代法計(jì)算系統(tǒng)固有特性的關(guān)鍵在于能否求出動(dòng)力矩陣，而動(dòng)力矩陣是由質(zhì)量矩陣和柔度矩陣進(jìn)行矩陣乘法運(yùn)算得到的[2]。對(duì)于多自由度單向串聯(lián)振動(dòng)系統(tǒng)，因?yàn)槿岫染仃囀枪战蔷仃嘯12]，而拐角矩陣具有沿主對(duì)角線分層的特點(diǎn)，如圖3所示，所以在程序語(yǔ)言的編寫過(guò)程中，根據(jù)單層循環(huán)結(jié)構(gòu)很難完成矩陣建立。

圖3 柔度矩陣的特有形式拐角矩陣

通過(guò)對(duì)多自由度單向串聯(lián)振動(dòng)系統(tǒng)的研究發(fā)現(xiàn)，柔度矩陣主對(duì)角線上元素為各自由度上剛度倒數(shù)值的累加和，因此通過(guò)雙層循環(huán)嵌套和矩陣元素的行列轉(zhuǎn)換，可以準(zhǔn)確建立任意階數(shù)的柔度矩陣。其程序指令如下：

i=1;w=0;

while i<=n %層間循環(huán)

j=i;w=w+1/k(i,1); %剛度倒數(shù)求和

while j<=n %層內(nèi)循環(huán)

K0(j,i)=w; K0(i,j)=w; %主對(duì)角線對(duì)稱元素賦值

j=j+1;

end

i=i+1;

end

由于矩陣迭代法的主要計(jì)算方式就是歸一化后根據(jù)程序設(shè)計(jì)精度反復(fù)求積和不斷賦值，所以在MATLAB程序語(yǔ)言的編寫過(guò)程中，通過(guò)矩陣中列與列的對(duì)比、運(yùn)算和變換，可以在提高運(yùn)算效率的同時(shí)減少計(jì)算空間，避免程序語(yǔ)句過(guò)于復(fù)雜而引起的變量重復(fù)調(diào)用。其對(duì)應(yīng)程序段如下：

%初始假設(shè)振型歸一化

a=A0(end,1);A0(:,1)=A0(:,1)./a;

D=K0*M0; %初始動(dòng)力矩陣

t=1;

while t<=s %求解第s階固有頻率和主振型

A(:,1)=A0(:,1); %初始假設(shè)振型賦值

A(:,2)=D*A(:,1);a=A(end,2);A(:,2)=A(:,2)./a;

%矩陣迭代法

while abs(sum(A(:,1)-A(:,2)))>=0.000001

A(:,1)=A(:,2);A(:,2)=D*A(:,1);

a=A(end,2);A(:,2)=A(:,2)./a;

end

M=A(:,2)'*M0*A(:,2); %第s階主質(zhì)量

%清型變換后的動(dòng)力矩陣

D=D-(A(:,2)*A(:,2)'*M0).*a./M;

t=t+1;

end

a=sqrt(1/a);Ai=A(:,2); %第s階固有特性

3 算例分析

在如圖4所示的三自由度簡(jiǎn)化振動(dòng)系統(tǒng)中，假設(shè)不考慮單位之間的相互換算，即各振動(dòng)單元的剛度k和質(zhì)量m均為1，并且各單元在起始時(shí)刻都處于零位[3]。根據(jù)MATLAB程序獲得的系統(tǒng)第三階固有特性見(jiàn)表1～3。

圖4 三自由度簡(jiǎn)化振動(dòng)系統(tǒng)

表1 設(shè)計(jì)精度為10-3時(shí)的計(jì)算值與迭代值對(duì)比表

表2 設(shè)計(jì)精度為10-4時(shí)的計(jì)算值與迭代值對(duì)比表

表3 設(shè)計(jì)精度為10-5時(shí)的計(jì)算值與迭代值對(duì)比表

通過(guò)以上數(shù)據(jù)表的對(duì)比可以發(fā)現(xiàn)，矩陣迭代法對(duì)于計(jì)算多自由度單向串聯(lián)振動(dòng)系統(tǒng)固有頻率的影響很小，但是對(duì)于主振型的影響很大。隨著系統(tǒng)階數(shù)的不斷提高，誤差呈現(xiàn)出逐級(jí)放大的趨勢(shì)。這可能是因?yàn)樵诰仃嚨ㄖ?，無(wú)論求解第幾階主振型都需要其對(duì)應(yīng)的上一階主振型參與計(jì)算，從而造成誤差累積。因此為了提高計(jì)算精度，需將程序的設(shè)計(jì)精度至少提高兩個(gè)數(shù)量級(jí)，例如，三自由度振動(dòng)系統(tǒng)的設(shè)計(jì)精度應(yīng)為10-5。但是，這在一定程度上增加了計(jì)算量，不過(guò)現(xiàn)有的高性能計(jì)算機(jī)可以很好地彌補(bǔ)矩陣迭代法的不足并發(fā)揮其長(zhǎng)處。

根據(jù)MATLAB/GUI軟件程序獲得的三自由度振動(dòng)系統(tǒng)的各單元振動(dòng)圖像如圖5所示。通過(guò)觀察可以發(fā)現(xiàn)，由于各單元之間的相互影響，運(yùn)動(dòng)形式不再是正弦函數(shù)，并且單元排列階數(shù)越高，其波形越不穩(wěn)定，如圖5中第三階主振型(單元M3密集虛線)所示，從而驗(yàn)證了隨自由度增加系統(tǒng)復(fù)雜性提高的結(jié)論[7，9]。

圖5 各單元振動(dòng)模擬圖像

4 結(jié)束語(yǔ)

本文采用MATLAB軟件，根據(jù)矩陣迭代法，計(jì)算出多自由度單向串聯(lián)振動(dòng)系統(tǒng)的固有特性，提出了針對(duì)特殊矩陣以及系統(tǒng)任意階固有特性的編程方法。通過(guò)對(duì)三自由度振動(dòng)系統(tǒng)的算例分析，對(duì)比了直接計(jì)算法和矩陣迭代法的區(qū)別以及各單元的模擬振動(dòng)響應(yīng)曲線。計(jì)算結(jié)果表明，隨著系統(tǒng)階數(shù)的提高，矩陣迭代法的誤差呈現(xiàn)出逐級(jí)放大的趨勢(shì)，但是通過(guò)提高設(shè)計(jì)精度，可以有效降低誤差。同時(shí)，根據(jù)圖像分析可以發(fā)現(xiàn)，隨著單元階數(shù)等級(jí)的提高，其對(duì)應(yīng)單元的波形越不穩(wěn)定，從而驗(yàn)證了隨自由度增加系統(tǒng)復(fù)雜性提高的結(jié)論，為設(shè)計(jì)人員分析多自由度系統(tǒng)自振特性及受迫振動(dòng)提供了參考依據(jù)，具有一定的實(shí)際應(yīng)用價(jià)值。

[1] 劉習(xí)軍. 工程振動(dòng)與測(cè)試技術(shù)[M]. 天津：天津大學(xué)出版社，1999：1-57.

[2] 文濤，胡青春. 基于MATLAB語(yǔ)言的多自由度振動(dòng)系統(tǒng)的固有頻率及主振型計(jì)算分析[J]. 中國(guó)制造業(yè)信息化，2007，36(1)：2-19.

[3] 潘延亮，董大偉，閆兵，等. 發(fā)動(dòng)機(jī)隔振系統(tǒng)振動(dòng)固有特性的優(yōu)化計(jì)算及影響因素分析[J]. 車用發(fā)動(dòng)機(jī)，2012(2)：17-25.

[4] 周斌. 內(nèi)燃機(jī)彈性基礎(chǔ)隔振系統(tǒng)固有頻率偏差分析[J]. 西南交通大學(xué)學(xué)報(bào)，1988，33(5)：503-507.

[5] 王秀，劉初升，彭利平，等. 含橫向?qū)еФ薊uler-Bernoulli梁彎曲振動(dòng)的固有特性分析[J]. 機(jī)械設(shè)計(jì)，2014，31(2)：81-85.

[6] 張鵬，安子軍. 擺線鋼球行星傳動(dòng)動(dòng)力學(xué)建模與固有特性分析[J]. 中國(guó)機(jī)械工程，2014，25(2)：157-162.

[7] 趙克剛，李健，姚偉浩. 車輛傳動(dòng)系統(tǒng)固有特性的識(shí)別[J]. 重慶理工大學(xué)學(xué)報(bào)：自然科學(xué)，2014，28(4)：5-9.

[8] 唐貴基，陳秀娟. 旋轉(zhuǎn)機(jī)械軸系扭振固有特性分析[J]. 汽輪機(jī)技術(shù)，2010，52(5)：339-344.

[9] 張憲，何洋，鐘江，等. 疲勞振動(dòng)試驗(yàn)臺(tái)的模態(tài)與諧響應(yīng)分析[J]. 機(jī)械設(shè)計(jì)與制造，2008(4)：12-14.

[10] 黃廷祝. 特殊矩陣分析及應(yīng)用[M]. 北京：科學(xué)出版社，2007：89-95.

[11] Balakumar B，Jame H D. MATLAB原理與工程應(yīng)用[M].2版. 北京：電子工業(yè)出版社，2006：292-341.

[12] 邱金蕙，王矞輝，李振全. 基于MATLAB/GUI的新型界面開(kāi)發(fā)方式[J]. 河北工業(yè)科技，2008，25(4)：16-23.

InherentCharacteristicsoftheOne-wayTandemVibrationSystemwithMultipleDegreesofFreedomBasedonMATLAB/GUIInterface

YU Xiang, ZHOU Song

(Shenyang Aerospace University, Liaoning Shenyang, 110136, China)

Inherent characteristics such as natural frequencies and mode shapes are the main parameters of vibration system in engineering vibration field. In this paper, it utilizes the MATLAB/GUI interface design platform and proposes a calculation method for one-way tandem vibration system with multiple degrees of freedom based on matrix iteration. In addition, it validates the method's rationality by solving practical examples and analyzing errors. The method lays a foundation for preliminary decoupling the response of complex vibration systems.

MATLAB; Multiple Segrees of Freedom; Vibration System; Matrix Iteration; Inherent Characteristics

10.3969/j.issn.2095-509X.2014.09.015

2013-11-07

于翔(1988—)，男，山東煙臺(tái)人，沈陽(yáng)航空航天大學(xué)碩士研究生，主要研究方向?yàn)檎駝?dòng)載荷下材料疲勞壽命預(yù)測(cè)。

O327；TB115

A

2095-509X(2014)09-0062-04