资源描述
第七章第七章 动力学与振动动力学与振动7.1轨迹轨迹7.2单自由度系统单自由度系统7.3多自由度系统多自由度系统第七章 动力学与振动7.1轨迹17.1轨迹轨迹举例说明举例说明:重力场中有两个物体重力场中有两个物体,其中质量为其中质量为m2的物体固定的物体固定,而质量为而质量为m1的物体绕的物体绕m2做平做平面圆周运动面圆周运动.做圆周运动的做圆周运动的m1物体的轨道半径物体的轨道半径用变量用变量r表示表示,角度用变量角度用变量a表示表示.m1 r a m2 两物体系统两物体系统7.1轨迹2卫星绕地球转动时，卫星绕地球转动时，m2等于地球的质量，等于地球的质量，m1等于卫星的质量，等于卫星的质量，r为卫星球心与地球为卫星球心与地球球心间的距离。其运动轨迹由下列方程组球心间的距离。其运动轨迹由下列方程组决定：决定：式中：式中：，其中，其中t是时间变量，是时间变量，p为为物体在地球表面做圆周运动物体在地球表面做圆周运动的周期。在地球表面，r=6.373x106 m。卫星绕地球转动时，m2等于地球的质量，m1等于卫星的质量，r3用龙格用龙格库塔法可以实现求解：库塔法可以实现求解：引入新状态变量：引入新状态变量：带入前面的微分方程组，可得四个一阶微分方程。带入前面的微分方程组，可得四个一阶微分方程。用龙格库塔法可以实现求解：4建立函数文件建立函数文件orbit.mfunction 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)：建立函数文件orbit.m5由初始条件建立执行文件由初始条件建立执行文件menu71.minitcond=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,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,双曲线轨迹双曲线轨迹);由初始条件建立执行文件menu71.m6程序运行结果程序运行结果程序运行结果77.2单自由度系统单自由度系统7.2.1概述概述一一.力学模型力学模型 弹簧弹簧质量质量阻尼系统阻尼系统其中：振体质量为其中：振体质量为m，弹簧的线性系数为，弹簧的线性系数为k，非线，非线 性系数为性系数为a，阻尼系数为，阻尼系数为c，外力，外力F（t）。）。mcK，aX（t）F（t）=X（0）kf（t）7.2单自由度系统mcK，aX（t）F（t）=X（0）kf（8二二.运动微分方程运动微分方程用用x表示系统的位移，则运动微分方程为：表示系统的位移，则运动微分方程为：式中：式中：固有频率固有频率 非线性系数非线性系数 阻尼因子阻尼因子引入新变量转化状态空间方程形式：引入新变量转化状态空间方程形式：二.运动微分方程97.2.2 线性系统的自由振动线性系统的自由振动一一.运动微分方程运动微分方程当当 时，得到线性振动系统的自由时，得到线性振动系统的自由振动方程。振动方程。二二.MATLAB求解求解对应的函数文件对应的函数文件FreeOcillation.mfunction xdot=FreeOcillation(t,x,dummy,zeta)xdot=x(2);-2.0*zeta*x(2)-x(1);三种阻尼系数（三种阻尼系数（1）阻尼系数为）阻尼系数为0.1时是欠阻尼时是欠阻尼情况（情况（2）阻尼系数为）阻尼系数为1时是临界阻尼情况时是临界阻尼情况（3）阻尼系数为）阻尼系数为5时是过阻尼情况时是过阻尼情况7.2.2 线性系统的自由振动10由初始条件（位移和速度均为由初始条件（位移和速度均为1时）建立执行文件时）建立执行文件menu72.mzeta=0.1 1.0 5.0;tspan=linspace(0,40,400);%生成生成0-40的四百个线的四百个线性点性点lintype=-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 onendsubplot(2,1,1);由初始条件（位移和速度均为1时）建立执行文件menu72.m11xlabel(Time(tau);ylabel(Displacement x(tau);title(Displacement as a function of(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);xlabel(Time(tau);12程序运行结果程序运行结果程序运行结果137.2.3 线性系统的强迫振动线性系统的强迫振动一一.运动微分方程运动微分方程二二.MATLAB求解求解若若对应的函数文件对应的函数文件ForceOcillation.mfunction xdot=ForceOcillation(t,x,dummy,zeta,Omega,x0)xdot=x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t);7.2.3 线性系统的强迫振动14为了获得频谱图为了获得频谱图,建立函数文件建立函数文件AmplitudeSpectrum.mfunctionf,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=211=2048。fft的应用见的应用见Help编制执行文件编制执行文件menu72f.m为了获得频谱图,建立函数文件AmplitudeSpectru15zeta=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(:,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);zeta=0.4;Omega=3.0;x0=50;funct16yy=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(Response spectrum of a linear system);hold onsubplot(2,1,1);xlabel(Time(tau);ylabel(Displacement x(tau);title(Response of a linear system);hold onyy=x(:,1);17程序运行结果程序运行结果程序运行结果187.2.4 线性系统的频率响应、阶跃响应及线性系统的频率响应、阶跃响应及脉冲响应脉冲响应单自由度振动系统的强迫振动微分方程可为：单自由度振动系统的强迫振动微分方程可为：通过通过LAPLACE变换，得传递函数：变换，得传递函数：其中：其中：Logspace,rad2deg,loglog,semilogx see help7.2.4 线性系统的频率响应、阶跃响应及脉冲响应Logs19一一.编制执行文件编制执行文件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;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;endfor cnt=1:length(zeta)figure(1)loglog(w,mag(cnt,:),k-)一.编制执行文件frequency72.m，求频率响应。20title(SDOF frequency response magnitudes for zeta=0.2 to 1.0 in steps of 0.2)xlabel(Frequency(rad/sec)ylabel(Magnitude)gridhold onendhold offfor cnt=1:length(zeta)figure(2)semilogx(w,phs(cnt,:),k-)title(SDOF frequency response phases for zeta=0.2 to 1.0 in steps of 0.2)xlabel(Frequency(rad/sec)ylabel(phase)gridhold onendhold offtitle(SDOF frequency response21程序运行结果程序运行结果幅频曲线幅频曲线程序运行结果幅频曲线22程序运行结果程序运行结果相频曲线相频曲线程序运行结果相频曲线23二二.求时频响应的基本函数命令求时频响应的基本函数命令可以通过上述命令求线性系统的波得图、乃奎可以通过上述命令求线性系统的波得图、乃奎斯特图、斯特图、阶跃响应、脉冲响应、初始条件响应、阶跃响应、脉冲响应、初始条件响应、输入输入u的响应的响应二.求时频响应的基本函数命令可以通过上述命令求线性系统的波得24例：博得图例：博得图编制执行文件编制执行文件bode72.mm=1zeta=0.1k=1wn=sqrt(k/m)den=1 2*zeta*wn wn2num=1/msys=tf(num,den)bode(sys)Bode,nyquist see help例：博得图Bode,nyquist see help25例：求正弦输入激励响应例：求正弦输入激励响应编制执行文件编制执行文件sin72.mt=0:0.01:50;u=sin(t);lsim(sys,u,t)Lsim see help例：求正弦输入激励响应Lsim see help26INITIAL的使用的使用对方程对方程建立状态向量建立状态向量方程写为方程写为其中其中INITIAL的使用对方程建立状态向量方程写为其中27则输出为：编制执行文件编制执行文件initial72.mm=1;c=0.1;k=1;A=0 1;-k/m-c/m;C=1 0;sys=ss(A,C,);x0=10,0;initial(sys,x0)C、D确定要输出的量确定要输出的量initial see help则输出为：编制执行文件initial72.mC、D确定要输出28
展开阅读全文