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

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

2、.016kg 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级杆组的动力学，再分别对这三个模型进行相应参数的求解。21图1-2 AB构件受力模型如上图1-2对于曲柄AB由理论力学可以列出表达式：R xa - rxbFxi = mi Re si*RyA -RyBFyi = mi Im s

3、iRxa IRyAM 一Mi M f RxaQSin 弓 - RyAG costRxb（A-r/si-RyB（ - r cost=Ji二i由运动学知识可以推得:Resi = ReA rc样 i cos 二 /2） rci cos二）* * Im si = Im A rci sin（弓二 /2） rc2;i sin（冇 二）将上述各式合并成矩阵形式有,（i-2i）migmi ReA gQ 力 cos（弓二 /2） mirci cos（刊 二）FXi RXB龄*#2mi Im A mirci sin（片二 /2）mirci sin（齐二）FyiRyB如图i-3,对构件BC的约束反力推导如下,i1图

4、1-3 BC构件受力模型RXBFX2RXC = m2 Re S 2RyB Fy2 RyeIm S2M 2 RxbQ Sin 二2 - RyBc2 cos Rxc （2 - Q）sin 出 - Rye （2 - rc2）cos2 = J2 二 2 如图1-4,对构件BC的约束反力推导如下，RXD FX3 1 RXC = m3 Re S 3RyDFy3 Rye msg =m3lm S3M3Rxdrc3 sin-RyDQ cos - Rxc（D- s）sin 岂 -Rye（D - g）cos= J3 3由运动学可以推导得,i1Im S22W*=Im B rc2 二2 sin（二2 二 /2）rc2

5、sin（r2二）2 Res2 二 ReB rc2 讥 cos（r 二 /2） rc2 cosd 二）Res3 二 ReD a 3 cos（k 二 / 2） g cos3 二）科材2Im S3 = Im D rc3 3 sin（二3 什第/2） rc3 sin（）3 出第）将上述BC构件，CD构件各式合并成矩阵形式有，-101000 I _RXB |010100RyBrc2 si n$_C2 CO羽2一（2 rc2）si2（2 -rc2）C0200Rxc00-1010RyC000-101Rxd:.00（3 -rc3）si3-（3 - rc3）CO3R3S in&-rc3 cod 一iRyD- 1

6、*m2 ReB + m2rc2 日2 cos2 +兀 /2） +m2rc20 2 cos2 十兀）_ FX200*2（1-22）m2 Im B + rni2c2。2 sin（日2 + 兀 /2） + rni2c2 日 2 sin（2 + 兀）一 Fy2 + mgJ2 日 2 - M 2 .2m3 ReD + mhQ 日3 cosd +71 /2） +口3心日3 cos3 +巧Fx3 * *2m3 Im D + m3rc3&3sin（爲 + 兀 /2） +口3心日 3 sin但3 + 巧 一 Fy3 +J3 8 3 M31如图1-5对构件5进行约束反力的推导如下,图1-5 CE杆件受力模型12*

7、Rxc F x5 + Rxe = m5 Re sRycFy5RyE _m5g =m5lm sM 5 - RxC rc5 si n T5 - RyC rc5 cosi5 RxE （5 - rc5 ） Si 门応，RyE （5 _ rc5 ） COS = J 5 T 5 如图1-6对滑块进行受力分析如下，-RyEFyERfF x6F y6/FxE -RxE滑块受力模型e-RxE -Rf sin v6 = m6 Re E-RyE Rf cos6 -m6g =m6由运动学可推，ImE*2Res ReC rc5 5cos（二5 二 /2） rc5 5 cos5二）冷冷.2Im s = ReC rc5 珀

8、 sin（屯 : /2） rc5 usin（v5二）* *Re E = S6 cos 玉 Im E = S6 s i n6-10 sin00-01-C5 Cos K0010（C5 -5）sin 丸-10一 （Tc501-5）cos 屯0-1000 sin=6COS 二6EcRyCRxERyE JLR（1-23）25c5 71 52COS（二5 工）-Fx5Re C m5c5 715 cos（ % 二 / 2） m5* * m5 Re C +m5c5 日5 sin（日5 + 兀 / 2） + m5c5 日5 sin（日5 + 兀）一 Fy5 + msgJ 5 日5 M 5m6 S6 cos 日6

9、 - Fx6vm6 S6 sin 日6 Fy6 +m6gm51二编程与仿真利用MATLAB进行仿真分析，主要包括两个步骤：首先是编制计算所需要的函数模块，然后利用其仿真工具箱Simulink建立仿真系统框图，设定初始参数进行仿真分析。针对建立完成的数学模型，为了进行矩阵运算，根据以上式子编制M函数文件chengcrank.m ，chengrrr.m、chengcrankdy.m、chengrrrdy.m、chengrrp.m 禾口 chengrrpdy.m 女口下：曲柄原动件M函数文件chengcrank.m：function y=che ngcra nk(x)%F unction to co

10、mpute the accleratio n of crank%ln put 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)A2*cos(x(1)+pi);r1*x(3)*s in (x(1)+pi/2)+r1*x (2) A2*si n(x(1)+pi)；y=ddB;RRR II级杆组M函数文件chengrrr.m :function y=che ngrrr(x)%fun c

