机械模态分析作业:如图1所示是一个单自由系统附件一个减振器形成的的两自由振动系统,已知m 1=105kg ，m 2=7kg ，k 1=10000N/m ，k 2=410N/m ，c 2=1.15N ·m-1·s ，F 1(t)=F 1e j ωt 。

求：（简化为粘性比例阻尼进行实模态分析）1． 物理坐标下的振动微分方程； 2． 频响函数矩阵；3． 频响函数的模态展式矩阵； 4． 脉冲相应函数；5． 画出H 11（ω）的幅频特性曲线，相频特性曲线，实频特性曲线，虚频特性曲线，Nyquist 图，Bode 图； 6． 固有频率，阻尼固有频率； 7． 画出振型图；8． 模态坐标系下的振动微分方程；9． 模态参数：复模态质量，复模态刚度，复模态阻尼。

10.按实模态系统，给出灵敏度分析。

11.集全班同学的数据（必要的话再补做不同m 2，k 2，c 2参数下的数据，画出x1的最大振幅与m 2，k 2，c 2，的变化曲线，从而分析出减振器的最佳参数。

解：1．振动微分方程对质量m 1、m 2绘分离体图（如图1-1），用牛二定律列分离体在铅垂方向的力平衡方程得1221221111122122122()()()()F c x x k x x k x m x c x x k x x m x ∙∙∙∙∙∙∙∙+-+--=----= （1.1）将（1.1）整理可得：112211221122222222000m x c c x k k k x F m c c k k x x x ∙∙∙∙∙∙⎡⎤⎡⎤-+-⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥++=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦（1.2）且m 1=105、m 2=7、k 1=10000、k 2=410、c 2=1.15，代入（1.2）得： ⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎦⎤⎢⎣⎡+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎦⎤⎢⎣⎡∙∙∙∙∙∙0 410 410-410- 104101.15 1.15- 1.15- 15.17 0 0 1051212121F x x x x x x （1.3）可以得出此二自由度系统振动微分方程为：()M x C x Kx f t ∙∙∙++= 其中M=⎥⎦⎤⎢⎣⎡7 0 0 105；C=⎥⎦⎤⎢⎣⎡ 1.15 1.15- 1.15- 15.1；K=⎥⎦⎤⎢⎣⎡410 410- 410 - 10410；f(t)=⎥⎦⎤⎢⎣⎡0 1F 图1-1、系统的分离体图 2．频响函数矩阵由书P25（1.4-58）公式可知，此二自由度系统频响函数矩阵为一2×2方阵，其表达式为：图1 两自由度振动系统21()()H K M j C ωωω-=-+，其中M=⎥⎦⎤⎢⎣⎡7 0 0 105；C=⎥⎦⎤⎢⎣⎡ 1.15 1.15- 1.15-15.1；K=⎥⎦⎤⎢⎣⎡410 410- 410 - 10410； （2.1）写成矩阵形式：（2.2）3．频响函数的模态展式矩阵1）求解瑞利阻尼矩阵由于粘性阻尼矩阵C 无法进行正交性对角化，故不能直接应用坐标变换将（1.3）解耦。

由于在该题中，粘性阻尼相对很小，对于小阻尼振动系统，可以利用瑞利比例阻尼来代替粘性阻尼，以获得可对角化的阻尼矩阵。

（1）瑞利比例阻尼系数的确定瑞利比例阻尼：K +M =βαC ，其中⎥⎦⎤⎢⎣⎡410 410- 410 - 10410；a 、β 为瑞利比例阻尼系数瑞利比例阻尼系数存在以下关系：1112222222βωαξωβωαξω⎧+=⎪⎪⎨⎪+=⎪⎩，其中i ω为圆频率i i f πω2=（i f 为系统固有频率，书中表示为i 0ω）；i ξ为阻尼比i ii 0=ωσξ将上式写为矩阵形式：111222122122ωωξαξωβω⎡⎤⎢⎥⎡⎤⎡⎤⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎣⎦可得：1111222122122ωωξαξβωω-⎡⎤⎢⎥⎡⎤⎡⎤⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎣⎦，其中、i i f πω2=，i ii 0=ωσξ （3.1）由此可知，只要我们确定了一个系统任意两阶的固有频率及其阻尼比，就可以确定出瑞利比例阻尼系数，从而得到瑞利比例阻尼矩阵。

（2）求该二阶系统的一、二阶固有频率及其阻尼比利用求解该系统振动微分方程()M x C x Kx f t ∙∙∙++=的特征值i λ来确定固有频率及其阻尼比。

