基于MATLAB的振动模态分析

来源：世旅网

摘要

振动系统是研究机械振动的运动学和动力学，研究单自由系统的振动有着实际意义，因为工程上有许多问题通过简化，用单自由度系统的振动理论就能得到满意的结果。模态是振动系统的一种固有振动特性，模态一般包含频率、振型、阻尼。

振动系统问题是个比较虚拟的问题，比较抽象的理论分析，对于问题的分析可以实体化建立数学模型，通过MATLAB可以转化成为图像。单自由度频率、阻尼、振型的分析，我们可以建立数学模型，最后通过利用MATLAB编程实现数据图形；多自由度主要研究矩阵的迭代求解，我们在分析抽象的理论的同时根据MATLAB编程实现数据的迭代最后可以得到所要的数据，使我们的计算更加简便。

利用MATLAB编程并验证程序的正确性。通过程序的运行，能快速获得多自由度振动系统的固有频率以及主振型，为设计人员提供了防止系统共振的理论依据，也为初步分析各构件的振动情况以及解耦分析系统响应奠定了基础。关键词：振动系统；单自由度；MATLAB；多自由度

Abstract

Vibration system is to study the kinematics and dynamics of mechanical vibration, the vibration of a single free system has practical significance, because there are many engineering problems by simplifying, using the vibration theory of a single degree of freedom system can be satisfied with the results.

Vibration system problems is a relatively virtual problems, more abstract and theoretical analysis, problem analysis for a mathematical model can be materialized by MATLAB can be converted into images. Single degree of freedom frequency, damping, mode shape analysis, we can create mathematical models, the final program data through the use of MATLAB graphics; many degrees of freedom main matrix iterative solution, our analysis based on abstract theory, while MATLAB programming The last iteration of data can be the desired data, so our calculations easier

Using MATLAB programming and verify the correctness of the program.Through the process of operation, can quickly obtain multiple degrees of freedom vibration system and the main vibration mode natural frequency for the design to prevent resonance provide the theoretical basis for the preliminary analysis of the vibration of each component, and laid the decoupling of system response basis.

Key words:vibrating system; Single Degree of Freedom ;MATLAB; multiple degree of

freedom

辽宁工程技术大学毕业设计（论文）

1 绪论

1.1问题的提出

机械振动是一门既古老又年轻的科学，随着人类科学技术的不断进步振动理论得到不断的发展和完善。机械振动在许多情况下是有害的，人们想方设法避免它：另一方面，人们利用机械振动原理制造了各种机械或仪表来为人类服务。振动机械是20世纪后半期得到迅速发展的一类机械，它是利用振动原理来完成各种工艺过程的机械设备。其中，Mathorks公司推出的MATLAB以其强大的功能和易用性受到越来越多科技工作者的欢迎。它把计算、可视化、程序设计融合到了一个交互的工作环境中，可以实现工程计算、算法研究、建模和仿真、数据分析及可视化、科学和工程绘图、应用程序开发（包括图形用户界面程序设计）等功能。它在美国等发达国家的大学里已经成为一种必须掌握的基本编程语言，而在国外的研究设计单位和工业部门，更是早己成为研究和解决工程计算问题的一种标准软件。在国内也有越来越多的科学技术工作者参加到学习和倡导这种语言的行列中来。应用MATLAB软件对选矿用振动筛的振动特性进行研究，可以充分发挥计算机技术的优势，为选矿用振动筛振动特性研究探索新的途径。

在工程振动中，确定系统固有频率与主振型是非常重要的。固有频率是决定系统振动特性的重要物理量，它既是防止系统共振的依据，又是多自由度系统解耦分析（模态分析）的前提，因此研究某系统振动时，首先要求出系统的固有频率。主振型则为初步分析各构件的振动情况以及解耦分析奠定了基础。对于多自由度振动系统，计算系统固有频率与主振型主要有2种方法[1]：(1)利用特征矩阵方程式与特征方程式求解；(2)矩阵迭代法求解。2种方法各有各的特色。对于低自由度的振动系统，方法一容易、快捷。但是在实际工程中，大多数振动

系统都是自由度较多，用特征矩阵方程式与特征方程式求解系统固有频率与主振型这种传缆的计算方法虽然从原则上可行，但当自由度增加时，惯性、刚度阵的阶数增高，计算量也急剧加大，这显然很不方便。但采用矩阵迭代法，即使是自由度很大的振动系统，计算量也只不过是多进行矩阵迭代而已，而且假设的初始矩阵愈接近实际状况，迭代的次数愈少，相应的计算量也愈少。

王超：基于MATLAB的振动系统编程分析

1.2 国内外研究现状

1.2.1机械振动理论的发展状况及应用现状

振动理论是力学的一个重要组成部分[2]，人类对振动现象的认识有悠久的历史。振动力学的物理基础在17世纪已经奠定，到了18世纪，振动力学已从物理学中独立出来。最主要的成就为线性振动理论的形成，它是与数学中的常微分方程和偏微分方程同步发展的。目前，振动及系统按运动微分方程的形式分为以下两种。

线性振动：描述其运动的方程为线性微分方程，相应的系统称为线性系统。线性振动的一个重要特性是线性叠加原理成立。非线性振动[3]:描述其运动的方程为非绒性微分方程，相应的系统称为非线性系统。对于非线性振动叠加原理不再成立。在实际的振动机械或振动系统中，严格的讲，都是非线性的。但是，建立振动系统的非线性力学模型难度大，求解困难，有些问题甚至无解可求。在实际的工程应用中，很多情况下在误差允许的范围之内用线性的方法解决复杂的近线性问题。线性振动有确定的力学模型一一线性微分方程，可以求得准确的解，能够描述出振动系统的主要特征。由于用线性振动的方法能够解决众多的工程实际问题，线性振动的理论一直倍受关注，并且在理论和实验方面已经得到很大的发展和成熟。特别是多自由度系统的振动的理论，可以说既是振动力学的核心又是应用得最广泛的振动理论。线性振动在当今不仅是作为基础科学的力学的一个重要组成部分，而且正走上向工程科学发展的道路，它在航空、机械、船舶、车辆、建筑、水利等工业技术部门中占有愈来愈重要的地位。线性振动的应用可分为两个方面：一个方面是减少由于振动而造成的危害，目的在于减振甚至于避免有害的振动；另一个方面利用振动，如工业上常采用的振动筛选、振动沉桩、振动输送以及按振动理论设计的测量传感器、地震仪等等就是这方面的典型例子。选矿用振动筛是振动筛选设各中的—种，线性振动理论在选矿用振动筛的设计制造及生产运行中有着广泛的应用，有关这方面的内容将在下一节中详细介绍。线性振动的理论在发展过程中产生了一个重要分支，那就是模态分析理论。在对选矿用振动筛进行分析时，需要通过实验来验证理论的正确性，振动实验则需要用到模态分析技术。模态分析技术从20世纪60年代后期发展至今已趋成熟[4]。它和有限元分析技术一起，已成为结构动力学中两大支柱。模态分析是结构动力学中的～种“逆问题”分析方法，它与传统的“正问题”方法（主要是指有限元方法）不同，是建立在实验（或实测）的基础上，采用实验与理论相结合的方法来处理工程中的振动问题。目前这一技术已发展成为解决工程中振动问题的重要手段，在机械、航空、航天、土木、建筑、造船、化工等

辽宁工程技术大学毕业设计（论文）

工程领域被广泛应用[5]。近十年来，模态分析理论吸取了振动理论、信号处理、信号分析、数据处理、数理统计及自动控制理论中的有关“营养”，结合自身内容的发展，形成了一套独特的理论为模态分析及参数识别技术的发展奠定了理论基础。模态分析的基础理论概念主要包括；机械阻抗、导纳、传递函数（或频响函数）、实模态、复模态等。模态测试技术主要采用同时测量输入及输出的方法，对一个振动系统来说，可以表示成图l-1所示的框图

