MATLAB系统仿真报告 - 有阻尼受迫振动系统

一、 问题描述

有阻尼受迫振动的结构及基本原理

图一 有阻尼的受迫振动系统

图1为有阻尼的受迫振动系统，质量为M，摩擦系数为B， 弹簧倔强系数为K。拉力、摩擦力和弹簧力三都影响质量为M的物体的加速度。如果系统的能量守恒，且振动一旦发生，它就会持久的、等幅的一直进行下去。但是，实际上所遇到的自由振动都是逐渐衰减直至最终停止，即系统存在阻尼。阻尼有相对运动表面的摩擦力、液体与气体的介质阻力、电磁阻力以及材料变形时的内阻力等作用。物体在驱动力作用下的振动是受迫振动。

二、 模型分析与建立

利用牛顿运动定律，建立系统的力平衡微分方程如下：

Mdxdt22?Bdxdt?Kx?f(t)

（1）

式中的f (t)是一个外加的激励力，如果 f (t) =F0 sinωt，则称为谐激励力，其中ω为外施激励频率，t是持续时间。故（1）式又可写成：

Mdxdt22

?Bdxdt?Kx?F0sinwt （2）

wn2（2）式是一个线性非齐次方程。令B/M = 2n（n为阻尼系数）），K/M=

n （

wn为固有振动频率），ξ = wn为相对阻尼系数或阻尼比，则（2）式可写为：

dx22 dt?2ndxdt?wnx?hsin(wt)2 （3）

根据阻尼对系统振动的影响，振动响应分为弱阻尼(ξ＜1)、（强阻尼ξ＞1）和临界阻尼（ξ=1）三种情况。这里仅讨论弱阻尼的情况。在弱阻尼情况下的振动为响应：x=Ae-ξwnt sin ( 1－ξ2wn t +φ ) +A1 sin (wt+θ) （4） 谐迫振动的主要特性有：

（1）式(4)包括瞬态与稳态响应两部分，其中瞬态响应是一个有阻尼的谐振。振动频率为系统固有频率wn，振幅A与初相位角?决定初始条件，振幅的衰减按

e??wnt规律，因此，振动持续时间决定于系统的阻尼比?。（2）谐振的稳态响

应是一个简谐振动，其频率比等于激励力的频率w，振幅为A1，相位角为?。（3）当外施激励频率等于系统固有频率时，系统发生位移共振，即振动位移最大。 将（1）式进行Laplace变换得：

MX(s)s?BX(s)s?KX(s)?F(s) 这样，该系统的传递函数为：

H(s)?X(s)F(s)?1Ms22?Bs?K?1/Ms?(B/M)s?(K/M)2 （5）

模型假设及模型建立

如图1所示，已知K= 43.8N/S，M=18.2kg ，B=1.49N/S，n=B/2M=0.0819，h为一常数，固有频率wn=1.5513。在MATLAB的simulink平台下，根据（3）式建立仿真模型文件，如图2所示，其中scope显示振动响应曲线，scope1显示加速度曲线，scope2显示速度曲线。

图2

三、 仿真结果及分析

根据图2，当输入激励力的频率ω作如下变化时，将分别得到其振动响应。

（ 1） 当激励力的频率w=1rad/s时 （小于系统固有频率wn=1.5513rad/s）其响应曲线如图4所示。还可以得到其速度和加速度的响应曲线，如图5、图6。

10.5 10

0-0.5-10

20

30

40

50

60

70

80

图4 w=1rad/s时的振动响应曲线

1 0.5 0 -0.5 -1 -1.5 0 1.510

20

30

40

50

60

70

80

图5 w=1rad/s时的速度曲线

1 0.5 0 -0.5 -1 -1.5 1.50

10

20

30

40

50

60

70

80

图6 w=1rad/s时的加速度曲线

（2）当激励力的频率ω=10rad/s（大于系统固有频率响应曲线如图7所示。

0.060.040.020wn=1.5513rad/s）时，其

0 10

-0.02-0.04 20

30

40

50

60

70

80

图7 w=15rad/s时的振动响应曲线

（3）当激励力的频率等于系统固有频率时，其响应曲线如图8所示，此时为共振

42 10

0-2-40 20

30

40

50

60

70

80

现象。

图8 w=wn=1.5513rad/s时的振动响应曲线

四、 仿真中遇到的问题

仿真之前，没有对仿真参数进行适当的设置，导致在仿真的时候是输出的图形不够光滑，没有达到预期的效果，同时也要对增益进行适当的设置，才能得到好的效果。对系统的积分误差、最大仿真步长与起始仿真步长进行合适设置：最大仿真步长为0.1；初始仿真步长为0.01；绝对误差改为1e-6。然后再进行仿真，就能得到更加符合实际的仿真图，从前面的仿真图中可以看出，减小系统仿真积分误差可以有效地提高系统的仿真性能，使仿真输出波峰与波谷的曲线变得比较光滑。

五、 总结

