2007年6月15日 星期五

機動學 作業十三

作業十三 b94611042 王志豪 (Due date: 12pm June 16,2007)

我的blog

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

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

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

由題目可知,我們欲求之轉速比為125。
首先,從講義或是課本中,可知道轉速比之公式:





轉速比公式
註:算出來答案的正負號代表著平面第一個齒輪和最後一個齒輪轉的方向,其中,正為相同,負為相反。

至於普通串聯齒輪的轉速比,其推導則為:













普通串聯齒輪轉速比推導

然而,題目所要求的為複式齒輪,因為有共軸齒輪,故須修正!













修正版公式推導
註:當中打「'」者為共軸齒輪之第二個齒輪。

欲求得最終轉速比......
是利用主動輪和驅動輪透過中間許多組齒列......
將各組減速比相乘而得。

按題目要求,轉速比要為1:125(若為放大)。
由於任兩齒輪間之轉速比以不超過10為原則,考慮:

125^(1/2)=11.18
125^(1/3)=5

由於125開平方後得到約11.18......
意即需兩組轉速比為11.18齒列作組合,不和我們的要求。
故此設計不適合僅用兩組齒列。
而將125開三次方後則是得到5這個答案,比10小。
符合要求下,因此我們可以設計三組轉速比為5的齒列來組合。

也就是使用組合串聯三組!
可使用1:5的齒輪組串連三次,就可得到125的轉速比。

整個組合的齒數順序為12:60;12:60;12:60
轉速比為(60/12)*(60/12)*(60/12)= 125
誤差值為0 %

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

這個嘛......其實這些作業裡面,我並不會有哪一次特別得意耶。說句實在話,每一次的作業,學到的東西都很多,我個人其實不是那種非常厲害可以瞬間就舉一反三的那種人物,每一次的作業我時常都問題一大堆,得叼擾同學們,不斷請教對方,這才逐漸能一步步地將作業打造成型,完成每一週的任務。在此我非常感謝那些曾經幫我解惑過的同學們。在一次次的作業中成長,真的學到了不少東西。

當然,如果真的真的要嚴格講起來最得意的作品,應該是作業八之後的作品吧。在此之前,不少地方一直抓不住要領,以作業八為分隔線,終於比較能切中要點,在各方面的資料彙整以及程式方面也完善很多。唯一比較缺憾的是,作業十個人費了不少心思做了許多的動畫,然而由於第一版的解釋不夠完善,以致於造成教授在批覽時的困惑,在此感到十分抱歉。其他的,我個人都還蠻滿意的啦,雖然花了很多時間,但我有一些外系的朋友看過後都覺得蠻不錯的,覺得很神奇(雖然他們說看不懂程式.....),辛苦的成績被肯定,還是挺高興的。

大致就是這樣吧,機動學在此告一段落,他會是我大學生涯中一段值得回首的記憶。

2007年6月8日 星期五

機動學 作業十二

作業十二 b94611042 王志豪 (Due date: 12pm June 6,2007)
我5/31日曾全程來上課

我的blog

題目:
1. 請聲明5/31日曾全程來上課。
2. 一組標準全齒輪齒輪之徑節為8(亦可使用自設值),齒數分別為30T與48T,其工作壓力角為20度(可為14.5或25度,自選)。

* 試求其接觸線長度,與接觸比。
* 兩齒輪之節圓、基圓直徑各為如何?請列式計算其結果。
* 此組齒輪是否會產生干涉現象?試列式證明之。
* 可否利用draw_gear.m繪出其接合情形,並繪出其動畫效果。

12.1
我5/31日曾全程來上課



12.2
==>一組標準全齒輪齒輪之徑節為8(亦可使用自設值),齒數分別為30T與48T,其工作壓力角為20度(可為14.5或25度,自選)。

* 試求其接觸線長度,與接觸比。
* 兩齒輪之節圓、基圓直徑各為如何?請列式計算其結果。
* 此組齒輪是否會產生干涉現象?試列式證明之。
* 可否利用draw_gear.m繪出其接合情形,並繪出其動畫效果。
Ans:

1.
假設齒輪徑節為8,齒數為30T跟48T,工作壓力角為20度.
此時,我們可以利用教學網站上的一個function來解決這個問題.
也就是function [c_ratio, c_length, ad, pc, pb, d2, d3, ag] = contact_ratio(pd, n2,n3, phi)

其中輸入的參數分別如下:
Pd:徑節
n2, n3:兩齒輪之齒數
phi:壓力角

輸出參數:
cr_ratio:接觸比
cr_length:接觸長度
ad:齒冠
pc, pb:周節及基周節
d2, d3:兩齒輪節圓直徑。
ag:兩齒輪之接近角、遠退角及作用角

************************************
function [c_ratio,c_length,ad,pc,pb,d2,d3,ag]=contact_ratio(pd,n2,n3, phi)
%
%Find the contact ratios
% Inputs:
% Pd: Diametrial pitch;
% n2,n4:number of both gears;
% phi: pressure angle, degrees
% Outputs:
% c_ratio, c_length: contact ratio and length
% ad:addendium
%   pc,pb: circular and basic circular pitches
%   r2, r3: radii of pitch circles
%  ag: angles of action, in matrix of
% [alpha2 beta2 theta2 alpha3 beta3 theta3]
% Example: [c_r,c_l,ad,pc,pb,d2,d3,ag]
% =contact_ratio(6,24,48,20)
% Revised: March 9, 2006
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;
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];
**************************************

主程式
==>
function [c_ratio, c_length, ad, pc, pb, d2, d3, ag] =contact_ratio(8,30,48,20)

結果如下:

接觸比 =1.7005
接觸長度 =0.6275
齒冠 =0.1250
周節=0.3927
基周節=0.3690
齒輪一節圓直徑=3.7500
齒輪二節圓直徑=6
兩齒輪的接近角、遠退角及作用角
==>(1)10.4850 9.9211 20.4061
(2)6.5532 6.2007 12.7538

因此
==>
*接觸比:c_ratio = 1.7005
*接觸長度:c_length =0.6275

所謂接觸比,即同時有幾對齒相互嚙合,數值越大力量分散較平均,因此最好大於1,這樣比較可以避免有未嚙合的狀況,因此,這個狀況是符合的.

若欲手動計算,可用課本上的公式
==>
接觸比 mc = 接觸路徑長度/基周節




2.
節圓
==>
由公式: Pd=N/D (徑節=齒數/節圓直徑)

30T的齒輪節圓直徑: 8=30/D D=3.75 (cm)
48T的齒輪節圓直徑: 8=48/D D=6 (cm)

和第一小題用程式執行的結果一樣.

基圓
==>
利用課本的公式,將節圓直徑乘上工作壓力角的餘弦值即是答案:

30T的齒輪基圓直徑: 3.75*cos(20)=3.524 (cm)
48T的齒輪基圓直徑: 6*cos(20)=5.638 (cm)


3.
關於這個,可以分為接近角以及退遠角兩方面來討論.

接近角
==>
必須滿足MP>=AP
條件可以變成(N2^2+2*N2*N3)sin(壓力角)^2>=4+4*N3
是以N2=30,N3=48,壓力角=20度,代入.
(900+2880)*sin(20)^2=442.176 > 4+4*N3=196
結論,不會干涉.

退遠角
==>
必須滿足NP>=BP
條件可以變成(N3^2+2*N3*N2)sin(壓力角)^2>=4+4*N2
(2304+2880)sin(20)^2=606.41 > 4+4*N2=124
結論,不會干涉.

所以總結,不會產生干涉現象!

以下即為齒輪干涉的概念圖==>














另外,假如不想要每次那麼麻煩計算.
網路講義中有提供了一個程式可供使用判斷呦.
這個程式就是function [x]=isinterf(phi,N1,N2)

***********************************************
function [x]=isinterf(phi,N1,N2)
%
% Test if the gear set exists an interference
% phi:pressure angle, in degrees
% N1,N2:teeth of both gears
% x=0:no interference; x=1 interence exists
x=0;
sinx=sin(phi*pi/180);
if N2if N1*(N1+2*N2)*sinx*sinx<4*(1+N2), x=1;end
************************************************
phi為壓力角
N1為齒輪一齒數
N2為齒輪二齒數
x=1表示產生干涉
x=0表示不產生干涉

主程式
==>
function [x]=isinterf(20,30,48)

ans=0
所以是不產生干涉,與計算結果符合的!


4.
function [coords]=draw_gear(Dp,N,phi,range,x0,y0)

關於這個函式,它的基本參數定義如下:
Dp: 節矩
N: 齒數
phi: 壓力角
range: 繪出之部份
x0,y0: 齒輪中心座標

*************************************************
function [coords]=draw_gear(Dp,N,phi,range,x0,y0)
% [coords]=draw_gear(Dp,N,phi,range,x0,y0)
% To draw a whole gear
% Inputs:
% Dp: Diametrical pitch
% N: no of teeth in a gear
% phi: pressure angle, degrees
% range: the section range to be drawn
% x0,y0: the location of the gear center
% Example [coords]=draw_gear(10,15,20,360,0,0)
[coord,theta,rp,rb]=tooth(Dp,N,phi);
coords=[];i=0;
while i 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;
*************************************************

[coords]=draw_gear(8,48,20,360,9.75/2,0,12.7538)

圖形如下==>















由於我們的目標是期望它能繪出動畫效果.
而這個程式本身只是會製齒輪外型使其成為一齒輪組.
因此,勢必要加入變數theta,使其擁有轉動的動畫效果.
方式是將程式裡頭的while迴圈中......
把coord1=rotate2D(coord,-i,x0,y0)中的-i改成-i+theta2.
一開始可供輸入的變數則變為如下:

function [coords]=draw_gear(Dp,N,phi,range,x0,y0,theta2)

改造程式完成後,便可開始進行動畫繪製工作.
在此,利用第一小題程式求得的作用角~~
以及第二小題求得的節圓直徑.
將這兩變數代入程式後,多個for迴圈.
如此一來,便大功告成!

主程式
==>
for theta=0:0.5:24
clf
[coords]=draw_gear(8,30,20,360,0,0,20.4061+theta)
[coords]=draw_gear(8,48,20,360,9.75/2,0, 12.7538-theta*.625)
pause(0.2)
end
axis equal




後來在看講義的時候,發現有個好用的function
叫做move2_gear(Dpitch,nn1,nn2,phi,omega1)
做出來的效果很好......
因此,在此弄上這個漂亮許多的版本哩~~
(p.s.1上面那個自製似乎轉動速度慢了點...=.=)
(p.s.2 ctrl+C 是好物......good!)

***************************************************
function move2_gear(Dpitch,nn1,nn2,phi,omega1)
% move2_gear(Dpitch,nn1,nn2,phi,omega1)
% To draw a whole gear
% Inputs:
% Inputs:
% Dpitch:dimetral pitch
% nn1,nn2: no. of teeth for both gears
% phi:pressure angle, degrees
% omega1: angular velocity of gear 1
% Example move2_gear(10,15,20,20,10)
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
******************************************
Dpitch:節矩
nn1,nn2:兩齒輪之齒數
phi:壓力角,degrees
omega1: 齒輪1之角速度,rad/s


主程式
==>
move2_gear(8,30,48,20,10)



2007年6月2日 星期六

機動學 作業十一

作業十一 b94611042 王志豪 (Due date: 12pm May 30, 2007)
我有上本週5/24的課

