Homework 8

Homework 8 b94611001 歐陽太閒
本週四　（4／26）曾來上課。

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

一 ,設桿2角度theta2=45度時,運動學分析

採用老師機動學教學部落格中提供的f4bar程式(f4bar程式之輸出之第七行夾角並未輸出)註解後版本如附件一
其中:
r(1:4) = [4 3 3 5] %各桿之長度
theta1 = 0 deg %第一桿水平角
theta2 = 45 deg %驅動桿水平夾角
td2 = 10 %驅動桿(第二桿)角速度
tdd2 = 0 %驅動桿(第二桿)角加速度
mode = -1 %閉合型
linkdrive = 0 %驅動桿第二桿


>> r(1:4) = [4 3 3 5]
t1=0
t2=45
ang=10
aacc=0
f4bar(r,t1,t2,ang,aacc,-1,0)

ans =

1.0e+003 *

0.0040 0 0 0
0.0021 + 0.0021i 0.0450 0.0100 0
0.0011 + 0.0028i 0.0695 -0.0163 0.4914
-0.0008 + 0.0049i 0.0995 -0.0050 0.3836

0.0212 + 0.0212i 0.0021 + 0.0021i
0.0041 - 0.0245i 0.0032 + 0.0049i
-0.2121 - 0.2121i 0
-1.8712 - 0.4391i 0

可得以下分析:

即各桿之(*1000)
位置 角度(度數) 角速度 角加速度
桿1 0.0040 0 0 0
桿2 0.0021 + 0.0021i 0.0450 0.0100 0
桿3 0.0011 + 0.0028i 0.0695 -0.0163 0.4914
桿4 -0.0008 + 0.0049i 0.0995 -0.0050 0.3836

Q速度 0.0212 + 0.0212i Q位置 0.0021 + 0.0021i
P速度 0.0041 - 0.0245i P位置 0.0032 + 0.0049i
Q加速 -0.2121 - 0.2121i
P加速 -1.8712 - 0.4391i

f4bar之程式似乎與網頁之說明有所不同,程式中註解data (1,6),data (2,6)分別為Q,P位置,但網頁為VQ,AQ,VP,AP之絕對值

二,四連桿位置及各點速度

由於需繪出一四連桿系統並撰寫一程式繪製之,故參考機動學教學平台主題6.1至6.5關於四連桿系統的位置,速度,加速度之分析,並據此撰寫一

函式以利日後使用.

由於教學平台之範例程式為以直線來表現連桿,故試寫一函式嘗試以連桿圖形來顯示實際的運動狀況,基礎架構仍同於drawlinks.m,概述如下:
1.輸入防呆處理
2.提取f4bar分析的結果
3.進行計算,產生繪圖用數據
4.進行繪圖
5.數據標示及顯示(P,Q桿端速度)

執行結果如下:

script如下:
title('作業8.2 作者:歐陽太閒 b94611001');%資訊
xlabel('x-axis');
ylabel('y-axis');
dw4shape([4 3 3 5],0,45,10,0,-1,0);


dw4shape程式碼如下:
function [values]=dw4shape(r,th1,th2,td,tdd,mode,linkdrive)
%歐陽太閒 b94611001
%本函數改寫自老師的drawlinks.m
%[values]=dw4shape(四桿長,第一桿角度,第二桿角度,角速度,角加速度,模式,驅動桿)
%[values]=dw4shape([4 3 3 5],0,45,10,0,-1,0)
%本程式需要f4bar.m及linkshape.m方能正常執行
r2d=180/pi;
axis equal;
if nargin<7,linkdrive=0;end %防呆
if nargin<6,mode=1;end %防呆
[values b]=f4bar(r,th1,th2,td,tdd,mode,linkdrive); %提取f4bar中的值
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;

