function [x]=isinterf(phi,N1,N2) % N1,N2= 兩齒輪齒數 % x=0:無干涉產生 ; x=1:有干涉產生 x=0; sinx=sin(phi*pi/180); if N2 < N1,nn=N1;N1=N2;N2=nn;end if N1*(N1+2*N2)*sinx*sinx<4*(1+N2), x=1; end
for k = 0:360; %以迴圈改變繪製起始角度以達動畫效果 title('b94611040-模擬兩齒輪運轉動畫') %標註 pause(0.08) %畫面停留時間 clf; %清除畫面 move_two_gears(8,30,48,20,k) end;
得到以下結果:
另一方法 直接使用講義之move2_gear.m function move2_gear(Dpitch,nn1,nn2,phi,omega1) % Dpitch= 徑節 % nn1,nn2= 兩輪齒數 % phi= 壓力角(deg) % omega1= 角速度 clf; d2r=pi/180;delt=0.01; [coord1,r1,rb1]=one_tooth(Dpitch,nn1,phi,360,0,0); [coord2,r2,rb2]=one_tooth(Dpitch,nn2,phi,360,0,0); st=180/nn2;if nn1+nn2>2*fix((nn1+nn2)/2),st=0;end coord2=rotate2D(coord2,180+st,0,0); xc1=coord1(:,1);yc1=coord1(:,2); xc2=coord2(:,1);yc2=coord2(:,2); height=max(r1,r2)*1.2; ar=min(abs(r1),abs(r2)); coord=bushing(ar/5,0,0); % Get the coordinates of 1st bushing xb1=coord(:,1)-r1;yb1=coord(:,2); xb2=coord(:,1)+r2;yb2=coord(:,2); coord=bushing(-r1,-r1,0);%Get the 1st pitch circle xp1=coord(:,1);yp1=coord(:,2); coord=bushing(-r2,r2,0);% Get the 2nd pitch circle xp2=coord(:,1);yp2=coord(:,2); plot(xb1,yb1,'r-');hold on; plot(xb2,yb2,'k-'); plot(xp1,yp1,'r:'); plot(xp2,yp2,'k:'); plot([-r1,r2]',[0,0]','r:'); xx1=min([r1,r2])/2;phir=(90-phi)*d2r; plot([0,0]',[-xx1*2,xx1*2]','b:'); plot([-xx1 xx1]',[-xx1*tan(phir), xx1*tan(phir)]','b:');
cir1=line('xdata',[],'ydata',[],'erasemode','xor','linewidth',1,'color','r'); cir2=line('xdata',[],'ydata',[],'erasemode','xor','linewidth',1,'color','k'); line1=line('xdata',[],'ydata',[],'erasemode','xor','linewidth',2,'color','r'); line2=line('xdata',[],'ydata',[],'erasemode','xor','linewidth',2,'color','k'); lx1=[0 -r1]';ly1=[0,0]'; lx2=[r2,0]';ly2=[0,0]'; axis([-2.5*r1 2.5*r2 -height height]); axis equal; title('Press Ctl-C to stop'); theta1=180;theta2=180;s1=omega1*delt/d2r; while 1, z1=rotate2D([xc1,yc1],theta1,-r1,0); z2=rotate2D([xc2,yc2],theta2,r2,0); L1=rotate2D([lx1,ly1],theta1,-r1,0); L2=rotate2D([lx2,ly2],theta2,r2,0); set(cir1,'xdata',z1(:,1),'ydata',z1(:,2)); % For 1st circle moving set(cir2,'xdata',z2(:,1),'ydata',z2(:,2)); % For 2nd circle moving set(line1,'xdata',L1(:,1),'ydata',L1(:,2)); % For 1st line set(line2,'xdata',L2(:,1),'ydata',L2(:,2)); % For 2nd line drawnow; pause(1/s1); %Stop for a while so we can see the graph theta1=theta1+s1; theta2=theta2-s1*r1/r2; if theta1 > 360, theta1=theta1-360;end; %Reverse the direction at bondary line if theta2 > 360,theta2=theta2-360;end; end
function [y,yy,yyy]=dwell(ctheta,range,pattern) %ctheta= 需要計算之凸輪角度 %range= 升程及返程之範圍 %pattern= 運動型式 % 1:等速運動uniform 2:抛物線parabolic 3:簡諧simple harmonic % 4:擺線cycloidal 5:多項式polynomial motion d2r=pi/180; theta=ctheta*d2r;range=range*d2r; dim=length(ctheta); y=zeros(size(ctheta));yy=y;yyy=y; for i=1:dim if theta(i)>=range(3) %for the last motion(downward) mode=pattern(2);betax=2*pi-range(3); switch mode, case 1, [y(i),yy(i),yyy(i)]=uniform(theta(i), range(3),betax,-1); case 2, [y(i),yy(i),yyy(i)]=parabolicm(theta(i), range(3),betax,-1); case 3, [y(i),yy(i),yyy(i)]=harmonicm(theta(i), range(3),betax,-1); case 4, [y(i),yy(i),yyy(i)]=cycloidm(theta(i), range(3),betax,-1); case 5, [y(i),yy(i),yyy(i)]=polynorm(theta(i), range(3),betax,-1); end;
elseif theta(i)>=range(2) % dewell on the top y(i)=1; elseif theta(i)>=range(1) % for the 1st motion(upward) mode=pattern(1);betax=range(2)-range(1); switch mode, case 1, [y(i), yy(i), yyy(i)]=uniform(theta(i), range(1),betax,+1); case 2, [y(i), yy(i), yyy(i)]=parabolicm(theta(i), range(1),betax,+1); case 3, [y(i), yy(i), yyy(i)]=harmonicm(theta(i), range(1),betax,+1); case 4, [y(i), yy(i), yyy(i)]=cycloidm(theta(i), range(1),betax,+1); case 5, [y(i), yy(i), yyy(i)]=polynorm(theta(i), range(1),betax,+1); end end end
hold on; plot(draw_inner(:,1),draw_inner(:,2),'k-'); plot(draw_inner(:,1),draw_inner(:,2),'bo'); plot(draw_outer(:,1),draw_outer(:,2),'r:'); plot(draw_outer(:,1),draw_outer(:,2),'bo'); plot(draw_block(:,1),draw_block(:,2),'k-'); hold off;
function [data,form] = f4bar(r,theta1,theta2,td2,tdd2,mode,linkdrive) %輸入參數為 %r= [r1 r2 r3 r4]: 各桿之長度 %theta1= 桿一之水平角度 %theta2= 驅動桿(可為桿2或3)之水平角度 %td2= 驅動桿之角速度(rad/sec) %tdd2= 驅動桿之角加速度(rad/sec^2) %mode= +1 or -1 組合模數,負值表閉合型;正值表分支型 %linkdrive = 0表示驅動桿為第二桿; 1表示驅動桿為第三桿 if nargin<7,linkdrive=0;end if nargin<6,mode=1;end data=zeros(4,6); % if coupler is the driver, interchange the vetor 3 & 2 if linkdrive==1,r=[r(1) r(3) r(2) r(4)];end rr=r.*r;d2g=pi/180; [theta,td,tdd]=deal(zeros(4,1)); theta(1:2)=[theta1 theta2]*d2g; t1=theta(1);tx=theta(2); s1=sin(t1);c1=cos(t1); sx=sin(tx);cx=cos(tx); % position calculations A=2*r(1)*r(4)*c1-2*r(2)*r(4)*cx; C=rr(1)+rr(2)+rr(4)-rr(3)-2*r(1)*r(2)*(c1*cx+s1*sx); B=2*r(1)*r(4)*s1-2*r(2)*r(4)*sx; pos=B*B-C*C+A*A; if pos>=0, form=1; % Check for the denominator equal to zero if abs(C-A)>=1e-5 t4=2*atan((-B+mode*sqrt(pos))/(C-A)); s4=sin(t4);c4=cos(t4); t3=atan2((r(1)*s1+r(4)*s4-r(2)*sx),(r(1)*c1+r(4)*c4-r(2)*cx)); s3=sin(t3);c3=cos(t3); else % If the denominator is zero, compute theta(3) first A=-2*r(1)*r(3)*c1+2*r(2)*r(3)*cx; B=-2*r(1)*r(3)*s1+2*r(2)*r(3)*sx; C=rr(1)+rr(2)+rr(3)-rr(4)-2*r(1)*r(2)*(c1*cx+s1*sx); pos=B*B-C*C+A*A; if pos>=0, t3=2*atan((-B-mode*sqrt(pos))/(C-A)); s3=sin(t3); c3=cos(t3); t4=atan2((-r(1)*s1+r(3)*s3+r(2)*sx),... (-r(1)*c1+r(3)*c3+r(2)*cx)); s4=sin(t4);c4=cos(t4); end end theta(3)=t3;theta(4)=t4; %velocity calculation td(2)=td2; AM=[-r(3)*s3, r(4)*s4; -r(3)*c3, r(4)*c4]; BM=[r(2)*td(2)*sx;r(2)*td(2)*cx]; CM=AM\BM; td(3)=CM(1);td(4)=CM(2);
%acceleration calculation tdd(2)=tdd2; BM=[r(2)*tdd(2)*sx+r(2)*td(2)*td(2)*cx+r(3)*td(3)*td(3)*c3-... r(4)*td(4)*td(4)*c4;r(2)*tdd(2)*cx-r(2)*td(2)*td(2)*sx-... r(3)*td(3)*td(3)*s3+r(4)*td(4)*td(4)*s4]; CM=AM\BM; tdd(3)=CM(1);tdd(4)=CM(2); %store results in array data % coordinates of P and Q if linkdrive==1, c2=c3;c3=cx;s2=s3;s3=sx; r(2:3)=[r(3) r(2)];theta(2:3)=[theta(3) theta(2)]; td(2:3)=[td(3) td(2)];tdd(2:3)=[tdd(3) tdd(2)]; else c2=cx;s2=sx; end for j=1:4, data(j,1:4)=[r(j)*exp(i*theta(j)) theta(j)/d2g td(j) tdd(j)] ; end % position vectors data(1,5)=r(2)*td(2)*exp(i*theta(2));%velocity for point Q data(2,5)=r(4)*td(4)*exp(i*theta(4));%velocity for point P data(3,5)=r(2)*(i*tdd(2)-td(2)*td(2))*exp(i*theta(2));%acc of Q data(4,5)=r(4)*(i*tdd(4)-td(4)*td(4))*exp(i*theta(4));%acc of P data(1,6)=data(2,1);%position of Q, again data(2,6)=data(1,1)+data(4,1);% position of P
%find the accelerations else form=0; if linkdrive==1, r=[r(1) r(3) r(2) r(4)]; for j=1:4, data(j,1)=r(j).*exp(i*theta(j));end % positions end end
>> for i=0:20:360,drawlinks([4 3 3 5],0,i,-1,0); end
Combination of links fail at degrees 0.0 Combination of links fail at degrees 20.0 Combination of links fail at degrees 340.0 Combination of links fail at degrees 360.0
function [val]=body2(r6,alph,values,link) %輸入參數為: % r6 = length of point from link's origin,Q %alph= angle in deg. between the base link and r %values= vetor output from the call of f4bar.m %link =the link that the point lays on, 1-4 %outputs:(in complex form) %val(1):initial position,(say Q, for link 2) %val(3):middle position,(say P, for link 2) %val(2):final position,or position of A %val(4):velocity of A %val(5):acceleration of A % r2=values(link,1); d2g=pi/180; switch link case 2 % on the crank val(1)=0;v0=0; a0=0; val(3)=values(2,1); case 3 % on the coupler val(1)=values(2,1); val(3)=values(1,1)+values(4,1); v0=values(1,5); a0=values(3,5); case 4 % on the rocker val(1)=values(1,1); val(3)=values(1,1)+values(4,1); v0=0; a0=0; case 1 % on the frame val(1)=0; v0=0; a0=0; val(3)=values(1,1); end th6=(values(link,2)+alph)*d2g; efact=exp(j*th6); w6=values(link,3); a6=values(link,4); val(2)=val(1)+r6*efact; val(4)=v0+j*r6*w6*efact; val(5)=a0+(j*r6*a6-r6*w6*w6)*efact;
function h=draw4link(values,linkdrive) %function draw4link(values,b,linkdrive) %draw the positions of four-bar links % values are data from f4bar.m %linkdrive: 0 for crank, 1 for coupler %clf; if nargin<2, linkdrive=0;end rr=values(:,1); rr(3,1)=rr(1,1)+rr(4,1); rx=real(rr(:,1));rx(4)=0; ry=imag(rr(:,1));ry(4)=0; h(1)=line([0 rx(1)],[0 ry(1)],'color','k','LineWidth',4); h(4)=line([rx(1) rx(3)],[ry(1) ry(3)],'color','g','LineWidth',1.5); if linkdrive==0 h(2)=line([0 rx(2)],[0 ry(2)],'color','b','LineWidth',2); h(3)=line([rx(2) rx(3)],[ry(2) ry(3)],'color','r',... 'LineWidth',1.5,'marker','o'); else h(2)=line([0 rx(2)],[0 ry(2)],'color','r','LineWidth',1.5); h(3)=line([rx(2) rx(3)],[ry(2) ry(3)],'color','b',... 'LineWidth',2,'marker','o'); end text(0,0,' O');text(rx(1),ry(1),' R'); text(rx(2),ry(2),' P');text(rx(3),ry(3),' Q'); axis equal grid on
最後以move_4paths將四連桿運動情形繪製成動畫 做動態模擬
function move_4paths(r,r6,th6,nlink,th1,td2,tdd2,sigma,driver,ntimes,npts) %Inputs: % r: row vector for four links % th1: frame angle % th2: crank angle or couple angle % td2,tdd2:angular velocity and acceleration of the driving link. % sigma: assembly mode % driver: 0 for crank, 1 for coupler % ntimes: no. of cycles % npts: number of points divided % r6,rh6,nlink:additional length and angle for nlink link. %example: % move_4paths([4 2 3 4],2,-30,3,0,10,0,1,0,4,100) % %clf; if nargin<10, ntimes=3;npts=100;end; figure(1); [Qstart, Qstop]=fb_angle_limits(r,th1,driver); npoint=abs(npts); th2=linspace(Qstart,Qstop,npoint); val=zeros(6,npoint); for i=1:npoint, [vr b]=f4bar(r,th1,th2(i),td2,tdd2,sigma,driver); [para]=body2(r6,th6,vr,nlink); val(1:3,i)=[vr(1,1);vr(2,1);vr(1,1)+vr(4,1)]; val(4:6,i)=[para(1);para(3);para(2)]; end x=real(val);y=imag(val); h=drawlimits(r,th1,sigma,driver); line(x(5,:)',y(5,:)','color','r','linestyle',':'); line(x(4,:)',y(4,:)','color','b','linestyle','-.'); line(x(6,:)',y(6,:)','color','k','linestyle',':'); range=1.2*([min(min(x)) max(max(x)) min(min(y)) max(max(y))]); axis(range);axis equal;grid off; for i=2:4,set(h(i),'erasemode','xor');end h0=patch('xdata',[],'ydata',[],'erasemode','xor','facecolor','r',... 'marker','o'); i=0;s=-1;axis off; for m=1:ntimes s=-s; if abs(Qstop-Qstart-360)<1,i=0;s=1;end; while 1, i=i+s; if i>npoint|i==0,break;end; set(h(2),'xdata',[0 x(2,i)], 'ydata',[0 y(2,i)]);%crank set(h(3),'xdata',[x(2,i) x(3,i)], 'ydata',[y(2,i) y(3,i)]);%coupler set(h(4),'xdata',[x(1,i) x(3,i)], 'ydata',[y(1,i) y(3,i)]);%Rocker set(h0,'xdata',[x(4:6,i)], 'ydata',[y(4:6,i)]); drawnow; %flush the draw buffer pause(0.1); end end % for m loop
function h=drawlimits(r,th1,sigma,driver) %function drawlmits(r,th1,sigma,driver) %draw the positions of four-bar links %call f4bar funcion %r: row vector for four links %th1: frame angle %sigma: assembly mode %driver: 0 for crank, 1 for coupler % Example:h=f4limits([4 2 3 4],0,1,0) [Qstart, Qstop]=fb_angle_limits(r,th1,driver) [values b]=f4bar(r,th1,Qstart,0,0,sigma,driver); if b==1, h=drawlinks(values,driver); else fprintf('Combination of links fails at degrees %6.1f\n',Qstart); end [values b]=f4bar(r,th1,Qstop,0,0,sigma,driver); if b==1, h=drawlinks(values,driver); else fprintf('Combination of links fails at degrees %6.1f\n',Qstart); end axis equal grid on
