2007年6月16日

6/16- 第十三次作業(last one)

作業十三 廖婉婷 b94611040



Q1.試設計一組複式齒輪,使其轉速比為125(請說明思考步驟及結果)。

由於每組之轉速比以維持在10以內為佳,超過此值時則需考慮增加齒輪組數。
因此轉速比為125時必須用較多的組數才成達到目標。
首先,將125開方,其值為11.18,仍然比10大,故使用兩組組合仍嫌不足,
如將125開立方,其值為5,比10小,且剛好為整數,故使用三組齒輪之組合。
在此設驅動之小齒輪數最小為12齒,
在此須注意轉速比要求125
轉速比定義為
VR = W(驅動輪) / W(被動輪)

設第一齒輪為驅動輪
也就是
VR = W1/W6 = (N2/N1)*(N4/N3)*(N6/N5)

假設為12:60;12:60;12:60的組合
則VR = (60/12)*(60/12)*(60/12)= 125,與題目要求相符
故其齒數順序為12:60;12:60;12:60三組。

模擬此齒列機構接合狀態
修改十二次作業之draw_2gears函式
以下為修正後內容:
function A=draw_gear_train(Dp,n2,n3,x0,y0,phi)
%修改字作業十二draw_2gears
%Inputs:
%Pd= 徑節
%N1,N2= 兩輪齒數
%x0,y0= 齒輪中心
%phi= 壓力角(deg)
[coords]=draw_gears(Dp,n2,phi,360,x0,y0,0);
[c_ratio,c_length,ad,pc,pb,r2,r3,d2,d3,ag]=contact_ratio(Dp,n2,n3,phi);
%這一行是為了取出d2和d3兩個節圓直徑
x1=x0+d2/2+d3/2;
%兩齒輪節圓半徑合 為了作為第二個齒輪的中心座標
[coords2]=draw_gears(Dp,n3,phi,360,x1,y0,360/n3/2);
%第二個齒輪旋轉一起始角度(半個齒的角度)使得兩齒輪能契合
[coords]=draw_gears(Dp,n2,phi,360,x1,y0,0);
%繪製第二組齒輪,比照以上方法
x2=x1+d2/2+d3/2;
%x座標以累加方式
[coords2]=draw_gears(Dp,n3,phi,360,x2,y0,360/n3/2);
[coords]=draw_gears(Dp,n2,phi,360,x2,y0,0);
%繪製第三組齒輪
x3=x2+d2/2+d3/2;
[coords2]=draw_gears(Dp,n3,phi,360,x3,y0,360/n3/2);

執行程式:
draw_gear_train(8,12,60,0,0,20)
title('b94611040-模擬齒列接合狀態') %標註

得到下圖:



若要模擬齒列動畫
修改作業十二move_two_gears函式
在此須要特別留意各齒輪的速度差異
也就是由繪製起始角k來控制
以下為修正內容:
function A=move_gear_train(Dp,n2,n3,x0,y0,phi,k)
%Inputs:
%Pd= 徑節;
%N1,N2= 兩輪齒數;
%x0,y0= 齒輪中心座標;
%phi= 壓力角(deg);
%新增函數k= 起始繪製角度;
[coords]=move_gears(Dp,n2,phi,360,x0,y0,k*25);
%設起始角為25*k
[c_ratio,c_length,ad,pc,pb,r2,r3,d2,d3,ag]=contact_ratio(Dp,n2,n3,phi);
x1=d2/2+d3/2;
[coords2]=move_gears(Dp,n3,phi,360,x1,y0,360/n3/2-5*k);
%第二個齒輪旋轉一起始角度(半個齒的角度)使得兩齒輪能契合
%又速度為第一齒輪1/5倍且為反向
[coords]=move_gears(Dp,n2,phi,360,x1,y0,360/n3/2-5*k);
x2=x1+d2/2+d3/2;
%第三齒與第二齒共軸,角速度相等且同向
%故k值相等
[coords2]=move_gears(Dp,n3,phi,360,x2,y0,360/n3/2+k);
%第四輪速度為第三輪1/5倍又反向
[coords]=move_gears(Dp,n2,phi,360,x2,y0,360/n3/2+k);
x3=x2+d2/2+d3/2;
%第五齒與第四齒共軸,角速度相等且同向
%故k值相等
[coords2]=move_gears(Dp,n3,phi,360,x3,y0,360/n3/2-k/5);
%第六輪速度為第五輪1/5倍又反向


執行程式:
for k = 0:360;%以迴圈改變繪製起始角度以達動畫效果
title('b94611040-模擬齒列運轉動畫') %標註
pause(0.08) %畫面停留時間
clf; %清除畫面
move_gear_train(8,12,60,0,0,20,k)
end;

得以下動畫:

我們可以看到由左輪開始驅動
且最右側齒輪轉速比最左側齒輪慢很多





Q2.請指出本學期中你自己最感得意的一次作業(請說明其原因,且該作業必須在自己的部落格內)。

最滿意的應該是作業五吧 連結在此
因為第一題就要求繪出整支手的外形
雖然之前作業就有現成的繪製連桿function
但畢竟仍有段落差
故花了一整個週末的下午,就為了繪製手的外形
還記得那時在動工前,還先在小畫家"擬稿"(下圖)

接著遇到第二個障礙就是在第二題
要模擬手指伸展的動作
雖然角度計算好後,function看似順利寫出
但實際跑程式時相對角度卻一直出問題
使得指尖在伸展途中一直會出現很不符合"人體工學"的扭轉
還好最後在作業時限之內終於改正問題
看到手指"順利"展開那刻真的是非常激動


至於第三小題要分析投手在投出球時各指的速度及加速度
一開始完全摸不著頭緒
不太了解題目希望我們作怎樣的分析
不過我還是試著自己去理解程式跑出的圖加以解釋(如下圖)

說實在最後能整理出這些論點我自己也蠻訝異的(從一開始思緒是零)
即使最後的答案可能不完全是老師期望的

總之這次(第五次)作業從週末開始動工
花的時間應該是我在全部裡面數一數二多吧
但最後除了分數外 在blog也得到老師的肯定
這才是最令我開心的!

2007年6月6日

6/6- 第十二次作業

作業十二 廖婉婷 b94611040


*5/31日曾全程來上課*


Q1:一組標準全齒輪齒輪之徑節為8(亦可使用自設值),齒數分別為30T與48T,
其工作壓力角為20度(可為14.5或25度,自選)。
試求其接觸線長度,與接觸比。


這題使用講義上function contact_ratio
此函式之輸入值為徑節、兩齒輪之齒數及壓力角
輸出參數為接觸比,接觸長度,節圓直徑等
須留意此函式中未輸出"基圓直徑"
因此在函式輸出參數中加入D2,D3代表兩齒輪之基圓直徑
其值可由以下式子求出:

基圓直徑 = 齒數x基周節/pi