qv=abs(values(1,5)); %取得點Q速度
if real(real(values(1,5)))>0;dirqv=atan(imag(values(1,5))/real(values(1,5)))*r2d;else;dirqv=pi/2;end %方向
pv=abs(values(2,5)); %取得點P速度
if real(real(values(2,5)))>0;dirpv=atan(imag(values(2,5))/real(values(2,5)))*r2d;else;dirpv=pi/2;end %方向

if b==1
linkshape([0 0],[rx(1),ry(1)],0.5); %第一桿
hold on;
if linkdrive==0 %第二桿為驅動桿時
linkshape([0 0],[rx(2),ry(2)],0.5); %繪製第二桿
linkshape([rx(2) ry(2)],[rx(3),ry(3)],0.5); %繪製第三桿
else
linkshape([0 0],[rx(2),ry(2)],1); %繪製第二桿
linkshape([rx(2) ry(2)],[rx(3),ry(3)],0.5);%繪製第三桿
end
linkshape([rx(3) ry(3)],[rx(1),ry(1)],0.5);%繪製第四桿
end
text(0,0,' O');
text(rx(1),ry(1),' R');
text(rx(2),ry(2),' P');
text(rx(2)/2,ry(2)/2, ['Pv=' num2str(pv) 'θ= ' num2str(dirpv)]);%標示P速度
text(rx(3),ry(3),' Q');
text((rx(3)-rx(2))/2+rx(2),(ry(3)-ry(2))/2+ry(2),['Qv=' num2str(qv) 'θ=' num2str(dirqv)]);%標示Q速度
%程式結束

另外,若以圖形方式呈現(由於圖形速度加速度過大,故未完整顯示),可在dw4shape末端加入以下程式碼
藍色為Q點速度,紅色為P點速度,黃色為Q點加速,粉紅色為P點加速


可利用以下程式碼:
qvx=real(values(1,5))
qvy=imag(values(1,5))
line([rx(2) rx(2)+qvx],[ry(2) ry(2)+qvy],'color','b','linewidth',2) %藍色為Q點速度
pvx=real(values(2,5))
pvy=imag(values(2,5))
line([rx(3) rx(3)+pvx],[ry(3) ry(3)+pvy],'color','r','linewidth',2) %紅色為P點速度
qax=real(values(3,5))
qay=imag(values(3,5))
line([rx(2) rx(2)+qax],[ry(2) ry(2)+qay],'color','y','linewidth',2) %黃色為Q點加速
pax=real(values(4,5))
pay=imag(values(4,5))
line([rx(3) rx(3)+pax],[ry(3) ry(3)+pay],'color','m','linewidth',2) 粉紅色為P點加速
axis equal;grid on
axis([-30 30 -30 30]);



若採用老師的drawlinks

title('作業8.2 作者:歐陽太閒 b94611001');%資訊
xlabel('x-axis');
ylabel('y-axis');
drawlinks([4 3 3 5],0,45,-1,0);
則其結果為

三,四連桿之限制角度

由於四連桿運動機構中有死點存在,故需撰寫一程式進行分析找尋驅動桿可運動的範圍
而本程式中利用機動學教學部落格中的fb_angle_limits.m進行角度分析
其結構如下:
1.程式資訊
2.自fb_angle_limits取得角度分析的資訊
3.以此資訊進行繪製
4.數據標示


script如下:
>>dw4lim([4 3 3 5],0,-1,0)
angs =
28.9550 (最小角)
ange =
331.0450 (最大角)

執行結果:


function dw4lim(r,th1,sigma,driver)
%歐陽太閒 b94611001
%dw4lim(四桿長,第一桿角度,模式,驅動桿)
%本程式用於進行尋找四連桿系統的驅動桿最大角及最小角
%本程式需有 dw4shape.m 及 fb_angle_limits.m方可正確執行
title('作業8.3 作者:歐陽太閒 b94611001');%資訊
xlabel('x-axis');
ylabel('y-axis');
[angs, ange]=fb_angle_limits(r,th1,driver) %利用fb_angle_limits進行分析
dw4shape(r,th1,Qstart,0,0,sigma,driver); %程式碼再利用
dw4shape(r,th1,Qstop,0,0,sigma,driver);
text(r(1),r(1),['s1=' num2str(angs,'%6.1f')]); %取到小數下第一位並轉成字串
text(r(1),-r(1),['s1=' num2str(ange,'%6.1f')]);
%程式結束