输入（激励）系统输出（响应）图1.1 模态分析框图 Fig. 1.1 Modal Analysis Diagram

通过测量激励和响应，进行模念分析可以确定系统。自从FFT问世以来，目前广泛采用宽频带激振技术。其中主要有脉冲、阶跃激励，快速正弦扫描等瞬态激励和纯随机、伪随机、周期随机、瞬态随机等激励方法。此外，由于F弦慢扫描技术测试精度高，它仍不失为重要激励手段。模态参数辨识的频域方法有：分量分析法、导纳圆辨识方法、正交多项式曲线拟合、非线性优化辨识方法等。模态参数辨识的时域方法与模态参数辨识的频域方法不同，它无需将所测得的响应与激励的时间历程信号变换到频域中去，而是直接在时域中进行参数辨识。它与频域法相比，两者所采取的分析路线不同，如图1.2所示。

时域信号 FFT 频率信号数估计传递函传递函数参数识别模态参数图1.2 模态参数辨识分析路线框图

Fig.1.2 Modal parameter identification of line diagram

时域法比频域法发展较晚，但近几年来有长足的进展。自70年代以来主要有： Ibrahim时域法（简称LTD法）、最小二乘复指数法(LSCE法)、多参考点复指数法(PRCE法)、特征系统实现算法(ERA)。模态分析技术在动态载荷识别、模型修正与结构动力修改中有广泛的应用，结构动态特征灵敏度分析是非常重要的方法之一。模态综合技术主要有组合系统法和模态综合法。随着电子技术与计算机技术的迅速发展，模态分析已成为解决复杂结构振动问题的主要工具，并与计算机辅助设计(CAD)．计算机辅助实验(CAT)相结合，

王超：基于MATLAB的振动系统编程分析

进入产品设计阶段，作为计算机辅助工程中的重要环节，有着广泛的应用[6]。

1.2.2 MATLAB软件的发展状况及应用现状

MATLAB软件概述：MATLAB的名称源自Matrix Laboratory，是一门计算语言口[7]。在工程计算领域，计算机技术的应用正逐步将科技人员从繁重的计算工作中解放出来。在科学计算和工程应用的过程中，一些技术人员尝试用Basic，Fortran以及C语言编制程序来减轻计算的工作量，但编制程序不仅要掌握所用语-的语法，还要对有关算法进行深入分析。为了满足用户对工程数学计算的要求，MATLAB的功能、特点、应用范围：MATLAB越来越广地被人们应用是源于它在求解方程、数值计算、程序编写上的优点，而它的这些优点是由它的功能和特点决定的。MATLAB的主要功能:(1)数值计算功能，一条MATLAB语句相当于几十条C语言或Fortran语言的语句。(2)符号计算功能，利用MATLAB的符号计算功能可以清晰地获得解的表达式，对于避免出错和提高程序的可读性均有很大的帮助。(3)数据分析和可视化功能，在科学计算和研究工作中，技术人员经常会遇到大量的原始数据，而对数据的分析往往难于入手。MATLAB能将这些数据以图形的方式显示出来，不仅使数据间的关系清晰明了，而且对于揭示其内在本质往往有着非常重要的作用。MATLAB提供了良好的用户界面，许多函数本身会自动绘制出图形，而且会自动选取坐标刻度，绘制出直角坐标、极坐标、对数坐标下的二维和三维图形，以及条形图、直方图、等高线图、饼形图、离散数据图和瀑布图等专用图形。(4)文字处理功能。MATLAB的主要特点： (1)功能强大，MATLAB不但在数值计算和符号计算方面具有强大的功能，而且在计算结果的分析和数据可视化方面也有着其它类似软件难以匹敌的优势[9]。Notebook，Simulink功能以及各种专业工具箱将MATLAB的应用扩展到非常广的领域。(2)界面友好、编程效率高，MATLAB的指令表达方式与标准教科书的数学表达式非常相近，用户不需要有较高的计算机编程基础，只要按照计算要求输入表达式，

MATLAB将为用户计算出结果。同时使用MATLAB语言设计的程序，其编译和执行速度都超过了传统c和Fortran语言设计的程序，在工程计算方面的编程效率也高于其它编程语言。(3)扩展性强，MATLAB的最重要特点之一就是它的可扩展性。这个特点使得用户能够自由地开发自己的应用程序。这些年来，许多使用MATLAB的科学家、工程师和技术人员已经开发出相当多的不同领域的应用程序。MATLAB的应用范围：MATLAB由主包和各种工具箱组成。主包是MATLAB的核心，工具箱是扩展的有专门功能的函数。例如，控制系统工具箱应用于连续和离散系统设计、

辽宁工程技术大学毕业设计（论文）

频域和时域响应等控制领域；信号处理工具箱应用于自适应去噪和压缩、谱分析和估计等信号处理领域；通信工具箱应用于信号编码、调制解调等通信领域。应用MATLAB的各种工具箱可以在很大程度上减小用户编程时的复杂度，因此MATLAB在很广的领域内得到了应用，其典型应用有；自动控制、图像信号处理、生物医学工程、语音处理、雷达工程、信号分析、振动理论、时序分析与建模、化学、统计学、经济学等

1.3 MATLAB语言的优点

MATLAB作为一个以矩阵和数组为核心计算的软件，对矩阵迭代法中的矩阵迭代计算尤其适合[10]。就所查的资料看，以前的学者和研究人员迭代求解系统固有频率与主振型时，大部分都是用Visiual Basic或Fortran语言来编写程序[11]。限于Visiual Basic或Fortran本身语句以及语法的局限性，用这种高级语言编写的程序涉及到选择合适的算法和编写冗长的语言代码以及键入和调试等一系列问题。即使有现成的标准予程序可供调用，要在一些较复杂的、科研问题中编写一个完整的程序仍然是一个复杂的、技巧性很强的工作。因此，用高级语言编写的程序一般代码段较长，需要调用的子程序较多，整个程序的通读性较差。相反，MATLAB则有简洁、可读性强等优点。

1.4本文研究的内容

振动机械在国民经济中占有重要的位置，振动筛是振动机械中的重要一员。一直以来有许多人对振动筛进行设计和研究，但是，振动筛的动态设计和计算机辅助设计近年来刚刚起步。振动特性是振动筛非常重要的有别于非振动机械的一个本质特点，却往往被设计者和制造者简单化。客观的说，一般的振动都是非线性的，但在许多情况下可以近似看作线性来处理。线性振动理论不论从基础理论还是实验技术方面近年来都有很大的发展，特别是应用现代化振动测试仪器测量振动信号以及应用计算机软件来分析处理振动信号，为从事振动研究的科技人员带来了极大的方便。

把振动的理论应用到工程实际中去，切实解决工程中遇到的实际的振动问题是研究振动理论的根本目的。需要对该力学模型进行深入的分析（借助MATLAB软件进行仿真分析）。

本文主要利用MATLAB对振动系统进行模拟分析对于虚拟抽象的理论图像化，处理单自由度振动的3个阻尼和强迫单自由度阻尼振动，多自由度系统振动矩阵迭代求解。

王超：基于MATLAB的振动系统编程分析

2单自由度系统的振动

2.1 单自由度振动系统数学模型的建立

标，由牛顿定律得运动方程为[13]:

[12]

建立和分析有粘性阻尼时的自由度振动微分方程。以静平衡位置为原点建立如图坐

m c kx m x m x

cx

cxkx0 (2-1) mx令

2nc2k,n mm其中n称为衰减系数，单位为1s；n是相应的无阻尼时的固有频率，式(2-1)可以写为：

22nxnxx0 (2-2)

