机械模态分析作业例题.doc

1、1题目：完成一个综合作业(What I hear, I forgot. What I see, I remember. What I do, I understand.)作业:如图所示的两自由振动系统,已知m1=100kg，m 2=5kg，k 1=10000N/m，k 2=500N/m，c 2=1Nm-1s，F 1(t)=F1ejt 。求：1 物理坐标下的振动微分方程； 2 频响函数矩阵；3 频响函数的模态展式矩阵；4 脉冲相应函数；5 画出 H11（）的幅频特性曲线，相频特性曲线，实频特性曲线，虚频特性曲线，Nyquist 图，Bode 图；6 固有频率，阻尼固有频率；7 画出振型图；8 模

2、态坐标系下的振动微分方程；9 模态参数：复模态质量，复模态刚度，复模态阻尼。解：1振动微分方程对质量 m1、m 2绘分离体图（如图 1-1） ，用牛二定律列分离体在铅垂方向的力平衡方程得（1.1）1211122()()Fcxkxkmx 将（1.1）整理可得：（1.2）1121212 220 0mxcxkxF且 m1=100、 m2=5、k 1=10000、k 2=500、c 2=1，代入（1.2）得：12220055 0xxx（1.3）可以得出此二自由度系统振动微分方程为： ()MxCKft其中 ； ； ； 图 1-1、系统的分离体图105M1C105K10Fx2频响函数矩阵由书 P25（1.

3、4-58）公式可知，此二自由度系统频响函数矩阵为一 22方阵，其表达式为：，其中 ； ； 21()HKjC05M1C05K（2.1）2写成矩阵形式：（2.2）12112 22()()05150()Hjj 3频响函数的模态展式矩阵1）求解瑞利阻尼矩阵由于粘性阻尼矩阵 C 无法进行正交性对角化，故不能直接应用坐标变换将（1.3）解耦。由于在该题中，粘性阻尼相对很小，对于小阻尼振动系统，可以利用瑞利比例阻尼来代替粘性阻尼，以获得可对角化的阻尼矩阵。（1）瑞利比例阻尼系数的确定瑞利比例阻尼： ，其中 ； ； 、 为瑞利比MK105105K例阻尼系数瑞利比例阻尼系数存在以下关系：，其中 为圆频率 （ 为

4、系统固有频率，书中表示为 ） ； 为阻尼12i 2iifi 0ii比 0ii将上式写为矩阵形式：1122可得： ，其中1122、 （3.1）2iif0ii由此可知，只要我们确定了一个系统任意两阶的固有频率及其阻尼比，就可以确定出瑞利比例阻尼系数，从而得到瑞利比例阻尼矩阵。（2）求该二阶系统的一、二阶固有频率及其阻尼比利用求解该系统振动微分方程的特征值 来确定固有频率及其阻尼比。由书 P23（1.4-43）-（1.4-46 ）公式为求()MxCKfti3解步骤，下面利用 Matlab 来计算固有频率 和阻尼比 ：0ii编写 Matlab 程序 polynomial.m 求特征方程，程序如下：sy

5、ms x;m1=100; m2=5; k1=10000; k2=500; c2=1;M=m1 0;0 m2;C=c2 -c2;-c2 c2;K=k1+k2 -k2;-k2 k2;y=det(M*x2+C*x+K)解以上求得的多项式： p=500 105 102500 10000 5000000; x0=roots(p)由特征值可得： 、2201.3568.94.5310.3560.489、2202.694.78.1020.1.（3）求瑞利比例阻尼系数及瑞利比例阻尼矩阵根据公式（3.1）编写 Matlab 程序 rayleigh.m求解特征方程，程序如下：function Cr=rayleigh

