MATLAB与工程应用-第7章-动力学与振动.ppt

1、第七章 动力学与振动,7.1轨迹 7.2单自由度系统 7.3多自由度系统,7.1轨迹 举例说明:重力场中有两个物体,其中质量为m2的物体固定,而质量为m1的物体绕m2做平面圆周运动.做圆周运动的m1物体的轨道半径用变量r表示,角度用变量a表示. m1 r a m2 两物体系统,卫星绕地球转动时，m2等于地球的质量，m1等于卫星的质量，r为卫星球心与地球球心间的距离。其运动轨迹由下列方程组决定： 式中： ，其中t是时间变量，p为物体在地球表面做圆周运动的周期。在地球表面，r=6.373x106 m。,用龙格库塔法可以实现求解： 引入新状态变量： 带入前面的微分方程组，可得四个一阶微分方程。,建立

2、函数文件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); 三组初始条件(t=0)：,由初始条件建立执行文件menu71.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=-. -. -.; for i=1:3 t,x=ode45(orbit,tspan,i

3、nitcond(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,双曲线轨迹);,程序运行结果,7.2单自由度系统 7.2.1概述 一.力学模型 弹簧质量阻尼系统 其中：振体质量为m，弹簧的线性系数为k，非线 性系数为a，阻尼系数为c，外力F（t）。,二.运动微分方程 用x表示系统的位移，则运动微分方程为： 式中： 固有频率 非线性系数 阻尼因子 引入新变量转化状态空间方程形式：,7.2.2 线性系统的

4、自由振动 一.运动微分方程 当 时，得到线性振动系统的自由振动方程。 二.MATLAB求解 对应的函数文件FreeOcillation.m function xdot=FreeOcillation(t,x,dummy,zeta) xdot=x(2);-2.0*zeta*x(2)-x(1); 三种阻尼系数（1）阻尼系数为0.1时是欠阻尼情况（2）阻尼系数为1时是临界阻尼情况 （3）阻尼系数为5时是过阻尼情况,由初始条件（位移和速度均为1时）建立执行文件menu72.m zeta=0.1 1.0 5.0; tspan=linspace(0,40,400);%生成0-40的四百个线性点 lintyp

5、e=-b -r -r; for i=1:3 t,x=ode45(FreeOcillation,tspan,1 1,zeta(i); subplot(2,1,1); plot(t,x(:,1),lintype(2*(i-1)+1:2*i); hold on subplot(2,1,2); plot(x(:,1),x(:,2),lintype(2*(i-1)+1:2*i); hold on end subplot(2,1,1);,xlabel(Time( tau); ylabel(Displacement x( tau); title(Displacement as a function of(

6、tau); axis(0 40 -2.0 2.0); text(2.7,-1.3,阻尼系数=0.1); text(3.6,-0.1,1.0); text(3.6,1.0,5.0); subplot(2,1,2); xlabel(Displacement); ylabel(Velocity); title(Phase portrait); axis(-2.0 2.0 -2.0 2.0); text(0.7,-1.25,阻尼系数=0.1); text(0.8,-0.65,1.0); text(0.8,0.1,5.0);,程序运行结果,7.2.3 线性系统的强迫振动 一.运动微分方程 二.MATLA

7、B求解 若 对应的函数文件ForceOcillation.m function xdot=ForceOcillation(t,x,dummy,zeta,Omega,x0) xdot=x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t);,为了获得频谱图,建立函数文件AmplitudeSpectrum.m functionf,amplitude=AmplitudeSpectrum(yy,Fs,Nstart,N); f=(Fs*(0:N-1)/N)*2.0*pi; amplitude=abs(fft(yy(Nstart:Nstart+N),N)/N; 采样速率30/60

8、00=0.005,则采样频率1/0.005=200,这个频率远远超出了必须达到的采样频率,结果显示截短频谱图,需设置Nstart=3200,N=211=2048。 fft的应用见Help 编制执行文件menu72f.m,zeta=0.4;Omega=3.0;x0=50; tspan=linspace(0,30,6000); options=odeset(RelTol,1e-8,AbsTol,1e-8); lintype=-b; t,x=ode45(ForceOcillation,tspan,0 0,options,zeta,Omega,x0); subplot(2,1,1); plot(t,x

9、(:,1); axis(0 30 -8 8); hold on subplot(2,1,2);,function xdot=ForceOcillation(t,x,dummy,zeta,Omega,x0) xdot=x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t);,yy=x(:,1); N=2048;Nstart=3200;Fs=200; f,Amplitude=AmplitudeSpectrum(yy,Fs,Nstart,N); semilogy(f(1:40),2*Amplitude(1:40); xlabel(Frequency); ylabel(Amp

10、litude); title(Response spectrum of a linear system); hold on subplot(2,1,1); xlabel(Time( tau); ylabel(Displacement x( tau); title(Response of a linear system); hold on,程序运行结果,7.2.4 线性系统的频率响应、阶跃响应及脉冲响应 单自由度振动系统的强迫振动微分方程可为： 通过LAPLACE变换，得传递函数： 其中：,Logspace, rad2deg, loglog, semilogx see help,一.编制执行文件

11、frequency72.m，求频率响应。 m = 1; zeta = 0.1:0.1:1; k = 1; wn = sqrt(k/m); 10 w = logspace(-1,1,400); rad2deg = 180/pi; s = j*w; for cnt = 1:length(zeta) xfer(cnt,:)=(1/m) ./ (s.2 + 2*zeta(cnt)*wn*s + wn2); mag(cnt,:) = abs(xfer(cnt,:); phs(cnt,:) = angle(xfer(cnt,:)*rad2deg; end for cnt = 1:length(zeta)

12、figure(1) loglog(w,mag(cnt,:),k-),title(SDOF frequency response magnitudes for zeta = 0.2 to 1.0 in steps of 0.2) xlabel(Frequency(rad/sec) ylabel(Magnitude) grid hold on end hold off for cnt = 1:length(zeta) figure(2) semilogx(w,phs(cnt,:),k-) title(SDOF frequency response phases for zeta = 0.2 to

13、1.0 in steps of 0.2) xlabel(Frequency(rad/sec) ylabel(phase) grid hold on end hold off,程序运行结果幅频曲线,程序运行结果相频曲线,二.求时频响应的基本函数命令,可以通过上述命令求线性系统的波得图、乃奎斯特图、阶跃响应、脉冲响应、初始条件响应、输入u的响应,例：博得图 编制执行文件 bode72.m m=1 zeta=0.1 k=1 wn=sqrt(k/m) den=1 2*zeta*wn wn2 num=1/m sys=tf(num,den) bode(sys),Bode,nyquist see help,例：求正弦输入激励响应 编制执行文件sin72.m t=0:0.