单摆运动过程的MATLA&Simulink建模与仿真
作者:王军 Email:wj820420@126.com
本文章为一次学习总结，发到网上供大家参考，希望大家转载的时候不要匿名篡改，保持良好的学术作风。
在高中物理学习过程中,我们接触了单摆.当时的单摆定义是:
细线一端固定在悬点,另一端系一个小球,如果细线的质量与小球相比可以忽略,球的直接与线的长度相比也可以忽略,这样的装置就叫做单摆.
根据上面的定义,得出以下试验结论. (1)当摆角很小时，周期与振幅无关; (2)周期与摆球质量无关;
(3)单摆振动的周期与摆长有关;单摆周期的平方与摆长成正比.
以上结论是在理想条件下得到的结论,现对这个理想条件下的单摆进行分析与仿真,将仿真结果与以上结论进行对比验证.
1 理想模式下单摆的数学模型.
首先根据理想条件,摆线质量忽略不计,空气阻力忽略不计.
设摆线长度为l,摆球质量为m,重力加速度为g,系统的初始时刻为t=0,在任意t≥0时刻摆球的线速度为v(t),角速度为ω(t),角位移?(t),以单摆的固定位置为坐标原点建立直角坐标系,水平方向为x轴方向.示意图如下所示：
在 t 时刻，摆锤所受切向力ft(t)是重力mg在其运动圆弧切线方向上的分力，即
ft(t) ??mg sin?(t)
完全理想条件下，根据牛顿第二运动定律，切向加速度为： a(t) ??g sin?(t)
因此得到单摆的运动微分方程组：
dv(t)dt??g sin?(t) (1)
1
d?(t)dt???(t)??v(t)l (2)
使用欧拉算法求解：将dv(t) ??v(t ??dt) ??v(t)和d?(t) ???(t ??dt) ??(t)代入 式（1）及式（2）中，并以仿真步进量Δ作为dt 的近似，得到基于时间的递推方程：
v (t +Δ)= v (t)+gsin?( t)Δ （3）
?(t +Δ)=?(t)?
v(t)lΔ （4）
注:本递推方程仅适合于摆角???2,也就是要求无论初始速度多少,摆角的最大
幅度不能超过90度,如果超过90度比并且初始速度为0时放手小球会自由下落一段时间才能摆动,本递推方程无法描述. 据此编写仿真程序：
在MALAB命令窗口输入以下命令： dt=0.0001; %仿真步进 T=16; %仿真时间长度 t=0:dt:T;%仿真计算时间序列 g=9.8; L=1.5;
th0=1.5; %初始摆角设置,不能超过π/2,即要求球摆动开始时绳子就要有拉力,如果初始摆角超过π/2,则球会经过一阵自由落体后才能进行摆动,上面的递推方程不能满足该情形. v0=0; %初始摆速设置
v=zeros(size(t)); %程序存储变量预先初始化，可提高执行速度 th=zeros(size(t)); v(1)=v0; th(1)=th0;
for i=1:length(t) %仿真求解开始 v(i+1)=v(i)+g*sin(th(i)).*dt; th(i+1)=th(i)-1./L.*v(i).*dt; end %使用双坐标系统来作图
[AX,B1,B2]=plotyy(t,v(1:length(t)),t,th(1:length(t)),'plot'); set(B1,'LineStyle','--'); %设置图线型 set(B2,'LineStyle',':');
set(get(AX(1),'Ylabel'),'String','线速度v(t)m/s');%作标注 set(get(AX(2),'Ylabel'),'String','角位移\\th(t)/rad'); xlabel('时间t/s'); legend(B1,'线速度v(t)',2);
2
legend(B2,'角位移\\th(t)',1);
在以上假设条件下得到仿真图形如下:
在其他条件不变的情况下,仅仅改变细线的长度L=3,再次进行仿真,仿真图象如下:
对比两幅图象可以看出,在理想条件下,同样摆角下,单摆的摆臂变化,影响单摆的最大线速度以及单摆的周期,当摆臂增加时,最大线速度增加,同时单摆的周期也增加.此结论正好与最初单摆理想条件下的试验结论一致.
实际情况中,当摆幅很小(?<50),单摆的摆动可以看作是简谐运动,现在更改单摆的初始摆角为50 情况,也就是令th0=0.085,L依然为1.5,仿真结果如下:
3
2 现实情况下单摆的数学模型.
现实情况下,绳子的质量,摆球的半径,空气的阻力等等都对单摆的摆动有影响,这些影响的主要作用就是阻止单摆的摆动,为简单起见, 可设单摆在摆动中受到阻力fz，显然阻力与摆锤的运动速度有关，即阻力是单摆线速度的函数：fz = f(v)，fz (t) ???kv(t)
上式中，k>0为阻力比例系数，式中的负号表示阻力方向与摆锤运动方向相反。切向加速度由切向合力ft?fz产生，根据牛顿第二运动定律，有
a(t) ??g sin?(t)??
kv(t)m
因此得到修正后的单摆运动微分方程组:
dv(t)dtd?(t)dt??g sin?(t)??
??v(t)lkv(t)m （5）
（6）
仍然使用欧拉算法求解：将dv(t) ??v(t ??dt) ??v(t)和d?(t) ???(t ??dt) ??(t)代入 式（5）及式（6）中，并以仿真步进量Δ作为dt 的近似，得到基于时间的递推方程：
v (t +Δ)= v (t)+(gsin?( t)??
kv(t)m)Δ
?(t +Δ)=?(t)?
v(t)lΔ
据此编写仿真程序：
在MALAB命令窗口输入以下命令： subplot(2,1,1)
4
dt=0.0001; %仿真步进 T=16; %仿真时间长度 t=0:dt:T;%仿真计算时间序列 g=9.8; L=1.5; m=8; k=3;
th0=1.5; %初始摆角设置,不能超过π/2 v0=0; %初始摆速设置
v=zeros(size(t)); %程序存储变量预先初始化，可提高执行速度 th=zeros(size(t)); v(1)=v0; th(1)=th0;
for i=1:length(t) %仿真求解开始 v(i+1)=v(i)+(g*sin(th(i))-k./m.*v(i)).*dt; th(i+1)=th(i)-1./L.*v(i).*dt; end %使用双坐标系统来作图
[AX,B1,B2]=plotyy(t,v(1:length(t)),t,th(1:length(t)),'plot'); set(B1,'LineStyle','-'); %设置图线型 set(B2,'LineStyle',':');
set(get(AX(1),'Ylabel'),'String','线速度v(t)m/s');%作标注 set(get(AX(2),'Ylabel'),'String','角位移\\th(t)/rad'); xlabel('时间t/s'); legend(B1,'线速度v(t)',2); legend(B2,'角位移\\th(t)',1);
其它条件不变,令阻力比例系数k=1.5,将两次仿真的图形放在一起,结果如下:
5