而其中齒數,基周節皆為函式輸出函數
因此函式改為:
function [c_ratio,c_length,ad,pc,pb,d2,d3,D2,D3,ag]=contact_ratio(pd,n2,n3,phi)
% Inputs:
% Pd= 徑節
% n2,n3= 兩齒輪之齒數
% phi= 壓力角(deg)
% Outputs:
% c_ratio, c_length= 接觸比,接觸長度
% ad= 齒冠
% pc,pb= 周節及基周節
% d2,d3= 兩齒輪節圓直徑
% D2,D3= 兩齒輪基圓直徑(增加之輸出參數)
% ag= 兩齒輪之接近角,遠退角及作用角
% = [alpha2 beta2 theta2 alpha3 beta3 theta3]
d2g=pi/180;
pangle=phi*d2g;
cosx=cos(pangle);sinx=sin(pangle);
ad=1./pd;pc=pi./pd;
pb=pc.*cosx;
r2=n2./(2*pd);
r3=n3./(2*pd);
d2=2*r2;
d3=2*r3;
rb2=r2.*cosx;rb3=r3.*cosx;
D2=n2.*pb./pi;
D3=n3.*pb./pi;
%此處為增加之基圓直徑參數運算
ax=sqrt((r3+ad).^2-(r3.*cosx).^2)-r3.*sinx;
xb=sqrt((r2+ad).^2-(r2.*cosx).^2)-r2.*sinx;
c_length=ax+xb;
c_ratio=c_length./pb;
ag1=[ax./rb2 xb./rb2 c_length./rb2]/d2g;
ag2=[ax./rb3 xb./rb3 c_length./rb3]/d2g;
ag=[ag1;ag2];

由題目給定:
徑節Pd= 8;齒數n2= 30,n3= 48;壓力角phi= 20
>> [c_ratio,c_length,ad,pc,pb,r2,r3,d2,d3,D2,D3,ag]=contact_ratio(8,30,48, 20)

c_ratio =
1.7005

c_length =
0.6275

ad =
0.1250

pc =
0.3927

pb =
0.3690

r2 =
1.8750

r3 =
3

d2 =
3.7500

d3 =
6

D2 =
3.5238

D3 =
5.6382

ag =
10.4850 9.9211 20.4061
6.5532 6.2007 12.7538
得到接觸比 = 1.7005; 接觸長度 = 0.6275




Q2:兩齒輪之節圓、基圓直徑各為如何?請列式計算其結果。

依課本所列式子:
Pd = N/d (徑節=齒數/節徑)
d = N/Pd

又由題目給定徑節Pd = 8;齒數分別為30,48
可以得到d2 = 30/8 = 3.75
d3 = 48/8 = 6

至於基圓直徑則由式:
Pb = pi*D/N (基周節 = pi*基圓直徑/齒數)
D = Pb*N/pi

由程式得到Pb = 0.369
代入式子得D2 = 0.369*30/pi = 3.52
D3 = 0.369*48/pi = 5.64


最後由第一題function contact_ratio驗證之結果:

d2 =
3.7500

d3 =
6

D2 =
3.5238

D3 =
5.6382

發現與計算結果相符!




Q3:此組齒輪是否會產生干涉現象?試列式證明之。


借用講義的圖
當兩輪齒冠圓與作用線交點A,B均落在作用線MN內,無干涉產生
因此要避免干涉
在接近角部分,線段MP應大於AP;退遠角部分NP應大於BP
因此由兩式MP>=AP;NP>=BP
詳細計算如下(參考自課本)



推導出
(N3^2 + 2*N3N2)sin^2φ >= 4 + 4*N2
(N2^2 + 2*N2N3)sin^2φ >= 4 + 4*N3


代入此題得:
左項 = (48^2+2*48*30)sin^2(20) = 606.41
右項 = 4 + 4*30 = 124
左項 > 右項

又由第二式
左項 = (30^2+2*30*48)sin^2(20) = 442.18
右項 = 4 + 4*48 = 196
左項 > 右項

兩式皆成立,我們得到此二齒輪無干涉產生


接著,由講義上isinterf函式作驗證
以下是函式內容:
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

題目給定壓力角phi=20,兩齒輪齒數=30,48
執行程式:
>> isinterf(20,30,48)

ans =

0
無干涉產生,與計算結果相符!




Q4:可否利用draw_gear.m繪出其接合情形,並繪出其動畫效果。

本題將講義中draw_gear函式稍作修改
因為使用draw_gear時發現兩齒輪角度上無法契合
以下為修改後函式 draw_gears:
function [coords]=draw_gears(Dp,N,phi,range,x0,y0,i) 
%利用講義之draw_gear.m作修改
%Dp= 徑節
%N= 齒數
%range= 繪製範圍
%x0,y0= 繪製中心座標
%新增變數i= 齒輪的起始角度
[coord,theta,rp,rb]=tooth(Dp,N,phi);
coords=[];
%去掉原本的i=0;將i值改為自行輸入
%也就是自行設定齒輪繪製起始角度
%以達到讓尺輪旋轉之效果
while i < range
coord1=rotate2D(coord,-i,x0,y0);
coords=[coords;coord1];
i=i+theta;
end
plot(coords(:,1),coords(:,2));hold on;
[coord]=bushing(rp/8,x0,y0);
plot(coord(:,1),coord(:,2),'b-');
[coord]=bushing(-rp,x0,y0);
plot(coord(:,1),coord(:,2),'r:');
[coord]=bushing(-rb,x0,y0);
plot(coord(:,1),coord(:,2),'b:');
axis equal;

接著再另設draw_2gears.m呼叫上面所列函式
以達繪製兩齒輪
以下為函式內容:
function A=draw_2gears(Dp,N1,N2,phi)
%Inputs:
%Pd= 徑節
%N1,N2= 兩輪齒數
%phi= 壓力角(deg)
[coords]=draw_gears(Dp,N1,phi,360,0,0,0);
%設以原點為第一齒輪中心
[c_ratio,c_length,ad,pc,pb,r2,r3,d2,d3,ag]=contact_ratio(Dp,N1,N2,phi);
%為了取d2和d3兩個節圓直徑參數值
x=d2/2+d3/2; %兩齒輪節圓半徑合 為了作為第二個齒輪的中心座標
[coords2]=draw_gears(Dp,N2,phi,360,x,0,360/N2/2);
%第二個齒輪設定中心在(x,0)
%並且旋轉一起始角度(半個齒的角度)使得兩齒輪能契合

執行程式:
draw_2gears(8,30,48,20)
title('兩齒輪接合情形')





