function plot_dwell(ctheta,s,pattern,range) %ctheta = cam angle (deg)--can be a matrix %pattern = denote the type of motion used(a 3 element-row matrix) % 1:uniform 2:parabolic 3:simple harmonic 4: cycloidal % 5:polynomial motion % example [4 3] %range =the degrees the specific motion starts % Output: y is for displacement, yy is the derivative of the displacement with % respect to theta, and yyy the second derivative with respect % to theta. % Example plot_dwell(0:10:360,2,[4 3],[90 180 240]); figure(1);clf; [y,yy,yyy]=dwell(ctheta,range,pattern) h1=plot(ctheta,y*s,'b-',ctheta,yy*s,'k-',ctheta,yyy*s,'r-') legend('Displacement','Velocity','Acceleration',3) xlabel('Elapsed Angle, degrees') grid
function [x,y]=pincam(cth,r0,s,e,L,range,pattern,cw) %Find the pin type cam with an offsect e %Inputs: % cth:angle of cam, degrees % r0:radius of base circle % e:offset % s:stroke % L:length of pin % cw:rotation direction of cam(-counterclockwise,+clockwise %pattern = denote the type of motion used(a 3 element-row matrix) % 1:uniform 2:parabolic 3:simple harmonic 4: cycloidal % 5:polynomial motion % example [4 3] %range =the degrees the specific motion starts, eg.[90 180 240] % Example: [x y]=pincam([10 60],5,2,1,10,[90 180 240],[4 3],-1) figure(1); clf; th=cth*pi/180; s0=sqrt(r0*r0-e*e); for i=1:length(cth) t=th(i)*cw; A=[cos(t) -sin(t);sin(t) cos(t)]; [ym,yy,yyy]=dwell(cth(i),range,pattern); x0=s0+ym*s; Sx=[0 x0 x0+L;e e e]; X=A\Sx; x(i)=X(1,2);y(i)=X(2,2); line(X(1,1:2),X(2,1:2)); line(X(1,2:3),X(2,2:3),'linewidth',3,'color','red') end hold on; plot([0 x],[0 y],'ro',x,y,'k-') axis equal
function [x,y]=pincam3(cth,r0,s,e,L,range,pattern,cw) %Find the pin type cam with an offsect e %Inputs: % cth:angle of cam, degrees % r0:radius of base circle % e:offset % s:stroke % L:length of pin % cw:rotation direction of cam(-counterclockwise,+clockwise %pattern = denote the type of motion used(a 3 element-row matrix) % 1:uniform 2:parabolic 3:simple harmonic 4: cycloidal % 5:polynomial motion % example [4 3] %range =the degrees the specific motion starts, eg.[90 180 240] % Example: [x y]=pincam([10 60],5,2,1,10,[90 180 240],[4 3],-1) figure(4); clf; th=cth*pi/180; s0=sqrt(r0*r0-e*e); for i=1:length(cth) t=th(i)*cw; A=[cos(t) -sin(t);sin(t) cos(t)]; [ym,yy,yyy]=dwell(cth(i),range,pattern); x0=s0+ym*s; Sx=[0 x0 x0+L;e e e]; X=A\Sx; x(i)=X(1,2);y(i)=X(2,2); %line(X(1,1:2),X(2,1:2));
function drawsldpaths(r6,th6,r,th1,td2,tdd2,sigma,npts,driver,mode) clf; figure(1); warning off; r(abs(r)＜eps)=eps; [Qstart, Qstop]=sld_angle_limits(r,th1,driver) npoint=abs(npts); th2=linspace(Qstart,Qstop,npoint); val=zeros(11,npoint); for i=1:npoint, if driver==2, r(1)=th2(i);end [vr b]=sldlink(r,th1,th2(i),td2,tdd2,sigma,driver); [para]=body(r6,th6,vr,3); if mod(i,5)==0|i==1|i==npoint, drawsldlinks(r,th1,th2(i),sigma,driver);
end val(1:3,i)=[vr(1,1)+vr(4,1);vr(2,1);para(2)]; switch driver case 0 val(4:7,i)=[abs(vr(1,1));vr(3,2);vr(3,3);vr(3,4)]; case 1 val(4:7,i)=[abs(vr(1,1));vr(2,2);vr(2,3);vr(2,4)]; case 2 val(4:7,i)=[abs(vr(2,2));vr(3,2);vr(2,3);vr(3,3)]; end val(8:11,i)=[vr(1,5);para(4);vr(4,6);para(5);]; end warning on; plot(val(1,:),'k-','LineWidth',1.5,'linestyle',':');% path of Q plot(val(2,:),'k-','LineWidth',1.5);% path of P plot(val(3,:),'g-','LineWidth',1.5);% path of A axis equal if mode==0, return;end; th2=th2(3:end-3);val=val(:,3:end-3); title0={'Crank Angle','Coupler Angle','Slider Pos'}; title1={'\Theta3(r) & r1(k)', '\Theta2(r) & r1(k)',... '\Theta2(r) & \Theta3(k)' }; title2={'Vel of A (r) & Slider(k)',... 'Acc of A(r) & Slider(k)' }; title3={'\omega(r) & \alpha(b) of Coupler',... '\omega(r) & \alpha(b) of Crank',... '\omega of Crank(r) & Coupler(b)'}; intitle=title0(driver+1); val(abs(val)＞10e+5)=NaN; val(8:11,:)=abs(val(8:11,:)); figure(2); clf; subplot(2,2,1); plot(th2,val(4,:),'k-'); hold on;fact=round(max(val(5,:))/max(val(4,:))*10)/10; plot(th2,val(5,:)/fact,'r-');% crank or coupler angle xlabel(intitle);ylabel(title1(driver+1)); grid on subplot(2,2,2); plot(th2,val(6,:),'r-'); fact=round(max(val(7,:))/max(val(6,:))*10)/10; hold on;plot(th2,val(7,:)/fact,'b-'); xlabel(intitle);ylabel(title3(driver+1)); grid on; subplot(2,2,3); plot(th2,val(8,:),'k-'); hold on;plot(th2,val(9,:),'r-'); xlabel(intitle);ylabel(title2(1)); grid on; subplot(2,2,4); plot(th2,val(10,:),'k-'); hold on;plot(th2,val(11,:),'r-'); xlabel(intitle);ylabel(title2(2)); grid on;
function [data,form] = f4bar(r,theta1,theta2,td2,tdd2,mode,linkdrive) %%function [data,form] = f4bar(r,theta1,theta2,td2,tdd2,mode,linkdrive) % This function analyzes a four-bar linkage when the driving link is % crank or coupler. The input data are: % theta1,theta2 are angles in degrees % r=[r1 r2 r3 r4]= lengths of links(1=frame) %td2 = crank or coupler angular velocity (rad/sec)% tdd2 = crank or coupler angular acceleration (rad/sec^2) %mode = +1 or -1. Identifies assembly mode %linkdrive = 0 for crank as driver; 1 for coupler as driver %data (1:4,1) = link positions for 4 links %data (1:4,2) = link angles in degrees %data (1:4,3) = link angular velocities %data (1:4,4) = link angular accelerations %data (1,5) = velocity of point Q %data (2,5) = velocity of point P %data (3,5) = acceleration of point Q %data (4,5) = acceleration of point P %data (1,6) = position of Q %data (2,6) = position of P %form = assembly status. form = 0, mechanism fails to % assemble. % program revised from Waldron's Textbook % Revised:DSFON, BIME, NTU. Date: Feb. 7, 2007 if nargin<7,linkdrive=0;end mode="1;end" data="zeros(4,6);" linkdrive="=" r="[r(1)" rr="r.*r;d2g=" t1="theta(1);tx=" s1="sin(t1);c1=" sx="sin(tx);cx=" 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;">=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 f orm=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
function [values]=drawlinks(r,th1,th2,mode,linkdrive) %function drawlinks(r,th1,th2,mode,linkdrive) %draw the positions of four-bar links %call f4bar funcion %r: row vector for four links %th1: frame angle %th2: crank angle or couple angle %mode: assembly mode %linkdrive: 0 for crank, 1 for coupler %clf; if nargin<5,linkdrive=0;end mode="1;end" rr="values(:,1);rr(3,1)=" rx="real(rr(:,1));rx(4)=" ry="imag(rr(:,1));ry(4)=" b="=" linkdrive="=" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjs_XAqnlzgujB0onb96zWoaaf3oMVZBRNlALHQ5d_x-tziETX1heySgyYNs2obGpwUqOybYeTqcd5J4jZmHYzXoS2BZTtzOj4Vgo7teXZmLcSJNEM-_9ahABQpdVEyLgkTU-FWQawPIpPu/s1600-h/b94611048-matlabpic-HW8-1-2.jpg"> 第三題 當桿２迴轉時，求出此組四連桿之限制角度，並繪出其位置 function [Ang1, Ang2]=fb_angle_limits(r,theta1,linkdrive) %%function [Ang1, Ang2]=fb_angle_limits(r,theta1,linkdrive) % Find initital & final angles for the driving link % linkdrive= (0 for crank; 1 for coupler as the driver). % Variables: % r=linkage row vector (cm) % theta1=frame angle(degree); % Ang1,Ang2=initial & final angles of the driving link(deg) %Program if nargin<3,linkdrive=0;end theta1="0;end" linkdrive="=" r="[r(1)" r1="r(1);r2=" r3="r(3);r4=" rmin="min(r);rmax=" rtotal="sum(r);Ang1=" ang2="2*pi;">(r3+r4)& abs(r1-r2)(r3+r4)& abs(r1-r2)>=abs(r3-r4) Ang1=acos((r2^2-(r4+r3)^2+r1^2)/(2*r1*r2)); Ang2=-Ang1; end % if (r1+r2)<=(r3+r4)&abs(r1-r2)>=abs(r3-r4) % Ang1=0; % Ang2=2*pi; % end if (r1+r2)<=(r3+r4)&abs(r1-r2)< ang1="acos((r2^2-(r4-r3)^2+r1^2)/(2*r1*r2));" ang2="2*pi-Ang1;" adjst =" (Ang2-Ang1)*1e-8;" ang1="theta1+(Ang1+adjst)*180/pi;" ang2="theta1+(Ang2-adjst)*180/pi;" tt="Ang1;Ang1=" ang2="TT;end" qstart =" 28.9550" qstop =" 331.0450" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjKjNMCI-_92A8-6ernZ2tyBvuqmCOXld5LerL2BNJOzwIH8RfxmxOaFppb9oPdO4BXs7md4y6G_-uYLw7GnzfJcKVm3jZmKy6HCMXrbN6_mDLq3S96o1N3MfQ_NJkDJuz2sd3JRvwGIpwI/s1600-h/b94611048-matlabpic-HW8-1-3.jpg">
第四題 設theta2=[0:20:360]，試繪出此組四連桿之重疊影像，解釋為何有些沒有值 clf; 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
b94611048 張志鵬 本次作業有用到課本裡的function gruebler 及 function grashof function [df]=gruebler(nlink,jointype) %% [df]=gruebler(nlink,jointype) % nlink:no. of total links % jointype:row matrix for number of joints for each type, % the order of elements is: % 1 R-joint 2 slider 3 compound joint(sliding & rolling) % 4 ball 5 cylinder 6 planar 7 cylinder rolling % 8 cam 9 helix 10 ball & 11 point contact % Example: df=gruebler(4,[4]) % Author:D.S.Fon Bime,NTU. Date:Jan. 30, 2007 code=[1 1 2 3 2 3 1 2 1 3 5]; n=length(jointype); dim=3; if n>3, dim=6; end; ff=0; njoint=0; for i=1:n, njoint=njoint+jointype(i); ff=ff+jointype(i)*code(i); end; df=dim*(nlink-njoint-1)+ff;
function ans=grashof(ground_no,linkage) % Function to test the Grashof linkage % Inputs: % ground_no:the ground link number in the order % linkage: row matrix for lengths of the 4 links % in original assigned order. % Example:ans=grashof(4,[4 4.2 2.6 2]) % Revised: March 4, 2006 ground=linkage(ground_no); link=sort(linkage);% sorting the links ig=find(linkage==link(1)); if link(1)+link(4)>link(3)+link(2), ans='Non-Grashof Linkage'; elseif link(1)+link(4)==link(3)+link(2) ans='Neutral Linkage'; elseif link(1)==ground, ans='Double-Crank Linkage'; else switch ig case 1 im=3; case 2 im=4; case 3 im=1; case 4 im=2; end
if ground==linkage(im) ans='Double-Rocker Linkage'; else ans='Crank-Rocker Linkage'; end end
