Matlab应用之动力学与振动.ppt

上一页目录下一页返回 2 F教学目标 介绍Matlab在动力学与振动中的应用，分 别用于轨迹，单自由度和多自由度线性与非线 性系统的自由振动和强迫振动的分析。 F学习要求 能够运用Matlab基本原理，对物体的运动轨迹和 单自由度系统进行简单的动力学分析。 上一页目录下一页返回 3 目录 6.1 轨迹 6.2 单自由度系统 6.3 多自由度系统 习题 上一页目录下一页返回 4 6.1 轨迹 举例说明:重力场中有两个物体,其中质量为m2的物体固定, 而质量为m1的物体绕m2做平面圆周运动.做圆周运动的m1物 体的轨道半径用变量r表示,角度用变量a表示. m2 m1 a r 上一页目录下一页返回 5 6.1 轨迹 例6.1：卫星绕地球转动时，m2等于地球的质量，m1等于 卫星的质量，r为卫星球心与地球球心间的距离。其运动轨迹 由下列方程组决定： 式中： ，其中t是时间变量，p为物体在地球表面做 圆周运动的周期。在地球表面，r=6.373x106 m。 上一页目录下一页返回 6 6.1 轨迹 用龙格库塔法可以实现求解： 引入新状态变量： 上一页目录下一页返回 7 建立函数文件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初始 轨轨迹类类 型 12001.5椭圆椭圆 21002pi圆圆 32004双曲线线 三组初始条件(t=0)： 上一页目录下一页返回 8 由初始条件建立执行文件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 轨迹 常微分方程的 数值求解函数 上一页目录下一页返回 9 程序运行结果 6.1 轨迹 上一页目录下一页返回 10 6.2 单自由度系统 6.2.1 概述 一.力学模型 m c K，a X（t） F（t）=X（0）kf（t ） 弹簧质量阻尼系统 其中：振体质量为m，弹簧的线性系数为k，非线性系数为a， 阻尼系数为c，外力F（t）。 上一页目录下一页返回 11 6.2 单自由度系统 二. 运动微分方程 用x表示系统的位移，则运动微分方程为： 式中： 固有频率: 非线性系数:阻尼因子: 上一页目录下一页返回 12 6.2 单自由度系统 引入新变量转化状态空间方程形式： 上一页目录下一页返回 13 6.2 单自由度系统 6.2.2 线性系统的自由振动 一.运动微分方程 当 时，得到线性振动系统的自由振动方程。 上一页目录下一页返回 14 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 上一页目录下一页返回 15 6.2 单自由度系统 由初始条件（位移和速度均为1时, ）建立执行 文件(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(:,2)为速度 hold on end 上一页目录下一页返回 16 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); 续上 ： 上一页目录下一页返回 17 6.2 单自由度系统 程序运行结果 上一页目录下一页返回 18 6.2 单自由度系统 6.2.3 非线性系统的自由振动 1、运动微分方程 一.非线性弹簧系统 上一页目录下一页返回 19 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的函数文件 上一页目录下一页返回 20 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); 上一页目录下一页返回 21 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 续上 ： 上一页目录下一页返回 22 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,:); 续上 ： 上一页目录下一页返回 23 6.2 单自由度系统 程序运行结果 上一页目录下一页返回 24 6.2 单自由度系统 二、非线性阻尼系统 1、运动微分方程 式中，常量 为摩擦系数， 为物体的重 量，k为线性弹簧的系数。干摩擦力是速度的分段函数，用 signum表示。速度为正时， signum取+1，速度为负时， signum取-1. 如果弹簧的弹性力不能克服干摩擦力，系统将 停止振动。即当 上一页目录下一页返回 25 6.2 单自由度系统 引入新变量将方程转化一阶方程形式： 两边同时求导 两边同时求导 上一页目录下一页返回 26 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 上一页目录下一页返回 27 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 上一页目录下一页返回 28 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); 续上 ： 上一页目录下一页返回 29 6.2 单自由度系统 程序运行结果 上一页目录下一页返回 30 6.3 多自由度系统 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,50 m=diag(100,100,100) VibrationMode,EigenValue=eig(k,m) 上一页目录下一页返回 32 附录：ode45函数 如果系数矩阵A的特征值连乘积小于零，且绝对值 最大和最小的特征值之比（刚性比）很大，则称此类 方程为刚性方程 ode是Matlab专门用于解微分方程的功能函数。该求解器 有变步长（variable-step）和定步长（fixed-step）两种类型。 不同类型有着不同的求解器，其中ode45求解器属于变步长的 一种，采用Runge-Kutta算法； 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为返回时间列向量；解矩阵Y的每一行对应于T的一个元素，列 数与求解变量数相等。 vfun为函数句柄，为根据待求解的ODE方程所编写的ode文件（odefile） ； vTSPANT0 TFINAL是微分系统yF(t,y)的积分区间；Y0为初始条件 voptions用于设置一些可选的参数值，缺省时，相对于第一种调用格式。P1 ，P2，的作用是传递附加参数P1，P2，到ode文件。当options缺省时 ，应在相应位置保留，以便正确传递参数。 上一页目录下一页返回 34 附录：ode45函数 所谓的odefile实际上是一个Matlab函数文件，一般作为整 个求解程序的一个子函数，表示ode求解问题 Z ode文件的最简单格式必须有一个自变量t和函数y作为输 入变量，一个y的导函数作为输出变量。其中自变量t不论在 ode文件中是否使用都必须作为第一输入变量，y则必须作为 第二输入变量，位置不能颠倒。 Z 可以向ode文件中传递参数，数目不受限制 odefile 上一页目录下一页返回 35 附录：ode45函数 为了能够解出方程，要用指令odeset确定求解的条件和要求。 在MATLAB中，求解方程组的指令都有默认的求解的条件和要求 （由结构数组options表示），但可以用odese