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');gridon;subplot(1,2,2);plot(t,ab);title('曲柄滑块机构的滑块a-t图');xlabel('时间t(曲柄旋转一周)');ylabel('滑块加速度a');gridon;%下面黄金分割法求滑块的速度与加速度最大值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);whileabs(a-b)=epsilonify1=y2b=a2;a2=a1;y2=y1;a1=b-0.618*(b-a);y1=-R.*omiga.*sin(omiga.*a1)+R.*omiga.*cos(omiga.*a1).*tan(alpha1);elsea=a1;a1=a2;y1=y2;a2=a+0.618*(b-a);y2=-R.*omiga.*sin(omiga.*a2)+R.*omiga.*cos(omiga.*a2).*tan(alpha1);endn1=n1+1;endvbm1=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));whileabs(a-b)=epsilonify1=y2b=a2;a2=a1;y2=y1;a1=b-0.618*(b-a);y1=-(-R.*omiga.*sin(omiga.*a1)+R.*omiga.*cos(omiga.*a1).*tan(alpha1));elsea=a1;a1=a2;y1=y2;a2=a+0.618*(b-a);y2=-(-R.*omiga.*sin(omiga.*a2)+R.*omiga.*cos(omiga.*a2).*tan(alpha1));endn1=n1+1;endvbm2=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));whileabs(a-b)=epsilonify1=y2b=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));elsea=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));endn1=n1+1;enddisp(['经过',num2str(n1),'次计算,找到的加速度最大值对应的w是:',num2str(omiga*(a+b)/2)])disp(['速度最大值点与加速度最大值相差的wt是:',num2str((omiga*(a+b)/2)-vbm1),'弧度。'])曲柄连杆机构动力学分析图clear;clc;closeall;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值forn=1:100alpha1(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,MXT(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);endfork=1:8subplot(3,3,k);plot(t,XT(:,k));gridon;endsubplot(3,3,9);plot(t,Q);title('Q_B-t曲线');gridon;Qm=max(Q);tm=t(find(Q==Qm));wtm=omega*tm;%找出滑块最大惯性力Qm及对应的时间t和弧度wtdisp([‘当t=’,num2str(tm),‘即wt=’,num2str(wtm),‘时，滑块B达到最大惯性力Qm=',num2str(Qm),'。'])