如果进一步令

nn (2-3)

其中无量纲的称为相对阻尼系数，则式(2-2)可写为：

22nxnxx0 (2-4) 为了求解，令

辽宁工程技术大学毕业设计（论文）

xest (2-6)

代入(2-4)后得到特征方程：

2s22nsn0 (2-7)

他的两个特征根为：

s1,2nn21 (2-8)

根据相对阻尼系数的不同大小，可以将阻尼分为三种状态：1时为过阻尼，1时为临界阻尼，01时为欠阻尼。

1）过阻尼状态

1，s1与s2是两个不等的负实根，令

*n21 （2-9）

初始条件

(0)x0 （2-10） x(0)x0,x系统初始条件响应为

x(t)ent(x0ch*t临界阻尼状态

0nx0x*sh*t) （2-11）

1,sn是二重根，方程（2-4）的通解为系统对式（2-10）的初始条件的响应为

0nx0)t] （2-12） x(t)ent[x0(x欠阻尼状态

1，其中

dn12 （2-13）

初始条件响应

x(t)ent(x0cosdt0nx0xdsindt) （2-14）

2.2 参数设定与求解

阻尼比分别取；应用Matlab对式(2-11)和(2-12),（2-14）求解。程序如下：

王超：基于MATLAB的振动系统编程分析

clear,format compact;

a=0.5;t=0:0.1:18;;w0=1;

k=1;x0=1;

wd=w0.*sqrt(1-a*a);x1=wd

y=exp(-a*w0.*t).*(x0.*cos(wd.*t)+((x1+a*wd*x0)./wd)*sin(wd.*t)) figure(1),plot(t,y,'r');hold on a=1.0;t=0:0.1:18; w0=1;wd=1;x1=wd;

y=exp(-wd.*t).*(x0+(x1+wd*x0).*t); figure(1),plot(t,y,'d');hold on

a=2.0;t=0:0.1:18;w0=1;wd=w0*sqrt(a*a-1);

y=exp(-a*w0.*t).*(x0.*cosh(wd.*t)+(x1+a*w0*x0)/w0.*sinh(t));

figure(1),plot(t,y,'v');hold on 结论：

图2-2为Matlab计算后给出的响应曲线，从中可以得到一些重要的结论[14]：在01的情况下，阶跃信号输入时，输出信号为衰减振荡，其振荡角频率（阻尼振荡角频率）为d，幅值按指数衰减越大，阻尼越大，衰减越快。

1时，振荡系统等同于两个一阶系统串联。此时虽然不产生振荡，但也需要经过较长时间才能达到稳态。

在一定的之下，欠阻尼系统能够更快地达到稳态值；而过阻尼系统反应迟饨，动作缓慢，所以系统通常设计成欠阻尼系统，取值为2

辽宁工程技术大学毕业设计（论文）

x(t)1.41.210.80.6a=1.00.4a=2.00.20a=0.5-0.2024681012141618t(s)

图2-2

算例绘制无阻尼单自由度系统的固有频率和周期随静变形的变化曲线。固有频率n和周期n

ngst，n2stg

取g9.81m/s2。可以利用下列MATLAB程序画出st在0~0.5范围内n和n的变换曲线：

%Ex2_17.m g=9.81;

for i=1:101 t(i)=0.01+(0.5-0.01)*(i-1)/100;w(i)=(g/t(i))^0.5; tao(i)=2*pi*(t(i)/g)^0.5; end

plot(t,w);gtext('w_n'); hold on;plot(t,tao);gtext('T_n');

王超：基于MATLAB的振动系统编程分析

xlabel('delta_s_t'); title('Example2.1');

Example2.13530252015wn105Tn000.050.10.150.20.250.30.350.40.450.5deltast

2.3单自由度系统的强迫振动

[15]

简谐激励是激励形式中最简单的一种，虽然它在实际中存在的场合比较少但掌握系统对于简谐激励的响应的规律，是理解系统对周期激励或更一般形式激励的响应基础。图所示的弹簧质量系统中，质量块上作用有简谐激振力

P(t)P0sint (2-15)

辽宁工程技术大学毕业设计（论文）

P(t) x P(t) mx m kx cx k c

其中P0为激振力幅，为激振频率。以静平衡位置为坐标原点建立图示的坐标系。从图的受力分析，得到运动微分方程为：

mxcxkxp0sint (2-16) 由常微分方程理论知道，方程(3.2)的通解x由相应的齐次方程的通解xh和非齐次方程的任意特解xp两部分组成，即

x(t)xh(t)xp(t) (2-17)

当欠阻尼时，式中xh(t)为有阻尼自由振动，它的特点是振动频率为阻尼固有频率，振幅按指数规律衰减，称为瞬态振动或瞬态响应；xp(t)是一种持续的等幅振动，它是由于简谐激励振力的持续作用而产生的，称之为稳态强迫振动或稳态振动，在间隔充分长时间考虑的振动就是这种稳态振动，而在刚受到外界激励时，系统的响应则是上述两种振动之和。可见，系统受简谐激励后的响应可以分为两个阶段，一开始的过程称为过渡阶段，经充分长时间后，瞬态响应消失这时进入过渡阶段，经充分长时间后，瞬态响应消失，这是进入稳态阶段。

将方程（2-15）的两端同除以质量m，并且令

cm22n （2-18）其中为相对阻尼系数，n为相应的无阻尼系统的固有频率，则方程（2-15）成为

王超：基于MATLAB的振动系统编程分析

2nxn2xxp0sint （2-19） m上述方程特解可以通过xBsin(t)或者xAcostBsint来求得，这里介绍用复数方法求式（2-19）的特解。先将式（2-19）写为下列的复数形式

2nxn2xxp0ite （2-20） m其中x是复数设复数形式的特解为

xBeit （2-21）

其中B称为复振幅，其意义是包含有相位的振幅。将式（2-21）代入（2-20），解得

BP01m22 （2-22） ni2n记为频率比，它定义为

 （2-23） n则式（2-22）可以写成

BP01p01ik12i2k(12)2(2)2eBei （2-24）

式中

Bp01k(12)2(2)2 （2-25）

tg1212 （2-26）将式（2-24）代入（2-21），得到复数形式的特解为

xBei(t) （2-27）

比较方程（2-17）与（2-18），可知（2-19）中的位移x是（2-20）中复数x的虚部，因此（2-25）的虚部就是方程（2-12）的特解，即有

xBsin(t) （2-28）

其中B为振幅，为相位差。由式（2-26）、2-23）及（2-24）得出稳态强迫振动有如下的基本特点：

辽宁工程技术大学毕业设计（论文）

1)线性系统对简谐激励的稳态响应是频率等同于激励频率而相位滞后与激振力的简谐振动；

2)稳态响应的振幅及相位差只取决于系统本身的物理性质和激振的频率及力幅，而与系统本身进入运动的方式无关。

无阻尼系统对简谐激励的稳态响应可以从式（2-26）得出。当n时，得到1，0，这时

xP01sint （2-29）

k12当n时，得到1，，这时

P01sin(t) （2-23） 2k1P1式（2-21）也可以写成（2-22）的形式，这时相位差反映在振幅0的符号中。

k12x上述结果也可以由直接设xBsint并代入下列方程而得到：

kxP0sint （2-24） mx为了具体讨论影响稳定响应的振幅和相位差的各种因素，记

B0P0 （2-25） kB0实际是质量块在激振力幅静作用下的最大位移。再引入无量纲的振幅放大因子

，它定义为

B1 （2-26） 222B0(1)(2)由式（2-26）和（2-19）可以分别画出以相对阻尼系数为参数的曲线——曲线与曲线，前者称为幅频响应曲线，后者称为相频响应曲线如图所示

程序如下

