机械振动简谐振动仿真

摘要

机械振动主要有简谐振动，阻尼振动，受迫振动三种。对三种振动建立模型，列出振动方程，再对三种振动给定初始条件，就可以利用Matlab Simulink功能对三种振动进行仿真模拟，得出振动的位移，速度，加速度，动能，势能，机械能随时间的变化关系图像。另外，我们对振动方程求解，得出振子位移关于时间的函数，再分别对其求一阶、二阶导数，就可以得出速度、加速度函数，再经过简单运算就可以得到动能、势能、机械能函数。我们再通过分析函数来分析其图像，再对比仿真模拟出的图像，就可以确定我们的仿真研究方法的可信度。

关键词：简谐振动；阻尼振动；受迫振动；共振

1引言——机械振动的仿真原理

1.1 Matlab Simulink功能简述

Simulink是基于Matlab的框图设计环境，可以用来对各种动态系统进行建模、分析和仿真，它的建模范围广泛，可以针对任何能用数学来描述的系统进行建模，例如航空航天动力学系统、卫星控制制导系统、通信系统、船舶及汽车等，其中包括了连续、离散，条件执行，事件驱动，单速率、多速率和混杂系统等。Simulink提供了利用鼠标拖放的方法来建立系统框图模型的图形界面，而且还提供了丰富的功能块以及不同的专业模块集合，利用Simulink几乎可以做到不书写一行代码即完成整个动态系统的建模工作。除此之外，Simulink还支持Stateflow，用来仿真事件驱动过程。

Simulink是从底层开发的一个完整的仿真环境和图形界面，是模块化了的编程工具，它把Matlab的许多功能都设计成一个个直观的功能模块，把需要的功能模块用连线连起来就可以实现需要的仿真功能了。也可以根据自己的需要设计自己的功能模块，Simulink功能强大，界面友好，是一种很不错的仿真工具[1]。

1.2机械振动的物理模型

物理学中的机械振动主要分为简谐振动、阻尼振动、受迫振动三种。下面我们根据这三种类型的振动建立物理模型来分别研究。

1.2.1简谐振动的物理模型

图1 弹簧振子做简谐振动物理实验模型

如上图所示，弹簧振子在O 附近做简谐振动。已知弹簧振子质量为m ，所受合力为F ，弹簧劲度系数为k ，则有：F kx =-。又由牛顿第二定律有： 22d x F ma m dt

== (1) 于是可以得到： 220d x k x dt m

+= (2) 令m k =2

ω，则可得： 2220d x x dt ω+= (3)

方程(3)的解x 即为弹簧振子在时刻t 时的振动位移，一阶导数x 即为弹簧振子在时刻t 时振动速度，其二阶导数x 即为弹簧振子在时刻t 时的加速度。

1.2.2阻尼振动的物理模型

如图1，若弹簧振子在x 轴上受到粘滞阻尼的作用力，则弹簧振子做阻尼的振动。设弹簧振子受到的阻尼力为： dx f v dt

γγ=-=- (4) 式中γ 为阻尼系数，与物体的形状以及周围性质有关。

弹簧振子受到的弹力为F kx =-，则对弹簧振子，有牛顿第二定律有：

22d d d d x x m kx t t

γ=-- (5) 整理后得： 22d d d d x k x x t m m t

γ=-- (6) 令20k m ω=，2n m

γ=，则有： 2202d d 20d d x x n x t t

ω++=

(7)

这就是阻尼振动的振动方程。其解x 即为弹簧振子在时刻t 时的振动位移，

一阶导数x

即为弹簧振子在时刻t 时振动速度，其二阶导数x 即为弹簧振子在时刻t 时的加速度。

1.2.3受迫振动的物理模型

如图2，弹簧振子在O 附近做阻尼振动。已知弹簧振子质量为m ,弹簧劲度系数为k 。平行于x 轴的平面对弹簧振子有阻尼力的作用。对弹簧振子施加一外加激励力()f t ,设0()sin f t F t ω=，则称为谐激励力，其中ω为外施激励频率，t 是持续时间。 对弹簧振子受力分析，其所受弹力为：F kx =-。由于阻尼振动是振幅(或能量)随时间不断减少的振动。能量减少的原因是有粘滞阻尼和辐射阻尼。为方便，均视为粘滞阻尼。则弹簧

图2 弹簧振子在外加激励力作用下做阻尼受迫振动

振子所受阻尼力为：

dx f v dt

γγ=-=- (8) 式中γ 为阻尼系数，与物体的形状以及周围性质有关。则对弹簧振子，由牛顿第二定律有：

202sin d x dx m kx F t dt dt

γω=--+ (9) 对(9)式变形可得:

202sin F d x k dx x t dt m m dt m

γω=--+ (10) 令00,2,F k n h m m m

