2007年5月30日 星期三

機動學第十一次作業

第十一次作業(Due date: 12pm May 30, 2007)

請聲明本週(5/24)有來上課。
本人有上本週的課


某凸輪開始時先在0-100∘區間滯留,然後提升後在200至260∘區間滯留,其高度(衝程)為5公分,
其餘l由260∘至360∘則為返程。升程採用等加速度運動,返程之運動型式自定。設刻度區間為10∘,
試繪出其高度、速度及加速度與凸輪迴轉角度間之關係。

下面plot_dwell.m為一執行Dwell函數之範例。其指令如下:

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


plot_dwell(0:10:360,5,[2 2],[100 200 260])

dwell內的parabolicm函數有小錯誤
上面的程式是有修改過的

而結果如下





上面的圖可以看到(紅色圈圈),速度值為斜直線,但是加速度竟是斜直線(應該是水平直線)。這是因為我們取間隔10度做變化,於是連線變為斜直線。

如果我們改成間隔1度的話
就好了



設凸輪之半徑為15公分,以順時針方向旋轉,其從動件為梢型,垂直接觸,長為10公分,從動件之運動
係依照第二項之運動型式。試繪出此凸輪之工作曲線。


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


pincam([0:10:360],15,5,0,10,[100 200 260],[2 2],-1)

突輪與從動件的關係




你能讓此凸輪迴轉嗎?

此題最有挑戰性的
是要找出凸輪的相對位置
然後加以旋轉

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));

%line(X(1,2:3),X(2,2:3),'linewidth',3,'color','red')
%pause(0.5)
%clf

end

for i=1:length(cth)
clf ;
t=th(i)*cw;
B=[cosd(90) -sind(90);sind(90) cosd(90)];
A=[cos(t) -sin(t);sin(t) cos(t)];
[ym,yy,yyy]=dwell(cth(i),range,pattern);
o=x*cos(-t)+y*sin(-t)
p=x*sin(-t)-y*cos(-t)
m=o*cosd(90)+p*sind(90)
n=o*sind(90)-p*cosd(90)
plot([0 m],[0 n],'ro',m,n,'k-')
axis([-50 50 -50 50])
grid on;
x0=s0+ym*s;
Sx=[0 x0 x0+L;e e e];
X=B*Sx;
%line(X(1,1:2),X(2,1:2));
line(X(1,2:3),X(2,2:3),'linewidth',3,'color','red')
pause(0.5)
end
hold on;

grid on


動畫~~~ 成功!

2007年5月26日 星期六

機動學第十次作業

本人有上本週的課


請思考速度與加速度的問題,當一桿以某特定點M等角速度迴轉時,其端點P之速度方向如何?
其加速度方向如何?

在動力學的分析中
P的速度為 w x r
而其方向為切線方向




但是因為這是一個 rotaional motion
所以P點一定有加速度
而加速度為 w x w x r





若該特定點M復以等速水平運動,則同一端點P之速度與加速度方向會變為如何?
若M點同時也有加速度,則點P會有何變化?

如果M點只是水平運動時
那在動力學裡也有學到
此時的P點速度為 (Vm + w x r)
而加速度仍是為w x w x r

如果M點也有加速度時
那P點的加速度為 Ap + w x w x r


若以此推理四連桿的運動,則點P與Q之速度與加速度方向會與桿一(固定桿)之兩端點之關係如何?
與我們前面的作業分析結果有無共通之處?(參看第六章之四連桿機構之運動分析)



固定桿兩端為O點和R點,此兩點由於接地,其速度與加速度皆為0,
此時我們若以第二桿為驅動桿,
則P點速度為
VP=Wop x Rp/o  方向與桿二垂直

AP=[αx Rp/o]+[Wop x (Wop x Rp/o)] 前項為切線加速度與桿2垂直,後項為法線加速度指向O點



