134第十章(附录)常用MATLAB电算编程10.1连杆机构的电算程序10.1.1如图10.1所示是一铰链四杆机构,试按以下给定的两连架杆对应位置用解析法编写设计四杆机构的程序:1)对应主动构件转角f的三个位置f1、f2、f3,满足从动件转角的三个对应位置p1、p2、p32)对应主动构件转角f的三个位置f1、f2、f3、f4、f5,满足从动件转角的三个对应位置p1、p2、p3、p4、p5图10.1连杆设计(1)图10.2连杆设计(2)10.1.2若已知图10.2机构的最大传动角max和最小传动角min以及尺寸a,d,试编制电算程序求b,c10.1.3若已知图10.2机构的极位夹角f和曲柄长a以及摇杆的摆角p和长度c,试编制电算程序求连杆长b和机架长d10.1.4若已知图10.1机构中连杆在三个位置与X轴方向的夹角和连杆上一点m的三个位置坐标,并且已知固定转动副A,D的坐标,试编制电算程序求连杆上B,C点的一组坐标,从而设计出该机构。10.2凸轮机构的电算程序10.2.1图10.3是滚子直动从动件盘形凸轮在工作中的一个位置,试根据此图编制电算程序求出凸轮轮廓。图10.3滚子直动从动件盘形凸轮13510.2.2试编制求凸轮最大压力角的程序10.3齿轮范成实验演示程序10.3.1编制一个加工齿轮的范成实验演示程序第10章常用MATLAB电算编程题解答与分析10.1连杆机构的电算程序10.1.1解:1)根据教材中的推导,令acbdcaR222221cdR2adR3得:)]()cos[()cos()cos(0003021pfpfppRffRR(10-1)把对应位置的已知条件代入以上算式解方程组即可。为方便起见取f0=p0=0以方便编程。functionlinkage3(d,f1,p1,f2,p2,f3,p3)f1=f1*pi/180;f2=f2*pi/180;f3=f3*pi/180;p1=p1*pi/180;p2=p2*pi/180;p3=p3*pi/180;fc='r1-r2*cos(f)+r3*cos(p)=cos(f-p)';s1=subs(fc,'f',f1);s1=subs(s1,'p',p1);s1=vpa(s1);s2=subs(fc,'f',f2);s2=subs(s2,'p',p2);s2=vpa(s2);s3=subs(fc,'f',f3);s3=subs(s3,'p',p3);s3=vpa(s3);[r1,r2,r3]=solve(s1,s2,s3,'r1','r2','r3');a=d/(r3)c=d/(r2)b=sqrt(a^2+c^2+d^2-2*a*c*r1)例如:若已知图10.1中d=50,a杆转角f为45o、90o、135o对应c杆转角p为52o、82o、112o可调用:linkage3(50,45,52,90,82,135,112)a=27.629285658965426760636076548711c=41.110355468665232376141549605063b=57.2362894665213494750616056630382)对于5个位置要求的情况,在程序中我们设定:对应图10.1,已知d,对应主动件a的转角f1、f2、f3、f4、f5,从动件有转角p1、p2、p3、p4、p5。因不是线性方程,要采用非线形方程解答的指令,即首先要选一组初选向量r:)]5(),4(),3(),2(),1([rrrrrr,各分量依次代表方程(10-1)中的R1、R2、R3和转角初136值f0、p0。计算结果是它们的值,存放在向量y中,)]5(),4(),3(),2(),1([yyyyyy,显然也就是:R1=y(1);R2=y(2);R3=y(3);f0=y(4);p0=y(5)。程序文件如下:functiony=linkage5(r)f1=35*pi/180;p1=5.5*pi/180;;f2=80*pi/180;;p2=34*pi/180;;f3=110*pi/180;;p3=54.2*pi/180;;f4=130*pi/180;;p4=66.8*pi/180;;f5=150*pi/180;;p5=77*pi/180;;y(1)=r(1)-r(2)*cos(f1+r(4))+r(3)*cos(p1+r(5))-cos((f1-p1)+(r(4)-r(5)));y(2)=r(1)-r(2)*cos(f2+r(4))+r(3)*cos(p2+r(5))-cos((f2-p2)+(r(4)-r(5)));y(3)=r(1)-r(2)*cos(f3+r(4))+r(3)*cos(p3+r(5))-cos((f3-p3)+(r(4)-r(5)));y(4)=r(1)-r(2)*cos(f4+r(4))+r(3)*cos(p4+r(5))-cos((f4-p4)+(r(4)-r(5)));y(5)=r(1)-r(2)*cos(f5+r(4))+r(3)*cos(p5+r(5))-cos((f5-p5)+(r(4)-r(5)));y=[y(1)y(2)y(3)y(4)y(5)]例如:图10.1中,已知原动件a和从动件c的五个对应位置为:f1=35o,p1=5.5o;f2=80o,p2=34o;f3=110o,p3=54.2o;f4=130o,p4=66.8o;f5=150o,p5=77o。机架d=50mm调用上述程序方法:r=[0.91.3215*pi/18045*pi/180]fsolve(@linkage5,r,optimset(‘display’,’off’))ans=0.80521.27781.91760.16870.8734继续计算:d=50;c=d/y(2)c=39.1285a=d/y(3)a=26.0738b=sqrt(a^2+c^2+d^2-2*a*c*y(1))b=55.3887f0=y(4)*180/pif0=9.67p0=y(5)*180/pip0=50.067510.1.2解:不难推得:bcdacbf2)(cos222maxbcadcbf2)(cos222min对照图10.2得计算程序:functionlinkage1(fmax,fmin,a,d)fmax=fmax*pi/180;137fmin=fmin*pi/180;s1='b^2+c^2-2*b*c*cos(fmax)=(a+d)^2';s2='b^2+c^2-2*b*c*cos(fmin)=(d-a)^2';s1=subs(s1,'a',a);s1=subs(s1,'d',d);s1=subs(s1,'fmax',fmax);s2=subs(s2,'a',a);s2=subs(s2,'d',d);s2=subs(s2,'fmin',fmin);[b,c]=solve(s1,s2,'b','c');b=vpa(b)c=vpa(c)例如:已知图10.2中最大传动角100o,最小传动角40o,a=30,d=80,求b,c长度调用如下:linkage1(100,40,30,80)b=[66.421270063616906454764670997951][-66.421270063616906454764670997951][76.903878877835238731601215219829][-76.903878877835238731601215219829]c=[76.903878877835238731601215219825][-76.903878877835238731601215219825][66.421270063616906454764670997940][-66.421270063616906454764670997940]舍弃不合理的负号答案。10.1.3解:参照图10.4运用平面几何的推导可知如下结论(为节省篇幅,不再推导):cos1)cos1(2sin2222acb)2sin2sin)(cos(2caba221以f表示极位夹角,k表示行程速比系数,以p表示摇杆摆角,以g2代表角2,以g1代表角1得运算程序:functionlinkage2(a,c,k,p)p=p*pi/180;f=180*((k-1)/(k+1))*pi/180;b=sqrt((2*c^2*(sin(p/2))^2-a^2*(1+cos(f)))/(1-cos(f)));g2=acos((b-a)*sin(f)/(2*c*sin(p/2)));g1=g2-f+p/2;d=sqrt((b-a)^2+c^2-2*(b-a)*c*cos(g1));b138d例如,已知曲柄a长75,摇杆c长290,摆角32o;行程速比系数k为1.25,求连杆和机架长度可调用:linkage2(75,290,1.25,32)b=176.0143d=278.7168图10.4连杆设计(3)10.1.4解:具体算法请详见本书参考文献[1]58页的推导,为节省篇幅这里仅介绍解题算式。对2,3位置都和位置1建立以下算式并求得参数:fyfxxximimmii111101sincosxfyfxfyfxAAiAiAiiiii111011011sincossincosyfyfxfxfyBAiAiAiiiii111011011cossinsincos)(2120120101011yxyyxxciiAiAii把计算结果代入以下算式:011111cyBxAiBiBi共有两个这样的方程,就可求得B点坐标(yxBB11,)同理,把所有以上算式的A换成D就解得C点坐标。以下是已知连杆三个位置与选定的X轴间的夹角f1,f2,f3及连杆上一点m在这三个位置时的坐标(xm1,ym1),(xm2ym2),(xm3,ym3),固定转动副A的坐标(xa,ya),求B点坐标(xb1,yb1)的程序:functionliangan3(f1,f2,f3,xa,ya)139symsxb1;symsyb1;xm1=1;ym1=1;xm2=2;ym2=0.5;xm3=3;ym3=1.5;f12=(f2-f1)*pi/180;f13=(f3-f1)*pi/180;x012=xm2-xm1*cos(f12)+ym1*sin(f12);y012=ym2-xm1*sin(f12)-ym1*cos(f12);x013=xm3-xm1*cos(f13)+ym1*sin(f13);y013=ym3-xm1*sin(f13)-ym1*cos(f13);a12=x012*cos(f12)+y012*sin(f12)-xa*cos(f12)-ya*sin(f12)+xa;b12=y012*cos(f12)-x012*sin(f12)+xa*cos(f12)-ya*sin(f12)+ya;c12=-x012*xa-y012*ya+0.5*(x012^2+y012^2);a13=x013*cos(f13)+y013*sin(f13)-xa*cos(f13)-ya*sin(f13)+0;b13=y013*cos(f13)-x013*sin(f13)+xa*cos(f13)-ya*sin(f13)+0;c13=-x013*xa-y013*ya+0.5*(x013^2+y013^2);s12='a12*xb1+b12*yb1+c12=0';s12=subs(s12,'a12',a12);s12=subs(s12,'b12',b12);s12=subs(s12,'c12',c12)s13='a13*xb1+b13*yb1+c13=0';s13=subs(s13,'a13',a13);s13=subs(s13,'b13',b13);s13=subs(s13,'c13',c13);[xb1,yb1]=solve(s12,s13,xb1,yb1);xb1yb110.2凸轮机