而若以老師的drawlimits進行分析:
drawlimits([4 3 3 5],0,-1,0)
Qstart =
28.9550 (最小角)
Qstop =
331.0450 (最大角)

執行結果:


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

本程式繪出四連桿的連續運動並顯示其驅動桿之旋轉有一範圍,亦即沒有值處乃是來自於其在f4bar的分析結果中乃是不存在的(虛數)
是故,即如老師在課堂上所說的,其沒有值處是由於其解為虛數,在實數座標系中不能顯現.

script如下:

title('作業8.4 作者:歐陽太閒 b94611001');%資訊
xlabel('x-axis');
ylabel('y-axis');
for k=[0:20]
dw4shape([4 3 3 5],0,k*20,10,0,-1,0);
end

執行結果:


而若以drawlinks為之
script如下:
title('作業8.4 作者:歐陽太閒 b94611001');%資訊
xlabel('x-axis');
ylabel('y-axis');
for k=[0:20]
drawlinks([4 3 3 5],0,k*20,-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
Combination of links fail at degrees 380.0

即為fbar中傳回b不為1,不能組合的無效結果,也就是在驅動桿運動範圍外的結果.


五,動畫

根據老師的drawlimit執行結果:
drawlimits([4 3 3 5],0,-1,0)
Qstart =
28.9550
Qstop =
331.0450
故可知其可旋轉範圍為28.9550 2到331.0450之間

可寫一程式使其迴轉:
其結構如下
1.資訊
2.for-loop進行繪製

程式碼如下:
title('作業8.5 作者:歐陽太閒 b94611001');%資訊
xlabel('x-axis');
ylabel('y-axis');
r=[4 3 3 5]
for k=[29:331] %簡化程式碼
pause(0.01);
clf;
axis([-8 8 -5 5]); %固定座標系
dw4shape(r,0,k,0,0,-1,0) %採用先前撰寫的dw4shape
end

結果如下:

影片連結
經老師指導,本影片已經改進為坐標不移動


心得:
本次作業雖可參考教學平台之資料,但仍在實際處理上遇到很多困難,其大多來自於角度分析時定義的限制,所以必須要使用複雜的if-then-else

結構來進行處理,而四連桿系統的運動本就複雜難以分析,故前人累積的智慧,使我們不必重新從頭想起,而可以以現成的公式,定理等進行處理.

而本次作業也已經進入到分析實際運用於生活中的結構,表面上看來簡單,實際上卻極為複雜.

附件一
%歐陽太閒 b94611001 加註
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 %防呆裝置 補足驅動
if nargin<6,mode=1;end %防呆裝置 補足模數
data=zeros(4,6); %建立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 %驅動桿為第三桿時,r的順序交換
rr=r.*r;d2g=pi/180;
[theta,td,tdd]=deal(zeros(4,1)); %存到theta,td,tdd
theta(1:2)=[theta1 theta2]*d2g; %theta1,theta2的輸入為角度量非徑度量
t1=theta(1);tx=theta(2);
s1=sin(t1);c1=cos(t1); %分別取t1,tx的水平和垂直分量
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, %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 分母為零先計算theta(3)
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),theta(4)計算完畢

theta(3)=t3;theta(4)=t4;


%velocity calculation %速度分析計算
td(2)=td2;
AM=[-r(3)*s3, r(4)*s4; -r(3)*c3, r(4)*c4]; %r3,r4分別乘上其對應角度的sin,cos化為分量
BM=[r(2)*td(2)*sx;r(2)*td(2)*cx]; %乘上角速度和方向分量
CM=AM\BM; %求解([RM] [AM] = [CM]公式,6.2)
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; %求解([DM]=[RM]\[FM])
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; %處理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

