北航机电系统仿真实验报告Word文档下载推荐.doc
《北航机电系统仿真实验报告Word文档下载推荐.doc》由会员分享，可在线阅读，更多相关《北航机电系统仿真实验报告Word文档下载推荐.doc（12页珍藏版）》请在冰点文库上搜索。
m
l
-
图1移动高架吊车模型
1、具有两个自由度的移动高架吊车模型
为此，使用如下定义的拉格朗日方程
其中：
q x（t）和θ（t）的自由度
D 由于摩擦而消耗的能量
Fq 由自由度q产生的力
Ec和Eq 系统的动能和势能
1．1系统移动时的动能
组件速度由下式决定：
于是
1．2系统的势能
1．3在q（t）=θ（t）自由度下的Lagrange方程
如果D=0，则自由度相等，因此不考虑总摩擦损失
根据自由度θ施加的力为零
简化后给出第一Lagrange方程
1．4在q（t）=x（t）自由度下的Lagrange方程
根据自由度x给吊车施力：
简化后给出Lagrange第二方程
1．5操作点附近的线性模型
得到的模型是非线性的，且不同自由度下有不同的方程。
如果只考虑在操作点θ0=0附近只有很小的θ变化，可考虑如下的简化
另外，有
这就给出了如下的线性微分方程
根据x（t）和θ（t）的微分系统给出
2、系统的数学模型
对x（t）和θ（t）的微分系统方程组进行拉氏变换得
将M=10kg，m=5kg，l=1m，g=9.8≈10，F0=1代入上式并进一步变换得
3、利用Simulink对系统进行仿真
3.1系统开环过程仿真
打开一空白模型编辑窗口，根据该系统的数学模型，创建、设置、连接模块，如图2所示。
图2仿真模型结构图
由于在模型文件中使用了Tofile模块，输出结果已经保存到MATLAB工作空间当中，可以使用MATLAB命令来绘制结果，代码如下，结果如图4所示。
●>
>
load（'
D:
\MATLAB7\out1'
）
plot（out（1,:
）,out（2,:
））
gridon
title（'
X（t）-t'
）;
xlabel（'
Time（sec）'
ylabel（'
Distance（m）'
◆>
\MATLAB7\out2'
θ（t）-t'
时间（秒）'
角度（弧度）'
运行仿真，结果如图3所示。
图3Simulink窗口Scope显示状态仿真结果图
图4MATLAB命令窗口绘制状态仿真结果图
比较图3和图4可以看出，结果一模一样。
3.2吊车停止过程仿真
为让吊车停止，最好在θ（t）经过零的瞬间取消力F（t）。
否则，悬挂物将继续振荡。
为此，在块振荡的2个周期内在移动高架吊车上施加1N的力，在θ（t）经过零点的时候取消力F（t），避免了残余振荡。
因此，对式进行拉氏反变换得，进一步得T==1.622，即将模块Step1的steptime设置为2T=3.245。
根据系统的数学模型和以上推导，搭建仿真模型图（如图5所示）。
图5系统仿真模型图
设定参数，运行仿真，结果如图6所示。
图6系统仿真结果
二、Maple自主实验—滑块摆
滑块摆由一置于光滑杆上的质量为m的滑块A和一质量为M的小球B和长度为L，质量不计的刚性杆铰接而成，如图1所示。
ψ
A
B
L
图1滑块摆模型
1.系统的数学模型
为此，使用如下定义的拉格朗日方程：
Ec和Ep 系统的动能和势能
不计各处摩擦，系统有两个自由度，以x和为广义坐标，以过A点的水平面为零势能面，系统的动能和势能分别为：
系统的Lagrange方程为：
计算出诸导数
带入Lagrange方程，得到系统的运动微分方程：
设定初始条件为：
m=1Kg，M=1Kg，g=9.8，L=2m
φ（0）=0rad，x（0）=0m，φ’（0）=-1.3rad/s，x’（0）=1m/s
2.利用Maple对系统进行仿真
在交互窗口中输入如下的程序语句，并“运行整个工作表”，即可播放（画）出所需的图像。
restart;
with（DEtools）:
#用微分方程求解包
with（plots）:
#用图形软件包
with（plottools）:
#用图形工具包
m:
=1:
M:
g:
=9.8:
l:
=2:
eq1:
=（m+M）*diff（x（t）,t$2）+M*l*diff（phi（t）,t$2）*cos（phi（t））-M*l*（diff（phi（t）,t））^2*sin（
phi（t））=0;
eq2:
=M*（l^2）*diff（phi（t）,t$2）+M*l*diff（x（t）,`$`（t,2））*cos（phi（t））+M*g*l*sin（phi（t））=0;
sys:
={eq1,eq2};
Ini:
={phi（0）=0,x（0）=0,D（phi）（0）=-1.3,D（x）（0）=1};
var:
={phi（t）,x（t）}:
val:
=array（1..100）:
forito100doval[i]:
=i/10enddo:
S:
=dsolve（`union`（sys,Ini）,var,type=numeric,method=rkf45,output=val）:
eval（S）:
forito100do
pos[i]:
=S[2,1][i,4];
ang[i]:
=S[2,1][i,2];
posY[i]:
=-cos（ang[i]）*l;
posX[i]:
=sin（ang[i]）*l;
enddo:
minx:
=10000:
maxx:
=-10000:
forito100do
temmin:
=`if`（minx<
pos[i],minx,pos[i]）:
minx:
=temmin:
temmax:
=`if`（maxx>
pos[i],maxx,pos[i]）:
=temmax:
fortto100do:
wall:
=curve（[[maxx,0],[minx,0]]）:
pospendx:
=posX[t]+pos[t]:
pendulum:
=disk（[pospendx,posY[t]],0.15,color=tan）:
mass:
=rectangle（[pos[t]-0.18,0.1],[pos[t]+0.18,-0.1],color=violet）:
#mass:
=disk（[pos[t],0],0.15,color=tan）:
lineM_P:
=curve（[[pos[t],0],[pospendx,posY[t]]]）:
G1[t]:
=display（pendulum,mass,lineM_P,wall）:
display（[seq（G1[t],t=1..100）],insequence=true,scaling=constrained,axes=none,title=`滑块摆的运动动画`）;
oldx:
=pos[1]:
oldcur:
=curve（[[pos[1],1],[pos[1],1]]）:
oldplot:
=display（oldcur）:
nwcur:
=curve（[[i-1,oldx],[i,pos[i]]]）:
plotcur（i）:
=display（oldplot,nwcur）:
=plotcur（i）:
oldx:
=pos[i]:
display（[seq（plotcur（t）,t=1..100）],insequence=true,title=`Phi-t`）;
oldphi:
=ang[1]:
oldcur1:
=curve（[[ang[1],1],[ang[1],1]]）:
oldplot1:
=display（oldcur1）:
nwcur1:
=curve（[[i-1,oldang1],[i,ang[i]]]）:
plotcur1（i）:
=display（oldplot1,nwcur1）:
=plotcur1（i）:
oldang1:
=ang[i]:
display（[seq（plotcur1（t）,t=1..100）],insequence=true,title=`X-t`）;
2.1滑块A的位移x随时间t的变化曲线
2.2角度φ随时间t的变化曲线
2.3滑块摆的运动动画
-11-
升级会员 