机械动态仿真实验指导

1、机械动态仿真实验指导书一课程简介机械动态仿真是为机械工程等相关专业开设的一门技术基础课，主要介绍基于MATLAB软件平台的机械运动仿真模块SimMechanics和Simulink的应用，目的在于培养学生了解和掌握机械动态仿真的基本理论、仿真方法及技术，为将来从事机械产品动态特性分析奠定基础。二课程实验目的实验目的：熟悉和掌握MATLAB软件的函数编程，掌握基于MATLAB软件平台的机械仿真模块SimMechanics和Simulink的建模方法、仿真技术。通过实验教学，使学生学会把实际问题抽象为数学模型、进而建立仿真模型的方法，增强分析和解决实际问题的能力。三实验方式与基本要求实验方式：学生

2、独立在计算机上进行MATLAB软件编程，建立四连杆机构、曲柄滑块机构、二连杆平面机器人等机构的数学模型，进而建立基于Simulink和SimMechanics动态仿真模型，设置仿真运行参数，观察不同参数下的仿真结果。基本要求：1为保证实验顺利进行，要求学生掌握MATLAB语言的基础知识；2由指导教师讲清实验的基本要求、目的、建模要点及注意事项；3要求学生严格遵守实验课守则，认真实验，按时完成实验报告。四实验报告每个实验均应撰写实验报告。实验报告内容应包括：实验名称、目的、内容、实验结果分析等。鼓励学生以小论文形式整理和撰写实验报告。五考试（考核）方法与规定本课程的实验考核采取理论测试与操作考核

3、相结合方式，并考虑平时实验及实验报告完成情况，按百分制综合平定成绩。实验1-1 小车弹簧系统的运动仿真参考书图6-17 质量-阻尼-弹簧机械系统的仿真。应用牛顿定律可以用数学模型描述：式中，y是质量m的位移；u为作用质量m上的力。采用现代控制理论的状态方程，可以将单自由度的二阶微分方程写成两个一阶微分方程：设 x1=y, x2=dy/dx, 则有 （1-2）根据（1-2）的数学模型即可编制仿真程序，之后即可上机调试运行。步骤：1在MATLAB快捷菜单栏中，打开仿真库，弹出Simulink Library Browser窗口。2在Simulink Library Browser中点击File-&

4、gt;New->Model，以创建仿真模型。（注意与MATLAB命令窗口建立的模型不同。）3将输入信号f视为正弦波，根据数学模型创建Simulink仿真模型如下：4在MATLAB命令窗口对质量m、刚度k、阻尼系数B赋值后，即可对模型仿真。 实验1-2 子系统模块在创建模型时，如果需要一个子系统也可直接在子系统窗口创建。这样可以省去了上面压缩子系统和重新安排窗口的步骤。x3要使用子系统模块创建新的子系统，先从Function & Table模块库中拖一个子系统模块到模型窗口中。双击子系统模块就会出现一个子系统编辑窗口。单个小车的运动方程如下：注意符号的写法。先建立单个小车的系统的子

5、系统。使用子系统模块创建如图所示的子系统。此子系统用来模拟一个小车的运动。子系统的输入为小车的左距xn-1和右距xn+1，输出为小车的当前位置xn。 子系统完成后，关闭子系统窗口。复制两次此子系统模块，并如图所示连接起来。 为了可以对每个小车的参数赋值，要做以下设置：建立文件set_k_m.m% set the spring constants and block mass valuesk1=1; k2=2; k3=4; k0=0;m1=1; m2=3; m3=2;仿真开始之前，在MATLAB命令窗口中运行此M文件。然后，指定示波器模块把显示数据保存到工作空间中。并设置仿真的起止时间为（010

6、0）。仿真结束后，在MATLAB命令窗口把所得到的小车3的显示数据绘制成图。实验2-1 通过运动学仿真求解曲柄滑块的速度一、 曲柄滑块的数学模型，采用闭环矢量方程yR2R3R1q3q2x （3-1）将此矢量方程分解到x和y坐标轴上，得到 （3-2） （3-3）将上边两式对时间求导，有 （3-4） （3-5） 其中是矢量R1大小的变化率，也是滑块相对于地面的平移速度。若曲柄的角速度已知，则方程改写为 （3-6）二、曲柄滑块运动学的Simulink仿真将w2视为仿真的输入，可以用数值积分从速度中计算出q2、q3和r1。显然，匀速时w2为常量，作为仿真的输入。选择constant模块。另外两个速度可