γω===，0ω为固有频率，n 为阻尼因数，则(10)可变为： 22022sin d x dx n x h t dt dt

ωω++= (11) 方程(11)的解就是时刻t 时弹簧振子的位移，其一阶导数x 即为弹簧振子在时刻t 时振动速度，其二阶导数x 即为弹簧振子在时刻t 时的加速度[2]。

我们记0/n ξω=为相对阻尼系数或阻尼比。根据阻尼对系统振动的影响，振

动响应分为弱阻尼(ξ＜1)、强阻尼(ξ＞1)和临界阻尼(ξ=1)三种情况，这里仅讨论弱阻尼的情况。

1.3 Matlab Simulink 仿真原理简述

在得到弹簧振子的简谐振动、阻尼振动和受迫振动方程后，通过这三个方程，我们可以用高等数学的方法求出这三个方程的通解。同时，我们可以用Matlab 的计算功能求出它们的通解。这三个方程的通解表示振子位移随时间的变化情况。

我们得到的这三个方程，前两个为二阶常系数线性齐次微分方程，第三个为二阶常系数非齐次微分方程。根据这三个方程，我们可以通过Matlab Simulink 中的各种模块模拟弹簧振子的位移、速度、加速度，再添加一个平方模块，设置好系数，就可以模拟振子动能、势能、机械能，用线连接各模块，这样流程图就做好了。设置好各模块的参数后，再设置好系统环境变量，点击运行，通过示波器模块就可以模拟出相应的图像曲线[3]。图像的横坐标均表示时间，纵坐标相应为位移、速度、加速度、动能、势能、机械能。图像表示这些物理量随时间变化关系。

通过这三种情况方程的通解，我们可以分析振子位移随时间变化情况，再和模拟出的图像对比分析。对方程通解求一阶导，就可以得到振子速度随时间变化关系，分析出速度随时间变化情况，再和模拟出的图像对比分析。同样我们可以求出方程通解的二阶导数，这就是振子加速度随时间变化关系，分析函数特征，再和模拟出的函数图像对比分析，就可以分析模拟出的图像是否正确，及其与理

论符合情况。

2简谐振动方程的解及其模拟仿真

2.1简谐振动方程的求解

这里，我们设系统初始条件为0t =s 时，04m x =，00m/s v =。通过高等数学方法解这个齐次微分方程可得：

12cos sin )cos()x C t C t t A t ωωω?ω?=+=+=+ (12)

式中21

tan ,C A C ?=-。则速度表达式为:sin()v x A t ωω?==-+，将初始条件代入(12)式，可得：

4cos()x t ω?=+ (13)

这就是满足初始条件的简谐振动方程的解。由(13)式我们可以得出弹簧振子位移随时间的变化情况。振子周期为2T π

ω=。0s t =时，振子位移正向最大位移

出，即图1中的A 位置，此时振子速度为0，加速度最大；经

4T ，振子向负方向运动到平衡位置，此时振子速度最大，加速度为0；再经4

T ，振子继续向负方向运动到负的最大位移处，此时速度为0，加速度最大；再经过

4T ，振子向正向运动到平衡位置，此时速度最大，加速度为0；最后经过

4

T ，振子回到初始位置，即正的最大位移处，完成一个周期的振动。 通过matlab7.0符号运算，可以得出该微分方程的解，相关程序见附录程序1。

2.2简谐振动模型的仿真研究

2.2.1基本模型的建立

我们设系统的固有频率1ω-=，则2200.5s ω-=。于是(3)式变为：

22

0.50d x x dt +=。打开Simulink Library Browser ，选择新建按钮，根据所

需要模拟的运动方程选取模块，其中包括Subtract、Intergrator、Gain以及Scope 模块，需要注意的是将Subtract模块中的List of signs改为--，以便让x前面的符号为负，为了使x前的系数为0.5，将Gainl中的值设为0.5，速度项系数Gain设为0。设置位移模块的初始值为4，速度模块的初始值设为0，加速度模块的初始值设为0。这样，几个关键模块的属性就根据方程的需要设置好了。

(1)运用Gain1将x和0.5相乘。

(2)运用Subtract使x前面的符号为负。

(3)运用Intergrator将x积分为x，将x积分为x。

Scope为示波器输出模块。最后，将各个模块按照方程的需要逐一连接，如

图3所示[4]。

图3简谐振动位移仿真模拟流程图

图4简谐振动位移仿真模拟图像

单击右键，选择Configuration Parameters设置系统的运行环境，初始运行时间设为0s，停止时间取为80s,最大步长设为0.1，初始步长设为0.01，设好后，点击OK。再点击图3“Scope”输出模块得到振子位移仿真曲线如图4所示。图中横坐标表示时间，单位为s，纵坐标表示位移，单位为m。

2.2.2 速度、加速度的监测

