第6章 Matlab應(yīng)用之動力學(xué)與振動課件

第6章動力學(xué)與振動第6章Matlab應(yīng)用之動力學(xué)與振動教學(xué)目標(biāo)介紹Matlab在動力學(xué)與振動中的應(yīng)用，分別用于軌跡，單自由度和多自由度線性與非線性系統(tǒng)的自由振動和強(qiáng)迫振動的分析。學(xué)習(xí)要求

能夠運用Matlab基本原理，對物體的運動軌跡和單自由度系統(tǒng)進(jìn)行簡單的動力學(xué)分析。第6章Matlab應(yīng)用之動力學(xué)與振動目錄6.1軌跡6.2單自由度系統(tǒng)6.3多自由度系統(tǒng)習(xí)題第6章Matlab應(yīng)用之動力學(xué)與振動6.1軌跡舉例說明:重力場中有兩個物體,其中質(zhì)量為m2的物體固定,而質(zhì)量為m1的物體繞m2做平面圓周運動.做圓周運動的m1物體的軌道半徑用變量r表示,角度用變量a表示.m2m1ar第6章Matlab應(yīng)用之動力學(xué)與振動6.1軌跡例6.1：衛(wèi)星繞地球轉(zhuǎn)動時，m2等于地球的質(zhì)量，m1等于衛(wèi)星的質(zhì)量，r為衛(wèi)星球心與地球球心間的距離。其運動軌跡由下列方程組決定：式中：，其中t是時間變量，p為物體在地球表面做圓周運動的周期。在地球表面，r=6.373x106m。第6章Matlab應(yīng)用之動力學(xué)與振動6.1軌跡用龍格—庫塔法可以實現(xiàn)求解：引入新狀態(tài)變量：第6章Matlab應(yīng)用之動力學(xué)與振動建立函數(shù)文件Orbit.mfunctionxd=Orbit(t,x)xd=[x(2)x(1)*x(4)^2-4.0*pi^2/x(1)^2x(4)-2.0*x(2)*x(4)/x(1)];6.1軌跡組X1初始X2初始X3初始X4初始軌跡類型12001.5橢圓21002pi圓32004雙曲線三組初始條件(t=0)：第6章Matlab應(yīng)用之動力學(xué)與振動由初始條件建立執(zhí)行文件execute_61.minitcond=[2001.5;1002*pi;2004];tspan=linspace(0,5,1000);options=odeset('RelTol',1e-6,'AbsTol',[1e-61e-61e-61e-6]);lintype=[‘k-'‘b-.'‘r--'];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,'雙曲線軌跡');6.1軌跡常微分方程的數(shù)值求解函數(shù)第6章Matlab應(yīng)用之動力學(xué)與振動程序運行結(jié)果6.1軌跡第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)6.2.1概述一.力學(xué)模型mcK，aX（t）F（t）=X（0）kf（t）彈簧—質(zhì)量—阻尼系統(tǒng)其中：振體質(zhì)量為m，彈簧的線性系數(shù)為k，非線性系數(shù)為a，阻尼系數(shù)為c，外力F（t）。第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)二.運動微分方程用x表示系統(tǒng)的位移，則運動微分方程為：式中：固有頻率:非線性系數(shù):阻尼因子:第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)引入新變量轉(zhuǎn)化狀態(tài)空間方程形式：第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)6.2.2線性系統(tǒng)的自由振動一.運動微分方程當(dāng)時，得到線性振動系統(tǒng)的自由振動方程。第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)二.MATLAB求解編寫方程對應(yīng)的函數(shù)文件FreeOscillation.m0三種阻尼系數(shù)()（1）阻尼系數(shù)為0.1時是欠阻尼情況（2）阻尼系數(shù)為1時是臨界阻尼情況（3）阻尼系數(shù)為5時是過阻尼情況functionxdot=FreeOscillation(t,x,zeta,Alpha)xdot=[x(2);-2.0*zeta*x(2)-x(1)-Alpha*x(1)^3];end第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)由初始條件（位移和速度均為1時,）建立執(zhí)行文件(execute_62.m)zeta=[0.11.05.0];Alpha=[0.0,0.0,0.0];tspan=linspace(0,40,400);%生成0-40的四百個線性點lintype=char('-k','--k','-.k');fori=1:3[t,x]=ode45(@FreeOscillation,tspan,[11],[],zeta(i),Alpha(i));figure(1);plot(t,x(:,1),lintype(i,:));%x(:,1)為位移holdonfigure(2);plot(x(:,1),x(:,2),lintype(i,:));%x(:,2)為速度holdonend第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)figure(1);xlabel('Time(\tau)');ylabel('Displacementx(\tau)');title('Displacementasafunctionof(\tau)');axis([040-1.51.5]);plot([0,40],[0,0],'k-')legend('\zeta=0.1','\zeta=1.0','\zeta=5.0')figure(2);xlabel('Displacementx(\tau)');ylabel('Velocity');title('Phaseportrait');axis([-2.02.0-2.02.0]);legend('\zeta=0.1','\zeta=1.0','\zeta=5.0');續(xù)上：第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)程序運行結(jié)果第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)6.2.3非線性系統(tǒng)的自由振動1、運動微分方程一.非線性彈簧系統(tǒng)第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)2、Matlab求解編寫常微分方程對應(yīng)的函數(shù)文件FreeOscillation.mfunctionxdot=FreeOscillation(t,x,zeta,Alpha)xdot=[x(2);-2.0*zeta*x(2)-x(1)-Alpha*x(1)^3];end與例6.2相同，只是改變了Alpha的值，可以直接借用例6.2的函數(shù)文件第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)由初始條件建立執(zhí)行文件(execute_63.m)程序如下zeta=0.2;Alpha=[0.00,-0.25,-0.25];x0=[-2.00,-2.00,-2.00];v0=[2.00,2.00,2.31];tspan=linspace(0.0,30.0,401);lintyp=char('-k','--k','-.k');options=odeset('RelTol',1e-8,'AbsTol',[1e-81e-8]);d=char('Linear:x_0=-2v_0=2\alpha=0',...'Nonlinear:x_0=-2v_0=2\alpha=-0.25',...'Nonlinear:x_0=-2v_0=2.31\alpha=-0.25');第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)fori=1:3[t,x]=ode45(@FreeOscillation,tspan,[x0(i)v0(i)]',options,zeta,Alpha(i));figure(1)plot(t,x(:,1),lintyp(i,:));holdonfigure(2)plot(x(:,1),x(:,2),lintyp(i,:));holdonend續(xù)上：第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)figure(1)xlabel('\tau');ylabel('x(\tau)');axis([0.0,30.0,-3.0,3.0]);legend(d(1,:),d(2,:),d(3,:));figure(2)xlabel('x(\tau)');ylabel('dx/d\tau');axis([-2.0,3.0,-2.0,3.0]);legend(d(1,:),d(2,:),d(3,:));續(xù)上：第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)程序運行結(jié)果第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)二、非線性阻尼系統(tǒng)1、運動微分方程式中，常量為摩擦系數(shù)，為物體的重量，k為線性彈簧的系數(shù)。干摩擦力是速度的分段函數(shù)，用signum表示。速度為正時，signum取+1，速度為負(fù)時，signum取-1.如果彈簧的彈性力不能克服干摩擦力，系統(tǒng)將停止振動。即當(dāng)?shù)?章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)引入新變量將方程轉(zhuǎn)化一階方程形式：兩邊同時求導(dǎo)兩邊同時求導(dǎo)第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)functionxdot=FrictionOscillation(t,x,d)%非線性阻尼系統(tǒng)ode文件ifabs(x(1))<=d&&x(2)==0.0;xdot=[0;0];elsexdot=[x(2);-d*sign(x(2))-x(1)];end2、Matlab求解編寫常微分方程對應(yīng)的函數(shù)文件FrictionOscillation.m第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)由初始條件(d=0.86,初始條件a（3.0，0.0），b（5.0，0.0）)建立執(zhí)行文件(execute_64.m)，求數(shù)值解d=0.86;x0=[3.0,5.0];v0=[0.0,0.0];tspan=linspace(0,12,120);options=odeset('AbsTol',[1e-3,1e-3]);lintyp=char('--k','-k');fori=1:2;[t,x]=ode45(@FrictionOscillation,tspan,[x0(i),v0(i)]',options,d);figure(1);plot(t,x(:,1),lintyp(i,:));holdonfigure(2)plot(x(:,1),x(:,2),lintyp(i,:));holdonend第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)figure(1)xlabel('\tau');ylabel('x(\tau)');axis([0.0,12.0,-4.0,6.0]);plot([0,12],[0,0],'k-');legend('x_0=3.0,v_0=0.0','x_0=5.0,v_0=0.0');figure(2)xlabel('x(\tau)');ylabel('dx/d\tau');text(2.5,0.5,'(3.0,0.0)');text(4.5,0.5,'(5.0,0.0)');plot([-4,6],[0,0],'k-',[0,0],[-6,4],'k-');axis([-4.0,6.0,-6.0,4.0]);續(xù)上：第6章Matlab應(yīng)用之動力學(xué)與振動6.2單自由度系統(tǒng)程序運行結(jié)果第6章Matlab應(yīng)用之動力學(xué)與振動6.3多自由度系統(tǒng)6.3.1多自由系統(tǒng)的固有頻率問題一、力學(xué)模型二、運動微分方程第6章Matlab應(yīng)用之動力學(xué)與振動三、Matlab求解例6.5三自由系統(tǒng)的振動模態(tài)及固有頻率設(shè)k1=100N/m，k2=50N/m，m1=m2=m3=100kg。求特征值與特征向量的程序如下：k=[100,-100,0;-100,150,-50;0,-50,50]m=diag([100,100,100])[VibrationMode,EigenValue]=eig(k,m)第6章Matlab應(yīng)用之動力學(xué)與振動附錄：ode45函數(shù)如果系數(shù)矩陣A的特征值連乘積小于零，且絕對值最大和最小的特征值之比（剛性比）很大，則稱此類方程為剛性方程ode是Matlab專門用于解微分方程的功能函數(shù)。該求解器有變步長（variable-step）和定步長（fixed-step）兩種類型。不同類型有著不同的求解器，其中ode45求解器屬于變步長的一種，采用Runge-Kutta算法；ode45表示采用四階，五階Runge-Kutta單步算法,截斷誤差為(Δx)^3。解決的是Nonstiff(非剛性)常微分方程。第6章Matlab應(yīng)用之動力學(xué)與振動附錄：ode45函數(shù)[T,Y]=ode45(@fun,TSPAN,Y0)[T,Y]=ode45(@fun,TSPAN,Y0,options)[T,Y]=ode45(@fun,TSPAN,Y0,options,P1,P2,…)[T,Y,TE,YE,IE]=ode45(@fun,TSPAN,Y0,options,P1,P2,…)調(diào)用格式：說明：輸出變量T為返回時間列向量；解矩陣Y的每一行對應(yīng)于T的一個元素，列數(shù)與求解變量數(shù)相等。@fun為函數(shù)句柄，為根據(jù)待求解的ODE方程所編寫的ode文件（odefile）；TSPAN＝[T0TFINAL]是微分系統(tǒng)y'＝F(t,y)的積分區(qū)間；Y0為初始條件options用于設(shè)置一些可選的參數(shù)值，缺省時，相對于第一種調(diào)用格式。P1，P2，…的作用是傳遞附加參數(shù)P1，P2，…到ode文件。當(dāng)options缺省時，應(yīng)在相應(yīng)位置保留[]，以便正確傳遞參數(shù)。第6章Matlab應(yīng)用之動力學(xué)與振動附錄：ode45函數(shù)所謂的odefile實際上是一個Matlab函數(shù)文件，一般作為整個求解程序的一個子函數(shù)，表示ode求解問題ode文件的最簡單格式必須有一個自變量t和函數(shù)y作為輸入變量，一個y的導(dǎo)函數(shù)作為輸出變量。其中自變量t不論在ode文件中是否使用都必須作為第一輸入變量，y則必須作為第二輸入變量，位置不能顛倒?？梢韵騩de文件中傳遞參數(shù)，數(shù)目不受限制odefile第6章Matlab應(yīng)用之動力學(xué)與振動附錄：ode45函數(shù)為了能夠解出方程，要用指令odeset確定求解的條件和要求。在MATLAB中，求解方程組的指令都有默認(rèn)的求解的條件和要求（由結(jié)構(gòu)數(shù)組options表示），但可以用odeset修改或重新建立，odesetoptions=odeset(‘name1’,value1,’name2’,value2,