1、上一页目录下一页返回2F教学目标教学目标介绍介绍Matlab在动力学与振动中的应用，分别用在动力学与振动中的应用，分别用于轨迹，单自由度和多自由度线性与非线性系于轨迹，单自由度和多自由度线性与非线性系统的自由振动和强迫振动的分析。统的自由振动和强迫振动的分析。F学习要求学习要求 能够运用能够运用Matlab基本原理，对物体的运动轨迹和基本原理，对物体的运动轨迹和单自由度系统进行简单的动力学分析。单自由度系统进行简单的动力学分析。上一页目录下一页返回3目录目录6.1 轨迹轨迹6.2 单自由度系统单自由度系统6.3 多自由度系统多自由度系统习题习题上一页目录下一页返回46.1 轨迹轨迹举例说明举例
2、说明:重力场中有两个物体重力场中有两个物体,其中质量为其中质量为m2的物体固定的物体固定,而而质量为质量为m1的物体绕的物体绕m2做平面圆周运动做平面圆周运动.做圆周运动的做圆周运动的m1物体物体的轨道半径用变量的轨道半径用变量r表示表示,角度用变量角度用变量a表示表示.m2m1ar上一页目录下一页返回56.1 轨迹轨迹例例6.1：卫星绕地球转动时，：卫星绕地球转动时，m2等于地球的质量，等于地球的质量，m1等于卫等于卫星的质量，星的质量，r为卫星球心与地球球心间的距离。其运动轨迹由为卫星球心与地球球心间的距离。其运动轨迹由下列方程组决定：下列方程组决定：0242222222ddaddrdad
3、rrddardrd式中：式中：，其中，其中t是时间变量，是时间变量，p为物体在地球表面做为物体在地球表面做圆周运动的周期。在地球表面，圆周运动的周期。在地球表面，r=6.373x106 m。pt/上一页目录下一页返回66.1 轨迹轨迹用龙格用龙格库塔法可以实现求解：库塔法可以实现求解：引入新状态变量：引入新状态变量：ddaxaxddrxrx432114244321224122124xxxddxxddxxxxddxxddx0242222222ddaddrdadrrddardrd上一页目录下一页返回7建立函数文件建立函数文件Orbit.mfunction xd=Orbit(t,x)xd=x(2)x
4、(1)*x(4)2-4.0*pi2/x(1)2 x(4)-2.0*x(2)*x(4)/x(1);14244321224122124xxxddxxddxxxxddxxddx6.1 轨迹轨迹组组X1初始初始X2初始初始X3初始初始X4初始初始轨迹类型轨迹类型12001.5椭圆椭圆21002pi圆圆32004双曲线双曲线三组初始条件三组初始条件(t=0)：上一页目录下一页返回8由初始条件建立执行文件由初始条件建立执行文件execute_61.minitcond=2 0 0 1.5;1 0 0 2*pi;2 0 0 4;tspan=linspace(0,5,1000);options=odeset(R
5、elTol,1e-6,AbsTol,1e-6 1e-6 1e-6 1e-6);lintype=k-b-.r-;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,双曲线轨迹双曲线轨迹);6.1 轨迹轨迹常微分方程的常微分方程的数值求解函数数值求解函数上一页目录下一页返回9程序运行结果程序运行结果6.1 轨迹轨迹上一页目录下
6、一页返回106.2 单自由度系统单自由度系统6.2.1 概述概述一一.力学模型力学模型mcK，aX（t）F（t）=X（0）kf（t）弹簧弹簧质量质量阻尼系统阻尼系统其中：振体质量为其中：振体质量为m，弹簧的线性系数为，弹簧的线性系数为k，非线性系数为，非线性系数为a，阻尼系数为阻尼系数为c，外力，外力F（t）。）。上一页目录下一页返回116.2 单自由度系统单自由度系统二二.运动微分方程运动微分方程)(20322fXxxddxdxd用用x表示系统的位移，则运动微分方程为：表示系统的位移，则运动微分方程为：式中：式中：tn固有频率固有频率:非线性系数非线性系数:阻尼因子阻尼因子:mkn2nman
7、mc2上一页目录下一页返回126.2 单自由度系统单自由度系统引入新变量转化状态空间方程形式：引入新变量转化状态空间方程形式：ddxxxx21)(20322fXxxddxdxd)(203112221 fXxxxddxxddx上一页目录下一页返回136.2 单自由度系统单自由度系统6.2.2 线性系统的自由振动线性系统的自由振动一一.运动微分方程运动微分方程当当 时，得到线性振动系统的自由振动方程。时，得到线性振动系统的自由振动方程。0)(,0F0222xddxdxd 上一页目录下一页返回146.2 单自由度系统单自由度系统二二.MATLAB求解求解编写方程对应的函数文件编写方程对应的函数文件F
8、reeOscillation.m)(203112221 fXxxxddxxddx0三种阻尼系数三种阻尼系数()（1）阻尼系数为）阻尼系数为0.1时是欠阻尼情况（时是欠阻尼情况（2）阻尼系数为阻尼系数为1时是临界阻尼情况（时是临界阻尼情况（3）阻尼系数为）阻尼系数为5时是过阻尼时是过阻尼情况情况function xdot=FreeOscillation(t,x,zeta,Alpha)xdot=x(2);-2.0*zeta*x(2)-x(1)-Alpha*x(1)3;end上一页目录下一页返回156.2 单自由度系统单自由度系统 由初始条件（位移和速度均为由初始条件（位移和速度均为1时时,）建立执
9、行）建立执行文件文件(execute_62.m)zeta=0.1 1.0 5.0;Alpha=0.0,0.0,0.0;tspan=linspace(0,40,400);%生成生成0-40的四百个线性点的四百个线性点lintype=char(-k,-k,-.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(:,
10、2)为速度为速度 hold onend 040ddxxxx21上一页目录下一页返回166.2 单自由度系统单自由度系统figure(1);xlabel(Time(tau);ylabel(Displacement x(tau);title(Displacement 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);
11、axis(-2.0 2.0-2.0 2.0);legend(zeta=0.1,zeta=1.0,zeta=5.0);续上：续上：上一页目录下一页返回176.2 单自由度系统单自由度系统程序运行结果程序运行结果上一页目录下一页返回186.2 单自由度系统单自由度系统6.2.3 非线性系统的自由振动非线性系统的自由振动1 1、运动微分方程、运动微分方程23220d xdxxxdd一一.非线性弹簧系统非线性弹簧系统上一页目录下一页返回196.2 单自由度系统单自由度系统2、Matlab求解求解编写常微分方程对应的函数文件编写常微分方程对应的函数文件FreeOscillation.mfunction
12、xdot=FreeOscillation(t,x,zeta,Alpha)xdot=x(2);-2.0*zeta*x(2)-x(1)-Alpha*x(1)3;end23220d xdxxxdd与例与例6.2相同，只是相同，只是改变了改变了Alpha的值，的值，可以直接借用例可以直接借用例6.2的函数文件的函数文件上一页目录下一页返回206.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;
13、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);上一页目录下一页返回216.2 单自由度系统单自由度系统for i=1:3 t,x=ode45(FreeOscillation,tspan,x0(i)v0(i),options,zet
14、a,Alpha(i);figure(1)plot(t,x(:,1),lintyp(i,:);hold on figure(2)plot(x(:,1),x(:,2),lintyp(i,:);hold onend续上：续上：上一页目录下一页返回226.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/dtau);axis(-2.0,3.0,-2.0,3.0);legend(d(1,:
15、),d(2,:),d(3,:);续上：续上：上一页目录下一页返回236.2 单自由度系统单自由度系统程序运行结果程序运行结果上一页目录下一页返回246.2 单自由度系统单自由度系统二、非线性阻尼系统二、非线性阻尼系统1 1、运动微分方程、运动微分方程22()0d xdxxd signumdd式中，常量式中，常量 为摩擦系数，为摩擦系数，为物体的重为物体的重量，量，k为线性弹簧的系数。干摩擦力是速度的分段函数，用为线性弹簧的系数。干摩擦力是速度的分段函数，用signum表示。速度为正时，表示。速度为正时，signum取取+1，速度为负时，速度为负时，signum取取-1.如果弹簧的弹性力不能克服
16、干摩擦力，系统将如果弹簧的弹性力不能克服干摩擦力，系统将停止振动。即当停止振动。即当/dmg k，mg0=dxxdd，且上一页目录下一页返回256.2 单自由度系统单自由度系统引入新变量将方程转化一阶方程形式：引入新变量将方程转化一阶方程形式：ddxxxx2112212()dxxddxxdsignum xd 22()d xdxxdsignumdd 两边同时求导两边同时求导1dxdxdd两边同时求导两边同时求导222dxd xdd上一页目录下一页返回266.2 单自由度系统单自由度系统function xdot=FrictionOscillation(t,x,d)%非线性阻尼系统非线性阻尼系统o
17、de文件文件if abs(x(1)=d&x(2)=0.0;xdot=0;0;else xdot=x(2);-d*sign(x(2)-x(1);end2、Matlab求解求解编写常微分方程对应的函数文件编写常微分方程对应的函数文件FrictionOscillation.m12212()dxxddxxdsignum xd 0=dxxdd，且上一页目录下一页返回276.2 单自由度系统单自由度系统 由初始条件由初始条件(d=0.86,初始条件初始条件a（3.0，0.0），），b（5.0，0.0）)建立执行文件建立执行文件(execute_64.m)，求数值解，求数值解d=0.86;x0=3.0,5.
18、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),lintyp(i,:);hold on figure(2)plot(x(:,1),x(:,2),lintyp(i,:);hold onend上一页目录下一页返回286.2 单自由度系统单自由度系统figure(1)xlabel(tau
19、);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,(5.0,0.0);plot(-4,6,0,0,k-,0,0,-6,4,k-);axis(-4.0,6.0,-6.0,4.0);续上：续上：上一页目录下一页返回296.2 单自由度系统单自由度系统程序运行结果程序运行结果上一页目录下一页返回306.3 多
20、自由度系统多自由度系统6.3.1 多自由系统的固有频率问题多自由系统的固有频率问题一、力学模型一、力学模型二、运动微分方程二、运动微分方程上一页目录下一页返回31三、三、Matlab求解求解例例6.5 三自由系统的振动模态及固有频率三自由系统的振动模态及固有频率设设k1=100N/m，k2=50N/m，m1=m2=m3=100kg。求特征值与特征向量的程序如下：求特征值与特征向量的程序如下：k=100,-100,0;-100,150,-50;0,-50,50m=diag(100,100,100)VibrationMode,EigenValue=eig(k,m)上一页目录下一页返回32附录：附录
21、：ode45函数函数如果系数矩阵如果系数矩阵A的特征值连乘积小于零，且绝对值最大的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程和最小的特征值之比（刚性比）很大，则称此类方程为为刚性方程刚性方程 ode是是Matlab专门用于解微分方程的功能函数。该求解器专门用于解微分方程的功能函数。该求解器有变步长（有变步长（variable-step）和定步长（）和定步长（fixed-step）两种类型。）两种类型。不同类型有着不同的求解器，其中不同类型有着不同的求解器，其中ode45求解器属于变步长的求解器属于变步长的一种，采用一种，采用Runge-Kutta算法；算法；
22、ode45表示采用四阶，五阶表示采用四阶，五阶Runge-Kutta单步算法单步算法,截断误差为截断误差为(x)3。解决的是。解决的是Nonstiff(非刚性非刚性)常微分方程常微分方程。上一页目录下一页返回33附录：附录：ode45函数函数1.T,Y=ode45(fun,TSPAN,Y0)2.T,Y=ode45(fun,TSPAN,Y0,options)3.T,Y=ode45(fun,TSPAN,Y0,options,P1,P2,)4.T,Y,TE,YE,IE=ode45(fun,TSPAN,Y0,options,P1,P2,)调用格式：调用格式：说明：说明：v输出变量输出变量T为返回时间列
23、向量；解矩阵为返回时间列向量；解矩阵Y的每一行对应于的每一行对应于T的一个元素，列的一个元素，列数与求解变量数相等。数与求解变量数相等。vfun为函数句柄，为根据待求解的为函数句柄，为根据待求解的ODE方程所编写的方程所编写的ode文件（文件（odefile）；）；vTSPANT0 TFINAL是微分系统是微分系统yF(t,y)的积分区间；的积分区间；Y0为初始条件为初始条件voptions用于设置一些可选的参数值，缺省时，相对于第一种调用格式。用于设置一些可选的参数值，缺省时，相对于第一种调用格式。P1，P2，的作用是传递附加参数的作用是传递附加参数P1，P2，到到ode文件。当文件。当op
24、tions缺省时，缺省时，应在相应位置保留应在相应位置保留，以便正确传递参数。以便正确传递参数。上一页目录下一页返回34附录：附录：ode45函数函数1.所谓的所谓的odefile实际上是一个实际上是一个Matlab函数文件，一般作为函数文件，一般作为整个求解程序的一个子函数，表示整个求解程序的一个子函数，表示ode求解问题求解问题2.ode文件的最简单格式必须有一个自变量文件的最简单格式必须有一个自变量t和函数和函数y作为输作为输入变量，一个入变量，一个y的导函数作为输出变量。的导函数作为输出变量。其中自变量其中自变量t不不论在论在ode文件中是否使用都必须作为第一输入变量，文件中是否使用都
25、必须作为第一输入变量，y则则必须作为第二输入变量，位置不能颠倒。必须作为第二输入变量，位置不能颠倒。3.可以向可以向ode文件中传递参数，数目不受限制文件中传递参数，数目不受限制odefile上一页目录下一页返回35附录：附录：ode45函数函数为了能够解出方程，要用指令为了能够解出方程，要用指令odeset确定求解的条件和要求。确定求解的条件和要求。在在MATLAB中，求解方程组的指令都有默认的求解的条件和要求中，求解方程组的指令都有默认的求解的条件和要求（由结构数组（由结构数组options表示），但可以用表示），但可以用odeset修改或重新建立，修改或重新建立，odesetoption
26、s=odeset(name1,value1,name2,value2,)options=odeset(oldopts,name1,value1,name2,value2,)options=odeset(oldopts,newopts)odeset语句格式如下：语句格式如下：上一页目录下一页返回36附录：附录：ode45函数函数 第一种调用格式是指定各个参数的取值，对不指定取值的参第一种调用格式是指定各个参数的取值，对不指定取值的参数，取默认值。在不引起混淆的情况下，参数名可以只键入前数，取默认值。在不引起混淆的情况下，参数名可以只键入前面的几个字母，也不必区分大小写，如用面的几个字母，也不必区
27、分大小写，如用“abst”表示表示AbsTol.但但数值的输入必须格式正确，否则仍采用默认值。数值的输入必须格式正确，否则仍采用默认值。第二种格式使用了原来的优化选项，但对其中的参数第二种格式使用了原来的优化选项，但对其中的参数1等指等指定了新值。定了新值。第三种格式合并了两个优化选项第三种格式合并了两个优化选项oldopts newopts，重复部分，重复部分取取newopts的指定值）：的指定值）：第四种格式可在屏幕上显示如下全部可设置的参数及其默认第四种格式可在屏幕上显示如下全部可设置的参数及其默认值。值。上一页目录下一页返回37附录：附录：ode45函数函数键入键入help odeset可查看全部参数的说明，下面对其中几个参数可查看全部参数的说明，下面对其中几个参数举例说明。举例说明。RelTol 相对误差，默认值为相对误差，默认值为1e-3AbsTol 绝对误差，默认值为绝对误差，默认值为1e-6OutputFcn 输出方式，输出方式，默认值为默认值为odeplot,其它选项有：其它选项有：odeplot 按时间顺序画出全部变量的解按时间顺序画出全部变量的解 odephase2 二维相空间中两个变量的图形二维相空间中两个变量的图形 odephase2三维相空间中三个变量的图形三维相空间中三个变量的图形 odeprint 打印输出打印输出参看课本参看课本P129