我的blog

題目:
1. 請聲明本週(5/24)有來上課。

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

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

4. 你能讓此凸輪迴轉嗎?

11.1
我有上本週5/24的課



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

關於這題,基本上,我們可以有兩種作法.

分別是可以利用講義中的function [y, yy, yyy]=parabol_cam(phi, phi_in, beta_range, direct, travel,rpm),或者,使用另外一個程式function plot_dwell(ctheta,s,pattern,range),這兩個程式同樣都可以解決這個問題.

關於parabol_cam的程式碼,我就不多詳述了.(貼上blog的時候,似乎系統對裡面一些字有意見,認為語法有問題,所以就不再貼上,詳細可見於講義中.),由於都是做差不多的事情,為了避免在分析時造成混淆,個人決定不採用parabol_cam,而是使用更加方便的plot_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
*********************************************************************************

關於這個function plot_dwell,我簡單說明一下:

ctheta:代表我們打算繪製哪個角度範圍內的圖形.
s:衝程.
pattern:代表升程以及返程時的運動方式.
range:代表的是升程是從哪個角度開始還有結束,以及返程是從哪裡開始.

有了這些條件,我們便可以開始行動解決問題了!!!

從題目所給的條件看來
第一項表示的是角度,故選擇0-360度之區間做圖,輸入0:10:360
第二項則是動件衝程:題目設定為5cm,所以輸入5
第三項表示形式,由於為等加速度運動,故寫入2
第四項表示運動始末的角度,故前兩項各為100.200,返程始於260度,因此第三項為360

而關於模式方面,由於返程運動模式自訂,故會分成五種狀況:
1:等速運動uniform
2:抛物線parabolic
3:簡諧simple harmonic
4:擺線cycloidal 
5:多項式polynomial motion

主程式1(返程等速度運動)
==>
plot_dwell(0:10:360,5,[2 1],[100 200 260])















等速度返程的運動情況


主程式2(返程等加速度運動,也就是拋物線!)
==>
plot_dwell(0:10:360,5,[2 2],[100 200 260])















等加速度返程的運動情況


主程式3(返程簡諧運動)
==>
plot_dwell(0:10:360,5,[2 3],[100 200 260])















簡諧運動返程的運動情況



主程式4(返程擺線運動)
==>
plot_dwell(0:10:360,5,[2 4],[100 200 260])















擺線運動返程的運動情況


主程式5(返程多項式運動)
==>
plot_dwell(0:10:360,5,[2 5],[100 200 260])















多項式運動返程的運動情況



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

由於題目說的是梢形元件,所以我們可以直接套用講義上的一個function [x,y]=pincam(cth,r0,s,e,L,range,pattern,cw),利用這個function來幫助我們解決問題.

******************************************************************
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 pincam,我簡單說明一下:
cth:凸輪角度,度數
r0:凸輪基圓半徑
e:偏置量
s:從動件衝程
L:從動件長度
cw:凸輪轉動方向(反時鐘為正,順時鐘為負)
pattern=運動的型式,二元素之列矩陣.
range=升程及返程之範圍,三元素列矩陣

從題目所給的條件看來
第一項表示角度,選擇0-360度之區間做圖,故輸入0:10:360
第二項則是表示突輪基圓半徑,題目為15cm,所以寫入15
第三項表示從動件衝程,題目設為5cm,因此寫5
第四項乃偏置量,因為題目沒有要求,所以設為0
第五項是從動件長度,題目為10cm,故輸入10
第六項表示運動始末的角度,故前兩項各為100.200,返程始於260度,故第三項為360
第七項為形式,因為皆為等加速度運動,所以選擇2
第八項表示凸輪轉動方向,題目所求為順時針,所以很明顯是-1

故主程式
==>
[x y]=pincam([0:10:360],15,5,0,10,[100 200 260],[2 2],-1);
















11.4
==>你能讓此凸輪迴轉嗎?
Ans:

在此我們用到11.3的function pincam,把這個程式作修改即可.
首先對裡頭一開始的凸輪繪製保留.
然後原先line出圖形的地方,刪除,以利於動畫的執行.
最後,把裡面在多加入一個旋轉迴圈,使這個凸輪迴轉就大功告成了!

改造後的新function
==>
function [x,y]=pincam2(cth,r0,s,e,L,range,pattern,cw)
%Inputs:
%cth:凸輪的角度
%r0:基圓半徑
%e:偏置量
%s:衝程
%L:梢長
%cw:順時或逆時
%range:運動範圍
figure(1);
%先叫出空白圖
pause(2);
%停兩秒,以利準備拍攝的程式
clf;
th=cth*pi/180;
s0=sqrt(r0*r0-e*e);
oo=37;
for nn=1:10:360 mm=nn*pi/180;
oo=oo-1;
clf;
axis([-30 30 -30 30])
for i=1:length(cth)t=th(i)*cw;A=[cos(t+mm) -sin(t+mm);sin(t+mm) cos(t+mm)];
[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);
end[yw,yww,ywww]=dwell(cth,range,pattern);
y1=yw*s+r0+L;
y2=yw*s+r0 ;
line([0 0],[y1(oo),y2(oo)],'linewidth',10,'color','blue');
line([0 30],[y1(oo),20],'linewidth',10,'color','red');
axis([-30 30 -30 30])
hold on;
plot([0 x],[0 y],'ro',x,y,'k-')
pause(0.03);
end


主程式
==>
[x y]=pincam2([0:10:360],15,5,0,10,[100 200 260],[2 1],-1)

2007年5月25日 星期五

機動學 作業十

作業十 b94611042 王志豪 (Due date: 12pm May 23, 2007)
我有上本週5/17的課

我的blog


題目:

1.請聲明本週(5/17)有來上課。

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

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

10.1
我有上本週5/17的課



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

