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

第七章动力学与振动7.1轨迹7.2单自由度系统7.3多自由度系统7.1轨迹举例说明:重力场中有两个物体,其中质量为m2的物体固定,而质量为m1的物体绕m2做平面圆周运动.做圆周运动的m1物体的轨道半径用变量r表示,角度用变量a表示.m1ram2两物体系统卫星绕地球转动时，m2等于地球的质量，m1等于卫星的质量，r为卫星球心与地球球心间的距离。其运动轨迹由下列方程组决定：式中：，其中t是时间变量，p为物体在地球表面做圆周运动的周期。在地球表面，r=6.373x106m。0242222222ddaddrdadrrddardrdpt/用龙格—库塔法可以实现求解：引入新状态变量：带入前面的微分方程组，可得四个一阶微分方程。ddaxaxddrxrx432114244321224122124xxxddxxddxxxxddxxddx0242222222ddaddrdadrrddardrd建立函数文件orbit.mfunctionxd=orbit(t,x)xd=[x(2);x(1)*x(4)^2-4.0*pi^2/x(1)^2;x(4);-2.0*x(2)*x(4)/x(1)];三组初始条件(t=0)：组X1初始X2初始X3初始X4初始轨迹类型12001.5椭圆21002pi圆32004双曲线由初始条件建立执行文件menu71.minitcond=[2001.5;1002*pi;2004];tspan=linspace(0,5,1000);options=odeset('RelTol',1e-6,'AbsTol',[1e-61e-61e-61e-6]);lintype=['-.''-.''-.'];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,'双曲线轨迹');程序运行结果7.2单自由度系统7.2.1概述一.力学模型弹簧—质量—阻尼系统其中：振体质量为m，弹簧的线性系数为k，非线性系数为a，阻尼系数为c，外力F（t）。mcK，aX（t）F（t）=X（0）kf（t）二.运动微分方程用x表示系统的位移，则运动微分方程为：式中：固有频率非线性系数阻尼因子引入新变量转化状态空间方程形式：)(20322fXxxddxdxdtnnmc2mkn2nmaddxxxx21)(203112221fXxxxddxxddx7.2.2线性系统的自由振动一.运动微分方程当时，得到线性振动系统的自由振动方程。二.MATLAB求解对应的函数文件FreeOcillation.mfunctionxdot=FreeOcillation(t,x,dummy,zeta)xdot=[x(2);-2.0*zeta*x(2)-x(1)];三种阻尼系数（1）阻尼系数为0.1时是欠阻尼情况（2）阻尼系数为1时是临界阻尼情况（3）阻尼系数为5时是过阻尼情况0)(F0222xddxdxd由初始条件（位移和速度均为1时）建立执行文件menu72.mzeta=[0.11.05.0];tspan=linspace(0,40,400);%生成0-40的四百个线性点lintype=['-b''--r''---r'];fori=1:3[t,x]=ode45('FreeOcillation',tspan,[11],[],zeta(i));subplot(2,1,1);plot(t,x(:,1),lintype(2*(i-1)+1:2*i));holdonsubplot(2,1,2);plot(x(:,1),x(:,2),lintype(2*(i-1)+1:2*i));holdonendsubplot(2,1,1);xlabel('Time(\tau)');ylabel('Displacementx(\tau)');title('Displacementasafunctionof(\tau)');axis([040-2.02.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('Phaseportrait');axis([-2.02.0-2.02.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线性系统的强迫振动一.运动微分方程二.MATLAB求解若对应的函数文件ForceOcillation.mfunctionxdot=ForceOcillation(t,x,dummy,zeta,Omega,x0)xdot=[x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t)];)(2022fXxddxdxd)cos()(0XF为了获得频谱图,建立函数文件AmplitudeSpectrum.mfunction[f,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/6000=0.005,则采样频率1/0.005=200,这个频率远远超出了必须达到的采样频率,结果显示截短频谱图,需设置Nstart=3200,N=2^11=2048。fft的应用见Help编制执行文件menu72f.mzeta=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,[00],options,zeta,Omega,x0);subplot(2,1,1);plot(t,x(:,1));axis([030-88]);holdonsubplot(2,1,2);functionxdot=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('Amplitude');title('Responsespectrumofalinearsystem');holdonsubplot(2,1,1);xlabel('Time(\tau)');ylabel('Displacementx(\tau)');title('Responseofalinearsystem');holdon程序运行结果7.2.4线性系统的频率响应、阶跃响应及脉冲响应单自由度振动系统的强迫振动微分方程可为：通过LAPLACE变换，得传递函数：其中：mknnmc2Logspace,rad2deg,loglog,semilogxseehelp一.编制执行文件frequency72.m，求频率响应。m=1;zeta=0.1:0.1:1;k=1;wn=sqrt(k/m);10w=logspace(-1,1,400);rad2deg=180/pi;s=j*w;forcnt=1:length(zeta)xfer(cnt,:)=(1/m)./(s.^2+2*zeta(cnt)*wn*s+wn^2);mag(cnt,:)=abs(xfer(cnt,:));phs(cnt,:)=angle(xfer(cnt,:))*rad2deg;endforcnt=1:length(zeta)figure(1)loglog(w,mag(cnt,:),'k-')title('SDOFfrequencyresponsemagnitudesforzeta=0.2to1.0instepsof0.2')xlabel('Frequency(rad/sec)')ylabel('Magnitude')gridholdonendholdoffforcnt=1:length(zeta)figure(2)semilogx(w,phs(cnt,:),'k-')title('SDOFfrequencyresponsephasesforzeta=0.2to1.0instepsof0.2')xlabel('Frequency(rad/sec)')ylabel('phase')gridholdonendholdoff程序运行结果幅频曲线程序运行结果相频曲线二.求时频响应的基本函数命令可以通过上述命令求线性系统的波得图、乃奎斯特图、阶跃响应、脉冲响应、初始条件响应、输入u的响应例：博得图编制执行文件bode72.mm=1zeta=0.1k=1wn=sqrt(k/m)den=[12*zeta*wnwn^2]num=[1/m]sys=tf(num,den)bode(sys)Bode,nyquistseehelp例：求正弦输入激励响应编制执行文件sin72.mt=0:0.01:50;u=sin(t);lsim(sys,u,t)LsimseehelpINITIAL的使用对方程建立状态向量方程写为其中则输出为：编制执行文件initial72.mm=1;c=0.1;k=1;A=[01;-k/m-c/m];C=[10];sys=ss(A,[],C,[]);x0=[10,0];initial(sys,x0)C、D确定要输出的量initialseehelp