%第八次作業結束

Homework 7

%homework 7 作者:歐陽太閒 b94611001

7.1. 繪製位置

根據題目要求:
r=[a a+5 a-5]cm，a=1+10
長度起始角度 angle=[0 0 0]度
角速度 w=[0.2 0.5 0.3]rad/s(原題目的0.3多一空格)
角加速度 wa=[0, 0.1, 0.2]rad/s^2

而當時間改變時角速度變量(w0為原角速度,wa為角加速度):
delw=wa*t+w0
故將其對時間積分得
deltheta=0.5*w*(t)^2+w0*t
即為角度改變量

因此利用老師的dyad_draw列成以下script file(老師的dyad.m最後一行的vect似乎應為vec)
由於本題只要求位置,故dyad_draw中繪製速度,加速度的line部份暫以%使其成為註解
直接將各桿角速度,角加速度之參考點視為原點時:
結果如圖
程式碼如下:

%作業7.1 考慮角速度,角加速度參考點.歐陽太閒 b94611001
axis equal; %等比例
clf;
title('作業7.1 作者:歐陽太閒 b94611001 單位:cm,rad/s,rad/s^2'); %資訊
xlabel('x-axis');
ylabel('y-axis');
a=1+10; %1+10
r=[a a+5 a-5]; %各桿長
ang=[0 0 0]; %預設角度
w1=[0.2 0.5 0.3]; %角速度
aa=[0 0.1 0.2]; %角加速度
r2d=180/(pi); %角徑轉換
for t=[1 2 3 4 5] %迴圈
% for t=[1 1.5 2 2.5 3 3.5 4 4.5 5]
delw=w1 + t.*aa; %新角速度
deltheta=0.5.*aa.*(t.^2)+ w1.*t %各桿分別的角度增加量
ang2=deltheta.*r2d; %角度增加量 = 各桿對原點所轉角度
dyad_draw(r,ang2,delw,aa);
end


以上的結果有些奇怪,因此若考量dyad_draw的各桿角度是指對原點的絕對角度,而非對前一桿的角度,
但角速度造成的角位移則是對角速度對應之原點造成角位移(否則其桿的對應原點將移動),因此在第二三桿的位移增加量上需考慮其上一桿的位

置,但角速度仍可直接計算.
而由於各桿之角速度和角加速度皆為定值,故可寫一script file如下:
結果如圖
若啟用速度,加速度繪製
結果如圖
程式碼如下:

%作業7.1 考慮角速度,角加速度參考點.歐陽太閒 b94611001
axis equal; %等比例
clf;
title('作業7.3 作者:歐陽太閒 b94611001 單位:cm,rad/s,rad/s^2'); %資訊
xlabel('x-axis');
ylabel('y-axis');
a=1+10; %1+10
r=[a a+5 a-5]; %各桿長
ang=[0 0 0]; %預設角度
w1=[0.2 0.5 0.3]; %角速度
aa=[0 0.1 0.2]; %角加速度
r2d=180/(2*pi); %角徑轉換

for t=[1 2 3 4 5]
%固定座標系
delw=w1 + t.*aa; %各桿角速度
deltheta=0.5.*aa.*(t^2)+w1*t; %各桿分別的角度增加量
ang2(1)=(deltheta(1))*r2d; %1桿實際絕對角度
ang2(2)=ang2(1)+(deltheta(2))*r2d; %2桿實際絕對角度
ang2(3)=ang2(2)+(deltheta(3))*r2d; %3桿實際絕對角度
dyad_draw(r,ang2,delw,aa);
end


以上程式架構皆為:
1.座標定義和資訊
2.定義參數
3.以for-loop給予時間
a.分別以公式計算角速度和位移角度
b.以dyad_draw畫出
4.結束

7.2繪出三桿之結點之速度，加速度與時間t之關係