for kesai=[0.05,0.10,0.15,0.25,0.375,0.50,1.0] lamda=0:0.01:5.0;

beta=1./(sqrt((1-lamda.^2).^2+(2*kesai*lamda).^2)); plot(lamda,beta)

王超：基于MATLAB的振动系统编程分析

hold on end

axis([0 5 0 3]);

300.052.50.1520.25子因0.375大放1.5幅振10.50.51.0000.511.522.533.544.55频率比偏心质量引起的强迫振动振幅放大因子

MB2me(12)2(2)2 (程序如下：

for kesai=[0.05,0.10,0.15,0.25,0.375,0.50] lamda=0:0.01:5.0;

beta=lamda./(sqrt((1-lamda.^2).^2+(2*kesai*lamda).^2)); plot(lamda,beta) hold on end

axis([0 5 0 3]);

2-27)

辽宁工程技术大学毕业设计（论文）

30.52.50.100.1520.25MB/me0.3751.50.5011.00.5000.511.522.5频率比33.544.55

支撑运动引起的强迫振动振幅放大因子

B1(2)2 (2-28) a(12)2(2)2程序如下：

for kesai=[0.05,0.10,0.15,0.25,0.375,0.50,1.0] lamda=0:0.01:5.0;

beta=sqrt((1+(2*kesai*lamda).^2)./((1-lamda.^2).^2+(2*kesai*lamda).^2)); plot(lamda,beta) hold on end

axis([0 5 0 3]);

王超：基于MATLAB的振动系统编程分析

30.052.50.100.152振幅放大因子0.251.50.3750.5011.00.51.141000.511.522.5频率比33.544.55

算例利用MATLAB，绘制弹簧-质量系统在简谐力作用下的响应曲线。已知数据

00.1m/s。如下：m5kg,k2000N/m,F(t)100cos30tN,x00.1m,x系统全解形式如下：

x(t)0xnsinnt(x0f0f0)costcost n2222nn式中，

f0F010020,nm5k20rad/s,30rad/s m利用MATLAB绘制解曲线上式的程序如下： %Ex3_11.m F0=100; wn=20; m=5; w=30;

辽宁工程技术大学毕业设计（论文）

x0=0.1; x0_dot=0.1; f_0=F0/m; for i=1:101

t(i)=2*(i-1)/100;

x(i)=x0_dot*sin(wn*t(i))/wn+(x0-f_0/(wn^2-w^2))*cos(wn*t(i)) +f_0/(wn^2-w^2)*cos(w*t(i)); end plot(t,x); xlabel('t'); ylabel('x(t)'); title('Ex3.11')

Ex3.110.150.10.050x(t)-0.05-0.1-0.15-0.200.20.40.60.81t1.21.41.61.82

2.4本章小结

基于MATLAB对单自由度自由振动绘制振动图像，进行粘性阻尼，强迫振动振幅放大因子绘图进行数据分析，使振动数据更加明显。

王超：基于MATLAB的振动系统编程分析

3 基于MATLAB的多自由度系统编程分析

3.1 多自由度系统

[16]

多自由度振动系统的数学模型：

MxCxKxf （3-1）

其中M、C、K、f和x分别为质量矩阵、阻尼矩阵、刚度矩阵、力向量和响应向量。把这个时域矩阵方程变换到拉氏域（变数为p），并假定初始位移和初始速度为零，则得：

(p2MpCK)X(p)F(p) （3-2）

或 Z(p)X(p)F(p) （3-3）

式中 Z(p)：动刚度矩阵。

由式3-2）或（3-3）可以得出传递函数矩阵H(p)：

X(p)H(p)F(p) （3-4）

借助矩阵相关理论计算出来：

H(p)Z(p)式中 adj(Z(p))：为伴随矩阵； Z(p)：为Z(p)的行列式。

1adj(Z(p))Z(p) （3-5）

式（3-5）的分母，叫做系统的特征方程。类似单自由度系统，特征方程的根，即系统极点，决定系统的共振频率。应为C一般粘性阻尼矩阵，不一定满足矩阵可对角化条件，为了把系统方程（3-2）转化为一般特征值问题公式，需引入恒等式：

(pMpM)X0 （3-6）

将此式与（3-2）合并：

(pAB)YF （3-7）

其中

辽宁工程技术大学毕业设计（论文）

0 AMM, BM0,

0KCpX0Y, F

XF令F=0, 则（3-7）的特征值满足下列方程：

pAB0 （3-8）

对于N自由度系统，此方程有2N个复共轭对出现的特征根：

iiji 其中i阻尼因子；i为阻尼固有频率。 jiii3.2 第一阶固有频率及主振型[17,18]

在求解系统动力响应时，系统较低的前几阶固有频率及相应的主振型占有重要的地位，为计算它们而采用下面的矩阵迭代法是比较简单的。

将i和i带入公式中，得

Aiii （3-9）

若将上式左端看作新列阵，上式表示：对于精确的主振型。新列阵(Ai) 与原来的列阵i的各个对应元素之间都相差同一常倍数，这个常倍数即特征值1。

记X1为初始迭代列阵，由展开定理，X1可以表示为

X1a11a22ann （3-10）

对上式左乘矩阵A，由式（3-9）得知第一次迭代后所得的列阵为

n2X2AX1a111a222annn=1aaan（3-11） 22n1111如果特征值1不是特征方程的重根，那么上式中的

23、、、n都小于1，因此比起111其他主振型1在X2内占的比重相对地比在X1中占的比重大，换句话说，用矩阵A迭代计算一次后，扩大了迭代列阵中第一阶主振型的优势。经第二次迭代后，得

王超：基于MATLAB的振动系统编程分析

2n22a11a2X3AX21an2112n 同理第（r-1）次迭代后的结果为

XrAXr1r1nr12a11a21a2n11r1n （3-12）

可见随着次数的增加，第一阶主振型的优势越来越扩大，当迭代次数充分大时，由上式近似地得

r1Xr1a11 （3-13）

这时再迭代一次，得出

Xr1AXr1Xr （3-14）

由此看到迭代后的新列阵Xr1与原来列阵Xr的各个对应元素之间都仅相差一倍数1，所以Xr或Xr1就是对应于1的第一阶主振型，而特征值1可由下式算出

Xr1l 1Xrl2，，n) （3-15） l(1，其中Xrl表示列阵Xr的第l个元素。为防止迭代过程中迭代列阵的元素变得过大或过小，每次迭代后需要使列阵归一化，例如使它最后一个元素成为1。下面是实用的矩阵迭代法的计算步骤：

1)选取初始迭代列阵X1使其最后一个元素为1；

2)对Xi作矩阵迭代，并使新列阵Y1归一化，即

Y1AX1， X21Y （3-16） Y1n13)重复步骤（2），第r次的迭代结果为

YrAXr， Xr11Yr （3-17） Yrn4)若在允许的误差范围内有Xr1=Xr，则将Xr1或Xr取作第一阶主振型1，由式（3-15）得知

辽宁工程技术大学毕业设计（论文）

YrnYrn （3-18） 1Xrn1因而第一阶固有频率为

11Yrn （3-19）

211由式（3-12）看出，矩阵迭代法计算1及1收敛速度取决于比值2，越小，122收敛的越快，上述的矩阵迭代法又称为逆迭代法。

m k k 图4-2 分析图

m 2k 2m Fig.4-2 analysis chart

用矩阵迭代法求解过程如下：

解用影响系数法求得系统的质量矩阵和刚度矩阵为

2Km00，KKM0m0002m0K3K2K02K 2K算出K的逆阵及系统的动力矩阵为

K11121111，AK1Mm124 122kk122.5125若X1111，第一次迭代后得到

T0.5000004m1 Y1AX17， X2Y10.857000kY1381.000000