7、以从闭环矢量方程中求得。实现这个过程需要三个积分模块作为仿真的开始。第一个输入是w2，这里设为常量并作为仿真的输入，选择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; 三、参数设置（建立初始条件）初始条件必须是机构在某个真实位置的角度和长度。假定曲柄的初始位置q2=0，此时连杆与

8、曲柄处于同一条线上。显然，q3=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('Piston Displacment (in)')类似地，连杆的速度变化可以用下列命令画出：plot(simout(:,1),simout(:,3)xlabel('Time(se

9、c) ')ylabel('Connectin Rod Speed (rad/s) ')实验2-2 通过运动学仿真求解曲柄滑块的加速度速度方程是对闭环矢量方程求关于时间的一阶导数而得到的。加速度方程要对时间求二阶导数，这时要特别注意对时间变量乘积的求导过程。闭环矢量方程对时间的二阶导数方程由下式给出： （4-7） （4-8）在仿真时，位移（q2、q3和r1）被视为已知量（因为它们是积分的结果）。类似地，如果将仿真扩展并计入加速度，则速度（w2、w3和r¢1）也被视为已知量。 （4-9）比较加速度方程和速度方程左端矩阵，可以看出左端2´2矩阵完全一样。加

10、速度仿真的大部分模块结构和前面集合的速度仿真结构一样，可以修改前面的模型。求解上面的加速度方程（4-9）需要一个新的函数文件compacc.m。该函数需要额外的输入a2、（仿真输入）w2、w3。Mux模块被扩展，仿真系统增加了三个新的积分模块。修改后的仿真系统如图。输出矢量也相应地要加以改变，变量simout现在包括6列数据：t、q2、w3、q3、r1和r1。仿真的初始条件theta-2=0; theta-3=0; r1=5.0; omega-2=0; omega-3=0; r1-dot=0输入加速度为10rad/s2; 仿真运行时间为4秒。当4秒结束时，曲柄转速应为40 rad/s，曲柄共旋

11、转80 rad （大约13圈）下图给出了连杆的角速度所经历的加速过程仿真曲线。注意到图中给出的周期数和所确定的曲柄圈数13相吻合。加速度不是按正弦规律变化。纵坐标是角速度omega-3实验3-1 四杆机构的位置问题q2q3q4杆长r1=12, r2=4, r3=10, r4=7在机构分析过程中，首先要进行位置分析。就单自由度机构而言，需要回答以下问题：若已知机构中某一根连杆的位置（相对于地面），那么机构中其他连杆的位置应如何确定？例如：若给定连杆2的转角q2和所有连杆的长度，则q3和q4可完全由方程求出。然而这组方程是关于q3和q4的非线性超越方程，非常难以求解。这里介绍非线性方程的数值解法。

12、采用牛顿辛普森方法是求解非线性方程的一种迭代法，它从某一给定的初始向量开始不断地给增量直到所得结果足够接近精确解。在前面推导得到了四连杆机构的闭环矢量方程以及其在x方向和y方向的投影方程， （2-2） （2-3）将方程重新组合得到如下方程：因此位置问题可表述为：对给定的一组连杆长度和它的，寻求适当的，使得函数f1和f2等于零。由于f1和f2是非线性和超越的（超越函数中包含未知量），故线性的矩阵分析方法已不适用，牛顿辛普森方法正适合于求解这类非线性方程。首先，以名义解的形式重新定义变量，认为名义解接近精确解，其间差值由以下修正因子描述： 其中：q3和q4代表问题的解；和为接近解的名义解；Dq3和

13、Dq4即为修正因子。运用上述泰勒级数，将结果表达为如下矩阵方程：四连杆机构位置问题的MATLAB求解求位置的方程function th3, th4=possol4(th,rs)th2=th(1);th3bar=th(2); %(guess)th4bar=th(3); %(guess)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);while norm(f)>epsilon J=-rs(3)

14、*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);rs(3)*sin(th3bar)-rs(4)*sin(th4bar)+rs(2)*sin(th2); norm(f);end;th3=th3bar;th4=th4bar;输入初始条件：rs=12,4,10,7th