由於本程式以dyad為發展基礎,故直接採取其內對於速度加速度之定義,以以下程式求得值帶入dyad並繪製.
以下程式架構皆為:
1.座標定義和資訊
2.定義參數
3.以for-loop給予時間
a.分別以公式計算角速度和位移角度
b.以dyad求得速度,加速度,合速度,合加速度等
c.存入變數
4.以for-loop給予陣列序
a.繪製圖形
5.結束

1.利用以下程式求得各桿端速度與t之關係:
結果如圖
以下為程式碼:

%作業7.2.1 各桿加速度 歐陽太閒 b94611001
a=1+10; %1+10
r=[a a+5 a-5]; %各桿長
ang=[0 0 0]; %預設角度
w1=[0.2 0.5 0.3]; %角速度
aa=[0 0.1 0.2]; %角加速度
r2d=180/(pi); %角徑轉換
title('作業7.2 作者:歐陽太閒 b94611001 綠,紅,粉紅分別為桿1,2,3速度'); %資訊
xlabel('x-axis');
ylabel('y-axis');

for j=[1:1:41]
t=j/10+1
delw=w1 + t.*aa; %各桿角速度
deltheta=0.5.*aa.*(t^2)+w1*t; %各桿分別的角度增加量
ang2(1)=(deltheta(1))*r2d; %1桿實際絕對角度
ang2(2)=ang2(1)+(deltheta(2))*r2d; %2桿實際絕對角度
ang2(3)=ang2(2)+(deltheta(3))*r2d; %3桿實際絕對角度
[vec,data] = dyad(r,ang2,delw,aa);%利用dyad
x=[0;cumsum(real(data(:,1)))]; %取值
y=[0;cumsum(imag(data(:,1)))];
sdata=sum(data); %將三桿值和
for k=1:length(r)
vl(k,j)=abs(data(k,2));
al(k,j)=abs(data(k,3));
%line([x0 vx],[y0 vy],'marker','s','linewidth',2);
%line([x0 ax],[y0 ay],'marker','s','color','r','linewidth',3)
end
end
for k=[1:1:40] %繪製
line([k/10+1 (k+1)/10+1],[vl(1,k) vl(1,k+1)],'color','g','linewidth',2) %桿一速度
line([k/10+1 (k+1)/10+1],[vl(2,k) vl(2,k+1)],'color','m','linewidth',2) %桿二速度
line([k/10+1 (k+1)/10+1],[vl(3,k) vl(3,k+1)],'color','r','linewidth',2) %桿三速度
end


2.各桿端加速度與t之關係
結果如圖
以下為程式碼:

%作業7.2.2 各桿加速度 歐陽太閒 b94611001
a=1+10; %1+10
r=[a a+5 a-5]; %各桿長
ang=[0 0 0]; %預設角度
w1=[0.2 0.5 0.3]; %角速度
aa=[0 0.1 0.2]; %角加速度
r2d=180/(pi); %角徑轉換
title('作業7.2 作者:歐陽太閒 b94611001 綠,紅,粉紅分別為桿1,2,3加速度'); %資訊
xlabel('x-axis');
ylabel('y-axis');

for j=[1:1:41]
t=j/10+1
delw=w1 + t.*aa; %各桿角速度
deltheta=0.5.*aa.*(t^2)+w1*t; %各桿分別的角度增加量
ang2(1)=(deltheta(1))*r2d; %1桿實際絕對角度
ang2(2)=ang2(1)+(deltheta(2))*r2d; %2桿實際絕對角度
ang2(3)=ang2(2)+(deltheta(3))*r2d; %3桿實際絕對角度
[vec,data] = dyad(r,ang2,delw,aa);%利用dyad
x=[0;cumsum(real(data(:,1)))]; %取值
y=[0;cumsum(imag(data(:,1)))];
sdata=sum(data); %將三桿值和
for k=1:length(r)
vl(k,j)=abs(data(k,2)); %取速度
al(k,j)=abs(data(k,3)); %取加速度
end
end
for k=[1:1:40] %繪製
line([k/10+1 (k+1)/10+1],[al(1,k) al(1,k+1)],'color','g','linewidth',2) %桿一加速度
line([k/10+1 (k+1)/10+1],[al(2,k) al(2,k+1)],'color','m','linewidth',2) %桿二加速度
line([k/10+1 (k+1)/10+1],[al(3,k) al(3,k+1)],'color','r','linewidth',2) %桿三加速度
end