如果要繪製動畫
講義之move2_gear函式可輸入函數直接繪製(稍後有執行範例)
不過也可直接由draw_gear函式直接改寫
承接上一題,
以下為改寫內容:
function [coords]=move_gears(Dp,N,phi,range,x0,y0,k) 
%利用老師的draw_gear.m修改成move_gears.m
%Dp= 節徑
%N= 齒輪數
%phi= 壓力角
%range= 繪製範圍
%x0,y0= 齒輪中心座標
%新增變數k= 齒輪的起始繪製角度
[coord,theta,rp,rb]=tooth(Dp,N,phi);
coords=[];
range2 = range+k;
%因起始角度增加,range也要跟著往後移
%不然會出現在360度內齒數越來越少之情況
while k < range2
coord1=rotate2D(coord,k,x0,y0);
%如果將k改為-k;旋轉方向會相反
coords=[coords;coord1];
k=k+theta;
end
plot(coords(:,1),coords(:,2));hold on;
[coord]=bushing(rp/8,x0,y0);
plot(coord(:,1),coord(:,2),'b-');
[coord]=bushing(-rp,x0,y0);
plot(coord(:,1),coord(:,2),'r:');
[coord]=bushing(-rb,x0,y0);
plot(coord(:,1),coord(:,2),'b:');
axis equal;

接著以move_two_gears函式呼叫
以繪製兩個齒輪
function A=move_two_gears(Dp,n2,n3,phi,k)
%Inputs:
%Pd= 徑節;
%N1,N2= 兩輪齒數;
%phi= 壓力角(deg);
%新增函數k= 起始繪製角度;
[coords]=move_gears(Dp,n2,phi,360,0,0,1.6*k);%第一個齒輪設定中心在(0,0)
[c_ratio,c_length,ad,pc,pb,r2,r3,d2,d3,ag]=contact_ratio(Dp,n2,n3,phi);
%這一行是為了取出d2和d3兩個節圓直徑
x=d2/2+d3/2;
%兩齒輪節圓半徑合 為了作為第二個齒輪的中心座標
[coords2]=move_gears(Dp,n3,phi,360,x,-0,360/n3/2-k);
%第二個齒輪設定中心在(x,0)
%並且旋轉一起始角度(半個齒的角度)使得兩齒輪能契合
%需留意此處起始角度加上一負號
%因此兩輪才會是反向旋轉


重點在新增一輸入函數k做為齒輪繪製起點角度
如此,就可以得到跟講義move2_gear.m相似動畫了!

執行程式:
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

執行程式:
move2_gear(8,30,48,20,20)
結果如下:

2007年5月30日

5/30- 第十一次作業

作業十一 廖婉婷 b94611040


*本週(5/24)有來上課*


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

本題我們先採用講義上function dwell
計算出隨角度改變之高度,速度,加速度值
再利用plot_dwell函式繪出關係圖
以下為函式內容:
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

由題目給定開始時先在0-100°區間滯留,然後提升後在200至260°區間滯留,
其餘由260°至360°則為返程。
也就是升程為100~200°,返程為260~360°
注意ctheta表現的方式
[x y]分別表示升程與返程之運動型式
由於題目限定升程為等加速度運動
而返程運動自訂
故分別列出[2 1] [2 2] [2 3] [2 4] [2 5]各項情況之結果
function plot_dwell(ctheta,s,pattern,range)
%ctheta= 凸輪角度
%pattern= 運動型式
% 1:等速運動uniform 2:抛物線parabolic 3:簡諧simple harmonic
%     4:擺線cycloidal 5:多項式polynomial motion
%range= 升程及返程範圍
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

1.升程:拋物線 返程:等速運動
axis([0 360 -180 180])
plot_dwell(0:10:360,8,[2 1],[100 200 260])
title('升程:拋物線 返程:等速運動')


2.升程:拋物線 返程:拋物線
axis([0 360 -180 180])
plot_dwell(0:10:360,8,[2 2],[100 200 260])
title('升程:拋物線 返程:拋物線')


3.升程:拋物線 返程:簡諧運動
axis([0 360 -180 180])
plot_dwell(0:10:360,8,[2 3],[100 200 260])
title('升程:拋物線 返程:簡諧運動')


4.升程:拋物線 返程:擺線運動
axis([0 360 -180 180])
plot_dwell(0:10:360,8,[2 4],[100 200 260])
title('升程:拋物線 返程:擺線運動')


5.升程:拋物線 返程:多項式
axis([0 360 -180 180])
plot_dwell(0:10:360,8,[2 5],[100 200 260])
title('升程:拋物線 返程:多項式')



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

本題使用講義pincam函式
另角度為0~360度 每10度作間隔
半徑r0 = 15
衝程s = 5
由於為垂直接觸 偏置量e = 0
從動件長度L = 10
凸輪順時針轉動故cw = -1
運動型式在此以上題之第一情況為代表
亦即升程為加速度運動 返程為等速運動
因此pattern = [2 1]
function [x,y]=pincam(cth,r0,s,e,L,range,pattern,cw)
%cth: 凸輪角度(deg)
%r0: 凸輪基圓半徑
%e: 偏置量
%s: 從動件衝程
%L: 從動件長度
%cw: 凸輪轉動方向(反時鐘為正,順時鐘為負)
%pattern: 運動型式
% 1:等速運動uniform 2:抛物線parabolic 3:簡諧simple harmonic
%     4:擺線cycloidal 5:多項式polynomial motion
%range: 升程及返程範圍
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]=pincam([0:10:360],15,5,0,10,[100 200 260],[2 1],-1)
title('從動件為梢形之凸輪工作曲線')

以下為結果:


Q:你能讓此凸輪迴轉嗎?

for i = 1:10:360;
[x y]=pincam(i,15,5,0,10,[120 180 210],[2 4],-1)
axis equal;
grid on;
end;

2007年5月23日

5/23- 第十次作業

作業十 b94611040 廖婉婷




*本週(5/17)有來上課*


Q1:請思考速度與加速度的問題,當一桿以某特定點M等角速度迴轉時,其端點P之速度方向如何?
其加速度方向如何?若該特定點M復以等速水平運動,則同一端點P之速度與加速度方向會變為如何?
若M點同時也有加速度,則點P會有何變化?
若以此推理四連桿的運動,則點P與Q之速度與加速度方向會與桿一(固定桿)之兩端點之關係如何?
與我們前面的作業分析結果有無共通之處?