经过这几天的努力，我终于把simulink这个软件包的基础部分学会了，而且在规定的时间内完成了仿真作业。先前，曾经学习过MATLAB，当时就喜欢上了这个软件，感叹其功能之强大。自己也在平时没事的时候学一点，但没有接触过simulink这个软件包。经过本次simulink仿真的学习，最大的收获是我又学到了MATLAB的其中一个功能，为我以后的学习提供了很好的实验手段和方法。刚一开始，首先英文给了我一个下马威，有一种恐惧的感觉，但是通过老师的精讲细解和自身的努力，我终于一步一步的踏上探索它的征程，其功能之强大吸引着我不断的努力与尝试，虽然辛苦，但苦中有乐。虽然完成了此次作业，但我不会放弃对它的深入学习和研究，真正的把他学透，为我所用成为自己在以后学习与研究中的一个法宝。

附录

附录一：参考资料

【1】matlab程序设计与仿真 刘卫国. 高等教育出版社

【2】matlab7基础与提高 飞思科技产品研发中心. 电子工业出版社 【3】数学模型 姜启源. 高等教育出版社

附录二：习题

1. 利用simulink仿真来实现摄氏温度到华氏温度的转换 Tf?95Tc?32

仿真结果如下：

400 200 0 100 50 0 0

10 20 30 40 50 60 70 80 90 100

上图为华氏温度，下图为摄氏温度

?y??x?y2. 设系统微分方程为?，试建立系统模型并仿真

y(1)?2?

1.2

10.8

0.6

0.4

0.2

0-0.2

0

1

-0.4

-0.6

-0.8

2

3

4

5

6

7

8

9

10

电流图

543210-1-2-3-4-5 0

1

2

3

4

5

6

7

8

9

1

2、脚本文件的方法 定义函数文件hou.m

function xdot=hou(t,x,flag,R,L,C) xdot=zeros(2,1);

xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t); xdot(2)=1/C*x(1); function in=f(t) in=0;

脚本文件：

L=0.5; C=0.02;

for R=[1 2 3 4 5 6 7 8 9 10 11 12 13] [t,x]=ode45('hou',[0 7],[0;1],[],R,L,C); figure(1);

plot(t,x(:,1),'g-');

hold on; xlabel('time'); text(3,0.1,'L(t)'); title('电流零输入响应'); grid; figure(2);

plot(t,x(:,2),'g-'); hold on; xlabel('time'); text(3,0.8,'C(t)'); title('电压零输入响应'); grid; end

电流零输入相应0.150.1L(t)0.050-0.05-0.1-0.15-0.20123time4567

电流零输入相应

电压零输入相应10.8C(t)0.60.40.20-0.2-0.4-0.6-0.80123time4567

电压零输入相应

8. 一池中有水2000m3，含盐 2 kg，以 6m3/ 分 的速率向池中注入浓度为 0.5 kg / m3 的盐水，又以 4 m3 / 分的速率从池中流出混合后的盐水，问欲使池中盐水浓度达到 0.2 kg / 【附加：试画出浓度vsm，需要多长时间？（1用simlink的方法，2用脚本文件的方法）时间的曲线】

1、simulink的方法

解：设t时刻盐水浓度为Ckg/m3 根据题意得方程

C?2?3t?4?Cdt2000?6t?4t3

仿真结果如下图所示：

0.2

0.200

0.200

0.

0.200

0.200

0.1

0.200

0.

0.

0.199

0.199

0.0

0.199

0.199

0

0

510152025

184.184.184.184.184.184.18185.185.185.

有图可以看出盐水浓度达到 0.2 kg / m3，需要184.82分钟。 2、 脚本文件的方法： t=0:250;

v1=6; v2=4; c1=0.5; c2=0.2;

c=(2+3*t)./(2000+6*t); plot(t,c); xlabel('t'); ylabel('C');

title('浓度随时间的变化曲线'); grid on

11、搭建特定的信号源，建立SIMULINK仿真模型、显示仿真结果。

建立sfunction文件

Sfunction的脚本文件：

function [sys,x0,str,ts]=sfunction(t,x,u,flag,m) switch flag,

case 0,

[sys,x0,str,ts]=mdlInitializeSizes; case 3,

sys=mdlOutputs(t,x,u,m); case {1,2,4,9} sys=[]; otherwise

error(['Unhandled flag=',numstr(flag)]); end

function[sys,x0,str,ts]=mdlInitializeSizes(T) sizes=simsizes; sizes.NumContStates=0; sizes.NumDiscStates=0; sizes.NumOutputs=1; sizes.NumInputs=1; sizes.DirFeedthrough=1; sizes.NumSampleTimes=1; sys=simsizes(sizes); x0=[]; str=[]; ts=[0 0];

function sys=mdlOutputs(t,x,u,m) if u>=0&u=m&u<2*m sys=(2*m-u)*(2*m-u); else sys=0; end

在user-defined sfuntion模型库中找到S-Function模型，修改名为sfunction,建立模型如下：

5 4 3 2 1 0 -1- 2- 3- 4- 0 5

1 2 3 4 5 6 7 8 9

1

本文来源：https://www.bwwdw.com/article/005t.html