3.合速度及合加速度與t之關係:
結果如圖
以下為程式碼:

%作業7.2.3 各桿加速度 歐陽太閒 b94611001
a=1+10; %1+10
r=[a a+5 a-5]; %各桿長
ang=[0 0 0]; %預設角度
w1=[0.2 0.5 0.3]; %角速度
aa=[0 0.1 0.2]; %角加速度
r2d=180/(pi); %角徑轉換
title('作業7.2 作者:歐陽太閒 b94611001 綠色為合速度'); %資訊
xlabel('x-axis');
ylabel('y-axis');
for j=[1:1:401]
t=j/100+1
delw=w1 + t.*aa; %各桿角速度
deltheta=0.5.*aa.*(t^2)+w1*t; %各桿分別的角度增加量
ang2(1)=(deltheta(1))*r2d; %1桿實際絕對角度
ang2(2)=ang2(1)+(deltheta(2))*r2d; %2桿實際絕對角度
ang2(3)=ang2(2)+(deltheta(3))*r2d; %3桿實際絕對角度
[vec,data] = dyad(r,ang2,delw,aa);%利用dyad
x=[0;cumsum(real(data(:,1)))]; %取值
y=[0;cumsum(imag(data(:,1)))];
sdata=sum(data); %將三桿值和
vv(j)=[abs(sdata(2))]; %取速度加速度絕對值
a2(j)=[abs(sdata(3))];
end
for k=[1:1:400] %繪製
line([k/100+1 (k+1)/100+1],[vv(k) vv(k+1)], 'marker','d','color','g','linewidth',3) %合速度
line([k/100+1 (k+1)/100+1],[a2(k) a2(k+1)], 'marker','d','color','m','linewidth',3) %合加速度
end


7.3動畫
本動畫利用7.1的程式修改而成,提高程式碼的重複利用率.
連結於此
採取1~5秒動作

以下為程式碼:

%作業7.3 動畫 歐陽太閒 b94611001
axis equal; %等比例
clf;
title('作業7.3 作者:歐陽太閒 b94611001 單位:cm,rad/s,rad/s^2'); %資訊
xlabel('x-axis');
ylabel('y-axis');
a=1+10; %1+10
r=[a a+5 a-5]; %各桿長
ang=[0 0 0]; %預設角度
w1=[0.2 0.5 0.3]; %角速度
aa=[0 0.1 0.2]; %角加速度
r2d=180/(pi); %角徑轉換

for t=[1:0.1:5]
pause(0.1);%時間間隔
clf;
axis([-40 40 -40 40]); %固定座標系
delw=w1 + t.*aa; %各桿角速度
deltheta=0.5.*aa.*(t^2)+w1*t; %各桿分別的角度增加量
ang2(1)=(deltheta(1))*r2d; %1桿實際絕對角度
ang2(2)=ang2(1)+(deltheta(2))*r2d; %2桿實際絕對角度
ang2(3)=ang2(2)+(deltheta(3))*r2d; %3桿實際絕對角度
dyad_draw(r,ang2,delw,aa);
end


以上程式架構為:
1.座標定義和資訊
2.定義參數
3.以for-loop給予時間
a.延時清除
b.分別以公式計算角速度和位移角度
c.以dyad_draw畫出
4.結束

心得:本次作業在期中考期間,時間較緊迫,但是在程式寫作和分析時題目仍相當有內涵,故花費了較以往作業更長的時間,而本次作業讓我們了解
到關於具有角速度,角加速度的自由擺動桿之運動分析,而利用老師提供的dyad和dyad_draw乃節省了大量的時間,免去重複開發和除錯的痛苦.

