基于MATLAB六杆机构动力学分析与仿真Word下载.docx

1、如图1-4，对构件BC的约束反力推导如下，图 1-4 CD构件受力模型由运动学可以推导得，将上述BC构件，CD构件各式合并成矩阵形式有，= （1-22）如图1-5 对构件5进行约束反力的推导如下，图1-5 CE杆件受力模型如图1-6 对滑块进行受力分析如下，滑块受力模型由运动学可推，（1-23）二 编程与仿真利用MATLAB进行仿真分析，主要包括两个步骤：首先是编制计算所需要的函数模块，然后利用其仿真工具箱Simulink建立仿真系统框图，设定初始参数进行仿真分析。针对建立完成的数学模型，为了进行矩阵运算，根据以上式子编制M函数文件chengcrank.m ，chengrrr.m、chengc

2、rankdy.m、chengrrrdy.m、chengrrp.m和chengrrpdy.m如下：曲柄原动件M函数文件chengcrank.m：function y=chengcrank（x）%Function to compute the accleration of crank%Input parameters%x（1）=theta-1%x（2）=dtheta-1%x（3）=ddtheta-1%0utput parameters%y（1）=ReddB%y（2）=ImddBr1=0.4;ddB=r1*x（3）*cos（x（1）+pi/2）+r1*x（2）2*cos（x（1）+pi）;r1*x（

3、3）*sin（x（1）+pi/2）+r1*x（2）2*sin（x（1）+pi）;y=ddB;RRR II级杆组M函数文件chengrrr.m：function y=chengrrr（x）%function to compute the acceleration for RRR bar group%x（1）=theta-2%x（2）=theta-3%x（3）=dtheta-2%x（4）=dtheta-3%x（5）=ReddB%x（6）=ImddB%Output parameters%y（1）=ddtheta-2%y（2）=ddtheta-3%y（3）=ReddC%y（4）=ImddCr2=1.2

4、; r3=0.8; ReddD=0; ImddD=0;a=r2*cos（x（1）+pi/2） -r3*cos（x（2）+pi/2）; r2*sin（x（1）+pi/2） -r3*sin（x（2）+pi/2）;b=-r2*cos（x（1）+pi） r3*cos（x（2）+pi）;-r2*sin（x（1）+pi） r3*sin（x（2）+pi）*x（3）2;x（4）2+ReddD-x（5）;ImddD-x（6）;ddth=inv（a）*b;y（1）=ddth（1）;y（2）=ddth（2）;y（3）=x（5）+r2*ddth（1）*cos（x（1）+pi/2）+r2*x（3）2*cos（x（1）+

5、pi）;y（4）=x（6）+r2*ddth（1）*sin（x（1）+pi/2）+r2*x（3）2*sin（x（1）+pi）;曲柄原动件动力学M函数文件 chengcrankdy.m：function y=chengcrankdy（x）%Function for Dyanmic analysis of crank%Input parameters%x（4）=RxB%x（5）=RyB%0utput parameters%y（1）=RxA%y（2）=RyA%y（3）=M1 g=9.8; %重力加速度 %曲柄长度rc1=0.2;%质心离铰链A的距离m1=1.2;%曲柄质量J1=0.016; %绕质心转

6、动惯量Fx1=0; Fy1=0; MF=0;%作用于质心的外力和外力矩ReddA=0; ImddA=0;%铰链A的加速度y（1）=m1*ReddA+m1*rc1*x（3）*cos（x（1）+pi/2）+m1*rc1*x（2）2*cos（x（1）+pi）-Fx1+x（4）;y（2）=m1*ImddA+m1*rc1*x（3）*sin（x（1）+pi/2）+m1*rc1*x（2）2*sin（x（1）+pi）-Fy1+x（5）+m1*g;y（3）=J1*x（3）-y（1）*rc1*sin（x（1）+y（2）*rc1*cos（x（1）-x（4）*（r1-rc1）*sin（x（1）+x（5）*（r1-rc

7、1）*cos（x（1）-MF;RRR II级杆组动力学M函数文件 chengrrrdy.m：function y=chengrrrdy（x）%Function for Dyanmic analysis of RRR dayard group%x（1）=theta-2 %x（3）=dtheta-2 %x（5）=ddtheta-2 %x（6）=ddtheta-3%x（7）=ReddB %x（8）=ImddB%x（9）=Fx3 %x（10）=Fy3 %x（11）=M3%y（1）=RxB%Y（2）=RyB%y（3）=RxC%y（4）=RyC%y（5）=RxD%y（6）=RyD %两杆的长度rc2=0.

8、6;rc3=0.4; %质心到铰链B的距离 %质心到铰链D的距离 m2=3; m3=2.2; %两杆的质量J2=0.25;J3=0.09;%两杆的转动惯量ReddD=0;ImddD=0;Fx2=0; Fy2=0;M2=0; %2杆的外力和外力矩a=zeros（6）;a（1,1）=1;a（1,3）=1;a（2,2）=1;a（2,4）=1;a（3,1）=rc2*sin（x（1）;a（3,2）=-rc2*cos（x（1）;a（3,3）=-（r2-rc2）*sin（x（1）;a（3,4）=（r2-rc2）*cos（x（1）;a（4,3）=-1;a（4,5）=1;a（5,4）=-1;a（5,6）=1;a

