matlab曲柄连杆机构分析

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),'弧度。'])

曲柄连杆机构动力学分析图

clear;clc;close all;

n=750;h=0.2;f=0.15;g=9.8;l=0.975;m2=3.8;R=0.0381;m1=2.1;m3=2.1*2.3;

omega=n*pi/30;tmax=2*pi/omega;

t=linspace(0,tmax,100);%计算曲柄转一周的总t值

for n=1:100

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

alpha1p(n)=-(R.*omega.*cos(omega*t(n)))./(l.*cos(alpha1(n)));%alpha1的一阶导数

alpha1pp(n)=(R.*omega.^2.*sin(omega.*t(n))+l.*sin(alpha1(n)).*alpha1p(n).^2)./

(l.*cos(alpha1(n)));%alpha1的二阶导数

xA=R.*cos(omega*t(n)); yA=R.*sin(omega*t(n));

xApp=-R.*omega.^2.*cos(omega.*t(n));

yApp=-R.*omega.^2.*sin(omega.*t(n)); %曲柄末端A的加速度

x1=xA/2; y1=yA/2;

x1pp=xApp./2;y1pp=yApp./2; %曲柄质心的加速度

xBpp=xApp-l*(alpha1pp(n).*sin(alpha1(n))+alpha1p(n).^2.*cos(alpha1(n))); yBpp=0; %刀头(滑块)B的加速度

x2=xA+l.*cos(alpha1(n))./2;y2=yA+l.*sin(alpha1(n))./2;

x2pp=xApp-l*(alpha1pp(n).*sin(alpha1(n))+alpha1p(n).^2.*cos(alpha1(n)))./2; y2pp=yApp+l*(alpha1pp(n).*cos(alpha1(n))-alpha1p(n).^2.*sin(alpha1(n)))./2;

%连杆质心的加速度

%对于曲柄m1有方程

% Ox-m1.*x1pp+Ax=0

% Oy-m1.*y1pp+Ay-m1.*g=0

% M-Ay.*xA-Ax.*yA+m1.*g.*x1+m1.*y1pp.*x1+m1.*x1pp*y1=0

%对于连杆m2有方程

% -Ax+Bn.*cos(alpha1)+Bt.*cos(alpha1+pi/2)-m2.*x2pp=0

% -Ay+Bn.*sin(alpha1)+Bt.*sin(alpha1+pi/2)-m2.*y2pp-m2.*g=0

% Bt.*l+m2.*(g+y2pp).*(x2-xA)+m2.*x2pp.*(y2-yA)=0

%对于刀头(滑块)m3有方程

% -f.*N+Bn.*cos(alpha1-pi)-Bt.*cos(alpha1-pi/2)-m3.*xBpp=0

% N-m3.*g+Bn.*sin(alpha1-pi)+Bt.*sin(alpha1-pi/2)=0

% 8个方程，未知数也有Ox,Oy,Ax,Ay,Bn,Bt,N,M共8个

A=[1,0,1,0,0,0,0,0;

0,1,0,1,0,0,0,0;

0,0,-yA,-xA,0,0,0,1;

0,0,-1,0,cos(alpha1(n)),cos(alpha1(n)+pi/2),0,0;

0,0,0,-1,sin(alpha1(n)),sin(alpha1(n)+pi/2),0,0;

0,0,0,0,0,l,0,0;

0,0,0,0,cos(alpha1(n)-pi),-cos(alpha1(n)-pi/2),-f,0;

0,0,0,0,sin(alpha1(n)-pi),sin(alpha1(n)-pi/2),1,0];

B=[m1.*x1pp;

m1.*y1pp+m1.*g;

-m1.*g.*x1-m1.*y1pp.*x1-m1.*x1pp*y1;

m2.*x2pp;

m2.*y2pp+m2.*g;

-m2.*(g+y2pp).*(x2-xA)-m2.*x2pp.*(y2-yA);

m3.*xBpp;

m3.*g];

X=A\B; %X向量的顺序为Ox,Oy,Ax,Ay,Bn,Bt,N,M

XT(n,:)=X‘; %将N转置后以行向量形式存入矩阵XT中

Bn(n)=X(5);Bt(n)=X(6);N(n)=X(7);

%为求滑块最大惯性力Qm从XT中取出的Bn、Bt及N值Q(n)=Bn(n)*cos(alpha1(n)-pi)-Bt(n)*cos(alpha1(n)-pi/2)-f*N(n);

end

for k=1:8

subplot(3,3,k);plot(t,XT(:,k));grid on;

end

subplot(3,3,9);plot(t,Q);title('Q_B-t曲线');grid on;

Qm=max(Q);tm=t(find(Q==Qm));wtm=omega*tm;%找出滑块最大惯性力

Qm及对应的时间t和弧度wt

disp([‘当t=’,num2str(tm),‘即wt=’,num2str(wtm),‘时，滑块B达到最大惯性力

Qm=',num2str(Qm),'。'])