重复上述步骤，各次的迭代结果列于表4-1。由表可见，经过6次迭代后已有X7X8，所

王超：基于MATLAB的振动系统编程分析

0.462598以第一阶主振型及基频取为1X70.860806， 1.0000001110.373087k17.184210mm k表4-1 第一阶主振型的迭代

Tab.4-1 The first order iteration of principal mode

2 3 4 5 6 7 1 0.500000 0.465517 0.462830 0.462617 0.462598 0.462598 Xr

1 0.875000 0.862069 0.860911 0.860814 0.860806 0.860806 1 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000

km1 8

7.250000

7.189655

7.184652

7.184245

7.184210

算例分析 clc; clear;

%建立质量矩阵M，刚度矩阵K\\ syms k m; M=[m 0 0 0 m 0 0 0 2*m]; K=[2*k -k 0 -k 3*k -2*k 0 -2*k 2*k];

%*********迭代第n阶主阵型******************** %n为计数器 %

辽宁工程技术大学毕业设计（论文）

n=1; while n<=3

%计算系统动力矩阵A if n==1

A(:,:,n)=K\\M; elseif n~=1

A(:,:,n)=A(:,:,n-1)-(MP(:,n-1)\(:,n-1))*f(:,n-1)*f(:,n-1)'*M; end

%定义初始迭代向量X(1) if n==1

X(:,1,n)=[1 1 1]'; elseif n==2

X(:,1,n)=[1 1 -1]'; elseif n==3

X(:,1,n)=[1 -1 1]'; end

%迭代过程，Y为中间矩阵,i为迭代次数 i=1;

Y(:,1,n)=A(:,:,n)*X(:,1,n); X(:,2,n)=Y(3,1,n).\\Y(:,1,n);

while abs(X(1,i,n)-X(1,i+1,n))>=0.000001&&abs(X(2,i,n)-X(2,i+1,n))>0.000001 Y(:,i+1,n)=A(:,:,n)*X(:,i+1,n); X(:,i+2,n)=Y(3,i+1,n).\\Y(:,i+1,n); X(:,2); i=i+1; end % X(:,:,n) f(:,n)=X(:,i,n); t(n)=Y(3,i,n);

MP(:,n)=f(:,n)'*M*f(:,n);

王超：基于MATLAB的振动系统编程分析

n=n+1; end

%输出数据过程

disp('第一阶主阵型的迭代结果'); X(:,:,1)

disp('第二阶主阵型的迭代结果'); X(:,:,2)

disp('第三阶主阵型的迭代结果'); X(:,:,3)

disp('φi的计算结果，矩阵的每列分别是1,2 ，3阶的'); f

disp('######******注意***********主阵型的迭代结果后面的0是系统的占位符号，不算计算结果');

3.3 本章小结

在工程振动中，确定系统固有频率与主振型时是非常重要的。我们在计算系统固有频率与主振型可以采用矩阵迭代的方法求解，相对于以前所采用的算法语言（Fortran、C语言等），本文采用MATLAB语言，不仅相对语句少，可读性强，还可以利用MATLAB的绘图功能对结果进行直观地分析。只要输入必要的数据，就可以快速地获得振动系统的固有频率以及主振型，对设计人员计算复杂多自由度系统固有频率具有参考意义，并为初步分析各构件的振动情况以及解耦分析系统响应奠定了基础。

辽宁工程技术大学毕业设计（论文）

4 连续系统的振动

4.1 运动微分方程

[19]

考虑长为l的弹性弦或索，每单位长度受大小为f(x,t)的横向力作用，作用在单位上的外力等于作用在单元上的惯性力，即

2(PdP)sin(d)fdxPsindx2 （4-1）

t其中，P是张力；为每单位长度的质量；为弦相对于x轴偏离的角度。对微长度dx，有

Pdx （4-2） xsintan （4-3）

xdP和

2sin(d)tan(d)dx （4-4）

xx2因此非均匀弦受强迫振动的运动微分方程式（4-1）可以简化为

(x,t)2(x,t)[P]f(x,t)(x) （4-5） 2xxt如果弦是均匀的，且张力为常力，则式（4-5）简化为

2(x,t)2(x,t)Pf(x,t)x2t2 （4-6）

如果f(x,t)0，则得自由振动方程为

2(,t)2(x,t)P （4-7）

x2t2或

22c2 （4-8） 2xt2其中

Pc()12 （4-9）

式（4-8）即为著名的波动方程。

王超：基于MATLAB的振动系统编程分析

4.2梁的横向振动

两端固定梁边界条件[20]

Y(0)Y'(0)0Y(l)Y'(l)0

频率方程和主振型

(1)coslchl1

Y(x)Ci[cosixchixri(sinixshix)](2)

rilchilsinilshiicossinlshlliicosilchilclear; clc;

b(1)=4.730041; b(2)=7.853205; b(3)=10.995608; b(4)=14.137166; for i=1:3 c(i)=b(i)./3;

r(i)=((sin(b(i))+sinh(b(i)))./(cos(b(i))-cosh(b(i))));

yx=@(x) cos(c(i)*x)-cosh(c(i)*x)+r(i)*(sin(c(i)*x)-sinh(c(i)*x)); fplot(yx,[0,3]);

title('两端固定，Ci取1，梁长取3。'); hold on; end hold off;

辽宁工程技术大学毕业设计（论文）

两端固定，Ci取1，梁长取3。21.51第二阶振型0.5第三阶振型0-0.5-1第一阶振型-1.5-200.511.522.53两端自由梁边界条件：

Y''(0)Y'''(0)0Y''(l)Y'''(l)0

频率和主振型

(1)coslchl1

Yi(x)Ci[cosixchixri(sinixshix)](2)

rilchil

icossinsinilshililshilcosilshilclear;

clc;

b(1)=4.730041; b(2)=7.853205; b(3)=10.995608; b(4)=14.137166; for i=1:3

c(i)=b(i)./3;

r(i)=((sin(b(i))+sinh(b(i)))./(cos(b(i))-cosh(b(i))));

yx=@(x) cos(c(i)*x)+cosh(c(i)*x)+r(i)*(sin(c(i)*x)+sinh(c(i)*x)); fplot(yx,[0,3]);

title('两端自由，Ci取1，梁长取3。');

王超：基于MATLAB的振动系统编程分析

hold on; end

hold off;

两端自由，Ci取1，梁长取3。21.51第三阶振型0.50-0.5-1-1.5-2第二阶振型第一阶振型00.511.522.53

初始条件[21,22]：

由于运动微分方程涉及对时间t的二阶导数与对x的四阶导数，因而为得到唯一确定得解

(x,t)，需要2个初始条件与4个边界条件为

(x,t0)0(x),自由振动

0(x) （4-10） (x,t0)t可以利用分离变量法求自由振动得解，即令

(x,t)W(x)T(t) （4-11）

将式（4-11）代入式运动微分方程经整理后有

c2d4W(x)1d2T(t)a2 （4-12） 42W(x)dxT(t)dt其中，a2为正的常量式，式（4-12）可以表示为两个式子：

d4W(x)4W(x)0 （4-13） 4dx 28

辽宁工程技术大学毕业设计（论文）

d2T(t)2T(t)0 （4-14） 2dt其中

式（4-14）的解可以表示为

42c2A2EI （4-15）

T(t)AcostBsint （4-16）

其中，A与B为常量，可以根据初始条件确定。为求式（4-13）的解，假定

W(x)Cesx （4-17）

其中，C与s为常量。将式（4-17）代入（4-13）后得

s440 （4-18）

该方程的根为

s1,2，s3,4i（4-19）因此方程（4-13）的解为

W(x)C1exC2exC3eixC4eix（4-20）

其中，C1，C2，C3与C4为常量。式（4-20）也可以表示为