(1)
P點指的是桿三以及桿四相交的位置.
在此,我以r1,r2,r3,r4來分別表示各桿.
這時後,即可用老師上課所用的方法.首先,我們找出r1、r3的瞬時中心M13.
其位置會在r2、r4的延長線上,所以我們可以很清楚的看出來P點的速度以及加速度方向.
P的速度方向應該是在垂直r2的方向,而加速度是在r2的方向.














基本桿圖示意

若要比較科學化的分析,也可以如此看待:
一開始,先假設PM間距離為r(單位m).
而此時M就以等角速度ω(單位rad/s)轉動.
則P點此時有速度r*ω(m/s),其方向與r為垂直方向,且有向心加速度r*ω^2(單位m/s^2)










精簡化分析


(2)
那麼,假若該特定點M復以等速水平運動,則同一端點P之速度與加速度方向又會變為如何呢?
簡單來說,假若M複以V(m/s)等速水平移動,則P點此時之速度為:

==> V+ω*r ,其中方向為其向量合.

至於加速度方面則不受影響,依舊為r*ω^2,維持不變.














向量合力圖

(3)
而談到若M點同時也有加速度,則關於點P的變化,茲說明如下:
一言以蔽之,由於M點有了加速度,故此時P點的加速度將不再只有法線加速度而已.
它的切線方向理論上也會有一個切線加速度.
不過,在此和M有水平速度不一樣的地方在於M有無加速度這項條件.
其實並不會影響到P點的速度以及加速度方向!!!
它所能影響的,不外乎只有速度以及加速度的大小而已.
這其實就是物理上所學的,關於力的方向可以各自獨立計算.
加速度大致會變成:V(m/s)+r*ω^2*cos(θ)(m/s) i + r*ω^2*sin(θ)(m/s) j


(4)
藉由上面這些條件,我們便可以開始分析四連桿的運動了.
首先,關於速度的部份,可知P的速度是在垂直r2方向,Q則是在垂直r4的方向.
加速度的部份,P在r2的延長線上,Q則是位於r4的延長線上面.
P之速度方向是P與桿一端點連線的垂直方向,加速度方向則是P與桿一端點的連線方向.
Q的速度方向,則就是Q與桿一端點連線垂直,加速度方向則為Q與桿一端點的同一連線方向.

除此之外,若是桿三驅動的話,可以看出P的速度方向就是在垂直r3的方向,P的加速度就是在r3的方向上.

而當桿三有了加速度,假設為a......
此時桿三除了角速度W1以外還會有另一個隨時間變動的角速度at.
結果跟我們分析一開始桿二驅動的時候有點類似:
都是只會影響速度以及加速度的大小,而不影響速度與加速度的方向.














參考用的基本示意圖,可與解說對照了解

另外,以桿三為驅動桿,較為特別之處,在於我們要利用Q.R的瞬時中心M24!
有了M24,我們可以發現P和Q之間的速度以及加速度方向是一樣的.
只不過,由於兩者和瞬時中心間的距離不同,所以速度跟加速度的大小也會不同.



10.3

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

在此,我是參照課本上的一個function slider_solve 來改寫的.
以下是我變形過的new function

主程式
==>
function slider_draw(R,L,e)
%R=桿一長
%L=桿二長
%e=偏置量
the1=slider_limit(R,L,e)
the2=90
ang=linspace(the1,the2,100);
d=slider_solve(ang,R,L,e,1);

x=R*cosd(ang);
y=R*sind(ang);
for n=1:100
hold on
line([0,x(n),d(n)],[0,y(n),e]);
line([d(n)-3,d(n)+3,d(n)+3,d(n)-3,d(n)-3],[e-2,e-2,e+2,e+2,e-2]);
plot(0,0,'go')
plot(x(n),y(n),'go')
plot(d(n),e,'go')
plot([0,0],[0,e-d(n)*(y(n)-e)/(x(n)-d(n))],'go:')
plot([x(n),0],[y(n),e-d(n)*(y(n)-e)/(x(n)-d(n))],'go:')
plot([x(n),d(n)],[y(n),y(n)*d(n)/x(n)],'go:')
plot([d(n),d(n)],[0,y(n)*d(n)/x(n)],'go:')
axis equal
axis ([-110 110 -110 110]);
pause(0.08)
clf
end

%以上,程式的主要概念在於先畫出外型,然後再做出瞬心位置,最後plot上去.

以下即分別是
R > L
R < L
R = L
三種不同情況下曲柄滑塊連桿的運動狀況:









以上是關於曲柄滑塊連桿組合運動狀況的部份.
接下來,即便是關於計算瞬心的部份了.
基本上,瞬心以甘迺迪定理可定出.....
瞬心13可由12及23的連線14及34的連線交點獲得.
24則為23及34連線與12及14連線之交點獲得.
在此,我寫了一個專門關於計算瞬心位置的function

function [M12,M13,M14,M23,M24,M34]=Centor(r2,th2,r3,th3)
M12=[0,0]
M13=M23-[r2*sin(th2)/tan(th2),r2*sin(th2)]
M14=[r2*cos(th2)+r3*cos(th3),0]
M23=[r2*cos(th2),r2*sin(th2)]
M24=M23+[r3*cos(th3),r3*cos(th3)*tan(th3)]
M34=M23+[r3*cos(th3),r3*sin(th3)]

這個程式基本上,是可以幫我們找出六個瞬心的!
而為了將他們畫出來,由上面這個瞬心function,可以延伸出一個新function.