在分析速度,加速度的問題時
我們可以使用之前的function dyad_draw
以下是函式內容:
function dyad_draw(rho,theta,td,tdd)
% Inputs: rho: 連桿長度
% theta: 連桿水平角度(deg)
% td: 角速度(rad/s)
% tdd: 角加速度(rad/s^2)
clf;
[vec,data] = dyad(rho,theta,td,tdd);
%先由function dyad計算出速度 加速度對應值再繪出
x=[0;cumsum(real(data(:,1)))];y=[0;cumsum(imag(data(:,1)))];
for i=1:length(x)-1
linkshape([x(i) y(i)],[x(i+1) y(i+1)],1);
end
for k=1:length(rho)
x0=x(k+1);y0=y(k+1);
vx=x0+real(data(k,2));vy=y0+imag(data(k,2));
ax=x0+real(data(k,3));ay=y0+imag(data(k,3));
line([x0 vx],[y0 vy],'marker','s','linewidth',2);
line([x0 ax],[y0 ay],'marker','s','color','r','linewidth',3)
end
sdata=sum(data);
ss=[real(sdata(1)) imag(sdata(1))];
vv=[real(sdata(2)) imag(sdata(2))];
aa=[real(sdata(3)) imag(sdata(3))];
line([0 ss(1)],[0 ss(2)],'linestyle',':','linewidth',2)
line([ss(1) ss(1)+vv(1)],[ss(2) ss(2)+vv(2)],...
'marker','d','color','g','linewidth',3)
line([ss(1) ss(1)+aa(1)],[ss(2) ss(2)+aa(2)],...
'marker','d','color','m','linewidth',3)

axis equal;grid on

今設一桿以特定點M作等角速度 ω(rad/s)轉動
則另一端P點有速度v = ωr(m/s)
其方向與M至P之向量垂直,與軌跡路徑相切
用於表示p點之速度快慢
至於加速度一般可分為切線及法線加速度
由於本題作等角速度迴轉,故角加速度為零,切線加速度亦為零
因此加速度僅為法線加速度之值: rω^2 (m/s^2), 其方向為向心
今假設桿長4, 水平角度30deg, 角速度0.5rad/s
結果如下圖所示:
dyad_draw(4,30,.5,.2)


若M為等角加速度0.2rad/s^2旋轉
執行程式
dyad_draw(4,30,.5,.2)
gtext ('切線加速度')
gtext ('法線加速度')
gtext ('合加速度')
%標明加速度及其切線,法線分量之位置

應該如下圖所示:


若該特定點M復以等速水平運動,則需考慮到相對運動之觀念
設如果M點是以Vm(m/s)的速度移動
那P點的速度為Vp = Vm + w x r ,
其中w x r為w和r之外積,r為m至p之相對位置
由於M點速度不影響該桿之角加速度 故加速度值不變
Vp和Vm的關係可由下圖作圖解:

若M點同時也有加速度設為A(m/s^2) 則要考慮到相對加速度
由上式相對速度關係式再對時間t微分
可以得到:Ap = Am + (Ap/m)t + (Ap/m)n
其中(Ap/m)t和(Ap/m)n 為相對加速度之切線及法線分量
因此式子變為Ap = Am + αxr - ω^2r
其關係見下圖:



若推至四連桿運動,因為桿一為固定
因此P與Q點可分別看作一端為固定點之桿作迴轉運動,也就是我們一開始討論的情形
可發現與前面作業所討論相同
首先見下圖

我們假設以下情況皆以桿二驅動
由於桿一固定 也就是0點固定點
若r1作等角速度ω迴轉
則p點有一速度Vp = ωr1, 方向與r垂直
因r1沒有角加速度 p點僅有法線加速度(Ap)n

若r1作等角加速度α迴轉
則p點產生了切線加速度(Ap)t
因此加速度合Ap = (Ap)t + (Ap)n

若要分析Q點速度,加速度
由於是第二桿驅動 p點已具有速度/加速度
因此討論到Q點則需採用相對速度/加速度的觀念
也就是VQ = Vp + w x r ; AQ = Ap + αxr - ω^2r





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


根據機構的定義,任二連桿在相同的平面上運動,應可找到一個瞬心
今天一曲柄滑塊機構,可看作一四連桿
又由甘迺迪定理可推得此機構共有4(4-1)/2 = 6個瞬心
首先由於四連桿之旋轉結即為瞬心 我們可以很快找到四個瞬心
另兩個瞬心則是彼此相對之桿造成的

因此要找出另兩個瞬心
由課本上瞬心定位法我們發現:
點13可由12 23之連線與14 34之連線交點獲得
同理
點24則為23 34連線與12 24之連線的交點

由以上找出了六個瞬心位置
便可由function繪出機構瞬心示意圖
以下為函示內容:
function drawcenter(r2,th2,r3)
%輸入函數為
%r2 = 曲桿長度
%th2 = 曲桿水平角度
%r3 = 連結桿長度
theta2 = th2*pi / 180;
%首先須記得將th2換算為弧度
X12 = 0; Y12 = 0;
X23 = r2*cos(theta2); Y23 = r2*sin(theta2);
X34 = X23+r3*cos(Y23/r3); Y34 = 0;
X13 = X34; Y13 = X34*tan(theta2);
X24 = 0; Y24 = X34*tan(Y23/r3);

% 繪出各部份
draw_inner = [X12 Y12;X23 Y23;X34 Y34;X12 Y12];
draw_outer = [X34 Y34;X13 Y13;X23 Y23;X24 Y24;X12 Y12];
draw_block = [X34-0.5 Y34-0.5;X34-0.5 Y34+0.5;X34+0.5 Y34+0.5;X34+0.5 Y34-0.5;X34-0.5 Y34-0.5];

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;

%標註
text(X12+0.2,Y12-0.2,'12');
text(X23+0.2,Y23-0.2,'23');
text(X34+0.2,Y34-0.2,'34');
text(X13+0.2,Y13-0.2,'13');
text(X24+0.2,Y24-0.2,'24');

今假設r2 = 60, th2 = 45, r3 = 50
執行程式
>> axis([0 100 0 100])
axis equal;
drawcenter (50,45,60)

得到以下結果:

須特別留意在曲柄滑塊組合中
點14因為滑塊的關係,與接觸面(視為另一桿)運動方向相同
其方向應與地面垂直
亦即點14存在於無窮遠的位置


若要分析瞬心位置與曲柄角間之關係
我們試著將上一題繪出瞬心之function作角度迴圈製成動畫
再製作動畫前
我們先利用之前作業之function sld_angle_limits
算出給定桿長之極限角度 方便跑動畫時之角度範圍
以下為函式內容:
function [Qstart, Qstop]=sld_angle_limits(r,theta1,linkdrive)
%r= 各桿桿長
%theta1= 第一桿水平角度
%linkdrive= 0 第二桿驅動
%linkdrive= 1 第三桿驅動
%linkdrive= 2 滑塊驅動
%例如本題給定條件r = [0 60 50 0]
%r1隨角度改變因此可設為0
%又假設為滑塊驅動情況 故linkdrive = 2
r1=r(1);r2=r(2);r3=r(3);r4=r(4);
g2d=180/pi;
switch linkdrive
case 0 %crank
if r3+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=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;
%Qstart = 驅動桿(二或三桿)之最低限角度 (deg)
%Qstop = 驅動桿(二或三桿)之最高限角度 (deg)

執行程式:
[Qstart, Qstop]=sld_angle_limits([0 60 10 0],0,2)

Qstart =

10.0000


Qstop =

110.0000

