基于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.016kg

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

3、1m1 Im s1M 1M FRXA rc1 sin 1RyA rc1 cos 1RXB (r1rc1 ) sin 1RyB ( r1rc1 ) cos 1J11由运动学知识可以推得：Re s1Re Arc11cos( 1Im s1Im Arc11sin( 1将上述各式合并成矩阵形式有，RXAm1Re Am1rc11 cos( 1/ 2)m1rc1RyAM 1m1Im Am1 rc11 sin( 1/ 2)m1rc1如图 1-3，对构件BC 的约束反力推导如下，/ 2)rc1/ 2)rc121 cos(121 sin( 121 cos(1)21 sin(1)FX1RXB(1-21)Fy1RyB

4、m1 g图 1-3 BC 构件受力模型RXBFX 2RXCm2 Re s 2RyBFy2RyCm2 gm2 Im s 2M 2RXB rc2 sin 2RyB rc2 cos 2RXC ( r2rc2 ) sin2RyC (r2rc 2 ) cos 2J22如图 1-4，对构件BC 的约束反力推导如下，图 1-4CD 构件受力模型RXDFX 3R XCm3 Re s 3RyDFy3RyCm3 gm3 Im s3M 3RXD rc3 sin 3RyD rc3 cos 3RXC (r3rc 3 ) sin 3RyC (r3rc 3 ) cos 3J 33由运动学可以推导得，2Im s 2Im Br

5、c22 sin( 2/ 2)rc 22 sin( 2)2Re s 2Re Brc22cos( 2/ 2)rc22cos( 2)2Re s 3Re Drc33cos( 3/ 2)rc33cos( 3)2Im s3Im Drc33 sin( 3/ 2)rc 33 sin( 3)将上述 BC 构件， CD 构件各式合并成矩阵形式有，101000010100rc2 sin 2rc2 cos 2(r2rc2 )sin 2(r2rc2 ) cos 20000101000010100(r3rc3 ) sin 3(r3rc3 )cos 3rc3 sin 3rc3 cos 32m2ReBm2 rc 22cos(

6、 2/ 2)m2rc 22cos( 2)FX 22m2Im Bm2 rc22sin( 2/ 2)m2 rc 22 sin( 2)F y2m2 g (1-22)J22M 22m3 Re D m3rc33 cos( 3/ 2) m3 rc33 cos( 3)FX 32m3Im Dm3 rc 33sin( 3/ 2)m3 rc33sin( 3)F y3m3 gJ33M 3如图 1-5 对构件 5 进行约束反力的推导如下，RXBRyBRXC=RyCRXDRyD图 1-5 CE 杆件受力模型RxCFx5RxEm5 Re sRyCFy5RyEm5 g m5 Im sM 5RxC rc5 sin 5RyC

7、rc5 cos 5 RxE (r5 rc5 ) sin 5 RyE (r5 rc5 ) cos 5 J 5 5如图 1-6 对滑块进行受力分析如下，滑块受力模型Fx6RxERF sin6m6 Re EFy 6RyERF cos 6 m6 gm6 Im E由运动学可推，2Re sReCrc 55 cos( 5/ 2)rc55 cos( 5)2Im sRe Crc55 sin( 5/ 2)rc55 sin( 5)Re Es6 cos6Im Es6 s i n 61010001010rc 5 sin5rc 5 cos 5(r c5r5 ) sin 5( rc 5r5 ) cos 500010sin0

8、001cos2RxCRyCRxE6RyE6R F（ 1-23）m5 Re C m5 r c 55 cos(5/ 2)m5 r c55 cos(5)Fx 52m5 Re C m5 rc5 5sin( 5/ 2)m5 rc 55sin( 5)F y 5m5 gJ 55M 5m 6s 6cos6F x 6m6 s 6 sin6F y 6m 6 g二 编程与仿真利用 MATLAB进行仿真分析， 主要包括两个步骤：首先是编制计算所需要的函数模块，然后利用其仿真工具箱Simulink建立仿真系统框图，设定初始参数进行仿真分析。针对建立完成的数学模型，为了进行矩阵运算，根据以上式子编制M 函数文件cheng

9、crank.m ，chengrrr.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)

10、+pi/2)+r1*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

11、)=ddtheta-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);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(

12、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 paramet

13、ers%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)+m1*rc1*x(2)2*si

14、n(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)=dtheta-2%x(4)=dtheta-

15、3%x(5)=ddtheta-2%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;% 两杆的转动惯量ReddD=0;ImddD=0;Fx2=0; Fy2=

16、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,6)=-rc3*cos(x(2);b=zeros(6,1)

17、;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)+pi/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*si

18、n(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-5%y(2)=ddsr5=1.2; th6=0; ReddD=0; ImddD=0;a=r

19、5*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 parameters%x(1)=theta-5%x(2)=dtheta-5%x(3)=ddtheta-5%x

20、(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;Fx6=1000; Fy6=0;M5=0;th6=0;a=zeros(5);a(1,1)=1;a(1,3)=1;a(2,2)=1;a(2,4)=1;a(3,1)=rc5*sin(x(1);a(3,2)=-rc5*cos

21、(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)+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 命令根据需要绘制曲线。图 3-1四 仿真的实现再设计完成仿真框图之后，为了进行仿真还必须设定初始参数值。连杆机构杆长已经在simulink框图中给定，如果设定1 初始夹角为62 ，1 =10rad/