W(x)C1cosxC2sinxC3coshxC4sinhx （4-21）

或

W(x)C1(cosxcoshx)C2(cosxcoshx)C3(sinxsinx)C4(sinxsinhx)（4-22）

在每种不同的形式下，C1，C2，C3，C4为不同的常量。可以由边界条件确定。亮的固有频率可由式（4-15）计算，即

2EIEI （4-23） (l)2AAl4函数W(x)称为梁的固有振型函数，为振动的固有频率。式（4-21）或式（4-22）中的位置常量C1，C2，C3与C4以及式（4-23）中的值可以根据梁的边界条件确定。常见梁的边界条件

王超：基于MATLAB的振动系统编程分析

1)自由端：弯矩和剪力分别为零，即

22 M(x)EI20,V(x)(EI2)0 （4-24）

xxx 2)简支端：挠度和弯矩分别为零，即

2 0,M(x)EI20 （4-25）

x 3)固定端：挠度和转角分别为零，即 0,0 （4-25） x4)梁的两端与弹簧、阻尼器和质量块相连。当梁的末端产生横向位移、转角x、

t 速度t与加速度2t2时，由于弹簧、阻尼器以及质量块所受的阻力分别与，

与2t2成比例，而在该末端阻力由简历来平衡，于是

22(EI2)a[kxcm2] （4-26） xxtt5)梁的末端与扭转弹簧、扭转阻尼器与转动惯性元件相连、这种情况下，边界条件为

223EI2a[ktctI0] （4-27） 2xxxtxt强迫振动[23]

可以运用振型叠加法求梁的强迫振动解。假设梁的挠度为

(x,t)Wn(x)qn(t) （4-28）

n1其中，Wn(x)为第n阶固有振型函数或满足微分方程的特征函数：

d4Wn(x)2EIAWn(x)0，n1,2 （4-29） n4dxqn(t)为对应的广义坐标。将式（4-28）代入强迫振动方程中，得

d4Wn(x)d2qn(t)EIqn(t)AWn(x)f(x,t) （4-30） 42dxdtn1n1根据式（4-29），式（4-30）可以表示为

辽宁工程技术大学毕业设计（论文）

d2qn(t)1W(x)q(t)W(x)f(x,t) （4-30） nnnn2dtAn1n12用Wm(x)乘以式，再从0积分到l，利用正交条件，得

d2qn(t)12q(t)Qn(t) （4-31） nn2dtAb其中，Q(t)称为相对于qn(t)的广义力，其值为

Qn(t)f(x,t)Wn(x)dx （4-32）

0l常量b为

bWn(x)dx （4-33）

0l本质上，式（4-31）可以视为无阻尼单自由度系统的运动方程。运用杜哈美积分，式（4-31）的解可以表示为

qn(t)AncosntBnsinnt1d （4-34）

Abn(t)2f01nanxsinsinsint （4-35）算例作图表示(x,t)22Aln1nll%Ex8_14.m x=20; f0=100; a=10; A=1; l=40;

ro=0.283/386.4; w=100; n=1;

wn=(n^2)*360.393674; for i=1:1001 t(i)=3*(i-1)/100;

w1(i)=(2*f0/(ro*A*1))*sin(n*pi*a/1)*sin(n*pi*x/1)*sin(w*t(i))/(wn^2-w^2);

王超：基于MATLAB的振动系统编程分析

end n=2; for i=1:1001

t(i)=3*(i-1)/1000;

w2(i)=(2*f0/(ro*A*l))*(sin(pi*a/l)*sin(pi*x/l)*sin(w*t(i))/(360.393674^2-w^2)+sin(2*pi*a/l)*sin(2*pi*x/l)*sin(w*t(i))/((2*360.39674)^2-w^2)); end for i=1:1001

t(i)=3*(i-1)/1000;

w3(i)=(2*f0/(ro*A*l))*(sin(pi*a/l)*sin(pi*x/l)*sin(w*t(i))/(360.393674^2-w^2)+sin(2*pi*a/l)*sin(2*pi*x/l)*sin(w*t(i))/((2*360.393674)^2-w^2)+sin(5*pi*a/1)*sin(5*pi*x/1)*sin(w*t(i))/((5*360.393674)^2-w^2)); end

subplot('311'); plot(t,w1); ylabel('w(x,t)'); title('x=20,n=1'); subplot('312'); plot(t,w2); ylabel('w(x,t)'); title('x=20,n=1,2'); subplot('313'); plot(t,w3); xlabel('t'); ylabel('w(x,t)'); title('x=20,n=1,2,5');

辽宁工程技术大学毕业设计（论文）

1x 10-29x=20,n=1w(x,t)0-100.511.5x=20,n=1,222.530.05w(x,t)0-0.050.0500.511.5x=20,n=1,2,522.53w(x,t)0-0.0500.511.5t22.53

4.3小结

本章通过队连续体的固有振型计算建立了数学模型。在计算系统固有频率与主振型后可以采用MATLAB的方法把抽象的理论图像数据化分析，相对于以前所采用的算法语言（Fortran、C语言等），本文采用MATLAB语言，不仅相对语句少，可读性强，还可以利用MATLAB的绘图功能对结果进行直观地分析。本程序对于杆、梁、等弹性振动系统的固有频率以及主振型的计算具有通用性，只要输入必要的数据，就可以快速地获得振动系统的固有频率以及主振型，对设计人员计算复杂多自由度系统固有频率具有参考意义，并为初步分析各构件的振动情况以及解耦分析系统响应奠定了基础。

王超：基于MATLAB的振动系统编程分析

5结论

本文以振动力学为基础对于抽象的振动系统进行了分析和数学模拟，再引用了MATLAB编程软件对单自由度的自由振动数学模型进行图像数据模拟，使振动数据更加明显，对粘性阻尼自由振动的图像实现充分体现了MATLAB的优点是真的数据更加明显清晰，对于单自由的系统的MATLAB实现了振幅放大因子的三个有关图像，而多自由度的振动系统中MATLAB的应用时主振型的计算更加简单速度，连续体的振动同时也应用了MATLAB的公式图像实现使公式更加明显。其中以梁得振动为代表，所以以MATLAB为基础的振动分析使振动系统更见容易理解易算。

辽宁工程技术大学毕业设计（论文）

致谢

四年的本科大学学习生活即将结束，在完成毕业论文之际，深感毕业课题的顺利完成是与许多人对笔者的帮助和关心分不开的，在此写下感激之情。

本课题是在张智慧老师的精心指导和亲切关怀下完成的，从论文的选题到课题的结束都凝聚了导师的心血。导师渊博的知识、严谨的治学作风、精益求精的钻研精神和忘

我工作的精神都值得笔者不断的学习。不但在学术上，在做人上导师正直坦荡的品德也为我树立了榜样。对此，笔者对方老师表示最诚挚的敬意和最衷心的感谢！在课题的实验研究过程中，得到了力学与工程学院各位领导和老师的大力指导和热情帮助。同时感谢力学院的张智慧老师在振动书籍使用上对笔者的大力指导和热情帮助。对理论力学的所有老师和工程力学系的所有老师在理论研究过程中提供的帮助表示感谢。在论文排版过程中，还得到了同学孟村影、赵显伟、丁学成、余海潮的帮助，在此表示感谢。衷心感谢张智慧老师对笔者的大力帮助！

几年来，父母在生活上和事业上给予了笔者莫大的支持和鼓励，付出了辛勤的劳动，谨以此文献给二老。

感谢打扫校园和教室的保洁员阿姨和公寓服务室阿姨，感谢你们给我提供优良的学习和生活环境。

感谢帮助过我的同学们，是你们让我的大学生活丰富多彩。虽然以后见不到了，我还是会想你们的。