function [M12,M13,M14,M23,M24,M34]=PrintCenter(r2,theta2,r3,theta3)
M12=[0,0]
M13=M23-[r2*sin(th2)/tan(th2),r2*sin(th2)]
M14=[r2*cos(th2)+r3*cos(th3),0]
M23=[r2*cos(th2),r2*sin(th2)]
M24=M23+[r3*cos(th3),r3*cos(th3)*tan(th3)]
M34=M23+[r3*cos(th3),r3*sin(th3)]
x1=0;
x2=r2*cos(th2);
x3=x2+r3*cos(th3);
x4=r2*cos(th2)+r3*cos(th3);
x5=x2-(r2*sin(th2)/tan(th3))
x6=x2+r3*cos(th3)
y1=0;
y2=r2*sin(th2)
y3=y2+r3*sin(th3)
y4=0
y5=y2-r2*sin(th2)
y6=y2+(r3*cos(th3)*tan(th2))
Print1=[x1 y1;x5 y5;x3 y3;x4 y4;x1 y1;x6 y6;x4 y4;x1 y1]
line(Print1(:,1),Print1(:,2))
text(x1,y1,'M12')
text(x2,y2,'M23')
text(x3,y3,'M34')
text(x4,y4,'M14')
text(x5,y5,'M13')
text(x6,y6,'M24')

而此時,又可以再度利用上面這個瞬心程式,做出另一個動畫哩.

主程式
==>
AXIS([-40 40 -40 40]);
for th1=0:5:360
AXIS([-40 40 -40 40]);
PrintCenter(20,th1,15,85)
pause(0.2)
clf
end














瞬心圖舉例



2007年5月18日 星期五

機動學 作業九

作業九 b94611042 王志豪 (Due date: 12pm May 16, 2007)
我有上本週5/3的課

機動學分組:C組:5/17(四) 8:00am-9:30am
(蔡篤明,劉昶志,許惠善,張鈞崴,張延瑋,王志豪,范詠晴,林詠舜)

我的blog

題目:

•請聲明本週(5/3)有來上課。
•請閱讀機動論壇中有關下週分組上課的情形,並登記組別。

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

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

首先,題目所說的滑塊機構,其長相如下:














基本模式

由此,我們的分析正式展開.
一開始,我們得先討論題目所述之連桿,其運轉範圍之限制.
而在這邊,首先,我們就必須要先討論其限制的角度.
為了討論這些連桿的限制角度,網路講義第七章中,提供了一個很好用的函數.
這個函數,即為function 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(1) 固定桿(第一桿)
*r(2) 曲桿 =52(because 42+10)
*r(3) 結合(連)桿=57(because 52+5)
*r(4) 垂直偏置量=10
*theta1 : 第一桿之水平角
*theta2 : 驅動桿之水平夾角
*td2 : 驅動桿角速度
*tdd2 : 驅動桿角加速度
*sigma : 組合模數
*driver : 驅動桿

在這裡,基本上我們是不用設定第一桿的長度.
這主要是因為第一桿的方向雖然是固定的,但其長度卻會隨機構運動而不停變動伸縮.
其中,對於r的設定基本如下:
r=[第一桿、曲桿、連桿、偏置量]
因此,透過這些的條件,我們便可以把數值代入主程式中:

1.這是以第二桿作驅動桿的(driver為0),其結果如下:
***主程式==>
r=[0 52 57 10];
theta1=0;
linkdrive=0;
sld_angle_limits(r,theta1,linkdrive);

Qstart =

-64.6683


Qstop =

244.6683















第二桿驅動限制角度位置圖


2.接下來我們便以第三桿為驅動桿(driver為1)
***主程式==>
r=[0 52 57 10];
theta1=0;
linkdrive=1;
sld_angle_limits(r,theta1,linkdrive);

Qstart =

-47.4631


Qstop =

227.4631















第三桿驅動限制角度位置圖


3.接下來我們便以滑塊為驅動桿(driver為2)
***主程式==>
r=[0 52 57 10];
theta1=0;
linkdrive=2;
sld_angle_limits(r,theta1,linkdrive);

Qstart =

-108.5403


Qstop =

108.5403















滑塊驅動限制角度位置圖


由上面的結果,便可知道機構運動時的基本限制!!!



在了解極限角度的狀況後,我們便可以開始討論其運動狀況.
關於這個部分,我們可以利用一個講義第七章中的還是來幫助我們運算.
即是function drawsldlimits(r,th1,sigma,driver)!!!
這個函式正是能幫助我們算出運動界限.
桿長r為:(第一桿、曲桿、連桿、偏置量)
第一桿的水平角度為th1=0
Sigma如同第八次作業中的mode.
Driver則依照我們討論的情形有所更改.
當中,0代表驅動桿為第二桿,1代表驅動桿為第三桿,2代表由滑塊驅動.

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

有了這個程式後,我們便可以開始討論各項狀況了.



1.以第二桿為驅動桿的狀況(driver為0)

在此,由機動學的講義中,可以知道這部分應該會有四種可能性

1.是當r3≧r4≧0,且 r3+r4 ≦ r2
2.是當r4≧r3≧0,且 r3+r4 ≧ r2
3.是當r3≧-r4≧0,且 r3-r4 ≦ r2
4.是當-r4≧r3≧0,且 r3-r4 ≧ r2

從條件[0 52 57 10]裡面,可以看出沒有符合我們上面四個條件.
縱使把偏置量改成-10cm,[0 52 57 -10],依然無法符合上面四個條件的任何一種.
所以可得結論,這個滑塊機構是無法以第二桿為驅動桿的.



2.以第三桿為驅動桿的狀況(driver為1)

討論類似上面,分成四部份討論.

1.是當0≦r4,且0≦r2-r4≦r3
2.是當0≦r4,且r2+r4≦r3
3.是當r4≦0,且0≦r2+r4<r3
4.是當r4≦0,且r2-r4<r3 時

從條件[0 52 57 10]我們可以看出我們符合第一個條件.
所以這個機構是可行的.
此時,欲畫基本圖型,程式如下:

drawsldlinks([0 52 57 10],0,i,-1,0)

結果參照上面角度分析的圖.
即是"第三桿驅動限制角度位置圖"

這時候,為了方便觀察分析,可以製作成動畫型式來呈現.
***主程式==>
for i=-47.4631:2:227.4631
pause(0.1)
clf
drawsldlinks([0 52 57 10],0,i,-1,0)
axis([-130 130 -130 130])
end





3.以滑塊為驅動桿的狀況(driver為2)

在這個以滑塊為驅動桿的部份,基本上只會有兩種情形產生.
1.|r2-r3|>r4時,其條件限制於r1所給之長短
2.|r2-r3|< r4時
由我們的桿長[0 52 57 10],可知是第二種情形
所以這個機構是可行的.
此時,欲畫基本圖型,程式如下:

drawsldlinks([0 52 57 10],0,i,-1,2)

結果參照上面角度分析的圖.
即是"滑塊驅動限制角度位置圖"

這時候,為了方便觀察分析,可以製作成動畫型式來呈現.
(這個動畫向右運動的情形和向左是一樣的)
而此時的Qstart跟Qstop指的是第一桿的起點位置以及終點位置.
如果想要使用drawsldlinks來做成動畫的話,我們必須先算出角度.

θmin =sin-1(r4/(r2+r3)) = 5.264
θmax =π- sin-1(r4/(r2+r3)) = 174.736

***主程式==>
for i=5.264:(174.736-5.264)/60:174.736
pause(0.1)
clf
drawsldlinks([0 52 57 10],0,i,-1,1)
axis([-130 130 -130 130])
end




(註:本來以前一直解決不了動畫執行的時候,其範圍大小會不斷變化,造成視野縮放不停,使得觀賞上較不順眼.而在5/17的時候,由於一個不小心睡過頭,錯過了小組討論,然而上課時,感謝教授的指正與教誨,讓我終於明白原來問題是出在於for迴圈以及clf身上,經修正之後,視野果然不復過去不停縮放的窘境,非常高興,十分感謝!!!)




在大致上的討論告一段落以後......
接下來便可以來討論滑塊連桿之軌跡.
此時,可以利用講義上的function drawsldpaths.
這是個蠻厲害的函式.
在前面,任何桿之延伸點上之位移、速度及加速度等東西......
均可利用sldlink()、body()、fb_sld_linits()及drawlinks()等函數進行分析!!!
然而,利用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)];%Sq,Sp,Sa %找q,p,a
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);];
%Vs, Va, As, Aa
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;
*****************************************************************

關於這個函式,其基本使用方式如下:
* r : 四連桿之長度,一列向量表示.
* th1: 第一桿的水平角度(通常為零).
* th2: 驅動桿或聯結桿的水平角度.
* td2, tdd2:驅動桿之角速度及角加速度.
* Sigma:連桿之組合模式選擇(需為+1或-1)
* Driver:驅動桿之設定( 0 表示苗桿驅動,1 表示由聯結桿驅動,3滑塊驅動).
* Npts: 設定分割的點數或位置.
* r6,rh6,nlink:桿上特定點之位置,包括桿長,與桿之夾角及附於何桿.
* Mode:mode=0,畫簡單位置圖;=1 畫所有圖表;=2畫所有圖表,但用簡單位置圖.


***主程式==>
drawsldpaths(0,0,[0 52 57 10],0,0.01,0.01,1,50,0,1)















軌跡分析

















軌跡各項要素分析

2007年5月9日 星期三

機動學 作業八

作業八 b94611042 王志豪 (Due date: 12pm May 9, 2007)
我有上本週4/26的課
我的blog

題目:

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

* 問題一:設桿2角度theta2=45度時,求各點之位置、速度與加速度為何?
* 問題二:繪出此四連桿之相關位置及標明各點之速度方向及大小(以程式為之)。
* 問題三:當桿2迴轉時,求出此組四連桿之限制角度,並繪出其位置(以程式為之)。
* 問題四:設theta2=[0:20:360],試繪出此組四連桿之重疊影像,解釋為何有些沒有值。
* 問題五:若將問題三考慮在內,只在可迴轉的範圍內迴轉,請問你能讓此組四連桿作成動畫方式迴轉嗎?

8.1
==>設桿2角度theta2=45度時,求各點之位置、速度與加速度為何?
Ans:

首先,按題目所示,我引用了網路講義第六章的一個函式function 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
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(r,theta1,theta2,td2,tdd2,mode,linkdrive)
* r(1:4) = 各桿之長度,r(1)為固定桿,其餘分別為曲桿、結合桿及被動桿,即r=[r1 r2 r3 r4]。
* theta1 = 第一桿之水平角,或為四連桿之架構角,以角度表示。
* theta2 = 驅動桿之水平夾角,以角度表示。一般為曲桿角,但若為結合桿驅動,則為結合桿之水平夾角。
* td2 = 驅動桿(第二桿或第三桿)之角速度(rad/sec)。
* tdd2 = 驅動桿(第二桿或第三桿)之角加速度(rad/sec^2)。
* mode = +1 or -1. 組合模數,負值表示閉合型,正值為分支型,但有時需視實際情況而定。
* linkdrive = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿) 。
* form = 組合狀態, 0 :表示無法組合; 1:可以正確組合
* data = 輸出矩陣,其大小為 4 X 7,各行之資料分配如下:

**********************************************************
    1   2(deg) 3(rad/s) 4(rad/s2) 5 6   7
I   桿1位置   θ1 ω1    α1    VQ  |VQ|  ∠VQ
II  桿2位置 θ2   ω2    α2    VP  |VP|  ∠VP
III 桿3位置  θ3   ω3    α3    AQ  |AQ|  ∠AQ
IV 桿4位置   θ4   ω4    α4    AP  |AP|   ∠AP
**********************************************************

