應(yīng)用之動(dòng)力學(xué)和振動(dòng)

1、關(guān)于應(yīng)用之動(dòng)力學(xué)與振動(dòng)第一張，PPT共三十七頁，創(chuàng)作于2022年6月教學(xué)目標(biāo)介紹Matlab在動(dòng)力學(xué)與振動(dòng)中的應(yīng)用，分別用于軌跡，單自由度和多自由度線性與非線性系統(tǒng)的自由振動(dòng)和強(qiáng)迫振動(dòng)的分析。學(xué)習(xí)要求 能夠運(yùn)用Matlab基本原理，對物體的運(yùn)動(dòng)軌跡和單自由度系統(tǒng)進(jìn)行簡單的動(dòng)力學(xué)分析。第二張，PPT共三十七頁，創(chuàng)作于2022年6月目錄6.1 軌跡6.2 單自由度系統(tǒng)6.3 多自由度系統(tǒng)習(xí)題第三張，PPT共三十七頁，創(chuàng)作于2022年6月6.1 軌跡舉例說明:重力場中有兩個(gè)物體,其中質(zhì)量為m2的物體固定,而質(zhì)量為m1的物體繞m2做平面圓周運(yùn)動(dòng).做圓周運(yùn)動(dòng)的m1物體的軌道半徑用變量r表示,角度用變量a

2、表示.m2m1ar第四張，PPT共三十七頁，創(chuàng)作于2022年6月6.1 軌跡例6.1：衛(wèi)星繞地球轉(zhuǎn)動(dòng)時(shí)，m2等于地球的質(zhì)量，m1等于衛(wèi)星的質(zhì)量，r為衛(wèi)星球心與地球球心間的距離。其運(yùn)動(dòng)軌跡由下列方程組決定：式中： ，其中t是時(shí)間變量，p為物體在地球表面做圓周運(yùn)動(dòng)的周期。在地球表面，r=6.373x106 m。第五張，PPT共三十七頁，創(chuàng)作于2022年6月6.1 軌跡用龍格庫塔法可以實(shí)現(xiàn)求解：引入新狀態(tài)變量：第六張，PPT共三十七頁，創(chuàng)作于2022年6月建立函數(shù)文件Orbit.mfunction xd=Orbit(t,x)xd=x(2)x(1)*x(4)2-4.0*pi2/x(1)2 x(4)-2

3、.0*x(2)*x(4)/x(1);6.1 軌跡組X1初始X2初始X3初始X4初始軌跡類型12001.5橢圓21002pi圓32004雙曲線三組初始條件(t=0)：第七張，PPT共三十七頁，創(chuàng)作于2022年6月由初始條件建立執(zhí)行文件execute_61.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 1e-6 1e-6);lintype=k- b-. r-;for i=1:3 t,x=ode45(Orbit,tspan,init

4、cond(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,雙曲線軌跡);6.1 軌跡常微分方程的數(shù)值求解函數(shù)第八張，PPT共三十七頁，創(chuàng)作于2022年6月程序運(yùn)行結(jié)果6.1 軌跡第九張，PPT共三十七頁，創(chuàng)作于2022年6月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（

5、t）。第十張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)二. 運(yùn)動(dòng)微分方程用x表示系統(tǒng)的位移，則運(yùn)動(dòng)微分方程為：式中：固有頻率:非線性系數(shù):阻尼因子:第十一張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)引入新變量轉(zhuǎn)化狀態(tài)空間方程形式：第十二張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)6.2.2 線性系統(tǒng)的自由振動(dòng)一.運(yùn)動(dòng)微分方程當(dāng) 時(shí)，得到線性振動(dòng)系統(tǒng)的自由振動(dòng)方程。第十三張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)二.MATLAB求解編寫方程對應(yīng)的函數(shù)文件FreeOscillation.m0三種阻尼系數(shù)( )（1）阻尼系數(shù)為

