Matlab应用之动力学与振动

,无忧,PPT,整理发布,单击此处编辑母版标题,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,单击此处编辑母版标题,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,上一页,目录,下一页,返回,*,无忧,PPT,整理发布,单击此处编辑母版标题,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,第6章 动力学与振动,教学目标,介绍,Matlab,在动力学与振动中的应用，分别用于轨迹，单自由度和多自由度线性与非线性系统的自由振动和强迫振动的分析。,学习要求,能够运用,Matlab,基本原理，对物体的运动轨迹和单自由度系统进行简单的动力学分析。,目录,6.1 轨迹,6.2 单自由度系统,6.3 多自由度系统,习题,6.1,轨迹,举例说明:重力场中有两个物体,其中质量为m2的物体固定,而质量为m1的物体绕m2做平面圆周运动.做圆周运动的m1物体的轨道半径用变量r表示,角度用变量a表示.,m2,m1,a,r,6.1,轨迹,例6.1：卫星绕地球转动时，m2等于地球的质量，m1等于卫星的质量，r为卫星球心与地球球心间的距离。其运动轨迹由下列方程组决定：,式中：，其中t是时间变量，p为物体在地球表面做圆周运动的周期。在地球表面，r=6.373x106 m。,6.1,轨迹,用龙格库塔法可以实现求解：,引入新状态变量：,建立函数文件Orbit.m,function xd=Orbit(t,x),xd=x(2),x(1)*x(4)2-4.0*pi2/x(1)2,x(4),-2.0*x(2)*x(4)/x(1);,6.1 轨迹,组,X1初始,X2初始,X3初始,X4初始,轨迹类型,1,2,0,0,1.5,椭圆,2,1,0,0,2pi,圆,3,2,0,0,4,双曲线,三组初始条件(t=0)：,由初始条件建立执行文件execute_61.m,initcond=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,initcond(i,:),options);,polar(x(:,3),x(:,1),lintype(2*(i-1)+1:2*i);,hold on,end,text(0.5,-1.2,椭圆轨迹);,text(-1.2,1,圆轨迹);,text(1.75,2,双曲线轨迹);,6.1 轨迹,常微分方程的数值求解函数,程序运行结果,6.1 轨迹,6.2 单自由度系统,6.2.1 概述,一.力学模型,m,c,K,，a,X（t）,F（t）=X（0）kf（t）,弹簧质量阻尼系统,其中：振体质量为m，弹簧的线性系数为k，非线性系数为a，阻尼系数为c，外力F（t）。,6.2 单自由度系统,二.运动微分方程,用x表示系统的位移，则运动微分方程为：,式中：,固有频率:,非线性系数:,阻尼因子:,6.2 单自由度系统,引入新变量转化状态空间方程形式：,6.2 单自由度系统,6.2.2 线性系统的自由振动,一.运动微分方程,当 时，得到线性振动系统的自由振动方程。,6.2 单自由度系统,二.MATLAB求解,编写方程对应的函数文件,FreeOscillation.m,0,三种阻尼系数()（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,6.2 单自由度系统,由初始条件（位移和速度均为1时,）建立执行文件,(execute_6,2.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(:,2)为速度,hold on,end,6.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);,axis(-2.0 2.0-2.0 2.0);,legend(zeta=0.1,zeta=1.0,zeta=5.0);,续上：,6.2 单自由度系统,程序运行结果,6.2 单自由度系统,6.2.3 非线性系统的自由振动,1、运动微分方程,一.非线性弹簧系统,6.2 单自由度系统,2、Matlab求解,编写常微分方程对应的函数文件,FreeOscillation.m,function 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的函数文件,6.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-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);,6.2 单自由度系统,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 on,end,续上：,6.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,:),d(2,:),d(3,:);,续上：,6.2 单自由度系统,程序运行结果,6.2 单自由度系统,二、非线性阻尼系统,1、运动微分方程,式中，常量 为摩擦系数，为物体的重量，k为线性弹簧的系数。干摩擦力是速度的分段函数，用signum表示。速度为正时，signum取+1，速度为负时，signum取-1.如果弹簧的弹性力不能克服干摩擦力，系统将停止振动。即当,6.2 单自由度系统,引入新变量将方程转化一阶方程形式：,两边同时求导,两边同时求导,6.2 单自由度系统,function xdot=FrictionOscillation(t,x,d),%非线性阻尼系统ode文件,if abs(x(1)=d,xdot=0;0;,else,xdot=x(2);-d*sign(x(2)-x(1);,end,2、Matlab求解,编写常微分方程对应的函数文件,FrictionOscillation.m,6.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);,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 on,end,6.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/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);,续上：,6.2 单自由度系统,程序运行结果,6.3 多自由度系统,6.3.1 多自由系统的固有频率问题,一、力学模型,二、运动微分方程,三、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),附录：ode45函数,如果系数矩阵A的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程为,刚性方程,ode是Matlab专门用于解微分方程的功能函数。该求解器有变步长（variable-step）和定步长（fixed-step）两种类型。不同类型有着不同的求解器，其中ode45求解器属于变步长的一种，采用Runge-Kutta算法；ode45表示采用四阶，五阶Runge-Kutta单步算法,截断误差为(x)3。解决的是,Nonstiff(非刚性)常微分方程,。,附录：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,方程所编