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

一 系统模型成立

为了对机构进行仿真分析，第一必需成立机构数学模型，即位置方程，然后利用MATLAB仿真分析工具箱Simulink对其进行仿真分析。图3．24所示是由原动件(曲柄1)和RRR—RRP六杆机构。各构件的尺寸为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·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-2关于曲柄AB由理论力学能够列出表达式：

RXARXBFX1m1Re•s•1

R••yARyBFy1m1Ims1

M1MFRXArc1sin1RyArc1cos1RXB(r1rc1)sin1RyB(r1rc1)cos1J••11由运动学知识能够推得：

Re•s••••••21ReArc11cos(1/2)rc11cos(1) Im•s••••••21ImArc11sin(1/2)rc11sin(1)

将上述各式归并成矩阵形式有，

RXA•••••2Rm1ReAm1rc11cos(1/2)m1rc11cos(1)FX1R1XByAm••m•••2(1-21) M1ImA1rc11sin(1/2)m1rc11sin(1)Fy1RyBm1g

如图1-3，对构件BC的约束反力推导如下，

1-2 AB

图

图1-3 BC构件受力模型

••RXBFX2RXCm2Res2

RyBFy2RyCm2gm2Ims2

••M2RXBrc2sin2RyBrc2cos2RXC(r2rc2)sin2RyC(r2rc2)cos2J22

如图1-4，对构件BC的约束反力推导如下，

••

图 1-4 CD构件受力模型

••RXDFX3RXCm3Res3

RyDFy3RyCm3gm3Ims3

M3RXDrc3sin3RyDrc3cos3RXC(r3rc3)sin3RyC(r3rc3)cos3J33

由运动学能够推导得，

••••Ims2ImBrc22sin(2/2)rc22sin(2) Res2ReBrc22cos(2/2)rc22cos(2)•••••••2•••••••2Re•s••••••23ReDrc33cos(3/2)rc33cos(3) Im•s••••••23ImDrc33sin(3/2)rc33sin(3)

将上述BC构件，CD构件各式归并成矩阵形式有，

101000010100rrrcosc2sin2c2cos2(2rc2)sin2(r2rc2)20000101000010100(r3rc3)sin3(r3rc3)cos3rc3sin3rc3cos3•••••2mReBmr/2)m22c22cos(22rc22cos(2)FX2•••••2m2ImBm2rc22sin(2/2)m2rc22sin(2)Fy2m2g (1-22)

J2••2M2m3ReD••2m•••3rc33cos(3/2)m3rc33cos(3)FX3mImD••mr••sin(/2)mr•2sin()F33c3333c333y3m3gJ3••3M3如图1-5 对构件5进行约束反力的推导如下，

图1-5 CE杆件受力模型

RR••xCFx5xEm5Res

RyCFy5RyEm••5gm5Ims

M5RxCrc5sin5RyCrc5cos5RxE(r5rc5)sin5RyE(r5rc5)cos5J5••5如图1-6 对滑块进行受力分析如下，

RXBRyBRXCR=

yCRXDRyD

滑块受力模型

Fx6RxERFsin6m6ReE

••Fy6RyERFcos6m6gm6ImE

由运动学可推，

•••••••2••ResReCrc55cos(5/2)rc55cos(5) ImsReCrc55sin(5/2)rc55sin(5)

ReEs6cos6 ImEs6sin6

10rc5sin500•••••••••2••••••••01rc5cos500••10(rc5r5)sin510•201(rc5r5)cos501RxCRyC0RxEsin6RyEcos6RF00mReCmrcos(/2)mrcos()F5555c555c55x5•••••2mReCmrsin(/2)mrsin()Fmg5555c555c55y55••JM555••m6s6cos6Fx6••m6s6sin6Fy6m6g（1-23）

二 编程与仿真

利用MATLAB进行仿真分析，要紧包括两个步骤：第一是编制计算所需要的函数模块，然后利用其仿真工具箱Simulink成立仿真系统框图，设定初始参数进行仿真分析。针对成立完成的数学模型，为了进行矩阵运算，依照以上式子编制M函数文件chengcrank.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)=Re[ddB] %y(2)=Im[ddB] r1=0.4;

ddB=[r1*x(3)*cos(x(1)+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)=Re[ddB] %x(6)=Im[ddB]

%Output parameters %y(1)=ddtheta-2 %y(2)=ddtheta-3 %y(3)=Re[ddC] %y(4)=Im[ddC]

r2=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(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 parameters %y(1)=RxA %y(2)=RyA %y(3)=M1

g=9.8; %重力加速度 r1=0.4; %曲柄长度

rc1=0.2;%质心离铰链A的距离 m1=1.2;%曲柄质量

J1=0.016; %绕质心转动惯量

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-rc1)*cos(x(1))-MF;

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-3 %x(5)=ddtheta-2 %x(6)=ddtheta-3 %x(7)=Re[ddB] %x(8)=Im[ddB]

%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)=RyD

g=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*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);

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*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)=Re[ddC] %x(4)=Im[ddC] %x(5)=ds

%Output parameters %y(1)=ddtheta-5 %y(2)=dds

r5=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 parameters %x(1)=theta-5 %x(2)=dtheta-5 %x(3)=ddtheta-5 %x(4)=dds-6 %x(5)=Re[ddC] %x(6)=Im[ddC] %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(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;

三 系统仿真框图

进入MATLAB，在命令栏中键入Simulink进入仿真界面，依照信息传递的逻辑关系，成立仿真系统框图如图3-1.

然后设定各环节的初始参数，即能够对机构进行运动学仿真分析，再利用MATLAB的plot命令依照需要绘制曲线。

图3-1

四 仿真的实现

再设计完成仿真框图以后，为了进行仿真还必需设定初始参数值。连杆机构杆长已经在simulink框图中给定，若是设定1初始夹角为62，1=10 rad/s，曲柄1作匀速转动（即10），接下来要确信杆2，3的角位移和角速度，杆5的角位移和角速度，滑块的速度。

能够利用辛普森方式（在MATLAB命令框中输入M函数为rrrposi）求得2=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.3612 1.8101 -2.2345 3.3250 -41*pi/180 0.6962 -3.1323];其中初始数值别离对应：theta-一、omega-一、theta-二、omega-二、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

移动副F的约束反力

曲柄力矩所做的功W1