再一次向曾经培养、教育、关心和帮助过笔者的领导、前辈、老师和朋友们致以最诚挚、最衷心的感谢！

王超：基于MATLAB的振动系统编程分析

参考文献

[1] 许本文,焦群英.机械振动与模态分析基础[M].北京:机械工业出版社,1988.

[2] 方治华,贾宏玉,沈利.应用Matlab对有死区的振动系统进行仿真分析[J],包头钢铁学院学报,2004,12. [3] 陈予怒,非线性振动,天津科学技术出版社[M],1982.

[4] 傅志方,华宏星,模态分析理论与应用[M],上海交通大学出版社,2002. [5] 周传荣,赵淳生,机械振动参数识别及其应用[M],科学出版社,1989.

[6] 席平原,张海涛,机械动态系统可视化仿真及应用[J],煤矿机械,2004,25（8）:3l-33.

[7] Edward B.Magrab,高会生,李新叶,胡智奇等译,Matlab原理与工程应用[J],电子工业出版社,2002. [8] 张廷中,MATLAB使用指南[M],北京:科学技术文献出版社,1998. [9] 张志涌等,精通MATLAB[M],北京:北京航空航天大学出版社,2000.

[10] Hanselman D,Littlefield B.精通Matlab[M].张航（译）[M]北京:清华大学出版社,2002. [11] 夏永源,张阿舟机械振动问题的计算机解法[M],北京国防工业出版社,1993. [12] 刘延柱等,振动力学[M],北京:高等教学出版社,1998.

[13] 孟立凡,郑宾传感器原理及技术[M],北京:国防工业出版社,2005.

[14] 付巍,基于Madah的单自由度振动系统的数学仿真实验[J],机械工程与自动化,2005(6):23-4. [15] 汤姆逊,达利,振动力学[M],清华大学出版社,2005.

