机械动态仿真实验指导Word格式.docx

机械动态仿真实验指导Word格式.docx

《机械动态仿真实验指导Word格式.docx》由会员分享，可在线阅读，更多相关《机械动态仿真实验指导Word格式.docx（24页珍藏版）》请在冰豆网上搜索。

在创建模型时，如果需要一个子系统也可直接在子系统窗口创建。

这样可以省去了上面压缩子系统和重新安排窗口的步骤。

要使用子系统模块创建新的子系统，先从Function&

Table模块库中拖一个子系统模块到模型窗口中。

双击子系统模块就会出现一个子系统编辑窗口。

单个小车的运动方程如下：

注意符号的写法。

先建立单个小车的系统的子系统。

使用子系统模块创建如图所示的子系统。

此子系统用来模拟一个小车的运动。

子系统的输入为小车的左距xn-1和右距xn+1，输出为小车的当前位置xn。

子系统完成后，关闭子系统窗口。

复制两次此子系统模块，并如图所示连接起来。

为了可以对每个小车的参数赋值，要做以下设置：

建立文件set_k_m.m

%setthespringconstantsandblockmassvalues

k1=1;

k2=2;

k3=4;

k0=0;

m1=1;

m2=3;

m3=2;

仿真开始之前，在MATLAB命令窗口中运行此M文件。

然后，指定示波器模块把显示数据保存到工作空间中。

并设置仿真的起止时间为（0－100）。

仿真结束后，在MATLAB命令窗口把所得到的小车3的显示数据绘制成图。

实验2-1通过运动学仿真求解曲柄滑块的速度

一、曲柄滑块的数学模型，采用闭环矢量方程

y

（3-1）

将此矢量方程分解到x和y坐标轴上，得到

（3-2）

（3-3）

将上边两式对时间求导，有

（3-4）

（3-5）

其中

是矢量R1大小的变化率，也是滑块相对于地面的平移速度。

若曲柄的角速度已知，则方程改写为

（3-6）

二、曲柄滑块运动学的Simulink仿真

将2视为仿真的输入，可以用数值积分从速度中计算出2、3和r1。

显然，匀速时2为常量，作为仿真的输入。

选择constant模块。

另外两个速度可以从闭环矢量方程中求得。

实现这个过程需要三个积分模块作为仿真的开始。

第一个输入是2，这里设为常量并作为仿真的输入，选择constant模块。

为此编写一个MATLAB函数来求解该方程。

函数名为compvel.m

内容如下：

function[x]=compvel（u）

r2=1.0;

r3=4.0;

a=[r3*sin（u（3））1;

-r3*cos（u（3））0];

b=[-r2*u

（1）*sin（u

（2））;

r2*u

（1）*cos（u

（2））];

x=inv（a）*b;

三、参数设置（建立初始条件）

初始条件必须是机构在某个真实位置的角度和长度。

假定曲柄的初始位置2=0，此时连杆与曲柄处于同一条线上。

显然，3=0，滑块距原点的距离为r2+r3的位置，即r1=5.0（in）。

在MATLAB命令窗口键入：

>

theta-2=0;

theta-3=0;

r1=5.0;

为了显示滑块的位移r1随时间的变化图形，可以用下列MATLAB命令：

plot（simout（:

1）,simout（:

6））

xlabel（'

Time（sec）'

ylabel（'

PistonDisplacment（in）'

类似地，连杆的速度变化可以用下列命令画出：

3））

ConnectinRodSpeed（rad/s）'

实验2-2通过运动学仿真求解曲柄滑块的加速度

速度方程是对闭环矢量方程求关于时间的一阶导数而得到的。

加速度方程要对时间求二阶导数，这时要特别注意对时间变量乘积的求导过程。

闭环矢量方程对时间的二阶导数方程由下式给出：

（4-7）

（4-8）

在仿真时，位移（2、3和r1）被视为已知量（因为它们是积分的结果）。

类似地，如果将仿真扩展并计入加速度，则速度（2、3和r1）也被视为已知量。

（4-9）

比较加速度方程和速度方程左端矩阵，可以看出左端22矩阵完全一样。

加速度仿真的大部分模块结构和前面集合的速度仿真结构一样，可以修改前面的模型。

求解上面的加速度方程（4-9）需要一个新的函数文件compacc.m。

该函数需要额外的输入2、（仿真输入）2、3。

Mux模块被扩展，仿真系统增加了三个新的积分模块。

修改后的仿真系统如图。

输出矢量也相应地要加以改变，变量simout现在包括6列数据：

t、2、3、3、r’1和r1。

仿真的初始条件theta-2=0;

omega-2=0;

omega-3=0;

r1-dot=0

