MATLAB第章 振動

1、7.1 7.1 軌跡軌跡7.2 7.2 單自由度系統(tǒng)單自由度系統(tǒng)7.3 7.3 多自由度系統(tǒng)多自由度系統(tǒng) 7.1 7.1 軌跡軌跡舉例說明舉例說明:重力場中有兩個物體重力場中有兩個物體,其中質(zhì)量為其中質(zhì)量為m2的物體固定的物體固定,而質(zhì)而質(zhì)量為量為m1的物體繞的物體繞m2做平面圓周運動。做圓周運動的做平面圓周運動。做圓周運動的m1物體的軌物體的軌道半徑用變量道半徑用變量r表示表示,角度用變量角度用變量表示。表示。兩物體系統(tǒng)兩物體系統(tǒng) 衛(wèi)星繞地球轉(zhuǎn)動時，衛(wèi)星繞地球轉(zhuǎn)動時，m2等于地球的質(zhì)量，等于地球的質(zhì)量，m1等于衛(wèi)星的質(zhì)等于衛(wèi)星的質(zhì)量，量，r為衛(wèi)星球心與地球球心間的距離。其運動軌跡由下列方程為

2、衛(wèi)星球心與地球球心間的距離。其運動軌跡由下列方程組決定：組決定：引入狀態(tài)變量引入狀態(tài)變量:建立函數(shù)文件建立函數(shù)文件orbit.mfunction xd=orbit(t,x)xd=x(2);x(1)*x(4)2-4.0*pi2/x(1)2; x(4);-2.0*x(2)*x(4)/x(1);三組初始條件三組初始條件(t=0)：由初始條件建立執(zhí)行文件由初始條件建立執(zhí)行文件menu71.minitcond=2 0 0 1.5;1 0 0 2*pi;2 0 0 4;tspan=linspace(0,5,1000);options=odeset(RelTol,1e-6,AbsTol,1e-6 1e-6

3、1e-6 1e-6);lintype=-. -. -.;for i=1:3 t,x=ode45(orbit,tspan,initcond(i,:),options); polar(x(:,3),x(:,1),lintype(2*(i-1)+1:2*i); hold onendtext(0.5,-1.2,橢圓軌跡橢圓軌跡);text(-1.2,1,圓軌跡圓軌跡);text(1.75,2,雙曲線軌跡雙曲線軌跡);運行結(jié)果運行結(jié)果 7.2 7.2 單自由度系統(tǒng)單自由度系統(tǒng)一、概述一、概述1、力學(xué)模型、力學(xué)模型其中：振體質(zhì)量為其中：振體質(zhì)量為m，彈簧的線性系數(shù)為，彈簧的線性系數(shù)為k，非線性系數(shù)為，非線

4、性系數(shù)為，阻尼系數(shù)為阻尼系數(shù)為c，外力，外力F(t）。）。彈簧彈簧質(zhì)量質(zhì)量阻尼系統(tǒng)阻尼系統(tǒng)2、運動微分方程、運動微分方程用用x表示系統(tǒng)的位移，則運動微分方程為：表示系統(tǒng)的位移，則運動微分方程為：引入新變量轉(zhuǎn)化狀態(tài)空間方程形式：引入新變量轉(zhuǎn)化狀態(tài)空間方程形式：二、線性系統(tǒng)的自由振動二、線性系統(tǒng)的自由振動1、運動微分方程、運動微分方程當(dāng)線性項為當(dāng)線性項為0時，得到線性振動系統(tǒng)的自由振動方程。時，得到線性振動系統(tǒng)的自由振動方程。2、MATLAB求解求解對應(yīng)的函數(shù)文件對應(yīng)的函數(shù)文件FreeOcillation.mfunction xdot=FreeOcillation(t,x,dummy,zeta)x

5、dot=x(2);-2.0*zeta*x(2)-x(1);三種阻尼系數(shù)三種阻尼系數(shù)（1）阻尼系數(shù)為）阻尼系數(shù)為0.1時是欠阻尼情況；（時是欠阻尼情況；（2）阻）阻尼系數(shù)為尼系數(shù)為1時是臨界阻尼情況；（時是臨界阻尼情況；（3）阻尼系數(shù)為）阻尼系數(shù)為5時是過阻尼時是過阻尼情況。情況。由初始條件（位移和速度均為由初始條件（位移和速度均為1時）建立執(zhí)行文件時）建立執(zhí)行文件menu72.mzeta=0.1 1.0 5.0;tspan=linspace(0,40,400);lintype=-b -r -r;for i=1:3 t,x=ode45(FreeOcillation,tspan,1 1,zeta(

6、i); subplot(2,1,1); plot(t,x(:,1),lintype(2*(i-1)+1:2*i); hold on subplot(2,1,2); plot(x(:,1),x(:,2),lintype(2*(i-1)+1:2*i); hold onendsubplot(2,1,1);xlabel(Time( tau);ylabel(Displacement x( tau);title(Displacement as a function of( tau);axis(0 40 -2.0 2.0);text(2.7,-1.3,阻尼系數(shù)阻尼系數(shù)=0.1);text(3.6,-0.1,

7、1.0);text(3.6,1.0,5.0);subplot(2,1,2);xlabel(Displacement);ylabel(Velocity);title(Phase portrait);axis(-2.0 2.0 -2.0 2.0);text(0.7,-1.25,阻尼系數(shù)阻尼系數(shù)=0.1);text(0.8,-0.65,1.0);text(0.8,0.1,5.0);運行結(jié)果運行結(jié)果三、線性系統(tǒng)的強迫振動三、線性系統(tǒng)的強迫振動1、運動微分方程、運動微分方程2、MATLAB求解求解對應(yīng)的函數(shù)文件對應(yīng)的函數(shù)文件ForceOcillation.mfunction xdot=ForceOcil

8、lation(t,x,dummy,zeta,Omega,x0)xdot=x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t);為了獲得頻譜圖為了獲得頻譜圖,建立函數(shù)文件建立函數(shù)文件AmplitudeSpectrum.mfunctionf,amplitude=AmplitudeSpectrum(yy,Fs,Nstart,N);f=(Fs*(0:N-1)/N)*2.0*pi;amplitude=abs(fft(yy(Nstart:Nstart+N),N)/N;采樣速率采樣速率30/6000=0.005,則采樣頻率則采樣頻率1/0.005=200,這個頻率遠(yuǎn)遠(yuǎn)這個頻率遠(yuǎn)遠(yuǎn)

9、超出了必須達到的采樣頻率超出了必須達到的采樣頻率,結(jié)果顯示截短頻譜圖結(jié)果顯示截短頻譜圖,需設(shè)置需設(shè)置Nstart=3200,N=211=2048。 t = 0:0.001:0.6; x = sin(2*pi*50*t)+sin(2*pi*120*t); y = x + 2*randn(size(t); plot(y(1:50) title(Signal Corrupted with Zero-Mean Random Noise) xlabel(time (seconds)fft的應(yīng)用的應(yīng)用Fft變換程序變換程序Y = fft(y,512);Pyy = Y.* conj(Y) / 512;f =

10、 1000*(0:256)/512;plot(f,Pyy(1:257)title(Frequency content of y)xlabel(frequency (Hz)時域運行結(jié)果時域運行結(jié)果頻域運行結(jié)果頻域運行結(jié)果編制執(zhí)行文件編制執(zhí)行文件menu72f.mzeta=0.4;Omega=3.0;x0=50;tspan=linspace(0,30,6000);options=odeset(RelTol,1e-8,AbsTol,1e-8);lintype=-b; t,x=ode45(ForceOcillation,tspan,0 0,options,zeta,Omega,x0); subplot

11、(2,1,1); plot(t,x(:,1); axis(0 30 -8 8); hold on subplot(2,1,2); yy=x(:,1); N=2048;Nstart=3200;Fs=200;f,Amplitude=AmplitudeSpectrum(yy,Fs,Nstart,N); semilogy(f(1:40),2*Amplitude(1:40);xlabel(Frequency);ylabel(Amplitude);title(Response spectrum of a linear system); hold onsubplot(2,1,1);xlabel(Time(

12、tau);ylabel(Displacement x( tau);title(Response of a linear system);hold on運行結(jié)果運行結(jié)果頻率響應(yīng)頻率響應(yīng) 階躍響應(yīng)階躍響應(yīng) 脈沖響應(yīng)脈沖響應(yīng)Bode函數(shù)函數(shù):g = tf(1 0.1 7.5,1 0.12 9 0 0);bode(g)2342912. 05 . 71 . 0)(ssssssHbode(g,0.1 , 100)nyquist(g)impulse(g) 脈沖響應(yīng)脈沖響應(yīng)step(g) 階躍響應(yīng)階躍響應(yīng)一、建立系統(tǒng)的運動微分方程一、建立系統(tǒng)的運動微分方程1 1. .用剛度影響系數(shù)法建立振動微分方程用剛度影響

13、系數(shù)法建立振動微分方程000000000212122221112112121 nnnnnnnnnyyykkkkkkkkkyyymmm0ykyM 可簡寫為可簡寫為 7.3 7.3 自由振動模態(tài)及固有頻率自由振動模態(tài)及固有頻率例：例：三質(zhì)量三質(zhì)量m1，m2，m3串聯(lián)于彈簧串聯(lián)于彈簧k1，k2，k3上，試列出自上，試列出自由振動微分方程。由振動微分方程。 解：解：求剛度系數(shù)。求剛度系數(shù)。 給給x1以單位位移，以單位位移，x2與與x3保持保持不動，即不動，即x1=1，x2=x3=0； 要產(chǎn)生這樣的位移狀態(tài)，在各要產(chǎn)生這樣的位移狀態(tài)，在各點加的力就是點加的力就是k11，k21，k31，顯，顯然有：然有：

14、 k11=k1+k2，k21=k2，k31=0（得到第一列）（得到第一列）再令再令x2=1x2=1，x1=x3=0 x1=x3=0，則有則有 k12=k12=k2k2，k22=k2+k3k22=k2+k3，k32=k32=k3k3mmmmmm1x2x3x1k2k3k11x11k12kmmm12k22k12x32kmmm32k13x33k于是得到剛度矩陣：于是得到剛度矩陣： 33332222133323123222113121100kkkkkkkkkkkkkkkkkkk 由于系統(tǒng)有對應(yīng)于廣義坐標(biāo)由于系統(tǒng)有對應(yīng)于廣義坐標(biāo)x1，x2，x3的集中質(zhì)量，故質(zhì)的集中質(zhì)量，故質(zhì)量矩陣是對角陣，于是得到自由振

15、動微分方程：量矩陣是對角陣，于是得到自由振動微分方程： 00000000000321333322221321321xxxkkkkkkkkkxxxmmm 2.2.用柔度影響系數(shù)法建立振動微分方程用柔度影響系數(shù)法建立振動微分方程 柔度影響系數(shù)法就是力法，柔度影響系數(shù)法就是力法，它所建立方程是各點的位移協(xié)調(diào)它所建立方程是各點的位移協(xié)調(diào)方程。柔度影響系數(shù)方程。柔度影響系數(shù)cij定義：在定義：在j點作用有單位力，而其他各點沒點作用有單位力，而其他各點沒有力，所引起有力，所引起i點的位移。點的位移。 iijcjjcj點（力） i點（位移） 1 cij Pj cjj1P2PiPjPnP 從右圖看出，當(dāng)從右圖

16、看出，當(dāng)j點作用單位力點作用單位力時，在梁上各點產(chǎn)生的位移時，在梁上各點產(chǎn)生的位移c1j，c2j，cij，cnj。我們暫時不管其他。我們暫時不管其他點，只研究點，只研究i點引起位移點引起位移cijPj。 同樣，在同樣，在1點加力點加力P1，在，在i點產(chǎn)生位移點產(chǎn)生位移ci1P1； 這就是這就是i點的位移協(xié)調(diào)方程，用同樣方法可以點的位移協(xié)調(diào)方程，用同樣方法可以得到梁上各點的位移協(xié)調(diào)方程：得到梁上各點的位移協(xié)調(diào)方程： 設(shè)在力系設(shè)在力系P1，P2Pn作用下作用下i點的位移是點的位移是yi，那么必有：那么必有： 同樣，在同樣，在2點加力點加力P2，在，在i點產(chǎn)生位移點產(chǎn)生位移ci2P2； 同樣，在同樣

17、，在j點加力點加力Pj，在，在i點產(chǎn)生位移點產(chǎn)生位移cijPj； 同樣，在同樣，在n點加力點加力Pn，在，在i點產(chǎn)生位移點產(chǎn)生位移cinPn； niniiiiiPcPcPcPcy22111 寫成矩陣形式：寫成矩陣形式： nnPcPcPcy12121111nnPcPcPcy22221212nnnnnnPcPcPcy2211 Pcy（1） 其中其中 nnnnnncccccccccc212222111211 稱為柔度矩陣。由功的互等定理可知稱為柔度矩陣。由功的互等定理可知cij=cji，因此因此c是對稱矩陣。我們前邊講過，自由振動是對稱矩陣。我們前邊講過，自由振動時時 ，即，即 iiiymP 于是（

18、于是（1）式可寫成：）式可寫成： 此即用柔度影響系數(shù)法建立的多自由度系統(tǒng)自此即用柔度影響系數(shù)法建立的多自由度系統(tǒng)自由振動微分方程。由振動微分方程。 nnnnnnnnnyyymmmcccccccccyyy 212121222211121121000000或簡寫成：或簡寫成： yMcy （2） 討論討論 ：就是說，剛度矩陣和柔度矩陣是互為逆矩陣。就是說，剛度矩陣和柔度矩陣是互為逆矩陣。 1、（、（2）式可寫成）式可寫成 0yMcy 等式兩邊同乘等式兩邊同乘 ，則，則 1c 011ycyMcc 即即 01ycyM 與前邊得到的與前邊得到的 比較，可見：比較，可見： 0ykyM 1 ck按逆陣性質(zhì)：按

19、逆陣性質(zhì)： 1 kc 2、方程（、方程（2）與建立的方程形式不同，實質(zhì)是）與建立的方程形式不同，實質(zhì)是一樣的。一樣的。 解：首先計算柔度系數(shù)，解：首先計算柔度系數(shù)，用莫爾法：用莫爾法： 例：一懸臂梁，固定三個集中質(zhì)量，梁的抗彎例：一懸臂梁，固定三個集中質(zhì)量，梁的抗彎剛度剛度EI為常數(shù)，梁質(zhì)量不計，求：試用柔度系數(shù)法為常數(shù)，梁質(zhì)量不計，求：試用柔度系數(shù)法列微分方程。列微分方程。 dxEIMMcljiij0(單位：長度單位：長度/力力) 注：注：1、Mi在在i加單位力的彎矩方程，加單位力的彎矩方程，Mj在在j加單位力的彎矩方程；加單位力的彎矩方程；2、原公式中、原公式中 號號指彎矩方程不連續(xù)時，分

20、段積分后求和。指彎矩方程不連續(xù)時，分段積分后求和。 lll123mmm于是：于是： EIldxEIxldxEIMcll3)(30202111EIldxEIxldxEIMcll38)2(3202202222EIldxEIxldxEIMcll39)3(3302302333EIldxEIxlxldxEIMMccll65)2)(300212112EIldxEIxlxldxEIMMccll43)3)(300313113EIldxEIxlxldxEIMMccll314)3)(2(32020323223故柔度矩陣為故柔度矩陣為 即即 自由振動微分方程：自由振動微分方程： 542882816585263EIl

21、c321332100000054288281658526yyymmmEIlyyy 321332154288281658526yyyEImlyyy 關(guān)于固有頻率及固有振型的概念，我們已經(jīng)學(xué)關(guān)于固有頻率及固有振型的概念，我們已經(jīng)學(xué)習(xí)過。研究多自由度系統(tǒng)的這些固有特征時習(xí)過。研究多自由度系統(tǒng)的這些固有特征時, , 采用采用矩陣為數(shù)學(xué)工具。矩陣為數(shù)學(xué)工具。 二、固有頻率與固有振型二、固有頻率與固有振型 注意注意: :對于一個正定系統(tǒng)（即系統(tǒng)沒有剛體位對于一個正定系統(tǒng)（即系統(tǒng)沒有剛體位移），用上述各種方法建立微分方程都可以，有時移），用上述各種方法建立微分方程都可以，有時用柔度系統(tǒng)法更為方便；但對于半正

22、定系統(tǒng)（系統(tǒng)用柔度系統(tǒng)法更為方便；但對于半正定系統(tǒng)（系統(tǒng)有剛體位移），則柔度系統(tǒng)沒有意義（加單位力后有剛體位移），則柔度系統(tǒng)沒有意義（加單位力后系統(tǒng)不能維持平衡而產(chǎn)生剛體運動），只能用形如系統(tǒng)不能維持平衡而產(chǎn)生剛體運動），只能用形如My+ky=0那樣的方程。那樣的方程。 y 對于有對于有n個自由度個自由度q1，q2qn的系統(tǒng)，其自由的系統(tǒng)，其自由振動微分方程的一般形式是：振動微分方程的一般形式是： 或?qū)懗删仃囆问交驅(qū)懗删仃囆问?這是一組二階常系數(shù)、線性常微分方程。設(shè)系這是一組二階常系數(shù)、線性常微分方程。設(shè)系統(tǒng)偏離平衡位置作自由振動時，存在各統(tǒng)偏離平衡位置作自由振動時，存在各qi按同一頻按同一頻

23、率率，同一相位，同一相位作簡諧振動的特解：作簡諧振動的特解： 0qkqM （1） )(11tieAq)(22tieAq )(tinneAq)(tieAq 這是一組關(guān)于振幅這是一組關(guān)于振幅A1，A2An的線性齊次的線性齊次代數(shù)方程組，要使代數(shù)方程組，要使A1，A2An有非零解（零解對有非零解（零解對應(yīng)不振動），必有系數(shù)行列式為零，即：應(yīng)不振動），必有系數(shù)行列式為零，即： 代入方程組（代入方程組（1），消去公因子），消去公因子 后，得后，得 )(tie0)()()(121212212111211nnnAmkAmkAmk0)()()2(222222222111221nnnAmkAmkAmk0)()2

24、()(222212222121221nnnAmkAmkAmk（2） 或或 0)(2AkM02kM或或 此即多自由度系統(tǒng)的特征方程（或稱頻率此即多自由度系統(tǒng)的特征方程（或稱頻率方程）方程） 022221211211221112211nnnnnnnnnnmkmkmkmkmkmk 在這個方程中，在這個方程中，kij，mij都是已知的，未知的只都是已知的，未知的只有有，展開后是一個關(guān)于，展開后是一個關(guān)于2 2的的n n次代數(shù)方程，從中次代數(shù)方程，從中可以求得可以求得2 2的的n n個根。這些根就是系統(tǒng)的固有頻率。個根。這些根就是系統(tǒng)的固有頻率。由小到大排列：由小到大排列：1 1，2 2n n叫做叫做基

25、階，二基階，二階階nn階固有頻率。階固有頻率。 把每一個固有頻率例如把每一個固有頻率例如代回方程組（代回方程組（2 2），就可），就可以求得該階固有振型：以求得該階固有振型： )1()1(3)1(2)1(1:nAAAA 這里解釋一下。在二自由度系統(tǒng)時，求得了固有這里解釋一下。在二自由度系統(tǒng)時，求得了固有頻率頻率1 1，2 2之后，把之后，把1 1（或（或2 2）代回振型方程組求）代回振型方程組求該階振型時，是代回兩個方程中的任意一個，求得該階振型時，是代回兩個方程中的任意一個，求得 ，及，及 。就是說，振型方程組的兩個方程不。就是說，振型方程組的兩個方程不獨立，而是線性相關(guān)。對于多自由度系統(tǒng)，

26、也是一樣，獨立，而是線性相關(guān)。對于多自由度系統(tǒng)，也是一樣，因為振型方程組（因為振型方程組（2 2）的系統(tǒng)行列式為零，即）的系統(tǒng)行列式為零，即，就說明矩陣，就說明矩陣 的秩的秩rnrn，至少有一個方程與，至少有一個方程與其他方程線性相關(guān)。我們在求固有振型時，應(yīng)劃出一其他方程線性相關(guān)。我們在求固有振型時，應(yīng)劃出一個方程，例如最后一個，而把剩下的個方程，例如最后一個，而把剩下的n-1n-1個方程中某一個方程中某一項移到等式右邊，例如把項移到等式右邊，例如把A A1 1項移過去。那么解出的項移過去。那么解出的A A2 2，A A3 3AAn n都包含都包含A A1 1，也就是說，求出了，也就是說，求出

27、了A A2 2，A A3 3AAn n各各自對自對A A1 1的比值，那么它們之間的比值自然也就確定了。的比值，那么它們之間的比值自然也就確定了。 1)1(1)1(2AA2)2(1)2(2AA02kM)(2kM 即該階固有振型：即該階固有振型： 顯然：顯然： 代回我們設(shè)的特解代回我們設(shè)的特解 )1()1(3)1(2)1(1:nAAAA注：上標(biāo)（注：上標(biāo)（1）表示第一階振型）表示第一階振型 )()1(1)1(1tieAq)()1(2)1(2tieAq )()1()1(tinneAq或或 第第1階主振動階主振動 )()1()1(tieAq)1()1(2)1(1)1()1(2)1(1:nnAAAqq

28、q第第1 1階振型：階振型： 同樣，當(dāng)把同樣，當(dāng)把2 2，3 3n n分別代入振型方程分別代入振型方程（2 2）時，可求得其他各階固有振型。即：）時，可求得其他各階固有振型。即： 就是說，當(dāng)系統(tǒng)按第一階固有頻率作振動時，就是說，當(dāng)系統(tǒng)按第一階固有頻率作振動時，各點振幅各點振幅 之間具有確定的比值，而且之間具有確定的比值，而且系統(tǒng)各點之間在任一瞬時也具有同樣確定的比值，系統(tǒng)各點之間在任一瞬時也具有同樣確定的比值，這說明系統(tǒng)有一定的振動形態(tài)，稱為第一階固有振這說明系統(tǒng)有一定的振動形態(tài)，稱為第一階固有振型。型。 )1()1(2)1(1,nAAA)1()1(2)1(1)1()1(2)1(1:nnAAA

29、qqq第第2 2階振型：階振型： )2()2(2)2(1)2()2(2)2(1:nnAAAqqq 第第n n階振型：階振型： )()(2)(1)()(2)(1:nnnnnnnnAAAqqq于是也可寫出各階主振動。于是也可寫出各階主振動。 方程（方程（1）的通解是這）的通解是這n個主振動的疊加：個主振動的疊加： 寫成矩陣形式：寫成矩陣形式： )()(1)()2(1)()1(1)1(121ntintitieAeAeAq)()(2)()2(2)()1(2)1(221ntintitieAeAeAq )()()()2()()1()1(21ntinntintinneAeAeAq )(tieAq（3） 式中

30、共有式中共有2n個特定常數(shù)個特定常數(shù) ：及：及，注：對于每一階振型來說，各階振幅之間比值，注：對于每一階振型來說，各階振幅之間比值已定，只要知道其中一個就可求得其他振幅，所以，已定，只要知道其中一個就可求得其他振幅，所以，每一振型中只有一個待定振幅。這些常數(shù)由初始每一振型中只有一個待定振幅。這些常數(shù)由初始條件（條件（t=0時的時的 及及 ）來確定。）來確定。 ) 1 () 1 (2) 1 (1,nAAAn21,nqqq21,nqqq21, 由（由（3）式可見，具有）式可見，具有n個自由度的無阻尼系統(tǒng)個自由度的無阻尼系統(tǒng)自由振動一般由自由振動一般由n個不同頻率的主振動組成。在每個不同頻率的主振動

31、組成。在每一主振動中，系統(tǒng)各點的位移始終保持一定的比例一主振動中，系統(tǒng)各點的位移始終保持一定的比例關(guān)系（等于各點振幅之比），即系統(tǒng)有一定的振動關(guān)系（等于各點振幅之比），即系統(tǒng)有一定的振動形態(tài)，稱為主振型（或固有振型）?？梢匀稳〕銎湫螒B(tài)，稱為主振型（或固有振型）?？梢匀稳〕銎渲幸稽c的振幅（例如中一點的振幅（例如A1）等于）等于1，使其標(biāo)準(zhǔn)化，這，使其標(biāo)準(zhǔn)化，這樣的振型，稱為標(biāo)準(zhǔn)化了的固有振型。樣的振型，稱為標(biāo)準(zhǔn)化了的固有振型。 解：取三個扭轉(zhuǎn)角解：取三個扭轉(zhuǎn)角1 1，2 2，3 3為廣義坐為廣義坐標(biāo)。標(biāo)。 例一例一、懸臂軸抗扭剛度、懸臂軸抗扭剛度GJ，質(zhì)量不計，其上，質(zhì)量不計，其上固定三個圓盤，

32、轉(zhuǎn)動慣量是固定三個圓盤，轉(zhuǎn)動慣量是I。求：系統(tǒng)的固有頻。求：系統(tǒng)的固有頻率，固有振型。率，固有振型。 123代入拉氏方程代入拉氏方程 232221212121IIIT22321221)(21)(2121lGJlGJlGJU0)(iiUTdtd得系統(tǒng)自由振動微分方程：得系統(tǒng)自由振動微分方程： 0)2(211lGJI 0)2(3212lGJI 0)(313lGJI 或?qū)懗删仃囆问剑夯驅(qū)懗删仃囆问剑?把得到的把得到的M、k矩陣代入振型方程矩陣代入振型方程 000110121012100010001321321lGJI 0)(2AkM注：或設(shè)注：或設(shè) 2211,AA000)11012101210001

33、0001(3212AAAlGJI設(shè)設(shè) ，則上式成為：，則上式成為：2GJIl注：可以直接代入頻率方程注：可以直接代入頻率方程 02kM頻率方程。頻率方程。 由一元三次方程求根公式解得：由一元三次方程求根公式解得： 0210121012展開為：展開為： 0165231=0.198，2=1.555，3=3.247代入代入 式中，求得三個固有頻率式中，求得三個固有頻率 2GJIlIlGJ445. 01IlGJ247. 12IlGJ802. 13 代回振型方程中的前二個（用代回振型方程中的前二個（用），把），把A3移到移到方程右端：方程右端： 將上面求得的將上面求得的1，2，3分別代入分別代入，就得，

34、就得到三個固有振型：到三個固有振型： 取取A1=1（即對標(biāo)準(zhǔn)化），得（即對標(biāo)準(zhǔn)化），得 0)2(21AA321)2(AAA11A 22A1)2(23A即：即： )34( : )2( :1:2321AAA247. 2:802. 1:1:)1(3)1(2)1(1AAA)802. 0( :445. 0:1:)1(3)1(2)1(1AAA555. 0: )247. 1( :1:)1(3)1(2)1(1AAA1802. 1247. 21445. 0802. 01247. 1555. 0編制編制m文件文件.k=2 -1 0;-1 2 -1;0 -1 1;m=1 0 0;0 1 0;0 0 1;vibrat

35、ionmode,eigenvalues=eig(k,m)v1=vibrationmode(:,1)v2=vibrationmode(:,2)v3=vibrationmode(:,3)b1=v1./v3b2=v2./v3b3=v3./v3for i=1:3subplot(3,1,i) plot(0 1 2 3,0 b3(i) b2(i) b1(i)hold onend運行結(jié)果運行結(jié)果應(yīng)用應(yīng)用MATLAB求解求解建立系統(tǒng)的運動微分方程建立系統(tǒng)的運動微分方程 在前一節(jié)看到，為了求多自由度系統(tǒng)的固有頻在前一節(jié)看到，為了求多自由度系統(tǒng)的固有頻率和固有振型，必須解高次代數(shù)方程（頻率方程或率和固有振型，必須

36、解高次代數(shù)方程（頻率方程或稱特征方程），其次數(shù)與系統(tǒng)自由度數(shù)目一致。當(dāng)稱特征方程），其次數(shù)與系統(tǒng)自由度數(shù)目一致。當(dāng)自由度數(shù)增加時，解次數(shù)很高的代數(shù)方程，是很困自由度數(shù)增加時，解次數(shù)很高的代數(shù)方程，是很困難的。所以一般用近似的數(shù)值方法。這里介紹的是難的。所以一般用近似的數(shù)值方法。這里介紹的是一種數(shù)值方法一種數(shù)值方法“矩陣迭代法矩陣迭代法”。 多自由度相同自由振動微分方程一般形式：多自由度相同自由振動微分方程一般形式： 0qkqM 設(shè)特解為：設(shè)特解為： tieAq三、用矩陣迭代法求固有頻率和固有振型三、用矩陣迭代法求固有頻率和固有振型代入方程得振型方程代入方程得振型方程 用用k-1前乘各項，移項后

37、得到前乘各項，移項后得到 02AkAM121AAMk引入符號引入符號 1MkD動力矩陣，動力矩陣， 就得到迭代格式（振型方程）：就得到迭代格式（振型方程）： 12AAD 對于一個給定的系統(tǒng)，動力矩陣對于一個給定的系統(tǒng)，動力矩陣D是可知的。是可知的。矩陣迭代法就是用逐次迭代、逐次逼近的方法求解矩陣迭代法就是用逐次迭代、逐次逼近的方法求解振型方程，以確定特征值振型方程，以確定特征值 和特征向量（固有振型）和特征向量（固有振型）A。 21集中質(zhì)量法概念m/4m/4m/4m/3m/3m/2l/4l/4l/4l/4l/3l/2l/3l/3l/2例例:集中質(zhì)量法求簡支梁的模態(tài)集中質(zhì)量法求簡支梁的模態(tài)力學(xué)模

38、型:有一簡支梁，其上有一個附加質(zhì)量（位置和尺寸均在圖上標(biāo)注），梁長L=2m ,截面尺寸b*h為0.04m*0.02m ,截面慣性矩為2.67* ，密度為7920kg/，質(zhì)量塊的質(zhì)量M為3.0kg 彈性模量E=210GPa 離散810一一.力學(xué)模型及集中質(zhì)量法簡化力學(xué)模型及集中質(zhì)量法簡化質(zhì)量矩陣0043.168000006.16804003.168004mmMmM3591171.7552.1451.36511 16112.1453.122.1451076871191.3652.1451.755lEIF柔度矩陣矩陣 = 稱作系統(tǒng)的動力矩陣 = DFM二、二、MATLAB計算計算m=1 0 0;0

39、1.94 0;0 0 1;m=1 0 0;0 1.94 0;0 0 1;h=0.00000195h=0.00000195* *3.186;3.186;b=9 11 7;11 16 11;7 11 9;b=9 11 7;11 16 11;7 11 9;f=h.f=h.* *b;b;a=fa=f* *m;m;d=1 1 1;d=1 1 1;for i=1:16;for i=1:16;y=ay=a* *d;d;d=(1/y(3,1)d=(1/y(3,1)* *y;y;endendv=dv=ds=ms=m* *y;y;w1=1/sqrt(s(3,1)subplot(3,1,1)plot(0 2,0 0

40、)hold onplot(0 0.5 1 1.5 2,0 v(1,1) v(2,1) v(3,1) 0)p=v*m*v;q=v*v*m;j=a-(s(3,1)/p)*q;d=1 1 -1;% x(1)=d;for i=1:16;y=j*d;% l=y(i)d=(1/y(3,1)*y;endv=ds=m*y;w2=1/sqrt(s(3,1)subplot(3,1,2)plot(0 2,0 0)hold onplot(0 0.5 1 1.5 2,0 v(1,1) v(2,1) v(3,1) 0)q=v*v*m;r=j-(s(3,1)/p)*q;d=1 -1 1;% x(1)=d;for i=1:1

41、6;y=r*d;% l=y(i)d=(1/y(3,1)*y;endv=ds=m*y;w3=1/sqrt(s(3,1)subplot(3,1,3)plot(0 2,0 0)hold onplot(0 0.5 1 1.5 2,0 v(1,1) v(2,1) v(3,1) 0)運行結(jié)果四、四、 強迫振動及減振器強迫振動及減振器m=50 10;k=200 40;c=10 6;N=m(2) c(2) k(2);c(2) k(2);D=m(1)*m(2) (c(1)+c(2)*m(2)+c(2)*m(1) . (k(1)+k(2)*m(2)+k(2)*m(1)+c(1)*c(2) . k(1)*c(2)+

42、c(1)*k(2) k(1)*k(2);sy=tf(N,D);% y,t=impulse(transferab(m,k,c),20)y,t=impulse(sy,20);subplot(2,1,1);plot(t,y(:,1,1);ylabel(x_1(t);title(impulse response of m_1)subplot(2,1,2);plot(t,y(:,2,1);xlabel(time t);ylabel(x_2(t)title(impulse response of m_2);編制編制MATLAB文件文件運運行行結(jié)結(jié)果果減振器：求減振器：求m1上施加一上施加一外力時，物體外力

43、時，物體m1和和m2位移的頻率響應(yīng)函數(shù)。位移的頻率響應(yīng)函數(shù)。程序如上程序如上k=200 40;c=10 6;omega=0.0:0.005:4for i=1:2 if i=1 sys=tf(1,m(1) c(1) k(1); mag,phas=bode(sys,omega); plot(omega,mag(1,:),-); hold on; elseN=m(2) c(2) k(2);c(2) k(2);D=m(1)*m(2) (c(1)+c(2)*m(2)+c(2)*m(1) . (k(1)+k(2)*m(2)+k(2)*m(1)+c(1)*c(2) . k(1)*c(2)+c(1)*k(2) k(1)*k(2);sys=tf(N,D); sys=tf(N,D); mag,phas=bode(sys,omega); plot