MATLAB机构运动仿真

机械原理第一次作业（matlab7.0）：

求：r1旋转360°时，θ2，θ3，ω2，ω3，α2，α3和C点的加速度. 设r1=1000,r2=2000,r3=1800,r4=2500

1、角位移的M函数：

function y=diyiwen(x)

% Input parameters

% x(1)=theta-1

% x(2)=theta-2 guess value

% x(3)=theta-3 guess value

% x(4)=r1

% x(5)=r2

% x(6)=r3

% x(7)=r4

% Output parameters

% y(1)=theta-2

% y(2)=theta-3

theta2=x(2);

theta3=x(3);

%

epsilon=1.0E-6;

%

f=[x(4)*cos(x(1))+x(5)*cos(theta2)-x(7)-x(6)*cos(theta3);

x(4)*sin(x(1))+x(5)*sin(theta2)-x(6)*sin(theta3)];

%

while norm(f)>epsilon

J=[-x(5)*sin(theta2) x(6)*sin(theta3);x(5)*cos(theta2) -x(6)*cos(theta3)];

dth=inv(J)*(-1.0*f);

theta2=theta2+dth(1);

theta3=theta3+dth(2);

f=[x(4)*cos(x(1))+x(5)*cos(theta2)-x(7)-x(6)*cos(theta3);

x(4)*sin(x(1))+x(5)*sin(theta2)-x(6)*sin(theta3)];

norm(f);

end;

y(1)=theta2;

y(2)=theta3;

r1旋转360°时，θ2，θ3的M文件程序：

r(1)=1000;

r(2)=2000;

r(3)=1800;

r(4)=2500;

dr=pi/180;

th(1)=0;

th(2)=59.8897 *dr;

th(3)=106.0172*dr;

y=diyiwen([th(1),th(2),th(3),r(1),r(2),r(3),r(4)])

dth=1*dr;

for i=1:360

y=diyiwen([th(1),th(2),th(3),r(1),r(2),r(3),r(4)]);

th23(i,:)=[th(1)/dr,th(2)/dr,th(3)/dr];

th(1)=th(1)+dth;

th(2)=y(1);

th(3)=y(2);

end

figure(1)

plot(th23(:,1),th23(:,2),th23(:,1),th23(:,3),th23(:,1),th23(:,1)) axis([0,360,0,360])

grid on

title('角位移线图')

xlabel('曲柄转角th(1)')

ylabel('从动件转角th(2),th(3)')

text(300,50,'角th(2)')

text(300,150,'角th(3)')

text(200,200,'角th(1)')

2、角速度的M函数：

function y=dierwen(x)

% Iput parameters

%

% x(1)=theta-1

% x(2)=theta-2

% x(3)=theta-3

% x(4)=dtheta-1

% x(5)=r1

% x(6)=r2

% x(7)=r3

%

% Output parameters

%

% y(1)=dtheta-2

% y(2)=dtheta-3

%

A=[-x(6)*sin(x(2)) x(7)*sin(x(3));

x(6)*cos(x(2)) -x(7)*cos(x(3))];

B=[x(5)*sin(x(1));-x(5)*cos(x(1))]*x(4);

y=inv(A)*B;

ω2，ω3

r1旋转360°时ω2，ω3的M文件程序：

th(1)=0;

dr=pi/180;

dth(1)=10;

dth(2)=59.8897*dr;

dth(3)=106.0172*dr;

dth1=1*dr;

r(1)=1000;

r(2)=2000;

r(3)=1800;

r(4)=2500;

for j=1:360

x=dierwen([th(1),dth(2),dth(3),dth(1),r(1),r(2),r(3)]);

y=diyiwen([th(1),dth(2),dth(3),r(1),r(2),r(3),r(4)]);

dth(2)=y(1);

dth(3)=y(2);

th(1)=th(1)+dth1;

dth23(j,:)=[th(1)/dr,x(1)/dr,x(2)/dr];

end

figure(2)

plot(dth23(:,1),dth23(:,2),dth23(:,1),dth23(:,3))

axis([0,360,-400,400])

grid on

title('角速度线图')