因此,我們便可以利用這個程式來解決問題.
首先,便是把題目所需條件代入程式中.

***主程式==>
[val,form]=f4bar([4 3 3 5],0,45,10,0,-1,0)
abs(val(:,1))
abs(val(:,2))
abs(val(:,3))
abs(val(:,4))
%我是以mode=-1來執行,也就是討論閉合型的狀況

此時,便可得解,而由於位置他是以向量的模式來表現,所以實際位置必須直接用手計算之.


各桿之速度(rad/s):
第一桿 : 0
第二桿 : 10.0000
第三桿 : 16.2681
第四桿 : 4.9677


各桿之加速度(rad/s^2):
第一桿 : 0
第二桿 : 0
第三桿 : 491.4428
第四桿 : 383.6120


P點的
位置 : (2.12 , 2.12)
速度 (rad/s) : 0.0212 + 0.0212i (10)
加速度 (rad/s^2) : -0.2121 - 0.2121i (0)


Q點的
位置 : (3.17 , 4.93)
速度 (rad/s) : 0.0041 - 0.0245i (16.2681)
加速度 (rad/s^2) : -1.8712 - 0.4391i (491.4428)


8.2

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

按題目所示,我引用了網路講義第六章的一個函式function drawlinks:

***********************************************************
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
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),' P');text(rx(3),ry(3),' Q');
else
fprintf('Combination of links fail at degrees %6.1f\n',th2);
end
axis equal
grid on
******************************************************************

這個程式的基本使用大致如下:
function drawlinks(r,th1,th2,mode,linkdrive)
* r(1:4) = 各桿之長度,r(1)為固定桿,其餘分別為曲桿、結合桿及被動桿。
* theta1 = 第一桿之水平角,或為四連桿之架構角,以角度表示。
* theta2 = 驅動桿之水平夾角,以角度表示。一般為曲桿角,但若為結合桿驅動,則為結合桿之水平夾角。
* mode = +1 or -1. 組合模式,負值表示閉合型,正值為分支型,但有時需視實際情況而定。
* linkdrive = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿)


以及引用之前作業用過的function dyad_draw
*******************************************************************
function dyad_draw(rho,theta,td,tdd)
% Call function dyad(rho,theta,td,tdd) to
% analyze a dyad linkage composing cranks and a dyad with drawing.
% Inputs: rho:length of links
% theta:incling angles, deg.
% td:angular velocity, rad/s
% tdd:angular acceleration, rad/s^2
% Example:
% dyad_draw([15 10 5],[0 60 90],[.5 .8 .5],0)
clf;
[vec,data] = dyad(rho,theta,td,tdd);
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
*********************************************************************

因此,我們便可以利用這兩個程式來解決問題.
於是,便把題目所需條件代入程式中.

***主程式==>
drawlinks([4 3 3 5],0,45,-1,0)















四連桿相關位置

drawlinks([4 3 3 5],0,45,-1,0)
dyad_draw([4 5],[0 99.5246],[0 4.9677]*0.05,[0 383]*0.001)
dyad_draw([3 3],[45 69.4856],[10 16.2681]*0.05,[0 491]*0.001)















連桿之速度方向與大小


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

關於這題,在此則是使用了極限函數fb_angle_limits:

*************************************************************
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="0;" 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, Ang2]=fb_angle_limits([4 3 3 5],0,0)

此時,我們便可以得到:
Ang1(即驅動桿與x軸之間最小夾角,start)= 28.9550度
Ang2(即驅動桿與x軸之間最大夾角,stop )= 331.0450度

因此,我們便可以利用這兩個角度來繪圖,同樣是用到function drawlinks

***主程式==>
[Ang1, Ang2]=fb_angle_limits([4 3 3 5],0,0)
drawlinks([4 3 3 5],0,28.9551,-1,0)
drawlinks([4 3 3 5],0,331.0449,-1,0)















四連桿之極限角度位置圖


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

在此,我將theta2設為i,用function drawlinks直接寫一個for迴圈即可

***主程式==>
for i=0:20:360
drawlinks([4 3 3 5],0,i,-1,0)
end















四連桿之重疊影像

跑完程式後,Matlab顯示了這些訊息:
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間才可做旋轉,因此其他角度是不會發生的.換個方式來講,沒有值的部份,其實就是沒有符合網路上第六章講義起先提到的:我們的四連桿必要滿足r2ejθ2+r3ejθ3=r1ejθ1+r4ejθ4,當四連桿滿足條件就代表實際情形中,它是可以轉到該角度的,故可得知我們上一題所算的限制角度,其實也就是運用這個閉合關係,完全是息息相關的.


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

當然能的.
在此,我們可以再次利用網路上的資源,也就是function move_4paths.
這是一個頗讚的動態模擬function,可以處理很多種軌跡.