得到給定桿長之曲柄滑塊極限角度約為10~110(deg)

讓程式跑角度迴圈製成動畫之指令:
for i = 10:110,
clf;
axis([0 100 -10 90])
axis equal
%在此必須固定座標軸
%不然畫面會隨著瞬心位置遠近而改變
%很難看出與滑塊機構之關係
drawcenter(60,i,50);
pause(0.1);
end;

結果如下:

2007年5月20日

5/8討論課心得

5/8是採老師帶領小組討論的方式
一組大約八人左右

在這次討論課中
我題出的問題是:
每次的作業老師在講義中都會提供相關程式碼
我們知道那個function input什麼,output什麼
但是如果問到函式內部的程式碼架構
也就是要深入瞭解其中程式語言 就變得不太能掌握
其實主要在想要將程式結果做成自己希望的形式時就有點吃力


不過老師的回答非常有道理:
我們這門課修的是機動學
主要要學的是一些機構的組合和運動方式
matlab在其中只是個輔助工具
我們是要知道如何利用那個程式碼 求得我們需要的數據另外作分析
並不是如何從頭去創造一個完整的程式指令
這也就是每個科目所著重的點不同
如果想要更加瞭解其中組成 或許可往matlab應用的相關課程深入

真的很感謝老師
我覺得這樣的回答 完全解開了我目前的疑惑!

2007年5月16日

5/16- 第九次作業

作業九 廖婉婷 B94611040


*本週(5/3)有來上課*

由課本上定義,滑塊之滑動方向若與水平方向有一距離e時
稱為偏置機構,而其中e稱之偏置量
而此偏置機構可視為一個四連桿結構
其中偏置量將影響滑桿機構之活動範圍
<本題著重在極限位置之分析>

首先依題目給定R=50,L=55,e=10
由課本之slider_limit函式
計算滑塊水平滑行最大距離(即衝程)和界限角度
以下為程式內容:
function [ s theta21,theta22 ] = slider_limit( R,L,e )
%R= 曲桿長
%L= 連桿長
%e= 偏置量
d2g = pi/180;
th1 = asin(e./(R+L));
th2 = asin(e./abs(R-L));
s = (R+L).*cos(th1)-abs(R-L).*cos(th2);
theta21 = th1/d2g;
theta22 = th2/d2g;

執行後得到以下結果:
 slider_limit( 50,55,10 )
abs(s)
abs(theta21)
abs(theta22)

ans =

1.0452e+002 -8.6603e+000i


ans =

104.8809


ans =

5.4650


ans =

117.4463

直接對各項取絕對值
我們得到衝程為104.88
界限角度為5.47和117.45

討論:
在這邊我先發現了個問題
由於篇置量固定
最大衝程應該出現在曲桿和連桿為一直線時
由畢氏定理計算出值應為104.52?






接著利用老師講義上所提供繪製極限位置之函式(drawsldlimits)
此程式除呼叫sld_angle_limits獲得迴轉之範圍外
呼叫drawsldlink函數繪出各桿之位置。
以下為函式內容:
function [Qstart, Qstop]=sld_angle_limits(r,theta1,linkdrive)
%輸入函數:
%r = 四連桿之長度向量
%theta1 = 第一桿之夾角,角度表示(deg)
%linkdrive = 驅動模式(=0 第二桿驅動; =1 第三桿驅動 =2 滑塊驅動)
%輸出(兩個角度)
%Qstart = 驅動桿(第二桿或第三桿)之最低限角度 (deg)
%Qstop = 驅動桿(第二桿或第三桿)之最高限角度 (deg)
r1=r(1);r2=r(2);r3=r(3);r4=r(4);
g2d=180/pi;
switch linkdrive
case 0 %crank
if r3+r4=0 & r3>r4
Qstart=asin((r4-r3)/r2);Qstop=asin((r4+r3)/r2);
elseif r3+r4>=r2 & r3>=r4 & 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=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;

function [values]=drawlinks(r,th1,th2,mode,linkdrive)
%r = 四桿之向量位置
%th1 = 第一桿水平角度
%th2 = 驅動桿水平角度
%mode = +1 or -1 組合模數,負值表閉合型;正值表分支型
%linkdriver = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿)
if nargin<5,linkdrive=0;end
if nargin<4,mode=1;end
[values b]=f4bar(r,th1,th2,0,0,mode,linkdrive);
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;
if b==1
plot([0 rx(1)],[0 ry(1)],'k-','LineWidth',4);
hold on;
if linkdrive==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)],'g-','LineWidth',1.5);
plot(rx,ry,'bo');
text(0,0,' O');text(rx(1),ry(1),' R');
text(rx(2),ry(2),' Q');text(rx(3),ry(3),' P');
else
fprintf('Combination of links fail at degrees %6.1f\n',th2);
end
axis equal
grid on

function drawsldlimits(r,th1,sigma,driver)
%r = 四桿之向量位置
%th1 = 第一桿水平角度
%sigma = +1 or -1 組合模數,負值表閉合型;正值表分支型
%driver = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿) ;2(滑塊驅動)
g2d=180/pi;
[Qstart, Qstop]=sld_angle_limits(r,th1,driver)
if Qstart==0 & Qstop==0, return; end
switch driver
case 0 %Link 2 as driving link
[rr]=drawsldlinks(r,th1,Qstart,sigma,driver);
text(real(rr(2,1)/2),imag(rr(2,1)/2),['s1=' num2str(Qstart,'%6.1f')]);
[rr]=drawsldlinks(r,th1,Qstop,sigma,driver);
text(real(rr(2,1)/2),imag(rr(2,1)/2),['s2=' num2str(Qstop,'%6.1f')]);
case 1 %Link 3 as driving link
[rr]=drawsldlinks(r,th1,Qstart,sigma,driver);
text(real(rr(2,1)+rr(3,1)/2),imag(rr(2,1)+rr(3,1)/2),...
['s1=' num2str(Qstart,'%6.1f')]);
[rr]=drawsldlinks(r,th1,Qstop,sigma,driver);
text(real(rr(2,1)+rr(3,1)/2),imag(rr(2,1)+rr(3,1)/2),...
['s1=' num2str(Qstop,'%6.1f')]);
case 2 %Slide as driving link
if abs(r(2)-r(3)) < r(4);
angx=asin((r(2)+r(3))/r(4))*g2d;
tmax=th1+angx;
tmin=th1+180-angx;
else
tmin=th1+asin((r(2)+r(3))/r(4))*g2d;
tmax=th1+asin((r(2)-r(3))/r(4))*g2d;
end
r(1)=Qstart;offr=0.1*Qstart;
[rr]=drawsldlinks(r,th1,tmin,sigma,driver);
text(real(rr(1,1)/2),imag(rr(1,1)/2)+offr,['r1=' num2str(Qstart,'%6.1f')]);
r(1)=Qstop;
[rr]=drawsldlinks(r,th1,tmax,sigma,driver);
text(real(rr(1,1)/2),imag(rr(1,1)/2)-offr,['r1=' num2str(Qstop,'%6.1f')]);
end
axis equal
grid on