6、0.1時(shí)是欠阻尼情況（2）阻尼系數(shù)為1時(shí)是臨界阻尼情況（3）阻尼系數(shù)為5時(shí)是過阻尼情況function xdot=FreeOscillation(t,x,zeta,Alpha)xdot=x(2);-2.0*zeta*x(2)-x(1)-Alpha*x(1)3; end第十四張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)由初始條件（位移和速度均為1時(shí), ）建立執(zhí)行文件(execute_62.m)zeta=0.1 1.0 5.0; Alpha=0.0,0.0,0.0;tspan=linspace(0,40,400); %生成0-40的四百個(gè)線性點(diǎn)lintype=char(-k,-k

7、,-.k);for i=1:3 t,x=ode45(FreeOscillation,tspan,1 1,zeta(i),Alpha(i); figure(1); plot(t,x(:,1), lintype(i,:); % x(:,1)為位移 hold on figure(2); plot(x(:,1),x(:,2),lintype(i,:); % x(:,2)為速度 hold onend 第十五張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)figure(1);xlabel(Time( tau);ylabel(Displacement x( tau);title(Displac

8、ement as a function of( tau);axis(0 40 -1.5 1.5);plot(0,40,0,0, k-)legend(zeta=0.1, zeta=1.0, zeta=5.0)figure(2);xlabel(Displacement x(tau);ylabel(Velocity);title(Phase portrait);axis(-2.0 2.0 -2.0 2.0);legend(zeta=0.1, zeta=1.0, zeta=5.0);續(xù)上：第十六張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)程序運(yùn)行結(jié)果第十七張，PPT共三十七頁，創(chuàng)作于

9、2022年6月6.2 單自由度系統(tǒng)6.2.3 非線性系統(tǒng)的自由振動(dòng)1、運(yùn)動(dòng)微分方程一.非線性彈簧系統(tǒng)第十八張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)2、Matlab求解編寫常微分方程對應(yīng)的函數(shù)文件FreeOscillation.mfunction xdot=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ù)文件第十九張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)由初始條件建立執(zhí)行文件(ex

10、ecute_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-8 1e-8);d=char(Linear:x_0=-2 v_0=2 alpha=0,. Nonlinear:x_0=-2 v_0=2 alpha=-0.25,. Nonlinear:x_0=-2 v_0=2.31 alpha=-0.25);第二十張

11、，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)for i=1:3 t,x=ode45(FreeOscillation,tspan,x0(i) v0(i),options,zeta,Alpha(i); figure(1) plot(t,x(:,1),lintyp(i,:); hold on figure(2) plot(x(:,1),x(:,2),lintyp(i,:); hold onend續(xù)上：第二十一張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)figure(1)xlabel(tau);ylabel(x(tau);axis(0.0,30.0,-3.0,3.0)

12、;legend(d(1,:),d(2,:),d(3,:);figure(2)xlabel(x(tau);ylabel(dx/dtau);axis(-2.0,3.0,-2.0,3.0);legend(d(1,:),d(2,:),d(3,:);續(xù)上：第二十二張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)程序運(yùn)行結(jié)果第二十三張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)二、非線性阻尼系統(tǒng)1、運(yùn)動(dòng)微分方程式中，常量 為摩擦系數(shù)， 為物體的重量，k為線性彈簧的系數(shù)。干摩擦力是速度的分段函數(shù)，用signum表示。速度為正時(shí)， signum取+1，速度為負(fù)時(shí)，signum取-

13、1. 如果彈簧的彈性力不能克服干摩擦力，系統(tǒng)將停止振動(dòng)。即當(dāng)?shù)诙膹?，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)引入新變量將方程轉(zhuǎn)化一階方程形式：兩邊同時(shí)求導(dǎo)兩邊同時(shí)求導(dǎo)第二十五張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)function xdot=FrictionOscillation(t,x,d)% 非線性阻尼系統(tǒng)ode文件if abs(x(1)=d & x(2)=0.0; xdot=0;0;else xdot=x(2);-d*sign(x(2)-x(1);end2、Matlab求解編寫常微分方程對應(yīng)的函數(shù)文件FrictionOscillation.m第

14、二十六張，PPT共三十七頁，創(chuàng)作于2022年6月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);for i=1:2; t,x=ode45(FrictionOscillation,tspan,x0(i),v0(i),options,d); figure(1); plot(t,x(:,1),

15、lintyp(i,:); hold on figure(2) plot(x(:,1),x(:,2),lintyp(i,:); hold onend第二十七張，PPT共三十七頁，創(chuàng)作于2022年6月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/dtau);text(2.5,0.5,(3.0,0.0);text(4.5,0.5,

16、(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ù)上：第二十八張，PPT共三十七頁，創(chuàng)作于2022年6月6.2 單自由度系統(tǒng)程序運(yùn)行結(jié)果第二十九張，PPT共三十七頁，創(chuàng)作于2022年6月6.3 多自由度系統(tǒng)6.3.1 多自由系統(tǒng)的固有頻率問題一、力學(xué)模型二、運(yùn)動(dòng)微分方程第三十張，PPT共三十七頁，創(chuàng)作于2022年6月三、Matlab求解例6.5 三自由系統(tǒng)的振動(dòng)模態(tài)及固有頻率設(shè)k1=100N/m， k2=50N/m， m1=m2=m3=100kg。求特征值與特征向量的程序如下：k=100,-100,0;-100,1

17、50,-50;0,-50,50m=diag(100,100,100)VibrationMode,EigenValue=eig(k,m)第三十一張，PPT共三十七頁，創(chuàng)作于2022年6月附錄：ode45函數(shù)如果系數(shù)矩陣A的特征值連乘積小于零，且絕對值最大和最小的特征值之比（剛性比）很大，則稱此類方程為剛性方程 ode是Matlab專門用于解微分方程的功能函數(shù)。該求解器有變步長（variable-step）和定步長（fixed-step）兩種類型。不同類型有著不同的求解器，其中ode45求解器屬于變步長的一種，采用Runge-Kutta算法； ode45表示采用四階，五階Runge-Kutta單步

18、算法,截?cái)嗾`差為(x)3。解決的是Nonstiff(非剛性)常微分方程。 第三十二張，PPT共三十七頁，創(chuàng)作于2022年6月附錄：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為返回時(shí)間列向量；解矩陣Y的每一行對應(yīng)于T的一個(gè)元素，列數(shù)與求解變量數(shù)相等。fun為函數(shù)句柄，為根據(jù)待求解的ODE方程所編寫的ode文件（

19、odefile）；TSPANT0 TFINAL是微分系統(tǒng)yF(t,y)的積分區(qū)間；Y0為初始條件options用于設(shè)置一些可選的參數(shù)值，缺省時(shí)，相對于第一種調(diào)用格式。P1，P2，的作用是傳遞附加參數(shù)P1，P2，到ode文件。當(dāng)options缺省時(shí)，應(yīng)在相應(yīng)位置保留，以便正確傳遞參數(shù)。第三十三張，PPT共三十七頁，創(chuàng)作于2022年6月附錄：ode45函數(shù)所謂的odefile實(shí)際上是一個(gè)Matlab函數(shù)文件，一般作為整個(gè)求解程序的一個(gè)子函數(shù)，表示ode求解問題ode文件的最簡單格式必須有一個(gè)自變量t和函數(shù)y作為輸入變量，一個(gè)y的導(dǎo)函數(shù)作為輸出變量。其中自變量t不論在ode文件中是否使用都必須作為第一輸