MATLAB第7章振动.

第七章动力学与振动7.1轨迹7.2单自由度系统7.3多自由度系统§7.1轨迹举例说明:重力场中有两个物体,其中质量为m2的物体固定,而质量为m1的物体绕m2做平面圆周运动。做圆周运动的m1物体的轨道半径用变量r表示,角度用变量θ表示。两物体系统卫星绕地球转动时，m2等于地球的质量，m1等于卫星的质量，r为卫星球心与地球球心间的距离。其运动轨迹由下列方程组决定：引入状态变量:建立函数文件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)：由初始条件建立执行文件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单自由度系统一、概述1、力学模型其中：振体质量为m，弹簧的线性系数为k，非线性系数为α，阻尼系数为c，外力F(t）。弹簧—质量—阻尼系统2、运动微分方程用x表示系统的位移，则运动微分方程为：引入新变量转化状态空间方程形式：二、线性系统的自由振动1、运动微分方程当线性项为0时，得到线性振动系统的自由振动方程。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时是过阻尼情况。•由初始条件（位移和速度均为1时）建立执行文件menu72.m•zeta=[0.11.05.0];•tspan=linspace(0,40,400);•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));•holdon•subplot(2,1,2);•plot(x(:,1),x(:,2),lintype(2*(i-1)+1:2*i));•holdon•end•subplot(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');运行结果三、线性系统的强迫振动1、运动微分方程2、MATLAB求解对应的函数文件ForceOcillation.mfunctionxdot=ForceOcillation(t,x,dummy,zeta,Omega,x0)xdot=[x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t)];为了获得频谱图,建立函数文件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。•t=0:0.001:0.6;•x=sin(2*pi*50*t)+sin(2*pi*120*t);•y=x+2*randn(size(t));•plot(y(1:50))•title('SignalCorruptedwithZero-MeanRandomNoise')•xlabel('time(seconds)')fft的应用Fft变换程序Y=fft(y,512);Pyy=Y.*conj(Y)/512;f=1000*(0:256)/512;plot(f,Pyy(1:257))title('Frequencycontentofy')xlabel('frequency(Hz)')时域运行结果频域运行结果编制执行文件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);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运行结果频率响应阶跃响应脉冲响应Bode函数:g=tf([10.17.5],[10.12900]);bode(g)2342912.05.71.0)(ssssssHbode(g,{0.1,100})nyquist(g)impulse(g)脉冲响应step(g)阶跃响应一、建立系统的运动微分方程1.用刚度影响系数法建立振动微分方程000000000212122221112112121nnnnnnnnnyyykkkkkkkkkyyymmm}0{}]{[}]{[ykyM可简写为§7.3自由振动模态及固有频率例：三质量m1，m2，m3串联于弹簧k1，k2，k3上，试列出自由振动微分方程。解：求刚度系数。给x1以单位位移，x2与x3保持不动，即x1=1，x2=x3=0；要产生这样的位移状态，在各点加的力就是k11，k21，k31，显然有：k11=k1+k2，k21=－k2，k31=0（得到第一列）再令x2=1，x1=x3=0，则有k12=－k2，k22=k2+k3，k32=－k3mmmmmm1x2x3x1k2k3k11x11k12kmmm12k22k12x32kmmm32k13x33k于是得到刚度矩阵：33332222133323123222113121100kkkkkkkkkkkkkkkkkkk由于系统有对应于广义坐标x1，x2，x3的集中质量，故质量矩阵是对角阵，于是得到自由振动微分方程：00000000000321333322221321321xxxkkkkkkkkkxxxmmm2.用柔度影响系数法建立振动微分方程柔度影响系数法就是力法，它所建立方程是各点的位移协调方程。柔度影响系数cij定义：在j点作用有单位力，而其他各点没有力，所引起i点的位移。iijcjjcj点（力）i点（位移）1→cijPj→cjj1P2PiPjPnP从右图看出，当j点作用单位力时，在梁上各点产生的位移c1j，c2j，…cij，…cnj。我们暂时不管其他点，只研究i点引起位移cijPj。同样，在1点加力P1，在i点产生位移ci1·P1；这就是i点的位移协调方程，用同样方法可以得到梁上各点的位移协调方程：设在力系P1，P2……Pn作用下i点的位移是yi，那么必有：同样，在2点加力P2，在i点产生位移ci2·P2；……同样，在j点加力Pj，在i点产生位移cij·Pj；……同样，在n点加力Pn，在i点产生位移cin·Pn；niniiiiiPcPcPcPcy22111……写成矩阵形式：nnPcPcPcy12121111nnPcPcPcy22221212nnnnnnPcPcPcy2211Pcy][（1）其中nnnnnncccccccccc212222111211称为柔度矩阵。由功的互等定理可知cij=cji，因此〔c〕是对称矩阵。我们前边讲过，自由振动时，即iiiymP于是（1）式可写成：此即用柔度影响系数法建立的多自由度系统自由振动微分方程。nnnnnnnnnyyymmmcccccccccyyy212121222211121121000000或简写成：yMcy（2）讨论：就是说，刚度矩阵和柔度矩阵是互为逆矩阵。1、（2）式可写成0yMcy等式两边同乘，则1][c0][][][11ycyMcc即0][1ycyM与前边得到的比较，可见：0][ykyM1][][ck按逆阵性质：1][][kc2、方程（2）与建立的方程形式不同，实质是一样的。解：首先计算柔度系数，用莫尔法：例：一悬臂梁，固定三个集中质量，梁的抗弯刚度EI为常数，梁质量不计，求：试用柔度系数法列微分方程。dxEIMMcljiij0(单位：〔长度〕/〔力〕)〔注：1、Mi——在i加单位力的弯矩方程，Mj——在j加单位力的弯矩方程；2、原公式中号指弯矩方程不连续时，分段积分后求和。〕lll123mmm于是：EIldxEIxldxEIMcll3)(30202111EIldxEIxldxEIMcll38)2(3202202222EIldxEIxldxEIMcll39)3(33023023