Homework 6

homework6 歐陽太閒 b94611001

6.1平面組合機構

6.1圖

1.此圖共有12桿(L1(j2於其上不運動應視為一桿),L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12),其中L10為固定(

接地)桿
且運動結中j1共1個運動結為旋轉結(R)DOF=1,
j2,j3,j4共3個運動結為複式結,其中j4為滑塊結(DOF=2),j2滑動結(DOF=2),j3滑槽結視為稜柱結DOF=1.
而M1,M2,M3,M4,M5,M6為多桿共結,故其所佔結數=(相連桿數-1),
故M1=(3-1)=2,M2=(3-1)=2,M3=(3-1)=2,M4=(3-1)=2,M5=(3-1)=2,M6=(3-1)=2
因此總運動結數:
J=1+1+2+6*2=16
Σf=1+1+2*2+6*2=18
2.古魯伯公式:
M=3*(N-1)-(3*J-Σf)
M=可動度
N=系統的連桿(link)總數
J=系統的運動結(joint)總數
fi=i號運動結的連結度
故利用古魯伯公式代入(N=12,J=16,f=18):
M=3*(N-1)-(3*J-Σf)=3*(12-1)-(3*16-18)=33-30=3

3.題目要求採用老師撰寫的function[df]=gruebler()
故可得知其呼叫指令為
>> gruebler(12,[14 0 2])
其結果為
ans =
3


另外一種算法為視M6非三桿共結,則可得一答案為(J=15,f=3*2)
M=3*(N-1)-(3*J-Σf)=3*(12-1)-(3*15-17)=33-28=5
其呼叫指令為
>> gruebler(12,[13 0 2])
ans =
5

4.至於j4滑塊結(外力閉合),j2滑動結(外型閉合),j3滑槽結(外型閉合)
對於本系統的可動度影響為由於前二個運動結皆屬於DOF=2之複式結(轉動,沿特定方向滑動)
而第三個滑槽結根據12日請教老師的結果,應為由於接了兩個桿而可視為一個稜柱對(DOF=1)和一個三桿相連旋轉節,共使得DOF=3故使得本系統之可動度較純使用旋轉對的系統多了數個自由度,且使其可以受和系統可動度(3)等量的外部輸入而作動,而非一可能有預載負荷的靜不定不可動結構(DOF=-1).

6.2立體組合機構

6.2圖

1.本系統共有6桿(L1,L2,L3,L4,L5,L6),其中L1上套有R1但仍接地為固定桿,而有6個運動結,R1,R2為旋轉結

(DOF=1),而S1,S2,S3為球結(DOF=3),C1為筒結(DOF=2)
故其Σf=2*1+3*3+1*2=13

2.若由空間古魯伯公式計算
N=6
J=6
Σf=13
M=6*(N-1)-(6*J-Σf)=6*(6-1)-(6*6-13)=7
為一可動結構,但若考慮惰性自由度則為5
3.若利用老師的function[df]=gruebler()函式
則其呼叫指令為:
>>gruebler(6,[2 0 0 3 1])
結果:
ans =
7
4.根據定義:惰性自由度(Idle Degree of Freedom)
是由於球狀結及圓筒結而使得連桿形成自轉運動(例如筒結的軸轉動)而抵銷1個自由度稱之(可由觀察而得).
因此在S1,S2間的L5會產生自轉,而S3和C1之間的L4也會產生自轉,依照定義共會少2個自由度,是為惰性自由度

,故實際機構上的可動度為7-2=5,此惰性自由度造成系統上的實際可動度較古魯伯公式所推得者為少,使得機

構之可動度較所推得者少.

6.3葛拉索機構

1.
a.根據葛拉索準則,葛拉索機構為一四連桿機構,當一機構中的最長桿與最短桿的和小於另兩桿的和,則其至少