雖然題目只給定R,L,e之值 也就是相當於四連桿之r2,r3,r4
由於前面slider_limit已求得衝程為104.88 也就是相當於四連桿之r1
因此在drawsldlimits函式中
即可直接代入四桿之值

(1)曲桿驅動(以第二桿為驅動桿)
因為曲桿驅動
driver為0
>> axis([-110 110 -40 40])
drawsldlimits([104.88 50 55 10],0,-1,0)

Qstart =

3.6000e-006


Qstop =

360.0000

討論:
在這邊我發現了一個問題
照程式結果跑出來的角度可作360度旋轉
但是根據課本上所述:
"若曲柄須能完成360度的迴轉,則其相關連桿之長度須依其所能夠成三角形關係而定
而若有偏置的情形時,則連結桿與曲柄之長度差應大於偏置距離e"
與本題條件所給R-L<0條件不符,顯然有錯誤
在請教過歐陽太閒同學後
發現是sld_angle_limits程式碼中少了r3>=r4>0這種情況(也就是講義所列第二情形)
因此將原程式碼中:
   elseif r3+r4>=r2 & r4>=r3 & r3>=0
Qstart=asin((r4-r3)/r2);Qstop=pi-asin((r4-r3)/r2);

修正為:
   elseif r3+r4>=r2 & r3>=r4 & r3>=0
Qstart=asin((r4-r3)/r2);Qstop=pi-asin((r4-r3)/r2);


重新跑過一次後:
axis([-110 110 -40 40])
drawsldlimits([104.88 50 55 10],0,-1,0)
title('第二桿驅動極限位置 -b94611040')

Qstart =

-64.1581


Qstop =

244.1581


得到了正確的結果





(2)連桿驅動(以第三桿為驅動桿)
%drive為1
axis([-110 110 -40 40])
drawsldlimits([104.88 50 55 10],0,-1,1)
title('第三桿驅動極限位置 -b94611040')

Qstart =

-46.6582


Qstop =

226.6582






(3)滑塊驅動
%driver為2
>> axis([-110 110 -40 40])
drawsldlimits([104.88 50 55 10],0,-1,2)
title('滑塊驅動極限位置 -b94611040')

Qstart =

-104.5227


Qstop =

104.5227


需特別注意
由於滑塊驅動屬於直線運動
故在程式中之Qstart及Qstop值為第一桿位移之最大與最小量。

討論
做到這邊 一直在思考r1的長度問題
原本是使用slider_limit計算出之最大衝程直接當作r1
但是後來考慮到此為滑塊驅動情況
當為曲稈驅動時會有角度的限制
因此試的代入其他值至r1 發現結果一樣
後來觀察sld_angle_limits程式碼
發現限制角度運算都是使用r2,r3和r4
r1會隨運動角度不同而改變 因此設任何值似乎並不影響



接著分別就第二及第三桿驅動之運動狀態作成動畫:
(1)曲桿驅動(以第二桿為驅動桿)


for i=[-64.1581:244.1581]
title('第二桿驅動極限位置 -b94611040')
pause(0.005);
clf;
drawsldlinks([0 50 55 10],0,i,1,0)
axis([-60 120 -80 80]);
axis equal;
end;


(2)連桿驅動(以第三桿為驅動桿)


for i=[-46.6582:226.6582]
title('第三桿驅動極限位置 -b94611040')
pause(0.005);
clf;
drawsldlinks([0 50 55 10],0,i,1,1)
axis([-100 100 -100 100]);
axis equal;
end;



如果要分析各桿之位置,速度,加速度
我們使用講義中的function sldlinks作分析
以下為程式碼內容:
由前面 已知r2=50, r3=55, r4=10,
r1值隨滑塊運動作變動,在此先設為1
另設驅動桿角速度為5, 角加速度為0
function [data,form] = sldlink(r,theta1,theta2,td2,tdd2,sigma,driver)
% r= [r1 r2 r3 r4],四桿之長
% td2= 驅動桿之角速度(rad/sec)
% tdd2= 驅動桿之角加速度(rad/sec^2)
% sigma= +1 or -1. 組合模數, 負值表示閉合型,正值為分支型
% driver= 0(驅動桿為第二桿); 1 (驅動桿為第三桿) ;2(滑塊驅動)
f=@(num,ndg) round(num*10^ndg)/10^ndg;
data=zeros(4,8);
if driver==1, r(2:3)=[r(3) r(2)];end
rr=r.*r;d2g=pi/180;
theta=[theta1 theta2 0 theta1+90]*d2g;
t1=theta(1);t4=theta(4);
[td,tdd,rd,rdd]=deal(zeros(4,1));
s1=sin(t1);c1=cos(t1);
s4=sin(t4);c4=cos(t4);
switch driver
case 2 % for the sliding block as driver
rd(1)=td2;
rdd(1)=tdd2;
A=-2*r(1)*r(2)*c1-2*r(2)*r(4)*c4;
B=-2*r(1)*r(2)*s1-2*r(2)*r(4)*s4;
C=rr(1)+rr(2)+rr(4)-rr(3)+2*r(1)*r(4)*(c1*c4+s1*s4);
arg=B*B-C*C+A*A;
arg=f(arg,5);
if (arg>=0)
form=1; %assembly flag
t2=2*atan((-B+sigma*sqrt(arg))/(C-A));
s2=sin(t2); c2=cos(t2);
t3=atan2((r(1)*s1+r(4)*s4-r(2)*s2),(r(1)*c1+r(4)*c4-r(2)*c2));
theta(2)=t2; theta(3)=t3;
s3=sin(t3); c3=cos(t3);
%velocity calculation
AM=[-r(2)*s2, -r(3)*s3; r(2)*c2, r(3)*c3];
BM=[rd(1)*c1;rd(1)*s1];
CM=AM\BM;
td(2)=CM(1); td(3)=CM(2);
%acceleration calculation
BM=[r(2)*td(2)*td(2)*c2+r(3)*td(3)*td(3)*c3+rdd(1)*c1;...
r(2)*td(2)*td(2)*s2+r(3)*td(3)*td(3)*s3+rdd(1)*s1];
CM=AM\BM;
tdd(2)=CM(1); tdd(3)=CM(2);
else
form=0; return;
end
case {0,1} % for link 2 or 3 as a crank
td(2)=td2; tdd(2)=tdd2; tx=theta(2);
sx=sin(tx); cx=cos(tx);
% position calculations
A=2*r(4)*(c1*c4+s1*s4)-2*r(2)*(c1*cx+s1*sx);
B=rr(2)+rr(4)-rr(3)-2*r(2)*r(4)*(cx*c4+sx*s4);
arg=A*A-4*B;
arg=f(arg,5);
if (arg>=0)
form=1; %assembly flag
r(1)=(-A+sigma*sqrt(arg))/2;
t3=atan2((r(1)*s1+r(4)*s4-r(2)*sx),(r(1)*c1+r(4)*c4-r(2)*cx));
theta(3)=t3;
s3=sin(t3); c3=cos(t3);
%velocity calculation
AM=[c1, r(3)*s3; s1, -r(3)*c3];
BM=[-r(2)*td(2)*sx;r(2)*td(2)*cx];
CM=AM\BM;
rd(1)=CM(1); td(3)=CM(2);
%acceleration calculation
BM=[-r(2)*tdd(2)*sx-r(2)*td(2)*td(2)*cx-r(3)*td(3)*td(3)*c3;...
r(2)*tdd(2)*cx-r(2)*td(2)*td(2)*sx-r(3)*td(3)*td(3)*s3];
CM=AM\BM;
rdd(1)=CM(1); tdd(3)=CM(2);
else
form=0;
if driver==1,
r(2:3)=[r(3) r(2)];
for j=1:4, data(j,1)=r(j).*exp(i*theta(j));end % position vectors
end
return;
end
if driver==1,
r(2:3)=[r(3) r(2)];
c2=c3;c3=cx;s2=s3;s3=sx;
td(2:3)=[td(3) td(2)];
tdd(2:3)=[tdd(3) tdd(2)];
theta(2:3)=[theta(3) theta(2)];
elseif driver==0
c2=cx;s2=sx;
end
end % end of switch
for j=1:4,
data(j,1)=r(j).*exp(i*theta(j));
data(j,2:4)=[theta(j)/d2g,td(j),tdd(j)];
end
data(:,5)=[rd(1);rdd(1);data(2,1);data(2,1)+data(3,1)];
data(1,6)=i*r(2)*td(2)*exp(i*theta(2));
data(2,6)=data(1,6)+i*r(3)*td(3)*exp(i*theta(3));
data(3,6)=r(2)*(i*tdd(2)-td(2)*td(2))*exp(i*theta(2));
data(4,6)=data(3,6)+r(3)*(i*tdd(3)-td(3)*td(3))*exp(i*theta(3));
%[Vq;Vp;Aq;Ap]
data(1:4,7)=abs(data(1:4,6));data(1:4,8)=angle(data(1:4,6))/d2g;