11、ti on to compute the accelerati on for RRR bar group%ln put parameters%x(1)=theta-2%x(2)=theta-3%x(3)=dtheta-2%x(4)=dtheta-3%x(5)=ReddB%x(6)=lmddB%Output parameters%y(1)=ddtheta-2%y(2)=ddtheta-3%y(3)=ReddC%y(4)=lmddCr2=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)

12、+pi/2) -r3*sin(x(2)+pi/2);11b=-r2*cos(x(1)+pi)r3*cos(x(2)+pi);-r2*si n(x(1)+pi)r3*si n(x(2)+pi)*x (3)A2;x (4) A2+ReddD-x(5);ImddD-x (6);ddth=i nv(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)A 2*cos(x(1)+pi); y(4)=x(6)+r2*ddth(1)*si n(x(1)+pi/2)+r2*x(3)A2*si n(x(1)+pi);曲

13、柄原动件动力学M函数文件 chengcrankdy.m：function y=che ngcra nkdy(x)%Function for Dyanmic an alysis of crank%l nput parameters%x(1)=theta-1%x(2)=dtheta-1%x(3)=ddtheta-1%x(4)=RxB%x(5)=RyB%0utput parameters%y(1)=RxA%y(2)=RyA%y(3)=M1g=9.8; %重力加速度r1=0.4; %曲柄长度rc仁0.2;%质心离铰链 A的距离m1=1.2;%曲柄质量J1=0.016; %绕质心转动惯量Fx仁0; Fy

14、仁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) A2*cos(x(1)+pi)-Fx1+x(4); y(2)=m1*lmddA+m1*rc1*x (3)*si n(x(1)+pi/2)+m1*rc1*x (2) A2*si n(x(1)+pi)-Fy1+x(5)+m1*g;y(3)=J1*x (3)-y(1)*rc1*si n(x(1)+y(2)*rc1*cos(x(1)-x(4)*(r1-rc1)*si n(x(1)+x (5)*(r1

15、-rc1)*cos(x (1)-Mf;RRR II级杆组动力学 M函数文件 chengrrrdy.m :function y=che ngrrrdy(x)%Function for Dyanmic an alysis of RRR dayard group%ln put parameters%x(1)=theta-2%x(2)=theta-3%x(3)=dtheta-2%x(4)=dtheta-3%x(5)=ddtheta-221%x(6)=ddtheta-3%x(7)=ReddB%x(8)=lmddB%x(9)=Fx3%x(10)=Fy3%x(11)=M3%Output parameters

16、%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=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*si n(x(1);a(3,2)=-rc2*cos(x(

17、1);a(3,3)=-(r2-rc2)*si n(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)*si n(x (2);a(6,4)=-(r3-rc3)*cos(x(2);a(6,5)=rc3*si n(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)A2*cos(x(1)+pi)-Fx2;b(2,1)=m2*rc2*x(5)*s in (

18、x(1)+pi/2)+m2*x(8)+m2*rc2*x (3)A2*si n(x(1)+pi)-Fy2+m2*g; b(3,1)=J2*x(5)-M2;21b(4,1)=m3*rc3*x(6)*cos(x(2)+pi/2)+m3*ReddD+m3*rc3*x(4)A2*cos(x(2)+pi)-x(9); b(5,1)=m3*rc3*x(6)*s in(x(2)+pi/2)+m3*ImddD+m3*rc3*x(4)A2*si n(x(2)+pi)-x(10)+m3*g; b(6,1)=J3*x(6)-x(11);y=i nv(a)*b;RRP II级杆组M函数文件：fun ctio n y=c

19、he ngrrp(x)%fun cti on to compute the accelerati on for RRP bar group%ln put parameters%x(1)=theta-5%x(2)=dtheta-5%x(3)=ReddC%x(4)=lmddC%x(5)=ds%Output parameters%y(1)=ddtheta-5%y(2)=ddsr5=1.2; th6=0; ReddD=0; ImddD=0;a=r5*cos(x(1)+pi/2)-cos(th6); r5*si n(x(1)+pi/2)-si n( th6);b=-r5*cos(x(1)+pi) 0;

20、-r5*sin(x(1)+pi) 0*卜(2)人2; x(5)+ReddD-x(3); ImddD-x(4);y=i nv(a)*b;RRP II级杆组动力学M函数文件：function y=che ngrrpdy(x)%Function for Dyanm5c an alysis of RRP dayard group%ln put parameters%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%

21、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;11Fx6=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*si n(x(1);a(3,2)=-rc5*cos(x(1);a(3,3)=-(r5-rc5)*si n(x(1);a(3,4)=(r5-rc5)*cos(x(1);a(4,3)=-1;a(4,5)=-s

22、in (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)A2*cos(x(1)+pi)-Fx5;b(2,1)=m5*x(6)+m5*rc5*x(3)*s in (x(1)+pi/2)+m5*rc5*x (2)A2*si n(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)*si n(th6)-Fx6+m6*g;y=i nv(a)*b;三系统仿

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