机械模态分析作业例题

题目：完成一个综合作业(What I hear, I forgot. What I see, I remember. What I do, I understand.)作业:如图所示的两自由振动系统,已知m 1=100kg ，m 2=5kg ，k 1=10000N/m ，k 2=500N/m ，c 2=1N ·m-1·s ，F 1(t)=F 1e j ωt。

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

解：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 cc x k k k x F m c c k k x x x ∙∙∙∙∙∙⎡⎤⎡⎤-+-⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥++=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦（1.2） 且m 1=100、m 2=5、k 1=10000、k 2=500、c 2=1，代入（1.2）得：111122************50005115005000x x x F x x x ∙∙∙∙∙∙⎡⎤⎡⎤--⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥++=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦（1.3） 可以得出此二自由度系统振动微分方程为：()M x C x K x f t ∙∙∙++=其中100005M ⎡⎤=⎢⎥⎣⎦；1111C -⎡⎤=⎢⎥-⎣⎦；10500500500500K -⎡⎤=⎢⎥-⎣⎦；1()0F f x ⎡⎤=⎢⎥⎣⎦图1-1、系统的分离体图2．频响函数矩阵由书P25（1.4-58）公式可知，此二自由度系统频响函数矩阵为一2×2方阵，其表达式为：21()()H K M j C ωωω-=-+，其中100005M ⎡⎤=⎢⎥⎣⎦；1111C -⎡⎤=⎢⎥-⎣⎦；10500500500500K -⎡⎤=⎢⎥-⎣⎦（2.1） 写成矩阵形式：12111222122()()10500100500()()()5005005H H j j H H H j j ωωωωωωωωωωω-⎡⎤-+--⎡⎤==⎢⎥⎢⎥---+⎣⎦⎣⎦（2.2）3．频响函数的模态展式矩阵1）求解瑞利阻尼矩阵由于粘性阻尼矩阵C 无法进行正交性对角化，故不能直接应用坐标变换将（1.3）解耦。

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

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

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

由书P23（1.4-43）-（1.4-46）公式为求解步骤，下面利用Matlab 来计算固有频率0i ω和阻尼比i ξ：编写Matlab 程序polynomial.m 求特征方程，程序如下：syms x;m 1=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*x^2+C*x+K)解以上求得的多项式：>> p=[500 105 102500 100005000000];>> x 0=roots(p)由特征值可得：20118.94528.9453ωλ==≈、11010.03560.00408.9453σξω===02211.1791ωλ==≈、22020.06940.006211.1791σξω===（3）求瑞利比例阻尼系数及瑞利比例阻尼矩阵 根据公式（3.1）编写Matlab 程序rayleigh.m 求解特征方程，程序如下：function Cr=rayleigh() %--计算瑞利阻尼系数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];x=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;可知：瑞利比例阻尼系数-0.3004α=、4=2.374110β-⨯瑞利比例阻尼矩阵-27.5422-0.1187-0.1187-1.3830C ⎡⎤=⎢⎥⎣⎦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=[500 -102500 5000000]; >> x0=roots(p)可知：20180ω=、202125ω=。

将其分别代入回2()0K M ωϕ-=，可得：112125005000500100ϕϕ-⎡⎤⎡⎤=⎢⎥⎢⎥-⎣⎦⎣⎦ 、122220005000500125ϕϕ--⎡⎤⎡⎤=⎢⎥⎢⎥--⎣⎦⎣⎦（3.2）求得模态矩阵[]11121221221154ϕϕϕϕϕϕϕ⎡⎤⎡⎤===⎢⎥⎢⎥-⎣⎦⎣⎦ 3）求解频响函数的模态展式矩阵（1）求模态质量矩阵、模态刚度矩阵和模态阻尼矩阵151000112250[]1405540180Ti diag m M ϕϕ⎡⎤⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦⎣⎦ 151********1180000[]1450050054022500Ti diag k K ϕϕ-⎡⎤⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦⎣⎦15-27.5422-0.11871163.30420.000963.3042[]14-0.1187-1.3830540.000948.7206048.7206Ti diag c C ϕϕ---⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤===≈⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-----⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦（2）由此可得频响函数的模态展式为：221()Ti ii i i iH k m j c ϕϕωωω==-+∑（3.3）写成矩阵形式为：22111221112212222211122211122211122221222111221221222222111222111222()()()()()k m j c k m j c k 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 代入：111221212122()()1514()(1800022563.3042)(2250018048.7206)()()525416H H H j j H H ωωωωωωωωω---⎡⎤⎡⎤⎡⎤==--⨯+--⨯⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦4．脉冲响应函数对（3.3）作傅立叶逆变换，得到脉冲响应函数矩阵：21()sin i Tti idi i i dih t et m σϕϕωω-==∑（4.1）5．11()H ω的幅频、相频、实频、虚频特性曲线以及导纳图和博德图 1）11()H ω的幅频特性曲线：11()H ω与ω的关系11()H ω=，其中01ωωΩ=，i ξ为阻尼比。

代入可得：1111()H ω=+（5.1）编写Matlab 程序figure1.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./(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ξξϕ⎛⎫⎛⎫-Ω-Ω=+⎪ ⎪-Ω-Ω⎝⎭⎝⎭。