用MATALAB分析阻尼运动

1、 Matlab语言课程论文用MATALAB分析阻尼运动 姓名：刘莎学号：12010245231专业：2010级通信专业指导老师：汤全武 所在学院：物理电气信息学院完成日期：2011年11月28号 用MATLAB分析阻尼运动（刘莎 12010245231 2010级通信专业） 【摘要】：大学物理中的阻尼运动是一个比较繁琐的问题，物体在运动过程中受各种阻力（如摩擦力、空气阻力等）的影响，出现能量逐渐衰减而导致运动减弱的现象，这种运动被称为阻尼运动。用MATLAB语言解决是个很好的工具，可以解决方程、绘制图像、处理数据，这大大简化运算过程与方法，在实验处理方面给予了很大的帮助，节约了时间，同时也节约

2、了计算量。也可以把MATALAB应用到实际问题当中，达到了学以致用，加深对阻尼运动的理解。【关键词】：力学 阻尼运动 matlab语言 绘图一 、问题提出MATLAB语言是当今科学界最具影响力也是最具活力的软件。是一种面向与工程的高级语言，它提供了强有力的科学用算，灵活的设计流程、高质量的图形可视化以及界面设计的便捷和其他语言接口等功能，MATLAB语言在各国高校与研究中发挥着重要的作用，它是一种集数据运算、符号运算，可视化建模、仿真和图形处理多方面与一体的语言，现在用MATLAB解决物理中的阻尼运动也是很方便的事情。可以用传统的方法通过MATLAB实现其振动规律曲线，并改变其初始条件，实现一

3、组曲线的脉冲过度函数。二、 模拟小球受阻尼运动通过把阻尼振动的小球的运动做投影，更好地观察在阻尼振动过程中，小球的振动衰弱情况，用ma tlab把阻尼振动简化成较容易的程序设计。小球受到阻力运动状态就会减弱，利用这种情况就可以分析阻尼运动的过程 。模拟弹簧阻尼振动的过程包括以下内容：（） 更新小球和弹簧的位置。（）更新水平线的位置和长度。（）实时地画出弹簧高度与时间的关系曲线。（）小球受到各种阻尼时，合外力和加速度都会减小，运动会受到阻碍，用MATLAB进行模拟。下面来模拟弹簧振子在阻尼力的作用下进行阻尼运动的过程，同时画出相应圆球质心随时间变化的曲线。相应的MATLAB程序如下：close

4、all;clear;clc; % 阻尼振动rectangle('position',12,8,2,0.3,'FaceColor',0.1,0.3,0.4); % 生成固定块axis(0,15,-1,10);hold on; % 设置坐标轴范围plot(13,13,7,8,'r','linewidth',2); % 画与弹簧连接的线y=2:.2:7; % 得到弹簧对应的纵坐标数据M=length(y); % 获取数据的长度x=12+mod(1:M,2)*2; % 生成弹簧的横坐标数据x(1)=13;x(end-3:end)=13;

