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

于 翔，周 松

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

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

于 翔，周 松

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

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

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

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

1 固有特性的求解方法

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

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

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

(1)

(2)

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

[D]=[δ][M]

(3)

通過矩陣迭代，在規(guī)定的有效位數(shù)內(nèi)，當(dāng)發(fā)現(xiàn){A}k≈{A}k-1時，運算結(jié)束。此時，{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)階固有特性時，每次迭代前都乘以清型變換后的動力矩陣[D]s。根據(jù)式(2)，則有：

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

(6)

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

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

2 程序設(shè)計

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

2.1GUI界面操作

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

圖2 軟件操作界面

2.2MATLAB程序

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

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

通過對多自由度單向串聯(lián)振動系統(tǒng)的研究發(fā)現(xiàn)，柔度矩陣主對角線上元素為各自由度上剛度倒數(shù)值的累加和，因此通過雙層循環(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; %主對角線對稱元素賦值

j=j+1;

end

i=i+1;

end

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

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

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

D=K0*M0; %初始動力矩陣

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=D-(A(:,2)*A(:,2)'*M0).*a./M;

t=t+1;

end

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

3 算例分析

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

圖4 三自由度簡化振動系統(tǒng)

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

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

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

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

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

圖5 各單元振動模擬圖像

4 結(jié)束語

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

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

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

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

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

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

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

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

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

[9] 張憲，何洋，鐘江，等. 疲勞振動試驗臺的模態(tài)與諧響應(yīng)分析[J]. 機(jī)械設(shè)計與制造，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的新型界面開發(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—)，男，山東煙臺人，沈陽航空航天大學(xué)碩士研究生，主要研究方向為振動載荷下材料疲勞壽命預(yù)測。

O327；TB115

A

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