一桿是可旋轉桿(可完全迴轉).若表示成不等式則為:s+g＜p+q
b.至於非葛拉索型機構,則是最長桿與最短桿的和大於另兩桿的和,則除驅動桿外三活動連桿皆為搖桿,即為非

葛拉索型機構,表成數學式為s+g＞p+q

2.(圖形繪製用程式碼參考附件一)
第一組：桿1-桿4分別為7,4,6,5cm
g=7,s=4,p=5,q=6
s+g=11=p+q
為中立連桿組(Neutral Linkage),葛拉索變點機構,可能有兩種運動形式(平行四邊形機構和平行四邊形反置

機構,需利用另一個平行四邊形機構方可避免)
6.3.1圖


第二組：桿1-桿4分別為5,7,6,3cm
g=7,s=3,p=5,q=6
s+g=10＜p+q=11
為葛拉索一型機構,第一桿為固定桿,為一Crank-Rocker Linkage機構
6.3.2圖



第三組：桿1-桿4分別為5,4,4,3cm
s+g=8=p+q,為一中立連桿組(Neutral Linkage),葛拉索變點機構,有兩種運動形式(如第一組)
6.3.2圖



3.利用grashof()函式檢驗
第一組
>> grashof(1,[7 4 6 5])
ans =
Neutral Linkage
第二組
>> grashof(1,[5 7 6 3])
ans =
Crank-Rocker Linkage
第三組
>>grashof(1,[5 4 4 3])
ans =
Neutral Linkage
以上三組皆和前述分析同,檢驗成功

4.若以以上三組要成為葛拉索機構之改善方法如下:
第一組:為中立連桿組,故可以加長桿3或桿4使s+g＜p+q其成為一葛拉索機構,
或者加上另一平行四邊形運動機構,而可以避免態點的存在,如圖(摘自老師所著課本)

第二組:為一Crank-Rocker的葛拉索機構.

第三組:為中立連桿組,故可以加長桿2或桿3使s+g＜p+q其成為一葛拉索機構,
或者加上另一平行四邊形運動機構.

心得:本次作業在使用古魯伯公式時體會,在工程的分析上不能一切照公式硬帶,須臨機應變,例如多桿共結或者是複式結的處理皆須另外考慮,而在處理葛拉索機構時,方知簡單的四連桿機構分析起來竟然如此複雜,而且在解位置需要大量的計算,但最後能推得各種組合的運動模式,由此可見前人的智慧.

附件1:link4(l)程式碼
function link4(l)
%作業6.3.2 顯示用程式碼 作者 歐陽太閒
%本程式為便利繪出三種四連桿之機構,故參考老師所著之程式範例4.7-function four_link2

(theta2,r,mode)
%並嘗試進行簡化改寫,以利使用.
clf; %清除畫面
axis equal; %軸等比例
title('作業6.3 作者:歐陽太閒 b94611001'); %資訊
xlabel('x-axis');
ylabel('y-axis');
d2r=pi/180; %角徑轉換因子
ang=(0:15:360)'*d2r; %形成一組軌跡所需角
l2=l.^2; %開始建立常用之多項式,根據老師所著課本之解四連桿位置及參考程式範例4.7
Bx=l(2)*cos(ang);%解B結之x位置
By=l(2)*sin(ang);%解B結之y位置
m=(l2(4)-l2(3)-l2(1)+l2(2))./(Bx-l(1))/2;m2=(m-l(1)).^2;
p=(By./(Bx-l(1)));p2=p.^2;
Cy=((m-l(1)).*p+sqrt(m2.*p2-(p2+1).*(m2-l2(4))))./(p2+1); %解Cy之位置
Cx=m-(By./(Bx-l(1))).*Cy; %解Cx之位置
Cx=l(1)+sqrt(l2(4)-Cy.^2); %解Cx之位置
for i=1:length(ang)
if isreal(Cy(i)) %去除虛數解
line([0 Bx(i) Cx(i) l(1) 0]',[0 By(i) Cy(i) 0 0]'); %實際繪製
end
end
%程式結束