机械模态分析

1 题目 完成一个综合作业 What I hear I forgot What I see I remember What I do I understand 作业 如图所示的两自由振动系统 已知 m1 100kg m2 3 5kg k1 10000N m k2 600N m c2 1 2N m 1 s F1 t F1ej t 求 1 物理坐标下的振动微分方程 2 频响函数矩阵 3 频响函数的模态展式矩阵 4 脉冲相应函数 5 画出 H11 的幅频特性曲线 相频特性曲线 实频特性曲线 虚频特性曲线 Nyquist 图 Bode 图 6 固有频率 阻尼固有频率 7 画出振型图 8 模态坐标系下的振动微分方程 9 模态参数 复模态质量 复模态刚度 复模态阻尼 解 1 振动微分方程 振动微分方程 对质量 m1 m2绘分离体图 如图 1 1 用牛二定律列分离体在铅垂方向的力平衡方程得 1 1 12212211 111 22122122 Fc xxkxxk xm x c xxkxxm x 将 1 1 整理可得 1 2 1122112211 222222 22 0 00 mxccxkkkxF mcckkx xx 且 m1 100 m2 3 5 k1 10000 k2 600 c2 1 2 代入 1 2 得 1 3 1111 2 22 10001 21 210600600 03 51 21 26006000 xxxF x xx 可以得出此二自由度系统振动微分方程为 M x C xKxf t 其中 图 1 1 系统的分 1000 03 5 M 1 21 2 1 21 2 C 10600600 600600 K 1 0 F f x 离体图 2 频响函数矩阵 频响函数矩阵 由书 P25 1 4 58 公式可知 此二自由度系统频响函数矩阵为一 2 2 方阵 其表达式为 其中 21 HKMj C 1000 03 5 M 1 21 2 1 21 2 C 10600600 600600 K 2 2 1 写成矩阵形式 2 2 1 2 1112 2 2122 10600 100600 6006003 5 HHjj HHjj 3 频响函数的模态展式矩阵 频响函数的模态展式矩阵 1 求解瑞利阻尼矩阵 由于粘性阻尼矩阵 C 无法进行正交性对角化 故不能直接应用坐标变换将 1 3 解耦 由于在该题 中 粘性阻尼相对很小 对于小阻尼振动系统 可以利用瑞利比例阻尼来代替粘性阻尼 以获得可对角 化的阻尼矩阵 1 瑞利比例阻尼系数的确定 瑞利比例阻尼 其中 为瑞利CMK 1000 03 5 M 10600600 600600 K 比例阻尼系数 瑞利比例阻尼系数存在以下关系 其中为圆频率 为系统固有频率 书中表示为 为阻尼 1 1 1 2 2 2 22 22 i 2 ii f i f 0i i 比 0 i i i 将上式写为矩阵形式 1 11 22 2 1 22 1 22 可得 其中 3 1 1 1 11 22 2 1 22 1 22 2 ii f 0 i i i 由此可知 只要我们确定了一个系统任意两阶的固有频率及其阻尼比 就可以确定出瑞利比例阻尼 系数 从而得到瑞利比例阻尼矩阵 2 求该二阶系统的一 二阶固有频率及其阻尼比 利用求解该系统振动微分方程的特征值来确定固有频率及其阻尼比 由书 M x C xKxf t i P23 1 4 43 1 4 46 公式为求解步骤 下面利用 Matlab 来计算固有频率和阻尼比 0i i 编写 Matlab 程序 polynomial m 求特征方程 程序如下 syms x m1 100 m2 4 k1 10000 k2 600 c2 1 2 M m1 0 0 m2 3 C c2 c2 c2 c2 K k1 k2 k2 k2 k2 y det M x 2 C x K 解以上求得的多项式 p 350 124 2 97100 10000 6000000 x0 roots p 由特征值可得 22 011 0 02399 63759 6375 1 1 01 0 0239 0 0025 9 6375 22 022 0 201413 584013 5855 2 2 02 0 0694 0 0062 11 1791 3 求瑞利比例阻尼系数及瑞利比例阻尼矩阵 根据公式 3 1 编写 Matlab 程序 rayleigh m 求解特征方程 程序如下 function Cr rayleigh 计算瑞利阻尼系数alpha和beta 4 xi1 0 0025 xi2 0 0062 f1 9 6375 f2 13 5855 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 3 5 k1 10000 k2 600 M m1 0 0 m2 K k1 k2 k2 k2 k2 Cr alpha M beta K 可知 瑞利比例阻尼系数 0 4628 4 2 0878 10 瑞利比例阻尼矩阵 44 0660 0 1253 0 1253 1 4945 C 2 求解模态矩阵 及特征矢量矩阵 书P23已说明根据粘性比例阻尼振动系统的微分方程所求得的特征矢量与该系统无阻尼振动下求得的 特征矢量相等 因此 我们可以利用求此二阶系统在无阻尼振动下的微分方程的特征矢量更简单的得出 模态矩阵 改写Matlab程序polynomial m求解此二阶系统在无阻尼振动下的微分方程的特征方程 程序如下 syms x m1 100 m2 3 5 k1 10000 k2 600 M m1 0 0 m2 K k1 k2 k2 k2 k2 y det K x 2 M 解以上求得的多项式 p 350 97100 6000000 x0 roots p 5 可知 将其分别代入回 可得 2 01 92 9017 2 02 184 5269 2 0KM 3 2 11 21 1309 83600 0 600274 84 12 22 7852 69600 0 60045 84415 求得模态矩阵 1112 12 2122 11 2 1813 088 3 求解频响函数的模态展式矩阵 1 求模态质量矩阵 模态刚度矩阵和模态阻尼矩阵 1263340 1386 113 08803 52 1813 0880 1386699 5351 T i diag mM 12 181060060011108400300 113 0886006002 1813 088300129080 T i diag kK 12 18 44 0660 0 125311 51 1685 3 3382 113 088 0 1253 1 49452 1813 088 0 4879 300 0675 T i diag cC 2 由此可得频响函数的模态展式为 3 3 2 2 1 T ii i iii H kmj c 写成矩阵形式为 6 22 111221112212 2222 1112221112221112 22 2122 211122122122 2222 111222111222 kmj ckmj ckmj ckmj cHH H HH kmj ckmj ckmj ckmj c 将所求 代入 i diag m i diag k i diag c 11122121 2122 12 18113 088 108400 129 0851 1685 129080 108 4300 0675 2 184 752413 088 171 296 HH Hjj HH 4 脉冲响应函数 脉冲响应函数 对 3 3 作傅立叶逆变换 得到脉冲响应函数矩阵 4 1 2 1 sin i T t ii di i idi h tet m 5 的幅频 相频 实频 虚频特性曲线以及导纳图和博德图的幅频 相频 实频 虚频特性曲线以及导纳图和博德图 11 H 1 的幅频特性曲线 与的关系 11 H 11 H 其中 为阻尼比 代入可 11 22222222 11112222 11 1 4 1 4 H kk 01 i 得 11 2222 2222 11 108400 1 4 0 0025 129080 1 4 9 6375 92 901792 9017184 5269184 5296 H 5 1 编写 Matlab 程序 figure1 m 画图 程序如下 omega 5 01 15 o1 omega sqrt 92 9017 o2 omega sqrt 184 5296 k1 108400 k2 129080 xi1 0 0025 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 on xlabel 频率 Hz ylabel 幅值 mm title m1的一阶幅频特性 输出图形 7 2 的相频特性曲线 与的关系 11 H 11 代入可得 1122 11 22 12 22 arctanarctan 11 5 2 1122 2 0 00252 0 0062 9 637513 584 arctanarctan 11 92 9017184 5269 编写 Matlab 程序 figure2 m 画图 程序如下 omega 5 01 15 o1 omega sqrt 92 9017 o2 omega sqrt 184 5269 k1 108400 k2 129080 xi1 0 0025 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 on xlabel 频率 Hz ylabel 相位角 title m1的一阶相频特性 输出图形 8 3 的实频特性曲线 与的关系 11 H 11 R H 代入可得 22 12 11 22 222222 11112222 11 1414 R H kk 5 3 22 11 22 2222 22 11 92 9017184 5269 10840014 0 0025 12908014 0 0062 92 901780184 526913 584 R H 编写 Matlab 程序 figure3 m 画图 程序如下 omega 5 01 15 o1 omega sqrt 92 9017 o2 omega sqrt 184 5269 k1 108400 k2 129080 xi1 0 0025 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 plot omega y11 LineWidth 2 grid on xlabel 频率 Hz ylabel 幅值 mm title m1的一阶实频特性 输出图形 9 4 的虚频特性曲线 与的关系 11 H 11 I H 代入可得 1122 11 22 222222 11112222 22 1414 I H kk 11 22 2222 22 2 0 0252 0 0062 9 637513 584 10840014 0 0025 12908014 0 0062 92 901792 9017184 5296184 5296 I H 5 4 编写 Matlab 程序 figure4 m 画图 程序如下 omega 5 01 15 o1 omega sqrt 92 9017 o2 omega sqrt 184 5269 k1 108400 k2 129080 xi1 0 0025 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 plot omega y11 LineWidth 2 axis 5 20 0 0035 0 005 grid on xlabel 频率 Hz ylabel 幅值 mm title m1的一阶虚频特性 输出图形 10 5 的导纳图 与的关系 11 H 11 R H 11 I H 编写 Matlab 程序 figure5 m 画图 程序如下 omega 0 01 15 o1 omega sqrt 92 9017 o2 omega sqrt 184 5269 k1 108400 k2 129080 xi1 0 0025 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 2 xi2 o2 k2 1 o2 2 2 2 xi2 o2 2 plot yR yI LineWidth 2 axis square grid on xlabel 实频幅值 mm ylabel 虚频幅值 mm title m1的一阶Nyquist图 输出图形 11 6 的博德图 与的关系 11 H 11 lg H lg 编写 Matlab 程序 figure6 m 画图 程序如下 omega 5 01 15 o1 omega sqrt 92 9017 o2 omega sqrt 184 5269 k1 108400 k2 129080 xi1 0 0025 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 grid on xlabel 频率取对数 ylabel 幅值取对数 title m1的一阶Bode图 输出图形 12 6 固有频率和阻尼固有频率 固有频率和阻尼固有频率 1 固有频率 特征值的模 i Hz 22 0111 0 02399 63759 6375 Hz 22 0222 0 201413 58413 585 2 阻尼固有频率 特征值的虚部 i Hz Hz 1 9 6375 d 2 13 585 d 7 振型图 振型图 根据 3 2 求得的模态矩阵 可以画出系统一 二阶 1112 12 2122 11 2 1813 088 主振型图为 13 一阶主振型图 二阶主振