MATLAB第章 振动

1、7.1 7.1 轨迹轨迹7.2 7.2 单自由度系统单自由度系统7.3 7.3 多自由度系统多自由度系统 7.1 7.1 轨迹轨迹举例说明举例说明:重力场中有两个物体重力场中有两个物体,其中质量为其中质量为m2的物体固定的物体固定,而质而质量为量为m1的物体绕的物体绕m2做平面圆周运动。做圆周运动的做平面圆周运动。做圆周运动的m1物体的轨物体的轨道半径用变量道半径用变量r表示表示,角度用变量角度用变量表示。表示。两物体系统两物体系统 卫星绕地球转动时，卫星绕地球转动时，m2等于地球的质量，等于地球的质量，m1等于卫星的质等于卫星的质量，量，r为卫星球心与地球球心间的距离。其运动轨迹由下列方程为

2、卫星球心与地球球心间的距离。其运动轨迹由下列方程组决定：组决定：引入状态变量引入状态变量:建立函数文件建立函数文件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)：由初始条件建立执行文件由初始条件建立执行文件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

3、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,双曲线轨迹双曲线轨迹);运行结果运行结果 7.2 7.2 单自由度系统单自由度系统一、概述一、概述1、力学模型、力学模型其中：振体质量为其中：振体质量为m，弹簧的线性系数为，弹簧的线性系数为k，非线性系数为，非线

