機械原理大作業(yè)1(連桿機構)27題matlab程序

附錄：程序清單yy=(0:0.1:360);%桿AB的角位移，每隔0.1度計數(shù)yy1=yy/180*pi;%轉化為弧度xb=280*cos(yy1);%點B的x坐標yb=280*sin(yy1);%點B的y坐標w=10;%桿AB的角速度vxb=-w*yb;%點B的速度在x方向的分量vyb=w*xb;%點B的速度在y方向的分量axb=-w*w*xb;%點B的加速度在x方向的分量ayb=-w*w*yb;%點B的加速度在y方向的分量xd=0;%點D的x坐標yd=160;%點D的y坐標vxd=0;%點D的速度在x方向的分量vyd=0;%點D的速度在y方向的分量axd=0;%點D的加速度在x方向的分量ayd=0;%點D的加速度在y方向的分量jbcd=ones(1,3601);%給角BCD賦初值fdb=ones(1,3601);%?給角BD賦初值li=350;%桿BC的長度lj=320;%桿CD的長度lbd=ones(1,3601);%給BD賦初值fi=ones(1,3601);%給桿BC的角位移賦初值fj=ones(1,3601);%給桿CD的角位移賦初值xc=ones(1,3601);%給點Cx坐標賦初值yc=ones(1,3601);%給點Cy坐標賦初值ci=ones(1,3601);%給中間變量賦初值cj=ones(1,3601);%給中間變量賦初值wi=ones(1,3601);%給桿BC的角速度賦初值wj=ones(1,3601);%給桿CD的角速度賦初值ss=ones(1,3601);%給ss賦初值ffg=ones(1,3601);%給構件5的角位移賦初值xg=-25;%點G的x坐標yg=80;%點G的y坐標vxg=0;%點G的速度在x方向的分量vyg=0;%點G的速度在y方向的分量axg=0;%點G的加速度在x方向的分量ayg=0;%點G的加速度在y方向的分量wgf=ones(1,3601);%給桿GF的角速度賦初值%求角BCD，角BDform=1:3601lbd(1,m)=sqrt((xd-xb(1,m))^2+(yd-yb(1,m))^2);if(lbd(1,m)<(li+lj)&&lbd(1,m)>abs(lj-li))jbcd(1,m)=acos((li*li+lbd(1,m)*lbd(1,m)-lj*lj)/(2*li*lbd(1,m)));elseiflbd(1,m)==(li+lj)jbcd(1,m)=0;elseif(lbd(1,m)==abs(lj-li)&&(li>lj))jbcd(1,m)=0;elseif(lbd(1,m)==abs(lj-li)&&(li<lj))jbcd(1,m)=pi;endif(xd>xb(1,m)&&yd>=yb(1,m))fdb(1,m)=atan((yd-yb(1,m))/(xd-xb(1,m)));elseif(xd==xb(1,m)&&yd>yb(1,m))fdb(1,m)=pi/2;elseif(xd<xb(1,m)&&yd>=yb(1,m))fdb(1,m)=atan((yd-yb(1,m))/(xd-xb(1,m)))+pi;elseif(xd==xb(1,m)&&yd<yb(1,m))fdb(1,m)=3*pi/2;elseif(xd>xb(1,m)&&yd<yb(1,m))fdb(1,m)=atan((yd-yb(1,m))/(xd-xb(1,m)))+2*pi;elseif(xd<xb(1,m)&&yd<yb(1,m))fdb(1,m)=atan((yd-yb(1,m))/(xd-xb(1,m)))+pi;endfi(1,m)=fdb(1,m)-jbcd(1,m);%桿BC的角位移iffi(1,m)<0fi(1,m)=fi(1,m)+2*pi;endend%求點C的坐標xc=xb+li*cos(fi);yc=yb+li*sin(fi);forn=1:3601%求桿CD的角位移if(xc(1,n)>xd&&yc(1,n)>=yd)fj(1,n)=atan((yc(1,n)-yd)/(xc(1,n)-xd));elseif(xc(1,n)==xd&&yc(1,n)>yd)fj(1,n)=pi/2;elseif(xc(1,n)<xd&&yc(1,n)>=yd)fj(1,n)=atan((yc(1,n)-yd)/(xc(1,n)-xd))+pi;elseif(xc(1,n)<xd&&yc(1,n)<yd)fj(1,n)=atan((yc(1,n)-yd)/(xc(1,n)-xd))+pi;elseif(xc(1,n)==xd&&yc(1,n)<yd)fj(1,n)=pi/2*3;elseif(xc(1,n)>xd&&yc(1,n)<=yd)fj(1,n)=atan((yc(1,n)-yd)/(xc(1,n)-xd))+2*pi;endendci=li*cos(fi);cj=lj*cos(fj);si=li*sin(fi);sj=lj*sin(fj);g1=ci.*sj-cj.*si;%求桿BC、CD的角速度wi=(cj*vxd-cj.*vxb+sj*vxd-sj.*vyb)./g1;wj=(ci*vxd-ci.*vxb+si*vxd-si.*vyb)./g1;g2=-axb+wi.^2.*ci-wj.^2.*cj;g3=-ayb+wi.^2.*si-wj.^2.*sj;%求桿BC、CD的角加速度ei=(g2.*cj+g3.*sj)./g1;ej=(g2.*ci+g3.*si)./g1;lbf=281.113856;ai=51.499/180*pi;fii=fi+ai*ones(1,3601);%求點F的坐標、速度、加速度xf=xb+lbf*cos(fii);yf=yb+lbf*sin(fii);vxf=vxb-lbf*wi.*sin(fii);vyf=vyb+lbf*wi.*cos(fii);axf=axb-lbf*wi.^2.*cos(fii)-lbf*ei.*sin(fii);ayf=ayb-lbf*wi.^2.*sin(fii)-lbf*ei.*cos(fii);%求桿GF的角位移fori=1:3601ss(1,i)=sqrt((xg-xf(1,i))^2+(yg-yf(1,i))^2);ifxf(1,i)>xg&&yf(1,i)>=ygffg(1,i)=atan((yf(1,i)-yg)/(xf(1,i)-xg));elseifxf(1,i)==xg&&yf(1,i)>ygffg(1,i)=pi/2;elseifxf(1,i)<xg&&yf(1,i)>=ygffg(1,i)=atan((yf(1,i)-yg)/(xf(1,i)-xg))+pi;elseifxf(1,i)<xg&&yf(1,i)<ygffg(1,i)=atan((yf(1,i)-yg)/(xf(1,i)-xg))+pi;elseifxf(1,i)==xg&&yf(1,i)<ygffg(1,i)=3*pi/2;elseifxf(1,i)>xg&&yf(1,i)<ygffg(1,i)=atan((yf(1,i)-yg)/(xf(1,i)-xg))+2*pi;endend%求桿GF的角速度fori=1:3601ifss(1,i)==0wgf(1,i)=0;elsewgf(1,i)=(vyf(1,i)*cos(ffg(1,i))-vxf(1,i)*sin(ffg(1,i)))/ss(1,i);endend%求桿GF的角加速度vss=vxf.*cos(ffg)+vyf.*sin(ffg);egf=(ayf.*cos(ffg)-axf.*sin(ffg)-2*