xlabel('曲柄旋转角度')

ylabel('从动件角速度dth(2),dth(3)(rad/s)')

text(250,200,'连杆2角速度')

text(50,100,'摇杆3角速度')

3、r1旋转360°时，α2，α3程序M文件：

r(1)=1000;

r(2)=2000;

r(3)=1800;

r(4)=2500;

th(1)=0;

dr=pi/180;

dth(2)=59.8897 *dr;

dth(3)=106.0172*dr;

dth(1)=10;

dth2=1*dr;

for k=1:360

x=dierwen([th(1),dth(2),dth(3),dth(1),r(1),r(2),r(3)]);

y=diyiwen([th(1),dth(2),dth(3),r(1),r(2),r(3),r(4)]);

dth(2)=y(1);

dth(3)=y(2);

w1=x(1);

w2=x(2);

C=[-r(2)*sin(dth(2)),r(3)*sin(dth(3));r(2)*cos(dth(2)),-r(3)*cos(dth(3))];

D=-[-w1*r(2)*cos(dth(2)),w2*r(3)*cos(dth(3));-w2*r(2)*sin(dth(2)),w2*r(3)*sin(dth(3))]…*[w1;w2]+10*[10*r(1)*cos(th(1));10*r(1)*sin(th(1))];

a=inv(C)*D;

a2(k)=a(1);

a3(k)=a(2);

jiao(k)=th(1);

th(1)=th(1)+dth2;

end

figure(3)

plot(jiao/dr,a2,jiao/dr,a3)

axis([0,360,-100,100])

grid on

title('角加速度线图')

xlabel('曲柄旋转转角')

ylabel('从动件角加速度（rad/s^2）')

text(200,40,'连杆2的角加速度')

text(150,-40,'摇杆3的角加速度')

4、C点的加速度函数：

C点的M函数文件：

function y=cdjsd(x)

% x(1)=r(1);

% x(2)=r(2);

% x(3)=th(1);

% x(4)=th(2);

% x(5)=w1;

% x(6)=w2;

% x(7)=a(1);

y=[-x(1)*sin(x(3)),-x(2)*sin(x(4));x(1)*cos(x(3)),x(2)*cos(x(4))]*[0;x(7)]+[-x(1)*x(5)*cos(x(3)), -x(2)*x(6)*cos(x(4));-x(1)*x(5)*sin(x(3)),-x(2)*x(6)*sin(x(3))]*[x(5)^2;x(6)^2];

r1旋转360°时，C点的加速度的M文件程序：

r(1)=1000;

r(2)=2000;

r(3)=1800;

r(4)=2500;

th(1)=0;

dr=pi/180;

dth(2)=59.8897*dr;

dth(3)=106.0172*dr;

dth(1)=10;

dth2=1*dr;

for k=1:360

x=dierwen([th(1),dth(2),dth(3),dth(1),r(1),r(2),r(3)]);

y=diyiwen([th(1),dth(2),dth(3),r(1),r(2),r(3),r(4)]);

dth(2)=y(1);

dth(3)=y(2);

w1=x(1);

w2=x(2);

C=[-r(2)*sin(dth(2)),r(3)*sin(dth(3));r(2)*cos(dth(2)),-r(3)*cos(dth(3))];

D=-[-w1*r(2)*cos(dth(2)),w2*r(3)*cos(dth(3));-w2*r(2)*sin(dth(2)),w2*r(3)*sin(dth(3))]*[w1;w 2]+10*[10*r(1)*cos(th(1));10*r(1)*sin(th(1))];

a=inv(C)*D;

a2=a(1);

ac=cdjsd([r(1),r(2),th(1),dth(2),w1,w2,a2]);

th(1)=th(1)+dth2;

acx(k)=ac(1);

acy(k)=ac(2);

jiao(k)=th(1);

end

plot(jiao/dr,acx,jiao/dr,acy)

axis([0,360,-300000,300000])

grid on

title('C点加速度线图')

xlabel('曲柄旋转角度')

ylabel('C点x,y方向的加速度') text(300,100000,'x方向')

text(300,-100000,'y方向')