要得到速度与加速度的实时振动曲线只需要在图3的基础上加入两个Scope 模块，如图5所示。

图5弹簧振子速度、加速度仿真模拟流程图

运行Scope1得到弹簧振子速度图像，如图6所示。图中横坐标表示时间，单位为s，纵坐标表示振子速度，单位为m/s。

图6 弹簧振子速度仿真模拟图像 图7 弹簧振子加速度仿真模拟图像

运行Scope2，就得到弹簧振子加速度监测图像，如图7所示。图中横坐标表示时间，单位为s ，纵坐标表示振子加速度，单位为2m/s 。

2.2.3 动能、势能、机械能监测

系统动能、势能、机械能的定义如下： 212k E mv =

(14) 212

p E kx = (15) p k t E E E += (16) 根据动能与势能的公式在原有的简谐振动模拟流程图中加入Product 模块(实

现2x 和2v 运算)和增益模块Gain 以及Sum 模块将两输入信号进行叠加便可将动能与势能及机械能波形输出出来。先对各个模块名进行编辑，设置好字体大小，再进行各个模块的属性设定。

前面我们已经设定了x 的系数即Gain1参数为0.5，即：0.5k m =，这里，我们取1k g m =,则0.5N/m k =。我们由此可得到动能中速度的平方项系数

20.5

m =，即Gain2参数为0.5；势能中x 平方项系数20.25k =,即Gain3的参数为0.25。最后用仿真信号线将各个模块连接起来，如下图8所示[5]。

图8 简谐振动弹簧振子动能、势能、机械能流程图

单击运行后，点击“t E ”输出模块得到总能量曲线，点击“k E ”输出模块得到动能曲线，点击“p E ”输出模块得到势能曲线如图9，图10，图11所示,图

中横坐标代表时间，单位为s ，纵坐标分别代表动能、势能、机械能，单位为J 。

图9简谐振动弹簧振子动能图像图10简谐振动弹簧振子势能图像

图11简谐振动弹簧振子机械能图像

2.3简谐振动的图像分析

由简谐振动方程的解(13）式知：位移随时间的变化关系为余弦函数。即：4cos()x t ω?=+。由初始条件知其振幅4m A =，初始位移为4m 。周期28.88s

T π

ω===。正如图4所示。这就是弹簧振子做简谐振动的位移随时

间的变化关系。

我们对(13)式求一阶导数，有：

4sin()))x t ωω???=-+=-+=-+ (17) 这就是弹簧振子的速度随时间变化的关系。它的图像是正弦函数图像，周期为28.88s

T π

ω===。根据初始条件，其初始速度为0，正如图6所示。理

论与图像相符合。

我们再对(13)式求二阶导数，有：

24cos())x t ωω??=-+=-+ (18) 这就是弹簧振子加速度随时间变化关系。它的图像是余弦函数，周期仍为28.88s

T π

ω===。根据初始条件，其初始加速度为22m/s -，正如图7所示。

理论与图像是相符合的。

我们将(14)式与(17)式联合，可以得到振子动能随时间变化关系如下：

2214sin )2

k E mx ?==+ (19) 其图像是将正弦函数负半轴部分沿t 轴对折上去后得到的。很容易看出上式必为非负，故图像在时间轴上方。其周期变为原来周期的一半，即：4.44s

T πω===。由于初始时刻速度为0，故初始动能为0，正如图9所示。理论与图像是相符合的。

同理，我们将(15)式与(13)式联合，可以得到弹簧振子势能随时间变化关系如下：

2214cos )2

p E kx ?==+ (20) 由上式可知，振子势能图像是将余弦函数负半轴部分沿t 轴对折上去得到的，上式比为非负，图像在时间轴上方。其周期与动能周期一样，均为4.44s 。由初

始条件知，其初始势能最大，为4J ，如图10所示。可见，理论与图像是相符合的。

我们将(19)式和(20)式相加，就可以得到振子的机械能为：4t k p E E E =+=J 。其图像为平行于时间轴的一条直线，该直线在纵轴上的截距为4J ，如图10所示。图像与理论是符合的。我们得出振子的机械能为一定值。从能量角度分析，做简谐振动的振子只受弹力作用，系统机械能守恒。

3阻尼振动方程的求解和仿真模拟

3.1弹簧振子做阻尼振动方程的求解

方程(7)是齐次方程。其特征方程为：

22020n λλω++= (21) 这里我们只讨论相对阻尼系数(或阻尼比)01n ξω=<的情况，即为弱阻尼振动。于是可以求得特征根为：

1n λ=-+

(22)

2n λ=--(23)

于是我们可以得到方程(7)的齐次方程通解为：

12()[)sin ]nt x t e C C -=+

(24)

我们用也可以用Matlab 的数学计算功能编写程序求解方程(7)，相关程序见附录程序2。

