基于MATLAB六杆机构动力学分析与仿真

1、六杆机构的动力学分析仿真一 系统模型建立为了对机构进行仿真分析，首先必须建立机构数学模型，即位置方程，然后利用MATLAB 仿真分析工具箱 Simulink 对其进行仿真分析。图 3 24 所示是由原动件 （曲柄 1） 和 RRRRRP 六杆机构。 各构件的尺寸为 r1=400mm ，r2=1200mm ，r3=800mm,r4=1500mm ， r5=1200mm ；各构件的质心为 rc1=200mm,rc2 600mm,rc3=400mm ， rc5=600mm ；质量为 m1=1 2kg ， m2 3kg， m3 2 2kg ； m5=3 6kg ，m6=6kg; 转动惯量为 J1=0.

2、016kg m2， J2=0.25kgm2；J3=0.09kgm2，J5=0.45kgm2；构件 6 的工作阻力 F6=1000N ，其他构件 所受外力和外力矩均为零，构件 1 以等角速度 10 rad/s逆时针方向回转，试求不计摩擦时， 转动副 A 的约束反力、驱动力矩、移动副 F 的约束反力。此机构模型可以分为曲柄的动力学、图 1-1RRR II 级杆组的动力学和 RRP II 级杆组的动力学，再分别对这三个模型进行相应参数的求解。1 / 12图 1-2 AB 构件受力模型如上图 1-2 对于曲柄 AB由理论力学可以列出表达式： ?R XARXBFX1 m1Res1R yARyB?Fy1

3、m1 Im s1M1MFRXArc1 sin 1 RyArc1 cos 1 RXB(r1 rc1 ) sin 1 RyB(r1 rc1) cos 1 J1由运动学知识可以推得：?2Res1ReArc1 1 cos( 1/2)rc1 1 cos( 1 )?2Im s1Im Arc1 1 sin( 1/2)rc1 1 sin( 1 )将上述各式合并成矩阵形式有，RXA?ReA?/2)?2) FX1 RXBm1m1rc11 cos( 1m1rc11 cos( 1RyA?2(1-21)M1m1Im Am1rc11 sin( 1/2)m1rc11 sin( 1) Fy1 RyB m1g如图 1-3，对构

4、件 BC 的约束反力推导如下，2 / 12图 1-3 BC 构件受力模型RXBFX2RXCm2 Res2RyBFy2RyCm2g m2 ImM2RXBrc2 sin2 RyBrc2 cos如图1-4，对构件BC 的约束反力推?2?s2?RXC(r2 rc2)sin 2 RyC (r2 rc2)cos 2 J2 2FX3RXC?m3 Res3Fy3RyC?m3g m3 Im s3M3RXD rc3 sin3 RyD rc3 cos 3RXC (r3?rc3)sin 3 RyC (r3 rc3)cos 3 J3 3由运动学可以推导得，3 / 12?Im s2?Im B?rc2?2 sin( 2/

5、2)?2rc2 2 sin( 2)?2Res2ReBrc22 cos( 2/ 2)rc2 2 cos( 2)?2Res3ReDrc33 cos( 3/2)rc3 3 cos( 3)?Im s3?Im D?rc3?3 sin( 3/2)?2rc3 3 sin( 3)将上述BC 构件， CD 构件各式合并成矩阵形式有，101000RXB010100RyBrc2 sin2 rc 2 cos2 (r2rc2 ) sin 2(r2 rc2 ) cos200RXC001010RyC000101RXD00(r3rc3 ) sin 3(r3 rc3 )cos3rc3 sin 3rc3 cos 3RyD?2m2

6、ReBm2rc2 2cos( 2/2)m2rc2 2 cos( 2 )FX 2?2m2Im Bm2rc2 2sin( 2/2)m2rc2 2 sin( 2 )Fy2m2g?J2 2M2?2m3 ReD m3rc33 cos(3 /2) m3rc3 3 cos( 3)FX3?2m3ImDm3rc3 3sin( 3/2)m3rc3 3 sin( 3 )Fy3m3g?(1-22)J3 3 M 3 如图 1-5 对构件 5 进行约束反力的推导如下，图 1-5 CE 杆件受力模型4 / 12RxCFx5RxE?m5 ResRyCFy5RyE?m5g m5Im sM5RxCrc5 sin?5 RyCrc5

7、 cos 5 RxE (r5 rc5 )sin 5 RyE (r5 rc5)cos 5 J5 5如图1-6 对滑块进行受力分析如下，滑块受力模型?Fx6RxERF sin 6m6 ReE?Fy6RyERF cos 6m6g m6Im Em6 s6 cos 6 Fx 6?m6 s6 sin 6 Fy6 m6g由运动学可推，?2ResReCrc5 5cos(5/ 2) rc55 cos(5)?2Im sReCrc5 5sin(5/2) rc55 sin( 5)?ReEs6 cos6Im Es6 sin610100RxC01010R yCrc5 sin5rc5 cos 5(rc5r5 )sin5(r

8、c5r5)cos50RxE0010sin 6RyE0?0?0?21cos 6RF（1-23）m5Re Cm5rc55 cos( 5/ 2)m5rc5 5cos( 5)Fx5?2m5 ReC m5rc5 5sin(5/ 2) m5rc5 5 sin(5)Fy5m5g?J 5 5 M 5 ?5 / 12二 编程与仿真利用 MATLAB 进行仿真分析， 主要包括两个步骤： 首先是编制计算所需要的函数模块， 然后利用其仿真工具箱 Simulink 建立仿真系统框图，设定初始参数进行仿真分析。针对建 立完成的数学模型，为了进行矩阵运算，根据以上式子编制 M 函数文件 chengcrank.m ， che

9、ngrrr.m、 chengcrankdy.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

10、*x(2)2*cos(x(1)+pi);r1*x(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%Input parameters%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)=ddth

11、eta-2%y(2)=ddtheta-3%y(3)=ReddC%y(4)=ImddCr2=1.2; 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);6 / 12 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

12、);y(3)=x(5)+r2*ddth(1)*cos(x(1)+pi/2)+r2*x(3)2*cos(x(1)+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(1)=theta-1%x(2)=dtheta-1%x(3)=ddtheta-1%x(4)=RxB%x(5)=RyB%0utput

13、 parameters%y(1)=RxA%y(2)=RyA%y(3)=M1g=9.8; % 重力加速度r1=0.4; % 曲柄长度rc1=0.2;% 质心离铰链 A 的距离m1=1.2;% 曲柄质量J1=0.016; % 绕质心转动惯量Fx1=0; Fy1=0; M F=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)+m