5、% 计算出弹簧上下端点的横坐标值D=plot(x,y); % 画出弹簧C=0:.1:2*pi; % 生成圆球的角度数据r=0.3; % 圆球的半径t1=r*sin(C); % 计算出圆球对中心的纵坐标数据F1=fill(13+r*cos(C),2+t1,'r'); % 画出圆球set(gca,'ytick',0:2:9); % 设置y轴的刻度set(gca,'yticklabels',num2str(-1:3'); % 重新设置y轴的刻度值plot(0,15,2,2,'black'); % 画出平衡位置H1=plot(0,

6、13,2,2,'g'); % 球心的跟踪线Q=plot(0,2.5,'color','r'); % 画出运动曲线;td=; % 记录时间的变量yd= ; % 记录y轴位置的变量T=0; % 设置初始时的时间值text(2,8,'damp osillation','fontsize',24); % 添加标注文字set(gcf,'doublebuffer','on'); % 设置渲染效果while T<12; % 利用循环处理阻尼运动过程的模拟pause(0.2); % 暂停一下,

7、显示动画效果Dy=1-0.5*exp(-T/4)*cos(pi*T); % 计算T时刻弹簧对平衡位置的位移Y=-(y-2)*Dy+7; % 计算弹簧的纵坐标数值Yf=Y(end)+t1; % 计算圆球的纵坐标数值td=td,T;yd=yd,Y(end); % 更新运动曲线的数据set(D,'ydata',Y); % 更新弹簧的位置数据set(F1,'ydata',Yf,'facecolor',rand(1,3); % 更新圆球的位置数据set(H1,'xdata',T,13,'ydata',Y(end),Y(end

8、); % 更新跟踪线的数据set(Q,'xdata',td,'ydata',yd); % 更新圆球的数据T=T+0.1; % 更新时间end %结束Kd=find(diff(sign(diff(yd)=-2)+1; % 计算极大值的位置X=td(Kd); % 得到极大值处的横坐标数值Y=yd(Kd); % 得到极大值处的纵坐标数值X=0,X,td(end); % 得到上侧包络线的横坐标数据Y=yd(1),Y,yd(end); % 得到上侧包络线的纵坐标数据plot(X,Y,':'); % 画出上侧包络线Kx=find(diff(sign(diff

9、(yd)=2)+1; % 计算极小值的位置X=td(Kx); % 得到极小值处的横坐标数值Y=yd(Kx); % 得到极小值处的纵坐标数值X=0,X,td(end); % 得到下侧包络线的横坐标数据Y=-(yd(1)-4),Y,-(yd(end)-4); % 得到下侧包络线的纵坐标数据plot(X,Y,':'); % 画出下侧包络线运行动画截图如图所示： 图（一）运行程序：通过MATLAB就可以清楚的看到小球受到阻尼运动时的运动情况，图像成衰减型，利用matlab就可以简单、省时的完成所需要的内容。三、小球受阻尼力作用下的圆周运动小球受阻尼力作用下的圆周运动，同时小球的大小随着

10、时间以e指数的规律减小。模拟这个运动过程需要含以下三方面的内容：计算小球在不同时刻的位置记录小球运动形成的轨迹实时更新小球的大小通过上述三方面的模拟可以得到这个过程的动画显示，通过MATLAB编写小球受到阻尼运动时的圆周运动程序如下：close all;clc;clear;n=1; % 阻尼系数w=6; % 圆频率ww=sqrt(w2-n2); % 等效的圆频率B1=80; % 系数B2=-20; % 系数axis(-105,105,-105,105);hold on; % 设置坐标轴范围set(gcf,'doublebuffer','on'); % 设置图形窗

11、口的渲染效果G=plot(100,0,'r'); % 轨迹曲线x=100; % 初始位置的横坐标y=0; % 初始位置的纵坐标r=2; % 小球半径t1=r*sin(0:.1:2*pi); % 计算小球的边界数据t2=r*cos(0:.1:2*pi); % 计算小球的边界数据X=x+r*t1; % 得到小球的横坐标数据Y=y+r*t2; % 得到小球的纵坐标数据Q=fill(X,Y,rand(1,3); % 画出小球t=0; % 时间的初始数值plot(-100,100,0,0); % 画出格线plot(0,0,-100,100); % 画出格线k=1; % 控制循环的参数xl

12、abel('Please press "space" key and stop this program!',.'fontsize',14,'color','r'); % 添加提示字符title('damp turbination ','B1=',num2str(B1),.' B2=',num2str(B2),'fontsize',14); % 加注图题while k; % 循环处理s=get(gcf,'currentkey');

13、% 获取当前键入的按键名称if strcmp(s,'space'); % 判断按键是否为空格键clc;k=0; % 把参数k设置为0，从而可以结束循环endpause(0.3); % 暂停一下显示动画效果z=exp(-n*t)*(B1*exp(i*ww*t)+B2*exp(-i*ww*t); % 计算当前时刻的位置xp=real(z); % 分离出实部数值yp=imag(z); % 分离出虚部数值r=0.5+1.5*exp(-n*t); % 计算当前小球的半径X=xp+t1*r; % 计算小球的横坐标数据Y=yp+t2*r; % 计算小球的纵坐标数据set(Q,'xda

14、ta',X,'ydata',Y); % 更新小球的数据x=x,xp; % 存储轨迹的横坐标数据y=y,yp; % 存储轨迹的纵坐标数据set(G,'xdata',x,'ydata',y); % 更新轨迹曲线的数据t=t+0.02; % 更新时间数值end 小球受阻尼做圆周运动动画截图如图所示。 图（二）四、阻尼运动的运动方程分析单自由度阻尼系统的阻尼系统对其固有振动模态的影响，用传统方法通过matlab实现其振动规律曲线，并改变其初始条件，实现一组曲线的脉冲过度函数，改变传统运算算法，使运算简便明了。（a）动力学分析：现象：振幅随时间减小

15、原因：阻尼假设：振动速度减小时，摩擦力正比于质点的速率。即：阻尼力Fr=-Cv 阻尼系数：C对物块应用牛顿第二定律： （1） (2)化解为二阶线性常系数齐次方程，即阻尼振动的动力学方程： (3)其中： 固有角频率 (4) 阻尼系数 （5）解得： （6）其中： 振幅 角频率 (7)振幅为随时间的推移，呈指数递减，越大，振动衰减越快；越小，振幅衰减越慢。随时间的推移，质点坐标单调地趋于零。质点运动是非周期的，甚至不是往复的。将质点移开平衡位置后释放，质点便慢慢回到平衡位置停下来，即过阻尼状态。（b）种阻尼的比较欠阻尼运动 过阻尼 临界阻尼 对以上传统公式算法，用matlab进行编程计算，实现其振动

16、规律曲线，并通过曲线对阻尼振动规律进行分析。所绘图如下图：曲线对阻尼振动规律进行分析的图如图所示： 图（三） 图(四)(c)把上述算法改进改进上述算法公式较为繁琐，计算速度慢，编程复杂，是因为要避免复数运算。而matlab本身就具备复数运算功能，通过改变其中一些算法，重新编程，可以使求解变得简单，程序也直观简单。设原始微分方程为： （1）先用roots函数求p1,p2，语句为： （2）然后不必管他们是否为复数， 把 （3） 对t求导，得 （4）带入初始条件可得: 将这两个线行方程联立，得解： （6）这样，只要给出t数组，就求出了x。根据以上改进，假如其它参数不变，求出=0.3一种情况下此系统的

17、脉冲响应，用matlab编程实现其运动规律曲线，如下图：运动规律曲线图如图所示 图（五）五、 用matlab仿真有阻尼的受迫振动 图（六） 有阻尼的受迫振动系统 图六为有阻尼的受迫振动系统，质量为M，摩擦系数为B， 弹簧倔强系数为K。拉力、摩擦力和弹簧力三都影响质量为M的物体的加速度。如果系统的能量守恒，且振动一旦发生，它就会持久的、等幅的一直进行下去。但是，实际上所遇到的自由振动都是逐渐衰减直至最终停止，即系统存在阻尼。阻尼有相对运动表面的摩擦力、液体与气体的介质阻力、电磁阻力以及材料变形时的内阻力等作用。物体在驱动力作用下的振动是受迫振动。（a）模型分析与建立 利用牛顿运动定律，建立系统的

18、力平衡微分方程如下： （7） 式中的f (t)是一个外加的激励力，如果 f (t) =F0 sint，则称为谐激励力，其中为外施激励频率，t是持续时间。故（1）式又可写成： （8）（8）式是一个线性非齐次方程。令B/M = 2n（n为阻尼系数），K/M= （ 为固有振动频率）， = 为相对阻尼系数或阻尼比，则（8）式可写为： （9）根据阻尼对系统振动的影响，振动响应分为弱阻尼(1)、（强阻尼1）和临界阻尼（=1）三种情况。这里仅讨论弱阻尼的情况。在弱阻尼情况下的振动为响应：x=Ae-wnt sin ( 12wn t + ) +A1 sin (wt+) （10）谐迫振动的主要特性有：（1）七式包

19、括瞬态与稳态响应两部分，其中瞬态响应是一个有阻尼的谐振。振动频率为系统固有频率，振幅A与初相位角决定初始条件，振幅的衰减按规律，因此，振动持续时间决定于系统的阻尼比。（2）谐振的稳态响应是一个简谐振动，其频率比等于激励力的频率w，振幅为，相位角为。（3）当外施激励频率等于系统固有频率时，系统发生位移共振，即振动位移最大。 将（7）式进行Laplace变换得： （11） 这样，该系统的传递函数为： （12）模型假设及模型建立 如图1所示，已知K= 43.8N/S，M=18.2kg ，B=1.49N/S，n=B/2M=0.0819，h为一常数，固有频率=1.5513。在MATLAB的simulin

20、k平台下，根据（3）式建立仿真模型文件，如图2所示，其中scope显示振动响应曲线，scope1显示加速度曲线，scope2显示速度曲线。图(七）(b)仿真结果及分析 根据图2，当输入激励力的频率作如下变化时，将分别得到其振动响应。 (1） 当激励力的频率w=1rad/s时 （小于系统固有频率=1.5513rad/s）其响应曲线如图4所示。还可以得到其速度和加速度的响应曲线，如图5、图6。01020304050607080-1-0.500.51图(八) w=1rad/s时的振动响应曲线01020304050607080-1.5-1-0.500.511.5 图（九） w=1rad/s时的速度曲线01020304050607080-1.5-1-0.500.511.5 图（十） w=1rad/s时的加速度曲线（2）当激励力的频率=10rad/s（大于系统固有频率=1.5513rad/s）时，其响应曲线如图7所示。01020304050607080-0.04-0.0200.020.040.06 图（十一） w=15rad/s