由书P23（1.4-43）-（1.4-46）公式为求解步骤，下面利用Matlab 来计算固有频率i 0ω和阻尼比i ξ： 编写Matlab 程序polynomial.m 求特征方程，程序如下：syms x;m 1=105; m2=7; k1=10000; k2=410; c2=1.15; M=[m1 0;0 m2];C=[c2 -c2;-c2 c2];K=[k1+k2 -k2;-k2 k2];y=det(M*x^2+C*x+K)解以上求得的多项式：>> p=[735 644 115920 11500 4100000]; >> x 0=roots(p)由特征值可得：2650.72593.72872.02211≈+==0λω、0396.02593.72872.011===01ωσζ 2882.102626.107253.022≈+==202λω、0707.02626.107253.0222===0ωσζ（3）求瑞利比例阻尼系数及瑞利比例阻尼矩阵根据公式（3.1）编写Matlab 程序rayleigh.m 求解特征方程，程序如下：function Cr=rayleigh()%--计算瑞利阻尼系数alpha和beta--xi1=0.0396; xi2=0.0707; f1=7.2650; f2=10.2882;omega1=2*pi*f1;omega2=2*pi*f2;A=[1/(2*omega1) omega1/2;1/(2*omega2) omega2/2];xi=[xi1;xi2];x=inv(A)*xi;alpha=x(1,1)beta=x(2,1)%--计算瑞利阻尼矩阵Cr(2*2)--alpha=-1.8801;beta=0.0026;m1=105; m2=7; k1=10000; k2=410;M=[m1 0;0 m2];K=[k1+k2 -k2;-k2 k2];Cr=alpha*M+beta*K;可知：瑞利比例阻尼系数-1.8801=α、0026.0=β瑞利比例阻尼矩阵⎥⎦⎤⎢⎣⎡=12.0947-1.0660- 1.0660- 3445.170-C 2）求解模态矩阵（及特征矢量矩阵）书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-x^2*M) 解以上求得的多项式：>> p=[735 -115920 4100000]; >> x0=roots(p)可知：554653201.=ω、1597.104=202ω。

将其分别代入回0)(=-2ϕωM K ，可得：035.1178 410- 410- 767.47862111=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡ϕϕ、0319.1179-410- 410-7685.526-2212=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡ϕϕ求得模态矩阵[]⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡==11112 1 1.2848 - 11.67501 122 22 ϕϕϕϕϕϕϕ3）求解频响函数的模态展式矩阵（1）求模态质量矩阵、模态刚度矩阵和模态阻尼矩阵[]⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡==116.5550 0 0 1394.10591.2848- 11.67501 1 7 00 1051.2848-1 11.6750 1 ϕϕM m diag T i[]⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡==121340 0 0 6721.806251.2848- 11.67501 1 410 410-410- 104101.2848- 1 11.67501 ϕϕK k diag Ti[]⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-⎥⎦⎤⎢⎣⎡==187.6- 0 0 1843.8-1.2848- 11.67501 1 12.0947- 1.0660- 1.0660- 3445.1701.2848- 1 11.6750 1 ϕϕC c diag T i（2）由此可得频响函数的模态展式为：221()Ti i i i i i H k m j c ϕϕωωω==-+∑ （3.3）写成矩阵形式为：22111221112212222211122211122211122221222111221221222222111222111222()()()()()k m j c k m j ck m j c k m j c H H H H H k m j c k m j c k m j c k m j c ϕϕϕϕϕϕωωωωωωωωωωωωωϕϕϕϕϕϕωωωωωωωω⎡⎤++⎢⎥-+-+-+-+⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦++⎢⎥-+-+-+-+⎣⎦将所求ϕ、[]i diag m 、[]i diag k 、[]i diag c 代入：⎥⎦⎤⎢⎣⎡⨯-+⎥⎦⎤⎢⎣⎡⨯-=⎥⎦⎤⎢⎣⎡=-2-2 1.6507 2848.1- 1.2848- 1)6.1875550.116-121340( 305613667501167501118.18431059.1394-8062.56721)( )()( )()(1122211211ωωωωωωωωωj . .. j H H H H H ）（4．脉冲响应函数对（3.3）作傅立叶逆变换，得到脉冲响应函数矩阵：21()sin iT ti i di i i di h t e tm σϕϕωω-==∑ （4.1）5．11()H ω的幅频、相频、实频、虚频特性曲线以及导纳图和博德图 1）11()H ω的幅频特性曲线：11()H ω与ω的关系11222222221111222211()(1)4(1)4H k k ωξξ=+-Ω+Ω-Ω+Ω，其中01ωωΩ=，i ξ为阻尼比。

代入可得：8471)x(0.07074)1(1213401)x(0.03964)7802.521(8062.567211)(22211105.⨯+105.8471-+52.7802⨯+-=222ωωωωωH （5.1）编写Matlab 程序figure1.m 画图，程序如下：omega=5:.01:15;o1=omega/sqrt(52.7802); o2=omega/sqrt(105.8471); k1=56721.8062; k2=121340; xi1=0.0396; xi2=0.0707;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）11()H ω的相频特性曲线：11ϕ与ω的关系112211221222arctan arctan 11ξξϕ⎛⎫⎛⎫-Ω-Ω=+ ⎪ ⎪-Ω-Ω⎝⎭⎝⎭。