1、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(

2、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;%下面黃金分割法求滑塊的速度與加速度最大值epsi

3、lon=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=epsilon if y1=epsilo

4、n 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.*omig

5、a.2.*sin(omiga.*a2).*tan(alpha1); end n1=n1+1;enddisp(經過,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

6、=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.*

7、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.*s

8、in(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(

9、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,

10、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+

11、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=AB; %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);endfor k=1:8 subplot(3,3,k);plot(t,XT(:,k);grid on;endsubpl