[16] 哈尔滨工业大学理论力学教研室,理论力学（下）[M],北京:高等教学出版社,1997 [17] 程伟,李敏,工程振动基础[Ml.北京:北京航空航天大学出版社,2004. [18] 蒋伟,机械动力学分析[M],北京:中国传媒大学出版社,2004.

[19] (美)Singiresu S. Rao,著李欣业,张明路,编译[M],北京:清华大学出版社,2009. [20] 郑兆昌,机械振动[M],北京:机械工业出版社,1980.

[21] 于长官,现代控制理论[M].哈尔滨:哈尔滨工业大学出版设,2004.

[22] 师汉民,机械振动系统——分析测试建模对策[M],华中科技大学出版社,2004. [23] 程耀东,机械振动学[J],杭州:浙江大学出版社,1988.37~66.

辽宁工程技术大学毕业设计（论文）

附录A译文

第三章简谐激励振动 3.1引导

间歇激励是工程系统中经常遇见的。它通常是由于旋转机械失衡而产生的。尽管纯简单间歇激励似乎比其他周期激励要较少遇见，但为了更好的理解在更一般的激励下系统的相应，明白在间歇激励下系统的特性是必要的。简谐激励可以是系统中的力或系统中某点的位移。

3.2间歇强迫振动

我们将首先考虑受简谐力F0sinl激励的具有粘性阻尼的单自由度系统，如图3.21-1.他的运动微分方程可以根据分离体图求得：

ckxF0sint （3.2-1） mxx

图3.2-1具有简谐激励的阻尼系统图3.2-1有阻尼的强迫振动的矢量关系

这个方程的解分为两部分，一部分是所谓余函数，即齐次方程的解，另一部分是特殊积分。余函数，在本种情况下就是在第二章讨论过的阻尼自由振动。

上述方程的特解是频率与激励频率相同的稳态振动。可以假定特解具有如下形式：

xXsin(t) （3.2-2）

式中X是振荡的振幅，是位移相对于激励F0的相位角。

把方程（3.2-2）代入微分方程（3.2-1），就可以找到上述方程中的振幅和相位角。回想起在简谐运动中，速度和加速度分别比位移导前3.2-2那样图解表示。从图很容易看出：

和，因此方程各项可以同样按图2王超：基于MATLAB的振动系统编程分析

XF0(km2)(c)2 （3.2-3）

tan1c （3.2-4）

km2现在我们把公式（3.2-3）和（3.2-4）写成无因次的形式，这样就能够把结果用简明的图像表示。把公式（3.2-3）和（3.2-4）的分子和分母都除以看，得：

X(1F0/Km2c2)()kk2 （3.2-5）

tanc/k （3.2-6） 2m1k上面的公式还可以进一步表示为下列量的函数：

nk无阻尼振荡的固有频率； mcc2mn临界阻尼系数；

c阻率； cccccc2。 kcckn于是得到了振幅和相位角的无因次表示式：

XkF01[1(22)][2()]2nn2( （3.2-7）

)n （3.2-8） tan21()n这些公式表明，无因次的振幅Xk/F0及相位角仅仅是频率比n和阻率的函数，同时可以图示于图3.2-3.这些曲线指出，在靠近共振的频率范围内，阻率对振幅和相位角

n1和大值n有很大的影响。研究与图3.2-2相应的力图的三种情况，即小值n，（见图3.2-4），可以使我们对振动特性有进一步的理解。

辽宁工程技术大学毕业设计（论文）

对于n1的情况，惯性力和阻尼力都小，这使得相位角也小，外加力的值接近等于弹簧力，如图3.2-4a。

n1时，相位角是

，力图如图3.2-4b。这时惯性力比较大，他与弹簧力相平衡。2外加力克服了阻尼力。共振时的振幅可以根据公式（3.2-5）或公式（3.2-7），或者从图3.2-4b求得：

XF0F0 （3.2-9） cn2k当n1时，接近，外力几乎全花费在克服惯性力上，如图3.2-4c。总括起来，我们可以写出微分方程和它的完整解，包括瞬变项：

F2nxn2x0sintx（3.2-10） m(t)xF0ksin(t)[1(22)][2]2nnX1ensin(12nt1) （3.2-11）

图3.2-3公式（3.27）和（3.28）的图解

王超：基于MATLAB的振动系统编程分析

图3.2-4 强迫振动中的矢量关系

3.3旋转失衡

在旋转机械中，失衡式振动激励的总根源。这里考虑限制在垂直方向移动的弹簧-质量系统，收到旋转部分失衡质量的激励，如图3.3-1.失衡用偏心质量m表示，它的偏心距为e并以角速度旋转。令x代表非旋转部分质量（Mm）离开静平衡位置的位移，则m的位移是：

xesint

图3.3-1由于旋转失衡引起的简谐扰动

因而运动方程为：

d2m2(xesint)kxcx (Mm)xdt上式可重新改写为：

cxkx(me2)sint （3.3-1） Mx显然上式与方程（3.2-1）相同，但在这里以me2代替F0因此稳态解可改写为：

Xme2(kM)(c)222 （3.3-2）

辽宁工程技术大学毕业设计（论文）

tanc （3.3-3） 2kM还可以进一步简化为无因次的形式：

MXme([1(2)n22)][2]2nn2( （3.3-4）

2)n （3.3-5） tan21()n这些公式图解表示于图3.3-2.完整解是：

x(t)X1ensin(1nt1)2me2(kM)(c)222sin(t)（3.3-6）

图3.3-2旋转失衡强迫振动方程（3.3-4）和（3.3-5）的图示

一对反转的偏心重量激振器用来产生弹簧支撑的质量的强迫振动，如图3.3-3。改变转速曾记录下0.60英寸的共振振幅。当转速增加到大大超过共振频率时，振幅出现在靠近不变的振幅0.08英寸处，求系统的阻率。

解：根据公式（3.3-4），共振振幅是：

王超：基于MATLAB的振动系统编程分析

meXM0.60英寸

2当大大地超过n时，同一公式变为：

Xme0.08英寸M

解这两个方程，得系统的阻率：

0.080.066620.60

2meem我们曾指出，偏心距为的质量形式了离心力。这种力根据它们在转子上的分

布而产生静失衡和动失衡。

静失衡当不平衡的质量全部处于一个平面内，如薄圆盘，则不平衡形成的是同一平面的径向力。如图3.3-4所示，这种失衡可以用静力试验法来找；将轮轴安放在一对水平轨道上当轮子比较重的点滚到轮轴最下方一点时轮子就停止了。因此，用不到快速旋转轮子就可找到失衡点，这叫做静平衡法。

图3.3-3 图3.3-4 具有静失衡的系统

动失衡当不平衡质量不只在一个平面内时，形成了一个力和一个偶，这属于动失衡。根据前述方法，静试验法可以找到合理，但如果不快速旋转转子的话，力偶边找不到。例如，我们考虑一个有两个圆盘的轴（图3.3-5），如果两个不平衡质量相等且相对于，则转子将相对轴心处于稳定平衡。但是当转子快速旋转时每一失衡圆盘将产生旋转离心力，使轴相对支撑产生摆动倾向。

一般地说，长转子，如电动机转子和汽车发动机的曲线，可以看成是由一系列个圆盘组成的，每一个都可能有某种失衡。为了找到失衡。这种转子必须用快速旋转法（Spin）。检测和改变转子不平衡的机器叫做平衡机。从本质上看，平衡机是由支撑座构成的，支承

辽宁工程技术大学毕业设计（论文）

座用弹簧安装，以便通过运动来检测不平衡力，如图3.3-6。知道了每一支承的振幅及其相对相位，就能找到转子不平衡并消除它。

图3.3-5 有动失衡的系统图3.3-6 转子平衡机

尽管薄圆盘可以用静力法平衡，但同样可以用动平衡法。我们描述一个可以用简单方法实现的试验。

圆盘被支撑在受弹簧约束的轴承上，它可以像图3.3-7那样沿水平移动。以任意一预定速度转动，注意最大偏斜时的振幅X0以及轮的位置“a”。装在轴承上的加速度计及闭光测频仪可以做这些检测。由初始失衡0而引起的振幅X0被刻在轮子沿o到a的方向。

图3.3-7 进行薄圆盘平衡试验

接着将实验重量1加到轮子的任一点，用同样的速度重复上述过程。由原始失衡0和试验重量1由原始失衡0和试验重量1产生的新的振幅X1和轮位“b”用矢量ob表示。矢量差ab就是试验重量w1的数值增加到w0（oa/ab），矢量ab与矢量oa大小相等，方向相反，这时因为X1为零，因此轮便被平衡。

王超：基于MATLAB的振动系统编程分析

图3.4-1由质量不平衡产生的对称旋转曲线

3.4转轴的弓状旋曲

转轴在一定的速度下会趋向于弓状而且会以复杂的形状旋曲。旋曲的定义为由弯曲轴线和支撑中心线组成的平面的转动。由于各种原因，诸如质量不平衡，轴种的滞后阻尼，陀螺力，轴承中液体摩擦等，都会形成这种现象。轴的旋曲可以发生在轴的转向的同向或反向。旋曲的速度等于或不等于轴的速度。

轴的旋曲是一个难以捉摸的课题。它的一般运动是属于自激运动一类的，在这种中，激励力诱发运动而又由运动本身来控制。轴旋曲的一般运动的研究超出本课程的范围。有兴趣的读者可参阅小冈特完成的有关本课题的优秀报告。

在这一节中我们将考虑同步旋曲这一最简单情况，在这中情况下，旋曲速度等于轴的转速。为此我们假设一个理想的系统：它是一个质量为m的圆盘对称地安装在轴上，轴支撑在两个轴承上，如图3.4-1.盘得质量中心G是在从几何中心S引出的径向距离e处。轴承中心线穿过盘平面的O点，轴中性的变位为OS。在同步旋曲时，O,S和G间保证固定的关系，盘和轴以不变的速度旋转。用xs和ys代表轴心的S的位置，质量中心G的坐标是(xsesint)和（ysesint）。假定粘性阻尼正比于S的速度，则在x，y向的运动微分方程：

d2s m2(xsecost)kxscxdt

辽宁工程技术大学毕业设计（论文）

d2s m2(ysesint)kyscydtscxskxsme2cost （3.4-1） x或 mscxskxsme2cost x m这些方程与方程（3.3-1）相似，我们用观察法可以写成它的解：

xsme2cos(t)(km)(c)222 （3.4-2）

ysme2sin(t)(km)(c)2s2s222

me2OSrxytan(km)(c)222 （3.4-3）

c （3.4-4）

km2显然，线段SGe比位移线段OSr导前一相位角，它与阻尼值和转动速度有关。当旋转速度等于临界速度nk/m（或轴的横向振动的固有频率）时，出现共振的条件，这时的振幅只是由于阻尼才受到限制。图3.4-2表示在三种不同速度条件下的盘-轴系统。

图3.4-2具有粘性阻尼的同步旋曲的振幅和相位关系

3.5 支撑运动

在许多情况下，动力系统是由支撑点的运动而激励，如图3.5-1。我们令y代表支撑点的谐位移，而质量m的坐标x是相对于固定坐标的。

在位移的位置，不平衡力是由阻尼器和弹簧产生，运动微分方程为：

王超：基于MATLAB的振动系统编程分析

k(xy)c(xy) （3.5-1） mx上式可改成为：

cxkxkycy （3.5-2） mx利用复数代数学，我们令

yYejt

xXei(t)Xeieit （3.5-2）

利用复数代数学，我们令

yYeiwxXei(t)Xeeiit （3.5-3）

这样可以使位移x和位移y相差一相位角。把这些公式代入方程（3.5-2），得；

(m2ick)Xei(kic)Y

Xeikic或（3.5-4） 2Y(km)ic因而振幅比的绝对值是：

k2(c)2222(km)(c)1(2XY2)n2[1()2]2[2]nn (3.5-5)

为了求相位角，我们使公式（3.5-4）两边的实数部分和虚数部分各自相等便可确定

sin和cos，然后取两者的比值便得相位角的计算公式：

tanmc32(k[1()2](c)2n221()2(2)nn3)n （3.5-6）

确定稳态振幅和相位的方程式（3.5-5）和（3.5-6）画在图3.5-2。从频率曲线图可以看出，在频率比n2时，不管阻尼如何不同，XY总是等于1。

辽宁工程技术大学毕业设计（论文）

图3.5-2方程（3.5-5）和（3.5-6）的图解

王超：基于MATLAB的振动系统编程分析

附录B

辽宁工程技术大学毕业设计（论文）

王超：基于MATLAB的振动系统编程分析

辽宁工程技术大学毕业设计（论文）

王超：基于MATLAB的振动系统编程分析

辽宁工程技术大学毕业设计（论文）

王超：基于MATLAB的振动系统编程分析

辽宁工程技术大学毕业设计（论文）

王超：基于MATLAB的振动系统编程分析

辽宁工程技术大学毕业设计（论文）

王超：基于MATLAB的振动系统编程分析

辽宁工程技术大学毕业设计（论文）

王超：基于MATLAB的振动系统编程分析

辽宁工程技术大学毕业设计（论文）

王超：基于MATLAB的振动系统编程分析

辽宁工程技术大学毕业设计（论文）

王超：基于MATLAB的振动系统编程分析

中文题目：基于MATLAB的振动系统编程分析

外文题目：ANALYSIS OF VIBRATION SYSTEM BASED ON MATLAB

ROGRAMMING

毕业设计（论文）共 63 页（其中：外文文献及译文完成日期 2011年6月

27页）图纸共0张答辩日期 2011年6月

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文