9、（6,3）=（r3-rc3）*sin（x（2）;a（6,4）=-（r3-rc3）*cos（x（2）;a（6,5）=rc3*sin（x（2）;a（6,6）=-rc3*cos（x（2）;b=zeros（6,1）;b（1,1）=m2*rc2*x（5）*cos（x（1）+pi/2）+m2*x（7）+m2*rc2*x（3）2*cos（x（1）+pi）-Fx2;b（2,1）=m2*rc2*x（5）*sin（x（1）+pi/2）+m2*x（8）+m2*rc2*x（3）2*sin（x（1）+pi）-Fy2+m2*g;b（3,1）=J2*x（5）-M2;b（4,1）=m3*rc3*x（6）*cos（x（2）+p

10、i/2）+m3*ReddD+m3*rc3*x（4）2*cos（x（2）+pi）-x（9）;b（5,1）=m3*rc3*x（6）*sin（x（2）+pi/2）+m3*ImddD+m3*rc3*x（4）2*sin（x（2）+pi）-x（10）+m3*g;b（6,1）=J3*x（6）-x（11）;y=inv（a）*b;RRP II级杆组M函数文件：function y=chengrrp（x）%function to compute the acceleration for RRP bar group%x（1）=theta-5%x（2）=dtheta-5%x（3）=ReddC%x（4）=ImddC%x

11、（5）=ds%y（1）=ddtheta-5%y（2）=ddsr5=1.2; th6=0;a=r5*cos（x（1）+pi/2） -cos（th6）; r5*sin（x（1）+pi/2） -sin（th6）;b=-r5*cos（x（1）+pi） 0; -r5*sin（x（1）+pi） 0*x（2）2; x（5）+ReddD-x（3）; ImddD-x（4）;RRP II级杆组动力学M函数文件：function y=chengrrpdy（x）%Function for Dyanm5c analysis of RRP dayard group%x（1）=theta-5 %x（2）=dtheta-5

12、%x（3）=ddtheta-5%x（4）=dds-6 %x（5）=ReddC %x（6）=ImddC%y（1）=RxC %Y（2）=RyC%y（3）=RxE %y（4）=RyE%y（5）=RF %移动副的约束反力 %杆的长度rc5=0.6; %质心到铰链B的距离 m5=3.6; m6=6; %杆、块的质量J5=0.45;Fx5=0; Fy5=0;Fx6=1000; Fy6=0;M5=0;th6=0;a=zeros（5）;a（3,1）=rc5*sin（x（1）;a（3,2）=-rc5*cos（x（1）;a（3,3）=-（r5-rc5）*sin（x（1）;a（3,4）=（r5-rc5）*cos（x

13、（1）;a（4,5）=-sin（th6）;a（5,5）=cos（th6）;b=zeros（5,1）;b（1,1）=m5*x（5）+m5*rc5*x（3）*cos（x（1）+pi/2）+m5*rc5*x（2）2*cos（x（1）+pi）-Fx5;b（2,1）=m5*x（6）+m5*rc5*x（3）*sin（x（1）+pi/2）+m5*rc5*x（2）2*sin（x（1）+pi）-Fy5+m5*g;b（3,1）=J5*x（3）-M5;b（4,1）=m6*x（4）*cos（th6）-Fx6;b（5,1）=m6*x（4）*sin（th6）-Fx6+m6*g;三 系统仿真框图进入MATLAB，在命令栏中

14、键入Simulink进入仿真界面，根据信息传递的逻辑关系，建立仿真系统框图如图3-1. 然后设定各环节的初始参数，即可以对机构进行运动学仿真分析，再利用MATLAB的plot命令根据需要绘制曲线。图3-1四 仿真的实现再设计完成仿真框图之后，为了进行仿真还必须设定初始参数值。连杆机构杆长已经在simulink框图中给定，如果设定初始夹角为62，=10 rad/s，曲柄1作匀速转动（即），接下来要确定杆2，3的角位移和角速度，杆5的角位移和角速度，滑块的速度。可以利用辛普森方法（在MATLAB命令框中输入M函数为rrrposi）求得=0.3612rad/s，=1.8101rad/s，再利用MAT

15、LAB（在命令框输入rrrvel）求出W2=-2.2345, W3=3.3250,再利用杆3的角位移和角速度、杆5的角位移求得（在MATLAB命令框中输入M函数为compvel）W5=0.6962，ds=-3.1323。对仿真框图中各积分器设定参数变量x并在matlab命令框输入变量 x=62*pi/180 10 0.3612 1.8101 -2.2345 3.3250 -41*pi/180 0.6962 -3.1323;其中初始数值分别对应：theta-1、omega-1、theta-2、omega-2、omega-3、theta-5,omega-5 ds，以及仿真时间为1s，后进行仿真，利用MATLAB中的plot绘图命令把角速度曲线分别绘制出来。在MATLAB命令中键入：plot（tout,simout（:,1），plot（tout,simout（:,2），plot（tout,simout（:,3） ，plot（tout,simout（:,4），plot（tout,simout2（:,5），即可得到点A的水平方向、垂直方向的约束反力、驱动力矩M1及其所作功W1的变化曲线，如图所示。转动副A水平方向力 转动副A垂直方向力曲柄上作用的力矩M1 曲柄力矩所做的功W1 移动副F的约束反力