6、()%-计算瑞利阻尼系数alpha和beta-xi1=0.0040; xi2=0.0062; f1=8.9453; f2=11.1791;omega1=2*pi*f1;omega2=2*pi*f2;A=1/(2*omega1) omega1/2;1/(2*omega2) omega2/2;xi=xi1;xi2;4x=inv(A)*xi;alpha=x(1,1)beta=x(2,1)%-计算瑞利阻尼矩阵Cr(2*2)-m1=100; m2=5; k1=10000; k2=500;M=m1 0;0 m2;K=k1+k2 -k2;-k2 k2;Cr=alpha*M+beta*K;可知：瑞利比例阻尼系

7、数 、-0.344=2.710瑞利比例阻尼矩阵 2758-.1-.C2）求解模态矩阵（及特征矢量矩阵）书P23已说明根据粘性比例阻尼振动系统的微分方程所求得的特征矢量与该系统无阻尼振动下求得的特征矢量相等。因此，我们可以利用求此二阶系统在无阻尼振动下的微分方程的特征矢量更简单的得出模态矩阵改写Matlab 程序polynomial.m求解此二阶系统在无阻尼振动下的微分方程的特征方程，程序如下：syms x;m1=100; m2=5; k1=10000; k2=500;M=m1 0;0 m2;K=k1+k2 -k2;-k2 k2;y=det(K-x2*M)解以上求得的多项式： p=500 -10

8、2500 5000000; x0=roots(p)可知： 、 。将其分别代入回 ，可得：201820152()0KM5、 （3.2）122500 1250求得模态矩阵 1122543）求解频响函数的模态展式矩阵（1）求模态质量矩阵、模态刚度矩阵和模态阻尼矩阵 10125045418TidagmM510025TikK15-27.4-.87163.4.963.0403878.726TidagcC （2）由此可得频响函数的模态展式为：（3.3）221()Tii iHkmjc写成矩阵形式为： 2211212121122 222211211()()()kjckjckjckmjcHm 将所求 、 、 、

9、代入：idagikidagc11221 212() 5 4()(80563.04)(0184.706)2Hj j 4脉冲响应函数对（3.3）作傅立叶逆变换，得到脉冲响应函数矩阵：（4.1）21()sinTtidiidhtem5 的幅频、相频、实频、虚频特性曲线以及导纳图和博德图1()H1） 的幅频特性曲线： 与 的关系1()H6，其中 ， 为阻尼比。代入可1222211121()()4()4Hkk01i得：1 2222()80(1)4(0.)50(1)4(0.6)8515（5.1）编写 Matlab 程序 figure1.m 画图，程序如下：omega=5:.01:15;o1=omega/sq

10、rt(80); o2=omega/sqrt(125);k1=18000; k2=22500; xi1=0.0040; xi2=0.0062;y11=1./(k1.*sqrt(1-o1.2).2+(2.*xi1.*o1).2)+1./(k2.*sqrt(1-o2.2).2+(2.*xi2.*o2).2);plot(omega,y11,LineWidth,2);grid onxlabel(频率 Hz)ylabel(幅值 mm)title(m1的一阶幅频特性)输出图形：2） 的相频特性曲线： 与 的关系1()H1。代入可得：1212arctnarctn7（5.2）1 2 20.40.68.9531.

11、79arctnarctn15编写 Matlab 程序 figure2.m 画图，程序如下：omega=5:.01:15;o1=omega/sqrt(80); o2=omega/sqrt(125);k1=18000; k2=22500; xi1=0.0040; xi2=0.0062;y11=atan(-2.*xi1.*o1)./(1-o1.2)+atan(-2.*xi2.*o2)./(1-o2.2);plot(omega,y11,LineWidth,2);axis(5 15 -2 2)grid onxlabel(频率 Hz)ylabel(相位角 )title(m1的一阶相频特性)输出图形：3）

12、的实频特性曲线： 与 的关系1()H1()RH。代入可得：2 2111244Rkk8（5.3）2 21 2 211805()804(.)504(0.6)815RH 编写 Matlab 程序 figure3.m 画图，程序如下：omega=5:.01:15;o1=omega/sqrt(80); o2=omega/sqrt(125);k1=18000; k2=22500; xi1=0.0040; xi2=0.0062;y11=(1-o1.2)./(k1.*(1-o1.2).2+(2.*xi1.*o1).2)+(1-o2.2)./(k2.*(1-o2.2).2+(2.*xi2.*o2).2);plo

