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])