Q點的速度與加速度
由於P點有速度
所以Q點的速度為
Vq = Vp + Wpq x Rq/p
但是由於Vq 是硍R點連接 而R點為接地點
所Q點的速度也為 Wq r x Rq/r  方向與桿四垂直

同樣加速度也是
Aq = Ap+α x Rq/p + Wpq x (Wpq x Rq/p)
也等於 α x Rq/r + Wqr x (Wqr x Rq/r) 
切線加速度與桿四垂直,法線加速度指向R點


設有一運動之曲柄滑塊連桿組合,設滑塊之偏置量為零,且在水平方向移動,
試以此機構之曲桿長度及角度,以及連結桿之長度為輸入項,利用matlab寫出一程式計算在不同曲柄角度時,
六點瞬心之對應位置。可順便探討六點瞬心與曲柄角間之關係。

在四連桿的時候
會有六個IC
但是有兩點固定
所以有兩點會在無限遠

現在用老師的程式改一下
用裡面的數據 把IC的置找出來


function [values]=drawsldlinksic(R1,R2,th2,sigma,driver)
clf;
r=[10,R1,R2,0]
th1=0;
AXIS([-5 5 -5 5]);
[values b]=sldlink(r,th1,th2,10,0,sigma,driver);
rr=values(:,1);rr(3)=rr(3)+rr(2);
rx=real(rr);
rx(4)=0;
ry=imag(rr);
ry(4)=0;
t=linspace(pi/2,2.5*pi,361);
C1=0.2*exp(j*t');
x1=real(C1);
y1=imag(C1);
line(x1,y1,'LineWidth',1.5);
C2=0.4*exp(j*t');
x2=real(C2);
y2=imag(C2);
line(x2,y2,'LineWidth',2);
line([x2(91) x2(91)],[y2(91) y2(91)-0.5],'LineWidth',1.5);
line([x2(271) x2(271)],[y2(271) y2(271)-0.5],'LineWidth',1.5);
line([x2(91)-1.5 x2(271)+1.5],[y2(271)-0.5 y2(271)-0.5],'LineWidth',1.5,'linestyle',':');
hold on;
if rx(1)>0,
the2=th2;
else
the2=180-th2;
end;
the3=-values(3,2)
ic13x=rx(1);
ic13y=rx(1)*tand(the2);
ic24x=0;
ic24y=rx(1)*tand(the3);
if b==1
plot([0 rx(1)],[0 0],'k-','LineWidth',5);
hold on;
plot([0 rx(1)],[0 ry(1)],'g-','LineWidth',1.5);
plot([rx(1) ic13x],[ry(1) ic13y],'k:','LineWidth',1.5);
plot([rx(2) ic13x],[ry(2) ic13y],'k:','LineWidth',1.5);
plot([0 ic24x],[0 ic24y],'k:','LineWidth',1.5);
plot([rx(2) ic24x],[ry(2) ic24y],'k:','LineWidth',1.5);

if driver==0
plot([0 rx(2)],[0 ry(2)],'b-','LineWidth',1.5);
plot([rx(2) rx(3)],[ry(2) ry(3)],'r-','LineWidth',2);

else
plot([0 rx(2)],[0 ry(2)],'r-','LineWidth',2);
plot([rx(2) rx(3)],[ry(2) ry(3)],'b-','LineWidth',1.5);
end
plot([rx(1) rx(3)],[ry(1) ry(3)],'k-');
plot(rx,ry,'bo');
plot(ic13x,ic13y,'mo');
plot(ic24x,ic24y,'mo');
text(0,0,' O 12');text(rx(1),ry(1),' R');
text(rx(2),ry(2),' P 23');
text(rx(3),ry(3),' Q 34');
text(ic13x,ic13y,'13');
text(ic24x,ic24y,'24');
length=max(abs(values(2:3,1)));
len=.20*length;ww=.15*length;

[coords] = sldbox(len,ww,rx(3),ry(3),th1);
plot(coords(:,1),coords(:,2),'r-','LineWidth',2);
[coords] = sldbox(len*3,0,rx(3),ry(3)-ww/2,th1);
plot(coords(:,1),coords(:,2),'r:','LineWidth',1.5);
else
fprintf('Combination of links fails at degrees %6.1f\n',th2);
end

axis equal;
grid on;


下面是動畫



我們可以觀一下IC的移動
何時會在上面 何時會在下面
有時會在無限遠哦~~

2007年5月17日 星期四

機動學第九次作業

B94611048 張志鵬

題目
請就教科書中第四章第五節之偏置機構作另類分析,分析過程可採你所知的方式(包括講義中所列的方法)。
運動中分以曲桿驅動及滑塊驅動的方式,並說明運動的界限或範圍。
設此機構之曲桿長Rcm,連桿Lcm,滑塊之偏置量為10cm等數據作分析。
其中,R=10+(學號末二碼),L=R+5 。


r(1) 固定桿
r(2) 曲桿 =58
r(3) 結合(連)桿=63
r(4) 垂直偏置量=10
theta1 第一桿之水平角
theta2 驅動桿之水平夾角
td2 驅動桿角速度
tdd2 驅動桿角加速度
sigma 組合模數
driver 驅動桿


限制角度
老師在網路上有寫出sld_angle_limits
可以做限制角度的分析

function [Qstart, Qstop]=sld_angle_limits(r,theta1,linkdrive)
%
%function [Qstart, Qstop]=sld_angle_limits(r,theta1,linkdrive)%
Find initital & final angles of driving link for slider
% linkdrive= 0 for crank; 1 for coupler as the driver, with
% Qstart & Qstop as the limit angles of driving link.
% linkdrive=2 for slider as input, with Qstart & Qstop as its
% initial & final positions of r1% Variables
% r= linkage row vector (cm)
% Examples:[Qst,Qsp]=sld_angle_limits([1 5 3 1],30,0)
r1=r(1);r2=r(2);r3=r(3);r4=r(4);g2d=180/pi;
switch linkdrive
case 0 %crank
if r3+r4<r2 & r4>=0 & r3>r4
Qstart=asin((r4-r3)/r2);Qstop=asin((r4+r3)/r2);
elseif r3+r4>=r2 & r4>=r3 & r3>=0
Qstart=asin((r4-r3)/r2);Qstop=pi-asin((r4-r3)/r2);
elseif r3-r4<=r2 & r4<0 & r3>=-r4
Qstart=asin((r4-r3)/r2);Qstop=asin((r4+r3)/r2);
elseif r3-r4>=r2 & r3>=0 & -r4>=r3
Qstart=-pi-asin((r4+r3)/r2);Qstop=asin((r4+r3)/r2);
else
Qstart=0;Qstop=2*pi;
end
case 1 %coupler
if r2-r4<=r3 & r4>=0 & r2>=r4
Qstart=asin((r4-r2)/r3);Qstop=pi-asin((r4-r2)/r3);
elseif r2+r4<r3 & r4>=0
Qstart=asin((r4-r2)/r3);Qstop=asin((r4+r2)/r3);
elseif r2+r4<=r3 & r4<=0 & r2+r4>=0
Qstart=-pi-asin((r4+r2)/r3);
Qstop=asin((r4+r2)/r3);
elseif r2-r4<r3 & r4<=0
Qstart=asin((r4-r2)/r3);Qstop=asin((r4+r2)/r3);
else %r2>=(r3+abs(r4))
Qstart=0;Qstop=2*pi;
end
case 2 %slider displacement
Qstart=0;Qstop=0;
arg2=(r2+r3)^2-r4^2;
if abs(r2-r3)>=r4
arg1=(r2-r3)^2-r4^2;
if arg1>0,Qstart=sqrt(arg1);end;
Qstop=sqrt(arg2);
else
if arg2<0, return; end
Qstart=sqrt(arg2);Qstop=-sqrt(arg2);
end
theta1=0;g2d=1;end %case
if Qstop<Qstart,TT=Qstart;Qstart=Qstop;Qstop=TT;end
adjust=(Qstop-Qstart)*1e-8;
Qstart=theta1+(Qstart+adjust)*g2d;
Qstop=theta1+(Qstop-adjust)*g2d;


以第二桿為驅動桿

r=[0 58 63 10];
theta1=0;
linkdrive=0;
sld_angle_limits(r,theta1,linkdrive);


結果如下

Qstart =

-66.0349

Qstop =

246.0349


限制角度位置圖

以第二桿為驅動桿


以第三桿為驅動桿


所以只有在這範圍裡面的值才有效

在了解角度限制後 就可以做動畫了

function [values]=drawsldlinks(r,th1,th2,sigma,driver)
%function drawlinks(r,th1,th2,sigma,driver)
%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
%sigma: assembly mode
%driver: 0 for crank, 1 for coupler
%% Example:
% drawsldlinks([1 2 3 4],0,60,1,0)
%clf;
[values b]=sldlink(r,th1,th2,10,0,sigma,driver)
;rr=values(:,1);rr(3)=rr(3)+rr(2);
rx=real(rr);rx(4)=0;
ry=imag(rr);ry(4)=0;
if b==1
plot([0 rx(1)],[0 0],'k-','LineWidth',4);%ground line
hold on;
plot([0 rx(1)],[0 ry(1)],'g-','LineWidth',1.5);
if driver==0
plot([0 rx(2)],[0 ry(2)],'b-','LineWidth',1.5);
plot([rx(2) rx(3)],[ry(2) ry(3)],'r-','LineWidth',2);
else
plot([0 rx(2)],[0 ry(2)],'r-','LineWidth',2);
plot([rx(2) rx(3)],[ry(2) ry(3)],'b-','LineWidth',1.5);
end
plot([rx(1) rx(3)],[ry(1) ry(3)],'k-');
plot(rx,ry,'bo');
text(0,0,' O');text(rx(1),ry(1),' R');
text(rx(2),ry(2),' P');text(rx(3),ry(3),' Q');
length=max(abs(values(2:3,1)));
len=.20*length;ww=.15*length;
[coords] = sldbox(len,ww,rx(3),ry(3),th1);
plot(coords(:,1),coords(:,2),'r-','LineWidth',2);
[coords] = sldbox(len*3,0,rx(3),ry(3)-ww/2,th1);
plot(coords(:,1),coords(:,2),'r:','LineWidth',1.5);else
fprintf('Combination of links fails at degrees %6.1f\n',th2);
axis equal;
grid on;
end
axis equal;
grid on;



以第二桿為驅動桿
運動動畫


for i=-66.0349:1:246.0349;
drawsldlinks([0 58 63 10],0,i,1,0)
axis([-80 80 -80 125]);
end;




位置分析
利用以上的程式 可以得到各點的位置、速度和加速度

一桿位置 0.0775 一桿角度 0 一桿角速度 0 一桿角加速度 0
二桿位置 0.0290 + 0.0502i 二桿角度 0.0600 二桿角速度 0.0100 二桿角加速度 0
三桿位置 0.0485 - 0.0402i 三桿角度 -0.0397 三桿角速度 -0.0060 三桿角加速度 0.0739
四桿位置 0.0000 + 0.0100i 四桿角度 0.0900 四桿角速度 0 四桿角加速度 0

Q位置向量 0.0290 + 0.0502i q速度 -2.9000 - 5.0229i q加速度-0.1200
P位置向量 0.0775 + 0.0100i p速度 -1.6611 p加速度 0.1800

0.0775 0 0 0 -0.7429
0.0290 + 0.0502i 0.0600 0.0100 0 -1.6611
0.0485 - 0.0402i -0.0397 -0.0060 0.0739 0.0290 + 0.0502i
0.0000 + 0.0100i 0.0900 0 0 0.0775 + 0.0100i


-0.5023 + 0.2900i 0.5800 0.1500
-0.7429 0.7429 0.1800
-2.9000 - 5.0229i 5.8000 -0.1200
-1.6611 1.6611 0.1800

其他的驅動桿的做法也大同小異



路徑分析
在老師的講議中
有一個程式叫drawsldpaths
我看了很久也看不太懂 只能說老師太利害了!!
相當的精彩 值得去了解



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;


drawsldpaths(0,0,[0 58 63 10],0,10,10,1,50,0,1)



2007年5月8日 星期二

第八次作業





B94611048 張志鵬本人4月26日曾來上課


題目
有一組四連桿,其桿長分別為r=[4 3 3 5],由桿2驅動設第一固定桿角度theta1=0度角速度 td2=10rad/s角加速度tdd2=0 rad/s^2

問題一 設桿2角度theta2=45度時,求各點之位置、速度與加速度為何?
下面為f4bar()函數之內容可以幫助我們分析四連桿

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

[val,form]=f4bar([4 3 3 5],0,45,10,0,-1,0)
四連桿之角度為:
第一桿 第二桿 第三桿 第四桿
0 45.0000 69.4856 99.5246
各桿之速度則為:
第一桿 第二桿 第三桿 第四桿
0 10.0000 16.2681 4.9677
各桿之加速度為:
第一桿 第二桿 第三桿 第四桿
0 0 491.4428 383.6120
P點的 位置 速度 加速度
2.1213 + 2.1213i 0.0212 + 0.0212i -0.2121 - 0.2121i
Q點的 位置 速度 加速度
3.1726 + 4.9311i 0.0041 - 0.0245i -1.8712 - 0.4391i

問題二 繪出此四連桿之相關位置及標明各點之速度方向及大小
程式drawlinks之目的即是利用MATLAB繪製四連桿之相關位置。
所以程式本身會呼叫f4bar.m函數以計算四連桿之向量位置,然後加以繪圖。

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="http://4.bp.blogspot.com/_SiVOqZP1-5c/RkMekeidzFI/AAAAAAAAAGY/mVcsi9DJ0lg/s1600-h/b94611048-matlabpic-HW8-1-2.jpg">
第三題 當桿2迴轉時,求出此組四連桿之限制角度,並繪出其位置
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="http://2.bp.blogspot.com/_SiVOqZP1-5c/RkHH-eidzDI/AAAAAAAAAGI/gKsRRewiEag/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

在第三題時 我們求出四連桿之限制角度
從28.955到331.045間才可做旋轉其他角度是不可能發生的



第五題 若將問題三考慮在內,只在可迴轉的範圍內迴轉,請問你能讓此組四連桿作成動畫方式迴轉嗎
下面這程式真的很精彩!!! 可以繪出所有的軌蹪

function move_4paths(r,r6,th6,nlink,th1,td2,tdd2,sigma,driver,ntimes,npts)
%%function move_4paths(r,r6,th6,nlink,th1,td2,tdd2,sigma,driver,ntimes,npts)
%%draw the positions of four-bar links
%call f4bar.m funcion, f4limits.m, fb_angle_limits.m, body.m
%%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=" npoint="abs(npts);" th2="linspace(Qstart,Qstop,npoint);" val="zeros(6,npoint);" i="1:npoint," x="real(val);y=" h="f4limits(r,th1,sigma,driver);" range="1.2*([min(min(x))" i="2:4,set(h(i),'erasemode','xor');end" h0="patch('xdata',[],'ydata',[],'erasemode','xor','facecolor','r',..." i="0;s=" m="1:ntimes" s="-s;" i="0;s=" i="i+s;">npointi==0,break;end;
set(h(2),'xdata',[0 x(2,i)], 'ydata',[ (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


move_4paths([4 3 3 5],0,0,3,0,10,0,1,0,4,150)



下面是我的影片