3.2弹簧振子做阻尼振动的模拟仿真研究

如图1所示，我们假设有粘滞阻尼力时，k =43.8N/m ,18.2kg m =,-11.49N/(ms )γ=，10.0819N/(ms kg)2n m γ

-==，0ω=1.5513-1s 。由此我们可计算

得到：120.1638N/(ms kg)n -=，20 2.4065ω=-2s 。与简谐振动的仿真模拟流程图类似，只是x 的系数不再为0，而应是0.1638。x 的系数为2.4065。初始条件设为004m,0m/s x v ==，连接好个模块后，如图12所示：

图12做阻尼振动的弹簧振子仿真模拟流程图

弹簧振子的位移，速度，加速度仿真模拟图像分别如图13、14、15所示，图中横坐标表示时间，单位为s，纵坐标分别表示位移、速度、加速度，单位分别为m、m/s、m/s2。

图13弹簧振子做阻尼振动的位移仿真图像

图14弹簧振子做阻尼振动的速度仿真图像图15弹簧振子做阻尼振动的加速度仿真图像

类似于简谐振动的能量仿真模拟流程图，只需修改相关系数即可得出弹簧振子做阻尼振动的动能、势能、机械能模拟图像。动能表达式中速度平方项的系数

为1

0.518.29.1

2

m=?=，即Gain4的Gain值设为9.1；势能表达式中位移平方项

的系数1

0.543.821.9

2

K=?=，即Gain3中Gain值设为21.9。其流程图如图16

所示，动能、势能、机械能分别如图17、18、19所示。图中横坐标均表示时间，单位为s，纵坐标分别表示动能、势能、机械能，单位均为J。

图16弹簧振子做阻尼振动的能量仿真模拟流程

图17阻尼振动的动能仿真模拟图像图18阻尼振动的势能仿真图像

图19 阻尼振动的机械能仿真模拟图像

3.3阻尼振动的图像分析

根据图像模拟过程，我们已知的条件有：k =43.8N /m ，18.2kg m =，-11.49N/(ms )γ=，10.0819N/(ms kg)2n m γ

-==，0ω=1.5513-1s ，

120.1638N/(ms kg)n -=，20 2.4065ω=-2s ，004m,0m/s x v ==。

阻尼振动方程的解即(24)式是振子位移随时间变化关系。我们将它稍作变形为：

()))nt x t e Ae ??--=+=+

(25)

式中：A =21

tan C C ?=-。上式可以分为两部分，第一部分为nt Ae -，表征阻尼振动振子的振幅；

第二部分为)?+，表示阻尼振动位移随时间呈余弦函数变化。很显然，第一部分是随时间逐渐减小的，这说明振子的振幅在逐渐减小。由初始条件可知：振子初始位移为4m

。振子的位移变化周期为4.054s T ==。由于振子振幅随时间增大逐渐减小，因此初始位移为最

大位移。以后，随时间增大振子的振幅逐渐减小。这与阻尼振动的概念是相符的。图13表示阻尼振动的位移随时间变化的图像，这与理论是相符的。

我们再对(25)式求一阶导数，有：

'0())nt v x t A e ω?-==-+ (26) 式中，'0???=+

，0?=-。上式表示振子速度随时间变化关系。

与振子位移与时间关系类似，(26)式也可以分为两部分。第一部分为0nt A e ω--，表征振子速度极值；

第二部分为')?+，表征振子速度随时间呈正弦函数变化。我们容易看出，第一部分的绝对值是随时间变化逐渐减小的，即速度大小的极值是随时间逐渐减小的。由初始条件知，其初速度为0。由第二部分可

知，振子的周期为 4.054s T ==。如图14所示，这就是振子做阻尼振

动的速度随时间变化的图像，图像与理论是相符合的。

我们将(26)式再对时间求一阶导数，就可以得到振子加速度随时间变化关系为：

2"0())nt a x t A e ω?-==-+ (27)

式中，"02???=+。同样可以将(27)式分为两部分。第一部分为20nt A e ω--，表征

加速度的极值情况；第二部分为")?+，表征振子加速度随时间呈余弦函数。容易看出，第一部分的绝对值随时间增大逐渐减小，即加速度大小的极值是逐渐减小的。初始时刻，其加速度大小的极值为最大值。从第二部分可以

看出，其周期仍为 4.054s T ==。图15是我们模拟出的振子加速度随

时间变化图像，这与理论是相符合的。

由于振子做阻尼振动的位移，速度的极值均随时间的增大而减小的，故其势能，动能，机械能的极值亦随时间增大而减小。一定时间后，它们都趋于0。分别如图17、18、19所示。我们从能量角度来分析振子的机械能。振子的弹力做功不改变其机械能，但振子所受的阻尼力一直做负功，所以振子的机械能不断减小，最后为0。