输入加速度为10rad/s2;

仿真运行时间为4秒。

当4秒结束时，曲柄转速应为40rad/s，曲柄共旋转80rad（大约13圈）

下图给出了连杆的角速度所经历的加速过程仿真曲线。

注意到图中给出的周期数和所确定的曲柄圈数13相吻合。

加速度不是按正弦规律变化。

纵坐标是角速度omega-3

实验3-1四杆机构的位置问题

杆长r1=12,r2=4,r3=10,r4=7

在机构分析过程中，首先要进行位置分析。

就单自由度机构而言，需要回答以下问题：

若已知机构中某一根连杆的位置（相对于地面），那么机构中其他连杆的位置应如何确定？

例如：

若给定连杆2的转角2和所有连杆的长度，则3和4可完全由方程求出。

然而这组方程是关于3和4的非线性超越方程，非常难以求解。

这里介绍非线性方程的数值解法。

采用牛顿－辛普森方法是求解非线性方程的一种迭代法，它从某一给定的初始向量开始不断地给增量直到所得结果足够接近精确解。

在前面推导得到了四连杆机构的闭环矢量方程以及其在x方向和y方向的投影方程，

（2-2）

（2-3）

将方程重新组合得到如下方程：

因此位置问题可表述为：

对给定的一组连杆长度和它的，寻求适当的，使得函数f1和f2等于零。

由于f1和f2是非线性和超越的（超越函数中包含未知量），故线性的矩阵分析方法已不适用，牛顿－辛普森方法正适合于求解这类非线性方程。

首先，以名义解的形式重新定义变量，认为名义解接近精确解，其间差值由以下修正因子描述：

其中：

3和4代表问题的解；

和

为接近解的名义解；

3和4即为修正因子。

运用上述泰勒级数，将结果表达为如下矩阵方程：

四连杆机构位置问题的MATLAB求解

求位置的方程

function[th3,th4]=possol4（th,rs）

th2=th

（1）;

th3bar=th

（2）;

%（guess）

th4bar=th（3）;

epsilon=1.0e-6;

f=[rs（3）*cos（th3bar）-rs（4）*cos（th4bar）+rs

（2）*cos（th2）-rs

（1）;

rs（3）*sin（th3bar）-rs（4）*sin（th4bar）+rs

（2）*sin（th2）];

whilenorm（f）>

epsilon

J=[-rs（3）*sin（th3bar）rs（4）*sin（th4bar）;

rs（3）*cos（th3bar）-rs（4）*cos（th4bar）];

dth=inv（J）*（-1.0*f）;

th3bar=th3bar+dth

（1）;

th4bar=th4bar+dth

（2）;

