单摆运动过程的MATLABSimulink建模与仿真.docx
《单摆运动过程的MATLABSimulink建模与仿真.docx》由会员分享，可在线阅读，更多相关《单摆运动过程的MATLABSimulink建模与仿真.docx（14页珍藏版）》请在冰豆网上搜索。
单摆运动过程的MATLABSimulink建模与仿真
单摆运动过程的MATLAB&Simulink建模与仿真
单摆运动过程的MATLA&Simulink建模与仿真
作者:
王军Email:
wj820420@
本文章为一次学习总结，发到网上供大家参考，希望大家转载的时候不要匿名篡改，保持良好的学术作风。
在高中物理学习过程中,我们接触了单摆.当时的单摆定义是:
细线一端固定在悬点,另一端系一个小球,如果细线的质量与小球相比可以忽略,球的直接与线的长度相比也可以忽略,这样的装置就叫做单摆.
根据上面的定义,得出以下试验结论.
（1）当摆角很小时，周期与振幅无关;
（2）周期与摆球质量无关;
（3）单摆振动的周期与摆长有关;单摆周期的平方与摆长成正比.
以上结论是在理想条件下得到的结论,现对这个理想条件下的单摆进行分析与仿真,将仿真结果与以上结论进行对比验证.
1理想模式下单摆的数学模型.
首先根据理想条件,摆线质量忽略不计,空气阻力忽略不计.
设摆线长度为l,摆球质量为m,重力加速度为g,系统的初始时刻为t=0,在任意t≥0时刻摆球的线速度为v（t）,角速度为ω（t）,角位移
（t）,以单摆的固定位置为坐标原点建立直角坐标系,水平方向为x轴方向.示意图如下所示：
在t时刻，摆锤所受切向力ft（t）是重力mg在其运动圆弧切线方向上的分力，即
（t）=mgsin
（t）
完全理想条件下，根据牛顿第二运动定律，切向加速度为：
a（t）=gsin
（t）
因此得到单摆的运动微分方程组：
gsin
（t）
（1）
（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）-
Δ（4）
注:
本递推方程仅适合于摆角
也就是要求无论初始速度多少,摆角的最大幅度不能超过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;
fori=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）;
legend（B2,'角位移\th（t）',1）;
在以上假设条件下得到仿真图形如下:
在其他条件不变的情况下,仅仅改变细线的长度L=3,再次进行仿真,仿真图象如下:
对比两幅图象可以看出,在理想条件下,同样摆角下,单摆的摆臂变化,影响单摆的最大线速度以及单摆的周期,当摆臂增加时,最大线速度增加,同时单摆的周期也增加.此结论正好与最初单摆理想条件下的试验结论一致.
实际情况中,当摆幅很小（
<50
）,单摆的摆动可以看作是简谐运动,现在更改单摆的初始摆角为50情况,也就是令th0=0.085,L依然为1.5,仿真结果如下:
2现实情况下单摆的数学模型.
现实情况下,绳子的质量,摆球的半径,空气的阻力等等都对单摆的摆动有影响,这些影响的主要作用就是阻止单摆的摆动,为简单起见,可设单摆在摆动中受到阻力fz，显然阻力与摆锤的运动速度有关，即阻力是单摆线速度的函数：
fz=f（v），fz（t）=-kv（t）
上式中，k>0为阻力比例系数，式中的负号表示阻力方向与摆锤运动方向相反。
切向加速度由切向合力ft-fz产生，根据牛顿第二运动定律，有
a（t）=gsin
（t）-
因此得到修正后的单摆运动微分方程组:
gsin
（t）-
（5）
（6）
仍然使用欧拉算法求解：
将dv（t）=v（t+dt）-v（t）和d
（t）=
（t+dt）-
（t）代入
式（5）及式（6）中，并以仿真步进量Δ作为dt的近似，得到基于时间的递推方程：
v（t+Δ）=v（t）+（gsin
（t）-
）Δ
（t+Δ）=
（t）-
Δ
据此编写仿真程序：
在MALAB命令窗口输入以下命令：
subplot（2,1,1）
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;
fori=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,将两次仿真的图形放在一起,结果如下:
其它条件不变,令阻力比例系数k=0,将将两次仿真的图形放在一起,结果如下:
由上图可以看出,在理想状态下,当线速度最大的时候,角位移为0,当线速度为0时角位移最大,两个同时发生没有滞后.非理想状态下,角位移与线速度0值和最大值发生会有一定的延迟,这是由于各方面阻力综合原因造成的.
3理想条件下特殊摆动的数学模型.
这里所说的特殊摆动,指的是在理想条件下,
=
的情况下,摆球有一个很大的初始速度,可以使摆球围绕悬点不停的旋转（垂直方向,重力能够起到作用）,此时已经不能称之为单摆了,但是在理想条件下,小球的转动依然符合一定的规律,先对其进行数学物理分析,再进行MATLAB仿真.
理力学原理的分析同理想状态下的单摆分析完全一样,只是单摆是来回摆动,而此时是围绕悬点来回旋转,摆动的时候角位移是在一个范围内变化,而旋转的情况下角位移是不停的增加的,如果时间无限长,则角位移无线长,图象显示就是一条上升的曲线,没有很直观的感受,因此我们只考虑线速度的变化规律.
基于时间的递推方程依然如下:
v（t+Δ）=v（t）+（gsin
（t）-
）Δ
（t+Δ）=
（t）-
Δ
只是在MATLAB仿真中不需要图象显示角位移
的图象,仅仅显示线速度v的图象就好.MATLAB仿真程序如下:
dt=0.0001;%仿真步进
T=16;%仿真时间长度
t=0:
dt:
T;%仿真计算时间序列
g=9.8;
L=1.5;
m=8;
k=0;%空气阻力比例系数
th0=pi;%初始摆角设置
v0=10;%初始摆速设置
v=zeros（size（t））;%程序存储变量预先初始化，可提高执行速度
th=zeros（size（t））;
v
（1）=v0;
th
（1）=th0;
fori=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
plot（t,v（1:
length（t）））
xlabel（'时间t/s'）;
ylabel（'线速度v（t）'）;
仿真图象如下所示:
由上图可以看出，单摆围绕悬点做圆周运动，速度在最高点的时候最低，在最低点的时候最大，现在验证最大速度值:
带入仿真初始值：
v0=10g=9.8L=1.5从而得出
与仿真结果明显一致。
由此还可以引出许多特殊情形,比如用硬制木棍代替细摆线,力学原理就发生了很大的改变,仿真的曲线也会有所不同.
4使用Simulink重新仿真理想模式下的单摆运动.
本次仿真以单变量
建立微分方程,用示波器观察
的变化情况.
由之前的分析易知道有:
单摆回复力
而单摆的切向加速度
上式为一个二阶非线性常微分方程.为了便于比较,将所有参数设置的与之前一致,也就是g=9.8,L=1.5,初始
则
Simulink模型模块如下所示:
模块搭建好后,进行参数设置,积分器Integrator为初始角速度设置为0，积分器Integrator1为初始
=1.5,三角函数sin为默认设置,两个增益
SimulationConfigurationParameters设置为起始时间0.终止时间16,采样周期auto,最大间隔0.2,最小0.1,解法为ode45（4阶5阶龙格库塔算法）.
全部设置好后,启动Simulink开始仿真,两个示波器的输出图象如下:
与之前MATLAB程序仿真结果完全一致.
至此,本次针对单摆运动的仿真全部结束.明显可以看出Simulink在系统结构明确的情况下,只需进行简单的参数设置,就可以观测任意环节的变量状态,十分简便.
关于本次仿真的一些感受：
1.仿真的条件要求很严格，不同条件（理想非理想状态，初始状态等等）下建立的模型不一致，需要严格界定。
2.仿真的关键在于充分了解系统的物理化学等特性，并在此基础上建立系统的数学模型。
3.仿真的结果需要认真验证，对于一些明显的结论，需要能够在仿真曲线上看的清楚。
4.MATLAB作为一种辅助工具，功能十分强大，但是使用起来很简单，需要深入学习。
在系统结构参数明确情况下,使用Simulink仿真更简便一些.
5.对于连续系统,其图象数据也是用极小的采样间隔得到的,因此我理解所有的连续系统都可以用微小间隔的离散采样进行仿真.
升级会员 