13、t(omega,y11,LineWidth,2);grid onxlabel(频率 Hz)ylabel(幅值 mm)title(m1的一阶实频特性)输出图形：4） 的虚频特性曲线： 与 的关系1()H1()IH。代入可得：1 222 21144Ikk91 2 22 20.40.68.9531.79()81()54()0 5IH （5.4）编写 Matlab 程序 figure4.m 画图，程序如下：omega=5:.01:15;o1=omega/sqrt(80); o2=omega/sqrt(125);k1=18000; k2=22500; xi1=0.0040; xi2=0.0062;y11

14、=(1-o1.2)./(k1.*(1-o1.2).2+(2.*xi1.*o1).2)+(1-o2.2)./(k2.*(1-o2.2).2+(2.*xi2.*o2).2);plot(omega,y11,LineWidth,2);axis (5 15 -0.0035 0.0005)grid onxlabel(频率 Hz)ylabel(幅值 mm)title(m1的一阶虚频特性)输出图形：5） 的导纳图： 与 的关系1()H1()R1()IH编写 Matlab 程序 figure5.m 画图，程序如下：omega=0:.01:15;o1=omega/sqrt(80); o2=omega/sqrt(1

15、25);k1=18000; k2=22500; xi1=0.0040; xi2=0.0062;yR=(1-o1.2)./(k1.*(1-o1.2).2+(2.*xi1.*o1).2)+(1-o2.2)./(k2.*(1-o2.2).2+(2.*xi2.*o2).2);yI=(-2.*xi1.*o1)./(k1.*(1-o1.2).2+(2.*xi1.*o1).2)+(-102.*xi2.*o2)./(k2.*(1-o2.2).2+(2.*xi2.*o2).2);plot(yR,yI,LineWidth,2);axis squaregrid onxlabel(实频幅值 mm)ylabel(虚频幅

16、值 mm)title(m1的一阶Nyquist图)输出图形：6） 的博德图： 与 的关系1()H1lg()Hl编写 Matlab 程序 figure6.m 画图，程序如下：omega=5:.01:15;o1=omega/sqrt(80); o2=omega/sqrt(125);k1=18000; k2=22500; xi1=0.0040; xi2=0.0062;y11=log(1./(k1.*sqrt(1-o1.2).2+(2.*xi1.*o1).2)+1./(k2.*sqrt(1-o2.2).2+(2.*xi2.*o2).2);plot(log(omega),y11,LineWidth,2)

17、;grid onxlabel(频率取对数 )ylabel(幅值取对数 )title(m1的一阶Bode图)输出图形：116固有频率和阻尼固有频率1）固有频率：特征值 的模。iHz*220110.3568.94.53Hz2217192）阻尼固有频率：特征值 的虚部。iHz Hz18.945d2.8d7振型图根据（3.2）求得的模态矩阵 ，可以画出系统一、二阶主振型图1212154为：一阶主振型图 二阶主振型图其中：x 表示系统各点的静平衡位置12v1表示系统在一阶模态下各点振幅比v2表示系统在二阶模态下各点振幅比8模态坐标下振动微分方程由书 P25（1.4-59）公式可以得出此二自由度系统在模态坐标下的振动微分方程为：（8.1）()Tiiidagmycdagkyft其中 ， ， ，25018idag18025ik63.0428.76idagc14T代入可写为一个解耦的方程组：（8.2）111222563.048()8750yyFt9复模态质量、复模态刚度、复模态阻尼由书 P30（1.5-31）公式可知，复模态质量、刚度、阻尼为：，其中 为 i 阶模态矢量。 HmiiiiiiMkKcC12iii（9.1）所以复模态质量、刚度、阻尼与实模态时的一致。可得：1）复模态质量： 12025018midag2）复模态刚度： 1200imkk3）复模态阻尼： 1263.428.706icdagc