14、1*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-rc1)*cos(x (1)-M F;RRR II 级杆组动力学 M 函数文件 chengrrrdy.m ：function y=chengrrrdy(x)%Function for Dyanmic analysis of RRR dayard group%Input parameters%x(1)=theta-2%x(2)=theta-3%x(3)=dt

15、heta-2%x(4)=dtheta-3%x(5)=ddtheta-27 / 12%x(6)=ddtheta-3%x(7)=ReddB%x(8)=ImddB%x(9)=Fx3%x(10)=Fy3%x(11)=M3%0utput parameters%y(1)=RxB%Y(2)=RyB%y(3)=RxC%y(4)=RyC%y(5)=RxD%y(6)=RyDg=9.8; % 重力加速度r2=1.2; r3=0.8; % 两杆的长度rc2=0.6;rc3=0.4; %质心到铰链 B 的距离 %质心到铰链 D 的距离 m2=3; m3=2.2; % 两杆的质量J2=0.25;J3=0.09;%两杆的转

16、动惯量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(6,3)=(r3-rc3)*sin(x(2);a(6,4)=-(r3-rc3)*cos(x(2);a(6,5)=rc3*sin(x(2);a(6

17、,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;8 / 12 b(4,1)=m3*rc3*x(6)*cos(x(2)+pi/2)+m3*ReddD+m3*rc3*x(4)2*cos(x(2)+pi)-x(9); b(5,1)=m3*rc3*x(6)

18、*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%Input parameters %x(1)=theta-5 %x(2)=dtheta-5 %x(3)=ReddC%x(4)=ImddC%x(5)=ds%Output parameters%y(1)=ddtheta-

19、5%y(2)=ddsr5=1.2; th6=0; ReddD=0; ImddD=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); y=inv(a)*b;RRP II 级杆组动力学 M 函数文件：function y=chengrrpdy(x)%Function for Dyanm5c analysis of RRP dayard group%Input paramete

20、rs %x(1)=theta-5%x(2)=dtheta-5 %x(3)=ddtheta-5%x(4)=dds-6 %x(5)=ReddC%x(6)=ImddC%0utput parameters%y(1)=RxC %Y(2)=RyC%y(3)=RxE%y(4)=RyE%y(5)=RF% 移动副的约束反力g=9.8; % 重力加速度 r5=1.2; % 杆的长度rc5=0.6; % 质心到铰链 B 的距离m5=3.6; m6=6; %杆、块的质量 J5=0.45;Fx5=0; Fy5=0;9 / 12Fx6=1000; Fy6=0;M5=0;th6=0;a=zeros(5);a(1,1)=1;

21、a(1,3)=1;a(2,2)=1;a(2,4)=1;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(1);a(4,3)=-1;a(4,5)=-sin(th6);a(5,4)=-1;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)+

22、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;y=inv(a)*b;三 系统仿真框图进入 MA TLAB ，在命令栏中键入 Simulink 进入仿真界面，根据信息传递的逻辑关系， 建立仿真系统框图如图 3-1. 然后设定各环节的初始参数，即可以对机构进行运动学仿真分 析，再利用 MA TLAB 的 plot 命令根据需要绘制曲线。10 / 12四 仿真的实现再设计完成仿真框图之后， 为了进行仿真还必须设定初始参数值。 连杆机构杆长已经在 simulink 框图中给定，如果设定 1 初始夹角为 62 ， 1=10 rad/s，曲柄 1 作匀速转动(即 1 0 )，接下来要确