4、性系数为，阻尼系数为阻尼系数为c，外力，外力F(t）。）。弹簧弹簧质量质量阻尼系统阻尼系统2、运动微分方程、运动微分方程用用x表示系统的位移，则运动微分方程为：表示系统的位移，则运动微分方程为：引入新变量转化状态空间方程形式：引入新变量转化状态空间方程形式：二、线性系统的自由振动二、线性系统的自由振动1、运动微分方程、运动微分方程当线性项为当线性项为0时，得到线性振动系统的自由振动方程。时，得到线性振动系统的自由振动方程。2、MATLAB求解求解对应的函数文件对应的函数文件FreeOcillation.mfunction xdot=FreeOcillation(t,x,dummy,zeta)x

5、dot=x(2);-2.0*zeta*x(2)-x(1);三种阻尼系数三种阻尼系数（1）阻尼系数为）阻尼系数为0.1时是欠阻尼情况；（时是欠阻尼情况；（2）阻）阻尼系数为尼系数为1时是临界阻尼情况；（时是临界阻尼情况；（3）阻尼系数为）阻尼系数为5时是过阻尼时是过阻尼情况。情况。由初始条件（位移和速度均为由初始条件（位移和速度均为1时）建立执行文件时）建立执行文件menu72.mzeta=0.1 1.0 5.0;tspan=linspace(0,40,400);lintype=-b -r -r;for i=1:3 t,x=ode45(FreeOcillation,tspan,1 1,zeta(

6、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);xlabel(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,

7、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);运行结果运行结果三、线性系统的强迫振动三、线性系统的强迫振动1、运动微分方程、运动微分方程2、MATLAB求解求解对应的函数文件对应的函数文件ForceOcillation.mfunction xdot=ForceOcil

8、lation(t,x,dummy,zeta,Omega,x0)xdot=x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t);为了获得频谱图为了获得频谱图,建立函数文件建立函数文件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,这个频率远远这个频率远远

9、超出了必须达到的采样频率超出了必须达到的采样频率,结果显示截短频谱图结果显示截短频谱图,需设置需设置Nstart=3200,N=211=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(Signal Corrupted with Zero-Mean Random Noise) xlabel(time (seconds)fft的应用的应用Fft变换程序变换程序Y = fft(y,512);Pyy = Y.* conj(Y) / 512;f =

10、 1000*(0:256)/512;plot(f,Pyy(1:257)title(Frequency content of y)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,0 0,options,zeta,Omega,x0); subplot

11、(2,1,1); plot(t,x(:,1); axis(0 30 -8 8); hold on subplot(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(Response spectrum of a linear system); hold onsubplot(2,1,1);xlabel(Time(

12、tau);ylabel(Displacement x( tau);title(Response of a linear system);hold on运行结果运行结果频率响应频率响应 阶跃响应阶跃响应 脉冲响应脉冲响应Bode函数函数:g = tf(1 0.1 7.5,1 0.12 9 0 0);bode(g)2342912. 05 . 71 . 0)(ssssssHbode(g,0.1 , 100)nyquist(g)impulse(g) 脉冲响应脉冲响应step(g) 阶跃响应阶跃响应一、建立系统的运动微分方程一、建立系统的运动微分方程1 1. .用刚度影响系数法建立振动微分方程用刚度影响

13、系数法建立振动微分方程000000000212122221112112121 nnnnnnnnnyyykkkkkkkkkyyymmm0ykyM 可简写为可简写为 7.3 7.3 自由振动模态及固有频率自由振动模态及固有频率例：例：三质量三质量m1，m2，m3串联于弹簧串联于弹簧k1，k2，k3上，试列出自上，试列出自由振动微分方程。由振动微分方程。 解：解：求刚度系数。求刚度系数。 给给x1以单位位移，以单位位移，x2与与x3保持保持不动，即不动，即x1=1，x2=x3=0； 要产生这样的位移状态，在各要产生这样的位移状态，在各点加的力就是点加的力就是k11，k21，k31，显，显然有：然有：

14、 k11=k1+k2，k21=k2，k31=0（得到第一列）（得到第一列）再令再令x2=1x2=1，x1=x3=0 x1=x3=0，则有则有 k12=k12=k2k2，k22=k2+k3k22=k2+k3，k32=k32=k3k3mmmmmm1x2x3x1k2k3k11x11k12kmmm12k22k12x32kmmm32k13x33k于是得到刚度矩阵：于是得到刚度矩阵： 33332222133323123222113121100kkkkkkkkkkkkkkkkkkk 由于系统有对应于广义坐标由于系统有对应于广义坐标x1，x2，x3的集中质量，故质的集中质量，故质量矩阵是对角阵，于是得到自由振

15、动微分方程：量矩阵是对角阵，于是得到自由振动微分方程： 00000000000321333322221321321xxxkkkkkkkkkxxxmmm 2.2.用柔度影响系数法建立振动微分方程用柔度影响系数法建立振动微分方程 柔度影响系数法就是力法，柔度影响系数法就是力法，它所建立方程是各点的位移协调它所建立方程是各点的位移协调方程。柔度影响系数方程。柔度影响系数cij定义：在定义：在j点作用有单位力，而其他各点没点作用有单位力，而其他各点没有力，所引起有力，所引起i点的位移。点的位移。 iijcjjcj点（力） i点（位移） 1 cij Pj cjj1P2PiPjPnP 从右图看出，当从右图

16、看出，当j点作用单位力点作用单位力时，在梁上各点产生的位移时，在梁上各点产生的位移c1j，c2j，cij，cnj。我们暂时不管其他。我们暂时不管其他点，只研究点，只研究i点引起位移点引起位移cijPj。 同样，在同样，在1点加力点加力P1，在，在i点产生位移点产生位移ci1P1； 这就是这就是i点的位移协调方程，用同样方法可以点的位移协调方程，用同样方法可以得到梁上各点的位移协调方程：得到梁上各点的位移协调方程： 设在力系设在力系P1，P2Pn作用下作用下i点的位移是点的位移是yi，那么必有：那么必有： 同样，在同样，在2点加力点加力P2，在，在i点产生位移点产生位移ci2P2； 同样，在同样

17、，在j点加力点加力Pj，在，在i点产生位移点产生位移cijPj； 同样，在同样，在n点加力点加力Pn，在，在i点产生位移点产生位移cinPn； niniiiiiPcPcPcPcy22111 写成矩阵形式：写成矩阵形式： nnPcPcPcy12121111nnPcPcPcy22221212nnnnnnPcPcPcy2211 Pcy（1） 其中其中 nnnnnncccccccccc212222111211 称为柔度矩阵。由功的互等定理可知称为柔度矩阵。由功的互等定理可知cij=cji，因此因此c是对称矩阵。我们前边讲过，自由振动是对称矩阵。我们前边讲过，自由振动时时 ，即，即 iiiymP 于是（

18、于是（1）式可写成：）式可写成： 此即用柔度影响系数法建立的多自由度系统自此即用柔度影响系数法建立的多自由度系统自由振动微分方程。由振动微分方程。 nnnnnnnnnyyymmmcccccccccyyy 212121222211121121000000或简写成：或简写成： yMcy （2） 讨论讨论 ：就是说，刚度矩阵和柔度矩阵是互为逆矩阵。就是说，刚度矩阵和柔度矩阵是互为逆矩阵。 1、（、（2）式可写成）式可写成 0yMcy 等式两边同乘等式两边同乘 ，则，则 1c 011ycyMcc 即即 01ycyM 与前边得到的与前边得到的 比较，可见：比较，可见： 0ykyM 1 ck按逆阵性质：按

19、逆阵性质： 1 kc 2、方程（、方程（2）与建立的方程形式不同，实质是）与建立的方程形式不同，实质是一样的。一样的。 解：首先计算柔度系数，解：首先计算柔度系数，用莫尔法：用莫尔法： 例：一悬臂梁，固定三个集中质量，梁的抗弯例：一悬臂梁，固定三个集中质量，梁的抗弯刚度刚度EI为常数，梁质量不计，求：试用柔度系数法为常数，梁质量不计，求：试用柔度系数法列微分方程。列微分方程。 dxEIMMcljiij0(单位：长度单位：长度/力力) 注：注：1、Mi在在i加单位力的弯矩方程，加单位力的弯矩方程，Mj在在j加单位力的弯矩方程；加单位力的弯矩方程；2、原公式中、原公式中 号号指弯矩方程不连续时，分

20、段积分后求和。指弯矩方程不连续时，分段积分后求和。 lll123mmm于是：于是： EIldxEIxldxEIMcll3)(30202111EIldxEIxldxEIMcll38)2(3202202222EIldxEIxldxEIMcll39)3(3302302333EIldxEIxlxldxEIMMccll65)2)(300212112EIldxEIxlxldxEIMMccll43)3)(300313113EIldxEIxlxldxEIMMccll314)3)(2(32020323223故柔度矩阵为故柔度矩阵为 即即 自由振动微分方程：自由振动微分方程： 542882816585263EIl

21、c321332100000054288281658526yyymmmEIlyyy 321332154288281658526yyyEImlyyy 关于固有频率及固有振型的概念，我们已经学关于固有频率及固有振型的概念，我们已经学习过。研究多自由度系统的这些固有特征时习过。研究多自由度系统的这些固有特征时, , 采用采用矩阵为数学工具。矩阵为数学工具。 二、固有频率与固有振型二、固有频率与固有振型 注意注意: :对于一个正定系统（即系统没有刚体位对于一个正定系统（即系统没有刚体位移），用上述各种方法建立微分方程都可以，有时移），用上述各种方法建立微分方程都可以，有时用柔度系统法更为方便；但对于半正

22、定系统（系统用柔度系统法更为方便；但对于半正定系统（系统有刚体位移），则柔度系统没有意义（加单位力后有刚体位移），则柔度系统没有意义（加单位力后系统不能维持平衡而产生刚体运动），只能用形如系统不能维持平衡而产生刚体运动），只能用形如My+ky=0那样的方程。那样的方程。 y 对于有对于有n个自由度个自由度q1，q2qn的系统，其自由的系统，其自由振动微分方程的一般形式是：振动微分方程的一般形式是： 或写成矩阵形式或写成矩阵形式 这是一组二阶常系数、线性常微分方程。设系这是一组二阶常系数、线性常微分方程。设系统偏离平衡位置作自由振动时，存在各统偏离平衡位置作自由振动时，存在各qi按同一频按同一频

23、率率，同一相位，同一相位作简谐振动的特解：作简谐振动的特解： 0qkqM （1） )(11tieAq)(22tieAq )(tinneAq)(tieAq 这是一组关于振幅这是一组关于振幅A1，A2An的线性齐次的线性齐次代数方程组，要使代数方程组，要使A1，A2An有非零解（零解对有非零解（零解对应不振动），必有系数行列式为零，即：应不振动），必有系数行列式为零，即： 代入方程组（代入方程组（1），消去公因子），消去公因子 后，得后，得 )(tie0)()()(121212212111211nnnAmkAmkAmk0)()()2(222222222111221nnnAmkAmkAmk0)()2

24、()(222212222121221nnnAmkAmkAmk（2） 或或 0)(2AkM02kM或或 此即多自由度系统的特征方程（或称频率此即多自由度系统的特征方程（或称频率方程）方程） 022221211211221112211nnnnnnnnnnmkmkmkmkmkmk 在这个方程中，在这个方程中，kij，mij都是已知的，未知的只都是已知的，未知的只有有，展开后是一个关于，展开后是一个关于2 2的的n n次代数方程，从中次代数方程，从中可以求得可以求得2 2的的n n个根。这些根就是系统的固有频率。个根。这些根就是系统的固有频率。由小到大排列：由小到大排列：1 1，2 2n n叫做叫做基

25、阶，二基阶，二阶阶nn阶固有频率。阶固有频率。 把每一个固有频率例如把每一个固有频率例如代回方程组（代回方程组（2 2），就可），就可以求得该阶固有振型：以求得该阶固有振型： )1()1(3)1(2)1(1:nAAAA 这里解释一下。在二自由度系统时，求得了固有这里解释一下。在二自由度系统时，求得了固有频率频率1 1，2 2之后，把之后，把1 1（或（或2 2）代回振型方程组求）代回振型方程组求该阶振型时，是代回两个方程中的任意一个，求得该阶振型时，是代回两个方程中的任意一个，求得 ，及，及 。就是说，振型方程组的两个方程不。就是说，振型方程组的两个方程不独立，而是线性相关。对于多自由度系统，

26、也是一样，独立，而是线性相关。对于多自由度系统，也是一样，因为振型方程组（因为振型方程组（2 2）的系统行列式为零，即）的系统行列式为零，即，就说明矩阵，就说明矩阵 的秩的秩rnrn，至少有一个方程与，至少有一个方程与其他方程线性相关。我们在求固有振型时，应划出一其他方程线性相关。我们在求固有振型时，应划出一个方程，例如最后一个，而把剩下的个方程，例如最后一个，而把剩下的n-1n-1个方程中某一个方程中某一项移到等式右边，例如把项移到等式右边，例如把A A1 1项移过去。那么解出的项移过去。那么解出的A A2 2，A A3 3AAn n都包含都包含A A1 1，也就是说，求出了，也就是说，求出

27、了A A2 2，A A3 3AAn n各各自对自对A A1 1的比值，那么它们之间的比值自然也就确定了。的比值，那么它们之间的比值自然也就确定了。 1)1(1)1(2AA2)2(1)2(2AA02kM)(2kM 即该阶固有振型：即该阶固有振型： 显然：显然： 代回我们设的特解代回我们设的特解 )1()1(3)1(2)1(1:nAAAA注：上标（注：上标（1）表示第一阶振型）表示第一阶振型 )()1(1)1(1tieAq)()1(2)1(2tieAq )()1()1(tinneAq或或 第第1阶主振动阶主振动 )()1()1(tieAq)1()1(2)1(1)1()1(2)1(1:nnAAAqq

28、q第第1 1阶振型：阶振型： 同样，当把同样，当把2 2，3 3n n分别代入振型方程分别代入振型方程（2 2）时，可求得其他各阶固有振型。即：）时，可求得其他各阶固有振型。即： 就是说，当系统按第一阶固有频率作振动时，就是说，当系统按第一阶固有频率作振动时，各点振幅各点振幅 之间具有确定的比值，而且之间具有确定的比值，而且系统各点之间在任一瞬时也具有同样确定的比值，系统各点之间在任一瞬时也具有同样确定的比值，这说明系统有一定的振动形态，称为第一阶固有振这说明系统有一定的振动形态，称为第一阶固有振型。型。 )1()1(2)1(1,nAAA)1()1(2)1(1)1()1(2)1(1:nnAAA

29、qqq第第2 2阶振型：阶振型： )2()2(2)2(1)2()2(2)2(1:nnAAAqqq 第第n n阶振型：阶振型： )()(2)(1)()(2)(1:nnnnnnnnAAAqqq于是也可写出各阶主振动。于是也可写出各阶主振动。 方程（方程（1）的通解是这）的通解是这n个主振动的叠加：个主振动的叠加： 写成矩阵形式：写成矩阵形式： )()(1)()2(1)()1(1)1(121ntintitieAeAeAq)()(2)()2(2)()1(2)1(221ntintitieAeAeAq )()()()2()()1()1(21ntinntintinneAeAeAq )(tieAq（3） 式中

30、共有式中共有2n个特定常数个特定常数 ：及：及，注：对于每一阶振型来说，各阶振幅之间比值，注：对于每一阶振型来说，各阶振幅之间比值已定，只要知道其中一个就可求得其他振幅，所以，已定，只要知道其中一个就可求得其他振幅，所以，每一振型中只有一个待定振幅。这些常数由初始每一振型中只有一个待定振幅。这些常数由初始条件（条件（t=0时的时的 及及 ）来确定。）来确定。 ) 1 () 1 (2) 1 (1,nAAAn21,nqqq21,nqqq21, 由（由（3）式可见，具有）式可见，具有n个自由度的无阻尼系统个自由度的无阻尼系统自由振动一般由自由振动一般由n个不同频率的主振动组成。在每个不同频率的主振动

31、组成。在每一主振动中，系统各点的位移始终保持一定的比例一主振动中，系统各点的位移始终保持一定的比例关系（等于各点振幅之比），即系统有一定的振动关系（等于各点振幅之比），即系统有一定的振动形态，称为主振型（或固有振型）。可以任取出其形态，称为主振型（或固有振型）。可以任取出其中一点的振幅（例如中一点的振幅（例如A1）等于）等于1，使其标准化，这，使其标准化，这样的振型，称为标准化了的固有振型。样的振型，称为标准化了的固有振型。 解：取三个扭转角解：取三个扭转角1 1，2 2，3 3为广义坐为广义坐标。标。 例一例一、悬臂轴抗扭刚度、悬臂轴抗扭刚度GJ，质量不计，其上，质量不计，其上固定三个圆盘，

32、转动惯量是固定三个圆盘，转动惯量是I。求：系统的固有频。求：系统的固有频率，固有振型。率，固有振型。 123代入拉氏方程代入拉氏方程 232221212121IIIT22321221)(21)(2121lGJlGJlGJU0)(iiUTdtd得系统自由振动微分方程：得系统自由振动微分方程： 0)2(211lGJI 0)2(3212lGJI 0)(313lGJI 或写成矩阵形式：或写成矩阵形式： 把得到的把得到的M、k矩阵代入振型方程矩阵代入振型方程 000110121012100010001321321lGJI 0)(2AkM注：或设注：或设 2211,AA000)11012101210001

33、0001(3212AAAlGJI设设 ，则上式成为：，则上式成为：2GJIl注：可以直接代入频率方程注：可以直接代入频率方程 02kM频率方程。频率方程。 由一元三次方程求根公式解得：由一元三次方程求根公式解得： 0210121012展开为：展开为： 0165231=0.198，2=1.555，3=3.247代入代入 式中，求得三个固有频率式中，求得三个固有频率 2GJIlIlGJ445. 01IlGJ247. 12IlGJ802. 13 代回振型方程中的前二个（用代回振型方程中的前二个（用），把），把A3移到移到方程右端：方程右端： 将上面求得的将上面求得的1，2，3分别代入分别代入，就得，

34、就得到三个固有振型：到三个固有振型： 取取A1=1（即对标准化），得（即对标准化），得 0)2(21AA321)2(AAA11A 22A1)2(23A即：即： )34( : )2( :1:2321AAA247. 2:802. 1:1:)1(3)1(2)1(1AAA)802. 0( :445. 0:1:)1(3)1(2)1(1AAA555. 0: )247. 1( :1:)1(3)1(2)1(1AAA1802. 1247. 21445. 0802. 01247. 1555. 0编制编制m文件文件.k=2 -1 0;-1 2 -1;0 -1 1;m=1 0 0;0 1 0;0 0 1;vibrat

35、ionmode,eigenvalues=eig(k,m)v1=vibrationmode(:,1)v2=vibrationmode(:,2)v3=vibrationmode(:,3)b1=v1./v3b2=v2./v3b3=v3./v3for i=1:3subplot(3,1,i) plot(0 1 2 3,0 b3(i) b2(i) b1(i)hold onend运行结果运行结果应用应用MATLAB求解求解建立系统的运动微分方程建立系统的运动微分方程 在前一节看到，为了求多自由度系统的固有频在前一节看到，为了求多自由度系统的固有频率和固有振型，必须解高次代数方程（频率方程或率和固有振型，必须

36、解高次代数方程（频率方程或称特征方程），其次数与系统自由度数目一致。当称特征方程），其次数与系统自由度数目一致。当自由度数增加时，解次数很高的代数方程，是很困自由度数增加时，解次数很高的代数方程，是很困难的。所以一般用近似的数值方法。这里介绍的是难的。所以一般用近似的数值方法。这里介绍的是一种数值方法一种数值方法“矩阵迭代法矩阵迭代法”。 多自由度相同自由振动微分方程一般形式：多自由度相同自由振动微分方程一般形式： 0qkqM 设特解为：设特解为： tieAq三、用矩阵迭代法求固有频率和固有振型三、用矩阵迭代法求固有频率和固有振型代入方程得振型方程代入方程得振型方程 用用k-1前乘各项，移项后

37、得到前乘各项，移项后得到 02AkAM121AAMk引入符号引入符号 1MkD动力矩阵，动力矩阵， 就得到迭代格式（振型方程）：就得到迭代格式（振型方程）： 12AAD 对于一个给定的系统，动力矩阵对于一个给定的系统，动力矩阵D是可知的。是可知的。矩阵迭代法就是用逐次迭代、逐次逼近的方法求解矩阵迭代法就是用逐次迭代、逐次逼近的方法求解振型方程，以确定特征值振型方程，以确定特征值 和特征向量（固有振型）和特征向量（固有振型）A。 21集中质量法概念m/4m/4m/4m/3m/3m/2l/4l/4l/4l/4l/3l/2l/3l/3l/2例例:集中质量法求简支梁的模态集中质量法求简支梁的模态力学模

38、型:有一简支梁，其上有一个附加质量（位置和尺寸均在图上标注），梁长L=2m ,截面尺寸b*h为0.04m*0.02m ,截面惯性矩为2.67* ，密度为7920kg/，质量块的质量M为3.0kg 弹性模量E=210GPa 离散810一一.力学模型及集中质量法简化力学模型及集中质量法简化质量矩阵0043.168000006.16804003.168004mmMmM3591171.7552.1451.36511 16112.1453.122.1451076871191.3652.1451.755lEIF柔度矩阵矩阵 = 称作系统的动力矩阵 = DFM二、二、MATLAB计算计算m=1 0 0;0

39、1.94 0;0 0 1;m=1 0 0;0 1.94 0;0 0 1;h=0.00000195h=0.00000195* *3.186;3.186;b=9 11 7;11 16 11;7 11 9;b=9 11 7;11 16 11;7 11 9;f=h.f=h.* *b;b;a=fa=f* *m;m;d=1 1 1;d=1 1 1;for i=1:16;for i=1:16;y=ay=a* *d;d;d=(1/y(3,1)d=(1/y(3,1)* *y;y;endendv=dv=ds=ms=m* *y;y;w1=1/sqrt(s(3,1)subplot(3,1,1)plot(0 2,0 0

40、)hold onplot(0 0.5 1 1.5 2,0 v(1,1) v(2,1) v(3,1) 0)p=v*m*v;q=v*v*m;j=a-(s(3,1)/p)*q;d=1 1 -1;% x(1)=d;for i=1:16;y=j*d;% l=y(i)d=(1/y(3,1)*y;endv=ds=m*y;w2=1/sqrt(s(3,1)subplot(3,1,2)plot(0 2,0 0)hold onplot(0 0.5 1 1.5 2,0 v(1,1) v(2,1) v(3,1) 0)q=v*v*m;r=j-(s(3,1)/p)*q;d=1 -1 1;% x(1)=d;for i=1:1

41、6;y=r*d;% l=y(i)d=(1/y(3,1)*y;endv=ds=m*y;w3=1/sqrt(s(3,1)subplot(3,1,3)plot(0 2,0 0)hold onplot(0 0.5 1 1.5 2,0 v(1,1) v(2,1) v(3,1) 0)运行结果四、四、 强迫振动及减振器强迫振动及减振器m=50 10;k=200 40;c=10 6;N=m(2) c(2) k(2);c(2) k(2);D=m(1)*m(2) (c(1)+c(2)*m(2)+c(2)*m(1) . (k(1)+k(2)*m(2)+k(2)*m(1)+c(1)*c(2) . k(1)*c(2)+

42、c(1)*k(2) k(1)*k(2);sy=tf(N,D);% y,t=impulse(transferab(m,k,c),20)y,t=impulse(sy,20);subplot(2,1,1);plot(t,y(:,1,1);ylabel(x_1(t);title(impulse response of m_1)subplot(2,1,2);plot(t,y(:,2,1);xlabel(time t);ylabel(x_2(t)title(impulse response of m_2);编制编制MATLAB文件文件运运行行结结果果减振器：求减振器：求m1上施加一上施加一外力时，物体外力

43、时，物体m1和和m2位移的频率响应函数。位移的频率响应函数。程序如上程序如上k=200 40;c=10 6;omega=0.0:0.005:4for i=1:2 if i=1 sys=tf(1,m(1) c(1) k(1); mag,phas=bode(sys,omega); plot(omega,mag(1,:),-); hold on; elseN=m(2) c(2) k(2);c(2) k(2);D=m(1)*m(2) (c(1)+c(2)*m(2)+c(2)*m(1) . (k(1)+k(2)*m(2)+k(2)*m(1)+c(1)*c(2) . k(1)*c(2)+c(1)*k(2) k(1)*k(2);sys=tf(N,D); sys=tf(N,D); mag,phas=bode(sys,omega); plot