《机械动态仿真》实验指导书一．课程简介机械动态仿真是为机械工程等相关专业开设的一门技术基础课，主要介绍基于MATLAB软件平台的机械运动仿真模块SimMechanics和Simulink的应用，目的在于培养学生了解和掌握机械动态仿真的基本理论、仿真方法及技术，为将来从事机械产品动态特性分析奠定基础。
二．课程实验目的实验目的：熟悉和掌握MATLAB软件的函数编程，掌握基于MA TLAB软件平台的机械仿真模块SimMechanics和Simulink的建模方法、仿真技术。
通过实验教学，使学生学会把实际问题抽象为数学模型、进而建立仿真模型的方法，增强分析和解决实际问题的能力。
三．实验方式与基本要求实验方式：学生独立在计算机上进行MATLAB软件编程，建立四连杆机构、曲柄滑块机构、二连杆平面机器人等机构的数学模型，进而建立基于Simulink和SimMechanics 动态仿真模型，设置仿真运行参数，观察不同参数下的仿真结果。
基本要求：1．为保证实验顺利进行，要求学生掌握MATLAB语言的基础知识；2．由指导教师讲清实验的基本要求、目的、建模要点及注意事项；3．要求学生严格遵守实验课守则，认真实验，按时完成实验报告。
四．实验报告每个实验均应撰写实验报告。
实验报告内容应包括：实验名称、目的、内容、实验结果分析等。
鼓励学生以小论文形式整理和撰写实验报告。
五．考试（考核）方法与规定本课程的实验考核采取理论测试与操作考核相结合方式，并考虑平时实验及实验报告完成情况，按百分制综合平定成绩。
实验1-1小车弹簧系统的运动仿真参考书图6-17 质量-阻尼-弹簧机械系统的仿真。
应用牛顿定律可以用数学模型描述：式中，y 是质量m 的位移；u 为作用质量m 上的力。
采用现代控制理论的状态方程，可以将单自由度的二阶微分方程写成两个一阶微分方程： 设 x1=y, x2=dy/dx, 则有⎪⎭⎪⎬⎫+--==m u x m f x m k dt dx x dt dx 21221 （1-2）根据（1-2）的数学模型即可编制仿真程序，之后即可上机调试运行。
步骤：1．在MATLAB 快捷菜单栏中，打开仿真库，弹出Simulink Library Browser 窗口。
2．在Simulink Library Browser 中点击File->New->Model ，以创建仿真模型。
（注意与MATLAB 命令窗口建立的模型不同。
）3．将输入信号f 视为正弦波，根据数学模型创建Simulink 仿真模型如下：4．在MA TLAB 命令窗口对质量m 、刚度k 、阻尼系数B 赋值后，即可对模型仿真。
实验1-2子系统模块在创建模型时，如果需要一个子系统也可直接在子系统窗口创建。
这样可以省去了上面压缩子系统和重新安排窗口的步骤。
要使用子系统模块创建新的子系统，先从Function & Table 模块库中拖一个子系统模块到模型窗口中。
双击子系统模块就会出现一个子系统编辑窗口。
单个小车的运动方程如下：)]()([1111n n n n n n nn x x k x x k m x -+-=++- 注意符号的写法。
先建立单个小车的系统的子系统。
使用子系统模块创建如图所示的子系统。
此子系统用来模拟一个小车的运动。
子系统的输入为小车的左距x n-1和右距x n+1，输出为小车的当前位置x n 。
子系统完成后，关闭子系统窗口。
复制两次此子系统模块，并如图所示连接起来。
为了可以对每个小车的参数赋值，要做以下设置： 建立文件set_k_m.m% set the spring constants and block mass values k1=1;k2=2;k3=4;k0=0; m1=1;m2=3;m3=2;仿真开始之前，在MATLAB 命令窗口中运行此M 文件。
然后，指定示波器模块把显示数据保存到工作空间中。
并设置仿真的起止时间为（0－100）。
仿真结束后，在MA TLAB 命令窗口把所得到的小车3的显示数据绘制成图。
实验2-1 通过运动学仿真求解曲柄滑块的速度x3一、 曲柄滑块的数学模型，采用闭环矢量方程132R R R =+ （3-1）将此矢量方程分解到x 和y 坐标轴上，得到13322r r r =+θθcos cos （3-2）03322=+θθsin sin r r （3-3）将上边两式对时间求导，有1333222rr r =--θωθωsin sin （3-4） 0333222=+θωθωcos cos r r （3-5）其中1r是矢量R 1大小的变化率，也是滑块相对于地面的平移速度。
若曲柄的角速度已知，则方程改写为⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-22222213333301θωθωωθθcos sin cos sin r r r r r （3-6） 二、曲柄滑块运动学的Simulink 仿真将ω2视为仿真的输入，可以用数值积分从速度中计算出θ2、θ3和r 1。
显然，匀速时ω2为常量，作为仿真的输入。
选择constant 模块。
另外两个速度可以从闭环矢量方程中求得。
实现这个过程需要三个积分模块作为仿真的开始。
第一个输入是ω2，这里设为常量并作为仿真的输入，选择constant 模块。
另外两个速度可以从闭环矢量方程中求得。
为此编写一个MA TLAB 函数来求解该方程。
函数名为compvel.m 内容如下：function [x]=compvel(u) r2=1.0; r3=4.0;R2R3R1θ3θ2xya=[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，滑块距原点的距离为r 2+r 3的位置，即r1=5.0（in ）。
在MA TLAB 命令窗口键入： >>theta-2=0; theta-3=0; r1=5.0;为了显示滑块的位移r 1随时间的变化图形，可以用下列MATLAB 命令： plot(simout(:,1),simout(:,6)) xlabel('Time(sec) ')ylabel('Piston Displacment (in)')类似地，连杆的速度变化可以用下列命令画出： plot(simout(:,1),simout(:,3)) xlabel('Time(sec) ')ylabel('Connectin Rod Speed (rad/s) ')实验2-2通过运动学仿真求解曲柄滑块的加速度速度方程是对闭环矢量方程求关于时间的一阶导数而得到的。
加速度方程要对时间求二阶导数，这时要特别注意对时间变量乘积的求导过程。
闭环矢量方程对时间的二阶导数方程由下式给出：132333332222222r r r r r =----θωθωθωθωcos sin cos sin （4-7）032333332222222=-+-θωθωθωθωsin cos sin cos r r r r （4-8）在仿真时，位移（θ2、θ3和r 1）被视为已知量（因为它们是积分的结果）。
类似地，如果将仿真扩展并计入加速度，则速度（ω2、ω3和r '1）也被视为已知量。
⎥⎦⎤⎢⎣⎡--++-=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-323322222223233222222213333301θωθωθαθωθωθααθθsin sin cos )cos cos sin (cos sin r r r r r r r r r （4-9） 比较加速度方程和速度方程左端矩阵，可以看出左端2⨯2矩阵完全一样。
加速度仿真的大部分模块结构和前面集合的速度仿真结构一样，可以修改前面的模型。
求解上面的加速度方程（4-9）需要一个新的函数文件compacc.m 。
该函数需要额外的输入α2、（仿真输入）ω2、ω3。
Mux 模块被扩展，仿真系统增加了三个新的积分模块。
修改后的仿真系统如图。
输出矢量也相应地要加以改变，变量simout 现在包括6列数据：t 、θ2、ω3、θ3、r’1和r1。
仿真的初始条件theta-2=0; theta-3=0; r1=5.0; omega-2=0; omega-3=0; r1-dot=0输入加速度为10rad/s 2; 仿真运行时间为4秒。
当4秒结束时，曲柄转速应为40 rad/s ，曲柄共旋转80 rad （大约13圈）下图给出了连杆的角速度所经历的加速过程仿真曲线。
注意到图中给出的周期数和所确定的曲柄圈数13相吻合。
加速度不是按正弦规律变化。
纵坐标是角速度omega-3实验3-1四杆机构的位置问题杆长r1=12, r2=4, r3=10, r4=7在机构分析过程中，首先要进行位置分析。
就单自由度机构而言，需要回答以下问题：若已知机构中某一根连杆的位置（相对于地面），那么机构中其他连杆的位置应如何确定？例如：若给定连杆2的转角θ2和所有连杆的长度，则θ3和θ4可完全由方程求出。
然而这组方程是关于θ3和θ4的非线性超越方程，非常难以求解。
这里介绍非线性方程的数值解法。
采用牛顿－辛普森方法是求解非线性方程的一种迭代法，它从某一给定的初始向量开始不断地给增量直到所得结果足够接近精确解。
在前面推导得到了四连杆机构的闭环矢量方程以及其在x 方向和y 方向的投影方程，4413322cos cos cos θθθr r r r +=+ （2-2） 443322sin sin sin θθθr r r =+ （2-3）将方程重新组合得到如下方程：因此位置问题可表述为：对给定的一组连杆长度和它的，寻求适当的，使得函数f1和f2等于零。
由于f1和f2是非线性和超越的（超越函数中包含未知量），故线性的矩阵分析方法已不适用，牛顿－辛普森方法正适合于求解这类非线性方程。
首先，以名义解的形式重新定义变量，认为名义解接近精确解，其间差值由以下修正因子描述：其中：θ3和θ4代表问题的解；3θ和4θ为接近解的名义解；∆θ3和∆θ4即为修正因子。
运用上述泰勒级数，将结果表达为如下矩阵方程：四连杆机构位置问题的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)>epsilonJ=[-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);rs(3)*sin(th3bar)-rs(4)*sin(th4ba r)+rs(2)*sin(th2)]; 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) 上面求解的特点是只要给定连杆的长度，输入角度以及初始估计，即可得到四连杆机构中其它两个连杆的角度。