15、=0,45*pi/180,135*pi/180观察输出结果，正确输出方法：th3,th4=possol4(th,rs)上面求解的特点是只要给定连杆的长度，输入角度以及初始估计，即可得到四连杆机构中其它两个连杆的角度。fourbar.mD2R=pi/180.0;rs=12.0,4.0,10.0,7.0;th=0.0,45*D2R,100*D2R;dth=5*D2R;for i=1:72 ths=possol4(th,rs) angles(i,:)=th(1)/D2R ths(1)/D2R ths(2)/D2R; th(1)=th(1)+dth; th(2)=ths(1); th(3)=ths(2

16、);end脚本文件将结果存储于名为angles的矩阵中（72行3列），则可由plot命令绘图。plot(angles(:,1), angles(:,2), angles(:,1), angles(:,3)axis(0 360 0 160)>> text(50,35,'Theta-3')>> text(110,110,'Theta-4')>> ylabel('Solved angles (degrees)')>> xlabel('Input angles (degrees)') c=-

17、rs(3)*sin(th3) rs(4)*sin(th4);rs(3)*cos(th3) -rs(4)*cos(th4)c = -69.5269 69.5269 71.8750 8.1250d=-250*rs(2)*sin(0);250*rs(2)*cos(0)d = 0 10000inv(c)*dans = 125 125这些初始条件是很容易被实现的。若曲柄以250rad/s的转速旋转，机构在0.025s之内旋转一周。如果仿真运行0.1s，则大致模拟4圈。实验3-2 有曲柄条件的仿真求解。ADCBBDr1r2r3r4ACC2C2q四杆机构中有的连架杆能作整周回转而成为曲柄，有的则不能。在机械

18、原理的教科书中，导出铰链四杆机构有曲柄的条件为：连架杆和机架中必有一杆是最短杆。1 最短杆与最长杆的长度之和小于其它两杆长度之和。上述两个条件必须同时满足，否则机构中就不存在曲柄，只能是双摇杆机构。这里假定铰链四杆机构中存在曲柄，进行可动性分析所关心的问题是：机构输入构件与输出构件之间起连接作用的连杆长度所许可的变化。对于一系列输入曲柄的位置来说，最关心的是能够装配的极小和极大连杆长度。分析：对于问题（1），可以用解析法求得以输入和输出角度位置的函数表示的连杆长度的极小或极大值。用上面的具有典型性的四杆机构尺寸，曲柄为最短杆，摇杆小于机架长度：r1=12, r2=4, r4=7，求连杆r3的极

19、小和极大值。建立一个MATLAB命令的脚本文件，起名为ganminmax.m，% calculate the condition of fourbar having crankD2R=pi/180.0;rs=12.0,4.0,7.0;th=0.0;dth=5*D2R;for i=1:72 c=(rs(2)*sin(th)2+(rs(1)-rs(2)*cos(th)2; a=c0.5-rs(3); b=c0.5+rs(3); angles(i,:)=th/D2R a b; th(1)=th(1)+dth;endangles(37,2),angles(1,3)plot(angles(:,1), a

20、ngles(:,2), angles(:,1), angles(:,3)axis(0 360 0 25)在MATLAB窗口运行该命令得到下图和>> ganminmaxans =9 15这两个值代表杆r2由0转过360度，能够保证有曲柄的连杆取值范围。 采用同样的方法能够确定，r4从0到360度时，容许连杆取值的范围。 脚本文件取名为lminmax.m% calculate the condition of fourbar having crankD2R=pi/180.0;rs=12.0,4.0,7.0;th=0.0;dth=5*D2R;for i=1:72 c=(rs(3)*sin

21、(th)2+(rs(1)+rs(3)*cos(th)2; a=c0.5-rs(2); b=c0.5+rs(2); angles(i,:)=th/D2R a b; th(1)=th(1)+dth;endangles(37,2),angles(1,3)plot(angles(:,1), angles(:,2), angles(:,1), angles(:,3)axis(0 360 0 25)仅有个别参数需要修改，总体上和前面的程序大致相同。实验4 四杆机构速度和加速度仿真为进一步说明运动学仿真的概念，将对曲柄滑块机构的仿真进行修改，用来对四连杆机构进行运动仿真。这里只需对Simulink模型作少许

22、修改即可完成。两个函数compacc.m和conchk.m必须作修改以反映适当的闭环矢量方程。写出四连杆的闭环矢量方程： （4-1）投影到x轴和y轴： （4-2） （4-3）对时间求一次导数 （4-4） （4-5）假定连杆2与一台电机相连，该电机能够提供足够大的驱动力矩使得w2相对保持常量。在此条件下，w2就称为机构输入，重写（4-4）和（4-5） （4-6） （4-7）在传统的动力学分析中，确定角速度在某个时刻的大小，是在求解位置之后进行的。也就是说，在某个时刻所有连杆转角的角度是已知的。方程可写成矩阵的形式： （4-8）再次对时间求导并整理后：（4-11）*利用式（4-8）求解四连杆机构在

23、图示位置时连杆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=1

24、0.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的建模方法）本章讨论一种在机器人学文献中

25、经常遇到的开链式机构的动力学仿真。平面连杆机构是一种简单的两自由度的机械装置，具有一定的复杂动力特性。下面导出这种机构的矢量方程和动力学方程。最终的动力学仿真要求具有两个输入参数，即由两个电动机所产生的扭矩。矢量方程平面两连杆机器人的矢量方程为： （5-1）注意到方程和我们前面所有遇到的方程形式有所不同，原因在于矢量角的定义方法不同。所给出的角度是相对于前一个连杆的方位，而不是相对于整体坐标系的x轴。这种习惯的表示方法在机器人学中很普遍，其源于安装在机器人手臂上的传感器所测得的是连杆的相对转角，而不是绝对转角。将矢量方程投影到x轴和y轴上 （5-2） （5-3）对上式求导数有 （5-4） （5

26、-5）令C12=cos(q1+q2)，S12=sin(q1+q2)，将上两式重新整理并写成如下矩阵形式 （5-6）这个矩阵是人们所熟知的雅可比矩阵，它是一种有效的控制算法的基础。对矩阵方程再求导，得到加速度方程 （5-7）上述推导出的方程构成了进行动力学仿真的基础，它们表明了有效负荷的加速度与两节点处电动机的角速度和角加速度之间的关系。确定两连杆质心处的加速度和节点变量之间的关系是十分重要的，而这些关系推导如下 （5-9） （5-10） （5-11）（5-12） 动力学方程按照前面所讲的方法分别取两个连杆的隔离体进行分析。从第一根连杆可以推导出下面三个运动方程 （5-13） （5-14） （5

27、-15）同样，第二根连杆的受力如图，有下述三个运动方程 （5-16） （5-17） （5-18）最后研究作用在机器人上的有效载荷。注意到有效载荷的变化是直接与第二根连杆相联系的，所以可以与连杆集中在一起考虑。由于有效载荷的质量是随着机器人举起的不同物体而变化的，并且我们关心的是维持有效载荷所需要的力的大小，因而，再增加几个方程对研究问题更加方便。将加载物视为集中质量，得到两个方程 （5-19） （5-20）综上所述，从矢量方程中得到六个标量方程，从动力学方程中得到八个标量方程，如果再将电动机的扭矩考虑进来，就有下面的14个未知量 联立约束方程联立六个运动方程和八个动力学方程，即可得到一个关于1

28、4个未知量的线性方程组，如下这个矩阵方程嵌入到MATLAB函数中，随后也可以嵌入到Simulink仿真中。上图为两连杆机器人的Simulink仿真模型，未连接的输出代表冗余的加速度和约束力。为了在仿真中获得较高的可信度，需要做一个简单的实验。回想一下机器人在垂直平面工作时在重力作用下的受力图。因此，如果让机器人从任何初始位置开始运动，将输入的扭矩值设置为，那么机器人将在自身重力的作用下下落，最后达到两个连长都在一条铅垂线上的位置。参照图7-1中节点角的定义，此时相应的节点角度是q1p/2和q20 rad。对于初始条件，选择q10和q2p/2 rad。这时对应于加载物的位置（就是机器人操纵装置的

29、末点位置）xpl=1.0和ypl=0.8。如同所有的仿真一样，积分求解器的初始条件必须是相容的。求解矩阵方程的函数程序如下：function out=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