f=[rs（3）*cos（th3bar）-rs（4）*cos（th4bar）+rs

（2）*cos（th2）-rs

（1）;

norm（f）;

end;

th3=th3bar;

th4=th4bar;

输入初始条件：

rs=[12,4,10,7]

th=[0,45*pi/180,135*pi/180]

观察输出结果，正确输出方法：

[th3,th4]=possol4（th,rs）

上面求解的特点是只要给定连杆的长度，输入角度以及初始估计，即可得到四连杆机构中其它两个连杆的角度。

fourbar.m

D2R=pi/180.0;

rs=[12.0,4.0,10.0,7.0];

th=[0.0,45*D2R,100*D2R];

dth=5*D2R;

fori=1:

72

[ths

（1）,ths

（2）]=possol4（th,rs）

angles（i,:

）=[th

（1）/D2Rths

（1）/D2Rths

（2）/D2R];

th

（1）=th

（1）+dth;

th

（2）=ths

（1）;

th（3）=ths

（2）;

end

脚本文件将结果存储于名为angles的矩阵中（72行3列），则可由plot命令绘图。

plot（angles（:

1）,angles（:

2）,angles（:

axis（[03600160]）

text（50,35,'

Theta-3'

text（110,110,'

Theta-4'

ylabel（'

Solvedangles（degrees）'

xlabel（'

Inputangles（degrees）'

c=[-rs（3）*sin（th3）rs（4）*sin（th4）;

rs（3）*cos（th3）-rs（4）*cos（th4）]

c=-69.526969.5269

71.87508.1250

d=[-250*rs

（2）*sin（0）;

250*rs

（2）*cos（0）]

d=0

10000

inv（c）*d

ans=125

125

这些初始条件是很容易被实现的。

若曲柄以250rad/s的转速旋转，机构在0.025s之内旋转一周。

如果仿真运行0.1s，则大致模拟4圈。

实验3-2有曲柄条件的仿真求解。

四杆机构中有的连架杆能作整周回转而成为曲柄，有的则不能。

在机械原理的教科书中，导出铰链四杆机构有曲柄的条件为：

连架杆和机架中必有一杆是最短杆。

1．最短杆与最长杆的长度之和小于其它两杆长度之和。

上述两个条件必须同时满足，否则机构中就不存在曲柄，只能是双摇杆机构。

这里假定铰链四杆机构中存在曲柄，进行可动性分析所关心的问题是：

机构输入构件与输出构件之间起连接作用的连杆长度所许可的变化。

对于一系列输入曲柄的位置来说，最关心的是能够装配的极小和极大连杆长度。

分析：

对于问题

（1），可以用解析法求得以输入和输出角度位置的函数表示的连杆长度的极小或极大值。

用上面的具有典型性的四杆机构尺寸，曲柄为最短杆，摇杆小于机架长度：

r1=12,r2=4,r4=7，求连杆r3的极小和极大值。

建立一个MATLAB命令的脚本文件，起名为ganminmax.m，

%calculatetheconditionoffourbarhavingcrank

rs=[12.0,4.0,7.0];

th=0.0;

c=（rs

（2）*sin（th））^2+（rs

（1）-rs

（2）*cos（th））^2;

a=c^0.5-rs（3）;

b=c^0.5+rs（3）;

）=[th/D2Rab];

[angles（37,2）,angles（1,3）]

axis（[0360025]）

在MATLAB窗口运行该命令得到下图和

ganminmax

ans=

915

这两个值代表杆r2由0转过360度，能够保证有曲柄的连杆取值范围。

采用同样的方法能够确定，r4从0到360度时，容许连杆取值的范围。

脚本文件取名为lminmax.m

c=（rs（3）*sin（th））^2+（rs

（1）+rs（3）*cos（th））^2;

a=c^0.5-rs

（2）;

b=c^0.5+rs

（2）;

仅有个别参数需要修改，总体上和前面的程序大致相同。

实验4四杆机构速度和加速度仿真

为进一步说明运动学仿真的概念，将对曲柄滑块机构的仿真进行修改，用来对四连杆机构进行运动仿真。

这里只需对Simulink模型作少许修改即可完成。

两个函数compacc.m和conchk.m必须作修改以反映适当的闭环矢量方程。

写出四连杆的闭环矢量方程：

（4-1）

投影到x轴和y轴：

（4-2）

（4-3）

对时间求一次导数

（4-4）

（4-5）

假定连杆2与一台电机相连，该电机能够提供足够大的驱动力矩使得2相对保持常量。

在此条件下，2就称为机构输入，重写（4-4）和（4-5）

（4-6）

在传统的动力学分析中，确定角速度在某个时刻的大小，是在求解位置之后进行的。

也就是说，在某个时刻所有连杆转角的角度是已知的。

方程可写成矩阵的形式：

再次对时间求导并整理后：

（4-11）

*********************************************************

利用式（4-8）求解四连杆机构在图示位置时连杆3、4的角速度。

设输入连杆2的角速度为100rad/s。

各连杆长度及在图示位置的转角如下：

r1=12.0;

r2=4.0;

r3=10.0;

r4=7.0;

th2=45*pi/180;

th3=24.652*pi/180;

th4=90.6794*pi/180;

j=[-r3*sin（th3）r4*sin（th4）;

r3*cos（th3）-r4*cos（th4）]

b=[100*r2*sin（th2）;

-100*r2*cos（th2）]

om34=inv（j）*b

**************************************************************

可见速度和加速度方程左端的2×

2矩阵完全相同。

function[x]=compacc4（u）

r1=12;

r2=4;

r3=10.0;

r4=7.0;

a=[-r3*sin（u（6））r4*sin（u（6））;

r3*cos（u（6））-r4*sin（u（7））];

b

（1）=r2*u

（1）*sin（u（5））+r2*u

（2）^2*cos（u（5））+r3*u（3）^2*cos（u（6））-r4*u（4）^2*cos（u（7））;

b

（2）=-r2*u

（1）*cos（u（5））+r2*u

（2）^2*sin（u（5））+r3*u（3）^2*sin（u（6））-r4*u（4）^2*sin（u（7））;

x=inv（a）*b'

;

四连杆仿真框图只需作少许修改，具体如图：

实验5-1平面铰接二连杆的运动仿真（基于Simulink的建模方法）

本章讨论一种在机器人学文献中经常遇到的开链式机构的动力学仿真。

平面连杆机构是一种简单的两自由度的机械装置，具有一定的复杂动力特性。

下面导出这种机构的矢量方程和动力学方程。

最终的动力学仿真要求具有两个输入参数，即由两个电动机所产生的扭矩。

矢量方程

平面两连杆机器人的矢量方程为：

（5-1）

注意到方程和我们前面所有遇到的方程形式有所不同，原因在于矢量角的定义方法不同。

所给出的角度是相对于前一个连杆的方位，而不是相对于整体坐标系的x轴。

这种习惯的表示方法在机器人学中很普遍，其源于安装在机器人手臂上的传感器所测得的是连杆的相对转角，而不是绝对转角。

将矢量方程投影到x轴和y轴上

（5-2）

（5-3）

对上式求导数有

（5-4）

（5-5）

令C12=cos（1+2），S12=sin（1+2），将上两式重新整理并写成如下矩阵形式

（5-6）

这个矩阵是人们所熟知的雅可比矩阵，它是一种有效的控制算法的基础。

对矩阵方程再求导，得到加速度方程

（5-7）

上述推导出的方程构成了进行动力学仿真的基础，它们表明了有效负荷的加速度与两节点处电动机的角速度和角加速度之间的关系。

确定两连杆质心处的加速度和节点变量之间的关系是十分重要的，而这些关系推导如下

（5-9）

（5-10）

（5-11）

（5-12）

动力学方程

按照前面所讲的方法分别取两个连杆的隔离体进行分析。

从第一根连杆可以推导出下面三个运动方程

（5-13）

（5-14）

（5-15）

同样，第二根连杆的受力如图，有下述三个运动方程

（5-16）

（5-17）

（5-18）

最后研究作用在机器人上的有效载荷。

注意到有效载荷的变化是直接与第二根连杆相联系的，所以可以与连杆集中在一起考虑。

由于有效载荷的质量是随着机器人举起的不同物体而变化的，并且我们关心的是维持有效载荷所需要的力的大小，因而，再增加几个方程对研究问题更加方便。

将加载物视为集中质量，得到两个方程

（5-19）

（5-20）

综上所述，从矢量方程中得到六个标量方程，从动力学方程中得到八个标量方程，如果再将电动机的扭矩考虑进来，就有下面的14个未知量

联立约束方程

联立六个运动方程和八个动力学方程，即可得到一个关于14个未知量的线性方程组，如下

这个矩阵方程嵌入到MATLAB函数中，随后也可以嵌入到Simulink仿真中。

上图为两连杆机器人的Simulink仿真模型，未连接的输出代表冗余的加速度和约束力。

为了在仿真中获得较高的可信度，需要做一个简单的实验。

回想一下机器人在垂直平面工作时在重力作用下的受力图。

因此，如果让机器人从任何初始位置开始运动，将输入的扭矩值设置为０，那么机器人将在自身重力的作用下下落，最后达到两个连长都在一条铅垂线上的位置。

参照图7-1中节点角的定义，此时相应的节点角度是1＝－/2和2＝0rad。

对于初始条件，选择1＝0和2＝/2rad。

这时对应于加载物的位置（就是机器人操纵装置的末点位置）xpl=1.0和ypl=0.8。

如同所有的仿真一样，积分求解器的初始条件必须是相容的。

求解矩阵方程的函数程序如下：

functionout=robot（u）

%

%u

（1）=omege-1

%u

（2）=theta-1

%u（3）=omege-2

%u（4）=theta-2

%u（5）=Torque-1

%u（6）=Torque-2

g=9.8067;

r1=1.0;

rc1=0.5;

r2=0.8;

rc2=0.1;

m1=2.5;

m2=1.8;

I1=0.15;

I2=0.05;

mpl=2.0;

S1=sin（u

（2））;

S12=sin（u

（2）+u（4））;

C1=cos（u

（2））;

C12=cos（u

（2）+u（4））;

a=zeros（14,14）;

b=zeros（14,1）;

a（1,1）=r1*S1+r2*S12;

a（1,2）=r2*S12;

a（1,7）=1;

a（2,1）=-r1*C1-r2*C12;

a（2,2）=-r2*C12;

a（2,8）=1;

a（3,1）=rc1*S1;

a（3,3）=1;

a（4,1）=-rc1*C1;

a（4,4）=1;

a（5,1）=r1*S1+rc2*S12;

a（5,2）=rc2*S12;

a（5,5）=1;

a（6,1）=-r1*C1-rc2*C12;

a（6,2）=-rc2*C12;

a（6,6）=1;

a（7,3）=-m1;

a（7,9）=1;

a（7,11）=1;

a（8,4）=-m1;

a（8,10）=1;

a（8,12）=1;

a（9,1）=I1;

a（9,11）=r1*S1;

a（9,12）=-r1*C1;

a（10,5）=-m2;

a（10,11）=-1;

a（10,13）=1;

a（11,6）=-m2;

a（11,12）