求助!!!急 曲柄四杆机构的程序!!!
谁能帮我修改一下程序啊,非常感谢!!!!程序如下:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输入已知数据
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
l1=40;
l2=50;
l3=75;
l4=35;
l6=70;
l5=60;
omega1=10;
omega6=0;
hd=pi/180;
du=180/pi;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 四杆机构曲线设计
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n1=1:asin(6/7)*du
theta4(n1)=n1+asin((l6/l5)*sin(n1*hd))*du; % 杆DF转过的角度
l7(n1)=l5*sin(theta4(n1)*hd)/sin(n1*hd); % 杆AF的长度
l8(n1)=sqrt((l6-l4*cos(theta4(n1)*hd)+l1*cos((n1)*hd))^2+(-l4*sin(theta4(n1)*hd)+l1*sin((n1)*hd))^2); % BD的长度
theta2(n1)=(n1*hd-acos(((l7(n1)+l1)^2+l8(n1)^2-(l4+l5)^2)/(2*(l7(n1)+l1)*l8(n1)))+acos((l2^2+l8(n1)^2-l3^2)/(2*l2*l8(n1))))*du; % 杆BC转过的角度
theta3(n1)=(theta4(n1)*hd+acos((l8(n1)^2+(l4+l5)^2-(l1+l7(n1))^2)/(2*l8(n1)*(l4+l5)))-acos((l8(n1)^2+l3^2-l2^2)/(2*l3*l8(n1))))*du; % 杆CD转过的角度
A=[-l2*sin(theta2(n1)*hd),l3*sin(theta3(n1)*hd),l4*sin(theta4(n1)*hd+pi),0;
l2*cos(theta2(n1)*hd),-l3*cos(theta3(n1)*hd),-l4*cos(theta4(n1)*hd+pi),0;
0,0,l5*sin(theta4(n1)*hd),cos(n1*hd);
0,0,-l5*cos(theta4(n1)*hd),-sin(n1*hd)];
B=[l1*sin(n1*hd+pi);-l1*cos(n1*hd+pi);l7(n1)*sin(n1*hd); -l7(n1)*cos(n1*hd)];
P1=(A\B )*omega1 ;
omega2(n1)=P1(1);
omega3(n1)=P1(2);
omega4(n1)=P1(3);
v7=P1(4);
vcx(n1)=-l1*omega1*sin(n1*hd+pi)-l2*omega2(n1)*sin(theta2(n1)*hd);
vcy(n1)=l1*omega1*cos(n1*hd+pi)+l2*omega2(n1)*cos(theta2(n1)*hd);
vc(n1)=sqrt(vcx(n1)^2+vcy(n1)^2); % C点的速度
C=[-l2*sin(theta2(n1)*hd),l3*sin(theta3(n1)*hd),-l4*sin(theta4(n1)*hd),0;
-l2*cos(theta2(n1)*hd),-l3*cos(theta3(n1)*hd),l4*cos(theta4(n1)*hd),0;
0,0,l5*sin(theta4(n1)*hd),cos(n1*hd);
0,0,-l5*cos(theta4(n1)*hd),sin(n1*hd)];
D=[-l2*omega2(n1)*cos(theta2(n1)*hd),l3*omega3(n1)*cos(theta3(n1)*hd),-l4*omega4(n1)*cos(theta4(n1)*hd),0;
-l2*omega2(n1)*sin(theta2(n1)*hd),l3*omega3(n1)*sin(theta3(n1)*hd),-l4*omega4(n1)*sin(theta4(n1)*hd),0;
0,0,l5*omega4(n1)*cos(theta4(n1)*hd),-omega1*sin(n1*hd);
0,0,l5*omega4(n1)*sin(theta4(n1)*hd),omega1*cos(n1*hd)];
E=[l1*omega1*cos(n1*hd);
l1*omega1*sin(n1*hd);
l7*omega1*cos(n1*hd)+v7*sin(n1*hd);
l7*omega1*sin(n1*hd)+v7*cos(n1*hd);];
P2=C\(E*omega1-D*P1);
alpha2(n1)=P2(1);
alpha3(n1)=P2(2);
alpha4(n1)=P2(3);
acx(n1)=l1*omega1^2*cos(n1*hd)-l2*omega2(n1)^2*cos(theta2(n1)*hd)-l2*alpha2(n1)*sin(theta2(n1)*hd);
acy(n1)=l1*omega1^2*sin(n1*hd)-l2*omega2(n1)^2*sin(theta2(n1)*hd)-l2*alpha2(n1)*cos(theta2(n1)*hd);
ac(n1)=sqrt(acx(n1)^2+acy(n1)^2); % C点的加速度
end
disp 'vc(50)=',vc(50)
disp 'ac(50)=',ac(50)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 位移,角速度,角加速度和四杆机构图形输出
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
n1=1:asin(6/7)*du;
subplot(2,2,1); % 绘制C点的速度曲线
plot(n1,vc)
title('c点速度线图');
xlabel('l7转角 \theta_1 / \circ')
ylabel('C点速度 / mm\cdots^{-1}')
grid on;hold on;
text(100,-300,'v_c')
subplot(2,2,2); % 绘制C点的加速度曲线
plot(n1,ac)
title('c点加速度线图');
xlabel('l7转角 \theta_1 / \circ')
ylabel('C点加速度 / mm\cdots^{-2}')
grid on;hold on;
text(100,-300,'a_c')
subplot(2,2,3); % 设计铰链四杆机构
x(1)=0;
y(1)=0;
x(2)=l7(50)*cos(50*hd);
y(2)=l7(50)*sin(50*hd);
x(3)=(l5+l6+10)*cos(50*hd);
y(3)=(l5+l6+10)*sin(50*hd);
x(4)=l7(50)*cos(50*hd);
y(4)=l7(50)*sin(50*hd);
x(5)=l7(50)*cos(50*hd)+4*cos(pi/2-50*hd);
y(5)=l7(50)*sin(50*hd)-4*sin(pi/2-50*hd);
x(6)=(l7(50)+8)*cos(50*hd)+4*cos(pi/2-50*hd);
y(6)=(l7(50)+8)*sin(50*hd)-4*sin(pi/2-50*hd);
x(7)=(l7(50)+8)*cos(50*hd)-4*cos(pi/2-50*hd);
y(7)=(l7(50)+8)*sin(50*hd)+4*sin(pi/2-50*hd);
x(8)=(l7(50)-8)*cos(50*hd)-4*cos(pi/2-50*hd);
y(8)=(l7(50)-8)*sin(50*hd)+4*sin(pi/2-50*hd);
x(9)=(l7(50)-8)*cos(50*hd)+4*cos(pi/2-50*hd);
y(9)=(l7(50)-8)*sin(50*hd)-4*sin(pi/2-50*hd);
x(10)=l7(50)*cos(50*hd)+4*cos(pi/2-50*hd);
y(10)=l7(50)*sin(50*hd)-4*sin(pi/2-50*hd);
x(11)=l7(50)*cos(50*hd);
y(11)=l7(50)*sin(50*hd);
x(12)=70;
y(12)=0;
x(13)=l6-l4*cos(theta4(50)*hd);
y(13)=-l4*sin(theta4(50)*hd);
x(14)=l2*cos(theta2(50)*hd)-l1*cos(50*hd);
y(14)=l2*sin(theta2(50)*hd)-l1*sin(50*hd);
x(15)=-l1*cos(50*hd);
y(15)=-l1*sin(50*hd);
x(16)=0;
y(16)=0;
plot(x,y);
title('铰链四杆机构');
grid on;
hold on;
xlabel('mm')
ylabel('mm')
axis([-40 140 -40 140]);
plot(x(1),y(1),'o');
plot(x(2),y(2),'o');
plot(x(11),y(11),'o');
plot(x(12),y(12),'o');
plot(x(13),y(13),'o');
plot(x(14),y(14),'o');
plot(x(15),y(15),'o');
plot(x(16),y(16),'o');
text(-20,-15,'l_1 ');
text(0,-25,'l_2 ');
text(50,-30,'l_3 ');
text(80,-15,'l_4 ');
text(65,30,'l_5');
text(40,10,'l_6 ');
text(15,30,'l_7');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 四杆机构运动仿真
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
m=moviein(20);
j=0;
for n1=1:asin(6/7)*du
j=j+1;
clf;
x(1)=0;
y(1)=0;
x(2)=l7(n1)*cos(n1*hd);
y(2)=l7(n1)*sin(n1*hd);
x(3)=(l5+l6+10)*cos(n1*hd);
y(3)=(l5+l6+10)*sin(n1*hd);
x(4)=l7(n1)*cos(n1*hd);
y(4)=l7(n1)*sin(n1*hd);
x(5)=l7(n1)*cos(n1*hd)+4*cos(pi/2-n1*hd);
y(5)=l7(n1)*sin(n1*hd)-4*sin(pi/2-n1*hd);
x(6)=(l7(n1)+8)*cos(n1*hd)+4*cos(pi/2-n1*hd);
y(6)=(l7(n1)+8)*sin(n1*hd)-4*sin(pi/2-n1*hd);
x(7)=(l7(n1)+8)*cos(n1*hd)-4*cos(pi/2-n1*hd);
y(7)=(l7(n1)+8)*sin(n1*hd)+4*sin(pi/2-n1*hd);
x(8)=(l7(n1)-8)*cos(n1*hd)-4*cos(pi/2-n1*hd);
y(8)=(l7(n1)-8)*sin(n1*hd)+4*sin(pi/2-n1*hd);
x(9)=(l7(n1)-8)*cos(n1*hd)+4*cos(pi/2-n1*hd);
y(9)=(l7(n1)-8)*sin(n1*hd)-4*sin(pi/2-n1*hd);
x(10)=l7(n1)*cos(n1*hd)+4*cos(pi/2-n1*hd);
y(10)=l7(n1)*sin(n1*hd)-4*sin(pi/2-n1*hd);
x(11)=l7(n1)*cos(n1*hd);
y(11)=l7(n1)*sin(n1*hd);
x(12)=70;
y(12)=0;
x(13)=l6-l4*cos(theta4(n1)*hd);
y(13)=-l4*sin(theta4(n1)*hd);
x(14)=l2*cos(theta2(n1)*hd)-l1*cos(n1*hd);
y(14)=l2*sin(theta2(n1)*hd)-l1*sin(n1*hd);
x(15)=-l1*cos(n1*hd);
y(15)=-l1*sin(n1*hd);
x(16)=0;
y(16)=0;
plot(x,y);
title('铰链四杆机构');
grid on;
hold on;
xlabel('mm')
ylabel('mm')
axis([-100 150 -100 150]);
plot(x(1),y(1),'o');
plot(x(2),y(2),'o');
plot(x(11),y(11),'o');
plot(x(12),y(12),'o');
plot(x(13),y(13),'o');
plot(x(14),y(14),'o');
plot(x(15),y(15),'o');
plot(x(16),y(16),'o');
m(j)=getframe;
end
axis([-40 140 -40 140]);
xlabel('mm')
ylabel('mm')
movie(m)
[[it] 本帖最后由 lvmeng_2008 于 2008-6-20 19:36 编辑 [/it]]