*****************************************************************************
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;">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=f4limits(r,th1,sigma,driver)
%function f4lmits(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=draw4link(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=draw4link(values,driver);
else
fprintf('Combination of links fails at degrees %6.1f\n',Qstart);
end
axis equal
grid on
**************************************************************************

***主程式1==>
move_4paths([4 3 3 5],0,0,3,0,10,0,1,0,4,150)


使用function move_4paths做出的動畫


當然,假如不依賴這個強大的function move_4paths的話,也是可以用兩個for迴圈,然後把最大極限角度最小極限角度當作上下限,然後區分成50等份(等分隨意,任人高興,純粹是動畫流暢度的問題),然後再把各角度代入function drawlinks中,運用for迴圈每次清圖(使用clf),便可以做出動畫了!只不過,這個做出來的品質明顯沒有使用function move_4paths來得好,因此這個程式,顯然地過於簡陋,還有頗多進步的空間......(汗).

***主程式2==>
AXIS([-5 5 -5 5]);
for theta=28.9551:(331.045-28.9551)/50:331.045
clf
axis equal
AXIS([-5 5 -5 5]);
drawlinks([4 3 3 5],0,theta,-1,0)
pause(0.2)
clf
end
grid on
for theta2=331.045:-(331.045-28.9551)/50:28.9551
clf
axis equal
AXIS([-5 5 -5 5]);
drawlinks([4 3 3 5],0,theta2,-1,0)
pause(0.2)
end
grid on


使用兩個for迴圈所製成的動畫

2007年4月26日 星期四

機動學 作業七

作業七 b94611042 王志豪
我有上本週4/19的課

7.題目
設有三桿相連結成端桿(dyad)之型式。請利用網路講義中之函數dyad及dyad_draw作分析。

function dyad_draw(rho,theta,td,tdd)

其中:
* 各桿之對應長度rho=[a, a+5, a-5]cm,a=(你的學號末一碼)+10;
* 各桿之對應起始角度theta=[0, 0, 0]度;
* 各桿之對應角速度為td=[0.2, 0.5, 0 .3]rad/s;
* 各桿之對應角加速度為tdd=[0, 0.1, 0.2]rad/s^2;
問:
1. 當t=[1 2 3 4 5]秒時,此端桿之對應方位如何?
2. 繪出三桿之結點之速度,加速度與時間t之關係。
3. 能作出此時間區段之動畫嗎?

Ans:
(1)
首先,按題目所示,引用網路講義裡的兩個函式如下:

*******************************************************************
function [vec, dyadata] = dyad(rho,theta,td,tdd)
%
% function [vec, th, dyadata] = dyad(rho,theta,td,tdd)
% Analyzes a dyad linkage composing a crank and a dyad.
% Inputs: rho:length of links
% theta:incling angles, deg.
% td:angular velocity, rad/s
% tdd:angular acceleration, rad/s^2
% Outputs: vec:absoute length of links
% th:angles of links, velocities & acc, deg
% dyadata:original data, in complex forms
% Example:[vec,th,dyadata] = dyad([5 10],[30 50],[2 4])
theta=theta(:); rho=rho(:);
n=length(rho);
if nargin<4, tdd="zeros(size(rho));" nargin="=" td="ones(size(rho));" td="ones(size(rho))*td;end" tdd="ones(size(rho))*tdd;end" td="td(:);tdd=" d2g="pi/180;" tt="exp(i*theta*d2g);" pp="rho.*tt;" vv="i*td.*pp;" aa="-pp.*td.^2+i*pp.*tdd;" dyadata="[pp" vec="[abs(sum(dyadata));angle(sum(dyadata))/d2g];" x="[0;cumsum(real(data(:,1)))];y=" i="1:length(x)-1" k="1:length(rho)" x0="x(k+1);y0=" vx="x0+real(data(k,2));vy=" ax="x0+real(data(k,3));ay=" sdata="sum(data);" ss="[real(sdata(1))" vv="[real(sdata(2))" aa="[real(sdata(3))" nargin="=" dd="1;end;" d="abs(dd);" ab="(B(1)+j*B(2))-(A(1)+j*A(2));" d="abs(AB);th=" t="linspace(pi/2,2.5*pi,20);" cout="max(d/2,0.2)*exp(j*t');Cin=">0,
P=[0;Cin;Cout(1:10);D+Cout(11:20);D+Cin;D+Cout(20);Cout(1)];
else
P=[Cin;0;D;D+Cin];
end
xx=real(P);yy=imag(P);
x=xx*cos(th)-yy*sin(th)+A(1);
y=xx*sin(th)+yy*cos(th)+A(2);
line(x,y)
axis equal
************************************************************************
(註:個人決定再引用網路講義中另一個函數linkshape,這可使圖形美觀.)

dyad_draw([12,17,7],[0,0,0],[0.2,0.5,0.3],[0,0.1,0.2]);

for n=0:1:5;
dyad_draw ([12,17,7],[0.2*n,0.5*n+(1/2)*0.1*n^2,0.2*n+(1/2)*0.2*n^2],[0.2,0.5+n*0.1,0.3+0.2*n],[0,0.1,0.2]);
pause(1);
end;
%在此,我使用了一些數學運算代入dyad_draw這個function裡面.
%由於角速度以及角加速度隨時間而變,所以設一n變數進去,使其模仿時間.
%利用for迴圈,使他重複執行五遍,便可得下列圖片之結果.

當t=[1 2 3 4 5]秒時,此端桿之對應方位:














第1秒















第2秒















第3秒















第4秒















第5秒



(2)
第一桿之速度
第一桿之加速度




第二桿之速度
第二桿之加速度




第三桿之速度
第三桿之加速度



%以上繪圖來自於dyad_draw,利用其角速度以及角加速度的值去跑出來的.
%簡單來說,就是利用第一題的程式,把其得到的值代入另一個繪圖小程式.
%這個程式是用plot畫出來的,順著桿長最大極限而畫出的.
%類似這個:
for a=1:length(rho)
figure
plot(aa(j,:))
figure
plot(bb(j,:))
end
%當初禮拜六畫這個圖的時候,代數好像不是用這個(忘了),不過概念上大概長這個模樣就是了.

(3)
axis equal;
dyad_draw([12,17,7],[0,0,0],[0.2,0.5,0.3],[0,0.1,0.2]);
pause(2)
clf
for n=1:5
axis equal;
dyad_draw([12,17,7],[0.2*t,0.5*t+0.5*0.1*t^2,0.3*t+0.5*0.2*t^2],[0.2,0.5+t*0.1,0.3+t*0.2],[0,0.1,0.2]);
pause(1)
clf
end;
%與第一題的作法雷同,加了些暫停,使其運作更流暢.


hw7動畫