由預設條件代入函式執行得:
[values,form] = sldlink([0 50 55 10],0,100,5,0,-1,0)
abs(values)

ans =
% column1 column2 column3 column4
0.0472 0 0 0
0.0500 0.1000 0.0050 0
0.0550 0.1345 0.0011 0.0307
0.0100 0.0900 0 0

column5 column6 column7 column8
0.2904 0.2500 0.2500 0.1700
0.9368 0.2904 0.2904 0.1800
0.0500 1.2500 1.2500 0.0800
0.0483 0.9368 0.9368 0.1800

其中第一行為各連桿位置向量。
第二行為各桿之水平夾角,第三及第四行為各桿之角度速度及角加速度。
第五行為第一桿之軸向速度與加速度、P及Q點之位置向量。
第六至八行則為P點與Q點之速度與加速度量,第六行為向量,第七行為絕對量,第八行為夾角。

此為以第二桿為軀動桿 且為預設角度之情況
因此若要以第三桿,連桿驅動, 則可比照此作法執行
另外,在嘗試將r1代入不同值時
我們發現在此處r1的值也不影響執行結果



此外,若要作連續性的速度/加速度分析
則可利用之前講義上的draw_4paths函式
由於在此不需作桿上特定點分析,因此r6,th6之值設為0
function draw_4paths(r6,th6,nlink,r,th1,theta,td2,tdd2,sigma,npts,driver,mode)
%Inputs:
% r= 四連桿之長度,以向量表示
% th1,theta: 第一桿及驅動桿之水平角度
% td2,tdd2: 驅動桿之角速度及角加速度
% sigma: 組合模數
% driver: 0為第二桿驅動, 1為第三桿驅動
% npts: 設定分割的點數或位置
% r6,rh6,nlink: 桿上特定點之位置,包括桿長,與桿之夾角及附於何桿
% mode=0 畫簡單位置圖
% =1 畫所有圖表
% =2 畫簡單位置圖+速度圖
% =3 畫所有圖表,但用簡單位置圖
代入預設值執行程式後:
>> draw_4paths(0,0,3,[0 50 55 10],0,30,2,5,1,100,0,1)


不過由於尚無法完全掌握此函式之內容
因此在瞭解清楚之前 還不敢亂加猜測其中意義

2007年5月6日

5/7- 第八次作業

作業八 廖婉婷 B94611040



本週四 (4/26)曾來上課


Q1.設桿2角度theta2=45度時,求各點之位置、速度與加速度為何?

本題採用講義上之f4bar函式
先分析出四連趕各桿/接點之位置,速度,加速度值 再取其所需
而題目一開始給定之條件為
r=[4 3 3 5], theta1=0, td2=10rad/s, tdd2=0rad/s^2
以下為函式內容:

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




%由題目給定條件執行程式
>> [data,form]=f4bar([4 3 3 5],0,45,10,0,-1,0)

data =

1.0e+003 *

Columns 1 through 5

0.0040 0 0 0 0.0212 + 0.0212i
0.0021 + 0.0021i 0.0450 0.0100 0 0.0041 - 0.0245i
0.0011 + 0.0028i 0.0695 -0.0163 0.4914 -0.2121 - 0.2121i
-0.0008 + 0.0049i 0.0995 -0.0050 0.3836 -1.8712 - 0.4391i

Column 6

0.0021 + 0.0021i
0.0032 + 0.0049i
0
0

%column1: 各桿之位置向量(複數型式)
%column2: 各桿之水平夾角
%column3: 各桿之角速度
%column4: 各桿之角加速度
%column5-1: Q點之速度
%column5-2: P點之速度
%column5-3: Q點之加速度
%column5-4: P點之加速度
%column6-1: Q點之位置向量
%column6-2: P點之位置向量

form =

1

%form=1 表示此四連桿可構成,若為0 則表不行



由程式結果得到
O點位置: 0 + 0i
Q點位置: 0.0021 + 0.0021i
P點位置: 0.0032 + 0.0049i
R點位置: 0.0040 + 0i(同桿一)


>> abs(data(:,5))
%對column5取絕對值
ans =

1.0e+003 *

0.0300
0.0248
0.3000
1.9220


又對column5取絕對值後
Q點速度為0.03
P點速度為0.0248
R點和O點速度為0(因桿1為固定桿)
(速度單位為: 長度/時間)

同理
Q點加速度為0.3
p點加速度為1.9220
R點和O點加速度為0
(加速度為: 長度/時間^2)




Q2.繪出此四連桿之相關位置及標明各點之速度方向及大小(以程式為之)。

