matlab曲柄连杆机构分析讲课讲稿

m a t l a b曲柄连杆机构

分析

clear;clc;

n=750;l=0.975;R=0.0381;h=0.2;omiga=n.*pi/30;tmax=2.*pi/omiga;

t=0:0.001:tmax; %计算曲柄转一圈的总t值

alpha1=atan((h+R.*sin(omiga.*t))./sqrt(l.*l-(h+R.*sin(omiga.*t))))+pi;

alpha1p=-(R.*omiga.*cos(omiga.*t))./(l.*cos(alpha1));

vb=-R.*omiga.*sin(omiga.*t)+R.*omiga.*cos(omiga.*t).*tan(alpha1);

ab=-R.*omiga.^2.*cos(omiga.*t)-

(R.*omiga.*cos(omiga.*t)).^2./(l.*(cos(alpha1)).^3)

-R.*omiga.^2.*sin(omiga.*t).*tan(alpha1);

subplot(1,2,1);plot(t,vb);title('曲柄滑块机构的滑块v-t图');

xlabel('时间t(曲柄旋转一周)');ylabel('滑块速度v');grid on;

subplot(1,2,2);plot(t,ab);title('曲柄滑块机构的滑块a-t图');

xlabel('时间t(曲柄旋转一周)');ylabel('滑块加速度a');grid on;

%下面黄金分割法求滑块的速度与加速度最大值

epsilon=input('根据曲线初始区间已确定，请输入计算精度epsilon(如输入

0.001):');

a=0;b=0.04; %初始区间

n1=0; %n1用于计算次数

a1=b-0.618*(b-a);y1=-R.*omiga.*sin(omiga.*a1)

+R.*omiga.*cos(omiga.*a1).*tan(alpha1);

a2=a+0.618*(b-a);y2=-R.*omiga.*sin(omiga.*a2)

+R.*omiga.*cos(omiga.*a2).*tan(alpha1);

while abs(a-b)>=epsilon

if y1<=y2

b=a2;a2=a1;y2=y1;

a1=b-0.618*(b-a);y1=-

R.*omiga.*sin(omiga.*a1)+R.*omiga.*cos(omiga.*a1).*tan(alpha1);

else

a=a1;a1=a2;y1=y2;

a2=a+0.618*(b-a);y2=-

R.*omiga.*sin(omiga.*a2)+R.*omiga.*cos(omiga.*a2).*tan(alpha1);

end

n1=n1+1;

end

vbm1=omiga*(a+b)/2;

disp(['经过',num2str(n1),'次计算,用黄金分割法找到速度最大值对应的wt是:', num2str(vbm1),'弧度。'])

a=0.04;b=0.08;%初始区间变化，对应的函数取负即可

n1=0; %n1用于计算迭代次数

a1=b-0.618*(b-a);y1=-(-

R.*omiga.*sin(omiga.*a1)+R.*omiga.*cos(omiga.*a1).*tan(alpha1));

a2=a+0.618*(b-a);y2=-(-

R.*omiga.*sin(omiga.*a2)+R.*omiga.*cos(omiga.*a2).*tan(alpha1));

while abs(a-b)>=epsilon

if y1<=y2

b=a2;a2=a1;y2=y1;

a1=b-0.618*(b-a);

y1=-(-R.*omiga.*sin(omiga.*a1)+R.*omiga.*cos(omiga.*a1).*tan(alpha1));

else

a=a1;a1=a2;y1=y2;

a2=a+0.618*(b-a);

y2=-(-R.*omiga.*sin(omiga.*a2)+R.*omiga.*cos(omiga.*a2).*tan(alpha1));

end

n1=n1+1;

end

vbm2=omiga*(a+b)/2;

disp(['经过',num2str(n1),'次计算,找到的另一个速度最大值对应的w

是:',num2str(vbm2),'弧度。'])

disp(['两个速度最大值点相差是:',num2str(vbm2-vbm1),'弧度(即',

num2str(180*(vbm2- vbm1)/pi),'度)。'])

a=0.02;b=0.06;%对应的函数取负即可

n1=0; a1=b-0.618*(b-a);

y1=-(-R.*omiga.^2.*cos(omiga.*a1)-(R.*omiga.*cos(omiga.*a1)).^2./

(l.*(cos(alpha1)).^3)-R.*omiga.^2.*sin(omiga.*a1).*tan(alpha1));

a2=a+0.618*(b-a);

y2=-(-R.*omiga.^2.*cos(omiga.*a2)-(R.*omiga.*cos(omiga.*a2)).^2./

(l.*(cos(alpha1)).^3)-R.*omiga.^2.*sin(omiga.*a2).*tan(alpha1));

while abs(a-b)>=epsilon

if y1<=y2

b=a2;a2=a1;y2=y1;a1=b-0.618*(b-a);

y1=-(-R.*omiga.^2.*cos(omiga.*a1)-(R.*omiga.*cos(omiga.*a1)).^2./ (l.*(cos(alpha1)).^3)-R.*omiga.^2.*sin(omiga.*a1).*tan(alpha1));

else

a=a1;a1=a2;y1=y2;a2=a+0.618*(b-a);

y2=-(-R.*omiga.^2.*cos(omiga.*a2)-(R.*omiga.*cos(omiga.*a2)).^2./ (l.*(cos(alpha1)).^3)-R.*omiga.^2.*sin(omiga.*a2).*tan(alpha1));

end

n1=n1+1;

end

disp(['经过',num2str(n1),'次计算,找到的加速度最大值对应的w

是:',num2str(omiga*(a+b)/2)])

disp(['速度最大值点与加速度最大值相差的wt

是:',num2str((omiga*(a+b)/2)-vbm1),'弧度。'])