MATLAB與工程應(yīng)用-第7章-動力學(xué)與振動課件

第七章動力學(xué)與振動7.1軌跡7.2單自由度系統(tǒng)7.3多自由度系統(tǒng)第七章動力學(xué)與振動7.1軌跡17.1軌跡舉例說明:重力場中有兩個物體,其中質(zhì)量為m2的物體固定,而質(zhì)量為m1的物體繞m2做平面圓周運(yùn)動.做圓周運(yùn)動的m1物體的軌道半徑用變量r表示,角度用變量a表示.m1ram2

兩物體系統(tǒng)7.1軌跡2衛(wèi)星繞地球轉(zhuǎn)動時(shí)，m2等于地球的質(zhì)量，m1等于衛(wèi)星的質(zhì)量，r為衛(wèi)星球心與地球球心間的距離。其運(yùn)動軌跡由下列方程組決定：式中：，其中t是時(shí)間變量，p為物體在地球表面做圓周運(yùn)動的周期。在地球表面，r=6.373x106m。衛(wèi)星繞地球轉(zhuǎn)動時(shí)，m2等于地球的質(zhì)量，m1等于衛(wèi)星的質(zhì)量，r3用龍格—庫塔法可以實(shí)現(xiàn)求解：引入新狀態(tài)變量：帶入前面的微分方程組，可得四個一階微分方程。用龍格—庫塔法可以實(shí)現(xiàn)求解：4建立函數(shù)文件orbit.mfunctionxd=orbit(t,x)xd=[x(2);x(1)*x(4)^2-4.0*pi^2/x(1)^2;x(4);-2.0*x(2)*x(4)/x(1)];三組初始條件(t=0)：建立函數(shù)文件orbit.m5由初始條件建立執(zhí)行文件menu71.minitcond=[2001.5;1002*pi;2004];tspan=linspace(0,5,1000);options=odeset('RelTol',1e-6,'AbsTol',[1e-61e-61e-61e-6]);lintype=['-.''-.''-.'];fori=1:3[t,x]=ode45('orbit',tspan,[initcond(i,:)],options);polar(x(:,3),x(:,1),lintype(2*(i-1)+1:2*i));holdonendtext(0.5,-1.2,'橢圓軌跡');text(-1.2,1,'圓軌跡');text(1.75,2,'雙曲線軌跡');由初始條件建立執(zhí)行文件menu71.m6程序運(yùn)行結(jié)果程序運(yùn)行結(jié)果77.2單自由度系統(tǒng)7.2.1概述一.力學(xué)模型彈簧—質(zhì)量—阻尼系統(tǒng)其中：振體質(zhì)量為m，彈簧的線性系數(shù)為k，非線性系數(shù)為a，阻尼系數(shù)為c，外力F（t）。mcK，aX（t）F（t）=X（0）kf（t）7.2單自由度系統(tǒng)mcK，aX（t）F（t）=X（0）kf（8二.運(yùn)動微分方程用x表示系統(tǒng)的位移，則運(yùn)動微分方程為：式中：固有頻率非線性系數(shù)

阻尼因子引入新變量轉(zhuǎn)化狀態(tài)空間方程形式：二.運(yùn)動微分方程97.2.2線性系統(tǒng)的自由振動一.運(yùn)動微分方程當(dāng)時(shí)，得到線性振動系統(tǒng)的自由振動方程。二.MATLAB求解對應(yīng)的函數(shù)文件FreeOcillation.mfunctionxdot=FreeOcillation(t,x,dummy,zeta)xdot=[x(2);-2.0*zeta*x(2)-x(1)];三種阻尼系數(shù)（1）阻尼系數(shù)為0.1時(shí)是欠阻尼情況（2）阻尼系數(shù)為1時(shí)是臨界阻尼情況（3）阻尼系數(shù)為5時(shí)是過阻尼情況7.2.2線性系統(tǒng)的自由振動10由初始條件（位移和速度均為1時(shí)）建立執(zhí)行文件menu72.mzeta=[0.11.05.0];tspan=linspace(0,40,400);%生成0-40的四百個線性點(diǎn)lintype=['-b''--r''---r'];fori=1:3