本題使用講義上drawlinks函式
繪出由第一題給定條件之四連桿接合情形
drawlinks函式由呼叫第一題f4bar函式,只取其column1的數值
即各桿之位置向量,就能構成基本O,P,Q,R四點,進而完成連桿接合圖
以下為函式內容:
function [values]=drawlinks(r,th1,th2,mode,linkdrive)
%r: 各桿之長度
%th1: 第一桿之水平角
%th2: 驅動桿之水平角
%mode: +1 or -1,組合模式,負值表示閉合型;正值為分支型
%linkdrive: 0是驅動桿為第二桿(crank);1是驅動桿為第三桿(coupler)
if nargin<5,linkdrive=0;end
if nargin<4,mode=1;end
[values b]=f4bar(r,th1,th2,0,0,mode,linkdrive);
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;
if b==1
plot([0 rx(1)],[0 ry(1)],'k-','LineWidth',4);
hold on;
if linkdrive==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)],'g-','LineWidth',1.5);
plot(rx,ry,'bo');
text(0,0,' O');text(rx(1),ry(1),' R');
text(rx(2),ry(2),' Q');text(rx(3),ry(3),' P');
else
fprintf('Combination of links fail at degrees %6.1f\n',th2);
end
axis equal
grid on




>> axis([-1 5 -1 5]);
axis equal;
drawlinks([4 3 3 5],0,45,-1,0)
%在顯示畫面上標明三點速度方向及大小
gtext ('速度方向0.0212 + 0.0212i 大小0.03')
gtext ('速度方向0.0041 - 0.0245i 大小0.0248')
gtext ('速度方向0 + 0 大小0')
gtext ('速度方向0 + 0 大小0')


ans =

Columns 1 through 5

4.0000 0 0 0 0
2.1213 + 2.1213i 45.0000 0 0 0
1.0513 + 2.8098i 69.4856 0 0 0
-0.8274 + 4.9311i 99.5246 0 0 0

Column 6

2.1213 + 2.1213i
3.1726 + 4.9311i
0
0


%Q速度方向0.0212 + 0.0212i 大小0.03
%P速度方向0.0041 - 0.0245i 大小0.0248
%R及0速度方向0 + 0 大小0


<圖> 四連桿之相關位置及各點之速度方向及大小(點圖可放大)



Q3.當桿2迴轉時,求出此組四連桿之限制角度,並繪出其位置(以程式為之)。

此題以講義中drawlimits函式
繪出四連桿之極限位置及限制角度
以下為函式內容:

function drawlimits(r,th1,sigma,driver)
%輸入參數為
%r: 四桿之位置向量
%th1: 第一桿之水平角度
%sigma: +1 or -1 組合模數,負值表閉合型;正值表分支型
%driver: 0表示驅動桿為第二桿; 1表示驅動桿為第三桿
[Qstart, Qstop]=fb_angle_limits(r,th1,driver)
clf;
[values b]=f4bar(r,th1,Qstart,0,0,sigma,driver);
rr=values(:,1);
rr(3)=rr(1)+rr(4);
rx=real(rr);rx(4)=0;
ry=imag(rr);ry(4)=0;
if b==1
plot([0 rx(1)],[0 ry(1)],'k-','LineWidth',4);
hold on;
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)],'-g');
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');
text(rx(2)/2,ry(2)/2,['s1=' num2str(Qstart,'%6.1f')]);
else
fprintf('Combination of links fails at degrees %6.1f\n',Qstart);
end
[values b]=f4bar(r,th1,Qstop,0,0,sigma,driver);
rr=values(:,1);
rr(3)=rr(1)+rr(4);
rx=real(rr);rx(4)=0;
ry=imag(rr);ry(4)=0;
if b==1
if driver==0
plot([0 rx(2)],[0 ry(2)],'b-','LineWidth',1);
plot([rx(2) rx(3)],[ry(2) ry(3)],'r-','LineWidth',1.5);
else
plot([0 rx(2)],[0 ry(2)],'r-','LineWidth',1.5);
plot([rx(2) rx(3)],[ry(2) ry(3)],'b-','LineWidth',1);
end
plot([rx(1) rx(3)],[ry(1) ry(3)],'g-');
plot(rx,ry,'bo');
text(rx(2),ry(2),' p');text(rx(3),ry(3),' q');
text(rx(2)/2,ry(2)/2,[' s2=' num2str(Qstop,'%6.1f')]);
else
fprintf('Combination of links fail at degrees %6.1f\n',Qstop);
end
axis equal
grid on




%同1,2題之條件執行程式
>> drawlimits([4 3 3 5],0,-1,0)

Qstart =

28.9550


Qstop =

331.0450



<圖> 四連桿之極限位置(點圖可放大)

因為r1+r2 < r3+r4
且|r1-r2| < |r3-r4|
故屬於第四類型,角度限制在右側(28.9度和331度)
又由葛拉索準則,最短及最長桿之和大於其他兩桿
故應屬雙曲柄機構





Q4.設theta2=[0:20:360],試繪出此組四連桿之重疊影像,解釋為何有些沒有值。

藉由loop讓程式繪出以20度為桿2移動間隔之四連桿組成方式
以下為程式執行結果:
>>  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


<圖> theta2=[0:20:360],四連桿之重疊影像(點圖可放大)

由程式結果可知theta2為0,20,340,360這四個角度之連桿組合無法顯示
結合問題三我們注意到長度4 3 3 5之連桿組合
其限制角度為28.9及331度,未達及超過該角就無法組合
故本題以20度為間隔 顯示角度只會在40~320度間





Q5.若將問題三考慮在內,只在可迴轉的範圍內迴轉,請問你能讓此組四連桿作成動畫方式迴轉嗎?


首先以講義body函式
作桿上任一點位置,速度,加速度之運算:
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;



接著函式draw4link
與之前drawlinks類似
繪出四連桿連接情形,主動桿使用藍色表示,搖桿為綠色,地桿為黑色
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




%執行程式
%可以看到驅動桿及搖桿之運動軌跡
%還有第三桿上任一點之路徑
move_4paths([4 3 3 5],2,-30,3,0,10,0,1,0,4,100)



另外,
若只需將連桿運動情形製成動畫模式
不需其軌跡路徑,亦可由以下指令執行:
clear;
frames = 200;
LENGTH = [4 3 3 5];
[start, stop] = fb_angle_limits(LENGTH, 0, 0);
for repeat = 1:3,
for i = 0:frames,
clf;
drawlinks(LENGTH, 0, start + (i*((stop-start)/frames)));
axis([-3.5 4.5 -5.5 4]);
pause(0.01);
end;
for i = 0:frames,
clf;
drawlinks(LENGTH, 0, stop - (i*((stop-start)/frames)));
axis([-3.5 4.5 -5.5 4]);
pause(0.008);
end;
end;