第6章 Matlab应用之动力学与振动概要课件

第6章动力学与振动第6章动力学与振动教学目标介绍Matlab在动力学与振动中的应用，分别用于轨迹，单自由度和多自由度线性与非线性系统的自由振动和强迫振动的分析。学习要求

能够运用Matlab基本原理，对物体的运动轨迹和单自由度系统进行简单的动力学分析。教学目标学习要求2目录6.1轨迹6.2单自由度系统6.3多自由度系统习题目录6.1轨迹36.1轨迹举例说明:重力场中有两个物体,其中质量为m2的物体固定,而质量为m1的物体绕m2做平面圆周运动.做圆周运动的m1物体的轨道半径用变量r表示,角度用变量a表示.m2m1ar6.1轨迹举例说明:重力场中有两个物体,其中质量为m2的46.1轨迹例6.1：卫星绕地球转动时，m2等于地球的质量，m1等于卫星的质量，r为卫星球心与地球球心间的距离。其运动轨迹由下列方程组决定：式中：，其中t是时间变量，p为物体在地球表面做圆周运动的周期。在地球表面，r=6.373x106m。6.1轨迹例6.1：卫星绕地球转动时，m2等于地球的质量56.1轨迹用龙格—库塔法可以实现求解：引入新状态变量：6.1轨迹用龙格—库塔法可以实现求解：6建立函数文件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)：建立函数文件Orbit.m6.1轨迹组X1初始X2初始X7由初始条件建立执行文件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轨迹常微分方程的数值求解函数由初始条件建立执行文件execute_61.m6.1轨迹8程序运行结果6.1轨迹程序运行结果6.1轨迹96.2单自由度系统6.2.1概述一.力学模型mcK，aX（t）F（t）=X（0）kf（t）弹簧—质量—阻尼系统其中：振体质量为m，弹簧的线性系数为k，非线性系数为a，阻尼系数为c，外力F（t）。6.2单自由度系统6.2.1概述一.力学模型mcK，106.2单自由度系统二.运动微分方程用x表示系统的位移，则运动微分方程为：式中：固有频率:非线性系数:阻尼因子:6.2单自由度系统二.运动微分方程用x表示系统的位移，116.2单自由度系统引入新变量转化状态空间方程形式：6.2单自由度系统引入新变量转化状态空间方程形式：126.2单自由度系统6.2.2线性系统的自由振动一.运动微分方程当时，得到线性振动系统的自由振动方程。6.2单自由度系统6.2.2线性系统的自由振动一.运136.2单自由度系统二.MATLAB求解编写方程对应的函数文件FreeOscillation.m0三种阻尼系数()（1）阻尼系数为0.1时是欠阻尼情况（2）阻尼系数为1时是临界阻尼情况（3）阻尼系数为5时是过阻尼情况functionxdot=FreeOscillation(t,x,zeta,Alpha)xdot=[x(2);-2.0*zeta*x(2)-x(1)-Alpha*x(1)^3];end6.2单自由度系统二.MATLAB求解编写方程对应的函数146.2单自由度系统由初始条件（位移和速度均为1时,）建立执行文件(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)为速度holdonend6.2单自由度系统由初始条件（位移和速度均为1时,156.2单自由度系统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');续上：6.2单自由度系统figure(1);续上：166.2单自由度系统程序运行结果6.2单自由度系统程序运行结果176.2单自由度系统6.2.3非线性系统的自由振动1、运动微分方程一.非线性弹簧系统6.2单自由度系统6.2.3非线性系统的自由振动1、186.2单自由度系统2、Matlab求解编写常微分方程对应的函数文件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的函数文件6.2单自由度系统2、Matlab求解编写常微分方程对应196.2单自由度系统由初始条件建立执行文件(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.2单自由度系统由初始条件建立执行文件(execute206.2单自由度系统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续上：6.2单自由度系统fori=1:3续上：216.2单自由度系统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,:));续上：6.2单自由度系统figure(1)续上：226.2单自由度系统程序运行结果6.2单自由度系统程序运行结果236.2单自由度系统二、非线性阻尼系统1、运动微分方程式中，常量为摩擦系数，为物体的重量，k为线性弹簧的系数。干摩擦力是速度的分段函数，用signum表示。速度为正时，signum取+1，速度为负时，signum取-1.如果弹簧的弹性力不能克服干摩擦力，系统将停止振动。即当6.2单自由度系统二、非线性阻尼系统1、运动微分方程式中246.2单自由度系统引入新变量将方程转化一阶方程形式：两边同时求导两边同时求导6.2单自由度系统引入新变量将方程转化一阶方程形式：两边256.2单自由度系统functionxdot=FrictionOscillation(t,x,d)%非线性阻尼系统ode文件ifabs(x(1))<=d&&x(2)==0.0;xdot=[0;0];elsexdot=[x(2);-d*sign(x(2))-x(1)];end2、Matlab求解编写常微分方程对应的函数文件FrictionOscillation.m6.2单自由度系统functionxdot=Frict266.2单自由度系统由初始条件(d=0.86,初始条件a（3.0，0.0），b（5.0，0.0）)建立执行文件(execute_64.m)，求数值解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,:));holdonend6.2单自由度系统由初始条件(d=0.86,初始条件a（276.2单自由度系统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]);续上：6.2单自由度系统figure(1)续上：286.2单自由度系统程序运行结果6.2单自由度系统程序运行结果296.3多自由度系统6.3.1多自由系统的固有频率问题一、力学模型二、运动微分方程6.3多自由度系统6.3.1多自由系统的固有频率问题一、30三、Matlab求解例6.5三自由系统的振动模态及固有频率设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)三、Matlab求解例6.5三自由系统的振动模态及固有频率31附录：ode45函数如果系数矩阵A的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程为刚性方程ode是Matlab专门用于解微分方程的功能函数。该求解器有变步长（variable-step）和定步长（fixed-step）两种类型。不同类型有着不同的求解器，其中ode45求解器属于变步长的一种，采用Runge-Kutta算法；ode45表示采用四阶，五阶Runge-Kutta单步算法,截断误差为(Δx)^3。解决的是Nonstiff(非刚性)常微分方程。附录：ode45函数如果系数矩阵A的特征值连乘积小于零，且绝32附录：ode45函数[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,…)调用格式：说明：输出变量T为返回时间列向量；解矩阵Y的每一行对应于T的一个元素，列数与求解变量数相等。@fun为函数句柄，为根据待求解的ODE方程所编写的ode文件（odefile）；TSPAN＝[T0TFINAL]是微分系统y'＝F(t,y)的积分区间；Y0为初始条件options用于设置一些可选的参数值，缺省时，相对于第一种调用格式。P1，P2，…的作用是传递附加参数P1，P2，…到ode文件。当options缺省时，应在相应位置保留[]，以便正确传递参数。附录：ode45函数[T,Y]=ode45(@fun,TS33附录：ode45函数所谓的odefile实际上是一个Matlab函数文件，一般作为整个求解程序的一个子函数，表示ode求解问题ode文件的最简单格式必须有一个自变量t和函数y作为输入变量，一个y的导函数作为输出变量。其中自变量t不论在ode文件中是否使用都必须作为第一输入变量，y则必须作为第二输入变量，位置不能颠倒。可以向ode文件中传递参数，数目不受限制odefile附录：ode45函数所谓的odefile实际上是一个Matl34附录：ode45函数为了能够解出方程，要用指令odeset确定求解的条件和要求。在MATLAB中，求解方程组的指令都有默认的求解的条件和要求（由结构数组options表示），但可以用odeset修改或重新建立，odesetoptions=odeset(‘name1’,value1,’name2’,value2,…)

options=odeset(oldopts,‘name1’,value1,’name2’,value2,…)

options=odeset(oldopts,newopts)

odeset语句格式如下：附录：ode45函数为了能够解出方程，要用指令odeset确35附录：ode45函数第一种调用格式是指定各个参数的取值，对不指定取值的参数，取默认值。在不引起混淆的情况下，参数名可以只键入前面的几个字母，也不必区分大小写，如用“abst”表示AbsTol.但数值的输入必须格式正确，否则仍采用默认值。第二种格式使用了原来的优化选项，但对其中的参数1等指定了新值。第三种格式合并了两个优化选项oldoptsnewopts，重复部分取newopts的指定值）：第四种格式可在屏幕上显示如下全部可设置的参数及其默认值。附录：ode45函数第一种调用格式是指定各个参数的取值，36附录：ode45函数键入helpodeset可查看全部参数的说明，下面对其中几个参数举例说明。

RelTol

相对误差，默认值为1e-3

AbsTol

绝对误差，默认值为1e-6

OutputFcn

输出方式，

默认值为‘odeplot’,其它选项有：

odeplot

按时间顺序画出全部变量的解

odephase2

二维相空间中两个变量的图形

odephase2三维相空间中三个变量的图形

odeprint

打印输出

参看课本P129附录：ode45函数键入helpodeset可查看全部参数37第6章动力学与振动第6章动力学与振动教学目标介绍Matlab在动力学与振动中的应用，分别用于轨迹，单自由度和多自由度线性与非线性系统的自由振动和强迫振动的分析。学习要求

能够运用Matlab基本原理，对物体的运动轨迹和单自由度系统进行简单的动力学分析。教学目标学习要求39目录6.1轨迹6.2单自由度系统6.3多自由度系统习题目录6.1轨迹406.1轨迹举例说明:重力场中有两个物体,其中质量为m2的物体固定,而质量为m1的物体绕m2做平面圆周运动.做圆周运动的m1物体的轨道半径用变量r表示,角度用变量a表示.m2m1ar6.1轨迹举例说明:重力场中有两个物体,其中质量为m2的416.1轨迹例6.1：卫星绕地球转动时，m2等于地球的质量，m1等于卫星的质量，r为卫星球心与地球球心间的距离。其运动轨迹由下列方程组决定：式中：，其中t是时间变量，p为物体在地球表面做圆周运动的周期。在地球表面，r=6.373x106m。6.1轨迹例6.1：卫星绕地球转动时，m2等于地球的质量426.1轨迹用龙格—库塔法可以实现求解：引入新状态变量：6.1轨迹用龙格—库塔法可以实现求解：43建立函数文件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)：建立函数文件Orbit.m6.1轨迹组X1初始X2初始X44由初始条件建立执行文件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轨迹常微分方程的数值求解函数由初始条件建立执行文件execute_61.m6.1轨迹45程序运行结果6.1轨迹程序运行结果6.1轨迹466.2单自由度系统6.2.1概述一.力学模型mcK，aX（t）F（t）=X（0）kf（t）弹簧—质量—阻尼系统其中：振体质量为m，弹簧的线性系数为k，非线性系数为a，阻尼系数为c，外力F（t）。6.2单自由度系统6.2.1概述一.力学模型mcK，476.2单自由度系统二.运动微分方程用x表示系统的位移，则运动微分方程为：式中：固有频率:非线性系数:阻尼因子:6.2单自由度系统二.运动微分方程用x表示系统的位移，486.2单自由度系统引入新变量转化状态空间方程形式：6.2单自由度系统引入新变量转化状态空间方程形式：496.2单自由度系统6.2.2线性系统的自由振动一.运动微分方程当时，得到线性振动系统的自由振动方程。6.2单自由度系统6.2.2线性系统的自由振动一.运506.2单自由度系统二.MATLAB求解编写方程对应的函数文件FreeOscillation.m0三种阻尼系数()（1）阻尼系数为0.1时是欠阻尼情况（2）阻尼系数为1时是临界阻尼情况（3）阻尼系数为5时是过阻尼情况functionxdot=FreeOscillation(t,x,zeta,Alpha)xdot=[x(2);-2.0*zeta*x(2)-x(1)-Alpha*x(1)^3];end6.2单自由度系统二.MATLAB求解编写方程对应的函数516.2单自由度系统由初始条件（位移和速度均为1时,）建立执行文件(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)为速度holdonend6.2单自由度系统由初始条件（位移和速度均为1时,526.2单自由度系统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');续上：6.2单自由度系统figure(1);续上：536.2单自由度系统程序运行结果6.2单自由度系统程序运行结果546.2单自由度系统6.2.3非线性系统的自由振动1、运动微分方程一.非线性弹簧系统6.2单自由度系统6.2.3非线性系统的自由振动1、556.2单自由度系统2、Matlab求解编写常微分方程对应的函数文件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的函数文件6.2单自由度系统2、Matlab求解编写常微分方程对应566.2单自由度系统由初始条件建立执行文件(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.2单自由度系统由初始条件建立执行文件(execute576.2单自由度系统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续上：6.2单自由度系统fori=1:3续上：586.2单自由度系统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,:));续上：6.2单自由度系统figure(1)续上：596.2单自由度系统程序运行结果6.2单自由度系统程序运行结果606.2单自由度系统二、非线性阻尼系统1、运动微分方程式中，常量为摩擦系数，为物体的重量，k为线性弹簧的系数。干摩擦力是速度的分段函数，用signum表示。速度为正时，signum取+1，速度为负时，signum取-1.如果弹簧的弹性力不能克服干摩擦力，系统将停止振动。即当6.2单自由度系统二、非线性阻尼系统1、运动微分方程式中616.2单自由度系统引入新变量将方程转化一阶方程形式：两边同时求导两边同时求导6.2单自由度系统引入新变量将方程转化一阶方程形式：两边626.2单自由度系统functionxdot=FrictionOscillation(t,x,d)%非线性阻尼系统ode文件ifabs(x(1))<=d&&x(2)==0.0;xdot=[0;0];elsexdot=[x(2);-d*sign(x(2))-x(1)];end2、Matlab求解编写常微分方程对应的函数文件FrictionOscillation.m6.2单自由度系统functionxdot=Frict636.2单自由度系统由初始条件(d=0.86,初始条件a（3.0，0.0），b（5.0，0.0）)建立执行文件(execute_64.m)，求数值解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,:));holdonend6.2单自由度系统由初始条件(d=0.86,初始条件a（646.2单自由度系统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]);续上：6.2单自由度系统figure(1)续上：656.2单自由度系统程序运行结果6.2单自由度系统程序运行结果666.3多自由度系统6.3.1多自由系统的固有频率问题一、力学模型二、运动微分方程6.3多自由度系统6.3.1多自由系统的固有频率问题一、67三、Matlab求解例6.5三自由系统的振动模态及固有频率设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)三、Matlab求解例6.5三自由系统的振动模态及固有频率68附录：ode45函数如果系数矩阵A的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程为刚性方程ode是Matlab专门用于解微分方程的功能函数。该求解器有变步长（variable-step）和定步长（fixed-step）两种类型。不同类型有着不同的求解器，其中ode45求解器属于变步长的一种，采用Runge-Kutta算法；ode45表示采用四阶，五阶Runge-Kutta单步算法,截断误差为(Δx)^3。解决的是Nonstiff(非刚性)常微分方程。附录：ode45函数如果系数矩阵A的特征值连乘积小于零，且绝69附录：ode45函数[T,Y]=ode45(@fun,TSPAN,Y0)[T,Y]=ode45(@fun