[t,x]=ode45('FreeOcillation',tspan,[11],[],zeta(i));subplot(2,1,1);plot(t,x(:,1),lintype(2*(i-1)+1:2*i));holdonsubplot(2,1,2);plot(x(:,1),x(:,2),lintype(2*(i-1)+1:2*i));holdonendsubplot(2,1,1);由初始條件（位移和速度均為1時(shí)）建立執(zhí)行文件menu72.m11xlabel('Time(\tau)');ylabel('Displacementx(\tau)');title('Displacementasafunctionof(\tau)');axis([040-2.02.0]);text(2.7,-1.3,'阻尼系數(shù)=0.1');text(3.6,-0.1,'1.0');text(3.6,1.0,'5.0');subplot(2,1,2);xlabel('Displacement');ylabel('Velocity');title('Phaseportrait');axis([-2.02.0-2.02.0]);text(0.7,-1.25,'阻尼系數(shù)=0.1');text(0.8,-0.65,'1.0');text(0.8,0.1,'5.0');xlabel('Time(\tau)');12程序運(yùn)行結(jié)果程序運(yùn)行結(jié)果137.2.3線性系統(tǒng)的強(qiáng)迫振動一.運(yùn)動微分方程二.MATLAB求解若對應(yīng)的函數(shù)文件ForceOcillation.mfunctionxdot=ForceOcillation(t,x,dummy,zeta,Omega,x0)xdot=[x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t)];

7.2.3線性系統(tǒng)的強(qiáng)迫振動14為了獲得頻譜圖,建立函數(shù)文件AmplitudeSpectrum.mfunction[f,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)超出了必須達(dá)到的采樣頻率,結(jié)果顯示截短頻譜圖,需設(shè)置Nstart=3200,N=2^11=2048。fft的應(yīng)用見Help編制執(zhí)行文件menu72f.m為了獲得頻譜圖,建立函數(shù)文件AmplitudeSpectru15zeta=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,[00],options,zeta,Omega,x0);subplot(2,1,1);plot(t,x(:,1));axis([030-88]);holdonsubplot(2,1,2);functionxdot=ForceOcillation(t,x,dummy,zeta,Omega,x0)xdot=[x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t)];zeta=0.4;Omega=3.0;x0=50;funct16yy=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('Responsespectrumofalinearsystem');holdonsubplot(2,1,1);xlabel('Time(\tau)');ylabel('Displacementx(\tau)');title('Responseofalinearsystem');holdonyy=x(:,1);17程序運(yùn)行結(jié)果程序運(yùn)行結(jié)果187.2.4線性系統(tǒng)的頻率響應(yīng)、階躍響應(yīng)及脈沖響應(yīng)單自由度振動系統(tǒng)的強(qiáng)迫振動微分方程可為：通過LAPLACE變換，得傳遞函數(shù)：其中：Logspace,rad2deg,loglog,semilogxseehelp7.2.4線性系統(tǒng)的頻率響應(yīng)、階躍響應(yīng)及脈沖響應(yīng)Logs19一.編制執(zhí)行文件frequency72.m，求頻率響應(yīng)。m=1;zeta=0.1:0.1:1;k=1;wn=sqrt(k/m);10w=logspace(-1,1,400);rad2deg=180/pi;s=j*w;forcnt=1:length(zeta)xfer(cnt,:)=(1/m)./(s.^2+2*zeta(cnt)*wn*s+wn^2);mag(cnt,:)=abs(xfer(cnt,:));phs(cnt,:)=angle(xfer(cnt,:))*rad2deg;endforcnt=1:length(zeta)figure(1)loglog(w,mag(cnt,:),'k-')一.編制執(zhí)行文件frequency72.m，求頻率響應(yīng)。20title('SDOFfrequencyresponsemagnitudesforzeta=0.2to1.0instepsof0.2')xlabel('Frequency(rad/sec)')ylabel('Magnitude')gridholdonendholdoffforcnt=1:length(zeta)figure(2)semilogx(w,phs(cnt,:),'k-')title('SDOFfrequencyresponsephasesforzeta=0.2to1.0instepsof0.2')xlabel('Frequency(rad/sec)')ylabel('phase')gridholdonendholdofftitle('SDOFfrequencyresponse21程序運(yùn)行結(jié)果幅頻曲線程序運(yùn)行結(jié)果幅頻曲線22程序運(yùn)行結(jié)果相頻曲線程序運(yùn)行結(jié)果相頻曲線23二.求時(shí)頻響應(yīng)的基本函數(shù)命令可以通過上述命令求線性系統(tǒng)的波得圖、乃奎斯特圖、階躍響應(yīng)、脈沖響應(yīng)、初始條件響應(yīng)、輸入u的響應(yīng)二.求時(shí)頻響應(yīng)的基本函數(shù)命令可以通過上述命令求線性系統(tǒng)的波得24例：博得圖編制執(zhí)行文件bode72.mm=1zeta=0.1k=1wn=sqrt(k/m)den=[12*zeta*wnwn^2]num=[1/m]sys=tf(num,den)bode(sys)Bode,nyquistseehelp例：博得圖Bode,nyquistseehelp25例：求正弦輸入激勵響應(yīng)編制執(zhí)行文件sin72.mt=0:0.01:50;u=sin(t);lsim(sys,u,t)Lsim